浅水波可以由Saint-Venant[1]方程组来描述, 该方程组被广泛地应用于河流和海洋流动的数值模拟. 一维的Saint-Venant方程组具有如下形式：

$ \begin{equation} \begin{cases} h_t + (hu)_x = 0, \\ (hu)_t + (hu^2+\frac12 gh^2 )_x = -ghB_x. \end{cases} \end{equation} $

(1.1)

其中$ h(x, t) $表示水深, $ u(x, t) $表示速度, $ g $表示重力加速度, 与时间无关的函数$ B(x) $表示底部地形.

我们将许多实际问题描述成一种“静稳态”附近的小扰动. 这里“静稳态”指的是水流速度等于0, 且水面保持水平的一种状态, 即

$ \begin{equation} u = 0, \quad w: = h+B = const. \end{equation} $

(1.2)

这里的$ w $也称作自由面. 特别的, $ h = 0 $表示干底情形. 在数值上精准捕捉这种静稳态附近的小扰动, 是设计格式的一大挑战. “平衡性”（Well-balanced property）的提出可以很好地刻画数值格式捕捉小扰动解的能力. 本文讨论的平衡性是指在数值上能够精确求解稳态解. 同时由于物理的限制, 水深需要保持恒非负$ h\geq0 $. 我们称能够保持水深恒非负的数值格式为“保正性”（Positivity-preserving property）格式. 对于本文研究的求解一维浅水波方程的二阶格式, “带底部”与“不带底部”在数值方法上都需要保持高阶、无震荡以及保持水深非负. 对于这两种情形, 人们常常关心静稳态及其扰动时的波动现象. “不带底部”时各物理量成为常数（包括水深和水速）, 而“带底部”的浅水方程只有速度是常数而水深却不再是常数, 这给数值格式的设计带来很大的困难, 尤其这种困难和保正性困难相互影响. 我们在设计格式时需要保证二阶、无震荡、保正和平衡性.

为方便讨论, 我们将(1.1)式写成向量形式：

$ \begin{equation} U_t + \left(F(U, B)\right)_x = S(U, B). \end{equation} $

(1.3)

其中

$ \begin{gathered} U = \begin{pmatrix} w \\ hu \end{pmatrix}, \quad F(U, B) = \begin{pmatrix} hu \\ \frac{(hu)^2}{w-B}+\frac{g}{2}(w-B)^2 \end{pmatrix}, \quad S(U, B) = \begin{pmatrix} 0 \\ -g(w-B)B_x \end{pmatrix}. \end{gathered} $

浅水波方程组是一类双曲守恒律方程组, 由于这类方程组在实际问题中的广泛应用, 如何求它的数值解一直是一个热点问题. 近年来出现了许多高分辨率重构方法（如WENO[2], CWENO[3, 4]等）. 自1990年Nessyahu和Tadmor[5]提出全离散中心格式之后, 中心格式引起了人们的广泛关注, Kurganov等[6]提出了半离散形式的中心迎风格式. 近年来, Chen和Noelle[7], Bollermann等[8]对半离散的中心迎风格式进行了改进, 使其具有更良好的性质. MUSCL格式[9]是一种求解双曲守恒律方程组的高阶格式, 它在提高了数值精度和分辨率的同时, 也能够保持本质无振荡. Van Leer和Hancock[10]改进了MUSCL格式, 并提出了MUSCL-Hancock方法.

对于带底部地形的浅水波方程组, Kurganov和Petrova[11]提出了一种具有平衡性和保正性的中心迎风格式. 该格式在时间和空间上均达到二阶精度, 但在每个时间步里, 需要在同一个单元边界处求解两次Riemann问题. 本文运用MUSCL-Hancock方法对其进行了改进, 不仅能保证数值精度, 而且在每个时间步里只求解一次Riemann问题, 从而减小计算量. 同时, 在合理的CFL条件下, 本文格式依然具有平衡性和保正性. 经过理论分析和数值实验, 我们验证了本文格式具备平衡性和保正性, 并且是稳定高效可信的；格式在时间和空间都达到了二阶精度, 并且能有效捕捉间断, 具备高分辨率.

