A REGULARIZED PHASE-FIELD LATTICE BOLTZMANN MODEL FOR TWO-PHASE FLOWS
-
摘要: 目前, 基于相场理论的格子玻尔兹曼 (lattice Boltzmann, LB)模型已经被广泛应用于气-液两相流流动问题中. 为了提高已有相场LB模型的数值稳定性, 提出了一种新型的正则化相场格子玻尔兹曼模型(regularized lattice Boltzmann model, RLBM)用于模拟大密度比和大黏度比的气液两相流动. 该模型由界面追踪和流场求解两个核心模块构成, 其中界面变化采用守恒型Allen-Cahn (A-C)相场方程控制, 流场演化则通过不可压Navier-Stokes (N-S)方程控制. 首先, 构建了两个正则化格子玻尔兹曼方程(lattice Boltzmann equation, LBE)分别获取流场信息和相场信息. 与标准的单松弛模型不同的是, 提出的模型在演化方程的碰撞项中引入了非平衡态的预前碰撞函数, 且该预前碰撞项仅与宏观量有关. 通过Chapman-Enskog (C-E)多尺度展开分析, 证实了所提出的模型能够准确地恢复到宏观流场控制方程和相场控制方程. 进一步地, 为了验证本模型的有效性, 模拟了4个两相流典型算例, 包括静态液滴、Rayleigh-Taylor (R-T)不稳定性问题、气泡上升和单个液滴撞击液膜. 数值结果证实了提出的模型能够准确地模拟大密度比、大黏度比、大雷诺数下的气液两相流流动问题. 更重要的是, 相较于传统的相场单松弛模型在小的迁移率$( {\theta _M} < 2.0 \times {10^{ - 2}}) $下就会诱发数值方法不收敛的问题, 提出的模型在模拟迁移率较小$( {\theta _M} = 1.0 \times {10^{ - 6}} )$的复杂两相流动时表现出更好的稳定性, 能够更准确地刻画界面流动, 捕捉界面形态.
-
关键词:
- 正则化相场模型 /
- 气-液两相流 /
- Allen-Cahn方程 /
- 格子玻尔兹曼方法
Abstract: At present, the lattice Boltzmann (LB) model based on the phase field theory has been widely applied in gas-liquid two-phase flow problems. In order to improve the numerical stability of the existing phase field LB models, a new regularized phase-field lattice Boltzmann model (RLBM) is proposed for simulating gas-liquid two-phase flows with large density ratio and high viscosity ratio in this work. The proposed model consists of two core modules, namely interface tracking and flow field solver, where the interface is governed by the conservating Allen-Cahn (A-C) phase-field equation, and the flow field is governed by the incompressible Navier-Stokes (N-S) equations. Firstly, two regularized lattice Boltzmann equations (LBE) have been constructed to obtain flow field and phase field information, respectively. Unlike the standard single-relation-time (SRT) model, a non-equilibrium pre-collision function, which is only related to the macroscopic variables, have been introduced into the collision term of the evolution equation in the proposed model. It has been confirmed that this model can accurately recover to the macroscopic flow field and phase filed governing equations by the multi-scale Chapman-Enskog (C-E) analysis. Furthermore, to verify the effectiveness of the present model, four typical two-phase flow cases were simulated in this paper, including static droplet, Rayleigh-Taylor (R-T) instability problems, bubble rising and a single droplet impacting on the liquid film. The obtained numerical results of these typical examples demonstrate that the proposed model can accurately simulate gas-liquid two-phase flow problems under large density ratio, high viscosity ratio and high Reynolds number. More importantly, compared to the traditional phase field SRT model, which could cause non-convergence issues at low mobility $( {\theta _M} < 2.0 \times {10^{ - 2}} )$, it has been found that the model proposed in this paper exhibits better stability in simulating complex two-phase flow with low mobility $( {\theta _M} = 1.0 \times {10^{ - 6}}) $, and can more accurately characterize interface flow and capture interface morphology. -
引言
两相流问题流动现象丰富, 广泛存在于自然界和工程领域, 是能源、环境、微流控和生物医学等领域重点考虑的共性科学问题. 随着高性能计算设备和计算理论的发展, 数值模拟逐渐成为了研究两相流动问题的主要方法.
两相流动中界面存在变形、破裂与融合等现象, 相与相之间存在复杂的相互作用. 准确捕捉界面变化和描述不同相之间的相互作用是两相流动数值模拟的关键问题之一[1]. 传统的两相流数值方法以离散Navier-Stokes (N-S)方程为基础, 并结合相应的界面捕捉方法(interface-capturing method)捕获流体界面变化的详细信息, 主要包括流体体积法[2](volume of fluid)、水平集法[3](level set method)和界面追踪法[4](front tracking method)等. 但这些方法都是基于宏观的控制方程进行分析, 属于宏观方法, 难以刻画不同流体之间的微观相互作用[5-6]. 格子Boltzmann方法(lattice Boltzmann method, LBM)是近30年发展起来的一种介观数值方法, 相对于传统宏观方法在求解两相流动问题有非常明显优势[7]. LBM的气体动理学背景使其不受连续介质假设的限制, 保证了该方法能够描述两相流质量及能量的传输过程, 同时兼具微观本质和粒子特性, 能够直观地刻画不同相之间的微观相互作用[8-10].
在多相流领域, 主要的LB模型有颜色梯度模型[11]、伪势模型[12]、自由能模型[13]和相场模型[6]. 其中, 基于扩散界面方法的相场模型能够通过求解界面控制方程来显式地描述界面变化, 在研究两相流动问题中表现出较好的研究前景[14]. 相场LB模型采用两套方程来分别追踪界面和演化流场. 控制流场流动的N-S方程组采用LBM求解, 描述界面位置的相场方程即Cahn-Hilliard (C-H)方程[15]或Allen-Cahn (A-C)方程[16]则可以选择LBM或是有限差分等方法求解. He等[6]在1999首次提出了一种双分布函数模型求解两相流动, 准确模拟了界面变形、破碎和卷起等复杂界面动力学行为. He等[6]和Lee等[17]提出了一种稳定的混合差分离散格式来刻画. 该模型实现了大密度比不可压缩两相流动问题的模拟. 但是, 上述两种模型所恢复的界面方程相较于C-H方程而言包含额外的人工虚假项. 为了准确地恢复到C-H方程, Zheng等[18]基于源项的思想修正了界面演化方程右端的空间作用力项, Zu等[14]将平衡态分布函数的空间作用力项引入到演化方程的右端, Liang等[19]将时间导数作为源项引入到界面LBE中. 这些行之有效的措施使得所恢复的界面方程能够和C-H方程保持一致, 并提高了模型的稳定性和界面追踪精度.
相较于包含4阶序参数导数的C-H方程, A-C方程中的扩散项只涉及序参数的2阶导数, 更容易求解且具有较好的稳定性[20-23]. Geier等[24]首次将LBM用于求解A-C方程, 并获得了比C-H方程更好的界面追踪精度. 随后Fakhari等[25-26]借鉴Lee等[17]的工作提出了一个有限差分LBE求解A-C方程, 结合多松弛格式(multi-relation-time, MRT)提高稳定性, 并应用自适应格式优化界面追踪, 成功模拟一系列大密度比多相流动问题. 但是以上这两个模型在通过C-E展开恢复到界面控制方程含有多余项, 不能够准确地恢复到A-C方程, 虽然这些多余的人工项对于一般的界面追踪来说影响不大, 但是在模拟涉及表面张力等复杂多相流动时则存在较大影响[27]. 针对这一问题, 国内两个研究团队(Ren等[28]和Wang等[29])独立地提出了两种类似LBE模型的方案来求解A-C方程, 二者的关键点都是在LBE中引入源项处理扩散项的一部分, 以精确恢复A-C方程. 不同的是Wang等[29]的单松弛 (single-relation-time, SRT)模型中空间导数项通过分布函数的非平衡部分计算得到, 且碰撞过程可以局部实现, 而Ren等[28]提出的MRT模型中空间导数的计算是基于一般的2阶差分格式. Liang等[30]结合Wang等[29]的工作提出了一个适用于模拟大密度比两相流动的单松弛LB模型, 该模型更为简洁且稳定性较好. Zu等[14]将源项基于平衡态分布函数实现, 提出了一个界面LB模型用以演化界面, 并能够准确地恢复到A-C控制方程. Xu等[31]基于MRT格式提出了一种非对角松弛矩阵用于将LBE完全恢复到守恒型A-C方程, 并进一步分析了多余的人工项对不同问题的影响, 结果表明该额外项对于简单问题特别是对小的马赫数问题影响不大, 但是对于复杂流动问题将会引起非物理扩散或界面波动不稳定.
在相场LB模型中, SRT模型凭借其形式简单、编程方便、计算高效等优势, 在多相流动问题研究中应用较为广泛. 然而, 该模型在松弛参数接近0.5时容易出现数值不稳定性问题[32]. 为此, 大多数改进型相场模型都是采用MRT格式. 显然, 相比于SRT模型, MRT模型能够在牺牲一定的计算效率下获得很好的计算精度, 但是MRT模型的参数选择没有规范标准, 往往需要依靠多次尝试来确定最佳参数, 这就给实际应用造成了不便[33]. 正则化格式是LBM体系内另一种具有良好稳定性的碰撞模型, 正则化格式的SRT模型仅引入一些只与宏观量相关的预前碰撞项, 并消除所有的高阶平衡态信息[33-36]. 目前, RLBM已经被用于求解对流扩散方程[35]、传热[37]、纳米流动[38]和电流体流动[39]等复杂问题, 数值结果表明RLBM在不引入复杂参数的条件下, 相较于传统SRT模型有更好的数值准确度, 并改善了模型的稳定性. 最近, RLBM在多相流领域也取得了一些进展. Ba等[40]将RLBM应用到颜色梯度模型中, 其模拟结果表明该模型得到的数值结果的准确性和稳定性能够媲美MRT模型获得的结果, 同时计算效率更高. 杨旭光等[41]基于由Peng-Robinson自由能模型, 发展了一个正则化自由能两相流模型, 获得了较好的稳定性. Otomo等[42]则在伪势模型的基础上实现了对多组分多相流动的模拟, 其模型在稳定性和鲁棒性方面都较原伪势模型更好. 但是这些模型都还缺乏模拟大密度比高雷诺数两相流动的能力, 在实际的应用中存在着局限.
鉴于RLBM在颜色梯度模型以及伪势模型中表现出来的优势, 本文将正则化碰撞模型应用到相场模型中, 发展了一种适用于大密度比两相流动的正则化相场LB模型. 正则化碰撞模型数值稳定性较好, 且算法实现简单, 计算效率较高. 本文针对A-C相场方程和含有外力项的不可压N-S方程组, 分别提出了相应的RLBM模型, 并借助Chapman-Enskog (C-E)多尺度展开分析梳理了该模型与宏观控制方程的关系, 证实了该模型能够准确恢复到宏观控制方程, 最后通过相应的算例验证了模型的优势.
1. 正则化相场LB两相流模型
1.1 控制方程
基于相场理论的非混溶两相流数值仿真需求解连续性方程和动量方程以描述流体动力学行为, 并结合相场方程追踪界面变化. 流场的演化由不可压N-S方程组描述[43]
$$ \nabla \cdot {\boldsymbol{u}} = 0 $$ (1) $$ \frac{{\partial \rho {\boldsymbol{u}}}}{{\partial t}} + \nabla \cdot (\rho {\boldsymbol{uu}}) = - \nabla p + \nabla \cdot \left\{ {\mu \left[ {\nabla {\boldsymbol{u}} + {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}}} \right]} \right\} + {\boldsymbol{F}} $$ (2) 其中$ {\boldsymbol{u}}, \rho, p $和$ \mu $分别代表流体的速度、密度、压力和动力黏度. $ {\boldsymbol{F}} $为力项, 其构成包括体积力$ {\boldsymbol{G}} $和表面张力$ {{\boldsymbol{F}}_s} $. 在本文中, 表面张力的计算采用$ {{\boldsymbol{F}}_s} = {\mu _\phi }\nabla \phi $. $\phi $为序参数, $ {\mu _\phi } $为化学势. 对两相流系统的自由能泛函做变分运算, 便可以得到化学势$ {\mu _\phi } $的表达式
$$ {\mu _\phi } = 4\beta \left(\phi - \frac{{{\phi _l} + {\phi _g}}}{2}\right)(\phi - {\phi _l})(\phi - {\phi _g}) - \kappa {\nabla ^2}\phi $$ (3) 式中, $ {\phi _l} $和$ {\phi _g} $分别是标记液相和气相的参数. $ \beta $和$ \kappa $是与界面厚度$W$和表面张力系数$\sigma $相关的参数, 二者可以通过以下关系式计算
$$ \beta \text{ = }\frac{12\sigma }{W}\text{, }\qquad \kappa = \frac{3W\sigma }{2} $$ (4) 序参数$\phi $的演化由A-C相场方程控制. 依据相场理论, 控制流体之间界面演化的A-C方程的守恒形式为[20, 24]
$$ \frac{{\partial \phi }}{{\partial t}} + \nabla \cdot (\phi {\boldsymbol{u}}) = \nabla \cdot \left[ {{\theta _M}(\nabla \phi - \lambda {\boldsymbol{n}})} \right] $$ (5) 式中, $ {\theta _M} $表示迁移率, $\lambda $为与序参数$\phi $和界面厚度$W$相关的函数
$$ \lambda = \frac{{1 - 4{{(\phi - {\phi _0})}^2}}}{W} $$ (6) 式中, ${\phi _0} = ({\phi _l} + {\phi _g})/2$表示两相之间的界面. $ {\boldsymbol{n}} $为界面处的单位梯度方向, 其定义为
$$ {\boldsymbol{n}} = \frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}} $$ (7) 1.2 正则化两相流LB模型
正则化LB模型是一类具有良好数值稳定性的碰撞模型, 其核心思想为在演化方程中引入仅与宏观量相关的预碰撞项, 且摒弃高阶非平衡态信息. 基于此思想, 针对A-C方程的正则化LB模型为
$$\begin{split} &{g}_{\alpha }(x + {{{\boldsymbol{e}}}}_{\alpha }{\delta }_{t},t + {\delta }_{t}) = {g}_{\alpha }^{{\mathrm{eq}}}({{\boldsymbol{x}}},t) + \left(1-\frac{1}{{\tau }_{g}}\right){g}_{\alpha }^{{\mathrm{neq}}}({{\boldsymbol{x}}},t) +\\ &\qquad {\delta }_{t}\left(1-\frac{1}{2{\tau }_{g}}\right){G}_{\alpha }\text{, }\quad \alpha \text{ = }0,1,2,\cdots,8\end{split} $$ (8) 其中
$$ g_\alpha ^{{\mathrm{eq}}} = {\omega _\alpha }\phi \left(1 + \frac{{{{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}}}}{{c_s^2}}\right)$$ (9) $$ g_\alpha ^{{\mathrm{neq}}}({{\boldsymbol{x}}},t) = \frac{{{\omega _\alpha }{{{\boldsymbol{e}}}_\alpha } \cdot \displaystyle\sum\nolimits_\beta {{{{\boldsymbol{e}}}_\beta }\left[ {g_\beta ^{}({{\boldsymbol{x}}},t) - g_\beta ^{{\mathrm{eq}}}({{\boldsymbol{x}}},t)} \right]} }}{{c_s^2}} $$ (10) $$ {G_\alpha } = \frac{{{\omega _\alpha }{{\boldsymbol{e}}_\alpha } \cdot [{\partial _t}(\phi {\boldsymbol{u}}) + c_s^2\lambda {\boldsymbol{n}}]}}{{c_s^2}} $$ (11) 式中, $ g_\alpha ^{}({{\boldsymbol{x}}},t) $和$ g_\alpha ^{{\mathrm{eq}}}({{\boldsymbol{x}}},t) $分别代表当前$t$时刻$ {{\boldsymbol{x}}} $位置处粒子沿着速度$ {{{\boldsymbol{e}}}_\alpha } $方向运动的序参数局部分布函数和平衡态分布函数, ${\delta _t}$为相邻演化步的时间间隔, $ {\tau _g} $表示无量纲松弛时间, $ {\omega _\alpha } $表示格子方向的权系数, $ c_s^{} $为格子声速. 在离散的格子模型, 为满足各向同性, 不同的格子速度模型中$ {\omega _\alpha } $和$ c_s^{} $的选取各不相同. 权衡计算效率和计算精度, 在二维模拟中一般采用D2Q9模型[44]. 在等温模型中, 格子声速${c_s} = c/\sqrt 3 $. 无量纲松弛参数可以通过迁移率计算获得
$$ {\theta _M} = c_s^2({\tau _g} - 0.5){\delta _t} $$ (12) 其中$c = {\delta _x}/{\delta _t}$为格子特征速度, ${\delta _x}$表示空间步长. 序参数则可以简单通过分布函数$ g_\alpha ^{}({{\boldsymbol{x}}},t) $的0阶矩计算
$$ \phi = \sum\limits_\alpha {g_\alpha ^{}} $$ (13) 在相场模型中, 密度和黏度等物理量皆可以表示成序参数的线性函数
$$ \rho = {\rho _g} + ({\rho _l} - {\rho _g})\phi ,\quad \eta = {\eta _g} + ({\eta _l} - {\eta _g})\phi$$ (14) 其中下标$l$和$g$分别表示液相和气相流体. 附录A中的C-E多尺度展开分析, 清晰地展示了本模型能够准确地恢复到A-C方程的守恒形式.
在模型中, 动量方程同样采用正则化LB模型求解, 其演化方程为
$$ \begin{split} &{f_\alpha }(x + {{{\boldsymbol{e}}}_\alpha }{\delta _t},t + {\delta _t}) = f_\alpha ^{{\mathrm{eq}}}({{\boldsymbol{x}}},t) +\\ &\qquad \left(1 - \frac{1}{{{\tau _f}}}\right)f_\alpha ^{{\mathrm{neq}}}({{\boldsymbol{x}}},t) + {\delta _t}\left(1 - \frac{1}{{2{\tau _f}}}\right){\boldsymbol{F}}_\alpha ^{}\end{split} $$ (15) 平衡态分布函数$ f_\alpha ^{{\mathrm{eq}}}({{\boldsymbol{x}}},t) $的表达式为[30]
$$ f_\alpha ^{{\mathrm{eq}}} = \left\{\begin{split} & ({\omega _\alpha } - 1)p + {\omega _\alpha }\rho \left[ {\frac{{({{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}})}}{{c_s^2}} + \frac{{{{({{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}})}^2}}}{{2c_s^4}} - \frac{{{{\left| {\boldsymbol{u}} \right|}^2}}}{{2c_s^2}}} \right],{\text{ }}\alpha {\text{ = 0}} \\ & {\omega _\alpha }\left\{ {p + \rho \left[ {\frac{{({{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}})}}{{c_s^2}} + \frac{{{{({{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}})}^2}}}{{2c_s^4}} - \frac{{{{\left| {\boldsymbol{u}} \right|}^2}}}{{2c_s^2}}} \right]} \right\},{\text{ }}\alpha {\text{ = 1,}}2,\cdots{\text{,8}} \end{split} \right. $$ (16) 分布函数的非平衡部分$ f_\alpha ^{{\mathrm{neq}}}({{\boldsymbol{x}}},t) $的表达式为
$$\begin{split} &f_\alpha ^{{\mathrm{neq}}}({{\boldsymbol{x}}},t) = \frac{{{\omega _\alpha }{{\boldsymbol{Q}}_\alpha }:\displaystyle\sum\nolimits_\beta {{{{\boldsymbol{e}}}_\beta }{{{\boldsymbol{e}}}_\beta }\left[ {f_\beta ^{}({{\boldsymbol{x}}},t) - f_\beta ^{{\mathrm{eq}}}({{\boldsymbol{x}}},t)} \right]} }}{{2c_s^4}} +\\ &\qquad \frac{{{\delta _t}}}{2}\left( {\frac{{{\omega _\alpha }{{\boldsymbol{Q}}_\alpha }:\displaystyle\sum\nolimits_\beta {{{{\boldsymbol{e}}}_\beta }{{{\boldsymbol{e}}}_\beta }{\boldsymbol{F}}_\beta ^{}} }}{{2c_s^4}} - {\boldsymbol{F}}_\alpha ^{}} \right)\end{split} $$ (17) 力项的表达式为
$$ {\boldsymbol{F}}_\alpha ^{}{\text{ = }}{\delta _t}\frac{{\left( {{{\boldsymbol{e}}_\alpha } - {\boldsymbol{u}}} \right)}}{{c_s^2}} \cdot {\varGamma _\alpha }{{\boldsymbol{F}}_s} + {\delta _t}\frac{{\left( {{{\boldsymbol{e}}_\alpha } - {\boldsymbol{u}}} \right)}}{{c_s^2}} \cdot \nabla \psi (\rho )({\omega _\alpha } - {\varGamma _\alpha }) $$ (18) 其中
$$ \psi (\rho ) = p - \rho c_s^2 ,\quad \varGamma _\alpha ^{} = {\omega _\alpha }\left[ {1 + \frac{{({{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}})}}{{c_s^2}} + \frac{{{{({{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}})}^2}}}{{2c_s^4}} - \frac{{{{\left| {\boldsymbol{u}} \right|}^2}}}{{2c_s^2}}} \right] $$ (19) 流场压力$p$和速度$ {\boldsymbol{u}} $则可以分别通过分布函数的零阶和一阶矩计算得到
$$ p = \frac{1}{{1 - {\omega _0}}}\left(\sum\limits_\alpha {f_\alpha } + 0.5{\delta _t}u \cdot \nabla \rho - {\omega _0}\rho \frac{{{u^2}}}{{2c_s^2}}\right), \;\; \alpha {\text{ = 1,}}2,\cdots{\text{,8}} $$ (20) $$ \rho {\boldsymbol{u}} = \frac{{\displaystyle\sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }f_\alpha ^{}} }}{{c_s^2}} + \frac{{{\delta _t}{{\boldsymbol{F}}_s}}}{2} $$ (21) 附录B中的C-E多尺度展开分析可以清晰地看到, 本文提出的正则化流场LB模型能够准确地恢复到不可压N-S方程组.
2. 正则化相场LB两相流模型的数值分析
2.1 静态液滴
静态液滴是验证两相流模型的基本问题, 能够检验模型界面演化和表面张力实现的准确性. 在该问题中, 一个直径为$R$的液滴浸润在${N_x} \times {N_y} = 200 \times 200$的正方形区域中心, 该区域内注满着其他流体. 在平衡状态下, 液滴内外压差$ \Delta P $的大小等于表面张力系数$\sigma $与半径$R$的比值, 即$ \Delta P = {P_{{\mathrm{in}}}} - {P_{{\mathrm{out}}}}{\text{ = }}\sigma /R $, 其中$ {P_{{\mathrm{in}}}} $和$ {P_{{\mathrm{out}}}} $分别表示平衡状态下液滴内部和外部的压力. 在模拟中, 我们采用周期性边界格式, 密度比设置为$ {\rho _l}/{\rho _g} = {100} $, 界面厚度采用$W = 4$, 液体的黏度相同, 即${\nu _l} = {\nu _g} = 0.1$.
首先, 设定迁移率$ {\theta _M} = 0.05 $, 模拟在$\sigma = 0.005$和$\sigma = 0.001$两种表面张力下不同半径静态液滴的内外压差. 图1中展示了不同表面张力下在半径$20 \leqslant R \leqslant 60$范围内液滴的内外压差随$1/R$的变化情况, 实线为理论分析值, 圆圈为通过本文发展模型模拟得到的结果. 从图中能够看到模拟结果与理论分析值的预测是一致的, 表明本文方法能够很好地表述两相流动问题的表面张力. 迁移率是影响相场模型数值稳定性的重要因素, 在传统的单松弛相场LB两相流模型中迁移率的取值受限比较大, 一般都需$ {\theta _M} \geqslant 0.02 $, 当$ {\theta _M} < 0.02 $时可能会诱发数值不稳定性问题. 为了进一步验证本文模型在迁移率方面的鲁棒性, 我们选取了不同的迁移率进行测试. 图2中展示了本模型在迁移率$ {\theta _M} =1.0 \times {10^{ - 6}} $和Liang等[30]模型在迁移率$ {\theta _M} = 1.0 \times{10^{ - 2}} $的速度场结果对比图. 从图中可以看到, Liang等[30]模型在迁移率$ {\theta _M} = 1.0 \times{10^{ - 2}} $时就会导致基于标准SRT格式相场模型计算得到的流场变得紊乱, 使得计算极不稳定. 而本文发展的模型在迁移率为$ {\theta _M} = 1.0 \times {10^{ - 6}} $这样一个很小的量级时都能够保证流场的稳定演化, 捕捉到清晰的流场信息, 表明了该模型在较小的迁移率下仍然具有较好的稳定性.
![]()
图 1 不同表面张力下内外压差随半径的变化Figure 1. Variation of pressure difference between inner and outer pressure with radius under different surface tension2.2 R-T不稳定性
静态液滴算例的流动过程较为简单, 不涉及界面变化. 在本章中, 将采用二维管道中的Rayleigh-Taylor (R-T)不稳定性这一问题来检验模型在处理复杂界面变化两相流动问题的可靠性. R-T不稳定性表述的是由密度差所诱发的不同流体间界面不稳定现象, 这种现象在天体物理学和海洋学等领域都十分常见. 该问题的物理背景是: 在一个高度和宽度分别为$4{L_0}$和${L_0}$的微管道中, 较重的液相流体置于管道的上方, 较轻的气相流体则放置在液相流体的下方, 当界面之间存在一种微小的扰动$\zeta \cos (2\text{π} x/\lambda )$时, 受到重力的影响该扰动会逐渐发展, 这两种流体将会在重力的加速下相互扩张. $\zeta = 0.1{L_0}$即为扰动的振幅, $\lambda = {L_0}$表示波长. 无量纲雷诺数${{Re}} = {L_0}\sqrt {{L_0}g} /\nu $和阿特伍德数${{At}} = ({\rho _l} - {\rho _g})/({\rho _l} + {\rho _g})$用于表征R-T不稳定性的特征, $g$表示重力加速度, $\nu $为液相流体的运动黏性系数. 参考速度的定义为${U_0} = \sqrt {{L_0}g} $. 界面初始位置${y_0}$设定为
$$ {y_0} = 2{L_0} + \zeta \cos (2\text{π} x/{L_0}) $$ (22) 为了保证序参数在扰动的界面光滑连续分布, 其初始分布设定为
$$ \phi (x,y) = 0.5 + 0.5\tanh \left[\frac{{2(y - {y_0})}}{W}\right] $$ (23) 其中, $W$为两界面的厚度. 在模拟中采用的网格为$200 \times 800$, 即${L_0} = 200$. 为更好地控制边界行为, 并模拟更真实的情况, 上下边界应用无滑移边界条件, 左右边界则被定义为周期性边界条件. 表面张力系数和界面厚度分别设置为$ \sigma = 5.0 \times {10^{ - 5}} $和$W = 4$.
图3中展示了${Re} = 3000$和${{At}} = 0.5$下$t = 0.5,{\text{ }}1.0, 1.5, 2.0,{\text{ }}2.5$和$3.0$时刻流体密度云图, 其中时间$t$是由特征时间$t = \sqrt {{L_0}/(g \cdot At)} $标准化的无量纲时间. 从图3中可以清晰地观察到R-T不稳定性的发展过程. 在开始阶段, 由于流体的运动和相互扩张作用, 在两个反向旋转涡旋产生之前重流体开始向下方扩张, 同时轻流体气泡开始向上运动; 随着时间的推移, 反向旋转的涡旋逐渐消失, 并伴随着向上的次级尖峰的产生, 这些现象展示了R-T不稳定性中流体的复杂运动和相互作用的演变. 为了更直观地说明本文提出模型的有效性和准确性, 图4给出了定量数据, 详细描绘了气泡和尖端在发展过程中的位置随时间的变化. 这些数据与已有的文献中的基准结果[28, 31, 45]进行了比较, 结果表明, 计算结果与文献中的数据[28, 31, 45]相吻合, 证明了该模型能够准确地预测和描述流体的运动和相互作用.
![]()
图 3 在At = 0.5和Re = 3000时R-T不稳定性界面随时间演化过程Figure 3. The evolution of R-T instability interface with time at At = 0.5 and Re = 3000![]()
图 4 在At = 0.5和Re = 3000时界面位置随时间变化Figure 4. The changes of the interface position with time at At = 0.5 and Re = 3000在相场模拟中, 迁移率是一个影响计算稳定和界面追踪精度的关键因素. 我们进一步分析了本模型对迁移率的敏感性, 其结果展示在图5中. 可以看到, 本模型在较小的迁移率下即$ {\theta _M} = 0.002 $都能够准确地模拟此类复杂流动问题. 此外, 相较于较大的迁移率$ {\theta _M} = 0.02 $, 较小的迁移率下$ {\theta _M} = 0.002 $界面扩散更小, 能够捕捉到更为精细的界面.
![]()
图 5 在Re = 3000, $t = 3.0$时不同迁移率的界面形态Figure 5. The interface morphology with different mobilities at Re = 3000 and $t = 3.0$2.3 气泡上升
在上面的两个算例中, 分别验证了模型在静态以及复杂流动状态下的可靠性, 但是这两个算例都是基于小密度比的工况下. 本节中, 将通过单气泡在微通道中受浮力影响上升这一算例验证模型在大密度比条件下的稳定性. 气泡上升是一种司空见惯的流体流动现象, 存在于生活的方方面面. 在物理、化学以及生物医学等多个领域, 气泡运动问题一直受到广泛的关注, 它作为基础性问题的重要性不可忽视. 气泡的上升过程不仅对液体中的声音、光线和热量等产生影响, 还能促进气体和液滴之间进行质量交换. 这个过程是浮力、表面张力和黏性力等多种力耦合作用的结果, 同时伴随着复杂的界面变化. 在初始时刻, 一个大小为${L_0} \times 2{L_0}$的方形管道底部放置了一个直径为$D = 0.5{L_0}$的圆形气泡, 气泡的位置为$(0.5{L_0},0.5{L_0})$. 管道的左右两边采用周期性边界条件处理, 上下壁面应用无滑移边界条件. 在该问题中, 较为关键的无量纲数为雷诺数${Re} = {\rho _l}D\sqrt {gD} /{\mu _l}$和Eotvos数${{Eo}} = {\rho _l}g{D^2}/\sigma $.
在模拟中液相和气相的密度比与黏度比分别设置为${\rho _l}/{\rho _g} = 1000$和${\mu _l}/{\mu _g} = 100$, 气液之间的界面厚度则设置为$W = 4$. 特征长度取作${L_0} = 300$. 迁移率设定为$ {\theta _M} = 0.005 $. 模拟中固定雷诺数${Re} = 35$, 我们分析了不同${{Eo}}$数对气泡上升的影响.
图6 ~ 图8中展示了气泡上升过程中界面形态在${{Eo}} = 10$, $50$和$ 125$下的演化情况. 当${{Eo}} = 10$时, 气泡的形状在很大程度上受到表面张力的限制作用, 导致气泡的变形程度较小. 在整个上升过程中, 气泡的形状都表现出较高的规则性, 并在后期维持一个半圆的形状. 然而, 随着${{Eo}}$数的增加, 表面张力对气泡形状的影响逐渐变得有限. 此时, 气泡的变形变得更加剧烈, 两侧各形成一个拖尾, 并随着时间的推移逐渐变得更加细长. 本文所模拟的气泡上升过程与已有文献中记录和描述的现象[46-47]高度一致, 这表明我们的模型能够准确地捕捉和模拟复杂两相流动中的关键物理现象, 验证了模拟结果的准确性和可靠性. 为了更深入地定量评估本文所提出的模型的准确性, 图9详细描绘了气泡在上升过程中质心位置的变化过程, 并与已有的文献结果进行了对比. 通过观察图9, 可以清晰地看到气泡质心位置的变化轨迹整体和时间近似线性正相关. 该图中的曲线与文献中的结果[48-49]非常接近, 整体趋势和主要特征保持高度一致, 进一步表明了我们的模型能够准确地模拟气泡上升过程中的动态变化, 即证实了本文提出的模型的可靠性.
![]()
图 6 ${{Eo}} = 10$时气泡上升过程界面形态比较Figure 6. Comparison of interface morphology during bubble rise process at Eo = 10![]()
图 7 ${{Eo}} = 50$时气泡上升过程界面形态比较Figure 7. Comparison of interface morphology during bubble rise process at Eo = 50![]()
图 8 ${{Eo}} = 125$时气泡上升过程界面形态比较Figure 8. Comparison of interface morphology during bubble rise process at Eo = 1252.4 液滴撞击液膜
最后, 将通过液滴撞击液膜这一算例来评估模型在模拟大密度比和复杂界面变化问题方面的性能. 液滴撞击液膜是一个具有重要实际意义的两相流动问题, 广泛存在于生活和工业领域, 如雨滴撞击水面、喷淋系统中的液滴撞击目标表面等. 液滴撞击液膜是一个复杂的多相流动问题, 涉及到多个物理过程的相互作用, 如表面张力和黏性力等. 因此, 这一算例广泛地用于验证模型在模拟复杂多相流动方面的能力, 通过模拟这一现象, 可以全面评估模型的性能和能力.
问题模型如图10所示, 一个直径为$D$的圆形液滴位于厚度为$h = 0.1{L_0}$的液膜上端, 并以大小为$U$的初始速度撞击液膜. 液膜位于$\left[ {0,{L_0}} \right] \times \left[ {0,3{L_0}} \right]$方形区域底部, 液滴的中心坐标为$ ({x_c},{y_c}) = (3{L_0}/2,h + D/2) $. 在本算例的计算中, 液滴的直径设为$D = 200$, 选择大小为$500 \times 1500$的格子作为计算域. 液滴撞击液膜的主要运动特征主要受密度比${\rho _l}/{\rho _g}$, 雷诺数${{Re}} = {\rho _l}UD/{\mu _l}$和韦伯数${{We}} = {\rho _l}{U^2}D/\sigma $三者影响, 其中$U$为液滴的撞击速度. 周期边界条件应用于左右边界, 上下壁面则采用无滑移边界条件. 无量纲时间设定为$t = nU/D$, 其中$n$为模拟计算步数. 其他相关参数分别设置为$U = 0.025$, $D = 200$, $h = 50$, $ {\theta _M} = 0.02 $和${{We}} = 8000$.
![]()
图 10 液滴撞击液膜问题示意图Figure 10. Schematic diagram of droplet impact on liquid film problem首先, 我们固定密度比为${\rho _l}/{\rho _g} = 1000$, 分析不同雷诺数下液滴撞击液膜的演化过程. 雷诺数在${{Re}} = 20,100$和$500$时液滴撞击液膜的演化过程呈现在图11中. 可以看到, 在雷诺数较小(${{Re}} = 20$)时, 液滴撞击液膜的过程分为接触和扩散两个阶段. 在接触阶段, 液滴与液膜相互接触并形成一个清晰的接触线. 在扩散阶段, 液滴的能量被传递到液膜上, 液滴的撞击速度在接触壁面后转化为水平方向的速度, 导致液膜向左右两边扩散. 当雷诺数增大(${{Re}} \geqslant 100$)后, 液体的黏性束缚作用减小, 出现了明显的溅射现象. 在液滴撞击液膜后, 液滴与液膜的交汇处产生射流, 这些射流随着时间的推移逐渐演化成美丽的水花. 图12为不同雷诺数下液滴撞击液膜在$t = 1.6$的瞬时状态对比. 从图中可以清晰地观察到雷诺数对液滴飞溅的影响. 有趣的是, 当雷诺数处于较低水平时, 我们可以观察到液滴飞溅的高度随着雷诺数的增大而显著增加. 然而, 当雷诺数达到一定水平后, 雷诺数变化的影响变得微乎其微. 这一现象表明, 在特定的物理条件下, 雷诺数对液滴飞溅的效应存在一个阈值. 本文所观察到的这些现象都和文献中已有的结果[50]相吻合. 图13中描绘了液滴在不同雷诺数下的铺展半径. Josserand等[51]的研究指出, 液滴的铺展半径$r/2 R$与无量纲演化时间成正比, 且与流体黏度无明显关联. 从图中可以清晰看到, 模拟数据与理论预测完全一致, 且雷诺数对结果的影响微乎其微, 这进一步证实了模型在模拟严苛条件下复杂流体流动的准确性和稳定性.
![]()
图 11 不同雷诺数下液滴撞击液膜界面的演化过程Figure 11. The evolution of droplet impact on liquid film interface at different Reynolds numbers![]()
图 12 不同雷诺数下液滴撞击液膜界面形态在$t = 1.6$时的对比Figure 12. Comparison of interface morphology for droplet impacting on liquid film at $t = 1.6$ under different Reynolds numbers![]()
图 13 液滴撞击液膜过程中无量纲半径的扩散过程Figure 13. The dimensionless radius diffusion process during droplet impacting on liquid film最后, 我们固定雷诺数为${{Re}} = 500$, 分析不同密度比对液滴飞溅现象的影响. 图14中展示了不同密度比下液滴撞击液膜的演化过程. 从图中可以看到, 当密度比较小(${\rho _l}/{\rho _g} = 10$)时, 液滴撞击液膜之后所形成的水花需要更大的动能才能够带动附近气体运动, 导致其缺乏足够的能量克服表面张力和黏性的束缚, 因此溅起水花的高度较低, 且向内弯曲, 水花臂较厚. 随着密度比的增加, 液体带动其周围气体运动所需耗费的动能逐渐减少, 因而溅起的水花高度明显增加, 水花臂也明显变薄. 在密度比为${\rho _l}/{\rho _g} = 100$和${\rho _l}/{\rho _g} = 1000$时, 液滴飞溅的整体形态就较为相似, 只是在${\rho _l}/{\rho _g} = 100$时液滴水花的端部会有明显的内卷, 在${\rho _l}/{\rho _g} = 1000$的情况下水花则向外飞溅.
![]()
图 14 不同密度比下液滴撞击液膜界面形态演化过程对比Figure 14. Comparison of interface morphology evolution processes for droplet impacting on liquid film at different density ratios3. 总结
RLBM模型在碰撞项中引入非平衡态矩, 能够显著地改善标准LB模型的计算稳定性. 本文的主要工作是基于A-C相场理论提出相应的RLBM两相流模型. 模型通过求解A-C相场方程追踪界面变化, 流场的演化则采用N-S方程控制. 通过C-E多尺度展开分析, 可验证本文提出的模型能够准确恢复到宏观控制方程. 针对静态问题, 采用静态液滴验证了模型在表述表面张力方面的准确性, 同时发现相较于传统的相场LB两相流模型在$ {\theta _M} < 2.0 \times {10^{ - 2}} $时就会诱发计算紊乱, 当前模型在很小的迁移率($ {\theta _M} = 1.0 \times {10^{ - 6}} $)下仍能够保持较好的稳定性. 进一步采用R-T不稳定性验证了本模型在较小迁移率下模拟复杂两相流动时界面变化的准确性. 最后, 采用单个气泡上升和液滴撞击液膜等具备复杂两相流动的基准算例检验了模型在模拟大密度比复杂两相流动问题时的适用范围和可靠性, 发现本文发展的模型在迁移率较小时仍能够准确地模拟各类复杂两相流动问题, 同时还具备处理大密度比大雷诺数等具备挑战性的复杂两相流动问题的能力.
附录A 守恒型A-C方程的正则化LB模型C-E多尺度分析
C-E多尺度展开分析为评判所构建的LBE方程和相对应的宏观控制方程组是否具备一致性提供了分析工具. 基于C-E多尺度展开分析理论, 将分布函数、空间导数和时间导数展开为小参数$\varepsilon $的集合, 之后将该集合式分别代入式(8)的Taylor展开式中, 归纳整理后可以得到流场分布函数和序参数分布函数在${\varepsilon ^0}$, ${\varepsilon ^1}$和${\varepsilon ^2}$尺度上的方程
$$ O({\varepsilon }^{0}): {g}_{\alpha }^{(0)} = {g}_{\alpha }^{({\mathrm{eq}})} + \left(1-\frac{1}{{\tau }_{g}}\right)\frac{{\omega }_{\alpha }{{\boldsymbol{e}}}_{\alpha }\cdot {\displaystyle \sum {}_{\beta }{{\boldsymbol{e}}}_{\beta }\left({g}_{\beta }^{(0)}-{g}_{\beta }^{({\mathrm{eq}})}\right)}}{{c}_{s}^{2}} $$ (A1) $$\begin{split} &O({\varepsilon }^{1}): {g}_{\alpha }^{(1)} + {\delta }_{t}{D}_{1}{g}_{\alpha }^{(0)} = \left(1-\frac{1}{{\tau }_{g}}\right)\frac{{\omega }_{\alpha }{{\boldsymbol{e}}}_{\alpha }\cdot {\displaystyle \sum {}_{\beta }{{\boldsymbol{e}}}_{\beta }{g}_{\beta }^{(1)}}}{{c}_{s}^{2}} + \\ &\qquad {\delta }_{t}\left(1-\frac{1}{2{\tau }_{g}}\right){G}_{\alpha }^{(1)} \end{split} $$ (A2) $$ \begin{split} &O({\varepsilon }^{2}): {g}_{\alpha }^{(2)} + {\delta }_{t}\left(\frac{\partial }{\partial {t}_{2}}{g}_{\alpha }^{(0)} + {D}_{1}{g}_{\alpha }^{(1)}\right) + \frac{{\delta }_{t}{}^{2}}{2}{D}_{1}{}^{2}{g}_{\alpha }^{(0)} = \\ &\qquad \left(1-\frac{1}{{\tau }_{g}}\right)\frac{{\omega }_{\alpha }{{\boldsymbol{e}}}_{\alpha }\cdot {\displaystyle \sum {}_{\beta }{{\boldsymbol{e}}}_{\beta }{g}_{\beta }^{(2)}}}{{c}_{s}^{2}}\end{split} $$ (A3) 其中${D_1} = \dfrac{\partial }{{\partial {t_1}}} + {{\boldsymbol{e}}_\alpha } \cdot {\nabla _{x1}}$. 将式(A2)代入式(A3)中可得
$$\begin{split} & g_\alpha ^{(2)} + {\delta _t}\frac{\partial }{{\partial {t_2}}}g_\alpha ^{(0)} + \frac{{{\delta _t}}}{2}{D_1}g_\alpha ^{(1)} + \frac{{{\delta _t}}}{2}{D_1}\left(1 - \frac{1}{{{\tau _g}}}\right)\frac{{{\omega _\alpha }{{\boldsymbol{e}}_\alpha } \cdot \displaystyle\sum {_\beta {{\boldsymbol{e}}_\beta }g_\beta ^{(1)}} }}{{c_s^2}}+ \\ &\qquad\frac{{{\delta _t}^2}}{2}{D_1}\left(1 - \frac{1}{{2{\tau _g}}}\right)G_\alpha ^{(1)}{\text{ = }}\left(1 - \frac{1}{{{\tau _g}}}\right)\frac{{{\omega _\alpha }{{\boldsymbol{e}}_\alpha } \cdot \displaystyle\sum {_\beta {{\boldsymbol{e}}_\beta }g_\beta ^{(2)}} }}{{c_s^2}} \end{split} $$ (A4) 联立式(A2)和式(A4)的0阶速度矩, 并进行简单变换, 可得${t_1}$和${t_2}$时间尺度的宏观控制方程
$$ \frac{{\partial \phi }}{{\partial {t_1}}} + {\nabla _{x1}}(\phi {\boldsymbol{u}}) = 0 $$ (A5) $$ \begin{split} &\frac{{\partial \phi }}{{\partial {t_2}}} + {\nabla _{x1}} \cdot \left[ {\left(1 - \frac{1}{{2{\tau _g}}}\right)\sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }g_\alpha ^{(1)}} } \right] +\\ &\qquad \frac{{{\delta _t}}}{2}{\nabla _{x1}} \cdot \left(1 - \frac{1}{{2{\tau _g}}}\right)\left(\frac{{\partial \phi {\boldsymbol{u}}}}{{\partial {t_0}}} + c_s^2\lambda {\boldsymbol{n}}\right) = 0 \end{split}$$ (A6) 根据式(A2), 可得
$$ \begin{split} & \sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }g_\alpha ^{(1)}} =- {\tau _g}{\delta _t}\left[ \frac{\partial }{{\partial {t_0}}}\sum\nolimits_\alpha {{{\boldsymbol{e}}_\alpha }g_\alpha ^{(0)}} + {\nabla _{x1}} \cdot \sum\nolimits_\alpha {{{\boldsymbol{e}}_\alpha }{{\boldsymbol{e}}_\alpha }g_\alpha ^{(0)}} -\right.\\ &\qquad\left. \left(1 - \frac{1}{{2{\tau _g}}}\right)\sum\nolimits_\alpha {{{\boldsymbol{e}}_\alpha }G_\alpha ^{(1)}} \right]= \\ &\qquad - {\tau _g}{\delta _t}\left[ {c_s^2{\nabla _{x1}} \cdot \phi + \frac{1}{{2{\tau _g}}}\frac{{\partial \phi {\boldsymbol{u}}}}{{\partial {t_1}}} - \left(1 - \frac{1}{{2{\tau _g}}}\right)c_s^2\lambda {\boldsymbol{n}}} \right] \end{split} $$ (A7) 将式(A7)代入式(A6), 可得
$$ \frac{{\partial \phi }}{{\partial {t_2}}} + {\nabla _{x1}} \cdot \left[ {c_s^2\left({\tau _g} - \frac{1}{2}\right){\delta _t}} \right]\left({\nabla _{x1}}\phi - \lambda {\boldsymbol{n}}\right) = 0 $$ (A8) 进一步耦合式(A8)和式(A5), 并令
$$ {\theta _M} = c_s^2\left({\tau _g} - \frac{1}{2}\right){\delta _t} $$ (A9) 显然, 本文发展的模型能够准确恢复到A-C相场方程(5).
附录B 不可压N-S方程组的正则化LB模型C-E多尺度分析
借助于C-E多尺度展开分析可完美地呈现从介观LBE恢复到宏观N-S方程的过程, 现在引用C-E多尺度展开分析检验本文所提出的针对连续性方程的RLB模型. 同样, 基于C-E多尺度展开分析理论与Taylor展开, 归纳整理后可以得到流场分布函数在${\varepsilon ^0}$, ${\varepsilon ^1}$和${\varepsilon ^2}$尺度上的方程
$$ O({\varepsilon }^{0}): {f}_{\alpha }^{(0)} = {f}_{\alpha }^{({\mathrm{eq}})} + \left(1-\frac{1}{{\tau }_{f}}\right)\frac{{\omega }_{\alpha }{{\boldsymbol{Q}}}_{\alpha }:{\displaystyle {\sum }_{\beta }{{\boldsymbol{e}}}_{\beta }{{\boldsymbol{e}}}_{\beta }\left({f}_{\beta }^{}-{f}_{\beta }^{{(\mathrm{eq})}}\right)}}{2{c}_{s}^{4}} $$ (B1) $$ \begin{split} &O({\varepsilon }^{1}): {f}_{\alpha }^{(1)} + {\delta }_{t}{D}_{1}{f}_{\alpha }^{(0)} = \frac{{\delta }_{t}}{2}{{\boldsymbol{F}}}_{\alpha }^{(1)} + \\ &\qquad \left(1-\frac{1}{{\tau }_{f}}\right)\left(\frac{{\omega }_{\alpha }{{\boldsymbol{Q}}}_{\alpha }:{\displaystyle {\sum }_{\beta }{{\boldsymbol{e}}}_{\beta }{{\boldsymbol{e}}}_{\beta }{f}_{\beta }^{(1)}}}{2{c}_{s}^{4}} + \frac{{\delta }_{t}}{2}\frac{{\omega }_{\alpha }{{\boldsymbol{Q}}}_{\alpha }:{\displaystyle {\sum }_{\beta }{{\boldsymbol{e}}}_{\beta }{{\boldsymbol{e}}}_{\beta }{{\boldsymbol{F}}}_{\beta }^{(1)}}}{2{c}_{s}^{4}}\right)\end{split} $$ (B2) $$ \begin{split} &O({\varepsilon }^{2}): {f}_{\alpha }^{(2)} + \frac{\partial }{\partial {t}_{2}}{f}_{\alpha }^{(0)} + {\delta }_{t}{D}_{1}{f}_{\alpha }^{(1)} + \frac{{\delta }_{t}{}^{2}}{2}{D}_{1}{}^{2}{f}_{\alpha }^{(0)} =\\ &\qquad \left(1-\frac{1}{{\tau }_{f}}\right)\frac{{\omega }_{\alpha }{{\boldsymbol{Q}}}_{\alpha }:{\displaystyle {\sum }_{\beta }{{\boldsymbol{e}}}_{\beta }{{\boldsymbol{e}}}_{\beta }{f}_{\beta }^{(2)}}}{2{c}_{s}^{4}}\end{split} $$ (B3) 将式(B5)代入式(B6)可得
$$ \begin{split} & f_\alpha ^{(2)} + {\delta _t}\frac{\partial }{{\partial {t_2}}}f_\alpha ^{(0)} + \frac{{{\delta _t}}}{2}{D_1}\left( {f_\alpha ^{(1)} + \frac{{{\delta _t}}}{2}{\boldsymbol{F}}_\alpha ^{(1)}} \right){\text{ }} - \\ &\qquad \left(1 - \frac{1}{{{\tau _f}}}\right) {\frac{{{\omega _\alpha }{{\boldsymbol{Q}}_\alpha }:\displaystyle\sum\nolimits_\beta {{{\boldsymbol{e}}_\beta }{{\boldsymbol{e}}_\beta }f_\beta ^{(2)}} }}{{2c_s^4}}} =- \left(1 - \frac{1}{{{\tau _f}}}\right)\frac{{{\delta _t}}}{2}{D_1}\cdot\\ &\qquad \left( {\frac{{{\omega _\alpha }{{\boldsymbol{Q}}_\alpha }:\displaystyle\sum\nolimits_\beta {{{\boldsymbol{e}}_\beta }{{\boldsymbol{e}}_\beta }f_\beta ^{(1)}} }}{{2c_s^4}} + \frac{{{\delta _t}}}{2}\frac{{{\omega _\alpha }{{\boldsymbol{Q}}_\alpha }:\displaystyle\sum\nolimits_\beta {{{\boldsymbol{e}}_\beta }{{\boldsymbol{e}}_\beta }{\boldsymbol{F}}_\beta ^{(1)}} }}{{2c_s^4}}} \right) \\ {\text{ }} \end{split} $$ (B4) 将式(B1)乘以${{\boldsymbol{e}}_\alpha }{{\boldsymbol{e}}_\alpha }$, 可以得到
$$ \sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }{{\boldsymbol{e}}_\alpha }f_\alpha ^{{(\mathrm{eq})}}} = \sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }{{\boldsymbol{e}}_\alpha }f_\alpha ^{(0)}} $$ (B5) 再将上式代入式(B1)中, 有
$$ f_\alpha ^{{(\mathrm{eq})}} = f_\alpha ^{(0)} $$ (B6) 结合式(24)和分布函数的各阶矩, 能够得到
$$ \sum\limits_\alpha {f_\alpha ^{(1)}} = - \frac{{{\delta _t}{\boldsymbol{u}} \cdot \nabla \rho }}{2} \text{, } \quad \sum\limits_\alpha {f_\alpha ^{(k)}} = 0, \;\;k > 1 $$ (B7a) $$ \sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }f_\alpha ^{(1)}} = - \frac{{{\delta _t} \cdot {{\boldsymbol{F}}_s}}}{2} \text{, }\quad \sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }f_\alpha ^{(k)}} = 0,{\text{ }}k > 1 $$ (B7b) 对式(B2)和$ {{\boldsymbol{e}}_\alpha } \cdot ({\mathrm{B2}}) $分别在$\alpha $方向上求和, 可得
$$\qquad\qquad \nabla \cdot {\boldsymbol{u}} = 0 $$ (B8a) $$\qquad\qquad \frac{{\partial \rho {\boldsymbol{u}}}}{{\partial {t_0}}} + {\nabla _{x1}} \cdot (p + \rho {\boldsymbol{uu}}) = {\boldsymbol{F}} $$ (B8b) 同理, 对式(B4)和$ {{\boldsymbol{e}}_\alpha } \cdot ({\mathrm{B4}}) $分别在$\alpha $方向上求和, 可得
$$ \begin{split} &\frac{{{\delta _t}}}{2}{D_1}\left( {f_\alpha ^{(1)} + \frac{{{\delta _t}}}{2}{\boldsymbol{F}}_\alpha ^{(1)}} \right) = \frac{{{\delta _t}}}{2}\frac{\partial }{{\partial {t_0}}}\left( { - \frac{{{\delta _t} \cdot {\boldsymbol{u}}\nabla \rho }}{2} + \frac{{{\delta _t} \cdot {\boldsymbol{u}}\nabla \rho }}{2}} \right) +\\ &\qquad \frac{{{\delta _t}}}{2}{\nabla _{x1}}\left( { - \frac{{{\delta _t} \cdot {{\boldsymbol{F}}_s}}}{2} + \frac{{{\delta _t} \cdot {{\boldsymbol{F}}_s}}}{2}} \right) = 0\end{split} $$ (B9) $$ \frac{{\partial \rho {\boldsymbol{u}}c_s^2}}{{\partial {t_2}}} + c_s^2\left(\frac{1}{2} - {\tau _f}\right){\delta _t}{\nabla _{x1}} \cdot \left\{ {\rho \left[{\nabla _{x1}} \cdot {\boldsymbol{u}} + {{\left( {{\nabla _{x1}} \cdot {\boldsymbol{u}}} \right)}^{\mathrm{T}}}\right]} \right\} = 0 $$ (B10) 令$ \mu = \rho c_s^2\left({\tau _f} - \dfrac{1}{2}\right){\delta _t} $, 并耦合不同时间尺度, 便可以得到不可压N-S方程
$$ \nabla \cdot {\boldsymbol{u}} = 0 $$ (B11a) $$ \frac{{\partial \rho {\boldsymbol{u}}}}{{\partial t}} + \nabla \cdot (\rho {\boldsymbol{uu}}) = - \nabla p + \nabla \cdot \left\{ {\mu \left[ {\nabla {\boldsymbol{u}} + {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}}} \right]} \right\} + {\boldsymbol{F}} $$ (B11b) -
![]()
图 1 不同表面张力下内外压差随半径的变化
Figure 1. Variation of pressure difference between inner and outer pressure with radius under different surface tension
![]()
图 3 在At = 0.5和Re = 3000时R-T不稳定性界面随时间演化过程
Figure 3. The evolution of R-T instability interface with time at At = 0.5 and Re = 3000
![]()
图 4 在At = 0.5和Re = 3000时界面位置随时间变化
Figure 4. The changes of the interface position with time at At = 0.5 and Re = 3000
![]()
图 5 在Re = 3000, $t = 3.0$时不同迁移率的界面形态
Figure 5. The interface morphology with different mobilities at Re = 3000 and $t = 3.0$
![]()
图 6 ${{Eo}} = 10$时气泡上升过程界面形态比较
Figure 6. Comparison of interface morphology during bubble rise process at Eo = 10
![]()
图 7 ${{Eo}} = 50$时气泡上升过程界面形态比较
Figure 7. Comparison of interface morphology during bubble rise process at Eo = 50
![]()
图 8 ${{Eo}} = 125$时气泡上升过程界面形态比较
Figure 8. Comparison of interface morphology during bubble rise process at Eo = 125
![]()
图 10 液滴撞击液膜问题示意图
Figure 10. Schematic diagram of droplet impact on liquid film problem
![]()
图 11 不同雷诺数下液滴撞击液膜界面的演化过程
Figure 11. The evolution of droplet impact on liquid film interface at different Reynolds numbers
![]()
图 12 不同雷诺数下液滴撞击液膜界面形态在$t = 1.6$时的对比
Figure 12. Comparison of interface morphology for droplet impacting on liquid film at $t = 1.6$ under different Reynolds numbers
![]()
图 13 液滴撞击液膜过程中无量纲半径的扩散过程
Figure 13. The dimensionless radius diffusion process during droplet impacting on liquid film
-
[1] Mirjalili S, Jain SS, Dodd M. Interface-capturing methods for two-phase flows: An overview and recent developments. Center for Turbulence Research Annual Research Briefs, 2017, 13: 117-135
[2] Rajendran S, Manglik RM, Jog M. New property averaging scheme for volume of fluid method for two-phase flows with large viscosity ratios. Journal of Fluids Engineering, 2022, 144(6): 061101 doi: 10.1115/1.4053548
[3] Doherty W, Phillips TN, Xie Z. A stabilised finite element framework for viscoelastic multiphase flows using a conservative level-set method. Journal of Computational Physics, 2023, 477: 111936 doi: 10.1016/j.jcp.2023.111936
[4] Inguva V, Kenig EY, Perot JB. A front-tracking method for two-phase flow simulation with no spurious currents. Journal of Computational Physics, 2022, 456: 111006 doi: 10.1016/j.jcp.2022.111006
[5] He XY, Chen SY, Doolen GD. A novel thermal model for the lattice Boltzmann method in incompressible limit. Journal of Computational Physics, 1998, 146(1): 282-300 doi: 10.1006/jcph.1998.6057
[6] He XY, Chen SY, Zhang RY. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh–Taylor instability. Journal of Computational Physics, 1999, 152(2): 642-663 doi: 10.1006/jcph.1999.6257
[7] Petersen KJ, Brinkerhoff JR. On the lattice Boltzmann method and its application to turbulent, multiphase flows of various fluids including cryogens: A review. Physics of Fluids, 2021, 33(4): 041302 doi: 10.1063/5.0046938
[8] Chen SY, Doolen GD. Lattice Boltmzann emthod for fluid flows. Annual Review of Fluid Mechanics, 1998, 30(1): 329-364 doi: 10.1146/annurev.fluid.30.1.329
[9] Dudek M, Øye G. Microfluidic study on the attachment of crude oil droplets to gas bubbles. Energy & Fuels, 2018, 32(10): 10513-10521
[10] Liu HH, Kang QJ, Leonardi CR, et. al. Multiphase lattice Boltzmann simulations for porous media applications. Computational Geosciences, 2016, 20: 777-805
[11] Wang NN, Ni WL, Liu HH. Geometrical wetting boundary condition for complex geometries in lattice Boltzmann color-gradient model. Physics of Fluids, 2024, 36(1): 012109 doi: 10.1063/5.0180592
[12] Qin ZR, Zhu JF, Chen WB, et al. An effective pseudo-potential lattice Boltzmann model with extremely large density ratio and adjustable surface tension. Physics of Fluids, 2022, 34(11): 113328 doi: 10.1063/5.0123727
[13] Zhang CH, Guo ZL, Wang LP. Improved well-balanced free-energy lattice Boltzmann model for two-phase flow with high Reynolds number and large viscosity ratio. Physics of Fluids, 2022, 34: 012110 doi: 10.1063/5.0072221
[14] Zu YQ, He SS. Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts. Physical Review E, 2013, 87(4): 043301 doi: 10.1103/PhysRevE.87.043301
[15] Huang JJ, Zhang LQ. Simplified method for wetting on curved boundaries in conservative phase-field lattice-Boltzmann simulation of two-phase flows with large density ratios. Physics of Fluids, 2022, 34(8): 82101 doi: 10.1063/5.0101291
[16] Schwarzmeier C, Holzer M, Mitchell T, et al. Comparison of free-surface and conservative Allen-Cahn phase-field lattice Boltzmann method. Journal of Computational Physics, 2023, 475(15): 111753
[17] Lee T, Lin CL. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. Journal of Computational Physics, 2005, 206(1): 16-47 doi: 10.1016/j.jcp.2004.12.001
[18] Zheng HW, Shu C, Chew YT. Lattice Boltzmann interface capturing method for incompressible flows. Physical Review E, 2005, 72(5): 056705 doi: 10.1103/PhysRevE.72.056705
[19] Liang H, Shi BC, Guo ZL, et al. Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows. Physical Review E, 2014, 89(5): 053320 doi: 10.1103/PhysRevE.89.053320
[20] Wang K, Xia YC, Li ZY. A phase-field Lattice Boltzmann method for liquid-vapor phase change problems based on conservative Allen-Cahn equation and adaptive treegrid. Computers & Fluids, 2023, 264: 105973
[21] Hu Y, Li DC, He Q. Generalized conservative phase field model and its lattice Boltzmann scheme for multicomponent multiphase flows. International Journal of Multiphase Flow, 2020, 132: 103432 doi: 10.1016/j.ijmultiphaseflow.2020.103432
[22] Chen XW, Null XQ, Song SH. Fourth-order structure-preserving method for the conservative Allen-Cahn equation. Advances in Applied Mathematics and Mechanics, 2023, 15(1): 159-181 doi: 10.4208/aamm.OA-2021-0325
[23] 章诗婷, 肖鸿威, 周锦翔等. 广义守恒相场简化多相流格子Boltzmann方法. 空气动力学学, 2022, 40(3): 75-86 (Zhang Shiting, Xiao Hongwei, Zhou Jinxiang. Generalized conservative phase-field model based lattice Boltzmann method for multiphase flows. Acta Aerodynamica Sinica, 2022, 40(3): 75-86 (in Chinese) Zhang Shiting, Xiao Hongwei, Zhou Jinxiang. Generalized conservative phase-field model based lattice Boltzmann method for multiphase flows. Acta Aerodynamica Sinica, 2022, 40(3): 75-86 (in Chinese)
[24] Geier M, Fakhari A, Lee T. Conservative phase-field lattice Boltzmann model for interface tracking equation. Physical Review E, 2015, 91(6): 063309 doi: 10.1103/PhysRevE.91.063309
[25] Fakhari A, Bolster D, Luo LS. A weighted multiple-relaxation-time lattice Boltzmann method for multiphase flows and its application to partial coalescence cascades. Journal of Computational Physics, 2017, 341: 22-43 doi: 10.1016/j.jcp.2017.03.062
[26] Fakhari A, Geier M, Lee T. A mass-conserving lattice Boltzmann method with dynamic grid refinement for immiscible two-phase flows. Journal of Computational Physics, 2016, 315: 434-457 doi: 10.1016/j.jcp.2016.03.058
[27] Liang H, Wang RL, Wei YK, et al. Lattice Boltzmann method for interface capturing. Physical Review E, 2023, 107(2): 025302 doi: 10.1103/PhysRevE.107.025302
[28] Ren F, Song BW, Sukop MC, et al. Improved lattice Boltzmann modeling of binary flow based on the conservative Allen-Cahn equation. Physical Review E, 2016, 94(2): 023311 doi: 10.1103/PhysRevE.94.023311
[29] Wang HL, Chai ZH, Shi BC, et al. Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations. Physical Review E, 2016, 94(3): 033304 doi: 10.1103/PhysRevE.94.033304
[30] Liang H, Xu JR, Chen JX, et al. Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows. Physical Review E, 2018, 97(3): 033309 doi: 10.1103/PhysRevE.97.033309
[31] Xu XC, Hu YW, Dai Bing, et al. Modified phase-field-based lattice Boltzmann model for incompressible multiphase flows. Physical Review E, 2021, 104(3): 035305 doi: 10.1103/PhysRevE.104.035305
[32] Sterling JD, Chen SY. Stability analysis of lattice Boltzmann methods. Journal of Computational Physics, 1996, 123(1): 196-206 doi: 10.1006/jcph.1996.0016
[33] Latt J, Chopard B. Lattice Boltzmann method with regularized pre-collision distribution functions. Mathematics and Computers in Simulation, 2006, 72(2): 165-173
[34] Kawaguchi M, Fukui T, Morinishi K. Comparative study of the virtual flux method and immersed boundary method coupled with regularized lattice Boltzmann method for suspension flow simulations. Computers & Fluids, 2022, 246: 105615
[35] Wang L, Yang XG, Wang HL, et al. A modified regularized lattice Boltzmann model for convection-diffusion equation with a source term. Applied Mathematics Letters, 2021, 112: 106766 doi: 10.1016/j.aml.2020.106766
[36] 李旭晖, 单肖文, 段文洋. 格子玻尔兹曼正则化碰撞模型的理论进展. 空气动力学学报, 2022, 40(3): 49-64 (Li Xuhui, Shan Xiaowen, Duan Wenyang. Theoretical progress on regularized lattice Boltzmann collision models. Acta Aerodynamica Sinica, 2022, 40(3): 49-64 (in Chinese) Li Xuhui, Shan Xiaowen, Duan Wenyang. Theoretical progress on regularized lattice Boltzmann collision models. Acta Aerodynamica Sinica, 2022, 40(3): 49-64 (in Chinese)
[37] Feng YL, Guo SL, Tao WQ, et al. Regularized thermal lattice Boltzmann method for natural convection with large temperature differences. International Journal of Heat & Mass Transfer, 2018, 125: 1379-1391
[38] Wang L, Chai ZH, Shi BC. Regularized lattice Boltzmann simulation of double-diffusive convection of power-law nanofluids in rectangular enclosures. International Journal of Heat and Mass Transfer, 2016, 102: 381-395 doi: 10.1016/j.ijheatmasstransfer.2016.06.041
[39] Wang L, Wei ZC, Li TF, et al. A lattice Boltzmann modelling of electrohydrodynamic conduction phenomenon in dielectric liquids. Applied Mathematical Modelling, 2021, 95: 361-378 doi: 10.1016/j.apm.2021.01.054
[40] Ba Y, Wang NN, Liu HH, et al. Regularized lattice Boltzmann model for immiscible two-phase flows with power-law rheology. Physical Review E, 2018, 97(3): 033307 doi: 10.1103/PhysRevE.97.033307
[41] 杨旭光, 汪垒. 多相多组分Peng-Robinson流体的正则化格子Boltzmann方法研究及模拟. 力学学报, 2023, 55(8): 1649-1661 (Yang Xuguang, Wang Lei. Regularized lattice Boltzmann method for multi-component and multi-phase Peng-Robinson fluids. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(8): 1649-1661 (in Chinese) Yang Xuguang, Wang Lei. Regularized lattice Boltzmann method for multi-component and multi-phase Peng-Robinson fluids. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(8): 1649-1661 (in Chinese)
[42] Otomo H, Crouse B, Dressler M, et al. Multi-component lattice Boltzmann models for accurate simulation of flows with wide viscosity variation. Computers & Fluids, 2018, 172: 674-682
[43] Jacqmin D. Calculation of two-phase Navier–Stokes flows using phase-field modeling. Journal of Computational Physics, 1999, 155(1): 96-127 doi: 10.1006/jcph.1999.6332
[44] Qian YH, Humières DD, Lallemand P. Lattice BGK models for Navier-Stokes equation. Europhysics Letters, 1992, 17(6): 479-496 doi: 10.1209/0295-5075/17/6/001
[45] Fakhari A, Mitchell T, Leonardi C, et al. Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios. Physical Review E, 2017, 96(5): 053301 doi: 10.1103/PhysRevE.96.053301
[46] Li X, Dong ZQ, Li Y, et al. A fractional-step lattice Boltzmann method for multiphase flows with complex interfacial behavior and large density contrast. International Journal of Multiphase Flow, 2022, 149: 103982 doi: 10.1016/j.ijmultiphaseflow.2022.103982
[47] Li QZ, Lu ZL, Chen Z, et al. An efficient simplified phase-field lattice Boltzmann method for super-large-density-ratio multiphase flow. International Journal of Multiphase Flow, 2023, 160: 104368 doi: 10.1016/j.ijmultiphaseflow.2022.104368
[48] Yuan HZ, Chen Z, Shu C, et al. A free energy-based surface tension force model for simulation of multiphase flows by level-set method. Journal of Computational Physics, 2017, 345: 404-426 doi: 10.1016/j.jcp.2017.05.020
[49] Hysing S, Turek S, Kuzmin D, et al. Quantitative benchmark computations of two-dimensional bubble dynamics. International Journal for Numerical Methods in Fluids, 2009, 60(11): 1259-1288 doi: 10.1002/fld.1934
[50] Chen Z, Shu C, Tan DS, et al. Simplified multiphase lattice Boltzmann method for simulating multiphase flows with large density ratios and complex interfaces. Physical Review E, 2018, 98(6): 063314 doi: 10.1103/PhysRevE.98.063314
[51] Josserand C, Zaleski S. Droplet splashing on a thin liquid film. Physics of Fluids, 2003, 15(6): 1650-1657 doi: 10.1063/1.1572815
下载:
计量
- 文章访问数: 207
- HTML全文浏览量: 173
- PDF下载量: 76
