DOI: 10.12677/AAM.2023.128367, PDF, HTML, XML,  被引量    科研立项经费支持
作者: 周翔宇, 张志壮, 钱守国, 李 刚^*：青岛大学数学与统计学院，山东 青岛
关键词: 浅水波方程；ADER方法；间断Galerkin格式；保证格式；微分变换过程；全离散；Shallow Water Wave Equations； ADER Approach； Discontinuous Galerkin Scheme； Positivity-Preserving Method； Differential Transformation Procedure； Fully-Discrete

文章引用：周翔宇, 张志壮, 钱守国, 李刚. 浅水波方程的高阶保正Well-Balanced ADER间断Galerkin格式[J]. 应用数学进展, 2023, 12(8): 3728-3743. https://doi.org/10.12677/AAM.2023.128367

1. 引言

在本文中，我们主要的工作是构造非平坦底部地形上的浅水波方程(SWEs)的保正well-balanced性质的间断Galerkin格式，该方程具体形式如下 [1] [2] [3] [4] ：

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

在上式中，h为水深，u为水深平均速度，b为水底地形。g是重力加速度常数。

浅水波方程(1)属于双曲守恒律的一种，并且其静水稳态解由下式给出：

$u=0,h+b=\text{const},$ (2)

其中源项与通量梯度完全平衡。这种平衡律的数值模拟，包括非平坦底部地形上的非线性浅水波方程，有一个主要的问题：标准数值方法可能无法在离散水平上保持稳态下通量梯度和源项之间的平衡，并且可能在稳态附近产生伪振荡。Well-balanced方法旨在克服这一问题，并且在稳态（或接近稳态）解决格式下表现良好。对于非线性SWEs，已经有多种方法被提出并应用，其中包括 [5] - [12] 及其参考文献。在相关问题的数值模拟中，另一个经常遇到的问题是对少水或无水地区的润湿和干燥情况的处理。典型的问题包括溃坝问题、洪波和爬高现象。数值上，在计算过程中可能会产生负值，这可能会带来额外的困难。有许多现有的保正方法可以克服这一困难，关于这些方法我们参考了文献 [13] - [18] 。

上述对浅水波方程的所有工作都是在有限差分法或有限体积法的基础上进行的。在过去的几十年中，高阶有限元Galerkin(DG)方法在求解包含双曲型守恒律的偏微分方程方面得到了广泛的关注。

近十多年来，高阶时间和空间任意导数方法在求解对流占优问题中得到了广泛的关注。Toro和合作者最初在双曲守恒律的有限体积格式背景下提出了ADER格式 [19] [20] [21] 。格式 [19] [20] [21] 是单步、单级和全离散的。后来，Dumbser和Munz [22] [23] [24] 提出了一种采用单阶段ADER方法 [19] [20] [21] 的DG方法，用于时间离散化。随后，Dumbser等人对该方法进行了一系列研究 [25] [26] [27] [28] 。

DG方法 [25] [26] [27] [28] 的关键是构造高阶数值通量。网格间状态以时间泰勒级数展开的形式表示，其中时间导数由使用控制PDE的空间导数表示。该程序被称为Cauchy-Kowalewski程序(在文献中也称为Lax-Wendroff程序 [29] )。此外，还需要求解网格间的广义黎曼问题，以获得状态本身及其空间导数。然而，Cauchy-Kowalewski方法非常繁琐，尤其是在求解高维问题的时候。因此，CauchyKowalewski程序的替代或简化是非常有意义的。值得注意的是，Dumbser和Munz [30] 提出了一种高效的Cauchy-Kowalewski过程算法。最近，Dumbser等人 [31] ，Duan和Tang [32] 采用局部连续时空Galerkin预测方法 [33] 来取代DG方法框架下的Cauchy-Kovalewski过程。

Galerkin (DG)方法是一类以高阶间断分段多项式空间作为测试函数空间的有限元方法 [34] [35] [36] 。半离散DG方法(也称为RKDG方法)使用间断的分段多项式空间作为数值解和测试函数空间，并通过TVD Runge-Kutta格式进行稳定的高阶时间积分 [37] - [42] 。Qian [43] 设计了几何形状不规则的明渠浅水波方程的高阶保正平衡方法。

最近，ADER方法由Toro和Titarev引入到有限体积格式中，该方法通过完全离散的一步公式获得空间和时间上任意高阶精度。原始版本的ADER方法在空间和时间上使用截断的Taylor展开，并通过Cauchy-Kovalevskaya (C-K)过程 [44] 将时间导数替换为空间导数，但是该过程不能处理刚性源项。为了克服这些问题，ADER-DG方法用局部时空DG预测器取代C-K过程 [45] 。因此，单步公式只需要逐点计算通量和源项，这样就极大地减少了空间元素之间的模板依赖，并允许几乎完美的并行化。

在本文中，我们构造了湖处于静止稳态的SWEs的高阶良好平衡的ADER-DG方法，并用微分变换代替C-K过程。为了达到良好的平衡性能，我们建议结合基于流体静力重建思想的数值通量，以及一种新的源项近似和分解算法。所得到的DG方法在保持湖泊静止稳态的well-balanced特性的同时，又保持了光滑解的高阶精度。

本文结构如下：在第2节中，我们提出了保正平衡的ADER-DG方法。第3节给出了一些数值算例来验证所提出的ADER-DG方法的性能。最后，第4节给出结论。

2. ADER-DG方法

在本节中，我们构造了浅水波方程的高阶良好平衡的ADER-DG方法(1)。为了准确地保持静水稳态(2)，我们基于流体静力重建技术构建了well-balanced的数值通量，并提出了一种新的基于分解算法的源项近似。

2.1. 符号记法

首先将空间区间 $\left[a,b\right]$ 离散成N个单元 $I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right],j=1,2,\cdots ,N$ 。每个单元格的中心是 $x_j=\frac{1}{2}\left(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}}\right)$ ，单元格 $I_j$ 的大小记为 ${\tau }_j=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$ ，最大网格大小记为 $\tau =\underset{1\le j\le N}{\max }{\tau }_j$ 。符号 ${\Omega }_j=I_j\times \left[t^n,t^{n+1}\right]$ 表示时空控制体积， $\left[t^n,t^{n+1}\right]$ 为时间间隔。在这里，我们取近似解空间如下：

