一种适用于非结构网格的间断Galerkin有限元LU-SGS隐式方法

马明生, 龚小权, 邓有奇, 赵 辉

1.西北工业大学 航空学院, 陕西 西安 710072 2.中国空气动力研究与发展中心计算空气动力研究所, 四川 绵阳 621000



一种适用于非结构网格的间断Galerkin有限元LU-SGS隐式方法

马明生1,2, 龚小权1, 邓有奇1, 赵 辉1

1.西北工业大学 航空学院, 陕西 西安 710072 2.中国空气动力研究与发展中心计算空气动力研究所, 四川 绵阳 621000

具有TVD性质的显式Runge-Kutta间断Galerkin(RKDG)格式在CFD领域得到广泛应用,但是显式计算稳定性差、计算效率低。为改善时间推进效率,基于高阶间断Galerkin有限元方法,采用欧拉一阶后差(BDF1),发展了一套高效的隐式LU-SGS(lower upper-symmetric Gauss-Seidel)求解方法,方法基于MPI并行实现,适合于不同计算精度。针对非线性系统左端项矩阵,对比了简化前后LU-SGS的计算效率。建立的间断Galerkin有限元方法基于非结构网格,采用Taylor基函数,计算精度最高达到四阶精度。通过NACA0012翼型以及M6机翼算例对发展的LU-SGS方法进行了考察,与显式算法相比,隐式格式的迭代步数和CPU时间均较大程度减小,效率能够提高1个量级以上。最后将隐式算法用于复杂外形翼身组合体F4的流场计算,结果表明所发展的隐式方法具有较好的鲁棒性,能够用于复杂外形计算。

间断Galerkin有限元;欧拉方程;Taylor基函数;LU-SGS;计算效率;非结构网格

随着CFD和计算机技术的发展,CFD工作者希望得到更精细的流场结构以及更准确的气动数据,比如阻力、力矩、热流。精细化的网格或者高精度的计算格式是实现这些要求的途径之一。在规则网格上进行光滑函数的数值近似,n阶精度的计算格式,计算误差与hn(h为网格尺寸)成正比。高精度格式能够给出一些复杂流动现象的精细流场结构,相比传统二阶格式具有更高的分辨率。而非结构网格相比结构网格具有更大的灵活性,更适应复杂外形和边界条件。因此发展适用于非结构网格的高精度格式一直是计算流体力学(CFD)研究的热点。

间断Galerkin有限元方法(discontinuous galerkin finite element method,DGM)是一种适合于结构网格以及非结构网格的高精度格式。它最早出现在1973年Reed和Hill[1]求解中子输运方程问题的论文中。DGM结合了有限体积法和有限元法的优点。作为一种有限元方法,它采用有限元法中基函数思想来构造单元内的变量分布,容易实现高精度。同时它又具有有限体积法积分形式的计算形式且融入了有限体积法中数值通量、Riemann间断、限制器等思想。DGM具有提高计算精度容易,高精度时计算模板紧致,适合非结构网格,适合处理复杂外形及复杂边界、初值问题,并行能力强,自适应计算方便等优点。间断Galerkin有限元方法[2]在求解守恒方程组问题中得到广泛应用。经过Cockburn和Shu等的持续研究,一种具有TVD性质的显式Runge-Kutta间断Galerkin(RKDG)格式得以逐步完善。

RKDG格式单步计算时间少,内存需求量相对较小,程序容易实现,但是随着计算精度的提高, RKDG格式对应的稳定性条件将越来越严格,时间推进步长受到明显限制,计算效率低下,特别是黏性计算。一些加速技术如当地时间步长、残值光顺、多重网格等,在一定程度上加速了收敛过程。即便如此,最大的时间步长依然受到当地稳定性条件的限制。与显式的时间格式相比,隐式方法基本不受时间推进步长的限制,因此发展高精度间断Galerkin有限元方法的时间隐式离散方法,有着重要的工程应用价值。

