Richtmyer-Meshkov（RM）不稳定性^[1-2]指的是激波与界面作用导致界面上扰动增长并失稳的现象。界面的形成可以是由于两侧的介质具有不同的密度^[3-5]、不同的温度^[6-7]或者不同的其他热力学参数（如比热比）^[8]。激波与单模小扰动界面相互作用问题是最简单、最基础的RM不稳定性问题，也是RM不稳定性研究领域最为关心的问题之一，因为复杂的扰动可以通过傅里叶分析分解为不同波长的单模波的叠加，因此单模RM不稳定性的研究可以为更复杂界面的RM不稳定性研究^[9-10]提供帮助。RM不稳定性在许多领域都有广泛的应用，包括惯性约束核聚变（inertial confinement fusion, ICF）^[11]、超新星爆发^[12]、超燃冲压发动机^[13]和气云爆炸^[14]等领域。前人通过理论、实验和数值模拟的方式不断推进对单模RM不稳定性的研究，下文将从这3个方面简要介绍相关进展。

1960年，Richtmyer^[1]将RM不稳定性看作冲击加载的Reyleigh-Taylor (RT)不稳定性^[15-16]，理论上得到了界面演化的脉冲模型，该模型于1969年在Meshkov^[2]的实验中被验证。后来，Meyer等^[17]以及Vandenboomgaerde等^[18]在脉冲模型^[1]的基础上，结合实验和数值模拟结果提出了修正的线性模型。该线性模型^[17]虽然能预测扰动增长的早期行为，但难以描述大振幅下的非线性增长规律。后来，学者们提出了大量的非线性模型，如通过Padé近似拓展后的Zhang-Sohn (ZS)模型^[19]、拟合线性和渐近阶段的Sadot模型^[20]、适用于大密度比的Dimonte-Ramaprabhu (DR)模型^[21]和能归一化气泡与尖钉增长率的Zhang-Guo (ZG)模型^[22]等。这些非线性模型进一步提升了模型对实验结果的预测能力。但是，目前还没有一种模型适用于所有的工况^[23-24]。

对单模RM不稳定性的实验研究主要通过激波管开展，实验技术中主要涉及界面的形成和流场的测量。单模界面的形成方法有高分子膜片技术^[25-26]和通过震荡形成驻波的无膜技术^[27-28]等，近年来，Guo等^[10]、罗喜胜等^[29]、Zhai等^[30-31]、Li等^[32]、王宏辉等^[33]、Liu等^[34]和马迪等^[35]提出了利用线约束的可控界面肥皂膜技术。流场测量主要通过高速摄影^[36]、粒子图像测速（particle image velocimetry, PIV）技术^[37-38]和PIV-PLIF (planar laser induced fluorescence)同步测试技术^[39]等进行。利用这些技术，前人获得了丰富的实验结果，例如：Niederhaus等^[40]和Collins等^[41]首次采用无膜技术给出了单模界面演化过程，获得了排除膜片影响的实验结果；近期Liu等^[34]通过线约束肥皂膜的新技术研究了二维小振幅air/SF₆单模界面演化，由于其界面形状精确可控，并且肥皂膜对结果影响较小，因此这些结果可以为理论和数值结果提供基础的参考数据。但是目前界面形成方法无法做到对结果无任何影响，而流场显示技术只能获得界面演化的宏观结构，难以在一次实验中同时获得流动的速度场和压力场等精确的分布信息。

近年来，随着高精度计算方法和计算机计算能力的发展，数值模拟越来越被广泛应用在单模RM不稳定性的研究中。数值模拟可以同实验数据和理论模型相互印证^[42]，共同推进对单模RM不稳定性的认识。例如，Bai等^[43]采用基于切割网格的MuSiC+和CCGF程序对广泛使用的线性和非线性模型进行了对比研究，对各种模型的准确性从数值模拟角度进行了评估。数值模拟能精细地给出流场信息，并可采用控制变量的方式揭示不同因素对RM不稳定性的影响规律，这些因素包含黏性^{[42, 44-45]}、马赫数^[46]、表面张力^[47]、湍流^{[44, 48-49]}、化学反应^[50-51]、真实气体模型^[52]、磁场^[53]和强度^[54-55]等。而本文中关注的高马赫数（或者说强激波）下的单模RM不稳定性现象在实验上较难开展研究，目前主要通过数值模拟进行研究。例如，Sun等^[56]通过分子动力学研究了马赫数高达10的微观尺度的RM不稳定性问题，考虑了黏性的影响。Dell等^[57]通过SPH (smoothed particle hydrodynamics)模拟了冻结流中马赫数高达10的单模扰动增长情况，发现初始扰动增长率是初始扰动的非单调函数，并且马赫数会影响初始扰动增长率的峰值。Rikanati等^[58]采用界面追踪方法对铍泡沫介质中激波马赫数高达15.3的实验进行了数值模拟，确定了初始大振幅和高激波马赫数都会导致初始增长率低于线性理论预测。Samtaney等^[59]在热化学平衡流中对激波马赫数为18的激波冲击氢气/氦气界面的RM不稳定性现象进行了数值模拟。Zanotti等^[60]采用高阶自适应程序研究了轻/重流体和重/轻流体中激波马赫数对扰动增长的影响，最大激波马赫为32.2，其研究发现，激波马赫数的增大会导致扰动增长的减缓，并且在激波马赫数2以内存在一个增长最快的激波马赫数。本文中主要考虑气体介质，当激波很强时，双原子气体介质中会出现热化学反应，产生高温非平衡现象，而考虑热化学非平衡效应的数值模拟工作还较缺乏，本文中将研究热化学非平衡效应下的单模界面RM不稳定性现象，并对比冻结流气体、热非平衡气体和热化学非平衡气体3种气体模式的不同，在本文中这些气体模式的定义如下。

