平面激波与气泡界面相互作用过程的数值研究

叶 斌, 王 超, 施红辉, 张 珂, 熊红平

(浙江理工大学机械与自动控制学院, 杭州 310018)

平面激波与气泡界面相互作用过程的数值研究

叶 斌, 王 超, 施红辉, 张 珂, 熊红平

(浙江理工大学机械与自动控制学院, 杭州 310018)

对在空气中的平面激波与He气泡和SF6气泡相互作用过程进行了数值计算。利用大涡模拟(LES)方法,使用具有三阶精度的Third-Order MUSCL算法对气泡的发展与变形过程进行了模拟,然后与文献中的实验结果进行了比较。结果发现:无论激波通过低密度界面还是高密度界面,都会有射流结构和涡的产生,但产生的机理并不相同;对于密度低的界面,气泡会演化为两个对称的“气泡”结构,并产生两个独立的涡和持续时间较短的射流;而对于密度高的界面,气泡会演化为“尖钉”结构。

激波; 数值计算; 大涡模拟; 气泡和尖钉; 射流

0 引 言

激波与液体或气体界面相互作用时会出现不稳定现象,并加速流体间的相互混合,这种现象称为Richtmyer-Meshkov(RM)不稳定性[1-2]。RM不稳定性在天体物理、惯性约束核聚变(ICF)、水中炸药爆炸、航天火箭发动机中燃料的混合和燃烧等方面都有重要的应用价值。

对于RM不稳定性,已经有理论、实验、数值三方面的研究成果。其中理论模型方面的研究主要有Richtmyer[1]提出的脉冲模型理论、可压缩线性理论以及Zhang等[2]提出的非线性理论等。在实验研究方面,最早由Meshkov[3]在1969年完成了验证性的实验;其后Hass等[4]于1987年利用平面弱激波轰击球形界面及柱形界面,进行了一系列的实验研究,通过放射性照相技术,记录了气泡迅速弯曲变形直至破裂的整个过程。2003年Layers等[5]分别吹制了由氦、氮、氪气体材料充填的肥皂泡,研究了平面激波和球形界面的相互作用情况;得到的放射性照片显示,三种气体实验的球形界面变形差异较大,但在不稳定性发展的后期,气泡的运动速度会不断增加,直至一个有限值。2007年施红辉等[6]利用竖直放置的矩形截面激波管,通过CCD高速摄影,研究了气/液界面多元扰动RM不稳定性后期阶段形成的流体混合现象,得出了混合区域宽度与时间成线性关系的结论。2011年司廷等[7]采用高速纹影法,开展了平面激波以及反射激波作用下球形气体界面演化的实验研究,发现对于He气泡,当反射激波作用后,左右界面的运动会发生明显变化,右界面的速度迅速降低(53 m/s),产生右涡环的右界面速度(69 m/s)小于无反射时的下游界面速度(134 m/s),而左界面的位置几乎保持不变。对于SF6气泡,当壁面与气泡的距离较小时将产生两个射流,而且能够观察到界面振动;当距离增大后气体界面上将产生小尺度的涡结构。在数值模拟方面,1994年胡光初等[8]成功将二阶迎风TVD数值格式运用到弱电离、化学非平衡斜激波反射流场的数值计算中;2003年,宗文刚等[9]重新研究了双重加权实质无波动(DWENO)激波捕捉格式的构造过程,并对两重加权值的计算方法进行了改进,将改进后的DWENO(r=3)格式进一步应用于一维、二维典型问题的数值模拟;2009年李平等[10]将适用于可压多介质黏性流动和湍流的大涡模拟(LES)方法以及改进的MVFT代码,应用到了对低密度流体界面不稳定性及其诱发的湍流混合问题的数值模拟中;2010年翟金明等[11]采用二阶精度的单调迎风中心格式(MUSCL)和非结构自适应网格技术与有限体积形式,对轴对称变截面管道中基于主一次的爆轰波传播过程进行了数值模拟;2014年吴宇等[12]采用Fluent中的压力速度耦合的PISO算法,利用VOF方法研究了二维和三维下的液柱在激波马赫数为1.1的情况,成功的计算出了气液界面上RM不稳定性的演化过程,结果表明:初始扰动数目对RM不稳定性的影响显著,液柱会在横截面平面(Z方向)发生变形失稳,导致沿液柱轴向出现变形失稳,数值计算的结果与实验结果吻合较好,说明应用VOF方法于气液界面上RM不稳定性的研究,也能有较满意的结果。

