董方方，朱谭谭，屈子健

(长安大学 公路学院，陕西 西安 710064)

寒冷地区岩体冻融常引发冻融滑坡、隧道冻害等工程灾害[1]。冻融循环作用下岩石的力学行为一直是岩土工程领域的研究热点。国内外研究人员采用不同岩石对冻融循环作用下岩石力学行为开展了一系列研究工作。杨鸿锐等[2]研究了不同温度循环范围内麦积山石窟饱水砂砾岩波速损失率、质量损失率和抗拉强度的变化规律及特征。王鲁男等[3]研究了化学腐蚀与冻胀荷载耦合作用下粉砂岩强度长期衰减规律和预测模型。杨秀荣[4]研究发现岩样在长期荷载作用下的破坏模式均为剪切破坏，且岩石破坏程度随着冻融次数的增加逐渐严重。张峰瑞等[5]提出了冻融岩石损伤黏性元件，构建了花岗岩冻融剪切蠕变本构模型。崔凯等[6]通过研究认为，冻融、溶蚀和盐分结晶等条件不同会导致孔隙的类型、数量和联通情况产生差异，进而引起岩石宏观性质出现不同程度的累积性损伤。张慧梅等[7]研究发现岩石损伤受缺陷随机分布和应力状态的影响，抗拉强度对岩石缺陷更为敏感。闻磊等[8]通过研究得到了抗拉强度与冻融循环次数的关系方程。然而，上述多数研究基于室内试验开展，数值模拟在岩石冻融力学行为研究中应用较少，且大多采用有限元法或有限差分法[9-12]。有限元法和有限差分法无法模拟孔隙水相变膨胀物理过程，且在岩石裂隙扩展等方面也具有一定局限性。基于此，本文提出采用孔隙水颗粒和颗粒膨胀的概念模拟低温条件下岩石孔隙水低温相变膨胀过程，并通过与室内试验结果对比验证了该方法的可行性和可靠性。采用该方法开展了岩石冻融循环后巴西劈裂试验，研究了岩石冻融循环后宏、细观拉伸力学行为。

1 基于颗粒流的冻融循环数值模拟方法

1.1 岩石冻融损伤机理

如图1所示，作为一种多孔介质，岩石由矿物颗粒和孔隙组成。在富水条件下，环境中的水会进入岩石填充岩石内部孔隙。温度大于0 ℃时，岩石内的水主要以液体形式存在。当温度低于0 ℃时，岩石内的孔隙水逐渐冻结成冰。由于冰的密度小于水的密度，孔隙水冻结成冰后体积增大挤压周围矿物颗粒，岩石内部产生冻胀力(图1)。冻胀力大于矿物颗粒间粘结强度或矿物颗粒自身强度时，岩石内部产生微裂隙，岩石孔隙增大，岩石体积膨胀。如果岩石处于富水环境，环境水会继续进入岩石内部充填新增孔隙。在长期冻融循环作用下，岩石内部损伤累积，岩石力学性能衰退。Winkler[13]通过分析岩石内部水分相变膨胀规律发现，孔隙冰在-5 ℃、-10 ℃和-22 ℃时膨胀压力可分别达到61、113和211.5 MPa，并认为冻融条件下岩石的抗拉强度比抗压强度更值得关注。

图1 岩石细观结构及冻胀原理示意图Fig.1 Schematic diagram of rock meso-structure and frost heave principle

1.2 基于颗粒流的冻融数值方法

富水环境下岩石冻融损伤的主要原因是孔隙水冻结成冰后体积膨胀，挤压矿物颗粒导致岩石产生宏、细观裂隙。基于此，采用以下方法可以模拟岩石冻胀物理过程：

(1)将低温条件下的富水岩石简化为岩石矿物颗粒和孔隙水；

(2)将岩石矿物颗粒和孔隙水分别简化成球形颗粒，即矿物颗粒和孔隙水颗粒；

(3)忽略矿物颗粒热胀冷缩引起的体积变化；

