流场非均匀性对非平面激波诱导的Richtmyer-Meshkov不稳定性影响的数值研究*

王 震，王 涛，柏劲松,2，肖佳欣

（1.中国工程物理研究院流体物理研究所，四川 绵阳 621999；2.中国工程物理研究院流体物理研究所冲击波物理与爆轰物理实验室，四川 绵阳 621999；3.北京临近空间飞行器系统工程研究所，北京 100076）

RM不稳定性[1]及其诱导产生的湍流混合现象，由于其广泛的应用如大特征尺度的超新星爆炸[2]、中等特征尺度的超燃冲压发动机中超声速燃烧[3]及小特征尺度的惯性约束聚变[4]等，而备受关注。其产生的关键机制为斜压涡量沉积，即激波加载两种密度不同的物质界面时，界面处密度梯度与压力梯度不共线而导致斜压涡量沉积（即 ∇ ρ×∇p/ρ2≠0），进而界面处的初始扰动会被放大而产生失稳现象。

半个多世纪以来，学者们对RM不稳定性进行了广泛深入的研究，结果表明界面、流场及激波是影响RM不稳定性发展的重要因素。Yang等[3]通过48组数据数值研究，发现激波强度、密度比、椭圆界面形状等9种因素将影响界面的不稳定性发展，他的研究为后续RM不稳定性影响因素的研究起了指导作用。Zou等[5]、Bai等[6]分别实验、数值模拟研究，表明平面激波加载长/短轴比值不同的椭圆形气柱界面时，界面处沉积的涡量不同，因而椭圆形气柱界面的形态存在差异。Bai等[7]以高斯分布函数描述初始流场密度的非均匀分布，数值研究表明初始流场非均匀性对平面激波诱导的RM不稳定性演化影响显著。Xiao等[8]、肖佳欣等[9]在此基础上更加细致地分析了初始均匀流场与非均匀流场中的涡量及环量分布，揭示了非均匀流场中RM不稳定性产生及演化机理。同时，Bai等[10]通过数值模拟两种非均匀流场中平面激波诱导的RM不稳定性现象，发现：在线性阶段或弱非线性阶段，界面扰动增长强烈依赖于初始流场非均匀性；而在非线性阶段，该依赖性逐渐减小。

可以发现，Bai等[7,10]、Xiao等[8]、肖佳欣等[9]侧重研究流场非均匀性对平面激波诱导产生的RM不稳定性演化的影响。然而，Ishizaki等[11]、Kane等[12]、Zou等[13]和Zhai等[14]分别通过理论分析、实验研究及数值模拟手段发现，均匀流场中非平面激波冲击密度不同的平界面时也会产生类RM不稳定性现象。实际上，在工程应用中非平面激波广泛存在，如在RM不稳定性重要应用领域之一的ICF中，由于非均匀激光的辐照或者靶丸外表面粗糙度的影响，均会不可避免的产生非平面激波。同时，受制造工艺精度的影响，靶丸内部不同物质界面可能存在微小的扰动，并且物质的密度并非绝对均匀分布。因此，在内爆压缩阶段，非均匀流场中非平面激波冲击扰动界面引起的不稳定性就有可能出现，进而影响靶丸的压缩效果，甚至造成点火失败。目前，通过实验手段实现稳定、可靠的非平面激波生成及其参数测定仍然是一大难点，而且物质内部的密度分布很难精确测量。因此，通过数值模拟手段研究非均匀流场中非平面激波加载扰动界面的RM不稳定性演化规律，就具有重要的学术价值，这对于拓展对界面不稳定性现象的认识具有重要的帮助。

在文献[6-10, 15-18]的基础上，本文中，通过大涡模拟方法，数值模拟δ1=0.616 2 m、δ2=0.496 1 m非均匀流场中非平面激波加载扰动Air/SF6界面时的RM不稳定性现象。通过界面形态、扰动振幅、环量及脉动速度相关量，重点分析初始流场非均匀性对非平面激波诱导的RM不稳定性演化的影响，以期为系统地开展非平面激波诱导的RM不稳定性影响因素理论分析及实验研究提供参考数据。

1 数值方法与计算模型

1.1 数值方法

采用可压缩多介质黏性流动和湍流大涡模拟程序MVFT数值求解经Favre滤波后的可压缩多介质黏性流动Navier-Stokes方程组，如下：

