基于ICEM CFD和AUTODYN的榴弹破片场仿真研究

刘立，鲁冬林，曾拥华，王小龙，柏林元

基于ICEM CFD和AUTODYN的榴弹破片场仿真研究

刘立1，2，鲁冬林1，曾拥华1，王小龙1，柏林元1

（1.解放军理工大学，野战工程学院，江苏南京210007；2.69241部队，新疆昌吉831707）

常规研究榴弹破片的途径是采用空中静爆或水井爆炸装置，收集并分析破片的数质量情况。运用ICEM CFD前处理程序及AUTODYN软件对某型榴弹的爆炸过程进行了仿真模拟，通过Gurney公式及Mott公式对榴弹自然破片的速度、质量分布情况进行分析验证。结果表明，仿真模拟同实际爆炸过程一致，破片场速度、质量与工程计算结果吻合较好。

ICEM CFD；AUTODYN；破片场

榴弹是弹丸内装有炸药，利用爆炸时产生的壳体破片和炸药爆炸的能量以形成杀伤和爆破作用的弹药的总称，在军事实践中应用非常广泛[1]。榴弹作为一种较为常见的杀伤手段，是工程机械等防护性能较弱的装备在战时面临的主要威胁。通过工程计算及实装试验可以得知，大部分防护能力较弱的装备在被榴弹直接命中时将造成重度损伤乃至完全报废。本文重点研究榴弹爆炸形成的破片场的分布规律及特性。

ICEM CFD是一款功能强大的前处理软件，不仅可以为主流的CFD软件提供高质量的网格，而且还可以完成多种CAE软件的前处理工作，能够利用ICEM CFD方便快捷地创建有限元网格模型[2－3]。AUTODYN是一款显式有限元分析程序，它提供很多高级功能，用来解决固体、流体、气体及其相互作用的高度非线性动力学问题，尤其是在解决军工问题方面有其独特优势[4]。

针对榴弹爆炸破片的形成过程，当前主要是采用战斗部空中静态爆炸试验或采用水井爆炸装置的方法，试验成本高危险大，不便进行定量分析。利用仿真软件进行数值模拟的方法不但能够有效节约成本，而且可以更加直观地反映爆炸过程，更加快速地获得破片的质量分布和速度变化等参数，是一种重要的装备易损性分析技术手段。

本文利用ICEM CFD对榴弹壳体模型进行网格化，而后导入AUTODYN中进行爆炸过程仿真模拟，对破片场参数及工程计算结果进行对比分析。对于杀爆榴弹，其主要杀伤手段为爆炸产生的冲击波及破片群，破片的杀伤威力主要由生成的破片数量、破片质量分布、破片速度等因素决定。本文研究破片场分布规律，不考虑破片在空气中速度衰减、杀伤面积及对目标的侵彻能力，仅对破片速度及质量分布情况进行讨论。

1 榴弹爆炸破片场工程计算

1.1 破片速度

榴弹爆炸后，破片初速可应用Gurney公示进行计算[5]：

式中，V0为破片平均初速（m/s）；为Gurney常数（m/s）；C为炸药质量（kg）；m为壳体质量（kg）；其中C/m称为弹药的爆炸载荷系数。

根据炸药理论，E与爆速D存在以下关系：

式中K为爆轰产物的绝热系数，通常取K=3，则：

该型榴弹有效装药质量C=3.27 kg，壳体质量m=14.8 kg，战斗部长度L=0.452 m，TNT炸药爆速D=6 930 m/s，由式（3）得其Gurney常数为2 440.58 m/s.由式（1）可算得破片平均初速为1 086 m/s.

1.2 破片质量

大量研究资料表明，对于非预制破片榴弹，Mott模型在预测每一质量等级破片的分布方面与试验数据的吻合程度令人满意[6]，其表达式为：

式中，N（M＞m）表示质量大于m的破片数量；N0为破片总数量；M0为所有破片的总质量；m¯为破片平均质量；μ为特征质量参数，标志弹体的破碎特性，取决于弹体的外壳结构、材料以及内部装药的性能。经查有关资料，该榴弹的μ值取1.766×10-3kg.

由式（4）计算，该型榴弹共生成4 270个大小不一的破片。破片质量分布如表所示。

表1 工程计算破片质量分布

经分析，随着破片质量的增大，其数量呈逐渐下降趋势，87%的破片质量小于8 g.

2 榴弹爆炸破片场有限元仿真

2.1 战斗部壳体网格划分

保险、引信、电线等对仿真结果影响微小的部分予以忽略，仅考虑榴弹壳体。该型榴弹战斗部具有对称性，为降低计算量、便于查看所建模型网格情况及战斗部爆炸时装药的膨胀、溢出等过程，不必建立整个战斗部模型，采用四分之一模型。

六面体网格在计算精度、变形特性、划分网格数量、抗畸变程度及再划分次数等方面具有明显优势[7]，是三维模型网格划分的首选。因此选用六面体网格对榴弹战斗部壳体进行有限元网格划分。

