多孔浅水方程的高阶平衡保正有限差分AWENO格式

High-Order Well-Balanced and Positivity-Preserving Finite-Difference AWENO Scheme with Hydrostatic Reconstruction for Porous Shallow Water Equations

摘要: 本文提出了求解具有不连续孔隙度和底部地形的一维多孔浅水方程的高阶平衡保正有限差分AWENO格式，所提出的格式保持了静水稳态的良好平衡特性。在这个数值框架中，采用静水重构方法具有两个主要优点：1) 使用任意单调通量和重新表述源项的方法获得了良好平衡特性。2) 采用Lax-Friedrichs (LF)通量的一阶方案在适当的时间步长内保持了水高保正性。通过大量的数值算例验证了该格式具有高阶精度和良好平衡特性，所有算例的数值结果与解析解一致。

Abstract: In this paper, we propose a higher-order well-balanced and positivity-preserving finite-difference AWENO scheme for solving the one-dimensional porous shallow water equation with discontinuous porosity and bottom topography. The proposed format maintains the well-balanced property of the hydrostatic steady state. In this numerical framework, the hydrostatic reconstruction (HR) method is employed with two main advantages: 1) The method using arbitrary monotone fluxes and reformulated source terms obtains the well-balanced property. 2) The first-order scheme using Lax-Friedrichs (LF) fluxes and the HR method maintains the water height preserving properties with an appropriate time step. We verify that the scheme has high-order accuracy and well-balanced properties through a large number of numerical examples, and the numerical results of all cases agree with the analytical solutions.

1. 引言

多孔浅水方程已被数值离散化，用于模拟城市地区的淹没[1]-[5]和植被地区的流动[6]-[10]等问题。忽略能量损失的一维多孔浅水方程的质量和动量表示为

$\begin{array}{l}\frac{\partial \left(\varphi h\right)}{\partial t}+\frac{\partial \left(\varphi hu\right)}{\partial x}=0,\\ \frac{\partial \left(\varphi hu\right)}{\partial t}+\frac{\partial \left(\varphi h{u}^2+\frac{1}{2}g\varphi h^2\right)}{\partial x}=\frac{1}{2}g{h}^2\frac{\partial \varphi }{\partial x}-g\varphi h\frac{\partial b}{\partial x},\end{array}$ (1.1)

其中，t表示时间，x为水流坐标。$h=h\left(x,t\right)$ 表示水深，$u=u\left(x,t\right)$ 为沿x方向的水平速度，$b=b\left(x\right)$ 代表底部地形，g是重力加速度，$\varphi =\varphi \left(x\right)$ 是通道横截面孔隙度且$\varphi \in \left[0,1\right]$ 。当$\varphi =1$ 时，孔隙度最大，水流完全通畅，这时就变为传统浅水方程

$\begin{array}{l}\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}=0,\\ \frac{\partial hu}{\partial t}+\frac{\partial \left(h{u}^2+\frac{1}{2}g{h}^2\right)}{\partial x}=-gh\frac{\partial b}{\partial x}.\end{array}$ (1.2)

方程(1.1)是带源项的双曲守恒律，满足静水定常解

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

求解这个方程格式的一个主要困难是：在离散水平上标准数值格式破坏了静水定常解下的通量梯度和源项的精确平衡，并且还会导致定常解附近出现伪振荡。此方程已经通过多种数值方法进行了解决[11] [12]。代表性的研究工作包括：间断Galerkin (RKDG)方法[13]-[15]，谱元方法[16]，无单元有限元方法[17]。本文针对具有不连续孔隙度和底部地形的一维多孔浅水方程开发了有限差分AWENO格式。在过去的几十年里，已经开发出很多高阶精确的数值格式，包括有限差分/体积加权基本非振荡(WENO)格式[11] [18]-[20]。有限差分WENO格式具有高分辨率，高精度，适应性强等优点[21] [22]，但是在静水定常解的条件下有数值振荡。因此，在有限差分WENO格式的基础上，我们开发有限差分AWENO格式，该格式具有以下优点：1) 实现简单，2) 计算效率高，3) 容易扩展多维模型，4) 单元界面处计算数值通量灵活性高[23]。在WENO格式构造数值通量和源项完成后，我们会加入修正项来保证格式不会出现数值振荡，参考[15]我们改写原方程，并把源项重新改写为

$U_t+F\left(U,\varphi \right)=S\left(U,\varphi ,\varphi b,\varphi b^2\right),$ (1.4)

其中

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

其中，质量流量$Q=hu$ ，$H=\phi h$ 为湿截面。

