求解多孔浅水波方程组的高精度和谐数值方法

High-Accuracy Numerical Methods for Solving the System of Porous Shallow Water Wave Equations

摘要: 本文提出了一种高阶平衡(Runge-Kutta Discontinuous Galerkin, RKDG)有限元方法，用于解决具有不连续孔隙度和复杂底部地形的多孔浅水波方程组。该方法能够在静水稳态解的条件下保持良好的平衡特性，即当数值通量梯度与源项在连续意义上完全平衡时，离散格式也能准确保持稳态解。与传统的恒定横截面浅水波方程相比，由于孔隙度和底部地形的非连续性，构造具有良好平衡性质的高阶数值格式更加复杂。为应对这一挑战，我们首先对源项进行了重构，并结合基于Lax-Friedrichs通量的静水重构方法，同时合理设计离散格式，确保该方法在处理不连续区域时仍能保持高阶精度和稳定性。通过一系列标准数值算例验证了所提方法的有效性，结果表明，该方法在保持静水平衡、实现高阶精度以及分辨不连续解方面具有优异的表现。

Abstract: In this paper, we propose a high-order balanced Runge-Kutta Discontinuous Galerkin (RKDG) finite element method for solving the porous shallow water wave equations with discontinuous porosity and complex bottom topography. The proposed method preserves the good balance properties of the steady-state solution, meaning that when the numerical flux gradient and the source terms are fully balanced in a continuous sense, the discretization scheme can also accurately maintain the steady-state solution. Compared to traditional shallow water wave equations with constant cross-sections, the discontinuities in porosity and bottom topography make it more challenging to construct a high-order numerical scheme with good balance properties. To overcome this challenge, we first reconstruct the source terms and apply a hydrostatic reconstruction method based on the Lax-Friedrichs flux, while also designing the discretization scheme to ensure high-order accuracy and stability when handling discontinuous regions. A series of standard numerical test cases validate the effectiveness of the proposed method. The results demonstrate that the method performs exceptionally well in maintaining hydrostatic balance, achieving high-order accuracy, and resolving discontinuous solutions.

1. 引言

多孔浅水波方程已经通过数值方法进行离散化处理，并广泛用于城市洪水模拟[1] [2]和植被区域水流研究[3] [4]等领域。忽略能量损失的多孔浅水波方程的质量和动量方程可表示为：

$\{\begin{array}{l}{\left(\varphi h\right)}_t+{\left(\varphi hu\right)}_x=0\hfill \\ {\left(\varphi hu\right)}_t+{\left(\varphi h{u}^2+\frac{1}{2}g\varphi h^2\right)}_x=\frac{1}{2}g{h}^2{\varphi }_x-g\varphi h{b}_x\hfill \end{array}$ (1.1)

其中$\varphi$ 表示孔隙度，$b$ 为底部地形，$h$ 为水面高度，$u$ 为速度，$g$ 为重力常数，$Q=\varphi hu$ 为质量流量。

模型(1.1)是带有源项的双曲守恒律(即双曲平衡律)，并允许存在如下的静水定常解：

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

在流体力学中，源项与通量之间的平衡关系在数值模拟中尤为重要。此模型主要的困难是在非守恒源项中孔隙度和底部地形的非连续性可以同时存在于一个位置，模型(1.1)的x导数项无法按经典方式定义。然而，为了反映不连续的孔隙度和底部地形的影响，必须正确地确定不连续几何形状对流动产生的力。在处理非连续孔隙度与非平坦底部地形的情况下，传统数值方法在离散化过程中可能破坏稳态解下源项与通量梯度之间的精确平衡。这种不平衡，即使在细网格下，也可能导致不稳定的附加振动，进而引起较大误差。其根本原因在于截断误差破坏了源项与通量之间的微妙平衡。为此，我们提出了一种具有良好平衡特性的格式，能够更好地保持静态水平衡，并有效捕捉到可能出现的小幅度误差。

本研究中，我们重点探讨高阶RKDG方法[5] [6]，该方案具有以下优点。首先，它的实现相对简单，便于应用。其次，计算效率较高，能够有效减少计算时间。此外，该方法允许近似解在单元交界处出现间断，这为处理复杂问题提供了灵活性。最后，与有限体积方法相比，RKDG方法避免了由于扩展节点模板以提高精度而带来的计算复杂性和非局部依赖性，因而在处理大规模问题时更具优势。

2. 标准RKDG方法

本节将回顾并简要介绍标准RKDG方法。首先我们介绍RKDG方法的一些标准符号。区间$I=\left[a,b\right]$ 被分成$N$ 个子区间并且我们用$I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right],j=1,\cdots ,N$ 来表示单元。首先假设网格均匀分布，单元的大小为$\Delta x=x_{j-\frac{1}{2}}-x_{j+\frac{1}{2}}$ ，并且我们使用$x_j=\frac{1}{2}\left(x_{j+\frac{1}{2}}+x_{j-\frac{1}{2}}\right)$ 来表示单元中心。分段多项式空间$V_h^K$ 是指由每个单元$I_j$ 上次数不超过指定值$K$ 的多项式所构成的集合：

