吳宗鐸，嚴 謹，蔣 頡，董美余，黃 技

(1. 廣東海洋大學 海洋工程學院，廣東 湛江 524088；2. 中國船舶研究設計中心，上海 201108；3. 上海振華重工集團股份有限公司，上海 200125)

基于Mie-Grüneisen狀態(tài)方程的水下爆炸數值模擬

吳宗鐸1，嚴 謹1，蔣 頡2，董美余3，黃 技1

(1. 廣東海洋大學 海洋工程學院，廣東 湛江 524088；2. 中國船舶研究設計中心，上海 201108；3. 上海振華重工集團股份有限公司，上海 200125)

由于Mie-Grüneisen狀態(tài)方程形式比較復雜，給界面的處理帶來很大難度。本文通過對Euler方程做分離變量和引入質量分數，完成了Mie-Grüneisen狀態(tài)方程下的多介質可壓縮流動的數值模擬，使計算過程得到簡化，并通過算例驗證了該方法的可靠性。在多項式形式的Mie-Grüneisen狀態(tài)方程下，利用該方法模擬了球形炸藥在水中爆炸后，氣體和水相互作用的近場情況。在計算模型中引入的氣體質量分數，很好反映了流場中不同區(qū)域內氣體、水及水氣并存3種不同的狀態(tài)。

水下爆炸；數值計算；Mie-Grüneisen狀態(tài)方程；多介質

0 引 言

由于水下爆炸在軍事和國防領域有著極大的重要性，因而一直受到大家的關注。多數時候，炸藥都以球形裝藥的形式出現。對于球形炸藥來說，得到其沖擊波壓力并不容易。早期對沖擊波壓力的研究以實驗為主[1 – 2]。其后，雖然出現了許多數值模擬，但水下爆炸的實驗仍未中斷[3 – 4]。但考慮保密性的要求，很多結果未能發(fā)表。

然而，由于實驗條件成本大，實施難，水下爆炸已經越來越多的被數值模擬所代替。但球形炸藥的三維特性使得它在普通的數值模擬中很難得到廣泛應用。其中很大一個原因就是普通的數值模擬中設置三維計算網格將消耗更多的計算機資源。因此，關于水下爆炸的研究工作多基于商業(yè)軟件，如Autodyn[5 – 6]，LS-DYNA[7]，Abaqus[8]，MSC-Dytran[9]等。而目前的水下爆炸數值模擬仍然多停留在二維平面上，如張阿曼的無網格SPH方法[10]、Daramizadeh[11]的五方程方法。除此以外，還有一個原因也對球形炸藥的數值模擬起著限制作用，那就是涉及到復雜狀態(tài)方程下的Riemann問題求解。目前對于水下爆炸問題，采用的狀態(tài)方程多為形式復雜Mie-Grüneisen狀態(tài)方程[5, 7, 10]，給沖擊波的模擬造成很大的難度。

本文將求解坐標系由傳統的二維網格體系轉到球坐標體系下。由于藥球爆炸后沖擊波傳播的球面對稱形式，沖擊波的參數僅和傳播距離有關[12]。這樣，可將沖擊波物理參數表示成球坐標系半徑的函數，無需建立三維的直角坐標網格來求解，為計算帶來很大的方便。對于Riemann問題的處理，本文在Mie-Grüneisen狀態(tài)方程下對能量守恒Euler方程做變量分離，通過近似處理把參考壓力、參考單位內能等復雜的物理量分離出來，使Euler方程組簡單易于計算。利用這種算法，對水下爆炸的近場情況進行數值計算。在激波的傳播過程中，整個流場被視為一個理想氣體和水的混合體，當中的氣體和水用質量分數加以區(qū)分。由于藥球爆炸后，水和水蒸氣共存的狀態(tài)大量存在于流場周圍，該方法能很好地反映這一狀態(tài)，同時也能較好地計算出水下爆炸周圍流場的情況。

1 控制方程

水下爆炸問題中，無論是炸藥還是介質，其物理參數密度ρ、壓力p和單位體積內能ρe都可寫成Mie-Grüneisen狀態(tài)方程的形式[12]，其表達式如下：

其中，Г為Mie-Grüneisen系數，可以表示成

因此，Г可以寫成密度ρ的函數。pref和eref為參考點的壓力和單位體積內能，與密度ρ有關。假設流場無粘，那么根據質量、動量和能量守恒關系，可以得出歐拉方程：

式中：ρ，p，u，E分別為密度，壓力，速度和單位體積能量，并且對于單位體積能量E，可以寫成內能和動能的總和，于是有

在炸藥和水的接觸間斷面上，大多數物理參數為間斷狀態(tài)但流體的壓力和速度仍然保持連續(xù)，間斷的物理參數導數可以視為無窮大，因此u和p的導數相對其他導數來說可以忽略不計，那么式（3）中質量守恒和能量守恒可以近似表示成：

將方程（1）代入式（3）中的能量守恒方程，有

由于Г，pref和eref都可以寫成只和密度ρ有關的數，再忽略p的變化，因此，可以對式（5）作變量分離：

