预裂爆破是岩体开挖工程中广泛使用的一种控制爆破技术。在主爆区爆破前沿设计轮廓线预先生成一条贯通裂缝，从而改善主爆区岩体破碎效果并降低围岩爆破振动^[1-2]。预裂爆破效果不仅与岩体力学性质和布孔方式相关，更依赖于岩体赋存的地质环境。随着地下工程的大量开发与利用，岩体开挖深度不断增加。浅部岩体所受地应力较低，地压对岩体预裂爆破影响较小。然而，深部岩体所受地应力较高，高地压强约束条件下，岩体预裂爆破会表现出与浅部不同的力学行为^[3-4]。

20世纪60年代，Langefors等^[5]首次开展了光面爆破试验，并取得了良好的控制爆破效果。在光面爆破技术的基础上，进一步发展出了预裂爆破技术。这2种技术都被归类为周边爆破技术，或者称为轮廓爆破技术^[6]。Nicholls等^[7]考虑了地应力的影响，率先开展了露天花岗岩原位预裂爆破试验，研究发现爆破主裂纹的扩展方向始终与水平应力方向保持一致。受到这一研究的启发，Yang等^[8]开展了单向压力下岩体深孔预裂爆破模型试验，发现当压力方向与炮孔连线垂直时，预裂成缝效果较差，而当压力方向与炮孔连线平行时，裂纹在孔间容易贯通。在工程实践方面，Lu等^[9]总结了国内大型地下厂房的岩体开挖程序及轮廓爆破方式，讨论了地应力对深埋岩体预裂爆破的影响，并针对水平应力较高的工程岩体提出了先拉槽再预裂的爆破开挖方案。Hu等^[10]通过数学模型分析了深谷高边坡预裂爆破损伤机理，发现地应力对岩体爆破裂纹扩展有明显的抑制作用，并指出传统的预裂爆破方法无法实现岩体理想的开挖效果。针对这一问题，众多学者又相继开展了预应力下聚能爆破^[11]和切缝爆破^[12]等新型控制爆破技术的研究，并取得了丰硕的成果。

数值模拟在岩体爆炸动力学研究领域具有显著的优势，它不仅能够丰富研究数据和信息，还能为深入解析岩体爆破损伤机理提供有力的理论支持^[13]。白羽^[14]在COMSOL Multiphysics软件环境下进行有限元编程，建立了地应力下岩体双孔爆破损伤力学模型，并模拟了爆破过程中的裂纹扩展行为，发现爆炸应力波在裂纹萌生阶段起主导作用，而爆生气体压力则会促使裂纹进一步扩展和贯通。杨建华等^[15]将多孔爆破简化为双孔爆破模型，采用FEM (finite element method)-SPH（smoothed particle hydrodynamics）耦合法模拟了孔间裂纹传播及贯通过程，并通过分析炮孔周围环向应力场的演化和分布特征，揭示了地应力下多孔爆破的破岩机制。Li等^[16]探讨了地应力作用下深埋隧洞轮廓爆破的超/欠挖问题，通过对隧洞顶板和侧壁进行受力分析，总结了爆破损伤特征，结合数值模拟结果，讨论了不同压力条件下隧洞周边爆破裂纹的扩展行为，并提出了相应的优化方案。Lu等^[17]利用自主开发的数值分析软件FDEM (finite-discrete element method)，研究了初始应力、孔间距和径向不耦合系数对爆炸应力叠加和孔间开裂特征的影响，为地应力影响下的轮廓爆破提供了思路，并通过锦屏二级水电站隧洞爆破开挖工程进行了有效性验证。

本文中，在前人研究的基础上，进一步探究地应力对岩体预裂爆破裂纹扩展行为的影响。基于弹性力学中平面应变问题假设，建立地应力下岩体预裂爆破的理论模型。通过Laplace变换和数值反演的方法计算岩体爆炸应力的传播和衰减过程，并分析初始应力对岩体预裂爆破压力场分布的影响。结合应力波叠加破坏理论推导岩体预裂爆破孔间成缝时的判别条件，即不同地应力条件下，孔间成缝时炮孔间距与装药不耦合系数之间的关系。此外，依据室内小型岩石爆破实验，结合爆破裂纹分布特征与压力衰减曲线验证有限元模型的合理性。通过验证的模型进一步模拟静水压力和非静水压力条件下岩体预裂爆破的损伤过程，并依据矢量演化揭示岩体预裂爆破时的成缝机制。结合理论计算结果对岩体预裂爆破提供参数优化，优化方案不仅可以保证岩体预裂成缝，同时，在很大程度上减少对围岩的破坏，以期为深部岩体预裂爆破提供理论参考。

1.   理论分析

地应力下岩体预裂爆破可简化为初始压力下双孔爆破引起的平面应变问题，简化的理论几何模型见图1。模型包括岩体、炸药和空气3部分，炮孔和装药直径分别为D_b和D_e，其比定义为装药不耦合系数K_dec，K_dec=D_b/D_e。两孔间距为S_hol，将炮孔中心连线方向上的中点定义为原点O。假设含两圆孔的无限大平板为均匀连续且各向同性的弹性介质，平板外侧分别受到水平方向压力σ_hor和竖直方向压力σ_v的作用，孔壁受到爆炸压力的作用。因此，地应力下岩体预裂爆破可以看作是静态荷载与动态荷载共同作用的结果。

