定向水耦合破岩聚能装药结构优化设计及数值模拟

Optimal Design and Numerical Simulation of Structure Optimization of Directional Water Coupled Rock-Breaking Polygonal Charge

摘要: 本研究针对定向破岩聚能装药，通过对固定装药半径下随着锥顶角和锥顶距变化下，有效装药率的变化趋势进行计算，开展了双向聚能装药结构参数的优化设计。并通过ALE算法对不同耦合填充的PVC聚能装药柱进行研究，分析聚能与非聚能方向的裂纹成型规律以及孔壁所受有效应力峰值，得到结论：在水耦合介质的作用下，非聚能方向的炮孔壁所受有效峰值应力更小，起到了更好的定向断裂控制效果。

Abstract: In this study, the optimized design of bi-directional polyenergy charge structure parameters was carried out for directional rock breaking polyenergy charge by calculating the trend of effective charge rate under the change of cone top angle and cone top distance with fixed charge radius. And through the ALE algorithm to study the PVC polyenergy charge column with different coupling filling, analyze the crack formation law and the effective peak stress on the hole wall in the polyenergy and non-polyenergy directions, and get the conclusion: under the action of water coupling medium, the effective peak stress on the gun hole wall in the non-polyenergy direction is smaller, which plays a better effect of directional fracture control.

1. 引言

水压定向断裂控制爆破技术是在孔壁受到冲击初始裂缝生成后，利用水具有的密度、粘度很大并且其不可压缩的特征，获得相较于空气耦合爆破时更佳的炸药能量利用率。同时水压定向断裂控制爆破的水耦合介质填充在装药爆轰后形成的水准静态压力相较空气准静态压力来说，具有更大的时间和压力水平作用于裂缝使其不断扩张[1]。

现阶段，李裕春等人[2]在线性聚能射流形成这一课题上开创性的开展了关于数值的相关模拟分析与讨论，对逆向的射流速度梯度有关的作用机理作出了分析。除此之外，例如Feng D.L.等许多学者[3]还将SPH应用到数值模拟分析在线性聚能射流中可能的侵蚀作用这一论题，最终得出了比原有传统有限元法更稳定与精准的效果。纪冲等学者[4]基于线型聚能装药切割器技术在某爆破工程中应用的实例，对线型聚能射流切割目标的原理及影响因素进行分析，提出了对于一般常用钢构件材料条件下工程爆破采用线型聚能切割器的参数设计方法，将线型聚能装药切割器技术在以工程建设为目的爆破施工中进一步地推广。学者孙飞、武双章与和发波等多人[5]-[7]共同在线性聚能装药的理论方面作出了研究，其中就以数值模拟的方法在聚能切割器上作出了正交优化的相关设置，而其最后也如愿以偿地得出了射流侵彻性能实现跨越性提高的极佳成果，除此之外，其因为对经济社会存在较好的效用，能为经济作出突出贡献。这就进一步使得线性聚能爆破技术朝着社会化方向应用和发展了。

通过分析以上研究成果可知，现阶段聚能装药和水耦合装药作为定向断裂控制爆破技术的重要组成，从理论、模拟到实践均经历了多重检验，为本文的研究奠定了基础。

本文通过求解爆轰作用下沿不同方向飞散的爆轰产物曲线，并通过蒙特卡洛打点法求得向各方向散射的装药量。利用MATLAB设计聚能装药有效部分及利用率的程序通过线性规划得到圆形装药下各参数的最优解，从而对双向线性聚能切割装置进行参数设计，在此基础上利用ALE算法对PVC材料聚能装药管下的定向断裂破岩模型进行数值模拟研究，探究不同耦合介质下聚能装药破岩发展过程和规律，探讨单孔聚能水压爆破破岩机理与岩石侵彻规律。

2. 双向线性聚能结构参数优化设计

2.1. 装药有效量及装药利用率计算

受限于光面爆破中周边孔承载破岩体积所需要的线装药密度，聚能装药不耦合系数的选择空间较小，因而本章建立以圆柱形装药截面的半径，药柱中心为原点的直角坐标系，其具体情况如下图1所示，可以看出α则是锥顶角，R是设装药的装药半径，c是其锥顶角到装药中心的距离长度。根据图1中AF = Ac求解该边界函数可求得聚能方向的装药有效部分。图给出了装药有效部分计算直角坐标图(第一象限)。

