DOI: 10.12677/aam.2024.1312522, PDF, HTML, XML,  被引量    科研立项经费支持
作者: 贺 婧, 高 巍^*：内蒙古大学数学科学学院，内蒙古 呼和浩特
关键词: 双曲守恒律方程；有限体积法；Hermite插值；杂交因子；杂交格式；Hyperbolic Conservation Law Equations； Finite Volume Method； Hermite Interpolation； Hybrid Indicator； Hybrid Scheme

文章引用：贺婧, 高巍. 一维双曲守恒律方程的保极值无振荡杂交有限体积格式[J]. 应用数学进展, 2024, 13(12): 5406-5419. https://doi.org/10.12677/aam.2024.1312522

1. 引言

双曲守恒律方程是一类重要的偏微分方程，其非线性对流特征导致光滑初值问题也会产生间断解，它常常导致数值计算误差过大而无法有效数值逼近物理解，从而计算失败。因此，研究者们一直在尝试构造无振荡、高分辨的高效数值方法，希望数值解能在消除非物理振荡的同时，在光滑区域具有高阶精度，在间断区域保持较高的分辨率。

一阶线性格式，如FOU (First-Order Upwind) [1]，格式构造形式相对比较简单，可以保持计算稳定，不会产生虚假数值振荡，但在间断区域附近其数值解会产生较大的数值耗散。减小数值耗散带来的模拟精度损失，沿着相继提出二阶及以上的高阶线性格式(如二阶迎风格式SOU (Second-Order Upwind) [2]、Lax-Wendroff格式[3]、QUICK (Quadratic Upstream Interpolation for Convective Kinematics) [4]、三阶迎风插值格式CUI (Cubic Upwind Interpolation) [5]。但是根据Godunov [6]定理，高阶线性格式虽然可以很好地逼近光滑解但在间断附近无法避免非物理振荡。接着人们提出了高分辨率格式(High Resolution Scheme) [7] [8]以保证解不会出现非物理振荡。早期高分辨率对流格式的相关研究是由Harten [7]提出的TVD (Total Variational Diminishing Constraint)格式，TVD性质无法区分光滑极值点和间断，从而导致数值解在非单调光滑区域出现精度损失，降为一阶逼近精度。之后Shu对于TVD格式进行了改进，放松了TVD性质的要求，得到一致二阶的TVB格式[9]。这两个格式都能够抑制间断解处和大梯度附近产生的非物理振荡，但TVB格式的不足之处是其中存在与问题相关的一个调节参数。另外一个重要的构造高分辨率格式的方式是结合Leonard提出的NVF (Normalized Variable Formulation)方法[10]与Gaskell和Lau提出的CBC (Convection Boundness Criterion)准则，它使得对流项离散数值格式满足局部有界性，从而保证数值计算的稳定性。CBC格式的特点是在NVF形式下的构造过程简单。与TVD格式类似，典型CBC类数值格式也由于对流有界准则的限制，使得数值解虽然可以避免数值振荡，但是在非单调光滑区域依然出现精度损失，无法达到真正一致二阶精度。

为了解决上述问题，本文结合上述高分辨率格式和QUICK格式来构建一种新的杂交格式。为此引入杂交因子来判断光滑与间断区域。杂交因子决定了杂交格式的性能，它是一个开关函数，其功能是尽可能准确识别解的复杂结构，使其在光滑区域具有高阶精度，从而使得杂交格式在光滑区域呈现线性高阶格式的特点，在解的间断区域又可充分利用高分辨率格式去抑制非物理振荡。在过往的研究中，对于解的局部特征，WENO-Z格式的光滑因子表现出良好的分辨效果。受此启发，本文构造了一个新的杂交因子，并得到相应的杂交格式。数值结果显示杂交格式不仅能克服光滑极值点处的降阶现象，还能降低格式的耗散。

2. 空间离散

一维双曲守恒律方程

$u_t+f{\left(u\right)}_x=0$ (1)

其中 $f\left(u\right)$ 表示流函数。

下面用FVM (Finite Volume Methods)对方程(1)进行空间离散。

下面首先将计算区间 $\left[a,b\right]$ 剖分成N等分

$a=x_{1/2}\le x_{3/2}\le \cdots \le x_{N+1/2}=b$ (2)

这里 $x_{j+1/2}\left(j=0,1,\cdots ,N\right)$ 表示剖分节点， $I_j=\left(x_{j-1/2},x_{j+1/2}\right)$ 表示第j个控制单元，其长度为 $\text{Δ}{x}_j=x_{j+1/2}-x_{j-1/2}$ ，在区间 $I_j$ 上积分可得

${\displaystyle {\int }_{I_j}{u}_t\text{d}x}+{\displaystyle {\int }_{I_j}f{\left(u\right)}_x\text{d}x}=0$ (3)

函数 $u\left(x_j,t\right)$ 在控制单元上 $I_j$ 的单元积分平均值为

$\bar{u}\left(x_j,t\right)=\frac{1}{\text{Δ}x}{\displaystyle {\int }_{I_j}u\left(x_j,t\right)\text{d}x}$ (4)

这里 $x_j=\frac{1}{2}\left(x_{j-1/2}+x_{j+1/2}\right)$ ，因此可以得到

$\frac{\text{d}\bar{u}\left(x_j,t\right)}{\text{d}t}=-\frac{\text{1}}{\text{Δ}x}\left[f\left(u\left(x_{j+1/2},t\right)\right)-f\left(u\left(x_{j-1/2},t\right)\right)\right]$ (5)

用以下格式近似方程(5)

$\frac{\text{d}{\bar{u}}_j\left(t\right)}{\text{d}t}=-\frac{1}{\text{Δ}x}\left({\hat{f}}_{j+1/2}-{\hat{f}}_{j-1/2}\right)$ (6)

这里 ${\bar{u}}_j\left(t\right)$ 是单元积分平均值 $\bar{u}\left(x_j,t\right)$ 的数值近似值， ${\hat{f}}_{j+1/2}$ 是数值流通量。

本文采用Lax-Friedrichs数值流通量[11]，其形式如下

${\hat{f}}_{j+\frac{1}{2}}=\hat{f}\left(u_{j+\frac{1}{2}}^{-},u_{j+\frac{1}{2}}^{+}\right)=\frac{1}{2}\left[f\left(u_{j+\frac{1}{2}}^{+}\right)+f\left(u_{j+\frac{1}{2}}^{-}\right)-\lambda \left(u_{j+\frac{1}{2}}^{+}-u_{j+\frac{1}{2}}^{-}\right)\right]$ (7)

这里 $\lambda ={\max }_{u\in \left(u^{-},u^{+}\right)}\left|f^{\prime }\left(u\right)\right|$ 。

2.1. 高分辨率格式的构造

归一化变量公式(NVF)是描述和分析高分辨率(HR)格式的一个框架，其由Leonard [10]引入。NVF是基于局部归一化变量的面公式，用于构造面f处的变量 ${\varphi }_f$ ，该方法会用到 ${\varphi }_U$ 上游点(upwind)， ${\varphi }_C$ 中点(central)， ${\varphi }_D$ 下游点(downwind)的节点值，如图1所示。

Figure 1. Three neighboring mesh points and the mesh face

图1. 三个相邻节点和单位边界

归一化变量表达式为

$\hat{\varphi }=\frac{\varphi -{\varphi }_U}{{\varphi }_D-{\varphi }_U}$ (8)

根据归一化关系，可得

${\varphi }_f=f\left({\varphi }_U,{\varphi }_C,{\varphi }_D\right)$ (9)

易知 ${\hat{\varphi }}_U=0$ ， ${\hat{\varphi }}_D=1$ ，上述公式可化简为

${\hat{\varphi }}_f=f\left({\hat{\varphi }}_C\right)$ (10)

可知 ${\hat{\varphi }}_f$ 的值只依赖于 ${\varphi }_C$ 的归一化值 ${\hat{\varphi }}_C$ 。

常见的线性高阶格式可以统一写为

${\varphi }_f={\varphi }_C+\left[\frac{1+K}{4}\left({\varphi }_D-{\varphi }_C\right)+\frac{1-K}{4}\left({\varphi }_C-{\varphi }_U\right)\right]$ (11)

对式(11)中的K赋不同的值就可得到不同的高阶格式，其中 $K\in \left[-1,1\right]$ 。

将式(11)正则化得

${\hat{\varphi }}_f=\left(1-\frac{K}{2}\right){\hat{\varphi }}_C+\frac{1}{4}\left(1+K\right)$ (12)

Godunov-Ryabenki [6]定理表明，高阶线性格式不能保持单调性，会产生非物理振荡。因此要想保证解不会出现非物理振荡，需对高阶线性格式施加有界性从而产生高分辨率格式。Lau和Gaskell基于NVF提出了对流有界准则CBC [1]，其后Yu等人[12]和Huo等人[13]进一步提出NVF形式下的BAIR (Boundedness Accuracy and Interpolative Reasonableness)准则。

$\{\begin{array}{l}\frac{3}{2}{\hat{\varphi }}_C\le f\left({\hat{\varphi }}_C\right)\le \frac{1}{2}\left(1+{\hat{\varphi }}_C\right),0<{\hat{\varphi }}_C<\frac{1}{2}\\ \frac{1}{2}\left(1+{\hat{\varphi }}_C\right)\le f\left({\hat{\varphi }}_C\right)\le \frac{3}{2}{\hat{\varphi }}_C\text{and}f\left({\hat{\varphi }}_C\right)\le 1,\frac{1}{2}\le {\hat{\varphi }}_C<1\\ {\hat{\varphi }}_f={\hat{\varphi }}_C,{\hat{\varphi }}_C\le 0\\ {\hat{\varphi }}_f={\hat{\varphi }}_C,{\hat{\varphi }}_C\ge 1\end{array}$ (13)

此外，TVD性质的表达式可以在NV公式中改写为

$\{\begin{array}{l}{\hat{\varphi }}_f\le 1,{\hat{\varphi }}_f\in \left[{\hat{\varphi }}_C,2{\hat{\varphi }}_C\right],0<{\hat{\varphi }}_C<1\\ {\hat{\varphi }}_f={\hat{\varphi }}_C,{\hat{\varphi }}_C\le 0\\ {\hat{\varphi }}_f={\hat{\varphi }}_C,{\hat{\varphi }}_C\ge 1\end{array}$ (14)

图2是CBC-BAIR条件和TVD约束准则。

Figure 2. TVD Constraint Guidelines and CBC-BAIR Conditions

图2. TVD约束准则和CBC-BAIR条件

在NV公式中，HR格式的一般形式可以写为： ${\hat{\varphi }}_f=f\left({\hat{\varphi }}_C\right)$ 。图2中的阴影区域描述了TVD约束和BAIR条件之间的关系。从图2中可以看出，所有已知的线性高阶方案的特征线都超出了NV图中的CBC (BAIR)和TVD区域，要想使用NVD来构造高分辨率有限体积格式，要求其的特征线必须落在NV图中的CBC (BAIR)和TVD区域内，才能满足对流有界性[14]。则有

1) 特征线在NV图中必须过点 $\left(0,0\right),\left(\frac{1}{2},\frac{3}{4}\right),\left(1,1\right)$

2) 在点 $\left(\frac{1}{2},\frac{3}{4}\right)$ 处有二阶局部截断误差

3) 满足 $f^{\prime }\left(\frac{1}{2}\right)=1-\frac{1}{2}K$ ，同时与相应的K格式具有相同的精度

