冰碛体三维细观结构随机重构及大型三轴试验离散元模拟研究

宋玉香, 钱均益, 刘 勇

(1.石家庄铁道大学 土木工程学院，河北 石家庄 050043；2.石家庄铁道大学 道路与铁道工程安全保障省部共建教育部重点实验室，河北 石家庄 050043)

冰碛体是一种冰川融化后冰水对宽级配碎屑物质进行搬运堆积而成的产物，通常只存在于特定的气候和特定的地区[1]。随着近年来中国高海拔地区的开发，科学经济地解决工程中冰碛体引起的问题成为迫切需要。对冰碛体力学特性研究通常采用室内大型直剪试验[2]、室内大型三轴试验[3]、原位剪切试验[4]并剔除或代替粗颗粒、巨颗粒，但改变粒径组成的试验成果并不能准确表征其性质，由此进行的工程特性分析往往与现场施工结果存在偏差。此外，受限于试验观察技术，试验难以从细观块体尺度上对冰碛体的变形破坏进行深入分析，如冰碛体内部结构损伤及剪切带的形成与演化。

随着近年来不连续介质数值模拟方法的发展[5-7]，目前已有不少学者采用离散元对冰碛体或土石混合体进行了大型剪切试验模拟，数值试验不仅花费代价小、具有可重复性，还可以突破仪器尺寸对试样粒径的限制。徐文杰等[7]利用三维扫描获取真实石块形态进而构建土石混合体的离散元模型，并利用离散元剪切试验研究了土石混合体的细观力学特性，但此类方法仅能代表特定情况下的冰碛体或土石混合体，可重复性和代表性仍不足。张强等[8]基于不规则块石随机生成技术和颗粒离散元模拟技术建立了土石混合体三维模型，并利用三轴试验对不同含石量土石混合体的力学特性和变形破坏机理进行了探讨研究，但颗粒流离散元模拟技术需要标定土体颗粒间的接触参数，标定依据的标准不同就会造成最终模拟结果的差异；同时，颗粒簇组合虽能模拟出不规则块石形状，但并不能准确表征冰碛体内部块石未经风化磨圆的次棱角性。因此，目前仍需要一种更准确的、具有可重复性的大型试验模型模拟方法。

基于块体离散元的方法，并考虑实际块石尺寸形状，建立可用于大尺度工程分析的冰碛体三轴试验模型，利用数值模拟软件3DEC模拟室内大型三轴试验，对试样变形破坏特征及内部剪切带进行分析，为深入分析冰碛体力学特性提供科学方法与依据。

1 冰碛体三维模型生成

1.1 三维随机块石集合生成

冰碛体内部块石形状极为复杂，多为棱角或次棱角状，但在粒度组成、块石形态及排列特征等方面具有分形自相似性特征[9]，其中块石长宽比主要集中在1.4～1.5。在进行三维细观结构随机重构时，需要对块石形状作适当的简化处理，选用不规则凸多面体来近似代替实际复杂形状的块石。冰碛体的三维随机块石集合生成具体步骤为：①根据Medley[10]建议的土/石阈值的取值并结合后续试验试样尺寸，确定块石粒径尺寸为40～100 mm,随机选取区间内尺寸作为块石椭球体基元长轴，以1.4的长宽比确定椭球体基元大小;②确定多面体顶点，根据凸多面体顶点与面之间所满足的空间拓扑关系确定凸多面体外表面集(即对于任一个面而言，除该面所包含的顶点以外，其余顶点均位于该面的同一侧)，连接三角形面集构成凸多面体;③考虑块石空间方位的随机性，在将生成的凸多面体做一个随机转动;④检查生成凸多面体质量并计算多面体体积，将符合实际冰碛体内部块石形态的凸多面体存入三维随机块石集合;⑤重复步骤①～步骤④，直至块石总体积大于当前块石含量所需体积值。

1.2 三维随机块石投放

块石投放是将块石集合中的块石由大到小逐个投放到指定的区域内，并保证预投放块石与已投放块石间不存在相互入侵情形。块石的相互入侵分2种：第1种是一个块石的部分入侵到另一块石内部；第2种是小块石全部入侵到大块石内部(如图1所示)。

图1 块石间入侵类型

冰碛体的三维随机块石投放具体步骤如图2所示。基于以上原理，利用python语言编制三维块石生成程序。基于三维块石生成程序，根据试验尺寸，建立50 cm的立方体冰碛体随机细观结构模型。考虑块石含量及空间分布2个重要指标，分别建立了含石量为5%、10%、15%、20%、25%、30%、35%随机细观结构模型(见图3)，在每一含石量下，又分别生成了3个不同块石空间分布结构模型。

图2 冰碛体三维随机块石投放流程图

图3 不同含石量模型

2 冰碛体大型三轴试验模拟

2.1 冰碛体离散元模型建立

