基于PFC^2D的冻融循环作用下冰碛土微观损伤研究

Microscopic Damage Evolution of Moraine Soils under Freeze-Thaw Cycles Based on PFC^2D Simulation

摘要

为探究冻融循环作用下冰碛土性质劣化的微观损伤机制，基于离散元理论提出一种通过水颗粒膨胀实现土体冻融损伤模拟的方法.利用颗粒流软件PFC^2D模拟三轴压缩试验，结合室内试验结果对比分析，在模拟冰碛土的力学性质变化方面显示出高度的准确性和可靠性，揭示了冻融冰碛土受载时微裂隙、位移场和力链场的演化过程与破裂特征.结果表明：（1）试样冻融过程中微裂隙呈由四周产生并逐渐向中间扩展的“累计演化”趋势，其中张拉微裂隙占据主导地位，在冻融前期（2~5次）颗粒间以水平挤压为主从而大量发育偏90°倾向张拉微裂隙；（2）冻融作用引起的冰碛土性质劣化在冻融前期（2~5次）尤为明显，粘聚力c随冻融次数N的增加呈负指数函数递减规律，而内摩擦角φ呈小幅波动态势；（3）试样受载过程中剪切微裂隙占据主导地位，微裂隙发育呈“慢→陡→缓”趋势演化，根据其裂隙演化特点，将加载过程的应力-应变曲线划分为4个变形阶段；（4）冻融后试样受载时减速斜率转换点B移动到峰值应力点C之前，说明B点在微裂隙扩展-贯通-形成破坏过程中可以作为“前兆特征”；冻融20次试样受载时破坏程度更剧烈且形成明显剪切破坏带.

Abstract

In order to investigate the microscopic damage mechanism of the degradation of the properties of freezing-thaw (F-T) damaged moraine soils, a method of simulating F-T damage of soils through water particle expansion is proposed based on the discrete element theory. Using the particle flow software PFC^2D to simulate the triaxial compression test, combined with the comparative analysis of the test results, this method is accurate and reliable in modeling the changes in mechanical properties of moraine soils and reveals the evolution of microcrack; displacement field; force chain field and rupture characteristics of F-T moraine soils during the loading process. The results show follows (1) The microcracks in the F-T process show a trend of “cumulative evolution” that arises from the surrounding area and gradually expands to the middle, and the tensile microcracks are dominant; at 2-5 times of F-T processes, horizontal compression between particles dominated and a large number of tension microcracks inclined at 90° were developed. (2) The deterioration of moraine properties caused by F-T is particularly obvious in the early freeze-thaw period (2-5 times), the cohesion c decreases as a negative exponential function with the number of F-T cycles, while the internal friction angle φ shows a small fluctuation. (3) The specimen loaded process is dominated by shear cleavage, with the trend of “first slow, then steep, and finally slow” evolution, and the stress-strain curve is divided into four deformation stages according to the evolution characteristics. (4) When the sample is loaded after freeze-thaw, the transition point B of deceleration slope moves before the peak stress point C, indicating that point B can be used as a “precursor feature” in the process of microcrack expansion-through-formation failure; The specimen with 20 F-T cycles were more severely damaged when loaded and formed distinct shear zones.

Graphical abstract

关键词

冰碛土 / 离散元 / PFC^2D / 冻融循环 / 微观损伤 / 破裂演化 / 三轴压缩 / 工程地质学.

Key words

moraine soil / discrete element method / PFC^2D / freezing-thaw cycle / microscopic damage / rupture evolution / triaxial compression / engineering geology

冰碛土是第四纪冰川行进并融化过程中在冰碛地貌中形成的一种特殊土壤，具有无分选性、级配不良和混杂堆积的特点（Rossi et al.， 2019； Wang et al.， 2024）.由于季节轮换和昼夜交替引起温度急剧变化，土体内水分在冻结和融化状态之间周期性转化，这直接导致土体孔隙结构变化、强度减弱及土体力学性质劣化（彭建兵等，2023；陈志敏等， 2024），一旦土体承载能力达到极限，可能引起山体滑坡在内的各种地质灾害（Wang et al.， 2023），如冰湖溃决和冰川泥石流等（Deng et al.， 2017； Veettil and Kamp， 2021），对基础设施和居民的生命财产安全构成严重威胁.目前，已经有众多学者从多角度对冰碛土这类特殊岩土体开展了广泛研究，对冰碛土基本物理力学性质已经取得基本认识.Zhou et al. （2019）对不同含水率冰碛土进行排水剪切试验，发现含水率主要影响冰碛土剪胀性.Wang et al. （2024）研究了含冰率与温度对冰碛土力学性质的影响，通过热力学三轴试验，其研究指出冰碛土粘聚力和内摩擦角随温度升高而减小；随着冰含量的增加，冰碛土的峰值强度增加.蒋婷婷等（2024）对冰碛土开展不固结不排水三轴试验，综合考虑了冻融循环次数、初始含水率和围压对冰碛土力学参数的影响规律.江强强等（2020）发现含石量较大时（20%~80%），块石形成骨架结构，强度由块石和细粒土共同决定.

