对流占优问题的非线性QUICK格式

A Nonlinear QUICK Scheme for Convection-Dominated Equations

摘要: 本文构造了一种新的非线性加权格式，利用这种加权格式对对流项进行离散，在加权系数中构造了新的间断阈值来判断光滑与间断，使得此格式在间断的时候利用一阶迎风格式，这样可以避免振荡，而在光滑解部分用三阶QUICK格式保证了三阶精度，然后利用三阶Runge-Kutta方法对时间进行离散，进而保证整体精度，使得数值解达到比较好的逼近效果。

Abstract: This paper constructed a nonlinear weighted scheme. The convective term is discretized by this new nonlinear weighted scheme. A new discontinuous threshold in the weighted coefficient is constructed to judge the smoothness and discontinuity, so that the scheme uses the first-order upwind scheme near the discontinuity which the oscillation can be avoided, and the third-order QUICK scheme is used to ensure the third-order accuracy in the smooth region. Time discretization is fulfilled by using the third order Runge-Kutta scheme. The new scheme achieves optimal order accuracy. It makes the numerical solution achieve a better approximation effect.

1. 引言

对流占优方程的间断初值问题一直以来都是数值计算的重要内容。线性格式有中心差分格式CD (Central Difference)，二阶迎风格式SOU (Second-order Upwind)，QUICK [1] (Quadratic Upstream Interpolation for Convective Kinematics)对于光滑初值都能达到二阶及以上的精度。对于间断初值问题，线性格式由于不满足对流有界性，容易引起非物理振荡。为此Gaskell和Lau提出了CBC (Convection Boundedness Criterion) [2] ，可以保证数值格式的有界性，但是未能保证数值格式的精度。为了克服这个缺点，Hou提出了BAIR [3] 准则。1983年Harten [4] 提出了TVD (Total Variation Diminishing)格式，满足了有界性的要求且具有高分辨率的优点，在间断处和极值点处却达不到三阶及以上的精度，尤其在极值点附近由于差分格式的解只具有一阶精度。考虑到前面几个格式的不足，本文提出了格式一种非线性QUICK格式，使得在间断处利用一阶迎风格式，在连续处利用三阶QUICK格式。于是如何判断间断与光滑成了关键问题。本文分析了 [5] 中的光滑因子的特性，考虑到 [6] 中关于平滑度度量的不足，即对于间断和光滑极值不能够明确的界定，于是本文设计了一种能够有效捕捉间断位置的方法，使得对间断和光滑极值的判断清晰而准确。通过算例分析，对于光滑解，达到一致高精度。对于间断解，可由有效抑制数值振荡，逼近效果良好。

本文整体结构安排如下：第一节本文用有限体积格式对守恒律方程进行空间离散，用Lax-Friedrichs [7] 型数值流通量对通量进行数值逼近。然后通过间断阈值对解的间断与光滑处进行判定，构造了非线性加权新格式来离散对流项。第二节对时间进行离散。第三节为一些数值算例及其结果的分析与讨论。第四节得出结论。

2. 空间离散格式

2.1. 有限体积离散

一维守恒律方程为

其中 $f\left(u\right)$ 为流通量。用有限体积法对(1)进行离散，首先将计算区间剖分为N等分， $a=x_{\frac{1}{2}}\le x_{\frac{3}{2}}\le \cdots \le x_{N+\frac{1}{2}}=b$ 且每个控制单元为 $I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right]$ ，其长度为 $\Delta x=\frac{b-a}{N}$ ，对(1)在区间 $I_j$ 上进行积分，则有：

${\displaystyle {\int }_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}{u}_t\text{d}x}+{\displaystyle {\int }_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}f{\left(u\right)}_x\text{d}x}=0$

令 $\frac{{\displaystyle \int u\left(x_j,t\right)\text{d}x}}{\Delta x}=\bar{u}\left(x_j,t\right)$ ，称其为积分平均值，于是可以得到

$\frac{\text{d}\bar{u}\left(x_j,t\right)}{\Delta x}=-\frac{f\left(u_{j+\frac{1}{2}},t\right)-f\left(u_{j-\frac{1}{2}},t\right)}{\Delta x}$

利用数值通量逼近上式如下：

$\frac{\text{d}\bar{u}\left(x_j,t\right)}{\Delta x}=-\frac{{\hat{f}}_{j+\frac{1}{2}}-{\hat{f}}_{j-\frac{1}{2}}}{\Delta x}$

其中 ${\bar{u}}_j\left(t\right)$ 是 $$ 的数值近似，数值流通量 ${\hat{f}}_{j+\frac{1}{2}}$ 来近似逼近 $f\left(u_{j+\frac{1}{2}},t\right)$ ，其中

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

Lax-Friedrichs型的数值流通量为：

$\hat{f}\left(u^{+},u^{-}\right)=\frac{1}{2}\left[f\left(u^{+}\right)+f\left(u^{-}\right)-\alpha \left(u^{+}-u^{-}\right)\right]$ (1)

