爆轟產(chǎn)物狀態(tài)方程的水下爆炸反演理論研究*

楊晨琛，李曉杰，閆鴻浩，王小紅，王宇新

（大連理工大學(xué)工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023）

炸藥爆轟產(chǎn)物狀態(tài)方程是爆轟物理研究中的基礎(chǔ)問題[1]。目前，在確定狀態(tài)方程的實(shí)驗(yàn)方法中，較為常用的是圓筒試驗(yàn)法(Cylinder Test)。圓筒試驗(yàn)法[2]源自對半球試驗(yàn)法[3]的改進(jìn)，經(jīng)過Lee 等[4-5]的發(fā)展，已成為一套測定爆轟產(chǎn)物驅(qū)動能力、標(biāo)定JWL 狀態(tài)方程的成熟方法，并納入我國軍用標(biāo)準(zhǔn)(GJB772A-1997)。圓筒試驗(yàn)法所標(biāo)定的JWL 方程參數(shù)，對于理想炸藥，在圓筒脹裂之前的壓力范圍內(nèi)（從爆壓降至約0.1 GPa）具有較高的精度，但當(dāng)藥量較大或非理想成分較多時，圓筒試驗(yàn)法存在一定的局限性。一是圓筒試驗(yàn)要求炸藥能制成藥柱、并無縫裝入標(biāo)準(zhǔn)無氧銅管，但大多數(shù)工業(yè)炸藥、含鋁炸藥的反應(yīng)區(qū)較寬，爆轟參數(shù)受裝藥尺寸影響較大[6-7]，標(biāo)準(zhǔn)直徑下的標(biāo)定結(jié)果可能無法適用于更大直徑的情況；二是圓筒試驗(yàn)需要使用高速攝影或激光干涉儀等測出筒壁的膨脹軌跡和速度曲線，其較高的成本、復(fù)雜的光路和電路，以及金屬飛散物，限制了其在大藥量爆轟測試中的應(yīng)用。然而，在工程爆破、武器研制、海洋開發(fā)等應(yīng)用領(lǐng)域中，大藥量的非理想炸藥恰恰存在著廣泛的應(yīng)用，在相關(guān)的爆炸問題研究、計算與設(shè)計中，急需通過實(shí)驗(yàn)手段確定在要求范圍內(nèi)準(zhǔn)確度較高的狀態(tài)方程。

水下爆炸試驗(yàn)(aquarium test)，或稱水箱法，為測定爆轟產(chǎn)物狀態(tài)方程提供了另一種途徑。水下爆炸試驗(yàn)最早見于Holton[8]的研究，起初只是作為測量炸藥爆壓的一種手段，經(jīng)過Cook[9]和Rigdon 等[10]的發(fā)展和改進(jìn)，由Bjarnholt 總結(jié)成一套成體系的測試方法[11]，隨后被納入美軍軍用標(biāo)準(zhǔn)(MIL-STD-1751)。Bjarnholt 總結(jié)了水下爆炸試驗(yàn)的三個優(yōu)點(diǎn)[12]：一是實(shí)驗(yàn)裝置結(jié)構(gòu)簡單，成本較低；二是被測炸藥的裝藥尺寸可以和現(xiàn)場尺寸一致，只要水域足夠大；三是被測炸藥和外部介質(zhì)（水）充分接觸，使得爆轟產(chǎn)物可以充分膨脹至很低的壓力。簡而言之，水下爆炸試驗(yàn)提供了一種成本低、尺寸限制少、壓力范圍廣的爆轟產(chǎn)物測試途徑。相對于制造等直徑的無氧銅管，尋找可用于爆炸測試的水域要簡單的多，約5 kg 藥量以下的測試可選爆炸水池，更大藥量時可考慮在積水礦坑、湖泊中進(jìn)行。此外，在常規(guī)圓筒試驗(yàn)中，爆轟產(chǎn)物在低壓區(qū)（0.1 GPa 以下）的信息隨著銅管破裂而丟失[1]，而在水下爆炸試驗(yàn)中，其信息將繼續(xù)以壓縮波形式向外傳遞，最終體現(xiàn)在水中流場質(zhì)點(diǎn)的壓力波動中。因此，通過水下爆炸試驗(yàn)測定的狀態(tài)方程（并不限定是JWL 方程）在低壓區(qū)也能保持一定的精度，因而不僅能用于由產(chǎn)物膨脹主導(dǎo)的接觸爆炸問題，在由沖擊波及其波后流場主導(dǎo)的空中、水下非接觸爆炸問題中，原理上也能保證一定的準(zhǔn)確度。

