模拟气泡上升运动的一种改进的Level Set方法

张一夫，苑宇，李维重

(1.辽宁工程技术大学 安全科学与工程学院，辽宁 葫芦岛 125105;2.大连交通大学 机械工程学院，辽宁 大连 116020;3，大连理工大学 能源与动力学院，辽宁 大连 116023)

0 引言

气泡现象普遍存在于自然界和工程实际中［1-2］，广泛用于环境、能源、化工及水利工程等领域，如液－气燃料燃烧问题中的液－气界面、水轮机和水泵的空化空蚀、船舶螺旋桨水流、气液化学反应、废水处理、破碎的波浪在水面下产生的大量气泡等.气泡在液体中的上升运动过程是一个非线性、复杂、不稳定的动力过程.气泡在变形过程中还常常伴随着融合、破碎、翻转、摆动等一系列的复杂运动，并且气泡的运动、变形以及尾涡还会对周围流场的速度、压力场产生影响.因此，对竖直流道内气泡上升过程的研究是一项具有挑战性和现实意义的研究课题.

目前，常用的追踪气泡的方法主要有:VOF、Level Set、Front Tracking、Phase Field 等方法［3-8］.相比较而言，Level Set方法不仅能够比较准确地追踪运动界面，而且无须进行复杂繁琐的界面重构技术，容易编程，具有较大的通用性，另一方面，可以采用高精度高分辨率的差分格式，对物质界面定位更加准确.但由于Level Set方法是一个非守恒的算法，将Level Set方法应用于不可压缩两相流问题的模拟中，由于数值方法引起的数值耗散问题，必然会引起气泡质量的变化(质量损失或增加)，使相界面逐渐偏离正确的位置，进而使流场失真;另一方面，考虑到Level Set方法在全场进行重新初始化时，会由于边界初网格非正交性引起的数值耗散.因此，本文作者引入体积修正法和局部重新初始化方法对临近并排双气泡的上升运动过程进行了数值模拟，成功地解决了上述问题.

1 数学模型

1.1 基本假设

本模型基于如下假设:

(1)流体模型为低雷诺数层流模型;

(2)用均一的单流体模型计算变密度、变黏度系数的不相混两相流;

(3)液相和气相都视为不可压缩流体.

1.2 控制方程组

本文采用的无量纲控制方程组包括:

(1)不可压缩流体连续性方程

(2)考虑表面张力和浮力的动量方程

(3)采用Level Set方法追踪两相界面的输运方程

式中，u为速度矢量;ρ为密度;μ为动力黏度;p为压力;B为浮力项，B=(1－ρ)/Fr;F为表面张力项，F=－ kδs(φ)n/We;n 为界面单位法向矢量，n=－▽φ/|▽φ|;k为界面的曲率k=▽·n=▽·(－▽φ/|▽φ|);δs为与界面有关的Dirac分布;φ为Level Set函数.无量纲数

1.3 界面追踪方法及两相物理量参数处理

气-液两相界面的区分是采用一个光滑的Level Set函数φ(x，t)来实现的.当光滑函数φ等于零时，该光滑函数的变量x就代表了分界面上的点，当φ值大于零为界面外部，小于零为界面内部.

为了避免在界面附近由于大密度比及大粘度系数比引起的数值不稳定，采用修正的Heaviside函数来光滑两相参数.因此，无量纲密度和黏度系数由如下公式得出

式中，λρ= ρg/ρl和λμ= μg/μl分别为密度比和黏度系数比.

光滑的Heaviside函数H(φ)定义如下

式中，ε为数值弥散界面的可调带宽.

由于数值方法的内在因素，随着计算的进行，Level Set函数往往无法保持符号距离函数(Signed Distance Function，SDF)的特性.因此，为了使φ(x，t)保持这种特性，需要采用一种方法来将Level Set函数重置为符号距离函数.利用符号距离函数的如下性质:

若φ为符号距离函数，则有:

重新初始化的偏微分方程为

式中，需将sign(φ0)光滑化为为重新初始化之前，上一时层的最后计算值.

