JWL狀態(tài)方程及其等效多方狀態(tài)方程在內(nèi)爆炸計(jì)算中的應(yīng)用分析

徐維錚，吳衛(wèi)國

1武漢理工大學(xué)高性能艦船技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北武漢430063

2武漢理工大學(xué)交通學(xué)院，湖北武漢430063

0 引 言

從戰(zhàn)例統(tǒng)計(jì)和試驗(yàn)資料來看，反艦導(dǎo)彈已成為艦船生命力的主要威脅，艦船在艙內(nèi)爆炸沖擊載荷作用下的毀傷效應(yīng)分析及防護(hù)結(jié)構(gòu)設(shè)計(jì)是目前艦船工程領(lǐng)域的研究熱點(diǎn)。對(duì)艙內(nèi)爆炸沖擊載荷的準(zhǔn)確分析是結(jié)構(gòu)響應(yīng)評(píng)估和防護(hù)結(jié)構(gòu)設(shè)計(jì)載荷輸入的前提，因此起著至關(guān)重要的作用。數(shù)值模擬作為一種研究內(nèi)爆炸載荷的有效方式，得到了國內(nèi)外學(xué)者的普遍重視。

針對(duì)約束空間內(nèi)爆炸過程的數(shù)值模擬，國內(nèi)外學(xué)者開展了許多研究。曹玉忠等［1］采用瞬時(shí)爆轟后的高溫高壓氣體作為爆炸初場(chǎng)，其與周圍空氣均采用理想氣體狀態(tài)方程來描述，數(shù)值模擬研究了半球頂圓柱筒密閉式抗爆容器的內(nèi)部爆炸流場(chǎng)。Benselama等［2］基于瞬時(shí)爆轟假定，采用理想氣體狀態(tài)方程或JWL狀態(tài)方程（Jones-Wilkins-Lee Equation of State，JWL-EOS）描述爆轟產(chǎn)物，對(duì)比分析了數(shù)值模擬結(jié)果與空中爆炸（空爆）的超壓實(shí)測(cè)結(jié)果，結(jié)果表明，爆轟產(chǎn)物的狀態(tài)方程對(duì)于空爆超壓載荷的影響較小。任會(huì)蘭等［3］采用基于消息傳遞接口（MPI）模式開發(fā)的內(nèi)爆炸過程三維并行程序，對(duì)不同藥量的炸藥在典型工事中的爆炸過程進(jìn)行了數(shù)值研究。在模擬過程中，假定炸藥瞬時(shí)爆轟，其爆轟產(chǎn)物采用可變指數(shù)多方氣體狀態(tài)方程進(jìn)行描述。Ma等［4］基于瞬時(shí)爆轟假定，采用可變指數(shù)多方氣體狀態(tài)方程來描述爆轟產(chǎn)物，通過求解非守恒形式的可壓縮歐拉方程，數(shù)值模擬研究了爆炸波在地鐵站內(nèi)的傳播過程。Sugiyama等［5］基于瞬時(shí)爆轟假定，采用理想的氣體狀態(tài)方程描述了爆轟產(chǎn)物及周圍空氣，全場(chǎng)具有統(tǒng)一的氣體絕熱指數(shù)1.4，在研究中，對(duì)二維軸對(duì)稱可壓縮歐拉方程進(jìn)行數(shù)值離散，研究了房間內(nèi)爆炸泄出過程中超壓峰值隨傳播距離衰減的特性。近期，徐維錚等［6-8］基于瞬時(shí)爆轟假定，同樣采用理想氣體狀態(tài)方程描述爆轟產(chǎn)物及周圍空氣，使用自主開發(fā)的三維內(nèi)爆炸波高精度數(shù)值計(jì)算程序，探討了泄壓口大小、WENO格式精度、炸藥質(zhì)量等因素對(duì)約束空間內(nèi)爆炸準(zhǔn)靜態(tài)超壓載荷的影響規(guī)律。

