三相流体的轴对称格子 Boltzmann 模型及其在 Rayleigh-Plateau 不稳定性的应用*

刘程 梁宏

(杭州电子科技大学物理系,杭州 310018)

基于多组分相场理论提出了一类模拟三相流体流动的轴对称格子Boltzmann 模型.该模型利用两个粒子分布函数来捕捉三种不同流体之间的相界面,另一个粒子分布函数来求解流体动力学方程以获得流场信息.为了刻画坐标变换引起的轴对称效应,巧妙地设计了演化方程中平衡态分布函数和外力项分布函数,从理论上保证本文模型可以正确恢复三相流体系统的宏观控制方程,并且轴对称效应产生的源项中不包含任何复杂的梯度项,从而比现有的轴对称格子Boltzmann 模型更加简单高效.首先通过模拟一系列轴对称多相流的基准算例,包括静态的双液滴、液体透镜的扩展和二元流体Rayleigh-Plateau 不稳定性,来验证本文模型的有效性与正确性.接下来,利用该模型研究了三相流体的Rayleigh-Plateau 不稳定性的增长过程,定量分析了波数和液柱半径比对复合液体线程破裂过程中界面动力学行为、界面破裂时间以及生成子液滴尺寸的影响.可以发现复合的液体线程在波数较大时破裂生成一个复合主液滴和卫星液滴,而在波数较小时可以生成更多数量的卫星液滴,这导致复合主液滴和卫星液滴的尺寸随着波数的增加呈现先增大而后减少的趋势.另外,我们发现内部流体优先于中间流体发生界面破裂,并且流体界面的破裂时间均随着波数的减小而增加.最后发现,增大液柱半径比可以促进内部液体线程的破裂,而延缓中间液体线程发生破裂,并且复合主液滴的尺寸随着液柱半径比的增加而增大,而复合卫星液滴的尺寸对液柱半径比的变化不显著.

1 引言

三相流体系统广泛存在于自然界和工程应用中,如水污染物的输运、强化采油中油水气的传输过程、微流控装置中复合液滴的生成技术等.由于涉及多种流体界面的迁徙、变形、破裂及合并等现象,三相流体的界面动力学行为十分复杂.随着计算机技术的快速发展,数值模拟已成为研究三相流体流动问题的有效手段.它可以方便地提供实验研究中不易测量的物理量,如在相界面变化过程中某一时刻的界面形状、界面处流体的速度、压力及密度分布等.根据不同的数值方法,一些学者对三相流体系统进行了研究,并在笛卡尔坐标系下建立了多种数值模型,如Smith 等[1]开发了一种水平投影方法,用于在水平集公式中处理多个连接的运动,并对三相流系统进行了研究;Bonhomme 等[2]发展了三相流的流体体积函数法,用于研究气泡穿过两相流体界面的动力学行为;Kim[3]基于相场理论,提出了包含表面张力的三组分不混溶流体流动模型.

格子Boltzmann (lattice Boltzmann,LB) 方法[4]是近三十年来发展起来的流体系统模拟方法,通过描述粒子分布函数的演化再现流体的宏观行为.由于其微观本质和介观特点,LB 方法相比传统数值方法具有算法简单、天然并行、易于处理流体间的微观相互作用等优点,被广泛用于模拟多相流体[5,6]以及三相流体输运问题.Spencer 等[7]基于两相流LB 方法中颜色模型,重新设计了三相流体间的颜色梯度,从而建立了三组分流体的颜色模型.Leclaire 等[8]在颜色梯度的LB 模型中引入了三个碰撞算子,并使用高阶离散算子计算颜色梯度,从而提出了数值稳定性更好的三相流LB 模型.Yu 等[9]基于Leclaire 等[8]开发的微扰算子,推导出了描述不同流体之间相互作用的界面力公式,又应用了Spencer 等[7]提出的重着色算法来维持界面并确保流体的不混溶性,提出了用于模拟具有表面张力的不混溶三相流体的颜色梯度LB 模型.Semprebon 等[10]基于自由能泛函理论,推导出流体间和流固间的相互作用,从而建立了表面张力可以独立调节的三相自由能LB 模型.Wöhrwag 等[11]将非理性流体状态方程耦合至三相自由能泛函中,并采用了熵的碰撞算子以提高模型的稳定性,从而发展了满足热力学一致性、稳定性好的三相流LB 模型.Liang 等[12,13]基于多组分的相场理论,建立了模拟非混相三相流体输运的LB 模型.该模型巧妙地设计了平衡态分布函数和外力项分布函数,从而理论上保证可以正确地恢复多组分的Cahn-Hilliard 方程和不可压Navier-Stokes 方程.

