非守恒双曲方程组的路径守恒ADER间断Galerkin方法：在浅水方程中的应用

A Path-Conservative ADER Discontinuous Galerkin Method for Non-Conservative Hyperbolic Equations: Applications in Shallow Water Equations

摘要: 本文提出了求解非守恒双曲型偏微分方程的一种新的路径守恒间断Galerkin (DG)方法。特别地，这里的方法采用了一级ADER (在空间和时间的任意导数)方法来实现时间离散化。此外，该方法采用微分变换(DT)过程而不是Cauchy-Kowalewski (C-K)过程来实现局部时间演化。与经典的ADER方法相比，该方法不需要求解内部单元的广义黎曼问题。与RKDG (Runge-Kutta DG)方法相比，该方法不需要中间步骤，因此需要较少的计算机存储空间。简而言之，当前的方法是一步一步完全离散的。而且，该方法在空间和时间上都容易获得高阶精度。浅水方程的数值结果表明，该方法具有较高的阶精度，对间断解具有较好的分辨率。

Abstract: In this article, we propose a new path-conservative discontinuous Galerkin (DG) method to solve the non-conservative hyperbolic partial differential equations (PDE). In particular, the method here applies the one-stage ADER (Arbitrary DERivatives in space and time) approach to fulfill the tem-poral discretization. In addition, this method uses the differential transformation (DT) procedure other than the Cauchy-Kowalewski (C-K) procedure to achieve the local temporal evolution. Com-pared with the classical ADER methods, the current method is free of solving generalized Riemann problems at inter-cells. In comparison with the Runge-Kutta DG (RKDG) methods, the proposed method needs less computer storage thanks to no intermediate stages. In brief, this current method is one-step, one-stage, and fully-discrete. Moreover, this method can easily obtain arbitrary high-order accuracy in space and in time. Numerical results for shallow water equations (SWEs) show that the method enjoys high-order accuracy and keeps good resolutions for discontinuous so-lutions.

1. 引言

物理学和工程领域的许多流体问题都可以根据第一原理描述为守恒定律

$\frac{\partial W}{\partial t}+\nabla \cdot \underset{\_}{\underset{\_}{F}}\left(W\right)=0,$ (1)

其中 $W$ 为守恒变量向量， $\underset{\_}{\underset{\_}{F}}\left(W\right)=\left(f\left(W\right),g\left(W\right)\right)$ 为物理通量张量。到目前为止，我们可以通过守恒原

理深刻地理解自然界中的大多数物理运动。然而，在空气动力学、天体物理学、航空航天、石油工业等领域对可压缩多相流/介质流的建模中，由于不同相(介质)之间复杂的相互作用，出现了非守恒积项(即未知解的空间导数)。因此，相关的数学模型不能以守恒形式表示。然而，上述流体问题可以用拟线性非守恒双曲方程组表示如下

$\frac{\partial W}{\partial t}+\underset{\_}{\underset{\_}{{\rm A}}}\left(W\right)\cdot \nabla W=\text{0},$ (2)

其中 $\underset{\_}{\underset{\_}{{\rm A}}}\left(W\right)=\left(A\left(W\right),B\left(W\right)\right)$ 为方程组矩阵。这里使用分块矩阵语法给出矩阵 $A\left(W\right)$ 和 $B\left(W\right)$ 的紧凑符号。其中，假设方程(2)为双曲型，即 $A\left(W\right)$ 和 $B\left(W\right)$ 分别具有m个实特征值和m个线性无关的特征向量的完整集合。特别地，如果 $A\left(W\right)$ 和 $B\left(W\right)$ 是 $f\left(W\right)$ 和 $g\left(W\right)$ 的雅可比矩阵，则方程(2)将简化为守恒定

律(1)的双曲方程。这一观察结果表明，方程(2)适合同时表示守恒定律和非守恒方程组。

然而，方程(2)的主要问题在于间断情况下弱解的经典定义的不足。直到Dal Maso、Le Floch和Murat的理论(DLM理论) [1] 的出现，这一重大问题才有了很大的进展。DLM理论给出了方程组(2)弱解的定义，通过引入路径 $\Psi =\Psi \left(\xi ;W^{-},W^{+}\right)$ 连接相空间中的两个状态 $W^{-}$ 和 $W^{+}$ 。随后，针对DLM理论 [1] ，Castro [2] 和Pares [3] 根据非守恒双曲方程组发展了路径守恒方法(2)。实际上，路径守恒方法 [2] [3] [4] 也可以被称为Toumi [5] 对Roe方法弱表述的扩展。

后来，继最初的成果 [1] [2] 之后，许多研究者对路径守恒方法进行了许多尝试。代表性研究主要有：ADER格式 [6] [7] [8] 、FORCE格式 [9] [10] 、HLLC Riemann解算方法 [11] [12] 、Osher Riemann解算方法 [13] [14] 、中心格式 [15] 、中心迎风格式 [4] [16] 、ADER-DG方法 [17] 等。关于最新进展和简要的历史回顾，我们参考文献 [18] 和 [19] 。