(4)孔隙水颗粒的体积随温度的变化而变化。

根据上述方法，采用PFC3D颗粒流软件，建立低温条件下的富水岩石数值模型如图2所示。其中，矿物颗粒用于模拟岩石的矿物，孔隙水颗粒用于模拟低温条件下岩石内的孔隙水。岩-岩粘结是岩石矿物颗粒之间的粘结，岩-冰粘结是岩石矿物颗粒和孔隙水颗粒之间的粘结，冰-冰粘结是孔隙水颗粒和孔隙水颗粒之间的粘结。在数值计算模型中，给三种不同颗粒和不同粘结赋予不同的力学参数。

图2 低温富水岩石数值模型细观结构Fig.2 Meso-structure of numerical model for water rich rock under low temperature

当温度低于0 ℃时，岩石中的孔隙水不会全部冻结成冰，部分孔隙水仍然以液态形式存在，这部分水即为未冻水。未冻水占总孔隙水的比例为未冻水含量。刘泉声等[14]通过分析累计孔隙体积的分布规律，得到岩石中未冻水含量与温度的关系表达式为

(1)

T=f(t)

(2)

式中，wu为未冻水含量，T为岩石的温度，t为冻融循环所持续的时间。岩石中的孔隙水冻结成冰后体积增大。数值计算过程中孔隙水颗粒体积随温度的变化可用下式表示：

(3)

其中，V为任意温度下孔隙水颗粒的体积，V0为0 ℃时孔隙水颗粒的体积(此时孔隙水颗粒为液态形式，不产生冻胀力)，ρw和ρi分别为水和冰的密度。岩石在冻融循环过程中会发生变形，其中部分变形是不可恢复变形。假定每次冻融循环后不可恢复变形占总变形的比例为η，则经过i次冻融循环作用后任意孔隙水颗粒的体积可以表示为

(4)

(1)建立低温富水岩石数值模型，颗粒大小在一定范围内随机分布，颗粒位置随机分布；

(2)给数值模型赋参数，并记录每个孔隙水颗粒的初始体积；

(3)获取当前冻融循环持续时间，采用式(2)计算当前时间点的温度；

(4)采用式(1)计算当前温度下未冻水含量；

(5)采用式(3)和式(4)计算每个孔隙水颗粒的当前条件下的体积目标值；

(6)改变孔隙水颗粒体积至目标值，并循环程序达到力学平衡；

(7)重复上述步骤，直到达到目标冻融循环次数。

1.3 冻融循环数值方案

在采用PFC3D开展岩石力学数值模拟时，常用的粘结模型有线性平行粘结模型和线性接触粘结模型。其中，线性平行粘结模型能够在颗粒之间传递力和力矩，且强度服从摩尔库伦准则。而线性接触粘结模型只能在颗粒之间传递力，不能传递力矩。所以本文在数值模拟中采用线性平行粘结数值模型。在冻融循环数值计算中，粘结和颗粒所采用的细观力学参数分别如表1和表2所示。需要注意的是，巴西劈裂数值计算时为常温环境，此时试样内只有矿物颗粒及其之间的粘结，孔隙水颗粒已经被删除。所以孔隙水颗粒之间和孔隙水颗粒与矿物颗粒之间粘结的破坏对岩石试样自身的力学行为的影响可以忽略不计。因此，为了消除这些粘结的破坏对数值模拟结果的影响，本文将这些粘结的抗拉强度和粘聚力设置为足够大(100 MPa)。冻融循环数值模拟中，忽略室温条件下孔隙水对岩石力学行为的影响，只考虑温度小于0 ℃的情况。由式(1)和式(3)可知，当温度低于-25 ℃时，岩石内未冻水含量接近于0，孔隙水颗粒体积基本不再发生变化。因此，将冻融循环温度变化范围设为-25 ℃～0 ℃。基于自然界中气温随时间的变化规律[15-16]，温度与冻融循环持续时间的关系采用正弦函数表示，循环周期为60 min，即温度由0 ℃降低到-25 ℃的时间需要30 min。温度随时间的变化可用下式表示：

