不规则沟渠浅水流动的高阶平衡有限差分AWENO格式

High-Order Well-Balanced Finite Difference AWENO Schemes for Shallow Water Flows along Channels with Irregular Geometry

摘要: 在本文中，我们提出了具有不规则几何形状和非平坦底部地形的开放沟渠中浅水方程的高阶有限差分加权本质无振荡(Alternative Weighted Essentially Non-Oscillatory) AWENO格式。所提出的格式保持了静水稳态解的良好平衡性质，即在通量梯度和源项之间存在精确平衡时，在离散状态保持稳态。与具有恒定横截面的传统浅水方程相比，由于沟渠不规则几何形状引起的影响，构建良好平衡格式并不是一项简单的工作。为了保持良好平衡性质，我们首先重新构造源项，然后使用具有Lax-Friedrichs通量的静水重构方法，最后借助一种新颖的源项逼近方法离散源项。基准数值示例被应用来验证所得格式的良好性能：平衡性能，高阶精度，以及对不连续解的高分辨率。

Abstract: In this paper, we present high-order finite alternative weighted essentially non-oscillatory (AWENO) schemes for shallow water flows along open channels with irregular geometry and over a non-flat bottom topography. The proposed schemes maintain the well-balanced property for the still water steady-state solutions, namely preserving steady state at the discrete level, when there is an exact balance between the flux gradient and the source term. Compared with the traditional shallow water equations with constant cross-section, the construction of the well-balanced schemes is not trivial work due to the effect induced by the irregular geometry of the channels. To preserve the well-balanced property, we first reformulate the source term, then use the hydrostatic reconstruction (HR) method with Lax-Friedrichs (LF) flux, and finally discrete the source term with the help of a novel source term approximation. Benchmark numerical examples are applied to validate the good performances of the resulting schemes: well-balanced property, high order accuracy, and high resolution for the discontinuous solutions.

1. 引言

浅水方程作为“浅水”自由水面流动的数学模型，在河流、河道等实际水动力问题的模拟中具有重要作用。浅水流动的数值模拟在水利和海岸工程中有着广泛的应用[1] [2]。通过对Navier-Stokes方程的深度积分[3]，可以推导出沿不规则几何形状的沟渠非平坦底部地形的一维浅水方程[4]，其形式为：

$\begin{array}{l}{H}_t+Q_x=0,\\ Q_t+{\left(\frac{Q^2}{H}+\frac{1}{2}g\sigma h^2\right)}_x=\frac{1}{2}g{h}^2{\sigma }_x-g\sigma h{b}_x,\end{array}$ (1.1)

其中${\sigma }_x$ 表示沟渠横截面的宽度，b为河底地形，h为水面高度，$H=\sigma h$ 为湿截面，$Q=Hu$ 为质量流量，u为速度，g为重力常数。

模型(1.1)属于带有源项的双曲守恒律(也称双曲平衡律)，并且允许给出的静水定常解：

$u=0,h+b=const,$ (1.2)

其中源项与通量的梯度完全平衡。包括非平坦地形上的非线性浅水方程和引力场作用下的欧拉方程在内的双曲平衡律的数值模拟面临的一个主要挑战是，标准数值格式在离散水平上破坏了稳态下通量梯度和源项之间的精确平衡，并且即使在极细的网格上也会导致稳态附近的伪振荡，因为截断误差占据了小扰动的主要部分，从而导致了解的大误差。为了解决这些问题，我们提出了平衡格式，在离散水平上保持静水定常解直至舍入误差，并且在相对粗糙的网格上也能解决定常解下的小扰动。

在本次工作中，我们重点研究了高阶有限差分AWENO格式[5] [6]，该方案具有以下优点：1) 实现简单；2) 计算效率高；3) 容易拓展到多维模型；4) 任意单调通量在计算单元边界处数值通量有较好的适用性。

2. 高阶有限差分AWENO格式

在本节中，我们将简要回顾有限差分AWENO格式[5] [6]。首先我们介绍AWENO格式下的一些标准符号。区间$I=\left[a,b\right]$ 被分成N个子区间，并且我们用$I_j=\left[x_{j-1/2},x_{j+1/2}\right],j=1,\cdots ,N$ 来表示单元。为简单起见，我们首先假设网格均匀分布，单元格的大小为$\Delta x=x_{j+1}-x_j$ ，并且我们使用$x_{j+1/2}=x_j+\Delta x/2$ 来表示单元$I_i$ 的中心。在不丧失一般性的情况下，我们首先考虑一维标量的情况。

$u_t+f{\left(u\right)}_x=0.$

半离散格式为

$\frac{\text{d}}{\text{d}t}{u}_j\left(t\right)=-\frac{1}{\Delta x}\left({\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}\right),$

其中${\widehat{f}}_{i+\frac{1}{2}}$ 为数值通量。

