水下爆炸柱形高壓氣泡膨脹過程的RGFM和高精度格式數值模擬*

史汝超,張亞軍,徐勝利

(1.中國科學技術大學近代力學系，安徽 合肥 230026; 2.中國工程物理研究院流體物理研究所，四川 綿陽 621999)

柱形高壓氣泡膨脹問題常見于柱形裝藥水下爆炸。對此類問題進行數值模擬，精確處理和追蹤高壓力比、高密度比的氣水界面一直是難點。Rayleigh-Plesset運動方程[1]為眾多氣泡運動理論研究提供了參考；J.R.Blake等[2]和A.Pearson等[3]采用邊界元方法研究了水下爆炸氣泡演化過程；T.A.Vernon[4]和Z.Zong[5]采用拉普拉斯方程描述速度場，研究高壓氣泡運動。

而對圓柱形高壓氣泡膨脹問題進行三維高精度數值模擬的文獻較少，是因為柱形高壓氣泡的三維精確建模和數值模擬較復雜。如果使用level set重新初始化技術對整個流場的level set初值直接進行定義，會造成界面的非物理移動[6]。因此，本文中，對流場的level set初值進行分區定義，以實現精確建模，同時用RGFM(real ghost fluid method)[7]處理氣水界面附近網格節點，并結合高精度格式，對柱形高壓氣泡膨脹問題進行高精度數值模擬。

1 控制方程和數值方法

1.1 Level set方法

使用level set方程[8]追蹤氣水界面：

(1)

(2)

若柱形氣泡傾斜放置，則將坐標系旋轉，使柱形氣泡中軸線與旋轉坐標系的z軸平行。因為坐標系旋轉不會改變φ，所以在旋轉坐標系中，仍可使用式(2)定義φ的初值。對于多爆源，用式(2)分別求解φ1、φ2、…、φn。在計算中，每隔一定的時間步，還需要對φ重新初始化[9]。

1.2 RGFM

如圖2所示：(1)先用ARPS(approximate Riemann problem solver)[10]求解界面參數pI、uI、ρIL、ρIR；(2)將界面參數pI、uI、ρIL賦給真實節點i；(3)將pI、uI、ρIL賦給虛擬節點i+1；(4)在虛擬區域中，將節點i+1的值沿界面法線方向外推，具體外推方程見文獻[11]。在三維計算中，用界面兩側法向偏轉角度最小的兩個節點代入一維模型求解界面參數。

圖1 柱形氣泡φ初值定義示意圖Fig.1 Schematics of defining initial φ of cylindrical bubble

圖2 RGFM氣水界面附近節點賦值示意圖Fig.2 Schematics of updating the nodes next to interface using RGFM

1.3 控制方程及狀態方程

忽略黏性項和重力項，用Euler方程描述流場；使用瞬時爆轟模型，描述水下爆炸的高壓氣體膨脹過程和沖擊波傳播過程，高壓氣體的初始密度為裝藥密度，初始壓力為給定值，并用理想氣體狀態方程描述；對于水，用Tait方程[12]描述。

1.4 數值格式

2 計算結果與分析

2.1 算例1

柱形高壓氣體的初始參數為:pg=21 GPa，ρg=1.6 t/m3，γg=3.0。水的初始參數為：pl=0.1 MPa，ρl=1.0 t/m3。如圖3(a)所示，計算區域取[0,30]×[0,30]×[0,30]。炸藥的中心坐標為(15,15,10)。柱形氣泡長10.0，底面直徑6.0,與面DCGH呈夾角45°。邊界條件取:計算域底部為無滑移壁面，其余5個面為流出的。網格數量為121×121×121。由圖4(a)，柱形高壓氣泡開始膨脹時，首先在水中產生一道激波，同時氣泡內部產生一道稀疏波，并向氣泡中心傳播。由于氣泡在底面方向和徑向的膨脹，兩個底面附近均產生一圈低壓區。由圖4(b)，激波到達固壁，并開始反射。由圖4(c)，隨著計算時間的推進，反射激波開始打向氣泡，抑制氣泡膨脹。

圖3 高壓氣泡初始位置示意圖Fig.3 Schematics of initial location of high pressure bubble

圖4 算例1流場壓力云圖Fig.4 Pressure cloud pictures of flow field of case 1

2.2 算例2

圖5 算例2流場壓力云圖Fig.5 Pressure cloud pictures of flow field of case 2