本文内容安排如下：第2节, 简单介绍中心迎风格式；第3节, 详细介绍本文设计的新格式, 即基于MUSCL-Hancock方法的中心迎风格式, 该格式在时间和空间上均能达到二阶精度, 同时对水深进行保正重构；第4节是理论分析部分, 证明了基于MUSCL-Hancock方法的中心迎风格式具有保正性和平衡性；第5节, 通过数值算例验证了该格式的鲁棒性.

中心迎风格式是一种求解双曲型方程组的有限体积方法, 该方法具有高分辨率. 该格式利用分片线性重构使其具有空间二阶精度, 通过增加恰当的数值粘性使其具有了迎风的性质. 本节简要介绍其计算方法.

我们将定义域等分为有限体积单元$ I_j = [x_{j-\frac12}, x_{j+\frac12}], (j = 1, 2, ..., N) $, 其单元长度为$ \Delta x : = x_{j+\frac12}-x_{j-\frac12} $. 对(1.3)式在单元$ I_j $上做积分平均, 得到如下半离散格式：

$ \begin{equation} \frac{d U_j(t)}{dt} = -\frac{H_{j+\frac12}(t)-H_{j-\frac12}(t)}{\Delta x}+ S_j(t). \end{equation} $

(2.1)

其中$ S_j $为源项的单元平均的一个合适离散,

$ \begin{equation} S_j(t) \approx \frac{1}{\Delta x}\int_{I_j}S(U(x, t), B(x))dx = \left(0, -g\left(\frac{w_{j+\frac12}^- + w_{j-\frac12}^+ }2-B_j\right) \frac{B_{j+\frac12}-B_{j-\frac12}}{\Delta x}\right)^\top. \end{equation} $

(2.2)

$ H_{j+\frac12} $是中心迎风数值通量：

$ \begin{equation} H_{j+\frac12}(t) = \frac{a^{+}_{j+\frac12}F(U^{-}_{j+\frac12}, B_{j+\frac12})-a^{-}_{j+\frac12}F(U^{+}_{j+\frac12}, B_{j+\frac12}) +a^{+}_{j+\frac12}a^{-}_{j+\frac12}(U^{+}_{j+\frac12}-U^{-}_{j+\frac12})}{a^{+}_{j+\frac12}-a^{-}_{j+\frac12}}. \end{equation} $

(2.3)

这里的$ U^{\pm}_{j+\frac12} $为$ U $在单元边界$ x_{j+\frac12} $处的左右极限值（如图 所示）, 并通过如下线性重构给出

$ \begin{equation} \widetilde q : = q_j + (q_x)_j(x-x_j), \quad x_{j-\frac12}<x<x_{j+\frac12}. \end{equation} $

(2.4)

其中$ q $分别代表$ w $和$ u $. 对于$ U^{\pm}_{j+\frac12} $, 我们令$ h^{\pm}_{j+\frac12} = w^{\pm}_{j+\frac12}-B_{j+\frac12}, (hu)^{\pm}_{j+\frac12} = h^{\pm}_{j+\frac12} u^{\pm}_{j+\frac12} $. 为避免分母过小, 单元平均速度定义为

$ \begin{equation*} u_{j}: = \begin{cases} \frac{({hu})_{j}}{ h_j}, & h_{j} \ge \epsilon, \\ 0, &\rm{其它}. \end{cases} \end{equation*} $

其中$ \epsilon $是一个很小的正数, 本文取$ \epsilon = 10^{-9} $.

为消除数值振荡, 我们对斜率采用如下方法重构：

$ \begin{equation} (q_x)_j = \rm{minmod}\left( \theta \frac{ q_j - q_{j-1}}{\Delta x}, \frac{ q_{j+1}- q_{j-1}}{2\Delta x}, \theta\frac{ q_{j+1} - q_{j}}{\Delta x} \right), \quad \theta \in[1, 2]. \end{equation} $

(2.5)

其中minmod函数定义为:

