燃料空气炸药爆轰产物JWL状态方程参数计算

赵星宇， 白春华， 姚箭， 孙彬峰

(北京理工大学 爆炸科学与技术国家重点实验室， 北京 100081)

0 引言

液体燃料[1]、固体燃料[2-3]或固体与液体(简称固液)混合燃料[4-5]在爆炸载荷驱动[4,6-7]作用下向周围抛撒，与空气混合后能够形成无约束多相燃料空气云雾，这种可爆云雾被称之为燃料空气炸药(FAE)[8]。FAE云雾在起爆能量的引爆下能够达到爆轰状态，可造成大范围面目标毁伤[4]。

国内外学者针对不同种类、不同相态混合而成的FAE云雾爆轰开展了大量试验与理论研究，得到了云雾区爆轰参数[9]，拟合了云雾爆轰超压场随比距离的变化规律[10-12]。基于实验数据，云雾爆轰数值仿真方面的研究也相继开展，对于爆轰过程化学反应的处理目前有不同的方法。徐胜利等[13]、王晔等[14]完全忽略爆轰过程，将爆轰产物等效为高压气体，其压力、温度等参数通过理论计算或试验数据推算而得。田园等[15]、昝文涛等[16]则考虑爆轰过程的真实化学反应，该方法精确程度较高但模型计算量较大，常用于二维模型。陈明生等[17-18]利用有限元分析软件LS-DYNA建立了云雾爆轰三维模型，云雾爆轰过程采用高能炸药燃烧模型进行描述，该方法精确度和计算量适中，但需要附加爆轰产物的状态方程。

Jones-Wilkins-Lee(JWL)状态方程是描述炸药爆轰产物压力- 密度- 比能关系的方程式之一[19]，体现了爆轰产物高温高压气团的膨胀做功过程[20]。固态炸药的JWL状态方程参数一般由圆筒试验[21]确定，但由于FAE在宏观上呈云雾状态，难以固定在试验标准圆筒中，因此圆筒试验法不适用于FAE云雾。

本文采用外场云雾爆轰试验的峰值超压- 距离数据来替代圆筒试验的筒壁位移- 时间数据，通过反向传播(BP)神经网络联合遗传算法(BPNN-GA)[22-23]进行寻优计算，确定适用于无约束多相FAE云雾的JWL状态方程参数，并采用同种燃料不同布场方式的云雾爆轰试验测试数据进行仿真模型与参数的验证。

1 云雾爆轰JWL状态方程待定参数模型

JWL状态方程参数的确定一般需要基于试验建立数值仿真模型，并重复“凑参数- 数值计算- 结果对比”这一循环过程[19,24]，直到计算结果接近试验数据。因此本文依照文献[17]中的外场云雾爆轰试验进行数值建模与结果对比过程。

1.1 外场云雾爆轰数值模型

外场云雾爆轰试验使用扇形装药结构，该扇形装药结构内部装填有85 kg的液态碳氢燃料与铝粉颗粒混合物[17]。试验测试系统布置如图1所示，沿0°、90°、135°、180°这4个方向，距离扇形装药结构中心5 m、8 m、10 m、15 m、20 m、30 m、40 m和50 m的地面每环布设4个德国Kilster公司生产的Kilster-ICP型压力传感器，每个方向上布置一台美国RedLake公司生产的HG-100K型高速摄像机(拍摄速率2 000帧/s)。

图1 85 kg FAE云雾爆轰试验测试系统布置图Fig.1 Test system for 85 kg FAE cloud detonation experiment

根据数据处理结果，起爆时刻云雾形态为扁平四棱柱，其截面近似为等腰梯形。以扇形装药结构所在位置作为数值模型原点，并指定地面所在平面为计算域Oxy平面，则云雾底面4点坐标依次为(-6.0 m, 10.5 m, 0 m)、(-13.0 m, -14.5 m, 0 m)、(13.0 m, -14.5 m, 0 m)、(6.0 m, 10.5 m, 0 m)。云雾高度为3.0 m，起爆点高度为2.5 m，水平距离起爆点与原点较为接近，在此模型中可忽略，则起爆点坐标为(0 m, 0 m, 2.5 m)。考虑云雾区尺寸及试验布场中压力测点位置(最远点50 m)，建立半径R=55 m、高度H=20 m的圆柱形计算域。

