摘要:采用多相流的相场格子Boltzmann方法数值研究了微通道内高雷诺数单模Rayleigh-Taylor (RT)不稳定性的后期演化规律, 重点分析表面张力对相界面动力学行为以及气泡与尖钉增长的影响. 数值实验表明, 随着界面张力的增大, 可以有效降低演化过程中相界面结构的复杂程度, 并抑制不稳定性后期相界面破裂形成离散液滴. 另外, 增大表面张力可以先促进后抑制气泡振幅的增长, 而当表面张力较小时, 尖钉振幅增长曲线之间并无明显差别, 当表面张力增大到一定值后, 它对尖钉振幅的抑制效果可明显地被观察到. 进一步, 根据不稳定性速度增长曲线, 将高雷诺数单模RT不稳定性的演化划分为线性增长、饱和速度增长、重加速、混沌混合四个发展阶段. 数值计算获取气泡与尖钉的饱和速度符合包含界面张力效应的势流理论模型. 另外还统计了不同表面张力和Atwood数下表征RT不稳定性后期演化的气泡与尖钉增长率, 结果显示气泡与尖钉后期增长率随着表面张力的增大总体上呈现出先促进后抑制的规律. 最后, 从数值计算和理论分析两方面研究了不同Atwood数下RT 不稳定性发生的临界表面张力, 发现两者结果符合得很好, 并且临界表面张力随着流体Atwood数的增大而增大.
关键词: Rayleigh-Taylor不稳定性/
格子Boltzmann方法/
表面张力/
后期增长/
临界值

English Abstract

全文HTML

--> --> -->

瑞利-泰勒(Rayleigh-Taylor, RT)不稳定性现象是一种常见的流体不稳定性问题, 它发生于重力场作用下, 两种不同密度的流体界面受到微小扰动后. 这类问题不仅存在于自然现象中, 例如卷云的形成, 还在天体物理学、惯性约束聚变、地球气候学等领域有着广泛且重要的应用, 也是多相流体界面、不稳定性、湍流混合等基本问题的理论基础, 因而吸引着许多****对其开展相关研究^[1]. RT不稳定性的先驱工作归功于英国流体力学大师 Rayleigh和Taylor^[2,3], 他们创造性地提出了线性稳定性理论, 并证明了无黏流体的初始扰动呈现指数增长的规律. 后来, Lewis^[4]实验验证了线性增长理论, 还报告了RT不稳定性随后进入非线性发展阶段. Sharp^[5]对RT不稳定性的早期研究工作做了很好的总结, 并定性地将RT不稳定性的增长划分为四个发展阶段. 根据界面初始扰动的不同, RT不稳定性可以分为单模和多模问题, 而本文主要关注单模RT不稳定性的研究. 2001年, Waddell 等^[6]实验研究了二维单模 RT不稳定性问题, 指出扰动振幅在初始阶段以指数形式增长, 即 $h = a_1 {\rm e}^{\gamma t}+ $$ a_2 {\rm e}^{-\gamma t}$, 其中 $ \gamma $ 是线性增长因子, 定义为

$ \gamma = \sqrt{A_t g k-\frac{\sigma k^3}{\rho_{\rm h}+\rho_{\rm l}}}, $

