饱和砂岩冻融循环损伤演化特性

Damage Evolution Characteristics of Saturated Sandstone in Freeze-Thaw Cycles

摘要

为揭示寒区岩体工程中岩石在冻融循环作用下变形损伤机制，通过室内冻融循环试验、低场核磁共振试验以及声发射监测试验对饱和砂岩孔隙度变化和宏观强度力学特性进行了分析，建立了饱和岩石冻融循环温度-渗流-应力-损伤耦合模型并进行了验证，开展了饱和砂岩在不同冻融循环次数下孔隙度变化及损伤演化数值模拟.研究结果表明：随着冻融次数增加，砂岩孔隙度增长速率变快，单轴抗压强度降低，且下降速率逐渐加快.孔隙尺寸与数量的变化引起岩石中砂岩强度降低.冻融受载荷作用下砂岩损伤是冻融损伤和载荷损伤共同作用的结果，且随着应变增加，砂岩的损伤变量最终趋于1.研究结果可为寒区岩石力学特性研究提供了理论参考和试验依据.

Abstract

To reveal the deformation and damage mechanism of rock in cold region rock engineering under freeze-thaw cycles， the porosity changes and macroscopic strength mechanical properties of saturated sandstone were analyzed through laboratory freeze-thaw cycle tests， low-field nuclear magnetic resonance experiments， and acoustic emission （AE） monitoring. A coupled temperature， permeability， stress， and damage model for saturated rock under freeze-thaw cycles was developed and validated. Numerical simulations were carried out to investigate the porosity variation and damage evolution in saturated sandstone subjected to different freeze-thaw cycles. The results indicate that as freeze-thaw cycles increase， the porosity growth rate of sandstone increases， while the uniaxial compressive strength decreases， and the rate of strength reduction gradually accelerates. The changes in pore size and quantity lead to deterioration in the strength of sandstone. Freeze-thaw-induced damage in sandstone arises from the combined effects of freeze-thaw cycles and loading， with the damage variable eventually approaching 1 as strain increases. The findings offer theoretical insights and experimental data for understanding the mechanical characteristics of rocks in cold regions.

Graphical abstract

关键词

冻融循环 / 核磁共振 / 孔隙度 / 损伤演化 / 耦合模型

Key words

freeze-thaw cycle / nuclear magnetic resonance / porosity / damage evolution / coupled model

我国寒区面积广泛，寒区工程岩体经历冻融循环的作用后，其内部产生孔隙水的相变、体积膨胀、应力集中等现象，从而导致岩体的损伤劣化和强度衰减，引发冻融剥蚀、冻融滑坡和隧道冻害等一系列工程灾害^［1-3］.因此，开展岩石冻融循环力学损伤机制和特性的研究，对于保障寒区岩体工程的可持续发展和长期稳定具有重要意义.

近年来，为了揭示冻融条件下岩石的损伤机制和力学特性，国内外学者在室内冻融实验和冻融理论方面开展了大量研究.在冻融实验研究方面：Song等^［4］介绍了核磁共振技术在页岩储层应用中的可行性.Gong等^［5］通过核磁共振和力学试验，揭示了冻融后花岗岩孔隙特征与力学性质.吴志军等^［6］利用低场核磁共振，通过T₂谱曲线和核磁共振图像等方法，揭示了裂隙岩体的渗透规律和机制.高峰等^［7］利用核磁共振技术测量了不同初始孔隙度的饱水砂岩的孔隙度变化量，建立了基于孔隙度变化量的相对剩余峰值强度劣化模型.Jiang等^［8］通过核磁共振和CT技术研究，表明岩石的微孔隙结构与冻融循环次数存在显著的相关性.这些研究表明，核磁共振技术可对岩石孔隙结构、孔径分布等物理特性进行快速、无损且准确的测量，在研究岩石冻融损伤中的孔隙特征和物理特性方面具有广泛的应用前景.通过核磁共振技术能够获得对岩石孔隙结构的深入理解，并进一步揭示冻融过程中孔隙特征与岩石力学性质损伤演化之间的关系.

在冻融理论研究方面：Yahaghi等^［9］建立了理论损伤模型，定量描述了砂岩性质的劣化规律，并进行了数值模拟，揭示了砂岩的损伤机制和破坏过程.朱谭谭等^［10］基于冻融循环试验，提出了一种基于颗粒流和颗粒膨胀的模拟方法，研究了轴向应力下岩石冻融变形和损伤规律.Yin等^［11］探讨了不同含水率下红砂岩的三轴蠕变特性，提出了冻融损伤和水分劣化模型，该模型准确描述了红砂岩的全阶段蠕变特性.寒区岩石通常受到温度、水和应力等因素的共同作用，在不同因素的耦合分析研究中，肖鹏等^［12］建立了一种考虑冻融与围压耦合作用的砂岩损伤本构模型，并通过三轴压缩试验验证了其可行性和准确性.王震等^［13］考虑了冻胀过程中水分迁移作用，建立了岩石力-热-水耦合弹塑性本构模型，发现了水分迁移等因素对冻融后的岩石损伤阈值等力学参数的影响.Liu等^［14］提出了一种简化的岩石冻融含相变THM（thermo-hydro-mechanical）耦合模型.刘泉声等^［15］运用双重孔隙介质模型等理论，推导出裂隙岩体在冻融条件下的THM耦合控制方程，研究了耦合条件下孔隙压力场和应力场等的分布规律.

