0 引　言

在自然界和工程技术领域中存在三种常见的流体不稳定性：瑞利-泰勒（Rayleigh-Taylor，RT）不稳定性、瑞奇迈尔-莫西科夫（Richtmyer-Meshkov，RM）不稳定性和开尔文-亥姆霍兹（Kelvin-Helmholtz，KH）不稳定性。外力场中当重流体被轻流体加速或支撑时，两流体界面处的扰动随时间发展起来的一种流体力学不稳定性现象，被称为RT不稳定性^[1-2]。RT不稳定性作为一种基本的流体不稳定现象，在科学工程的许多领域都具有重要意义。例如在惯性约束聚变点火过程中，靶丸在压缩发生内爆时会产生RT不稳定性，其诱导的湍流混合会直接影响靶丸的能量增益并导致点火失败^[3]。这时发生的RT不稳定性是具有破坏作用的，我们应当尽可能减弱它。而对于超燃冲压发动机来说，RT不稳定性的存在则会使得燃料充分混合，提高其燃烧效率。RT不稳定性还广泛存在于天文领域，如旋转恒星^[4]、极光斑^[5]、行星状星云^[6]、超新星爆炸^[7]等。因此，对RT不稳定性的研究不仅具有重要的理论意义，还具有实际的应用价值。

一般的，RT不稳定性的研究方法主要有实验分析^[8-9]、理论研究^[10-13]和数值模拟^[14-21]三种。然而，自然界的RT不稳定性现象相当复杂，本身的物理系统具有较强的非线性，绝大多数实际问题的计算均超出了现有解析求解的能力范围，理论研究需要各种简化及线性假设，所得到的结果往往较为有限。实验研究则常常会存在仪器昂贵、实验危险和无法观察短时间行为等缺点。因而对RT不稳定性的研究在理论求解和实验方法上都会遇到这样或那样的不足。数值模拟作为一种当今新兴的研究工具已然成为重要的科学研究手段。相较于实验与理论研究，数值模拟可通过计算机模拟得到复杂物理系统演化过程中的精细结构和物理机制。对于RT不稳定性现象，研究人员通过数值模拟已取得了较为丰硕的成果。例如，Gallis等^[14]利用分子气体动力学的直接模拟蒙特卡罗方法，数值模拟再现了RT不稳定性混合层增长的许多定性特征，并且在线性、非线性和自相似状态下与理论和经验模型具有定量上的一致性。Wei等^[15]提出了模拟二维不可压缩 RT 不稳定性的耦合格子玻尔兹曼模型和一种修正的平衡分布函数, 并利用一个简单的标度讨论了混合RT 不稳定性的区域。Liang等^[16]利用格子玻尔兹曼方法研究了雷诺数（Re）对多模、互不相溶的RT不稳定性中的演化界面动力学和气泡/尖钉振幅的影响。Lyubimova等^[17]在相场方法的基础上，研究了限制在水平平面内的RT不稳定性中的扩散与对流演化。Scott等^[18]利用基于小波的二维可压缩直接数值模拟研究了单模RT不稳定性对涡旋动力学的影响。Li等^[19]以欧拉（Euler）模拟为基础对双界面的RT不稳定性弱非线性区域内的界面耦合效应进行了数值研究。结果表明，当下界面的阿特伍德数（At）较小时，下界面的微扰增长幅度与上界面的At呈正相关；但当下界面上的At较大时，二者呈负相关。Luo等^[20]利用二元混合流体模型研究了可压缩性对RT不稳定性的影响，发现可压缩性引起的初始密度分层起稳定作用，而膨胀-压缩效应起失稳作用。Livescu等^[21]利用直接数值模拟研究了在重力为零或反转时RT不稳定性的演化。他们发现，当加速度改变时，一些湍流量发生了显著变化。以上列举的是不同学者利用不同数值模拟方法对RT不稳定性做出的一系列研究，他们的结论丰富了我们对RT不稳定性的认知。

在RT不稳定性的数值模拟中，根据初始扰动界面结构的设置，可以分为单模RT不稳定性（流体界面被一个具有单个频率的余弦或者正弦所扰动）和多模RT不稳定性（界面被两个及以上不同频率的余弦或者正弦所扰动）。此前，多模RT不稳定性已被许多专家学者研究，并得到了一些有意义的结论。例如，Banerjee等^[22]结合不可压缩Euler方程和隐式大涡模拟方法研究初始条件对不可压缩流体的多模RT不稳定性的影响，发现RT混合的整体增长强烈依赖于初始条件。Burton等^[23]使用非线性大涡模拟方法对具有超高At的多模RT不稳定性进行研究，发现尖钉的高度和混合层生长速率受到初始密度比的强烈影响。Zhang等^[24]对多模烧蚀可压缩RT不稳定性的自相似非线性演化进行了二维和三维数值研究，发现气泡发展速度受到初始条件和烧蚀速度的影响。Liang等^[16]使用改进的相场格子玻尔兹曼方法研究了Re对不可压缩多模非混相RT不稳定性的影响，研究表明多模RT不稳定性在不同Re下表现出不同的界面动力学情形。Yilmaz等^[25]利用大涡模拟模拟了高At下的多模三维RT不稳定性，结果表明高At下RT不稳定性快速发展，尖钉的生长速率和速度增加，混合区域的不对称性增大。Hamzehloo等^[26]利用相场方法研究了At、Re、表面张力和初始扰动幅值的不同组合对不可压缩流体的多模非混相RT不稳定性的影响，发现三维多模RT不稳定性在初始阶段表现出与单模RT不稳定性相似的指数界面增长率。Ding等^[27]利用分子动力学方法研究了强加速度下的双模微观RT不稳定性，发现微观RT不稳定性表现出较弱的非线性性质，并且双模RT不稳定性在微观尺度上的模耦合行为与宏观尺度上的模耦合行为有明显差异。

