基于VOF法的近自由面水下爆炸气泡运动数值模拟

李健， 黄红生， 林贤坤， 荣吉利， 项大林

(1.广西科技大学 广西汽车零部件与整车技术重点实验室， 广西，柳州 545006；2.北京理工大学 宇航学院， 北京 100081)



基于VOF法的近自由面水下爆炸气泡运动数值模拟

李健1， 黄红生1， 林贤坤1， 荣吉利2， 项大林2

(1.广西科技大学 广西汽车零部件与整车技术重点实验室， 广西，柳州 545006；2.北京理工大学 宇航学院， 北京 100081)

对于近水面水下爆炸气泡的运动特性，在结合SOLA算法差分格式，利用投影算法对气泡的控制方程进行数值求解的基础上，基于VOF法进行了数值模拟，并利用MSC.DYTRAN软件进行仿真分析，通过对两者的分析结果对比研究. 结果表明，基于结合SOLA算法差分格式的投影算法，利用VOF法可对水下爆炸气泡各时刻的运动界面进行较精确的捕捉和重构，从而实现精确追踪水下爆炸气泡脉动全物理过程的运动界面. 同时，也研究了距离参数对气泡坍塌、射流及环状气泡的影响，掌握了距离参数对近水面水下爆炸气泡运动特性的影响.

投影算法；SOLA算法；气泡；VOF法；水下爆炸；界面追踪

流体中脉动的流体同相振荡时彼此吸引，反相振荡时相互排斥，这种在流体中脉动物体间的相互作用力称为Bjerknes力. 在自由面的Bjerknes力和浮力两者作用下，近水面水下爆炸的气泡运动是一个复杂的能量释放与传递过程，包括膨胀、收缩、变形及气泡向自由表面漂移等一系列运动，可能形成向上、向下或相对射流，射流击穿初始球形气泡而形成环形气泡[1]. 随着气泡脉动的进行，自由面还会产生连续的水冢，因此气泡在近水面附近的运动状态会相当复杂. 为了揭示近水面气泡运动的规律，需实时掌握气泡运动、漂移、坍塌以及射流特性. 因此，追踪各时刻气泡的运动界面，研究近水面水下爆炸气泡的运动特性，对掌握弹药毁伤威力及气泡运动对目标的损伤具有重要的意义.

在水下爆炸气泡数值仿真研究方面，目前国内外主要借助成熟商业软件[2-4]或边界积分法(BIM)[5-6]. 商业软件往往涉及的面较宽泛但不够专业，如MSC.DYTRAN软件，虽可实现两相流的数值计算，但因其中的Euler算法的基础是FVM法，而无法直接捕捉运动界面. BIM法虽可较准确描述气泡边界变形，但其本质上还是属于Lagrange方法，在处理气泡严重变形、破碎、融合方面缺点较明显，精度难以保证，特别是当水下爆炸气泡由于连续脉动出现的气泡多次破碎、融合现象时，BIM法的局限性更加明显. 由Hirt C W等[7]首先提出的VOF法，可灵活处理介质面、自由面、间断面或各种内部的运动界面. 特别以其存储小、应用简单，追踪的界面锐利性好、精细，可处理自由面重入等强非线性现象，而深受广大CFD工作者的喜爱[8]. 随着VOF法不断完善，目前已有多种成熟的界面重构算法，如FLAIR和Youngs界面重构技术. Ashgriz N等[9]通过对任意运动界面的两个相邻网格，构造一条斜线段，作为跨过该边界网格的界面近似，该方法只适用于均匀网格；Seifollahi M等[10]采用基于PLIC重构技术的VOF方法，考虑表面张力情况下模拟了气泡水中上浮过程. 因此，对于水下爆炸气泡运动这类问题的数值计算，VOF法结合界面重构技术较其它方法具有较明显的优势.

为了对近水面水下爆炸气泡脉动全物理过程进行数值模拟，研究气泡的运动特性，本文首先基于N-S(Navier-Stokes)方程建立水下爆炸二维气泡的数学模型，结合SOLA算法差分格式，利用投影算法对气泡的控制方程进行数值求解，并利用VOF法中的Youngs界面重构技术重构各时刻气泡的运动界面. 通过与商业软件的对比，验证了数值计算方法与程序设计的正确性，以此为基础，基于VOF方法研究近自由面水下爆炸气泡的运动特性.

1 气泡运动的数学模型

1.1 控制方程

