浅水波方程组的熵稳定有限体积格式

Entropy Stable Finite Volume Scheme for Shallow Water Equations

摘要: 本文针对非平底地形上的浅水波方程组，提出了一种高精度熵稳定有限体积格式。首先，我们构造了一个具有二阶精度的well-balanced的半离散熵守恒格式，该格式满足给定熵对的熵恒等式，并精确地保持静水稳态。本文的关键点是使通量梯度和源项的离散化相互匹配；其次，二阶熵守恒格式的可负担熵守恒通量对于最终的高阶格式也是至关重要的。然后，以二阶熵守恒格式为基本模块，实现了高阶well-balanced的半离散熵守恒格式。第三，通过在现有的熵守恒格式中添加适当的耗散项，提出了高阶精度的well-balanced的半离散熵稳定格式，该耗散项基于标度熵变量的加权本质非振荡重构，以克服熵守恒的数值振荡。最后，使用Runge-Kutta方法对半离散格式进行时间积分，以实现最终格式。用大量的数值结果说明所提出的格式满足离散熵不等式，具有良好的平衡性，保持了光滑解的真正高阶精度，并且能够很好地捕获稳态下的小扰动。

Abstract: This article develops a high-order accurate entropy stable finite volume scheme for the shallow wa-ter equations over non-flat bottom topography. Firstly, we construct a second-order accurate well-balanced semi-discrete entropy conservative scheme, which satisfies the entropy identity for the given entropy pair and maintains the still water steady states exactly. The key idea is to make both discretizations for the flux gradient and the source term match each other. Secondly, afforda-ble entropy conservative fluxes of the second-order entropy conservative schemes are also critical for the ultimate high-order scheme. Then, high-order well-balanced semi-discrete entropy con-servative scheme is achieved by taking the second-order entropy conservative schemes as a build-ing block. Thirdly, the high-order accurate well-balanced semi-discrete entropy stable schemes are proposed by adding a suitable dissipation term to the existing entropy conservative scheme based on the weighted essentially non-oscillatory reconstruction of the scaled entropy variables to sup-press the numerical oscillations of the entropy conservative scheme. Finally, the semi-discrete scheme is integrated in time using the Runge-Kutta approach to achieve the eventual scheme. Ex-tensive numerical results strongly indicate that the proposed scheme satisfies the discrete entropy inequality, is well-balanced, keeps the genuine high-order accuracy for smooth solutions, and can well capture small perturbations according to the steady state.

1. 引言

浅水波方程组(SWE)广泛用于描述在重力和底部地形影响下的浅水流(如河流和海岸流、湖泊流、潮汐等)，其垂直维度远小于水平维度。其一维形式为

$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}=-S\left(U,b\right),$ (1)

其中守恒变量 $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}}$ ，以及源项 $S\left(U,b\right)={\left(0,gh{b}_x\right)}^{\text{T}}$ 。其中，

$h\left(x,t\right)$ 表示流体的高度， $u\left(x,t\right)$ 表示流体速度的深度平均， $b\left(x\right)$ 表示底部地形， $g=9.812$ 表示重力常数。特别的，方程组(1)满足下面的静水稳态方程

$u=0,h+b=\text{Constant}\text{.}$ (2)

实际上，由于控制方程的高度非线性，方程组(1)的理论分析非常困难。因此，高阶格式的数值求解是解决这些问题的有效方法，并在计算流体力学领域引起了广泛关注 [1] [2] [3] 。SWE高阶格式的代表性研究(1)主要包括：动理学格式 [4] 、气体动理学格式 [5] 、中心迎风格式 [6] 、加权基本非振荡(WENO)格式 [7] - [16] 、Hermite WENO格式 [17] 、中心格式 [18] [19] [20] [21] 、Runge-Kutta间断Galerkin (RKDG)方法 [22] [23] [24] 、ADER (Arbitrary DERivatives in space and time)格式 [25] [26] 、谱元法 [27] 、Godunov方法 [28] [29] 、无网格伽辽金法 [30] 、中心–迎风格式 [31] [32] 、ADER间断伽辽金(ADER-DG)法 [33] 。

此外，非线性双曲守恒律(如方程组(1))一般没有精确解。因此，需要定义弱解。不幸的是，这种弱解可能不是唯一的。事实上，熵条件从所有可能的弱解中筛选物理相关解方面起着重要作用。此外，物理相关解也称为熵解。

定义1. (熵函数 [34] )严格凸函数 $\eta \left(U\right)$ 称为方程组(1)的熵函数，如果存在相关的熵通量 $q\left(U\right)$ ，使得

$q^{\prime }\left(U\right)=V^{\text{T}}{F}^{\prime }\left(U\right),$ (3)

其中 $V={\eta }^{\prime }{\left(U\right)}^{\text{T}}$ 称为熵变量， $\left(\eta ,q\right)$ 称为熵对。

对于给定的熵对 $\left(\eta ,q\right)$ ，将(1)的两端乘以 $V^{\text{T}}$ ，其光滑解满足以下熵恒等式

$\frac{\partial \eta \left(U\right)}{\partial t}+\frac{\partial q\left(U\right)}{\partial x}=-V^{\text{T}}S\left(U,b\right),$ (4)

然而，在出现间断解的情况下，上述恒等式(4)不再成立。特别地，如果对于所有给定的熵对 $\left(\eta ,q\right)$ ，能使

$\frac{\partial \eta \left(U\right)}{\partial t}+\frac{\partial q\left(U\right)}{\partial x}\le -V^{\text{T}}S\left(U,b\right),$ (5)

在分布意义上成立，则(1)的弱解 $U$ 被称为熵解(或物理相关解) [34] 。

通常，熵恒等式(4)和熵不等式(5)都被称为熵条件。理论上，熵条件对于双曲平衡定律(如(1))的适定性非常重要；数值上，满足(5)所述离散熵不等式的高阶格式可以提高其自身的鲁棒性。因此，对于给定的熵对，构造满足离散熵不等式(5)的数值格式具有重要意义。满足离散熵不等式的格式也称为熵稳定格式。近年来，熵稳定格式受到了广泛的关注。浅水波方程组的高阶熵稳定(ES)格式主要包括：间断Galerkin谱元法 [35] ，节点间断Galerkin方法 [36] ，间断Gallerkin (DG)方法 [37] [38] ，有限差分格式 [39] [40] 。

