0.   引言

随着科学技术的进步,深部岩体中的资源开采与工程建设越来越多。由于深部岩体中的高地应力,开挖/开采中的瞬时卸载容易引发动力效应,可能造成与静力问题不同的破坏模式^[1-2]。

国内外学者对于高地应力下隧道开挖所产生的动力效应进行了比较全面的研究。Miklowitz^[3]采用拉普拉斯变换的方法,求解了隧道瞬时开挖的弹性解析解。Cai^[4]分别采用FLAC以及Phase²模拟隧道开挖过程,经过对比后认为开挖过程的模拟必须考虑动力效应。Li等^[5]将不同卸载路径下的理论解与二维离散元模拟结果进行了对比,并研究了不同侧压力系数、卸载路径下裂纹数量以及能量的变化。Fan等^[6]通过基于能量的损伤判断方式,对比了准静态开挖和动态开挖所产生的损伤区大小,认为由于远场应力做功使得近场能量瞬间集聚,瞬时卸载的损伤区增大。Chen等^[7]分别计算了考虑钻爆法开挖的动力效应以及只考虑准静态应力重分布时,围岩损伤范围,并与锦屏二期工程隧道实测数据进行了对比,认为高地应力下的动态开挖效应不可忽视。Yan等^[8]研究发现当埋深超过800 m后,地应力迅速增加,原始地应力的卸载效应会超过爆破荷载的效应,在损伤区的形成中占主导地位。

现有的研究中大多只考虑了完整岩体中的开挖动力效应,但实际岩体中通常含有节理,而节理一般对围岩的破坏有很大影响。此外,目前研究中多采用有限元方法进行数值模拟,只能关注能量、强度、应变等宏观参量,需要借助宏观本构体现岩石动力特性,且难以模拟动力问题中因应力波动造成的微裂隙的渐进产生和扩展。因此,本文利用离散元从宏微观角度研究节理岩体中洞室开挖卸荷引发的动力破坏过程,分析连通率对于动力响应及破坏机制的影响。以期为高地应力下的节理岩体开挖分析提供有价值的参考。

1.   模型建立及模拟步骤

因隧道的轴向尺寸远大于截面尺寸,且沿轴向受力近似均匀,可简化为平面应变问题。随着深度的增大,地应力逐渐趋向于静水压力状态^[9],因此本文针对静水应力下的平面应变问题进行模拟。岩体中缺陷对工程的稳定性有显著影响,其尺度范围十分广泛,从微观（岩石晶体之间的裂纹、裂隙等）到宏观（节理、断层等）各种尺度。本文关注与隧道几何尺寸相同尺度的节理,且认为节理有一定走向。

1.1   微观接触模型

本文采用两种微观模型分别模拟岩石和节理的力学特性：其一为考虑胶结厚度与宽度的微观本构模型^[10],这一微观模型能再现无节理岩石真实破坏特性；其二为Ivars等^[11]提出的光滑节理模型,该模型通过生成节理表面以避免颗粒受到接触方向的影响,两个颗粒可以沿着接触面相互覆盖而产生相对滑移。节理可以根据是否设置黏结强度而分为无黏结节理与有黏结节理两种形式。本文采用无黏结的光滑节理模型模拟岩体中的节理。

1.2   微观参数的确定

砂岩是一种广泛分布在地层中的沉积岩,本文选用某砂岩并根据其宏观力学特性、粒度、孔隙率进行离散元制样和微观参数标定。最终微观参数的标定结果如表1所示,其相应的宏观参数如表2所示,由此可以计算得到岩体的横波和纵波的波速分别为1366,2477 m/s^[18]。为了模拟动力过程,还需要标定微观局部阻尼,本文利用PFC模拟弯曲激振的试验并按照宏观阻尼比为2%^[12]进行标定,得到微观局部阻尼为0.05。光滑节理模型的刚度采用其默认值,即颗粒之间的弹簧元件串联之后的刚度（法向刚度7.5×10⁹ N/m,切向刚度5×10⁹ N/m）,而摩擦系数与颗粒的相同。

表  1  DEM岩石试样模型参数

Table  1.  Material parameters of DEM rock specimens

类型	微观参数	数值
颗粒部分	颗粒密度$\rho$/(kg·m-3)	2700
	颗粒法向刚度$k_{\text{n}}$/(N·m-1)	1.5×1010
	颗粒切向刚度$k_{\text{s}}$/(N·m-1)	1×1010
	颗粒摩擦系数$\mu$	0.7
	颗粒抗转动系数$\beta$	1.2
