基于CDEM的隧道卸压爆破及岩爆抑制效应模拟*

陈柯竹 李天斌 马春驰 张 航 高美奔

(①成都理工大学环境与土木工程学院，成都 610059，中国) (②地质灾害防治与地质环境保护国家重点实验室(成都理工大学)，成都 610059，中国)

0 引 言

岩爆是地下工程开挖或扰动过程中，高储能岩体通过脆性破裂快速释放应变能，导致围岩爆裂、弹射的动力现象(伍法权等，2010；冯夏庭等，2019)。岩爆具很强的突发性、随机性和危害性，严重影响地下工程施工的安全，可通过改变围岩应力状态、改变围岩物理力学性质、合理支护及施工等方式降低岩爆可能性(李天斌等，2016；王鹰等，2016)。改变围岩应力状态是指通过应力释放和转移，降低应力梯度或将应力集中区向深部推移，经过大量工程的检验(欧阳振华，2014；顾合龙等，2016；李俊平等，2018)，卸压爆破对围岩能量集中区进行能量释放来改变围岩应力分布，有效减缓围岩体应力集中程度，达到降低岩爆灾害的效果。

李鹏(2017)通过数值模拟得出了卸压爆破缓解高应力区围岩内部应力集中的计算结果，分析了卸压爆破后围岩应力转移现象；窦林名(2001)研究了围岩的岩石力学特性对应力分布的影响，以及通过卸压爆破改变围岩应力状态的方法；魏明尧等(2011)分析了不同卸压条件下围岩应力分布转移规律。吴世勇等(2011)通过FLAC3D分析了双排倾斜辐射孔和垂直超前孔等两种不同岩爆防治方案，验证了爆破卸压的重点部位；李俊平等(2015)应用FLAC3D静态模拟巷道掘进中巷帮钻孔卸压的钻孔参数的卸压效果，并采用ANSYS/LS-DYNA程序和FLAC3D程序相结合，动态模拟了硬岩巷道端面的钻孔爆破卸压过程，分析了不同卸压方案的钻孔爆破卸压效果。宋希贤等(2014)采用RFPA2D模拟和分析动载荷作用下巷道围岩动态损伤破坏发展以及围岩中应力场重新分布的过程，揭示动力扰动下深部巷道卸压孔与锚杆联合支护的作用机理。

综上所述，前人利用多种数值模拟方法实现了如下方面的研究，总结了较为丰富的成果：(1)卸压爆破前后应力转移情况研究；(2)卸压爆破效果的影响因素研究(如岩性、爆破位置等)；(3)不同卸压方案的对比分析，优化爆破参数，获得工程的最优爆破方案。

然而，前人多采用连续介质范畴的理论和方法进行连续-非连续的动态过程分析，且尚未进行隧道动态开挖过程的卸压爆破模拟的效果评价，采用基于FEM/DEM耦合的数值模拟方法CDEM研究冲击荷载作用下岩体破裂扩展-贯通的连续-非连续地质力学过程，以及卸压爆破对岩爆发育过程的抑制效应。

本文以巴陕高速公路米仓山隧道卸压爆破段为例，运用了基于连续-非连续CDEM方法的数值模拟软件，考虑了岩石的非均质性，施加对应的地应力，动(卸压爆破)、静(开挖)结合的方式模拟了隧道开挖、卸压爆破全过程。通过两组对比模拟，研究了卸压爆破抑制岩爆效应的机理，通过裂纹的萌生、扩展、交汇和贯通全过程，以及应力向径向深部转移的现象得出了未卸压条件下在强破损区岩爆可能性大。结合卸压爆破探讨从隧道横向纵向切应力峰值的减小以及损伤区域的变化以及损伤值的减小得出了卸压爆破有降低岩爆的效果，同时对比得出了卸压爆破对隧道轮廓非开挖岩体损伤区域的影响。

1 CDEM方法简介

本文采用基于CDEM方法(冯春等，2014)的基于GPU加速的块体动力学仿真软件(GDEM Block-Dyna)，引入应变强度分布模型实现动静荷载作用下岩石的损伤破裂行为的描述，引入Landau爆源模型实现卸压爆破过程的模拟。

