CPR 方法与高分辨率二阶格式的混合算法

郭 嘉，朱华君，石国权，燕振国，*，宋松和，唐玲艳

(1. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000；2. 国防科技大学 理学院 数学系，长沙 410073)

0 引 言

高阶精度格式相比于低阶格式在精细结构模拟中更具优势，随着工程问题对于精细结构模拟要求的逐步提升，高阶精度格式成为计算流体力学中的一个重点研究方向[1-3]。在高阶精度格式中，重构修正方法（correction procedure via reconstruction, CPR）具有格式紧致、高效、可应用于复杂非结构网格的优点。它最初由Huynh[4]提出，并由王志坚等[5]推广至非结构网格上， 现被较多应用到科学计算和工程预测之中。但CPR 在捕捉激波时容易产生虚假的数值振荡[1]，故CPR 的激波捕捉技术成为研究焦点[1-2,4-5]。

舒其望等[6]在CPR 方法中引入WENO 限制器，保持了光滑区域上的高阶精度特性，并在间断附近有效抑制了数值振荡。李万爱等[7]通过提出一个p加权的过程，在间断Galerkin 方法（discontinuous Galerkin,DG）上构造了新的光滑因子，提高了激波捕捉性能，然而间断附近的振荡现象未被完全消除。具有子单元分辨率的人工黏性法[8]也被提出，但是它依赖于参数的人工选择。

一个可行的解决办法是构造有限元方法（如CPR 或DG）和有限差分法（ finite difference method,FD）/有限体积法（finite volume method, FV）的混合格式，使得在间断附近具有更好的局部激波捕捉能力。这种方法的思路是，在光滑区域使用高阶有限元方法保持紧致和高分辨率的性质，在间断附近采用FD 或FV 提供鲁棒的激波捕捉能力。Dumbser 等提出DGFV 混合格式[9]并使用二阶FV 来捕捉激波。刘铁钢等[10]提出具有良好几何灵活性和高计算效率的Runge-Kutta 间 断 Galerkin 方 法（ the Runge-Kutta discontinuous Galerkin method, RKDG）和加权本质无振荡格式（weighted essentially non-oscillatory schemes,WENO）的区域混合格式。这些格式避免了高阶有限元方法直接处理激波，而是使用具有良好激波捕捉能力的格式进行激波捕捉，从而使得混合格式具有良好的激波捕捉特性。

朱华君等[11]构造了基于区域混合思路的混合WCNS-CPR 格式，解决了混合格式的空间精度和几何守恒律等基础问题，该格式计算效率较高且易于应用到复杂曲网格中。最近，基于单元动态混合思想，郭嘉和朱华君等[12]构造了WCNS-CPR 格式的激波捕捉技术。本文，我们在混合WCNS-CPR 格式基础上，针对高阶CPR 方法，通过引入激波侦测器和基于高阶WCNS 插值的二阶格式来模拟带有激波的流动问题，进一步提高混合算法的计算效率并增强混合算法的激波捕捉能力。

1 数值方法

以一维守恒律方程（1）为例，对三阶CPR 方法和高分辨率二阶格式作简要的介绍。

1.1 三阶CPR 方法

1.2 高分辨率二阶格式

WCNS 格式[13-14]的构造过程主要包含三部分：非线性插值、数值通量以及通量导数的求解。文献[15]中指出，基于对称守恒型网格导数算法（ symmetrical conservative metric method， SCMM） 的WCNS 可分解为多个二阶格式。本文采用了基于高阶WCNS 插值的二阶格式，此格式为三阶WCNS 的分解格式之一，即非线性插值为高阶WCNS 插值，而通量导数使用二阶差分算子的格式，这里简称为高分辨率二阶格式（ high-resolution second order finite difference scheme，HRFD2）。

1.2.1 三阶WCNS 非线性插值

1.3 混合格式

本文采用与高阶WCNS-CPR 混合算法[12]中相同的混合策略，关键在于将WCNS 格式替换为基于WCNS 插值的二阶格式HRFD2，进一步提高计算效率并增强混合算法的激波捕捉能力。下面对混合格式中的格式切换策略和界面处理方法进行介绍。

1.3.1 格式切换策略

混合格式在不光滑区域使用具有激波捕捉能力的HRFD2，而在光滑区域使用具有较高计算效率和较高分辨率的CPR 格式。因此需要激波侦测器来侦测激波所在位置。本文使用NI 激波侦测器[18]来判断激波所在单元，并将其标记为问题单元（Troubled cell）。NI的公式为：

其中：r+1 是非线性加权的候选模板个数，对于HRFD2，r=1。为了保证激波不会移动到光滑单元（Smooth cell），在问题单元附近引入缓冲单元（Buffer cell），如图1 所示。

