理想磁流体方程的MUSCL-型熵相容格式研究

The Research of MUSCL-Type Entropy Consistent Schemes for Ideal Magnetohydrodynamics Equations

摘要: 针对磁流体动力学方程，本文提出了一种基于MUSCL-Hancock方法求解MHD方程的熵相容格式(EC-MHM格式)，获得了一种求解理想磁流体动力学方程的高分辨率熵相容格式。该格式在解的光滑区域具有高精度；在解的间断区域可以合理地控制耗散，使抹平现象得到改善，还可有效避免非物理现象的产生。文中还证明了熵相容格式的收敛性。采用熵稳定格式、熵相容格式和新的高分辨率熵相容格式对理想磁流体动力学方程进行数值模拟。结果表明：新格式能准确地捕捉解的结构，且具有无振荡、高分辨、鲁棒等特性。

Abstract: Focusing on the idea magnetohydrodynamic (MHD) equations, this paper presents an entropy-consistent scheme based on the MUSCL-Hancock method for solving MHD equations, termed the Entropy Consistent MUSCL-Hancock (EC-MHM) scheme, and thus achieving a high-resolution entropy-consistent formulation for solving ideal MHD equations. This scheme exhibits high accuracy in smooth regions of the solution and effectively controls dissipation in discontinuous zones, leading to an improvement in the smearing phenomenon and efficiently preventing the emergence of non-physical oscillations. The convergence of the entropy consistent scheme is also proved. The ideal MHD equations are numerically simulated by entropy stable scheme, entropy consistent scheme and the new high resolution entropy consistent scheme. The results show that the new scheme can accurately capture the structure of the solution, and has the characteristics of no oscillation, high resolution and robustness.

1. 引言

磁流体动力学(Magnetohydrodynamics, MHD)是将经典流体力学和电动力学的方法结合起，针对导电流体和磁场之间的相互作用进行研究的一门学科。常应用于等离子体物理学、受控热核聚变、天体物理研究和电磁推进等领域。虽然MHD方程并不是严格的双曲型方程，但其波的特征结构与流体力学方程类似，因此可以将双曲型守恒律方程组的数值求解方法推广应用于MHD方程组的数值求解当中，以获得求解理想磁流体动力学方程的高效、稳定的数值方法。

对于一般双曲守恒律方程，“熵”一直是一个备受重视的量，它在解的光滑区域保持不变，但在激波等间断区域会有所增加(数学模型中表现为熵减，即熵耗散)，物理上对应着发生熵增，因此对于“熵”的研究是一个非常重要的课题。Lax在1954年提出了弱解的概念[1]，允许间断解存在，然而弱解并不惟一。随后又在文献[2]中提出熵稳定条件，并证明了满足熵稳定条件的Cauchy问题的数值解在任何时刻惟一且具有物理意义。1987年，Tadmor [3]构造了一类二阶的熵守恒格式，同时提出一种精确量化数值黏性的方法，即：比熵守恒格式含有更多黏性的三点格式是熵稳定的。2006年，Roe [4]提出在熵守恒格式的基础上添加Roe格式的数值黏性项，得到了一类熵稳定格式。该格式对激波具有良好的捕捉效果。近年来，国内外发展了一系列求解MHD方程的数值方法[5]-[8]，但使用熵稳定格式求解的方法很少。直至2016年，Winters等[9]提出了一种针对理想MHD方程的熵守恒通量，并加入适当的耗散得到熵稳定格式。然而该格式对于间断的捕捉效果并不理想，整体精度只有一阶。

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

2. 一维理想磁流体动力学方程

2.1. 一维理想磁流体动力学方程

磁流体动力学(Magnetohydrodynamics, MHD)方程是通过流体力学中的Navier-Stokes方程和电动力学中的Maxwell方程耦合而成，其写成守恒律的形式为

$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}=0$ , (1)

该式中，

$U=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho e\\ B_1\\ B_2\\ B_3\end{array}\right]$ , $F\left(U\right)=\left[\begin{array}{c}\rho u\\ \rho u^2+p+\frac{1}{2}{\Vert B\Vert }^2-B_1^2\\ \rho uv-B_1{B}_2\\ \rho uw-B_1{B}_3\\ u\left(\rho e+p+\frac{1}{2}{\Vert B\Vert }^2\right)-B_1\left(u\cdot B\right)\\ 0\\ u{B}_2-v{B}_1\\ u{B}_3-w{B}_1\end{array}\right]$ . (2)

$\rho$ 、p分别为流体密度和压力，$u=\left(u,v,w\right)$ 为流体速度矢量，$B=\left(B_1,B_2,B_3\right)$ 为磁感应强度矢量，$\rho e$ 表示总能量，即

