一种无数值积分的间断Galerkin有限元方法

马小乐, 曹 伟

(天津大学 机械工程学院力学系, 天津 300072)

0 引 言

间断Galerkin有限元方法(DG-FEM)最早是由Reed和Hill在1973年提出并应用于求解中子输运方程[1]，且由Lesaint和Raviart在1974年首次对该方法进行了收敛阶分析，并应用于水平对流方程的数值计算中[2]。此后，Wellford和Oden在1975-1976年间，首先使用DG方法进行了波在弹性介质中的传播研究[3]，Delfour和Trochu也在1978年将此方法运用到了最佳控制问题[4]，但均局限求解线性方程。Chavent等最早将间断Galerkin有限元方法推广到非线性守恒律问题[5]，并通过引入梯度限制器使得计算格式稳定[6]，但在时间和空间方向上均只能实现低阶收敛。20世纪80年代末到90年代初，Cockburn和Shu通过引进稳定的Runge-Kutta方法[7]并配合改进的梯度限制器[8]，克服了这一局限，保证了在光滑区域上获得正常的数值精度以及清晰的间断解[9]，继而使用广义梯度限制器和具有高阶精度的RK方法使得数值结果重新获得了高阶精度[10]并将其推广到一维方程组[11]、高维标量守恒律[12]及高维方程组[13]，即著名的RKDG方法。

间断Galerkin有限元方法具有良好的守恒性、收敛性等数学特性，并且结合了有限体积法与连续有限元法的诸多优点。首先，该方法只需要简单地增加单元自由度即可实现高阶精度，且每个单元都是相互独立的，对一般非规则网格具有更强的适应性，利于处理复杂的区域边界和具有复杂边界条件的问题；另外，该方法的质量矩阵是单元的而非整体的，对该矩阵求逆相对容易，继而可以方便的构造显示半离散格式；同时，数值通量具体形式的定义具备很强的灵活性，不同的定义方法可以反映具体问题的物理特性，从而保证波占优问题的稳定性。

尽管间断Galerkin有限元方法拥有众多的优点，但是该方法在每个单元需要求解的变量数都有所增加，而且随着精度的提高，变量数是非线性增长的，在复杂外形的大型数值模拟过程中，所需的计算量和存储量已无法满足实际工程问题的数值模拟需求。构造更加高效的隐式算法[14]、引入并行计算[15]等成为缓解上述问题的有效办法。Aktins和Shu则从另一个角度出发，针对使用间断Galerkin有限元方法求解过程中数值积分带来的计算量和存储量大的不利影响，提出了无数值积分的DG方法[16-17]以解决这个问题。

不同于Aktins和Shu以简单的单项式作为基函数的无数值积分DG方法，本文在单元基函数为Lagrange插值多项式的DG格式的基础上，又引入了基于Jacobi正交多项式的单元近似函数表达式，并通过应用该多项式的正交与求导等性质，无需数值积分的运算过程，方便、直接地构造出了一维和二维守恒律间断Galerkin有限元方法的显式半离散格式。基于该格式相关思想，将有利于构造更为高效的高阶间断Galerkin有限元计算方法。

1 数值方法

1.1 一维守恒律

考虑一维守恒律

在每一个单元上应用Galerkin加权余量法的基本思想，使单元方程余量正交于该单元内的Np个单元权(基)函数，即

对方程(4)进行分部积分，并引入数值通量[18]f*代替单元边界处的通量项，之后再做一次分部积分，即可构造出间断Galerkin有限元方法下一维守恒律的积分表达式

将单元多项式近似解函数(2)和单元多项式近似通量函数(3)带入方程(5)中，整理可得相应矩阵形式的积分表达式

以积分表达式(6)为数值求解的出发点，一般需要通过数值积分的方法确定每一个单元上的质量矩阵和刚度矩阵，为了避免这一运算过程，首先引入参考变量r∈I=-1,1及线性变换

继而单元多项式近似解函数(2)可以写为参考变量的形式

在参考区间I上的Np个插值节点rj(j=1,2,…,Np)选取为多项式

在选定了单元插值节点以后，引入单元多项式近似解函数的另一种表达形式

由于(8)式和(13)式逼近的是同一函数的同阶多项式，因此一定存在恒等关系

