裂隙饱水砂岩冻融损伤演化与受载破裂特征细观模拟研究

Meso-simulation Study on Freeze‒Thaw Damage Evolution and Loading Fracture Characteristics of Fractured Water‒Saturated Sandstone

摘要

寒区岩土工程通常涉及含有节理和裂隙的岩体，其内部结构面饱水后在冻融循环的作用下易引发一系列工程问题。本文基于离散元水冰颗粒相变膨胀方法，开展裂隙饱水砂岩冻融循环及单轴压缩数值模拟，并结合室内力学试验及CT（computed tomography）进行对比验证，从细观上探讨了裂隙饱水砂岩的冻融损伤演化规律和细观受载破裂特征。结果表明：裂隙饱水砂岩冻融损伤程度随倾角不同而异，损伤程度从大到小对应的倾角依次为0°、15°、45°、60°、90°、75°，伴随冻融作用进入后期阶段，次生裂纹在既有裂隙周边萌生与扩展；冻融循环过程中张拉微裂纹占据主导地位，裂隙周围因水冰相变的冻胀扩容作用出现应力集中现象，在冻融后期，低倾角裂隙左侧的次生裂纹与边界贯通应力得以释放；随着冻融循环次数的增加，试样的峰值应力逐渐降低，试样受载破裂过程中微裂纹数量呈“慢—缓—陡”的趋势演化，张拉微裂纹数目远多于剪切微裂纹，试样张拉破坏模式显著；冻融前后试样破裂过程与形态明显不同，试样受载诱发冻融产生的次生裂纹互相贯通，进而发育形成翼裂纹，随着冻融循环次数的增加，裂纹扩展更加剧烈，最终破碎程度更为完全。研究可为探索裂隙岩体冻融损伤机制及破坏预测方法提供参考。

Abstract

Objective Geotechnical engineering in cold regions typically involves rock masses containing joints and fractures. During seasonal freeze-thaw cycles, water trapped within these internal structural planes undergoes freezing and thawing, resulting in freeze-thaw damage to the rock mass. This damage leads to the deterioration of mechanical properties and overall structural stability, significantly increasing the risk of engineering accidents. As engineering activities in cold regions continue to expand, the problem of freeze-thaw degradation in fractured rock masses becomes increasingly significant. Therefore, this study investigates the degradation effects and damage mechanisms of saturated fractured rock masses during cyclic freeze-thaw processes at the micromechanical level. The objective is to accurately assess how primary fractures influence the development of secondary fractures during freeze-thaw cycles and to determine whether micromechanical damage alters the loading-fracture characteristics of rock masses. This research provides new insights into the mechanisms of freeze-thaw damage and failure prediction methods for rock slopes in cold regions and advances the interdisciplinary field of rock mechanics and numerical simulation. Methods Fractures within rock masses were connected through rock bridge zones. When external loads caused rock bridge failure, step-like sliding occurred on slopes with discontinuous joints of varying steepness. Therefore, in laboratory experiments, three carefully selected rock specimens with fractures (high dip angle of 75°, low dip angle of 15°, rock bridge length of 15 mm, dip angle of 90°) were subjected to freeze-thaw cycles and uniaxial compression tests after saturation. CT scans were performed on the specimens after they reached predetermined freeze-thaw cycles (0, 15, 30 cycles) or after loading. Using Avizo software, three-dimensional reconstructed images enabled analysis of micro-parameter evolution and quantitative fractal geometric characterization. In numerical simulations, a pore-water-saturated rock model was established in Particle Flow Code (PFC) based on the discrete element method to represent the phase-change expansion of water and ice particles. Parallel-bonding contact models were selected to represent particle cohesion, and the micro-parameters were calibrated using laboratory mechanical test results and trial-and-error methods. The entire freeze-thaw cycle was controlled using Fish language to simulate the coupled expansion of pore water and ice particles during phase transitions. Microcrack evolution diagrams, crack rose diagrams, displacement plots, and force chain diagrams were recorded and processed for different freeze-thaw cycle counts (0, 10, 15, 20, 25, 30 cycles). The extent to which internal structural degradation influenced local strain development during loading deformation and determined failure characteristics remained unclear. Therefore, subsequent uniaxial compression simulations were conducted on freeze-thaw-cracked saturated sandstone. These simulations employed deformation-controlled loading at a rate of 0.05 mm/min until specimen failure occurred. Microcrack evolution diagrams, displacement plots, and force chain diagrams were recorded and processed for different loading stages. Results and Discussions In indoor experiments, the processed CT images underwent three-dimensional reconstruction to evaluate the effects of freeze-thaw damage on distinct fracture zones independently. The reconstructed 3D model was segmented into Fracture Zone I, the rock bridge region, and Fracture Zone Ⅱ. As the number of freeze-thaw cycles increased, both the increment in Euler number and the increment in fracture volume of the specimens gradually increased, indicating an accelerated progression of freeze-thaw damage during the later stages. Second, the study calculated the surface porosity of each layer based on the CT scan images. The increment in surface porosity during the later freeze-thaw stage exceeded that observed in the earlier stage, which further validated the accelerated progression of freeze-thaw damage in the later phase. In addition, the study revealed that damage was concentrated in Fracture Zone I during both the early and late freeze-thaw stages, while the rock bridge region exhibited minimal damage. In numerical simulations, microcrack development during the early freeze-thaw phase primarily occurred on specimen surfaces and around joints. As the number of freeze-thaw cycles increased, microcracks progressively expanded inward from the specimen periphery and joint margins. Then, crack clustering emerged around the joints and exhibited localized cumulative effects. The intensity of crack clustering was significantly higher around low-angle joints than around high-angle joints. After 20 freeze-thaw cycles, secondary cracks induced by freeze-thaw damage became observable around the low-angle joints, which caused severe specimen degradation. This finding further confirmed that low-angle joints sustained greater damage than high-angle joints during freeze-thaw cycles. At any freeze-thaw cycle count, tensile microcracks significantly outnumbered shear cracks and dominated the damage pattern. Further analysis of the displacement and force chain diagrams revealed that particle displacement and freeze-induced swelling forces were primarily concentrated around the fractures. After 25 freeze-thaw cycles, a distinct low-force chain line appeared from the left tip of the low-angle fractures to the specimen boundary. This occurred because secondary fractures generated by freeze-thaw damage propagated to the boundary and caused stress release. Regarding the evolution of microcracks and fracture mechanisms during the simulated uniaxial compression of freeze-thaw fractured sandstone, the study identified four distinct loading stages: (Ⅰ) crack-free stage, (Ⅱ) slow crack development stage, (Ⅲ) rapid crack propagation stage, and (Ⅳ) through-crack failure stage. A comprehensive evaluation of fracture evolution under loading for both unfrozen specimens and specimens subjected to 30 freeze-thaw cycles revealed several similarities. Fractures initiated on the right side of the low-angle fractures, propagated along the rock bridge regions, and connected with the high-angle fractures. Both conditions exhibited block detachment on the left side of the low-angle fractures. However, clear differences also emerged. Multiple inner wing cracks were observed around the low-angle fractures in the specimens subjected to 30 freeze-thaw cycles during failure. Crack propagation was more intense, and failure was more complete in these specimens, which indicated a greater destructive potential in freeze-thaw-cycled specimens. Conclusions This study employs the discrete element method to simulate freeze-thaw cycles and uniaxial compression in fractured, water-saturated sandstone by modeling the phase change expansion of water-ice particles. The results are validated through laboratory mechanical tests and computed tomography (CT) analysis. At the microscale, the evolution of freeze-thaw damage and the characteristics of microfractures in fractured, water-saturated sandstone are investigated. The results indicate that freeze-thaw damage primarily concentrates near low-angle fractures, which trigger secondary crack propagation during the later stages of the cycles. Tensile microcracks dominate throughout the freeze-thaw process, while stress concentration occurs around fractures due to freeze expansion caused by water-ice phase transitions. During loading-induced failure, microcracks exhibit a "slow-moderate-rapid" evolution trend. Tensile microcracks significantly outnumber shear microcracks, indicating a pronounced tensile failure mode. Secondary cracks induced by freeze-thaw cycles become interconnected under loading. As the number of freeze-thaw cycles increases, crack propagation intensifies, ultimately leading to more complete fragmentation. This study provides a theoretical foundation for improving operational safety in geotechnical engineering in cold regions.

