切缝药包爆炸作用机理数值模拟*

申 涛，罗 宁,3，向俊庠，高祥涛

(1.中国矿业大学深部岩土力学与地下工程国家重点实验室,江苏 徐州 221116; 2.中国矿业大学力学与土木工程学院,江苏 徐州 221116; 3.中国科学技术大学近代力学系,安徽 合肥 230027; 4.贵州省公安厅治安总队,贵州 贵阳 550001)

近年来，在深部岩土工程、隧道开挖、边坡开挖、贵重石材开采等工程领域中切缝药包定向控制爆破技术应用广泛且效果显著。切缝药包的实质是在具有一定密度和强度的管状炸药包装上开有不同形状、不同角度和不同数目的切缝来同时控制爆炸应力场的分布和爆生气体对介质的准静态作用与尖劈作用，达到定向控制爆破的目的。因此,切缝药包定向控制爆破效果显著，切缝管发挥了重要作用。Fourney等[1]在炮孔中使用轴向不同方向、不同数量间隙切缝的管来控制爆炸后的定向裂缝方向和数量，并通过脆性聚酯材料、花岗岩、树脂玻璃系列模型实验证明了使用该装药方式可以获得指定方向和数量的径向裂纹，达到控制断裂的目的。王树仁等[2]、Yang等[3]开展动光弹爆炸加载实验和铅铸体爆炸实验获得了切缝药包爆炸时动态应力场、爆生裂纹产生、扩展及止裂过程的定性记录，证明了药包切缝方向炮孔壁最先产生裂隙，且承受的爆轰波作用强度大于其他方向。王汉军等[4]、罗勇等[5]、谢华刚等[6]综合运用弹塑性力学、爆炸力学与岩石断裂力学分析了切缝药包在岩石中爆炸过程中切缝管的作用和爆炸后岩石中裂纹形成、扩展条件及方向。魏晨慧等[7]、Wang等[8]、Ma等[9]对切缝药包岩石定向断裂控制爆破初始裂纹形成过程进行了数值模拟，并对影响爆破效果的不耦合系数、切缝宽度以及切缝管壁厚等主要因素进行了分析。然而，前人的研究中对于“炸药-切缝管-空气”三者共同作用下切缝药包爆轰过程[10]、爆生气体的动力学行为以及爆炸后压力时空分布问题的研究还比较欠缺。由于爆炸过程的瞬时性以及高温、高压、强破坏等极端条件限制，因此采用数值模拟方法对切缝药包爆炸机理研究是一个非常好的选择。基于高速纹影切缝药包的实验研究[10]，采用LS-DYNA软件中流固耦合算法模拟切缝药包爆炸过程，通过“炸药-切缝管-空气”三者在爆炸作用下动力学行为特征的数值模拟，揭示切缝药包爆炸作用机理。

1 计算模型

为研究“炸药-切缝管-空气”三者共同作用下切缝管爆炸过程，建立如图1所示3种模型。双缝/耦合切缝药包准二维模型，如图1(a)所示，此处双缝指沿管轴线方向开2条缝，耦合指药包直径与切缝管内径一致。在模拟过程中考虑炸药、切缝管和空气3种物质。模型尺寸参数与高速纹影实验相同，模拟中空气域大小为40 cm×40 cm。计算过程中炸药采用中心点起爆方式，空气边界设置为非反射边界。模型有限元网格和局部网格放大图如图2(a)所示，共549 00个单元。为避免计算过程中炸药和空气单元变形过大导致的求解异常，采用适用于处理大变形问题的多物质ALE算法。切缝管采用Lagrange单元算法，同时设置切缝管与炸药、空气之间基于罚函数的流固耦合约束[11]。为对比有/无切缝管时爆炸压力场的时空分布，建立无切缝管爆炸准二维模型，如图1(b)所示，其中炸药和空气单元采用多物质ALE算法。在数值模拟中，一般采用准二维模型近似模拟切缝药包爆炸过程，即在模型厚度方向取单层三维单元，同时施加厚度方向的位移约束来近似为二维问题进行模拟[8-9]。因此为探讨准二维模拟方法的可行性，以及切缝药包爆炸过程中切缝管的形态变化，建立双缝/耦合切缝药包三维模型，如图1(c)所示。考虑炸药、切缝管、堵塞和空气4种物质，模型中切缝管截面参数与准二维模拟中一致，切缝管长10 cm，切缝长8 cm，空气域大小为40 cm×40 cm×20 cm，其有限元模型如图2(b)所示，模型共3 706 000个单元。

