基于高阶物面近似的自适应间断有限元法欧拉方程数值模拟

孙 强，吕宏强，伍贻兆

(南京航空航天大学，航空宇航学院，江苏南京 210016)

基于高阶物面近似的自适应间断有限元法欧拉方程数值模拟

孙 强，吕宏强*，伍贻兆

(南京航空航天大学，航空宇航学院，江苏南京 210016)

将高阶间断有限元与网格自适应相结合，于非结构网格上求解二维Euler方程。将数值解多项式的高阶项贡献用人工粘性项系数的形式进行量化，网格自适应过程中以人工粘性项系数作为网格自适应的指示器。在系数达到设定的上限阀值的区域进行网格加密，在系数达到设定的下限阀值的区域将迭代过程中加密过的网格稀疏以减少网格量。所有自适应均在高阶曲线逼近真实物面的基础上进行，以保证数值结果的精度。典型数值算例结果与实验结果进行了对比，表明采用该自适应间断有限元法可以保证以尽可能小的计算量得到高精度结果。

高阶间断有限元;自适应方法;Euler方程;人工粘性;物面高阶近似

0 引言

间断Galerkin有限元法(DGM)发展于20世纪70年代［1］。由于该方法能够构造高精度的格式，并且易于实现自适应计算和并行计算，已成为计算流体力学领域研究热点之一。Bassi等在1997年用DGM成功求解了二维Euler方程和N-S方程［2-3］;Cockburn等提出Runge-Kutta间断Galerkin方法(RKDG)［4-5］; Nguyen等将S-A湍流模型成功应用于DGM［6］。国内吕宏强等在稀疏的非结构网格上求解了二维Euler方程和线化Euler方程［7-9］;张来平等发展了静动态混合重构的DG/FV混合格式［10］;于剑等对人工粘性和NS方程进行了研究［11-12］。

自适应方法根据研究对象的不同和解的性质自动调整网格分布或者数值解的逼近精度。在流体力学计算中，解通常在局部变化较大，其他大部分区域是光滑的［12］，此时使用自适应方法可以有效降低计算量。自适应方法通常被分为三类:对局部网格加密或稀疏的h型自适应;调整局部单元内逼近多项式次数的p型自适应;调整网格单元位置保持单元总数不变的r型自适应。

由于间断有限元方法通常将未知数表达为高阶形式，在相同的网格上计算量与有限体积法和有限差分法相比成倍增加，故而自适应方法与间断有限元方法相结合势在必行。Yang等提出了一种于非结构网格上模拟界面的自适应方法［13］;Flaherty和Remacle研究小组则首次将网格度量场引入自适应间断有限元方法中［14］。徐云、吴迪、蔚喜军［15-16］等则根据后误差估计提出了不同的h型网格自适应方法。

在高阶情况下(p≥2)，间断有限元对于物面形状的表达十分敏感［3］，若物面边界以直线表达，则无法得到收敛解。本文设计了一种h型自适应高阶间断有限元法，物面边均采用高阶曲线逼近真实物面以保证数值解的精度。将数值解多项式的高阶项贡献以人工粘性［17］系数ε的形式量化，以系数ε作为网格自适应的指示器。数值模拟过程中产生的非线性方程组采用牛顿法进行迭代，每个迭代步产生的大型线性系统采用块高斯-赛德尔方法求解［18］。

1 间断有限元数值离散

二维守恒形式的Euler方程为

其中，U为守恒变量，对流通量F(U)=(f(U)，g (U))，表达式分别为:

式中，ρ为密度;v1和v2为x、y方向的速度分量;E为单位总能;P为压强;h为单位总焓。

式(1)两边乘检验函数V，在计算域运用分部积分后，弱解形式为:

将计算域划分为网格{e}，其中e表示网格单元。在每个单元e内:

式(4)中，φi(x，y)为基函数。将式(4)代入式(3)中，得

式(6)称为p阶间断有限元离散［7］，p为式(4)中基函数的阶数。

