求解带源项的浅水波方程的高分辨率熵相容格式

High Resolution Entropy Consistent Scheme for Solving Shallow Water Wave Equations with Source Terms

摘要: 浅水波方程对湖泊、河流等波动问题的研究具有重要意义。源项是对底部地势的描述。带源项的浅水波方程可以归结为非线性的双曲守恒律问题。本文采用结构网格，构造了一种熵相容格式求解带源项的浅水波方程，并对熵守恒变量使用基于MUSCL格式的斜率限制器重构，构造具有2阶精度的熵相容格式。在数值实验中证明了该格式有效地避免了非物理现象的产生，并且可以准确地捕捉激波，具有良好的稳健性。

Abstract: Shallow water wave equation is of great significance to the study of wave problems in lakes and riv-ers. The source term is a description of the topography at the bottom. The shallow water wave equation with source term can be reduced to a nonlinear hyperbolic conservation law problem. In this paper, an entropy consistent scheme is constructed to solve the shallow water wave equation with the source term by using the structural grid. The entropy conservation variables are recon-structed by using a slope limiter which is based on the MUSCL scheme to construct an entropy con-sistent scheme with second-order accuracy. Numerical experiments show that the scheme can ef-fectively avoid the occurrence of non-physical phenomena, and can accurately capture shock waves, which has good robustness.

1. 引言

浅水波方程是对各种环境流体运动的数学描述，源项是底部地形。它对海洋，河流等波动问题的研究具有重要意义。带有源项的浅水波方程具有如下形式：

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

其中

$U=\left(\begin{array}{c}h\\ hu\end{array}\right),f\left(U\right)=\left(\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\end{array}\right),S=\left(\begin{array}{c}0\\ ghb\end{array}\right),$

h表示水深，u表示水流速度，g表示重力加速度， $b=b\left(x\right)$ 表示底部地势，是关于x的函数。

当源项为0时，浅水波方程可以归结为双曲守恒律方程。即便双曲守恒律方程的初值函数是光滑的，它的解也可能会产生激波、接触间断等，这种情况的出现是不符合古典解有关理论的。为解决该类问题，1954年，Lax提出了弱解 [1] 这一概念，为数值方法的发展打下了良好的理论基础。Lax在1973年证明了熵稳定条件等价于“粘性消失” [2]。文献 [3] 把熵稳定条件引入到了带源项的浅水波方程。1987年Tadmor提出了熵变量和熵势的概念，构造了二阶精度的熵守恒格式，可以适用于各种网格 [4]。在熵守恒格式的基础上适当地加入一些粘性，就可以达到熵稳定。2006年，Roe提出在熵守恒格式上增加一个数值粘性项，称为ERoe格式，这是一类熵稳定格式 [5]，该格式保持了Roe格式中，对激波的捕捉效果较好这一特点。熵守恒和熵稳定格式是获得具有物理意义的，且数值上稳定的解的一种高效途径，具有更实用的指导意义。2009年，Roe和Ismail进一步分析得到：当解在跨过激波时，会产生激波强度的立方量级熵增，在这基础上继续发展了求解Burgers方程与Euler方程时，将数值粘性项发展为更精准量化的熵相容格式 [6]。熵相容格式对数值粘性的控制比熵稳定格式的数值粘性更为精细，其数值结果更加符合物理现象，保证了激波稳定性的同时能得到更加精确的数值解。2018年，张海军，封建湖等构造了带源项的浅水波方程的熵稳定格式，模拟效果较好 [7]。2021年，任璇，封建湖等，提出了一种基于MUSCL重构的斜率限制器，有效地提高了双曲方程熵相容格式的分辨率 [8]。

本文采用有限体积法对带源项的浅水波方程求解，对源项使用中心差分格式离散，在熵相容格式的基础上，使用基于MUSCL重构的斜率限制器，对方程求解。通过数值实验证明该格式具有良好的性质，有效地避免了产生伪振荡等非物理现象。

2. 有限体积法

对式(1)在 $\left[x_{i-1/2},x_{i+1/2}\right]$ 上积分，得到半离散格式

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

其中 $U_i$ 是 $\left[x_{i-1/2},x_{i+1/2}\right]$ 上的单元平均， $F$ 是数值通量， ${\tilde{S}}_i=\left(\begin{array}{c}0\\ \frac{g}{\Delta x}\left({\bar{h}}_{i+\frac{1}{2}}{\left[b\right]}_{i+\frac{1}{2}}+{\bar{h}}_{i-\frac{1}{2}}{\left[b\right]}_{i-\frac{1}{2}}\right)\end{array}\right)$ 是地势的中心差分。

