0.   引言

岩石裂隙是地下岩体中普遍存在的薄弱部位,深刻影响着岩体的力学与渗流特性,分析含裂隙岩石在外力作用下的受压变形与损伤破坏规律是进行各类与岩石裂隙力学行为相关研究的基础。自然界的岩石处于干燥、部分或完全饱和状态,在流体的物理化学作用下,裂隙的微观结构及黏结特性将发生变化,从而改变岩石裂隙的宏观力学性质。因此,研究不同干湿状态下的岩石裂隙受压变形破坏特征具有重要意义。

岩石裂隙的变形破坏特征首先取决于岩石基质的力学特性。在以往的研究中,学者们通过大量的试验定量分析了水对各种岩石力学特性的影响。例如,Erguler等^[1]测试了含黏土岩石在不同含水量下的力学特性,结果表明相比于干燥的含黏土岩石,饱和岩石的单轴抗压强度（UCS）、弹性模量和抗拉强度分别降低了90%,93%和90%；Rajabzadeh等^[2]以碳酸盐岩为试验对象进行了单轴抗压试验,发现在饱和状态下这3项力学参数相较于干燥状态下降了约70%；花岗岩等结晶岩则因其良好的抗水降解特性,力学性质受水的影响较小,在饱和状态下的抗压强度降低率一般小于10%^[3]。由此可见,由于各类岩石之间矿物组分和胶结状态的不同,其力学特性受水的影响程度也存在显著的差异。

在岩石裂隙受压变形行为方面,以往研究通过大量的力学试验,发现并定量描述了在法向应力作用下裂隙的非线性变形行为,建立了经验模型描述该行为^[4]。然而,通常的物理试验只能获取裂隙宏观的应力-变形关系,缺乏对裂隙内部细观接触变形破坏过程的精确描述。相比于物理试验,数值方法能宏细观结合地分析岩石裂隙的受压变形过程。目前针对该问题,直接有限元法^[5-6]和基于半空间法的边界元法^[7-9]是两种较为常用的数值模拟方法,且半空间法具有计算效率更高的优势。Hopkins^[8]提出了一种基于力平衡的数值模拟计算方法,该方法使用不同高度的圆柱体模拟粗糙的裂隙表面,同时考虑了裂隙上下表面以及圆柱形微凸体的弹性变形。相较于早期的Greenwood- Williamson(GW)模型^[10],该方法考虑了粗糙裂隙中相邻接触体之间的相互影响,从而具有更高的计算精度。另一种基于Boussinesq解的半空间方法将裂隙表面视作连续面,可以更为真实地反映裂隙面的几何特征,并且能够实现对裂隙受压过程中局部塑性变形的模拟^[7]。到目前为止,大量的数值模拟研究只考虑了接触体的弹性变形,且一般采用数值生成的裂隙面。为精确描述真实岩石裂隙的受压变形破坏行为,Kling等^[11]基于花岗闪长岩岩芯裂隙面在0～10 MPa的法向应力范围内的压缩试验数据,进行了弹塑性接触模型的校准和验证；Zou等^[12]选取汨罗麻岩为对象,将裂隙压缩试验结果与使用半空间弹塑性接触模型计算所得的数值模拟结果进行了对比分析,从法向应力-位移曲线与裂隙粗糙面损伤区等角度分析并验证了数值模型的准确性。

大量试验研究已经证明岩石的破坏产生声发射,通过检测和统计声发射事件可以定量描述岩石内部的断裂破坏过程。对岩石声发射的研究主要集中在3个方面：①对岩石损伤过程中声发射事件率、能量率和振铃计数等声发射基本参数的研究。Rudajev等^[13]通过花岗岩单轴压缩试验,研究了岩石破坏损伤过程的声发射规律与参数；陈国庆等^[14]利用声发射事件率和能量率分析了开放型岩桥裂纹贯通机理及脆性破坏特征；Liu等^[15]研究不同温度下煤岩的力学特性、声发射能量和振铃计数,分析了岩石的损伤破坏行为。②利用地震震级理论、信息熵理论、分形理论和频谱分析等数学方法研究岩石损伤过程中的声发射特征参数。龚囱等^[16]利用地震震级理论研究了红砂岩蠕变破坏过程中的声发射b值特征；王春来等^[17]利用信息熵理论获得了花岗岩破坏过程中声发射事件空间熵值和时间熵值；Li等^[18]通过大理岩拉伸试验获得了声发射信号的双主频分布特征,讨论了声发射信号频率与岩石微破裂之间的关系。③利用岩石声发射事件进行岩石损伤定位。赵兴东等^[19]采用Geiger定位算法研究了含预制裂纹花岗岩损伤过程中微裂纹的三维空间演化规律；刘建坡等^[20]采用单纯形定位算法研究了预制孔岩石破裂过程中裂纹的时空演化规律；Dong等^[21]提出一种无测速三维声发射定位方法,并通过试验验证了算法的可靠性和精度。目前岩石声发射研究多集中于完整试样的声发射基本参数、特征参数和定位算法等方面,对岩石裂隙面受压破坏声发射响应规律的研究尚不深入。