這樣Mie-Gruneisen狀態(tài)方程下的m種介質的控制方程可以表示成包含ρ，p，u，E和Г，pref，eref以及質量分數αi的如下的歐拉關系式。方程（3）中的變量U可以寫成：

通量F則表示成：

2 計算方法

本文采用MUSCL-TVD格式的有限體積法進行計算。假設問題為一維Euler方程，其有限體積形式有：

3 計算實例

半徑為R0裸藥球在無限水域中爆炸，假設流體無粘且藥球質量分布均勻，t=0時刻，炸藥爆炸并形成高壓氣體，氣體密度ρ0等于炸藥球的密度，即1 630 kg/m3，爆壓pH=2.1×1010Pa，氣團半徑等于藥球半徑[1]。假設TNT炸藥爆炸后形成的氣體為理想氣體，狀態(tài)方程如下：

式中：μ=ρ/ρ0–1，其系數如表1所示。

這樣，轉化成如（1）所示的Mie-Grüneisen的形式以后，相關參數為：

表 1 炸藥和水的狀態(tài)方程參數Tab. 1 The EOS parameters for explosive and water

圖1為壓力峰值與文獻[2]的對比結果。從圖中可以看出，2條壓力峰值曲線與大致保持一致，且在距爆點約12R0處吻合得比較好。當沖擊波過了12R0以后，沖擊波的壓力衰減速度降低，曲線將有所緩和，與本文的計算存在少量差距。而在距爆點很近的地方（6R0以內），誤差會有所增加。由于這部分區(qū)域的爆炸沖擊波的作用很強，壓力變化非常劇烈，有些實驗結果甚至未能在6R0以內給出合理的經驗公式[1]。

圖2為距爆心不同距離處的壓力峰值曲線，橫坐標是測點與爆心的距離R和初始半徑R0的比值。本文計算結果與文獻[2]的經驗公式相比，近場壓力值總體有些偏大。但基本反映了壓力在沖擊波前后的變化情況。在沖擊波傳播到某點時，該點壓力瞬間上升到峰值；隨后沖擊波繼續(xù)前進，此時由于沖擊波作用減弱該點壓力則迅速下降。隨后，向內收縮的稀疏波聚集到爆心處又重新擴散開來，造成壓力的脈動變化，形成第2次峰值。

圖3為不同時刻壓力的分布圖，從圖中可以看出，隨著時間的推移，壓力曲線逐漸變得平穩(wěn)，并在最大峰值過后一定距離出現了壓力脈動，形成二次脈動壓力。圖4給出了不同位置的二次壓力峰值，可以看出其數值比較小。值得注意的是，由于本文的計算是基于能量守恒的假設，未考慮能量的損耗，因此本文的計算結果比實際情況要大，由脈動壓力形成的第2次的壓力峰值也比較明顯。實際情況下，第2次壓力峰值很難在曲線上清楚顯現出來[1]。

圖5理想氣體的質量分數α和質量ρα在不同時刻的變化情況，從圖6中可以看出，沖擊波傳播開來以后，中央的氣團逐漸由單一的氣體狀態(tài)變?yōu)闅怏w和液態(tài)水并存的一種狀態(tài)。這情況與實際的氣水模型存在一定差別，實際的氣水模型中并未有爆轟氣體溶解在水中，因此本計算模型仍然存在著一些數值上的不準確，但并不影響沖擊波壓力的計算。

4 結 語

利用本文的分離變量方法能很好地解決Mie-Grüneisen狀態(tài)方程下的多介質可壓縮流動的激波問題，并且能有效地簡化計算過程。該方法在處理水下爆炸近場問題的時候能取得很好的效果。按本文的計算方法處理得到的計算結果和文獻[1]對比，取得了不錯的效果，并對計算第2次的脈動壓力有著很好的幫助。

炸藥在水中爆炸以后，水的狀態(tài)包括液態(tài)水、水蒸氣和液氣混合態(tài)。再加上爆炸產生的氣體、液態(tài)水、氣體和氣水混合體共存于流場中。利用質量分數可以區(qū)分開流場中的各種成分。炸藥爆炸后，早期形成高密度氣團，隨后氣團擴散很快，在氣團擴散過程中，氣體的質量分數和密度都急劇下降，最后形成一個中心大部分區(qū)域稀薄，兩邊小部分區(qū)域稠密的狀態(tài)。

由于本文僅考慮了動量和能量守恒的理想情況，而實際情況比較復雜，能量損耗不可避免，因此，本文所計算的壓力值比實際數值大，沖擊波的傳播速度也比實際要快。

[ 1 ]COLE R. H. Underwater explosions[M]. New York: Dover Publications, INC. 1965.

[ 2 ]ZAMYSHLYAYEV B V, YAKOVLEV Y S. Dynamic loads in underwater explosion, AD-757183 [R]. Naval Intelligence Support Center, 1978.