图1为CDEM方法中块体的组成(冯春等，2010)。通过虚拟裂缝可表征多裂纹的萌生、扩展及交汇贯通过程。通过对罚弹簧上设置相适应的本构和断裂准则，可在虚拟界面上实现拉伸、剪切断裂过程；断裂发生后，虚拟界面即转化为真实的接触界面，通过赋予相应的接触模型及接触参数，即可对接触面的力学行为进行准确刻画。

图 1 CDEM的块体及界面Fig. 1 Block and interface of CDEM

图 2 米仓山隧道地形及地质纵断面图Fig. 2 Terrain and geological profile of Mountain Micang tunnel

1.1 单元边界本构模型

CDEM中岩石渐进破坏过程是通过单元边界的相继断裂实现的。本文在单元边界上采用的本构模型为应变强度分布模型(冯春等，2014)，该模型假设岩石代表性体积单元是由大量的细长纤维组成，岩石的破裂是内部强度不一的纤维随着拉伸应变的增大而相继出现断裂，从纤维开始断裂到全部断裂的过程，模拟了岩石的非均质性，刻画了岩石渐进破坏的过程。

1.2 爆源模型

本文中在模拟卸压爆破时，起爆点火方式采用到时起爆，爆生气体膨胀压力采用基于朗道-斯坦纽科维奇公式，在爆源模型中假定爆生气体的膨胀为等熵膨胀即爆生气体在膨胀过程中与围岩无热交换(陈保君等，2008)。

2 卸压爆破数值模拟

2.1 工程地质概况

米仓山特长公路隧道全长13.81ikm，为双线分离式隧道，左右线间距为40im，设计时为国内第二，世界第三长的公路隧道。隧址区处于秦岭东西向复杂构造带与龙门山华夏系构造的结合部，属于中深切割的侵蚀-构造剥蚀地貌区域，最大埋深约1000im，围岩以石英闪长岩为主，夹花岗岩、闪长岩等透镜体，属于典型的硬岩隧道，地质纵断面图如图 2所示。

卸压爆破区间段位于隧道中部(约K46+052)围岩主要为石英闪长岩，为Ⅲ级围岩，呈巨块状镶嵌结构或块状整体结构，基本无水，根据区段内测得最大水平主应力δ1为25.2iMPa，应力方向与洞轴线方向基本一致。石英闪长岩单轴饱和抗压强度为64.69iMPa，由相关《公路隧道设计规范》中地应力判别δ1=2.6<3，可知此区段为极高地应力区段。并且在现场施工过程中时常发生岩爆，图 3所示为岩爆引起初期支护开裂、变形、侵限甚至坍塌。岩爆严重威胁到施工人员的安全与施工进度以及永久性结构安全，故采取了卸压爆破的方式进行岩爆防治。

图 3 岩爆引起的初支开裂破坏Fig. 3 Initial cracking damage induced by rock burst

图 4 掌子面卸压孔布置示意图Fig. 4 Schematic diagram of the pressure relief hole arrangement on the face

2.2 岩爆段卸压爆破方案

为了改善掌子面附近围岩的应力集中情况，米仓山隧道右线桩号约K46+052采用了掌子面径向钻孔。掌子面卸压孔共6个，卸压孔与隧道轴线夹角为15°(图4)。卸压爆破单孔取装药量为2.2ikg，单根炸药长度为20icm，重量为0.1ikg，采用两端装药，底部连续装药，装药量为1.5ikg，长度为3im，其余部分为间隔装药，装药量为0.7ikg，长度为1.4im，堵孔长度为1.8im。具体卸压参数如表 1所示。

表 1 右线掌子面(K46+052)前方卸压参数Table1 Pressure relief parameter in front of the right line tunnel face(K46+052)

孔直径/mm炮孔深度/m装药量/kg孔间排距/m炮孔堵塞长度/m85102.21.1～1.51.8

表 2 岩石及衬砌力学参数表

Table2 Rock and lining mechanics parameter table

力学参数材料密度/kg·m-3弹性模量/Pa泊松比黏聚力/Pa抗拉强度/Pa内摩擦角/(°)剪胀角/(°)厚度/cm石英闪长岩28845.5e100.2672e71.2e755.015/衬砌24002.3e100.2////15

