燃料空氣炸藥爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)計算

趙星宇， 白春華， 姚箭， 孫彬峰

(北京理工大學 爆炸科學與技術(shù)國家重點實驗室， 北京 100081)

0 引言

液體燃料[1]、固體燃料[2-3]或固體與液體(簡稱固液)混合燃料[4-5]在爆炸載荷驅(qū)動[4,6-7]作用下向周圍拋撒，與空氣混合后能夠形成無約束多相燃料空氣云霧，這種可爆云霧被稱之為燃料空氣炸藥(FAE)[8]。FAE云霧在起爆能量的引爆下能夠達到爆轟狀態(tài)，可造成大范圍面目標毀傷[4]。

國內(nèi)外學者針對不同種類、不同相態(tài)混合而成的FAE云霧爆轟開展了大量試驗與理論研究，得到了云霧區(qū)爆轟參數(shù)[9]，擬合了云霧爆轟超壓場隨比距離的變化規(guī)律[10-12]。基于實驗數(shù)據(jù)，云霧爆轟數(shù)值仿真方面的研究也相繼開展，對于爆轟過程化學反應(yīng)的處理目前有不同的方法。徐勝利等[13]、王曄等[14]完全忽略爆轟過程，將爆轟產(chǎn)物等效為高壓氣體，其壓力、溫度等參數(shù)通過理論計算或試驗數(shù)據(jù)推算而得。田園等[15]、昝文濤等[16]則考慮爆轟過程的真實化學反應(yīng)，該方法精確程度較高但模型計算量較大，常用于二維模型。陳明生等[17-18]利用有限元分析軟件LS-DYNA建立了云霧爆轟三維模型，云霧爆轟過程采用高能炸藥燃燒模型進行描述，該方法精確度和計算量適中，但需要附加爆轟產(chǎn)物的狀態(tài)方程。

Jones-Wilkins-Lee(JWL)狀態(tài)方程是描述炸藥爆轟產(chǎn)物壓力- 密度- 比能關(guān)系的方程式之一[19]，體現(xiàn)了爆轟產(chǎn)物高溫高壓氣團的膨脹做功過程[20]。固態(tài)炸藥的JWL狀態(tài)方程參數(shù)一般由圓筒試驗[21]確定，但由于FAE在宏觀上呈云霧狀態(tài)，難以固定在試驗標準圓筒中，因此圓筒試驗法不適用于FAE云霧。

本文采用外場云霧爆轟試驗的峰值超壓- 距離數(shù)據(jù)來替代圓筒試驗的筒壁位移- 時間數(shù)據(jù)，通過反向傳播(BP)神經(jīng)網(wǎng)絡(luò)聯(lián)合遺傳算法(BPNN-GA)[22-23]進行尋優(yōu)計算，確定適用于無約束多相FAE云霧的JWL狀態(tài)方程參數(shù)，并采用同種燃料不同布場方式的云霧爆轟試驗測試數(shù)據(jù)進行仿真模型與參數(shù)的驗證。

1 云霧爆轟JWL狀態(tài)方程待定參數(shù)模型

JWL狀態(tài)方程參數(shù)的確定一般需要基于試驗建立數(shù)值仿真模型，并重復“湊參數(shù)- 數(shù)值計算- 結(jié)果對比”這一循環(huán)過程[19,24]，直到計算結(jié)果接近試驗數(shù)據(jù)。因此本文依照文獻[17]中的外場云霧爆轟試驗進行數(shù)值建模與結(jié)果對比過程。

1.1 外場云霧爆轟數(shù)值模型

外場云霧爆轟試驗使用扇形裝藥結(jié)構(gòu)，該扇形裝藥結(jié)構(gòu)內(nèi)部裝填有85 kg的液態(tài)碳氫燃料與鋁粉顆粒混合物[17]。試驗測試系統(tǒng)布置如圖1所示，沿0°、90°、135°、180°這4個方向，距離扇形裝藥結(jié)構(gòu)中心5 m、8 m、10 m、15 m、20 m、30 m、40 m和50 m的地面每環(huán)布設(shè)4個德國Kilster公司生產(chǎn)的Kilster-ICP型壓力傳感器，每個方向上布置一臺美國RedLake公司生產(chǎn)的HG-100K型高速攝像機(拍攝速率2 000幀/s)。

