高階DG格式多重網格方法構造中的不同隱式策略影響研究

李廣佳，郝海兵，陳智，張強

(1.中國航天空氣動力技術研究院 第十一總體設計部，北京　100074) (2.中國航空計算技術研究所 第七研究室，西安　710068) (3.中航飛機研發中心 動力燃油所，西安　710089) (4.西北工業大學 航空學院，西安　710072)

?

高階DG格式多重網格方法構造中的不同隱式策略影響研究

李廣佳1，郝海兵2，陳智3，張強4

(1.中國航天空氣動力技術研究院 第十一總體設計部，北京100074) (2.中國航空計算技術研究所 第七研究室，西安710068) (3.中航飛機研發中心 動力燃油所，西安710089) (4.西北工業大學 航空學院，西安710072)

DG方法是一種非常具有潛力的高精度方法，但其在對復雜外形的數值模擬方面仍存在內存需求量大、計算量巨大等不足。為了進一步提高DG方法求解Euler方程的效率，在傳統p型多重網格的基礎上，結合LU-SGS和GMRES兩種隱式迭代方法，研究其整體加速性能。p型多重網格方法通過對不同階次多項式近似解進行遞歸迭代求解，來達到加速收斂的目的。高階近似 (p>0)使用顯式龍格庫塔格式，最低階近似 (p=0)使用隱式格式。對NACA0012翼型和ONERA M6機翼跨音速無粘流動進行數值模擬，結果表明：與顯式TVD-RKDG時間格式相比，DG(p0)層上采用LU-SGS和GMRES的p型多重網格方法收斂速度均得到明顯提高，且GMRES迭代法性能最佳，LU-SGS迭代法次之。

間斷Galerkin方法；p型多重網格；歐拉方程；LU-SGS；GMRES

0　引　言

間斷Galerkin方法(DG方法)由于集合了有限元法(FEM)和有限體積法(FVM)的優點，近年來得到了廣泛關注。該方法最早由W.H.Reed等[1]于1973年在求解中子輸運方程時引入，但在隨后的很長一段時間內并未被人們所重視，直到20世紀90年代，B.Cockburn等[2-3]提出了Runge-Kutta DG方法(TVD-RKDG)，并給出了部分收斂性和穩定性證明，該方法才被廣泛地應用于計算流體力學。DG方法只需在單元內部通過提高逼近多項式的階次就能實現高階精度，具有很好的緊致特性，且易于實行并行計算。此外，由于解在單元邊界上保持間斷，該方法非常適用于求解間斷問題。目前，DG方法作為一種最具潛力的高精度方法[4]，將會對計算流體力學的發展產生積極的推動作用。

本文在文獻[10]的基礎上，對DG(p0)層上的不同隱式策略在p型多重網格方法中的應用展開研究，并考察其整體加速性能。其中，DG(p0)層隱式求解主要選取FVM中的兩種常用的隱式迭代方法，即LU-SGS和GMRES。

1　控制方程

非定常Euler方程在直角坐標系下的守恒形式為

(1)

對于氣體動力學方程：

(2)

式中：γ為比熱比，取γ=1.4。

2　間斷Galerkin方法

利用間斷Galerkin方法求解式(1)，首先需要將計算區域劃分成互不重疊的子域。子域可以選取為任意形狀，對于二維空間，本文采用三角形非結構網格；對于三維空間，采用四面體非結構網格。

定義有限元空間Vh：

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

(3)

式中：τh為子域空間；V(K)為局部函數空間，取作p(p=1,2,L)次多項式的集合。

假設間斷函數在有限元空間中的近似解為Uh，將式(1)兩邊同乘以試驗函數v，寫成變分形式，再采用分部積分，得到守恒方程組的弱形式：

(4)

式中：Ω為計算域；?Ω為Ω的邊界。

在每個單元內：

(5)

式中：φj(x)為基函數。

將式(5)帶入式(4)，得到半離散形式的守恒方程：

(6)

式(6)稱為p階間斷Galerkin有限元離散，p為式(5)中基函數的最大階數。可以看出Uh、vh在整個計算區域內不再連續，在單元邊界上存在間斷，且整個流場的自由度已經轉變成求解插值系數u，而非流場守恒變量U。