表1 冻融循环数值模拟中粘结的细观力学参数

表2 冻融循环数值模拟中颗粒细观力学参数

(5)

采用上述方法建立巴西圆盘数值模型如图3所示，圆盘数值模型为圆柱形，直径和高度均为50 mm。数值模型共包含74 423个矿物颗粒，33 527个孔隙水颗粒，347 483个岩-岩粘结，157 700个岩-冰粘结和9 174个冰-冰粘结。数值模拟过程中首先开展冻融循环试验，冻融循环5、10、15、20、25和30次后分别开展巴西劈裂试验，巴西劈裂荷载施加采用位移控制模式，位移速率为0.15 m/s。根据Zhang等[17]和Potyondy等[18]的研究成果，这一加载速率在岩石力学颗粒流数值模拟中能够满足静态力学试验的要求。

图3 巴西圆盘数值模型Fig.3 Numerical model of Brazilian disc

2 结果与分析

2.1 变形和强度特征

不同冻融循环次数后巴西劈裂荷载-位移曲线见图4。由图可知，冻融后巴西劈裂荷载-位移曲线整体上分为线性、非线性和峰后三个阶段。线性段荷载随位移的增大基本呈线性变化，且随着冻融循环次数的增大，荷载-位移曲线线性段的斜率逐渐减小。这是因为随着冻融循环次数的增大，试样在冻融过程中产生的裂隙增多，试样抵抗变形的能力减弱。荷载-位移曲线线性段斜率随冻融循环次数的增大而减小这一变化规律，在刘慧等[19]开展的冻融后巴西劈裂室内试验中也得到了验证(图4(b))。在非线性段，竖向荷载随位移的增大呈非线性变化，此时试样在竖向荷载作用下产生了较多的微裂隙，试样抵抗变形的能力进一步降低，试样进入塑性变形阶段。随着冻融循环次数的增大，试样的峰后特征逐渐明显，这与图4(b)所示室内试验结果相一致。这是因为冻融循环次数增大导致试样内微裂隙增多，试样逐渐由脆性破坏向塑性破坏转变。值得注意的是，巴西劈裂室内试验中，荷载-位移曲线在初始加载阶段常表现出一个下凹型非线性阶段，例如尹乾等[20]和Zhu等[21]开展的室内巴西劈裂试验荷载-位移曲线均有一个下凹非线性段。然而，如图4(a)所示，数值模拟中荷载-位移曲线没有出现下凹非线性段。这是因为数值模拟中无法模拟试验仪器各构件之间间隙压缩和试样内微裂隙的张开和闭合。

图4 冻融循环后巴西劈裂荷载-位移曲线Fig.4 Load-displacement curve of rock samples after freeze-thaw cycles in Brazil split test

冻融循环后试样间接抗拉强度随冻融循环次数的变化规律见图5，其中室内试验抗拉强度均为平均值。由图可知，数值模拟和室内试验岩石间接抗拉强度随冻融循环次数的增大均表现为逐渐减小。这是因为孔隙水低温相变膨胀导致试样产生裂隙损伤，随着冻融循环次数的增大，试样内冻胀损伤累积(图6)，试样承载能力降低。如图5所示，数值模拟中试样抗拉强度随冻融循环次数的增大基本呈线性减小。冻融循环次数由0增大到30次，试样抗拉强度由6.36 MPa减小到3.30 MPa，减小了48.11%。根据闻磊等[8]和Zhang等[22]的室内试验结果，红砂岩、灰岩、石英砂岩和花岗岩的抗拉强度均随冻融循环次数的增大逐渐减小。其中，红砂岩的抗拉强度随冻融循环次数的增大基本呈线性减小，这与数值模拟结果一致。然而，花岗岩、灰岩和石英砂岩的抗拉强度随冻融循环次数的增大呈非线性减小。所以，岩石抗拉强度随冻融循环次数的变化规律与岩石自身的性质有关。