$V_h^K=\left\{\nu \left(x\right):{\nu |}_{I_j}\in P^K\left(I_j\right),j=1,2,\cdots ,N\right\}$

空间$V_h^K$ 允许其中的元素在单元界面$x_{j+1/2}$ 处存在间断，所以我们在这里所了解的方法也被称为间断方法。因为$V_h^K$ 为有限维空间，故其包含有限个基函数。

DG方法

我们首先考虑了一维标量的双曲守衡方程[7]，标准形式为：

$\{\begin{array}{l}{u}_t+f{\left(u\right)}_x=0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),\hfill \end{array}x\in \left[a,b\right].$ (2.1)

接下来，使用Galerkin方法，我们使试探函数$\nu \left(x\right)\in V_h^K$ ，并在单元$I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right]$ 上进行积分：

我们在每个单元$I_j$ 上寻找未知量$u\left(x,t\right)$ 的DG近似：

${u_h\left(x,t\right)|}_{x\in I_j}={\displaystyle \sum _{l=0}^K{u}_j^l}\left(t\right)v_j^{\left(l\right)}\left(x\right)\in V_h^K,$

并且对于所有试探函数$\nu \left(x\right)\in V_h^K$ ，满足以下等式：

其中，

${\widehat{f}}_{j+\frac{1}{2}}=\widehat{f}\left(u_h\left(x_{j+\frac{1}{2}}^{-},t\right),u_h\left(x_{j+\frac{1}{2}}^{+},t\right)\right),$

这里$\widehat{f}\left(a_1,a_2\right)$ 是数值通量，用来逼近$f\left(u\left(x_{j+\frac{1}{2}},t\right)\right)$ 。我们通常用到的Lax-Friedrichs流量：

$\widehat{f}\left(a_1,a_2\right)=\frac{1}{2}\left(f\left(a_1\right)+f\left(a_2\right)-\alpha \left(a_2-a_1\right)\right),$

这里$\alpha =\max \left|f^{\prime }\left(u\right)\right|$ ，对于双曲守衡方程，我们取$\alpha =\max \left|\lambda \left(u\right)\right|$ ，其中$\lambda \left(u\right)$ 为雅可比矩阵$\frac{\partial f\left(u\right)}{\partial u}$ 的特征值。

3. 良好平衡的RKDG方法

3.1. 改进RKDG方法

本文目标提出改进后的RKDG方法，以提高在模拟物理现象时的计算精度和数值稳定性。主要目的在于证明，只需对传统的RKDG方法的初始值或通量进行适度的修改，即可实现良好的平衡性质，为构建高阶平衡格式提供了一种极为简洁的途径。与传统的RKDG方法相比，改进后的格式仅需很少的额外计算成本便能实现很好的平衡特性，并获得高精度的数值模拟结果。

对于多孔浅水波方程(1.1)，我们关注的是保持静水稳态解(1.2)。对于此类静水解，第一个方程${\left(hu\right)}_x=0$ 对任何格式都能精确满足。我们将集中讨论第二个方程，它可以表示为：

${\left(hu\right)}_t+f{\left(U\right)}_x=s\left(\varphi ,h,b\right),$

其中$U={\left(h,hu\right)}^{\text{T}}$ ，T代表转置。

在RKDG方法[8] [9] [10]中，$U$ 通过在分段多项式空间$V_h^K$ 近似为$U_h$ 。我们将底部函数$b$ 投影到同样的空间$V_h^K$ 中获得近似值$b_h$ 。这意味着如果$h+b=const$ ，则$h_h+b_h=const$ 。我们的数值格式具有如下形式：

(3.1)

与传统RKDG比较，我们可以看到，通量${\widehat{f}}_{j+\frac{1}{2}}$ 和${\widehat{f}}_{j-\frac{1}{2}}$ 被左侧通量${\widehat{f}}_{j+\frac{1}{2}}^l$ 和右侧通量${\widehat{f}}_{j-\frac{1}{2}}^r$ 替代。我们把上述格式改写为：

(3.2)

