青岛大学，数学与统计学院，山东 青岛

收稿日期：2022年8月9日；录用日期：2022年9月7日；发布日期：2022年9月14日

摘要

在本文研究中，我们针对流体力学中的可压缩欧拉方程组，以熵守恒数值通量为基础，并在此基础上通过添加数值粘性的方式，建立了熵稳定数值通量，以满足熵不等式，最终得到了熵稳定有限差分格式。严格的理论分析以及广泛的数值结果均验证了本格式保持高分辨率和无伪振荡的良好特性。我们相信该方法在流体力学领域会有着相当广泛的应用前景。

关键词

熵守恒，熵稳定，熵不等式，欧拉方程组

Entropy Stable Finite Difference Schemes for Compressible Euler Equations

Xiangyu Zhou, Zhizhuang Zhang, Jinmei Gao, Shouguo Qian^*, Gang Li^*

School of Mathematics and Statistics, Qingdao University, Qingdao Shandong

Received: Aug. 9^th, 2022; accepted: Sep. 7^th, 2022; published: Sep. 14^th, 2022

ABSTRACT

In this article, we aim at the compressible Euler equations model in fluid mechanics, based on the entropy conservation numerical flux, and on this basis, by adding numerical viscosity, we achieve the satisfaction of entropy inequality, and establish the entropy stable numerical flux. Rigorous theoretical analysis and extensive numerical results verify the good characteristics of this method to maintain high resolution and no pseudo oscillation. We believe that this method will have a wide application prospect in the field of fluid mechanics.

Keywords:Entropy Conservation, Entropy Stability, Entropy Inequality, Euler Equations

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

在本文中，我们针对空气动力学中的可压缩欧拉方程组构造高精度的熵稳定有限差分格式。当忽略掉粘性和热通量时，该方程组可以用来描述可压缩流体的流动。对于欧拉方程组，熵条件是热力学第二定律的直接表示。也就是说，在光滑解中守恒的熵应该在冲击时耗散 [1]。不满足熵条件的数值方法可能导致非物理解，例如负密度和负压力。因此提出了各种熵固定技巧来修改Roe格式。这些熵的修正过程涉及到一个人工参数，这个参数在具体问题中有不同的表现形式 [2] [3]。

为了实现熵稳定性，Tadmor引入了熵变量和熵势的概念，设计了一类二阶精度的熵守恒格式 [4] ，并为格式是否熵稳定提供了一个判别准则，即：如果三点格式比熵守恒格式包含更多粘性，则它们是熵稳定的 [5] [6]。根据该准则，可以通过向熵守恒格式添加合适的扩散算子来设计熵稳定的格式。我们通过利用物理粘性和热传导的耗散通量将这一思想应用于欧拉方程组。然而，只有在非常精细的网格上，才会有比较满意的结果。此外，守恒律系统的熵守恒数值通量的计算复杂且代价高昂。为了克服这一障碍，Ismail和Roe为欧拉方程组提出了一种简单且计算高效的熵守恒数值通量。Fjordholm等人提出了一种构造高阶(二阶以上)熵稳定格式的方法 [7] [8]。该方法将高阶熵守恒数值通量和适当的数值扩散算子相结合，涉及满足符号特性的无重构 [9] [10] [11]。

此外，在过去几十年中，欧拉方程组的其他数值格式也得到了很好的发展，例如加权基本无振荡格式，间断伽辽金格式，变分迭代法和无网格方法 [12] - [20]。在各种数值方法中，加权基本无振荡格式和间断伽辽金格式已被研究并得到广泛普及。从成功的基本无振荡格式发展而来，加权基本无振荡格式使用非线性权重在解的平滑区域自动实现高阶精度，并在不连续处实现基本无振荡分辨率 [21]。近年来设计了许多加权基本无振荡格式类型方法 [10] [22]。间断伽辽金格式由Reed和Hill提出 [23] ，并由Cockburn和Shu进一步研究使用高阶龙格–库塔方式进行推进，所得格式称为龙格–库塔间断伽辽金格式 [24]。近年来出现了许多加权基本无振荡格式器 [4] [25] [26] [27] ，构造了不同的限制以抑制冲击附近的振荡。但是这些方法具有严格稳定性的结果仍远未完成，尤其是非线性问题。