图5 数值模拟和室内试验间接抗拉强度随冻融循环次数的变化规律Fig.5 Simulated and experimental indirect tensile strength of rock sample after freeze-thaw cycles

2.2 损伤演化过程

冻融循环过程中试样损伤演化特征如图6(a)—(f)所示，Khanlari等[23]采用砂岩开展室内冻融试验得到的试样破裂形态如图6(g)所示，冻融循环数值模拟过程中试样内微裂隙数量变化如图6(h)所示。如图6(a)所示，当冻融循环次数n达到5次时，试样内部已经开始产生少量的裂隙。随着冻融循环次数的增大，试样内微裂隙数量逐渐增多。随着裂隙数量的增多，裂隙在试样表面局部区域产生聚集，见图6(b)—(f)，椭圆形标记范围内尤为明显。这主要是因为冻融循环过程中试样内部受到外部的约束作用，能够承受更大的冻胀力。而试样表面不受约束作用，冻融循环过程中容易发生破坏。这与室内试验(图6(g))结果相似，证明本文所采用的数值方法能够较为真实地反映岩石在冻融循环过程中的破裂形态。如图6(h)所示，随着冻融循环次数增大，裂隙数量的增长速率也逐渐增大，试样内裂隙数量随冻融循环次数增大基本呈指数形式变化。

图6 冻融循环过程中试样损伤演化特征及其与室内试验结果对比Fig.6 Damage evolution characteristics of numerical model during freeze-thaw cycles and the comparation between the results from numerical simulation and laboratory

不同冻融循环次数后的圆盘在巴西劈裂实验中的破坏特征如图7所示。需要注意的是，图7中的所有裂隙均为巴西劈裂试验中产生的裂隙，冻融循环过程中产生的裂隙未显示。由图可知，当冻融循环次数较小时，圆盘主要表现为中心劈裂破坏，中心裂隙的宽度较小，裂隙主要集中在圆盘竖向对称轴附近。随着冻融循环次数的增大，圆盘的中心裂隙宽度明显增大，在远离圆盘竖向对称轴的位置出现了大量的裂隙(图7(g))。随着冻融循环次数的增大，劈裂试验开始前圆盘内已经存在大量裂隙。劈裂试验中，圆盘在竖向荷载作用下已存在裂隙容易产生起裂扩展。所以，随着冻融循环次数的增大，圆盘偏离中心线位置的裂隙逐渐增多。此外，圆盘中心裂隙上下端部呈分叉形态。这一现象与Zhu等[21]室内巴西劈裂试验中圆盘破裂特征相同(图8)。在竖向荷载作用下，圆盘上下与墙体(数值计算)或压头(室内试验)为小面积面接触。圆盘侧向变形受到墙体或压头的约束作用，导致圆盘在上下两端很难沿水平方向开裂。这导致中心裂隙扩展到上下端部时向两个方向分别扩展，中心裂隙在上下端部表现为分叉形态(图7、图8)。

图7 不同冻融循环次数后巴西圆盘破裂特征Fig.7 Failure characteristics of Brazilian disc after different freeze-thaw cycles

图8 巴西劈裂室内试验圆盘破坏形态[21]Fig.8 Failure patterns of Brazilian discs in laboratory tests[21]

定义巴西劈裂试验中试样中产生的裂隙面法线方向与加载方向的夹角为裂隙倾角，裂隙倾角变化范围为0到90°。将裂隙根据倾角大小每隔5°分为一组，如倾角为0～5°的裂隙为一组。不同倾角范围内裂隙数量占总裂隙的变化如图9所示。其中图9只统计巴西劈裂试验中产生的裂隙，剔除了冻融循环过程中产生的裂隙。巴西劈裂试验中，试样在水平方向承受拉应力，试样产生拉伸破坏，裂隙容易沿加载方向起裂扩展。所以，如图9所示，在任意冻融循环次数下，裂隙所占比例随倾角的增大逐渐增大。以冻融循环次数为10的试样为例，倾角在0～90°范围内的裂隙(每间隔5°为一组)的比例分别为0.22%、0.58%、1.05%、1.45%、2.03%、2.74%、3.35%、4.05%、4.63%、5.68%、6.42%、7.23%、8.15%、8.95%、9.75%、10.66%、11.29%、11.78%。随着冻融循环次数的增大，倾角较小的裂隙所占比例逐渐增大，倾角较大的裂隙所占比例逐渐减小。例如，冻融循环次数由0增加到30次，倾角在85°～90°范围内的裂隙比例由11.81%减小到10.02%，减小了15.16%。

