爆轟產物狀態方程的水下爆炸反演理論研究*

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

（大連理工大學工程力學系工業裝備結構分析國家重點實驗室，遼寧 大連 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%以下。如果暫不考慮實驗測試本身的精度丟失問題，利用本文的逆特征線反演算法和遺傳優化算法，可得到壓力范圍較寬、準確性較高的爆轟產物狀態方程。