$\rho e=\frac{p}{\gamma -1}+\frac{\rho }{2}{\Vert u\Vert }^2+\frac{1}{2}{\Vert B\Vert }^2$ , (3)

式中$\gamma$ 为比热比，${\Vert B\Vert }^2=B_1^2+B_2^2+B_3^2$ 。对于理想MHD方程而言，一个额外的无散度条件$\nabla \cdot B=0$ 需要被满足，否则在进行数值模拟时会出现数值不稳定[13]；又由于$B=\left(B_1,B_2,B_3\right)$ ，则根据散度的向量形式以及空间变量x可知，此时有$\partial B_1/\partial x=0$ 成立，即$B_1$ 为常数。

方程(1)的Jacobi矩阵$\partial F/\partial U$ 的特征值系统在文献[14]中有详细叙述：包含一条散度波，波速为${\lambda }_D=u$ ；一条熵波，波速为${\lambda }_E=u$ ；两条快波，波速为${\lambda }_{\pm f}=u\pm c_f$ ；两条慢波，波速为${\lambda }_{\pm s}=u\pm c_s$ ；两条Alfvén波，波速为${\lambda }_{\pm a}=u\pm c_a$ ，其中

$c_a^2=\frac{B_1^2}{\rho }$ , $c_{f,s}^2=\frac{1}{2}\left({\frac{\gamma \rho +\Vert B\Vert }{\rho }}^2\pm \sqrt{{\left({\frac{\gamma \rho +\Vert B\Vert }{\rho }}^2\right)}^2-4\frac{\gamma \rho B_1^2}{{\rho }^2}}\right)$ . (4)

与之对应的Jacobi矩阵的右特征向量矩阵为

$R=\left[\begin{array}{cccccccc}\rho & 0& \rho & 1& 0& \rho & 0& \rho \\ -c_f& 0& -c_s& 0& 0& c_s& 0& c_f\\ \frac{c_f{b}_1{b}_2}{c_f^2-b_1^2}& B_3& \frac{c_s{b}_1{b}_2}{c_S^2-b_1^2}& 0& 0& -\frac{c_s{b}_1{b}_2}{c_S^2-b_1^2}& -B_3& \frac{c_f{b}_1{b}_2}{c_f^2-b_1^2}\\ \frac{c_f{b}_1{b}_3}{c_f^2-b_1^2}& -B_2& \frac{c_s{b}_1{b}_3}{c_s^2-b_1^2}& 0& 0& -\frac{c_s{b}_1{b}_3}{c_s^2-b_1^2}& B_2& -\frac{c_f{b}_1{b}_3}{c_f^2-b_1^2}\\ \rho a^2& 0& \rho a^2& 0& 0& \rho a^2& 0& \rho a^2\\ 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\sqrt{\rho }{c}_f^2{b}_2}{c_f^2-b_1^2}& -\sqrt{\rho }{B}_3& \frac{\sqrt{\rho }{c}_s^2{b}_2}{c_s^2-b_1^2}& 0& 0& \frac{\sqrt{\rho }{c}_s^2{b}_2}{c_s^2-b_1^2}& -\sqrt{\rho }{B}_3& \frac{\sqrt{\rho }{c}_f^2{b}_2}{c_f^2-b_1^2}\\ \frac{\sqrt{\rho }{c}_f^2{b}_3}{c_f^2-b_1^2}& \sqrt{\rho }{B}_2& \frac{\sqrt{\rho }{c}_s^2{b}_3}{c_s^2-b_1^2}& 0& 0& \frac{\sqrt{\rho }{c}_s^2{b}_3}{c_s^2-b_1^2}& \sqrt{\rho }{B}_2& \frac{\sqrt{\rho }{c}_f^2{b}_3}{c_f^2-b_1^2}\end{array}\right]$ , (5)

式中$b=\left(b_1,b_2,b_3\right)=B/\sqrt{\rho }$ ，$a^2=\gamma p/\rho$ 。

根据文献[3]的理论，可以通过定义熵函数、熵通量函数和熵变量等，构造出保持熵守恒的数值格式；再根据[3] [4] [17]的理论，在熵守恒格式的基础上，适当加上一定量的数值耗散，即可得到熵稳定格式，从而更准确地模拟这些物理过程[3]。因此，为了方便一维MHD方程熵相容格式的构造，我们对一维MHD方程的熵函数、熵通量函数和熵变量进行定义。

定义熵函数$E\left(U\right)=-\rho s/\left(\gamma -1\right)$ ，熵通量函数$Q\left(U\right)=-\rho us/\left(\gamma -1\right)$ ，其中物理熵$s=\ln \left(p\right)-\gamma \ln \left(\rho \right)$ ，相应的熵变量$\nu$ 和熵势$\phi$ 为

