基于遗传算法辨识炸药JWL状态方程参数的研究

2022-05-16顾晓辉邢柏阳

振动与冲击 2022年9期

崔浩, 郭锐，宋浦，顾晓辉，周昊，邢柏阳

(1.南京理工大学机械工程学院，南京 210094;2.西安近代化学研究所燃烧与爆炸技术重点实验室，西安 710065)

炸药爆轰产物的状态方程是描述爆轰波阵面过后爆轰产物的压力、温度和密度等物理量之间关系的复杂函数，目前主要使用经验或半经验的状态方程式。其中，JWL(Jones-Wilkins-Lee)状态方程[1]可以比较精确地描述炸药爆轰产物膨胀驱动的过程，Kury等[2]提出的圆筒试验可帮助确定该方程中的未知参数。炸药爆炸产物JWL状态方程参数的确立具有重要意义，是开展爆效毁伤、武器设计和安全防护等工作的基础。

关于确定炸药爆炸产物JWL状态方程参数的方法学者们已进行了大量的研究，大多采用圆筒试验和数值仿真相结合的方法，通常需要不断手动调整参数组合直到仿真和试验数据误差小于1%[3-6]，但此种方法具有试验成本高、参数获取效率较低等缺陷。此外，还有将圆筒试验数据代入理论模型来确定状态方程参数的方法[7-9]，首先根据试验数据计算出炸药的一些爆轰参数，随后代入守恒方程中依次计算JWL状态方程中的低压、中压和高压贡献项的参数，此操作一定程度上提高了炸药JWL参数的获取效率。但由于圆筒试验存在成本高、操作繁杂、数据读取费时等缺陷，于是便有学者提出了理论拟合JWL参数的方法[10]。其中，赵峥等[11]提出了一种确定炸药JWL参数的γ拟合法，并采用差分进化法得到了4种炸药的JWL状态方程参数，发现拟合JWL状态方程和圆筒试验法得到的JWL状态方程的P-V曲线相差较小；王成等[12]基于γ律状态方程和遗传算法拟合了5种炸药的JWL参数，并将其代入数值计算中进行检验，发现拟合JWL参数具有较高的精度；沈飞等[13]基于圆筒试验中的能量转换关系，建立了一种确定炸药JWL参数的简易计算流程，发现此种方法确定的JWL状态方程的P-V曲线能够满足工程应用的需要。

本文结合遗传算法优秀的搜索能力和圆筒能量模型，提出了一种仅需已知炸药的初始密度和爆速，便可通过自编程遗传算法程序辨识炸药JWL状态方程参数的方法，而炸药的初始密度和爆速通过简单测量即可获得。同时对炸药驱动圆筒壁的过程进行数值计算，以检验辨识得到的JWL状态方程参数的有效性。

1 JWL状态方程及圆筒能量模型

1.1 JWL状态方程

等熵形式下的JWL状态方程为

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

(1)

式中：P为爆轰产物压力；V为爆轰产物相对比容，其值为爆轰产物体积和炸药初始体积之比；A、B、C、R1、R2、ω为6个待定参数。式(1)中右端依次是爆轰产物膨胀过程中高压、中压和低压3个阶段的贡献项。

由热力学关系对P积分，可得等熵线上的内能为

(2)

炸药理想爆轰时，根据C-J条件和Hugoniot关系式可得到如下3个相容方程

根据C-J条件-(∂ps/∂V)VCJ=ρ0D2可得

(3)

由爆轰产物状态的Hugoniot关系式可得

(11)

因C-J等熵线通过C-J点，则

(4)

式中：ρ0为炸药初始密度；D为炸药爆速；VCJ为C-J点处爆轰产物的相对比容；PCJ为炸药爆压；E0为单位体积炸药初始内能。相容方程组中VCJ和PCJ的值可由下式确定

(5)

(6)

式中,γ为炸药的绝热指数(或多方指数)，其值可由炸药初始密度ρ0、爆速D、爆热Q或炸药初始体积能量E0[14-15]确定

(7)

1.2 圆筒能量模型

Gurney模型[16-17]描述了受限炸药在无限长圆筒内爆炸时的能量分布，但以前提出的Gurney模型假设圆筒壁在膨胀过程中像刚体一样向外移动，圆筒壁厚度保持不变且移动速度相同，忽视了圆筒壁的塑性形变和厚度变化。为此，本文采用改进的Gurney模型[18-19]，此模型描述圆筒能量密度的精确性经过了多种乳化炸药的圆筒试验验证[20]，示意图如图1所示。

