Ma数预处理间断Galerkin算法
谭勤学, 任静, 蒋洪德
清华大学 热能工程系, 燃气轮机研究所, 北京 100084
通讯作者:任静,教授, E-mail:renj@mail.tsinghua.edu.cn

作者简介: 谭勤学(1986—), 男(汉), 湖南, 博士研究生。

摘要

该文为研究间断Galerkin方法对低 Ma数流动计算的实用性,将有限体积预处理矩阵方法引入间断Galerkin框架,针对低 Ma数问题发展了三维粘性流动求解方法。行波算例(traveling wave)表明: 在间断Galerkin框架下引入预处理矩阵方法可用于低 Ma数粘性流动的计算,且能保持间断Galerkin方法原有的离散精度。顶盖驱动流动、层流边界层、后台阶湍流流动和方腔内自然对流4个经典算例,检验了预处理间断Galerkin方法求解低 Ma数流动的可行性及程序的可靠性。不同 Ma数绕NACA0012无粘流动算例进一步表明,该文所用预处理间断Galerkin方法计算收敛速度几乎与 Ma数无关。

关键词: 间断Galerkin方法; 预处理矩阵方法; 低Ma数; 无积分; 矩阵运算
中图分类号:O354 文献标志码:A 文章编号:1000-0054(2015)01-0134-07
Preconditioning discontinuous Galerkin method for low Mach number flows
Qinxue TAN, Jing REN, Hongde JIANG
Gas Turbine Research Center, Department of Thermal Engineering,Tsinghua University, Beijing 100084, China
Abstract

The finite volume preconditioning method is applied to the discontinuous Galerkin method for three-dimensional viscid low Mach number flows. The traveling wave case is used to verify that the preconditioning discontinuous Galerkin method is suitable for viscid low Mach number flows and that the method retains the original discrete accuracy of the discontinuous Galerkin method. The results for four classical cases (lid driven incompressible flow, a Blasius boundary layer, a backward facing step with turbulent flow and natural convection in a square enclosure) validate the applicability of the preconditioning discontinuous Galerkin method and the reliability of the program. The flow around a NACA0012 airfoil for three different Mach numbers show that the convergence speed is independent of the Mach number.

Keyword: discontinuous Galerkin method; precondition method; low Mach number; Quadrature free; matrix operation

间断 Galerkin方法( discontinuous galerkin method, DGM)最早是由 Reed Hill在1973年提出的,用于中子输运方程的求解[1]。该方法结合了有限体积方法( FVM)和有限元方法( FEM)的基本思想,特别易于处理复杂边界及边值问题,同时具有灵活处理间断的能力,且可以通过适当选取基函数,提高单元插值多项式的次数来实现任意高阶精度[2]。由于这些优异性质,在欧洲用于探讨未来工业设计高阶 CFD( computational fluid dynamies)程序的 ADIGMA计划中[3], 间断 Galerkin方法被认为是最有前途的算法之一[4]

近年来,间断Galerkin方法得到了迅速发展,其中Cockburn和Shu发展的Runge Kutta间断Galerkin方法应用较为广泛[5]。此方法对守恒变量采用时间推进法求解,主要用于可压缩问题的计算。对低 Ma数( Ma<0.3,不可压缩流动可认为 Ma=0)流动问题的求解,由于流体流动速度远小于音速 c, 会导致时间推进计算收敛速度较慢,另外一方面还可能导致计算得到非物理解,给方程的求解带来很大难度[6]。在有限体积方法求解中可采用预处理技术来解决该问题,该类方法通过矩阵变换使系统特征值处于同一数量级。目前,比较具有代表性的预处理矩阵有Turkel预处理矩阵、 Weiss和Smith预处理矩阵等[7]。国际上对于间断Galerkin方法采用预处理矩阵技术求解低 Ma数流动有了一定的研究: Nigro在间断Galerkin方法中使用预处理矩阵,显式Runge Kutta方法及隐式GMRES方法求解了低 Ma数的无粘流动[8]; 但对低 Ma数粘性流动的研究相对较少。

本文基于非结构网格,改进了一种用矩阵运算替代数值积分的间断Galerkin实现方法,并将有限体积预处理矩阵方法引入间断Galerkin方法用于低 Ma数粘性及无粘流动的计算。

1 控制方程

守恒形式的三维 Navier-Stokes方程

Ut+Fxx+Fyy+Fzz=S,

Fx=Fc x-Fv x,

Fy=Fc y-Fv y,

Fz=Fc z-Fv z.(1)

其中: Fc Fv S分别为无粘通量、粘性通量、源项。对一般问题重力可忽略,源项 S=0; 但对于旋转坐标系下的计算需引入离心力及Coriolis源项,对自然对流等问题需引入重力源项。

2 预处理矩阵方法