在本文中，我们将为可压缩欧拉方程组构造一个熵稳定格式。新的数值通量包含两部分：熵守恒数值通量和数值扩散项。熵守恒部分准确地保留了熵，并且在平滑流中表现良好。然而，它将在冲击附近引起振荡。为了消除这种振荡，我们设计了一个合适的数值扩散项，将其添加到了熵守恒格式上，最终实现了熵稳定格式。

本文的其余部分结构如下。在第2节中，我们简要概述了双曲型守恒律的熵条件。第3节描述了用于求解可压缩欧拉方程组的熵稳定格式的详细设计过程。第4节介绍了我们实施的数值实验，并展示了其结果，旨在评估该格式在具体实践中的良好表现。最后的第5节则是结论。

2. 熵条件

对于拟线性双曲型守恒律方程/方程组的初值问题

$u_t+f{\left(u\right)}_x=0,u\left(x,0\right)=u_0\left(x\right),x\in \Re ,t>0,$ (2.1)

而言，即使初始数据 $u_0\left(x\right)$ 是十分光滑的，解 $u\left(x,t\right)$ 也可能会出现间断。因此，必须从弱解的意义上来考虑这个初值问题的解，也就是说要允许存在更为一般的弱解类。但是，弱解类则会过于宽泛，这就导致了该问题不存在唯一的解。这就需要我们额外的增加条件，来找出物理相关的解。这个条件就被称为熵条件，而据此条件选择出来的解则被称为熵解。熵解u可以看作是具有某种耗散的解 $u^{\epsilon }$ 的极限，即 $u={\lim }_{\epsilon \to 0}{u}^{\epsilon }$，其中 $u^{\epsilon }$ 是下列方程的解

$u_t^{\epsilon }+f{\left(u^{\epsilon }\right)}_x=\epsilon {\left[\beta \left(u^{\epsilon }\right)u_x^{\epsilon }\right]}_x,\beta \left(u^{\epsilon }\right)>0,0<\epsilon <1.$ (2.2)

现以标量守恒律方程 $\left(m=1\right)$ 为例，介绍一般形式的熵条件。如果标量函数对 $\left(\eta \left(u\right),q\left(u\right)\right)$ 满足

$q^{\prime }\left(u\right)={\eta }^{\prime }\left(u\right)f^{\prime }\left(u\right),$ (2.3)

则 $\eta \left(u\right),q\left(u\right)$ 分别被称为熵函数和熵通量，同时 $\left(\eta \left(u\right),q\left(u\right)\right)$ 被称为熵对，其中 $\eta \left(u\right)$ 是u的凸函数。当解u光滑的时候，u就会满足以下等式

$\eta {\left(u\right)}_t+q{\left(u\right)}_x=0.$ (2.4)

如果解u间断，那么u就会在弱意义下满足以下不等式

$\eta {\left(u\right)}_t+q{\left(u\right)}_x\le 0.$ (2.5)

(2.4)和(2.5)式分别被称为熵等式和熵不等式，两者又被合称为熵条件。对于标量守恒律方程，可以找到无穷多的标量函数对 $\left(\eta \left(u\right),q\left(u\right)\right)$，使之满足条件 $q^{\prime }\left(u\right)={\eta }^{\prime }\left(u\right)f^{\prime }\left(u\right)$，成为熵对，即使取定了 $\eta \left(u\right)$ 之后，也可以找到相当一部分能够使 $q^{\prime }\left(u\right)={\eta }^{\prime }\left(u\right)f^{\prime }\left(u\right)$ 成立的 $q\left(u\right)$，并且这些 $q\left(u\right)$ 之间仅仅相差了一个与u没有关系的“常数”。对于 $m\ge 3$ 的守恒律方程组，一般无法根据(2.3)确定出熵对 $\left(\eta \left(u\right),q\left(u\right)\right)$。