本研究的主要目的是针对非守恒双曲型方程组提出一种新的路径守恒DG方法(2)。该方法采用单阶段ADER方法实现高阶时间离散化，因此称为路径守恒ADER-DG方法。其基本思想是采用DT过程 [20] [21] [22] 代替C-K过程，通过低阶空间展开系数来表示解的时间展开系数。此外，DT过程可以实现任意高阶精度的局部时间演化。就我们而言，这将是第一次尝试将DG方法与DT过程一起应用于非守恒方程组。然后，我们将该方法推广到以非守恒形式处理SWEs。特别是Li等 [17] 提出了一种以双曲平衡律形式求解SWEs的ADER-DG方法。此外，Li [23] 等人将ADER-DG方法推广到求解气体动力学中的欧拉方程。在此，对于非保守形式的SWEs，ADER-DG方法的成功将说明ADER-DG方法具有通用性。

本文的结构如下：第2节阐述了路径守恒的ADER-DG方法的一般框架，然后将所得方法以非守恒形式应用于一维浅水波方程。第3节实现了典型算例来验证所提出方法的性能。第4节给出了一些结论。

2. 路径守恒的ADER-DG法的一般公式

本文针对非守恒方程组(2)给出了一种一般的ADER-DG方法框架。实际上，我们首先将方程组(2)在

给定的近似空间中乘以一个测试函数 $\varphi \left(\stackrel{\to }{x}\right)$ ，然后在一个时间单元 $C_j\times \left[t^n,t^{n+1}\right)$ 上积分，得到

${\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{C_j}{\overset{}{\int }}\frac{\partial W}{\partial t}\varphi \text{d}V\text{d}t}}+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{C_j}{\overset{}{\int }}\left(\underset{\_}{\underset{\_}{{\rm A}}}\left(W\right)\cdot \nabla W\right)\varphi \text{d}V\text{d}t}}=0,$

在这里 $C_j$ 是空间单元。然后，我们用 $W_{\tau }\left(\stackrel{\to }{x},t\right)$ 作为时间多项式近似 $W\left(\stackrel{\to }{x},t\right)$。由于 $W_{\tau }$ 通常在单元间出现跳跃，我们提出了一种路径守恒方法来处理跳跃。

接下来，我们得到PDE(6)的DG方法：对于任意测试函数 $\varphi$ ，数值解 $W_{\tau }$ 满足以下等式

$\begin{array}{l}{\displaystyle \underset{C_j}{\int }{W}_{\tau }\left(\stackrel{\to }{x},t^{n+1}\right)\varphi \text{d}V}-{\displaystyle \underset{C_j}{\int }{W}_{\tau }\left(\stackrel{\to }{x},t^n\right)\varphi \text{d}V}+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{C_j\backslash \partial C_j}{\overset{}{\int }}\left(\underset{\_}{\underset{\_}{{\rm A}}}\left(W_{\tau }\right)\cdot \nabla W_{\tau }\right)\varphi \text{d}V\text{d}t}}\\ +{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{\partial C_j}{\overset{}{\int }}{D}^{-}\left(W_{\tau }^{-},W_{\tau }^{+},\stackrel{\to }{n}\right)\varphi \text{d}S\text{d}t}}=0,\end{array}$ (3)

其中 $\stackrel{\to }{n}$ 为单元 $C_j$ 边界处的单位向外法向量。

DG方法(3)采用ADER方法实现时间离散化。因此，DG方法(3)是一层一层完全离散的。因此，我

们将提出的方法(3)称为ADER-DG方法。此外，为了实现时间离散，我们需要通过DT过程实现 $W_{\tau }\left(\stackrel{\to }{x},t^n\right)$ 在时间单元 $C_j\times \left[t^n,t^{n+1}\right)$ 中的局部时间演化。DT过程将在第2.1.3节中详细描述。这里， $W_{\tau }^{-}$ 来自单元内部 $C_j$ 和 $W_{\tau }^{+}$ 来自单元 $C_j$ 外部。此外，符号 $D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)$ 表示单元间的跳跃项，满足以下要求 [2] [3] [6] ：

· 对于每一个 $W_{\tau }$ 和 $\stackrel{\to }{n}$ ，

$D^{-}\left(W_{\tau },W_{\tau }\right)=0.$

· 对于每一个 $W_{\tau }^{-}$ ， $W_{\tau }^{+}$ 和 $\stackrel{\to }{n}$ ，

$D^{-}\left(W_{\tau }^{-},W_{\tau }^{+},\stackrel{\to }{n}\right)+D^{-}\left(W_{\tau }^{-},W_{\tau }^{+},-\stackrel{\to }{n}\right)={\displaystyle {\int }_0^1\left(\underset{\_}{\underset{\_}{{\rm A}}}\left(\Psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)\stackrel{\to }{n}\right)\frac{\partial \Psi }{\partial \xi }\text{d}\xi }.$

另外， $\Psi =\Psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)$ 表示连接状态 $W_{\tau }^{-}$ 和 $W_{\tau }^{+}$ 的足够光滑的路径。满足上述条件的方法称为路

径守恒方法。

此外，根据 [17] [23] ，所提出的求解守恒定律(1)双曲方程组的ADER-DG方法(3)的守恒等价如下

