真实虚拟流方法在多介质可压缩流动模拟中的应用＊

陈 宏，朱卫兵，张小彬，孙润鹏，郭金鑫

（1.哈尔滨工程大学航天与建筑工程学院，黑龙江 哈尔滨 150001；2.中国航天科工集团三十一研究所高超声速冲压发动机技术重点实验室，北京 100074）

多介质可压缩流应用范围较广泛，如水下爆炸、航空发动机及超音速冲压发动机燃料的混合和燃烧效率的提高、惯性约束热核聚变的控制等［1－2］，而多介质可压缩流中激波及大密度比介质界面附近的计算精度问题一直是计算流体力学研究中的难点和热点。由于在界面两侧流体的状态方程不同，流体密度在界面处出现阶跃，采用传统的单相流方式处理此类问题时，界面处会出现非物理振荡。当界面变化剧烈时，会使计算无法持续进行。因此，界面附近流场的处理就显得尤为重要。

为了解决多介质流计算中产生的振荡现象，R.P.Fedkiw等［3］提出了采用GFM处理多介质流动问题。该方法在抑制界面附近非物理振荡方面取得了较好效果，但是在解决强间断和大密度比问题（如气－水界面问题）时会引起虚假的物理解。R.Caiden等［4］用GFM处理气－水界面问题，通过对水介质中的速度进行插值和对空气侧的压力进行插值定义了ghost流体的状态。该方法尽管处理气－水界面表现优异，但不适合用于求解气－气界面。T.G.Liu等［5］详细分析了原始GFM的缺陷，认为其处理波系与界面作用时不够精确，在此基础上他们提出了修正虚拟流法（modified ghost fluid method，MGFM），利用特征关系和激波跳跃条件求解当地Riemann解问题，进而更新对应的ghost点。C.M.Wang等［6］对MGFM进行了改进，提出了新的处理方式，即RGFM。该方法对靠近界面的真实流体状态采用Riemann解进行了修正，从而减小了计算误差。在此基础上，C.M.Wang等［6］初步计算了气－水界面变化问题。陈荣三等［7］通过在界面处构造Riemann问题对原始GFM进行了改进，对一维强激波与气－气、气－液界面的相互作用问题以及射流问题进行了计算。本文中，将level set方法［8］、HLLC方法［9］结合GFM及RGFM运用于气－气、气－水界面两相流的数值计算，并比较用2种ghost方法得到的计算结果。在气－气界面两相流中，用RGFM捕捉到的波系和界面分布更准确，同时采用RGFM处理气－水界面复杂变化时能得到较好的结果，而采用GFM处理气－水界面问题时遇到了困难，程序无法进行下去。

1 计算方法

1.1 基本方程

守恒方程为二维Euler方程，表示如下：

式中：x和y分别为横向和纵向坐标，ρ和p分别为密度和压力，u和v为分别为x和y方向上的速度，体积总能量

由于HLLC格式在HLL（Harten－Lax－Van Leer）格式的基础上考虑了接触间断和剪切波问题，因此，本文中采用HLLC格式离散数值通量

式中：Sl和Sr分别为左、右波速，Sm为中间波速，详见文献［9］。中间状态变量可由Rankine－Hugonoit关系式得到：

式中：k＝l，r；ρ＊k为左右2束波之间的物质密度；q≡unx＋vny。

中间波速和最大、最小波速分别估计为：

式中：c为声速，u＊、v＊为Roe平均变量。

1.2 level set方法

level set方程为基本思想为：构造函数φ（x，t），使得在任意时刻，运动界面Γ（t）恰是φ（x，t）的零等值面，一般可取φ（x，t）为点x到界面Γ（0）的符号距离。φ（x，t）以适当的速度移动，在任意时刻，只要求出φ（x，t）的值，物质界面的位置也就确定了。

level set方法通过距离函数φ确定物质界面的位置。x到Γ（0）的符号距离

level set方程的数值求解：方程（12）空间项采用五阶 WENO格式［10］求解，时间项采用三阶TVDRunge－Kutta方法［11］求解。level set重构方程表示为

1.3 真实虚拟流法（RGFM）

与原始的GFM［3］相比，真实虚拟流法（RGFM）允许修改界面处真实流体的压力、速度和密度，使计算结果更精确。通过求解Riemann解得出界面处的参数，并将其赋予虚拟网格节点，如图1所示。假设物质界面位于节点i和i＋1之间，首先在界面处定义Riemann问题，界面左侧状态Ul＝Ui，右侧状态Ur＝Ui＋1，然后在此基础上求解Riemann解（本文在此处采用精确Riemann解［9］），可以得到界面处的速度ui和压力pi以及界面左侧的密度ρi，l和右侧的密度ρi，r，将所求得的Riemann解重新赋值（ui、pi和ρi，l）给ghost点（i＋1和i＋2），同时将（ui、pi和ρi，l）赋给真实流体中紧邻界面的节点i。

图1 真实虚拟流方法Fig.1Real ghost fluid method

2 数值模拟

2.1 一维算例—Sod激波管问题

