0.   引言

在工程建设中,岩体经常会与水发生接触,如坝基岩体、库岸边坡、深埋引水隧道等,都会不同程度地受到地下水或地表水对其的软化作用,导致岩体的力学性质发生明显的降低,如峰值强度、黏聚力、弹性模量等^[1-2]。因此,研究水对岩石的软化作用具有重要的实际意义。

水对岩石软化的作用机制较为复杂,至今已有许多学者通过各种技术手段对水软化过程中的不同问题进行了研究。杨春和等^[3]研究了板岩遇水软化后的过程和机制：随着浸泡时间的延长,板岩的吸水率增加,板岩内部结构变得松散,颗粒间的黏结力降低,从而使得岩石的力学强度减小；He等^[4]研究了泥岩力学性质随含水率的变化规律；闫章程等^[5]等采用声发射技术,研究了干燥和饱水灰岩的声发射特征差异,证明水对声发射信号具有明显的弱化效应。朱宝龙等^[6]研究表明,水岩作用使得页岩黏土化趋势显著,由矿物间致密的片层状结构逐渐变为松散的絮状结构,孔隙度增加。上述研究均表明,岩石吸水后宏观力学性质劣化的重要原因之一是微细观结构的变化^[7]。

随着计算机技术的发展,数值模拟已成为研究岩石力学性质的一种重要手段。近年来,基于离散元法的颗粒流软件（particle flow code, PFC）已在各个岩土工程领域取得了成功的应用,成为从细观角度研究岩石力学问题的重要方法^[8]。如柳万里等^[9]运用颗粒流软件分析了岩石在天然和饱和两种状态下的能量演化规律；邓华锋等^[10]发现水岩作用导致岩石内部颗粒强度逐渐弱化,非均匀性增强,承载力降低；卞康等^[11]的结果表明随吸水时间增加,颗粒流模型的剪切带宽度变大,主剪切面的贯通程度提高,岩样的破碎程度明显提升；Zhao等^[12]通过对摩擦系数和黏结强度进行折减,模拟水对岩石力学性质的弱化作用。Potyondy^[13]在应力腐蚀理论基础上提出了平行黏结应力腐蚀模型,以模拟应力或水溶液等对岩石胶结的应力腐蚀作用。上述研究^[9-13]已经充分证明颗粒流软件在水岩作用领域的适用性和正确性。但是这些研究大部分是不同含水率下采用不同的细观参数,且参数的选取和室内试验联系不明确。

成岩过程中,岩石经受了漫长的地质历史时期的作用,其内部胶结物的强度存在很大的差异。因此胶结物处力学性质存在一定的不同,所以水对胶结的腐蚀程度也不同。在本文的数值建模当中,将胶结物强度的非均质性考虑进去,将其假定为服从Weibull分布。

本文通过构建基于平行黏结半径乘子的损伤因子,将页岩的室内试验^[14]和数值模拟联系起来,模拟不同泡水时间下的页岩力学性质（如单轴压缩强度、弹性模量等）弱化过程,提出了平行黏结水弱化模型（parallel- bonded water-weakening, PWW）。该页岩中含有35%左右的绿泥石,绿泥石在水的作用下膨胀、溶解,极大地影响着页岩的水理性质^[14]。由于在室内试验当中很难定量研究岩石胶结物的强度分布对岩石力学性质的影响,因此本文所提PWW模型耦合了Weibull分布,来表征胶结物的强度分布,以更加接近岩石存在非均质性的实际情况。通过分析单轴压缩试验中,不同黏结面积系数（不同泡水时间）下模型的胶结物非均质性、能量耗散、微裂纹的倾角和破坏模式等,揭示水对页岩的弱化作用的细观机制。

1.   平行黏结水弱化模型（PWW模型）

在平行黏结接触模型中,平行黏结半径乘子$\lambda$是一个极其重要的细观参数,对模型的宏观力学响应产生很大的影响^[15]。