根据对称性，采用1/2模型进行计算与六面体网格划分，如图2所示。

图2 85 kg FAE云雾爆轰数值仿真模型Fig.2 Numerical simulation model of 85 kg FAE cloud detonation

1.2 计算方法与边界条件

考虑不同网格尺寸对计算精度的影响，对网格单元边长为1.0 m、0.8 m、0.6 m、0.4 m这4种条件进行对比计算，取90°方向50 m测点的峰值超压Δp50对比结果见表1所示。网格尺寸的差异会导致模型同一测点计算结果的差异，随着网格尺寸的细化，网格尺寸对计算精度的影响逐渐降低。当网格尺寸从0.6 m细化至0.4 m时，计算值仅改变4.2%，综合计算精度和计算成本，选取网格尺寸为0.6 m. 模型六面体单元总数约50万。

表1 不同网格尺寸峰值超压对比Tab.1 Overpressure comparison of different mesh sizes

为简化模型，假设FAE云雾符合Chapman-Jouguet(CJ)理论[8]，即云雾的化学反应在爆轰波通过后瞬间完成并释放反应热。空气简化为符合γ定律的理想气体[18],初始压强为1×105Pa. 模型采用任意拉格朗日- 欧拉(ALE)多物质算法进行计算求解。云雾起爆时刻定义为0 ms，计算总时长为150 ms，数据保存时间间隔为0.5 ms.

计算域地面简化为刚性壁面，能够全反射冲击波，设为刚性壁面边界条件。模型Oyz平面为对称平面，对该平面节点x轴方向运动施加约束。该模型空气区边界采用无反射边界条件，并在边界单元沿内法向施加1×105Pa的压力，以模拟外界大气对计算区域的压力。

1.3 材料模型与参数

材料包括空气区材料与云雾区材料两部分。

空气采用MAT_NULL材料模型和EOS_LINEAR_POLYNOMIAL状态方程，其详细参数[18]如表2所示。状态方程表达形式为

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E,

(1)

式中：p为压力(Pa)；μ=V-1，V为比容,V=ρ/ρ0，ρ为密度(kg/m3)，ρ0为初始密度(kg/m3)；C0～C6为与材料性质相关的常数；E为能量密度(J/m3)，即单位体积内能。

表2 空气区材料参数表Tab.2 Parameters of air zone

云雾区采用高能炸药燃烧模型(MAT_HIGH_EXPLOSIVE_BURN)，模型所需参数由实验中压力传感器数据计算所得，爆速D=1 270 m/s，爆压pCJ=3.32 MPa，初始密度ρ0由(2)式估算：

(2)

式中：K为爆轰产物多方指数，此处近似取1.2. 由(2)式可得初始密度ρ0=4.53 kg/m3.

云雾状态方程采用JWL状态方程。JWL状态方程表达形式为

(3)

式中：E0为FAE云雾的初始能量密度(J/m3)，根据单位体积内云爆燃料所释放的化学能可得初始能量密度E0=2×106J/m3；A、B、R1、R2、ω为与炸药爆轰性能相关的参数，是数值模型中需要待定的参数。

2 JWL状态方程参数寻优计算

2.1 JWL参数相容关系

等熵条件下，JWL方程[19]为

ps=Ae-R1V+Be-R2V+CV-(ω+1),

(4)

式中：ps为等熵压力(Pa)；C为方程参数。根据CJ理论，等熵线、波速线、Hugoniot线相切于一点，据此得到JWL参数相容方程组[19,21]为

(5)