$V_{\tau }^k=\left\{{\phi \left(x,t\right):\phi \left(x,t\right)|}_{{\Omega }_j}\in P^k\left({\Omega }_j\right)\right\},$ (3)

其中 $P^k\left({\Omega }_j\right)$ 表示 ${\Omega }_j$ 上k次以下的时空多项式集合。注意， $V_{\tau }^k$ 中的函数允许在元件界面上具有间断。此外，我们使用 ${\left\{v_{\tau }\right\}}_{j+\frac{1}{2}}=\left(v_{\tau ,j+\frac{1}{2}}^{+}+v_{\tau ,j+\frac{1}{2}}^{-}\right)/2$ 表示函数 $v_{\tau }$ 在元素界面 $x_{j+\frac{1}{2}}$ 处的算术平均值。

2.2. 微分变换过程

在本节中，我们将介绍一种基于将微分方程转化为其时空序列展开解的系数递推关系的方法，称为微分变换法(DTM) [46] 。值得注意的是，该方法在评估C-K过程时更加简洁高效 [47] 。为方便起见，我们首先考虑一维平流方程

$\begin{array}{l}{q}_t+q_x=0,\\ q\left(x,0\right)=q_0\left(x\right).\end{array}$ (4)

将 $q\left(x,0\right)$ 的初始条件展开为泰勒多项式

$q_{\tau }\left(x,0\right)={\displaystyle \underset{k_x=0}{\overset{\infty }{\sum }}\tilde{q}\left(k_x,0\right){\left(x-x_j\right)}^{k_x}}.$ (5)

接下来，我们假设解可以表示为每个时空域内的以下多项式

$q_{\tau }\left(x,t\right)={\displaystyle \underset{k_x=0}{\overset{\infty }{\sum }}{\displaystyle \underset{k_t=0}{\overset{\infty }{\sum }}\tilde{q}\left(k_x,k_t\right){\left(x-x_j\right)}^{k_x}{\left(t-t^n\right)}^{k_t},\left(x,t\right)\in {\Omega }_j}.}$ (6)

微分变换 $\tilde{q}$ 的形式是这样的

$\tilde{q}\left(k_x,k_t\right)=\frac{1}{k_x!k_t!}{\frac{{\partial }^{k_x+k_t}q\left(x,t\right)}{\partial x^{k_x}\partial t^{k_t}}|}_{x=x_j,t=t^n},k_x,k_t=0,1,2,\cdots .$ (7)

将(6)代入平流式(4)，可得递归关系如下：

$\tilde{q}\left(k_x,k_t+1\right)=-\frac{k_x+1}{k_t+1}\tilde{q}\left(k_x+1,k_t\right).$ (8)

因此，从式(5)中的初始微分变换 $\tilde{q}\left(k_x,0\right)$ 开始，反复重复式(8)，得到任意时刻的微分变换。因此，我们可以通过泰勒级数展开(6)得到平流方程(4)的解，即任意高阶解。

随后，我们考虑了浅水波方程的微分变换方法(13)。我们将 $U\left(x,0\right)$ 的初始条件投影到分段多项式空间 $V_{\tau }^k$ 中，得到 $U_{\tau }\left(x,0\right)$ ，其形式为

$U_{\tau }\left(x,0\right)={\displaystyle \underset{k_x=0}{\overset{k}{\sum }}\tilde{U}\left(k_x,0\right){\left(x-x_j\right)}^{k_x}}.$ (9)

与Runge-Kutta DG方法不同，我们通过k次的时空多项式在每个时空域 ${\Omega }_j$ 内寻找浅水波方程(13)的数值解，如下所示

$U_{\tau }\left(x,t\right)={\displaystyle \underset{k_x=0}{\overset{k}{\sum }}{\displaystyle \underset{k_t=0}{\overset{k}{\sum }}\tilde{U}\left(k_x,k_t\right){\left(x-x_j\right)}^{k_x}{\left(t-t^n\right)}^{k_t}.}}$ (10)

与2.2节类似，我们将(10)代入(45)，得到浅水波方程(12)的递推关系如下：

$\begin{array}{l}\tilde{h}\left(k_x,k_t+1\right)=-\frac{k_x+1}{k_t+1}\left(\tilde{h}\tilde{u}\right)\left(k_x+1,k_t\right),\\ \left(\tilde{h}\tilde{u}\right)\left(k_x,k_t+1\right)=-\frac{k_x+1}{k_t+1}{\tilde{G}}_1\left(k_x+1,k_t\right)-\frac{k_x+1}{k_t+1}{\tilde{G}}_2\left(k_x+1,k_t\right)-\frac{k_x+1}{k_t+1}{\tilde{G}}_3\left(k_x+1,k_t\right).\end{array}$ (11)

其中

$\begin{array}{l}{\tilde{G}}_1\left(k_x,k_t\right)={\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}{\tilde{G}}_{1,a}\left(r,s\right){\tilde{G}}_{1,b}}}\left(k_x-r,k_t-s\right),\\ {\tilde{G}}_2\left(k_x,k_t\right)=\frac{g}{2}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}\tilde{h}\left(r,s\right)\tilde{h}\left(k_x-r,k_t-s\right)}},\\ {\tilde{G}}_3\left(k_x,k_t\right)={\displaystyle \underset{r=0}{\overset{k_x}{\sum }}\tilde{h}\left(r,k_t\right){\tilde{G}}_{3,a}\left(k_x-r\right)},\\ {\tilde{G}}_{1,a}\left(k_x,k_t\right)={\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}\left(\tilde{h}\tilde{u}\right)\left(r,s\right)}}\left(\tilde{h}\tilde{u}\right)\left(k_x-r,k_t-s\right),\\ {\tilde{G}}_{2,b}\left(k_x,k_t\right)=-\frac{1}{\tilde{h}\left(0,0\right)}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}\tilde{h}\left(r,s\right){\tilde{G}}_{1,b}\left(k_x-r,k_t-s\right)}},\\ {\tilde{G}}_{3,a}\left(k_x\right)=g\left(k_x+1\right)\tilde{b}\left(k_x+1\right).\end{array}$