综上所述,对岩石裂隙受压变形破坏特征的研究缺乏宏观应力-位移关系与细观接触变形破坏相结合的深入探讨,较少考虑岩石干湿状态造成的影响。针对这些问题,本文基于自主开发的裂隙面接触弹塑性分析方法^[12],开展干燥与饱和状态下含裂隙花岗岩和细砂岩的单轴压缩试验与弹塑性接触数值计算,结合不同状态下两种岩石的声发射响应特征和数值计算结果,分析裂隙受压变形过程中塑性变形的增长趋势和声发射响应规律,揭示塑性变形均值增量与应力之间的定量关系。

1.   研究方法

1.1   试验方法与流程

准备4个尺寸为φ50 mm×100 mm的圆柱形花岗岩与红砂岩标准试样。在试样半高位置切一道宽2 mm,深1 mm的环形凹槽,使用万能材料试验机沿环形凹槽将岩样横向劈裂形成新鲜的裂隙。将岩样沿中轴线旋转一定角度后直立放置,形成具有一定开度的单裂隙岩石试样,如图1所示。由于劈裂产生的裂隙面咬合较好,法向加载时裂隙充分闭合,其力学行为与完整试样差别较小。因此,这里采用旋转试样的方式让裂隙面形成一定程度的错位,从而更充分地研究接触体的变形破坏行为。分别制备干燥与饱和状态下的花岗岩和红砂岩试样。干燥试样采用烘干法制备；制备饱和试样时,将试样浸置于纯净水中,利用真空泵连续施加24 h的负压,以达到完全饱和状态,随后立即开展压缩试验。另准备同批次的花岗岩和红砂岩标准试样各6个,按干燥和饱和状态分为2组,分别进行P波波速测试和单轴压缩试验,获得干湿状态下两种岩石的平均单轴抗压强度UCS、弹性模量E、泊松比$\nu$、P波波速等基本参数,如表1所示。

图  1  岩石试样和裂隙错位情况

Figure  1.  Rock samples and dislocated fractures

下载: 全尺寸图片 幻灯片

表  1  岩样的基本性质

Table  1.  Fundamental properties of rock samples

试样编号	UCS/MPa	弹性模量/GPa	泊松比	P波波速/(m/s)	旋转角度/(°)
G-D	148.94	47.33	0.24	3179	6.11
G-S	133.20	44.99	0.23	4292	6.80
R-D	71.45	15.10	0.40	2689	3.49
R-S	39.98	11.20	0.35	3022	14.77
注：G为花岗岩,R为红砂岩,D为干燥,S为饱和。

下载: 导出CSV 

| 显示表格

试验系统主要包括轮廓测量系统、加载系统和声发射信号采集系统。采用分辨率为1 μm的Keyence VR3000高精度三维轮廓仪测量仪对压缩试验前后的试样裂隙表面进行扫描,获得表面几何形貌数据。采用MTS815电液伺服试验机在恒定应力速率控制模式下对裂隙试样进行分级加载,具体加卸载方案如表2所示。由于饱和红砂岩的单轴抗压强度较低,为保证试验后岩样的完整性,将其最大应力值降低为5 MPa,同步的加载速率也相应降低。声发射系统包括美国物理声学公司（PAC）生产的Express8采集卡（最大采样率10 MHz）,前置放大器（增益40 dB）和R6a声发射传感器,设置系统触发阈值为50 dB,采样率3 MHz,采样长度1024个点。6个声发射探头位置坐标如图2（b）与表3所示。安装完成后的试样如图2（a）所示,在加卸载过程中,采用柔性密封装置保证饱和试样中的裂隙内充满水,并实时同步采集试样的荷载、位移和声发射数据。