3. 数值格式的构造

3.1. 熵守恒格式

现在对初值问题进行守恒型空间离散。考虑 $m=1$ 和均匀的空间网格 $x_j=jh$，其中 $j\in {\rm Z}$，h是空间网格步长，那么半离散格式为

$\frac{\text{d}}{\text{d}t}{u}_j=-\frac{1}{h}\left({\hat{f}}_{j+\frac{1}{2}}-{\hat{f}}_{j-\frac{1}{2}}\right),$ (3.1.1)

其中数值通量 ${\hat{f}}_{j+\frac{1}{2}}$ 为物理通量 $f\left(u\right)$ 的一个适当的近似，并且满足下面的条件。例如，当 ${\hat{f}}_{j+\frac{1}{2}}$ 是两点通量，也就是当 ${\hat{f}}_{j+\frac{1}{2}}={\hat{f}}_{j+\frac{1}{2}}\left(u_j+u_{j+1}\right)$ 时，有 ${\hat{f}}_{j+\frac{1}{2}}\left(u,u\right)=f\left(u\right)$。此时该格式也被称为三点格式，这个条件被称

为相容性条件。

下面将论述以上的半离散格式满足熵等式的条件，或者说熵守恒充分条件。取定熵对 $\left(\eta \left(u\right),q\left(u\right)\right)$。

定义熵变量 $v={\eta }^{\prime }\left(u\right)$，熵势 $\varphi =vu-\eta \left(u\right)$ 和熵通量势 $\psi =vf\left(u\right)-q\left(u\right)$。如果数值通量 ${\hat{f}}_{j+\frac{1}{2}}$ 满足

$\left(v\left(u_j\right)-v\left(u_{j+1}\right)\right){\hat{f}}_{j+\frac{1}{2}}=\psi \left(u_j\right)-\psi \left(u_{j+1}\right),$ (3.1.2)

则半离散格式(3.1.1)是熵守恒的，(3.1.2)就被称为Tadmor的熵守恒数值通量的充分条件 [4]。事实上，

将 $v_j=v\left(u_j\right)$ 左乘以(3.1.1)式，再把 $v_j$ 写成 $v_j=\frac{1}{2}\left(v_j+v_{j+1}\right)+\frac{1}{2}\left(v_j+v_{j-1}\right)$ 的形式，则半离散格式的解满足下式

$\begin{array}{c}{v}_j\frac{\text{d}}{\text{d}t}{u}_j=\frac{\text{d}}{\text{d}t}\eta \left(u_j\right)=-\frac{1}{h}{v}_j\left({\hat{f}}_{j+\frac{1}{2}}-{\hat{f}}_{j-\frac{1}{2}}\right)\\ =-\frac{1}{h}\left(\frac{1}{2}\left(v_j+v_{j+1}\right){\hat{f}}_{j+\frac{1}{2}}+\frac{1}{2}\left(v_j-v_{j+1}\right){\hat{f}}_{j+\frac{1}{2}}\right)\\ +\frac{1}{h}\left(\frac{1}{2}\left(v_j+v_{j-1}\right){\hat{f}}_{j-\frac{1}{2}}+\frac{1}{2}\left(v_j-v_{j-1}\right){\hat{f}}_{j-\frac{1}{2}}\right)\\ =-\frac{1}{h}\left(\frac{1}{2}\left(v_j+v_{j+1}\right){\hat{f}}_{j+\frac{1}{2}}+\frac{1}{2}\left(\psi \left(u_j\right)-\psi \left(u_{j+1}\right)\right)\right)\\ +\frac{1}{h}\left(\frac{1}{2}\left(v_j+v_{j-1}\right){\hat{f}}_{j-\frac{1}{2}}+\frac{1}{2}\left(\psi \left(u_j\right)-\psi \left(u_{j-1}\right)\right)\right)\end{array}$

