随着地表浅层矿产资源的逐渐枯竭，通过深部开采获取资源将逐渐常态化^[1-2]。与传统浅部开采相比，深部开采需要更加先进的技术和设备，从而保证资源的有效利用与合理开发。其中，钻爆法仍是深部硬岩开挖的主要手段，广泛应用于巷道开挖、井筒掘进和采场爆破等工程领域。浅部开采所受地应力较小，地压对岩体爆破碎裂的影响可忽略不计，而深部开采所受地应力较大，千米深井实测地应力结果可达30～50 MPa^[3]，深部岩体在高地应力下爆破时会表现出与浅部不同的力学行为。

关于地应力下岩体爆破，国内外学者已开展诸多试验、理论及数值计算方面的研究。其中，最早的研究可追溯至1966年由Nicholls等^[4]开展的地应力下预裂爆破试验，首次提出爆破径向裂纹更倾向于沿最大主应力方向扩展。基于该试验，Kutter等^[5]通过不同材料小型爆破试验发现单向压力对爆破裂纹的萌生和扩展具有导向作用，从而印证了Nicholls等的结论。国内方面，肖正学等^[6]基于前人的研究，结合室内与现场爆破试验数据，发现了类似的爆破裂纹扩展行为规律，并指出初始应力场会改变爆轰波的传播规律。此后，刘殿书等^[7]基于动光弹实验记录的动态等差条纹图，分析了不同初始应力下爆破应力波的变化过程，验证了肖正学等^[6]的观点，并将该结论进一步具体化。在爆破试验测量仪器方面，杨仁树等^[8-9]将动态焦散线方法与高速摄影技术相结合，建立了一套适用于爆炸、侵彻等超动态问题研究的数字激光动态焦散线试验系统，通过该系统解决了岩体爆破领域难以实时观测爆破过程的一些问题。岳中文等^[10]通过该系统分析了单向初始压力作用下切缝药包爆破裂纹的扩展行为，指出切缝药包可以有效控制爆炸能量分布，使能量沿切缝方向集中释放，进而实现控制爆破中材料的定向断裂。随后，Yang等^[11]进一步开展了切缝药包控制爆破相关研究，基于弹性与断裂力学理论，分析了初始应力对切缝药包爆破裂纹扩展行为的影响，讨论了不同切缝角度下的爆破损伤特征，揭示了深部岩体切缝药包爆破的破裂机理。

此外，数值模拟以其简便/普适性，被广泛应用于岩体爆破过程研究。常用的数值方法主要包括有限元、离散元两类，其中使用最为广泛的软件当属LS-DYNA。Ma等^[12]基于该软件中的Johnson-Holmquist-2 (JH-2) 模型模拟了单孔爆破裂纹扩展过程，讨论了应力加载速率、自由面距离、节理位置和初始压力等因素对爆破裂纹扩展行为的影响，并用于控制爆破工程指导。Xie等^[13-14]基于该软件中的拉压损伤Riedel-Hiermaier-Thoma (RHT) 模型模拟了地应力下岩体掏槽爆破裂纹扩展过程，发现了不同压力条件下掏槽爆破效果不佳的原因，即静水压力抑制裂纹扩展及非静水压力影响损伤分布，进而提出改变掏槽孔距离及布置方位的优化设计方法。Li等^[15-16]通过数值模拟与图像处理相结合的方法研究了深部岩体爆破的块度特征，分析了地应力对岩体爆破的块度尺寸及形状影响，研究结果表明，随着地应力的增加，岩体爆破块度会更大更圆，且块度粒径分布的范围会更广。

上述研究多集中在初始应力作用下岩体爆破动态力学响应现象分析，有必要深入了解初始应力对岩体爆破破裂机制的影响。此外，如何定量分析岩体爆破裂纹扩展行为和损伤特征与初始应力的关系也有待深入挖掘。本文基于弹性力学建立初始应力作用下单孔爆破动静组合理论模型，分析静载模型应力分布及动载模型应力演化过程，讨论初始应力对岩体爆破损伤的影响，探索破裂机理；采用显式动力学分析软件LS-DYNA，介绍RHT模型中岩石材料参数的标定及修正方法，并通过爆破实验进行模型验证；使用验证的模型分析不同初始应力作用下岩体爆破应力波演化过程及其对裂纹扩展行为的影响，并通过分形维数对爆破损伤量化分析，为深部岩体爆破工程提供理论指导。

1.   动静组合理论模型

初始应力作用下岩体爆破问题可简化为平面应变问题，如图1所示。假定含圆孔无限大平板为均匀连续且各向同性的弹性介质，外侧分别受水平初始应力σ_x和竖直初始应力σ_y作用，圆孔内部受动态载荷p(t)作用，初始应力下岩体单孔爆破模型可看作静态荷载和动态荷载共同作用下的组合模型。

图  1  理论分析几何模型

Figure  1.  Geometrical model for theoretical analysis

下载: 全尺寸图片 幻灯片

1.1   静态应力分布特征

根据弹性力学理论，可分别计算不同初始应力条件下炮孔近场(r=a)与远场(r=10a)的岩体径向应力 ${\sigma }_{r}^{\mathrm{s}}$ 、环向应力 ${\sigma }_{\theta }^{\mathrm{s}}$ 和剪切应力 ${\tau }_{r\theta }^{\mathrm{s}}$ ^[17]

式中：a和r分别为炮孔半径和距炮孔中心的距离；定义竖直方向应力σ_y为压力p，取值20 MPa；K定义为侧压系数，K=σ_x/σ_y。

