经济社会发展巨大的资源消耗使得地球浅部资源已趋于枯竭，目前矿产和能源资源开采已逐步进入1 000～2 000 m以深水平，未来深部资源开采将成为常态。以金属矿为例，目前全国金属矿千米深井约32座，如安徽铜陵的冬瓜山铜矿开采深度大于1 000 m，云南会泽的铅矿开采深度1 300 m，辽宁抚顺的红透山铜矿开采深度达到1 600 m^[1]。据统计，我国1/3以上地下金属矿山未来10年内开采深度将超过1 000 m，部分可到2 000～3 000 m^[2]。地应力实测结果表明，1 000～5 000 m开采深度地应力将达到50～135 MPa^[1]。钻孔爆破目前仍然是深部硬岩开挖的主要手段，深部岩体爆破破岩是爆炸荷载与高地应力共同作用的结果^[3]。

对于高地应力作用下的岩体爆破，已经开展了一些试验与数值模拟研究^[3-7]。研究结果表明，高地应力的存在改变了爆生裂纹的传播方向和扩展长度，裂纹首先呈辐射状从炮孔壁向外传播，随后逐渐平行于最大主应力方向向外扩展；地应力抑制了裂纹扩展的长度，地应力场水平越高，裂纹扩展长度越小。然而，目前相关的研究主要针对单孔爆破方面，而实际工程的爆破是多个炮孔同时起爆，需要考虑炮孔间的相互作用。在多孔爆破研究方面，He等^[8]通过数值模拟研究发现，双孔爆破引起的应力波叠加会在炮孔之间形成径向裂纹而使岩石断裂，并指出随着不耦合系数的增大，垂直于炮孔连线的裂纹的长度呈线性减小；Zhao等^[9]研究发现，增加相邻炮孔间微差起爆时间，可以在炮孔之间形成更好的裂纹网络，并使径向裂纹贯通；李洪伟等^[10]通过现场试验和数值分析发现，相邻炮孔的间距越大，主裂纹相互贯通的位置越靠近炮孔连线方向。上述研究虽考虑到了炮孔间的相互影响，但针对的是无地应力或者地应力水平很低的浅部岩体，有关高地应力条件下炮孔间的相互作用及破岩机理还缺乏深入系统的研究。因此，有必要更深入地研究高地应力作用下多孔爆破的破岩机制。

本文中，采用SPH（smoothed particle hydrodynamics）-FEM（finite element method）耦合数值模拟方法，研究不同地应力条件下双孔爆破岩体裂纹起裂、传播及贯通过程，分析高地应力作用下炮孔周围应力场动态调整过程、分布特征及对岩体开裂的影响，探讨高地应力岩体爆破的炮孔布置方式，以期研究成果可为深部岩体工程爆破设计优化提供理论依据。

1.   分析模型

岩体爆破是一个非常复杂的三维动力学过程，涉及到多个柱状炮孔内炸药的爆轰、岩体结构面和岩体各向异性特征等各种复杂问题。为了揭示高地应力对炮孔间裂纹传播及贯通过程的影响，需要对复杂的问题进行简化。本文中假定：（1）岩体为各向同性均质材料；（2）同一微差段内的炮孔同时起爆；（3）炸药爆轰波速无限大，柱状炮孔内的炸药同时起爆。这样，由于柱状炮孔的长度远大于炮孔直径，从而可以采用平面应变模型来研究这一复杂问题而不失一般性。建立如图1所示的双孔爆破计算模型，模型尺寸为8 m×8 m，2个炮孔位于模型的中部，炮孔直径为42 mm，炮孔间距为0.8 m。炮孔采用耦合装药结构，药卷直径等于炮孔直径。

图  1  双孔爆破分析模型

Figure  1.  The analysis model for double-hole blasting

下载: 全尺寸图片 幻灯片

炮孔周围的岩体赋存于地应力环境之中，水平向地应力为σ_x，竖直向地应力为σ_y。为了对比研究地应力大小和主应力方向对岩石爆破开裂过程的影响，本文中研究10种地应力工况，如表1所示，其中工况1为无地应力场σ_x=σ_y=0 MPa，工况2～5为不同水平的静水地应力场，工况6～10是侧压力因数λ为2、3、4的非静水地应力场。