胶结部分	最大胶结厚度$h_{\max }$/m	1×10-5
	胶结弹性模量$E_{\text{b}}$/N	2.3×108
	胶结抗压强度${\sigma }_{\text{c}}$/Pa	2.7×109
	胶结抗拉强度${\sigma }_{\text{t}}$/Pa	2.7×108
	胶结延伸率$s_{\text{p}}$	0.15

下载: 导出CSV 

| 显示表格

表  2  无节理岩石宏观参数

Table  2.  Macro-parameters of intact rock

对比	弹性模量[13]/ GPa	泊松比[13]	单轴抗压强度[14]/MPa	巴西劈裂强度[15]/MPa	内摩擦角[13]/(°)
试验	10	0.152	38	2.55	40.1
DEM	10.2	0.15	36	2.552	39.7
误差/%	2	1.32	5	0.07	1.00

下载: 导出CSV 

| 显示表格

1.3   模型试样生成

本文假设颗粒级配均匀分布,且根据砂岩粒度选择粒径为0.3～0.5 mm。通过分层欠压法^[16]生成均匀试样,试样内含有20万颗粒。由于砂岩的孔隙率通常在10%～20%之间^{[13, 17-18]},因此本文通过墙体伺服的方式施加围压将颗粒压至孔隙率15.3%,生成胶结并计算至平衡状态,得到完好的岩石试样。随后将岩石试样压至指定的地应力水平,生成节理并计算至稳定状态,得到DEM节理岩体试样。

采用岩桥距离$L_{\text{b}}$节理长度$L_j$、角度$\theta$、连通率$R_j=L_j/(L_j+L_b)$、间距$L_I$描述节理几何特征,如图1所示。

图  1  DEM模型示意图

Figure  1.  Schematic diagram of DEM specimen

下载: 全尺寸图片 幻灯片

1.4   模拟步骤

为了利用离散元准确模拟节理岩体的开挖动力响应过程,需要解决计算域边界设置的问题。

由于实际开挖过程中所产生的应力波能在无限域中向外传播,而数值模模拟范围有限,所以需要在人工截断处的边界上模拟实际应力波向外辐射的过程,根据简单的一维波动理论,任何一个截断处的动应力大小为^[18]

式中 ${\sigma }_{nd}$,${\tau }_s$分别为黏性边界所对应的法向,切向动应力；$\rho$为岩体的密度；$c_p$为纵波波速；$c_s$为横波波速；$v$为岩体在边界相应位置处的运动速度。

在考虑初始地应力的情况下,可以得到截断边界处任意一个颗粒上所受到的法向外力为

式中,${\sigma }_s$为初始的地应力,$d$为颗粒的粒径。当物体的变化速率十分缓慢的时候,则边界上的退化为准静态条件。

在模拟过程中,设置初始地应力为25 MPa的静水压力,并采用瞬间删除中心颗粒的方式近似模拟快速开挖过程,同时利用PFC内建的fish语言,在循环计算过程中插入人工黏性边界的计算过程以模拟应力波的辐射。

由于篇幅有限,因此本文将一部分变量设为定值（$\theta =0$,$L_t=R$以及$L_I=1.5R$）,分别对不同连通率（0.4,0.6,0.8）下的工况进行分析。

2.   数值模拟结果

2.1   动力与静力的计算结果对比

为了研究动力效应对于计算结果的影响,本文对比了动力和静力两种模式下的计算结果。静力模拟中,将局部阻尼设置为0.7,并对开挖部分的材料刚度分多个阶段进行软化^[4]来控制卸荷速率,当最后一个阶段的计算达到稳定后对其进行挖除,以此来模拟准静态变化过程。在动态效应研究时,设置局部阻尼为0.05,并且开启黏性边界。

图2展示了动力与静力问题的胶结破坏分布。对于动力问题,本文给出了胶结破坏的时空分布来分析破坏演化的过程（不同颜色的破坏发生在不同时刻）。而对准静态仅给出最后破坏情况。当考虑动力效应时,本文通过监测边界上颗粒的速度变化计算了应力波到达边界的时间点,在此时刻内（3×10^-5 s以内,下文简称阶段一）,卸载所产生的应力波在岩体内传播,并主导了胶结破坏的产生。由图2（a）可见,瞬时卸载所产生的动力效应（阶段一）,导致了围岩迅速产生损伤。对于连通率较大的情况（大于等于0.6）,应力波传播过程中可能导致节理端部产生裂纹,并扩展形成破裂区。这一破裂区可以看做是伪自由面,在破裂区以及其内部的岩体依旧具有一定的承载能力,但是主要依靠外部比较完好岩石进行承载；随后在距离洞室较远的地方,由于内层围岩发生破裂弱化以及节理尖端损伤,在经过一段时间后也将形成破裂区。不考虑动力效应时,由图2（b）可见,围岩仅产生局部少量的损伤,原有节理的尖端会发生裂纹的扩展,但是较少发生聚合并形成破坏区。