${\displaystyle \underset{C_j}{\int }{W}_{\tau }\left(\stackrel{\to }{x},t^{n+1}\right)\varphi }\text{d}V-{\displaystyle \underset{C_j}{\int }{W}_{\tau }\left(\stackrel{\to }{x},t^n\right)\varphi }\text{d}V+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{C_j\backslash \partial C_j}{\overset{}{\int }}\underset{\_}{\underset{\_}{F}}\left(W_{\tau }\right)\cdot \nabla }}\varphi \text{d}V\text{d}t+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{F}_{j+\frac{1}{2}}}\varphi \text{d}S\text{d}t=0,$ (4)

以 $F_{j+\frac{1}{2}}:=F_{j+\frac{1}{2}}\left(W_{\tau }^{-},W_{\tau }^{+},\stackrel{\to }{n}\right)$ 为数值通量近似单元间物理通量 $\underset{\_}{\underset{\_}{F}}$ 。实际上，当非守恒方程组(2)满足守恒

律(1)，即A和B分别是f和g的雅可比矩阵时，上述非守恒形式的路径守恒式ADER-DG方法(3)可简化为保守的ADER-DG方法(4)。Dumbser等人在有限体积格式的框架下对这种等价性进行了详细的证明 [6] 。

2.1. 一维SWEs的应用

基于路径保守的ADER-DG方法的一般框架，我们取以下一维SWEs为例

$\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}$ (5)

来说明该方法(3)的具体实现步骤。式中， $h\left(x,t\right)$ 、 $u\left(x,t\right)$ 分别为水深、流体速度。符号 $b\left(x\right)$ 代表底部地形，字母“g”代表重力加速度。

按照 [2] 中的方法，我们将(5)中的几何源项合并到项 $A\left(W\right)\frac{\partial W}{\partial x}$ 中，得到如下的非守恒形式

$\frac{\partial W}{\partial t}+A\left(W\right)\frac{\partial W}{\partial x}=0,$ (6)

在这里

$W=\left(\begin{array}{l}h\\ hu\\ b\end{array}\right),A\left(W\right)=\left(\begin{array}{ccc}0& 1& 0\\ c^2-u^2& 2u& c^2\\ 0& 0& 0\end{array}\right),$

用 $c=\sqrt{gh}$ 表示声速。这种操作使方程组(6)更易于获得Well-balanced方法。

从数学的观点来看，非守恒方程组(6)保留稳态解，满足

$A\left(W\right)\frac{\partial W}{\partial x}=0.$

特别地，静止水域的稳态解有以下几种形式

$h+b=\text{Constantand}u=0.$

传统的方法不能准确地保持这种稳定状态，并导致非物理振荡。Well-balanced (W-B)方法 [24] [25] 可以在离散水平上保持达到机器精度的稳态，并且即使在相对粗糙的网格上也可以解决稳态的小扰动 [26] ，从而相应地提高计算效率。

2.1.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)$ 和 ${\tau }_j=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$ 作为网格中心和单元 $I_j$ 的大小。最大网格尺寸定义为 $\tau =\underset{1\le j\le N}{\max }{\tau }_j$ 。这里，我们应用 ${\Omega }_j=I_j\times \left[t^n,t^{n+1}\right)$ 作为时间单元，并令

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

作为近似空间，其中 $P^k\left({\Omega }_j\right)$ 表示时间单元 ${\Omega }_j$ 上的一组时间多项式，其阶数最高为k。

2.1.2. 一维路径保守的ADER-DG方法的构建

对于以为方程组(6)，ADER-DG方法(3)如下：对于 $\forall \varphi \left(x\right)\in V_{\tau }^k$ ，解 $W_{\tau }\left(x,t\right)$ 满足下列等式

$\begin{array}{l}{\displaystyle \underset{I_j}{\int }{W}_{\tau }\left(x,t^{n+1}\right)\varphi \left(x\right)\text{d}x}-{\displaystyle \underset{I_j}{\int }{W}_{\tau }\left(x,t^n\right)\varphi \left(x\right)\text{d}x}+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{I_j\backslash \partial I_j}{\int }\left(A\left(W_{\tau }\right)\frac{\partial W_{\tau }}{\partial x}\right)\varphi \left(x\right)\text{d}x\text{d}t}}\\ +{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{\partial I_j}{\int }{D}^{-}\left(W_{\tau }^{-},W_{\tau }^{+},\stackrel{\to }{n}\right)\varphi \left(x\right)\text{d}S\text{d}t}}=0,\text{for}j=1,2,3,\cdots ,N.\end{array}$ (8)

进而，方程(8)也有如下的等式形式

$\begin{array}{l}{\displaystyle \underset{I_j}{\int }{W}_{\tau }\left(x,t^{n+1}\right)\varphi \left(x\right)\text{d}x}-{\displaystyle \underset{I_j}{\int }{W}_{\tau }\left(x,t^n\right)\varphi \left(x\right)\text{d}x}+{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{\displaystyle \underset{I_j\backslash \partial I_j}{\int }\left(A\left(W_{\tau }\right)\frac{\partial W_{\tau }}{\partial x}\right)\varphi \left(x\right)\text{d}x\text{d}t}}\\ +{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{D}^{-}\left(W_{j+\frac{1}{2}}^{-},W_{j+\frac{1}{2}}^{+}\right)\varphi \left(x\right)\text{d}t}-{\displaystyle \underset{t^n}{\overset{t^{n+1}}{\int }}{D}^{-}\left(W_{j-\frac{1}{2}}^{-},W_{j-\frac{1}{2}}^{+}\right)\varphi \left(x\right)\text{d}t}=0,\text{for}j=1,2,3,\cdots ,N.\end{array}$ (9)

