球形破片侵彻明胶运动模型及破片参数敏感性分析

蒋明飞,陈 莉,刘 坤,赵 磊,吴志林

(南京理工大学 机械工程学院, 江苏 南京 210094)

破片是轻武器弹药中常见的杀伤元,可按需预制成球形、菱形、小箭形等形状,并装入套筒式弹体中,通过炸药或发射药获得必要的初始速度以杀伤目标[1]。弹道明胶因与人体肌肉组织相似的密度、黏弹性及良好的均匀性、透明性,在生物医学和创伤弹道中被广泛用作软组织模拟物[2-4]。

众多学者开展了破片与明胶靶标的实验与仿真研究[5-15]。刘飞等[5]在实验研究的基础上,结合瞬态动力分析软件,对球形杀伤元在明胶介质中的运动规律进行了分析研究,并提出了明胶靶标的标定方法。温垚珂等[6]分别使用Lagrange法、任意拉格朗日-欧拉(arbitrary Lagrange-Euler, ALE)法和SPH-Lagrange耦合算法对钢球高速侵彻明胶的动力学过程进行了模拟,发现Lagrange法的计算精度和效率最高;并在此基础上使用Lagrange法模拟了三种球形破片中等速度侵彻明胶的作用过程[7]。文献[8-9]利用LS-DYNA软件模拟了球形破片与明胶的相互作用过程,研究了速度衰减、能量变化和空腔演化规律。李金明等[10]对4.8 mm球形破片侵彻明胶靶标的过程进行了数值仿真研究,发现破片质量是影响明胶靶标瞬时空腔最大直径及侵彻深度的主要因素。苑大威等[11]采用仿真方法研究了菱形破片以不同速度、姿态角侵彻明胶的作用过程,总结菱形破片翻转规律,得出菱形破片侵彻明胶深度的经验公式。Swain等[12]、Ye等[13]开展了不同直径的钢球以不同速度侵彻明胶的实验,建立了侵彻深度的经验公式。张志倩等[14]通过实验和数值模拟方法研究了低速陶瓷球侵彻明胶的过程,发现陶瓷球直径是影响侵彻深度的主要因素。田浩成等[15]利用遗传算法对不同工况的仿真数据进行拟合,得到了计算瞬时空腔的简化数学模型。

由于数值模拟方法计算效率较低,很多学者通过理论建模的方法分析了球形破片侵彻明胶的运动规律。Seglets[16]认为球形破片侵彻明胶时会受到惯性阻力和材料阻力的作用,其中材料阻力是应变率的指数函数,并建立了球形破片侵彻明胶的运动模型。Liu等[17]考虑球形破片分离角和速度的关系,建立了球形破片惯性阻力系数和最大瞬时空腔的联系;后又基于准静态柱形空腔膨胀理论,提出了包含惯性项和率相关强度项的球形破片侵彻明胶的阻力模型[18]。莫根林等[19-21]建立了长方体破片侵彻明胶的六自由度运动模型,并根据实验数据试算得到长方体破片运动模型中的动态阻力系数;利用球形破片在明胶中的运动模型和空腔实验数据,建立了空腔振幅、角频率与能量的经验公式。刘坤等[22-24]在考虑了惯性阻力和黏性阻力的模型基础上引入明胶静态强度项;后又提出了考虑明胶应变率效应的理论模型,使用最小二乘法确定了最佳惯性阻力系数和黏性阻力系数。以上这些理论模型均建立在破片已经完全进入明胶中的假设,忽略了球形破片由空气进入明胶这个跨介质的阶段。

本文基于动态空腔膨胀理论,考虑了球形破片由空气进入明胶这个阶段的速度衰减,建立了球形破片侵彻明胶的分段运动理论模型,通过实验确定了最优阻力系数。利用该模型分析了理论计算过程中的误差来源;并通过理论推导和Sobol′法去分析球形破片参数对侵彻深度影响的敏感性。

1 理论模型