2.1. 五阶WENO插值

我们简要回顾守恒变量u的五阶WENO插值过程的必要组成部分。单元边界处的值$u_{i+\frac{1}{2}}^{-}$ 在目标值$x_{i+\frac{1}{2}}$ 的左侧改写成三阶插值多项式$q^{\left(k\right)}\left(x\right)$ 的凸组合，其中在每个子模板中的非线性权为${\omega }_k$ ，子模板取$S_k=\left[x_{i-2+k},\cdots ,x_{i+k}\right],k=0,1,2$ ，在全局模板$S^5={\displaystyle {\cup }_{k=0}^2}$ 中。例如

$u_{i+\frac{1}{2}}^{-}={ℐ}_{i+\frac{1}{2}}^{-}\left[u\right]:={\displaystyle \sum _{k=0}^2{\omega }_k}{u}_{i+\frac{1}{2}}^{\left(k\right)},$

并且

$u_{i+\frac{1}{2}}^{\left(0\right)}=\frac{3}{8}{u}_{i-2}-\frac{5}{4}{u}_{i+1}+\frac{15}{8}{u}_i,u_{i+\frac{1}{2}}^{\left(1\right)}=-\frac{1}{8}{u}_{i-1}+\frac{3}{4}{u}_i+\frac{3}{8}{u}_{i+1},u_{i+\frac{1}{2}}^{\left(2\right)}=\frac{3}{8}{u}_i+\frac{3}{4}{u}_{i+1}-\frac{1}{8}{u}_{i+2},$

其中$\left\{u_i\right\}$ 是单元中心的点值，并且$u_{i+\frac{1}{2}}^{\left(k\right)}=u^{\left(k\right)}\left(x_{i+\frac{1}{2}}\right)$ 。改进的非线性权${\omega }_k$ [6]被定义为

${\omega }_k=\frac{{\alpha }_k}{{\displaystyle \sum _{s=0}^2{\alpha }_s}},{\alpha }_k=d_k\left(1+{\left(\frac{{\tau }_5}{{\beta }_k+\epsilon }\right)}^p\right),$

其中$\{d_0=\frac{1}{16},d_1=\frac{5}{8},d_2=\frac{5}{16}\}$ 是线性权。p是权参数并且灵敏参数$\epsilon$ 被用来避免0的除法。在这项研究之中，$p=2$ 和$\epsilon ={10}^{-12}$ 被选择进行数值算例。子模板$S_k$ 的局部低阶光滑指示子${\beta }_k$ 为

${\beta }_k={\displaystyle \sum _{\ell =1}^2\Delta }{x}^{2\ell -1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\left(\frac{{\partial }^{\ell }}{\partial x^{\ell }}{u}^{\left(k\right)}\left(x\right)\right)}^2\text{d}x},$

其中${\tau }_5=\left|{\beta }_0-{\beta }_2\right|$ 。

在上面的讨论中，由于全局模板$S_5$ 偏向于目标点$x_{j{+}_2^1}$ 的左侧，右单元边界值$u_{i+\frac{1}{2}}^{+}={ℒ}_{i+\frac{1}{2}}^{+}\left[u\right]$ 关于目标值$x_{j+\frac{1}{2}}$ 镜像对称于$u_{j+\frac{1}{2}}^{-}$ 。

2.2. 生成数值通量

对于五阶AWENO格式，数值通量${\widehat{f}}_{i+\frac{1}{2}}$ 可以表示为：

${\widehat{f}}_{i+\frac{1}{2}}=ℱ\left(u_{i+\frac{1}{2}}^{-},u_{i+\frac{1}{2}}^{+}\right)+{\mathcal{D}}_{i+\frac{1}{2}}\left[f\right],$

其中$ℱ\left(.,.\right)$ 是单调通量，${\mathcal{D}}_{i+\frac{1}{2}}$ 是高阶中心差分算子，$u_{i+\frac{1}{2}}^{\pm }$ 通过五阶WENO插值获得。在这次研究中，我们使用Lax-Friedrichs通量求解方程组。为保证数值稳定，我们采用简单的中心有限差分法近似高阶修正导数项，全局阶为$O\left(\Delta x^6\right)$ ，如下：

${\mathcal{D}}_{i+\frac{1}{2}}\left[f\right]=-\frac{1}{24}{\Delta x^2{f}_{xx}|}_{i+\frac{1}{2}}+\frac{7}{5760}{\Delta x^4{f}_{xxxx}|}_{i+\frac{1}{2}},$

其中