图  2  不同连通率胶结破坏分布图

Figure  2.  Distribution of cracks under different joint persistencies

下载: 全尺寸图片 幻灯片

为了进一步量化胶结破坏的数量并分析对比两种计算模式下不同工况的胶结破坏差异性,图3给出了不同工况下的胶结破坏数目。由图3可以看出两种计算模式下胶结破坏数量均以拉剪破坏为主,但是当连通率从0.6增大至0.8后,动力计算结果中的胶结破坏数量增长速度远大于准静态计算的结果。这意味着高地应力、高连通率的岩体中,卸载动力效应会导致后续产生严重的破坏。

图  3  不同连通率下胶结破坏数量

Figure  3.  Number of cracks under different joint persistencies

下载: 全尺寸图片 幻灯片

2.2   微观破坏模式分析

本节通过胶结破坏位置以及其周边的颗粒位移变化对连通率0.8的工况破坏模式进行微观分析。当考虑动力效应时,根据胶结破坏所处位置的不同可以大致将破坏分成两种形式,见图4。第一种为靠近洞口处的破坏（1.5R范围内）,如图4（a）,底板和顶板与节理之间岩体（位置一）发生大面积的损伤,这部分岩体受到周围岩体的挤压,并且可能向洞室中心飞溅；而洞室两帮（位置二）通常会由于节理端部的应力集中而发生小部分破坏。第二种为远离洞口处的破坏（1.5R范围外）,如图4（b）,其通常表现为层间节理首尾翼型裂纹的贯通,最终可能形成多层破裂区。

图  4  破坏模式示意图

Figure  4.  Schematic diagram of failure modes

下载: 全尺寸图片 幻灯片

图5,6详细展示了洞周和其他地方的破坏。其中小的黑色箭头表示了颗粒位移方向,大的红色箭头表示了颗粒整体位移方向。红色、绿色、蓝色线段分别对应了拉剪、纯拉、压剪胶结破坏,黄色线段代表了节理。图5（a）为洞周底板位置,由于节理面的存在使得周围岩体沿着节理面发生滑动,导致了中间的岩体受到了挤压从而产生严重的破坏。图5（b）为洞周两帮,由于节理端部的应力集中,使得与节理共面的围岩直接被压碎。

图  5  靠近洞口的破坏（考虑动力效应）

Figure  5.  Failures around tunnel (dynamic solution)

下载: 全尺寸图片 幻灯片

图  6  远离洞口的破坏（考虑动力效应）

Figure  6.  Failures away from tunnel (dynamic solution)

下载: 全尺寸图片 幻灯片

在远离洞口的位置处,相对较近处（图6（a））,颗粒均向着洞室中心移动,由此可以推断出此处的宏观裂纹是岩体之间的相对拉伸导致的,而微观上颗粒的胶结破坏模式大多为拉剪破坏。在相对较远处（图6（b））,通过裂纹两侧颗粒位移方向可以看出,左侧岩体向洞室中心移动,而右侧节理之间的块体发展为旋转块体,两者之间产生相对的拉伸与剪切,因此裂纹表现为拉剪裂纹,同时微观上的胶结破坏也是由拉剪主导的。

当采用静力计算时,其破坏的范围与程度均显著降低。图7为静力问题的洞口结果,与图5动力结果相比,可以发现准静态模式下,图7（a）中的岩块不再受到两侧的挤压而产生显著破坏,图7（b）中节理端部的挤压更为明显。图8为静力问题远离洞口的结果,与图6相比,准静态模式下,最终仅在节理端部产生少量破坏,颗粒整体位移呈现挤压的状态,胶结破坏中压剪百分比增大（图8（a））；颗粒位移仅在端部呈现出一定的拉剪趋势（图8（b））。

图  7  靠近洞口的破坏（不考虑动力效应）

Figure  7.  Failures around tunnel (static solution)

下载: 全尺寸图片 幻灯片

图  8  远离洞口的破坏（不考虑动力效应）

Figure  8.  Failures away from tunnel (static solution)

下载: 全尺寸图片 幻灯片

2.3   特定位置颗粒速度的时程变化分析

波在传播过程中遇到节理之后,会产生透射、反射,因此本文分别选取了多个观测点,分别观测不同连通率节理附近的颗粒速度相对于无节理时的速度变化,以分析节理对于动力效应传播的影响。