因此，将初始微分变换 $\tilde{U}\left(k_x,0\right)$ 和 $\tilde{b}\left(k_x\right),k_x=0,\cdots ,k$ 代入(11)，得到微分变换 $\tilde{U}\left(k_x,k_t\right)$ ，并利用(10)得到在时空域 ${\Omega }_j$ 上的近似 $U\left(x,t\right)$。

2.3. ADER-DG方法的构建

在本节中，我们将利用2.2节中的微分变换过程构建SWEs(1)的ADER-DG方法。

首先，我们将原控制方程(1)改写为

$\begin{array}{l}{h}_t+{\left(hu\right)}_x=0,\\ {\left(hu\right)}_t+{\left(\frac{{\left(hu\right)}^2}{h}+\frac{1}{2}g{h}^2\right)}_x=\frac{1}{2}g{\left(b^2\right)}_x-g\left(h+b\right)x_x.\end{array}$ (12)

我们将源项 $-gh{b}_x$ 按照 [48] [49] 中的方法改写为动量方程中的等效形式 $\frac{1}{2}g{\left(b^2\right)}_x-g\left(h+b\right)b_x$ 。

为了便于表示，我们引入了可以将(12)改写为紧致向量形式的符号

$U_t+F\left(U\right)=S\left(U\right),$ (13)

其中 $U={\left(h,hu\right)}^{\text{T}},F\left(U\right)={\left(hu,h{u}^2+\frac{1}{2}g{h}^2\right)}^{\text{T}}$ 和 $S\left(U\right)={\left(0,-g\left(h+b\right)b_x+\frac{1}{2}g{b}_x^2\right)}^{\text{T}}$ 分别表示保守变量、物理通量和几何源项。

式(13)与任意空间测试函数 $\phi \left(x\right)$ 的乘法并使用 ${\Omega }_j$ 上的部分积分产生以下弱形式

$\begin{array}{l}{\displaystyle \underset{I_j}{\int }U\left(x,t^{n+1}\right)\varphi \left(x\right)\text{d}x}-{\displaystyle \underset{I_j}{\int }U\left(x,t^n\right)\varphi \left(x\right)\text{d}x}+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}F\left(x_{j+\frac{1}{2}},t\right)\text{d}t\cdot \varphi \left(x_{j+\frac{1}{2}}^{-}\right)}\\ -{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}F\left(x_{j-\frac{1}{2}},t\right)\text{d}t\cdot \varphi \left(x_{j-\frac{1}{2}}^{+}\right)}-{\displaystyle \underset{{\Omega }_j}{\iint }F\left(x,t\right)\varphi \left(x\right)\text{d}x\text{d}t}\\ ={\displaystyle \underset{{\Omega }_j}{\iint }S\left(U_{\tau }\left(t\right)\right)\varphi \left(x\right)\text{d}x\text{d}t}\text{for}j=1,2,\cdots ,N\text{.}\end{array}$

这里，为了便于表示，我们将物理通量写作一个二元函数 $F\left(x,t\right)=F\left(U\left(x,t\right)\right)$ 。

求解(13)的全离散ADER-DG方法定义如下：对于所有测试函数 $\varphi \left(x\right)={\left(x-x_j\right)}^{k_x}$ 且 $k_x=0,1,\cdots ,k$ ，从 $V_{\tau }^k$ 出发，解 $U_{\tau }\left(x,t\right)\in V_{\tau }^k$ 满足

$\begin{array}{l}{\displaystyle \underset{I_j}{\int }{U}_{\tau }\left(x,t^{n+1}\right)\varphi \left(x\right)\text{d}x}-{\displaystyle \underset{I_j}{\int }{U}_{\tau }\left(x,t^n\right)\varphi \left(x\right)\text{d}x}-{\hat{F}}_{j+\frac{1}{2}}\cdot \varphi \left(x_{j+\frac{1}{2}}^{-}\right)+{\hat{F}}_{j-\frac{1}{2}}\cdot \varphi \left(x_{j-\frac{1}{2}}^{+}\right)\\ +{\displaystyle \underset{{\Omega }_j}{\iint }{F}_{\tau }\left(x,t\right)\varphi \left(x\right)\text{d}x\text{d}t}+{\displaystyle \underset{{\Omega }_j}{\iint }S\left(U_{\tau }\left(t\right)\right)\varphi \left(x\right)\text{d}x\text{d}t}j=1,2,\cdots ,N.\end{array}$ (14)

这里用数值通量 ${\hat{F}}_{j+\frac{1}{2}}$ 来近似网格之间物理通量的时间积分，即 ${\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}F\left(x_{j+\frac{1}{2}},t\right)\text{d}t}$ 。其中 $\hat{F}\left(U_{\tau ,j+\frac{1}{2}}\left(t\right)\right)=f\left(U_{\tau ,j+\frac{1}{2}}^{-}\left(t\right),U_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\right)$ 数值通量与简单的Lax-Friedrichs通量一致

$f\left(U_{\tau ,j+\frac{1}{2}}^{-}\left(t\right),U_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\right)=\frac{1}{2}\left[F\left(U_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)\right)+F\left(U_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\right)-\alpha \left(U_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)-U_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)\right)\right].$ (15)

其中

$\alpha =\underset{x\in I}{\max }\left|\frac{1}{\Delta t^n}{\displaystyle {\int }_{t^n}^{t^{n+1}}\lambda \left(U\left(x,t\right)\right)\text{d}t}\right|$ . (16)

其中 $\lambda \left(U\right)$ 为雅可比矩阵 $F^{\prime }\left(U\right)$ 的特征值，其具体形式为

