歐拉方程的隱式間斷有限元算法研究

段治健

1.咸陽師范學(xué)院數(shù)學(xué)與信息科學(xué)學(xué)院，陜西咸陽 712000

2.西北工業(yè)大學(xué)翼型葉柵空氣動力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，西安 710072

歐拉方程的隱式間斷有限元算法研究

段治健1，2

1.咸陽師范學(xué)院數(shù)學(xué)與信息科學(xué)學(xué)院，陜西咸陽 712000

2.西北工業(yè)大學(xué)翼型葉柵空氣動力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，西安 710072

針對Euler方程，設(shè)計了適合間斷Galerkin有限元方法的LU-SGS、GMRES以及修正LU-SGS隱式算法。采用Roe通量以及Van A lbada限制器技術(shù)實(shí)現(xiàn)了經(jīng)典LU-SGS、GMRES算法，引入高階項誤差補(bǔ)償，發(fā)展了修正LU-SGS算法。以NACA 0012、RAE2822翼型為例驗(yàn)證分析了算法的可靠性和高效性。結(jié)果表明修正LU-SGS算法存儲量較少，程序?qū)崿F(xiàn)方便，而且計算效率是LU-SGS算法的2.5倍以上，接近于循環(huán)GMRES算法。

間斷有限元方法；LU-SGS算法；GMRES算法；限制器

1 引言

近年來，間斷有限元方法（Discontinuous Galerkin FEM，DGM）[1]在解決含有間斷現(xiàn)象的問題中發(fā)揮著越來越大的作用，由于DGM保持了有限元法和有限差分法的優(yōu)點(diǎn)，可以處理復(fù)雜區(qū)域邊界、復(fù)雜邊界條件，易于實(shí)現(xiàn)自適應(yīng)計算；可以得到任意階精度的格式；容易實(shí)現(xiàn)并行計算；具有良好的穩(wěn)定性和間斷捕捉能力，所以被廣泛應(yīng)用于水動力學(xué)、氣動力學(xué)、波傳播等問題[1-4]。盡管DGM具有高精度、收斂快的優(yōu)點(diǎn)，但對大型問題而言計算效率是一個瓶頸，然而隱式算法是提高計算效率最重要的途徑。鑒于隱式格式的高效率和較好的穩(wěn)定性，以及近年來非結(jié)構(gòu)網(wǎng)格技術(shù)的蓬勃發(fā)展，許多高性能的隱式算法在計算空氣動力學(xué)領(lǐng)域中脫穎而出。隱式時間推進(jìn)格式往往都是無條件穩(wěn)定的，相對于顯式格式，時間步長約束小、總體計算時間短。著名的隱式方法有：ADI方法、LU-SGS方法[5]、GMRES方法[6]等，并得到了進(jìn)一步發(fā)展[7-8]。相比之下，LU-SGS格式占用內(nèi)存低，計算量小，每步時間推進(jìn)只需要從前到后、從后到前掃描計算即可。

本文研究了LU-SGS、GMRES、修正LU-SGS隱式時間推進(jìn)格式，針對NACA 0012、RAE2822翼型比較了三種算法的優(yōu)劣性，發(fā)展了一套高效的雙曲守恒律隱式間斷有限元方法。

2 間斷有限元方法

不考慮體積力和外部熱源，笛卡爾坐標(biāo)系下的Euler方程形式如下：

其中，U為守恒變量，F(xiàn)j為對應(yīng)無粘通量，ρ，e，p分別為氣體的密度、單位體積總內(nèi)能及壓強(qiáng)。ui是xi方向的速度，δij是K ronecker函數(shù)?？紤]到理想氣體的熱力學(xué)特性，狀態(tài)方程如下：

取γ=1.4為比熱比。

在式（1）兩邊乘以測試函數(shù)W，在計算域內(nèi)運(yùn)用分部積分公式，得到弱形式：

這里Γ是Ω的邊界，nj是邊界外法向量。

在計算域劃分網(wǎng)格上，用近似解Uh，Wh代替解析解U，W，每個單元內(nèi)：

因?yàn)殚g斷有限元方法在區(qū)域邊界不要求連續(xù)，這樣數(shù)值通量處理方法與有限體積方法相同，如Godunov、Engquist-Osher、HLL/HLLC、LF、LLF、Roe數(shù)值通量等。本文采用Van A lbada限制器[9]，數(shù)值通量采用Roe格式，該格式相對簡單，實(shí)際計算會更加方便，數(shù)值效果良好。

3 GMRES方法