对于单元间的跳转项 $D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)$ 有不同的选择，例如

· Osher跳跃项：

$D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)=\frac{1}{2}{\displaystyle {\int }_0^1\left(A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)-\left|A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)\right|\right)}\frac{\partial \psi }{\partial \xi }\text{d}\xi .$ (10)

· Roe跳跃项：

$D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)={\left({\tilde{A}}_{\psi }\left(W_{\tau }^{-},W_{\tau }^{+}\right)\right)}^{-}\cdot \left(W_{\tau }^{+}-W_{\tau }^{-}\right),$ (11)

其中， ${\tilde{A}}_{\psi }\left(W_{\tau }^{-},W_{\tau }^{+}\right)$ 为Toumi等在 [5] 中定义的某种意义上 $A\left(W_{\tau }\right)$ 的Roe线性化矩阵，即函数 ${\tilde{A}}_{\psi }$ 满足以

下性质：

- 对每个 $W_{\tau }^{-},W_{\tau }^{+}$ ，矩阵 ${\tilde{A}}_{\psi }\left(W_{\tau }^{-},W_{\tau }^{+}\right)$ 具有m个不同的实特征值

${\lambda }_1<{\lambda }_2<\cdots <{\lambda }_m.$

- 兼容性属性

${\tilde{A}}_{\psi }\left(W_{\tau },W_{\tau }\right)=A\left(W_{\tau }\right),\text{foreach}{W}_{\tau }.$

- 对于任意 $W_{\tau }^{-},W_{\tau }^{+}$ ，矩阵 ${\tilde{A}}_{\psi }$ 满足下式要求

${\tilde{A}}_{\psi }\left(W_{\tau }^{-},W_{\tau }^{+}\right)\cdot \left(W_{\tau }^{+}-W_{\tau }^{-}\right)={\displaystyle {\int }_0^{1}A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)}\frac{\partial \psi }{\partial \xi }\text{d}\xi ,$

根据广义Roe的性质。

另外，对于式(10)和式(11)中矩阵的绝对值算子，我们通常使用以下符号

$\Lambda =diag\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m\right),$

$\left|\Lambda \right|=diag\left(\left|{\lambda }_1\right|,\left|{\lambda }_2\right|,\cdots ,\left|{\lambda }_m\right|\right),$

${\Lambda }^{-}=\frac{1}{2}\left(\Lambda -\left|\Lambda \right|\right),$

$A^{-}=R{\Lambda }^{-}{R}^{-1},$

其中 $R$ 是矩阵 $A$ 的右特征向量的矩阵， $R^{-1}$ 表示它的逆矩阵。

目前为止，可以将当前方法(9)视为给定路径 $\psi \left(\xi \right)$ 的函数，形式如下

$\psi \left(\xi \right)=\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right),\text{0}\le \xi \le \text{1}\text{.}$

并且，函数 $\psi \left(\xi \right)=\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)$ 是Lipschitz连续的，满足一定的正则性和相容条件

$\psi \left(0;W_{\tau }^{-},W_{\tau }^{+}\right)=W_{\tau }^{-},$

$\psi \left(1;W_{\tau }^{-},W_{\tau }^{+}\right)=W_{\tau }^{+},$

$\psi \left(\xi ;W_{\tau },W_{\tau }\right)=W_{\tau }.$

在本文中，我们应用简单的分段路径

$\psi \left(\xi \right)=\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)=W_{\tau }^{-}+\xi \left(W_{\tau }^{+}-W_{\tau }^{-}\right),\text{0}\le \xi \le 1,$ (12)

比如 [2] [3] [6] 。

这种路径的选择十分有用的，它保证了前文所提出的方法对于SWEs [2] [3] [6] 是W-B的。此外，由路径(12)，Osher跳跃项(10)可简化为以下形式

$D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)=\frac{1}{2}{\displaystyle {\int }_0^1\left(A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)-\left|A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)\right|\right)\text{d}\xi }\cdot \left(W_{\tau }^{+}-W_{\tau }^{-}\right).$ (13)

同时，Roe跳跃项(11)可简化为以下形式

$D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)=\frac{1}{2}\left({\displaystyle {\int }_0^{1}A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)\text{d}\xi }-\left|{\displaystyle {\int }_0^{1}A\left(\psi \left(\xi ;W_{\tau }^{-},W_{\tau }^{+}\right)\right)\text{d}\xi }\right|\right)\cdot \left(W_{\tau }^{+}-W_{\tau }^{-}\right).$ (14)

对于式(13)和式(14)中路径积分的计算，我们采用了具有适当高阶精度的高斯正交规则。

2.1.3. DT算法

