多次爆破作用下巷道围岩累积损伤演化与破裂特征

黄江，陈士海，2，曾凡福，苏松

(1.华侨大学 土木工程学院，福建 厦门 361021；2.华侨大学 福建省隧道与城市地下空间工程技术研究中心，福建 厦门 361021；3.中铁一局集团有限公司，陕西 西安 710054；4.中铁十四局集团有限公司，山东 济南 250014)

地球浅部资源日益枯竭[1]，我国将深入地下1 000～2 000 m开采矿产资源，甚至延伸到更深的地方[2].钻爆法因经济、高效和适应性强等优点已经成为开采矿产资源的主要方式.使用钻爆法开挖地下巷道时，预测爆破引起围岩的损伤区域对岩体长期稳定具有重要的意义[3-4].深部岩体处于复杂的力学环境，较大的地应力是影响爆破开挖的重要因素[5].岩体中的地应力对围岩损伤的演化和爆破裂纹的发展都有很大影响[6-8].在爆破掘进过程中，爆破冲击波是瞬间产生并作用于岩体，深部岩体的原位应力受到爆破荷载的影响，可能导致一些难以预料的灾害[9].岩体的初始应力状态对地下深部工程结构稳定性评估起着重要作用[10].爆破扰动与地应力耦合作用下岩体受损，其力学特性也受到影响.特别是多次爆破扰动作用下岩体物理力学性质不断劣化，围岩损伤程度逐渐累积，岩体损伤区域也在不断地扩大[11-12].因此，研究多次爆破荷载作用下深部巷道围岩损伤演化与破裂特征具有重要意义.

深部巷道爆破开挖过程中，围岩累积损伤演化过程十分复杂，在深部工程进行现场爆破试验比较困难.所以，数值模拟成为了分析和评估岩石力学中复杂问题的一种有效方式[13-15].陈俊桦等[16]建立了考虑初始损伤的弹塑性损伤本构模型，将数值模拟结果与工程实例结合，提出评价围岩受爆破影响的损伤判据.李新平等[17]把损伤本构嵌入FLAC3D软件中，通过数值模拟方法预测了地下洞室爆破开挖损伤范围.刘闽龙等[18]建立了各向异性动态损伤本构模型并将其用于隧道爆破开挖的数值模拟中，研究了单孔爆破的损伤范围.Yi等[19]通过理论模型和数值模拟技术研究了地应力条件下单孔爆破裂纹发展趋势.数值模拟技术已经成为岩石力学研究中的利器.因此，本文对多次爆破作用下巷道围岩累积损伤演化与破裂特征进行研究.

1 考虑扰动效应的时效损伤本构模型

1.1 一维时效损伤本构模型

在多次爆破开挖后，深部岩体工程可能发生渐进性变形.时效性损伤本构模型能有效描述非线性时效变形特征.经典的一维Bingham流变模型是由线性弹性元件和理想粘塑性元件串联而成，能较好地描述岩石的稳态蠕变阶段，但是不能描述减速蠕变阶段和加速蠕变阶段.深部巷道爆破掘进是一个动态扰动过程，一般在加速蠕变阶段前，围岩就已经破坏.

经典的一维Bingham流变模型经过改进后，可描述岩石减速蠕变阶段.采用非线性函数F(t)描述岩石减速蠕变阶段[20]，其表达式为

F(t)=1-e-at.

(1)

式(1)中：a为系数；t为时间；e为自然常数.

在常应力作用下，非线性弹性原件的弹性模量函数E(t)[21]的表达式为

图1 扰动荷载与无扰动荷载作用下岩石应变曲线Fig.1 Rock strain curves under disturbance load effect and non-disturbance load

.

(2)

式(2)中：E0为非线性弹性原件的初始弹性模量.

采用非线性函数，按式(2)改进经典的一维Bingham流变模型中的弹性元件，得到一个描述岩石减速蠕变阶段的模型.蠕变方程为

(3)

式(3)中：ε为应变；σ为应力；σs为屈服应力；η1为粘滞系数.

假设扰动荷载作用产生应变增量，作用期间忽略蠕变变形.扰动荷载与无扰动荷载作用下岩石应变曲线，如图1所示.由图1可知：在扰动荷载作用期间[22]，岩石产生了较大的无法恢复的变形，扰动荷载作用引起的应变增量远大于岩石的时效变形.这意味着岩石内部衍生了微裂纹，存在损伤.

