2. 中国空气动力研究与发展中心 计算空气动力研究所, 四川 绵阳 621000
2. Computational Aerodynamics Institute, China Aerodynamics Research and Development Center, Mianyang, Sichuan 621000, China
随着计算机技术和计算流体力学(CFD)的发展，近几十年来网格生成技术及数值计算方法取得了飞速的进步。在网格生成技术中，非结构网格具有适合复杂外形、方便网格自适应等突出优点，进一步发展的混合网格技术部分克服了非结构网格需要大量计算资源的不足，代表了网格技术的发展趋势[1]。而对数值计算方法来说，目前已知的绝大多数CFD商业软件和in-house工业应用软件均以二阶精度的有限体积方法(FVM)为基础，尽管其已成功解决了大量的工程实际问题，但在CFD领域内仍有许多问题需要高阶精度(三阶或以上精度)方法才能较好地解决。这主要是因为低阶方法具有较大的数值耗散与色散，对一些非常复杂的流动现象，如旋涡主导的流动、分离、湍流等问题，其通常难以给出精细的流场结构。尤其对于湍流的大涡模拟(LES)、直接数值模拟(DNS)等问题，低阶方法固有的数值耗散可能大到掩盖这些真实的物理粘性，采用高阶方法已是CFD界的共识[2]。在计算气动声学、计算电磁学等领域，需要对长时间历程的波传播进行高精度的数值模拟，低阶方法在网格规模受计算机资源限制时往往无法准确模拟流场的特性，带来无法接受的计算结果，此时亦有必要使用高阶方法。此外分析表明，对于误差水平要求不高的问题，低阶方法可以用较小的计算代价来解决；而对误差水平要求很高的问题，使用高阶方法计算代价更小，效率更高[2]。
近二十年来，基于非结构/混合网格的高阶精度计算方法发展迅速，CFD工作者已提出大量的计算方法，主要包括k-exact FVM[3, 4]、间断Galerkin方法(DGM)[5, 6]、有限谱体积(SV)方法[7]、有限谱差分 (SD)方法[8]、及将DG、SV、SD等方法统一在一个框架之内的CPR(Correction Procedure via Reconstruction)方法[9]等等。更多相关内容可以参 考Ekaterinaris[10]、Wang[2]和张来平等[11]的综述文章。
高阶k-exact FV方法和以DG方法为代表的DG/SV/SD/CPR等方法都各具优点，但仍有可以改 进的空间。如高阶FVM需要扩充模板来提高重构精度，在非结构网格上模板的搜寻扩展很不方便并且方法不紧致，而DGM的计算量和存储量非常大。由此构造结合FV方法和DG方法优点的混合方法是一种较好的选择，其已受到许多学者的关注，并提出了多种混合方案。这些混合方法的核心思想是使用本单元和相邻单元的多个自由度重构一个更高阶的分布，使用这个更高阶的分布来得到本单元的较低阶自由度信息的高阶更新。
Cockburn等[12]为了进一步提高DG方法精度最早提出对解使用重构算法的思想，后来Ryan等[13]做了进一步的发展。上述工作仅仅是在最后输出解时使用重构算法，因而可以认为是DG方法的后处理技术。Dumbser和Munz首次提出从DG方法计算开始时刻即使用线性重构算子，并构造出一类重构的DG方法，命名为PNPM方法[14, 15]；基于类似思想，Luo等提出了RDG[16, 17](Reconstructed DG)方法，构造了3阶RDG(P1P2)格式；张来平等构造出一类DG/FV混合方法[18, 19, 20, 21, 22]；在原始CPR格式的基础上，王志坚等借鉴PNPM方法和DG/FV混合方法的思想，构造了一系列PNPM-CPR格式[23, 24]。这几种混合方法均显示出了很好的性能，具有每个自由度比FV方法和DG方法更高效、方法紧致、隐式时间离散内存需求更低等很多优点[11]。
张来平等在构造DG/FV混合方法[18, 19, 20, 21, 22]时提出了“静态重构”和“动态重构”的概念，并基于静动态“混合重构”的思想构造了一类三阶以上精度的DG/FV混合格式，其已应用于三角形/四边形混合网格下的标量方程和Euler/N-S方程的数值模拟，大量的数值算例表明该方法达到了设计精度，且相比同阶精度DG方法具有更高的计算效率[25]。
本文对DG/FV混合方法研究进展进行了简要综述，重点介绍了“混合重构”的基本思想、针对RANS方程的求解方法、隐式时间离散格式、色散耗散特性及稳定性条件、计算量理论分析等。在此基础上，给出了系列典型算例的层流和湍流计算结果，并与相应的精确解和文献结果进行了比较，以验证DG/FV方法对粘性流动问题模拟的适用性。 1 DG/FV混合方法 1.1 控制方程
在CFD中控制方程一般为Navier-Stokes方程，它的守恒形式可简写为：
其中 U 为守恒变量； F c为无粘(对流)通量； F v为粘性通量，它们的具体形式可参考文献[26]。
对于湍流计算我们使用雷诺平均N-S(RANS)方程，本文使用Spalart-Allmaras(S-A)一方程湍流模型。S-A模型的湍流涡粘性表达式为νt=vfv1，其中 v为模型因变量，其满足的微分方程为：
上式中的各种系数可参考文献[26]。本文中对式(1)和式(2)使用解耦的方式求解，并且目前仅初步使用低阶方式求解式(2)。 1.2 静态重构和动态重构
本质上说，数值格式的构造过程就是离散和重构的过程。离散即利用网格技术将计算域分解为离散网格单元；重构可以看成把解耦的信息进行耦合，把丢失的信息“还原”。不同的“还原”方式导致了不同的格式，因此重构技术对格式的构造起到至关重要的作用。在DG方法中，解一般认为是跨单元间断的多项式。解多项式系数的约束关系是时间相关的，其随“时间”同步推进计算，因此我们称之为“动态重构”(或“时间相关重构”)。解多项式的信息由初边值条件和控制方程经过Galerkin有限元方法“提炼”而来。而在k-exact FV方法中，解只有单元平均值随时间推进更新，且更新时使用一个重构的解高阶多项式分布来计算数值通量。解高阶分布的系数由相邻单元的单元平均值插值而来，不与单元平均值一起推进计算，可以看成是时间上的一种“后处理”过程。重构多项式的高阶信息来源于邻近单元，因此我们称之为“静态重构”(或“网格相关重构”)。
在本文DG/FV混合方法中，对每个单元构造其模板，并使用模板中每个单元上物理量初始较低的PDG阶多项式分布重构出更高阶的PFV阶多项式分布(PFV>PDG)，此过程称为“静态重构”；使用此PFV阶高阶分布来计算数值通量和数值积分，从而对单元中物理量的PDG阶多项式分布系数使用常规的DG方法存储和时间推进(此称为“动态重构”)，得到下一时刻的PDG阶分布系数。理论分析和数值计算表明，此静动态“混合重构”算法能够明显减少计算量和存储量。上述静态重构具体策略可以采用类似k-exact的最小二乘重构[14, 15, 25]、强插值重构[16, 17]、Gauss-Green公式重构[18, 19, 20, 21, 22]等算法。 1.3 混合重构及DG/FV空间离散
引入辅助变量 Z ，将式(1)重写为如下一阶方程组：
将计算域Ω划分成互不重合的单元集合Ωh= {Ωe}，其中Ωe表示第e个控制单元。令解空间为 V hPDG={(v)d+2: v e∈PolyPDG(Ωe)Ωe∈Ωh }，其中d为问题空间维数、PDG为DG方法的解多项式阶数、PolyPDG表示不超过PDG阶的多项式构成的线性空间；令重构空间为 V hPFV，其中PFV≥PDG为重构多项式阶数，具体定义类似 V hPDG。
单元e中，变量 U h和重构变量 W h可以分别写成如下的形式：
其中 c l和cl是多项式系数，bl( ξ )取参考单元内的正交多项式基函数，ξ 是参考单元的局部坐标。式(4)中变量 U h即为需要计算和存储的自由度，它的时间推进过程即为“动态重构”。另外可以看到，重构变量 W h可以分为两部分，它的低阶部分为了保证守恒性直接取为 U h，而高阶部分即为需要通过网格模板“静态重构”的部分。
我们数值求解式(3)即是寻求 U h∈ V hPDG，满足式(3)的弱形式。设由某种重构算法得到重构变量 W h∈ V hPFV(具体参考文献[14, 15, 16, 17, 18, 19, 20, 21, 22])，并令辅助变量 Z h∈( V hPFV)d。将式(3)乘以检验函数bj，然后在单元Ωe上使用分部积分，可得其弱形式为：
上式中为了消除 W h、 F c和 F v在边界Ωe上的二义性，分别使用了其数值通量 、 H c和 H v。 W h+和 Z h+分别表示相邻单元的重构变量和辅助变量，n 为边界单位外法向。无粘数值通量 H c可以使用常规的欧拉黎曼求解器得到[18, 19, 20, 21, 22]，本文使用Roe格式[27]。与粘性相关的数值通量、 H v如何计算是DG/FV混合方法求解N-S方程的关键内容。本文使用由Bassi和Rebay提出的BR格式[28, 29]。由于BR1格式需要求解和存储辅助变量 Z h，其计算量和存储量大，且BR1格式需要用到邻居的邻居单元信息，其不紧致，所以此处使用其改进后的BR2格式。具体如下，由式(5)再经过一次分部积分，得：
依此定义提升算子R h= Z h-
W h，则在每个单元中 R h可由式(7)方便的算出(由于使用正交基函数，不用单元体积分)，而 W h可以直接得到，因此可得 Z h= W h+ R h。进一步定义修正提升算子 r f为：
其中f是单元e的一个边界面。
粘性数值通量 H v的计算公式为：
其中ηf为稳定性因子，一般取3。
式(9)定义的关于面f的修正提升算子 r f仅依赖本单元和面f对应的相邻单元，因此式(10)中面f的粘性通量只与此两个单元相关，从而具有紧致性。
最终的半离散方程可以写为：
1.4 时间离散
随着空间精度阶的提高，使用显式时间离散时DGM及DG/FV方法对应的稳定性条件将越来越严格，时间推进步长受到明显的限制，而隐式时间离散方法的时间推进步长受计算精度阶和计算网格尺度的限制较小，因此本文使用之前发展的基于Newton/Gauss-Seidel迭代的时间隐式离散方法[22]。以如下一阶Euler后差时间离散为例：
定义关于 U n+1的非定常残差Run( U n+1)为： 若能求解Run( U n+1)=0则可得到满足一阶精度的 U n+1，但上式关于 U n+1是非线性的，不能直接求解。对其使用Newton迭代方法求解，可得： 其中l是Newton迭代指标，l=0时 U n+1,l取 U n，l→∞时 U n+1,l+1= U n+1,l= U n+1。
式(14)是关于Δ U n+1,l的大型线性方程组，但注意到RHS只与目标单元以及若干邻近的单元相关，从而式(14)左端矩阵是一个块稀疏矩阵。本文采用Gauss-Seidel迭代来求解式(14)，具体方法见文献[21, 22]。
以上求解过程可能还需要对曲边界进行修正，其它如边界条件和限制策略[30]等限于篇幅在此不再赘述。 2 数值特性分析 2.1 数值色散耗散特性
流体力学方程是复杂的非线性方程组，一般来说此方程组的数学性质，如解的存在性、唯一性、数学提法的适定性等等都还是正在研究中的问题，且很难找到一般情况下方程的解析解或精确解。所以这里我们仅针对一维线性对流方程分析DG/FV混合方法的色散、耗散特性，得到几种典型精度阶的DG/FV格式的数值特性，并和典型精度阶DG格式的数值特性进行比较。
图 1给出了几种典型格式的数值色散关系曲线，此处及下文中各格式括号中的数字表示其设计精度阶。从图 1可以看出，在高波数区，DG/FV(3)的色散特性接近DGM(2)，在低波数区，DG/FV(3)的色散特性介于DGM(2)和DGM(3)之间，表现出令人满意的色散特性。图 2给出了几种典型格式的数值耗散关系曲线，可以看到其数值耗散均小于零，说明格式是稳定的。DG/FV(2)格式的耗散较大，DG/FV(3)的耗散特性同样介于DGM(2)和DGM(3)之间，既表现出低波数区有尽量小的耗散，在高波数区又具有较大的耗散，耗散特性比较理想。
 |
