多介质流体力学的模型可以用扩展欧拉方程组描述。早期，使用守恒格式直接求解，但在接触间断(物质界面)附近，会产生非物理的压力和速度振荡。随后，很多学者对此进行分析，并提出各种补救措施^[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, ρu²+p, (ρE+p)u)^T。刚性气体状态方程为p=(γ-1)ρe-γp_∞，其中γ为比热比，p_∞为参考压力。由于流体流动时状态量γ随着流动保持不变，因此满足^[10]：∂γ/∂t+u∂γ/∂x=0。定义质量比分Y=ρ₁/(ρ₁+ρ₂)，建立γ与Y的关系:$\frac{1}{{\gamma - 1}} = \frac{Y}{{{\gamma _1} - 1}} + \frac{{1 - Y}}{{{\gamma _2} - 1}}$，则质量比分同样满足:

写成守恒律形式∂(ρY)/∂t+∂(ρuY)/∂x=0。这相当于在Euler方程的基础上增加了一个新的守恒量及其守恒律方程，得到了一个完整的多介质Euler方程组，方便起见，仍记为:

式中：W=(ρ, ρu, ρE, ρY)^T，F(W)=(ρu, ρu²+p, (ρE+p)u, ρuY)^T。

2.   守恒谱体积方法

直接将谱体积方法推广应用于多介质Euler方程组(2)。下面给出格式设计的基本框架：将计算区域Ω划分为N个子区间{S_i=[x_i-1/2, x_i+1/2]Ω, i∈Z}，每个S_i称为谱体积元，h_i=x_i+1/2-x_i-1/2表示其区间长度。假设要设计一个k阶精度格式，则将每个谱体积元S_i进一步分化成k个控制体，即C_{i, j}=(x_{i, j-1/2}, x_{i, j+1/2}),j=1, …,k。在每个控制体上定义控制体单元平均值$\overline W _{i, j}^n = \frac{1}{{{h_{i, j}}}}\int_{{C_{i, j}}} {W\left( {x, {t_n}} \right){\rm{d}}x, j = 1, \cdots, k}$。将控制方程(2)在C_{i, j}上积分，可以得到如下半离散格式：

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

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

步骤2：计算数值通量。根据重构多项式的性质，在谱体积单元内部状态连续的控制体单元界面处，直接使用连续的数值通量${{\hat F}_{i, j + 1/2}} = F\left( {{W_{i, j + 1/2}}} \right)$, 而在谱体积单元界面处以及谱体积单元内部经过限制后的控制体单元界面处，由于两侧的状态值不同，使用间断数值通量${{\hat F}_{i, j + 1/2}} = F\left( {W_{_{i, j + 1/2}}^{n - \left( {\bmod } \right)}, W_{_{i, j + 1/2}}^{n + \left( {\bmod } \right)}} \right)$。常用的数值通量计算方法，如Lax-Friedrichs、HLLC、Roe、BGK通量等都可以使用，本文中选择最简便的Lax-Friedrichs数值通量。

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

3.   稳定性分析

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

式中：$\Delta _j^0\rho = \bar \rho _{i, j + 1}^n - \bar \rho _{i, j - 1}^n, \Delta _j^ + \rho = \bar \rho _{i, j + 1}^n - \bar \rho _{i, j}^n, \Delta _j^ - \rho = \bar \rho _{i, j}^n - \bar \rho _{i, j - 1}^n$。

若n时刻速度u=u₀和压力p=p₀为常数，则根据(5)结合(4)，可得：(ρu)_{i, j}ⁿ⁺¹=ρ_{i, j}ⁿ⁺¹u₀，即在n+1时刻速度可以保持常数：u_{i, j}ⁿ⁺¹=u₀。进一步，将状态方程代入(6)，则n+1时刻满足：$\frac{{p_{\mathit{i}, j}^{n + 1}}}{{\left( {\gamma - 1} \right)_{i, j}^{n + 1}}} = \frac{{p_{\mathit{i}, j}^n}}{{\left( {\gamma - 1} \right)_{i, j}^n}} - \frac{{\Delta t}}{h}\left( {\frac{1}{2}\Delta _j^0\frac{{pu}}{{\gamma - 1}} - \frac{\alpha }{2}\left( {\Delta _j^ + \frac{p}{{\gamma - 1}} - \Delta _j^ - \frac{p}{{\gamma - 1}}} \right)} \right)$。

若压力在n+1时刻也保持常数，则需要满足：$\frac{1}{{\left( {\gamma - 1} \right)_{i, j}^{n + 1}}}{\rm{ = }}\frac{1}{{\left( {\gamma - 1} \right)_{i, j}^n}} - \frac{{\Delta t}}{h}\left( {\frac{1}{2}\Delta _j^0\frac{u}{{\gamma - 1}} - \frac{\alpha }{2}\left( {\Delta _j^ + \frac{1}{{\gamma - 1}} - \Delta _j^ - \frac{1}{{\gamma - 1}}} \right)} \right)$，上式等价于：

式(7)可以看做偏微分方程∂Y/∂t+∂(uY)/∂x=0的离散。显然，这与质量比分满足的方程(1)不同，即第2节中的谱体积格式无法满足压力无振荡的必要条件(7)，并因此造成了物质界面处压力、速度的非物理振荡。

4.   拟守恒谱体积方法

根据第3节的稳定性分析，对第2节中的数值格式做一些修正。将质量比分方程(1)写成如下形式∂Y/∂t+∂f/∂x-Y∂g/∂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

N	Lerr1	order		Lerr1	order		Lerr1	order		Lerr1	order
	k=2			k=3			k=4			k=4
10	4.80×10-3	-		3.91×10-4	-		7.22×10-6	-		4.69×10-7	-
20	1.17×10-4	2.04		6.77×10-5	2.53		4.34×10-7	4.06		1.93×10-8	4.60
40	2.93×10-5	2.00		9.64×10-6	2.81		2.47×10-8	4.14		6.79×10-10	4.83
80	7.34×10-5	2.00		1.26×10-6	2.94		1.56×10-9	3.98		2.21×10-11	4.94
160	1.84×10-5	2.00		1.59×10-7	2.99		9.77×10-11	4.00		6.85×10-13	5.01

下载: 导出CSV 

| 显示表格

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

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

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

图  1  守恒谱体积格式的一维多介质运动界面问题

Figure  1.  One-dimensional moving interface problem for conservative spectral volume scheme

下载: 全尺寸图片 幻灯片

图  2  拟守恒谱体积格式的一维多介质运动界面问题

Figure  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  拟守恒谱体积格式的高压力比气液激波管问题

Figure  3.  High pressure ratio gas-liquid shock tube problem for quasi-conservative spectral volume scheme

下载: 全尺寸图片 幻灯片

图  4  拟守恒谱体积格式的高压力比气液激波管问题

Figure  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  拟守恒谱体积格式的三点问题密度等值线图

Figure  5.  Density contours of triple point problem for quasi-conservative spectral volume scheme

下载: 全尺寸图片 幻灯片

6.   结论

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