目前关于岩石冻融循环的研究主要体现在对其宏观力学特性的影响规律方面，针对微观孔隙度及孔径分布变化的研究较少^［4-6］.同时现有岩石的损伤本构模型大多很难考虑微观孔隙度的变化，而是将其作为不变量，这与冻融循环作用下岩石损伤演化过程中孔隙度的实际变化不符.

因此，本文选取内蒙古扎哈淖尔露天煤矿南帮砂岩作为研究对象，探讨砂岩在冻融循环过程中以及冻融后受载荷影响的损伤行为.首先利用核磁共振技术得到冻融次数下砂岩样品的T2谱、孔隙度和孔径结构等参数，并基于力学试验，建立考虑损伤场与应力场、渗流场和温度场相互作用的多场耦合数学模型，开展冻融循环作用下砂岩力学特性演化数值模拟，分析细观结构演化特征与砂岩力学特性之间的关系.最后依据宏观损伤力学理论，揭示饱和砂岩在冻融循环和载荷作用下损伤演化规律和机制.

1 试验方法与结果分析

1.1 试样制备

试验样品取自内蒙古扎哈淖尔露天煤矿南帮砂岩岩芯，将岩芯用线切割机加工成25 mm\times 25 mm\times50 mm的立方体，采用双面磨石机将试样的端面磨平，试样高度误差控制在0.3 mm以内，两端面平整度控制在0.1 mm以内.试样的矿物成分见表1，初始平均物理参数见表2，加工后的试样选取完整性好、波速接近的岩块以减小试验的离散性.

1.2 试验流程

试验方案根据《工程岩体试验方案标准》（GB/T50266—2013）进行设计，通过冻融循环试验来模拟内蒙古寒冷地区岩石受季节变化和昼夜温差影响而产生的损伤现象.具体试验设备及流程如图1所示.

1） 岩样制取、烘干及饱和试验：通过电热风鼓风箱对筛选并分组的试样烘干，温度设定为105 ℃，时间为24 h，直至试样质量不再发生变化.烘干后，试样在室温状态冷却1 h后记录其干燥质量.将烘干后的试样放入真空饱水装置中进行饱水处理，先用真空泵抽真空2 h，然后注水再抽2 h，最后保持真空状态静置24 h，直至试样的质量不再发生变化.

2） 冻融循环流程：采用可程式恒温恒湿冻融试验机.考虑试样尺寸及当地冬季气温，制定如图2所示冻融循环方案：将冻结温度设定在-20 ℃，融化温度设定在20 ℃，冻结和融化过程的时间都为1 h，且到达设定温度时保持当前温度2 h，即一个冻融循环周期为6 h.最终设定0，10，20，30，40和50次的冻融循环，共6组，每组3个试件.

3） 低场核磁共振试验：采用MesoMR23-060H-I型低场核磁共振试验机对冻融前后的砂岩岩样开展核磁共振扫描试验，获得饱和砂岩岩样的孔隙度和T2谱，通过砂岩岩样的孔隙度和T2谱的分析，以研究砂岩的内部细观结构的变化趋势.

4） 冻融循环试验：将处理后的砂岩试样置于可程式恒温恒湿试验机中，第1次冻融循环周期先放入要冻融50次的砂岩试样，每完成10次冻融循环，对所有试样进行一次重新饱和处理，进行一次核磁共振试验测试，以分析冻融循环对砂岩孔隙结构的影响.

5） 单轴压缩试验：采用YAW-3000型微机控制电液伺服压力试验机，并采用PCI-Ⅱ声发射监测系统进行监测，加载速率为0.001 mm/s，待砂岩完成冻融循环试验和核磁共振试验后，分别对不同冻融循环次数的岩样进行单轴压缩试验，获得饱和砂岩岩样的宏观力学指标.

1.3 核磁共振测试

低场核磁共振是探测岩石内不同大小孔隙中的流体信号累积来揭示孔隙结构，通过对完全饱水岩石进行CPMG脉冲序列检测，对岩石孔隙中的水分子¹H含量进行检测，T₂弛豫时间谱揭示了孔隙尺寸与孔隙中的水信号之间的关系.T₂谱的不同峰值代表了不同尺寸孔隙中的流体量，而峰值大小和面积则表明了相应孔隙大小的分布和数量.通过这种方式量化岩石内部孔隙尺寸和分布的信息.

对砂岩进行核磁共振测试，分别获得了岩石0，10，20，30，40和50次的核磁共振T₂谱，以B1，B2岩样进行分析.

如图3所示，完全饱水环境下砂岩的T₂谱呈双峰分布，第1个波的范围在0.1~75 ms，第2个波的范围在75~1 000 ms.随着冻融循环次数的增加，砂岩岩样T₂谱的幅值及峰面积随着冻融循环的作用逐渐增大，说明岩石整体孔隙度随着冻融循环次数的增加而增大.

对于岩石孔隙中存在的流体，其弛豫行为可分为3种主要类型：表面弛豫、扩散弛豫和体积弛豫.当岩石孔隙中的流体是单一的，表面弛豫的强度显著大于扩散弛豫和体积弛豫.因此，通常将核磁共振中观测到的整体横向弛豫速率简化为表面弛豫效应，横向弛豫时间与孔径大小关系如下：

\frac{\mathrm{1}}{T_{\mathrm{2}}}={\rho }_{\mathrm{2}}\left(\frac{S}{V}\right)={\rho }_{\mathrm{2}}\left(\frac{F_{\mathrm{s}}}{R}\right).(1)

