浅水磁流体方程的高分辨率熵稳定格式

High-Resolution Entropy Stable Scheme for Shallow Water Magnetohydrodynamics Equations

摘要: 本文构造了一种求解浅水磁流体方程(Shallow Water Magnetohydrodynamics Equations，简称为SWMHD方程)的高分辨率熵稳定格式(ES2)。主要构造过程是将一种基于MUSCL (Monotone Up-stream-centred Scheme for Conservation Laws)型重构方法的新型斜率限制器应用于SWMHD方程的熵稳定格式(ES)中，得到了具有二阶精度的高分辨率熵稳定格式，该格式在间断区域通过合理控制耗散的方式有效地抑制了非物理现象的产生。将新构造的格式(ES2)和熵稳定格式(ES)，熵守恒格式(C)应用于一维和二维SWMHD方程问题求解中，数值结果表明，该格式能够准确地捕捉解的结构，在光滑区域达到了高阶精度，且在间断区域没有非物理振荡。

Abstract: In this paper, a high-resolution entropy stable scheme (ES2) is constructed to solve the shallow wa-ter magnetohydrodynamics equations (SWMHD equation). The main construction process is to ap-ply a new slope limiter based on MUSCL (Monotone Upstream-centred Scheme for Conservation Laws) type reconstruction method to the entropy stable scheme (ES) of the SWMHD equation, and obtain a high-resolution entropy stable scheme with second-order accuracy, which effectively sup-presses the occurrence of non-physical phenomena by reasonably controlling dissipation in the discontinuous region. The newly constructed scheme (ES2), entropy stable scheme (ES), and entro-py conservation scheme (C) are applied to solve the one-dimensional and two-dimensional SWMHD equation problems, and the numerical results show that the scheme can accurately capture the structure of the solutions, achieve high-order accuracy in the smooth region, and have no non-physical oscillations in the discontinuous region.

1. 引言

浅水波方程广泛应用于大气流动、潮汐、风暴潮、河流和海岸流、湖泊流、海啸等。浅水磁流体方程(SWMHD)是在浅水波方程的基础上考虑了磁场的影响，最早由Gilman [1] 提出用于描述太阳胸廓线 [2] (太阳内部的辐射区和对流区被一层薄层分隔开，这层薄层被称为太阳胸廓线)的现象。同时SWMHD方程也是理想磁流体力学(MHD)方程在自由表面、浅层和导电流体的情况下的近似值，与MHD方程一样，SWMHD方程也需要对无散约束条件进行特殊处理。虽然SWMHD方程不是严格的双曲型方程，但其波的结构与流体力学方程相似，因而我们可以将求解一般双曲守恒律的方法推广到SWMHD方程上。

双曲守恒律方程的理论研究及相应求解的数值方法研究已有一百多年的历史。有关双曲守恒律方程数值解的研究发现：即使其初值函数光滑，其解在某个时刻也可能会出现间断，该问题的出现违反了古典解的有关理论。为解决该类问题，Lax在1954年提出了弱解 [3] 的概念，允许间断解存在；然而弱解并不唯一，由此产生了守恒型差分格式的解是否收敛于物理解的问题，这就需要额外的标准来寻找满足相关物理特性的唯一解。在数学上，满足“粘性消失”的Cauchy问题的解在任何时刻都是唯一且有物理意义，这种解称为熵解。1973年Lax在文 [4] 中证明了熵稳定条件等价于这种“粘性消失”的机制，并由一个单元熵不等式来表示。在实际物理问题中，跨过激波时熵是增加的，但反映在数学模型中，恰好差一个负号，表现为熵减，所以要求所采用的数值方法必须引起熵的耗散，才能达到熵稳定的要求。1987年Tadmor引入熵变量和熵势的概念，构造了一类二阶的熵守恒格式，其数值通量保持总熵不变 [5] 。2006年，Roe提出在熵守恒格式的基础上添加Roe格式的数值粘性项，得到了一类熵稳定格式 [6] ，称为ERoe格式。2009年Roe和Ismail通过进一步分析得到：解在跨过激波时产生了激波强度立方量级的熵增，在此基础上发展了对数值粘性项更精确量化的熵相容格式 [7] 。2015年，刘友琼、封建湖、任炯 [8] 提出了一种复合型通量限制器，得到了一类基于通量限制器的高分辨率熵相容格式。2021年，任璇、封建湖提出了一种以MUSCL方法为基础的新的斜率限制器 [9] ，有效地提高了熵相容格式的精度和分辨率。随后沈亚玲 [10] ，高凡琪 [11] 等将此类新型斜率限制器运用于磁流体方程和相对论流体力学方程的数值求解中，并通过大量算例证实该方法是一种效果颇佳的求解方法。