4) $\{\begin{array}{l}{f}^{\prime }\left(0\right)={\theta }_1,\frac{3}{2}\le {\theta }_1\le 2\\ f^{\prime }\left(1\right)={\theta }_2,0\le {\theta }_1\le \frac{1}{2}\end{array}$

上述构造过程使得满足BAIR准则的CBC格式也同时可以具备TVD性质，使得数值计算稳定性得

以加强。本文中令 ${\theta }_1=\frac{17}{10}$ ， ${\theta }_2=\frac{3}{10}$ ， $K=\frac{1}{3}$ ，根据以上的条件，利用Hermite插值设计一个四阶多项式

作为该方案的特征线。

${\hat{\varphi }}_f=\{\begin{array}{l}\frac{1}{10}{\hat{\varphi }}_C\left(12{\hat{\varphi }}_C^3-24{\hat{\varphi }}_C^2+5{\hat{\varphi }}_C+17\right),0<{\hat{\varphi }}_C<1\\ {\hat{\varphi }}_C,\text{ elsewhere}\end{array}$ (15)

图3中红色的曲线为BAIR区域内的高分辨率格式的NV线。通过图3可以看出，构造的高分辨率格式既满足BAIR条件，也满足TVD约束。

Figure 3. NV lines in high resolution format within the BAIR area

图3. BAIR区域内的高分辨率格式的NV线

2.2. 杂交因子

杂交格式构造的关键是设法得到一个有效的杂交因子，它可以有效分辨解的光滑区域与间断区域。在间断附近的计算，杂交格式转换为可以已知数值振荡的高分辨率格式，在光滑解区域，杂交格式保持线性高阶格式的特征。具体而言，本文的高阶线性格式使用三阶QUICK，高分辨率格式使用2.1节构造的高分辨率格式。

${\bar{u}}_j$ 为单元 $I_j$ 上的积分平均值，在模板 $S_0=\left\{I_{j-2},I_{j-1},I_j\right\},S_1=\left\{I_{j-1},I_j,I_{j+1}\right\},S_2=\left\{I_j,I_{j+1},I_{j+2}\right\}$ 上构造二次重构多项式分别记为 $P_k\left(x\right),k=0,1,2$ ，满足

$\frac{1}{\text{Δ}x}{\displaystyle {\int }_{I_{j+k+i-3}}{P}_k\left(x\right)\text{d}x}={\bar{u}}_{j+k+i-3},k=0,1,2,i=0,1,2$ (16)

为了在解的光滑区域获得更高的数值精度，选定大模板 $S_5=\left\{I_{j-2},I_{j-1},I_j,I_{j+1},I_{j+2}\right\}$ ，在5点模板 $S_5$ 上构造四次重构多项式 $P\left(x\right)$ ，且满足

$\frac{1}{\text{Δ}x}{\displaystyle {\int }_{I_{j+i}}P\left(x\right)\text{d}x}={\bar{u}}_{j+i},i=-2,-1,0,1,2,$ (16)

上式在 $I_j$ 边界处满足

$p\left(x_{j+\frac{1}{2}}\right)={\displaystyle \underset{k=0}{\overset{2}{\sum }}{c}_k{P}_k}\left(x_{j+\frac{1}{2}}\right)$ (18)

其中 $c_k$ 为线性权，光滑因子 ${\beta }_k^{}$ [15]取为

${\beta }_k={\displaystyle \underset{l=1}{\overset{2}{\sum }}\text{Δ}{x}^{2l-1}}{\displaystyle {\int }_{I_j}{\left(\frac{{\text{d}}^l}{\text{d}{x}^l}{P}_k\left(x\right)\right)}^2\text{d}x}$ (19)

对于均匀网格，3个光滑因子的计算结果为

$\begin{array}{l}{\beta }_0^{}=\frac{13}{12}{\left({\bar{u}}_{j-2}-2{\bar{u}}_{j-1}+{\bar{u}}_j\right)}^2+\frac{1}{4}{\left({\bar{u}}_{j-2}-4{\bar{u}}_{j-1}+3{\bar{u}}_j\right)}^2\\ {\beta }_1^{}=\frac{13}{12}{\left({\bar{u}}_{j-1}-2{\bar{u}}_j+{\bar{u}}_{j+1}\right)}^2+\frac{1}{4}{\left({\bar{u}}_{j-1}-{\bar{u}}_{j+1}\right)}^2\\ {\beta }_2^{}=\frac{13}{12}{\left({\bar{u}}_j-2{\bar{u}}_{j+1}+{\bar{u}}_{j+2}\right)}^2+\frac{1}{4}{\left(3{\bar{u}}_j-4{\bar{u}}_{j+1}+{\bar{u}}_{j+2}\right)}^2\end{array}$ (20)

下面定义新的计算权重的方法如下

${\alpha }_k=\frac{{\tau }_5}{{\beta }_k+\text{epsilon}},k=\left\{0,1,2\right\}$ (21)

其中 ${\tau }_5$ 由Borges等人[16]提出的一个新的光滑因子，定义为 ${\beta }_0$ 和 ${\beta }_2$ 的绝对差值

${\tau }_5^{}=\left|{\beta }_0^{}-{\beta }_2^{}\right|$ (22)

满足

1) 如果5点模板不包含间断点，则 ${\tau }_5\ll {\beta }_k$

2) 如果解在 $S_k=\left(x_{j-2+k},x_{j-1+k},x_{j+k}\right),k=0,1,2$ 是连续的，但在整个 $S_5$ 中是不连续的，则 ${\beta }_k\ll {\tau }_5$