式中：V为岩石孔隙体积，cm³；S为岩石孔隙表面积，cm²；{\rho }_{\mathrm{2}}为横向表面弛豫强度，μm/ms；R为岩石孔隙半径；F_{\mathrm{s}}为岩石孔隙的形状因子，与岩石的高径比有关，本文F_{\mathrm{s}}为2.计算得到T₂谱分布相对应的岩石孔径，如图4所示.

将砂岩的孔隙按照孔径的大小分为3类^［16］：小孔（r≤0.1 μm）、中孔（0.1 μm<r<1μm）、大孔（r≥1 μm）.由图4可见，砂岩孔径分布主要集中在0.1 μm以下，说明砂岩内部结构中小孔孔隙占主导，随着冻融循环次数的增加，第1波峰的强度降低，说明岩石中的小孔数量在减少，第2波峰的强度上升，意味着岩石中的小孔在冻融作用下逐渐贯穿形成中孔或大孔.砂岩内原有小尺寸孔隙逐渐扩展，小孔占比随之减少，中孔和大孔分布比例逐渐上升，说明冻胀作用对岩石孔径增加起主导作用，岩石内部孔隙度逐渐增大.

1.4 砂岩宏观强度

声发射能够有效地揭示岩石破裂失稳过程中损伤演化的特性，通过单轴压缩试验过程中声发射监测，研究载荷和冻融循环对砂岩的孔隙结构及其力学性质的影响.在试验中采用了PCI-Ⅱ声发射监测系统，传感器均匀分布在样品表面，每个样品上布置3个传感器，以确保全面覆盖声发射信号.图5展示了通过冻融循环后的单轴压缩试验，砂岩经历0，10，20，30，40和50次冻融循环后的声发射单轴压缩曲线.

由图5可见，在砂岩加载初期，声发射信号相对较少.随着应力增大，砂岩内部孔隙开始扩张，微小裂缝数量逐渐增加，进而声发射信号显著增加.岩石在失稳时，砂岩内部颗粒因持续摩擦而产生更多声发射信号.可以观察到随着冻融次数的增加，岩石达到应力峰值的时间更长，冻融导致岩石内部孔隙破坏，降低了岩石脆性，使其在外力作用下表现出明显的延性特征.岩石在经历30，40，50次冻融循环后，可以明显地观察到强度峰值附近曲线出现振荡现象，这种现象是由于冻融循环次数增加导致砂岩内部结构不稳定，部分区域微裂缝扩展贯通并最终断裂，岩石在受压状态下这些区域局部破坏，但在整体上仍保持一定的承压能力，从而引起宏观力学行为的非连续性.此外，冻融作用下的孔隙结构破坏也会影响声发射的累计振铃计数，随着冻融次数的增加，累计振铃计数减少.这一变化表明，冻融作用下岩石的损伤是不可逆的，损伤程度随冻融次数增加而加剧.

图6为不同冻融次数砂岩单轴压缩应力-应变曲线.冻融环境下，砂岩力学特性与其细观结构参数变化特征基本一致，冻融次数越多，岩石内部孔隙越发育，中孔和大孔的比例上升，其承压强度也越低.具体而言，冻融循环次数增加，其单轴抗压强度下降，应变增大，砂岩由脆性破坏向延性破坏转变.经过10，20，30，40，50次冻融循环，砂岩试样单轴压缩强度分别降低了6.46%，14.45%，19.39%，36.88%，47.91%；试样弹性模量分别降低了6.35%，16.01%，23.64%，32.86%，51.08%.这表明在冻融循环过程中，岩样的孔隙和微裂纹不断发育并且扩展，导致其损伤进一步加剧.

岩石单轴压缩曲线可以分为4个阶段：压密阶段，这个阶段应力较低，岩石内部的初始微裂纹被压密；线弹性阶段，随着应力的增加，岩石的应力、应变之间存在着线性关系；裂纹扩展阶段，当应力到达一定程度时，岩石内部的微裂纹开始扩展，导致岩石的强度降低；失稳破裂阶段，微裂纹扩展并且贯穿形成宏观裂纹，导致岩石最终发生破坏.

图7展示了岩样不同冻融循环次数的单轴抗压强度.通过拟合可以得到砂岩冻融循环与单轴抗压强度的关系式为

\sigma =\mathrm{61.97}-\mathrm{11.6}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{N}{\mathrm{44.1}}\right).(2)

式中：\sigma为抗压强度，MPa；N为冻融循环次数.

随着冻融次数增加，砂岩单轴抗压强度呈现下降趋势，且下降速率逐渐上升，表明岩石中孔隙尺寸与数量的变化是导致砂岩强度降低的本质原因.

2 岩石冻融损伤多场耦合模型

岩石在冻融作用下的损伤过程是一个复杂的多场耦合过程，受到应力场、渗流场和温度场的共同影响，这些场之间的相互作用直接影响着岩石的损伤演化过程.然而，现有研究多集中于损伤场与3个物理场中的某一个或两个场的耦合关系，忽略了它们之间的多因素正交作用.因此，建立一个综合考虑损伤场、应力场、渗流场和温度场多场耦合相互作用的数学模型，对于深入理解和准确预测岩石冻融损伤过程具有重要的意义.

非均质性是岩石的本质特性.引入Weibull^［17］统计分布来表征岩石材料细观结构的非均质特性，将岩石材料细观弹性模量、渗透率等力学参数进行Weibull统计分布赋值：

\phi (i)=\frac{m}{i_{\mathrm{0}}}{\left(\frac{i}{i_{\mathrm{0}}}\right)}^{m-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\left(\frac{i}{i_{\mathrm{0}}}\right)}^m\right].(3)