然而，水下爆炸法遇到的一個很大困難來自于算法。由于水下爆炸試驗(yàn)通常測的是近場沖擊波軌跡和中遠(yuǎn)場時程壓力曲線，因此相比于只需計算爆轟產(chǎn)物流場的圓筒試驗(yàn)程序，水下爆炸程序還需計算水中沖擊波流場，致使后者的單次計算量遠(yuǎn)遠(yuǎn)超過前者，結(jié)果是在相同精度下的迭代尋優(yōu)時間很長。有些研究者曾試圖簡化模型以加快單次計算，例如Itoh 等根據(jù)高速攝影測得的近場水氣界面位置和沖擊波軌跡，使用一維流體動力學(xué)程序[13]確定了一種乳化炸藥的JWL 方程參數(shù)[14]，但結(jié)果與圓筒試驗(yàn)的標(biāo)定結(jié)果相差較大，說明更準(zhǔn)確的水下爆炸程序至少應(yīng)采用二維模型。由于根據(jù)試驗(yàn)數(shù)據(jù)確定狀態(tài)方程本質(zhì)上是一個反問題，即在已測數(shù)據(jù)的邊界條件和模型結(jié)構(gòu)的幾何約束下，給出滿足流場控制方程的爆轟產(chǎn)物狀態(tài)方程。因此，本文基于以往對水的高壓狀態(tài)方程[15]、爆炸近場測試技術(shù)[16]和特征線正演算法[17-18]的研究，試圖根據(jù)水下爆炸實(shí)驗(yàn)數(shù)據(jù)反演爆轟產(chǎn)物內(nèi)部的狀態(tài)。首先用特征線正演程序計算出水中沖擊波及其波后流場，再以此作為初始數(shù)據(jù)，用逆特征線法復(fù)現(xiàn)出水氣界面，最后結(jié)合遺傳算法實(shí)現(xiàn)JWL 方程參數(shù)的尋優(yōu)。

1 水下爆炸試驗(yàn)的測試模型

對于大藥量的非理想炸藥的水下爆炸測試，可行性較高的方案是連續(xù)壓導(dǎo)探針[16]與壓電傳感器[19]的聯(lián)合測量，分別用于采集近場的沖擊波軌跡和中遠(yuǎn)場的波后壓力時程曲線。盡管大藥量水下爆炸是球形裝藥為主，但長圓柱形炸藥的并聯(lián)使用也廣泛存在，因此本文同時考慮了這兩種最常用的裝藥形式，圖1 展示了球形裝藥與無限長圓柱裝藥的半模型，可以看出兩者的流場具有很大的相似性。在爆炸氣泡膨脹到極限直徑之前，流場中的粘性輸運(yùn)、熱傳導(dǎo)等擴(kuò)散項(xiàng)可以忽略，因此球形模型對應(yīng)的是一維非定常無粘流，而柱形模型對應(yīng)的是二維定常無粘超聲速流（本文暫不考慮亞聲速情形），流場的控制方程都是雙曲型偏微分方程。

圖 1 球形與柱形裝藥水下爆炸模型Fig. 1 Model of an underwater explosion of a spherical or cylindrical geometry explosive

在進(jìn)行反演計算之前，首先需將被測數(shù)據(jù)還原成流場中的數(shù)據(jù)點(diǎn)。對于近場沖擊波軌跡，通過斜率提取可得波速，再由沖擊絕熱方程可補(bǔ)全沖擊波的其余參數(shù)，其中球形模型對應(yīng)正沖擊波，而柱形模型對應(yīng)斜沖擊波。對于波后壓力時程曲線，可直接還原為流場中的壓力分布數(shù)據(jù)。一旦獲得足夠的流場數(shù)據(jù)點(diǎn)，使用下文的逆特征線法便可復(fù)現(xiàn)所覆蓋區(qū)域內(nèi)所有水中質(zhì)點(diǎn)（包括水氣界面）的流動狀態(tài)。

