基于多参量判据的深地下工程岩爆倾向性研究

2022-04-20张恒源范俊奇郭佳奇石晓燕孙飞跃

高压物理学报 2022年2期

张恒源，范俊奇，郭佳奇,3，石晓燕，孙飞跃

（1. 河南理工大学土木工程学院，河南焦作 454000；2. 军事科学院国防工程研究院，河南洛阳 471023；3. 河南省地下工程与灾变防控重点实验室，河南焦作 454000）

自1738 年英国South Stafford 锡矿岩爆被首次报道以来，世界上已有南非、美国、中国等众多国家发生了不同程度和规模的岩爆[1–4]。岩爆是高地应力条件下地下工程开挖过程中产生的一种突发性动力失稳地质灾害，在水利水电、交通隧道等工程施工过程中时有发生[5–6]。随着国内外交通和经济的快速发展，隧道以及矿山等深部地下工程的开挖深度逐步增加，“三高一扰动”施工环境日益恶化，岩爆已成为深地下工程长远发展过程中难以避免的地质灾害之一[2,7]。因此，正确认识岩爆灾害机制、准确预测岩爆发生可能性及强度，对保证施工安全和确保工程进度意义重大。

国内外专家学者从多种角度对岩爆及其判据进行了研究，提出了应力判据、能量判据、脆性判据以及综合判据等岩爆预测方法，并将判据合理地应用在工程案例中，为地下工程减少了损失。现有的应力判据主要根据围岩应力状态对岩爆进行预测，如Turchaninov 判据[1]、Russense 判据[8]等，并且此类判据主要利用岩石最大切应力、最大主应力、单轴抗压强度等指标，因此此类判据在形式上较为相似，主要区别体现在不同工程下岩爆分级临界值的不同，如Barton 判据[9]和陶振宇判据[1]的表达式相同但是取值并不相同。脆性是影响岩石在荷载作用下力学行为和破坏特征的重要特性[10]，岩石在发生岩爆破坏的过程中，脆性破坏是最主要的特征之一，一些学者在此基础上提出了岩石脆性指标[11]、脆性系数[12]、变形脆性系数[13]等脆性判据。然而岩爆并不是由单一的某个因素所决定的，单独将某种因素作为判断标准显然不够合理，一些学者基于此提出了考虑多因素影响的岩爆判据，如谷明成等[14]结合秦岭隧道岩爆具体情况，提出了多因素综合判据；Zhang 等[15]提出了五因素岩爆判据；邱士利等[16]通过分析锦屏二级水电站岩爆案例，提出了多因素岩爆判据RVI。相对于单一因素判据而言，综合判据考虑因素更全面，具有明显的优势。此外，Kidybiński[17]提出了岩石弹性变形能指数Wet；Wiles[18]针对矿柱岩爆问题提出了局部能量释放密度（local energy release density，LERD），Beck 等[19]、苏国韶等[20]、邱士利等[21]、徐婕[22]分别在LERD 的基础上加以改进，提出了模拟地层功（model groud work，MGW）、局部能量释放率（local energy release rate，LERR）、相对能量释放指数（relative energy release index，RERI）、岩爆能量释放率（rock burst energy release rate，RBERR）等能量判据，并在工程中加以应用；陈卫忠等[23]通过室内试验提出了基于岩石实际储存能量与极限能量之比的岩爆判据；郭建强等[24]基于弹性应变能建立了岩爆烈度分级预测模型；宫凤强等[25]提出了剩余弹性能指数岩爆判据指标；孙飞跃等[26]基于能量原理提出了岩爆倾向性指标。围岩内部积聚能量的迅速释放引起岩爆，因此从能量方面考虑能真实反映岩爆本质。对以上判据进行深入分析发现，现有岩爆判据大多考虑工程岩体处于三向受压的情况，忽略了岩体可能存在的双向受压单向受拉的情况。根据谢和平等[27]提出的单元整体破坏准则可知，岩体单元可释放应变能与表面能临界值的关系是岩体发生破坏的依据，但目前文献鲜有与岩体表面能相关的判据。综合判据和能量判据都可反映岩爆过程中部分因素造成的影响并具有一定的适用性，因而，若能将能量判据与综合判据相结合，则判据所考虑的岩爆因素会更加全面。

鉴于此，本研究依据岩石破坏过程中的能量转化机制和单元整体破坏准则，在充分考虑应力状态的基础上，提出一种多参量岩爆判据；基于三维离散元平台，对上述判据进行二次开发，分析不同工况下围岩应力、能量和岩爆判据指标的分布特征，探究深地下工程在开挖扰动作用下围岩的动态响应规律；将建立的岩爆预测方法应用于锦屏二级水电站4#引水隧洞岩爆实际工程案例，验证其可靠性和适用性。

