寒区工程岩体处于恶劣的环境及复杂的应力状态下，通常会伴随冻融和长期荷载的共同作用。在寒区露天岩体工程施工过程中则表现为冻融、长期加载和卸载的综合作用，为保证工程建设的安全运营，必须同时考虑冻融与加卸载蠕变共同作用引起的时效损伤对工程岩体强度的影响。

目前，在岩石力学特性和冻融环境下岩石破坏机理的研究方面，国内外学者都有诸多重要贡献。例如，Liu等^[1]对红砂岩进行室内单轴抗压强度和纵波速度测试，探讨了冻融环境下含水饱和度对岩石力学性质的影响；Shi等^[2]通过冻融循环和三轴压缩试验，分析了冻融次数和围压对红砂岩力学性能的影响；Huang等^[3]将岩石的冻融损伤用弹性模量来表示，建立了冻融和荷载共同作用下的损伤本构模型，并将其应用于寒冷地区热-水-机械耦合条件下隧道的稳定性分析；张峰瑞等^[4]建立了花岗岩的冻融剪切蠕变本构模型，通过对蠕变参数的敏感性分析，研究了蠕变参数对花岗岩蠕变变形的影响。上述工作主要是对冻融岩石物理力学特性的研究，但从本质上看，岩石的宏观形变和破坏是许多细观、晶粒尺度的裂缝贯通和孔洞坍塌造成的^{[5 − 6]}。为更加清晰地认识冻融岩石蠕变特性，必须分析其细观孔隙结构，从细观层面来阐释冻融岩石的蠕变特性。

核磁共振（nuclear magnetic resonance，NMR）作为一种无损测试技术，是测定岩体内部结构演化的主要方法之一，可系统获得流体在样品中赋存及迁移的相关参数，例如渗透率、孔隙度、T₂谱分布（T₂为横向弛豫时间）等。目前这一技术已被广泛用于冻融岩石细观结构评估^[7]，周科平等^[8]、李杰林等^[9]对不同冻融岩石进行核磁共振测试和成像分析，对冻融岩石的细观损伤与宏观力学特性之间的联系进行了研究，得出了岩石在冻融作用下孔隙结构的演化及扩展规律；陈国庆等^[10]对不同温差冻融后的石英砂岩、变质砂岩及红砂岩进行核磁共振检测，探讨了各类砂岩在不同温差冻融作用下的孔隙结构发展规律。以上研究为冻融岩石细观结构表征积累了丰富的资料，但将核磁共振技术应用于冻融岩石加卸载蠕变试验中，从细观角度描述弹性后效特征与蠕变力学特性联系的研究却十分罕见。岩石内部结构具有较强的分形特征^{[11 − 13]}，因此，利用分形理论可以定量描述岩石细观结构的复杂性^[14]。另外，利用分形维数来描述加卸载蠕变过程中冻融岩石孔隙结构复杂程度的变化，探讨其与细观参数及宏观特性之间关系的研究同样罕见。

基于此，本文通过对不同冻融次数的红砂岩岩样进行分级加卸载蠕变试验，并配合核磁共振检测，同时引入分形理论，探究冻融后红砂岩在蠕变过程中细观参数、孔隙结构复杂程度与蠕变力学特性之间的关系，并建立了考虑冻融对岩石蠕变损伤影响的本构模型。

1.   试验方法

1.1   试样制备

试验所用岩石采自陕北某工程岩体，是以长石、石英、蒙脱石等成分为主的中细粒红砂岩。首先，按国际岩石力学学会试验规程^[15]将岩块加工成直径50 mm、高度100 mm、端面不平整度不超过0.02 mm的圆柱体（图1）；然后烘干岩样，待其冷却后测得其干密度和纵波波速；最后将岩样饱和，测得其饱和质量，进而获得岩样的饱和密度、饱和含水率和孔隙度。岩样物理性质具有良好的相似性，表1为岩样基本物理参数的平均值。依据纵波波速和干密度相近的原则，将岩样分为Ⅰ、Ⅱ两组（表2），每组3个岩样对应不同的冻融次数。将组Ⅰ岩样用于抗压强度试验，组Ⅱ岩样用于蠕变试验。