图1 几何模型参数与边界设置Fig.1 Geometric model parameters and boundary settings

图2 有限元网格与局部放大网格Fig.2 Finite element meshes and local amplified meshes

2 数值模拟材料模型与参数

2.1 炸药JWL状态方程

炸药采用高能炸药模型，其中爆轰产物膨胀所产生的压力peos通过JWL状态方程[12]计算:

(1)

式中:V为相对体积；E为初始体积内能；A、B、R1、R2、ω为跟炸药性质相关的常数。模拟中采用的是低密度的太安(PETN)炸药，具体参数[13]为：ρ=0.88 g/cm3,A=348.62 GPa,B=11.288 GPa,R1=7,R2=2,ω=0.24,D=5 170 m/s,E0=5 025 MJ/m3,pCJ=6.2 GPa。

2.2 切缝管材料模型与状态方程

切缝管类似于薄壁壳结构，在爆炸的强冲击载荷作用下，切缝管发生大变形甚至破坏。Steinberg等[14]、Cochran等[15]通过实验研究发现，在爆炸高压下，材料应变率很大(>105s-1)，畸变律中最基本的2个材料特性参数：表征材料抗弹性畸变特性的剪切模量G和表征材料抗塑性畸变特性的屈服强度Y都依赖于压力(应力球量)和温度，相比之下，应变率效应则允许忽略。因此该模型假设应变率饱和，材料剪切模量随着压力的增大而增大，随着温度的增大而减小。这时可以把固体本构畸变律归结为确定G=G(p,T)和Y=Y(p,T)的问题，当Y0[1+β(εεi)]n≤Ymax时，该模型表达式[12]为:

(2)

(3)

固体高压状态方程通常使用Grüneisen状态方程,加载和卸载时分别为：

(4)

p=ρ0c2μ+(γ0+aμ)E

(5)

式中:c为us-up曲线的截距；S1、S2、S3为us-up曲线斜率系数；0为Grüneisen系数；a为0的一阶体积修正系数；E为体积初始内能；μ为体积参数，μ=ρρ，其中ρ为材料当前密度，ρ为材料初始密度。模拟中切缝管材料选用铜，具体参数[16]见表2。

表1 铜管参数[16]Table 1 Parameters of copper pipe[16]

2.3 理想气体状态方程

为求解空气中的爆轰压力，将空气视为理想气体，此时可采用满足γ律关系的线性多项式状态方程。在LS-DYNA中线性多项式状态方程通过下式计算压力[12]：

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E

(6)

当为理想气体时，取C0=C1=C2=C3=C6，C4=C5=-1，且γ=cp-cV，其中cp为气体定压比容，cV为气体定热比容，具体参数[13]为：ρ=1.29×10-3g/cm3,=1.4,C4=0.4,C5=0.4,E0=2.5×10-4MJ/m3。

3 结果与分析

3.1 冲击波相互作用

建立如图所示炸药-切缝管-空气计算模型来分析切缝药包爆炸作用下冲击波相互作用机理。图3(a)即双缝/耦合切缝药包截面图。定义材料初始密度ρ0和声波速度c0乘积表示材料的阻抗值，根据冲击波的相互作用和反射计算公式与阻抗匹配技术可探讨切缝管作用机理[17]:

(7)

式中:pin、pt、pr分别为入射、透射、反射压力，ρ1、ρ2和c1、c2分别为不同材料的初始密度和声波速度。