$\begin{array}{l}=-\frac{1}{h}\left(\frac{1}{2}\left(v_j+v_{j+1}\right){\hat{f}}_{j+\frac{1}{2}}-\frac{1}{2}\left(\psi \left(u_j\right)+\psi \left(u_{j+1}\right)\right)\right)\\ +\frac{1}{h}\left(\frac{1}{2}\left(v_j+v_{j-1}\right){\hat{f}}_{j-\frac{1}{2}}-\frac{1}{2}\left(\psi \left(u_j\right)+\psi \left(u_{j-1}\right)\right)\right)\\ =:-\frac{1}{h}\left({\hat{q}}_{j+\frac{1}{2}}-{\hat{q}}_{j-\frac{1}{2}}\right),\end{array}$ (3.1.3)

其中数值熵通量为

${\hat{q}}_{j+\frac{1}{2}}=\frac{1}{2}\left({\eta }^{\prime }\left(u_j\right)+{\eta }^{\prime }\left(u_{j+1}\right)\right){\hat{f}}_{j+\frac{1}{2}}-\frac{1}{2}\left(\psi \left(u_j\right)+\psi \left(u_{j+1}\right)\right).$ (3.1.4)

此外，令 $u_j=u_{j+1}=u$，利用 ${\hat{f}}_{j+\frac{1}{2}}\left(u,u\right)=f\left(u\right)$ 和 $\psi \left(u\right)$ 的定义，得

${\hat{q}}_{j+\frac{1}{2}}\left(u,u\right)={\eta }^{\prime }\left(u\right)f\left(u\right)-\psi \left(u\right)=q\left(u\right).$ (3.1.5)

所得的数值熵通量是与 $q\left(u\right)$ 相容的。

3.2. 熵稳定格式

当半离散格式(3.1.1)中的数值熵通量 ${\hat{f}}_{j+\frac{1}{2}}$ 满足

$\left({\eta }^{\prime }\left(u_{j+1}\right)-{\eta }^{\prime }\left(u_j\right)\right){\hat{f}}_{j+\frac{1}{2}}\le \psi \left(u_{j+1}\right)-\psi \left(u_j\right),$ (3.2.1)

其中 $\eta$ 为u的凸函数，则该格式满足以下半离散的熵不等式

$\frac{\text{d}}{\text{d}t}\eta \left({\mu }_j\right)\le -\frac{1}{h}\left({\hat{q}}_{j+\frac{1}{2}}-{\hat{q}}_{j-\frac{1}{2}}\right),$ (3.2.2)

并称此时的半离散格式(3.1.1)是熵稳定的。

在熵守恒数值通量后添加一个合适的熵耗散项，是目前构造熵稳定数值通量的最简单的做法。为了

方便区分，将满足熵守恒充分条件的数值通量记为 ${\hat{f}}_{j+\frac{1}{2}}^{EC}={\hat{f}}^{EC}\left(u_j,u_{j+1}\right)$，相应的数值熵通量记为 ${\hat{q}}_{j+\frac{1}{2}}^{EC}$，例如式(3.1.4)中的 ${\hat{q}}_{j+\frac{1}{2}}$。如果取

${\hat{f}}_{j+\frac{1}{2}}^{ES}={\hat{f}}_{j+\frac{1}{2}}^{EC}-\frac{1}{2}{d}_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right),$ (3.2.3)

其中 $d_{j+\frac{1}{2}}\ge 0$，但是不同时取0 (否则就会退化为熵守恒格式)，那么半离散格式就可以写成以下形式

$\frac{\text{d}}{\text{d}t}{u}_j=-\frac{1}{h}\left({\hat{f}}_{j+\frac{1}{2}}^{ES}-{\hat{f}}_{j-\frac{1}{2}}^{ES}\right).$ (3.2.4)

类似前面的做法，有