$\nu =E^{\prime }\left(U\right)={\left[\frac{\gamma -s}{\gamma -1}-\frac{\rho {\Vert u\Vert }^2}{2p},\frac{\rho u}{p},\frac{\rho v}{p},\frac{\rho w}{p},-\frac{\rho }{p},\frac{\rho B_1}{p},\frac{\rho B_2}{p},\frac{\rho B_3}{p}\right]}^{\text{T}}$ , (6)

$\phi \left(\nu \right)=\nu \cdot F-Q=\rho u+\frac{\rho u{\Vert B\Vert }^2}{2p}-\frac{\rho B_1\left(u\cdot B\right)}{p}$ . (7)

2.2. 熵稳定条件

熵是热力学中表征物质状态的参量之一，其物理意义是系统混乱程度的度量。在双曲守恒律系统中，当解连续时，熵在整个系统中保持不变，称为熵守恒；但当出现间断解时，熵会增加，从而熵守恒格式的数值解会出现伪振荡，因此需要对总熵进行耗散。如果熵耗散过程满足熵稳定条件，则称它是熵稳定的。

设$\left(E,Q\right)$ 为双曲守恒律方程(1)的一组熵对，将(1)式的两端同时乘$E^{\prime }{\left(U\right)}^T$ ，得到

$E{\left(U\right)}_t+E^{\prime }{\left(U\right)}^T{F}^{\prime }\left(U\right)U_x=0$ , (8)

如果$Q^{\prime }{\left(U\right)}^T=E^{\prime }{\left(U\right)}^T{F}^{\prime }\left(U\right)$ ，则有

$E{\left(U\right)}_t+Q{\left(U\right)}_x=0$ , (9)

我们称(9)式为熵等式。

当间断发生时，熵会有所增加，熵等式(9)不再成立。为此，我们从弱解的角度出发，将$U^{\epsilon }$ 作为带有粘性机制的双曲守恒律方程[12]：

$U_t^{\epsilon }+F{\left(U^{\epsilon }\right)}_x=\epsilon U_{xx}^{\epsilon }$ , (10)

的极限解。假设$\left(E,Q\right)$ 为双曲守恒律方程(1)的一组熵对，对方程(10)的两端同时左乘$E^{\prime }{\left(U^{\epsilon }\right)}^T$ ，由于熵函数$E\left(U\right)$ 为凸函数[15]可以得到：

$\begin{array}{c}E{\left(U^{\epsilon }\right)}_t+Q{\left(U^{\epsilon }\right)}_x=\epsilon E^{\prime }\left(U^{\epsilon }\right)U_{xx}^{\epsilon }\\ =\epsilon \left(E{\left(U^{\epsilon }\right)}_{xx}-{\left(U_x^{\epsilon }\right)}^T{E}^{″}\left(U^{\epsilon }\right)U_x^{\epsilon }\right)\\ \le \epsilon E{\left(U^{\epsilon }\right)}_{xx}\end{array}$ , (11)

当$\epsilon \to 0$ 时，得到熵不等式

$E{\left(U\right)}_t+Q{\left(U\right)}_x\le 0$ , (12)

我们称上述不等式为熵稳定条件。

3. 熵相容格式

3.1. 熵相容格式

熵守恒格式，即保持总熵不变，在离散的情况下满足离散熵等式：

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

式中，$E\left(U_i\left(t\right)\right)$ 为凸的熵函数，$Q_{i+1/2}$ 为与熵通量函数Q相容的数值熵通量函数[3]。2015年，Winters等[16]在Tadmor提出的逐段分解显式构造方法的基础上，通过引入参数向量

$z=\left[z_1,z_2,\cdots ,z_8\right]=\left[\sqrt{\frac{\rho }{p}},\sqrt{\frac{\rho }{p}}u,\sqrt{\frac{\rho }{p}}v,\sqrt{\frac{\rho }{p}}w,\sqrt{\rho p},B_1,B_2,B_3\right]$ , (14)

和Ismail公式

$\widehat{p}=\left\{\left\{z_1\right\}\right\}{z}_5^{\ln }$ , ${\widehat{p}}_1=\frac{\left\{\left\{z_5\right\}\right\}}{\left\{\left\{z_1\right\}\right\}}$ , ${\widehat{p}}_2=\frac{\gamma +1}{2\gamma }\frac{z_5^{\ln }}{z_1^{\ln }}+\frac{\gamma -1}{2\gamma }\frac{\left\{\left\{z_5\right\}\right\}}{\left\{\left\{z_1\right\}\right\}}$ ,

