杀爆战斗部破碎形态及破片飞散特性数值模拟研究

战斗部威力与毁伤效能指标一直是武器研制、鉴定和使用部门聚焦的核心内容。杀伤爆破型战斗部主要靠冲击波和高速破片对目标实施毁伤,战斗部壳体在爆轰压力驱动下膨胀断裂并形成破片飞散是一个复杂的物理过程,Taylor[1]在拉伸断裂假设的基础上研究了金属柱壳的破碎问题,并提出了Taylor断裂判据,但Taylor判据仅适用于加载压力较低时金属壳体的断裂行为,当加载压力较高时,断裂模式由拉伸断裂向拉剪混合或纯剪断裂过渡;Hoggatt[2]提出的绝热剪切断裂模式认为在爆炸载荷驱动下内壁面首先发生绝热剪切失稳,形成绝热剪切带向外扩展;Gurney[3]在爆轰气体密度均匀假设的基础上提出了预测金属圆管破片速度的经验公式,其认为破片初始速度主要由壳体质量与装药质量的比值确定,但是Gurney公式没有将轴向稀疏波和壳体结构考虑进去,只是对二维模型的推算。后来Hennequin[4]通过定量地分析预制破片和自然破片之间速度的区别,在Gurney公式中引入了能量修正系数F,提高计算精度;范亚夫等[5]采用量纲分析法,得到了考虑壳体强度对膨胀初速的影响的Gurney公式修正式;Hiroe等[6]基于PENT炸药线起爆技术,对低碳钢等材料柱壳的膨胀破碎进行了试验研究;汤铁刚等[7]采用高速摄影方法研究了45#钢柱壳在爆炸加载下的膨胀断裂,证实了断裂应变的应变率效应。

随着计算机技术的发展,国内外相关学者开始尝试采用数值模拟的方法来研究爆炸载荷驱动下柱状金属圆管壳体破碎及破片飞散问题[8-15]。但是对于壳体结构以及曲面母线更为复杂的战斗部,其不同部位的破碎特性、破碎形态以及破片飞散速度都鲜有研究。因此本文利用非线性动力学程序LS-DYNA,在结构尺寸和材料参数相同的情况下,采用FEM(有限元)和SPH(光滑粒子)2种方法构建了某型号杀爆战斗部的数值模型,研究了战斗部爆炸瞬间裂纹的生成与扩展规律以及不同部位的破碎形态,分析了破片沿战斗部轴线的飞散速度,并在此基础上计算了爆炸能量在各物质之间的传递及分配规律。

1 计算模型与参数

1.1 FEM模型

根据某型样弹战斗部结构尺寸,采用TrueGrid软件建立1/4三维数值计算模型,模型如图1所示。模型中炸药和空气采用ALE算法,编为一个ALE多物质组;战斗部壳体采用Lagrange算法,在两者之间使用流固耦合算法来模拟爆炸加载过程及壳体的大变形。空气网格为边长100 cm×100 cm×200 cm的长方体区域并对其施加透射边界。为保证填充炸药几何形状的精确可靠以及炸药爆轰阶段的计算精度,战斗部初始位置附近空气网格尽量致密。距离战斗部较远区域的网格较为稀疏,可避免模型体积过大导致计算时间过长这一问题。

利用非线性动力学程序LS-DYNA计算空气中爆炸的过程及现象。战斗部的网格划分情况以及装药结构如图2所示,网格尺寸为0.3 cm,网格尺寸尽量均匀且径向尺寸略大于轴向尺寸,以保证网格在壳体膨胀的过程中不会出现畸变。为避免流固耦合过程中发生初始渗漏,对耦合算法的参数进行调整,增加Lagrange壳体单元参与耦合的节点数量。

1.2 SPH模型

光滑粒子动力学方法(SPH)是一种自适应、无网格粒子形式,具有拉格朗日性质的数值模拟方法,避免了网格畸变等问题,适用于处理大变形及冲击载荷等问题。SPH模型如图3所示。首先采用TrueGrid软件将装药和战斗部壳体划分为Lagrange六面体单元网格,并根据战斗部结构的对称性建立1/4模型。然后利用前处理软件LS-PREPOST将有限元模型中的六面体单元离散成光滑粒子,最后通过修改K文件添加SPH粒子的相关属性。

