空氣自由場(chǎng)中強(qiáng)爆炸沖擊波傳播二維數(shù)值模擬

姚成寶，李 若，田 宙，郭永輝

(1.北京大學(xué)數(shù)學(xué)科學(xué)學(xué)院，北京 100871； 2.西北核技術(shù)研究所，陜西 西安 710024； 3.北京大學(xué)應(yīng)用物理與技術(shù)研究中心，北京 100871)

?

空氣自由場(chǎng)中強(qiáng)爆炸沖擊波傳播二維數(shù)值模擬

姚成寶1,2，李 若1,3，田 宙1，郭永輝1

(1.北京大學(xué)數(shù)學(xué)科學(xué)學(xué)院，北京 100871； 2.西北核技術(shù)研究所，陜西 西安 710024； 3.北京大學(xué)應(yīng)用物理與技術(shù)研究中心，北京 100871)

編寫了適用于模擬具有高密度比、高壓力比的強(qiáng)激波問題的二維柱對(duì)稱多介質(zhì)流體計(jì)算程序。利用有限體積方法求解流體的Euler方程組，采用level set方法捕捉爆炸產(chǎn)物與空氣的運(yùn)動(dòng)界面，并通過求解物質(zhì)界面兩側(cè)Riemann問題的精確解來計(jì)算爆炸產(chǎn)物與空氣之間的數(shù)值通量。研制了三角形網(wǎng)格自適應(yīng)技術(shù)來實(shí)現(xiàn)網(wǎng)格的自動(dòng)加密和粗化，在保證捕捉激波峰值的前提下有效地提高了計(jì)算效率。利用計(jì)算程序?qū)? kt TNT當(dāng)量的空氣自由場(chǎng)強(qiáng)爆炸問題進(jìn)行數(shù)值模擬，計(jì)算得到的峰值超壓、沖擊波到達(dá)時(shí)間等物理參數(shù)與點(diǎn)爆炸理論結(jié)果基本一致。

爆炸力學(xué)；空中強(qiáng)爆炸；level set；Riemann問題；網(wǎng)格自適應(yīng)；自由場(chǎng)；沖擊波

空中爆炸[1]是指炸藥、炮彈等在距地面或水面一定高度處的爆炸。隨著計(jì)算機(jī)技術(shù)和計(jì)算理論的快速發(fā)展，數(shù)值模擬在該領(lǐng)域內(nèi)的研究起到越來越重要的作用[2-3]。和常規(guī)炸藥相比，強(qiáng)爆炸產(chǎn)生的初始?jí)毫Ω撸S傳播距離衰減更慢，殺傷破壞距離更遠(yuǎn)，數(shù)值模擬的難度更大。

針對(duì)空中強(qiáng)爆炸問題的強(qiáng)對(duì)流性，本文中選擇基于歐拉坐標(biāo)系下的多介質(zhì)流體計(jì)算模型模擬沖擊波的傳播過程。文獻(xiàn)[4]中認(rèn)為爆炸產(chǎn)物與空氣的界面最初是分開的，只是隨后由于爆炸產(chǎn)物的脈動(dòng)作用，特別是由于界面周圍產(chǎn)生渦流作用，使得界面逐漸模糊，最后與空氣混合在一起。文獻(xiàn)[4]的實(shí)驗(yàn)結(jié)果表明，對(duì)爆炸破壞作用具有實(shí)際意義的只是第一次的膨脹和壓縮過程。基于上述觀點(diǎn)，本文中采用level set方法捕捉爆炸產(chǎn)物與空氣的界面，并認(rèn)為界面兩側(cè)的物質(zhì)不發(fā)生相互融合。大量的數(shù)值結(jié)果表明[5-7]，在物質(zhì)界面兩側(cè)，由于爆炸產(chǎn)物與空氣具有不同的狀態(tài)方程，且密度、壓力相差巨大，需采取合適的方法來處理界面上的接觸間斷，否則會(huì)產(chǎn)生強(qiáng)烈的非物理振蕩。本文中通過求解Riemann精確解來獲得界面上接觸間斷兩側(cè)的物質(zhì)狀態(tài)，并以此求解界面上的數(shù)值通量，以避免非物理振蕩。

在空中強(qiáng)爆炸沖擊波傳播數(shù)值模擬過程中，為準(zhǔn)確捕捉?jīng)_擊波壓力峰值，本文中采用網(wǎng)格自適應(yīng)技術(shù)來實(shí)現(xiàn)網(wǎng)格的自動(dòng)加密和稀疏；利用計(jì)算程序?qū)? kt TNT當(dāng)量空中強(qiáng)爆炸產(chǎn)生的沖擊波傳播過程進(jìn)行數(shù)值模擬，將計(jì)算得到的沖擊波參數(shù)與點(diǎn)爆炸理論結(jié)果[8]進(jìn)行對(duì)比。

