1. 引言

一维双曲守恒律方程组的一般形式为

${\partial }_{t}u+{\partial }_{x}f\left(u\right)=0$ ， (1)

其中 $u=u\left(x,t\right)$ 是守恒型变量， $\left(x,t\right)\in R\times R^{+}$ ， $R\times R^{+}\to R^N\left(N\ge 1\right)$ ， $f\left(u\right)$ 是光滑的通量函数。双曲守恒律方程是计算流体力学中一类反映物理现象及其规律的重要方程，如气体动力学中的Euler方程、浅水动力学中的浅水波方程等，这些方程所涉及的物理模型在航天航空、气象、海洋等领域应用广泛，因此对该类方程的求解具有重要意义。研究非线性双曲守恒律的主要困难是：即使初值函数光滑，其解在时间推进的某个时刻也可能会出现间断，产生激波、接触间断等。间断解的出现使双曲守恒律方程的解违背了古典解理论，于是Lax在1954年提出了弱解 [1] 的概念，允许间断解存在，这为数值方法更好的发展打下了良好的基础。但是由于弱解不惟一，很多数值方法无法保证这种弱解具有物理意义。我们知道，任何物理的问题都存在粘性机制，而在数学上求解时可以认为粘性是接近于零的。当“粘性消失”时得到的Cauchy问题的解在任何时刻惟一且有物理意义，这种解称为熵解。1973年Lax从物理学中的热力学第二定律出发，在文献 [2] 中证明了熵稳定条件等价于这种“粘性消失”的机制，并由一个单元熵不等式来表示：

$U{\left(u\right)}_t+F{\left(u\right)}_x\le 0$ ， (2)

其中 $u$ 是惟一的具有物理意义的解， $U\left(u\right)$ 是 $u$ 的凸函数，称为熵函数，满足 ${U^{\prime }}^{\prime }\left(u\right)>0$ ， $F\left(u\right)$ 是熵通量函数。满足式(2)的 $U\left(u\right)$ ， $F\left(u\right)$ 被称为熵对。为了满足熵稳定条件，同时避免非物理解的产生，长期以来人们的做法是引入额外的耗散来控制熵的变化。

为合理地研究熵稳定与数值粘性之间的关系，1987年Tadmor在文献 [3] 里引入熵变量和熵势的概念，构造了一类二阶的熵守恒格式，其数值通量保持总熵不变。他还同时提出一种关于数值粘性更精确量化的比较原则，即：如果一个三点守恒型格式含有比熵守恒格式更多的数值粘性，则该三点守恒型格式是熵稳定的。但是具体应该添加多少数值粘性，一直是研究人员探索的问题。2004年，Tadmor提出一种逐段分解的显式构造方法来得到守恒型方程组的熵守恒数值通量，该方法可用于各种守恒型系统 [4] 。结合熵稳定格式的定义，Tadmor和Zhong把该方法应用于Navier-Stokes方程组 [5] 和浅水波方程 [6] ；数值实验表明这种方法只有在网格划分非常细密时才有效。2006年，Roe提出在熵守恒格式的基础上加上他本人提出的Roe格式的数值粘性的方法，得到了一类熵稳定格式，称为ERoe格式 [7] 。随后Roe和Ismail通过进一步分析得到：对于Euler方程，解在跨过激波时产生了激波强度立方量级的熵增，并在此基础上发展了对数值粘性项更精确量化的熵相容格式 [8] 。熵相容格式是目前对熵的变化估计得最精确的一种熵稳定格式，但Roe和Ismail的熵相容格式只有一阶精度，对间断的分辨率较低。

Van Leer提出的MUSCL数据重构方法 [9] 是求解双曲守恒律方程的一种高阶方法，通过MUSCL数据重构得到了许多高阶格式，如MUSCL-Hancock方法，广义黎曼问题法和斜率限制中心格式等 [9] ；2021年，任璇、封建湖提出了一种以MUSCL数据重构为基础的斜率限制器方法 [10] ，有效地提高了熵相容格式的精度和分辨率；随后沈亚玲等 [11] 将该斜率限制器方法应用到理想磁流体动力学(MHD)方程中，并通过数值模拟验证了该格式是求解MHD方程的一种较为理想的方法。本文进一步研究熵相容格式，利用MUSCL-Hancock方法构造一种在时间和空间上全离散的二阶高分辨率熵相容格式，数值结果表明新格式具有稳健性、高精度性和无振荡性等良好特征并且能够提高计算效率。