本文的主要目的是提出非平底地形上一维和二维SWE的高阶ES格式。通过对非平底地形进行适当的离散化，提出了两点熵守恒(EC)通量来构造满足熵恒等式的半离散二阶精确平衡EC格式。然后，以上述两点EC格式为构建块，相应地实现了高阶well-balanced的EC格式。然而，EC格式可能会在间断性周围产生数值振荡。为了克服这一缺点，考虑到尺度熵变量，使用WENO重构将一些合适的耗散项添加到EC通量中，以实现满足半离散熵不等式的高阶精确平衡ES格式。半离散EC格式和ES格式都与高阶Runge-Kutta方法相结合，用于时间离散。

本文的组织如下：在第2节中，我们给出了一维SWE的平衡EC有限体积格式。随后，我们在第3节中将所提出的格式扩展到二维领域。在第4节中，我们给出了大量的数值示例，以证明所得格式的性质。最后，第5节得出了相关结论。

2. 一维EC格式的构造

为方便起见，原始控制方程(1)可以改写如下

$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}=-G\frac{\partial b}{\partial x},$ (6)

其中 $G={\left(0,gh\right)}^{\text{T}}$ 。这里，对于一维SWE (6)，我们将总能量作为数学熵函数 [34] [40]

$\eta =\eta \left(U\right)=\frac{1}{2}\left(h{u}^2+g{h}^2\right)+ghb,$

以及由下式给出的熵变量

$V={\eta }^{\prime }{\left(U\right)}^{\text{T}}={\left(g\left(h+b\right)-\frac{1}{2}{u}^2,u\right)}^{\text{T}}.$

因此，考虑到一致条件(3)，熵通量为

$q=q\left(U\right)=\frac{1}{2}h{u}^3+g{h}^{2}u+ghub.$

此外，我们还定义了熵势 $\psi$ 如下

$\psi =V^{\text{T}}F\left(U\right)-q\left(U\right)=\frac{1}{2}g{h}^{2}u.$

2.1. EC格式

首先，我们将空间域划分为均匀网格，其中 $x_0<x_1<\cdots <x_N$ 作为网格节点以及空间步长取为 $\Delta x=x_i-x_{i-1}$ 。我们首先用数值通量 ${\tilde{F}}_{i+\frac{1}{2}}$ 构造EC格式。(6)的半离散有限体积格式具有以下形式

$\frac{\text{d}}{\text{d}t}{U}_i=-\frac{1}{\Delta x}\left({\tilde{F}}_{i+\frac{1}{2}}-{\tilde{F}}_{i-\frac{1}{2}}\right)+S_{\text{i}},$ (7)

$S_i$ 是源项的近似值，其中

$S_i=S_{i+\frac{1}{2}}+S_{i+\frac{1}{2}}:=-\frac{1}{2}{\left\{\left\{G\right\}\right\}}_{i+\frac{1}{2}}\frac{{\left[\left[b\right]\right]}_{i+\frac{1}{2}}}{△x}-\frac{1}{2}{\left\{\left\{G\right\}\right\}}_{i-\frac{1}{2}}\frac{{\left[\left[b\right]\right]}_{i-\frac{1}{2}}}{△x},$ (8)

其中 ${\left\{\left\{a\right\}\right\}}_{i+\frac{1}{2}}=\frac{1}{2}\left(a_i+a_{i+1}\right)$ 表示a在点 $x_{i+\frac{1}{2}}=x_i+\Delta x/2$ 的平均值。此外， ${\tilde{F}}_{i+\frac{1}{2}}$ 是在点 $x_{i+\frac{1}{2}}=x_i+\Delta x/2$ 处近似物理通量 $F\left(U\right)$ 的数值通量。

特别地，对于一些与给定熵通量q一致的数值熵通量 ${\tilde{q}}_{i+\frac{1}{2}}$ ，如果(7)的解满足以下半离散熵恒等式，则格式(7)被称为EC

$\frac{\text{d}}{\text{d}t}\eta \left(U_i\right)+\frac{1}{\Delta x}\left({\tilde{q}}_{i+\frac{1}{2}}-{\tilde{q}}_{i-\frac{1}{2}}\right)=0$ (9)

并且数值通量 ${\tilde{F}}_{i+\frac{1}{2}}$ 被相应地称为EC通量。事实上，格式(8)依赖于等式(4)表示离散状态。

下面的引理引入了给定的半离散格式(7)为EC的充分条件，以及源项的离散化。

引理1. 半离散格式(7)是EC并具有二阶精度，如果(7)中使用的一致两点数值通量 ${\tilde{F}}_{i+\frac{1}{2}}:=F\left(U_i,U_{i+1}\right)$ 满足

${\left[\left[V\right]\right]}^{\text{T}}\tilde{F}=\left[\left[\psi \right]\right]+g\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left[\left[b\right]\right],$ (10)

其中 $\left\{\left\{a\right\}\right\}=\frac{1}{2}\left(a_i+a_{i+1}\right),\left[\left[a\right]\right]=a_{i+1}-a_i$ 。(8)中的数值熵通量 ${\tilde{q}}_{i+\frac{1}{2}}$ 具有以下形式

${\tilde{q}}_{i+\frac{1}{2}}={\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}-{\left\{\left\{\psi \right\}\right\}}_{i+\frac{1}{2}}-\frac{1}{4}g{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[u\right]\right]}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}.$ (11)

证明：将(7)的两端左乘以 $V_i^{\text{T}}$ 得到

