摘要:本文采用相场格子Boltzmann方法研究了竖直微通道内中等Atwoods数流体的单模Rayleigh-Taylor不稳定性问题, 系统分析了雷诺数对相界面动力学行为以及扰动在各发展阶段演化规律的影响. 数值结果表明高雷诺数条件下, 不稳定性界面扰动的增长经历了四个不同的发展阶段, 包括线性增长阶段、饱和速度阶段、重加速阶段及混沌混合阶段. 在线性增长阶段, 我们计算获得的气泡与尖钉振幅符合线性稳定性理论, 并且线性增长率随着雷诺数的增加而增大. 在第二个阶段, 我们观察到气泡与尖钉将以恒定的速度增长, 获得的尖钉饱和速度略高于Goncharov经典势能模型的解析解[Phys. Rev. Lett. 2002 88 134502], 这归因于系统中产生了多个尺度的旋涡, 而涡之间的相互作用促进了尖钉的增长. 随着横向速度和纵向速度的差异扩大, 气泡和尖钉界面演化诱导产生的Kelvin–Helmholtz不稳定性逐渐增强, 从而流体混合区域出现许多不同层次的涡结构, 加速了气泡与尖钉振幅的演化速度, 并在演化后期阶段, 导致界面发生多层次卷起、剧烈变形、混沌破裂等行为, 最终形成了非常复杂的拓扑结构. 此外, 我们还统计了演化后期气泡与尖钉的无量纲加速度, 发现气泡和尖钉的振幅在后期呈现二次增长规律, 其增长率系数分别为0.045与0.233. 而在低雷诺条件下, 重流体在不稳定性后期以尖钉的形式向下运动而轻流体以气泡的形式向上升起. 在整个演化过程中, 界面变得足够光滑, 气泡与尖钉在后期的演化速度接近于常数, 未观察到后期的重加速与混沌混合阶段.
关键词: 格子Boltzmann方法/
Rayleigh-Taylor不稳定性/
相场方法/
雷诺数/
二次增长

English Abstract

全文HTML

--> --> -->

瑞利-泰勒(Rayleigh-Taylor, RT)不稳定性是一种经典而又古老的流体界面不稳定性现象. 当重流体置于轻流体之上并在相界面处施加一个微小的扰动, 两相流体的界面是不稳定的, RT不稳定性现象将会发生. RT不稳定性问题广泛存在于天体物理(卷云、超星系爆炸、蟹状星云等)^[1]、地球物理(盐丘、油气矿的形成)^[2]、工程界(混合过程、惯性约束聚变)^[3], 也在核能约束聚变热能利用中扮演着重要角色. 另外, 我国著名核物理学家贺贤土院士指出有效抑制RT不稳定性的后期湍流混合是惯性约束聚变过程中实现点火成功的关键^[4]. 因此, RT不稳定性问题在流体力学领域具有重要的理论价值和广泛的工程应用, 从而长期以来吸引许多****对其开展了相关研究. 早在1883年, 著名****Rayleigh^[5]在研究云的分层问题中第一次描述了RT不稳定性现象. 后来, Taylor^[6]在原子能武器的研究中进一步发现了RT不稳定性现象. 他们对RT不稳定性进行理论分析并提出了描述扰动发展的经典线性增长理论: 当界面初始扰动的振幅小于它的波长时, 扰动将变得不稳定, 其振幅将以指数形式增长. 后来, 一些****相继分析了流体压缩性、流体黏性、界面张力对扰动增长率的影响, 并提出了修正的线性增长理论^[7–9]. Lewis^[10]第一次通过实验研究了RT不稳定性问题, 验证了Taylor的线性增长理论, 并初步而定性地描述了不稳定性中扰动的发展阶段. 后来, Sharp^[11]进一步研究了RT不稳定性问题, 并较为完善地划分了不稳定性发展的四个阶段: 初始阶段, 扰动振幅符合经典的线性增长理论, 将以指数形式增长, 直至增长到初始扰动波长的$ 10- 40\% $; 紧接着, 界面扰动增长转变为非线性发展阶段, 表现为重流体渗透到轻流体中呈现出蘑菇的形状, 称为尖钉, 而轻流体上浮至重流体中, 称为气泡; 在轻重流体相互渗透的过程中, 两相系统的非线性强度逐渐增强, 蘑菇头部发生挤压而尾部出现卷吸现象, 这表明产生了第二类不稳定性现象, 即开尔文-赫姆霍兹(Kelvin-Helmholtz, KH)不稳定性; 最后, 受两类不稳定性的双重影响, 气泡和尖钉发生破裂行为, 流体界面表现为后期的混沌状态.
自Sharp^[11]的研究工作以来, 许多****对RT不稳定性的各个发展阶段开展了研究, 包括单模态RT不稳定性和随机扰动模态的RT不稳定性^[12–15]. 本文主要关注相界面初始扰动为单模态的RT不稳定性问题. Waddell等^[16]实验研究了低Atwood数下二维单模RT不稳定性, 其结果表明扰动振幅在初始阶段以指数的形式增长, 获得的线性增长率与经典的线性理论分析相一致, 还进一步发现尖钉和气泡的平均速度在演化后期接近于常数. Goncharov^[17]通过理论分析证实了这一点, 并定量地给出理想流体的尖钉和气泡恒定的饱和速度