$\begin{array}{l}{\Delta x^2{f}_{xx}|}_{i+\frac{1}{2}}=\frac{1}{48}\left(-5f_{i-2}+39f_{i-1}-34f_i-34f_{i+1}+39f_{i+2}-5f_{i+3}\right)+O\left(\Delta x^6\right),\\ {\Delta x^4{f}_{xxxx}|}_{i+\frac{1}{2}}=\frac{1}{2}\left(f_{i-2}-3f_{i-1}+2f_i+2f_{i+1}-3f_{i+2}+f_{i+3}\right)+O\left(\Delta x^6\right).\end{array}$

2.3. 推广到方程组

我们对方程(1.1)的齐次形式使用AWENO格式：

$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}=0.$

首先，守恒变量$U_{i+\frac{1}{2}}$ 通过Roe平均左特征向量$L_{i+\frac{1}{2}}$ 投影到特征场中：

${\tilde{U}}_{i+l}=L_{i+\frac{1}{2}}{U}_{i+l},l=-2,\cdots ,3.$

然后，我们使用WENO插值程序计算边界值${\tilde{U}}_{i+\frac{1}{2}}^{-}$ ：

${\tilde{U}}_{i+\frac{1}{2}}^{-}={\displaystyle \sum _{k=0}^2{\omega }_k}{\tilde{U}}^{\left(k\right)}\left(x_{i+\frac{1}{2}}\right).$

最后，插值${\tilde{U}}_{i+\frac{1}{2}}^{-}$ 通过右特征向量$R_{i+\frac{1}{2}}$ 投影回物理空间：

$U_{i+\frac{1}{2}}^{-}=R_{i+\frac{1}{2}}{\tilde{U}}_{i+\frac{1}{2}}^{-}.$

右单元边界值$U_{i+\frac{1}{2}}^{+}$ 可以通过与标量类似的情况计算。数值通量${\widehat{F}}_{i+\frac{1}{2}}$ 可以被近似为：

${\widehat{F}}_{i+\frac{1}{2}}=ℱ\left(U_{i+\frac{1}{2}}^{-},U_{i+\frac{1}{2}}^{+}\right)+{\mathcal{D}}_{i+\frac{1}{2}}\left[F\right].$

3. 良好平衡的AWENO格式

3.1. 改写方程

在本节中，我们提出了沿静水稳态沟渠的一维浅水方程的良好平衡有限差分AWENO格式(1.2)。然而，原始系统(1.1)的标准AWENO格式并不能直接导致良好平衡格式。在此，我们将原方程重新表述为：

$\begin{array}{l}{H}_t+Q_x=0,\\ Q_t+{\left(\frac{Q^2}{H}+\frac{1}{2}g\sigma h^2\right)}_x=\frac{1}{2}g{\left(h+b\right)}^2{\sigma }_x-g\left(h+b\right){\left(\sigma b\right)}_x+\frac{1}{2}g{\left(\sigma b^2\right)}_x,\end{array}$ (3.1)

其中我们将源项$\frac{1}{2}g{h}^2{\sigma }_x-g\sigma h{b}_x$ 改写为$\frac{1}{2}g{\left(h+b\right)}^2{\sigma }_x-g\left(h+b\right){\left(\sigma b\right)}_x+\frac{1}{2}g{\left(\sigma b^2\right)}_x$ 。定义

$\begin{array}{l}U={\left(H,Q\right)}^{\top },F\left(U,\sigma \right)={\left(\begin{array}{c}Q\\ \frac{Q^2}{H}+\frac{1}{2}g\sigma h^2\end{array}\right)}_x,\\ S\left(U,\sigma ,\sigma b,\sigma b^2\right)=\frac{1}{2}g{\left(h+b\right)}^2{\left(\begin{array}{c}0\\ \sigma \end{array}\right)}_x-g\left(h+b\right){\left(\begin{array}{c}0\\ \sigma b\end{array}\right)}_x+\frac{1}{2}g{\left(\begin{array}{c}0\\ \sigma b^2\end{array}\right)}_x.\end{array}$

上述方程的紧致形式为：

$U_t+F{\left(U,\sigma \right)}_x=S\left(U,\sigma ,\sigma b,\sigma b^2\right).$

一维守恒方程组的半离散形式，使用有限差分AWENO格式，可以总结为：

${\frac{\text{d}U}{\text{d}t}|}_{x=x_i}+\frac{{\widehat{F}}_{i+\frac{1}{2}}-{\widehat{F}}_{i-\frac{1}{2}}}{\Delta x}=\frac{{\widehat{S}}_{i+\frac{1}{2}}-{\widehat{S}}_{i-\frac{1}{2}}}{\Delta x},$(3.2)

其中${\widehat{F}}_{i+\frac{1}{2}}$ 为数值通量。

3.2. 应用静水重构方法构造数值通量

定义一个新的向量$B={\left(0,b\left(x\right)\right)}^{\top }$ ，$S_1={\left(0,\sigma \left(x\right)\right)}^{\top }$ ，$Q={\left(h,hu\right)}^{\top }$ 其插值计算公式为：