Graphical abstract

关键词

冻融循环 / 裂隙岩体 / 离散元 / 细观损伤 / CT扫描 / 破裂特征

Key words

freeze‒thaw cycle / fractured rock mass / discrete element / mesoscopic damage / CT scan / rupture characteristics

本刊网刊

09 ‒ 18

网络出版日期:2025 ‒ 09 ‒ 28

近年来，随着“一带一路”倡议的深入推进，中国西部寒区基础设施建设快速发展，隧道、矿山、水电等大型岩体工程相继展开^[1]。昼夜交替和季节性变化引起的自然冻融循环，使得寒区岩石物理力学性能发生劣化，极易引发如岩崩^[2]、衬砌开裂^[3‒4]、边坡滑移失稳^[5]等工程地质问题，这不仅危及工程安全，同时给后期运营埋下重大隐患。因此，深入研究冻融环境下岩石的损伤劣化规律具有重要的工程实践意义。

岩体作为一种天然地质结构，内部常存有褶皱、节理、断层等结构面。在冻融、渗流、载荷等作用下，结构面间相互影响并通过次生结构面扩展贯通，从而进一步降低岩体的力学特性^[6]。此外，寒区岩质边坡的宏观裂隙通常赋存一定水分，其在冻结过程中较微裂隙可产生更为严重的冻融损伤进而驱动岩体破坏，因此冻害研究的关键在于裂隙冰^[7]。另一方面，冻融损伤评价还需考虑裂隙饱和度的影响^[7‒8]，饱水冻融条件（冻结时裂隙充冰，融解时裂隙充水）更易引发岩体起裂扩展^[9]。因而，在高海拔区域，含有饱水节理的岩质边坡愈发不稳定^[10‒11]，这将对寒区道路的选线设计与安全运营产生显著影响。虽然已有研究通过现场监测^[12]和理论模型^[13]探讨了冻结岩石的破裂特征，但却难以直观地分析岩石材料损伤的内禀机制。

岩石裂隙的无序分布是其宏观非线性变形的主要原因，而宏观变形破坏是细观裂纹萌生、扩展、贯通的综合效应。因此，从细观角度探究岩石的结构演化对其宏观力学特性影响具有重要意义^[14]。目前，CT扫描技术凭借其卓越的三维高精度成像能力，已经成为研究岩石细观结构演化的重要工具^[15]。Maji等^[16]采用CT设备观测到微裂隙向宏观裂隙的转化过程，并对试件破坏阶段的裂纹扩展进行定量评估；De Kock等^[14]使用CT扫描技术观测到岩石在冻结期间的裂隙萌生、冰结晶放热，以及融解期间的裂隙闭合现象，发现当冰晶含量达裂隙体积的4.3%时，结晶压力超过岩石强度引发起裂；Liu等^[17]利用高精度CT扫描仪获得不同冻融次数下砂岩的图像，并基于图像中的细观参数定义损伤变量，建立了该损伤变量与冻融循环之间的关系；宋勇军等^[18]对不同饱和度的砂岩进行冻融试验，并结合CT技术分析了砂岩渗透率、孔喉参数、分形维数等的变化。尽管针对冻融岩石细观损伤机制的研究已取得了丰硕成果，但如何精准观测和定量表征岩石在循环冻融过程中内部细观微裂纹的萌生、扩展及局部演化机制，仍是当前研究亟待突破的瓶颈。

数值技术是研究岩石复杂力学行为的一个重要手段，并为解决岩石破裂等问题提供了有效途径^[19‒20]。现有数值方法主要为有限元法和离散元法。有限元法受限于连续性假设，不能直接捕捉岩石裂纹的萌生、扩展等过程^[21]，故而部分学者将研究重点转向能够重现岩石内部微裂纹演化全过程的离散元方法^[22‒23]。颗粒流法作为一种基于离散元的方法，已广泛应用于岩石力学研究领域。Lin等^[24]建立了循环冻融节理岩石的冻胀压力计算模型，采用颗粒流法对岩石节理内部和尖端冻胀应力的分布进行研究；宋勇军等^[25]利用PFC^3D模拟冻融岩石力学行为，分析了岩石的冻融受载变形过程和破裂演化特征；刘佳诺等^[26]利用颗粒流软件PFC^2D对冰碛土试样进行冻融循环及三轴压缩数值模拟，揭示了冰碛土在冻融作用下的微观结构变化、力学性能劣化规律及损伤累积过程。综上所述，虽然针对冻融岩石细观损伤特性的研究已有大量成果，但鲜有研究考虑原生裂隙的影响。裂隙水相变产生的冻胀力和冻融变形会极大程度影响岩石的完整性和力学特性，而冻融中原生裂隙饱水如何影响次生裂隙发展，细观损伤是否改变岩石受载破坏特征等问题尚缺乏相关研究。

基于以上考虑，本文借助水冰颗粒相变耦合的方法，开展裂隙饱水岩石冻融循环及单轴压缩模拟，并通过室内力学试验结果校准模型细观参数，同时结合CT扫描进行对比验证，旨在研究裂隙饱水岩石冻融过程中微裂纹演化机制，并对其受载时内部微裂纹衍生及破裂演化特征进行深入探究。研究方法和试验结果可为寒区道路边坡等岩土工程的施工与安全运营提供参考依据。

1 室内试验

1.1 试样制备

为校准数值模型细观参数并评价数值模拟可靠性，开展室内试验。试验所用白砂岩采自中国西部某边坡工程，通过X射线衍射（XRD）分析可知该砂岩主要由石英（质量占比54.1%）、高岭石（质量占比37.2%）、钾长石（质量占比5.2%）和其他矿物（质量占比3.5%）组成，XRD图谱如图1所示。为确保试样均质性，选取质地均匀、无明显层理的新鲜岩块，沿同一方向钻取岩芯。而后，基于ISRM建议的方法，将岩块加工成直径φ50 mm×高H100 mm的标准岩样。

天然岩体中的节理通常具有不连续性，而工程选址通常会避开大型结构面，因此短小断续节理成为工程实践中的常见问题^[27]。岩体中的裂隙通过岩桥区域连接，当外部荷载导致岩桥失效时，缓‒陡型断续节理边坡可能发生阶梯状滑移^[28]。这种破坏模式在含有低或高倾角节理组的岩质边坡中尤为常见。为此，选取15°与75°的倾角，模拟实际工程中断续节理岩体的典型工况。试验中：裂隙迹长统一设置为20 mm，宽度为1 mm；岩桥长度为15 mm，倾角为90°。试件通过钻孔和线切割技术预制裂隙，其空间几何特征及试样成品如图2所示。