表  1  地应力加载工况

Table  1.  In-situ stress conditions used in numerical calculations

应力场	工况	σx/MPa	σy/MPa	λ=σx/σy
无初始地应力场	1	0	0
静水地应力场	2	10	10	1
	3	20	20	1
	4	40	40	1
	5	80	80	1
非静水地应力场	6	40	20	2
	7	30	10	3
	8	60	20	3
	9	120	40	3
	10	80	20	4

下载: 导出CSV 

| 显示表格

2.   数值计算方法与模型

2.1   SPH-FEM耦合数值模拟方法

传统的有限元法(FEM)难以处理岩体爆破等大变形问题。光滑粒子流体力学方法(SPH)是一种无网格的拉格朗日数值计算方法，与传统的基于网格的数值计算方法相比，SPH方法省去了网格初始划分和重构的步骤，且粒子之间没有固定连接，不会出现网格畸变现象，适用于模拟大变形的问题。但是，SPH方法存在计算效率低、难以解决边界条件的问题。基于此，本文中采用SPH-FEM耦合数值计算方法来模拟岩体爆破，炮孔近区采用SPH方法模拟炸药爆轰和岩体破碎等极端变形问题，中远区采用FEM方法处理爆炸地震波传播及边界约束问题，以提高计算效率。在耦合界面处，SPH粒子和FEM网格以点-面胶结的方式来保证两者间的位移协调，如图2所示。为保障计算精度和提高计算效率，模型中心1.6 m×1.6 m范围内采用SPH模拟，其余范围采用FEM网格模拟，如图3所示，FEM网格最小尺寸为8 mm，最大尺寸为40 mm。有关SPH-FEM耦合数值模拟方法的理论原理详见文献[11-12]。

图  2  SPH-FEM耦合示意图

Figure  2.  Illustration of the coupled SPH-FEM algorithm

下载: 全尺寸图片 幻灯片

图  3  数值计算模型

Figure  3.  The numerical model used in the calculations

下载: 全尺寸图片 幻灯片

2.2   材料模型

Johnson等^[13]提出的JH-2本构模型描述了冲击载荷作用下脆性材料的力学行为，被广泛应用于岩石爆破数值模拟^{[3, 8]}。该材料模型主要包括3部分：多项式状态方程、强度模型和损伤模型。在JH-2模型中，多项式状态方程描述了流体静水压力p与体积应变μ之间的关系：

式中：$\mu = {\rho }/{{\rho _0}} - 1$，ρ为材料的当前密度，ρ₀为材料的初始密度；K₁、K₂和K₃均为材料常数；Δp为流体静水压力的增量。

JH-2强度模型同时考虑了完整材料的强度和断裂状态下的残余强度，材料从完整状态到断裂状态的变化过程由损伤变量D控制。

当材料未发生损伤时（D=0），归一化的材料强度为：

当材料完全开裂时（D=1），归一化的材料残余强度为：

在材料从初始状态向完全开裂状态转化过程中（0＜D＜1），归一化的材料损伤强度为：

式中：A、B、C、M和N均为材料常数；上标星号*表示归一化的量；${p^*}$为归一化的静水压力；${T^*}$为归一化的抗拉强度；$\sigma _{\rm{i}}^*$和$\sigma _{\rm{f}}^*$分别为材料在完整状态和完全开裂状态下的归一化强度；${\sigma ^*}$为归一化的等效应力；${\dot \varepsilon ^*}$为归一化的等效应变率。以上归一化的量分别采用下式计算：

式中：p为实际静水压力；T为材料所能承受的实际最大静水拉应力；σ_i和σ_f分别为材料在完整状态和完全开裂状态下的实际强度；σ为实际等效应力；${\sigma _{{\rm{HEL}}}}$为材料的Hugoniot弹性极限；${p_{{\rm{HEL}}}}$为Hugoniot弹性极限的静水压应力；$\dot \varepsilon$为实际等效应变率；${\dot \varepsilon _0}$为参考应变率，取${\dot \varepsilon _0}$=1.0 s⁻¹。

