基于BP神经网络预测炸药JWL状态方程参数对EFP速度的影响

郝礼楷,顾文彬,谢兴博,张亚栋,钟明寿,邹绍昕,康耕新

(1.陆军工程大学 野战工程学院, 南京 210000; 2.中国人民解放军31539部队, 北京 100000;3.东南大学 爆炸安全防护教育部工程研究中心, 南京 210000;4.陆军工程大学 国防工程学院, 南京 210000)

0 引言

在高温高压的爆轰产物作用下,金属药型罩翻转闭合形成具有一定形状、较高速度的弹丸,该弹丸被称为爆炸成型弹丸(explosively formed projectile,EFP)。EFP的速度越大,动能也越大,侵彻能力也越强[1],这是评价毁伤威力的重要指标,对战斗部设计与威力评估具有十分重要意义。因此,学者们对EFP速度进行了大量研究,提出了多种工程测量方法[2]、理论计算模型[3]和数值仿真模型[4]。炸药作为聚能装药的关键部件,对EFP速度的影响至关重要。然而,以往工作[5-7]的关注重点集中在炸药密度、爆速和爆压等因素,而对炸药状态方程参数的影响规律研究较少。

爆轰产物作用及状态变化实际上是由炸药状态方程所控制的,对参数取值较为敏感。Jones-Wilkins-Lee(JWL)状态方程形式简单,能够准确反映爆轰产物的膨胀驱动做功能力,广泛运用于武器设计、爆轰工程。JWL状态方程中包含多个未知参数,需要结合圆筒试验,经过大量修正计算并基于经验进行人为“调参数”,取值过程复杂,人为因素影响较大。另外有学者通过贝叶斯分析方法[8]、遗传优化算法[9]、非嵌入式多项式混沌法[10]等方法确定JWL状态方程参数,但是并未得到广泛运用。总体而言,目前尚未形成操作性强的参数确定方法,无法给出精确数值,拟合结果因人而异。因此,需要对JWL状态方程参数的快速鉴别校正和合理取值提出新的解决方法,满足EFP速度的研究需求。机器学习具有强大的数据处理与分析能力,能够避免繁重的计算任务,具有较高准确度,引入BP神经网络可以有效解决上述问题。

为此,本文中研究JWL状态方程各特征参数对EFP速度的影响,建立考虑JWL状态方程参数的EFP速度计算公式。在此基础上,以JWL状态方程各特征参数和EFP速度分别作为BP神经网络的输入和输出进行反复训练与测试,最终得到精度较高的EFP速度预测网络模型。

1 JWL状态方程

JWL状态方程的标准压力形式为[11]:

(1)

JWL状态方程由3项构成,方程式右端第1项主要控制高压段作用,第2项主要控制中压段作用,第3项在低压段作用明显,即状态方程参数控制爆轰产物的高、中、低压3个阶段。P-V曲线如图1所示[12]。

图1 P-V曲线Fig.1 P-V curve

等熵条件下,JWL状态方程为[9]:

(2)

理想爆轰时,根据C-J条件和Hugoniot关系式,爆轰产物满足以下3个约束守恒方程[9],即:

(3)

2 EFP成型数值仿真

2.1 装药结构

选用顾文彬等[13]优化设计得到的Φ65 mm EFP聚能装药,装药结构如图2所示,主装药为压装JH-2炸药,装药长径比为1,采用船尾形装药结构。药型罩由紫铜制成,采用中间厚、边缘薄的变壁厚球缺型结构。药型罩的顶部厚度δ=2.1 mm,内曲率半径r1=67 mm,外曲率半径r2=62 mm。

图2 Φ 65 mm EFP结构图Fig.2 Structure diagram of Φ 65 mm EFP

2.2 数值模型

数值模型由药型罩、炸药和空气构成,考虑到模型对称性和减小计算量,使用LS-DYNA软件分别建立如图3所示的1/2二维有限元模型和1/4三维有限元模型,单元均采用多物质ALE算法。起爆点位主装药上端面中心处,在模型的对称面上设置对称约束,对空气施加透射边界条件,计算采用cm-g-μs单位制。紫铜药型罩采用MAT_JOHNSON-COOK模型和EOS_Gruneisen状态方程描述;JH-2炸药采用MAT_HIGH_EXPLOSIVE_BURN模型和EOS_JWL状态方程描述;空气采用MAT_NULL模型和EOS_LINEAR_POLYNOMIAL状态方程描述。材料参数取值如表1—表3所示[14-16]。