為了便于計算，通常還需要進行坐標變換。針對三角形單元，本文采用面積坐標；針對四面體單元，則采用體積坐標。將總體笛卡爾坐標轉換成局部自然坐標，將空間中的任意單元轉變成局部坐標系中的標準單元，則式(6)變為

(7)

式中：|J|為體坐標變換雅可比矩陣；|Js|為面坐標變換雅可比矩陣。

對于式(7)中的積分，采用高斯數值積分：

(8)

(9)

其中數值通量選取HLLC格式。HLLC近似黎曼解方法是一種具有實際應用價值的、能夠精確模擬接觸間斷和剪力波的最簡單的平均狀態方法。采用HLLC近似黎曼解方法得到的結果與采用精確黎曼解方法得到的結果基本相同。此外，DG方法在計算過程中會產生數值振蕩，嚴重時可能導致無法求解，因此，本文采用間斷探測器和斜率限制器[11]相結合的方法來抑制間斷處的數值振蕩。

3　p型多重網格

目前，基于非結構網格的有限體積法在求解Euler方程和N-S方程時，廣泛使用h型多重網格方法進行加速收斂。h型多重網格方法通過粗網格上的解修正細網格上的解以達到加速收斂的目的。p型多重網格方法可以看作是h型多重網格方法在高階有限元空間的一種推廣，即方程組的解在不同階次多項式近似層上進行迭代求解。通過將高階近似上的低頻誤差轉換成低階近似上的高頻誤差，達到快速消除高階近似上的低頻誤差的目的。針對DG(p1)方法，本文使用兩層V循環，主要包括以下步驟：

(2) 將p1層上的解和殘值限制到低一層近似解p0上：

(10)

(3) 在p0層上計算強迫項：

Fp0=Rp0-R(up0)

(11)

(4) 在p0層上計算一個時間步，計算殘值：

R=R(up0)+Fp0

(12)

(13)

從初始解開始，低階近似多項式上的解都會對高階近似多項式上的解進行修正，從而加速收斂。

(14)

不同階次層之間的限制和延拓算子均可在母單元上計算。當選取正交基函數時，式(14)比較容易計算。

4　GMRES算法

在最低層p0近似上，DG方法即為一階精度中心格式的有限體積法，其具有相對成熟的隱式算法，例如ADI方法、點隱式方法、LU-SGS方法、GMRES方法等。本文選取LU-SGS[12]和GMRES[13]方法，其中LU-SGS方法詳見文獻[10]，限于篇幅，本文不再贅述；GMRES以Galerkin原理為基礎，建立Krylov子空間的一組標準正交基，再在Krylov子空間中利用最小二乘方法得到解向量。

將p0層上的Euler方程空間離散并進行線化，得

(15)

式中：V為單元體積；R為p0層上解的殘值。

將網格面上的無粘通量通過雅可比矩陣最大特征值分裂，得到整個流場空間的一個N維線性方程組(N為總網格數)：

AΔUn=Rn

(16)

GMRES算法的具體過程如下：

(2) 開始進行內迭代，通過Arnoldi方法構造Krylov子空間的標準規范正交基：

①對于j=1,2,…,m，

b.對于i=1,2,…,j，計算

(17)

②Krylov子空間的標準正交基可以寫成矩陣Vm：

(18)

由Gram-Schmidt正交化過程可知，式(18)可表示為如下矩陣形式：

(19)

式中：Hm∈Rm,m為上Hessenberg矩陣，

(20)

滿足Galerkin條件：

(21)

式(21)滿足最小二乘解條件：

(22)

若令z(m)=Vmy(m)，則式(22)右端極小化泛函可表示為

J(y)=‖‖r(0)‖2v(1)-AVmy‖2

(23)

利用式(19)可得

(24)

由此式(24)的解為

(25)

式中：y(m)為極小化式(24)的J(y)。

(4) 重開始

根據數學上收斂條件‖r(m)‖2來判定是否終止迭代計算。若‖r(m)‖2適當小則終止迭代，否則

(26)

回到步驟(2)進行重新迭代，直至滿足步驟(4)中的收斂條件。

5　算例與分析

