基于Level Set方法的水中气泡上升过程数值模拟

田　辉，房　媛，王文成 ，邹克武，叶阳辉

(1.承德石油高等专科学校 机械工程系，承德　067000；2.承德石油高等专科学校 学生处，承德　067000；3.西安交通大学 能源与动力工程学院，西安　710049)



基于Level Set方法的水中气泡上升过程数值模拟

田辉1，房媛2，王文成1，邹克武1，叶阳辉3

(1.承德石油高等专科学校 机械工程系，承德067000；2.承德石油高等专科学校 学生处，承德067000；3.西安交通大学 能源与动力工程学院，西安710049)

摘要：基于Level Set方法发展了一种高分辨率的求解气液相界面迁移特性的数值方法。通过三步Runge-Kutta Crank-Nicholson投影方法求解流场，Level Set方法捕捉气液交界面位置。在验证算法的捕捉性能的基础上，研究了气泡在水中上升过程的变形以及运动特性。分析结果显示气泡上升过程中的变形是气泡顶部/底部压差、两侧涡及表面张力综合作用的结果；气泡在接近壁面时由于两侧不对称涡及壁面排挤作用的影响，上升过程不断向中心区域靠近并伴随有逆时针旋转；气泡初始位置距壁面越近，壁面的影响越显著。

关键词：Level Set;投影法;Zalesak问题;气泡上升

气泡运动大量存在于能源、动力、石化、航空航天等各个领域，气液相界面迁移对流体系统的宏观流动、传热及传质特性有着显著的影响。掌握气泡的运动规律和变形机理，可为设计高效节能并可以持续稳定运行的两相流设备提供理论依据[1]。为此国内外学者发展了各种数值模拟方法用于研究多相流相界面迁移特性，如当前被普遍应用的：Front Tracking Method[2]、VOF[3]、Level Set[4]等系列算法。

Front Tracking Method通过拉格朗日思想跟踪表征界面位置的虚拟标记粒子来获取相界面的运动规律。其计算精度较高，一些典型算例已将此方法的求解结果作为基准解。但此方法会消耗大量的计算资源，并且在求解拓扑结构发生变化的问题时虚拟粒子的增减机制尚需完善。因而也制约了Front Tracking Method在工程问题上的应用。VOF方法通过欧拉思想定义计算网格内某相的体积含率来表征流场内相的分布，体积含率F∈(0,1)的区域即为界面位置。VOF方法可以很好的保证系统的质量守恒，但不能准确地求解当地曲率，进而不能准确描述表面张力的影响。同时VOF方法还必须根据特定的界面重构方法来获得界面的位置，限制了界面的捕捉精度。同样基于欧拉思想的Level Set方法通过定义计算域内各点到界面的距离为相函数，定义距离函数φ=0描述相界面。由于距离函数为连续函数便于进行微分、求导等操作，所以可以较准确地考虑表面张力对相界面迁移的影响，同时可以获得更锐利的相界面分布。但Level Set方法随着计算的进行存在质量亏损问题，Sussman[5]提出了一种距离函数重新初始化方法在获得锐利的界面分布的同时有效地保证系统质量守恒。陈斌[6-8]基于Front Tracking Method研究了壁面对气泡上升的影响，计算结果显示壁面对气泡的作用主要体现在对气泡周围流场的抑制上，由此而形成的非对称流场导致气泡逐渐偏离壁面。然而文中并未给出气泡由静止达到稳定状态的过程内气泡变形与其运动轨迹的相互影响。

本文通过自编程序实现了基于投影法及Level Set方法对气液两相界面迁移特性的捕捉。在介绍数值方法的基础上，通过求解Zalesak问题验证算法在求解曲率变化较大问题中的准确性、可行性。在此基础上对水中气泡自由上升过程的变形及运动规律进行了数值研究。

1数值方法

1.1控制方程

描述不可压缩气液两相流问题的无量纲Navier-Stokes方程可表示为：

(1)

·U=0

(2)

通过Level Set方法捕捉气液交界面。相函数φ被定义为计算域内各点到界面的距离。φ>0表示液体区域，φ<0表示气体区域，而φ=0给出了气液相界面的位置。Osher[4]等推导了速度场U作用下函数φ随时间变化所满足的相函数附加方程(3)，为Level Set方法的提出奠定了基础。

