【基础理论与应用研究】
浅埋炸药冲击试验台架数值模拟方法研究
现代战争中防护型车辆时常面临着简易爆炸装置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粒子算法中考虑到土壤与结构的耦合效应，土壤冲击时域上晚于冲击波，所以在试件经历了冲击波之后，土壤继续对试件产生能量传递，造成泡沫铝吸收能量继续增长。
图15 ALE计算结果
图16 土壤能量曲线
图17 背板能量曲线
图18 泡沫铝内能曲线
图19 台架运动量曲线
4) 台架运动量上对比中可以看出，PBM-DEM粒子算法台架最大跳动量为68.2 mm，ALE算法台架最大跳动量为64.1 mm，时域上峰值出现时间PBM-DEM粒子算法较ALE算法滞后。
5) 土壤能量曲线对比中可以看出，峰值上差异性在4%，总体趋势类似，但是由于ALE算法在土壤飞溅超出计算域时会被删除，而PBM-DEM粒子算法中土壤粒子不会被删除，只是由于粒子与结构碰撞引起土壤能量降低，所以其土壤剩余能量大于ALE算法。
6) 两种算法计算相同的时长上，PBM-DEM粒子算法需要的时间更短，计算效率更高。
4 台架试验
4.1 试验内容
本次台架试验目的是测试某型钢板-泡沫铝复合材料抗爆炸冲击性能，同时验证前期爆炸仿真精度。
根据试验要求制作台架如图20所示。外廓支架均采用10 mm厚的Q235钢材制造；试验台架由四根圆柱导杆及Q235实心扁钢制作；配重块按照要求配至1 t；试件两端采用M24螺栓连接在活动支架上。
根据前期仿真边界，将试件迎爆面调节至离地高300 mm处，地雷埋深100 mm，实现中心引爆。
图20 台架试验装置
图21-图22中显示了爆炸后台架及试件的整体状况。根据试验前对背板的网格处理，针对各点变形进行测量，最终拟合背板残余变形曲面，如图23所示。试验前在导杆上涂抹黄油，根据测量试验后黄油在导杆上的痕迹长度，来描述试验台架的垂向运动量，如图24所示。
图21 台架试验后状态
图22 试件试验后状态
图23 变形曲面拟合
图24 台架运动量测量
试验结果表明，背板最大残余变形量为29.3 mm，坐标为(200,250)，台架活动部件最大垂直位移为66.7 mm。
4.2 仿真试验对比分析
根据试验结果测出的各点变形量，拟合变形曲面，发现最大变形区域偏向于试件无约束段，有螺栓约束段变形较小。对比两种算法计算结果与试验实际测量结果如表6所示。
表6 仿真结果与试验结果
PBM-DEM粒子算法计算结果最大变形点与试验值差别在30 mm左右，最大变形量误差为3%，台架运动量差别为2.5%；ALE算法计算结果最大变形点与试验值差别在133 mm左右，最大变形量误差为3%，台架运动量差别在4%；对比分析得出背板变形情况与PBM-DEM粒子算法计算结果更接近。
5 结论
1) 采用了PBM-DEM粒子算法模拟了浅埋地雷爆炸冲击试验台架试验，分析了DEM粒子半径对仿真计算的影响，发现粒子半径和耦合结构网格尺寸一致时，DEM粒子与结构耦合效果较好，且未产生粒子渗透的现象。
2) 采用ALE算法模拟了浅埋地雷爆炸冲击试验台架试验，分析了两种不同算法计算结果中的背板残余变形量、背板的总能量、泡沫铝内能、土壤动能，台架运动量及计算时长。发现两种算法在前15 ms内的计算结果吻合较好，15 ms后由于PBM-DEM粒子算法中土壤粒子与结构耦合效应，导致结果差生差异。
3) 和ALE算法相比，PBM-DEM算法仿真计算结果与实爆试验结果更接近，拟合度更高，且计算效率更高。
[1] 彭兵,王显会,王磊,等.某轻型车辆底部结构爆炸仿真与试验研究[J].科学技术与工程,2017(12):178-183.
[2] 魏然.基于爆炸防护的车身结构多学科优化研究[D].南京：南京理工大学,2017.
[3] 李利莎,谢清粮,郑全平,等.基于langrange、ALE 和SPH算法的接触爆炸模拟计算[J].爆破,2011(3):18-22.
[4] TRAJKOVSKI J,KUNC R,PREBIL I.Blast response of centrally and eccentrically loaded flat-,U-,and V-shaped armored plates comparative study[J].Shock Waves,2017,27(4):583-591.
[5] OLOVSSONA L,HANSSENB C A G,B?RVIKC C T，et al.A particle-based approach to close-range blast loading[J].European Journal of Mechanics A/Solids,2010(29):1-6.
[6] TENG H,WANG J.Particle blast method (PBM) for the simulation of blast loading[C]//Paper Presented at the 13th International LS-DYNA Conference,Detroit,USA.2014.
[7] B?RVIK T,OLOVSSON L,HANSSEN A G,et al.A discrete particle approach to simulate the combined effect of blast and sand impact loading of steel plates[J].Mech.Phys.Solids,2011,59(5):940-958.
[8] JENSEN M R,SMITH W.Discrete particle method is a predictive tool for simulation of mine blast—a parameter study of the process and approach[C]//Paper Presented at the Proceedings of the 2015 NDIA Ground Vehicle Systems Engineering and Technology Symposium (GVSETS),Michigan.2015.
[9] BARANOWSKI P，MALACHOWSKI J.Numerical study of selected military vehicle chassis subjected to blast loading in terms of tire strength improving[J].Bull.Pol.Acad.Sci.Tech.Sci,2013,63(4):867-878.
[10] GUO Q,ZHOU Y,WANG X,et al.Numerical simulations and experimental analysis of a vehicle cabin and its occupants subjected to a mine blast[J].Proc.Inst.Mech.Eng.D.Automob.Eng,2016,230(5):623-631.
[11] KARAJAN N,HAN Z,TEN H,et al.Interaction possibilities of bonded and loose particles in LS-DYNA[C]//Paper Presented at the 9th European LS-DYNA Conference,Manchester,UK.2013.
[12] HALLQUIST J.LS-DYNA Keyword User’s Manual[M].Version:R10.0.Livermore,California.2017.
[13] 王飞,陈卫东.爆炸冲击载荷作用下板壳结构数值仿真分析[J].强度与环境,2010,37(4):37-39.
[14] 曾爱,王显会,周云波.防爆油箱爆炸冲击响应数值模拟方法研究[J],兵器装备工程学报,2019(02):55-59.
Research on Numerical Simulation Method of Impact Test Bench for Shallow Buried Explosives
Citation format:WU Mengyang, ZHOU Yunbo, ZHAO Feng, et al.Research on Numerical Simulation Method of Impact Test Bench for Shallow Buried Explosives[J].Journal of Ordnance Equipment Engineering,2020,41(10):117-124.