不耦合装药结构是爆破开采过程中常用的施工手段。大量实验和工程实际表明^[1-3]，选取合适的不耦合系数能够提升炸药能量的利用率和岩石的破碎效率。因此，研究装药不耦合系数与岩石爆破损伤特性的关系，对于指导爆破施工设计等工作具有重要的现实意义。

常见的爆破破岩理论^[4]认为，炸药在岩体内部的爆炸过程可分为爆炸应力波预裂破岩的动态冲击作用和爆生气体膨胀扩腔的准静态作用。叶志伟等^[5]通过轮廓爆破方法，研究了不耦合装药结构下爆破孔壁压力的变化规律，实验结果表明，薄壁钢管内爆生气体的准静态压力峰值随不耦合系数的增大近似呈线性增长。徐颖等^[6]利用不耦合装药结构对爆生裂纹数量进行了有效控制，并通过实验获取了具有最佳破岩效应的不耦合系数。费鸿禄等^[7]较早地论证了在不耦合装药结构的光面爆破过程中爆炸冲击波压力和爆生气体压力的计算方法。杨仁树等^[8]通过实验证明了不耦合偏心装药的有机玻璃爆破损伤的分形效应，建立了损伤变量与分形维数之间的预测模型。Yuan等^[9]利用颗粒流程序（particle flow code，PFC），对不同耦合系数下的岩体爆破损伤过程进行了有效模拟，论证了爆炸冲击波和爆生气体对岩体爆破损伤的破坏模式。

本研究根据已有理论，在考虑爆炸冲击波和爆生气体共同作用的基础上，提出一种“动”、“静”荷载混合施加的PFC爆破模拟方法，并以此进行花岗岩爆破数值模拟研究，建立花岗岩爆破损伤程度与装药不耦合系数之间的预测模型，以期为指导爆破施工设计提供理论支持。

1.   理论分析与计算

1.1   炮孔压力计算

耦合装药结构下，球形药包在岩体内部爆炸时，炮孔压力${p_{{\text{m }}}}$为^[10]

式中：${\,\rho _{\text{c}}}$为炸药密度，${v_{\text{c}}}$为炸药爆速，$\kappa$为绝热指数，${\,\rho _{\text{r}}}$为岩石密度，${v_{\text{r}}}$为岩石纵波波速。

空气不耦合装药结构下，设炮孔半径为${r_{\text{p}}}$，装药半径为${r_{\text{t}}}$，则装药不耦合系数$\zeta = {r_{\text{p}}}/{r_{\text{t}}}$。由于存在空气间隔，因此爆轰压力${p_{\text{t}}}$到达孔壁时为

式中：$\eta$为压力增大倍数，一般取8～11，本研究取8。

考虑理想情况下岩石与炸药的阻抗匹配时，有${\,\rho _{\text{r}}}{v_{\text{r}}}/\left( {{\rho _{\text{r}}}{v_{\text{r}}} + {\rho _{\text{c}}}{v_{\text{c}}}} \right){\text{ = 1/2}}$；将爆生气体视为等熵膨胀，且取$\kappa$=3。此时，耦合装药下，炮孔压力为爆轰波阵面压力；不耦合装药结构下，由于爆炸冲击波在炮孔内的多次反射，压力峰值发生变化，则爆炸冲击压力p_d可表示为

爆生气体压力p_g的计算公式^[11]为

式中：${p_\kappa }$为临界压力，${p_\kappa }{\text{ = }}200\;{\rm{MPa}}$。

由式(3)和式(4)可知，装药不耦合系数 $\zeta$ 是影响爆炸冲击压力和爆生气体压力的主要因素。

1.2   爆生气体扩腔作用起始时间的确定

不耦合装药结构下，爆生气体扩腔作用的起始时间是关系到爆破参量计算结果是否符合实际的重要因素。爆炸冲击波和爆生气体对炮孔壁产生作用的时间间隔Δt可采用如下近似计算方法^[9]

式中：${D_{\text{s}}}$为粉碎区半径；${v_{{\text{cp}}}}$为爆炸冲击作用下的裂纹扩展速度，可由分离式霍普金森压杆（split Hopkinson pressure bar，SHPB）实验获得。