(3)

采用文献[5]中Sussman提出的距离函数重新初始化方法保证系统质量守恒。重新初始化方程表示为：

φt=signε(φ0)(1-|φ|)

(4)

(5)

(6)

1.2求解方法

采用三步Runge-Kutta Crank-Nicholson投影方法求解Navier-Stokes方程，在时间上保证二阶精度。空间上扩散相采用二阶中心差分格式离散，对流相采用二阶迎风格式离散。为增加计算稳定性采用积分平均法求解相函数附加方程(3)，五阶WENO格式离散相函数重新初始化方程(4)。离散方程可表示为：

(7)

(8)

(9)

(10)

(11)

φm+1=φm-Δt(αm+1Um·φm+βm+1Um-1·φm-1)

(12)

其中：

(13)

(14)

2算法验证

以Zalesak[9]问题为测试算例对本文发展的界面捕捉算法进行有效性验证。计算域为1×1的正方形区域，四个边界对应为周期性边界条件。计算域被一开了矩形槽的圆划分为内、外两部分。圆心位于(0.5，0.75)，半径为0.15，矩形槽长0.25，宽0.05，如图1所示。给定速度场为：

(15)

本文沿用参考文献[9]中对计算误差的定义：

(16)

其中：φe(t)、φc(t)分别为t时刻精确解及数值解，L为界面长度。

t=6.28时刻带缺口的圆形沿计算域中心旋转完一周。表1给出了本文算法与文献[9]中结果的误差比较。表中结果显示本文发展算法在较粗网格内计算精度稍差，随着网格加密计算精度与ELVIRA方法相当。图2给出了四个典型时刻的数值解及初始位置精确解。可见旋转一周后数值解除在尖角部分与精确解有所偏差外，其余部分吻合良好。

表1　测试问题计算误差E(t=6.28)

网格CLSVOFELVIRA本文算法100×1000.005720.005670.00591200×2000.002520.002620.00259400×4000.001060.001210.00117

3气泡上升过程的数值模拟

本部分以静止水中上升的三维气泡为研究对象。记气泡初始半径为R0，计算域为11R0×11R0×22R0。计算域上边界设置自由面边界条件，∂w/∂z=0；其余五个边界设为固壁，边界无滑移条件。计算采用多重网格技术，四重均分网格，网格数为41×41×81。

3.1上升气泡的稳定过程

图3(a)中每隔单位时间给出气泡的变形情况，展示了气泡在水中自由上升的过程。可以看出最初阶段上升速度、变形均较小，由于惯性造成气泡顶部与底部速度差异使得气泡有被拉长的趋势。但随着气泡与周围水的相对运动加剧，在气泡两侧各形成一个涡。其表现为推动气泡底部加速上升并使气泡不断向两侧伸展。当上升到一定阶段，气泡顶部与底部的速度差、两侧涡的强度对气泡的作用与气泡表面张力作用相平衡时，气泡形状将不再变化保持稳定上升。图3(b)给出了气泡达到稳定时两侧涡的分布情况。由于周围水的惯性较大，使气泡在上升的启动过程以及达到稳定上升前受到较大的阻力。当浮力与阻力达到平衡时气泡便以恒定的速度上升。上升过程中气泡平均速度变化情况如图4所示。

3.2壁面对气泡运动的影响

通过考察初始位置距壁面不同距离的气泡上升过程来分析壁面对气泡上升过程的影响。图5中所示分别为初始位置距壁面3.5R0及2R0的气泡自由上升过程。图中可见气泡在上升过程不断向中心区域靠近。并在此过程中伴随有逆时针旋转，且气泡初始位置越靠近壁面旋转越明显。

如图6(a)所示，上升初期气泡左侧(远离壁面的一侧)已经形成了较明显的涡，而右侧涡强度较小。受左侧逆时针旋转涡的影响气泡则呈现出明显的旋转特性。同时在壁面排挤作用影响下气泡逐渐向区域中心运动。如图6(b)所示，随着气泡的上升右侧的涡并没有发展成熟，而表现为明显的壁面排挤作用。使气泡及左侧的涡向区域中心运动。此时由于气泡与左侧涡中心相对位置的变化，旋转作用有所减弱。