$\begin{array}{l}\frac{\text{d}{\eta }_i}{\text{d}t}=-\frac{1}{\Delta x}\left[V_i^{\text{T}}\left({\tilde{F}}_{i+\frac{1}{2}}-{\tilde{F}}_{i-\frac{1}{2}}\right)-\frac{1}{2}g{u}_i{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}\frac{{\left[\left[b\right]\right]}_{i+\frac{1}{2}}}{\Delta x}-\frac{1}{2}g{u}_i{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}\frac{{\left[\left[b\right]\right]}_{i-\frac{1}{2}}}{\Delta x}\right]\\ =-\frac{1}{\Delta x}\left[V_i^{\text{T}}\left({\tilde{F}}_{i+\frac{1}{2}}-{\tilde{F}}_{i-\frac{1}{2}}\right)+\frac{1}{2}g{u}_i{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}g{u}_i{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\right].\end{array}$

右侧可以进一步重新排列如下

$\begin{array}{l}{V}_i^{\text{T}}\left({\tilde{F}}_{i+\frac{1}{2}}-{\tilde{F}}_{i-\frac{1}{2}}\right)+\frac{1}{2}g{u}_i{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}g{u}_i{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\\ ={\left({\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}-\frac{1}{2}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}\right)}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}-{\left({\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}-\frac{1}{2}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}\right)}^{\text{T}}{\tilde{F}}_{i-\frac{1}{2}}\\ +\frac{1}{2}g\left({\left\{\left\{u\right\}\right\}}_{i+\frac{1}{2}}-{\left[\left[u\right]\right]}_{i+\frac{1}{2}}\right){\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}g\left({\left\{\left\{u\right\}\right\}}_{i-\frac{1}{2}}-{\left[\left[u\right]\right]}_{i-\frac{1}{2}}\right){\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\\ ={\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}-\frac{1}{2}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}+\frac{1}{2}g\left({\left\{\left\{u\right\}\right\}}_{i+\frac{1}{2}}-{\left[\left[u\right]\right]}_{i+\frac{1}{2}}\right){\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}\\ -{\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i-\frac{1}{2}}-\frac{1}{2}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i-\frac{1}{2}}+\frac{1}{2}g\left({\left\{\left\{u\right\}\right\}}_{i-\frac{1}{2}}-{\left[\left[u\right]\right]}_{i-\frac{1}{2}}\right){\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\end{array}$

$\begin{array}{l}={\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}-\frac{1}{2}{\left[\left[\psi \right]\right]}_{i+\frac{1}{2}}-\frac{1}{2}g{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left\{\left\{u\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}\\ +\frac{1}{2}g\left({\left\{\left\{u\right\}\right\}}_{i+\frac{1}{2}}-{\left[\left[u\right]\right]}_{i+\frac{1}{2}}\right){\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}\\ -{\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i-\frac{1}{2}}-\frac{1}{2}{\left[\left[\psi \right]\right]}_{i-\frac{1}{2}}-\frac{1}{2}g{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left\{\left\{u\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\\ +\frac{1}{2}g\left({\left\{\left\{u\right\}\right\}}_{i-\frac{1}{2}}-{\left[\left[u\right]\right]}_{i-\frac{1}{2}}\right){\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\end{array}$

$\begin{array}{l}={\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}-\frac{1}{2}{\left[\left[\psi \right]\right]}_{i+\frac{1}{2}}-\frac{1}{4}g{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left\{\left\{u\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}\\ -{\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i-\frac{1}{2}}-\frac{1}{2}{\left[\left[\psi \right]\right]}_{i-\frac{1}{2}}-\frac{1}{4}g{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left\{\left\{u\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\\ =\left({\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{\tilde{F}}_{i+\frac{1}{2}}-{\left[\left[\psi \right]\right]}_{i+\frac{1}{2}}-\frac{1}{4}g{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left\{\left\{u\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}\right)\\ -\left({\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^T{\tilde{F}}_{i-\frac{1}{2}}-{\left[\left[\psi \right]\right]}_{i-\frac{1}{2}}-\frac{1}{4}g{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left\{\left\{u\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\right)\\ ={\tilde{q}}_{i+\frac{1}{2}}-{\tilde{q}}_{i-\frac{1}{2}},\end{array}$

其中 $a_i={\left\{\left\{a\right\}\right\}}_{i+\frac{1}{2}}-\frac{1}{2}{\left[\left[a\right]\right]}_{i+\frac{1}{2}}$ 和 $a_i={\left\{\left\{a\right\}\right\}}_{i-\frac{1}{2}}-\frac{1}{2}{\left[\left[a\right]\right]}_{i-\frac{1}{2}}$ 用于第一个等式；条件(9)用于第二个等式； $\frac{1}{2}{\left[\left[a\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}{\left[\left[a\right]\right]}_{i-\frac{1}{2}}={\left\{\left\{a\right\}\right\}}_{i+\frac{1}{2}}-{\left\{\left\{a\right\}\right\}}_{i-\frac{1}{2}}$ 用于第三等式。因此，格式(7)是EC，并有以 ${\tilde{F}}_{i+\frac{1}{2}}={\tilde{F}}_{i+\frac{1}{2}}\left(U_i+U_{i+1}\right)$ 为数值通量形式。

此外，由于使用了二阶中心差分，源项的离散化也是二阶精确的。此外， [41] 中的结果表明，通量梯度的离散化也是二阶精度的。因此，格式(7)以及 ${\tilde{F}}_{i+\frac{1}{2}}={\tilde{F}}_{i+\frac{1}{2}}\left(U_i+U_{i+1}\right)$ 是二阶精确的。

接下来，我们给出满足条件(9)的EC通量。

定理2. 对于一维SWE(1)，以下数值通量

$\tilde{F}=\left(\begin{array}{c}\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\\ \left\{\left\{h\right\}\right\}{\left\{\left\{u\right\}\right\}}^2+\frac{g}{2}\left\{\left\{h^2\right\}\right\}\end{array}\right)$ (12)

是EC通量，并且与(1)中的物理通量 $F\left(U\right)$ 一致。

证明：这里，关键在于应用以下恒等式

$\left[\left[ab\right]\right]=\left\{\left\{a\right\}\right\}\left[\left[b\right]\right]+\left[\left[a\right]\right]\left\{\left\{b\right\}\right\},$

并用其等效形式重写下面的跳转

$\begin{array}{l}{\left[\left[V\right]\right]}^{\text{T}}={\left(g\left[\left[h\right]\right]+g\left[\left[b\right]\right]-\left\{\left\{u\right\}\right\}\left[\left[u\right]\right],\left[\left[u\right]\right]\right)}^{\text{T}},\\ \left[\left[\psi \right]\right]=g\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left[\left[h\right]\right]+\frac{1}{2}g\left\{\left\{h^2\right\}\right\}\left[\left[u\right]\right].\end{array}$

接下来，我们将上述表达式代入(9)，并得到

$\begin{array}{l}\left(g\left[\left[h\right]\right]+g\left[\left[b\right]\right]-\left\{\left\{u\right\}\right\}\left[\left[u\right]\right]\right)\cdot {\tilde{F}}^1+\left[\left[u\right]\right]\cdot {\tilde{F}}^2\\ =g\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left[\left[h\right]\right]+\frac{1}{2}g\left\{\left\{h^2\right\}\right\}\left[\left[u\right]\right]+g\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left[\left[b\right]\right]\\ -\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left\{\left\{u\right\}\right\}\left[\left[u\right]\right]+\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left\{\left\{u\right\}\right\}\left[\left[u\right]\right]\\ =\left(g\left[\left[h\right]\right]+g\left[\left[b\right]\right]-\left\{\left\{u\right\}\right\}\left[\left[u\right]\right]\right)\cdot \left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}+\left[\left[u\right]\right]\cdot \left(\left\{\left\{h\right\}\right\}{\left\{\left\{u\right\}\right\}}^2+\frac{1}{2}g\left\{\left\{h^2\right\}\right\}\right).\end{array}$

然后，通过比较相同跳跃的系数，我们得到以下等式

$\left(g\left[\left[h\right]\right]+g\left[\left[b\right]\right]-\left\{\left\{u\right\}\right\}\left[\left[u\right]\right]\right)\cdot {\tilde{F}}^1=\left(g\left[\left[h\right]\right]+g\left[\left[b\right]\right]-\left\{\left\{u\right\}\right\}\left[\left[u\right]\right]\right)\cdot \left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\},$

$\left[\left[u\right]\right]{\tilde{F}}^2=\left[\left[u\right]\right]\left(\left\{\left\{h\right\}\right\}{\left\{\left\{u\right\}\right\}}^2+\frac{1}{2}g\left\{\left\{h^2\right\}\right\}\right).$

并有

${\tilde{F}}^1=\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\},$

${\tilde{F}}^2=\left\{\left\{h\right\}\right\}{\left\{\left\{u\right\}\right\}}^2+\frac{1}{2}g\left\{\left\{h^2\right\}\right\}.$

最终获得了(11)中的数值通量。

备注1：此外，还容易验证EC通量(11)与物理通量 $F\left(U\right)$ 一致。实际上，在(11)中，设 $\left(h_L,u_L,b_L\right)=\left(h_R,u_R,b_R\right)=\left(h,u,b\right)$ 并有

${\tilde{F}}^1=hu,{\tilde{F}}^2=h{u}^2+\frac{1}{2}g{h}^2.$

2.2. Well-Balanced性质分析

在这里，我们简要总结了当前格式的well-balanced性质，并得出以下定理。

定理2. 具有EC通量(11)的当前格式(7)是well-balanced的，因为它保持了静水稳定状态(2)。换句话说，对于给定的初始条件 $u_i\equiv 0,h_i+b_i\equiv \text{Constant},\forall i$ ，(7)的解满足

$\frac{\text{d}}{\text{d}t}{h}_i\equiv 0,\frac{\text{d}}{\text{d}t}{\left(hu\right)}_i\equiv 0.$

证明：首先，在假设 $u_i\equiv 0,h_i+b_i\equiv \text{Constant}$ ，我们可以从(11)中发现 ${\tilde{F}}^1=0$ ，并根据(7)验证

$\frac{\text{d}}{\text{d}t}{h}_i\equiv 0.$

之后得到

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{\left(hu\right)}_i=-\frac{1}{\Delta x}\left({\tilde{F}}_{i+\frac{1}{2}}^2-{\tilde{F}}_{i-\frac{1}{2}}^2\right)-G_i\frac{{\left\{\left\{b\right\}\right\}}_{i+\frac{1}{2}}-{\left\{\left\{b\right\}\right\}}_{i-\frac{1}{2}}}{\Delta x}\\ =-\frac{g}{\Delta x}[\left(\frac{1}{2}{\left\{\left\{h^2\right\}\right\}}_{i+\frac{1}{2}}-\frac{1}{2}{\left\{\left\{h^2\right\}\right\}}_{i-\frac{1}{2}}\right)+\left({\left\{\left\{hb\right\}\right\}}_{i+\frac{1}{2}}-{\left\{\left\{hb\right\}\right\}}_{i-\frac{1}{2}}\right)\\ -\left({\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left\{\left\{b\right\}\right\}}_{i+\frac{1}{2}}-{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left\{\left\{b\right\}\right\}}_{i-\frac{1}{2}}\right)+h_i\left({\left\{\left\{b\right\}\right\}}_{i+\frac{1}{2}}-{\left\{\left\{b\right\}\right\}}_{i-\frac{1}{2}}\right)]\\ =-\frac{g}{\Delta x}[\left(\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[h\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[h\right]\right]}_{i-\frac{1}{2}}\right)+\left({\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{b}_{i+1}-{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{b}_{i-1}\right)\\ -\left({\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left\{\left\{b\right\}\right\}}_{i+\frac{1}{2}}-{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left\{\left\{b\right\}\right\}}_{i-\frac{1}{2}}\right)]\end{array}$

$\begin{array}{c}=-\frac{g}{\Delta x}\left[\left(\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[h\right]\right]}_{i-\frac{1}{2}}\right)+\frac{1}{2}\left({\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[b\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[b\right]\right]}_{i-\frac{1}{2}}\right)\right]\\ =-\frac{g}{2\Delta x}\left[\left({\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[h+b\right]\right]}_{i+\frac{1}{2}}+{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[h+b\right]\right]}_{i-\frac{1}{2}}\right)\right]\\ \equiv 0,\end{array}$

其中第二个等式中使用了 $\frac{1}{2}{\left\{\left\{h^2\right\}\right\}}_{i+\frac{1}{2}}-\frac{1}{2}{\left\{\left\{h^2\right\}\right\}}_{i-\frac{1}{2}}=\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i+\frac{1}{2}}{\left[\left[h\right]\right]}_{i+\frac{1}{2}}+\frac{1}{2}{\left\{\left\{h\right\}\right\}}_{i-\frac{1}{2}}{\left[\left[h\right]\right]}_{i-\frac{1}{2}}$ 。最后一个等式是由于 $\left[\left[h+b\right]\right]\equiv 0$ ，因为 $h_i+b_i=\text{Constant,}\forall \text{i}$ 。因此，该格式如所要求的那样是well-balanced的。

2.3. 熵稳定格式

我们知道，对于双曲守恒定律或平衡守恒定律，熵恒等式只能适用于光滑解。换句话说，由于激波等间断的存在，熵是不守恒的。更重要的是，熵守恒格式可能在间断处产生非物理震荡。这些现象促使我们提出满足给定熵对的离散熵不等式的熵稳定格式。特别是熵稳定格式可以抑制熵守恒格式在强间断处产生的数值振荡。根据 [41] 的想法，在此采用的主要措施是在熵守恒通量 ${\tilde{F}}_{i+\frac{1}{2}}$ 上增加一个耗散项，然后采用如下形式得到熵稳定通量

${\hat{F}}_{i+\frac{1}{2}}={\tilde{F}}_{i+\frac{1}{2}}-\frac{1}{2}{D}_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}},$ (13)

满足

${\left[\left[V\right]\right]}^{\text{T}}\hat{F}-\left[\left[\psi \right]\right]-g\left(\left[\left[hbu\right]\right]-\left[\left[hb\right]\right]\left\{\left\{b\right\}\right\}\right)\le 0,$ (14)

其中 $D$ 是一个半正定矩阵。

在此，我们采用

$D=R\left|\Lambda \right|R^{\text{T}},$ (15)

其中，R表示雅可比矩阵 $\frac{\partial F}{\partial U}$ 的右特征向量的矩阵，且当 ${\lambda }_1<{\lambda }_2<\cdots <{\lambda }_m$ 时， $\left|\Lambda \right|=diag\left(\left|{\lambda }_1\right|,\left|{\lambda }_2\right|,\cdots ,\left|{\lambda }_m\right|\right)$ 是的特征值。

结合(13)中的数值通量 $\hat{F}$ 和(8)中源项的离散化，得到最终格式

$\frac{\text{d}}{\text{d}t}{U}_i=-\frac{1}{\Delta x}\left({\hat{F}}_{i+\frac{1}{2}}-{\hat{F}}_{i-\frac{1}{2}}\right)+S_i.$ (16)

以下定理表明，高阶精确格式(16)满足熵离散不等式，因此是熵稳定的。

定理4. 高阶精确格式(16)是满足熵稳定的，即，其满足以下熵离散不等式

$\frac{\text{d}}{\text{d}t}\eta \left(U_i\right)+\frac{1}{\Delta x}\left({\hat{q}}_{i+\frac{1}{2}}-{\hat{q}}_{i-\frac{1}{2}}\right)\le 0,$

这里的数值熵通量是由

${\hat{q}}_{i+\frac{1}{2}}={\tilde{q}}_{i+\frac{1}{2}}-\frac{1}{2}{\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{D}_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}.$

证明：使用 $V_i^{\text{T}}$ 左乘(16)并运用(11)后，得到

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\eta \left(U_i\right)=-\frac{1}{\Delta x}\left({\tilde{q}}_{i+\frac{1}{2}}-{\tilde{q}}_{i-\frac{1}{2}}\right)+\frac{1}{2\Delta x}\left(V_i^{\text{T}}{D}_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}-V_i^{\text{T}}{D}_{i-\frac{1}{2}}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}\right)\\ =-\frac{1}{\Delta x}\left({\tilde{q}}_{i+\frac{1}{2}}-{\tilde{q}}_{i-\frac{1}{2}}\right)\\ +\frac{1}{2\Delta x}\left[\left({\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}-\frac{1}{2}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}^{\text{T}}\right)D_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}-\left({\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^{\text{T}}+\frac{1}{2}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}^{\text{T}}\right)D_{i-\frac{1}{2}}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}\right]\\ =-\frac{1}{\Delta x}\left[\left({\tilde{q}}_{i+\frac{1}{2}}-\frac{1}{2}{\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}{D}_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}\right)-\left({\tilde{q}}_{i-\frac{1}{2}}-\frac{1}{2}{\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^{\text{T}}{D}_{i-\frac{1}{2}}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}\right)\right]\\ -\frac{1}{4\Delta x}\left({\left[\left[V\right]\right]}_{i+\frac{1}{2}}^{\text{T}}{D}_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}+{\left[\left[V\right]\right]}_{i-\frac{1}{2}}^{\text{T}}{D}_{i-\frac{1}{2}}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}\right)\\ \text{:}=-\frac{1}{\Delta x}\left({\hat{q}}_{i+\frac{1}{2}}-{\hat{q}}_{i-\frac{1}{2}}\right)-\frac{1}{4\Delta x}\left({\left[\left[V\right]\right]}_{i+\frac{1}{2}}^{\text{T}}{D}_{i+\frac{1}{2}}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}+{\left[\left[V\right]\right]}_{i-\frac{1}{2}}^{\text{T}}{D}_{i-\frac{1}{2}}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}\right),\end{array}$ (17)

其中 $V_i^{\text{T}}={\left\{\left\{V\right\}\right\}}_{i+\frac{1}{2}}^{\text{T}}-\frac{1}{2}{\left[\left[V\right]\right]}_{i+\frac{1}{2}}^{\text{T}}$ 和 $V_i^{\text{T}}={\left\{\left\{V\right\}\right\}}_{i-\frac{1}{2}}^{\text{T}}+\frac{1}{2}{\left[\left[V\right]\right]}_{i-\frac{1}{2}}^{\text{T}}$ 已经在第二个等式中使用。由(15)中矩阵D的正定性，可得如下不等式

这里的数值熵通量是由

$\frac{\text{d}}{\text{d}t}\eta \left(U_i\right)+\frac{1}{\Delta x}\left({\hat{q}}_{i+\frac{1}{2}}-{\hat{q}}_{i-\frac{1}{2}}\right)\le 0.$

2.4. 时间离散

半离散格式(7)可以表示为

$\frac{\text{d}}{\text{d}t}{U}_i=F\left(U_i\right):=-\frac{1}{\Delta x}\left({\hat{F}}_{i+\frac{1}{2}}-{\hat{F}}_{i-\frac{1}{2}}\right)+S_i.$ (18)

另外，时间离散采用传统的Runge-Kutta方式，

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

2.5. 格式总结

这里的数值熵通量是由对于一维方程组(1)，所提格式(16)的具体步骤总结如下：

1) 首先，我们获取点 ${\left\{U_i\right\}}_{i+1}^N$ 的值。

2) 构造如(12)所示的高阶熵守恒通量。