图  1  红砂岩岩样尺寸及安装示意

Figure  1.  Rock sample size and installation instructions

下载: 全尺寸图片 幻灯片

表  1  各项物理参数的平均值

Table  1.  Mean physical parameters of the rock samples

物理参数	纵波波速/ （m∙s−1）	干密度/ （g·cm−3）	饱和密度/ （g·cm−3）	饱和含 水率/%	孔隙 度/%
平均值	2 781	2.21	2.67	7.05	15.59

下载: 导出CSV 

| 显示表格

表  2  岩样分组

Table  2.  Grouping of the rock samples

组别	冻融0次	冻融30次	冻融60次
Ⅰ	A-1	B-1	C-1
Ⅱ	A-R	B-R	C-R

下载: 导出CSV 

| 显示表格

1.2   试验方案

试验主要包括３个方面：冻融循环试验、核磁共振检测、单轴分级加卸载蠕变试验。

（1）冻融循环试验

首先在压力值为−0.1 MPa的真空饱和仪中放入岩样，使岩样饱和，持续时间24 h，然后将饱和红砂岩岩样放入冻融循环试验箱中，进行冻融循环试验。冻融循环温度及高低温持续时间按照《工程岩体试验方法标准》（GBT 50266—2013）^[16]设定，冻结温度设定为−20 ℃，融解温度设定为20 ℃，冻结和融解时间均为12 h，每24 h为一个冻融循环周期，按表2分组进行冻融循环试验。

（2）核磁共振检测

采用磁感应强度为（0.3±0.05）T的NMC12–010V核磁共振测试分析系统进行检测。岩样在进行核磁共振检测前，通过真空抽气的方式进行强制饱和处理，完毕后，再从水中取出，将表面水分擦干，包上保鲜膜后检测。

（3）单轴分级加卸载蠕变试验

采用TAW-1000型岩石力学试验机，首先完成组Ⅰ中岩样的单轴抗压强度试验，获得不同冻融岩样的单轴抗压强度σ_ucs，然后根据组Ⅰ岩样的试验结果确定组Ⅱ岩样的各级加载应力值。用加载应力水平来表示预设加载应力与抗压强度的比值，不同冻融次数岩样的各级应力水平加载值如表3所示。加卸载方式采用应力控制，加载速率为0.01 MPa/s，每级荷载持续时间12 h后完全卸载至0，卸载速率为0.01 MPa/s，待变形趋于稳定后（约5 h）卸除岩样并进行核磁共振检测，得到该级加载比条件下岩样的核磁共振弛豫T₂谱和孔隙度等数据，然后进行下一级加载，依次循环，加载方式如图2所示。

表  3  单轴蠕变试验中不同冻融岩样各级应力加载值

Table  3.  Stress loading values at various levels in the uniaxial creep test of sandstone under different freezing-thawing conditions

冻融循环 次数/次	岩样 编号	冻融后T2 谱面积	加载应力/MPa
			第一级	第二级	第三级	第四级	第五级
0	A-R	10 668.82	10.01	13.34	16.68	20.01	23.35
30	B-R	11 584.31	9.71	12.95	16.19	19.42	22.66
60	C-R	11 348.65	8.40	11.20	14.00	16.79	19.59

下载: 导出CSV 

| 显示表格

图  2  冻融岩石分级加卸载蠕变应力加载方式

Figure  2.  Graded loading and unloading creep stress loading method of freezing-thawing rock

下载: 全尺寸图片 幻灯片

2.   试验结果

2.1   冻融岩石蠕变力学性质

2.1.1   单轴强度试验

图3为单轴抗压强度试验的应力-应变曲线，表4为试验结果。曲线中，OA是压实阶段，AB是线弹性阶段，BC是塑性发展阶段^[17]。据表4可知，抗压强度和弹性模量随冻融次数的增加而减小：冻融次数从0次增加到30次时，降幅分别为2.99%和4.69%；冻融次数从30次增加到60次时，降幅分别为13.53%和17.38%，较前者变化明显。这表明，岩样的力学性能受冻融循环作用的影响较大，其主要原因是岩样孔隙结构在受冻胀应力的反复作用时，内部孔隙的连续扩张和微裂纹的不断扩展。