2.3 数值模型及方案设计

2.3.1 本构及相关参数选取

在数值模拟中，每个有限单元采用考虑应变软化效应的莫尔-库仑模型(Soften-MC)，单元之间的虚拟界面上采用应变强度分布模型，岩石的拉应变与剪应变取值服从Weibull概率分布，卸压爆破采用Landau爆源模型。

2.3.2 模拟方案设计

根据前期勘察资料，以及现场地质资料得出卸压爆破段围岩岩性单一，为花岗闪长岩，计算模型岩石力学参数见表 2。模型分组如图 5所示，尺寸为100im×100im×80im，含节点数为29302，总块体数为54900。模型中材料由围岩组、开挖组以及材料性质变化组3个部分组成，在开挖掘进时，对左右两侧边界施加水平位移约束，对底面边界施加垂直方向位移约束，模型中未能包含的上覆岩层所产生的压应力以体力的形式施加在每个有限单元上，在运用爆炸模块进行卸压爆破动态计算时，将模型边界条件更改为黏性边界条件以消除人工边界引起的虚假反射。

图 5 模型数值计算分组开挖示意图Fig. 5 Numerical simulation diagram of model

为了探讨现场卸压爆破之后岩爆减少的现象是否由卸压爆破造成，本文设计了两组数值模拟方案：A组为不进行卸压爆破的情况下的连续开挖；B组为先开挖一段围岩之后进行卸压爆破再继挖下一段围岩，这两组数值模拟所采用的所有参数均相同。因分组开挖单元在开挖过程中被移除，为减少不必要的计算量，故只对围岩部分开启虚拟界面计算，因此在损伤云图中未包含开挖组的损伤，以及损伤体积增长也未包含开挖组部分。

(1)A组数值模拟：模型共计17组，其中编号1～8组为待开挖岩体，为了减少开挖原因导致的结果差异，故每组开挖长度保持一致，9～16组为材料性质可变动的组，例如开挖第1组岩体后，与之对应的第11组材料从岩石的力学参数变为衬砌的力学参数。

(2)B组数值模拟：B组数值模拟选取的材料力学参数以及地应力参数均与A组相同，卸压爆破施加在开挖第1组与第2组围岩之间进行。卸压爆破区间位于2～3组开挖岩体之中。

2.3.3 初始地应力数据选取

地应力场作为影响岩爆的一个重要因素，本次数值模拟所选取的原岩应力场数据为K46+725断面用应力解除法测定地应力数据。该测点与卸压爆破区段相距约700im，高程差仅有约20～50im，且区间内无断层褶皱等地质构造现象，参考贾金晓(2016)米仓山隧道隧址区初始地应力场反演结果，卸压爆破区间段地应力与K46+725处差别很小，故选取此点地应力数据作为数值计算所使用的地应力值。将所测定的该区段的地应力测试结果(表3)转换为与模型对应的坐标系，具体应力分量如表 4所示。

表 3 地应力测试结果表Table3 Test result of ground stress

测点埋深/m主应力参数主应力大小/MPa方位角/(°)倾角/(°)K46+725550δ125.2251-1.5δ213.1-21.666.8δ312.9149-75.6

表 4 应力分量计算结果

Table4 Calculation result of the stress component

δx/MPaδy/MPaδz/MPaδxy/MPaδyz/MPaδxz/MPa25.11313.64.20.240.45

3 数值模拟结果分析

3.1 未卸压条件的数值模拟结果分析

天然状态下岩石自身存在着损伤，在隧道开挖所诱发的压应力下，围岩体内局部形成了损伤拉应力状态。局部微破裂将会随着开挖卸荷导致的围岩应力增加而扩展，从而在隧道轮廓的某些区域裂纹聚集和贯通，随着进一步的开挖掘进，这些严重损伤区域在高地应力作用下就可能成为岩爆的高发区。

伴随着开挖掘进，应力重分布所导致的应力集中现象，当含有接触面的区域某个块体与目标块体相互连接的“半弹簧”超过其特征长度时，弹簧失效则代表该块体与目标块体发生分离。

