基于LU-SGS迭代的DGM隱式方法研究

郝海兵, 張強(qiáng), 楊永, 梁益華

間斷Galerkin方法[1](DGM)由于具備易于實(shí)現(xiàn)高精度以及求解間斷問(wèn)題等優(yōu)點(diǎn)而受到廣泛關(guān)注。該方法在任意單元內(nèi)部通過(guò)提高逼近多項(xiàng)式的階次實(shí)現(xiàn)高精度;空間離散時(shí)允許單元邊界上存在間斷,具有靈活處理復(fù)雜區(qū)域邊值問(wèn)題的能力;在單元邊界處可以采用有限體積法中的思路,基于Riemann問(wèn)題來(lái)構(gòu)造數(shù)值通量函數(shù),實(shí)現(xiàn)逆風(fēng)格式,從而更有利于求解間斷問(wèn)題。目前,國(guó)外對(duì)該方法已經(jīng)展開(kāi)了近十多年的研究,并取得了較大的進(jìn)展[2],而國(guó)內(nèi)相關(guān)研究起步較晚,在近幾年才逐漸受到重視。國(guó)內(nèi)主要進(jìn)展包括邱建賢等[3]將WENO格式作為限制器引入到DGM中;張來(lái)平等[4]提出了靜動(dòng)態(tài)混合重構(gòu)的DG/FV混合格式等。當(dāng)前DGM主要應(yīng)用于空間離散,而時(shí)間方向推進(jìn)多沿用了Shu等[2]構(gòu)造的TVD-RKDG方法。因此,當(dāng)求解更高精度或大規(guī)模計(jì)算問(wèn)題時(shí),穩(wěn)定性條件變得越來(lái)越苛刻,從而導(dǎo)致當(dāng)?shù)貢r(shí)間步長(zhǎng)取值越來(lái)越小,收斂率會(huì)急劇下降,計(jì)算非常耗時(shí)。針對(duì)該問(wèn)題,目前主要采用p型多重網(wǎng)格方法[2]和隱式時(shí)間離散2種加速收斂技術(shù)來(lái)提高計(jì)算效率。本文主要對(duì)隱式格式離散進(jìn)行研究,p型多重網(wǎng)格方法見(jiàn)參考文獻(xiàn)[2,5]。目前,國(guó)外很多學(xué)者對(duì)間斷Galerkin有限元隱式算法進(jìn)行了多種試探性的研究,不過(guò)目前還沒(méi)有形成成熟穩(wěn)定的算法。主要進(jìn)展包括如下:Bassi和Rebay首次將GMRES隱式算法應(yīng)用于DG求解二維可壓縮N-S方程[6];Hartmann等人使用GMRES-Newton算法求解二維Euler和N-S方程[7]等。

本文基于非結(jié)構(gòu)網(wǎng)格,針對(duì)定常Euler方程組,發(fā)展相應(yīng)的DGM隱式算法?？紤]到LU-SGS迭代法不需要存儲(chǔ)和處理大型矩陣,并且非常適合于非結(jié)構(gòu)網(wǎng)格的無(wú)序性,具有顯著的計(jì)算效率,并且在有限體積法中得到了廣泛應(yīng)用。本文借鑒其思想,將該方法推廣到高精度DGM隱式格式中,構(gòu)造適合于DGM隱式求解的LU-SGS迭代格式。通過(guò)對(duì)NACA0012翼型和ONERA M6機(jī)翼跨聲速無(wú)粘流動(dòng)進(jìn)行數(shù)值模擬驗(yàn)證其加速收斂性能。

1 控制方程

考慮非定常Euler方程在直角坐標(biāo)系的守恒形式:

(1)

式中:

ρ、e、p分別表示氣體的密度、單位體積總內(nèi)能及壓力,ui表示xi方向的速度。對(duì)于氣體動(dòng)力學(xué)方程。

γ=1.4是比熱比。

2 間斷Galerkin方法

間斷Galerkin方法求解方程組(1),首先需要將計(jì)算區(qū)域劃分成互不重疊的子域。子域可以選取為任意形狀,對(duì)于二維空間,本文采用三角形非結(jié)構(gòu)網(wǎng)格,三維空間,采用四面體非結(jié)構(gòu)網(wǎng)格。定義有限元空間Vh:

Vh={vh∈L2(Ω):vh|K∈V(K);?K∈τh}

式中：τh為子域空間,V(K)為局部函數(shù)空間,取作p(p=1,2,…)次多項(xiàng)式的集合。

假設(shè)間斷函數(shù)在有限元空間中的近似解為Uh,將Euler方程(1)兩邊同乘以試驗(yàn)函數(shù)v,寫(xiě)成變分形式,再由Green公式,得到守恒方程組的弱形式:

(2)

式中:Ω為計(jì)算域;?Ω為Ω的邊界。在每個(gè)單元內(nèi):

(3)

式中:φj(x)是基函數(shù)。將(3)式帶入(2)式,得到半離散形式的守恒方程。

(4)