冻结气体（frozen gas, FG）：介质为经典的量热完全气体，不考虑分子振动能激发，气体的比热比保持恒定为1.4，同时不考虑化学反应，气体组分保持为初始组分比例，即氮气摩尔分数为0.79，氧气摩尔分数为0.21。

热非平衡气体（thermal nonequilibrium gas, TNG）：气体的比热比随温度变化，考虑分子振动能的激发，但不考虑化学反应，气体组分同样保持为初始组分比例。该气体模式不是本文的重点，主要用来同冻结流和热化学非平衡流进行对比，以揭示热效应和化学效应的影响程度。

热化学非平衡气体（thermo-chemical nonequilibrium gas, TCNG）：气体的比热比不是常数，考虑分子振动能的激发，同时考虑化学反应，但没有电离发生。

本文中，将研究高马赫数激波作用下单模界面RM不稳定性问题，涵盖轻/重界面和重/轻界面2种情况，涉及的激波马赫数范围分别为6～9和8～11。本文中不考虑后期的湍流混合阶段。下文将首先介绍数值算法和程序，然后给出网格收敛性测试。结果讨论部分首先给出以上3种气体模式下波系和界面的演化过程，然后对比数值模拟与不同理论模型得到的扰动增长规律，详细讨论不同气体模式、激波马赫数及初始扰动尺度对扰动振幅及其增长率的影响，最后分析涡量和环量的演化特征。

1.   数值方法

1.1   控制方程和热化学模型

本文中仅考虑无黏气体，不考虑黏性、传热和扩散，控制方程在二维Euler方程的基础上增加了组分方程和振动能方程：

式中：U为守恒型变量，t为时间，F_j为j (j=1, 2)方向对流通量，S为化学源项和振动能源项。U、F_j和S的具体表达式分别为：

式中：ρ_s为组分s(s = 1～5)的密度，s依次代表N₂、O₂、NO、N和O；ρ为总密度；u₁为x方向速度；u₂为y方向速度；E为总能，表达式为$E =\rho e +\text{ }0.5\rho ({u}_{1}^{2}+{u}_{2}^{2})$；e为比内能；E_v为振动内能；ω_s为由化学反应造成的单位体积质量生产率；ω_v为单位体积的振动能生成率。e的具体表达式为：

式中：等号右边依次为平动内能、转动内能、振动内能和等效零点能；X_s为组分s的质量分数；R_s为组分s的气体常数，其值为R_s=R_u/W_s，R_u为通用气体常数，W_s为组分s的相对原子量；T_tr为转动-平动温度；${e_{{\text{v}},s}}$为组分s的振动能；$h_s^{(0)}$为生成焓，每种组分的生成焓可以在文献[61]中找到。需要注意的是，转动内能和振动内能只包含双原子分子。组分s的振动能表达式为：

式中：h为普朗克常数，k为波尔兹曼常数，v_s为振动特征频率，T_v为振动温度，θ_v,s为组分s的特征振动温度。总的振动内能为${E_{\text{v}}} =\displaystyle \sum\limits_{s = 1}^5 {{\rho _s}{e_{{\text{v}},s}}}$，振动能源项为：

式中：松弛时间τ_s在T_tr低于8000 K时通过Millikan等^[62]提出的方法进行计算，高于8000 K时采用Park^[63]提出的修正公式进行计算；化学反应引起的振动能变化项系数$\tilde c$在本文中取1。在热非平衡气体中，无化学反应源项，而在冻结气体中，无振动能和化学反应源项。

本文中采用的化学反应模型为十七步反应模型^[61]：

式中：Q_s为N₂、O₂、NO、N和O中的任意一种。反应源项：

式中：n_r为涉及的反应个数；α_s,r和β_s,r为第r个反应的反应物和生成物的化学计量系数；K_f,r为正反应速率；K_b,r为逆反应速率。K_f,r由Arrhenius公式得到：

式中：C_0r、C_1r和C_2r为Arrhenius常数，不同化学反应值不一样，本文中采用Park的值^[64]；对于离解反应，${T_{\text{f}}} = \sqrt {{T_{\text{v}}}{T_{{\text{tr}}}}}$，对于其他反应，${T_{\text{f}}} = {T_{{\text{tr}}}}$。逆反应速率K_b,r由平衡常数^[64]得到：

平衡常数$K_r^{{\text{eq}}}\left( {{T_{\text{b}}}} \right)$为：