根据上述研究方法的不同，可以将流体不稳定性的模拟工具分为以下三个层次：微观分子动力学模型、介观动理学模型和宏观流体力学模型。宏观流体力学模型主要是指基于Euler方程组或N-S方程组的各类模型，可以用于描述较大时空尺度的缓变行为。这类模型基于宏观连续性假设，无法描述小尺度上粒子随机运动对应的热涨落行为，缺乏描述微介观结构和快模式的能力。为了从更底层的基础上深入理解流体系统特征机制，需要利用其它数值模型手段，比如微观分子动力学方法^[28]。虽然这种方法在纳米流等问题中得到了成功应用，但由于其对计算资源的高需求，本质上局限于小尺度问题。为了解决物理精度与计算效率的矛盾问题，可以借助动理学理论构建介观层次的数值模拟方法。

近十年，基于非平衡统计物理的离散玻尔兹曼方法（discrete Boltzmann method，DBM）被成功发展应用于研究复杂流场中热力学非平衡效应（thermodynamic non-equilibrium effect，TNE）和流体力学非平衡效应（hydrodynamic non-equilibrium effect，HNE）。作为一种粗粒化的物理建模方法，DBM可以从以下两方面弥补NS模型的不足：1）适用于模拟研究具有锐利界面的流体流动，如冲击波；2）可以描述和刻画更底层的热力学非平衡信息。事实上，我们通过查普曼-恩斯柯（Chapman-Enskog，CE）多尺度分析展开，可以得到连续极限假设条件下，介观层次的分布函数演化与宏观层次的物理量演化之间的内在联系。根据CE多尺度分析，利用克努森数（Kn）中的各阶项，可以构造不同层次的描述热力学非平衡信息的DBM。在分布函数需要满足的动理学矩中存在守恒矩和非守恒矩两种类型，利用守恒矩可得到流体的宏观物理量信息，同时利用 $(f - {f^{{\rm{eq}}}})$ 的动理学矩可描述系统偏离热力学平衡态的具体信息，其中 ${f^{{\rm{eq}}}}$ 是平衡态分布函数^[29-33]。通过这些分析方法，一些以前无法提取的信息可被分层、定量地研究。

目前，DBM已被广泛用于模拟研究各类复杂流体流动，包括流体不稳定性^[34-43]、多相流^[44-45]、反应流^[46-48]、激波^[49]和爆轰^[50-56]等。这里，我们重点关注并简要介绍一下DBM研究流体不稳定性的成果。Lai等^[34]研究了可压缩性对RT不稳定性的影响，发现可压缩性在演化早期抑制了RT不稳定性的发生，在后期则促进了RT不稳定性。Li等^[35]利用DBM模拟了可压缩流体系统的多模RT不稳定性，探讨了其演化机理。Chen等^[36]利用多松弛时间DBM研究了黏性、热传导和普朗特数对RT不稳定性的影响。他们还模拟了二维RM不稳定性和RT不稳定性共存系统，讨论了RT与RM不稳定系统的相似和不同之处；研究了两种不稳定性的协作和竞争机制，并探讨了重力场 $g$ 和马赫数对非平衡的影响，发现在重力加速度的共同作用下，RT和RM不稳定共存系统中扰动的增长速度可能会增加^[37]。Gan等^[38]利用DBM研究了黏性和热传导对KH不稳定性的影响，发现黏性效应稳定了KH不稳定性，并提高了局部和全局TNE强度；而热传导效应则表现为先抑制后增强KH不稳定性。Lin等^[39]研究了KH不稳定性的动态非平衡过程，研究表明由于物理梯度和非平衡区域的共同作用，在KH不稳定性的整个演化周期中，热力学非平衡强度先增大后减小，而混合熵的增长速率则呈现先减小后增大最后减小的趋势，混合自由焓的变化趋势与混合熵的变化趋势相反。Chen等^[40]采用多松弛时间DBM模拟了耦合Rayleigh-Taylor-Kelvin-Helmholtz不稳定系统，并且引入形态边界长度和TNE强度研究了该系统的复杂流体结构和动力学过程。Ye等^[41]研究了Kn效应对可压缩RT不稳定性的影响，发现Kn增大会抑制RT不稳定性的发展，但会增强全局HNE和TNE效应。Chen等^[42]研究了比热比效应对可压缩RT不稳定性的影响，发现由于宏观物理量梯度与热力学非平衡区域之间的竞争，TNE强度先增大后减小，并随着比热比的减小而增大。Zhang等^[43]引入示踪粒子作为可压缩RT不稳定流动离散玻尔兹曼模拟的补充，并研究了混相双流体系统界面附近RT不稳定性流动的精细结构和TNE行为，同时讨论了可压缩性和黏性对RT混合的影响。DBM结果还得到了分子动力学^[28]、直接模拟蒙特卡罗方法^[57-58]等的证实和补充。