为了构建方法(9)，我们需要提前实现从 $t^n$ 开始的时间单元的局部时间演化。这个操作的原因是我们需要在每个时间单元中以多项式的形式得到数值解。然后，我们可以在(9)中以高阶精度计算时间积分和时间积分。实际上，为了实现这一目标，ADER方法 [27] [28] [29] [30] [31] 使用C-K过程反复微分控制PDE，并通过空间导数得到时间导数。为了获得高阶时间精度，我们需要高阶时间导数。此时，由于链式规则的使用，C-K程序将变得非常繁琐。Dumbser和Munz [31] 利用Leibnize规则提出了一种高效的算法。最近，Dumbser等人 [6] [7] ，Tang等人 [32] 应用局部DG预测方法 [29] 来代替C-K程序。最近，Li等人通过DT程序开发了一种针对SWEs的ADER-DG方法 [17] 。

在本研究中，我们采用DT程序而不是C-K程序。实际上，DT过程最初是针对非线性初值问题发展起来的 [33] [34] 。随后，Ayaz将DT过程推广到二维情况 [20] 以及方程组情况 [21] 。Kurnaza等 [22] 将其推广到更一般的n维情况。此外，Norman和Finkel [33] 应用该方法构建了一维SWEs的多矩有限体积格式。

下面我们给出DT过程的具体定义，如 [20] [21] [22] 所示。假设函数 $q\left(x,t\right)$ 在单元格 $I_j\text{at}{t}^n$ 是已知的，则DT定义如下

$\tilde{q}\left(k_x,k_t\right)={\frac{1}{k_x!k_t!}\frac{{\partial }^{k_x+k_t}u\left(x,t\right)}{\partial x^{k_x}\partial t^{k_t}}|}_{x=x_j,t=t^n},\text{with}q\left(x,t\right)={\displaystyle \underset{k_t=0}{\overset{k}{\sum }}{\displaystyle \underset{k_x=0}{\overset{k-k_t}{\sum }}\tilde{q}\left(k_x,k_t\right){\left(x-x_j\right)}^{k_x}}{\left(t-t^n\right)}^{k_t}}.$ (15)

其中， $\tilde{q}\left(k_x,k_t\right)$ 表示根据原函数 $q\left(x,t\right)$ 变换后的函数。

实际上， $\tilde{q}\left(k_x,k_t\right)$ 表示关于 $q\left(x,t\right)$ 的展开系数用截断泰勒级数的形式表示。表1显示了这里使用的

一些转换后的函数。

Table 1. Transformed functions of some functions

表1. 部分函数的变换函数

然后，具体说明了DT过程的实现步骤。一开始，我们得到

$b_r\left(x\right)={\displaystyle \underset{k_x=0}{\overset{k}{\sum }}\tilde{b}\left(k_x\right){\left(x-x_j\right)}^{k_x}},$

$h_{\tau }\left(x,0\right)={\displaystyle \underset{k_x=0}{\overset{k}{\sum }}\tilde{h}\left(k_x,0\right){\left(x-x_j\right)}^{k_x}},$

${\left(hu\right)}_{\tau }\left(x,0\right)={\displaystyle \underset{k_x=0}{\overset{k}{\sum }}\left(\widetilde{hu}\right)\left(k_x,0\right){\left(x-x_j\right)}^{k_x}},$

使用 $L^2$ 投影来近似单元格 $I_j$ 中的 $b\left(x\right),h\left(x,0\right)$ 以及 $hu\left(x,0\right)$ 。

随后，我们对式(5)的两端施加DT过程，得到如下递归式

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

在这里，我们使用以下辅助变量

${\tilde{G}}_1\left(k_x,k_t\right)={\displaystyle \underset{s=0}{\overset{k_t}{\sum }}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\tilde{G}}_{1,a}}}\left(r,s\right){\tilde{G}}_{1,b}\left(k_x-r\right)\left(k_t-s\right),$

${\tilde{G}}_2\left(k_x,k_t\right)=\frac{g}{2}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}{\displaystyle \underset{r=0}{\overset{k_x}{\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{s=0}{\overset{k_t}{\sum }}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}\left(\widetilde{hu}\right)\left(r,s\right)\left(\widetilde{hu}\right)}}\left(k_x-r,k_x-s\right),$

${\tilde{G}}_{1,b}\left(k_x,k_t\right)=\frac{1}{\tilde{h}\left(0,0\right)}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\tilde{G}}_{1,b}\left(k_x-r,k_t-s\right)\tilde{h}}}\left(r,s\right),$

${\tilde{G}}_{3,a}\left(k_x\right)=g\left(k_x+1\right)\tilde{b}\left(k_x+1\right).$

然后，将 $\tilde{h}\left(k_x,0\right)$ ， $\left(\widetilde{hu}\right)\left(k_x,0\right)$ 和 $\tilde{b}\left(k_x\right)$ ， $k_x=0,1,\cdots ,k$ 带入(16)，即可递归得到

$\tilde{h}\left(k_x,k_t\right),\left(\widetilde{hu}\right)\left(k_x,k_t\right),\text{for}{k}_t=0,1,\cdots ,k;k_x=0,1,\cdots ,k-k_t,$