岩样制备后，将其置于真空干燥箱中在105 ℃下烘干48 h；随后，测定干密度和纵波波速；接着，利用真空饱和装置（真空度-0.1 MPa）进行24 h饱水处理；最后，基于孔隙率及其他物理参数，筛选出3个裂隙试样用于裂隙饱水冻融试验。岩样基本物理参数见表1。

1.2 试验流程

试验分为冻融循环和单轴压缩，分别在达到设定冻融次数或加载至峰值应力后对试样进行CT扫描，试验岩样及设备如图3所示。试验步骤具体为：1）冻融循环。设定-20 ℃为冻结温度，20 ℃为融解温度，单次循环时间为24 h（图4）。循环次数为0、15、30次，在达到预设的冻融次数后进行CT扫描。为实现裂隙饱冰，将试样装入等效直径为50.93 mm的矩形塑料自封袋内冻结。由于塑料袋紧密贴合且上部开口，试件表面难以形成有效冻胀压力。2）单轴压缩。以0.05 mm/min的速率进行加载，达到峰值应力后卸载并进行CT扫描。CT扫描采用天津三英精密仪器有限公司生产的MultiscaleVoxel‒2000工业CT进行，其曝光时间为0.5 s，电流为1.3 mA，电压为385 kV。单次扫描可获得1 200像素×1 200像素的16 bit灰度图2 400张，空间分辨率约为41.67 μm。

1.3 细观损伤特征

对经历不同冻融循环次数的试样CT扫描图像进行后期处理。首先，基于Avizo软件，采用非局部均值滤波降噪。非局部均值滤波作为一种线性滤波，可去除外界环境及岩样本身等因素造成的噪音及环状伪影，以利于试验数据的处理。随后，进行边缘裁剪并使用射束硬化校正技术以减少由射束硬化引起的伪影。最后，使用交互法阈值和交互式顶帽算法对CT扫描图像进行分割，通过实际测量的孔隙率反演确定每个试件的特定阈值^[7]，具体处理流程如图5（a）所示。由于二维CT图像仅能反映岩石某一特定截面的结构特征，无法直观地展示出岩石内部空间结构的整体情况。因此，进一步地，将处理过的CT图像进行三维重建。为对不同倾角裂隙的冻融损伤效果进行单独评估，将重建的三维模型分为裂隙Ⅰ区、岩桥区域、裂隙Ⅱ区，如图5（b）所示。

基于CT扫描数据，可以提取多种细观结构参数。在三维分析中，常采用欧拉数和裂隙体积作为评价指标，用于表征岩石受冻融作用后的损伤程度。欧拉数作为多孔介质连通性的量化指标，能够间接反映裂隙网络中水的迁移能力；裂隙体积则直接表征了冻融循环导致的裂隙扩展程度。欧拉数E_{\mathrm{u}}可由下式^[29]计算：

E_{\mathrm{u}}={\beta }_{\mathrm{0}}-{\beta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}（1）

式中，{\beta }_{\mathrm{0}}为孤立孔隙数量，{\beta }_{\mathrm{1}}为连通孔隙数量，{\beta }_{\mathrm{2}}为封闭孔隙数量。E_{\mathrm{u}}值越高，表明结构连通性越差；喉道充分发育的裂隙网络的E_{\mathrm{u}}值通常为负。

表2为冻融砂岩每百层切片的细观参数变化，{\Delta }_{\mathrm{1}}和{\Delta }_{\mathrm{2}}分别表示第2次和第3次CT扫描获得的参数相较于前一次扫描的变化量。由表2可知，试件不同区域欧拉数与裂隙体积的{\Delta }_{\mathrm{2}}均大于{\Delta }_{\mathrm{1}}，这表明冻融损伤进入后期呈加速发展的趋势。此外，裂隙Ⅰ区的裂隙体积增量明显大于其他区域，说明低倾角裂隙周边受冻融影响而引起的损伤更加剧烈。

面孔隙率n_{\mathrm{s}}为材料内部某一剖切面上孔隙的面积占所在切面的面积之比，可由下式计算：

n_{\mathrm{s}}=\frac{S_{\mathrm{p}}}{S}\times \mathrm{100}\mathrm{\%}（2）

式中，S_{\mathrm{p}}指某剖切面上孔隙的面积，S为该剖切面的总面积。

为了深入研究预制裂隙对试样孔隙结构的影响，基于CT扫描图像计算了每一层的面孔隙率，图6为面孔隙率增量随切片层数增加的变化趋势及次生裂隙形貌特征。由图6（a）可知，冻融后期的面孔隙率增量大于冻融前期面孔隙率增量，这进一步验证了冻融损伤在后期加速发展。此外，可以注意到，不论是冻融前期还是冻融后期，损伤均集中于裂隙Ⅰ区，而岩桥区域损伤最小，这主要是因为裂隙内发生水冰相变，对周围岩石颗粒造成冻胀挤压，导致裂隙周围出现次生裂隙的发展，而岩桥区域的裂隙扩张因裂隙水迁移的滞后性而发展缓慢，故面孔隙率增量较低。由图6（b）可知，在试样预制裂隙周围可以观察到由冻融损伤引起的次生裂隙，这些次生裂隙的扩展与连通导致试样内部结构的完整性遭到破坏，进而导致试样宏观力学性能的显著退化。

为进一步探究不同倾角裂隙饱水冻融受到的损伤存在差异的具体机制，从断裂力学角度进行研究。根据断裂力学，对于一条等宽、扁平、完全冻结的裂隙，初始裂隙在产生宏观裂纹前，可简化为竖直裂隙面处的恒张拉力模型^[30]。故在宏观裂隙产生之前，内部冻结模型可以简化为垂直裂隙面处的法向拉伸模型。Ⅰ型开裂情况对应的应力强度因子K_{\mathrm{I}}的表达式如下：

K_{\mathrm{I}}={\sigma }_{\mathrm{z}}\sqrt[]{\frac{\mathrm{\pi }a}{\mathrm{2}}}（3）

式中，{\sigma }_{\mathrm{z}}为裂隙中的冻结应力，a为裂隙长度。

K_{\mathrm{I}\mathrm{C}}为材料的平面应变断裂韧性，在冻融循环作用下，当K_{\mathrm{I}}>K_{\mathrm{I}\mathrm{C}}时^[31]，裂隙尖端开始开裂。裂缝长度r_{\mathrm{f}}可以由下式计算：

r_{\mathrm{f}}=\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}{(\frac{K_{\mathrm{I}\mathrm{C}}}{{\sigma }_{\mathrm{t}}})}^{\mathrm{2}}（4）

式中，{\sigma }_{\mathrm{t}}为岩石抗拉强度。

裂隙区初始冻融循环过程可视为纯Ⅰ型裂纹断裂过程，其初始扩展方向应与断裂方向大致相同，扩展长度应与岩石的断裂韧性和抗拉强度有关。由于裂隙呈扁平状，可以认为冻融作用主要造成裂隙端部产生拉应力集中，并沿裂隙走向^[32]产生二维裂隙扩展。因此，沿裂隙倾角的切割剖面可认为是裂隙^[33]不同倾角的岩体应力弱面，同时将抗冻裂范围X_{\mathrm{f}\mathrm{t}}视为冻融局部化损伤抗力范围。图7为二维剖面中抵抗裂缝扩展的最大范围。图7中，X_{\mathrm{f}\mathrm{t}\mathrm{1}}为低倾角裂隙的抗冻裂范围，X_{\mathrm{f}\mathrm{t}\mathrm{2}}为高倾角裂隙的抗冻裂范围，O_{\mathrm{1}}为低倾角裂隙的中心，O_{\mathrm{2}}为高倾角裂隙的中心，{\theta }_{\mathrm{1}}为低倾角裂隙的水平夹角角度，{\theta }_{\mathrm{2}}为高倾角裂隙的水平夹角角度，{\theta }_{\mathrm{3}}为低倾角裂隙右侧至试样右上角所连直线与基准水平线所成夹角的角度，L_{\mathrm{1}}为试样的半径，L_{\mathrm{2}}为高倾角裂隙上端至试样下平面的垂直距离。根据Shen等^[34]的研究成果，抗冻裂范围X_{\mathrm{f}\mathrm{t}}可按下式计算：

