馬小樂, 曹 偉

(天津大學 機械工程學院力學系, 天津 300072)

0 引 言

間斷Galerkin有限元方法(DG-FEM)最早是由Reed和Hill在1973年提出并應用于求解中子輸運方程[1]，且由Lesaint和Raviart在1974年首次對該方法進行了收斂階分析，并應用于水平對流方程的數(shù)值計算中[2]。此后，Wellford和Oden在1975-1976年間，首先使用DG方法進行了波在彈性介質中的傳播研究[3]，Delfour和Trochu也在1978年將此方法運用到了最佳控制問題[4]，但均局限求解線性方程。Chavent等最早將間斷Galerkin有限元方法推廣到非線性守恒律問題[5]，并通過引入梯度限制器使得計算格式穩(wěn)定[6]，但在時間和空間方向上均只能實現(xiàn)低階收斂。20世紀80年代末到90年代初，Cockburn和Shu通過引進穩(wěn)定的Runge-Kutta方法[7]并配合改進的梯度限制器[8]，克服了這一局限，保證了在光滑區(qū)域上獲得正常的數(shù)值精度以及清晰的間斷解[9]，繼而使用廣義梯度限制器和具有高階精度的RK方法使得數(shù)值結果重新獲得了高階精度[10]并將其推廣到一維方程組[11]、高維標量守恒律[12]及高維方程組[13]，即著名的RKDG方法。

間斷Galerkin有限元方法具有良好的守恒性、收斂性等數(shù)學特性，并且結合了有限體積法與連續(xù)有限元法的諸多優(yōu)點。首先，該方法只需要簡單地增加單元自由度即可實現(xiàn)高階精度，且每個單元都是相互獨立的，對一般非規(guī)則網(wǎng)格具有更強的適應性，利于處理復雜的區(qū)域邊界和具有復雜邊界條件的問題；另外，該方法的質量矩陣是單元的而非整體的，對該矩陣求逆相對容易，繼而可以方便的構造顯示半離散格式；同時，數(shù)值通量具體形式的定義具備很強的靈活性，不同的定義方法可以反映具體問題的物理特性，從而保證波占優(yōu)問題的穩(wěn)定性。

盡管間斷Galerkin有限元方法擁有眾多的優(yōu)點，但是該方法在每個單元需要求解的變量數(shù)都有所增加，而且隨著精度的提高，變量數(shù)是非線性增長的，在復雜外形的大型數(shù)值模擬過程中，所需的計算量和存儲量已無法滿足實際工程問題的數(shù)值模擬需求。構造更加高效的隱式算法[14]、引入并行計算[15]等成為緩解上述問題的有效辦法。Aktins和Shu則從另一個角度出發(fā)，針對使用間斷Galerkin有限元方法求解過程中數(shù)值積分帶來的計算量和存儲量大的不利影響，提出了無數(shù)值積分的DG方法[16-17]以解決這個問題。

不同于Aktins和Shu以簡單的單項式作為基函數(shù)的無數(shù)值積分DG方法，本文在單元基函數(shù)為Lagrange插值多項式的DG格式的基礎上，又引入了基于Jacobi正交多項式的單元近似函數(shù)表達式，并通過應用該多項式的正交與求導等性質，無需數(shù)值積分的運算過程，方便、直接地構造出了一維和二維守恒律間斷Galerkin有限元方法的顯式半離散格式?；谠摳袷较嚓P思想，將有利于構造更為高效的高階間斷Galerkin有限元計算方法。

1 數(shù)值方法

1.1 一維守恒律

考慮一維守恒律

在每一個單元上應用Galerkin加權余量法的基本思想，使單元方程余量正交于該單元內的Np個單元權(基)函數(shù)，即

對方程(4)進行分部積分，并引入數(shù)值通量[18]f*代替單元邊界處的通量項，之后再做一次分部積分，即可構造出間斷Galerkin有限元方法下一維守恒律的積分表達式

將單元多項式近似解函數(shù)(2)和單元多項式近似通量函數(shù)(3)帶入方程(5)中，整理可得相應矩陣形式的積分表達式

以積分表達式(6)為數(shù)值求解的出發(fā)點，一般需要通過數(shù)值積分的方法確定每一個單元上的質量矩陣和剛度矩陣，為了避免這一運算過程，首先引入?yún)⒖甲兞縭∈I=-1,1及線性變換

