馬明生, 龔小權(quán), 鄧有奇, 趙 輝

1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000

?

一種適用于非結(jié)構(gòu)網(wǎng)格的間斷Galerkin有限元LU-SGS隱式方法

馬明生1,2, 龔小權(quán)1, 鄧有奇1, 趙 輝1

1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000

具有TVD性質(zhì)的顯式Runge-Kutta間斷Galerkin(RKDG)格式在CFD領(lǐng)域得到廣泛應(yīng)用,但是顯式計(jì)算穩(wěn)定性差、計(jì)算效率低。為改善時(shí)間推進(jìn)效率,基于高階間斷Galerkin有限元方法,采用歐拉一階后差(BDF1),發(fā)展了一套高效的隱式LU-SGS(lower upper-symmetric Gauss-Seidel)求解方法,方法基于MPI并行實(shí)現(xiàn),適合于不同計(jì)算精度。針對(duì)非線性系統(tǒng)左端項(xiàng)矩陣,對(duì)比了簡(jiǎn)化前后LU-SGS的計(jì)算效率。建立的間斷Galerkin有限元方法基于非結(jié)構(gòu)網(wǎng)格,采用Taylor基函數(shù),計(jì)算精度最高達(dá)到四階精度。通過(guò)NACA0012翼型以及M6機(jī)翼算例對(duì)發(fā)展的LU-SGS方法進(jìn)行了考察,與顯式算法相比,隱式格式的迭代步數(shù)和CPU時(shí)間均較大程度減小,效率能夠提高1個(gè)量級(jí)以上。最后將隱式算法用于復(fù)雜外形翼身組合體F4的流場(chǎng)計(jì)算,結(jié)果表明所發(fā)展的隱式方法具有較好的魯棒性,能夠用于復(fù)雜外形計(jì)算。

間斷Galerkin有限元;歐拉方程;Taylor基函數(shù);LU-SGS;計(jì)算效率;非結(jié)構(gòu)網(wǎng)格

隨著CFD和計(jì)算機(jī)技術(shù)的發(fā)展,CFD工作者希望得到更精細(xì)的流場(chǎng)結(jié)構(gòu)以及更準(zhǔn)確的氣動(dòng)數(shù)據(jù),比如阻力、力矩、熱流。精細(xì)化的網(wǎng)格或者高精度的計(jì)算格式是實(shí)現(xiàn)這些要求的途徑之一。在規(guī)則網(wǎng)格上進(jìn)行光滑函數(shù)的數(shù)值近似,n階精度的計(jì)算格式,計(jì)算誤差與hn(h為網(wǎng)格尺寸)成正比。高精度格式能夠給出一些復(fù)雜流動(dòng)現(xiàn)象的精細(xì)流場(chǎng)結(jié)構(gòu),相比傳統(tǒng)二階格式具有更高的分辨率。而非結(jié)構(gòu)網(wǎng)格相比結(jié)構(gòu)網(wǎng)格具有更大的靈活性,更適應(yīng)復(fù)雜外形和邊界條件。因此發(fā)展適用于非結(jié)構(gòu)網(wǎng)格的高精度格式一直是計(jì)算流體力學(xué)(CFD)研究的熱點(diǎn)。