雙圓柱形裝藥，爆炸氣體產物的初始參數為：pg=21 GPa，ρg=1.6 t/m3，γg=3.0。水的初始參數為：pl=0.1 MPa，ρl=1.0 t/m3。如圖3(b)所示，計算域取[0,30]×[0,10]×[0,30]，炸藥形狀為兩個相同尺寸的圓柱體，高6.0，底面直徑3.0，同計算域底面DCGH平行。炸藥的中心坐標分別為(9,5,12)、(21,5,18)。邊界條件取面ABCD和面GHEF為無滑移壁面，其余4個面為流出的。面GHEF上點1、2的坐標分別為(6,0,12)、(21,0,18)。網格數量為151×51×151。由圖5(a)，柱形氣泡首先膨脹，在水中產生一道壓縮波，在氣泡內部產生一道稀疏波，兩個截面分別為y=5和x=21。由圖5(b)，大約在t=0.90 ms時，壓縮波到達壁面，并產生反射波，兩個截面分別為x=9和x=21。由圖5(c)，大約在t=1.10ms時，固壁的反射波與氣泡發生作用，開始抑制氣泡在徑向上的膨脹；同時兩個氣泡膨脹產生的激波發生疊加，彼此抑制對方的膨脹。由圖6(a)，激波在到達固壁后再反射的這個過程中，對固壁的作用區域為橢圓形。由圖6(b)，柱形高壓氣泡在膨脹過程中，形狀從圓柱逐漸向橢球變化。由圖6(c)，兩個指定點的壓力峰值分別為約4.11和6.56 MPa。

圖6 算例2的計算結果Fig.6 The calculation results of case 2

3 結 論

對柱形裝藥水下爆炸高壓氣泡膨脹過程進行了高精度數值模擬，給出了不同時刻壓力云圖、氣泡形狀的變化趨勢以及指定點的壓力曲線。在對流場的level set初值精確定義的前提下，用RGFM結合高精度格式(五階WENO、四階R-K法)可以很好地處理柱形高壓氣泡膨脹問題以及追蹤高密度比、高壓力比的氣水界面。計算結果表明，柱形高壓氣泡在膨脹過程中，其形狀逐漸向橢球形變化。受固壁反射波的影響，高壓氣泡在固壁法線方向上的膨脹受到抑制；雙圓柱形高壓氣泡膨脹產生的激波，可以彼此抑制另一個氣泡的膨脹。

[1] Plesset M S. The dynamics of cavitation bubbles[J]. Journal of Applied Mechanics, 1949,16(16):277-282.

[2] Blake J R, Gibson D C. Growth and collapse of a vapour cavity near a free surface[J]. Journal of Fluid Mechanics, 1981,111:123-140.

[3] Pearson A, Cox E, Blake J R, et al. Bubble interactions near a free surface[J]. Engineering Analysis with Boundary Elements, 2004,28(4):295-313.

[4] Vernon T A. Whipping response of ship hulls from underwater explosion bubble loading[R]. Dartmouth, Nova Scotia: Defence Research Estabishment, 1986.

[5] Zong Z. A hydroplastic analysis of a free-free beam floating on water subjected to an underwater[J]. Journal of Fluids and Structures, 2005,20(3):359-372.

[6] Russo G, Smereka P. A remark on computing distance functions[J]. Journal of Computational Physics, 2000,163(1):51-67.

[7] Wang C W, Liu T G, Khoo B C. A real ghost fluid method for the simulation of multimedium compressible flow[J]. SIAM Journal on Scientific Computing, 2006,28(1):278-232.

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

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

[10] 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.

[11] Aslam T. A partial differential equation approach to multidimensional extrapolation[J]. Journal of Computational Physics, 2004,193(1):349-355.

[12] Wardlaw Jr A B. Underwater explosion test cases[R]. Indian Head Division: Naval Surface Warfare Center, 1998.

[13] Shu C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[C]∥Advanced numerical approximation of nonlinear hyperbolic equations. Springer, 1997:325-432.

[14] Shu C W, Osher S. Efficient implementation of essentially non-oscillatory shock capturing schemes[J]. Journal of Computational Physics, 1988,77(2):439-471.

[15] Jiang G S, Peng D. Weighted ENO schemes for Hamilton-Jacobi equations[J]. SIAM Journal on Scientific Computing, 2000,21(6):2126-2143.