其中$ g $ 是重力加速度, $ A_t $是Atwood数, $ k $是波数, $ \sigma $是流体间的表面张力, $ \rho_{\rm h} $和$ \rho_{\rm l} $分别为重轻流体的密度. 进一步, 他们还观察到在不稳定性演化后期气泡与尖钉的平均速度接近于常数. Glimm等^[7]利用前追踪方法模拟了二维单模RT不稳定性问题, 发现不稳定性的增长经历饱和速度阶段后会出现重加速阶段. 后来, Wilkinson和Jacobs^[8]对三维单模RT不稳定性进行了实验研究, 指出由于界面涡结构间的相互作用, 气泡与尖钉在后期阶段被重加速, 导致演化速度高于Goncharov势能模型^[9]的理论值. 2012 年, Ramaprabhu等^[10]大规模模拟了不同Atwood数和雷诺数下三维单模RT不稳定性的后期演化问题, 首次将高雷诺数RT不稳定性的发展归结为线性增长、饱和速度增长、重加速、混沌混合四个发展阶段, 并进一步观察到气泡与尖钉在重加速阶段的演化速度超过经典势流模型^[9]的解析解, 而在后期的混沌混合阶段, 发现气泡与尖钉的增长会受到抑制. Wei和Livescu^[11]利用直接数值模拟方法高精度研究了低Atwood数下二维单模RT不稳定性问题, 观察到不稳定性在高雷诺数条件下同样经历了四个发展阶段, 但不同于Ramaprabhu等^[10]的结果, 他们发现气泡的平均振幅在后期混沌阶段具有二次增长的规律. Lai等^[12]采用离散 Boltzmann 方法研究了RT不稳定性的热动非平衡效应. 近期, Liang等^[13-15]采用介观格子Boltzmann (lattice Boltzmann, LB)方法模拟了二维及三维微通道内单模RT不稳定性后期演化问题, 分析了Atwood数和雷诺数对不稳定性的增长阶段和相界面动力学行为的影响.
上述研究均忽略了表面张力对RT不稳定性的影响, 已有研究表明RT不稳定性在表面张力作用下可以显示毛细波、收缩、破裂等独特的界面动力学行为^[16]. 目前, 绝大多数的工作^[16-20]集中于分析表面张力对初始扰动为多模态的RT不稳定性增长, 而较少****考虑表面张力对单模RT不稳定性演化的影响. Daly^[21]采用有限差分法数值研究了表面张力对单模RT不稳定性早期增长的影响, 获取的不同界面张力下线性增长因子与修正的线性稳定性理论一致, 并且发现存在表面张力时可以抑制界面的卷吸行为. Zhang等^[22]利用多相流LB方法模拟了单模RT不稳定性问题, 发现表面张力作用下可以减弱相界面的卷吸程度, 在演化后期诱导相界面发生夹断生成离散液滴, 并且观察到界面张力存在时可以减小气泡与尖钉的演化速度. 基于Goncharov势流理论^[9], Young和Ham^[23]分析界面张力对无黏流体的单模RT不稳定性增长的影响, 并提出了包含表面张力效应的气泡渐进速度的解析模型,

$ {u_{\rm b}} = \sqrt{\frac{2A_t g}{3k(1+A_t)}-\frac{\sigma k}{9\rho_{\rm h}}}.$

进一步, 他们还利用水平集数值方法模拟了表面张力作用下单模RT不稳定性在饱和速度阶段的增长, 发现获得的气泡渐进速度与提出的解析模型相一致. 同样从势流理论出发, Sohn^[24]分析了流体黏性和表面张力对单模RT不稳定性发展的影响, 发现黏性和表面张力均抑制不稳定性的发展, 减小气泡的饱和增长速度. 结合黏性效应, 气泡在饱和速度阶段的渐进速度修正为

$ {u_{\rm b}} = \sqrt{\frac{2A_t g}{3k(1+A_t)}+\frac{4}{9}k^2 \nu_{\rm h}^2-\frac{\sigma k}{9\rho_{\rm h}}}-\frac{2}{3}k\nu_{\rm h}, $

其中$ \nu_{\rm h} $ 表示重流体的运动学黏性. 夏同军等^[25]分别基于势流模型和Zufiria模型^[26]分析了表面张力对RT不稳定性非线性增长阶段的影响, 推导出两种模型下气泡渐进速度和渐进曲率的解析表达式, 结果表明, 界面张力降低了气泡的速度, 但对曲率没有影响. 另外发现, 基于势流模型所得气泡渐进速度大于Zufiria模型的预测值. 进一步, 基于Zufiria理论模型, Li等^[27]分析了流体黏性和表面张力对不稳定性饱和速度阶段的影响, 提出了包含两种效应的气泡渐进速度的解析式. Guo等^[28]建立了含表面张力的单模RT不稳定性的弱非线性模型, 分析了RT不稳定性在表面张力作用下从线性增长到非线性增长的转变.
综上所述, 已有****分析了表面张力对单模RT不稳定性增长的影响, 但这些研究聚焦于不稳定性中气泡增长的中前期阶段, 而对表面张力作用下不稳定性的演化后期及尖钉增长的描述还未有报道. 另外, 上述研究所考虑的雷诺数均较小. 本文将采用基于相场理论的多松弛LB方法研究冗长微通道内高雷诺数RT不稳定性的演化规律, 着重分析表面张力对相界面动力学行为及气泡与尖钉后期增长的影响.

LB方法^[29]是近几十年发展起来的一种流体系统建模和模拟的介观方法, 因其具有易于刻画流体间与流固间的微观相互作用、高性能并行计算及自动追踪相界面等优点, 因而被广泛用来研究多相流界面动力学问题. 目前, 从流体间相互作用力的不同物理背景, 已提出了多种不同类别的多相流LB模型, 包括颜色模型、伪势模型、自由能模型和相场LB模型. 相比前三类模型, 相场LB模型因在相界面追踪方面具有坚实的物理学基础而受到一些****的广泛关注, 感兴趣的读者可以参考近期关于多相流体的相场LB方法的综述论文^[30]. 相场模型中界面追踪方程包括Cahn-Hilliard方程和Allen-Cahn方程, 基于Allen-Cahn模型的LB方法在求解序参数时具有更低的数值耗散, 适用于模拟大密度比两相流问题^[31], 而基于Cahn-Hilliard模型的LB方法在模拟界面大拓扑变化时则更具有优势. 本文采用Liang等^[13]提出的基于Cahn-Hilliard方程的相场LB模型, 该模型克服了已有模型中相界面求解精度低以及数值稳定性差等问题, 提高了LB方法在相界面捕获方面的能力, 可显式计算流场中的宏观量, 并在求解复杂的流体界面RT不稳定性问题上显示出较大的潜力^[14,32]. 在相场理论中, 多相流体系统的总自由能泛函可表述为^[33]

