ABSTRACT

Referring to the construction idea of the MUSCL scheme, we present a construction method for high-order numerical fluxes. The higher order numerical flux is applied to the finite volume (difference) ENO, WENO and DG schemes to obtain single step higher order Semi-Lagrangian scheme. For the one-dimensional Euler equations, this paper presents a new feature line processing scheme in the feature space, which solves the difficulty that the Semi-Lagrangian method is difficult to generalize to multidimensional. Numerical experiments show that the new format is smaller than the original format, the efficiency is higher, and the simulation effect on the shock wave is also greatly improved.

Keywords:MUSCL, Semi-Lagrangian, High-Order Numerical Flux, One-Step High Order Semi-Lagrangian Scheme, Euler Equations, Sod Shock Tube Problem

一种高阶数值通量的探讨

温燕静¹，杨玉月^1*，魏雁霞²

¹湘潭大学数学系，湖南 湘潭

²中科软科技股份有限公司，北京

收稿日期：2019年5月9日；录用日期：2019年5月24日；发布日期：2019年5月31日

摘 要

参考MUSCL格式的构造思想，我们给出了一种高阶数值通量的构造方法。将高阶数值通量应用于有限体积(差分)ENO，WENO和DG等格式得到单步高阶Semi-Lagrangian格式。针对一维Euler方程组，本文在特征空间给出了一种新的特征线的处理方案，解决了Semi-Lagrangian方法难以推广到多维的难点。数值实验表明，新格式比原格式误差更小，效率更高，对激波的模拟效果也有较大提升。

关键词 :MUSCL，Semi-Lagrangian，高阶数值通量，单步高阶Semi-Lagrangian格式，Euler方程组， Sod激波管问题

Copyright © 2019 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

1959年，Godunov首次提出以Riemann问题解为基础来构造网格物理量均值的Godunov格式，它有效地解决了所有内部状态变量的局部黎曼问题 [1] [2]。随后，Van Leer首次讨论了二阶Godunov格式的构造得到了MUSCL (Monotone Upwind Scheme for Conservation Laws)格式，这个格式包括了保单调性的插值和利用Riemann解得到的迎风通量，提高了数值精度和分辨率 [2]。MUSCL格式已经被许多人研究 [3] [4] [5]，并被广泛用于科学和工程问题的模拟 [5]。

1959年，Wiin-Nielsen提出了一种修正轨迹方法 [6]，随后Sawyer将其发展为Semi-Lagrangian方法 [7]。该方法具有固定的(欧拉)数值网格，根据Riemann不变量沿特征线不变的性质，从而推进每个时间步，使得时间步长不受CFL条件的限制，减少计算工作量和数值扩散 [8]。Semi-Lagrangian方法成功地推广到了有限差分，有限元和谱方法，并广泛应用于天气预报，海洋学模型的数值模拟和线性对流的建模 [9] [10]。1982年，Robert将半隐式(SI)和半拉格朗日(SL)这两种方法相结合，在NWP模型中实现了显著的效率增益，为许多环境应用产生了许多有效的SISL模型 [11]。1999年，Strain将水平集方法和半拉格朗日方法相结合，能有效、自适应和模块化求解几何运动界面问题 [10] [12]。2006年，Bonaventura和Restelli等人将DG方法与半拉格朗日方法相结合为SLDG，能有效地融合了两种方法优点，又避免了缺点 [13]，对于shallow water方程模拟，SLDG具有较高的精度和有效性，并可自动选择局部逼近程度 [11]。2010年，Qiu等人将Strang时间分裂与空间高阶WENO重构相结合，应用于Semi-Lagrangian方法来求解Vlasov方程，比使用显式时间步长的Vlasov-Poisson系统的欧拉公式计算更高效，更灵活 [14]。

然而，对于一维或高维Euler方程组，求解特征线方程仍然是Semi-Lagrangian方法的难点。2011年，Lentine等人基于文 [15] 提出的分裂方法，用保守Semi-Lagrangian对流代替Meno对流格式求解方程组的对流项，但局部产生虚假振荡 [16]。目前尚未在现有文献中发现更好的处理方案。2015年，吴浪根据文 [17] 提出的方法，针对守恒律方程给出了几种高阶Semi-Lagrangian格式，通过数值模拟验证了该方法的高阶精度和无振荡性质 [18]。处理过程用到了一维Euler方程组的Riemann不变量，并反解出原物理量。这个过程增加了计算量，针对不同维数Euler方程组需要给出不同的Riemann不变量，不方便推广到多维问题，即便是通过交替方向法也不行。