在细观层面，何迎红和屈智炯（1990）、屈智炯等（1992）在20世纪90年代采用扫描电镜研究了冰碛土颗粒的微观结构并进行了力学性状的研究，发现其应力-应变曲线特征与围压有密切关系.吕士展等（2014）对冰碛土原状样进行了CT扫描，获得了冰碛土内部真实细观组构图像与孔隙特征.Qiu et al. （2023）利用声发射（AE）分析了冻融循环后冰碛土试样的单轴压缩变形响应，指出冰碛土的变形曲线表现为4个变形阶段.随着非连续介质数值模拟技术的显著进步（周先齐等，2007），众多学者现已采纳离散元方法深入探索岩土体冻融循环过程以及冰碛体剪切过程的内部动态变化.这类数值实验的明显优势体现在其能够细致揭示冻融循环如何细微地影响岩土结构，尤其是观察到裂纹的形成、颗粒的重新排列等微观级别的结构变动，这些微观变化通过传统的实验手段往往是无法直观捕捉到的.例如孙渊等（2023）利用PFC^3D软件，揭示了黄土冻融损伤机理与根系固土作用；徐文杰等（2016）利用三维扫描技术构建土石混合体数字模型，并开展三维数值直剪试验研究了土石混合体的力学性质；宋勇军等（2023）基于离散元提出水冰颗粒相变耦合膨胀方法，实现了岩石冻融过程的裂隙演化过程.不同于岩石，土颗粒排列方式相对松散，其孔隙率更大且孔隙结构更为复杂（蒋德旺等， 2019）.冻融循环作用会导致土体内部的细微孔隙逐步合并，形成大中孔隙（Xie et al.， 2015； Han et al.， 2023； Tao et al.， 2023），这些孔隙可以形成毛细粘聚力并容纳水分，影响土体的物理强度指标与水分的运移（Leuther and Schlüter， 2021）.目前针对冰碛土在冻融循环过程中微裂隙场、位移场和力链场的演化过程还罕有报告.

综合上述考虑，本研究基于PFC^2D软件，通过水颗粒膨胀法对冰碛土试样进行冻融循环及三轴压缩模拟，模型考虑了毛细粘聚力并结合室内三轴试验结果进一步校准模型的合理性.不仅从宏观层面评估冰碛土的力学性能劣化特征，更从微观上洞察到冻融过程中土体颗粒的运动、重新排列和破裂等行为，揭示冰碛土在冻融作用下的微观结构变化、力学性能劣化规律以及损伤累积过程，为冰碛土冻融循环下的性质劣化提供新的视角.

1 试验材料及方案

1.1　试验材料

青藏高原山南市努确沟流域内共分布有6处冰川泥石流沟、8处冰湖，广泛分布有季节性冻融冰碛土.试验物料取自埋深30 cm的冰碛堤区域现场.为了减小对冰碛土的扰动，参照吕士展等（2014）的做法，人工小心挖出原状土柱后将其装入边长为20 cm的立方体盒，然后用保鲜膜包裹带回实验室.测得其湿密度后将土样在烘箱中烘干24 h，对其中一部分土样进行筛分试验，试验步骤见图1.平均颗粒级配曲线见图2，不均匀系数Cu=2.2，曲率系数Cc=0.94，属级配不均匀砂土.剩余土样过2 mm直径筛孔以备后用，由于土样粗粒部分（粒径大于2 mm）质量占比24.1%，低于粗粒控制土体力学性质的临界值40%（Kuenza et al.， 2004）.因此，移除粒径大于2 mm的粗粒对冰碛土的物理力学性质影响很小，可有效反映其本身的特性.根据《GB/T50123-2019土工试验方法标准》，试验测得的土样含水率、比重等基本物理性质指标如表1所示.采用Bruker D8 Advance X射线粉晶衍射仪分析，角度范围为5°~40°，扫描速度为5°/min，得到的XRD 衍射谱图如图3所示，冰碛土矿物成分为石英、钠长石、透长石、角闪石、白云母和少许绿泥石.

将过筛后的土样配水并采用分层压实法制得18个标准圆柱体试样（底面半径为25 mm、高度100 mm），含水率均为10.71%.为了防止水分转移，盖上保鲜膜密封24 h.

1.2　试验方案

对18个土样进行0、2、5、10、15、20次冻融循环试验，结合山南市近5 a夏季与冬季月最高（最低）温度的平均值 （图4），确定每一个循环周期持续24 h（汤明高等， 2022），即：冰碛土在-16 ℃的温度下冻结12 h，在28 ℃温度下静置12 h进行解冻，完成1次冻融循环后将试样放入养护箱中静置12 h，湿度设置为60%（与野外监测资料一致），从而实现补水的过程，冻融过程温度曲线如图5.考虑到冰碛堤冰碛土的埋深浅，同时冰碛土在实际冻融环境中经历周期性的性质变化，对不同冻融循环次数后的试样进行三级围压（2、10、20 kPa）固结不排水试验（CU），三轴压缩时轴向加载速度设置为0.5 mm/min，具体试验方案如表2所示.