图9 不同倾角裂隙比例变化规律Fig.9 Change in the proportion of cracks with different dip angle

根据相关研究发现，拉伸破坏率或剪切破坏率能够反映岩石试样的破坏形式及破坏诱因。定义冻融循环过程和巴西劈裂试验中拉伸裂隙占总裂隙的比例为拉伸破坏率，其计算见下式：

(6)

其中，Rt为拉伸破坏率，Nt和N分别为试样内拉伸裂隙数量和总裂隙数量。在冻融循环过程和巴西劈裂试验中，试样拉伸破坏率随冻融循环次数的变化规律如图10所示。其中，在计算巴西劈裂试验中拉伸破坏率时，只统计了劈裂试验过程中产生裂隙，剔除了劈裂试验前即冻融循环过程中产生的裂隙。在冻融循环过程中，孔隙水相变膨胀，部分矿物颗粒之间承受拉应力。而岩石的抗拉强度远小于其抗压强度。因此，在冻融循环过程中，试样内的裂隙主要为拉伸裂隙。冻融循环过程中试样的拉伸破坏率的最大值和最小值分别为98.15%(n=5次)和92.53%(n=30次)。这也证明了Winkler[13]认为冻融作用下岩石抗拉强度比抗压强度更值得关注这一观点的正确性。随着冻融循环次数的增大，试样内裂隙增多导致试样损伤累积，力学性能衰退。试样抵抗剪应力的能力减弱。所以，随着冻融循环次数的增大，试样的拉伸破坏率减小，即剪切破坏率增大。

在劈裂试验中，试样在水平方向承受拉应力，试样的破坏仍以拉伸破坏为主。劈裂试验中试样的拉伸破坏率最大值和最小值分别为96.88%(n=30次)和83.37%(n=5次)。如图10所示，在劈裂试验中，拉伸破坏率随着冻融循环次数的增大逐渐呈非线性增大。这是因为巴西劈裂过程中试样水平方向以拉应力为主，随着冻融循环次数的增大，试样力学性能衰退，抗拉性能减弱，拉伸破坏率增大。

图10 冻融循环和巴西劈裂中试样拉伸破坏率变化规律Fig.10 Tensile failure rate of the sample during the processes of freeze-thaw and Brazil splitting test

3 结论

本文提出了基于孔隙水颗粒和颗粒膨胀概念的颗粒流方法模拟孔隙水低温相变膨胀物理过程，并采用该方法开展了岩石冻融循环试验和冻融后巴西劈裂试验，研究了岩石冻融循环作用下的损伤演化过程和冻融后拉伸力学行为。可以得出以下主要结论：

1)采用本文所提出的方法开展冻融循环试验和冻融后巴西劈裂试验得到的结果与室内试验结果相一致，本文提出的方法在岩石冻融损伤和冻融后拉伸力学行为研究领域具有可行性和可靠性。

2)冻融循环过程中，裂隙更容易在试样表面起裂聚集，导致试样产生表面剥落破坏。冻融损伤累积会导致试样在巴西劈裂试验中产生大量远离试样中心线的裂隙。

3)巴西劈裂试验中，倾角较大的裂隙数量较多。冻融循环会导致劈裂试验中大倾角裂隙比例减小，小倾角裂隙比例增大。冻融循环和巴西劈裂过程中试样拉伸破坏率均远大于剪切破坏率。