（5）

式中，\theta为裂隙的水平夹角角度，a为裂隙的长度，L为试样半径。

不同裂隙倾角下的抗冻裂范围如图8所示。由图8可知：当裂缝倾角小于59.5°时，抗冻裂范围X_{\mathrm{f}\mathrm{t}}的大小随裂隙倾角的增大而增大；当裂隙倾角大于59.5°时，抗冻裂范围X_{\mathrm{f}\mathrm{t}}随裂隙倾角的增大而减小。这就解释了为什么15°裂隙倾角时受到的冻融损伤大于75°裂隙倾角时受到的损伤。由图8还可知，裂隙饱水冻融损伤程度由大到小对应的裂隙倾角依次为0°、15°、45°、60°、90°、75°，Tan等^[35]也得到了相同的结论。

图9为不同冻融次数后砂岩试样的CT二维纵向切片图像。由图9可知：冻融循环15次后，低倾角裂隙右侧尖端开始出现次生裂隙；当冻融次数达到30次时，低倾角裂隙两侧尖端、高倾角裂隙上方尖端及试件右侧区域均出现了较明显的次生裂隙，而高倾角裂隙下方仍未见明显次生裂隙。根据前文断裂力学分析，高倾角裂隙因其抗冻裂范围较大、抗冻融损伤能力较强，该区域内裂隙发育不完全。

为更准确评估试样在不同冻融次数后内部细观损伤的演化过程，图10给出了CT二维横向切片图像。

由图10可知：冻融20次后，试样上、中、下部的边缘区域均已出现裂纹；随着冻融次数增加，裂纹持续扩展并相互贯通，导致岩石颗粒大量脱落，甚至在试样外围出现片状或块状剥落；冻融30次后，试样劣化程度显著加剧，且损伤主要集中于外围区域。持续的冻融作用促使劣化自表及里、由外向内进一步发展。

1.4 力学特征

图11为不同冻融循环次数下砂岩单轴压缩应力‒应变曲线。由于试样达到峰值应力后即刻卸载进行CT扫描，故曲线无峰后阶段。由图11可知，随着冻融次数的增加，试样峰值应力与弹性模量逐渐降低。材料的宏观力学行为从根本上取决于其细观裂隙网络的演化特征；结合CT扫描可知，冻融循环导致试样内部裂隙与孔隙增多，从而严重弱化了岩石抵抗变形的能力。

综上所述，实验室测试手段多聚焦于岩石宏观力学响应的获取，受观测尺度与精度的限制，细观层面上微结构演化、孔隙分布特征及裂纹扩展路径等关键信息仍难以被系统捕捉和量化分析。CT扫描技术固然可对岩石冻融循环特定阶段的细观损伤进行一定的研究，但由于室内试验的局限性，以及仅靠无损检测手段难以直观地分析岩石材料损伤的内禀机制，更多时候只能间接得出相关结论。因此，本文借助数值方法辅助室内试验，深入分析裂隙饱水岩石在冻融循环条件内部微裂纹演化过程，并对其受载破裂特征进行细观模拟研究。

2 冻融过程数值模拟

2.1 模拟原理

冻融循环导致岩石性能退化的核心机制在于其微观层面的结构变化。当水发生相变结冰时，其体积膨胀会对矿物基质产生膨胀压力，这种压力削弱了颗粒间的结合力，同时促使孔隙结构重组和微裂隙形成。在解冻过程中，融水又会填充进新形成的孔隙空间。随着冻融过程的反复进行，这种结构性损伤会持续累积，进而引发岩石力学强度下降和渗透性变化等一系列退化现象，最终严重影响岩体的工程稳定性。

本文选择颗粒流软件（PFC）模拟这一过程。PFC将粒子系统中的每个粒子视为具有质量、形状和机械特性的独立离散元素。粒子之间由一系列节点和连接点连接，形成网络结构。在细观层面上，PFC可以模拟水/冰颗粒和岩土颗粒，有利于建立考虑冻结的裂隙砂岩模型，同时也是观察损伤裂纹演化的有效途径。此外，与PFC^3D相比，选择PFC^2D模型可以更直观地观察岩石在冻融和加载过程中的损伤情况。

平行黏结接触模型由Potyondy和Cundall^[36]提出，该模型基于线性模型，可在颗粒与颗粒之间或颗粒与墙之间黏结，同时传递力与力矩。黏结破坏后，该模型退化为线性接触模型。基于上述原理，平行黏结模型可以模拟颗粒系统中颗粒膨胀、碰撞、断裂等现象，更适合反映冰颗粒冻胀对岩石骨架造成的破坏。因此，本文选择平行黏结模型用于样品键合，其原理如附录A图A1所示^[37]。

通过体应变等效的思想，利用水颗粒体积膨胀产生的力来表征岩石中的有效冻胀应力，假设如下：

1）将饱和砂岩试样简化为由岩石矿物颗粒和孔隙水颗粒组成。

2）冻融过程中不考虑岩石内部孔隙水分迁移。

3）随温度变化（升温或降温），孔隙水颗粒达到目标温度的时间是一致的。

4）岩石矿物颗粒和冰粒的变形均假定为弹性变形。

在一次冻融循环中，单个孔隙水颗粒初始面积S_{\mathrm{0}}可由下式表示：

S_{\mathrm{0}}=\mathrm{\pi }{r}_{\mathrm{0}}^{\mathrm{2}}（6）

式中，r_{\mathrm{0}}为孔隙水颗粒初始半径。单个孔隙水颗粒相变冻结后面积S_F可由下式计算：

S_{\mathrm{F}}=\mathrm{\pi }{\left(r_{\mathrm{0}}+u_{\mathrm{v}}\right)}^{\mathrm{2}}=S_{\mathrm{0}}+\mathrm{\Delta }{S}_{{}_{\mathrm{W}}}（7）

式中，\mathrm{\Delta }{S}_{{}_{\mathrm{W}}}为孔隙水颗粒冻结后面积增量，u_{\mathrm{v}}为孔隙水颗粒半径膨胀量。

根据刘泉声等^[38]对低温下岩石孔隙冻胀模型的研究成果，u_{\mathrm{v}}由下式计算：

u_{\mathrm{v}}=r_{\mathrm{0}}\cdot \frac{p_{\mathrm{i}}}{E_{\mathrm{m}}}\cdot \frac{[\mathrm{1}+v_{\mathrm{m}}+\mathrm{2}(\mathrm{1}-\mathrm{2}{v}_{\mathrm{m}})n]}{\mathrm{2}(\mathrm{1}-n)}（8）

式中，E_{\mathrm{m}}、{\nu }_{\mathrm{m}}分别为岩石基质的弹性模量与泊松比，n为岩石孔隙率，p_{\mathrm{i}}为岩石孔隙冰压力。

根据Dana等^[39]的研究成果，岩石孔隙冰压力大小可以表示为：