$ E = \int_{\varOmega}\left[F(\phi)+\frac{\alpha}{2} {\left| {\nabla\phi} \right|}^2\right]\mathrm{d}\varOmega, $

其中$ \phi $是序参数, 用于捕捉相界面; $ F(\phi) $是体相区的自由能密度泛函, 具有双势阱形式$F(\phi) = \beta\phi^2(1- $$ \phi)^2$, 参数$ \alpha $和$ \beta $依赖于表面张力$ \sigma = \sqrt{2\alpha \beta/6} $和界面厚度$ D = \sqrt{8\alpha/\beta} $. 对总的自由能泛函求变分运算, 可以获得化学势的表达式为

$ \mu = 4\beta\phi(\phi-1)(\phi-0.5)-\alpha\nabla^2\phi. $

进一步, 假设序参数$ \phi $由流体速度$ {u} $发生对流, 和化学势$ \mu $产生扩散, 从而得到描述相界面动力学的经典Cahn-Hilliard相场方程:

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

其中, $ M $是迁移率. 为了描述多相流体输运过程, 相界面追踪的Cahn-Hilliard方程需要耦合含外力的不可压Navier-Stokes方程:

$ \nabla \cdot {{u}} = 0, \tag{7a}$

$\begin{split}\rho\left({{\partial {{u}}} \over {\partial t}} + {{u}}\cdot\nabla {{u}}\right) = \;&-\nabla p + \nabla\cdot\left[ {\nu\rho(\nabla {{u}} + \nabla{{{u}}^{\rm T}})} \right]\\&+ {{{F}}_{\rm s}} + {{G}},\end{split} \tag{7b}$

其中$ \rho $是流体密度, $ p $是流体动力学压力, $ \nu $是运动学黏性, 界面张力$ {{F}}_{\rm s} $取势形式为$ {F}_{\rm s} = \mu\nabla{\phi} $, $ {G} $是系统中施加的外力.
本文的LB模型利用两个分布函数$ f_i $和$ g_i $分别来求解Cahn-Hilliard和Navier-Stokes相互耦合的多相系统控制方程, 对应的多松弛碰撞算子的演化方程可表示为^[29]

$\begin{split}&{f}_i({{x}} +{{{c}}_i}{\delta _t},t + {\delta _t}) -{f}_i({{x}},t) \\=\;& -({{M}}^{-1}{{S}}^{f}{{M}})_{ij}[{{f}_j({{x}},t) -f_j^{\rm eq}({{x}},t)}]\\&+ {\delta _t}\left[{M}^{-1}\left({I}-\frac{{S}^f}{2}\right){M}\right]_{ij}F_j({{x}},t),\end{split} \tag{8a}$

$\begin{split}&{g}_i({{x}} +{{{c}}_i}{\delta _t},t + {\delta _t}) -{g}_i({{x}},t) \\=\;& -({{M}}^{-1}{{S}}^{g}{{M}})_{ij}[{{g}_j({{x}},t) -g_j^{\rm eq}({{x}},t)}]\\&+ {\delta_t}\left[\mathbf{M}^{-1}\left({I}-\frac{{S}^g}{2}\right){M}\right]_{ij}G_j({{x}},t),\end{split} \tag{8b}$

其中, $ {f}_i({{x}}, t) $, $ {g}_i({{x}}, t) $是粒子分布函数, $ f_i^{\rm eq}({{x}}, t) $和$ g_i^{\rm eq}({{x}}, t) $是分布函数对应的平衡态分布函数, $ {{c}}_i $是离散速度, $ {{M}} $表示碰撞矩阵, $ {{S}}^f $和$ {{S}}^g $是对角松弛矩阵, $ F_i({{x}}, t) $和$ G_i({{x}}, t) $分别为源项和外力项的分布函数. 为了恢复正确的宏观方程, 平衡态分布函数可设定为^[13]