在非切缝方向，炸药爆炸后产生的冲击波与切缝管相互作用，由于炸药阻抗低于切缝管，因此当炸药中压力p1波阵面达到界面时，压力上升为p2，如图3(b)中t2时刻所示。从t3到t4时刻之间，有反向压力p1波阵面向炸药中传播，另一压力p2波阵面向管中传播。当炸药中压力波阵面和初始的冲击波卸载部分相遇时压力减小到p2-p1。压力从切缝管传向空气时，由于切缝管阻抗大于空气阻抗，因此当管中压力p2波阵面达到界面时，压力下降为p3，同时在切缝管中产生一个稀疏脉冲向切缝管中传播，当这个稀疏波传播过程中遇到初始脉冲的卸载部分，在t4时刻将形成拉伸脉冲且在2个方向传播。在切缝方向，炸药爆炸后产生的冲击波则直接作用于低阻抗的空气。综上所述，切缝药包爆炸过程中，由于切缝管阻抗大于炸药阻抗，爆炸冲击波直接作用于管，除产生透射波外，尚有向爆炸中心反射的压缩波，王树仁等[2]指出该反射波能量约为总爆轰能量的(10～13)%。

图3 冲击波相互作用机理模型Fig.3 Shock wave interaction mechanism model

3.2 切缝药包爆轰过程

高速纹影实验利用气流对光波的扰动，将肉眼不可见的气流变化转换成图像，同时结合高速摄像机将爆炸过程记录下来[18]。图4中同时给出了双缝/耦合切缝药包爆轰过程纹影实验与数值模拟结果。

图4 切缝药包爆轰过程高速纹影实验和数值模拟结果Fig.4 Experimental and numerical simulation results for the detonation process of a split-tube charge holder

图4中黄色为切缝管，灰色为炸药及爆炸后爆生气体产物。初始时刻，炸药截面为圆形，炸药爆炸后爆轰波由炸药内部传播到炸药边界处并作用于切缝管，由于切缝管的约束作用，爆生气体首先从切缝方向开始溢出。该阶段切缝管由于惯性效应并无较大变形，爆轰气体与产物沿切缝方向开始溢出且未发生绕流，其形态类似纺锤，如图4中2.25 μs时刻所示。随后爆生气体持续从切缝处往外膨胀，在爆炸冲击波和爆生气体的作用下切缝管发生膨胀变形，其形态类似哑铃，如图4中4.75 μs时刻所示。之后爆生气体继续膨胀，并不断作用于切缝管，切缝管曲率不断变小。此时切缝管对爆生气体的控制作用不断减弱，爆生气体继续膨胀并发生绕流，即爆生气体开始向非切缝方向膨胀，如图4中7.25和9.75 μs时刻所示，随后爆生气体以相似的形态不断膨胀。

炸药在爆炸过程中，冲击波的传播速度高于爆生气体的膨胀速度，图4中显示了爆炸过程节点速率等位线。爆炸初始阶段，非切缝方向冲击波瞬间从炸药内部传到切缝管内壁，而后冲击波沿切缝管传播，切缝方向在空气中形成冲击波，爆生气体也开始沿切缝溢出。随后，当冲击波传到切缝管外时，冲击波阵面与爆生气体出现界面分离，由于冲击波速度大于爆生气体的膨胀速度，两者之间的距离不断加大，界面分离越来越明显。整个爆炸过程中冲击波阵面的形态变化与爆生气体与产物形态变化类似，不过由于冲击波速度大于爆生气体膨胀速度，沿垂直切缝方向的两冲击波阵面之间相互交汇。

由于数值计算所采用炸药与实验存在差异，数值模拟所选用的炸药爆速约为实验中炸药爆速的4倍，因此两者在时间尺度上无法进行直接对比。但是从整个爆炸发展过程来看，两者结果吻合。切缝药包对爆炸能量释放及爆生气体动态力学行为实现了定向控制作用。

3.3 爆炸压力时空分布