由于间断有限元法允许解在单元的边界处不连续，因而数值通量的定义不是唯一的。此处的数值通量直接借用传统有限体积法的处理方法，即将式(6)中的通量函数F(Uh)·n用数值通量代替，其中和分别表示该单元和其相邻单元在该边界处的值。本文所有计算都采用LLF格式数值通量:

在物面和远场分别采用无穿透和无反射边界条件。式(6)最终的离散形式可写为:

其中，M是质量矩阵，如果间断有限元采用正交基函数，则是对角矩阵;R(u)是残值;u包含Uh中所有的未知系数。对式(8)的求解采用牛顿-块高斯赛德尔方法求解［18］。

2 人工粘性

人工粘性法［18］的思想就是在原Euler方程中加入一个耗散修正项，从而抑制激波附近的非物理振荡。方程(1)变为如下形式:

式(9)中，参数ε控制粘性项的大小。

式中，(·，·)e表示L2(Ωe)空间的标准内积。这里预估Se的值能达到1/p4数量级。

单元的光顺指示器确定后，粘性项系数通过以下函数定义:

式(11)中，se= log10Se，参数 ε0= De/p，s0= log10(1/p4)，κ是人为定义常数。其中，De是计算元单元e的特征尺度。

由式(10)和式(11)可以看出，se是由U中高阶项决定，表征了U的高阶项的贡献大小。解的光顺程度越差se越大。εe只是作为se的分段函数，在数值表征上意义相同，它抹平光顺区域的较小的se，突出非光顺区域，使计算和表达更简便。

3 网格自适应

3.1 网格加密

在计算过程中单元的ε达到设定好的上限阀值则认为该单元需要进行加密，采用三角形单元四分法对网格单元进行加密，如图1所示。

图1 单元加密Fig.1 Element refinement

图中实线单元为非物面单元，可以直接取三条边的中点连线剖分加密。

Bassi指出高阶间断有限元对于物面形的表达精度十分敏感［3］，因而本文中的物面单元的物面边均采用高阶曲线逼近真实物面:在计算过程中设计物面边

的高阶曲线(本文中为6阶)，即使网格十分稀疏，构造好的物面曲边依然可以十分精确的表达真实物面，如图1中左图虚线所示。

在物面单元需要进行剖分时，利用已构造好的高阶曲线找到逼近真实物面的物面边的中点，然后连接另外两条边的中点生成四个新的网格单元，并且新生成的两个物面单元的物面边也用高阶曲线拟合，如图1中右图虚线所示。

网格加密需要遵守以下原则:

1)相邻单元的被剖分次数相差最多为1。

2)当单元的几何尺寸细化到设定好的阀值之后就不可继续剖分。

3.2 网格稀疏

在计算过程中，部分单元被加密后可能在进一步的迭代过程中解变得十分光滑，需要对其进行粗化以减少网格量。当加密过的单元的四个子单元ε均达到设定好的下限阀值时，则这四个子单元进行合并。如图2所示。

图2 单元的粗化Fig.2 Element coarsening

图2中，实线表示非物面单元的合并过程。若粗化后的单元为物面单元，物面边依旧采用高阶曲线逼近真实物面，如图中虚线所示。

网格粗化是网格加密的逆向操作，需要遵守以下原则:

1)初始单元不可被粗化。

2)只有被加密过的单元才可以被粗化。

3)当邻单元已经比判断需要粗化的四个单元剖分次数大1时，不可进行粗化。

3.3 网格数据结构

为使程序的可移植性更高，将网格自适应作为模块进行操作。虽然二叉树结构存储加密和粗化网格直观高效，但自适应过程中增加与减少的点和线编号复杂且对于程序的结构不利。因而本文所有单元、点、线的存储结构均为链表结构，加密过程中单元的增加操作如图3所示。

