徐春光，白曉征，劉瑜，劉君

(1.國防科學技術大學 航天與材料工程學院，湖南 長沙410073;2.中國人民解放軍總參謀部 陸航研究所，北京101121;3.大連理工大學 航空航天學院，遼寧 大連116024)

0 引言

凝聚態(tài)炸藥爆轟及其與空氣的相互作用問題具有廣泛的應用背景，相應的實驗及理論研究已經有了上百年的歷史，例如對于常見的球形裝藥空中爆炸問題，目前已經積累了大量的實驗數(shù)據(jù)及經驗公式，理論研究也相對完善。對于結構所受爆炸波的沖擊載荷，相關研究也很豐富［1－2］。但是實驗研究具有成本高、危險性大的缺點，而數(shù)值模擬具有低成本、低危險性、高效率的優(yōu)點，而且能夠模擬無法進行實驗，或者實際中極少出現(xiàn)的故障情況，近年來隨著數(shù)值模擬技術的發(fā)展，數(shù)值模擬受到了越來越多的重視，在實驗設計、實驗結果推廣、故障診斷、工程實際等方面都得到了越來越多的應用。

數(shù)值模擬爆炸波的遠場效應時，可不考慮爆轟反應的細節(jié)而采用瞬時爆轟模型，認為炸藥爆轟瞬時完成，形成一團高溫、高壓的均勻燃氣;也可采用點爆炸模型或者經驗公式給出計算的初場，這些方法在爆炸波與物體相互作用的模擬方面都取得了很好的效果［3－4］。如果爆炸初場的半徑取得較大，則爆轟產物密度較小，可用完全氣體狀態(tài)方程來描述，這樣將進一步降低數(shù)值計算的難度。

但當模擬爆炸近區(qū)的流場時，尤其當此流場與爆轟波發(fā)展過程密切相關時，就必須同時模擬爆轟反應及其與空氣的相互作用，具體的講，必須同時考慮炸藥、爆轟產物、空氣3 種組分，以及前兩種組分間的反應動力學模型。對這種典型的多介質流動，有若干種算法可以選擇。拉格朗日型方法可以得到無耗散的介質界面，但對大變形問題，需要進行網格重分。歐拉型方法中，也有一些解決方法，如采用較為成熟的Level Set 方法捕捉界面，但同樣難以處理界面拓撲的大幅變化，而且界面處物理量不守恒，對于物理量變化劇烈的爆炸流場而言，這一問題可能是十分嚴重的;采用改進的虛擬流體方法［5］(MGFM)可很大程度上改善這一問題，但仍未根本解決。近年來有學者陸續(xù)提出了一種新的多流體計算模型［6－7］，對每種相(氣、液、固)采用獨立的連續(xù)、動量、能量方程來描述，同時求解一個附加的體積分數(shù)的方程，這一方法也在含能材料的爆轟模擬中得到了應用［8］，其主要存在的問題是求解方程多，計算量大，而且方程是非守恒的。

在某型爆炸激波管的設計中，采用TNT 裝藥的爆炸來驅動產生空氣中沖擊波，圓柱形的TNT 藥柱置于驅動管內，藥柱與驅動管內壁間為靜止空氣。在驅動管出口處引爆藥柱后，爆轟波沿藥柱向管底傳播，爆轟產物同時向側面的空氣中飛散，在管壁和管底上也引起很大的壓力載荷。這一問題中流場的發(fā)展與爆轟波的傳播過程密切相關，爆轟產物與空氣間存在復雜的相互作用，二者間界面也經歷了很大的變化。經過調研，目前的各種算法解決這一問題均存在一定程度的不足或缺陷，為此本文提出一種炸藥、爆轟產物與空氣3 種介質的混合模型，通過求解固定網格下的多組分Euler 方程組來解決這一問題。此混合模型在固相與氣相間采用等壓假設，滿足體積可加性;在氣相之間采用等溫假設。炸藥與爆轟產物均采用JWL 狀態(tài)方程，反應模型為文獻［9］提出的“點火—生長”模型，采用時空兩階精度的有限體積方法離散，采用AUSM+－up 格式計算通量。計算了空氣中球形TNT 裝藥的爆炸問題，并將爆炸近區(qū)的超壓、沖量等與實驗結果進行了比較，驗證了計算方法的有效性和精度，為爆炸激波管的模擬奠定了基礎。

1 計算方法

1.1 控制方程

考慮化學反應的三維非定常Euler 方程組為

