一種無(wú)數(shù)值積分的間斷Galerkin有限元方法

馬小樂(lè), 曹 偉

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

0 引 言

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

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

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

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

1 數(shù)值方法

1.1 一維守恒律

考慮一維守恒律

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

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

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

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

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

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

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

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

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

則有

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

=JkMij(17)

Jk為線(xiàn)性坐標(biāo)變換下的雅克比行列式。

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

Mij

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

M=(VVT)-1(19)

對(duì)于單元?jiǎng)偠染仃?，由線(xiàn)性坐標(biāo)變換可直接轉(zhuǎn)化為參考區(qū)間上的剛度矩陣，即

=Sij(20)

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

則可得微分矩陣Dr的表達(dá)式為

Dr=VrV-1(22)

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

即參考單元?jiǎng)偠染仃嚨谋磉_(dá)式為

S=MDr(24)

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

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

1.2 二維守恒律

考慮二維守恒律

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

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

引入?yún)⒖紗卧狪={r=(r,s)|(r,s)≥-1,r+s≤0}及坐標(biāo)變換

其中，v1、v2、v3表示單元Dk上三個(gè)頂點(diǎn)的物理坐標(biāo)，則可得積分表達(dá)式的矩陣形式

式中，Jk為坐標(biāo)變換Jacobi行列式，M,Sr和Ss為參考單元上的質(zhì)量矩陣及剛度矩陣，且分別定義為

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

引入增量函數(shù)

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

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

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

Vij=φj(ri) (38)

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

M=(VVT)-1(40)

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

Sr=MDr，Ss=MDs(42)

與一維守恒律稍有不同，需要計(jì)算面積分項(xiàng)

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

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

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

2 算例與分析

2.1 一維線(xiàn)性波動(dòng)方程

考慮具體的一維線(xiàn)性波動(dòng)方程

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

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

分析表中的結(jié)果可知，首先，無(wú)論是增加單元個(gè)數(shù)K還是增加單元多項(xiàng)式近似函數(shù)階數(shù)N，整體誤差均是逐漸減小的，即該格式是明顯收斂的；此外，當(dāng)單元多項(xiàng)式近似函數(shù)階數(shù)N保持不變時(shí)，隨著單元個(gè)數(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)問(wèn)題，其控制方程為一維守恒型可壓縮Euler方程組

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

初始條件給定為

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

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

(a)

(b)

(c)

(d)

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

(a)

(b)

(c)

(d)

2.3 二維可壓縮Euler方程組

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

其中，u,v,p,ρ分別為流向速度，法向速度，壓力和密度，E是單位體積動(dòng)能和內(nèi)能的和，且有E=ρ[e+(u2+v2)/2]，以及完全氣體假設(shè)下的狀態(tài)方程p=(γ-1)[E-ρ(u2+v2)/2]。對(duì)此方程可以根據(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ù)值精度存在一定的波動(dòng)，但總體上仍趨向于間斷Galerkin有限元方法的最優(yōu)收斂階N+1階。

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

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

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

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

(a) ρ

(b) p

(c) Ma

(a) ρ

(b) p

(c) Ma

3 結(jié) 論

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