根据动态空腔膨胀理论,假设侵彻靶材的弹头表面微元的受力由材料阻力和惯性阻力组成,其中材料阻力与靶材的屈服强度相关,惯性阻力与靶材密度和法向膨胀速度平方的乘积相关[25]。借鉴上述分析方法,假设球形破片为刚性球体,侵彻明胶时迎风面上各微元面仅受垂直于微元面的惯性阻力和材料阻力的作用,则微元所受阻力f可表示为:

f=[AY+Bρt·(v·n)2]·n·dS

(1)

式中,Y和ρt分别是明胶材料的屈服强度和密度,A和B分别是无量纲材料阻力系数、惯性阻力系数,v为微元的速度,dS为微元的面积,n为微元的外法线方向。

为描述球形破片在明胶中的运动规律,建立固定坐标系O′x′y′z′,使得O′y′z′平面和明胶入射平面重合。在球形破片的质心建立局部坐标系Oxyz,使得x轴和x′轴平行、y轴和y′轴平行。坐标系与明胶的相对位置如图1所示。

图1 模型坐标系Fig.1 Coordinate system of the model

仅考虑球形破片在x′方向的平动,不考虑球形破片的空间转动和在其他方向的平动,则球形破片表面微元的速度和质心的速度相同均为v。在局部坐标系Oxyz中,令微元面的直角坐标为(x,y,z)、柱坐标为(r,θ,x),两者的关系可表示为:

(2)

由于微元处于球形破片的表面,其位置坐标满足球面几何条件:

x2+y2+z2=R2

(3)

式中,R为球形破片半径。将式(2)代入式(3)可得:

x2+r2=R2

(4)

对式(3)等号两侧进行微分运算,可得微元的单位外法线n为:

(5)

将式(2)和式(4)代入式(5),得到单位外法线n的柱坐标形式为:

(6)

柱坐标下微元的面积大小dS可表示为:

dS=R·dθ·dx

(7)

若微元处在明胶内,在固定坐标系中,需满足以下关系:

x′=xc+x>0

(8)

式中,x′、xc分别为微元和球形破片质心在固定坐标系O′x′y′z′中的坐标。整理式(8)可得:

x>-xc

(9)

假设球形破片侵彻明胶过程中,仅破片的前半球面与明胶接触作用,则球形破片侵彻明胶的过程可分为两个阶段:未完全侵入阶段和完全侵入阶段。在未完全侵入阶段,-R

将式(6)～(8)代入式(1)并积分,得到球形破片在两个阶段的运动阻力F为:

(10)

式中,vn为微元的法向速度,vn=v·n=v·x/R。

根据牛顿第二定律,球形破片的运动方程可表示为:

(11)

其中,m是球形破片的质量。

第一阶段:球形破片所受的阻力与位移相互关联耦合,难以求出位移的解析解。令ve为球形破片的前半球面刚好完全进入明胶时的速度,第一阶段结束时,球形破片速度由入靶速度v0衰减为ve,侵彻明胶深度Xc1=R。

第二阶段:球形破片的前半球面完全进入明胶后,其阻力形式相对简单,联立式(10)和式(11),有如下关系:

(12)

式中,m=ρp4πR3/3,R=D/2,ρp为球形破片密度,D为球形破片直径。将m代入式(12)可得:

(13)

将式(13)移项,并对两边同时积分可得:

(14)

对式(14)进行化简,可得球形破片在第二阶段位移xc2的表示式为:

(15)

特别地,当球形破片速度v衰减为0时,即第二阶段结束时,球形破片在明胶中的侵彻深度Xc2为:

(16)

球形破片在明胶中总的侵彻深度Xc=Xc1+Xc2,引入无量纲侵彻深度Xc/D:

(17)

式(17)给出了刚性球侵彻软介质类靶标的一般表达式,具有普适性。若只考虑球形破片在第二阶段的动能损失,定义符号ψ为球形破片在单位侵彻距离上的动能损失,其表达式为:

(18)

定义符号χ为第二阶段球形破片在单位侵彻距离上的比动能损失,其表达式为:

(19)

式中,ΔEsk2代表球形破片在第二阶段的比动能损失,ΔEsk2=ΔEk2/Smax,其中Smax代表球形破片的最大横截面积,Smax=πD2/4。

2 实验方法与结果