其中${\widehat{f}}_{j+\frac{1}{2}}=\widehat{F}\left(U_h\left(x_{j+\frac{1}{2}}^{-},t\right),U_h\left(x_{j+\frac{1}{2}}^{+},t\right)\right)$ 。(3.2)的左侧是传统的RKDG格式，而右侧是本文对源项的逼近方法。左通量${\widehat{f}}_{j+\frac{1}{2}}^l$ 与右通量${\widehat{f}}_{j-\frac{1}{2}}^r$ 的具体设计将在后文阐述，此处${\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j+\frac{1}{2}}^l$ 和${\widehat{f}}_{j-\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}^r$ 是$O\left(\Delta x^{k+1}\right)$ 的高阶修正项。因此，格式(3.1)是一个k + 1阶守恒格式，并且将收敛于弱解。

孔隙度的空间变化可能导致复杂的数值计算，尤其是在不连续或快速变化的情况下。为了准确捕捉孔隙度对流体流动的影响，必须对孔隙度进行合理的数值处理。我们将孔隙度定义为：${\varphi }_j\left(x\right)=\frac{1}{\text{Δ}x}{\displaystyle {\int }_{I_j}\varphi \left(x\right)\text{d}x}$ 。在边界处定义唯一孔隙度${\varphi }_{h,j+\frac{1}{2}}={\varphi }_{j+\frac{1}{2}}=\varphi \left(x_{{}_{j+\frac{1}{2}}}\right)$ ，避免因孔隙度不连续或变化剧烈所引发的数值不稳定。能够有效消除不连续带来的影响，从而提高源项积分的精度。

为了获得良好的平衡性质，我们需要计算残差：

(3.3)

为了使静水状态到零误差残差，需满足以下三个条件：

• 公式(2.2)中的所有积分都应使用适当的高斯求积去精确计算。显然，因为$h_h$ 、$b_h$ 和$v_h$ 在每个单元$I_j$ 中都是多项式，因此$f$ 和$s$ 都是多项式。

• 我们假设：

${\widehat{f}}_{j+\frac{1}{2}}^l=f\left(U_h\left(x_{j+\frac{1}{2}}^{-},t\right)\right),{\widehat{f}}_{j-\frac{1}{2}}^r=f\left(U_h\left(x_{j-\frac{1}{2}}^{+},t\right)\right).$ (3.4)

对于静水问题这个条件并不是显而易见的，稍后我们将讨论如何使其在RKDG方法中成立。

• 我们假设$U$ 的数值近似解$U_h$ ，是方程${\left(hu\right)}_t+f{\left(U\right)}_x=s\left(\varphi ,h,b\right)$ 的定常解。这意味着：

${\left(\frac{1}{2}g{\varphi }_h{h}_h^2\right)}_x=\frac{1}{2}g{h}_h^2{\left({\varphi }_h\right)}_x-g{\varphi }_h{h}_h{\left(b_h\right)}_x\text{or}{\partial }_{x}f\left(U_h\right)=s\left({\varphi }_h,h_h,b_h\right).$ (3.5)

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

我们定义一个新的向量：

$h_{h,j+\frac{1}{2}}^{\ast ,+}=\max \left(0,h_{h,j+\frac{1}{2}}^{+}+b_{h,j+\frac{1}{2}}^{+}-\max \left(b_{h,j+\frac{1}{2}}^{-},b_{h,j+\frac{1}{2}}^{+}\right)\right),$ (3.6)

$h_{h,j+\frac{1}{2}}^{\ast ,-}=\max \left(0,h_{h,j+\frac{1}{2}}^{-}+b_{h,j+\frac{1}{2}}^{-}-\max \left(b_{h,j+\frac{1}{2}}^{-},b_{h,j+\frac{1}{2}}^{+}\right)\right),$ (3.7)

并重新定义$U$ 为：

$U_{h,j+\frac{1}{2}}^{\ast ,\pm }=\left(\begin{array}{c}{\varphi }_{h,j+\frac{1}{2}}^{}{h}_{h,j+\frac{1}{2}}^{\ast ,\pm }\\ {\left(\varphi hu\right)}_{h,j+\frac{1}{2}}^{\pm }\end{array}\right).$ (3.8)

左通量${\widehat{f}}_{j+\frac{1}{2}}^l$和右通量${\widehat{f}}_{j-\frac{1}{2}}^r$ 修正为：

${\widehat{f}}_{j+\frac{1}{2}}^l=\widehat{F}\left(U_{h,j+\frac{1}{2}}^{\ast ,-},U_{h,j+\frac{1}{2}}^{\ast ,+}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{-}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{\ast ,-}\right)}^2\end{array}\right),$ (3.9)