圖1 85 kg FAE云霧爆轟試驗測試系統(tǒng)布置圖Fig.1 Test system for 85 kg FAE cloud detonation experiment

根據(jù)數(shù)據(jù)處理結(jié)果，起爆時刻云霧形態(tài)為扁平四棱柱，其截面近似為等腰梯形。以扇形裝藥結(jié)構(gòu)所在位置作為數(shù)值模型原點，并指定地面所在平面為計算域Oxy平面，則云霧底面4點坐標依次為(-6.0 m, 10.5 m, 0 m)、(-13.0 m, -14.5 m, 0 m)、(13.0 m, -14.5 m, 0 m)、(6.0 m, 10.5 m, 0 m)。云霧高度為3.0 m，起爆點高度為2.5 m，水平距離起爆點與原點較為接近，在此模型中可忽略，則起爆點坐標為(0 m, 0 m, 2.5 m)。考慮云霧區(qū)尺寸及試驗布場中壓力測點位置(最遠點50 m)，建立半徑R=55 m、高度H=20 m的圓柱形計算域。

根據(jù)對稱性，采用1/2模型進行計算與六面體網(wǎng)格劃分，如圖2所示。

圖2 85 kg FAE云霧爆轟數(shù)值仿真模型Fig.2 Numerical simulation model of 85 kg FAE cloud detonation

1.2 計算方法與邊界條件

考慮不同網(wǎng)格尺寸對計算精度的影響，對網(wǎng)格單元邊長為1.0 m、0.8 m、0.6 m、0.4 m這4種條件進行對比計算，取90°方向50 m測點的峰值超壓Δp50對比結(jié)果見表1所示。網(wǎng)格尺寸的差異會導致模型同一測點計算結(jié)果的差異，隨著網(wǎng)格尺寸的細化，網(wǎng)格尺寸對計算精度的影響逐漸降低。當網(wǎng)格尺寸從0.6 m細化至0.4 m時，計算值僅改變4.2%，綜合計算精度和計算成本，選取網(wǎng)格尺寸為0.6 m. 模型六面體單元總數(shù)約50萬。

表1 不同網(wǎng)格尺寸峰值超壓對比Tab.1 Overpressure comparison of different mesh sizes

為簡化模型，假設(shè)FAE云霧符合Chapman-Jouguet(CJ)理論[8]，即云霧的化學反應(yīng)在爆轟波通過后瞬間完成并釋放反應(yīng)熱。空氣簡化為符合γ定律的理想氣體[18],初始壓強為1×105Pa. 模型采用任意拉格朗日- 歐拉(ALE)多物質(zhì)算法進行計算求解。云霧起爆時刻定義為0 ms，計算總時長為150 ms，數(shù)據(jù)保存時間間隔為0.5 ms.

計算域地面簡化為剛性壁面，能夠全反射沖擊波，設(shè)為剛性壁面邊界條件。模型Oyz平面為對稱平面，對該平面節(jié)點x軸方向運動施加約束。該模型空氣區(qū)邊界采用無反射邊界條件，并在邊界單元沿內(nèi)法向施加1×105Pa的壓力，以模擬外界大氣對計算區(qū)域的壓力。

1.3 材料模型與參數(shù)

材料包括空氣區(qū)材料與云霧區(qū)材料兩部分。

空氣采用MAT_NULL材料模型和EOS_LINEAR_POLYNOMIAL狀態(tài)方程，其詳細參數(shù)[18]如表2所示。狀態(tài)方程表達形式為

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E,

(1)

式中：p為壓力(Pa)；μ=V-1，V為比容,V=ρ/ρ0，ρ為密度(kg/m3)，ρ0為初始密度(kg/m3)；C0～C6為與材料性質(zhì)相關(guān)的常數(shù)；E為能量密度(J/m3)，即單位體積內(nèi)能。

表2 空氣區(qū)材料參數(shù)表Tab.2 Parameters of air zone