初始时刻压力比p2／p1＝10的Sod激波管问题初始条件为：当x＜0.6时，ρ2＝1.000，p2＝10.0，u2＝0.0，γ＝1.4；当x＞0.6时，ρ1＝0.125，p1＝1.0，u1＝0.0，γ＝1.4。计算结果如图2所示，计算时间为0.035s，图中给出了用GFM和RGFM得到的密度、速度和压力的数值解。从图2中可以看出，与解析解相比，RGFM捕捉到的激波、膨胀波及界面的分布明显好于GFM的计算结果。

压力比p2／p1＝1 000时，由于初始时刻流场中存在较大的间断（压力和密度），随着时间的延长，在流场中产生了向右运动的强激波和接触间断，同时还产生了向左运动的膨胀波。由于初期接触间断的速度较高，其紧随激波向右运动。图3为t＝0.005s时计算结果与解析解的比较。由图3可以看出，在这种极端的情况下采用RGFM得到的计算结果与解析解仍然符合较好，而采用GFM得到的结果则在界面处出现了振荡现象，得到了非物理解。

图2 p2／p1＝10时的Sod激波管问题Fig.2Sod’s shock tube problem at p2／p1＝10

图3 p2／p1＝1 000时的Sod激波管问题Fig.3Sod’s shock tube problem at p2／p1＝1 000

2.2 二维算例

2.2.1 激波打柱形氦气泡

激波与气泡之间的相互作用应用背景广泛，从天体物理到人体组织的空泡损伤（如对肾结石的震波碎石等）。J.L.Hass等［13］实验研究了激波与圆柱形氦气泡的相互作用。本算例即以该实验为基础，研究了激波由重气体进入轻气体的物理过程，本算例也在该实验的基础上对马赫数为1.22的激波在空气中与氦气泡的相互作用过程进行了数值模拟。由于计算域对称，只计算了上半区域。图4为计算域示意图，对几何尺寸进行了量纲一化，以下各算例与此相同。气泡内区域为Ⅰ，气泡外激波前区为Ⅱ，激波后区域为Ⅲ，入射激波速度为1.22Ma，从右向左运动。计算网格为1 600×200，其量纲一初始条件［3］为：气泡内，ρ＝0.138 0，p＝1.0000，u＝0.000，v＝0.0，γ＝1.67；激波前，ρ＝1.0000，p＝1.000，u＝0.000，v＝0.0，γ＝1.4；当x＞225时，ρ＝1.376 4，p＝1.569 8，u＝－0.394，v＝0.0，γ＝1.4。

图5给出了对该物理过程采用RGFM得到的计算结果和实验图像［13］，为了能清晰地观察到波系的变化，采用纹影技术［14］进行了后处理。入射激波在空气中从右向左运动，首先打到气泡右侧界面上，右界面受此作用开始向左运动，同时在右界面上形成Richtmyer－Meshkov（RM）不稳定性。入射激波打到界面上产生的反射波为膨胀波，而折射波为激波；由于激波在氦气泡中的传播速度比在空气中的快，所以可以看到折射激波在气泡中形成弧形激波走在了入射激波的前面，这样就在气泡上下界面形成了速度差，使得气泡上下界面开始出现Kelvin－Helmholtz（KH）不稳定性，如图5（a）所示。

当氦气泡中的弧形折射激波遇到气泡左界面时，同样会发生反射和折射，反射波和折射波均为激波，并且在左界面也产生了RM不稳定性。由图5（b）可以看出，此时反射波较弱，在实验纹影图中已经很难观察到；而折射波以圆弧形从氦气泡中进入空气，在图中可以看到圆弧形激波遇到激波管上下壁面发生反射，且该反射波再次与气泡的上下界面作用，形成新的RM和KH不稳定性；此时气泡在波系的作用下呈半圆形。

由图5（c）和（d）可以看出，气泡右侧界面受气流和不稳定性影响向内塌陷，并逐渐形成丁状结构。随着丁状结构的发展，最终导致气泡断裂，并在断裂处形成射流。

图4 激波管内激波打气泡示意图Fig.4Sketch of the dimensionless computational domain for shock－bubble interaction in a shock tube

图5 采用RGFM得到的计算结果与实验图像［13］的比较Fig.5Comparison between computational results by RFGM and experimental images［13］

图6为激波打气泡中后期气泡断裂后的形态图像，图中同时给出了对应时刻采用GFM得到的计算结果。从图6中可以看出，采用RGFM得到的计算结果与实验结果符合更好，而采用GFM得到的计算结果误差稍大。目前关于激波与氦气泡相互作用问题的数值研究中，大部分文献只给出了气泡变化的初期形态，很少给出气泡大变形甚至断裂后的形态。本文算法给出了整个气泡变化过程的图像，而且均与实验结果［13］符合较好，验证了本文算法的适用性。

图6 用不同方法得到的激波打氦气泡中后期气泡断裂后的形态Fig.6Shapes of helium bubble after fracture in the midanaphase of the interaction between shock wave and helium bubble by different methods