$ \begin{equation*} \rm{minmod}(z_1, z_2, ...): = \begin{cases} \min\limits_j\lbrace z_j \rbrace, & z_j>0, \forall j\\ \max\limits_j\lbrace z_j \rbrace, & z_j<0, \forall j\\ 0, &\rm{其它}. \end{cases} \end{equation*} $

在数值通量(2.3)式的求解中, 波速估计如下

$ \begin{equation} \begin{aligned} a^{+}_{j+\frac12} = \max\lbrace u^{+}_{j+\frac12}+\sqrt {gh^{+}_{j+\frac12}}, u^{-}_{j+\frac12}+\sqrt{gh^{-}_{j+\frac12}}, 0 \rbrace, \\ a^{-}_{j+\frac12} = \min \lbrace u^{+}_{j+\frac12}- \sqrt{gh^{+}_{j+\frac12}}, u^{-}_{j+\frac12}-\sqrt{gh^{-}_{j+\frac12}}, 0 \rbrace. \end{aligned} \end{equation} $

(2.6)

为了使得重构后的水深非负, 即$ \widetilde w(x)\geq B(x) $, 我们对自由面梯度作如下“保正性修正”[11]：

●  如果$ w^{-}_{j+\frac12}<B_{j+\frac12} $, 则令$ (w_x)_j : = {(B_{j+\frac12}- w_j)}/{\frac{\Delta x}{2}} $

$ \Rightarrow w^{-}_{j+\frac12} = B_{j+\frac12}, \quad w^{+}_{j-\frac12} = 2 w_j - w^{-}_{j+\frac12}; $

●  如果$ w^{+}_{j-\frac12}<B_{j-\frac12} $, 则令$ (w_x)_j : = { ( w_j - B_{j-\frac12})}/{\frac{\Delta x}{2}} $

$ \Rightarrow w^{+}_{j-\frac12} = B_{j-\frac12}, \quad w^{-}_{j+\frac12} = 2 w_j-w^{+}_{j-\frac12}. $

注  如下对半离散格式(2.1)的时间方向离散采用具有二阶精度的强稳定的两步Runge-Kutta法(SSP-RK)[12]：

$ \begin{equation} \begin{aligned} &U_j^{n+1, (1)} = U_j^n + \Delta t R_j^n, \\ &U_j^{n+1, (2)} = U_j^{n+1, (1)} + \Delta t R_j^{n+1, (1)}, \\ &U_j^{n+1} = \frac{U_j^n + U_j^{n+1, (2)}}{2}. \end{aligned} \end{equation} $

(2.7)

在第2节中描述的传统中心迎风格式, 为了在时间方向达到二阶精度, 对半离散格式(2.1)采用两步Runge-Kutta方法[12]更新物理量在新时刻的值. 从而在每一个时间步里, 在单元边界处求解两次Riemann问题. 为了减少Riemann问题的求解次数, 我们采用MUSCL-Hancock方法来更新物理量. 本文提出的数值格式在时间和空间上均达到二阶精度.

MUSCL-Hancock方法的一般步骤为:

1. 重构$ U(x, t_n) $在单元边界处的左右极限值$ U^{n, +}_{j-\frac12}, U^{n, -}_{j+\frac12} $；

2. 计算$ U(x, t_{n+\frac12}) $在单元边界处的左右极限值$ U^{n+\frac12, +}_{j-\frac12}, U^{n+\frac12, -}_{j+\frac12} $；

3. 利用第2步得到的单元边界值计算数值通量$ H^{n+\frac12}_{j\pm\frac12} $和源项$ S_j^{n+\frac12} $, 从而更新$ t_{n+1} $时刻的单元平均值$ U^{n+1}_j $.

上述方法第1步对自由面的重构与传统中心迎风格式相同, 下文3.1节和3.2节详细说明第2步和第3步.

为了计算在$ t_{n+\frac12} $时刻的数值源项$ S_j^{n+\frac12} $(见(2.2))和数值通量$ H_{j+\frac12}^{n+\frac12} $(见(2.3)), 我们在这一节估计$ U(x, t_{n+\frac12}) $在$ x_{j+\frac12} $处的左右极限值（如图 2所示）.