式中:ρ 為混合氣體的密度;u、v、w 分別為x，y，z 方向速度;p 為混合氣體壓力;E 為氣體的總能，E =，e 為混合氣體的比內能，ρi為第i 種組分的密度;Qc為化學反應的放熱項;σi為第i 種組分的生成率。在本文中，下標1 代表未反應炸藥，下標2 代表產物氣體，下標3 代表空氣，于是有

1.2 狀態(tài)方程

空氣采用量熱完全氣體模型，狀態(tài)方程為

式中:T 為氣體溫度;R 為氣體常數(shù);γ 為比熱比。

炸藥及爆炸產物均采用JWL 狀態(tài)方程描述，其形式為

或

表1 炸藥的狀態(tài)方程參數(shù)Tab.1 Coefficients of JWL state equations

1.3 反應模型

反應模型采用Lee 等提出的“點火—生長”反應模型［9］，具體形式為

式中:ρ0為炸藥的初始密度;λ 為反應度，含義為爆炸產物的質量分數(shù);I、G、x、y、z、γ 是與炸藥有關的6個可調系數(shù)，Lee 等選擇x=2/9，y =2/3，γ =4，其余3 個系數(shù)用擬合實驗數(shù)據(jù)得到。Lee－Tarver 反應率模型描述了熱點的形成和生長，(7)式中的第1 項表示熱點的形成，第2 項表述反應的增長。其中(ρ/ρ0－1)是未反應炸藥在沖擊波作用下的相對壓縮度。壓縮度越大，沖擊波所產生的熱點就越多，溫度也越高，反應也就越強。本文采用的模型系數(shù)如表2所示。

表2 “點火—生長”模型中的系數(shù)Tab.2 Coefficients for ignition and growth model

引入反應速率模型后，(2)式中源項的各分量就可寫為:

式中Q0為單位質量炸藥爆轟釋放的能量。

1.4 離散方法

本文采用格心格式的有限體積方法離散控制方程(1)式，將無粘通量與反應源項同步求解，也即由n 時刻的流動變量Qn，根據(jù)(2)式分別計算無粘通量與反應源項，然后代入(1)式，更新流動變量為Qn+1.采用Green 公式計算單元內變量梯度，獲得空間二階精度;為了抑制振蕩，采用了Venkatakrishnan 限制器;采用AUSM+－up［10］格式計算單元邊界上的通量。時間推進采用二階精度的4 步Runge－Kutta 方法。

AUSM+－up 是根據(jù)界面馬赫數(shù)對通量進行分裂的，因此用來計算不同狀態(tài)方程的流體時，只需要計算界面聲速即可，無需考慮通量的Jacobi 矩陣，很適合計算狀態(tài)方程形式復雜的問題。

2 混合物狀態(tài)方程

為了推導混合物的狀態(tài)方程，首先給出Mie－Grüneisen 狀態(tài)方程對應的聲速公式。對于滿足Mie－Grüneisen 狀態(tài)方程的流體模型，其壓力可以由密度和內能顯式地給出，狀態(tài)方程的一般形式為p=p(ρ，e).比較常見的Mie－Grüneisen 狀態(tài)方程包括理想和剛性氣體狀態(tài)方程、Van der Waals 狀態(tài)方程、JWL 狀態(tài)方程、Cochran－Chan 狀態(tài)方程和BKW、HOM 狀態(tài)方程等。

由熵的定義及熱力學第一定律可以推導得出，對于狀態(tài)方程為p =p(ρ，e)的物質，其聲速可表示為

這樣，就建立了聲速與狀態(tài)方程之間的關系，若已知物質的狀態(tài)方程，就可根據(jù)(10)式求得其聲速。

2.1 炸藥與產物的混合物狀態(tài)方程

文獻［9］指出，使用JWL 狀態(tài)方程計算爆轟問題時，在炸藥與爆轟產物間應當采用等壓假設，即同一個單元內的炸藥和爆轟產物具有相同的壓力，而溫度可能不同［11］。根據(jù)這兩個假設即可得到炸藥與產物的混合物狀態(tài)方程，過程較為簡單［11］，這里僅給出最后結果:

2.2 產物與空氣的混合物的狀態(tài)方程

爆轟產物與空氣同為氣態(tài)，在構造兩者混合模型時，借鑒了氣態(tài)爆轟計算中的模型，認為滿足局部動量平衡和熱力學平衡假設，也即在爆炸流場的任意一個流體微元內，這兩種組分的溫度、速度相同，壓力滿足分壓定理。