2 水中流場的特征線算法

2.1 逆特征線法的基本原理

逆特征線法（inverse MOC），或稱逆序特征線法，是一種在由雙曲型偏微分方程控制的流場中沿時間或空間逆序推進(jìn)的特征線算法，可用于根據(jù)邊界條件、出流條件反推流場初始條件的反演計算，原理上是利用了雙曲型方程關(guān)于時間變量或等位變量可逆的性質(zhì)。與相應(yīng)的正特征線法相比，逆特征線法的區(qū)別僅在于計算順序的不同，而在精度、穩(wěn)定性方面沒有本質(zhì)區(qū)別。與其他反演算法相比，逆特征線法繼承了特征線算法的優(yōu)點(diǎn)，即方程可降維、計算簡單和幾何處理能力強(qiáng)等。因此，在流場控制方程可看作雙曲型方程的前提條件下，采用逆特征線法可求解一些困難的反演問題，如地震波數(shù)據(jù)分析中的不均勻聲學(xué)介質(zhì)中波系重建問題[20]，以及乘波體外形設(shè)計中的給定激波的物面生成問題[21]。

在特征線理論中，特征線是變量空間中帶有黎曼不變量的積分曲線，通過異族特征線的相交可以解出交點(diǎn)的參數(shù)。只要確定了特征線上的黎曼不變量，那么沿下游計算和沿上游計算是等價的。如圖2所示，點(diǎn)M、N 為已知數(shù)據(jù)層上兩點(diǎn)，區(qū)域Σ1、Σ2分別是點(diǎn)M 的依賴域和影響域，區(qū)域Σ3、Σ4分別是點(diǎn)N 的依賴域和影響域。那么從MN 之間的已知點(diǎn)出發(fā)，除了無法計算區(qū)域Σ1、Σ2、Σ3和Σ4的流態(tài)之外，既可沿下游求共同的影響域，即正演區(qū)MNP（Forward area），也可沿上游求解共同的依賴域，即反演區(qū)MNQ（Inverse area）。因此，為了將更長的水氣界面納入已知水下爆炸數(shù)據(jù)的反演區(qū)中，不僅需要沖擊波數(shù)據(jù)，還需要增加波后流場的數(shù)據(jù)，此即本文采用聯(lián)合測量方案的原因。

2.2 逆特征線法的計算過程

以一維球?qū)ΨQ非定常無粘流為例，其特征線方程[17]為：

圖 2 正特征線法求解的正演區(qū)與逆特征線法求解的反演區(qū)Fig. 2 A forward area calculated by MOC and an inverse area calculated by inverse MOC

其中：D/Dt 代表沿特征線的微分；p，ρ 和T 分別為質(zhì)點(diǎn)壓力、密度和溫度；u 為質(zhì)點(diǎn)速度，u=dr/dt ；c 為質(zhì)點(diǎn)聲速， c2=(?p/?ρ)s；γ 與Grüneisen 系數(shù)Γ 相關(guān)，γ =(?p/?e)ρ=ρΓ ；代表曲率影響，=u/r=dr/(rdt) ，當(dāng)r 很大時趨于0；為熵變率，=ds/dt 。具體推導(dǎo)過程見文獻(xiàn)[17]。相較于以往的逆特征線法[22-23]，該特征線方程用沿跡線的熵變代替了沿特征線的熵變，即將熵信息隱含于其余熱力學(xué)量中，從而避免了對熵梯度的直接計算。

特征線方程本質(zhì)上是連續(xù)性方程和動量方程的組合，而求解流場還需結(jié)合能量方程：