损伤模型中，损伤变量D描述材料从完整强度到剩余强度的转变，其表达式为：

式中：$\Delta {\varepsilon _{\rm{p}}}$为一个积分循环中形成的等效塑性应变增量；$\varepsilon _{\rm{p}}^{\rm{f}}$是在应力作用下材料发生损伤的等效塑性应变。$\varepsilon _{\rm{p}}^{\rm{f}}$的表达式为：

式中：D₁和D₂均为损伤常数。

根据Banadaki等^[14]的研究，本文中所取的岩石材料主要物理力学参数如表2所示。

表  2  岩石材料物理力学参数

Table  2.  Physical and mechanical parameters of the rock material

ρ/(g·cm-3)	K/GPa	G/GPa	σHEL/GPa	A	N	C	B	M	D1	D2
2.66	25.7	21.9	4.5	0.76	0.62	0.005	0.25	0.62	0.005	0.7

下载: 导出CSV 

| 显示表格

采用Jones-Wilkins-Lee (JWL)状态方程^[15]模拟炸药爆轰过程，JWL状态方程描述了炸药爆轰产物压力、体积和能量之间的关系：

式中：p_d为爆轰产物的压力；E₀为爆轰产物初始体积内能；V为爆轰产物的相对体积；A₁、B₁、R₁、R₂、ω为炸药常数。参考文献[14]，有关的炸药参数取值为：初始爆轰压力p_d0=16.0 GPa，E₀=7.38 GJ/m³，密度ρ_e=1 320 kg/m³，爆轰速度c_d=6 690 m/s，A₁=586 GPa，B₁=21.6 GPa，R₁=5.81，R₂=1.17，ω=0.28。计算得到的炮孔壁上的爆炸荷载压力时程曲线如图4所示，荷载峰值为3.16 GPa。

图  4  爆炸荷载压力时程曲线

Figure  4.  Blasting pressure varying with time

下载: 全尺寸图片 幻灯片

2.3   数值计算模型验证

无地应力条件下(σ_x=σ_y=0 MPa)单个炮孔爆破时，数值模拟计算得到的炮孔周围岩体裂纹分布如图5(a)所示。本文计算条件下，炮孔周围形成了粉碎区和破碎区，粉碎区半径约等于5倍炮孔半径，破碎区半径约为19倍炮孔半径。Banadaki等^[14]针对上述物理力学性质的花岗岩开展了室内爆破模型试验，并采用染料浸染和高强度紫外线下数字摄影的方式记录了炮孔周围的裂纹分布，如图5(b)所示。该试验中，粉碎区半径为5～6倍炮孔半径，破碎区半径为20～22倍炮孔半径。对比图5(a)和5(b)，本文数值模拟的岩体爆破开裂模式与Banadaki等^[14]的室内试验结果基本一致，这表明本文中采用的数值模型和材料参数是合理的。此外，Banadaki等^[14]也采用JH-2模型、JWL状态方程和上述材料参数对室内爆破模型试验进行了数值模拟，验证了该模型和参数选择的合理性。

图  5  本文数值模拟与Banadaki等^[14]的试验结果的对比

Figure  5.  Comparison between the numerical simulation and the experimental result by Banadaki, et al^[14]

下载: 全尺寸图片 幻灯片

3.   地应力对岩石爆破开裂的影响

3.1   静水地应力场条件下的岩石爆破开裂

