0.   引言

磁流体动力学(magnetohydrodynamic, MHD)方程作为天体物理学和工程中重要的物理模型，由Navier-Stokes方程和Maxwell方程耦合而成，可以很好地模拟空间等离子体的运动规律.MHD方程数值计算的难点之一是必须满足一个附加的磁场无散度条件，否则会导致在计算中出现负密度和负压.Godunov^[1]从方程的对称性出发，提出了通过添加一个非守恒源项来控制磁场散度，最终得到了修正后的MHD方程.

MHD方程具有双曲守恒律方程的性质，因此考虑将求解双曲守恒律方程的方法应用到MHD方程中.双曲守恒律方程的数值计算离不开熵.Tadmor^[2-3]首次引入熵守恒通量，添加合适的耗散算子来构造求解非线性双曲守恒律的熵稳定格式.熵稳定格式满足离散熵不等式，在间断区域可以保持熵的耗散，保证不会出现伪振荡.Roe^[4]提出将Roe数值格式中的黏性项作为熵稳定格式中的耗散项，进而得到Roe型的熵稳定格式.郑素佩等^[5-7]利用旋转不变性得到混合数值通量并用来求解浅水波方程以及Euler方程，在保持熵稳定性质的同时又提高了数值格式的分辨率.Fjordholm等^[8]为双曲守恒律方程设计了一种基于ENO(essentially non-oscillatory)重构的任意高阶精确的熵稳定格式，称之为TeCNO格式. 为了构建保持熵稳定并且精度更高的重构方式，Biswas等^[9]通过在耗散项部分添加开关算子，得到了一种基于WENO重构的熵稳定通量.Duan和Tang^[10-11]构造了基于该重构的高阶精度熵稳定有限差分格式，用来求解浅水磁流体方程和相对论流体方程.郑素佩等^[12]利用Lagrange多项式证明了一种满足符号性质的WENO重构，并构造了熵稳定格式.Duan和Tang^[13]将自适应移动网格与高阶熵稳定格式结合用来求解相对论流体方程.

国内外有许多学者研究将熵稳定格式应用于MHD方程求解中.Winters和Gassner^[14]针对于理想MHD方程构造了一组方便快捷的熵守恒通量，并引入了Janhunen源项，最终得到了求解MHD方程的一阶熵稳定格式.为了提高精度，Chandrashekar和Klingenberg^[15]基于TeCNO格式的构造思想，利用满足符号性质的Minmod限制器对特征熵变量进行重构，得到二阶熵稳定格式.翟梦情等^[16]将移动网格与熵稳定格式相结合，并用来求解一维理想MHD方程.Liu和Shu等^[17]将满足熵稳定性质的DG(discontinuous Galerkin)方法推广至MHD方程，并达到三阶精度.熵稳定格式的精度主要与耗散项中熵变量值的跳跃有关.通过对耗散项进行WENO重构，可以显著提高熵稳定格式的精度.但目前还没有将WENO重构与熵稳定格式结合用来求解理想磁流体方程的研究.

本文考虑利用五阶WENO重构的思想对特征熵变量进行重构，并且将高阶熵守恒通量与离散源项部分结合，通过在耗散项中加入开关函数，使得耗散项部分仅在间断位置被激活，从而得到一种低耗散、高分辨率的用来求解理想MHD方程的四阶熵稳定格式.

1.   控制方程

理想磁流体方程是由质量、动量、能量以及磁场所构成的一组双曲守恒律方程组，可表示为如下形式:

其中ρ为质量，$\rho \boldsymbol{u}=\left[\begin{array}{lll}\rho u & \rho v & \rho w\end{array}\right]^{\mathrm{T}}$为动量，ρe为总能量，B=[B₁  B₂  B₃]^T为磁场，p为压力，符号⊗表示张量积.压力p可由其他变量表示为

其中γ为绝热指数，在数值计算中被考虑为常数.记守恒变量U=[ρ  ρu  ρe  B]^T.对于理想MHD方程而言，一个额外的无散度条件