1 控制方程和狀態(tài)方程[9]

由于空中強(qiáng)爆炸問題本身具有的對(duì)稱性，本文采用二維柱對(duì)稱計(jì)算模型來進(jìn)行計(jì)算。控制方程為：

(1)

式中：ρ表示密度，u和v表示方向r和z方向的速度，E表示總能量密度，p表示壓力。

爆炸產(chǎn)物和空氣均采用γ律狀態(tài)方程描述：

(2)

式中：γ為比熱比；對(duì)于強(qiáng)爆炸產(chǎn)物，γ=1.2；對(duì)于空氣，γ=1.4。e為比內(nèi)能，且滿足：

(3)

2 計(jì)算方法

2.1 level set方法[10-12]

level set方法的主要思想是引入level set距離函數(shù)φ(x,t)，把隨時(shí)間運(yùn)動(dòng)的物質(zhì)界面Γ(t)定義為距離函數(shù)的零等值面。在每個(gè)時(shí)刻，只要求出距離函數(shù)φ(x,t)的值，就可以知道物質(zhì)界面的位置。level set函數(shù)的演化方程為：

(4)

式中：φt、φx和φy分別為φ對(duì)t、x和y的偏導(dǎo)數(shù)。在計(jì)算過程中，為保持φ(x,t)始終是x點(diǎn)到界面Γ(t)的符號(hào)距離，需要對(duì)φ(x,t)進(jìn)行重新初始化：

(5)

式中：φ0為重新初始化前的φ(x,t)的值，且

(6)

式中：ε為一個(gè)任意小的正數(shù)。

2.2 Riemann問題求解[13-14]

圖1 Riemann問題解的類型Fig.1 Solution type of the Riemann problem

律狀態(tài)方程的Riemann問題為求解式(7)所示的多介質(zhì)Euler方程組的初值問題：

(7)

式中：下角標(biāo)1和2分別表示接觸間斷左側(cè)和右側(cè)。

這種初始時(shí)刻的間斷是不穩(wěn)定的，在t>0時(shí)會(huì)立即分解為圖1所示的波系結(jié)構(gòu)。其中,w表示線性退化的接觸間斷波；w1和w2分別表示接觸間斷兩側(cè)的非線性波(根據(jù)初值條件的不同，可能是稀疏波或者激波)；ρI1和ρI2分別為接觸間斷上w左側(cè)和右側(cè)的密度。在接觸間斷兩側(cè)，壓力和速度保持不變，而密度發(fā)生間斷。記W=(ρ,u,p)T，則壓力p滿足：

(8)

其中,

(9)

(10)

(11)

利用牛頓迭代法對(duì)式(8)進(jìn)行迭代求解，可得到界面上的最終壓力p。界面上的速度u的表達(dá)式為

(12)

3 數(shù)值離散方法

3.1 守恒方程組的離散求解

采用有限體積方法來離散求解守恒方程組，在每個(gè)網(wǎng)格單元τ上從tn到tn+1時(shí)刻對(duì)式(1)進(jìn)行積分，可得：

(13)

(14)

式中：Ut、Fx、Gy分別為U、F、G對(duì)t、x、y的偏導(dǎo)數(shù)，Ω為網(wǎng)格單元的體積。利用散度定理，可得：

(15)

式中：?τ為τ的表面，S為表面積。對(duì)于包含2種流體的混合單元，除了需要計(jì)算單元邊界上的數(shù)值通量外，還需要在物質(zhì)界面上利用Riemann問題的解來計(jì)算2種流體在界面上的數(shù)值通量。

3.2 level set方程求解

本文中利用特征線方法來數(shù)值求解level set函數(shù)的演化方程，將式(4)轉(zhuǎn)換成:

(16)

根據(jù)特征線方法的思想，只需知道tn+1時(shí)刻網(wǎng)格頂點(diǎn)Xi上的流體質(zhì)點(diǎn)在tn時(shí)刻的位置 ，就可以得到tn+1時(shí)刻的level set函數(shù)在Xi處的值φn+1(Xi)。具體求解過為：

(17)

(18)

式中：Vn[X(tn)]表示tn時(shí)刻X(tn)處流體質(zhì)點(diǎn)的速度。

圖2 網(wǎng)格自適應(yīng)示意圖Fig.2 Mesh adaption

利用文獻(xiàn)[15]中的顯式正系數(shù)格式來進(jìn)行求解level set函數(shù)的重新初始化方程(5)，保持φ(x,t)始終是x點(diǎn)到界面Γ(t)的符號(hào)距離。

3.3 網(wǎng)格自適應(yīng)技術(shù)