間斷Galerkin有限元方法(discontinuous galerkin finite element method,DGM)是一種適合于結(jié)構(gòu)網(wǎng)格以及非結(jié)構(gòu)網(wǎng)格的高精度格式。它最早出現(xiàn)在1973年Reed和Hill[1]求解中子輸運(yùn)方程問(wèn)題的論文中。DGM結(jié)合了有限體積法和有限元法的優(yōu)點(diǎn)。作為一種有限元方法,它采用有限元法中基函數(shù)思想來(lái)構(gòu)造單元內(nèi)的變量分布,容易實(shí)現(xiàn)高精度。同時(shí)它又具有有限體積法積分形式的計(jì)算形式且融入了有限體積法中數(shù)值通量、Riemann間斷、限制器等思想。DGM具有提高計(jì)算精度容易,高精度時(shí)計(jì)算模板緊致,適合非結(jié)構(gòu)網(wǎng)格,適合處理復(fù)雜外形及復(fù)雜邊界、初值問(wèn)題,并行能力強(qiáng),自適應(yīng)計(jì)算方便等優(yōu)點(diǎn)。間斷Galerkin有限元方法[2]在求解守恒方程組問(wèn)題中得到廣泛應(yīng)用。經(jīng)過(guò)Cockburn和Shu等的持續(xù)研究,一種具有TVD性質(zhì)的顯式Runge-Kutta間斷Galerkin(RKDG)格式得以逐步完善。

RKDG格式單步計(jì)算時(shí)間少,內(nèi)存需求量相對(duì)較小,程序容易實(shí)現(xiàn),但是隨著計(jì)算精度的提高, RKDG格式對(duì)應(yīng)的穩(wěn)定性條件將越來(lái)越嚴(yán)格,時(shí)間推進(jìn)步長(zhǎng)受到明顯限制,計(jì)算效率低下,特別是黏性計(jì)算。一些加速技術(shù)如當(dāng)?shù)貢r(shí)間步長(zhǎng)、殘值光順、多重網(wǎng)格等,在一定程度上加速了收斂過(guò)程。即便如此,最大的時(shí)間步長(zhǎng)依然受到當(dāng)?shù)胤€(wěn)定性條件的限制。與顯式的時(shí)間格式相比,隱式方法基本不受時(shí)間推進(jìn)步長(zhǎng)的限制,因此發(fā)展高精度間斷Galerkin有限元方法的時(shí)間隱式離散方法,有著重要的工程應(yīng)用價(jià)值。

目前間斷Galerkin有限元方法已有多種高效的隱式離散方法,比如帶預(yù)處理的GMRES[3]、二階的Crank-Nicholson(CN2)格式[4]、四階隱式的Runge-Kutta(IRK4)等。LU-SGS隱式方法因?yàn)槠湓谟?jì)算效率、魯棒性及內(nèi)存占用方面的優(yōu)勢(shì),在有限體積方法(FVM)中被大量采用。Yoon與Jameson在1988年[5]首先將其作為迭代方法用于方程求解,隨后Rieger和Jameson將其用于求解三維黏性流場(chǎng),從此各種改進(jìn)的LU-SGS隱式方法被用來(lái)求解基于結(jié)構(gòu)及非結(jié)構(gòu)網(wǎng)格的流場(chǎng)。詳細(xì)介紹LU-SGS在高精度間斷Galerkin有限元方法中應(yīng)用的文獻(xiàn)較少,Haga和Wang將該方法用于基于非結(jié)構(gòu)四面體網(wǎng)格的譜方法[6]。郝海兵、張強(qiáng)等在DGM四面體網(wǎng)格上實(shí)現(xiàn)了p=1階的LU-SGS迭代[7]。郭永恒在其博士論文中[8]研究了不同LU-SGS改進(jìn)方法對(duì)收斂的影響。劉偉在文獻(xiàn)[9]中發(fā)展了基于Newton/Gauss-Seidel迭代的DGM隱式方法。

本文基于建立在非結(jié)構(gòu)網(wǎng)格上的高精度間斷Galerkin有限元方法,發(fā)展了LU-SGS隱式計(jì)算方法,比較了不同計(jì)算精度下的收斂曲線,通過(guò)流場(chǎng)熵增云圖給出了不同精度下的計(jì)算結(jié)果。同時(shí)對(duì)左端項(xiàng)進(jìn)行了簡(jiǎn)化,給出了左端項(xiàng)不同簡(jiǎn)化方法對(duì)收斂性的影響。給出了一系列二維、三維算例不同精度下的收斂曲線,比較了二階精度下有限體積和DGM下的收斂速度。