本文采用 Weiss & Smith的预处理矩阵[9]: 将 Navier-Stokes方程中的守恒变量 U用原始变量 Q代替,并引入预处理矩阵 Г, 控制方程为

(2)

其中:

Q=[p,u,v,w,T]T,Θ=1Ur2-ρTρCp,

Γ= Θ000ρTΘvxρ00ρTvxΘvy0ρ0ρTvyΘvz00ρρTvzΘH-1ρvxρvyρvzρTH+ρCp.

计算过程中 Ur是保证计算精度及稳定性的关键因素,对于可压缩流动,采用

Ur=min( c,max( uref, u, μρh)).(3)

对于不可压流动,采用

Ur=max( uref, u, μρh).(4)

其中: h表示网格尺度; uref为流动中的参考速度。由式(3)可知对于超音速、跨音速等问题 Ur=c, 方程与原始Navier-Stokes方程等价。

2.1 预处理矩阵方法数值通量

在使用预处理矩阵方法求解式(2)时,需要对无粘通量格式进行相应的修正,本文使用预处理 Roe格式,具体形式为

Fcn= 12( FR n+FL n- FcnQ ΔU)≈ 12FRn+FLn-ΓΓ-1FcnQΔQ.(5)

其中: 下标中 n表示法向方向。对于粘性通量,采用中心格式,不针对预处理进行修正。

2.2 预处理矩阵方法特征向量及特征值

预处理 Roe格式或隐式计算时都需要用到无粘通量的 Jacobi矩阵

AΓ-1 (Fc·n)Q =RΓΛΓ RΓ-1.(6)

其中: ΛΓ为对角阵,其特征值、特征向量为

RΓ= 000un-λ2un-λ3ny0-nz-nxρ-nxρ-nxnz0-nyρ-nyρ0-nynx-nzρ-nzρnznxnyun-λ2ρCpun-λ3ρCp,

λ123 =u, λ4 =u'-c', λ5 =u'+c';

u'=u(1), c'= α2u2+Ur2;

α=(1 -βUr2) /2, β=( ρP+ ρTρCP).

3 预处理间断Galerkin方法

将式(2)与测试函数ϕ相乘,并对其在任意单元上分部积分,得到变分形式

(7)

使用有限元空间近似函数 Qh代替函数 Q, DGM重构基函数空间和测试函数空间都取为标准单元上 p阶多项式函数空间 Pp

(8)

上述离散过程中的面积分和体积分采用标准单元上的Gauss数值积分计算: 将积分单元映射到标准单元,然后在标准单元上采用式(9)所示Gauss积分公式, φ表示Lagrange插值多项式, n表示Gauss积分点数量,另外对于 p+1阶精度的DGM, 其体积分精度需要2 p阶,而面积分要求为2 p+1阶,即

(9)

其中:

1) 时间项离散

结合式(8)及式(9),由于其为虚拟时间积分,解的精度仅依赖于空间项的离散,与时间项无关。为简化计算过程,单元内部 Γ的值可用单元内平均值替代而不会影响到计算精度,此时有

(10)

2) 单元体积分离散

结合式(8)及式(9),可得单元体积分过程有

(11)

其中: i=1, ..., N; j=1, ..., Nq; Nq为体数值积分点数量, A B C表达式为

Ai, j= ϕi,jξ Wj, F^x = JξxFx+ξyFy+ξzFz;

Bi, j= ϕi,jη Wj, F^y = JηxFx+ηyFy+ηzFz;

Ci, j= ϕi,jζ Wj, F^z = JζxFx+ζyFy+ζzFz.

3) 单元面积分离散

结合上述式(8)及式(9),单元面积分过程为

(12)

其中: i=1, ..., N; j=1, ..., Ns; Ns为面数值积分点数量; Fn, j为第 j个面数值积分点的数值通量。

总结上述离散过程,采用预处理间断Galerkin方法所得半离散形式方程为

(13)

本文采用3步Runge Kutta方法求解该方程,在此不赘述该方法。另外从上述公式可看到,此种以矩阵向量乘积替代积分过程的方法,在前处理阶段预先存储相关矩阵,可有效减少计算量、提高计算效率,并可简化程序编写过程。经本文程序实测,此种实现方法相对直接数值积分,计算效率有一定提高。

4 预处理计算方法验证

在此算例中,计算域为(0,1)×(0,1)控制方程为不可压缩 Navier-Stokes方程, x和 y方向均为周期性边界 10其初始有2个行波,行波在粘性作用下逐渐耗散,其解析解为

u(x,y,t)=1+2 cos(2π(x-t)) sin(2π(y-t)) e-8π2vt, (14)

v(x,y,t)=1-2 sin(2π(x-t)) cos(2π(y-t)) e-8π2vt, (15)