為了有效地捕捉激波峰值，并盡可能減少計(jì)算量，本文中采用網(wǎng)格自適應(yīng)技術(shù)。根據(jù)相鄰網(wǎng)格的壓力梯度來構(gòu)造自適應(yīng)指示因子，在壓力變化劇烈的地方加密網(wǎng)格，在變化平緩的地方對(duì)網(wǎng)格進(jìn)行粗化。典型的網(wǎng)格自適應(yīng)如圖2所示。

由圖2可知，在沖擊波波陣面附近，為捕捉激波峰值，網(wǎng)格被劇烈的加密；而在沖擊波波陣面后以及沖擊波尚未達(dá)到的地方，網(wǎng)格非常稀疏。

4 算 例

對(duì)1 kt TNT當(dāng)量空中強(qiáng)爆炸產(chǎn)生的沖擊波傳播過程進(jìn)行數(shù)值模擬，并將計(jì)算得到的沖擊波參數(shù)與點(diǎn)爆炸理論結(jié)果[8]進(jìn)行對(duì)比。

(1) 沖擊波壓力等值線圖。

計(jì)算得到不同時(shí)刻的沖擊波壓力等值線如圖3所示。

圖3 不同時(shí)刻的沖擊波壓力等值線圖Fig.3 Pressure contours at different time

(2) 沖擊波參數(shù)。

將計(jì)算得到的不同爆心距(r)的沖擊波超壓峰值(pop)、正壓沖量(I)、沖擊波到時(shí)(ts)和正壓作用時(shí)間(Δt)等沖擊波參數(shù)與點(diǎn)爆炸理論結(jié)果進(jìn)行比對(duì)，如圖4～7所示。

圖4 峰值超壓Fig.4 Peak overpressure

圖5 正壓沖量Fig.5 Overpressure impulse

圖6 沖擊波到時(shí)Fig.6 Arrival time

圖7 正壓作用時(shí)間Fig.7 Positive time duration

從圖4～7可知，程序計(jì)算得到的峰值超壓、沖擊波到時(shí)與點(diǎn)爆炸理論結(jié)果基本一致，正壓作用時(shí)間略有差別。由于點(diǎn)爆炸理論中沒有給出正壓沖量的結(jié)果，這里只給出程序的計(jì)算結(jié)果。

(3) 沖擊波超壓時(shí)間歷程曲線。

計(jì)算得到不同爆心距的沖擊波超壓(po)歷程曲線，如圖8所示。從圖8可知，對(duì)離爆心一定距離的觀察點(diǎn)，當(dāng)沖擊波到達(dá)時(shí)，空氣壓力迅速上升至峰值；沖擊波過后，空氣壓力逐漸降低；當(dāng)爆炸產(chǎn)物過度膨脹形成的負(fù)壓稀疏波到達(dá)時(shí)，空氣呈稀疏狀態(tài)，空氣壓力逐漸下降至負(fù)壓。隨后由于爆炸產(chǎn)物的脈動(dòng)作用，該處再次形成二次激波。如此反復(fù)，直至最后恢復(fù)至初始?jí)毫Α膱D8中給出的距爆心為200 m和300 m處的沖擊波超壓歷程來看，計(jì)算得到的沖擊波超壓曲線能夠正確地描述空中強(qiáng)爆炸時(shí)某固定位置處沖擊波的正壓、負(fù)壓以及二次激波，表明程序的計(jì)算合理，結(jié)果可信。

5 結(jié) 論

針對(duì)空中強(qiáng)爆炸問題的特點(diǎn)，研究了適合計(jì)算強(qiáng)對(duì)流問題的數(shù)值計(jì)算方法，在此基礎(chǔ)上編寫了相應(yīng)的二維多介質(zhì)流體計(jì)算程序。其中，采用有限體積方法求解流體控制方程，利用level set方法捕捉流體的運(yùn)動(dòng)界面，并通過求解物質(zhì)界面兩側(cè)的Riemann問題精確解來計(jì)算物質(zhì)界面兩側(cè)的數(shù)值通量。計(jì)算過程中采用了網(wǎng)格自適應(yīng)技術(shù)，在捕捉激波峰值的前提下，有效地減小了計(jì)算量，提高了計(jì)算效率。

計(jì)算了1 kt TNT當(dāng)量的空中強(qiáng)爆炸問題，模擬了沖擊波在空氣自由場(chǎng)中的傳播過程，給出了從爆炸源區(qū)附近到1 km距離處各爆心距的超壓峰值、正壓沖量、沖擊波到時(shí)和正壓作用時(shí)間等沖擊波參數(shù)，計(jì)算結(jié)果與點(diǎn)爆炸理論結(jié)果基本一致，說明本文程序采用的計(jì)算方法能夠應(yīng)用于空中爆炸等強(qiáng)對(duì)流問題的數(shù)值模擬。