将LGL点ri(i=1,2,…,Np)依次带入关系式(14)中并定义广义Vandermonde矩阵

则有

基于以上关系表达式，可以将任意单元Dk上的质量矩阵变换至参考单元上的质量矩阵Mij

=JkMij(17)

Jk为线性坐标变换下的雅克比行列式。

将式(16)代入Mij的表达式中，并应用正交关系(12)式(取α=0,β=0)，则有

Mij

即参考单元质量矩阵的表达式为

M=(VVT)-1(19)

对于单元刚度矩阵，由线性坐标变换可直接转化为参考区间上的刚度矩阵，即

=Sij(20)

将(16)式两端同时对r求导，并定义

则可得微分矩阵Dr的表达式为

Dr=VrV-1(22)

将参考单元上的质量矩阵与微分矩阵相乘

即参考单元刚度矩阵的表达式为

S=MDr(24)

至此，将式(17)、(20)代入积分表达式(6)中，方程两端对M求逆并应用矩阵关系式(24)，整理即可得显式半离散格式

观察该显示半离散格式，M和Dr可通过(19)和(22)式获得，且矩阵V和Vr可由LGL点以及Jacobi正交多项式及其关系式(10)和(11)唯一确定，该过程简单有效，且无需引入数值积分的运算过程。

1.2 二维守恒律

考虑二维守恒律

及相应的初始条件和边界条件，其中，x=(x,y),f=(f1,f2)。

使用K个直边三角形单元Dk对Ω进行相容的几何剖分，并在每个单元上定义N阶单元多项式近似解函数和单元多项式近似通量函数

引入参考单元I={r=(r,s)|(r,s)≥-1,r+s≤0}及坐标变换

其中，v1、v2、v3表示单元Dk上三个顶点的物理坐标，则可得积分表达式的矩阵形式

式中，Jk为坐标变换Jacobi行列式，M,Sr和Ss为参考单元上的质量矩阵及刚度矩阵，且分别定义为

其中，0≤λ1,λ2,λ3≤1, (i,j)≥0,i+j≤N，继而可以得到这些节点的物理坐标xe(λ1,λ2,λ3)。

引入增量函数

在确定二维插值节点之后，同样需要引入另一种单元多项式近似函数表达形式

且φn(r)的具体表达式为

仿照一维守恒律构造相关矩阵的思想，由式(37)及插值节点可以定义矩阵

Vij=φj(ri) (38)

继而可以得到二维格式的质量矩阵、微分矩阵和刚度矩阵的表达式

M=(VVT)-1(40)

Dr=VrV-1，Ds=VsV-1(41)

Sr=MDr，Ss=MDs(42)

与一维守恒律稍有不同，需要计算面积分项

由于在直边三角形单元上，当单元近视解函数为N阶多项式时，每条边上将刚好分布有N+1个节点，以其中一条边为例可计算积分如下

其中1≤i≤Np,1≤j≤N+1，且

为定义在三角单元边上的质量矩阵，该矩阵的元素可以由边上的节点按照一维问题处理得到。综上，即可完整的构造出二维守恒律DG方法的显式半离散格式。

2 算例与分析

2.1 一维线性波动方程

考虑具体的一维线性波动方程

且受制于如下初始条件和边界条件

该方程的精确解为u(x,t)=sin(x-2πt)，选取Lax-Friedrichs数值通量，并使用Runge-Kutta方法进行时间上的迭代求解。表1给出了T=10时刻，N和K取不同值时L1范数定义下的整体误差及相应条件下的数值精度。

分析表中的结果可知，首先，无论是增加单元个数K还是增加单元多项式近似函数阶数N，整体误差均是逐渐减小的，即该格式是明显收敛的；此外，当单元多项式近似函数阶数N保持不变时，随着单元个数K的增加，由整体误差所确定的数值精度均逐渐逼近于间断Galerkin有限元方法的理论精度N+1阶。

表1 求解单波方程得到的L1误差及数值精度Table 1 L1 error and numerical precision of wave equation

2.2 一维可压缩Euler方程组

考虑经典激波管(SOD)问题，其控制方程为一维守恒型可压缩Euler方程组

其中，u、ρ、p分别为流向速度，密度和压力，E是单位体积动能和内能的和，即E=ρ(e+u2/2)，此外，由完全气体假设下的状态方程可得p=(γ-1)(E-ρu2/2)。