SWMHD方程的现有数值研究包括演化Galerkin格式 [12] 、时空守恒单元/求解单元(CE/SE)方法 [13] 、中心逆风格式 [14] 、Roe型格式 [15] 、高阶CE/SE格式 [16] 等，现有的数值方法奠定了数值求解SWMHD方程的基础，但是使用熵稳定格式求解SWMHD方程的方法很少。近来Winters和Gassner [17] 在2016年提出了SWMHD方程在平坦地形下的二阶熵稳定有限体积格式(满足半离散熵不等式)，但是其格式精度较低，在间断处出现了一定程度的抹平现象，随后在2021年段俊明和汤华中 [18] 等构造出了非平坦地形下SWMHD方程的高阶熵稳定格式。本文在上述工作的基础上将一种基于MUSCL方法的新型斜率限制器 [9] 加入熵稳定格式中，得到了一种新型的高分辨率熵稳定格式，该格式可以识别间断区域，并在间断区域自动添加耗散项，格式具有高分辨率且达到了二阶精度，显著改良了熵稳定格式在间断处抹平严重的现象，更加贴近精确解。

2. 一维浅水磁流体方程

平坦地形下的SWMHD方程可以写成如下的守恒形式

$\frac{\partial q}{\partial t}+\frac{\partial f\left(q\right)}{\partial x}=0$ (1)

式中，

$q=\left[\begin{array}{l}h\\ h{u}_1\\ h{u}_2\\ h{b}_1\\ h{b}_2\end{array}\right],f\left(q\right)=\left[\begin{array}{c}h{u}_1\\ h{u}_1^2+\frac{1}{2}g{h}^2-h{b}_1^2\\ h{u}_1{u}_2-h{b}_1{b}_2\\ 0\\ h{u}_1{b}_2-h{u}_2{b}_1\end{array}\right]$ . (2)

h为流体层深度， $u=\left(u_1,u_2\right)$ 为流体速度矢量， $b=\left(b_1,b_2\right)$ 为磁感应强度矢量，常数g为重力加速度。此外，与理想磁流体力学方程相似，方程(1)还必须满足无散条件： $\nabla \cdot hb=0$ ，在一维情况下，无散条件可以简化为： $\partial \left(h{b}_1\right)/\partial x=0$ 。

由于 $\nabla \cdot hb=0$ ，方程(1)为退化的，为了缓解这种退化本文采用文献 [17] 中的方法，将无散条件放宽为 $\nabla \cdot hb\approx 0$ 并添加了与无散条件成正比的源项 $s$ ，此时 $hb$ 为非常量。由此得到方程(1)的改进版本：

$\frac{\partial q}{\partial t}+\frac{\partial f}{\partial x}=\frac{\partial }{\partial t}\left[\begin{array}{l}h\\ h{u}_1\\ h{u}_2\\ h{b}_1\\ h{b}_2\end{array}\right]+\frac{\partial }{\partial x}\left[\begin{array}{c}h{u}_1\\ h{u}_1^2+\frac{1}{2}g{h}^2-h{b}_1^2\\ h{u}_1{u}_2-h{b}_1{b}_2\\ 0\\ h{u}_1{b}_2-h{u}_2{b}_1\end{array}\right]=-\frac{\partial \left(h{b}_1\right)}{\partial x}\left[\begin{array}{l}0\\ 0\\ 0\\ u_1\\ u_2\end{array}\right]=s$ . (3)