定义

${\bar{\phi }}_{i+\frac{1}{2}}=\frac{{\phi }_{i+1}+{\phi }_i}{2},{\left[\phi \right]}_{i+\frac{1}{2}}={\phi }_{i+1}-{\phi }_i.$

时间上采用使用三阶Runge-Kutta时间离散

$\{\begin{array}{l}{U}_i^{*}=U_i^n+\Delta tL\left(U_i^n\right),\\ U_i^{**}=\frac{3}{4}{U}_i^n+\frac{1}{4}{U}_i^{*}+\frac{1}{4}\Delta tL\left(U_i^{*}\right),\\ U_i^{\left(n+1\right)}=\frac{1}{3}{U}_i^n+\frac{2}{3}{U}_i^{**}+\frac{2}{3}\Delta tL\left(U_i^{**}\right).\end{array}$ (3)

其中 $L\left(U\right)=-\frac{1}{\Delta x}\left(F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}}\right)-{\tilde{S}}_i$。

3. 带源项浅水波方程的求解

3.1. 熵相容格式

令熵函数为 $E\left(U\right)=\frac{\left(h{u}^2+g{h}^2\right)}{2}+ghb$，熵通量函数 $H\left(U\right)=gu{h}^2+\frac{h{u}^3}{2}$，定义熵变量 $V=E^{\prime }\left(U\right)={\left(g\left(h+b\right)-\frac{u^2}{2},u\right)}^{\text{T}}$，熵势 $\Phi =\frac{1}{2}gu{h}^2$，从而可得熵守恒格式的数值通量

$F_{i+\frac{1}{2}}^c=\left(\begin{array}{c}{\bar{h}}_{i+\frac{1}{2}}{\bar{u}}_{i+\frac{1}{2}}\\ \frac{g}{2}{\left({\bar{h}}^2\right)}_{i+\frac{1}{2}}+{\bar{h}}_{i+\frac{1}{2}}{\left({\bar{u}}_{i+\frac{1}{2}}\right)}^2\end{array}\right).$ (4)

双曲守恒律方程的解在光滑处才是守恒的，根据熵条件

${\partial }_{t}E\left(U\right)+{\partial }_{x}H\left(U\right)\le 0,$ (5)

在激波等不连续点附近需要产生耗散，因此，需要额外添加耗散项，构建熵稳定格式。将熵变量的数值耗散项添加到熵守恒格式的数值通量，有

$F_{i+\frac{1}{2}}^{es}=F_{i+\frac{1}{2}}^c-\frac{1}{2}{D}_{i+\frac{1}{2}}\left(V_{i+1}-V_i\right),$ (6)

其中 $D_{i+\frac{1}{2}}$ 是对称半正定矩阵。

Roe格式的数值通量具有如下形式 [6]

$F_{i+\frac{1}{2}}^{Roe}=\frac{1}{2}\left(f\left(U_{i+1}\right)+f\left(U_i\right)\right)-\frac{1}{2}{R}_{i+\frac{1}{2}}{\Lambda }_{i+\frac{1}{2}}{R}_{i+\frac{1}{2}}^{-1}\left(U_{i+1}-U_i\right),$ (7)

其中 $R_{i+\frac{1}{2}}$ 是通量雅可比矩阵 ${\partial }_{U}f\left(U\right)$ 的特征向量， ${\Lambda }_{i+\frac{1}{2}}=diag\left(\left|{\lambda }_1\right|,\left|{\lambda }_2\right|,\cdots ,\left|{\lambda }_n\right|\right)$ 是以通量雅可比矩阵的特征值为对角线元素的非负对角矩阵。这些矩阵和特征值都是Roe平均状态。

因为 $\Delta U={\partial }_{V}U\Delta V$，其中 ${\partial }_{V}U$ 是对称正定的。由Barth特征向量缩放定理 [9]，有 ${\partial }_{V}U=\tilde{R}{\tilde{R}}^{\text{T}}$，Roe通量改写为

$F_{i+\frac{1}{2}}^{Roe}=\frac{1}{2}\left(f\left(U_{i+1}\right)+f\left(U_i\right)\right)-\frac{1}{2}{\tilde{R}}_{i+\frac{1}{2}}{\Lambda }_{i+\frac{1}{2}}{\tilde{R}}_{i+\frac{1}{2}}^{\text{T}}{\left[V\right]}_{i+\frac{1}{2}}.$ (8)

