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

收稿日期：2022年7月11日；录用日期：2022年8月4日；发布日期：2022年8月15日

摘要

在本文中，我们针对流体力学中的Burgers方程以及浅水波方程，构造了高精度熵稳定格式。我们首先以熵守恒数值通量为基础，通过添加适当的数值熵粘性的方式，构造了熵稳定数值通量，实现了熵不等式，最终建立了熵稳定数值格式。广泛的数值结果均验证了本格式保持高分辨率和无伪振荡的良好特性。我们相信该格式在流体力学领域会有着相当广泛的应用前景。

关键词

熵守恒，熵稳定，熵不等式，浅水波方程，Burgers方程

An Entropy Stable Scheme for Shallow Water Equations Based on Entropy Conservative Scheme

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

School of Mathematics and Statistics, Qingdao University, Qingdao Shandong

Received: Jul. 11^th, 2022; accepted: Aug. 4^th, 2022; published: Aug. 15^th, 2022

ABSTRACT

In this article, we aim at building high-order entropy stable scheme for the Burgers equation and shallow water equation in fluid mechanics. Based on the entropy conservative numerical flux, we first construct the entropy stable numerical flux by adding appropriate numerical entropy viscosity, then achieve the entropy inequality, and achieve the entropy stable finite difference scheme eventually. Extensive numerical results illustrate that the resulting scheme keeps high resolution and is free of spurious oscillations. We believe that this scheme will have a wide application prospect in the field of fluid mechanics.

Keywords:Entropy Conservation, Entropy Stable, Entropy Inequality, Shallow Water Equation, Burgers Equation

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. 引言

在本文中，我们主要通过高精度数值格式来求解拟线性双曲型守恒律方程/方程组的初值问题

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

该初值问题具有以下特点，即使给出了充分光滑的初始条件 $u_0\left(x\right)$，相应的解 $u\left(x,t\right)$ 也会出现间断，这就导致了无法得到古典解。因此，必须从弱解的意义上来考虑该初值问题的解，也就是说要允许存在更为一般的弱解类。但是，弱解类又过于宽泛，这就导致了该初值问题的弱解可能不唯一。因此，我们需要增加额外的条件，来筛选出物理相关的解 [1]。此处增加的条件就被称为熵条件，而据此条件选择出来的物理相关的解也被称为熵解。熵解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.$ (1.2)

浅水波方程是用来描述浅水环境下各种类型流体运动规律的重要数学模型。它可以用于刻画湖泊、河流、灌溉渠和沿海地区的流体流动。在处理垂直运动尺度明显小于水平运动尺度这样的近似浅水环境下，通过简化流体动力学中的Navier-Stokes方程，可以得到浅水流体的平均深度的连续性方程和动量方程 [2] [3]。本文研究的目的在于设计一个有效的熵稳定格式，用于求解一维浅水波方程

$\frac{\partial }{\partial t}\left[\begin{array}{c}h\\ hu\end{array}\right]+\frac{\partial }{\partial x}\left[\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\end{array}\right]=0,$ (1.3)

这里， $h=h\left(x\right)$ 表示流体的高度，u表示速度，g表示重力加速度。该方程的一般紧致向量形式为

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

其中， $U={\left[h,hu\right]}^{\text{T}}$ 是质量和动量的守恒变量， $F\left(U\right)={\left[hu,h{u}^2+\frac{1}{2}g{h}^2\right]}^{\text{T}}$ 是流通量。

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

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

本文针对浅水波方程的空间导数离散过程，提出了一种具有高分辨的熵稳定数值通量表达式；并将该格式分别应用于一维浅水波方程的数值算例中，数值结果表明该格式具有高分辨率和基本无震荡的良好特点，是模拟浅水波方程的理想格式。

本文的其余部分结构如下。在第2节中，我们简要概述了双曲型守恒定律的相关理论。第3节描述了用于求解Burgers方程和浅水波方程的熵稳定格式的详细设计过程。第4节介绍了我们所做的数值实验及结果，旨在评估该格式在实践中的表现。最后的第5节则是结论。

2. 熵条件

