斜激波与物质界面相互作用的数值模拟

2023年12月20日发(作者：沙雕幼儿园数学试卷)

CNIC-01770

IAPCM-0043

斜激波与物质界面相互作用的数值模拟

唐维军王晨星

(北京应用物理与计算数学研究所，北京，100088)

摘要

利用多介质PPM方法研究了斜激波与物质交界面的相互作用。

采用与体积分数耦合的Euler方程组作为计算模型，用双波近似来

求解一般刚性气体状态方程Riemarm问题，通过体积分数的计算来

获得界面的位置，在整个流场采用统一的高阶PPM格式进行计算。

文中对斜激波与不同物质界面相互作用进行了数值模拟，并给出了

交界面上由于斜压效应产生的涡列的演化过程，特别是强斜激波与

不同物质界面的相互作用情况。

关键词：斜激波PPM方法

75

Numerical Simulation of Oblique Shock Wave and Interface

(In

TANG Weijun

Chinese)

WANG Chenxing

(Institute of Applied Physics and Computational Mathematics, Beijing, 100088)

ABSTRACT

The results of numerical computations of the refraction of a plane shock

wave at the different interfaces are presented. The unsteady，two-dimensional，

compressible, Euler equations numerically assuming stiffened gas equation of

state are solved. The applied numerical method is PPM (Piecewise Parabolic

Method) method for multi-component. The Riemann problem using two-shock

approximate are solved. By computing the volume-fraction the interface is ob-tained ? and the high-order PPM method is used in the whole computational

field. Applying this method，some cases of the inter-action when the shock

wave meets the interfaces of the different fluids are computed. The evolvement

of the interfaces because of the baroclinic effects on the interface is presented，

especially the case of the strong shock striking at the interfaces.

Key words： Oblique shock wave，PPM

引言

当斜激波与多种物质分界面相互作用时，会出现Richtmyer-Meshkov不稳定性现

象[1]，以及由于斜压性导致产生的涡量。研究激波与物质界面相互作用在燃料混合和激光

惯性约束聚变方面有着重要意义。研究平面界面扰动演化的实验如文献[2〜6]所述，研究

柱面和球面界面的实验如文献[7]所述。数值模拟方面的研究如文献[8〜13]所述。激波与

界面相互作用的理论很不完善，如Richtmyei[14〕提出的不稳定性线性理论仅对激波和在穿

过界面很短的时间内有效，StUrtevantraS现实验结果与线性理论相差很大。

Hawley和Zabusky™对弱激波作用下的空气-氦（Air-He)以及空气-氟里昂（Air-Freon)界面演变过程的使用进行了数值模拟，并对两种慢-快界面和快-慢界面生成涡量的

大小和方向进行了比较。严长林[13]利用基于有限体积法的二阶Godimov格式对空气-氟利

昂（Air-Freon)的情形数值模拟了强激波与交界面之间的相互作用过程。Pkone和BoriS[1(>]

的数值模拟只在边际上求解，并且关于涡量生成的公式不易在一般情形应用。Winkler

等[11〕及Chalmers等[8]的数值模拟发现了新的现象，但没有给出充分的解释。上述数值模

拟都是采用二阶精度的格式，且是使用理想气体的状态方程，绝热指数也一样，模拟的是两

种物质。

本文采用三阶多介质PPM方法[15]进行了数值模拟，且对含理想气体的状态方程及刚

性气体状态方程进行数值模拟，界面两边状态方程参数不同，并给出了交界面演化的详细数

值图象。对慢-快界面和快慢界面及弱激波和强激波的差别进行了分析。本文使用La-grange 步再加上remapping步的方法，结合了 Lagrange方法和Euler方法的优点；同时增

加了体积分数的计箅，对界面的描述比上述数值模拟更精确。

1 数值方法

1.1 多介质流体体积方程

多介质可压缩流动中各介质的状态方程一般可以采用一般的刚性气体状态方程，可表

示为

Pip, e) = (r-l)pe-yPo. (1)

假设在一个网格区域内含有m种不同的介质，它们各自所占的体积分数为,且各体

m

积分数满足方程$^>=1。假设已知各组分的物理量广>，(i)M，PG 7(i)，力2，一般地，网

格平均守恒物理量可以写成各组分介质物理量的体积分数加权形式：

m m m -P = ，^ = ，戽=I] /YCO^O^O + 1 Y(,)^(,)M2N