2. 高阶有限差分AWENO格式

简要回顾求解双曲型守恒律的五阶保守有限差分AWENO格式的框架。在不失一般性的前提下，考虑一维标量的情况

$\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial t}=0.$ (2.1)

考虑由单元中心$x_i=a+\left(i-\frac{1}{2}\right)\Delta x$ 定义的均匀网格，单元边界点由$x_{i\pm \frac{1}{2}}=x_i\pm \frac{1}{2}\Delta x$ 得到，其中$\Delta x=\frac{b-a}{N}$ 为均匀网格间距，计算区域为$\left[a,b\right]$ ，考虑的半离散形式

${\frac{\partial u}{\partial t}|}_{x=x_i}+\frac{{\widehat{f}}_{i+\frac{1}{2}}-{\widehat{f}}_{i-\frac{1}{2}}}{\partial t}=0,i=1,\cdots ,N,$

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

2.1. 五阶WENO重构

使用三阶插值多项式$q^{\left(k\right)}\left(x\right)$ 和非线性权值${\omega }_k$ 通过全局模板$S^5={\displaystyle \cup {}_{k=0}^2}$ 的每个子模板 $S_k=\left[x_{i-2+k},\cdots ,x_{i+k}\right],k=0,1,2$ 得到基于点$x_{i+\frac{1}{2}}$ 的单元左侧边界值$u_{i+\frac{1}{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}$ ，其中，$u_{i+\frac{1}{2}}^{\left(k\right)}=u^{\left(k\right)}\left(x_{i+\frac{1}{2}}\right)$ ，改进的非线性权定义为${\omega }_k=\frac{{\alpha }_k}{{\displaystyle \sum _{s=0}^2{\alpha }_s}}$ ，${\alpha }_k$ 可由以下公式得到${\alpha }_k=d_k\left(1+{\left(\frac{{\tau }_5}{{\beta }_k+ϵ}\right)}^p\right)$ ，其中，$d_0=\frac{1}{16}$ ，$d_1=\frac{5}{8}$ ，$d_2=\frac{5}{16}$ 是理想的线性权，在数值算例中取$p=2$ ，$ϵ={10}^{-12}$ 。取${\beta }_k$ 为${\beta }_k={\displaystyle \sum _{l=1}^2\Delta }{x}^{2l-1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}^{x_{i+\frac{1}{2}}}}\left(\frac{{\partial }^l}{\partial x^l}{u}^{\left(k\right)}\left(x\right)\right)}$ ，取${\tau }_5=\left|{\beta }_0-{\beta }_2\right|$ 。

以上讨论中，由于全局模板偏向目标点$x_{i+\frac{1}{2}}$ 左侧，因此单元左侧边界值$u_{i+\frac{1}{2}}^{-}$ 已经被插值求出。右单元边值$u_{i+\frac{1}{2}}^{+}$ 与左单元边界值选用模板相对于目标点$x_{i+\frac{1}{2}}$ 镜像对称。

2.2. 构造数值通量${\widehat{f}}_{i+\frac{1}{2}}$

对于五阶AWENO格式，数值通量表示为

${\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],$ (2.2.1)

其中$ℱ\left(\cdot ,\cdot \right)$ 表示单调通量，${\mathcal{D}}_{i+\frac{1}{2}}\left[f\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}},$ (2.2.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.4)，我们可以对方程的齐次形式使用AWENO格式

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

首先，通过雅可比矩阵$\frac{\partial F\left(U\right)}{\partial U}$ ，得到平均左特征向量$L_{i+\frac{1}{2}}$ ，然后将$U_{i+l}$ 投影到特征空间${\tilde{U}}_{i+l}=L_{i+\frac{1}{2}}{U}_{i+l}$ ，$l=-2,\cdots ,3$ 。 然后，通过五阶AWENO插值计算边界值${\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)$ 。最后，通过右特征向量$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}}^{+}$ 通过相似的过程得到，并且使用的模板${\tilde{U}}_{i+\frac{1}{2}}^{-}$ 和${\tilde{U}}_{i+\frac{1}{2}}^{+}$ 关于目标点$x_{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].$ (2.3.2)

注1. 我们将这个五阶WENO插值过程记作$U_{i+\frac{1}{2}}^{\pm }={ℒ}_{i+\frac{1}{2}}^{\pm }\left[U\right]$ ，同时，算子$ℒ[\cdot ]$ 是一个线性算子。修正项${\mathcal{D}}_{i+\frac{1}{2}}\left[F\right]$ 通过得到。

3. 平衡AWENO格式的构造