其中 $\alpha =\underset{u\in \left(u^{+},u^{-}\right)}{\max }\left|f^{\prime }\left(u\right)\right|$ 。

2.2. 对流项离散

本文构造了的非线性权的QUICK格式来逼近对流项 $$ ，有如下格式

$u_{j+\frac{1}{2}}^{-}=\left(-\frac{1}{6}+\frac{1}{6}\omega \right){\bar{u}}_{j-1}+\left(\frac{5}{6}+\frac{1}{6}\omega \right){\bar{u}}_j+\left(\frac{2}{6}-\frac{2}{6}\omega \right){\bar{u}}_{j+1}$ (2)

$$ (3)

其中 $\omega$ 为权，当解在光滑区域时，此时，数值格式变为三阶QUICK格式；当解在间断区域时，此时 $\omega =1$ ，数值格式变为一阶迎风格式，一阶格式可以抑制振荡，于是如何判断光滑与间断成为了关键问题。

为了有效的判断光滑与间断，使得格式在光滑处利用三阶QUICK格式，在间断处利用一阶迎风格式。根据对文献 [8] [9] [10] 里k个模板下的光滑因子 ${\beta }_r$ 进行分析，

${\beta }_r=\underset{l=1}{\overset{k-1}{{\displaystyle \sum }}}\Delta x^{2l-1}{\displaystyle {\int }_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}{\left(\frac{{\text{d}}^l}{\text{d}{x}^l}{p}_r\left(x\right)\right)}^2\text{d}x}$ (4)

其中 $$ 是在k个模板 $s_r\left(i\right)=\left\{x_{i-r},\cdots ,x_{i-r+k-1}\right\},r=0,1,2\cdots ,k-1$ 下构建的多项式。(4)右侧的积分是在 $I_j$ 上插值多项式的所有导数的在 $L_2$ 范数下的平方和。可知当 $k=3$ 时，根据文献 [11] [12] [13] 中对光滑度的测量 ${\beta }_r$ ，我们引入如下形式的光滑因子，

$r\left(j\right)={\left({\bar{u}}_{j+1}-2{\bar{u}}_j+{\bar{u}}_{j-1}\right)}^2+{\left({\bar{u}}_{j+1}-{\bar{u}}_{j-1}\right)}^2$

通过分析计算可知，在连续时 $r\left(j\right)$ 之间几乎相同，而在间断处 $r\left(j\right)$ 之间存在较大的变化， $r\left(j\right)$ 对数值解的间断部分比较敏感，通过分析文献 [6] ，我们给出了的如下光滑参数：

可以很好的衡量数值解的陡度和光滑程度，在光滑处通过泰勒展开可知 ${\lim }_{\Delta x\to 0}\varnothing =1$ ，而在间断处 $\varnothing$ 的值变化比较大，由于 $\varnothing$ 值在间断处的变化不能够明确的界定间断与光滑，为了更加明确的判断，分析了以上的光滑因子与光滑参数的特性，我们给出 $\varnothing$ 的界定值，构造出了如下间断阈值：

$ratio=\frac{{\left(\max \left\{{\bar{u}}_j\right\}-\min \left\{{\bar{u}}_j\right\}\right)}^2}{\max \left\{r\left(j-1\right),r\left(j\right),r\left(j+1\right)\right\}+dd{r}_{\max }},j=1,2,\cdots ,N$

其中 $dd{r}_{\max }=\max \left\{{\left({\bar{u}}_{j+1}-2{\bar{u}}_j+{\bar{u}}_{j-1}\right)}^2\right\}$ ，通过大量实验发现在间断处有 $\varnothing >ratio$ ，而在光滑处有 $$ ，根据算例发现通过 与 的比较确实能够有效的判断间断与光滑区域。

3. 时间离散

对空间进行离散后得到关于时间的常微分方程组，为了保持整体格式的三阶精度，对时间利用三阶Runge-Kutta方法来进行离散，并且可以有效避免振荡，三阶Runge-Kutta格式如下：

$\begin{array}{l}{u}^{\left(1\right)}=u^n+\Delta tL\left(u^n\right)\\ u^{\left(2\right)}=u^n+\frac{1}{4}\left[u^{\left(1\right)}+\Delta tL\left(u^{\left(1\right)}\right)\right]\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}\left[u^{\left(2\right)}+\Delta tL\left(u^{\left(2\right)}\right)\right]\end{array}$

4. 数值算例

4.1. 线性对流方程