图1 三种单元类型示意图Fig. 1 Illustration of three different types of cells

所采用的格式切换策略如下：当激波侦测器侦测到激波时，在激波附近的单元格式切换为HRFD2。当激波离开后，这些单元的格式由HRFD2 切换为CPR。综上，即根据激波侦测器的值，将所有单元分为不光滑单元（问题单元和缓冲单元）和光滑单元，在不光滑单元采用HRFD2，而在光滑单元采用CPR。

1.3.2 界面处理

为了避免两种格式间出现错误的信息交换，我们需要仔细研究界面处理，主要的困难在于HRFD2-CPR 界面的处理。下面以一维情况为例，给出界面处理方案。首先，如图2 所示，第i+1 单元的HRFD2 需要其相邻的第i单元（CPR 单元）的虚拟点信息{gi,k,k=1,2,3}。HRFD2 单元在虚拟点的函数值是通过使用CPR 的信息 {ci,k,k=1,2,3}插值得到的，即：

图2 界面处理Fig. 2 Schematic of interface treatment

此外，我们也要考虑格式切换时的信息传递。当第i个单元在第tn时刻是光滑单元，而在下一时刻为不光滑单元时，HRFD2 的求解点值可以通过CPR 的拉格朗日多项式插值得到。而对于相反的情形，即第tn时刻为不光滑单元，下一时刻切换为光滑单元时，信息应该由HRFD2 传递给CPR。这个过程涉及到了跨界面非线性插值，我们以第i个单元的第1 个求解点为例：

2 数值算例

本节中我们将会测试CPR 与HRFD2 的混合格式（简记为CPR-HRFD2）的空间精度、激波捕捉能力和计算效率。同时会进行CPR-HRFD2 与HRFD2 的对比。

我们首先在二维等熵涡问题中测试CPRHRFD2 的精度，简单对比了HRFD2 和CPR 的计算效率，并在一维激波管问题中验证格式的激波捕捉性能，在二维黎曼问题和双马赫反射问题测试激波捕捉能力以及对比计算效率，随后用高马赫数问题来测试格式的强激波捕捉能力，最后在激波-旋涡干扰问题中研究对于复杂流动问题的模拟。时间推进格式为显式三阶龙格库塔格式。L∞和L2误差分别用来度量误差的局部和全局值。L2误差计算公式为：

2.1 精度测试

我们在二维等熵涡问题中对构造的混合格式进行精度测试。下面考虑二维欧拉方程：

图3、图4、图5、图6 分别展示了CPR-HRFD2 的不同单元类型分布、x方向NI值、y方向NI值以及误差分布图。从图4、图5 可看出，当NI值大于 δNI时，单元会被判定为问题单元，单元分布图中这个单元显示为红色，表示此单元由HRFD2 计算。缓冲单元也显示为红色，由HRFD2 计算。反之，当NI值小于δNI时，单元会被判定为光滑单元，单元分布图中这个单元则显示为蓝色，此单元由CPR 计算。从表1可知，CPR-HRFD2 具有二阶精度。

表1 二维精度测试Table 1 2D accuracy test

图3 不同单元分布(自由度为120)Fig. 3 Distribution of different cells, DOFs = 120

图4 x 方向NI 侦测结果分布(自由度为120)Fig. 4 NI detection result in x-direction, DOFs = 120

图5 y 方向NI 侦测结果分布(自由度为120)Fig. 5 NI detection result in y-direction, DOFs = 120

图6 误差分布(自由度为120)Fig. 6 Error distribution, DOFs = 120

此外，针对二维等熵涡问题，进行了HRFD2 与CPR 的计算时间对比测试，取相同的自由度，时间步长ΔT=0.000 1， 计算到T=1。由表1 误差分布和表2计算时间可知，在相同计算误差下CPR 计算效率比HRFD2 高。

表2 二维等熵涡问题HRFD2 与CPR 计算时间对比Table 2 CPU time of HRFD2 and CPR in isentropic vortex problem

2.2 激波捕捉能力测试

2.2.1 激波管问题

考虑一维欧拉方程：

这一算例，采用基于解析解的Dirichlet 边界条件。δNI=0.2 ， 时间步长为 ΔT=0.000 1， 计算到T=3.0。参考解由三阶WCNS 在2100 自由度网格上计算得到。

由图7 可知，激波和接触间断都可以被高分辨率地捕捉到。在接触间断和激波附近，没有明显的振荡，这展现了CPR-HRFD2 良好的激波捕捉能力。此外，大部分单元由CPR 计算，保证具有更高的计算效率。

图7 Sod 激波管问题，自由度为120Fig. 7 Sod shock tube problem, DOFs = 120

2.2.1.2 Shu-Osher 激波管问题[20]

