• ISSN 1001-1455  CN 51-1148/O3
  • EI、Scopus、CA、JST收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊

爆轰产物状态方程的水下爆炸反演理论研究

杨晨琛 李晓杰 闫鸿浩 王小红 王宇新

杨晨琛, 李晓杰, 闫鸿浩, 王小红, 王宇新. 爆轰产物状态方程的水下爆炸反演理论研究[J]. 爆炸与冲击, 2019, 39(9): 092201. doi: 10.11883/bzycj-2018-0210
引用本文: 杨晨琛, 李晓杰, 闫鸿浩, 王小红, 王宇新. 爆轰产物状态方程的水下爆炸反演理论研究[J]. 爆炸与冲击, 2019, 39(9): 092201. doi: 10.11883/bzycj-2018-0210
YANG Chenchen, LI Xiaojie, YAN Honghao, WANG Xiaohong, WANG Yuxin. An inverse method for the equation of state of detonation products from underwater explosion tests[J]. Explosion And Shock Waves, 2019, 39(9): 092201. doi: 10.11883/bzycj-2018-0210
Citation: YANG Chenchen, LI Xiaojie, YAN Honghao, WANG Xiaohong, WANG Yuxin. An inverse method for the equation of state of detonation products from underwater explosion tests[J]. Explosion And Shock Waves, 2019, 39(9): 092201. doi: 10.11883/bzycj-2018-0210

爆轰产物状态方程的水下爆炸反演理论研究

doi: 10.11883/bzycj-2018-0210
基金项目: 国家自然科学基金(11272081,11672067)
详细信息
    作者简介:

    杨晨琛(1993- ),男,博士研究生,yccdut@126.com

    通讯作者:

    李晓杰(1963- ),男,教授,博导,从事爆炸与冲击动力学研究,dalian03@qq.com

  • 中图分类号: O383