為了研究p型多重網格方法中，DG(p0)層上分別采用LU-SGS、GMRES不同隱式迭代方法的整體加速性能，分別對繞NACA0012翼型和ONERAM6機翼的跨音速無粘流動進行數值模擬，并和TVD-RKDG方法進行比較，來驗證DG(p0)層上的不同隱式策略的p型多重網格方法加速收斂效果和精度。算例中使用的計算機基本配置CPU為酷睿I5 3.2GHz、內存為4G，所有網格均采用Delaunay方法生成。

NACA0012翼型的流動計算狀態為：Ma∞=0.8，攻角α=1.25°。計算的整個流場包含3 531個網格節點，6 789個網格單元，計算網格如圖1所示。

圖1　NACA0012翼型計算網格Fig.1　Computational grid of NACA0012 airfoil

關于CPU時間和迭代步數的密度最大殘值收斂歷程對比如圖2所示。

(a) CPU時間

(b) 迭代步數 圖2　NACA0012翼型殘值收斂歷程比較Fig.2　Comparison of convergence historyversus of NACA0012 airfoil

從圖2可以看出：當采用TVD-RKDG方法時，計算耗時約70min；當采用LU-SGS隱式迭代方法時，計算耗時約12min，計算效率提高6倍左右；而采用GMRES隱式迭代方法時，計算耗時約6min，計算效率提高12倍左右。

ONERAM6機翼的流動計算狀態為：Ma∞=0.84，攻角α=3.06°。計算的整個流場包含21 019個網格節點和99 350個網格單元。機翼表面網格如圖3所示。

圖3　M6機翼表面網格Fig.3　Surface mesh of M6 wing

采用p型多重網格方法計算后的機翼上表面壓力等值線如圖4所示，可以看出，λ狀激波結構非常清晰，外弦和內弦激波大約在展長87%處相交，在94%處又分開，和實驗結果[14]基本吻合。

圖4　M6機翼上表面壓力等值線Fig.4　Pressure isolines on upper surface of M6 wing

關于CPU時間和迭代步數的密度最大殘值收斂歷程對比如圖5所示。

(a) CPU時間

(b) 迭代步數 圖5　M6機翼殘值收斂歷程比較Fig.5　Comparison of convergence historyversus of M6 wing

從圖5可以看出：當采用TVD-RKDG方法時，計算耗時約900min；當采用LU-SGS隱式迭代方法時，計算耗時約200min，計算效率提高4.5倍左右；而采用GMRES隱式迭代方法時，計算耗時約100min，計算效率提高9倍左右。

6　結　論

(1) 本文研究了DG(p0)層上采用不同的隱式迭代方法對p型多重網格方法整體加速性能的影響，其中隱式迭代方法選擇了在FVM中被廣泛使用的LU-SGS和GMRES方法。通過對NACA0012翼型和M6機翼跨音速無粘流動的數值模擬，驗證了本文方法的加速收斂效果。

(2) DG(p0)層上采用不同的隱式迭代方法對p型多重網格方法整體加速性能影響比較明顯；與顯式TVD-RKDG時間格式相比，LU-SGS隱式迭代方法和GMRES隱式迭代方法均可顯著提高收斂速度，其中GMRES隱式迭代方法的性能最佳，計算效率提高約10倍；LU-SGS隱式迭代方法次之，計算效率提高約5倍。

[1]ReedWH,HillTR.TrianglemeshmethodsfortheNeutrontransportequation[R].LA-UR-73-479, 1973.

[2]CockburnB,HouS,ShuCW.TheRunge-KuttalocalprojectiondiscontinuousGalerkinfiniteelementmethodforconservationlaws.IV.themultidimensionalcase[J].MathematicsofComputation, 1990, 54: 545-581.

[3]CockburnB,ShuCW.Runge-KuttadiscontinuousGalerkinmethodsforconvection-dominatedproblems[J].JournalofScientificComputing, 2001, 16(3): 173-261.

[4]WangZJ.High-ordermethodsfortheEulerandNavier-Stokesequationsonunstructuredgrids[J].ProgressinAerospaceSciences, 2007, 43(1-3): 1-41.

[5]RonquistEM,PateraAT.SpectralelementmultigridI:formulationandnumericalresults[J].JournalofScientificComputing, 1987, 2(4): 389-406.

[6]BassiF,GhidoniA,RebayS,etal.High-orderaccuratep-multigriddiscontinuousGalerkinsolutionoftheEulerequations[J].InternationalJournalforNumericalMethodsinFluids, 2009, 60(8): 847-865.