图 6为隧道开挖未卸压模型中轴线纵剖面损伤分布云图，由于篇幅限制，文中后续出现的横剖面图均在强破损区选取的典型剖面。图 7为未卸压条件下强损伤区围岩损伤破坏裂纹扩展变化图，在隧道附近围岩相继出现破损，最终在拱顶偏右的部位形成严重破损区。

图 6 未卸压开挖强损伤区典型纵剖面损伤分布云图Fig. 6 Nephograms of longitudinal section damage distribution of unrelieved excavation model

在开挖过程中裂纹从萌生到扩展，到最后的局部裂纹区域贯通形成严重破损区，开挖掘进以及应力集中导致了部分岩石内部破裂，强度降低。图 8为强破损区典型剖面最大主应力云图，破损区应力的释放导致应力峰值向径向一定深度方向发生转移，应力峰值达到了37.4iMPa，围岩处于高围压状态。高围压卸荷比低围压卸荷有更高的张拉断裂比重，岩石体现出更明显的脆性破裂性质(马春驰等，2015)，该部分完整性较好的岩体在短时间内蓄积了大量的弹性应变能，当这部分弹性应变能大于块体单元表面能时，右边墙靠近拱顶强破损区的岩石可能在这部分能量的作用下发生岩爆。其余部位处围岩虽然出现了破损，产生了少量裂纹尚未贯通形成破碎区，围岩完整性稳定性较高。

图 7 未卸压条件下裂纹萌生扩展变化图Fig. 7 Crack initiation and propagation diagram under unreduced pressurea. 裂纹萌生；b. 裂纹扩展；c. 裂纹贯通形成严重破损区

图 8 未卸压开挖强损伤区典型剖面最大主应力云图Fig. 8 Maximum principal stress nephograms of typical section of strong damage zone

3.2 卸压条件的数值模拟结果及其对比分析

从数值模拟结果来看，在进行卸压之后的隧道开挖模型发生岩爆的可能性大幅降低，达到抑制岩爆的作用。下文通过损伤变化、裂纹扩展规律、地应力变化3个方面进行对比分析。

3.2.1 卸压对围岩损伤的影响

在爆破应力波和爆生气体共同作用下势必会对围岩产生损伤或者破坏(图 9)，将开挖掘进过程中损伤因子大于0.1的块体计入统计，围岩损伤率Dr=VD/V，式中，VD表示损伤因子大于0.1的块体的体积，V表示模型总体积(未包含已开挖部分)。如图 9所示，在未卸压开挖第1组时，围岩开始出现破损，损伤率Dr为0.025%，掘进至20im时，损伤率增加到0.26%。从损伤面积折线图体现出了围岩损伤的萌生到扩展再到损伤区域贯通的过程，而在B组卸压条件下，由于卸压爆破的影响，围岩损伤体积明显高于A组未卸压条件下。从曲线斜率可以看出卸压爆破虽然导致了损伤区域短时间内快速增长，但损伤区域增长率出现递减的趋势，而在未卸压条件下，每步开挖都会导致损伤区域的增长且由于应力的集中程度较高，围岩损伤破坏率在逐步提高的过程，并且最终可能导致岩爆的形成。

图 9 损伤区域体积占比图Fig. 9 Ratio of damage zone volume

图 10 强损伤区典型损伤云图对比图(未卸压/卸压)Fig. 10 Nephograms of damage comparisona. 未卸压条件下典型剖面损伤云图；b. 卸压条件下典型剖面损伤云图

但如图 10所示，对比同一截面的损伤值以及损伤区域可以看出，由于在爆破荷载作用下围岩出现全方位的扰动，区别于应力集中条件下的局部区域扰动，卸压条件下围岩破损区域虽明显大于未卸压条件下的围岩损伤(影响范围约隧道洞径一倍范围)，但未出现损伤因子为1的区域，受损围岩仍然具有一定的强度。从图10b中可以看出在拱顶以及拱底虽然出现大范围的破损，但损伤因子范围在0.1～0.5内，且由于“半弹簧”以及应变强度分布模型的原因，在其部分区域的目标弹簧的应变并未达到产生破损所需要的应变值，所以该部分区域在计算结果云图中并未显示出现损伤。

