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等提出DG-FV混合格式^[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方法和高分辨率二阶格式作简要的介绍。

$ \frac{{\partial {\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}}({\boldsymbol{U}})}}{{\partial x}} = 0 $

(1)

1.1 三阶CPR方法

对于三阶CPR方法，第i个单元的三个求解点处的守恒变量 ${\boldsymbol{U}}_{i,k}^{}$ （ $ k = 1,2,3 $ ），可通过求解以下微分方程得到：

$ \frac{{\partial {{\boldsymbol{U}}_{i,k}}}}{{\partial t}} + {\left. {\frac{{\partial {{\boldsymbol{F}}_i}(\xi )}}{{\partial \xi }}} \right|_{\xi = {\xi _k}}} + {\left. {\frac{{\partial {{\boldsymbol{\sigma }}_i}(\xi )}}{{\partial \xi }}} \right|_{\xi = {\xi _k}}} = 0,\;\;\xi \in [ - 1,1] $

(2)

式中： ${{\boldsymbol{F}}_i}(\xi )$ 为通量多项式， $ {{\boldsymbol{\sigma }}_i}(\xi ) $ 是通量修正多项式，

$ \begin{split}& {{\boldsymbol{F}}_i}(\xi ) = \sum\limits_{l = 1}^3 {{{\boldsymbol{F}}_{i,l}}} {\tilde L_l}(\xi ) \\& {{\boldsymbol{\sigma }}_i}({\xi _l}) = [{{\tilde {\boldsymbol{F}}}}( - 1) - {{\boldsymbol{F}}_i}( - 1)]{g_{\rm{L}}}({\xi _l}) + [{{\tilde {\boldsymbol{F}}}}(1) - {{\boldsymbol{F}}_i}(1)]{g_{\rm{R}}}({\xi _l}) \\ \end{split} $

(3)

其中， $ {\tilde L_l}(\xi ) $ 是基于CPR求解点的一维Lagrange多项式， ${g_{\rm{L}}}(\xi )$ 和 ${g_{\rm{R}}}(\xi )$ 是修正函数， ${{\tilde {\boldsymbol{F}}}}$ 是数值通量，

$ \begin{split} {{\tilde {\boldsymbol{F}}}}( - 1) = {{\tilde {\boldsymbol{F}}}}({{\boldsymbol{U}}_{i - 1}}(1),{{\boldsymbol{U}}_i}( - 1)) \\ {{\tilde {\boldsymbol{F}}}}(1) = {{\tilde {\boldsymbol{F}}}}({{\boldsymbol{U}}_i}(1),{{\boldsymbol{U}}_{i + 1}}( - 1)) \\ \end{split}$

(4)

这里 ${{\boldsymbol{U}}_i}(\xi )$ 为第i个单元的变量分布函数，且有 ${{\boldsymbol{U}}_i}(\xi ) = \displaystyle\sum\limits_{k = 1}^3 {{{\boldsymbol{U}}_{i,k}}} {\tilde L_k}(\xi )。$

本文使用的求解点是高斯点： ${{\xi _1} = - {\sqrt {15} }}/{5}, {{\xi _2} = 0,}{\:\:{\xi _3} = {\sqrt {15} }/{5}}$ 。修正函数为：

$ \begin{split}& {g_{\rm{L}}}(\xi ) = {g_{{\rm{DG,L}}}} = {{{( - 1)}^3}} \cdot ({\tilde L_3}(\xi ) - {\tilde L_2}(\xi ))/2 \\& {g_{\rm{R}}}(\xi ) = {g_{{\rm{DG,R}}}} =({\tilde L_3}(\xi ) + {\tilde L_2}(\xi )) /2 \end{split}$

其中： ${\tilde L_3}(\xi ) = (5{\xi ^3} - 3\xi )/2，{\tilde{L}}_{2}(\xi ) = {\scriptstyle }(3{\xi }^{2}-1)/2$ 。

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非线性插值