将照相明胶颗粒与自来水按1 ∶9重量配比,置于恒温60 ℃的水浴炉中保温1 h,然后将澄清透明的明胶溶液注入尺寸为350 mm×250 mm×200 mm的不锈钢模具中,待冷却至室温后放置在4 ℃的恒温箱中保温24 h,最后脱模后再保温2 h,即可进行球形破片侵彻明胶的实验。

为了便于发射球形破片,需要将其放置在塑料弹托中,如图2所示,并用黄油填满弹托的缝隙,然后将弹托装入带发射药的7.62 mm枪弹药筒中,使用弹道枪瞄准明胶中心位置进行垂直射击。射击时,通过光电靶测量球形破片进入明胶前的速度,光电靶的中心与明胶的距离为1 m,光电靶之间的距离为1.2 m。高速摄像机用于拍摄球形破片在明胶中的运动过程,拍摄帧率为10 000帧/s,分辨率为640×320 px。实验系统示意图如图3所示,其中明胶侧面的光源用于提高拍摄画面的质量。

图2 用于发射球形破片的弹托Fig.2 Sabots used to fire spherical fragments

图3 实验系统示意图Fig.3 Schematic diagram of the experimental system

实验中各发射了三种不同直径的钢质球形破片,破片编号为#1～#3,其直径D、质量m和通过光电靶的时间差Δt如表1所示。三种球形破片侵彻明胶的运动过程如图4所示,可见球形破片在明胶中的运动可视为水平直线运动,其在竖直方向的运动可忽略不计。使用图像后处理软件的测量功能,可获得球形破片的位移随时间的变化规律。

表1 球形破片的参数

(a) #1球形破片的侵彻过程(a) Penetration process of the #1 spherical fragment

(b) #2球形破片的侵彻过程(b) Penetration process of the #2 spherical fragment

(c) #3球形破片的侵彻过程(c) Penetration process of the #3 spherical fragment

3 结果讨论

3.1 阻力系数确定

根据表1中球形破片通过光电靶的时间和光电靶之间的距离,求得#1～#3球形破片入射明胶的入靶速度分别为651 m/s、712 m/s和931 m/s。

令实验中#i球形破片进入明胶中第1帧对应的时刻为Ti1,则第j帧对应的时刻Tij为:

Tij=Ti1+dt×(j-1)

(20)

式中,dt为高速摄像机的拍摄间隔。

令Tij时刻球形破片位移的实验值为pij,理论计算值为yij,则位移的实验值和理论值的均方根误差总和σ为:

(21)

式中,ni为#i球形破片的有效实验数据个数。

为使运动模型能够较好地模拟所有球形破片的运动规律,阻力系数A和B的取值要使σ最小。已知明胶材料的密度ρt=1 030 kg/m3,屈服强度Y=2.2×105Pa[7],使用式(10)～(11)和位移的实验数据,借助四阶龙格-库塔法和直接结果搜索法,编程求解阻力系数。经计算,最优的阻力系数为A=2.77,B=0.34。

球形破片位移的实验值和理论值的比较如图5所示,可以看出实验值与理论模型计算值误差很小,说明分段运动模型能较好地模拟不同直径的钢质球形破片侵彻明胶的运动规律。

图5 位移的实验值和理论值的比较Fig.5 Comparison between theoretical and experimental displacements

3.2 模型验证及误差分析

以同样的实验方法测试直径5.16 mm、质量1.39 g的钨球侵彻明胶的运动位移,光电靶测得钨球入射明胶的入靶速度为800 m/s。钨球破片位移的实验值和理论值的比较如图6所示,说明分段运动模型能够有效预测不同密度的球形破片侵彻明胶的运动规律,具有普适性。

图6 钨球的运动位移Fig.6 Motion displacements of the tungsten ball

本文使用的#1～#3球形破片入射明胶的入靶速度v0是通过光电靶测得的,速度ve是通过式(10)～(11)计算得到的,表2比较了两种速度之间的差距。速度v0和ve之间的差距是由于球形破片在未完全侵入阶段的速度衰减,当前实验条件下球形破片在未完全侵入阶段的速度衰减为6～10 m/s。多数学者认为的球形破片在明胶中运动的初始速度实际是ve,而他们大都采用光电靶测得的入靶速度v0来替代,这样就在初始速度的取值上产生了误差,导致位移的计算结果往往不够准确。

