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

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

（大连理工大学工程力学系工业装备结构分析国家重点实验室，辽宁 大连 116023）

炸药爆轰产物状态方程是爆轰物理研究中的基础问题[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 方程参数的寻优。

1 水下爆炸试验的测试模型

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

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

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

2 水中流场的特征线算法

2.1 逆特征线法的基本原理

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

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

2.2 逆特征线法的计算过程

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

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

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

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

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

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

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

2.3 球形水下爆炸的特征线网格

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

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

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

3 爆轰产物状态方程的确定方法

3.1 基于水气界面边界条件的优化目标

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

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

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

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

3.2 JWL 方程参数优化的遗传算法

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

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

本文遗传算法的基本过程包含三部分：首先是生成初始种群，将三个待定参数R1、R2、ω 看作不同的基因，将每一组参数{R1，R2，ω}看作基因的组合即染色体，以染色体作为个体的特征生成初始种群。为了增强算法的全局搜索能力，本文对基因进行二进制编码操作，即将基因的实数解空间映射到二进制编码空间，相应的编码规则见表1。一般来说，R1，R2，ω 的取值范围是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

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

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

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

表 2 交叉操作与突变操作的示意过程Table 2 Schematic process of crossover operation and mutation operation

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

4 算例与结果

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

表 3 几种炸药的爆轰参数Table 3 Detonation parameters of several explosives

4.1 水下爆炸测试结果

图 6 连续压导探针的测试结果示意图Fig. 6 Schematic diagram of test results of a continuous pressure-induced conduction probe

图 7 压电传感器的测试结果示意图Fig. 7 Schematic diagram of test results of a piezoelectric sensor

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

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

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

4.2 水气界面反演结果

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

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

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

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

4.3 状态方程优化结果

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

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

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

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

表 4 JWL 方程参数的反演结果与标准数据对比Table 4 Comparison between inversed results and reference parameters of JWL EOS

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

5 结 论

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

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

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

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