觀察JWL 狀態(tài)方程(5)式或(6)式還可發(fā)現(xiàn)，如果令A=B=0，則JWL 狀態(tài)方程就退化為量熱完全氣體的狀態(tài)方程，JWL 狀態(tài)方程中的ω 相當于完全氣體中的(γ－1)，γ 為比熱比。因此，可以將完全氣體狀態(tài)方程寫為JWL 狀態(tài)方程的形式，下文均采用這一方法以簡化公式的推導和書寫。

采用上面的混合模型后，爆轟產物與空氣的混合氣體的狀態(tài)方程可寫為

2.3 炸藥、產物與空氣的混合物狀態(tài)方程

當流體微元中同時含有炸藥、爆轟產物與空氣3 種組分時，炸藥為固相，而產物與空氣為氣相，它們在微元內的分布情況可如圖1所示。為了建立3組分的混合模型，將上述兩種兩組分混合模型結合起來，假設兩種氣體組分(爆轟產物與空氣間滿足動力學、熱力學平衡和分壓定理，而氣體組分與固體組分(炸藥)間具有體積可加性且壓力平衡。具體地說，有如下4 點:1)每個流體微元的體積可分為兩部分，分別為固體組分和氣體組分所占據(jù);2)氣體組分間熱力學平衡:T23=T2=T3，壓力滿足分壓定理:p23=p2+p3;3)氣體與固體組分間壓力平衡:p1=p23，溫度不等:T1≠T23;4)所有組分的速度都相等。

圖1 炸藥、產物與空氣間的流體混合模型Fig.1 Fluid－mixture model of explosives，detonation products and air

微元內混合物的分布如圖1所示。這種模型在物理上較為直觀，并且在ρ1=0 時退化為爆轟產物與空氣的混合模型，ρ3=0 時退化為炸藥與爆轟產物的混合模型。下面推導此混合模型的狀態(tài)方程:

首先考慮氣相。在微元dV23內，爆轟產物和空氣的質量分別為ρ2dV 和ρ3dV，在此微元內其密度分別為

根據(jù)(12)式可得微元dV23內壓力與內能的關系式:

然后考慮固相，在微元dV1內，炸藥的密度ρ'1=根據(jù)狀態(tài)方程(5)式，有

聯(lián)立(13)式與(14)式，根據(jù)壓力平衡條件p =p1=p23，可得到包含炸藥、爆轟產物、空氣的混合物的狀態(tài)方程:

式中:

可以證明，當ρ1=0 時，(15)式便退化為(12)式，即僅含爆轟產物與空氣的混合物狀態(tài)方程;當ρ3=0時，(15)式便退化為(11)式，即炸藥與爆轟產物的混合物的狀態(tài)方程。將(15)式分別對ρ 和e 求導，得到再根據(jù)(10)式便可得到此時的聲速。推導過程比較繁瑣，這里僅給出最后結果:

式中:

由上述推導可知，該流體混合模型可以處理包含炸藥、爆轟產物、空氣的流動問題，并且混合物的狀態(tài)方程(15)式及聲速公式(16)式均是顯函數(shù)形式，無需迭代求解，計算效率很高。該方法同時具有流體混合型方法的固有優(yōu)點，如:適合計算界面的大變形問題，無需進行任何特殊處理，給計算帶來了極大的便利，但在流體界面處也不可避免地引入了數(shù)值耗散，導致捕捉到的流體界面是一個連續(xù)過渡的區(qū)域，界面附近的若干個單元內將產生人為的流體混合。此外，本文方法使用JWL 狀態(tài)方程描述炸藥和爆轟產物，使用完全氣體模型描述空氣，如何拓展至其他形式的狀態(tài)方程仍有待研究。

3 空氣中球形裝藥的爆炸

如引言所述，模擬爆炸波在炸藥內的傳播過程時，僅需考慮炸藥與爆轟產物兩種物質，模擬爆炸波在遠場的傳播過程時，可忽略炸藥的影響，僅模擬爆轟產物與空氣，甚至可以將爆轟產物簡化為空氣處理，目前對空中爆炸問題的數(shù)值模擬也絕大多數(shù)集中于這兩種情況［3－4］，但對于爆炸近區(qū)的流場，要同時考慮炸藥、爆轟產物和空氣3 種物質，相關的數(shù)值研究很少。事實上，這類問題可供對比的實驗數(shù)據(jù)也極為有限，尤其是非球形裝藥的爆炸問題。為了校驗本文數(shù)值方法，本節(jié)計算了球形裝藥的爆炸問題，并與實驗結果進行比較。