(4)式稱為p階間斷Galerkin有限元離散,p為(3)式中基函數(shù)的最大階數(shù)。不難看出Uh、vh在整個(gè)計(jì)算區(qū)域上不再連續(xù),在單元邊界上是間斷的,而且整個(gè)流場(chǎng)的自由度已經(jīng)轉(zhuǎn)變成求解插值系數(shù)Ui,不再是流場(chǎng)守恒變量。

為了便于計(jì)算,通常還需要進(jìn)行坐標(biāo)變換。針對(duì)三角型單元,我們采用面積坐標(biāo),針對(duì)四面體單元,采用體積坐標(biāo)。將總體笛卡爾坐標(biāo)轉(zhuǎn)換成局部的自然坐標(biāo),將空間中任意單元轉(zhuǎn)變成局部坐標(biāo)系中的標(biāo)準(zhǔn)單元。(4)式變?yōu)?

(5)

式中：|J|為體坐標(biāo)變換雅可比矩陣,|Js|面坐標(biāo)變換雅可比矩陣。

對(duì)于(5)式中的積分,我們選用高斯數(shù)值積分。

3 隱式時(shí)間推進(jìn)格式

TVD-RKDG時(shí)間推進(jìn)格式中,由于受穩(wěn)定性條件限制,計(jì)算中CFL數(shù)過(guò)小,導(dǎo)致效率比較低,而隱式推進(jìn)格式一般是無(wú)條件穩(wěn)定的,計(jì)算步長(zhǎng)可以取得較大,從而整體效率較高。本文在借鑒非結(jié)構(gòu)網(wǎng)格有限體積LU-SGS隱式計(jì)算格式的基礎(chǔ)上,構(gòu)造出了適用于間斷Galerkin有限元法的LU-SGS隱式計(jì)算格式。

LU-SGS方法基本思想是運(yùn)用通量線化假設(shè)和最大特征值方法進(jìn)行雅可比矩陣分裂,把塊對(duì)角矩陣分解為上、下2個(gè)三角矩陣,這種分解可以避免繁雜的矩陣求逆運(yùn)算,極大地提高了計(jì)算效率。該方法本質(zhì)上仍然是一種近似因式分解方法,對(duì)于結(jié)構(gòu)網(wǎng)格比較適合,而在非結(jié)構(gòu)網(wǎng)格下,由于網(wǎng)格存儲(chǔ)的不規(guī)則性,該方法的計(jì)算效率要差一些。

對(duì)公式(5)在時(shí)間方向上采用一階向后歐拉積分:

(6)

將其展開(kāi):

(7)

由于方程(6)為一非線性系統(tǒng),直接求解非常麻煩,對(duì)(7)式中的二、三項(xiàng)分別進(jìn)行一階泰勒展開(kāi),進(jìn)行線化處理:

(8)

式中:

(9)

為了簡(jiǎn)化處理,本文將方程(6)中時(shí)間項(xiàng)和空間項(xiàng)的數(shù)值通量分開(kāi)單獨(dú)處理,在時(shí)間項(xiàng)上無(wú)粘通量Fnum采用LF數(shù)值通量格式,可以分裂為

(10)

式中：λij為雅克比矩陣的譜半徑。

則(9)式可以簡(jiǎn)化成:

將公式(7)、(8)式帶入(6)式,省略高階項(xiàng):

(11)

式中:

將方程(6)代入上式,得到等價(jià)的矩陣形式

(12)

(13)

(14)

式中：L為嚴(yán)格下三角矩陣,D為純對(duì)角線矩陣,U為嚴(yán)格上三角矩陣

(15)

將(15)式變形為

(D+L)D-1(D+U)Δu=Rn+(LD-1U)Δu

(16)

忽略上式右端第二項(xiàng)高階小量,可得:

(D+L)D-1(D+U)Δu=Rn

(17)

上式只需要通過(guò)向前掃描和向后掃描兩個(gè)過(guò)程即可求解自由度,其中

向前掃描:

(18)

向后掃描:

(19)

具體實(shí)現(xiàn)過(guò)程為:

(20)

(21)

4 算例與分析

為了驗(yàn)證本文基于LU-SGS迭代法構(gòu)造的DGM隱式算法的正確性和有效性,本文分別對(duì)繞NACA0012翼型和ONERA M6機(jī)翼的跨聲速無(wú)粘流動(dòng)進(jìn)行數(shù)值模擬,并和TVD-RKDG計(jì)算結(jié)果進(jìn)行比較,來(lái)驗(yàn)證隱式算法的加速收斂效果和精度。本文算例中使用的計(jì)算機(jī)基本配置CPU為酷睿I5 3.2GHz、內(nèi)存為4G,所有網(wǎng)格均采用Delaunay方法生成。