${\widehat{f}}_{j-\frac{1}{2}}^r=\widehat{F}\left(U_{h,j-\frac{1}{2}}^{\ast ,-},U_{h,j-\frac{1}{2}}^{\ast ,+}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{+}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{\ast ,+}\right)}^2\end{array}\right),$ (3.10)

在这里，我们常用的数值通量是Lax-Friedrichs通量：

$\widehat{F}\left(a,b\right)=\frac{1}{2}\left(f\left(a\right)+f\left(b\right)-\alpha \left(b-a\right)\right),$ (3.11)

其中，$\alpha ={\max }_u\left|F^{\prime }\left(u\right)\right|$ 。很容易得到${\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j+\frac{1}{2}}^l$ 和${\widehat{f}}_{j-\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}^r$ 在一般解的情况下达到$O\left(\Delta x^{k+1}\right)$ 的精度，因此最初的k + 1阶精度得以保持。在静水条件下$h_h+b_h=const$ ，所以

${\widehat{f}}_{j+\frac{1}{2}}^l=\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{\ast ,-}\right)}^2\end{array}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{-}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{\ast ,-}\right)}^2\end{array}\right)=\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{-}\right)}^2\end{array}\right)=f\left(U_{h,j+\frac{1}{2}}^{-}\right),$ (3.12)

同理可得：

${\widehat{f}}_{j-\frac{1}{2}}^r=f\left(U_{h,j-\frac{1}{2}}^{+}\right).$ (3.13)

3.3. 证明良好平衡性质

通过使用修正的数值通量，结合源项的精确近似，基于RKDG方法我们能够重新构造数值格式(3.1)。核心目标是确保在通量和源项间保持精确平衡，从而在静水条件下实现稳定解。数值格式重新写为：

命题1. 满足上述三个条件的RKDG方法，对于多孔浅水波方程的静水解问题是精确的。

证明.

数值通量${\widehat{f}}_{j+\frac{1}{2}}^l$ 为：

$\begin{array}{c}{\widehat{f}}_{j+\frac{1}{2}}^l=\widehat{F}\left(U_{h,j+\frac{1}{2}}^{\ast ,-},U_{h,j+\frac{1}{2}}^{\ast ,+}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{-}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{\ast ,-}\right)}^2\end{array}\right)\\ =\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{\ast ,-}\right)}^2\end{array}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{-}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{\ast ,-}\right)}^2\end{array}\right)=\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j+\frac{1}{2}}^{-}\right)}^2\end{array}\right)\\ =f\left(U_h\left(x_{j+\frac{1}{2}}^{-},t\right)\right),\end{array}$

数值通量${\widehat{f}}_{j-\frac{1}{2}}^r$ 为：

$\begin{array}{c}{\widehat{f}}_{j-\frac{1}{2}}^r=\widehat{F}\left(U_{h,j-\frac{1}{2}}^{\ast ,-},U_{h,j-\frac{1}{2}}^{\ast ,+}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{+}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{\ast ,+}\right)}^2\end{array}\right)\\ =\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{\ast ,+}\right)}^2\end{array}\right)+\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{+}\right)}^2-\frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{\ast ,+}\right)}^2\end{array}\right)=\left(\begin{array}{c}0\\ \frac{g}{2}\varphi {\left(h_{h,j-\frac{1}{2}}^{+}\right)}^2\end{array}\right)\\ =f\left(U_h\left(x_{j-\frac{1}{2}}^{+},t\right)\right).\end{array}$

则$R$ 对于静水问题将简化为：

因此，我们验证了该格式在数值解法中具有良好的平衡性。

3.4. 平衡格式的总结

我们对求解多孔浅水波方程(1.1)进行了总结。采用良好平衡RKDG方法的完整过程具体步骤如下：

1. 通过对源项的重新构造，将方程转化为形式(3.1)，以便与RKDG方法兼容并保证其良好的平衡性。

2. 利用RKDG方法构造数值通量${\widehat{f}}_{j+\frac{1}{2}}^{\pm }$ ，通过对传统RKDG方法的修改确保稳态解能够被准确捕捉。

3. 我们将数值通量的余项与源项进行近似求和，并基于这一近似采用龙格–库塔(Runge-Kutta)方法推进时间，并通过迭代过程逐步提高解的精度，最终实现平衡。

4. 数值结果

本节展示了通过一系列基础数值实验验证所提多孔浅水波方程高阶RKDG方法的效果。实验采用了RKDG方法与三阶TVD龙格–库塔方法[11]的耦合，CFL值设为0.6，重力常数g取9.812。实验涵盖了连续和不连续孔隙度的情形，并且所有测试结果均已得到验证。

