相对论流体力学方程的MUSCL-型熵相容格式

The MUSCL-Type Entropy Consistent Schemes for Relativistic Hydrodynamics Equations

摘要: 本文基于MUSCL-Hancock数据重构方法提出了一种求解相对论流体力学方程的高分辨率熵相容格式(EC-MHM格式)。首先将一个修正的斜率限制器应用到MUSCL-Hancock方法的数据重构中，使之与熵相容格式结合，从而得到高分辨率的熵相容通量；在时间上，主要利用双曲守恒律方程的守恒型差分形式来更新下一时间层，从而提高了格式的计算效率。文中还证明了熵相容格式的收敛性。新格式在解的光滑区域具有高精度的特性，然而在间断区域，EC-MHM格式可以有效抑制非物理现象的发生；最后通过一系列数值算例验证了新格式具有无振荡、高分辨率等良好性能。

Abstract: This paper presents a high-resolution entropy-consistent scheme, termed the EC-MHM (Entropy-consistent, MUSCL-type High-resolution Method), for solving relativistic hydrodynamics equations, based on the MUSCL-Hancock data reconstruction methodology. Firstly, a modified slope limiter is applied to the data reconstruction of MUSCL-Hancock method, and combine it with the entropy consistent scheme, so as to obtain the high resolution entropy consistent flux. For the discretization of time derivative, the conservative finite difference scheme of hyperbolic conservation laws is adopted to update the solution at the next time level. The convergence of the entropy consistent scheme is also proved. In regions where the solution is smooth, the EC-MHM scheme exhibits high precision characteristics. Conversely, in discontinuous zones, the EC-MHM can effectively prevent the occurrence of non-physical phenomena. Finally, a series of numerical examples are simulated, and the new scheme is proved to have good properties such as no oscillation and high resolution.

1. 引言

相对论流体力学(Relativistic Hydrodynamics, RHD)是计算流体力学的一个分支，在天体物理学、宇宙学等多种学科中占有重要地位。它是一门牵涉多个理论与应用方面的综合性学科，与现代物理学、天体物理学、宇宙学等多种学科有着密切的联系。与计算流体力学一样，RHD数值模拟是采用数值计算的方法来求解相对论流体力学的控制方程，得到流体的运动情况。此方式能够打破采用观测方法的局限性，以及简化复杂的理论分析；由于RHD方程组可以向双曲守恒律方程组转化，于是可以将求解两类方程组的数值方法相互联系，一些经典的数值格式也成功应用于RHD方程中。

在对双曲守恒律方程组的求解中，即使初值条件光滑，数值解也可能出现间断现象。在实际问题中，这些现象表现为熵增；那么为了解决间断带来的熵增问题，1954年，Lax提出了弱解[1]的概念，但弱解不唯一，还需添加其他条件得到唯一的弱解，于是Lax又提出了与“粘性消失”机制等价的熵稳定条件[2]，为后续研究物理解提供了一种直接有效的途径。Tadmor在1987年提出了一类熵守恒格式[3]和满足熵稳定格式的条件：一个守恒型的三点差分格式，只要含有比熵守恒格式更多的数值粘性，则该格式是熵稳定的。这为构造熵稳定格式提供了理论保证。2006年，Roe将数值粘性项添加到熵守恒格式上，得到了一类熵稳定格式[4]，称之为ERoe格式。2009年，Roe、Ismail两个人在对熵增现象进行更深入的分析后，得到了更加精确量化的熵相容格式[5]，并且该格式的数值结果更为精确。2019年，Tang等人将“熵”的相关理论推广应用到相对论流体力学方程中[6]，给出其熵对及熵守恒格式、高阶熵稳定格式的数值通量，虽然该格式在处理复杂边界问题时具有高精度特性，但是计算量较大。

Van Leer提出的MUSCL数据重构方法[7]也称为变量外推法，是一种高阶方法，近年来基于MUSCL数据重构得到了许多高阶数值格式，如MUSCL-Hancock方法，广义黎曼问题法和斜率限制中心格式等[7]。随着近四十年的深入研究，在求解双曲守恒律方程时，相对论流体力学方程数值求解方法得到了快速的发展。在国内的汤华中教授也构造了一系列有效的格式[6] [8]-[10]。2021年，任璇、封建湖得到了一种新方法，即基于MUSCL数据重构提出的斜率限制器方法[11]，提高了格式的精度和分辨率；之后高凡琪等[12]将该方法应用到相对论流体力学方程中，并通过不同数值算例验证了该格式的一些良好特性。本文将在之前的基础上进一步研究RHD方程的熵相容格式，将基于MUSCL方法的斜率限制器进行修正，并将它应用到MUSCL-Hancock方法的数据重构中，使之与熵相容格式结合，构造一种在时间和空间上全离散的二阶高分辨率熵相容格式，最后通过RHD方程的一系列数值算例证明了新格式具有无振荡、高分辨率等良好性能。