图  3  应力-应变曲线

Figure  3.  Stress-strain curves

下载: 全尺寸图片 幻灯片

表  4  组Ⅰ中不同冻融岩样的力学参数

Table  4.  Mechanical parameters of different freezing-thawing rock samples in group I

参数	冻融次数/次
	0	30	60
抗压强度/MPa	33.37	32.37	27.99
弹性模量/GPa	39.65	37.79	31.22

下载: 导出CSV 

| 显示表格

2.1.2   蠕变行为

图4为不同冻融岩样的单轴分级加卸载蠕变试验曲线，根据图中数据可得到除最后一级外不同加载条件下的稳态蠕变速率和塑性应变，如图5所示。当应力水平在0.3σ_ucs～0.5σ_ucs时，二者均随应力水平的提升近似线性增长，当应力水平增高至0.6σ_ucs时，二者均急剧增大。此外，冻融作用使得岩样内部颗粒间的黏结作用减弱^[18]，故冻融次数越多，这一现象更为显著。例如，当应力水平从0.5σ_ucs提升至0.6σ_ucs时，对于冻融0次的岩样，塑性应变量从0.9421×10⁻³增加到0.1149×10⁻²，增幅为22.01%；然而，冻融30次岩样的塑性应变量从0.3419×10⁻²增加到0.5974×10⁻²，增幅为74.75%，冻融60次岩样的塑性应变量从0.6307×10⁻²急剧增加到1.2043×10⁻²，增幅为90.95%。

图  4  不同冻融岩样的蠕变曲线

Figure  4.  Creep curves of different freezing-thawing rock samples

下载: 全尺寸图片 幻灯片

图  5  蠕变速率、塑性应变分别与加载比的关系

Figure  5.  Relationship between the creep rate, plastic strain and loading ratio

下载: 全尺寸图片 幻灯片

2.2   蠕变过程中冻融岩石孔隙结构细观参数演化

2.2.1   孔径分布变化特征

孔径分布可由T₂曲线的核磁共振信号强度来表征^[19]。具体来说，孔隙数量随着信号强度的增加而增加^[20]，孔隙越大其弛豫时间T₂越长，且其表面积（S）与体积（V）之比越小，反之亦然^[21]。根据核磁共振理论^[22]，本文将T₂<10 ms的孔隙视为小孔，T₂>10 ms的孔隙视为大孔。

冻融0次，30次和60次后的岩样在不同应力水平的T₂分布曲线如图6所示，孔径分布占比如图7所示。可见，T₂谱面积随应力水平的提高而增长，且冻融次数越多趋势越明显，这是冻融循环削弱了内部颗粒间的黏聚力和摩擦力导致的。试验后岩样大孔（T₂>10 ms）占比约为小孔（T₂<10 ms）的1.5倍。当应力水平在0.3σ_ucs～0.5σ_ucs之间变化时，应力水平对孔径分布的变化没有明显影响，孔径分布在一个较小的范围内波动。由于岩样内部微小孔隙的萌生、扩展和连通是导致其孔隙结构变化的重要因素^[23]，所以孔径分布的变化是随机的。然而，当应力水平在0.5σ_ucs～0.6σ_ucs之间变化时，小孔占比均突然下降，大孔占比均突然上升。其原因在于高应力水平下，岩样内部损伤愈发严重，原来的小孔演化为大孔，导致大孔占比增加，小孔占比减少。

图  6  冻融红砂岩蠕变过程中T₂分布曲线

Figure  6.  T₂ distribution during the creep process of rock samples with different freezing-thawing times

下载: 全尺寸图片 幻灯片

图  7  蠕变过程中红砂岩的孔径分布占比

Figure  7.  Proportion of pore size distribution of rock samples during creep

下载: 全尺寸图片 幻灯片