3) ${\tau }_5>\min \left({\beta }_0,{\beta }_1,{\beta }_2\right)$

因此，杂交因子定义如下

$\theta =\left|\frac{1}{1+{\left({\displaystyle \sum {\alpha }_k}\right)}^2}-1\right|$ (23)

对于包含不连续面的模板，我们可以得出

${\displaystyle \underset{k}{\sum }{\alpha }_k}\gg 1$ (24)

此时

$\theta \approx 1$ (25)

而对于那些光滑模版，光滑度量因子 ${\beta }_k$ 和 ${\tau }_5$ 在 $x_i$ 处进行泰勒展开，可得

$\begin{array}{l}{\beta }_0={u^{\prime }}_j^2\text{Δ}{x}^2+\left(\frac{13}{12}{{u^{\prime }}^{\prime }}_j^2-\frac{2}{3}{u^{\prime }}_j{{{u^{\prime }}^{\prime }}^{\prime }}_j\right)\text{Δ}{x}^4+\left(-\frac{13}{6}{{u^{\prime }}^{\prime }}_j{{{u^{\prime }}^{\prime }}^{\prime }}_j+\frac{\text{1}}{\text{2}}{u^{\prime }}_j{{{{u^{\prime }}^{\prime }}^{\prime }}^{\prime }}_j\right)\text{Δ}{x}^5+{\rm O}\left(\text{Δ}{x}^6\right)\\ {\beta }_1={u^{\prime }}_j^2\text{Δ}{x}^2+\left(\frac{13}{12}{{u^{\prime }}^{\prime }}_j^2-\frac{2}{3}{u^{\prime }}_j{{{u^{\prime }}^{\prime }}^{\prime }}_j\right)\text{Δ}{x}^4+{\rm O}\left(\text{Δ}{x}^6\right)\\ {\beta }_2={u^{\prime }}_j^2\text{Δ}{x}^2+\left(\frac{13}{12}{{u^{\prime }}^{\prime }}_j^2-\frac{2}{3}{u^{\prime }}_j{{{u^{\prime }}^{\prime }}^{\prime }}_j\right)\text{Δ}{x}^4+\left(\frac{13}{6}{{u^{\prime }}^{\prime }}_j{{{u^{\prime }}^{\prime }}^{\prime }}_j-\frac{\text{1}}{\text{2}}{u^{\prime }}_j{{{{u^{\prime }}^{\prime }}^{\prime }}^{\prime }}_j\right)\text{Δ}{x}^5+{\rm O}\left(\text{Δ}{x}^6\right)\end{array}$ (26)