2. 熵相容格式

2.1. 一维双曲守恒律方程

2.1.1. 熵守恒格式

一维双曲守恒律方程(1)的半离散格式为

$\frac{\text{d}}{\text{d}t}{u}_i\left(t\right)=-\frac{1}{\Delta x}\left(f_{i+1/2}-f_{i-1/2}\right)$ ， (3)

熵守恒格式能够保持总熵不变，在离散条件下满足单元熵等式

$\frac{\text{d}}{\text{d}t}U\left(u_i\left(t\right)\right)+\frac{1}{\Delta x}\left(F_{i+1/2}-F_{i-1/2}\right)=0$ ， (4)

其中 $U\left(u_i\left(t\right)\right)$ 为凸的熵函数， $F_{i+1/2}$ 为数值熵通量函数。Tadmor在文献 [3] 中给出满足上述熵守恒条件的数值通量 $f_{i+1/2}^C$ ，该数值通量与熵变量和熵势的关系为：

${\left[v\right]}_{i+1/2}^T{f}_{i+1/2}^C={\left[\psi \right]}_{i+1/2}$ ， (5)

式(4)中的数值熵通量为 $F_{i+1/2}=v_{i+1/2}^T{f}_{i+1/2}^C-{\bar{\psi }}_{i+1/2}$ ， ${[\cdot ]}_{i+1/2}={(\cdot )}_{i+1}-{(\cdot )}_i$ ， ${\overline{[\cdot ]}}_{i+1/2}=\frac{{(\cdot )}_{i+1}-{(\cdot )}_i}{2}$ ，熵势 $\psi =v^{T}f-F$ ，则称 $f_{i+1/2}^C$ 是熵守恒格式的数值通量。

当计算标量情况下的双曲守恒律方程时，用公式(5)进行计算，可得到熵守恒通量 $f_{i+1/2}^C$ 的惟一表达形式

$f_{i+1/2}^C=\frac{{\left[\psi \right]}_{i+1/2}}{{\left[v\right]}_{i+1/2}}=\frac{{\psi }_{i+1}-{\psi }_i}{v_{i+1}-v_i}$ ， (6)

2.1.2. 熵相容格式

熵守恒格式能够保持总熵不变，在光滑区域具有良好的表现且具有二阶精度，但是当间断解出现时熵会增加，使得熵守恒格式的数值解出现伪振荡。2006年Roe提出在熵守恒的基础上添加一个Roe-耗散项得到熵稳定格式，1999年Barth证明了 $u_v\approx R{R}^T$ [12] ，将Roe-耗散项进一步调整为如下形式

$D_{i+1/2}^{ES}=\frac{1}{2}{R}_{i+1/2}\left|{\Lambda }_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ ， (7)

其中 $\left|{\Lambda }_{i+1/2}\right|=diag\left(\left|{\lambda }_{i+1/2}^1\right|,\cdots ,\left|{\lambda }_{i+1/2}^N\right|\right)$ 是 $A\left(u\right)=\frac{\text{d}f}{\text{d}u}$ 的特征值矩阵， $R$ 是 $A\left(u\right)=\frac{\text{d}f}{\text{d}u}$ 的右特征向量构成的矩阵，得到熵稳定格式的数值通量 $f_{i+1/2}^{ES}$ ：

$f_{i+1/2}^{ES}=f_{i+1/2}^C-D_{i+1/2}^{ES}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left|{\Lambda }_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ 。 (8)

与熵守恒格式相比，熵稳定格式只多了耗散项 $D_{i+1/2}^{ES}$ ，但是这个耗散项有时不足以抵消跨越间断所产生的熵增，这使得在数值实验中，该耗散项在激波附近仍然可能出现伪振荡。对于Euler方程，Ismail和Roe曾经证明了在有激波时系统的熵增约为密度跳跃的立方 [8] ，因此需要在上述的熵稳定通量中再添加一项使之产生 $O\left({\left[u\right]}^3\right)$ 的熵增，本文将采用文献 [10] 中的熵相容通量