本文针对两种不同密度的气体混合的情况,建立了多组分流动的二维轴对称模型,采用具有三阶精度的MUSCL算法,对气泡界面RM不稳定性现象进行大涡模拟(LES),分析了扰动界面的发展和演化规律及射流结构和涡结构的分布和发展,并且将数值结果与实验结果进行了比较。本文在前期计算中发现,VOF方法不适用于气体间的混合,而k-ε模型由于湍流度太大使得对混合后期的描述不理想,因此在本文中采用了大涡模拟。

1 数值计算模型

1.1 控制方程

本文利用大涡模拟(LES)的方法对可压缩多组分流进行模拟,并利用Fave技术对控制方程进行滤波,滤波后的控制方程如下:

1.2 计算模型与数值方法

本文选取的计算模型是参照文献[4,6]中实验模型而建立的一个二维矩形的管道(由于本文重点研究气泡横截面的变形情况,因此将模型简化为二维结构),模型上下对称,计算模型如图1所示。假设入射激波由管道的左侧进入并向右传播,整个管道除了球以外都为空气,其中激波左侧为波后空气,右侧为波前空气(静止),波后物理参数查表得出,球中为He或者SF6气体。计算区域长为8 cm×3 cm,球直径D=3 mm,计算区域足够大,可以保证上下边界不会产生反射来影响计算结果。上下壁面为绝热光滑壁面,入口为压力入口,出口为压力出口,本文不考虑粘性影响,故没有边界剪切层。

本文选用非化学反应的组分输运模型来实现多组分可压缩的流动和湍流的模拟,选择LES模型描述非稳态的湍流场,压力与速度耦合采用Coupled算法求解,控制体压力的求解采用Body-Force-Weighted算法进行计算,连续性方程、动量方程、能量方程采用Third-Order MUSCL格式进行离散,组分方程采用quick格式进行离散。

2 计算结果及分析

2.1 He球形界面的变形和发展情况

He气泡变形与发展情况的数值计算结果和文献[4,6]中的实验结果对比情况如表1所示。从表中可以看出,两者的界面变形与发展情况能够很好的吻合,并且数值计算结果能够很好地反映气泡与空气的混合情况,到气泡发展的后期本文数值计算的运动边界更光滑,结构更清晰,基于图中气泡变形图片的细致对比,可以验证本文的多组分无化学反应数值计算的算法是精确可靠的。当入射激波到达Air/He的上游(左)界面且与气泡相互作用后,上游界面首先开始向气泡里凹,气泡由球形逐渐变为椭圆形,随着时间的推移变形加剧,下游(右)界面基本保持不变,形成类似于月牙形的形状(如表1(a)—(d)中气泡的发展过程)。上游界面不断向下游界面运动,最终与气泡的下游界面相碰,整个气泡如同两片肺叶,从上下游界面相碰撞处分开,形成上下两个相互对称的气泡(如表1(e)—(g)中气泡的发展过程)。随着进一步的发展,两个气泡的尾部都回出现涡,此时尾部已经是空气与He的混合物(如表1(g)所示),空气不断被卷入涡中,使得尾部的空气的含量越来越多。随着界面的发展,涡的体积不断增大,整个气泡的长度在不断的增长,即上下游的距离越来越远,上游界面和下游界面之间连接的部分越来越细,直至断裂,分开为上下相互独立两部分,气泡中He含量会越来越少(如表1(g)—(n)中气泡的发展过程)。从数值计算与实验的对比可知,在气泡发展的前后期吻合度都非常高。