2. 一维相对论流体力学方程

一维相对论流体力学(RHD)方程的守恒律形式为

$U_t+F{\left(U\right)}_x=0$ (1)

其中$U={\left[D,m,E\right]}^T$ 为守恒型变量，$F\left(U\right)={\left[Du,mu+p,m\right]}^{\text{T}}$ 为守恒型通量，质量密度，动量密度，能量密度分别为：

$D=\rho W,m=DhWu,E=DhW-p.$ (2)

其中$\rho$ 表示物质的固有密度，p为压力，u是归一化(取光速$c=1$ )的流体速度；$W=\frac{1}{\sqrt{1-u^2}}$ 为Lorentz因

子，$h=1+e+p/\rho$ 表示比焓，e是比内能。对于理想气体，可添加状态方程

$p=\left(\gamma -1\right)\rho e$ (3)

使方程(1)封闭，其中$\gamma \in \left(1,2\right]$ 是绝热指数。

方程(1)的雅可比矩阵为[13]

$A=\frac{\partial F\left(U\right)}{\partial U}=\frac{1}{N}\left(\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ 0& A_{32}& 0\end{array}\right)$ (4)

其中$A_{11}=\frac{u^{2}h\left(1-c_s^2\right)}{\gamma -1}+\frac{uh}{W^2}$ ，$A_{12}=-\frac{1}{W^3}+\frac{\gamma }{W\left(\gamma -1\right)}$ ，$A_{13}=-\frac{u\gamma }{W\left(\gamma -1\right)}$ ，$A_{21}=-\frac{h^2}{W^3}\left(1-\frac{c_s^2}{\gamma -1}\right)$ ，$A_{22}=-\frac{uh}{W^3}\left(1-\frac{c_s^2}{\gamma -1}\right)+\frac{2uh\left(1-c_s^2\right)}{\gamma -1}$ ，$A_{23}=-\frac{u^{2}h\left(1-c_s^2\right)}{\gamma }+\frac{h}{W^2}$ ，$A_{32}=\frac{h\left(1-c_s^2{u}^2\right)}{\gamma -1}$ ，$N=\frac{h\left(1-c_s^2{u}^2\right)}{\gamma -1}$ 。

雅可比矩阵的特征值分别为${\lambda }_{-}=\frac{u-c_s}{1-u{c}_s}$ ，${\lambda }_0=u$ ，${\lambda }_{+}=\frac{u+c_s}{1+u{c}_s}$ ，$c_s=\sqrt{\left(\gamma p\right)/\left(\rho h\right)}$ 表示声速，h是比焓。

雅可比矩阵的特征值所对应的右特征向量矩阵为[10]

$R=\left(\begin{array}{ccc}1& 1& 1\\ \left(u-c_s\right)Wh& uW& \left(u+c_s\right)Wh\\ \left(1-u{c}_s\right)Wh& W& \left(1+u{c}_s\right)Wh\end{array}\right)\left(\begin{array}{ccc}\sqrt{\frac{\rho W\left(1-u{c}_s\right)}{2\gamma }}& 0& 0\\ 0& \sqrt{\frac{\left(\gamma -1\right)\rho W}{\gamma }}& 0\\ 0& 0& \sqrt{\frac{\rho W\left(1+u{c}_s\right)}{2\gamma }}\end{array}\right)$ (5)

已知雅可比矩阵$A$ 的特征值都是实数，并且特征值对应的特征向量线性无关，可知一维相对论流体力学方程是双曲型方程组。因此双曲守恒律方程组的数值求解方法可推广到RHD方程组中。但由于RHD方程有Lorentz因子，使得求解该方程变得较为复杂。联系表达式式(2)、(3)、和$h=1+e+p/\rho$ ，有