i = l ； = 1 V /

通过接触界面两侧压力平衡和速度一致的假设导出流体体积分数等参数的演化方程，

结合Euler方程组，得出了基于准守恒型控制方程组的计算模型[6’16~1S]：

77

+ 悬,+

+

+

p) = 0

+ p) = 0 (2)

ft(pE) + ^[(^E + p)u~] + ^[(^E + p)v] - 0

I— H u — h v = 0 z = 2，\"•， m

at ox oy

该模型保证了界面两侧能量方程的一致近似并避免了压力振荡，可求解高密度比、含强

激波的多介质可压缩流动，并可方便地推广到含任意数目交界面和多维的情况。

为了更有效地捕捉界面，本文采用Lagrange-Remapping两步计算；对于二维问题，采

用维数分裂的方法，在每个方向上分别进行求解，即在每个方向上都要进行Lagrange步和

remapping步的计算。以广〉，7(\'〉分别表示z和y方向分裂的计算结果，其上标t表示时间

步长，则由i\"时刻的值可以得到：

j^JrH-2 ___ CAT)

A

Y(A/) y(Ax)

A

(U\")

即下一时刻的物理量的值。

—维Lagrange-Remapping型PPM方法可分为如下四步：

(1)物理量的分段抛物插值；

(2)近似Riemann问题求解网格边界时间平均的压力和速度；

(3)Lagrange方程组推进求解；

(4)最后将物理量变回静止的Euler网格上。

以工方向分裂后的一维多介质Lagrange方程组为例。

1.2 Lagrange 步

at

+

am

=0 (3a)

(3b)

(3c)

at

at

+

dm

dm

= 0

学=0

dt

(3d)

方程组中r，u, E, K分别为比容、z方向的速度、比总能、i组分的体积分数，a在平

面、轴对称和球对称时分别取0，1，2。M为质量坐标，r为空间坐标。

物理量的插值与文献[19]所描述的过程一样，这时需要插值的物理量为片，<，p]’

(Y(i)；且无需利用接触间断的检测，因为在拉格朗日步的计算中，接触间断是自动保持的。

为获得网格边界上时间平均的压力和速度，需要在网格边界上求解近似Riemarm问题。利

用双波近似求解Riemarm问题，可得出时间平均的网格物理量,fJ+m ：

uR(p) — uL(p) = 0

_

78

=处-碌y _ = ««

其中’紹㈤=叩+(穿仏餘-1)]，。”\'^。

对Lagrange方程组（3a)〜（3d)进行差分近似，得到下一时刻的网格平均值

xJ+l/2 ~ .^fl/2 ~t\" A^M 汁 1/2

(4a)

(4b)

x =(碏;/2 严一«+1+1/2严

n+l __ (^jtl/2 )° 一

Ti («+l)Am)

)(4c)

ufl = + y (A>+1/2 + AJ-1/2

告 iPi-YZ — Pmn) (4d)

Ef1 = E]+ -^-(A^u^nVj-m

l7Ylj

(y<0-Aj+muj+mp^m) (4e)

)f

1= (Y(i))； i = 1, 2, N-l (4f)

1.3 Remapping 步

考虑鄉=9义为Lagnmge坐标下网格的守恒质量，=(〜,„),¥]为Euler坐标下

网格的守恒质量，AmA^-p/VwAl^w为网格边界通过的流体质量，体积分数最后的变换

公式为

AE^1/2 =

(p^r)j = [Am, -

+ —z^iphm + 丫 PL ) AV;+1/2

(pJ+^AVJ+u, - pJ-m^yj-mWV0(5a)

(5b) ,

(Meui«)； = LujAmj — iuj+m Am/+1/2 — u*-m -in )]/im]