2.2.2 水下激波打气泡

本算例主要分析水下一入射激波扫过圆柱形气泡后的物理变化过程，计算域如图7所示，各边界均为自由边界。初始条件为：气泡内，ρ＝1.00，p＝1.0，u＝0.0000，v＝1.0，γ＝1.4；激波后，ρ＝1 220.75，p＝9 120.0，u＝1.284 1，v＝0.0，γ＝7.0；其余为，ρ＝1 000.0，p＝1.00，u＝0.0000，v＝0.0，γ＝7.0。

采用GFM处理气－水界面时，程序不能顺利计算下去，因此本部分只给出了采用RGFM得到的计算结果。图8给出了网格数为840×720的情况下不同时刻气泡界面的位置，从图中可以看出气泡从变形到破碎的变化过程。

图7 计算域Fig.7Computational domain

图8 不同时刻的界面分布Fig.8Distribution of the interface at different times

图9给出了不同时刻流场密度变化的纹影图。激波在水中从左向右运动扫过柱形气泡物质界面，在空气侧产生一道折射激波，而在水中产生反射膨胀波，如图9（a）所示，由于界面两侧声阻不同，且差异较大，气泡内折射激波速度滞后于入射激波，气泡受到波系作用而产生变形，呈现出“肾”形。随着时间的延长，气泡内折射激波在气泡右侧界再次发生反射和折射，同时在气泡凹陷处形成水射流，如图9（b）所示。图9（c）～（d）中水射流已经将气泡截断，并产生了冲击波。受该冲击波的影响，断裂的小气泡中形成了二次水射流现象。该算例验证了本文算法处理多介质流中激波与大密度比界面相互作用问题的适用性。

图9 不同时刻的密度纹影Fig.9Density schlieren images at different times

3 结 论

（1）将RGFM引入到多介质可压缩流的数值模拟中，研究激波与大密度比介质界面的相互作用问题，实现了界面不稳定的精确模拟。研究和开发了适用于可压缩多介质流体界面（气－气界面、气－水界面）不稳定性发展演化的计算方法和程序。

（2）一维和二维算例表明，RGFM在处理气－气界面时优于原始GFM，能得到更精确的计算结果。

（3）在水下激波打气泡的算例中，采用RGFM成功的捕捉到了大密度比（1 000∶1）界面从变形到破碎的变化过程，证明该方法可以处理大密度比多介质流动中的界面大变形问题。

［1］Brouillette M.The Richtmyer－Meshkov instability［J］.Annual Review of Fluid Mechanics，2002，34：445－468.

［2］Waitz I A，Marble F E，Zukoski E E.Investigation of a contoured wall injector for hyper－velocity mixing augmentation［J］.AIAA Journal，1993，31（6）：1014－1021.

［3］Fedkiw R P，Aslam T，Merriman B，et al.A non－oscillatory Eulerian approach to interfaces in multimaterial flows：The ghost fluid method［J］.Journal of Computational Physics，1999，152（2）：457－492.

［4］Caiden R，Fedkiw R P，Anderson C.A numerical method for two－phase flow consisting of separate compressible and incompressible regions［J］.Journal of Computational Physics，2001，166（1）：1－27.

［5］Liu T G，Khoo B C，Yeo K S.Ghost fluid method for strong shock impacting on material interface［J］.Journal of Computational Physics，2003，190（2）：651－681.

［6］Wang C W，Liu T G，Khoo B C.A real ghost fluid method for the simulation of multimedium compressible flow［J］.SIAM Journal on Scientific Computing，2006，28（1）：278－302.

［7］陈荣三，盛万成.用改进的ghost fluid方法模拟强激波与界面的相互作用［J］.爆炸与冲击，2008，28（2）：144－148.

Chen Rong－san，Sheng Wan－cheng.A modified ghost fluid method for strong shock wave impacting on material interface［J］.Explosion and Shock Waves，2008，28（2）：144－148.

［8］Osher S，Sethian J A.Fronts propagating with curvature－dependent speed：Algorithms based on Hamilton－Jacobi formulations［J］.Journal of Computational Physics，1988，79（1）：12－49.

［9］Toro E F.Riemann solvers and numerical methods for fluid dynamics：A pratical introduction［M］.2nd ed.Springer，1999：115－162.

［10］Jiang G S，Shu C W.Efficient implementation of weighted ENO schemes［J］.Journal of Computational Physics，1996，126（1）：202－228.

［11］Chen S，Merriman B，Osher S，et al.A simple level set method for solving Stefan problem［J］.Journal of Computational Physics，1997，135（1）：8－29.

［12］Tang H S，Huang D.A second－order accurate capturing scheme for 1Dinviscid flows of gas and water with vacuum zones［J］.Journal of Computational Physics，1996，128（2）：301－318.

［13］Haas J F，Sturtevant B.Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities［J］.Journal of Fluid Mechanics，1987，181（1）：41－76.

［14］Quirk J J，Karni S.On the dynamics of a shock－bubble interaction［J］.Journal of Fluid Mechanics，1996，318：129－163.