4结论

本文发展了一种基于Level Set Method求解气液相界面迁移特性的数值方法。通过求解Zalesak问题验证了算法在捕捉曲率变化较大问题中的准确性、可行性。在此基础上以水中的气泡为研究对象，模拟了气泡初始位置在计算域中央以及靠近壁面情况下的上升过程。结果表明：1)气泡受浮力驱动向上运动，在其顶部与底部间的压差、两侧涡及表面张力共同作用下发生变形，当以上各种作用达到平衡时气泡达到稳定上升状态；2)靠近壁面的气泡受到其两侧不等强度涡的影响及壁面排挤作用的影响，在上升的过程中伴随着逆时针的旋转及向中心区域的运动。气泡初始位置距壁面越近，壁面的影响越显著。

参考文献：

[1]田辉,房媛,王文成，等. 离心泵内变工况流动特性的数值研究[J]. 承德石油高等专科学校学报, 2015(2):28-31.

[2]Tryggvason, G. A front-tracking method for the computations of multiphase flow[J]. Journal of Computational Physics, 2001. 169(2): 708-759.

[3]Rider, W.J. and D.B. Kothe. Reconstructing volume tracking[J]. Journal of Computational Physics, 1998. 141(2): 112-152.

[4]Osher, S. and J.A. Sethian. Fronts Propagating with Curvature-Dependent Speed-Algorithms Based on Hamilton-Jacobi Formulations[J]. Journal of Computational Physics, 1988. 79(1): 12-49.

[5]Sussman, P. Smereka, S. Osher. A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[J]. Journal of Computational Physics, 1994. 114(1): 146-159.

[6]陈斌, T.Kawamura, Y.Kodama. 静止水中单个上升气泡的直接数值模拟[J]. 工程热物理学报, 2005(6): 82-84.

[7]陈斌. 倾斜壁面附近上升气泡的直接数值模拟[J]. 工程热物理学报, 2007(6): 965-967.

[8]陈斌. 高粘度流体中上升气泡的直接数值模拟[J]. 工程热物理学报, 2006,02: 255-258.

[9]Sussman, E.G. Puckett. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J]. Journal of Computational Physics, 2000. 162(2): 301-337.

Numerical Simulation of Bubble Rising in Water Based on Level Set Method

TIAN Hui1, FANG Yuan2, WANG Wen-cheng1, ZOU Ke-wu1, YE Yang-hui3

(1.Department of Mechanical Engineering, Chengde Petroleum College,Chengde 067000, Hebei, China;2.Department of Student Affairs, Chengde Petroleum College, Chengde 067000, Hebei, China;3.School of Energy and Power Engineering, Xi’an Jiaotong University,Xi’an 710049, Shaanxi, China)

Abstract:Based on level set method, a high resolution numerical method was developed to study the migration characteristics of gas-liquid interface. Three step Runge-Kutta Crank-Nicholson projection method was used to solve the Navier-Stokes equations. Level set method was applied to capturing the phase interface between air and water. Zalesak’s problem was solved to validate the numerical method. The numerical results of bubble rising in water show that the deformation of bubble during the rising process is subject to (1) the pressure difference between the top and bottom of bubble, (2) the vortex, and (3) surface tension. When close to the wall, due to the wall repulsion and the asymmetric vortex, the deformed bubble moves to the centerline with clockwise rotation. The closer to the wall, the wall effect on the bubble will be more significant.

Key words:level set; projection method; Zalesak’s problem; bubble rising

基金项目：国家自然科学基金资助项目(两相非定常空化空泡时空演化机理及其对高效离心泵水动力性能影响的研究)：51076126

收稿日期：2015-10-30

作者简介：田辉(1982-)，男，河北承德人，承德石油高等专科学校机械工程系讲师，博士，主要从事气液两相界面迁移特性数值研究。

中图分类号：O359

文献标识码：A

文章编号：1008-9446(2016)02-0028-05