$\begin{array}{c}{f}_{i+1/2}^{EC}=f_{i+1/2}^C-D_{i+1/2}^{ES}-D_{i+1/2}^{EP}\\ \approx f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle \underset{k=1}{\overset{q}{\sum }}{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}\end{array}$ ， (9)

其中， $s_k$ 是高斯点， $w_k$ 是与高斯点对应的系数，积分区间为[0, 1]，在这里采用Gauss-Legendre积分公式，此时 $s_{1,3}=\frac{1}{2}\pm \frac{\sqrt{15}}{10}$ ， $s_2=\frac{1}{2}$ ， $w_{1,3}=\frac{5}{9}$ ， $w_2=\frac{8}{9}$ 。

记 $D_{i+1/2}=\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle \underset{k=1}{\overset{q}{\sum }}{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|$ ，则上式可简记为

$f_{i+1/2}^{EC}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left|D_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ 。 (10)

在标量情况下，熵相容格式的数值通量为

$\begin{array}{c}{f}_{i+1/2}^{EC}=f_{i+1/2}^C-D_{i+1/2}^{ES}-D_{i+1/2}^{EP}\\ =f_{i+1/2}^C-\frac{1}{2}\left|a\left(u_{i+1/2}\right)\right|{\left[v\right]}_{i+1/2}-\frac{1}{2}\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)a\left(u\left(\xi \right)\right)U_{uu}\left(\xi \right)\text{d}\xi }\right|{\left[v\right]}_{i+1/2}\end{array}$ 。 (11)

其中 $a\left(u\right)$ 是特征速度，记 $D_{i+1/2}=\left|a\left(u_{i+1/2}\right)\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)a\left(u\left(\xi \right)\right)U_{uu}\left(\xi \right)\text{d}\xi }\right|$ ，则上式可简

化为

$f_{i+1/2}^{EC}=f_{i+1/2}^C-\frac{1}{2}{D}_{i+1/2}{\left[v\right]}_{i+1/2}$ 。 (12)

在文献 [10] 中已证明该熵相容格式是收敛的。

对于一维无粘Burgers方程

$u_t+f{\left(u\right)}_x=0$ ， (13)

其中， $f\left(u\right)=u^2/2$ ，熵函数为 $U\left(u\right)=u^2/2$ ，熵通量函数为 $F\left(u\right)=2u^2/3$ ，特征值为 $a\left(u\right)=u$ ，进而得到熵变量 $v\left(u\right)=U^{\prime }\left(u\right)=u$ ，熵势 $\psi \left(u\right)=v\left(u\right)f\left(u\right)-F\left(u\right)=u^3/6$ ，由式(6)可得出Burgers方程的熵守恒数值通量为 $f_{j+1/2}^C=\left(u_j^2+u_j{u}_{j+1}+u_{j+1}^2\right)/6$ ，最终得到Burgers方程的熵稳定数值通量和熵相容数值通量分别为

$f_{i+1/2}^{ES}=f_{i+1/2}^C-\frac{1}{4}\left|u_i+u_{i+1}\right|\left(u_{i+1}-u_i\right)$ ， (14)

$f_{i+1/2}^{EC}=f_{i+1/2}^C-\frac{1}{2}\left(\frac{1}{2}\left|u_i+u_{i+1}\right|+\frac{1}{6}\left|u_{i+1}-u_i\right|\right)\left(u_{i+1}-u_i\right)$ 。 (15)

2.2. 一维浅水波方程

考虑描述浅水情况下流体运动规律的浅水波方程，具有以下守恒格式：

$\frac{\partial }{\partial t}\left[\begin{array}{c}h\\ hu\end{array}\right]+\frac{\partial }{\partial x}\left[\begin{array}{c}hu\\ g{h}^2/2+h{u}^2\end{array}\right]=0$ ， (16)

其中 $u={\left(h,hu\right)}^{\text{T}}$ ，h是水深，u是速度，g是重力加速度，取熵函数 $U\left(u\right)=\left(g{h}^2+u^{2}h\right)/2$ ，熵通量为 $F=h{u}^3/2+gu{h}^2$ ，熵变量为 $v=U^{\prime }\left(u\right)={\left[gh-u^2/2,u\right]}^{\text{T}}$ 。熵守恒的数值通量采用Fjordholm在文献 [13] 所提出的浅水波方程的显式熵守恒格式的数值通量：