| 图 1 几种典型格式的数值色散关系Fig. 1 Dispersive behavior of some schemes |
 |
| 图 2 几种典型格式的数值耗散关系Fig. 2 Dissipative behavior of some schemes |
图 3和图 4分别为DG/FV(3)格式的非物理模态的数值色散和耗散曲线。从图 3可以看出，非物理模态波的传播方向与物理模态波相反，不过从图 4可以看出非物理模态波数值耗散很大，在短距离、短时间内就会被耗散掉，这与Hu和Atkins[31]研究所得DGM的数值色散耗散特性是类似的。
 |
| 图 3 DG/FV(3)格式的非物理模态数值色散关系Fig. 3 Dispersive of DG/FV(3) in spurious mode |
 |
| 图 4 DG/FV(3)格式的非物理模态的数值耗散关系Fig. 4 Dissipative of DG/FV(3) in spurious mode |
使用Fourier方法对模型方程进行稳定性分析，得到了一些典型阶数DG和DG/FV格式的稳定性条件(见表 1)，其中时间离散使用1-5阶的显式RK方法。由表 1可以看出DG/FV格式的CFL数与低 一阶DG格式的CFL数相当，比同阶DG格式的CFL数有大幅提高，说明了此DG/FV混合格式在时间推进步长方面具有显著优势。
| 显式RK | 1阶 | 2阶 | 3阶 | 4阶 | 5阶 |
| DGM(1) | 1.00 | 1.000 | 1.2564 | 1.3927 | 1.6085 |
| DG/FV(2) | × | 1.000 | 1.1758 | 1.3846 | 1.5796 |
| DGM(2) | × | 0.333 | 0.4096 | 0.4642 | 0.5347 |
| DG/FV(3) | × | 0.100 | 0.4188 | 0.4642 | 0.5361 |
| DGM(3) | × | × | 0.2122 | 0.2352 | 0.2717 |
文献[32]在对数值通量的计算量作了合理的假设后给出了几种精度DG格式的计算量的理论分析结果。本文参考其方法，理论分析了常用的3-4阶精度DG和DG/FV混合格式的计算量，其中二维标量方程三角形网格下所 得结果见表 2。从中我们可以看出：1)无论是DG还是DG/FV格式，单元高斯积分的计算量都占总计算量很大比重(40%~50%)，这说明了减少单元积分点数目的重要性； 2)边界高斯积分的计算量所占比重随着方法精度阶数提高逐渐减小；3)DG/FV格式中重构运算计算量占总计算量约15%，数值计算也表明重构计算量比重较小；4)DG/FV格式计算量只有同阶DG格式计算量的50%左右。
| 积分点上 物理量 | 积分点 上通量 | 单元 积分 | 边界 积分 | 重构 | 总 计算量 | 相比同阶 DGM | |
| DGM(2) | 45 | 27 | 108 | 54 | 0 | 234 | — |
| DG/FV(3) | 110 | 31 | 144 | 54 | 54 | 393 | 49% |
| DGM(3) | 165 | 46.5 | 432 | 162 | 0 | 805.5 | — |
| DG/FV(4) | 304 | 50.5 | 504 | 162 | 144 | 1164.5 | 50% |
| DGM(4) | 456 | 78 | 1440 | 360 | 0 | 2334 | — |
表 3给出了本文使用的几种DG和DG/FV格式所用的单元高斯积分和边界积分节点数。从中可以看出，相比同阶精度的DG格式，DG/FV混合格式所用积分点数大幅减少。这正是DG/FV格式较同阶精度DG格式计算量小的主要原因。
| 分布多 项式阶 | 重构多 项式阶 | 单元积 分点数 | 边积分 点数 | |
| DGM(2) | 1 | — | 3 | 2 |
| DG/FV(3) | 1 | 2 | 4 | 2 |
| DGM(3) | 2 | — | 6 | 3 |
| DG/FV(4) | 2 | 3 | 7 | 3 |
| DGM(4) | 3 | — | 12 | 4 |
表 4给出了二维标量方程DG和DG/FV格式计算效率的数值计算结果。与三阶DGM(3)相比，同阶精度的DG/FV(3)混合方法计算量节省了约40%~50%，存储量也节省了约30%~40%。
| 网格 | CPU(s) | 相比同阶 DGM | 内存(M) | |
| DGM(2) | 10×10×2 | 1.65×10-3 | — | 0.31 |
| 20×20×2 | 6.52×10-3 | — | 0.80 | |
| 40×40×2 | 2.80×10-2 | — | 2.65 | |
| 80×80×2 | 1.13×10-1 | — | 8.16 | |
| DG/FV(3) | 10×10×2 | 2.57×10-3 | 57% | 0.34 |
| 20×20×2 | 9.98×10-3 | 58% | 0.95 | |
| 40×40×2 | 4.07×10-2 | 56% | 3.21 | |
| 80×80×2 | 1.58×10-1 | 54% | 11.1 | |
| DGM(3) | 10×10×2 | 4.50×10-3 | — | 0.43 |
| 20×20×2 | 1.73×10-2 | — | 1.30 | |
| 40×40×2 | 7.22×10-2 | — | 4.50 | |
| 80×80×2 | 2.95×10-1 | — | 17.1 |
为了验证DG/FV混合方法的数值精度阶，首先计算了库埃特(Couette)流动问题。此问题描述的是在两个平板y=0和y=H中，由上平板匀速运动所导致的槽道流动。取上平板速度U=1，温度T1=0.85；下平板温度T0=0.8；粘性系数μ=0.01为常数；槽道高度H=2。计算区域取为[0,4]×[0,2]，左右边界使用周期边界条件，上下边界使用等温壁条件。计算使用四套网格，粗网格为120个三角形单元，其余网格由粗网格逐次加倍而得。此问题的精确解如下[8]:
表 5给出了计算所得密度与其精确解的误差的L2模及数值精度阶。可以看出各阶DG/FV混合格式的数值精度阶数均达到了设计精度，同时DG/FV格式的计算误差位于低一阶DGM和同阶DGM结果之间，且随着网格加密更接近于同阶DGM结果。
| 单元数 | L2 Error | Order | L2 Error | Order | L2 Error | Order | L2 Error | Order |
| DG/FV(2) | DGM(2) | DG/FV(3) | DGM(3) | |||||
| 120 | 6.30×10-3 | — | 3.01×10-3 | — | 1.89×10-3 | — | 5.56×10-4 | — |
| 480 | 1.50×10-3 | 2.07 | 7.51×10-4 | 2.00 | 2.59×10-4 | 2.87 | 7.23×10-5 | 2.94 |
| 1920 | 3.51×10-4 | 2.10 | 1.76×10-4 | 2.09 | 2.96×10-5 | 3.13 | 8.53×10-6 | 3.08 |
| 7680 | 7.23×10-5 | 2.28 | 4.21×10-5 | 2.06 | 3.44×10-6 | 3.10 | 9.84×10-7 | 3.12 |
| DG/FV(4) | DGM(4) | DG/FV(5) | DGM(5) | |||||
| 120 | 3.39×10-4 | — | 9.89×10-5 | — | 6.42×10-5 | — | 3.98×10-5 | — |
| 480 | 1.73×10-5 | 4.29 | 7.28×10-6 | 3.76 | 2.57×10-6 | 4.64 | 1.98×10-6 | 4.33 |
| 1920 | 9.06×10-7 | 4.26 | 4.40×10-7 | 4.05 | 1.11×10-7 | 4.53 | 7.53×10-8 | 4.72 |
| 7680 | 4.28×10-8 | 4.40 | 2.67×10-8 | 4.04 | 3.58×10-9 | 4.95 | 3.01×10-9 | 4.64 |
方腔流动是一个层流计算的经典算例，本文采用的是边长为1的正方形空腔，计算使用混合网格，其中包含1440个四边形单元和1110个三角形单元，如图 5所示。计算条件是Ma=0.1，Re=1×104。
 |