其中：d/dt 代表物質(zhì)導(dǎo)數(shù)，e 是質(zhì)點(diǎn)比內(nèi)能，q 是沿跡線的吸熱。聯(lián)立式(1)～(3)，通過特征線的交叉尋找解點(diǎn)，便可逐層推進(jìn)地求解雙曲型流場。如前所述，正特征線法與逆特征線法在原理上并無二致，只不過求解時推進(jìn)的方向相反。

圖 3 水中流場的兩類反演計算格式示意圖Fig. 3 Schematic diagram of two types of numerical scheme for the inversion of the water field

對于給定沖擊波反求物面的問題，節(jié)點(diǎn)的計算格式主要分為沖擊波邊界和內(nèi)流場兩類。如圖3 所示，兩者都是先從已知點(diǎn)A 和B 預(yù)估出點(diǎn)D 以及點(diǎn)C 的位置，再按點(diǎn)C 的位置插值流動參數(shù)，然后根據(jù)點(diǎn)A、B、C 的參數(shù)求解點(diǎn)D 的參數(shù)，最后進(jìn)行多次校正便可獲得點(diǎn)D 的解，這種預(yù)估-校正過程在理論上具有二階精度。

2.3 球形水下爆炸的特征線網(wǎng)格

盡管特征線算法中也存在計算網(wǎng)格，但不同于有限元法、有限差分中的預(yù)設(shè)網(wǎng)格，特征線法的網(wǎng)格是在計算過程中形成的，因而便于根據(jù)流場局部變化設(shè)計出具有自適應(yīng)性的網(wǎng)格。對于本文的球形模型和柱形模型，水中流動都是涉及兩個坐標(biāo)的有旋流，計算每個未知點(diǎn)都需要兩條特征線和一條跡線，為了提高網(wǎng)格的自適應(yīng)性，計算中只保留同一族的特征線作為連續(xù)的數(shù)據(jù)層，而另一族特征線與跡線用于輪流與第一族特征線交叉而生成新的網(wǎng)格點(diǎn)[24]。由于新生網(wǎng)格點(diǎn)位于同一族特征線上，因而避免了以往的同族特征線交叉問題[25]。

圖 4 從沖擊波及其波后流場到水氣界面的自適應(yīng)特征線網(wǎng)格Fig. 4 An adaptive characteristic grid from known shock wave and after-shock flow to gas-water interface

如圖4 所示，球形模型的逆特征線法求解過程為：從已知數(shù)據(jù)層如沖擊波上的一點(diǎn)出發(fā)，沿著一條特征線向內(nèi)連續(xù)延伸，通過另一條特征線定位新的交點(diǎn)，結(jié)合跡線計算新交點(diǎn)的參數(shù)，直到延伸至水氣界面并形成新的界面點(diǎn)，然后開始新一輪循環(huán)，同時配合沖擊波局部網(wǎng)格加密來控制新生水氣界面點(diǎn)之間的間距。當(dāng)反演計算還原了水氣界面上足夠的數(shù)據(jù)點(diǎn)，便可將水氣界面的位置、壓力等作為擬合目標(biāo)，進(jìn)行爆轟產(chǎn)物狀態(tài)方程的參數(shù)優(yōu)化。

3 爆轟產(chǎn)物狀態(tài)方程的確定方法

3.1 基于水氣界面邊界條件的優(yōu)化目標(biāo)

JWL 狀態(tài)方程在參考等熵線(principle isentrope)上的壓力函數(shù)為：

式中：V 為無量綱的相對比容，A，B，C，R1，R2和ω 為待定參數(shù)。如果炸藥爆轟再滿足CJ 條件，可以得到狀態(tài)方程另需滿足的三個相容方程[1]：

式中：pJ、D 和E0分別為爆壓、爆速和初始體積內(nèi)能，其中E0可通過量熱計實(shí)驗(yàn)或熱化學(xué)計算獲得。因此，為了減少計算量，對于滿足CJ 假設(shè)的理想炸藥，考慮這三個約束條件，待定參數(shù)可以減少為R1，R2和ω 這三個參數(shù)。對于非理想炸藥，如銨油炸藥、乳化炸藥和含鋁炸藥等，由于反應(yīng)區(qū)較寬、能量釋放較慢，第三個相容方程需進(jìn)行部分修正。