通常按照调整几何拓扑结构、进行分块、建立映射关系、设置尺寸生成网格、检查网格质量的步骤进行ICEM CFD网格划分。为防止不够规则的四分之一壳体自动生成的网格出现乱序，导致网格质量下降甚至“孔洞”问题的产生，壳体要进行合理分块。对于两个壳体前后两端等扭曲度较大的位置，甚至要进行多次分块。同时，合理增大节点数量可以有效地改善网格质量[8]。共生成65 526个六面体网格，如图1所示。

图1 四分之一榴弹壳体模型网格

使用ICEM CFD提供Pre-mesh Quality功能对划分的网格进行质量考察，如图2、图3所示，最小雅克比矩阵与最大雅克比矩阵行列式的比值Determinant 2×2×2为0.765～0.999，介于0.5～1之间；每个网格单元的最小内角Angle的值为37.08～89.91，介于30°～90°之间，表明网格质量符合要求。

图2 四分之一榴弹壳体Determinant 2×2×2值

图3 四分之一榴弹壳体Angle值

2.2 算法选择及材料定义

拉格朗日算法可以很好地描述固体材料的行为，且相对于其他算法计算速度更快，但在固体有较大变形的时候会导致网格扭曲，通过采用附加侵蚀的方法可以克服避免网格发生畸变，适用于榴弹壳体及TNT装药。欧拉算法适合于描述液体和气体的行为，自由边界面和材料的交界面可以通过固定的欧拉网格来表达，由于网格是固定的，大变形或者有流动的情形并不会导致网格畸变，适用于本模型中的空气。

榴弹壳体选择30CrMnSi，装药选择TNT，破片在空气场中飞行。三种材料最后通过定义拉格朗日-欧拉耦合作用进行计算。拉格朗日单元同欧拉单元耦合类型选择Fully Coupled.由于建立的空气场范围有限且空气阻力对破片初速影响较小可忽略不计，因此边界条件设置为无限边界，Boundaries选择Flow_out.

空气及TNT参数取自ATUODYN材料库，空气采用Ideal Gas状态方程，TNT采用JWL状态方程。JWL状态方程可以用来计算由化学能转化成的压力，其压力和能量的关系可表示为：

式中A1、B1、R1、R2、ω为材料常数；p为压力；V为相对体积；E0为初始比内能。

TNT的材料参数如表2所示。

表2 TNT材料参数

壳体选用Shock状态方程、von Mises屈服准则和Principal-Strain失效模型[9]。其中von-Mises屈服准则是von Mises于1913年提出的一个屈服准则，即最大弹性畸变能准则，即当弹性畸变能达到某一定值时，材料发生屈服。在描述榴弹爆炸时壳体的膨胀破裂过程时应用von Mises屈服准则是十分有效的[10]。Von Mises应力在三维空间中被定义为：

式中，sv为von-Mises应力，sx、sy、sz为正应力，sxy，syz，sxz为剪应力[11]。

30CrMnSi的主要参数如表3所示。

表3 壳体材料参数

2.3 模型建立及运算

将利用ICEM CFD建立好的榴弹壳体网格模型导入AUTODYN中，空气场模型与榴弹中心线同轴。由于榴弹壳体具有不规则轮廓，在其包围的内侧填充TNT时，要采用分段填充处理的方法将壳体内部的空气用TNT替换。空气和TNT模型网格化时，由于榴弹壳体完全破裂后，网格疏密程度对破片在空气场中的飞行速度影响较小，因此对空气场模型网格采取由心至两侧稀疏化处理。通过测试，可以显著提高仿真速度。起爆点设置在装药最前端，起爆路径选择直接起爆。最终建立的网格化仿真模型如图4所示。

图4 网格化仿真模型

根据建立的仿真模型，采用AUTODYN软件在工作站上进行数值模拟，经过286 229个循环，战斗部壳体完全破碎为934个破片。由于是四分之一模型，根据对称原理，可得该榴弹战斗部共计产生3 736个破片。图5反映了榴弹爆炸时壳体的变化过程。TNT引爆后，榴弹前端首先出现变形，50 μs时，榴弹前端已明显膨大，一较大的破片脱离壳体。随着时间增加，爆轰波由前端向后端移动，其释放的爆轰压力使钢质壳体发生膨胀，榴弹直径明显扩大。壳体表面出现轴向轻微裂纹并不断增大，而后出现径向裂纹直至完全破裂，爆轰产物从裂纹向外流动并喷出。当气体产物开始衰减时，大小不一的破片就从爆炸云团中飞出。破片形成过程的仿真结果与实际爆炸过程基本一致[12]。

图5 榴弹壳体膨胀破裂过程

3 对比分析

3.1 速度对比分析