$B_{i+\frac{1}{2}}^{\pm }=P^{-1}{ℐ}_{i+\frac{1}{2}}^{\pm }\left[PB\right],$

${\left(S_1\right)}_{i+\frac{1}{2}}^{\pm }=P^{-1}{ℐ}_{i+\frac{1}{2}}^{\pm }\left[P{S}_1\right],$

其中线性WENO插值算子${ℐ}_{i+\frac{1}{2}}^{\pm }\left[PB\right]$ 和${ℐ}_{i+\frac{1}{2}}^{\pm }\left[P{S}_1\right]$ 使用了${ℐ}_{i+\frac{1}{2}}^{\pm }\left[U\right]$ 的一般非线性权，其中P是一个常数排列矩阵：

$P=(\begin{array}{cc}0& 1\\ 1& 0\end{array}).$

显然有

$Q_{i+\frac{1}{2}}^{\pm }+P{B}_{i+\frac{1}{2}}^{\pm }={ℐ}_{i+\frac{1}{2}}^{\pm }\left[Q+PB\right].$

使用HR方法[7] [8] [9]，我们设

$h_{i+\frac{1}{2}}^{*,\pm }=\max \left(0,h_{i+\frac{1}{2}}^{\pm }+b_{i+\frac{1}{2}}^{\pm }-\max \left(b_{i+\frac{1}{2}}^{+},b_{i+\frac{1}{2}}^{-}\right)\right),$

${\sigma }_{i+\frac{1}{2}}^{*,\pm }=\max \left(0,\max \left({\sigma }_{i+\frac{1}{2}}^{+},{\sigma }_{i+\frac{1}{2}}^{-}\right)\right),$

左右插值的值$U_{i+\frac{1}{2}}^{\pm }$ 修正为：

$U_{i+\frac{1}{2}}^{*,\pm }=(\begin{array}{c}{\sigma }_{i+\frac{1}{2}}^{*,\pm }{h}_{i+\frac{1}{2}}^{*,\pm }\\ {\sigma }_{i+\frac{1}{2}}^{*,\pm }{h}_{i+\frac{1}{2}}^{*,\pm }{u}_{i+\frac{1}{2}}^{\pm }\end{array}).$

数值通量${\widehat{F}}_{i+\frac{1}{2}}$ 修正为${\widehat{F}}_{i+\frac{1}{2}}^{-}$ 和${\widehat{F}}_{i+\frac{1}{2}}^{+}$ ，

${\widehat{F}}_{i+\frac{1}{2}}^{\pm }=ℱ\left(U_{i+\frac{1}{2}}^{*,-},U_{i+\frac{1}{2}}^{*,+}\right)+{\mathscr{H}}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[F\right],$

其中静水修正项为：

${\mathscr{H}}_{i+\frac{1}{2}}^{\pm }=(\begin{array}{c}0\\ \frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2-\frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{*,\pm }{\left(h_{i+\frac{1}{2}}^{*,\pm }\right)}^2\end{array}).$

3.3. 源项离散

源项进一步修改为：

$S\left(U,\sigma ,\sigma b,\sigma b^2\right)=\frac{1}{2}g{\left(h+b\right)}^2\frac{\partial S_1}{\partial x}-g\left(h+b\right)\frac{\partial \left(B\circ S_1\right)}{\partial x}+\frac{1}{2}g\frac{\partial \left(B\circ B\circ S_1\right)}{\partial x},$

其中符号$\circ$ 定义为向量中对应分量相乘。为了与数值通量${\widehat{F}}_{{}_{i+\frac{1}{2}}}$ 保持一致，我们在单元边界处$x_{i+\frac{1}{2}}$ 近似源项[10]。

${\widehat{S}}_{i+\frac{1}{2}}^{\pm }=\frac{1}{2}g{\left(h_i+b_i\right)}^2{\left({\widehat{S}}_1\right)}_{i+\frac{1}{2}}^{\pm }-g\left(h_i+b_i\right){\left(\widehat{B°S_1}\right)}_{i+\frac{1}{2}}^{\pm }+\frac{1}{2}g{\left(\widehat{B°B°S_1}\right)}_{i+\frac{1}{2}}^{\pm },$

定义其中的量

$\begin{array}{l}{\left({\widehat{S}}_1\right)}_{i+\frac{1}{2}}^{\pm }={\left({\text{S}}_{\text{1}}\right)}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[{\text{S}}_{\text{1}}\right],\\ {\left(\widehat{\text{B}°{\text{S}}_1}\right)}_{i+\frac{1}{2}}^{\pm }={\text{B}}_{i+\frac{1}{2}}^{\pm }°{\left({\text{S}}_1\right)}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[\text{B}°{\text{S}}_1\right],\\ {\left(\widehat{\text{B}°\text{B}°{\text{S}}_1}\right)}_{i+\frac{1}{2}}^{\pm }={\text{B}}_{i+\frac{1}{2}}^{\pm }°{\text{B}}_{i+\frac{1}{2}}^{\pm }°{\left({\text{S}}_1\right)}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[\text{B}°\text{B}°{\text{S}}_1\right].\end{array}$

