蒙特卡洛剖分投影法的爆震弹破片平均比动能评估*

刘星雨，战仁军，欧阳的华，李永利

(1. 武警工程大学 装备管理与保障学院， 陕西 西安 710086； 2. 上海交通大学 海洋智能装备与系统教育部重点实验室， 上海 200240；3. 西安建筑科技大学 资源工程学院， 陕西 西安 710055)

爆震弹是利用含能材料爆炸时产生的强烈声响和眩目闪光，使有生目标暂时致聋致盲的非致命性弹种[1]。此类弹种采用手投、枪发或发射器方式进行发射，被广泛用于部队反恐、处突、海上维权执法等任务中。其封装壳体多为塑料或铝制材料，在爆炸过程中形成的较大破片易对人造成不可逆转的伤害。为提升其安全性，国内外学者目前主要采用的是改善壳体材料或外置高强度开孔金属壳的方法[2]。而对于采用壳体预刻槽方法控制破片的质量分布、飞散速度、比动能，进而减小破片场的致伤与致死半径的研究则较少。特别是基于上述指标下破片“成型、飞散、着靶”全时域的理论模型与仿真方法，还少有研究。为此，建立破片全时域的仿真模型，自动化地求解破片质量分布、速度、平均比动能、致伤半径、致死半径等安全性指标，是高效地指导武器使用、优化武器性能、节约研发成本的客观需求。

目前，国内外学者对杀伤性弹药破片的研究相对较多，特别是杀伤性预制破片战斗部威力场的数值仿真研究，提出了很多成熟的仿真模型与计算方法，如文献[3-10]基于LS-DYNA、AUTODYN平台对各类战斗部破片场进行了建模、仿真评估、效能优化、快速计算等。但对非致命弹药破片的安全特性研究，特别是全时域下的破片比动能研究还相对较少。上述研究中依托LS-DYNA、AUTODYN等通用有限元平台，通常只能获得爆炸初始时间段内的破片场飞散特性及冲击波压力，难以得到全时域全破片的平均比动能。而在建立全时域的飞散模型及平均比动能计算过程中，其重要参数平均迎风面积目前主要采用均匀取向理论的试验方法获得，而基于有限元结果的任意形状破片平均迎风面积求解还未见相关报道。

故本文根据爆震弹破片的比动能定义、形成机理及飞散原理，建立了平均比动能计算模型，并结合有限元壳体破碎模型及破片飞散模型，对同一装药下的自然破片和半预制破片全时域平均比动能及安全半径进行了评估。其中，平均比动能计算模型通过蒙特卡洛剖分投影(Monte Carlo Subdivision Projection Simulation，MCSPS)方法建立，可对任意形状有限元破片的平均迎风面积及比动能进行求解；破碎模型采用了相容流固耦合方法进行数值模拟；飞散模型是经有限元质量损耗修正后的质点外弹道模型。通过模型求解，可得到距爆心任意位置上的破片垂直靶分布、平均比动能阈值，进而对自然破片与半预制破片的致伤、致死半径进行了评估。

1 模型

1.1 基于MCSPS法的平均比动能计算模型

1.1.1 计算模型

根据破片的平均比动能定义，其表达式[11]为：

(1)

MCSPS方法通过对有限元破片的坐标点集进行N次平移、随机旋转、投影、三角剖分、凹域边界搜索排序、叉乘法求面积并取平均，可得到N次投影面积Aij的平均值。当N值趋于无穷大时，所得均值结果可近似为该破片的平均迎风面积。即

(2)

各计算步骤如式(3)～(9)所示。

(3)

(4)

(5)

(6)

(7)

(8)

(9)

其中：f1，f2，f3分别为对坐标向量的平移、旋转和投影运算；f4为平面坐标向量的Delaunay凸包三角剖分[12]，P为三角剖分凸包点集矩阵，Clist为凸包三角形矩阵；f5为带最大边长约束的三角剖分凹域边界搜索与排序，Pb为排序后的边界顶点向量；f6为叉乘法求多边形面积，Aij为第i枚破片第j次MCSPS后的迎风面积；f7为进行N次MCSPS后所得迎风面积取平均，N为MCSPS总次数。