图6表示榴弹壳体完全破碎后，各破片的速度分布情况。由图可以看出破片初速主要分布在800 m/s～1 300 m/s之间，速度相对较低的破片考虑是壳体前后两段及壳体壁厚较大处产生。根据仿真模拟结果，忽略由壳体前后两端产生的质量明显过大的破片，所得破片的平均速度为1 041.34 m/s，二者速度差仅为44.66 m/s，误差为4.12%.由此可见，仿真运算同工程计算结果吻合度较高。

图6 破片速度分布

由于Gurney公式是假设炸药势能直接转变为金属材料的动能和爆轰产物的膨胀，破片初速是在理想状态下计算得到的，没有考虑爆轰产物溢出等情况，因此理论计算值比仿真结果偏大是合理的。破片初速的仿真结果与工程计算值在可允许的误差范围内，因此利用该模型来描述自然破片战斗部的破片飞散规律是可信的。

3.2 质量对比分析

由表4的仿真模拟结果可以看出，榴弹战斗部壳体完全破裂为大小、形状不同的破片后，不同质量区间内破片的数量分布情况。可以看出破片主要分布在0 g～8 g之间，当破片质量大于16 g时，其数量已明显减少。

表4 仿真模拟破片质量分布

由图7可知，当破片质量小于1 g时，工程计算同仿真结果误差较大，工程计算结果明显大于仿真模拟结果，考虑是由于部分破片质量过小，杀伤破坏意义可以忽略不计。当破片质量大于1 g时，二者吻合度较高，证明了仿真模拟与Mott公式对该型榴弹破片质量分布预测的参考意义。

图7 工程计算同仿真模拟破片数质量分布

4 结束语

本文应用ICEM CFD及AUTODYN软件对榴弹壳体膨胀、破裂、自然破片形成、飞散过程进行了仿真模拟计算，对破片的质量、速度分布情况进行了分析。通过运用工程算法进行理论分析，表明此方法可以清晰地反映榴弹爆炸过程，对破片速度、质量等参数具有较准确的预测效果，可以在一定程度上弥补榴弹实爆试验中破片参数不易测量、难以量化分析的缺陷。由于文中工程计算结果源自相关经验公式及材料参数，仿真模拟结果的准确程度有待进一步试验验证。

[1]安振涛，王超，甄建伟.常规弹药爆炸破片和冲击作用规律理论研究[J].爆破，2012，29（3）：15-16.

[2]纪兵兵，陈金瓶.ICEM网格划分与技术实例详解[M].北京：中国水利水电出版社，2012.

[3]许蕾，罗会信.基于ANASYS ICEM CFD和CFX数值仿真技术[J].机械工程师，2008，31（12）：65-66.

[4]石少卿，汪敏，孙波，等.AUTODYN工程动力分析应用实例[M].北京：中国建筑工业出版社，2011.

[5]姚武文.飞机战伤模式与机理[M].北京：航空工业出版社，2006.

[6]黄经伟.破片式战斗部破片规律及破片对步兵战车的侵彻作用研究[D].南京：南京理工大学，2014（10）：21-22.

[7]王纯，刘艳梅，周涛，等.基于ICEM CFD对汽轮机末级三围叶片流场网格划分方法的优化[J].汽轮机技术，2012，54（10）：324.

[8]尚伟，计欣华，杨晓婧.光弹实验在有限元建模中的应用[J].实验室研究与探索，2014，33（10）：16-17.

[9]AUTODYN Users Manual[M].California：Century Dynamics Corporation，2005.

[10]杨海军，张爱林.基于Von Mises应力的预应力钢结构拓扑优化设计[J].北京工业大学学报，2010（4）：475-477.

[11]沙云东，郭小鹏.随机声在和作用下的复杂薄壁结构Von Mises应力概论分布研究[J].震动与冲击，2011，30（01）：138

[12]王凤英，刘天生.毁伤理论与技术[M].北京：北京理工大学出版社，2009.

Simulation Study of Grenade Fragment Field Based on ICEM CFD and AUTODYN

LIU Li1，2，LU Dong-lin1，ZENG Yong-hua1，WANG Xiao-long1，BAI Lin-yuan1
（1.College of Field Engineer，PLA Univ.of Sic.&Tech.，Nanjing Jiangsu 210007，China；2.69241 Forces，Changji Xinjiang 831707，China）

The conventional way of studying grenade fragments is to collect and analyze the number quality of fragments by aerial static explosion or water well exploding device.In this paper，ICEM CFD preprocessing program and AUTODYN software were used to simulate the explosion process of a grenade.The velocity and mass distribution of grenade natural fragments were analyzed and verified by Gurney formula and Mott formula.The results show that the simulation is consistent with the actual explosion process，and the velocity and mass of the fragment field agree well with the engineering calculation results.

ICEM CFD；AUTODYN；fragment field

E932.4

A

1672-545X（2017）07-0060-04

2017-04-03

刘立（1986-），男，江苏邳州人，硕士研究生，主要研究方向为工程装备保障；鲁冬林（1969-），男，湖北荆州人，硕士，副教授，主要研究方向为工程装备保障。