波动方程的高阶间断有限元方法

High Order Discontinuous Finite Element Method for Wave Equation

摘要: 近年来，间断有限元方法在求解双曲型方程中获得广泛应用，与传统有限元和差分方法相比也有许多优势。在本文中，对双曲型的波动方程应用间断有限元分析，通过构造间断有限元的强形式，选择合适的数值通量，与精确解相比得到了较小的误差以及较高的收敛阶，达到了间断有限元的理论收敛阶。

Abstract: In recent years, discontinuous finite element method has been widely used in solving hyperbolic equations. Compared with traditional finite element method and difference method, discontinuous finite element method has many advantages. In this paper, the discontinuous finite element analysis is applied to the hyperbolic wave equation. By constructing the strong form of the discon-tinuous finite element and selecting the appropriate numerical flux, the smaller error and higher convergence order are obtained compared with the exact solution, and the theoretical convergence order of the discontinuous finite element is reached.

1. 引言

近几年，间断有限元方法(DG方法)在诸多领域获得广泛应用。DG方法最早是由Reed和Hill在求解中子输运方程时提出 [1]，之后其他学者进行了DG方法的误差分析并获得了最优收敛阶。经过数十年的发展，目前已出现了多种DG方法。这些DG方法大致可以分为两大类：基于数值通量形式的DG方法和基于内部罚函数的DG方法。第一类基于数值通量形式的DG方法主要用于求解双曲方程，以Cockburn和Shu提出的求解双曲守恒律方程的Runge-Kutta DG (RK-DG)有限元方法具有代表性 [2] [3]。Cockburn和Shu提出了几种常用的数值通量，适用于各种不同情况。在实际应用中，学者可以针对不同的方程构造独特的数值通量以满足实际情况 [4] [5]，和其他差分格式相比可以达到更高收敛阶 [6]，所以应用的更为广泛。

2. 波动方程及其间断有限元方法

2.1. 波动方程

以带有初值条件的一维无源声波方程为例，其方程形式如下：