这里的源项 $s$ 与理想磁流体力学方程的Janhunen源项类似，源项 $s$ 的添加保持了浅水磁流体系统的质量守恒和动量守恒，也保证了黎曼问题的正定性。

改进后的一维浅水磁流体力学方程的通量函数 $f\left(q\right)$ 的Jacobi矩阵 $\frac{\partial f}{\partial q}$ 的特征值系统可参考文献 [17] 。SWMHD方程共有五条波，包含一条散度波，波速为 ${\lambda }_1=u$ ；两条Alfvén波，波速为 ${\lambda }_{2,3}=u\pm b_1$ ；两条磁重力波，波速为 ${\lambda }_{4,5}=u\pm c_g$ ；其中 $c_g=\sqrt{gh+b_1{}^2}$ 。相对应的Jacobi矩阵的右特征向量矩阵如下

$R=\left[\begin{array}{ccccc}1& 0& 1& 0& 1\\ u_1-c_g& 0& u_1& 0& u_1+c_g\\ u_2& 1& u_2& 1& u_2\\ 0& 0& c_g^2/b_1& 0& 0\\ b_2& 1& b_2& -1& b_2\end{array}\right]$ . (4)

同时，我们定义一维SWMHD方程的熵函数

$E\left(q\right)=\frac{1}{2}\left(g{h}^2+h{u}_1^2+h{u}_2^2+h{b}_1^2+h{b}_2^2\right)$ , (5)

和熵通量函数

$F\left(q\right)=g{h}^2{u}_1+\frac{1}{2}\left(h{u}_1^2+h{u}_2^2+h{b}_1^2+h{b}_2^2\right)-h{u}_2{b}_1{b}_2$ , (6)

熵变量

$v={\left(gh-\frac{1}{2}\left(u_1^2+u_2^2+b_1^2+b_2^2\right),u_1,u_2,b_1,b_2\right)}^{\text{T}}$ , (7)

熵势

$\phi =v^{\text{T}}f-F=\frac{1}{2}g{h}^2{u}_1-h{b}_1\left(u_1{b}_1+u_2{b}_2\right)$ . (8)

本文在空间方向上采用均匀网格上的半离散守恒型有限体积格式

$\frac{\text{d}}{\text{d}t}{q}_i\left(t\right)+\frac{1}{\Delta x}\left(f_{i+1/2}-f_{i-1/2}\right)=\frac{1}{2}\left(s_{i-1/2}+s_{i+1/2}\right)$ , (9)

时间方向上采用三阶Runge-Kutta [5] 方法：