從上述研究現(xiàn)狀可知，描述炸藥瞬時(shí)爆轟后的爆轟產(chǎn)物狀態(tài)的方程主要有3種：JWL狀態(tài)方程、可變指數(shù)多方氣體狀態(tài)方程和理想氣體狀態(tài)方程。其中，JWL狀態(tài)方程既可用來描述爆炸沖擊載荷的高壓段，也可描述低壓段，并且已被目前的商用程序（如AUTODYN，LSDYNA等）所采用；可變指數(shù)多方氣體狀態(tài)方程是以一種隨密度改變絕熱指數(shù)的方式來近似描述爆轟產(chǎn)物狀態(tài)的方程，目前已被相關(guān)學(xué)者采用［3-4］；雖然上述研究并未給出確定其中參數(shù)的方法，且理想氣體狀態(tài)方程在描述高壓段時(shí)也有一定的誤差，但作為一種簡單的狀態(tài)方程，它便于程序編寫，因此有學(xué)者也將其應(yīng)用到了內(nèi)爆炸場(chǎng)的數(shù)值研究中。

在真實(shí)內(nèi)爆炸過程中，針對(duì)不同炸藥的JWL狀態(tài)方程，其ω參數(shù)存在一定的差異，并對(duì)爆炸過程會(huì)產(chǎn)生一定的影響。當(dāng)爆炸發(fā)生在不同氣體環(huán)境時(shí)，初始?xì)怏w絕熱指數(shù)γ0的差異也會(huì)對(duì)爆炸過程產(chǎn)生影響。因此，可變指數(shù)多方氣體狀態(tài)方程作為描述爆轟產(chǎn)物狀態(tài)的一種等效方式，具備理想氣體狀態(tài)方程的簡潔性，且便于程序編寫及數(shù)值計(jì)算，是一種較有優(yōu)勢(shì)的狀態(tài)方程。

鑒此，本文將基于Fortran平臺(tái)，采用五階WENO有限差分格式，利用自主開發(fā)的艙室內(nèi)炸藥爆炸過程二維數(shù)值計(jì)算程序，研究參數(shù)ω及γ0在不同取值下對(duì)內(nèi)爆炸參數(shù)計(jì)算的影響規(guī)律，以此提出用于確定等效多方氣體狀態(tài)方程中參數(shù)的方法。

1 控制方程及數(shù)值求解

本文將整個(gè)爆炸場(chǎng)視作由爆轟產(chǎn)物及空氣兩種成分組成，將爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)方程與可壓縮歐拉方程進(jìn)行耦合求解，具體形式如下［9］：

其中，

式中：U為守恒變量；G，H分別為x,y方向的數(shù)值通量；ρ為氣體密度；u,v分別為x,y方向氣體的速度分量；p為氣體壓力；E為單位體積流體的總能量；e為氣體內(nèi)能；λ為爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)。

由于式（1）中含有ρ，u，v，p，e，λ這6個(gè)未知變量，所以還需補(bǔ)充狀態(tài)方程才能將方程組封閉。在本文計(jì)算中，爆轟產(chǎn)物的狀態(tài)采用JWL狀態(tài)方程［2］描述：

式中：A，B，R1，R2，ω均為JWL狀態(tài)方程中的參數(shù)，其中A=3.712×1011Pa，B=3.232×109Pa，R1=4.15，R2=0.95，ω=0.3，ρE為炸藥密度，ρE=1 630 kg/m3。

此外，周圍氣體采用理想氣體狀態(tài)方程描述：

式中：γ0=1.4。此時(shí)初始?xì)怏w密度ρ0=1.0 kg/m3；初始大氣壓力p0=1.0×105Pa。

本文分別采用五階WENO有限差分格式［10］和三階TVD龍格—庫塔（TVD-RK）法對(duì)上述方程中的空間項(xiàng)及時(shí)間項(xiàng)進(jìn)行數(shù)值離散［7］。

2 內(nèi)爆炸數(shù)值計(jì)算

2.1 艙室尺寸及測(cè)點(diǎn)布置