而在装药后的爆炸过程中，聚能装药会沿着图2所示的样子逐层脱离，沿着这两个脱落方向飞散的爆轰产物会在中间相交，形成一个明显的交界，求解这一交界面在截面方向上形成的曲线便可计算得到向两个不同方向散射的装药面积比。

Figure 1. Diagram of calculation parameters of effective charge

图1. 有效装药计算参数几何示意图

Figure 2. Diagram of layer by layer shedding of shaped charge

图2. 聚能装药逐层脱落示意图

李必红博士[8]通过最小二乘法将n个观测值代入求得f(x)理论的回归模型$Y={\beta }_0+{\beta }_{1}X+{\beta }_2{X}^2+e$ 的各项系数，从而得到$Y=f\left(x\right)$ 这一边界，其中自变量为X，因变量为Y，表示在二维空间中的坐标；e为随机误差，${\beta }_0$ ，${\beta }_1$ ，${\beta }_2$ 为回归模型通过观测数据拟合得到的系数，用于描述边界曲线f(x)。

由于这一方法求解f(x)的原函数较为复杂容易出错，本章选用蒙特卡洛打点法对f(x)所包围的装药面积进行计算。

蒙特卡洛计算积分的平均值法求定积分的过程是在f(x)上抽样求出一个近似面积值，$f\left(x_i\right)\ast \left(b-a\right)$ 然后计算其数学期望，用数学公式表述这一过程为：

$\begin{array}{c}S=\frac{1}{4}\left[f\left(x_1\right)\left(b-a\right)+f\left(x_2\right)\left(b-a\right)+f\left(x_3\right)\left(b-a\right)+f\left(x_4\right)\left(b-a\right)\right]\\ =\frac{1}{4}\left(b-a\right)\left(f\left(x_1\right)+f\left(x_2\right)+f\left(x_3\right)+f\left(x_4\right)\right)=\frac{1}{4}\left(b-a\right){\displaystyle \sum }_{i=1}^{4}f\left(x_i\right)\end{array}$ (1-1)

其中公式中各物理量的含义：f(x)：被积函数，即需要积分的函数；a和b：积分区间的下限和上限；x_i为在积分区间[a, b]内随机抽取的样本点；S为蒙特卡洛方法计算的积分估计值。采样n个满足均匀分布随机样本，则有：

$S=\frac{b-a}{n}{\displaystyle \sum }_{i=1}^{n}f\left(a_i\right)\simeq \theta$ (1-2)

采样点越多，估计值也就越来越接近。n为抽样的样本数量；$\theta$ 为表示期望值的符号，用于表示积分的数学期望。

这一方法假定了x在[a, b]间上均匀分布，而大多时候x在[a, b]上的分布并不均匀，因此上面方法就会存在很大的误差。

因此假设x在[a, b]上的概率密度函数为p(x)，加入到$\theta ={\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}$ 中变换为：