对高速碰撞和非线性大变形问题进行计算时,SPH方法通常采用变光滑长度的时间积分格式[15]:

式(1)中:h(t)为随时间变化的光滑长度; v为粒子速度。当材料发生拉伸或者断裂引起粒子相互分离时,光滑长度会增加,当材料受到撞击或挤压引起粒子间距压缩时,光滑长度会减少。通过改变光滑长度来保持邻域内粒子数量的恒定。在计算时必须设定好光滑长度的最小放缩系数和最大放缩系数:初始光滑长度一定,如果放缩系数过小,可能引起数值断裂,导致弹体壳体的金属粒子在爆轰波的驱动下呈飞散状;放缩系数太大,则核函数影响域内的粒子数太多,会造成计算量的急剧增加,可能产生压缩不稳定性。在LS-DYNA中,系统默认选择所有粒子两两最小间距的最大值为初始光滑长度[16]。本模型采用的最小放缩系数为0.1,最大放缩系数为0.67。

1.3 材料模型与参数

TNT炸药采用Mat_High_Explosive_Burn材料模型和JWL状态方程,具体参数见表1[13]。其中,ρ为密度;D为爆速;PCJ为炸药爆炸产生的爆轰波C-J面压力;A、B、R1、R2、ω为试验确定的常数;E0单位体积的初始内能。壳体为高强度合金钢材料,采用Plastic-Kinematic模型,具体参数见表2[13]。其中,E为弹性模量;λ为泊松比;σs为屈服强度;C、P为应变率相关的常数;β为随动硬化系数;fs为失效体应变[14]通常取爆轰产物泄漏时刻壳体材料对应的应变值[17]。空气采用空物质材料本构模型和线性多项式状态方程。

2 计算结果和分析

2.1 爆炸过程及破碎形态分析

弹体爆炸初始阶段破裂的过程如图4所示,图中显示为Mises(v-m)应力云图。爆轰从顶端开始,60 μs爆轰结束,90～110 μs弹体从顶端自上而下出现破裂,破裂时的壳体膨胀为初始直径的1.5～1.6倍左右。由于本模型中弹体壳体采用拉格朗日算法,对于形变量过大的网格须进行删除,因此很多小尺寸破片实际上已经被删。数值计算结果显示,内层壳体的应变率在105 s-1量级,内层壳体单元在爆轰阶段迅速达到失效阈值,并被大量删除,在后续的弹体膨胀阶段,随着应变率的下降,弹体壳体开始出现“绝热剪切”破裂模式。

裂纹产生后,在钢质弹体中的扩展速度约为1 550 m/s[18]。裂纹传播的准则为裂纹出现后不能传入压缩区,裂纹的两边立即出现卸载波,在卸载区内不会出现新的裂纹。但是目前在FEM方法中尚没有计算模型可以包含绝热剪切、裂纹生成及扩展等诸多效应。

SPH算法每执行一个运算步,都会对所有粒子的邻域进行重新搜索,重新确认。并计算邻域内所有粒子的速度、应力应变状态。如果邻域内某个粒子的塑性应变满足Plastic-Kinematic模型中设定的失效准则,该粒子会被认定脱离接触,从而开始裂纹的演变。

图5给出了SPH模型壳体破裂过程的塑性应变云图以及局部壳体的破碎形态。T=76 μs时,炸药爆轰基本结束,弹体顶端产生较大塑性变形,局部发生破裂,最大塑性形变量为195%,弹体圆弧段和圆柱段发生径向膨胀,弹体表面塑性应变区域纵横交错。T=120 μs时,塑性应变区域进一步发展形成裂纹,由图中可以看出,纵向裂纹要多于横向裂纹,破片的初始形态多为竖条状,圆弧段和圆柱段结合部位以及圆柱段中部有明显横向裂纹出现。为了验证该型杀爆战斗部裂纹扩展形态,采用土中爆炸的方式,限制破片生成和飞散,获得裂纹扩展阶段的圆柱部大块破片如图5(c)所示,可以看出弹体裂纹主要以绝热剪切的形式沿纵横扩展,绝热剪切面与径向呈40°～60°夹角,横向裂纹较少,主要呈现拉伸断裂模式,与SPH模型仿真结果基本吻合。T=398 μs时刻,爆炸过程基本结束,可以看出弹体圆弧段和圆柱段末端破碎性较好,这是由于弧度和折角的存在,爆轰波对壳体材料微元的抛掷方向发生较大变化,从而导致这些部位更容易产生横向裂纹。