记$ \lambda = \frac{\Delta t}{\Delta x} $. 利用(2.4)式重构$ U(x, t_n) $在单元边界处的左右极限值$ U_{i+\frac12}^{n, \pm} $之后, 我们用下式计算在$ t_{n+\frac12} $时刻的单元边界值

$ \begin{equation} \begin{cases} w^{n+\frac12, \pm}_{j+\frac12} = w^{n, \pm}_{j+\frac12} - \frac{\lambda}{2}\left( (hu)^{n, -}_{j+\frac12}-(h u)^{n, +}_{j-\frac12} \right), \\ ( {hu})^{n+\frac12, \pm}_{j+\frac12} = (h u)^{n, \pm}_{j+\frac12}- \frac{\lambda}{2} \left((\frac12 gh^2+h{ u}^2)^{n, -}_{j+\frac12}-(\frac12 gh^2+h{ u}^2)^{n, +}_{j-\frac12} \right) \\ \quad\quad\quad\qquad\quad -\frac{\lambda}{2}g(w_j^n-B_j)(B_{j+\frac12}-B_{j-\frac12}). \end{cases} \end{equation} $

(3.1)

并记

$ \begin{equation} w_j^{n+\frac12} = \frac{ w^{n+\frac12, -}_{j+\frac12} + w^{n+\frac12, +}_{j-\frac12} }{2}. \end{equation} $

(3.2)

虽然在进行“保正性修正”之后, 有$ w_{j+\frac12}^{n, \pm}\geq B_{j+\frac12} $成立, 但是我们仍然不能保证在半时间步的保正性, 即$ w_{j+\frac12}^{n+\frac12, \pm}\geq B_{j+\frac12} $可能被破坏. 然而, 在一定的CFL条件限制下$ w_j^{n+\frac12}\geq B_j $成立(其证明见引理4.1), 于是我们可以类似地在半时间步作如下的“保正性修正”[11]:

●  如果$ w^{n+\frac12, -}_{j+\frac12}<B_{j+\frac12} $, 则令

$ \begin{equation} w^{n+\frac12, -}_{j+\frac12} = B_{j+\frac12}, \quad w^{n+\frac12, +}_{j-\frac12} = 2 w_j^{n+\frac12}-w^{n+\frac12, -}_{j+\frac12}; \end{equation} $

(3.3)

●  如果$ w^{n+\frac12, +}_{j-\frac12}<B_{j-\frac12} $, 则令

$ \begin{equation} w^{n+\frac12, +}_{j-\frac12} = B_{j-\frac12}, \quad w^{n+\frac12, -}_{j+\frac12} = 2 w_j^{n+\frac12}-w^{n+\frac12, +}_{j-\frac12}. \end{equation} $

(3.4)

利用已经估算出来的守恒量$ U^{n+\frac12, \pm}_{j+\frac12} $, 我们可以更新其在$ t_{n+1} $时刻的单元平均值

$ \begin{equation} U^{n+1}_j = U^n_j - \lambda\left(H^{n+\frac12}_{j+\frac12}-H^{n+\frac12}_{j-\frac12}\right)+\Delta t S^{n+\frac12}_j. \end{equation} $

(3.5)

其中数值通量$ H^{n+\frac12}_{j+\frac12} $与数值源项$ S^{n+\frac12}_{j} $分别由(2.3)式和(2.2) 式在$ t_{n+\frac12} $时刻求得, 即利用$ U^{n+\frac12, \pm}_{j+\frac12} $去实现.

我们在设计求解带底部的一维浅水波方程组的数值格式时, 常常关注两种重要性质: 一是格式要能保证水深非负$ h\geq0 $, 即“保正性”；二是要能在静稳态条件下, 维持稳态解, 即“平衡性”. 本节我们理论证明了本文设计的MUSCL-Hancock中心迎风格式具有这两种性质.

保正性(Positivity-preserving property)证明包含两个步骤, 首先是考虑在$ t_{n+\frac12} $时刻的保正性, 其次考虑在$ t_{n+1} $时刻的保正性.