2 冻融损伤方法与模型建立

2.1　模拟方法的实现

冻融循环造成土体劣化的本质在于其对土体微观结构的影响.当土体经历冻融循环时，水分在土体中的冻结和融化过程导致体积的反复膨胀与收缩，这导致土粒重新排列（Konrad， 1989； 肖东辉等， 2014； Wang et al.， 2020），造成孔隙重分布并发育微裂隙，由此削弱颗粒之间的粘聚力.重复的冻融循环会不断累积这种结构性的损伤，导致土体强度下降、渗透性变化等一系列劣化现象，最终影响土体的工程性能.

结合上述土体经历冻融后的微观破坏情况分析，提出一种基于颗粒离散元模型的水冰相变膨胀颗粒方法.利用PFC^2D软件来模拟冰碛土的冻融循环过程，并采用以下假设：

（1）将冰碛土试样简化为由土颗粒骨架与孔隙水颗粒组成.

（2）不考虑土体中的非冻结水（谭龙等， 2015），试样内的水颗粒总是同时冻结与解冻.

（3）冻融过程中不考虑孔隙水渗流.

每次冻结过程中，遍历所有水颗粒膨胀其体积至原体积的1.09倍（Zhu et al.， 2021），使得水颗粒与周围颗粒重叠，从而产生膨胀力的作用.膨胀力由颗粒接触处产生，向土骨架内部传递，当力的大小超过颗粒之间的粘结强度时，粘结键断裂并生成裂隙，造成冻融损伤.随后冰颗粒融化，水颗粒体积恢复到原始体积，同时提高土颗粒间的粘结强度以模拟毛细粘聚力的形成.最后再次遍历所有水颗粒，将其体积膨胀至当前体积的1.01倍以模拟试验中的补水过程并作为下一次冻融的初始水体积，以上完成一次冻融循环.具体冻融过程水颗粒体积变化与裂隙的形成过程如图6.

2.2　PB粘结模型与数值模型建立

本试验的的冰碛土由筛分后细颗粒胶结形成，具有较强结构性（冯俊德等， 2008）.为了模拟颗粒间的粘结特性（张永双等， 2007），同徐安花等（2025）一致，颗粒间的接触模型采用平行粘结接触模型（parallel bonded contact model）.

平行粘结接触模型由Potyondy and Cundall（2004）提出，是一种基于线性模型的模型，可在颗粒与颗粒之间或颗粒与墙之间接触，可以同时传递力与力矩，建立后因平行粘结键的存在，接触处的相对运动会致使粘结材料内产生力和力矩.未粘结或者粘结破坏后，该模型退化为线性接触模型.粘结破坏模式取决于抗剪强度以及抗拉强度，当拉力超过抗拉强度时则产生拉伸破坏；当剪力超过抗剪强度时，则产生剪切破坏.

模型设置为高100 mm、宽度50 mm，高纵比为2︰1的标准圆柱试样.水颗粒最小粒径为0.2 mm，最大粒径为0.25 mm，土颗粒级配特征与试验土样保持一致.采用分层压缩法分5层逐层夯实压缩颗粒至目标孔隙率0.183，最终生成冰碛土颗粒9 443个.土体颗粒密度1 700 kg/m³，水颗粒密度1 000 kg/m³，含水率用水颗粒的面积占比表征，设定为10.7%，图7为生成的模型.

2.3　细观参数标定

初始状态试样通过进行无侧限单轴压缩模拟试验来进行验证，结合室内单轴压缩实验结果，对未经冻融循环影响的试样进行应力-应变曲线的参数校准.在进行单轴及三轴压缩数值模拟时，鉴于实际实验的加载过程，解冻后试样的强度主要由土颗粒骨架支撑.因此，在单轴压缩数值模拟过程中，需要剔除模型中的孔隙水颗粒.

平行粘结接接触模型无法很好地体现试样受载时的压密阶段，基于“试错法”连续迭代调整微观参数（吴禄源等， 2023），使得单轴抗压强度和弹性模量的数值结果与宏观力学参数较为吻合，并且试样的最终破坏形态相似（图8）.将水-水颗粒间粘结和水-土颗粒间粘结的强度设置足够大，以避免微裂隙在土骨架以外的地方产生从而干扰结果.初始状态下数值模型的接触参数列于表3中.

3 模拟验证与结果分析

3.1　模拟结果与室内结果对比分析

为验证上述数值模拟结果的可靠性，将不同冻融次数后的模拟试样进行三轴压缩模拟.基于伺服原理控制左右墙体移动，模拟围压的作用，其取值与试验围压保持一致，即σ₃取值2、5、10 kPa.随后控制上下墙体向模型中心移动，模拟加压过程，应变速率与试验一致，设定为0.5 mm/min.以每次加载时轴向峰值应力σ为纵轴，σ₃为横轴，绘制最佳关系曲线，并按公式（1）与公式（2）求得粘聚力c与内摩擦角φ：

c=\frac{{\sigma }_c\left(\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)}{\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\phi }，（1）

\phi =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\frac{m-\mathrm{1}}{m+\mathrm{1}}，（2）