有研究认为，爆生气体对炮孔的准静态作用完全发生在爆炸冲击波的动态作用之后。然而，实际上，随着炸药的不断消耗，高温高压的爆生产物不断累积，爆生气体的扩腔作用随着炸药的爆轰同步进行。为符合实际爆破过程，本研究将爆生气体扩腔作用的起始时间设置为炸药完全消耗时刻，即爆炸冲击压力达到峰值时开始对炮孔壁施加爆生气体的准静态作用。

1.3   PFC爆破“动”、“静”荷载混合施加原理

为进一步研究不同耦合系数下岩体爆破损伤特性，利用PFC对岩体爆破过程进行模拟分析，并提出爆破“动”、“静”荷载混合施加方法来模拟爆炸冲击波的动态作用和爆生气体的准静态作用。

颗粒膨胀加载法^[12-13]是采用PFC模拟爆破过程的主要方法之一。如图1所示，在炮孔中心位置设置动态荷载施加颗粒和准静态荷载施加颗粒。假定初始动、静荷载施加颗粒半径均为${r_0}$，随着荷载施加，颗粒逐渐膨胀，当爆轰压力为p时，荷载施加颗粒半径增量为

图  1  PFC中动静荷载共同施加原理

Figure  1.  Sketch of dynamic and quasi-static loads application in PFC

下载: 全尺寸图片 幻灯片

式中：${d_{\text{r}}}$为荷载施加颗粒半径增量，${S_{\rm n}}$为颗粒间法向接触刚度。

根据装药结构的不同，可由式(3)计算出爆炸冲击压力峰值p_dmax，则动态荷载颗粒对爆破空腔壁产生的最大径向接触力之和为

式中：F_dmax为爆炸冲击的动荷载作用下岩石颗粒受到的最大推力，$S_{\text{n}}^{\rm d}$为动荷载颗粒刚度，$d_{\text{r}}^{\rm d}$为动荷载颗粒半径增量。

同理，根据式(4)可计算出不同不耦合系数下的爆生气体峰值压力p_gmax，则静态荷载颗粒对爆破空腔产生的最大径向接触力为

式中：F_gmax为爆炸冲击的静荷载作用下岩石颗粒受到的最大推力，$S_{\text{n}}^{\rm g}$为静荷载颗粒刚度，$d_{\text{r}}^{\rm g}$为静荷载颗粒半径增量。

由式(6)可知，荷载施加颗粒的法向接触刚度${S_{\text{n}}}$是影响动静荷载大小的关键参量。由于两种颗粒始终处于重叠状态，如果同时将它们的刚度设定为较大值，则会导致荷载施加颗粒之间的接触力过大，进而影响通过“监测圆”获取的炮孔压力大小。为了解决这一问题，在动态冲击阶段将静载施加颗粒的刚度$S_{\text{n}}^{\text{g}}$和半径增量$d_{\text{r}}^{\text{g}}$设置为很低的数值（本研究取10⁻¹⁰），从而无法对动态冲击阶段的炮孔荷载产生影响。爆炸冲击压力达到峰值p_dmax后，恢复$S_{\text{n}}^{\text{g}}$和$d_{\text{r}}^{\text{g}}$至正常增量大小，此时动载压力开始降低，静载压力逐渐上升。当岩体内部动载残存应力$\sigma$_d与静载压力$\sigma$_g相等时，爆生气体的扩腔运动开始，经峰值p_gmax后，爆生气体压力缓慢下降至零。这也与实际炸药爆炸过程中爆炸冲击压力和爆生气体对孔壁的作用过程基本相同。