由图5(e)所示的破片形态可以看出,靠近弹体顶端的破片比较短小,分布方向有轴向的、横向的也有斜向的。而圆弧段下部和圆柱段中下部破碎性较差,且多为纵向分布的长条状破片,在水介质破碎性试验中收集的质量在50 g以上的破片多集中在该区域。

2.2 破片速度分析

为了定量分析破片速度的分布规律,从弹体顶端到底部中心沿轴线方向等距离选取观测点。图6—图9给出了2种方法在弹体壳体不同位置计算得到的速度对比曲线。由图6可知:SPH模型弹体顶端破片速度在80 μs的时间内加速到310 m/s的平台,而FEM模型加速到350 m/s的平台却用了240 μs,并且可以明显的区分出爆轰波加载与爆轰产物加载2个过程。SPH方法将爆炸加载的流体力学过程简化为粒子应力应变状态的传递[19],这种应力传递比流固耦合算法更加简洁高效,将爆轰波加载与爆轰产物加载过程合并。这也导致SPH方法计算的破片速度普遍偏低,弹体顶端破片最终速度较FEM方法偏低了7.3%,圆弧段偏低了17.4%,圆柱段偏低了3.2%,弹底中心偏低了19.2%。

圆弧段和弹底中心的计算偏差最大,主要原因是这2个位置的弹体结构变化较大,轴向飞散破片会受到径向飞散破片的速度牵连,导致破片速度降低,并发生方向偏转。从图9可以看出SPH模型弹底中心破片速度从峰值1 290 m/s下降到945 m/s,下降幅度达到26.7%。

图7、图8、图9所示FEM模型没有出现速度牵连导致速度明显下降的现象,主要是因为受牵连位置,很容易达到应力应变失效阈值,从而导致单元删除,失去牵连作用,单元被删除的位置更加容易产生应力集中,从而导致材料的断裂。

为了获得破片速度沿弹体轴向的分布规律,在SPH模型中,从弹体顶端开始沿轴线方向等距离选取选择11个观测点,其中最后一个观测点为弹体底部中心位置。FEM模型由于部分位置数据缺失,选择了非等间距的9个有效观测点。其中弹体顶端的位置是L=0 mm,速度均选取300 μs之后的稳定速度。由此得到各模型破片速度沿壳体轴向分布曲线如图10所示。由于弹体末端和弹体底部结构非线性过度,速度也出现了跳变,除此之外,对其他观测点进行多项式拟合,得到了弹体破片速度沿壳体轴向的分布函数:

其中: V表示破片速度,m/s;L表示观测点的轴向位置,cm,该分布规律可用于杀伤效应评估。

为了验证仿真计算结果的准确性,采用基于迎风面积的单靶距测速法,对某型样弹战斗部进行速度测量。根据破片速度矢量的分布情况,靶板所测速度为弹体壳体圆柱部分的速度,所测破片飞散角基本都在90°左右,对应数值模型中25 cm<L<40 cm的壳体区域。测试结果如表3所示,其中K值为综合了破片形状因子和阻力因子的符合系数,是通过20块典型破片的质量和迎风面积计算出来的平均数据,再根据固定靶距处破片的平均速度,推算得到破片初始速度平均初速为1 143.95 m/s。

FEM方法计算的圆柱段最大速度为1 094 m/s(除去弹底),相对试验数据偏低了4.4%,SPH方法计算圆柱段的最大速度为1 058 m/s,相对试验数据偏低了7.5%。

杀爆弹战斗部壳体在爆轰波的驱动下,前期会出现塑性拉伸,随后会从壳体内部产生绝热剪切裂纹,通常呈斜45°向外扩展。爆炸驱动下壳体材料失效断裂是一种复杂的物理力学机制现象,目前还没有成熟的材料模型和状态方程可以准确描述该过程,在塑性拉伸阶段,Plastic-Kinematic模型和Johnson-Cook模型都可以描述,但是绝热剪切带的产生与扩展就很难进行计算了。由于无法计算绝热剪切现象,仿真计算得到的塑性变形能会偏高,从而导致壳体动能偏低,因此计算结果偏低4.4%和7.5%是可以接受的。