式中：i、j代表x、y、z的3个方向分别代表可解尺度流体密度、速度、压强及单位质量总能量；N为流体种类数；为第s种流体体积分数，满足为扩散系数，ν为流体运动黏性系数，Sc为Schmidt数；为牛顿流体黏性应力张量，δij为Kronecker符号为亚格子应力张量；分别为单位时空可解尺度、亚格子能量流， λl= µlcp/Prl、 λt= µtcp/Prt分别为可解尺度、亚格子热传导系数，µl、 µt为流体动力黏性系数，cp为比定压热容，为流体温度，Prl、Prt为Prandtl数。采用Vreman亚格子模型[19]模化处理，实现控制方程组封闭。

采用理想气体状态方程，如下：

式中：γ代表气体比热比，对于空气γ=1.40，e代表比内能。

采用算子分裂技术将方程组（1）描述的物理过程分解成无黏通量、黏性通量及热通量3个部分进行计算。无黏通量采用多介质高精度分段抛物线方法（PPM）进行计算，两步Lagrange-Remapping型PPM方法由以下4个步骤组成：（1）物理量分段抛物插值；（2）近似Riemann问题求解；（3）Lagrange方程组推进求解；（4）将物理量映射回静止的欧拉网格上。在无黏通量的基础上，采用二阶空间中心差分及两步Runge-Kutta时间推进方法求解黏性通量、热通量及标量输运通量，详细步骤见文献[20]。

1.2 计算模型

采用图1所示计算模型，计算域为[-0.02 m,0.25 m]×[0.00 m, 0.20 m]，离散网格数为 540×400。左边界为自由来流，右、上、下边界为固壁。扰动air/SF6界面呈现正弦函数分布：

式中：初始平衡位置xi0=0.03 m，波长λi=0.05 m，初始振幅Ai=0.005 m。

马赫数为1.25的非平面激波初始波阵面呈现正弦函数分布：

式中：波阵面平衡位置xs0=0.01 m，波长λs=0.05 m，初始振幅As=0.005 m，初始相位为φ。

本文中，SF6气体初始密度呈现高斯函数分布：

式中：yc=0 m；δi(i=1, 2)为流场非均匀系数，反映了流场密度分布的非均匀程度，δi越小，流场的非均匀性越强，密度分布越不均匀；选取δ1=0.616 2 m、δ2=0.496 1 m，即在计算域上部气体密度为 ρSF6，而下部气体密度则分别为 0.9ρSF6、 0.85ρSF6，如图2所示。初始时刻气体参数见表1。

图1 计算模型Fig.1 Simplified computational model

图2 均匀流场及非均匀流场沿y方向的密度分布Fig.2 Density distributions along y direction in uniform and non-uniform flows

表1 Air/SF6初始参数（20 ℃, 101.325 kPa）Table 1 Properties of air and SF6 gases

2 结果与分析

讨论的重点在于初始流场密度的非均匀分布对非平面激波导致的RM不稳定性演化的影响，对于非平面激波与平面激波引起的RM不稳定性的差别仅进行简单对比，其产生机制差异更加深入的分析将在后续工作中开展。

2.1 流场可视化

图3给出了均匀流场及δ1=0.616 2 m、δ2=0.496 1 m非均匀流场中，非平面激波同相及反相加载扰动Air/SF6界面时的流场密度云图。同时，还给出了均匀流场中，相同初始位置的平面激波冲击同一物质界面的演化过程作为对比。当入射激波初始冲击界面后，流动迅速进入线性阶段，界面扰动振幅逐渐增大而形状变化较小。约0.2 ms后非线性效应逐步增强，流动进入弱非线性阶段，界面上形成典型的“气泡”及“尖钉”结构，如图3(b)～(d)所示。约1.6 ms反射激波再次冲击已经变形的界面，非线性效应占据主导作用，流动进入强非线性阶段，界面开始出现多尺度涡结构，流动呈现湍流混合状态，如图3(e)～(f)所示。可以看到，非平面激波波后流场是非均匀分布的，且其波阵面形状也随时间发生变化，而平面激波波后则为均匀流场，如图3(a)所示。均匀流场中，透射激波波阵面为垂直于激波运动方向的平面，“气泡”及“尖钉”结构呈轴对称分布；非均匀流场中，透射激波波阵面则为斜面，“气泡”及“尖钉”结构非对称增长，且δi越小，其对称性越差。同时，无论非平面激波同相加载还是反相加载，非均匀流场中的界面形态在反射激波加载前与均匀流场差异明显，如图3(b)～(d)所示；而反射激波加载后这种差异减小，如图3(e)～(f)所示。此外，反射激波加载前，平面激波加载引起的“气泡”及“尖钉”结构流向尺寸小于非平面激波同相加载时的尺寸，但大于其反相加载时的数值。