近孔口的节理及观测点的布置形式如图9所示。其中在y轴方向上布置了2组观测点,每组观测点包含了3个连续且间隔为0.1R的观测点。第一组位于两条相邻节理的岩桥中心,第二组观测点位于节理的中位线上。在x方向布置第三组观测点,其位置分别位于岩桥中心以及节理中心。以下用G代表组,P代表点,例如G1P1代表第一组第一个点（x=0,y=1.4R）,G2P2代表第二组第二个点（x=0.5R,y=1.5R）。

图  9  近孔口的节理布置示意图

Figure  9.  Schematic diagram of joints around tunnel

下载: 全尺寸图片 幻灯片

图10～12展示了3组观测点在不同连通率下的颗粒速度时程变化。第一组的观测点的位置所得到的颗粒速度如图10所示。在G1P1位置处,阶段一内由于开挖卸载所产生应力波向外传递,因此可以看到颗粒速度瞬间增大,并且有节理岩体的颗粒速度大于无节理岩体的颗粒速度,说明节理的存在放大了岩桥中心的动力效应；另一方面,该点的动力效应随着连通率的增大而增大。随后两个点的初始速度随连通率、时间的变化规律与第一个点相似,但是峰值速度随距离的增大而呈现递减。

图  10  第一组观测点颗粒速度时程变化

Figure  10.  Change of velocity of first group of particles

下载: 全尺寸图片 幻灯片

图  11  第二组观测点颗粒速度时程变化

Figure  11.  Change of velocity of second group of pariticles

下载: 全尺寸图片 幻灯片

图  12  第三组观测点颗粒速度时程变化

Figure  12.  Change of velocity of third group of particles

下载: 全尺寸图片 幻灯片

对比3个点的颗粒速度随连通率的变化可以发现,当连通率小于等于0.6时,其速度随时间快速收敛至0。而对于连通率为0.8的工况,由于出现了多个阶段的破坏,因此破坏所释放的能量导致其相应产生了多个波峰。

第二组观测点位于节理的中位线上（x=0.5R）。在G2P1位置处,含节理岩体的颗粒初始峰值随着连通率的增大而减小,但是均大于无节理的情况。进一步对比有节理与无节理的情况可以看出,无节理岩体的速度曲线经过了一个波峰后快速递减至零,而有节理岩体中的颗粒速度则维持了一段时间的振动,体现了节理的反射效应。

G2P2位于节理的中央,相比于无节理的情况,此处颗粒起始振动产生了明显的滞后现象。另一方面,颗粒起始峰值速度小于无节理的情况,并且随着连通率的减小而减小。由于光滑节理模型的刚度较小,因此颗粒振动的频率显著降低。

在G2P3处,其初始的波峰与G2P2类似,都具有滞后的特性。同时与G2P1相比可以发现,在无节理的情况下,其颗粒峰值速度随空间的衰减幅度较小；而对于存在节理的情况,其振动的峰值均发生了显著的衰减,说明节理面起到一个缓冲的作用,阻碍了应力波的传播。

在第三组观测点中,G3P2（x=2R,y=0）处的颗粒速度在卸载动力效应结束之后产生了一个巨大的波峰,这可能是由于颗粒自身的胶结发生破坏所导致的,由此可见胶结破坏所释放的能量会造成十分剧烈的影响。G3P2（x=1.5R,y=0）的变化趋势与G2P2类似,不再赘述。

3.   结论

本文采用离散元研究了高地应力下节理岩体中动力破坏机制。通过改变节理连通率,研究了节理密度对于胶结破坏数量以及颗粒速度的影响,并分析了宏微观破坏模式。得到了以下3点结论。

（1）在25 MPa的地应力水平下,准静态计算结果仅产生少量破坏,而动力模式计算所得到的岩体破坏范围与程度均显著增大。准静态计算结果仅在节理端部产生破坏,而动力模式计算结果产生了多条聚合形成的破裂带。

（2）当节理连通率为0.8时,静、动力计算所得结果仅在洞口两帮的微观破坏较为类似。动力问题中,底板受到两侧岩体的挤压破坏；远离洞口产生了翼型裂纹,并扩展贯通为宏观的拉伸、拉剪裂纹,微观上均为胶结拉剪破坏。

（3）节理面对传播中的应力波起到反射与阻碍的作用。因此,节理面与孔口之间的岩块将受到反射后的压应力波以及岩块两侧岩体的挤压作用,从而产生较大的破坏。而节理面后的岩体受到了节理面的缓冲作用,在卸载动力效应下的破坏较少。