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

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 100081)

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

1.2 計(jì)算方法與邊界條件

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

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

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

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

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

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

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

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í)驗(yàn)中壓力傳感器數(shù)據(jù)計(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)方程表達(dá)形式為

(3)

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

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

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

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

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

(4)

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

(5)

式中：VCJ為CJ比容，VCJ=K/(K+1).A、B、C為表征壓力的參數(shù)，單位是Pa；R1、R2、ω?zé)o量綱。ρ0、D、E0、pCJ、VCJ均已計(jì)算得到，因此A、B、C、R1、R2、ω這6個(gè)參數(shù)僅有3個(gè)是獨(dú)立的，當(dāng)給定一組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)造仿真計(jì)算值與實(shí)測(cè)值偏差函數(shù)

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

(6)

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

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

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

偏差函數(shù)Err為三元函數(shù)，且形式未知，僅能通過一組R1、R2、ω值，經(jīng)過數(shù)值仿真后得到峰值壓力仿真值pn，進(jìn)而算得Err值。智能算法[22]在解決此類復(fù)雜、多元、非線性函數(shù)的優(yōu)化問題有著出色的表現(xiàn)，本文便采用BP神經(jīng)網(wǎng)絡(luò)根據(jù)有限組輸入- 輸出數(shù)據(jù)訓(xùn)練神經(jīng)網(wǎng)絡(luò)模型，訓(xùn)練完成的BP神經(jīng)網(wǎng)絡(luò)便能夠建立R1、R2、ω值(輸入)與Err值(輸出)之間的映射關(guān)系，完成函數(shù)擬合，進(jìn)一步地，采用遺傳算法對(duì)訓(xùn)練好的BP神經(jīng)網(wǎng)絡(luò)模型尋找最小值，最小值所對(duì)應(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ù)來訓(xùn)練，當(dāng)輸出層結(jié)果不滿足期望輸出時(shí)，則通過BP不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，以逼近訓(xùn)練所用的輸出數(shù)據(jù)[22]。利用此特性，本文選取36組R1、R2、ω值(作為輸入層)，根據(jù)JWL參數(shù)相容關(guān)系計(jì)算出其余3個(gè)參數(shù)A、B、C，再通過云霧爆轟數(shù)值模型計(jì)算相應(yīng)的Err值(作為輸出層)，用于訓(xùn)練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ò)拓?fù)浣Y(jié)構(gòu)如圖3所示，采用3-7-1的3層結(jié)構(gòu)，其中隱含層節(jié)點(diǎn)個(gè)數(shù)由2N+1規(guī)則[23]確定，N為輸入層節(jié)點(diǎn)數(shù)，因此隱含層為7節(jié)點(diǎn)。

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

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

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

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

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

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

運(yùn)行遺傳算法進(jìn)化50代左右得到誤差函數(shù)Err最小值的近似解為20.76，如圖5所示，此時(shí)對(duì)應(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ù)對(duì)比驗(yàn)證

3.1 單爆源云霧爆轟實(shí)驗(yàn)驗(yàn)證

單爆源云霧爆轟實(shí)驗(yàn)所用裝藥結(jié)構(gòu)橫截面亦為扇形，裝填有300 kg固液混合云爆燃料，燃料配比和裝填密度與1.1節(jié)的扇形裝藥結(jié)構(gòu)完全一致。從正上方拍攝的燃料分散形成云霧圖像如圖6(a)所示，為方便建模，將云霧截面近似成如白色虛線所示的不規(guī)則多邊形，并將裝藥結(jié)構(gòu)所在位置設(shè)為原點(diǎn)，經(jīng)過圖像處理計(jì)算，可得云霧邊界形狀特征坐標(biāo)依次為(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，起爆點(diǎn)高度為4.50 m，忽略起爆點(diǎn)與裝藥結(jié)構(gòu)間的水平距離。

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

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

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

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

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

3.2 多爆源云霧爆轟試驗(yàn)驗(yàn)證

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

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

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

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

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

4 結(jié)論

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

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

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

3)對(duì)于多爆源云霧爆轟，沿45°方向的峰值超壓仿真值與試驗(yàn)值在變化趨勢(shì)上一致，在50 m測(cè)點(diǎn)處，計(jì)算值與實(shí)測(cè)值偏差為11.1%.