上述模型的提出促进了LB 方法在三相流领域的发展.然而,需要指出的是,这些模型均基于笛卡尔直角坐标系建立的.而在三相流体系统中,存在着相界面显示为轴对称特性的问题,如两种不同液滴的对头碰撞、气泡穿过液-液两相界面的上升动力学过程、三相Rayleigh-Plateau 不稳定性问题等.当采用LB 方法来模拟此类轴对称三相流动问题,最直接的方式是利用三维的三相流LB 模型和合适的曲边界处理格式.然而,这种方式没有利用轴对称流动的特点: 具有轴对称特性的流体系统可以将三维流动转化为子午面上的二维问题,大大减少了计算量和数据储存.为了充分利用轴对称流动的特点,一些学者考虑圆柱坐标系下建立了有效的多相流LB 模型.然而,绝大多数的轴对称多相模型是基于两相流情形[14],针对圆柱坐标系下的三相流LB 模型的研究鲜有报道.Haghani 等[15]从Cahn-Hilliard 相场理论出发,建立了轴对称坐标下的三相流LB 模型,但该模型中由轴对称效应引入的源项包含许多复杂的梯度项,这导致模型的算法复杂且梯度项的离散也会带来额外的数值误差.另一方面,针对三相流体Rayleigh-Plateau 不稳定性的研究工作也较少.Yang 等[16]利用轴对称的相场算法,对三相流体的Rayleigh-Plateau 不稳定性问题开展了初步研究,定性分析了黏性比、界面张力和初始扰动对复合线程的界面动力学的影响.

鉴于上述分析,本文将基于多组分的相场理论,提出一个简单高效的轴对称三相流LB 模型.该模型利用了合适的平衡态分布函数和外力项分布函数,从理论上保证可以正确地描述轴对称三相流体流动,并且由轴对称效应引入的源项不包含任何的梯度项,从而比已有轴对称三相流LB 模型更加简单.通过模拟静态的双液滴、液体透镜的扩展和二元流体Rayleigh-Plateau 不稳定性算例来验证本文LB 模型的准确性,并以此模拟了三相流体Rayleigh-Plateau 不稳定性的增长过程,详细分析了波数和内外液柱半径比对复合液体线程的界面动力学行为、界面破裂时间及生成子液滴尺寸的影响机制.

2 数值模型

2.1 宏观控制方程

本文模型建立在一个含有三种不互溶流体组成的不可压缩系统,并引入ϕ1,ϕ2,ϕ3三个序参数来描述三种流体在混合物所占的体积分数,且满足如下的约束条件:

根据多组分流体的相场理论[17,18],三相流体系统的总的自由能泛函可以表述为

其中,Ω为三相系统所处的流域,F(ϕ1,ϕ2,ϕ3) 为体相区的自由能,D为界面厚度的特征尺度,λi为与表面张力有关的物理参数.在相场理论[19]中,化学势定义为自由能泛函对序参数的变分.然而,针对三相流体系统,为了满足序参数的守恒条件,需要在自由能的变分项中引入一个拉格朗日乘数β,从而化学势µi可以表述为

序参数ϕi的时间演化由流体速度的对流和化学势的扩散来描述,其宏观控制方程可以写成多组分Cahn-Hilliard 方程的形式:

其中,u为流体速度,Mi为迁移率.假设迁移率满足约束条件M1λ1=M2λ2=M3λ3=M0,并对方程(4)中所有i进行求和,可以得到关于变量S的方程:

其中S=ϕ1+ϕ2+ϕ3.根据方程(1)给出的守恒条件,S=1 显然是方程(5)的解,可以推导出

其中,λT定义为

将β的表达式(6)代入方程(3),可以得到化学势µi的表达式为

为了使相场模型满足减少一致性条件以及具有很好的适定性属性,Boyer 等[17,18]理论构造了如下的体相区自由能:

其中,λi是与表面张力有关的参数,

其中,σ12,σ13,σ23分别表示三相体系中两种流体之间的表面张力.将(9)式代入(8)式,最终可以获得化学势的表达式:

其中,p是压力,ρ是流体密度,ν是运动学黏性系数,G是外力,Fs是表面张力.另外,为了减少相界面处的虚假速度,本文选用势能形式的表面张力

为了获得轴对称三相流体动力学的宏观控制方程,进行了从笛卡尔坐标系到柱坐标系的变换

其中,x=rcosθ,y=rsinθ,r,z,θ分别表示径向坐标,轴向坐标和方位角.经过一系列的代数计算,可以得到柱坐标系下的宏观控制方程:

2.2 轴对称多相流的格子Boltzmann 模型

本文需要利用两个粒子分布函数的LB 演化方程来捕捉三相流体间的相界面,另一个粒子分布函数的LB 方程来求解流场动力学特性.采用最简单的单松弛碰撞模型,追踪三相流体相界面的LB 演化方程可以表述为[12]

其中,(x,t) 是粒子在t时刻在x位置的序参数分布函数,为对应的平衡态分布函数,τi是无量纲松弛时间,δt是时间增量,(x,t) 是外力项分布函数.将(16)式的最右端项吸收至对流项,从而巧妙设计了序参数的平衡态分布函数:

其中,ωk和ek分别是权重系数和离散速度,cs是声速,η是引入的一个自由参数,用来调整迁移率系数.采用最具有代表性的D2Q9 离散速度模型[20],其对应权系数的取值为ω0=4/9,ω1-4=1/9,,离散速度ek定义为

式中c=δx/δt表示格子速度,δx代表格子步长.为了能够准确地恢复轴对称的多组分Cahn-Hilliard方程,还构造了全新的外力项分布函数:

在本模型中,通过计算粒子分布函数的零阶矩,来获得序参数ϕi的宏观量

另外,序参数ϕ3的值可以通过守恒条件(1)得到,即ϕ3=1-ϕ1-ϕ2,并且流体的密度ρ可以看成序参数ϕi的线性插值函数:

其中,ρ1,ρ2,ρ3分别表示三相流体系统所对应的相密度.利用Chapman-Enskog 多尺度分析技术[21],可以证明本模型可以准确地恢复到三相流体界面追踪的轴对称Cahn-Hilliard 相场方程以及迁移率Mi和松弛因子之间的关系:

由(24)式可知,当迁移率非常小的情况下,无量纲松弛时间τi将会接近 0.5,这可能会导致LB 方法的数值不稳定,因此引入自由参数η是非常有必要的,且不会增加任何额外的计算成本.

对于三相不可压缩流体的建模,需要将多组分Cahn-Hilliard 相场方程与Navier-Stoke 流动方程相耦合.基于单松弛碰撞模型,描述流体动力学的的LB 演化方程可以写成[12]

其中,gk(x,t) 为描述流场的粒子分布函数,(x,t)为平衡态分布函数,τg为与流体黏度有关的松弛因子,Gk(x,t) 为外力项分布函数.考虑到轴对称效应带来的影响,为了正确的恢复连续性方程,本文采用了如下修正的平衡态分布函数[14]:

另外,需要考虑将外力项引入介观LB 方法中所产生的格子离散效应.Liang 等[22]近期考虑了格子离散效应,设计了一类准确并且简化的外力项分布函数.在轴对称坐标系下,流场的外力项分布函数可以表示为