當不考慮空氣組分時，如凝聚態(tài)炸藥中的爆轟波傳播問題，(15)式自動退化為(11)式，此時計算程序的精度和有效性已得到了驗證。但對于同時包含炸藥、爆轟產物、空氣3 種組分的問題，可供驗證的實驗數(shù)據(jù)尤其是非球形爆炸的近區(qū)實驗數(shù)據(jù)幾乎沒有，為了校驗本文方法，模擬了空氣中球形裝藥的爆炸問題，重點考察爆炸近區(qū)的計算結果與實驗值的比較。

采用一維球對稱程序模擬了空氣中球形裝藥的爆炸問題。裝藥為TNT，裝藥半徑為0.1 m，裝藥量約為6.8 kg，在球心處點火。計算采用的反應模型為“點火—生長”模型，狀態(tài)方程及反應模型中的系數(shù)來自文獻［12］.如圖2所示給出了不同時刻的壓力波形，可看出，15 μs 時，爆轟波仍在炸藥中傳播，未到達炸藥與空氣界面，爆轟波的峰值壓力為46 GPa，爆速為6 810 m/s，略低于理論爆速6 849 m/s.C－J 點壓力為18.86 GPa，也略低于理論值19.1 GPa.24 μs 時，沖擊波已傳至空氣中，由于波后無反應區(qū)提供能量，沖擊波波后壓力迅速下降，遠低于爆轟波波后壓力。從圖上可以看出，隨著傳播距離的增大，沖擊波的波后壓力逐漸下降，波速也在變慢。

圖2 不同時刻的爆轟波壓力波形Fig.2 Distributions of pressure at different times

如圖3所示給出了42 μs 時沖擊波附近的壓力和密度曲線，在此圖上可更清楚地看到爆轟波到達空氣界面后的波系發(fā)展情況。爆轟波到達產物與空氣的界面后，透射波為沖擊波，在空氣中傳播;反射波為反射膨脹波，向爆心傳播。爆轟波后的強稀疏波到達產物與空氣的界面后，透射波為透射稀疏波，反射波為反射激波。透射稀疏波不停地追上空氣中沖擊波，使得波后峰值壓力不斷被削弱，沖擊波波速下降。同樣的，向爆心運動的反射激波不停追上反射膨脹波，在激波后產生反向運動的稀疏波。由此可見，計算得到的各種波系及流體界面結構清晰，反映了計算方法具有較高的精度和分辨率。

計算過程中沿徑向每隔0.05 m 設置一個壓力監(jiān)測點，如圖4所示給出了標尺距離為1.0、1.5、2.0 的3 個測點的超壓時程曲線。標尺距離的定義為Rs=r/W1/3，r 為測點到爆心的距離(m)，W 為裝藥量(kg)。

圖3 42 μs 時的壓力和密度曲線Fig.3 Distribution curves of pressure and density at 42 microseconds

關于空中爆炸的超壓峰值隨距離的變化，有相當多的實驗數(shù)據(jù)與經驗公式可供比較，但適用于爆炸近區(qū)的數(shù)據(jù)較少。本文主要采用了文獻［13－16］中的一些數(shù)據(jù)，這些數(shù)據(jù)在離爆心較遠的距離上能給出彼此接近的結果，但在裝藥近區(qū)有所差異。其中Baker 的數(shù)據(jù)是對許多理論和實驗結果的匯編，具有很寬的適用范圍;Henrych［14］和Brode［16］的經驗公式在爆炸近區(qū)會給出偏大的結果;Henrych的經驗公式在直至裝藥表面的范圍內都是足夠準確的。如圖5所示給出了計算得到的超壓峰值與實驗結果的比較，可看出，計算結果與各實驗值均符合很好，即使標尺距離小于時，計算結果仍然與Baker、Henrych、Kinney 等［15］的結果高度吻合，曲線中標尺距離最小的點距裝藥表面僅5 cm.

圖4 標尺距離1.0、1.5、2.0 的3 個監(jiān)測點的超壓時間曲線Fig.4 Over－pressure waveforms at scaled distances of 1.0，1.5，2.0

圖5 超壓峰值與實驗結果、經驗公式的比較Fig.5 Comparison of peak over－pressure between computation and experiment and formula