2.3 能量传递分析

通过对杀爆战斗部的能量分配和传递过程进行定量分析,可以快速评估破碎形态和飞散初速,炸药释放的能量一部分消耗于壳体的变形、破碎和破片的飞散,另一部分消耗于爆轰产物的膨胀和形成空气冲击波。壳体的变形与破碎消耗的能量约占炸药总能量的3%,对于破碎性较好的杀爆战斗部,这个比重可能达到6%[19]。简单来说,破片内能即塑性变形能越高,其破碎性越好。

图11所示为不同模型炸药及其爆轰产物内能和破片动能的变化曲线。FEM模型的内能曲线实际上是炸药和空气所构成的欧拉场的内能曲线,其中包含了一部分空气内能,单位体积空气的初始内能为2.5×10-7 MJ/ cm3[20],但所占比重较少,而SPH模型不包含空气域,因此FEM模型的内能曲线在初始时刻略高于SPH模型。2条内能曲线的拐点出现在T=70 μs时刻,该时刻为爆轰结束时刻,在这之前炸药内能迅速释放,之后的内能残余主要保存在爆轰产物里面,2个模型的破片动能曲线也是在这个时刻之后达到峰值。由于FEM模型壳体单元不断删除,导致后续破片动能曲线的逐渐下降。由图中显示的能量变化过程可以看出,破片动能峰值为5.6 MJ,占爆炸总能量的36%;用于生成空气冲击波的能量主要是爆轰产物的残余内能,占爆炸总能量的53%;另外11%的能量主要分配给爆轰产物的动能以及壳体拉伸、破裂所产生的塑性变形能。结果符合杀伤爆破型战斗部的经验数据[21]。

图12 给出了不同模型破片内能即塑性变形能的变化曲线,FEM模型计算的塑性变形能要明显高于SPH模型,其峰值到达1.5 MJ,占到了爆炸总能量的9.6%,随后由于失效单元大量删除,塑性变形能的计算结果显著下降。SPH模型的塑性变形能最终稳定在0.9 MJ,占爆炸总能量的5.7%。

3 结论

本文采用FEM和SPH方法构建杀伤爆破型战斗部的数值模型,对其破碎形态和破片飞散特性进行了计算分析,研究结果表明:

1) SPH方法和FEM方法相比,可以计算弹体不同部位的破碎形态,较好的模拟裂纹扩展和破片生成。

2) SPH方法忽略了爆轰产物加载过程,而FEM方法无法计算破片速度牵连作用。FEM和SPH数值计算结果与试验的偏差分别为4.4%和7.5%,基本满足精度要求。

3) SPH方法对爆炸能量传递过程的计算较FEM更加准确。破片动能占爆炸总能量的36%;爆轰产物的残余内能用于生成空气冲击波,占爆炸总能量的53%;另外11%的能量主要分配给爆轰产物的动能以及壳体拉伸、破裂所产生的塑性变形能,其中塑性变形能占爆炸总能量的5.7%。

[1] TAYLOR G I.The fragmentation of tubular bombs[J].Advisory Council on Scientific Research and Technical Development,1963,5(1):202-320.

[2] HOGGATT C R,RECHT R F.Fracture behavior of tubular bombs[J].Journal of Applied Physics,1968,39(3):1856-1862.

[3] GURNEY R W.The initial velocity of fragments form bombs,shells,and grenades[J].BBL Report,1943,405.

[4] HENNEQUIN E.Influence of the edge effects on the initial velocity of fragments from a warhead[J].NASA STI/Recon Technical Report N,1986,86:32593.

[5] 范亚夫,段祝平.关于Gurney公式的强度效应修正[J].弹道学报,2009,21(1):27-30.FAN Yafu,DUAN Zhuping.Strength effect correction of Gurney formula[J].Journal of Ballistics,2009,21(1):27-30.

[6] HIROE H,FUJIWARA K,HATA H,et al.Deformation and fragmentation behaviour of exploded metal cylinders and the effects of wall materials,configuration,explosive energy and initiated locations[J].International Journal of Impact Engineering,2008,35:1578-1586.