式中：σ_c 为最佳关系曲线在纵坐标上的截距，单位为kPa；m为最佳关系曲线的斜率.

试验与模拟的最佳关系曲线的拟合效果较好（图9），轴向应力均随着围压的增大而变大，这说明围压的围箍作用有助于提高土体的承载能力.室内试验与数值模拟的最佳关系曲线较为贴近，认为数值模拟可以还原真实室内试验.

试验与模拟得到结果对比见图10，两种方法测得的同一指标与冻融次数N的关系曲线规律相似，即随着冻融次数N增加，粘聚力c呈负指数函数递减，而内摩擦角呈小幅波动态势，这一规律同已有相关研究一致（刘耀辉等， 2023），进一步验证了数值模拟的结果.

对于结果的差异性做以下说明：

（1）天然冰碛土作为一种松散颗粒材料，往往包含不同粒径和形状的颗粒，建模时需要使用异形颗粒来构建形状更复杂的颗粒团.本文试验所采用的冰碛土分布在冰碛堤上，常年受到强烈冰缘寒冻的风化作用，磨圆度较好，形状差异性小，故采用圆形颗粒建模.然而，模拟未能充分考虑土颗粒形状差异的影响，这导致试验结果与模拟结果有些许误差.

（2）模拟过程忽略了水分渗流对土体的软化作用，导致模拟结果与试验结果出现差异.

（3）在实际的冻融循环过程中，即便是在土体结构遭到破坏后，土颗粒间仍然由于范德华力、分子力以及由毛细孔作用产生的基质吸力而相互作用.而平行粘结接触模型未能充分考虑这一问题，由此导致结果误差，可进一步研究适用于冰碛土颗粒的本构模型.

3.2　不同冻融次数下微观损伤特性分析

冰碛土试样在不同冻融循环次数下微裂隙的演化过程见图11，图中绿色部分代表张拉微裂隙，蓝色部分代表剪切微裂隙.可以看到微裂隙逐步从试样周围开始发育并向内部扩展，这是由于试样四周的水颗粒受冰碛土土颗粒的约束相较中心位置更少，并且其直接接触四周刚性墙体边界，更容易在冻结体积膨胀过程中使土骨架断裂损伤，由此在四周形成大范围微裂隙.

随着冻融循环次数不断增加，尽管土颗粒间的法向粘结强度逐渐增大，但中间土骨架逐渐抵抗不住水颗粒发生水冰相变产生的冻胀应力的作用而发生损伤，整体表观为微裂隙从四周向内部逐渐扩展的“累积演化”趋势.

为了验证这一结果，对冻融5次后试样距离中心不同位置处取样进行电镜扫描.取样位置与扫描结果如图12所示：远离中心区域（a位置）孔隙结构更为松散，孔隙直径更大（134~310 μm），在孔隙周围可见明显的微裂隙发育；中心区域（b位置）的孔隙分布均匀，直径在60 μm左右，直径更小，微裂隙发育不明显.冻胀作用会导致孔隙的扩张，进一步引起土体骨架微裂隙的发育（Liu et al.， 2020），更大的孔隙直径对应更强的冻胀作用历史，由此证明冻融初期（2~5次）土体微裂隙主要在四周发育.

图13对不同冻融循环次数试样的微裂隙进行统计.可以看出，在任一冻融循环次数下，张拉微裂隙数量都远大于剪切裂隙数量，占据主导地位.这是因为孔隙水颗粒体积发生膨胀并产生冻胀应力，从而对周围的土骨架产生挤压，当冻结应力超过冰碛土土骨架的抗拉强度时，土骨架断裂产生张拉微裂隙；其次，张拉微裂隙占比在冻融次数2~5时涨幅较大，说明张拉作用在冻融初期尤为强烈.

为了探究冻融过程中微裂隙发育的形态，统计了微裂隙的倾向和数量，以微裂隙数量为径向、倾向为圆周方向绘制裂隙玫瑰花图（图14）.发现在冻融2~5次时，张拉微裂隙在偏竖直方向发育明显，这是因为在冻融初期颗粒间以水平挤压为主，粘结破坏后产生了垂直于粘结的竖向裂隙.后续冻融过程中张拉微裂隙沿四周均匀发展.剪切微裂隙在整个冻融过程无明显发育倾向，表现为无序性.

图15绘制不同冻融次数下试样的力链场与位移场的演化过程，力链指受载时力在颗粒间的传递路径，力链粗细代表力的大小.由力链场云图显示（图15b），冻融产生的膨胀力先在试样四周集中，力链方向以水平向为主，随着冻融循环的进行，力链向试样内部扩散，传递方向改变为向四周扩散.

从位移云图来看（图15a），在冻融次数为2~5次时，试样中仅有部分位置发生了较小位移，位移主要集中在表面，最大值为1.95×10^-4 m，随冻融循环次数的不断增加，位移整体增大，当冻融循环次数为20次时，试样达到最大位移量为4.73×10^-4 m，反映出冻融作用对冰碛土试样损伤进一步加重，冻融后期四周位移大于内部位移，上端部位移大于下端部位移.