考虑到Level Set方法在全场进行重新初始化时，会由于边界处网格非正交性引起的数值耗散，因此，将Level Set重新初始化方法局限在了界面附近一个很小的范围内，可以大大提高Level Set计算的精度.因此，本文采用类似于Peng等［9］提出的窄带计算的局部重新初始化方法，将两相界面定义在一个有限的区间(－β≤φ≤β)内.界面仍采用符号距离函数的定义方法(如图1所示)，对φ值进行重新初始化计算.

图1 界面函数初始值定义(当β=0.15时)

为保证气泡在计算过程中具有很好的质量守恒性，采用体积校正方法［10］对气泡体积进行修正.这种方法的主要思想是通过不连续相(气泡或液滴)体积的盈亏来调整界面位置:

假设气泡是球形的，半径的增量为

其中，C1为模拟体系的维数.由于分散相体积的增加(或减少)对应于从界面到分散相中心距离的增加(或减少)，所以可取Level Set函数的校正量正比于Δr.因此，Level Set函数的校正量可表达为

其中，C2是一个经验参数，一般在0.01～0.1之间的范围内取值，C2较大可能导致重新初始化方程不熟练，较小会降低计算效率.本文中C2取为0.01 .

这样，Level Set方法的体积校正方程为

1.4 表面张力项处理

表面张力项采用连续表面张力模型(Contin-uous Surface Force，CSF)［11］模型计算来处理.对于Level Set方法，CSF模型的表面张力表达式为

其中，δ函数表示如下:

2 控制方程的离散和求解

分别采用三阶精度的QUICK格式以及二阶精度的中心差分格式，对N-S方程中对流项与扩散项进行离散.主场的压力与速度耦合采用经典的SIMPLE［12-13］算法求解.物理量及参数定义采用同位网格［14］，所有物理量都定义在网格节点的中心处.界面上的逆变速度采用非交错网格上的Rhie-Chow动量插值方法［15］获得，在计算界面压力梯度时，可以有效抑制不合理的压力场产生.

Level Set输运方程采用二阶上风和带有Superbee限制算子的平均积分型TVD格式进行离散求解［16］，以保证非线性对流项数值稳定性.重新初始化方程采用修正的Godunov方法对五阶WENO差分格式进行求解，得到重新初始化值.时间导数项采用三阶的TVD-Runge-Kutta方法，以保证每一时间步Level Set方程都满足符号距离函数的特性.

采用C++语言编制代码求解控制方程.计算流程如下:

(1)初始化函数φ为到界面的符号距离函数;

(2)采用SIMPLE算法求解动力学方程，从tn时间步的速度和压力场获得tn+1的值;

(3)求解Level Set方程和局部重新初始化方程以保证φ为到界面的符号距离;

(4)采用体积校正法修正气泡体积;

(5)更新tn+1时层上的流体特性参数;

(6)时间递进一步，重复(2～5)步直至最后时间步.

3 模拟结果与讨论

本文用未采用体积校正法的全场重新初始化方法(方法1)和采用体积校正法的局部重新初始化方法(方法2)分别对竖直流道内并排双气泡的上升过程进行了数值模拟.左右边界和底部边界均采用无滑移边界条件，上端为开口边界，采用外推边界条件.初始时，气泡形状为圆形，流场静止，两气泡圆心间距为Δd=2.4R，气泡初始圆心位置距底面为 3R0.模拟工况:Re=10，We=20，Fr=1.0，ρb/ρl=1/1 000 和 μb/μl=1/100.计算网格为120×240.利用两种方法分别模拟出的气泡形态对比如图2所示.

图2 两种方法模拟出的气泡形态演变对比