目前间断Galerkin有限元方法已有多种高效的隐式离散方法,比如带预处理的GMRES[3]、二阶的Crank-Nicholson(CN2)格式[4]、四阶隐式的Runge-Kutta(IRK4)等。LU-SGS隐式方法因为其在计算效率、鲁棒性及内存占用方面的优势,在有限体积方法(FVM)中被大量采用。Yoon与Jameson在1988年[5]首先将其作为迭代方法用于方程求解,随后Rieger和Jameson将其用于求解三维黏性流场,从此各种改进的LU-SGS隐式方法被用来求解基于结构及非结构网格的流场。详细介绍LU-SGS在高精度间断Galerkin有限元方法中应用的文献较少,Haga和Wang将该方法用于基于非结构四面体网格的谱方法[6]。郝海兵、张强等在DGM四面体网格上实现了p=1阶的LU-SGS迭代[7]。郭永恒在其博士论文中[8]研究了不同LU-SGS改进方法对收敛的影响。刘伟在文献[9]中发展了基于Newton/Gauss-Seidel迭代的DGM隐式方法。

本文基于建立在非结构网格上的高精度间断Galerkin有限元方法,发展了LU-SGS隐式计算方法,比较了不同计算精度下的收敛曲线,通过流场熵增云图给出了不同精度下的计算结果。同时对左端项进行了简化,给出了左端项不同简化方法对收敛性的影响。给出了一系列二维、三维算例不同精度下的收敛曲线,比较了二阶精度下有限体积和DGM下的收敛速度。

1 数值方法

1.1 控制方程

守恒形式的无黏可压缩欧拉方程组可写为

(1)

其中守恒变量、无黏对流通量为

(2)

ρ、u、v、w、p、E、H分别代表密度、3个方向速度、压力、总能和总焓。

(3)

γ为比热比。

1.2 DGM空间离散

对方程组(1)采用间断Galerkin有限元离散。将其两端同时乘以测试函数φi,并对整个求解域积分Ω,再采用分部积分得到如下公式

(4)

式中，∂Ω为求解域Ω的边界,n为边界的法向量。将求解域Ω剖分为互不重叠的非结构网格单元,并对每个单元采用上面的弱形式。每个计算单元的每个守恒变量在单元内的函数分布为

(5)

uj为单元守恒变量的自由度,即为每时间步需要求解的量。N=N(p,d)由插值多项式阶数p和求解的空间维数d决定。φj是阶数为p的基函数。由于Taylor基函数[10]对所有单元形式都一样且为等级基,对非结构混合网格较为合适,因此本文采用Taylor基函数。将方程(5)带入方程(4)得

(6)

最后整个求解域方程可写为如下形式

(7)

Res为无黏通量残差项,由面积分和体积分构成,二者均采用高斯数值积分计算。M为全局的质量矩阵,由各单元的质量矩阵构成。无黏通量残差项具体计算为

(8)

wv(l)、|Jl|、φil分别为体积分对应的单元高斯积分点的权系数、雅克比矩阵行列式值以及第i个基函数的梯度。ws(k)、|Jsk|、φis分别为单元面面积分对应的高斯积分点k的权系数、雅克比矩阵行列式值以及第i个基函数的值。Wak、Wbk分别对应边界面上当前单元和相邻单元在积分点k的变量值,Wal为当前单元体积分点l的变量值。

1.3 隐式LU-SGS方法

对方程(7)第一式采用一阶欧拉后差得到

(9)

方程(9)为一大型的线性系统Ax=B,其中矩阵A为一大型稀疏块矩阵。以三维p=3(四阶精度)DGM为例,矩阵A的维数为n-cell×(5×20),n-cell为计算单元总数。如果直接存储该矩阵并求逆,其存储量和计算量太大,目前大多采用迭代法来求解该线性系统。

将面积分和体积分带入左端项的Resn。

(10)

将左端项矩阵写为

(11)

式中，D为对角块矩阵,L为严格下三角矩阵,U为严格上三角矩阵。

为方便推导,进一步得到矩阵D、L、U,本文只写出左端项矩阵中与当前计算单元有关的块,将当前单元编号用a表示,相邻单元用b表示,

对(10)式采用链式法则

(12)

(13)

进一步得到面积分通量对守恒变量的偏导数

(14)

将方程(8)的第一项以及方程(14)带入方程(12)得

(15)

(16)

(17)

2 数值模拟结果及分析