$F_{i+\frac{1}{2}}^C=\left(\begin{array}{c}{\bar{h}}_{i+\frac{1}{2}}{\bar{u}}_{i+\frac{1}{2}}\\ \frac{g}{2}{\left(\overline{h^2}\right)}_{i+\frac{1}{2}}+{\bar{h}}_{i+\frac{1}{2}}{\left({\bar{u}}_{i+\frac{1}{2}}\right)}^2\end{array}\right)$ ， (17)

其中 $\overline{(\cdot )}=\left({(\cdot )}_i+{(\cdot )}_{i+1}\right)/2$ ， $\overline{{(\cdot )}^2}=\left({(\cdot )}_i^2+{(\cdot )}_{i+1}^2\right)/2$ 。

最终得到浅水波方程的熵稳定格式和熵相容格式的数值通量分别为

$f_{i+1/2}^{ES}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left|{\Lambda }_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ ， (18)

$f_{i+1/2}^{EC}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle \underset{k=1}{\overset{q}{\sum }}{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ ， (19)

其中 $s_{1,3}=\frac{1}{2}\pm \frac{\sqrt{15}}{10}$ ， $s_2=\frac{1}{2}$ ， $w_{1,3}=\frac{5}{9}$ ， $w_2=\frac{8}{9}$ ，对应雅可比矩阵的特征值矩阵为

$\Lambda =\left[\begin{array}{cc}u-\sqrt{gh}& \\ & u+\sqrt{gh}\end{array}\right]$ ， (20)

使用Roe平均值计算的平均特征值和特征向量矩阵为

${\Lambda }_{i+1/2}=\left[\begin{array}{cc}\bar{u}-\sqrt{g\bar{h}}& \\ & \bar{u}+\sqrt{g\bar{h}}\end{array}\right]$ ， $R_{i+1/2}=\frac{1}{\sqrt{2g}}\left[\begin{array}{cc}1& 1\\ \bar{u}-\sqrt{g\bar{h}}& \bar{u}+\sqrt{g\bar{h}}\end{array}\right]$ 。 (21)

2.3. 一维欧拉方程

考虑无粘可压缩气体动力学中的Euler方程

$\frac{\partial }{\partial t}\left[\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right]+\frac{\partial }{\partial x}\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \left(E+p\right)u\end{array}\right]=0$ ， (22)

其中 $\rho ,u,p$ 分别代表气体的密度，速度，压力。总熵定义为 $E=p/\left(\gamma -1\right)+\rho u^2/2$ ，总焓为 $H=\left(E+p\right)/\rho$ ， $\gamma$ 为理想气体数字，取 $\gamma =1.4$ 。Euler方程的熵对为

$U=\frac{-\rho S}{\gamma -1}$ ， $F=\frac{-\rho uS}{\gamma -1}$ ， (23)

其中S为物理熵。由熵变量 $v=U_u$ 可知

$v={\left[\frac{\gamma -S}{\gamma -1}-\frac{\rho u^2}{2p},\frac{\rho u}{p},-\frac{\rho }{p}\right]}^{\text{T}}$ 。 (24)

本文采用Ismail和Roe所提出的Euler方程的熵守恒格式 [8] ，具体的熵守恒格式的数值通量为

$f_{i+1/2}^C=\left[\begin{array}{c}\hat{\rho }\hat{u}\\ \hat{\rho }{\hat{u}}^2+{\hat{p}}_1\\ \hat{\rho }\hat{u}\hat{H}\end{array}\right]$ ， (25)

其中参数变量Z为

$Z={\left[z_1,z_2,z_3\right]}^T=\sqrt{\rho /p}{\left[1,u,p\right]}^T$ ， (26)

$\widehat{(\cdot )}$ 表示在单元交界面处的平均值，满足：

$\hat{\rho }=\overline{z_1}{z}_3^{\ln },\hat{u}=\overline{z_2}/\overline{z_1},{\hat{p}}_1=\overline{z_3}/\overline{z_1}\hat{a}=\sqrt[2]{\gamma {\hat{p}}_2/\hat{\rho }}$ ， (27)