炮孔近场和远场的应力分布特征如图2所示，近场为孔壁自由面，径向应力与剪切应力均为0。由图2(a)和2(b)可知：水平方向，炮孔近场环向应力随侧压系数的增大呈线性减小的趋势，而炮孔远场几乎不受侧压系数影响；竖直方向，炮孔近场与远场环向应力均随侧压系数的增大呈线性增长趋势，且近场的增长速率更加显著。此外，侧压系数K=0时，孔壁在水平与竖直方向分别出现了压应力和拉应力集中，考虑到岩体抗拉/压强度特性，初始拉应力存在处岩体破坏将更为严重，根据环向应力变化规律可知，拉应力转化为压应力的临界侧压系数K=1/3。炮孔远场径向应力与环向应力分布规律相似，仅在方向上发生90°偏转，如图2(c)所示，即竖直方向上径向应力的大小与侧压系数无关，而水平方向上随侧压系数的增大径向应力也随之增加。剪切应力在水平与竖直方向上均为0，其他方向随侧压系数增大呈先减后增的规律，其方向会在K=1.0前后发生变化，如图2(d)所示。

图  2  炮孔周边静态应力分布

Figure  2.  Stress distribution around the borehole under static load

下载: 全尺寸图片 幻灯片

距炮孔中心不同位置各应力变化过程如图3所示，水平与竖直方向上各应力随距离的增加均趋于平缓变化，不同的是，环向与径向应力主导方向会在侧压系数K=1.0前后发生变化。K>1.0时，水平方向以环向应力为主导，如图3(a)所示，K<1.0时，竖直方向以环向应力为主导，如图3(b)所示。通常认为环向应力是产生径向裂纹的主要原因，因此，各向异性压力条件下径向裂纹在不同方向上的分布有所差异，且随着水平、竖直方向初始应力差的增加，径向裂纹扩展各向异性也会更加明显。

图  3  不同位置处岩体静态应力的变化

Figure  3.  Stress variation of rock mass at different distances under static load

下载: 全尺寸图片 幻灯片

1.2   动态应力演化规律

将炸药产生的爆炸荷载等效为动态荷载p(t)作用于炮孔壁处， $p\left(t\right)={p}_{\mathrm{V}\mathrm{N}}{\left(\mathrm{e}\gamma /n\right)}^{n}{t}^{n}{\mathrm{e}}^{-\gamma t}$ ，其中，p_VN为炸药起爆后炮孔壁处峰值压力，n为控制应力波演化的模型参数，通常取整数， $\gamma$ 为压力衰减指数。根据上述基本假设，极坐标系下岩体中爆炸应力波的控制方程可表示为^[18]：

式中： $\varphi \left(r,t\right)$ 为位移势函数，t为应力波传播时间，v_p为岩体的纵波波速， ${\sigma }_{r}^{\mathrm{d}}(r,t)$ 为爆炸荷载下岩体的径向应力，其动态径向应力和环向应力随时间变化的关系可由下式求得：

式中： ${\sigma }_{\theta }^{{\rm{d}}}$ 为动态环向应力；λ和μ为拉梅常数， $\lambda =\nu E/\left[\left(1+\nu \right)\left(1-2\nu \right)\right]$ ， $\mu =E/\left(2+2\nu \right)$ ，E为弹性模量，v为泊松比。应力波控制方程经过拉普拉斯变换后，可求得式(2)偏微分方程的通解：

式中： ${s}_{1}=p\left(m\right)/\left(\lambda +2\mu \right)$ ， ${s}_{2}=m/{v}_{\mathrm{p}}$ ， ${t}{'}=(r-a)/{v}_{\mathrm{p}}$ ， $w={v}_{\mathrm{p}}/{v}_{\mathrm{s}}$ ，m为拉普拉斯变换参数，v_s为岩体横波波速， $H(t-{t}')$ 为Heaviside函数，K₀和K₁分别为第二类零阶和一阶修正的Bessel函数。

式(4)经过拉普拉斯逆变换后并不能求得其具体的解析表达式，因此需要结合数值反演的方法进行求解^[19]。目前关于数值反演的算法有很多，如Durbin算法、Crump算法和Stehfest算法等，其中Durbin算法稳定可靠性较好^[20-21]，通过该算法可求得其爆炸荷载下岩体中应力波随时间的变化规律：

式中：α可取任意实数，0≤α≤Re(m)；T为求解时间间隔，0≤t≤T/2。

距炮孔中心不同位置处岩体动态应力波演化规律如图4所示。爆炸荷载下动态径向应力先是急剧增加，增长至峰值应力后迅速衰减，且孔壁处应力最大，如图4(a)所示。较大的径向压应力会使得炮孔近区岩体破坏程度加剧，最终形成压缩粉碎区。相较于孔壁处而言，其他位置的径向应力会出现由压缩应力转变为拉伸应力的现象。需要说明的是，径向拉应力的存在会进一步诱导环向裂纹形成，随着距离的增加，峰值径向压应力逐渐减小，且衰减速度相对较快，而峰值径向拉应力几乎没有发生显著变化，整体呈先增后减，最后趋于平缓的变化趋势，如图4(b)所示。动态环向应力与径向应力的变化规律有所区别，如图4(c)所示，其应力时程曲线出现三处应力峰值：爆炸压应力达到峰值后迅速衰减转变为拉应力，达到拉应力峰值后又继续衰减转变为压应力，相较于压应力峰值而言，拉应力峰值相对较大，这也是径向裂纹较为发育的原因之一。峰值环向压应力与拉应力随距离的变化趋势则刚好与峰值径向应力相反，其拉应力峰值相较于压应力较高，压应力随距离的增加无明显的变化，且距炮孔中心5倍的炮孔半径之外，拉伸与压缩应力变化趋势几乎相同，如图4(d)所示。