3.4. 证明良好平衡性质

所有上述程序一起组成了沿沟渠的浅水方程的有限差分AWENO格式，我们可以使用修正后的数值通量${\widehat{F}}_{i+\frac{1}{2}}^{\pm }$ 和源项${\widehat{S}}_{i+\frac{1}{2}}^{\pm }$ 来重写数值格式(3.1)，重写后的格式和命题如下所述：

$\frac{\partial U}{\partial t}+\frac{{\tilde{F}}_{i+\frac{1}{2}}^{-}-{\tilde{F}}_{i-\frac{1}{2}}^{+}}{\Delta x}=0,{\tilde{F}}_{i+\frac{1}{2}}^{\pm }={\widehat{F}}_{i+\frac{1}{2}}^{\pm }-{\widehat{S}}_{i+\frac{1}{2}}^{\pm }.$

命题1. 上述沿沟渠的浅水方程有限差分AWENO格式对于定常解来说是良好平衡的。

证明

$\begin{array}{l}{\widehat{F}}_{i+\frac{1}{2}}^{\pm }=\left(\begin{array}{c}0\\ \frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{*,\pm }{\left(h_{i+\frac{1}{2}}^{*,\pm }\right)}^2\end{array}\right)+\left(\begin{array}{c}0\\ \frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2-\frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{*,\pm }{\left(h_{i+\frac{1}{2}}^{*,\pm }\right)}^2\end{array}\right)+{\mathcal{D}}_{i+\frac{1}{2}}\left[F\right]\\ =\left(\begin{array}{c}0\\ \frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2\end{array}\right)+{\mathcal{D}}_{i+\frac{1}{2}}\left[F\right],\end{array}$

边界值${\widehat{S}}_{i+\frac{1}{2}}^{\pm }$ 为：

$\begin{array}{l}{\widehat{S}}_{i+\frac{1}{2}}^{\pm }=\frac{1}{2}g{\left(h_i+b_i\right)}^2{\left({\widehat{S}}_1\right)}_{i+\frac{1}{2}}^{\pm }-g\left(h_i+b_i\right){\left(\widehat{B°S_1}\right)}_{i+\frac{1}{2}}^{\pm }+\frac{1}{2}g{\left(\widehat{B°B°S_1}\right)}_{i+\frac{1}{2}}^{\pm }\\ =\frac{1}{2}g{\left(h_i+b_i\right)}^2\{(\begin{array}{c}0\\ {\sigma }_{i+\frac{1}{2}}^{\pm }\end{array})+{\mathcal{D}}_{i+\frac{1}{2}}\left[S_1\right]\}-g\left(h_i+b_i\right)\left\{\left(\begin{array}{c}0\\ {\sigma }_{i+\frac{1}{2}}^{\pm }{b}_{i+\frac{1}{2}}^{\pm }\end{array}\right)+{\mathcal{D}}_{i+\frac{1}{2}}\left[B°S_1\right]\right\}\\ +\frac{1}{2}g\left\{\left(\begin{array}{c}0\\ {\sigma }_{i+\frac{1}{2}}^{\pm }{\left(b_{i+\frac{1}{2}}^{\pm }\right)}^2\end{array}\right)+{\mathcal{D}}_{i+\frac{1}{2}}\left[B°B°S_1\right]\right\}.\end{array}$

由于算子${\mathcal{D}}_{i+\frac{1}{2}}[\cdot ]$ 是线性的，所以我们可以得到：

$\begin{array}{l}{\mathcal{D}}_{i+\frac{1}{2}}\left[F\right]-\left\{\frac{1}{2}g{\left(h_i+b_i\right)}^2{\mathcal{D}}_{i+\frac{1}{2}}\left[S_1\right]-g\left(h_i+b_i\right){\mathcal{D}}_{i+\frac{1}{2}}\left[B\circ S_1\right]+\frac{1}{2}g{\mathcal{D}}_{i+\frac{1}{2}}\left[B\circ B\circ S_1\right]\right\}\\ ={\mathcal{D}}_{i+\frac{1}{2}}\left[F-\frac{1}{2}g{\left(h_i+b_i\right)}^2{S}_1+g\left(h_i+b_i\right)\left(B\circ S_1\right)-\frac{1}{2}g\left(B\circ B\circ S_1\right)\right]\\ ={\mathcal{D}}_{i+\frac{1}{2}}\left[\left(\begin{array}{c}0\\ 0\end{array}\right)\right]=0.\end{array}$