$ f_i^{\rm eq} = \left\{ \begin{aligned}&\phi+({\omega_i} - 1)\eta\mu, \quad { i = 0}, \\&{\omega_i}\eta\mu+{\omega_i}\phi{{{{{c}}_i} \cdot {{u}}}\over {c_{\rm s}^2}}, \quad { i\neq0}, \end{aligned} \right. $

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

且

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

其中, $ \eta $是调节迁移率的自由参数, $ c_{\rm s} $是声速, $ \omega_i $ 是权系数. 对于二维流动, 本文采用流行D2Q9格子模型, 其对应的权系数$ \omega_i $为$ \omega_0 = 4/9 $, $\omega_{1-4} = $$ 1/9$, $ \omega_{5-8} = 1/36 $, 离散速度 $ {c}_{i} $ 定义为^[34,35]

$ { c}_{i} = \left\{ \begin{aligned}& (0,0)c, \quad\;\qquad\qquad\qquad\qquad\qquad\quad {i = 0}, \\&(\cos [(i-1)\pi /2],\sin [(i-1)\pi /2])c,\quad { i = 1—4}, \\&\sqrt{2}\Big(\cos [(i-5)\pi /2+\pi /4],\\&\;\sin [(i-5)\pi /2+\pi /4]\Big)c, \qquad\qquad \quad{ i = 5—8}, \end{aligned} \right. $

其中, $ c = \delta_x/\delta_t $, $ \delta_x $和$ \delta_t $分别代表格子长度和时间步长, $ c_{\rm s} = c/\sqrt{3} $, 对于多相流LB模型, 通常设定$ c = \delta_x = \delta_t = 1 $. 另外, D2Q9 模型所对应的碰撞矩阵 $ {M} $ 定义为^[36]

$\begin{array}{l} {{M}} \!=\!\! \left(\!\!\!\! \begin{array}{rrrrrrrrr} 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}$

以及松弛矩阵$ {{S}}^f $和$ {{S}}^g $可表示为

$ {{S}}^f \!=\! {\rm diag}\left(s_0^f,\,s_1^f,\,s_2^f,\,s_3^f,\,s_4^f,\,s_5^f,\,s_6^f,\,s_7^f,\,s_8^f\right), $

$ {{S}}^g = {\rm diag}\left(s_0^g,\;s_1^g,\;s_2^g,\;s_3^g,\;s_4^g,\;s_5^g,\;s_6^g,\;s_7^g,\;s_8^g\right), $

其中, $ 0 < s_i^f < 2, \; 0 < s_i^g < 2 $. 为了准确地恢复宏观控制方程, 离散源项$ F_i $和外力项分布函数$ G_i $分别定义为^[13]

$ F_i = \frac{\omega_i{c}_{i}\cdot\partial_t(\phi{u})}{c_{\rm s}^2}, $

$ {{G}_i} \!=\! \frac{({c}_i-{u})}{c_{\rm s}^2}\cdot[s_i({u})\nabla(\rho{c_{\rm s}^2})+({F}_{\rm s}+{F}_a+{G})(s_i({u})+\omega_i)], $

其中, ${F}_a = \dfrac{{\rm d}\rho}{{\rm d} \phi}M\nabla^2\mu {u}$是为了消除恢复的动量方程中误差项所引入的界面力. 通过计算分布函数的零阶矩和一阶矩, 可以得到宏观量$ \phi $, $ {u} $和$ p $, 具体计算公式为

$ \phi = \sum\limits_i {{f_i}}, \tag{17a}$

$ \rho{{u}} = {{\sum\limits_i {{{{c}}_i}{g_i}}+ 0.5{\delta_t}({{{F}_{\rm s}}+{G}})}}, \tag{17b}$

$ p = {{c_{\rm s}^2} \over {(1 - {\omega_0})}}\left[ {\sum\limits_{i \ne 0} {{g_i}} +{{{\delta _t}} \over2}{{u}} \cdot \nabla \rho + \rho {s_0}({ u})} \right].\tag{17c} $

另外, 流体密度$ \rho\; $可看成序参数$ \phi $的线性函数:

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

通过Chapman-Enskog多尺度分析^[13], 可以证明本文采用的相场LB模型能够正确地恢复到Cahn-Hilliard方程和Navier-Stokes方程, 并且迁移率$ M $和运动学黏度$ \nu $的数学表达式为

$ M = \eta c_{\rm s}^2(\tau_f-0.5)\delta_t,\tag{19a}$

$ \nu = c_{\rm s}^2(\tau_g-0.5)\delta_t, \tag{19b} $

其中, $ 1/\tau_f = {s_3^f} = {s_5^f} $, $ 1/\tau_g = {s_7^g} = {s_8^g} $. 已有研究表明^[36,37], 松弛参数的选取影响多松弛LB方法的数值稳定性和精度. 在实际的模拟中, 松弛因子$ {s_3^f} $和$ {s_7^g} $可根据迁移率和流体黏性给定, $ {s_4^g} = {s_6^g} = $$ 1.7 $, 余下的松弛参数均取 1. 另外, 数值实验表明, 松弛参数在一定温和的范围内取值, 本文的多松弛LB方法可以给出相当的计算结果. 演化方程中时间导数采用显式欧拉格式进行离散, 即$\partial_t(\phi{u}) = $$ [\phi(t){u}(t)-\phi(t-\delta_t){u}(t-\delta_t)]/\delta_t$, 而梯度和拉普拉斯算子的离散方法均采用二阶中心各向同性差分格式^[38]:

$ \nabla \phi({{x}}) = \sum\limits_{i \ne 0}{\frac{\omega_i{c}_i\phi({{x}} +{{{c}}_i}{\delta_t})}{c_{\rm s}^2 \delta_t}}, $

$ \nabla^2\phi({{x}}) = \sum\limits_{i \ne 0} {\frac{2\omega_i[\phi({{x}} +{{{c}}_i}{\delta_t})-\phi({{x}})]}{c_{\rm s}^2 \delta_t^2}}, $

而密度梯度的计算公式为$\nabla\rho = \dfrac{{\rm d}\rho}{{\rm d}\phi}\nabla\phi$.

本文关注冗长的微管道内多相流体界面单模扰动的RT不稳定性, 网格划分为 $L\times W = 8160 \times $$ 544$. 初始时刻, 在微管道的中心截面处施加波长为$ W $的微小扰动:

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

其中, 波数$ k $定义为$ k = 2\pi/W $, 并设定序参量的初始分布为

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

其中, 界面厚度$ D $固定为$ 4 $个格子网格. 需要指出的是, RT不稳定性的界面演化可以用两个无量纲参数来描述, 即雷诺数(Reynolds number, $ Re $)和阿特伍德数(Atwood number, $ A_t $), 分别定义为:

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

本文着重分析高雷诺数RT不稳定性的后期演化规律, 因而在数值模拟中雷诺数均设定为 $ Re = 10^4 $. 为了表征重力效应, 对微管道中轻重流体均施加一个竖直方向上的浮力:

$ {{G}} = \left[0, - \left(\rho-\frac{\rho_{\rm h}+\rho_{\rm l}}{2}\right)g\right], $

其中, 重流体密度设定为$ \rho_{\rm h} = 1 $, 轻流体的密度可由给定的阿特伍德数确定, 重力加速度满足关系式$ \sqrt{gW} = 0.04 $. 另外, Peclet数 ($ Pe $)设定为50, 其定义式为$ Pe = \sqrt{gW}D/M\beta $. 边界条件设置如下: 左右边界采用周期性边界条件, 而上下边界均采用无滑移半反弹边界格式. 半反弹格式将壁面设定在格点中线上, 假定粒子与壁面碰撞后速度发生逆转, 即临近壁面流体点$ {x}_f $的粒子分布函数设定为$ f_{\bar{i}}({x}_f, t+\delta_t) = f'_i({x}_f, t) $, 其中$ {c}_i = -{c}_{\bar{i}} $是指向壁面的粒子速度, $ f'_i $表示碰撞后的粒子分布函数^[39]. 具体而言, 当$ {{x}_f} $位于上壁面时, 半反弹边界格式实现的数学表达式为

$\begin{split} f_4({x}_f,t+\delta_t) =\;& f'_2({x}_f,t),\; f_7({x}_f,t+\delta_t) \\=\;& f'_5({x}_f,t),\; f_8({x}_f,t+\delta_t) \\= \;&f'_6({x}_f,t), \end{split}$

而当$ {{x}_f} $为下壁面时, 半反弹边界条件可设置为

$\begin{split}f_2({x}_f,t+\delta_t) =\;& f'_4({x}_f,t),\; f_5({x}_f,t+\delta_t)\\=\;& f'_7({x}_f,t),\; f_6({x}_f,t+\delta_t) \\=\;& f'_8({x}_f,t). \end{split}$

在本文的研究中, 选取特征长度为扰动的初始波长$ W $, 特征速度和特征时间分别给定为$\sqrt{gW}$和$ \sqrt{W/g} $, 下文统计的相关物理量均已被相应特征值所无量纲化.
图1给出了$ A_t = 0.7 $时, 四种不同表面张力下高雷诺数RT不稳定性的相界面演化过程. 可以发现, 对于不同的表面张力, 相界面在初始阶段表现出相似的动力学行为: 重流体和轻流体相互渗入形成了尖钉和气泡结构. 随后在不同表面张力下, 尖钉和气泡的增长表现出明显的差异. 在表面张力很小, 仅为$ 5\times10^{-6} $时, 尖钉持续下降并出现非线性效应, 这是由于气泡和尖钉界面上存在的速度差诱导产生了非线性 Kelvin-Helmholtz (KH)不稳定性, 进而使尖钉前端卷起形成一对涡. 随着两个旋涡的持续增长, 在卷起的尾端出现了二级涡, 并且混合区域内的非线性效应也越来越剧烈, 导致界面在多个位置发生卷起, 形成了复杂的界面拓扑结构. 在不稳定性继续发展过程中也伴随着流体界面混沌的破裂, 最终导致混合区域存在着大量的离散液滴, 同时也观察到相界面的增长失去了左右对称性. 相界面的非对称性发展同样在高雷诺数下混相RT不稳定性现象中被观察到^[40]. 随着表面张力的增加, KH不稳定性的出现开始延后, 一级涡的产生相应地延缓, 随后产生的涡的数量也相应地减少, 后期的结构复杂度下降. 不同于$ \sigma = 5\times10^{-6} $情形, RT不稳定性在较大界面张力的增长过程中, 相界面图案始终保持着左右对称性, 整个演化过程中界面破裂的现象减少, 系统中离散的液滴数量也相应减少, 尤其当$ \sigma = 1\times10^{-2} $时, 不稳定性的发展过程中未观察到离散液滴. 另外, 我们通过数值实验进一步发现, 继续增大表面张力, 扰动会随时间发生振荡并逐渐恢复平稳, RT不稳定现象并未发生, 这表明存在一个临界的表面张力的值, 超过临界值后, RT不稳定现象不会发生, 而小于该临界值, RT 不稳定性将会发生. 下文将细致讨论临界表面张力的影响因素以及关系表达式.
图 1 $ Re = 10^4 $, $ A_t = 0.7 $时, 表面张力对高雷诺数非混相RT不稳定相界面演化的影响　(a) $ \sigma = 5\times10^{-6} $; (b) $\sigma = $$ 5\times10^{-4}$; (c) $ \sigma = 5\times10^{-3} $; (d) $ \sigma = 1\times10^{-2} $
Figure1. Effect of the surface tension on the evolution of phase interface in the immiscible RT instability with $ Re = 10^4 $, $ A_t = 0.7 $: (a) $ \sigma = 5\times10^{-6} $; (b) $ \sigma = 5\times10^{-4} $; (c) $ \sigma = 5\times10^{-3} $; (d) $ \sigma = 1\times10^{-2} $.

气泡和尖钉的振幅是单模RT不稳定研究中两个重要的物理量, 因而, 本文统计了不同表面张力条件下随时间变化的气泡与尖钉振幅, 以定量比较界面张力的影响. 图2为气泡与尖钉振幅在不同表面张力下随时间的演化曲线. 可以看出, 对于气泡而言, 随着表面张力的增加, 气泡振幅的增长并不是一直减缓的, 在表面张力很小的一段范围内, 表面张力反而促进了气泡的增长. 当$ \sigma = 5\times10^{-4}\; $时, 气泡振幅的增长在发展前期受表面张力影响不大, 但是在后期, 增长明显快于$ \sigma = 5\times10^{-6} $和$ \sigma = 5\times10^{-5} $时的气泡振幅. 而表面张力增加到$ \sigma = 5\times10^{-3} $后, 气泡振幅的增长无论是前期还是后期均慢于$ \sigma = 5\times10^{-4} $时的增长. 继续增大表面张力到$ \sigma = 1\times10^{-2} $, 在前中期的抑制效果更加明显. 上述结果与文献[23,24]中普遍认为表面张力抑制气泡的增长有所不同, 原因可能是前人所考虑的界面张力范围有限. 而对于尖钉而言, 当表面张力相对较小时, 振幅曲线没有明显的差异, 这表明当表面张力足够小时, 其对尖钉振幅的影响不再重要. 但是当表面张力增大到一定值后, 它对尖钉振幅的抑制效果则可以被明显观察到. 我们还进一步模拟了不同表面张力条件下中低 Atwood数$(A_t = $$ 0.5, \; 0.1)$的RT不稳定性问题, 结果表明, 随着界面张力的增大, 表面张力对于气泡振幅的影响同样是先促进后抑制的, 说明这是一个普遍的规律, 而表面张力可以减缓尖钉的增长.
图 2 表面张力对随时间变化的气泡和尖钉振幅的影响
Figure2. Influence of surface tension on the time variations of bubble and spike amplitudes.

接着对气泡和尖钉的增长速度进行定量分析. 图3给出了在不同表面张力下气泡与尖钉速度随时间变化的演化曲线. 根据文献[11,14], 高雷诺数下的单模RT不稳定性的发展经历了四个不同的阶段, 包括线性增长阶段、饱和速度阶段、重加速阶段和混沌阶段. 当表面张力足够小时, 气泡和尖钉速度在线性阶段的发展十分相似, 增长曲线基本重合, 即表面张力的效应不显著. 而当表面张力增大到一定值后, 线性阶段的气泡与尖钉速度减缓明显. 上述数值结果与线性稳定性理论^[6]分析一致: 扰动振幅在初始阶段的增长具有指数形式, 其线性增长因子如(1)式所示, 这从理论上表明当表面张力非常小时, 其值对初始阶段扰动增长的影响可以忽略, 而当表面张力较大时, 表面张力的效应逐渐显著, 增大其值可以有效地减弱扰动的线性发展. 紧接着线性增长阶段, 从图3可以发现气泡与尖钉的增长将以近似恒定的速度增长, 这表明不稳定性的发展进入饱和速度阶段. Goncharov^[9]基于经典的势流理论模型预测了无黏流体在忽略表面张力条件的单模RT不稳定性中气泡与尖钉的饱和速度, 其无量纲形式分别表示为
图 3 表面张力对随时间演化的气泡和尖钉增长速度的影响. 蓝色和虚线分别表示势能模型预测气泡与尖钉速度在 $ \sigma = 5\times10^{-3} $ 和 $ \sigma = 1\times10^{-2} $ 时的解析解
Figure3. Influence of surface tension on the time evolutions of bubble and spike growth velocities. The blue and yellow dotted lines represent the analytical solutions of the bubble and spike velocities from potential flow model at $ \sigma = 5\times10^{-3} $ and $ \sigma = 1\times10^{-2} $.

$ {u_{\rm b}} = \sqrt{\frac{A_t}{3\pi(1+A_t)}},\; \; \; {u_{\rm s}} = \sqrt{\frac{A_t}{3\pi(1-A_t)}}. $

Sohn^[24]分析流体黏性和表面张力对不稳定性增长的影响, 并基于Goncharov的解析势流模型, 提出了修正的气泡渐进速度的数学表达式((3)式). 受Goncharov^[9]工作的启发, 耦合黏性和表面张力效应的尖钉恒定速度的无量纲形式可以写为

$ {u_{\rm s}} = \sqrt{\frac{A_t}{3\pi(1-A_t)}-\frac{2\pi}{9Bo}+\left(\frac{4\pi}{3Re}\right)^2}-\frac{4\pi}{3Re}, $

其中Bo为邦德数, 定义为$ Bo = \rho_{h}gW^2/\sigma $. 根据 Sohn^[24]的解析模型, 我们发现气泡与尖钉满足 $ \sigma\leqslant5\times 10^{-3} $时渐进速度的理论解几乎差别不大, 因而在图3中仅给出$ \sigma = 5\times10^{-3} $和$ \sigma = 1\times10^{-2} $条件下气泡与尖钉饱和速度的理论值. 可以看出, 气泡渐进速度的数值计算结果与理论解符合较好, 并且发现在界面张力较大的情形下, 增大其值可以有效地降低气泡的渐进速度. 而对于尖钉而言, 界面张力对尖钉饱和速度的解析解影响不大, 数值计算的尖钉渐进速度在表面张力较小时与理论解一致, 而在表面张力较大时, 数值预测的结果略低于解析解, 这表明界面张力较大时可以轻微地减弱尖钉在饱和速度阶段的增长. 随着不稳定性的非线性强度持续增强, 气泡与尖钉的演化速度会超过各自对应的饱和速度理论值, 这说明不稳定性的发展进入了重加速阶段. 从图3可以观察到, 气泡的再加速阶段在小表面张力($ \sigma = 5\times $$ 10^{-6} $)时并不明显, 在增大表面张力至$ \sigma = 1\times10^{-2} $后, 重加速阶段则十分清晰. 最后, 气泡与尖钉速度随时间上下波动, 显示不稳定的行为, 其值会被反复加速和减速, RT不稳定性的演化进入混沌混合阶段. 从图3可以看出, 增大表面张力可以有效地抑制气泡与尖钉速度在混沌混合阶段的波动程度. 另外, 已有研究表明, 高雷诺数的单模RT不稳定性的增长在后期的混沌混合阶段呈现出平均二次增长的规律^[11,14], 即气泡与尖钉演化振幅随时间变化的函数可表示为$ {{h}}_{\rm {s, b}} = \alpha_{\rm {s, b}}{A_t}gt^2 $, 其中$ \alpha_{\rm {s, b}} $分别表示气泡与尖钉的后期增长率系数. 值得注意的是, 气泡与尖钉增长率是RT不稳定性研究中最为关心且重要的参数, 目前已提出了多种不同方法用于测量后期阶段的气泡与尖钉增长率. 为了显示良好的收敛性, 本文采用Cabot等^[41]提出的测量气泡和尖钉的后期增长率的方法, $\alpha_{\rm {s, b}} = $$ {\dfrac{{{\dot{h}}_{\rm {s, b}}}^2}{4{A_t}g{h_{\rm {s, b}}}}}$, 其中$ \dot{h}_{\rm {s, b}} $表示气泡与尖钉振幅对时间的一阶导数. 图4给出了不同Atwood数下气泡与尖钉的后期增长率与表面张力的关系曲线. 可以看出, 相同的Atwood数下, 尖钉的后期增长率远大于对应的气泡后期增长率, 其效应随着Atwood数的增长而越发显著. 而对于相同的界面张力, 尖钉后期增长率会随着Atwood数的增大而显著递增. 另外还发现, 尖钉与气泡的增长系数随着表面张力的增大总体上呈现出先递增而后减小的趋势.
图 4 不同Atwood数下气泡与尖钉的后期增长率系数 $ \alpha_{\rm {b, s}} $ 随表面张力的变化
Figure4. Variations of bubble and spike late-time growth coefficients $ \alpha_{\rm {b, s}} $ with respect to the surface tension under different Atwood numbers

最后讨论RT不稳定性现象发生的临界表面张力. 通过大规模的数值实验, 发现当表面张力小于某个临界阈值时, RT不稳定性现象将会发生, 而表面张力大于该临界值时, 相界面是稳定的, RT不稳定性受到抑制. 另外还发现, 临界表面张力的值依赖于流体Atwood数. 临界表面张力的范围可以通过半分法来确定. 对于固定流体Atwood数, 考虑一个较大界面张力的范围$ [\sigma_0\; \; \sigma_1] $, 其中$ \sigma = \sigma_0 $时, RT不稳定性能发生, 而当$ \sigma = \sigma_1 $时, RT不稳定性不能发生. 接着, 考虑$\sigma = ({\sigma_0+\sigma_1})/{2}$的情形, 如果通过数值模拟发现 RT不稳定性能发生, 则临界表面张力坐落在$\left[({\sigma_0+\sigma_1})/{2}\; \; \sigma_1\right]$范围内, 否则临界表面张力的范围为$\left[\sigma_0\; \; ({\sigma_0+\sigma_1})/{2}\right]$. 反复重复上述的数值实验, 最终可以确定临界表面张力的数值预测的变化范围. 另外, 根据线性稳定性理论, RT不稳定性在初始阶段的线性增长因子如(1)式所示. 要使RT不稳定性能够发生, 必须保证线性增长因子中开根号的值不小于 0, 因而可以推导出临界表面张力的理论关系式为

$ \sigma_c = \frac{(\rho_{\rm h}-\rho_{\rm l})g}{k^2}, $

这与利用文献[24,27]提出的模型推导出的结果是相同的. 图5给出了LB方法数值计算获得的不同Atwood数下的临界表面张力以及理论模型所预测的临界值. 从图5可以发现, 临界表面张力会随着Atwood数的增加而增大, 并且理论给出临界值基本落在数值模拟所预测的范围内, 即数值计算结果与理论模型结果一致, 这也表明利用(30)式对临界表面张力进行计算是合理的.
图 5 不同 Atwood 数下的临界表面张力
Figure5. Critical surface tensions at various Atwood numbers.

本文基于相场理论的LB方法模拟了冗长微通道内非混相流体的单模RT不稳定性问题, 考察了高雷诺数条件下表面张力对相界面动态过程及扰动增长的影响. 数值结果表明, Atwood数为0.7的单模RT不稳定性在不同的表面张力下表现出明显不同的界面动力学特征: 低表面张力时, 不稳定性演化后期流体界面出现混沌的破裂, 诱导系统中产生大量的离散液滴, 同时也观察到相界面在后期演化过程中显示非对称性特征; 而随着表面张力逐渐增大, 界面区域中产生离散液滴的数量减少, 在极大的表面张力作用下相界面甚至未出现破裂情形, 在整个演化过程中界面始终保持关于中轴线对称. 另外还分析了表面张力对气泡与尖钉振幅增长和演化速度的影响. 在表面张力较小时, 表面张力对尖钉增长的作用不显著, 继续增大表面张力可以明显地抑制尖钉的增长, 而表面张力对气泡增长的影响则呈现出先促进后抑制的规律. 进一步统计了不同表面张力和Atwood数下气泡与尖钉演化后期的二次增长系数, 可以发现对于固定的Atwood数, 尖钉的后期增长率远大于相应的气泡增长率, 并且其效应随着Atwood数的增大而更加显著. 此外还发现, 气泡与尖钉的后期增长率随着表面张力的增大总体上呈现先增大而后递减的规律. 最后, 通过理论分析建立了RT不稳定性现象发生的临界表面张力的数学表达式, 揭示了临界表面张力与流体Atwood数呈现递增关系, 并且通过LB方法计算获取的临界表面张力与理论模型结果相符合.