$\{\begin{array}{l}W={\left(1-m^2/{\left(E+p\right)}^2\right)}^{-\frac{1}{2}}\\ E+p=DhW=W\left(D+\frac{\gamma }{\gamma -1}pW\right)\end{array}$ (6)

为了得到原始变量$\left(\rho ,u,p\right)$ 的值，需对式(6)进行求解。此时式(6)变成关于变量p的非线性方程，在已知守恒型变量$U=\left(D,m,E\right)$ 的数值时，利用牛顿迭代法求解压力p，然后由(6)的第一式得到W，再根据

关系式$\rho =\frac{D}{W}$ ，$h=1+\frac{\gamma p}{\left(\gamma -1\right)\rho }$ ，$u=\frac{m}{DhW}$ ，求出密度$\rho$ 和速度u。

另外，为构造一维RHD方程的熵相容格式，定义其熵函数为$M\left(U\right)=-\frac{\rho Ws}{\gamma -1}$ ，熵通量函数为$N\left(U\right)=-\frac{\rho uWs}{\gamma -1}$ ，其中热力学熵$s=\ln \left(p\right)-\gamma \ln \left(\rho \right)$ ，对应的熵变量为：$V={\left[\frac{\gamma -s}{\gamma -1}+\frac{\rho }{p},\frac{\rho Wu}{p},-\frac{\rho W}{p}\right]}^{\text{T}}$ 。

3. 熵相容格式

3.1. 熵守恒格式

称在离散情况下满足以下熵等式的数值格式为熵守恒格式：

$\frac{\text{d}}{\text{d}x}M\left(U_i\left(t\right)\right)+\frac{1}{\Delta x}\left(N_{i+1/2}-N_{i-1/2}\right)=0$ (7)

其中M为凸的熵函数，$N_{i+1/2}$ 为数值熵通量函数[2]，且$N_{i+1/2}$ 与熵通量函数N相容。2019年，Tang在文献[9]中提出了求解RHD方程的熵守恒格式。通过引入参数向量$z={\left(z_1,z_2,z_3\right)}^{\text{T}}={\left(\rho ,\rho /p,u\right)}^{\text{T}}$ ，得到求解一维相对论流体力学方程的离散熵守恒通量为$F_{i+1/2}^C={\left(F_1,F_2,F_3\right)}^{\text{T}}$ ，其中

$\{\begin{array}{l}{F}_1=z_1^{\ln }\overline{\left(uW\right)},\\ F_2=q^{-1}\left[\left(1+\frac{1}{\left(\gamma -1\right)z_2^{\ln }}\right)\overline{z_2}{z}_3^{\text{Lor}}{F}_1+\overline{z_1}{\overline{W}}^2+\overline{z_1}\overline{z_3}\overline{W}{z}_3^{\text{Lor}}\right],\\ F_3=q^{-1}\left[\overline{z_1}\overline{W}\overline{\left(uW\right)}+\overline{z_1}\overline{z_3}\overline{\left(uW\right)}{z}_3^{\text{Lor}}+\left(1+\frac{1}{\left(\gamma -1\right)z_2^{\ln }}\right)F_1\left(\overline{z_2}\overline{W}+\overline{z_2}\overline{z_3}{z}_3^{\text{Lor}}\right)\right].\end{array}$ (8)

公式中$q=\overline{z_2}{\overline{W}}^2+\overline{z_2}\overline{z_3}\overline{W}{z}_3^{\text{Lor}}-\overline{z_2}\overline{\left(uW\right)}{z}_3^{\text{Lor}}$；$\overline{z_{i+1/2}}=\frac{z_i+z_{i+1}}{2}$ ；$\overline{\left(uW\right)}$ 表示u和W的乘积取平均值；对数平均[5] $z_{i+1/2}^{\ln }=\frac{z_i-z_{i+1}}{\ln \left(z_i\right)-\ln \left(z_{i+1}\right)}$ ；$z_3^{\text{Lor}}$ 是Lorentz平均，$z_3^{\text{Lor}}=u^{\text{Lor}}=\frac{u_i+u_{i+1}}{\sqrt{1-u_i^2}\sqrt{1-u_{i+1}^2}\left(\sqrt{1-u_i^2}+\sqrt{1-u_{i+1}^2}\right)}$ 。

3.2. 熵稳定格式

在解的光滑区域熵守恒格式有着良好的性能，但是在间断处会出现非物理现象。因此为了很好地捕捉到解的各个结构，需要总熵达到一定的耗散。Tadmor在文献[3]中提出满足熵稳定格式的条件：一个守恒型三点差分格式，若比熵守恒格式含有更多的数值粘性，且满足以下离散熵不等式，即为熵稳定格式：

$\frac{\text{d}}{\text{d}x}M\left(U_i\left(t\right)\right)+\frac{1}{\Delta x}\left(N_{i+1/2}-N_{i-1/2}\right)<0$ (9)

因此在本文将Roe格式的数值粘性项[4]添加到数值通量(8)中，得到求解一维RHD方程的熵稳定格式的数值通量为：

$F^{ES}=F^C-\frac{1}{2}R\left|\Lambda \right|R^{-1}\left[U\right].$ (10)

其中$\Lambda$ 是方程(1)的雅可比矩阵$A$ 的特征值构成的对角矩阵，$R$ 是$A$ 的特征值对应的右特征向量构成的矩阵，表达式为公式(5)；$[·]={(·)}_R-{(·)}_L$ ，下标L、R分别表示用于计算矩阵中的变量的左、右状态，后续将不再提及。关于上式中$\left|\Lambda \right|$ 的取法在Roe格式和Lax-Friedrichs格式两种中任取其一即可：

$\left|\Lambda \right|=\{\begin{array}{ll}\text{diag}\left(\left|{\lambda }_1\right|,\cdots ,\left|{\lambda }_n\right|\right),\hfill & Roescheme\hfill \\ \max \left(\left|{\lambda }_1\right|,\cdots ,\left|{\lambda }_n\right|\right)I,\hfill & Lax\text{-}Friedrichsscheme\hfill \end{array}$ (11)

其中$I$ 是单位矩阵，${\lambda }_1,\cdots ,{\lambda }_n$ 是雅可比矩阵$A$ 的特征值。

1999年，Barth [14]证明了$U_V\approx R{R}^{\text{T}}$ ，故有：

$\frac{1}{2}R\left|\Lambda \right|R^{-1}\left[U\right]\approx \frac{1}{2}R\left|\Lambda \right|R^{-1}{U}_V\left[V\right]=\frac{1}{2}R\left|\Lambda \right|R^{-1}R{R}^{\text{T}}\left[V\right]=\frac{1}{2}R\left|\Lambda \right|R^{\text{T}}\left[V\right]$ (12)

因此得到在离散情况下，一维相对论流体力学方程的熵稳定格式(ES)的数值通量为：

$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}$ (13)

上述式子中，$F_{i+1/2}^C$ 的表达式为公式(8)，${\Lambda }_{i+1/2}=\text{diag}\left(\frac{\overline{u_{i+1/2}}-{\left({\widehat{c}}_s\right)}_{i+1/2}}{1-\overline{u_{i+1/2}}{\left({\widehat{c}}_s\right)}_{i+1/2}},\overline{u_{i+1/2}},\frac{\overline{u_{i+1/2}}+{\left({\widehat{c}}_s\right)}_{i+1/2}}{1+\overline{u_{i+1/2}}{\left({\widehat{c}}_s\right)}_{i+1/2}}\right)$，其中$\overline{u_{i+1/2}}=\frac{u_i+u_{i+1}}{2}$ ，$c_s=\sqrt{\frac{\gamma p}{\rho h}}$ 为音速，上式中${\left({\widehat{c}}_s\right)}_{i+1/2}=\sqrt{\frac{\gamma {\widehat{p}}_{i+1/2}}{{\widehat{\rho }}_{i+1/2}{\widehat{h}}_{i+1/2}}}$ ，${\widehat{\rho }}_{i+1/2}={\rho }_{i+1/2}^{\ln }$ ，${\widehat{p}}_{i+1/2}=\frac{{\rho }_{i+1/2}^{\ln }}{{\left(p/\rho \right)}_{i+1/2}^{\ln }}$ ，${\widehat{h}}_{i+1/2}=1+\frac{\gamma }{\gamma -1}\frac{{\widehat{p}}_{i+1/2}}{{\widehat{\rho }}_{i+1/2}}$ ，${\left[V\right]}_{i+1/2}=V_{i+1}-V_i$ 。

3.3. 熵相容格式

尽管熵稳定(ES)格式比熵守恒格式含有更多的数值粘性，但当跨过激波时产生的熵增不够大时，熵稳定格式的耗散机制不能完全抵消间断区域处的熵增量，因此仍然不能避免出现伪振荡现象。只有继续增加耗散量对间断区域处的熵增进行充分地耗散，才有可能得到对数值粘性项更精确量化的熵稳定格式。Ismail和Roe详细研究了Euler方程组的解穿过激波时的熵增[5]，提出了一种对熵增的估计更为精确地熵稳定格式，并将其称之为熵相容格式。因此在本文中将文献[15]中的熵相容格式的表达式推广到相对论流体力学方程中，其数值通量为：

$F^{EC}=F^C-\frac{1}{2}R\left(\left|\Lambda \right|+{\left(\frac{1}{6}\left|\left[\Lambda \right]\right|\right)}^{+}\right)R^{-1}\left[U\right]$ (14)

利用关系式(12)，将式中的数值粘性项写成熵变量的差分形式，可得到熵相容格式数值通量的另一种表示，从而得到一维相对论流体力学方程的熵相容格式的数值通量为：

$F^{EC}=F^C-\frac{1}{2}R\left(\left|\Lambda \right|+{\left(\frac{1}{6}\left|\left[\Lambda \right]\right|\right)}^{+}\right)R^{\text{T}}\left[V\right]$ (15)

其中，${\left(\frac{1}{6}\left|\left[\Lambda \right]\right|\right)}^{+}$ 表示矩阵中的负数替换为0，其余数值不变的矩阵[16]。

因此得到在离散情况下，一维RHD方程的熵相容格式(EC)的数值通量为：

$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(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}\right)R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}$ (16)