$ u_{\rm b} = \sqrt{\frac{2A_tg}{C_gk(1+A_{\rm t})}},\; \; \; u_{\rm s} = \sqrt{\frac{2A_tg}{C_gk(1-A_{\rm t})}}, $

其中$ u_{\rm b} $是气泡的速度, $ u_{\rm s} $是尖钉的速度, g是重力加速度, $ A_{\rm t} $是Atwood数, k是波数, $ C_g $对于三维情形取1, 而对于二维情形则取为3. Wilkinson和Jacobs^[18]对三维单模RT不稳定性问题进行了实验研究, 观察到气泡与尖钉在演化后期的增长速度均高于Goncharov^[17]的势能理论模型的解析解, 并将这种与理论预测不一致现象归因于涡结构间的相互作用. 后来, 一些****还分析了涡结构相互作用、流体黏性、表面张力对第二阶段中气泡演化速度的影响, 并提出了包含这些物理因素的理论表达式.
单模态RT不稳定性的上述研究还仅限于扰动发展的前三个阶段, 而对于不稳定性的后期阶段, 非线性效应往往越来越强烈, 流体界面则会显示出非常复杂而剧烈的拓扑变化, 因而难以通过实验方法测量演化后期的相关物理统计量, 也难以通过理论分析方法来求解. 近几十年来, 随着计算机计算技术与数值方法的快速发展, 数值模拟已经成为一种基本的研究手段, 在RT不稳定性研究方面发挥了越来越重要的作用. He等^[19]采用一个双分布等温格子Boltzmann (Lattice Boltzmann, LB)方法模拟了三维单模RT不稳定性问题, 分析了雷诺数和Atoowd数对不稳定性界面结构的影响. Glimm等^[20]基于欧拉方程的前追踪方法研究了二维单模RT不稳定性问题, 发现不稳定性的发展经历了饱和恒定速度阶段后会出现一个加速阶段, 这个阶段后来被Wilkinson和Jacobs^[18]实验所发现. Celani^[21]利用相场方法模拟了二维单模RT不稳定性问题, 并在初始阶段获得了与线性理论相符合的扰动增长率. 需要指出的是, 大多数的数值研究是用于验证所发展的计算流体力学方法的正确性与有效性的, 研究阶段也局限于不稳定性发展的中前期. Ramaprabhu等^[22]模拟了三维混相流体的单模RT不稳定性现象, 考察了Atwood数和雷诺数对气泡和尖钉振幅在后期阶段的演化规律. 他们发现低Atwood数的不稳定性在高雷诺数情形下经历了指数增长、饱和速度增长、重加速、混沌混合等四个发展阶段, 并进一步指出在重加速阶段, 气泡和尖钉速度已经超过经典势能模型所预测的理论解, 而在后期阶段, 气泡和尖钉速度会出现急剧的下降. Wei等^[23]采用直接数值模拟方法研究了低Atwoods数下二维混相流体的单模RT不稳定性问题, 分析了雷诺数对气泡和尖钉振幅增长的影响. 他们同样地观察到混相不稳定性在高雷诺数时经历了一系列的发展阶段, 但不同于Ramaprabhu等^[22]的结果, 气泡演化速度在后期的混沌阶段出现了随时间波动的现象, 并且其平均加速度接近于常数, 这表明气泡在演化后期具有二次增长的规律. Liang等^[24,25]利用基于相场理论的多相流LB方法模拟二维及三维长微管道内低Atwoods数下非混相流体的单模RT不稳定性问题, 详细地考察了雷诺数对不稳定性各发展阶段扰动增长和流体界面动态行为的影响. 数值结果表明高雷诺数时, 不稳定性也经历了四个发展阶段, 包括线性增长阶段、饱和速度增长阶段、重加速阶段及混沌混合阶段. 在后期阶段, 相界面发生多层次卷起、剧烈变形、混沌破裂等行为, 形成了非常复杂的拓扑结构, 并且气泡振幅显示二次增长的规律. 而低雷诺数时, 相界面的演化变得相对光滑, 后期的发展阶段也相继难以到达. 最近, Hu等^[26]数值研究了低Atwoods数混相流体的单模RT不稳定性问题, 发现气泡速度在重加速阶段后会出现不停地加速与减速, 并将这种现象归因于旋涡强度的变化.
综上所述, 已有许多****^[16–25]对单模RT不稳定性问题开展了研究, 丰富了人们对不稳定性发展规律的认识. 然而, 绝大多数的工作^[16–21]局限于不稳定性发展的前期阶段, 并且所考虑的雷诺数均较小. 尽管一些****^[22–26]对高雷诺数下不稳定性演化的后期阶段开展了部分研究, 但考虑的两种流体为相互混溶的^[22,23,26], 所关注的流体间Atwood数也太小^[23–26]. 鉴于此, 本文将采用基于相场理论的多松弛LB方法^[24]研究中等Atwood数下非混相流体RT不稳定性的演化规律, 着重分析雷诺数对相界面动力学行为和扰动后期增长的影响.