表  2  试样加卸载方案

Table  2.  Loading and unloading plans

试样编号	加载速率/(MPa·s-1)	目标应力/MPa	卸载速率/(MPa·s-1)	保载时间/s
G-D	0.1	10	0.1	45
G-W	0.1	10	0.1	45
R-D	0.1	10	0.1	45
R-W	0.05	5	0.05	45

下载: 导出CSV 

| 显示表格

图  2  试验设备及传感器设置

Figure  2.  Diagram of experimental equipment and settings of sensors

下载: 全尺寸图片 幻灯片

1.2   声发射定位算法

Geiger算法是目前用于计算岩石声发射源位置的有效方法之一^{[19, 22-23]}。该算法利用给定的初始点,通过多次迭代得到最终结果,每次迭代中,利用最小二乘法计算一个修正向量,再将修正向量加到上次迭代的计算结果中,从而获得一个修正后的新试验点,并判断这个试验点是否满足设定要求。若满足,认为该点即为声发射源坐标,若不满足则继续进行迭代。每次迭代的结果由下式获得：

式中,x,y,z为试验点坐标,t为事件发生时间,(x_i,y_i, z_i)为第i个传感器的位置,t_i为纵波到达第i个传感器的时间,v_p为纵波波速。纵波到达各传感器的时间t_oi,可用试验点坐标计算出的到达时间的一阶泰勒展开式表示：

式中,t_oi为第i个传感器获得的纵波到达时间,t_ci试验点坐标计算出的纵波到达第i个传感器的时间。式（2）中,$\frac{\partial t_i}{\partial x}=\frac{(x_i-x)}{vR}，\frac{\partial t_i}{\partial y}=\frac{(y_i-y)}{vR}，\frac{\partial t_i}{\partial z}=\frac{(z_i-z)}{vR}，$$\frac{\partial t_i}{\partial t}=1$,$R={\left[{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2\right]}^{\frac{1}{2}}$。

对于N个传感器,可以得到N个方程,用矩阵形式可表示为

式中,$A=\left[\begin{array}{cccc}\frac{\partial t_1}{\partial x}& \frac{\partial t_1}{\partial y}& \frac{\partial t_1}{\partial z}& 1\\ \frac{\partial t_2}{\partial x}& \frac{\partial t_2}{\partial y}& \frac{\partial t_2}{\partial z}& 1\\ ⋮& ⋮& ⋮& ⋮\\ \frac{\partial t_n}{\partial x}& \frac{\partial t_n}{\partial y}& \frac{\partial t_n}{\partial z}& 1\end{array}\right]，\theta =\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\\ \Delta t\end{array}\right]，B=\left[\begin{array}{c}{t}_{\text{o1}}-t_{\text{c1}}\\ t_{\text{o2}}-t_{\text{c2}}\\ ⋮\\ t_{\text{on}}-t_{\text{cn}}\end{array}\right]。$利用最小二乘法求解式（3）可得

通过方程（4）求出修正向量,以(θ+Δθ)作为新的试验点,继续迭代,直至满足误差要求。

1.3   数值模拟方法

基于接触力学中的变分原理,以接触系统最小势能的概念为基础,提出了裂隙表面接触的弹塑性计算模型^[12]。该模型通过Boussinesq解^[24]将裂隙表面接触应力与位移相关联,可实现裂隙表面接触体的弹塑性变形计算。对于两个粗糙裂隙表面接触的弹塑性变形问题,系统的总互补势能V^*[7]的表达式为

式中,Ω为接触域,${\sigma }_{\text{n}}$为法向应力,ū_z(x,y)为接触区域内的复合表面法向位移,${\overline{u}}_{\text{z}}^{\ast }(x,y)$为接触区域内两个接触微凸体基于几何干涉的总位移,$\Delta u_{\text{z}}^{\text{p}}$为接触应力达到岩石基质硬度时,塑性变形区域的复合增量表面位移,即${\sigma }_{\text{n}}$>H,H的取值根据文献中基于UCS的预测公式得到^[25]。利用数值方法求解方程（5）时,需将接触面离散成矩形单元结构网格,该网格与轮廓扫描获得的裂隙面形貌数据一致。为解决线性互补问题,使用约束共轭梯度算法^[26]迭代求解接触压力值,设置收敛条件${\sigma }_{\text{n}}$>0且Σ${\sigma }_{\text{n}}$=F。