现有研究中多将爆炸应力波简化为三角波^[14]或简谐波^[15]。由于三角波荷载的变化趋势单一，难以表达爆炸冲击作用后爆生气体扩腔作用对炮孔压力的影响，因此，为了符合实际爆破过程，本研究采用半正弦波作为爆炸应力波形，实现爆破荷载的施加。爆生气体的准静态作用对孔壁荷载的变化较复杂。张玉磊等^[16]指出，爆生气体的膨胀做功过程与正弦函数的半周期变化相似。为此，本研究采用半正弦波作为爆生气体作用下炮孔压力的变化曲线，然而其半正弦周期比爆炸应力波的衰减时间更长。图2为本研究采用的爆炸冲击波与爆生气体共同作用下的炮孔压力时程曲线。

图  2  炮孔压力时程曲线

Figure  2.  Curve of borehole pressure-time history

下载: 全尺寸图片 幻灯片

2.   数值模拟

2.1   模型建立与细观参数标定

如图3所示，在长20 m、宽20 m的范围内共生成颗粒76445个，颗粒半径在0.0280～0.0465 m之间服从均匀分布，将模型伺服至稳定状态。为模拟无限岩体中的爆破过程，需在岩体周边设置无反射边界^[17]，以吸收爆炸应力波在边界处的反射拉伸作用。然而，在PFC中无法直接设置无反射边界，为此选取模型边界0.1 m范围内的颗粒作为无反射边界，在PFC中通过自编FISH语言，时刻监测爆破过程中边界颗粒的应力变化，施加与其大小相同、方向相反的外力，使其受力始终处于恒定状态，以此达到设置无反射边界的目的。

图  3  PFC岩体数值模型

Figure  3.  Numerical model of rock mass in PFC

下载: 全尺寸图片 幻灯片

颗粒间的黏结采用PFC中模拟岩土类材料常用的平行黏结模型（parallel bond model，PBM）^[18]，其最大的优点是既可以承受力又可以承受力矩。PBM由黏结元件和弹簧元件组成，黏结状态（bonded）下，颗粒间的接触具有一定的法向刚度和切向刚度。黏结键强度符合莫尔-库仑准则，当外界作用力超过其黏结强度时，黏结键发生破坏，颗粒间的接触变为无黏结状态（unbonded），此时颗粒间的接触退化为完全弹性接触。PBM黏结元件的基本结构如图4所示。图4中，g_s为黏结元件法向累积相对位移，${\bar k_{\text{n}}}$为线性接触法向刚度，${\bar k_{\text{s}}}$为线性接触切向刚度，${\bar K_{\text{n}}}$为平行黏结法向刚度，${\bar K_{\text{s}}}$为平行黏结切向刚度，${\,\mu _{\text{c}}}$为摩擦系数，${\bar \sigma _{\text{n}}}$为法向黏结强度，${\bar \sigma _{\text{s}}}$为切向黏结强度，${\varphi _{\text{c}}}$为内摩擦角。

图  4  平行黏结模型的基本结构

Figure  4.  Structure of parallel bond model

下载: 全尺寸图片 幻灯片

细观参数的合理选取是保证PFC数值模拟结果符合实际的重要前提。以弹性模量、单轴抗压强度和泊松比作为标准，通过PFC静力学模拟实验标定数值模型的细观参数，使数值模型与真实岩石材料具有相近的宏观力学性质。在只关注岩石破坏模式和损伤状态的前提下，以该套细观参数进行动力学实验或大变形过程模拟同样具有一定的合理性^{[9, 12, 18]}。

张学朋等^[19]的研究表明，当PFC数值模拟的试样颗粒数较少时，模型无法获得稳定的力学性质。为了避免颗粒数目对标定结果产生影响，将待试标定样大小增加至宽4 m、长8 m，共生成6067个颗粒，颗粒体系与上述岩石介质的数值模型相同。岩石材质选取花岗岩，根据室内实验研究结果^[20]，采用“试错法”^[21]在PFC中进行细观参数标定工作，最终确定了花岗岩试样的应力-应变曲线和破坏模式，如图5所示。试样的主要细观参数如表1所示，其中：R_min为最小颗粒半径，R_max为最大颗粒半径，R_max/R_min为粒径比，${\bar E_{\rm l}}$为有效模量，${\bar E_{\rm p}}$为线性接触模量，${\bar k_{\text{n}}}$/${\bar k_{\text{s}}}$为线性接触刚度比，${\bar K_{\text{n}}}$/${\bar K_{\text{s}}}$为平行黏结刚度比。表2列出了数值模拟和室内实验得到的试样抗拉强度${\sigma _{\text{t}}}$、弹性模量E_t和泊松比${v_{\text{t}}}$。