基于MCSPS的平均迎风面积仿真算法流程如图1所示。其中，Pe为破片单元数据库，Pf为破片汇总数据库，i2、i分为Pe库行搜索标记的首端和末端，q为破片序号，mtkl为蒙特卡洛次数变量，N为蒙特卡洛总次数，设置为50 000次，area_mtkl为蒙特卡洛投影面积矩阵。

图1 基于MCSPS的平均迎风面积仿真算法流程Fig.1 Flow chart of simulation algorithm of average windward area based on MCSPS

1.1.2 凹域三角剖分边界搜索排序

对于投影点集为凸包的情况，可直接根据计算机图形学原理，采用Delaunay三角剖分进行求解。而对于凹多边形围成图形区域(简称凹域)，则无法直接进行求解。提出采用最大边长约束的原则，将凸包剖分算法所得大于约束值R的三角边去除，对凹域边界重新进行搜索，并用叉乘法求解凹多边形面积，可大大提高平均迎风面积的精度及计算效率。R取值为：

(10)

式中：R为边长约束值；λscale为放大因子，取8；li为剖分三角形第i条边的长度；ne为三角边总数。

带最大边长约束的凹域三角剖分边界搜索排序算法流程如图2所示。图2中：T为无重复节点的投影点集坐标矩阵；dm为凸包剖分边界矩阵；dt_Clist为凸包三角剖分三角形节点编号矩阵；edges为凸包三角剖分升序排列的边界矩阵；I为edges中不包含重复边的行号向量；IN为edges的行号向量；IO为边界边行号向量；sort_vertices为用于存放边界顶点序列的向量；polygon为按顺(逆)时针排序的边界边节点坐标矩阵。

图2 带最大边长约束的三角剖分边界搜索算法流程图Fig.2 Flow chart of boundary search-sort algorithm with maximum edge length constrain of triangulation

凹域三角剖分边界搜索的差集运算表达式为：

EIO=EI-EIN-I

(11)

式中：EIO为凹多边形区域边界边集合，有且仅有一条公共边；EIN-I为凹多边形区域的内部边集合，存在两个公共边；EI为不包含重复边的凹多边形边集合。

通过搜索算法得到边界顶点坐标矩阵polygon，采用叉乘法求解多边形面积可得到Aij，将结果代入式(2)和式(1)可得到破片的平均迎风面积及平均比动能。

1.2 壳体有限元破碎模型

1.2.1 有限元模型

某型爆震弹其结构原理如图3(a)所示。全弹长度为110 mm(不含击发机构)，直径为38 mm；主装药为一种含铝烟火药，质量为60 g；壳体为ABS塑料材料。对击发机构及壳体结合部作适当简化，无刻槽与内刻槽壳体几何模型如图3(b)、图3(c)所示。其中，1/4内刻槽壳体采用4条周向槽及2条轴向槽进行划分，槽深2 mm，槽宽2 mm。

图3 爆震弹结构原理及壳体几何模型Fig.3 Structure principle and shell geometry model

通过相容流固耦合方法对壳体破碎过程进行1/4对称建模，如图4所示，爆震剂及空气采用共节点Euler单元，连接座、壳体采用共节点的Lagrange单元。在对称边界面设置位移约束条件；在空气外表面上设置无反射边界。为了验证爆震剂参数的可靠性，建立6.5 cm×20 cm×306 cm的空气单元网格，用于计算3 m内的冲击波压力，以便与试验进行对比，网格大小为0.1 cm×0.1 cm×0.2 cm。为了确保Euler单元应力应变输运的精度，爆震剂与空气单元网格采用相同大小进行划分。为了提高壳体失效过程中的应力应变输运精度，将靠近壳体附近单元网格尺寸采用0.1 cm×0.1 cm×0.1 cm进行划分。

图4 爆震弹流固耦合有限元模型Fig.4 Finite element model of fluid-solid coupling

弹体采用立式放置，炸高距地面1 m。爆震剂采用单点起爆，炸点位置距连接座下端面1.5 cm。