图  4  不同距离处岩体的动态应力波演化

Figure  4.  Dynamic stress waves evolution in rock mass at various distances

下载: 全尺寸图片 幻灯片

2.   数值计算模型及验证

理论模型将岩石视为弹性介质是为了便于获取数学解析表达式，而在实际工程中，岩体爆破时的力学行为是复杂且非线性的。因此，在弹性理论模型基础上，通过非线性数值模型可弥补理论模型的不足并论证理论模型的合理性。数值模拟计算结果的可靠性很大程度取决于材料本构方程，LS-DYNA软件中常用于模拟岩石材料的有RHT模型、JH-2模型、CSCM模型等。通过数值计算结果发现，RHT材料模型可以更准确地描述深部岩体受到爆炸荷载所表现的力学响应和破坏行为^[22]。

2.1   RHT模型参数标定

RHT模型参数标定需要岩石的基本物理力学参数，通过试验测得花岗岩的力学参数包括：密度（2620 kg·m⁻³）、单轴抗压强度（162 MPa）、泊松比（0.18）、剪切模量（21.9 GPa）、纵波波速（4600 m/s）、横波波速（2820 m/s）。岩石的抗压/拉强度具有明显的应变率效应，在RHT岩石材料模型中，材料的应变强化效应由应变率增强因子F_r为:

式中：p为材料压力；f_c和f_t分别为单轴抗压、抗拉强度； ${\dot{\varepsilon }}_{0}^{\mathrm{c}}$ 和 ${\dot{\varepsilon }}_{0}^{\mathrm{t}}$ 分别为压缩、拉伸荷载下的参考应变率， ${\dot{\varepsilon }}_{0}^{\mathrm{c}}=3.0\times 10^{-5}\;{{\rm{s}}}^{-1}$ ， ${\dot{\varepsilon }}_{0}^{\mathrm{t}}={3.0\times 10}^{-6}\;{\mathrm{s}}^{-1}$ ；β_c和β_t分别表示压缩和拉伸应变率指数， ${\beta }_{\mathrm{c}}=4/(20+3{f}_{\mathrm{c}})$ ， ${\beta }_{\mathrm{t}}=2/(20+{f}_{\mathrm{c}})$ ，根据试验测得岩石单轴抗压强度可知，β_c=0.008，β_t=0.011。

A和N为材料的破坏面参数，用于描述RHT模型中材料的破坏面强度。通过归一化压力满足 $3{p}_{0}^{*} {\text{＞}} {F}_{\mathrm{r}}$ ，得到：

式中： ${\sigma }_{\mathrm{f}}^{*}$ 和 ${p}_{0}^{*}$ 分别为材料归一化压缩强度和压力， ${\sigma }_{\mathrm{f}}^{*}={\sigma }_{\mathrm{f}}/{f}_{{\rm{c}}}$ ， ${p}_{0}^{*}=p/{f}_{\mathrm{c}}$ 。材料处于准静态时，F_r=1，通过Hoek-Brown经验方程可计算不同围压下岩石三轴抗压强度（ ${\sigma }_{1}={\sigma }_{3}+{f}_{\mathrm{c}}\sqrt{\eta {\sigma }_{3}/{\sigma }_{1}+1}$ ，其中σ₁和σ₃分别为最大和最小主应力，η为材料经验系数，η=24），进而求得一系列归一化压缩强度及压力的数据点，再根据式(7)拟合表1中数据点，求得A=2.48，N=0.79，拟合曲线如图5(a)所示。

表  1  不同围压下花岗岩力学参数

Table  1.  Mechanical parameters of granite sample under various confining pressure

σ2/MPa	σ3/MPa	σ1/MPa	${p}_{0}^{*}$	${\sigma }_{\mathrm{f}}^{*}$		σ2/MPa	σ3/MPa	σ1/MPa	${p}_{0}^{*}$	${\sigma }_{\mathrm{f}}^{*}$
0	0	162	0.33	1.00		60	60	569	1.42	3.14
10	10	265	0.59	1.58		70	70	616	1.56	3.37
20	20	342	0.79	1.99		80	80	660	1.69	3.58
30	30	408	0.96	2.33		90	90	703	1.82	3.79
40	40	466	1.12	2.63		100	100	466	1.94	3.98
50	50	520	1.28	2.90

下载: 导出CSV 

| 显示表格

图  5  RHT模型破坏面的相关参数拟合

Figure  5.  Fitting curve of failure surface parameters for RHT model

下载: 全尺寸图片 幻灯片

RHT材料模型相较于其他模型的优点之一是引入了偏应力张量第三不变量及罗德角，其中罗德角因子可以描述失效面子午线失效压缩强度的折减，折减值Q与相对压力的函数关系可表示为：

式中：Q₀为拉伸与压缩荷载下材料子午线半径之比，B为岩石材料罗德角相关系数，参考文献[23]中的试验数据，进行线性回归求得：Q₀=0.68，B=0.05，如图5(b)所示。