引理[3]：若与$F$ 相容的数值通量函数$F_{i+1/2}=F\left(U_i,U_{i+1}\right)$ 关于$U_i$ ，$U_{i+1}$ 是Lipschitz连续的，那么存在一个数值粘性矩阵$G_{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}{G}_{i+1/2}\left(U_{i+1}-U_i\right)$ (17)

定理1：数值通量为(16)的熵相容格式(EC)是收敛的。

证明：前文提到的熵格式均是守恒型三点格式，将EC格式的数值通量$F_{i+1/2}^{EC}$ 写成带数值粘性的形式：

$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{G}}_{i+1/2}^{EC}{\left[U\right]}_{i+1/2}$

其中${\tilde{G}}_{i+1/2}^{EC}=G_{\text{i}+1/2}^C+R_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+{\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}\right)R_{i+1/2}^{-1}$，$G_{i+1/2}^C$ 是熵守恒格式对应的数值粘性矩阵。

从关系式(11)可知，熵相容格式的数值通量$F_{i+1/2}^{EC}$ 可写成带数值粘性的形式：

$F_{i+1/2}^{EC}=\frac{1}{2}\left(F\left(U_i\right)+F\left(U_{i+1}\right)\right)-\frac{1}{2}{G}_{i+1/2}^{EC}{\left[V\right]}_{i+1/2}$