3) 通过将耗散项加到原始熵守恒通量(12)上，构造高阶熵守恒通量(13)。

4) 通过公式(8)源项离散。

5) 得到半离散格式(18)。

6) 使用Runge-Kutta方法更新到下一个时间层(19)。

7) 重复步骤(2)-(6)。

3. 二维方程组

在本节中，我们将一维ES格式(16)拓展到以下二维SWES

$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}+\frac{\partial G\left(U\right)}{\partial y}=-S^1\frac{\partial b}{\partial x}-S^2\frac{\partial b}{\partial y},$ (20)

同时有

$U=\left(\begin{array}{c}h\\ hu\\ hv\end{array}\right),F\left(U\right)=\left(\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\\ huv\end{array}\right),G\left(U\right)=\left(\begin{array}{c}hu\\ huv\\ h{v}^2+\frac{1}{2}g{h}^2\end{array}\right),$

$S^1\left(U,b\right)=\left(\begin{array}{c}0\\ gh\\ 0\end{array}\right),S^2\left(U,b\right)=\left(\begin{array}{c}0\\ 0\\ gh\end{array}\right),$

其中v表示y方向上的速度。特别的，二位方程组(20)保持了下面静水稳态解

$u=v=0,h+b=\text{Constant}\text{.}$ (21)