E表示需要加密的单元，生成的四个子单元中间单元依旧为编号E，其余三个则排在总单元后三个，点与边的操作类似。在这种操作下，网格模块与程序流场计算模块的接口不受自适应操作的影响，从而提高程序的可移植性。

同样的，网格稀疏的过程单元减少过程如图4所示。

图3 单元增加Fig.3 Elements increase

图4 单元减少Fig.4 Elements decrease

其中E和三个fi表示由加密操作生成的四个子单元，合并后生成的单元依旧编号为E，而fi后的单元编号依次减3。

在网格加密过程中，每个子单元均设置两个标记，用于记录所在的初始单元号和本身的加密次数。在后续的网格稀疏过程中，根据这两个信息和单元相邻关系进行单元合并。

3.4 数据传递

网格加密和粗化的过程中，使各物理量在子单元内的积分值等于合并单元内的积分值，从而保证物理量的守恒。假设合并单元为E，各子单元为ei，uE(x，y)为合并单元的守恒变量的函数，uei(x，y)为子单元守恒变量的函数，则需满足以下方程:

其中m为合并单元剖分后的子单元个数。

4 数值结果

4.1 圆柱绕流测试算例

为验证物面构造与网格自适应结合的效果，本文首先对圆柱绕流(Ma=0.38、α=0)这一经典算例进行了数值模拟。

前人［3，8］做圆柱绕流计算为获得精度较高的结果，网格均人为设计为规则分布。本文为验证物面构造和自适应间断有限元的优点，初始网格使用了分布较随意、单元尺寸较极端的计算网格。

图5给出了圆柱绕流网格，该网格共有80个单元，物面仅有8个点。

图5 圆柱绕流计算网格Fig.5 Mesh around a circle

图6 ε分布和自适应后的网格Fig.6 ε distribution and mesh after three times adaption

图6(a)为在原始网格上p=4时ε分布，可以看出主要分布在流场数值可能出现数值解变化比较剧烈的区域。对这一区域进行网格自适应，2次自适应后的网格如图6(b)所示。可以看出网格主要在ε集中区域自适应，从而仅需对这一区域的少量网格进行自适应操作。自适应前后的网格对比，发现自适应操作确实基于真实物面曲线，自适应后的网格更加逼近真实物面。

图7给出了p=4情况下的圆柱绕流等马赫线图和压强系数分布，可以看出压强系数分布光滑并且对称。

图7 圆柱绕流等马赫线图和Cp分布Fig.7 Mach isolines and Cpdistribution

4.2 跨声速NACA0012翼型绕流

采用上述自适应间断有限元法给出了经典算例跨声速(Ma=0.8、α=1.25)下NACA0012翼型绕流的数值结果。

图8给出的NACA0012翼型的网格共有448个单元，而物面单元只有32个。

图8 NACA0012翼型计算网格Fig.8 Mesh around NACA0012 airfoil

图9给出了引入人工粘性未加入网格自适应的粘性项系数ε的分布与马赫云图。

图9 未加自适应的ε分布及等马赫线图Fig.9 ε distribution and Mach isolines without adaption

从图中可以看出虽然网格十分稀疏，但依然可以在一个网格尺度内捕捉到激波，并且ε的分布集中在激波区域。

图10给出了1次、2次自适应后的网格变化。图11给出了网格自适应过程中局部网格的细化和粗化的动态变化情况。

从图中可以看出，自适应过程中激波区域的网格加密，而远离激波区域在精度低的时候加密了的单元则合并降低网格量。

图12给出了自适应结束时的网格分布和ε分布。

图13给出了添加自适应后的五阶精度(p=4)的等马赫线图和表面压强系数分布。

图10 1次和2次自适应后的网格Fig.10 Mesh after once and twice adaption

图11 1次和2次自适应后的局部网格Fig.11 Local mesh after once and twice adaption

图12 自适应结束的网格和ε分布Fig.12 Mesh and ε distribution after adaption

图13 等马赫线图和Cp分布曲线Fig.13 Mach isolines and Cpdistribution