为了描述动态扰动荷载作用过程中的岩石损伤程度，建立一个扰动损伤元件，该元件只在扰动荷载作用时间td内起作用，产生的应变εR即为扰动荷载作用引起的应变增量.扰动荷载每作用一次，岩石内部就会产生相应的损伤，损伤元件本构模型方程为

σ=EdεR(1-D).

(4)

式(4)中：Ed为扰动损伤元件的弹性模量；D为岩石损伤.

不考虑单次扰动施加过程中的损伤变化，损伤变量的微分形式为

(5)

设ε0是第一次扰动开始加载时的等效应变，此时，D=0；εd是最后一次扰动加载结束后的等效应变，此时，D=1.D的表达式为

D=[(εR-ε0)/(εd-ε0)](εd/εR).

(6)

考虑到爆破损伤只会在爆破过程中产生，因此引入算子K(σd)，保证损伤元件只在爆破荷载扰动过程中产生作用.算子K(σd)表达式为

(7)

式(7)中：σd为爆破荷载.

每爆破一次都会产生一个应变增量，记第n次爆破产生的应变增量为εR(n)，引入一个与爆破次数有关的衰减函数F(n)，其表达式为

F(n)=e-bn.

(8)

式(8)中：n为扰动作用次数；b为拟合系数.

损伤元件的本构模型方程为

(9)

考虑扰动效应的时效损伤元件模型，如图2所示.考虑扰动效应的一维时效损伤本构模型方程为

图2 考虑扰动效应的时效损伤元件模型Fig.2 Time-dependent damage element model considering disturbance effect

(10)

1.2 三维时效损伤本构模型

岩石多处于三维的应力状态，一维时效损伤本构模型难以描述岩石在多向应力状态下的蠕变行为，因此，需将式(10)拓展为三维情形.假设岩石为各向同性体，根据广义Hooke定律，弹性元件三维本构关系为

(11)

式(11)中：Si，j为偏应力；G为剪切模量；ei，j为偏应变；σm为球应力；K为体积模量；εm为球应变.

分解岩石内部的应力，得到非线性弹性元件A的应变方程为

(12)

在三维应力下，岩石体积蠕变近乎为0，三维蠕变方程中不考虑体积蠕变，仅考虑瞬时体积变形，并引入Drucker-Prager准则描述岩石的屈服，采用相关联的流动法则，即塑性势函数(Q)与屈服函数(F)关系为Q=F.则F的表达式为

(13)

式(13)中：I1为应力的第1不变量；J2为偏应力的第2不变量；ρ，k为Drucker-Prager准则材料常数.

粘塑性体B的三维本构关系为

(14)

设扰动损伤元件C中体积变化是弹性的，流变性质主要表现在剪切变形方面，即

(15)

三维形式下的损伤(Dt)表达式为

Dt=[(εRS-ε0S)/(εdS-ε0S)](εdS/εRS).

(16)

式(16)中：εRS，ε0S和εdS均为三维形式下的等效应变.

考虑扰动效应的三维时效损伤本构模型表达式为

(17)

式(17)中：a为衰减系数，t为时间.

(a) 岩石试样 (b) 数值模型图3 岩石试样与数值模型Fig.3 Rock samples and numerical model

采用LS-DYNA软件提供材料本构模型，二次开发材料子程序接口User Defined Material.采用Visual Studio 2013软件，通过Fortran语言编写材料本构模型源程序.采用Parallel Studio XE 2010软件，将自定义本构模型代码与材料本构模型源程序联合编译成动态链接库文件.为了验证三维时效损伤本构模型的合理性和正确性，在ANSYS/LS-DYNA软件中进行仿真模拟.数值模型以实际岩样规则为准，岩石试样的直径为50 mm，高为100 mm.岩石试样与数值模型，如图3所示.图3中：三维时效损伤本构模型共有3 146个节点和2 500个单元.

根据岩石三轴动态扰动多级加载实验的结果[22]，得到红砂岩蠕变损伤本构模型的参数.考虑扰动效应的时效损伤模型参数如下：非线性弹性元件A的初始弹模为8.85 GPa，泊松比为0.2；粘弹性元件B的弹性元件衰减系数为35.8，粘弹性元件粘滞系数为150；损伤元件C的损伤原件剪切模量为10.8 GPa，损伤原件衰减系数为1.1，初始损伤应变与破坏时损伤应变差值为0.025.