本文给出了一种高阶数值通量的构造方法，将其应用于有限体积(差分) ENO，WENO和DG，即相应的Semi-Lagrangian方法。针对于一维Euler方程组，在特征空间上给出了一种新的处理特征线的方法。数值实验表明：高阶数值通量能应用于有限体积(差分) ENO，WENO和DG等格式，针对Euler方程组的特征线处理方法能保持空间离散格式本身的本质无震荡性，能获得更高的分辨率且具备更高的计算效率。

2. 二阶数值通量的推导

为了推导双曲型方程数值计算格式中的二阶数值通量，我们考察如下方程

$u_t+f{\left(u\right)}_x=0,\left(x,t\right)\in \left[a,b\right]\times \left[0,T\right],$ (1)

$u\left(x,0\right)=u_0\left(x\right).$ (2)

已知 $u\left(x,t_n\right)$ 在第j个单元( $I_j:=\left[x_{j-1/2},x_{j+1/2}\right]$ )端点的单侧极限分别为 $u_{j+1/2}^{-}\left(t_n\right)$ 和 $u_{j-1/2}^{+}\left(t_n\right)$ ，容易给出它在 $I_j$ 内的插值多项式

$p_j\left(x\right)=u\left(x,t_n\right)=u_{j+1/2}^{-}\left(t_n\right)+\frac{u_{j+1/2}^{-}\left(t_n\right)-u_{j-1/2}^{+}\left(t_n\right)}{h_j}\left(x-x_{j+1/2}\right),$ (3)

其中 $h_j=x_{j+1/2}-x_{j-1/2}$ 。

过点 $\left(x_{j+1/2},t_{n+1}\right)$ 的特征线与直线 $t=t_n$ 相交于一点X

$\frac{\text{d}}{\text{d}s}X\left(x,t;t+s\right)=\frac{\partial }{\partial u}f\left(u\left(X\left(x,t;t+s\right),t+s\right)\right)=h\left(X\left(x,t;t+s\right),t+s\right),$ (4)

其中 $X\left(x,t;t\right)=x$ 。利用一阶Runge-Kutta方法，可得

$X=x_{j+1/2}-\Delta sh\left(X\left(x_{j+1/2},t_n;t_n+s\right),t_n+s\right),$ (5)

将X代入 $u\left(x,t_n\right)$ ，得

$u\left(X,t_n\right)=u_{j+1/2}^{-}\left(t_n\right)+\frac{u_{j+1/2}^{-}\left(t_n\right)-u_{j-1/2}^{+}\left(t_n\right)}{h_j}\Delta sh\left(X\left(x_{j+1/2},t_n;t_n+s\right),t_n+s\right).$ (6)

其中 $u\left(X,t_n\right)$ 是近似代替 $u\left(x_{j+1/2},t_{n+1}\right)$ 。

最后，利用梯形公式可计算出 $\left[t_n,t_{n+1}\right]$ 内 $x_{j+1/2}$ 处的二阶数值通量

${\hat{f}}_{j+1/2}=\frac{f\left(u\left(X,t_n\right)\right)+f\left(u\left(x_{j+1/2},t_n\right)\right)}{2}.$ (7)

3. 二阶数值通量在DG格式的应用

为了将二阶数值通量应用于DG格式，我们将方程(1)两边同时乘 $\phi \left(x\right)$ ，并在时空单元 ${\Omega }_j$ 上求积分，

${\displaystyle {\iint }_{{\Omega }_j}\left({\partial }_t{u}_h+{\partial }_{x}f\left(u_h\right)\right)\phi \left(x\right)\text{d}x\text{d}t}=0,$ (8)

其中 $\phi \left(x\right)$ 为区间 $I_j$ 上的光滑函数。由格林公式，可得

$\begin{array}{l}{\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}\left(u_h\left(x,t_{n+1}\right)-u_h\left(x,t_n\right)\right)\phi \left(x\right)\text{d}x}+{\displaystyle {\int }_{t_n}^{t_{n+1}}f\left(u_h\left(x_{j+1/2},t\right)\right)\phi \left(x_{j+1/2}\right)\text{d}t}\\ -{\displaystyle {\int }_{t_n}^{t_{n+1}}f\left(u_h\left(x_{j-1/2},t\right)\right)\phi \left(x_{j-1/2}\right)\text{d}t}-{\displaystyle {\iint }_{{\Omega }_j}f(u_h){\partial }_x\phi \left(x\right)\text{d}x\text{d}t}=0,\end{array}$ (9)

我们再将方程(2)两边同乘 $\phi \left(x\right)$ 并在单元 $I_j$ 上积分得

${\displaystyle {\int }_{I_j}{u}_h}\left(x,0\right)\phi \left(x\right)\text{d}x={\displaystyle {\int }_{I_j}{u}_0}\left(x\right)\phi \left(x\right)\text{d}x.$ (10)

这样就得到了方程(1-2)的弱形式(9-10)。

