不同不耦合系数下花岗岩爆破损伤特性的离散元模拟

袁增森，徐振洋，潘 博，李广尚

（1. 辽宁科技大学矿业工程学院，辽宁 鞍山 114051；2. 北京科技大学土木与资源工程学院，北京 100083）

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

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

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

1 理论分析与计算

1.1 炮孔压力计算

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

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

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

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

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

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

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

图1 PFC 中动静荷载共同施加原理Fig. 1 Sketch of dynamic and quasi-static loads application in PFC

图2 炮孔压力时程曲线Fig. 2 Curve of borehole pressure-time history

2 数值模拟

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

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

图3 PFC 岩体数值模型Fig. 3 Numerical model of rock mass in PFC

图4 平行黏结模型的基本结构Fig. 4 Structure of parallel bond model

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

图5 数值模拟结果与室内实验[20]结果对比Fig. 5 Comparison between numerical and experimental results[20]

表1 主要细观参数Table 1 Main microscopic parameters

表2 试样的主要力学参数Table 2 Main mechanical parameters of the sample

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

根据Yuan 等[9]的研究，炸药密度取 ρc= 1 000 kg/m3，爆速vc= 2 700 m/s。数值模拟中共取6 种不耦合系数，不耦合系数与炮孔半径的对应关系见表3，其中： ζ=1.0 对应耦合装药结构，其余为不耦合装药结构。由式(3)、式(4)和式(5)，可计算出不同不耦合系数下的爆炸冲击压力和爆生气体压力，分别如图6和图7 所示。

图6 爆炸冲击应力对孔壁的压力随不耦合系数的变化Fig. 6 Curve of the pressure on blasthole wall under explosion stress with different decoupling coefficients

图7 爆生气体对孔壁的压力随不耦合系数的变化Fig. 7 Curve of the pressure on blasthole wall under detonation gas stress with different decoupling coefficients

表3 不同不耦合系数对应的炮孔半径Table 3 Blasthole radius vs. decoupling coefficient

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

图8 炮孔压力时程曲线Fig. 8 Curves of blasthole pressure-time history

3 计算结果分析

3.1 爆生裂纹扩展模式

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

图9 ζ=1.4 时炮孔压力曲线Fig. 9 Curve of blasthole pressure ( ζ=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 ζ=1.4 时爆生裂纹的扩展过程Fig. 10 Propagation of blasting induced cracks ( ζ=1.4)

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

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

图11 不同不耦合系数下花岗岩的爆破损伤情况Fig. 11 Blasting induced damage of rock mass under different decoupling coefficients

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

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

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

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

图12 爆生裂纹数随不耦合系数的变化曲线Fig. 12 Number of blasting induced crack vs.decoupling coefficient

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

从图13 可以看出，e 指数曲线对小不耦合系数拟合度更加敏感。为了提高预测模型的精度，通过模拟实验获得了ζ为1.25、1.30、1.35 时的爆生裂纹数，依次为14 916、12 135 和9 427，实验过程不再赘述。由此得到拟合曲线的拟合精度R2为0.980 8。该预测模型对于爆破施工设计具有较高的参考价值。

图13 指数拟合结果Fig. 13 Exponential fitting result

4 结 论

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

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

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

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