平行黏结水弱化模型除颗粒流模型现有的假设^[15]之外,另有4个假设：

（1）将吸水弱化过程中岩石孔隙率增加等带来的不利影响均等效为平行黏结半径乘子的下降。

（2）水弱化效应仅影响颗粒间的胶结接触,对颗粒本身的几何形态不产生影响,因此每个平行黏结接触都是发生水弱化效应的潜在位置。

（3）水弱化只发生在接触表面,不论浸水时间多长,折减平行黏结半径乘子都瞬间完成。

（4）模型中不同接触的平行黏结半径乘子不同,代表岩石中胶结物的非均质性。

接下来详细介绍平行黏结水弱化模型的原理。

1.1   宏观损伤因子

当岩石吸水软化以后,力学性质劣化。随着泡水时间的增加,岩石内部的损伤逐渐累积,由损伤力学可知,岩石的宏观损伤因子可由室内试验测得的弹性模量得到^[14]

式中 $D_t$是岩石泡水t天以后的宏观损伤因子；E₀为岩石干燥状态下室内试验测得的弹性模量；$E_t$为岩石泡水t天以后室内试验测得的弹性模量。依据室内试验^[14]结果绘制宏观损伤因子随浸水时间的变化曲线如图1所示。如图所示,两者符合幂函数关系（如式（2）所示）,即随着浸水时间的增加,岩石含水率慢慢趋于稳定,弹性模量的降低也趋于稳定,不会无限制地降低。

图  1  宏观损伤因子与浸水时间的关系

Figure  1.  Relationship between macro-damage factor and water immersion time

下载: 全尺寸图片 幻灯片

1.2   胶结物力学参数的非均质性

许多学者在理论推导^{[14, 16]}或者数值模拟^[17-18]中均采用力学参数服从Weibull分布来表征岩石力学参数的非均质性,结果表明,非均质性对岩石的强度和裂纹扩展有很大影响。非均质性对岩石强度的最大影响在于会引起局部应力集中现象,从而加快微裂纹的萌生与扩展^[19]。

由于在室内试验当中很难定量研究岩石胶结物的强度对于岩石力学性质的影响,特别是吸水软化后的岩石力学性质,因此本文中采用数值模拟具有一定的优势。为了表征岩石胶结物的非均质性,本文将颗粒流模型中的细观参数平行黏结半径乘子$\lambda$的取值设定为服从Weibull分布。平行黏结半径乘子$\lambda$的概率密度函数为

累积分布函数为

式中,m为均质性因子（当m越小时,表示模型非均质性越大）,$\overline{\lambda }$为平行黏结半径乘子特征值。

颗粒流模型中第i个接触的平行黏结半径乘子${\lambda }_i$的取值为

式中,${\xi }_i$为第i个接触上服从(0,1]的均匀分布的随机数。

接下来描述二维情况下,如何计算当平行黏结半径乘子$\lambda$服从Weibull分布时两颗粒之间的力和力矩。

单个平行黏结的面积：

式中 ${\overline{R}}_i$为第i个接触的平行黏结半径；R_i为第i个接触两端的颗粒半径中的更小值,R_i=min(R⁽¹⁾,R⁽²⁾),R⁽¹⁾,R⁽²⁾分别为接触两端的颗粒半径；${\overline{A}}_i$为第i个接触的平行黏结面积。

当颗粒受到外部荷载时,颗粒流软件在计算时会进行接触的力和力矩的传递与更新。平行黏结半径乘子$\lambda$服从Weibull分布时两颗粒之间的接触力和力矩计算公式为

式中 ${\overline{F}}_i^{\text{n}},{\overline{F}}_i^{\text{s}},{\overline{M}}_i^{\text{n}}$分别为第i个接触的时步更新前的平行黏结法向接触力、切向接触力、平行黏结法向力矩；$F_i^{\text{n}},F_i^{\text{s}},M_i^{\text{n}}$分别为第i个接触的时步更新后的平行黏结法向接触力、切向接触力、平行黏结法向力矩；$\Delta {\delta }_i^{\text{n}},\Delta {\delta }_i^{\text{s}},\Delta {\theta }_i^{\text{n}}$平行黏结法向和切向位移增量、平行黏结的相对转角。