下面引理给出$ t_{n+\frac12} $时刻的保正性.

引理4.1   如果时间步长满足CFL条件限制, $ \lambda \max\limits_j |u^{n, \pm}_{j\pm\frac12}|\leq \frac13 $, 则由公式(3.1)计算所得在$ t_{n+\frac12} $时刻的单元平均水深恒非负, 即$ w^{n+\frac12}_{j} \ge B_j $.

证   利用(3.2)式可以得到

$ \begin{equation} \begin{aligned} w^{n+\frac12}_{j} - B_j = & \frac{w^{n+\frac12, -}_{j+\frac12} + w^{n+\frac12, +}_{j-\frac12}}{2} - \frac{B_{j+\frac12} + B_{j-\frac12}}{2} \\ = & \frac{w^{n, -}_{j+\frac12} + w^{n, +}_{j-\frac12} - \lambda ((hu)^{n, -}_{j+\frac12}-(h u)^{n, +}_{j-\frac12})}{2} - \frac{B_{j+\frac12} + B_{j-\frac12}}{2} \\ = & \frac{h^{n, -}_{j+\frac12} + h^{n, +}_{j-\frac12} - \lambda ((hu)^{n, -}_{j+\frac12}-(h u)^{n, +}_{j-\frac12})}{2} \\ = & \frac12\left(1 - \lambda u^{n, -}_{j+\frac12} \right) h^{n, -}_{j+\frac12} + \frac12\left(1 + \lambda u^{n, +}_{j-\frac12}\right)h^{n, +}_{j-\frac12} \\ \geq& 0 \end{aligned} \end{equation} $

(4.1)

其中第二个等式中用到了(3.1)式；由CFL条件, 有$ 1 - \lambda u^{n, -}_{j+\frac12}>0 $, $ 1 + \lambda u^{n, +}_{j-\frac12}> 0 $, 从而最后的不等式成立. 引理得证.

下面讨论$ t_{n+1} $时刻的保正性.

定理4.1   如果时间步长满足CFL条件限制, $ \lambda \max\limits_j |a^{n+\frac12, \pm}_{j\pm\frac12}|\leq \frac13 $, 则由公式(3.5)更新的在$ t_{n+1} $时刻的单元平均水深恒非负, 即$ w^{n+1}_{j} \ge B_j $.

证   我们仅考虑(3.5)式的第一个分量式, 首先将其减去底部$ B_j $, 类似(4.1)式, 再代入由(2.3)式计算出的数值通量, 可以得到

$ \begin{equation} \begin{aligned} w^{n+1}_j - B_j = & h_j^n - \lambda \left(H^{n+\frac12, (1)}_{j+\frac12}-H^{n+\frac12, (1)}_{j-\frac12}\right) \\ = & h_j^n - \lambda a^{n+\frac12, +}_{j+\frac12} \frac{u^{n+\frac12, -}_{j+\frac12}-a^{n+\frac12, -}_{j+\frac12}}{a^{n+\frac12, +}_{j+\frac12}-a^{n+\frac12, -}_{j+\frac12}} h^{n+\frac12, -}_{j+\frac12} -\lambda a^{n+\frac12, -}_{j+\frac12}\frac{a^{n+\frac12, +}_{j+\frac12}-u^{n+\frac12, +}_{j+\frac12}}{a^{n+\frac12, +}_{j+\frac12}-a^{n+\frac12, -}_{j+\frac12}} h^{n+\frac12, +}_{j+\frac12} \\ &+ \lambda a^{n+\frac12, +}_{j-\frac12} \frac{u^{n+\frac12, -}_{j-\frac12}-a^{n+\frac12, -}_{j-\frac12}}{a^{n+\frac12, +}_{j-\frac12}-a^{n+\frac12, -}_{j-\frac12}} h^{n+\frac12, -}_{j-\frac12} +\lambda a^{n+\frac12, -}_{j-\frac12}\frac{a^{n+\frac12, +}_{j-\frac12}-u^{n+\frac12, +}_{j-\frac12}}{a^{n+\frac12, +}_{j-\frac12}-a^{n+\frac12, -}_{j+\frac12}} h^{n+\frac12, +}_{j-\frac12}. \end{aligned} \end{equation} $