${\hat{p}}_2=\frac{\gamma +1}{\gamma }\frac{z_3^{\ln }}{z_1^{\ln }}+\frac{\gamma -1}{2\gamma }\frac{{\bar{z}}_3}{{\bar{z}}_1},\hat{H}=\frac{{\hat{a}}^2}{\gamma -1}+\frac{{\hat{u}}^2}{2}$ ， (28)

在上式中， $\bar{a}=\frac{a_l+a_r}{2}$ 表示算术平均， $a^{\ln }=\frac{a_l-a_r}{\ln a_i-\ln a_r}\left(a_l\ne a_r\right)$ ， $a^{\ln }=a_l,\left(a_l=a_r\right)$ ，表示对数平均。在这里选用文献 [10] 中的熵稳定通量和熵相容通量，分别为

$f_{i+1/2}^{ES}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}{\Lambda }_{i+1/2}{S}_{i+1/2}{R}_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ ， (29)

$f_{i+1/2}^{EC}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}{D}_{i+1/2}{S}_{i+1/2}{R}_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}$ ， (30)

其中 $D_{i+1/2}=\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\sum }_{k=1}^q{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)$ ， $S_{i+1/2}=diag\left(\frac{\hat{\rho }}{2\gamma },\frac{\left(\gamma -1\right)\hat{\rho }}{\gamma },\frac{\hat{\rho }}{2\gamma }\right)$ ， $s_{1，3}=\frac{1}{2}\pm \frac{\sqrt{15}}{10}$ ， $s_2=\frac{1}{2}$ ， $w_{1,3}=\frac{5}{9}$ ， $w_2=\frac{8}{9}$ ， $\Lambda$ 就是雅可比矩阵A的特征值矩阵，有

$\Lambda =diag\left(u-a,u,u+a\right)$ 。 (31)

3. 基于MUSCL方法的熵相容格式

本文所涉及的熵稳定格式、熵相容格式和基于MUSCL数据重构的熵相容格式，为了在时间方向达到二阶精度，对半离散格式均采用二阶精度的强稳定的两步Runge-Kutta方法 [14]