其中$G_{i+1/2}^{EC}=G_{i+1/2}^C+R_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+{\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}\right)R_{i+1/2}^{\text{T}}$ ，

令$P_{i+1/2}^{EC}=G_{i+1/2}^{EC}-G_{i+1/2}^C=R_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+{\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}\right)R_{i+1/2}^{\text{T}}$ ，则

$\begin{array}{c}{\left[V\right]}_{i+1/2}^{\text{T}}{P}_{i+1/2}^{EC}{\left[V\right]}_{i+1/2}={\left[V\right]}_{i+1/2}^{\text{T}}{G}_{i+1/2}^{EC}{\left[V\right]}_{i+1/2}-{\left[V\right]}_{i+1/2}^{\text{T}}{G}_{i+1/2}^C{\left[V\right]}_{i+1/2}\\ ={\left[V\right]}_{i+1/2}^{\text{T}}{R}_{i+1/2}\left(\left|{\Lambda }_{i+1/2}\right|+{\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}\right)R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}\\ ={\left(R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}\right)}^{\text{T}}\left|{\Lambda }_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}+{\left[V\right]}_{i+1/2}^{\text{T}}{R}_{i+1/2}{\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}{R}_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}\end{array}$

其中${\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}$ 是半正定矩阵，$\left|{\Lambda }_{i+1/2}\right|$ 是正定矩阵，$R_{i+1/2}^{\text{T}}$ 是可逆矩阵，${\left[V\right]}_{i+1/2}$ 为非零向量，从而$R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}$ 是非零向量，有：

${\left(R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}\right)}^{\text{T}}\left|{\Lambda }_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}>0$ ，${\left[V\right]}_{i+1/2}^{\text{T}}{R}_{i+1/2}{\left(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}{R}_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}\ge 0$ ，

故${\left[V\right]}_{i+1/2}^T{P}_{i+1/2}^{EC}{\left[V\right]}_{i+1/2}>0$ ，即与熵守恒格式相比EC格式含有更多的数值粘性。由比较原则[3]可得，EC格式满足熵稳定条件，因此数值解收敛。

4. 高分辨率熵相容格式

在每个控制单元$I_i=\left[x_{i-1/2},x_{i+1/2}\right]$ 上，将双曲守恒律方程的精确解$u\left(x\right)$ 表示为分段的精确解$u_i\left(x\right)$ ，然后对$u_i\left(x\right)$ 进行重构，记重构之后的函数为${\tilde{u}}_i\left(x\right)$ 。如果有不等式$\left|u_i\left(x\right)-{\tilde{u}}_i\left(x\right)\right|\le O\left({\left(\Delta x\right)}^{k+1}\right)$ 成立，则称重构函数${\tilde{u}}_i\left(x\right)$ 具有k阶精度。

$u_i^n$ 表示在时间层$t^n$ 上单元$I_i$ 的平均值，表达式为：$u_i^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\tilde{u}}_i\left(x\right)\text{d}x}=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{u}_i\left(x\right)\text{d}x}$，之后对一个控制体单元上的$u_i^n$ 进行分段线性重构，可将其表示为：

${\tilde{u}}_i\left(x\right)=u_i^n+\left(x-x_i\right)\frac{{\Delta }_i}{{\Delta }_x},x\in \left[x_{i-1/2},x_{i+1/2}\right]$ (18)