${\widehat{u}}_1=\frac{\left\{\left\{z_2\right\}\right\}}{\left\{\left\{z_1\right\}\right\}}$ , ${\widehat{v}}_1=\frac{\left\{\left\{z_3\right\}\right\}}{\left\{\left\{z_1\right\}\right\}}$ , ${\widehat{w}}_1=\frac{\left\{\left\{z_4\right\}\right\}}{\left\{\left\{z_1\right\}\right\}}$ , ${\widehat{u}}_2=\frac{\left\{\left\{z_1{z}_2\right\}\right\}}{\left\{\left\{z_1^2\right\}\right\}}$ , ${\widehat{v}}_2=\frac{\left\{\left\{z_1{z}_3\right\}\right\}}{\left\{\left\{z_1^2\right\}\right\}}$ , ${\widehat{w}}_2=\frac{\left\{\left\{z_1{z}_4\right\}\right\}}{\left\{\left\{z_1^2\right\}\right\}}$ ,

${\widehat{B}}_1=\left\{\left\{z_6\right\}\right\}$ , ${\widehat{B}}_2=\left\{\left\{z_7\right\}\right\}$ , ${\widehat{B}}_3=\left\{\left\{z_8\right\}\right\}$ , ${\dot{B}}_1=\left\{\left\{z_6^2\right\}\right\}$ , ${\dot{B}}_2=\left\{\left\{z_7^2\right\}\right\}$ , ${\dot{B}}_3=\left\{\left\{z_8^2\right\}\right\}$ ,

$\widehat{B_1{B}_2}=\left\{\left\{z_6{z}_7\right\}\right\}$ , $\widehat{B_1{B}_3}=\left\{\left\{z_6{z}_8\right\}\right\}$ ,

推导出MHD方程的熵守恒通量：

$F_{i+1/2}^C=\left[\begin{array}{c}\widehat{\rho }{\widehat{u}}_1\\ {\widehat{p}}_1+\widehat{\rho }{\widehat{u}}_1^2+\frac{1}{2}\left({\dot{B}}_1+{\dot{B}}_2+{\dot{B}}_3\right)-{\dot{B}}_1\\ \widehat{\rho }{\widehat{u}}_1{\widehat{v}}_1-\widehat{B_1{B}_2}\\ \widehat{\rho }{\widehat{u}}_1{\widehat{w}}_1-\widehat{B_1{B}_3}\\ \frac{\gamma {\widehat{u}}_1{\widehat{p}}_2}{\gamma -1}+\frac{\widehat{\rho }{\widehat{u}}_1}{2}\left({\widehat{u}}_1^2+{\widehat{v}}_1^2+{\widehat{w}}_1^2\right)+{\widehat{u}}_2\left({\widehat{B}}_2^2+{\widehat{B}}_3^2\right)-{\widehat{B}}_1\left({\widehat{v}}_2{\widehat{B}}_2+{\widehat{w}}_2{\widehat{B}}_3\right)\\ 0\\ {\widehat{u}}_2{\widehat{B}}_2-{\widehat{v}}_2{\widehat{B}}_1\\ {\widehat{u}}_2{\widehat{B}}_3-{\widehat{w}}_2{\widehat{B}}_1\end{array}\right]$ , (15)

其中$\left\{\{·\}\right\}=\frac{{(·)}_{i+1}+{(·)}_i}{2}$ ，${(·)}^{\ln }=\frac{{(·)}_{i+1}-{(·)}_i}{\ln {(·)}_{i+1}-\ln {(·)}_i}$ 。

由于Euler方程的解在跨过激波时产生密度跳跃的立方级的熵增[17]，熵守恒格式在激波处会出现伪振荡现象。因此需要构造满足熵稳定条件的数值格式。首先用熵增推导出的耗散项去修正熵守恒格式，得到熵稳定格式；对熵稳定格式的熵增进一步精确估计，即可得到Ismail和Roe定义的熵相容格式[17]，即对熵增估计更精确的熵稳定格式格式。

对于双曲守恒律方程组来说，选择路径积分

$U=\left(1-\xi \right)U_i+\xi U_{i+1},0\le \xi \le 1$ , (16)

来连接左状态$U_i$ 和右状态$U_{i+1}$ 时，解在跨越激波时产生的熵增[17]为：

$\dot{E}=\frac{1}{2}{\left[\nu \right]}_{i+1/2}^{\text{T}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{-1}\text{d}\xi }{\left[U\right]}_{i+1/2}$, (17)

我们知道，熵增总是正的，即$\dot{E}>0$ 。由于熵变量$\nu =E^{\prime }\left(U\right)$ 且E为凸函数，则有

${\left[\nu \right]}_{i+1/2}=E^{″}\left(\xi \right){\left[U\right]}_{i+1/2}$ , (18)

式中$E^{″}\left(\xi \right)>0$ ，表明${\left[\nu \right]}_{i+1/2}$ 和${\left[U\right]}_{i+1/2}$ 具有相同的正负号，因此要使式(17)满足$\dot{E}>0$ ，只需保持中间矩阵${\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{-1}\text{d}\xi }$ 的正定性，为此本文做出如下修正：