式中：${\textit{z}} = 10\;000/{T_{\text{b}}}$，${T_{\text{b}}} = {T_{{\text{tr}}}}$。

为显示界面的运动和变形，本文中采用颜色函数f进行标记，初始时f在界面一侧为0，在另一侧为1。其控制方程为对流方程：

1.2   程序简介

本文中使用VAS2D程序进行数值模拟，该程序为基于二阶有限体积方法的二维程序，采用自适应非结构网格。在数值求解过程中，采用了分裂方法，将对流项和源项的求解分开。对流项的求解采用了MUSCL-HANCOCK方法进行时空二阶重构，使用了HLL (Harten-Lax-van Leer)格式求解数值通量^[65-66]，源项的求解采用了具有A-stability单步隐式时间格式^[67]。在界面颜色函数f的求解采用了类似于Adapt-HLLC (Harten-Lax-van Leer contact)格式^[68]的Adapt-HLL格式。VAS2D已在不稳定性问题中得到广泛应用^[30-31]，Li等^[32]对该程序模拟热化学非平衡气体的适用性进行了大量验证，限于篇幅的原因，本文中不再给出。

2.   问题描述

本文中将研究右行激波与单模界面的相互作用过程，初始界面两侧流体分轻/重和重/轻2种分布方式（轻/重界面表示初始激波从左侧轻流体进入右侧重流体，重/轻则相反），轻/重界面流场初始分布如图1所示。单模界面两侧的温度不同，波前轻流体的温度为600 K，重流体的温度为293 K。在激波到达界面之前，界面两侧处于压力平衡状态，设置初始压力p₀为1 kPa。初始时设置较低压力，一方面考虑到在高马赫数激波管实验中多采用抽真空的方式实现高马赫数的流动条件；另一方面可以防止较高马赫数的初始激波作用下出现温度过高产生电离从而超过本文程序的计算范围（一般保持在10000 K以下）。本文中研究的激波马赫数M_s和界面初始振幅a₀范围见表1，单模界面波长λ = 20a₀。计算域为[0, 200λ] × [0, 0.5λ]，初始界面中心位置处于x = 100λ处，数值模拟只计算半个波长。初始激波从左向右运动，其初始位置距离界面上游为0.025λ。上下边界为无黏固壁，左右边界为开口边界。在本文的研究中，主要关注扰动的早中期增长过程，并未涉及后期的湍流状态。在扰动早中期演化过程中，正弦界面具有对称性，通过测试表明，在本文研究的最长时间段内，多个波长的结果并未丧失对称性，因此选择半个波长计算是合理的且极大减小了计算量。在无黏流动中，固体壁面边界条件在数值上等同于对称边界。在无黏固壁处，法向速度采用无穿透边界条件，切向速度、温度、压力沿着壁面法向均为零梯度边界条件。由于计算量的限制和非结构数据传递的困难，本文中并未进一步考察周期边界条件和多个波长的影响。计算采用自适应非结构网格进行。下文中结果以初始激波接触上游界面的时刻为零时刻。由于考察时间段内，界面将移动较远的距离，因此采用匀速移动的惯性坐标系，坐标系的移动速度为透射激波速度。由于计算域足够大，在考察的时间段，不会有从计算域左、右边侧的反射波影响到界面扰动增长区域。为进行对比研究，在冻结流、热非平衡流和热化学非平衡流中，M_s取值一致，为此需要对非平衡流中初始激波波后区域进行特定的初始化。对于热化学非平衡流，初始激波波后的压力、密度、温度和组元分布可以通过一维定常控制方程^[32]得到，并插值到二维初始网格上，热非平衡流的初始化方法同热化学非平衡流类似。以轻/重界面M_s = 8的热化学非平衡流为例，激波波后压力和密度分布情况如图2所示，x = 0处为激波面的位置。由于非平衡气体中从非平衡到平衡状态的弛豫过程具有有限速率，使得激波波后存在弛豫距离，也即该问题本身具有內秉尺度。因此，本文结果主要采用有量纲形式表述，以便对比各种不同气体模式的结果。

图  1  激波冲击轻/重单模界面初始示意图

Figure  1.  Schematic of shock wave impact on a light/heavy single-mode interface

下载: 全尺寸图片 幻灯片

表  1  初始条件

Table  1.  Initial conditions

界面类型	气体模式	Ms	a0/mm		界面类型	气体模式	Ms	a0/mm
轻/重界面	冻结流	6, 7, 8, 9	0.75, 7.5, 75		重/轻界面	冻结流	8, 9, 10, 11	0.75, 7.5, 75
	热非平衡流					热非平衡流
	热化学非平衡流					热化学非平衡流

下载: 导出CSV 

| 显示表格

图  2  正激波后的压力和密度分布

Figure  2.  Pressure and density distributions after the normal shock wave

下载: 全尺寸图片 幻灯片