于是使用Roe的数值耗散项

$D_{i+\frac{1}{2}}={\tilde{R}}_{i+\frac{1}{2}}{\Lambda }_{i+\frac{1}{2}}{\tilde{R}}_{i+\frac{1}{2}}^{\text{T}},$ (9)

为了简便，下文将 ${\tilde{R}}_{i+\frac{1}{2}}$ 写为 $R_{i+\frac{1}{2}}$。

带源项的浅水波方程有

$R_{i+\frac{1}{2}}=\frac{1}{\sqrt{2g}}\left(\begin{array}{cc}1& 1\\ {\bar{u}}_{i+\frac{1}{2}}-\sqrt{g{\bar{h}}_{i+\frac{1}{2}}}& {\bar{u}}_{i+\frac{1}{2}}+\sqrt{g{\bar{h}}_{i+\frac{1}{2}}}\end{array}\right),$ (10)

${\Lambda }_{i+\frac{1}{2}}=diag\left(\left|{\bar{u}}_{i+\frac{1}{2}}-\sqrt{g{\bar{h}}_{i+\frac{1}{2}}}\right|,\left|{\bar{u}}_{i+\frac{1}{2}}+\sqrt{g{\bar{h}}_{i+\frac{1}{2}}}\right|\right).$ (11)

Roe和Ismail指出，解在跨越激波时会产生立方级的熵增 [6]，因此，需要更多的耗散项修正熵守恒格式，从而抑制振荡。

在二维双曲守恒律方程中，跨越激波产生熵增为 [6]

${\tilde{U}}_{i+\frac{1}{2}}=\frac{1}{2}{\left[V\right]}_{i+\frac{1}{2}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{-1}\left(U\left(\xi \right)\right)\text{d}\xi {\left[U\right]}_{i+\frac{1}{2}}}.$ (12)

选择积分路径为 $U=\left(1-\xi \right)U_i+\xi U_{i+1}$， $0<\xi <1$。熵增总是正的，即有 $\tilde{U}>0$。

熵增 $\tilde{U}$ 有

$\begin{array}{c}{\tilde{U}}_{i+\frac{1}{2}}=\frac{1}{2}{\left[V\right]}_{i+\frac{1}{2}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{-1}\left(U\left(\xi \right)\right)\text{d}\xi {\left[U\right]}_{i+\frac{1}{2}}}\\ =\frac{1}{2}{\left[V\right]}_{i+\frac{1}{2}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{-1}\left(U\left(\xi \right)\right){\left[U\right]}_{i+\frac{1}{2}}\text{d}\xi }\\ \approx \frac{1}{2}{\left[V\right]}_{i+\frac{1}{2}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{-1}\left(U\left(\xi \right)\right)R\left(U\left(\xi \right)\right)R^{\text{T}}\left(U\left(\xi \right)\right){\left[V\right]}_{i+\frac{1}{2}}\text{d}\xi }\\ =\frac{1}{2}{\left[V\right]}_{i+\frac{1}{2}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{\text{T}}\left(U\left(\xi \right)\right)\text{d}\xi {\left[V\right]}_{i+\frac{1}{2}}}.\end{array}$ (13)

为了保证(13)的符号性，只需保证积分项正定，因此，作如下修正：

${\tilde{U}}_{i+\frac{1}{2}}=\frac{1}{2}{\left[V\right]}_{i+\frac{1}{2}}\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{\text{T}}\left(U\left(\xi \right)\right)\text{d}\xi }\right|{\left[V\right]}_{i+\frac{1}{2}},$ (14)

对应耗散项为

$D_{i+\frac{1}{2}}^{ec}=\frac{1}{2}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\left(U\left(\xi \right)\right)\Lambda \left(U\left(\xi \right)\right)R^{\text{T}}\left(U\left(\xi \right)\right)\text{d}\xi {\left[V\right]}_{i+\frac{1}{2}}},$ (15)

对上式的积分项使用Gauss-Legendre求积公式近似，于是有

$D_{i+\frac{1}{2}}^{ec}=\frac{1}{2}{R}_{i+\frac{1}{2}}\left|{\displaystyle \underset{k=1}{\sum }{\omega }_k\left(2s_k-1\right)\Lambda \left(U\left(s_k\right)\right)}\right|R_{i+\frac{1}{2}}^{\text{T}}{\left[V\right]}_{i+\frac{1}{2}},$ (16)