为了进一步从动理学角度深入研究可压缩流体RT不稳定性的演化规律和物理机制，本文使用Bhatnagar-Gross-Krook（BGK）DBM模拟分析多模RT系统中的HNE和TNE的动态演化过程和宏观表征，并通过分析流体系统中温度梯度与非平衡面积占比解释非平衡分量和全局非平衡效应的演化规律，进而从介观与宏观两方面了解多模RT不稳定性。

1 离散玻尔兹曼模型

为方便起见，本文统一对物理量进行了无量纲化。包含外力作用的玻尔兹曼方程形式如下：

$ {\partial _t}f + {\boldsymbol{v}} \cdot {\nabla _{\boldsymbol{r}}}f + {\boldsymbol{a}} \cdot {\nabla _{\boldsymbol{v}}}f = \varOmega $

(1)

其中， ${\boldsymbol{a}}$ 为加速度向量， ${\boldsymbol{v}}$ 为粒子速度， ${\boldsymbol{r}}$ 为空间坐标， $ f $ 是关于空间、时间、粒子速度的分布函数， $\varOmega$ 为碰撞项。由于式（1）右端碰撞项过于复杂，很难对其进行直接计算。人们提出了一系列简化碰撞模型，最常用的模型为BGK玻尔兹曼模型，其离散形式为：

$ {\partial _t}{f_i} + {{\boldsymbol{v}}_i} \cdot {\nabla _{\boldsymbol{r}}}{f_i} + {\boldsymbol{a}} \cdot {\nabla _{{{\boldsymbol{v}}_i}}}{f_i} = - \frac{1}{\tau }({f_i} - f_i^{{\rm{eq}}}) $

(2)

其中， $ i $ 表示粒子速度个数， $f_i^{{\rm{eq}}}$ 为平衡态分布函数。

DBM的构建过程是粗粒化物理建模的过程。通过CE多尺度展开，可以将上述离散玻尔兹曼方程恢复到可压缩N-S方程组：

$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho {\boldsymbol{u}}) = 0 $

(3)

$ \frac{\partial }{{\partial t}}(\rho {\boldsymbol{u}}) + \nabla \cdot (\rho {\boldsymbol{uu}} + p{\boldsymbol{I}} + {\boldsymbol{P}}) = \rho {\boldsymbol{a}} $

(4)

$ \begin{split}& \hfill \frac{\partial }{{\partial t}}\big[\rho (E + \frac{{{{\left| {\boldsymbol{u}} \right|}^2}}}{2})\big] + \nabla \cdot \big[\rho {\boldsymbol{u}}(E + T + \frac{{{{\left| {\boldsymbol{u}} \right|}^2}}}{2}) + {\boldsymbol{J}} + {\boldsymbol{P}} \cdot {\boldsymbol{u}}\big] = \rho {\boldsymbol{a}} \cdot {\boldsymbol{u}} \end{split} $

(5)

其中， $\; \rho$ 是密度， ${\boldsymbol{u}}$ 是流体速度， $ T $ 是温度， $E = (D + n)T/2$ 为系统内能， $ D $ 为空间维数， $ n $ 为额外自由度数， $ p $ 为压强， $\kappa = (D + n + 2)\mu /2$ 是热传导系数， $\mu = \tau \rho T = \tau p$ 为动力黏性系数， ${\boldsymbol{P}}$ 和 ${\boldsymbol{J}}$ 分别是与黏性应力张量和热通量有关的量，其具体形式分别为：

$ {\boldsymbol{P}} = - \mu \left[ {\nabla {\boldsymbol{u}} + {{(\nabla {\boldsymbol{u}})}^{\text{T}}} - \frac{2}{{D + n}}(\nabla \cdot {\boldsymbol{u}}){\boldsymbol{I}}} \right] $

(6)

$ {\boldsymbol{J}} = - \kappa \nabla T $

(7)

为了研究方便，我们可通过CE多尺度分析，恢复出对应的N-S层次的可压缩宏观流体力学方程组。此时，需要保证离散的平衡态分布函数 $f_i^{{\rm{eq}}}$ 满足下述7个动理学矩关系：

$ {\boldsymbol{M}}_0^{{\rm{eq}}} = \sum\limits_i {f_i^{{\rm{eq}}}} = \rho $

(8)

$ {\boldsymbol{M}}_1^{{\rm{eq}}} = \sum\limits_i {f_i^{{\rm{eq}}}} {{\boldsymbol{v}}_i} = \rho {\boldsymbol{u}} $

(9)

$ {\boldsymbol{M}}_{2,0}^{{\rm{eq}}} = \sum\limits_i {f_i^{{\rm{eq}}}\left( {{{\boldsymbol{v}}_i} \cdot {{\boldsymbol{v}}_i} + \eta _i^2} \right)} = \rho \big[(D + n)T + {\boldsymbol{u}} \cdot {\boldsymbol{u}}\big] $

(10)

$ {\boldsymbol{M}}_2^{{\rm{eq}}} = \sum\limits_i {f_i^{{\rm{eq}}}{{\boldsymbol{v}}_i}{{\boldsymbol{v}}_i}} = \rho \left( {T{\boldsymbol{I}} + {\boldsymbol{uu}}} \right) $