云霧區(qū)采用高能炸藥燃燒模型(MAT_HIGH_EXPLOSIVE_BURN)，模型所需參數(shù)由實驗中壓力傳感器數(shù)據(jù)計算所得，爆速D=1 270 m/s，爆壓pCJ=3.32 MPa，初始密度ρ0由(2)式估算：

(2)

式中：K為爆轟產(chǎn)物多方指數(shù)，此處近似取1.2. 由(2)式可得初始密度ρ0=4.53 kg/m3.

云霧狀態(tài)方程采用JWL狀態(tài)方程。JWL狀態(tài)方程表達形式為

(3)

式中：E0為FAE云霧的初始能量密度(J/m3)，根據(jù)單位體積內(nèi)云爆燃料所釋放的化學能可得初始能量密度E0=2×106J/m3；A、B、R1、R2、ω為與炸藥爆轟性能相關(guān)的參數(shù)，是數(shù)值模型中需要待定的參數(shù)。

2 JWL狀態(tài)方程參數(shù)尋優(yōu)計算

2.1 JWL參數(shù)相容關(guān)系

等熵條件下，JWL方程[19]為

ps=Ae-R1V+Be-R2V+CV-(ω+1),

(4)

式中：ps為等熵壓力(Pa)；C為方程參數(shù)。根據(jù)CJ理論，等熵線、波速線、Hugoniot線相切于一點，據(jù)此得到JWL參數(shù)相容方程組[19,21]為

(5)

式中：VCJ為CJ比容，VCJ=K/(K+1).A、B、C為表征壓力的參數(shù)，單位是Pa；R1、R2、ω無量綱。ρ0、D、E0、pCJ、VCJ均已計算得到，因此A、B、C、R1、R2、ω這6個參數(shù)僅有3個是獨立的，當給定一組R1、R2、ω值，A、B、C可通過(5)式解出。

R1、R2、ω的取值范圍[21]是4.0≤R1≤7.0、0.8≤R2≤2、0.2≤ω≤0.3，由(5)式解出A、B、C為正值的解。

2.2 構(gòu)造仿真計算值與實測值偏差函數(shù)

由2.1節(jié)可知，每給定一組R1、R2、ω值，就能求解其余JWL參數(shù)，進而通過云霧爆轟數(shù)值模型進行計算，得到地面各測點處峰值壓力計算值，將計算值與實驗值對比，偏差越小則此JWL參數(shù)越優(yōu)。因此構(gòu)造如(6)式偏差函數(shù)：

(6)

式中：M為測點數(shù)；qi為權(quán)重系數(shù)；pn為數(shù)值仿真得到的第i個測點峰值壓力值；pe為實測的第i個測點峰值壓力值。實驗共布設(shè)32個壓力測點(見圖1)，以爆心距離15 m為分界線，將測點分為兩種情況，每種情況有16個測點。當爆心距離小于等于15 m時，測點位于云霧覆蓋區(qū)內(nèi)部或云霧邊界處，測點壓力一般為兆帕量級；當爆心距離大于15 m時，測點位于云霧覆蓋區(qū)外部，超壓迅速下降至10-1兆帕量級甚至更低。因此將位于5 m、8 m、10 m、15 m的測點其權(quán)重系數(shù)qi設(shè)為1，將位于20 m、30 m、40 m、50 m的測點其權(quán)重系數(shù)qi設(shè)為10，通過權(quán)重系數(shù)將云霧區(qū)內(nèi)外的誤差保持在同一數(shù)量級。

至此，F(xiàn)AE爆轟產(chǎn)物的JWL參數(shù)計算轉(zhuǎn)化為一個數(shù)學上的多元優(yōu)化問題，即在4.0≤R1≤7.0、0.8≤R2≤2.0、0.2≤ω≤0.3范圍內(nèi)尋找一組R1、R2、ω值，使得偏差函數(shù)Err達到最小值。

2.3 BP神經(jīng)網(wǎng)絡(luò)聯(lián)合遺傳算法尋優(yōu)