对于水下爆炸气泡运动问题，实际上是气液两相流问题，气泡的运动由控制方程决定，它主要包括连续性方程和动量方程. 由于水可认为是不可压缩流体，若不考虑黏性和假设液体的密度为常数，则在笛卡儿坐标系下，二维气泡运动的数学模型可表示为如下的N-S方程组：

(1)

式中：ρ为密度；u，v为x，y方向的流速；t为时间；p为压力；Hx，Hy为x，y方向单位质量流体所受外部体力，包括重力和表面张力的作用.

1.2 运输方程

设在二维计算区域Ω中，流体A所在的区域为Ω1，流体B所在的区域记为Ω2. 采用正交网格对计算区域Ω进行网格划分，对于网格，定义

(2)

各网格的C值构成VOF函数. 对于两种不相溶的流体组成的流场，设(u，v)是流体速度向量，则VOF函数C满足则C(x，y，t)满足

(3)

2 控制方程的数值解法

本文采用投影算法结合SOLA算法差分格式求解控制方程. 式(1)的动量方程可表示为

(4)

式中A(vn)为动量方程中的对流项. 投影算法分两步计算第n+1步的速度vn+1和压力pn+1，即

(5)

(6)

因vn+1需满足散度为0，即满足连续方程·vn+1=0. 因此，对等式(6)两边取散度，有

(7)

因此，由式(5)可得中间速度v*，再由式(7)可得压力pn+1，最后由式(6)可得速度vn+1.

2.1 网格边界中心速度的中间值求解

对于网格(i，j)，为了保证计算稳定和易于收敛，对式(1)的连续方程采用如下隐式差分格式,

(8)

利用SOLA算法差分格式，对方程(5)进行离散，有

(9)

(10)

2.2 网格中心的压力求解

对式(7)的压力方程，采用5点差分格式的超松驰迭代法(SOR)求解，其迭代形式为

(11)

(12)

3 VOF界面追踪的理论与方法

3.1 运输方程的数值解法

对方程(3)沿一个单元网格空间和单元时间积分

(13)

式中Δt为单位时间间隔. 对方程(11)积分后有

(14)

方程左端的3项分别为在时间内通过网格单元内流量的总变化量、垂直边的变化量和水平边的变化量. 其中ΔCi±1/2，j和ΔCi，j±1/2的积分形式为

(15)

(16)

3.2 交界面的单位法向量计算

采用差分方法，通过网格(i，j)附近的8个网格，如图1所示，组成一个求解中心交界面网格中的界面法向的差分格式.

(17)

(18)

(19)

(20)

(21)

4 仿真实例

4.1 对比研究

为了验证本文算法在近水面气泡运动数值计算的准确性，基于文献[2]中的计算工况，利用VOF法结合Youngs界面重构技术，对水域为12 m×8 m，装药量为3.0 kg的爆炸在水下2.0 m处爆炸气泡与自由面水冢运动全过程进行数值计算，并将计算结果与MSC.DYTRAN软件计算结果进行对比，原因是文献[2]已经验证了仿真分析结果的正确性. 气泡的状态方程为

(22)

p=pc+p0V0Vλ， p0=1.39×105WV0λ,

式中：p为气泡压力；pc为可冷凝气体的饱和蒸汽压，量级与大气压相当，一般可忽略；p0和V0分别为气泡形成时的初始压力和初始体积，W为装药量，λ为气体的比热比，对于TNT水中爆炸的爆轰产物而言，λ取1.25.

空气的状态方程采用理想气体状态方程：

(23)

式中：p为气泡压力；e为比热力学能；γ为常数；ρ为材料密度；对于空气，取ρ=1.184 8 kg/m3，γ=1.4，e=2.14×10-7m2/s2.

针对气泡与自由面运动过程，基于MSC.DYTRAN进行仿真分析，结果如图3所示，采用本文所提算法进行数值计算，结果如图4所示.

其中图3(a)、图4(a)为气泡初始状态，图3(b)、图4(b)为气泡第一膨胀至体积最大状态. 根据图3(c)、3(d)与图4(c)、4(d)，可以明显看到在Bjerknes力作用下气泡的坍塌并形成射流. 从图3(e)、3(f)与图4(e)、4(f)可以看出，在射流作用下，气泡被击穿，并以此继续膨胀，通过对比可看出，VOF运动界面追踪法与MSC.DYTRAN的分析结果吻合得较好. 由此可见，本文所提出算法可很好地模拟气泡膨胀、脉动、坍塌、破碎以及自由面连续水冢等一系列复杂的物理现象.