在有限元模型底部施加Z方向的约束，在有限元模型顶部施加轴向均布荷载43.41 MPa，在侧面施加与三轴压缩蠕变试验相同的均布径向荷载0 .待试样变形趋于平稳后，开始施加动态扰动荷载，扰动荷载的幅值为10 MPa，应力条件与实际三轴压缩蠕变试验保持一致.

岩石试样加载模型，如图4所示.图4中：σw为围压应力；σax为轴压应力.轴向应变数值模拟与试验结果对比，如图5所示.

图4 岩石试样加载模型 图5 轴向应变数值模拟与试验结果对比Fig.4 Loading model of rock sample Fig.5 Comparison between axial strain numerical simulation and test results

由图5可知：轴向应变数值模拟曲线与试验应变结果曲线较为吻合；蠕变模型参数具有正确性.因此，考虑扰动效应的时效损伤本构模型在LS-DYNA软件中二次开发成功.

2 数值模拟

2.1 数值模型的建立

某地区矿山巷道岩层主要以红砂岩为主，致密坚硬，岩层性质较稳定.矿山巷道为圆形断面，其半径为2.5 m.大多数地下工程都涉及无限域或半无限域，而有限元法处理这类问题通常是在有限区域里进行离散.为了使有限元法处理不产生过大的误差，离散区域必须有足够的范围，并使区域外边界条件尽可能接近实际状态.考虑巷道爆破影响范围及“边界效应”的影响，建立的有限元模型尺寸(长×宽×高)为40 m×40 m×50 m，有限元模型采用3D Solid 164实体单元进行网格划分.

图6 数值模型Fig.6 Numerical model

由于距离爆源越近的岩体受到的冲击作用越强烈，所以爆破开挖处的网格加密的模型单元数为162 475，节点数为172 096，最小单元网格尺寸为0.03 m，最大的单元网格尺寸为3.76 m.圆形巷道处于模型中心位置，数值模型，如图6所示.图6中：为了模拟岩体对爆破冲击波的削弱作用，在模型周围施加无反射边界条件.

2.2 地应力工况

地应力的存在对巷道爆破开挖有很大影响，所以当深部岩体进行多次爆破模拟时，应考虑岩石所受竖直地应力的自重应力(σv)及水平地应力的构造应力(σh)[23]，构造应力与自重应力比值大多为0.8～1.2[24]，因此，深部岩体的侧压力系数λ=1.未来10 a，中国矿产资源多数在地下2 000 m进行开采，地应力大约为50 MPa[25].据此拟定了4种不同地应力爆破工况:工况1～4，其地应力分别为0，15，30，45 MPa.

(a) 地应力 (b)爆破荷载图7 地应力与爆破荷载的施加方式Fig.7 Application method of in-situ stress and blasting load

2.3 爆破荷载的施加

为探明在多次爆破冲击作用下，不同地应力岩体累积损伤演化过程与爆破损伤区域的扩展规律，首先，对模型施加地应力以形成初始应力场；然后，在巷道轮廓上施加爆破荷载.地应力与爆破荷载的施加方式，如图7所示.

图8 爆破荷载示意图Fig.8 Schematic diagram of blasting load

采用LS-DYNA软件计算爆破荷载，爆破荷载形式、荷载峰值和作用时间是计算的关键.在爆破开挖时，爆破荷载对巷道轮廓的压应力急剧上升，达到峰值后又立刻衰减，整个爆破过程都非常迅速、剧烈.在复杂的巷道爆破工程中，时程曲线一般采用半经验半理论的计算公式.

为了方便使用，将Sharpe-Duvall爆破理论中的双指数函数型爆破冲击荷载进一步简化为三角形荷载，其具有线性上升段和线性下降段，因此具有良好的适用性[26].数值模拟的加载时间为0.05 s，卸载时间为0.20 s，垂直施加于巷道轮廓面的等效爆破荷载峰值为137 MPa.爆破荷载示意图，如图8所示.

2.4 数值计算结果及分析

通过重启动技术实现多次爆破冲击波的施加，采用关键字“*STRESS_INITIALIZATION”实现应力、应变、损伤等参数的继承，计算5次爆破荷载作用下深部巷道围岩累积损伤效应.