繼而單元多項式近似解函數(shù)(2)可以寫為參考變量的形式

在參考區(qū)間I上的Np個插值節(jié)點rj(j=1,2,…,Np)選取為多項式

在選定了單元插值節(jié)點以后，引入單元多項式近似解函數(shù)的另一種表達形式

由于(8)式和(13)式逼近的是同一函數(shù)的同階多項式，因此一定存在恒等關系

將LGL點ri(i=1,2,…,Np)依次帶入關系式(14)中并定義廣義Vandermonde矩陣

則有

基于以上關系表達式，可以將任意單元Dk上的質量矩陣變換至參考單元上的質量矩陣Mij

=JkMij(17)

Jk為線性坐標變換下的雅克比行列式。

將式(16)代入Mij的表達式中，并應用正交關系(12)式(取α=0,β=0)，則有

Mij

即參考單元質量矩陣的表達式為

M=(VVT)-1(19)

對于單元剛度矩陣，由線性坐標變換可直接轉化為參考區(qū)間上的剛度矩陣，即

=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)唯一確定，該過程簡單有效，且無需引入數(shù)值積分的運算過程。

1.2 二維守恒律

考慮二維守恒律

及相應的初始條件和邊界條件，其中，x=(x,y),f=(f1,f2)。

使用K個直邊三角形單元Dk對Ω進行相容的幾何剖分，并在每個單元上定義N階單元多項式近似解函數(shù)和單元多項式近似通量函數(shù)