An inverse method for the equation of state of detonation products from underwater explosion tests

  • 摘要: 由于水下爆炸过程中爆轰产物的信息以水中压缩波形式向外传递,本文旨在探讨如何利用水下爆炸试验数据确定爆轰产物的状态方程。相较于常规圆筒试验,水下爆炸试验具有装置简单成本低、装药尺寸限制少、测定压力范围更广的特点,更适用于大药量或非理想炸药的现场测试。本文从水中冲击波轨迹和波后压力时程曲线出发,发展了由冲击波及其波后流场还原水气界面的逆特征线算法,以及根据水气界面确定爆轰产物状态方程的遗传算法。与水下爆炸正演结果对比,发现逆特征线法可以较准确地还原水气界面的位置和压力参数,有效压力下限可达2 MPa,远低于圆筒试验的测试下限0.1 GPa。根据水气界面的反演结果,遗传算法也能稳定地优化JWL方程参数,8种常用炸药的等熵衰减压力误差在爆压至0.01 GPa的区间内都小于3%。结果表明,利用本文的逆特征线算法和遗传算法,理论上可以反演出压力范围较宽、准确性较高的爆轰产物状态方程。
  • 炸药爆轰产物状态方程是爆轰物理研究中的基础问题[1]。目前,在确定状态方程的实验方法中,较为常用的是圆筒试验法(Cylinder Test)。圆筒试验法[2]源自对半球试验法[3]的改进,经过Lee等[4-5]的发展,已成为一套测定爆轰产物驱动能力、标定JWL状态方程的成熟方法,并纳入我国军用标准(GJB772A-1997)。圆筒试验法所标定的JWL方程参数,对于理想炸药,在圆筒胀裂之前的压力范围内(从爆压降至约0.1 GPa)具有较高的精度,但当药量较大或非理想成分较多时,圆筒试验法存在一定的局限性。一是圆筒试验要求炸药能制成药柱、并无缝装入标准无氧铜管,但大多数工业炸药、含铝炸药的反应区较宽,爆轰参数受装药尺寸影响较大[6-7],标准直径下的标定结果可能无法适用于更大直径的情况;二是圆筒试验需要使用高速摄影或激光干涉仪等测出筒壁的膨胀轨迹和速度曲线,其较高的成本、复杂的光路和电路,以及金属飞散物,限制了其在大药量爆轰测试中的应用。然而,在工程爆破、武器研制、海洋开发等应用领域中,大药量的非理想炸药恰恰存在着广泛的应用,在相关的爆炸问题研究、计算与设计中,急需通过实验手段确定在要求范围内准确度较高的状态方程。

    水下爆炸试验(aquarium test),或称水箱法,为测定爆轰产物状态方程提供了另一种途径。水下爆炸试验最早见于Holton[8]的研究,起初只是作为测量炸药爆压的一种手段,经过Cook[9]和Rigdon等[10]的发展和改进,由Bjarnholt总结成一套成体系的测试方法[11],随后被纳入美军军用标准(MIL-STD-1751)。Bjarnholt总结了水下爆炸试验的三个优点[12]:一是实验装置结构简单,成本较低;二是被测炸药的装药尺寸可以和现场尺寸一致,只要水域足够大;三是被测炸药和外部介质(水)充分接触,使得爆轰产物可以充分膨胀至很低的压力。简而言之,水下爆炸试验提供了一种成本低、尺寸限制少、压力范围广的爆轰产物测试途径。相对于制造等直径的无氧铜管,寻找可用于爆炸测试的水域要简单的多,约5 kg药量以下的测试可选爆炸水池,更大药量时可考虑在积水矿坑、湖泊中进行。此外,在常规圆筒试验中,爆轰产物在低压区(0.1 GPa以下)的信息随着铜管破裂而丢失[1],而在水下爆炸试验中,其信息将继续以压缩波形式向外传递,最终体现在水中流场质点的压力波动中。因此,通过水下爆炸试验测定的状态方程(并不限定是JWL方程)在低压区也能保持一定的精度,因而不仅能用于由产物膨胀主导的接触爆炸问题,在由冲击波及其波后流场主导的空中、水下非接触爆炸问题中,原理上也能保证一定的准确度。

    然而,水下爆炸法遇到的一个很大困难来自于算法。由于水下爆炸试验通常测的是近场冲击波轨迹和中远场时程压力曲线,因此相比于只需计算爆轰产物流场的圆筒试验程序,水下爆炸程序还需计算水中冲击波流场,致使后者的单次计算量远远超过前者,结果是在相同精度下的迭代寻优时间很长。有些研究者曾试图简化模型以加快单次计算,例如Itoh等根据高速摄影测得的近场水气界面位置和冲击波轨迹,使用一维流体动力学程序[13]确定了一种乳化炸药的JWL方程参数[14],但结果与圆筒试验的标定结果相差较大,说明更准确的水下爆炸程序至少应采用二维模型。由于根据试验数据确定状态方程本质上是一个反问题,即在已测数据的边界条件和模型结构的几何约束下,给出满足流场控制方程的爆轰产物状态方程。因此,本文基于以往对水的高压状态方程[15]、爆炸近场测试技术[16]和特征线正演算法[17-18]的研究,试图根据水下爆炸实验数据反演爆轰产物内部的状态。首先用特征线正演程序计算出水中冲击波及其波后流场,再以此作为初始数据,用逆特征线法复现出水气界面,最后结合遗传算法实现JWL方程参数的寻优。

    对于大药量的非理想炸药的水下爆炸测试,可行性较高的方案是连续压导探针[16]与压电传感器[19]的联合测量,分别用于采集近场的冲击波轨迹和中远场的波后压力时程曲线。尽管大药量水下爆炸是球形装药为主,但长圆柱形炸药的并联使用也广泛存在,因此本文同时考虑了这两种最常用的装药形式,图1展示了球形装药与无限长圆柱装药的半模型,可以看出两者的流场具有很大的相似性。在爆炸气泡膨胀到极限直径之前,流场中的粘性输运、热传导等扩散项可以忽略,因此球形模型对应的是一维非定常无粘流,而柱形模型对应的是二维定常无粘超声速流(本文暂不考虑亚声速情形),流场的控制方程都是双曲型偏微分方程。

    图  1  球形与柱形装药水下爆炸模型
    Figure  1.  Model of an underwater explosion of a spherical or cylindrical geometry explosive

    在进行反演计算之前,首先需将被测数据还原成流场中的数据点。对于近场冲击波轨迹,通过斜率提取可得波速,再由冲击绝热方程可补全冲击波的其余参数,其中球形模型对应正冲击波,而柱形模型对应斜冲击波。对于波后压力时程曲线,可直接还原为流场中的压力分布数据。一旦获得足够的流场数据点,使用下文的逆特征线法便可复现所覆盖区域内所有水中质点(包括水气界面)的流动状态。

    逆特征线法(inverse MOC),或称逆序特征线法,是一种在由双曲型偏微分方程控制的流场中沿时间或空间逆序推进的特征线算法,可用于根据边界条件、出流条件反推流场初始条件的反演计算,原理上是利用了双曲型方程关于时间变量或等位变量可逆的性质。与相应的正特征线法相比,逆特征线法的区别仅在于计算顺序的不同,而在精度、稳定性方面没有本质区别。与其他反演算法相比,逆特征线法继承了特征线算法的优点,即方程可降维、计算简单和几何处理能力强等。因此,在流场控制方程可看作双曲型方程的前提条件下,采用逆特征线法可求解一些困难的反演问题,如地震波数据分析中的不均匀声学介质中波系重建问题[20],以及乘波体外形设计中的给定激波的物面生成问题[21]

    在特征线理论中,特征线是变量空间中带有黎曼不变量的积分曲线,通过异族特征线的相交可以解出交点的参数。只要确定了特征线上的黎曼不变量,那么沿下游计算和沿上游计算是等价的。如图2所示,点MN为已知数据层上两点,区域Σ1Σ2分别是点M的依赖域和影响域,区域Σ3Σ4分别是点N的依赖域和影响域。那么从MN之间的已知点出发,除了无法计算区域Σ1Σ2Σ3Σ4的流态之外,既可沿下游求共同的影响域,即正演区MNP(Forward area),也可沿上游求解共同的依赖域,即反演区MNQ(Inverse area)。因此,为了将更长的水气界面纳入已知水下爆炸数据的反演区中,不仅需要冲击波数据,还需要增加波后流场的数据,此即本文采用联合测量方案的原因。

    图  2  正特征线法求解的正演区与逆特征线法求解的反演区
    Figure  2.  A forward area calculated by MOC and an inverse area calculated by inverse MOC

    以一维球对称非定常无粘流为例,其特征线方程[17]为:

    DpDt+ρcDuDt+ρc2˙εγT˙s=0(alongDrDt=u+c)
    (1)
    DpDtρcDuDt+ρc2˙εγT˙s=0(alongDrDt=uc)
    (2)

    其中:D/Dt代表沿特征线的微分;pρT分别为质点压力、密度和温度;u为质点速度,u=dr/dtc为质点声速,c2=(p/ρ)sγ与Grüneisen系数Γ相关,γ=(p/e)ρ=ρΓ˙ε代表曲率影响,˙ε=u/r=dr/(rdt),当r很大时趋于0;˙s为熵变率,˙s=ds/dt。具体推导过程见文献[17]。相较于以往的逆特征线法[22-23],该特征线方程用沿迹线的熵变代替了沿特征线的熵变,即将熵信息隐含于其余热力学量中,从而避免了对熵梯度的直接计算。

    特征线方程本质上是连续性方程和动量方程的组合,而求解流场还需结合能量方程:

    de+pd(1/ρ)=δq(alongdrdt=u)
    (3)

    其中:d/dt代表物质导数,e是质点比内能,q是沿迹线的吸热。联立式(1)~(3),通过特征线的交叉寻找解点,便可逐层推进地求解双曲型流场。如前所述,正特征线法与逆特征线法在原理上并无二致,只不过求解时推进的方向相反。

    对于给定冲击波反求物面的问题,节点的计算格式主要分为冲击波边界和内流场两类。如图3所示,两者都是先从已知点AB预估出点D以及点C的位置,再按点C的位置插值流动参数,然后根据点ABC的参数求解点D的参数,最后进行多次校正便可获得点D的解,这种预估-校正过程在理论上具有二阶精度。

    图  3  水中流场的两类反演计算格式示意图
    Figure  3.  Schematic diagram of two types of numerical scheme for the inversion of the water field

    尽管特征线算法中也存在计算网格,但不同于有限元法、有限差分中的预设网格,特征线法的网格是在计算过程中形成的,因而便于根据流场局部变化设计出具有自适应性的网格。对于本文的球形模型和柱形模型,水中流动都是涉及两个坐标的有旋流,计算每个未知点都需要两条特征线和一条迹线,为了提高网格的自适应性,计算中只保留同一族的特征线作为连续的数据层,而另一族特征线与迹线用于轮流与第一族特征线交叉而生成新的网格点[24]。由于新生网格点位于同一族特征线上,因而避免了以往的同族特征线交叉问题[25]

    图4所示,球形模型的逆特征线法求解过程为:从已知数据层如冲击波上的一点出发,沿着一条特征线向内连续延伸,通过另一条特征线定位新的交点,结合迹线计算新交点的参数,直到延伸至水气界面并形成新的界面点,然后开始新一轮循环,同时配合冲击波局部网格加密来控制新生水气界面点之间的间距。当反演计算还原了水气界面上足够的数据点,便可将水气界面的位置、压力等作为拟合目标,进行爆轰产物状态方程的参数优化。

    图  4  从冲击波及其波后流场到水气界面的自适应特征线网格
    Figure  4.  An adaptive characteristic grid from known shock wave and after-shock flow to gas-water interface

    JWL状态方程在参考等熵线(principle isentrope)上的压力函数为:

    p=AeR1V+BeR2V+CV(1+ω)
    (4)

    式中:V为无量纲的相对比容,ABCR1R2ω为待定参数。如果炸药爆轰再满足CJ条件,可以得到状态方程另需满足的三个相容方程[1]

    AeR1VJ+BeR2VJ+CV(1+ω)J=pJ
    (5)
    AR1eR1VJ+BR2eR2VJ+C(1+ω)V(2+ω)J=ρ0D2
    (6)
    AeR1VJ/R1+BeR2VJ/R2+CVωJ/ω=E0+1/2pJ(1VJ)
    (7)

    式中:pJDE0分别为爆压、爆速和初始体积内能,其中E0可通过量热计实验或热化学计算获得。因此,为了减少计算量,对于满足CJ假设的理想炸药,考虑这三个约束条件,待定参数可以减少为R1R2ω这三个参数。对于非理想炸药,如铵油炸药、乳化炸药和含铝炸药等,由于反应区较宽、能量释放较慢,第三个相容方程需进行部分修正。

    由于水中反演已经确定了水气界面在水一侧的流动参数,若将水气界面看作不掺混型两相界面,即其边界条件为压力与法向速度连续,则单次迭代计算只需涉及炸药流场即可。因此本文的优化问题可表述为:对于爆轰产物的定型膨胀问题,在产物边界的位置、法向速度已给定的约束条件下,从R1R2ω构成的三维空间中选择一点,使产物边界的压力最符合所要求的分布函数,如误差上限最小、方差最小,或其他范数准则。图5是C4炸药(R1R2ω分别为4.5,1.4和0.25)的目标函数在三个参数截面上的穷举搜索结果,可以看出该优化问题的单峰性较好。

    图  5  本文优化问题的穷举搜索结果(30万数据点)
    Figure  5.  Exhaustive search results of this optimization problem (300 000 data points)

    遗传算法(genetic algorithm)是一种求解非线性优化问题的仿生算法,基于生物进化和遗传学原理,通过模拟生物的遗传、变异和自然选择过程,在种群更替中实现对全局最优解的启发式搜索[26]。遗传算法对目标函数的限制较少,只需提供计算点对应的函数值,因而具有较高的泛用性和鲁棒性。对于JWL方程参数优化问题,由于难以构造出连续可导的目标函数,遗传算法作为一种成熟可行的方案而陆续得到应用[27-29]

    本文遗传算法的基本过程包含三部分:首先是生成初始种群,将三个待定参数R1R2ω看作不同的基因,将每一组参数{R1R2ω}看作基因的组合即染色体,以染色体作为个体的特征生成初始种群。为了增强算法的全局搜索能力,本文对基因进行二进制编码操作,即将基因的实数解空间映射到二进制编码空间,相应的编码规则见表1。一般来说,R1R2ω的取值范围是4~5,1~2和0.2~0.4,本文配合算例拓宽到3.0~8.0,0.5~3.0和0.2~0.5。

    表  1  状态方程参数的二进制编码
    Table  1.  Binary encoding of JWL EOS parameters
    参数 范围 网格数 分辨率 样本
    R1 3.00~8.11 29 (512) 0.010 010010110 (4.5)
    R2 0.500~3.05 28 (256) 0.010 01011010 (1.4)
    ω 0.200~0.515 26 (64) 0.005 001010 (0.25)
    下载: 导出CSV 
    | 显示表格

    然后是对种群进行选择,即以适应度函数为依据对种群进行优胜劣汰。适应度函数主要基于目标函数而构造,由于本文是目标函数最小化的问题,适应度函数构造为:

    f(X)=1c+Max(p/pW1)
    (6)

    式中:ppw分别对应爆轰产物正演的和水中反演的界面压力分布,Max函数的作用是取两者的最大相对误差(即本文目标函数),c为防止分母为0的一个常数。当计算了所有个体的适应度值后,需要选择用于产生后代的父本和母本,本文采用个体被选中的概率正比于适应度值的概率选择机制(轮盘赌),同时为了避免最优解丢失,采用保留最高适应度值的策略(精英保留)。

    最后是更新种群,对父本和母本按照一定的规则进行染色体交叉、基因突变,生成新的个体以恢复种群规模。交叉规则采用随机单点交叉,按照交叉概率选择基因交换点的位置,交叉概率越大,交换点位置越高,染色体变化越大。突变规则采用随机多点突变,按照突变概率反转子代基因中的每个二进制位的值。由于选取原则目前尚无统一的理论指导,本文参考Stoffa[30]等的研究,采用适中的交叉概率(0.6)与较小的突变概率(0.01)。表2为以基因R1为例的交叉和变异结果,可以看出二进制编码增强了种群的多样性,如交叉结果并不一定在父本和母本之间,而突变结果更加灵活。

    表  2  交叉操作与突变操作的示意过程
    Table  2.  Schematic process of crossover operation and mutation operation
    R1 交叉 突变
    Samlpe 1 001010101×010101010→001101010+010010101 010(0)10101→010(1)10101
    (3.85)×(4.70)→(4.06)+(4.49) (4.49)→(4.81)
    Samlpe 2 001111111×010000000→001000000+010111111 0101(1)111(1)→0101(0)111(0)
    (4.27)×(4.28)→(3.64)+(4.91) (4.91)→(4.74)
    下载: 导出CSV 
    | 显示表格

    本文遗传算法的具体流程如下:首先产生100个个体作为初始种群,然后计算所有个体的适应度值并按比例分配被选择概率以及挑选父本和母本,最后经过编码、交叉、突变和解码操作得到99个子代,加上保留的最佳个体,构成含100个新个体的新种群。重复最后两个过程,直到进化出适应度最高的个体且10代无提高。

    为了验证反演算法的有效性,先用水下爆炸正演算法[17-18]建立流场数据库,然后从中提取与实验测试结果相对应的信息,包括水中冲击波轨迹和固定位置测法的波后压力时程曲线,以此作为初始数据进行反演计算,最后比较反演结果与正演结果的差异。正演算法模型参考图1,同样基于CJ假设和近场绝热假设,采用爆轰产物的JWL状态方程以及水的Polynomial型状态方程。水气界面初始参数由水的绝热冲击方程与爆轰产物的膨胀波方程联立确定,不考虑早期的热效应,即所有水质点的卸载过程为等熵过程。由于特征线网格随计算生成,网格密度主要由CFL条件控制,在水气边界、冲击波边界附近适当加密。冲击波头最远位置为30倍的初始装药半径,最终节点数约800~1 000万。表3为本文所有炸药算例(爆压范围约5~30 GPa)的基本信息[31]

    表  3  几种炸药的爆轰参数
    Table  3.  Detonation parameters of several explosives
    炸药 ρ0/(g·cm−3) D /(km·s−1) pJ/GPa E0/(GJ·m−3)
    ANFO 0.931 4.160 5.15 2.48
    C4 1.601 8.193 28.0 9.0
    Comp B 1.717 7.980 29.5 8.5
    HNS 1.40 1.400 6.340 14.5 6.0
    LX-01 1.230 6.840 15.5 6.1
    PETN 1.50 1.500 7.450 22.0 8.56
    Tetryl 1.730 7.910 28.5 8.2
    TNT 1.630 6.930 21.0 6.0
    下载: 导出CSV 
    | 显示表格

    以装药半径0.5 m的球形C4炸药(约838.3 kg)为算例,图6为连续压导探针的测试结果示意图,包含了水中冲击波的传播轨迹,以及一小段炸药中的爆轰波传播轨迹,水中曲线上每一点的斜率对应冲击波扫过的瞬时速度,通过滤波去噪、拟合和求导处理可获得冲击波速度衰减曲线,结合冲击绝热方程可复现出冲击波阵面上的其他物理参数。图7为压电传感器的测试结果示意图,包括距离药球球心10倍、15倍、20倍半径的三个固定位置的压力时程曲线。之所以考虑这三个位置含有多方面的原因,首先,常用的电气石压电传感器的测压上限约为200~300 MPa,这三处的峰压均在此之下,都存在被测可行性;其次,冲击波强度随距离衰减,压电传感器布置越远,计算模型的绝热无粘假设越难成立,且信息丢失、外界干扰带来的误差影响越大;最后,球形炸药水下爆炸火球的极限半径约为15倍装药半径(基于Plesset & Prosperetti公式[32]),在此范围内的传感器将被损毁而增加测试成本,大多数水下爆炸测试布点常在此之外。

    图  6  连续压导探针的测试结果示意图
    Figure  6.  Schematic diagram of test results of a continuous pressure-induced conduction probe
    图  7  压电传感器的测试结果示意图
    Figure  7.  Schematic diagram of test results of a piezoelectric sensor

    图8是以上测试结果复现出的数据节点,三处测点对应的三种测试方案各有利弊,如Plan A的测试精度更高,但压力传感器是一次性的,而Plan C的传感器可重复使用,但精度会有下降,设计实验时应综合考虑。在反演计算方面,Plan C的计算规模最大、反演范围最广,误差累计比另外两者更严重。为了评估反演算法的各项性能,本文采用难度最大的Plan C的结果作为数据来源。

    图  8  水下爆炸测试的数据节点示意图
    Figure  8.  Schematic diagram of data nodes of an underwater explosion test

    表3中的基本爆轰参数,和从正演结果中提取的虚拟实验数据作为输入,从初始节点沿特征线向内反演,以逐点累积的方式求解水气界面。仍以球形C4炸药为例,图9为反演获得的水气界面结果,以及作为对比的冲击波、波后流场输入数据,其中特征线用于划分影响域和依赖域。可以看出,冲击波只对应很小的水气界面范围(R/R0=1~1.6),这说明20倍半径内的冲击波只能反映炸药膨胀过程中非常早期的信息,因而所能标定的状态方程有效范围也不大,而如果增加一条压力时程曲线的数据(本文联合测试方案),则可以明显地扩宽对炸药膨胀信息的获取范围。

    图  9  水气界面的冲击波反演区域与波后流场反演区域
    Figure  9.  Inversed area of the shock wave and after-shock flow on the gas-water interface

    图10为球形C4炸药的水气界面的迹线、压力的正演结果和反演结果,对比可以发现两者吻合度很高,不仅界面位置较为准确,而且界面上很小的压力波动也能还原出来,如球心反射造成二次峰压(2nd peak, 15.8 MPa)和三次峰压(3th peak, 3.2 MPa),表明逆特征线法可以较准确地还原水气界面的位置和压力参数。

    图  10  水气界面的位置、压力的正演结果与反演结果的比较
    Figure  10.  Comparison between forward results and inversed results of position and pressure on gas-water interface

    图11为所确定的C4炸药的爆轰产物等熵线,可以看出与JWL方程对应的标准等熵线符合程度较高,其中根据冲击波数据所确定的范围较小(v/v0<3),而根据压力时程曲线所确定的范围很宽(3<v/v0<140)。相比于圆筒试验的标定压力下限0.1 GPa,使用本文方法可以轻易突破0.01 GPa,且随着压力时程曲线的测时延长可进一步缩小。例如,当本文测时取为15 ms时,对应的压力下限已达0.002 GPa。

    图  11  C4炸药的反演等熵线与JWL参考等熵线的比较
    Figure  11.  Comparison between inverse isentrope and JWL principle isentrope of C4 explosive

    为了进一步检验所确定状态方程的精确程度。考虑到在球形炸药水下爆炸过程中,球心是压力变化最剧烈的位置:在爆轰波未入水前球心压力恒定,而当稀疏波进入爆轰产物后压力快速衰减,接下来稀疏波在球心汇聚反射又造成压力反复上升,因此球心压力适合作为比较两组状态方程的参考指标。图12为分别用所确定的C4炸药JWL方程与原方程计算的球心压力时程曲线,可以看出两者吻合良好,无论是二次、三次峰压的位置和强度,还是低压力下(小于0.01 GPa)的衰减,都具有较高的还原精度。

    图  12  用球心的时程压力曲线验证所反演的JWL方程的有效性
    Figure  12.  Configuration of inverse JWL EOS by pressure-time history curve of sphere center

    本文利用上述方法求解了8种常见炸药的JWL方程参数,如表4所示,其中标准参数取自LLNL炸药手册[31]。可以看出,在爆轰产物的压力从爆压衰减到0.01 GPa的范围内,反演结果的相对误差都在3%以下,总体上还原了爆轰产物在较宽的压力范围内的衰减情况。

    表  4  JWL方程参数的反演结果与标准数据对比
    Table  4.  Comparison between inversed results and reference parameters of JWL EOS
    实验 A/GPa B/GPa C/GPa R1 R1 ω 最大误差(p>0.01 GPa)/%
    ANFO (1) 49.46 1.89 0.474 3.90 1.10 0.333 3
    ANFO (2) 50.82 2.101 0.437 3.96 1.11 0.32 2.8
    C4 (1) 609.77 12.95 1.043 4.50 1.40 0.25
    C4 (2) 610.64 12.71 1.033 4.50 1.38 0.25 1.3
    Comp B (1) 524.23 7.678 1.082 4.20 1.10 0.34
    Comp B (2) 528.29 8.021 1.087 4.22 1.10 0.35 2.9
    HNS 1.40 (1) 366.5 6.75 1.163 4.80 1.40 0.32
    HNS 1.40 (2) 365.60 6.482 1.157 4.79 1.37 0.32 1.0
    LX-01 (1) 311.04 4.761 1.042 4.50 1.00 0.35
    LX-01 (2) 305.98 4.244 1.007 4.45 0.95 0.34 2.0
    PETN 1.50 (1) 625.3 23.29 1.149 5.25 1.60 0.28
    PETN 1.50 (2) 630.48 23.52 1.111 5.27 1.59 0.275 1.3
    Tetryl (1) 586.83 10.67 0.774 4.40 1.20 0.275
    Tetryl (2) 584.47 10.51 0.768 4.39 1.20 0.27 1.1
    TNT (1) 373.77 3.7471 0.735 4.15 0.90 0.35
    TNT (2) 376.61 4.011 0.769 4.17 0.92 0.36 1.7
     注:(1) Reference parameters;(2) Inversed results
    下载: 导出CSV 
    | 显示表格

    本文提出了一种通过水下爆炸试验反演炸药状态方程的方法,并基于逆特征线法和遗传算法开发了相应的二维计算程序,其中逆特征线法用于还原难以测量的水气界面参数,同时大幅减少后续优化模块的计算量,而遗传算法用于JWL方程参数的迭代寻优,最后通过算例验证了该方法的可行性和合理性。主要结论如下。

    (1)用水下爆炸试验反演爆轰产物膨胀过程中的信息,可行的测试对象是近场的冲击波轨迹和中远场的波后压力时程曲线,其中压力传感器的接入位置推荐在10~20倍装药半径之间,以获得较显著的压力波动范围和较小的外界影响误差。

    (2)用逆特征线算法从冲击波及其波后流场反演水气界面,可以较为准确地还原界面的位移和压力波动,包括界面后期的二次、三次峰压等细节变化。更重要的是,有效界面压力最低可达2 MPa,远小于圆筒试验的测压下限0.1 GPa,因而更适合测定爆轰产物低压区的膨胀规律。

    (3)用遗传算法从水气界面边界条件确定爆轰产物状态方程,能在合适时间内至少获得近似全局最优解。从所确定的若干种炸药的JWL状态方程来看,在0.01 GPa之前与原方程的误差都在3%以下。如果暂不考虑实验测试本身的精度丢失问题,利用本文的逆特征线反演算法和遗传优化算法,可得到压力范围较宽、准确性较高的爆轰产物状态方程。

  • 图  1  球形与柱形装药水下爆炸模型

    Figure  1.  Model of an underwater explosion of a spherical or cylindrical geometry explosive

    图  2  正特征线法求解的正演区与逆特征线法求解的反演区

    Figure  2.  A forward area calculated by MOC and an inverse area calculated by inverse MOC

    图  3  水中流场的两类反演计算格式示意图

    Figure  3.  Schematic diagram of two types of numerical scheme for the inversion of the water field

    图  4  从冲击波及其波后流场到水气界面的自适应特征线网格

    Figure  4.  An adaptive characteristic grid from known shock wave and after-shock flow to gas-water interface

    图  5  本文优化问题的穷举搜索结果(30万数据点)

    Figure  5.  Exhaustive search results of this optimization problem (300 000 data points)

    图  6  连续压导探针的测试结果示意图

    Figure  6.  Schematic diagram of test results of a continuous pressure-induced conduction probe

    图  7  压电传感器的测试结果示意图

    Figure  7.  Schematic diagram of test results of a piezoelectric sensor

    图  8  水下爆炸测试的数据节点示意图

    Figure  8.  Schematic diagram of data nodes of an underwater explosion test

    图  9  水气界面的冲击波反演区域与波后流场反演区域

    Figure  9.  Inversed area of the shock wave and after-shock flow on the gas-water interface

    图  10  水气界面的位置、压力的正演结果与反演结果的比较

    Figure  10.  Comparison between forward results and inversed results of position and pressure on gas-water interface

    图  11  C4炸药的反演等熵线与JWL参考等熵线的比较

    Figure  11.  Comparison between inverse isentrope and JWL principle isentrope of C4 explosive

    图  12  用球心的时程压力曲线验证所反演的JWL方程的有效性

    Figure  12.  Configuration of inverse JWL EOS by pressure-time history curve of sphere center

    表  1  状态方程参数的二进制编码

    Table  1.   Binary encoding of JWL EOS parameters

    参数 范围 网格数 分辨率 样本
    R1 3.00~8.11 29 (512) 0.010 010010110 (4.5)
    R2 0.500~3.05 28 (256) 0.010 01011010 (1.4)
    ω 0.200~0.515 26 (64) 0.005 001010 (0.25)
    下载: 导出CSV

    表  2  交叉操作与突变操作的示意过程

    Table  2.   Schematic process of crossover operation and mutation operation

    R1 交叉 突变
    Samlpe 1 001010101×010101010→001101010+010010101 010(0)10101→010(1)10101
    (3.85)×(4.70)→(4.06)+(4.49) (4.49)→(4.81)
    Samlpe 2 001111111×010000000→001000000+010111111 0101(1)111(1)→0101(0)111(0)
    (4.27)×(4.28)→(3.64)+(4.91) (4.91)→(4.74)
    下载: 导出CSV

    表  3  几种炸药的爆轰参数

    Table  3.   Detonation parameters of several explosives

    炸药 ρ0/(g·cm−3) D /(km·s−1) pJ/GPa E0/(GJ·m−3)
    ANFO 0.931 4.160 5.15 2.48
    C4 1.601 8.193 28.0 9.0
    Comp B 1.717 7.980 29.5 8.5
    HNS 1.40 1.400 6.340 14.5 6.0
    LX-01 1.230 6.840 15.5 6.1
    PETN 1.50 1.500 7.450 22.0 8.56
    Tetryl 1.730 7.910 28.5 8.2
    TNT 1.630 6.930 21.0 6.0
    下载: 导出CSV

    表  4  JWL方程参数的反演结果与标准数据对比

    Table  4.   Comparison between inversed results and reference parameters of JWL EOS

    实验 A/GPa B/GPa C/GPa R1 R1 ω 最大误差(p>0.01 GPa)/%
    ANFO (1) 49.46 1.89 0.474 3.90 1.10 0.333 3
    ANFO (2) 50.82 2.101 0.437 3.96 1.11 0.32 2.8
    C4 (1) 609.77 12.95 1.043 4.50 1.40 0.25
    C4 (2) 610.64 12.71 1.033 4.50 1.38 0.25 1.3
    Comp B (1) 524.23 7.678 1.082 4.20 1.10 0.34
    Comp B (2) 528.29 8.021 1.087 4.22 1.10 0.35 2.9
    HNS 1.40 (1) 366.5 6.75 1.163 4.80 1.40 0.32
    HNS 1.40 (2) 365.60 6.482 1.157 4.79 1.37 0.32 1.0
    LX-01 (1) 311.04 4.761 1.042 4.50 1.00 0.35
    LX-01 (2) 305.98 4.244 1.007 4.45 0.95 0.34 2.0
    PETN 1.50 (1) 625.3 23.29 1.149 5.25 1.60 0.28
    PETN 1.50 (2) 630.48 23.52 1.111 5.27 1.59 0.275 1.3
    Tetryl (1) 586.83 10.67 0.774 4.40 1.20 0.275
    Tetryl (2) 584.47 10.51 0.768 4.39 1.20 0.27 1.1
    TNT (1) 373.77 3.7471 0.735 4.15 0.90 0.35
    TNT (2) 376.61 4.011 0.769 4.17 0.92 0.36 1.7
     注:(1) Reference parameters;(2) Inversed results
    下载: 导出CSV
  • [1] 孙承纬. 应用爆轰物理 [M]. 北京: 国防工业出版社, 2000: 272−273.
    [2] WILKINS M, SQUIER B, HALPERIN B. The Equation of State of PBX 9404 and LX 04-01: UCRL-7797 [R]. USA: Lawrence Radiation Laboratory, 1964.
    [3] KURY J W, HORNIG H C, LEE E L, et al. Metal acceleration by chemical explosives [C] // Fourth International Symposium on Detonation. Arlington, Texas: Office of Naval Research, 1965: 3−12.
    [4] LEE E L, HORNIG H C, KURY J W. Adiabatic expansion of high explosive detonation products: UCRL-50422[R]. University of California Radiation Laboratory, 1968. DOI: 10.2172/4783904
    [5] LEE E L, FINGER M, COLLINS W. JWL equation of state coefficients for high explosives: UCID-16189 [R]. Lawrence Livermore Laboratory, 1973.
    [6] FORBES J W, LEMAR E R. Detonation wave velocity and curvature of a plastic-bonded, nonideal explosive PBXN-111 as a function of diameter and confinement [J]. Journal of Applied Physics, 1998, 84(12): 6600–6605. DOI: 10.1063/1.369033.
    [7] 韩勇, 黄辉, 黄毅民, 等. 不同直径含铝炸药的作功能力 [J]. 火炸药学报, 2008, 31(6): 5–7. DOI: 10.3969/j.issn.1007-7812.2008.06.002.

    HAN Yong, HUANG Hui, HUANG Yimin, et al. Power of aluminized explosives with different diameters [J]. Chinese Journal of Explosives & Propellants, 2008, 31(6): 5–7. DOI: 10.3969/j.issn.1007-7812.2008.06.002.
    [8] HOLTON W C. The detonation pressures in explosives as measured by transmitted shocks in water [R]. Naval Ordnance Laboratory, 1954.
    [9] COOK M A, PACK D H, MCEWAN W S. Promotion of shock initiation of detonation by metallic surfaces [J]. Transactions of the Faraday Society, 1960, 56: 1028–1038. DOI: 10.1039/tf9605601028.
    [10] RIGDON J K. Explosive performance: SANL-712-004 [R]. Office of Scientific & Technical Information Technical Reports, 1969. DOI: 10.2172/532483.
    [11] BJARNHOLT G. Suggestions on standards for measurement and data evaluation in the underwater explosion test [J]. Propellants, Explosives, Pyrotechnics, 1980, 5(2-3): 67–74. DOI: 10.1002/prep.19800050213.
    [12] BJARNHOLT G. Strength testing of explosives by underwater detonation [J]. Propellants, Explosives, Pyrotechnics, 1978, 3(1-2): 70–71. DOI: 10.1002/prep.19780030117.
    [13] MADER C L, GAGE W R. Fortran sin: a one-dimensional hydrodynamic code for problems which include chemical reactions, elastic-plastic flow, spalling, and phase transitions: LA-3720[R]. USA: Los Alamos Scientific Laboratory, 1967.
    [14] HAMASHIMA H, KATO Y, ITOH S. Determination of JWL Parameters for Non-Ideal Explosive [C] // AIP conference proceedings. AIP, 2004, 706(1): 331−334. DOI: 10.1063/1.1780246
    [15] 李晓杰, 张程娇, 王小红, 闫鸿浩. 水的状态方程对水下爆炸影响的研究 [J]. 工程力学, 2014, 31(8): 46–52. DOI: 10.6052/j.issn.1000-4750.2013.03.0180.

    LI Xiaojie, ZHANG Chengjiao, WANG Xiaohong, YAN Honghao. Numerical study on the effect of equations of state of water on underwater explosions [J]. Engineering Mechanics, 2014, 31(8): 46–52. DOI: 10.6052/j.issn.1000-4750.2013.03.0180.
    [16] 李科斌, 李晓杰, 闫鸿浩, 王小红. 一种可实现水下爆炸参数连续测量的新型电测方法 [J]. 兵工学报, 2017, 38(S1): 108–112.

    LI Kebin, LI Xiaojie, YAN Honghao, WANG Xiaohong. New electrometric method for the continuous measurement of underwater explosion parameters [J]. Acta Armamentarii, 2017, 38(S1): 108–112.
    [17] 李晓杰, 张程娇, 闫鸿浩, 等. 水下爆炸近场非均熵流的特征线差分解法 [J]. 爆炸与冲击, 2012, 32(6): 604–608. DOI: 10.3969/j.issn.1001-1455.2012.06.008.

    LI Xiaojie, ZHANG Chengjiao, YAN Honghao, et al. Difference method of characteristics in isentropic flow of underwater explosion in near-field region [J]. Explosion and Shock Waves, 2012, 32(6): 604–608. DOI: 10.3969/j.issn.1001-1455.2012.06.008.
    [18] 李晓杰, 杨晨琛, 张程娇, 闫鸿浩, 王小红. 水下爆炸非均熵二维定常流的三族特征线解法 [J]. 爆炸与冲击, 2017: 1–8. DOI: 10.11883/1001-1455(2017)01-0001-09.

    LI Xiaojie, Yang Chenchen, YAN Honghao, WANG Xiaohong. A Fdm of three characteristic lines of two-dimensional non-isentropic steady flow of cylindrical explosive underwater explosion [J]. Explosion and Shock Waves, 2017: 1–8. DOI: 10.11883/1001-1455(2017)01-0001-09.
    [19] 池家春, 马冰. TNT/RDX(40/ 60)炸药球水中爆炸波研究 [J]. 高压物理学报, 1999, 13(3): 199–204. DOI: 10.3969/j.issn.1000-5773.1999.03.008.

    CHI Jiachun, MA Bing. Underwater explosion wave by a spherical charge of composition B-3 [J]. Chinese Journal of High Pressure Physics, 1999, 13(3): 199–204. DOI: 10.3969/j.issn.1000-5773.1999.03.008.
    [20] CARRION P M. Computation of velocity and density profiles of acoustic media with vertical inhomogeneities using the method of characteristics applied to the slant stacked data [J]. The Journal of the Acoustical Society of America, 1985, 77(4): 1370–1376. DOI: 10.1121/1.392028.
    [21] SOBIECZKY H, DOUGHERTY F C, JONES K. Hypersonic waverider design from given shock waves [C] // Proceedings of the first international hypersonic waverider symposium. University of Maryland College Park, 1990: 17−19.
    [22] 钱翼稷. 超音速轴对称有旋流特征线法的计算程序 [J]. 北京航空航天大学学报, 1996, 22(4): 454–459.

    QIAN Yiji. Computer program of supersonic axisymmetric rotational characteristic method [J]. Journal of Beijing University of Aeronautics and Astronautics, 1996, 22(4): 454–459.
    [23] 乔文友, 黄国平, 夏晨, 汪明生. 发展用于高速飞行器前体/进气道匹配设计的逆特征线法 [J]. 航空动力学报, 2014, 29(06): 1444–1452.

    QIAO Wenyou, HUANG Guoping, XIA Chen, WANG Mingsheng. Development of inverse characteristic method for matching design of high-speed aircraft forebody/inlet [J]. Journal of Aerospace Power, 2014, 29(06): 1444–1452.
    [24] YANG C, LI X, ZHANG C. Numerical study of two-dimensional cylindrical underwater explosion by a modified method of characteristics [J]. Journal of Applied Physics, 2017, 122(10): 105903. DOI: 10.1063/1.4986881.
    [25] 张程娇. 炸药爆轰产物参数的特征线差分反演方法研究[D]. 大连理工大学, 2016: 101−102.
    [26] 师学明, 王家映. 地球物理资料非线性反演方法讲座(四)遗传算法 [J]. 工程地球物理学报, 2008(02): 129–140. DOI: 10.3969/j.issn.1672-7940.2008.02.001.

    SHI Xueming, WANG Jiaying. Lecture on non-linear inverse methods in geophysics (4) Genetic Algorithm Method [J]. Chinese Journal of Engineering Geophysics, 2008(02): 129–140. DOI: 10.3969/j.issn.1672-7940.2008.02.001.
    [27] 江厚满, 张若棋, 张寿齐. 用遗传算法确定材料物态方程参数 [J]. 高压物理学报, 1998(01): 48–54.

    JIANG Houman, JIANG Ruoqi, ZHANG Shouqi. Applying genetic algorithm to determine parameters in equation of state [J]. Chinese Journal of High Pressure Physics, 1998(01): 48–54.
    [28] 温丽晶, 段卓平, 张震宇, 欧卓成, 黄风雷. 采用遗传算法确定炸药爆轰产物JWL状态方程参数 [J]. 爆炸与冲击, 2013, 33(S1): 130–134.

    WEN Lijing, DUAN Zhuoping, ZHANG Zhengyu, OU Zhuocheng, HUANG Fenglei. Determination of JWL-EOS parameters for explosive detonation products using genetic algorithm [J]. Explosion and Shock Waves, 2013, 33(S1): 130–134.
    [29] 王成, 徐文龙, 郭宇飞. 基于基因遗传算法和γ律状态方程的JWL状态方程参数计算 [J]. 兵工学报, 2017, 38(S1): 167–173.

    WANG Cheng, XU Wenlong, GUO Yufei. Calculation of JWL equation of state parameters based on genetic algorithm and γ equation of state [J]. Acta Armamentarii, 2017, 38(S1): 167–173.
    [30] STOFFA P L, SEN M K. Nonlinear multiparameter optimization using genetic algorithms: Inversion of plane-wave seismograms [J]. Geophysics, 1991, 56(11): 1794–1810. DOI: 10.1190/1.1442992.
    [31] DOBRATZ B, CRAWFORD P. LLNL handbook of explosives: UCRL-52997[R]. Lawrence Livermore National Laboratory, 1985.
    [32] PLESSET M S, PROSPERETTI A. Bubble dynamics and cavitation [J]. Annual review of fluid mechanics, 1977, 9(1): 145–185. DOI: 10.1146/annurev.fl.09.010177.001045.
  • 期刊类型引用(3)

    1. 崔浩,武军安,周昊,杨永亮,邢柏阳,陈雄,郭锐. 探针式圆筒试验方法及TNT炸药爆轰产物JWL状态方程参数标定. 含能材料. 2023(07): 707-713 . 百度学术
    2. 贾新昆,卢邦飞. 工程扰动下矿山高陡边坡力学响应规律研究. 金属矿山. 2021(09): 51-59 . 百度学术
    3. 崔浩,郭锐,顾晓辉,宋浦,杨永亮,江琳,俞旸晖. BP神经网络和圆筒能量模型标定炸药的JWL参数. 火炸药学报. 2021(05): 665-673 . 百度学术

    其他类型引用(5)

  • 加载中
图(12) / 表(4)
计量
  • 文章访问数:  4994
  • HTML全文浏览量:  1588
  • PDF下载量:  86
  • 被引次数: 8
出版历程
  • 收稿日期:  2018-06-19
  • 修回日期:  2018-10-05
  • 网络出版日期:  2019-08-25
  • 刊出日期:  2019-09-01

目录

/

返回文章
返回