当冲击作用于不同密度或不同可压缩性2种物质的扰动界面时产生Richtmyer-Meshkov不稳定性(RMI)。这种不稳定性理论上由Richtmyer发现并描述^[1], 由Meshkov从实验中证实^[2]。该界面不稳定性问题在许多自然现象及科学和工程领域中起着重要作用^[3-8], 如超新星爆炸、磁化等离子体、磁约束、太阳磁化层、地下盐矿、火山岛及外壳与内部流体混合导致中子收益降低的激光驱动惯性约束聚变和冲击波与火焰相互作用导致的爆燃转爆轰等。此外, RMI也可能从受冲击的金属表面产生喷射物。

RMI的演化通常经历由不稳定模式的振幅h_k和波长λ=2π/k描绘的若干阶段。对于kh_k≪1, 扰动随kUt呈正比例增长, U为激波作用后的界面运动速度。当kh_k达到某一值, 非线性使增长率降低, 驱动模式耦合, 且增长率随着扰动谱宽的增大而减小。然后, 由于尖钉下落(重流体进入轻流体)比气泡上升(轻流体进入重流体)快, 界面变得不对称。对于宽的不稳定谱, 非线性最终导致产生湍流混合层。RMI的脉冲性质令问题复杂, 使得RMI定性上与常见的Rayleigh-Taylor不稳定性(RTI)不同。由于冲击的可压缩性、复杂的物质特性以及后期的非线性运动直至湍流混合, RMI演化的计算是困难的。当然, 随着计算机技术的迅猛发展, 这可以采用多维高分辨率流体力学模拟来进行, 但它们计算强度大, 无法用于工程设计优化研究。因此, 目前实际应用中通常采用捕捉较低分辨率时不稳定流动主要特征的简化“混合模型”^[9]。杨玟等^[10-11]对此进行了尝试, 将传统的k-ε模型应用于界面不稳定性引起的混合, 取得了令人满意的结果。

但是, 由于与RMI相关的其它物理过程非常复杂, 较复杂的混合模型(如k-ε)也难以直接应用到工程设计中。目前, 很多实际应用中对RMI诱发混合现象的处理都非常简单, 假设混合层宽度以指数形式t^θ_i增长。而大量实验研究表明该比例关系仍不确定^[3-6], 因为考虑压缩性的计算是困难的, 它们与实验不符。即使指数律粗略满足, 但不同工况下θ_i的差别也很大, 它显著依赖于初始扰动谱。由此可见, 工程设计中对RMI诱发混合现象的处理过于粗糙。

本文中, 在简单介绍描述作用于混合层中产生的气泡和尖钉的浮阻力模型基础上, 采用该模型对激波管低压缩情况和激光加载高压缩情况下的RMI诱发混合层宽度(气泡与尖钉宽度之和)进行计算, 验证模型和选取参数的有效性。

1.   模型介绍

目前, 典型的浮阻力模型可写为如下形式^[12]:

式中：下角标i, j表示2种不同的流体，下角标为1时表示重流体(尖钉)，为2时表示轻流体(气泡)；ρ_i为重流体/轻流体的密度；v_i是尖钉/气泡的渗透速度，且v_i=dh_i/dt，h_i表示尖钉/气泡的瞬时宽度；C_a是附加质量力系数；β是浮力产生的模型常数；C_d是阻力系数；a(t)为激波脉冲加速度；V_i为尖钉/气泡的体积，A_i为尖钉/气泡的截面积。方程左端第一项为惯性力，第二项为浮力，右端为阻力。关于模型的详细论述可参考文献[13-14]，这里不再重复。对于Richtmyer-Meshkov不稳定性，通常认为冲击简单地给予界面上的气泡和尖钉一个脉冲，则它们随后的运动可以由惯性力和阻力相等来得到(加速度为零)。因此脉冲加速度情况是有启发性的，可以用来研究不稳定性的惯性特性。

本文所求解的模型方程是一组二阶常微分方程, 将它们简化为一阶微分方程:dh_i/dt=v_i; dv_i/dt= -f_iC_dv_i|v_i|/h_i。采用四阶Runge-Kutta方法进行求解。

2.   结果分析与讨论

采用上述模型和数值方法, 对关注的激波管低压缩情况和激光加载高压缩情况下模型的性能进行了考察。这2种工况下RMI产生的机理不同:对于弱冲击, 主要贡献来自于压力梯度和密度梯度不重合引起的旋涡沉积; 对于强冲击, 存在激波在经折射后产生了显著的反射, 这产生增长率的振荡, 但它们最终衰减。

2.1   激波管低压缩的情况下的模拟

首先采用上述模型对4种不同激波脉冲加速度情况下气泡和尖钉宽度进行了计算。图 1给出了所采用的4种加速度曲线，g为重力加速度。脉冲加速度a约为150g，持续时间t₀约为10 ms。这些曲线为LANL的Dimonte等LEM(Linear Electric Motor)实验的测量曲线^[15]。实验中流体和脉冲加速度的性质参数见表 1，其中R为密度比，R=(1+A)(1-A)，A为Atwood数，A=(ρ₂-ρ₁)(ρ₂+ρ₁)，We为韦伯数，Re为雷诺数。对于每一种情况，通过调整阻力系数C_d和初始振幅h_i0来使随时间变化的解与实验数据相符。但是，数值实验发现：在大多数情况下h_i0对结果的影响远小于C_d的影响。