(Eeuler); =

((y(i) )euler), =

A^ _ (AEUi/2 一 厶均*-1/2)]/么4

— ((Y(\') )A1/2AV^1/2 — )A\" AVA/2)]/V?

(5c)

(5d)

(5e)

对于y方向的维数分裂，与X方向的类似，只是对.y方向的速度《进行计算，代替了 ,0：

方向维数分裂时对:C方向的速度《进行计算，而其他物理量的计算与.Z方向进行时是一样

的。

2 数值结果分析

在下面的数值模拟中，激波从计算区域的左侧进人，交界面为与激波成60°的斜面，交

界面的左边为空气（密度为灼），右边为He或Fre0n(密度为作）。设鸬/外，则对于Air-He情形巧=0. 14<1，即激波从重气体区进入轻气体区；对于Air-Freon情形，77=3〉1，即

激波从轻气体区进入重气体区。

计算区域及激波和交界面的初始位置如图1所示，波前为静止气体，波后的状态可以按

照正激波的关系式给出。对于上下壁面，按固壁边界条件给出；上下游边界按流入流出边界

条件给出。

2.1

2.. 1.1

弱激波作用下交界面的演化

在弱激波的情况下，假定激波的马赫数均为1.2。

弱激波作用下Air-He交界面的演化

Air:(o=l. 2, P=l, 013, u=0, 7=1.4, H

计算区域为[0，38,. 75] X[0,6,, 2]，采用1 250X200的网格进行计算。初始条件为：

79

p=0. 167, P=L 013, u = Q,He： 7=

1. 63, P„=0

图2给出了一定时间序列的密度等值线图。激

波和交界面初次相遇的图形如图2(a)所示；激波与

交界面相互作用如图2(b)所示，被人射激波扫过的

界面产生弯曲变形。交界面的位置转换是由于

Richtmyer-Meshkov不稳定性引起的，这是当激波

由重气体向轻气体入射时发生的典型特征。在界面发生转化的过程中，下端点靠近下壁面

处开始有一薄层涡量产生（如图2(c)所示）。涡量是由于激波所产生的压力梯度，与界面所

产生的密度梯度不平行，即斜压项存在，导致涡量产生。除靠近计算区域下边界的部分界面

外，界面的其他部分基本上是一条直线（如图2(d)所示）。在随后的发展过程中（图2(e)〜

(h))，除靠近计算区域下边界的部分界面外，其他界面始终是一条较为光滑的，曲率半径很

大的曲线，界面两边的流体并没有混合。在靠近计算区域下边界的部分界面附近在开始不

久就形成了涡卷。在其后的发展过程中，这个涡卷变大，将周围的涡卷逐渐卷吸进去。但靠

近计箅域下边界的涡卷并不是失稳所致，而是“壁面效应”所导致的。

图2 交界面（Air_-He)的演化

2. 1. 2 弱激波作用下Air-Freon交界面的演化

计算区域为[0，38,.75]乂[0，6.2]，采用625X 100的网格进行计算。初始条件为

Air： ,0=1, P = l, u = 0, y=1.4, Poo=0

Freon:jo=3, P=l, u=0, 7=1.4, PTO=0

图3是按一定时间序列排列的密度等值线图。图3(a)是激波与界面初次相遇的图形。

随后，激波与界面相互作用，且靠近下边界的界面已经有涡卷起（如图3(b)所示）。但除靠

近下边界处的界面外，被扰动的界面几乎呈一条直线。由于界面上的速度大于界面下的速

度，界面被拉长。这进一步加大了初始扰动的幅值。随着时间的发展，激波完全穿过界面，

此时已有涡形成(如图3(c)所示），且大小交替出现。在图3(b)中靠近下边界处卷起的涡已

经离开了下边界，不再附着在下边界处。此时，界面的失稳问题已经变成Kelvin-Helmholtz

失稳问题。随着时间的继续发展，小涡变成大涡，界面继续被拉长（如图3(d)所示）。最后

开始出现涡的并对（如图3(e)，图3(f)所示）。

..

图3交界面（Air-#reon)的演化

^丄丄又o.

f 2.2

2. 2.1

强激波作用下交界面的演化

强激波作用下Air-Freon交界面的演化

激狭马赫数为3，计算区域为[0，124]X[0，6.2]，计箅网格数为4 000X200。初始条

|

件如下：

Air： p=l, P=l, m = 0, y=1.4,

Freon:(o=3, P = 1, m=0, y—1. 4 , =0

强激波的作用过程和弱激波的情况很相似。交界面在斜压效应作用下也会有涡量产

生；涡层都出现了并对的现象。但在强激波的情况下，在激波穿越界面时，界面会发生很大

的变形；交界面的运动速度也要比弱激波的情况快；交界面上的涡层的数目与弱激波情形也

有不同，比弱激波时要少一些。但涡量比弱激波情形要大。涡旋占据的空间宽度也比弱激

波情形要宽。图4给出不同时间序列的界面附近区域的密度等值线图。

图4 交界面（Air-Freon)的演化

2. 2. 2

数：

强激波作用下Air-He交界面的演化

激波马赫数为10，计箅区域为[0,35] X[0，5]，计算网格数采用1 400X200。初始参

Air:1o=1.2, P=l, u=0, y=l. 4, Pre=0

He： p=0. 167，P=l, m=0, 7=1. 63, Po»=0

从图5可见，强激波的作用过程与弱激波的十分相似：当激波与交界面相互作用时，被

人射激波扫过的界面都会产生弯曲变形，这也是

由于Richtmyer-Meshkov不稳定性引起的；在

界面发生转化的过程中，下端点靠近下壁面处都

会有一薄层涡量产生；在随后的发展过程中，除

靠近计算区域下边界的部分界面外,其他界面始

终是较为光滑的；在靠近计算区域下边界的部分

界面附近都形成了涡卷。但是也有一些不同：在

强激波的情况下，界面运动的距离要比弱激波的

长；另外，由于壁面效应更加明显，形成的涡卷也

要比弱激波的大。

2.3 气体与液体交界面的情况

液气界面的情况

在下面的数值模拟中考虑如下的情况：假定

一个马赫数为1,. 95的激波从计算区域的左侧进

人，交界面为与激波成60°的斜面，交界面的左边为液体，右边为气体，此时激波从重流体区

进人轻流体区。边界条件同上。

初始条件采用文献[8]中的数据：

静止气体…=1，P=l, « = 0，y=1.4, Poo=0

静止液体:(o=5，P=l, m = 0, 7=4, P«o = l

激波后液体:(o=7.093, P=10, « = 0.728 8，y=4, P„ = l

计算区域为[0，37.,5]父[0，5]，采用1 500X200的网格进行计算。

2.3.1

(e)

图5 交界面（Air-He)的演化

⑷

从图6可以看出其作用过程与前面气体的情况很相似：当激波与交界面相互作用时，由

于Richtmyer-Meshkov的不稳定性，被人射激波扫过的界面都会产生弯曲变形；靠近下壁

面处都会有一薄层涡量产生，随后都会形成涡卷；除靠近下边界的部分外，其他界面始终是

较为光滑的。但是也有一些不同：由于壁面效应更加明显，形成的涡卷也要比弱激波的大。

2.3.2 气液界面的情况

在下面的数值模拟中，假定一个马赫数为1. 95的激波从计算区域的左侧进入，交界面

为与激波成60°的斜面，交界面的左边为气体，右边为液体，此时激波从轻流体区进入重流

体区。边界条件同上。

计箅区域为[0，37. 5]X[0，5]，计箅网格数为1 500X200。初始条件如下：

静止气体:lo=l，P=l, u=0, y=l. 4, P«,=0

静止液体屮=5，P=l, m=0, y=4,

激波后的状态可按照正激波的关系式计算。

由图7可以看到，这种情况与前面的两种气体分界面的情况也同样很相似。交界面的

角度没有发生反转；交界面在斜压效应的作用下有涡量产生；涡量不断扩散，出现涡的并对

现象。但形成的涡卷要更明显；涡的并对也比前面的弱激波情形要快。

3结论

本文使用多介质体积分数的欧拉方程组的计算模型和经典的单介质PPM方法，结合

了数值模拟斜激波与两种密度不同的流体的交界面作用的情况。从上面的算例中可以看

出，这种方法可以较好地捕捉到物质的交界面，并且可以较好地模拟高密度比、水气两相介

质和带有强激波的可压缩流动。

从计算结果中可以看出，在激波与两种不同密度流体的交界面相互作用时：

(1)当激波从重流体进入轻流体时，交界面将发生反转，这是Rkhtmyer-Meshkov不稳

定性的典型特点；交界面的下端会产生涡；但随着激波强度的增加，壁面效应会更明显，形成

的涡卷也会随之变大。

(2)当激波从轻流体进人重流体时，交界面不会发生反转，但交界面上会产生一系列的

涡，并且这些涡会不断扩张，出现并对现象；当激波强度加大时，并对现象会更明显，从而加

速了两种流体的混合；而且，并对现象会随着两种流体的密度比的加大而愈发明显。

参考文献

1 Yang Yumin，Zhang Qiang„ Small amplitude theory of Richtymer-Meshkov instability, Phys. Fluids，

1994,6(5)： 1856〜1873

2 Andronov V，Bakhrakh S M, Meshkov E： E，Mokhov V N„ Sov„ Phys.，1976，JETP 44： 424〜427

3 Benjamin R E, Trease H E, Shaner J W. Phys. Fluids, 1984, 27： 2390

4 Meshkov E E. Sov„ Fluid Dyn.. 1969,4(5)： 101

5 Sturtevant B„ in Shcok Tubes and Waves，edited by Gionig，H,，(VCH Verlag, Berlin, 1987)，89〜

100

6 Zaitsev, S G, Lazareva E V, Chernukha V V，et al. DokL, Akad,, Nauk SSR.t 1985, 283:94 [Sov,.

Phys.. Dokl.，1985，30:579

7 Haas J-F，Sturtevant B. Interaction of weak shock waves with cylindrical and spherical gas inhomogene-ities. J.. Fluid Mech„，1987，181： 41 〜76

8 Chalmers J W, Hodson S W，Winkler K_H A，Woodward P R，Zabusky N J.. Fluid Dyn. Res,,，1988’

3：392

83

9 Hawley J F，Zabusky N J. Vortex Paradigm for Shock-Accelerated Density-Satisfied Interfaces. Physi-cal Review Letter, 1989, 63 (12)： 1241 〜1244

10 Picone J M, Boris J P„ JL Fluid Mech.,，1988,189： 23

11 Winkler K-H A, Chalmers J W, Hodson S W，Woodward P R, Zabusky N J. Phys,. Today, 1987,40

(10)： 28

12 Youngs D L, Numerical simulation of tuibulent mixing by Rayleigh-Taylor instability. Physical, 1984,

D12： 32〜44

13严长林，孙德军，尹协远，童秉纲,,激波加速的密度分层流体交界面的演化,,空气动力学学报，2001，19

(2):228〜233

14 Richtmyer R D, Taylor instability in shock acceleration of compressible fluids.. Comm.. Pure Appl.

Math.，1960，13:297〜319

15马东军，孙德军，尹协远,，高密度比多介质可压缩流动的PPM方法.计箅物理，2001，18(6):517〜522

16 Saurel R，Abgrall R. A simple method for compressible multicomponent flows:. SI AM J. Sci., Com-put.，1999，21(3)： 1115〜1145

17 Shyue K M.. An efficient shock-captuxing algorithm for compressible multicomponent problems,, J.

Comput. Phys..，1998，142： 208〜242

18柏劲松可压缩多介质流体高精度数值计箅方法和网格自适应技术研究[博士论文]中国工程物理

研究院，2003

19 Collela P，Woodward P. The Piecewise Parabolic Method (PPM) for gas-dynamical simulations. J.

Comput,. Phys..，1984，54:174〜201.,