RHT模型中材料的破坏积累量由损伤参数D_r确定，其损伤变量和累积塑性应变表达式如下：

式中：Δε_p为塑性应变积累量； ${\varepsilon }_{\mathrm{p}}^{\mathrm{d}}$ 为破坏时的极限应变， ${\varepsilon }_{\mathrm{p}}^{\mathrm{d}}={{D}_{1}\left[{p}_{0}^{*}-(1-D){p}_{\mathrm{t}}^{*}\right]}^{{D}_{2}}$ ，D₁和D₂为损伤变量，D₁=0.04，D₂=1.0， ${p}_{\mathrm{t}}^{*}$ 为破坏截止压力。

RHT模型中剪切应力分量与压应力分量相耦合，压力与应力的关系由Mie-Grüneisen形式的多项式曲线连同p-α孔隙压实曲线关系描述，其状态方程为：

式中：B₀和B₁为材料常数，B₀=B₁=1.22；α₀和ρ₀分别为初始孔隙度与密度； $\kappa$ 为体积应变；e为单位质量的内能；多项式系数A₁，A₂，A₃分别为 ${A}_{1}={\rho }_{0}{c}_{0}^{2}$ ， ${A}_{2}={\rho }_{0}{c}_{0}^{2}(2s-1)$ ， ${A}_{3}={\rho }_{0}{c}_{0}^{2}(3{s}^{2}-4s+1)$ ，其中c₀为材料声速，s为经验系数，计算求得A₁=33.95 GPa, A₂=41.42 GPa, A₃=8.70 GPa。

根据上述经验公式可确定RHT模型大部分参数，其余难以确定且对数值计算结果比较敏感参数，如屈服面参数 ${G}_{\mathrm{c}}^{*}$ 和 ${G}_{\mathrm{t}}^{*}$ ，残余面参数A_f和N_f等，可在RHT模型原始文献^[24]推荐取值的基础上，结合分离式霍普金森压杆(SHPB)动态力学试验进行修正和优化^[25]。SHPB系统中入射杆和透射杆的长度分别为2 000和1500 mm，巴西圆盘岩样的直径和厚度分别为50和25 mm，通过试错法不断调整RHT模型中的参数，直至数值模拟结果与试验结果匹配^[26]，其比较结果如图6和图7所示。

图  6  数值模拟与试验中岩体冲击破坏过程比较

Figure  6.  Comparison of failure process between numerical simulation and test

下载: 全尺寸图片 幻灯片

图  7  数值模拟与试验中岩体应力时程曲线的对比

Figure  7.  Comparison of stress with time between numerical simulation and test

下载: 全尺寸图片 幻灯片

圆盘试样在冲击荷载下的破坏过程及特征如图6所示，试样与杆件接触端裂纹萌生起裂，此后裂纹沿弱面不断扩展延伸，进而两端的裂纹相遇聚结，最终岩样被裂纹贯穿产生宏观大面积裂隙。此外，在接触面会产生较大的破碎区，试样整体为拉伸劈裂破坏。在对比岩样动态破坏过程及特征的基础上，进一步比较了其应力时程曲线，如图7所示。首先进行动态应力平衡验证，如图7(a)所示，试样两端处于较好的受力平衡状态。然后基于三波法求得应力时程曲线，如图7(b)所示，由图可知，岩样的应力时程曲线变化趋势基本一致，且峰值强度相同。值得注意的是，试验中的应力时程曲线相对粗糙，且峰后表现的延性更为显著，这可能与试样的均质性有关。为说明破坏过程与应力时程曲线的对应关系，将破坏过程中的特征点对应地标注在应力时程曲线上（其中t_i、t_g、t_c和t_p分别代表裂纹的起裂、增长、凝聚和贯穿的时间）。总体而言，数值模拟与实验结果整体吻合较好，该模型能够较准确地描述岩样损伤区域范围及动态破坏过程，具体的模型参数如表2所示。

表  2  岩石RHT模型材料参数

Table  2.  RHT model parameters for rock mass

参数名称	符号	值	来源		参数名称	符号	值	来源
密度	ρ0	2620 kg/m3	试验测定		损伤因子	D1	0.04	参考文献[14]
初始孔隙度	α	1.00			状态方程参数	B0	1.22
抗压强度	fc	162 MPa			状态方程参数	B1	1.22
参考压缩应变率	${\dot{\varepsilon }}_{0}^{\rm{c}}$	3.0×10−5 s−1	模型给定		拉伸体积塑性应变分数	${P}_{\mathrm{t}}^{\mathrm{f}}$	0.001
参考拉伸应变率	$\dot{\varepsilon }_{0}^{\mathrm{t}}$	3.0×10−6 s−1			孔隙压实压力	Pcomp	6.00 GPa	理论计算
破坏压缩应变率	${\dot{\varepsilon }}_{\mathrm{c}}$	3.0×1025 s−1			压缩应变率指数	βc	0.008
破坏拉伸应变率	${\dot{\varepsilon }}_{\mathrm{t}}$	3.0×1025 s−1			拉伸应变率指数	βt	0.011
损伤因子	D2	1.00			状态方程参数	T1	33.95 GPa
侵蚀塑性应变	${\varepsilon }_{\mathrm{p}}^{\mathrm{f}}$	2.00	默认取值		状态方程参数	T2	0.00 GPa
孔隙度指数	NP	3.0			Hugoniot多项式系数	A1	33.95 GPa
最小损伤残余应变	${\varepsilon }_{\mathrm{p}}^{m}$	0.012			Hugoniot多项式系数	A2	41.42 GPa
剪切模量减小因子	ξ	0.50			Hugoniot多项式系数	A3	8.71 GPa
相对抗剪强度	${F}_{\mathrm{s}}^{*}$	0.18	试验优化		孔隙坍塌压力	pcrush	108
相对抗拉强度	${F}_{\mathrm{t}}^{*}$	0.06			洛德角相关因子	Q0	0.68
压缩屈服面参数	${G}_{\mathrm{c}}^{*}$	0.50			洛德角相关因子	B	0.05
拉伸屈服面参数	${G}_{\mathrm{t}}^{*}$	0.70			破坏面参数	A	2.48
残余面参数	Af	1.62			破坏面参数	N	0.79
残余面参数	Nf	0.62			弹性剪切模量	G	21.9 GPa