由于水中反演已經(jīng)確定了水氣界面在水一側(cè)的流動參數(shù)，若將水氣界面看作不摻混型兩相界面，即其邊界條件為壓力與法向速度連續(xù)，則單次迭代計算只需涉及炸藥流場即可。因此本文的優(yōu)化問題可表述為：對于爆轟產(chǎn)物的定型膨脹問題，在產(chǎn)物邊界的位置、法向速度已給定的約束條件下，從R1、R2和ω 構(gòu)成的三維空間中選擇一點(diǎn)，使產(chǎn)物邊界的壓力最符合所要求的分布函數(shù)，如誤差上限最小、方差最小，或其他范數(shù)準(zhǔn)則。圖5 是C4 炸藥（R1，R2和ω 分別為4.5，1.4 和0.25）的目標(biāo)函數(shù)在三個參數(shù)截面上的窮舉搜索結(jié)果，可以看出該優(yōu)化問題的單峰性較好。

3.2 JWL 方程參數(shù)優(yōu)化的遺傳算法

遺傳算法(genetic algorithm)是一種求解非線性優(yōu)化問題的仿生算法，基于生物進(jìn)化和遺傳學(xué)原理，通過模擬生物的遺傳、變異和自然選擇過程，在種群更替中實(shí)現(xiàn)對全局最優(yōu)解的啟發(fā)式搜索[26]。遺傳算法對目標(biāo)函數(shù)的限制較少，只需提供計算點(diǎn)對應(yīng)的函數(shù)值，因而具有較高的泛用性和魯棒性。對于JWL 方程參數(shù)優(yōu)化問題，由于難以構(gòu)造出連續(xù)可導(dǎo)的目標(biāo)函數(shù)，遺傳算法作為一種成熟可行的方案而陸續(xù)得到應(yīng)用[27-29]。

圖 5 本文優(yōu)化問題的窮舉搜索結(jié)果（30 萬數(shù)據(jù)點(diǎn)）Fig. 5 Exhaustive search results of this optimization problem(300 000 data points)

本文遺傳算法的基本過程包含三部分：首先是生成初始種群，將三個待定參數(shù)R1、R2、ω 看作不同的基因，將每一組參數(shù){R1，R2，ω}看作基因的組合即染色體，以染色體作為個體的特征生成初始種群。為了增強(qiáng)算法的全局搜索能力，本文對基因進(jìn)行二進(jìn)制編碼操作，即將基因的實(shí)數(shù)解空間映射到二進(jìn)制編碼空間，相應(yīng)的編碼規(guī)則見表1。一般來說，R1，R2，ω 的取值范圍是4～5，1～2 和0.2～0.4，本文配合算例拓寬到3.0～8.0，0.5～3.0 和0.2～0.5。

表 1 狀態(tài)方程參數(shù)的二進(jìn)制編碼Table 1 Binary encoding of JWL EOS parameters

然后是對種群進(jìn)行選擇，即以適應(yīng)度函數(shù)為依據(jù)對種群進(jìn)行優(yōu)勝劣汰。適應(yīng)度函數(shù)主要基于目標(biāo)函數(shù)而構(gòu)造，由于本文是目標(biāo)函數(shù)最小化的問題，適應(yīng)度函數(shù)構(gòu)造為：

式中：p 和pw分別對應(yīng)爆轟產(chǎn)物正演的和水中反演的界面壓力分布，Max 函數(shù)的作用是取兩者的最大相對誤差（即本文目標(biāo)函數(shù)），c 為防止分母為0 的一個常數(shù)。當(dāng)計算了所有個體的適應(yīng)度值后，需要選擇用于產(chǎn)生后代的父本和母本，本文采用個體被選中的概率正比于適應(yīng)度值的概率選擇機(jī)制（輪盤賭），同時為了避免最優(yōu)解丟失，采用保留最高適應(yīng)度值的策略（精英保留）。