3.3.1切缝管内部压力

图5 5个特征单元的位置Fig.5 Positions of five elements

图5所示在切缝管内炸药模型水平和垂直方向选取5个特征单元，相邻单元间距为0.15 cm。水平方向3个单元压力时程曲线如图6(a)所示，对比有/无切缝管时压力曲线和局部放大图可以看出，两者在升时阶段基本完全重合，对于中心处单元峰值略有差异，有切缝管时峰值为3.6 GPa，无切缝管时为3.4 GPa。而降时阶段两者差异非常明显，无切缝管时压力下降很快且无波动，3 μs左右均降到最低点处；有切缝管时压力波动现象非常明显，压力在下降过程中出现3个峰值过程，压力下降时间增大，7 μs左右才降为最低处。垂直方向3个单元压力时程曲线如图6(b)所示。对比有/无切缝管时压力曲线可以看出，两者在升时阶段基本完全重合，但靠近管处单元压力峰值大幅增加，有切缝管时峰值为5.8 GPa，无切缝管时为3.4 GPa。在降时阶段，有切缝管时压力下降依旧缓慢且明显波动。结合前文分析可知，由于切缝管的存在，非切缝方向爆炸后产生的冲击波遇到管后发生反射，靠近切缝管界面处压力峰值增大，使压力降时阶段出现波动，延长了冲击波作用时间；切缝方向上则对压力峰值影响较小，但同时也延长了冲击波作用时间。

图6 有无切缝管时药包内部单元压力时程曲线Fig.6 Pressure-time curves in the inside elements of the charge with and without split-tube

图7 特征单元的位置Fig.7 Positions of elements

3.3.2切缝管外部径向压力

图7所示在切缝管外空气模型水平和垂直方向选取8个特征单元，相邻单元之间间距为1.4 cm。图8(a)为水平方向单元压力时程曲线，压力峰值与增长速率较大，同时压力峰值出现的时刻也比较早。由于在空气中，爆炸产生的压力会以指数式递减，因此压力峰值随着单元距离爆炸中心距离的增大而减小。图8(b)为垂直方向单元压力时程曲线，与水平方向相比压力非常小，并且压力峰值出现的时间比较晚。通过对比不同方向上单元峰值压力与峰值压力时刻可以发现，切缝方向单元压力峰值分别是垂直切缝方向的27.2、7.5、6.5倍，且切缝方向单元峰值压力时刻相较于垂直切缝方向超前14.6、24.6、32.6 μs。由此可以看出，切缝管在爆炸过程中有效地控制爆炸能量的分布，使炸药能量重新分配，从而达到控制爆破的作用。此外，图6(c)中垂直方向随着单元距离的增大，单元压力呈现增大趋势，结合前文分析，此时压力的增大与爆生气体的绕流和切缝管运动产生的扰动有很大的关系。由于5.6 cm处距爆炸中心较远，因此在计算时间内没有采集到压力信号。

图8 有切缝管时药包外部单元压力时程曲线Fig.8 Pressure-time curves in the outside elements of the charge with split-tube

3.3.3切缝管外部环向压力

图9 特征单元的位置Fig.9 Positions of elements