1 數(shù)值方法

1.1 控制方程

守恒形式的無(wú)黏可壓縮歐拉方程組可寫(xiě)為

(1)

其中守恒變量、無(wú)黏對(duì)流通量為

(2)

ρ、u、v、w、p、E、H分別代表密度、3個(gè)方向速度、壓力、總能和總焓。

(3)

γ為比熱比。

1.2 DGM空間離散

對(duì)方程組(1)采用間斷Galerkin有限元離散。將其兩端同時(shí)乘以測(cè)試函數(shù)φi,并對(duì)整個(gè)求解域積分Ω,再采用分部積分得到如下公式

(4)

式中，?Ω為求解域Ω的邊界,n為邊界的法向量。將求解域Ω剖分為互不重疊的非結(jié)構(gòu)網(wǎng)格單元,并對(duì)每個(gè)單元采用上面的弱形式。每個(gè)計(jì)算單元的每個(gè)守恒變量在單元內(nèi)的函數(shù)分布為

(5)

uj為單元守恒變量的自由度,即為每時(shí)間步需要求解的量。N=N(p,d)由插值多項(xiàng)式階數(shù)p和求解的空間維數(shù)d決定。φj是階數(shù)為p的基函數(shù)。由于Taylor基函數(shù)[10]對(duì)所有單元形式都一樣且為等級(jí)基,對(duì)非結(jié)構(gòu)混合網(wǎng)格較為合適,因此本文采用Taylor基函數(shù)。將方程(5)帶入方程(4)得

(6)

最后整個(gè)求解域方程可寫(xiě)為如下形式

(7)

Res為無(wú)黏通量殘差項(xiàng),由面積分和體積分構(gòu)成,二者均采用高斯數(shù)值積分計(jì)算。M為全局的質(zhì)量矩陣,由各單元的質(zhì)量矩陣構(gòu)成。無(wú)黏通量殘差項(xiàng)具體計(jì)算為

(8)

wv(l)、|Jl|、φil分別為體積分對(duì)應(yīng)的單元高斯積分點(diǎn)的權(quán)系數(shù)、雅克比矩陣行列式值以及第i個(gè)基函數(shù)的梯度。ws(k)、|Jsk|、φis分別為單元面面積分對(duì)應(yīng)的高斯積分點(diǎn)k的權(quán)系數(shù)、雅克比矩陣行列式值以及第i個(gè)基函數(shù)的值。Wak、Wbk分別對(duì)應(yīng)邊界面上當(dāng)前單元和相鄰單元在積分點(diǎn)k的變量值,Wal為當(dāng)前單元體積分點(diǎn)l的變量值。

1.3 隱式LU-SGS方法

對(duì)方程(7)第一式采用一階歐拉后差得到

(9)

方程(9)為一大型的線性系統(tǒng)Ax=B,其中矩陣A為一大型稀疏塊矩陣。以三維p=3(四階精度)DGM為例,矩陣A的維數(shù)為n-cell×(5×20),n-cell為計(jì)算單元總數(shù)。如果直接存儲(chǔ)該矩陣并求逆,其存儲(chǔ)量和計(jì)算量太大,目前大多采用迭代法來(lái)求解該線性系統(tǒng)。

將面積分和體積分帶入左端項(xiàng)的Resn。

(10)

將左端項(xiàng)矩陣寫(xiě)為

(11)

式中，D為對(duì)角塊矩陣,L為嚴(yán)格下三角矩陣,U為嚴(yán)格上三角矩陣。

為方便推導(dǎo),進(jìn)一步得到矩陣D、L、U,本文只寫(xiě)出左端項(xiàng)矩陣中與當(dāng)前計(jì)算單元有關(guān)的塊,將當(dāng)前單元編號(hào)用a表示,相鄰單元用b表示,

對(duì)(10)式采用鏈?zhǔn)椒▌t