最后是更新種群，對父本和母本按照一定的規(guī)則進(jìn)行染色體交叉、基因突變，生成新的個體以恢復(fù)種群規(guī)模。交叉規(guī)則采用隨機(jī)單點(diǎn)交叉，按照交叉概率選擇基因交換點(diǎn)的位置，交叉概率越大，交換點(diǎn)位置越高，染色體變化越大。突變規(guī)則采用隨機(jī)多點(diǎn)突變，按照突變概率反轉(zhuǎn)子代基因中的每個二進(jìn)制位的值。由于選取原則目前尚無統(tǒng)一的理論指導(dǎo)，本文參考Stoffa[30]等的研究，采用適中的交叉概率（0.6）與較小的突變概率（0.01）。表2 為以基因R1為例的交叉和變異結(jié)果，可以看出二進(jìn)制編碼增強(qiáng)了種群的多樣性，如交叉結(jié)果并不一定在父本和母本之間，而突變結(jié)果更加靈活。

表 2 交叉操作與突變操作的示意過程Table 2 Schematic process of crossover operation and mutation operation

本文遺傳算法的具體流程如下：首先產(chǎn)生100 個個體作為初始種群，然后計算所有個體的適應(yīng)度值并按比例分配被選擇概率以及挑選父本和母本，最后經(jīng)過編碼、交叉、突變和解碼操作得到99 個子代，加上保留的最佳個體，構(gòu)成含100 個新個體的新種群。重復(fù)最后兩個過程，直到進(jìn)化出適應(yīng)度最高的個體且10 代無提高。

4 算例與結(jié)果

為了驗(yàn)證反演算法的有效性，先用水下爆炸正演算法[17-18]建立流場數(shù)據(jù)庫，然后從中提取與實(shí)驗(yàn)測試結(jié)果相對應(yīng)的信息，包括水中沖擊波軌跡和固定位置測法的波后壓力時程曲線，以此作為初始數(shù)據(jù)進(jìn)行反演計算，最后比較反演結(jié)果與正演結(jié)果的差異。正演算法模型參考圖1，同樣基于CJ 假設(shè)和近場絕熱假設(shè)，采用爆轟產(chǎn)物的JWL 狀態(tài)方程以及水的Polynomial 型狀態(tài)方程。水氣界面初始參數(shù)由水的絕熱沖擊方程與爆轟產(chǎn)物的膨脹波方程聯(lián)立確定，不考慮早期的熱效應(yīng)，即所有水質(zhì)點(diǎn)的卸載過程為等熵過程。由于特征線網(wǎng)格隨計算生成，網(wǎng)格密度主要由CFL 條件控制，在水氣邊界、沖擊波邊界附近適當(dāng)加密。沖擊波頭最遠(yuǎn)位置為30 倍的初始裝藥半徑，最終節(jié)點(diǎn)數(shù)約800～1 000 萬。表3 為本文所有炸藥算例（爆壓范圍約5～30 GPa）的基本信息[31]。

表 3 幾種炸藥的爆轟參數(shù)Table 3 Detonation parameters of several explosives

4.1 水下爆炸測試結(jié)果

圖 6 連續(xù)壓導(dǎo)探針的測試結(jié)果示意圖Fig. 6 Schematic diagram of test results of a continuous pressure-induced conduction probe

圖 7 壓電傳感器的測試結(jié)果示意圖Fig. 7 Schematic diagram of test results of a piezoelectric sensor

以裝藥半徑0.5 m 的球形C4 炸藥（約838.3 kg）為算例，圖6 為連續(xù)壓導(dǎo)探針的測試結(jié)果示意圖，包含了水中沖擊波的傳播軌跡，以及一小段炸藥中的爆轟波傳播軌跡，水中曲線上每一點(diǎn)的斜率對應(yīng)沖擊波掃過的瞬時速度，通過濾波去噪、擬合和求導(dǎo)處理可獲得沖擊波速度衰減曲線，結(jié)合沖擊絕熱方程可復(fù)現(xiàn)出沖擊波陣面上的其他物理參數(shù)。圖7 為壓電傳感器的測試結(jié)果示意圖，包括距離藥球球心10 倍、15 倍、20 倍半徑的三個固定位置的壓力時程曲線。之所以考慮這三個位置含有多方面的原因，首先，常用的電氣石壓電傳感器的測壓上限約為200～300 MPa，這三處的峰壓均在此之下，都存在被測可行性；其次，沖擊波強(qiáng)度隨距離衰減，壓電傳感器布置越遠(yuǎn)，計算模型的絕熱無粘假設(shè)越難成立，且信息丟失、外界干擾帶來的誤差影響越大；最后，球形炸藥水下爆炸火球的極限半徑約為15 倍裝藥半徑（基于Plesset& Prosperetti 公式[32]），在此范圍內(nèi)的傳感器將被損毀而增加測試成本，大多數(shù)水下爆炸測試布點(diǎn)常在此之外。