$\alpha =\underset{x\in I}{\max }\left(\frac{1}{\Delta t^n}{\displaystyle {\int }_{t^n}^{t^{n+1}}\left(\left|u_{\tau }\left(x,t\right)\right|+\sqrt{g{h}_{\tau }\left(x,t\right)}\right)\text{d}t}\right)$ . (17)

2.4. 平衡良好的数值通量

请注意，在前一小节中提出的ADER-DG方法(45)不具有良好的平衡特性。在本节中，我们提出了对数值通量的修正，目的是准确地保持稳态(2)。

为了保持正性(如第2.14节所述)，我们还在时间步长 $t^n$ 上引入了新的网格边界值：

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

其中

$h_{\tau ,j+\frac{1}{2}}^{*,+}=h_{\tau ,j+\frac{1}{2}}^{*,-}$

在(2)中满足稳定性。

现在，让我们讨论平衡的数值通量。在(15)定义的Lax-Friedrichs数值通量 ${\hat{F}}_{j+\frac{1}{2}}$ 中，附加项 $-\alpha \left(U_{\tau ,j+\frac{1}{2}}^{n,+}-U_{\tau ,j+\frac{1}{2}}^{n,-}\right)$ 对数值黏度起到了作用，这是非线性守恒律所必需的。然而，它们可能破坏稳态的良好平衡特性。因此，我们将该通量(15)修改为

$\begin{array}{l}f\left(U_{\tau ,j+\frac{1}{2}}^{-}\left(t\right),U_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\right)=\frac{1}{2}[\left(\begin{array}{c}{h}_{\tau ,j+\frac{1}{2}}^{*,+}{u}_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\\ \frac{{\left({\left(hu\right)}_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\right)}^2}{u_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)}+\frac{1}{2}g{\left(h_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\right)}^2\end{array}\right)+\left(\begin{array}{c}{h}_{\tau ,j+\frac{1}{2}}^{*,-}{u}_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)\\ \frac{{\left({\left(hu\right)}_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)\right)}^2}{u_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)}+\frac{1}{2}g{\left(h_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)\right)}^2\end{array}\right)\\ \begin{array}{c}\stackrel{}{}\\ \\ \\ \end{array}-\alpha \left(\left(\begin{array}{c}{h}_{\tau ,j+\frac{1}{2}}^{*,+}\\ {\left(hu\right)}_{\tau ,j+\frac{1}{2}}^{+}\left(t\right)\end{array}\right)-\left(\begin{array}{c}{h}_{\tau ,j+\frac{1}{2}}^{*,-}\\ {\left(hu\right)}_{\tau ,j+\frac{1}{2}}^{-}\left(t\right)\end{array}\right)\right)].\end{array}$ (19)

这里我们将通量项的第一个分量从 ${\left(hu\right)}_{\tau ,j+\frac{1}{2}}^{\pm }\left(t\right)$ 修改为 $h_{\tau ,j+\frac{1}{2}}^{*,\pm }\left(t\right)u_{\tau ,j+\frac{1}{2}}^{*,\pm }\left(t\right)$ 以便在第2.14节中解释保正性。

2.5. 源项近似

接下来，我们给出源项积分的高阶平衡近似，其形式为

${\displaystyle \underset{{\Omega }_j}{\iint }{S}^2\left(t\right)v\text{d}x\text{d}t}=-g{\displaystyle \underset{{\Omega }_j}{\iint }\left(h+b\right)\left(t\right)b_{x}v\text{d}x\text{d}t}+\frac{1}{2}g{\displaystyle \underset{{\Omega }_j}{\iint }{b}_x^{2}v\text{d}x\text{d}t},$ (20)

其中 $S^2$ 表示源项的第二分量。我们可以进一步把这个积分分解为

$\begin{array}{l}{\displaystyle \underset{{\Omega }_j}{\iint }{S}^2\left(t\right)v\text{d}x\text{d}t}=-g{\displaystyle \underset{{\Omega }_j}{\iint }\left[\left(h+b\right)\left(t\right)-\overline{{\left(h+b\right)}_j}\left(t\right)\right]b_{x}v\text{d}x\text{d}t}\\ -g{\displaystyle \underset{{\Omega }_j}{\iint }\overline{{\left(h+b\right)}_j}\left(t\right)b_{x}v\text{d}x\text{d}t}+\frac{1}{2}g{\displaystyle \underset{{\Omega }_j}{\iint }{b}_x^{2}v\text{d}x\text{d}t}\\ =-g{\displaystyle \underset{{\Omega }_j}{\iint }\left[\left(h+b\right)\left(t\right)-\overline{{\left(h+b\right)}_j}\left(t\right)\right]b_{x}v\text{d}x\text{d}t}\\ -g{\displaystyle {\int }_{t^n}^{t^{n+1}}\overline{{\left(h+b\right)}_j}\left(t\right)\text{d}t\cdot \left(\left(bv\right)\left(x_{j+\frac{1}{2}}^{-}\right)-\left(bv\right)\left(x_{j-\frac{1}{2}}^{+}\right)-{\displaystyle {\int }_{I_j}b{v}_x\text{d}x}\right)}\\ +\frac{1}{2}g{\displaystyle {\int }_{t^n}^{t^{n+1}}\text{d}t\cdot \left(\left(bv\right)\left(x_{j+\frac{1}{2}}^{-}\right)-\left(bv\right)\left(x_{j-\frac{1}{2}}^{+}\right)-{\displaystyle {\int }_{I_j}b{v}_x\text{d}x}\right)},\end{array}$ (21)

其中项 $b,b^2,\left(h+b\right)\left(t\right),{\overline{\left(h+b\right)}}_j\left(t\right)$ 分别用它们的数值近似值 $b_{\tau },b_{\tau }^2,\overline{{\left(h+b\right)}_{\tau ,j}}\left(t\right)$ 代替。此外，将 $b_{\tau },b_{\tau }^2$ 的边界值替换为它们在网格面处的平均值，记为 $\left\{b_{\tau }\right\},\left\{b_{\tau }^2\right\}$ ，与稳态时的数值通量(51)一致。