通过计算粒子分布函数的一阶矩和零阶矩,可以获得流体流动的速度

和流体动力学压力p为

利用多尺度理论分析手段[22],可以证明本模型可以恢复正确的轴对称不可压Navier-Stokes 方程,并且运动性黏性与松弛因子之间的关系为

在数值模拟中,需要采用合适的差分格式对模型中时间导数项和空间导数项进行离散运算,本模型采用显式的欧拉差分格式计算时间导数项,

以及各向同性的二阶中心差分格式计算梯度算子和拉普拉斯算子,

其中χ表示任意变量.另外,为了避免轴对称LB 模型由于涉及r-1的项而在r=0 处产生奇异现象,将第一层网格放在r=0.5δx处,并对轴向边界处应用对称边界条件[14].

3 数值结果和讨论

本节首先通过模拟几个典型的轴对称多相流问题来验证本文提出的介观LB 模型,包括静态的双液滴、液体透镜和二元流体Rayleigh-Plateau 不稳定性问题.接着,利用该模型来研究三相流体的Rayleigh-Plateau 不稳定性问题,详细分析波数和内外液柱的半径比对液线破裂过程中界面运动过程以及产生的复合液滴尺寸的影响.

3.1 静态的双液滴

静态的双液滴作为三相流的基准算例,首先被用来验证本文所建立的轴对称三相流 LB 模型.在静态的双液滴测试中,将计算网格设置为Nz×Nr=400×100,其中尺寸相同的两个液滴静止放置于计算区域中,液滴的半径为R=50,两液滴的中心坐标分别是 (z1,r0)=(100,0) 和 (z2,r0)=(300,0).初始的序参数分布设定为

而序参数ϕ3可以由方程(1)得到,其他物理参数设置为: 三相系统密度ρ1=10,ρ2=5,ρ3=1,界面厚度D=4,运动学黏性ν=0.1,界面张力σ12=σ13=σ23=0.1,松弛时间τ1=τ2=0.8,迁移率M0=0.01.粒子分布函数fi和gi的边界条件设置如下: 左右边界运用周期边界条件,上壁面采用无滑移反弹边界条件,下边界则采用对称性边界条件.图1 给出了稳态时三个序参数在计算区域中的分布.可以看出,序参数的稳态分布符合初始条件给定的双曲分布,避免了三相流的数值方法可能会产生非物理的虚假相.进一步统计了序参数沿着界面方向的稳态分布以及相应的解析解,并将结果列于图2.从图中可以观察到数值预测的序参数分布与解析解能很好地符合,这表明本文提出轴对称三相流的LB 模型能够准确捕获界面的分布.

图1 静态双液滴测试中三个序参数的稳态分布 (a) ϕ1 的分布;(b) ϕ2 的分布;(c) ϕ3 的分布Fig.1.Steady-state distributions of three order parameters in the static double-droplet test: (a) Distribution of ϕ1 ;(b) distribution of ϕ2 ;(c) distribution of ϕ3.

图2 静态双液滴测试中序参数沿水平线的稳态分布Fig.2.Steady distributions of the order parameters across the horizontal line in static double-droplet test.

3.2 液体透镜的扩展

液体透镜的扩展问题被广泛用于测试所提出的三相流LB 模型,然而绝大多数工作限制于二维笛卡尔坐标系[9,10,12],针对三维或柱坐标空间的模拟鲜有报道.本节将模拟三维液体透镜的扩展来验证所提出的三相流轴对称LB 模型.如图3 所示,圆形的液体透镜初始时刻放置于另外两种不混溶的流体之间,在表面张力作用下会发生变形,并最终达到稳定的形状.根据Neumann定律[23],液体透镜在平衡时的形状由三个表面张力系数决定:

图3 液体透镜相关几何参数示意图Fig.3.Schematic diagram of relevant geometric parameters of liquid lens.

