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

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

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

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

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

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

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

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

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

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

2 水中流場(chǎng)的特征線算法

2.1 逆特征線法的基本原理

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

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

2.2 逆特征線法的計(jì)算過(guò)程

以一維球?qū)ΨQ(chēng)非定常無(wú)粘流為例，其特征線方程[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 很大時(shí)趨于0；為熵變率，=ds/dt 。具體推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[17]。相較于以往的逆特征線法[22-23]，該特征線方程用沿跡線的熵變代替了沿特征線的熵變，即將熵信息隱含于其余熱力學(xué)量中，從而避免了對(duì)熵梯度的直接計(jì)算。

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

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

圖 3 水中流場(chǎng)的兩類(lèi)反演計(jì)算格式示意圖Fig. 3 Schematic diagram of two types of numerical scheme for the inversion of the water field

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

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

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

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

如圖4 所示，球形模型的逆特征線法求解過(guò)程為：從已知數(shù)據(jù)層如沖擊波上的一點(diǎn)出發(fā)，沿著一條特征線向內(nèi)連續(xù)延伸，通過(guò)另一條特征線定位新的交點(diǎn)，結(jié)合跡線計(jì)算新交點(diǎn)的參數(shù)，直到延伸至水氣界面并形成新的界面點(diǎn)，然后開(kāi)始新一輪循環(huán)，同時(shí)配合沖擊波局部網(wǎng)格加密來(lái)控制新生水氣界面點(diǎn)之間的間距。當(dāng)反演計(jì)算還原了水氣界面上足夠的數(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 為無(wú)量綱的相對(duì)比容，A，B，C，R1，R2和ω 為待定參數(shù)。如果炸藥爆轟再滿足CJ 條件，可以得到狀態(tài)方程另需滿足的三個(gè)相容方程[1]：

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

由于水中反演已經(jīng)確定了水氣界面在水一側(cè)的流動(dòng)參數(shù)，若將水氣界面看作不摻混型兩相界面，即其邊界條件為壓力與法向速度連續(xù)，則單次迭代計(jì)算只需涉及炸藥流場(chǎng)即可。因此本文的優(yōu)化問(wèn)題可表述為：對(duì)于爆轟產(chǎn)物的定型膨脹問(wè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ù)在三個(gè)參數(shù)截面上的窮舉搜索結(jié)果，可以看出該優(yōu)化問(wèn)題的單峰性較好。

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

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

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

本文遺傳算法的基本過(guò)程包含三部分：首先是生成初始種群，將三個(gè)待定參數(shù)R1、R2、ω 看作不同的基因，將每一組參數(shù){R1，R2，ω}看作基因的組合即染色體，以染色體作為個(gè)體的特征生成初始種群。為了增強(qiáng)算法的全局搜索能力，本文對(duì)基因進(jìn)行二進(jìn)制編碼操作，即將基因的實(shí)數(shù)解空間映射到二進(jìn)制編碼空間，相應(yīng)的編碼規(guī)則見(jiàn)表1。一般來(lái)說(shuō)，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

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

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

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

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

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

4 算例與結(jié)果

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

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

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

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

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

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

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

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

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

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

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

圖 9 水氣界面的沖擊波反演區(qū)域與波后流場(chǎng)反演區(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 方程對(duì)應(yīng)的標(biāo)準(zhǔn)等熵線符合程度較高，其中根據(jù)沖擊波數(shù)據(jù)所確定的范圍較?。╲/v0＜3），而根據(jù)壓力時(shí)程曲線所確定的范圍很寬（3＜v/v0＜140）。相比于圓筒試驗(yàn)的標(biāo)定壓力下限0.1 GPa，使用本文方法可以輕易突破0.01 GPa，且隨著壓力時(shí)程曲線的測(cè)時(shí)延長(zhǎng)可進(jìn)一步縮小。例如，當(dāng)本文測(cè)時(shí)取為15 ms 時(shí)，對(duì)應(yīng)的壓力下限已達(dá)0.002 GPa。

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

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

圖 12 用球心的時(shí)程壓力曲線驗(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ù)對(duì)比Table 4 Comparison between inversed results and reference parameters of JWL EOS

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

5 結(jié) 論

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

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

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

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