如图9所示，分别选择有/无切缝管模型中距离炸药中心8 mm处切缝管外空气单元分析其压力峰值和峰值时刻。图10(a)为压力峰值极坐标图，切缝药包爆炸后，切缝方向压力峰值远大于非切缝方向。其中0°方向压力峰值为104 MPa，3.7°方向压力峰值最大为172 MPa，90°方向压力峰值最小，为11 MPa，切缝方向压力约为垂直切缝方向压力的9.5～15.6倍。无切缝管时炸药爆炸后沿环向压力峰值均为27 MPa。计算可知，切缝方向压力是无切缝管时的3.9～6.4倍，而非切缝方向压力是其0.4倍。根据前文分析可知，非切缝方向爆炸冲击波经过切缝管传播后压力会有所衰减；切缝方向爆炸冲击波和爆生气体会在该方向形成聚能作用。压力峰值出现的时刻也是一个影响爆炸效果的重要因素，如图10(b)所示，垂直切缝方向压力峰值时刻晚于切缝方向。无切缝管时爆炸压力峰值时刻沿环向一致，均为4.0 μs；有切缝管时，切缝方向压力峰值最早出现时刻为3.7 μs，沿切缝方向到非切缝方向，随着角度增大，峰值压力出现时刻不断延后，垂直切缝方向压力峰值最晚出现，为38.2 μs。综上所述，切缝管的存在，增大了切缝方向爆轰波传播速度，其最主要的作用是可以控制爆炸后不同方向的爆炸压力分布，增大切缝方向且降低非切缝方向压力峰值，同时提前切缝方向且延缓非切缝方向压力峰值出现时刻。

图10 有/无切缝管时药包环向压力峰值与峰值时刻Fig.10 Circumferential peak pressures and their corresponding times of the charge with and without split-tube

3.4 切缝管形态变化

图11 切缝管形态变化过程Fig.11 Morphological change process of split-tube

图11给出了切缝药包爆炸过程准二维模拟和三维模拟切缝管形态变化过程。图11(a)为准二维方法模拟结果，图11(b)为三维方法模拟结果俯视图。对比2图可以看出2种不同方法模拟所得结果基本一致，在爆炸冲击波和爆生气体的共同作用下，切缝管曲率不断变小，由圆形截面变为椭圆截面，随后基本变为近似矩形截面。图11(c)为三维方法模拟结果轴测视图，可以看出切药包爆炸后，切缝管的膨胀过程与起爆点位置有关，在爆炸冲击波及爆生气体的共同作用下，从起爆点处以相同的变形特征沿切缝管轴向发展。从以上定性分析过程可以看出，对于切缝药包爆炸过程的简化力学模型研究，采用准二维方法模拟是合理可行的。切缝药包爆炸过程中，切缝管与爆炸冲击波和爆生气体相互作用，可以有效控制爆炸能量的释放，在切缝方向增大爆炸冲击波、爆生气体对介质的作用，而非切缝方向产生减弱作用。

4 结 论

基于双缝/耦合切缝药包爆轰行为高速纹影实验，建立“炸药-切缝管-空气”模型，采用数值模拟方法初步研究了其定向控制爆破机理。具体分析了爆炸过程中冲击波相互作用，爆生气体的动力学行为和压力时空分布与切缝管形态变化过程。主要结论如下：(1)由于炸药、切缝管、空气三者的阻抗差异，切缝药包爆炸后爆炸冲击波与切缝管相互作用，除产生透射外，还有向爆炸中心反射的压缩波；切缝方向爆炸冲击波直接作用于空气，非切缝方向冲击波经过切缝管衰减后再作用于空气。(2)切缝药包能够在爆炸过程中有效地控制爆生气体的扩散行为，使之在爆炸过程中按照一定形态变化；爆炸冲击波与切缝管相互作用，增大了切缝管内部爆炸冲击波压力且延长了冲击波作用时间；在切缝管外部增大了切缝方向爆轰波传播速度和爆炸冲击波压力，使得切缝方向压力超前且高于非切缝方向，这在定向断裂控制爆破中有着非常积极的作用。(3)切缝药包爆炸过程中切缝管变形与起爆点位置有关，在爆炸冲击波和爆生气体的共同作用下，切缝管曲率不断变小且从起爆点处以相同的变形特征沿切缝管轴向发展；对于切缝药包爆炸过程的简化力学模型研究，采用准二维方法与三维方法可以获得一致的结果，因此在数值模拟过程中准二维方法是合理可行的。(4)采用流固耦合算法可以有效模拟爆炸过程中固体大变形问题，避免求解过程中网格的严重扭曲以及负体积造成程序求解出错问题；基于流固耦合算法的双缝/耦合切缝药包数值模拟结果与高速纹影实验结果吻合。