下载: 导出CSV 

| 显示表格

2.2   炸药及空气材料

LS-DYNA软件采用高能炸药模型*MAT_HIGH_EXPLOSIVE_BURN模拟试验中的炸药材料，并结合Jones-Wilkens-Lee (JWL)状态方程来描述其爆炸压力、体积与能量之间的关系：

式中：p_e为爆轰压力，V为相对体积， ${E}_{0}^{\mathrm{e}}$ 为炸药单位体积内能，A_e、B_e、R₁、R₂、ω为材料常数。具体炸药模型的材料参数如表3所示^[27] ，其中v_d和p_CJ分别表示爆轰速度和初始压力。

表  3  炸药模型材料参数

Table  3.  Parameters for the explosive material

ρe/(kg·m−3)	vd /(m·s−1)	pCJ/GPa	${E}_{0}^{\mathrm{e}}$/(J·m−3)	Ae/GPa	Be/GPa	R1	R2	ω
1320	6690	16	7.38×109	586	21.6	5.81	1.77	0.282

下载: 导出CSV 

| 显示表格

爆破时常采用不耦合装药，LS-DYNA中空气介质使用典型材料模型*MAT_NULL，并结合线性多项式状态方程*EOS_LINEAR_POLYNOMIAL定义气体压力、密度和内部能量之间的关系：

式中：p_a为气体压力， ${E}_{0}^{\mathrm{a}}$ 为单位气体体积的内部能量，u为动态黏度系数， $u =\left(\rho /{\rho }_{0}\right)-1$ ，ρ和ρ₀分别是材料的密度和初始密度， ${C}_{4}={C}_{5}=\gamma -1$ ，γ为比热比，通常C₀=C₁=C₂=C₃=C₆=0。具体空气模型的材料参数如表4所示^[27]。

表  4  空气模型材料参数

Table  4.  Mateiral parameters for the air

ρa/(kg·m−3)	C0	C1	C2	C3	C4	C5	C6	${E}_{0}^{\mathrm{a}}$/(kJ·m−3)	V0
1.29	0	0	0	0	0.4	0.4	0	250	1.0

下载: 导出CSV 

| 显示表格

2.3   数值模型验证

数值模型参数确定后，可通过相关试验^[27]验证选取的岩石、炸药及空气材料的合理性。试验中圆柱岩样的高度和直径分别为150和144 mm，炮孔和药卷直径分别为6.45和1.70 mm，药卷周围包裹聚乙烯，并在其外侧嵌套铜管，确保试样爆裂后碎块不会飞溅，便于后期扫描完整断面的裂纹形态。在ANSYS软件中建立相同尺寸的几何模型并划分网格，模型共包含约332万个单元和338万个节点，其详细结构尺寸及横截面局部网格如图8所示。

图  8  数值模型及局部网格

Figure  8.  Configuration of numerical model and local mesh

下载: 全尺寸图片 幻灯片

炸药起爆后产生大量冲击波作用于孔壁，其应力峰值远高于岩体的动态抗压强度，在冲击波作用下，炮孔周围形成压剪粉碎区（Zone Ⅰ）。冲击波穿过耦合介质消耗许多能量进而衰减转变成应力波，其应力峰值低于岩体抗压强度，因而压剪粉碎区的范围不再继续扩展。然而岩石动态抗拉强度相比抗压强度较低，因此在应力波作用下形成了拉剪破碎区（Zone Ⅱ）。应力波在向自由面传播的同时，径向裂纹沿裂隙尖端继续扩展，最终形成边界处的裂隙发育区（Zone Ⅲ）。除径向裂纹外，在拉剪破碎区还有明显的环向裂纹产生，这与实验和理论也是相符的，其比较结果如图9所示。从裂纹区域分布及扩展特征分析，数值模拟与试验和理论结果具有高度一致性。

图  9  岩体爆破裂纹形态理论、试验与模拟结果比较

Figure  9.  Comparison of rock blasting crack patterns among theoretical, experimental and simulated results

下载: 全尺寸图片 幻灯片

为进一步验证数值模型的可靠性及合理性，在比较爆破裂纹扩展过程及分布特征的基础上，进一步对比分析了爆炸压力分布。其结果数据主要来源于理论计算、试验测量及数值模拟，其中，理论计算结果可由经验公式^[22]求得：

式中：p_max为炸药的起爆峰值压力，ρ_s和ρ_e分别为岩体及炸药的密度， v_d为炸药的爆速，v_p为岩体的波速，γ_e为爆轰产物的绝热膨胀系数。压力衰减规律可由下式求得：

式中：p(r)为爆炸压力，ζ为爆炸压力衰减指数。