$\{\begin{array}{l}{q}_i^{\ast }=q_i^k+\Delta tL\left(q_i^k\right)\hfill \\ q_i^{\ast \ast }=\frac{3}{4}{q}_i^k+\frac{1}{4}{q}_i^{\ast }+\frac{1}{4}\Delta tL\left(q_i^{\ast }\right)\hfill \\ q_i^{k+1}=\frac{1}{3}{q}_i^k+\frac{2}{3}{q}_i^{\ast \ast }+\frac{2}{3}\Delta tL\left(q_i^{\ast \ast }\right)\hfill \end{array}$ , (10)

其中 $L\left(q\right)=-\frac{1}{\Delta x}\left(f_{i+1/2}-f_{i-1/2}\right)+\frac{1}{2}\left(s_{i+1/2}-s_{i-1/2}\right)$ 。

3. 熵稳定格式

熵守恒格式要求保持总熵不变且满足离散熵等式 [19] ：

$\frac{\text{d}}{\text{d}t}E\left(q_i\left(t\right)\right)+\frac{1}{\Delta x}\left(F_{i+1/2}-F_{i-1/2}\right)=0$ , (11)

其中， $E\left(q_i\left(t\right)\right)$ 为凸的熵函数， $F_{i+1/2}$ 为熵通量函数，本文采用Winters和Gassner在文献 [17] 中的熵守恒格式数值通量，即

$f_{i+1/2}^c=\left[\begin{array}{c}\left\{h\right\}\left\{u_1\right\}\\ \left\{h\right\}{\left\{u_1\right\}}^2+\frac{1}{2}g{\left\{h\right\}}^2-\left\{h{b}_1\right\}\left\{b_1\right\}\\ \left\{h\right\}\left\{u_1\right\}\left\{u_2\right\}-\left\{h{b}_1\right\}\left\{b_2\right\}\\ \left\{h\right\}\left\{u_1\right\}\left\{b_1\right\}-\left\{h{b}_1\right\}\left\{u_1\right\}\\ \left\{h\right\}\left\{u_1\right\}\left\{b_2\right\}-\left\{h{b}_1\right\}\left\{u_2\right\}\end{array}\right]$ , (12)

与之对应的源项 $s$ 有以下的离散形式

$s_i=\frac{1}{2}\left(s_{i+1/2}+s_{i-1/2}\right)=-\frac{1}{2}\left({\left[h{b}_1\right]}_{i+1/2}{\left[\begin{array}{c}0\\ 0\\ 0\\ \begin{array}{l}\frac{\left\{u_1{b}_1\right\}}{\left\{\Delta x{b}_1\right\}}\\ \frac{\left\{u_2{b}_2\right\}}{\left\{\Delta x{b}_2\right\}}\end{array}\end{array}\right]}_{i+1/2}+{\left[h{b}_1\right]}_{i-1/2}{\left[\begin{array}{c}0\\ 0\\ 0\\ \begin{array}{l}\frac{\left\{u_1{b}_1\right\}}{\left\{\Delta x{b}_1\right\}}\\ \frac{\left\{u_2{b}_2\right\}}{\left\{\Delta x{b}_2\right\}}\end{array}\end{array}\right]}_{i-1/2}\right)$ . (13)

其中 ${\{\cdot \}}_{i+1/2}=\frac{{(\cdot )}_{i+1}+{(\cdot )}_i}{2}$ ， ${[\cdot ]}_{i+1/2}={(\cdot )}_{i+1}-{(\cdot )}_i$ ， ${\left[a^2\right]}_{i+1/2}=2{\left\{a\right\}}_{i+1/2}{\left[a\right]}_{i+1/2}$ ， ${\left[ab\right]}_{i+1/2}={\left\{a\right\}}_{i+1/2}{\left[b\right]}_{i+1/2}-{\left\{b\right\}}_{i+1/2}{\left[a\right]}_{i+1/2}$ 。

熵守恒格式在解的间断处表现出了严重的不稳定性，可以通过在熵守恒格式数值通量里添加适当的耗散项来控制这些不稳定因素，添加耗散项后的格式为熵稳定格式，本文采用文献 [17] 中的熵稳定格式数值通量

$f^{ES}=f^C-\frac{1}{2}R\left|\Lambda \right|S{R}^{\text{T}}\left[v\right]$ , (14)

其中， $R^{\text{T}}$ 为右特征向量矩阵 $R$ 的转置，特征值对角矩阵

$\Lambda =diag\left(u_1-c_g,u_1-b_1,u_1,u_1+b_1,u_1+c_g\right)$ , (15)

对角缩放矩阵

$T=diag\left(\frac{c}{c_g\sqrt{2g}},\frac{c}{\sqrt{2g}},\frac{b_1}{c_g\sqrt{g}},\frac{c}{\sqrt{2g}},\frac{c}{c_g\sqrt{2g}}\right)$ , (16)

为了方便书写引入了对角标度矩阵 $S=T^2$ ，上式中 $c=\sqrt{gh}$ 。

在离散情况下，(14)式还可以写成

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

4. 高分辨率熵稳定格式

熵守恒格式在解的光滑区域变现良好，但是在间断处有严重的伪振荡现象，而熵稳定格式在解的间断处能够有效地避免伪振荡的产生，但是一阶耗散项的加入使格式整体只有一阶精度。本节我们通过添加斜率限制器将熵守恒格式的数值通量(12)和熵稳定格式的数值通量(14)组合在一起，得到了一个新的熵稳定格式数值通量，所得到的新数值通量可以自动识别间断区域，在间断区域自动添加数值耗散项，有效地提高了格式的精度和分辨率。

4.1. 斜率限制器

已知 $q\left(x\right)$ 为SWMHD方程的精确解，在每个控制单元 $I_i=\left[x_{i-1/2},x_{i+1/2}\right]$ 上，精确解表示为 $q_i\left(x\right)$ ，我们对每个控制单元上的 $q_i\left(x\right)$ 进行重构。

不失一般性，假设对标量函数 $h_i\left(x\right)$ 进行了重构，将重构后的函数记为 $R_i\left(x\right)$ 。如果满足不等式 $\left|h_i\left(x\right)-R_i\left(x\right)\right|\le O\left({\left(\Delta x\right)}^{k+1}\right)$ ，那么就称这个重构函数 $R_i\left(x\right)$ 具有k-阶精度。在这里本文采用文献 [9] 中的MUSCL数据重构方法，选取

$R_i\left(x\right)=h_i^n+H\left(x-x_i\right)$ , (18)

其中

$H\left(x-x_i\right)={\frac{\partial h_i}{\partial x}|}_{x_i}\left(x-x_i\right)+\frac{1}{2}{\frac{{\partial }^2{h}_i}{\partial x^2}|}_{x_i}\left({\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}\frac{{\left(x-x_i\right)}^2}{2}\text{d}x}\right)$ , (19)