参照经典二阶DG格式，取Legendre多项式 $P_l\left(x\right)$ 为局部基函数，并设 $u_h\left(x,t\right)={\displaystyle \underset{l=0}{\overset{k}{\sum }}{u}_j^l}\left(t\right){\phi }_l\left(x\right)$ ，其中 ${\phi }_l\left(x\right)=P_l\left(\xi \right)=P_l\left(\frac{x-x_j}{\Delta x/2}\right)$ 。将方程(9)中的 $\phi \left(x\right)$ 替换成 ${\phi }_l\left(x\right)$ ，则有

当l = 0时,(9)式简写为

$u_j^0\left(t_{n+1}\right)=u_j^0\left(t_n\right)-\frac{\Delta t}{\Delta x}\left({\hat{f}}_{j+1/2}-{\hat{f}}_{j-1/2}\right),$ (11)

当l = 1时，(9)式简写为

$u_j^1\left(t_{n+1}\right)=u_j^1\left(t_n\right)-\frac{3\Delta t}{\Delta x}\left({\hat{f}}_{j+1/2}+{\hat{f}}_{j-1/2}\right)+\frac{3}{\Delta x}{\displaystyle {\iint }_{{\Omega }_j}f\left(u_h\right){\partial }_x{\phi }_1}\left(x\right)\text{d}x\text{d}t.$ (12)

将(12)式中的二重积分项简记为K，利用如下梯形公式计算

$K=\frac{\Delta t}{\Delta x}\left({\displaystyle {\int }_{I_j}f\left(u_h\left(x,t_{n+1}\right)\right)\text{d}x}+{\displaystyle {\int }_{I_j}f\left(u_h\left(x,t_n\right)\right)\text{d}x}\right),$ (13)

其中 ${\displaystyle {\int }_{I_j}f\left(u_h\left(x,t_n\right)\right)\text{d}x}$ 继续利用梯形公式(或其他积分公式)， $u_h\left(x,t_{n+1}\right)$ 由方程(7,11)得到的单元均值 $u_j^0\left(t_{n+1}\right)$ 来近似代替，得

$K=\Delta t\left[f\left(u_j^0\left(t_{n+1}\right)\right)+\frac{1}{2}\left(f\left(u_h\left(x_{j+1/2},t_n\right)\right)+f\left(u_h\left(x_{j-1/2},t_n\right)\right)\right)\right].$ (14)

再利用上一节所介绍的二阶数值通量的推导方法，可求解(12)式中的 ${\hat{f}}_{j+1/2}$ 这样，我们得到相应的单步二阶DG格式(10)~(12)。

4. 一维Euler方程组上的二阶数值通量的推导

一维Euler方程组可以表示为

$u_t+{\left[F\left(u\right)\right]}_x=0,$ (15)

$u={\left(\rho ,\rho u,\rho E\right)}^{\text{T}},F\left(u\right)={\left(\rho u,\rho u^2+p,\left(\rho E+p\right)u\right)}^{\text{T}}$

其中压力 $p=\left(\gamma -1\right)\left(\rho E+\frac{1}{2}\rho u^2\right)$ ，比热比 $\gamma =1.4$ 。Jacobian矩阵 $A=\frac{\partial F\left(u\right)}{\partial \stackrel{\to }{u}}$ ，A的特征值为

${\lambda }_1=u-c,{\lambda }_2=u,{\lambda }_3=u+c,$

方程组(15)化成特征空间形式

$u_t+{\left(R\Lambda R^{-1}u\right)}_x=0,$

其中R为右特征矩阵， $R^{-1}$ 为左特征矩阵， $\Lambda =diag\left(u-c,u,u+c\right)$ 。

已知 $u\left(x,t_n\right)$ 在第j个单元( $I_j=\left[x_{j-1/2},x_{j+1/2}\right]$ )端点的单侧极限分别为 $u_{j+1/2}^{-}\left(t_n\right)$ 和 $u_{j-1/2}^{+}\left(t_n\right)$ 可以给出它在 $I_j$ 内的插值多项式

$p_j\left(x,t_n\right)=u\left(x,t_n\right)=u_{j+1/2}^{-}\left(t_n\right)+\frac{u_{j+1/2}^{-}\left(t_n\right)-u_{j-1/2}^{+}\left(t_n\right)}{h_j}\left(x-x_{j+1/2}\right),$

特征线方程 $\frac{\text{d}X}{\text{d}t}=-\lambda$ 与直线 $t_n$ 相交，有

$\frac{\text{d}}{\text{d}s}{X}_k\left(x,t;t+s\right)=-{\lambda }_k\left(X_k\left(x,t;t+s\right),t+s\right),k=1,2,3$ (16)

利用二阶Runge-Kutta法求解常微分方程(16)，再将 $X_k$ 代入 $u\left(x,t_n\right)$ 得到 $u_k\left(X_k,t_n\right)$ 。令

