下载:
-
三相流体系统广泛存在于自然界和工程应用中, 如水污染物的输运、强化采油中油水气的传输过程、微流控装置中复合液滴的生成技术等. 由于涉及多种流体界面的迁徙、变形、破裂及合并等现象, 三相流体的界面动力学行为十分复杂. 随着计算机技术的快速发展, 数值模拟已成为研究三相流体流动问题的有效手段. 它可以方便地提供实验研究中不易测量的物理量, 如在相界面变化过程中某一时刻的界面形状、界面处流体的速度、压力及密度分布等. 根据不同的数值方法, 一些学者对三相流体系统进行了研究, 并在笛卡尔坐标系下建立了多种数值模型, 如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不稳定性的增长过程, 详细分析了波数和内外液柱半径比对复合液体线程的界面动力学行为、界面破裂时间及生成子液滴尺寸的影响机制.
-
本文模型建立在一个含有三种不互溶流体组成的不可压缩系统, 并引入
$ \phi_1 $ ,$ \phi_2 $ ,$ \phi_3 $ 三个序参数来描述三种流体在混合物所占的体积分数, 且满足如下的约束条件:$ \sum\limits_{{{i}} = 1}^3 {{\phi _i} = 1},\; \; \; 0 \leqslant {\phi _i} \leqslant 1. $ 根据多组分流体的相场理论[17,18], 三相流体系统的总的自由能泛函可以表述为
$ \varPsi = \int\nolimits_\varOmega \bigg[\frac{{12}}{D}F({\phi _1},{\phi _2},{\phi _3}) + \sum\limits_{i = 1}^3 {\frac{3}{8}D{\lambda _i}|} \nabla {\phi _i}| \bigg] {\rm{d}}\varOmega, $ 其中,
$\varOmega $ 为三相系统所处的流域,$ F({{\phi _1}, \; {\phi _2}, \; {\phi _3}}) $ 为体相区的自由能,$ D $ 为界面厚度的特征尺度,$ \lambda_i $ 为与表面张力有关的物理参数. 在相场理论[19]中, 化学势定义为自由能泛函对序参数的变分. 然而, 针对三相流体系统, 为了满足序参数的守恒条件, 需要在自由能的变分项中引入一个拉格朗日乘数 β, 从而化学势$ \mu_i $ 可以表述为$ \begin{split} &{\mu _i} = \frac{{\partial \varPsi }}{{\partial {\phi _i}}} + \beta = \frac{{12}}{D}\frac{{\partial F}}{{\partial {\phi _i}}} - \frac{3}{4}D{\lambda _i}{\nabla ^2}{\phi _i} + \beta,\; \;\\ & \qquad ~~~~~~~~i = 1,2,3. \end{split} $ 序参数
$\phi_i$ 的时间演化由流体速度的对流和化学势的扩散来描述, 其宏观控制方程可以写成多组分Cahn-Hilliard方程的形式:$ \frac{{\partial {\phi_i}}}{{\partial t}} + \nabla\cdot({\phi _i}{\boldsymbol{u}}) = \nabla \cdot ({M_i}\nabla {\mu _i}),\; \; i = 1,2,3, $ 其中, u为流体速度,
$ M_i $ 为迁移率. 假设迁移率满足约束条件$ M_1\lambda_1=M_2\lambda_2=M_3\lambda_3=M_0 $ , 并对方程(4)中所有i进行求和, 可以得到关于变量S的方程:$ \begin{split} & \frac{{\partial S}}{{\partial t}} + \nabla \cdot (S{\boldsymbol{u}}) = \\ & \nabla \cdot {M_0}\nabla \bigg( - \frac{3}{4}D{\nabla ^2}S + \sum\limits_{i = 1}^3 {\frac{{12}}{{D{\lambda _i}}}\frac{{\partial F}}{{\partial {\phi _i}}} + \beta \sum\limits_{i = 1}^3 {\frac{1}{{{\lambda _i}}}} } \bigg), \end{split} $ 其中
$ S=\phi_1+\phi_2+\phi_3 $ . 根据方程(1)给出的守恒条件,$ S=1 $ 显然是方程(5)的解, 可以推导出$ \beta = - \sum\limits_{i = 1}^3 {\frac{{4{\lambda _{\rm{T}}}}}{{D{\lambda _i}}}\frac{{\partial F}}{{\partial {\phi _i}}}}, $ 其中,
$ \lambda_{\rm{T}} $ 定义为$ \frac{3}{{{\lambda _{\rm{T}}}}} = \sum\limits_{i = 1}^3 {\frac{1}{{{\lambda _i}}}}. $ 将 β 的表达式(6)代入方程(3), 可以得到化学势
$ \mu_i $ 的表达式为$ \begin{split} &{\mu _i} = \frac{{4{\lambda _{\rm{T}}}}}{D}\sum\limits_{j \ne i} {\left[\frac{1}{{{\lambda _j}}} \bigg(\frac{{\partial F}}{{\partial {\phi _i}}} - \frac{{\partial F}}{{\partial {\phi _j}}}\bigg)\right] - \frac{3}{4}D{\lambda _i}{\nabla ^2}{\phi _i}},\; \; \;\\ &~~~~~~~~~~~~ i = 1,2,3. \\[-10pt]\end{split} $ 为了使相场模型满足减少一致性条件以及具有很好的适定性属性, Boyer等[17,18]理论构造了如下的体相区自由能:
$ \begin{split} F({\phi _1},{\phi _2},{\phi _3}) = \;&\frac{{{\lambda _1}}}{2}\phi _1^2{(1 - {\phi _1})^2}+ \frac{{{\lambda _2}}}{2}\phi _2^2{(1 - {\phi _2})^2} \\ &+ \frac{{{\lambda _3}}}{2}\phi _3^2{(1 - {\phi _3})^2},\\[-15pt] \end{split} $ 其中,
$ \lambda_i $ 是与表面张力有关的参数,$ \begin{split} &{\lambda _1} = {\sigma _{12}} + {\sigma _{13}} - {\sigma _{23}},\\ &{\lambda _2} = {\sigma _{12}} + {\sigma _{23}} - {\sigma _{13}},\\ &{\lambda _3} = {\sigma _{13}} + {\sigma _{23}} - {\sigma _{12}}, \end{split} $ 其中,
$ \sigma_{12} $ ,$ \sigma _{13} $ ,$ \sigma _{23} $ 分别表示三相体系中两种流体之间的表面张力. 将(9)式代入(8)式, 最终可以获得化学势的表达式:$ \begin{split} {\mu _i} = \;&\frac{{12}}{D}[{\lambda _i}{\phi _i}(1 - {\phi _i})(1 - 2{\phi _i})\\ &- 2{\lambda _{\rm{T}}}{\phi _1}{\phi _2}(1 - {\phi _1} - {\phi _2})] - \frac{3}{4}D{\lambda _i}{\nabla ^2}{\phi _i},\; \; \\ &i = 1,2, \\[-10pt]\end{split} $ 其中,
$ \mu_3 $ 的值可以由方程$\displaystyle\sum\nolimits_{i = 1}^3 {\frac{{{\mu _i}}}{{{\lambda _i}}} = 0}$ 确定. 为了描述流体流动, 多组分Cahn-Hillard相场方程(4)需要与不可压Navier-Stokes方程相耦合,$ \begin{split}& \nabla \cdot {\boldsymbol{u}} = 0,\\ &\frac{{\partial (\rho {\boldsymbol{u}}})}{{\partial t}} + \nabla \cdot (\rho {\boldsymbol{u}}{\boldsymbol{u}}) \\ =\;& - \nabla p+ \nabla \cdot [\rho \nu (\nabla {\boldsymbol{u}} + \nabla {{\boldsymbol{u}}^{\rm{T}}})] + {{\boldsymbol{F}}_{\rm{s}}} + \boldsymbol{G}, \end{split} $ 其中, p是压力, ρ是流体密度, ν是运动学黏性系数, G是外力,
$ {\boldsymbol{F}}_{\rm{s}} $ 是表面张力. 另外, 为了减少相界面处的虚假速度, 本文选用势能形式的表面张力${{\boldsymbol{F}}_{\rm{s}}}=\displaystyle\sum\nolimits_{i = 1}^3 {{\mu _i}\nabla {\phi _i}}$ [12].为了获得轴对称三相流体动力学的宏观控制方程, 进行了从笛卡尔坐标系到柱坐标系的变换
$ \begin{array}{l} (x,y,z) \to (r,\theta,z), \end{array} $ 其中,
$ x=r{\rm{cos}}\theta $ ,$ y=r{\rm{sin}}\theta $ , r, z, θ分别表示径向坐标, 轴向坐标和方位角. 经过一系列的代数计算, 可以得到柱坐标系下的宏观控制方程:$ {\partial _t}{\phi _i} + {\partial _\alpha }({\phi _i}{u_\alpha }) = {\partial _\alpha }({M_i}{\partial _\alpha }{\tilde \mu _i}) + \frac{{{M_i}{\partial _r}{{\tilde \mu }_i}}}{r} - \frac{{{\phi _i}{u_r}}}{r}, $ 和
$ \begin{split} &{\partial _\alpha }{u_\alpha } + \frac{{{u_r}}}{r} = 0,\\ &{\partial _t}(\rho {u_\beta }) + {\partial _\alpha }(\rho {u_\alpha }{u_\beta }) = - {\partial _\beta }p\\ &+ {\partial _\alpha }[\nu \rho ({\partial _\alpha }{u_\beta } + {\partial _\beta }{u_\alpha })] + {{\tilde F}_{{\rm{s}}\beta }} + {G_\beta } \\ &+ \frac{{\nu \rho ({\partial _r}{u_\beta } + {\partial _\beta }{u_r})}}{r} - \frac{{2\rho \nu {u_r}{\delta _{\beta r}}}}{{{r^2}}} - \frac{{\rho {u_r}{u_\beta }}}{r}, \end{split} $ 式中,
$ \alpha, \beta =\{ r, z\} $ , δ是Kronecker函数,$ {\tilde \mu _i} $ 经过坐标变换后修正的化学势${\tilde \mu _i} = \mu_i - \dfrac{3}{4}D{\lambda _i}{\partial _r}{\phi _i}/r$ , 表面张力${\tilde F_{{\rm{s}}\beta }} = \displaystyle\sum\nolimits_{i = 1}^3 {{{\tilde \mu }_i}\nabla {\phi _i}}$ . 本文考虑的多相流是无涡流, 因此上述方程均已省略了与方位角速度相关的项. 另外需要指出的是, 现有三相流的轴对称LB模型[15]是基于上述宏观控制方程而提出的, 这会导致模型中外力项包含了许多复杂的梯度项, 增加了模型实现的复杂性. 我们尝试把径向坐标r吸收至偏导数运算符合内, 通过数学上等价变形, 可以得到形式上更加简洁的三相流体的宏观控制方程:$ {\partial _t}(r{\phi _i}) + {\partial _\alpha }(r{\phi _i}{u_\alpha }) = {\partial _\alpha }({M_i}{\partial _\alpha }r{\tilde \mu _i}) - {M_i}{\partial _r}{\tilde \mu _i}, $ 和
$ \begin{split} &{\partial _\alpha }(r{u_\alpha }) = 0,\\ &{\partial _t}(r\rho {u_\beta }) + {\partial _\alpha }(r\rho {u_\alpha }{u_\beta }) = - {\partial _\beta }(rp) \\ &+ {\partial _\alpha }[r\rho \nu ({\partial _\alpha }{u_\beta }) + ({\partial _\beta }{u_\alpha })] + r({{\tilde F}_{{\rm{s}}\beta }} + {G_\beta }) \\ &+ (p - \frac{{2\rho \nu }}{r}{u_r}){\delta _{\beta r}}. \end{split} $ -
本文需要利用两个粒子分布函数的LB演化方程来捕捉三相流体间的相界面, 另一个粒子分布函数的LB方程来求解流场动力学特性. 采用最简单的单松弛碰撞模型, 追踪三相流体相界面的LB演化方程可以表述为[12]
$ \begin{split}& f_k^i({\boldsymbol{x}} + {{\boldsymbol{e}}_k}{\delta _t},t + {\delta _t}) - f_k^i({\boldsymbol{x}},t) \\ =\;& - \frac{1}{{{\tau _i}}}[f_k^i({\boldsymbol{x}},t) - f_k^{i,{\rm{eq}}}({\boldsymbol{x}},t)] \\ &+ \Big(1 - \frac{1}{{2{\tau _i}}}\Big) {\delta _t}F_k^i({\boldsymbol{x}},t),\; \; \; i=1,\; 2, \end{split} $ 其中,
$ f_k^i({\boldsymbol{x}}, t) $ 是粒子在$ t $ 时刻在$ {\boldsymbol{x}} $ 位置的序参数分布函数,$ f_k^{i, {\rm{eq}}}({\boldsymbol{x}}, t) $ 为对应的平衡态分布函数,$ \tau_i $ 是无量纲松弛时间,$ \delta_t $ 是时间增量,$F_k^i({\boldsymbol{x}}, t)$ 是外力项分布函数. 将(16)式的最右端项吸收至对流项, 从而巧妙设计了序参数的平衡态分布函数:$ \begin{split} &f_k^{i,{\rm{eq}}}({\boldsymbol{x}},t) = \\ &\left\{ \begin{aligned} &r\phi_i + r({\omega_k} - 1)\eta{{\tilde \mu }_i}, & { k=0 } \; \\ &{\omega_k}r\eta{{\tilde \mu }_i} + {\omega _k}\frac{{{e_{k\alpha }}r{\phi _i}{u_\alpha } + {e_{kr}}{M_i}{{\tilde \mu }_i}}}{{c_{\rm{s}}^2}}, & { k\neq0 }, \end{aligned} \right. \end{split} $ 其中,
$ \omega _k $ 和$ {\boldsymbol{e}}_k $ 分别是权重系数和离散速度,$ c_{\rm{s}} $ 是声速, η是引入的一个自由参数, 用来调整迁移率系数. 采用最具有代表性的$ D2 Q9 $ 离散速度模型[20], 其对应权系数的取值为$ {\omega_0}=4/9 $ ,${\omega_{1-4}}= 1/9$ ,$ {\omega_{5-8}}=1/36 $ ,$ {c_{\rm{s}}}=c/\sqrt 3 $ , 离散速度$ {\boldsymbol{e}}_k $ 定义为$ {{\boldsymbol{e}}_k} = \left\{ {\begin{aligned} & (0,0)c,\; \;& k=0\;\;\;\;\;\;\\ &(\cos [(k - 1)\pi /2],\sin [(k - 1)\pi /2])c,\; \;& k=1—4\\ &\sqrt 2 (\cos [(k - 5)\pi /2+\pi/4],\sin [(k - 5)\pi /2+\pi/4])c,\; \;& k=5—8 \end{aligned}}, \right. $ 式中
$ c={\delta_x/\delta_t} $ 表示格子速度,$ \delta_x $ 代表格子步长. 为了能够准确地恢复轴对称的多组分Cahn-Hilliard 方程, 还构造了全新的外力项分布函数:$ F_k^i = \frac{{{\omega _k}{e_{k\alpha }}{\partial _t}(r{\phi _i}{u_\alpha } + {M_i}{{\tilde \mu }_i}{\delta _{\alpha r}})}}{{c_{\rm{s}}^2}}. $ 在本模型中, 通过计算粒子分布函数的零阶矩, 来获得序参数
$ \phi_i $ 的宏观量$ {\phi _i} = \frac{1}{r}\sum\limits_k {f_k^i}. $ 另外, 序参数
$ \phi_3 $ 的值可以通过守恒条件(1)得到, 即$ \phi_3=1-\phi_1-\phi_2 $ , 并且流体的密度$ \rho $ 可以看成序参数$ \phi_i $ 的线性插值函数:$ \begin{array}{l} \rho = {\phi _1}{\rho _1} + {\phi _2}{\rho _2} + {\phi _3}{\rho _3}, \end{array} $ 其中,
$ \rho_1 $ ,$ \rho_2 $ ,$ \rho_3 $ 分别表示三相流体系统所对应的相密度. 利用Chapman-Enskog多尺度分析技术[21], 可以证明本模型可以准确地恢复到三相流体界面追踪的轴对称Cahn-Hilliard相场方程以及迁移率$ M_i $ 和松弛因子之间的关系:$ \begin{array}{l} {M_i} =\eta c_{\rm{s}}^2({\tau _i} - 0.5). \end{array} $ 由(24)式可知, 当迁移率非常小的情况下, 无量纲松弛时间
$ \tau_i $ 将会接近 0.5, 这可能会导致LB方法的数值不稳定, 因此引入自由参数η是非常有必要的, 且不会增加任何额外的计算成本.对于三相不可压缩流体的建模, 需要将多组分Cahn-Hilliard相场方程与Navier-Stoke流动方程相耦合. 基于单松弛碰撞模型, 描述流体动力学的的LB演化方程可以写成[12]
$ \begin{split}& {g_k}({\boldsymbol{x}} + {{\boldsymbol{e}}_k}{\delta_t},t + {\delta _t}) - {g_k}({\boldsymbol{x}},t) \\ =\;& - \frac{1}{{{\tau _g}}}[{g_k}({\boldsymbol{x}},t) - g_k^{{\rm{eq}}}({\boldsymbol{x}},t)] + \left(1 - \frac{1}{{2{\tau _g}}}\right){\delta _t}{G_k}({\boldsymbol{x}},t), \end{split} $ 其中,
$ {g_k}({\boldsymbol{x}}, t) $ 为描述流场的粒子分布函数,$ g_k^{{\rm{eq}}}({\boldsymbol{x}}, t) $ 为平衡态分布函数,$ \tau_g $ 为与流体黏度有关的松弛因子,$ {G_k}({\boldsymbol{x}}, t) $ 为外力项分布函数. 考虑到轴对称效应带来的影响, 为了正确的恢复连续性方程, 本文采用了如下修正的平衡态分布函数[14]:$ g_k^{{\rm{eq}}} = \left\{ {\begin{aligned} &\frac{{r\rho }}{{c_{\rm{s}}^2}}({\omega _k} - 1) + r\rho {s_k}({\boldsymbol{u}}),\; \; \;& k = 0\;\\ &\frac{{r\rho }}{{c_{\rm{s}}^2}}{\omega _k} + r\rho {s_k}({\boldsymbol{u}}),\; \; \;& k \ne 0, \end{aligned}} \right. $ 其中,
${s_k}({\boldsymbol{u}}) = {\omega _k}\left[\frac{{{{\boldsymbol{e}}_k} \cdot {\boldsymbol{u}}}}{{c_{\rm{s}}^2}} + \frac{{{{({{\boldsymbol{e}}_k} \cdot {\boldsymbol{u}})}^2}}}{{2c_{\rm{s}}^4}} - \frac{{{\boldsymbol{u}} \cdot {\boldsymbol{u}}}}{{2c_s^2}}\right]. $ 另外, 需要考虑将外力项引入介观LB方法中所产生的格子离散效应. Liang等[22]近期考虑了格子离散效应, 设计了一类准确并且简化的外力项分布函数. 在轴对称坐标系下, 流场的外力项分布函数可以表示为
$ {G_k} = {\omega _k}\left[\frac{{{e_{k\alpha }}{F_\alpha }}}{{c_{\rm{s}}^2}} + \frac{{{u_\alpha }{\partial _\beta }(r\rho ){e_{k\alpha }}{e_{k\beta }}}}{{c_{\rm{s}}^2}} - \rho {u_r}\right], $ 其中,
$ F_\alpha $ 定义为$ {F_\alpha } = r({\tilde F_{{\rm{s}}\alpha }} + {G_\alpha }) + \left(p - \frac{{2\rho \nu }}{r}{u_r}\right){\delta _{\alpha r}}. $ 通过计算粒子分布函数的一阶矩和零阶矩, 可以获得流体流动的速度
$ {u_\alpha } = \frac{{\displaystyle\sum\nolimits_k {{e_{k\alpha }}{g_k} + 0.5{\delta _t}[r({{\tilde F}_{{\rm{s}}\alpha }} + {G_\alpha }) + p{\delta _{\alpha r}}]} }}{{r\rho + {\delta _t}{r^{ - 1}}\rho \nu {\delta _{\alpha r}}}}, $ 和流体动力学压力
$ p $ 为$\begin{split} p = \;&\frac{{c_s^2}}{{1 - {\omega _0}}}\Bigg[\frac{1}{r}\sum\limits_{k \ne 0} {g_k} + \frac{{{\delta _t}}}{2}{\boldsymbol{u}} \cdot \nabla \rho \\ &+ \rho {s_0}({\boldsymbol{u}}) - \frac{{{\omega _0}\rho \nu {u_r}}}{{rc_{\rm{s}}^2}} \Bigg].\end{split} $ 利用多尺度理论分析手段[22], 可以证明本模型可以恢复正确的轴对称不可压Navier-Stokes方程, 并且运动性黏性与松弛因子之间的关系为
$ \begin{array}{l} \nu = c_{\rm{s}}^2({\tau _g} - 0.5){\delta _t}. \end{array} $ 在数值模拟中, 需要采用合适的差分格式对模型中时间导数项和空间导数项进行离散运算, 本模型采用显式的欧拉差分格式计算时间导数项,
$ {\partial _t}\chi ({\boldsymbol{x}},t) = \frac{{\chi ({\boldsymbol{x}},t) - \chi ({\boldsymbol{x}},t - {\delta _t})}}{{{\delta _t}}}, $ 以及各向同性的二阶中心差分格式计算梯度算子和拉普拉斯算子,
$ \begin{split} &\nabla \chi ({\boldsymbol{x}},t) = \sum\limits_{k \ne 0} {\frac{{{\omega _k}{{\boldsymbol{e}}_k}\chi ({\boldsymbol{x}} + {{\boldsymbol{e}}_k}{\delta _t},t)}}{{c_{\rm{s}}^2{\delta _t}}},} \\ &{\nabla ^2}\chi ({\boldsymbol{x}},t) = \sum\limits_{k \ne 0} {\frac{{2{\omega_k}[\chi ({\boldsymbol{x}} + {{\boldsymbol{e}}_k}{\delta_t},t) - \chi ({\boldsymbol{x}},t)]}}{{c_{\rm{s}}^2\delta_t^2}}}, \end{split} $ 其中χ表示任意变量. 另外, 为了避免轴对称LB模型由于涉及
$ r^{-1} $ 的项而在$ r=0 $ 处产生奇异现象, 将第一层网格放在$ r=0.5\delta_x $ 处, 并对轴向边界处应用对称边界条件[14]. -
本节首先通过模拟几个典型的轴对称多相流问题来验证本文提出的介观LB模型, 包括静态的双液滴、液体透镜和二元流体Rayleigh-Plateau不稳定性问题. 接着, 利用该模型来研究三相流体的Rayleigh-Plateau不稳定性问题, 详细分析波数和内外液柱的半径比对液线破裂过程中界面运动过程以及产生的复合液滴尺寸的影响.
-
静态的双液滴作为三相流的基准算例, 首先被用来验证本文所建立的轴对称三相流 LB 模型. 在静态的双液滴测试中, 将计算网格设置为
${N_z}\times {N_r}=400\times100$ , 其中尺寸相同的两个液滴静止放置于计算区域中, 液滴的半径为$ R=50 $ , 两液滴的中心坐标分别是$ ({z_1}, {r_0}) = (100, 0) $ 和$ ({z_2}, r_0) = (300, 0) $ . 初始的序参数分布设定为$ \begin{split}& {\phi _1}(z,r)= \\ & 0.5 + 0.5\tanh \left[2\frac{{R - \sqrt {{{(z - {z_1})}^2} + {{(r - {r_0})}^2}} }}{D}\right],\\ &{\phi _2}(z,r)= \\ & 0.5 + 0.5\tanh \left[2\frac{{R - \sqrt {{{(z - {z_2})}^2} + {{(r - {r_0})}^2}} }}{D}\right], \end{split} $ 而序参数
$ \phi_3 $ 可以由方程(1)得到, 其他物理参数设置为: 三相系统密度$ \rho_1=10 $ ,$ \rho_2=5 $ ,$ \rho_3=1 $ , 界面厚度$ D=4 $ , 运动学黏性$ \nu=0.1 $ , 界面张力$ \sigma_{12}=\sigma_{13}=\sigma_{23}=0.1 $ , 松弛时间$ \tau_1=\tau_2=0.8 $ , 迁移率$ M_0=0.01 $ . 粒子分布函数$ f_i $ 和$ g_i $ 的边界条件设置如下: 左右边界运用周期边界条件, 上壁面采用无滑移反弹边界条件, 下边界则采用对称性边界条件. 图1 给出了稳态时三个序参数在计算区域中的分布. 可以看出, 序参数的稳态分布符合初始条件给定的双曲分布, 避免了三相流的数值方法可能会产生非物理的虚假相. 进一步统计了序参数沿着界面方向的稳态分布以及相应的解析解, 并将结果列于图2. 从图中可以观察到数值预测的序参数分布与解析解能很好地符合, 这表明本文提出轴对称三相流的LB模型能够准确捕获界面的分布.
图 1 静态双液滴测试中三个序参数的稳态分布 (a)
$ \phi_1 $ 的分布; (b)$ \phi_2 $ 的分布; (c)$ \phi_3 $ 的分布Figure 1. Steady-state distributions of three order parameters in the static double-droplet test: (a) Distribution of
$ \phi_1 $ ; (b) distribution of$ \phi_2 $ ; (c) distribution of$ \phi_3 $ 
图 2 静态双液滴测试中序参数沿水平线的稳态分布
Figure 2. Steady distributions of the order parameters across the horizontal line in static double-droplet test
-
液体透镜的扩展问题被广泛用于测试所提出的三相流LB模型, 然而绝大多数工作限制于二维笛卡尔坐标系[9,10,12], 针对三维或柱坐标空间的模拟鲜有报道. 本节将模拟三维液体透镜的扩展来验证所提出的三相流轴对称LB模型. 如图3所示, 圆形的液体透镜初始时刻放置于另外两种不混溶的流体之间, 在表面张力作用下会发生变形, 并最终达到稳定的形状. 根据Neumann定律[23], 液体透镜在平衡时的形状由三个表面张力系数决定:
图 3 液体透镜相关几何参数示意图
Figure 3. Schematic diagram of relevant geometric parameters of liquid lens
$ \begin{split} &\cos ({\theta _1}) = \frac{{\sigma _{23}^2 + \sigma _{12}^2 - \sigma _{13}^2}}{{2{\sigma _{23}}{\sigma _{12}}}},\; \; \;\\ &\cos ({\theta _2}) = \frac{{\sigma _{23}^2 + \sigma _{13}^2 - \sigma _{12}^2}}{{2{\sigma _{23}}{\sigma _{13}}}}, \end{split}$ 其中,
$ \theta_1 $ ,$ \theta_2 $ 表示接触角. 另外, 透镜的理论长度d可以表示为$ {\left(\frac{d}{2}\right)^2}\sum\limits_{j = 1}^2 {\frac{1}{{\sin ({\theta _j})}}\left[\frac{{{\theta _j}}}{{\sin ({\theta _j})}} - \cos ({\theta _j})\right] = A,} $ 其中, A 是初始时刻圆形透镜的面积. 液体透镜的高度可以解析地表示为
$ {h_j} = \frac{{1 - \cos ({\theta _j})}}{{\sin ({\theta _j})}}\left(\frac{d}{2}\right),\; \; \; j = 1,2. $ 该问题的初始条件设置如下: 半径
$ R=100 $ 的圆形透镜坐落在$ {N_z}\times{N_r}=500\times250 $ 计算区域内, 其中心点的坐标为$(z_{\rm{c}}, \; r_{\rm{c}})=(250, \; 0)$ , 从而序参数的初始分布可以设定为$\begin{split}& {\phi _1}(z,r)= \\ & 0.5 + 0.5\tanh \Bigg[2\frac{{R - \sqrt {{{(z -{z_{\rm{c}}})}^2} + {{(r - {r_{\rm{c}}})}^2}} }}{D}\;\Bigg],\\ &{\phi _2}(z,r)= \\ & \max \left\{ \left[0.5 + 0.5\tanh \frac{{2(z - {z_{\rm{c}}})}}{D}\right] - {\phi _1},0\right\}, \\[-18pt]\end{split} $ 其他相关的物理参数给定为
$ \rho_1=10 $ ,$ \rho_2=1 $ ,$ \rho_3= 5 $ , 界面厚度D = 4, 运动学黏度$ \nu=0.1 $ , 松弛时间$ \tau_1=\tau_2=0.8 $ , 迁移率$ M_0=0.01 $ . 边界条件采用与静态的双液滴测试相同的边界处理格式. 图4给出了三种不同表面张力比$\sigma_{12}:\sigma_{13}:\sigma_{23}=0.6:0.6: 1, 1:1:1, \; 1:4/3:1$ 时, 液体透镜到达系统平衡状态具有的界面形状. 从图中可以发现, 液体透镜在不同的表面张力比时具有不同稳态界面形状, 其界面形状与文献报道的数值结果[12]定性上一致. 进一步测量了液体透镜在平衡时的接触角$ \theta_1 $ 和$ \theta_2 $ , 并将数值预测的计算结果与相应的解析解总结于表1中. 从表1可以发现, LB方法所预测接触角的数值结果与解析解的最大相对误差不超过1.7%, 即数值结果与理论解相符合. 本文还计算了透镜的长度d和高度$ (h_1, h_2) $ 以及相应的解析解, 并将结果列于表2中. 从表2可以发现, 当表面张力比为$ (0.6:0.6:1) $ 时, 长度d的数值解与解析解存在着较大误差, 而当表面张力比取其他情形, 数值预测的结果与解析解相比大体上可以接受. 这种误差较高的结果也存在于前人的轴对称或三维的LB方法模拟[15,24]中. 可以通过采用更精细的网格或更高阶的差分格式来计算梯度和拉普拉斯项以降低计算误差.
图 4 不同界面张力比液体透镜的平衡形态 (a)
$\sigma_{12}:\sigma_{13}: $ $ \sigma_{23}=0.6:0.6:1$ ; (b)$ \sigma_{12}:\sigma_{13}:\sigma_{23}=1:1:1 $ ; (c)$\sigma_{12}: $ $ \sigma_{13}:\sigma_{23}=1:4/3:1$ .Figure 4. Equilibrium morphology of liquid lens at different interfacial tension ratios: (a)
$ \sigma_{12}:\sigma_{13}:\sigma_{23}=0.6:0.6:1 $ ; (b)$ \sigma_{12}:\sigma_{13}:\sigma_{23}=1:1:1 $ ; (c)$\sigma_{12}:\sigma_{13}:\sigma_{23}= 1:4/ $ $ 3:1$ .表面张力比
$(\sigma_{12}:\sigma_{13}:\sigma_{23})$数值解 解析解 相对误差 $\theta_1$ $\theta_2$ $\theta_1$ $\theta_2$ $\theta_1$ $\theta_2$ $0.6:0.6:1$ 34.0 34.0 33.6 33.6 1.2% 1.2% $1:1:1$ 60.1 60.1 60 60 0.2% 0.2% $1:4/3:1$ 84.4 49.0 83.6 48.2 1.0% 1.7% 表 1 不同表面张力比的液体透镜的平衡接触角
$\theta_1$ 和$\theta_2$ Table 1. Equilibrium contact angles
$\theta_1$ and$\theta_2$ of liquid lens at different surface tension ratios表面张力比$(\sigma_{12}:\sigma_{13}:\sigma_{23})$ 数值解 解析解 相对误差 d $h_1$ $h_2$ d $h_1$ $h_2$ d $h_1$ $h_2$ $0.6:0.6:1$ 322.1 49.2 49.2 391.8 59.1 59.1 17.8% 16.8% 16.8% $1:1:1$ 253.0 73.2 73.2 277.0 80.0 80.0 8.7% 8.5% 8.5% $1:4/3:1$ 234.5 106.1 53.4 251.6 112.5 56.3 6.8% 5.7% 5.2% 表 2 不同表面张力比液体透镜的长度d和高度
$h_1$ ,$h_2$ Table 2. Length d and heights
$h_1$ ,$h_2$ of liquid lens at different surface tensions -
Rayleigh-Plateau不稳定性问题广泛存在于日常生活和实际工程中, 如水龙头中流出的射流、喷墨打印和基因芯片排列等. 在不稳定性的演化过程中, 其界面动力学行为呈现出轴对称的性质, 因此, 本文提出的轴对称多相流LB模型适合用于模拟经典的Rayleigh-Plateau不稳定性问题. 根据文献调研可知[14], 已有的研究多数着眼于两相情形, 而本文的轴对称三相流模型LB模型理论上也满足减少一致性条件, 因而, 我们首先模拟二元流体的Rayleigh-Plateau不稳定性问题来验证数值模型的正确性. 计算的网格分辨率设定为
${N_{\rm{z}}} \times {N_r} = \lambda \times 200$ , 其中$ \lambda $ 表示波长. 序参数的初始分布设置如下:$ \begin{split} &{\phi _1}(z,r) = 0.5 + 0.5\tanh \frac{{2(R + d - r)}}{D},\\ &{\phi _2}(z,r) = 0,\\ &{\phi _3}(z,r) = 1 - {\phi _1}(z,r) - {\phi _2}(z,r). \end{split} $ 其中,
$ R=60 $ 为液柱的半径,$ d=0.1 R\cos(2\pi z/\lambda ) $ 为两相界面处施加的微小扰动. 为了考察波数$ k $ ($ k=2\pi R/\lambda $ ) 对不稳定性界面动力学行为的影响, 本文考虑三种不同的波长, 分别设定为$ \lambda=600 $ , 1200和1800, 其他物理参数为$ \rho_1=10 $ ,$ \rho_3=1 $ ,$ \sigma_{13}=0.1 $ ,$ D=4 $ ,$ M_0=0.1 $ ,$ \tau_1=\tau_2=1.0 $ ,$ \nu=0.1 $ . 另外选择毛细时间$ \sqrt {{R^3}{\rho _1}/{\sigma _{13}}} $ 为特征时间, 且如下的演化时间均被特征值无量纲化. 图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 $ .Figure 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不稳定性问题, 重点分析波数和内部液体与中间液体半径比对三相流体界面动力学行为、界面破裂时间及生成子液滴尺寸的影响. 根据三相系统中流体所处的位置, 将三种流体分别称为内部液体、中间液体和外部液体. 计算网格设定为
$ {N_{{z}}}\times{N_r}=\lambda\times 200 $ , 边界条件与上一个算例相同. 序参数的初始分布给定如下:$ \begin{split}& {\phi _1}(z,r) \\ =\;& 0.5 + 0.5\tanh \frac{{2[RI +0.1RM\cos (2\pi z/\lambda ) - r]}}{D},\\ &{\phi _2}(z,r) \\ =\;& 0.5 + 0.5\tanh \frac{{2[RM +0.1RM\cos (2\pi z/\lambda ) - r]}}{D}\\ &-\phi_1(z,r),\\[-15pt] \end{split} $ 其中,
$ RI $ ,$ RM $ 表示内部液体和中间液体的初始半径, 分别取$ RI=40 $ ,$ RM=80 $ . 其他物理参数为$ \rho_1=\rho_2=\rho_3=1 $ ,$ \sigma_{12}=\sigma_{13}=\sigma_{23}=0.1 $ ,$D= 4$ ,$ \tau_1=\tau_2=0.8 $ ,$ \nu=0.1 $ ,$ M_0=0.01 $ . 另外, 选取内部液柱的半径$ RI $ 为特征长度和毛细时间$ \sqrt {{{RI}^3}{\rho_1}/{\sigma_{13}}} $ 为特征时间, 且如下的统计量均被这些特征值无量纲化. -
首先考察波数对三相流体Rayleigh-Plateau不稳定性演化的影响, 并通过调节初始扰动的波长λ来获得不同的波数
$ k=2\pi RI/\lambda $ , 其中波长$ \lambda $ 的值分别为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.14 $ Figure 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 不同波数下复合液体线程破裂生成主液滴和卫星液滴的尺寸
Figure 7. The sizes of main and satellite droplets originating from the breakup of compound liquid threads with different wavenumbers
-
本节研究内部液体与中间液体半径比对三相流体Rayleigh-Plateau不稳定性演化的影响. 计算区域设置为
$ {N_{{z}}}\times {N_r}=1200\times 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.5 $ Figure 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 不同液柱半径比下复合液体线程破裂生成复合主液滴和卫星液滴的尺寸
Figure 9. The sizes of the compound main and satellite droplets originating from the breakup of compound liquid threads with different radius ratios
-
本文基于多组分Cahn-Hilliard相场理论, 建立了一类简单高效准确的轴对称三相流LB模型. 该模型利用了三个粒子分布函数: 两个序参数分布函数追踪三相流体的相界面, 另一个密度分布函数求解流体的速度场和压力场分布. 通过引入修正的平衡态分布函数和合适的源项, 本文模型可以正确地描述轴对称三相流体系统, 并且不同于前人的轴对称三相LB模型, 本文模型中由轴对称效应引起的源项不包含任何复杂的梯度项, 从而在算法实现上更加简单高效. 为了验证所提出的轴对称三相流LB模型, 模拟了一系列经典的轴对称多相流算例, 包括静态的双液滴、液体透镜的扩展和二元流体Rayleigh-Plateau不稳定性. 数值结果显示本文模型可以准确地捕捉三相流体的相界面, 预测的不同界面张力比下液体透镜的稳态形状和Rayleigh-Plateau不稳定性中界面动力学过程与解析解或文献报道数据一致. 进一步采用提出的LB算法研究了三相流体Rayleigh-Plateau不稳定性问题, 详细分析了波数和液柱半径比对不稳定性的界面演化过程、线程破裂时间以及子液滴尺寸的影响. 可以发现随着波数的减少, 不稳定性可以诱导生成更多的复合子液滴. 不同于两相Rayleigh-Plateau不稳定性情形, 液体线程破裂形成子液滴的尺寸随着波数的增加先增大而后减少. 另外, 增大液柱半径比可以促进复合主液滴尺寸的增长, 而对复合卫星液滴尺寸的影响较小.
-
基于多组分相场理论提出了一类模拟三相流体流动的轴对称格子Boltzmann模型. 该模型利用两个粒子分布函数来捕捉三种不同流体之间的相界面, 另一个粒子分布函数来求解流体动力学方程以获得流场信息. 为了刻画坐标变换引起的轴对称效应, 巧妙地设计了演化方程中平衡态分布函数和外力项分布函数, 从理论上保证本文模型可以正确恢复三相流体系统的宏观控制方程, 并且轴对称效应产生的源项中不包含任何复杂的梯度项, 从而比现有的轴对称格子Boltzmann模型更加简单高效. 首先通过模拟一系列轴对称多相流的基准算例, 包括静态的双液滴、液体透镜的扩展和二元流体Rayleigh-Plateau不稳定性, 来验证本文模型的有效性与正确性. 接下来, 利用该模型研究了三相流体的Rayleigh-Plateau不稳定性的增长过程, 定量分析了波数和液柱半径比对复合液体线程破裂过程中界面动力学行为、界面破裂时间以及生成子液滴尺寸的影响. 可以发现复合的液体线程在波数较大时破裂生成一个复合主液滴和卫星液滴, 而在波数较小时可以生成更多数量的卫星液滴, 这导致复合主液滴和卫星液滴的尺寸随着波数的增加呈现先增大而后减少的趋势. 另外, 我们发现内部流体优先于中间流体发生界面破裂, 并且流体界面的破裂时间均随着波数的减小而增加. 最后发现, 增大液柱半径比可以促进内部液体线程的破裂, 而延缓中间液体线程发生破裂, 并且复合主液滴的尺寸随着液柱半径比的增加而增大, 而复合卫星液滴的尺寸对液柱半径比的变化不显著.
-
关键词:
- 格子Boltzmann方法 /
- 轴对称 /
- 三相流 /
- Rayleigh-Plateau不稳定性
Based on the multi-component phase field theory, in this paper we propose an axisymmetric lattice Boltzmann model for three-phase fluids. The proposed model takes advantage of two particle distribution functions for capturing phase interface among three different fluids, and another particle distribution function for solving the hydrodynamic equations for flow field. In order to describe the axisymmetric effect arising from the coordinate transformation, we elaborately design the equilibrium distribution function and forcing distribution function in the evolution equation, which ensures that the model can accurately recover the macroscopic governing equation for three-phase fluids. Also, the introduced source terms accounting for the axisymmetric effect contain no additional gradient term, which makes it be simpler than the existing lattice Boltzmann model for axisymmetric three-phase fluids. To validate the proposed model, a series of axisymmetric multiphase benchmark examples are performed, including the static double droplets, the spreading of liquid lens, and the binary-fluid Rayleigh-Plateau instability. It is reported that the present model can accurately capture the phase interface, and the predicted steady shapes of the liquid lens agree well with the analytical profiles. Then, the proposed model is used to study the three-phase Rayleigh-Plateau instability and the effects of the wavenumber and the radius ratio of liquid column on the interfacial dynamic behaviour, the breakup time of liquid threads and the size of daughter droplet are investigated in detail. It can be found that the compound liquid thread at a high wavenumber could break up into one main droplet and one satellite droplet, but the multiple satellite droplets can be produced at a low wavenumber, which leads to that the sizes of main and satellite droplets increase with the wavenumber at first and then decrease with it. Besides, we can observe that the inner fluid undergoes the breakup at earlier time than the middle fluid, and the breakup time for both inner and middle fluids increases with the decrease of the wavenumber. Finally, we can find that increasing the radius ratio of liquid column accelerates the breakup of inner-fluid thread, but prevents the breakup of the middle-fluid thread. In addition, the size of the compound main droplet increases with the radius ratio of liquid column, while the size of the compound satellite droplet doest not change much with it.[1] Smith K A, Solis F J, Chopp D L 2002 Interfaces Free Bound. 4 263
[2] Bonhomme R, Magnaudet J, Duval F, Piar B 2012 J. Fluid Mech. 707 405
[3] Kim J 2007 Comput. Methods Appl. Mech. Eng. 196 4779
Google Scholar
[4] Guo Z L, Shu C 2013 Lattice Boltzmann Method and Its Applications in Engineering (Singapore: World Scientific) pp10–32
[5] 胡晓亮, 梁宏, 王会利 2020 物理学报 69 044701
Google Scholar
Hu X L, Liang H, Wang H L 2020 Acta Phys. Sin. 69 044701
Google Scholar
[6] 贺传晖, 刘高洁, 娄钦 2021 物理学报 70 244701
Google Scholar
He C H, Liu G J, Lou Q 2021 Acta Phys. Sin. 70 244701
Google Scholar
[7] Spencer T J, Halliday I, Care M C 2010 Phys. Rev. E 82 066701
Google Scholar
[8] Leclaire S, Reggio M, Trepanier J 2013 J. Comput. Phys. 246 318
Google Scholar
[9] Yu Y, Liu H H, Liang D, Zhang Y H 2019 Phys. Fluids 31 012108
Google Scholar
[10] Semprebon C, Kruger T, Kusumaatmaja H 2016 Phys. Rev. E 93 033305
Google Scholar
[11] Wöhrwag M, Semprebon C, Moqaddam M A, Karlin I 2018 Phys. Rev. Lett. 120 234501
Google Scholar
[12] Liang H, Shi B C, Chai Z H 2016 Phys. Rev. E 93 013308
Google Scholar
[13] Liang H, Xu J R, Chen J X, Chai Z H, Shi B C 2019 Appl. Math. Model. 73 487
Google Scholar
[14] Liang H, Chai Z H, Shi B C, Guo Z L, Zhang T 2014 Phys. Rev. E 90 063311
Google Scholar
[15] Haghani-Hassan-Abadi R, Rahimian M 2020 Acta Mech. 231 2323
Google Scholar
[16] Yang J X, Li Y B, Lee C Y, Kim J 2021 Eur. J. Mech.-B/Fluid 89 203
Google Scholar
[17] Boyer F, Lapuerta C 2006 ESAIM: Math. Model. Numer. Anal. 40 653
Google Scholar
[18] Boyer F, Lapuerta C, Minjeaud S, Piar B, Quintard M 2010 Transp. Porous Media 82 463
Google Scholar
[19] Jacqmin D 1999 J. Comput. Phys. 155 96
Google Scholar
[20] Qian Y, d’Humieres D, Lallemand P 1992 Europhys. Lett. 17 479
Google Scholar
[21] Liang H, Shi B C, Guo Z L, Chai Z H 2014 Phys. Rev. E 89 053320
Google Scholar
[22] Liang H, Li Y, Chen J X, Xu J R 2019 Int. J. Heat Mass Transfer 130 1189
Google Scholar
[23] Rawlinson J S, Widom B 1982 The Molecular Theory of Capillarity (Oxford: Clarendon Press) pp208–214
[24] Jiang F, Tsuji T 2017 Water Resour. Res. 53 11
Google Scholar
[25] Zheng L, Zheng S, Zhai Q L 2021 Comput. Fluid 218 104857
Google Scholar
[26] Tjahjadi M, Stone H A, Ottino J M 1992 J. Fluid Mech. 243 297
Google Scholar
-
图 1 静态双液滴测试中三个序参数的稳态分布 (a)
$ \phi_1 $ 的分布; (b)$ \phi_2 $ 的分布; (c)$ \phi_3 $ 的分布Fig. 1. Steady-state distributions of three order parameters in the static double-droplet test: (a) Distribution of
$ \phi_1 $ ; (b) distribution of$ \phi_2 $ ; (c) distribution of$ \phi_3 $ 
图 2 静态双液滴测试中序参数沿水平线的稳态分布
Fig. 2. Steady distributions of the order parameters across the horizontal line in static double-droplet test
图 3 液体透镜相关几何参数示意图
Fig. 3. Schematic diagram of relevant geometric parameters of liquid lens
图 4 不同界面张力比液体透镜的平衡形态 (a)
$\sigma_{12}:\sigma_{13}: $ $ \sigma_{23}=0.6:0.6:1$ ; (b)$ \sigma_{12}:\sigma_{13}:\sigma_{23}=1:1:1 $ ; (c)$\sigma_{12}: $ $ \sigma_{13}:\sigma_{23}=1:4/3:1$ .Fig. 4. Equilibrium morphology of liquid lens at different interfacial tension ratios: (a)
$ \sigma_{12}:\sigma_{13}:\sigma_{23}=0.6:0.6:1 $ ; (b)$ \sigma_{12}:\sigma_{13}:\sigma_{23}=1:1:1 $ ; (c)$\sigma_{12}:\sigma_{13}:\sigma_{23}= 1:4/ $ $ 3:1$ .
图 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 $ 
图 6 不同波数下三相Rayleigh-Plateau不稳定性的界面演化过程 (a)
$ k=0.42 $ ; (b)$ k=0.21 $ ; (c)$ k=0.14 $ Fig. 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 不同波数下复合液体线程破裂生成主液滴和卫星液滴的尺寸
Fig. 7. The sizes of main and satellite droplets originating from the breakup of compound liquid threads with different wavenumbers
图 8 不同内外液柱半径比下三相Rayleigh-Plateau不稳定性的界面演化过程 (a)
$ RM/RI=1.5 $ ; (b)$RM/RI= $ $ 2$ ; (c)$ RM/RI=2.5 $ Fig. 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 不同液柱半径比下复合液体线程破裂生成复合主液滴和卫星液滴的尺寸
Fig. 9. The sizes of the compound main and satellite droplets originating from the breakup of compound liquid threads with different radius ratios
表 1 不同表面张力比的液体透镜的平衡接触角
$\theta_1$ 和$\theta_2$ Table 1. Equilibrium contact angles
$\theta_1$ and$\theta_2$ of liquid lens at different surface tension ratios表面张力比
$(\sigma_{12}:\sigma_{13}:\sigma_{23})$数值解 解析解 相对误差 $\theta_1$ $\theta_2$ $\theta_1$ $\theta_2$ $\theta_1$ $\theta_2$ $0.6:0.6:1$ 34.0 34.0 33.6 33.6 1.2% 1.2% $1:1:1$ 60.1 60.1 60 60 0.2% 0.2% $1:4/3:1$ 84.4 49.0 83.6 48.2 1.0% 1.7% 
下载: 导出CSV
表 2 不同表面张力比液体透镜的长度d和高度
$h_1$ ,$h_2$ Table 2. Length d and heights
$h_1$ ,$h_2$ of liquid lens at different surface tensions表面张力比$(\sigma_{12}:\sigma_{13}:\sigma_{23})$ 数值解 解析解 相对误差 d $h_1$ $h_2$ d $h_1$ $h_2$ d $h_1$ $h_2$ $0.6:0.6:1$ 322.1 49.2 49.2 391.8 59.1 59.1 17.8% 16.8% 16.8% $1:1:1$ 253.0 73.2 73.2 277.0 80.0 80.0 8.7% 8.5% 8.5% $1:4/3:1$ 234.5 106.1 53.4 251.6 112.5 56.3 6.8% 5.7% 5.2%
下载: 导出CSV
-
[1] Smith K A, Solis F J, Chopp D L 2002 Interfaces Free Bound. 4 263
[2] Bonhomme R, Magnaudet J, Duval F, Piar B 2012 J. Fluid Mech. 707 405
[3] Kim J 2007 Comput. Methods Appl. Mech. Eng. 196 4779
Google Scholar
[4] Guo Z L, Shu C 2013 Lattice Boltzmann Method and Its Applications in Engineering (Singapore: World Scientific) pp10–32
[5] 胡晓亮, 梁宏, 王会利 2020 物理学报 69 044701
Google Scholar
Hu X L, Liang H, Wang H L 2020 Acta Phys. Sin. 69 044701
Google Scholar
[6] 贺传晖, 刘高洁, 娄钦 2021 物理学报 70 244701
Google Scholar
He C H, Liu G J, Lou Q 2021 Acta Phys. Sin. 70 244701
Google Scholar
[7] Spencer T J, Halliday I, Care M C 2010 Phys. Rev. E 82 066701
Google Scholar
[8] Leclaire S, Reggio M, Trepanier J 2013 J. Comput. Phys. 246 318
Google Scholar
[9] Yu Y, Liu H H, Liang D, Zhang Y H 2019 Phys. Fluids 31 012108
Google Scholar
[10] Semprebon C, Kruger T, Kusumaatmaja H 2016 Phys. Rev. E 93 033305
Google Scholar
[11] Wöhrwag M, Semprebon C, Moqaddam M A, Karlin I 2018 Phys. Rev. Lett. 120 234501
Google Scholar
[12] Liang H, Shi B C, Chai Z H 2016 Phys. Rev. E 93 013308
Google Scholar
[13] Liang H, Xu J R, Chen J X, Chai Z H, Shi B C 2019 Appl. Math. Model. 73 487
Google Scholar
[14] Liang H, Chai Z H, Shi B C, Guo Z L, Zhang T 2014 Phys. Rev. E 90 063311
Google Scholar
[15] Haghani-Hassan-Abadi R, Rahimian M 2020 Acta Mech. 231 2323
Google Scholar
[16] Yang J X, Li Y B, Lee C Y, Kim J 2021 Eur. J. Mech.-B/Fluid 89 203
Google Scholar
[17] Boyer F, Lapuerta C 2006 ESAIM: Math. Model. Numer. Anal. 40 653
Google Scholar
[18] Boyer F, Lapuerta C, Minjeaud S, Piar B, Quintard M 2010 Transp. Porous Media 82 463
Google Scholar
[19] Jacqmin D 1999 J. Comput. Phys. 155 96
Google Scholar
[20] Qian Y, d’Humieres D, Lallemand P 1992 Europhys. Lett. 17 479
Google Scholar
[21] Liang H, Shi B C, Guo Z L, Chai Z H 2014 Phys. Rev. E 89 053320
Google Scholar
[22] Liang H, Li Y, Chen J X, Xu J R 2019 Int. J. Heat Mass Transfer 130 1189
Google Scholar
[23] Rawlinson J S, Widom B 1982 The Molecular Theory of Capillarity (Oxford: Clarendon Press) pp208–214
[24] Jiang F, Tsuji T 2017 Water Resour. Res. 53 11
Google Scholar
[25] Zheng L, Zheng S, Zhai Q L 2021 Comput. Fluid 218 104857
Google Scholar
[26] Tjahjadi M, Stone H A, Ottino J M 1992 J. Fluid Mech. 243 297
Google Scholar
计量
- 文章访问数: 680
- PDF下载量: 21
- 被引次数: 0