网格收敛性验证如图3所示，对冻结流中M_s = 11和a₀ = 7.5 mm条件下的重/轻界面数值模拟分别采用了4层、5层和6层加密，分别对应一个波长320、640和1280个网格。由于固定了波长对应的网格数，当波长变化时，最小网格尺寸等比例缩放。可以看到，5层加密网格已经收敛，因此后文中采用5层加密进行计算。如无特殊说明，本文中后续分析的扰动参数均为a₀ = 7.5 mm，λ = 150 mm，轻/重界面M_s = 8，重/轻界面M_s = 9。

图  3  不同网格尺度下重/轻界面的振幅增长率

Figure  3.  Amplitude growth rates at different grid resolutions for the heavy/light interface

下载: 全尺寸图片 幻灯片

3.   结果分析

3.1   流动基本特征

激波马赫数为8的激波冲击轻/重界面时的界面和压力场演化如图4所示，可以看到，由于初始激波从低波阻抗介质进入高波阻抗介质，在与界面作用后形成透射激波和反射激波（见图4(a)）。初始界面上的小扰动在激波经过后逐渐增强，进入非线性阶段后轻流体进入重流体形成气泡，而重流体进入轻流体形成尖钉结构，尖钉结构逐渐出现典型的蘑菇状卷起（见图4(d)）。对比3种气体模式下的流场演化表明，透射激波、反射激波和界面的速度有较明显的差异，三者的运动速度均是热化学非平衡流最快，热非平衡流次之，冻结流最慢。透射激波和反射激波在与界面相互作用早期为非平直形态，随着时间的演化逐渐趋于平直。非平直的激波与上下壁面作用会形成横向传播的反射波，这些反射波又会相互作用形成棋盘状的压力分布形态（见图4(c)），该横波在前人的研究^[69-70]中也有发现。

图  4  轻/重界面时的界面和压力场演化

Figure  4.  Evolution of interface and pressure field for the light/heavy interface

Top: FG; middle: TNG; bottom: TCNG.

下载: 全尺寸图片 幻灯片

激波马赫数为9的激波冲击重/轻界面的界面和压力场演化如图5所示，由于初始激波从高波阻抗介质进入低波阻抗介质，与界面作用后形成透射激波和反射膨胀波（见图5(a)）。单模界面振幅在经历压缩后反向增长，从线性阶段逐渐进入非线性阶段，同样形成气泡和尖钉结构，尖钉结构最后形成典型的蘑菇状卷起（见图5(d)）。对比3种气体模式的流场发展可以发现，透射激波、反射稀疏波和界面的速度有较明显的差异。透射激波在冻结流中最快，热非平衡流次之，热化学非平衡流最慢。反射稀疏波及界面向右移动速度均为热化学非平衡流最快，热非平衡流次之，冻结流最慢。在重/轻界面情况下也存在由于横向波相互作用而形成的棋盘状压力分布形态（见图5(c)）。

图  5  重/轻界面时的界面和压力场演化

Figure  5.  Evolution of interface and pressure field for the heavy/light interface

Top: FG; middle: TNG; bottom: TCNG.

下载: 全尺寸图片 幻灯片

图6给出了轻/重界面和重/轻界面时透射激波、界面和反射波（激波或者稀疏波）的数值模拟结果与一维平衡流理论的时空位置的对比。理论解计算方法可参考一般气体模型中的黎曼问题计算方法^[71]。理论计算时，采用初始激波波后无穷远处的流动状态和下游界面后方气体无穷远处的气体状态为黎曼问题两侧的初始状态，得到了冻结气体、热平衡气体（thermal equilibirum gas, TEG）和热化学平衡气体（thermo-chemical equilibirum gas, TCEG）的理论解。对比数值结果和理论结果发现，数值结果中透射波、界面和反射波的相对快慢与平衡流理论相同，冻结流的结果与理论吻合很好。由于非平衡效应的存在，热非平衡气体和热化学非平衡气体的结果同平衡流理论解相差相对冻结流较大，体现了非平衡效应的影响。

图  6  激波与轻/重和重/轻界面作用界面及波系随时间的变化

Figure  6.  Positions of the interfaces and wave systems of the shock wave-light/heavy and -heavy/light interface interaction at different times

The positions of the interface, transmitted and reflected shock waves are indicated in red, black and blue, respectively.

下载: 全尺寸图片 幻灯片

3.2   扰动增长与增长率

本节中，首先简要回顾典型的单模界面RM不稳定性的扰动增长理论模型，并在高马赫数下对比数值模拟结果和多种理论模型，随后分析了不同气体模式、初始扰动尺度和激波马赫数对扰动增长的影响规律。

3.2.1   理论模型与数值结果的对比

对于冻结流单模RM不稳定性，目前已发表大量理论模型^{[17-23, 72]}。这些模型在低马赫数情况下，得到了较多的实验和数值模拟结果支持^{[28, 33-34]}，但在高马赫数情况下，这些模型的预测能力还没有得到详细的检验。下文中将对比高马赫数情况下这些理论模型和本文的数值模拟结果，为理论模型的评估提供数值支撑。

单模界面RM不稳定性的线性理论最早由Richtmyer^[1]提出。在不可压缩假设下，RM不稳定性可以看作瞬时加速情况下的RT不稳定性。RT理论中对于一个单模小扰动界面，在重力作用下，界面扰动幅度a满足：