不同静水地应力水平下，炮孔内炸药起爆后的裂纹传播过程如图6所示。在0.05 ms内，不同地应力水平下裂纹扩展速度基本一致。这是由于在炸药起爆之初，爆炸荷载远高于地应力，地应力对爆炸致裂过程几乎没有影响。随着爆炸荷载压力衰减，地应力对岩石爆破开裂的影响得以体现。在起爆后0.2 ms，σ_x=σ_y=0 MPa条件下2个炮孔产生的裂纹已经连接贯通，而此刻其他地应力水平下炮孔间的裂纹并无贯通。在起爆后0.5 ms，所有工况下的裂纹扩展过程均已停止，裂纹长度不再增加，如图7所示。最后，σ_x=σ_y=0 MPa，σ_x=σ_y=10 MPa和σ_x=σ_y=20 MPa等3种地应力水平下炮孔间的裂纹能完全贯通，形成新的自由面；而σ_x=σ_y=40 MPa和σ_x=σ_y=80 MPa地应力水平下炮孔间的裂纹不能贯通，达不到形成良好的爆破开挖面效果。由此可见，地应力对岩石爆炸致裂起抑制作用，随着地应力水平的提高，裂纹扩展长度减小。因此，对于高地应力岩体爆破，为了形成良好的爆破开挖面，宜适当减小炮孔间的间距。高地应力条件下爆生裂纹的扩展长度除了与地应力水平有关外，还与岩石性质、岩体结构面、炸药种类、装药结构、起爆方式等众多因素有关，现场爆破过程中炮孔间距的减小量需要通过现场试验确定。

图  6  不同静水地应力水平下的岩石爆破开裂过程

Figure  6.  Blast-induced rock fracture processes under different hydrostatic in-situ stress levels

下载: 全尺寸图片 幻灯片

图  7  不同静水地应力水平下裂纹扩展长度随时间的变化

Figure  7.  Variation of crack length with time under different hydrostatic in-situ stress levels

下载: 全尺寸图片 幻灯片

对于高地应力岩体爆破，由于地应力对爆炸致裂的抑制作用，若不减小炮孔间距而又要实现炮孔间裂纹的连接贯通，势必要采用高密度、高爆轰波速的炸药来提高爆炸荷载压力。对于σ_x=σ_y=40 MPa和炮孔间距S=0.8 m的计算工况，当炮孔壁上的爆炸荷载压力提高至3.58 GPa时，炮孔间的裂纹得以连接贯通，相比于图6(d)，炮孔壁上的爆炸荷载压力提高了13%；而当地应力达到80 MPa时，炮孔壁上的爆炸荷载压力提高至4.12 GPa才可实现炮孔间裂纹的贯通，如图8所示，相比于图6(e)，炮孔壁上的爆炸荷载压力提高了30%。爆炸荷载压力提高后尽管实现了爆孔间裂纹的贯通，但不可避免地会对保留岩体造成更大的损伤破坏。对于σ_x=σ_y=80 MPa，爆炸荷载压力未提高时，在垂直于炮孔连线方向上裂纹最大扩展长度为0.17 m，0～0.1 m范围内的裂纹密度(定义为单位面积内裂纹的长度^[16])为90 m⁻¹，＞0.1～0.3 m范围内的裂纹密度为10 m⁻¹，如图6(e)所示。当爆炸荷载压力提高30%后，该方向上的裂纹延伸到了0.26 m，0～0.1 m范围内的裂纹密度基本不变，为100 m⁻¹，而＞0.1～0.3 m范围内的裂纹密度显著增大，变为26 m⁻¹，如图8所示。可见，爆炸荷载压力提高30%后，保留岩体的损伤深度和0.1 m以外岩体的损伤程度都显著增大。因此，对于高地应力岩体爆破，通过提高爆炸荷载压力来实现炮孔间的裂纹贯通不可取，宜通过减小炮孔间距来实现形成良好的爆破开挖面，同时达到保护保留岩体的目的。

图  8  地应力σ_x=σ_y=80 MPa条件下炮孔间裂纹贯通及相应的爆炸载荷压力随时间的变化

Figure  8.  The crack connection between the blastholes under σ_x=σ_y=80 MPa and the corresponding blasting pressure varying with time

下载: 全尺寸图片 幻灯片

3.2   非静水地应力场条件下的岩石爆破开裂