p_{\mathrm{i}}=\frac{\mathrm{0.029}}{\frac{\mathrm{1}}{E_{\mathrm{m}}}\cdot \frac{\mathrm{1}+\mathrm{2}n+\left(\mathrm{1}-\mathrm{4}n\right){\nu }_{\mathrm{m}}}{\mathrm{2}\left(\mathrm{1}-n\right)}+\mathrm{1.029}\cdot \frac{\mathrm{1}-\mathrm{2}{\nu }_{\mathrm{i}}}{E_{\mathrm{i}}}}（9）

式中，E_{\mathrm{i}}、{\nu }_{\mathrm{i}}分别为孔隙冰的弹性模量与泊松比，分别为9 GPa、0.35。

由于岩体基质中的表面吸附及孔隙毛细作用，始终保持一定含量的未冻水，而未冻水不影响岩石的冻结膨胀，因此只考虑冻结孔隙中的冰压力。当孔隙水完成相变时，孔隙冰的体积占据全部孔隙体积（不考虑非自由未冻水的影响），孔隙冰压力作用于孔壁，导致岩石孔隙空间膨胀。在任何温度下，孔隙水/冰粒的面积可以由下式计算：

S=\left\{\begin{array}{l}\begin{array}{l}{S}_{\mathrm{0}}+\mathrm{\Delta }{S}_{\mathrm{w}}\left[\mathrm{1}+\mathrm{0.139}{\left(\frac{\mathrm{1}}{\mathrm{\Delta }T}\right)}^{\mathrm{1}/\mathrm{3}}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{1}+{\mathrm{e}}^{-\mathrm{0.268}\mathrm{\Delta }T}}{\mathrm{2}}\right)\right]·\\ \left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{0.268}\mathrm{\Delta }T}\right),T⩽\mathrm{0}°\mathrm{C};\end{array}\\ S_{\mathrm{0}},T>\mathrm{0}°\mathrm{C}\end{array}\right.（10）

式中：S为孔隙水/冰颗粒在任意温度下的面积；\mathrm{\Delta }T=T_{\mathrm{m}}-T，T_{\mathrm{m}}为常压下自由体积水冻结点（标准大气压下其值为273.15 K），T为当前获取的温度。

在冻融循环模拟中，由于采用颗粒膨胀法表征冻融循环可以忽略实验室试验中恒温段的影响，结合自然界中气温-时间的变化规律，将试样温度随时间的变化简化为温度变化的余弦函数^[40]，以控制温度的快速和缓慢变化。

整个试样冻融循环通过编写Fish语言控制，孔隙水颗粒随温度变化发生水冰颗粒相变耦合膨胀，计算路线如图12所示。附录A图A2给出了PFC模拟冻融循环过程的原理图，描述了在冻融循环期间建模的水颗粒的冻结和膨胀过程，其中冻结和膨胀力来自由相变引起的水颗粒的体积膨胀。

2.2 模型构建及参数标定

图13为裂隙岩石数值模型。模型高度为100 mm，宽度为50 mm，裂隙尺寸与砂岩试样保持一致（图2）。模型由岩石颗粒与水颗粒两种介质组成，裂隙内充满水颗粒以模拟饱水冻融，颗粒间黏结选择平行黏结模型。为使模型尽量精确且更加符合室内效果，通过调整颗粒半径与孔隙率进行多次试验，最终确定颗粒细观参数，见表3。

岩石在冻融循环过程中的劣化是由于岩石骨架的损坏造成的，因此，应将样品中水‒水颗粒间和岩‒水颗粒间的黏结强度设置得足够高，以防止样品中的水颗粒影响测试结果。选取未经历冻融循环的岩样的单轴压缩参数，采用试错法对数值模型中的岩石颗粒力学参数进行标定。模型相关参数见表4。

模型校正结果如图14所示，蓝色虚线的斜率为室内试验岩样的弹性模量。由于数值模型难以完全模拟天然岩样的初始缺陷，故应变结果存在一定偏差。然而，忽略压密阶段^[23]，试样的峰值强度和弹性模量与试验结果高度吻合。此外，数值结果较好地反映了曲线的峰后阶段，故使用表4中细观参数进行数值模拟合理可行。

3 结果与分析

3.1 不同冻融次数下的细观损伤特性

裂隙试样在不同冻融循环次数下的微裂纹演化过程如图15所示。图15中，紫色为试样产生的拉伸微裂纹，黄色为剪切微裂纹。由图15可知：在冻融前期，微裂纹的发育主要集中于试样表面和裂隙周围；随着冻融循环次数的增加，微裂纹逐步从试样周围及裂隙周边向内部扩展，这是因为试样四周的水颗粒受岩石颗粒的约束相较中心位置更少，更容易在冻结膨胀过程中造成岩石颗粒的断裂损伤；而裂隙内水颗粒在冻结状态下体积膨胀，会对周围的岩石颗粒产生冻胀力，导致岩石颗粒间的接触破坏，试样四周及裂隙周边产生大范围微裂隙。

进一步地，裂隙周边出现裂纹聚集，局部表现出累积效应，低倾角裂隙周边的裂纹聚集程度明显高于高倾角裂隙周边。当冻融循环次数达到20次以后，低倾角裂隙周边可以观察到由冻融损伤引起的次生裂纹，造成试样劣化严重。由此可以说明，在冻融循环过程中，低倾角裂隙受到的损伤大于高倾角裂隙。总体来看，微裂隙呈现先慢后快的发展趋势，表明冻融损伤进入后期呈加速发展的趋势。进一步验证了冻融损伤是不断累积的过程，当局部损伤累积到一定程度后会引起岩石受载能力下降甚至是发生破坏。这与第1节所得结论一致，同时验证了模型的可靠性。

为进一步评估冻融作用下微裂纹演变过程，对试样中微裂纹数量和方向进行统计，结果如图16所示，径向方向表示裂纹数量，圆周方向表示裂纹方向。

由图16可知，在任意冻融循环次数下，张拉微裂纹数量都远大于剪切裂纹数量。这与Shen等^[34]通过试验发现岩石冻融易发生拉伸破坏的结论一致。这是由于孔隙中的水颗粒在冻结过程中体积膨胀，产生了冻胀应力，进而对周边的岩石骨架施加压力。当这种冻结应力超越了岩石颗粒黏结所能承受的抗拉强度时，岩石骨架会发生断裂，从而形成张拉微裂纹。进一步观察可知，在冻融前期，张拉微裂纹在竖直方向发育比较明显，而在冻融后期张拉微裂纹沿四周均匀发展。这是因为在冻融初期，颗粒间以水平挤压为主，黏结破坏后产生了垂直于黏结方向的竖向裂隙。而剪切微裂纹在整个冻融过程中无明显发育倾向，表现为无序性。然而，冻融过程中裂纹的演变情况并不足以全面揭示冻融损伤的演化机制。

进一步地，图17给出了不同冻融次数下试样的位移场与力链场的演化过程。由图17（a）可知：冻融前期颗粒位移由端部向中部扩展，这主要是由于靠近试样中心的岩石颗粒受周围颗粒约束较多，故发生的位移相对更小，而试样表面附近的岩石颗粒受到的约束较少；随着冻融循环次数的增加，裂隙周边出现较大的位移，主要集中于低倾角裂隙尖端，向试样边界发展，最大位移达到1.532×10^-4 m。这归因于在冻融初期，裂隙周边岩石颗粒被周围颗粒紧密约束，所以位移不明显，而随着冻融次数的增加，裂隙周边次生裂纹扩展，较小的岩石颗粒在冻胀力的作用下出现较大的位移。

