An automatic identification method for width of shear band of sand in PFC simulations
-
摘要: 从细观尺度研究颗粒材料的应变局部化有助于深刻理解其物理力学机制，而剪切带的出现和演化是应变局部化的典型特征。目前在离散元法（DEM）模拟中对砂土应变局部化特征（如剪切带宽度、倾角、分布类型等性质）的观测方法还主要依赖手动、目测和估计等方式。针对这一问题，采用PFC2D开展了大量密砂双轴压缩模拟试验，提出了一种基于“α-shape”和“Delaunay三角化”算法的剪切带及其宽度的自动识别方法。进而，基于多组砂土试样开展双轴压缩试验，验证了所提出的剪切带及其宽度自动识别方法的有效性。Abstract: Studying the strain localization of granular materials from the meso-scale is helpful for deeply understanding the physical and mechanical natures of strain localization, and the occurrence and evolution of shear bands are the typical characteristics of strain localization. In the numerical simulations of discrete element method (DEM), observing and measuring the evolution of strain localization of sand (such as the width, inclination angle and distribution type of shear bands, etc.) greatly rely on the manual, visual and estimation ways. Towards this problem, many numerical experiments on the dense sands in biaxial compression are carried out by using PFC2D. Based on the numerical experiments, an identification method based on "α-shape" and "Delaunay triangulation" algorithms is developed to automatically identify the shear bands and their widths. Furthermore, based on some groups of biaxial compression tests on the dense sands, the effectiveness of the proposed automatic identification method is validated.
-
0. 引言
位于寒冷地区的水利工程，地基易受极端低温影响而发生冻胀现象，进而造成对水工建筑物结构的破坏。寒区工程超重力试验平台是利用离心力场模拟超重力场，具有克服小比尺模型的应力差异过大，避免了原位试验复杂耗时的优势，且具备热交换、干/湿环境模拟的功能，是开展寒区水利工程模型试验的重要设备。采用该设备进行寒区工程地基冻胀问题的模拟研究时，模型上方的热交换系统可在一定时间内持续为模型提供冷/热边界，利用热交换系统与模型表面间的之间的空气夹层上下表面温差进行热交换。
开展超重力场下寒区模型试验时，试验设计需了解土体与热交换系统间夹层空腔温度变化。而对空气夹层的换热过程可简化为不同尺寸的空腔及土体的耦合对流换热问题。对于空腔结构内部重力导致的对流问题，石峰等[1-2]研究了矩形腔体中Rayleigh-Benard对流结构，Ning等[3]研究了具有水平流的流体对流的周期行为问题，胡宇鹏等[4]研究了腔体内对流的稳定性问题，传统对此类的问题研究以常重力状态下为主[5-7]，对于超重力环境对低温边界下空气夹层-结构对流换热问题研究较少。由于土工模型箱内空气夹层换热过程较为复杂，易受到离心机工作状态的影响，获得超重力对低温对流换热的影响规律是设计试验面临的挑战问题。
本文基于数值模拟开展土体-夹层在超重力条件下对流换热仿真，通过分析土体界面的温度变化以及空腔内换热过程，获得模型表面温度随时间的变化规律，以及空腔内的温度场分布，获得了对换热过程对流、传导相关无量纲数变化的认识，可为开展寒区土体模拟试验提供参考。
1. 夹层传热计算模型
寒区工程超重力试验设备可简化描述为封闭空腔内自然对流传热的物理模型：模型箱内盛放一土体模型，此时模型加速度为Ng（常重力的倍数），方向由离心机主轴指向转臂。热交换板为冷表面，初始温度为T0；模型表面为光滑平面（热面），初始温度为T′0；模型表面与模型箱下表面的距离为L（夹层间隙厚度），该距离内形成一矩形密闭空腔。超重力场下的密闭空腔换热模型简化如图 1所示。
图 1中土体除上表面外均为绝热，空气夹层上表面为冷边界，土体与夹层界面通过对流及热传导进行热交换。可知空腔内流动包含两种运动：沿固壁的边界层运动以及围绕空腔中心的旋转运动，后者具有典型的Rayleigh-Benard（简称RB）对流的特征。RB对流是研究非平衡开放系统中对流运动的稳定性、时空结构、涡环（又称斑图）形成及非线性动力学特性的典型模型之一[1-4]，指在一个封闭的空腔内，其上表面温度恒定，下表面加热，由于底部流体被加热膨胀，相对于上部温度较低的流体其密度较小，下部流体上升，在上升中接触其他流体部分时消耗能量，温度变低密度变大。当底部温度继续升高，某些底部流体的温度足够大，密度足够小，能够上升到顶部而能量不至于全部耗散，从而形成温度差导致腔体内流体运动的流动现象[6]。需要指出的是，存在临界瑞利数（Ra），当空腔内局部Ra高于临界值时，才会出现RB对流[5-7]，如图 3所示，否则仅为热传导状态。
![]()
图 2 空腔界面温度在不同Ng、不同间隙时随时间变化Figure 2. Variation of temperature at interface with time under different Ng and clearances![]()
图 3 20 cm空腔50g状态换热末态时刻温度分布Figure 3. Temperature nephogram of 20 cm thick clearance at overload of 50g and last time step该物理问题与常重力场下的Rayleigh-Benard对流相似，仅重力加速度不同。在Boussinesq近似假设（认为浮力诱导的流体运动中，当温度足够小时，仅在浮力项中考虑密度变化）下，在Ng重力场下的流体力学基本控制方程如式（1）所示：
∇U=0 ，∂U/∂t+(U⋅∇)U=υ∇2U−∇P/ρ0+ρ/ρ0Ngg ，(U⋅∇)T+∂T/∂t=λ/ρ0CP∇2T 。} (1) 式中U为流体速度场；T为温度场；ρ为空气密度；g为重力加速度；υ为空气运动黏度；Cp为空气比热；t为时间；λ为导热系数。由于直接求解物理模型较为复杂，一般采用基于CFD（计算流体动力学）数值方法开展研究。
根据以上分析，模型箱的边界条件可以简化为：矩形方腔左右两侧绝热，上部为空气，下层为土体（密度为1.80 g/cm3，导热系数为2 W/m·℃，热容为4 MJ/m3·℃，热扩散系数为2.5×10-7 m2/s）。问题可简化为在Ng的超重力场下，研究矩形密闭空腔内的温度变化及温度场分布。考虑到空腔温度变化和从土体交换的热量直接相关，土体的热量与体积、质量等参数密切相关，因此建立考虑空气夹层与土体的温度分析模型，以模拟全场温度变化。同时，考虑到模型箱为长方体结构，取二维模型（侧向剖面）以表征真实三维结构，以降低计算量，提高仿真效率。计算模型的边界条件为：土体与夹层界面为换热边界，夹层上表面为-30℃，其它界面均为绝热边界，全场初始温度为10℃；即T0= -30℃；T′0=10℃，侧壁及底面均为绝热边界。根据试验状态取常重力及两种超重力状态，即Ng =1，10和50；根据土体情况取夹层厚度L =1，10和20 cm；针对两种参数组合共9种工况开展分析。
2. 夹层-土体对流换热分析
2.1 模型表面温度随时间变化
给出不同空气夹层厚度时界面温度在不同Ng值及间隙下随时间变化，如图 2所示。
图 2结果温度变化趋势与文献[8]一致。分析表明：随着Ng增大，界面降温加速，高Ng下降温速率大于低Ng工况，以20 cm为例，1g（常重力）、10g和50g状态下前2×105 s的平均降温速率分别为13.5，15，16.5℃/105 s；随着换热时间变长，降温速度随之降低，换热过程减速；降温速率与Ng并非线性关系：1，10，50g状态下200000~400000 s的平均降温速率分别为3.5，3，2℃/105 s；原因在于随着土体温度变低，土体界面与空气温差变小，导致温差驱动的对流换热的热量变小，导致降温过程随之减缓。
其次，在超重力状态下，随着夹层间隙减小，降温速率提升，小间隙下降温速率大于大间隙：以1 cm为例，1g，10g和50g状态下前2×105 s平均降温速率分别为12.5，15.5，17.5℃/105 s。空气与土体界面的换热热流率随间隙增大而减小，随Ng上升而上升：如10cm间隙时，初始阶段1g，10g和50g状态下，提取界面换热热流分别为78，165，325 W/m2；1 cm间隙时，10g和50g状态下界面热流率分别为520，630 W/m2，大于间隙较大时。界面对流换热系数变化趋势也一致。
需要指出的是，由于模型箱总容积一定，当夹层厚度较大时（20 cm），土体质量较小（总热容小），夹层厚度较小时，土体质量大（总热容大），导致通过对流及热传导将土体温度降低所需热量不同。
2.2 温度场分布特性
提取典型时刻全场温度分布及空腔流动速度，取具有代表性的20 cm空腔厚度和50g工况为例，如图 3所示，图 4给出了典型状态空腔内涡环的状态。
![]()
图 4 换热末态10 cm空腔50g状态涡环Figure 4. Vortex ring of 10 cm-clearance and overload 50g at last time step结合数值仿真结果及图 3，4可以看出：①20 cm间隙时，空腔内呈现典型的RB流动特征，形成两个对称的涡环，涡环外部流速较快，内部流速越来越小；出现涡环，超重力使得流动速度变快，1g，10g和50g下涡环中最大速度分别为0.028，0.65，1.02 m/s；在靠近壁面角落部位，会形成小的涡环，随着空腔流场的发展而变化。②10 cm间隙时，不同Ng下均会出现2对涡环；靠近侧边壁面处涡环较细长，中部涡环长宽比趋于1。③1 cm间隙时，1g状态有明显涡环成对出现；10g和50g状态下无明显的涡环发展过程，仅在靠近壁面交角附近出现较小涡环。④局部涡环会随着换热过程而发展。例如，在1 cm工况时，在初始阶段，先在靠近侧边壁面附近形成小的涡环，随之向空腔中部扩展；等到一定程度，涡环发展稳定；随着换热过程继续，涡环消失，仅在壁面附近形成较小的涡流，空腔气体温度随厚度方向变化明显，体现出传热为主的特征。涡环的稳定分布与流场速度、压力及界面换热系数、温度等直接相关。
2.3 空气夹层换热规律
常用Nu数（努塞尔数）和Ra数（瑞利数）描述换热过程中的参数变化[8]，Nu数表征换热系统从流场中吸热的强弱，Ra数表征流场中自然对流的强弱。温度在空气和土体界面上并非完全均匀分布，因此界面上同一时刻对流换热系数、Nu数等并非均匀分布，计算时取界面中点数据。给出土体与空气界面上Nu数随时间的变化如图 5所示。
不同间隙及Ng下，Nu数均随时间降低，在夹层流场启动阶段，Nu数变化较大，表示初始状态流场对流换热较为显著，后续慢慢稳定降低；说明起始阶段（温差最大时刻）换热系统工作效率较高。
3. 结论
通过建立超重力下模型箱土体夹层换热模型，分析了不同厚度及Ng对换热过程的影响规律，研究获得以下结论：
（1）随着超重力Ng值增大，界面降温加速，高Ng值下降温速率大于低Ng值；随着夹层间隙厚度减小，降温速率提升，小间隙下降温速率大于大间隙。
（2）空气与土体界面的换热热量在一定时刻内随间隙增大而减小，随Ng上升而上升，与夹层空腔中涡环发展有关。
（3）超重力状态会使得传热规律发生变化。涡环的稳定分布与流场速度、压力及界面换热系数、温度等直接相关；20 cm间隙时，空腔内呈现典型的RB流动特征，形成两个对称的涡环，超重力使得流动速度变快；10 cm间隙时，不同Ng值下均会出现2对涡环；靠近侧边壁面处涡环较细长，中部涡环长细比小；1 cm间隙时，1g状态有明显涡环成对出现；高Ng值状态下无明显的涡环发展过程，热量传递主要以传热为主。
（4）不同间隙及Ng值下，Nu数均随时间降低，在流场启动阶段，Nu数变化较大；在不同间隙状态下，随着超重力Ng值变大，界面对流换热降低。
-
![]()
图 2 通过滚球法实现的“α-shape”算法
Figure 2. α-shape algorithm accomplished by rolling ball method
![]()
图 4 删除转动角较小的颗粒后进行多次“降噪”
Figure 4. Multiple times of "denoising" after deleting soil particles with relatively small rotation
![]()
图 7 基于“α-shape”算法的剪切带边界的识别
Figure 7. Identification of shear band boundary based on 'α-shape' algorithm
![]()
图 9 5组不同中值粒径试样的剪切带及识别结果
Figure 9. Shear bands and identificated results from five groups of sand samples with different median diameters
![]()
图 10 5组不同中值粒径试样的单一型剪切带及识别结果
Figure 10. Single shear bands and identificated results from five groups of sand samples with different median diameters
表 1 双轴压缩模型的参数取值
Table 1 Parameters for the biaxial compression model
物理意义 数值 物理意义 数值 试样大小 5 cm×10 cm 接触刚度kn/ks 2×107 N/m 颗粒密度 2.6×103 kg/m3 轴向应变速率 5%/min 孔隙率 0.15 目标轴向应变 10% 围压 100 kPa 阻尼系数 0.7 颗粒摩擦系数 0.5 墙壁摩擦系数 0 
下载: 导出CSV
表 2 5组试样PFC2D双轴压缩试验及结果
Table 2 Identificated results of five biaxial loading tests
dmed
/mm识别宽度/mm 相对宽度 网格法宽度/mm 相对
误差/%孔隙率场宽度/mm 相对
误差/%6 8.0 13.3 7.2 -10 11.7 47 7 8.8 12.5 9.1 4 12.6 44 9 9.8 10.9 9.6 -2 13.3 36 12 11.7 9.8 12.6 7 16.2 38 15 13.5 9.0 15.0 11 17.8 32
下载: 导出CSV
-
[1] ALSHIBLI K A, STURE S. Shear band formation in plane strain experiments of sand[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(6): 495–503. doi: 10.1061/(ASCE)1090-0241(2000)126:6(495)
[2] 王冬勇. 基于二阶锥规划有限元法的岩土体稳定性及应变局部化研究[D]. 北京: 北京交通大学, 2020. WANG Dong-yong. Stability and Strain Localization Analysis of Geotechnical Problems Based on Second-Order Cone Programming Optimized Finite Element Method[D]. Beijing: Beijing Jiaotong University, 2020. (in Chinese)
[3] 蒋明镜, 李秀梅. 双轴压缩试验中砂土剪切带形成的离散元模拟分析[J]. 山东大学学报(工学版), 2010, 40(2): 52–58. https://www.cnki.com.cn/Article/CJFDTOTAL-SDGY201002011.htm JIANG Ming-jing, LI Xiu-mei. DEM simulation of the shear band of sands in biaxial tests[J]. Journal of Shandong University (Engineering Science), 2010, 40(2): 52–58. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SDGY201002011.htm
[4] TANG H X, DONG Y, WANG T, et al. Simulation of strain localization with discrete element-Cosserat continuum finite element two scale method for granular materials[J]. Journal of the Mechanics and Physics of Solids, 2019, 122: 450–471. doi: 10.1016/j.jmps.2018.09.029
[5] TIAN J Q, LIU E L, HE C. Shear band analysis of granular materials considering effects of particle shape[J]. Acta Mechanica, 2020, 231(11): 4445–4461. doi: 10.1007/s00707-020-02771-y
-
期刊类型引用(17)
1. 姚贵，张志明，袁勇，禹海涛，吴彦霖. 基于振动台试验的中庭式地铁车站地震响应数值模拟. 山东理工大学学报(自然科学版). 2025(03): 57-63 . 
百度学术
2. 崔春义，许民泽，许成顺，赵经彤，刘海龙，孟坤. 考虑地震需求统计不确定性的地铁车站结构地震易损性分析. 岩土工程学报. 2025(03): 453-462 . 本站查看
3. 张启良. 大跨度无柱平顶结构地铁车站抗震分析. 智能城市. 2025(01): 148-151 . 百度学术
4. 侯本伟，游丹，范世杰，许成顺，钟紫蓝. 基于网络效率的城市轨道交通网络抗震韧性评估. 清华大学学报(自然科学版). 2024(03): 509-519 . 百度学术
5. 丘志杰，龙慧，汪博豪，刘璐瑶. 相邻框架结构与地铁车站的间距比对地铁车站地震响应的影响. 南华大学学报(自然科学版). 2024(04): 26-36 . 百度学术
6. 李维沈，李文婷，徐昊哲. 地下结构抗震性能指标限值的影响因素研究. 自然灾害学报. 2024(06): 178-191 . 百度学术
7. 姚凡夫，杨帆，庄海洋. CFRP布加固的地铁地下车站结构中柱抗震性能研究. 震灾防御技术. 2024(04): 754-762 . 百度学术
8. 苗晗，蒋录珍，安军海，李莎，马晓明. 叠合装配式管廊结构抗震性能水平与评价方法研究. 震灾防御技术. 2023(01): 53-64 . 百度学术
9. 庄海洋，李晟，王伟，陈国兴. 采用不同隔震形式的双层地铁地下车站结构地震反应分析. 振动工程学报. 2023(02): 379-388 . 百度学术
10. 白立广. 地下双层车站盖挖段主体结构施工技术. 建筑机械. 2023(08): 63-67 . 百度学术
11. 王立新，范飞飞，汪珂，李储军，姚崇凯，甘露. 地铁车站不同减震层的减震机理及性能分析. 铁道标准设计. 2022(05): 131-139 . 百度学术
12. 游裕鑫，邵国建，李昂，刘旭. 异跨框架式地铁地下车站结构三维非线性地震响应分析. 河南科学. 2022(04): 610-617 . 百度学术
13. 罗永鸿，张梓鸿，许成顺，李洋. 基于Pushover分析方法的多层地铁车站地震反应研究. 震灾防御技术. 2022(01): 143-153 . 百度学术
14. 张椿民. 地下多层地铁车站结构抗震设计对比分析. 市政技术. 2022(12): 50-58 . 百度学术
15. 范世杰，游丹，侯本伟，许成顺. 基于网络效率的城市轨道交通网络震后性能评估. 防灾减灾工程学报. 2022(06): 1165-1173+1190 . 百度学术
16. 陈文斌，庄海洋，李晟，陈苏. 基于柱顶隔震的3层3跨地铁地下车站结构抗震性能研究. 震灾防御技术. 2021(01): 146-156 . 百度学术
17. 丁录董，徐军林，庄海洋，陈文斌. 预制+现浇装配式地铁地下车站结构地震反应的三维有限元分析. 世界地震工程. 2021(04): 157-166 . 百度学术
其他类型引用(16)
计量
- 文章访问数: 85
- HTML全文浏览量: 21
- PDF下载量: 25
- 被引次数: 33