初值为这一算例，采用Dirichlet 边界条件。时间步长为ΔT=0.000 1， 参数 δNI=0.2， 计算到T=1.8。参考解由三阶WCNS 在9 000 自由度网格上计算得到。

由图8 可知，在激波侦测器判别的充分光滑区域内使用了CPR。此外，CPR-HRFD2 可以很好地捕捉激波和间断。密度的数值解贴近参考解，展现了高分辨率特点。

图8 Shu-Osher 激波管问题，自由度为600Fig. 8 Shu-Osher shock tube problem, DOFs = 600

2.2.1.3 Lax 激波管问题[21]

由图9 可知，在x≈2.3和x≈3.8处，激波和间断被高分辨率地捕捉到，整个流场没有观察到明显振荡。

图9 Lax 激波管问题，自由度为300Fig. 9 Lax shock tube problem, DOFs = 300

2.2.2 二维黎曼问题

计算区域(x,y)∈[0,1]×[0,1]，初值为

计算到T=0.8， 取参数 δNI=0.1。

图10 为HRFD2 结果。由图11、图12 和图13 可知，在光滑区域，流场由CPR 精确计算，且无明显振荡。此外，小尺度结构也被CPR-HRFD2 高分辨率地捕捉到了，结果图（图11）和HRFD2 的结果图（图10）相近。大部分区域由CPR 计算，减少了计算代价。

图10 HRFD2 密度计算结果，自由度为960×960Fig. 10 Density result of HRFD2, DOFs = 960×960

图11 CPR-HRFD2 的计算密度结果，自由度为960×960Fig. 11 Density contour calculated by CPR-HRFD2 method，DOFs = 960×960

图12 CPR-HRFD2 的x 方向NI 侦测Fig. 12 NI detection contour by CPR-HRFD2 method in x direction

图13 CPR-HRFD2 的y 方向NI 侦测Fig. 13 NI detection contour by CPR-HRFD2 method in y direction

2.2.3 双马赫反射问题[22]

由图14、图15 和图16 知，CPR-HRFD2 的结果图和HRFD2 以及三阶WCNS 的结果图基本相同，此外，重要的小尺度和大尺度结构都被CPR-HRFD2 捕捉到了，其他区域也由CPR 精确计算。由图17 和图18知，整个区域CPR 单元个数远超过HRFD2 个数。另外，由表3 知，CPR-HRFD2 计算效率高于HRFD2（计算时间为18 480 s）和三阶WCNS（计算时间为18 840 s），当 δNI=0.9时，相比于HRFD2 提高了34.7%的计算效率，相比于三阶WCNS 提高了36.0%的计算效率。

表3 双马赫反射问题的CPU 时间对比（自由度为1 200×300）Table 3 CPU time comparison of double Mach problem, DOFs = 1 200×300

图14 三阶WCNS 格式密度计算结果，自由度为1 200×300Fig. 14 Density contour by WCNS3 method, DOFs = 1 200×300

图15 HRFD2 密度计算结果，自由度为1 200×300Fig. 15 Density contour by HRFD2 method, DOFs = 1 200×300

图16 CPR-HRFD2 密度计算结果，自由度为1 200×300Fig. 16 Density contour by CPR-HRFD2 method,DOFs = 1 200×300

图17 CPR-HRFD2 的x 方向NI 侦测Fig. 17 NI detection contour by CPR-HRFD2 method in x direction

图18 CPR-HRFD2 的y 方向NI 侦测Fig. 18 NI detection contour by CPR-HRFD2 method in y direction

2.2.4 马赫数 2000 射流问题

我们研究马赫数2 000 射流问题，计算域为[0,1]×[-0.25,0.25]， 参 数 γ=5/3。 初 值 为(ρ,u,v,p)=(0.5,0,0,0.412 7)。右边界、上边界和下边界是出流条件，左边 界 的 入 流 条 件 为：若y∈[-0.05,0.05]，(ρ,u,v,p)=(5,800,0,0.412 7)； 否则 (ρ,u,v,p)=(5,0,0,0.412 7)。对于这个问题，喷流速度为800，相对于射流气体中的声速来说，大约是马赫数2 100。我们用自由度1 710×855，计算到T=0.001， 时间步长 ΔT=5.0×10-4Δx,画出密度和压力的对数结果图。文献[23]中使用带保正限制器的五阶WCNS 格式对这个问题进行了模拟。如图19 和图20 所示，数值模拟结果与文献[23]、文献[24]相符，展现了CPR-HRFD2 格式的强激波捕捉能力。此外由图21 和图22 可知，大部分区域使用CPR 计算，保持了高计算效率。

图19 CPR-HRFD2 密度计算结果，自由度1 710×855Fig. 19 Density contour by CPR-HRFD2 method, DOFs = 1 710×855