(11)

$ \begin{split} {\boldsymbol{M}}_{3,1}^{{\rm{eq}}} =& \sum\limits_i {f_i^{{\rm{eq}}}\left( {{{\boldsymbol{v}}_i} \cdot {{\boldsymbol{v}}_i} + \eta _i^2} \right)} {{\boldsymbol{v}}_i} \\ =& \rho {\boldsymbol{u}}\big[(D + n + 2)T + {\boldsymbol{u}} \cdot {\boldsymbol{u}}\big] \end{split}$

(12)

$ \begin{split} {\boldsymbol{M}}_3^{{\rm{eq}}} =& \sum\limits_i {f_i^{{\rm{eq}}}{{\boldsymbol{v}}_i}{{\boldsymbol{v}}_i}} {{\boldsymbol{v}}_i} = \rho \big[T({{\boldsymbol{u}}_{\alpha} }{{\boldsymbol{e}}_{\beta} }{{\boldsymbol{e}}_{\chi} }{\delta _{{\beta} {\chi} }} +\\& {{\boldsymbol{u}}_{\beta} }{{\boldsymbol{e}}_{\alpha} }{{\boldsymbol{e}}_{\chi} }{\delta _{{\alpha} {\chi} }} + {{\boldsymbol{u}}_{\chi} }{{\boldsymbol{e}}_{\beta} }{{\boldsymbol{e}}_{\alpha} }{\delta _{{\alpha} {\beta} }}) + {\boldsymbol{uuu}}\big] \end{split} $

(13)

$ \begin{split} {\boldsymbol{M}}_{4,2}^{{\rm{eq}}} = &\sum\limits_i {f_i^{{\rm{eq}}}\left( {{{\boldsymbol{v}}_i} \cdot {{\boldsymbol{v}}_i} + \eta _i^2} \right)} {{\boldsymbol{v}}_i}{{\boldsymbol{v}}_i} \\ =&\rho T\big[(D + n + 2)T + {\boldsymbol{u}} \cdot {\boldsymbol{u}}\big]{{\boldsymbol{e}}_{\alpha} }{{\boldsymbol{e}}_{\beta} }{\delta _{{\alpha} {\beta} }} +\\&\rho {\boldsymbol{u}} {\boldsymbol{u}} \big[(D + n + 4)T + {\boldsymbol{u}} \cdot {\boldsymbol{u}}\big] \end{split} $

(14)

其中， ${{\boldsymbol{e}}_{\alpha} }$ 为 $ {\alpha} $ 方向上的单位向量， $ {\delta _{{\alpha} {\beta} }} $ 是克罗内克函数。以上方程组可以写成一个矩阵形式：

$ {\boldsymbol{C}} \cdot {{\boldsymbol{f}}^{{\rm{eq}}}} = {{\boldsymbol{\hat f}}^{{\rm{eq}}}} $

(15)

其中：

$ {{\boldsymbol{f}}^{{\rm{eq}}}} = {\Big[f_1^{{\rm{eq}}},f_2^{{\rm{eq}}}, \cdots ,f_N^{{\rm{eq}}}\Big]^{\text{T}}} $

(16)

$ {{\boldsymbol{\hat f}}^{{\rm{eq}}}} = {\Big[\hat f_1^{{\rm{eq}}},\hat f_2^{{\rm{{\rm{eq}}}}}, \cdots ,\hat f_N^{{\rm{eq}}}\Big]^{\text{T}}} $

(17)

${\boldsymbol{C}}$ 是关于离散速度 ${{\boldsymbol{v}}_i}$ 和自由参数 ${\eta _i}$ 的系数矩阵， ${{\boldsymbol{\hat f}}^{{\rm{eq}}}}$ 的元素为式（8）～式（14）右边的宏观量。要保证以上方程组有解，则离散速度的个数要不少于需要用到的动理学矩数目。由于式（8）～式（14）包含16个独立的分量，因此离散速度个数至少要有16个。如图1所示，我们给出了二维十六速度模型（D2V16）。离散速度值由式（18）取得：