基于分子动理学理论的LB方法^[27]是近十几年发展起来的介观数值方法, 它不再基于宏观连续介质模型, 而通过描述流体粒子分布函数的演化再现复杂流动的宏观行为, 从而相比传统数值方法具有一些独特的优势, 比如易于处理复杂物理边界、程序实现简单且天然并行、直观刻画多相流体间微观相互作用. 目前, 从流体组分间相互作用力的不同物理背景出发, 已经提出了多种描述多相流体输运的LB模型, 包括颜色模型、伪势模型、自由能模型、基于相场理论的LB方法. 而基于相场理论的LB方法在描述多相流体界面动力学方面具有清晰的物理机制, 可用于模拟多相系统中界面具有拓扑变化的流动, 因而受到了许多****的广泛关注, 并已成功地应用于轴对称多相流^[28,29]、三相流^[30,31]、复杂微通道内液滴动力学^[32–34]、液滴润湿固体表面^[35]等复杂多相流问题. 本文将采用Liang等^[24]提出的相场LB模型研究长微通道中非混相RT不稳定性的后期演化规律, 该模型相比前人相场LB模型可正确地恢复Cahn-Hilliard和Navier-Stokes的耦合方程, 并且速度和压力场可通过分布函数进行显示计算. 另外, 为了提高模型在高雷诺数时的数值稳定性, 我们还在相场LB模型的演化过程中采用多松弛的碰撞算子^[36]. 本文相场LB模型利用两个独立的分布函数$ f_i $和$ g_i $, 其对应的多松弛演化方程可表述为:

$\begin{split} &{f_i}({{x}} + {{{e}}_i}\delta_x,t + \delta_t) - {f_i}({{x}},t) \\ =\, & - ({{{M}}^{ - 1}}{{{S}}^f}{{M}})_{ij}\left[ {{f_j}({{x}},t) \!-\! f_j^{\rm {eq}}({{x}},t)} \right] \\&\!+\! {\delta_t}{F_i}({{x}},t), \end{split}$

$\begin{split} & {g_i}({{x}} + {{{e}}_i}\delta_x,t + \delta_t) - {g_i}({{x}},t) \\ =\, & - ({{{M}}^{ - 1}}{{{S}}^g}{{M}})_{ij}\left[ {{g_j}({{x}},t) - g_j^{\rm {eq}}({{x}},t)} \right]\\& + {\delta_t}{G_i}({{x}},t), \end{split}$

其中$ {f_i}({{x}}, t) $是描述粒子在空间$ { x} $和时间t时的序参数分布函数, $ {g_i}({{x}}, t) $是密度分布函数, $ f_i^{\rm {eq}}({{x}}, t) $,$ g_i^{\rm {eq}}({{x}}, t) $为两分布函数所对应的平衡态分布函数, $ {M} $是分布函数空间到矩空间的变换矩阵, $ {{{S}}^f} $, $ {{{S}}^g} $为松弛矩阵, $ {F_i}({{x}}, t) $、$ {G_i}({{x}}, t) $分别为源项和外力项的分布函数, $ {\delta_x} $和$ {\delta_t} $为空间和时间步长. 为了匹配宏观控制方程, 平衡态分布函数$ f_i^{\rm {eq}}({\mathit{\boldsymbol{ x}}}, t) $和$ g_i^{\rm {eq}}({\mathit{\boldsymbol{ x}}}, t) $可设定为^[24]

