基于变密度流体层下气泡浮升运动的数值模拟分析

　　摘要：本文开发了一种耦合气体和两层不同密度液体的直接模拟数值模型,以研究变密度流体层中气泡浮升的运动特性.基于高精度空间离散的5阶加权本质无振荡格式和时间离散的3阶Runge-Kutta格式,提出一种多水平集函数耦合的数值方法,实现了气液相界面和各液体表面的捕捉以及交界面处密度和黏度的计算.模型用已有文献的计算结果进行了验证,并应用于水、油两层流体中气泡运动的研究,有效地捕捉了复杂条件下气泡浮升运动的精细结构.研究发现:相较于均匀密度流体层中的运动,气泡在变密度流体层中运动时,其下方的压差液柱更容易在底部发生断裂分离;气泡上方液体的厚度较大,则气泡处于水下变形阶段的时间就越长,并使气泡底部的形态更趋向于扁平.
　　关键词： 力学 变密度流体层 数值模拟 气泡 水平集法
　　气体在液体中的运动广泛存在于自然界和工程领域中,如发动机水下排气、水下爆炸引起的气泡运动、石油开采.其中的动力学过程在一个多世纪以来一直受到人们的关注[1,2].
　　数值模拟是研究气体在液体中运动的一种有效方法[3,4].由于液体与气体之间的密度比一般远大于重液体与轻液体之间的密度比,所以气液相界面的捕捉一直是数值模拟的重点与难点.目前,捕捉气液相界面最有效且最具有代表性的方法有两种,即流体体积(VOF)法和水平集法.VOF法由Hirt等[5]在1981年最先提出,其将供体-受体结合方程用于求解对流F函数,以确保数值解.但该方法若不进行几何重构,可能会产生较大的数值耗散,进而影响计算精度.
　　而在水平集法中,气液两相流中的相界面被设为零水平集.水平集函数是符号距离函数,相界面的一侧符号为正,另一侧符号为负,函数值的大小为空间位置与零水平集(相界面)的距离.据此,气液界面的形状和曲率可以较容易地通过求解水平集函数直接获得.水平集法由Osher等[6]在1988年最先提出,并将之用于捕捉锋面运动.该方法在形状拓扑上具有较大的优势,但可能伴有质量损失.Sussman等[7]首次将水平集法运用于计算不可压缩二相流.在此之后,水平集法在计算流体力学中的应用得到了很好的发展.Adalsteinsson等[8]提出水平集函数可以由一个扩展速度推进,而非流体速度.Jiang等[9]提出用高精度的加权本质无振荡(WENO)有限差分格式来求解流场可有效减少质量损失.现阶段的水平集法已经能够较好地捕捉变形剧烈的流体交界面[10,11,12].针对多相问题,Starinshak等[12]提出了当计算域中含多种材料时水平集函数的处理方法,但应用该方法时,一旦各水平集函数内部(正数部分)出现重叠,可能在重叠处出现严重的质量损失问题,甚至导致程序的崩溃.
　　当气体在层化的海洋环境等非均匀流体中运动时,由于周围液体的密度不恒定,其运动过程应当被视为一种在变密度流体层中的运动.此时,气泡周围的变密度流体层可能在此种复杂运动过程中扮演着十分重要的作用.
　　近年来,已有不少关于气泡动力学的研究论文,如梅登飞等[13]的黏性与非黏性颗粒在流化床中的气泡行为模拟.然而,虽然有一些气液耦合的研究,如李少白等[14]的非牛顿流体中线双气泡相互作用的数值模拟,但以上研究主要考虑两种流体间的相互作用,并没有考虑多种液体的存在,或液体内部间的物质交换对气泡上升运动的影响.此外,目前已有一些研究分析了流体黏性对气泡运动的影响[2],但关于气泡在变密度流体层中浮升运动的研究仍较少,尤其缺乏对气体运动穿过多种液体并进入大气环境这种复杂过程下的气泡浮升运动精细结构的捕捉.因此,为实现气泡在变密度流体层中浮升运动的有效模拟,并考虑到变密度流体层可以等效为若干双层流体的叠加,本文开发了一种耦合气体和两层不同密度液体的直接模拟模型,对气泡在变密度流体层中的运动进行计算分析.基于5阶加权本质无振荡和3阶Runge-Kutta高精度格式[4]提出一种多水平集函数耦合数值(MLS)法,用以捕捉气体表面以及各液体表面;使用投影法计算压力,并通过求解多物质流动系统的平滑化方程计算交界面处的密度和黏度.通过与已有文献结果的比较,验证了所建立的MLS法的有效性.采用该模型对水面薄油层的两层流体中的气泡运动进行研究,实现了对复杂条件下气泡浮升运动精细结构的有效捕捉.
　　1、数值模型
　　耦合气泡和不同密度流体的数值模型在计算中采用水和油作为典型变密度流体层,并将水作为底层流体,油作为上层流体,油的上方则为空气.模型设置如图1所示.其中:d1为气泡上方水层的厚度;d2为油层的厚度.
　　图1模型设置示意图
　　1.1 多水平集函数的构建及多物质流动系统的平滑化
　　首先,以气液相界面和各液体表面构建多个水平集函数ϕm(X)(下文简写为ϕm).其中:m为流体序号;X为点的空间坐标,X=(x,y,z).当m=1时,流体的种类为气体,则有
　　当m 2时,流体的种类为液体,则有
　　式中:Ωm为m号流体所在的区域;ψm为流体的表面微元集,即流体界面的水平零集;s为坐标点到流体界面的距离.
　　由于对控制方程进行离散化后,交界面处微元的密度和黏度不再能够简单地考虑为某一流体的密度和黏度,所以需要进行平滑化处理.因此,引入Heaviside函数H(ϕm),则有
　　式中:ε为界面厚度,即平滑层厚度的一半.引入H(ϕm)后,便可实现两相流的平滑化,如水气两相流的密度及黏度可以表示为
　　式中:ρw为水的密度;ρg为气体的密度;μw为水的黏度;μg为气体的黏度;ϕwg为以水气交界面为零水平集且水相内符号为正的水平集函数.
　　据此,先对气液相界面进行平滑化,而后按液体序号依此完成各液体表面的平滑化,最终实现多物质流动系统的平滑化.
　　式中:n为流体种类总数;ρm、μm分别为m号流体的密度和黏度.
　　1.2 控制方程及水平集演化
　　采用三维不可压缩Navier-Stokes方程组作为控制方程.该方程组包括连续性方程
　　和动量方程
　　式中: ˆ  ^ 为梯度算子;uˆu^为流体速度;pˆp^为压力;FˆF^为体积力;ρˆρ^为流体密度;μˆμ^为流体黏度;tˆt^为时间;控制方程中变量的上标＂^＂符号表示有量纲的物理量,不带上标则表示无量纲的物理量.由于主要考虑薄油层的密度和黏度对气泡运动的影响,且水、油间的表面张力对气泡的形态影响不大,水、气间的表面张力系数和油气间的表面张力系数相近[15],故统一取表面张力系数为7.20 10-2N/m.此时,结合水平集函数求解的密度及黏度形式,对式(8)和(9)进行无量纲化处理[16]可得
　　式中: 为无量纲梯度算子;u为无量纲流体速度;p为无量纲压力;Fsf为无量纲表面张力;ρ为无量纲流体密度,ρ=ρ(X);μ为无量纲流体黏度,μ=μ(X);t为无量纲时间;ez为重力方向无量纲单位矢量,ez=(0,0,1);Re为雷诺数;Fr为弗劳德数;We为韦伯数;uˆstu^st为特征速度;ρˆstρ^st为特征密度;μˆstμ^st为特征黏度;lˆstl^st为特征长度;gˆstg^st为重力加速度;σˆstσ^st为特征表面张力系数.Fsf的计算公式为
　　式中:ϕ为以气液相界面为零水平集的水平集函数,等效于ϕm=1;κ(ϕ)为曲率;δ(ϕ)为Diracδ函数.使用投影法求解压力项和下一时刻的速度.根据水平集演化方程[4]得到
　　以及加入质量修正的重距离化方程[4]
　　求解下一时刻的水平集函数,并依此计算下一时刻流场的密度和黏度分布.式中:τ为人工时间;ϕm0为τ=0时即求解式(18)时的ϕm;S(ϕm0)为平滑符号函数;λ为质量补偿系数;Ω为计算域.
　　1.3 数值计算步骤
　　MLS法的主要计算流程如下:
　　(1)构建多个水平集函数;
　　(2)耦合所有水平集函数,求解整个流场的密度及黏度分布;
　　(3)根据整个流场的密度及黏度分布和现有速度场,计算压力场和下一时刻的速度场;
　　(4)根据新的速度场,计算新的各水平集函数;
　　(5)重复步骤(2) (4)直到模拟时间结束.
　　计算中的各流体间互不混溶,并可近似为不可压缩流体[16].在此条件下,气泡内的初始压力为1个标准大气压(1.013 105Pa),对应的气体初始密度为1kg/m3,水、气间的密度比为1000 1.
　　计算采用结构网格,通过有限差分法对控制方程中的各项进行离散.同时,对于水平集法,采用具有5阶精度的WENO格式对空间项进行离散,并采用具有3阶精度的Runge-Kutta格式对时间项进行离散.计算中需要满足以下Courant-Friedrich-Levy(CFL)条件,即可保证模型计算时的收敛性[4]:
　　式中:u、v、w为空间直角坐标系的无量纲速度分量.
　　2模拟验证
　　通过与文献的模拟气泡在自由表面爆裂和水中气泡上浮结果的对比,以验证MLS法的有效性.
　　2.1 气泡在自由表面爆裂算例
　　以Gu等[4]的气泡在自由表面爆裂算例作为验证算例1(C1),该算例为在[-3,3] [-3,3] [-6,6]的计算域中,设置一个半径为1,球心坐标(x,y,z)=(0,0,-3.2)的气泡,水深区间为[-6,-2.0],水的上方为气体.C1示意图如图2所示.
　　图2C1示意图
　　针对该算例,选用计算网格数为110 110 220,采用无滑移边界条件,取Δt=0.005Δx,Re=474,Fr=0.64,We=1.0,ρg=0.001ρw,μg=0.01μw.
　　图3水气交界面的演变过程
　　本文模型(M1)的计算结果与C1的模拟结果对比如图3所示.由图3可知,M1与C1的结果一致表明整个流场的演化过程可定性地分为3个阶段.第1阶段为水下变形阶段,由于气泡各处的压力不同,气泡会在上升的同时发生形变;因底部受到的压力较大,表面张力和黏性力难以维持其原有形态,气泡底部将发生凹陷.第2阶段为融合阶段,气泡在浮力的作用下,穿透液体表面,在液体中部形成开口;此后随气泡继续上升,开口将进一步扩大.第3阶段为惯性阶段,此时气泡完全与上方气体融合,伴随上升的液体形成液柱结构;后续在重力和黏性力的主导作用下,由于液柱内各处速度不同,液柱进一步发生分离现象.另一方面,对比结果表明,M1精确地捕捉到了液柱二次分离过程中产生的细小流体结构,具有更好的连续性和质量守恒性.
　　2.2 水中气泡上浮算例
　　以Premlata等[2]的水中气泡上浮算例作为验证算例2(C2),该算例的参数为Re=100,Fr=1.0,We=200,ρg=0.001ρw,μg=0.01μw.针对该算例选用的计算网格数为100 100 200,采用无滑移边界条件,并取Δt=0.005Δx.本文模型(M2)的计算结果与C2的模拟结果的对比如图4所示.由图4可知,直到t=0.4时,气泡都几乎是球形的;当t>0.4时,气泡底部持续发生凹陷;当t=1.6时,气泡发生显著的变形.此后,气泡继续向上移动,形成环形的＂甜甜圈状结构＂.M2与和C2所获得的气泡形态演变过程的模拟结果具有良好的一致性.
　　图4气泡形状演变过程
　　2.3 质量守恒性分析
　　为进一步精确地检验MLS法的有效性,计算每一时刻的相对质量变化率ηmc,则有
　　式中:mc为计算质量.M1和M2中的ηmc随时间的变化如图5所示.气泡在自由表面爆裂时,水气相界面的变化明显比在水中上浮时水气相界面的变化更剧烈.这一点与前者的质量损失较大相符合.同时,两算例中的最大相对质量变化率均小于0.6%[4],说明了MLS法已具有良好的质量守恒性,已能够以较小的质量损失较为精确地捕捉气泡的运动.
　　图5相对质量变化率
　　3、变密度层中气泡浮升运动的数值模拟
　　气泡在自由表面爆裂时,液气交界面的变化十分剧烈,流场中的水气混合过程相当复杂.为对该过程进行研究,设置了4组算例(MF1 MF4),各算例的主要参数如表1所示.4组算例中,气泡的半径为1,球心为(x,y,z)=(0,0,-3),水深区间为[-6,-2.0].4组算例的网格数,无量纲数以及离散方法等均与C1相同.在MF2和MF3中,取油的密度ρo=0.8ρw及黏度μo=2.0μw.
　　图6MF1的模拟结果
　　图7MF2的模拟结果
　　图8MF3的模拟结果
　　图9MF4的模拟结果
　　MF1和MF2的模拟结果分别如图6和7所示,MF3和MF4的模拟结果分别如图8和9所示.由图6 9可知,无论是均匀流体层中的气泡运动,还是变密度流体层中的气泡运动,整个流场的运动过程均可定性地分为3个阶段:水下变形阶段、融合阶段以及惯性阶段.
　　对比图6(b)与图8(b)、图7(b)与图9(b)可知,在水下变形阶段时,增大气泡上方液体层的厚度会增加气泡处于该阶段的时间.例如,当气泡上方的液体层厚度d(d=d1+d2)从0.2增大到0.4时,该阶段的持续时间约从0.3增大到0.6.同时,厚度的增加会改变气泡与上方气体融合时的形状和运动状态,进而影响到后续两个阶段的变化.另外,气泡上方液体层的种类对水下变形阶段的气泡形态影响较小,这主要是由于油层的厚度较薄,变密度层中的气泡和均匀流体中的气泡在液体下所受到的压力基本相似.
　　在融合阶段,当气泡在上方较厚的液体层作用下,产生了较大的形变而变得足够扁平时,气泡上下的压力差将会减少,而压力作用的面积则会增大,使得压差液柱的横截面积增加.另一方面,虽然厚度相同的均匀水体与变密度流体层中的气泡运动具有相似的发展过程,但是由于变密度层的非均匀性,各流体间的相互作用更为复杂.具体而言在该阶段中,油和水体的上部分会因受气泡给予的推力而向两侧边界运动,而其下方的水体则会因受压力而向上运动形成压差液柱.之后,在四周边界的约束下,部分液体回流.由于油与水互不混溶,油体将主要在水体的表面流动.
　　在惯性阶段,压差液柱的发展是流场变化的主要特点.此时,流场主要受惯性、黏性力以及重力的控制.在该阶段早期,虽然受到向下的有效重力,但在惯性的作用下,压差液柱将继续上升.然而,相较于均匀流体层中的运动,由于油层和水体是互不混溶的,当回流的流量较大时,油层回流会冲击压差液柱的底部,对压差液柱产生剪切和挤压的作用,使得黏性力难以保持压差液柱和下方水体的连续,最终产生分离现象.
　　图10zh演化图
　　4个算例(MF1 MF4)中液体能达到的最大高度随着时间变化的结果如图10所示,其中zh为液体最高点在z轴方向上的无量纲坐标值.由图10可知,当d=0.2时,MF1和MF2中水下变形阶段的时间均约为0.3,同时最大zh均大于3.0.在随后的融合阶段中,由于气泡两侧的液体向气泡底部流动,液面的最大高度会有所下降.在惯性阶段中,由于气泡已经与上方气体融合,液体内部间的相互作用成为影响液气界面的主要因素.因此,MF1与MF2中液气分布的区别将进一步增大,两者zh的最大差值可达0.5.另一方面,当d=0.4时,由于气泡处于水下变形阶段的时间较长,气泡形态趋于扁平,气泡上下两侧的压力差较小,故最大zh将比d=0.2时减小约2.8.
　　4、结论
　　本文针对气泡在变密度流体层中的运动开发了一种耦合气体和两层不同密度液体的直接模拟模型,并基于水平集法提出了一种求解多水平集函数耦合的数值方法.通过与已有数值模拟结果的对比验证了MLS法的有效性.采用该模型对水、油两层流体中气泡的运动进行了研究,有效地捕捉了复杂条件下气泡浮升运动的精细结构.研究结果表明,当液体和液体中的气泡可视为不可压缩流体,且不考虑各液体间的扩散作用时,主要结论如下.
　　(1)若边界对流体的运动具有一定的约束,相较于均匀密度流体层中的运动,由于边界的约束,油层将会回流,从而对水体产生剪切挤压的作用,使得气泡在变密度流体层中运动时,其下方的压差液柱将更容易在底部发生断裂分离.
　　(2)增大气泡上方初始液体厚度会使气泡处于水下变形阶段的时间更长,同时气泡底部的形态也将更为扁平.