高階DG格式多重網(wǎng)格方法構(gòu)造中的不同隱式策略影響研究

李廣佳，郝海兵，陳智，張強(qiáng)

(1.中國航天空氣動(dòng)力技術(shù)研究院 第十一總體設(shè)計(jì)部，北京　100074) (2.中國航空計(jì)算技術(shù)研究所 第七研究室，西安　710068) (3.中航飛機(jī)研發(fā)中心 動(dòng)力燃油所，西安　710089) (4.西北工業(yè)大學(xué) 航空學(xué)院，西安　710072)

?

高階DG格式多重網(wǎng)格方法構(gòu)造中的不同隱式策略影響研究

李廣佳1，郝海兵2，陳智3，張強(qiáng)4

(1.中國航天空氣動(dòng)力技術(shù)研究院 第十一總體設(shè)計(jì)部，北京100074) (2.中國航空計(jì)算技術(shù)研究所 第七研究室，西安710068) (3.中航飛機(jī)研發(fā)中心 動(dòng)力燃油所，西安710089) (4.西北工業(yè)大學(xué) 航空學(xué)院，西安710072)

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

間斷Galerkin方法；p型多重網(wǎng)格；歐拉方程；LU-SGS；GMRES

0　引　言

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

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

1　控制方程

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

(1)

對(duì)于氣體動(dòng)力學(xué)方程：

(2)

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

2　間斷Galerkin方法

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

定義有限元空間Vh：

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

(3)

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

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

(4)

式中：Ω為計(jì)算域；?Ω為Ω的邊界。

在每個(gè)單元內(nèi)：

(5)

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

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

(6)

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

為了便于計(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)單元，則式(6)變?yōu)?/p>

(7)

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

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

(8)

(9)

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

3　p型多重網(wǎng)格

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

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

(10)

(3) 在p0層上計(jì)算強(qiáng)迫項(xiàng)：

Fp0=Rp0-R(up0)

(11)

(4) 在p0層上計(jì)算一個(gè)時(shí)間步，計(jì)算殘值：

R=R(up0)+Fp0

(12)

(13)

從初始解開始，低階近似多項(xiàng)式上的解都會(huì)對(duì)高階近似多項(xiàng)式上的解進(jìn)行修正，從而加速收斂。

(14)

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

4　GMRES算法

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

將p0層上的Euler方程空間離散并進(jìn)行線化，得

(15)

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

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

AΔUn=Rn

(16)

GMRES算法的具體過程如下：

(2) 開始進(jìn)行內(nèi)迭代，通過Arnoldi方法構(gòu)造Krylov子空間的標(biāo)準(zhǔn)規(guī)范正交基：

①對(duì)于j=1,2,…,m，

b.對(duì)于i=1,2,…,j，計(jì)算

(17)

②Krylov子空間的標(biāo)準(zhǔn)正交基可以寫成矩陣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) 重開始

根據(jù)數(shù)學(xué)上收斂條件‖r(m)‖2來判定是否終止迭代計(jì)算。若‖r(m)‖2適當(dāng)小則終止迭代，否則

(26)

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

5　算例與分析

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

NACA0012翼型的流動(dòng)計(jì)算狀態(tài)為：Ma∞=0.8，攻角α=1.25°。計(jì)算的整個(gè)流場(chǎng)包含3 531個(gè)網(wǎng)格節(jié)點(diǎn)，6 789個(gè)網(wǎng)格單元，計(jì)算網(wǎng)格如圖1所示。

圖1　NACA0012翼型計(jì)算網(wǎng)格Fig.1　Computational grid of NACA0012 airfoil

關(guān)于CPU時(shí)間和迭代步數(shù)的密度最大殘值收斂歷程對(duì)比如圖2所示。

(a) CPU時(shí)間

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

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

ONERAM6機(jī)翼的流動(dòng)計(jì)算狀態(tài)為：Ma∞=0.84，攻角α=3.06°。計(jì)算的整個(gè)流場(chǎng)包含21 019個(gè)網(wǎng)格節(jié)點(diǎn)和99 350個(gè)網(wǎng)格單元。機(jī)翼表面網(wǎng)格如圖3所示。

圖3　M6機(jī)翼表面網(wǎng)格Fig.3　Surface mesh of M6 wing

采用p型多重網(wǎng)格方法計(jì)算后的機(jī)翼上表面壓力等值線如圖4所示，可以看出，λ狀激波結(jié)構(gòu)非常清晰，外弦和內(nèi)弦激波大約在展長(zhǎng)87%處相交，在94%處又分開，和實(shí)驗(yàn)結(jié)果[14]基本吻合。

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

關(guān)于CPU時(shí)間和迭代步數(shù)的密度最大殘值收斂歷程對(duì)比如圖5所示。

(a) CPU時(shí)間

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

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

6　結(jié)　論

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

(2) DG(p0)層上采用不同的隱式迭代方法對(duì)p型多重網(wǎng)格方法整體加速性能影響比較明顯；與顯式TVD-RKDG時(shí)間格式相比，LU-SGS隱式迭代方法和GMRES隱式迭代方法均可顯著提高收斂速度，其中GMRES隱式迭代方法的性能最佳，計(jì)算效率提高約10倍；LU-SGS隱式迭代方法次之，計(jì)算效率提高約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] 郝海兵, 楊永. 非結(jié)構(gòu)網(wǎng)格上p型多重網(wǎng)格法流場(chǎng)數(shù)值模擬[J]. 計(jì)算力學(xué)學(xué)報(bào), 2011, 28(3): 360-365.

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

[11] 郝海兵, 楊永, 左歲寒. 高階間斷Galerkin方法求解三維歐拉方程的研究[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 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”計(jì)算項(xiàng)目(2014CB744803)

郝海兵,hao_hb@163.com

1674-8190(2016)03-279-07

V211

A

10.16615/j.cnki.1674-8190.2016.03.003

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

郝海兵(1981-)，男，博士，高級(jí)工程師。主要研究方向：理論與計(jì)算流體力學(xué)、氣動(dòng)優(yōu)化設(shè)計(jì)。

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

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