式中：t为时间，g为重力加速度，k = 2π/λ为扰动波数，a₀为界面的初始振幅，A为Atwood数。A的定义为：

式中：ρ₁和ρ₂分别为激波作用前界面上、下游流体密度。Richtmyer^[1]将RT理论中的重力加速度替换为脉冲形式g = ∆vδ(t)后，其扰动增长率可以写为：

式中：∆v为界面受激波冲击获得的阶跃速度，可以通过一维黎曼问题的求解^[71]得到。

为提高线性脉冲模型对实验的预测能力，Richtmyer^[1]进一步提出在线性理论中使用波后Atwood数及波后振幅，此时界面振幅的线性增长率为：

式中：A⁺和a⁺分别为激波刚离开界面后的Atwood数和振幅，波后Atwood数A⁺通过求解黎曼问题^[71]得到。波后振幅a⁺采取以下方式计算：

式中：U_t为透射激波速度。下文结果中RM脉冲模型指的是做了波后量修正的模型^[34]。

由RM脉冲模型可以看出，界面扰动幅度随时间线性增长，且增长速率是一个常数。在轻/重流体中，界面扰动幅度从与激波作用之后开始增长；而在重/轻流体中，界面扰动幅度先减小然后反向增长。为了进一步提高理论模型的预测能力，近年来有多种修正方法被提出，例如Meyer等^[17]提出用激波作用前后的平均振幅取代波后振幅（Meyer-Blewett (MB)模型）：

Vandenboomgaerde等^[18]假设在激波穿越界面的过程中Atwood数和振幅等参数线性变化（Van模型）：

Mikaelian^[72]借鉴了Layzer势流理论的思想，提出了适用于增长中期的非线性模型（Mik模型），其形式如下：

式中：线性增长率v₀由线性阶段模型得到，b为气泡，s对应尖钉；E_b/s与波后Atwood数有关^[73-74]。且：

Zhang等^[19]通过级数高阶展开来预测轻/重界面的扰动增长（Zhang-Sohn (ZS)模型），并通过Padé近似将模型有效范围扩展到了非线性阶段，其振幅增长率为：

式中：v_lin参考Liu等^[34]的做法取为线性模型增长率v₀。

Sadot等提出的非线性模型（Sad模型）^[20]将可压缩流中的ZS模型^[19]与Alon等^[75]提出的势流模型相结合，得到的扰动增长速率为：

式中：C的值通过匹配势流模型和实验结果进行确定，当A⁺ ≥ 0.5时，C = 1/(3π)；当A⁺<0.5时，C = 1/(2π)。

Dimonte等^[21]在A⁺比较大的情况提出了一个非线性理论（Dimonte-Ramaprabhu (DR)模型），该模型气泡或尖钉的扰动增长率为：

式中：系数${C_{{\text{b/s}}}} = \dfrac{{4.5 \pm \left| A \right| + (2 \mp \left| A \right|)\left| {k{a_0}} \right|}}{4}$，系数${F_{{\text{b/s}}}} = 1 \pm \left| A \right|$，下标b对应式中正负号的负号，下标s对应式中正负号的正号。

最近，Zhang等^[22]在有限密度比时统一了尖钉结构和气泡结构的增长率（Zhang-Guo (ZG)模型），其气泡和尖钉的增长率为：

式中：正Atwood数对应气泡，负Atwood数对应尖钉。本文中对比了冻结流中高马赫数激波冲击轻/重界面和重/轻界面时的数值模拟结果和理论模型结果（图7～8）。这里参考前人的做法，将理论曲线平移时刻t^*和振幅a^*，再同数值结果进行对比^[34]，t^*为界面被冲击到界面振幅线性增长所需要的时间，a^*为t^*时刻的数值振幅。轻/重界面t^*估计值为：

图  7  轻/重界面振幅的数值模拟结果与理论解的对比

Figure  7.  Comparison of the amplitudes of the light/heavy interface between numerical simulation results and theoretical solutions

下载: 全尺寸图片 幻灯片

图  8  重/轻界面振幅的数值模拟结果与理论解的对比

Figure  8.  Comparison of the amplitudes of the heavy/light interface between numerical simulation results and theoretical solutions

下载: 全尺寸图片 幻灯片

式中：U_f为反射激波速度。重/轻界面时t^*则直接取为扰动幅度降低到最小值的时刻^[35]。

图7(a)和图8(a)中的理论解基于波后振幅，图7(b)和图8(b)中的理论解基于平均振幅。从图7(a)可以看到，对于冻结流中轻/重界面，基于波后振幅的ZS模型能相对较好地预测高马赫数激波作用下的振幅增长规律，优于基于平均振幅的ZS模型（图7(b)）。但对于重/轻界面，基于波后振幅的理论结果与数值结果相差较大（图8(a)），而基于平均振幅的ZS模型能相对较好地预测振幅增长规律（图8(b)）。

3.2.2   不同气体模式的影响