此时， $R_i\left(x\right)$ 具有2阶精度。为了计算方便，采用简单的中心差分来近似式(19)中的一阶和二阶偏导数，即

${\frac{\partial h_i}{\partial x}|}_{x_i}=\frac{h_{i+1}^n-h_{i-1}^n}{2\Delta x},{\frac{{\partial }^2{h}_i}{\partial x^2}|}_{x_i}=\frac{h_{i+1}^n-2h_i^n+h_{i-1}^n}{\Delta x^2}$ , (20)

再对重构后的 $R_i\left(x\right)$ 进行进一步的限制，

$R_i\left(x\right)=h_i^n+{\phi }_{i}H\left(x-x_i\right)$ , (21)

其中， ${\phi }_i\in \left[0,1\right]$ 称之为斜率限制器。本文选取文献 [9] 中的斜率限制器，即 ${\phi }_i=\max \left(\min \left({\phi }_L,{\phi }_R\right),0\right)$ 。其中，左斜率限制器

${\phi }_L=\min \left(1,\frac{h_{i-1}^n-h_i^n}{H\left(x_{i-1/2}-x_i\right)}\right)$ , (22)

右斜率限制器

${\phi }_R=\min \left(1,\frac{h_{i+1}^n-h_i^n}{H\left(x_{i+1/2}-x_i\right)}\right)$ , (23)

对于 ${\phi }_L$ 和 ${\phi }_R$ ，一般情况下，为了避免分母 $H\left(x_{i\pm 1/2}-x_i\right)$ 为0，可以在分母上加上一个很小的常数 $\epsilon ={10}^{-6}$ 。

4.2. 高分辨率熵稳定格式

本小节将3.1节提出的斜率限制器运用到第2节提出的熵稳定格式(14)上，得到了具有良好特性的高分辨率熵稳定格式。

按照3.1节的方法进行数据重构时，要先将守恒型变量 $q$ 转化为原始变量 $\tilde{q}$ ，即

$q=\left[\begin{array}{ccccc}h& h{u}_1& h{u}_2& h{b}_1& h{b}_2\end{array}\right]\to \left[\begin{array}{ccccc}h& u_1& u_2& b_1& b_2\end{array}\right]=\tilde{q}$ , (24)

再对原始变量 $\tilde{q}$ 进行上述的数据重构，重构完成后还要恢复变量值 $q$ ，将重构后的变量值代入熵稳定数值通量(14)就得到了具有高分辨率的熵稳定数值通量 $f_{i+1/2}^{ES2}$ ，即

$f^{ES2}={\tilde{f}}^C-\frac{1}{2}\tilde{R}\left|\tilde{\Lambda }\right|\tilde{S}{\tilde{R}}^{\text{T}}\left[\tilde{v}\right]$ , (25)

