不同流体间的扰动界面在冲击波作用下会发生不稳定性增长，称为Richtmyer-Meshkov（RM）不稳定性^[1-2]。当RM不稳定性发展进入非线性阶段后，扰动界面会发展为“蘑菇”状的尖钉/气泡结构；之后由于界面两侧切向速度不匹配而产生剪切效应，导致界面上又发生Kelvin-Helmholtz（KH）不稳定性^[3]，并诱发更多小尺度涡。KH不稳定性的出现会破坏“蘑菇”状界面结构的对称性，使之破碎，并加速湍流混合。流体界面不稳定性和湍流混合涉及许多重要领域，如武器物理、惯性约束聚变、超音速燃烧、天体物理等，因此具有重大的研究意义。

较早之前，对于RM界面不稳定性和湍流混合，研究人员关注的重点是湍流混合区宽度的增长规律及影响因素^[4-16]。随着实验诊断技术、高精度数值算法及并行计算技术的发展，研究重点逐渐深入到湍流混合的内在机制和统计规律^[17-23]以及以此为基础的湍流混合建模^[24]。然而，由于湍流混合的多尺度特性，实验中所能获取的流场信息依然十分有限，无法给出足够多的湍流混合区瞬态数据和精细复杂的混合区流场结构；而数值模拟基本上能够给出所需要的有关湍流混合区流场的所有信息，特别是借助超大规模的高性能计算，可以在不同尺度上研究湍流混合的物理机制及其演化的统计特性。Cabot等^[25]就曾在IBM的BlueGene/L超级计算机上利用6万多个核实现了Rayleigh-Taylor（RT）不稳定性及湍流混合的近300亿网格（3 072³）的超大规模直接数值模拟，研究了湍流混合区增长的标度率问题和雷诺数效应。

AWE开展的“反尖端”RM不稳定性激波管实验中初始的SF₆/空气右界面呈“∨”形，整体上类似于反的士兵臂章，称为“反尖端”界面。本研究采用自研的可压缩多介质黏性流动和湍流大涡模拟程序MVFT^{[6, 21]}，对“反尖端”界面不稳定性及湍流混合进行大规模并行数值模拟，在高分辨率的网格下细致地研究“反尖端”界面不稳定性的发展过程和规律，包括波系作用、湍流混合区结构的演化、湍流混合区统计量的变化等。

1.   计算方法

可压缩多介质黏性流动和湍流大涡模拟程序MVFT所求解的流场控制方程是经过Favre滤波后的可压缩多介质黏性流动Navier-Stokes方程组

式中：下标i、j分别代表x、y、z 3个方向，遵循张量运算法则； $\bar \rho$ 、 ${\tilde u_k}(k = i,j)$ 、 $\bar p$ 、 $\tilde E$ 分别为可解尺度流体的密度、速度、压强和单位质量总能量；N为流体种类； ${\tilde Y^{(s)}}$ 表示第s种流体的体积分数（s= $1,2,\cdots\!,$ N–1），且满足 $\displaystyle\sum\nolimits_{s=1}^N {{{\tilde Y}^{(s)}}} = 1$ ； $\tilde D$ 为扩散系数， $\tilde D = {\nu / {S\!c}}$ ， $\nu$ 为流体运动黏性，Sc为Schmidt数； $q_j^{\text{l}} = - {\lambda _{\text{l}}}{{\partial \tilde T} / {\partial {x_j}}}$ 是单位时间单位空间的可解尺度能量流， ${\lambda _{\rm{l}}} = {\mu _{\rm{l}}}{c_{{p}}}/P{r_{\rm{l}}}$ 是可解尺度热传导系数，c_p为比定压热容，Pr_l是Prandtl数， $\tilde T$ 为流体温度； ${\tau _{ij}},\;Q_i^{\text{T}},\;Q_i^{\text{Y}}$ 为滤波后产生的亚格子输运项，分别表示亚格子应力张量、热流量和体积分数输运通量，需要通过亚格子尺度应力模型进行模化处理以使控制方程组封闭，本研究采用Vreman亚格子模型^[26]。牛顿流体黏性应力张量 ${\sigma _{ij}}$ 表示为

