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

段治健

1.咸陽師范學院數學與信息科學學院，陜西咸陽 712000

2.西北工業大學翼型葉柵空氣動力學國防科技重點實驗室，西安 710072

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

段治健1，2

1.咸陽師范學院數學與信息科學學院，陜西咸陽 712000

2.西北工業大學翼型葉柵空氣動力學國防科技重點實驗室，西安 710072

針對Euler方程，設計了適合間斷Galerkin有限元方法的LU-SGS、GMRES以及修正LU-SGS隱式算法。采用Roe通量以及Van A lbada限制器技術實現了經典LU-SGS、GMRES算法，引入高階項誤差補償，發展了修正LU-SGS算法。以NACA 0012、RAE2822翼型為例驗證分析了算法的可靠性和高效性。結果表明修正LU-SGS算法存儲量較少，程序實現方便，而且計算效率是LU-SGS算法的2.5倍以上，接近于循環GMRES算法。

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

1 引言

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

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

2 間斷有限元方法

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

其中，U為守恒變量，Fj為對應無粘通量，ρ，e，p分別為氣體的密度、單位體積總內能及壓強。ui是xi方向的速度，δij是K ronecker函數。考慮到理想氣體的熱力學特性，狀態方程如下：

取γ=1.4為比熱比。

在式（1）兩邊乘以測試函數W，在計算域內運用分部積分公式，得到弱形式：

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

在計算域劃分網格上，用近似解Uh，Wh代替解析解U，W，每個單元內：

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

3 GMRES方法

具有平方斂速的共軛梯度法[10]是New ton迭代隱式算法的一種，GMRES算法屬于其中一種，因而倍受關注，循環GMRES算法相對于GMRES算法收斂性通常差一些，但是存儲量相對較少。因此本文采用循環GMRES算法。對于給定的非對稱線性方程組Ax=f，循環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方法

當式（4）中取為常數時，則方程變為：

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

從而得到一個N維的線性方程組（N為網格總數）：

將系數矩陣分解為A=D+L+U，式（7）變為：

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

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

5 數值算例及結果分析

算例1 NACA 0012翼型跨音速無粘繞流，計算網格節點2 270個，單元數為4 355個。計算狀態為Ma=0.8，α=1.25°。

算例2 RAE2822翼型跨音速無粘繞流，計算網格節點2 270個，單元數為4 355個。計算狀態為Ma=0.725，α=2.54°。

圖1和圖2分別為NACA 0012和RAE2822翼型的壓力系數分布曲線，數值計算的結果與實驗值基本吻合。

圖1 NACA 0012翼型表面壓力系數對比

圖2 RAE2822翼型表面壓力系數對比

圖3和圖4分別為兩個算例時間收斂曲線。可以看出，循環GMRES算法計算效率遠高于LU-SGS格式。計算效率是LU-SGS格式的3倍左右。修正LU-SGS格式效率接近于循環GMRES算法，且所需存儲量較少，單步迭代時間短，計算效率是LU-SGS格式的2.5倍以上。

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

圖4 RAE2822翼型計算時間曲線

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

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

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

6 結束語

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

[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]李勁杰，楊青，楊永年.三維非結構網格Euler方程的LU-SGS算法及其改進[J].計算物理，2006，23（6）：748-752.

[8]李春娜，葉正寅.基于二維非結構網格的GMRES隱式算法[J].西北工業大學學報，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].計算機科學，2010，37（3）：242-244.

[12]曹芳芳，呂全義.解非對稱塊三對角線性方程組的并行算法[J].西北工業大學學報，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.

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

段治健（1980—），男，博士，講師，研究領域為理論與計算流體力學、信息處理中的并行算法。E-mail：zhijian_duan@126.com

2013-12-01

2014-01-28

1002-8331（2014）16-0021-04

CNKI網絡優先出版：2014-03-03,http://www.cnki.net/kcms/doi/10.3778/j.issn.1002-8331.1311-0474.htm l