1.2.2 材料本构及失效准则

爆震剂：采用HIGH_EXPLOSIVE_BURN高能炸药本构，其琼斯-威尔金斯-李(Jones-Wilkins-Lee, JWL)状态方程为：

(12)

式中：p为压力；Vg为爆轰产物体积和炸药初始体积之比；E0为爆轰产物单位体积的初始内能；A，B，R1，R2，ω为特征参数。

密度ρ根据装药质量及体积计算得到，爆速D根据文献[1]得到，JWL方程各参数借鉴含铝炸药文献[13]，并作以适当调整，可使1 m、2 m、3 m处的冲击波超压与试验基本吻合。各参数取值见表1。

空气：采用Null材料本构和线性多项式状态方程Linear-Polynominal进行建立，即

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

(13)

式中，C0～C6为与气体性质有关的状态方程参数，C0、C1、C2、C3为压强量纲，C4、C5、C6无量纲。地面空气μ=0，排除大气相对压力，由此得出的空气材料模型参数见表2。

表1 爆震剂材料模型参数

表 2 空气材料模型参数

壳体：采用PLASTIC_KINEMATIC塑性随动材料本构，该模型应变率采用了Cowper-Symonds模型。动态屈服应力表达式为：

(14)

(15)

式中，Etan为切线模量，E为弹性模量。其材料模型参数见表3[2]，ρ为材料密度，V为泊松比，σ0为屈服极限，FS为最大有效塑性应变。

表3 壳体材料模型参数

LS-DYNA中add_erosion关键字提供了最大主应变失效、最大主应力失效、等效塑性应变失效、最大静水拉压失效等失效准则。经过多次试算对比，并结合材料本构中最大有效塑性应变FS取值，发现采用最大主应变失效准则可较好地解决壳体裂纹附近单元出现大面积失效的情况。

1.3 破片飞散模型

当壳体完全破碎后，冲击波对破片的速度影响逐渐减弱，直至某一时刻速度不再显著增加。此后，各破片在空气阻力FR、自身重力G及旋转力矩M作用下做自由飞行。根据质点外弹道学原理[14]并考虑有限单元失效法中质量、平均迎风面积的损耗因素，引入破碎前壳体质量mr进行修正。可推得经修正后的第i枚破片在二维平面内速度vi、水平夹角θi、水平位移xi、垂直位移yi和比动能edi随时间t变化的常微分方程组。

(16)

设置迭代终止条件为yin≤0或tin≥Tmax(Tmax为时间阈值，取10 s)。迭代初始值yi0从破碎模型中直接读取，vi0、θi0、xi0需从破碎模型当前帧及前一帧破片单元的速度、坐标分量数据中进行计算，其表达式为：

(17)

根据空间坐标矩阵变换关系，将各破片二维轨迹矩阵[xi,yi]按各破片初始方位角绕坐标轴y旋转，可获得各破片的三维飞散轨迹矩阵[x′i,y′i,z′i]以及垂直靶分布情况。

2 结果分析与讨论

2.1 冲击波超压对比

为了验证爆震剂参数的可靠性，进行自由场冲击波超压测试试验，测得了距爆心1 m、2 m、3 m处的冲击波超压峰值并与仿真超压进行了对比，不同距离上的仿真与试验结果对比如图5所示。图5中，无刻槽与内刻槽壳体在1 m、2 m、3 m处的冲击波超压峰值仿真结果基本相同；两类壳体在1 m、2 m处的超压峰值与试验结果基本吻合，但在3 m处试验结果比仿真结果更低。3 m处仿真结果与试验结果产生误差的原因主要有两个方面：一是JWL状态方程及高能炸药本构模型主要是基于TNT炸药而建立，对于含铝烟火药的爆炸数值模拟会存在一定误差；二是LS-DYNA等商用有限元软件无反射边界面附近仍会存在冲击波超压反射的情况，使得边界附近单元压力值偏大。

2.2 质量分布对比