注:图中每幅图对应的数字,只代表时间的先后顺序,并不代表准确时间,数值与实验的对比结果不是相同发展时间的对比。

激波与He气泡球形界面相互作用过程形成的涡量云图如图2所示。从图中可以看出,在不稳定性发展的初期(如图2(a)—(c)所示),都是在气泡的表面产生涡量,这是因为气泡表面密度梯度和压力梯度不平行导致的。这也表明了激波与球形界面的相互作用的前期中,斜压机制是导致涡量产生的主要因素。当入射激波与气泡作用后在气泡的内部会产生折射激波以及在气泡的表面会产生绕射激波,由于He声速较大、声阻抗较小,折射激波的速度会较入射激波快,在激波与气泡上游界面接触时,上游界面的表面产生涡,说明入射激波对涡的产生贡献较大,而气泡内部的折射激波贡献较小(如图2(a)所示)。随着激波的进一步运动,射流结构的产生,涡量会向下游界面聚集,在射流撞向下游界面后,在下游界面处产生一对称的涡环(如图2(d)—(g)所示),此时上游的涡量较少,涡主要集中在涡环的下游,而此时在气泡的尾部有“尖钉”结构的产生(如表1(d)—(f)所示),说明“尖钉”结构的产生与涡量有关,在下文会进行定量分析说明“尖钉”结构的产生与涡量的具体关系。随着进一步的发展,上游界面的涡量会越来越少,上下游之间的连线会越来越长越来越细最后断裂,这是因为涡环中的涡向右运动的速度较快的缘故。由此可见涡量的产生和分布对于界面的变形和失稳起着非常重要的作用。

马赫数为1.23与马赫数为1.30的上游界面与下游界面的位移与时间的关系曲线如图3所示。激波与气泡开始作用后,上游界面获得一定的速度向右运动,由于He的密度较小,气泡内的折射激波很快运动到下游,使得下游界面也获得一个较大的速度,很明显在激波与气泡作用的中后期上下游界面的运动速度与时间成线性关系,从图中可以看出,下游界面的运动速度要比上游界面的大,导致涡量主要集中在气泡下游尾部,产生“尖钉”结构。经计算,马赫数为1.23的上游界面的运动速度为117 m/s,下游界面的运动速度为172 m/s,而参考文献[7](具有与文献[4,6]相同的实验模型)中马赫数为1.23的实验结果是,上游界面的运动速度为114 m/s,下游界面的运动速度为151 m/s,实验与数值的下游界面运动速度存在21 m/s的差,氦气实验和计算的结果吻合较差一些,但从界面演化的趋势来看,二者还是一致的。马赫数为1.3的上游界面的运动速度为147 m/s,下游界面的运动速度为205 m/s,通过比较可以发现,马赫数增大,上下游界面的运动速度都会增大,下游界面比上游界面的速度的增加会略微快些。

2.2 SF6球形界面的变形和发展情况

SF6气泡的变形与发展的数值计算结果与文献中的实验结果的对比情况如图4所示。从图4中可看出,气泡的发展过程可分为3个阶段:a)气泡被压缩成椭圆过程。此过程气泡不断被压缩,体积减小,气泡逐渐在x轴方向的半径变得比y轴方向的半径小,整个气泡呈现椭圆型,且气泡的下游界面有向上游界面移动的趋势,使得下游界面的曲率变大,而上游界面并没有多大变化(如图4(a)—(f)显示的气泡发展过程),这主要是由于在气泡内的折射激波运动的速度比气泡外的入射激波慢的原因。b)射流以及尖钉结构的形成。随着界面的进一步的发展,在下游界面会有射流结构的产生,这种结构是在气泡内部的折射激波聚焦所产生的高压驱动作用下成的,随着射流结构的发展,气泡也发生了很大的变形(此时已不是椭圆形状),并且射流会越来越长,但是射流尾部的密度会变得越来越小(如图4(j)—(n)所示),这主要是因为SF6气体与空气的剪切和混合以及气体的扩散作用引起的。在射流结构的发展过程中,会伴随有尖钉结构的出现,尖钉会很快翻转成蘑菇状。c)气体全面混合过程。在这个过程,蘑菇状不断向里翻转,空气不断进入SF6中,导致气体密度越来越小,而射流也不断与空气混合,实验中的射流结构较数值中消失早,原因可能为实际流动为三维,而计算采用两维,必然有定量差异。同时在这个过程会产生很多漩涡,在扩散效应以及漩涡的诱导下,界面得到了进一步的发展,界面不稳定性得到了进一步的加剧,最终两种气体会充分混合(如图4(n)—(p)所示)。