图  1  预裂爆破理论模型示意图

Figure  1.  Schematic diagram of theoretical model for presplit blasting

下载: 全尺寸图片 幻灯片

1.1   静态应力分布

根据弹性力学理论，炮孔周边岩体中的应力分布状态可由吉尔西公式获得，其中计算圆孔周边弹性体中任意一点的应力解析表达式如下^[18]：

式中：${\sigma }_{r}^{\mathrm{s}}$、${\sigma }_{\theta }^{\mathrm{s}}$和${\tau }_{r\theta }^{\mathrm{s}}$分别为径向应力、切向应力和剪切应力，其中，上标“s”表示静态应力场，下标“r”和“θ”分别表示极径和极角。k和ψ分别为侧向压力系数和距离比，定义式如下：

式中：a为炮孔半径，L为距炮孔连线中心O的距离。将竖直方向的压力σ_v定义为地应力P，根据地应力随岩体埋藏深度的变化规律^[19]，750 m左右埋深的地应力约为20 MPa。

单个炮孔周边的切向应力分布如图2所示。由图2(a)可知，静水压力条件下，切向应力随初始压力的升高呈线性增高的趋势。当初始压力唯一时，孔壁周边的切向应力在各方向上均相等。非静水压力条件下，切向应力分布如图2(b)所示。水平方向上，切向应力随侧压系数的增大呈线性降低的趋势，当侧压系数为0和2时，切向应力分别取最大值6P和最小值2P。竖直方向上，切向应力表现出与水平方向相反的变化规律，且该方向上的变化幅度相对较大。值得注意的是，单向压力时，在孔壁竖直方向出现了局部拉应力集中现象。考虑到岩体的抗拉强度远低于抗压强度，因此，拉应力集中区附近的岩体更容易发生破裂。根据切向应力随侧压系数的变化规律可知，竖直方向上拉应力转换至压应力的临界侧压系数为1/3，水平方向上的临界侧压系数为3。

图  2  地应力作用下炮孔周边切向应力分布

Figure  2.  Tangential stress distribution around the borehole under in-situ stress

下载: 全尺寸图片 幻灯片

2个炮孔连线方向上应力分布如图3所示。不同压力下中点O附近处的径向与切向应力均处于稳定状态，而在孔壁周边区域，其应力分布差异显著。由图3(a)可知，静水压力下，径向应力和切向应力远离炮孔时分别呈提高和降低趋势。当初始压力一定时，中点附近处的径向应力与切向应力基本相等。随着初始压力的升高，中点附近处的切向应力也随之升高，较高的切向压应力会削弱爆炸后叠加的切向拉应力，进而阻碍孔间裂纹贯通。非静水压力条件下静态应力的分布如图3(b)所示，径向应力和切向应力均在远离炮孔时趋于平缓变化，但主导应力会在侧压系数k=1.0前后发生变化。k<1.0时，中点附近以切向应力为主导；k>1.0时，径向应力在中点附近占主导地位。

图  3  炮孔中心连线方向上静态应力分布

Figure  3.  Static stress distribution along the center line of two boreholes

下载: 全尺寸图片 幻灯片

1.2   动态应力的传播和衰减

根据理论模型的基本假设，爆炸应力波在岩体中的传播过程可简化为平面应变问题，极坐标系下爆炸荷载所引起的应力变化满足以下运动方程^[20]：

式中：${\sigma }_{r}^{\mathrm{d}}$和${\sigma }_{\theta }^{\mathrm{d}}$分别为爆炸荷载下炮孔周边岩体中的动态径向应力和切向应力，u_d为岩体内某一质点的径向位移，ρ为材料密度，t为应力波传播时间。

线弹性材料在平面应变问题下，某一质点的应力状态可由广义的胡克定律表示：

式中：λ和G为拉梅常数，可根据材料的弹性模量E和泊松比μ分别计算获得：$\lambda =E\mu /\left(1-\mu -2{\mu }^{2}\right)$，$G=E/\left(2+2\mu \right)$。将式(4)代入式(3)可以得到控制方程^[21]：

式中：c_p为纵波波速，可根据材料密度和拉梅常数计算，${c}_{\mathrm{p}}={\left[\left(\lambda +2G\right)/\rho \right]}^{0.5}$。

控制方程的边界条件如下：

为便于微分方程进行Laplace变换，将指数函数形式的爆炸压力曲线近似为多段线性函数^[22]：

式中：p(t)为等效爆炸荷载，t_r和t_con分别为爆炸压力的上升和持续时间，p_max为爆炸压力峰值。p_max的表达式为：

式中：ρ_e和v_de分别为炸药的密度和爆轰速度，且ρ_e=1320 kg/m³，v_de=6690 m/s；γ为等熵指数，通常取3.0；η为绝热膨胀系数，通常取1.5。

应力与位移关系式(4)经过一系列Laplace变换后可以转换得到相应的应力数学表达式^[23]：