力链是受载时力在颗粒间的传递路径，力链颜色深浅代表力的大小。由图17（b）可知，冻结产生的冻胀力主要集中于裂隙周边，随着冻融次数的增加，冻胀力逐渐增大。值得注意的是，当冻融循环达到25次以后，低倾角裂隙左侧尖端到试样边界存在明显的低力链线，这是因为冻融损伤产生的次生裂隙扩展至边界，引起应力释放。

3.2 冻融裂隙饱水砂岩受载破裂特征

3.2.1 应力-应变与微裂纹演变

众所周知，岩石因循环冻胀/融缩作用产生疲劳损伤，进而影响物理力学性质。但结构内部的劣化在受载变形过程中如何影响局部应变发展，又在何种程度上决定失效破坏特征尚未明了。因此，本文对模拟冻融裂隙砂岩单轴压缩过程中的微裂纹与应力-应变曲线衍生关系进行分析并对冻融裂隙砂岩破裂演化机制进行评估。模拟采用变形控制加载，以0.05 mm/min速率加载至试样破坏。

图18为不同冻融循环次数下试样单轴压缩模拟过程中微裂纹数量随应力-应变变化的曲线。由图18可知，随着冻融循环次数的增加，试样的峰值应力逐渐降低，冻融循环30次后试样的峰值应力相较于冻融0次下降34.9%。由此可以说明，冻融损伤造成岩体内部劣化，进而影响其宏观力学特性，模拟结果较为合理地揭示了冻融作用下岩石内部损伤累积过程。值得注意的是，整个加载过程张拉微裂纹数目均高于剪切微裂纹，说明加载过程试样的失效破坏以张拉作用为主。

分析图18可以明显得出，整个加载过程大致分为4个阶段：1）无裂纹阶段。此时试样被压缩，颗粒间接触力不断增大，但并未达到张拉起裂应力或剪切起裂应力，因此未形成微裂纹。2）裂纹缓慢发育阶段。随着压力不断增大，颗粒间黏结开始断裂，从而形成张拉和剪切微裂纹。3）裂纹快速扩展阶段。微裂纹进入快速增长期，小裂纹逐渐贯通形成宏观裂隙。4）裂纹贯通破坏阶段。裂纹快速贯通，大型宏观破坏带已经形成，试样呈现破碎状态。

试样在进行单轴压缩时，起裂对应的轴向应力为起裂应力阈值。其中，{\sigma }_{\mathrm{N}}^{\mathrm{s}}为剪切微裂纹起裂应力阈值，{\sigma }_{\mathrm{N}}^{\mathrm{t}}为张拉微裂隙起裂应力阈值，{\sigma }_{\mathrm{f}}为峰值轴向应力。将图15中试样受载后拉伸、剪切起裂应力及峰值应力统计于表5。由表5可知，随着冻融循环次数的增加，以上指标均逐渐降低。这表明岩石在经历多次冻融循环作用后，其骨架间的胶结性能会明显减弱；在受到外力荷载时，即便应力水平相对较低，也足以诱发微裂纹的萌生，进而导致岩石的整体承载能力大幅度降低。

3.2.2 受载破裂演化过程

以冻融0次和冻融30次的两个试样作为研究对象，探究试样单次受载破坏演化的全过程，如图19所示。结合图18分析如下：

1）A点前阶段。该阶段为无裂纹阶段，受冻融循环影响，冻融30次后试样初始损伤大于冻融0次的试样，其张拉起裂应力（4.82 MPa）低于未冻融试样的起裂应力（14.78 MPa）。由位移场可知，试样均处于被压缩的状态。由力链场可知，低倾角裂隙两侧出现应力集中现象。由微裂纹演化可知：当加载到A点时，冻融0次试样在低倾角裂隙尖端出现一定数量的微裂纹，但未观察到明显的裂纹聚集现象；而冻融30次后试样会在低倾角裂隙两侧观察到1^#、2^#内翼裂纹。由此可以推断，裂隙饱冰冻融会使低倾角裂隙周围产生次生裂纹，在加载时，次生裂隙互相贯通迅速扩展为翼裂纹，这也是冻融循环引起岩石宏观力学特性降低的根源。

2）AB阶段。此阶段为裂纹缓慢发展阶段，随着加载应力的增大，低倾角裂隙周围颗粒位移明显，表明低倾角裂隙正被压缩，这与室内试验观察到的现象一致。当加载至B点时，冻融0次与30次试样裂纹均从低倾角裂隙右侧开始，沿着岩桥区域不断扩展，逐渐接近高倾角裂隙；由力链场同样可以得到这一结论，可以明显观察到低倾角裂隙右侧存在应力集中，并向下延伸。值得注意的是，冻融30次试样的1^#、2^#内翼裂纹继续向远离裂隙的方向发育。

3）BC阶段。此阶段为裂纹快速扩展阶段，可以观察到，低倾角裂隙已经几乎被压密。冻融0次与30次试样应力依旧集中于低倾角裂隙尖端，并向周围扩散。观察微裂纹分布图19可以发现，冻融0次与30次试样岩桥区域已经基本贯通。此外，冻融0次试样在低倾角左侧产生2^#外翼裂纹，并向试样边界不断扩展；而冻融30次试样损伤进一步加剧，在低倾角裂隙左侧及高倾角裂隙右侧均出现明显的翼裂纹。

4）CD阶段。此阶段为裂纹贯通破坏阶段，此时试样已经破坏。观察位移场可知，冻融0次与30次试样低倾角裂隙左侧均存在块状脱落，冻融30次试样高倾角右侧也存在颗粒剥落。由微裂纹分布图19可知：冻融0次试样高倾角裂隙左侧出现3^#外翼裂纹，2^#外翼裂纹与边界贯通，试样整体已经破碎；而冻融30次试样在低倾角裂隙左侧及高倾角裂隙右下侧均出现大量翼裂纹聚集，试样内部的翼裂纹数量也大于冻融循环0次试样，这说明经过冻融损伤的试样裂纹扩展更加剧烈，破碎更为完全。

综上所述，冻融0次与30次试样受载破裂过程表现出一定的相似性，两者破碎均从低倾角裂隙右侧开始，沿着岩桥区域扩展并与高倾角裂隙贯通，且均在低倾角裂隙左侧产生块状脱落。但两者又存在明显的区别，破碎时，冻融30次试样低倾角裂隙周边存在数条内翼裂纹，且裂纹扩展更加剧烈，破碎更为完全。这归因于裂隙饱水冻融引起低倾角裂隙附近次生裂纹的发育，在低应力状态下，这些次生裂隙会迅速贯通，形成肉眼可见的翼裂纹，对岩样内部结构造成破坏，这也是冻融损伤引起宏观力学性能降低的根源。

为了进一步评估不同饱水冻融次数对试样破坏形态的影响，室内试验中，当试样达到峰值应力后对其卸载并进行CT扫描，图20为裂隙饱水冻融0次与30次试样最终破坏形态。

由图20可知，冻融0次试样最终破坏模式由1^#外翼裂纹、2^#内翼裂纹、3^#反翼裂纹构成，冻融30次试样破坏模式由1^#外翼裂纹、2^#、3^#、4^#内翼裂纹及5^#外翼裂纹构成。值得注意的是，数值模拟试样虽未出现3^#内翼裂纹，但在图19中可观察到低倾角裂隙下侧裂纹扩展现象，而最终破坏形态未能观察到可能是3^#内翼裂纹发育不完全所致。综合来看，冻融0次与30次试样均为拉剪混合破坏，但冻融30次试样破碎更为完全，说明冻融试样破坏危害更大。数值结果与实验室试样破裂形态吻合较好，进一步验证了上述数值模拟试样受载破裂过程的合理性。