需要被满足，否则在进行数值模拟时会出现数值不稳定.给定第i个维度的熵对，即熵函数和熵通量函数:

其中s=ln(pρ^－γ), u_i为u的第i个分量.基于文献[12]中给出的熵变量的定义，可以得到

其中β=ρ/(2p).Godunov^[1]针对于理想MHD方程的无散度条件，提出通过添加一个额外的非守恒源项来控制磁场散度.最终得到修正后的MHD方程为

其中$\phi^{\prime}(\boldsymbol{V})=\left[\begin{array}{llll}0 & \boldsymbol{B}^{\mathrm{T}} & \boldsymbol{u} \cdot \boldsymbol{B} & \boldsymbol{u}^{\mathrm{T}}\end{array}\right]$.

2.   理想MHD方程四阶熵稳定格式

本节考虑在高阶熵守恒通量以及高阶离散源项的基础上，通过添加具有开关函数并且基于WENO重构的耗散项，以确保仅在满足符号性质的位置进行对特征熵变量的WENO重构，从而得到一种耗散更低并且满足熵稳定性质的高阶数值通量.

在一维情况下，式(6)可以表示为如下守恒形式:

其中守恒变量U=[ρ  ρu  ρe  B]^T，B₁(U)为磁场B的第一个分量，通量函数为

考虑在均匀划分的空间网格上对式(7)进行离散，得如下半离散有限体积格式：

其中$\boldsymbol{U}_{i}=\frac{1}{\left|C_{i}\right|} \int_{C_{i}} \boldsymbol{U}(x, t) \mathrm{d} x$为守恒变量U在第i个单元$C_{i}=\left[x_{x-1 / 2}, x_{x+1 / 2}\right]$上的单元平均值，$\hat{\boldsymbol{F}}_{x+1 / 2}$为与通量函数相容的数值通量.采用中心格式的磁场散度部分满足

2.1   熵稳定格式

若数值格式(8)满足离散熵等式则为熵守恒格式，即

其中$\hat{q}_{x+1 / 2}$为数值熵通量函数.Chandrashekar等^[15]提出了数值通量为熵守恒通量的一个充分条件，并给出了如下熵守恒通量：

其中

式中$\bar{a}_{x+1 / 2}=\frac{1}{2}\left(a_{i}+a_{i+1}\right), \tilde{a}_{x+1 / 2}=\frac{a_{i+1}-a_{i}}{\ln \left(a_{i+1}\right)-\ln \left(a_{i}\right)}$.对离散源项部分，取$\hat{B}_{1, i+1 / 2}=\bar{B}_{1, i+1 / 2}, \hat{B}_{1, i-1 / 2}=\bar{B}_{1, i-1 / 2}$.

熵守恒数值格式具有二阶精度，为了提高熵守恒格式的精度，Fjordholm等^[8]通过对二阶熵守恒通量进行线性组合，得到任意偶数阶精度的熵守恒通量.对于离散源项部分，为了获得与高阶熵守恒通量一致的精度，同样考虑线性组合，得到任意偶数阶精度的离散磁场分量.本文考虑四阶精度的数值通量和离散磁场分量，具体表示为

熵守恒格式在光滑区域可以保持高精度，但是当解出现间断时会产生熵耗散，仅用熵守恒格式会产生伪振荡.因此需要在熵守恒格式中添加额外的数值耗散机制以保证格式的熵稳定性.若数值格式满足离散熵不等式：

即为熵稳定格式.本文考虑将Roe格式中的数值黏性作为耗散项添加至熵守恒通量来构造熵稳定数值通量：

其中R_x+1/2为Jacobi矩阵$\frac{\partial}{\partial \boldsymbol{U}}\left(\boldsymbol{F}(\boldsymbol{U})+\phi^{\prime}(\boldsymbol{V}) B_{1}(\boldsymbol{U})\right)$的右特征向量的离散形式.|Λ_x+1/2|可以选择Roe型的耗散项

或者Lax-Friedrichs型的耗散项

