浅埋炸药冲击试验台架数值模拟方法研究

2020-11-05吴梦阳周云波王利辉李进军

兵器装备工程学报 2020年10期

吴梦阳，周云波，赵丰，王利辉，李进军，张明，周迪

(1.南京理工大学机械工程学院，南京 210094；2.陆军某项目管理中心，北京 100073；3.中国人民解放军32379部队，北京 100073； 4.中国人民解放军32381部队，北京 100073)

现代战争中防护型车辆时常面临着简易爆炸装置IED(Improvise Explosive Device)的威胁，其巨大的爆炸冲击有可能击穿车辆底部，给车辆及车内乘员带来毁灭性的损伤。因此，近年来针对军用防护型车辆底部抗冲击技术的研究逐步发展，通过分析爆炸冲击应力波在车辆结构及乘员系统中的传播路径，逐层降低应力波强度，以降低乘员损伤指标作为最终目的，系统化提升整车抗爆炸冲击能力，致力于降低车辆及乘员的损伤[1]。

自北约STANAG 4569标准发布后，各国逐步开展车辆抗爆炸试验研究，试验技术也逐渐走向成熟。整车爆炸试验是评估车辆抗爆炸冲击能力最直接的手段，为车辆抗爆炸防护设计提供最有效的素材。但是，由于整车爆炸试验周期长、可重复性不高等原因，利用CAE仿真技术，结合试验结果，可以有效预测爆炸环境下车辆及乘员的损伤效应，为后期车辆防护能力优化设计提供方向[2]。

近年来，用任意拉格朗日-欧拉算法(ALE)和流固耦合算法模拟爆炸源、冲击波传递以及结构响应是公认的爆炸数值模拟算法之一。但是由于ALE算法需要建立大量的空气和土壤单元，造成计算效率不高，并且欧拉算法与结构耦合的鲁棒性较低，可能会造成计算结果不稳定。

在此基础上，基于粒子的无网格法逐渐被开发，如SPH、PBM(Particle-Blast Method)、DEM(Discrete-Element Method)等，用来协调计算的鲁棒性。李利莎[3]使用Lagrangian算法、ALE算法和SPH算法分别模拟了炸药冲击混凝土过程，对比了三种不同算法的优劣性。Trajkovski[4]使用了ALE算法和SPH算法模拟爆炸冲击靶板，比较两种算法的差异，发现SPH算法可有效避免冲击波渗漏问题。Olovsson等[5]首次使用PBM算法进行爆炸模拟，之后，滕海龙[6]对算法进行了改进，使用PBM算法进行了一系列的近场爆炸模拟计算，计算结果与试验结果具有较好的一致性，同时发现PBM算法在计算经济性上优于ALE方法。Børvik等[7]用离散单元法DEM模拟了土壤介质，用来研究球形炸药在土壤中的爆炸特性。Jensen等[8]使用PBM-DEM粒子算法对地雷爆炸进行了模拟研究，旨在了解不同土壤密度、爆炸物形状等变量下车体的结构响应。

本文主要研究内容是根据某浅埋爆炸试验台架试验要求，分别使用PBM-DEM爆炸算法以及ALE算法进行仿真分析；比较了两种算法结果中试件背板最大残余变形量、台架的运动量、试件背板的总能量、泡沫铝内能以及土壤动能；最终进行实爆试验，验证两种算法的精度，对比分析得出PBM-DEM粒子算法计算结果与试验结果更接近，其计算时间也比ALE算法计算时间短。

1 PBM-DEM数值模拟算法理论

1.1 PBM粒子爆炸算法

PBM粒子算法是针对基于分子动力学理论KMT(Kinetic Molecular Theory)的微粒法CPM(Corpuscular Method)进行了改进。CPM主要是用于模拟气囊的爆炸驱动，通过将大量的分子整合为一个粒子，整合比例接近1018∶1。CPM避免了时间尺度与长度尺度的矛盾性，但CPM不适用于热平衡假设失效的爆炸模拟。

CPM是建立在KMT基础上的一种粒子法，允许对气体体积进行数值处理时存在一定偏差。在微观层面上动力学分子理论表现为研究气体分子及其相互作用，宏观层面上表现为理想气体定律。该理论基于以下假设[9]：

1) 分子之间的平均距离比它们自身尺寸要大。

2) 存在热力学平衡，即分子处于随机运动状态。

3) 分子遵循牛顿运动定律。

4) 分子-分子和分子-结构的相互作用是完全弹性的。

动力学分子理论最早可追溯于1738年，当时Daniel Bernoulli提出了一个理论，即空气对活塞的压力是通过离散的分子碰撞建立起来的。以动力学理论为出发点，James Clerk Maxwell导出了热平衡时分子速度分布的一个非常精确的表达式[10]:

式中:f(v)是分子速度的概率密度函数;M是摩尔质量;R是通用气体常数；T是温度。