图20 CPR-HRFD2 压力计算结果，自由度1 710×855Fig. 20 Pressure contour by CPR-HRFD2 method,DOFs = 1 710×855

图21 CPR-HRFD2 的x 方向NI 侦测Fig. 21 NI detection contour by CPR-HRFD2 method in x direction

图22 CPR-HRFD2 的y 方向NI 侦测Fig. 22 NI detection contour by CPR-HRFD2 method in y direction

2.2.5 激波-旋涡干扰问题

在最后，我们测试了激波-旋涡干扰问题[25]。该问题包含了一个运动的旋涡和一道静止正激波，可以发展出具有平滑特征和间断的复杂流动结构。如图23 所示，该问题的计算域 (x,y)∈[0,2]×[0,1]。初始条 件 由 位 于 (0.5,y)的 静 态 激 波Mas和 中 心 位 于(xc,yc)=(0.25,0.5)的 复合旋涡Mav给出，初始条件的具体设置参考文献[25]。

图23 激波-旋涡干扰的初始条件和边界条件Fig. 23 Initial and boundary conditions of shock-vortex

左右边界分别设为流入、流出边界条件，上下边界设为滑移固壁边界条件。时间步长 ΔT=0.000 1，计算停止时间T=0.7， 参数 δNI=0.9。

如图24、图25 所示，CPR-HRFD2 在保持激波捕捉能力的同时，对平滑旋涡产生的小尺度旋涡结构也有很好的捕捉能力，计算结果与文献[25]中的计算结果相符合。由图26 和图27 可知，大部分区域使用CPR 计算，保持了较高的计算效率。

图24 HRFD2 密度计算结果，自由度960×480Fig. 24 Density contour by HRFD2 method, DOFs = 960×480

图25 CPR-HRFD2 密度计算结果，自由度960×480Fig. 25 Density contour by CPR-HRFD2 method,DOFs = 960×480

图26 CPR-HRFD2 的x 方向NI 侦测Fig. 26 NI detection contour by CPR-HRFD2method in x direction

激波-旋涡干扰问题包含复杂的流动结构。按照如图28 所示的切线位置，可以对不同结构特征进行定量分析[26]。不同参考线位置对应着不同的主要特征，见表4。

图28 参考线位置Fig. 28 Illustration of the position of the reference lines

表4 参考线Table 4 Reference lines

为对比CPR-HRFD2 和HRFD2 的激波捕捉能力和分辨率特性，在此选取①②③参考线进行分析。

从图29 中①可看出，CPR-HRFD2 和HRFD2 都有较强的激波捕捉能力，激波附近和波后区域流场光滑，无明显的振荡。从图29 中①③则可以明显看出，对于旋涡的捕捉，HRFD2 的分辨率不如CPR-HRFD2。从图26 和图27 侦测的结果来看，CPR-HRFD2 在旋涡处也有使用HRFD2 计算，但是与完全使用HRFD2相比，旋涡结构附近更多地使用CPR 求解，且相比于HRFD2，CPR 更高的精度避免了更大的耗散。如果增加 δNI，减少或避免侦测到旋涡结构，分辨率特性差异将会更加明显。

图29 HRFD2 和CPR-HRFD2 的参考线结果Fig. 29 The reference lines calculated by HRFD2 and CPR-HRFD2

另外，文献[26]中采用基于无参数梯度限制器的三阶和四阶CPR 方法对该问题进行了模拟，数值结果中观察到数值振荡，在定激波的右侧产生许多高频数值振荡。而在本文的数值结果中不管是密度图还是切线图都没有看到明显的数值振荡现象，说明了本文对高阶CPR 方法所发展的基于混合算法的激波捕捉格式对激波旋涡问题具有更好的模拟能力。

3 结 论

本文发展了三阶CPR 与HRFD2 的混合算法，使用基于非线性权的激波侦测器来判断问题单元和缓冲单元，提出了格式切换和界面处理策略。在问题单元和缓冲单元内，HRFD2 用来捕捉激波，而在光滑单元内使用CPR，用于提高计算效率。文中进行了精度测试，激波捕捉能力测试，以及计算效率测试，得到了如下结论：

1）二维的精度测试结果表明CPR-HRFD2 保持了精度；

2）通过对一维激波管问题，二维黎曼问题，双马赫反射问题，高马赫数射流问题和激波-旋涡干扰问题进行数值模拟，取得了比较好的结果，表明该混合格式具有强激波捕捉能力，可以应用于高超声速的高效数值模拟；

3）CPR-HRFD2 是一种高分辨率激波捕捉格式，同时具有较高计算效率。相比于HRFD2，CPR-HRFD2计算量更小，计算效率更高。

本文的构造思想可以直接推广到更高阶的CPR方法与基于更高阶WCNS 插值的二阶格式的混合算法中。