$\dot{E}=\frac{1}{2}{\left[\nu \right]}_{i+1/2}^{\text{T}}\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{-1}\text{d}\xi }\right|{\left[U\right]}_{i+1/2}$, (19)

其对应的数值耗散为

$D^{EP}=\frac{1}{2}\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{-1}\text{d}\xi }\right|{\left[U\right]}_{i+1/2}$ , (20)

在计算过程中，由于矩阵积分${\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{-1}\text{d}\xi }$ 计算量大且形式复杂，为了简便计算，使用Gauss-Legendre积分公式进行近似，得到

$D^{EP}=\frac{1}{2}{R}_{i+1/2}\left|{\displaystyle \sum _{k=1}^n{\omega }_k\left(2S_k-1\right)\Lambda \left(S_k{U}_{i+1}+\left(1-S_k\right)U_i\right)}\right|R_{i+1/2}^{-1}{\left[U\right]}_{i+1/2}$ , (21)

其中$S_k$ 是高斯点，${\omega }_k$ 是与高斯点对应的系数，整个积分区间为$\left[0,1\right]$ 。通常情况下我们采用三点Gauss-Legendre积分公式进行积分，此时

$S_{1,3}=\frac{1}{2}\pm \frac{\sqrt{15}}{10}$ , $S_2=\frac{1}{2}$ , ${\omega }_{1,3}=\frac{5}{9}$ , ${\omega }_2=\frac{8}{9}$ . (22)

在求解理想MHD方程的过程中，将熵耗散项$D^{EP}$ 添加到熵稳定格式的数值通量[4]

$F_{i+1/2}^{ES}=F_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left|\Lambda \right|R_{i+1/2}^{-1}{\left[U\right]}_{i+1/2}$ , (23)

中，从而产生足够的熵耗散来抵消熵增，得到一种MHD方程的新的熵相容格式的数值通量：

$F_{i+1/2}^{EC}=F_{i+1/2}^C-\frac{1}{2}{R}_{i+1/2}\left(\left|\Lambda \right|+\left|{\displaystyle \sum _{k=1}^n{\omega }_k\left(2S_k-1\right)\Lambda \left(U\left(S_k\right)\right)}\right|\right)R_{i+1/2}^{-1}{\left[U\right]}_{i+1/2}$ , (24)

式中，特征值对角矩阵

$\Lambda =diag\left({\widehat{u}}_1-c_f,{\widehat{u}}_1-c_a,{\widehat{u}}_1-c_s,{\widehat{u}}_1,{\widehat{u}}_1,{\widehat{u}}_1+c_s,{\widehat{u}}_1+c_a,{\widehat{u}}_1+c_f\right)$ . (25)

$\left|\Lambda \right|=diag\left(\left|{\widehat{u}}_1-c_f\right|,\left|{\widehat{u}}_1-c_a\right|,\left|{\widehat{u}}_1-c_s\right|,\left|{\widehat{u}}_1\right|,\left|{\widehat{u}}_1\right|,\left|{\widehat{u}}_1+c_s\right|,\left|{\widehat{u}}_1+c_a\right|,\left|{\widehat{u}}_1+c_f\right|\right)$ ，$R_{i+1/2}$ 是与之对应的右特征向量矩阵，$R_{i+1/2}^{-1}$ 是矩阵$R_{i+1/2}$ 的逆矩阵。

3.2. 熵相容格式的稳定性论证

首先引入关于守恒型差分格式的如下的定理3.1。

定理3.1 [3] 若与$F$ 相容的数值通量$F_{i+1/2}\left(U_i,U_{i+1}\right)$ 关于$U_i$ ，$U_{i+1}$ Lipschits连续，那么一定存在一个矩阵$K_{i+1/2}$ 使得$F_{i+1/2}$ 可以写成如下形式：

$F_{i+1/2}=\frac{1}{2}\left(F\left(U_i\right)+F\left(U_{i+1}\right)\right)-\frac{1}{2}{K}_{i+1/2}\left(U_{i+1}-U_i\right)$ . (26)

由上述定理可知，守恒型熵守恒格式的数值通量$F_{i+1/2}^C$ 为

$F_{i+1/2}^C=\frac{1}{2}\left(F\left(U_i\right)+F\left(U_{i+1}\right)\right)-\frac{1}{2}{K}_{i+1/2}^C\left(U_{i+1}-U_i\right)$ , (27)

其中，$K_{i+1/2}^C$ 是与之相对应的数值粘性矩阵。

定理3.2熵相容格式满足熵稳定条件。