[1] 奧爾連科 Л П.爆炸物理學(xué)[M].孫承緯,譯.北京:科學(xué)出版社,2013,456-463.

[2] 張洪武,何揚(yáng),張昌權(quán).空中爆炸沖擊波地面載荷的數(shù)值模擬[J].爆炸與沖擊,1992,12(2):188-198. Zhang Hong-wu, He Yang, Zhang Chang-quan. Numerical simulation on ground surface loading of shock wave from air explosion[J]. Explosion and Shock Waves, 1992,12(2):188-198.

[3] 岳鵬濤,徐勝利,彭金華.FAE爆炸波對(duì)地面目標(biāo)作用的三維數(shù)值模擬[J].爆炸與沖擊,2000,20(2):195-201. Yue Peng-tao, Xu Sheng-li, Peng Jin-hua. 3D numerical simulations on the interaction between FAE blast waves and ground targets[J]. Explosion and Shock Waves, 2000,20(2):195-201.

[4] 北京工業(yè)學(xué)院.爆炸及其作用[M].北京:國防工業(yè)出版社,1979.

[5] Liu T G, Khoo B C, Yeo K S. Ghost fluid method for strong shock impactingon material interface[J]. Journal of Computational Physics, 2003,190(2):651-681.

[6] Fedkiw R P, Aslam T, Merriman B,et al. A non-oscillatory eulerian approach to interfaces in multimaterial flows: The ghost fluid method[J]. Journal of Computational Physics, 1999,152(2):457-492.

[7] Colella P, Glaz H M. Efficient solution algorithms for the Riemann problem for real gases[J]. Journal of Computational Physics,1985,59(2):264-289.

[8] 喬登江.強(qiáng)爆炸物理概論(上冊(cè))[M].北京:國防工業(yè)出版社,2003:51-55.

[9] 楊秀敏.爆炸沖擊現(xiàn)象數(shù)值模擬[M].合肥:中國科學(xué)技術(shù)大學(xué)出版社,2010:122-125.

[10] Osher S, Sethian J A. Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-jaccobi formulations[J]. Journal of Computational Physics, 1988,79(1):12-49.

[11] Sussman M, Semerka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow[J]. Journal of Computational Physics, 1994,114(4):146-159.

[12] 張軍,趙寧,任登鳳,等.level set方法在自適應(yīng)Catersian網(wǎng)格上的應(yīng)用[J].爆炸與沖擊,2008,28(5):156-161. Zhang Jun, Zhao Ning, Ren Deng-feng, et al. Application of the level set method on adaptive cartesian grids[J].Explosion and Shock Waves, 2008,28(5):156-161.

[13] Eleuteuo F T. Riemann solvers and numerical methods for fluid dynamics[M]. First Edition. Springer, 2009:16-123

[14] 劉儒勛,舒其望.計(jì)算流體的若干新方法[M].北京:科學(xué)出版社,2003:121-130

[15] Di Y N, Li R, Tang T, et al. level set calculations for incompressible two-phase flows on a dynamically adaptive grid[J]. Journal of Scientific Computing, 2007,31(1/2):75-98.

(責(zé)任編輯 王小飛)

Two dimensional simulation for shock wave produced by strong explosion in free air

Yao Cheng-bao1, 2, Li Ruo2, Tian Zhou1, Guo Yong-hui1

(1.SchoolofMathematicalSciences,PekingUniversity,Beijing100871,China; 2.NorthwestInstituteofNuclearTechnology,Xi’an710024,Shaanxi,China; 3.CenterofAppliedPhysicsandTechnology,PekingUniversity,Beijing100871,China)

Aimed to simulate the propagation of blast wave with high density ratio and high pressure ratio produced by strong explosion in the air, a two dimensional numerical program is written in which the problem is treated as a two-medium compressible flow with sharp material interface in Eulerian grids. In this method, the finite volume method is used to solve the Euler equations, level set method is used to capture the moving interface, and the numerical flux across the interface is calculated by exactly solving the Riemann problem. Mesh adaption technique in triangle meshes is adopted to refine or coarsen the meshes which can both capture the peak overpressure and improve the computational efficiency. One kiloton nuclear charge of strong explosion in free air is simulated. The shock wave parameters, including peak overpressure, shock arrival time and so on, are in consistence with the point explosion theory, which show the accuracy and efficiency of the numerical methods.

mechanics of explosion; strong explosion in air; level set; Riemann problem; mesh adaption; free air; shockwave

10.11883/1001-1455(2015)04-0585-06

2013-12-04；

2014-07-09

姚成寶(1984- )，男，博士研究生，助理研究員，yaocheng@pku.edu.cn。

O381 國標(biāo)學(xué)科代碼： 13035

A