一种高精度欧拉算法在成形装药金属射流问题中的应用＊

陈千一

（中国商用飞机有限责任公司上海飞机设计研究院，上海201210）

空心炸药能将炸药爆炸产生的能量汇聚到一起，在局部产生极强的作用力；如果给空心炸药的镂空部分加上衬垫，炸药爆炸产生的高压气体可以将衬垫压溃，形成一股高速射流，这就是成形装药。

成形装药金属射流属于爆炸加载问题，涉及到炸药的爆轰过程、不同相态物质之间的相互作用以及固体的大变形甚至破碎等高度非线性问题，目前主要靠实验手段研究此类问题，数值模拟暂时只能起到辅助作用。但是实验成本昂贵，又具有一定的危险性，而数值模拟成本低廉、高效安全，且可以获取一些实验技术无法获取的信息。近年来，一些研究人员尝试用现代计算流体力学的高精度格式来模拟爆炸问题。刘凯欣课题组将一种改进的CE／SE格式与粒子水平集方法相结合，成功模拟了爆炸焊接［1］和金属射流［2］等爆炸加载问题，计算结果与实验及他人的计算结果吻合得较好，证实了计算方法的可靠性。

材料在高速冲击载荷作用下的温度计算一直是一个棘手的问题，温度对材料的动力学行为有着重要的影响，从塑性流动到结构性相变、甚至液化或者汽化，其中都有温度的作用。大多数文献只考虑由塑性变形引起的温度变化，而忽略了冲击温升［3－6］。对于冲击温升的计算已有一些教科书给出了详细的讨论［7－8］，传统的冲击温升计算方法在欧拉型数值程序中会遇到2方面的问题：首先是程序中难以判断加载过程与卸载过程，且在冲击动力学问题中会出现反复的加载与卸载的循环过程，使得实际操作起来更加困难；其次，由于欧拉程序记录的是空间点的物理量，而不是跟踪某一物质点的始末，因此无法通过状态量直接求出温度的变化。

本文中，构建一套高精度的求解成形装药金属射流问题的数值方案，提出一种适用于欧拉程序的同时考虑塑性温升和冲击温升的计算方法。首先使用固体炸药爆炸的经典算例对自行开发的数值程序进行验证，然后对成形装药金属射流问题进行模拟，给出炸药的起爆过程、金属衬垫的失效以及射流的形成过程。讨论不同衬垫形状对射流的形状、速度和温度的影响。

1 算法描述

1.1 控制方程

使用流体弹塑性模型求解金属射流，一般将应力张量拆开表示，

式中：σij为应力张量的分量，sij为偏应力张量的分量，p为静水压。当固体物质受到强冲击载荷时，其力学行为接近于流体，此时可以忽略应力偏量。

将轴对称问题简化到柱坐标（r，z，θ）的r－z平面上，不考虑外力、热传导和热扩散，方程组的守恒形式可以写成

式中：Q为守恒量矢量，E和F分别为r和z方向上的流通量矢量，S为源项矢量，它们的表达式为

式中：ρ为物质的密度，u和v分别为r和z方向上的速度分量，E＝ρe＋（1／2）ρ（u2＋v2）为物质的体积能量，e为物质的质量内能，σrr、σzz、σθθ和σrz为应力张量的分量，srr、szz、sθθ和srz为偏应力张量的分量，εp为等效塑性应变，T为物质的温度。对于固体炸药的气体产物，只需要求解前4个方程。衬垫材料失效后，开始塑性流动并形成射流。随着温度的上升，应力偏量趋近于0，金属材料表现得更像是无黏流体。使用改进的时－空守恒元解元（CE／SE）格式［9－10］求解上述控制方程组。

1.2 温度的计算

考虑塑性温升和冲击温升对射流温度的影响。塑性温升的计算公式为

由于金属衬垫为固体材料，而固体材料的冲击绝热线与卸载等熵线十分接近［7］，假设加载过程和卸载过程都是等熵的，因此有