证明：本文中提到的熵格式都是守恒型三点格式，因此熵相容格式的数值通量可以写为带有粘性的形式：

$F_{i+1/2}^{EC}=\frac{1}{2}\left(F\left(U_i\right)+F\left(U_{i+1}\right)\right)-\frac{1}{2}{K}_{i+1/2}^{EC}\left(U_{i+1}-U_i\right)$ , (28)

其中$K_{i+1/2}^{EC}=K_{i+1/2}^C+R_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \left(U\left(\xi \right)\right)\text{d}\xi }\right|\right)R_{i+1/2}^{-1}$ ，根据Barth在文献[18]中证明的$U_v\approx R{R}^T$ ，$v$ 为熵变量$v=E^{\prime }\left(U\right)$ ，(28)式可以写为

$F_{i+1/2}^{EC}=\frac{1}{2}\left(F\left(U_i\right)+F\left(U_{i+1}\right)\right)-\frac{1}{2}{\tilde{K}}_{i+1/2}^{EC}\left(v_{i+1}-v_i\right)$ , (29)

其中${\tilde{K}}_{i+1/2}^{EC}=K_{i+1/2}^C+R_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \left(U\left(\xi \right)\right)\text{d}\xi }\right|\right)R_{i+1/2}^T$。

令$P_{i+1/2}^{EC}={\tilde{K}}_{i+1/2}^{EC}-K_{i+1/2}^C=R_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \left(U\left(\xi \right)\right)\text{d}\xi }\right|\right)R_{i+1/2}^T$，则有：

$\begin{array}{l}{\left(v_{i+1}-v_i\right)}_{i+1/2}^{{\rm T}}{P}_{i+1/2}^{EC}\left(v_{i+1}-v_i\right)\\ ={\left(v_{i+1}-v_i\right)}_{i+1/2}^{{\rm T}}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \left(U\left(\xi \right)\right)\text{d}\xi }\right|\right)R_{i+1/2}^T\left(v_{i+1}-v_i\right)\\ ={\left(R_{i+1/2}^T\left(v_{i+1}-v_i\right)\right)}^{{\rm T}}\left(\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \left(U\left(\xi \right)\right)\text{d}\xi }\right|\right)R_{i+1/2}^T\left(v_{i+1}-v_i\right)\end{array}$ , (30)

我们知道，$R_{i+1/2}^T$ 为可逆矩阵，$\left(v_{i+1}-v_i\right)$ 为非零向量，$\left|{\Lambda }_{i+1/2}\right|+\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \left(U\left(\xi \right)\right)\text{d}\xi }\right|$ 为对角矩阵。

因此，由(30)式可以得到${\left(v_{i+1}-v_i\right)}_{i+1/2}^{{\rm T}}{P}_{i+1/2}^{EC}\left(v_{i+1}-v_i\right)>0$ ，也就是${\left(v_{i+1}-v_i\right)}_{i+1/2}^{{\rm T}}{\tilde{K}}_{i+1/2}^{EC}\left(v_{i+1}-v_i\right)>{\left(v_{i+1}-v_i\right)}_{i+1/2}^{{\rm T}}{K}_{i+1/2}^C\left(v_{i+1}-v_i\right)$ ，这表明熵相容格式比熵守恒格式包含更多的数值粘性。由定理3.1可知，熵相容格式满足熵稳定条件。

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

采用MUSCL-Hancock方法[12]来更新下一时间层的物理量，是一个全离散过程。计算过程如下：

第一步：在$t_n$ 层利用MUSCL重构方法，重构$U_i\left(x\right)$ 在单元交界面处的左右极限值，详见[12]，结果为$U_i^L$ ，$U_i^R$ ；

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

${\overline{U}}_i^L=U_i^L-\frac{1}{2}\frac{\Delta t}{\Delta x}\left[f\left(U_i^R\right)-f\left(U_i^L\right)\right]$ , ${\overline{U}}_i^R=U_i^R-\frac{1}{2}\frac{\Delta t}{\Delta x}\left[f\left(U_i^R\right)-f\left(U_i^L\right)\right]$ . (31)

第三步：利用上述所得边界外推值，得出熵相容数值通量$F_{i+1/2}^{EC-MHM}$ ：

$F_{i+1/2}^{EC-MHM}={\overline{F}}_{i+1/2}^C-\frac{1}{2}{\overline{R}}_{i+1/2}\left(\left|\overline{\Lambda }\right|+\left|{\displaystyle \sum _{k=1}^n{\omega }_k\left(2S_k-1\right)\overline{\Lambda }\left(\overline{U}\left(S_k\right)\right)}\right|\right){\overline{R}}_{i+1/2}^{-1}{\left[\overline{U}\right]}_{i+1/2}$ , (32)

其中“$\overline{·}$ ”是重构后的取值。