中国新疆、西南地区的冰碛体中块石成分主要由花岗岩、玄武岩、花岗闪长岩、千枚岩、板岩、灰岩中的一种或几种组成；土体成分以粉质黏土为主，其成分及含量很大程度上决定了冰碛体的胶结程度及水敏性。为简化分析，试验模拟中将冰碛体视为特殊的土石混合体，块石视为花岗岩，土视为具有一定胶结的粉质黏土。基于已有研究可知[8]，冰碛体内部块石在加载破坏过程中一般处于弹性受力阶段，因此为更有效地分析冰碛体的工程力学特性，对冰碛体内部块石采用弹性本构模型进行模拟。根据工程地质手册，具体参数见表1。

表1 土石模型参数

2.2 大型三轴试验模拟

采用块石离散元分析软件3DEC进行大型三轴试验模拟，对于冰碛体模型的四周节点施加围压，底部节点固定。通过设置应变率/计算步进行加载，根据《铁路工程土工试验规程》(TB 10102—2010)，选取0.5%/min轴向应变速率进行加载，但通过数值模拟试验发现按此加载速度试样会发生局部破坏，因此通过反复试算，最终确定加载速率为0.1%/min，加载至轴向应变达到20%时结束。

3 冰碛体大型三轴试验模拟分析

3.1 空间分布影响性分析

选取具有代表性的同一含石量不同分布的冰碛体试样的偏应力-应变曲线如图4所示。

图4 同一含石量不同分布的冰碛体试样的偏应力-应变曲线

由图4可以看出，在低含石量(小于20%)情况下，相同含石量不同块石空间分布的冰碛体试样三轴试验模拟获得的偏应力-应变曲线重合，曲线差异率在1%以内；当含石量在20%～25%之间，偏应力-应变曲线基本重合，曲线差异率在1%～3%之间；当含石量大于25%时，内部块石空间分布不同的冰碛体，其偏应力-应变曲线表现出一定的差异性，差异率在5%以上，但总体趋势相同。这表明当含石量超过一定值时，块石分布会对冰碛体的力学特性产生影响。

3.2 围压和块石含量影响性分析

不同围压、不同含石量冰碛体试样的偏应力-应变曲线如图5所示。

图5 不同冰碛体试样的偏应力-应变曲线

由图5可知，冰碛体试样三轴荷载作用下，其偏应力-应变曲线一开始表现出近似弹性变形的特点，随着轴向荷载的不断增大，试样又迅速屈服。随着围压的增大，冰碛体的强度也增大。在高围压作用下(800 kPa)，冰碛体偏应力-应变曲线在屈服后表现为应变硬化型。当围压相同时，冰碛体的含石量越高，其强度也越大且峰值应力对应的应变也逐渐增大，冰碛体的残余强度也随含石量的增加而增大。块石含量在20%～25%之间,冰碛体偏应力-应变曲线发生突变，强度显著增加，且在低围压作用下峰值明显；当块石含量在5%～15%或30%～35%区间内，冰碛体峰值强度及残余强度相近。

从冰碛体试样三轴试验偏应力-应变曲线中，取其峰值偏应力(若没有明显峰值时取轴向应变为15%对应的轴向偏应力值)，将数值与对应围压相加得到冰碛体三轴试验峰值强度。如表2所示。

表2 冰碛体三轴试验峰值强度 kPa

根据表2绘制同一围压下冰碛体三轴试验峰值应力随含石量变化图，并记录不同围压下冰碛体峰值应力随含石量的增长率，如图6所示。从图6可知，冰碛体峰值应力随含石量增大不断增大，含石量35%的冰碛体较含石量5%的冰碛体峰值应力提升超过100%，且围压越大增幅越大，但不同围压间增长率相差不大。由图6可见，冰碛体峰值应力在含石量10%～30%之间发生显著增长，增长率超过10%，其中含石量15%～25%峰值应力增长率最大，峰值应力变化剧烈。

图6 冰碛体三轴试验峰值应力及其增长率随含石量变化图

3.3 强度特性分析

依据表2冰碛体峰值应力绘制试样在不同围压条件下的极限莫尔应力圆和强度包络线。通过强度包络线拟合摩尔库仑强度方程，求得各含石量试样对应的黏聚力、内摩擦角如表3所示。

表3 各含石量试样抗剪强度参数

由表3可知，冰碛体黏聚力和内摩擦角随着含石量的增大而增大。含石量的变化对黏聚力影响更显著，含石量35%冰碛体较含石量5%冰碛体，黏聚力提高了68.21 kPa，内摩擦角提高了8.63°。依据变化趋势可分为3段：当含石量小于15%时，黏聚力和内摩擦角增长速率较缓；当在含石量10%～30%之间时，黏聚力和内摩擦角剧烈变化，近似线性增长；其中，黏聚力在冰碛体含石量15%～20%之间增长率最大，为38.95%，内摩擦角冰碛体含石量20%～25%之间增长率最大，为8%；当冰碛体含石量在大于30%阶段，黏聚力和内摩擦角变化幅度不大，增速变缓，力学特性相近。

