史汝超,張亞軍,徐勝利

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

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

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

1 控制方程和數(shù)值方法

1.1 Level set方法

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

(1)

(2)

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

1.2 RGFM

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

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

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

1.3 控制方程及狀態(tài)方程

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

1.4 數(shù)值格式

2 計算結(jié)果與分析

2.1 算例1

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

圖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

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

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

3 結(jié) 論

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

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