本節(jié)主要討論狀態(tài)方程（4）及式（5）中ω和γ0參數(shù)的不同取值對(duì)內(nèi)爆炸參數(shù)計(jì)算的影響規(guī)律。圖1（a）所示為用于計(jì)算的長方形艙室尺寸，其中艙內(nèi)設(shè)置有3個(gè)測(cè)點(diǎn)，通過測(cè)點(diǎn)輸出爆炸參數(shù)的時(shí)程曲線。艙室壁面邊界條件設(shè)為剛性邊界，不考慮艙室結(jié)構(gòu)的變形［11］。均勻網(wǎng)格尺寸取為10 mm，如圖1（b）所示。

圖1 艙室尺寸及網(wǎng)格劃分Fig.1 Sizes of cabin and mesh partition

2.2 爆炸初場(chǎng)

基于瞬時(shí)爆轟假定［5］，本文將初始方形炸藥（邊長l=40 mm）等效為均勻的高壓氣團(tuán)，如圖2所示。具體參數(shù)按照J(rèn)WL狀態(tài)方程計(jì)算如下：，其中氣體初始內(nèi)能e0等于炸藥爆熱值Q，e0=4.69×106J/kg。炸藥周圍氣體參數(shù)由理想氣體狀態(tài)方程計(jì)算，即ρ=ρ0,p=p0。

圖2 爆炸初場(chǎng)（ω=0.3,γ0=1.4）Fig.2 Initial condition for the explosion（ω =0.3,γ0=1.4）

2.3 參數(shù)ω的影響

由于炸藥類型不同，ω的值也會(huì)有所不同，所以探討ω的取值對(duì)爆炸過程計(jì)算的影響規(guī)律具有一定的工程意義。本文將式（4）右端前兩個(gè)指數(shù)高壓項(xiàng)簡記為A+B項(xiàng)，其隨氣體密度的變化曲線如圖3所示。

圖3 指數(shù)高壓項(xiàng)隨氣體密度衰減曲線Fig.3 Curves of the exponential high pressure term attenuated with gas density

由圖3可知，隨著爆轟產(chǎn)物密度的降低，高壓項(xiàng)數(shù)值迅速減小，當(dāng)爆轟產(chǎn)物密度膨脹為初始密度的1/12時(shí)，高壓項(xiàng)數(shù)值趨于0。即爆轟產(chǎn)物膨脹到一定程度后，爆轟產(chǎn)物JWL狀態(tài)方程演化為理想氣體狀態(tài)方程形式：

且等效氣體絕熱指數(shù)γe=ω+1。

由式（4）及式（6）可知，在 JWL狀態(tài)方程中A，B，R1，R2，ρE這5個(gè)參數(shù)存在于前兩個(gè)指數(shù)高壓項(xiàng)中，因在爆炸過程中高壓項(xiàng)迅速衰減，故本文不對(duì)這5個(gè)參數(shù)的影響進(jìn)行探討。鑒于參數(shù)ω既可影響爆炸沖擊載荷高壓段的計(jì)算，也可影響其低壓段的計(jì)算，且不同炸藥類型有不同的ω值［12］，因此，為探討參數(shù)ω對(duì)整個(gè)爆炸過程計(jì)算的影響規(guī)律，本文算例選取了3種ω值（ω=0.2，0.3，0.4）進(jìn)行計(jì)算，同時(shí)保持γ0=1.4 。

圖4所示為不同ω取值下炸藥中心測(cè)點(diǎn)3處的超壓衰減時(shí)程曲線。由圖可知，隨著ω值的增大，炸藥內(nèi)部的壓力衰減速率減緩。

圖5所示為不同ω取值下測(cè)點(diǎn)1和測(cè)點(diǎn)2處爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)λ的時(shí)程曲線。由圖可知，ω的取值對(duì)于不同測(cè)點(diǎn)處的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)時(shí)程曲線有較大影響，且隨著ω值的增大，爆轟產(chǎn)物到達(dá)艙壁面的初始時(shí)間提前，亦即ω值的增大加快了爆炸初期爆轟產(chǎn)物的膨脹速率。

圖5 不同ω取值下測(cè)點(diǎn)處的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)時(shí)程曲線Fig.5 Time history curves of mass fraction of detonation products at gauging points for different values ofω

