基于微分变换策略的全离散间断Galerkin方法

High Order Discontinuous Galerkin Methods with Differential Transform Strategy

摘要: 在本文中，针对一维双曲守恒律方程，我们基于微分变换策略提出了一种高阶全离散间断Galerkin方法。与传统间断Galerkin方法相比较，该方法的主要特点为存储量较小，在时间上能够到达任意高阶精度。与ADER方法相比较，避免了复杂的Cauchy-Kowalewski步骤，同时程序编写简单。数值结果表明该方法具有高阶精度，针对间断解保持高分辨率。

Abstract: In this paper, based on the differential transform strategy, we introduce high order full discrete discontinuous Galerkin methods for one-dimensional hyperbolic conservation laws. Compared with the standard discontinuous Galerkin methods, the current methods enjoy the following advantages including low storage as well as keeping arbitrary high order accuracy in time. In comparison with the ADER methods, the resulting methods avoid the complicated Cauchy-Kowalewski procedure and the coding is easy at the same time. In addition, numerical results also indicate that the resulting methods enjoy genuine high order accuracy for smooth solutions, and keep steep discontinuity transition.

1. 引言

在本文中，考虑双曲守恒律

$\{\begin{array}{l}{u}_t+f{\left(u\right)}_x=0\\ u\left(x,0\right)=u_0(x)\end{array}$

的高阶数值方法。间断Galerkin (DG)方法作为一种高阶方法首先由Reed和Hill [1] 于1973年在求解中子输运问题时提出。上个世纪90年代，Cockburn，Chi-Wang Shu等人 [2] [3] [4] [5] 提出了基于Runge-Kutta时间离散方法的Runge-Kutta DG (RKDG)方法，并将该方法成功应用于双曲守恒律方程。DG方法使用分段间断多项式空间作为试探和检验函数空间，并结合了有限元方法和有限体积格式的优点，特别易于处理复杂边界和边值问题，拥有灵活处理间断的能力，同时又适合自适应算法与并行策略，因此在计算流体力学领域得到了广泛应用 [6] - [17] 。关于DG方法的历史回顾与最新进展见文献 [18] [19] [20] 。传统的RKDG方法为显格式，其优点在于每个时间步长计算中不需要求解大规模线性代数方程组。但是，传统RKDG方法作为一种高精度方法，由于中间步的存在，导致存储量较高。此外，RKDG方法因为使用Gauss数值求积，所以计算成本非常高。

近年来以ADER方法 [21] 为代表的时–空任意高阶方法取得了长足进步。ADER方法是一类基于有限体积积分的具有高阶精度的数值方法。在ADER方法框架内，在间断面处对时间进行Taylor展开，人们通常使用Cauchu-Kowalewski步骤(又称为Lax-Wendroff步骤) [22] 将时间导数项转换为时间导数项，然后分别求解由展开式中的各个部分所构造的黎曼问题。但是Cauchu-Kowalewski步骤非常复杂，计算成本较高，同时编程复杂。

在本文中，针对双曲守恒律基于微分变换策略提出了高阶全离散DG方法。该DG方法是单步、显式的和全离散的，从而避免了存储量过高的缺点。此外微分变换策略较之Cauchu-Kowalewski过程更加简洁，程序编写简单。更进一步，因为数值解关于时间、空间皆表示为多项式，故可以使用精确积分，从而避免了计算成本较高的Gauss积分，节约了计算成本，因此更加适应实际问题的需求。

2. 微分变换简介

基于微分方程本身以得到方程解的时–空Taylor级数系数的递推关系，这种策略被称为微分变换(differential transformation) [23] 。与ADER方法中所用到的Cauchy-Kowalewski (C-K)技巧或者Lax-Wendforr技巧相比较，微分变换策略显得更简单，程序编写更加简洁。但是，并非所有方程都能够使用微分变换策略。目前，可以确定的是针对浅水波方程、欧拉方程可以使用微分变换策略 [24] 。

以一维Burgers方程的初值问题

$\{\begin{array}{l}{q}_t+{\left(\frac{1}{2}{q}^2\right)}_x=0\\ q\left(x,0\right)=q_0(x)\end{array}$

为例解释如何借助微分变换策略求解微分方程。

记Burgers方程为

$q_t+f{\left(q\right)}_x=0$ (1)

其中 $f\left(q\right)=\frac{1}{2}{q}^2:=f\left(x,t\right)$ 为流通量。首先得到 $Q\left(k,0\right)$ ，即把初始条件表示为 $q\left(x,0\right)={\displaystyle {\sum }_0^{\infty }Q\left(k,0\right){\left(x-x_j\right)}^k}$ 。然后，我们把未知函数 $q\left(x,t\right)$ 以及流通量 $f\left(x,t\right)$ 都展成时–空Taylor级数

$q\left(x,t\right)=\underset{k=0}{\overset{\infty }{{\displaystyle \sum }}}\underset{h=0}{\overset{\infty }{{\displaystyle \sum }}}Q\left(k,h\right){\left(x-x_j\right)}^k{\left(t-t^n\right)}^h,$

$f\left(x,t\right)=\underset{k=0}{\overset{\infty }{{\displaystyle \sum }}}\underset{h=0}{\overset{\infty }{{\displaystyle \sum }}}F\left(k,h\right){\left(x-x_j\right)}^k{\left(t-t^n\right)}^h,$ (2)

其中 $F\left(k,h\right)=\frac{1}{2}{\displaystyle {\sum }_{r=0}^h{\displaystyle {\sum }_{s=0}^{h}Q\left(r,s\right)Q\left(k-r,h-s\right)}}$ 。

将(2)中两个等式代入方程(1)，得到

$\underset{k=0}{\overset{\infty }{{\displaystyle \sum }}}\underset{h=0}{\overset{\infty }{{\displaystyle \sum }}}hQ\left(k,h\right){\left(x-x_j\right)}^k{\left(t-t^n\right)}^{h-1}=-\underset{k=1}{\overset{\propto }{{\displaystyle \sum }}}\underset{h=0}{\overset{\propto }{{\displaystyle \sum }}}kF\left(k,h\right){\left(x-x_j\right)}^{k-1}{\left(t-t^n\right)}^h$

调整下求和符号中的下标，得到

$\underset{k=0}{\overset{\infty }{{\displaystyle \sum }}}\underset{h=0}{\overset{\infty }{{\displaystyle \sum }}}\left(h+1\right)Q\left(k,h+1\right){\left(x-x_j\right)}^k{\left(t-t^n\right)}^h=-\underset{k=0}{\overset{\infty }{{\displaystyle \sum }}}\underset{h=0}{\overset{\infty }{{\displaystyle \sum }}}\left(k+1\right)F\left(k+1,h\right){\left(x-x_j\right)}^k{\left(t-t^n\right)}^h$

比较上式两端，得到如下递推关系

$Q\left(k,k+1\right)=-\frac{k+1}{h+1}F\left(k+1,h\right)$ (3)

借助于已有的 $Q\left(k,0\right)$ ，并利用递推公式(3)，可以得到解 $q\left(x,t\right)$ 的时–空Taylor级数展开式。

针对方程组，可以逐方程地借助微分变换策略来求解微分方程。

3. 一维问题的全离散DG方法

考虑非线性双曲守恒律方程

$U_t+F{\left(U\right)}_x=0$ (4)

首先将方程(4)两端同乘以检验函数 $v\left(x\right)\in v_h^k$ ，并在时–空单元 $I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right]\times \left[t^n,t^{n+1}\right]$ 上积分，然后利用分部积分得到：

$\begin{array}{l}{\displaystyle {\int }_{I_j}{U}^{n+1}v\text{d}x\text{d}t}-{\displaystyle {\int }_{I_j}{U}^{n}v\text{d}x\text{d}t}+{\displaystyle {\int }_{t^n}^{t^{n+1}}F\left(U\left(x_{j+\frac{1}{2},}t\right)\right)\text{d}t}\cdot v\left(x_{j+\frac{1}{2}}^{-}\right)\\ -{\displaystyle {\int }_{t^n}^{t^{n+1}}F\left(U\left(x_{j-\frac{1}{2}},t\right)\right)\text{d}t}\cdot v\left(x_{j-\frac{1}{2}}^{+}\right)-{\displaystyle {\iint }_{I_j}F\left(U\right)v_x\text{d}x\text{d}t}=0\end{array}$

对于积分 ${\displaystyle {\int }_{t^n}^{t^{n+1}}F\left(U\left(x_{j+\frac{1}{2}},t\right)\right)\text{d}t}={\displaystyle {\int }_{t^n}^{t^{n+1}}F\left(x_{j+\frac{1}{2}},t\right)\text{d}t}$ ，采用简单的Lax-Friedrichs数值流通量来逼近，即

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

这里 $\alpha =\max \left|\lambda \left(U\right)\right|$ ，其中 $\lambda \left(U\right)$ 为雅克比矩阵 $F\text{'}\left(U\right)$ 的特征值，并且最大值是在整个计算区域来取。

最后，得到全离散DG方法：对于 $\forall v\left(x\right)\in V_h^k$ ， $U_h\in V_h^k$ 满足

${\displaystyle {\int }_{I_j}{U}_h^{n+1}v\text{d}x\text{d}t}-{\displaystyle {\int }_{I_j}{U}_h^{n}v\text{d}x\text{d}t}+{\hat{F}}_{j+\frac{1}{2}}\cdot v\left(x_{j+\frac{1}{2}}^{-}\right)-{\hat{F}}_{j-\frac{1}{2}}\cdot v\left(x_{j-\frac{1}{2}}^{+}\right)-{\displaystyle {\iint }_{I_j}{F}_h{v}_x\text{d}x\text{d}t}=0$

由于 $U_h$ ， $F_n$ 均表示为关于x，t的多项式，所以对于积分采用精确积分，从而避免了计算成本较高的Gauss积分。

4. 数值结果

在本节中，通过数值实验来展示所提出的DG方法的良好性能。针对所有数值算例，我们分别将时间与空间取为二次多项式与四次多项式(即分别取 $\mathcal{K}=\mathcal{H}=2$ ，此时取 $\text{CFL}=0.18$ 。或者 $\mathcal{K}=\mathcal{H}=4$ ，同时取 $\text{CFL}=0.1$ )作为基函数。重力常数g取为9.812。

4.1. 精度测试

借助具有光滑解的算例来测试方法的精度阶。考虑线性对流方程：

$\{\begin{array}{l}{u}_t+u_x=0,0\le x\le 2\\ u\left(x,0\right)=\sin \left(\text{π}x\right),\text{periodic}\end{array}$

将该算例计算至 $t=2$ ，同时将方法的数值误差与精度阶分别展示于表1。从表中，可以很明显地发现该方法拥有预期的高阶精度。

4.2. 线性对流方程

借助带有间断解的算例来测试方法能否保持高分辨率

$\{\begin{array}{l}{u}_t+u_x=0,-1\le x\le 1\\ u\left(x,0\right)=u_0\left(x\right),\text{periodic}\end{array}$

Table 1. Third-order method error and precision order

表1. 三阶方法的误差和精度阶

其中

$u_0\left(x\right)=\{\begin{array}{l}\frac{1}{6}\left(G\left(x,\beta ,z-\delta \right)+G\left(x,\beta ,z+\delta \right)+4G\left(x,\beta ,\delta \right)\right),-0.8\le x\le -0.6\\ 1,-0.4\le x\le -0.2\\ 1-\left|10\left(x-0.1\right)\right|,0\le x\le 0.2\\ \frac{1}{6}\left(F\left(x,\alpha ,a-\Delta \right)+F\left(x,\alpha ,a+\Delta \right)+4F\left(x,\alpha ,a\right)\right),0.4\le x\le 0.6\\ 0,\text{otherwise}\end{array}$

这里，取 $G\left(x,\beta ,z\right)={\text{e}}^{-\beta {\left(x-z\right)}^2}$ ， $F\left(x,\alpha ,a\right)=\sqrt{\max \left(1-{\alpha }^2{\left(x-a\right)}^2,0\right)}$ 。参数取作 $a=0.5$ ， $z=-0.7$ ， $\delta =0.005$ ， $\alpha =10$ ， $\beta =\frac{\log 2}{36{\delta }^2}$ 。

基于带有100个单元的网格，将该算例计算至 $t=8$ ，同时将数值结果展示于图1。我们从图中发现数值结果能够同精确解保持很好的一致性，且保持高分辨率。

4.3. Burgers方程

接下来，测试该方法针对强间断问题能否保持高分辨率，考虑一维Burgers方程：

$u+{\left(\frac{u^2}{2}\right)}_x=0$

考虑初始条件

$u\left(x,0\right)=0.5+\sin \left(\text{π}x\right),x\in \left[0,2\right]$

该问题在 $t=1.5/\text{π}$ 时刻会产生激波。接下来考虑另一种初始条件

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

该初值问题会在左侧产生稀疏波，右侧产生激波。我们把以上两种初值问题的数值结果展示于图2并与精确解相比较。很显然数值结果保持陡峭的间断过渡，从而说明该方法对于间断解同样保持高分辨率。

4.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=0\end{array}$

为例来考虑方程组情形。

其中h，u分别表示水深，深度平均速度。g表示重力常数。

考虑浅水波方程中经典的Riemann问题，又称溃坝问题 [25] ，其初始条件具有以下情形：

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(h_L,u_L\right)\text{if}x\le x_0\\ \left(h_R,u_R\right)\text{otherwise}\end{array}$

其中 $x_0$ 为大坝位置， $\left(h_L,u_L\right)$ 与 $\left(h_R,u_R\right)$ 分别表示大坝左右两侧的水体深度与速度。

Figure 1. The numerical result of the linear convection equation at time t = 2. Grid n = 100

图1. 线性对流方程在t = 2时刻的数值结果。网格n = 100

Figure 2. Numerical results for the Burgers equation, grid n = 100. Left: t = 1.5/π, right: t = 0.5

图2. Burgers方程的数值结果，网格n = 100。左：t = 1.5/π，右：t = 0.5

我们考虑以下几种不同的初始条件：

4.4.1. 情形一

初始条件为

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(1,2.5\right)\text{if}x\le 10\\ \left(0.1,0\right)\text{otherwise}\end{array}$

计算区域为区间 $\left[0,50\right]$ 。计算该算例至 $t=7$ ，同时将计算结果展示于图3并与精确解相比较。数值结果表明数值结果与精确解保持一致，同时保持非振荡性质。

4.4.2. 情形二

初始条件如下：

$\left(h,u\right)\left(x,0\right)=\{\begin{array}{l}\left(1,-5\right)\text{if}x\le 25\\ \left(1,5\right)\text{otherwise}\end{array}$

计算区域为 $\left[0,50\right]$ 。将 $t=2.5$ 时刻的数值结果展示于图4。数值结果表明数值结果与精确解保持一致。

Figure 3. At the time t = 7, the result of the grid is n = 100. Water depth (left) and speed (right)

图3. 在t = 7时刻的数值结果，网格n = 100。水体深度(左)与速度(右)

Figure 4. At the time t = 2.5, the result of the grid is n = 100. Water depth (left) and speed (right)

图4. 在t = 2.5时刻的数值结果，网格n = 100。水体深度(左)与速度(右)

4.5. 欧拉方程

考虑空气动力学中的可压缩欧拉方程

$\{\begin{array}{l}{\rho }_t+{\left(\rho u\right)}_x=0\\ {\left(\rho u\right)}_t+{\left(\rho u^2+p\right)}_x=0\\ E_t+{\left(\left(E+p\right)u\right)}_x=0\end{array}$

其中ρ为密度，u为速度，E表示总能量，p是压力，且与总能量有如下关系

$E=\frac{p}{\gamma -1}+\frac{1}{2}\rho u^2$

其中γ为绝热比常数，在本文中取作 $\gamma =1.4$ 。

4.5.1. Sod问题

初始条件如下

Figure 5. The numerical result of Sod problem at t = 1.3, grid n = 200. From left to right: density, speed, pressure. nx = 200

图5. Sod问题在t = 1.3时刻的数值结果，网格n = 200。从左向右：密度，速度，压力。nx = 200

Figure 6. Lax problem numerical results at t = 1.3. Grid n = 200. From left to right: density, speed, pressure

图6. Lax问题在t = 1.3时刻的数值结果。网格n = 200。从左向右：密度，速度，压力

$\left(\rho ,u,p\right)\left(x,0\right)=\{\begin{array}{l}\left(1,0,1\right)\text{if}x\le 0\\ \left(0.5,0,0.571\right)\text{otherwise}\end{array}$

计算区域为[−5, 5]。将 $t=2$ 时刻的数值结果展示于图5。数值结果保持很好的分辨率。

4.5.2. Lax问题

初始条件如下

$\left(\rho ,u,p\right)\left(x,0\right)=\{\begin{array}{l}\left(0.445,0.698,3.528\right)\text{if}x\le 0\\ \left(0.5,0,0.571\right)\text{otherwise}\end{array}$

计算区域为[−5, 5]。将时刻 $t=1.3$ 的计算结果展示于图6。数值结果保持很好的间断过度。

5. 结论

在本文中，基于微分变换策略提出了一种高阶全离散间断Galerkin方法。该方法的存储量较小，程序编写简单。针对单个方程以及方程组，开展了广泛数值实验。数值结果表明该方法针对光滑问题具有高阶精度，同时针对间断解拥有高分辨率。

致谢

本研究得到国家自然科学基金项目(11771228)的资助。