2.2.2   孔隙度的演化规律

图8为随应力水平提升，不同冻融岩石孔隙度增量的变化情况。为量化孔隙度增量与应力水平间的关系，对图8中数据进行拟合，发现孔隙度增量与应力水平间存在良好的指数关系，指数方程可以表示为$y={y}_{0}+A{\mathrm{e}}^{-x/t}$。

图  8  孔隙度增量受应力水平的影响

Figure  8.  Effect of loading ratio on porosity increment

下载: 全尺寸图片 幻灯片

对于冻融后的岩样，当应力水平在0.3σ_ucs～0.5σ_ucs时，岩样孔隙度增量逐渐增大，但增幅较小，岩样内部仍以弹性变形为主，在力的作用下，内部颗粒之间胶结较弱的区域因失衡而破坏，产生新的裂隙^[24]。孔隙度增量在应力水平增长到0.6σ_ucs时迅速增加，冻融0次、30次和60次的岩样分别增加10.65%、17.98%和23.75%，岩样进入以塑形变形为主的阶段，胶结较弱的区域基本破坏，此时在外力作用下，原有的孔隙会继续扩展或多个孔隙发生贯通，表现为裂隙的扩展^[25]。以上现象表明，随着应力水平提高，岩样内部孔隙结构不断劣化，在应力水平达到并超过0.6σ_ucs时，蠕变阶段岩石内部损伤愈发严重。

图8中应力水平为0.6σ_ucs时，冻融30次和60次岩样孔隙度增量的增幅明显高于冻融0次岩样，为进一步探讨冻融循环次数与孔隙度增量的关系，对蠕变阶段不同冻融岩样孔隙度增量的变化情况进行分析，如图9所示。可以看出，随着冻融循环次数增多，孔隙度增量近似线性增长，说明在扩大孔隙和裂隙方面，冻融循环发挥了一定的推动作用；此外，高应力水平下冻融作用对孔隙度增量的影响更为明显，这是由于冻融过程中水冰相变不仅增加了孔隙体积，同时也使得岩石基质之间的胶结性质降低。

图  9  冻融次数对孔隙度增量的影响

Figure  9.  Effect of freezing-thawing times on porosity increment

下载: 全尺寸图片 幻灯片

2.2.3   蠕变过程中的分形特征

根据Frenkel-Halsey-Hill分形理论^{[26 − 27]}可得：

式中：r——孔隙半径/μm；

r_max——孔隙半径最大值/μm；

r_min——孔隙半径最小值/μm；

S_V——半径小于r的孔的累积体积分数；

D——分形维数。

由于r_min<<r_max，故式（1）可由式（2）表示为：

在核磁共振测试中，横向弛豫时间T₂可用来反映孔隙大小，因此可以将式（2）转换成式（3）：

式中：T_2max——最大横向弛豫时间/ms。

对式（3）两边取对数，可得：

本文利用式（4）来计算分形维数。

在孔径分类的基础上，利用式（4）在相应范围内线性拟合计算出小孔和大孔的分形维数D₁和D₂。此外，通过对整个孔径范围的线性拟合，得到总孔隙的分形维数D_T。由于篇幅限制，这里仅给出冻融0次岩样第一级应力加载后分形维数的线性拟合计算过程，如图10所示。

图  10  分形维数计算过程

Figure  10.  Fractal dimension calculation process

下载: 全尺寸图片 幻灯片

图11是随应力水平变化的不同冻融岩样分形维数的变化趋势。随着应力水平的提升，小孔的分形维数D₁在2附近上下波动。这表明小孔在蠕变过程中无分形特征。

图  11  冻融红砂岩分形维数随应力水平的演化

Figure  11.  Evolution of fractal dimension of freezing-thawing red sandstone with stress level

下载: 全尺寸图片 幻灯片