3.3　不同冻融次数后冰碛土受载损伤分析

仅依靠室内标准三轴实验很难了解到加载时微裂隙的“起裂-发育”演化过程，图16绘制了当围压为10 kPa时不同冻融循环次数试样三轴压缩模拟过程中应力应变曲线与微裂隙的关系曲线.

与冻融循环模拟过程不同，整个加载过程剪切微裂隙数目均高于张拉微裂隙，说明加载过程土体的失效破坏以剪切作用为主.整个加载过程中，微裂隙发育呈“慢-陡-缓”趋势.从局部放大图可知，在裂隙发育初期，剪切微裂隙提前出现，张拉微裂隙出现时间略微滞后，整体缓慢发展.随后剪切微裂隙增长速率突然变快，将此增速斜率转换点定义为A点，做为微裂隙的“起裂点”.随着轴向应力的不断增加，微裂隙数量快速上升后逐渐变缓，这是由于裂隙迅速发育贯通，土体强度即将失效，分别将剪切微裂隙迅速增长段曲线段和增速变缓曲线段做线性拟合，两直线交点垂直投影到应力-应变曲线上的B点定义为减速斜率转换点.将试样的峰值轴向应力点定为C点，应变为3.5%时定义为D点.

综上所述，冻融后试样应力-应变与微裂隙发育曲线可大致分为以下4个阶段：微裂隙产生并缓慢发育OA阶段（Ⅰ）、微裂隙起裂并迅速发育AB阶段（Ⅱ）、微裂隙贯通并形成破坏带至试样破坏BC阶段（Ⅲ）、峰值应力后期微裂隙缓慢扩展CD阶段（Ⅳ）.

试样在进行三轴压缩时，将起裂时对应的轴向应力称为“起裂应力阈值”，其中剪切微裂隙起裂应力阈值为σ_N^s（应变为剪切起裂应变ε^s），张拉微裂隙起裂应力阈值为σ_N^t，峰值轴向应力为σ_f.表4统计了上述各指标，随冻融循环次数N的变大，这些指标均逐渐降低，在冻融前期（2~5次）降低幅度较大.以上结果表明冰碛土土颗粒骨架经反复冻融循环作用后胶结变弱，受载时起裂应力变低，从而导致强度性质发生劣化，尤其在前5次冻融过程中强度性质劣化明显.

3.4　试样三轴受载多场演化过程

以冻融0次和冻融20次的两个试样作为研究对象，对比在σ₃=10 kPa条件下冻融前后试样在单次受载破坏演化的全过程.图17和图18分别记录为冻融循环0次和20次试样在三轴压缩模拟时受载破坏演化全过程（包括力链场、位移场、及微裂隙演化图），结合并参考图16a和图16f来分析试样受载破坏演化的全过程.

Ⅰ阶段（A点前阶段）：该阶段为微裂隙产生并缓慢发育阶段，历时很短，由于冻融循环的作用冻融20次试样初始损伤大于冻融0次的试样（图17a，图18a），其剪切起裂应力（28.6 kPa）低于未冻融试样的起裂应力（44.4 kPa）.由位移场得知（图17b，图18b），冻融0次与20次试样都处于被压缩状态，但冻融20次试的位移云图呈上下两端内凹的形态，这是由于冻融20次试样经冻融循环作用，颗粒重新排列后造成颗粒间空隙变大，加载时中间颗粒往两侧边界移动，造成左右两侧颗粒的位移更小，中间部分的颗粒位移相对较大；从力链云图（图17c，图18c）可知，冻融0次试样受载时力链在上下边界产生，并竖直向试样中心传递，分布较均匀；而冻融20次试样的力链产生在上下边界，并未向试样中心传递，可见冻融损伤导致颗粒间力的传递效应减弱.

Ⅱ阶段（AB阶段）：该阶段为微裂隙起裂并迅速发育阶段，在B点时候，冻融0次试样（σ=98 kPa）与20次试样（σ=86.5 kPa）微裂隙均迅速增多，由上一阶段上下端部产生的较大位移区域逐渐向内部扩展，由于加载的压密作用，位移场两端内凹现象减弱.冻融0次试样上下端部的最大位移量为5.21×10^-5~6.41×10^-5 m，最大力链值为128.65 N；冻融20次试样上下端部的最大位移量为4.82×10^-5~4.91×10^-5 m，最大力链值为133.35 N，结果表明，冻融后试样在轴向应力较低时（σ=86.5 kPa），较小位移量就可造成较大的力链，进一步导致裂隙的快速发育.

Ⅲ阶段（BC阶段）：该阶段为微裂隙贯通直至试样破坏阶段.由于B点与C点的时间间隔很短，两点处微裂隙演化图、位移场、力链场无明显差异.在冻融循环次数为0次时，C点位于B点之前（图16a），冻融之后，C点均位于B点之后（图16b~16f）.这说明试样在未冻融情况下，首先达到轴向应力峰值，随后微裂隙贯通形成破坏带；试样在经受冻融作用后，微裂隙先贯通并形成破坏带随后试样达到峰值应力.这表明，冻融作用导致了试样内部微观结构的变化，使得微裂隙在冻融循环中更容易发育和贯通.