(12)

(13)

進(jìn)一步得到面積分通量對(duì)守恒變量的偏導(dǎo)數(shù)

(14)

將方程(8)的第一項(xiàng)以及方程(14)帶入方程(12)得

(15)

(16)

(17)

2 數(shù)值模擬結(jié)果及分析

本文首先通過(guò)等熵渦算例驗(yàn)證了DGM的計(jì)算精度。為驗(yàn)證本文發(fā)展的LU-SGS方法的收斂性,計(jì)算了幾類典型的測(cè)試算例。測(cè)試過(guò)程中,右端對(duì)流通量項(xiàng)采用Roe′s格式,并采用三階TVD-Runge-Kutta方法來(lái)對(duì)比分析。利用NACA0012翼型亞聲速算例考察了不同精度的LU-SGS收斂曲線,研究了CFL數(shù)對(duì)收斂歷程的影響,給出了左端項(xiàng)簡(jiǎn)化前后的收斂曲線對(duì)比。計(jì)算了三維M6機(jī)翼亞聲速繞流,給出了不同精度的收斂曲線。最后通過(guò)F4翼身組合體給出了發(fā)展的LU-SGS隱式方法在復(fù)雜外形中的應(yīng)用。在分析過(guò)程中,本文采用成熟且經(jīng)過(guò)驗(yàn)證的有限體積方法程序作為對(duì)比分析的參考。

2.1 等熵渦問(wèn)題

等熵渦問(wèn)題有精確解,用來(lái)對(duì)二維Euler問(wèn)題格式的精度進(jìn)行驗(yàn)證,這里把問(wèn)題擴(kuò)展到三維問(wèn)題。等熵渦問(wèn)題描述的是在初始的平均流(ρ,u,v,w,p)=(1,1,1,0,1)上加入一個(gè)等熵渦,渦的大小按照(18)式給出。

(18)

式中，T=p/ρ,熵S=p/ρr,ε為渦的強(qiáng)度ε=5.0;計(jì)算區(qū)域取為[0,10]×[0,10]×[0,10],邊界條件按照精確解給定,時(shí)間上采用顯式三階Runge-Kutta法推進(jìn)至2 s。計(jì)算采用的網(wǎng)格為N×N×N(N=10,20,40)六面體網(wǎng)格。表1給出了不同精度下的L1 Error和L2 Error,可以看到方法達(dá)到了設(shè)計(jì)的精度。

表1 等熵渦問(wèn)題精度驗(yàn)證(六面體網(wǎng)格)

2.2 NACA0012翼型亞聲速繞流

為考察DGM的精度以及發(fā)展的隱式算法的收斂性,本文計(jì)算了NACA0012翼型在來(lái)流馬赫數(shù)Ma=0.5時(shí)亞聲速繞流。計(jì)算網(wǎng)格如圖1所示。

圖1 NACA0012計(jì)算網(wǎng)格

從圖3a)的曲線1和5,CFL=1,簡(jiǎn)化左端項(xiàng)矩陣對(duì)收斂速度有一定影響,導(dǎo)致迭代步數(shù)略有增加。從圖3b)看到,對(duì)于CFL=1的情況,由于簡(jiǎn)化左端項(xiàng)矩陣后,單步計(jì)算時(shí)間減小,總的收斂時(shí)間減小。由于簡(jiǎn)化左端項(xiàng)矩陣后CFL數(shù)不能取大,因此未簡(jiǎn)化左端項(xiàng)矩陣當(dāng)CFL數(shù)取大時(shí),其收斂速度快于簡(jiǎn)化后的狀態(tài)。

圖2 DGM不同計(jì)算精度及二階FVM收斂歷程 圖3 DGM二階精度不同CFL數(shù)收斂歷程

2.3 ONERA M6機(jī)翼亞聲速繞流