2.6. 斜率限制器

对于双曲守恒律系统，当解包含间断点时，ADER-DG方法通常与斜率限制程序相结合，以压缩间断点附近可能出现的振荡。文献中有许多不同的斜率限制器选择，本文我们考虑经典的全变分有界(TVB)限制器 [50] 。

标准限制器可能与良好平衡性质相冲突，我们根据 [51] 中提出的思想提出了一个良好平衡限制器程序。未知 $U_{\tau }$ 上的标准TVB限制器包括两个步骤。第一步是根据网格平均值 $\overline{U_{\tau ,j}},\overline{U_{\tau ,j\pm 1}}$ 和网格边界值 $U_{\tau ,j+1/2}^{-},U_{\tau ,j-1/2}^{+}$ 来检查单元 $I_j$ 中是否需要任何限制。第二步是对该单元 $I_j$ 中的变量 $U_{\tau }$ 应用TVB限制器。在提出的平衡斜率限制程序中，我们将首先根据单元平均值 $\overline{{\left(h+b\right)}_{\tau ,j}},\overline{{\left(h+b\right)}_{\tau ,j\pm 1}}$ 和 ${\left(h+b\right)}_{\tau ,j+1/2}^{-},{\left(h+b\right)}_{\tau ,j-1/2}^{+}$ 来检查是否需要限制。如果单元 $I_j$ 被标记为需要限制，则实际的TVB限制器像往常一样应用于 $U_{\tau }$ 。注意 ${\left(h+b\right)}_{\tau }$ 在静水稳态(2)时变为常数，因此在达到稳态时不施加限位器，保持了良好的平衡性质。

2.7. 均衡格式概要

我们提出的浅水波方程(1)的良好平衡ADER-DG方法由式(45)给出，其中的数值通量定义在式(51)中，源项近似在上文中已经给出。综合前面各小节的结果，可以很直观地证明如下结果：

命题1 上述浅水波方程(1)的ADER-DG方法可以保持稳态解(2)的well-balanced性。

证明 在稳态(2)时，我们有

${\left(h+b\right)}_{\tau }=C,u_{\tau }=0.$

我们只需要证明包含源项的动量方程的平衡性。将方程中的通量项用 $F^{\left(2\right)}$ 表示。对于式(51)中给出的数值通量，式(45)左侧的数值近似为

$\begin{array}{l}{\displaystyle {\int }_{t^n}^{t^{n+1}}{\hat{F}}^{\left(2\right)}\left(U_{j+\frac{1}{2}}\left(t\right)\right)\text{d}t}\cdot v_{j+\frac{1}{2}}^{-}-{\displaystyle {\int }_{t^n}^{t^{n+1}}{\hat{F}}^{\left(2\right)}\left(U_{j-\frac{1}{2}}\left(t\right)\right)\text{d}t}\cdot v_{j-\frac{1}{2}}^{+}-{\displaystyle {\iint }_{{\Omega }_j}{F}^{\left(2\right)}\left(U\left(t\right)\right)v_x\text{d}x\text{d}t}\\ ={\displaystyle {\int }_{t^n}^{t^{n+1}}\frac{1}{2}g\left\{h_{\tau ,j+\frac{1}{2}}^2\right\}\left(t\right)\text{d}t\cdot v\left(x_{j+\frac{1}{2}}^{-}\right)}-{\displaystyle {\int }_{t^n}^{t^{n+1}}\frac{1}{2}g\left\{h_{\tau ,j-\frac{1}{2}}^2\right\}\left(t\right)\text{d}t\cdot v\left(x_{j-\frac{1}{2}}^{+}\right)}-{\displaystyle {\iint }_{{\Omega }_j}\frac{1}{2}g{h}_{\tau }^2\left(t\right)v_x\text{d}x\text{d}t.}\end{array}$ (22)

由于 ${\left(h+b\right)}_{\tau }\left(t\right)=\overline{{\left(h+b\right)}_{\tau ,j}}\left(t\right)\equiv C$ ，源项近似变成

$\begin{array}{l}{\displaystyle {\iint }_{{\Omega }_j}{S}_{\tau }^{\left(2\right)}\left(t\right)v\text{d}x\text{d}t}=-gC{\displaystyle {\int }_{t^n}^{t^{n+1}}\text{d}t}\cdot \left({\left\{b_{\tau }\right\}}_{j+\frac{1}{2}}{v}_{j+\frac{1}{2}}^{-}-{\left\{b_{\tau }\right\}}_{j-\frac{1}{2}}{v}_{j-\frac{1}{2}}^{+}-{\displaystyle {\int }_{I_j}{b}_{\tau }{v}_x\text{d}x}\right)\\ +\frac{1}{2}g{\displaystyle {\int }_{t^n}^{t^{n+1}}\text{d}t}\cdot \left({\left\{b_{\tau }^2\right\}}_{j+\frac{1}{2}}{v}_{j+\frac{1}{2}}^{-}-{\left\{b_{\tau }^2\right\}}_{j-\frac{1}{2}}{v}_{j-\frac{1}{2}}^{+}-{\displaystyle {\int }_{I_j}{b}_{\tau }^2{v}_x\text{d}x}\right).\end{array}$

注意这个等式

${\left(\frac{1}{2}g{h}_{\tau }^2\left(t\right)\right)}_x={\left(\frac{1}{2}g{\left(h+b\right)}_{\tau }^2\left(t\right)-g{\left(h+b\right)}_{\tau }\left(t\right)b_{\tau }+\frac{1}{2}g{b}_{\tau }^2\right)}_x=-gC{\left(b_{\tau }\right)}_x+\frac{1}{2}g{\left(b_{\tau }^2\right)}_x.$

并由此我们得出结论，通量和源项近似在静水稳态时相互平衡，从而得到期望的良好平衡特性。

2.8. Positivity-Preserving方法