式中：上划线“⁻”表示Laplace变换，m为Laplace变换参数，p(m)为爆炸荷载p(t)的Laplace变换式，K₀和K₁分别表示第二类0阶和1阶的Bessel函数。

式(9)即为爆炸荷载下炮孔周边岩体动态应力场的Laplace空间解，结合Stehfest算法进行数值反演，便可求得其具体的解析解^[24]。此处炮孔直径D_b和装药直径D_e分别设定为40和20 mm，孔距S_hol设定为15倍的炮孔直径，即600 mm。根据应力波叠加理论，可求得爆破引起的动态切向应力的传播和衰减过程曲线，如图4所示。单孔爆破时，应力时程曲线出现3处峰值，如图4(a)所示。炸药起爆后，孔壁上的切向应力迅速达到压力峰值，此后压应力转变为拉应力，拉应力达到峰值后又继续衰减转变为压应力，并逐渐趋于稳定。双孔爆破时，动态切向应力的传播和衰减过程与单孔爆破时基本一致。尤其在孔壁附近，其应力时程曲线几乎相同，如图4(b)所示。

图  4  不同爆破荷载下动态切向应力演化

Figure  4.  Evolution of the tangential stress under different blasting loads

下载: 全尺寸图片 幻灯片

需要说明的是，在中点附近，单孔与双孔爆破所引起的切向应力差别显著，如图5(a)所示。当距离比r_i≤12时，双孔爆破产生的应力波无明显叠加现象，而在中点O处，其叠加效应最显著，叠加后的切向拉应力约为单孔爆破时的2倍。动态切向拉应力对孔间裂缝的形成起到关键的控制作用，下面分析炮孔连线上切向拉应力峰值分布特征。由图5(b)可知，切向拉应力曲线总体以O为中心，左右呈对称分布，中点附近呈W形变化，最小值出现在距中点2.5倍的炮孔半径处。因此，岩体预裂爆破时，最难形成裂缝的不是中点，而是中点附近的某一点，求得该点处的切向拉应力，并与岩石的抗拉强度比较，便可判定预裂爆破时孔间裂缝能否贯通。

图  5  动态切向应力的比较与孔间峰值应力分布

Figure  5.  Dynamic tangential stresses and distribution of maximum stress in adjacent holes

下载: 全尺寸图片 幻灯片

1.3   预裂爆破孔间成缝准则

相邻炮孔间应力波传播和叠加过程如图6所示。应力波叠加类型分为平行波与斜波2种，其波阵面可分解为径向压缩应力分量σ_c和切向拉伸应力分量σ_t，如图6(a)所示。当波阵面在中点O附近相遇时，炮孔连线上的平行波首先发生叠加。由于切向应力分量方向相同，应力波叠加后拉应力显著增强，在拉应力作用下，岩体会沿相反方向发生变形。斜波叠加略滞后于平行波，由于其波阵面切向分量不在同一方向，此时会与另一个波阵面的法向分量相互抵消，进而使得拉应力降低，如图6(b)所示。值得注意的是，地应力在竖直方向会产生压应力，根据应力波叠加破坏理论^[25]，给出预裂爆破孔间成缝准则：当竖直方向上叠加的合外力大于岩体抗拉强度时，便会发生拉伸破裂，进而孔间裂纹扩展连接，最终贯通成缝：

图  6  相邻炮孔间应力波传播和叠加示意图

Figure  6.  Diagrams of stress wave propagation and superposition in adjacent holes

下载: 全尺寸图片 幻灯片

式中：${\sigma }_{\theta ,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{d}}$为切向拉应力峰值的最小值；f_t为岩体抗拉强度；σ_g为岩体受到的地应力，可以根据σ_g=[(σ_max+σ_min)+(σ_min−σ_max)cos(2α)]/2确定，其中σ_max和σ_min分别为最大主应力和最小主应力，α为相邻2个炮孔的连线与水平方向的夹角。

基于上述理论分析，将孔径固定为40 mm，改变不耦合系数和孔距，进而计算不同地应力下岩体预裂爆破孔间成缝时孔距与装药直径的关系，如图7所示。当地应力一定时，孔距越大，装药直径越大，说明孔距扩大时，需要更多的药量产生更高的爆炸压力来保证预裂爆破的成功。此外，孔距与装药直径间近似呈线性变化关系，随着地应力的升高，曲线斜率逐渐减小，说明地应力越高，孔距增大相同的距离所需增加的装药量也越多。类似地，当孔距一定时，地应力水平越高，装药直径越大，说明地应力的存在会对爆破裂纹产生抑制作用，地应力越高，就需要更多的爆炸能量，从而保证爆破裂纹能够顺利贯通。对于不同孔距而言，地应力水平较低时，装药直径对爆破裂纹扩展行为的影响较显著。当装药直径一定时，地应力水平越高，孔距越小。说明在较高地应力水平中开展预裂爆破时，在不改变开挖轮廓线的前提下，又要保证预裂缝能够顺利贯通，相邻2个炮孔间的距离应当减小，意味着需要开凿更多的预裂孔来实现预期的爆破效果。

图  7  不同地应力下孔距与装药直径之间的关系

Figure  7.  Relationship between hole spacing and charge diameter under different in-situ stresses

下载: 全尺寸图片 幻灯片