為進(jìn)一步分析爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)的最終演變趨勢(shì)，通過質(zhì)量加權(quán)，可得到充分混合后的艙內(nèi)爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)λf為

式中：m為炸藥質(zhì)量；mkong為空氣質(zhì)量；V為艙室體積。爆炸初期，爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)為1.0，周圍環(huán)境氣體對(duì)應(yīng)的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)為0。

由式（7）可知，λf僅與炸藥質(zhì)量、艙室體積相關(guān)，而與ω的取值無關(guān)，即艙內(nèi)爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)將最終趨于一致。從圖5也可見這種趨勢(shì)。

圖6所示為不同ω取值下測(cè)點(diǎn)1和測(cè)點(diǎn)2處氣體絕熱指數(shù)γ的時(shí)程曲線。由圖可知，ω的取值對(duì)氣體絕熱指數(shù)的時(shí)程有著顯著影響，即隨著ω的增大，最終混合物的氣體絕熱指數(shù)也增大。

圖6 不同ω取值下測(cè)點(diǎn)處的氣體絕熱指數(shù)時(shí)程曲線Fig.6 Time history curves of gas specific heat ratio at gauging points for different values ofω

根據(jù)式（7）的分析思路，可得到艙內(nèi)最終混合物的氣體絕熱指數(shù)γf為

式（8）的推導(dǎo)用到了式（6）的結(jié)論。由式（6）可知，當(dāng)爆轟產(chǎn)物膨脹到后期，JWL狀態(tài)方程演化為理想的氣體狀態(tài)方程形式，此時(shí)爆轟產(chǎn)物的絕熱指數(shù)為γe，艙室內(nèi)最終混合物的氣體絕熱指數(shù)γf為爆轟產(chǎn)物與空氣絕熱指數(shù)的質(zhì)量加權(quán)之和。因此，由式（8）可知，當(dāng)其他參數(shù)不變時(shí)，γf會(huì)隨參數(shù)ω的增大而增大。

圖7所示為不同ω取值下測(cè)點(diǎn)處的超壓時(shí)程曲線。由圖可知，隨著ω的增大，沖擊波峰值增大，沖擊波到達(dá)艙壁面的初始時(shí)間提前，且最終形成的準(zhǔn)靜態(tài)壓力峰值也隨之增大。

圖7 不同ω取值下測(cè)點(diǎn)處的超壓時(shí)程曲線Fig.7 Time history curvs of overpressure at gauging points for different values ofω

為了說明準(zhǔn)靜態(tài)壓力峰值增大的原因，本文采用了文獻(xiàn)［6］中的密閉空間內(nèi)爆炸準(zhǔn)靜態(tài)超壓峰值Δps計(jì)算公式：

然后，進(jìn)一步將式（8）代入式（9），可得

由式（10）可知，當(dāng)其他參數(shù)不變時(shí)，Δps將隨著參數(shù)ω的增大而增大。

2.4 參數(shù) γ0的影響

當(dāng)爆炸發(fā)生在不同氣體環(huán)境（例如氮?dú)?、氧氣、空氣）中時(shí)，不同氣體有著不同的初始絕熱指數(shù)γ0［13］，并對(duì)爆炸過程產(chǎn)生一定的影響，故探討γ0的取值對(duì)于研究其對(duì)爆炸參數(shù)計(jì)算的影響規(guī)律具有一定的工程意義。本文選取3種γ0值（γ0=1.2，1.3，1.4），同時(shí)保持ω=0.3。

圖8所示為不同γ0取值下測(cè)點(diǎn)處的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)時(shí)程曲線。由圖可知，γ0的取值對(duì)不同測(cè)點(diǎn)處爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)時(shí)程曲線有著較大影響，且隨著γ0值的減小，爆轟產(chǎn)物到達(dá)艙壁面的初始時(shí)間提前，即γ0值的增大可使爆炸初期爆轟產(chǎn)物的膨脹速率減緩。

圖8 不同γ0取值下測(cè)點(diǎn)處的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)時(shí)程曲線Fig.8 Time history curves of mass fraction of detonation products at gauging points for different values ofγ0