圖 8 水下爆炸測試的數(shù)據(jù)節(jié)點(diǎn)示意圖Fig. 8 Schematic diagram of data nodes of an underwater explosion test

圖8 是以上測試結(jié)果復(fù)現(xiàn)出的數(shù)據(jù)節(jié)點(diǎn)，三處測點(diǎn)對應(yīng)的三種測試方案各有利弊，如Plan A 的測試精度更高，但壓力傳感器是一次性的，而Plan C 的傳感器可重復(fù)使用，但精度會有下降，設(shè)計實(shí)驗(yàn)時應(yīng)綜合考慮。在反演計算方面，Plan C 的計算規(guī)模最大、反演范圍最廣，誤差累計比另外兩者更嚴(yán)重。為了評估反演算法的各項(xiàng)性能，本文采用難度最大的Plan C 的結(jié)果作為數(shù)據(jù)來源。

4.2 水氣界面反演結(jié)果

以表3 中的基本爆轟參數(shù)，和從正演結(jié)果中提取的虛擬實(shí)驗(yàn)數(shù)據(jù)作為輸入，從初始節(jié)點(diǎn)沿特征線向內(nèi)反演，以逐點(diǎn)累積的方式求解水氣界面。仍以球形C4 炸藥為例，圖9 為反演獲得的水氣界面結(jié)果，以及作為對比的沖擊波、波后流場輸入數(shù)據(jù)，其中特征線用于劃分影響域和依賴域。可以看出，沖擊波只對應(yīng)很小的水氣界面范圍（R/R0=1～1.6），這說明20 倍半徑內(nèi)的沖擊波只能反映炸藥膨脹過程中非常早期的信息，因而所能標(biāo)定的狀態(tài)方程有效范圍也不大，而如果增加一條壓力時程曲線的數(shù)據(jù)（本文聯(lián)合測試方案），則可以明顯地擴(kuò)寬對炸藥膨脹信息的獲取范圍。

圖10 為球形C4 炸藥的水氣界面的跡線、壓力的正演結(jié)果和反演結(jié)果，對比可以發(fā)現(xiàn)兩者吻合度很高，不僅界面位置較為準(zhǔn)確，而且界面上很小的壓力波動也能還原出來，如球心反射造成二次峰壓（2nd peak, 15.8 MPa）和三次峰壓（3th peak, 3.2 MPa），表明逆特征線法可以較準(zhǔn)確地還原水氣界面的位置和壓力參數(shù)。

圖 9 水氣界面的沖擊波反演區(qū)域與波后流場反演區(qū)域Fig. 9 Inversed area of the shock wave and after-shock flow on the gas-water interface

圖 10 水氣界面的位置、壓力的正演結(jié)果與反演結(jié)果的比較Fig. 10 Comparison between forward results and inversed results of position and pressure on gas-water interface

4.3 狀態(tài)方程優(yōu)化結(jié)果

圖11 為所確定的C4 炸藥的爆轟產(chǎn)物等熵線，可以看出與JWL 方程對應(yīng)的標(biāo)準(zhǔn)等熵線符合程度較高，其中根據(jù)沖擊波數(shù)據(jù)所確定的范圍較小（v/v0＜3），而根據(jù)壓力時程曲線所確定的范圍很寬（3＜v/v0＜140）。相比于圓筒試驗(yàn)的標(biāo)定壓力下限0.1 GPa，使用本文方法可以輕易突破0.01 GPa，且隨著壓力時程曲線的測時延長可進(jìn)一步縮小。例如，當(dāng)本文測時取為15 ms 時，對應(yīng)的壓力下限已達(dá)0.002 GPa。