3.2.2 卸压对围岩裂纹扩展的影响

在卸压条件下，由于爆破对围岩的扰动，导致应力重分布，峰值应力向后推移，围岩破损起裂情况也较未卸压条件下出现了不同。卸压爆破后掘进过程中典型剖面裂纹的发展情况如图 11所示，对比未卸压条件下裂纹扩展过程，裂纹反应了破损情况，从破损规模上卸压条件下的开挖破损区域大于未卸压条件下的开挖，但从裂纹扩展贯通上的角度来看，卸压条件下由于爆破的扰动导致围岩出现破损以及应力的调整，此时虽有裂纹的产生，并且损伤半径比值(卸压/未卸压)分别为3.4倍(拱顶)、1.47倍(拱底)，损伤体积比值(卸压/未卸压)为5.36倍，但是并未能如未卸压开挖条件下的裂纹贯通形成强破损区，故降低了岩爆发生的可能性。

图 11 卸压条件下裂纹扩展变化图Fig. 11 Crack propagation change under pressure reliefa. 卸压条件裂纹萌生；b. 卸压条件裂纹扩展最终形态

3.2.3 卸压对地应力的影响

图 12 卸压孔口前方最大主应力峰值图Fig. 12 Peak diagram of maximum principal stress before pressure relief hole

由图 12可以看出，在卸压孔口前方卸压与未卸压两种条件下的最大主应力分布规律，未进行卸压的情况下，最大主应力为39.7iMPa，为原岩初始最大主应力的1.6倍，峰值区处于强破损区，最大应力梯度为4.8iMPa·m-1，此时应力集中程度高，岩爆危险性大；进行卸压之后，应力峰值明显降低，应力峰值区向后推移了约4im，最大应力梯度降低为4.0iMPa·m-1，应力集中程度明显降低。相对于未卸压条件下，爆破对围岩体提前进行了损伤破坏，且损伤范围较之未卸压条件下范围更广，降低了围岩体积聚的弹性应变能，对前方的围岩体也起到了松弛作用，从而缓解了岩爆产生所需要的强度和能量条件，卸压效果明显。

4 微震监测及卸压效果检验

微震监测技术是基于弹性波的解译来分析岩体内部微裂纹扩展以及岩体稳定性的监测方法(马春驰，2017)。已有的研究成果表明(谢鹏飞等，2012；乔小龙，2017；马春驰等，2018)，微震监测可以接收围岩破损时产生的声发射信号，并定位其声源位置。

4.1 微震监测结果对比分析

米仓山隧道微震监测采用的是加拿大ESG微震监测系统。通过微震监测接受的声发射数据的变化来反映防治岩爆措施的实施效果。张航(2015)采用该套设备对汶川县303省道紫荆隧道进行微震监测定位的精度验证，实验结果表明定位误差约13.5im。

为了验证数值模拟结果的正确性，选取未卸压区段K46+120—K46+80区段，与卸压区段K46+60—K46+30两者对比验证。

4.1.1 未卸压区段微震监测结果分析

K46+120—K46+90区间段施工的每日微震事件监测数据汇总如图 13所示，围岩的破裂伴随着微震信号产生。对应数值模拟分析的损伤的变化可以看出，在掘进开始时，裂纹开始萌生，随着开挖的进行，损伤区域的快速增加，微震后续监测中的微震事件数快速增加。

图 13 未卸压区段微震事件统计图Fig. 13 Statistics of microseismic events in unrelieved sections

图 14 12月29日岩爆破坏示意图Fig. 14 Tunnel destroy by rock burst

图 15为2017年12月22～29日微震监测微震事件分布图，由图可以看出，微震事件主要分布在隧道右边墙和拱底，事件总数为108次，最大地震矩1.2，地震能量达1+e7焦耳，释放能量加大，对应未卸压组数值模拟结果中岩体部分区域所产生的强破损区，由于微震监测的误差以及数值模拟的局限性，强破损区的位置没有呈现较好的一致性，但如图 14所示，经历了破裂损伤的快速扩展后，围岩损伤率降低，最终由于强破损区的形成以及高地应力的作用下在2017年12月29日发生了岩爆灾害，且岩爆位置与数值模拟结果具有较好的对应性。结合微震监测数据以及现场岩爆实际位置可以验证未卸压组数值模拟结果的正确性。