岩体不同位置处的理论峰值压力可由式(13)与式(14)联合计算求得，而数值模拟结果则是通过沿径向布置若干单元测点输出的峰值压力，试验结果为Banadaki室内试验测量的压力，如图10所示。由图可知，理论计算与数值中峰值压力的衰减趋势基本一致，且试验中测得的压力值也基本分布在其压力衰减曲线上。通过裂纹分布特征及压力衰减规律来看，数值模拟与试验和理论结果较为吻合，故该模型适用于本文的研究工作。

图  10  压力衰减曲线比较

Figure  10.  Comparison of attenuation curve of peak pressure

下载: 全尺寸图片 幻灯片

3.   单孔爆破数值模拟

在数值模型验证的基础上，进行单孔爆破数值模拟计算，分析不同初始压力下岩体静态应力初始化分布及压力演化过程对裂纹扩展行为及损伤特征的影响。模拟主要分为2组进行：各向同性初始压力与各向异性初始压力条件下的单孔岩体爆破，不同工况下的应力加载条件如表5所示。

表  5  初始应力加载条件

Table  5.  Initial stress conditions in numerical simulation

应力状态	工况	埋深/m	σx/MPa	σy/MPa		应力状态	工况	埋深/m	σx/MPa	σy/MPa
各向同性初始压力	E-1	0	0	0		各向异性初始压力	A-1	750	5	20
	E-2	375	10	10			A-2	750	10	20
	E-3	750	20	20			A-3	1125	30	20
	E-4	1125	30	30			A-4	1500	40	20

下载: 导出CSV 

| 显示表格

为验证理论模型的合理性，数值模型采用平面单层网格，约束垂直平面位移确保计算在平面应变条件下进行。模型几何尺寸的边长为4 m，其中炮孔半径为21 mm，不耦合系数为1.3。模型中单元数量共计约108万个，单元的平均尺寸为4 mm。为保证计算效率，且要解决爆炸大变形可能引起的网格畸变问题，数值模型中岩体单元采用拉格朗日算法，炸药与空气单元则采用任意拉格朗日欧拉算法，通过流固耦合关键字与多物质组关键字实现爆炸荷载的有效传递。在模型四周施加围压来模拟初始应力，并设置无反射边界条件以模拟无限区域岩体，在距炮孔中心不同位置处布置若干测点，记录压力变化过程，数值模型及压力测点布置如图11所示。

图  11  数值模型及压力测点布置

Figure  11.  Numerical model and arrangement of the measurement points

下载: 全尺寸图片 幻灯片

3.1   静态应力初始化

初始应力下岩体爆破数值模拟分两步计算：第一步在边界施加初始应力，待应力初始化稳定后再进行爆破计算。应力初始化方法有很多，如动力松弛法、显-隐式转化法和dynain文件法等，综合考虑模拟效果和计算效率，采用dynain文件法更稳定，该方法可用显/隐式求解功能实现应力初始化^[28]。使用关键字*INTERFACE_SPRINGBACK_LSDYNA将应力初始化结果输出，然后通过关键字*INCLUDE将第一步的计算结果文件包含进来，再进行第二步计算，具体的计算流程如图12所示。

图  12  初始应力下岩体爆破计算流程

Figure  12.  Flow chart for the rock blast under initial stress

下载: 全尺寸图片 幻灯片

为说明施加于岩体上的初始应力处于稳定状态，以工况E-3为例进行分析。选取炮孔周边的四个测点，输出其压力时程曲线，如图13所示。初始应力时程曲线主要包括两个阶段——应力增长阶段与稳定阶段，加载初期应力以一定速度稳定增长，增长至预设的压力目标值时便趋于稳定状态。此外，不同测点的应力时程曲线基本重合，因此施加于岩体上的初始应力是稳定可靠的。

图  13  应力初始化验证

Figure  13.  Verification of stress initialization

下载: 全尺寸图片 幻灯片

在后处理软件LS-PREPOST中建立局部柱坐标系，绘制岩体径向、环向及剪切应力的分布云图，如图14~图18所示。云图可直观地显示炮孔周围岩体的应力分布，如侧压系数K=0时，竖直方向上径向与环向应力出现局部拉应力集中，剪切应力在一、三象限为压应力，二、四象限为拉应力，如图14所示；侧压系数K=1.0时，径向与环向应力在炮孔壁处成同心圆状分布，且不存在剪切应力，如图16所示。上述的模拟结果与1.1节中理论模型分析的结果高度一致，再次验证了数值模型的合理性及应力初始化方法的可靠性。

图  14  应力初始化模拟结果 (K=0)

Figure  14.  Numerical results of the stress initialization (K=0)

下载: 全尺寸图片 幻灯片

图  15  应力初始化模拟结果 (K=0.5)

Figure  15.  Numerical results of the stress initialization (K=0.5)

下载: 全尺寸图片 幻灯片

图  16  应力初始化模拟结果 (K=1.0)

Figure  16.  Numerical results of the stress initialization (K=1.0)

下载: 全尺寸图片 幻灯片

图  17  应力初始化模拟结果 (K=1.5)

Figure  17.  Numerical results of the stress initialization (K=1.5)

下载: 全尺寸图片 幻灯片

图  18  应力初始化模拟结果 (K=2.0)

Figure  18.  Numerical results of the stress initialization (K=2.0)

下载: 全尺寸图片 幻灯片

3.2   压力演化过程