圖9所示為不同γ0取值下測(cè)點(diǎn)處的氣體絕熱指數(shù)時(shí)程曲線。由圖可知，γ0的取值對(duì)氣體絕熱指數(shù)γ的時(shí)程有著顯著影響，隨著γ0的增大，最終混合物的氣體絕熱指數(shù)γf增大。由式（8）也可知，當(dāng)其他參數(shù)不變時(shí)，γf隨著參數(shù)γ0的增大而增大。

圖9 不同γ0取值下測(cè)點(diǎn)處的氣體絕熱指數(shù)時(shí)程曲線Fig.9 Time history curves of gas specific heat ratio at gauging points for different values ofγ0

圖10所示為不同γ0取值下測(cè)點(diǎn)處的超壓時(shí)程曲線。由圖可知，隨著γ0的增大，沖擊波峰值增大，沖擊波到達(dá)艙室壁面的初始時(shí)間提前。由圖還可看出，γ0的取值對(duì)最終形成的準(zhǔn)靜態(tài)壓力峰值的影響很小，具體可由式（10）判斷。由式（10）可知，當(dāng)其他參數(shù)不變時(shí)，Δps與γ0呈非線性的關(guān)系，乘號(hào)兩端均包含γ0且具有相反的變化特征，從而使得準(zhǔn)靜態(tài)壓力峰值的變化并不顯著。

3 等效多方氣體狀態(tài)方程

圖10 不同γ0取值下測(cè)點(diǎn)處的超壓時(shí)程曲線Fig.10 Time history curves of overpressure at gauging points for different values ofγ0

由于JWL狀態(tài)方程中的參數(shù)較多且包含指數(shù)的計(jì)算，因此使程序變量增多、計(jì)算時(shí)間增加。為提高工程計(jì)算的效率，作為一種可靠的方式，可以采用一種等效的多方氣體狀態(tài)方程來描述爆轟產(chǎn)物的狀態(tài)。例如，文獻(xiàn)［3-4］在研究炸藥爆炸數(shù)值計(jì)算時(shí)，采用了可變絕熱指數(shù)多方氣體狀態(tài)方程來描述炸藥爆轟產(chǎn)物，具體形式如下：

式中：γb0為爆轟產(chǎn)物的初始等效絕熱指數(shù)；γbf為爆轟產(chǎn)物膨脹到最后的絕熱指數(shù)；下標(biāo)b為爆轟產(chǎn)物絕熱指數(shù)隨膨脹過程的衰減系數(shù)。

3.1 參數(shù)取值

盡管文獻(xiàn)［3-4］提出了采用可變絕熱指數(shù)多方氣體狀態(tài)方程來描述炸藥爆轟產(chǎn)物，但卻未對(duì)炸藥爆炸過程影響較大的參數(shù)γb0，γbf及b的取值問題進(jìn)行充分探討。實(shí)際上，如何合理地確定這3個(gè)參數(shù)值十分重要。因此，本文擬基于JWL狀態(tài)方程來探討上述3個(gè)參數(shù)的取值方法，提出一種與JWL狀態(tài)方程近似等效的多方氣體狀態(tài)方程，并通過相關(guān)數(shù)值算例來驗(yàn)證其可行性及可靠性。

本文研究時(shí)，γb0的取值由式（12）計(jì)算。假設(shè)炸藥瞬時(shí)爆轟，則ρ=ρE，e=Q，將其代入式（4）和式（11）中，按照壓力等效原則，可得

式中，γbf的取值由式（13）計(jì)算。當(dāng)爆轟產(chǎn)物得到充分膨脹后，產(chǎn)物的狀態(tài)方程退化為式（6），并按照壓力等效原則，可得

式中，b的取值直接影響到氣體絕熱指數(shù)隨密度變化的關(guān)系（主要在高壓段起作用），目前，暫無一種解析方式可直接確定其數(shù)值，但可通過數(shù)值計(jì)算結(jié)果與JWL狀態(tài)方程的計(jì)算結(jié)果進(jìn)行對(duì)比，來近似確定其數(shù)值。

3.2 數(shù)值計(jì)算結(jié)果分析