式中：μ_l为流体黏性。气体介质的状态方程采用理想气体状态方程。

具体计算过程为：首先采用算子分裂技术将（1）式描述的物理过程分解为3个子过程进行计算，即整个通量分解为无黏通量、黏性通量和热通量3部分。无黏通量的计算采用多介质高精度PPM（Piecewise Parabolic Method）方法，两步Lagrange-Remapping型的PPM方法分4个步骤进行：（1）物理量分段抛物插值，（2）近似Riemann问题求解，（3）Lagrange方程组推进求解，（4）将物理量映射到静止的欧拉网格上。然后在无黏通量的基础上，采用二阶空间中心差分方法和两步Runge-Kutta时间推进方法求解黏性通量、热通量及标量输运通量。

2.   结果分析

本研究所采用的模型是英国AWE开展的“反尖端”RM不稳定性激波管实验^[27]，计算模型如图1所示，由细金属丝制成的网状支架支撑微米量级的硝化纤维薄膜将实验气体隔开，形成预设的初始界面。激波管内截面尺寸为20 cm×10 cm，“反尖端”区气体为SF₆，其宽度为15 cm，两侧气体是空气，“反尖端”界面平衡位置距离激波管尾端20 cm。另外，在初始界面上设置振幅小于0.01 cm、波长为0.5～5 cm的小尺度随机扰动，以模拟支架的扰动效应。初始入射冲击波马赫数为1.26。计算域：(x, y, z)∈[0, 80 cm]×[–10 cm, 10 cm]×[0, 10 cm]。实验得到的结果十分有限，只是不同时刻混合区的演化图像及壁面气泡和尖钉位置（距激波管末端的距离），不足以深入分析界面不稳定性及湍流混合发展规律。为此，本研究借助大规模数值模拟，利用1 024个CPU，采用粗（网格大小1.00 mm）、中（网格大小0.50 mm）、细（网格大小0.25 mm）3种不同分辨率的网格，详细研究“反尖端”界面不稳定性及湍流混合的演化特性。细网格下的网格总量为1.024×10⁹。气体初始参数如表1所示，其中γ为气体比热比。

图  1  计算模型和“反尖端”界面

Figure  1.  Computational model and inverse chevron interface

下载: 全尺寸图片 幻灯片

表  1  空气和SF₆的初始参数

Table  1.  Initial properties of air and SF₆

Gas	ρ/（kg·m–3）	p/MPa	γ	μl/（Pa·s）	Diffusion coefficient/（m2·s–1）
SF6	5.97	0.1	1.09	1.474 6×10–5	0.97×10-5
Air	1.18	0.1	1.40	1.852 6×10–5	2.04×10–5

下载: 导出CSV 

| 显示表格

冲击波从低阻抗介质（如空气）向高阻抗介质（如SF₆）方向加载物质界面时，会产生一个透射冲击波和一个反射冲击波；冲击波从高阻抗介质向低阻抗介质方向加载物质界面，会产生一个透射冲击波和一个反射稀疏波。图2给出了以密度显示的“反尖端”计算模型轴线上的一维近似波谱，可以看出复杂的波与界面以及波与波之间的相互作用，冲击波加载使流场密度增大，稀疏波卸载使流场密度减小。首先，入射冲击波从左侧加载空气/SF₆左界面，产生一个左行反射冲击波和一个右行透射冲击波；该右行透射冲击波加载SF₆/空气右界面，在SF₆气体中产生一个左行反射稀疏波，向右边空气中透射一个右行冲击波；右边空气中的右行冲击波会从激波管尾端反射回来再次与SF₆/空气右界面作用，向SF₆中产生一个左行透射冲击波（会在左界面上发生分解，向SF₆中反射一个右行稀疏波），并向右边空气中产生一个右行反射冲击波，该右行反射冲击波会反复从激波管尾端反射回来加载SF₆/空气右界面，在界面上发生分解；前面第一次在SF₆/空气右界面上分解向SF₆中产生的左行反射稀疏波和空气/SF₆左界面作用后向SF₆中产生一个反射压缩波，该压缩波和SF₆/空气右界面作用分解产生一个左行反射稀疏波，之后在SF₆左右界面之间反复分解成压缩波和稀疏波。如图2所示，在SF₆气体区域的胞格结构显示了冲击波、压缩波、稀疏波、界面之间的复杂相互作用，使得界面不稳定性发展更复杂，而波和界面的每一次作用都会促使混合区快速发展^[22]。