偏差函數(shù)Err為三元函數(shù)，且形式未知，僅能通過一組R1、R2、ω值，經(jīng)過數(shù)值仿真后得到峰值壓力仿真值pn，進而算得Err值。智能算法[22]在解決此類復雜、多元、非線性函數(shù)的優(yōu)化問題有著出色的表現(xiàn)，本文便采用BP神經(jīng)網(wǎng)絡(luò)根據(jù)有限組輸入- 輸出數(shù)據(jù)訓練神經(jīng)網(wǎng)絡(luò)模型，訓練完成的BP神經(jīng)網(wǎng)絡(luò)便能夠建立R1、R2、ω值(輸入)與Err值(輸出)之間的映射關(guān)系，完成函數(shù)擬合，進一步地，采用遺傳算法對訓練好的BP神經(jīng)網(wǎng)絡(luò)模型尋找最小值，最小值所對應(yīng)的R1、R2、ω值即為本文所求參數(shù)。

2.3.1 采用BP神經(jīng)網(wǎng)絡(luò)擬合偏差函數(shù)Err

BP神經(jīng)網(wǎng)絡(luò)是一種3層或3層以上的前饋網(wǎng)絡(luò)，它由輸入層、隱含層、輸出層組成，該神經(jīng)網(wǎng)絡(luò)需要一定數(shù)量的輸入- 輸出數(shù)據(jù)來訓練，當輸出層結(jié)果不滿足期望輸出時，則通過BP不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，以逼近訓練所用的輸出數(shù)據(jù)[22]。利用此特性，本文選取36組R1、R2、ω值(作為輸入層)，根據(jù)JWL參數(shù)相容關(guān)系計算出其余3個參數(shù)A、B、C，再通過云霧爆轟數(shù)值模型計算相應(yīng)的Err值(作為輸出層)，用于訓練BP神經(jīng)網(wǎng)絡(luò)。其中R1取4.0、5.0、6.0、7.0；R2取1.1、1.4、1.7；ω取0.2、0.3、0.4. 上述取值正交化成36組數(shù)據(jù)。本文所用BP神經(jīng)網(wǎng)絡(luò)拓撲結(jié)構(gòu)如圖3所示，采用3-7-1的3層結(jié)構(gòu)，其中隱含層節(jié)點個數(shù)由2N+1規(guī)則[23]確定，N為輸入層節(jié)點數(shù)，因此隱含層為7節(jié)點。

圖3 BP神經(jīng)網(wǎng)絡(luò)拓撲結(jié)構(gòu)Fig.3 Topology of back propagation neural network

采用數(shù)學分析軟件MATLAB神經(jīng)網(wǎng)絡(luò)工具箱實現(xiàn)上述BP神經(jīng)網(wǎng)絡(luò)，算法及參數(shù)選擇參考自文獻[22]，簡述如下：隱含層和輸出層的傳遞函數(shù)為tan-Sigmoid，訓練算法為Levenberg-Marquardt算法，最大訓練次數(shù)為100，學習率為0.1，訓練目標精度為0.000 000 4. 從上述36組數(shù)據(jù)中隨機選擇32組用于訓練，其余4組用于預測輸出，測試該網(wǎng)絡(luò)的擬合性能。

由于BP神經(jīng)網(wǎng)絡(luò)選取訓練組數(shù)據(jù)的隨機性，可能導致每次計算所得的JWL狀態(tài)方程參數(shù)有細微波動，因此本文經(jīng)過多次重復訓練，挑選其中預測輸出與期望輸出最為接近的網(wǎng)絡(luò)模型如圖4所示。測試結(jié)果表明，此神經(jīng)網(wǎng)絡(luò)能夠預測Err函數(shù)且綜合誤差為0.266 0.

圖4 BP神經(jīng)網(wǎng)絡(luò)測試結(jié)果Fig.4 Test result of back propagation neural network

2.3.2 遺傳算法尋優(yōu)

遺傳算法是一類借鑒生物進化規(guī)律，即優(yōu)勝劣汰遺傳機制演化而來的搜索算法[22]，它采用概率化的尋優(yōu)方法，不需要特定規(guī)則便能自適應(yīng)地調(diào)整搜索方向，具有很好的全局尋優(yōu)能力，擅長解決復雜多元函數(shù)的極值尋優(yōu)問題。將2.3.1節(jié)中訓練好的BP神經(jīng)網(wǎng)絡(luò)作為待尋優(yōu)函數(shù)，調(diào)用MATLAB軟件優(yōu)化工具箱中的遺傳算法尋找其最小值，其參數(shù)設(shè)置[21-22]簡述如下：變量數(shù)為3(R1、R2、ω)，變量邊界為[4.0, 0.8, 0.2]至[7.0, 2.0, 0.4]，種群類型為實數(shù)向量，種群數(shù)量為20，交叉概率為0.8，變異概率為0.2，最大進化代數(shù)為300.