(4.2)

由引理4.1和第3.1节的保正性修正, 可以得到$ h^{n+\frac12, \pm}_{j \pm \frac12}\geq0 $. 又由波速估计(2.6)式可知$ a^{n+\frac12, +}_{j\pm\frac12}\ge 0 $, $ a^{n+\frac12, -}_{j\pm\frac12}\le 0 $, 且$ a^{n+\frac12, +}_{j\pm\frac12}-a^{n+\frac12, -}_{j\pm\frac12} \geq a^{n+\frac12, +}_{j\pm\frac12}-u^{n+\frac12, +}_{j\pm\frac12} \ge 0 $及$ a^{n+\frac12, +}_{j\pm\frac12}-a^{n+\frac12, -}_{j\pm\frac12} \geq u^{n+\frac12, -}_{j\pm\frac12}-a^{n+\frac12, -}_{j\pm\frac12} \ge 0 $. 因此

$ \begin{equation} w^{n+1}_j - B_j \geq h_j^n - \frac13 ( h^{n+\frac12, -}_{j+\frac12} + h^{n+\frac12, +}_{j-\frac12}) = h_j^n - \frac23 h_j^{n+\frac12}. \end{equation} $

(4.3)

(4.3)中的等号利用了第3.1节中(3.1)式和修正步的守恒性. 由(4.1)式知$ h_j^{n+\frac12} \leq \frac43 h^{n}_j $, 于是$ w^{n+1}_j - B_j \geq \frac19 h_j^n\geq0 $. 定理得证.

“平衡性”（Well-balanced property）可以很好的刻画数值格式捕捉稳态附近的小扰动解的能力. 本节讨论格式的平衡性. 如同保正性分析, 这里同样也包含两个步骤：首先是考虑在单元边界处的不变性, 即$ U_{j+\frac12}^{n+\frac12} = U_{j+\frac12}^{n} $；然后考虑单元中心处的不变性, 即$ U_j^{n+1} = U_j^n $.

定理4.2   如果初值满足$ w_j = const, u_j = 0, \forall j $, 那么可以得到$ U_j^{n+1} = U_j^{n}, \forall j $.

证   首先考虑在$ t_{n+\frac12} $时刻物理量的计算. 由静稳态初始条件和初始数据重构可以保证$ w^{n, \pm}_{j\pm\frac12} = w_j^n = const $且$ u^{n, \pm}_{j\pm\frac12} = 0 $, 从而$ h^{n, \pm}_{j\pm\frac12} = w_j^n-B_{j+\frac12} $. 利用(3.1)式可得

$ \begin{equation} w^{n+\frac12, \pm}_{j+\frac12} = w^{n, \pm}_{j+\frac12} = w_j^n, \end{equation} $

(4.4)

以及

$ \begin{equation} \begin{aligned} ( {hu})^{n+\frac12, \pm}_{j+\frac12} = & - \frac{\lambda}{2} \left((\frac12 gh^2)^{n, -}_{j+\frac12}-(\frac12 gh^2)^{n, +}_{j-\frac12} + g(w^n_j-B_j)(B_{j+\frac12}-B_{j-\frac12}) \right) \\ = & - \frac{\lambda}{2} gh^n_j\left(h^{n, -}_{j+\frac12} + B_{j+\frac12}- h^{n, +}_{j-\frac12} - B_{j-\frac12}\right) = 0, \end{aligned} \end{equation} $

(4.5)

于是$ u^{n+\frac12, \pm}_{j\pm\frac12} = 0 $和$ h^{n+\frac12, \pm}_{j\pm\frac12} = w_j^n-B_{j+\frac12} $.

下面考虑$ t_{n+1} $时刻物理量的更新. 将$ U_{j+\frac12}^{n+\frac12, \pm} $代入(2.2)式和(2.3)式计算可得数值源项和数值通量分别为