$\{\begin{array}{l}\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0\\ u\left(x,0\right)=u_0\left(x\right)\end{array},x\in \left[a,b\right],t>0$ (4)

4.1.1. 情形1

在给定初始条件 $u\left(x,0\right)=\sin \left(2\text{π}x\right),x\in \left[0,1\right]$ ， $$ ，利用新格式计算当 $\Delta x=20,40,80,160$ 时，时间t = 0.1时的数值解，并求出格式的 $L_1,L_{\infty }$ 误差及其精度阶(如表1)，计算公式如下：

$E_{L_1}={\Vert u^{num}-u^{exact}\Vert }_{L_i}=\Delta x{\displaystyle {\sum }_{j=1}^n\left|u_j^{num}-u_j^{exact}\right|}$

$E_{L_{\infty }}={\Vert u^{num}-u^{exact}\Vert }_{L_{\infty }}=\underset{1\le j\le N}{\max }\left|u_j^{num}-u_j^{exact}\right|$

精度阶的计算公式： $order=\frac{\ln \left(E_N/E_{2N}\right)}{\ln 2}$

Table 1. Errors and orders for schemes

表1. 格式误差与数值精度阶对比

由表可知此格式对于光滑初值条件下，新格式的数值解能够达到三阶精度，具有较高的分辨率。

4.1.2. 情形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\text{π}\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}$

利用新格式计算当网格数 $N=200$ 时， $\text{CFL}=0.\text{5}$ ， $x\in \left[0,20\right]$ ，时间 $t=0.5$ 时的数值解并与精确解在图1中进行比较，精确解与QUICK格式下的数值解在图2中进行比较。

Figure 1. Numerical solution of the above two schemes and the exact solution with initial condition of composite discontinuity

图1. 间断复合初始条件下的非线性QUICK格式和QUICK格式的数值解与精确解

由图1可知，QUICK格式与非线性QUICK格式相比，QUICK格式在间断处要产生的振荡，而新格式在间断处容易避免振荡，在光滑处的逼近效果也比较好。

4.1.3. 情形3

在给定间断初始条件

$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<x\le 0.6\\ 1,0.6<x\le 0.8\\ 0,\text{elsewhere}\end{array}$ , $x\in \left[-1,2\right],\text{CFL}=0.1,T=0.2$

当网格数 $$ ，新格式的数值解及QUICK格式与精确解图像如图2。

Figure 2. Numerical solution of the above two schemes and the exact solution with W-shaped initial condition

图2. 间断初始条件下的非线性QUICK格式和QUICK格式的数值解与精确解

由图2可知QUICK的数值解在间断处会产生振荡，而新格式的数值解与精确解相比较，在间断处没有振荡，相比较新格式达到了比较好的逼近效果。

4.1.4. 情形4

对于有间断及相邻近处有极值的初始条件而言，对于格式的要求比较高，对于前面提到的格式，会在间断处产生振荡，不能够很好的区分光滑极值和间断，容易混淆这两种情况，从而在极值处的逼近效果不是很理想，下面我们给出在新格式下的数值解。

在给定有间断和极值的初始条件如下

$u\left(x,0\right)=\{\begin{array}{l}\sqrt{1-{\left(\frac{10x}{3}\right)}^2},x\in \left[-0.3,0.3\right]\\ 0,其他\end{array},x\in \left[-1,1\right],\text{CFL}=0.\text{5}$

利用新格式计算当 $N=600$ 时，时间 $$ 时的数值解和QUICK格式的数值解与精确解的比较如图3进行比较。

由图3可知，对于有光滑极值和极值邻近处有间断的初值问题，由数值解图像可知新格式无论是在间断处还是在极值点附近其逼近效果都非常好，而QUICK格式在间断处却会产生振荡，可知非线性QUICK格式比QUICK格式的逼近效果要好。

4.2. 无粘Burgers方程

无粘burgers方程解会随着时间的推移产生间断，给数值求解带来了很大困难，本文利用新格式求解此方程如下，

Figure 3. Numerical solution of the above two schemes and the exact solution with complicated structures

图3. 复杂结构初始条件下的非线性QUICK格式和QUICK格式的数值解与精确解

$\frac{\partial u}{\partial t}+\frac{\partial \left(\frac{u^2}{2}\right)}{\partial x}=0,x\in \left[0,2\right],t=\frac{1.5}{\text{π}}$

其初始条件为 $u\left(x,0\right)=0.5+\text{sin}\left(2\text{π}x\right)$ ，新格式的数值解图像及QUICK格式的数值解的图像如图4：

Figure 4. Numerical solution of the above two schemes and the exact solution of Burgers equation

图4. Burgers方程的非线性QUICK格式和QUICK格式的数值解与精确解

由新格式与QUICK格式的数值解与精确解的图像比较可发现，在随着时间的推移，在逼近解的间断处时，QUICK格式下的数值解会有明显的振荡，而在非线性QUICK格式下的数值解会避免振荡，有很好的逼近效果。

5. 结论

本文基于有限体积下对双曲守恒方程进行数值离散，利用三阶QUICK格式与一阶迎风格式(FOU)加权后得到新的非线性QUICK格式来逼近对流项，本文通过引入间断阈值使得对流项的离散格式在间断处利用一阶迎风格式并且在光滑处使用三阶QUICK格式，通过算例分析，对于光滑初值能够达到三阶以上精度，对于有间断和极值的初值问题也可以避免振荡，使得数值解达到了比较好的逼近效果。

基金项目

感谢内蒙古自治区研究生科研创新项目(1402020201-46)及内蒙古自然科学基金项目(2015MS0101)和内蒙古自治区人才开发基金项目(12000-1300020240)的支持。