结果就是 $W_{\tau }\left(x,t\right):={\left(h_{\tau }\left(x,t\right),{\left(hu\right)}_{\tau }\left(x,t\right),b_{\tau }\left(x\right)\right)}^{\text{T}}={\displaystyle \underset{k_t=0}{\overset{k}{\sum }}{\displaystyle \underset{k_x=0}{\overset{k-k_t}{\sum }}\tilde{W}\left(k_x,k_t\right){\left(x-x_i\right)}^{k_x}{\left(t-t^n\right)}^{k_t}}\in V_{\tau }^k}$ 。

每个空间单元 ${\Omega }_j$ 中的 $\tilde{W}\left(k_x,k_t\right):={\left(\tilde{h},\left(\widetilde{hu}\right),\tilde{b}\right)}^{\text{T}}\left(k_x,k_t\right)$ 。此外，为了更好地理解这一过程，我们

给出了DT过程的详细算法。

备注1. 总而言之，DT过程的关键功能是根据现有解 $W_{\tau }\left(x,t^n\right)$ 为每个空间单元局部提供高阶时间演化。

备注2. C-K过程直接使用控制PDE的符号展开，并且需要重新计算许多项。因此，C-K过程根据复杂度导致指数增长。然而，DT过程相对简单，但复杂性是可以预测的。

备注3. 在实践中，不需要根据底部地形b应用DT过程，因为底部b只依赖于空间变量x。

2.1.4. 斜率限制器

一般来说，对于间断问题，斜率限制器是必不可少的。在这里，我们使用总变分界(TVB)限制器 [35]

来控制非物理振荡。事实上，我们只根据数值解 $U_{\tau }:={\left(h_{\tau },{\left(hu\right)}_{\tau }\right)}^{\text{T}}$ 来实现TVB限制器步骤，不包括与时间t无关的底部形貌 $b_{\tau }\left(x\right)$ 。具体来说，我们需要在单元平均的基础上识别出“坏单元”(即包含间断的单元)从 $W_{\tau }\left(x,t^{n+1}\right)$ 在 $t^{n+1}$ 上的 ${\bar{U}}_{\tau ,j}\left(t^{n+1}\right),{\bar{U}}_{\tau ,j\pm 1}\left(t^{n+1}\right)$ 和单元平均值 ${\bar{U}}_{\tau ,j}\left(x_{j+\frac{1}{2}}^{\pm },t^{n+1}\right)$ 。

为了便于说明，我们首先给出一些符号

$\begin{array}{l}{\tilde{U}}_j=U_{\tau ,j+\frac{1}{2}}^{-}\left(t^{n+1}\right)-{\bar{U}}_{\tau ,j}\left(t^{n+1}\right),\\ {\widetilde{\stackrel{˜}{U}}}_j=U_{\tau ,j}^{-}\left(t^{n+1}\right)-{\bar{U}}^{+}{}_{\tau ,j-\frac{1}{2}}\left(t^{n+1}\right),\\ {\Delta }_{+}{\bar{U}}_j=U_{\tau ,j+1}^{-}\left(t^{n+1}\right)-{\bar{U}}_{\tau ,j}\left(t^{n+1}\right),\\ {\Delta }_{-}{\bar{U}}_j=U_{\tau ,j}^{-}\left(t^{n+1}\right)-{\bar{U}}_{\tau ,j-1}\left(t^{n+1}\right).\end{array}$ (17)

然后，我们得到下面更新后的值

$\begin{array}{l}{\tilde{U}}_j^{\left(\mathrm{mod}\right)}=m\left({\tilde{U}}_j,{\Delta }_{+}{\bar{U}}_j,{\Delta }_{-}{\bar{U}}_j\right),\\ {\widetilde{\stackrel{˜}{U}}}_j^{\left(\mathrm{mod}\right)}=m\left({\widetilde{\stackrel{˜}{U}}}_j,{\Delta }_{+}{\bar{U}}_j,{\Delta }_{-}{\bar{U}}_j\right),\end{array}$

对(17)中的变量使用TVB限制器 [35] 。这里， $m\left(a_1,a_2\text{,}{a}_3\right)$ 是一个最小模函数

$m\left(a_1,a_2,a_3\right)=\{\begin{array}{l}{a}_1,\text{if}\left|a_1\right|\le M\Delta x^2,\\ s\cdot \underset{1\le i\le 3}{\min }\left|a_i\right|,\text{if}\left|a_1\right|>M\Delta x^2\text{and}s=sign\left(a_1\right)=sign\left(a_2\right)=sign\left(a_3\right),\\ 0,\text{otherwise},\end{array}$

$M>0$ 。单元 $I_j$ 被识别为坏单元，只要

${\tilde{U}}_j^{\left(\mathrm{mod}\right)}\ne {\tilde{U}}_j\text{or}{\widetilde{\stackrel{˜}{U}}}_j^{\left(\mathrm{mod}\right)}\ne {\widetilde{\stackrel{˜}{U}}}_j.$

随后，单元极限值 $U_{\tau }\left(x_{j+\frac{1}{2}}^{\pm },t_{n+1}\right)$ 被定义为

$\begin{array}{l}{U}_{\tau }^{\left(\mathrm{mod}\right)}\left(x_{j+\frac{1}{2}}^{-},t_{n+1}\right)={\bar{U}}_j+{\tilde{U}}_j^{\left(\mathrm{mod}\right)},\\ U_{\tau }^{\left(\mathrm{mod}\right)}\left(x_{j-\frac{1}{2}}^{+},t_{n+1}\right)={\bar{U}}_j-{\widetilde{\stackrel{˜}{U}}}_j^{\left(\mathrm{mod}\right)}.\end{array}$ (18)

最后，根据有限的单元值(18)和单元平均值确定多项式 $U_{\tau }\left(x,t^{n+1}\right)$ 。

2.1.5. 一维路径守恒的ADER-DG方法的实施步骤

对于以为方程组(6)，在一个时间区间 $\left[t^n,t^{n+1}\right)$ 内，本文方法的具体步骤如下：

1) 最初，由 $W\left(x,0\right)$ 得到