在 [52] 中，Zhang和Shu提出了设计双曲守恒律的高阶最大原理保持方法的框架。从那时起，该方法得到了许多关注，并被广泛应用，包括 [53] 中的浅水波方程。结果表明，在合适的CFL条件下，能保持水高的非负性而不影响质量守恒，并能保持通解的高阶精度。在这里，我们将探讨这种方法在第2节中为浅水波方程提出的平衡良好的ADER-DG方法中的应用。

如 [54] 所述，实现这一目标的关键组成部分是以下两项：该方法一阶版本的正性，以及与高阶方法耦合的简单保正限制器。在建立之后，我们只考虑简单的欧拉时间离散化，同样的结果可以推广到多步和TVD高阶龙格–库塔方法。考虑平衡良好的ADER-DG方法(45)与(51)中定义的数值通量，取测试函数 $v=1$ 导致湿截面的单元平均更新如下

${\bar{h}}_j^{n+1}={\bar{h}}_j^n-{\lambda }_1\left[{\hat{F}}_{j+\frac{1}{2}}^{\left(1\right)}-{\hat{F}}_{j-\frac{1}{2}}^{\left(1\right)}\right],$ (23)

其中

${\hat{F}}_{j+\frac{1}{2}}^{\left(1\right)}=\frac{1}{2}{\displaystyle {\int }_{t^n}^{t^{n+1}}\left[h_{j+\frac{1}{2}}^{*,-}{u}_{j+\frac{1}{2}}^{-}\left(t\right)+h_{j+\frac{1}{2}}^{*,+}{u}_{j+\frac{1}{2}}^{+}\left(t\right)-\alpha \left(h_{j+\frac{1}{2}}^{*,+}-h_{j+\frac{1}{2}}^{*,-}\right)\right]\text{d}t.}$ (24)

并且有 ${\lambda }_1=1/\Delta x$ 。为方便起见，设 $\lambda ={\lambda }_1\Delta t^n$ 。为了不丧失一般性，我们在本节中忽略下标 $\tau$ 。该方法的一阶形式为

$h_j^{n+1}=h_j^n-{\lambda }_1\left[{\hat{F}}_{j+\frac{1}{2}}^{\left(1\right)}-{\hat{F}}_{j-\frac{1}{2}}^{\left(1\right)}\right].$ (25)

其数值通量定义为

${\hat{F}}_{j+\frac{1}{2}}^{\left(1\right)}=\frac{1}{2}\Delta t^n\left[h_{j+\frac{1}{2}}^{*,-}{u}_{j+\frac{1}{2}}^{-}\left(t\right)+h_{j+\frac{1}{2}}^{*,+}{u}_{j+\frac{1}{2}}^{+}\left(t\right)-\alpha \left(h_{j+\frac{1}{2}}^{*,+}-h_{j+\frac{1}{2}}^{*,-}\right)\right],$ (26)

其中 $\alpha =\max \left(\left|u\right|+\sqrt{gh}\right)$ 。

关于这种一阶方法的保正性，我们有如下引理：

引理2 在CFL条件 $\lambda \alpha \le 1$ 下，考虑具有数值通量(57)的一阶格式(56)。如果 $h_j^n,h_{j\pm 1}^n$ 是非负的，那么 $h_j^{n+1}$ 也是非负的。

证明 格式(56)可写成

${\displaystyle \underset{r=1}{\overset{N}{\sum }}{\hat{w}}_r=1}$

$\begin{array}{l}{h}_j^{n+1}=\left[1-\frac{1}{2}\lambda \left(\alpha +u_j^n\right)\frac{h_{j+\frac{1}{2}}^{*,-}}{h_j^n}-\frac{1}{2}\lambda \left(\alpha -u_j^n\right)\frac{h_{j-\frac{1}{2}}^{*,+}}{h_j^n}\right]h_j^n+\frac{1}{2}\lambda \left(\alpha +u_{j-1}^n\right)h_{j-\frac{1}{2}}^{*,-}+\frac{1}{2}\lambda \left(\alpha -u_{j+1}^n\right)h_{j+\frac{1}{2}}^{*,+}\\ \ge \left(1-\lambda \alpha \right)h_j^n+\frac{1}{2}\lambda \left(\alpha +u_{j-1}^n\right)h_{j-\frac{1}{2}}^{*,-}+\frac{1}{2}\lambda \left(\alpha -u_{j+1}^n\right)h_{j+\frac{1}{2}}^{*,+}.\end{array}$

因此， $h_j^{n+1}$ 是 $h_{j-1}^n$ ， $h_j^n$ 和 $h_{j+1}^n$ 的线性组合，并且所有系数都是非负的，因为 $0\le h_{j+\frac{1}{2}}^{*,-}$ ， $h_{j-\frac{1}{2}}^{*,+}\le h_j^n$ 。因此， $h_j^{n+1}\ge 0$ 。

接下来，我们讨论高阶格式。我们参考 [54] 了解细节，这里只介绍主要思想。我们在区间 $I_j$ 上引入N点( $2N-3\ge k$ ) Legendre gaas-lobatto正交规则，用 $S_j=\left\{x_{j-\frac{1}{2}}={\hat{x}}_j^1,{\hat{x}}_j^2,\cdots ,{\hat{x}}_j^{N-1},{\hat{x}}_j^N=x_{j+\frac{1}{2}}\right\}$ ，表示这些正交点，在区间 $\left[-1/2,1/2\right]$ 上对应的正交权值 ${\hat{w}}_r$ 满足。由于求积对于k次多项式是精确的，我们有

${\bar{h}}_j^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{h}_j^n\text{d}x}={\displaystyle \underset{r=1}{\overset{N}{\sum }}{\hat{w}}_r{h}_j^n\left({\hat{x}}_j^r\right)}={\displaystyle \underset{r=2}{\overset{N-1}{\sum }}{\hat{w}}_r{h}_j^n\left({\hat{x}}_j^r\right)}+{\hat{w}}_r{h}_{j-\frac{1}{2}}^{+}+{\hat{w}}_N{h}_{j+\frac{1}{2}}^{-}.$ (27)

按照 [55] 的方法，我们得到的结果是：