图 15 未卸压区段岩爆微震事件分布及震源参数Fig. 15 Distribution and source parameters of rockburst microseismic events in unrelieved sections

4.1.2 卸压区段微震监测结果分析

图 16为K46+60—K46+30区间段施工的每日微震事件监测数据汇总，在进行了卸压爆破的后续2id内，微震事件数较多，在卸压爆破的作用下岩体的损伤破裂较多。从图 17微震事件分布图对比数值计算结果，可以看出卸压爆破导致围岩破损范围较未卸压条件增大，但最大地震矩降低为0.529，结合微震事件统计和数值模拟结果可以看出，后续开挖岩体损伤率有逐渐降低的趋势，虽然围岩的破损范围更广，但有利于能量提前释放，减缓了应力集中程度。数值模拟中连续开挖岩体破裂情况较好地吻合了微震监测结果，从而验证了数值模拟结果的正确性。

图 16 卸压区段微震事件监测统计图Fig. 16 Statistics of microseismic events in relieved sections

图 17 卸压区段岩爆微震事件分布及震源参数Fig. 17 Distribution and source parameters of rockburst microseismic events in unloading section

微震监测情况同数值模拟结果综合分析可以看出，卸压爆破通过提前对隧道轮廓周围围岩的扰动，导致围岩破损，能量提前释放，使得应力重分布，降低应力峰值的方法降低了岩爆发生的可能性。从现场施工情况来看，卸压区段未发生岩爆，卸压效果明显。

4.2 小 结

微震监测结果与两组数值模拟对比中卸压爆破段与未卸压爆破段破损区域的范围变化相一致；以及现场岩爆实例与数值模拟未卸压段强破损区位置相对应。两者结合验证了数值模拟结果的正确性。

在岩爆段施工过程中可采用CDEM法分别模拟两种开挖工况，模拟结果结合微震监测的综合分析，判断岩爆灾害的位置和等级，结合施工经验，选取钻孔卸压或者爆破卸压的方式从而降低岩爆灾害的发生。

5 结果与结论

本文运用连续-非连续的计算方法进行了隧道连续开挖诱发岩爆的数值模拟，并将卸压爆破作为唯一变量进行了两组模拟，结合现场微震监测，对卸压爆破抑制岩爆效应机理进行了分析，并得出以下结论：

(1)卸压爆破作为一种“以爆治爆”的岩爆防治方法，通过爆破产生的应力波破碎一定范围的岩体进而衰减成弹性波扰动围岩，使围岩蓄积的弹性能得到释放，避免了围岩弹性能的继续蓄积，能量的转移降低了岩爆发生的风险。

(2)CDEM适合围岩卸压爆破的连续-非连续全过程分析，为其他工程采用卸压爆破抑制岩爆效果和分析爆破对围岩体影响程度提供参考，可为实际工程中在岩爆段采取卸压爆破后围岩支护方案提供一些的参考价值；通过模拟结合微震监测为岩爆段施工提供是否卸压以及卸压强度的选择。

(3)围岩连续开挖过程中伴随着裂纹的起裂、裂纹的扩展进而在某些部位产生贯通的裂纹即形成局部强破损区域，这些局部的强破损区在高地应力和应力梯度较大的条件下易导致岩爆灾害，得出了在未卸压条件下持续开挖诱发岩爆灾害的可能性较大。

(4)对比分析了隧道围岩地应力场、损伤值以及损伤区域变化规律，由于卸压爆破对围岩的扰动，并使其损伤，虽损伤区域较未卸压分别扩大2.4倍(拱底)、0.47倍(拱底)，但因提前释放了隧道轮廓附近大部分围岩储存的弹性能，使得围岩整体损伤因子处于0.1～0.5范围内，同时也导致了应力峰值向后转移了4im。

(5)分析了卸压爆破对隧道轮廓附近围岩稳定性的影响，由于围岩整体损伤因子处于一个较低的水平，裂纹未能交汇贯通，围岩仍具有一定的稳定性。