将轮廓测量仪获得的裂隙表面三维几何模型进行等间距（0.1 mm）点云划分,构建三维数值裂隙模型。如图3所示,所建立的高精度三维数值裂隙模型较好地反映了实际裂隙的表面特征。在数值计算中根据表1设置不同状态下岩石的力学参数,固定下表面,按照$\Delta {\sigma }_{\text{n}}$=0.1 MPa逐步将应力加载到上表面,记录加载过程中裂隙面法向位移、法向应力、塑性破坏区域与面积等数据。

图  3  裂隙面几何形貌

Figure  3.  Morphological characteristics of fracture samples

下载: 全尺寸图片 幻灯片

2.   法向应力作用下裂隙变形破坏行为

2.1   裂隙变形曲线和破坏区域分析

图4为单轴压缩试验和数值计算得到的花岗岩与红砂岩试样的法向应力-法向位移关系曲线。

图  4  试验与数值计算得到的法向应力-法向位移曲线对比

Figure  4.  Comparison of experimentally and numerically obtained normal stress-normal displacement curves

下载: 全尺寸图片 幻灯片

结果表明,裂隙受压闭合过程中,法向位移u与法向应力呈非线性增长关系,曲线斜率随法向应力的增大逐渐减小,与以往研究得到的基本规律一致^{[4, 27]}。以往的研究中常将岩石定义为纯弹性体,只考虑岩石裂隙受压过程中的弹性变形,忽略接触体的塑性变形,导致纯弹性计算（EM）得到的法向位移相较于试验结果偏小。本研究采用的弹塑性数值方法（EPM）考虑了接触体塑性变形,较弹性接触模型得到了更为准确的结果。

如表1所示,干燥状态下两种岩石的单轴抗压强度、弹性模量和泊松比均大于饱和状态。饱和花岗岩的单轴抗压强度、弹性模量和泊松比分别是干燥试样对应参数的89.4%,95.1%和96.3%,饱和红砂岩的这3项基本参数是干燥状态的56.0%,74.2%和87.8%。

经检测,所研究的花岗岩为致密块状构造,孔隙率低（0.9%）,岩石颗粒亲水性弱,所以岩石受水化作用的影响较小；而红砂岩属于颗粒支撑类型,孔隙-接触式胶结,孔隙率较大（9.8%）,且含有一定的碳酸钙矿物（约为10%）,因而水的软化作用较为明显。如图4所示,在相同应力条件下,细砂岩的位移量大于花岗岩,饱和状态下的位移量大于干燥状态时的位移量。使用基于标准试样测得的力学参数（表1）进行数值计算,得到的应力-位移曲线与试验曲线吻合良好,说明精确测试和获取岩石在不同干湿状态下的力学性质以及裂隙表面形貌是准确评价裂隙变形破坏特征的关键。

将扫描获得的单轴试验前后裂隙表面形貌点云数据相减,可得到表面的破坏区域,如图5所示；导出数值结果中的裂隙表面塑性变形量,将其绘制成高程图像,即为数值计算得到的破坏区域,如图5所示。试验测量和数值计算得到的裂隙表面破坏区域吻合较好,进一步验证了该数值方法对裂隙受压变形破坏行为模拟的有效性。图5中分别列出了各岩样在目标应力作用下裂隙面的接触面积比C_a,即裂隙上下面接触面积与其总表面积之比,该参数反映了接触破坏区域面积的大小。强度偏低的砂岩接触面积比更大,而水的软化作用可进一步降低岩石的各项基本力学参数,使岩石的脆性减弱,塑性增强,所以饱和岩样的接触面积比更大。脆性材料破坏容易引起周边颗粒的开裂、松动或脱落,试验后需将这些部分完全去掉再进行表面测量；相比之下,数值计算仅考虑接触面上的塑形破坏,因此得到的破坏区域相对试验测量略小。

图  5  试验测量与数值计算得到的表面破坏区域对比