$\begin{array}{c}{v}_j\frac{\text{d}}{\text{d}t}{u}_j=-\frac{1}{h}{v}_j\left({\hat{f}}_{j+\frac{1}{2}}^{ES}-{\hat{f}}_{j-\frac{1}{2}}^{ES}\right)=-\frac{1}{h}{v}_j\left({\hat{f}}_{j+\frac{1}{2}}^{EC}-{\hat{f}}_{j-\frac{1}{2}}^{EC}\right)\\ +\frac{1}{2h}{v}_j\left(d_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right)-d_{j-\frac{1}{2}}\left(v_j-v_{j-1}\right)\right)\\ =-\frac{1}{h}\left({\hat{q}}_{j+\frac{1}{2}}^{EC}-{\hat{q}}_{j-\frac{1}{2}}^{EC}\right)+S_j,\end{array}$ (3.2.5)

并由加一项和减一项方法有

$\begin{array}{c}{S}_j=\frac{1}{2h}\left(v_j{d}_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right)-v_j{d}_{j-\frac{1}{2}}\left(v_j-v_{j-1}\right)\right)\\ =\frac{1}{4h}\left(\left(v_j-v_{j+1}\right)d_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right)-\left(v_j-v_{j-1}\right)d_{j-\frac{1}{2}}\left(v_j-v_{j-1}\right)\right)\\ +\frac{1}{4h}\left(\left(v_j+v_{j+1}\right)d_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right)-\left(v_j+v_{j-1}\right)d_{j-\frac{1}{2}}\left(v_j-v_{j-1}\right)\right)\\ \le \frac{1}{4h}\left(\left(v_j+v_{j+1}\right)d_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right)-\left(v_j+v_{j-1}\right)d_{j-\frac{1}{2}}\left(v_j-v_{j-1}\right)\right),\end{array}$ (3.2.6)

由此可得，该格式满足半离散的熵不等式，其中

${\hat{q}}_{j+\frac{1}{2}}={\hat{q}}_{j+\frac{1}{2}}^{EC}-\frac{1}{4}\left(v_j+v_{j+1}\right)d_{j+\frac{1}{2}}\left(v_{j+1}-v_j\right).$ (3.2.7)

以上给出了一种半离散的熵稳定格式的构造方式。上述 $d_{j+\frac{1}{2}}$ 可以取为 $\left|{f^{\prime }}_{j+\frac{1}{2}}\right|$ 等。

3.3. 欧拉方程组的熵稳定格式

基于质量、动量、能量守恒定律，空气动力学中的可压缩欧拉方程组具有如下形式

${\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right)}_t+{\left(\begin{array}{c}\rho u\\ \rho u^2+\rho \\ u\left(E+\rho \right)\end{array}\right)}_x=0,$ (3.3.1)

其紧致形式为

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

其中

$U=\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right),F\left(U\right)=\left(\begin{array}{c}\rho u\\ \rho u^2+\rho \\ u\left(E+\rho \right)\end{array}\right).$ (3.3.3)

这里 $\rho ,u,P$ 分别表示流体密度，速度和压力。总能量定义为 $E=\frac{P}{\gamma -1}+\frac{1}{2}\rho u^2$，其中 $\gamma$ 为气体常量。本文中，我们取熵函数 $\eta \left(U\right)=\frac{-\rho s}{\gamma -1}$，熵通量 $q\left(U\right)=\frac{-\rho us}{\gamma -1}$，其中 $s=\ln P-\gamma \ln \rho$ 是物理意义上的熵 [7]。则可求出相应的熵变量为 $V=\frac{\partial \eta \left(U\right)}{\partial U}={\left[\frac{\gamma -s}{\gamma -1}-\frac{\rho u^2}{2P},\frac{\rho u}{P},\frac{-\rho }{P}\right]}^{\text{T}}$，进而计算出熵势为

$\psi \left(U\right)=VU-q\left(U\right)=\rho u$。

Ismail和Roe为欧拉方程组构造了一个熵守恒数值通量 [6] ，其具体步骤如下：首先定义向量

$Z=\left[z_1,z_2,z_3\right]=\left[\sqrt{\frac{\rho }{P}},\sqrt{\frac{\rho }{P}}u,\sqrt{\rho P}\right],$ (3.3.4)

熵守恒数值通量为