1.3   细观损伤因子

构建基于平行黏结半径乘子特征值的细观损伤因子：

式中 $D_{\text{p}}$为细观损伤因子；$\overline{\lambda }$为平行黏结完整状态下的半径乘子特征值（表示岩石干燥状态时的颗粒流模型,此时$\overline{\lambda }$=1）；${\overline{\lambda }}_{\text{w}}$为平行黏结接触折减以后的半径乘子特征值（表示岩石浸水t天以后的颗粒流模型）。

由式（10）可知,${\overline{\lambda }}_{\text{w}}$可由细观损伤因子$D_{\text{p}}$和岩石干燥状态时的颗粒流模型的平行黏结半径乘子特征值$\overline{\lambda }$计算得出

由于宏观损伤因子和颗粒流模型中定义的细观损伤因子属于不同维度、不同原理的两种损伤因子,无法确定具体的数学函数表达式,因此也不能直接从宏观损伤因子推导出细观损伤因子的取值。但是由于两者均是描述同一岩石受水软化作用后的损伤程度,所以宏观损伤因子与浸水时间的函数关系式对推测细观损伤因子与浸水时间的函数关系式有重要的指示作用。现假定两者是服从同一函数形式的,即细观损伤因子与浸水时间的关系应当如下式所示：

式中,A,B,C为材料参数,通过参数标定确定数值。

2.   离散元模型的建立

所建立的页岩颗粒流模型如图2所示,建立尺寸大小为50 mm×100 mm的矩形颗粒流模型。颗粒流模型的孔隙率为0.1,其颗粒半径服从0.15～0.249 mm之间的均匀分布,采用半径扩大法共生成35262个颗粒,平行黏结接触个数N=86163。通过赋予上下墙体一定的速度,达到准静态加载,模拟室内单轴压缩试验。

图  2  颗粒流模型示意图

Figure  2.  Schematic diagram of particle flow code model

下载: 全尺寸图片 幻灯片

当确定浸水t天以后颗粒流模型的平行黏结半径乘子特征值${\overline{\lambda }}_{\text{w}}$后,进而根据式（6）可以确定模型的总平行黏结面积A。在二维离散元模拟中,常用饱和度来表征岩石受水的影响^[20]。由于室内试验中很难对颗粒胶结物在吸水后的溶解进行测量,因此假定溶解的胶结物含量等于平行黏结模型中接触胶结面积的减小量。在微观上,由水所引起的力学性质弱化现象是由于颗粒之间的胶结物含量不断减少的过程。由于数值模拟中的饱和度^[20]和室内试验中的饱和度容易引起混淆,所以本文仿照文献[20]的定义,在本文PWW模型中,定义黏结面积系数$S_r$,其计算公式为

式中,A_max,A_w,A_min分别为干燥状态下模型的总平行黏结面积、浸水t天以后模型的总平行黏结面积和饱和状态下模型的总平行黏结面积,N为平行黏结接触总数。当黏结面积系数越大时,表征岩石含水率越高,其所受到的水软化作用越明显。

PWW模型参数标定步骤为：首先试算式（12）中的材料参数A,B,C,然后通过式（1）,（11）和（12）就可确定不同浸水时间下颗粒流模型的平行黏结半径乘子的特征值；最后再根据室内试验测得的不同浸水时间下岩石的单轴抗压强度标定颗粒流模型中的其它细观参数,如平行黏结抗拉强度,平行黏结弹性模量等。通过反复试算所有的细观参数,最终使得所建立的颗粒流模型和室内试验所获得的宏观力学性质一致。

针对不同岩石,式（12）中的材料参数不同,根据参数标定的结果取值。通过对页岩的室内试验^[14]进行参数标定,材料参数A,B,C为-0.597,26.847,0.599时较为合适。细观损伤因子随浸水时间的关系如图3所示,和图1中的宏观损伤因子随浸水时间的关系类似。