$W_{\tau }\left(x,0\right)={\displaystyle \underset{k_x=0}{\overset{k}{\sum }}\tilde{W}\left(k_x,0\right){\left(x-x_j\right)}^{k_x}}\text{in}{I}_j,\text{for}j=1,2,\cdots ,N$ .

2) 使用递归步(16)，根据 $\tilde{W}\left(k_x,0\right)$ 从 $W_{\tau }\left(x,t^n\right)$ 在 $t^n$ 处得到 $\tilde{W}\left(k_x,k_t\right)$ ，并在每个 ${\Omega }_j$ 上得到 $W_{\tau }\left(x,t\right),\text{for}j=1,2,\cdots ,N$ 。

3) 根据(13)和(14)构造跳跃项 $D^{-}\left(W_{\tau }^{-},W_{\tau }^{+}\right)$ 。

4) 用一级公式(9)更新 $W_{\tau }\left(x,t^{n+1}\right)$ 。

5) 需要时，在 $W_{\tau }\left(x,t^{n+1}\right)$ 上使用TVB斜率限制器。

6) 重复步骤2)~5)。

3. 数值结果

在此，我们用几个典型的算例来证明所提出的方法。在所有的算例中，我们都在近似空间中使用了二阶以上的时间多项式(即 $k=2$ )。为了保证数值的稳定性，我们取Courant-Friedrichs-Levy (CFL)数值为0.18。另外，我们设重力常数 $g=9.812$ 。为了节省空间，由于Roe类型方法保持了W-B性质，因此我们只使用带有Roe类型跳跃项的方法来呈现数值结果。

3.1. W-B性能测试

首先，我们用一个算例数值证明了W-B性质，如 [36] 。初始条件为

$u=0\text{and}h+b=10,x\in \left[0,10\right].$

数值证明了Well-balanced性质。特别是，我们在以下两种不同的底部上进行计算：