图  5  数值模拟结果与室内实验^[20]结果对比

Figure  5.  Comparison between numerical and experimental results^[20]

下载: 全尺寸图片 幻灯片

表  1  主要细观参数

Table  1.  Main microscopic parameters

Linear contact parameters
Rmin/m	Rmax/Rmin	${\bar E{_{\rm l}} }$/GPa	${\bar k{_{\text{n} }} }$/${\bar k{_{\text{s}} } }$	${\mu {_{\text{c} }} }$
0.028	1.66	45	2.5	1.0
Parallel bond parameters
${\bar \sigma {_{\text{n} }} }$/MPa	${\bar \sigma {_{\text{s} } }}$/MPa	${\bar E{_{\rm p}} }$/GPa	${\bar K{_{\text{n} } } }$/${\bar K{_{\text{s} } } }$	${\varphi {_{\text{c} } } }$/(°)
72	150	45	2.5	30

下载: 导出CSV 

| 显示表格

表  2  试样的主要力学参数

Table  2.  Main mechanical parameters of the sample

Method	${E{_{\text{t} } } }$/GPa	${\sigma {_{\text{t} } } }$/MPa	${v{_{\text{t} } } }$	$\,\rho$/(kg·m−3)
Experiment	21.4 – 63.7	137.7 – 163.8	0.23 – 0.26	2650 – 2740
Numerical simulation	27.2	150.93	0.24	2700

下载: 导出CSV 

| 显示表格

2.2   不同不耦合系数下炮孔压力变化

根据Yuan等^[9]的研究，炸药密度取${\,\rho _{\text{c}}}$= 1000 kg/m³，爆速${v_{\text{c}}}$= 2700 m/s。数值模拟中共取6种不耦合系数，不耦合系数与炮孔半径的对应关系见表3，其中：$\zeta$=1.0对应耦合装药结构，其余为不耦合装药结构。由式(3)、式(4)和式(5)，可计算出不同不耦合系数下的爆炸冲击压力和爆生气体压力，分别如图6和图7所示。

表  3  不同不耦合系数对应的炮孔半径

Table  3.  Blasthole radius vs. decoupling coefficient

$\zeta$	rp		$\zeta$	rp
1.0	0.060		1.6	0.096
1.2	0.072		1.8	0.108
1.4	0.084		2.0	0.120

下载: 导出CSV 

| 显示表格

图  6  爆炸冲击应力对孔壁的压力随不耦合系数的变化

Figure  6.  Curve of the pressure on blasthole wall under explosion stress with different decoupling coefficients

下载: 全尺寸图片 幻灯片

图  7  爆生气体对孔壁的压力随不耦合系数的变化

Figure  7.  Curve of the pressure on blasthole wall under detonation gas stress with different decoupling coefficients

下载: 全尺寸图片 幻灯片

根据计算结果，利用“动”、“静”荷载混合施加方法，在PFC中进行不同不耦合系数下花岗岩爆破模拟实验，得到的炮孔压力时程曲线见图8。

图  8  炮孔压力时程曲线

Figure  8.  Curves of blasthole pressure-time history

下载: 全尺寸图片 幻灯片

3.   计算结果分析

3.1   爆生裂纹扩展模式

为研究爆炸冲击波的“动”荷载和爆生气体的“静”荷载对爆生裂纹扩展模式的影响，以$\zeta$=1.4时的爆生裂纹扩展过程为例进行说明。图9显示了$\zeta$=1.4时的炮孔压力时程曲线和爆生裂纹扩展过程。

图  9  $\zeta$=1.4时炮孔压力曲线

Figure  9.  Curve of blasthole pressure ($\zeta$=1.4)

下载: 全尺寸图片 幻灯片