p(x,y,t)=-( cos(4π(x-t))+ cos(4π(y-t))) e-16π2vt. (16)

本文计算v=0.01, 为考察本文所采用预处理矩阵方法计算低Ma数及不可压流动问题对计算精度的保持性,采用一阶多项式重构,对流通量格式为预处理 Roe格式,计算到时刻t=1.0 s图1图2图3表明对于压力和速度, DGM计算结果与理论结果定性上吻合很好,证明文中所采用预处理矩阵方法可较好用于低Ma数问题的计算。

图1 理论压力 p分布与计算结果对比

图2 理论压力 x方向速度 u分布与计算结果对比

图3 理论压力 y方向速度 v分布与计算结果对比

表1定量给出了L2误差分析结果,此时速度、压力计算精度均能达到 DGM一阶插值多项式原有的2阶离散精度。此算例定性及定量上表明: 在 DGM框架下引入预处理矩阵方法,可较好计算低Ma数问题,且能保持 DGM的离散精度。

表1 一阶插值多项式 L2误差分析
5 算例

为进一步检验预处理矩阵方法应用于高精度 DGM求解低Ma数问题的可行性及程序的可靠性: 本文计算了顶盖驱动方腔不可压缩粘性流动、 Rayleigh-Benard自然对流、 Blusius层流边界层、后台阶湍流流动以及绕 NACA0012流动等经典算例。

5.1 顶盖驱动不可压粘性流动

计算中Re=uL/v=1 000, 网格数量为60×60的非均匀六面体单元,最小尺寸0.0002m, 采用2阶多项式插值。边界条件: 顶盖为切向速度u的壁面,其它为静止壁面,计算得到方腔内流线如图4所示,中间水平面上的V速度随x变化分布及中间垂直面上的U速度随y变化分布与文[11]的比较如图5所示,可看到本文计算结果和经典结果吻合很好。

图5 无量纲 U, V速度分布与参考值比较[10]

5.2 层流平板边界层

本算例计算长度L=1.0 m, 计算域为(0,0.1)×(0,1)的矩形区域,其中基于平板长度的Re数为106, 网格数为60×100, 对壁面边界进行加密,最小网格尺度为0.00005 m, 1阶多项式重构; 给定进口总温总压,出口静压。定义无量纲速度

η=y Rex12/x, Rex=ux/v. (17)

图6给出了x=0.9 m处,边界层内速度分布与经典 Blasius速度分布的比较,可看到数值解与理论解一致。

图6 边界层内速度分布与经典Blasius比较

5.3 后台阶湍流流动

后台阶流动为测试湍流计算的经典算例。本文所用计算模型及边界条件如图7, 进口Ma数约0.128, 1阶多项式重构,湍流模型有 Spalart-Allmaras一方程模型及 Menter SST双方程模型。实验中分离再附点位置为x/H=6.26±0.1, 本文 Spalart-Allmaras模型计算结果为6.19, Menter SST模型计算结果为6.39; x方向速度在y方向的分布与实验对比如图8所示(图中U/U ref的x坐标取为物理位置起点x值)[12]。由计算结果可知,分离再附位置及速度分布都和实验吻合较好,验证了本文预处理 DGM对该低Ma数湍流流动计算的可行性及程序湍流计算的正确性。

图7 后台阶计算模型及边界条件

图8 x方向速度分布与实验对比图

5.4 封闭方腔内自然对流

本算例H=0.1 m, 计算域为(0,0.1)×(0,0.1), 网格为50×50的壁面加密网格,使用理想气体模型并考虑重力, 1阶多项式重构。边界条件: 左右壁面温度分别为T C和T H, 上下壁面绝热,计算中所用到的无量参数定义为

Ra= gβPr(TH-TC)H3v2,U= uHν/Pr,V= vHv/Pr,Nu=H Θx,Θ= T-TcTH-Tc.

计算得到无量纲温度及无量纲速度分布如图9图10图11所示,与 Davis[13]结果比较可知,本文中预处理 DGM计算得到的分布情况与其一致。表2定量上给出了误差情况,其速度最大值及最小值的误差在1.0 %左右,而Nu数的误差小于 0.5 %,而Ra=106误差较大,其原因是其流动趋向湍流。

图9 自然对流不同 Ra数无量纲温度 Θ分布

图10 自然对流不同 Ra数无量纲x方向速度 U/URef分布

图11 自然对流不同 Ra数无量纲 y方向速度 V/VRef分布

表2 间断Galerkin误差分析
5.5 绕NACA0012翼型无粘计算

本算例为绕 NACA0012翼型无粘流动,计算所用网格单元数量为112×32全六面体网格,进口攻角0°, 一阶多项式重构。远场来流 Mach数为:Ma=10-1、 10-2、 10-3, 计算所得Ma分布如图12所示,与文[8]中相关计算结果一致,验证了计算的准确性; 其对应的压力方程收敛历史如图13: 此时低Ma数预处理 DGM收敛速度几乎与Ma数无关。