4.2 距离参数的影响

气泡坍塌、射流以及环状气泡等物理现象的产生与爆心初始位置有密切关系，因此需研究距离参数对气泡坍塌、射流及环状气泡的影响. 以气泡脉动过程中的最大半径Rm为特征长度、爆心处静水压力为特征压力，将各参数量纲一化，则特征时间为t=Rm(ρ/pamb)0.5. 以h=H/Rm量纲一的参数表示爆心初始位置和自由面之间的距离H与气泡最大半径Rm之间的比值，研究球形TNT炸药在h变化时气泡及水冢运动的特性. 计算结果如图5所示，以气泡对称轴上剖视图表示.

由图5(a)可知，由于h=0.81较小，气泡受自由面Bjerknes力影响明显，气泡在坍塌阶段形成背离自由面的凹型射流，射流穿透气泡致使其形成环状气泡；当h=0.98时，自由面Bjerknes力效应有所减弱，但仍占优势，气泡产生对射流，但向下方向射流较强；当h=1.16时，气泡形成对射流，但浮力逐渐占优势，最终形成指向自由面的射流；当h=1.35时，自由面Bjerknes力效应逐渐减弱，浮力作用效果明显，形成明显的向上射流.

图6为相同装药量在不同水深处气泡体积的时程曲线图. 从图中可以看到，气泡体积会随着水深的增加而减小. 图7为气泡壁上端点与下端点随时间的变化规律.

从图中可以看到，当气泡与自由面距离较近时，在气泡收缩阶段，其上端点会产生背离自由面方向的速度，下端点会产生指向自由面的速度，但上端点的速度要远高于下端点的速度，最终气泡会被射流击穿形成环状气泡，且击穿点的位置远低于气泡初始位置. 随着气泡与自由面距离的增加，其被射流击穿的位置会逐渐上移，表现为气泡上端点的速度会随着气泡与自由面距离的增加而减小，气泡下端的速度会随着气泡与自由面距离的增加而增加，其主要原因就是Bjerknes力会随着特征距离的增加而减弱，特别是当特征距离h>3.0，即当爆心距自由面距离超过3倍气泡最大半径时，Bjerknes力对气泡运动的影响非常小，浮力占绝对优势，此时不仅气泡会形成明显的指向自由面的射流，且自由面基本无水冢产生，可以不考虑自由界面的影响.

5 结 论

基于N-S方程建立了近自由面水下爆炸二维气泡运动的数学模型，采用SOLA算法差分格式，利用投影算法对气泡的控制方程进行了数值求解，基于VOF法结合Youngs界面重构技术对近水面水下爆炸气泡脉动全物理过程进行了数值模拟. 通过计算得到如下一些主要结论：

① 采用投影算法结合SOLA算法差分格式对控制方程进行求解，可较准确获得水下爆炸气泡各时刻的速度和压力，为采用VOF法的相关理论与方法精确追踪气泡的运动界面奠定了良好的基础.

② 采用VOF法结合Youngs界面重构技术，可对水下爆炸气泡各时刻的运动界面进行较精确的捕捉和重构，从而实现精确追踪水下爆炸气泡脉动全物理过程的运动界面.

③ 近水面水下爆炸中，特征距离对气泡运动影响明显. 当爆心与自由面的初始距离小于气泡膨胀最大半径时，自由面水柱高而细长，自由面水柱高度会随着特征距离的增加而减小，随着爆心距自由面距离的增加，当特征距离h>3.0时，自由面Bjerknes作用力对气泡运动的影响基本可忽略，此时气泡形成指向自由面的射流，且自由面无水冢产生.

④ 由于Bjerknes作用力的存在，近自由面水下爆炸气泡上端受自由面影响，会产生背离自由面的射流，而气泡下端在浮力的作用下产生指向自由面的射流，且上端点的射流速度会随着特征距离的增加而减小，下端点的射流速度会随着特征距离的增加而增加.

[1] 李健，荣吉利，项大林.近自由面水下爆炸气泡运动的数值计算研究[J].工程力学，2011，27(6):200-205．

Li Jian， Rong Jili， Xiang Dalin. Numerical study of bubble motion by underwater explosion near free surface[J]. Engineering Mechanics， 2011，27(6):200-205. (in Chinese)