上式中，$\frac{{\Delta }_i}{{\Delta }_x}$ 是单元$I_i$ 上选择的合适斜率。${\tilde{u}}_i\left(x\right)$ 在单元$I_i$ 的极值分别表示为：

$u_i^L={\tilde{u}}_i\left(x_{i-1/2}\right)=u_i^n-\frac{1}{2}{\Delta }_i$ ，$u_i^R={\tilde{u}}_i\left(x_{i+1/2}\right)=u_i^n+\frac{1}{2}{\Delta }_i$ (19)

通常地，我们称$u_i^L$ 与$u_i^R$ 为边界外推值，称整个重构方法为边界外推法[7]。另外，$u_i^n$ 和重构函数${\tilde{u}}_i\left(x\right)$ 在控制体单元$I_i$ 上的积分值相等，因此确保重构过程中的守恒性。

4.1. 斜率限制器

首先将${\tilde{u}}_i\left(x\right)$ 在$x_i$ 点处进行泰勒展开，得：

${\tilde{u}}_i\left(x\right)=u_i\left(x_i\right)+{\frac{\partial u_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}\frac{{\left(x-x_i\right)}^2}{2}+\cdots$ (20)

同样地将$u_i\left(x\right)$ 在点$x_i$ 处进行泰勒级数展开，得：

$\begin{array}{l}{u}_i^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{u}_i\left(x\right)\text{d}x}=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}\left(u_i\left(x_i\right)+{\frac{\partial u_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}\frac{{\left(x-x_i\right)}^2}{2}+\cdots \right)\text{d}x}\\ =u_i\left(x_i\right)+0+\frac{1}{\Delta x}{\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}{\displaystyle {\int }_{I_i}\frac{{\left(x-x_i\right)}^2}{2}\text{d}x}+\cdots \end{array}$ (21)

联合得到：${\tilde{u}}_i\left(x\right)=u_i^n+{\frac{\partial u_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{1}{2}\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}\left[{\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\left(x-x_i\right)}^2\text{d}x}\right]+\cdots$

为了方便后续计算，将${\tilde{u}}_i\left(x\right)$ 的表达式记为：

${\tilde{u}}_i\left(x\right)=u_i^n+{\frac{\partial u_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{1}{2}\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}\left[{\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\left(x-x_i\right)}^2\text{d}x}\right]$ (22)

满足不等式$\left|{\tilde{u}}_i\left(x\right)-u_i^n\right|\le O\left({\left(\Delta x\right)}^3\right)$ ，即重构后的格式具有二阶精度。

可以将其写成：

${\tilde{u}}_i\left(x\right)=u_i^n+H\left(x-x_i\right)$ (23)

其中，$H\left(x-x_i\right)={\frac{\partial u_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{1}{2}\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}\left[{\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\left(x-x_i\right)}^2\text{d}x}\right]$ 。

本文采用中心差商对一阶、二阶偏导数进行处理，即：

${\frac{\partial u_i}{\partial x}|}_{x_i}=\frac{u_{i+1}^n-u_{i-1}^n}{2\Delta x}$ ，${\frac{{\partial }^2{u}_i}{\partial x^2}|}_{x_i}=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}$ (24)

对重构式(23)进行限制，有：

${\tilde{u}}_i\left(x\right)=u_i^n+{\phi }_{i}H\left(x-x_i\right)$ (25)

其中${\phi }_i\in \left[\text{0,1}\right]$ 是一个斜率限制器，将其加到(23)式中，是为了使得解在光滑区域可以达到高精度的特性，而在间断区域能够选择原始数据，从而避免产生伪振荡。${\tilde{u}}_i\left(x\right)$ 在每个控制单元$I_i=\left[x_{i-1/2},x_{i+1/2}\right]$ 的左右边界外推值分别为：

$u_i^L=u_i^n+{\phi }_{i}H\left(x_{i-1/2}-x_i\right)$ ，$u_i^R=u_i^n+{\phi }_{i}H\left(x_{i+1/2}-x_i\right)$ (26)

为确保不会产生新的局部极值，这二者需满足下列不等式：