4 结 论

本文结合室内试验，通过颗粒流程序对裂隙饱水砂岩冻融细观损伤与破裂过程进行研究。基于水冰颗粒相变耦合膨胀方法，较好地实现了岩石裂隙饱水冻融循环过程模拟，揭示了冻融作用下裂隙饱水岩石细观损伤累积发展过程，深入分析了冻融后试样受载破裂演化行为。主要结论如下：

1）冻融损伤是一个持续累积的过程，微裂纹的演化呈现初期缓慢、后期加速的趋势。裂纹聚集区域主要位于试件边缘及低倾角裂隙周围。冻融循环15次后，在低倾角裂隙周边可明显观察到由冻融作用产生的次生裂纹。裂隙在饱水冻融条件下的损伤程度随倾角不同而异，损伤程度从大到小对应的倾角依次为0°、15°、45°、60°、90°、75°。

2）冻融循环过程中张拉微裂纹占据主导地位。在冻融前期，张拉微裂纹在竖直方向发育比较明显，而在冻融后期张拉微裂纹沿四周均匀发育。冻融过程中裂隙内出现应力集中，冻融25次以后，低倾角裂隙左侧出现低力链线，表明冻融损伤引起的次生裂纹扩展至边界，从而导致应力释放。

3）随着冻融循环次数的增加，试样的峰值应力逐渐降低，冻融循环30次后试样的峰值应力相较于冻融0次下降34.9%。加载过程中微裂纹演化特征呈“慢—缓—陡”的趋势演化，试样产生的张拉微裂纹数目高于剪切微裂纹，表明试样失效破坏以张拉作用为主。

4）试样受载后裂纹均由低倾角裂隙右侧开始发育，沿岩桥区域扩展并与高倾角裂隙贯通，后于裂隙尖端向试样边缘发育，直至破坏。但冻融30次试样在加载初期会引起冻融产生的次生裂纹互相贯通，进而发育形成翼裂纹；且冻融30次试样裂纹扩展更加剧烈，最终破碎程度更为完全。

附录见本刊网络版，扫描标题旁的二维码可阅读网络全文。

参考文献

[1] Tan Hao, Song Yongjun, Yang Huimin,et al.Energy evolution and strain localization in fractured sandstone under freeze-thaw cycling and uniaxial loading-unloading[J].International Journal of Rock Mechanics and Mining Sciences,2024,177:105746. doi:10.1016/j.ijrmms.2024.105746

[2] Yu Yang, Jiang Qincai, Wang Zehua,et al.Dynamic characteristics and damage model of freezing thawing red sandstone under impact load[J].Journal of China Coal Society,2025,50(5):2437‒2452.

[3] 于洋,江亲财,王泽华,等.冲击荷载作用下冻融红砂岩动力特性及损伤模型[J].煤炭学报,2025,50(5):2437‒2452.

[4] Lv Zhitao, Xia Caichu, Wang Yuesong,et al.Analytical el-asto‒plastic solution of frost heaving force in cold region tunnels considering transversely isotropic frost heave of surrounding rock[J].Cold Regions Science and Technology,2019,163:87‒97. doi:10.1016/j.coldregions.2019.04.008

[5] Yan Jian, He Chuan, Yan Qixiang,et al.In-situ test and calculational analysis on frost heaving force of moraine stratum in Que'er moutain tunnel[J].Rock and Soil Mechanics,2019,40(9):3593‒3602. doi:10.16285/j.rsm.2018.1169

[6] 严健,何川,晏启祥,等.雀儿山隧道冰碛地层冻胀力原位测试及计算分析[J].岩土力学,2019,40(9):3593‒3602. doi:10.16285/j.rsm.2018.1169

[7] Jiang Haiqiang, Liu Xingchao, Chen Yuwei,et al.Model test investigation on heat and deformation behaviors of canal slopes with protective layers caused by freeze‒thaw action[J].Cold Regions Science and Technology,2024,228:104338. doi:10.1016/j.coldregions.2024.104338

[8] Draebing D, Mayer T.Topographic and geologic controls on frost cracking in Alpine rockwalls[J].Journal of Geophysical Research:Earth Surface,2021,126(6):e2021JF006163. doi:10.1029/2021jf006163

[9] Tan Hao, Song Yongjun, Guo Xixi,et al.Research on meso-damage and strain localization of fractured sandstone after freeze‒thaw cycles[J].Chinese Journal of Rock Mechanics and Engineering,2022,41(12):2485‒2496.

[10] 谭皓,宋勇军,郭玺玺,等.冻融裂隙砂岩细观损伤与应变局部化研究[J].岩石力学与工程学报,2022,41(12):2485‒2496.

[11] Zhao Tao, Feng Yunfeng, Jia Hailiang,et al.Investigation of the mechanical behaviour of frozen fissured sandstone addressing the role of fissure ice[J].Geomechanics and Geophysics for Geo-Energy and Geo-Resources,2024,10(1):140. doi:10.1007/s40948-024-00860-z

[12] Tan Hao, Song Yongjun, Guo Xixi.Analysis of porosity,permeability,and anisotropy of sandstone in freeze‒thaw environments using computed tomography and fractal theory[J].Fractals,2021,29(8):2150239. doi:10.1142/s0218348x2150239x

[13] Krautblatter M, Funk D, Günzel F K.Why permafrost rocks become unstable:A rock‒ice-mechanical model in time and space[J].Earth Surface Processes and Landforms,2013,38(8):876‒887. doi:10.1002/esp.3374

[14] Mamot P, Weber S, Schröder T,et al.A temperature- and stress-controlled failure criterion for ice-filled permafrost rock joints[J].The Cryosphere,2018,12(10):3333‒3353. doi:10.5194/tc-12-3333-2018

[15] Draebing D, Krautblatter M, Hoffmann T.Thermo-cryogenic controls of fracture kinematics in permafrost rockwalls[J].Geophysical Research Letters,2017,44(8):3535‒3544. doi:10.1002/2016GL072050

[16] Huang Shibing, Liu Quansheng, Liu Yanzhang,et al.Frost heaving pressure and characteristics of frost cracking in elliptical cavity(crack) of rock mass under coupled thermal-mechanical condition at low temperature[J].Chinese Journal of Geotechnical Engineering,2018,40(3):459‒467. doi:10.11779/CJGE201803009

[17] 黄诗冰,刘泉声,刘艳章,等.低温热力耦合下岩体椭圆孔(裂)隙中冻胀力与冻胀开裂特征研究[J].岩土工程学报,2018,40(3):459‒467. doi:10.11779/CJGE201803009

[18] De Kock T, Boone M A, De Schryver T,et al.A pore-scale study of fracture dynamics in rock using X-ray micro-CT under ambient freeze‒thaw cycling[J].Environmental Science & Technology,2015,49(5):2867‒2874. doi:10.1021/es505738d

[19] Song Yongjun, Tan Hao, Yang Huimin,et al.Fracture evolution and failure characteristics of sandstone under freeze-thaw cycling by computed tomography[J].Engineering Geology,2021,294:106370. doi:10.1016/j.enggeo.2021.106370

[20] Maji V, Murton J B.Micro-computed tomography imaging and probabilistic modelling of rock fracture by freeze‒thaw[J].Earth Surface Processes and Landforms,2020,45(3):666‒680. doi:10.1002/esp.4764

[21] Liu Hui, Yang Gengshe, Yun Yehui,et al.Investigation of sandstone mesostructure damage caused by freeze‒thaw cycles via CT image enhancement technology[J].Adva-nces in Civil Engineering,2020,2020:8875814. doi:10.1155/2020/8875814