[7] 汤铁刚,李庆忠,孙学林,等.45钢柱壳膨胀断裂的应变率效应[J].爆炸与冲击,2006,26(2):129-133.TANG Tiegang,LI Qingzhong,SUN Xuelin,et al.Strain-rate effects of expanding fracture of 45 steel cylinder shells driven by detonation[J].Explosion and Shock Wave,2006,26(2):129-133.

[8] 孔祥韶,吴卫国,杜志鹏,等.圆柱形战斗部爆炸破片特性研究[J].工程力学,2014,31(1):243-249.KONG Xiangshao,WU Weiguo,DU Zhipeng,et al.Research on fragments characteristic of cylindrical warhead[M].Engineering Mechanics,2014,31(1):243-249.

[9] 蒋建伟,张谋,门建兵.小口径榴弹自然破片形成过程的数值模拟[J].弹箭与制导学报,2009,29(1):114-117.JIANG Jianwei,ZHANG Mou,MEN Jianbing.Numerical simulation of the formation of natural fragments from a small caliber shell[J].Journal of Projectiles,Rockets,Missiles and Guidance,2009,29(1):114-117.

[10] GLANVILLE J P,FAIRLIE G,HAYHURST C,et al.Numerical simulation of fragmentation using AUTODYN-2D and 3D in explosive ordnance safety assessment[C]//The 6th International Explosive Ordnance Symposium.Canberra.Australia.2003:29-31.

[11] TANAPORNRAWEEKIT G,KULSIRIKASEM W.Effects of material properties of warhead casing on natural fragmentation performance of high explosive (HE) warhead[J].World Academy of Science Engineering and Technology,2011,59:1275-1280.

[12] GUILKEY J E,HARMAN T B,BANERJEE B.An Eulerian-Lagrangian approach for simulating explosions of energetic devices[J].Computers &Structures,2007,85(11/14):660-674.

[13] HUANG Guangyan,LI Wei,FENG Shunshan.Axial distribution of fragment velocities from cylindrical casing under explosive loading[J].International Journal of Impact Engineering,2015,76(8):20-27.

[14] LIANG Minzu,LI Xiangyu,LU Fangyun.Effects of U-notches on the dynamic fracture and fragmentation of explosively driven cylinders[J].Theoretical and Applied Fracture Mechanics,2015,77(3):50-58.

[15] 王新征,张松林,邹广平.内部短药柱爆炸作用下钢筒破裂的数值分析[J].高压物理学报,2010,24(1):61-66.WANG Xinzheng,ZHANG Songlin,ZOU Guangping.Numerical analysis on fragmentation properties of the steel cylinder subjected to detonation of internal short cylinderical explosive charge[J].Chinese Journal of High Pressure Physics,2010,24(1):61-66.

[16] HALLQUIST J O.LS-DYNA keyword user’s manual[Z].Livermore:Livermore Software Technology Corporation,2001.

[17] 张志春,强洪夫,孙新利.杀伤战斗部破片飞散特性的数值模拟[J].解放军理工大学学报(自然科学版),2008,9(6):671- 675.ZHANG Zhichun,QIANG Hongfu,SUN Xinli.Numerical simulation for fragments flight characteristics of high explosive projectile[J].Journal of PLA University of Science and Technology,2008,9(6):671-675.

[18] FARAUD M.SPH simulations of debris impacts using two different computer codes[J].International Journal of Impact Engineering,1999,23(1):249-260.

[19] 许猛,杜忠华,刘宇珩,等.紫铜爆炸成型弹丸等效弹丸侵彻规律研究[J].弹道学报,2024,36(1):70-77.XU Meng,DU Zhonghua,LIU Yuhang,et al.Study on penetration law of equivalent projectile of copper exploively-formed projectile[J].Journal of Ballistics,2024,36(1):70-77.

[20] 梁晶,熊振宇,裴东兴,等.爆炸冲击波信号的降噪处理[J].火力与指挥控制,2022,47(8):79-84.LIANG Jing,XIONG Zhenyu,PEI Dongxing,et al.research on noise reduction process of explosion shock wave signal[J].Fire Control &Command Control,2022,47(8):79-84.

[21] 韩夫亮.预制破片战斗部爆炸威力场数值模拟分析[D].长春:吉林大学,2008.HAN Fuliang.The simulation analysis of prefabricate fragment warhead explosion power field[D].Changchun:Jilin University,2008.