其中 ${\tilde{f}}_{i+1/2}^C$ ， $\tilde{R}$ ， $\tilde{\Lambda }$ ， $\tilde{S}$ ， $\tilde{v}$ 为 $f_{i+1/2}^C$ ， $R$ ， $\Lambda$ ， $S$ ， $v$ 重构后的值。

在离散情况下，(25)式可以写成

$f_{i+1/2}^{ES2}={\tilde{f}}_{i+1/2}^C-\frac{1}{2}{\tilde{R}}_{i+1/2}\left|{\tilde{\Lambda }}_{i+1/2}\right|{\tilde{S}}_{i+1/2}{\tilde{R}}_{i+1/2}^{\text{T}}{\left[\tilde{v}\right]}_{i+1/2}$ , (26)

5. 数值算例

本节对所得到的熵守恒格式(C，数值通量表达式为(12))，熵稳定格式(ES，数值通量表达式为(14))，高分辨率熵稳定格式(ES2，数值通量表达式为(25))进行了一系列的数值模拟，进行数值模拟时采用的工具为MATLAB软件。以下的一维数值算例中，条件数CFL = 0.5，均采用200个均匀网格点进行计算。参考解(Reference)均由CFL = 0.1时取5000个均匀网格点的熵稳定格式(ES)得到。

算例1：精度测试

此算例主要是用于测试ES2格式的精度，在计算区域 $\Omega =\left[-1,1\right]$ 上考虑如下初值问题

$h\left(x,0\right)=1$ ， $u_1\left(x,0\right)=0$ ， $u_2\left(x,0\right)=\sin \left(2\pi \left(x\right)\right)$ ， $b_1\left(x,0\right)=1$ ， $b_2\left(x,0\right)=v_2\left(x,0\right)$ ，

常数 $g=1$ ，采用周期性边界条件，计算到 $t=1$ 时刻。该算例的精确解为

$h\left(x,t\right)=1$ ， $u_1\left(x,t\right)=0$ ， $u_2\left(x,t\right)=\sin \left(2\pi \left(x+t\right)\right)$ ， $b_1\left(x,t\right)=1$ ， $b_2\left(x,t\right)=v_2\left(x,t\right)$ 。

见表1，表1展示了不同网格数下变量 $u_2$ 的数值误差以及收敛阶。从表1结果来看，ES2格式在L¹和L²范数意义下收敛阶均超过二阶，具有高精度特性。

Table 1. Numerical errors and order of convergence for variables u 2

表1. 变量 $u_2$ 的数值误差和收敛阶

算例2：强黎曼问题

在区域 $\Omega =\left[-1,1\right]$ 上求解如下初值问题