${\tau }_5=\left|\frac{13}{3}{{u^{\prime }}^{\prime }}_j{{{u^{\prime }}^{\prime }}^{\prime }}_j-{u^{\prime }}_j{{{{u^{\prime }}^{\prime }}^{\prime }}^{\prime }}_j\right|\text{Δ}{x}^5+{\rm O}\left(\text{Δ}{x}^6\right)$ (27)

以及

${\alpha }_k=\frac{{\tau }_5}{{\beta }_k+eps}=\{\begin{array}{l}{\rm O}\left(\text{Δ}{x}^3\right),u^{\prime }\ne 0\\ {\rm O}\left(\text{Δ}x\right),u^{\prime }=0\end{array}$ (28)

我们可以得到 $\theta$ 截断误差为

$\theta =\left|\frac{1}{1+{\left({\displaystyle \underset{k}{\sum }{\alpha }_k}\right)}^2}-1\right|=\left(\frac{1}{1+{\rm O}\left(\text{Δ}{x}^p\right)}\right)-1,p=\{\begin{array}{l}\text{6, }{u}^{\prime }\ne 0\\ \text{2, }{u}^{\prime }=0\end{array}$ (29)

上式表明，在光滑模板上， $\theta$ 是0的低阶近似值。很明显

$0<\theta <1$ (30)

事实上，杂交因子可以写成一般形式

$\theta \left(x\right)=\left|\frac{1}{1+x^2}-1\right|$ (31)

其中 $x={\displaystyle \underset{k}{\sum }{\alpha }_k}$ 。

2.3. 杂交格式

基于杂交因子 $\theta$ ，我们采用以下方法来构造保极值的杂交格式来逼近界面值，其格式如下

$F_{Hybrid}=\theta F_{HR}+\left(1-\theta \right)F_{QUICK}$ (32)

其中 $F_{HR}$ 表示2.1节构造的高分辨率格式， $F_{QUICK}$ 表示高阶线性QUICK格式。

显然在光滑的地方的 $\theta$ 值更接近0，这时使用线性高阶格式，杂交格式为三阶精度。而在间断的地方 $\theta$ 的值更接近1，杂交格式使用高分辨率格式。此时 $\theta ={\rm O}\left(\text{Δ}{x}^2\right)$ 同样满足杂交格式的精度是三阶的。

3. 时间离散

空间离散后，得到一个关于时间变量的一阶常微分方程组

$\frac{\text{d}u}{\text{d}t}=L\left(u\right)$ (33)

为了保证杂交格式的三阶精度，本文使用Gottlieb和Shu提出的三阶SSP Runge-Kutta方法[17]对时间进行离散

$\{\begin{array}{l}{\varphi }^{\left(1\right)}={\varphi }^n+L\left({\varphi }^{\left(1\right)}\right)\\ {\varphi }^{\left(2\right)}=\frac{3}{4}{\varphi }^n+\frac{1}{4}\left[{\varphi }^{\left(1\right)}+\text{Δ}tL\left({\varphi }^{\left(1\right)}\right)\right]\\ {\varphi }^n=\frac{1}{3}{\varphi }^n+\frac{2}{3}\left[{\varphi }^{\left(2\right)}+\text{Δ}tL\left({\varphi }^{\left(2\right)}\right)\right]\end{array}$ (34)

4. 数值算例

4.1. 一维线性对流方程

$\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0$ (35)

在 $a=1$ 时，对不同的初始条件进行了数值实验并验证数值格式的精度及准确性。

4.1.1. 情形一

给定光滑初值条件 $u\left(x,0\right)=\sin \left(2\pi x\right)$ ，在区间 $\left[-1,1\right]$ 上具有周期边界条件。对于精度测试， $L_1$ 误差、 $L_2$ 误差、 $L_{\infty }$ 误差的计算时间为 $T=0.1$ ， $CFL=0.1$ 。将高分辨率格式、QUICK格式以及杂交格式的数值解进行比较，计算网格数分别取 $N=\text{10}$ 、20、40、80、160。

$L_1$ 误差、 $L_2$ 误差、 $L_{\infty }$ 误差以及精度阶(见表1)，格式的精度计算阶公式如下

${\Vert E\Vert }_{L_P}={\left(\text{Δ}x{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left|{\bar{u}}_i\left(\text{computed}\right)-{\bar{u}}_i\left(\text{exact}\right)\right|}^p}\right)}^{\frac{1}{P}},p=1,2$ (36)

${\Vert E\Vert }_{L{}_{\infty }}=\underset{1\le i\le N}{\max }\left|{\bar{u}}_i\left(\text{computed}\right)-{\bar{u}}_i\left(\text{exact}\right)\right|$ (37)

$\text{order}=\frac{\ln \frac{E_N}{E_{2N}}}{\ln 2}$ (38)

由表1，对于光滑解问题，三种格式的 $L_1\text{, }{L}_2,L_{\infty }$ 误差随着剖分网格的不断加密在不断的减小而且杂交格式可以达到三阶精度。与高分辨率格式及QUICK格式相比较可知，杂交格式在光滑极值点处有较好的近似效果。

Table 1. $L_1$ 、 $L_2$ and $L_{\infty }$ errors and accuracy rates

表1. $L_1$ 、 $L_2$ 和 $L_{\infty }$ 误差和精度阶

Scheme	Mesh	$L_1$ error	$L_1$ order	$L_2$ error	$L_2$ order	$L_{\infty }$ error	$L_{\infty }$ order
高分辨率	10	3.44E−02	……	4.00E−02	……	8.73E−02	……
	20	1.29E−02	1.42	1.61E−02	1.31	3.40E−02	1.36
	40	3.36E−03	1.82	5.03E−03	1.69	1.40E−02	1.28
	80	8.18E−04	2.16	1.55E−03	1.70	5.57E−03	1.28
	160	2.06E−04	1.99	4.66E−04	1.73	2.18E−03	1.40
QUICK	10	7.70E−03	……	8.68E−03	……	1.19E−02	……
	20	1.02E−03	2.92	1.13E−03	2.94	1.60E−03	2.89
	40	1.29E−04	2.98	1.43E−04	2.98	2.02E−04	2.99
	80	1.61E−05	3.00	1.79E−05	3.00	2.55E−05	2.99
	160	1.96E−06	3.04	2.17E−06	3.04	3.46E−06	2.88
杂交	10	1.93E−02	……	2.61E−02	……	6.88E−02	……
	20	9.98E−04	4.27	1.52E−03	4.10	3.18E−03	4.44
	40	1.25E−04	3.00	1.47E−04	3.37	3.09E−04	3.36
	80	1.59E−05	2.97	1.79E−05	3.03	2.88E−05	3.42
	160	2.00E−06	2.99	2.23E−06	3.00	3.58E−06	3.00

4.1.2. 情形二

给定间断复合初始条件

$u\left(x,0\right)=\{\begin{array}{l}1,1\le x\le 3\\ x-4,4\le x\le 5\\ -x+6,5\le x\le 6\\ \cos \left(0.5\pi \left(x-8\right)\right),7\le x\le 9\\ \exp \left(-4.5{\left(x-11\right)}^2\right),10\le x\le 12\\ 0,\text{ elsewhere}\end{array}$

其中，求解区域是 $\left[0,20\right]$ ， $CFL=0.1$ ， $T=0.5$ ， $N=400$ 。从图4(a)可知，QUICK格式间断解处仍存在一定的数值耗散，高分辨率格式在光滑极值点处存在降阶现象。而和图4(b)中精确解(红色实线)相比，可以看出，杂交格式在不连续点附近没有产生非物理振荡。

(a) (b)

Figure 4. Compare the numerical and accurate solutions in QUICK, high-resolution, and hybridization scheme

图4. 比较QUICK格式、高分辨率格式和杂交格式的数值解与准确解

4.1.3. 情形三

选取如下的W形间断复合初始条件

$u\left(x,0\right)=\{\begin{array}{l}1,0\le x\le 0.2\\ 4x-\frac{3}{5},0.2\le x\le 0.4\\ -4x+\frac{13}{5},0.4\le x\le 0.6\\ 1,0.6\le x\le 0.8\\ 0,\text{ elsewhere}\end{array}$

计算区间 $\left[-0.5,2.5\right]$ ， $N=800$ ， $T=1.0$ ，图5将QUICK、高分辨率格式和杂交格式的数值解与准确解进行了比较。

从图5(a)中蓝色可知，高分辨率格式在间断处的近似效果良好，图5(a)中绿色可知，QUICK格式在大梯度附近产生非物理振荡，而图5(b)中，杂交格式的结果表明，杂交因子可以有效地判断间断使得杂交格式在不连续点附近没有非物理振荡，与精确解的结果逼近效果良好。

(a) (b)

Figure 5. Non-smooth initial value condition (a) HR and QUICK; (b) Hybrid

图5. 非光滑初值条件(a) HR和QUICK；(b) 杂交

4.1.4. 情形四

$u\left(x,0\right)=\{\begin{array}{l}\sin \left(2\pi x\right),x\in \left[0.3,0.8\right]\\ \cos \left(2\pi x\right)-0.5,\text{ elsewhere}\end{array}$

求解区域是 $\left[0,1\right]$ ， $CFL=0.5$ ， $T=0.1$ ，在 $N=80$ 、 $N=160$ 对杂交格式的数值解与准确解进行了比较。图6可以看出杂交格式在不连续点附近不存在非物理振荡，且随着网格的加密和格式精度上升，精确解的吻合度越来越高。

Figure 6. The numerical solution is compared with the exact solution under the mesh number N = 80, 160

图6. 网格数N = 80, 160下数值解与准确解比较

由图6对比可以发现，杂交格式在计算Riemann初值问题时具有良好的效果，杂交格式下的数值解在间断区域没有非物理振荡，而且在光滑极值点处克服了降阶现象。

4.2. 非线性方程

4.2.1. 一维无粘Burgers方程

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial }{\partial x}\left(\frac{u^2}{2}\right)=0\\ u\left(x,0\right)=0.5+\sin \pi x\end{array}$ (39)

周期性边界条件 $0\le u\le 2$ ， $0\le t\le T$ 。

其中时间步长 $\text{Δ}t=CFL\frac{\text{Δ}x}{\underset{i}{\max }\left(\left|u_i\right|\right)}$ 。如图7所示计算域为 $\left[0，2\right]$ ，取 $CFL=0.5$ ， $N=400$ 进行计算。图7(b)中可以看出，杂交格式在不连续点附近没有产生非物理振荡。

(a) (b)

Figure 7. Compare the numerical and exact solutions in hybridization scheme (a) $T=\frac{1}{\pi }$ ; (b) $T=\frac{1.5}{\pi }$

图7. 比较杂交格式的数值解与精确解(a) $T=\frac{1}{\pi }$ ; (b) $T=\frac{1.5}{\pi }$

图7(a)中看出，在时刻 $T=\frac{1}{\pi }$ 时无间断，图7(b)中看出，在时刻 $T=\frac{1.5}{\pi }$ 处产生间断。杂交格式下的数

值解在光滑解处可以保持原有高精度，在间断解处未产生非物理振荡，且具有较好的分辨率。随着网格的细化，数值解近似于精确解。

4.2.2. 一维Buckely-Leverett问题

$\{\begin{array}{l}{u}_t+{\left(\frac{4u^2}{4u^2+{\left(1-u\right)}^2}\right)}_x=0\\ u\left(x,0\right)=\{\begin{array}{l}1,x\in \left[-0.5,0\right]\\ 0,\text{ otherwhere}\end{array}\end{array}$ (40)

时间步长 $\text{Δ}t=CFL\frac{\text{Δ}x}{\underset{i}{\max }\left(\left|u_i\right|\right)}$ ，取 $CFL=0.5$ ，并将计算求解区域 $\left[-1,1\right]$ 等分为200份。图8中的红色

实线为该方程在五阶WENO格式下的数值解。利用杂交格式数值计算该方程，时间 $T=0.4$ ，得到逼近结果如图8的蓝色圆圈。

Figure 8. Comparison of numerical and exact results of the Buckley-Leverett equation

图8. Buckley-Leverett 方程的数值和精确结果的比较

图8将Buckley-Leverett 方程在杂交格式下的数值解和精确解作比较，不难发现，杂交格式下的数值解能够很好的模拟精确解，在间断附近也可以很好地抑制非物理振荡，且杂交因子可以有效的分辨出间断区域与光滑区域。

5. 总结

本文首先构造了一种同时满足BAIR准则和TVD性质的CBC格式，并以此为基础设计一种无振荡保极值杂交格式。借鉴WENO-Z方法的光滑因子设计了一种开关函数，实现在间断区域使用高分辨率格式捕捉间断而在光滑区域使用高阶的迎风线性重构，从而达到在有效地消除非物理振荡及光滑极值点降阶现象的同时保持良好的数值逼近效果。

致 谢

作者感谢审稿专家给予本文宝贵的修改意见和建议以及期刊编辑对于本文的付出。

基金项目

内蒙古大学科研发展基金(21100-5187133)和内蒙古自治区人才开发基金项目(12000-1300020240)。

NOTES

^*通讯作者。

参考文献

[1]	Spalding, D.B. (1972) A Novel Finite Difference Formulation for Differential Expressions Involving Both First and Second Derivatives. International Journal for Numerical Methods in Engineering, 4, 551-559. https://doi.org/10.1002/nme.1620040409
[2]	Shyy, W. (1985) A Study of Finite Difference Approximations to Steady-State, Convection-Dominated Flow Problems. Journal of Computational Physics, 57, 415-438. https://doi.org/10.1016/0021-9991(85)90188-3
[3]	Lax, P. and Wendroff, B. (1960) Systems of Conservation Laws. Communications on Pure and Applied Mathematics, 13, 217-237. https://doi.org/10.1002/cpa.3160130205
[4]	Leonard, B.P. (1979) A Stable and Accurate Convective Modelling Procedure Based on Quadratic Upstream Interpolation. Computer Methods in Applied Mechanics and Engineering, 19, 59-98. https://doi.org/10.1016/0045-7825(79)90034-3
[5]	Agarwal, R. (1981) A Third-Order-Accurate Upwind Scheme for Navier-Stokes Solutions at High Reynolds Numbers. 19th Aerospace Sciences Meeting, St. Louis, 12-15 January 1981, 12-15. https://doi.org/10.2514/6.1981-112
[6]	Roe, P.L. (1986) Characteristic-Based Schemes for the Euler Equations. Annual Review of Fluid Mechanics, 18, 337-365. https://doi.org/10.1146/annurev.fl.18.010186.002005
[7]	Harten, A. (1983) High Resolution Schemes for Hyperbolic Conservation Laws. Journal of Computational Physics, 49, 357-393. https://doi.org/10.1016/0021-9991(83)90136-5
[8]	Sweby, P.K. (1984) High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws. SIAM Journal on Numerical Analysis, 21, 995-1011. https://doi.org/10.1137/0721062
[9]	Shu, C. (1987) TVB Uniformly High-Order Schemes for Conservation Laws. Mathematics of Computation, 49, 105-121. https://doi.org/10.1090/s0025-5718-1987-0890256-5
[10]	Leonard, B.P. (1988) Simple High‐accuracy Resolution Program for Convective Modelling of Discontinuities. International Journal for Numerical Methods in Fluids, 8, 1291-1318. https://doi.org/10.1002/fld.1650081013
[11]	Qiu, J., Khoo, B.C. and Shu, C. (2006) A Numerical Study for the Performance of the Runge-Kutta Discontinuous Galerkin Method Based on Different Numerical Fluxes. Journal of Computational Physics, 212, 540-565. https://doi.org/10.1016/j.jcp.2005.07.011
[12]	Yu, W. Q. Tao, D. S. Zhang, Q. W. W, B. (2001) Discussion on Numerical Stability and Boundedness of Convective Discretized Scheme. Numerical Heat Transfer, Part B: Fundamentals, 40, 343-365. https://doi.org/10.1080/104077901317091721
[13]	Hou, P.L., Tao, W.Q. and Yu, M.Z. (2003) Refinement of the Convective Boundedness Criterion of Gaskell and Lau. Engineering Computations, 20, 1023-1043. https://doi.org/10.1108/02644400310503008
[14]	Peng, J. and Shen, Y.Q. (2017) A Novel Weighting Switch Function for Uniformly High-Order Hybrid Shock-Capturing Schemes. International Journal for Numerical Methods in Fluids, 9, 681-703.
[15]	Blossey, P.N. and Durran, D.R. (2008) Selective Monotonicity Preservation in Scalar Advection. Journal of Computational Physics, 227, 5160-5183. https://doi.org/10.1016/j.jcp.2008.01.043
[16]	Borges, R., Carmona, M., Costa, B. and Don, W.S. (2008) An Improved Weighted Essentially Non-Oscillatory Scheme for Hyperbolic Conservation Laws. Journal of Computational Physics, 227, 3191-3211. https://doi.org/10.1016/j.jcp.2007.11.038
[17]	Gottlieb, S. and Shu, C. (1998) Total Variation Diminishing Runge-Kutta Schemes. Mathematics of Computation, 67, 73-85. https://doi.org/10.1090/s0025-5718-98-00913-2