SF6气泡涡量云图如图5所示。从图中可以看出在入射激波与SF6气泡相互作用的初期,在气泡的上游表面会有涡的产生(如图5(a)—(b)所示)。随着入射激波的进一步向右运动,在气泡的全部表面都有涡量的出现,而此时入射激波未到气泡的最右端面,此涡量主要是由已完全经过气泡最右端面的绕射激波引起的,从图5(c)—(d)可以看出,在气泡的上下顶点附近处的涡量最大,这可能是因为在顶点处的密度梯度与压力梯度的不平行度最大,而在气泡的对称轴处却没有涡量,可能是因为此处的密度梯度与压力梯度的平行,这也进一步说明了斜压效应是导致涡量产生的主要原因。随着射流结构的产生,在射流结构的表面有相应的涡量产生,并且涡量会向射流结构的尾部聚集,使得尾部的涡量越来越多(如图5(i)—(n)所示),这是由于射流不断向外运动与空气相互作用引起速度剪切的原因。大涡量的存在会导致射流尾部的速度变大,使得射流结构越来越长。到了后期,射流结构的涡量会越来越少,涡量主要会向下游界面的上下端面聚集,且也不仅仅聚集在气泡的表面,此阶段气泡与空气不断的混合,界面会发生大规模的破碎,进一步加剧界面的不稳定性。

马赫数为1.23与马赫数为1.30的上游界面与射流结构的位移与时间的关系曲线图如图6所示。激波与气泡相互作用的初期下游界面基本不动,直到绕射激波扫过时下游界面才会有略微向上移动,而后形成射流结构,射流结构的位移与时间成线性的关系,射流结构产生后的上游界面的位移也与时间成线性的关系。从图6中可以看出,射流结构的速度要比上游界面的速度大,这主要是由于上文所说的在射流结构尾部的涡量诱导所致,而上游结构的涡量少。经计算,马赫数为1.23的上游界面的变形速度为64 m/s,射流结构的速度为118 m/s,而参考文献[7](具有与文献[4,6]相同的模型)中马赫数为1.23的实验结果中,上游界面的变形速度为57.5 m/s,射流结构的速度为99.4 m/s,实验与数值的下射流结构的运动速度存在18.6 m/s的差,也就是说SF6气泡实验和计算的结果后期射流结构吻合较差一些,但从界面演化的趋势来看,二者还是一致的。马赫数为1.3的上游界面的变形速度为78 m/s,射流结构的速度为155 m/s,速度差为77 m/s马赫数的增加,上游界面与射流结构的速度都会增加,并且射流结构会增加的更快。

3 结 论

本文选用非化学反应的组分输运模型和LES模型来实现多组分可压缩的流动及湍流的模拟,对激波与氦气和SF6气泡相互作用的发展过程进行了数值研究,总结了平面激波作用下不同密度的球形界面变形和发展的物理规律。结论如下:

a) 在He气泡变形及发展的过程中,气泡的上游界面变形较大,并且会形成射流和涡环结构。通过涡量云图发现,由于涡的自旋效应使得下游涡环中所聚集的涡量会越来越多,上下游之间连线会越来越长、越来越细,最终两者分离开来。由此可知,涡量的产生和分布对He与环境中的空气混合、界面的变形和失稳都有重要的影响。同时通过位移时间曲线图,发现了涡量的分布导致上下游界面速度不同,从而对气泡变形产生影响。