基于网格点上的守恒变量 ${\boldsymbol{U}}_i^{}$ 通过插值可得到单元 $[ x_{i-\frac{1}{2}} , x_{i +\frac{1}{2}} ]$ 的左边界和右边界上的值 ${\boldsymbol{\tilde U}}_{i - \frac{1}{2}}^{\rm{R}}$ 和 ${\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{L}}$ 。下面我们只给出 ${\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{L}}$ 的求解过程， ${\boldsymbol{\tilde U}}_{i - \frac{1}{2}}^{\rm{R}}$ 求解过程类似。首先，对于三点模板 $ {S_{i + \frac{1}{2}}} = \{ {x_{i - 1}},{x_i},{x_{i + 1}}\} $ ，守恒变量 ${\boldsymbol{U}}_k^{}$ 转换为特征变量 ${\boldsymbol{Q}}_k^{}$ 的过程如下:

$ {\boldsymbol{Q}}_k^{} = \sum\limits_{m{\text{ = 1}}}^n {{{\boldsymbol{l}}_{i,m}}{{\boldsymbol{U}}_{k,m}}} ,\;\;k = i - 1,i,i + 1 $

(5)

式中：n为方程和变量数目， $ {\boldsymbol{Q}}_{k,m}^{} $ 为 $ {\boldsymbol{Q}}_k^{} $ 的第m个分量， $ {{\boldsymbol{l}}_{i,m}} $ 为第i个网格点处矩阵 $ {\boldsymbol{A}} = \partial {\boldsymbol{F}}/\partial {\boldsymbol{U}} $ 的第m个左特征向量。

模板 $ {S_{i + \frac{1}{2}}} = \{ {x_{i - 1}},{x_i},{x_{i + 1}}\} $ 可拆分成两个子模板 $ {S_{i + \frac{1}{2},1}} = \{ {x_{i - 1}},{x_i}\} $ 和 $ {S_{i + \frac{1}{2},2}} = \{ {x_i},{x_{i + 1}}\} $ 。在子模板上构造插值多项式 $ \begin{array}{*{20}{c}} {{{\boldsymbol{p}}^{(m)}}(x),}&{m = 1,2} \end{array} $ ，并分别得到在 $ x_{i + \frac{1}{2}}^{} $ 处的值：

$ \begin{gathered} {{\boldsymbol{p}}^{(1)}}({x_{i + \frac{1}{2}}}) = \dfrac{3}{2}{{\boldsymbol{Q}}_i} - \frac{1}{2}{{\boldsymbol{Q}}_{i - 1}} \\ {{\boldsymbol{p}}^{(2)}}({x_{i + \frac{1}{2}}}) = \dfrac{1}{2}{{\boldsymbol{Q}}_i} + \frac{1}{2}{{\boldsymbol{Q}}_{i + 1}} \\ \end{gathered} $

(6)

进而，基于 ${\text{\{ }}{{\boldsymbol{Q}}_{i - 1}},{{\boldsymbol{Q}}_i}{\text{, }}{{\boldsymbol{Q}}_{i + 1}}{\text{\} }}$ 可得到特征变量在 $ x_{i + \frac{1}{2}}^{} $ 处第m个变量的三阶非线性加权插值形式：

$ {\mkern 1mu} {\boldsymbol{Q}}_{i + \frac{1}{2},m}^{\rm{L}} = \sum\limits_{k = 1}^2 {\omega _{k,m}^{}{\boldsymbol{p}}_m^{(k)}({x_{i + \frac{1}{2}}})} $

(7)

式中 $ {\omega _{k,m}} = \dfrac{{{\beta _{k,m}}}}{{{\beta _{1,m}} + {\beta _{2,m}}}} $ 为优化的Z型权^[16-17]，