為進(jìn)一步考察LU-SGS方法在三維外形計(jì)算效率,本文計(jì)算了ONERA-M6亞聲速流場(chǎng),計(jì)算馬赫數(shù)Ma=0.5,迎角0°。網(wǎng)格如圖4所示,網(wǎng)格共約14萬(wàn)個(gè)四面體單元,機(jī)翼前后緣及空間分別采用各向異性和各向同性四面體網(wǎng)格。計(jì)算采用64個(gè)CPU。

圖4 M6機(jī)翼計(jì)算網(wǎng)格

圖5給出了DGM不同計(jì)算精度下LU-SGS殘差收斂曲線,同時(shí)給出了FVM以及RKDG顯式結(jié)果。隱式計(jì)算的CFL=1 000,顯式計(jì)算的CFL=0.3。從殘差隨迭代步數(shù)圖看,發(fā)展的LU-SGS隱式方法收斂效率遠(yuǎn)高于顯式方法。二階精度下,由于DGM計(jì)算的自由度為FVM的4倍且計(jì)算面通量和體積分通量時(shí)采用更多的高斯積分點(diǎn),二階DGM的單步計(jì)算時(shí)間遠(yuǎn)大于FVM。

圖5 不同計(jì)算精度DGM及二階FVM收斂曲線

2.4 DLR-F4翼身組合體亞聲速繞流

為驗(yàn)證發(fā)展的隱式算法的初步工程計(jì)算能力,數(shù)值模擬了DLR-F4翼身組合體亞聲速繞流。計(jì)算網(wǎng)格如圖6所示。

圖6 F4表面網(wǎng)格

四面體網(wǎng)格單元總數(shù)為637 126個(gè),計(jì)算馬赫數(shù)Ma=0.5,來(lái)流迎角為0°,計(jì)算采用128個(gè)CPU,圖7給出了殘差收斂曲線,相同CFL數(shù)以及相同的dt計(jì)算方法下,相對(duì)于有線體積方法,DGM殘差收斂速度較快,但是下降量級(jí)較少。這是由于在迭代過(guò)程中DGM一直修改迭代得到的自由度,造成限制單元?dú)埐畈唤?因此需要進(jìn)一步研究該方法的限制器。

圖7 不同方法及精度的殘差收斂曲線

3 結(jié) 論

本文發(fā)展了一種高效的LU-SGS隱式計(jì)算方法,方法適用于不同計(jì)算精度,通過(guò)一系列的算例表明:

1) 本文發(fā)展的隱式方法CFL數(shù)能夠取很大。相比顯式RKDG,其計(jì)算效率提高了1～2個(gè)量級(jí)。

2) 發(fā)展的隱式算法對(duì)左端項(xiàng)殘差比較敏感,需小心簡(jiǎn)化左端項(xiàng)殘差,簡(jiǎn)化可能限制CFL數(shù),導(dǎo)致整個(gè)計(jì)算時(shí)間增加,計(jì)算效率降低。

3) 隱式計(jì)算時(shí),單元左端項(xiàng)矩陣包含時(shí)間項(xiàng)、面積分項(xiàng)和體積分項(xiàng),單元對(duì)應(yīng)的矩陣D的矩陣元素幾乎都不為零,降低了CFL數(shù)對(duì)收斂性的影響。

4) 發(fā)展的隱式方法能夠初步應(yīng)用于復(fù)雜外形計(jì)算。

[1] Reed W H, Hill T R. Triangular Mesh Methods for the Neu-Tron Transport Equation[R]. Technical Report LA-UR-73-479, Los Alamos Scientific Lab, 1973

[2] Cockburn B, Shu Chiwang. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws Ⅳ: The Multidimensional Case[J]. Mathematics of Computation, 1990, 54: 545-581

[3] Persson Per-Olof, Peraire Jaime. An Efficient Low Memory Implicit DG Algorithm for Time Dependent Problems[C]∥44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2006: 9-12

