刘 娜，陈艺冰

（1.北京应用物理与计算数学研究所，北京100088；2.中国工程物理研究院高性能数值模拟软件中心，北京100088；3.北京应用物理与计算数学研究所计算物理实验室，北京100088）

多介质流体力学计算的谱体积方法＊

刘 娜1，2，3，陈艺冰1，3

（1.北京应用物理与计算数学研究所，北京100088；2.中国工程物理研究院高性能数值模拟软件中心，北京100088；3.北京应用物理与计算数学研究所计算物理实验室，北京100088）

针对高维及多物理耦合计算耗费大等困难，设计适合多介质流动模拟的模板紧致、易于并行、高阶精度、计算耗费小的谱体积方法。该方法是求解双曲型守恒率谱体积方法的直接推广，针对多介质流动物质界面捕捉的困难，利用拟守恒格式的思想避免物质界面处的非物理振荡。数值模拟结果表明，本方法具有高阶精度、高分辨率，且节约计算量，并且可以有效避免物质界面处非物理振荡。

流体力学；谱体积方法；高精度；刚性气体；多介质；模板紧致

多介质流体力学的模型可以用扩展欧拉方程组描述。早期，使用守恒格式直接求解，但在接触间断（物质界面）附近，会产生非物理的压力和速度振荡。随后，很多学者对此进行分析，并提出各种补救措施［1－2］，其中R.Abgrall等提出的拟守恒型格式，由于其简单、健壮得到了广泛应用。此后，K.M.Shyue［3－4］、Y.B.Chen等［5］进一步发展了这个方法。这些工作通常采用MUSCL限制器，精度不超过二阶。

当前，实际工程问题中对精细化设计提出越来越高的要求，不少更高阶（高于2阶）精度的格式被应用于求解多介质问题。其中，E.Johnsen等［6］将WENO格式推广于求解扩展欧拉方程组，取得成功。他的算法极大提高了对物质界面的分辨率，但WENO格式进行重构时，需要利用目标单元及其周围单元的信息，随着格式阶数的增大，需要不断扩展模板，在高维情形复杂度将急剧增加。J.Zhu等［7］采用DG格式来求解这个问题。与WENO相比，DG格式仅采用局部模板即可获得高精度。实际上，DG格式是近年来广受关注的高精度紧致格式的典型代表之一。除了DG方法之外，谱体积（SV）格式［8－9］也是一类成功的高精度紧致格式。它是谱方法和有限体积方法的结合体，其主要思想是根据细分的控制体的单元均值，重构出谱体积元中的近似解。本质上，DG方法和SV方法都是采用目标单元内部信息来获得网格单元内解的逼近，因此模板都较为紧致，且在并行计算时不需要通信周围单元的信息，是一种高并行度的格式。与DG方法相比，SV方法的一个优点在于节省了表面积分和体积分的计算，并且根据Z.J.Wang［8］的分析，SV方法对CFL数的限制也更低。

本文中，针对高维及多物理耦合计算耗费大等困难，设计适合多介质流动模拟的谱体积方法。

1 控制方程

非定常可压缩流体可以用欧拉方程来描述，以一维为例可以记为：W／t＋F（W）／x＝0，其中W＝（ρ，ρu，ρE）T，F（W）＝（ρu，ρu2＋p，（ρE＋p）u）T。刚性气体状态方程为p＝（γ－1）ρe－γp，其中γ为比热比，p为参考压力。由于流体流动时状态量γ随着流动保持不变，因此满足［10］：γ／t＋uγ／x＝0。定义质量比分Y＝ρ1／（ρ1＋ρ2），建立γ与Y的关系：则质量比分同样满足：

式中：W＝（ρ，ρu，ρE，ρY）T，F（W）＝（ρu，ρu2＋p，（ρE＋p）u，ρuY）T。

2 守恒谱体积方法

直接将谱体积方法推广应用于多介质Euler方程组（2）。下面给出格式设计的基本框架：将计算区域Ω划分为N个子区间｛Si＝［xi－1／2，xi＋1／2］Ω，i∈Z｝，每个Si称为谱体积元，hi＝xi＋1／2－xi－1／2表示其区间长度。假设要设计一个k阶精度格式，则将每个谱体积元Si进一步分化成k个控制体，即Ci，j＝（xi，j－1／2，xi，j＋1／2），j＝1，…，k。在每个控制体上定义控制体单元平均值tn）dx，j＝1，…，k。将控制方程（2）在Ci，j上积分，可以得到如下半离散格式：

要完成k阶谱体积格式，还需要3个步骤。