其中,θ1,θ2表示接触角.另外,透镜的理论长度d可以表示为

其中,A是初始时刻圆形透镜的面积.液体透镜的高度可以解析地表示为

该问题的初始条件设置如下: 半径R=100 的圆形透镜坐落在Nz ×Nr=500×250 计算区域内,其中心点的坐标为 (zc,rc)=(250, 0),从而序参数的初始分布可以设定为

其他相关的物理参数给定为ρ1=10,ρ2=1,ρ3=5,界面厚度D=4,运动学黏度ν=0.1,松弛时间τ1=τ2=0.8,迁移率M0=0.01.边界条件采用与静态的双液滴测试相同的边界处理格式.图4 给出了三种不同表面张力比σ12:σ13:σ23=0.6:0.6:1,1:1:1, 1:4/3:1 时,液体透镜到达系统平衡状态具有的界面形状.从图中可以发现,液体透镜在不同的表面张力比时具有不同稳态界面形状,其界面形状与文献报道的数值结果[12]定性上一致.进一步测量了液体透镜在平衡时的接触角θ1和θ2,并将数值预测的计算结果与相应的解析解总结于表1中.从表1 可以发现,LB 方法所预测接触角的数值结果与解析解的最大相对误差不超过1.7%,即数值结果与理论解相符合.本文还计算了透镜的长度d和高度 (h1,h2) 以及相应的解析解,并将结果列于表2 中.从表2 可以发现,当表面张力比为(0.6:0.6:1) 时,长度d的数值解与解析解存在着较大误差,而当表面张力比取其他情形,数值预测的结果与解析解相比大体上可以接受.这种误差较高的结果也存在于前人的轴对称或三维的LB 方法模拟[15,24]中.可以通过采用更精细的网格或更高阶的差分格式来计算梯度和拉普拉斯项以降低计算误差.

表1 不同表面张力比的液体透镜的平衡接触角θ1 和 θ2Table 1.Equilibrium contact angles θ1 and θ2 of liquid lens at different surface tension ratios.

表2 不同表面张力比液体透镜的长度d 和高度 h1,h2Table 2.Length d and heights h1,h2 of liquid lens at different surface tensions.

图4 不同界面张力比液体透镜的平衡形态 (a) σ12 :σ13 :σ23=0.6:0.6:1 ;(b) σ12 :σ13 :σ23=1:1:1 ;(c) σ12 :σ13 :σ23=1:4/3:1.Fig.4.Equilibrium morphology of liquid lens at different interfacial tension ratios: (a) σ12 :σ13 :σ23=0.6:0.6:1 ;(b) σ12 :σ13 :σ23=1:1:1 ;(c) σ12 :σ13 :σ23=1:4/3:1.

3.3 Rayleigh-Plateau 不稳定性

Rayleigh-Plateau 不稳定性问题广泛存在于日常生活和实际工程中,如水龙头中流出的射流、喷墨打印和基因芯片排列等.在不稳定性的演化过程中,其界面动力学行为呈现出轴对称的性质,因此,本文提出的轴对称多相流LB 模型适合用于模拟经典的Rayleigh-Plateau 不稳定性问题.根据文献调研可知[14],已有的研究多数着眼于两相情形,而本文的轴对称三相流模型LB 模型理论上也满足减少一致性条件,因而,我们首先模拟二元流体的Rayleigh-Plateau 不稳定性问题来验证数值模型的正确性.计算的网格分辨率设定为Nz×Nr=λ×200,其中λ表示波长.序参数的初始分布设置如下:

其中,R=60 为液柱的半径,d=0.1Rcos(2πz/λ)为两相界面处施加的微小扰动.为了考察波数k(k=2πR/λ) 对不稳定性界面动力学行为的影响,本文考虑三种不同的波长,分别设定为λ=600,1200 和1800,其他物理参数为ρ1=10,ρ3=1,σ13=0.1,D=4,M0=0.1,τ1=τ2=1.0,ν=0.1.另外选择毛细时间为特征时间,且如下的演化时间均被特征值无量纲化.图5 给出了不同波数下Rayleigh-Plateau 不稳定性的界面动力学行为.从图中可以发现,界面处不稳定性的扰动随着时间逐渐发展,导致中间区域的液体线程变得越来越纤细,而两端的液体尺寸逐渐增大.经过一段时间的演化,液柱界面在两个夹点发生破裂,形成了一段细长的液体灯丝和一个离散的主液滴.接下来,不稳定性在不同的波数下显示出不同的界面动力学特征.当波数k=0.63 或0.31 时,液体灯丝在表面张力的作用下逐渐收缩,并最终形成一个卫星液滴.而当波数较小k=0.21 时,液体灯丝最初也随时间而逐渐收缩,并伴随着灯丝两侧的体积逐渐增大,从而在界面表面又形成了二次Rayleigh-Plateau 不稳定性,导致了液体线程发生二次破裂,最终释放出三个卫星液滴.本文计算获得不同波数下Rayleigh-Plateau 不稳定性的界面动力学行为与前人的研究结果[14,25]在定性上是相符合的.

图5 不同波数下两相Rayleigh-Plateau 不稳定性的界面演化过 程 (a) k=0.63 ;(b) k=0.31 ;(c) k=0.21.Fig.5.Interfacial evolution processes of two-phase Rayleigh-Plateau instability with different wavenumbers:(a) k=0.63 ;(b) k=0.31 ;(c) k=0.21.

接下来采用本文提出的轴对称三相流LB 模型研究三相流体的Rayleigh-plateau 不稳定性问题,重点分析波数和内部液体与中间液体半径比对三相流体界面动力学行为、界面破裂时间及生成子液滴尺寸的影响.根据三相系统中流体所处的位置,将三种流体分别称为内部液体、中间液体和外部液体.计算网格设定为Nz ×Nr=λ×2 00,边界条件与上一个算例相同.序参数的初始分布给定如下:

3.3.1 波数的影响

首先考察波数对三相流体Rayleigh-Plateau不稳定性演化的影响,并通过调节初始扰动的波长λ来获得不同的波数k=2πRI/λ,其中波长λ的值分别为600,900,1200,1500,1800.图6给出了不同波数下三相Rayleigh-Plateau 不稳定性的界面演化过程.从图中可以观察到,不稳定性在初始阶段表现出相似的界面动力学特征: 界面随时间发生变形,其扰动的振幅逐渐增大,导致中间区域的液体线程变得越来越细薄,而两端的液体体积逐渐增大.另外可以发现,内部液体的表面扰动也随时间逐渐增长,并经历了更大的界面变形,从而在波数较大(k=0.42 或0.21)时,优先发生界面破裂形成了液体灯丝和主液滴结构,而在波数较小(k=0.14)时,在界面表面上形成了多模态的扰动.随后,当波数k=0.42 时,内部的液体灯丝在表面张力作用下首先逐渐收缩,在界面表面上形成了二次扰动,并在二次Rayleigh-Plateau 不稳定性的作用下,形成的液体灯丝会再次发生断裂,分离出两个卫星液滴.卫星液滴随着中间液体的断裂,分别向位于两端的主液滴移动并与主液滴的内部流体发生合并,从而最终形成了位于两端的复合主液滴和由中间液体断裂形成的卫星液滴,并且复合主液滴表现为中间液体包裹内部液体的形态.当波数k=0.21 时,内部的液体灯丝会在多个位置发生二次破裂,从而形成了更多数量的卫星液滴.紧接着,两测的卫星液滴向主液滴方向移动,与主液滴中内部流体接触而合并,最终形成了复合的主液滴.而中间液体界面断裂后形成的液体灯丝受表面张力作用逐渐收缩,致使余下的卫星液滴向中间聚拢并融合成单一液滴,并最终形成了复合的卫星液滴.当波数充分小(k=0.14)时,多模态的扰动随后逐渐发展,并使中间流体和内部液体的灯丝在多个位置发生断裂,最终生成了一个复合主液滴和三个复合的卫星液滴,这表明随着波数的减小,三相流体Rayleigh-Plateau 不稳定性可以诱导生成更多的复合子液滴,这与二元Rayleigh-Plateau 不稳定性的结果[26]在定性上是一致的.另外还统计了波数对界面破裂时间的影响,发现不同波数(k=0.42,0.28, 0.21, 0.17, 0.14)下内部流体线程破裂的无量纲时间分别为10.6,15.4,20.0,23.8,26.3,而相应的中间流体线程的无量纲破裂时间为21.9,27.5,34.4,38.8,40.0,这表明在三相Rayleigh-Plateau不稳定性的演化中,内部流体优先于中间流体发生界面破裂,并且内外流体界面的破裂时间均随着波数的减小而增加.