[7]FidkowskiKJ,OliverTA,LuJ,etal. p-multigridsolutionofhigh-orderdiscontinuousGalerkindiscretizationsofthecompressibleNavier-Stokesequations[J].JournalofComputationalPhysics, 2005, 207(1): 92-113.

[8]LuoH,BaumJD,L?hnerR.Ap-multigriddiscontinuousGalerkinmethodfortheEulerequationsonunstructuredgrids[J].JournalofComputationalPhysics, 2006, 211(2): 767-783.

[9]LuoH,BaumJD,L?hnerR.Fastp-multigriddiscontinuousGalerkinmethodforcompressibleflowsatallspeeds[C].AIAAJournal, 2008, 46(3): 635-652.

[10] 郝海兵, 楊永. 非結構網格上p型多重網格法流場數值模擬[J]. 計算力學學報, 2011, 28(3): 360-365.

HaoHaibing,YangYong.Theresearchofp-multigridsolutionfordiscontinuousGalerkinmethod[J].ChineseJournalofComputationalMechanics, 2011,28(3): 360-365.(inChinese)

[11] 郝海兵, 楊永, 左歲寒. 高階間斷Galerkin方法求解三維歐拉方程的研究[J]. 西北工業大學學報, 2011, 29(1): 128-132.

HaoHaibing,YangYong,ZuoSuihan.Effectivelyapplyinghigh-orderdiscontinuousGalerkinmethod(DGM)tosolving3-DEulerequationsonunstructuredgrids[J].JournalofNorthwesternPolytechnicalUniversity, 2011, 29(1): 128-132.(inChinese)

[12]JamesonA,YoonS.Lower-upperimplicitschemeswithmultiplesgridsfortheEulerequations[J].AIAAJournal, 1987, 25(7): 929-935.

[13]SaadY,SchultzMH.GMRES:ageneralizedminimalresidualalgorithmforsolvingnon-symmetriclinearsystems[J].SIAMJournalonScientificandStatisticalComputing, 1986, 7(3): 856-869.

[14]BarthTJ.A3-DupwindEulersolverforunstructuredmeshes[C].AIAA-91-1548-CP, 1991.

(編輯：馬文靜)

Study of Different Implicit Schemes for High Order DG Multigrid Method

Li Guangjia1， Hao Haibing2, Chen Zhi3, Zhang Qiang4

(1.The Eleventh Design Department, China Academy of Aerospace Aerodynamics, Beijing 100074, China) (2.No.7 Department, Aeronautical Computing Technique Research Institute, Xi’an 710068, China) (3.Institute of Power and Fuel, AVIC Aircraft Research and Development Center, Xi’an 710089, China) (4.School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China)

As a kind of high order method, the DG method has great potential in the numerical simulation of complex configuration, but there are a large amount of internal memory requirement and a huge amount of computation. In order to further improve the efficiency of the method for solving the Euler equation, the overall acceleration performance is studied by combining the two implicit iterative methods of LU-SGS and GMRES. The convergence acceleration ofp-multigrid method is achieved through the recursive iterative solving of different order polynomial approximations. In order to achieve better convergence effect, implicit scheme is implemented on the lowest-order approximation while explicit schemes are implemented on the higher-order approximations. The transonic inviscid flow around NACA0012 airfoil and ONERA M6 wing are simulated. The numerical results show that, compared with the explicit TVD-RKDG method, convergence rates for thep-type multigrid methods, which adopt LU-SGS or GMRES scheme onDG(p0) layer, are significantly improved, while the performance of GMRES iterative method is best and LU-SGS iteration is secondary.

discontinuous Galerkin methods;p-type multigrid; Euler equations; LU-SGS; GMRES

2016-04-20；

2016-06-06

國家“973”計算項目(2014CB744803)

郝海兵,hao_hb@163.com

1674-8190(2016)03-279-07

V211

A

10.16615/j.cnki.1674-8190.2016.03.003

李廣佳(1979-)，男，碩士，高級工程師。主要研究方向：計算流體力學、飛行器設計。

郝海兵(1981-)，男，博士，高級工程師。主要研究方向：理論與計算流體力學、氣動優化設計。

陳智(1986-)，女，工程師。主要研究方向：計算流體力學。

張強(1979-)，男，博士，副教授。主要研究方向：理論與計算流體力學。