从图中可以看出，自适应后的数值结果较原始网格上的数值精度有显著提高，并且激波捕捉更精确接近实验值。

4.3 超声速NACA0012翼型绕流

为验证方法的适应性，采用上述自适应间断有限元法给出了超声速(Ma=5.0、α=0)下NACA0012翼型绕流的数值结果(p=4)。

初始网格如图14，初始网格数为660。自适应计算结束后的马赫云图及网格、ε分布如图15所示。

从图中可以看出，自适应操作主要集中在头部激波区域。

图14 局部网格Fig.14 Local mesh

图15 马赫云图及ε分布Fig.15 Mach contour and ε distribution

5 结论

(1)将数值解的高阶项贡献量化后作为自适应指示器可行并直观有效。

(2)即使在很稀疏的初始网格基础上，使用基于真实物面逼近的网格自适应方法也可以获得高精度的数值解。

(3)本文设计的自适应方法提高了数值解精度的同时有效地控制了网格量。

对于非定常问题和三维问题，自适应方法在节省网格工作量和提高计算效率方面作用更为突出。同时具有加密和粗化的功能在非定常计算中会得到更有意义的应用，该部分研究正在进行中。在粘性流动的情况下可以在不引入人工粘性项的前提下计算人工粘性项系数，将其作为网格自适应判据，这部分工作也正在进行。

［1］ Reed W H，Hill T R.Triangular mesh methods for the neutron transport［R］.Los Alamos Scientific Laboratory Report LA-UR-73-479，1973.

［2］ Bassi F，Rebay S.A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations［J］.Journal of Computational Physics，1997，131 (2):267-279.

［3］ Bassi F，Rebay S.High-order accurate discontinuous finite element solution of the 2D Euler equations［J］.Journal of Computational Physics，1997，138(2):251-285.

［4］ Cockburn B，Shu C.The Runge-Kutta discontinuous Galerkin method for conservation laws V:multidimensional systems［J］.Journal of Computational Physics，1998，141(2):199-224.

［5］ Cockburn B，Shu C.The local discontinuous Galerkin method for time-dependent convection-diffusion systems［J］.SIAM Journal on Numerical Analysis，1998，35(6):2440-2463.

［6］ Nguyen N C，Persson P O，Peraire J.RANS solutions using high order discontinuous Galerkin methods［R］.AIAA 2007-0914.

［7］ Lu Hongqiang，Wu Yizhao，Zhou Chunhua，et al.High resolution of subsonic flows on coarse grids［J］.Acta Aeronautica et Astronautica Sinica，2009，30(2):200-204.(in Chinese).

吕宏强，伍贻兆，周春华，等.稀疏非结构网格上的亚音速流高精度数值模拟［J］.航空学报，2009，30(2):200-204.

［8］ Lu Hongqiang，Zhu Guoxiang，Song Jiangyong，et al.High-order discontinuous Galerkin solution of linearized Euler equations［J］.Chinese Journal of Theoretical and Applied Mechanics，2011，43 (3):621-624.(in Chinese).

吕宏强，朱国祥，宋江勇，等.线化欧拉方程的高阶间断有限元数值解法研究［J］.力学学报，2011，43(3):621-624.

［9］ Lu H.High-order discontinuous Galerkin solution of low-Re viscous flows［J］.Modern Physics Letters B，2009，23(03):309-312.

［10］Zhang Laiping，Liu Wei，He Lixin，et al.A class of discontinuous Galerkin/finite volume hybrid schemes based on the“static re-construction”and“dynamic re-construction”［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(6):1013-1022.(in Chinese).

张来平，刘伟，贺立新，等.基于静动态混合重构的DG/FV混合格式［J］.力学学报，2010，42(6):1013-1022.

［11］Yu J，Yan C.An artificial diffusivity discontinuous Galerkin scheme for discontinuous flows［J］.Computers＆Fluids，2013，75(20): 56-71.