步骤1：重构。利用给定的tn时刻的控制体单元平均值进行重构。在每个谱体积单元Si内部，由k个控制体的单元平均值，重构出在谱体积元Si内的一个至多k－1次多项式Wni（x）。这样得到的重构多项式可能会引起数值振荡，为了避免高阶重构带来的数值振荡，对重构多项式利用基于原始变量的修正TVB minmod斜率限制器的思想进行限制。详细过程可以参见文献［8］。经过限制，在整个计算区域的重构函数具有如下性质：（1）在谱体积单元内部不需要限制的控制体界面处，重构多项式函数是连续的；（2）在谱体积单元界面处或谱体积单元内部限制器起作用的控制体单元界面处，重构多项式函数是间断的。

步骤2：计算数值通量。根据重构多项式的性质，在谱体积单元内部状态连续的控制体单元界面处，直接使用连续的数值通量，而在谱体积单元界面处以及谱体积单元内部经过限制后的控制体单元界面处，由于两侧的状态值不同，使用间断数值通量。常用的数值通量计算方法，如Lax－Friedrichs、HLLC、Roe、BGK通量等都可以使用，本文中选择最简便的Lax－Friedrichs数值通量。

步骤3：时间方向离散。时间方向可以采用4阶Runge－Kutta方法离散。

3 稳定性分析

根据第2节谱体积方法在多介质守恒扩展欧拉方程组上的直接推广计算多介质问题，会产生物质界面处压力、速度的非物理振荡。在第5.2节将显示一维多介质运动界面问题的计算结果。R.Abgrall等［1－2］指出，这些非物理振荡是由于对组分方程的守恒离散造成的。下面以Lax－Friedrichs通量为例进行分析，为了简便起见，时间方向考虑Euler方法离散，则质量、动量、能量方程的离散可以表示为：

若n时刻速度u＝u0和压力p＝p0为常数，则根据（5）结合（4），可得：，即在n＋1时刻速度可以保持常数：进一步，将状态方程代入（6），则n＋1时刻满足

4 拟守恒谱体积方法

根据第3节的稳定性分析，对第2节中的数值格式做一些修正。将质量比分方程（1）写成如下形式Y／t＋f／x－Yg／x＝0，其中f＝uY，g＝u。质量比分方程的半离散格式修正为：

在物质界面处，速度u保持常数，则离散格式（8）满足压力无振荡的必要条件（7），且离散格式与质量比分方程（1）相容。

5 数值模拟

数值模拟用以检验格式的精度和性能，计算中一维算例CFL数取1.0，二维取0.5。

5.1 精度测试

在区域［0，2］上有两种介质组成的正弦波：

周期边界条件。表1中给出t＝1.0时各阶格式精度测试结果，结果显示2～5阶格式都可以达到相应收敛率。

表1 格式的数值精度Table 1 Numerical accuracy of present schemes

5.2 一维多介质运动界面问题

在区域［0，1］上有两种介质流体，以相同的速度匀速向右移动，界面匀速向右移动，初值条件为：

图1显示了采用守恒谱体积格式100个谱体积单元在t＝0.1时的压力和速度，可以看出，在物质界面处出现速度与压力的数值振荡，与第3节中的分析结果一致。为方便看清图像，图中只显示了奇数阶结果。图2显示了格式修正后的计算结果，结果显示，此时在物质界面处不会产生速度和压力的数值振荡，物质界面附近密度和比热比的局部放大图显示高阶格式具有更高的分辨率。

图1 守恒谱体积格式的一维多介质运动界面问题Fig.1 One－dimensional moving interface problem for conservative spectral volume scheme

图2 拟守恒谱体积格式的一维多介质运动界面问题Fig.2 One－dimensional moving interface problem for quasi－conservative spectral volume scheme

5.3 高压力比气液激波管问题

在区域［－0.2，1］上求解初值问题：

图3显式了t＝0.000 2时刻使用720个谱体积单元1、3、5阶拟守恒谱体积格式计算的结果，结果显示，此时在物质界面处不会产生速度和压力的数值振荡。图4显示了密度及其在物质界面和激波附近的局部放大图，结果表明高阶格式具有更高分辨率。

图3 拟守恒谱体积格式的高压力比气液激波管问题Fig.3 High pressure ratio gas－liquid shock tube problem for quasi－conservative spectral volume scheme

图4 拟守恒谱体积格式的高压力比气液激波管问题Fig.4 High pressure ratio gas－liquid shock tube problem for quasi－conservative spectral volume scheme

5.4 三点问题

初值条件为：