1 基于能量原理的岩爆倾向性判据

1.1 能量转化机制

式中：Ei为卸载弹性模量，µ为泊松比，E0为初始弹性模量， σ1为主应力， σ2为中间主应力， σ3为最小主应力。

图1 单位体积中的能量耗散Uid 和可释放应变能Uie 的量值关系Fig. 1 Quantitative relationship between released strain energy Uid and dissipated energy Uie per unit volume

1.2 岩爆倾向性判据

1.2.1 受压情况

上述岩爆判据意义明确，简单实用，能够较合理地定量判别深部地下工程的施工过程中岩爆地质灾害的发生程度以及位置范围，更重要的是该判据全面考虑了单元体受力的各种状况，因而，在进行深地下工程岩爆灾害的模拟预测分析时采用该判据具有十分重要的意义。

2 深地下工程岩爆灾害的三维离散元模拟

通过室内试验研究岩爆的费用高、耗时长，同时还可能存在不可控变量的干扰等缺点。相比于传统物理试验，数值模拟的发展为岩爆研究提供了极大的帮助。岩爆现象的发生受多种因素影响，即使完全相同的岩体，在不同工况下，其发生岩爆的可能性并不相同。因此，开展不同工况下开挖的洞室围岩能量及岩爆倾向性的变化研究将加深对岩爆孕育过程的认识。天然岩体中存在节理、断层和裂隙等，这些不连续面在加载过程中对岩体的力学性质、变形特征等有重要影响。离散元方法可以将岩体划分为离散岩块和节理，岩块之间可以分离并产生相对滑动，能够较为真实地模拟岩体离散性以及岩爆引起的不连续大变形特征，因此，本研究采用以非连续介质为基础的离散元法进行模拟。

FISH 语言是3DEC 自带的编程语言，很多程序命令无法实现的功能均可通过FISH 语言进行编程实现。本研究利用FISH 语言编写了相关程序，得到了应力、能量、岩爆判据Crs的分布云图，实现了岩爆倾向性范围和程度的可视化。本节通过对不同埋深（600、1200、1800 m）、不同侧压力系数（1.0、1.5 和2.0）下的洞室进行三维离散元模拟，分析地下工程围岩的应力、能量演化规律和岩爆倾向性分布规律。

2.1 离散元模型的构建

2.1.1 数值模型及边界约束条件

考虑到圣维南原理和洞室开挖的影响范围，建立的模型几何尺寸为50 m×60 m×20 m（高度×宽度×厚度），数值模型见图2。该计算模型的上边界为应力约束边界，下边界、前后和左右边界均为位移约束边界，如图2(a)所示。洞室开挖尺寸如图2(b)所示。

图2 隧道数值计算模型Fig. 2 Numerical calculation model for tunnel

2.1.2 本构模型及材料力学参数

本构模型主要是描述单元的应力-应变关系，对于弹塑性本构模型的选取则主要考虑单元的屈服准则和流动法则。为真实反映围岩的受力状况，本研究的本构关系采用描述硬岩力学行为的Mohr-Coulomb 屈服准则，该准则的破坏包络线与剪切屈服函数以及拉应力屈服函数相对应，是一个与拉伸破坏相关的流动法则[31]。在进行数值模拟时，不同埋深下采用的岩体力学参数由三轴试验获得，参数见表1，其中，h为埋深，E为弹性模量。

表1 岩体物理力学参数Table 1 Physical and mechanical parameters of rock mass

2.2 数值模拟方案和监测点的布置

在数值分析中，为模拟不同埋深及不同地应力条件下的围岩响应特征，洞室埋深h选取600、1200、1800 m，侧压力系数K选取1.0、1.5、2.0，模拟方案如表2 所示。

表2 数值模拟工况Table 2 Numerical simulation conditions

为揭示能量的时空分布规律，在拱顶、拱底与拱肩处布置了4 个监测点（A、B、C、D），监测点位置如图3 所示。

图3 数值计算模型中监测点位置Fig. 3 Locations of the monitoring points in the numerical model

2.3 模拟结果及分析

2.3.1 应力场演化分析

围岩应力状态的改变是发生岩爆破坏的诱因之一。在洞室开挖过程中，原有地应力平衡状态被打破，围岩应力重新分布。在工程设计和施工阶段，围岩应力重分布后的应力状态是预防岩爆发生所需考虑的重要因素。此外，研究不同埋深、不同侧压力系数下洞室开挖后围岩应力状态的重新分布情况对进一步了解岩爆的孕育与发展有重要意义。