［12］Yu Jian，Yan Chao.Study on discontinuous Galerkin method for navier-stokes equations［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(5):962-969.(in Chinese).

于剑，阎超.Navier-Stokes方程间断Galerkin有限元方法研究［J］.力学学报，2010，42(5):962-969.

［13］Yang X，James A J，Lowengrub J，et al.An adaptive coupled levelset/volume-of-fluid interface capturing method for unstructured triangular grids［J］.Journal of Computational Physics，2006，217 (2):364-394.

［14］Remacle J F，Li X，Shephard M S，et al.Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods［J］.International Journal for Numerical Methods in Engineering，2005，62(7):899-923.

［15］Xu Yun，Yu Xijun.Adaptive discontinuous Galerkin methods for hyperbolic conservation laws［J］.Chinese Journal of Computational Physics，2009，26(2):159-168.(in Chinese).

徐云，蔚喜军.自适应间断有限元方法求解双曲守恒律方程［J］.计算物理，2009，26(2):159-168.

［16］Wu Di，Yu Xijun.Adaptive discontinuous Galerkin method for euler equations［J］.Chinese Journal of Computational Physics，2010，27 (004):492-500(in Chinese).

吴迪，蔚喜军.自适应间断有限元方法求解三维欧拉方程［J］.计算物理，2010，27(004):492-500.

［17］VonNeumann J，Richtmyer R D.A method for the numerical calculation of hydrodynamic shocks［J］.Journal of Applied Physics，1950，21(3):232-237.

［18］Xia Y D，Wu Y Z，Lv H Q，et al.Parallel computation of a highorder discontinuous Galerkin method on unstructured grids［J］.Acta Aerodynamica Sinica，2011，29(5):537-541.(in Chinese).

夏轶栋，伍贻兆，吕宏强，等.高阶间断有限元法的并行计算研究［J］.空气动力学学报，2011，29(5):537-541.

［19］Lu H，Berzins M，Goodyer C E，et al.Adaptive high-order discontinuous Galerkin solution of elastohydrodynamic lubrication point contact problems［J］.Advances in Engineering Software，2012，45 (1):313-324.

Adaptive discontinuous Galerkin method to solve Euler equations based on high-order approximative boundary

Sun Qiang，Lu Hongqiang*，Wu Yizhao
(College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China)

A high-order discontinuous method(DGM)is integrated with adaptive method to solve Euler equations on unstructured mesh.Contribution of the polynomial’s highest-order terms is quantified in the form of artificial viscous coefficient.The coefficient is regarded as the indicator of h-adaptivity.Elements where the coefficients are greater than the upper limit are refined.Those where the coefficients are less than the lower limit are coarsened if they have been refined.A high-order geometric approximation of curved boundaries is adopted to ensure the convergence.Numerical results of test cases are consistent with corresponding experimental ones.High accurate numerical results can be obtained with the h-adaptive method at low expense.

high-order DGM;adaptive method;Euler equations;artificial viscosity;high-order boundary approximation

V211.3

Adoi:10.7638/kqdlxxb-2014.0027

0258-1825(2015)04-0446-08

2014-04-23;

2014-07-08

国家自然科学基金(11272152);航空科学基金(20101552018);江苏高校优势学科建设工程资助项目

孙强(1987-)，男，山东莱阳市人，博士研究生，研究方向:计算流体力学.E-mail:submarineultra@126.com

吕宏强*(1977-)，男，博士，副教授，E-mail:hongqiang.lu@nuaa.edu.cn

孙强，吕宏强，伍贻兆.基于高阶物面近似的自适应间断有限元法欧拉方程数值模拟［J］.空气动力学学报，2015，33(4):446-453.

10.7638/kqdlxxb-2014.0027 Sun Q，Lu H Q，Wu Y Z.Adaptive discontinuous Galerkin method to solve Euler equations based on high-order approximative boundary［J］.Acta Aerodynamica Sinica，2015，33(4):446-453.