為了進(jìn)一步檢驗(yàn)所確定狀態(tài)方程的精確程度?？紤]到在球形炸藥水下爆炸過程中，球心是壓力變化最劇烈的位置：在爆轟波未入水前球心壓力恒定，而當(dāng)稀疏波進(jìn)入爆轟產(chǎn)物后壓力快速衰減，接下來稀疏波在球心匯聚反射又造成壓力反復(fù)上升，因此球心壓力適合作為比較兩組狀態(tài)方程的參考指標(biāo)。圖12 為分別用所確定的C4 炸藥JWL 方程與原方程計算的球心壓力時程曲線，可以看出兩者吻合良好，無論是二次、三次峰壓的位置和強(qiáng)度，還是低壓力下(小于0.01 GPa)的衰減，都具有較高的還原精度。

圖 11 C4 炸藥的反演等熵線與JWL 參考等熵線的比較Fig. 11 Comparison between inverse isentrope and JWL principle isentrope of C4 explosive

圖 12 用球心的時程壓力曲線驗(yàn)證所反演的JWL方程的有效性Fig. 12 Configuration of inverse JWL EOS by pressure-time history curve of sphere center

表 4 JWL 方程參數(shù)的反演結(jié)果與標(biāo)準(zhǔn)數(shù)據(jù)對比Table 4 Comparison between inversed results and reference parameters of JWL EOS

本文利用上述方法求解了8 種常見炸藥的JWL 方程參數(shù)，如表4 所示，其中標(biāo)準(zhǔn)參數(shù)取自LLNL 炸藥手冊[31]。可以看出，在爆轟產(chǎn)物的壓力從爆壓衰減到0.01 GPa 的范圍內(nèi)，反演結(jié)果的相對誤差都在3%以下，總體上還原了爆轟產(chǎn)物在較寬的壓力范圍內(nèi)的衰減情況。

5 結(jié) 論

本文提出了一種通過水下爆炸試驗(yàn)反演炸藥狀態(tài)方程的方法，并基于逆特征線法和遺傳算法開發(fā)了相應(yīng)的二維計算程序，其中逆特征線法用于還原難以測量的水氣界面參數(shù)，同時大幅減少后續(xù)優(yōu)化模塊的計算量，而遺傳算法用于JWL 方程參數(shù)的迭代尋優(yōu)，最后通過算例驗(yàn)證了該方法的可行性和合理性。主要結(jié)論如下。

（1）用水下爆炸試驗(yàn)反演爆轟產(chǎn)物膨脹過程中的信息，可行的測試對象是近場的沖擊波軌跡和中遠(yuǎn)場的波后壓力時程曲線，其中壓力傳感器的接入位置推薦在10～20 倍裝藥半徑之間，以獲得較顯著的壓力波動范圍和較小的外界影響誤差。

（2）用逆特征線算法從沖擊波及其波后流場反演水氣界面，可以較為準(zhǔn)確地還原界面的位移和壓力波動，包括界面后期的二次、三次峰壓等細(xì)節(jié)變化。更重要的是，有效界面壓力最低可達(dá)2 MPa，遠(yuǎn)小于圓筒試驗(yàn)的測壓下限0.1 GPa，因而更適合測定爆轟產(chǎn)物低壓區(qū)的膨脹規(guī)律。

（3）用遺傳算法從水氣界面邊界條件確定爆轟產(chǎn)物狀態(tài)方程，能在合適時間內(nèi)至少獲得近似全局最優(yōu)解。從所確定的若干種炸藥的JWL 狀態(tài)方程來看，在0.01 GPa 之前與原方程的誤差都在3%以下。如果暫不考慮實(shí)驗(yàn)測試本身的精度丟失問題，利用本文的逆特征線反演算法和遺傳優(yōu)化算法，可得到壓力范圍較寬、準(zhǔn)確性較高的爆轟產(chǎn)物狀態(tài)方程。