式中：VCJ为CJ比容，VCJ=K/(K+1).A、B、C为表征压力的参数，单位是Pa；R1、R2、ω无量纲。ρ0、D、E0、pCJ、VCJ均已计算得到，因此A、B、C、R1、R2、ω这6个参数仅有3个是独立的，当给定一组R1、R2、ω值，A、B、C可通过(5)式解出。

R1、R2、ω的取值范围[21]是4.0≤R1≤7.0、0.8≤R2≤2、0.2≤ω≤0.3，由(5)式解出A、B、C为正值的解。

2.2 构造仿真计算值与实测值偏差函数

由2.1节可知，每给定一组R1、R2、ω值，就能求解其余JWL参数，进而通过云雾爆轰数值模型进行计算，得到地面各测点处峰值压力计算值，将计算值与实验值对比，偏差越小则此JWL参数越优。因此构造如(6)式偏差函数：

(6)

式中：M为测点数；qi为权重系数；pn为数值仿真得到的第i个测点峰值压力值；pe为实测的第i个测点峰值压力值。实验共布设32个压力测点(见图1)，以爆心距离15 m为分界线，将测点分为两种情况，每种情况有16个测点。当爆心距离小于等于15 m时，测点位于云雾覆盖区内部或云雾边界处，测点压力一般为兆帕量级；当爆心距离大于15 m时，测点位于云雾覆盖区外部，超压迅速下降至10-1兆帕量级甚至更低。因此将位于5 m、8 m、10 m、15 m的测点其权重系数qi设为1，将位于20 m、30 m、40 m、50 m的测点其权重系数qi设为10，通过权重系数将云雾区内外的误差保持在同一数量级。

至此，FAE爆轰产物的JWL参数计算转化为一个数学上的多元优化问题，即在4.0≤R1≤7.0、0.8≤R2≤2.0、0.2≤ω≤0.3范围内寻找一组R1、R2、ω值，使得偏差函数Err达到最小值。

2.3 BP神经网络联合遗传算法寻优

偏差函数Err为三元函数，且形式未知，仅能通过一组R1、R2、ω值，经过数值仿真后得到峰值压力仿真值pn，进而算得Err值。智能算法[22]在解决此类复杂、多元、非线性函数的优化问题有着出色的表现，本文便采用BP神经网络根据有限组输入- 输出数据训练神经网络模型，训练完成的BP神经网络便能够建立R1、R2、ω值(输入)与Err值(输出)之间的映射关系，完成函数拟合，进一步地，采用遗传算法对训练好的BP神经网络模型寻找最小值，最小值所对应的R1、R2、ω值即为本文所求参数。

2.3.1 采用BP神经网络拟合偏差函数Err

BP神经网络是一种3层或3层以上的前馈网络，它由输入层、隐含层、输出层组成，该神经网络需要一定数量的输入- 输出数据来训练，当输出层结果不满足期望输出时，则通过BP不断调整网络的权值和阈值，以逼近训练所用的输出数据[22]。利用此特性，本文选取36组R1、R2、ω值(作为输入层)，根据JWL参数相容关系计算出其余3个参数A、B、C，再通过云雾爆轰数值模型计算相应的Err值(作为输出层)，用于训练BP神经网络。其中R1取4.0、5.0、6.0、7.0；R2取1.1、1.4、1.7；ω取0.2、0.3、0.4. 上述取值正交化成36组数据。本文所用BP神经网络拓扑结构如图3所示，采用3-7-1的3层结构，其中隐含层节点个数由2N+1规则[23]确定，N为输入层节点数，因此隐含层为7节点。

图3 BP神经网络拓扑结构Fig.3 Topology of back propagation neural network

采用数学分析软件MATLAB神经网络工具箱实现上述BP神经网络，算法及参数选择参考自文献[22]，简述如下：隐含层和输出层的传递函数为tan-Sigmoid，训练算法为Levenberg-Marquardt算法，最大训练次数为100，学习率为0.1，训练目标精度为0.000 000 4. 从上述36组数据中随机选择32组用于训练，其余4组用于预测输出，测试该网络的拟合性能。