2.4.1 围岩累积损伤效应 根据SL 47-2020《水工建筑物岩石地基开挖施工技术规范》[27]，当岩石损伤达到0.19，即判定围岩受到爆破损伤.假设开挖时未爆破，围岩无损，即初始损伤(D0)为0.采用LS-PrePost软件提取多次爆破后围岩的损伤云图(图9)，并分析其损伤演化过程.

(a) 工况2第1次爆破 (b) 工况2第2次爆破 (c) 工况2第5次爆破

(d) 工况4第1次爆破 (e) 工况4第2次爆破 (f) 工况4第5次爆破图9 不同工况下爆破损伤云图Fig.9 Blasting damage cloud diagrams under different working conditions

由图9可知：在地应力与巷道爆破荷载耦合作用下，围岩产生了一定范围的损伤；同一地应力状态下，巷道爆破损伤范围随着爆破次数的增加而扩大；围岩损伤程度随爆破次数增加而增加；距离巷道轮廓越近，围岩损伤程度越大，这是因为初始的爆炸冲击波直接作用于巷道围岩导致岩体力学性质急速劣化，岩体内部裂纹快速萌生与发展，爆炸冲击波对中远区岩体作用减小.

为了探讨多次爆破荷载作用下爆破损伤范围变化特征，分别提取了每次爆破后围岩累积损伤值为0.19的位置，并测定了与巷道轮廓的距离.不同地应力下损伤范围，如图10所示.图10中：s为爆破损伤.由图10可知：在较大地应力、相同爆破次数下，巷道爆破损伤范围明显小于较小地应力的损伤范围，爆破损伤范围随着地应力增加而减小.

当侧压力系数λ=1时，在圆形巷道形成对称分布的环向压应力场，环向压应力场与爆破应力波方向相反，在爆破过程中压应力场消耗大量爆破能量，削弱了爆破应力波对围岩的冲击作用.由于形成微裂纹的能量大幅度减少，所以裂隙的扩展被抑制.地应力越大，对爆破荷载的削弱作用越强，围岩内部裂纹长度与数量也会相应减少，地应力对围岩损伤演化起到抑制作用，因此，爆破损伤范围也会减小.

不仅较大地应力的爆破损伤范围小于较小地应力的爆破损伤范围，而且每次爆破损伤范围的增量也更小.当地应力为0时，在相同爆破次数下，爆破损伤范围的增量明显大于地应力为15，30，45 MPa的爆破损伤范围增量；当地应力为30，45 MPa时，第2次爆破后损伤范围为0，其原因是地应力的抑制作用使围岩的损伤程度未达到0.19.由此可见，爆破损伤随着地应力的增加呈非线性减小.

为了研究不同地应力下的巷道围岩累积损伤变化特征，在数值模拟结果中提取了巷道轮廓顶部单元损伤值.不同地应力下爆破累积损伤结果，如图11所示.由图11可知：围岩损伤程度随着爆破次数的增加而增加，呈现一种阶梯式增加趋势，但是每级阶梯的高度(每次爆破损伤增量)随着爆破次数的增加而降低，表明在多次爆破荷载作用下，巷道围岩损伤具有非线性累积特性；由于地应力抑制岩体内部微裂纹的萌生、扩展与贯通，围岩损伤程度随着地应力的增加而减小，特别是在多次爆破荷载作用下，地应力的抑制作用越明显.

图10 不同地应力下损伤范围 图11 不同地应力下爆破累积损伤结果 Fig.10 Damage zone under different in-situ stresses Fig.11 Cumulative damage results of blasting under different in-situ stresses

爆破损伤增量，如图12所示.图12中：ΔD为岩石损伤增量.由图12可知：每次爆破产生的损伤增量并不一样，损伤增量随爆破次数n的增加而减小，与爆破扰动次数呈现出明显的非线性关系.抑制系数(ξ)与地应力的关系，如图13所示.

图12 爆破损伤增量 图13 抑制系数与地应力的关系 Fig.12 Blasting damage increment Fig.13 Relationship between restraining coefficient and in-situ stress

单次爆破损伤增量随着爆破扰动次数n呈指数衰减形式，引入地应力抑制系数ξσ，其表示地应力对爆破损伤的影响.损伤增量(ΔD)表达式为

ΔD(σ)=(1-ξσ)αe-βn

(18)

式(18)中：α，β为拟合参数.当地应力为0时，地应力对损伤无影响，ξσ=0；当地应力为无穷大时，地应力抑制系数ξσ=1，表示完全抑制爆破损伤.地应力抑制系数ξσ表达式为