具有平方斂速的共軛梯度法[10]是New ton迭代隱式算法的一種，GMRES算法屬于其中一種，因而倍受關(guān)注，循環(huán)GMRES算法相對于GMRES算法收斂性通常差一些，但是存儲量相對較少。因此本文采用循環(huán)GMRES算法。對于給定的非對稱線性方程組Ax=f，循環(huán)GMRES算法如下：

給定K rylov子空間：Km=Span{r0，A r0，…，Am-1r0}，而Lm=A Km=Span{A r0，A2r0，…，Amr0}，給定誤差ε。

（1）選擇初始解x0∈Rn，計算殘差r0=f-Ax0。

4 LU-SGS方法

當(dāng)式（4）中取為常數(shù)時，則方程變?yōu)椋?/p>

其中，Vi為網(wǎng)格單元i的體積，Δt為時間步長，為第n+1層上的殘值，F(xiàn)ij為單元i和單元j的公共網(wǎng)格面上的無粘通量，n為網(wǎng)格面的外法線單位矢量，ΔS為網(wǎng)格面的面積。將式（5）泰勒展開，線化為：

從而得到一個N維的線性方程組（N為網(wǎng)格總數(shù)）：

將系數(shù)矩陣分解為A=D+L+U，式（7）變?yōu)椋?/p>

其中，D為對角矩陣，U，L為嚴(yán)格上、下三角矩陣。忽略高階小量(LD-1U)ΔQn，方程只需兩次掃描便可完成求解。忽略了高階小量雖不會影響精度，但增加了截斷誤差會影響收斂速度，所以可以將高階小量進(jìn)行補(bǔ)償。具體計算過程如下：

上式計算過程同（1）、（2）。

5 數(shù)值算例及結(jié)果分析

算例1 NACA 0012翼型跨音速無粘繞流，計算網(wǎng)格節(jié)點(diǎn)2 270個，單元數(shù)為4 355個。計算狀態(tài)為Ma=0.8，α=1.25°。

算例2 RAE2822翼型跨音速無粘繞流，計算網(wǎng)格節(jié)點(diǎn)2 270個，單元數(shù)為4 355個。計算狀態(tài)為Ma=0.725，α=2.54°。

圖1和圖2分別為NACA 0012和RAE2822翼型的壓力系數(shù)分布曲線，數(shù)值計算的結(jié)果與實(shí)驗(yàn)值基本吻合。

圖1 NACA 0012翼型表面壓力系數(shù)對比

圖2 RAE2822翼型表面壓力系數(shù)對比

圖3和圖4分別為兩個算例時間收斂曲線?？梢钥闯觯h(huán)GMRES算法計算效率遠(yuǎn)高于LU-SGS格式。計算效率是LU-SGS格式的3倍左右。修正LU-SGS格式效率接近于循環(huán)GMRES算法，且所需存儲量較少，單步迭代時間短，計算效率是LU-SGS格式的2.5倍以上。

圖3 NACA 0012翼型計算時間曲線

圖4 RAE2822翼型計算時間曲線

圖5 NACA 0012翼型隨CFL數(shù)變化計算時間曲線

圖5和圖6中明顯可以看出，克服了顯式方法嚴(yán)格的條件限制，CFL數(shù)從1變到10時，殘值與時間效率提高了近2倍，從100到1 000時，幾乎吻合，此時CFL數(shù)對于計算效率的影響幾乎為0?？傮w來說，LU-SGS格式的單步計算量最少，收斂迭代步數(shù)最多。循環(huán)GMRES格式單步計算時間最長，收斂所需的步數(shù)最少。修正LU-SGS格式計算效率接近循環(huán)GMRES格式，且所需存儲量少。

圖6 NACA 0012翼型隨CFL數(shù)變化迭代步數(shù)收斂曲線

6 結(jié)束語

本文研究了歐拉方程的三種隱式間斷有限元算法，顯然GMRES格式效率最高，但是其算法復(fù)雜，編程實(shí)現(xiàn)比較困難，而且對內(nèi)存的需求較大。修正LU-SGS格式明顯優(yōu)于傳統(tǒng)LU-SGS格式，計算效率接近于GMRES算法，顯示出了其良好的穩(wěn)定性和求解效率。當(dāng)然，線性方程組的并行處理技術(shù)[11-12]，近年來發(fā)展的p、hp多重網(wǎng)格方法[13-14]，或者系數(shù)矩陣條件數(shù)較大時，采用ILU（0）、Block-ILU（0）[15]、Gauss-Seidel、SSOR預(yù)處理方法等，都可以很好地提高計算效率。