[2] Li Jian， Rong Jili. Bubble and free surface dynamics in shallow underwater explosion[J]. Ocean Engineering， 2011，38:1861-1868．

[3] Rajendran R. Numerical simulation of underwater explosion bulge test[J]. Materials and Design， 2009，30:4335-4341．

[4] 刘雨，黄风雷，吕中杰.蒸汽爆炸过程中压力波冲击气泡现象的数值模拟[J].北京理工大学学报，2013，33(4):339-342．

Liu Yu， Huang Fenglei， Lü Zhongjie. Numerical simulation of shock wave impacting bubble in vapor explosion[J]. Transactions of Beijing Institute of Technology， 2013，33(4):339-342. (in Chinese)

[5] Zhang A M， Yao X L， Yu X B. The dynamics of three-dimensional underwater explosion bubble[J]. Journal of Sound and Vibration， 2008，311:1196-1212.

[6] Wang C， Khoo B C， Yeo K S. Elastic mesh technique for 3D BIM simulation with an application to underwater explosion bubble dynamics[J]. Computers & Fluids， Volume， 2003，32(9):1195-1212.

[7] Hirt C W， Nichols B D. Volume of fluid method for the dynamics of free boundary [J]. Comput Phys， 1981，93:201-225．

[8] 赵大勇，李维仲.VOF方法中几种界面重构技术的比较[J].热科学与技术，2003，2(4):318-323．

Zhao Dayong， Li Weizhong. Comprison of three kinds of interface reconstruction methods applied to VOF[J]. Journal of Thermal Science and Technology, 2003，2(4):318-323. (in Chinese)

[9] Ashgriz N， Poo J Y. FLAIR: flux line-segment model for advection and interface reconstruction[J]. J Comput Phys， 1991，93(2):449-468.

[10] Seifollahi M， Shirani E， Ashgriz N. An improved method for calculation of interface pressure force in PLIC-VOF methods[J]. Eruopean Journal Mechanics B/Fluids， 2008，27:1-23．

[11] 万德成.用VOF法数值模拟孤立波迎撞及在后台阶上的演化[J].计算物理，1998，15(6):655-666．

Wan Decheng. Numerical simulations of the collisions between two solitary waves and the transformation of a solitary wave over backward step by VOF method[J]. Chinese Journal of Computational Physics， 1998，15(6):655-666. (in Chinese)

[12] 邹志利，邱大洪，王永学.VOF方法模拟波浪槽中二维非线性波[J].水动力学研究与进展:A辑，1996(1):93-103．

Zou Zhili， Qiu Dahong， Wang Yongxue. Numerical simulation of nonlinear wave generated in wave flume by VOF technique[J]. Journal of Hydrodynamics: Serial A， 1996(1):93-103. (in Chinese)

[13] Eugen F F Botta， Marcel H M Ellenbroek. A modified SOR method for the poisson equation in unsteady free-surface flow calculations[J]. Journal of Computational Physics， 1985，60(1):119-134.

(责任编辑：刘雨)

Numerical Study on Bubble Motion by Underwater Explosion Near Free Surface Based on VOF

LI Jian1， HUANG Hong-sheng1， LIN Xian-kun1， RONG Ji-li2， XIANG Da-lin2

(1.Guangxi Key Laboratory of Automobile Components and Vehicle Technology， Guangxi University of Science and Technology， Liuzhou， Guangxi 545006， China; 2.School of Aerospace Engineering，Beijing Institute of Technology， Beijing 100081， China)

The difference format of SOLA and projection algorithm were introduced to solve the control equation of bubble motion in underwater explosion. The dynamic behavior of bubble near free surface was simulated and analyzed with VOF and MSC.DYTRAN software. The computational results show that the VOF method can accurately track and reconstruct the motion of bubbles based on projection algorithm combined with the SOLA difference format. The interactions between bubble and free surface were studied systematically and summarized the dynamic behaviors of bubbles, including the influence of distance parameter on bubble collapse, eject and annularity bubble, that have close relation with distance parameter.

projection algorithm; SOLA algorithm; bubble; volume of fluid(VOF); underwater explosion; interface tracking

2014-07-13

国家自然科学基金资助项目(51209042，11272057)

李健(1980—)，男，博士，教授，E-mail:lijian0772@126.com.

林贤坤(1976—)，男，博士，副教授，E-mail:linxk0209@yahoo.cn.

O 352

A

1001-0645(2016)02-0122-06

10.15918/j.tbit1001-0645.2016.02.003