式中：dS为材料的熵增，cV为材料的比定容热容。由γ／V＝γ0／V0以及dV＝－1／ρ（d）2ρ，可以得到

这种近似的方法比较适合于欧拉型的数值程序，计算时不需要判断加载还是卸载。将塑性温升和冲击温升结合并写成欧拉形式

方程（3）中的ST就可以写为

2 算法验证

一端固壁一端自由的一维炸药棒是一个模拟固体炸药的经典算例，这个问题的爆轰过程的数值模拟已经十分成熟，本文中将它作为一个基准算例来验证我们的程序。一根长0.1m的TNT炸药棒从固壁端起爆，用燃烧函数［11］和JWL状态方程［12］描述炸药的爆轰过程和爆轰产物。TNT炸药的材料参数如表1所示，其中ρ0为炸药的密度，pCJ为炸药的C－J爆压，vd为炸药的爆速，e为炸药的内能，A、B、R1、R2和ω为JWL方程中的系数。

表1 TNT炸药和B炸药的材料参数Table 1 Material parameters of TNT and composition B

图1为起爆10μs后的压力分布曲线，分别使用了5种不同的网格尺寸：在0.1m的长度上划分200、500、1 000、2 000、4 000个网格。从图中可以看到，随着网格数n的增加，数值结果逐渐收敛。当网格数大于或等于1 000时，计算结果不受网格尺寸的影响。图2为炸药起爆后14μs内不同时刻的压力分布曲线，时间间隔为2μs。图中的虚线是由实验测得的TNT炸药的C－J爆压（21GPa），可以看到，随着爆轰过程的推进，压力峰值越来越接近C－J爆压。

图1 不同网格尺寸下起爆后10μs时的压力分布Fig.1 Pressure distribution at 10μs after detonation with different resolutions

图2 起爆后14μs内不同时刻的压力分布图，时间间隔2μsFig.2 Pressure profiles along the TNT slab at an interval of 2μs until 14μs after detonation

3 数值结果与分析

对4种不同形状的金属衬垫进行模拟，考察金属衬垫的形状对射流的影响。图3为初始时刻4种不同形状金属衬垫的构型图。图3（a）～（c）是圆锥形金属衬垫，它们的顶角θ依次为43.6°、77.3°、120°；图3（d）的衬垫是球面的一部分，球面半径为0.06m。4种构型的金属衬垫的装药半径均为0.04m。装药采用B炸药，其材料参数如表1所示。金属衬垫全部为OFHC铜，使用Johnson－Cook本构模型和 Mie－Gruneisen状态方程模拟这种材料。炸药的起爆方式可分为点起爆和面起爆，本文中全部采用面起爆。

图3 初始时刻4种不同形状衬垫的构型图Fig.3 Configurations of the shaped charges with different－shaped liners

图4为不同形状的金属衬垫在装药起爆后50μs时的变形图和密度分布等值云图，图4（a）～（c）是圆锥形衬垫的计算结果，图4（d）是球面衬垫的计算结果。考虑到球面衬垫的顶角为180°，可以看到，随着金属衬垫顶角的增大，射流的长度逐渐变短。这是因为顶角为43.6°的衬垫产生的射流速度最大，在同样的时间里射流的变形最大，被拉伸的长度最长。此外，还可以看到，随着衬垫顶角的增大，射流回流部分（即翅膀以上部分）所占的比重逐渐减小，也就是说，随着衬垫顶角的增大，金属衬垫的利用率有所增加。从密度云图可以看到，射流后部的表面，也就是与炸药的接触面附近密度最大，是由于受到炸药的爆炸载荷作用后，材料急剧压缩；而前部的射流密度较低，这是由于射流被拉长导致局部的体积膨胀。

图4 t＝50μs时不同形状的金属衬垫形成的射流的变形图与密度分布Fig.4 Deformation and density contour plots of the jets at 50μs from liners with different shapes