使用CPM做爆炸模拟时，每克爆轰产物的分子量大约在1023～1024之间。所以，在整个爆炸过程中，不可能做到对每个分子进行模拟。在权衡了数值计算精度要求和计算效率之后，针对分子动力学基础的微粒法进行了以下假设[11]：

1) 粒子均被定义为刚性球体；

2) 每个粒子包含了1015～1020个分子，具体数值取决于实际情况；

3) 对于每个单独的粒子，平动动能Wt和转动动能Ws之间是平衡的。

但在实际计算过程中，需要注意真实分子云的平动能和转动能的比值是一个统计值，不是对每个分子都有效。在热平衡假设的基础上，确定了Wt和Ws比例。对于大多数安全气囊应用来说，这是一个有效的假设，但在气流速度极高的爆炸模拟中则不是这样。故PBM算法针对CPM假设进行如下补充[12]：

1)Ws是每个粒子的集中标量变量。因此，粒子没有被赋予任何旋转/振动的自由度。

2) 当粒子撞击结构时，Wt值将发生变化，而Ws假定不受影响。因此，Ws和Wt的比值发生变化。

3) 当两个粒子碰撞时，它们的能量被重新分配，以恢复或保持Ws和Wt之间的正确比例。这样做的方式确保了总能量平衡和动量平衡。

1.2 DEM离散单元法

在模型前处理过程中，DEM方法是将模型离散为刚性球体，并通过罚函数定义接触。DEM是根据牛顿运动定律分别计算每个粒子的运动，采用六自由度系统来描述运动。图1描绘了DEM粒子间相互作用，其中包括法向和切向的刚度和阻尼力、静态摩擦力、滚动摩擦力和表面张力。由dint定义的两粒子相互作用距离。

dint=|x1-x2|-(r1+r2)

其中:xn是球体中心的坐标矢量;rn是球形粒子的半径。

图1 粒子相互作用示意图

粒子碰撞状态分为3种，如图2所示，当dint小于0时，表面张力和机械接触都存在，这可能会导致粒子之间相互渗透；当两个粒子不直接接触时，只存在表面张力；当两个dint足够大时，粒子间不再存在相互作用力。

图2 粒子不同接触状态示意图

法向弹簧力Fkn表达式为：

Fkn=Kndint

式中:Kn为法向接触刚度；KNorm为法向弹簧常数的比例因子；kn为根据弹性模量En和泊松比vn计算得出的压缩模量。

切向刚度表达式为：

Kt=KnKShear

法向阻尼力Fdn决定了两个粒子之间的阻尼特性，表达式为：

Fdn=Dndint

式中:mn为质点质量；Dn定义法向临界阻尼比；dnorm为法向阻尼系数。正像法向阻尼力一样，切向阻尼力由切向阻尼系数dtan定义。

2 仿真计算

为了验证某钢板复合泡沫铝夹心结构的抗爆炸性能，设计了一款试验台架，三维模型如图3所示。台架由外廓支架、试验台架以及配重块组成，均由Q235制造。试验台架部分与外廓支架间设立四根滑动导杆，以确保其可以垂向运动。试件安装在试验台架底部，两端由M24螺栓固定。

试件为钢板复合泡沫铝三明治夹芯结构，总体尺寸为600 mm×600 mm×40 mm。面板及背板为6 mm防弹钢板，芯层结构为28 mm泡沫铝，钢板与泡沫铝使用胶结。

图3 试验台架三维模型示意图

LS-DYNA是由Livermore Software Technology Corporation(LSTC)开发的高级通用多物理场仿真软件包。该软件包被用来计算许多复杂的现实世界问题，其起源和核心竞争力在于使用显式时间积分的高度非线性瞬态动态有限元分析(FEA)。LS-DYNA多用于汽车，航空航天，建筑和土木工程，军事，制造和生物工程行业。

2.1 PBM-DEM仿真计算

根据试验台架模型建立1∶1爆炸仿真模型，如图4所示。台架模型划分拉格朗日六面体实体单元网格，网格单元尺寸为10 mm，台架材料使用选择关键字Mat_Plastic_Motional来定义；试件建立的是5 mm尺寸的实体单元，由于爆炸冲击的响应时间短，需要考虑材料的应变率效应，因此，防弹钢材料使用Mat_Johnson_Cook材料模型来描述，防弹钢材料参数进行霍普金森杆进行动态拉伸试验获得，拉伸试样如图5所示；泡沫铝材料模型使用Mat_Crushable_Foam来表示，压缩曲线由准静态压缩试验获得，准静态压缩试验如图6所示。钢材料模型的主要参数如表1、表2所示(由实验数据拟合计算得出)。防弹钢应力应变曲线如图7，泡沫铝应力应变曲线如图8。