$ \begin{equation} S_j^{n+\frac12} = \begin{pmatrix} 0 \\ -g (w_j^n -B_j) \frac{B_{j+\frac12}-B_{j-\frac12}}{\Delta x} \end{pmatrix}, \quad H_{j+\frac12}^{n+\frac12} = \begin{pmatrix} 0 \\ \frac{g}{2}(w_j^n-B_{j+\frac12})^2 \end{pmatrix}, \end{equation} $

(4.6)

代入(3.5)式可得

$ \begin{equation} \begin{aligned} U^{n+1}_j - U^n_j = & - \lambda\left(H^{n+\frac12}_{j+\frac12}-H^{n+\frac12}_{j-\frac12}\right)+\Delta t S^{n+\frac12}_j \\ = & -\lambda \begin{pmatrix} 0 \\ \frac{g}{2}(w_j^n-B_{j+\frac12})^2 - \frac{g}{2}(w_j^n-B_{j-\frac12})^2 + g (w_j^n -B_j) (B_{j+\frac12}-B_{j-\frac12})\end{pmatrix}\\ = &\begin{pmatrix}0\\0\end{pmatrix}. \end{aligned} \end{equation} $

(4.7)

于是$ U^{n+1}_j = U^n_j $. 定理得证.

本文在第2节通过(2.4)式的分片线性重构使得格式在空间上达到二阶精度, 同时第2节的(2.7)式使用两步Runge-Kutta方法使得格式在时间上达到二阶精度. 本节我们通过一系列数值算例来验证本文格式的精度、有效性和鲁棒性. 第一个算例验证格式的数值精度；第二个算例验证格式捕捉静稳态解的小扰动的能力；第三个算例验证格式的保正性以及处理干界面的能力. 最后我们来验证格式对于间断底部上复合波的高分辨率. 在所有的计算中, 斜率限制器(2.5)中的参数$ \theta $均取为1.3. 由于这些算例缺少解析解, 对于以下所有数值算例, 我们均取在同等条件下的密网格（网格数为3200）上的数值结果作为参考解.

本算例选自参考文献[13], 旨在验证本文格式的数值精度. 本算例中使用的初始数据均为连续函数. 水深$ h $, 动量$ hu $的初始值分别为

$ \begin{equation*} h(x, 0) = 5+e^{\cos(2\pi x)}, \quad hu(x, 0) = \sin(\cos(2\pi x)), \quad \forall x\in [0, 1] \end{equation*} $

底部函数$ B(x) $为

$ \begin{equation*} B(x) = \sin^2(\pi x), \quad \forall x\in [0, 1]. \end{equation*} $

重力加速度$ g = 9.812 $. 为了比较精度, 模拟结束时间选在解出现间断之前, 这里取$ t_{\rm stop} = 0.1 $. 本算例采用周期边界条件. 我们分别用传统的MUSCL中心迎风格式和本文设计的MUSCL-Hancock中心迎风格式进行模拟, 当计算单元数分别为$ 25, 50, 100, 200, 400 $和$ 800 $时, 表 1显示的是的$ L^1 $误差和精度, 表 2显示的是两种格式程序计算时间. 结果表明这两种格式均能达到二阶数值精度, 同时我们还发现本文格式的$ L^1 $误差约为已有的MUSCL中心迎风格式的$ 1/3 $, 计算时间减少了约$ 1/3 $.

本算例选自文献[11], 我们将用它来验证格式的平衡性, 以及捕捉稳态解的微小扰动传播的能力.

这里底部定义为一个复杂的分段函数

$ \begin{equation} B(x) = \begin{cases} 10(x-3), &0.3 \le x \le 0.4, \\ 1-0.0025\sin^2(25\pi(x-0.4)), &0.4 \le x \le 0.6, \\ -10(x-0.7), &0.6 \le x\le 0.7, \\ 0, &\rm{其它}. \end{cases} \end{equation} $

(5.1)

初始条件为