Ⅳ阶段（CD阶段）：该阶段为峰值应力后期微裂隙缓慢扩展阶段，次生微裂隙不断产生，试样已失去承载能力.冻融0次与冻融20次试样云图演变趋势类似，但破坏形式有所不同（表5）.由于冻融作用导致了试样内部微观结构的变化，使得微裂隙在冻融循环中更容易发育和贯通.在未经冻融情况下，试样达到轴向应力峰值时，微裂隙的形成和扩展相对受限，因此破坏主要集中在试样表面或颗粒间的接触面.而经过冻融作用后，试样性质劣化明显，导致微裂隙的形成和扩展更加频繁和显著，微裂隙更容易贯通并形成破坏带.

综上所述：冻融循环0次与20次试样微裂隙、位移场、力链场云图的演化规律类似，但经过冻融循环后的试样由于颗粒的重新排列，竖直方向受载时表现出颗粒向水平两侧“逃逸”，力链不能快速传递到试样中心.经冻融循环后的试样在各阶段的微裂隙发育演进过程均明显快于未冻融试样，导致最终破坏形式不同，可将B点作为试样破坏前兆.

表5对比了数值模拟和室内试验在10 kPa围压条件下试样受载破坏的形态（蓝色为完好颗粒，红色为破坏颗粒）.发现在未经冻融情况下，破坏主要集中在试样表面.而冻融过后的试样均贯通形成明显剪切带，由此说明冻融循环作用对试样内部结构的“扰动”十分显著（邱恩喜等， 2024），其微裂隙更容易形成并扩展，导致试样在受载后破坏明显，进一步说明了模拟的合理性.

PFC数值模拟结果一方面揭示了冰碛土在冻融循环下的强度劣化特性，另一方面重现了冻融循环下试样受载微裂隙发育演化全过程.研究结果可为青藏高原冰湖溃决等链式灾害防控研究和工程建设设计提供参数指导和理论支撑.

4 结论

基于离散元软件PFC^2D，本文提出的水颗粒膨胀法实现了土的冻融损伤过程模拟，考虑了毛细粘聚力的作用并揭示了冰碛土在冻融循环作用下的强度劣化过程，主要得出以下结论：

（1）随着冻融循环次数的增加，粘聚力呈负指数函数减少，内摩擦角呈小幅波动态势；增大围压可以提高冰碛土土体的承载能力.

（2）冻融循环过程中张拉微裂隙占据主导地位，微裂隙呈从四周向内部逐渐扩展的“累积演化”趋势，在冻融循环2~5次时张拉作用尤其强烈，这一期间发育的张拉微裂隙的倾向集中于90°；膨胀力链首先在试样四周表面聚集，随后向试样内部传递，其传递方向先是水平传递随后向四周扩散；试样整体位移随着冻融次数的增大而增大，冻融15~20次时外部位移大于中心位移，上端部位移大于下端部位移.

（3）数值模拟三轴加载时微裂隙演化特征呈“慢-陡-缓”的趋势演化；随着冻融循环次数的增加，试样剪切起裂应力σ_N^s、剪切起裂应变ε^s、峰值应力σ_f在前5次冻融循环过程中快速降低，之后缓慢减小；冻融后的试样，微裂隙减速斜率转换点B移动到峰值应力点C之前，说明在微裂隙“扩展-贯通-形成破坏带”这一过程中，点B可作为试样破坏的“前兆特征”.

（4）未冻融试样受载时，微裂隙的形成和扩展相对受限，破坏集中在试样表面.冻融后试样由于内部微观结构变化，受载后微裂隙更快贯通并形成明显剪切带.

参考文献

[1] Chen, Z. M., Liu, Y. H., Li, N., 2024. Research on the Shear Characteristics of Moraine under Freeze-Thaw Action and the Surrounding Rock Pressure. Chinese Journal of Underground Space and Engineering, 20(2): 460-470 (in Chinese with English abstract).

[2] Deng, M. F., Chen, N. S., Liu, M., 2017. Meteorological Factors Driving Glacial till Variation and the Associated Periglacial Debris Flows in Tianmo Valley, South-Eastern Tibetan Plateau. Natural Hazards and Earth System Sciences, 17(3): 345-356. https://doi.org/10.5194/nhess-17-345-2017

[3] Feng, J. D., Li, J. G., Wang, R., et al., 2008. Large Scale Direct Shear Test on Strength Behavior of Railway Moraine Soils in Yunnan. Rock and Soil Mechanics, 29(12): 3205-3210 (in Chinese with English abstract).

[4] Han, M. X., Peng, W., Ma, B., et al., 2023. Micro-Composition Evolution of the Undisturbed Saline Soil Undergoing Different Freeze-Thaw Cycles. Cold Regions Science and Technology, 210: 103825. https://doi.org/10.1016/j.coldregions.2023.103825