首先以标量守恒律方程( $m=1$ )为例，介绍一般形式的熵条件。如果标量函数对 $\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.1)

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

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

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

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

式(2.2)与(2.3)又分别被称为熵等式和熵不等式，或者通称熵条件。对于标量守恒律方程( $m=1$ )，可以找到无穷多的标量函数对 $\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没有关系的“常数”。

3. 数值格式

3.1. 熵守恒格式

现在对初值问题(1.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)$，再用 $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).\end{array}$ (3.1.2)

由此可知，如果 ${\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.3)

这实际上是Tadmor的熵守恒通量的充分条件 [4]，则格式是熵守恒的，或者满足下面的半离散熵等式

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\eta \left(u_j\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)\\ =-\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.4)

其中数值熵通量

${\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.5)

是与 $q\left(u\right)$ 相容的。事实上，令 $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.6)

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}{l}{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}{l}{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. Burgers方程

一维情况下，Burgers方程为 $u_t+{\left(\frac{1}{2}{u}^2\right)}_x=0$，如果取熵函数为 $\eta =\frac{1}{2}{u}^2$，则可根据充分条件(3.1.3)，得到一维Burgers方程的熵守恒数值通量为

${\hat{F}}_{j+\frac{1}{2}}=\frac{1}{6}\left(u_{j+1}^2+u_j^2+u_{j+1}{u}_j\right).$ (3.3.1)

在(3.3.1)的基础上添加粘性项后，有

${\tilde{F}}_{j+\frac{1}{2}}=\frac{1}{6}\left(u_{j+1}^2+u_j^2+u_{j+1}{u}_j\right)-\frac{1}{2}\left|u_j\right|\left(u_{j+1}-u_j\right).$ (3.3.2)

类似地，二维情况下，控制方程为 $u_t+{\left(\frac{1}{2}{u}^2\right)}_x+{\left(\frac{1}{2}{u}^2\right)}_y=0$，取熵函数为 $\eta =\frac{1}{2}{u}^2$，根据充分条件(3.1.3)

可得到二维Burgers方程的熵守恒数值通量如下

${\hat{F}}_{i+\frac{1}{2},j}=\frac{1}{6}\left(u_{i+1,j}^2+u_{i,j}^2+u_{i+1,j}{u}_{i,j}\right),$ (3.3.3)

${\hat{G}}_{i,j+\frac{1}{2}}=\frac{1}{6}\left(u_{i,j+1}^2+u_{i,j}^2+u_{i,j+1}{u}_{i,j}\right).$ (3.3.4)

在(3.3.3)和(3.3.4)的基础上添加粘性项后，有

${\tilde{F}}_{i+\frac{1}{2},j}=\frac{1}{6}\left(u_{i+1,j}^2+u_{i,j}^2+u_{i+1,j}{u}_{i,j}\right)-\frac{1}{2}\left|u_{i,j}\right|\left(u_{i+1,j}-u_{i,j}\right),$ (3.3.5)

${\tilde{G}}_{i,j+\frac{1}{2}}=\frac{1}{6}\left(u_{i,j+1}^2+u_{i,j}^2+u_{i,j+1}{u}_{i,j}\right)-\frac{1}{2}\left|u_{i,j}\right|\left(u_{i,j+1}-u_{i,j}\right).$ (3.3.6)

其半离散格式为

$\frac{\text{d}}{\text{d}t}{u}_{i,j}=-\frac{1}{\text{d}x}\left({\tilde{f}}_{i+\frac{1}{2},j}-{\tilde{f}}_{i-\frac{1}{2},j}\right)-\frac{1}{\text{d}y}\left({\tilde{g}}_{i,j+\frac{1}{2}}-{\tilde{g}}_{i,j-\frac{1}{2}}\right).$ (3.3.7)

3.4. 浅水波方程

对于一维的浅水波方程(1.3)，先确定函数对 $\left(\eta \left(U\right),q\left(U\right)\right)$。首先，取定熵函数 $\eta \left(U\right)=\frac{1}{2}h{u}^2+\frac{1}{2}g{h}^2$，由一般形式的熵条件 $\nabla q^{\text{T}}=\nabla {\eta }^{\text{T}}\frac{\partial F\left(U\right)}{\partial U}$，其中

$\nabla \eta ={\left(\frac{\partial \eta }{\partial h},\frac{\partial \eta }{\partial \left(hu\right)}\right)}^{\text{T}}={\left(gh-\frac{u^2}{2},u\right)}^{\text{T}},$ (3.4.1)

$\frac{\partial F\left(U\right)}{\partial U}=\left(\begin{array}{cc}\frac{\partial \left(hu\right)}{\partial h}& \frac{\partial \left(hu\right)}{\partial \left(hu\right)}\\ \frac{\partial \left(h{u}^2+\frac{1}{2}g{h}^2\right)}{\partial h}& \frac{\partial \left(h{u}^2+\frac{1}{2}g{h}^2\right)}{\partial \left(hu\right)}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ gh-u^2& 2u\end{array}\right).$ (3.4.2)

经过计算可知， $\nabla q={\left(ghu-u^3,gh+\frac{3}{2}{u}^2\right)}^{\text{T}}$，则可取 $q\left(U\right)=\frac{1}{2}h{u}^3+g{h}^{2}u$，并且易知熵变量 $V=\nabla \eta ={\left(gh-\frac{1}{2}{u}^2,u\right)}^{\text{T}}$，则可根据熵变量计算出熵势 $\phi =V^{\text{T}}U-\eta \left(U\right)=\frac{1}{2}g{h}^2$，熵通量势 $\psi =V^{\text{T}}F\left(U\right)-q\left(U\right)=\left(gh-\frac{1}{2}{u}^2,u\right)\left(\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\end{array}\right)-\left(\frac{1}{2}h{u}^3+g{h}^{2}u\right)=\frac{1}{2}g{h}^{2}u$，

充分条件为

${\left(V_j-V_{j+1}\right)}^{\text{T}}{\hat{F}}_{j+\frac{1}{2}}=\left({\psi }_j-{\psi }_{j+1}\right),$ (3.4.3)

即

$\left(\left(g{h}_j-\frac{1}{2}{u}_j^2\right)-\left(g{h}_{j+1}-\frac{1}{2}{u}_{j+1}^2\right),u_j-u_{j+1}\right)\left(\begin{array}{l}{\hat{F}}_{j+\frac{1}{2}}^{\left(1\right)}\\ {\hat{F}}_{j+\frac{1}{2}}^{\left(2\right)}\end{array}\right)=\frac{1}{2}g{h}_j^2{u}_j-\frac{1}{2}g{h}_{j+1}^2{u}_{j+1}.$ (3.4.4)

通过计算，得到如下熵守恒数值通量

$F_{j+\frac{1}{2}}^{EC}=\left(\begin{array}{c}\overline{h_{j+\frac{1}{2}}}\overline{u_{j+\frac{1}{2}}}\\ \frac{1}{2}g\overline{h_{j+\frac{1}{2}}^2}+\overline{h_{j+\frac{1}{2}}}{\left(\overline{u_{j+\frac{1}{2}}}\right)}^2\end{array}\right),$ (3.4.5)

其中 $\overline{u_{j+\frac{1}{2}}}=\frac{1}{2}\left(u_j+u_{j+1}\right)$。

熵稳定数值通量为

${\tilde{F}}^{ES}{}_{j+\frac{1}{2}}={\hat{F}}_{j+\frac{1}{2}}-\frac{1}{2}{D}_{j+\frac{1}{2}}\left(V_{j+1}-V_j\right),$ (3.4.6)

其中 $D_{j+\frac{1}{2}}=R_{j+\frac{1}{2}}{\Lambda }_{j+\frac{1}{2}}{R}_{j+\frac{1}{2}}^{\text{T}}$，$\Lambda =diag\left(\left|{\lambda }_1\right|,\left|{\lambda }_2\right|\right)$，${\lambda }_1=u-c$，${\lambda }_2=u+c$ 为矩阵 $\frac{\partial F\left(U\right)}{\partial U}$ 的特征值，其中 $c=\sqrt{gh}$，$R=\left(R_1,R_2\right)$，$R_1,R_2$ 为分别对应于特征值 ${\lambda }_1,{\lambda }_2$ 的右特征向量， $R_1={\left(1,u-c\right)}^{\text{T}}$，$R_2={\left(1,u+c\right)}^{\text{T}}$，另外，以上均使用半点处 $x_{j+\frac{1}{2}}$ 的信息，

$u_{j+\frac{1}{2}}=\frac{\sqrt{h_j}{u}_j+\sqrt{h_{j+1}}{u}_{j+1}}{\sqrt{h_j}+\sqrt{h_{j+1}}},h_{j+\frac{1}{2}}=\frac{1}{2}\left(h_j,h_{j+1}\right).$ (3.4.7)

4. 数值结果

在本节中，我们将通过多个算例来测试所提出的熵稳定格式的良好性能。在整个计算过程中，我们取重力加速度g为9.812，为了保持数值稳定性取 $CFL=0.18$。关于时间离散，我们采用经典的三阶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}$