轻/重界面时不同气体模式下的振幅演化结果见图9(a)，振幅增长率结果见图9(b)。由图9(a)可知，初始激波的压缩导致振幅首先出现减小，在激波经过界面后，界面振幅才逐渐增大。由图9(b)可知，增长率先升高后逐渐降低。重/轻界面时，不同气体模式下振幅演化结果见图10(a)，振幅增长率结果见图10(b)。由图10(a)可知，由于初始激波的压缩，不同气体模式下振幅均出现减小，激波经过界面后界面振幅先反向，随后逐渐增大。图10(b)中振幅增长率先升高，在非线性效应影响下，后期增长率随后出现降低。

图  9  不同气体模式中轻/重界面振幅及其振幅增长率对比

Figure  9.  Comparison of amplitudes and amplitude growth rates of the light/heavy interface among different gas models

下载: 全尺寸图片 幻灯片

图  10  不同气体模式中重/轻界面振幅及其振幅增长率对比

Figure  10.  Comparison of amplitudes and amplitude growth rates of the heavy/light interface among different gas models

下载: 全尺寸图片 幻灯片

不论是在轻/重还是重/轻界面的结果中，振幅和振幅增长率最大的是冻结流，最小的是热化学非平衡流，而热非平衡流的结果介于两者之间，说明热效应和化学反应都会减弱扰动的增长。由脉冲理论可知，在波长和初始振幅一定的情况下，振幅增长率正比于由激波撞击导致的界面初始跳跃速度和Atwood数。由3.1节可知，热化学非平衡流、热非平衡流和冻结流的界面跳跃速度依次降低，而相应的Atwood数依次增大。以M_s = 8的轻/重界面为例，热化学非平衡流的界面跳跃速度为3 078.6 m/s，热非平衡流的界面跳跃速度为2 882.1 m/s，冻结流的界面跳跃速度为2733.4 m/s，依次降低；而热化学非平衡流的Atwood数为0.148，热非平衡流的Atwood数为0.204，冻结流的Atwood数为0.232，依次增大。可见，随着气体模式从冻结流过渡到热化学非平衡流，扰动增长变慢，根据脉冲理论，界面跳跃速度和Atwood数的联合效应将影响扰动增长，但两者具有相反的变化规律，因此，不同气体模式下密度的不同是造成扰动增长变化的主因，而不是跳跃速度。对比振幅增长率可以看出，热效应和化学反应对振幅增长率峰值影响较大，但3种气体模式的扰动振幅增长率在后期差距较小。

3.2.3   初始扰动尺度的影响

通过固定波长而改变振幅的方式来研究初始小振幅和大振幅对扰动增长规律的影响是常用的研究手段，而本文的出发点是对比热化学非平衡流动同经典冻结流的区别，初始振幅相对波长保持在小振幅范围内，且初始振幅与波长比例保持不变。为便于对比不同初始扰动尺度对结果的影响，采用以下方式对时间进行无量纲化^[34]：

为便于比较取冻结流理论解，振幅的无量纲化公式为：

由轻/重界面的冻结流结果可以看到，在无量纲情况下，不同初始扰动尺度具有相同的振幅曲线（见图11(a)），这同冻结流中欧拉方程无内秉尺度的事实相对应。但对于热化学非平衡流，初始扰动尺度越大，扰动振幅增长越慢（见图11(b)），该结论在无量纲或有量纲时均成立。热化学非平衡流与冻结流不同，在热化学非平衡流中入射激波波后流场不再是均匀的，在到达平衡状态前存在一段非平衡距离（见表2），这里波后非平衡区域的流向距离L_eq定义为从激波波后压力到达p_eq − ( p_eq − p_ps )/1000时距离激波面的距离，其中p_eq为平衡压力，p_ps为初始激波波后（post-shock）压力（即冻结流中波后压力）。在相同的非平衡距离下，当初始扰动尺度增加时，非平衡距离与初始扰动尺度比值下降，波后状态越接近波后无穷远的平衡流状态，与冻结流结果差距越大。重/轻界面冻结流中，不同初始扰动尺度下振幅曲线同样重合（见图12(a)），而热化学非平衡流也是初始扰动尺度越大，扰动振幅增长越慢（见图12(b)），与轻/重界面情况的分析类似。

图  11  不同初始扰动尺度下轻/重界面振幅增长对比

Figure  11.  Comparison of amplitude growths of the light/heavy interface among different initial disturbance scales

下载: 全尺寸图片 幻灯片

表  2  不同马赫数非平衡距离

Table  2.  Non-equilibrium distances at different Mach numbers

界面类型	Ms	Leq/mm		界面类型	Ms	Leq/mm
轻/重界面	6	582.5		重/轻界面	8	569.4
	7	129.6			9	133.0
	8	41.4			10	47.5
	9	20.0			11	20.6

下载: 导出CSV 

| 显示表格

图  12  不同初始扰动尺度下重/轻界面振幅增长对比

Figure  12.  Comparison of amplitude growths of the heavy/light interface among different initial disturbance scales

下载: 全尺寸图片 幻灯片