[22] Song Yongjun, Yang Huimin, Tan Hao,et al.Study on damage evolution characteristics of sandstone with different saturations in freeze-thaw environment[J].Chinese Journal of Rock Mechanics and Engineering,2021,40(8):1513‒1524.

[23] 宋勇军,杨慧敏,谭皓,等.冻融环境下不同饱和度砂岩损伤演化特征研究[J].岩石力学与工程学报,2021,40(8):1513‒1524.

[24] Gu Xinbao, Zhou Xiaoping.Numerical simulation of pro-pagation and coalescence of cracks using peridynamic theory[J].Rock and Soil Mechanics,2017,38(2):610‒616. doi:10.16285/j.rsm.2017.02.038

[25] 谷新保,周小平.裂纹扩展和连接过程的近场动力学数值模拟[J].岩土力学,2017,38(2):610‒616. doi:10.16285/j.rsm.2017.02.038

[26] Liu Quansheng, Wang Zhongwei.Research progress of ro-ck numerical simulation based on digital image processing[J].Chinese Journal of Rock Mechanics and Engineering,2020,39(S2):3286‒3296.

[27] 刘泉声,王中伟.基于数字图像处理的岩石数值模拟研究进展[J].岩石力学与工程学报,2020,39():3286‒3296.

[28] Jing L.A review of techniques,advances and outstanding issues in numerical modelling for rock mechanics and ro-ck engineering[J].International Journal of Rock Mechanics and Mining Sciences,2003,40(3):283‒353. doi:10.1016/s1365-1609(03)00013-3

[29] Ding Zeng, Feng Xiaojun, Wang Enyuan,et al.Acoustic em-ission response and evolution of precracked coal in the meta-instability stage under graded loading[J].Engineering Geology,2023,312:106930. doi:10.1016/j.enggeo.2022.106930

[30] Shi Zhanming, Li Jiangteng, Wang Ju,et al.Experimental and numerical study on fracture characteristics and constitutive model of sandstone under freeze-thaw-fatigue[J].International Journal of Fatigue,2023,166:107236. doi:10.1016/j.ijfatigue.2022.107236

[31] Lin Hang, Lei Daxing, Yong Rui,et al.Analytical and numerical analysis for frost heaving stress distribution wit-hin rock joints under freezing and thawing cycles[J].Environmental Earth Sciences,2020,79(12):305. doi:10.1007/s12665-020-09051-x

[32] Song Yongjun, Sun Yinwei, Li Chenjing,et al.Meso-fra-cture evolution characteristics of freeze-thawed sandstone based on discrete element method simulation[J].Rock and Soil Mechanics,2023,44(12):3602‒3616. doi:10.16285/j.rsm.2023.0448

[33] 宋勇军,孙银伟,李晨婧,等.基于离散元法模拟的冻融砂岩细观破裂演化特征研究[J].岩土力学,2023,44(12):3602‒3616. doi:10.16285/j.rsm.2023.0448

[34] Liu Jianuo, Li Mingli, Jiang Yuanjun,et al.Microscopic damage evolution of moraine soils under freeze-thaw cycles based on PFC^2D simulation[J].Earth Science,2025,50(10):4137‒4154.

[35] 刘佳诺,李明俐,姜元俊,等.基于PFC^2D的冻融循环作用下冰碛土微观损伤研究[J].地球科学,2025,50(10):4137‒4154.

[36] Wu Qiong, Kulatilake P H S W, Tang Huiming.Comparison of rock discontinuity mean trace length and density estimation methods using discontinuity data from an outcrop in Wenchuan area,China[J].Computers and Geotechnics,2011,38(2):258‒268. doi:10.1016/j.compgeo.2010.12.003

[37] Huang Da, Cen Duofeng, Ma Guowei,et al.Step-path failure of rock slopes with intermittent joints[J].Landslides,2015,12(5):911‒926. doi:10.1007/s10346-014-0517-6

[38] Wang Dengke, Zeng Fanchao, Wei Jianping,et al.Quantitative analysis of fracture dynamic evolution in coal subjected to uniaxial and triaxial compression loads based on industrial CT and fractal theory[J].Journal of Petroleum Science and Engineering,2021,196:108051. doi:10.1016/j.petrol.2020.108051

[39] Yu Meilu, Wang Zhongwen, Xu Ying,et al.The failure behavior of prefabricated fractured sandstone with different rock bridge inclination angles under freeze-thaw cycles[J].Unconventional Resources,2025,6:100152. doi:10.1016/j.uncres.2025.100152

[40] Zhou Yun, Tan Xianjun, Yang Diansen,et al.A numerical method to consider the interaction between multiple fractures in frozen rocks based on XFEM[J].Computers and Geotechnics,2024,169:106240. doi:10.1016/j.compgeo.2024.106240

[41] Liang Hehao, Guo Hongjie, Zhang Jianjing,et al.Influence of gradation composition on crack evolution of stone mastic asphalt based on digital image processing[J].Theoretical and Applied Fracture Mechanics,2025,136:104848. doi:10.1016/j.tafmec.2025.104848

[42] Cao Guorui, Longbang Qing, Duan Junbang,et al.Investigations on the crack resistance characteristics of asphalt concrete based on digital image correlation method[J].Engineering Fracture Mechanics,2024,307:110290. doi:10.1016/j.engfracmech.2024.110290

[43] Shen Fengqi, Qiu Wenliang, Xie Kai,et al.A discrete ele-ment-based study of freeze‒thaw damage in water-filled fractured rock[J].Computers and Geotechnics,2025,185:107343. doi:10.1016/j.compgeo.2025.107343

[44] Tan Tao, Zhang Chunyang, Li Wanru,et al.Evolution of freeze‒thaw damage characteristics and corresponding models of intact and fractured rocks under uniaxial compression[J].Rock Mechanics and Rock Engineering,2024,57(10):8013‒8033. doi:10.1007/s00603-024-03996-2

[45] Potyondy D O, Cundall P A.A bonded-particle model for rock[J].International Journal of Rock Mechanics and Mining Sciences,2004,41(8):1329‒1364. doi:10.1016/j.ijrmms.2004.09.011

[46] Itasca Consulting Group Inc.PFC 5.0 help manual[M].Mi-nneapolis:Itasca Consulting Group Inc,2014.

[47] Liu Quansheng, Huang Shibing, Kang Yongshui,et al.Study of unfrozen water content and frost heave model for saturated rock under low temperature[J].Chinese Journal of Rock Mechanics and Engineering,2016,35(10):2000‒2012. doi:10.13722/j.cnki.jrme.2015.1157

[48] 刘泉声,黄诗冰,康永水,等.低温饱和岩石未冻水含量与冻胀变形模型研究[J].岩石力学与工程学报,2016,35(10):2000‒2012. doi:10.13722/j.cnki.jrme.2015.1157

[49] Dana E, Skoczylas F.Gas relative permeability and pore structure of sandstones[J].International Journal of Rock Mechanics and Mining Sciences,1999,36(5):613‒625. doi:10.1016/s0148-9062(99)00037-6

[50] Zheng Zhidan, Shi Zhanming, Wang Mengxiang,et al.Fatigue fracture and energy evolution characteristics of sandstone under water‒ice phase transitions via the discrete element method[J].Computers and Geotechnics,2025,179:107023. doi:10.1016/j.compgeo.2024.107023

基金资助

国家自然科学基金项目(42277182)

国家自然科学基金项目(11972283)