4.1. 一维Burgers方程

初始条件为

$\begin{array}{l}{u}_t+{\left(\frac{1}{2}{u}^2\right)}_x=0,\\ u\left(x,0\right)=0.5+\sin \left(\pi x\right),\end{array}$

计算区域为 $\left[0,1\right]$，将该算例计算至 $t=\frac{1.5}{\pi }$，数值结果见图1。

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

4.2. 二维Burgers方程

初始条件为

$\begin{array}{l}{u}_t+{\left(\frac{1}{2}{u}^2\right)}_x+{\left(\frac{1}{2}{u}^2\right)}_y=0,\\ u\left(x,y,0\right)=0.5+\sin \left(\pi \frac{x+y}{2}\right),\end{array}$

计算区域为 $\left[-2,2\right]\times \left[-2,2\right]$，将该算例计算至 $t=\frac{1.5}{\pi }$，数值结果见图2。

通过Burgers方程的图像可以看出，熵稳定格式的数值结果对间断解保持陡峭的间断过渡，而且所得到的数值解与精确解拟合的比较好。

接下来，我们分别考虑浅水波方程的不同的Riemann问题。

图2. 在4.2节例子的数值结果。3D图(左)，沿着x = y的截面图(右)

4.3. 一维浅水波方程：Riemann问题1