由理论分析可知，环向应力对径向主裂纹的扩展行为具有决定性作用，为进一步探究裂纹的扩展机制，本节主要选取E-1（无初始压力）、E-4（各向同性压力）及A-4（各向异性压力）三种工况，分析环向应力的演化过程，如图19~图21所示。工况E-1及E-4下环向应力演化规律相似，如图19与图20所示。爆炸初期，炮孔周边产生较强的爆轰波，环向应力呈圆形辐射状传播，接着冲击波能量不断衰减，环向应力逐渐降低并趋于稳定状态。此外，相较于无初始压力而言，各向同性压力作用下环向拉应力波的扩展范围及扩展时间受到明显的抑制作用：无初始压力时，其扩展时间停留在0.5 ms前后，40 MPa初始压力时，其扩展时间仅停留在0.3 ms 前后，且其扩展范围也大幅度缩小。相较于各向同性压力条件，各向异性压力下环向应力演化过程有显著差别，尤其在0.2 ms后，环向拉应力波在水平与竖直方向分布出现差异，差异性随时间增加愈加明显，其传播方向逐渐集中于最大主应力方向，如图21所示。

图  19  无初始压力下环向应力的演化过程 (E-1)

Figure  19.  Evolution of hoop stress without initial pressure (E-1)

下载: 全尺寸图片 幻灯片

图  20  各向同性压力下环向应力的演化过程 (E-4)

Figure  20.  Evolution of hoop stress under equibiaxial pressure (E-4)

下载: 全尺寸图片 幻灯片

图  21  各向异性压力下环向应力的演化过程 (A-4)

Figure  21.  Evolution of hoop stress under anisotropic pressure (A-4)

下载: 全尺寸图片 幻灯片

根据监测点记录的数据，绘制相应的环向应力演化曲线，如图22所示。各向同性压力下水平与竖直方向监测点是中心对称的，因此其环向应力时程曲线基本重合，如图22(a)所示。而各向异性压力下水平与竖直方向上的环向应力时程曲线有所区别，主要表现在拉压应力峰值的差异性，最大主应力方向上的压应力削弱，而拉应力增强，如图22(b)所示。这可以很好地解释裂纹倾向于最大主应力方向扩展的原因。值得注意的是，不同初始压力下，其环向应力的变化规律具有相似性，即爆炸初期因其强烈的冲击作用产生较高的压应力，岩体发生压缩破坏，随后迅速衰减转变为拉应力，岩体发生拉伸破坏，最后趋于稳定状态，其演化规律与1.2节中理论模型分析结果一致。

图  22  不同压力条件下环向应力时程曲线

Figure  22.  Time histories of hoop stress under various pressure conditions

下载: 全尺寸图片 幻灯片

不同距离处峰值环向应力的变化特征如图23所示。各向同性压力下，峰值环向压应力的变化趋势一致，即随着距离的增加而降低，降低速率呈先增后减的规律。同样地，峰值环向拉应力随距离的增加也显著降低，且随着初始压力的增加，其降低速率变得更加明显，如图23(a)所示。各向异性压力下，峰值环向拉应力在水平与竖直方向有所差异，最大主应力方向峰值环向拉应力相对较大，且随着距离的增加，差异性更加明显，尤其在0.3 m外，如图23(b)所示。

图  23  不同距离峰值环向应力变化

Figure  23.  Peak hoop stress variation at various distances

下载: 全尺寸图片 幻灯片

3.3   裂纹扩展行为

各向同性压力下岩体爆破裂纹扩展特征如图24所示。径向主裂纹在各方向分布匀称，大致呈圆形裂纹带状分布，如图24(a)所示。其扩展半径随压力的增加而不断减小，且减小速率逐渐降低，当围压增加至40 MPa左右，裂隙扩展半径几乎不再发生明显变化，此外，径向主裂纹的数量随压力增大也会显著减少，与裂纹扩展半径有相似的变化规律，如图24(b)所示。经过指数函数拟合，可获得圆形裂纹带扩展半径L_c(m)与围压p(MPa)之间的关系：

图  24  各向同性压力下爆破裂纹扩展特征

Figure  24.  Characteristics of blasting crack propagation under different equibiaxial pressures

下载: 全尺寸图片 幻灯片

径向裂纹的扩展过程如图24(c)所示。在0.1 ms内，不同压力下裂纹扩展速度几乎相同，这是因为爆炸初期的爆轰压力远大于初始应力，初始应力对岩体爆破裂纹扩展的抑制作用被弱化。而0.1 ms后，裂纹的扩展速度发生了明显的变化，随初始压力的增加其扩展速度明显降低，且随着时间的推移，这种趋势愈加明显。初始压应力对裂纹的扩展具有明显的抑制作用，深部岩体在爆破时可能会因为高地应力挤压作用导致破碎困难，此时建议结合预裂爆破等手段进行卸压以削弱地应力的抑制作用，从而提高爆破破岩效率。

各向异性压力条件下岩体爆破裂纹扩展特征如图25所示。爆破裂纹在水平及竖直方向大致对称分布，径向主裂纹呈椭圆形裂纹带扩展，如图25(a)所示。椭圆裂纹带横轴长度L_x与竖轴长度L_y的比值定义为椭圆裂纹带轴长比c，即c=L_x/L_y。随侧压系数的增大，横轴裂纹长度无明显变化，竖轴裂纹长度显著减小，不同方向主裂纹长度的变化规律与初始应力具有良好的一致性。此外，椭圆裂纹带轴长比会随侧压系数的增加而增大，且增长的速率也逐渐提高，如图25(b)所示，也就是说，当初始水平与竖直应力差增大时，主裂纹带分布的各向异性会更加显著。经过指数函数拟合，可获得椭圆裂纹带的轴长比c与侧压系数K的关系：