非静水地应力水平下(λ=3)，炮孔内炸药起爆后的裂纹传播过程如图9所示。在0.05 ms内，不同非静水地应力水平下的裂纹扩展速度基本一致，裂纹呈辐射状向四周传播。在起爆后0.2 ms，裂纹扩展呈现方向性，岩体开裂范围呈椭圆形分布，椭圆的长轴与水平地应力（最大主应力）方向基本平行。这是因为在非静水地应力情况下，最小主应力方向区域内的环向压缩应力大于最大主应力方向区域内环向压缩应力，根据前面的分析可知，压应力越高，对爆生裂纹的抑制作用越强，因此裂纹在最小主应力方向上传播受到的地应力抑制作用更大，裂纹主要沿最大主应力方向传播。图10给出了竖直向地应力σ_y=20 MPa，侧压力因数λ分别为2、3和4条件下的爆生裂纹最终分布。λ=2时，椭圆形裂纹分布的长轴a和短轴b之比a/b=2.0；λ=3时，a/b=2.2；λ=4时，a/b=3.4。可见，随着侧压力因数λ的增加，爆生裂纹扩展表现出的方向性越来越显著。

图  9  不同非静水地应力水平下的岩石爆破开裂过程(λ=3)

Figure  9.  Blast-induced rock fracture processes under different non-hydrostatic in-situ stress levels for λ=3

下载: 全尺寸图片 幻灯片

图  10  不同侧压力因数下的爆生裂纹分布(σ_y=20 MPa)

Figure  10.  Distributions of blast-induced cracks under different lateral pressure coefficients for σ_y=20 MPa

下载: 全尺寸图片 幻灯片

σ_x=3σ_y=30 MPa和σ_x=3σ_y=60 MPa 2种地应力工况下炮孔间的裂纹能完全贯通，形成新的自由面，而σ_x=3σ_y=120 MPa地应力水平下，炮孔间的裂纹不能贯通，达不到良好的破岩效果。对于σ_x=3σ_y=60 MPa的地应力，若炮孔不沿最大主应力方向布置，而沿最小主应力方向布置（炮孔连线方向为竖直向），如图11所示，此时，炮孔间的裂纹不能贯通，达不到良好的破岩效果，而且在垂直于炮孔连线方向上，裂纹扩展的长度增大到了0.3 m，较炮孔沿最大主应力方向布置的情况增大了76%。可见，对于高地应力岩体爆破，在现场条件允许的前提下，炮孔宜沿最大主应力方向布置，以达到良好的裂纹贯通和破岩效果。

图  11  炮孔沿最小主应力方向布置时的爆生裂纹分布(σ_x=3σ_y=60 MPa)

Figure  11.  Distribution of blast-induced cracks for the blasthole arrangement along the minimum principal stress direction for σ_x=3σ_y=60 MPa

下载: 全尺寸图片 幻灯片

4.   炮孔周围应力场动态演化过程

为进一步揭示高地应力岩体多孔爆破的破岩力学机制，下文分析爆破过程中炮孔周围应力场的动态演化过程。选取炮孔连线中间位置的E点和垂直于炮孔连线方向上的F点进行分析，E点和F点距炮孔1的中心均为0.4 m，如图12所示。由于岩体的抗拉强度远低于其抗压强度，爆炸荷载作用下炮孔周围岩体的开裂主要是环向拉应力导致的径向裂纹开裂，因此，下文主要分析环向应力的动态调整过程。为了对比研究双孔爆破时炮孔间应力波相互作用对岩体开裂的影响，本节先分析单孔爆破时炮孔周围的动态应力场。图13给出了无地应力工况下炮孔1爆破时E、F两点的环向应力时程曲线，图中应力以受拉为正，受压为负（下同）。可以看到，单孔爆破时，由于E、F两点中心对称，其应力时程曲线完全重合。由于爆炸荷载的冲击特性，围岩环向应力在爆炸荷载作用的初始阶段为压应力，随之由于岩体的径向压缩转变为拉应力，在距炮孔中心0.4 m处该拉应力可达62 MPa，远高于岩体抗拉强度，岩体发生张拉破坏。图14给出了炮孔1和2同时起爆时E、F两点的环向应力时程曲线。对比图13可以看到，双孔爆破时，由于应力波的叠加，E点环向拉应力得到了加强，达到了98 MPa；而F点由于受到炮孔2产生的径向压缩波的作用，该位置的环向拉应力减弱，降低为43 MPa。可见，多孔爆破时，炮孔间应力波的相互作用促使裂纹优先向炮孔连线方向上传播。