| 图 5 计算使用混合网格，单元数为2550Fig. 5 Hybrid grid with 2550 cells |
计算所得流线和文献[33]结果如图 6所示，其中文献结果采用二阶精度有限差分方法和257×257的均匀结构网格。从流线图看，在如此稀疏的网格上，本文4阶DG/FV(4)格式和DGM(4)格式计算结果与文献都符合的很好，方腔内一个主涡在中心位置，在三个角附近分布了三个较小的二次涡。本文计算比较准确地分辨出了三个二次涡以及右下角的小涡。图 7给出了计算域内x=0.5直线位置的速度u剖面分布和y=0.5直线上的速度v剖面分布，并与文献结果进行了对比。可以看出，本文DG/FV(4)格式和DGM(4)格式计算结果均与文献结果符合的很好。从图 7还可以看出，4阶DG/FV(4)格式计算结果比3阶DGM(3)格式计算结果有很大改进。
 |
| 图 7 几中格式计算所得速度剖面Fig. 7 Velocity profile of several schemes 7 |
图 8给出了几种典型格式的显式和隐式时间离散的收敛历程比较，可以看出各阶格式采用隐式时间离散时计算CPU时间都大大减少，收敛性能提高约1~2个量级。另外可以看出，尽管迭代步数几乎相同，相比同阶精度DGM，DG/FV混合格式计算CPU时间显著减少，计算效率明显提高。
 |