图1 炸药爆轰产物驱动圆筒壁膨胀图

如图1所示，在圆筒试验中，从一端起爆铜管内的炸药，爆轰波以爆速D在炸药中传播，爆轰波过后炸药转换为爆轰产物并向外膨胀，圆筒壁在爆轰产物的驱动下向外移动，高速扫描照相机记录下观测点处的圆筒壁径向位移。图1中L为圆筒壁径向位移，um为圆筒外壁径向移动速度，θ为圆筒壁角(初始时θ为0°，随着圆筒壁的移动逐渐增大)。R0为圆筒初始内半径，x0为圆筒初始壁厚，随着圆筒壁向外移动，圆筒内半径变为R，圆筒壁厚度变为x。

在爆轰产物膨胀过程中，其内能不断转换成圆筒壁和爆轰产物的动能，两者的动能之和通过如下方程求出

(8)

式中,ρm为圆筒材料密度，在标准圆筒试验中其值为8.93 g/cm3。

对于单位长度炸药，初始体积是以R0为半径的薄切片，而转化为爆轰产物后的体积则是以Rθ为母线、R为半径的圆锥的侧面。此时爆轰产物的相对比容为

(9)

圆筒外径径向位移可表示为

L=(R+x)-(R0+x0)

(10)

假定圆筒密度在膨胀过程中保持不变，根据质量守恒可得

(11)

结合式(10)和式(11)，式(9)可表示为

(12)

在Φ25.4 mm圆筒试验[21]中R0=12.7 mm,x0=2.6 mm。此外，文献[22]给出了计算CHNO型炸药初始体积能量的经验公式

(13)

式中,Ed(7.0)为爆轰产物相对体积为7.0时的总动能。在圆筒试验中，炸药驱动圆筒壁移动是一个绝热的过程，随着圆筒壁的移动，爆轰产物的内能仅转换为圆筒壁和爆轰产物的动能，定义圆筒试验能量模型的能量差为

(14)

对于一组JWL参数，若在爆轰产物驱动圆筒壁移动的各个时期，f(V)的值均接近0，则说明此组JWL参数适用于该炸药。然而在Φ25.4 mm圆筒试验中，一般只关注圆筒壁径向位移为6 mm和19 mm时的圆筒能量，按式(12)计算出这两处的爆轰产物相对体积分别为2.4和7.0，因此在实际计算中只需判断f(2.4)和f(7.0)是否接近0即可。此外，Castedo等研究发现圆筒位移6 mm和19 mm时的圆筒壁角分别为10.0°和11.6°，沈飞等的研究和文献[23]给出了计算CHNO型炸药在这两个位置处的圆筒外壁径向移动速度经验公式

(15)

(16)

采用本节计算方法得到了4种炸药的爆轰参数，和Urtiew等研究中的试验数据进行对比，如表1所示。观察表1可发现，除TNT的E0值外，计算结果和试验数据相差较小，两者的最大误差不超过2%。

表1 炸药爆轰参数对比

2 遗传算法辨识JWL状态方程参数过程

遗传算法[24]于1975年由Holland教授等创立，是一种基于自然选择和基因遗传学原理的随机并行搜索算法，且而不需要任何初始化信息即可寻求到全局最优解[25]。由于遗传算法不受搜索空间的约束，对目标函数没有连续、可导或单峰的要求，早已广泛应用于确定材料物态方程参数的研究[26]。

不同JWL参数组合代表了不同的圆筒模型能量差，本文的目的是寻求能使能量差最小的JWL参数组合。由于JWL方程含有3个独立的自变量，仅参数组合的多样性就使得寻优过程变得困难，而遗传算法被认为是解决具有较大、复杂搜索空间问题的最适合的搜索方法之一。本文希望借助遗传算法优秀的搜索寻优能力寻求最优的JWL参数组合，为了解决遗传算法中的“早熟”现象及收敛速度慢的问题，对选择、交叉和变异算子进行了优化，同时采用了具有全局收敛性的精英保存策略。

2.1 染色体和适应度函数

每个染色体对应一个确定的解，本文采用求解精度更高、运算效率更快的浮点数编码编辑染色体的基因信息。由3个相容方程可知，JWL状态方程中的6个参数只有3个是独立的，因此每条染色体携带3条基因。为了方便计算，本文选取取值范围确定的R1、R2和ω作为自变量，并设置R1、R2和ω的取值范围分别是4～5、1～2和0.2～0.4[27]。