Figure  5.  Comparison of experimentally measured and numerically calculated asperity damage areas

下载: 全尺寸图片 幻灯片

2.2   裂隙变形曲线和破坏区域分析

将加载过程中每一应力增量对应的接触体塑形破坏单元体积除以裂隙的接触面积,得到塑性破坏变形量均值增量u_p,其与法向应力的关系如图6所示。塑性变形均值增量与应力整体呈负相关,且随着应力的不断增大,其下降梯度减小。加载初期的塑性变形均值增量普遍较大,这主要由于加载初期裂隙面的接触面积较小,接触体尖端上产生应力集中容易产生破坏所致。度过加载初期后裂隙面接触面积不断增大,且接触体尖端已被压平,因而塑性变形量的增长随之放缓。

图  6  塑性变形均值增量-法向应力关系图

Figure  6.  Relationship between mean plastic deformation increment and normal stress

下载: 全尺寸图片 幻灯片

对塑性变形均值增量-法向应力图进行曲线拟合：

式中,u_p为塑性变形均值增量,e₀为裂隙初始（即刚稳定接触时）平均开度,${\sigma }_{\text{n}}$为法向应力,E为弹性模量,R为试样半径,C_a为初始接触面积比,C_ai为根据数值模拟中每一个步长上的结果得到的在不同应力条件下的接触面积比,RSD为裂隙平均开度的相对标准差。u_p/e₀与力学及几何参数之间的关系如图7所示。该公式量纲统一,反映了塑形变形增量随着应力增长而逐渐降低的趋势,且充分考虑了裂隙力学和几何参数的控制作用,可实现对裂隙受压过程中塑形破坏的定量评价。

图  7  u_p/e₀与力学及几何参数之间的关系

Figure  7.  Relationship among u_p/e₀, mechanical and geometric parameters

下载: 全尺寸图片 幻灯片

3.   基于声发射的裂隙受压变形破坏分析

3.1   声发射定位分析

根据声发射P波到时数据、试样的P波波速和声发射探头位置坐标,采用Geiger算法可获得声发射源位置。声发射P波到时采用噪声阈值法^[28]获得,根据背景噪声信号幅值的平均值确定声发射门槛阈值,利用声发射信号首次越过门槛阈值的时间向前线性反推确定P纵波到时。

图8给出了试样的三维声发射源位置和试验破坏区域对比图。从三维声发射源位置的X-Z、Y-Z投影可以看出四个试样的定位点均集中在Z坐标为50 mm附近,这与裂隙位置一致,表明受压破坏主要发生在裂隙处。由声发射源位置在X-Y平面投影点与试验破坏区域对比可知,投影点均散布于裂隙表面的破坏区域,两者在空间位置上吻合较好。花岗岩裂隙面粗糙度大,接触体数量少且脆性较高,微凸体脆性破裂产生的声发射事件较多且能量较大,因而定位效果更好。红砂岩则接触体脆性偏低,受压破坏时产生声发射事件数较少,定位效果相对较差。

图  8  声发射定位点与试验测量得到的表面破坏区域对比图

注：红色为三维定位点，青色为三维定位点的X-Y投影，蓝色为三维定位点的X-Z投影，绿色为三维定位点的Y-Z投影。

Figure  8.  Comparison between AE positioning points and experimentally measured asperity damage areas

下载: 全尺寸图片 幻灯片

3.2   岩石裂隙声发射响应规律

图9给出了加卸载过程中两种试样的声发射振铃计数、累计计数、信号幅值随时间的响应关系。应力加载阶段（时间约为0～550 s）,两种岩石试样的声发射振铃计数明显高于卸载阶段（时间约为550～1055 s）。当法向应力增加时,花岗岩试样振铃计数明显增大,当应力进入保载阶段,振铃计数相较于加载阶段变小,法向应力由1 MPa增至2 MPa过程中,振铃计数最大,声发射信号强度和频度较高,裂隙面上微凸体的破坏较剧烈。声发射振铃计数主要集中在应力加载初期,说明接触体损伤破坏主要发生在此阶段,与图6所示计算规律相符。从能量的角度,在施加应力初期,应变能主要集中于接触体尖端,随着应力增大,接触体以塑性破坏的形式释放应变能,并产生声发射事件。