然后，${\tilde{F}}_{i+\frac{1}{2}}^{\pm }$ 简化为：

$\begin{array}{l}{\tilde{F}}_{i+\frac{1}{2}}^{\pm }={\widehat{F}}_{i+\frac{1}{2}}^{\pm }-{\widehat{S}}_{i+\frac{1}{2}}^{\pm }\\ =\left(\begin{array}{c}0\\ \frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2\end{array}\right)-\frac{1}{2}g\left\{{\left(h_i+b_i\right)}^2\left(\begin{array}{c}0\\ {\sigma }_{i+\frac{1}{2}}^{\pm }\end{array}\right)-2\left(h_i+b_i\right)\left(\begin{array}{c}0\\ {\sigma }_{i+\frac{1}{2}}^{\pm }{b}_{i+\frac{1}{2}}^{\pm }\end{array}\right)+\left(\begin{array}{c}0\\ {\sigma }_{i+\frac{1}{2}}^{\pm }{\left(b_{i+\frac{1}{2}}^{\pm }\right)}^2\end{array}\right)\right\}\\ =\left(\begin{array}{c}0\\ \frac{1}{2}g{\sigma }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2-\frac{1}{2}g{\left(h_i+b_i\right)}^2\left({\sigma }_{i+\frac{1}{2}}^{\pm }\right)+g\left(h_i+b_i\right)\left({\sigma }_{i+\frac{1}{2}}^{\pm }{b}_{i+\frac{1}{2}}^{\pm }\right)-\frac{1}{2}g\left({\sigma }_{i+\frac{1}{2}}^{\pm }{\left(b_{i+\frac{1}{2}}^{\pm }\right)}^2\right)\end{array}\right)\\ =\left(\begin{array}{c}0\\ 0\end{array}\right).\end{array}$

因此，有${\tilde{F}}_{i+\frac{1}{2}}^{-}={\tilde{F}}_{i+\frac{1}{2}}^{+}$ ，这意味着格式是良好平衡的。

3.5. 平衡格式的总结

最后，我们总结求解浅水沟渠方程(1.1)的良好平衡有限差分AWENO格式的完整过程。AWENO格式详细的算法如下所示：

重新表述源项，并将方程改写成如下形式(3.1)。

应用AWENO程序构造数值通量${\widehat{F}}_{i+\frac{1}{2}}^{\pm }$ 。

近似改写后源项中的导数项${\sigma }_x,{\left(\sigma b\right)}_x$ 和${\left(\sigma b^2\right)}_x$ 。

将数值通量的余数与源项近似相加，采用龙格–库塔方法即时推进时间。然后重复步骤2~4。

4. 数值结果

在本节中，我们通过大量的基准数值算例验证了所提出的浅水沿河道流动的良好平衡AWENO格式。在这些例子中实现了五阶有限差分AWENO格式(即k = 2)与三阶TVD龙格–库塔方法[11]的耦合。CFL值取0.6。重力常数g取值为9.812。具有连续和不连续宽度函数的沟渠已经进行了测试。

4.1. 精度测试

我们首先在一个光滑解的例子上测试所提格式的5阶精度。周期性底部地形与沟渠宽度函数如下：

$b\left(x\right)={\sin }^2\left(\pi x\right),\sigma \left(x\right)={\text{e}}^{\sin \left(2\pi x\right)},$

在本例中考虑。初始条件如下

$h\left(x,0\right)=3+{\text{e}}^{\cos \left(2\pi x\right)},Q\left(x,0\right)=\sin \left(\cos \left(2\pi x\right)\right),x\in \left[0,1\right],$

具有周期边界条件。我们运行测试直到停止时间$t=0.1$ ，此时解仍然是光滑的。

由于该非线性系统没有精确解，我们采用同样的格式，将N = 25,600个单元格细化为参考解，然后在计算误差和收敛率时将其作为精确解，能够发现达到预期的5阶精度。

4.2. 良好平衡测试

第二个测试问题[12] [13]的目的是验证当前AWENO格式的良好平衡性。我们考虑给出的底部地形：