| 图 8 几种典型格式残差收敛曲线Fig. 8 Convergence history of several schemes |
对于剪切层算例，计算条件为Ma1=0.5和Ma2=0.25，Re=ρ1U1δ/μ=500。计算域为[0,800] ×[-100,100]，使用25,272个单元的混合网格。计算初始条件和边界条件参考Colonius等的文献[34]，计算最终时间为t=1357.28=68Tf。图 9给出了两种格式计算所得涡量等值线云图并与文献[34]结果进行了比较。可以看出，本文4阶DG/FV(4)和5阶DG/FV(5)计算结果很一致并和文献结果类似。
 |
| 图 9 剪切层算例涡量等值线云图与文献结果比较Fig. 9 Vorticity contours at time t=68 Tf |
对于低速平板湍流算例，计算条件为来流Ma∞=0.2，Re=1.03×107，计算网格为76×61的结构网格，x方向平板端点附近第一层网格尺度为2×10-3，y方向第一层网格尺度为2×10-6。湍流模型为前述的SA模型。图 10给出了湍流平板某个站位下速度剖面分布的三阶DG/FV(3)格式计算结果，并和二阶FVM结果进行了比较，其中二阶FVM计算使用网格单元数量4倍于三阶格式计算所用网格。图 10可以看出，较粗网格下的三阶计算结果与密网格下的二阶FV计算结果相符良好。
 |