图  2  用流场密度显示的一维近似波谱图

Figure  2.  1D approximate wave visualized using flow filed density

下载: 全尺寸图片 幻灯片

图3给出了“反尖端”界面及流场演化的实验图像和数值模拟结果比较，左列是实验结果，右3列从左至右依次是网格为1.00、0.50、0.25 mm时数值模拟的沿展向平均的密度场图像，清晰地显示了波和界面的相互作用及波系在流场中的相互作用和演化过程。实验和数值模拟二者之间，包括界面的形状和位置及波阵面的形状和位置均吻合很好。当冲击波与SF₆区域左右界面作用后，界面上的小尺度随机扰动开始增长；而且在冲击波与SF₆区域的右界面作用时，由于冲击波是从高阻抗介质向低阻抗介质方向加载界面，所以右界面的大尺度“反尖端”界面发生反相，从中心轴处生长出大尺度尖钉，靠近激波管上下壁处形成大尺度气泡结构，并且尖钉和气泡逐渐长大。对比3种不同网格分辨率数值模拟的展向平均密度场可以看出，界面大尺度结构和波系的演化差别很小，只是在中后期界面小尺度结构之间有较小的差别，这是由于展向平均抹去了一些细节结构造成的。图4、图5和图6分别给出了2.0、3.0、4.0 ms时刻不同网格分辨率下以SF₆体积分数（Y_SF6）显示的湍流混合区三维图像。图中显示了复杂的三维湍流混合区演化，包括小尺度气泡/尖钉的长大、变形、融合等；不同分辨率网格计算所得的湍流混合区图像差别很大，其中高分辨率网格捕捉到了更精细的湍流混合区结构。图7给出了大尺度壁面气泡和中心尖钉位置随时间变化的曲线，该位置是相对于激波管尾端的距离，不同颜色曲线表示不同网格分辨率下的数值计算结果。可以看出，在不同网格分辨率下计算所得的大尺度壁面气泡和中心尖钉的位置差别较小，这是因为大尺度壁面气泡和尖钉的位置由湍流混合区的包络线确定，而湍流混合区的包络线主要由大尺度结构决定，湍流混合区大尺度结构的平均特征在不同网格分辨率下的差别很小；虽然差别较小，但是细网格下的模拟结果与实验结果吻合得更好。

图  3  “反尖端”界面演化的实验图像（左列）和以密度显示的数值模拟结果（三维计算的展向平均， 右3列从左至右计算网格尺寸依次为1.00、0.50和0.25 mm）比较

Figure  3.  Comparison of experimental (left column) and simulated density images (three right columns on different grid resolutions) of inverse chevron interface

下载: 全尺寸图片 幻灯片

图  4  2.0 ms时刻不同网格分辨率下以SF₆体积分数显示的湍流混合区三维图像

Figure  4.  3D images of turbulent mixing zone visualized using SF₆ volume fraction on different grid resolutions at 2.0 ms

下载: 全尺寸图片 幻灯片

图  5  3.0 ms时刻不同网格分辨率下以SF₆体积分数显示的的湍流混合区三维图像

Figure  5.  3D images of turbulent mixing zone visualized using SF₆ volume fraction on different grid resolutions at 3.0 ms

下载: 全尺寸图片 幻灯片

图  6  4.0 ms时刻不同网格分辨率下以SF₆体积分数显示的的湍流混合区三维图像