其中 $\omega$ 为求积系数， $s_k$ 为高斯点，在 $\left[0,1\right]$ 上积分。本文中采用三点求积公式，高斯点和求积系数分别为 $s_1=-\sqrt{\frac{3}{5}},s_2=0,s_3=\sqrt{\frac{3}{5}}$， ${\omega }_1=\frac{5}{9},{\omega }_2=\frac{8}{9},{\omega }_3=\frac{5}{9}$。因此，熵相容格式的数值通量为

$F_{i+\frac{1}{2}}^{ec}=F_{i+\frac{1}{2}}^c-\frac{1}{2}{R}_{i+\frac{1}{2}}\left({\Lambda }_{i+\frac{1}{2}}+\left|{\displaystyle \underset{k=1}{\sum }{\omega }_k\left(2s_k-1\right)\Lambda \left(U\left(s_k\right)\right)}\right|\right)R_{i+\frac{1}{2}}^{\text{T}}{\left[V\right]}_{i+\frac{1}{2}}.$ (17)

3.2. 基于MUSCL格式的斜率限制器

在每个控制体单元 $I_i=\left[x_{i-1/2},x_{i+1/2}\right]$ 上，单元长度为 $\Delta x$，对每个 $U_i$ 进行重构，记重构后的函数为 $L_i\left(x\right)$，如果有

$\left|L_i\left(x\right)-U_i\left(x\right)\right|\le O\left({\left(\Delta x\right)}^{k+1}\right),$ (18)

则称这个重构函数具有k-阶精度。 $U_i^n$ 表示在第n个时间层上，控制体单元 $I_i$ 的平均，有

$U_i^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{L}_i\left(x\right)\text{d}x}=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{U}_i\left(x\right)\text{d}x},$ (19)

于是有，

$L_i\left(x\right)=U_i^n+\left(x-x_i\right)\frac{{\Delta }_i}{\Delta x},x\in \left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]，$ (20)

其中 ${\Delta }_i$ 是选择的重构斜率。于是，点 $x_i$ 处的左右状态分别为

$\{\begin{array}{l}{L}_i^l=U_i^n-\frac{1}{2}{\Delta }_i,\\ L_i^r=U_i^n+\frac{1}{2}{\Delta }_i.\end{array}$ (21)

构造时，首先将重构函数 $L_i\left(x\right)$ 进行泰勒展开，有

$L_i\left(x\right)=U_i\left(x_i\right)+{\frac{\partial U_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}\frac{{\left(x-x_i\right)}^2}{2}+\cdots$ (22)

继续将 $U_i\left(x\right)$ 在 $x_i$ 处进行泰勒展开，可得

$\begin{array}{c}{U}_i^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{U}_i\left(x\right)\text{d}x}\\ =\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}\left(U_i\left(x_i\right)+{\frac{\partial U_i}{\partial x}|}_{x_i}\left(x-x_i\right)+{\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}\frac{{\left(x-x_i\right)}^2}{2}+\cdots \right)\text{d}x}\\ =U_i\left(x_i\right)+0+\frac{1}{\Delta x}\cdot {\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}\cdot {\displaystyle {\int }_{I_i}\frac{{\left(x-x_i\right)}^2}{2}\text{d}x+\cdots }\end{array}$ (23)

于是

$U_i\left(x_i\right)=U_i^n-\frac{1}{\Delta x}{\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}{\displaystyle {\int }_{I_i}\frac{{\left(x-x_i\right)}^2}{2}\text{d}x}+\cdots$ (24)

将(24)带入(22)，可得

$L_i\left(x\right)=U_i^n+{\frac{\partial U_i}{\partial x}|}_{x_i}\left(x-x_i\right)+\frac{1}{2}{\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}\left({\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\left(x-x_i\right)}^2\text{d}x}\right)+\cdots$ (25)

取

$L_i\left(x\right)=U_i^n+{\frac{\partial U_i}{\partial x}|}_{x_i}\left(x-x_i\right)+\frac{1}{2}{\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}\left({\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\left(x-x_i\right)}^2\text{d}x}\right).$ (26)

重构函数是二阶的，即 $\left|L_i\left(x\right)-U_i\left(x\right)\right|\le O\left({\left(\Delta x\right)}^3\right)$。

记