在程序中监测射流物质的最大速度和最大温度并每隔1μs记录1个值。图5（a）、（b）分别为4种不同形状的金属衬垫形成的射流的速度最大值和温度最大值的时间历程曲线。可以看到，随着衬垫顶角的增大，射流速度逐渐减小。温度的时间历程曲线略有不同，在受到爆炸载荷作用的初期，衬垫急剧压缩，塑性变形产生塑性温升，密度增大产生冲击温升，材料内部的温度迅速达到最大值。随着射流的形成，尽管塑性变形仍然使得温度上升，但是由于射流的伸长，材料密度快速下降，冲击温升部分开始贡献负温升。因此，在爆炸载荷作用的初期，材料的温度随着衬垫顶角的减小而增大，而当射流形成并趋于稳定之后，4种衬垫顶角的射流温度十分接近。

图5 不同形状衬垫下射流内部速度最大值和温度最大值的历史曲线Fig.5 Maximum－velocity histories and maximum－temperature histories of the jet from liners with different shapes

4 结 论

模拟了由4种不同形状衬垫构成的成形装药，通过分析发现，在爆炸载荷作用的初期，材料的温度随着衬垫顶角的减小而增大，而当射流形成并趋于稳定之后，4种构型的衬垫产生的射流温度十分接近。随着衬垫顶角的增大，射流的长度逐渐变短，速度逐渐减小，对金属衬垫的利用率有所增加。

本文中目前只完成了射流的形成过程，射流的穿甲问题是接下去需要开展的工作，这涉及到另一个内容十分丰富的领域——穿甲力学。此外，还可以考虑将现有的计算平台向三维扩展。

［1］Liu Wei－dong，Liu Kai－xin，Chen Qian－yi，et al.Metallic glass coating on metals plate by adjusted explosive welding technique［J］.Applied Surface Science，2009，255：9343－9347.

［2］Chen Qian－yi，Liu Kai－xin.A high－resolution Eulerian method for numerical simulation of shaped charge jet including solid－fluid coexistence and interaction［J］.Computers ＆ Fluids，2012，56：92－101.

［3］Udaykumar H S，Tran L，Belk D M，et al.An Eulerian method for computation of multimaterial impact with ENO shock－capturing and sharp interfaces［J］.Journal of Computational Physics，2003，186（1）：136－177.

［4］Tran L B，Udaykumar H S.A particle－level set－based sharp interface Cartesian grid method for impact，penetration，and void collapse［J］.Journal of Computational Physics，2004，193（2）：469－510.

［5］Ma S，Zhang X，Lian Y，et al.Simulation of high explosive explosion using adaptive material point method［J］.Computer Modeling in Engineering ＆ Sciences，2009，39（2）：101－123.

［6］Teng X，Wierzbicki T.Evaluation of six fracture models in high velocity perforation［J］.Engineering Fracture Mechanics，2006，73（12）：1653－1678.

［7］Meyers M A.Dynamic behavior of materials［M］.New York：Wiley，1994.

［8］汤文辉，张若棋.物态方程理论及计算概论［M］.长沙：国防科技大学出版社，1999.

［9］Wang Jing－tao，Liu Kai－xin，Zhang De－liang.An improved CE／SE scheme for multi－material elastic－plastic flows and its applications［J］.Computers ＆ Fluids，2009，38（3）：544－551.

［10］Chen Qian－yi，Wang Jing－tao，Liu Kai－xin.Improved CE／SE scheme with particle level set method for numerical simulation of spall fracture due to high－velocity impact［J］.Journal of Computational Physics，2010，229：7503－7519.

［11］李德元，徐国荣，水鸿寿，等.二维非定常流体力学数值方法［M］.北京：科学出版社，1998.

［12］Lee W H.Computer simulation of shaped charge problems［M］.World Scientific：Singapore，2006.