$b\left(x\right)=\{\begin{array}{l}0.25\left(1+\cos \left(10\pi \left(x-0.5\right)\right)\right),\text{if}0.4\le x\le 0.6,\hfill \\ 0,\text{otherwise}.\hfill \end{array}$ (4.1)

在区域[0, 1]中。随宽度$\sigma \left(x\right)$ 变化的沟渠形式为：

$\sigma \left(x\right)=\{\begin{array}{l}1-{\sigma }_0\left(1+\cos \left(2\pi \frac{x-\left(x_l+x_r\right)/2}{x_r-x_l}\right)\right),\text{if}x\in \left[x_l,x_r\right],\hfill \\ 1,\text{otherwise}.\hfill \end{array}$ (4.2)

其中$x_l$ 和$x_r$ 分别为收缩的左右边界，$1-2{\sigma }_0$ 表示在$\left(x_l+x_r\right)/2$ 点处的最小沟渠宽度。在本例中，我们选取$x_l=0.25,x_r=0.75$ 和${\sigma }_0=0.2$ 。给出了稳态解的初始条件：

$h+b=1,Q=\sigma hu=0.$

并施加周期边界条件。

我们在一个有200个均匀单元的网格上到$t=1$ 时刻解决这个问题。数值表面水平$h+b$ 和底部b如图1所示。

4.3. 小扰动测试

在本节中，我们模拟小扰动到稳态解的传播，以证明所提出的AWENO格式在这种具有挑战性的情况下的能力，该测试最早由[12]提出。设底部地形为(4.1)，初始条件为：

$h+b=\{\begin{array}{l}1+0.01,\text{if}0.1\le x\le 0.2,\hfill \\ 1,\text{otherwise},\hfill \end{array}Q=\sigma hu=0,$

在具有简单传输边界条件的计算区域[0, 1]中。对(4.2)中定义的两组沟渠$\sigma \left(x\right)$ 进行了检验，其中一组沟渠左移收缩$x_l=0.15,x_r=0.65,{\sigma }_0=0.2$ ，另一组沟渠右移收缩$x_l=0.\text{3}5,x_r=0.\text{8}5,{\sigma }_0=0.2$ 。对于这些

(a) (b)

Figure 1. Numerical results of the example in Section 4.2, the vertical axis represents the width of the channel, surface elevation, and bottom topography respectively, while the horizontal axis represents the position

图1. 在4.2节例子的数值结果，纵坐标分别表示沟渠宽度，表面水平和底部地形，横坐标表示位置

涉及如此小的稳态解扰动的测试，非平衡数值格式通常难以计算并导致振荡[14]。在图2中，我们展示了沟渠的三维视图和地形的垂直剖面。所提出的良好平衡有限差分AWENO格式在200个均匀计算单元上不同时间的数值结果，与精化的2000个单元“参考”格式的数值结果对比如图3所示。为了比较，我们还给出了200单元的非平衡AWENO格式的数值结果。我们可以清楚地看到，数值结果没有虚假的数值振荡；我们的格式可以在数值上，在相对粗糙的网格上很好地捕捉到如此小的扰动，而非平衡型则无法捕捉到小扰动。

(a) (b)

(c)

Figure 2. Numerical results of the example in Section 4.3, the vertical axis represents the bottom terrain, and the horizontal axis represents the position

图2. 在4.3节例子的数值结果，纵坐标表示底部地形，横坐标表示位置

(a) (b)

(c) (d)

(e) (f)

(g) (h)

(i) (j)

Figure 3. Numerical results of the example in Section 4.3, the vertical axis represents the width of the channel and the surface level, while the horizontal axis represents the position

图3. 在4.3节例子的数值结果，纵坐标分别表示沟渠宽度和表面水平，横坐标表示位置

4.4. 一个收敛发散的沟渠

这里我们考虑的是一个收敛发散沟渠中的经典跨临界定常流动，最初由García-Navarro等人在[15]中提出。该试验涉及许多实际问题，如桥墩之间的流动。底部设为平面(即b = 0)，收敛发散沟渠为：

$\sigma \left(x\right)=\{\begin{array}{l}5-0.7065\left(1+\cos \left(2\pi \frac{x-250}{300}\right)\right),\text{if}\left|x-250\right|\le 150,\hfill \\ 5,\text{otherwise}.\hfill \end{array}$

如图4所示，在计算区域[0, 500]中。初始条件为：

$h=2,Q=\sigma hu=20.$

边界条件为上游$Q=20$ ，下游$h=1.85$ 。宽度的变化不仅影响传播锋面的轮廓，也影响稳态流的轮廓。另外，由于边界条件，在亚临界流和超临界流之间出现了一个水力跃变。我们使用200个均匀单元计算这个测试直到t = 5000 (直到达到稳定状态)。我们在图5中给出了水高h和弗劳德数与参考解比较的数值结果。从图5中，我们观察到水在接近最大收缩点时加速(在我们的例子中σ(250) = 3.587)，并且流动在那里发展为临界。然后从亚临界流过渡到超临界流，形成一个静止的水力跃变，与下游条件所需的亚临界剖面连接。数值结果与参考解吻合较好，与文献[15] [16]吻合较好。

Figure 4. Numerical results of the example in Section 4.4, the vertical axis represents the width of the ditch, while the horizontal axis represents the position

图4. 在4.4节例子的数值结果，纵坐标表示沟渠宽度，横坐标表示位置

(a) (b)

Figure 5. Numerical results of the example in Section 4.4, the vertical axis represents surface elevation and water flow velocity respectively, while the horizontal axis represents position

图5. 在4.4节例子的数值结果，纵坐标分别表示表面水平和水流速度，横坐标表示位置

4.5. 驼峰上移动水的稳定状态

在最后一个数值例子[13]中，我们考虑了提出的AWENO格式对不同沟渠的稳定跨临界和亚临界流动的收敛性。在[17]中，使用相同的试验来检测浅水方程(即恒定沟渠宽度)的良好平衡AWENO格式的性能。计算区域设为[0, 25]，底部地形定义为

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

我们选择不同的可变沟渠宽度，在每种情况下定义，以展示沟渠对最终解决格式的影响。初始条件为

$h\left(x,0\right)=0.5-b\left(x\right),Q\left(x,0\right)=0.$

我们取200个均匀的计算单元，并将最终时间设为t = 200，当水流达到静水稳态时。根据不同的边界条件，流量可以是亚临界或跨临界的，有或没有稳定激波。为了便于比较，我们还提供了解析解，其计算方法如[15]所示。具体来说，为了得到沟渠中每个网格点的解析解，我们求解了一个三次方程，该方程来源于每个分支的水头的强制保护。此外，激波两侧的两个剖面通过跳跃条件连接。

1) 亚临界流

边界条件上游设为$hu=\text{4}\text{.42}$ ，下游设为$h=2$ 。在这个测试用例中，考虑了(4.2)中定义的两个不同的沟渠宽度$\sigma \left(x\right)$ 集，其中一个具有左移收缩$x_l=3.75,x_r=16.25,{\sigma }_0=0.05$ ，另一个具有右移收缩$x_l=8.75,x_r=21.25,{\sigma }_0=0.05$ 。在$t=200$ 时，水流将演化为流动的水稳态。聚合状态是亚临界流。我们在图6中给出了$t=200$ 时的表面水平$h+b$ 和质量流率Q，并将其中的解析解进行了比较。很明显，数值解和分析情况很一致。我们还可以观察到不同的$\sigma \left(x\right)$ 对收敛稳态解的影响。

2) 无激波跨临界流动。

边界条件设上游为$hu=\text{1}\text{.53}$ ，下游为$h=0.66$ 。考虑在(4.2)中定义的两个不同的沟渠集$\sigma \left(x\right)$ ，一个是左移收缩$x_l=3.75,x_r=16.25,{\sigma }_0=0.15$ ，另一个是右移收缩$x_l=8.75,x_r=21.25,{\sigma }_0=0.15$ 。收敛的稳态解是无激波的跨临界流。表面水平$h+b$ 和质量流量Q如图7所示，与解析解吻合较好。

3) 带激波的跨临界流动。

(a) (b)

(c) (d)

Figure 6. Numerical results of the example in Section 4.5, the vertical axis represents surface elevation and mass flow rate Q respectively, while the horizontal axis represents position

图6. 在4.5节例子的数值结果，纵坐标分别表示表面水平和质量流量Q，纵坐标表示位置

边界条件设为上游$hu=0.18$ ，下游$h=0.33$ 。考虑在(4.2)中定义的两个不同的沟渠集$\sigma \left(x\right)$ ，一个是左移收缩$x_l=3.75,x_r=16.25,{\sigma }_0=0.15$ ，另一个是右移收缩$x_l=8.75,x_r=21.25,{\sigma }_0=0.15$ 。收敛的稳态解是一个中间出现激波的跨临界流。我们在图8中给出了表面水平$h+b$ 和质量流量Q，它们与解析解吻合得很好。

(a) (b)

(c) (d)

Figure 7. Numerical results of the example in Section 4.5, the vertical axis represents surface elevation and mass flow rate Q respectively, while the horizontal axis represents position

图7. 在4.5节例子的数值结果，纵坐标分别表示表面水平和质量流量Q，纵坐标表示位置

(a) (b)

(c) (d)

Figure 8. Numerical results of the example in Section 4.5, the vertical axis represents surface elevation and mass flow rate Q respectively, while the horizontal axis represents position

图8. 在4.5节例子的数值结果，纵坐标分别表示表面水平和质量流量Q，纵坐标表示位置

5. 结论

在本文中，我们针对具有不规则几何形状和非平坦底部地形的开放沟渠浅水方程，提出了高阶有限差分AWENO格式。我们通过静水重构等方法近似数值通量，通过修改源项方法离散源项，提出的格式保持了良好平衡特性。即在通量梯度和源项之间存在精确平衡时，在离散状态保持稳态。严格的理论分析和大量的数值实验都证明了本文提出的高阶有限差分AWENO格式具有良好的均衡性，对解具有较高的分辨率，在光滑区域具有较高的精度阶。

基金项目

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

NOTES

^*通讯作者。