式中：\phi (i)为岩石材料细观力学参数的分布函数；m为均质度系数，反映材料的均质性，m的值越大表示材料单元力学性质越均匀；i为岩石材料细观力学参数；i_{\mathrm{0}}为所有岩石材料细观结构参数的近似统计平均值.

2.1 应力-应变场及损伤演化方程

对于孔隙岩石材料的控制方程为

{\mathit{\sigma }}_{ij}^{\mathrm{\text{'}}}={\mathit{\sigma }}_{ij}-b\overline{P}{\delta }_{ij}.(4)

式中：{\mathit{\sigma }}_{ij}^{\mathrm{\text{'}}}为岩石的有效应力张量，i,j=\mathrm{1,2},\mathrm{3}；{\mathit{\sigma }}_{ij}为岩石所受应力张量；\overline{P}为平均孔隙压力；b为Biot系数，与多孔体系的体积弹性模量有关；{\delta }_{ij}为克罗内克符号.

岩石材料结构单元的静力平衡微分方程为

{\mathit{\sigma }}_{ij}^{\mathrm{\text{'}}}+F_i=\mathrm{0}.(5)

式中F_i为净体力在i(x,y,z)方向上的分量.

岩石材料的总应变由弹性应变和由温度变化引起的应变构成：

{\mathit{\epsilon }}_{ij}={\mathit{\epsilon }}_{ij}^{\mathrm{e}}+{\mathit{\epsilon }}_{ij}^{\mathrm{h}}.(6)

式中：{\mathit{\epsilon }}_{ij}为应变张量；{\mathit{\epsilon }}_{ij}^{\mathrm{e}}为弹性应变张量；{\mathit{\epsilon }}_{ij}^{\mathrm{h}}为温度应变张量.

岩石所受有效应力与应变的关系为

\begin{array}{cc}{\mathit{\sigma }}_{ij}^{\mathrm{\text{'}}}=\mathit{H}{\mathit{\epsilon }}_{ij}^{\mathrm{e}}=\mathit{H}({\epsilon }_{ij}-{\epsilon }_{ij}^{\mathrm{h}}),& {\mathit{\epsilon }}_{ij}^{\mathrm{h}}={\alpha }_{\mathrm{t}}(T-T_{\mathrm{f}})\end{array}{\delta }_{ij}.(7)

式中：\mathit{H}为弹性刚度矩阵；T为实际温度；{\alpha }_{\mathrm{t}}为热膨胀系数；T_{\mathrm{f}}为温度应变为0时的临界温度.

同时，应变与位移的关系为

{\mathit{\epsilon }}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}(u_{i,j}+u_{j,i}).(8)

假设岩石满足弹性损伤理论，岩石细观单元弹性模量随损伤程度的变化如下：

E=E_{\mathrm{0}}(\mathrm{1}-D).(9)

式中：D为损伤变量；E_{\mathrm{0}}，E分别表示单元损伤前和损伤后的弹性模量.

当岩石材料结构单元的应力或应变状态到达一个损伤阈值，单元产生损伤.系统首先从最大拉应力准则进行判断，当单元的应力状态满足最大拉应力准则时，单元发生拉伸破坏；当单元的剪应力满足Mohr-Coulomb准则的损伤阈值时，单元发生剪切破坏.

f_{\mathrm{1}}=-{\sigma }_{\mathrm{3}}-{\sigma }_{\mathrm{t}\mathrm{0}},f_{\mathrm{2}}={\sigma }_{\mathrm{1}}-\frac{{\sigma }_{\mathrm{3}}(\mathrm{1}+\mathrm{s}\mathrm{i}\mathrm{n}\phi )}{(\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi )}-{\sigma }_{\mathrm{c}\mathrm{0}}.(10)

式中：{\sigma }_{\mathrm{t}\mathrm{0}}和{\sigma }_{\mathrm{c}\mathrm{0}}分别为单元的单轴抗拉强度和单轴抗压强度；{\sigma }_{\mathrm{1}}和{\sigma }_{\mathrm{3}}分别为单元的最大和最小主应力.当f_{\mathrm{1}}=\mathrm{0}或f_{\mathrm{2}}=\mathrm{0}时，材料单元分别发生拉伸损伤和剪切损伤.当f_{\mathrm{1}}>\mathrm{0}和f_{\mathrm{2}}>\mathrm{0}，表示材料单元会继续受损，即损伤将持续发展.相反，若和，则表示损伤保持不变，保持与上一个时间步相同的数值，而不会发生新的损伤.材料单元损伤变量D如式（11）和（12）所示.

细观单元发生拉伸损伤：

(11)

细观单元发生剪切损伤：

(12)

式中：{\epsilon }_{\mathrm{t}\mathrm{u}}为最大拉伸主应变；{\sigma }_{\mathrm{t}\mathrm{r}}和{\sigma }_{\mathrm{c}\mathrm{r}}分别为材料单元的拉伸损伤残余强度和剪切损伤残余强度；{\epsilon }_{\mathrm{t}\mathrm{0}}和{\epsilon }_{\mathrm{c}\mathrm{0}}分别对应材料单元发生拉伸破坏和压缩破坏的最大拉应变和最大压应变.

2.2 渗流场及渗透率损伤演化方程

根据孔隙结冰理论和静水压力理论的假设，基于达西定律和流体质量守恒方程，平衡方程的一般表达式可写成^［18］

\gamma \dot{p}=\mathrm{d}\mathrm{i}\mathrm{v}\left(\frac{K}{\eta }\nabla p_{\mathrm{w}}\right)+Q_{\mathrm{m}}-b{\dot{\epsilon }}_{\mathrm{v}}，(13)

\gamma =\frac{\mathrm{1}-n}{R_{\mathrm{m}}}+\frac{n{S}_{\mathrm{w}}}{R_{\mathrm{w}}}+\frac{n{S}_{\mathrm{i}}}{R_{\mathrm{i}}}，(14)

Q_{\mathrm{m}}=\left(\frac{\mathrm{1}}{{\rho }_{\mathrm{i}}}-\frac{\mathrm{1}}{{\rho }_{\mathrm{w}}}\right){\dot{w}}_{\mathrm{i}}+\overline{\alpha }\dot{T}.(15)

式中：\dot{p}为孔隙水压力；\eta为水的动力黏滞系数；Q_{\mathrm{m}}为质量源项；n为孔隙度；K为渗透率；{\epsilon }_{\mathrm{v}}为体应变；R_{\mathrm{m}}，R_{\mathrm{w}}和R_{\mathrm{i}}分别为岩石、水、冰的压缩模量；S_{\mathrm{w}}和S_{\mathrm{i}}分别为水和冰所占孔隙的比例；{\rho }_{\mathrm{i}}和{\rho }_{\mathrm{w}}分别是冰和水的密度；\overline{\alpha }为材料的体膨胀系数.

式（13）描述了岩石内部的渗流现象，并考虑了水和冰的占比变化以及热膨胀效应对渗流场的影响.通过求解这些方程，可以获得岩石内部水分渗流的变化情况，从而进一步研究损伤等因素对渗流场的影响.岩石在冻融过程中，产生的冻胀力和水压力会改变岩石的内部结构，即渗透率和孔隙度都会随着冻融次数的增加而改变，岩石的渗透率和孔隙度之间满足如下幂函数关系式^［19］：

n=n_{\mathrm{0}}{\left(K/K_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{3}}.(16)

式中：n_{\mathrm{0}}为初始孔隙度；K_{\mathrm{0}}为初始渗透率.

在岩石经历弹性变形过程中，渗透率会随着岩石的压实和膨胀而变化，在岩石压缩和拉伸过程中，可以使用以下公式^［20］描述岩石渗透率的变化，岩石在压缩状态下时：

K=\left\{\begin{array}{c}\begin{array}{l}\begin{array}{cc}{K}_{\mathrm{0}}{\mathrm{e}}^{-\beta \left({\sigma }_{\mathrm{1}}-\alpha p\right)},& D=\mathrm{0}\end{array},\\ \begin{array}{cc}\xi K_{\mathrm{0}}{\mathrm{e}}^{-\beta \left({\sigma }_{\mathrm{1}}-\alpha p\right)},& D\end{array}>\mathrm{0},\end{array}\end{array}\right.(17)

细观单元在拉伸状态下，渗透率与岩石损伤程度的关系式为

(18)

式中：\xi为加载单元被压缩时的渗透率变化系数；{\xi }^{\mathrm{\text{'}}}为加载单元受拉伸状态时的渗透率变化系数；\beta为孔隙压力和应力的耦合系数；\alpha为孔隙压力系数.

2.3 温度场及热传导损伤演化方程

在冻融循环过程中，岩石内部水和冰的含量时刻在变化，考虑水热相变潜热，岩石冻融循环中的温度场方程可以通过热传导定律来描述，假设岩石内部的固体和流体之间保持能量守恒及热平衡，即

\left(\rho c_{\mathrm{p}}\right)\frac{\partial T}{\partial t}=\rho c_{\mathrm{p}}{\gamma }_{\mathrm{u}}\nabla T+\nabla \left(l_{\mathrm{P}}\cdot \nabla T\right)+Q_{\mathrm{T}}，(19)

c_{\mathrm{P}}=n{S}_{\mathrm{w}}{c}_{\mathrm{w}}+n{S}_{\mathrm{i}}{c}_{\mathrm{i}}+\left(\mathrm{1}-n\right)c_{\mathrm{m}}，(20)

l_{\mathrm{P}}=n{S}_{\mathrm{w}}{l}_{\mathrm{w}}+n{S}_{\mathrm{w}}{l}_{\mathrm{w}}+\left(\mathrm{1}-n\right)l_{\mathrm{m}}.(21)

式中：\rho为岩石材料整体密度；Q_{\mathrm{T}}为热源项；t为时间；c_{\mathrm{p}}，c_{\mathrm{i}}，c_{\mathrm{w}}，c_{\mathrm{m}}分别为材料整体、冰、水、岩石的比热容；l_{\mathrm{p}}，l_{\mathrm{i}}，l_{\mathrm{w}}，l_{\mathrm{m}}分别为材料整体、冰、水、岩石的导热系数.

损伤必然会导致岩体的导热系数增加.关于损伤对导热系数和比热容的影响，假设损伤对导热系数的影响符合指数关系^［21］，即

l_{\mathrm{m}}\left(T,D\right)=l_{\mathrm{m}}\left(T\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(D/{\alpha }_{\mathrm{l}}\right).(22)

式中{\alpha }_{\mathrm{l}}表示损伤对导热性能的影响程度.

综上，方程（9），（11），（12），（17），（18）和（22）描述了岩石冻融循环损伤的演化过程，方程（4），（13）和（19）构成了耦合的非线性方程组，共同构建了渗流-应力-温度-损伤的控制方程.多场耦合关系示意图如图8所示.

3 冻融循环作用下砂岩损伤演化模拟

3.1 数值模型

本文基于Comsol Multiphysics软件的二次开发平台，将建立的耦合模型写入程序，并通过MATLAB软件实现，对砂岩冻融循环损伤过程作数值模拟分析，模拟中的模型尺寸与实际砂岩试样保持一致.模型各表面温度按照图1温度变化所设置.模型底面固定约束，其他面自由.使用的参数指标如表3所示，模型示意图如图9所示.

实验设置每0.5 h为一步，每完成一步计算后得到岩石的强度、弹性模量、渗透率、孔隙度和损伤等参数.将弹性模量等参数导入MATLAB进行损伤评估.若判定参数未达到损伤阈值就不记录数据，相反，若判定参数达到或超过设定的损伤阈值，记录损伤，同时对其他参数进行折减，得到的参数作为下一次计算的初始值从而实现参数的循环迭代更新.

3.2 砂岩孔隙度分析

在核磁共振测量中，¹H与水分子的关系使岩石孔隙中的水含量与核磁信号强度成正比.通过核磁共振试验可以得到岩石的孔隙度，饱水岩石的核磁信号越强，表示岩石的孔隙度越高.砂岩内部孔隙发育，可以达到完全饱水状态，所以同样可以用称重法来计算砂岩的孔隙度.根据构建的渗流-应力-温度-损伤的岩石冻融循环损伤模型，通过MATLAB软件，对B2号砂岩的冻融循环损伤过程进行数值模拟计算，得到砂岩的孔隙度，数值模拟的初始孔隙度取称重法和核磁共振法测得的平均值16.15%.

用3种方法计算得出砂岩不同冻融循环次数的孔隙度如图10所示，砂岩的孔隙度发展的趋势接近一致，呈上升趋势，随着冻融循环次数的增加而稳定升高.岩石冻结时，液态水转变为固态，引发冻胀现象使孔隙体积扩大.融化过程中，固态水转为液态，岩石因孔隙尺寸扩大而吸收更多水分，冻融循环导致孔隙度逐渐增加，从而使试件的饱水质量逐步上升.

同时从图10中可以看出，冻融初期孔隙度上升缓慢，冻融0~20次，增幅较小；随着冻融次数的增加，孔隙度增大的速率变快，这是因为在冻融初始阶段，砂岩岩样内部的骨架结构紧密，冻融作用后使得原本紧密的岩样骨架结构开始松散，颗粒间的黏结力减弱，随着颗粒间隙的增大，颗粒之间的相互作用也增强，孔隙逐步发育，砂岩内部的裂纹逐渐扩展并贯穿，裂纹尺寸随骨架的破裂而增长，使得岩石表面经过初期冻融后与初始状态显著不同.当冻融40~50次时，孔隙度的增大速率有所下降，这说明孔隙的发育不是无限增大的.

综上，数值模拟与试验结果的对比显示，本研究提出的岩石冻融循环过程中渗流-应力-温度-损伤耦合的数值模型能有效表征砂岩孔隙度和损伤的变化特征.

3.3 砂岩损伤模型验证

岩石的冻融循环是一个持续较长时间的过程，相比之下，岩石的单轴压缩发生在较短的时间内，所以需要分别模拟计算岩石冻融损伤和受载荷损伤.岩石冻融和受载荷产生的损伤是不可逆的，即二者的损伤量可以叠加，将经过冻融处理的岩石模拟数据如损伤量、弹性模量等参数代入单轴压缩程序中作为计算的初始值.在模拟之前先要处理原有单轴压缩应力-应变曲线^［22］，数据处理后最终得到结果如图11所示.

由图11可以看出，砂岩在弹性阶段的应力-应变曲线实验数据与数值模拟结果高度一致，二者的峰值强度也非常接近.然而，在经历40次和50次冻融循环后，试验的应力-应变曲线在强度峰值附近出现了明显的振荡现象.这一现象可以解释为：随着冻融循环次数的增加，砂岩内部结构逐渐失稳，局部区域出现微裂缝的形成与扩展，最终导致断裂.在受压过程中，这些区域首先发生局部破坏，但岩石整体仍保持一定的承载能力，随着外部载荷的持续作用，应力水平进一步上升，导致宏观力学行为呈现非连续性.然而，数值模拟并未捕捉到这一现象，因此未能反映出应力的上升和下降过程.总体而言，图11进一步验证了砂岩渗流-应力-温度-损伤耦合数值模型的合理性.

3.4 砂岩冻融受载荷损伤演化特征

岩石的冻融受载荷损伤可以通过微观和宏观两个层面来研究，岩石遭受的冻融和受载荷损伤导致其微观结构和力学性能出现明显劣化.根据宏观唯像损伤力学理论，可以通过宏观物理参数来揭示岩石的损伤和劣化程度.岩石的弹性模量在冻融循环周期中容易测定和分析，本研究采用弹性模量作为冻融损伤变量：

D_n=\mathrm{1}-\frac{E_n}{E_{\mathrm{0}}}.(23)

式中：E_{\mathrm{0}}为冻融前岩石的弹性模量；E_n为不同冻融次数下岩石的弹性模量.

岩石内部的微观结构是非均质性的，岩石在受外部载荷作用的过程中，为了从微观角度表征岩石微元体强度的分布，由式（3）基于Weibull分布函数，岩石中基本单元参数i和i_{\mathrm{0}}分别替换为岩石的应变\epsilon和应力-应变关系曲线应力峰值对应的应变值{\epsilon }_{\mathrm{0}}，因此，可以推导出受载荷作用下岩石的损伤演化方程：

D_{\mathrm{p}}={\int }_{\mathrm{0}}^{\epsilon }P(\epsilon )d\epsilon =\mathrm{1}-exp\left[-\frac{\mathrm{1}}{m}\left({(\frac{\epsilon }{{\epsilon }_{\mathrm{0}}})}^m\right)\right].(24)

式中，P(\epsilon )是岩石在受载荷作用中的强度分布函数，可以衡量岩石微元体的损伤率，还能反映微元体破坏的累计效应，这种累计效应最终导致岩石的宏观强度特性劣化.

冻融后的岩石受载荷作用下的本构方程^［23］可表示为

D_{\mathrm{m}}=D_{\mathrm{n}}+D_{\mathrm{p}}-D_{\mathrm{n}}{D}_{\mathrm{p}}.(25)

式中，D_{\mathrm{m}}为岩石在冻融和受载荷双重作用下的总损伤量，冻融与载荷之间的耦合作用实际上对总体损伤的程度产生了一定的减弱效果.

由式（23），（24），（25）可推导出冻融循环次数和应变为控制变量的砂岩冻融与载荷耦合的总损伤演化方程为

D_{\mathrm{m}}=\mathrm{1}-\frac{E_n}{E_{\mathrm{0}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{\mathrm{1}}{m}\left({\left(\frac{\epsilon }{{\epsilon }_{\mathrm{0}}}\right)}^m\right)\right].(26)

式（26）表明，在只考虑岩石受载荷作用造成的损伤时，E_{\mathrm{0}}=E_n，即D_{\mathrm{m}}=D_{\mathrm{P}}；而在仅考虑冻融损伤的情况下，\epsilon =\mathrm{0}，即D_{\mathrm{m}}=D_n.将不同冻融循环次数砂岩实验和模拟得到弹性模量和应变等参数代入式（26），可以计算出砂岩冻融受载荷作用下的损伤演化曲线.如图12所示，两组曲线都呈现一个先缓慢增加、再快速增长、最后趋于平缓的过程，当\epsilon =\mathrm{0}时，分别对应不同冻融循环次数下岩石的弹性模量，随着冻融循环次数的增加，损伤程度逐渐增大.

损伤演化曲线分为4个阶段：在岩石加载初期，岩样内部的微孔隙和微裂纹逐渐被压密，岩石表现出抗压性，强度提高，砂岩的损伤劣化变量初期增长较慢；随着应变增加，岩石抗压性逐渐降低，岩石损伤演化曲线呈线性趋势；当载荷持续施加时，岩体内部原有裂纹会逐渐闭合，促使岩石恢复一部分抗压性，因此，岩石损伤演化曲线斜率逐渐降低；在加载过程后期，无论是由冻融还是荷载引起的损伤，损伤变量增加的速度都会逐渐减缓，岩石内部微裂纹逐渐贯穿形成宏观裂纹，岩石损伤程度最终趋于1，即完全破坏.这个演化过程解释了在不同阶段砂岩的结构演化特征，说明砂岩冻融受载荷损伤演化曲线可以表征砂岩宏观力学特性.

4 结 论

1） 建立了一种能够描述岩石冻融循环过程中温度、应力和渗流与损伤耦合作用的本构模型，通过数值模拟计算出了砂岩孔隙度随冻融次数增加的变化趋势，冻融初期孔隙度上升缓慢，随后孔隙度上升速率变快，这一结果与称重法和核磁共振法测得的孔隙度的趋势基本吻合.此模型为深入探讨岩石冻融损伤机理提供了理论基础.

2） 砂岩的核磁共振T₂谱分布呈现出明显的“双峰”趋势，随着冻融次数的增加，第1峰值逐渐减小，第2峰值随之增大，直观地反映了岩石内部微观结构变化.砂岩的孔隙按照孔径的大小分为微孔、中孔和大孔.随着冻融次数的增加，微孔逐渐扩展为中孔和大孔，从而导致岩石的力学性能出现不同程度的削弱.

3） 冻融受载荷作用下砂岩的损伤演化模型可以用来揭示细观损伤与宏观力学参数之间的关系，两种因素的共同作用使砂岩的总损伤加剧，并且随着应变的增加，砂岩的损伤变量最终趋于1，有效展现了岩石的损伤程度，揭示了岩石冻融受载荷的损伤机理.

参考文献

[1] Lyu Z T， Xia C C， Liu W D. Analytical solution of frost heaving force and stress distribution in cold region tunnels under non-axisymmetric stress and transversely isotropic frost heave of surrounding rock ［J］. Cold Regions Science and Technology， 2020， 178：103117.

[2] Li Z G， Xu T， Zhao L C， et al. Enhancing stability analysis of open-pit slopes via integrated 3D numerical modeling and data monitoring ［J］. Engineering Failure Analysis， 2024， 163： 108495.

[3] 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.

[4] Song Y Q， Kausir R. NMR application in unconventional shale reservoirs-a new porous media research frontier ［J］. Progress in Nuclear Magnetic Resonance Spectroscopy， 2019， 112：17-33.

[5] Gong Y F， Song J X， Wu S Z， et al. Evolution of pore structure and analysis of freeze damage in granite during cyclic freeze-thaw using NMR technique［J］. Engineering Geology， 2024， 335： 107545.

[6] 吴志军， 卢槐， 翁磊， 等. 基于核磁共振实时成像技术的裂隙砂岩渗流特性研究 ［J］. 岩石力学与工程学报， 2021， 40（2）： 263-275.

[7] Wu Zhi-jun， Lu Huai， Weng Lei， et al. Investigations on the seepage characteristics of fractured sandstone based on NMR real-time imaging ［J］. Chinese Journal of Rock Mechanics and Engineering， 2021，40（2）：263–275.

[8] 高峰， 熊信， 周科平， 等. 冻融循环作用下饱水砂岩的强度劣化模型 ［J］. 岩土力学， 2019， 40（3）： 926-932.

[9] Gao Feng， Xiong Xin， Zhou Ke-ping， et al. Strength deterioration model of saturated sandstone under freeze-thaw cycle ［J］. Rock and Soil Mechanics， 2019， 40（3）： 926-932.

[10] Jiang H B， Li K N， Jin J. The variation characteristics of micro-pore structures of underground rocks in cold regions subject to freezing and thawing cycles ［J］. Arabian Journal of Geosciences， 2020， 13：1-7.

[11] Yanaghi J， Liu H Y， Chan A， et al. Experimental， theoretical and numerical modelling of the deterioration and failure process of sandstones subject to freeze-thaw cycles ［J］. Engineering Failure Analysis， 2022， 141（1）： 106686.

[12] 朱谭谭， 李昂， 黄达， 等. 应力-冻融耦合作用下砂岩变形与损伤特征研究［J］. 岩石力学与工程学报，2023，42（2）： 342-351.

[13] Zhu Tan-tan， Li Ang， Huang Da， et al. Deformation and damage characteristics of sandstone under the combined action of stress and freeze-thaw cycle ［J］. Chinese Journal of Rock Mechanics and Engineering， 2023， 42（2）： 342-351.

[14] Yin W， Wang X Y， Zheng S A， et al. Triaxial creep test and damage model study of layered red sandstone under freeze-thaw cycles［J］. Case Studies in Construction Materials， 2024， 21： e03785.

[15] 肖鹏， 陈有亮， 杜曦， 等. 冻融循环作用下砂岩的力学特性及细观损伤本构模型研究 ［J］. 岩土工程学报， 2023， 45（4）： 805-815.

[16] Xiao Peng， Chen You-liang， Du Xi， et al. Mechanical properties of sandstone under freeze-thaw cycles and the study of meso-damage constitutive model ［J］. Chinese Journal of Geotechnical Engineering， 2023，45（4）： 805-815.

[17] 王震， 朱珍德， 陈会官， 等. 冻融作用下岩石力-热-水耦合本构模型研究 ［J］. 岩土力学， 2019， 40（7）： 2608-2616.

[18] Wang Zhen， Zhu Zhen-de， Chen Hui-guan， et al. A thermo-hydro-mechanical coupled constitutive model for rocks under freeze-thaw cycles ［J］. Rock and Soil Mechanics， 2019， 40（7）： 2608-2616.

[19] Liu N F， Liang S H， Wang S J， et al. THM model of rock tunnels in cold regions and numerical simulation［J］. Scientific Reports， 2024， 14（1）： 3465.

[20] 刘泉声， 康永水， 刘滨， 等. 裂隙岩体水-冰相变及低温温度场-渗流场-应力场耦合研究 ［J］. 岩石力学与工程学报， 2011， 30（11）： 2181-2188.

[21] Liu Quan-sheng， Kang Yong-shui， Liu Bin， et al. Water-ice phase transition and thermo-hydro-mechanical coupling at low temperature in fractured rock［J］. Chinese Journal of Rock Mechanics and Engineering， 2011， 30（11）： 2181-2188.

[22] Nie R S， Zhou J， Chen Z X， et al. Pore structure characterization of tight sandstones via a novel integrated method： a case study of the Sulige gas field， Ordos Basin （Northern China） ［J］. Journal of Asian Earth Sciences， 2021， 213： 104739.

[23] Xu B， Xu T， Xue Y C， et al. Phase field modeling of mixed-mode crack in rocks incorporating heterogeneity and frictional damage ［J］. Engineering Fracture Mechanics， 2024， 298： 109936.

[24] Zuber B， Marchand J. Predicting the volume instability of hydrated cement systems upon freezing using poro-mechanics and local phase equilibria ［J］. Materials and Structures，2004， 37： 257-270.

[25] Zhu W C， Wei C H， Liu J， et al. A model of coal-gas interaction under variable temperatures ［J］. International Journal of Coal Geology， 2011， 86（2/3）： 213-221.

[26] Xu T， Tang C A， Yang T H， et al. Numerical investigation of coal and gas outbursts in underground collieries ［J］. International Journal of Rock Mechanics and Mining Sciences， 2006， 43（6）： 905-919.

[27] 朱万成， 魏晨慧， 田军， 等. 岩石损伤过程中的热-流-力耦合模型及其应用初探 ［J］. 岩土力学， 2009，30（12）： 3851-3857.

[28] Zhu Wan-cheng， Wei Chen-hui， Tian Jun， et al. Coupled thermal-hydraulic-mechanical model during rock damage and its preliminary application ［J］. Rock and Soil Mechanics， 2009， 30（12）： 3851-3857.

[29] Hao S W， Wang H Y， Xia M F， et al. Relationship between strain localization and catastrophic rupture ［J］. Theoretical and Applied Fracture Mechanics， 2007， 48（1）： 41-49.

[30] Huang S B， Liu Q S， Cheng A P， et al. A statistical damage constitutive model under freeze-thaw and loading for rock and its engineering application［J］. Cold Regions Science and Technology， 2018， 145： 142-150.

基金资助

国家自然科学基金资助项目(42172312)

国家自然科学基金资助项目(51950410595)

科技部外国专家重点支撑计划项目(ZDKY2023077L)