多孔浅水方程(1.1)的保守型特有限差分AWENO的半离散形式表示为

${\frac{\text{d}U}{\text{d}t}|}_{x=x_j}+\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.1)

定义新的向量$B={\left(0,b\right)}^{\top }$ ，$S_1={\left(0,\varphi \right)}^{\top }$ 同样地，向量$B$ 和$S_1$ 通过AWENO插值重构得到

$\begin{array}{l}{B}_{i+\frac{1}{2}}^{\pm }=P^{-1}{ℒ}_{i+\frac{1}{2}}^{\pm }\left[PB\right],\\ S_1{}_{i+\frac{1}{2}}^{\pm }=P^{-1}{ℒ}_{i+\frac{1}{2}}^{\pm }\left[P{S}_1\right],\end{array}$ (3.2)

其中，线性WENO插值算子${ℒ}_{i+\frac{1}{2}}^{\pm }\left[PB\right]$ 使用与${ℒ}_{i+\frac{1}{2}}^{\pm }\left[U\right]$ 相同的非线性权值，且$P$ 是置换矩阵$P=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$ 。

定义向量$Q={\left(h,hu\right)}^{\top }$ ，根据上述插值过程，显然有

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

使用HR方法，我们设$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)$ ，${\varphi }_{i+\frac{1}{2}}^{*,\pm }=\max \left(0,\max \left({\varphi }_{i+\frac{1}{2}}^{+},{\varphi }_{i+\frac{1}{2}}^{-}\right)\right)$ ，修改左右插值$U_{i+\frac{1}{2}}^{\pm }$ 为

