二维磁流体方程的高分辨率旋转通量格式

doi:10.3785/j.issn.1008-9497.2024.01.005

二维磁流体方程的高分辨率旋转通量格式

郑素佩^,, 翟梦情^,^,, 李琦, 建芒芒

长安大学理学院，陕西西安 710064

High resolution rotated flux scheme for two-dimensional magnetohydrodynamics equations

ZHENG Supei^,, ZHAI Mengqing^,^,, LI Qi, JIAN Mangmang

School of Science，Chang'an University，Xi'an 710064，China

通讯作者: ORCID：https：//orcid.org/0000-0001-6184-6788，E-mail：zhaimengqing2016@163.com.

收稿日期: 2022-09-27 修回日期: 2023-02-27 接受日期: 2023-03-06

基金资助:

国家自然科学基金资助项目. 11971075. 12101073
陕西省自然科学基金青年项目. 2020JQ-338

Received: 2022-09-27 Revised: 2023-02-27 Accepted: 2023-03-06

作者简介 About authors

郑素佩（1978—），ORCID：https：//orcid.org/0000-0003-2502-6998，女，博士，副教授，主要从事微分方程数值算法研究. 。

摘要

若待解方程满足旋转不变性，则可通过旋转通量法有效消除近似Riemann求解器的激波不稳定现象，抑制非物理现象的产生。针对二维理想磁流体（MHD）方程和浅水波磁流体（SWMHD）方程，构造了通量函数的类旋转矩阵，给出了方程的旋转不变性证明；根据该性质对控制方程做类一维处理，推导了方程的半离散旋转通量格式；利用通量限制器，将熵稳定通量和反扩散通量进行加权组合，得到能够自适应调整耗散量的高分辨率旋转通量格式。数值实验表明，此格式能精确捕捉解的结构，分辨率高、鲁棒性强，且易向高维推广。

关键词： 理想磁流体方程 ; 浅水波磁流体方程 ; 旋转不变性 ; 高分辨率熵稳定通量

Abstract

The rotated flux method could be used to effectively eliminate the shock instability of the approximate Riemann solver and suppress the generation of non-physical phenomenon if the equations to be solved meet the rotational invariance. For the 2D ideal magnetohydrodynamics (MHD) and shallow water magnetohydrodynamics (SWMHD) equations, the rotation-like matrix of flux function was constructed, and the corresponding rotational invariance theorem was given with proof, which was then used to deal with the governing equations applying quasi-1D method to derive the semi-discrete rotated flux scheme. Combining entropy stable flux and anti-diffusive flux by a flux limiter, a new flux that can adaptively adjust the dissipation term was obtained. Numerical experiments show that the new scheme can accurately capture the structure of solution, has high resolution, strong robustness and can be easily extended to higher dimensions.

Keywords： ideal magnetohydrodynamics equations ; shallow water magnetohydrodynamics equations ; rotational invariance ; high resolution entropy stable flux

