改进水平集法模拟不可压缩两相流动时质量守恒的体积校正法

中国科学B辑:化学2008年第38卷第7期:636~644改进水平集法模拟不可压缩两相流动时质量守恒的体积校正法李向阳①,王跃发①,禹耕之①,杨超①②,毛在砂①*①中国科学院过程工程研究所,绿色过程与工程重点实验室,北京100190②江苏省海洋资源开发研究院,连云港222005*联系人,E-mail:[e-mail]收稿日期:2007-10-19;接受日期:2008-03-26国家自然科学基金(批准号:20490206,50404009,20576133,20676134),中国石油天然气集团总公司和国家基础研究项目(批准号:2004CB217604,2007CB613507)资助摘要发展了一种体积校正法来改进水平集法模拟涉及相界面变形的不可压缩两相流动时的质量守恒.在传统的重新初始化程序之后,该体积校正法可以根据不连续相的质量亏盈,通过整体移动相界面的位置,使质量亏盈小于预先指定的误差范围.用3个算例验证了包含该体积校正法的水平集算法:二维单气泡或单液滴在液体中的运动,二维水滴通过空气介质落入水池,两个浮力推动的三维变形气泡相互影响下的运动.每个算例中质量守恒都符合要求的精度,数值模拟结果与实验数据吻合.关键词水平集算法质量守恒重新初始化体积校正气泡液滴1前言在自然界和各种工程领域中经常出现多相流动问题.由于计算机性能和计算流体力学方法的快速发展,近年来数值模拟方法越来越多地应用于多相流的研究中.在水平集算法中,两流体分界面用水平集函数(定义为距界面的代数距离)的零水平集隐含地给出,该水平集函数可以很容易地用一个简单的对流方程随着时间推进来更新[1,2].水平集方法现在已经广泛地应用于流体力学、燃烧、计算机显示和材料科学的研究中[3,4].然而,水平集函数的推进是一个非守恒的算法.随着计算时间增长,总是会出现气泡和液滴的质量亏盈的问题,使相界面逐渐偏离正确的位置,同时使流场失真.到目前为止,文献中报道了3类改善质量守恒的方法[5].第一类就是目前广泛采用的改进重新初始化方法[6~10].Sussman等[6~8]发展了一种类似于VOF(volumeoffluidmethod)方法的重新初始化程序.他们认为在重新初始化过程中,由于界面是不动的,那么每个网格中体积也是不变的,提出采用紧致的九点积分格式计算每一个网格单元的体积,根据网格单元体积与初始体积的差别来重新初始化水平集函数,这大大增加了计算负担和程序的复杂程度.与该方法不同,Chang等[9]注意到导致气泡液滴体积亏盈的数值扩散与局部曲率有密切关系,提出在每一个时间步还同时求解一个水平集函数的Hamilton-Jacobi方程直至稳态,得到了二维模拟的精确解.Wang和Zhou[10]认为重新初始化过程中体积亏盈有两个原因:一是数据之间的相互影响,二是界面两侧中国科学B辑:化学2008年第38卷第7期637信息的不对称性.他们通过锁定界面位置,引入了一种新的重新初始化水平集函数的算法.该方法针对的是界面而不是直接处理体积,必然会存在体积的亏盈或波动.除此之外,以上提到的改进的重新初始化方法都只是在二维的计算中得到了验证,在全三维问题模拟中达到同样的效果将是个更大的挑战[11].Oka和Ishii[11]在三维问题的模拟中曾设法很好地保持了质量守恒,但是在每一个时间步内同时使用三种重新初始化方法非常耗用计算机时间.第二类是几种“杂化”方法,比如把水平集方法和volume-tracking[12~14]或者front-tracking[15,16]方法结合起来.例如Enright等[15]的方法是把水平集函数和一些Lagrangian标记粒子结合起来,该方法能够保持界面几何形状的光滑描述,同时大大改进了under-resolved流动的质量守恒.但这些方法往往比较复杂,会大大增加计算负荷,尤其是当应用于三维计算时尤其突出.第三类,通过减小空间离散误差的方法来改进水平集方法中的质量守恒问题.能够更准确的求解界面结构的AMR(adaptivemeshrefinement)方法也曾被用来改进质量守恒问题[7,12,16~18].这类方法中一个典型例子包括两个协同的步骤:(1)用线性WENO(weightedessentiallynonoscillatory)格式来减小空间离散误差;(2)用SAMR(structuredadaptivemeshrefinement)方法提高局部解的精度.两个协同步骤大大增加了程序的复杂程度和计算的负担,尤其是在非并行的单机计算时.本文中提出了一种改善水平集方法应用中质量(体积)守恒的体积校正法,该方法不会增加计算代码的复杂程度,也不会显著地增大计算负担.2控制方程本研究中作了以下的合理假设:(1)两种流体是黏性的、牛顿型的、不可压缩的流体;(2)流动是等温的;(3)表面张力是一个定值[17].无因次形式的运动方程和连续性方程为:()()()T11,ptReFrWeμρρκφδφε∂⎛⎞⎛⎞+⋅∇=−∇+∇⋅∇+∇⎜⎟⎜⎟∂⎝⎠⎝⎠++UUUUUgn(1)0,∇⋅=U(2)其中,al/ρρρ=和al/μμμ=分别是分散相密度和黏度与连续相参数的比值,界面的单位法向量n定义为.φφ∇=−∇n(3)磨光的()δφε定义为:()()()11cos/if,20otherwise.φφδφε⎧+πεε⎪=ε⎨⎪⎩(4)()κφ是自由界面的曲率:().φκφφ⎛⎞∇=∇⋅=−∇⋅⎜⎟⎜⎟∇⎝⎠n(5)在水平集算法中,密度场和黏度场随着流场推进.为了避免在界面附近,尤其是大密度梯度条件下出现数值不稳定性,对密度和黏度函数作如下平滑处理:()alal()/1/(),Hρφρρρρφε=+−(6)()alal()/1/().Hμφμμμμφε=+−(7)其中的磨光的Heaviside函数定义为:()()0if,11sin//if2ifHφφφφεφφε−ε⎧⎪⎪⎛⎞=++ππε,⎨⎜⎟ε⎝⎠⎪⎪1ε.⎩≤(8)方程(1)中的无因次数Re,Fr和We分别定义为:ll,LVReρμ≡2,VFrLg≡2l,VLWeρσ≡(9)其中L和V分别为特征长度和特征速度.水平集方程为:0.tφφ∂+⋅∇=∂U(10)Sussman等[19]提出的重新初始化方法是积分下列方程至稳态:0sgn()(1||),φφφτ∂=−∇∂(11)方程中τ是重新初始化过程中迭代的虚拟时间,0φ是t时刻的水平集函数,0sgn()φ是为了避免数值不稳定而磨光的符号函数.为了获得更高的空间精度,采用五阶的weightedessentiallynonoscillatory(WENO)李向阳等:改进水平集法模拟不可压缩两相流动时质量守恒的体积校正法638格式来离散重新初始化方程[20].Chang等[9,21]提出的改进质量守恒的重新初始化方程:0(())(())0,AAPφτκφφτ∂+−−+∇=∂(12)其中0A表示初始时刻的总质量,()Aτ表示与()φτ对应的总质量.参数P是一个正的常数,一般取为1.0.Yang和Mao[22]在计算中发现,使用方程(12)来重新初始化,尽管两相流体的总质量守恒很好,然而流体颗粒(气泡或液滴)的质量却逐渐丢失.因此,他们将()Aτ的定义修改为:()(),ijijjArxyφτρφεε=ΔΔ∑≤(13)式中ijφε≤表示流体颗粒内部和虚拟界面厚度2ε内的节点,也就是说方程(13)中的()Aτ取的是气泡/液滴的质量而不是方程(12)中所取的体系总质量.3改进质量守恒的体积校正法由于已假设流体是不可压缩的,所以质量守恒和体积守恒是等价的.本文的重新初始化过程中,首先求解方程(11)使水平集函数成为代数距离函数,然后根据下边公式所示不连续相(气泡/液滴)体积(0()VHrxyφφ=ΔΔ∑≤)的亏盈来调整相界面的位置:()00.VtVVV−Δ=(14)假设气泡/液滴是球形的,半径的增量为:()110011,CrrrVr⎛⎞Δ=−=+Δ−⎜⎟⎝⎠(15)其中1C取模拟体系的维数.由于分散相体积的增加对应于从界面到分散相中心距离的增加,所以可取水平集函数的校正量正比于Δr.于是水平集函数的校正量可以表达为:()11211,CCVδφ⎛⎞=+Δ−⎜⎟⎝⎠(16)这里的2C是一个经验参数,较大的2C可能导致重新初始化方程不收敛,较小的2C降低计算效率.2C一般在0.01~0.1中按经验选取.这样,水平集函数的体积校正方程可以写为:0.φφδφ=+(17)本文采用控制体积法和Patankar[23]所描述的幂指数方案来求解方程(1)和(2).耦合了本文提出的体积校正的重新初始化程序的水平集方法的主要计算步骤如下:1.初始化流场(U和p)、各相的物性参数(ρ,μ和σ),初始化(,0)φX为相对于相界面的符号距离函数;2.使用压力耦合的半隐式SIMPLEC算法在每一个时间步内求解方程(1)和(2)得到速度场和压力场;3.在每一个时间步,利用已得到的速度场求解方程(10)来更新水平集函数φ;4.通过求解方程(11)至稳态,使φ更新为代数距离函数;5.判断分散相(气泡/液滴)体积的亏盈是否超出了(预先指定的)标准.如果没有超出标准,接受步骤4得到的φ为新的水平集函数,并根据它来更新物性;然后回到步骤2开始下一时间步的计算;否则,求解方程(17),然后回到步骤4.4结果和讨论在各种工业过程中,常常会遇见具有鲜明相界面的多相流动,实践证明很难通过数值方法模拟这类流动.因此本文从该类流动中选择3个典型的问题作为算例,来演示本文提出的体积校正的重新初始化程序:二维单气泡或单液滴在液体中的运动、二维水滴通过空气介质落入水池和两个浮力推动的三维变形气泡的相互作用运动.4.1二维轴对称单气泡或单液滴在液体中的运动在过去的几十年里,许多研究者通过实验和数值模拟的手段对气泡和液滴的运动进行了研究[11,22,24~29].本算例的重点放在确定气泡的稳态形状和上升速度上.和Yang和Mao[22]的模拟一样,计算的区域为Ω={(,)|010,08}xyxRyR≤≤≤≤,采用100×80非均匀网格进行计算,该设置可以得到足够的空间精度.首先模拟了空气泡分别在高黏度和低黏度的糖浆水溶液中自由上升运动.数值计算得出平衡态气泡的上升速度分别为0.181和0.317m/s,与实验结果0.190和0.306m/s[24]非常接近.图1和2所示的气泡稳定态形状也与实验可视化结果吻合较好[24].中国科学B辑:化学2008年第38卷第7期639为了验证气泡是否适用传统的零密度、无黏度的传统观点,模拟了气泡密度和黏度的量级对计算的影响.成功地模拟了零密度、无黏度的气泡在糖浆中的上升运动,如图3所示,与Yang和Mao[22]的模拟结果比较一致,但与图1中的气泡形状有明显差别.此外,目前的算法还能够用来捕捉运动中流体颗粒的拓扑变形.模拟得到的存在较大密度比和有拓扑变形的运动气泡的稳定态形状如图4所示,其中无因次时间的定义是//(2),tgRθ=R是气泡半径.以上模拟结果表明,在层流条件下,使用欧拉网格,耦合了本文提出的体积校正的重新初始化程序的水平集方法对这类相界面存在较大的密度和黏度梯度的稳态流体颗粒的运动是很有效的.4.2水滴通过空气介质落入水池水滴通过空气介质落入水池中的运动是有趣的自由界面问题[1].在这一过程中,会有一系列复杂的界面拓扑变形发生,对数值模拟提出了挑战.图1高黏度糖溶液中9.3cm3气泡稳态运动时的形状与实验的比较(a)数值模拟的气泡;(b)实验拍的照片(Bhaga和Weber,图3(b))[24](R=0.013m,ρl=1390kg·m−3,ρg=1.226kg·m−3,μl=2.786Pa·s,μg=1.78×10−5Pa·s,σ=0.0794N·m−1)图2低黏度糖溶液中9.3cm3气泡稳态运动时的形状与实验的比较(a)数值模拟的气泡;(b)实验拍的照片