图  3  细观损伤因子随浸水时间的关系

Figure  3.  Relationship between mesoscopic damage factor and immersion time

下载: 全尺寸图片 幻灯片

表1列出根据PWW模型计算得出的不同浸水时间下颗粒流模型的平行黏结半径乘子特征值,其他的细观参数取值如表2所示。采用表1,2中的细观参数,共建立7个颗粒流模型。

表  1  PWW模型计算参数

Table  1.  Parameters of PWW model

时间/d	弹性模量/GPa	宏观损伤因子	细观损伤因子	平行黏结半径乘子特征值	黏结面积/m2	黏结面积系数
0	27.38	0.000	0.000	1.000	30.14	0.000
2	25.67	0.062	0.051	0.949	28.62	0.087
15	20.90	0.236	0.251	0.749	22.61	0.432
30	17.59	0.357	0.412	0.588	17.73	0.712
50	13.55	0.505	0.514	0.486	14.65	0.888
70	12.36	0.548	0.550	0.450	13.57	0.950
90	11.40	0.583	0.579	0.421	12.70	1.000

下载: 导出CSV 

| 显示表格

表  2  颗粒流模型细观参数

Table  2.  Mesoscopic parameters of particle flow code model

参数	取值
颗粒密度/(kg·m-3)	2650.0
最小粒径/mm	0.15
粒径比	1.66
颗粒接触模量/GPa	13.5
颗粒刚度比	3.0
颗粒摩擦系数	0.5
均质性因子	10.0
平行黏结模量/GPa	13.5
平行黏结刚度比	3.0
平行黏结抗拉强度/MPa	20.55
平行黏结黏聚力/MPa	40.55
平行黏结内摩擦角/(°)	30.0

下载: 导出CSV 

| 显示表格

室内试验的应力–应变曲线如图4（a）所示,单轴抗压强度、弹性模量对比如图4（b）,（c）所示。从图4（b）,（c）可知,7个模型的数值模拟值和室内试验值^[14]均较为接近,证明PWW模型模拟水对岩石强度和变形特性弱化作用的适用性和准确性。

图  4  数值模拟与室内试验结果对比

Figure  4.  Comparison between numerical simulation and laboratory test results

下载: 全尺寸图片 幻灯片

3.   结果分析

3.1   胶结物非均质性分析

在室内试验当中很难定量研究岩石胶结物的强度分布对岩石力学性质的影响,因此本文中采用数值模拟具有一定的优势。

为了定量评价某一黏结面积系数下,岩石胶结物的非均质性对强度造成的影响,接下来以S_r=0.888为例,分析平行黏结半径乘子的非均质性对模型宏观力学性质的影响。图5为S_r=0.888,$\overline{\lambda }$=0.486时,平行黏结半径乘子的概率密度图。图6为不同均质性因子m下,模型的应力–应变曲线和宏观力学性质变化关系。

图  5  平行黏结半径乘子概率密度（S_r=0.888, $\overline{\lambda }$=0.486）

Figure  5.  Probability density of parallel bond radius multipliers

下载: 全尺寸图片 幻灯片

图  6  宏观力学性质与均质性因子的关系

Figure  6.  Relationship between macro-mechanical properties and homogeneity factors

下载: 全尺寸图片 幻灯片

由图6可知,平行黏结半径乘子的非均质性对模型的单轴抗压强度和弹性模量均有很大影响。当均质性因子由1增加至40时,模型均质性增加,单轴抗压强度和弹性模量由21.7 MPa,13.9 GPa增加至25.8 MPa,15.2 GPa,分别增加了18.9%和9.4%。单轴抗压强度和弹性模量随均质性因子m的变化关系符合指数函数关系（R²分别为0.98,0.94）。该结论和其他学者所做的数值模拟^[21]或者理论推导^[22]的结果相一致。