[5] He, Y. H., Qu, Z.J.,1990. Mechanical Property and Microstructure Research of Glacial Till. Advanced Engineering Sciences, 22(5): 57-62 (in Chinese).

[6] Jiang, D. W., Cui, P., Wang, J., et al., 2019. Experimental Study on the Effect of Shear Strength of Moraine Soil with Fine Grain Content. Journal of Glaciology and Geocryology, 41(1): 129-139 (in Chinese with English abstract).

[7] Jiang, Q. Q., Xu, Y. Q., Wang, H., 2020. Research on Shear Deformation Characteristics of Soilrock Mixtures under Different Stone Contents. Journal of Engineering Geology, 28(5): 951-958 (in Chinese with English abstract).

[8] Jiang, T. T., Pan, H. L., Ai, Y. F., et al., 2024. Effect of Freeze-Thaw Cycles and Water Content on the Mechanical Properties of Moraine Soil. Bulletin of Geological Science and Technology, 43(2): 238-252 (in Chinese with English abstract).

[9] Konrad, J. M., 1989. Physical Processes during Freeze-Thaw Cycles in Clayey Silts. Cold Regions Science and Technology, 16(3): 291-303. https://doi.org/10.1016/0165-232X(89)90029-3

[10] Kuenza, K., Towhata, I., Orense, R. P., et al., 2004. Undrained Torsional Shear Tests on Gravelly Soils. Landslides, 1(3): 185-194. https://doi.org/10.1007/s10346-004-0023-3

[11] Leuther, F., Schlüter, S., 2021. Impact of Freeze-Thaw Cycles on Soil Structure and Soil Hydraulic Properties. SOIL, 7(1): 179-191. https://doi.org/10.5194/soil-7-179-2021

[12] Liu, J. J., Zha, F. S., Xu, L., et al., 2020. Strength and Microstructure Characteristics of Cement-Soda Residue Solidified/Stabilized Zinc Contaminated Soil Subjected to Freezing-Thawing Cycles. Cold Regions Science and Technology, 172: 102992. https://doi.org/10.1016/j.coldregions.2020.102992

[13] Liu, Y. H., Chen, Z. M., Guo, L. M., et al., 2023. Shear Characteristics of Frozen Moraine Soil under Freeze-Thaw Cycle. China Sciencepaper, 18(2): 166-171, 203 (in Chinese with English abstract).

[14] Lü, S. Z., Wang, R., Hu, M. J., et al., 2014. Computerized Tomography (CT) Scanning Test Research on Intact Moraine Soil on West Side of Yulong Snow Mountain. Rock and Soil Mechanics, 35(6): 1593-1599, 1622 (in Chinese with English abstract).

[15] Peng, J. B., Zhang, Y. S., Huang, D., et al., 2023. Interaction Disaster Effects of the Tectonic Deformation Sphere, Rock Mass Loosening Sphere, Surface Freeze-Thaw Sphere and Engineering Disturbance Sphere on the Tibetan Plateau. Earth Science, 48(8): 3099-3114 (in Chinese with English abstract).

[16] Potyondy, D. O., Cundall, P. A., 2004. A Bonded-Particle Model for Rock. International Journal of Rock Mechanics and Mining Sciences, 41(8): 1329-1364. https://doi.org/10.1016/j.ijrmms.2004.09.011

[17] Qiu, E. X., He, Q. L., Chen, Q. L., et al., 2023. Influence of Freeze-Thaw Cycles on Mechanical Properties of Moraine Soils. Transportation Geotechnics, 42: 101097. https://doi.org/10.1016/j.trgeo.2023.101097

[18] Qiu, E. X., Pan, H. Y., He, Q. L., et al., 2024. Tests on the Mechanical Properties of Moraine Soils under Freeze-Thaw Conditions and the Modified Duncan-Zhang Model. Journal of Engineering Geology, 32(3): 772-784 (in Chinese with English abstract).

[19] Qu, Z. J., Liu, K. M., Xiao, X. J., et al., 1992. Study of Microstructure, Stress-Strain Behavior and Constitutive Model of Till. Chinese Journal of Geotechnical Engineering, 14(6): 19-28 (in Chinese).

[20] Rossi, A. M., Kendrick, K. J., Graham, R. C., 2019. Pedogenic Evolution on the Arid Bishop Creek Moraines, Eastern Sierra Nevada, California. CATENA, 183: 104222. https://doi.org/10.1016/j.catena.2019.104222

[21] Song, Y. J., Sun, Y. W., Li, C. J., et al., 2023. Meso-Fracture Evolution Characteristics of Freeze-Thawed Sandstone Based on Discrete Element Method Simulation. Rock and Soil Mechanics, 44(12): 3602-3616 (in Chinese with English abstract).

[22] Sun, Y., Li, H., Chen, Z. F., et al., 2023. Numerical Simulation of Freeze-Thaw Damage of Root-Soil Composite Based on Discrete Element Method. Science Technology and Engineering, 23(16): 7025-7032 (in Chinese with English abstract).

[23] Tan, L., Wei, C. F., Tian, H. H., et al., 2015. Experimental Study of Unfrozen Water Content of Frozen Soils by Low-Field Nuclear Magnetic Resonance. Rock and Soil Mechanics, 36(6): 1566-1572 (in Chinese with English abstract).