Figure  6.  3D images of turbulent mixing zone visualized using SF₆ volume fraction on different grid resolutions at 4.0 ms

下载: 全尺寸图片 幻灯片

图  7  大尺度壁面气泡和中心尖钉的位置D随时间变化曲线

Figure  7.  Positions of wall-bubble and center-spike with large scale in time

下载: 全尺寸图片 幻灯片

图8和图9显示了1.0、2.0、3.0、4.0 ms时中心轴上流场密度和SF₆体积分数在不同网格分辨率下的分布。图中清晰地显示了冲击波的传播及其在界面上的分解，以及界面的发展演化，比如1.0 ms时刻的密度分布显示，冲击波加载使流场密度增大，SF₆气体中的右行透射冲击波与SF₆/空气界面作用后向SF₆气体中反射的左行稀疏波的卸载效应使流场密度减小。在前期（1.0～2.0 ms），不同网格分辨率计算的流场密度和SF₆体积分数的分布差别很小，说明混合区流场湍流发展不充分，小尺度脉动还未完全发展起来；后期（3.0～4.0 ms）不同网格分辨率的计算结果差别增大，表明在冲击波的多次复杂加载作用下混合区流场已经进入完全湍流阶段，更多小尺度脉动出现，细网格对这些小尺度含能涡结构有更高的分辨率，而粗网格将大量的小尺度含能涡结构耗散掉。

图  8  不同时刻中心轴上的流场密度分布

Figure  8.  Flow density distributions along the centerline at different times

下载: 全尺寸图片 幻灯片

图  9  不同时刻中心轴上SF₆体积分数分布

Figure  9.  SF₆volume fraction distributions along the centerline at different times

下载: 全尺寸图片 幻灯片

湍动能和拟涡能可以用来表征流场湍流脉动和湍涡的强度和发展。图10和图11分别为1.0、2.0、3.0、4.0 ms时刻无量纲化的流场湍动能K和拟涡能Ω在不同网格分辨率下沿冲击波运动方向的分布，无量纲化速度参数取初始入射冲击波后的流场速度133.6 m/s，无量纲化长度参数取计算模型宽度0.1 m。可以看出：随着混合区的发展，小尺度脉动逐渐发展起来；随着网格分辨率的增大，对这些小尺度脉动的捕捉能力也逐渐提高，捕捉到了更强的湍流脉动和湍涡。所以更小尺寸的网格对流场有更高的分辨率，可以捕捉到更精细的混合区湍流结构，这对于研究RM不稳定性及湍流混合演化的动力学行为非常重要。

图  10  不同时刻无量纲化湍动能沿冲击波运动方向的分布

Figure  10.  Dimensionless turbulent kinetic energy distributions along motion direction of shock wave at different times

下载: 全尺寸图片 幻灯片

图  11  不同时刻无量纲化拟涡能沿冲击波运动方向的分布

Figure  11.  Dimensionless enstrophy distributions along motion direction of shock wave at different times

下载: 全尺寸图片 幻灯片

3.   结　论

利用可压缩多介质黏性流动和湍流大涡模拟程序MVFT，对“反尖端”界面不稳定性及其诱发的湍流混合进行了大规模并行三维数值模拟分析，揭示了界面不稳定性及湍流混合的复杂发展过程和规律。冲击波和空气/SF₆/空气界面作用并发生分解而产生冲击波、稀疏波、压缩波，这些次生波在SF₆气体中运动并相互作用，而且多次加载物质界面以及波和界面的每一次作用都会加速湍流混合区的发展和物质混合。高/低阻抗构型的“反尖端”界面受冲击加载后发生反相而发展为大尺度的壁面气泡和中心轴尖钉结构，该大尺度气泡和尖钉结构基本确定了湍流混合区的平均几何特征和包络范围而不依赖计算网格。采用高分辨率的计算网格能够捕捉到更精细的混合区小尺度湍涡结构及更强的湍流脉动，显示出湍流混合区的复杂结构和特征。