運行遺傳算法進化50代左右得到誤差函數(shù)Err最小值的近似解為20.76，如圖5所示，此時對應(yīng)的JWL參數(shù)取值為R1=4.0、R2=1.774、ω=0.227. 再通過(5)式解得A=4.56 MPa、B=5.89 MPa、C=0.27 MPa.

圖5 遺傳算法尋優(yōu)結(jié)果Fig.5 Optimized result of genetic algorithm

3 JWL狀態(tài)方程參數(shù)對比驗證

3.1 單爆源云霧爆轟實驗驗證

單爆源云霧爆轟實驗所用裝藥結(jié)構(gòu)橫截面亦為扇形，裝填有300 kg固液混合云爆燃料，燃料配比和裝填密度與1.1節(jié)的扇形裝藥結(jié)構(gòu)完全一致。從正上方拍攝的燃料分散形成云霧圖像如圖6(a)所示，為方便建模，將云霧截面近似成如白色虛線所示的不規(guī)則多邊形，并將裝藥結(jié)構(gòu)所在位置設(shè)為原點，經(jīng)過圖像處理計算，可得云霧邊界形狀特征坐標依次為(10.73 m, 0.38 m)、(18.49 m, 6.56 m)、(16.14 m, 10.63 m)、(6.18 m, 7.86 m)、(1.44 m, 8.53 m)、(3.59 m, 22.99 m)、(0.19 m, 22.99 m)、(-2.63 m, 9.20 m)、(-7.95 m, 7.19 m)、(-16.05 m, 15.52 m)、(-19.35 m, 13.22 m)、(-9.10 m, 2.68 m)、(-13.99 m, 1.96 m)、(-4.84 m, 16.00 m)、(3.59 m, -15.95 m)。云霧高度為7.44 m，起爆點高度為4.50 m，忽略起爆點與裝藥結(jié)構(gòu)間的水平距離。

圖6 單爆源云霧爆轟試驗與數(shù)值仿真模型Fig.6 Experimental and numerical models of single-source cloud detonation

建立300 kg異形橫截面云霧爆轟有限元模型如圖6(b)所示，其計算域半徑R=70 m，高度H=20 m，并采用2.3節(jié)所計算的該固液混合FAE云霧爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)用于該模型進行數(shù)值仿真計算。仿真壓力云圖與實驗現(xiàn)場高速攝影(布設(shè)于90°方向)結(jié)果對比如圖7所示，多相云霧被起爆約20 ms左右，云霧已全部完成爆轟，爆轟波從云霧上邊界向外傳播形成明顯的空氣沖擊波波陣面(黑色虛線為沖擊波陣面輪廓示意線)，由于云霧的特殊形狀，從云霧邊界透射入空氣的沖擊波陣面左側(cè)略高于右側(cè)，數(shù)值仿真結(jié)果也反映出了這一特征。

圖7 起爆后20 ms實驗與仿真結(jié)果對比Fig.7 Comparison of experimental and numerical results at 20 ms after ignition

選取如圖6所示沿60°、120°、180° 3個方向布設(shè)在云霧爆轟區(qū)外的地面壓力傳感器所測量的峰值超壓值，與數(shù)值計算值的對比如圖8所示。由于云霧橫截面的不規(guī)則形狀，沿著3個方向的云霧覆蓋范圍各不相同，導致了沿3個方向上的超壓衰減規(guī)律有所差異，在120°方向上云霧覆蓋范圍較小，超壓衰減較快；而60°方向和180°方向云霧覆蓋范圍較大，超壓衰減稍慢。隨著距爆心距離的增加，計算值與實測值逐漸接近，在距爆心距離為50 m時，沿3個方向計算值與4發(fā)實驗實測值的平均值之間相對誤差分別為8.2%、9.0%、5.5%. 這表明所提出的FAE爆轟產(chǎn)物JWL方程參數(shù)計算方法與計算所得的參數(shù)能夠模擬外場單爆源FAE爆轟。