[4] Wang L, Mavriplis D J. Implicit Solution of the Unsteady Euler Equation for High-Order Accurate Discontinuous Galerkin Discretizations[J]. J Comput Phys, 2007, 225:1994-2005

[5] Yoon S, Jameson A. Lower-Upper Implicit Schemes with Multiple Grids for the Euler Equatins[R]. AIAA-1986-0105

[6] Hagal Takanori, Sawada Keisuke, Wang Z J. An Implicit LU-SGS Scheme for the Spectral Volume Method on Unstructured Tetrahedral Grids[J]. Commun Comput Phys, 2009,6(5):978-996

[7] 郝海兵,張強(qiáng),楊永. 基于LU-SGS迭代的DGM隱式算法研究[J]. 西北工業(yè)大學(xué)學(xué)報(bào),2014,32(3):346-350

Hao Haibing, Zhang Qiang, Yang Yong. An Implicit Scheme of Discontinuous Galerkin Method(DGM) Based on LU-SGS Iterative Method[J]. Journal of Northwestern Polytechnical University, 2014,32(3):346-350 (in Chinese)

[8] 郭永恒. 雙曲守恒律方程組的間斷Galerkin方法研究[D]. 西安:西北工業(yè)大學(xué),2013

Guo Yongheng. Research of Discontinuous Galerkin Method for Hyperbolic Conservation Equations[D]. Xi′an, Northwestern Polytechnical University,2013 (in Chinese)

[9] 劉偉,張來(lái)平,赫新,等. 基于Newton/Gauss-Seidel迭代的 DGM隱式方法[J]. 力學(xué)學(xué)報(bào),2012,44(4):792-796

Liu Wei, Zhang Laiping, He Xin, et al. An Implicit Algorithm for Discontinuous Galerkin Method Based on Newton/Gauss-Seidel Iterations[J]. Acta Mechanica Sinica, 2012,44(4): 792-796 (in Chinese)

[10] Hong L, Baum J D, Lohner R. A Discountinuous Galerkin Method Based on a Taylor Basis for The Compressible Flows on Arbitrary Grids[J]. J Comput Phys, 2008, 227: 8875-8893

An Implicit LU-SGS Scheme for the Discontinuous Galerkin Method on Unstructured Grids

Ma Mingsheng1,2,Gong Xiaoquan1, Deng Youqi1, Zhao Hui1

1.School of Aeronautics, Northwestern Polytechnical University,Xi′an 710072, China 2.China Aerodynamics Research and Development Center,Mianyang 621000, China

The TVD explicit Runge-Kutta discontinuous Galerkin(TVD-RKDG) are widely used in CFD, but it has poor stability property and low computational efficiency. To improve the time advancing efficiency, an efficient implicit lower-upper symmetric Gauss-Seidel (LU-SGS) scheme based on backwards differencing methods (BDF1) has been applied to the high-order discontinuous Galerkin method. The method is implemented parallelly through MPI library and is suit for different computational orders. The left matrix of nonlinear system is simplified and computational efficiency is compared. The DGM is based on the Taylor basis functions on unstructured grid and the accuracy is up to forth order. The developed LU-SGS scheme is verified through NACA0012 airfoil and ONERA M6 wing. Compared to the traditional explicit Runge-Kutta scheme, the implicit scheme can greatly reduce the iteration number and CPU time. The computational efficiency can be accelerated evidently with more than one order of magnitude speed. At last the flow field of DLR-F4 is simulated to verify the robustness of the developed implicit scheme. The result proves that the developed LU-SGS scheme can be applied to complex configuration simulation.

discontinuous Galerkin method; Euler equations; Taylor basis; LU-SGS; computational efficiency; unstructured grid

2016-03-20

馬明生(1962—)，西北工業(yè)大學(xué)兼職教授，主要從事飛行器設(shè)計(jì)與計(jì)算流體力學(xué)研究。

V211.3

A

1000-2758(2016)05-0754-07