$W_k\left(X_k,t_n\right)=R^{-1}{u}_k\left(X_k,t_n\right)={\left(W_{k1},W_{k2},W_{k3}\right)}^{\text{T}}$ ，则 $W\left(x_{j+1/2},t_{n+1}\right)={\left(W_{11},W_{22},W_{33}\right)}^{\text{T}}$ 。

最后，利用梯形公式可计算出 $\left[t_n,t_{n+1}\right]$ 内 $x_{j+1/2}$ 处的二阶数值通量

${\hat{F}}_{j+1/2}=\frac{F\left(u\left(x_{j+1/2},t_n\right)\right)+F\left(u\left(x_{j+1/2},t_{n+1}\right)\right)}{2}=\frac{R\Lambda \left(W\left(x_{j+1/2},t_n\right)+W\left(x_{j+1/2},t_{n+1}\right)\right)}{2}.$ (17)

5. 数值算例

求解方程(1)的守恒型有限体积(差分)格式通常可以写成

$u_j^{n+1}=u_j^n-\frac{\Delta t}{\Delta x}\left({\hat{f}}_{j+1/2}-{\hat{f}}_{j-1/2}\right),$ (18)

我们将由公式(7)或(17))计算所得的通量替代(18)式中的 $\hat{f}$ ，这样就得到单步二阶格式.

算例1：一维Burger方程

考虑方程(1-2)，其中 $f\left(u\right)=u^2/2,\left(x,t\right)\in \left[0,2\text{π}\right]\times \left[0,2\right],u_0\left(x\right)=\sin x+0.5$ ，采用周期边界条件。为了考

表1. 单步二阶ENO，DG格式与原格式在t = 0.5时刻的数值误差表

表2. t = 0.5时三阶WENO格式的误差表

察格式的截断误差阶，我们给出 $t=0.5$ (尚未产生间断)的误差表，见表1~表2。其中后缀“HNF”表示使用公式(7)的单步二阶格式，另一个表示使用相应的二阶SSP Runge-Kutta时间离散格式。

由表1~表2可以看出：1) 本文给出的二阶数值通量能应用于ENO，WENO和DG等格式，得到相应的单步二阶格式。2) 单步二阶格式的计算误差比相应原格式更小(三阶WENO格式当网格充分细密时出现超收敛,误差更小)。3) 单步二阶格式所耗CPU时间更少，当N = 6400时，FV-ENO2-HNF格式比原格式节约了48%，FD-ENO2-HNF格式比原格式节约了56%，DG2-HNF格式比原格式节省了88%。

图1. t = 1.5真解与数值解的局部对比

当t = 1.5时，该问题已产生激波，见图1。由图可知：单步二阶格式在处理激波间断时具有更高的分辨率，模拟效果更好。

算例2：Sod激波管问题

考虑方程(15)，其中 $x\in \left[0,1\right]$ 。采用紧支边界，初始条件：t=0时刻，

$\left(\rho ,u,p\right)=\{\begin{array}{l}\left(1,0,1\right),x\le 0.5\\ \left(0.125,0,0.1\right),x>0.5\end{array}$

真解N=1600，数值解N=100时，使用二阶数值通量单步ENO格式与标准二阶FD-ENO2格式的对比(见图2~图3)。

图2. t = 0.15时真解与数值解的整体对比

由图2~图3可以看出：使用二阶数值通量的单步格式比原格式具有更高的分辨率。当N = 1000时，

图3. t = 0.15时真解与数值解的局部对比

FD-ENO2-HNF格式所消耗的CPU时间为282毫秒(ms)，标准FD-ENO2格式所消耗的CPU时间为400毫秒(ms)，FD-ENO-HNF格式比标准FD-ENO2格式节约了29.5%左右。

6. 总结

本文以二阶数值通量为例，给出了一种高阶数值通量的构造方法，能应用于有限体积(差分) ENO，WENO和DG格式得到单步高阶Semi-Lagrangian格式。格式天然具备守恒性、紧致性、并行性，且并行效率比时间离散应用SSP Runge-Kutta方法的相应格式更高。本文还给出了一种在特征空间处理特征线的方法，能应用于Euler方程组，比原格式具备更好的模拟效果和更高的效率。同时，这种方法能利用交替方向法推广到多维问题，并无需推导其Riemann不变量，也为Semi-Lagrangian方法在多维Euler方程组中的研究与应用提供了一种有效的途径。

文章引用

温燕静,杨玉月,魏雁霞. 一种高阶数值通量的探讨
An Approach of High Order Numerical Flux[J]. 应用数学进展, 2019, 08(05): 990-997. https://doi.org/10.12677/AAM.2019.85113

参考文献