其中$\left\{\lambda_{i+1 / 2, k}\right\}_{k=1, 2, \cdots, 8}$为Jacobi矩阵$\frac{\partial}{\partial \boldsymbol{U}}\left(\boldsymbol{F}(\boldsymbol{U})+\boldsymbol{\phi}^{\prime}(\boldsymbol{V}) B_{1}(\boldsymbol{U})\right)$的特征值，具体表示参见附录A.

2.2   五阶WENO重构

熵稳定通量的耗散项部分包含熵变量的差分，使得熵稳定格式精度只有一阶.为了得到高阶熵稳定格式，本文考虑五阶WENO重构方法.基于2.1小节中提出的四阶熵守恒通量，添加重构后的耗散项得到高阶熵稳定数值通量：

其中v_x+1/2^－和v_x+1/2⁺分别为特征熵变量v=R_x+1/2^TV在单元交界面i+1/2左右的WENO重构值.

WENO重构后v_x+1/2^－的第j个分量表示为v_i+1/2，j^－, 表达式如下：

其中v_{k, j}(k=i－2, i－1, i, i+1, i+2)表示第k个单元上的特征熵变量v_k的第j个分量.式(15)中的加权系数ω_k具体可表示为

其中ε=10^－6，用来防止出现分母部分为0；$d_{1}=\frac{1}{10}, d_{2}=\frac{3}{5}, d_{1}=\frac{3}{10} ; \beta_{k}$为光滑因子，采用Shu^[18]提出的计算方法，即

v_{i+1/2, j}⁺与v_{i+1/2, j}^－计算方法相同，只需要取与v_{i+1/2, j}^－对称的系数即可.

由于WENO重构过程不一定满足符号性质，即

符号性质保证了重构后数值格式仍然具有熵稳定性.Biswas等^[9]针对于WENO重构提出并证明了通过在每个特征熵变量分量的跳跃处添加一个开关函数：

即可保证经WENO重构后的数值通量仍然保持符号性质.由于最终得到的高阶熵稳定数值通量为

其中〈v〉_x+1/2的第k个分量〈v_k〉_x+1/2=S_k(v_{i+1/2, k}⁺－v_{i+1/2, k}^－).由式(18)可以看出，在WENO重构不满足符号性质的位置，对应耗散项部分变为0.与基于ENO重构的TeCNO格式相比，由于上述格式仅在符号性质满足的区域进行WENO重构，因此具有更小的耗散.有关ES-WENO格式的熵稳定性证明参见附录B.

2.3   二维理想MHD方程四阶熵稳定格式

二维理想MHD方程在网格划分均匀的区域上的有限体积半离散格式为

其中$\boldsymbol{U}_{i j}=\frac{1}{\left|C_{i j}\right|} \iint_{C_{i j}} \boldsymbol{U}(x, y, t) \mathrm{d} x \mathrm{~d} y$表示守恒变量U在均匀四面体计算单元C_ij=[x_x-1/2, x_x+1/2]×[y_j－1/2, y_j+1/2]上的单元平均值.$\hat{\boldsymbol{F}}_{i \pm 1 / 2, j}, \hat{\boldsymbol{G}}_{i, j \pm 1 / 2}$分别为x, y方向上相容的数值通量，$\hat{B}_{1, i \pm 1 / 2, j}, \hat{B}_{2, i, j \pm 1 / 2}$为离散的x, y方向上的磁场分量，并且满足

鉴于x方向上的数值通量表示与一维情况类似，本小节只考虑y方向.首先给出熵守恒数值通量：

其中

对y方向上的离散源项部分，考虑$\hat{B}_{2, i, j+1 / 2}=\bar{B}_{2, i, j+1 / 2}, \hat{B}_{2, i, j-1 / 2}=\bar{B}_{2, i, j-1 / 2}$.四阶熵守恒通量和四阶离散源项部分采用与式(11)相同的表示方式.最终得到y方向上高阶熵稳定通量如下：

其中R_{i, j+1/2}与Λ_{i, j+1/2}具体表达参考附录A，并取n=[0, 1, 0]^T; 〈v〉_{i, j+1/2}为带有开关函数的特征熵变量的跳跃值.

3.   数值算例