表1 空气材料参数Table 1 Material parameters of air

表2 JH-2炸药材料参数Table 2 Material parameters of JH-2 explosive

表3 紫铜材料参数Table 3 Material parameters of red copper

图3 EFP成型的数值模型Fig.3 Numerical simulation model of EFP forming

2.3 仿真结果

将数值仿真得到的EFP形态与试验X光照片[13]进行对比,如图4所示,EFP长径比和轮廓形态吻合度较高,马鞍形收缩在仿真中得到了较好体现,表明数值模型能够有效模拟EFP的速度梯度。EFP的速度V、长度L、前部最大直径D1、中部最大直径D2和长径比L/D2的数值仿真结果与X光试验结果的定量对比如表4所示。

表4 EFP成型仿真与X光试验的结果对比Table 4 Quantitative comparison between numerical simulation and X-ray test

图4 仿真结果与X光试验对比Fig.4 Comparison between numerical simulation and X-ray test

仿真结果与试验结果较为一致,最大相对误差为6.69%,表明上述数值模型构建和材料参数选取较为合理,仿真结果可信。考虑到三维仿真较为耗时,后续仅进行二维仿真计算。

3 JWL状态方程参数对EFP速度的影响

根据表2中JH-2炸药的材料参数确定JWL状态方程参数的初始值,即A=6.1,B=0.068 01,R1=4.1,R2=1.3,ω=0.36,E0=0.1。采用控制变量法,改变某一特征参数的取值,其余5个参数的取值保持不变(即初始取值)。在此基础上,得到JWL状态方程不同取值条件下EFP速度的仿真结果,分析影响规律,拟合变化曲线。

将文献[5]中5种常见炸药的材料参数作为JWL状态方程参数的取值范围,具体取值范围如下,A:3.5～8.5,B:0.02～0.22,R1:4～6.5,R2:0.8～1.9,ω:0.2～0.375,E0:0.07～0.1。各特征参数的增量分别为A:0.5,B:0.02,R1:0.25,R2:0.1,ω:0.025,E0:0.002 5。

3.1 参数A的影响

当参数A的取值范围取3.5～8.5,变化增量取0.5时,EFP速度随参数A取值变化的变化规律如图5所示。

图5 EFP速度随参数A的变化曲线Fig.5 The change curve of EFP velocity with parameter A

随着参数A的增加,EFP速度呈指数增长趋势,增长趋势较为稳定,拟合得到基于参数A的EFP速度计算表达式为:

VA=4 470.76-2 522.08e(-A/20.41)

(4)

3.2 参数B的影响

当参数B的取值在0.02～0.22,变化增量取0.02时,EFP速度随参数B取值变化的变化规律如图6所示。

图6 EFP速度随参数B的变化曲线Fig.6 The change curve of EFP velocity with parameter B

随着参数B的增加,EFP速度呈指数增长趋势,增长趋势较为稳定,拟合得到基于参数B的EFP速度计算表达式为:

VB=9 118.25-6 683.33e(-B/2.408)

(5)

3.3 参数R1的影响

当参数R1的取值范围取4～6.5,变化增量取0.25时,EFP速度随参数R1取值变化的变化规律如图7所示。

图7 EFP速度随参数R1的变化曲线Fig.7 The change curve of EFP velocity with parameter R1

随着参数R1的增加,EFP速度呈指数减小趋势,减小趋势逐渐降低且趋于平稳,拟合得到基于参数R1的EFP速度计算表达式为:

VR1=1 932.22+36 516.08e(-R1/1.029)

(6)

3.4 参数R2的影响

当参数R2的取值范围取0.8～1.9,变化增量取0.1时,EFP速度随参数R2取值变化的变化规律如图8所示。

图8 EFP速度随参数R2的变化曲线Fig.8 The change curve of EFP velocity with parameter R2

随着参数R2的增加,EFP速度呈指数减小趋势,减小趋势逐渐降低且趋于平稳,拟合得到基于参数R2的EFP速度计算表达式为:

VR2=2 480.32+1 408.52e(-R2/0.537)

(7)

3.5 参数ω的影响

当参数ω的取值范围取0.2～0.375,变化增量取0.025时,EFP速度随参数ω取值变化的变化规律如图9所示。

图9 EFP速度随参数ω的变化曲线Fig.9 The change curve of EFP velocity with parameter ω

随着参数ω的增加,EFP速度呈指数增长趋势,增长趋势较为稳定,拟合得到基于参数ω的EFP速度计算表达式为:

Vω=2 708.62-611.70e(-ω/0.214)