[1]Reed W H，Hill T R.Triangular mesh methods for the Neutron Transport equation，LA-UR-73-479[R].Los Alamos Scientific Laboratory，1973.

[2]Cockburn B，Shu C-W.Foreword for the special issue on discontinuous Galerkin method[J].Journal of Scientific Computing，2005：22-23.

[3]Luo H，Beaum J D，Lohner R.On the computation of steadystate compressible flow s using a discontinuous Galerkin method[J].International Journal for Numerical Methods in Engineering，2008，73：597-623.

[4]Qiu J X，Liu T，Khoo B C.Runge-Kutta discontinuous Galerkin methods for compressible two-medium flow simulations：One-dimensional case[J].Journal of Computational Physics，2007，222：353-373.

[5]Yoon S，Jameson A.Lower-upper symmetric Gauss-Seidel method for the Euler and Navier-Stoker equations[J]. AIAA Journal，1988，26（9）：1025-1026.

[6]Saad Y，Schultz M H.A generalized m inimal residual algorithm for solving nonsymmetric linear system s[J]. SIAM Journal on Scientific and Statistical Computing，1986，7：856-869.

[7]李勁杰，楊青，楊永年.三維非結(jié)構(gòu)網(wǎng)格Euler方程的LU-SGS算法及其改進(jìn)[J].計算物理，2006，23（6）：748-752.

[8]李春娜，葉正寅.基于二維非結(jié)構(gòu)網(wǎng)格的GMRES隱式算法[J].西北工業(yè)大學(xué)學(xué)報，2007，25（5）：630-635.

[9]Jawahar P，Kamath H.A high-resolution procedure for Euler and Navier-Stokes computations on unstructured grids[J]. J Comput Phys，2000，164：165-203.

[10]Orkw is P D，George J H.A comparison of CGS preconditioning methods for New ton’s method solvers[C]//AIAA，1993.

[11]段治健，楊永，馬欣榮，等.求解帶狀線性方程組的一種并行算法[J].計算機(jī)科學(xué)，2010，37（3）：242-244.

[12]曹芳芳，呂全義.解非對稱塊三對角線性方程組的并行算法[J].西北工業(yè)大學(xué)學(xué)報，2011，29（2）：318-322.

[13]K rzysztof J F，Todd A O，James L，et al.P-Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations[J]. Journal of Computational Physics，2005，207：92-113.

[14]Cristian R N，Dim itri J M.H-order discontinuous Galerkin methods using an hp-multigrid approach[J].Journal of Computational Physics，2006，213：330-357.

[15]Laslo T D，David L D.Preconditioning methods for discontinuous Galerkin solutions of the Navier-Stokes equations[J].Journal of Computational Phsics，2009，228：3917-3935.

DUAN Zhijian1，2

1.Department of Mathematics,Xianyang Normal University,Xianyang,Shaanxi 712000,China
2.National Key Laboratory of Aerodynamic Design and Research,Northwestern Polytechnical University,Xi’an 710072,China

The implicit schemes for solving Euler equation are investigated on unstructured grids,including LU-SGS, GMRES and improved LU-SGS schemes.By Roe numerical flux and Van Albada typed limiter,the traditional LU-SGS and GMRES schemes are explored,and the improved LU-SGS scheme is developed by adding the error compensation of high order term.In addition,the transonic inviscid flow around NACA 0012 airfoil and RAE2822 airfoil as the examples are calculated.The numerical experiments indicate that the error compensation LU-SGS algorithm has the advantages of low storage requirements and easy programming,and the computational efficiency is close to GMRES algorithm and more than 2.5 times of LU-SGS one.

discontinuous Galerkin FEM;LU-SGS algorithm;GMRES algorithm;limiter

A

TP301

10.3778/j.issn.1002-8331.1311-0474

DUAN Zhijian.Efficient implicit algorithm s for discontinuous Galerkin finite elem ent method of Eu ler equation. Computer Engineering and Applications,2014,50（16）：21-24.

國家自然科學(xué)基金（No.11002117）；咸陽師范學(xué)院科研基金項目（No.09XSYK 204，No.09XSYK 209）。

段治?。?980—），男，博士，講師，研究領(lǐng)域?yàn)槔碚撆c計算流體力學(xué)、信息處理中的并行算法。E-mail：zhijian_duan@126.com

2013-12-01

2014-01-28

1002-8331（2014）16-0021-04

CNKI網(wǎng)絡(luò)優(yōu)先出版：2014-03-03,http://www.cnki.net/kcms/doi/10.3778/j.issn.1002-8331.1311-0474.htm l