对于不耦合装药，孔壁压力在10 μs左右开始施加。58 μs左右炮孔压力达到峰值，如图10(a)所示，此时炮孔周围岩石介质的损伤程度高，形成粉碎区^[4]。至87 μs，为裂纹扩展阶段，由图10(b)可知，粉碎区外的径向裂纹在各个方向均有发育。由炮孔压力时程曲线可知，87 μs时爆生气体的扩腔作用开始施加，至100 μs时，在爆炸冲击波的动态预裂破岩作用下，爆生气体对径向裂纹产生驱裂作用，如图10(c)所示。爆生气体对裂纹的驱裂作用主要体现在对径向主裂纹尖端的延展作用，而对岩石的粉碎没有过多的作用。图10(d)显示了200 μs时爆生裂纹的最终分布情况。可以看出，100～200 μs区间爆生裂纹的扩展几乎全部源于爆生气体对裂纹尖端的延展作用。

图  10  $\zeta$=1.4时爆生裂纹的扩展过程

Figure  10.  Propagation of blasting induced cracks ($\zeta$=1.4)

下载: 全尺寸图片 幻灯片

3.2   不同不耦合系数下花岗岩爆破损伤模式及分析

图11为计算至裂纹停止扩展时爆生裂纹的分布情况。$\zeta$=1.0时，即耦合装药结构下，爆炸冲击压力p_t=1.32 GPa，爆生气体压力p_g=911 MPa。炮孔周围各向裂纹发育较好，在爆生气体的扩腔作用和对裂纹的驱裂作用下，径向裂纹延展距离较远，个别径向裂纹扩展至边界处。岩石介质各部分弹性应变能的累积差异造成了部分裂纹的不连续传播。$\zeta$=1.2时，p_t=3.52 GPa，p_g=305 MPa。较强的爆炸冲击应力使炮孔近区岩石得到充分破碎，径向裂纹数量极大地提升，爆生气体对裂纹的驱裂作用非常明显。$\zeta$=1.4时，p_t=1.40 GPa，p_g=189 MPa，裂纹发育模式与$\zeta$=1.0时相似，但裂纹扩展距离明显较小。$\zeta$=1.6时，p_t=0.62 GPa，p_g=130 MPa，随着爆炸冲击压力的大幅降低，爆破损伤区域减小，裂纹数量减少，爆生气体的准静态压力降低至岩石的单轴抗压强度以下，对爆生裂纹的扩展作用明显减弱。$\zeta$=1.8时，p_t=0.38 GPa，p_g=94 MPa，此时只有炮孔近区很小区域的岩石发生破坏，爆生气体压力对爆生裂纹扩展只起到微弱作用。$\zeta$=2.0时，p_t=0.17 GPa，p_g=70 MPa，爆炸冲击应力比岩石的单轴抗压强度稍大，只能造成炮孔近区花岗岩介质的起裂，而较低的爆生气体压力在没有爆炸冲击应力预裂破碎的情况下对花岗岩介质的爆破损伤过程几乎没有影响。

图  11  不同不耦合系数下花岗岩的爆破损伤情况

Figure  11.  Blasting induced damage of rock mass under different decoupling coefficients

下载: 全尺寸图片 幻灯片