图3 密度云图Fig.3 Simulated density contours

2.2 扰动振幅分析

图4 界面扰动振幅演化历史Fig.4 Perturbation amplitude evolution histories of interface

为定量描述流场非均匀性对非平面激波诱导的RM不稳定性的影响，定义沿激波运动方向、SF6气体体积分数为1%～99%区域宽度的一半为界面扰动振幅。界面扰动振幅的演化过程可以宏观地表征RM不稳定性中物质界面的变形及后期湍流混合的程度的变化。入射激波初始加载及反射激波加载后界面扰动振幅的演化过程如图4所示，图4(b)中tre=0 ms即对应反射激波加载后界面扰动振幅再次增长的时刻。从图4(a)可以看出，由于非平面激波初始加载，Air/SF6界面扰动振幅开始增长，两种非均匀流场中界面扰动振幅均大于均匀流场中的数值，并且δi越小，其界面扰动振幅偏离均匀流场条件下的数值越远。约1.6 ms时刻反射激波到达并再次冲击变形的物质界面，1.9 ms时刻（即图4(b)中的0 ms时刻）反射激波与界面作用完成，界面扰动振幅再次增长，而且增长速度更快。反射激波加载后的湍流混合阶段，均匀流场与两种非均匀流场中界面扰动振幅增长规律性不强，但二者的差异总体上呈减小趋势。前述规律说明，非平面激波无论是同相加载还是反相加载，在反射激波加载前阶段，初始流场的非均匀性会显著影响界面失稳的发展；而反射激波加载后阶段，虽然界面扰动增长加剧了，但是流场逐步趋向于均匀分布，因而初始流场的非均匀性对界面扰动振幅增长的影响有所减小。此外，从图中还可以看到，平面激波加载引起的界面扰动振幅大小总体上介于非平面激波同相加载与反相加载时的数值之间。

2.3 环量分析

在RM不稳定性中，涡量的生成与分布决定了界面的演化，在二维流动中，忽略质量力、体积力、黏性效应及可压缩效应的涡量输运方程为：

这说明流场中密度梯度与压力梯度方向的不共线是流体微团运动过程中涡量变化的主要原因。二维流动中，涡量表达式为：

式中：u、v分别为x、y方向速度分量。涡量在整个计算域中的分布决定了环量的分布，环量是解释界面扰动演化快慢的统计量之一，其定义为速度矢量沿封闭曲线的积分，由斯托克斯定理可变换为封闭曲线所围成面积上涡量的积分。即：

式中：L、A分别代表计算域[-0.02 m, 0.25 m]×[0.00 m, 0.20 m]所围成的封闭曲线及面积。

图5给出了均匀流场及δ1=0.616 2 m、δ2=0.496 1 m非均匀流场中，非平面激波同相及反相加载时流场中正、负及总环量的变化规律。总环量反映了激波与界面相互作用过程中整个流场生成的涡量的多少，总环量为负表示整个流场中产生的正涡量小于负涡量，反之则表示正涡量大于负涡量。从图中可以看出，均匀流场中正、负环量高度对称增长，总环量为0；非均匀流场中，负环量主导了总环量的大小，总环量并不为0。同时，图5还展示了均匀流场中平面激波诱导下流场中的环量变化。可以看到，由于非平面激波x方向速度分量u存在沿y方向的速度梯度，因而在入射激波冲击界面前，非平面激波波阵面处存在涡量，从而流场中正、负环量不为0，而入射激波为平面激波时，流场中正、负环量则从0开始增长，这个特点是非平面激波与平面激波诱导的RM不稳定性最为明显的区别。此外，平面激波加载时，流场中的正、负环量明显小于非平面激波引起的数值。

图5 环量演化历史Fig.5 Evolution histories of circulations

为了进一步分析流场非均匀性对非平面激波诱导产生的RM不稳定性演化的影响，引入了非均匀流场与均匀流场中正、负环量的相对差值，表示非均匀流场中环量偏离均匀流场的程度。定义为：

图6 正、负环量相对差值Fig.6 Relative differences of positive and negative circulations

表2 正、负环量相对差值极大值Table 2 Maximum values of relative difference for positive and negative circulations

图7 沿x方向的脉动速度一点二阶相关量分布Fig.7 Correlation distribution of fluctuating velocities along x direction

