馬小樂(lè), 曹 偉
(天津大學(xué) 機(jī)械工程學(xué)院力學(xué)系, 天津 300072)
間斷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ì)算方法。
考慮一維守恒律
在每一個(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ò)程。
考慮二維守恒律
及相應(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方法的顯式半離散格式。
考慮具體的一維線(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
考慮經(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)
考慮控制方程為二維守恒型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
本文以單元基函數(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ì)算格式。