从图2中可以看出，利用上述两种方法模拟出的气泡形态演变过程有明显差异.方法1模拟出的结果来看，随着气泡的上升，气泡明显变小，气泡体积守恒性很差.而方法2模拟出的结果，气泡守恒性很好.另外，两种方法模拟出的气泡上升速度也有所不同(如图3所示)，随着时间的推移，方法1模拟出的结果中，气泡在经过开始的加速阶段后，速度达到一个峰值，但由于气泡体积的逐渐减小，气泡的速度开始下降.在整个过程中，气泡始终都达不到一个稳定的状态.而由方法2模拟的结果来看，当气泡上升到一定高度后，气泡形态不在变化，因此气泡的速度最终也维持恒定.图4分别对利用两种方法模拟并排双气泡时，气泡体积偏差随时间的变化规律进行了对比.从图中可以看出，利用方法1得出的结果中，在初始阶段气泡的体积守恒性保持的很好，但随着时间的推移，体积偏差逐渐偏离1.0.而利用方法2获得的结果来看，气泡的体积偏差始终在1.0附近一个很小的范围内波动.由此可见，采用方法2对气泡的模拟更加准确，结果明显好于方法1.

图3 两种方法模拟出的气泡上升速度随时间变化

图4 气泡体积盈亏随时间的变化曲线

4 结论

本文通过建立气液两相流模型，提出改进的贴体坐标下Level Set方法，并结合考虑表面张力项采用CSF模型，对竖直流道内并排临近双气泡的上升运动过程进行了数值模拟.速度与压力的解耦采用基于同位网格下经典的SIMPLE算法，准确地追踪了气泡上升过程，获得了气泡运动过程形状，位置及速度等特性参数变化规律.采用基于体积校正法和局部重新初始化的Level Set方法，解决了Level Set方法中由于数值耗散引起的气泡体积的非守恒问题，获得了满意的结果.

［1］MONTES F J，GALAN M A，CERRO R L.Mass transfer from oscillating bubbles in bioreactors［J］.Chemical Engineering Science，1999，54:3127-3136.

［2］DARMANA D，HENKET R L B，DEEN N G，et al.Detailed modelling of hydrodynamics，mass transfer and chemical reactions in a bubble column using a discrete bubble model:Chemisorption of CO2into NaOH solution，numerical and experimental study［J］.Chemical Engineering Science，2007，62(9):2556-2575.

［3］HIRT C W，NICHOLS B D.Volume of fluid(VOF)method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，39(1):201-225.

［4］YOUNGS D L.Time-dependent multi-material flow with large fluid distortion numerical method for fluid dynamics［M］.New York:Academic Press，1982:273-285.

［5］OSHER S，SETIAN J A.Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations［J］.Journal of Computational Physics，1988，79:12-49.

［6］CHEN S，MERRIMAN B，OSHER S，et al.A simples Level Set method for solving stefan problems［J］.Journal of Computational Physics，1997，135:8-29.

［7］CHARRIER P，TESSIERAS B.On front-tracking methods applied to hyperbolic systems of nonlinear conservation laws［J］.SIAM Journal on Numerical Analysis，1986，23(3):461-472.

［8］WANG S L，SEKERKA R F.Algorithms for phase field computation of the dendritic operating state at large supercoolings［J］.Journal of Computational Physics，1996，127(1):110-117.

［9］PENG D，MERRIMAN B，OSHER S，et al.A PDE-based fast local Level Set method［J］.Journal of Computational Physics，1999，155(2):410-438.

［10］李向阳，王跃发，禹耕之，等.改进水平集法模拟不可压缩两相流动时质量守恒的体积校正法［J］.中国科学 B 辑:(化学)，2008，38(7):636-644.

［11］BRACKBILL J U，KOTHE D B，ZEMACH C.A continuum method for modeling surface tension［J］.Journal of Computational Physics，1992，100:335-354.

［12］PATANKAR S V.Numerical Heat Transfer and Fluid Flow［M］.New York:Hemispher，1980.

［13］陶文铨.数值传热学［M］.西安:西安交通大学出版社，2001.

［14］PERICCM，KESSLER R，SCHEUERER G.Comparison of finite-volume numerical methods with staggered and colocated grids［J］.Computers and Fluids，1988，16(4):389-403.

［15］RHIE C M，CHOW W L.A numerical study of t he turbulent flow past an isolated airfoil wit h t railing edge separation［J］.AIAA Journal，1983，21(11):1525-1532.

［16］刘儒勋，王志峰.数值模拟方法和运动界面追踪［M］.合肥:中国科学技术大学出版社，2001.