图12 不同 Ma数远场来流计算所得相对 Ma数分布

6 结论

本文发展了用矩阵运算代替数值积分的间断 Galerkin实现方法,使程序易于编写,并能有效地提高程序效率。针对本文预处理矩阵方法给出了其特征向量,并将预处理矩阵方法引入 DGM求解框架,相关算例表明:

1) 行波算例验证了在间断 Galerkin框架下引入预处理矩阵方法,可合理计算低Ma数问题,且能保持间断 Galerkin的离散精度。

2) 顶盖驱动方腔不可压流动、层流平板边界层、后台阶湍流流动、低Ra数下的自然对流问题从不可压、低Ma数层流流动、低Ma数湍流流动、浮升力驱动低Ma数流动这4个层面的进一步检验了预处理矩阵方法用于间断 Galerkin求解不可压及低Ma数问题的可行性及程序的可靠性。

3) 不同Ma数绕 NACA0012无粘计算表明: 采用本文所述预处理矩阵方法计算低Ma数问题,计算收敛速度几乎与Ma数无关。

The authors have declared that no competing interests exist.

参考文献
[1] Reed W H, Hill T R. Triangular Mesh Methods for the Neutron Transport Equation, Report LA-UR-73-479 [R]. New Mexico, USA: Los Alamos Scientific Laboratory, 1973 [本文引用:1]
[2] 刘儒勋, 舒其望. 计算流体力学的若干新方法 [M]. 北京: 科学出版社, 2003.
LIU Ruxun, Shu C W. Kinds of New Methods in Computational Fluid Dynamics [M]. Beijing: Science Press, 2003. (in Chinese) [本文引用:1]
[3] Kroll N. ADIGMA—A European project on the development of adaptive higher-order variational methods for aerospace applications [C]//47th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Orlando, USA: The American Institute of Aeronautics and Astronautics (AIAA)2009. [本文引用:1]
[4] 谭勤学, 任静, 蒋洪德. 全矩阵无积分间断有限元实现[J]. 中国科技论文, 2013, 8(8): 776-779
TAN Qinxue, REN Jing, JIANG Hongde. Quadrature-free implementation of discontinuous Galerkin method using the matrix product formula[J]. China Science Paper, 2013, 8(8): 776-779 [本文引用:1]
[5] Cockburn B, Hou S, Shu C W, et al. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case[J]. Mathematics of Computation, 1990, 54(190): 545-581. [本文引用:1] [JCR: 1.409]
[6] 李雪松, 顾春伟, 徐建中. 新型高低速流统一算法及其叶轮机应用[J]. 航空学报, 2007, 28(6): 1334-1338.
LI Xuesong, GU Chunwei, XU Jianzhong. New general scheme for flows at all speeds and it's application to turbomachinery[J]. Acta Aeronau Tica ET Astronau Tica Sinica, 2007, 28(6): 1334-1338. (in Chinese) [本文引用:1]
[7] 杨果. 直接间断有限元方法求解奇异摄动问题 [D]. 长沙: 湖南师范大学, 2009
YANG Guo. The Direct Discontinuous Galerkin Methods for Singularly Perturbed Problems [D]. Changsha: Hunan Normal University, 2009 [本文引用:1]
[8] Nigro A, De Bartolo C, Hartmann R, et al. Discontinuous Galerkin solution of preconditioned Euler equations for very low Mach number flows[J]. International Journal for Numerical Methods in Fluids, 2010, 63(4): 449-467. [本文引用:2] [JCR: 1.329]
[9] Gustavson J, Weiss J M, Smith W A. Preconditioning applied to variable and constant density flows[J]. AIAA Journal, 1995, 33(11): 2050-2057. [本文引用:1] [JCR: 1.165]
[10] Brown D L. Performance of under-resolved two-dimensional incompressible flow simulations, II[J]. Journal of Computational Physics, 1997, 138(2): 734-765. [本文引用:1] [JCR: 2.485]
[11] Ghia U, Ghia K N, Shin C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. Journal of Computational Physics, 1982, 48(3): 387-411. [本文引用:1] [JCR: 2.485]
[12] Driver D M, Driver D M, Seegmiller H L. Features of a reattaching turbulent shear layer in divergent channel flow[J]. AIAA Journal, 1985, 23(2): 163-171. [本文引用:1] [JCR: 1.165]
[13] de Vahl Davis G. Natural convection of air in a square cavity: a bench mark numerical solution[J]. International Journal for Numerical Methods in Fluids, 1983, 3(3): 249-264. [本文引用:1] [JCR: 1.329]