表2 球形破片在不同时刻的速度

球形破片在未完全侵入阶段的速度衰减受破片直径D、入靶速度v0和质量m的影响。图7为#1和#2球形破片以不同入靶速度冲击时,在未完全侵入阶段的速度衰减。可见随着入靶速度的增加,不同直径的钢质球形破片的速度衰减均呈线性增加,约占入靶速度的1%。

图7 入靶速度对未完全侵入阶段速度衰减的影响Fig.7 Effect of entering-target velocity on velocity attenuation in the incomplete entering stage

不同质量的#1和#2球形破片,以入靶速度650 m/s冲击时,在未完全侵入阶段的速度衰减如图8所示。当直径相同时,质量越小,即密度越低,速度衰减越大,当密度特别低时,球形破片在未完全侵入阶段的速度衰减可达几十到上百米每秒;当质量相同时,大直径低密度的球形破片速度衰减大于小直径高密度的球形破片,但随着密度的增加,两者之间的差距在逐渐减小。结果表明:密度对球形破片在未完全侵入阶段的速度衰减影响较大,低密度的球形破片在未完全侵入阶段的速度衰减将不能忽略。

图8 质量对未完全侵入阶段速度衰减的影响Fig.8 Effect of mass on velocity attenuation in the incomplete entering stage

球形破片在未完全侵入阶段存在速度衰减,考虑和忽略此阶段的运动模型在计算侵彻深度上存在误差。未完全侵入阶段对侵彻深度的影响规律如图9所示。相同密度下,低入靶速度破片的侵彻深度误差大于高入靶速度破片;相同入靶速度下,密度较低的破片侵彻深度误差较大,随着破片密度的增加,侵彻深度误差逐渐减小。入靶速度为600 m/s时,密度为7 800 kg/m3(钢)的破片侵彻深度误差约为0.6%,而密度为1 000 kg/m3(高聚物)的破片侵彻深度误差约为4.5%。结果表明:低密度、低入靶速度的球形破片忽略未完全侵入阶段对侵彻深度的影响较大。

图9 未完全侵入阶段对侵彻深度的影响Fig.9 Effect of the incomplete entering stage on penetration depth

3.3 破片参数的敏感性分析

由上文的分析可知,球形破片的参数直径D、密度ρp和入靶速度v0均会对侵彻明胶的过程产生影响,首先用理论推导分析球形破片参数对侵彻深度等参数影响的敏感性。将已知的相关参数代入式(17)～(19),发现Xc/D是与球形破片密度ρp、初始速度ve相关的函数;ψ是与球形破片直径D、初始速度ve相关的函数;而χ是只与初始速度ve相关的函数。通过进一步的量纲分析可知,ψ反映了球形破片侵彻明胶过程中的平均阻力;χ则反映了球形破片表面微元侵彻明胶过程中的平均阻抗应力。

为了更直观地判断球形破片参数密度ρp、直径D和速度ve对参数Xc/D、ψ影响的敏感性,根据适用于侵彻明胶的轻武器弹药的相关性能指标[26],取密度ρp变化范围为1 000～20 000 kg/m3,直径D的变化范围为3～9 mm,速度ve的变化范围为0～1 000 m/s,作三维空间曲面图来分析。图10给出了Xc/D在(ρp,ve)平面上的变化,显然,无量纲侵彻深度对速度ve更敏感而非对密度ρp更敏感。图11给出了ψ在(D,ve)平面上的变化,可以看出球形破片侵彻明胶过程中的平均阻力也对速度ve更敏感。此外,图10～11中三维曲线顶部的等高线可用于指导球形破片的参数优化设计。

图10 无量纲侵彻深度Xc/D在(ρp,ve)平面上的变化Fig.10 Variation of dimensionless penetration depth Xc/D on (ρp,ve) plane

图11 参数ψ在(D,ve)平面上的变化Fig.11 Variation of parameter ψ on (D,ve) plane