ξσ=1-[(e-μσ+ν)/σ].

(19)

式(19)中：μ，ν为拟合参数.

利用MATLAB软件进行全局拟合，得到损伤增量函数为

ΔD(σ)=[28.604 1(1-e-0.034 96σ)/σ]×0.360 87e-0.467 78n.

(20)

此时，地应力抑制系数ξσ表达式为

ξσ=1-[28.604 1(1-e-0.034 96σ)]/σ.

(21)

2.4.2 围岩破裂分析 为了研究爆破过程中的巷道破裂特征，进行模型试验与数值模拟.通过地下工程三维模型试验系统模拟爆破开挖过程，地下工程三维模型的尺寸(长×宽×高)为1 500 mm×2 000 mm×1 500 mm；巷道半径为100 mm；填筑材料为砂浆混凝土；材料质量(m)比为m(细砂)∶m(硅酸盐水泥)∶m(石膏)∶m(重晶石粉)∶m(水)=8.0∶0.2∶0.5∶0.3∶1.0(m为质量)[28].在水中加入质量浓度为1%的缓凝剂硼砂.地下工程三维模型，如图14所示.

爆破进尺为0.1 m，在隧道掌子面预埋炮孔，采用非炸药柱状药卷模拟炸药，实施电火花爆破.采用LS-DYNA软件建立数值模型，在K文件中添加关键字“*MAT_ADD_EROSION”，模拟巷道爆破开挖过程中岩体破裂过程.电火花爆破装置，如图15所示.

图14 地下工程三维模型 图15 电火花爆破装置 Fig.14 3D model of underground engineering Fig.15 Electric spark blasting device

(a) 模型试验 (b) 数值模拟图16 模型试验与数值模拟Fig.16 Model test and numerical simulation

模型试验与数值模拟，如图16所示.由图16可知以下4点结论.

1) 在多次爆破荷载作用下，巷道顶部与两侧有明显岩体破裂痕迹，顶部有部分岩块剥落，岩体裂纹形成阶段主要有微观裂纹阶段和宏观裂纹阶段.

2) 在多次爆破荷载作用下，巷道围岩发生变形，岩体中的微小裂隙随机形成，彼此独立.微小裂隙对岩体力学性能影响微小，岩体损伤处于较低水平，因此，对岩体承载能力影响较小.

3) 随着围岩持续变形及围岩损伤程度的增加，裂纹数量和扩展速率都增加，岩体损伤累积效应显现，从微观破坏向宏观破坏过渡，裂纹的产生和演化从随机变为有序，巷道围岩上出现数条较大的主裂纹，围岩表面的开裂与剥落对巷道围岩力学性能的影响也在逐渐增大.

4) 围岩表面因受拉破坏而破裂[29-30].爆破荷载垂直作用于巷道壁上，爆破应力波是一个快速加载、卸载的过程.在多次爆破荷载作用下，围岩强烈地压缩，巷道周围岩体积蓄了部分弹性变形能.在卸载阶段，压力迅速降低，达到一定程度时，岩石内部积蓄的弹性变形能快速释放，转变为卸载波，爆破荷载对岩体的作用由压缩转变为拉伸，当拉伸超过一定的强度时，岩石拉裂形成裂纹[31].

模型试验与数值模拟结果共同揭示了在多次爆破荷载作用下的围岩破裂演化机制，这将有助于深部巷道爆破开挖施工安全理论的发展.

3 结论

建立了考虑扰动效应的时效损伤本构模型，在LS-DYNA软件中实现本构模型的二次开发，通过数值模拟研究了不同地应力下巷道围岩爆破损伤演化特征与围岩破裂特征，得到了以下3个结论.

1) 在地应力与多次爆破扰动荷载耦合作用下，距离巷道轮廓越近，围岩损伤程度越大；爆破损伤范围随着地应力的增加呈非线性减小；地应力越大，抑制作用越明显.

2) 在多次爆破荷载作用下，巷道围岩损伤具有非线性累积特性.随着爆破作用次数的增加，围岩损伤阶梯式增加，但每次爆破损伤增量减小，地应力对巷道围岩损伤的抑制作用更加明显.

3) 在多次爆破荷载作用下，巷道顶部与两侧有明显岩体破裂痕迹，顶部有部分岩块剥落，巷道围岩上出现数条较大的主裂纹，围岩表面的开裂与剥落对巷道围岩力学性能的影响也在逐渐增大.