本文選取2.2節(jié)的炸藥參數(shù)（邊長l=40 mm）進(jìn)行數(shù)值計(jì)算。圖11所示為不同b取值下測(cè)點(diǎn)處的超壓時(shí)程曲線與JWL狀態(tài)方程計(jì)算結(jié)果的對(duì)比。由圖可知，當(dāng)參數(shù)b取為2 000時(shí)，數(shù)值計(jì)算得到的超壓載荷與JWL狀態(tài)方程的計(jì)算結(jié)果吻合較好。

圖11 不同b取值下測(cè)點(diǎn)處的超壓時(shí)程曲線（l=40 mm）Fig.11 Time history curves of overpressure at gauging points for different values ofb（l=40 mm）

圖12所示為不同b取值下測(cè)點(diǎn)處爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)時(shí)程曲線，圖13所示為不同b取值下測(cè)點(diǎn)處氣體絕熱指數(shù)時(shí)程曲線（兩圖中方形炸藥邊長l=40 mm）。由圖12和13圖可知，b的取值對(duì)爆炸參數(shù)時(shí)程曲線影響較大，當(dāng)b=2 000時(shí)，其在高壓段（爆炸初期）與JWL狀態(tài)方程的計(jì)算結(jié)果更吻合。因此，針對(duì)該爆炸工況，b=2 000是一種合理的取值。

圖12 不同b取值下測(cè)點(diǎn)處爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)時(shí)程曲線Fig.12 Time history curves of mass fraction of detonation products at gauging points for different values ofb

圖13 不同b取值下測(cè)點(diǎn)處的氣體絕熱指數(shù)時(shí)程曲線Fig.13 Time history curves of gas specific heat ratio at gauging points for different values ofb

為了驗(yàn)證參數(shù)b=2 000同樣適用于其他藥量，本文在其他參數(shù)不變的情況下，選取了邊長l=160 mm的方形炸藥再次進(jìn)行數(shù)值計(jì)算。圖14給出了在此炸藥參數(shù)及不同b取值下兩個(gè)測(cè)點(diǎn)處的超壓時(shí)程曲線，圖15給出了不同b取值下測(cè)點(diǎn)2處的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)時(shí)程曲線及氣體絕熱指數(shù)時(shí)程曲線。由對(duì)圖14和圖15的分析可知，針對(duì)大質(zhì)量炸藥內(nèi)爆炸工況，參數(shù)b=2 000的數(shù)值計(jì)算結(jié)果與JWL狀態(tài)方程的計(jì)算結(jié)果同樣吻合較好。

圖14 不同b取值下測(cè)點(diǎn)處的超壓時(shí)程曲線（l=160 mm）Fig.14 Time history curves of overpressure at gauging points for different values ofb（l=160 mm）

圖15 不同b取值下測(cè)點(diǎn)2處的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)及氣體絕熱指數(shù)時(shí)程曲線Fig.15 Time history cures of mass fraction of detonation products and gas specific heat ratio at gauging point-2 for different values ofb

4 結(jié) 論

通過本文的研究，得到如下主要結(jié)論：

1）在JWL狀態(tài)方程中，參數(shù)ω值的增大可使爆炸初期爆轟產(chǎn)物的膨脹速率加快，最終增大了混合物氣體絕熱指數(shù)和沖擊波峰值，使沖擊波到達(dá)艙室壁面的初始時(shí)間提前，且最終形成的準(zhǔn)靜態(tài)壓力峰值也增大。

2）初始?xì)怏w絕熱指數(shù)γ0的增大可使爆炸初期爆轟產(chǎn)物的膨脹速率減緩，增大了最終混合物氣體絕熱指數(shù)和沖擊波峰值，且沖擊波到達(dá)艙室壁面的初始時(shí)間提前，但對(duì)最終形成的準(zhǔn)靜態(tài)壓力峰值影響很小。

3）本文用來確定等效多方氣體狀態(tài)方程參數(shù)的方法具有可行性和可靠性，通過數(shù)值模擬結(jié)果分析，初步確定了其中的參數(shù)b=2 000。