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]。国内吕宏强等在稀疏的非结构网格上求解了二维Eul-er方程和线化Euler方程^{[7, 8, 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]。

二维守恒形式的Euler方程为

点击浏览公式

点击浏览公式

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

点击浏览公式

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

点击浏览公式

点击浏览公式

点击浏览公式

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

点击浏览公式

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

点击浏览公式

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

点击浏览公式

在每个单元e内，定义为U中最高至p－1阶对应的项求和，定义光顺指示器^[19]形式如下：

点击浏览公式

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

点击浏览公式

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

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

图 1 单元加密 Fig. 1 Element refinement

图选项

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

Bassi指出高阶间断有限元对于物面形的表达精度十分敏感^[3]，因而本文中的物面单元的物面边均采用高阶曲线逼近真实物面：在计算过程中设计物面边的高阶曲线(本文中为6阶)，即使网格十分稀疏，构造好的物面曲边依然可以十分精确的表达真实物面，如图 1中左图虚线所示。

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

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

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

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

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

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

图选项

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

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

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

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

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

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

图 3 单元增加 Fig. 3 Elements increase

图选项

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

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

图 4 单元减少 Fig. 4 Elements decrease

图选项

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

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

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

点击浏览公式

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

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

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

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

图选项

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

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

图选项

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

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

图选项

采用上述自适应间断有限元法给出了经典算例跨声速(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次自适应后的网格变化。

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

图选项

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

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

图选项

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

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

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

图选项

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

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

图选项

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

为验证方法的适应性，采用上述自适应间断有限元法给出了超声速(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)本文设计的自适应方法提高了数值解精度的同时有效地控制了网格量。

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