图6 不同波数下三相Rayleigh-Plateau 不稳定性的界面演化过程 (a) k=0.42 ;(b) k=0.21 ;(c) k=0.14Fig.6.Interfacial evolution processes of three-phase Rayleigh-Plateau instability with different wavenumbers:(a) k=0.42 ;(b) k=0.21 ;(c) k=0.14.

进一步,研究了波数对复合液体线程破裂生成子液滴尺寸的影响,统计结果如图7 所示,其中R*表示生成子液滴的无量纲半径,其值由生成的复合液滴的半径与内部液体的初始半径RI的比值得到.由于在低波数情形,不稳定性可以诱导生成多个卫星液滴,因此本文仅考虑系统中尺寸最大的复合主液滴和尺寸最小的卫星液滴.从图7 可以发现,复合主液滴和卫星液滴的尺寸均随着波数的增大呈现出先增长而后减少的趋势.这是由于在波数非常小时,不稳定性的波长较大,液体线程发生破裂时产生了更多数目的子液滴,从而根据液体体积守恒,导致生成复合主液滴和卫星液滴的尺寸较小.当波数较大时,液体线程发生破裂生成一个主液滴和卫星液滴,其尺寸随着波数的增加而减小,该结果与二元流体Rayleigh-Plateau 不稳定性的结果[14,25,26]一致.

图7 不同波数下复合液体线程破裂生成主液滴和卫星液滴的尺寸Fig.7.The sizes of main and satellite droplets originating from the breakup of compound liquid threads with different wavenumbers.

3.3.2 内外液柱半径比的影响

本节研究内部液体与中间液体半径比对三相流体Rayleigh-Plateau 不稳定性演化的影响.计算区域设置为Nz ×Nr=1200×200,其中内部液体的半径RI的值固定为40,而中间液体的半径RM的值分别取 60,70,80,90,100.图8 给出了不同内外液柱半径比下三相流体Rayleigh-Plateau 不稳定性的界面演化过程.从图中可以发现,三元Rayleigh-Plateau 不稳定性在不同的液柱半径比下表现出不同的相界面动力学过程.当液柱半径比为RM/RI=1.5 时,不稳定性的界面扰动在初始阶段随时间逐渐发展,导致中间区域的线程变薄而两端部分不断增大.内部流体经历更大的界面变形,从而优先于中间流体发生破裂,形成了主液滴和液体灯丝.紧接着,液体灯丝的末尾也发生断裂释放出两个卫星液滴.卫星液滴向主液滴方向移动,与主液滴发生合并.与此同时,中间流体的线程也随时间的演化而越来越薄,并发生破裂形成了复合的主液滴和液体灯丝.在复合液体灯丝收缩的过程中,相界面处形成了微小的扰动,进而发生二次不稳定性现象,导致相界面发生破裂生成了三个卫星液滴.由于惯性作用,两侧的复合卫星液滴逐渐靠近中间位置的复合卫星液滴,在表面张力的作用下最终合并为一个体积更大的复合卫星液滴.当增大液柱半径比RM/RI至2 或2.5 时,可以发现内部流体仍能优先于中间流体发生破裂,形成了主液滴和液体灯丝.然而,不同于低液柱半径比情形,液体灯丝随后逐渐收缩,受到中间流体的界面变形的影响,发生破裂形成了被中间流体包裹的多个子液滴.最外侧的两个子液滴向两侧运动,并与主液滴接触后发生融合.紧接着,中间流体的界面也发生破裂,形成了复合的主液滴和包裹着多个子液滴的复合液体灯丝.在表面张力作用下,复合液滴灯丝逐渐收缩,多个子液滴同时也向中间位置移动并发生合并,最终形成了一个复合的卫星液滴.此外,还研究了内外液柱半径比对三相Rayleigh-Plateau 不稳定性界面破裂时间的影响,结果表明内部流体线程在不同半径比下(RM/RI=1.5, 1.75,2.0, 2.25, 2.5)发生破裂的无量纲时间分别为 21.0,20.5,20.0,19.4,18.1,而中间流体线程破裂所对应的无量纲时间为26.3,30.6,34.4,39.4,44.0.可以发现内部液体线程仍然比中间液体线程更早发生破裂,并且增大液柱半径比可以促进内部液体线程的破裂,而延缓中间液体线程发生破裂,这是由于中间流体在较大的液柱半径比下占比更多的体积,需要更大的界面变形直至破裂.