b) 在SF6气泡变形及发展的过程中,SF6气泡的下游界面变形较大,并且会形成射流及尖钉结构。通过涡量云图发现,涡量的大小和分布对与环境中的空气混合、尖钉结构发展为蘑菇状以及最后的失稳都有重要的影响。同时通过位移时间曲线图,发现了涡量的分布导致上游界面速度小于射流结构的速度,从而推动着射流的发展。

[1] Richtmyer R D. Taylor instability in shock acceleration of compressible fluids[J]. Commun Pure Appl Math, 1960, 13(2): 297-319.

[2] Zhang Q, Sohn S I. An analytical nonlinear theory of Richtmyer-Meshkov instability[J]. Physics Letters A, 1996, 212(3): 149-155.

[3] Meshkov E E. Instability of the interface of two gases accelerated by a shock wave[J]. Fluid Dynamics, 1969, 4(5): 101-104.

[4] Hass J F, Sturtevant B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987, 181: 41-76.

[5] Layers G, Jourdan G, Houas L. Distortion of a spherical gaseous interface accelerated by a plan shock wave[J]. Phys Rev Lett, 2003, 91: 174502.

[6] 施红辉, 卓启威. Richtmyer-Meshkov不稳定性流体混合区发展的实验研究[J]. 力学学报, 2007, 23(3): 417-421.

[7] 司 廷, 翟志刚, 罗喜胜, 等. 平面激波以及反射激波作用下球形气体界面演化的实验研究[J]. 气体物理, 2011, 6(3): 76-83.

[8] 胡光初, 乐嘉陵, 曹文祥, 等. TVD数值技术在弱电离化学非平衡斜激波反射流场计算中的应用[J]. 计算物理, 1994, 11(2): 195-202.

[9] 宗文刚, 邓小刚, 张涵信. 双重加权实质无波动激波捕捉格式的改进和应用[J]. 空气动力学学报, 2003, 21(4): 399-407.

[10] 李平, 柏劲松, 王 涛, 等. 激波作用下气柱不稳定性发展诱发湍流大涡数值模拟[J]. 中国科学: G 辑, 2009, 39(9): 1241-1247.

[11] 翟金明, 王世合, 李永池, 等. 基于主-次激波在变截面管中传播的数值模拟[J]. 计算力学学报, 2010, 27(5): 925-929.

[12] 吴 宇, 施红辉, 王 超, 等. 液柱在激波冲击下RM不稳定性和破裂过程的数值计算[J]. 浙江理工大学学报, 2014, 31(5): 502-506.

(责任编辑： 康 锋)

Numerical Study of Interaction Process of Planar Shock Waveand Bubble Interface

YE Bin, WANG Chao, SHI Hong-hui, ZHANG Ke, XIONG Hong-ping

(School of Mechanical Engineering and Automation, Zhejiang Sci-Tech University,Hangzhou 310018, China)

The numerical simulation was carried out for interaction process among planer shock wave and He or SF6bubbles in air. The large eddy simulation (LES) method and the Third-Order MUSCL algorithm with third-order accuracy were adopted to simulate bubble development and deformation process. Then, the results were compared with experimental results in literatures. The results show that no matter whether the shock wave passes through low-density interface or high-density interface, the jet and eddy structures will be generated, but the generating mechanisms are different. For the low-density interface, the buddle will evolve into two symmetrical smaller bubbles, and generate two separate vortices and jet flow which lasts in a short time. For the high-density interface, the bubble will evolve into the “spike” structure.

shock wave; numerical calculation; large eddy simulation (LES); bubble and spike; jet flow

1673- 3851 (2015) 05- 0682- 06

2014-10-30

国家自然科学基金项目(10802077);浙江理工大学流体工程技术创新团队资助项目(11132932611309)

叶 斌(1989-),男,江西上饶人,硕士研究生,主要从事汽车制造技术方面的研究。

王 超,E-mail:wangchao@zstu.edu.cn

O382

A