初始条件为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(1,2.5\right)x\le 10,\\ \left(0.1,0\right)\text{else},\end{array}$

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

图3. 在4.3节例子的数值结果。水体深度(左)，水体速度(右)

4.4. 一维浅水波方程：Riemann问题2

初始条件为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(1,0\right)x\le 10,\\ \left(0.5,0\right)\text{else},\end{array}$

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

图4. 在4.4节例子的数值结果。水体深度(左)，水体速度(右)

4.5. 一维浅水波方程：Riemann问题3

初始条件为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(1,0\right)x\le 0,\\ \left(0.5,0\right)\text{else},\end{array}$

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

图5. 在4.5节例子的数值结果。水体深度(左)，水体速度(右)

4.6. 一维浅水波方程：Riemann问题4

初始条件为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(1,-5\right)x\le 25,\\ \left(1,-1\right)\text{else},\end{array}$

计算区域为 $\left[0,50\right]$，将该算例计算至 $t=2.5$，数值结果见图6。

图6. 在4.6节例子的数值结果。水体深度(左)，水体速度(右)

通过上述四个算例的数值结果，可以清楚地发现熵稳定格式的数值结果能够很好地保持高分辨率，数值结果均与精确解很好地匹配，且不存在伪震荡。

5. 结论

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

基金项目

本研究得到了国家自然科学基金面上项目(No. 11771228)的资助。

文章引用

张志壮,周翔宇,高金梅,李 刚. 一种基于熵守恒格式的浅水波方程的熵稳定格式
An Entropy Stable Scheme for Shallow Water Equations Based on Entropy Conservative Scheme[J]. 应用数学进展, 2022, 11(08): 5648-5659. https://doi.org/10.12677/AAM.2022.118596

参考文献