$F_{j+\frac{1}{2}}^{EC}={\left[\hat{\rho }\hat{u},\hat{\rho }{\hat{u}}^2+{\hat{P}}_1,\hat{\rho }\hat{u}\hat{H}\right]}^{\text{T}},$ (3.3.5)

其中带有 $\left(\hat{\cdot }\right)$ 的量表示单元界面处满足下列关系的平均值

$\hat{\rho }={\bar{z}}_1{z}_3^{\ln },\hat{u}=\frac{{\bar{z}}_2}{{\bar{z}}_1},{\hat{P}}_1=\frac{{\hat{z}}_3}{{\hat{z}}_1},{\hat{P}}_2=\frac{\gamma +1}{2\gamma }\frac{z_3^{\ln }}{z_1^{\ln }}+\frac{\gamma -1}{2\gamma }\frac{z_3^{\ln }}{z_1^{\ln }},\hat{H}=\frac{1}{2}{\hat{u}}^2+\frac{\gamma }{\gamma -1}\frac{{\hat{P}}_2}{\hat{\rho }},$ (3.3.6)

其中 $a_{j+\frac{1}{2}}^{\ln }=\frac{a_{j+1}-a_j}{\ln \left(a_{j+1}\right)-\ln \left(a_j\right)}$。

应当注意，当 $a_{j+1}\to a_j$ 的时候，由上述对数平均值定义计算出来的数值并不适定。为了克服这个问

题，本文中我们用另一种形式表示对数平均值 [6]。设 $\zeta =\frac{a_{j+1}}{a_j}$，使得 $a_{j+\frac{1}{2}}^{\ln }=\frac{a_{j+1}+a_j}{\ln \zeta }\frac{\zeta -1}{\zeta +1}$，其中

$\ln \left(\zeta \right)=2\left(\frac{1-\zeta }{1+\zeta }+\frac{1}{3}\frac{{\left(1-\zeta \right)}^3}{{\left(1+\zeta \right)}^3}+\frac{1}{5}\frac{{\left(1-\zeta \right)}^5}{{\left(1+\zeta \right)}^5}+\frac{1}{7}\frac{{\left(1-\zeta \right)}^7}{{\left(1+\zeta \right)}^7}+O\left({\zeta }^9\right)\right)$。

由此取得数值形式良好的对数平均值。

对熵守恒数值通量适当的加入某些粘性项，可以使其满足离散的熵不等式，构造出熵稳定数值通量，来避免伪振荡的产生。熵稳定数值通量表示如下

$f_{j+\frac{1}{2}}^{ES}=f_{j+\frac{1}{2}}^{EC}-\frac{1}{2}R\left|\Lambda \right|R^{\text{T}}\left(V_{j+1}-V_j\right),$ (3.3.7)

其中 $\Lambda =diag\left({\lambda }_{j+1/2}^1,{\lambda }_{j+1/2}^2,{\lambda }_{j+1/2}^3\right)=diag\left(u-c,u,u+c\right)$，$c=\sqrt{\left(\gamma -1\right)\left(H-u^2/2\right)}$，

$R=\left(R_1,R_2,R_3\right)=\left(\begin{array}{ccc}1& 1& 1\\ u-c& u& u+c\\ H-uc& u^2/2& H+uc\end{array}\right).$

其中 $R_1,R_2,R_3$ 分别为对应于特征值 ${\lambda }^1,{\lambda }^2,{\lambda }^3$ 的右特征向量， $u,c$ 均取单元界面处即半点处的值。

4. 数值结果

在本节中，我们将通过计算不同算例来测试所提出的熵稳定格式的性能。对于时间的计算，我们采用经典的三阶Runge-Kutta方法，

$\begin{array}{l}{U}_{\tau }^{\left(1\right)}=U_{\tau }^n+\Delta t\mathcal{L}\left(U_{\tau }^n\right),\\ U_{\tau }^{\left(2\right)}=\frac{3}{4}{U}_{\tau }^n+\frac{1}{4}\left(U_{\tau }^{\left(1\right)}+\Delta t\mathcal{L}\left(U_{\tau }^{\left(1\right)}\right)\right),\\ U_{\tau }^{n+1}=\frac{1}{3}{U}_{\tau }^n+\frac{2}{3}\left(U_{\tau }^{\left(2\right)}+\Delta t\mathcal{L}\left(U_{\tau }^{\left(2\right)}\right)\right).\end{array}$