在理论公式特别复杂的情况下,有时无法求出位移的解析解,还可以借助Sobol′法来分析各因素的敏感性。Sobol′法是一种广泛应用于全局敏感性分析的算法[27],可用基于方差的蒙特卡罗法,通过采样计算模型响应的总方差及各项偏方差,从而求得敏感性[28]。Sobol′法的主要思想是:令输入变量的定义域为单位空间,即Ωk=(x|0≤xi≤1;i=1,…,k),将函数f(x1,…,xk)分解成2n项多级函数的总和,即:

f1,2,…,k(x1,…,xk)]

(22)

然后求解f(x1,…,xk)的总方差:

(23)

式(22)中其他项的偏方差表示为:

(24)

敏感度系数定义为式(24)与式(23)的比值:

(25)

对式(23)两边同时除以D,可以推导出:

(26)

其中,Si称为因素的一阶敏感度系数,它代表因素xi对输出的主要影响,即对方差的贡献大小,根据此定义有:

(27)

类似地,Sij(i≠j)被称为二阶敏感度系数,以评估xi和xj两因素耦合作用对总体方差的影响。 同理,评估给定参数的主要影响和有关该变量的所有相互作用的总敏感度系数可以表示为:

(28)

式中,～i表示不包含i。

利用式(27)～(28)配合四阶龙格-库塔法,再结合Sobol′法编程,求解输入参数直径D、密度ρp和速度ve对输出参数Xc、Xc/D、ψ和χ影响的敏感度系数。图12给出了各输入参数对输出参数影响的敏感度系数,系数值越大表明输出参数对该因素越敏感,系数值特别小代表此因素对输出参数影响微乎其微,即此因素与输出参数不相关。由图12可知,用Sobol′法计算的结果与上文理论分析得出的结论基本一致。敏感度系数的影响规律为:影响侵彻深度Xc的敏感度系数由大到小依次是速度、密度和直径;影响无量纲侵彻深度Xc/D的敏感度系数由大到小依次是速度和密度,与直径无关;影响平均阻力ψ的敏感度系数由大到小依次是速度和直径,与密度无关;而平均阻抗应力χ只与速度相关。对比图12(a)和图12(b)可知,以输出参数无量纲侵彻深度Xc/D为例,速度和密度的总敏感度系数均大于一阶敏感度系数,表明这两个因素关联耦合影响的程度大于单一因素。因此,引入无量纲碰撞函数I和无量纲密度比λ,其表达式为:

(29)

将式(29)代入式(17)可得:

(30)

显然,这两个无量纲参数I和λ,控制并决定着刚性球侵彻明胶的深度,式(30)的形式同样适用于刚性球侵彻其他软介质类的情况。

(a) 一阶敏感度系数(a) First-order sensitivity coefficients

(b) 总敏感度系数(b) Total sensitivity coefficients

4 结论

本文为揭示球形破片对人体组织致伤机理,通过理论分析的方法,开展球形破片侵彻明胶的研究。基于动态空腔膨胀理论,考虑球形破片未完全侵入阶段的速度衰减,建立了球形破片侵彻明胶的分段运动理论模型,并通过实验进行验证。分析了理论计算过程中的误差来源,并推导得到了无量纲侵彻深度的表达式。利用Sobol′法进行了球形破片参数(直径、密度和速度)对侵彻深度等参数影响的敏感性分析。得出以下结论:

1)理论模型和实验结果符合较好,表明本文所建立的运动模型能较好地描述球形破片侵彻明胶的运动规律;

2)未完全侵入阶段的速度衰减与破片的密度密切相关,钢质球形破片的速度衰减约为入靶速度的1%,破片密度较低时速度衰减较大,不可忽略;

3)理论方法和Sobol′法对破片参数的敏感性分析规律基本一致,破片参数对侵彻深度影响的敏感性由高到低依次是速度、密度和直径。

本文推导的侵彻深度公式,可为杀伤元参数优化设计和毁伤评估提供理论参考。

致谢

63961部队的陈川琳工程师,在Sobol′法的算法上提供了帮助和指导,谨致谢意!