当荷载作用较小时，冰碛体中土体未达到屈服状态，土体与块石都处于弹性变形阶段，但有部分土-石接触面可能产生应力集中进而滑动或者破坏，所以冰碛体在屈服前偏应力-应变曲线并不是完全的线性增长，而是一种近似线性的形态。当偏应力-应变曲线达到峰值后，土体破坏，但破坏后土体及内部块石仍有残余强度，因此表现出应变软化特性。而当试样处于高围压状态(800 kPa)时，土颗粒或块石咬合紧密，而块石仍处于弹性阶段，使得应力-应变曲线会在达到峰值后变为应变硬化型或直接变为持续硬化型。

当含石量较低时(小于15%)，冰碛体的强度特性主要取决于土体，由于内部块石在土体中处于“悬浮”状态，因此其力学特性近似纯土体，随着块石含量增大，抗剪强度参数缓慢增大；当含石量在15%～30%之间增长时，冰碛体力学特性受内部结构变化的影响，力学特性发生剧烈变化，抗剪强度不断增大。当含石量大于30%时，块石之间在荷载作用下接触增多，容易形成“咬合”作用，块石开始逐渐成为冰碛体的主要骨架，对内部土颗粒约束作用增强，抗剪强度参数缓慢增长。

3.4 变形破坏特征分析

不同冰碛体试样沿中心进行竖向切片后，破坏后试样内部单元的位移云图如图7所示。

图7 不同冰碛体试样的内部剖切位移云图

由图7可知，在围压作用下，冰碛体试样破坏后侧向均表现出了不均匀的鼓胀变形，与其他室内试验研究结果[4]相吻合。冰碛体试样破坏后形成的剪切带并非一个规则的条带，而是一个曲折的条带，有明显绕石现象，其位置基本上位于接近于试样的中部，且形态上近似一个非对称的扁平X形分布。由图7(a)、图7(b)可以看出，对于同一冰碛体试样，在不同围压下，试样破坏后的形态是不同的，围压大小影响冰碛体的变形破坏形态，由此也体现出了冰碛体的非均质、非均匀特性；由图7(a)、图7(c)可以看出，在相同围压下，不同含石量的冰碛体试样破坏后的形态也是不同的；图7(a)、图7(d)为在相同含石量不同空间分布冰碛体试样同围压下三轴试验破坏后内部剖切位移云图，可以发现二者形态不完全相同，这体现出了冰碛体变形破坏规律的复杂性。综上所述，冰碛体在三轴加载下其内部剪切带的大小和分布形态受围压、块石含量、空间分布共同影响。

冰碛体在三轴荷载作用下发生剪切破坏时大多表现出较明显的剪胀现象，并且随着围压的减小剪胀现象越明显。这是由于低荷载条件下组成冰碛体的岩土体强度较颗粒之间的咬合力和摩阻力大，尤其是内部块石强度，发生剪切作用时，土颗粒及块石沿着相互之间的接触面滑移，进而引起试样体积的增大，最终表现出剪胀特点。而随着块石含量和围压的增大，块石之间在荷载作用下接触增多，容易形成“咬合”作用，块石开始逐渐成为冰碛体的主要骨架，对内部土体产生一定的约束作用，内部单元运动的整体性变强，侧面鼓胀变形量变小。

4 结论

运用计算机三维随机离散块体建立了不同空间分布和块石含量的冰碛体结构模型，开展了不同围压下冰碛体的大型离散元三轴试验模拟，对冰碛体的力学特性与变形破坏规律进行了分析，得出如下结论：

(1)在低荷载作用下冰碛体试样的应力-应变曲线近似直线状，表现出弹性变形的特点；随着荷载的进一步增大，试样迅速屈服。冰碛体强度随围压的增大而增强；围压大于800 kPa时，冰碛体试样在破坏后表现出应变硬化。

(2)含石量对冰碛体抗剪强度影响显著，随含石量的增加，冰碛体的强度不断增大，表现在黏聚力和内摩擦角同时增大，含石量35%冰碛体较含石量5%冰碛体，黏聚力提高了68.21 kPa，内摩擦角提高了8.63°；其中含石量15%～30%阶段，冰碛体处于“块石悬浮”向“土石密实”转化阶段，抗剪强度参数急剧增大，块石对土颗粒约束增强，其余阶段黏聚力和内摩擦角缓慢增大。冰碛体在含石量大于25%时，内部块石空间分布会对其力学特性产生影响，差异率在5%以上，但总体趋势相同。

(3)在围压加载下，冰碛体试样表现出剪胀现象，且围压越低剪胀现象越明显；破坏后所形成的剪切带表现出一定的绕石现象，形态上呈非对称的X形分布，且受块石含量、空间分布、围压大小共同影响。