$\{\begin{array}{l}\frac{{\partial }^{2}U\left(x,t\right)}{\partial t^2}-c^2\frac{{\partial }^{2}U\left(x,t\right)}{\partial x^2}=0x\in \left[L,R\right],t>0\\ U\left(x,0\right)=f\left(x\right)x\in \left[L,R\right]\\ \frac{\partial U\left(x,0\right)}{\partial t}=g\left(x\right)x\in \left[L,R\right]\end{array}$

首先对方程进行降阶处理，使其化为满足守恒律的形式。引入位移变量 $P=\left(x,t\right)$，满足如下的关系：

$c\frac{\partial U}{\partial x}=\frac{\partial P}{\partial t}.$

将其带入原声波方程中，并进行积分得到满足双曲守恒律的方程组：

$\{\begin{array}{l}\frac{{\partial }^{2}U}{\partial t^2}-c\frac{\partial }{\partial x}\left(\frac{\partial P}{\partial t}\right)=0\\ c\frac{\partial U}{\partial x}=\frac{\partial P}{\partial t}\end{array}\Rightarrow \{\begin{array}{l}\frac{\partial U}{\partial t}-c\frac{\partial P}{\partial x}=0\\ \frac{\partial P}{\partial t}\text{=}c\frac{\partial U}{\partial x}\end{array}\text{.}$

为简化方程，引入以下记号：则方程可以写成 $\frac{\partial W}{\partial t}+V\frac{\partial W}{\partial x}=0.$

$W=\left[\begin{array}{l}U\\ P\end{array}\right]V=\left[\begin{array}{cc}0& -c\\ -c& 0\end{array}\right]$

至此，通过引入新的变量将二阶声波方程转化为一阶方程组，并满足守恒律的要求 [7]。

2.2. 空间域间断有限元方法

对一阶方程组 $\frac{\partial W}{\partial t}+V\frac{\partial W}{\partial x}=0$ 进行空间离散。首先将求解区域 $\Omega =\left[L,R\right]$ 划分为K个不重叠的单元 $D^k=\left[x_l^k,x_r^k\right]$ 去逼近计算区域 $\Omega =\left[L,R\right]$，在单元上选取 $N_p$ 个节点。区域划分如图1所示：

Figure 1. Element division of calculation area

图1. 计算区域单元剖分

引进线性映射如下：将每个计算单元变换到区间[−1, 1]上。

$x\in D^k:x\left(r\right)=x_l^k+\frac{1+r}{2}{h}^k,h^k=x_r^k-x_l^k.$

在求解区域上构建试验函数空间： $V_h:=\left\{{\phi }_h\in L^2\left(\Omega \right):{{\phi }_h|}_K\in P^k\left(K\right)\right\}$，在单元上把局部解用 $N=N_p-1$ 阶多项式表示： $\left[\begin{array}{c}{E}_h^k\left(x,t\right)\\ H_h^k\left(x,t\right)\end{array}\right]={\displaystyle \underset{i=1}{\overset{N_p}{\sum }}\left[\begin{array}{c}{E}_h^k\left(x_i^k,t\right)\\ H_h^k\left(x_i^k,t\right)\end{array}\right]}{\phi }_i^k\left(x\right)$。

构造残量 ${\mathcal{R}}_h\left(x,t\right)=\frac{\partial W_h}{\partial t}+V\frac{\partial W_h}{\partial x}$，现要求残量正交于 $V_h$ 中所有的试验函数 [8]，即

${\displaystyle {\int }_{D^k}{R}_h}\left(x,t\right){\phi }_n\left(x\right)\text{d}x=0,1\le n\le N_p.$

经过一次分部积分得到DG方法弱形式如下，其中 $\hat{n}$ 表示局部外法向量，在一维问题中是一个标量。

${\displaystyle {\int }_{{\text{D}}^k}\left(\frac{\partial W_h^k}{\partial t}{\phi }_n-V{W}_h^k\frac{\text{d}{\phi }_n}{\text{d}x}\right)}\text{d}x=-{\left[V{W}_h^k{\phi }_n\right]}_{x_i^k}^{x_r^k}=-{\displaystyle {\int }_{\partial D^k}\hat{n}}\cdot V{W}_h^k{\phi }_n\text{d}x,1\le n\le N_p.$

在图1中，如果对试验函数和局部解没有约束条件，在单元界面上就可能会出现多个值 $W_R^{k-1}、W_L^k$，必须选择一个合适的解来进行接下来的计算。这就是数值通量的选择，目前先将其记为 ${\left(W_h\right)}^{*}$ [6]。

通过再一次分部积分，可以得到DG方法的强形式如下，这是计算常用的形式。

${\displaystyle {\int }_{D^k}{\mathcal{R}}_h}\left(x,t\right){\phi }_n\text{d}x={\displaystyle {\int }_{\partial D^k}\hat{n}}\cdot \left(V{W}_h^k-{\left(V{W}_h\right)}^{*}\right){\phi }_n\text{d}x,1\le n\le N_p.$

将残量 ${\mathcal{R}}_h\left(x,t\right)$ 带入上式积分得到：

${\mathcal{M}}^k\frac{\text{d}{W}_h^k}{\text{d}t}+V\cdot \mathcal{S}{W}_h^k={\left[{\mathcal{l}}^k\left(x\right)\left(V{W}_h^k-{\left(VW\right)}^{*}\right)\right]}_{x_l^k}^{x_r^k}.$

其中：

$\begin{array}{l}{\mathcal{M}}_{ij}^k={\displaystyle {\int }_{x_l^k}^{x_r^k}{\phi }_i^k}\left(x\right){\phi }_j^k\left(x\right)\text{d}x=\frac{h^k}{2}{\displaystyle {\int }_{-1}^1{\phi }_i}\left(r\right){\phi }_j\left(r\right)\text{d}r\\ {\mathcal{S}}_{ij}^k={\displaystyle {\int }_{x_l^k}^{x_r^k}{\phi }_i^k}\left(x\right)\frac{\text{d}{\phi }_j^k\left(x\right)}{\text{d}x}\text{d}x={\displaystyle {\int }_{-1}^1{\phi }_i}\left(r\right)\frac{\text{d}{\phi }_j\left(r\right)}{\text{d}r}\text{d}r\end{array}$ .

接下来将局部解的表达式带入，即可得到一阶方程组的DG方法的强形式，如下：

$\frac{\text{d}{U}_h^k}{\text{d}t}+\frac{-c}{J^k}{\mathcal{D}}_r{P}_h^k=\frac{1}{J^k}{\mathcal{M}}^{-1}{\left[{\phi }^k\left(x\right)\left(\left(-c\right)H_h^k-\left(-c\right)H^{*}\right)\right]}_{x_l^k}^{x_r^k}=\frac{1}{J^k}{\mathcal{M}}^{-1}{{\displaystyle \oint }}_{x_l^k}^{x_r^k}\hat{n}\cdot \left(\left(-c\right)H_h^k-\left(-c\right)H^{*}\right){\phi }^k\left(x\right)\text{d}x.$

$\frac{\text{d}{P}_h^k}{\text{d}t}+\frac{-c}{J^k}{\mathcal{D}}_r{U}_h^k=\frac{1}{J^k}{\mathcal{M}}^{-1}{\left[{\phi }^k\left(x\right)\left(\left(-c\right)E_h^k-\left(-c\right)E^{*}\right)\right]}_{x_l^k}^{x_r^k}=\frac{1}{J^k}{\mathcal{M}}^{-1}{{\displaystyle \oint }}_{x_l^k}^{x_r^k}\hat{n}\cdot \left(\left(-c\right)E_h^k-\left(-c\right)E^{*}\right){\phi }^k\left(x\right)\text{d}x.$

其中： $J^k$ 为线性映射的Jacobi， ${\mathcal{D}}_{r,(i,j)}={\frac{\text{d}{\phi }_j}{\text{d}r}|}_{r_i}$。 ${\mathcal{D}}_r$ 和 ${\mathcal{S}}_{ij}$ 之间存在如下的计算关系：

${\left(\mathcal{M}{\mathcal{D}}_r\right)}_{(i,j)}={\displaystyle \underset{n=1}{\overset{N_p}{\sum }}{\mathcal{M}}_{in}}{\mathcal{D}}_{r,(n,j)}={{\displaystyle \underset{n=1}{\overset{N_p}{\sum }}{\displaystyle {\int }_{-1}^1{\phi }_i}}\left(r\right){\phi }_n\left(r\right)\frac{\text{d}{\phi }_j}{\text{d}r}|}_{r_n}\text{d}r={{\displaystyle {\int }_{-1}^1{\phi }_i}\left(r\right){\displaystyle \underset{n=1}{\overset{N_p}{\sum }}\frac{\text{d}{\phi }_j}{\text{d}r}}|}_{r_n}{\phi }_n\left(r\right)\text{d}r={\displaystyle {\int }_{-1}^1{\phi }_i}\left(r\right)\frac{\text{d}{\phi }_j\left(r\right)}{\text{d}r}\text{d}r={\mathcal{S}}_{ij}.$

至此即完成了DG方法的空间离散。

2.3. 时间离散

完成空间域离散后，可以得到半离散问题如下，其中 $W_h$ 为未知量。此时问题的求解为常微分方程 [9]。

$\frac{\text{d}{W}_h}{\text{d}t}={\mathcal{L}}_h\left(W_h,t\right).$

经典的基于数值通量的RK-DG方法采用Runge-Kutta进行时间离散，这是一种经典的单步法 [7]。一般的高阶Runge-Kutta时间离散格式表达式如下 [10]：

$\{\begin{array}{l}{u}^{(0)}=W^{(n)}\\ u^{(i)}={\displaystyle \underset{l=0}{\overset{i-1}{\sum }}\left[{\alpha }_{il}{u}^{(l)}+{\beta }_{il}\Delta tL\left(u^{(l)}\right)\right]i=1,\cdots ,r}\\ W^{(n+1)}=u^{(\; r\; )}\end{array}$

其中当 $r\le 4$ 时，上述时间离散格式的精度是 $r$ 阶。在这里选取变形的四阶显式RK方法(ERK)：其基本形式如下，相关系数可查阅文献 [5] 得到。

$\begin{array}{l}{p}^{(0)}=W^n\\ \{\begin{array}{l}{k}^{(i)}=a_i{k}^{(i-1)}+\Delta t{\mathcal{L}}_h\left(p^{(i-1)},t^n+c_i\Delta t\right),i\in [1,\cdots ,5]\\ p^{(i)}=p^{(i-1)}+b_i{k}^{(i)}\end{array}\\ W^{n+1}=p^{(\; 5\; )}\end{array}$

3. 计算方法的实现

在本节将介绍计算过程中选取的基函数、数值积分计算以及数值通量选择选择等。

3.1. 基函数

如果考虑基函数 ${\phi }_n\left(r\right)=r^{n-1}$，在计算过程中需要计算单元质量矩阵 $M_{ij}$，随着多项式次数的增加，得到的质量矩阵条件数很差，因为

${\mathcal{M}}_{ij}={\left({\phi }_i,{\phi }_j\right)}_I\text{=}{\displaystyle {\int }_{-1}^1{\phi }_i{\phi }_j}\text{d}x\text{=}\frac{1}{i+j-1}\left[1+{\left(-1\right)}^{i+j}\right].$

一个稳妥的办法是选择一组规范正交基来计算，可以通过对基函数进行Gram-Schmidt直交化来实现。其计算形式如下： $P_n\left(r\right)$ 为传统的n阶Lengendre多项式。由此基函数得到的质量矩阵要优于简单基函数。

${\phi }_n\left(r\right)={\tilde{P}}_{n-1}\left(r\right)=\frac{P_{n-1}\left(r\right)}{\sqrt{{\gamma }_{n-1}}}其中{\gamma }_n=\frac{2}{2n+1}.$

3.2. 数值积分

从DG方法的强形式可以看出，计算中需要大量的积分运算。如果使用精确积分，会导致巨大的计算量，因此可以考虑使用数值积分进行计算，同时计算格式要满足精度要求。根据已有经验来看，对于K次多项式积分运算，单元内部的积分应该使用精度至少为2K次的数值积分格式，对于单元边界上的积分，应使用精度至少为2K + 1次的格式。

目前使用较为广泛的是采用Legendre-Gauss-Lobatto (LGL)积分点进行积分 [11]。该积分点与等距节点积分和Gauss积分不同，等距积分形成的矩阵随着多项式次数的增加会出现列几乎线性相关的情况，不利于计算。Gauss积分通过张量积形式形成高维情形时，节点会集中在积分区域内部。而LGL积分点可以通过其他方法避免这些情况。查阅相关文献可以得到LGL积分点的计算方法，得到积分点与积分权重。

3.3. 数值通量

关于数值通量的选择，首选具有相容性、守恒性、单调性的数值通量。此种类型的数值通量最初来自于求解双曲守恒律方程的高分辨率格式的构造上，它不仅具有高精度，而且还能保持数值格式的稳定性和收敛性，具有很好的数学性质 [4] [5]。

首先定义单元边界处的平均值 $\left\{\left\{u\right\}\right\}$ 以及沿法向的跳跃 $\left[\left[u\right]\right]$ 如下，其中 $u^{-}$ 表示单元内部极限值， $u^{+}$ 表示单元外部极限值。

$\left\{\left\{u\right\}\right\}=\frac{u^{-}+u^{+}}{2},\left[\left[u\right]\right]={\hat{n}}^{-}{u}^{-}+{\hat{n}}^{+}{u}^{+}.$

目前常用的满足如上性质的数值通量函数有：

${\left(u\right)}^{*}=\left\{\left\{u\right\}\right\}+\frac{1-\alpha }{2}\left[\left[u\right]\right],0\le \alpha \le 1$， $\alpha =1$ 代表中心通量， $\alpha =0$ 代表迎风通量。

Lax-Friedrichs通量： ${\left(u\right)}^{*}=\left\{\left\{u\right\}\right\}+\frac{C}{2}\left(u^{+}-u^{-}\right)$，其中C与计算区域的情况有关。

目前，在声波方程的计算中比较常用的是Lax-Friedrichs通量，其计算公式相对简单，但不会影响计算的精度。

3.4. 时间步长

为了保持数值算法的稳定性，空间、时间步长需要满足一定的关系式，对该方法进行冯·诺伊曼稳定性分析，以便知道数值 ${\text{CFL}}_{L^{\text{2}}}$ 的值是什么。

由条件 $\left|c\right|\frac{\Delta t}{\Delta x}\le {\text{CFL}}_{L^{\text{2}}}$ 确保其 $L^2$ 稳定性，因为只有在这种条件下，舍入误差才不会被放大 [5]。

例如，对于使用k次多项式和k + 1阶达到k + 1次的RK方法(k + 1阶精度)的DG离散，实际中可以采取 $CF{L}_{L^2}=\frac{1}{2k+1}$。此外，对于 $k\ge 2$，数值 $\frac{1}{2k+1}$ 比 ${\text{CFL}}_{L^{\text{2}}}$ 数值估算值小5% [5]。故在本文中时间步长满足 $\Delta t\le \frac{{\text{CFL}}_{L^{\text{2}}}}{C}{\min }_{k,i}\Delta x_i^k$，其中空间步长为 $\Delta x_i^k=x_{i+1}^k-x_i^k$，C为传播中最大速度。

4. 数值实验

本节将通过两个例子来验证DG方法的计算效果。

4.1. 一阶单波方程

先来看一个一阶简单线性波动方程，该方程具有初值，形式如下：

$\{\begin{array}{l}\frac{\partial u}{\partial t}+2\text{π}\frac{\partial u}{\partial x}=0x\in \left[-2\text{π},2\text{π}\right],t>0\hfill \\ u\left(x,0\right)=\cos \left(-x\right)x\in \left[-2\text{π},2\text{π}\right]\hfill \\ u\left(-2\text{π},t\right)=\cos \left(2\text{π}t\right)t\ge 0\hfill \end{array}$ .

该问题的精确解为 $u\left(x,t\right)=\cos \left(-x+2\text{π}t\right)$。对于单个方程，其DG方法要比方程组更简单，这里便不再重复。图2为一阶波方程DG方法数值解与精确解的比较，其中单元基函数选取为4阶，区域剖分为4段(a)、8段(b)和16段(c)，数值通量选为中心通量，最终时间为 $t=1.5$，可以看出计算结果误差较小。

(a) (b) (c)

Figure 2. Comparison between numerical solution and exact solution of first-order wave equation

图2. 一阶波动方程数值解与精确解的比较

同时我们还计算了在不同单元个数K和多项式阶数N情况下的 $L^2$ 误差以及收敛阶，从表1中可以看出：随着单元剖分个数的增加和多项式阶数的提高，收敛阶始终趋于 $N+1$ 阶 [12]。

Table 1. Error and convergence order of numerical solution

表1. 数值解的误差与收敛阶

4.2. 二阶波动方程

对于二阶波动方程，带有初值条件如下：其中选取波速 $c=-1$，计算区间为 $\left[-10,10\right]$。

$\{\begin{array}{l}\frac{{\partial }^{2}U\left(x,t\right)}{\partial t^2}-c^2\frac{{\partial }^{2}U\left(x,t\right)}{\partial x^2}=0x\in \left[-10,10\right],t>0\\ U\left(x,0\right)=\sin \left(x\right)x\in \left[-10,10\right]\\ \frac{\partial U\left(x,0\right)}{\partial t}=c\sin \left(x\right)x\in \left[-10,10\right]\end{array}$ .

在实际计算中，选取3阶基函数进行插值，数值通量选为Lax-Friedrichs通量，最终计算时间为 $t=1$。结果如图3：其中a为计算空间10等分，b为计算空间20等分，c为计算空间30等分。

(a) (b) (c)

Figure 3. Comparison between numerical solution and exact solution of second-order wave equation

图3. 二阶波动方程数值解与精确解的比较

对其进行误差分析和收敛阶计算，即使用N阶多项式模拟时，随着剖分单元数的增加，收敛阶逐渐趋向于 $N+1$ 阶。

5. 结论

在本文中，通过对二阶无源波动方程空间域进行间断有限元离散，在时间上采用经典的Rutta-Kutta方法进行离散，得到了较高的收敛阶，这和DG方法的理论分析是一致的。在计算过程中，选取经Gram-Schmidt直交化的Lengendre多项式作为基函数，避免了随着多项式阶数的提高质量矩阵条件数变坏的情况。数值积分选取LGL积分点和权重，避免了等距积分形成行列式线性相关的情况，时间离散上采用采用经典的Rutta-Kutta单步法四阶变形形式，减少了计算中的存储用量。最后通过两个实例分析验证了间断有限元方法可以获得更高的收敛阶，并且提高多项式的阶数要比传统有限元更方便，在未来会获得更广泛的应用。

基金项目

中国石油大学(北京)科研基金项目资助(编号：2462020XKJS02)。