染色体的适应性或一组JWL参数所确定的圆筒能量差，通过式(17)所示的适应度函数进行评估

F=|f(2.4)|+|f(7.0)|

(17)

由式(17)可知，当两处的能量差绝对值之和越小时，适应度函数值越小，表明染色体的适应性越好，染色体所携带的遗传信息越符合炸药真实的JWL状态方程。为了计算适应度函数值的大小，需在适应度子程序中先求出每条染色体(R1、R2、ω)对应的A、B、C，随后代入式(14)中分别求出特定位置的能量差。

2.2 基因操作步骤

遗传算法从一个由随机生成的染色体组成的种群开始，通过应用一系列的选择、交叉和变异等基因操作，使种群朝着适应性更好的染色体方向发展，本文使用的遗传算法详细步骤如下。

步骤1种群初始化。在自变量取值范围值随机生成80个染色体(个体)，80个染色体组成了初始种群，染色体数量的选定根据以往关于遗传算法的工作经验。

步骤2选择算子。通过选择操作，80个父代染色体将被挑选入配对库中进行交叉、变异等基因操作。本文采用锦标赛选择机制，每次从种群中随机取出sn个染色体(sn值通常取2)，然后选择其中适应性最好的一个进入配对库，重复此操作直到挑选出80个父代染色体。

锦标赛选择机制可以使适应性较好的个体具有较大的“生存”机会，并且由于只使用适应性作为选择的标准，不受适应度函数值的影响，因此避免了超级个体的影响，能在一定程度上避免过早收敛和种群停滞现象[28]。

步骤3交叉算子。交叉操作即为配对库中的80个父代染色体两两配对并按照某种方式互换基因以产生下一代。本文的交叉机制采用算数交叉，其基本思想是两个父代染色体线性组合产生出两个新的个体，子代染色体的产生机制如式(18)所示

(18)

步骤4变异算子。为了减少种群多样性的损失，降低陷入局部最优解的风险，交叉生成的子代染色体将以一定的概率发生变异。本文采用的高斯变异机制的基本流程为：当随机数小于变异概率时，随机选取染色体的某一条基因进行变异，用符合均值为μ，方差为σ2的正态分布的随机数来替换选取的基因值[29]。那么，变异后的基因值为

(19)

式中,μ=(xmin+xmax)/2,σ=(xmax-xmin)/6,xmin和xmax分别为变异基因值的上下限。

步骤5精英保存策略。为了防止选择、交叉、变异等遗传操作的随机性将当前群体中最好的个体破坏掉，本文采用具有全局收敛性的精英保存策略。当前群体中适应性最好的个体不参与交叉运算和变异运算，而是用它来替换掉本代群体经过交叉、变异等基因操作后产生的适应性最差的个体。该策略的实施是保证遗传算法收敛的重要条件之一。

步骤6种群替换。执行步骤2～步骤5后种群进化了一代，形成了新的种群，重复此步骤直到种群连续30代在适应性上没有任何提高，连续30代无进化则认为种群进化成熟。

2.3 算法流程及进化过程分析

遗传算法具体的操作流程如图2，输入炸药初始密度和爆速后，可根据式(15)和式(16)计算出u6 mm和u19 mm，并代入式(13)中计算出E0，然后根据式(7)计算出γ的值，随后可通过式(5)、式(6)计算出VCJ和PCJ的值。最后通过基因操作种群进化成熟后挑选出最优个体，读取最优个体携带的遗传信息获得R1、R2、ω的值并代入相容方程中计算A、B、C的值，得到的参数组合即为炸药的JWL参数组合。

图2 基因遗传算法辨识JWL状态方程参数流程图

图3为TNT炸药爆轰产物内能随相对体积的变化曲线，观察图中可发现JWL爆轰产物内能方程中3项的占比和相互之间关系随相对体积而变化，为了确保JWL参数在物理上的合理性，对JWL参数添加如下限定条件：

图3 TNT炸药爆轰产物比内能分布

(1)当爆轰产物相对体积大于某一特定值(通常为V>7[20，22])时，前两项之和必须小于第三项，即

(20)

(2)当爆轰产物的相对体积大于2时，第一项要小于第二项

(21)

(3)同样地，高压时，JWL方程第一项必须大于第二项(通常选取C-J状态)