引入?yún)⒖紗卧狪={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，繼而可以得到這些節(jié)點的物理坐標xe(λ1,λ2,λ3)。

引入增量函數(shù)

在確定二維插值節(jié)點之后，同樣需要引入另一種單元多項式近似函數(shù)表達形式

且φn(r)的具體表達式為

仿照一維守恒律構造相關矩陣的思想，由式(37)及插值節(jié)點可以定義矩陣

Vij=φj(ri) (38)

繼而可以得到二維格式的質量矩陣、微分矩陣和剛度矩陣的表達式

M=(VVT)-1(40)

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

Sr=MDr，Ss=MDs(42)

與一維守恒律稍有不同，需要計算面積分項

由于在直邊三角形單元上，當單元近視解函數(shù)為N階多項式時，每條邊上將剛好分布有N+1個節(jié)點，以其中一條邊為例可計算積分如下

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

為定義在三角單元邊上的質量矩陣，該矩陣的元素可以由邊上的節(jié)點按照一維問題處理得到。綜上，即可完整的構造出二維守恒律DG方法的顯式半離散格式。

2 算例與分析

2.1 一維線性波動方程

考慮具體的一維線性波動方程

且受制于如下初始條件和邊界條件

該方程的精確解為u(x,t)=sin(x-2πt)，選取Lax-Friedrichs數(shù)值通量，并使用Runge-Kutta方法進行時間上的迭代求解。表1給出了T=10時刻，N和K取不同值時L1范數(shù)定義下的整體誤差及相應條件下的數(shù)值精度。

分析表中的結果可知，首先，無論是增加單元個數(shù)K還是增加單元多項式近似函數(shù)階數(shù)N，整體誤差均是逐漸減小的，即該格式是明顯收斂的；此外，當單元多項式近似函數(shù)階數(shù)N保持不變時，隨著單元個數(shù)K的增加，由整體誤差所確定的數(shù)值精度均逐漸逼近于間斷Galerkin有限元方法的理論精度N+1階。

表1 求解單波方程得到的L1誤差及數(shù)值精度Table 1 L1 error and numerical precision of wave equation

2.2 一維可壓縮Euler方程組

考慮經(jīng)典激波管(SOD)問題，其控制方程為一維守恒型可壓縮Euler方程組

其中，u、ρ、p分別為流向速度，密度和壓力，E是單位體積動能和內能的和，即E=ρ(e+u2/2)，此外，由完全氣體假設下的狀態(tài)方程可得p=(γ-1)(E-ρu2/2)。

初始條件給定為

對該問題應用與一維線性波動方程相同的構造方法進行求解，為消除非物理振蕩，需要在RK方法的每一步使用經(jīng)典MUSCL限制器[23]。

圖1和圖2分別給出了使用該方法在K=500,N=2和N=3時，T=0.2時刻求解激波管問題所獲得的數(shù)值解與精確解的比較結果。由圖可以看出，該方法可以有效的捕捉到激波的準確位置，在光滑區(qū)域獲得精度較高的解，同時，在限制器的作用下，可以完全消除非物理振蕩。

(a)

(b)

(c)

(d)

此外，對比圖1和圖2所示結果可以發(fā)現(xiàn)，兩種情況下所得的數(shù)值解差別不大，這是由于在限制器的作用下，不可避免的影響了該方法運用高階基函數(shù)時的優(yōu)勢，針對這一問題，可以通過增加單元個數(shù)或者使用高階基函數(shù)配合更先進的限制器使計算結果得到改進。

(a)

(b)

(c)

(d)

2.3 二維可壓縮Euler方程組

考慮控制方程為二維守恒型Euler方程組

其中，u,v,p,ρ分別為流向速度，法向速度，壓力和密度，E是單位體積動能和內能的和，且有E=ρ[e+(u2+v2)/2]，以及完全氣體假設下的狀態(tài)方程p=(γ-1)[E-ρ(u2+v2)/2]。對此方程可以根據(jù)二維守恒律的處理方法構造相應的顯式半離散格式。

考慮具有精確解

表2顯示了由初始網(wǎng)格逐步均勻加密后得到的L1范數(shù)定義下密度ρ的整體誤差以及相應條件下的數(shù)值精度，h代表初始網(wǎng)格中各單元的大小。

表2 密度ρ的L1誤差及數(shù)值精度Table 2 L1 error and numerical precision of density

從表中的結果可以看出，同樣地，二維模式下該格式仍是明顯收斂的，雖然由于非結構網(wǎng)格導致數(shù)值精度存在一定的波動，但總體上仍趨向于間斷Galerkin有限元方法的最優(yōu)收斂階N+1階。

考慮同樣由二維可壓縮Euler方程組控制的前臺階問題，初始條件給定為馬赫數(shù)為3的超聲速均勻流場

流入邊界條件由初始值給定，墻邊界條件取為如下反射邊界條件

因假設流出邊界處仍為超聲速流場，所以不需要強加邊界條件。求解過程中，在RK方法的每一時間步均使用了適用于二維可壓縮Euler方程組的梯度的限制器[24]。

圖3和圖4分別顯示了選取K=25330，N=2和N=3時，在T=4時刻整個求解區(qū)域上密度ρ，壓力p以及馬赫數(shù)Ma的等值線圖。將該方法的計算結果與文獻[25]中的結果進行比較，整體求解區(qū)域上均吻合的很好，數(shù)值解達到了理想的精度，且激波的位置被準確的捕捉到。同時，與一維情況相同，限制器的作用對選取高階基函數(shù)時的計算結果造成了一定程度的影響，同樣的，可以通過增加單元個數(shù)或使用更先進的限制器加以改進。

(a) ρ

(b) p

(c) Ma

(a) ρ

(b) p

(c) Ma

3 結 論

本文以單元基函數(shù)為Lagrange插值多項式，單元插值節(jié)點為LGL點(二維情形下插值節(jié)點為LGL點的擴展)的DG格式作為求解出發(fā)點，同時引入了基于Jacobi正交多項式的單元近似函數(shù)表達式，通過建立兩種不同基函數(shù)的聯(lián)系，構造了一維和二維守恒律無數(shù)值積分的間斷Galerkin有限元方法顯式半離散格式。應用該方法對具有連續(xù)解或間斷解的不同算例進行數(shù)值模擬，驗證了其所得的計算結果可以很好的達到間斷Galerkin有限元方法的理論高階精度，且該方法配合適當?shù)南拗破鳎梢杂行У奶幚砗虚g斷的問題。由于不再需要通過數(shù)值積分來計算各單元的體積分與面積分，在保證了高階精度的同時，理論上可以減少數(shù)值計算所需內存量及計算量。但是，在處理含間斷問題時，為了完全抑制數(shù)值振蕩，使用限制器是必須的，繼而不可避免的會對激波造成污染，同時會一定程度地破壞解在光滑區(qū)域上的高階精度，也即在一定程度上限制了該計算方法運用高階基函數(shù)的優(yōu)勢。進一步的研究，可以應用該思想構造高階以及三維問題的無數(shù)值積分DG求解格式，并具體地對比驗證其在內存使用量及計算量方面的優(yōu)勢，繼而考慮配合其他方法及更為先進的限制器構造更為高效的且可以達到高階精度的間斷Galerkin有限元方法計算格式。