热化学非平衡气体中不同初始扰动尺度的变化也会影响同理论解的对比情况。由于非平衡距离的存在，3.2.1节中的理论模型中的界面阶跃速度和Atwood数不再是一个常数。但考虑随时间变化的界面阶跃速度和Atwood数的理论研究还有待进一步开展，本文中只考虑热化学非平衡流在2种特殊情况下的界面阶跃速度和Atwood数，一种是通过紧靠初始激波面处的波后状态（也即冻结气体状态）计算的黎曼问题的界面状态，另一种是用初始激波波后无穷远处的平衡状态在热化学平衡气体假设下计算的界面状态（计算方法可参考一般气体模型中的黎曼问题计算方法^[71]）。对于M_s = 8的轻/重界面算例，数值模拟结果与基于冻结流状态和平衡流状态的理论解对比见图11(b)。可以看到，当λ = 150, 1500 mm时，基于平衡流状态的理论解比基于冻结流状态的理论解更靠近数值模拟结果，而当初始扰动尺度减小时（λ = 15 mm），数值模拟的结果逐渐靠近冻结流的理论解。产生这种现象是因为初始激波波后的非平衡距离是恒定的，初始扰动尺度变小会使得界面附近的流体对界面扰动的影响作用增强，而初始扰动尺度变大时，该尺度下影响界面的初始激波波后流体逐渐进入到平衡流状态，因此数值解更靠近平衡流理论解。对于M_s = 9的重/轻界面算例，基于冻结流状态和平衡流状态的理论解与数值模拟结果的对比见图12(b)。3种初始扰动尺度情况下都是基于平衡流状态的理论解比基于冻结流状态的理论解更靠近数值模拟结果，这与轻/重界面时不一致。注意到，在图7(b)中，对比了各种不同的理论模型，发现重/轻界面时均表现较差，而ZS模型相对较好。因此，在重/轻界面时，以上差异可能同理论模型的较差预测能力有关。

3.2.4   激波马赫数的影响

对于轻/重界面，不同M_s时冻结流中的振幅演化如图13(a)所示，而热化学非平衡流中的振幅演化在图13(b)中给出。Samtaney等^[59]研究发现，不同M_s时的振幅增长曲线在以下的时间尺度化方式下几乎重合：

图  13  不同激波马赫数时轻/重界面振幅增长对比

Figure  13.  Comparison of amplitude growths of the light/heavy interfaces at different Mach numbers of shock waves

下载: 全尺寸图片 幻灯片

在本文冻结流的数值模拟结果中，上述时间尺度化方法的确可以较好地归一化不同M_s时的振幅（见图13(a)），这同Samtaney等的结论^[59]一致。但在热化学非平衡流中，不同M_s下的振幅增长在该尺度化方法下并不重合（见图13(b)）。重/轻界面冻结流的振幅增长如图14(a)所示，热化学非平衡流的振幅增长如图14(b)所示，同样可以看到该尺度化方法也未能很好归一化热化学非平衡流的结果。

图  14  不同激波马赫数时重/轻界面振幅增长对比

Figure  14.  Comparison of amplitude growths of the heavy/light interface at different Mach numbers of shock waves

下载: 全尺寸图片 幻灯片

3.3   涡量场分布与环量演化规律

对于单模RM不稳定性的扰动增长机理主要有2种解释：一种是压力扰动机制^[76]，即界面与激波作用后，在非平直的界面附近产生压力扰动，压力扰动会影响界面的发展；另一种是斜压机制^[77-78]，认为由界面处涡量的产生和演化诱导界面失稳。因此，研究涡量场对认清高马赫数激波作用下单模RM不稳定性具有重要意义，本节将对涡量场演化和环量进行分析。

涡量${\boldsymbol{\varOmega }}$的定义为：

在二维可压缩无黏流动中，涡量${\boldsymbol{\varOmega }}$的生成仅与斜压效应相关^[79]，即：

而环量Г的定义为：

式中：L为区域边界，dl为区域边界的长度微元，u为速度矢量，S为积分区域，dS为面积微元，ω为垂直面积微元dS的涡量分量。

全场环量的计算采用全计算域的涡量积分，全场正（负）向环量为正（负）向涡量区域的积分，界面区域环量的积分域为标记函数f = 0.05～0.95的区域。不同气体模式下，涡量的形成和分布较类似，下文以冻结流为例，介绍其主要结果。轻/重界面时涡量场与涡量生成项（式(32)右端项）如图15所示。高马赫数激波与界面作用时涡量有2种生成机制，第1种是初始激波与密度不同的界面相互作用时产生的斜压效应（见图15(a)），轻/重界面生成负向涡量，重/轻界面形成正向涡量，这种机制是低马赫数激波与界面作用时的主要机制。第2种涡量生成机制是由透射激波波后的横波与透射激波相互作用而在透射激波面附近生成大量的涡量（见图15(b)），横波引起的涡量随后会形成明显的涡街（见图15(d)）。第2种机制形成的涡量在高马赫数情况下比较明显。本文中也做过低马赫数下的RM不稳定性数值模拟（未给出），发现横向扰动带来的涡量强度远小于界面区域。重/轻界面冻结流涡量场与涡量生成项如图16所示，同样存在这2种涡量生成机制（见图16(a)～(b)），并同样会形成涡街（见图16(d)）。