(8)

3.6 参数E0的影响

当参数E0的取值在0.07～0.1,变化增量取0.002 5时,EFP速度随参数E0取值变化的变化规律如图10所示。

图10 EFP速度随参数E0的变化曲线Fig.10 The change curve of EFP velocity with parameter E0

随着参数E0的增加,EFP速度呈指数增长趋势,增长趋势较为稳定,拟合得到基于参数E0的EFP速度计算表达式为:

VE0=-6 595.02+8 445.85e(E0/1.186)

(9)

4 基于拟合公式计算EFP速度

传统EFP速度计算公式[3,17]主要考虑炸药密度、爆速、爆压、药型罩材料及结构等因素,未能引入炸药状态方程参数的影响,计算结果存在偏差。如表5所示,对于JH-2炸药,在装药密度、爆速和爆压非常接近甚至相同的情况下,由于JWL状态方程参数取值受到人为因素的干扰,各特征参数均存在不同程度的偏差。仿真结果的精度很大程度上取决于JWL状态方程参数的取值,因此很有必要引入JWL状态方程的影响。

表5 不同文献中的JH-2炸药参数Table 5 Parameters of JH-2 explosive in different literatures

4.1 EFP速度计算表达式

JWL状态方程的A、B、R1、R2、ω、E0参数取值由最小值到最大值时,EFP速度的变化率分别为:ηA=(2 570-2 340)/2 570=0.165 7,同理ηB=0.164 8;ηR1=0.247 2;ηR2=0.096 7;ηω=0.046 9;ηE0=0.088 8。

总变化率为:η=ηA+ηB+ηR1+ηR2+ηω+ηE0=0.812 8。

各参数影响权重分别为:σA=ηA/η=0.203 9,同理σB=0.202 8;σR1=0.304 1;σR2=0.119 0;σω=0.061 0;ηE0=0.109 2。

基于数值仿真结果,分析JWL状态方程参数对EFP速度的影响权重。R1的影响最大,ω的影响最小,各特征参数对EFP速度影响的大小顺序是:R1、A、B、R2、E0、ω。从理论角度分析原因,R1和A控制爆轰产物在高压段的膨胀驱动做功,显著影响药型罩翻转变形。使用控制变量法分析各参数对EFP速度的影响规律,是基于R1取4.1进行的数值仿真,因此可以引入修正系数K(K=4.1/R1),消除R1取值偏差给计算结果带来的影响。由于JWL状态方程各参数对EFP速度影响的拟合公式均为指数形式,因此可以将各个拟合公式进行线性叠加,得到式(10)所示的基于JWL状态方程参数表示的EFP速度计算公式。

VEFP=σR1VR1+K(σAVA+σBVB+σR2VR2+

σωVω+σE0VE0)=587.59+11 104.54e(-R1/1.029)+

K(2 500.98-514.25e(-A/20.41)-

1 355.38e(-B/2.408)+167.61e(-R2/0.537)-

37.31e(-ω/0.214)+922.29e(E0/1.186))

(10)

4.2 公式计算与仿真结果的对比验证

保持JH-2炸药的密度为1.7 g/cm3、爆速为8 400 m/s和爆压为30 GPa不变,采用表5中不同文献给出的JH-2炸药的JWL状态方程参数,开展数值仿真,仿真结果与公式计算结果对比如表6所示。仿真结果与公式计算对比如图11所示。

表6 EFP速度的仿真结果与公式计算对比Table 6 Comparison between numerical simulation and formula calculation of EFP velocity

图11 仿真结果与公式计算对比Fig.11 Comparison between numerical simulation and formula calculation

由图11可知,仿真结果与公式计算结果较为接近,对于取值方案1、3、5、6、7、9,两者吻合度较高,最大误差为6.87%;对于取值方案2、4、8,两者偏差较大,最大误差为13.21%。由此得出结论,使用K作为修正系数得到的速度表达式计算EFP速度是可行的。部分计算结果与仿真结果仍存在一定偏差,需要针对此方法加以改进,因此下文引入基于BP神经网络的EFP速度预测方法。

5 基于BP神经网络预测EFP速度

5.1 BP神经网络构建

BP神经网络是一种按误差逆传播训练的多层前馈神经网络,可以实现信号正向反馈以及误差反向传播。正向传播中,信号从输入层经隐含层神经元进行逐层处理,最终传输至输出层[24]。每一层神经元状态只影响下一层神经元状态,如果输出层未取得期望输出,则自动转入反向传播,使用梯度下降法调整各层网络连接的权值和阈值,边传播边修正,交替进行并反复训练,直至结果误差趋向设定的极小值,最终实现预测输出值逼近期望输出值。