图7（a）,（b）分别为不同黏结面积系数、不同均质性因子下岩石的单轴抗压强度和弹性模量。

图  7  均质性因子m和黏结面积系数S_r耦合对岩石宏观力学性质的影响

Figure  7.  Effects of homogeneity factor and bonded area coefficient coupling on macro-mechanical properties of rock

下载: 全尺寸图片 幻灯片

从图7（a）可知,岩石的黏结面积系数和非均质性均对岩石的单轴压缩强度和弹性模量存在一定的影响。当黏结面积系数较高时,单轴压缩强度对非均质性更为敏感,如S_r=0时,均质性因子从1增加至40,单轴压缩强度增加了47.4%；S_r=1.0时,均质性因子从1增加至40,单轴压缩强度增加了55.6%。但是黏结面积系数较大时,提高均质性因子,单轴压缩强度的增加值较小,如S_r=0,1,均质性因子从1增加至40,单轴压缩强度分别增加了18.1,8.4 MPa。从图7（b）可知,当黏结面积系数较低时,弹性模量对非均质性更为敏感。例如S_r=0时,均质性因子从1增加至40,弹性模量增加了45.7%。

本文所提出的PWW模型,由于耦合了Weibull分布,能够方便地研究岩石胶结物的非均质性对宏观力学性质的影响,为研究黏结面积系数和非均质性共同影响下岩石的宏观力学性质提供了有力的数值模拟工具。

3.2   能量耗散

能量耗散用于岩石内部微裂纹的萌生和扩展,是引起岩石内部损伤的本质原因。在PFC2D中,通过FISH函数监测模型在加载过程中各能量（耗散能、边界能、总应变能、摩擦能和动能）的变化^[15]。接下来以黏结面积系数S_r=0.432,m=10为例,分析加载过程中模型的能量耗散规律,如图8所示。

图  8  模型能量耗散与应力–应变关系(S_r=0.432)

Figure  8.  Energy dissipation and stress-strain relationship of model

下载: 全尺寸图片 幻灯片

由图8可知,在加载初期,由于上下墙体移动输入模型的边界能全部转化为模型的总应变能,耗散能、颗粒的动能和摩擦能几乎为零。随着轴向应力–轴向应变曲线进入线弹性阶段,模型开始萌生微裂纹,裂纹扩展过程中接触所储存的应变能逐渐转化为耗散能。边界能和总应变能增加速率增加,模型能量逐渐积累。轴向应力–轴向应变曲线进入峰前塑性屈服阶段后,边界能和总应变能曲线逐渐分离,耗散能开始缓慢增加,裂纹处于不稳定扩展阶段。模型进入峰后加速破坏阶段后,总应变能达到峰值,此时模型内接触所储存的能量达到最大值,随后应变能迅速释放,耗散能则急剧增加直至模型破坏。整个加载过程中,总应变能先增加后减小,耗散能则先缓慢增加后急剧增加。

表1中7个模型的总应变能变化如图9所示。

图  9  总应变能与黏结面积系数的关系

Figure  9.  Relationship between total strain energy and bonded area coefficient

下载: 全尺寸图片 幻灯片

由图9可知,随着黏结面积系数的增加（即岩石浸水时间的增加）,模型所储存的总应变能的总量和增速逐渐降低,总应变能峰值随黏结面积系数的增加大致呈线性下降。干燥状态下（S_r=0.0）总应变能峰值为331.1 J,饱和状态下（S_r=1.0）总应变能峰值为90.9 J,较前者下降了72.5%。由此可见,水对岩石存储能量的能力影响非常大。水的软化作用会使得宏观裂纹易于贯通,且贯通断裂面抵抗外力变弱,更加有利于能量释放^[23]。随着黏结面积系数的增加,岩石吸收的总应变能减少,这也是在地下工程的建设过程中经常采用岩层注水预防冲击地压的内在机制^[24]。

3.3   微裂纹倾角

