运动激波冲击环境气体(如空气、氮气)中具有不同密度的气泡会引发Richtmyer-Meshkov(RM)不稳定，在气泡界面诱导产生不同尺度的涡，使界面失稳变形，进而促进两侧不同气体相互掺混，这种现象广泛存在于自然界和工业领域，如超新星爆炸^[1]、超燃推进^[2]等。由于气泡界面两侧气体物性参数的差异，激波在冲击气泡时还可出现透射、绕射、聚焦、反射等复杂的波系变化，这会对气泡界面造成压力扰动，影响气泡界面的失稳过程。因此，相关研究的开展具有非常重要的实际应用价值和理论意义。

国内外学者已对激波与气泡作用过程开展了较多实验与数值研究。Haas^[3]采用阴影照相技术获得了弱激波分别与He轻气泡和R22重气柱相互作用的图像，并将激波扫过轻/重气体的情形归类为发散类型(divergent case)和收敛类型(convergent case)，但由于一次只能拍摄一张结果，故他们的实验需要较好的重复性。Layes等^[4-7]研究了不同强度的激波与球形气泡作用，其中气泡内分别采用小于、接近和大于外部空气的气体(分别为He、N₂、Kr)，他们通过高速摄影技术追踪气泡的变形，获得了不同密度气泡结构清晰的演变过程。Giordano等^[8]对Layes等^[4]的实验进行了二维数值模拟，发现了Kr重气泡在激波聚焦作用下形成的射流结构，但未对射流引发机理进行分析。同一时期，Tomkins等^[9]采用平面激光诱导荧光技术(planar laser-induced fluorescence, PLIF)技术获得了激波与单个SF₆重气柱作用的结果，并测量了气柱内外不同气体之间的混合率。Haehn等^[10]实验研究了N₂中的激波分别与Ar和SF₆气泡作用的情形，考察了反射激波作用对气泡变形的影响。Zhai等^[11]采用高速纹影和二维VAS2D程序对入射激波与SF₆重气泡的作用过程进行了实验与数值研究，获得了详细的流场演化过程，由于采用的高速摄影相机拍摄频率较高(时间间隔低于30 μs，小于Layes等^[4-7]在实验中采用的70 μs)，发现了重气泡下游界面产生的射流。Zhu Yuejin等^[12]、朱跃进等^[13]采用大涡模拟方法，对入射及反射激波与He轻气泡和Kr重气泡的作用也进行了数值研究，分析了轻、重气泡形成涡环的三维变形失稳过程。

然而不难看出，以往研究主要考察气泡界面的形态演变、结构尺寸变化及界面内外气体混合等问题。事实上，由于重气泡内的声阻抗较大，进入重气泡内的透射激波传播速度比外部绕射激波慢，内部透射激波和外部绕射激波能够在气泡下游极点附近分别发生碰撞，并伴随有局部温度和压力的瞬时升高，这种现象即为激波聚焦，这种特殊的流动现象类似于点爆炸，有可能在气泡下游极点附近引发高速气体射流，使重气泡界面出现更复杂的失稳变形。但由于射流结构较小，已有实验或计算的时空分辨率不够，且射流结构在流场黏性作用下不断衰减等原因，较长时间内重气泡内激波聚焦引发射流的现象都没有得到重视。有部分学者开始讨论激波聚焦引发射流这一现象。邹立勇等^[14-15]采用二维无黏的VAS2D程序数值研究了平面激波冲击不同密度的球形重气泡(SF₆、R22、R12和Kr)的情形，考察了不同的气泡内外声阻抗比对激波聚焦和射流形成的影响。Zhai等^[11]、Zhai Zhigang等^[16]同样采用VAS2D程序研究了激波在SF₆重气泡中聚焦和引发射流的现象，认为射流形成主要受气泡处压力扰动影响。王革等^[17]用二维欧拉方程数值研究了R22重气柱在激波作用下形成射流的现象，从波系结构角度解释了射流形成的原因。沙莎等^[18-19]对激波与R22重气柱和SF₆气泡的作用进行了数值模拟研究，结果显示气柱/气泡的RM失稳过程，并对激波在重气体内聚焦引发射流的现象进行了分析。