本节给出几个理想MHD方程数值算例用以检验该数值格式的精度、低耗性以及捕捉间断的能力.时间离散部分采用四阶精确的Runge-Kutta方法^[19].空间离散考虑四阶精度的熵守恒通量以及五阶WENO重构的耗散项.格式理论上整体精度可以达到四阶.二维算例网格数均选择256×256.如无特殊声明，耗散项部分均选择Roe型耗散项.一维算例中ES-WENO、TeCNO、ES和ref分别代表基于WENO重构的熵稳定格式、基于ENO重构的熵稳定格式、一阶熵稳定格式和参考解，其中参考解由5 000个网格的一阶熵稳定格式求得.

3.1   一维光滑Alfvén波问题

光滑Alfvén波算例经常用来计算理想MHD方程数值格式的精度.本文考虑一维情况下的光滑Alfvén波，用来验证ES-WENO格式高阶精度的特性，选择周期性边界条件并引入垂直磁场B_⊥=0.1sin(2πx)，平行磁场B_//=1.0以及z方向上的磁场B_z=0.1cos(2πx).具体的初始条件数据由下式给出：

计算区间为[0, 1]，绝热指数γ=5/3.由于周期性，解会在1个单位时间内恢复到初始状态.表 1给出了计算该算例到T=5时刻的L¹, L^∞误差以及对应的收敛阶.可以发现，本文构造的ES-WENO格式基本可以达到理论上的四阶精度.

表  1  T=5时不同网格数下B₁的L¹, L^∞误差以及对应的收敛阶

Table  1.  L¹, L^∞ errors in B₁ at T=5 and corresponding convergence rates for different mesh numbers

N	L1 error	order	L∞ error	order
16	9.165E-4		1.477E-3
32	2.838E-5	5.013	4.514E-5	5.032
64	2.100E-6	3.756	3.325E-6	3.763
128	1.320E-7	3.992	2.076E-7	4.001
256	8.337E-9	3.985	1.312E-8	3.984

下载: 导出CSV 

| 显示表格

3.2   Ryu-Jones Riemman问题

该一维算例用来测试ES-WENO格式低耗散的特性.计算区间为[－1, 1].初始条件如下：

绝热指数γ=5/3.选择Neumann边界条件.在200个网格下分别用ES-WENO和TeCNO格式计算该算例至T=0.4, 如图 1(a)所示，其中ES-WENO和TeCNO格式均选择耗散更大的Lax-Friedrichs型的耗散项用来抑制间断处出现的局部振荡.由图 1(b)可知，ES-WENO格式比TeCNO格式具有更小的耗散量，因此具有更高的分辨率.图 2给出了ES-WENO格式对应的开关函数在初始时刻以及终止时刻在单元界面上的函数值，与TeCNO格式不同，ES-WENO格式的耗散项只有在间断处及其邻近区域不为0.

图  1  Ryu-Jones Riemman问题

Figure  1.  The Ryu-Jones Riemman problem

下载: 全尺寸图片 幻灯片

图  2  ES-WENO格式在T=0和T=0.4对应的开关函数值

Figure  2.  The switch function values corresponding to the ES-WENO scheme at T=0 and T=0.4

下载: 全尺寸图片 幻灯片

3.3   Torrilhon Riemman问题

该一维算例用来测试ES-WENO格式的高分辨率特性.计算区域为[-1, 1.5].初始条件如下：

绝热指数γ=5/3.选择Neumann边界条件.ES-WENO格式选择Lax-Friedrichs型的耗散项.图 3给出了在400个网格下，该算例在T=0.4时刻的计算结果.从图中可以看到，接触间断只有在密度场中出现，其余场中依次出现的波形分别为快稀疏波、慢稀疏波、快激波和慢激波.熵稳定格式在间断处出现严重的抹平现象，并且在部分极值处丢失精度，新格式精确识别出了每一种波形，在极值处依然保持高精度并且没有产生伪振荡.

图  3  Torrilhon Riemman问题

Figure  3.  The Torrilhon Riemman problem

下载: 全尺寸图片 幻灯片