(22)

图4为4种炸药最优个体的适应度函数值随进化次数的变化曲线。观察图4可知，种群进化成熟后，4种炸药最优个体的适应度值均接近0，说明进化后圆筒模型能量差很小。由于采用了精英保存策略，种群中最优个体在进化过程中无“退化”现象，遗传算法收敛。

图4 最优个体的适应度值随进化次数变化曲线

3 遗传算法辨识结果及数值检验

3.1 遗传算法辨识结果

在仅已知炸药密度和爆速的前提下，利用遗传算法程序辨识得到了4种常用炸药的JWL状态方程参数，表2列出了辨识JWL参数和标准JWL参数的对比。其中标准JWL参数取自AUTODYN材料库,均是通过圆筒试验标定，其确定的JWL状态方程可认为是炸药的标准JWL状态方程。4种炸药的JWL状态方程P-V曲线对比如图5所示，观察图5可发现，对于这4种炸药，两条P-V曲线偏差均较小，为了定量分析辨识曲线和标准曲线的相似度，引入决定系数R2

图5 辨识JWL状态方程和标准JWL状态方程的P-V曲线对比

(23)

式中：SSE为误差平方和；SST为试验数据和均值之差平方和；R2的取值范围为[0,1]，其值越接近1，表明辨识曲线越接近标准曲线。计算发现TNT、HMX、Comp B和PBX-9404的辨识JWL状态方程P-V曲线的R2值分别为0.997 7、0.993 8、0.999 3和0.998 6，其值均大于0.99，说明辨识JWL状态方程和标准JWL状态方程吻合较好。观察表2可发现，采用本文方法辨识得到的6个JWL参数均不同于标准JWL参数，但是两者定义的P-V曲线却非常接近，这是因为JWL状态方程的多项式和非线性允许不同参数组合可以同时高精度地定义同一个方程。

表2 炸药JWL状态方程参数对比

赵铮等的研究采用差分进化法拟合得到了4种炸药的JWL参数，其中有两种炸药与本文辨识的炸药相同，即TNT和HMX。将采用本文遗传算法辨识得到的JWL方程和赵铮等研究中的差分进化法拟合的JWL方程的P-V曲线进行对比，结果如图6所示。观察图6可发现，对于这两种炸药，采用遗传算法辨识得到的JWL方程P-V曲线更接近标准JWL方程P-V曲线。

(a)TNT

3.2 数值检验

为了进一步检验本文辨识方法的有效性和精度，利用显式动力学软件AUTODYN-2D对Φ25.4 mm标准圆筒试验过程进行数值计算，考虑到圆筒模型的对称性，建立二维X轴对称模型。从一端起爆装填在铜管内的炸药，起爆方式为面起爆，同时在距离起爆端200 mm圆筒外径处设置一个高斯点记录圆筒径向位移时程曲线。圆筒试验仿真涉及到的材料包括炸药、空气和铜管，其中铜管材料为OHFC无氧铜，采用Steinberg-Guinan强度模型和Mie-Grüneisen状态方程，材料参数均取自AUTODYN材料库，炸药JWL状态方程参数则分别采用辨识JWL参数和标准JWL参数。炸药和空气采用欧拉网格，空气域大小为405.0 mm×76.2 mm，铜管采用拉格朗日网格，并在仿真中设置欧拉和拉格朗日网格互相耦合。为了保证计算精度的同时提高计算速度，铜管网格大小划分成0.5 mm×0.5 mm，空气和炸药网格大小划分为0.500 mm×0.488 mm。为了防止冲击波在边界反射，空气域的边界条件设置成无反射流出Flow-out。

图7和图8分别为圆筒外壁径向位移时程曲线和圆筒外壁速度变化对比图。从图7可看出，采用辨识JWL参数和标准JWL参数计算得到的圆筒外壁径向位移曲线整体偏差较小，尤其是圆筒壁膨胀前期，实线和虚线基本重合。TNT、HMX、Comp B和PBX-9404的圆筒径向位移曲线的R2值分别为0.999 9、0.997 4、0.997 9和0.998 0，其值均接近1。图8显示对于所有炸药的圆筒外壁速度变化曲线，实线在波形、起跳点和速度值等方面较好地重现了虚线，且R2值均大于0.99。两组参数数值计算结果的高相似性证明了本文提出的利用遗传算法和圆筒能量模型辨识炸药JWL状态方程参数的有效性和高精度。