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

1武漢理工大學高性能艦船技術教育部重點實驗室，湖北武漢430063

2武漢理工大學交通學院，湖北武漢430063

0 引 言

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

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

從上述研究現(xiàn)狀可知，描述炸藥瞬時爆轟后的爆轟產(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)的方程，目前已被相關學者采用［3-4］；雖然上述研究并未給出確定其中參數(shù)的方法，且理想氣體狀態(tài)方程在描述高壓段時也有一定的誤差，但作為一種簡單的狀態(tài)方程，它便于程序編寫，因此有學者也將其應用到了內(nèi)爆炸場的數(shù)值研究中。

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

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

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

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

其中，

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

由于式（1）中含有ρ，u，v，p，e，λ這6個未知變量，所以還需補充狀態(tài)方程才能將方程組封閉。在本文計算中，爆轟產(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。此時初始氣體密度ρ0=1.0 kg/m3；初始大氣壓力p0=1.0×105Pa。

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

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

2.1 艙室尺寸及測點布置

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

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

2.2 爆炸初場

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

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

2.3 參數(shù)ω的影響

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

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

由圖3可知，隨著爆轟產(chǎn)物密度的降低，高壓項數(shù)值迅速減小，當爆轟產(chǎn)物密度膨脹為初始密度的1/12時，高壓項數(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個參數(shù)存在于前兩個指數(shù)高壓項中，因在爆炸過程中高壓項迅速衰減，故本文不對這5個參數(shù)的影響進行探討。鑒于參數(shù)ω既可影響爆炸沖擊載荷高壓段的計算，也可影響其低壓段的計算，且不同炸藥類型有不同的ω值［12］，因此，為探討參數(shù)ω對整個爆炸過程計算的影響規(guī)律，本文算例選取了3種ω值（ω=0.2，0.3，0.4）進行計算，同時保持γ0=1.4 。

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

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

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

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

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

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

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

圖6 不同ω取值下測點處的氣體絕熱指數(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）的推導用到了式（6）的結論。由式（6）可知，當爆轟產(chǎn)物膨脹到后期，JWL狀態(tài)方程演化為理想的氣體狀態(tài)方程形式，此時爆轟產(chǎn)物的絕熱指數(shù)為γe，艙室內(nèi)最終混合物的氣體絕熱指數(shù)γf為爆轟產(chǎn)物與空氣絕熱指數(shù)的質(zhì)量加權之和。因此，由式（8）可知，當其他參數(shù)不變時，γf會隨參數(shù)ω的增大而增大。

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

圖7 不同ω取值下測點處的超壓時程曲線Fig.7 Time history curvs of overpressure at gauging points for different values ofω

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

然后，進一步將式（8）代入式（9），可得

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

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

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

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

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

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

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

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

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

圖10 不同γ0取值下測點處的超壓時程曲線Fig.10 Time history curves of overpressure at gauging points for different values ofγ0

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

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

3.1 參數(shù)取值

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

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

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

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

3.2 數(shù)值計算結果分析

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

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

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

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

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

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

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

圖15 不同b取值下測點2處的爆轟產(chǎn)物質(zhì)量分數(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 結 論

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

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

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

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