取爆轟波到達裝藥表面的時刻為零時刻，如圖6所示給出了沖擊波到達的時間隨標尺距離的變化情況，與Baker 的結果符合很好。如圖7所示給出了超壓沖量隨標尺距離的變化情況，作為對比的分別是Baker 的匯編數(shù)據(jù)，Sadofskyi 的理論公式以及Henrych 的實驗結果。在這里，不同文獻的結果出現(xiàn)了數(shù)倍的差異，而本文計算結果落在了文獻結果之間，與Baker 的結果更為接近，且變化趨勢與Baker 的結果相同。

總的來說，計算結果與各種實驗結果符合較好，反映了本文計算方法具有相當?shù)臏蚀_性。

圖6 沖擊波到達時間的比較Fig.6 Comparison of arriving time

圖7 超壓沖量的比較Fig.7 Comparison of impulse

4 結論

提出了一種用于模擬爆炸近區(qū)中炸藥、爆轟產物、空氣等相互作用流場的數(shù)值方法。算法為固定網格下的歐拉型方法，引入了一種炸藥、爆轟產物、空氣3 種組分的混合模型，該模型中炸藥及爆轟產物采用JWL 狀態(tài)方程，空氣采用理想氣體模型，在固相(炸藥)與氣相(爆轟產物、空氣)間采用等壓假設，且體積可加;氣相間滿足等溫假設及分壓定理。該模型無需迭代求解，計算效率較高。計算了空氣中球形TNT 裝藥的爆轟問題，得到了清晰的波系結構及流體界面，界面附近無非物理振蕩，計算結果與現(xiàn)有實驗值符合良好，驗證了本文方法的有效性和計算精度，為下一步的應用打下了基礎。

References)

［1］ Zhou X Q，Hao H.Prediction of airblast loads on structures behind a protective barrier［J］.International Journal of Impact Engineering，2008，35(5):363－375.

［2］ Wu C，Hao H.Numerical simulation of structural response and damage to simultaneous ground shock and airblast loads［J］.International Journal of Impact Engineering，2007，34(3):556－572.

［3］ 寧建國，王仲琦，趙衡陽，等.爆炸沖擊波繞流的數(shù)值模擬研究［J］.北京理工大學學報，1999，19(5):543－547.NING Jian－guo，WANG Zhong－qi，ZHAO Heng－yang，et al.Study on the flow of the explosive shock wave around the wall numerical simulation［J］.Transactions of Beijing Institute of Technology，1999，19(5):543－547.(in Chinese)

［4］ 劉君，劉瑞朝，賈忠湖，等.爆炸波與物體干擾流場的數(shù)值模擬［J］.空氣動力學學報，2000，18(1):55－60.LIU Jun，LIU Rui－chao，JIA Zhong－h(huán)u，et al.Numerical simulation of blast interaction with bulidings［J］.Acta Aerodynamica Sinica，2000，18(1):55－60.(in Chinese)

［5］ Liu T G，Khoo B C，Yeo K S.Ghost fluid method for strong shock impacting on material interface［J］.Journal of Computational Physics，2003，190(2):651－681.

［6］ Saurel R，Abgrall R.A multiphase godunov method for compressible multifluid and multiphase flows［J］.Journal of Computational Physics，1999，150(2):425－467.

［7］ Abgrall R，Saurel R.Discrete equations for physical and numerical compressible multiphase mixtures［J］.Journal of Computational Physics，2003，186(2):361－396.

［8］ Chinnayya A，Daniel E，Saurel R.Modelling detonation waves in heterogeneous energetic materials［J］.Journal of Computational Physics，2004，196(2):490－538.

［9］ Lee E L，Tarver C M.Phenomenological model of shock initiation in heterogeneous explosives［J］.Physics of Fluids，1980，23(12):2362－2372.

［10］ Liou M－S.A sequel to AUSM，part II:AUSM+－up for all speeds［J］.Journal of Computational Physics，2006，214(1):137－170.

［11］ Clutter J K，Belk D.Simulation of detonation wave interaction using an ignition and growth model［J］.Shock Waves，2002，12(3):251－263.

［12］ Kury J W，Breithaupt R D，Tarver C M.Detonation waves in trinitrotoluene［J］.Shock Waves，1999，9(9):227－237.

［13］ Baker W E.Explosions in air:part 1［M］.Austin:University of Texas Press，1974.

［14］ Henrych J.The dynamics of explosion and its use［M］.Amsterdam:Elsevier Scientific Publishing Company，1979.

［15］ Kinney G F，Graham K J.Explosive shocks in air［M］.New York:Spring－Verlag，1985.

［16］ Brode H L.Numerical solution of spherical blast waves［J］.Journal of Applied Physics，1955，26(6):766－775.