3.4   二维Orszag-Tang漩涡问题

该算例首先由Orszag和Tang^[20]提出，并且用来验证ES-WENO格式的数值稳定性.光滑的初始条件由下式给出：

计算区域为[0, 2π]²，采用周期性边界条件，绝热指数γ=5/3.图 4给出了密度分别在T=0.5, 1, 2, 3, 4时的等高线图.观察发现，随着时间推移，解在计算区域形成了几个激波和一个漩涡结构，许多细微的复杂特征都被新格式捕捉到并且没有产生伪振荡.由2.2小节可以得知，熵稳定格式的总熵由于熵耗散的原因从而不是恒定的，而是随着时间递减.图 5(a)给出了总熵随时间变化曲线，可以发现总熵在开始的时间里是恒定的，这是因为间断尚未产生.在大约T=1之后，间断开始产生，总熵也随之开始递减，因此完全离散的数值格式满足熵稳定的性质并且比TeCNO格式的耗散小.

图  4  二维Orszag-Tang漩涡问题

Figure  4.  The 2D Orszag-Tang vortex problem

下载: 全尺寸图片 幻灯片

图  5  总熵随时间变化图

Figure  5.  The changes of the total entropy with time

下载: 全尺寸图片 幻灯片

注   为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.

3.5   第一转子问题

该算例由Balsara和Spicer^[21]提出，用来检验ES-WENO格式扭转Alfvén波的传播.计算区域为[0, 1]²，采用周期性边界条件.初始条件由下文给出，首先定义

其中$r_{0}=0.1, r_{1}=0.115, u_{0}=2, r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$.

若r＜r₀，则

若r₀＜r＜r₁, 则

若r>r₁, 则

其余初始条件均取常数，

绝热指数γ=1.4.计算该算例到T=0.15并在图 6中给出了密度ρ，压力p，Mach数Ma=‖u‖/2以及磁场压力‖B‖²/2的等高线图，其中$a=\sqrt{\gamma p / \rho}$.观察到得到的数值解具有高分辨率的特点，并且没有产生振荡.总熵随时间的变化曲线在图 5(b)中给出，同理论分析一致，熵随时间的推进而递减，证明新格式满足熵稳定的性质并且比TeCNO格式的耗散小.

图  6  第一转子问题

Figure  6.  The 1st rotor problem

下载: 全尺寸图片 幻灯片

3.6   二维云-激波相互作用

该问题模拟了一个由强激波向静态的高密度云传播的过程，并在激波穿过云时产生了非常复杂的结构.该算例用来验证ES-WENO格式高分辨率的特性.计算区域为[0, 1]²，绝热指数γ=5/3.首先考虑在x=0.05处设置激波，考虑初始条件如下：

之后以(0.25, 0.5)为圆心，半径为0.15，设置密度为1的静态云.左边界采用流入边界条件，其余边界均为流出边界条件.图 7给出了该算例在T=0.06时的密度对数ln ρ以及磁场范数‖B‖.从图中可以发现，激波与密度云的相互作用在前方产生了弓形激波，在后方产生了尾形激波，并且在内部产生了类似于湍流的复杂结构.

图  7  云-激波相互作用

Figure  7.  Cloud-shock interactions

下载: 全尺寸图片 幻灯片

4.   结论

针对理想MHD方程组，在熵稳定格式的基础上，通过对耗散项进行满足符合性质的高阶WENO重构，并且使用任意偶数阶的熵守恒通量以及对Godunov源项进行类似的高阶离散，得到了高阶熵稳定格式.新格式不仅满足熵稳定性质，并且有效抑制了熵稳定格式的抹平现象，具有更低的耗散.通过对一维、二维MHD方程标准算例进行测试，该格式具有低耗散、高精度、无振荡并且可以有效捕捉间断等特点.

附录A.   理想MHD方程特征向量以及特征值的离散形式

方便起见，定义如下符号:

其中$\boldsymbol{b}=\frac{\boldsymbol{B}}{\sqrt{\rho}}$.定义n^⊥为由b和n张成并且与法向量n正交的单位向量，本文选择将n^⊥表示为如下形式:

原始变量p_r=[ρ  u^T  p  B^T]^T的特征向量以及对应特征值的紧致形式由下式给出.为了得到守恒变量U=[ρ  ρu^T  ρe  B^T]^T的特征向量矩阵，需要对原始变量特征向量矩阵左乘Jacobi矩阵：

其中0为3×1的零向量，O和I分别为3×3的零矩阵和单位矩阵.

熵波和散度波的特征值为λ_{1, 2}=u·n，特征向量为

Alfvén波的特征值为λ_±a=u·n±b·n，特征向量为

快磁声波的特征值为λ_±f=u·n±c_f，特征向量为

慢磁声波的特征值为λ_±s=u·n±c_s，特征向量为

其中特征向量矩阵和特征值变量取值均考虑为以下平均状态：

附录B.   ES-WENO格式的熵稳定性证明

本节只考虑一维情况下ES-WENO格式的熵稳定证明，二维情况可直接推广.为了方便证明，引入熵势$\psi_{1, i}=\boldsymbol{V}_{i}^{\mathrm{T}} \boldsymbol{F}\left(\boldsymbol{U}_{i}\right)-$ $q_{1}\left(\boldsymbol{U}_{i}\right)+\phi\left(\boldsymbol{V}_{i}\right) B_{1, i}, \phi^{\prime}\left(\boldsymbol{V}_{i}\right)$与V_i的关系式为$\phi^{\prime}\left(\boldsymbol{V}_{i}\right) \cdot \boldsymbol{V}_{i}=\phi\left(\boldsymbol{V}_{i}\right)=\phi_{i}$, 以及符号$\Delta(\cdot)_{x+1 / 2}=(\cdot)_{i+1}-(\cdot)_{i}$.

首先考虑数值格式

的熵守恒性质，上式与V_i做点积可得

其中，α_{2, 1}=4/3, α_{2, 2}=－1/6，并且

将Chandrashekar等^[15]提出的MHD方程熵守恒通量满足的充分条件

代入上式，化简可得$H_{i, i+r}=\hat{q}^{\mathrm{EC}}\left(\boldsymbol{U}_{i}, \boldsymbol{U}_{i+r}\right)+\psi_{1, i}$, 其中

同理可得$H_{i, i+r}=\hat{q}^{\mathrm{EC}}\left(\boldsymbol{U}_{i}, \boldsymbol{U}_{i+r}\right)+\psi_{1, i}$.

因此离散熵不等式

成立，其中$\hat{q}_{x+1 / 2}^{\mathrm{EC}-4}=-\frac{1}{\Delta x} \sum_{r=1}^{2} \alpha_{2, r} \sum_{s=0}^{1} \hat{q}^{\mathrm{EC}}\left(\boldsymbol{U}_{i-s}, \boldsymbol{U}_{i-s+r}\right)$.

下面证明特征熵变量重构后的数值格式的熵稳定性，此时数值格式为

上式与V_i做点乘可得

其中$\Delta \boldsymbol{v}_{x+1 / 2}^{\mathrm{T}}\left|\boldsymbol{\varLambda}_{x+1 / 2}\right|\langle\boldsymbol{v}\rangle_{x+1 / 2}=\sum_{k=1}^{8}\left(v_{i+1, k}-v_{i, k}\right)\left|\lambda_{i+1 / 2, k}\right| S_{k}\left(v_{i+1 / 2, k}^{+}-v_{i+1 / 2, k}^{-}\right) \geqslant 0$, 因此离散熵不等式

成立，因此格式满足熵稳定性质.其中$\hat{q}_{x+1 / 2}^{\mathrm{ES}-4}=\hat{q}_{x+1 / 2}^{\mathrm{EC}-4}-\frac{1}{2} \overline{\boldsymbol{V}}_{x+1 / 2}^{\mathrm{T}} \boldsymbol{R}_{x+1 / 2}\left|\boldsymbol{\varLambda}_{x+1 / 2}\right|\langle\boldsymbol{v}\rangle_{x+1 / 2}$为对应的数值熵通量函数.