由于BP神经网络选取训练组数据的随机性，可能导致每次计算所得的JWL状态方程参数有细微波动，因此本文经过多次重复训练，挑选其中预测输出与期望输出最为接近的网络模型如图4所示。测试结果表明，此神经网络能够预测Err函数且综合误差为0.266 0.

图4 BP神经网络测试结果Fig.4 Test result of back propagation neural network

2.3.2 遗传算法寻优

遗传算法是一类借鉴生物进化规律，即优胜劣汰遗传机制演化而来的搜索算法[22]，它采用概率化的寻优方法，不需要特定规则便能自适应地调整搜索方向，具有很好的全局寻优能力，擅长解决复杂多元函数的极值寻优问题。将2.3.1节中训练好的BP神经网络作为待寻优函数，调用MATLAB软件优化工具箱中的遗传算法寻找其最小值，其参数设置[21-22]简述如下：变量数为3(R1、R2、ω)，变量边界为[4.0, 0.8, 0.2]至[7.0, 2.0, 0.4]，种群类型为实数向量，种群数量为20，交叉概率为0.8，变异概率为0.2，最大进化代数为300.

运行遗传算法进化50代左右得到误差函数Err最小值的近似解为20.76，如图5所示，此时对应的JWL参数取值为R1=4.0、R2=1.774、ω=0.227. 再通过(5)式解得A=4.56 MPa、B=5.89 MPa、C=0.27 MPa.

图5 遗传算法寻优结果Fig.5 Optimized result of genetic algorithm

3 JWL状态方程参数对比验证

3.1 单爆源云雾爆轰实验验证

单爆源云雾爆轰实验所用装药结构横截面亦为扇形，装填有300 kg固液混合云爆燃料，燃料配比和装填密度与1.1节的扇形装药结构完全一致。从正上方拍摄的燃料分散形成云雾图像如图6(a)所示，为方便建模，将云雾截面近似成如白色虚线所示的不规则多边形，并将装药结构所在位置设为原点，经过图像处理计算，可得云雾边界形状特征坐标依次为(10.73 m, 0.38 m)、(18.49 m, 6.56 m)、(16.14 m, 10.63 m)、(6.18 m, 7.86 m)、(1.44 m, 8.53 m)、(3.59 m, 22.99 m)、(0.19 m, 22.99 m)、(-2.63 m, 9.20 m)、(-7.95 m, 7.19 m)、(-16.05 m, 15.52 m)、(-19.35 m, 13.22 m)、(-9.10 m, 2.68 m)、(-13.99 m, 1.96 m)、(-4.84 m, 16.00 m)、(3.59 m, -15.95 m)。云雾高度为7.44 m，起爆点高度为4.50 m，忽略起爆点与装药结构间的水平距离。

图6 单爆源云雾爆轰试验与数值仿真模型Fig.6 Experimental and numerical models of single-source cloud detonation

建立300 kg异形横截面云雾爆轰有限元模型如图6(b)所示，其计算域半径R=70 m，高度H=20 m，并采用2.3节所计算的该固液混合FAE云雾爆轰产物JWL状态方程参数用于该模型进行数值仿真计算。仿真压力云图与实验现场高速摄影(布设于90°方向)结果对比如图7所示，多相云雾被起爆约20 ms左右，云雾已全部完成爆轰，爆轰波从云雾上边界向外传播形成明显的空气冲击波波阵面(黑色虚线为冲击波阵面轮廓示意线)，由于云雾的特殊形状，从云雾边界透射入空气的冲击波阵面左侧略高于右侧，数值仿真结果也反映出了这一特征。

图7 起爆后20 ms实验与仿真结果对比Fig.7 Comparison of experimental and numerical results at 20 ms after ignition

