馬小樂, 曹 偉

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

0 引 言

間斷Galerkin有限元方法(DG-FEM)最早是由Reed和Hill在1973年提出并應(yīng)用于求解中子輸運方程[1]，且由Lesaint和Raviart在1974年首次對該方法進行了收斂階分析，并應(yīng)用于水平對流方程的數(shù)值計算中[2]。此后，Wellford和Oden在1975-1976年間，首先使用DG方法進行了波在彈性介質(zhì)中的傳播研究[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ù)值結(jié)果重新獲得了高階精度[10]并將其推廣到一維方程組[11]、高維標量守恒律[12]及高維方程組[13]，即著名的RKDG方法。

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

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

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

1 數(shù)值方法

1.1 一維守恒律

考慮一維守恒律

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

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

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

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

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

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

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

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

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

則有

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

=JkMij(17)

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

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

Mij

即參考單元質(zhì)量矩陣的表達式為

M=(VVT)-1(19)

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

=Sij(20)

將(16)式兩端同時對r求導(dǎo)，并定義

則可得微分矩陣Dr的表達式為

Dr=VrV-1(22)

將參考單元上的質(zhì)量矩陣與微分矩陣相乘

即參考單元剛度矩陣的表達式為

S=MDr(24)

至此，將式(17)、(20)代入積分表達式(6)中，方程兩端對M求逆并應(yīng)用矩陣關(guān)系式(24)，整理即可得顯式半離散格式

觀察該顯示半離散格式，M和Dr可通過(19)和(22)式獲得，且矩陣V和Vr可由LGL點以及Jacobi正交多項式及其關(guān)系式(10)和(11)唯一確定，該過程簡單有效，且無需引入數(shù)值積分的運算過程。

1.2 二維守恒律

考慮二維守恒律

及相應(yīng)的初始條件和邊界條件，其中，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為參考單元上的質(zhì)量矩陣及剛度矩陣，且分別定義為

其中，0≤λ1,λ2,λ3≤1, (i,j)≥0,i+j≤N，繼而可以得到這些節(jié)點的物理坐標xe(λ1,λ2,λ3)。

引入增量函數(shù)

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

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

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

Vij=φj(ri) (38)

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

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，且

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

2 算例與分析

2.1 一維線性波動方程

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

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

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

分析表中的結(jié)果可知，首先，無論是增加單元個數(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是單位體積動能和內(nèi)能的和，即E=ρ(e+u2/2)，此外，由完全氣體假設(shè)下的狀態(tài)方程可得p=(γ-1)(E-ρu2/2)。

初始條件給定為

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

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

(a)

(b)

(c)

(d)

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

(a)

(b)

(c)

(d)

2.3 二維可壓縮Euler方程組

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

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

考慮具有精確解

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

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

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

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

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

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

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

(a) ρ

(b) p

(c) Ma

(a) ρ

(b) p

(c) Ma

3 結(jié) 論

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