$ {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \;\begin{array}{*{20}{c}} {{\beta _{k,m}} = {d_k}[{C_q} + {{(\frac{{\tau _5^m}}{{\varepsilon + {\boldsymbol{I}}{{\boldsymbol{S}}_{k,m}}}})}^{{p_z}}}],}&{k = 1,2} \end{array} $

(8)

其中： $ {d_1} = {1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}，{d_2} = {3 \mathord{\left/ {\vphantom {3 4}} \right. } 4} $ 为线性权， $ \varepsilon = {10^{ - 6}} $ ， $ {{\boldsymbol{\tau }}_5} = | {{\boldsymbol{I}}{{\boldsymbol{S}}_0} - } {{\boldsymbol{I}}{{\boldsymbol{S}}_1}} | $ ， $ {C_q} $ 和 $ {p_z} $ 是参数，本文取 ${C}_{q} = 4$ ， ${p}_{z} = 2$ ， $ {\boldsymbol{I}}{{\boldsymbol{S}}_k} $ 为优化的光滑因子：

$ \begin{gathered} {\boldsymbol{I}}{{\boldsymbol{S}}_0} = \frac{1}{4}{(3{{\boldsymbol{Q}}_i} - 4{{\boldsymbol{Q}}_{i - 1}} + {{\boldsymbol{Q}}_{i - 2}})^2} + {({{\boldsymbol{Q}}_i} - 2{{\boldsymbol{Q}}_{i - 1}} + {{\boldsymbol{Q}}_{i - 2}})^2} \\ {\boldsymbol{I}}{{\boldsymbol{S}}_1} = \frac{1}{4}{(3{{\boldsymbol{Q}}_{i + 2}} - 4{{\boldsymbol{Q}}_{i + 1}} + 3{{\boldsymbol{Q}}_i})^2} + {({{\boldsymbol{Q}}_{i + 2}} - 2{{\boldsymbol{Q}}_{i + 1}} + {{\boldsymbol{Q}}_i})^2} \\ \end{gathered} $

(9)

最后，特征变量 ${\boldsymbol{Q}}_{i + \frac{1}{2}}^{\rm{L}}$ 可以转化为守恒变量 ${\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{L}}$ ， ${{\tilde U}}_{i + \frac{1}{2}}^{\rm{L}} = \displaystyle\sum\limits_{m{{ = 1}}}^n {{\boldsymbol{Q}}_{i + \frac{1}{2},m}^{\rm{L}}{{\boldsymbol{r}}_{i,m}}}$ ， ${{\boldsymbol{r}}_{i,m}}$ 是第i个网格点处矩阵 $ {\boldsymbol{A}} = \partial {\boldsymbol{F}}/\partial {\boldsymbol{U}} $ 的第m个右特征向量。

1.2.2 数值通量

数值通量 $ {{\tilde {\boldsymbol{F}}}}_{i + \frac{1}{2}}^{} $ 可以基于 ${\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{L}}$ 和 ${\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{R}}$ 求出： ${{\tilde {\boldsymbol{F}}}}_{i + \frac{1}{2}}^{} = {\boldsymbol{F}}({\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{L}},{\boldsymbol{\tilde U}}_{i + \frac{1}{2}}^{\rm{R}})$ 为黎曼通量。本文使用Lax-Friedrich通量，具体表达式如下：

$ {{\tilde {\boldsymbol{F}}}}_{i + \frac{1}{2}}^n = \frac{1}{2}[{\boldsymbol{F}}({\boldsymbol{U}}_i^n) + {\boldsymbol{F}}({\boldsymbol{U}}_{i + 1}^n) - {\alpha _n}({\boldsymbol{U}}_{i + 1}^n - {\boldsymbol{U}}_i^n)] $

(10)

其中， $ {\alpha _n} = \mathop {\max }\limits_{{{\boldsymbol{U}}^n}} \left| {{\boldsymbol{F'}}({\boldsymbol{U}})} \right| $ 。

1.2.3 通量导数

本文采用二阶差分算子计算通量导数的近似值 $ {({{{\tilde {\boldsymbol{F}}}}_x})_i} $ ： $ {({{\boldsymbol{F}}_x})_i} = \dfrac{1}{h}({{{\tilde {\boldsymbol{F}}}}_{i + \frac{1}{2}}} - {{{\tilde {\boldsymbol{F}}}}_{i - \frac{1}{2}}}) $ ，其中 $ h $ 为网格步长。

1.3 混合格式

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

1.3.1 格式切换策略

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

$\begin{split}& NI = \frac{1}{{{{[r(r + 1)]}^{1/2}}}}\cdot \\&\;\;\;\;\;\;\;\;\;\; {\Bigg[\sum\limits_{k = 0}^r {{{\Bigg(\frac{{[1/(r + 1)] - \Big[({\omega _k}/{d_k})\Big/\displaystyle\sum\limits_{k = 0}^r {({\omega _k}/{d_k})} \Big]}}{{[1/(r + 1)]}}\Bigg)^2}}} \Bigg]^{1/2}} \end{split}$

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

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

1.3.2 界面处理

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

$ {\boldsymbol{U}}({g_{i,k}}) = \sum\limits_{k = 1}^3 {\boldsymbol{U}} ({c_{i,k}}){\tilde L_k}({g_{i,k}}) $

(11)

其中 $ {\tilde L_k} $ 是基于 ${c_{i,k}}$ 点的拉格朗日基函数。

其次，我们使用CPR单元的分布函数求得 ${\boldsymbol{U}}_{i + \frac{1}{2}}^{ - ({\rm{cpr}})}$ ，并由HRFD2虚拟点及HRFD2单元值进行非线性插值计算得到 ${\boldsymbol{U}}_{i + \frac{1}{2}}^{ + ({\rm{HRFD2}})}$ ，进而构造界面处的数值通量：

$ {{{\tilde {\boldsymbol{F}}}}_{i + \frac{1}{2}}} = {{\tilde {\boldsymbol{F}}}}\left( {{\boldsymbol{U}}_{i + \frac{1}{2}}^{ - ({\rm{cpr}})},{\boldsymbol{U}}_{i + \frac{1}{2}}^{ + ({\rm{HRFD2}})}} \right) $

(12)

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

1）将模板 $S = \left\{ {{w_{i - 1,3}},\;{w_{i,1}},\;{w_{i,2}}} \right\}$ 划分为两个子模板 ${S_1} = \left\{ {{w_{i - 1,3}},\; {w_{i,1}}} \right\},$ ${S_2} = \left\{ {{w_{i,1}},\; {w_{i,2}}} \right\}；$

2）分别在 $ {S_1} $ 和 $ {S_2} $ 上重构插值多项式 $ {{\boldsymbol{U}}_1}({c_{i,1}}) $ 和 $ {{\boldsymbol{U}}_2}({c_{i,1}}) $ ：

$ \begin{gathered} {{\boldsymbol{U}}_1}({c_{i,1}}) = {\boldsymbol{U}}({w_{i - 1,3}}){{\hat L}_1}({c_{i,1}}) + {\boldsymbol{U}}({w_{i,1}}){{\hat L}_2}({c_{i,1}}) \\ {{\boldsymbol{U}}_2}({c_{i,1}}) = {\boldsymbol{U}}({w_{i,1}}){{\hat L}_1}({c_{i,1}}) + {\boldsymbol{U}}({w_{i,2}}){{\hat L}_2}({c_{i,1}}) \\ \end{gathered} $

(13)

其中 $ {\hat L_{}} $ 是基于WCNS求解点的线性拉格朗日多项式。

3）计算非线性权。线性权可由

$ \begin{split}& {d_1}{{\boldsymbol{U}}_1}({c_{i,1}}) + {d_2}{{\boldsymbol{U}}_2}({c_{i,1}}) = {\boldsymbol{U}}({w_{i - 1,3}}){{\hat L}_1}({c_{i,1}}) +\\& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{U}}({w_{i,1}}){{\hat L}_2}({c_{i,1}}) + {\boldsymbol{U}}({w_{i,2}}){{\hat L}_3}({c_{i,1}}) \end{split} $

(14)

得到，则非线性权由 ${\omega _k} = \dfrac{{{\beta _k}}}{{{\beta _1} + {\beta _2}}}$ 和公式（8）得到。

4）最后得插值结果 $ {\boldsymbol{U}}({c_{i,1}}) = {\omega _1}{{\boldsymbol{U}}_1}({c_{i,1}}) + {\omega _2}{{\boldsymbol{U}}_2}({c_{i,1}}) $ 。

2 数值算例

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

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

$ {{\boldsymbol{e}}_{{L_2}}} = \sqrt {\frac{{\displaystyle\sum\limits_{i = 1}^N {({\boldsymbol{U}}_i^h - {\boldsymbol{U}}(} {r_i}){)^2}}}{N}} $

(15)

其中 $ {\boldsymbol{U}}_i^h $ 和 $ {\boldsymbol{U}}({r_i}) $ 分别是第i个求解点处的数值解和真解，N是求解点总个数。

2.1 精度测试

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

$ \frac{\partial }{{\partial t}} \left[ \begin{array}{l} \rho \\ \rho u\\ \rho v\\ E \end{array} \right] + \frac{\partial }{{\partial x}}\left[ {\begin{array}{*{20}{c}} {\rho u}\\ {\rho {u^2} + p}\\ {\rho uv}\\ {u\left( {E + p} \right)} \end{array}} \right] + \frac{\partial }{{\partial y}}\left[ {\begin{array}{*{20}{c}} {\rho v}\\ {\rho uv}\\ {\rho {v^2} + p}\\ {v\left( {E + p} \right)} \end{array}} \right] = 0$

(16)

其中： $ p = (\gamma - 1)[E - \rho {{({u^2} + {v^2})} \mathord{/ {\vphantom {{({u^2} + {v^2})} 2}} } 2}], $ 对理想气体， $ \gamma = 1.4 $ ； $\;\rho 、u、v、p、E$ 分别代表密度、x方向速度、y方向速度、压力和总能量。

初值是均匀流 $ \{ \rho ,u,v,p\} = \{ 1,1,1,1\} $ 。然后将等熵涡添加到平均流中，等熵涡表达式为：

$ \begin{split}& (\Delta u,\;\Delta v) = {\dfrac{\varepsilon }{{2\text {π} }}{{\rm{e}}^{0.5(1 - {r^2})}}( - y,x)}\\&{\Delta T = - \dfrac{{(\gamma - 1){\varepsilon ^2}}}{{8\gamma {{\text {π}} ^2}}}{{\rm{e}}^{(1 - {r^2})}}} \end{split}$

其中： $ r = \sqrt {{x^2} + {y^2}} $ ，涡强度 $ \varepsilon = 5 $ 。计算区域为 $ [ - 10,10] \times [ - 10,10] $ ，采用周期边界条件。时间步长 $ \Delta T = 0.000\;1 $ ，计算到 $ T = 1 $ ，参数 $ {\delta _{NI}} = 0.000\;01 $ 。在这个光滑的算例中我们将 $ {\delta _{NI}} $ 取为很小的值，来避免全流场都采用CPR来计算。如果全流场采用CPR的情况下，其误差分布同样在表1中给出。本文取时间步长足够小，使得时间离散误差不会影响空间精度的测试。

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

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

2.2 激波捕捉能力测试 2.2.1 激波管问题

考虑一维欧拉方程：

$\begin{split}& \qquad\qquad\qquad{{\boldsymbol{U}}_t} + {\boldsymbol{F}}{({\boldsymbol{U}})_x} = 0 \\& {\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} \rho \\ {\rho v} \\ E \end{array}} \right],{\mkern 1mu} {\mkern 1mu} {\boldsymbol{F}}{\mkern 1mu} ({\boldsymbol{U}}) = \left[ {\begin{array}{*{20}{c}} {\rho v} \\ {\rho {v^2} + p} \\ {v(E + p)} \end{array}} \right],E = \dfrac{p}{{\gamma - 1}} + \dfrac{1}{2}\rho {v^2} \\\end{split}$

(17)

其中： $ \;\rho $ 是密度， $ v $ 是速度， $ p $ 是压力， $ E $ 是总能量， $ \gamma = 1.4 $ 。

2.2.1.1 Sod激波管问题^[19]

初值为 $(\rho ,u,p) = \left\{ {\begin{array}{*{20}{l}}{\begin{array}{*{20}{c}}{(0.125,0,0.1)}\\{(1.0,0,1.0)}\end{array}}&{\begin{array}{*{20}{c}}{ - 5\leqslant x < 0}\\{0 \leqslant x \leqslant 5}\end{array}}\end{array}} \right. $

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

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

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

初值为

$ (\rho ,u,p) = \left\{ \begin{gathered} \begin{array}{*{20}{c}} {(3.857143,2.629369,10.33333),}&{x < 4} \end{array} \\ \begin{array}{*{20}{c}} {(\begin{array}{*{20}{c}} {1.0 + 0.2\sin (5x),}&{0,}&{1.0} \end{array}),}&{\begin{array}{*{20}{c}} {}&{ x \geqslant 4} \end{array}} \end{array} \\ \end{gathered} \right. $

这一算例，采用Dirichlet边界条件。时间步长为 $ \Delta T = 0.000\;1 $ ，参数 $ {\delta _{NI}} = 0.2 $ ，计算到 $ T = 1.8 $ 。参考解由三阶WCNS在9000自由度网格上计算得到。

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

2.2.1.3 Lax激波管问题^[21]

初值为 $ (\rho ,u,p) = \left\{ \begin{gathered} \begin{array}{*{20}{l}} {(0.445,0.698,3.528),}& { - 5 \leqslant x < 0} \end{array} \\ \begin{array}{*{20}{c}} {(\begin{array}{*{20}{c}} {0.5,}& {0,}& {0.571} \end{array}),}& {\begin{array}{*{20}{c}} {}& {0 \leqslant x \leqslant 5} \end{array}} \end{array} \\ \end{gathered} \right. $

这一算例，采用基于解析解的Dirichlet边界条件。时间步长为 $ \Delta T = 0.000\;1 $ ，参数 $ {\delta _{NI}} = 0.2 $ ，计算到 $ T = 1.5 $ 。参考解由自由度为2100的三阶WCNS得到。

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

2.2.2 二维黎曼问题

计算区域 $ (x,y) \in \left[ {0,1} \right] \times \left[ {0,1} \right] $ ，初值为

$ (\rho ,u,v,p)(x,y,0) = \left\{ {\begin{array}{*{20}{l}} {\begin{array}{*{20}{c}} {(1.5,\; 0,\; 0,\; 1.5),} \qquad\qquad\quad\\ \qquad\qquad {0.8 < x < {\text{1}}.{\text{0}},{\text{0}}.{\text{8}} < y < {\text{1}}.{\text{0}}} \\ {(0.5323,\; 1.206,\; 0,\; 0.3),} \qquad\qquad\quad\\ \;\;\; \qquad\qquad{0 < x \leqslant 0.8,0.8 < y < 1.0} \\ {(0.138,\; 1.206,\; 1.206,\; 0.029),} \qquad\qquad\quad\\ \;\;\;\;\;\;\; \qquad\qquad{0 < x \leqslant 0.8,0 < y \leqslant 0.8} \\ {(0.5323,\; 0,\; 1.206,\; 0.3),} \qquad\qquad\quad\\ \,\;\;\; \qquad\qquad{0.8 < x < {\text{1}}.{\text{0,0}} < y \leqslant {\text{0}}.{\text{8}}} \end{array}} \end{array}} \right. $

计算到 ${T = 0}{\text{.8}}$ ，取参数 $ {\delta _{NI}} = 0.1 $ 。

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

2.2.3 双马赫反射问题^[22]

这个算例由Woodward提出，目的在于测试高分辨率格式的鲁棒特性。一个马赫数为10的激波设定在下边界 $x = 1/6$ 处，和x轴正方向的夹角为60^o。左边界为入流条件，右边界为出流条件。下边界 $ y = 0 $ 和 $x \leqslant 1/6$ 处给定激波波后值，在 $ y = 0 $ 和 $x > 1/6$ 给定滑移固壁边界条件。上边界 $ y = 1 $ 设定为时间相关的边界条件。初值为： $(\rho ,u,v,p) = \left\{\begin{array}{ll}\begin{array}{c}(1.4,0.0,0.0,1.0),\\ (8.0,7.145,-4.125,116.5),\end{array} & \begin{array}{c}y < \sqrt{3}\Big(x-1/6\Big)\\ y\geqslant\sqrt{3}\Big(x-1/6\Big)\end{array} \end{array}\right.$ 取时间步长为 $ \Delta T = 0.000\;1 $ ，计算到 $ T = 0.2 $ ，参数 $ {\delta _{NI}} = 0.9 $ 。

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

2.2.4 马赫数 2000射流问题

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

2.2.5 激波-旋涡干扰问题

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

左右边界分别设为流入、流出边界条件，上下边界设为滑移固壁边界条件。时间步长 $ \Delta T = 0.000\;1 $ ，计算停止时间 $ T = 0.7 $ ，参数 $ {\delta _{NI}} = 0.9 $ 。

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

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

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

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

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

3 结　论

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

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

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

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

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