主应力差是最大主应力和最小主应力之差，即∆σ=σ1−σ3。图4 为不同侧压力系数、埋深下洞室开挖后的围岩主应力差云图。

图4 不同工况下数值模拟的主应力差云图Fig. 4 Contour maps of the principal stress difference under different conditions obtained by numerical simulation

由图4 可知，不同工况下主应力差的分布范围产生了较大的变化，侧压力系数越大则埋深越大，主应力差的影响范围越大。3 倍洞径内主应力差较小值多集中在洞室拱底及拱脚处，而主应力差较大值多集中在洞室拱顶。对比不同工况可知，埋深以及侧压力系数主要影响主应力差较小值（小于5 MPa）的分布范围。在埋深为600 m 的情况下，侧压力系数由1.0 增加到2.0 时，洞室附近主应力差较小值分布范围由1 倍洞径增加到2 倍洞径；当埋深为1200 和1800 m 时，随着侧压力系数的增加，主应力差较小值的分布范围变化相对较小。说明当埋深为600 m 时侧压力系数为主要影响因素，而当埋深大于1 200 m时侧压力系数的影响减弱。当侧压力系数确定时，主应力差较小值的分布范围随着埋深的增加而增大，且当侧压力系数从1.0 增加至2.0 时，埋深对主应力差较小值分布范围的影响效果没有明显减弱。

2.3.2 能量演化过程分析

为计算岩爆判据指标Crs的分布情况，采用FISH 编程语言编写计算函数，监测所有计算单元的弹性应变能演化过程，获得了不同埋深与侧压力系数下洞室开挖后围岩的弹性应变能密度分布状态（图5）和弹性应变能密度的时空演变过程（图6）。

由图5 可知，不同工况下围岩的弹性应变能密度较大值多集中在洞室拱顶、拱脚处，峰值多出现在洞室拱顶位置，最大值为48.566 J/m3。在相同侧压力系数下，随着埋深的增加，洞室拱肩、拱脚处应变能密度开始增加，尤其在埋深为1800 m、侧压力系数为2.0 时，拱顶、拱肩及拱脚将积聚大量的弹性应变能，且分布范围较广。

图5 弹性应变能密度分布Fig. 5 Distribution of the elastic strain energy density

图6 为不同工况下洞室拱顶、拱底及拱肩4 个监测点的弹性应变能密度随时间的变化。由图6 可知，4 个监测点的弹性应变能密度随着时间的增加而增大，在经过一定的时间后，拱顶（监测点D）的弹性应变能密度明显高于其他位置，因此可以认为拱顶更容易发生岩爆。在侧压力系数一定时，随着埋深的增加，4 个监测点的弹性能应变能密度增加明显：当侧压力系数为1.0 时，随着埋深由600 m 增加至1800 m，监测点D的最大弹性应变能密度由26 MJ/m3增加至38 MJ/m3；侧压力系数为1.5 时，监测点D的弹性应变能密度由600 m 时的28.5 MJ/m3增加至1800 m 时的40 MJ/m3；侧压力系数为2.0 时，监测点D的弹性应变能密度由600 m 时的29 MJ/m3增至1800 m 时的42.5 MJ/m3。

图6 弹性应变能密度的时空分布Fig. 6 Spatial and temporal distribution of the elastic strain energy density

2.3.3 岩爆倾向性指标分布特征

图7 为采用FISH 编程语言编写计算函数得到的岩爆判别指标Crs的分布云图。图7 显示，在工况1（洞室埋深600 m，侧压力系数1.0）中，洞室拱顶位置出现Crs>19，表明此时拱顶容易发生岩爆。对比图4和图5 可知，工况1 开挖后在3 倍洞径内洞室拱顶的应力差较大，且拱顶应变能密度较其他位置更大，因此可以认为此时洞室拱顶处的岩爆倾向性较大。而当侧压力系数确定时，随着埋深增加，洞室边墙、底板的Crs显著增大，拱顶、边墙与底板均有明显岩爆倾向性，拱肩处Crs增加不明显，岩爆倾向性较小。由于Crs具有多参量判据的性质，岩爆判别结果并不等同于单独主应力或能量场的变化特征。埋深确定时，不同侧压力系数下岩爆倾向性指标Crs的强弱分布位置大体一致，即岩爆发生的位置大致相同，但岩爆烈度存在差别。随着侧压力系数增加，洞室边墙、底板由中等岩爆转变为强烈岩爆的区域范围也在增加，发生强烈岩爆的可能性加大。洞室围岩的岩爆倾向性指标Crs与围岩的主应力差及能量场均呈对称分布，但主要分布位置以及分布特征仍有区别。

