2. 中国空气动力研究与发展中心计算空气动力研究所, 四川 绵阳 621000
随着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) |
对方程(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) |
方程(11)左端项的面积分通量Resfn, 通过简单的LF数值通量格式计算
(13) |
进一步得到面积分通量对守恒变量的偏导数
(14) |
由方程(8)看到, 体积分通量为对流项与基函数梯度的点乘, 这里将基函数梯度类似看为面的法向量, 从而能够得到
将方程(8)的第一项以及方程(14)带入方程(12)得
(15) |
以三维p=3(四阶精度)DGM为例, 每个单元的左端项矩阵是大小为(5×20)×(5×20)的矩阵。令
(16) |
由于
(17) |
从方程(17)看到, 一般情况下矩阵
本文首先通过等熵涡算例验证了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, 可以看到方法达到了设计的精度。
精度 | 网格 | L1误差 | 阶数 | L2误差 | 阶数 |
DG(p=1) | 10*10*10 | 3.03×10-3 | 8.52×10-3 | ||
20*20*20 | 6.36×10-4 | 2.254 | 2.18×10-3 | 1.961 | |
40*40*40 | 1.47×10-4 | 2.107 | 5.17×10-4 | 2.080 | |
DG(p=2) | 10*10*10 | 5.01×10-4 | 1.23×10-3 | ||
20*20*20 | 6.10×10-5 | 3.037 | 1.54×10-4 | 3.000 | |
40*40*40 | 8.33×10-6 | 2.872 | 1.89×10-5 | 3.025 | |
DG(p=3) | 10*10*10 | 2.01×10-4 | 5.13×10-4 | ||
20*20*20 | 9.99×10-6 | 4.331 | 3.00×10-5 | 4.095 | |
40*40*40 | 5.77×10-7 | 4.112 | 1.80×10-6 | 4.058 |
为考察DGM的精度以及发展的隐式算法的收敛性, 本文计算了NACA0012翼型在来流马赫数Ma=0.5时亚声速绕流。计算网格如图 1所示。
对称面共7 200个三角形单元, 通过展向拉伸对称面单元生成三棱柱单元。计算采用48个CPU并行, 给出了2~4阶DGM结果, 并与二阶的FVM进行了对比, 图 2给出了二阶RKDG、不同精度隐式DGM以及二阶FVM的残差收敛曲线, 此时显式的CFL=0.3, 隐式计算的CFL=10。从图看到, 相同二阶精度下, FVM和DGM残差收敛速度相当, 隐式收敛效率远高于显式。显示了二阶精度下, 不同CFL数的收敛曲线。从图看到, CFL数大于10后, 收敛速度基本不变, 从
从图 3a)的曲线1和5, CFL=1, 简化左端项矩阵对收敛速度有一定影响, 导致迭代步数略有增加。从图 3b)看到, 对于CFL=1的情况, 由于简化左端项矩阵后, 单步计算时间减小, 总的收敛时间减小。由于简化左端项矩阵后CFL数不能取大, 因此未简化左端项矩阵当CFL数取大时, 其收敛速度快于简化后的状态。
2.3 ONERA M6机翼亚声速绕流为进一步考察LU-SGS方法在三维外形计算效率, 本文计算了ONERA-M6亚声速流场, 计算马赫数Ma=0.5, 迎角0°。网格如图 4所示, 网格共约14万个四面体单元, 机翼前后缘及空间分别采用各向异性和各向同性四面体网格。计算采用64个CPU。
图 5给出了DGM不同计算精度下LU-SGS残差收敛曲线, 同时给出了FVM以及RKDG显式结果。隐式计算的CFL=1 000, 显式计算的CFL=0.3。从残差随迭代步数图看, 发展的LU-SGS隐式方法收敛效率远高于显式方法。二阶精度下, 由于DGM计算的自由度为FVM的4倍且计算面通量和体积分通量时采用更多的高斯积分点, 二阶DGM的单步计算时间远大于FVM。
2.4 DLR-F4翼身组合体亚声速绕流为验证发展的隐式算法的初步工程计算能力, 数值模拟了DLR-F4翼身组合体亚声速绕流。计算网格如图 6所示。
四面体网格单元总数为637 126个, 计算马赫数Ma=0.5, 来流迎角为0°, 计算采用128个CPU, 图 7给出了残差收敛曲线, 相同CFL数以及相同的dt计算方法下, 相对于有线体积方法, DGM残差收敛速度较快, 但是下降量级较少。这是由于在迭代过程中DGM一直修改迭代得到的自由度, 造成限制单元残差不降, 因此需要进一步研究该方法的限制器。
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. DOI:10.1016/j.jcp.2007.03.002 |
[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. DOI:10.4208/cicp |
[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. |
[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. |
[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. DOI:10.1016/j.jcp.2008.06.035 |
2. China Aerodynamics Research and Development Center, Mianyang 621000, China