通过破碎模型得到的未刻槽、内刻槽壳体破碎情况如图6(a)、图6(b)所示。两类壳体在40 μs 时刻已完全破碎。根据飞散方向不同，将破片分为顶部破片、底部破片和径向破片，并进行编号，如图6(c)、图6(d)。从图6可知，内刻槽半预制破片比无刻槽自然破片数量多，径向破片尺寸更小，质量分布更加集中。

(a) 未刻槽壳体破碎情况(a) Breaking situation of non-grooved shell fragments

通过同一破片搜索排序算法，可得到以破片为单位的数量、质量、单元数、体积等。其质量频数分布情况如图7所示。图7中，两类破片质量分布呈现指数衰减的规律，大多数破片主要集中于0～1 g之间；除顶部击发座破片超过3 g外，半预制破片1～3 g的破片数比自然破片更少；破片均值和标准差更小，质量分布更加集中。

图7 仿真破片质量分布折线图Fig.7 Line graph of mass-frequency distribution of simulation fragments

两类破片的质量统计见表4。表4中，半预制破片质量最大值、极差、均值、标准差均比自然破片小。

表4 破片质量统计

综上可知，两类壳体结构的不同导致了应力集中点位置的不同，使得各破片单元失效的扩展方向产生差异，进而影响了破片的形状及质量。仅从质量因素上看，半预制破片可在一定程度上降低破片比动能，但具体情况还需和破片速度及平均迎风面积一起进行分析。

另外，通过破片回收试验，得到了同一壳体参数下的自然破片与半预制破片的质量分布情况，如图8所示。从质量分布上看，两类壳体的试验破片质量分布情况与仿真破片结果较为相似，绝大多数破片质量在0～1 g区间内，且两类破片质量呈指数衰减规律。但从破片形状上看，内刻槽仿真破片与试验破片形状较为吻合，但无刻槽仿真破片No.7、No.9、No.11、No.13与试验破片形状仍存在一定差异。

图8 试验破片质量分布折线图Fig.8 Line graph of mass-frequency distribution of test fragments

自然破片仿真结果与试验结果存在差异的原因是塑形随动本构模型是基于唯象理论并结合单元失效法出的破片失效与断裂，难以真实模拟出材料微孔洞及应力分布不均匀导致的壳体裂纹扩展的不确定性，使得仿真自然破片与真实破片的形状、质量分布与试验结果存在一定的差异。在后续研究中，可引入细观尺度下的弱点随机分布因子，对现有本构进行二次开发，以使仿真自然破片与真实破片的形状、质量分布更加吻合。

2.3 初始速度对比

根据破碎模型结果，通过LS-PREPOST后处理平台对顶部、底部、径向的各破片建立跟踪单元，其初始速度的时间历程如图9所示。图9中：两类破片在50 μs之后速度逐渐趋于稳定，不再显著增加；底部方向破片单元受旋转力矩影响较大，速度曲线呈现周期振荡，需采用中值滤波算法对其进行平滑处理后再进行数据读取。

(a) 自然破片初始速度(a) Initial velocity of natural fragments

不同方向上的破片初始速度统计见表5。从表5可知，两类壳体底部破片速度均大于径向破片速度，而径向破片速度又大于顶部破片速度。造成不同方向上速度差异的原因主要有两个方面：一是顶部破片的卸载效应使其方向上的应力波显著降低；二是底部端面距起爆点较远，应力波输运至此处时会在端面附近产生叠加，使应力波显著增大。

表5 不同方向的破片初始速度统计

其中，顶部方向半预制破片与自然破片速度差异不大，而径向、底部方向半预制破片的速度均值、标准差均比自然破片更低。仅从初始速度上看，半预制破片速度阈值更低且更加集中，可在一定程度上降低破片比动能。

2.4 平均迎风面积等初始参数

将破碎模型有限元数据导出，通过MCSPS平均迎风面积仿真程序计算，可得两类破片的数量、质心坐标、质心速度、单元数、平均迎风面积等初始参数，其中各破片迎风面积数据见表6。表6中，自然破片4、6、8、10、14、15单元数不足3个，节点投影后会出现无法构建三角剖分的超级三角形情况。对于不足3个节点的破片飞屑采用了近似计算的方法。将表6中各破片初始参数代入飞散模型中进行求解，可得到各破片平均比动能随距离的衰减情况。