$\theta ={\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\theta ={\displaystyle {\int }_a^{b}f\left(x\right)\text{d}{z}_j}={\displaystyle {\int }_a^b\frac{f\left(x\right)}{p\left(a\right)}p\left(c\right)\text{d}{z}_j}\simeq \frac{1}{n}{\displaystyle \sum }_{i=1}^n\frac{f\left(a_i\right)}{p\left(x_i\right)}$ (1-3)

这就是蒙特卡洛期望法计算积分的一般形式。其中p(x)是在区间[a, b]上的概率密度函数，描述x在[a, b]上的分布情况。p(c)为概率密度函数p(x)在某个点c的取值，用于加权计算积分。

由此设计MATLAB程序。

2.2. 聚能装药爆破参数确定

根据工程实际用的装药结构，取R = 20 mm，讨论半锥角α和锥顶距c变化时固定装药半径下装药利用率η的变化

通过本节设计的计算软件进行分析得到关系如表1所示。

Table 1. Charge utilization ratio

表1. 装药利用率取值

c/cm η/% α/˚	0.75	0.8	0.9
30	45.93	45.57	44.19
40	48.74	48.08	46.68
45	54.55	49.68	48.11
50	52.35	51.34	49.8
55	50.41	53.66	49.17

从表可以看出：半锥角恒定时，装药利用率随着锥顶距c的减小整体呈现上升趋势关系。锥顶距c恒定时，半锥角与装药利用率二者呈开口向下的抛物线关系，存在着极大值，在α = 45˚时，η取得最大值54.55%。

所以，在半径为20 mm时，取α = 45˚，c = 7.5 mm，药柱的在聚能方向的利用率最大。

3. 水耦合填充双向聚能装药的数值仿真

3.1. 爆破数值模型参数设置

利用ANSYS/LS-DYNA软件对聚能装置的能量聚集作用以及在不同装药耦合介质中相同孔隙的爆炸效果进行了分析和比较，为了更精确地模拟出真实的工作环境，利用流固耦合的方法，将炸药与大气通过共节点相连，利用ALE算法将空气和爆炸物划分为液体，防止网格过度变形会影响计算的准确性。根据模型的对称性利用Lagrange网格对PVC聚能药柱和岩石建立了1/4的模型，在右侧和上方边界分别应用了透射边界条件，在左侧和下放边界分别应用了对称边条件。

在爆炸作用下，从爆炸点起，依次形成压碎区，裂隙区和弹性震动区，与装药半径相比，一般粉碎区半径其2~7倍，裂隙区是其8~150倍。夏祥等人通过对柱状装药进行了数值分析，得出裂隙区半径为装药半径的73倍。因此在装药直径为40 mm、炮孔半径R = 60 mm的前提下，为研究聚能水压装药条件下的裂缝扩展和岩体内部侵入等问题，选择了高H = 100 cm的岩体。

采用ALE算法对炸药、耦合介质材料进行了优化，对PVC管材、岩体进行了Lagrange算法的优化。选择与耦合介质材料相同的材料参数，建立ALE网格区域，组成双层网格。利用ALE多物料传输算法，实现了炸药在全网格间的流动，并利用流固耦合方法计算出了岩体的变形。该模型采用1/4对对称结构，在Z方向上设置法向约束，在对称轴线上分别设置相应的对称边界条件。由于实际岩体可看作无限大，在其外部均设置无反射边界条件。聚能装药结构如图3所示。

Figure 3. Schematic diagram of shaped charge structure

图3. 聚能装药结构示意图

装药结构断面如图所示。聚能装药结构，即将炸药放置于聚能管中，管外径Φ42 mm、内经Φ40 mm，钻孔Φ60 mm。由于炮眼直径在深度方向上要比横切面大得多，所以将其简化为一个平面应变问题，并在模型的厚度方向上设置了一个单一的网格。模型长度和宽度均为5 m，厚度为0.5 cm。本论文模拟研究所用的聚能管是一种双极型线性聚能管，其几何参数见图4。

3.2. 材料本构模型和参数选取

本次模拟中岩石材料选择RHT模型。RHT模型是基于HJC模型的而产生的一种模型，与后者相比，RHT模型不仅拥有HJC模型的压力依赖性、应变率敏感性和压缩损伤软化等特点，可以引入拉压损伤来模拟岩石等脆性材料在爆炸载荷作用下的动态响应过程。模型中引入了破坏面、弹性极限面和残余强度面，描述了失效应力、初始屈服应力和残余应力的变化规律，相关状态方程参数见表2岩石材料参数。

用关键字*MAT_HIGH_EXPLOSIVE_BURN来模拟炸药模型，同时爆炸过程用关键字*EOS_JWL状态方程来模拟。主要参数为：爆速D = 6720 m/s，爆压P_CJ = 0.097 × 10¹¹ Pa。相关状态方程参数见表3。

Figure 4. Geometric parameters of complete charge

图4. 完整装药几何参数

Table 2. Rock material parameters

表2. 岩石材料参数

密度(kg∙m^-3)	G (GPa)	A	B	C	N	Pcrush (MPa)
2600	16	0.79	1.6	0.007	0.61	16
evcrush	Plock (MPa)	evlock	Tc (MPa)	fc (MPa)	D1	D2
0.001	800	0.10	4	40	0.04	1.0
K1 (GPa)	K2 (GPa)	K3 (GPa)	rate0	EFmin	Smax
85	−171	208	1.0	0.01	7.0

Table 3. Explosive material parameters

表3. 炸药材料参数

密度(g/cm^3)	爆速(cm/μs)	PCJ (10^11 Pa)	A	B	R1	R2	ω
1.15	0.65	0.097	2.14	0.00182	4.2	0.9	0.15

3.3. 模拟结果分析

将模拟起爆后空气耦合及水耦合下的PVC管在聚能方向上的形成历程进行截图记录并分析。由模拟结果可知，5 μs时刻，水耦合下的炸药冲击波开始作用于PVC聚能管使得聚能方向上的应力出现集中；10~15 μs时刻，形成聚能射流的头部形成，但由于装药的耦合系数影响，使得此时未能完全成型的射流已经开始作用于炮孔壁，如图5未完全成型的射流作用于孔壁；25~30 μs时刻，射流在孔壁上发生侵彻，形成0.5倍孔径的初始裂纹。空气耦合下的PCV聚能管在5 μs时已经形成了仅带尖端头部的空气射流，在作用于孔壁岩体生成初始裂纹后快速发生膨胀，形成径向作用于裂纹面的静压使其延伸开展。

对比可以看到的是空气作为耦合填充介质时5~10 μs这段时间内仍然对非聚能方向的岩壁造成了较大粉碎破坏，消耗了较多能量，从而导致后续爆轰气体膨胀阶段，裂纹的生长较水耦合装药下的更慢，更短。同时对比聚能方向上的岩体初生裂纹可以看到的是，在与X轴夹角为30˚~45˚之间，空气耦合填充的工况下出现平行于聚能方向的裂纹，且随着爆轰气体的进一步作用裂纹的延伸角度逐渐增大，不利于定向断裂的控制需求。而PVC聚能管在水耦合填充的帮助下形成了更优质的初生裂纹，非聚能方向上的岩壁保存情况相较空气耦合有显著的提升。下面结合聚能以及非聚能方向上的岩壁有效应力的大小进一步分析两种耦合填充下的聚能破岩机理。

Figure 5. The imperfectly formed jet acts on the pore wall

图5. 未完全成型的射流作用于孔壁

Figure 6. Pressure cloud image under water-coupled shaped charge

图6. 水耦合聚能装药下压力云图

从图6和图7中可以看到，在两种不同介质的耦合下的，聚能装药呈现出两种风格迥异的孔壁破坏形式。由于该模拟在Z方向上设置了对称约束，所以装药的爆轰气体在以及作为耦合填充物的水，均可以有效地在孔壁初始裂纹形成后释放能量，进一步对生成地裂纹拓展延伸。因此，分析图中地爆生裂纹，可以对两种耦合介质下初始裂纹地质量以及孔内各方向的爆轰能量利用进行评估。

Figure 7. Pressure cloud image under air coupled shaped charge

图7. 空气耦合聚能装药下的压力云图

提取两种不同耦合介质下聚能方向和非聚能方向上的有效应力如图8，从数据上看，水耦合介质填充时聚能方向孔壁应力峰值均出现在15 μs，聚能方向的峰值为178 MPa，远大于非聚能方向的23.7 MPa。空气耦合填充的聚能方向的出现有效应力的峰值时间为20 μs，数值为135 MPa，非聚能方向则于45 μs时达到峰值119 MPa，两者相距离25 μs聚能方向de1峰值为非聚能方向1.13倍，可以看到单纯的PVC聚能管对岩壁有一定的保护作用，但不明显。水耦合的非聚能方向的有效应力随时间缓慢上升，在60 μs时达到峰值55.2 MPa，随后与空气耦合下的两个方向一同回落至36.5~40 MPa，如图8所示。

4. 结论

(1) 通过蒙特卡洛打点法求解沿两个方向飞散的爆轰产物的交界面曲线所包围的沿着聚能方向散射的装药量截面积，得到了固定装药半径下随着锥顶角和锥顶距变化下有效装药的变化趋势，半锥角恒定

Figure 8. Effective stress curves in different directions

图8. 不同方向上的有效应力曲线

时，装药利用率随着锥顶距c的减小整体呈现上升趋势关系。锥顶距c恒定时，半锥角与装药利用率二者总体呈开口向下的抛物线关系，存在着极大值，在α = 45˚时，η取得最大值54.55%。在半径为20 mm时，取α = 45˚，c = 7.5 mm，药柱在聚能方向的利用率最大。

(2) 水耦合下的PVC聚能装药能够相较于空气耦合填充下的装药能够形成质量更好的射流用以侵彻聚能方向的孔壁。

(3) 空气耦合填充下更容易出现在与聚能方向夹角的30˚~45˚之间的裂纹，且随着爆轰气体的进一步作用裂纹的延伸角度逐渐增大，不利于定向断裂的控制需求。而PVC聚能管在水耦合填充下形成了更集中于聚能方向上一点处的初生裂纹，非聚能方向上的岩壁保存情况相较空气耦合有显著的提升，定向断裂效果明显。

基金项目

江苏省科技厅基础研究项目(自然科学基金)面上项目(BK20211397)；江苏省地质局科技创新项目(2022KY10)。

NOTES

^*通讯作者。