综合上述文献可以发现，目前重气泡内激波聚焦引发射流的研究仍存在一些不足，如多数文献未考虑流场黏性，且激波聚焦和射流引发的具体过程尚未完全清楚。为了揭示相关流动现象的细节，流场结构的时空分辨率需进一步提高。因此，本文中将采用高精度计算格式和高网格分辨率对平面激波冲击SF₆重气泡的过程进行数值研究，详细分析激波在重气泡内的聚焦过程，并进一步澄清射流的形成机理，以期为后续研究提供参考。

1.   计算方法与模型

1.1   计算方法

计算采用单组份二维Navier-Stokes方程，具体表达为：

式中：$\mathit{\boldsymbol{F}}\rm{=}\left[\begin{matrix} \mathit{\rho u} \\ \mathit{\rho }{{\mathit{u}}^{\rm{2}}}\rm{+}\mathit{p}\rm{-}{{\mathit{\tau }}_{\mathit{xx}}} \\ \mathit{\rho uv}\rm{-}{{\mathit{\tau }}_{\mathit{yx}}} \\ \rm{(}\mathit{E}\rm{+}\mathit{p}\rm{)}\mathit{u}\rm{-}\mathit{k}\frac{\partial \mathit{T}}{\partial \mathit{x}}\rm{-}\mathit{u}{{\mathit{\tau }}_{\mathit{xx}}}\rm{-}\mathit{v}{{\mathit{\tau }}_{\mathit{yx}}} \\ \mathit{\rho uY}\rm{-}\mathit{D\rho }\frac{\partial \mathit{Y}}{\partial \mathit{x}} \\ \end{matrix} \right]$, $\begin{align} & \mathit{\boldsymbol{G}}\rm{=}\left[\begin{matrix} \mathit{\rho v} \\ \mathit{\rho uv}\rm{-}{{\mathit{\tau }}_{\mathit{xy}}} \\ \mathit{\rho }{{\mathit{v}}^{\rm{2}}}\rm{+}\mathit{p}\rm{-}{{\mathit{\tau }}_{\mathit{yy}}} \\ \rm{(}\mathit{E}\rm{+}\mathit{p}\rm{)}\mathit{v}\rm{-}\mathit{k}\frac{\partial \mathit{T}}{\partial \mathit{y}}\rm{-}\mathit{u}{{\mathit{\tau }}_{\mathit{xy}}}\rm{-}\mathit{v}{{\mathit{\tau }}_{\mathit{yy}}} \\ \mathit{\rho vY}\rm{-}\mathit{D\rho }\frac{\partial \mathit{Y}}{\partial \mathit{y}} \\ \end{matrix} \right] \\ & \\ \end{align}$, $\mathit{\boldsymbol{U}}\rm{=}\left[\begin{align} & \mathit{\rho } \\ & \mathit{\rho u} \\ & \mathit{\rho v} \\ & \mathit{E} \\ & \mathit{\rho Y} \\ \end{align} \right]$, $\mathit{\boldsymbol{W}}\rm{=}\frac{\mathit{v}}{\mathit{y}}\left[\begin{align} & \ \ \mathit{\rho } \\ & \ \mathit{\rho u} \\ & \ \mathit{\rho v} \\ & \mathit{p}\rm{+}\mathit{E} \\ & \ \mathit{\rho Y} \\ \end{align} \right]$。其中：ρ为密度，u、v分别为x和y方向的速度，E为单位体积总能量，p为压力，γ为质量热容比，Y为环境气体的质量分数，τ为黏性应力张量，热传导系数k、扩散系数D、运动学黏性ν按照环境气体性质给出，$\mathit{E}\rm{=}\frac{\mathit{p}}{\mathit{\gamma }\rm{-1}}\rm{+}\frac{1}{2}\mathit{\rho }\rm{(}{{\mathit{u}}^{\rm{2}}}\rm{+}{{\mathit{v}}^{\rm{2}}}\rm{)}$。在上述控制方程组中，空间导数项的黏性部分采用十阶空间中心差分计算，无黏部分采用高精度九阶WENO格式^[20]计算，时间推进过程采用三阶Runge-Kutta方法求解，计算过程保证满足CFL数值稳定性条件。

1.2   计算模型