图  15  不同时刻冻结流中轻/重界面涡量场和涡量生成项

Figure  15.  Vorticity field and generation terms of the light/heavy interface in frozen gas at different times

Top: vorticity field; bottom: vorticity generation term.

下载: 全尺寸图片 幻灯片

图  16  不同时刻冻结流中重/轻界面涡量场与涡量生成项

Figure  16.  Vorticity field and generation terms of the heavy/light interface in frozen gas at different times

Top: vorticity field; bottom: vorticity generation term.

下载: 全尺寸图片 幻灯片

对比图17(a)轻/重界面的全场环量、全场正向和负向环量可知，全场环量的振荡趋势与全场正向、负向环量变化趋势一致。正向环量增加引起全场环量增加，负向环量（幅值）增加引起全场环量减少。对比涡量场可以确认环量的震荡与透射激波波后的涡街对应。图17(a)中全场环量平均振荡周期为53 μs，而采用流场纵向距离75 mm除以透射激波波后流场声速（一维理论预测1 499.7 m/s）估计的横波传播周期约为50 μs。而重/轻界面时（见图17(b)），全场环量平均振荡周期为48.5 μs，用透射激波波后流场声速1659.9 m/s估计的横波传播周期为45.2 μs，两者时间同样接近，说明环量场的震荡同透射激波波后横向波传播相关。

图  17  冻结流中轻/重和重/轻界面环量的演化

Figure  17.  Evolution of circulations of the light/heavy and heavy/light interfaces in frozen gas

下载: 全尺寸图片 幻灯片

轻/重界面时不同气体模式的全场和界面区域环量演化如图18(a)所示。由图可知热化学非平衡流环量幅值最小，冻结流最大，热非平衡流介于之间，且热化学非平衡流全场环量震荡衰减较快。上述结论同样适用于重/轻界面（见图18(b)）。从单模RM不稳定性扰动增长的斜压机制角度而言，热化学非平衡流界面附近环量幅值比冻结流小，对应界面区域涡量较小，因此定性上可断定涡量诱导的扰动增长比冻结流较慢，这同3.2.1节中发现的热化学非平衡流的振幅增长比冻结流慢是对应的。

图  18  不同气体模式轻/重和重/轻界面环量演化

Figure  18.  Evolution of circulations of the light/heavy and heavy/light interfaces in different gases

下载: 全尺寸图片 幻灯片

4.   结　论

对高马赫数激波冲击空气介质中因温度分层形成的单模界面而产生的RM不稳定性现象进行了数值模拟研究，考虑了轻/重和重/轻界面2种情况。高马赫数激波冲击界面时，会产生高温热化学反应，为揭示其中的热化学非平衡效应，对比了冻结流、热非平衡流和热化学非平衡流的区别，发现在波系结构、界面运动、扰动增长、涡量和环量规律等方面，高温气体的热效应（振动能激发）和化学反应均会对结果产生影响，导致了很多与经典冻结流不同的规律。相关结论如下：

(1)冻结流、热非平衡流和热化学非平衡流中透射激波、反射波及界面速度显著不同。轻/重界面时三者运动速度均为热化学非平衡流最快，热非平衡流次之，冻结流最慢；重/轻界面时透射激波波速为冻结流最快，热非平衡流次之，热化学非平衡流最慢，而界面和反射膨胀波运动规律则与轻/重界面相同。非平衡流相较于冻结流，界面和波系的相对快慢同基于冻结流和平衡流的一维黎曼问题的理论解是一致的。

(2)相对于其他常用理论模型，ZS模型能更好预测高马赫数激波作用下的单模RM不稳定性的振幅增长规律。对于3种气体模式，冻结流中振幅及其增长率最大，热非平衡流次之，热化学非平衡流最小。当初始扰动尺度变化时，冻结流中的无量纲振幅增长曲线可归一化，但该结论在热化学非平衡流中不成立。而当激波马赫数改变时，由Samtaney等^[59]提出的经典的时间尺度化方法虽然可以较好地归一化冻结流中不同激波马赫数下的振幅增长曲线，但不适用于热化学非平衡流。

(3)高马赫数激波冲击时，在2个区域中有较强的涡量生成，一个位于界面区域（这与低马赫数时相同），另一个位于透射激波波后。波后区域形成的原因是透射激波波后存在较强的横波，其与透射激波相互作用产生较强涡量。横波的来回传播和相互作用形成涡街状的涡量分布，并且引起环量的震荡。不同气体模式中环量幅值也明显不同，冻结流的环量幅值最大，热非平衡流次之，热化学非平衡流最小。不同气体模式中环量幅值的大小规律定性上同上述振幅增长快慢规律对应，这种发现符合扰动振幅增长的涡诱导机制的解释。

由以上结论可以看到，在热化学非平衡流中，高马赫数激波作用下的单模RM不稳定性振幅增长及其影响因素的作用规律同经典冻结流的规律有很多区别，因此本文的研究可以补充对强激波与单模界面作用而形成的RM不稳定性的认识。