式中分别表示非均匀流场、均匀流场中的正、负环量。图6给出了非平面激波同相和反相加载时，δ1=0.616 2 m、δ2=0.496 1 m非均匀流场中正、负环量与均匀流场中相对差值的变化规律。同时，表2列出了反射激波加载前、过渡阶段及反射激波加载后两种非均匀流场中正、负环量与均匀流场相对差值的极大值。可以看出，反射激波加载前，负环量的相对差值大于正环量的相对差值，并且δi越小，差别越明显。这种现象产生的原因为，由于初始流场的非均匀分布，增加了物质界面处的密度梯度，导致界面区域各方向速度发生改变，使得流场中生成的负涡量大于正涡量，从而流场中负环量值大于正环量值[8]；而且，界面处密度梯度随着δi减小而增大（见图2），因而流场中负涡量值与正涡量值的差异相应增大，环量也呈现相同的规律。此外，δ2=0.496 1 m非均匀流场中的正、负环量相对差值大于δ1=0.616 2 m非均匀流场中的数值。在过渡阶段，两种非均匀流场与均匀流场的正、负环量相对差值急剧增大，仍然是负环量的相对差值大于正环量的相对差值。反射激波加载后，两种非均匀流场与均匀流场中正、负环量相对差值明显减小。这些说明，流场非均匀性引起的环量差别主要存在于反射激波加载前及过渡阶段，而反射激波加载后差异则显著减小。由于流场中涡量主要集中在变形的界面区域，环量的分布则由涡量分布决定。因此，反射激波加载后，两种非均匀流场与均匀流场正、负环量差别的缩小就可以解释反射激波加载后初始流场非均匀性对界面扰动振幅影响减小的现象。

2.4 脉动速度相关量分析

脉动速度相关量可以表征激波加载后流场脉动的剧烈程度，其实质是通过速度的脉动统计量来反映界面的演化。在二维流动中，脉动速度一点二阶相关量表示脉动运动的平均动量输运：

式中分别代表速度分量u、v的脉动量。图7给出了反射激波加载前、后流场中沿x方向的脉动速度一点二阶相关量分布，其中图7(b)、(d)分别是图7(a)、(c)中均匀流场条件下非平面激波及平面激波加载引起的脉动速度相关量的局部放大图。由图可以看出，反射激波加载前，非均匀流场中非平面激波导致的脉动速度一点二阶相关量比均匀流场中的数值高出约两个数量级，并且δi越小，非均匀流场中脉动速度一点二阶相关量越大，如图7(a)、(c)。反射激波加载后，非均匀流场与均匀流场中脉动速度一点二阶相关量差别缩小至一个数量级，而且两种非均匀流场中的脉动速度一点二阶相关量相对大小差异也较反射激波加载前有所减小，如图7(e)、(f)。这与前述界面扰动振幅及环量的演化规律相吻合，这从高阶统计量的角度分析了初始流场非均匀性对非平面激波诱导的RM不稳定性演化发展的影响。此外，平面激波引起的流场脉动速度一点二阶相关量总体上小于非平面激波诱导的速度脉动量数值，图7(b)、(d)。

3 结 论

运用MVFT程序，数值模拟了均匀流场及两种不同系数的非均匀流场中非平面激波同相及反相加载扰动界面时的RM不稳定性现象，并与平面激波诱导的RM不稳定性进行了简要对比分析。主要考察了界面扰动振幅、环量及流场脉动速度的演化规律，具体结论如下。

（1）无论非平面激波同相加载还是反相加载，入射激波加载后至反射激波加载前，流动进入线性与弱非线性阶段，界面的不稳定性增长对初始流场的非均匀性十分敏感；流场非均匀系数δi越小，初始流场密度分布越不均匀，非平面激波加载后，流场中脉动速度相关量及负环量数值越大，界面扰动振幅越大。反射激波加载后，非线性效应占据主导作用，流场逐步趋向于均匀分布，初始流场非均匀性对界面演化的影响逐渐减弱。

（2）非平面激波与平面激波诱导的RM不稳定性有着明显的区别：入射激波冲击界面前，非平面激波波后呈现非均匀流场，且波阵面形状随时间发生变化。此外，由于x方向速度分量u存在y方向的速度梯度，从而非平面激波波阵面区域存在涡量，该涡量与非平面激波加载物质界面时生成的涡量共同作用，使得流场中环量及速度脉动与平面激波加载时存在显著的差异，因而呈现出与平面激波诱导的RM不稳定性不一样的界面演化过程。