1 引言
军用地面车辆在战场上面临各式各样的威胁，如埋雷、简易爆炸装置(IED)和其他动能弹的威胁[1]。这些威胁通常携带装在金属外壳中的烈性炸药。当炸药爆炸时，产生的高压使金属外壳膨胀破碎，并形成自然破片以非常高的速度传播，最终产生的冲击波和高速自然破片撞击军用车辆，对车辆及乘员造成损伤。
整车实爆试验是验证车辆防护性能最直接有效的方法，但由于其危险性大、试验成本高、试验周期长、不可重复等原因，在进行车辆防护性能评估时，通常采用有限元仿真技术结合试验的方法，预测军用车辆在各式威胁下车辆的结构响应和车内乘员的损伤情况，并为后续车辆防护设计提供基础[2]。
在分析军用车辆抗爆炸冲击过程中，研究人员通常分析爆炸冲击波对车辆结构的损伤而忽略高速破片对车辆结构的破坏。而有研究表明爆炸冲击波与高速破片对于目标结构存在耦合作用，在耦合作用下目标结构的损伤情况要大于单独作用[3]。在爆炸数值仿真方法上，任意拉格朗日-欧拉算法(ALE)在模拟爆炸冲击波传递中有着较高的精度，被广泛应用在爆炸数值仿真中[4-6]。但该算法需要建立大量的空气和土壤网格，计算效率低，并且需要控制流固耦合参数防止流场泄漏，造成计算不稳定。随后基于无网格的光滑粒子法(Smoothed particle hydrodynamics，SPH)被提出，以解决爆炸环境下大变形造成网格畸形问题。胡建宇[7]通过ALE算法对某装甲车辆进行爆炸仿真分析，并对车辆底部梁结构进行优化；石秉良[8]采用SPH法模拟炸药爆炸对驾驶室底部结构响应，验证了SPH算法在处理车辆底部爆炸等问题的可行性。目前在研究破片类型上，主要分为预制破片[9]和自然破片[10]。预制破片分布比较规律，破片大小相似，无法模拟榴弹爆炸产生的大小不同的破片。本文研究的为榴弹爆炸产生的自然破片，相比预制破片更接近实际情况。
本研究首先通过爆炸冲击钢板台架试验，对比分析ALE算法和SPH算法中钢板的最大残余变形量、钢板动能与内能，并与试验结果对比，验证SPH法的精度。其次建立某军用车辆和大口径榴弹有限元模型，并对浅埋榴弹爆炸形成自然破片的过程进行仿真。最后以某军用车辆为例，分析了浅埋榴弹爆炸，冲击波和破片联合作用下的车身结构损伤情况。
2 光滑粒子算法(SPH)理论
在SPH法中，流体由一组有限的粒子表示，通过使用一种平滑函数，在区域Ω中，函数f(x)在x点处的值
近似为积分插值[11]：
(1)
式中：W(x-x′,h)为平滑函数，通常称为插值核；h表示核的搜索区域，在有限元分析中对应为网格大小的光滑长度；x为粒子所处的位置矢量。SPH算法中粒子搜索如图1所示。
图1 SPH粒子搜索示意图
Fig.1 Schematic diagram of SPH particle search
为每个粒子引入体积质量
，SPH公式的质量、动量和能量守恒方程的离散形式分别为[12]：
(2)
(3)
(4)
式中：Wij=W(xi-xj,h)；x表示直角坐标下粒子位置；vα 表示速度分量； ρ表示密度；t表示时间；σαβ表示应力张量分量；E表示比内能。
3 SPH算法验证
3.1 钢板台架试验
依照AEP-55 Vol试验标准[13]建立爆炸冲击钢板试验台架(见图2、图3)，台架由支撑横纵梁焊接而成，在2个横梁之间装配钢板安装平台，钢板安装平台由40 mm厚上下板和立板焊接组成，材料均为q235钢。在钢板安装平台均匀配重7.5 t模拟车重，配重后台架质量为9 t。高强钢通过螺栓与安装平台连接，高强钢的屈服强度为1 145 MPa，尺寸为1 500 mm×1 500 mm×20 mm，使用水平仪让钢板与地面保持平行，安装完成后，钢板下表面距地面400 mm。
本次试验为浅埋8 kg当量TNT下钢板底部正中心爆炸试验。试验中的TNT爆炸物为圆柱状，高径比为1∶3，爆炸物上表面距土壤表面100 mm，起爆点为圆柱体中心距上表面1/3处。试验后，高强钢板中心处出现凸起，四周向下翘曲，同时钢板与安装平台的螺栓连接全部断裂脱离，试验后钢板塑性变形如图4所示，使用3D扫描设备得到钢板最大塑性变形量为221.8 mm(见图5)。
图2 试验台架场景图
Fig.2 Test bench
图3 试验台架示意图
Fig.3 Schematic diagram of test bench
图4 试验后钢板变形情况
Fig.4 Deformation of steel plate after test
图5 钢板迎爆面3D扫描成像
Fig.5 3D scanning imaging of the bursting surface of the steel plate
3.2 钢板台架仿真
参照爆炸冲击钢板试验台架建立有限元模型，其中配重不考虑变形故采用材料模型MAT_RIGID的实体单元模拟，横纵梁、钢板安装台架以及钢板采用材料模型MAT_PLASTIC_KINEMATIC的壳单元模拟，钢板与安装平台之间的螺栓采用考虑失效的BEAM单元连接，设定其屈服强度为 1 500 MPa，当螺栓应力超过其屈服强度时认为螺栓失效。高强钢板与台架的材料参数见表1所示，台架的网格控制在10～20 mm，分别采用ALE算法和SPH算法模拟TNT爆炸物对钢板结构响应。
表1 试验台架中材料参数
Table 1 Material parameters in the test bench
材料ρ/(kg·m-3)E/GPaσs/MPaεfsQ2357 8502052350.25高强钢7 85021011450.12
在爆炸冲击钢板台架仿真中，ALE算法中通过关键字*CONSTRAINED_LAGRANGE_IN_SOLID完成拉格朗日网格与欧拉网格的耦合，其中炸药、土壤、空气为欧拉网格，试验台架及高强钢板为拉格朗日网格。并通过关键字*INITIAL_VOLUME_FRACTION_GEOMETRY定义TNT的形状和位置。
SPH算法的基本思想是将有限单元离散化，在爆炸冲击钢板台架模型中用SPH粒子模拟爆炸物与土壤，粒子之间距离设置为5 mm，粒子间光滑长度设置为1.2，添加控制卡片*CONTROL_SPH设置粒子作用范围。爆炸发生后，SPH粒子以极高的速度飞散，接触算法难以捕捉粒子，为了抑制粒子穿透钢板，模型中采用较小的粒子间距5 mm并调整接触刚度，增大接触惩罚系数。添加关键字*CONTACT_NODES_TO_SURFACE完成SPH粒子与拉格朗日单元的耦合。
仿真后得到在8 kg当量TNT下钢板的结构动态响应，图6为2种算法下钢板的残余变形云图，图7给出了1～3 ms两种算法的仿真动画，其中ALE算法中爆炸冲击波类似于流体在网格中流动；SPH算法中，冲击波以粒子形式运动传播。图8给出了2种算法下钢板中剖面变形曲线，可以看出2种算法下的钢板变形近似，SPH算法下的钢板最大塑性变形量为226 mm，与试验结果误差在2%之内。图9给出了2种算法下钢板能量曲线，2种算法下钢板能量变化趋势及最大值近似。图10为2种算法下钢板接触力曲线，SPH算法的钢板接触力出现时间略早于ALE算法。在计算时间上，SPH算法所用时间为ALE算法的75%，SPH算法计算效率更高。2种算法的仿真结果如表2所示，结果表明SPH法在模拟爆炸方面具有一定的准确性。
图6 钢板挠度云图
Fig.6 deflection cloud diagram of steel plate
图7 ALE算法(左)和SPH算法(右)仿真动画
Fig.7 Simulation comparison of ALE algorithm(left) and SPH algorithm(right)
图8 2种算法的钢板中剖面变形曲线
Fig.8 Deformation curve of center section of steel plate with two algorithms
图9 2种算法的钢板能量曲线
Fig.9 Steel plate energy curves of the two algorithms
图10 2种算法的钢板接触力曲线
Fig.10 Steel plate contact force curves of the two algorithms
表2 2种算法仿真结果对比
Table 2 Comparison of simulation results of the two algorithms
爆炸算法钢板残余变形量/mm误差/%钢板最大内能/kJ钢板最大动能/kJ网格单元数计算时间/sALE230.33.691 095.8855.2664 22618 146SPH226.01.891 048.1940.31 017 10113 623
4 榴弹爆炸破碎数值仿真
4.1 榴弹有限元模型
采用LS-DYNA建立浅埋155 mm榴弹模型，榴弹壳体外径D=155 mm，内径d=112 mm，榴弹长度L=640 mm，其中榴弹壳体质量为40 kg，填充TNT炸药，质量为7 kg。
榴弹壳体采用拉格朗日网格，网格大小为1～5 mm，材料为合金钢，采用带损伤失效的Johnson-Cook本构模型，材料参数取自文献[14]，榴弹炸药采用SPH粒子，粒子间距为10 mm，壳体与炸药的材料模型参数如表3～表5所示。建立的浅埋155 mm榴弹有限元模型如图11所示。
表3 榴弹壳体与炸药材料模型参数
Table 3 High explosive projectile shell and explosive material model
部件材料状态方程强度模型失效模型侵蚀准则壳体合金钢GRUNEISENJohnson-CookPrincipal Strain应用炸药TNTJWL无无应用
表4 榴弹壳体Johnson-Cook模型参数
Table 4 Johnson-Cook model of high explosive projectile shell
材料密度ρ/(g·cm-3)E/GPaG/GPa泊松比μA/MPaB/MPaCnm合金钢7.8321177.50.36186990.020 60.351.03材料D1D2D3D4D5合金钢0.10.302-0.785 30.016 80.826
表5 炸药材料模型以及JWL参数
Table 5 Explosive material model and JWL parameters
P/GPaD/(m·s-1)ρ/(g·cm-3)A/GPaB/GPaR1R2ωE0/(kJ·m-3)V216 9301.63373.773.734 74.1540.90.357 2001.0
图11 浅埋155 mm榴弹有限元模型示意图
Fig.11 Finite element model of shallow buried 155 mm high explosive projectile
设置榴弹起爆点为弹头，榴弹水平浅埋于土壤中，榴弹最高处距土壤表面100 mm。为了体现出榴弹破裂的随机性，在数值仿真中采用Stochastic随机破坏模型模拟壳体材料的随机破坏。在Stochastic模型中通过设定弱化网格来获得材料随机破化形式，获得不同大小、数量的破片。Stochastic随机破坏模型表达式为：
(5)
式中：P为应变为ε时的材料破坏概率；C和γ是与材料相关的参数，对于钢，通常取C=0.004 67，γ=10[15]。通过关键字*DEFINE_STOCHASTIC_VARITION定义材料的随机失效应变曲线。
整个浅埋155 mm榴弹模型共有单元670 650个，节点686 116个。炸药起爆后，炸药产生巨大的压力使榴弹壳体膨胀，局部破坏形成长条状自然破片，浅埋155 mm榴弹破片形成的过程如图12所示。
从起爆时刻到40 μs时，榴弹壳体从起爆处出现轴向破裂，随后在40 μs到160 μs之间，壳体裂纹从起爆处延伸至壳体尾部，至此榴弹自然破片呈现出长条状；在爆炸冲击波的作用下，长条状自然破片沿着径向破裂，t=400 μs时，破片与壳体尾部分离，长条状破片沿径向破裂成大小不一多条破片；t=600 μs时，破片继续沿径向膨胀，壳体最大直径达到220 mm，是原直径的1.4倍；在600 μs至1 000 μs破片继续沿径向膨胀，在t=1 000 μs时，破片分布呈现球状，最大直径达到550 mm，是原直径的3.5倍，最终破片形状以长条状为主，同时伴随着细小颗粒状破片。
图12 浅埋155 mm榴弹破片形成过程示意图
Fig.12 Formation process of shallowly buried 155 mm high explosive projectile fragment
沿着壳体轴向选取若干测试点，如图13所示。提取各个测试点的破片速度值如图14所示。从图14可以看出，壳体前部的破片速度较小，分布在50～100 m/s；壳体前中部的破片速度分布在100～250 m/s；壳体中部的破片占破片质量的大部分，破片速度分布在200～450 m/s，且破片末速度相对一致；壳体后部破片速度波动较大，破片速度分布在50～450 m/s。从破片速度大小来看，壳体中部的破片速度最大，且破片质量较大，因此壳体中部产生的破片威胁最大。图15给出了起爆后1.2 ms的破片速度分布，此时榴弹破片基本成形，破片速度最大值为427 m/s，集中在壳体的中部；壳体两端破片速度较低，破片速度分布呈现壳体中部破片速度高、两端破片速度低的特点。
图13 榴弹破片速度测试点示意图
Fig.13 High explosive projectile fragment velocity test point
图14 破片速度曲线
Fig.14 Fragment velocity curve
图15 t=1.2 ms时破片速度分布云图
Fig.15 Fragment velocity distribution at t=1.2 ms
4.2 浅埋榴弹对车辆结构响应
建立某车辆有限元模型，对其进行浅埋155 mm榴弹爆炸仿真，分析冲击波与破片耦合作用下车身的结构响应。车身底部防护组件材料为某防弹钢，主要参数见表6。
155 mm榴弹浅埋于土壤表面100 mm以下。榴弹起爆点位于弹头，榴弹轴向与车身纵向平行，且弹体位于车身正中心下，弹头朝向车身前部。图16给出了起爆后10 ms内榴弹爆炸冲击波和破片的传播图。在起爆后约2 ms，爆炸冲击波作用于车辆底部防护组件，并且爆炸冲击波要早于破片作用于防护组件(见图16)。随后在冲击波和破片的耦合作用下，在t=2.6 ms时防护组件面板应力云图如图17所示，其中面板最大应力达到2 592 MPa，远超面板材料的屈服强度，防护组件中心产生大变形，最大变形量达到345.3 mm，如图18所示。在整个爆炸响应过程中，爆炸冲击波对防护组件造成较大的变形，破片对面板造成较小的损伤，面板仅在部分螺栓连接处出现单元失效的情况(见图18)。
表6 防弹钢Johnson-Cook模型主要参数
Table 6 Parameters of Johnson-Cook model of bulletproof steel
AJC/MPaBJC/MPaCJCnD1D2D31 3421500.0630.2540.322.52-8.835
图16 冲击波与破片传播图
Fig.16 Shockwave and fragment propagation
图17 t=2.6 ms面板von Mises应力云图
Fig.17 Panel von Mises stress diagram t=2.6 ms
图18 防护组件最大变形云图
Fig.18 Maximum deformation diagram of protective assembly
5 结论
1) 本文采用ALE算法和SPH算法模拟浅埋地雷爆炸冲击钢板台架试验，对比2种算法的仿真动画和目标钢板的爆炸响应。SPH算法中可以清晰看到土壤粒子与目标结构的耦合效应，2种算法得到的钢板变形与试验结果近似，验证了SPH算法在处理浅埋爆炸问题的可行性。
2) 采用SPH粒子模拟爆炸物和土壤，分析了浅埋榴弹爆炸形成自然破片的过程。自然破片呈现两头速度低，中部速度大的特点，破片最大速度为427 m/s，且中部破片呈现长条状，两端破片呈现颗粒状，威胁较大的破片位于壳体中部。
3) 对某车辆进行底部浅埋榴弹爆炸仿真，结果表明爆炸冲击波要早于破片作用于防护组件。爆炸冲击波对防护组件作用效果表现为防护组件大变形，破片对防护组件作用效果表现为局部损伤。爆炸冲击波与破片联合作用下车底结构最大应力为2 592 MPa，远超材料屈服强度。
[1] Louis M French.Military traumatic brain injury:an examination of important differences[J].Annals of the New York Academy of Sciences,2010,1208(1).
[2] 吴梦阳,周云波,赵丰,等.浅埋炸药冲击试验台架数值模拟方法研究[J].兵器装备工程学报,2020,41(10):117-124.
Wu M Y,Zhou Y B,Zhao F,et al.Numerical simulation of shallow buried explosive impact test bench[J].Journal of Ordnance and Equipment Engineering,2020,41(10):117-124.
[3] 陈长海,侯海量,朱锡,等.破片式战斗部空中爆炸下冲击波与破片的耦合作用[J].高压物理学报,2018,32(01):148-156.
Chen C H,Hou H L,Zhu X,et al.Coupling effects of shock wave and fragments on air explosion of fragmented warhead[J].Chinese Journal of High Pressure Physics,2018,32(01):148-156.
[4] Zhu F,Zhao L M,Lu G X,et al.Structural response and energy absorption of sandwich panels with an aluminium foam core under blast loading[J].Advances in Structural Engineering,2009,11(5).
[5] 李利莎,谢清粮,郑全平,等.基于Lagrange、ALE和SPH算法的接触爆炸模拟计算[J].爆破,2011,28(01):18-22,27.
Li L S,Xie Q L,Zheng Q P,et al.Simulation of contact explosion based on lagrange,ALE and SPH algorithm[J].Blasting,2011,28(01):18-22,27.
[6] 魏然.爆炸冲击下车身结构防护机理及多学科优化研究[D].南京：南京理工大学,2017.
Wei R.Study on protection mechanism and multidisciplinary optimization of body structure under explosion impact[D].Nanjing:Nanjing University of Science and Technology,2017.
[7] 胡建宇,王显会,周云波,等.爆炸冲击下某车辆底部梁结构设计优化[J].兵器装备工程学报,2020,41(01):53-58.
Hu J Y,Wang X H,Zhou Y B,et al.Structural design optimization of vehicle bottom beam under explosion impact[J].Journal of Ordnance Engineering,2020,41(01):53-58.
[8] 石秉良,周孔亢,张云,等.基于SPH算法的驾驶室底部结构对爆炸冲击波响应数值仿真[J].机械工程学报,2016,52(16):132-139.
Shi B L,Zhou K K,Zhang Y,et al.Numerical simulation of response of cab bottom structure to blast shock wave based on SPH algorithm[J].Journal of Mechanical Engineering,2016,52(16):132-139.
[9] 孙新,白春华,王浩喆,等.预制破片战斗部安全距离试验与数值仿真研究[J].兵器装备工程学报,2021,42(06):27-33.
Sun X,Bai C H,Wang H Z,et al.Experimental and numerical simulation study on the safety distance of precast fragment warhead[J].Journal of Ordnance and Equipment Engineering,201,42(06):27-33.
[10] 孟德荣， 贾振铎， 赵铮.半穿甲弹侵爆复合防护结构仿真分析[J] .兵器装备工程学报，2020，41(11):46-50．
Meng D R， Jia Z D， Zha Z.Simulation analysis of composite protective structure of semi-piercing projectile penetrating and exploding[J] .Journal of Ordnance Equipment Engineering， 2020，41(11): 46-50．
[11] Wang W,Huang Y,Grujicic Mica,et al.Study of impact-induced mechanical effects in cell direct writing using smooth particle hydrodynamic method[J].Journal of Manufacturing Science and Engineering,2008,130(2).
[12] Monaghan J J,Gingold R A.Shock simulation by the particle method SPH[J].Academic Press,1983,52(2).
[13] NATO.Procedures for evaluating the protection level of logistic and light armored vehicles-mine threat[M].NATO Standardization Organization,2011:1-74.
[14] Babu V,Kankanalapalli S,Vunnam M.Effect of fragmenting buried high explosive projectile on military ground vehicles[C]//31st International Symposium on Ballistics.2019.
[15] 蒋建伟,张谋,门建兵,等.PELE弹侵彻过程壳体膨胀破裂的数值模拟[J].计算力学学报,2009,26(04):568-572.
Jiang J W,Zhang M,Meng J B,et al.Numerical simulation of shell expansion and fracture during PELE projectile penetration[J].Chinese Journal of Computational Mechanics,2009,26(04):568-572.