本文首先通过等熵涡算例验证了DGM的计算精度。为验证本文发展的LU-SGS方法的收敛性,计算了几类典型的测试算例。测试过程中,右端对流通量项采用Roe′s格式,并采用三阶TVD-Runge-Kutta方法来对比分析。利用NACA0012翼型亚声速算例考察了不同精度的LU-SGS收敛曲线,研究了CFL数对收敛历程的影响,给出了左端项简化前后的收敛曲线对比。计算了三维M6机翼亚声速绕流,给出了不同精度的收敛曲线。最后通过F4翼身组合体给出了发展的LU-SGS隐式方法在复杂外形中的应用。在分析过程中,本文采用成熟且经过验证的有限体积方法程序作为对比分析的参考。

2.1 等熵涡问题

等熵涡问题有精确解,用来对二维Euler问题格式的精度进行验证,这里把问题扩展到三维问题。等熵涡问题描述的是在初始的平均流(ρ,u,v,w,p)=(1,1,1,0,1)上加入一个等熵涡,涡的大小按照(18)式给出。

(18)

式中，T=p/ρ,熵S=p/ρr,ε为涡的强度ε=5.0;计算区域取为[0,10]×[0,10]×[0,10],边界条件按照精确解给定,时间上采用显式三阶Runge-Kutta法推进至2 s。计算采用的网格为N×N×N(N=10,20,40)六面体网格。表1给出了不同精度下的L1 Error和L2 Error,可以看到方法达到了设计的精度。

表1 等熵涡问题精度验证(六面体网格)

2.2 NACA0012翼型亚声速绕流

为考察DGM的精度以及发展的隐式算法的收敛性,本文计算了NACA0012翼型在来流马赫数Ma=0.5时亚声速绕流。计算网格如图1所示。

图1 NACA0012计算网格

从图3a)的曲线1和5,CFL=1,简化左端项矩阵对收敛速度有一定影响,导致迭代步数略有增加。从图3b)看到,对于CFL=1的情况,由于简化左端项矩阵后,单步计算时间减小,总的收敛时间减小。由于简化左端项矩阵后CFL数不能取大,因此未简化左端项矩阵当CFL数取大时,其收敛速度快于简化后的状态。

图2 DGM不同计算精度及二阶FVM收敛历程 图3 DGM二阶精度不同CFL数收敛历程

2.3 ONERA M6机翼亚声速绕流

为进一步考察LU-SGS方法在三维外形计算效率,本文计算了ONERA-M6亚声速流场,计算马赫数Ma=0.5,迎角0°。网格如图4所示,网格共约14万个四面体单元,机翼前后缘及空间分别采用各向异性和各向同性四面体网格。计算采用64个CPU。

图4 M6机翼计算网格

图5给出了DGM不同计算精度下LU-SGS残差收敛曲线,同时给出了FVM以及RKDG显式结果。隐式计算的CFL=1 000,显式计算的CFL=0.3。从残差随迭代步数图看,发展的LU-SGS隐式方法收敛效率远高于显式方法。二阶精度下,由于DGM计算的自由度为FVM的4倍且计算面通量和体积分通量时采用更多的高斯积分点,二阶DGM的单步计算时间远大于FVM。

图5 不同计算精度DGM及二阶FVM收敛曲线

2.4 DLR-F4翼身组合体亚声速绕流

为验证发展的隐式算法的初步工程计算能力,数值模拟了DLR-F4翼身组合体亚声速绕流。计算网格如图6所示。

图6 F4表面网格

四面体网格单元总数为637 126个,计算马赫数Ma=0.5,来流迎角为0°,计算采用128个CPU,图7给出了残差收敛曲线,相同CFL数以及相同的dt计算方法下,相对于有线体积方法,DGM残差收敛速度较快,但是下降量级较少。这是由于在迭代过程中DGM一直修改迭代得到的自由度,造成限制单元残差不降,因此需要进一步研究该方法的限制器。

图7 不同方法及精度的残差收敛曲线

3 结 论

本文发展了一种高效的LU-SGS隐式计算方法,方法适用于不同计算精度,通过一系列的算例表明:

1) 本文发展的隐式方法CFL数能够取很大。相比显式RKDG,其计算效率提高了1～2个量级。

2) 发展的隐式算法对左端项残差比较敏感,需小心简化左端项残差,简化可能限制CFL数,导致整个计算时间增加,计算效率降低。

3) 隐式计算时,单元左端项矩阵包含时间项、面积分项和体积分项,单元对应的矩阵D的矩阵元素几乎都不为零,降低了CFL数对收敛性的影响。

4) 发展的隐式方法能够初步应用于复杂外形计算。