图  1  计算采用的4种不同脉冲加速度曲线

Figure  1.  Four kinds of impulsive accelerations used in the calculation

下载: 全尺寸图片 幻灯片

表  1  实验中采用的流体和脉冲加速度性质参数

Table  1.  Fluid combinations and characteristics for impusive accerleration experiments

No.	流体1	流体2	ρ1/(g·cm-3)	ρ2/(g·cm-3)	R	A	We	Re
1	H2O	CCl2F2	1.000	1.57	1.57	0.22	4 000	2 600
2	SF6	C4H10	0.067	0.81	12.10	0.85	1 100	8 000
3	SF6	CCl2F2	0.067	1.57	23.40	0.92	11 000	23 000
4	SF6	CCl2F2	0.032	1.57	49.10	0.96	6 000	25 000

下载: 导出CSV 

| 显示表格

图 2给出了4种加速度驱动下气泡和尖钉宽度随位移Z的变化, $Z=\iint a \mathrm{~d} t^{\prime} \mathrm{d} t$, 激波作用时Z≈Ut。由图可见, 4种加速度情况下计算的气泡和尖钉宽度与实验基本吻合。计算中阻力系数C_d的取值为3.67±0.73, 与文献[16]中分析得到的C_d的不确定度1.2接近。从图中还可看出:气泡和尖钉的不对称性随着密度比R的增大而增大。此外, 本文中还对实验结果按指数律h_i=h_i0t^θ_i进行了拟合, 其中h_i0的取值范围为0.5~1.0 cm。R=49.1时, θ₁≈0.85, θ₂≈0.33; R=23.4时, 指数迅速下降, θ₁≈0.45, θ₂≈0.24; R=1.57时, θ₁≈0.28, θ₂≈0.22。由此可见, 指数θ_i随密度比变化而变化, 但具体变化规律还未从数值模拟和实验中最终确定, 这主要是由于θ_i对实验初始条件敏感, 需要计算和实验之间更直接的比较。

图  2  气泡和尖钉宽度随位移的变化

Figure  2.  The width of bubble and spike with displacement

下载: 全尺寸图片 幻灯片

2.2   激光加载高压缩情况下的模拟

为了考察模型在高压缩情况下的性能，我们进一步对Nova激光器上马赫数Ma>10的实验进行了模拟。实验采用一靶丸装置在Nova激光器上进行^[17]。流体1由厚度为125 μm、初始密度为1.7 g/cm³的铍烧蚀层组成。流体2是未压缩密度为0.12 g/cm³的泡沫。波速为46 km/s的入射冲击与界面相互作用产生反射稀疏波和速度为3 km/s的透射激波。界面经加速后速度为56 km/s，物质被压缩后，ρ₁=2 g/cm³，ρ₂=0.5 g/cm³，A=-0.6。这些参数通过对比热比γ₁=1.8和γ₂=1.45的流体求解理想的黎曼问题得到。

图 3给出了Nova实验中计算的加速度曲线。由图可见，激光驱动在4 ns后停止，这导致泡沫减压，由于A < 0而产生Rayleigh-Taylor(RT)分量，因此冲击压缩后流动是亚音速的，本文模型是适用的。图 4给出了混合区总宽度H随位移Z的变化(由于实验不能分辨气泡和尖钉，因此给出了总振幅H)。从图中可看出：混合区总宽度的计算值与实验值吻合，而且C_d=2.0和C_d=5.36的曲线之间包括了全部的实验数据。但是，阻力系数C_d的不确定度约为3.36，明显大于低压缩情况的值(约为1.46)。此外，拟合得到总的混合宽度以指数为0.5的指数律增长，这超过了激波管低压缩时得到的指数，推测其原因可能是：(1)激光驱动随时间减小，使得压力降低、界面减速，这导致扰动膨胀，并引入RT分量(因为Aa>0)。这些影响可能显著增加推测的指数；(2) A=0.6时Nova上的初始扰动比激波管上的更对称，如果指数对初始条件敏感，这可能导致不同的指数。

图  3  Nova实验中的加速度曲线

Figure  3.  Acceleration history for Nova experiment

下载: 全尺寸图片 幻灯片

图  4  混合区宽度随位移的变化

Figure  4.  Variaion of total width with displacement

下载: 全尺寸图片 幻灯片

3.   结论

采用浮阻力模型对激波管低压缩和激光加载高压缩情况下Richtmyer-Meshkov不稳定性诱发的物质渗透边界的演化过程进行了计算, 计算结果与实验吻合得较好。这表明本研究中模型参数的选取、方程中现象学比例因子的添加和模型假设是合适的。但是由于实验测量的局限性, 模型中的一些问题仍然是突出的, 包括阻力项的大小和形式、压缩的影响、“附加质量”的描述等。为了更好地评估模型, 需要一些实验上的完善。首先, 气泡和尖钉必须单独分辨, 因为它们的表现相当不同, 尤其在A较大的情况。其次, 实验持续时间应当延长至足以揭示模型的差别为止。尽管如此, 本文模型仍明显优于当前实际应用中所采用的经验公式(本研究也显示指数θ_i随工况的不同而显著变化)。