为了保持数值稳定性取 $CFL=0.18$。

4.1. Lax问题

初始条件为

$\left(p,u,\rho \right)\left(x,0\right)=\{\begin{array}{l}\left(0.445,0698,3.528\right)\text{if}x\le 0,\\ \left(0.5,0,0.571\right)\text{if}x>0,\end{array}$

其中计算区域为 $\left[-5,5\right]$，将该算例计算至 $t=1.3$，数值结果见图1。

图1. 在4.1节例子的数值结果

4.2. Sod问题

初始条件为

$\left(p,u,\rho \right)\left(x,0\right)=\{\begin{array}{l}\left(1,0,1\right)\text{if}x\le 0,\\ \left(0.125,0,0.1\right)\text{if}x>0,\end{array}$

其中计算区域为 $\left[-5,5\right]$，将该算例计算至 $t=2$，数值结果见图2。

图2. 在4.2节例子的数值结果

4.3. Shu-Osher问题

初始条件为

$\left(p,u,\rho \right)\left(x,0\right)=\{\begin{array}{l}\left(3.857143,2.629369,10.333333\right)\text{if}x<-4,\\ \left(1+0.2\sin \left(5x\right),0,1\right)\text{if}x\ge -4,\end{array}$

其中计算区域为 $\left[-5,5\right]$，将该算例计算至 $t=1.8$，数值结果见图3。

图3. 在4.3节例子的数值结果

4.4. 123问题

初始条件为

$\left(p,u,\rho \right)\left(x,0\right)=\{\begin{array}{l}\left(1,-2,0.4\right)\text{if}x<0.5,\\ \left(1,2,0.4\right)\text{if}x\ge 0.5,\end{array}$

其中计算区域为 $\left[0,1\right]$，将该算例计算至 $t=0.15$，数值结果见图4。

图4. 在4.4节例子的数值结果

4.5. 修正的冲击/湍流相互作用问题

初始条件为

$\left(p,u,\rho \right)\left(x,0\right)=\{\begin{array}{l}\left(1.515695,0.523346,1.805\right)\text{if}x\le 4.5,\\ \left(1+0.1\sin \left(20\pi x\right),0,1\right)\text{if}x>4.5,\end{array}$

其中计算区域为 $\left[-5,5\right]$，将该算例计算至 $t=5$，数值结果见图5。

图5. 在4.5节例子的数值结果

通过上述多个算例的数值结果，可以清楚地发现在通过添加粘性项的方式，所得到的熵稳定格式，应用于可压缩欧拉方程组的计算，均可以使数值结果很好地保持高分辨率，在间断区域保持了陡峭的间断过渡，而且很好地解决了接触间断。一些算例中产生的冲击，在使用熵稳定格式之后，也得到了很好的解决。最为关键的是，所得到的数值结果图像均与精确解图像很好地拟合在一起，并且没有伪振荡现象发生。

5. 结论

本文针对欧拉方程组，通过在熵守恒格式的基础上通过添加适当的数值熵粘性的方法构造了熵稳定格式，并将熵稳定格式运用在算例上来检验其特性。欧拉方程组的多个算例的数值结果表明，应用该格式计算得到的数值解能够与精确解很好地拟合，数值结果能够保持高分辨率，避免伪振荡的发生，并且对于间断解保持陡峭的间断过渡。

致谢

本研究得到了山东省自然科学基金面上项目(No. ZR2021MA072)的资助。

文章引用

周翔宇,张志壮,高金梅,钱守国,李 刚. 可压缩欧拉方程组的熵稳定有限差分格式
Entropy Stable Finite Difference Schemes for Compressible Euler Equations[J]. 应用数学进展, 2022, 11(09): 6331-6341. https://doi.org/10.12677/AAM.2022.119669

参考文献