计算以文献[11]中实验为依据，设计如图 1所示的初始流场，流向(x)长0.25 m，法向(y)高0.07 m。SF₆重气泡位于(x=0.024 m，y=0.035 m)处，半径R₀=0.015 m，密度ρ₀=6.14 kg/m³，气泡外环境气体为空气，密度ρ₁=1.19 kg/m³，SF₆和空气具体的物性参数可见表 1。初始时刻，入射激波沿x方向从左向右传播，马赫数为1.23，激波后气体状态由Rankine-Hugoniot关系给出。初始流场的上、下边界为固壁，采用无滑移刚性绝热壁面边界；x方向左端面为进口，采用零梯度边界；x方向右端面为出口条件。实际计算区域取初始流场上半部分，中心对称轴为对称边界，网格尺寸为0.05 mm×0.05 mm。

图  1  初始流场示意图

Figure  1.  Schematic of initial flow

下载: 全尺寸图片 幻灯片

表  1  气体参数

Table  1.  Gas parameters

气体	ρ/(kg·m-3)	γ	c/(m·s-1)	ρc/(kg·m-2·s-1)	M
空气	1.19	1.40	346.0	411.7	29
SF6	6.14	1.09	133.9	822.1	146

下载: 导出CSV 

| 显示表格

2.   结果与讨论

2.1   实验验证

为验证数值方法的可靠性，图 2给出了不同时刻的实验纹影与数值纹影对比图。通过对每组纹影图的对比可以看出，入射激波的波系变化、各时刻SF₆重气泡的界面变形均与实验结果基本一致，这验证了计算方法的可靠性。需要注意的是，由于数值方法与计算网格的精度均较高，数值纹影图中气泡界面处有明显的不同尺度的涡，且SF₆重气泡下游的射流结构长度比实验结果长。

图  2  实验与数值结果对比

Figure  2.  Comparisons between experiment and simulation

下载: 全尺寸图片 幻灯片

为进一步比较实验与数值计算的结果，图 3~4给出了气泡上游极点位置和射流长度随时间的变化曲线，图中x表示上游极点坐标，x₀表示初始时刻气泡上游极点位置，气泡直径D₀=0.03 m，L_t为气泡上游极点距射流头部的距离，特征时间t_c≈249 μs^[11]，气泡界面位置按照环境质量分数分别为0.01、0.50和0.99进行处理。从图中可以看出，在入射激波与SF₆重气泡作用及之后的一段时间内，不同Y值所对应的计算结果与实验结果均吻合较好；当时间进一步发展，不同Y值对应的计算结果彼此之间呈现差异，只有Y=0.01时的气泡上游极点位置和射流长度变化一直与实验较吻合。由于Y值表示的是环境气体的质量分数，故在初始时刻，SF₆重气泡内Y=0，外界环境Y=1，当气泡受激波冲击时，气泡界面处的斜压效应和速度剪切会使界面发生Richtmyer-Meshkov不稳定和Kelvin-Helmholtz不稳定，在界面处形成不同尺度的涡，从而促进气泡内外气体的混合，使流场中出现Y介于0~1之间的混合气体层，而Y=0.01则表示所选位置处只有少量环境的空气与SF₆掺混，其密度较接近纯SF₆气泡，可近似把此位置看作气泡边界，因此，实验结果与Y=0.01时的计算结果吻合一定程度上反映出实验拍摄的局限性，即对气泡内外不同密度气体混合较好区域的捕捉效果稍差，这一点在图 2的可视化结果对比中亦有所体现, 图 3~4的定量对比表明数值方法的可靠性和合理性。

图  3  气泡上游极点位置随时间的变化规律

Figure  3.  Variations of upstream pole with time

下载: 全尺寸图片 幻灯片

图  4  气泡上游射流长度随时间的变化规律

Figure  4.  Variations of jet length with time

下载: 全尺寸图片 幻灯片

2.2   波系与SF₆气泡界面演化