$\{\begin{array}{l}{u}^{*}=u^n+\Delta tL\left(u^n\right)\\ u^{**}=u^{*}+\Delta tL\left(u^{*}\right)\\ u^{n+1}=0.5\cdot \left(u^n+u^{**}\right)\end{array}$ ， (32)

在 $t_n$ 层的单元边界处进行重构， $u\left(x\right)$ 表示双曲守恒律方程的精确解，在每个控制单元 $I_i=\left[x_{i-1/2},x_{i+1/2}\right]$ 上可以表示为分段的精确解 $u_i\left(x\right)$ ，在这里对每个 $u_i\left(x\right)$ 进行线性重构，记重构后的函数为 $R_i\left(x\right)$ ，表示如下：

$R_i\left(x\right)=u_i^n+{\left(u_x\right)}_i\left(x-x_i\right),x_{j-1/2}<x<x_{j+1/2}$ ， (33)

$R_i\left(x\right)$ 在左右两个交界面处的极限值分别为

$R_i^L=R_i\left(x_{i-1/2}\right)=u_i^n-\frac{1}{2}\Delta x\cdot {\left(u_x\right)}_i$ ， $R_i^R=R_i\left(x_{i+1/2}\right)=u_i^n+\frac{1}{2}\Delta x\cdot {\left(u_x\right)}_i$ 。 (34)

$R_i^L$ 和 $R_i^R$ 被称为边界外推值。通过MUSCL数据重构，得到的基于MUSCL数据重构的熵相容格式在大梯度附近会产生杂散振荡，为此我们采用广义的minmod限制器，使得解在光滑区域选择高阶的数据重构，在间断区域选择原始数据，从而避免产生伪振荡。具体形式如下(在第5节的算例中 $\alpha$ 值取2)：

${\left(u_x\right)}_i=\min \mathrm{mod}\left(\alpha \frac{u_i-u_{i-1}}{\Delta x},\frac{u_{i+1}-u_{i-1}}{2\Delta x},\alpha \frac{u_{i+1}-u_i}{\Delta x}\right),1<\alpha <2$ ， (35)

其中，minmod函数定义为

$\min \mathrm{mod}\left(a_1,a_2,a_3,\cdots \right):=\{\begin{array}{l}\min \left\{a_i\right\},a_i>0,\forall i\\ \max \left\{a_i\right\},a_i<0,\forall i\\ 0,\text{otherwise}\end{array}$ ， (36)

将公式(34)分别代入Burgers方程、浅水波方程、Euler方程的熵相容通量，得到基于MUSCL数据重构的熵相容通量：

$f_{i+1/2}^{EC-MUSCL}=f_{i+1/2}^C-\frac{1}{2}\left(\frac{1}{2}\left|R_{i+1}^L+R_i^R\right|+\frac{1}{6}\left|R_{i+1}^L-R_i^R\right|\right)\left(R_{i+1}^L-R_i^R\right)$ ， (37)

$f_{i+1/2}^{EC-MUSCL}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle \underset{k=1}{\overset{q}{\sum }}{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)R_{i+1/2}^T{\left[v\right]}_{i+1/2}^{Limited1}$ ， (38)

$f_{i+1/2}^{EC-MUSCL}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\sum }_{k=1}^q{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)S_{i+1/2}{R}_{i+1/2}^T{\left[v\right]}_{i+1/2}^{Limited1}$ ， (39)

其中 ${\left[v\right]}_{i+1/2}^{Limited1}$ 是重构后的熵变量值。在双曲守恒律方程中，上述的构造方法是先将守恒变量 $u$ 转化为原始变量进行重构，再将重构后的值代回熵变量，最终得到基于MUSCL重构的高分辨率熵相容通量 $f_{i+1/2}^{EC-MUSCL}$ 。

4. 基于MUSCL-Hancock方法的熵相容格式

第3节提及的两步龙格库塔法需要在单元边界处需要求解两次数值通量，但采用MUSCL-Hancock方法来更新下一时间层的物理量，是一个全离散过程，在单元交界面处只需求解一次数值通量，从而可以有效地提高格式的计算效率。计算过程如下：

第一步：在 $t_n$ 层重构 $u_i\left(x\right)$ 在单元交界面处的左右极限值，与第3节的过程一致，结果为 $R_i^L$ ， $R_i^R$ ；

第二步：利用下式将 $t_n$ 层的边界外推值 $R_i^L$ ， $R_i^R$ 推进到 $t_{n+1/2}$ 层得到相应的边界外推值：

${\bar{R}}_i^L=R_i^L-\frac{1}{2}\frac{\Delta t}{\Delta x}\left[f\left(R_i^R\right)-f\left(R_i^L\right)\right]$ ， ${\bar{R}}_i^R=R_i^R-\frac{1}{2}\frac{\Delta t}{\Delta x}\left[f\left(R_i^R\right)-f\left(R_i^L\right)\right]$ ， (40)

第三步：利用4.2得到的边界外推值分别代入Burgers方程、浅水波方程、Euler方程的熵相容通量得到基于MUSCL-Hancock方法的熵相容通量：

$f_{i+1/2}^{EC-MHM}=f_{i+1/2}^C-\frac{1}{2}\left(\frac{1}{2}\left|{\bar{R}}_{i+1}^L+{\bar{R}}_i^R\right|+\frac{1}{6}\left|{\bar{R}}_{i+1}^L-{\bar{R}}_i^R\right|\right)\left({\bar{R}}_{i+1}^L-{\bar{R}}_i^R\right)$ ， (41)

$f_{i+1/2}^{EC-MHM}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle \underset{k=1}{\overset{q}{\sum }}{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)R_{i+1/2}^T{\left[v\right]}_{i+1/2}^{Limited2}$ ， (42)

$f_{i+1/2}^{EC-MHM}=f_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\sum }_{k=1}^q{w}_k\left(2s_k-1\right)\Lambda \left(u\left(s_k\right)\right)}\right|\right)S_{i+1/2}{R}_{i+1/2}^T{\left[v\right]}_{i+1/2}^{Limited2}$ ， (43)

在这里熵相容的通量的计算方法与第3节相同。

利用上述得到的基于MUSCL-Hancock方法的熵相容通量 $f_{i+1/2}^{EC-MHM}$ ，在这里我们可以更新 $t_{n+1}$ 时刻的单元平均值，公式如下：