$\left[h,u_1,u_2,b_1,b_2\right]=\{\begin{array}{l}\left[1,0,0,1,0\right],\hfill \\ \left[2,0,0,0.5,1\right],\hfill \end{array}\begin{array}{c}x\le 0,\\ x>0.\end{array}$

其中 $g=1$ ，计算时采用Neumann边界条件， $t=0.4$ 时的数值解和参考解见图1。从图1中能够观察到，C格式在间断区域有明显的伪振荡现象，而ES格式和ES2格式都有效改善了伪振荡现象，同时ES2格式由于采取了二阶数据重构，能够更精确地捕捉到波的结构，跨越间断所需的单元个数也明显减少，显著提高了解的分辨率。

Figure 1. Numerical results of example 2

图1. 算例2的数值结果

算例3：弱黎曼问题

在区域 $\Omega =\left[-1,1\right]$ 上求解如下初值问题

$\left[h,u_1,u_2,b_1,b_2\right]=\{\begin{array}{l}\left[1,0,1,1,1\right],\hfill \\ \left[1-{10}^{-4},0,1,\left(1+{10}^{-4}\right)/\left(1-{10}^{-4}\right),1.0002\right],\hfill \end{array}\begin{array}{c}x\le 0,\\ x>0.\end{array}$

其中 $g=1$ ，计算时采用Neumann边界条件。变量 $h$ ， $u_1$ ， $u_2$ ， $b_1$ ， $b_2$ ， $h{u}_2$ ， $h{b}_2$ 在 $t=0.4$ 时的数值解和参考解见图2。观察图2可以发现C格式仍然有明显的振荡，ES格式虽然消除了振荡但是在间断处抹平现象严重，而ES2格式由于添加了斜率限制器的原因，不仅消除了振荡而且在每个间断处都能够准确地捕捉，在间断处更加贴近精确解，具有高分辨率这一良好特性。

Figure 2. Numerical results of example 3

图2. 算例3的数值结果

算例4： ${\left(h{b}_1\right)}_L={\left(h{b}_1\right)}_R$ 的黎曼问题

在区域 $\Omega =\left[0,1\right]$ 上求解如下初值问题

$\left[h,u_1,u_2,b_1,b_2\right]=\{\begin{array}{l}\left[1,0.2,0.7,0.5,0.4\right],\hfill \\ \left[0.5,-0.1,0.3,1,0.1\right],\hfill \end{array}\begin{array}{c}x\le 0.5,\\ x>\mathrm{0.5.}\end{array}$

其中 $g=9.81$ ，计算时采用Neumann边界条件。 $t=0.1$ 时的数值解和参考解见图3，此问题中 $h{b}_1=1/2$ 保持不变。从图3中可以发现ES2格式有效地改善了C格式抹平严重的现象，而且跨越间断所需的单元个数也明显减少，同时在间断区域也没有明显的伪振荡，展现了其高分辨率的优良特性。

Figure 3. Numerical results of example 4

图3. 算例4的数值结果

以上三个算例良好的数值结果说明了ES2格式在一维问题上的优良特性，我们进一步将所得到的ES2格式逐维推广到二维情况，对Rotor-like问题进行了数值模拟，验证了新格式在二维问题上的优良特性。

算例5：Rotor-like问题

此算例来自文献 [14] ，类似于MHD方程的Rotor问题 [20] ，计算区域 $\Omega =\left[-1,1\right]\times \left[-1,1\right]$ ，初始条件为

$\left[h,u_1,u_2,b_1,b_2\right]=\{\begin{array}{l}\left[10,-y,x,0.1,0\right],\hfill \\ \left[1,0,0,1,0\right],\hfill \end{array}\begin{array}{c}r\le 0.2,\\ r>\mathrm{0.2.}\end{array}$

其中， $r=\sqrt{x^2+y^2}$ ， $g=1$ ，最终模拟时间为 $t=0.15$ ，计算域被离散为100 × 100个网格点，采用40条等值线绘出，CFL = 0.5，边界条件仍然采用Neumann边界条件，数值结果见图4。图4给出了ES格式(右边)和ES2格式(左边)的数值模拟结果，本文由于对无散条件采用了与文献 [14] 不同的处理方法，所以变量 $b_1$ 和 $b_2$ 的结果略有不同。从计算结果看，ES2格式的图形整体形状更接近于圆形，中间区域的圆形轮廓也更加清晰，说明其对于边界处理得更好，除此之外，在变量h的模拟解中，我们发现ES2格式的中间区域也没有明显的杂乱线条，这与文献 [14] 中结果相当，由此充分展示了高分辨率熵稳定格式的优良性能。

Figure 4. Numerical results of example 5

图4. 算例5的数值结果

6. 结论

对于SWMHD方程，由于其并不是严格的双曲型方程，所以在构造熵格式时需要添加额外的源项并对其无散条件进行弱化，由此增加了熵格式求解的难度，本文在熵格式的基础上添加了一种斜率限制器，最终所得到的高分辨率熵稳定格式在一维和二维问题的数值求解过程中都体现出了其高精度，高分辨率，无振荡，鲁棒的优良特性，不失为求解SWMHD方程的一种较为理想的方法。另外，本文的研究对象主要是平坦地形下的SWMHD方程，对于非平坦地形也可以用类似的方法进行求解。

NOTES

^*通讯作者。