水下爆炸柱形高压气泡膨胀过程的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.