第四步：利用上述得到的熵相容通量$F_{i+1/2}^{EC-MHM}$ 代入以下守恒型格式：

$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)$ , (33)

更新$t_{n+1}$ 时刻的单元平均值。我们把格式(33)称为高分辨率熵相容格式。

实际重构的过程是：将守恒型变量$U$ 转换为原始变量，即

$\left[\rho ,\rho u,\rho v,\rho w,\rho e,B_1,B_2,B_3\right]\to \left[\rho ,u,v,w,p,B_1,B_2,B_3\right]$ , (34)

对原始变量进行上述的数据重构，重构完成后，恢复变量值，将重构后的变量值代入数值通量中再进行计算，从而得到具有高分辨率的熵相容格式。

5. 数值算例

对熵稳定格式(ES) (数值通量表达式为式(23))，熵相容格式(EC) (数值通量表达式为式(24))和高分辨率熵相容格式(EC-MHM) (数值通量表达式为式(32))进行数值模拟，对比说明三种格式的表现。除算例5.1.和5.5.在计算时采用周期性边界条件外，其余的数值算例均采用Neumann边界条件，Ref.代表参考解，由选取5000个均匀网格的熵稳定格式[19]计算得到。

5.1. 一维光滑Alfvén波问题

光滑Alfvén波算例经常用来计算理想MHD方程数值格式的精度。本文考虑一维情况下的光滑Alfvén波的如下初值问题：

$\rho \left(x,0\right)=1$ , $u\left(x,0\right)=0$ , $v\left(x,0\right)=0.1\sin \left(2\pi x\right)$ , $w\left(x,0\right)=1$ , $p\left(x,0\right)=0.1$ , $B_1\left(x,0\right)=1$ , $B_2\left(x,0\right)=v\left(x,0\right)$ , $B_3\left(x,0\right)=0.1\cos \left(2\pi x\right)$ .

该算例的精确解为：

$\rho \left(x,t\right)=1$ , $u\left(x,t\right)=0$ , $v\left(x,t\right)=0.1\sin \left(2\pi \left(x+t\right)\right)$ , $w\left(x,t\right)=1$ , $p\left(x,t\right)=0.1$ , $B_1\left(x,t\right)=1$ , $B_2\left(x,t\right)=v\left(x,t\right)$ , $B_3\left(x,t\right)=0.1\cos \left(2\pi \left(x+t\right)\right)$ .

在计算时采用周期性边界条件，计算区间为$\left[0,1\right]$ ，绝热指数$\gamma =5/3$ 。

见表1展示了计算该算例到$t=5$ 时刻的不同网格数下变量v的数值误差以及收敛阶。见表1结果来看，EC-MHM格式在L¹、L²范数意义下收敛阶均超过二阶，具有高精度特性。

Table 1. L¹, L² errors in v at $t=5$ and corresponding convergence rates for different mesh numbers

表1. $t=5$ 时不同网格数下v的L¹，L²误差以及对应的收敛阶

网格数	L^1范数	收敛阶	L^2范数	收敛阶
25	3.6000e−2	−	1.3900e−2	−
50	7.7500e−3	2.2157	2.4365e−3	2.5122
100	1.6518e−3	2.2302	4.1581e−4	2.5508
200	3.5292e−4	2.2266	6.5096e−5	2.6753
400	6.8964e−5	2.3554	1.1032e−5	2.5609
800	1.2976e−5	2.4100	1.7013e−6	2.6970

5.2. Brio-Wu激波管问题

控制方程为式(1)、(2)，在计算区域$\Omega =\left[-1,1\right]$ 上给定初始条件