命题3 对于浅水波方程(1)的格式(45)，我们设 $h_j^n\left(x\right)$ 为单元格 $I_j$ 中的DG多项式。如果 $h_{j-\frac{1}{2}}^{+},h_{j-\frac{1}{2}}^{-}$ 和 $h_j^n\left({\hat{x}}_j^r\right)\left(r=2,\cdots ,N-1\right)$ 都是非负的，则在CFL条件下 ${\bar{h}}_j^{n+1}$ 也是非负的：

$\lambda \alpha \le {\hat{w}}_1.$ (28)

证明我们把(58)代入(54)，重写(54)通过加减 ${\displaystyle {\int }_{t^n}^{t^{n+1}}{f}^{\left(1\right)}\left(h_{j-\frac{1}{2}}^{*,+},u_{j-\frac{1}{2}}^{+}\left(t\right),h_{j+\frac{1}{2}}^{*,-},u_{j+\frac{1}{2}}^{-}\left(t\right)\right)\text{d}t}$ ：

${\bar{h}}_j^{n+1}={\displaystyle \underset{r=2}{\overset{N-1}{\sum }}{\hat{w}}_r{h}_j^n\left({\hat{x}}_j^r\right)}+{\hat{w}}_1{h}_1^{+}+{\hat{w}}_N{h}_N^{-},$ (29)

以及

${\hat{h}}_1=h_{j-\frac{1}{2}}^{+}-\frac{{\lambda }_1}{{\hat{w}}_r}{\displaystyle {\int }_{t^n}^{t^{n+1}}\left[f^{\left(1\right)}\left(h_{j-\frac{1}{2}}^{*,+},u_{j-\frac{1}{2}}^{+}\left(t\right),h_{j+\frac{1}{2}}^{*,-},u_{j+\frac{1}{2}}^{-}\left(t\right)\right)-f^{\left(1\right)}\left(h_{j-\frac{1}{2}}^{*,-},u_{j-\frac{1}{2}}^{-}\left(t\right),h_{j-\frac{1}{2}}^{*,+},u_{j-\frac{1}{2}}^{+}\left(t\right)\right)\right]\text{d}t,}$ (30)

${\hat{h}}_N=h_{j+\frac{1}{2}}^{+}-\frac{{\lambda }_1}{{\hat{w}}_N}{\displaystyle {\int }_{t^n}^{t^{n+1}}\left[f^{\left(1\right)}\left(h_{j+\frac{1}{2}}^{*,-},u_{j+\frac{1}{2}}^{-}\left(t\right),h_{j+\frac{1}{2}}^{*,+},u_{j+\frac{1}{2}}^{+}\left(t\right)\right)-f^{\left(1\right)}\left(h_{j-\frac{1}{2}}^{*,+},u_{j-\frac{1}{2}}^{+}\left(t\right),h_{j+\frac{1}{2}}^{*,-},u_{j+\frac{1}{2}}^{-}\left(t\right)\right)\right]\text{d}t.}$ (31)

注意到

$\begin{array}{c}{\hat{h}}_1=h_{j-\frac{1}{2}}^{+}-\frac{{\lambda }_1}{2{\hat{w}}_1}{\displaystyle {\int }_{t^n}^{t^{n+1}}[h_{j-\frac{1}{2}}^{*,+}{u}_{j-\frac{1}{2}}^{+}\left(t\right)+h_{j+\frac{1}{2}}^{*,-}{u}_{j+\frac{1}{2}}^{-}\left(t\right)-\alpha \left(h_{j+\frac{1}{2}}^{*,-}-h_{j-\frac{1}{2}}^{*,+}\right)}\\ -\left(h_{j-\frac{1}{2}}^{*,-}{u}_{j-\frac{1}{2}}^{-}\left(t\right)+h_{j-\frac{1}{2}}^{*,+}{u}_{j-\frac{1}{2}}^{+}\left(t\right)-\alpha \left(h_{j-\frac{1}{2}}^{*,+}-h_{j+\frac{1}{2}}^{*,-}\right)\right)]\text{d}t\\ =\left(1-\frac{\lambda \alpha }{{\hat{w}}_1}\frac{h_{j-\frac{1}{2}}^{*,+}}{h_{j-\frac{1}{2}}^{+}}\right)h_{j-\frac{1}{2}}^{+}+\frac{\lambda }{2{\hat{w}}_1}\left(\alpha +\frac{1}{\Delta t^n}{\displaystyle {\int }_{t^n}^{t^{n+1}}{u}_{j+\frac{1}{2}}^{-}\left(t\right)\text{d}t}\right)h_{j+\frac{1}{2}}^{*,-}+\frac{\lambda }{2{\hat{w}}_1}\left(\alpha -\frac{1}{\Delta t^n}{\displaystyle {\int }_{t^n}^{t^{n+1}}{u}_{j-\frac{1}{2}}^{-}\left(t\right)\text{d}t}\right)h_{j-\frac{1}{2}}^{*,-}\\ \ge \left(1-\frac{\lambda \alpha }{{\hat{w}}_1}\right)h_{j-\frac{1}{2}}^{+}+\frac{\lambda }{2{\hat{w}}_1}\left(\alpha +\frac{1}{\Delta t^n}{\displaystyle {\int }_{t^n}^{t^{n+1}}{u}_{j+\frac{1}{2}}^{-}\left(t\right)\text{d}t}\right)h_{j+\frac{1}{2}}^{*,-}+\frac{\lambda }{2{\hat{w}}_1}\left(\alpha -\frac{1}{\Delta t^n}{\displaystyle {\int }_{t^n}^{t^{n+1}}{u}_{j-\frac{1}{2}}^{-}\left(t\right)\text{d}t}\right)h_{j-\frac{1}{2}}^{*,-}.\end{array}$ (32)