根据试验内容要求，试件迎爆面离地高度为300 mm，炸药埋深100 mm。根据模型坐标，计算得出炸药位置，通过Define_Pblast_Geometry关键字定义炸药形状特征，并使用Particel_Blast关键字定义TNT材料模型，其定义的爆轰波参数如爆速、初始内能以及密度均来源于JWL状态方程，同时该关键字还定义了爆炸粒子与DEM粒子和结构的耦合。

图4 PBM-DEM模型

图5 霍普金森拉伸试样

图6 准静态压缩试验

图7 防弹钢应力应变曲线

图8 泡沫铝应力应变曲线

表1 Q235材料参数

表2 防弹钢材料参数

根据前期大量仿真经验发现，DEM粒子半径需要保持在合理的范围区间内。半径过大会导致粒子携带动能过大，造成击破模型现象；半径过小会导致粒子容易穿透拉格朗日网格从而造成渗漏，同时粒子数增加，影响计算效率。一般设置粒子半径在耦合的拉格朗日网格尺寸的1～1.5倍，由于耦合的试件网格尺寸为5 mm，故将土壤粒子半径设置为5/6/7/8 mm分别进行计算。通过Disc_Sphere_Generation关键字构建DEM粒子。使用Mat_Rigid材料模型定义DEM粒子，粒子间相互作用关系由关键字Define_discrete定义，粒子与结构间的相互作用关系由Define_De_To_Surface_Coupling定义，部分参数见表3[8]。

表3 Define_De_To_Surface_Coupling参数

2.2 ALE仿真计算

为了与PBM-DEM模型比较，根据试验内容要求同时建立了基于ALE算法的台架有限元模型，模型如图9所示。数值流场长为2 000 mm，空气场高1 800 mm，土壤场深400 mm，流场的网格尺寸为15 mm。空气场材料采用Mat_Null模型与线性多项式状态方程定义[13]：

P=C1+C2μ+C3μ2+(C4+C5μ+C6μ2)E

P表示气体压力，Ci(i=0、1、2、3、4、5、6)表示多项式系数，E表示单位体积空气的初始内能。C0=-0.1，C1=C2=C3=C6=0，C4=C5=0.4；E0=0.253 mJ/mm3；V0=1[14]。土壤的超塑性材料模型由Mat_Soil_And_Foamble_Fail定义。

采用关键词Mat_High_Explosive_Burn来描述TNT的材料模型。爆炸产物的压力、体积和内能之间的关系由EOS_JWL定义。

式中ω、AJWL、BJWL、R1和R2是与爆炸材料有关的参数;V是相对体积；E0初始内能。TNT的材料性能和JWL参数见表4[14]。

图9 ALE模型

表4 JWL参数

3 计算结果

3.1 PBM-DEM粒子算法计算结果

通过分析粒子半径在5/6/7/8 mm工况下的计算结果发现，7/8 mm粒子半径下依然出现击穿迎爆面的现象；6 mm粒子半径工况下未出现击穿情况，但是出现了网格不平滑现象；5 mm粒子半径工况下，迎爆面呈现了较平滑的耦合状态，同时未出现粒子穿透现象。具体如图10～图13所示。

图10 半径8 mm计算结果图11 半径7 mm计算结果

图12 半径6 mm计算结果图13 半径5 mm计算结果

当DEM土壤粒子半径为5 mm时，试件面板未出现被DEM粒子击穿现象，且粒子与试件有了较好的接触耦合，所以选用该工况仿真计算结果作为最终PBM-DEM粒子算法计算结果。

3.2 两种数值模拟方法比较

分析对比两种算法背板残余变形量X1、背板的总能量E1、泡沫铝内能E2、土壤动能Esoil，台架运动量X2及计算时长t。具体仿真工况结果见表5，其背板变形云图和有关曲线见图14～图19。

表5 仿真工况结果

图14 PBM-DEM计算结果

1) 背板最大残余变形量曲线对比可以看出，通过PBM-DEM粒子算法计算得出的背板最大残余变形量为28.5 mm，ALE算法计算得出的背板最大残余变形量为31.3 mm。

2) 背板总能量包括动能以及变形产生的内能，峰值上两种算法的差异在5%左右，但是最终能量保持阶段中，ALE算法中由于使用壳单元模拟背板，模型计算仍在震荡未平稳，所以残余能量较高，但是从曲线中可以看出残余能量正在缓慢下降，相比下PBM-DEM粒子算法中背板最终残余能量曲线已经基本水平。

3) 泡沫铝内能曲线对比中可以看出，当爆炸冲击来临时，泡沫铝急剧压缩吸收能量，内能曲线迅速升高，曲线拐点处能量差别在6%左右。由于ALE算法中忽略了土壤的动能效应，所以当爆炸冲击完毕后，ALE算法泡沫铝内能曲线基本水平；PBM-DEM粒子算法中考虑到土壤与结构的耦合效应，土壤冲击时域上晚于冲击波，所以在试件经历了冲击波之后，土壤继续对试件产生能量传递，造成泡沫铝吸收能量继续增长。