多模态瑞利-泰勒不稳定性的离散玻尔兹曼数值研究

陈 璐，赖惠林，*，林传栋，李德梅

(1. 福建师范大学 数学与统计学院，福州 350117；2. 福建师范大学 福建省分析数学及应用重点实验室，福州 350117；3. 福建师范大学 福建省应用数学研究中心，福州 350117；4. 中山大学 中法核工程与技术学院，珠海 519082)

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−feq)的动理学矩可描述系统偏离热力学平衡态的具体信息，其中feq是平衡态分布函数[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 离散玻尔兹曼模型

图1 离散速度示意图Fig. 1 Schematic diagram of discrete velocities

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

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

2 数值模拟

图2 t 分别为0 、0.5、1.0、1.5、2.0、2.5、3.0和 3.5时的温度轮廓图Fig. 2 Contours of the temperature at t =0,0.5,1.0,1.5,2.0,2.5,3.0 and 3 .5 respectively

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

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

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

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

由于多模随机初始扰动中扰动种子的不确定性，我们模拟了10 组不同随机扰动情况下的RT 不稳定性，并根据其统计特性进行分析。首先，我们分析了RT 不稳定系统中系统温度梯度的变化趋势。图3（a）

图3 在10 组不同随机扰动情况下全局平均温度梯度随时间的演化图Fig. 3 Time evolution of the global average temperature gradients under ten initial conditions with different random perturbations

图4 不同随机扰动下与热传导相关的非平衡量的全局平均强度演化图Fig. 4 Time evolution of the global non-organized energy flux with different initial random perturbations

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

图5 不同随机扰动下非平衡区域面积占比随时间演化图Fig. 5 Time evolution of the proportion of the non-equilibrium region with different initial random perturbations

图6 t 分别为0 .02、0.5、1.0、1.5、2.0、2.5、3.0 和 3.5时的局部非平衡强度轮廓图Fig. 6 Contours of the local non-equilibrium effects at t =0.02, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 and 3 .5 respectively

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

图7 不同随机扰动下全局平均TNE 强度随时间演化图Fig. 7 Evolution of the global average TNE intensity with different initial random perturbations

3 结 论

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