$L_i\left(x\right)=U_i^n+H\left(x-x_i\right),$ (27)

其中

$H\left(x-x_i\right)={\frac{\partial U_i}{\partial x}|}_{x_i}\left(x-x_i\right)+\frac{1}{2}{\frac{{\partial }^2{U}_i}{\partial x^2}|}_{x_i}\left({\left(x-x_i\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_i}{\left(x-x_i\right)}^2\text{d}x}\right).$ (28)

将一阶偏导和二阶偏导使用中心差商格式近似

${\frac{\partial U_i}{\partial x}|}_{x_i}=\frac{U_{i+1}^n-U_{i-1}^n}{2\Delta x},{\frac{{\partial }^2{U}_i}{\partial x}|}_{x_i}=\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{\Delta x^2},$

将重构函数(27)增加斜率限制器，进一步提高精度，

$L_i\left(x\right)=U_i^n+{\Phi }_{i}H\left(x-x_i\right),$ (29)

其中 ${\Phi }_i$ 为增加的斜率限制器 [8]。

带有斜率限制器的重构在点 $x_i$ 的左右状态，分别为

$\{\begin{array}{l}{L}_i^l=U_i^n+{\Phi }_{i}H\left(x_l-x_i\right),\\ L_i^r=U_i^n+{\Phi }_{i}H\left(x_r-x_i\right).\end{array}$ (30)

令

$\{\begin{array}{l}\min \left(U_{i-1}^n,U_i^n\right)\le L_i^l\le \max \left(U_{i-1}^n,U_i^n\right),\\ \min \left(U_{i+1}^n,U_i^n\right)\le L_i^r\le \max \left(U_{i+1}^n,U_i^n\right),\end{array}$ (31)

以确保不产生新局部极值，将(31)减去 $U_i^n$ 得

$\{\begin{array}{l}\min \left(U_{i-1}^n-U_i^n,0\right)\le L_i^l-U_i^n\le \max \left(U_{i-1}^n-U_i^n,0\right),\\ \min \left(U_{i+1}^n-U_i^n,0\right)\le L_i^r-U_i^n\le \max \left(U_{i+1}^n-U_i^n,0\right).\end{array}$ (32)

当 $U_{i-1}^n-U_i^n>0$， $U_{i+1}^n-U_i^n>0$ 时，有

$\{\begin{array}{l}0<{\Phi }_{i}H\left(x_l-x_i\right)<U_{i-1}^n-U_i^n,\\ 0<{\Phi }_{i}H\left(x_r-x_i\right)<U_{i+1}^n-U_i^n,\end{array}$ (33)

即

$\{\begin{array}{l}0<L_i^l-U_i^n<U_{i-1}^n-U_i^n,\\ 0<L_i^r-U_i^n<U_{i+1}^n-U_i^n.\end{array}$ (34)

因此，选择限制器

$\{\begin{array}{l}{\Phi }_i^l=\min \left(1,\frac{U_{i-1}^n-U_i^n}{H\left(x_l-x_i\right)}\right),\\ {\Phi }_i^r=\min \left(1,\frac{U_{i+1}^n-U_i^n}{H\left(x_r-x_i\right)}\right).\end{array}$ (35)

为了限制 ${\Phi }_i\in \left[0,1\right]$，增加限制条件

$\{\begin{array}{l}{\Phi }_i^l=\max \left(\min \left(1,\frac{U_{i-1}^n-U_i^n}{H\left(x_l-x_i\right)}\right),0\right),\\ {\Phi }_i^r=\max \left(\min \left(1,\frac{U_{i+1}^n-U_i^n}{H\left(x_r-x_i\right)}\right),0\right).\end{array}$ (36)

于是，重构后的左右状态分别为

$\{\begin{array}{l}{L}_i^l=U_i^n+{\Phi }_i^{l}H\left(x_l-x_i\right),\\ L_i^r=U_i^n+{\Phi }_i^{r}H\left(x_r-x_i\right).\end{array}$ (37)

3.3. 高分辨率熵相容格式

将3.2节中的斜率限制器应用到带源项的浅水波方程，得到了具有带源项浅水波方程的二阶精度的高分辨率熵相容格式。

带源项的浅水波方程高分辨率熵相容格式的数值通量为

$F_{i+\frac{1}{2}}^{ec}=F_{i+\frac{1}{2}}^c-\frac{1}{2}{R}_{i+\frac{1}{2}}\left({\Lambda }_{i+\frac{1}{2}}+\left|{\displaystyle \underset{k=1}{\sum }{\omega }_k\left(2s_k-1\right)\Lambda \left(U^{*}\left(s_k\right)\right)}\right|\right)R_{i+\frac{1}{2}}^{\text{T}}{\left[V^{*}\right]}_{i+\frac{1}{2}}.$ (38)

其中 $U^{*},V^{*}$ 是使用斜率限制器重构后的变量值。利用(37)式计算重构， $L_i^l$ 代替左状态 $U_i$， $L_i^r$ 代替右状态 $U_{i+1}$，根据重构后的 $U^{*}$ 值得到熵变量 $V^{*}$，最后将重构值带入(38)，得到带源项浅水波方程的具有二阶精度的熵相容格式数值通量。

4. 数值实验

算例中的参考解使用2000个加密网格的熵稳定格式。EC表示熵相容格式，EC-Limited表示使用斜率限制器的高分辨率熵相容格式。

算例1. 考虑在区间 $\left[-1,1\right]$ 上，连续的底部地势 $b\left(x\right)=1$ 上的大型溃坝问题，初值条件为

$\begin{array}{l}u\left(x,0\right)=\{\begin{array}{l}15x>0,\\ 1x<0,\end{array}\\ h\left(x,0\right)=1.\end{array}$

取100个网格点，计算到 $t=0.1$，取 $g=1$， $CFL=0.1$，使用Neumann边界条件，图1所示，高度h的变化情况，熵相容格式的抹平现象更严重，带有斜率限制器的熵相容格式激波捕捉效果更好。

Figure 1. Variation of height h of large dam break problem

图1. 大型溃坝问题高度h的变化情况

算例2. 在 $\left[0,150\right]$ 上计算溃坝问题，底部地势和初值条件分别为

$\begin{array}{l}b\left(x\right)=\{\begin{array}{l}0.8\left|x-75\right|\le 18.75,\\ 0\text{else},\end{array}\\ h\left(x,0\right)=\{\begin{array}{l}10-b\left(x\right)0\le x\le 75,\\ 2-b\left(x\right)\text{else},\end{array}\\ u\left(x,0\right)=0.\end{array}$

取200个网格点，计算到 $t=4.5$，取 $g=1$， $CFL=0.1$，使用Neumann边界条件，图2所示是 $h+b$ 的变化情况，熵相容格式有轻微的抹平现象，使用斜率限制器的熵相容格式效果更好，在捕捉到激波的情况下，抑制了伪振荡的产生。

算例3. 考虑在 $\left[0,1\right]\times \left[0,1\right]$ 上二维圆形溃坝问题，初值条件为

$\begin{array}{l}h\left(x,y,0\right)=\{\begin{array}{l}2\sqrt{x^2+y^2}<0.5,\\ 1\text{else},\end{array}\\ u\left(x,y,0\right)=0,v\left(x,y,0\right)=0,b\left(x,y\right)=0.\end{array}$

Figure 2. Variation of the total height $h+b$ of the one-dimensional dam break problem

图2. 一维溃坝问题总高度 $h+b$ 的变化情况

Figure 3. Variation of height h of two-dimensional circular dam break problem. (a) 3D graphics of height h in entropy consistent scheme (left) and high-resolution entropy consistent scheme (right); (b) The cross section $y=x$ and the change in height h

图3. 二维圆形溃坝问题高度h的变化情况。(a)熵相容格式(左)和高分辨率熵相容格式(右)高度h的三维图形；(b)截面 $y=x$ 高度h的变化

将空间离散为100 × 100个网格点，取 $CFL=0.1$，计算到 $t=0.2$，使用Neumann边界条件，如图3所示分别为熵相容格式和带有斜率限制器的高分辨率熵相容格式的关于高度h的三维图形和 $y=x$ 的截面，的变化，熵相容格式的抹平现象较严重，带有斜率限制器的熵相容格式改善了这一情况，并且没有伪振荡产生，模拟效果更好。

5. 结论

浅水波方程的研究在实际生活中具有重要意义，通过数值实验表明带有基于MUSCL重构的斜率限制器的熵相容格式在捕捉激波时具有较好的效果，且抹平现象较少，避免了伪振荡的产生。该格式明显优于熵相容格式。但本文未对多维情况进一步推导，多维情况还需要对限制器进行修正，对多个维度重构，得到其高分辨率的熵相容格式。

NOTES

^*通讯作者。