$b\left(x\right)=5\exp \left(-0.4{\left(x-5\right)}^2\right)\text{and}b\left(x\right)=\{\begin{array}{l}4\text{if}4\le x\le \text{8,}\\ 0\text{else}\text{.}\end{array}$

其中第一个底部是光滑的，而第二个底部是间断的。

表2和表3分别给出了在 $t=0.5$ 时根据两种不同的底部地形产生的误差。很明显，所有的误差都在机器精度的同一量级，故而预期的W-B性能是相应实现的。

Table 2. Errors according to the example over the first bottom

表2. 根据第一个底部的算例得出的误差

Table 3. Errors according to the example over the second bottom

表3. 根据第二个底部的算例得出的误差

3.2. 精度测试

接下来，我们用 [36] 中的一个算例，以及相应的底部地形和初始条件来证明其准确性。初始条件为

$h\left(x,0\right)=5+\exp \left(\cos \left(2\pi x\right)\right)\text{and}\left(hu\right)\left(x,0\right)=\sin \left(\cos \left(2\pi x\right)\right),x\in \left[0,1\right]$

在底部上 $b\left(x\right)={\sin }^2\left(\pi x\right)$ 。

由于解析解无法得到，我们首先需要在一个32,000单元的精细网格上计算这个例子直到t = 0.1 s，然后将得到的数值解作为参考。随后，我们计算t = 0.1 s时的L¹误差，并根据参考解得到精度阶数，如表4所示。得到了明显的三阶精度。

Table 4. Errors and accuracy orders for h and hu

表4. h和hu的误差和精度阶

3.3. 稳态水流的扰动

在这里，我们应用 [37] 中的一个算例来测试所提出的方法在凹凸形状底部地形捕捉小扰动的能力。初始数据为

$h\left(x,0\right)=\{\begin{array}{l}1-b\left(x\right)+\epsilon \text{if}\text{1}\text{.1}\le x\le 1.2,\\ 1-b\left(x\right)\text{else,}\end{array}\text{and}u\left(x,0\right)=0,$

其中 $\epsilon >0$ 作为凹凸的参数

$b\left(x\right)=\{\begin{array}{l}0.25\left(\cos \left(10\pi \left(x-1.5\right)\right)+1\right)\text{if}\text{1}\text{.4}\le x\le 1.6,\\ 0\text{else,}\end{array}x\in \left[0,2\right].$

随着时间的推移，初始扰动分解为两个脉冲，沿两个不同的方向运动，如图1所示，与文献 [37] [38] 中的数值结果保持高度一致。显然，这两种类型的脉冲都可以很好地分解。

Figure 1. Surface level (left) and water discharge (right) at t = 0.2

图1. t = 0.2时的地表水位(左)和水量(右)

3.4. 矩形底部的溃坝问题

接下来，我们处理一个溃坝问题 [38] [39] [40] ，并使用以下初始条件

在矩形底部上方有

$b\left(x\right)=\{\begin{array}{l}8\text{if}\left|x-750\right|\le \text{1500}/8\text{,}\\ 0\text{else,}\end{array}x\in \left[0,1500\right].$

图2分别给出了 $t=15\text{s}$ 和 $t=60\text{s}$ 时的数值解，表明该方法得到了很好的数值解果，与文献 [38] [39] [40] 中的数值解保持一致。

Figure 2. Surface level at t = 15 (left) and t = 60 (right)

图2. t = 15 s (左)和t = 60 s (右)时的表面水平

3.5. 经过驼峰的稳定水流

进而，我们使用一个广泛应用的算例来验证该方法 [41] 。实际上，这个算例是基于初始条件对跨临界和亚临界流进行建模的

$u\left(x,0\right)=0\text{and}h\left(x,0\right)=0.33,x\in \left[0,25\right]$

越过驼峰

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

随后，我们在空间区间的两端施加不同的边界条件，计算到t = 200 s。此外，我们还展示了精确解 [42] 加以比较，从而更直观地验证我们的方法。

· 算例A：无激波的跨临界流动

上游边界的流量 $hu=1.53{\text{m}}^{\text{2}}/\text{s}$ ；下游边界的水深 $h=0.66\text{m}$ 。图3中的Case A的结果表明，数值结果同实际结果相吻合。

· 算例B：带激波的跨临界流动

上游边界的流量 $hu=0.18{\text{m}}^{\text{2}}/\text{s}$ ；下游边界的水深 $h=0.33\text{m}$ 。图3中的Case B的结果表明，数值结果与精确解有良好的一致性，无伪震荡。

· 算例C：亚临界流

上游边界的流量 $hu=0.18{\text{m}}^{\text{2}}/\text{s}$ ；下游边界的水深 $h=2\text{m}$ 。图3中的Case C的结果表明，数值解和精确解的拟合效果良好。

Figure 3. Surface level h + b at t = 200 s

图3. t = 200 s时的表面水位h + b

3.6. 跨过台阶的溃坝问题

这里，我们利用一个经过台阶的溃坝算例来测试我们的数值方法，以下算例可以很好地测试所构造的方法的鲁棒性。

· 算例A

首先，我们使用以下初始数据实现 [43] 中的算例

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(4,0\right)\text{if}x\le 0,\\ \left(1,0\right)\text{else,}\end{array}x\in \left[-10,10\right]$

跨过一个台阶状底部

$b\left(x\right)=\{\begin{array}{l}0\text{if}x\le 0,\\ 1\text{else}\text{.}\end{array}$

随着时间的推移，这个算例产生了向左移动的稀疏波向右移动的激波。

· 算例B

初始条件为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(4,5\right)\text{if}x\le 0,\\ \left(1,-0.9\right)\text{else,}\end{array}x\in \left[-10,10\right]$

在上面相同的底部。随着时间的推移，这个算例发展出两个朝着不同方向移动的冲击。

· 算例C

本算例来自 [13] ，初始数据为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(0.75,-9.49365\right)\text{if}x\le 0,\\ \left(1.10594,-4.94074\right)\text{else,}\end{array}x\in \left[-15,5\right]$

越过一个台阶

$b\left(x\right)=\{\begin{array}{l}0\text{if}x\le 0,\\ 0.2\text{else}\text{.}\end{array}$

· 算例D

初始数据为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(0.75,-1.35624\right)\text{if}x\le 0,\\ \left(1.10594,-4.94074\right)\text{else,}\end{array}x\in \left[-10,4\right]$

在与算例C相同的底部。

图4为同一网格200个单元在t = 1时的数值结果，与精确的结果吻合较好。我们可以清楚地发现数值结果和精确解保持高度一致，并拥有陡峭的间断过渡。

Figure 4. Surface level at t = 1

图4. t = 1时的表面水平

4. 结论

本文针对非保守双曲方程组，提出了一种新的DG方法。该方法采用一阶ADER方法实现时间离散化。然后，我们将该方法应用于求解SWEs。为了实现局部时间的高阶演化，本文采用DT过程代替C-K过程，由空间膨胀系数递归得到时间膨胀系数。与C-K程序相比，DT程序更简洁，编程更方便。此外，由于该方法不需要中间阶段，因此需要较少的计算机存储空间，并且不需要求解单元间的广义黎曼问题。由于该方法具有明确的一步性和紧凑的模板结构，因此它是超级计算机并行计算的理想选择。我们可以很容易地在空间和时间上进行任意高阶精度，而不需要太多的编码工作。基于DT程序。总之，所提出的方法是一步一步完全离散的。此外，大量的数值结果表明，该方法具有高阶精度、良好的平衡性和对间断解的良好分辨率。

基金项目

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

NOTES

^*通讯作者。