初始条件给定为

对该问题应用与一维线性波动方程相同的构造方法进行求解，为消除非物理振荡，需要在RK方法的每一步使用经典MUSCL限制器[23]。

图1和图2分别给出了使用该方法在K=500,N=2和N=3时，T=0.2时刻求解激波管问题所获得的数值解与精确解的比较结果。由图可以看出，该方法可以有效的捕捉到激波的准确位置，在光滑区域获得精度较高的解，同时，在限制器的作用下，可以完全消除非物理振荡。

(a)

(b)

(c)

(d)

此外，对比图1和图2所示结果可以发现，两种情况下所得的数值解差别不大，这是由于在限制器的作用下，不可避免的影响了该方法运用高阶基函数时的优势，针对这一问题，可以通过增加单元个数或者使用高阶基函数配合更先进的限制器使计算结果得到改进。

(a)

(b)

(c)

(d)

2.3 二维可压缩Euler方程组

考虑控制方程为二维守恒型Euler方程组

其中，u,v,p,ρ分别为流向速度，法向速度，压力和密度，E是单位体积动能和内能的和，且有E=ρ[e+(u2+v2)/2]，以及完全气体假设下的状态方程p=(γ-1)[E-ρ(u2+v2)/2]。对此方程可以根据二维守恒律的处理方法构造相应的显式半离散格式。

考虑具有精确解

表2显示了由初始网格逐步均匀加密后得到的L1范数定义下密度ρ的整体误差以及相应条件下的数值精度，h代表初始网格中各单元的大小。

表2 密度ρ的L1误差及数值精度Table 2 L1 error and numerical precision of density

从表中的结果可以看出，同样地，二维模式下该格式仍是明显收敛的，虽然由于非结构网格导致数值精度存在一定的波动，但总体上仍趋向于间断Galerkin有限元方法的最优收敛阶N+1阶。

考虑同样由二维可压缩Euler方程组控制的前台阶问题，初始条件给定为马赫数为3的超声速均匀流场

流入边界条件由初始值给定，墙边界条件取为如下反射边界条件

因假设流出边界处仍为超声速流场，所以不需要强加边界条件。求解过程中，在RK方法的每一时间步均使用了适用于二维可压缩Euler方程组的梯度的限制器[24]。

图3和图4分别显示了选取K=25330，N=2和N=3时，在T=4时刻整个求解区域上密度ρ，压力p以及马赫数Ma的等值线图。将该方法的计算结果与文献[25]中的结果进行比较，整体求解区域上均吻合的很好，数值解达到了理想的精度，且激波的位置被准确的捕捉到。同时，与一维情况相同，限制器的作用对选取高阶基函数时的计算结果造成了一定程度的影响，同样的，可以通过增加单元个数或使用更先进的限制器加以改进。

(a) ρ

(b) p

(c) Ma

(a) ρ

(b) p

(c) Ma

3 结 论

本文以单元基函数为Lagrange插值多项式，单元插值节点为LGL点(二维情形下插值节点为LGL点的扩展)的DG格式作为求解出发点，同时引入了基于Jacobi正交多项式的单元近似函数表达式，通过建立两种不同基函数的联系，构造了一维和二维守恒律无数值积分的间断Galerkin有限元方法显式半离散格式。应用该方法对具有连续解或间断解的不同算例进行数值模拟，验证了其所得的计算结果可以很好的达到间断Galerkin有限元方法的理论高阶精度，且该方法配合适当的限制器，可以有效的处理含有间断的问题。由于不再需要通过数值积分来计算各单元的体积分与面积分，在保证了高阶精度的同时，理论上可以减少数值计算所需内存量及计算量。但是，在处理含间断问题时，为了完全抑制数值振荡，使用限制器是必须的，继而不可避免的会对激波造成污染，同时会一定程度地破坏解在光滑区域上的高阶精度，也即在一定程度上限制了该计算方法运用高阶基函数的优势。进一步的研究，可以应用该思想构造高阶以及三维问题的无数值积分DG求解格式，并具体地对比验证其在内存使用量及计算量方面的优势，继而考虑配合其他方法及更为先进的限制器构造更为高效的且可以达到高阶精度的间断Galerkin有限元方法计算格式。