图  9  岩石裂隙声发射响应特征

Figure  9.  AE response characteristics of rock fractures

下载: 全尺寸图片 幻灯片

试样累计计数随法向应力的增加逐步增大,当法向应力减小时,累计计数增长明显放缓,逐渐趋于平稳。红砂岩振铃计数和累计计数均小于花岗岩试样,这是由于红砂岩强度和脆性较弱,受压破坏过程中释放的能量较少,声发射事件较少所致。花岗岩声发射信号幅值主要集中于50～100 dB,红砂岩信号幅值集中于50～80 dB,由于花岗岩强度和脆性较大,接触体受压破坏产生的应力波能量较大,使得声发射信号幅值略高。

3.3   水对岩石裂隙面声发射响应特征的影响

岩石裂隙由劈裂产生,且旋转了不同的角度,因此干燥与饱和条件下同种岩石的裂隙面接触状态存在一定的差异,使得两种条件下岩石的振铃计数峰值分别出现于不同应力阶段,声发射响应规律表现出微小差异性。由图9可知,干燥状态两种岩石的声发射振铃计数、累计计数都大于饱和状态。这是由于当试样处于饱和状态时,裂隙面接触体受压破坏释放的应力波传播至固液界面时,会发生反射、折射,使得应力波能量被削弱,从而使声发射信号幅值不能达到声发射采集系统设置的门槛阈值,导致部分声发射事件无法被记录,使得振铃计数和累计计数降低^[29]。另一方面,水的软化作用使得饱和状态岩石刚度减小,屈服强度和脆性降低,岩石应变能释放和储能能力降低,导致声发射响应减弱。

3.4   岩石变形破坏过程中的声发射特性

图6中应力与塑性变形均值增量的关系和图9中干燥花岗岩声发射振铃计数具有较好的对应关系,说明在岩石裂隙受压过程中发生的变形破坏主要集中在加载初期,这与图4中法向应力-法向位移曲线反映的单裂隙试样变形过程一致。干燥花岗岩声发射振铃计数的峰值相较于塑性变形均值增量的峰值出现得稍晚,这是由于试验在单轴压缩试验开始阶段存在压密过程,此时位移量虽较大,但裂隙面接触面积较小,多是裂隙面接触体的尖端破坏,其破坏体积较小,对应的振铃计数较少。对于红砂岩试样,则因自身性质的原因导致其未表现出和花岗岩试样一样明显的规律性,但振铃计数的峰值也出现在加载初期且具有和花岗岩试样一样的变化趋势。

4.   结论

岩石裂隙受压后的塑性变形直接改变其几何形貌,从而对岩体的力学与渗流特征产生显著的影响,是造成工程中岩体渗透性降低,地热及油气资源开采效率下降的重要诱因。本研究对干燥与饱和状态下含裂隙花岗岩和红砂岩开展了单轴压缩试验与弹塑性接触数值计算,同步采集了受压变形过程中的声发射信号,综合分析数据得到3点结论。

（1）基于裂隙面接触的弹塑性数值计算与试验结果在法向应力-法向位移关系曲线上有较高匹配度；同时,试验测量得到的破坏区域位置、数值计算得到的破坏区域位置与声发射定位点三者之间基本一致,表明研究采用的计算和测量手段具有较高的精度。

（2）塑性变形均值增量与应力呈负相关,且其减小的幅度随着应力的增大而变小,表明破坏主要发生在加载初期较尖锐的接触体上,随着接触面积的增大,应力从集中变得分散,破坏量逐步趋于稳定。通过拟合,建立了两者之间的公式,定量表征了裂隙力学和几何特征对塑形破坏的控制作用。

（3）饱和状态下两类岩石的强度降低,声发射特征产生明显差异。水的软化作用使饱和状态岩石裂隙的振铃计数、累计计数等低于干燥状态,且花岗岩声发射响应强于红砂岩。振铃计数与计算得到的塑形变形量有较好的对应关系。

本文只对岩石裂隙受压过程中的塑性变形均值增量和声发射特性进行了研究分析,且选取的应力范围较小,所得到的经验公式虽符合试验及声发射响应规律,但其适用性还需进一步验证。今后将针对更多种类的岩石裂隙受压过程中的力学及声学特性进行更加深入的研究。