图  25  各向异性压力下爆破裂纹扩展特征

Figure  25.  Characteristics of blasting crack propagation under different anisotropic pressures

下载: 全尺寸图片 幻灯片

径向裂纹在水平与竖直方向的扩展过程如图25(c)所示，其扩展规律与各向同性压力下类似，在此不再赘述。需要说明的是，不管侧压系数如何变化，径向主裂纹总是会倾向于最大主应力方向扩展。初始应力对裂纹扩展长度抑制的同时，对其扩展的方向也有诱导作用，因此，在实际工程中，若要对岩体进行爆破卸压，炮孔应尽可能沿最大主应力方向布置，以便于促进炮孔间裂隙的贯通。

3.4   损伤分形维数

岩石损伤破裂过程具有明显的分形特征，因此，岩石的损伤演化可以看作是分形的演化过程^[29]。同样，爆破裂纹的扩展过程也可以用分形维数来表征。本文采用计盒分形维数法定量描述初始应力下爆破裂纹的分布形态，分形维数D的数学表达式为^[30]

式中：δ_i 为方盒的边长，N(δ_i)为覆盖爆破裂纹的方盒数量。

将爆破裂纹损伤图像转换成黑白二值图像，然后将二值图像划分成若干个边长为δ_i的方盒，最后再统计覆盖爆破裂纹的所有盒子数量N(δ_i)，如图26所示。当方盒长度δ_i趋于0时，方盒数量N(δ_i)与方盒长度δ_k分别取对数后的比值将趋于D_b，对数据点用最小二乘法拟合，其线性方程为：

图  26  爆破裂纹盒子划分

Figure  26.  Box division for blasting cracks

下载: 全尺寸图片 幻灯片

式中：线性方程斜率D_b为爆破裂纹的分形维数，S为线性拟合曲线的纵轴截距。

此外，不同初始应力下岩体爆破破裂的分形损伤ω可用公式ω=(D−D₀)/(D_max−D₀)来量化比较，其中，D₀为岩石爆破前的初始分形维数，D_max为岩石完全损伤时的最大分形维数。对于完整岩石而言，初始分形维数D₀=0。同时，本文理论模型及数值模型是基于二维平面应变条件下进行的，岩石破裂也是二维平面的破裂，因此最大分形维数D_max=2。也就是说，初始应力下岩体爆破破裂的分形损伤ω与分形维数D_b呈线性相关，因此仅选择一个分形参量的变化规律分析即可。

不同初始应力作用下岩体爆破裂纹分形维数拟合线如图27和图28所示，由图可知，所有的分形维数数据都可以用一条直线很好地进行拟合。各向同性压力下，随着围压的不断增加，爆破裂纹的分形维数不断减小，初始压力为0时分形维数最大：D_b=1.403，当围压增长至10、20和30 MPa时，对应的爆破裂纹分形维数分别为1.341、1.214和1.087，如图27所示。各向异性压力下，侧压系数K为0.25、0.5、1.5和2.0时，对应的爆破裂纹分形维数分别为1.384、1.247、1.131和1.078，随着侧压系数的增加，分形维数呈现降低的变化趋势，如图28所示。

图  27  各向同性压力下爆破裂纹分形维数拟合线

Figure  27.  Fractal dimension of blasting cracks and its fitting lines under equibiaxial pressure

下载: 全尺寸图片 幻灯片

图  28  各向异性压力下爆破裂纹分形维数拟合线

Figure  28.  Fractal dimension of blasting cracks and its fitting lines under anisotropic pressure

下载: 全尺寸图片 幻灯片

4.   结　论

本文基于弹性力学建立了初始应力下单孔爆破动静组合理论模型，研究了静态应力分布及动态应力演化特征，揭示了岩体爆破破裂机制。通过相应的数值模拟，讨论了初始应力对岩体爆破压力演化及裂纹扩展行为的影响，并引入分形维数定量表征岩体爆破损伤特性，结论如下：

(1) 动静组合理论模型可以较好地反应初始应力下岩体爆破时的静态应力分布及动态应力演化过程，解释其爆破碎裂机理，并在相应数值模型中得到验证；

(2) 初始应力下岩体爆破裂纹的扩展行为可以通过环向应力分析；各向同性压力下环向应力分布完全对称，各向异性压力下环向应力主要集中于最大主应力方向；初始应力对裂纹扩展范围抑制的同时，对裂纹扩展方向也具有导向作用，随着水平与竖直初始应力差的增加，裂纹分布差异性显著增强；合理调整环向应力分布，可以改善岩体爆破碎裂效果；

(3) 各向同性压力下爆破裂纹带近似呈圆形，裂纹扩展半径及数量随压力的增加显著减小；各向异性压力下爆破裂纹带近似呈椭圆形，裂纹带的轴长比随侧压系数的增加显著增大；高地应力下岩体爆破破碎效果不佳时，建议结合预裂技术先卸压再爆破；

(4) 初始应力下岩体爆破破裂过程具有明显的分形特征，各向同性压力下分形维数随围压的增加而减小，各向异性压力下呈现相同的变化规律，分形维数可以定量表征岩体爆破破裂的损伤特性。