NACA0012翼型的流動(dòng)計(jì)算狀態(tài)為Ma∞=0.8,攻角α=1.25°。整個(gè)流場(chǎng)包括3 531個(gè)網(wǎng)格節(jié)點(diǎn),6 789個(gè)網(wǎng)格單元。圖1給出了翼型表面壓力系數(shù)分布比較,從圖中可以看出,采用DGM隱式算法并沒(méi)有改變精度,并且翼型上下表面的激波均捕捉的很好。圖2分別給出了關(guān)于CPU時(shí)間和迭代步數(shù)的密度最大殘值收斂歷程比較,當(dāng)采用TVD-RKDG時(shí)計(jì)算花費(fèi)了116 min左右,采用隱式算法后,計(jì)算花費(fèi)了14 min左右,計(jì)算效率提高8倍左右。

圖1 翼型表面壓力系數(shù)分布比較圖 圖2 殘值收斂歷程比較

ONERA M6機(jī)翼的流動(dòng)計(jì)算狀態(tài)為Ma∞=0.84,攻角α=3.06°。整個(gè)流場(chǎng)包含65 115個(gè)網(wǎng)格節(jié)點(diǎn)和459 353個(gè)網(wǎng)格單元。圖3顯示了機(jī)翼上表面壓力等直線圖,圖中λ狀激波結(jié)構(gòu)捕捉非常清晰,外弦和內(nèi)弦激波大約在展長(zhǎng)87%處相交,在94%處又分開(kāi),和實(shí)驗(yàn)結(jié)果[15]基本吻合。圖4給出了M6機(jī)翼65%展長(zhǎng)處的剖面壓力系數(shù)分布比較,結(jié)果和實(shí)驗(yàn)值吻合較好,并且DG隱式算法和TVD-RKDG得到的壓力系數(shù)分布完全重合。圖5分別給出了關(guān)于CPU時(shí)間和密度最大殘值收斂歷程比較,從圖中可見(jiàn),當(dāng)采用TVD-RKDG時(shí)計(jì)算花費(fèi)了2 300 min左右,采用隱式算法后,計(jì)算花費(fèi)了500 min左右,計(jì)算效率提高5倍左右。

圖3 M6機(jī)翼上表面壓力等值線

圖4 65 %展向位置處剖面壓力系數(shù)分布的對(duì)比 圖5 殘值收斂歷程比較

5 結(jié) 論

本文主要構(gòu)造了適合于DGM隱式求解的LU-SGS迭代格式,并發(fā)展了一套高效、實(shí)用的非結(jié)構(gòu)網(wǎng)格Euler方程求解程序。通過(guò)數(shù)值模擬跨聲速Euler方程中來(lái)驗(yàn)證其效率和精度。數(shù)值算例表明:采用DGM隱式算法能夠很好的維持DGM的精度,并能明顯的提高流場(chǎng)的收斂速度以及降低計(jì)算成本。在本文的二維算例中,計(jì)算效率達(dá)到8倍左右,三維算例中,計(jì)算效率達(dá)到5倍左右。雖然,本文的工作是基于DG(P1),但由于公式推導(dǎo)具有普遍性,可以比較容易推廣到更高精度。因此,今后工作還需要進(jìn)一步考察本文構(gòu)造的隱式算法在更高精度DGM中的性能。

參考文獻(xiàn)：

[1] Reed N H, Hill T R. Triangle Mesh Methods for the Neutron Transport Equation[R]. Los Almos Scientific Laboratory, Report LA-UR-73-479, 1973

[2] Fidkowski K J, Oliver T A, Lu J, Darmofal D L. P-Multigrid Solution of High-Order Discontinuous Galerkin Discretizations of the Compressible Navier-Stokes Equations[J]. J Comp Phys, 2005, 207: 92-113

[3] Qiu J, Shu C W. Hermite WENO Schemes and Their Application as Limiters for Runge-Kutta Discontinuous Galerkin Method: One Dimensional Case [J]. J Comp Phys, 2004, 193: 115-135

[4] Zhang L P, Liu W, He L X, Deng X G. A New Class of DG/FV Hybrid Schemes for One-Dimensional Conservation Law[C]∥The 8th Asian Conference on Computational Fluid Dynamics, Hong Kong, 2010: 10-14

[5] 郝海兵,楊永. 非結(jié)構(gòu)網(wǎng)格上P型多重網(wǎng)格法流場(chǎng)數(shù)值模擬[J]. 計(jì)算力學(xué)學(xué)報(bào), 2011, 28(3): 360-365

Hao Haibing, Yang Yong. The Research of P-Multigrid Solution for Discontinuous Galerkin Method [J]. Chinese Journal of Computational Mechanics, 2011, 28(3): 360-365 (in Chinese)

[6] Bassi F, Rebay S. GMRES for Discontinuous Galerkin Solution of the Compressible Navier-Stokes Equations∥Cockburn B, Karniadakis G E, Shu C W, Discontinuous Galerkin Method: Theory, Computations and Applications[M], Springer-Verlag, 2000

[7] Hartmann R, Houston P. Symmetric Interior Penalty DG Methods for the Compressible Navier-Stokes Equations I: Method Formulation[J]. Internatimal Journal of Numerical Analysis and Modeling, 2006(3): 1-20