初始空气（轻介质）位于水（重介质）上方，两种介质流体保持静止，静止高密度高压气体位于其右侧分别与低压空气和水产生作用，分别产生向左运动的稀疏波、向右移动的激波，以及推动物质界面向右移动。由于右侧上方的空气较轻，因此右侧上方的激波和物质界面向右移动的速度更快，于是在3种状态的流体的交汇点产生漩涡。形成漩涡后，上层空气与下层水的交界面发生弯曲，并与激波发生作用，形成马赫杆。图5显示了t＝1.0时使用600×200个网格2～5阶格式计算的结果，可以看到，各阶格式都能计算出此时基本的波结构，并且随着格式精度的增加能够捕捉到更精细的结构。

图5 拟守恒谱体积格式的三点问题密度等值线图Fig.5 Density contours of triple point problem for quasi－conservative spectral volume scheme

6 结 论

以多介质拟守恒格式为基础，通过高精度谱体积重构，构造了一类求解多介质问题的高阶精度拟守恒谱体积格式，有效提高了求解多介质问题的精度，并且继承谱体积格式模板紧致、易于并行等优点。数值结果表明，新格式求解多介质问题具有高阶精度，且模板紧致，并且不会在物质界面处产生非物理振荡。

参考文献：

［1］Abgrall R.How to prevent pressure oscillations in multicomponent flow calculation：A quasi conservative approach［J］.Journal of Computational Physics，1996，125（1）：150－160.

［2］Abgrall R，Karni S.Computations of compressible multifluids［J］.Journal of Computational Physics，2001，169（2）：594－623.

［3］Shyue K M.An effcient shock－capturing algorithm for compressible multicomponent problems［J］.Journal of Computational Physics，1998，142（1）：208－242.

［4］Shyue K M.A fluid－mixture type algorithm for compressible multicomponent flow with Mie－Gruneisen equation of state［J］.Journal of Computational Physics，2001，171（2）：678－707.

［5］Chen Y B，Jiang S.A non－oscillatory kinetic scheme for multicomponent flows with the equation of state for a stiffned gas［J］.Journal of Computational Mathematics，2011，29（6）：661－683.

［6］Johnsen E，Colonius T.Implementation of WENO schemes in compressible multicomponent flow problems［J］.Journal of Computational Physics，2006，219（2）：715－732.

［7］Zhu J，Qiu J X，Liu T G，et al.RKDG methods with WENO type limiters and conservative interfacial procedure for one－dimensional compressible multi－medium flow simulations［J］.Applied Numerical Mathematics，2011，61（4）：554－580.

［8］Wang Z J.Spectral（finite）volume method for conservation laws on unstructured grids：Basic formulation［J］.Journal of Computational Physics，2002，178（1）：210－251.

［9］Wang Z J，Liu Y.Spectral（finite）volume method for conservation laws on unstructured gridsⅢ：One dimensional systems and partition optimization［J］.Journal of Scientific Computing，2004，20（1）：137－157.

［10］Karni S.Multicomponent flow calculations by a consistent primitive algorithm［J］.Journal of Computational Physics，1994，112（1）：31－43.

High order spectral volume method for multi－component flows

Liu Na1，2，3，Chen Yibing1，3
（1.Institute of Applied Physics and Computational Mathematics，Beijing100088，China；2.CAEP Software Center for High Performance Numerical Simulation，Beijing100088，China；3.Laboratory of Computational Physics，Institute of Applied Physics and Computational Mathematics，Beijing100088，China）

Numerical simulation of multi－material flows has been an important issues in CFD，and most CFD production codes used for multi－material flow simulation is of either first or second order accuracy，too inefficient and costly with its grid refinement for high accuracy required problems.In this paper，a high－order，efficient，compact method，called the spectral volume method，was developed for the simulation of the multi－material flow as an extension of the spectral volume method for the conservation laws.It has been pointed out that the conservative spectral volume method for the multi－material flow will cause oscillation，and the reason for this has been analyzed.So the idea of quasi－conservative scheme was borrowed to prevent the spurious oscillations in the vicinity of a material contact discontinuity.Several numerical experiments proved that there is no oscillation near the material interface and the result also demonstrates the accuracy，the efficiency and the high performance of the scheme for the multi－material flow simulation.

fluid mechanics；spectral volume method；high order；stiffened gas；multi－component flows；compact stencil

O357.4国标学科代码：1302511

A

10.11883／1001－1455（2017）01－0114－06

（责任编辑 丁 峰）

2015－05－11；

2015－09－20

国家自然科学基金项目（11101047，11501043，11671050，91430218，U1630247）；

中国工程物理研究院科学技术发展基金项目（2013A0202011）

刘 娜（1986— ），女，博士，助理研究员，liu＿na＠iapcm.ac.cn。