圖8 數(shù)值仿真與實測地面超壓結(jié)果對比Fig.8 Comparison of numerically simulated and measured results of ground peak overpressure

3.2 多爆源云霧爆轟試驗驗證

多爆源云霧爆轟試驗數(shù)據(jù)取自文獻[18]。4個扇形裝藥結(jié)構(gòu)同步分散所產(chǎn)生的云霧如圖9(a)所示， 每個裝藥結(jié)構(gòu)裝填有85 kg的固液混合云爆燃料，燃料配比、裝填密度與殼體結(jié)構(gòu)與1.1節(jié)所述裝藥結(jié)構(gòu)完全一致，裝藥結(jié)構(gòu)與試驗場地中心距離為22 m. 圖9(b)為根據(jù)試驗所建立的多爆源云霧爆轟數(shù)值仿真模型，云霧模型的尺寸與1.1節(jié)云霧尺寸相同，模型計算域半徑R=70 m，高度H=20 m.

圖9 多爆源云霧爆轟實驗與數(shù)值仿真模型Fig.9 Experimental and numerical models of multi-source cloud detonation

縱截面壓力云圖與試驗圖像結(jié)果對比如圖10所示。從圖10(b)中可以看出，多爆源云霧同步爆轟產(chǎn)生的沖擊波于模型中心處相互作用，產(chǎn)生局部高壓區(qū)域。圖11為沿兩兩云霧之間布設(shè)的地面壓力傳感器(沿45°方向)采集的峰值超壓數(shù)據(jù)與數(shù)值計算值之間的對比結(jié)果，計算值與實測值在變化趨勢上一致，均反映了沖擊波疊加后出現(xiàn)高壓區(qū)(15 m位置)再隨傳播而衰減的趨勢。在距離模型中心較近區(qū)域，計算值與實測值沒有完全吻合，這是由于在實驗條件下，燃料由爆炸拋撒形成云霧，每個云霧爆源的形狀尺寸及濃度分布都存在著差異，而數(shù)值計算模型則無法反映這些差異。在模型中心處(0 m位置)，計算結(jié)果反映出了沖擊波疊加使得附近的壓力顯著增加。隨著沖擊波傳播距離的增加，各云霧爆源差異對峰值超壓影響逐漸下降，計算值與實測值逐漸接近，在50 m測點處，計算值與實測值偏差為11.1%. 因此本文所提出的FAE爆轟產(chǎn)物JWL方程參數(shù)計算方法與計算所得的參數(shù)也能夠模擬外場多爆源FAE爆轟。

圖10 起爆后24 ms FAE云霧爆轟試驗與仿真結(jié)果對比Fig.10 Comparison of measured and numerically simulated results of FAE cloud detonation at 24 ms after ignition

圖11 45°方向數(shù)值仿真與實測地面超壓結(jié)果對比Fig.11 Comparison of numerically simulated and measured results of ground peak overpressure in 45° direction

4 結(jié)論

本文類比標準圓筒試驗的JWL狀態(tài)方程參數(shù)計算循環(huán)過程，使用外場云爆試驗的峰值壓力數(shù)據(jù)建立FAE爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)的確定方法。得出主要結(jié)論如下：

1)引入BP神經(jīng)網(wǎng)絡(luò)聯(lián)合遺傳算法簡化狀態(tài)方程參數(shù)優(yōu)化過程，增加了尋優(yōu)速度和精度。計算得到了固液混合FAE爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)，可模擬單爆源及多爆源云霧爆轟與沖擊波傳播過程。

2)對于單爆源云霧爆轟，計算所得的壓力云圖與云霧起爆后所拍攝的圖像在形貌上一致，在距爆心距離50 m測點處，峰值超壓的計算值與實測值偏差分別為8.2%、9.0%、5.5%.

3)對于多爆源云霧爆轟，沿45°方向的峰值超壓仿真值與試驗值在變化趨勢上一致，在50 m測點處，計算值與實測值偏差為11.1%.