$u_j^{n+1}=u_j^n-\frac{\Delta t}{\Delta x}\left(f_{i+1/2}^{EC-MHM}-f_{i-1/2}^{EC-MHM}\right)$ 。 (44)

5. 数值算例

Burgers方程、一维浅水波方程、一维欧拉方程的熵稳定通量分别为(14)，(18)，(29)，熵相容通量分别为(15)，(19)，(30)，基于MUSCL方法的熵相容通量分别为(37)，(38)，(39)，基于MUSCL-Hancock方法的熵相容通量分别为(41)，(42)，(43)，进行数值模拟，通过数值结果说明所构造的EC-MHM格式具有更好的性质。

符号约定：

ES：熵稳定格式；EC：熵相容格式；EC-MUSCL：基于MUSCL重构的熵相容格式；EC-MHM：基于MUSCL-Hancock方法的熵相容格式；Exact：精确解；Ref：参考解。

5.1. 一维Burgers方程

对于一维Burgers方程的不同初值问题，在下面的实验中均采取40网格进行计算。

5.1.1. 间断初值问题

在区间[−1, 1]上考虑初值问题

$u\left(x,0\right)=\{\begin{array}{l}-1,\text{if}1/3\le \left|x\right|\le 1\\ 1,\text{if}\left|x\right|<1/3\end{array}$ ，

该问题的解析解为

$u\left(x,t\right)=\{\begin{array}{ll}1,\hfill & \text{if}t-1/3<x<1/3\hfill \\ \left(x+\frac{1}{3}\right)/t,\hfill & \text{if}-t-1/3<x<t-1/3\hfill \\ -1,\hfill & \text{otherwise}\hfill \end{array}$ 。

使用Neumann边界条件，计算至t = 0.32，CFL = 0.3，计算结果如图1所示，ES格式在激波附近会产生振荡，相比之下，EC格式具有更高的分辨率，并且能够更准确的捕捉激波，但二者在稀疏波区域的音速点处仍会产生膨胀激波。而EC-MUSCL格式和EC-MHM格式则具有更好的数值结果，能够准确捕捉稀疏波和激波，并消除振荡。从熵耗散图可以看出EC-MUSCL格式和EC-MHM格式的熵耗散介于ES格式和EC格式之间，既提高了稳定性也提高了计算结果的准确性；二者相比，EC-MHM格式比EC-MUSCL格式耗散更少，同时也具有更高的分辨率。

5.1.2. 可压缩波问题

在区域[−1, 1]上考虑初值问题

$u\left(x,0\right)=\{\begin{array}{ll}{u}_0,\hfill & \text{if}\left|x\right|>1\hfill \\ u_0-u_1\ast \sin \left(\pi x\right),\hfill & \text{if}\left|x\right|\le 1\hfill \end{array}$ ，

其中振幅 $u_1=0.5$ ， $u_0=0$ 。它的解析解为

$u\left(x,t\right)=\frac{u_0}{u_1}-\sin \left(\pi x-\frac{u}{u_1}{u}_{1}t\right)$ 。

使用Neumann边界条件，计算至t₁ = 0.32，CFL = 0.3，计算结果如图2所示，在解的连续区域，ES格式和EC格式对解的捕捉效果差别甚微，而EC-MUSCL格式和EC-MHM格式的计算结果则更优；在表1中给出了EC-MUSCL格式和EC-MHM格式在t₁ = 0.32时的L₁误差和阶数，结果表明EC-MHM格式比EC-MUSCL格式具有更小的误差和更高的精度，因此EC-MHM格式的计算结果也是最佳的；在表2中给出了EC-MUSCL格式和EC-MHM格式在t₁ = 0.32时程序的计算时间及EC-MHM格式和EC-MUSCL格式计算时间的比值，结果表明EC-MHM格式具有更高的计算效率，这与前文的理论分析也是一致的。

5.2. 一维浅水波方程

考虑在不同初始条件的溃坝问题，都使用Neumann边界条件，这些问题的精确解都是由一个向左的稀疏波和一个向右的激波组成。在区域[−1, 1]上，设初始条件为

$h\left(x,0\right)=\{\begin{array}{ll}2\hfill & x<0\hfill \\ 1.5\hfill & x\ge 0\hfill \end{array}$ ， $u\left(x,0\right)=0$ ，

取100个网格点，CFL = 0.45，计算至t = 0.4。

该问题的参考解为1200网格点时Roe格式的计算结果。从图3中可以看出，ES格式和EC格式均可以抑制膨胀激波的产生，但是二者都有明显的抹平现象，而EC-MUSCL格式和EC-MHM格式具有更小的熵耗散，因此可以更好地捕捉稀疏波和激波；与EC-MUSCL格式相比，EC-MHM格式对稀疏波和激波的捕捉效果都要更好，具有更高的分辨率。

5.3. Euler方程

5.3.1. Lax激波管问题

考虑在区间[0, 1]上的初值问题

$\left(\rho ,u,p\right)=\{\begin{array}{ll}\left(0.445,0.698,3.528\right),\hfill & \text{if}x<0.5\hfill \\ \left(0.5,0,0.571\right),\hfill & \text{otherwise}\hfill \end{array}$ ，

使用Neumann边界条件，CFL = 0.15，取200个网格点，计算至t = 0.16。

该问题的精确解由一个左稀疏波，一个中间的接触间断和一个右激波。关于密度的计算结果和总熵随时间的变化图如图4所示，EC格式和ES格式在间断位置的抹平都比较严重，相比之下，EC-MUSCL格式和EC-MHM格式在激波和接触间断处则表现出更准确的捕捉效果；EC-MUSCL格式在接触间断结束的位置会出现轻微的振荡，而EC-MHM格式能够更好的捕捉接触间断而不会产生振荡。

5.3.2. Shu-Osher问题

该问题是关于激波和熵波相互作用的问题 [15] ，它的解既包含激波又包含复杂的连续解，考虑在区间[−5, 5]上的Shu-Osher问题，给定初始条件

$\left(\rho ,u,p\right)=\{\begin{array}{ll}\left(3.857143,2.629369,10.333333\right),\hfill & -5\le x<-4\hfill \\ \left(1+0.2\sin \left(5x\right),0,1\right),\hfill & -4\le x<5\hfill \end{array}$ ，

使用Neumann边界条件，CFL = 0.1分别取400个网格点和800个网格点，计算至t = 1.8。

以4000网格的高分辨率熵相容格式 [16] 的数值结果作为参考解。关于密度的计算结果如图5所示，ES格式和EC格式能够捕捉精确解的结构，但是抹平现象严重，而EC-MHM格式和EC-MUSCL格式则具有更高的分辨率，从密度图中可以明显看出，EC-MHM格式比EC-MUSCL格式能够更准确地捕捉精确解，这也说明EC-MHM具有一定的实用性。

5.3.3. 爆炸波问题

考虑在区域[−0.5, 0.5]上的爆炸波问题，给定初始条件

$\left(\rho ,u,p\right)=\{\begin{array}{l}\left(1.0,0,1000\right),-0.5<x<-0.4\\ \left(1.0,0,0.01\right),-0.4<x<0.4\\ \left(1.0,0,100\right),0.4<x<0.5\end{array}$ ，

使用反射边界条件，取CFL = 0.005，计算至0.038，空间网格数为N = 400和800。

爆炸波解的结构很复杂，常用来检验数值格式捕捉激波的能力。从图6中可以看出四种格式都可以用来计算这个问题，但是EC-MHM具有更好的分辨力。因此在解决复杂波问题时，EC-MHM格式既能够更好地模拟精确解的结构同时还可以提高计算效率，节省时间。

6. 结论

本文构造了一种基于MUSCL-Hancock方法的熵相容格式，相比其它的熵相容格式而言，新的EC-MHM格式具有二阶精度，同时它只需在每个时间层 $\left[t_n,t_{n+1}\right]$ 的中间时刻t_n+1计算一次数值通量，从而达到提高计算效率的目的；若干数值结果表明，新的EC-MHM格式具有高精度、高分辨率和鲁棒等特性。该格式是求解双曲守恒律方程的一种较好的方法，后续可以将这种方法进一步推广到二维甚至更高维。

NOTES

^*通讯作者。

参考文献