微裂纹的倾角反映了模型在外部荷载下应力场的分布。倾角以水平向右为0°,逆时针旋转,将倾角分为18个组。本文表1中建立的7个模型的微裂纹的倾角分布如图10所示。

图  10  微裂纹的倾角分布

Figure  10.  Distribution of dip angle of microcracks

下载: 全尺寸图片 幻灯片

由图10可知,水对微裂纹的倾角影响较大。不同黏结面积系数下,平行黏结模型的黏结面积不同,各个接触的作用力不同,其模型的应力场分布也不同,微裂纹倾角的分布也随之不同。模型在干燥状态下（S_r=0）,微裂纹倾角集中于80°～100°,微裂纹平行或近似平行于轴向加载方向,这与在室内试验^[25]中观察到的花岗岩的结果类似。随着黏结面积系数的增加,微裂纹倾角的分布范围逐渐增加,在饱和状态下（S_r=1.0）倾角分布在70°～120°,表明水的存在使得微裂纹倾角分布更加分散。

在干燥状态下,岩石内部结构致密,孔隙和原生微裂隙较少。在轴向荷载下,接触在拉应力作用下发生拉破坏,微裂纹倾角集中于90°附近。在PFC2D中,岩石浸水时间为90天对应的模型（S_r=1.0）为平行黏结半径乘子特征值${\overline{\lambda }}_{\text{w}}$是0.421,接触黏结面积为12.70 m²,比干燥状态下的模型（S_r=0,平行黏结半径乘子特征值为1.0）的黏结面积（30.14 m²）减少了57.9%。室内试验^[26]中也观察到页岩浸水以后矿物胶结作用减弱,矿物之间的接触形式从面-面接触逐渐转变为点-面接触。当岩石颗粒黏结松散,疏松多孔以后,岩石内部力学性质薄弱处增加。在外力作用下,微裂纹的空间分布和倾角分布也更加分散。

3.4   破坏模式

图11是表1中7个模型在轴向应力下降至70%时的微裂纹分布图（红色线为素描的破坏面）和接触力分布图。

图  11  破坏模式与接触力分布（接触力单位为N）

Figure  11.  Failure modes and distribution of contact force

下载: 全尺寸图片 幻灯片

由图11可知,不同黏结面积系数的模型呈现出不同的破坏模式。当黏结面积系数较小时（如S_r=0,0.432）,岩石发育单个主破裂面,从右上角贯通至左下角,次生破裂面较少；当黏结面积系数较大时（如S_r=1.0）,由多条次破裂面组合而成的宏观破裂面导致模型失去承载能力,模型破碎程度更高。由于水对岩石内部细观结构有一定的软化和润滑作用,使得孔隙增加,原生裂隙更易扩展,随着黏结面积系数的增加,破裂面更加密集,且主破裂面的贯通性增强。

4.   结论

（1）揭示了岩石胶结物的非均质性与岩石宏观力学性质之间的关系。随着均质性因子m的增加,模型均质性增加,单轴抗压强度和弹性模量也随之增加。单轴抗压强度和弹性模量随均质性因子m的变化关系符合指数函数关系。同一均质性因子下、不同黏结面积系数时,单轴抗压强度和弹性模量对黏结面积系数的敏感性不同。

（2）随着黏结面积系数的增加,岩石所储存的总应变能的总量和增速逐渐降低,总应变能峰值随黏结面积系数的增加大致呈线性下降。

（3）岩石在干燥状态下（S_r=0）,微裂纹倾角集中于80°～100°,随着黏结面积系数的增加,微裂纹倾角的分布范围逐渐增加,在饱和状态下（S_r=1.0）倾角主要分布在70°～120°,水的存在使得微裂纹倾角分布更加分散。

（4）随着黏结面积系数的增加,颗粒间接触的黏结面积逐渐减小,接触黏结面积小的地方更易产生应力集中现象。水逐渐填充原生微裂隙和孔隙,岩石储能能力下降,岩石萌生微裂纹所需的应力也相应降低,微裂纹更易扩展。