为了对入射激波冲击SF₆重气泡引发射流的过程进行详细分析，图 5给出了不同时刻流场内复杂波系和SF₆气泡界面的数值纹影图。在图 5(a)所示时刻，从左向右运动的入射激波(IS)恰好运动到气泡中心位置，由于SF₆的声阻抗大于外界空气的声阻抗，故进入SF₆气泡内部的透射激波(TS)运动速度较慢，此时流场中还存在一道向左传播的弧状反射激波(RS)；到图 5(b)时，入射激波(IS)继续向右运动，在此过程中于界面外侧产生绕射激波(DS)，而绕射激波(DS)则在气泡内形成了向气泡流向中心对称轴运动的入射绕射激波(DTS)，此时透射激波(TS)进一步弯曲，和DTS以及气泡下游界面围成一个未受波系扰动的局部区域(UZ)，故此处的压力、密度等物理量与初始时刻一致；随着时间推进，气泡外侧两道相向运动的绕射激波(DS)先发生碰撞，形成新的透射激波，紧接着气泡内侧上下对称的2道入射绕射激波(DTS)也发生碰撞，形成1道向左传播的小激波(SS1)，在这段时间内，未受波系扰动的局部区域(UZ)周围波系透射激波(TS)和入射绕射激波(DTS)处于汇聚过程，其强度均有所增加，故透射激波(TS)和入射绕射激波(DTS)的波后流场物理性质的差异逐渐增大，形成了1道新的小激波(SS2)，而未受波系扰动的局部区域(UZ)的面积则逐渐减小，见图 5(c)、(d)所示；图 5(e)~(h)时，由透射激波(TS)、入射绕射激波(DTS)和左传播的小激波(SS1)围绕的未扰动区(UZ)进一步缩小，波系继续靠拢并发生聚焦，出现类似点爆炸过程，在气泡下游界面内侧形成高压膨胀区(EZ)；随后，围绕高压膨胀区(EZ)的激波向四周传播，使高压膨胀区(EZ)不断扩大，其中靠近下游极点附近的激波冲击下游界面后发生透射，形成向左传播的反射稀疏波(RRS)和向右传播的新透射激波(TSS)，从而加强了对SF₆重气泡下游极点附近界面的扰动，使其附近出现小的凸起，见图 5(i)~(k)所示，到t=160.4 μs时，气泡下游极点附近已明显形成射流头部，如图 5(J)所示。

图  5  波系与SF₆气泡演化过程

Figure  5.  Evolution of waves and SF₆ gas bubble

下载: 全尺寸图片 幻灯片

2.3   激波聚焦和射流引发分析

2.3.1   密度与压力分布

为具体分析激波在SF₆重气泡内的聚焦过程，图 6给出了与图 5(e)~(h)所示的4个时刻对应的密度场云图，图中黑色实线为压力等值线。在图 6(a)所示时刻，气泡内侧未受激波扰动的区域(UZ)向内收缩呈竖狭长状，围绕未受激波扰动的区域(UZ)的波系在其法向方向的上下顶点处首先汇聚形成2个高压区(HZ)和1道新激波(SS2)(与图 5(d)中SS2一致)；到图 6(b)时，围绕未受激波扰动的区域(UZ)的波系在流向上进一步汇聚，使未受激波扰动的区域(UZ)不断收缩并最终消失，而高压区(HZ)的密度和压力值则迅速升高，此时，流场中仅有这对上下对称且相互靠近的高压区；随着高压区(HZ)的相向运动，到t=119.7μs时，上下对称的高压区(HZ)已经在气泡流向中心对称轴处发生碰撞汇聚，形成1个高压高密度的区域(HF)，通过分析流场发现高压高密度的区域(HF)的局部密度超过98 kg/m³，压力超过60个大气压，由于捕捉高压区(HZ)恰好碰撞形成焦点的时刻极为困难，且压力和密度峰值也会受计算步长和网格尺寸影响，故实际上激波在SF₆重气泡内聚焦时的密度和压力峰值比上述计算值更大，这表明SF₆重气泡能够对入射激波起到极强的汇聚作用，而不断扩大的高压高密度膨胀区(EZ)则会对气泡下游界面产生冲击和扰动，为产生向下游运动的射流提供条件。事实上，重气泡内外的声阻抗比值、入射激波强度都对气泡内的激波聚焦过程有重要影响，能够改变激波聚焦区域的压力峰值以及和气泡下游极点之间的相对位置，这都是关系到能否成功引发射流或者射流强度的重要影响因素，故有必要开展更多的研究，以揭示声阻抗和激波强度对激波聚焦和射流引发的影响。

图  6  密度场分布云图

Figure  6.  Contour of density field distribution

下载: 全尺寸图片 幻灯片

2.3.2   定量分析