$ \begin{equation} (w(x, 0), u(x, 0)) = \begin{cases} (1+\epsilon, 0), &0.1<x<0.2, \\ (1, 0), &\rm{其它}. \end{cases} \end{equation} $

(5.2)

常数$ g = 1 $, 计算区域为[-1, 1]. $ \epsilon $为一个小参数. 当取$ \epsilon = 0 $, 表示静稳态. 我们选择程序结束时间为$ t_{\rm stop} = 10 $, 并把在单元数为200时的水深$ h $和动量$ hu $的$ L^{\infty} $误差显示在表 3中. 结果显示无论是$ h $还是$ hu $, 数值误差都处于机器精度范围内, 这从数值上证明了这两种中心迎风格式都具有平衡性. 我们也模拟在稳态附近的小扰动问题. 在$ [0.1, 0.2] $上的微小扰动参数取为$ \epsilon = 0.001 $, 结束时间$ t_{\rm stop} = 1 $, 单元数为100和200时的计算结果显示在图 3中. 结果表明两种迎风格式都能很好的捕获这种小扰动的传播. 我们需要强调的是, 数值结果也显示了MUSCL-Hancock中心迎风格式比MUSCL中心迎风格式具有更好的分辨率.

本算例选自文献[11], 底部为阶梯函数：

$ \begin{equation} B(x) = \begin{cases} 2, & x \le 0.5 , \\ 0.1, & x>0.5. \end{cases} \end{equation} $

(5.3)

初始值：

$ \begin{equation} (w(x, 0), u(x, 0)) = \begin{cases} (2.222, -1), & x \le 0.5, \\ (0.8246, -1.6359), & x>0.5. \end{cases} \end{equation} $

(5.4)

定义域为[0, 1], 常数$ g = 2 $, 结束时间为$ t_{\rm stop} = 0.2 $. 该问题的精确解由一个稀疏波, 接触间断和激波组成. 使用本文设计的MUSCL-Hancock中心迎风格式分别在100个单元和200个单元上进行计算, 结果显示在图 4中. 计算结果表明本文格式能够很好的捕捉这几个波. 在底部跳跃点$ x = 0.5 $附近, 水深已经非常接近于0. 由于格式具有保正性, 所以对于这种几乎干底问题也能够进行很好的模拟.

本算例选自文献[11], 底部是阶梯函数：

$ \begin{equation} B(x) = \begin{cases} 1.5, & x \le 0.5, \\ 1.1, & x>0.5. \end{cases} \end{equation} $

(5.5)

初始Riemann条件：

$ \begin{equation} (w(x, 0), u(x, 0)) = \begin{cases} (5, 1), & x \le 0.5, \\ (1.6, -2), & x>0.5. \end{cases} \end{equation} $

(5.6)

取常数$ g = 2 $, 定义域为[0, 1], 结束时间为$ t_{\rm stop} = 0.15 $. 该问题的精确解包含一个复合三波, 它由声速稀疏波连接到一个接触间断, 并伴随着零速度激波和另一个接触间断. 我们可以从图 5中明显地看出, 本文格式具有非常高的整体分辨率. 数值模拟结果在底部间断处有一个突降和突升.

关于求解带底部地形的一维浅水波方程, 本文给出基于MUSCL-Hancock方法的中心迎风格式. 相对于传统的MUSCL中心迎风格式而言, 在每个时间步$ [t_n, t_{n+1}] $, 本文设计的格式只在$ t_{n+\frac12} $时刻求解Riemann问题, 并更新数值通量和数值源项, 从而达到减少计算量的目的. 为了格式的保正性, 在每个时间步, 我们分别在$ t_n $时刻和$ t_{n+\frac12} $时刻对自由面的斜率进行了修正. 这种修正对于数据重构具有保正性和守恒性. 然后我们严格证明了在一定的CFL条件限制下, 格式对于解的更新具有平衡性与保正性. 数值算例也验证了这些性质. 另外, 数值结果表明, 与MUSCL格式相比, 本文格式具有更高的分辨率. 我们下一步将考虑格式的高维推广, 以及更一般的稳态解的捕捉能力, 例如动稳态解.