这里我们要处理一个矩形空间域 $\Omega =\left[a,b\right]\times \left[c,d\right]$ 并在一个均匀的笛卡尔网络中 $\left(x_i,y_j\right),i=0,1,\cdot \cdot \cdot ,N_x;j=0,1,\cdot \cdot \cdot ,N_y$ 其空间步长为 $\Delta x=x_i-x_{i-1},\Delta y=y_i-y_{i-1}$ ，这相当于用 $a=x_0<x_1<\cdots <x_{N_x}=b$ 并且 $c=y_0<y_1<\cdots <y_{N_y}=d$ 。

与一维情况类似，我们取下列熵函数

$\eta \left(U,b\right)=\frac{1}{2}\left(h{u}^2+h{v}^2+g{h}^2\right)+ghb,$

还有下面的熵变量

$V=\eta \text{'}{\left(U\right)}^{\text{T}}={\left(g\left(h+b\right)-\frac{1}{2}{u}^2-\frac{1}{2}{v}^2,u,v\right)}^{\text{T}}.$

相应的在x-，y-方向上的熵变量定义为

$\begin{array}{l}{q}^x\left(U\right)=\frac{1}{2}h{u}^3+g{h}^{2}u+ghub+\frac{1}{2}\left(hu\right)v^2,\\ q^y\left(U\right)=\frac{1}{2}h{v}^3+g{h}^{2}v+ghvb+\frac{1}{2}\left(hv\right)u^2.\end{array}$ (22)

与一维情况类似，在x-和y-方向上的二阶EC通量为

$\begin{array}{l}\tilde{F}=\tilde{F}\left(U\right)=\left(\begin{array}{c}\left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\\ \left\{\left\{h\right\}\right\}{\left\{\left\{u\right\}\right\}}^2+\frac{g}{2}\left\{\left\{h^2\right\}\right\}\\ \left\{\left\{h\right\}\right\}\left\{\left\{u\right\}\right\}\left\{\left\{v\right\}\right\}\end{array}\right),\\ \tilde{G}=\tilde{G}\left(U\right)=\left(\begin{array}{c}\left\{\left\{h\right\}\right\}\left\{\left\{v\right\}\right\}\\ \left\{\left\{h\right\}\right\}\left\{\left\{v\right\}\right\}\left\{\left\{u\right\}\right\}\\ \left\{\left\{h\right\}\right\}{\left\{\left\{v\right\}\right\}}^2+\frac{g}{2}\left\{\left\{h^2\right\}\right\}\end{array}\right).\end{array}$

然后再EC通量上加上适当的耗散项产生高阶ES通量，如下所示，

$\begin{array}{l}{\hat{F}}_{i+\frac{1}{2},j}^{kth}={\hat{F}}_{i+\frac{1}{2},j}^{2pth}-\frac{1}{2}{\alpha }_{i+\frac{1}{2},j}{S}_{i+\frac{1}{2},j}{R}_{i+\frac{1}{2},j}{\langle \langle w\rangle \rangle }_{i+\frac{1}{2},j},\\ {\hat{G}}_{i,j+\frac{1}{2}}^{kth}={\hat{G}}_{i,j+\frac{1}{2}}^{2pth}-\frac{1}{2}{\alpha }_{i,j+\frac{1}{2}}{S}_{i,j+\frac{1}{2}}{R}_{i,j+\frac{1}{2}}{\langle \langle w\rangle \rangle }_{i,j+\frac{1}{2}},\end{array}$ (23)

当k为偶数时 $p=k/2$ ，当k为奇数时 $p=\left(k+1\right)/2$ 。这里跳跃 ${\langle \langle w\rangle \rangle }_{i+\frac{1}{2},j},{\langle \langle w\rangle \rangle }_{i,j+\frac{1}{2}}$ 分别在x-和y-方向上通过WENO重构得到，并且黏度 ${\alpha }_{i+\frac{1}{2},j}$ 和 ${\alpha }_{i,j+\frac{1}{2}}$ 分别在x-和y-方向上选取。

此外，对于源项，我们首先对原想中的倒数进行近似，如下所示，

$\begin{array}{l}\frac{1}{\Delta x}\left({\tilde{b}}_{i+\frac{1}{2},j}^{2pth}-{\tilde{b}}_{i-\frac{1}{2},j}^{2pth}\right)={\frac{\partial b}{\partial x}{\rm B}|}_{i,j}+O\left(\Delta x^{2p}\right),\\ \frac{1}{\Delta y}\left({\tilde{b}}_{i,j+\frac{1}{2}}^{2pth}-{\tilde{b}}_{i,j-\frac{1}{2}}^{2pth}\right)={\frac{\partial b}{\partial y}{\rm B}|}_{i,j}+O\left(\Delta y^{2p}\right),\end{array}$

并且

$\begin{array}{l}{\tilde{b}}_{i+\frac{1}{2},j}^{2pth}={\displaystyle \underset{r=1}{\overset{p}{\sum }}{\alpha }_r^p{\displaystyle \underset{s=0}{\overset{p}{\sum }}\left(b_{i-s,j}+b_{i-s+r,j}\right)}},\\ {\tilde{b}}_{i,j+\frac{1}{2}}^{2pth}={\displaystyle \underset{r=1}{\overset{p}{\sum }}{\alpha }_r^p{\displaystyle \underset{s=0}{\overset{p}{\sum }}\left(b_{i,j-s}+b_{i,j-s+r}\right)}}.\end{array}$

然后我们用下面的形式近似源项，

$\begin{array}{l}-S^1{\frac{\partial b}{\partial x}|}_{i,j}\approx -\frac{S_{ij}^1}{\Delta x}\left({\tilde{b}}_{i+\frac{1}{2},j}^{2pth}-{\tilde{b}}_{i-\frac{1}{2},j}^{2pth}\right),\\ -S^2{\frac{\partial b}{\partial x}|}_{i,j}\approx -\frac{S_{ij}^2}{\Delta x}\left({\tilde{b}}_{i,j+\frac{1}{2}}^{2pth}-{\tilde{b}}_{i,j-\frac{1}{2}}^{2pth}\right).\end{array}$ (24)

最后我们得到(20)的半离散ES有限体积格式，格式如下，

$\begin{array}{l}\frac{\text{d}}{\text{d}t}{U}_{ij}=-\frac{1}{\Delta x}\left({\hat{F}}_{i+\frac{1}{2},j}^{kth}-{\hat{F}}_{i-\frac{1}{2},j}^{kth}\right)-\frac{1}{\Delta y}\left({\hat{G}}_{i,j+\frac{1}{2}}^{kth}-{\hat{G}}_{i,j-\frac{1}{2}}^{kth}\right)\\ -\frac{S_{ij}^1}{\Delta x}\left({\tilde{b}}_{i+\frac{1}{2},j}^{2pth}-{\tilde{b}}_{i-\frac{1}{2},j}^{2pth}\right)-\frac{S_{ij}^2}{\Delta y}\left({\tilde{b}}_{i,j+\frac{1}{2}}^{2pth}-{\tilde{b}}_{i,j-\frac{1}{2}}^{2pth}\right),\end{array}$ (25)

结合(23)中的ES通量 ${\hat{F}}_{i+\frac{1}{2},j}^{kth}$ 和 ${\hat{G}}_{i,j+\frac{1}{2}}^{kth}$ 以及(24)中对源项的离散化。

此外，(25)对时间的离散化也使用了龙格-库塔方法(19)。

4. 数值结果

在本节中，我们将通过计算不同算例来测试当前格式的性能。在所有的计算中我们采用五阶WENO重建(即 $k=5$ )，采用 $p=3$ 。因此我们得到了五阶精度的ES通量。总之，我们得到了具有五阶精度的空间离散化的半离散格式。此外，为了保证数值的稳定性，我们取Courant-Friedrichs-Levy (CFL)数值为0.6。

4.1. Well-Balanced性质的测试

首先我们应用一个例子 [9] [42] [43] [44] 和下面初始条件，

$u=0\text{m}/\text{s}\text{and}h+b=10\text{m},x\in \left[\text{0,10}\right]\text{m}$

数值证明了良好的性质。我们考虑了两个不同的底部：第一个是光滑的

$b\left(x\right)=5\exp \left(-\frac{2}{5}{\left(x-5\right)}^2\right)\text{m},$

第二个底部是间断的

$b\left(x\right)=\{\begin{array}{l}4\text{mif}4\le x\le 8\text{m},\\ 0\text{motherwise}\text{.}\end{array}$

表1和表2分别给出了解在 $t=0\text{s}$ 和 $t=0.5\text{s}$ 之间的数值误差。很明显所有的误差都在机器精度的同一量级，因此实现了预期良好平衡的性能。

Table 1. Numerical errors of the example in Section 4.1 over a smooth bottom

表1. 算例4.1节在光滑底部上的数值误差

Table 2. Numerical errors of the example in section 4.1 over a discontinuous bottom

表2. 算例4.1在间断底部上的数值误差

4.2. 精度测试

接下来，我们使用 [9] [42] [43] [44] 中的一个例子，以及下面的底部地形和下面的初始条件

$\begin{array}{l}b\left(x\right)={\sin }^2\left(\pi x\right)\text{m},\\ h\left(x,0\right)=\left(5+\exp \left(\cos \left(2\pi x\right)\right)\right)\text{m},x\in \left[0,1\right]\text{m}\\ \left(hu\right)\left(x,0\right)=\sin \left(\cos \left(2\pi x\right)\right){\text{m}}^2/\text{s},\end{array}$

来测试精度的顺序。

由于解析解无法得到，我们首先需要在一个12,800单元的精细网格上计算这个例子直到 $t=0.1\text{s}$ ， 然后将得到的数值解作为参考。随后，我们计算 $t=0.1\text{s}$ 时的 $L^1$ 误差，并根据参考解得到精度阶数，如表3所示。

Table 3. L1 error and orders of accuracy for the example in Section 4.2

表3. 算例4.2的L¹误差和精度阶数

说明预期的五阶精度明显达到。

4.3. 定常流的扰动问题

在这里，本例用于数值测试当前格式根据给定稳态 [45] - [50] 在凹凸形状底部地形上捕捉小扰动的能力。

$b\left(x\right)=\{\begin{array}{l}0.25\left(\cos \left(10\pi \left(x-1.5\right)\right)+1\right)\text{mif}\text{1}\text{.4}\le x\le 1.6\text{m},\\ 0\text{motherwise},\end{array}x\in \left[0,2\right]\text{m}.$

认为初始条件是对静水稳态的一个小扰动，即

$h\left(x,0\right)=\{\begin{array}{l}\left(1-b\left(x\right)+\epsilon \right)\text{mif}\text{1}\text{.1}\le x\le 1.2\text{m},\\ \left(1-b\left(x\right)\right)\text{motherwise}，\end{array}\text{and}u\left(x,0\right)=0\text{m}/\text{s}\text{,}$

其中 $\epsilon >0$ 。在这里，我们考虑两种不同情况：对于大脉冲 $\epsilon =0.2\text{m}$ ；对于单独的小脉冲 $\epsilon =0.001\text{m}$ 。

随着时间的推移，初始扰动分解为两个向不同方向运动的脉冲。图1和图2给出了 $t=0.2\text{s}$ 时，在200哥单元网格上两种情况的数值解。这两个脉冲即使是在相对粗糙的网格上都得到了很好的解析，并且与参考文献 [9] [50] [51] 中的结果一致。

4.4. 矩形凸起上的溃坝问题

然后，我们来数值模拟一个矩形底部上的溃坝问题 [7] [9] [14] [51] 。其中，初始条件为

$\begin{array}{l}h\left(x,0\right)=\{\begin{array}{l}\left(20-b\left(x\right)\right)\text{mif}x\le 750\text{m},\\ \left(15-b\left(x\right)\right)\text{motherwise},\end{array}x\in \left[0,1500\right]\text{m},\\ u\left(x,0\right)=0\text{m}/\text{s},\end{array}$

Figure 1. Surface level h + b (left) and water discharge hu (right) of the examples in section 4.3 at t = 0.2 s

图1. 在t = 0.2 s时，4.3节算例的表面水位h + b (左)和排水(右)

Figure 2. Free surface level h + b of the example in section 4.4 at t = 15 s (left) and t = 60 s (right)

图2. 在t = 15 s (左)和t = 60 s (右)时，在4.4节算例中的自由表面能级h + b

底部地形呈矩形

图2分别给出了 $t=15\text{s}$ 和 $t=60\text{s}$ 时的数值解。显然目前格式得到了很好的数值解，与参考值和 [7] [9] [14] [51] 中的数值解一致。

4.5. 驼峰底部上的定常流

进一步，通过一个广泛应用的实例 [52] 对所提出的格式进行验证。实际上，本例基于给出的初始条件对跨临界流和次临界流进行数值模拟

$h\left(x,0\right)=0.33\text{mand}u\left(x,0\right)=0\text{m}/\text{s},x\in \left[0,25\right]\text{m}$

底部地形呈驼峰状

$b\left(x\right)=\{\begin{array}{l}\left(0.2-0.05{\left(x-10\right)}^2\right)\text{mif}\text{8}\le x\le \text{12m,}\\ 0\text{motherwise}\text{.}\end{array}$

随后，我们在具有200个单元的网格上计算该例子，并计算至 $t=200\text{s}$ ，施加不同的边界条件。此外，我们还展示了从 [53] 中得到的精确解，以便更好地进行比较。

l 算例A：无激波的跨临界流动

上游边界的流量 $hu=1.53{\text{m}}^{\text{2}}/\text{s}$ ，下游边界的水深 $h=0.66\text{m}$ 。图3中算例A结果表明，数值计算结果能够和实际结果吻合。

Figure 3. Surface level h + b in section 4.5 at t = 200

图3. 在t = 200时，第4.5节中的表面水位h + b

l 算例B：带激波的跨临界流动

上游边界的排水量 $hu=0.18{\text{m}}^{\text{2}}/\text{s}$ ，下游边界的水深 $h=0.33\text{m}$ 。图3中算例B结果表明，数值解与精确解具有良好的一致性，且无伪振荡。

l 算例C：亚临界流

其中，上游边界流量 $hu=0.18{\text{m}}^{\text{2}}/\text{s}$ ，下游边界水深 $h=2\text{m}$ 。图3中算例C表明，数值解与精确解拟合效果很好。

4.6. 台阶状底部上的溃坝问题

l 算例A：

此外，我们实现了一个来自 [54] [55] 的示例以及下面的初始条件

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(4,0\right)\text{if}x\le \text{0m,}\\ \left(1,0\right)\text{otherwise,}\end{array}x\in \left[-10,10\right]\text{m}$

并超过一个台阶形状的底部

$b\left(x\right)=\{\begin{array}{l}\text{0}\text{mif}x\le \text{0m},\\ \text{1}\text{motherwise}\text{.}\end{array}$

随着时间的推移，这个数值算例产生了一个向左移动的稀疏薄和一个向右移动的激波。

l 算例B：

此外，我们使用不同的初始条件

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(4,5\right)\text{mif}x\le \text{0m,}\\ \left(1,-0.9\right)\text{motherwise},\end{array}x\in \left[-10,10\right]\text{m}$

与上面相同的底部。随着时间推移，这个算例出现了两个向不同方向移动的冲击。

l 算例C：

初始条件由

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(0.75,-9.49365\right)\text{if}x\le \text{0m,}\\ \left(1.10594,-4.94074\right)\text{otherwise},\end{array}x\in \left[-15,5\right]\text{m}$

并超过一个台阶形状的底部

$b\left(x\right)=\{\begin{array}{l}0\text{mif}x\le \text{0m,}\\ \text{0}\text{.2}\text{motherwise,}\end{array}$ (26)

l 算例D：

初始条件如下

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(0.75,-1.35624\right)\text{if}x\le \text{0m,}\\ \left(1.10594,-4.94074\right)\text{otherwise},\end{array}x\in \left[-10,4\right]\text{m}$

水体底部为(26)式。

图4给出了上述情况在同一网格上400个单元与精确网格的数值解，很明显，两者吻合较好。

Figure 4. Surface level h + b of the examples in Section 4.6

图4. 4.6节算例的表面h + b层

5. 结论

本文提出了一、二维非平坦地形浅水波方程的高阶精确熵稳定有限体积格式。首先，我们构造了二阶精确的well-balanced的半离散熵守恒有限体积格式，同时匹配物理通量的离散化和与非平底地形有关的源项的离散化。其次，以二阶精确熵守恒为基础，实现高阶精确均衡熵守恒格式。但熵守恒格式可能在间断附近产生杂散振荡，因此我们在熵守恒通量中加入适当的耗散项并构造满足半离散熵不等式的高阶半离散良好平衡的熵守恒格式。为了构造高阶耗散项，本文采用了基于缩放的熵变量的WENO重构。然后，采用高阶显式Runge-Kutta方法实现了时间区域的离散化。最后，大量的算例表明，本文所提出的格式保持了高阶精度，具有良好的平衡性，并能很好地捕捉到稳态下的小扰动。

致谢

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

NOTES

^*通讯作者。