2.   数值模型验证

理论模型可提供具体的数学解析表达式，但与岩石材料所表现出的动态力学响应有所区别。为此，在线弹性理论模型的基础上，结合非线性有限元模型，进一步探究地应力下岩体预裂爆破损伤特征。基于Banadaki^[26]开展的花岗岩小型爆破实验，确定材料参数并验证数值模型。该实验选用材料较致密的柱状花岗岩，其高度H_roc和直径D_roc分别为150和144 mm，如图8(a)所示。在岩样中心钻孔，孔径D_b为6.45 mm，钻孔中心装有直径D_e为1.7 mm的药卷，药卷外侧包裹聚乙烯材料，其外侧直径D_pol为4.5 mm。为防止岩样爆炸后崩解，在炮孔外侧嵌套铜管，其外径D_cop和内径D_a分别为6.45和5.25 mm。对爆破后的岩样进行分段切片处理，切片位置分别布置在距岩样顶面25 mm (#T)、75 mm (#M)和125 mm (#B)处。切片后用染料浸染，将其放置在高强度紫外光下观察切片中的爆破裂纹分布形态。为尽可能与实验条件保持一致性，在ANSYS软件中建立了相同尺寸的三维立体模型，该模型划分了约332万个六面体网格和338万个节点。数值模型中包含岩石、铜管、空气、聚乙烯和炸药5种材料，模型的结构尺寸及局部网格如图8(b)所示。

图  8  爆破实验设置^[26]及相应的有限元模型

Figure  8.  Experimental setup for rock lab-scale blasting^[26] and the corresponding finite element model

下载: 全尺寸图片 幻灯片

2.1   材料参数

Riedel-Hiermaier-Thoma (RHT)模型中充分考虑了岩石材料在爆炸荷载下的拉压损伤特性，RHT模型表现出的围压、应变率及损伤软化效应，可以较好地反映岩体在地应力下的爆炸力学响应特征，因此，被广泛应用于岩石爆炸损伤破裂过程的模拟。RHT模型中考虑了材料孔隙度与压力之间的关系：

式中：p_roc为岩石材料所受的压力，B₀和B₁为材料常数，ρ_roc和μ_roc分别为岩石材料的密度和体积应变，e为体积内能，φ₁和φ₀分别为材料当前孔隙度和初始孔隙度，A₁、A₂和A₃为材料的Hugoniot多项式系数。

RHT模型通过应力极限面描述在不同加载阶段材料强度的变化特征。弹性屈服面前，材料表现为线弹性变形，弹性屈服面后，材料表现为塑性变形。当强度达到极限时，材料内部开始产生损伤并不断累积，累积到一定程度材料发生软化，持续加载后最终完全破坏。其损伤参数定义为：

式中：Δε_p为累积塑性应变；${\varepsilon }_{\mathrm{p}}^{\mathrm{f}}$为破坏时材料的塑性应变，且$\varepsilon_{\mathrm{p}}^{\mathrm{f}}=D_1\left[p^*-\left(1-D\right)p_{\mathrm{t}}^*\right]^{D_2}$，其中D₁和D₂为损伤因子，${p}_{}^{*}$为归一化压力，${p}_{\mathrm{t}}^{\mathrm{*}}$为破坏截止压力。

RHT模型部分参数可根据力学实验和经验公式确定，其余参数可以参考相关文献进行取值。本文使用的岩石RHT模型参数如表1^[27]所示。

表  1  岩石RHT模型材料参数^[27]

Table  1.  Parameters of RHT model for rock mass^[27]

参数	取值		参数	取值		参数	取值
密度ρroc/(kg·m−3)	2620		状态方程参数B0	1.22		剪切模量减小因子ξ	0.50
抗压强度fc/MPa	162		状态方程参数B1	1.22		参考压缩应变率${\dot{\varepsilon }}_{0}^{\mathrm{c}}$/s−1	3.0×10−5
初始孔隙度φ0	1.00		侵蚀塑性应变${\varepsilon }_{\mathrm{s}}^{\mathrm{f}}$	2.00		参考拉伸应变率${\dot{\varepsilon }}_{0}^{\mathrm{t}}$/s−1	3.0×10−6
损伤因子D1	0.04		压缩屈服面参数${G}_{\mathrm{c}}^{\mathrm{*}}$	0.50		破坏压缩应变率${\dot{\varepsilon }}_{\mathrm{c}}$/s−1	3.0×1025
损伤因子D2	1.00		拉伸屈服面参数${G}_{\mathrm{t}}^{\mathrm{*}}$	0.70		破坏拉伸应变率${\dot{\varepsilon }}_{\mathrm{t}}$/s−1	3.0×1025
破坏面参数A	2.48		洛德角相关因子B	0.05		最小损伤残余应变${\varepsilon }_{\mathrm{p}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$	0.012
破坏面参数N	0.79		洛德角相关因子Q0	0.68		孔隙坍塌压力pcrush/MPa	108
残余面参数Af	1.62		压缩应变率指数βc	0.008		孔隙压实压力pcomp/GPa	6.00
残余面参数Nf	0.62		拉伸应变率指数βt	0.011		拉伸体积塑性应变分数${f}_{\mathrm{t}}^{\mathrm{p}}$	0.001
相对抗剪强度${F}_{\mathrm{s}}^{\mathrm{*}}$	0.18		弹性剪切模量G/GPa	21.9		Hugoniot多项式系数A1/GPa	33.95
相对抗拉强度${F}_{\mathrm{t}}^{\mathrm{*}}$	0.06		状态方程参数T1/GPa	33.95		Hugoniot多项式系数A2/GPa	41.42
孔隙度指数Npor	3.0		状态方程参数T2/GPa	0.00		Hugoniot多项式系数A3/GPa	8.71

下载: 导出CSV 

| 显示表格

LS-DYNA程序使用*MAT_NULL表征聚乙烯材料，结合Grüneisen状态方程定义材料压力、密度与初始内能之间的关系，表达式如下：

式中：p_pol为压力，c_pol和E_pol分别为材料波速和体积内能，μ_pol为动态黏度系数。S₁、S₂、S₃、γ_pol和α_pol为材料常数，V₀为相对体积。聚乙烯的模型参数如表2^[28]所示。

表  2  聚乙烯材料参数^[28]

Table  2.  Parameters for polyethylene^[28]

ρpol/(kg·m−3)	cpol/(m·s−1)	S1	S2	S3	γpol	αpol	V0
915	2901	1.481	0.0	0.0	1.64	0.0	1.0

下载: 导出CSV 

| 显示表格

LS-DYNA程序使用模型*MAT_JOHNSON_COOK描述铜管的变形特征，该模型中材料的强度函数与应变率和温度相关，其具体参数如表3^[28]所示。

表  3  铜材料参数^[28]

Table  3.  Parameters for copper^[28]

参数	值		参数	值		参数	值
密度ρcop/(kg·m−3)	8330		材料常数ncop	0.31		EOS常数S1	1.49
杨氏模量Ecop/GPa	110		材料常数mcop	1.09		EOS常数S2	0.0
泊松比μcop	0.35		材料常数wcop	0.025		EOS常数S3	0.0
熔化温度T/K	1357		材料常数Acop/GPa	0.089		EOS常数αcop	0.47
初始内能E0/(m·s−1)	0.0		材料常数Bcop/GPa	0.292		EOS常数Ccop	4430

下载: 导出CSV 

| 显示表格

模型中对炸药和空气分别通过模型*MAT_HIGH_EXPLOSIVE_BURN和*MAT_NULL进行表征，并结合相应的状态方程描述其压力与能量之间的关系，详细参数参考文献[27]。

2.2   模型验证

为更好地描述爆破裂纹扩展特征，分别分析沿径向和轴向爆炸压力的演化过程，如图9所示。径向：炸药起爆，产生高能冲击波，作用于孔壁，形成小范围压剪粉碎区。冲击波穿过聚乙烯等耦合介质后消耗许多能量，衰减转变为应力波，其应力峰值低于岩体的动态抗压强度，压剪粉碎区的范围不再扩展，然而岩石的动态抗拉强度相对较低，因此，在应力波作用下形成了拉剪破碎区。应力波向自由面传播的同时，径向裂纹沿裂隙尖端不断扩展，传播至自由面时发生反射，边界处在反射拉伸应力波作用下形成了裂纹剥落区。轴向：炸药在顶端起爆后，孔壁受到冲击波作用，形成V形损伤区。冲击波呈箭头状向下不断传播，V形损伤区的范围也随之不断扩大。冲击波传播至底端时。发生发射，与径向反射波的作用效果相同，在底部也会产生相应的裂纹剥落。

图  9  爆破荷载下岩石压力和损伤过程

Figure  9.  Explosion pressure and damage evolution of rock mass under blasting load

下载: 全尺寸图片 幻灯片

数值模拟结果与实验结果^[26]的比较如图10所示，图10(a)和(b)分别为爆破裂纹扩展模式及爆炸峰值压力衰减曲线。由图10(a)可知，随着到顶端起爆点距离的增大，爆破裂纹的分布更密集，且边界处的损伤剥落更明显。这是因为，炸药自顶端起爆后，轴向上的爆炸应力波与反射应力波叠加，致使爆炸应力场逐渐增强，进而导致岩体的损伤加剧。值得注意的是，实验结果中爆破裂纹分布更离散，这可能是由于实验中所用的岩样存在天然裂隙，这些缺陷会影响爆破裂纹扩展的连续性。受单元尺寸的限制，难以在有限元模型中重构含微小缺陷的岩体进行计算。但从裂纹的区域分布及扩展行为来看,数值模型可以定性地表征岩体爆破破裂现象。在定性验证数值模型的基础上，又进一步定量对比了岩体的爆炸压力衰减曲线，如图10(b)所示。由图10(b)可知，理论计算与数值模拟结果中爆炸峰值压力的变化趋势基本一致，且实验中测得的压力数据点也基本分布在压力衰减曲线上。综上所述，数值模拟与实验和理论分析的结果吻合较好，该模型参数适用于本文的研究工作。

图  10  岩样爆破实验^[26]与模拟结果比较

Figure  10.  Comparison of experimental^[26] and numerical results for rock sample blasting

下载: 全尺寸图片 幻灯片

3.   预裂爆破数值计算

基于验证的材料模型，进行预裂爆破数值计算，分析地应力对岩体预裂爆破裂纹扩展行为及损伤分布特征的影响。模拟方案分2组：静水压力组和非静水压力组，并以无地应力工况(N0)作为对照组，不同工况下应力加载条件如表4所示。

表  4  地应力加载条件

Table  4.  In-situ stress conditions used in the numerical study

应力 状态	工况	σhor/MPa	σv/MPa		应力 状态	工况	σhor/MPa	σv/MPa
静水 压力	H1	10	10		非静水 压力	A1	0	20
	H2	20	20			A2	10	20
	H3	30	30			A3	30	20
	H4	40	40			A4	40	20

下载: 导出CSV 

| 显示表格

为确保与理论模型中平面应变问题条件的一致性，预裂爆破采用准三维模型，即在厚度方向设置单层网格，并约束垂直平面方向的位移。预裂爆破数值模型如图11所示，模型的长和宽分别为1.0和0.6 m，炮孔布置在模型中心，两炮孔连线中点为坐标原点(0, 0)。炮孔之间的距离为0.6 m，炮孔半径为20 mm，采用不耦合装药，不耦合系数为2.0。模型中的单元类型为六面体八节点实体单元，单元的平均尺寸约为2 mm，单元总计约30万个。为保证数值计算结果具有较好的收敛性，模型使用流固耦合方法，其中岩体采用拉格朗日单元算法，炸药和空气采用欧拉单元算法，并通过多物质组和流固耦合等关键字实现爆炸荷载的有效传递^[29]。模型周边的单元设置为无反射边界条件用于模拟无限区域岩体，并在这些单元上施加初始压力模拟地应力作用，采用dynain文件法进行应力初始化，该方法可用于显/隐式求解计算，运算过程稳定且效率较高。

图  11  预裂爆破数值模型及炮孔附近局部网格

Figure  11.  Numerical model of presplit blasting and local mesh near the borehole

下载: 全尺寸图片 幻灯片

3.1   应力初始化

LS-PREPOST后处理软件中通过设置局部柱坐标系，可直观体现地应力下岩体中的静态应力分布特征，如图12所示。不同侧压系数下，径向应力和切向应力的大小以中点O均呈现左右对称分布，且在炮孔连线方向中点附近趋于平稳变化。值得注意的是，单向压力时，竖直方向上的径向与切向均会出现局部拉应力集中现象。这些现象与1.1节中的理论分析结果具有良好的一致性。

图  12  静态应力初始化结果

Figure  12.  Results of static stress initialization

下载: 全尺寸图片 幻灯片

为进一步验证理论计算结果的可靠性，在定性对比岩体静态应力分布特征的基础上，选取炮孔连线上的若干单元，监测其应力的变化特征，如图13所示。炮孔连线方向上，径向应力和切向应力的理论解与相应的数值解基本吻合，但理论计算的曲线相较于数值计算结果更光滑。这可能是因为，理论分析中假设材料线弹性变形，而数值模拟所用的RHT模型中假设材料弹塑性变形。此外，数值模型中单元数量有限，能够选取的测点较少，因此数值计算曲线与理论计算曲线存在微小差异。

图  13  不同方法得到的静态应力分布的比较

Figure  13.  Comparison of static stress distributions obtained by different methods

下载: 全尺寸图片 幻灯片

3.2   爆炸压力演化

预裂爆破中爆炸压力演化过程如图14所示。预裂孔同时起爆后，爆炸应力波以炮孔为中心向四周呈环状扩散。t=60 μs时，左右两侧预裂孔的爆炸应力波相遇，在到达岩体边界后应力波被吸收。t=70 μs时，爆炸应力波发生叠加并持续增强，致使叠加区域的爆炸压力和塑性应变显著增加。随后，应力波叠加效应弱化，爆炸能量不断耗散，动态压力持续衰减并逐渐趋于稳定状态。

图  14  爆炸压力演化过程

Figure  14.  Explosion pressure evolution process

下载: 全尺寸图片 幻灯片

中点O处的动态切向应力变化曲线如图15所示。随静水压力的提高，切向压应力和拉应力分别呈现升高和降低的趋势，且拉应力持续时间也逐渐减少，如图15(a)所示。静水压力在30 MPa前，其对应的切向拉应力峰值均高于岩石的抗拉强度，因此，孔间裂纹能够贯穿。而当静水压力增加至30 MPa时，切向拉应力峰值仅为5.2 MPa，低于岩石的抗拉强度，且拉应力的作用时间仅有35 μs，爆破裂纹难以在孔间贯通。不同侧压系数条件下，炮孔连线方向上中点O的切向应力变化曲线如图15(b)所示。同样的，侧压系数为0.5、1.5和2.0时，其对应的切向拉应力峰值分别为28.8、25.0和22.1 MPa，均高于岩石的抗拉强度，因此岩体的预裂爆破效果较好。而侧压系数为0时，即单向受压时，点O处的切向应力峰值仅有6.5 MPa，低于岩石的抗拉强度，此时的预裂爆破效果最差。当竖直方向上的初始压力固定，水平方向上初始压力不断升高时，岩体预裂爆破效果并不是随之变得越来越差。因此，非静水压力条件下，预裂爆破效果不仅与初始压力值相关，还受水平方向与竖直方向上的主应力差影响，爆破裂纹更倾向于沿着最大主应力方向扩展。实际岩体爆破工程中，建议预裂孔应尽可能沿最大主应力方向布置，从而更好地发挥地应力对爆破裂纹的诱导作用。

图  15  不同工况下中点O处动态切向应力的时程曲线

Figure  15.  Dynamical tangential stress-time curves at the middle point O under different conditions

下载: 全尺寸图片 幻灯片

3.3   裂纹扩展机制

地应力下岩体预裂爆破裂纹扩展模式如图16所示。由图16可知，爆破裂纹扩展模式大致可以分为两类，即孔间裂纹贯穿型（Type Ⅰ）和孔间裂纹非贯穿型（Type Ⅱ）。静水压力下，岩体爆破裂纹的数量随着静水压力的增加而减少，且损伤程度逐渐降低，破裂模式逐渐由Ⅰ型向Ⅱ型转变。当静水压力达到30 MPa时，爆破裂纹未能在孔间贯穿。非静水压力条件下，随着侧压系数的增加，爆破裂纹与损伤程度也呈现降低减缓的趋势，但其破裂模式却呈现出由Ⅱ型向Ⅰ型的转变趋势，其原因可通过裂纹扩展机制解释。总的而言，岩体破裂模式的变化规律与上述分析结论具有较好的一致性。

图  16  地应力下岩体预裂爆破裂纹扩展模式

Figure  16.  Fracture patterns in presplit blasting under in-situ stress

下载: 全尺寸图片 幻灯片

选取工况N0和H3的计算结果，分别代表Ⅰ型和Ⅱ型2种破裂模式。通过分析爆破裂纹的扩展行为及相应的质点位移变化规律，从而揭示岩体预裂爆破的损伤机制^[30]，如图17所示。Ⅰ型破裂过程如图17(a)所示，t=20 μs时，爆破裂纹在孔壁处呈辐射状向四周传播，在爆炸压力作用下岩体质点沿径向产生塑性变形。径向裂纹不断向外扩展，应力波在孔间相遇叠加，t=100 μs时，岩体质点在叠加应力波作用下沿切向产生拉应力，促使孔间裂纹形成。t=200 μs时，孔间裂纹与径向裂纹相遇凝聚，最终形成宏观破裂面。Ⅱ型破裂过程如图17(b)所示，爆炸初期，Ⅱ型爆破裂纹的扩展行为与Ⅰ型基本一致，但压碎区与裂隙区的范围相对较小。t=100 μs时，从位移矢量图可以看出，初始压力对孔间中点附近岩体质点的抑制作用。t=200 μs时，即使质点有朝相反方向的运动趋势，但位移变化量较小，并不能促使裂纹的形成，最终裂纹扩展趋于稳定状态。

图  17  爆破裂纹及位移矢量演化

Figure  17.  Evolution of blasting cracks and corresponding displacement vectors

下载: 全尺寸图片 幻灯片

3.4   损伤分布特征

为更好地理解岩石预裂爆破损伤特征，通过Hough变换法^[31]量化表征爆破裂纹的分布形态。首先，导出数值计算模型的损伤云图，然后通过图像处理软件对损伤云图进行二值化处理。最后，将二值化图像导入到Matlab程序中进行二次图像处理，便可获得最终的爆破裂纹分布形态。

爆破裂纹分布特征如图18所示。正东方向，即水平方向正半轴定义为0°方向，并规定逆时针为正方向。从裂纹分布玫瑰图可以看出，预裂爆破中，裂纹的角度分布具有中心对称性，且大部分角度的频率分布在1%～5%，这说明裂纹角度分布是相对平均的。此外，相较于水平方向，竖直方向上的裂纹占比较少，这是由于炮孔连线方向上应力波叠加效应导致的。

图  18  爆破裂纹的分布特征

Figure  18.  Distribution characteristics of blasting cracks

下载: 全尺寸图片 幻灯片

从裂纹的长度统计直方图来看，裂纹分布频率具有上下对称性，且在竖直方向上的分布频率最高。通过正弦函数拟合，发现裂纹的分布频率随着角度的增大变化较明显，在静水压力时，波峰和波谷分别出现在竖直方向和水平方向。值得注意的是，在非静水压力条件下，且当侧压系数k<1时，裂纹的分布频率随着角度的增大变化较缓，正弦函数拟合曲线会出现更多的峰值。这说明，爆破裂纹更倾向于沿着最大主应力方向扩展，压力差过大时会影响预裂爆破效果

4.   预裂爆破方案优化

4.1   方案概况

基于水电站中的边坡爆破工程，结合上述的理论分析结果，对其预裂爆破参数进行优化，优化前后的方案分别命名为方案1和2，其简化的几何模型如图19所示。模型的长和宽分别为6.4和3.6 m，炮孔排距为1.6 m。模型分为两部分，即岩体保留区和爆破开挖区。岩体保留区长1.6 m，宽3.6 m，预裂孔直径和装药直径分别为40和20 mm，相邻孔距为0.6 m。参考1.3节中的计算数据，对其装药直径进行优化。爆破开挖区的长和宽分别为4.8和3.6 m，生产孔的直径和装药直径保持不变，分别设定为40和30 mm，相邻炮孔间距为1.2 m。方案2与方案1中生产孔的起爆方式有所区别。方案1中，起爆预裂孔后，同时起爆生产孔。而方案2中，预裂孔起爆后，首先起爆临近自由面的生产孔，然后再起爆下一排生产孔。需要说明的是，该简化模型与理论模型的边界条件略有不同。其中，自由边界条件会对应力分布及传播产生影响，也可能会诱导开挖区的爆破裂纹扩展方向发生改变，但预裂孔距自由边界较远，几乎可以忽略自由面处透射边界条件对孔间成缝效果的影响。

图  19  岩石预裂爆破优化设计

Figure  19.  Optimal design scheme for rock presplit blasting

下载: 全尺寸图片 幻灯片

4.2   结果分析

不同方案下岩石爆炸压力演化过程如图20所示。首先起爆预裂孔，其产生的应力波以炮孔为中心向四周扩散，随后相遇并发生叠加，进而导致炮孔连线方向上的裂纹优先扩展。由于采用不耦合装药结构，爆炸压力相对较低，导致其他方向上爆破裂纹相对较少，预裂缝的形成可以阻止生产孔中的爆破能量向岩体保留区传递。预裂孔之间形成贯通裂隙面后，随即同时起爆2排生产孔。由于生产孔装药量较大，可以观察到周边岩体的损伤更严重，如图20(a)所示。除生产孔间爆炸应力波叠加外，还可以看到生产孔与预裂孔之间的应力波也会发生叠加。在优化方案中，岩石爆炸应力演化过程如图20(b)所示。参考1.3节中计算的理论数据，将预裂孔中的药量减小到建议值后，可以观察到预裂孔产生的爆炸压力显著降低，尽管爆炸能量有所减小，但仍可以看到孔间的损伤凝聚。需要说明的是，方案2与方案1中生产孔的起爆方式有所区别。为了尽可能利用生产孔中炸药产生的爆炸能量，首先起爆临近自由面一排的炮孔，由于爆炸应力波在自由面处发生反射，会促进拉伸应力波的产生，进一步加剧自由面附近岩体的层裂。此外，靠近自由面的炮孔起爆后，会在很大程度上降低内部岩体的夹滞作用，从而为下一轮岩体爆破提供更多的自由补偿空间。

图  20  不同方案下岩石爆炸压力演化过程

Figure  20.  Evolution of explosion pressure in different design schemes for rock blasting

下载: 全尺寸图片 幻灯片

将岩石的爆破损伤图像导出，通过边缘检测算法将损伤分区，进而获得岩石的块体尺寸分布特征，如图21所示。岩石在爆破荷载下产生了破碎区和裂隙发育区，裂隙发育区中的径向裂纹交错分布，形成了形状各异、大小不一的碎块。方案1中，预裂孔间的裂纹虽然可以贯通，但由于装药量大，导致岩体保留区的裂纹发育程度较高，如图21(a)所示。而优化方案中，将装药量降低到理论计算值后，预裂孔间的裂纹仍可以贯穿，且岩体保留区的损伤程度较低，如图21(b)所示。优化方案中不仅降低了对岩体保留区的扰动，还可以观察到爆破开挖区中的岩体破碎程度较高，尤其是靠近自由面一侧的岩体，破碎程度更高。由以上分析可知，理论计算给出的建议值具有较好的应用价值。

图  21  不同方案中岩石爆破块体分布特征

Figure  21.  Distribution characteristics of rock fragmentation induced by blasting in different schemes

下载: 全尺寸图片 幻灯片

5.   结　论

基于弹性力学平面应变问题假设，建立了地应力下岩体预裂爆破理论模型。通过Laplace变换和数值反演的方法分析了爆炸应力波的演化规律，讨论了地应力对岩体预裂爆破应力场分布的影响。此外，依据爆破实验验证了RHT岩石材料模型，采用验证的模型模拟了岩体预裂爆破过程，讨论了不同地应力条件对岩体预裂爆破损伤的影响，得到的主要结论如下。

(1)静水压力条件下，地应力对爆破产生的径向与切向压应力有增强作用，而对其拉应力有削弱作用，预裂孔间的质点可能会因为切向位移受限而无法形成贯通裂隙。

(2)非静水压力条件下，地应力在抑制爆破裂纹扩展范围的同时，对裂纹的扩展方向也有一定的导向作用，当炮孔中心连线方向与最大主应力方向平行时，可以获得较好的预裂爆破效果。

(3)基于应力波叠加破坏理论，提出了地应力下岩体预裂爆破孔间成缝的判别准则，得到了不同地应力条件下炮孔间距与不耦合系数之间的关系，结合预裂爆破参数优化得到了较好的验证。

(4)通过数值模拟中预裂爆破裂纹的扩展行为，验证了孔间拉伸破裂成缝的判别准则，并结合环向拉应力演化过程及位移矢量的变化规律揭示了岩体预裂爆破的成缝机制。