[24] Tang, M. G., Xu, Q., Deng, W. F., et al., 2022. Degradation Law of Mechanical Properties of Typical Rock in Sichuan-Tibet Traffic Corridor under Freeze-Thaw and Unloading Conditions. Earth Science, 47(6): 1917-1931 (in Chinese with English abstract).

[25] Tao, Y., Yang, P., Li, L., et al., 2023. Characterizing Unfrozen Water Content of Saline Silty Clay during Freezing and Thawing Based on Superposition of Freezing Point Reduction. Cold Regions Science and Technology, 213: 103933. https://doi.org/10.1016/j.coldregions.2023.103933

[26] Veettil, B. K., Kamp, U., 2021. Glacial Lakes in the Andes under a Changing Climate: A Review. Journal of Earth Science, 32(6): 1575-1593. https://doi.org/10.1007/s12583-020-1118-z

[27] Wang, J. Q., Wang, Q., Kong, Y. Y., et al., 2020. Analysis of the Pore Structure Characteristics of Freeze-Thawed Saline Soil with Different Salinities Based on Mercury Intrusion Porosimetry. Environmental Earth Sciences, 79(7): 161. https://doi.org/10.1007/s12665-020-08903-w

[28] Wang, J. T., Zhou, W. J., Dong, G. C., et al., 2024. Repeated Glacial Fluctuations during the Last Glacial Maximum in the Southeastern Tibetan Plateau: ¹⁰Be Surface Exposure Dating of Moraines in the Lahaku Valley, Haizishan Plateau, China. Palaeogeography, Palaeoclimatology, Palaeoecology, 636: 111959. https://doi.org/10.1016/j.palaeo.2023.111959

[29] Wang, X. M., Yin, J., Luo, M. H., et al., 2023. Active High-Locality Landslides in Mao County: Early Identification and Deformational Rules. Journal of Earth Science, 34(5): 1596-1615. https://doi.org/10.1007/s12583-021-1505-0

[30] Wu, L. Y., Zhu, Y. H., Bai, H. B., et al., 2023. Study on the Correlation of Macro and Meso Parameters of Parallel Bond Model Sandstone. Journal of Mining Science and Technology, 8(4): 487-501 (in Chinese with English abstract).

[31] Xiao, D. H., Feng, W. J., Zhang, Z., 2014. The Changing Rule of Loess' s Porosity under Freezing-Thawing Cycles. Journal of Glaciology and Geocryology, 36(4): 907-912 (in Chinese with English abstract).

[32] Xie, S. B., Qu, J. J., Lai, Y. M., et al., 2015. Effects of Freeze-Thaw Cycles on Soil Mechanical and Physical Properties in the Qinghai-Tibet Plateau. Journal of Mountain Science, 12(4): 999-1009. https://doi.org/10.1007/s11629-014-3384-7

[33] Xu, A. H., Yan, Y. H., Chang, D., et al., 2025. Study on Macroscopic and Microscopic Mechanical Properties of Frozen Clay. Journal of Glaciology and Geocryology, 47(2): 372-381 (in Chinese with English abstract).

[34] Xu, W. J., Wang, S., 2016. Meso-Mechanics of Soil-Rock Mixture with Real Shape of Rock Blocks Based on 3D Numerical Direct Shear Test. Chinese Journal of Rock Mechanics and Engineering, 35(10): 2152-2160 (in Chinese with English abstract).

[35] Zhang, Y. S., Guo, C. B., Shi, J. S., et al., 2007. Research on the Engineering Geological Properties of Moraine/Outwash Rocks on the West Side of the Yulong Mountains. Geoscience, 21(1): 150-156 (in Chinese with English abstract).

[36] Zhou, G. G. D., Chen, L. L., Mu, Q. Y., et al., 2019. Effects of Water Content on the Shear Behavior and Critical State of Glacial Till in Tianmo Gully of Tibet, China. Journal of Mountain Science, 16(8): 1743-1759. https://doi.org/10.1007/s11629-019-5440-9

[37] Zhou, X. Q., Xu, W. Y., Niu, X. Q., et al., 2007. A Review of Distinct Element Method Researching Progress and Application. Rock and Soil Mechanics, 28(S1): 408-416 (in Chinese with English abstract).

[38] Zhu, T. T., Chen, J. X., Huang, D., et al., 2021. A DEM-Based Approach for Modeling the Damage of Rock under Freeze-Thaw Cycles. Rock Mechanics and Rock Engineering, 54(6): 2843-2858. https://doi.org/10.1007/s00603-021-02465-4

基金资助

四川省科技计划项目(2025YHNH0008)

四川省科技计划项目(YFNH0004)

地质灾害与环境保护国家重点实验室自由探索课题项目(SKLGP20233006)

成都理工大学2022年中青年骨干教师发展资助计划(10912-SJGG2022-07282)

珠峰大科学计划2.0：面向人地协调的青藏高原冰湖溃决链式灾变机理与防控研究