$\left[\rho ,u,v,w,p,B_1,B_2,B_3\right]=\{\begin{array}{l}\left[1,0,0,0,1,0.75,1,0\right],x\le 0.5,\\ \left[0.125,0,0,0,0.1,0.75,-1,0\right],x>\mathrm{0.5.}\end{array}$

其中$\gamma =5/3$ ，空间网格数为 500，计算至$t=0.12$ 。计算结果与参考解见图1。将EC-MHM格式与ES格式和EC格式进行比较，可发现EC-MHM格式对每个间断都能准确地捕捉，跨越间断所需的单元个数也明显减少，提高了解的分辨率。从图中可以看出，ES格式和EC格式的总熵耗散大致相同，EC-MHM格式相较于ES格式而言熵耗散更小，既有效地抑制抹平现象，还避免了间断附近的伪振荡，数值解更加贴近参考解。

Figure 1. Calculation results of the Brio-Wu shock tube problem

图1. Brio-Wu激波管问题的计算结果

5.3. Torrihon黎曼问题

控制方程为式(1)、(2)，在计算区域$\Omega =\left[-1,1.5\right]$ 上给定初始条件

$\left[\rho ,u,v,w,p,B_1,B_2,B_3\right]=\{\begin{array}{l}\left[3,0,0,0,1,1.5,1,0\right],x\le 0;\\ \left[1,0,0,0,1,1.5,\cos \left(1.5\right),\sin \left(1.5\right)\right],x>0.\end{array}$

其中$\gamma =5/3$ ，空间网格数为500，计算至$t=0.4$ 。计算结果与精确解见图2。可以看出，相较ES格式和EC格式而言，EC-MHM格式能够准确地捕捉到解的结构，解的分辨率明显提高，使得数值结果更加贴近参考解。通过总熵的变化来看，随着时间的推移，ES格式和EC格式的耗散较多，表明其数值结果在间断处抹平现象严重；相较而言，EC-MHM格式对总熵耗散控制地更加合理，说明其数值结果对解的捕捉更加锐利。

Figure 2. Calculation results of the Torrihon Riemann problem

图2. Torrihon黎曼问题的计算结果

5.4. Ryu and Jones黎曼问题

控制方程为式(1)、(2)，在计算区域$\Omega =\left[-1,1\right]$ 上给定初始条件

$\left[\rho ,u,v,w,p,B_1,B_2,B_3\right]=\{\begin{array}{l}\left[1,0,0,0,1,0.7,0,0\right],x\le 0,\\ \left[0.3,0,0,1,0.2,0.7,1,0\right],x>0.\end{array}$

其中$\gamma =2$ ，空间网格数为500，计算至$t=0.4$ 。计算结果与精确解见图3。此问题的解包括：稀疏波、慢速稀疏波、接触间断、慢速激波和快速稀疏波。结果表明，ES格式的计算结果和EC格式的计算结果几乎相同，在间断处的抹平现象比较严重；而EC-MHM格式对每个间断捕捉更加准确，跨越间断时所需的单元个数也明显减少。从总熵图来看，ES格式和EC格式的总熵耗散大致相同，EC-MHM格式的总熵耗散比EC格式小，对总熵耗散控制的更加合理，有效改善了严重的抹平现象，从侧面说明新构造的高分辨率熵相容格式的优良性质。

5.5. First Rotor问题

计算区域$\Omega =\left[0,1\right]\times \left[0,1\right]$ ，$\gamma =1.4$ ，记$r=\sqrt{{\left(x-0.5\right)}^2+{\left(y-0.5\right)}^2}$ ，$r_0=0.1$ ，$r_1=0.115$ ，$u_0=2$ ，初始条件为：

$\{\begin{array}{l}\rho =10,\left(u,v\right)=\frac{u_0}{r_0}\left(-\left[y-0.5\right],\left[x-0.5\right]\right),r<r_0,\\ \rho =1+9f,\left(u,v\right)=\frac{f{u}_0}{r_0}\left(-\left[y-0.5\right],\left[x-0.5\right]\right),f=\frac{r_1-r}{r_1-r_0},r_0<r<r_1,\\ \rho =1,\left(u,v\right)=\left(0,0\right),r>r_1,\\ w=0,p=1,B=5/\sqrt{4\pi }{\left(1,0,0\right)}^{{\rm T}}.\end{array}$

Figure 3. Calculation results of the Ryu and Jones Riemann problem

图3. Ryu and Jones黎曼问题的计算结果

在本算例中选用周期性边界条件，计算至$t=0.15$ ，空间网格数为200 × 200。采用逐维推广的方法将熵相容格式和高分辨率熵相容格式分别应用于First Rotor问题数值模拟中，计算结果的云图见图4和图5。其中$Ma=\Vert u\Vert /a$ 表示马赫数，$a=\sqrt{\gamma p/\rho }$ 。

Figure 4. The result of the entropy stable scheme calculation of the First Rotor problem

图4. First Rotor问题的熵稳定格式计算结果

Figure 5. The result of the high-resolution entropy consistent scheme calculation of the First Rotor problem

图5. First Rotor问题的高分辨率熵相容格式计算结果

见图4与图5比较可知，高分辨率熵相容格式对激波等间断的捕捉效果较为清晰锐利；验证了运用新型斜率限制器得到的高分辨率熵相容格式具有一定的优良特性。

6. 结论

由于理想磁流体动力学方程波的特征结构与流体力学方程类似，因此将求解一般双曲守恒律的方法推广应用于理想磁流体方程中，在此基础上加入新构造的MUSCL-型斜率限制器，并对其采用MUSCL-Hancock方法来更新下一时间层的物理量，获得了理想磁流体动力学方程的高分辨率熵相容格式。数值结果表明，该格式具有高精度、无振荡、高分辨率、鲁棒等特性，是模拟理想磁流体方程较为理想的方法。

NOTES

^*通讯作者。