$\{\begin{array}{c}\min \left(u_{i-1}^n,u_i^n\right)\le u_i^L\le \max \left(u_{i-1}^n,u_i^n\right)\\ \min \left(u_{i+1}^n,u_i^n\right)\le u_i^R\le \max \left(u_{i+1}^n,u_i^n\right)\end{array}$ (27)

上式中不等式两边同时减去$u_i^n$ 得到：

$\{\begin{array}{c}\min \left(u_{i-1}^n-u_i^n,0\right)\le u_i^L-u_i^n\le \max \left(u_{i-1}^n-u_i^n,0\right)\\ \min \left(u_{i+1}^n-u_i^n,0\right)\le u_i^R-u_i^n\le \max \left(u_{i+1}^n-u_i^n,0\right)\end{array}$ (28)

为保证上述不等式成立，文献[12]给出了修正的斜率限制器为：

$\{\begin{array}{c}{\phi }_L=\min \left(1,\frac{u_{i-1}^n-u_i^n}{H\left(x_L-x_i\right)}\right),\\ {\phi }_R=\min \left(1,\frac{u_{i+1}^n-u_i^n}{H\left(x_R-x_i\right)}\right),\\ {\phi }_i=\max \left(\min \left({\phi }_L,{\phi }_R\right),0\right).\end{array}$ (29)

4.2. MUSCL-Hancock重构方法

已知一维情况下，RHD方程的熵相容格式数值通量为

$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(\frac{1}{6}\left|{\left[\Lambda \right]}_{i+1/2}\right|\right)}^{+}\right)R_{i+1/2}^{\text{T}}{\left[V\right]}_{i+1/2}$ (30)

近年来基于MUSCL数据重构得到了许多高阶数值格式，如MUSCL-Hancock重构方法[7]，它在提高了数值精度和分辨率的同时，也能够保持本质无振荡。在本文中应用该重构方法在每个单元交界面处只需要求解一次数值通量，从而减小了计算量。

MUSCL-Hancock方法的计算步骤如下：

第一步：与MUSCL数据重构方法[7]类似，利用式(26)得到$t_n$ 层重构$U_i\left(x\right)$ 在单元交界面处的左右极限值，也叫做边界外推值$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]$ , (31)

${\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].$ (32)

第三步：将上述的边界外推值${\overline{U}}_i^L$ ，${\overline{U}}_i^R$ 代入到一维RHD方程的熵相容通量中，从而得到高分辨率熵相容格式(EC-MHM)的数值通量为：

$F_{i+1/2}^{EC\text{-}MHM}={\tilde{F}}_{i+1/2}^C-\frac{1}{2}{\tilde{R}}_{i+1/2}\left(\left|{\tilde{\Lambda }}_{i+1/2}\right|+{\left(\frac{1}{6}\left|{\left[\tilde{\Lambda }\right]}_{i+1/2}\right|\right)}^{+}\right){\tilde{R}}_{i+1/2}^{\text{T}}{\left[\tilde{V}\right]}_{i+1/2}$ (33)

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

第四步：利用以下守恒型差分格式更新$t_{n+1}$ 时刻的单元平均值：

$U_j^{n+1}=U_j^n-\frac{\Delta t}{\Delta x}\left(F_{i+1/2}^{EC\text{-}MHM}-F_{i-1/2}^{EC\text{-}MHM}\right)$ (34)

在实际重构中，首先用牛顿迭代法求解非线性方程(公式(6))中的p，进而利用公式(2)的关系式，将守恒型变量$U$ 变为原始变量，即：

$\left[\rho W,\rho h{W}^{2}u,\rho hW-p\right]\to \left[\rho ,u,p\right]$ (35)

对其进行MUSCL数据重构后再回到守恒型变量，其次利用上述中的第二步得到时间推进后的边界外推值，代入熵相容格式中得到近似的数值通量$F_{i+1/2}^{EC\text{-}MHM}$ 。我们把格式(33)称为MUSCL-Hancock型高分辨率熵相容格式。

5. 数值算例

以下Riemann问题中，控制方程均为式(1) (2)，边界条件均为Neumann边界条件，绝热系数$\gamma =5/3$ ，算例中的熵稳定格式(ES)、熵相容格式(EC)与高分辨率熵相容格式(EC-MHM)都是由400个网格点得到的，数值通量分别为公式(13)、公式(16)、公式(33)；算例中的每个参考解是由文献[11]中的高分辨率熵相容格式采用4000个网格点计算得到。

算例1 在区域$\left[0,1\right]$ 上考虑一维相对论流体力学方程的如下初始问题：

$\rho \left(x,0\right)=1+0.2\sin \left(2\pi x\right)$ ，$u\left(x,0\right)=0.2$ ，$p\left(x,0\right)=1$

在计算时采用周期性边界条件，计算到$t=2$ 时刻。该算例的精确解为：

$\rho \left(x,t\right)=1+0.2\sin \left(2\pi \left(x-ut\right)\right)$ ，$u\left(x,t\right)=0.2$ ，$p\left(x,t\right)=1$

表1给出了当网格数N分别取不同值时密度$\rho$ 的数值误差以及收敛阶，同时图1展示了相应的变量$\rho$ 的精确解与数值解的结果对比图。从表1给出的结果来看，新构造的高分辨率熵相容格式在$L^1$ 和$L^2$ 范数意义下的收敛阶均超过二阶，在$L^{\infty }$ 范数意义下收敛阶接近于二阶，因此该格式(EC-MHM格式)是具有高精度的特性的。