大孔的分形维数D₂与应力水平呈正相关，其数值范围在2～3之间，D₂在相同应力水平下，冻融次数越多其值越大。此外，当冻融次数从0次增加到60次时，D₂与应力水平的线性拟合曲线斜率随着冻融次数的增多，从0.056增加到0.194，增幅为255.6%。这可能与冻融和长期荷载对岩样孔隙结构的共同作用有关，受冻融作用越强的岩样内部骨架越松散，在长期荷载的作用下其大孔孔隙结构复杂程度变化更为剧烈。

总孔隙分形维数D_T随应力水平的变化趋势与孔隙度的变化趋势一致，呈现出先缓慢上升、应力水平到达0.6σ_ucs后急剧上升的趋势。在应力水平较低时，岩样抵抗变形的能力较强，长期荷载并没使得孔隙结构的复杂程度出现较大增长，但随着应力水平的提升，损伤不断积累，当作用于岩样的长期荷载接近其屈服强度时，大量裂缝产生并合并，岩样抵抗变形的能力变弱，孔隙体积和孔隙网络复杂性大大增加，导致分形维数D_T急剧增大。同时，冰的反复冻胀与融缩使得岩样内部颗粒间的黏聚力和胶结强度逐步降低，微小孔隙不断萌生和融合，进而导致内部缺陷和损伤增加，故受冻融损伤越严重的岩样在应力水平分界点（0.6σ_ucs）处D_T变化更为显著。

对图12中的数据进行线性拟合，发现冻融岩石的孔隙度与D_T呈正相关，且拟合曲线的相关系数（R²）均大于0.87，表明孔隙度与D_T之间高度相关。因此，随着分形维数的增加，冻融岩石的孔隙度增大，即分形维数越大，其孔隙网络的发育程度就越大。

图  12  D_T与孔隙度的关系

Figure  12.  Relationship between D_T and porosity

下载: 全尺寸图片 幻灯片

3.   冻融岩石细观结构演化对宏观特性的影响

为探究冻融岩石的孔隙度及孔隙结构复杂程度双因素共同作用对其宏观蠕变力学行为的影响，绘制塑性应变量及稳态蠕变速率分别与孔隙度和D_T的三维关系图（图13）。由于篇幅有限，这里仅给出冻融30次岩样的三维关系图。由图13（a）可知，岩样的塑性应变量随孔隙度的增大而增大，当孔隙度较小时，分形维数对塑性应变量的影响较小，当孔隙度较大时，随着分形维数的增加，塑性应变量增加。因此，冻融岩石孔隙结构的复杂程度只在孔隙度较大时才对塑性应变量产生明显影响，在孔隙度相同的情况下，分形维数越大，冻融岩石的塑性应变量越大。由图13（b）可知，冻融岩石稳态蠕变速率受孔隙度及分形维数的影响与塑性应变量相似。因此，对于受长期荷载作用的冻融岩石，在其蠕变阶段，随着时间的推移，次生裂纹逐渐萌生，孔隙度随之增长，孔隙结构逐渐复杂，进而影响了冻融岩石的宏观蠕变力学行为。

图  13  双因素共同作用对蠕变力学行为的影响

Figure  13.  Influence of two factors on creep mechanical behavior

下载: 全尺寸图片 幻灯片

4.   冻融岩石蠕变模型

4.1   冻融岩石损伤黏性元件

当施加的荷载达到或超越一定值时，岩石内部会有损伤产生，应考虑荷载产生的损伤对蠕变参数造成的影响^[28]，本文仅考虑加速蠕变阶段应力产生的损伤，因为减速蠕变阶段和稳定蠕变阶段岩石内部损伤程度较小^[29]。引入损伤变量（D_c）来描述加速蠕变阶段，受冻融作用影响岩样的黏性系数损伤劣化，图14为冻融条件下的损伤黏性元件。

图  14  冻融损伤黏性元件

Figure  14.  Freezing-thawing damage viscous element

下载: 全尺寸图片 幻灯片

Abel黏壶^[30]的本构关系为：

式中：σ——加载应力/MPa；

γ——阶数；

t——时间/h。

当γ=0时，上式表示理想弹性体；当γ=1时，上式表示牛顿体；当0≤γ≤1时，上式可表达介于理想弹性体和牛顿体之间的材料。