选取如图6所示沿60°、120°、180° 3个方向布设在云雾爆轰区外的地面压力传感器所测量的峰值超压值，与数值计算值的对比如图8所示。由于云雾横截面的不规则形状，沿着3个方向的云雾覆盖范围各不相同，导致了沿3个方向上的超压衰减规律有所差异，在120°方向上云雾覆盖范围较小，超压衰减较快；而60°方向和180°方向云雾覆盖范围较大，超压衰减稍慢。随着距爆心距离的增加，计算值与实测值逐渐接近，在距爆心距离为50 m时，沿3个方向计算值与4发实验实测值的平均值之间相对误差分别为8.2%、9.0%、5.5%. 这表明所提出的FAE爆轰产物JWL方程参数计算方法与计算所得的参数能够模拟外场单爆源FAE爆轰。

图8 数值仿真与实测地面超压结果对比Fig.8 Comparison of numerically simulated and measured results of ground peak overpressure

3.2 多爆源云雾爆轰试验验证

多爆源云雾爆轰试验数据取自文献[18]。4个扇形装药结构同步分散所产生的云雾如图9(a)所示， 每个装药结构装填有85 kg的固液混合云爆燃料，燃料配比、装填密度与壳体结构与1.1节所述装药结构完全一致，装药结构与试验场地中心距离为22 m. 图9(b)为根据试验所建立的多爆源云雾爆轰数值仿真模型，云雾模型的尺寸与1.1节云雾尺寸相同，模型计算域半径R=70 m，高度H=20 m.

图9 多爆源云雾爆轰实验与数值仿真模型Fig.9 Experimental and numerical models of multi-source cloud detonation

纵截面压力云图与试验图像结果对比如图10所示。从图10(b)中可以看出，多爆源云雾同步爆轰产生的冲击波于模型中心处相互作用，产生局部高压区域。图11为沿两两云雾之间布设的地面压力传感器(沿45°方向)采集的峰值超压数据与数值计算值之间的对比结果，计算值与实测值在变化趋势上一致，均反映了冲击波叠加后出现高压区(15 m位置)再随传播而衰减的趋势。在距离模型中心较近区域，计算值与实测值没有完全吻合，这是由于在实验条件下，燃料由爆炸抛撒形成云雾，每个云雾爆源的形状尺寸及浓度分布都存在着差异，而数值计算模型则无法反映这些差异。在模型中心处(0 m位置)，计算结果反映出了冲击波叠加使得附近的压力显著增加。随着冲击波传播距离的增加，各云雾爆源差异对峰值超压影响逐渐下降，计算值与实测值逐渐接近，在50 m测点处，计算值与实测值偏差为11.1%. 因此本文所提出的FAE爆轰产物JWL方程参数计算方法与计算所得的参数也能够模拟外场多爆源FAE爆轰。

图10 起爆后24 ms FAE云雾爆轰试验与仿真结果对比Fig.10 Comparison of measured and numerically simulated results of FAE cloud detonation at 24 ms after ignition

图11 45°方向数值仿真与实测地面超压结果对比Fig.11 Comparison of numerically simulated and measured results of ground peak overpressure in 45° direction

4 结论

本文类比标准圆筒试验的JWL状态方程参数计算循环过程，使用外场云爆试验的峰值压力数据建立FAE爆轰产物JWL状态方程参数的确定方法。得出主要结论如下：

1)引入BP神经网络联合遗传算法简化状态方程参数优化过程，增加了寻优速度和精度。计算得到了固液混合FAE爆轰产物JWL状态方程参数，可模拟单爆源及多爆源云雾爆轰与冲击波传播过程。

2)对于单爆源云雾爆轰，计算所得的压力云图与云雾起爆后所拍摄的图像在形貌上一致，在距爆心距离50 m测点处，峰值超压的计算值与实测值偏差分别为8.2%、9.0%、5.5%.

3)对于多爆源云雾爆轰，沿45°方向的峰值超压仿真值与试验值在变化趋势上一致，在50 m测点处，计算值与实测值偏差为11.1%.