$f_i^{\rm {eq}} = \left\{ {\begin{aligned}& {\phi + ({\omega _i} - 1)\eta \mu ,}\qquad {i = 0},\\& {{\omega _i}\eta \mu + {\omega _i}\frac{{{{{e}}_i} \cdot \phi {{u}}}}{{c_{\rm s}^2}},}\quad {i \ne 0},\end{aligned}} \right.$

$ g_i^{eq} = \left\{ {\begin{aligned}&{\frac{p}{{c_{\rm s}^2}}({\omega _i} - 1) + \rho {s_i}({{u}})},\quad {i = 0},\\& {\frac{p}{{c_{\rm s}^2}}{\omega _i} + \rho {s_i}({{u}})},\qquad~~~ {i \ne 0},\end{aligned}} \right.$

且

$ {s_i}({\bf{u}}) = {\omega_i}\left[ {{{{{{e}}_i} \cdot {{u}}} \over {c_{\rm s}^2}} + {{{{({{{e}}_i} \cdot {{u}})}^2}} \over {2c_{\rm s}^4}} - {{{{u}} \cdot {{u}}} \over {2c_{\rm s}^2}}} \right], $

其中$ \phi $是序参数, $ \mu $是化学势, $ \eta $为调节迁移率的自由参数, $ \rho $是密度, p是流体动力学压力, $ {u} $是流体速度, $ {\omega_i} $为权系数, $ c_{\rm s} $为声速. 根据相场理论, 范德瓦尔斯(van der Waals)流体的化学势$ \mu $可表述为序参数的相关函数, 具体定义为

$ \mu = 4\beta \phi(\phi - 1)(\phi - 0.5)-k{\nabla ^2}\phi, $

其中参数k,$ \beta $与界面厚度$ (D) $、表面张力$ (\sigma) $相关,

$ k = 1.5D\sigma,\; \; \; \beta = \frac{12\sigma}{D}. $

另外, $ {\omega_i} $和$ c_{\rm s}^2 $的选取依赖于格子速度的离散模型. 针对本文所研究的二维微通道流动, 我们采用最有代表性的D2 Q9格子离散速度模型^[37,38], 其权系数和声速为$ {\omega_0} = 4/9 $, $ {\omega_{1-4}} = 1/9 $, $ {\omega_{5-8}} = 1/36 $, $ c_{\rm s}^2 = c^2/3 $, 离散速度$ \mathit{\boldsymbol{e}}_{i} $为

$ {{{e}}_i} \!=\! \left\{ {\begin{aligned}& {(0,0),}\qquad\qquad \qquad \qquad \qquad\qquad {i = 0},\\& {(\cos [(i \!-\! 1){\text{π}} /2],\sin [(i \!- \!1){\text{π}} /2]),}~ {i = 1—4},\\& {\sqrt 2 (\cos [(i \!-\! 5){\text{π}} /2 \!+\! {\text{π}} /4],\sin [(i \!-\! 5){\text{π}} /2 \!+\! {\text{π}} /4]),}\\ &~~\qquad \qquad \qquad\qquad \qquad \qquad \qquad {i = 5—8},\end{aligned}} \right. $

且采用规则的正方形格子, 其对应的变换矩阵为^[36]

$\begin{array}{*{20}{c}}{{M}} \!=\! \left(\!\!\!\!{\begin{array}{*{20}{c}}1&1&1&1&1&1&1&1&1\\{ - 4}&{ - 1}&{ - 1}&{ - 1}&{ - 1}&2&2&2&2\\4&{ - 2}&{ - 2}&{ - 2}&{ - 2}&1&1&1&1\\0&1&0&{ - 1}&0&1&{ - 1}&{ - 1}&1\\0&{ - 2}&0&2&0&1&{ - 1}&{ - 1}&1\\0&0&1&0&{ - 1}&1&1&{ - 1}&{ - 1}\\0&0&{ - 2}&0&2&1&1&{ - 1}&{ - 1}\\0&1&{ - 1}&1&{ - 1}&0&0&0&0\\0&0&0&0&0&1&{ - 1}&1&{ - 1}.\end{array}}\!\!\!\right).\end{array}$

同样地, 在多松弛的D2Q9格子模型, $ f_i $, $ g_i $所对应的松弛矩阵分别定义为:

$ {{S}}^f = (1,\; s_1^f,\; s_1^f,\; s_2^f,\; s_3^f,\; s_2^f,\; s_3^f,\; 1,\; 1), \tag{10a}$

$ {{S}}^g = (1,\; 1,\; 1,\; 1,\; s_1^g,\; 1,\; s_1^g,\; s_2^g,\; s_2^g), \tag{10b}$

其中$ s_2^f $是与相场方程中迁移率有关的松弛因子, 松弛因子$ s_2^g $的值取决于流体黏性, $ s_1^f $, $ s_3^f $, $ s_1^g $为可调节的松弛因子. Lallemand和Luo^[36]通过理论分析表明, 调节这些自由的松弛参数可以有效地提高LB方法的数值稳定性. 另外, Liang等^[24,25]通过数值实验发现多松弛模型适用于研究高雷诺数的RT界面不稳定性问题. 进一步, 为了恢复正确的相界面追踪Cahn-Hilliard方程, LB演化方程中的源项分布函数$ {F_i}({{x}}, t) $和外力项分布函数$ {G_i}({{x}}, t) $分别定义为:^[24,39]

$ {{F}_i} = \left[{ M}^{-1}\left({ I}-\frac{S^f}{2}\right){ M}\right]_{ij}{{{\omega_j}{{{e}}_j} \cdot {\partial_t}\phi {{u}}} \over {c_{\rm s}^2}}, $

$\begin{split} {{G}_i} =\; & \left[{ M}^{-1}\left({ I}-\frac{S^g}{2}\right){ M}\right]_{ij}\\ & \times{\omega_j}\left[\frac{{ e}_j\cdot ({ F}_{\rm s}+{ G})}{c_{\rm s}^2}+\frac{{{u}}{\nabla \rho}:{{{e}}_j}{{{e}}_j}}{c_{\rm s}^2}\right],\end{split} $

其中$ { I} $是单位矩阵, $ {F}_{\rm s} $是界面张力, $ { G} $是外力项. 在相场模型中, 表面张力取为势能形式$ {F}_{\rm s} = \mu\nabla \phi $, 发现能够有效地减少相界面处的虚假速度.
在本模型中, 流体的宏观量可以通过求解粒子分布函数的各阶矩获得, 具体的计算式如下^[39]:

$ \phi = \sum\limits_i {{f_i}}, $

$ \rho{{u}} = {{\sum\limits_i {{{{e}}_i}{g_i}} + \frac{{\delta _t}}{2}({{{F}}_{\rm s}+{{G}}})}}, $

$ p = {{c_{\rm s}^2} \over {(1 - {\omega _0})}}\left[ {\sum\limits_{i \ne 0} {{g_i}} + {{{\delta _t}} \over 2}{{u}} \cdot \nabla \rho + \rho {s_0}({u})} \right], $

而宏观量密度($ \rho $)可以看成是关于序参数($ \phi $)的一个简单的线性插值函数,

$ \rho = \phi(\rho_{\rm h}-\rho_{\rm l})+\rho_{\rm l}, $

其中$ \rho_{\rm h} $和$ \rho_{\rm l} $分别为重质流体和轻质流体的密度.
最后, 通过Chapman-Enskog多尺度理论分析^[24,39], 可以证明本文采用的基于相场理论的多松弛LB模型可以正确地恢复界面追踪的Cahn-Hilliard控制方程,

$ {{\partial \phi } \over {\partial t}} + \nabla \cdot \phi {{u}} = \nabla \cdot {M}(\nabla {\mu} ), $

和描述流体动力学的不可压Navier-Stokes方程组

$ \nabla \cdot {{u}} = 0, $

$ {{\partial { \rho {u}}} \over {\partial t}} + \nabla \cdot({\rho {{u}}{u}}) = - \nabla p + \nabla \cdot \left[ {\nu\rho(\nabla {{u}} + \nabla {{{u}}^{\rm T}})} \right] + {{{F}}_{\rm s}} + {{G}}, $

其中 迁移率M、流体运动性黏性$ \nu $与松弛 因子关系可分别表示为:

$ M =\eta c_{\rm s}^2 \left(\frac{1}{s_2^f}-0.5\right)\delta_t, $

$ \nu = c_{\rm s}^2 \left(\frac{1}{s_2^g}-0.5\right)\delta_t. $

本文将采用基于相场理论的多松弛LB模型研究二维长管道内非混相流体的单模RT界面不稳定性问题, 并详细地考察雷诺数对界面不稳定性发展的影响, 以及定量分析气泡与尖钉的振幅随着雷诺数的变化规律. 为了研究不稳定性的后期演化规律, 本文考虑的物理问题为一个长方形的微通道, 其高度和宽度分别为H和W, 且$ H/W = 8 $. 初始时刻, 密度较大的重流体置于另一种轻质流体的上方, 且给流体界面处一个微小的扰动, 在重力场的作用下, 扰动会逐渐发展, 并最终达到混沌混合状态. 在长方形微通道的中心截面处, 施加一个波长为W具有余弦函数的微小初始扰动,

$ h(x,y) = 0.5H+0.05W\cos(kx), $

其中$ k = 2{\text{π}}/W $是扰动波数. 为了使序参数变量在相界面处光滑连续的变化, 设定序参数$ \phi(x, y) $的初始分布为

$ \phi(x,y) = 0.5+0.5\tanh \frac{2(y-h)}{D}. $

为了表征二维RT不稳定性的演化特征, 引了两个常用的无量纲参数, 即雷诺数($ Re $)和阿特伍德数($ A_t $), 分别定义为^[12,13]

$ Re = \frac{W\sqrt{gW}}{\nu},\; \; \; \; \; \; A_t = \frac{\rho_{\rm h}-\rho_{\rm l}}{\rho_{\rm h}+\rho_{\rm l}}, $

其中g是重力加速度, $ {\nu} $是流体运动性黏性. 在本文的研究中, 我们考虑中等的阿特伍德数下界面不稳定性的演化规律, 将轻重流体的密度分别设定为1和3, 其对应的Atwoods数为0.5, 其他的物理参数选取如下: $ W = 256 $, $ D = 4 $, $ \sigma = 5\times 10^{-5} $, $ \sqrt{gW} = 0.04 $. 为了演化计算, 需要选择合适的边界格式处理物理边界处的粒子分布函数, 本文对上下固壁采用无滑移的半反弹边界条件, 左右边界应用周期边界条件. 另外, 标记气泡和尖钉的振幅为$ H_{\rm b} $和$ H_{\rm s} $, 分别定义为气泡与尖钉的前端到对应初始位置的距离, 因此气泡和尖钉的振幅在初始时刻的值为0. 进一步, 我们将管道宽度W和$ 1/\sqrt{A_tgk} $选为特征长度和特征时间, 还定义了气泡和尖钉无量纲的演化速度, 也被称为气泡和尖钉的Froude数^[23],

$ Fr_{\rm b} = \frac{u_{\rm b}}{\sqrt{\dfrac{A_tgW}{1+A_t}}},\; \; \; \; \; Fr_{\rm s} = \frac{u_{\rm s}}{\sqrt{\dfrac{A_tgW}{1+A_t}}}, $

其中$ u_{\rm b} $和$ u_{\rm s} $为气泡和尖钉前端的速度, 可以由气泡和尖钉振幅计算获得. 除了特别声明, 本文接下来给出的长度、速度、时间($ \tau $)等相关物理量均已被相应的特征值所无量纲化.
图1描述了四种典型的不同雷诺数下非混相RT不稳定性中相界面的演化过程. 从图1中可以发现, 对于不同的雷诺数, 不稳定性的扰动在初始阶段显示相似的界面动力学特行为: 重流体在重力作用往下运动而轻流体向上浮起, 即轻重流体之间相互渗透, 从而形成了气泡和尖钉图案. 紧接着, 流体界面在不同的雷诺数下展示出显著不同的动力学特征. 对高Re情形(Re = 10000), 尖钉继续向下运动并逐渐地向上卷起, 形成了旋转方向相反的两个旋涡, 这是KH不稳定性出现并作用于相界面的结果. 随着时间的演化, 两个旋涡不断地增长, 在卷起的尾端处形成了一对二级旋涡. 随后, 在多个旋涡相互作用下, 不稳定性系统的非线性效应越来越剧烈, 尖钉卷起的长度也越来越长, 并逐渐地靠近进而接触中轴线附近的流体界面. 在高流体界面剪切力作用下, 中轴线附近的界面在多处位置出现卷起与变形行为. 最终, 流体界面发生了混沌的破裂, 在系统中形成了许多离散的小液滴. 另外, 我们还观察到流体界面在整个演化过程中始终保持关于中轴线对称. 当Re减少至2048时, 同样观察到尖钉发生卷起行为, 在尾端也形成了一对二级涡, 并最终导致相界面在多个位置发生卷吸、变形和破裂, 形成较为复杂的结构. 然而, 相比Re = 10000的情形, 在演化后期, 系统中相界面的混沌程度减弱了. 当Re数进一步降低至50时, 重流体的尖钉往下运动, 经过一段演化时间, 在尾端发生卷吸现象, 也形成了一对旋转方向相互相反的旋涡, 但与高雷诺数情形相比, 界面卷起发生的时刻推迟了, 旋涡的卷吸幅度也相应地减弱了. 最后, 形成的旋涡随着时间演化而不断地发展, 伴随着尾端卷起的部分也越来越长. 在整个演化过程中, 未观察到二次旋涡卷吸和界面后期发生破裂的现象. 当Re充分小(Re = 5)时, 卷吸现象不再发生, 重流体将以尖钉的方式不断地向下运动, 界面也变得足够光滑, 未出现混沌破裂等复杂拓扑现象, 这是由于强黏性作用使流体之间的剪切层保持稳定, 流动在整个过程表现为层流状态.
图 1 雷诺数对非混相RT不稳定性中相界面演化图案的影响　(a) $ Re = 10000 $; (b) $ Re = 2048 $; (c) $ Re = 50 $; (d) $ Re = 5 $
Figure1. The effect of the Reynolds number on the evolution of interfacial patterns in the immiscible RT instability: (a) $ Re = 10000 $; (b) $ Re = 2048 $; (c) $ Re = 50 $; (d) $ Re = 5 $.

上面讨论了雷诺数对单模RT不稳定性中相界面动力学行为的影响, 而气泡与尖钉振幅及演化速度是描述RT不稳定性问题中另外两个非常重要的物理量. 为了进一步显示雷诺数的效应, 我们定量地分析了不同雷诺数下气泡和尖钉振幅、运动速度随时间的演化规律. 图2分别给出了气泡与尖钉在不同雷诺数下随时间变化的演化曲线. 从图中可以发现, 对所有的Re数情形, 气泡和尖钉的振幅均随着时间演化而不断增大. 当Re逐渐增大时, 可以观察到同一时刻所获得的尖钉振幅也越大, 而当Re增大至足够大时, 雷诺数对尖钉振幅的影响将不再显著. 在不可压流体的RT不稳定性中, 尖钉的运动特征在理论上由单位质量的浮力和黏性耗散力之间的竞争关系决定^[40],
图 2 雷诺数对无量纲化的气泡与尖钉随时间演化振幅的影响
Figure2. The effect of the Reynolds number on the dimensionless bubble and spike amplitudes.

$ \frac{{\rm d} H_{\rm s}}{{\rm d}t} = u_{\rm s},\; \; \; \; \; \frac{{\rm d} u_{\rm s}}{{\rm d}t} = \frac{\text{δ} \rho}{\bar{\rho}}g+ {F_{\rm d}}, $

其中 $\text{δ} \rho = (\rho_{\rm h}-\rho_{\rm l})/2 $, $ \bar{\rho} = (\rho_{\rm h}+\rho_{\rm l})/2 $, $ F_{\rm d} $是黏性耗散力, 定义为重力方向上的动量耗散率$ F_{\rm d} = -\epsilon/\nu $, $ \epsilon $是能量耗散率. 对于无黏时间尺度, Sreenivasan^[41]理论分析给出了能量耗散率的数学表达式, $ \epsilon = C \nu^3/W $, 其中C是常数. 根据上述分析, 在雷诺数从小增大的过程中, 黏性耗散力在不断地减少, 从而理论上可以获得更大的尖钉运动速度, 以及更大的尖钉的振幅; 而当雷诺数足够大时, 黏性耗散力已充分小, 浮力相比黏性耗散力更占统治地位, 从而继续增大雷诺数对尖钉振幅增长的影响将不再显著. 在本文模拟中, 我们确实观察到与上述的理论分析相一致的数值结果. 另外, 从图2中还可以发现, 当雷诺数在一定范围内, 增大雷诺数可以有效地促进气泡振幅的增长, 而当雷诺数充分大时, 气泡振幅的增长在后期会随着雷诺数的增大而呈现出一种递减的趋势. 这是由于当雷诺数足够大时, 在浮力驱动的不稳定性演化中后期, 诱导产生了KH不稳定性. 受两类不稳定性的共同作用, 系统中出现了许多不同尺度的旋涡, 这些涡效应在一定程序上减缓了气泡向上运动.
进一步, 我们还定量地统计了不同雷诺数下气泡和尖钉随时间演化的运动速度$ u_{\rm b} $, $ u_{\rm s} $, 其大小是通过气泡和尖钉的振幅对演化时间差分计算获得的, 并根据(25)式无量纲化得到气泡与尖钉的Froude数. 图3描述了不同雷诺数下气泡与尖钉随时间演化的Froude数. 从图3中可以发现, 气泡与尖钉的运动速度在不同雷诺数下表现出显著不同的演化规律. 雷诺数越高, 尖钉往下运动的速度越快, 而当雷诺数充分大时, 尖钉的演化速度对雷诺数的变化不再敏感, 这与上述(26)式理论分析的结果是相统一的. 雷诺数对气泡运动速度的影响则显现出先促进而后抑制的规律, 这是由于当雷诺数充分高时, 在不稳定性的演化后期产生了许多不同尺度的旋涡, 这些旋涡的运动减弱了气泡往上运动的速率而使之发生旋转. Ramaprabhu等^[22]数值研究了混相流体的RT不稳定性的后期演化规律, 他们也发现当雷诺数足够大时, 继续增大雷诺数会减小气泡往上的演化速度. 另外, 我们进一步对气泡与尖钉速度演化图进行分析, 将高雷诺数下非混相流体的单模RT不稳定性的演化归结于四个发展阶段, 包括线性增长阶段、饱和速度增长阶段、重加速阶段、混沌混合阶段. 在初始阶段, 不稳定性的扰动发展符合线性稳定性理论, 其振幅的增长具有指数形式的规律^[16],
图 3 雷诺数对无量纲化的气泡和尖钉演化速度的影响
Figure3. The effect of the Reynolds number on the dimensionless bubble and spike velocities.

$ H = a_1 {\rm e}^{\gamma t}+a_2 {\rm e}^{-\gamma t}, $

其中H代表气泡或者尖钉的振幅, $ a_1 $与$ a_2 $是拟合参数, $ \gamma $是线性增长因子. 图4给出了不同雷诺数下随时间演化的气泡与尖钉振幅的数值模拟结果以及曲线拟合结果, 可以发现气泡和尖钉的初始增长确实符合线性稳定性理论, 并且获得的线性增长因子随着雷诺数的增大而增大. 紧接着线性增长阶段, 气泡与尖钉将以近似恒定的速度增长, 这表明不稳定性的发展进入饱和速度增长阶段. Goncharov^[17]解析分析了单模RT不稳定性的非线性增长区域, 提出了经典的势能理论模型以预测气泡与尖钉的饱和增长速度, 其表达式如(1)式所示. 进一步根据(1)式, 可以推导出气泡与尖钉在饱和速度阶段所对应的无量纲Froude数分别为0.325与0.564. 从图3可以发现, 高雷诺数下尖钉在饱和速度阶段的Froude数略高于势能模型的解析解, 这是由于在实际模拟中, 界面在该阶段发生卷吸行为, 产生了许多不同尺度的旋涡, 这些涡效应会促进尖钉的发展, 而势能模型的理论解未包含涡效应. Goncharov^[17]通过数值模拟同样验证了尖钉实际演化速度高于势能模型的解析解. 另外, 我们发现在高雷诺数条件下, 继续增大雷诺数会减少气泡演化速度, 从而导致气泡饱和速度小于势能模型的解析解. 接下来, 各尺度的涡结构之间相互作用逐渐增加, 使得气泡和尖钉Froude数高于经典势能模型的理论解, 这预示着不稳定发展进入了重加速阶段. 重加速阶段不能持续地发展下去, 在演化后期, 气泡与尖钉的Froude数变得不稳定, 开始随着时间波动, 这表明界面的演化进入了混沌混合阶段. 为了揭示混沌混合阶段不稳定性的发展规律, 我们计算了气泡与尖钉增长的无量纲加速度, $ \alpha_{\rm {b, s}} = \ddot{h}_{\rm {b, s}}/2{A_t}g $, 其中$ \ddot{h}_{\rm {b, s}} $表示气泡与尖钉振幅对时间的两阶导数, 实际通过对气泡与尖钉振幅关于时间的二阶差分计算获得. 图5给出了高雷诺下气泡与尖钉无量纲加速度随时间的演化曲线. 从图可以发现, 气泡与尖钉加速度在演化后期不稳定性, 分别绕着常数0.045与0.233上下波动, 预示着后期气泡与尖钉平均加速度是一个常数, 并表明RT不稳定性的后期发展呈现出二次增长的规律. 当雷诺数足够低时, 在不稳定性的整个演化过程中, 不能观察到重加速阶段与混沌混合阶段, 气泡与尖钉在后期阶段将以恒定的饱和速度增长. 另外, 我们发现低雷诺数下气泡与尖钉的饱和速度低于经典的势能模型的理论解, 其原因在于势能模型考虑的是理想无黏性流体的不稳定性现象, 未考虑流体黏性对演化速度的影响.
图 4 不同雷诺数下, 气泡和尖钉振幅在初始阶段的演化曲线, 其中数据点是统计结果, 实线则是拟合结果
Figure4. The curves of the early-time bubble and spike amplitudes with different Reynolds numbers, where the data points and solid lines are the statistical and fitting results.

图 5 高雷诺数下, 气泡和尖钉的无量纲加速度演化曲线, 红色和蓝色实线分别为0.045和0.233
Figure5. The curves of dimensionless bubble and spike accelerations at a high Reynolds number, and the red and blue solid lines are 0.045 and 0.233.

本文基于相场理论的LB方法模拟了一个长微管道内非混相流体的RT不稳定性问题, 该方法可以准确地追踪相界面动力学行为, 以及采用多松弛碰撞模型可以很好地处理高雷诺数流动. 本文着重分析雷诺数对中等Atwoods数的不稳定性演化中界面动力学行为和扰动增长规律的影响. 研究发现在高雷诺数情形下, RT不稳定性的发展先后经历线性增长阶段、饱和速度增长阶段、重加速阶段以及混沌混合阶段. 在线性阶段, 重流体与轻流体之间相互渗透, 扰动的增长符合经典的线性稳定性理论. 紧接着, 不稳定性的演化进入了饱和速度阶段. 在该阶段中, 轻重流体形成了气泡与尖钉的结构, 气泡与尖钉将以恒定的速度增长, 并随着时间演化在尖钉尾端处出现卷吸行为. 随着横向速度和纵向速度的差异逐渐扩大, 气泡与尖钉的演化诱导产生了KH不稳定性, 进而使系统中出现了许多不同尺度的涡结构, 从而加速了气泡与尖钉的演化速度. 重加速阶段不能一直持续下去, 气泡与尖钉的演化速度在后期阶段会发生波动现象, 这表明不稳定性的发展进入了混沌混合阶段. 在混沌混合阶段, 界面会发生多层次卷起、剧烈变形、混沌破裂, 形成了非常复杂的界面拓扑结构, 系统中也产生了许多离散的小液滴. 另外, 我们统计了演化后期的气泡与尖钉加速度, 其无量纲的值分别围绕着常数0.045与0.233上下波动, 这表明不稳定性在演化后期具有二次增长的规律. 而当雷诺数足够小, 扰动的增长变得非常缓慢, 流体相界面也变得足够光滑, 未出现卷吸和破裂行为, 也未观察到后期的重加速与混沌混合阶段, 气泡与尖钉将以恒定速度一直演化下去.