由上述花岗岩爆破损伤模拟结果可知，装药不耦合系数对爆生裂纹的产生、扩展和分布具有非常重要的影响。究其原因，是不耦合装药结构对爆炸冲击压力和爆生气体压力产生较大的影响，从而改变了岩石的爆破破坏模式。随着不耦合系数的增大，爆炸冲击压力先升高后降低；爆生气体压力随着不耦合系数的增大呈快速递减趋势。由爆破破岩理论^[22]可知，粉碎区半径为药包半径的3～7倍，爆炸粉碎区域过大，会消耗掉绝大部分炸药能量，导致传入岩体内部用于切割破坏岩石的能量较少，造成炸药能量的不合理分布。因此，盲目地增加炮孔的入射能量并不能提高岩石整体的破岩效率。根据$\zeta$=1.2时岩体的爆破损伤模式可知，与$\zeta$=1.0时相比，粉碎区域仅略微增大，但爆生裂纹宽度和扩展距离得到显著增加，提升了炸药能量的有效利用率。对比模拟的其他情况，$\zeta$=1.2时爆破破岩效率最佳。因此，选取合理的不耦合系数对于提升炸药能量利用效率和岩石破碎效率具有重要意义。$\zeta$=1.4时的爆炸冲击压力比$\zeta$=1.0时只增加了80 MPa，爆生气体压力却下降了722 MPa。由爆生裂纹的最终分布情况可知，两种装药结构下，炮孔近区爆破破坏模式相似，损伤范围大致相同，说明爆生气体的准静态压力对岩石爆破破碎作用很小。然而，对于裂纹扩展，$\zeta$=1.0时较大的爆生气体压力是爆生裂纹得以进一步扩展的主要原因。由此可知，爆生气体的准静态压力对于爆生裂纹的扩展具有重要作用。随着不耦合系数的进一步增大，$\zeta$为 1.6、1.8、2.0时，爆炸冲击波和爆生气体的压力大幅降低，爆破损伤程度均比耦合装药时小，说明较大的装药不耦合系数不利于爆破破岩工作的进行。

值得注意的是，$\zeta$=1.2时，多数径向裂纹扩展至模型边界，由于设置的岩石模型区域较小，模拟中无法获取裂纹完全扩展后的分布情况，因此所获取的裂纹数量比完全扩展后的数量少。若进一步增大模型区域，会导致计算时间显著增加。考虑到本研究只关注模型区域内岩石介质的爆破损伤情况，因此研究结果仍然具有一定的合理性。另外，实际岩石具有不均质性，存在原生裂隙、节理等力学弱势面，可能对岩石的爆破损伤模式产生影响，对此本研究均不予考虑。

3.3   不同不耦合系数下花岗岩爆破损伤程度的预测模型

爆生裂纹数是反映岩石介质爆破损伤程度的重要参数。通过PFC中的DFN模块统计出不同不耦合系数下的爆生裂纹数y，如图12所示。

图  12  爆生裂纹数随不耦合系数的变化曲线

Figure  12.  Number of blasting induced crack vs. decoupling coefficient

下载: 全尺寸图片 幻灯片

由于$\zeta$=1.0和$\zeta$=1.2时的数据量较少，这里仅给出$\zeta$≥1.2时的爆生裂纹数预测模型，如图13所示，拟合公式为

图  13  指数拟合结果

Figure  13.  Exponential fitting result

下载: 全尺寸图片 幻灯片

从图13可以看出，e指数曲线对小不耦合系数拟合度更加敏感。为了提高预测模型的精度，通过模拟实验获得了$\zeta$为1.25、1.30、1.35时的爆生裂纹数，依次为14916、12135和9427，实验过程不再赘述。由此得到拟合曲线的拟合精度R²为0.9808。该预测模型对于爆破施工设计具有较高的参考价值。

4.   结　论

通过理论分析，结合数值模拟方法，研究了径向装药不耦合系数对花岗岩爆破损伤特性的影响，得到以下结论：

(1) 根据实际爆破过程，基于“颗粒膨胀加载法”，提出了一种“动”、“静”荷载混合施加的PFC爆破模拟方法，从实验中炮孔压力的变化情况来看，该方法可以很好地模拟爆炸冲击波和爆生气体共同作用下炮孔压力的变化；

(2) 分析了不同不耦合系数下爆炸冲击作用和爆生气体作用对花岗岩爆破损伤程度的影响，不耦合系数为1.2时花岗岩的爆破破岩最佳，由不耦合系数为1.4时的岩体损伤情况可知，爆生气体对爆生裂纹扩展具有重要作用；

(3) 花岗岩爆生裂纹数随着装药不耦合系数的增大近似呈指数下降，根据不耦合系数大于或等于1.2时的爆生裂纹数量，建立了岩体损伤程度与不耦合系数的预测模型，为$y = 1.13\times {{{10}^7}} {{\rm e}^{ - 5.22\zeta}}$，拟合度达0.9808，预测模型对爆破施工设计具有较高的参考价值。