图8 不同内外液柱半径比下三相Rayleigh-Plateau 不稳定性的界面演化过程 (a) RM/RI=1.5 ;(b) RM/RI=2 ;(c) RM/RI=2.5Fig.8.Interfacial evolution of three-phase Rayleigh-Plateau instability with different radius ratios of the inner and outer liquid columns: (a) RM/RI=1.5 ;(b) RM/RI=2 ;(c) RM/RI=2.5.

最后,统计了不同液柱半径比下,液体线程破裂生成复合的主液滴和卫星液滴的尺寸,结果如图9 所示.从图中可以发现,复合主液滴的尺寸随着液柱半径比的增加而增大,而复合卫星液滴的尺寸随着液柱半径比的增大变化不显著,这表明中间流体体积增大的部分随着不稳定性界面的变形均往复合的主液滴方向运动.

图9 不同液柱半径比下复合液体线程破裂生成复合主液滴和卫星液滴的尺寸Fig.9.The sizes of the compound main and satellite droplets originating from the breakup of compound liquid threads with different radius ratios.

4 结论

本文基于多组分Cahn-Hilliard 相场理论,建立了一类简单高效准确的轴对称三相流LB 模型.该模型利用了三个粒子分布函数: 两个序参数分布函数追踪三相流体的相界面,另一个密度分布函数求解流体的速度场和压力场分布.通过引入修正的平衡态分布函数和合适的源项,本文模型可以正确地描述轴对称三相流体系统,并且不同于前人的轴对称三相LB 模型,本文模型中由轴对称效应引起的源项不包含任何复杂的梯度项,从而在算法实现上更加简单高效.为了验证所提出的轴对称三相流LB 模型,模拟了一系列经典的轴对称多相流算例,包括静态的双液滴、液体透镜的扩展和二元流体Rayleigh-Plateau 不稳定性.数值结果显示本文模型可以准确地捕捉三相流体的相界面,预测的不同界面张力比下液体透镜的稳态形状和Rayleigh-Plateau 不稳定性中界面动力学过程与解析解或文献报道数据一致.进一步采用提出的LB 算法研究了三相流体Rayleigh-Plateau 不稳定性问题,详细分析了波数和液柱半径比对不稳定性的界面演化过程、线程破裂时间以及子液滴尺寸的影响.可以发现随着波数的减少,不稳定性可以诱导生成更多的复合子液滴.不同于两相Rayleigh-Plateau 不稳定性情形,液体线程破裂形成子液滴的尺寸随着波数的增加先增大而后减少.另外,增大液柱半径比可以促进复合主液滴尺寸的增长,而对复合卫星液滴尺寸的影响较小.