4.1. 良好平衡测试

第二个测试问题[12] [13]旨在验证所提出的RKDG方法的平衡性。我们考虑在区域[0, 1]内给定的底部地形：

$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)

孔隙度$\varphi \left(x\right)$ 的空间变化形式为：

$\varphi \left(x\right)=\{\begin{array}{l}1-{\varphi }_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{\varphi }_0$ 表示在$\left(x_l+x_r\right)/2$ 点处的最小孔隙度。在本例中，我们选取$x_l=0.25,x_r=0.75$ 和${\varphi }_0=0.2$ 作为初始条件，施加周期边界条件，并给出如下稳态解：

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

采用包含200个均匀单元的网格计算到时刻$t=1$ 。图1显示了表面高度$h+b$ 和底部地形$b$ 的变化。

Figure 1. Numerical results of the example in Section 4.1

图1. 4.1节算例的数值结果

4.2. 小扰动测试

本节中我们通过测试验证了该模型在稳定性方面的表现。该测试的过程和结果展示了使用RKDG方法时，模型能够有效解决在具有挑战性边界条件下的问题。该测试最早由[12]提出，初始条件为：

$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=\varphi hu=0,$

在计算区域[0, 1]中，采用简单的传输边界条件。测试了两组孔隙度$\varphi \left(x\right)$ ，第一组为孔隙度左移收缩$x_l=0.15,x_r=0.65,{\varphi }_0=0.2$ ，第二组为孔隙度右移收缩$x_l=0.\text{3}5,x_r=0.\text{8}5,{\varphi }_0=0.2$ 。采用了小范围的非平衡数值方法来处理可能出现的数值不稳定问题[14]。通过对200个均匀划分的数据点进行数值计算，使用优化后的RKDG方法与2000个单元的参考数据进行对比(见图2)。

Figure 2. Numerical results of the example in Section 4.2

图2. 4.2节算例的数值结果

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

本数值算例中[13]，我们研究了RKDG方法在不同孔隙度条件下的稳定性。在[15]中，使用相同的试验验证了浅水波方程的良好平衡RKDG方法的性能。计算区域设为[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根据不同的边界条件进行调整。为便于对比，还提供了数值解的解析解[16]。对流速和水位进行分段处理，并应用水量守恒以保证数值解的精度。

(1) 亚临界流

在此测试中，上游边界条件设为$hu=\text{4}\text{.42}$ ，下游边界条件设为$h=2$ 。考虑(4.2)中定义的两组孔隙度分布$\varphi \left(x\right)$ ，第一组孔隙度为左移收缩$x_l=3.75,x_r=16.25,{\varphi }_0=0.05$ ，第二组孔隙度为右移收缩$x_l=8.75,x_r=21.25,{\varphi }_0=0.05$ 。随着时间推移，当$t=200$ 时，水流最终稳定，形成亚临界流。此时，流速和水位不再发生显著变化，达到平衡状态。图3中展示了数值解与解析解的对比。结果表明，数值解与解析解基本一致，且孔隙度的不同对稳态解及其收敛性产生了影响。

Figure 3. Numerical results of the example in Section 4.3

图3. 4.3节算例的数值结果

(2) 无激波跨临界流动

在此测试中，上游边界条件设为$hu=\text{1}\text{.53}$ ，下游边界条件设为$h=0.66$ 。考虑(4.2)中定义的两组不同的孔隙度分布$\varphi \left(x\right)$ ，第一组孔隙度为左移收缩$x_l=3.75,x_r=16.25,{\varphi }_0=0.15$ ，第二组孔隙度为右移收缩$x_l=8.75,x_r=21.25,{\varphi }_0=0.15$ 。该稳态解为无激波的跨临界流。图4展示了表面水平$h+b$ 和质量流量$Q$ 的数值解与解析解的比较。结果表明，数值解与解析解的高度一致。

Figure 4. Numerical results of the example in Section 4.3

图4. 4.3节算例的数值结果

5. 结论

本文针对具有不规则孔隙度和非平坦底部地形的多孔浅水波方程，提出了高阶RKDG方法。该方法采用静水重构近似数值通量，并合理离散源项，从而保持通量梯度与源项之间的精确平衡，实现良好的平衡特性与静水稳态保持。理论分析和数值实验表明，所提出的高阶RKDG方法在保持静水平衡性、分辨不连续结构以及在光滑区域实现高阶精度方面均表现出优异性能。

致 谢

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

NOTES

^*通讯作者。