[ 3 ]楊振, 沈曉樂. 爆破戰(zhàn)斗部水中兵器爆炸威力評定方法研究[J]. 爆破, 2015, 32(2): 51–53. YANG Zhen, SHEN Xiao-le. Research on evaluation method of underwater blast brisance of weapon's explosive[J]. Blasting, 2015, 32(2): 51–53.

[ 4 ]魯忠寶, 胡宏偉, 劉銳, 等. 典型裝藥水下爆炸的殉爆規(guī)律研究[J]. 魚雷技術, 2014, 22(3): 230–235. LU Zhong-bao, HU Hong-wei, LIU Riu, et al. Research on law of sympathetic detonation of typical charge[J]. Torpedo Technology, 2014, 22(3): 230–235.

[ 5 ]張社榮, 李宏璧, 王高輝, 等. 水下爆炸沖擊波數值模擬的網格尺寸確定方法[J]. 振動與沖擊, 2015, 34(8): 93–100. ZHANG She-rong, LI Hong-bi, WANG Gao-hui, et al. A method to determine mesh size in numerical simulation of shock wave of underwater explosion[J]. Journal of Vibration and Shock, 2015, 34(8): 93–100.

[ 6 ]姜濤, 王桂芹, 詹發(fā)民, 等. 基于AUTODYN的潛艇典型艙段水中爆炸沖擊損傷研究[J]. 爆破器材, 2015, 44(6): 61–64. JIANG Tao, WANG Gui-qin, ZHAN Fa-min, et al. Impact damage analysis of typical submarine compartment subjected to underwater blasting based on AUTODYN[J]. Explosive Materials, 2015, 44(6): 61–64.

[ 7 ]梅群, 侯中華, 朱俊鋒, 等. 水下爆炸沖擊波壓力時程的數值模擬[J]. 河南科技大學學報(自然科學版), 2010, 31(4): 57–59 MEI Qun, HOU Zhong-Hua, ZHU Jun-Feng, et al. Numerical simulation of underwater explosion shock wave[J]. Journal of Henan University of Science and Technology: Natural Science, 2010, 31(4): 57–59.

[ 8 ]賈則, 權琳, 張姝紅, 等. 艦艇結構在水中兵器靜態(tài)爆炸作用下的沖擊響應[J]. 艦船科學技術, 2015, 37(3): 55–58. JIA Ze, QUAN Lin, ZHANG Shu-hong, et al. The research in shock response of warship under static explosion by underwater weapon[J]. Ship Science and Technology, 2015, 37(3): 55–58.

[ 9 ]袁建紅, 朱錫, 張振華. 水下爆炸載荷數值模擬方法[J]. 艦船科學技術, 2011, 33(9): 18–23. YUAN Jian-hong, ZHU Xi, ZHANG Zhen-hua. Numerical simulaton method study of underwater explosion load[J]. Ship Science and Technology, 2011, 33(9): 18–23.

[10]ZHANG A-man, YANG Wen-shan, YAO Xiong-liang. Numerical simulation of underwater contact explosion[J]. Applied Ocean Research, 2012, 34(1): 10–20.

[11]A. DARAMIZADEH a, M. R. ANSARI. Numerical simulation of underwater explosion near air–water free surface using a five-equation reduced model[J]. Ocean Engineering, 2015, 110(12): 25–35.

[12]張寶平, 張慶明, 黃風雷. 爆轟物理學[M]. 北京: 兵器工業(yè)出版社, 2006.

The numerical simulation of underwater explosion based on Mie-Grüneisen equation of state

WU Zong-duo1, YAN Jin1, JIANG Jie2, DONG Mei-yu3, HUANG Ji1
(1. Guangdong Ocean University School, College of Ocean Engineering, Zhanjiang 524088, China; 2. China Ship Development and Design Institute, Shanghai 201108, China; 3. Shanghai Zhenghua Heavy Industries Co., Ltd, Shanghai 200125, China)

Due to the complicated form of Mie-Grüneisen equation of state, it is difficult to handle the interface well. We separate several variables from the Euler equation and introduce mass fractions, manage to achieve the numerical simulation of multi-medium compressible flows under Mie-Grüneisen equation of state, while the process of calculation is simplified here. The approach was verified well in our case. Then under the polynomial equation of state, which is a general form of Mie-Grüneisen equation of state, we used this approach to simulate the underwater explosion of a sphere bomb, and obtained the near-field character of gas-water interaction. By introducing the mass fraction, the three states of different region, which include gas, water and gas-water mixture, are revealed clearly.

underwater；numerical calculation；Mie-Grüneisen equation of state；multi-medium

O382.1

A

1672 – 7649(2017)07 – 0029 – 05

10.3404/j.issn.1672 – 7649.2017.07.006

2016 – 07 – 12；

2016 – 11 – 10

廣東青年創(chuàng)新人才類資助項目(2014KQNCX086)；漁船漁港設施裝備標準化與安全保障能力研究資助項目(GDOU2016050258)

吳宗鐸(1984 – )，男，博士，講師，主要從事激波間斷及水下爆炸沖擊波數值模擬。