$U_{i+\frac{1}{2}}^{*,\pm }=(\begin{array}{c}{\varphi }_{i+\frac{1}{2}}^{*,\pm }{h}_{i+\frac{1}{2}}^{*,\pm }\\ \\ {\varphi }_{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],$ (3.4)

加上静水修正项 ${\mathscr{H}}_{i+\frac{1}{2}}^{\pm }=(\begin{array}{c}0\\ \frac{1}{2}g{\varphi }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2-\frac{1}{2}g{\varphi }_{i+\frac{1}{2}}^{*,\pm }{\left(h_{i+\frac{1}{2}}^{*,\pm }\right)}^2\end{array}).$

3.1. 源项离散

我们把方程(1.4)的源项重新写成

$S\left(U,\varphi ,\varphi b,\varphi 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},$ (3.1.1)

其中符号$\circ$ 定义为对应分量相乘，与构造数值通量一致的方法，我们构造源项

${\widehat{S}}_{i+\frac{1}{2}}^{\pm }=\frac{1}{2}g{(h+b)}^2{\widehat{S_1}}_{i+\frac{1}{2}}^{\pm }-g\left(h+b\right){\left(\widehat{B\circ S_1}\right)}_{i+\frac{1}{2}}^{\pm }+\frac{1}{2}g{\left(\widehat{B\circ B\circ S_1}\right)}_{i+\frac{1}{2}}^{\pm },$ (3.1.2)

分别定义${\widehat{S_1}}_{i+\frac{1}{2}}^{\pm }$ ，${\left(\widehat{B\circ S_1}\right)}_{i+\frac{1}{2}}^{\pm }$ 和${\left(\widehat{B\circ B\circ S_1}\right)}_{i+\frac{1}{2}}^{\pm }$ 为

$\begin{array}{l}{\widehat{S_1}}_{i+\frac{1}{2}}^{\pm }=S_1{}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[S_1\right],\\ {\left(\widehat{B\circ S_1}\right)}_{i+\frac{1}{2}}^{\pm }=B_{i+\frac{1}{2}}^{\pm }\circ S_1{}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[B\circ S_1\right],\\ {\left(\widehat{B\circ B\circ S_1}\right)}_{i+\frac{1}{2}}^{\pm }=B_{i+\frac{1}{2}}^{\pm }\circ B_{i+\frac{1}{2}}^{\pm }\circ S_1{}_{i+\frac{1}{2}}^{\pm }+{\mathcal{D}}_{i+\frac{1}{2}}\left[B\circ B\circ S_1\right],\end{array}$ (3.1.3)

其中，$S_1{}_{i+\frac{1}{2}}^{\pm }$ 和$B_{i+\frac{1}{2}}^{\pm }$ 通过(3.3)得到。

3.2. 平衡格式的证明

我们使用修正过的${\widehat{F}}_{i+\frac{1}{2}}^{\pm }$ 和${\widehat{S}}_{i+\frac{1}{2}}^{\pm }$ 把通量改写为

${\tilde{F}}_{i+\frac{1}{2}}^{\pm }={\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.$ (3.2.1)

命题1 对于方程(1.1),半离散格式(3.2.1)是五阶光滑，并对于静水定常解(1.3)是平衡的。

证明 半离散格式五阶光滑解是命题1的推论。

${\widehat{F}}_{i+\frac{1}{2}}^{\pm }=\left(\begin{array}{c}0\\ \frac{1}{2}g{\varphi }_{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{\varphi }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2-\frac{1}{2}g{\varphi }_{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],$

然后${\widehat{S}}_{i+\frac{1}{2}}^{\pm }$ 写成

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

由于$D_{i+\frac{1}{2}}\left[F\right]$ 是线性算子，这里使用的${\varphi }_{i+\frac{1}{2}}^{\pm }$ 和$b_{i+\frac{1}{2}}^{\pm }$ 都是通过AWENO重构得到的，所以我们可以得到

$\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[B\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}B+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\\ \frac{1}{2}g\varphi h^2-\frac{1}{2}g{\left(h+b\right)}^2\varphi +g\left(h+b\right)\varphi b-\frac{1}{2}g\varphi b^2\end{array}\right)\right]=0.\end{array}$

在静水状态下，${\tilde{F}}_{i+\frac{1}{2}}^{\pm }$ 则可以化简为

$\begin{array}{c}{\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{\varphi }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2\end{array}\right)-\left\{-g\left(h_i+b_i\right)\left(\begin{array}{c}0\\ b_{i+\frac{1}{2}}^{\pm }{\varphi }_{i+\frac{1}{2}}^{\pm }\end{array}\right)+\frac{1}{2}g\left(\begin{array}{c}0\\ {\left(b_{i+\frac{1}{2}}^{\pm }\right)}^2{\varphi }_{i+\frac{1}{2}}^{\pm }\end{array}\right)+\frac{1}{2}g{\left(h_i+b_i\right)}^2\left(\begin{array}{c}0\\ {\varphi }_{i+\frac{1}{2}}^{\pm }\end{array}\right)\right\}\\ =\left(\begin{array}{c}0\\ \frac{1}{2}g{\varphi }_{i+\frac{1}{2}}^{\pm }{\left(h_{i+\frac{1}{2}}^{\pm }\right)}^2+g\left(h_i+b_i\right)\left(b_{i+\frac{1}{2}}^{\pm }{\varphi }_{i+\frac{1}{2}}^{\pm }\right)-\frac{1}{2}g{\varphi }_{i+\frac{1}{2}}^{\pm }{\left(b_{i+\frac{1}{2}}^{\pm }\right)}^2-\frac{1}{2}g{\left(h_i+b_i\right)}^2{\varphi }_{i+\frac{1}{2}}^{\pm }\end{array}\right)=0.\end{array}$

因此，格式(3.2.1)是平衡的。

4. 数值结果

在本节中，我们将通过计算不同算例来测试所提出的有限差分AWENO格式的性能。对于时间离散，我们采用经典的三阶Runge-Kutta方法。

4.1. 小扰动测试

以[24]引入的扰动波传播为例，检验这种格式是否可以解决小脉冲在不规则地形上的传播问题，地形与孔隙度函数为

$b\left(x\right)=\{\begin{array}{ll}0.025\left(1+\cos \left(10\pi \left(x-0.5\right)\right)\right),\hfill & \text{if}0.4\le x\le 0.6,\hfill \\ 0,\hfill & \text{otherwise,}\hfill \end{array}$ (4.1.1)

和

$\varphi \left(x\right)=\{\begin{array}{ll}1-{\varphi }_0\left(1+\cos \left(2\pi \frac{x-\left(x_r+x_l\right)/2}{x_r-x_l}\right)\right),\hfill & \text{if}0.4\le x\le 0.6,\hfill \\ 0,\hfill & \text{otherwise, }\hfill \end{array}$ (4.1.2)

其中，$x_l$ 和$x_r$ 分别表示孔隙度收缩的左右边界，${\varphi }_0$ 为决定孔隙度大小的形状参数，见图1。

由于孔隙度函数在$\left(x_l+x_r\right)/2$ 处孔隙度为$1-2{\varphi }_0$ 的最小值，因此${\varphi }_0$ 应小于$0.5$ ，初始状态由下式给出

$h\left(x,t=0\right)=\{\begin{array}{ll}1+0.01-b\left(x\right),\hfill & \text{if}0.1\le x\le 0.2,\hfill \\ 1-b\left(x\right),\hfill & \text{otherwise, }\hfill \end{array}\text{ and }u\left(x,t=0\right)=0,$

计算区域为$x\in \left[0,1\right]$ ，在区域的东西边界施加了零梯度的流出条件。一个箱形扰动波分裂成两个脉冲

Figure 1. Numerical results of the example in Section 4.1, the vertical axis represents the bottom terrain and the porosity, while the horizontal axis represents the position

图1. 4.1节算例的底部地形与孔隙度的初始条件，纵坐标分别表示底部地形和孔隙度，横坐标表示位置

波。一个向西移动，穿过西部边界离开区域，而另一个向东传播，与地形相互作用。与之前的研究相同[25] [26]，考虑了两种不同孔隙度$\varphi \left(x\right)$ 集合的情况：(i) 情况1，$\left(x_l,x_r,{\varphi }_0\right)=\left(0.15,0.65,0.20\right)$ ，(ii) 情况2，$\left(x_l,x_r,{\varphi }_0\right)=\left(0.35,0.85,0.20\right)$ 。分别用200和2000个均匀网格进行数值试验，结果较为吻合(图2)。这意味着通量和源项小到足以分辨微小脉冲的演化。因此，我们可以得出结论，该格式能够解决流体接近稳态的问题。

4.2. 经过驼峰的稳定水流

驼峰定常流动能够验证了数值格式对亚临界流动、无激波跨临界流动和有激波跨临界流动定常解的收敛性[27]-[29]，是浅水求解器的经典基准试验。为了测试求解Saint-Venant方程的数值格式，根据通道

Figure 2. Numerical results of the example in Section 4.1, the vertical axis represents the width of the surface level, while the horizontal axis represents the position

图2. 4.1节算例的数值结果，纵坐标表示表面水平，横坐标表示位置

宽度变化对经典格式进行了修正。按照[25]孔隙度和底部地形用式(4.1.1)和(4.1.2)表示，在计算区域$x\in \left[0,25\right]$ 中，考虑两组不同的孔隙度函数参数：(i) 一组为$\left(x_l,x_r,{\varphi }_0\right)=\left(3.75,16.25,{\varphi }_0\right)$ ，(ii) 另一组为$\left(x_l,x_r,{\varphi }_0\right)=\left(8.75,21.25,{\varphi }_0\right)$ 。初始条件取$h\left(x,t=0\right)=0.5-b\left(x\right)$ and $u\left(x,t=0\right)=0$ ,不同情况相关的边界条件如下

(1) 亚临界流：上游$\varphi hu=4.42$ ，下游$h=2$ ；

(2) 激波过渡流：上游$\varphi hu=0.18$ ，下游$h=0.33$ ；

(3) 无激波过渡流：上游$\varphi hu=1.53$ ，下游$h=0.66$ 。

为了得到相应的流动条件，在亚临界流动情况下，形状参数${\varphi }_0$ 设置为0.05，在有激波和无激波的过渡流动情况下，形状参数${\varphi }_0$ 设为0.15。在所有情况下，使用200个均匀电池和不同的CFL条件进行$t=200\text{s}$ 的数值模拟，结果与解析解吻合良好(图3~5)。

5. 结论

在本文中，我们提出了具有不连续孔隙度和底部地形的一维多孔浅水方程的高阶有限差分AWENO格式。所提出的格式保持了静水稳态的良好平衡特性，即在通量梯度和源项之间存在精确平衡时，在离散状态保持稳态。通过一种新颖的源项分裂和一种通量修正技术的良好平衡数值通量以及对每个分裂源项的良好平衡数值近似来实现良好平衡特性。大量的数值实验表明，本文提出的AWENO格式具有良好平衡特性，对不连续解具有较高分辨率，在光滑区域具有较高精度。

Figure 3. Numerical results of the example in Section 4.2, the vertical axis represents the width of the surface level and miss flow rate Q respectively, while the horizontal axis represents the position

图3. 4.2节算例的数值结果，纵坐标表示表面水平和质量流量Q，横坐标表示位置

Figure 4. Numerical results of the example in Section 4.2, the vertical axis represents the width of the surface level and miss flow rate Q respectively, while the horizontal axis represents the position

图4. 4.2节算例的数值结果，纵坐标表示表面水平和质量流量Q，横坐标表示位置

Figure 5. Numerical results of the example in Section 4.2, the vertical axis represents the width of the surface level and miss flow rate Q respectively, while the horizontal axis represents the position

图5. 4.2节算例的数值结果，纵坐标表示表面水平和质量流量Q，横坐标表示位置

致 谢

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