为进一步分析激波聚焦过程引起的流场物理性质变化，图 7~8给出了SF₆重气泡中心对称轴(y=0.035 m)在x=0.01~0.07 m之间的密度和压力分布曲线。在图 7中，左侧和中间的虚线框内分别为SF₆气泡上游界面和下游界面附近位置，而最右侧的虚线框内为入射激波附近位置，可以看出，入射激波在所选的六个时刻内向下游方向运动明显。本节重点关注气泡下游界面附近的密度和压力变化，从图中可以看出，在t=114.5 μs时，中间虚线框内所对应位置处的密度和压力曲线均呈现两边高中间低的分布，这是因为此时下游界面内侧存在未受激波扰动的局部区域(UZ)，见图 5(d)，随着时间发展，围绕未扰动的区域(UZ)的波系逐渐汇聚，使得未扰动的区域(UZ)在流向上逐渐被压缩，到t=118.4 μs时未扰动的区域(UZ)已消失，而随着激波的不断汇聚，密度和压力曲线在t=119.7 μs时达到峰值，表明此时激波已经汇聚完成，密度和压力大大提高。在激波聚焦之后，聚焦区域开始向四周膨胀，中心的密度和压力值均有所下降，由于聚焦区域在紧靠气泡下游界面的内侧，故向下游运动的激波很快便传播至气泡下游界面，在t=126.7 μs时，气泡内外存在较高的压力差，对下游界面有较强的冲击作用，有利于引发射流结构。

图  7  不同时刻SF₆重气泡中心对称轴上的密度分布曲线

Figure  7.  Density profiles of SF₆ bubble along symmetric axis at different times

下载: 全尺寸图片 幻灯片

图  8  不同时刻SF₆重气泡中心对称轴上的压力分布曲线

Figure  8.  Pressure profiles of SF₆ bubble along symmetric axis at different times

下载: 全尺寸图片 幻灯片

2.3.3   涡量分布

在入射激波与SF₆重气泡作用过程中，气泡界面处既存在由密度梯度和压力梯度不一致引起的斜压效应，也存在由气泡内外速度差异引起的剪切作用，故气泡界面能够形成大量涡量，从而影响气泡内外不同气体介质的混合，进一步加剧气泡界面扰动。为研究涡量生成与演化对气泡下游界面射流形成的影响，图 9给出了不同时刻SF₆重气泡的涡量分布云图，其中黑色表示顺时针旋转，为负涡量，白色表示逆时针旋转，为正涡量。从图中可以看出，在t=103.1 μs时刻，气泡上表面为负涡量，下表面为正涡量；到t=118.4 μs时，激波在气泡内侧发生汇聚，聚焦区域的上半部分出现白色正涡量，下半部分出现黑色负涡量；随着聚焦完成后的高压高密度区域向四周不断膨胀，当向下游运动的激波冲击气泡下游界面时，气泡的上半部分仍形成逆时针旋转的正涡，下半部分形成顺时针旋转的负涡，见图 9(c)，此时涡对的旋转在气泡流向中心对称轴处能够诱导产生向下游的速度，从而加速气泡内侧流体冲击下游界面，到t=160.4 μs时，气泡下游极点处形成射流，且其头部的涡量分布仍为上正下负，即此时涡的自旋仍会促进射流向下游运动。由于涡对周围流体的卷吸作用，在射流发展过程中，其内部的空气质量分数Y会逐渐增大，射流结构本身也会出现不稳定现象，产生新的细小结构，要想进一步分辨射流结构本身的发展和失稳过程，需要精度更高的实验研究手段和数值方法，这有待进一步研究。

图  9  不同时刻SF₆重气泡的涡量分布云图

Figure  9.  Vorticity contours of SF₆ heavy bubble at different times

下载: 全尺寸图片 幻灯片

3.   结论

基于二维Navier-Stokes方程，采用高精度数值方法和高网格分辨率对平面入射激波冲击SF₆重气泡的作用过程进行了详细数值研究，主要考察了激波在重气泡内部的聚焦和射流引发过程，并通过与实验结果的对比验证了计算的可靠性。结果表明：

(1) SF₆重气泡内的激波聚焦过程分2个阶段，首先是波系在重气泡内侧未扰动区域法向的上下顶点处先汇聚形成一对高压区，随后高压区相向运动在流向中心对称轴处发生碰撞汇聚，完成聚焦过程；

(2) SF₆重气泡具有极强的聚能效应，激波聚焦形成的压力和密度峰值远大于初始时刻值，该高压可促进气泡下游界面形成射流；

(3) 激波聚焦在气泡下游界面附近还引起流场涡量变化，涡对的旋转对SF₆重气泡下游界面的射流形成也具有重要的促进作用。