[1] Reed W H, Hill T R. Triangular Mesh Methods for the Neu-Tron Transport Equation[R]. Technical Report LA-UR-73-479, Los Alamos Scientific Lab, 1973

[2] Cockburn B, Shu Chiwang. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws Ⅳ: The Multidimensional Case[J]. Mathematics of Computation, 1990, 54: 545-581

[3] Persson Per-Olof, Peraire Jaime. An Efficient Low Memory Implicit DG Algorithm for Time Dependent Problems[C]∥44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2006: 9-12

[4] Wang L, Mavriplis D J. Implicit Solution of the Unsteady Euler Equation for High-Order Accurate Discontinuous Galerkin Discretizations[J]. J Comput Phys, 2007, 225:1994-2005

[5] Yoon S, Jameson A. Lower-Upper Implicit Schemes with Multiple Grids for the Euler Equatins[R]. AIAA-1986-0105

[6] Hagal Takanori, Sawada Keisuke, Wang Z J. An Implicit LU-SGS Scheme for the Spectral Volume Method on Unstructured Tetrahedral Grids[J]. Commun Comput Phys, 2009,6(5):978-996

[7] 郝海兵,张强,杨永. 基于LU-SGS迭代的DGM隐式算法研究[J]. 西北工业大学学报,2014,32(3):346-350

Hao Haibing, Zhang Qiang, Yang Yong. An Implicit Scheme of Discontinuous Galerkin Method(DGM) Based on LU-SGS Iterative Method[J]. Journal of Northwestern Polytechnical University, 2014,32(3):346-350 (in Chinese)

[8] 郭永恒. 双曲守恒律方程组的间断Galerkin方法研究[D]. 西安:西北工业大学,2013

Guo Yongheng. Research of Discontinuous Galerkin Method for Hyperbolic Conservation Equations[D]. Xi′an, Northwestern Polytechnical University,2013 (in Chinese)

[9] 刘伟,张来平,赫新,等. 基于Newton/Gauss-Seidel迭代的 DGM隐式方法[J]. 力学学报,2012,44(4):792-796

Liu Wei, Zhang Laiping, He Xin, et al. An Implicit Algorithm for Discontinuous Galerkin Method Based on Newton/Gauss-Seidel Iterations[J]. Acta Mechanica Sinica, 2012,44(4): 792-796 (in Chinese)

[10] Hong L, Baum J D, Lohner R. A Discountinuous Galerkin Method Based on a Taylor Basis for The Compressible Flows on Arbitrary Grids[J]. J Comput Phys, 2008, 227: 8875-8893

An Implicit LU-SGS Scheme for the Discontinuous Galerkin Method on Unstructured Grids

Ma Mingsheng1,2,Gong Xiaoquan1, Deng Youqi1, Zhao Hui1

1.School of Aeronautics, Northwestern Polytechnical University,Xi′an 710072, China 2.China Aerodynamics Research and Development Center,Mianyang 621000, China

The TVD explicit Runge-Kutta discontinuous Galerkin(TVD-RKDG) are widely used in CFD, but it has poor stability property and low computational efficiency. To improve the time advancing efficiency, an efficient implicit lower-upper symmetric Gauss-Seidel (LU-SGS) scheme based on backwards differencing methods (BDF1) has been applied to the high-order discontinuous Galerkin method. The method is implemented parallelly through MPI library and is suit for different computational orders. The left matrix of nonlinear system is simplified and computational efficiency is compared. The DGM is based on the Taylor basis functions on unstructured grid and the accuracy is up to forth order. The developed LU-SGS scheme is verified through NACA0012 airfoil and ONERA M6 wing. Compared to the traditional explicit Runge-Kutta scheme, the implicit scheme can greatly reduce the iteration number and CPU time. The computational efficiency can be accelerated evidently with more than one order of magnitude speed. At last the flow field of DLR-F4 is simulated to verify the robustness of the developed implicit scheme. The result proves that the developed LU-SGS scheme can be applied to complex configuration simulation.

discontinuous Galerkin method; Euler equations; Taylor basis; LU-SGS; computational efficiency; unstructured grid

2016-03-20

马明生(1962—)，西北工业大学兼职教授，主要从事飞行器设计与计算流体力学研究。

V211.3

A

1000-2758(2016)05-0754-07