变系数Abel黏壶的本构关系:

式中：D_c——黏性系数的劣化。

冻融岩样在蠕变过程中的损伤与时间呈负指数函数关系^[29],即：

式中：α——与红砂岩性质及冻融次数相关的系数/h^–1。

变系数Abel黏壶的本构关系，可结合式（5）—（7）得到：

式中：${\eta }_{2}^{\gamma }$e^–αt——与冻融次数相关且随时间演化的黏性 系数/（GPa·h）。

当应力保持不变时，有：

式中：K——累加次数；

$\Gamma$——伽马函数。

4.2   蠕变模型的建立

由上文可知，不同冻融岩样依次经历减速、稳态和加速蠕变阶段，且冻融次数、应力状态和加载时间均影响岩样的蠕变特性。经典的Poyting-Thomson模型能较好地反映岩样蠕变的减速蠕变和稳态蠕变阶段，却无法反映加速蠕变阶段。因此，本文在Poyting-Thomson模型基础上串联一个由变系数Abel黏壶和塑性元件并联组成的黏塑性体表示加速蠕变阶段。

图15为改进的Poyting-Thomson模型，其中E₁、E₂和${\eta }_{1}$是与冻融次数相关的参数，σ_s为长期强度。其本构方程为：

图  15  改进的Poyting-Thomson模型

Figure  15.  Improved saturated Poyting-Thomson model

下载: 全尺寸图片 幻灯片

当σ<σ_s时：

当σ≥σ_s时：

4.3   模型验证及参数辨识

依据上文蠕变试验数据，采用1stOpt软件内置的通用全局优化算法识别模型参数，得到各应力水平下的模型参数及拟合结果。由于篇幅限制，这里仅列出了冻融30次岩样的模型参数及拟合结果，分别如表5和图16所示。不同冻融岩样的模型拟合结果相似，均与试验结果吻合良好，能够同时描述瞬时变形、减速蠕变、稳态蠕变和加速蠕变过程。

表  5  蠕变模型参数

Table  5.  Creep model parameters

加载等级	${E_1}$ /GPa	${E_2}$ /GPa	${\eta _1}$ /（GPa·h）	${\eta _2}$ /（GPa·h）	$\gamma$	$\alpha$/h−1	${R^2}$
第一级	129.969	19.024	2.245	—	—	—	0.97
第二级	89.972	13.931	1.568	—	—	—	0.95
第三级	52.979	12.401	1.496	—	—	—	0.98
第四级	22.135	9.671	1.054	—	—	—	0.93
第五级	10.973	8.066	0.809	0.001 3	0.161	0.169	0.91
注：“—”表示此处为空。

下载: 导出CSV 

| 显示表格

图  16  蠕变模型验证曲线

Figure  16.  Creep model validation curve

下载: 全尺寸图片 幻灯片

5.   结论

（1）当应力水平在0.3σ_ucs～0.5σ_ucs间变化时，应力水平对孔径分布的变化没有明显影响，所有岩样的孔径分布在一个较小的范围内波动，孔隙度增幅较小；当应力水平在0.5σ_ucs～0.6σ_ucs间变化时，大孔占比均突然增加，小孔占比均突然降低，孔隙度增幅迅速增大，且冻融加剧了这一现象。

（2）小孔没有明显的分形特征，D₂与应力水平之间呈正相关，D_T随应力水平的提升先缓慢增大，然后在0.6σ_ucs处突增，且D_T与孔隙度间的关系可用正线性函数描述。

（3）孔隙结构复杂程度仅在孔隙度较大时才对冻融岩石蠕变特性产生明显影响，分形维数越大，蠕应变量和蠕变速率在相同孔隙度下越高；孔隙结构复杂程度和孔隙度越大，岩石受损越严重。

（4）在Poyting-Thomson模型的基础上，根据冻融岩石的蠕变特性，串联一个改进的黏塑性体，建立考虑了冻融影响的蠕变模型。经验证，蠕变模型的拟合结果与试验结果吻合良好，能准确反映不同冻融岩样各加载阶段的蠕变状态。