表6 MCSPS算法下的自然与半预制破片初始参数

为验证MCSPS算法求迎风面积的可靠性，选择长方体、圆柱体预制破片进行求解并与理论公式[15]结果进行对比。设长方体尺寸为1 cm×2 cm×4 cm，网格尺寸为2 mm×2 mm×2 mm；圆柱体为半径r=1 cm，高度h=4 cm，网格尺寸为2 mm×2 mm×1 mm。当循环次数N=16 384(214)时，MCSPS算法所得长方体破片、圆柱体破片迎风面积结果与理论公式对比见表7。从表7中可知， MCSPS算法结果与理论公式结果基本吻合，其中长方体破片误差仅为0.12%，而圆柱体破片误差为5.29%。

表7 MCSPS算法与理论公式平均迎风面积误差对比

MCSPS算法与理论公式结果存在误差的原因主要有两点：一是伪随机数误差，即由软件产生的均匀分布样本是基于线性同余算法建立的，与真实均匀分布样本不可避免地存在一定误差，但此类误差在16 384次循环下对结果影响较小，如长方体破片的误差结果；二是剖分精度误差，即有限元网格是一种对真实物体的剖分重构，网格尺寸的大小会直接影响物体几何特征的精度，特别是对于圆柱体这类具有曲线、曲面的几何体而言，网格尺寸过大会导致圆柱体被分割成正多棱柱体，进而所得平均迎风面积将比真实的理论值要小，如圆柱体的误差结果，达到了5.29%。

2.5 平均比动能及安全半径对比

根据表6数据及飞散模型，得到了与爆心不同距离的破片垂直靶分布和平均比动能阈值，其中两类破片1.5 m、5 m处垂直靶分布如图10所示。从分布上看，两类破片飞散角相差不大，分布均呈现不完全对称，其中半预制破片比自然破片横向分布更稀疏、纵向分布更密集。从比动能数值上看，在1.5 m处自然破片的比动能阈值比半预制破片阈值更大，其径向破片6、7、8、10、13、14均超过半预制破片的比动能阈值；而在5 m处其比动能阈值反而远低于半预制破片的比动能阈值。造成此现象的原因主要有两点：一是自然破片7、13的尺寸较大且形状更加不规则，单位质量的平均迎风面积相比于半预制破片更大，飞散过程中速度衰减更快，进而使得自然破片比动能衰减得更快；二是由于自然破片6、8、10、14初始射角为负，在不足5 m时就已着地。

(a) 1.5 m处的自然破片(a) Natural fragments at 1.5 m

两类破片1 m至5 m内的比动能最大值、最小值、均值、标准差等结果见表8。在小于等于2.5 m范围内，自然破片比动能最大值、均值均高于半预制破片；超过2.5 m后，半预制破片比动能最大值反而高于自然破片。

表8 不同距离全破片的平均比动能评估结果

根据弹丸、破片、小箭杀伤人员的最小比动能标准160 J/cm2以及擦伤皮肤的最小比动能标准9.8 J/cm2[11]，对表7数据进行插值，得到两类破片致死半径和致伤半径，见表9。从表9中可知，相比于自然破片，半预制破片的致死半径有显著减小，但致伤半径并未有显著改善。根据非致命防暴弹的安全性要求，壳体破片比动能不宜高于擦伤皮肤标准，即自然破片、半预制破片安全半径分别为4.7 m、4.6 m。

表9 自然与半预制破片致死、致伤半径

3 结论

1) 基于MCSPS方法建立的平均比动能计算模型，可求解任意形状的有限元破片平均迎风面积和平均比动能。

2)全时域下的破片比动能评估相对于初始时刻的比动能评估更全面。

3)相比于自然破片，采用半预制破片方案可显著减小爆震弹破片的致死半径，但并不能对其安全半径有显著改善。要进一步降低其致伤半径，还应考虑起爆点位置、壳体厚度、材料等因素。