$ {{\boldsymbol{v}}_i} = \left\{ {\begin{array}{*{20}{l}} {{v_a}\left[ {\cos \dfrac{{(i - 1)\text {π} }}{2},\sin \dfrac{{(i - 1)\text {π} }}{2}} \right],i = 1, \cdots ,4} \\ {{v_b}\left[ {\cos \dfrac{{(2i - 1)\text {π} }}{4},\sin \dfrac{{(2i - 1)\text {π} }}{4}} \right],i = 5, \cdots ,8} \\ {{v_c}\left[ {\cos \dfrac{{(i - 9)\text {π} }}{2},\sin \dfrac{{(i - 9)\text {π} }}{2}} \right],i = 9, \cdots ,12} \\ {{v_d}\left[ {\cos \dfrac{{(2i - 9)\text {π} }}{4},\sin \dfrac{{(2i - 9)\text {π} }}{4}} \right],i = 13, \cdots ,16} \end{array}} \right. $

(18)

另外，当 $i = 1, \cdots ,4$ 时， ${\eta _i} = {\eta _0}$ ； $i = 5, \cdots ,16$ 时， ${\eta _i} = 0$ 。

在上述7个矩关系中，只有在式（8）～式（10）中的平衡态分布函数 $f_i^{{\rm{eq}}}$ 可以用分布函数 ${f_i}$ 代替，这是由于其满足质量、动量和能量守恒。但如果用 ${f_i}$ 代替任何另外四个需要的矩关系中的 $f_i^{{\rm{eq}}}$ ，则式中右侧的值与左侧的值可能会有偏差，此时这种偏差可以作为描述系统偏离平衡态的一种度量。在这个动理学思想的基础上，我们定义包含宏观信息的非平衡量 ${{\boldsymbol{\varDelta }}_{m,n}} = {{\boldsymbol{M}}_{m,n}}({f_i}) - {{\boldsymbol{M}}_{m,n}}(f_i^{{\rm{eq}}})$ 。这里， ${{\boldsymbol{M}}_{m,n}}$ 是粒子速度矩。进一步，我们用 ${\boldsymbol{v}}_i^* = {{\boldsymbol{v}}_i} - {\boldsymbol{u}}$ 代替 ${{\boldsymbol{M}}_{m,n}}$ 中的 ${{\boldsymbol{v}}_i}$ ，便得到了描述分子热涨落信息的热力学非平衡量 ${\boldsymbol{\varDelta }}_{m,n}^* = {\boldsymbol{M}}_{m,n}^*({f_i}) - {\boldsymbol{M}}_{m,n}^*(f_i^{{\rm{eq}}})$ ，包括 $\varDelta _{2{\alpha} {\beta} }^*$ 、 $\varDelta _{3{\alpha} {\beta} {\chi} }^*$ 、 $\varDelta _{(3,1){\alpha} }^*$ 、 $\varDelta _{(4,2){\alpha} {\beta} }^*$ 12个分量，其中 ${\alpha} 、{\beta} 、{\chi} 等于 x$ 或 $y$ 。

为了定性分析系统全局热力学非平衡效应，我们具体定义以下几种热力学非平衡量：

$ \left| {{\boldsymbol{\varDelta }}_2^*} \right| = \left| {\varDelta _{2xx}^*} \right| + \left| {\varDelta _{2xy}^*} \right| + \left| {\varDelta _{2yy}^*} \right| $

(19)

$ \left| {{\boldsymbol{\varDelta }}_{3,1}^*} \right| = \left| {\varDelta _{3,1x}^*} \right| + \left| {\varDelta _{3,1y}^*} \right| $

(20)

$ \left| {{\boldsymbol{\varDelta }}_{4,2}^*} \right| = \left| {\varDelta _{4,2xx}^*} \right| + \left| {\varDelta _{4,2xy}^*} \right| + \left| {\varDelta _{4,2yy}^*} \right| $

(21)

$ \left| {{\boldsymbol{\varDelta }}_3^*} \right| = \left| {\varDelta _{3xxx}^*} \right| + \left| {\varDelta _{3xxy}^*} \right| + \left| {\varDelta _{3xyy}^*} \right| + \left| {\varDelta _{3yyy}^*} \right| $

(22)

基于以上几种非平衡量，我们可进一步定义一种可以描述流体系统整体偏离平衡态程度的全局热力学非平衡量：

$ \left| {{\boldsymbol{\varDelta }}_{}^*} \right| = \left| {{\boldsymbol{\varDelta }}_2^*} \right| + \left| {{\boldsymbol{\varDelta }}_{3,1}^*} \right| + \left| {{\boldsymbol{\varDelta }}_{4,2}^*} \right| + \left| {{\boldsymbol{\varDelta }}_3^*} \right| $

(23)

对流体系统中的热力学非平衡强度进行积分并求平均就可以得到以下几种热力学非平衡量的全局平均强度：

$ \bar D = \frac{1}{{{L_x}{L_y}}}\int_0^{{L_x}} {\int_0^{{L_y}} {\left| {{{\boldsymbol{\varDelta }}^*}} \right|} } {\text{d}}x{\text{d}}y $

(24)

$ {\bar D_2} = \frac{1}{{{L_x}{L_y}}}\int_0^{{L_x}} {\int_0^{{L_y}} {\left| {{\boldsymbol{\varDelta }}_2^*} \right|} } {\text{d}}x{\text{d}}y $

(25)

$ {\bar D_{3,1}} = \frac{1}{{{L_x}{L_y}}}\int_0^{{L_x}} {\int_0^{{L_y}} {\left| {{\boldsymbol{\varDelta }}_{3,1}^*} \right|} } {\text{d}}x{\text{d}}y $

(26)

其中 $L_x $ 和 $L_y $ 分别表示计算区域的长和宽。

另外，对于离散玻尔兹曼演化方程（2）中的时间导数采用一阶精度的向前Euler格式进行处理，空间导数采用二阶精度的无波动、无自由参数的耗散有限差分格式（nonoscillatory and nonfree-parameter dissipation finite difference scheme，NND格式）^[59]。

2 数值模拟

本文在 $\varOmega = [ - d/2,d/2] \times [ - d/2,d/2]$ 的二维方形区域中进行模拟，将具有不同密度的上下两层同种物质流体置于重力加速度为常数的重力场中，并在界面处施加多模初始扰动：

$ {y_c}(x) = \sum\limits_{n = 21}^{30} {[{a_n}\cos ({k_n}x) + {b_n}\sin ({k_n}x)]} $

(27)

其中， ${k_n} = 2n\text{π} /{L_x}$ ， $ {L_x} = d $ ， ${a_n}$ 和 ${b_n}$ 是区间（0，1）中均匀分布的随机数。初始时刻流体系统满足静力学平衡条件：

$ {\partial _y}p(y) = - g\rho (y) $

(28)

其中 $g$ 为重力加速度。流体系统的初始条件为：

$ \left\{\begin{array}{l}\begin{array}{l}T(y) = {T}_{u}\text{，}\rho (y) = \dfrac{{p}_{0}}{{T}_{u}}\mathrm{exp}\bigg\{\dfrac{g}{{T}_{u}}\Big[2d-{y}_{c}(x)\Big]\bigg\}\text{，}\\ \;\;\;\;\quad\qquad\qquad\qquad\qquad\qquad\qquad y > {y}_{c}(x)\end{array} \\ \begin{array}{l}T(y) = {T}_{b}\text{，}\rho (y) = \dfrac{{p}_{0}}{{T}_{b}}\mathrm{exp}\bigg\{\dfrac{g}{{T}_{u}}\Big[2d-{y}_{c}(x)\Big]-\\\qquad \dfrac{g}{{T}_{u}}\Big[y-{y}_{c}(x)\Big]\bigg\}\text{，}\\\;\;\;\;\quad\qquad\qquad\qquad\qquad\qquad\qquad y > {y}_{c}(x)\end{array} \end{array}\right. $

(29)

其中， ${y_c}(x)$ 表示初始扰动界面位置， ${p_0}$ 表示流体系统顶部的初始压强， ${T_u}$ 和 ${T_b}$ 分别表示上部流体和下部流体的初始温度。另外，在现实物理系统中，流体界面附近的物理量变化并非强间断，而是光滑过渡，为此引入tanh函数描述过渡层变化，即温度场设置为：

$ T(y) = \frac{{{T_b} + {T_u}}}{2} - \frac{{{T_b} - {T_u}}}{2} \cdot \tanh \left[ {\frac{{y - {y_c}(x)}}{W}} \right] $

(30)

其中 $ W $ 表示流体界面处过渡层宽度。

需要说明的是，本文模拟的是可压缩流体系统，即由于重力的作用，系统中的压强从上到下逐层增大。因此，当设定上下层部分流体温度分别为定值 ${T_u}$ 和 ${T_b}$ 时，那么各部分流体密度将随着压强的变化而变化（见式（29）），并且存在如下关系： $ p = \rho RT $ ，其中 $ R $ 为阿伏加德罗常数。同时，为使流体界面处的压强保持一致，根据理想气体的物态方程可知：

$ {\rho _u}{T_u} = {\rho _b}{T_b} $

(31)

其中 $\;{\rho _u}$ 和 $\;{\rho _b}$ 分别是界面两侧上下流体的密度。另外，在处理边界条件时， $x$ 方向采用周期边界条件， $y$ 方向采用固壁边界条件。

数值模拟中，在满足网格无关性条件下，将模拟区域划分为 ${N_x} \times {N_y}=$ $512 \times 512$ 的均匀网格，空间步长为 $\Delta x = \Delta y$ $ = 5 \times {10^{ - 4}}$ ，时间步长为 $\Delta t = 2 \times {10^{ - 5}}$ ，并且 $d = \Delta x \times {N_x}$ 。系统上下部分流体的温度分别为 ${T_u} = 1.0$ 和 ${T_b} = 4.0$ ，模拟区域顶部的初始压强为 ${p_0} = 1.0$ 。离散速度参数设置为（ ${v_a} $ ， ${v_b} $ ， ${v_c} $ ， ${v_d}$ ）=（0.9，0.9，4.2，4.2），其他参数如下： ${\eta _0} = 15.0$ 、 ${a_x} = 0.0$ 、 ${a_y}=\;$ $ - g = - 1.0$ 、 $\tau = 4 \times {10^{ - 5}}$ 。

为了直观感受流体系统随时间的演化情况，图2给出了系统中不同时刻的温度云图。可以看出，初始时刻界面处存在不同振幅不同波长的微小扰动，并且由于在界面附近使用了tanh函数，因此物理场的空间梯度变化较为光滑。演化开始后，伴随着系统耗散和热传导的发生，小波长随机扰动界面变得光滑，过渡层变宽，宏观物理量梯度从尖锐转向光滑。大约 $t = 1.5$ 时，在重力的作用下，随机扰动界面不断增大，流体界面发生了较大的形变，并逐渐形成典型的尖钉与气泡结构。在后期（约 $t > 2.5$ ），由于界面切向速度的作用，尖钉（气泡）两侧形成了向上（下）翻转的“涡”结构。随后，尖钉与气泡不断延伸，相互耦合、相互渗透，系统进入混合阶段。

本文的模拟结果与经典多模RT不稳定性的主要特征基本一致^[60-61]，但同时也存在些许不同，造成这种差异的原因有以下几点：

1）随机扰动模态不同。在程序中，通过随机函数在物质界面产生随机扰动模态，不同软件平台的随机函数产生的结果往往不一致。

2）流体可压性的影响。与不可压模型不同，本文使用的可压缩DBM可以用于准确预测和研究RT不稳定性的可压缩效应。

3）与Euler、N-S等传统流体力学模型不同，本文使用的DBM不仅包含了黏性和热传导的影响，还同时包含了其他重要的非平衡效应。

由于多模随机初始扰动中扰动种子的不确定性，我们模拟了10组不同随机扰动情况下的RT不稳定性，并根据其统计特性进行分析。首先，我们分析了RT不稳定系统中系统温度梯度的变化趋势。图3（a）为 $x$ 方向全局平均温度梯度 $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 随时间演化图。总体来看， $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 的发展趋势为：先下降、后上升、再下降。具体过程如下：

第一阶段，即初始阶段（大约 $0 < t < 0.1$ ）， $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 快速下降。物理原因如下：在多模RT不稳定性的初始界面中，小波长随机扰动很快被耗散掉，扰动模式快速减少，物理场变得光滑，从而宏观物理量梯度降低。

第二阶段， $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 先缓慢增长（大约 $0.1 < t < 0.9$ ）、后快速增长（大约 $0.9 < t < 3.5$ ），大致呈现出指数增长趋势。物理原因如下：当 $0.1 < t < 0.9$ 时，扰动的波长小、频率大，使得扰动界面发展缓慢，热传导起主要作用，此时过渡层变宽，物理量梯度变化缓慢。当 $0.9 < t < 3.5$ 时，界面快速拉长，重流体向下发展形成“尖钉”，轻流体不断上升形成“气泡”。并且由于界面剪切的作用，KH不稳定性产生作用，加速了RT不稳定性的发展，尖钉气泡顶部卷起形成了典型的“蘑菇”状的涡结构， $x$ 方向上的温度等温层结构逐渐变得复杂。同时，界面附近的物理量梯度下降，此时界面的拉伸作用占主导地位。因此，全局平均温度梯度快速增长。

第三阶段（ $t > 3.5$ ）： $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 整体呈现下降趋势。此阶段，扰动界面继续拉长，但系统逐渐趋于平衡，两流体进一步融合，耗散作用使得扰动界面附近的物理量梯度下降、小结构逐渐消失。由于在这个过程中后者起主导作用，因此物理量梯度下降。

图3（b）展示了 $y$ 方向全局平均温度梯度 $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ 的演化情形。如图所示， $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ 的发展趋势为：首先保持定值，然后增长，最后下降。具体分析如下：初始阶段（大约 $0 < t < 2.0$ ）， $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ 维持在一个定值附近。这是由于在这段时间里，界面扰动发展缓慢，系统中界面两侧的剪切作用较小，界面两侧并未明显出现“涡”结构，此时 $y$ 方向的温度变化是均匀的，因此 $y$ 方向的温度梯度保持在一个定值附近。增长阶段（大约 $2.0 < t < 3.5$ ），流体系统继续演化，出现了尖钉与气泡结构，尖钉（气泡）向下（上）延伸，流体界面不断拉长，混合区域面积不断增加，并且界面两侧的剪切作用引发了KH不稳定性的产生，导致尖钉与气泡两侧产生了涡结构，系统中的流场结构变得较为复杂， $y$ 方向的温度变化不再均匀，全局平均温度梯度增加。下降阶段（大约 $t > 3.5$ ），流体界面继续延展，此时混合程度进一步增大，在耗散作用下混合区域界面变得光滑，物理量梯度变小，小结构也随之慢慢消失，系统中全局平均温度梯度逐渐减小。值得一提的是， $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ 的下降曲线并不是一种逐渐下降的趋势，而是一种振荡交替下降。这是由于系统中同时存在界面拉伸与耗散的作用，这两种作用互相竞争，从而使之振荡下降。图3（c）展示了全局平均温度梯度 $\mathop {\left| {\nabla T} \right|}\limits^{\_\_\_\_\_} $ 的变化曲线图，事实上， $\mathop {\left| {\nabla T} \right|}\limits^{\_\_\_\_\_} $ 由 $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 和 $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ 共同决定，因此 $\mathop {\left| {\nabla T} \right|}\limits^{\_\_\_\_\_} $ 的演化趋势可以根据 $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ 与 $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ 的演化规律得到。

图4给出了与热传导相关的非平衡量的全局平均强度演化图，图4（a、b、c）分别对应 ${\bar D_{3,1x}}$ 、 ${\bar D_{3,1y}}$ 和 ${\bar D_{3,1}}$ 的演化图。首先，由式（7）可知， ${\bar D_{3,1x}}$ 、 ${\bar D_{3,1y}}$ 和 ${\bar D_{3,1}}$ 的变化与热传导系数和温度梯度相关。从图中可以看出，在10组不同实验中， ${\bar D_{3,1x}}$ （ ${\bar D_{3,1y}}$ ）的变化趋势与 $\mathop {\left| {{\nabla _x}T} \right|}\limits^{\_\_\_\_\_\_} $ （ $\mathop {\left| {{\nabla _y}T} \right|}\limits^{\_\_\_\_\_\_} $ ）的变化趋势基本相同。在其演化过程中，主要有以下几个物理机制相互作用、相互竞争：1）温度梯度的增加（减小）使得热交换增强（减弱）；2）演化过程中流体接触面积的增加（减小）导致热交换增加（减弱）。对于机制2）而言，造成流体接触面积增大（减小）的原因有：1）尖钉（气泡）下降（上升）使得流体界面被拉伸延长，接触面积增加；2）KH不稳定性的产生使得流体接触面增加，混合区的形态结构变得更加复杂；3）随着尖钉（气泡）接触到下（上）壁面，流体界面的增加空间受限，当流体的混合达到一定程度后，混合区面积不再持续增加。总体来说，在系统演化过程中，流体的流动使得温度接触界面具有拉长的趋势，而热传导具有降低温度梯度、消除温度梯度的趋势。需要说明的是，尖钉与气泡先后到达壁面，接触区域仍在不断增长，但温度梯度变化较为复杂，所以此时 ${\bar D_{3,1x}}$ 的起伏较大。最后，流体系统演化逐渐趋于饱和，系统向平衡态发展，温度梯度与接触面积都呈下降趋势，因此 ${\bar D_{3,1x}}$ 不断减小。我们用同样的思路分析 ${\bar D_{3,1y}}$ 的变化趋势：前期，温度梯度的影响对 ${\bar D_{3,1y}}$ 的发展占主导作用，此时 ${\bar D_{3,1y}}$ 保持一个稳定的值；中期，随着温度梯度与接触面积的不断增加， ${\bar D_{3,1x}}$ 随之上升；后期，接触面积不再增加，并且温度梯度逐渐减小， ${\bar D_{3,1x}}$ 表现为下降的趋势。对于 ${\bar D_{3,1}}$ 而言， ${\bar D_{3,1x}}$ 和 ${\bar D_{3,1y}}$ 的共同作用使其产生如图4（c）所示的起伏变化趋势。

进一步可以发现，全局平均TNE强度 $\bar D$ 与宏观物理量梯度（包括温度、密度和速度梯度）和非平衡区域面积密切相关。数值研究表明，在RT不稳定性发展的早期阶段，局部物理量梯度减小，局部TNE强度减小。随着RT不稳定性的发展，非平衡区域面积增加，使得全局平均TNE强度增强，即宏观物理梯度效应和非平衡区域效应是相互竞争的。进一步，我们给出了非平衡区域面积占比（ ${S _r}$ ）的演化图，见图5。为了直观了解非平衡强度的演化，图6展示了系统局部非平衡强度在不同时刻的演化云图（从蓝色到红色表示非平衡强度增加）。在前期（大约 $0 < t < 1.2$ ），非平衡区域先增加后减小。增加的原因是：接触界面过渡层变宽，非平衡区域增加；减小的原因是：界面处的非平衡强度会随着扩散作用逐渐减小，这导致了大于阈值的非平衡区域占比的减小。之后（大约 $t > 1.2$ ），非平衡区域先增加、后减小，其原因是：流体开始形成尖钉与气泡结构，并且随着尖钉（气泡）下降（上升），非平衡区域不断伸展， ${S _r}$ 不断增大，并在 $t = 3.5$ 左右到达峰值。随着尖钉气泡到达上下壁面，非平衡区域增加空间受限，流体的融合也更加充分，此时系统中的物理量梯度变得光滑，非平衡强度逐渐减小，使得 ${S _r}$ 下降，最后保持在一个定值附近。同时，系统中出现的涡结构使得流场更加复杂，并在后期呈现出不规则的振荡。

图7给出了系统全局TNE强度随时间演化图。 $\bar D$ 的演化情况与宏观物理量梯度（包括温度、密度和速度梯度）和非平衡区域面积均密切相关。从图中可以看到， $\bar D$ 在演化早期（大约 $t < 1.0$ ）首先急剧下降，然后是较长时期的缓慢变化，数值小于0.01。随着流体界面逐渐拉长，系统中出现了越来越多的小结构，非平衡区域的增大也使得 $\bar D$ 呈现出增长的趋势，在大概 $t = 3.5$ 时到达峰值。随着系统物理量梯度的减弱，非平衡区域也减小，此时系统逐渐向平衡态发展， $\bar D$ 逐渐减小，直至趋于稳定状态。

3 结　论

本文利用DBM研究了可压缩流体中的多模初始扰动的RT不稳定性。首先，分析了与热通量相关的热力学非平衡分量的演化趋势，有两个物理机制起主要作用：一是温度梯度的增大使非平衡强度增大；二是两种流体接触面积的增大促进热交换，从而使非平衡区域增大。两个因素相互竞争共同影响全局非平衡强度的发展趋势。其次，分析了非平衡面积占比的变化趋势。界面热扩散作用使得界面上的局部非平衡强度有减小的趋势，同时界面的拉伸增加非平衡的区域。两种竞争机制使得非平衡区域呈现出复杂的变化趋势。最后，分析了全局平均TNE强度，全局平均TNE强度先增后减最后趋于稳定。此过程也存在两个竞争机制：非平衡区域的面积增大（减小）会增强（减弱）非平衡强度，物质界面物理梯度的增大（减小）对全局平均热力学非平衡强度有相同的影响，二者相互竞争，使其呈现出先减、后增、再减的趋势。这些现象和规律有助于我们更好地理解和分析可压缩RT不稳定性背后的物理机制和机理，为今后相关的流体不稳定性研究提供介观尺度的物理参考。