设定该网络的输入为JWL状态方程A、B、R1、R2、ω、E0参数,输出为EFP速度。输入参数为6个,输出参数为1个,可以采用单隐含层的前馈神经网络结构。隐含层的节点数将显著影响网络结构的可靠性,根据试算法[25]确定隐含层神经元数目取4,网络结构为6-4-1,如图12所示。

图12 BP神经网络结构Fig.12 The structure of BP neural network

网络学习的样本数据来源于前文拟合的6个函数表达式,每个拟合函数生成10 000组训练样本,共计60 000组样本作为训练集,表6中9组结果数据作为测试集。在确定神经网络模型结构以及数据集来源之后,为提高网络学习效率,消除数据之间的数量级差别,采用mapminmax函数对训练样本进行归一化处理,所有数据归一化到区间[-1,1]内,函数形式为:

Y=(Ymax-Ymin)*(X-Xmin)/(Xmax-Xmin)+Ymin

(11)

式(11)中:Y为归一化后的值;X为输入值;Xmin和Xmax分别为输入值中的最小值和最大值;Ymin和Ymax为自定义常数,分别设定为-1和1。

使用newff函数构建BP神经网络,为保证预测结果的可靠性,需要选择合适的传递函数和学习函数。其中隐含层的传递函数为tansig,输出层的传递函数为purelin,学习算法为traingdx,设定网络目标误差为0.000 1,同时反复调整与训练有关的其他参数。当误差低于设定的目标误差值或达到最大训练次数,训练终止。网络达到相对稳定的状态之后对输出结果进行反归一化处理,得到输出结果的实际数据,经过多次训练后得到了精度较高的速度预测网络模型。

5.2 网络预测与公式计算、数值仿真的对比验证

测试集的9组预测结果如表7和图13所示,网络预测结果与仿真结果非常接近,误差不超过2.3%,远小于公式计算的误差。由此表明,网络训练效果良好,具有较高精度。

表7 网络预测与公式计算、数值仿真对比Table 7 Comparison of network prediction,formula calculation and numerical simulation

为了验证上述建立的BP神经网络的适用性和准确性,不再按照表5中不同文献给出的JH-2炸药参数进行取值,而是保持JH-2炸药的密度为1.7 g/cm3、爆速为8 400 m/s和爆压为30 GPa不变,对JWL状态方程参数进行随机取值,取值方案如表8所示。继续开展数值仿真,并与网络预测值、公式计算结果进行对比,对比结果如表9和图14所示。

表8 JWL状态方程的随机取值方案Table 8 The random value scheme for the JWL equation of state

表9 随机取值方案下EFP速度的网络预测与公式计算、数值仿真对比Table 9 Comparison of network prediction,formula calculation and numerical simulation of EFP velocity under random value scheme

图14 随机取值时网络预测与公式计算、数值仿真对比Fig.14 Comparison of network prediction,formula calculation and numerical simulation under random value scheme

由图14可知,网络预测、公式计算与数值仿真结果均吻合较好。但是,公式计算的部分结果误差依旧较大,最大误差为23.4%,网络预测的结果误差均未超过4.53%,网络模型的整体预测精度优于公式计算。由此得出结论,上述训练完成的网络模型的预测精度是比较高的,能够较为准确地反映EFP速度与炸药JWL状态方程参数的关系。除此以外,网络预测在保证EFP速度结果准确性的基础上,求解时间一般为几毫秒,而数值仿真时间通常需要消耗1～2 h,网络运行效率远远高于数值仿真效率。

6 结论

1) JWL状态方程各特征参数对EFP速度影响的大小顺序是:R1、A、B、R2、E0、ω。随着参数R1、R2增大,EFP速度呈指数减小趋势,减小趋势逐渐降低且趋于平稳。随着参数A、B、E0、ω增大,EFP速度呈指数增加趋势,增大趋势较为稳定。

2) 建立了考虑JWL状态方程参数的EFP速度计算公式,与仿真结果相比,绝大部分公式计算的结果误差能够控制在10%以内。

3) 搭建了考虑JWL状态方程参数的EFP速度网络预测模型,训练后的模型具有较高精度,结果误差能够控制在5%以内。网络模型求解时间仅为几毫秒,在保证高预测精度的基础上显著提高求解速度,可作为EFP速度预测、JWL状态方程参数检验校准、仿真参数调整的一种工具。