算例2 在区域$\left[0,1\right]$ 上给定如下初始条件[17]：

$\left(\rho ,u,p\right)=\{\begin{array}{l}\left(5,0,50\right),x<0.5,\\ \left(2+0.3\sin \left(50x\right),0,5\right),x>\mathrm{0.5.}\end{array}$

图2表示$t=0.35$ 时刻各格式的结果对比图，从图中可以看出EC和ES格式对间断处的捕捉效果不好，不能准确捕捉到波的峰值；然而EC-MHM格式不仅可以准确识别出激波和稀疏波的位置，而且该格式的总熵耗散相对较小；总之结果表明三个格式中EC-MHM格式既能有效抑制抹平现象，又能很好地捕捉到解的结构。

Table 1. The numerical error and convergence order of density $\rho$

表1. 密度$\rho$ 的数值误差和收敛阶

网格数	$L^1$ 范数	收敛阶	$L^2$ 范数	收敛阶	$L^{\infty }$ 范数	收敛阶
25	2.4000e−3	—	3.7000e−3	—	9.8000e−3	—
50	4.6270e−4	2.3749	7. 7894e−4	2.2479	2.4000e−3	1.9709
100	9.0770e−5	2.3498	1.8115e−4	2.1043	8.3405e−4	1.5837
200	1.5863e−5	2.5165	4.1936e−5	2.1109	2.5079e−5	1.7337
400	3.3313e−6	2.2515	9.7303e−6	2.1076	7.5464e−6	1.7326

Figure 1. Comparison of results under different grids

图1. 不同网格下的结果对比

Figure 2. Numerical results for Case 2

图2. 算例2的数值结果

算例3 在区域$\left[0,1\right]$ 上给定如下初始条件：

$\left(\rho ,u,p\right)=\{\begin{array}{l}\left(1,-0.7,20\right),x<0.5,\\ \left(1,0.7,20\right),x>\mathrm{0.5.}\end{array}$

图3给出了计算到$t=0.4$ 时刻的数值解和参考解，该算例的解是由两个反向的稀疏波组成。由图可知，在$x=0.5$ 附近，出现了“undershoot”现象，可以明显看出新构造的EC-MHM格式更为逼近参考解，并且相对于ES格式和EC格式该格式(EC-MHM)对稀疏波的捕捉效果更好。

算例4 在区域$\left[0,1\right]$ 上给定如下初始条件：

$\left(\rho ,u,p\right)=\{\begin{array}{l}\left(1,0,{10}^3\right),x<0.5,\\ \left(1,0,{10}^{-2}\right),x>\mathrm{0.5.}\end{array}$

图4给出了$t=0.4$ 时刻得到的参考解和数值解，该算例均采用800个网格点计算。由图可知稀疏波较为弯曲，接触间断与激波之间的区域比较窄。对于该算例来说，ES格式和EC格式的解差别不大，计算得到的激波位置相对不够准确，从图中可以看到EC-MHM格式对激波的捕捉效果明显更好，有效抑制

Figure 3. Numerical results for Case 3

图3. 算例3的数值结果

Figure 4. Numerical results for Case 4

图4. 算例4的数值结果

了过度抹平现象；并且，EC-MHM格式对稀疏波的捕捉也较为精确，相比于ES格式和EC格式而言该格式对稀疏波的波头和波尾的捕捉效果要好。

算例5 在区域$\left[0,1\right]$ 上给定如下初始条件：

$\left(\rho ,u,p\right)=\{\begin{array}{l}\left(10,0,\frac{40}{3}\right),x<0.5,\\ \left(1,0,{10}^{-6}\right),x>\mathrm{0.5.}\end{array}$

图5给出了各格式计算到$t=0.4$ 时刻的数值解和参考解。该算例的解是由向左的稀疏波、接触间断和向右的激波组成。从图中可以看出，与EC格式相比，新构造的格式(EC-MHM)产生的总熵耗散较小，有效地改善了数值解在间断附近的抹平现象。在跨越间断时，EC-MHM格式的单元数目明显比熵相容格式(EC)少，并且较为逼近参考解，提高了分辨率。

Figure 5. Numerical results for Case 4

图5. 算例5的数值结果

6. 结论

本文基于MUSCL-Hancock方法构造了一种求解相对论流体力学方程的高分辨率熵相容格式，并证明了熵相容格式的收敛性。新构造的EC-MHM格式在解的光滑区域能够达到高精度，且该格式在间断区域可以有效抑制伪震荡的发生；同时它是全离散守恒型格式，相比半离散Runge-Kutta方法计算效率更高；最后通过RHD方程的一系列数值算例验证了EC-MHM的格式具有无振荡、高分辨率等良好性能。

NOTES

^*通讯作者。