高阶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-)，男，博士，副教授。主要研究方向：理论与计算流体力学。