| 图 10 湍流平板计算速度剖面分布Fig. 10 Velocity profile of turbulent flat plate |
对于NACA0012翼型绕流算例，计算条件为来流Ma∞=0.7，α=1.49°，Re=9×106。湍流模型仍是前述的SA模型。计算使用网格分布为前缘0.0015，尾缘0.003，法向4.75×10-6，远场15倍弦长，网格量为265×132，如图 11所示。计算所得升力、阻力系数同文献[35]结果比较见表 6，从表中可以看出，3阶DG/FV(3)计算结果介于DGM(2)和DGM(3)之间，和文献结果也比较符合。
 |
| 图 11 NACA0012翼型计算所用网格Fig. 11 Grid of NACA 0012 airfoil |
 |
| 图 12 NACA0012翼型压力等值线分布(DG/FV(3))Fig. 12 Pressure contours of NACA0012 airfoil (DG/FV(3)) |
 |
| 图 13 NACA0012翼型计算所得湍流粘性系数分布(DG/FV(3))Fig. 13 Turbulent viscosity contours of NACA0012 airfoil (DG/FV(3)) |
本文对基于非结构/混合网格的高阶精度DG/FV混合格式进行了简要的综述，重点介绍了 “混合重构”基本思想、与BR2方法结合的RANS方程求解方法、Newton/Gauss-Seidel耦合隐式计算格式，理论分析了其色散耗散特性及稳定性条件，并给出DG/FV混合方法与DG方法计算量的定性理论分析和数值结果，表明相比同阶精度的DG方法，DG/FV混合方法计算量减少约40%。随后，利用多个典型算例考察了DG/FV的计算精度和计算效率，并与相应的DG方法进行了对比。对有解析解的Couette流动问题计算表明，本文方法达到了设计精度；方腔流动、剪切层算例和低速平板湍流、NACA0012翼型绕流等经典算例结果表明，随着精度的提高，在较粗的计算网格上亦能得到高精度的计算结果，其计算结果的精度与同阶的DG方法相当，但是计算效率大幅提高；而隐式方法进一步提高了本文方法的计算效率。由此可见高阶精度DG/FV混合格式将具有良好的工程应用前景。
下一步我们将针对三维复杂外形的湍流数值模拟进行研究，重点研究曲边界修正方法及高精度边界条件的实现、针对高速流动的限制器和间断侦测方法或人工粘性方法、基于GMRES的隐式计算方法、三维大规模分区并行计算技术、湍流模型的紧耦合高精度计算方法、基于动网格的高精度非定常计算方法及几何守恒问题、与DES和LES方法的有机结合等，期待早日实现在复杂工程问题中的成功应用。
| [1] | BAKER T J. Mesh generation: art or science[J]. Progress in Aerospace Sciences, 2005, 41: 29-63. |
| [2] | WANG Z J. High-order methods for the Euler and Navier-Stokes equations on unstructured grids[J]. Progress in Aerospace Sciences, 2007, 43: 1-41. |
| [3] | BARTH T J, FREDERICKSON P O. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction[R]. AIAA Paper 90-0013. |
| [4] | OLLIVIER-GOOCH C, NEJAT A, MICHALAK K. Obtaining and verifying high-order unstructured finite volume solutions to the Euler equations[J]. AIAA Journal, 2009, 47(9): 2105-2120. |
| [5] | COCKBURN B, SHU C W. The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems[J]. J. Comput. Phys., 1998, 141: 199-224. |
| [6] | HE L X, ZHANG L P, ZHANG H X. Discontinuous Galerkin finite element method on 3D arbitrary elements[J]. ACTA Aerodynamica Sinica, 2007, 25(2): 157-162. (in Chinese) 贺立新, 张来平, 张涵信. 任意单元间断Galerkin 有限元计算方法研究[J]. 空气动力学学报, 2007, 25(2): 157-162. |
| [7] | SUN Y, WANG Z J, LIU Y. Spectral (finite) volume method for conservation laws on unstructured grids VI: extension to viscous flow[J]. J. Comput. Phys., 2006, 215: 41-58. |
| [8] | SUN Y Z, WANG Z J, LIU Y. High-order multidomain spectral difference method for the Navier-Stokes Equations on unstructured hexahedral grids[J]. Commun. Comput. Phys., 2007, 2(2): 310-333. |
| [9] | WANG Z J, GAO H, HAGA T. A unifying discontinuous formulation for hybrid meshes[M]. In: Adaptive High-Order Methods in Computational Fluid Dynamics. Edited by WANG Z J. Singapore: World Scientific Publishing, 2011. |
| [10] | EKATERINARIS J A. High-order accurate, low numerical diffusion methods for aerodynamics[J]. Prog. Aero. Sci., 2005,41:192-300. |
| [11] | ZHANG L P, HE L X, LIU W, et al. Reviews of high-order methods on unstructured and hybrid grid[J]. Advances in Mechanics, 2013, 43(2): 202-236.(in Chinese) 张来平, 贺立新, 刘伟, 等. 基于非结构/混合网格的高阶精度格式研究进展[J]. 力学进展, 2013, 43(2): 202-236. |
| [12] | COCKBURN B, LUSKIN M, SHU C W, SULI E.Enhanced accuracy by post-processing for finite element methods for hyperbolic equations[J]. Mathematics of Computation, 2003, 72: 577-606. |
| [13] | RYAN J K, SHU C W, ATKINS H L.Extension of a post-processing technique for the discontinuous Galerkin method for hyperbolic equations with applications to an aeroacoustic problem[J]. SIAM Journal on Scientific Computing, 2005, 26: 821-843. |
| [14] | DUMBSER M, BALSARA D S, TORO E F. A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes[J]. J. Comput. Phys., 2008, 227: 8209-8253. |
| [15] | DUMBSER M. Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier-Stokes equations[J]. Computers and Fluids, 2010, 39: 60-76. |
| [16] | LUO H, LUO L P, NOURGALIEV R, et al. A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids[J]. J. Comput. Phys., 2010, 229: 6961-6978. |
| [17] | LUO H, LUO L P, ALI A, et al. A parallel, reconstructed discontinuous Galerkin method for the compressible flows on arbitrary grids[J]. Commun. Comput. Phys., 2011, 9(2): 363-389. |
| [18] | ZHANG L P, LIU W, HE L X, et al. A class of hybrid DG/FV methods for conservation laws I: basic formulation and one-dimensional systems[J]. J. Comput. Phys., 2012, 231: 1081-1103. |
| [19] | ZHANG L P, LIU W, HE L X, et al. A class of hybrid DG/FV methods for conservation laws II: two-dimensional cases[J]. J. Comput. Phys., 2012, 231: 1104-1120. |
| [20] | ZHANG L P, LIU W, HE L X, et al. A class of hybrid DG/FV methods for conservation laws III: two-dimensional Euler equations[J]. Commun. Comput. Phys., 2012, 12(1): 284-314. |
| [21] | ZHANG L P, LIU W, LI M, et al. A class of DG/FV hybrid schemes for conservation law IV: 2D viscous flows and implicit algorithm[J], Computers and Fluids, 2014, 97: 110-125. |
| [22] | ZHANG L P,LI M, LIU W, et al. An implicit algorithm of high-order DG method and hybrid DG/FV schemes based on Newton/Gauss-Seidel iteration for 2D inviscid flows on arbitrary grids[J]. Communication in Computational Physics, 2014(in press). |
| [23] | WANG Z J, SHI L, FU S, et al. A PNPM-CPR framework for hyperbolic conservation laws[R]. AIAA Paper 2011-3227. |
| [24] | SHI L, WANG Z J, FU S, et al. A PNPM-CPR method for Navier-Stokes equations[R]. AIAA Paper 2012-460. |
| [25] | LI M, LIU W, ZHANG L P, et al. A class of high order hybrid DG/FV methods for 2D viscous flows[J].Transactions of Nanjing University of Aeronautics and Astronautics, 2013, 30(S): 68-73. |
| [26] | PETERSON W, DAVID W Z. Three-dimensional aerodynamic computations on unstructured grids using a Newton-Krylov approach[J]. Computers and Fluids, 2008, 37(2): 107-120. |
| [27] | ROE P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J].J. Comput. Phys., 1981, 43: 357-372. |
| [28] | BASSI F, REBAY S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[J]. J. Comput. Phys., 1997, 131: 267-293. |
| [29] | BASSI F, REBAY S, MARIOTTI G, et al. A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows[A]. In: 2nd European Conference on Turbomachinery Fluid Dynamics and Thermodynamics[C]. Decuypere R, Dibelius G, editors. Technologisch Instituut, Antwerpen, Belgium, 1997: 99-108. |
| [30] | ZHANG L P, LIU W, HE L X, et al. A shock detection method and applications in DGM for hyperbolic conservation laws on unstructured grids[J]. ACTA Aerodynamica Sinica, 2011, 29(4): 401-406. (in Chinese) 张来平, 刘伟, 贺立新, 等. 一种新的间断侦测器及其在DGM中的应用[J]. 空气动力学学报, 2011, 29(4): 401-406. |
| [31] | HU F, ATKINS H. Eigensolution analysis of the discontinuous Galerkin method with nonuniform grids, I: one space dimension[J]. J. Comput. Phys., 2002, 182(2): 516-545. |
| [32] | SUN Y Z, WANG Z J. Evaluation of discontinuous Galerkin and spectral volume methods for scalar and system conservation laws on unstructured grids[J]. Int. J. Numer. Meth. Fluids, 2004, 45: 819-838. |
| [33] | GHIA U, GHIA K N, SHIN C T. High Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. J. Comput. Phys., 1982, 48: 387-411. |
| [34] | COLONIUS T, LELE S K, MOIN P. Sound generation in a mixing layer[J]. J. Fluid Mech., 1997, 330: 375-409. |
| [35] | MAKSYMIUK C M, PULLIAM T H. Viscous transonic airfoil workshop results using ARC2D[R]. AIAA Paper 87-0415. |