PDF (2770KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

郑素佩, 翟梦情, 李琦, 建芒芒. 二维磁流体方程的高分辨率旋转通量格式. 浙江大学学报(理学版)[J], 2024, 51(1): 29-40 doi:10.3785/j.issn.1008-9497.2024.01.005

ZHENG Supei, ZHAI Mengqing, LI Qi, JIAN Mangmang. High resolution rotated flux scheme for two-dimensional magnetohydrodynamics equations. Journal of Zhejiang University(Science Edition)[J], 2024, 51(1): 29-40 doi:10.3785/j.issn.1008-9497.2024.01.005

磁流体力学是将经典流体力学和电动力学相结合研究导电流体和电磁场相互作用的学科，在等离子体物理学、天体物理以及流动控制等领域有广泛应用。较为经典的数学模型有理想磁流体（MHD）方程和浅水波磁流体（SWMHD）方程等，当满足适当条件时，在垂直方向对三维理想MHD方程进行积分，即可获得SWMHD方程。鉴于这两类方程均具有类双曲型特点，考虑将双曲型方程的数值求解格式应用于此两类方程。近年来，此类推广已有一系列成果^［1-7］，其中不乏熵稳定格式的推广。熵稳定格式是从物理概念出发，较熵守恒格式含有更多黏性^［8］，其构造理念符合间断位置产生熵增的规律，因而能够有效避免伪振荡的发生，得到唯一且具有物理意义的解^［9］。2016年，WINTERS等^［4，6］以文献［10］为基础，针对此两种方程构造了一类熵稳定数值通量格式，将熵相关理论推广到MHD方程的数值求解，但该通量格式未能实现对间断的高精度捕捉。郑素佩等^［11-12］基于高阶重构陆续得到多种高分辨率熵稳定格式，将求解精度提升到五阶。2021年，DUAN等^［5］提出的高阶精确熵稳定有限差分格式，成功求解了SWMHD方程，验证了熵稳定格式具有良好的性能。沈亚玲等^［13］引入限制器机制构造了一类求解理想MHD方程的高分辨率熵相容格式，该格式能合理控制耗散，有效改善了抹平现象。

针对二维问题，即使常用的Roe近似Riemann求解器也可能出现包括红斑现象在内的非物理问题^［14］，造成出现此现象的主要原因是近似Riemann求解器存在激波不稳定性问题，该问题可通过旋转通量得以改善或解决^［14］。旋转通量法，尝试将边界通量分解为2个正交分量，利用通量函数的旋转不变性对二维方程做类一维处理，将其简化为一维形式，而前提条件是待解方程必须满足旋转不变性。DAVID等^［15］验证了旋转Riemann求解器在多个迎风角度下的求解效果与高阶MUSCL格式相当。NISHIKAWA等^［16］基于Riemann求解器、结合Roe格式和HLL格式得到一种混合格式，该混合格式能自动切换至少波（HLL）或全波（Roe）求解器，在抑制红斑现象的同时避免过度耗散。ZHANG等^［17］利用加权压强函数定义旋转角度，实现对耗散影响的自适应调整，从而有效消除了由动量扰动造成的激波不稳定。郑素佩等^［18-19］引入了HLL数值通量，构造了求解Euler方程和浅水波方程的混合旋转通量格式，该格式间断捕捉能力好，分辨率高。

鉴于此，本文重点研究求解理想MHD方程和SWMHD方程的旋转通量法，给出类旋转矩阵以及对应的旋转不变性命题，引入通量限制器和反扩散通量构造高分辨率旋转通量格式，最后用数值算例验证了所构造的高分辨率旋转通量格式的性能。

1　控制方程

二维理想MHD方程由流体力学中的Navier-Stokes方程和电动力学中的Maxwell方程耦合而成，其守恒律形式^［4］为

\frac{\partial q}{\partial t} + \frac{\partial f (q)}{\partial x} + \frac{\partial g (q)}{\partial y} = 0

，（1）

其中， $q = {[ρ, ρ u, ρ v, ρ w, ρ e, B_{1}, B_{2}, B_{3}]}^{T}$ ， $(x, y)$ 为空间坐标， $t$ 为时间，

f (q) = [\begin{matrix} ρ u \\ ρ u^{2} + p + \frac{1}{2} {‖B‖}^{2} - B_{1}^{2} \\ ρ u v - B_{1} B_{2} \\ ρ u w - B_{1} B_{3} \\ u (ρ e + p + \frac{1}{2} {‖B‖}^{2}) - B_{1} (u \cdot B) \\ 0 \\ u B_{2} - v B_{1} \\ u B_{3} - w B_{1} \end{matrix}]

，

g (q) = [\begin{matrix} ρ v \\ ρ u v - B_{1} B_{2} \\ ρ v^{2} + p + \frac{1}{2} {‖B‖}^{2} - B_{2}^{2} \\ ρ v w - B_{2} B_{3} \\ v (ρ e + p + \frac{1}{2} {‖B‖}^{2}) - B_{2} (u \cdot B) \\ v B_{1} - u B_{2} \\ 0 \\ v B_{3} - w B_{2} \end{matrix}]

，

其中， $u = {(u, v, w)}^{T}$ 为流速度矢量， $B = (B_{1}, B_{2}, B_{3})^{T}$ 为磁场强度矢量， $ρ$ 为密度， $p$ 为压力，与总能量 $ρ e$ 满足关系式

p = (γ - 1) (ρ e - \frac{ρ}{2} {‖u‖}^{2} - \frac{1}{2} {‖B‖}^{2})

，（2）

其中， $γ$ 为比热比，特征系统参见文献［20］。

二维SWMHD方程由浅水波方程和Maxwell方程组合而成^［6］，在不可压缩性、流动变量二维变化以及 $z$ 方向磁流体静力平衡假设下，在 $z$ 方向对三维理想MHD方程积分，即可得到二维SWMHD方程的守恒律形式（同式（1）），其中，

q = {[h, h u, h v, h B_{1}, h B_{2}]}^{T}

，

\begin{matrix} f (q) = [\begin{matrix} \begin{matrix} h u \\ h u^{2} + \frac{1}{2} g h^{2} - h B_{1}^{2} \end{matrix} \\ h u v - h B_{1} B_{2} \\ \begin{matrix} 0 \\ h u B_{2} - h v B_{1} \end{matrix} \end{matrix}] ， \\ g (q) = [\begin{matrix} \begin{matrix} h v \\ h u v - h B_{1} B_{2} \end{matrix} \\ h v^{2} + \frac{1}{2} g h^{2} - h B_{2}^{2} \\ \begin{matrix} h v B_{1} - h u B_{2} \\ 0 \end{matrix} \end{matrix}] ， \end{matrix}

$h$ 为水深， $g$ 为恒定重力加速度^［6］。

此外，MHD方程在磁场中还需满足无散度条件（又称高斯定理） $\nabla \cdot B = 0$ ，忽略此条件通常造成非线性数值不稳定，得到负压力、负密度等异常结果。在二维情况下，此条件可简化为

\frac{\partial B_{1}}{\partial x} + \frac{\partial B_{2}}{\partial y} = 0

。（3）

2　两类方程的旋转不变性

TORO^［21］提出了二维Euler方程的旋转不变性理论，即对于任意旋转角 $θ$ 和守恒变量 $U = {[ρ, ρ u, ρ v, E]}^{T}$ ，方程均具有旋转不变特性，满足

c o s θ f (U) + s i n θ g (U) = T^{- 1} f (T U)

，（4）

其中， $f (U)$ 和 $g (U)$ 分别为x和y方向的通量函数， $T = T (θ)$ 为旋转矩阵， $T^{- 1}$ 为 $T$ 的逆，有

T = [\begin{matrix} A_{1} & 0 & 0 \\ 0 & A_{2} & 0 \\ 0 & 0 & A_{3} \end{matrix}]

，（5）

其中，

A_{1} = A_{3} = 1

，

A_{2} = [\begin{matrix} c o s θ & s i n θ \\ - s i n θ & c o s θ \end{matrix}]

。

受此理论启发，假设二维理想MHD方程的类旋转矩阵 $T_{I M H D}$ 为 $8 \times 8$ 阶分块对角矩阵。不难发现，三角函数子块 $A_{2}$ 位于旋转矩阵 $T$ 的第2~3行和第2~3列，与变量 $U$ 的速度相应行（第2~3行）对应。由变量 $q$ 的表达式，可知类旋转矩阵 $T_{I M H D}$ 的第2~3行、第2~3列位置应为子块 $A_{2}$ 。考虑流速度 $u$ 和磁场强度 $B$ 均为矢量，不妨借鉴子块 $A_{2}$ 与速度分量的关系，将类旋转矩阵中与磁场强度分量 $B_{1}$ 、 $B_{2}$ 对应的位置记为子块 $A_{2}$ ，这样，类旋转矩阵为

T_{I M H D} = [\begin{matrix} 1 \\ c o s θ & s i n θ \\ - s i n θ & c o s θ \\ 1 \\ 1 \\ c o s θ & s i n θ \\ - s i n θ & c o s θ \\ 1 \end{matrix}]

。（6）

下文将以式（6）为控制方程的类旋转矩阵，给出二维理想MHD方程的旋转不变性定理及证明。

定理1　对于任意旋转角 $θ$ 和变量 $q$ ，二维理想MHD方程的通量函数均满足旋转不变性，即

c o s θ f （ q ） + s i n θ g （ q ） = T_{I M H D}^{- 1} f (T_{I M H D} q)

。（7）

证明　令

\tilde{q} = T_{I M H D} q = {[\tilde{ρ}, \tilde{ρ} \tilde{u}, \tilde{ρ} \tilde{v}, \tilde{ρ} \tilde{w}, \tilde{ρ} \tilde{e}, {\tilde{B}}_{1}, {\tilde{B}}_{2}, {\tilde{B}}_{3}]}^{T}

，

则有

\{\begin{array}{l} \begin{array}{l} \tilde{u} = u c o s θ + v s i n θ, \\ \tilde{v} = - u s i n θ + v c o s θ, \end{array} \\ {\tilde{B}}_{1} = B_{1} c o s θ + B_{2} s i n θ, \\ {\tilde{B}}_{2} = - B_{1} s i n θ + B_{2} c o s θ, \\ \tilde{ρ} = ρ, \tilde{e} = e, \tilde{w} = w, {\tilde{B}}_{3} = B_{3}, \end{array}

易得

\tilde{p} = p

，

{‖\tilde{B}‖}^{2} = {‖B‖}^{2}

，

{‖\tilde{u}‖}^{2} = {‖u‖}^{2}

，

\tilde{u} \cdot \tilde{B} = u \cdot B

。

因此，式（7）右端

\begin{array}{l} T_{I M H D}^{- 1} f (T_{I M H D} q) = \\ [\begin{matrix} ρ \tilde{u} \\ (ρ {\tilde{u}}^{2} + p + \frac{1}{2} {‖B‖}^{2} - {\tilde{B}}_{1}^{2}) c o s θ - (ρ \tilde{u} \tilde{v} - {\tilde{B}}_{1} {\tilde{B}}_{2}) s i n θ \\ (ρ {\tilde{u}}^{2} + p + \frac{1}{2} {‖B‖}^{2} - {\tilde{B}}_{1}^{2}) s i n θ + (ρ \tilde{u} \tilde{v} - {\tilde{B}}_{1} {\tilde{B}}_{2}) c o s θ \\ ρ \tilde{u} w - {\tilde{B}}_{1} B_{3} \\ \tilde{u} (ρ e + p + \frac{1}{2} {‖B‖}^{2}) - {\tilde{B}}_{1} (u \cdot B) \\ - (\tilde{u} {\tilde{B}}_{2} - \tilde{v} {\tilde{B}}_{1}) s i n θ \\ (\tilde{u} {\tilde{B}}_{2} - \tilde{v} {\tilde{B}}_{1}) c o s θ \\ \tilde{u} B_{3} - w {\tilde{B}}_{1} \end{matrix}] \end{array}

。

分析可知，上述列向量第1行

ρ \tilde{u} = ρ (u c o s θ + v s i n θ) = ρ u c o s θ + ρ v s i n θ

。

第2行的

\begin{array}{l} (ρ {\tilde{u}}^{2} + p + {‖B‖}^{2} / 2 - {\tilde{B}}_{1}^{2}) c o s θ = \\ [ρ {(u c o s θ + v s i n θ)}^{2} + p + \\ {‖B‖}^{2} / 2 - (B_{1} c o s θ + B_{2} {s i n θ)}^{2}] c o s θ = \\ ρ u^{2} c o s^{3} θ + ρ v^{2} s i n^{2} θ c o s θ + \\ 2 ρ u v s i n θ c o s^{2} θ + (p + {‖B‖}^{2} / 2) c o s θ - \\ {B_{1}}^{2} c o s^{3} θ - {B_{2}}^{2} s i n^{2} θ c o s θ - 2 B_{1} B_{2} s i n θ c o s^{2} θ ， \end{array}

\begin{array}{l} (ρ \tilde{u} \tilde{v} - {\tilde{B}}_{1} {\tilde{B}}_{2}) s i n θ = \\ [ρ (u c o s θ + v s i n θ) (- u s i n θ + v c o s θ) - \\ (B_{1} c o s θ + B_{2} s i n θ) (- B_{1} s i n θ + \\ B_{2} c o s θ)] s i n θ = \\ - ρ u^{2} s i n^{2} θ c o s θ + ρ v^{2} s i n^{2} θ c o s θ - \\ ρ u v s i n^{3} θ + ρ u v s i n θ c o s^{2} θ + \\ {B_{1}}^{2} s i n^{2} θ c o s θ - {B_{2}}^{2} s i n^{2} θ c o s θ + \\ B_{1} B_{2} s i n^{3} θ - B_{1} B_{2} s i n θ c o s^{2} θ ， \end{array}

整合后可得

\begin{array}{l} (ρ {\tilde{u}}^{2} + p + {‖B‖}^{2} / 2 - {\tilde{B}}_{1}^{2}) c o s θ - \\ (ρ \tilde{u} \tilde{v} - {\tilde{B}}_{1} {\tilde{B}}_{2}) s i n θ = \\ ρ u^{2} c o s θ + ρ u v s i n θ + (p + {‖B‖}^{2} / 2) c o s θ - \\ B_{1}^{2} c o s θ - B_{1} B_{2} s i n θ = \\ (ρ u^{2} + p + {‖B‖}^{2} / 2 - B_{1}^{2}) c o s θ + \\ (ρ u v - B_{1} B_{2}) s i n θ 。 \end{array}

第4行

\begin{array}{l} ρ \tilde{u} w - {\tilde{B}}_{1} B_{3} = ρ w (u c o s θ + v s i n θ) - \\ B_{3} (B_{1} c o s θ + B_{2} s i n θ) = \\ ρ u w c o s θ + ρ v w s i n θ - \\ B_{1} B_{3} c o s θ - B_{2} B_{3} s i n θ = \\ (ρ u w - B_{1} B_{3}) c o s θ + (ρ v w - B_{2} B_{3}) s i n θ 。 \end{array}

至此，列向量第1、第2、第4行的等式关系分别得证。继续验证剩余行，可得旋转不变性关系式恒成立。

证毕。

类似地，给出二维SWMHD方程的旋转不变性定理，此证略。

定理2　对于任意旋转角 $θ$ 和变量 $q$ ，二维SWMHD方程的通量函数均具有旋转不变性，即

c o s θ f （ q ） + s i n θ g （ q ） = T_{S W M H D}^{- 1} f (T_{S W M H D} q)

，（8）

其中，类旋转矩阵

T_{S W M H D} = [\begin{matrix} \begin{matrix} 1 \end{matrix} & \begin{matrix} c o s θ \end{matrix} & \begin{matrix} s i n θ \end{matrix} & \begin{matrix} \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} \end{matrix} \\ - s i n θ & c o s θ & \begin{matrix}  \end{matrix} \\ \begin{matrix} c o s θ & s i n θ \end{matrix} \\ \begin{matrix} - s i n θ & c o s θ \end{matrix} \end{matrix}]

。（9）

3　数值离散

利用旋转不变性式（8）对控制方程离散，嵌入高分辨率通量得到空间离散格式。时间方向采用三阶强稳定Runge-Kutta方法，详见文献［22］。

3.1　空间半离散格式

将计算区域离散为有限个四边形单元，以任一单元为控制体 $V$ ，在 $V$ 上对控制方程积分：

\frac{d}{d t} \iint_{V}^{} q d V + \int_{A} h \cdot n d A = 0

，（10）

其中， $A$ 为 $V$ 的边界，由4个线段 $A_{s} A_{s + 1}$ （s=1，2，3，4； $A_{5} = A_{1}$ ）构成， $h = (f, g)$ ， n 为边界单位外法向量，边界总通量为

\int_{A} h \cdot n d A = \sum_{s = 1}^{4} \int_{A_{s}}^{A_{s + 1}} h \cdot n_{s} d A

。（11）

定义 $θ_{s}$ 为 $x$ 轴正向与界面外法向量 $n_{s}$ 的夹角，有 $n_{s} = (c o s θ_{s}, s i n θ_{s})$ ，则式（11）可改写为

\sum_{s = 1}^{4} \int_{A_{s}}^{A_{s + 1}} h \cdot n_{s} d A = \sum_{s = 1}^{4} \int_{A_{s}}^{A_{s + 1}} [c o s θ_{s} f + s i n θ_{s} g] d A 。

（12）

将式（10）左端第1项看作控制体 $V$ 上 $q$ 的平均值的时间变化率，有

\frac{d}{d t} \iint_{V}^{} q d V = | V | \frac{d}{d t} q

。（13）

将式（11）~式（13）代入式（10），可得

\frac{d}{d t} q = - \frac{1}{| V |} \sum_{s = 1}^{4} \int_{A_{s}}^{A_{s + 1}} [c o s θ_{s} f + s i n θ_{s} g] d A

。（14）

利用控制方程的旋转不变性和式（8），将式（14）转化为

\frac{d}{d t} q = - \frac{1}{| V |} \sum_{s = 1}^{4} \int_{A_{s}}^{A_{s + 1}} {T_{s}}^{- 1} f (T_{s} q) d A

，（15）

右端求和项的每项可近似为

\int_{A_{s}}^{A_{s + 1}} {T_{s}}^{- 1} f (T_{s} q) d A \approx L_{s} {T_{s}}^{- 1} {\tilde{f}}_{s}

，（16）

其中， $L_{s}$ 为线段 $A_{s} A_{s + 1}$ 的长度， ${\tilde{f}}_{s} = f (T_{s} q)$ 为对应的边界通量。最后将式（16）代入式（15），即可得到方程的半离散格式

\frac{d}{d t} q = - \frac{1}{| V |} \sum_{s = 1}^{4} L_{s} {T_{s}}^{- 1} {\tilde{f}}_{s}

。（17）

3.2　高分辨率熵稳定通量

对二维MHD方程，现有熵稳定通量函数为

\{\begin{array}{l} f_{}^{E S} = f_{}^{E C} - R_{x} Λ_{x} S_{x} R_{x}^{T} {[v]}_{x} / 2 ， \\ g_{}^{E S} = g_{}^{E C} - R_{y} Λ_{y} S_{y} R_{y}^{T} {[v]}_{y} / 2 ， \end{array}

（18）

其中， $f_{}^{E C}$ ， $g_{}^{E C}$ 分别为x，y方向的熵守恒通量函数， $- R Λ S R_{}^{T} [v] / 2$ 为Roe型耗散项^［23］， $R$ 为右特征向量矩阵， $Λ$ 为以特征值绝对值为对角元素的对角矩阵， $S$ 为对角放缩矩阵， $R_{}^{T}$ 为 $R$ 的转置， $v$ 为熵变量，具体表达式可参见文献［4，6］。

二维理想MHD方程熵稳定通量中的熵守恒通量函数 $f_{}^{E C}$ 和 $g_{}^{E C}$ 分别为

f_{}^{E C} = [\begin{matrix} \hat{ρ} {\hat{u}}_{1} \\ {\hat{p}}_{1} + \hat{ρ} {\hat{u}}_{1}^{2} + \frac{1}{2} （ {\dot{B}}_{1} + {\dot{B}}_{2} + {\dot{B}}_{3} ） - {\dot{B}}_{1} \\ \hat{ρ} {\hat{u}}_{1} {\hat{v}}_{1} - B_{1} B_{2} \\ \hat{ρ} {\hat{u}}_{1} {\hat{w}}_{1} - B_{1} B_{3} \\ \begin{array}{l} \frac{γ {\hat{u}}_{1} {\hat{p}}_{2}}{γ - 1} + \frac{\hat{ρ} {\hat{u}}_{1}}{2} ({\hat{u}}_{1}^{2} + {\hat{v}}_{1}^{2} + {\hat{w}}_{1}^{2}) + \\ {\hat{u}}_{2} （ {\hat{B}}_{2}^{2} + {\hat{B}}_{3}^{2} ） - {\hat{B}}_{1} （ {\hat{v}}_{2} {\hat{B}}_{2} + {\hat{w}}_{2} {\hat{B}}_{3} ） \end{array} \\ 0 \\ {\hat{u}}_{2} {\hat{B}}_{2} - {\hat{v}}_{2} {\hat{B}}_{1} \\ {\hat{u}}_{2} {\hat{B}}_{3} - {\hat{w}}_{2} {\hat{B}}_{1} \end{matrix}]

，（19）

g_{}^{E C} = [\begin{matrix} \hat{ρ} {\hat{v}}_{1} \\ \hat{ρ} {\hat{u}}_{1} {\hat{v}}_{1} - B_{1} B_{2} \\ {\hat{p}}_{1} + \hat{ρ} {\hat{v}}_{1}^{2} + \frac{1}{2} （ {\dot{B}}_{1} + {\dot{B}}_{2} + {\dot{B}}_{3} ） - {\dot{B}}_{2} \\ \hat{ρ} {\hat{v}}_{1} {\hat{w}}_{1} - B_{2} B_{3} \\ \begin{array}{l} \frac{γ {\hat{v}}_{1} {\hat{p}}_{2}}{γ - 1} + \frac{\hat{ρ} {\hat{v}}_{1}}{2} ({\hat{u}}_{1}^{2} + {\hat{v}}_{1}^{2} + {\hat{w}}_{1}^{2}) + \\ {\hat{v}}_{2} （ {\hat{B}}_{1}^{2} + {\hat{B}}_{3}^{2} ） - {\hat{B}}_{2} （ {\hat{u}}_{2} {\hat{B}}_{1} + {\hat{w}}_{2} {\hat{B}}_{3} ） \end{array} \\ {\hat{v}}_{2} {\hat{B}}_{1} - {\hat{u}}_{2} {\hat{B}}_{2} \\ 0 \\ {\hat{v}}_{2} {\hat{B}}_{3} - {\hat{w}}_{2} {\hat{B}}_{2} \end{matrix}]

。（20）

采用Ismail均值公式^［10］

\{\begin{array}{l} \begin{array}{l} \hat{ρ} = \{\{z_{1}\}\} z_{5}^{l n}, {\hat{p}}_{1} = \frac{\{\{z_{5}\}\}}{\{\{z_{1}\}\}}, \\ {\hat{p}}_{2} = \frac{γ + 1}{2 γ} \frac{z_{5}^{l n}}{z_{1}^{l n}} + \frac{γ - 1}{2 γ} \frac{\{\{z_{5}\}\}}{\{\{z_{1}\}\}}, \end{array} \\ \begin{array}{l} {\hat{u}}_{1} = \frac{\{\{z_{2}\}\}}{\{\{z_{1}\}\}}, {\hat{v}}_{1} = \frac{\{\{z_{3}\}\}}{\{\{z_{1}\}\}}, {\hat{w}}_{1} = \frac{\{\{z_{4}\}\}}{\{\{z_{1}\}\}}, \\ {\hat{u}}_{2} = \frac{\{\{z_{1} z_{2}\}\}}{\{\{{z_{1}}^{2}\}\}}, {\hat{v}}_{2} = \frac{\{\{z_{1} z_{3}\}\}}{\{\{{z_{1}}^{2}\}\}}, {\hat{w}}_{2} = \frac{\{\{z_{1} z_{4}\}\}}{\{\{{z_{1}}^{2}\}\}}, \end{array} \\ \begin{array}{l} {\hat{B}}_{1} = \{\{z_{6}\}\}, {\hat{B}}_{2} = \{\{z_{7}\}\}, {\hat{B}}_{3} = \{\{z_{8}\}\}, \\ {\dot{B}}_{1} = \{\{{z_{6}}^{2}\}\}, {\dot{B}}_{2} = \{\{{z_{7}}^{2}\}\}, {\dot{B}}_{3} = \{\{{z_{8}}^{2}\}\}, \end{array} \\ \begin{array}{l} B_{1} B_{2} = \{\{z_{6} z_{7}\}\}, B_{1} B_{3} = \{\{z_{6} z_{8}\}\}, \\ B_{2} B_{3} = \{\{z_{7} z_{8}\}\} ， \end{array} \end{array}

（21）

其中，参数向量^［4］

\begin{array}{l} z = {[z_{1}, z_{2}, \dots, z_{8}]}^{T} = \\ {[\sqrt[]{\frac{ρ}{p}}, \sqrt[]{\frac{ρ}{p}} u, \sqrt[]{\frac{ρ}{p}} v, \sqrt[]{\frac{ρ}{p}} w, \sqrt[]{ρ p}, B_{1}, B_{2}, B_{3}]}^{T} 。 \end{array}

（22）

二维SWMHD方程的熵守恒通量函数分别为

f_{}^{E C} = [\begin{matrix} \begin{matrix} \{\{h\}\} \{\{u\}\} \\ \begin{array}{l} \{\{h\}\} {\{\{u\}\}}^{2} + \frac{1}{2} g {\{\{h\}\}}^{2} \\ - \{\{h B_{1}\}\} \{\{B_{1}\}\} \end{array} \end{matrix} \\ \{\{h\}\} \{\{u\}\} \{\{v\}\} - \{\{h B_{1}\}\} \{\{B_{2}\}\} \\ \begin{matrix} \{\{h\}\} \{\{u\}\} \{\{B_{1}\}\} - \{\{h B_{1}\}\} \{\{u\}\} \\ \{\{h\}\} \{\{u\}\} \{\{B_{2}\}\} - \{\{h B_{1}\}\} \{\{v\}\} \end{matrix} \end{matrix}]

，（23）

g_{}^{E C} = [\begin{matrix} \begin{matrix} \{\{h\}\} \{\{v\}\} \\ \{\{h\}\} \{\{u\}\} \{\{v\}\} - \{\{h B_{2}\}\} \{\{B_{1}\}\} \end{matrix} \\ \begin{array}{l} \{\{h\}\} {\{\{v\}\}}^{2} + \frac{1}{2} g {\{\{h\}\}}^{2} \\ - \{\{h B_{2}\}\} \{\{B_{2}\}\} \end{array} \\ \begin{matrix} \{\{h\}\} \{\{v\}\} \{\{B_{1}\}\} - \{\{h B_{2}\}\} \{\{u\}\} \\ \{\{h\}\} \{\{v\}\} \{\{B_{2}\}\} - \{\{h B_{2}\}\} \{\{v\}\} \end{matrix} \end{matrix}]

。（24）

然而，熵稳定数值通量在间断位置的过度耗散导致出现不可忽略的抹平现象，为提高格式的分辨率，考虑引入通量限制器。其基本思想是利用限制器函数将高阶通量与低阶通量结合，构造新的通量形式，新通量往往保留了两种通量的优点，分辨率高且鲁棒性强。

下面以 $x$ 方向上的通量 $f$ 为例，构造新的数值通量形式。

分别取 $f_{}^{E C}$ ， $f_{}^{E S}$ 为高阶、低阶数值通量，利用通量限制器函数 $ψ$ 对两者进行组合，得到高分辨率熵稳定通量

\begin{array}{l} f^{E S L} = f_{}^{E S} + ψ (f_{}^{E C} - f_{}^{E S}) = \\ f_{}^{E C} - (I - ψ) (f_{}^{E C} - f_{}^{E S}) = \\ f_{}^{E C} - \frac{1}{2} R_{x} Λ_{x} (I - ψ) S_{x} R_{x}^{T} {[v]}_{x}, \end{array}

（25）

其中， $I$ 为单位矩阵， $ψ = ψ (θ)$ 为对角矩阵， $ψ (θ) = d i a g (φ^{1} (θ), φ^{2} (θ), \dots, φ^{k} (θ))$ ， $φ^{k} (θ)$ 表示第 $k$ 个方程的限制器函数，光滑度指标向量 $θ = [θ^{1}, θ^{2}, \dots, θ^{k}]$ ，在节点 $(x_{i + 1 / 2}, y_{j + 1 / 2})$ 处，有

\begin{array}{l} θ_{i + 1 / 2, j + 1 / 2}^{k} = \{\begin{matrix} (θ^{-})^{k} \\ (θ^{+})^{k} \end{matrix} = \\ \{\begin{matrix} α_{i - 1 / 2, j + 1 / 2}^{k} / α_{i + 1 / 2, j + 1 / 2}^{k}, λ_{i + 1 / 2, j + 1 / 2}^{k} > 0, \\ α_{i + 3 / 2, j + 1 / 2}^{k} / α_{i + 1 / 2, j + 1 / 2}^{k}, λ_{i + 1 / 2, j + 1 / 2}^{k} < 0, \end{matrix} \end{array}

（26）

其中， $λ_{i + 1 / 2, j + 1 / 2}^{k}$ 为控制方程Jacobi矩阵对应的第 $k$ 个特征值， $α_{i + 1 / 2, j + 1 / 2}^{k}$ 为 $α_{i + 1 / 2, j + 1 / 2}$ 的第 $k$ 个分量，有

α_{i + 1 / 2, j + 1 / 2} = L_{i + 1 / 2, j + 1 / 2} Δ_{i + 1 / 2, j + 1 / 2} q

，

$L_{i + 1 / 2, j + 1 / 2}$ 为对应的左特征向量矩阵。显然，可将通量看作低阶通量 $f_{}^{E S}$ 与校正项 $ψ (f_{}^{E C} - f_{}^{E S})$ 的组合，且该校正项的添加在一定程度上抵消了低阶格式的部分耗散量，因此又称校正项为反扩散通量。

通量限制器函数种类繁多，本文选用Minmod限制器^［24］，形式为

φ_{M} (θ) = m a x \{0, m i n [θ, 1]\} = m i n m o d (θ, 1)

。（27）

注1　 $x$ 方向有变量差分算子 $[\cdot] = {(\cdot)}_{i + 1, j} - {(\cdot)}_{i, j}$ ，算术平均算子 $\{\{\cdot\}\} = [{(\cdot)}_{i + 1, j} + {(\cdot)}_{i, j}] / 2$ ，对数平均算子 ${(\cdot)}^{l n} = [{(\cdot)}_{i + 1, j} - {(\cdot)}_{i, j}] / [l n ({(\cdot)}_{i + 1, j}) - l n ({(\cdot)}_{i, j})]$ 。将 ${(\cdot)}_{i + 1, j}$ 替换为 ${(\cdot)}_{i, j + 1}$ ，即可得到 $y$ 方向的算子表达式。

注2　对数平均算子 ${(\cdot)}^{l n}$ 存在分母为零的情形，可参考文献［10］附录B的做法处理。

3.3　高分辨率旋转通量格式

将高分辨率熵稳定数值通量嵌入半离散格式，即可得到熵稳定旋转通量格式。为进一步提高格式的分辨率，在每个时间层上对单元均值 $q_{i + 1 / 2, j + 1 / 2}$ 进行二阶重构^［25］，将重构所得值选为单元界面的左右状态，数值通量函数由 $f_{i, j + 1 / 2} = f (q_{i - 1 / 2, j + 1 / 2}, q_{i + 1 / 2, j + 1 / 2})$ 转化为 $f_{i, j + 1 / 2} = f ({\hat{q}}_{i - 1 / 2, j + 1 / 2}, {\hat{q}}_{i + 1 / 2, j + 1 / 2})$ ，其中，

\begin{array}{l} {\hat{q}}_{i \pm 1 / 2, j + 1 / 2}^{} = q_{i \pm 1 / 2, j + 1 / 2}^{} + \frac{1}{4} (x_{i, j + 1} + x_{i, j} - \\ x_{i \pm 1, j + 1} - x_{i \pm 1, j}) {\tilde{S}}_{i \pm 1 / 2, j + 1 / 2} ， \end{array}

${\tilde{S}}_{i + 1 / 2, j + 1 / 2}$ 为单元均值 $q_{i + 1 / 2, j + 1 / 2}$ 的偏导数 $q_{x}$ 的近似值，有

{\tilde{S}}_{i + 1 / 2, j + 1 / 2} = s i g n (S_{i + 1 / 2, j + 1 / 2}) \cdot \frac{|S_{i + 1 / 2, j + 1 / 2}^{+} \cdot S_{i + 1 / 2, j + 1 / 2}^{-}|}{|S_{i + 1 / 2, j + 1 / 2}^{+}| + |S_{i + 1 / 2, j + 1 / 2}^{-}|}

，

其中，

s i g n (S_{i + 1 / 2, j + 1 / 2}) = s i g n (S_{i + 1 / 2, j + 1 / 2}^{+}) + s i g n (S_{i + 1 / 2, j + 1 / 2}^{-})

，

S_{i + 1 / 2, j + 1 / 2}^{+} = \frac{4 (q_{i + 3 / 2, j + 1 / 2}^{} - q_{i + 1 / 2, j + 1 / 2}^{})}{(x_{i + 2, j + 1} + x_{i + 2, j} - x_{i, j + 1} - x_{i, j})}

，

S_{i + 1 / 2, j + 1 / 2}^{-} = \frac{4 (q_{i + 1 / 2, j + 1 / 2}^{} - q_{i - 1 / 2, j + 1 / 2}^{})}{(x_{i + 1, j + 1} + x_{i + 1, j} - x_{i - 1, j + 1} - x_{i - 1, j})}

。

4　数值实验

用高分辨率旋转通量格式进行数值求解，并与非旋转熵稳定格式的数值结果进行比较，以验证算法的性能。所讨论的算例均采用 $150 \times 150$ 网格，courant friedrichs-leay条件数取0.1。鉴于理想MHD方程变量较多，文中仅展示其部分变量（如压力、密度、Mach数以及磁场强度等），并给出总熵随时间的变化图，包括熵稳定（ES）格式和高分辨率旋转通量（ESLR）格式。

4.1　二维Riemann问题

选定计算区域 $Ω = [0,2] \times [0,2]$ ，将其划分为4个子区域， $I : (1,2) \times (1,2)$ ， $I I : (0,1) \times (1,2)$ ， $I I I : (0,1) \times (0,1)$ ， $I V : (1,2) \times (0,1)$ ，初始时刻各变量的取值见表1。采用Neumann边界条件，得到 $T = 0.1$ 0时刻二维Riemann问题2种格式在区域内的压力和密度，见图1和图2。分析可知，随着时间的推进，子区域 $I$ 和 $I I$ 之间有稀疏波产生， $I I$ 和 $I I I$ 、 $I I I$ 和 $I V$ 之间则有激波产生。可以看到，熵稳定格式（图1）有效抑制了伪振荡的产生，能够准确捕捉解的结构。不足之处在于间断位置过渡带宽度较大，分辨率较低。高分辨率旋转通量格式（图2）保留了熵稳定格式的优点，且耗散量低于熵稳定格式（图3），从而减少了间断位置的抹平现象，过渡带明显比熵稳定格式窄，激波分辨率也得到了显著提高。

表1 二维Riemann问题初始条件

Table 1 The initial condition for 2D Riemann problem

变量	区域
变量	I	II	III	IV
$ρ$	0.930 8	1.030 4	1.000 0	1.888 7
$ρ u$	1.455 7	1.577 4	1.750 0	0.233 4
$ρ v$	-0.463 3	-1.045 5	-1.000 0	-1.742 2
$ρ w$	0.057 5	-0.101 6	0.000 0	0.073 3
$ρ e$	5.083 8	5.781 3	6.000 0	12.999 0
$B_{1}$	0.350 1	0.350 1	0.564 2	0.564 2
$B_{2}$	0.983 0	0.507 8	0.507 8	0.983 0
$B_{3}$	0.350 0	0.157 6	0.253 9	0.491 5

新窗口打开| 下载CSV

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 二维Riemann问题熵稳定格式的数值结果

Fig.1 Numerical results of ES scheme for 2D Riemann problem

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 二维Riemann问题高分辨率旋转通量格式的数值结果

Fig.2 Numerical results of ESLR scheme for 2D Riemann problem

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 二维Riemann问题总熵随时间的变化

Fig.3 Change of total entropy with time for 2D Riemann problem

4.2　First Rotor问题

选择计算区域 $Ω = [0,1] \times [0,1]$ ，取 $γ = 1.4$ ， $r = \sqrt[]{{(x - 0.5)}^{2} + {(y - 0.5)}^{2}}$ ， $u_{0} = 2$ ， $r_{0} = 0.1$ ， $r_{1} = 0.115$ ，初始条件有

\{\begin{array}{l} \begin{array}{l} ρ = 10, (u, v) = \frac{u_{0}}{r_{0}} (- [y - 0.5], [x - 0.5]), \\ r < r_{0}, \end{array} \\ \begin{array}{l} ρ = 1 + 9 f, (u, v) = \frac{f u_{0}}{r} (- [y - 0.5], [x - 0.5]), \\ f = \frac{r_{1} - r}{r_{1} - r_{0}}, r_{0} < r < r_{1}, \end{array} \\ ρ = 1, (u, v) = (0,0), r > r_{1}, \end{array}

$w = 0$ ， $p = 1$ ， $B = {(5 / \sqrt[]{4 π}, 0,0)}^{T}$ ，终止时间 $T = 0.15$ ， $M a = ‖u‖ / a$ ， $a = \sqrt[]{γ p / ρ}$ 。区域中心为一个稠密的、快速旋转的圆柱体，周围区域则是静止的低密度流体，有初始条件均匀的磁场穿过。本文选用周期性边界条件，计算First Rotor问题2种格式在区域内的压力、密度、磁场强度和Mach数，分别如图4和图5所示。可以看到，熵稳定格式的数值结果（图4）整体比较光滑，无振荡产生，但存在一定的抹平现象，涡流边界较为模糊。高分辨率旋转通量格式的数值结果（图5）较熵稳定格式有明显改进，耗散量更少（图6），对间断的捕捉更锐利，边界更清晰，进一步说明了高分辨率旋转通量格式具有良好的性能。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 First Rotor问题熵稳定格式的数值结果

Fig.4 Numerical results of ES scheme for First Rotor problem

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 First Rotor问题高分辨率旋转通量格式的数值结果

Fig.5 Numerical results of ESLR scheme for First Rotor problem

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 First Rotor 问题总熵随时间的变化

Fig.6 Change of total entropy with time for First Rotor problem

4.3　爆炸波问题

选择计算区域 $Ω = [0,1] \times [0,1]$ ， $γ = 5 / 3$ ，令 $r = \sqrt[]{{(x - 0.5)}^{2} + {(y - 0.5)}^{2}}$ ，当 $r < 0.1$ 时， $p = 10$ ，当 $r > 0.1$ 时， $p = 0.1$ 。其余初始条件有 $ρ = 1$ ， $B = {(1 / \sqrt[]{2}, 1 / \sqrt[]{2}, 0)}^{T}$ ， $u = 0$ ， $v = 0$ ， $w = 0$ ，终止时间 $T = 0.15$ 。初始时刻，在圆心为 $(0.5,0.5)$ ，半径为0.1的圆内压强远大于周围区域，导致一个极强的爆炸波快速向外传播，磁场被剧烈压缩，中心区域密度迅速下降。考虑终止时刻爆炸波并未传播到边界，因此边界形式并无过多限制，本文选用Neumann边界条件。图7和图8分别给出了终止时刻爆炸波问题2种格式在计算区域内的压力、密度、磁场强度和Mach数。不难发现，熵稳定格式在间断位置仍存在不可忽略的抹平现象（图7），而旋转不变性以及通量限制器的引入使得此现象得到很大改善（图8），求解效果远优于熵稳定格式，分辨率的提高尤为显著（图9）。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 爆炸波问题熵稳定格式的数值结果

Fig.7 Numerical results of ES scheme for blast wave problem

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 爆炸波问题高分辨率旋转通量格式的数值结果

Fig.8 Numerical results of ESLR scheme for blast wave problem

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 爆炸波问题总熵随时间的变化

Fig.9 Change of total entropy with time for blast wave problem

4.4　Orszag-Tang-like湍流问题

该问题包含MHD湍流的许多重要特征，涉及涡旋系统演化过程中产生的冲击波之间的相互作用。计算区域 $Ω = [0,2 π] \times [0,2 π]$ ，初始条件为

\begin{matrix} h = γ^{2}, u = - s i n y, v = s i n x, \\ B_{1} = - s i n y, B_{2} = s i n 2 x, \end{matrix}

Fig.13 and Fig.14 are the same.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 Orszag-Tang-like湍流问题的数值结果

注左为熵稳定格式，右为高分辨率旋转通量格式，

图13、图14同

Fig.10 Numerical results of Orszag-Tang-like turbulence problem

Note Left is ES scheme， right is ESLR scheme，

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 Orszag-Tang-like 湍流问题总熵随时间的变化

Fig.11 Change of total entropy with time for Orszag-Tang-like turbulence problem

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 2种分离的传导流体问题总熵随时间的变化

Fig.12 Change of total entropy with time for two separated conducting fluids problem

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 2种分离的传导流体问题的数值结果

Fig.13 Numerical results of two separated conducting fluids problem

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 Rotor-like问题的数值结果

Fig.14 Numerical results of Rotor-like problem

其中， $γ = 5 / 3$ ，周期边界， $T = 2$ .00时刻的计算结果见图10和图11。对比分析可知，2种格式均能精确求解此问题，但由于终止时刻间断与复杂流体的相互作用较为明显，高分辨率旋转通量格式的求解效果远超熵稳定格式（图10）。此外，由总熵变化曲线可知，高分辨率旋转通量格式的熵变化量明显低于熵稳定格式（图11）。前1 s内，区域无间断或间断较弱，由限制器机制，高分辨率旋转通量格式所含通量在该时段转化为熵守恒形式，熵增仅出现在强间断形成之后。因此，高分辨率旋转通量格式更符合物理规律，具有更强的稳定性。

4.5　2种分离的传导流体问题

计算区域 $Ω = [- 1,1] \times [- 1,1]$ ，包含以 $(0,0)$ 为圆心、 $r_{0} = 0.3$ 为半径的圆，取 $g = 1$ ，记 $r = \sqrt[]{x^{2} + y^{2}}$ ，初始条件为

[h, u, v, B_{1}, B_{2}]^{T} = \{\begin{array}{l} {[10,0, 0,0.1,0]}^{T}, r \leq r_{0}, \\ {[1,0, 0,1, 0]}^{T}, r > r_{0}, \end{array}

采用Neumann边界，终止时间 $T = 0.15$ ，该时刻区域内的总熵变化以及高度、速度、磁场强度分别如图12和图13所示。可知，高分辨率旋转通量格式在保留鲁棒性的同时，分辨率、激波捕捉能力等均有一定提升，能有效识别流体中相互作用形成的特定间断。

4.6　Rotor-like问题

此问题是理想MHD方程中典型转子问题的推广，在计算区域 $Ω = [- 1,1] \times [- 1,1]$ ，取 $g = 1$ ， $r = \sqrt[]{x^{2} + y^{2}}$ ， $r_{0} = 0.1$ ，初始条件为

[h, u, v, B_{1}, B_{2}]^{T} = \{\begin{array}{l} {[10, - y, x, 0.1,0]}^{T}, r \leq r_{0}, \\ {[1,0, 0,1, 0]}^{T}, r > r_{0}, \end{array}

选用流出边界，终止时刻 $T = 0.2$ 0。区域内高度、速度等数值结果以及总熵变化分别如图14和图15所示。与前文若干算例类似，高分辨率旋转通量格式表现出更好的间断捕捉能力，数值结果几乎等同于文献［5］的格式在 $400 \times 400$ 网格上的结果。

图15

新窗口打开| 下载原图ZIP| 生成PPT

图15 Rotor-like问题总熵随时间的变化

Fig.15 Change of total entropy with time for Rotor-like problem

5　结论

证明了理想MHD方程和SWMHD方程通量函数的旋转不变性，利用通量限制器构造了针对2类MHD方程的高分辨率旋转通量格式。数值实验表明，旋转不变性以及通量限制器的引入使得格式的性能大幅提升。高分辨率旋转通量格式保留了熵稳定格式的特性，减少了数值耗散，消除了激波不稳定性，对解结构（特别是激波）的捕捉更为锐利。后续可进一步引入加权压强函数或加权密度函数，实现旋转角的自适应定义，也可从通量限制器角度开展研究，以期获得更高的分辨率。

http://dx.doi.org/10.3785/j.issn.1008-9497.2024.01.001

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

， GAO

Z M

， DAI

Z H

A 3D staggered Lagrangi-an scheme for ideal magnetohydrodynamics on unstr-uctured meshes

［J］. International Journal for Numerical Methods in Fluids， 2019， 90（11）： 584-602. DOI：10.1002/fld.4736

[本文引用: 1]

[2]

， TANG

High-order low-dissipation targeted ENO schemes for ideal magnetohydrodynamics

［J］. Journal of Scientific Computing， 2019， 80（1）： 692-716. DOI：10.1007/s10915-019-00941-2

[3]

MATTIA

， MIGNONE

A comparison of approxi-mate non-linear Riemann solvers for relativistic MHD

［J］. Monthly Notices of the Royal Astronomical Society， 2022， 510（1）： 481-499. DOI：10.1093/mnras/stab3373

[4]

WINTERS

A R

， GASSNER

G J

Affordable， entropy conserving and entropy stable flux functions for the ideal MHD equations

［J］. Journal of Computational Physics， 2016， 304（1）： 72-108. DOI：10.1016/j.jcp.2015.09.055

[本文引用: 4]

[5]

DUAN

J M

， TANG

H Z

High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

［J］. Journal of Computational Physics， 2021， 431： 110136. DOI：10.1016/j.jcp.2021. 110136

[本文引用: 2]

[6]

WINTERS

A R

， GASSNER

G J

An entropy stable finite volume scheme for the equations of shallow water magnetohydrodynamics

［J］. Journal of Scientific Computing， 2016， 67： 514-539. DOI：10.1007/s10915-015-0092-6

[本文引用: 4]

[7]

KEMM

Roe-type schemes for shallow water magne-tohydrodynamics with hyperbolic divergence cleaning

［J］. Applied Mathematics and Computation， 2016， 272： 385-402. DOI：10.1016/j.amc.2015.05.079

[本文引用: 1]

[8]

TADMOR

The numerical viscosity of entropy stable schemes for systems of conservation laws I

［J］. Mathematics of Computation， 1987， 49（179）： 91-103. DOI：10.1090/s0025-5718-1987-0890255-3

[本文引用: 1]

[9]

LAX

P D

. Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves［M］. Philadelphia： Society for Industrial and Applied Mathematics， 1973： 1-48. doi:10.1137/1.9781611970562.ch1

[本文引用: 1]

[10]

ISMAIL

， ROE

P L

Affordable， entropy-consistent Euler flux functions II： Entropy production at shocks

［J］. Journal of Computational Physics， 2009， 228 （15）： 5410-5436. DOI：10.1016/j.jcp. 2009.04.021

[本文引用: 3]

[11]

郑素佩，建芒芒，封建湖，等.

保号WENO-AO型中心迎风格式

［J］. 计算物理，2022，39（6）：677-686. DOI：10.19596/j.cnki.1001-246x-8507 .

[本文引用: 1]

ZHENG

S P

， JIAN

M M

， FENG

J H

， et al.

Sign preserving WENO-AO-type central upwind schemes

［J］. Chinese Journal of Computational Physics， 2022，39（6）：677-686. DOI：10.19596/j.cnki.1001-246x-8507 .

[本文引用: 1]

[12]

郑素佩，徐霞，封建湖，等.

高阶保号熵稳定格式

［J］. 数学物理学报， 2021， 41（5）： 1296-1310. DOI：10.3969/j.issn.1003-3998.2021.05.005

[本文引用: 1]

ZHENG

S P

， XU

， FENG

J H

， et al.

High-order sign preserving entropy stable schemes

［J］. Acta Mathematica Scientia， 2021， 41（5）： 1296-1310. DOI：10.3969/j.issn.1003-3998.2021.05.005

[本文引用: 1]

[13]

沈亚玲，封建湖，郑素佩，等.

基于一种新型斜率限制器的理想磁流体方程的高分辨率熵相容格式

［J］. 计算物理， 2022， 39（3）： 297-308. DOI：10.19596/j.cnki.1001-246x.8405

[本文引用: 1]

SHEN

Y L

， FENG

J H

， ZHENG

S P

， et al.

High resolution entropy consistent scheme for ideal magne-tohydrodynamics equations based on a new slope limiter

［J］. Chinese Journal of Computational Physics， 2022， 39（3）： 297-308. DOI：10.19596/j.cnki.1001-246x.8405

[本文引用: 1]

[14]

QUIRK

J J

A contribution to the great Riemann solver debate

［J］. International Journal for Numerical Methods in Fluids， 1994， 18（6）： 555-574. DOI：10.1002/fld.1650180603

[本文引用: 2]

[15]

DAVID

， KENETH

， POWELL

Use of a rotated Riemann solver for the two-dimensional Euler equati-ons

［J］. Journal of Computational Physics， 1993， 106（2）： 201-214. DOI：10.1016/S0021-9991（83）71103-4

[本文引用: 1]

[16]

NISHIKAWA

， KITAMURA

Very simple， carbun-cle-free， boundary-layer-resolving， rotated-hybrid Riemann solver

［J］. Journal of Computational Physics， 2008， 227（4）： 2560-2581. DOI：10.1016/j.jcp.2007.11.003

[本文引用: 1]

[17]

ZHANG

， LIU

， CHEN

B S

， et al.

Evaluation of rotated upwind schemes for contact discontinuity and strong shock

［J］. Computers & Fluids， 2016， 134： 11-22. DOI：10.1016/j.compfluid.2016.05.010

[本文引用: 1]

[18]

贾豆，郑素佩.

求解二维Euler方程的旋转通量混合格式

［J］. 应用数学和力学， 2021， 42（2）： 170-179. DOI：10.21656/1000-0887.410196

[本文引用: 1]

JIA

， ZHENG

S P

A hybrid scheme of rotational flux for solving 2D Euler equations

［J］. Applied Mathematics and Mechanics， 2021， 42（2）： 170-179. DOI：10.21656/1000-0887.410196

[本文引用: 1]

[19]

郑素佩，李霄，赵青宇，等.

求解二维浅水波方程的旋转混合格式

［J］. 应用数学和力学， 2022， 43（2）： 176-186. DOI：10.21656/1000-0887.420063

[本文引用: 1]

ZHENG

S P

， LI

， ZHAO

Q Y

， et al.

A rotated mixed scheme for solving 2D shallow water equations

［J］. Applied Mathematics and Mechanics， 2022， 43（2）： 176-186. DOI：10.21656/1000-0887.420063

[本文引用: 1]

[20]

JEFFREY

， TANIUTI

. Non-linear Wave Propagation［M］. New York： Academic Press， 1964： 167-194. DOI：10.1002/zamm.19650450434

[本文引用: 1]

[21]

TORO

. Riemann Solvers and Numerical Methods for Fluid Dynamics［M］. Berlin： Springer， 2013： 105-106.

[本文引用: 1]

[22]

GOTTLIEB

， KETCHESON

D I

， SHU

C W

High order strong stability preserving time discretizations

［J］. Journal of Scientific Computing， 2009， 38（3）： 251-289. DOI：10.1007/s10915-008-9239-z

[本文引用: 1]

[23]

ROE

P L

Entropy Conservation Schemes for Euler Equations

［R］. Lyon： Talk at HYP，2006.

[本文引用: 1]

[24]

SWEBY

P K

， BAINES

M J

On convergence of Roe's scheme for the general non-linear scalar wave equation

［J］. Journal of Computational Physics， 1984， 56（1）： 135-148. DOI：10.1016/0021-9991（84）90087-1

[本文引用: 1]

[25]

TANG

H Z

， TANG

Adaptive mesh methods for one-and two-dimensional hyperbolic conservation laws

［J］. SIAM Journal on Numerical Analysis， 2003， 41（2）： 487-515. DOI：10.1137/s003614290138437x

[本文引用: 1]

A 3D staggered Lagrangi-an scheme for ideal magnetohydrodynamics on unstr-uctured meshes

2019

... 磁流体力学是将经典流体力学和电动力学相结合研究导电流体和电磁场相互作用的学科，在等离子体物理学、天体物理以及流动控制等领域有广泛应用.较为经典的数学模型有理想磁流体（MHD）方程和浅水波磁流体（SWMHD）方程等，当满足适当条件时，在垂直方向对三维理想MHD方程进行积分，即可获得SWMHD方程.鉴于这两类方程均具有类双曲型特点，考虑将双曲型方程的数值求解格式应用于此两类方程.近年来，此类推广已有一系列成果^［1-7］，其中不乏熵稳定格式的推广.熵稳定格式是从物理概念出发，较熵守恒格式含有更多黏性^［8］，其构造理念符合间断位置产生熵增的规律，因而能够有效避免伪振荡的发生，得到唯一且具有物理意义的解^［9］.2016年，WINTERS等^［4，6］以文献［10］为基础，针对此两种方程构造了一类熵稳定数值通量格式，将熵相关理论推广到MHD方程的数值求解，但该通量格式未能实现对间断的高精度捕捉.郑素佩等^［11-12］基于高阶重构陆续得到多种高分辨率熵稳定格式，将求解精度提升到五阶.2021年，DUAN等^［5］提出的高阶精确熵稳定有限差分格式，成功求解了SWMHD方程，验证了熵稳定格式具有良好的性能.沈亚玲等^［13］引入限制器机制构造了一类求解理想MHD方程的高分辨率熵相容格式，该格式能合理控制耗散，有效改善了抹平现象. ...

High-order low-dissipation targeted ENO schemes for ideal magnetohydrodynamics

2019

A comparison of approxi-mate non-linear Riemann solvers for relativistic MHD

2022

Affordable， entropy conserving and entropy stable flux functions for the ideal MHD equations

2016

... 二维理想MHD方程由流体力学中的Navier-Stokes方程和电动力学中的Maxwell方程耦合而成，其守恒律形式^［4］为 ...

... 其中，

f_{}^{E C}

，

g_{}^{E C}

分别为x，y方向的熵守恒通量函数，

- R Λ S R_{}^{T} [v] / 2

为Roe型耗散项^［23］，

R

为右特征向量矩阵，

Λ

为以特征值绝对值为对角元素的对角矩阵，

S

为对角放缩矩阵，

R_{}^{T}

为

R

的转置，

v

为熵变量，具体表达式可参见文献［4，6］. ...

... 其中，参数向量^［4］ ...

High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

2021

... 选用流出边界，终止时刻

T = 0.2

0.区域内高度、速度等数值结果以及总熵变化分别如图14和图15所示.与前文若干算例类似，高分辨率旋转通量格式表现出更好的间断捕捉能力，数值结果几乎等同于文献［5］的格式在

400 \times 400

网格上的结果. ...

An entropy stable finite volume scheme for the equations of shallow water magnetohydrodynamics

2016

... 二维SWMHD方程由浅水波方程和Maxwell方程组合而成^［6］，在不可压缩性、流动变量二维变化以及

z

方向磁流体静力平衡假设下，在

z

方向对三维理想MHD方程积分，即可得到二维SWMHD方程的守恒律形式（同式（1）），其中， ...

...

h

为水深，

g

为恒定重力加速度^［6］. ...

... 其中，

f_{}^{E C}

，

g_{}^{E C}

分别为x，y方向的熵守恒通量函数，

- R Λ S R_{}^{T} [v] / 2

为Roe型耗散项^［23］，

R

为右特征向量矩阵，

Λ

为以特征值绝对值为对角元素的对角矩阵，

S

为对角放缩矩阵，

R_{}^{T}

为

R

的转置，

v

为熵变量，具体表达式可参见文献［4，6］. ...

Roe-type schemes for shallow water magne-tohydrodynamics with hyperbolic divergence cleaning

2016

The numerical viscosity of entropy stable schemes for systems of conservation laws I

1987

1973

Affordable， entropy-consistent Euler flux functions II： Entropy production at shocks

2009

... 采用Ismail均值公式^［10］ ...

... 注2　对数平均算子

{(\cdot)}^{l n}

存在分母为零的情形，可参考文献［10］附录B的做法处理. ...

保号WENO-AO型中心迎风格式

2022

保号WENO-AO型中心迎风格式

2022

高阶保号熵稳定格式

2021

高阶保号熵稳定格式

2021

基于一种新型斜率限制器的理想磁流体方程的高分辨率熵相容格式

2022

基于一种新型斜率限制器的理想磁流体方程的高分辨率熵相容格式

2022

A contribution to the great Riemann solver debate

1994

... 针对二维问题，即使常用的Roe近似Riemann求解器也可能出现包括红斑现象在内的非物理问题^［14］，造成出现此现象的主要原因是近似Riemann求解器存在激波不稳定性问题，该问题可通过旋转通量得以改善或解决^［14］.旋转通量法，尝试将边界通量分解为2个正交分量，利用通量函数的旋转不变性对二维方程做类一维处理，将其简化为一维形式，而前提条件是待解方程必须满足旋转不变性.DAVID等^［15］验证了旋转Riemann求解器在多个迎风角度下的求解效果与高阶MUSCL格式相当.NISHIKAWA等^［16］基于Riemann求解器、结合Roe格式和HLL格式得到一种混合格式，该混合格式能自动切换至少波（HLL）或全波（Roe）求解器，在抑制红斑现象的同时避免过度耗散.ZHANG等^［17］利用加权压强函数定义旋转角度，实现对耗散影响的自适应调整，从而有效消除了由动量扰动造成的激波不稳定.郑素佩等^［18-19］引入了HLL数值通量，构造了求解Euler方程和浅水波方程的混合旋转通量格式，该格式间断捕捉能力好，分辨率高. ...

... ［14］.旋转通量法，尝试将边界通量分解为2个正交分量，利用通量函数的旋转不变性对二维方程做类一维处理，将其简化为一维形式，而前提条件是待解方程必须满足旋转不变性.DAVID等^［15］验证了旋转Riemann求解器在多个迎风角度下的求解效果与高阶MUSCL格式相当.NISHIKAWA等^［16］基于Riemann求解器、结合Roe格式和HLL格式得到一种混合格式，该混合格式能自动切换至少波（HLL）或全波（Roe）求解器，在抑制红斑现象的同时避免过度耗散.ZHANG等^［17］利用加权压强函数定义旋转角度，实现对耗散影响的自适应调整，从而有效消除了由动量扰动造成的激波不稳定.郑素佩等^［18-19］引入了HLL数值通量，构造了求解Euler方程和浅水波方程的混合旋转通量格式，该格式间断捕捉能力好，分辨率高. ...

Use of a rotated Riemann solver for the two-dimensional Euler equati-ons

1993

Very simple， carbun-cle-free， boundary-layer-resolving， rotated-hybrid Riemann solver

2008

Evaluation of rotated upwind schemes for contact discontinuity and strong shock

2016

求解二维Euler方程的旋转通量混合格式

2021

求解二维Euler方程的旋转通量混合格式

2021

求解二维浅水波方程的旋转混合格式

2022

求解二维浅水波方程的旋转混合格式

2022

1964

... 其中，

γ

为比热比，特征系统参见文献［20］. ...

2013

... TORO^［21］提出了二维Euler方程的旋转不变性理论，即对于任意旋转角

θ

和守恒变量

U = {[ρ, ρ u, ρ v, E]}^{T}

，方程均具有旋转不变特性，满足 ...

High order strong stability preserving time discretizations

2009

... 利用旋转不变性式（8）对控制方程离散，嵌入高分辨率通量得到空间离散格式.时间方向采用三阶强稳定Runge-Kutta方法，详见文献［22］. ...

Entropy Conservation Schemes for Euler Equations

2006

... 其中，

f_{}^{E C}

，

g_{}^{E C}

分别为x，y方向的熵守恒通量函数，

- R Λ S R_{}^{T} [v] / 2

为Roe型耗散项^［23］，

R

为右特征向量矩阵，

Λ

为以特征值绝对值为对角元素的对角矩阵，

S

为对角放缩矩阵，

R_{}^{T}

为

R

的转置，

v

为熵变量，具体表达式可参见文献［4，6］. ...

On convergence of Roe's scheme for the general non-linear scalar wave equation

1984

... 通量限制器函数种类繁多，本文选用Minmod限制器^［24］，形式为 ...

Adaptive mesh methods for one-and two-dimensional hyperbolic conservation laws

2003

... 将高分辨率熵稳定数值通量嵌入半离散格式，即可得到熵稳定旋转通量格式.为进一步提高格式的分辨率，在每个时间层上对单元均值

q_{i + 1 / 2, j + 1 / 2}

进行二阶重构^［25］，将重构所得值选为单元界面的左右状态，数值通量函数由

f_{i, j + 1 / 2} = f (q_{i - 1 / 2, j + 1 / 2}, q_{i + 1 / 2, j + 1 / 2})

转化为

f_{i, j + 1 / 2} = f ({\hat{q}}_{i - 1 / 2, j + 1 / 2}, {\hat{q}}_{i + 1 / 2, j + 1 / 2})

，其中， ...

〈

〉