图  12  应力观测点的布置

Figure  12.  Arrangement of the stress observation points

下载: 全尺寸图片 幻灯片

图  13  单孔爆破时的环向应力变化曲线(σ_x=σ_y=0 MPa)

Figure  13.  Circumferential stress histories under single-hole blasting for σ_x=σ_y=0 MPa

下载: 全尺寸图片 幻灯片

图  14  双孔爆破时的环向应力变化曲线(σ_x=σ_y=0 MPa)

Figure  14.  Circumferential stress histories under double-hole blasting for σ_x=σ_y=0 MPa

下载: 全尺寸图片 幻灯片

不同静水地应力水平下，双孔爆破时测点E和F处的环向应力变化过程如图15所示。随着地应力水平的提高，E点处的环向拉应力逐渐降低，且拉应力的作用时间也越来越短。σ_x=σ_y=10 MPa和σ_x=σ_y=20 MPa两种地应力水平下E点的环向拉应力峰值分别为89 MPa和82 MPa，均高于岩石的动态抗拉强度，炮孔间的裂纹能够贯通。当静水地应力场达到80 MPa时，E点处的环向拉应力峰值仅为4 MPa，低于岩石的动态抗拉强度，且拉应力的作用时间仅有0.01 ms，裂纹不能在炮孔间贯通。由此可见，炮孔间的环向拉应力随着地应力水平的提高而降低，地应力对爆炸荷载的环向拉伸效应起抑制作用。由于受多个炮孔爆炸应力波相互作用的影响，相比于炮孔连线上的E点，垂直于炮孔连线方向上F点的环向拉应力更低，作用时间也更短。这表明，静水地应力条件下多个炮孔同时起爆时，垂直于炮孔连线方向传播的裂纹更易受到地应力的抑制而停止扩展。

图  15  不同静水地应力水平下炮孔周围环向应力随时间的变化曲线

Figure  15.  Circumferential stress histories around the blasthole under different hydrostatic in-situ stress levels

下载: 全尺寸图片 幻灯片

不同非静水地应力水平下(λ=3)，观测点E和F处的环向应力变化过程如图16所示。可以看到，当σ_x=3σ_y=30 MPa时，最小主应力方向上F点的环向拉应力峰值仅为15 MPa；σ_x=3σ_y=60 MPa和σ_x=3σ_y=120 MPa两种地应力工况下，F点甚至没有出现拉应力。而在炮孔连线方向（最大主应力方向）上的E点，3种非静水地应力工况下的环向拉应力峰值分别为94、65和24 MPa，远高于相同距离处F点的环向拉应力。因此，非静水地应力条件下，裂纹优先向最大主地应力方向扩展。

图  16  不同非静水地应力条件下炮孔周围环向应力随时间的变化曲线

Figure  16.  Circumferential stress histories around the blasthole under different non-hydrostatic in-situ stress levels

下载: 全尺寸图片 幻灯片

5.   结　论

通过上述数值模拟分析，可以得到以下结论：

（1）地应力对岩石爆炸致裂起抑制作用，对于高地应力岩体爆破，宜减小炮孔间的间距，以实现炮孔间裂纹贯通，达到形成良好的爆破开挖面效果。

（2）非静水地应力场条件下，爆破产生的裂纹主要沿最大主地应力方向扩展，沿最大主地应力方向布置炮孔，有利于炮孔间裂纹的连接贯通。

（3）多孔爆破时，爆炸应力波的相互叠加使炮孔间部位的环向拉应力得到加强，而垂直于炮孔连线方向上的环向拉应力受到减弱，因而静水地应力条件下垂直于炮孔连线方向传播的爆生裂纹更易受到地应力的抑制。