图7 岩爆倾向性指标 Crs 的分布云图Fig. 7 Distribution cloud map of the rockburst judgement index Crs

通过分析深部地下洞室围岩应力场及能量场的分布演化规律，研究了深部地下工程岩爆的演化特征，并将基于能量原理建立的岩爆倾向性指标Crs在数值模拟中加以应用。结果表明：岩爆判据指标Crs为评估岩爆倾向性提供了一种新的途径，但其是否适用于工程实例还有待进一步验证。

3 工程验证

锦屏二级水电站位于四川省凉山彝族自治州的雅砻江干流锦屏大河湾上，共有7 条隧道横贯锦屏山，包括4 条从景峰桥至大水沟的引水隧洞（1#、2#、3#、4#）、2 条平行交通辅助洞以及1 条施工排水洞。其中，引水隧洞之间的间距为60 m；2 条辅助洞间隔35 m，平行布置在引水隧洞南侧，与引水隧洞之间的距离为65 m。锦屏二级水电站引水隧洞全洞平均埋深1610 m，该工程位于西南高地应力区，其埋深大、高地应力等特点决定了发生岩爆等灾害的概率高。4 条引水隧洞自施工开挖至贯通竣工，频繁发生岩爆现象，出现施工机械损坏、工人受伤、工期延误等情况，严重影响工程进度。为了验证岩爆判据指标Crs的可靠性，对锦屏二级水电站4#引水隧洞岩爆实例进行数值模拟。

本岩爆案例发生于隧洞桩号K9+728～766 范围内的边墙，爆坑达2 m 深，岩爆的剧烈冲击作用使锚杆拉断且严重扭曲变形，洞室断面以及岩爆现场如图8(a)、图8(b)所示，地应力状态以及岩石力学参数如表3 和表4 所示[32–33]。其中，σx、σy、σz分别为x、y、z方向的主应力，τxy为作用于x面上沿y方向的切应力，τxz为作用于x面上沿z方向的切应力，τyz为作用于y面上沿z方向的切应力，cp为黏聚力峰值，cr为黏聚力残余值， ϕi为摩擦角初始值， ϕp为摩擦角峰值，ψ 为剪胀角。

表3 岩爆洞段的地应力状态Table 3 In-situ stress state in the cross section of the rockburst cavern

表4 岩体力学参数Table 4 Mechanical parameters of rock mass

图8 锦屏二级水电站4#引水隧洞断面尺寸及岩爆洞段现场Fig. 8 Dimension of the tunel 4# of Jinping II Hydropower Station and the rockburst cavern

隧洞开挖后的弹性应变能密度分布、岩爆判别指标分布以及岩爆破坏的模拟结果分别如图9(a)、图9(b)、图9(c)所示。由图9(a)可知，隧洞开挖后，聚集在围岩内部的弹性应变能以不同形式释放，洞室左侧弹性应变能随着距洞室距离的增加逐渐减小，在隧洞左侧边墙积聚的弹性应变能密度达到了最大值22.6 MJ/m3，表明此处的弹性应变能有可能超过表面能临界值并剩余部分弹性应变能，剩余的弹性应变能转化为动能，使岩体发生脱落、弹射等动力现象。由图9(b)可知，4#引水隧洞Crs的最大值在60 左右，出现在掌子面左侧边墙。对比图9(b)、图9(c)可知，Crs较大部位出现了岩块弹射现象，与岩爆实际发生区域基本一致；由Crs的判别式可知，此部位发生强烈岩爆，而实际爆坑深度达2 m，模拟强度与实际岩爆强度相符合。由此可见，根据Crs指标对岩爆进行判断是合理的，能够很好地反映岩爆的位置以及强度，具有良好的工程实用性。

图9 4#引水隧洞模拟结果Fig. 9 Simulation results of 4# headrace tunnel

高地应力下地下洞室开挖卸荷过程中，硬岩破坏的本质特征是开裂[34]，因此抑制硬岩内部裂口发展是避免岩爆等硬岩破坏的关键。运用本研究提出的岩爆判据对工程进行岩爆倾向性预测，根据预测结果采取相应措施可降低岩爆灾害的影响。目前，避免岩爆灾害的具体举措有优化开挖方案、围岩加固、应力解除等，如高地应力下硬岩地下洞室稳定性设计优化的裂化抑制法[35]可通过综合优化开挖方案、支护参数、支护时机以减少围岩开裂，从而减小或规避高地应力下地下洞室岩爆灾害。