可以清楚地看到 ${\hat{h}}_1$ 是 $h_{j-\frac{1}{2}}^{+},h_{j+\frac{1}{2}}^{*,-}$ 和 $h_{j-\frac{1}{2}}^{*,-}$ 的线性组合，并且所有系数都是非负的，因为 $0\le h_{j-\frac{1}{2}}^{*,+}\le h_{j-\frac{1}{2}}^{+}$ (这也证明了我们在 $\frac{\lambda \alpha }{{\hat{w}}_1}\le 1$ 的合适CFL条件下分别定义α和为(49)和(50)的原因。因此，我们得到 ${\hat{h}}_1\ge 0$ ，同样有 ${\hat{h}}_N\ge 0$ 。由于 ${\hat{w}}_1={\hat{w}}_N$ ，这两个CFL条件相同，也就是(59)。因此， ${\bar{h}}_j^{n+1}\ge 0$ ，因为它是 ${\hat{h}}_1,{\hat{h}}_N$ 和 $h_j^n\left({\hat{x}}_j^r\right)\left(r=2,\cdots ,N-1\right)$ 的凸组合。

我们在DG多项式 $U_j^n\left(x\right)={\left(h_j^n\left(x\right),{\left(hu\right)}_j^n\left(x\right)\right)}^{\text{T}}$ 上采用了以下保正极限器 [56] [57] ，

${\tilde{U}}_j^n\left(x\right)=\theta \left(U_j^n\left(x\right)-{\bar{U}}_j^n\right)+{\bar{U}}_j^n,\theta =\min \left\{1,\frac{{\bar{h}}_j^n}{{\bar{h}}_j^n-m_j}\right\},$ (33)

并有

$m_j=\underset{x\in S_j}{\min }{h}_j^n\left(x\right)=\underset{r=1,\cdots ,N}{\min }{h}_j^n\left({\hat{x}}_j^r\right).$ (34)

很容易观察到 ${\tilde{h}}_j^n\left({\hat{x}}_j^r\right)\ge 0\left(r=1,\cdots ,N\right)$ ，该极限保持变量 $U_j^n\left(x\right)$ 的局部守恒。

我们计算修正的多项式 ${\tilde{U}}_j^n\left(x\right)$ ，并在均衡方法中使用 ${\tilde{U}}_j^n\left(x\right)$ 代替 $U_j^n\left(x\right)$ (45)。根据 [56] [57] 中的证明，我们可以验证在CFL条件下，与该保正限制器耦合的平衡方法是高阶精度、保正和质量守恒的

$\lambda \alpha \le {\hat{w}}_1.$ (35)

为了提高效率，只有当对下一个时间步长的初步计算产生负湿截面时，我们才实现时间步长限制(66)。

3. 数值结果

为了证明well-balanced的ADER-DG格式的良好性能，我们给出了几个标准的浅水波方程的数值算例的结果。在所有的算例中，CFL数和引力常数g分别取0.16和9.812。

3.1. 平坦底部的黎曼问题

在本节中，我们考虑包含平坦底部的干面积(即 $b\left(x\right)\equiv 0$ )的黎曼问题。这些例子在 [58] 中使用过，这

里选择它们是为了证明我们的方法的保正能力。该例的计算域设置为[−300, 300]，初始条件为

$hu\left(x,0\right)=0\text{and}u\left(x,0\right)=\{\begin{array}{l}10\text{if}x\le 0,\\ 0\text{otherwise}.\end{array}$

左侧为静水高度设为10，右侧为干区。该问题的解析解可在 [59] 中找到。我们用我们的平衡保正方法计算了这个问题，该方法具有简单的传输边界条件和200个均匀单元。 $t=4,8,12$ 时刻的解如图1所示。我们还在这些图中绘制了精确的解，以便进行比较。我们可以观察到，数值结果很好地反映了精确解。

3.2. 抛物线碗

Sampson等人 [60] 推导出了具有抛物型底地形的一维浅水波方程的解析解。这为我们的数值方法提供了一个很好的测试用例。这个例子已在 [61] 中用于具有摩擦源项的浅水波方程。我们取抛物线底部：

$b\left(x\right)=h_0{\left(x/a\right)}^2,$

常数 $h_0$ 和a下文定义。计算域设置为[−5000, 5000]。对于不含摩擦源项的浅水波方程，解析水面由式给出

$h\left(x,t\right)+b\left(x\right)=h_0-\frac{B^2}{4g}\cos \left(2\omega t\right)-\frac{B^2}{4g}-\frac{Bx}{2a}\sqrt{\frac{8g_0}{g}}\cos \left(\omega t\right),$

其中 $\omega =\sqrt{2g{h}_0}/a,B=\text{constant}$ 。干湿锋面的确切位置如图2所示：

$x_0=-\frac{B\omega a^2}{2g{h}_0}\cos \left(\omega t\right)\pm a.$

3.3. 干燥河床的生成问题

我们首先考虑 [62] 中提出的平底 $b\left(x\right)=0$ 的数值实验。初始条件为

$h\left(x,0\right)=0.1,hu\left(x,0\right)=\{\begin{array}{l}-0.3\text{if}x\le 5,\\ 0.3\text{otherwise},\end{array}x\in \left[0,10\right].$

在这个例子中，干燥的河床会在两个方向相反的稀有波的中间产生。干层的生成使得这个例子很难用数值计算。将t = 1时刻得到的结果进行比较，正如图3可以观察到的那样，目前的格式在两个稀疏波之间留下了一个小的湿区。通过图中绘制出的精确解我们可以清晰的观察到，数值结果具有良好的保正性，并且很好地反映了精确解。

4. 结论

本文设计了一种有效的保正高阶ADER-DG方法，用于求解具有非平坦底地形的一维浅水波方程。所提出的方法具有两个良好的特点：对静水稳态解具有良好的平衡性和在干湿锋附近保持正性。通过一种新颖的源项分裂和每个分裂源项的适当的良好平衡近似来实现良好平衡特性。采用 [63] 中的一个简单的保正限幅器，以确保所得方法保持截面湿面积的非负性。我们进行了相应的数值模拟实验，结果表明本文所提出的方法具有良好的平衡性，在稳态解附近的小扰动测试中有效，在干湿前沿附近保持正性以及高阶精度，并且在连续解和间断解上都表现良好。

基金项目

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