FVM in CFD 学习笔记_第11章_对流项离散

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 11 Discretization of the Convection Term

在前面的第8-9章，咱们详细讲解了在正交、非正交、结构、非结构网格上的通有稳态扩散方程的离散方法。本章咱们讲讲CFD的控制方程中另一个非常重要的项，对流项的离散方法。最初，跟扩散项中所采用的离散方式一样，对流项也是对物理量采用对称和线性分布（廓线profile）假设来离散的，然而，这种分布廓线有很大缺陷，促使人们提出了使用迎风廓线来修正其缺陷。尽管迎风廓线可以得到物理上说得通的结果，然而其被表明是高度diffusive（扩散？）的，导致结果只有一阶精度。为了提高精度，提出了偏迎风的高阶廓线。离散误差倒是降低了，但是高阶廓线却引出了另一种形式的误差，dispersion误差（色散？），关于这种误差的处理方法将放到下一章讲解。还有一点需要说明，我们假设驱动扩散项的流场，是事先已知的。那么流场的计算将放到后续章节里详解（可能你已经猜到了，因为离散方程的系数矩阵本身就是与流场量相关的，哪怕是定常流场，也需要迭代计算，而不可压缩中还有压力速度的耦合问题，还得构造解耦算法来折腾呢）。

1 引言

尽管对流项看起来很简单，然而其离散方法却非常麻烦，以至于人们研究了30来年。这些工作给影响其离散的问题照亮了道路，也催生了大量的离散格式。鉴于文献的数量是如此庞大，咱们不得不用两个章节来分析它们。本章将引入基本概念和高阶（High Order，缩写为HO）偏迎风格式。下一章会讲到对流项的界限（bounding）问题，其有助于发展高阶精度的无振荡对流格式，即，所谓的高分辨率（High Resolution（HR））格式。

为清晰起见，新概念的引入将首先在1维网格上进行，然后再扩展到多维非正交网格上。本章将用1维对流-扩散问题的离散作为开端，通过稳定性判据，表明中心差分算法的缺陷。随后，表明迎风格式可通过该稳定性测试。接着解释了低阶格式的数值diffusion（扩散）误差和高阶格式的数值dispersion（色散）误差。而对高分辨率格式的处理则放到下一章进行讲解。

2 稳态一维对流和扩散

为便于说明问题，给出一个非常简单的一维定常对流-扩散方程，其控制方程为
$$\frac{d(\rho u\varphi )}{dx}-\frac{d}{dx}\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\right)=0$$
万分幸运的是，这个问题存在着解析解，所以它通常用作基准算例来跟各种各样的数值格式做对照。

2.1 解析解

如果截面是不变的，那么这个稳态1维问题的连续方程是
$$\frac{d(\rho u)}{dx}=0$$
这表明$\rho u$是定值，那么对流扩散方程变成了
$$\rho u\frac{d\varphi }{dx}-\frac{d}{dx}\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\right)=0$$
对$x$积分，得
$$\rho u\varphi -{\Gamma }^{\varphi }\frac{d\varphi }{dx}=c_1$$
$c_1$为积分常数，由边界条件来定。上面这个方程进一步写成
$$\frac{d\varphi }{dx}=\frac{\rho u}{{\Gamma }^{\varphi }}\varphi -\frac{c_1}{{\Gamma }^{\varphi }}$$
引入变量替换，定义$\Phi$
$$\Phi =\frac{\rho u}{{\Gamma }^{\varphi }}\varphi -\frac{c_1}{{\Gamma }^{\varphi }}$$
则方程变为
$$\frac{d\Phi }{dx}=\frac{\rho u}{{\Gamma }^{\varphi }}\Phi$$
即
$$\frac{d\Phi }{\Phi }=\frac{\rho u}{{\Gamma }^{\varphi }}dx$$
两侧同时积分，得
$$Ln(\Phi )=\frac{\rho u}{{\Gamma }^{\varphi }}x+c_3$$
其中$c_3$为积分常数，进一步写为
$$\Phi =c_2{e}^{\frac{\rho u}{{\Gamma }^{\varphi }}x}$$
其中$c_2=e^{c_3}$为积分常数，再把$\Phi$替换回最初的变量$\varphi$，得到解析解
$$\varphi =\frac{c_2{\Gamma }^{\varphi }{e}^{\frac{\rho u}{{\Gamma }^{\varphi }}x}+c_1}{\rho u}$$
一维问题的网格剖分如下图所示

以上图的$W$和$E$两节点之间的部分为例，边界条件为
$$\{\begin{array}{llc}\varphi ={\varphi }_W& at& x=x_W\\ \varphi ={\varphi }_E& at& x=x_E\end{array}$$
可根据该边界条件求得积分常数$c_1$和$c_2$，从而获得解析解为
$$\frac{\varphi -{\varphi }_W}{{\varphi }_E-{\varphi }_W}=\frac{e^{P{e}_L\frac{x-x_W}{L}}-1}{e^{P{e}_L}-1}$$
其中$P{e}_L$为Peclet数（基于长度$L$），代表了变量$\varphi$的对流输运速率与扩散输运速率的比值，即
$$P{e}_L=\frac{\rho uL}{{\Gamma }^{\varphi }}L=x_E-x_W$$
不同的$P{e}_L$下所得到的不同结果如下图所示，可见，$P{e}_L=0$时对应的是纯扩散问题，$\varphi$从$W$到$E$的变化是一条直线，而随着$P{e}_L$的增大，对流输运速率所占比重越来越大，$\varphi$从$W$到$E$的变化也逐渐变为了更加陡峭的廓线形式。

2.2 数值解

在1维单元上对于1维对流扩散方程$\frac{d(\rho u\varphi )}{dx}-\frac{d}{dx}\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\right)=0$进行数值积分，得
$${\int }_{V_C}[\nabla \cdot (\rho \mathbf{v}\varphi )-\nabla \cdot ({\Gamma }^{\varphi }\nabla \varphi )]dV=0$$
其中$\mathbf{v}=u\mathbf{i}$为速度矢量，守恒方程可以写成对流通量和扩散通量的形式，即
$${\int }_{V_C}\nabla \cdot ({\mathbf{J}}^{\varphi ,C}+{\mathbf{J}}^{\varphi ,D})dV=0$$
其中对流通量${\mathbf{J}}^{\varphi ,C}=\rho \mathbf{v}\varphi$，扩散通量${\mathbf{J}}^{\varphi ,D}=-{\Gamma }^{\varphi }\nabla \varphi$。接下来使用散度定理，把体积分转化成面积分，即
$${\int }_{V_C}\nabla \cdot ({\mathbf{J}}^{\varphi ,C}+{\mathbf{J}}^{\varphi ,D})dV={\int }_{\partial V_C}({\mathbf{J}}^{\varphi ,C}+{\mathbf{J}}^{\varphi ,D})\cdot d\mathbf{S}={\int }_{\partial V_C}\left(\rho u\varphi \mathbf{i}-{\Gamma }^{\varphi }\frac{d\varphi }{dx}\mathbf{i}\right)\cdot d\mathbf{S}=0$$
将面积分用流过单元面的通量加和形式表示，即
$$\sum _{f\sim nb(C)}{\left(\rho u\varphi \mathbf{i}-{\Gamma }^{\varphi }\frac{d\varphi }{dx}\mathbf{i}\right)}_f\cdot {\mathbf{S}}_f=0$$
注意，如果面积矢量${\mathbf{S}}_f$是沿着反方向的（$-x$方向），那么它的符号是负的，对于恒定截面的情况，上式变为
$$\left[(\rho u\Delta y\varphi {)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e\right]-\left[(\rho u\Delta y\varphi {)}_w-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]=0$$
单元面上的速度$u$是已知的（$u_e,u_w$已知），那些个梯度项（$(d\varphi /dx{)}_e,(d\varphi /dx{)}_w$）用之前讲过的离散方法可以得到。现在的问题是如果继续离散的话，面上的值${\varphi }_e$和${\varphi }_w$需要用邻近的节点值来得到，确定这些面上值的方法就是文献中所谓的“对流格式”。

2.3 初步推导：中心差分格式（The Center Difference （CD） Scheme）

乍一看，好像没啥难的啊，就跟在扩散项离散中采用的线性插值廓线一样，来做这个对流离散不就好了么，即，在给定面比如右侧面$e$上的$\varphi$值的计算，可以假设变量的空间分布为
$$\varphi (x)=k_0+k_1(1-x_C)$$
其中的$k_0$和$k_1$为系数，使用$\varphi ={\varphi }_E,x=x_E$和$\varphi ={\varphi }_C,x=x_C$，可得
$${\varphi }_e={\varphi }_C+\frac{{\varphi }_E-{\varphi }_C}{x_E-x_C}(x_e-c_C)$$
这便是中心差分格式，其可通过Taylor展开并略去2阶以上的项来推出，即其有2阶精度。

在均匀网格上，该离散型式为
$${\varphi }_e=\frac{{\varphi }_C+{\varphi }_E}{2}$$
这样，对于上一小节最后得到的离散格式中，第1个中括号项$\left[(\rho u\Delta y\varphi {)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e\right]$变为
$$\begin{array}{rl}\left[(\rho u\Delta y\varphi {)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e\right]& =(\rho u\Delta y{)}_e\frac{{\varphi }_C+{\varphi }_E}{2}-{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)\\ & =Flux{C}_e{\varphi }_C+Flux{F}_e{\varphi }_E+Flux{V}_e\end{array}$$
其中
$$\begin{array}{rl}& Flux{C}_e={\Gamma }_e^{\varphi }\frac{\Delta y_e}{\delta x_e}+\frac{(\rho u\Delta y{)}_e}{2}\\ & Flux{F}_e=-{\Gamma }_e^{\varphi }\frac{\Delta y_e}{\delta x_e}+\frac{(\rho u\Delta y{)}_e}{2}\\ & Flux{V}_e=0\end{array}$$
同样，第2个中括号项$-\left[(\rho u\Delta y\varphi {)}_w-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]$变为
$$\begin{array}{rl}-\left[(\rho u\Delta y\varphi {)}_w-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]& =-\left[(\rho u\Delta y{)}_w\frac{{\varphi }_C+{\varphi }_W}{2}-{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]\\ & =Flux{C}_w{\varphi }_C+Flux{F}_w{\varphi }_W+Flux{V}_w\end{array}$$
其中
$$\begin{array}{rl}& Flux{C}_w={\Gamma }_w^{\varphi }\frac{\Delta y_w}{\delta x_w}-\frac{(\rho u\Delta y{)}_w}{2}\\ & Flux{F}_w=-{\Gamma }_w^{\varphi }\frac{\Delta y_w}{\delta x_w}-\frac{(\rho u\Delta y{)}_w}{2}\\ & Flux{V}_w=0\end{array}$$
从而推得对流扩散方程的离散形式为
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W=0$$
其中
$$\begin{array}{rl}& a_E=Flux{F}_e=-{\Gamma }_e^{\varphi }\frac{\Delta y_e}{\delta x_e}+\frac{(\rho u\Delta y{)}_e}{2}\\ & a_W=Flux{F}_w=-{\Gamma }_w^{\varphi }\frac{\Delta y_w}{\delta x_w}-\frac{(\rho u\Delta y{)}_w}{2}\\ & a_C=Flux{C}_e+Flux{C}_w={\Gamma }_e^{\varphi }\frac{\Delta y_e}{\delta x_e}+\frac{(\rho u\Delta y{)}_e}{2}+{\Gamma }_w^{\varphi }\frac{\Delta y_w}{\delta x_w}-\frac{(\rho u\Delta y{)}_w}{2}\end{array}$$
由于问题是一维的，有$\Delta y_e=\Delta y_w$，故可把它们都设成1。此外，从连续方程可得到$(\rho u\Delta y{)}_e-(\rho u\Delta y{)}_w=0$，且$u$为常数。假设扩散系数为定值，则离散方程系数变为
$$\begin{array}{rl}& a_E=-\frac{{\Gamma }^{\varphi }}{x_E-x_C}+\frac{(\rho u{)}_e}{2}\\ & a_W=-\frac{{\Gamma }^{\varphi }}{x_C-x_W}-\frac{(\rho u{)}_w}{2}\\ & a_C=-(a_E+a_W)\end{array}$$
回代到离散格式$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W=0$，可得
$$\frac{{\varphi }_C-{\varphi }_W}{{\varphi }_E-{\varphi }_W}=\frac{a_E}{a_E+a_W}$$
如果是均匀网格，可以写成$P{e}_L=\frac{\rho uL}{{\Gamma }^{\varphi }},L=x_E-x_W$的形式，单元$C$的数值解为（注意，书中给出的公式RHS为$\frac{1}{2}(1-\frac{P{e}_L}{2})$实际应为$\frac{1}{2}(1-\frac{P{e}_L}{4})$）
$$\frac{{\varphi }_C-{\varphi }_W}{{\varphi }_E-{\varphi }_W}=\frac{1}{2}\left(1-\frac{P{e}_L}{4}\right)$$
单元$C$的精确解，由上一小节的解析解给出
$$\frac{{\varphi }_C-{\varphi }_W}{{\varphi }_E-{\varphi }_W}=\frac{e^{P{e}_L/2}-1}{e^{P{e}_L}-1}$$
让$P{e}_L$从-10变化到10，数值解和解析解给出的结果展示在下图中，当$P{e}_L$值（绝对值）较小时，数值解和解析解非常接近，当$P{e}_L$的值（绝对值）逐渐增大并跨过阈值后，中心差分格式给出的数值解就和解析解产生了很大偏差，其变得不受约束（无界）呈现出不符合物理意义的特性。即，解析解对于正的和负的$P{e}_L$分别逐渐趋近于0和1，中心差分格式给出的数值解则随着$P{e}_L$从$-\infty$增大到$+\infty$而从$+\infty$线性减小到$-\infty$。这表明前面离散过程中采用的某些假设是不切实际的或者说是不符合物理意义的，到底是什么原因造成的呢？

如下图所示，在点$C$处的扩散所受到$C$点上游和下游条件的影响是等效的（下图a），而对流过程则是与方向高度相关的过程，其仅在流动方向上有输运特性（下图b）。因此，给上下游节点赋予了同样的权重的线性廓线假设，对于扩散项的近似是没问题的（下图a），却没办法描述对流项的方向特性，此时陡峭的廓线更为合适（下图b），这便是造成前述中心差分格式数值解和解析解之间偏差的原因。

下图c给出了对流-扩散联合效应下的影响区域以及更加实际的廓线，该影响区域在低Peclet数下趋近于扩散区域（下图a），在高Peclet数下则趋近于对流区域（下图b）。因此，只要扩散在传输特性中占优，那么使用线性廓线能产生符合物理实际的解。然而，一旦对流压倒了扩散（问题变得对流占优了），再用线性廓线的假设就会产生不符合物理意义的解。可以容易地估算出Peclet数在何值下会产生该现象，假设流动是沿着$x$正方向的，注意到$a_E$系数有可能会变成正值，这样便会导致非物理解（如果流动是沿着$x$负方向的话，那么$a_W$可能会变成正值），即

$$a_E=-{\Gamma }_e^{\varphi }\frac{\Delta y_e}{\delta x_e}+\frac{(\rho u\Delta y{)}_e}{2}>0\Rightarrow \frac{(\rho u{)}_e\delta x_e}{{\Gamma }_e^{\varphi }}>2$$
定义Peclet数为
$$Pe=\frac{\rho u\delta x}{{\Gamma }^{\varphi }}$$
这样，对于均匀网格来说，$Pe=P{e}_L/2$，那么$a_E>0$的条件变为
$$Pe>2$$
因此，单元的Peclet数（Pe）大于2的话，离散过程变得不相符，因为此时邻居值的增加会导致$C$处值的降低，这又会反过来进一步增加邻居值，导致错误的放大。

当然，通过减小网格尺度让单元Peclet数小于2，就可以避免这种情况。然而，对于很多实际情况，这样做的话将会显著增加存储量和计算量，以至超出计算机负荷。而且，对于纯对流的问题（即不考虑粘性扩散效应，例如，Euler流动），这么做的话也是不可行的，所以需要发展修正方法。

2.4 迎风格式（Upwind Scheme）

如上图，与对流过程更加契合的格式为迎风格式，迎风格式基本上是模仿了对流的物理特性，即，单元面上的值是依赖于迎风节点上的值的，换言之，依赖于流动方向的。图中所示的单元面上的值用下式给出
$${\varphi }_e=\{\begin{array}{llc}{\varphi }_C& if& {\dot{m}}_e>0\\ {\varphi }_E& if& {\dot{m}}_e<0\end{array}{\varphi }_w=\{\begin{array}{llc}{\varphi }_C& if& {\dot{m}}_w>0\\ {\varphi }_W& if& {\dot{m}}_w<0\end{array}$$
（注意，这里的${\varphi }_w$没有写反，因为质量流量$\dot{m}$跟速度$u$是不同的，看下面说明）

其中${\dot{m}}_e$和${\dot{m}}_w$是在面$e$和$w$的质量流量（面积矢量的正方向是朝外的，而非指向单元中心$C$的）
$$\begin{array}{rl}& {\dot{m}}_e=(\rho \mathbf{v}\cdot \mathbf{S}{)}_e=(\rho uS{)}_e=(\rho u\Delta y{)}_e\\ & {\dot{m}}_w=(\rho \mathbf{v}\cdot \mathbf{S}{)}_w=-(\rho uS{)}_w=-(\rho u\Delta y{)}_w\end{array}$$
如此，面$e$的对流通量可写成
$$\begin{array}{rl}{\dot{m}}_e{\varphi }_e& =||{\dot{m}}_e,0||{\varphi }_C-||-{\dot{m}}_e,0||{\varphi }_E\\ & =Flux{C}_e^{Conv}{\varphi }_C+Flux{F}_e^{Conv}{\varphi }_E+Flux{V}_e^{Conv}\end{array}$$
其中
$$\begin{array}{rl}& Flux{C}_e^{Conv}=||{\dot{m}}_e,0||\\ & Flux{F}_e^{Conv}=-||-{\dot{m}}_e,0||\\ & Flux{V}_e^{Conv}=0\end{array}$$
其中的$||a,b||=\max (a,b)$代表了$a$和$b$中的最大值。同样手法，也可推导面$w$的对流通量
$$\begin{array}{rl}{\dot{m}}_w{\varphi }_w& =||{\dot{m}}_w,0||{\varphi }_C-||-{\dot{m}}_w,0||{\varphi }_W\\ & =Flux{C}_w^{Conv}{\varphi }_C+Flux{F}_w^{Conv}{\varphi }_W+Flux{V}_w^{Conv}\end{array}$$
其中
$$\begin{array}{rl}& Flux{C}_w^{Conv}=||{\dot{m}}_w,0||\\ & Flux{F}_w^{Conv}=-||-{\dot{m}}_w,0||\\ & Flux{V}_w^{Conv}=0\end{array}$$
把扩散项用上标$Diff$表示，推得一维稳态对流扩散方程的离散格式为
$$\begin{array}{rl}(Flux{C}_e^{Conv}& +Flux{C}_e^{Diff}+Flux{C}_w^{Conv}+Flux{C}_w^{Diff}){\varphi }_C\\ & +(Flux{F}_e^{Conv}+Flux{F}_e^{Diff}){\varphi }_E+(Flux{F}_w^{Conv}+Flux{F}_w^{Diff}){\varphi }_W=0\end{array}$$
可写成下面形式
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W=b_C$$
其中
$$\begin{array}{rl}{a}_E& =Flux{F}_e^{Conv}+Flux{F}_e^{Diff}\\ & =-||-{\dot{m}}_e,0||-{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}\\ a_W& =Flux{F}_w^{Conv}+Flux{F}_w^{Diff}\\ & =-||-{\dot{m}}_w,0||-{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}\\ a_C& =\sum _f\left(Flux{C}_f^{Conv}+Flux{C}_f^{Diff}\right)\\ & =||{\dot{m}}_e,0||+||{\dot{m}}_w,0||+{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}+{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}\\ & =-(a_E+a_W)+\underset{=0}{\underbrace{({\dot{m}}_e+{\dot{m}}_w)}}\\ b_C& =-\sum _f\left(Flux{V}_f^{Conv}+Flux{V}_f^{Diff}\right)\\ & =0\end{array}$$
显而易见，迎风格式推出了负的邻居系数，而且在满足连续方程的情况下（即，${\dot{m}}_e+{\dot{m}}_w=0$），主节点上的系数为
$$a_C=-(a_W+a_E)$$
这确保了有界性。

在连续方程满足的前提下，假设为均匀网格和常扩散系数，用${\varphi }_E$和${\varphi }_W$所给出的${\varphi }_C$值为
$$\frac{{\varphi }_C-{\varphi }_W}{{\varphi }_E-{\varphi }_W}=\frac{2+||-P{e}_L,0||}{4+||-P{e}_L,0||+||P{e}_L,0||}=\frac{2+||-P{e}_L,0||}{4+|P{e}_L|}$$
将该迎风格式的数值解与精确的解析解、之前中心差分格式的数值解一并展示于下图中。当$P{e}_L$的值（绝对值）较小时，迎风格式并没有中心差分格式的结果准确，这是因为迎风廓线只有一阶精度，而线性廓线具有二阶精度；而当$P{e}_L$的值（绝对值）较大时，中心差分格式是不稳定的，其表现为无界性，已经失去了物理意义，而此时的迎风格式，尽管并不是特别准确，但是物理意义上是正确的。

如此一来，就有了精确性和稳定性之间的权衡。使用迎风格式，得到的解即便在高$P{e}_L$下依旧保持有界性和有物理意义，然而付出的代价就是精度较低。另一方面，二阶精度的中心差分格式当$P{e}_L$超出特定阈值的时候会变得不稳定，产生不符合物理意义的解。两种格式都受到了误差的影响，一个影响了精确性，另一个影响了稳定性。这些都是什么误差呢？在引入背风（downwind）格式后，咱们再讨论这个问题。

2.5 背风格式（Downwind Scheme）

如果把迎风格式反过来用，即，用背风格式的话，会发生什么有趣的事情呢？如上图所示，在这种插值廓线中，是把界面背风方向节点处的值视作界面上的值，这样，面$e$和面$w$处的值为
$${\varphi }_e=\{\begin{array}{llc}{\varphi }_E& if& {\dot{m}}_e>0\\ {\varphi }_C& if& {\dot{m}}_e<0\end{array}{\varphi }_w=\{\begin{array}{llc}{\varphi }_W& if& {\dot{m}}_w>0\\ {\varphi }_C& if& {\dot{m}}_w<0\end{array}$$
面上的对流通量可以写成
$$\begin{array}{rl}{\dot{m}}_e{\varphi }_e& =-||-{\dot{m}}_e,0||{\varphi }_C+||{\dot{m}}_e,0||{\varphi }_E\\ & =Flux{C}_e^{Conv}{\varphi }_C+Flux{F}_e^{Conv}{\varphi }_E+Flux{V}_e^{Conv}\\ {\dot{m}}_w{\varphi }_w& =-||-{\dot{m}}_w,0||{\varphi }_C+||{\dot{m}}_w,0||{\varphi }_W\\ & =Flux{C}_w^{Conv}{\varphi }_C+Flux{F}_w^{Conv}{\varphi }_W+Flux{V}_w^{Conv}\end{array}$$
可进一步推得一维稳态对流扩散方程的离散形式为
$$\begin{array}{rl}(Flux{C}_e^{Conv}& +Flux{C}_e^{Diff}+Flux{C}_w^{Conv}+Flux{C}_w^{Diff}){\varphi }_C\\ & +(Flux{F}_e^{Conv}+Flux{F}_e^{Diff}){\varphi }_E+(Flux{F}_w^{Conv}+Flux{F}_w^{Diff}){\varphi }_W=0\end{array}$$
可写成下面形式
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W=b_C$$
其中
$$\begin{array}{rl}{a}_E& =Flux{F}_e^{Conv}+Flux{F}_e^{Diff}\\ & =||{\dot{m}}_e,0||-{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}\\ a_W& =Flux{F}_w^{Conv}+Flux{F}_w^{Diff}\\ & =||{\dot{m}}_w,0||-{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}\\ a_C& =\sum _f\left(Flux{C}_f^{Conv}+Flux{C}_f^{Diff}\right)\\ & =-||-{\dot{m}}_e,0||-||-{\dot{m}}_w,0||+{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}+{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}\\ & =-(a_E+a_W)+\underset{=0}{\underbrace{({\dot{m}}_e+{\dot{m}}_w)}}\\ b_C& =-\sum _f\left(Flux{V}_f^{Conv}+Flux{V}_f^{Diff}\right)\\ & =0\end{array}$$
在连续方程满足的前提下，假设为均匀网格和常扩散系数，用${\varphi }_E$和${\varphi }_W$所给出的${\varphi }_C$值为
$$\frac{{\varphi }_C-{\varphi }_W}{{\varphi }_E-{\varphi }_W}=\frac{2-||P{e}_L,0||}{4-||-P{e}_L,0||-||P{e}_L,0||}=\frac{2-||P{e}_L,0||}{4-|P{e}_L|}$$
不用把上式的图形画出来，也可以看出，当$|P{e}_L|\to 4$的时候，解变得完全无界。背风格式与其它格式联合起来预测突变界面时是有用的，然而这里引入背风格式的目的主要是为了便于在下一节给出更好的稳定性分析。

3 截断误差：数值扩散和反扩散（Truncation Error: Numerical Diffusion and Anti-Diffusion）

截断误差是由离散过程的自然近似所造成的，对于笛卡尔网格上的一维情形较容易分析。下面将做迎风、背风、中心差分格式的扩散和反扩散分析。

3.1 迎风格式

针对通过迎风格式离散得到的方程，尝试将其恢复到最初的积分方程并考察其截断误差。假设流动是沿着$x$正方向的，将${\varphi }_C$和${\varphi }_W$分别写成${\varphi }_e$和${\varphi }_w$函数的形式，迎风格式为
$${\varphi }_e={\varphi }_C{\varphi }_w={\varphi }_W$$
这样，使用迎风格式的一维对流扩散方程离散形式为
$$(\rho u\Delta y{)}_e{\varphi }_C-(\rho u\Delta y{)}_w{\varphi }_W-\left[{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]=0$$
将${\varphi }_C$做一维Taylor展开，以单元面$e$的值为基准，得
$$\begin{array}{rl}{\varphi }_C& ={\varphi }_e+{\left(\frac{d\varphi }{dx}\right)}_e(x_C-x_e)+\frac{1}{2}{\left(\frac{d^2\varphi }{d{x}^2}\right)}_e(x_C-x_e{)}^2+...\\ & ={\varphi }_e-{\left(\frac{d\varphi }{dx}\right)}_e(x_e-x_C)+...\end{array}$$
对于均匀网格
$${\varphi }_C={\varphi }_e-{\left(\frac{d\varphi }{dx}\right)}_e\frac{(\delta x{)}_e}{2}+\frac{1}{2}{\left(\frac{d^2\varphi }{d{x}^2}\right)}_e{\left(\frac{(\delta x{)}_e}{2}\right)}^2+...$$
可获得${\varphi }_W$的相似表达式
$${\varphi }_W={\varphi }_w-{\left(\frac{d\varphi }{dx}\right)}_w\frac{(\delta x{)}_w}{2}+\frac{1}{2}{\left(\frac{d^2\varphi }{d{x}^2}\right)}_w{\left(\frac{(\delta x{)}_w}{2}\right)}^2+...$$
截去二阶项和高阶项，将其代入到对流项中，离散方程的左端项变为
$$\begin{array}{rl}& (\rho u\Delta y{)}_e{\varphi }_C-(\rho u\Delta y{)}_w{\varphi }_W-\left[{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]=\\ & (\rho u\Delta y{)}_e\left[{\varphi }_e-{\left(\frac{d\varphi }{dx}\right)}_e\frac{(\delta x{)}_e}{2}\right]-(\rho u\Delta y{)}_w\left[{\varphi }_w-{\left(\frac{d\varphi }{dx}\right)}_w\frac{(\delta x{)}_w}{2}\right]-\\ & \left[{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]\end{array}$$
可以改写为
$$\begin{array}{rl}& (\rho u\Delta y{)}_e{\varphi }_C-(\rho u\Delta y{)}_w{\varphi }_W-\left[{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]=\\ & (\rho u\Delta y{)}_e{\varphi }_e-(\rho u\Delta y{)}_w{\varphi }_w-\left[{\left({\Gamma }^{\varphi }+\rho u\frac{\delta x}{2}\right)}_e{\left(\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }+\rho u\frac{\delta x}{2}\right)}_w{\left(\frac{d\varphi }{dx}\Delta y\right)}_w\right]\end{array}$$
显而易见，待求解的方程中额外加入了扩散项分量，这常被称为截断误差。数值扩散的值为
$${\Gamma }_{truncation}^{\varphi }=\rho u\frac{\delta x}{2}$$
该截断误差，也称为流向扩散，通过更改扩散系数的值影响了待求解方程，从而降低了解的精确性。这样，相当于改变了对流扩散方程的扩散效应。另一方面，该额外的流向扩散效应也是有益的，因为其对解进行了限定，起到了稳定作用，可得到物理上有意义的解。

显然，要想减少这个流向数值扩散，需要针对对流项构造高阶精度的近似格式。然而，后面章节中会讲到，这些格式在构造的时候还需要保证解的有界性。

3.2 背风格式

跟上一小节类似，假设流动是沿着$x$正方向的，将${\varphi }_E$和${\varphi }_C$分别写成${\varphi }_e$和${\varphi }_w$函数的形式，背风格式为
$${\varphi }_e={\varphi }_E{\varphi }_w={\varphi }_C$$
这样，使用背风格式的一维对流扩散方程离散形式为
$$(\rho u\Delta y{)}_e{\varphi }_E-(\rho u\Delta y{)}_w{\varphi }_C-\left[{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]=0$$
在均匀网格上，将${\varphi }_E$和${\varphi }_C$做一维Taylor展开，分别以单元面$e$和$w$的值为基准，得
$$\begin{array}{rl}& {\varphi }_E={\varphi }_e+{\left(\frac{d\varphi }{dx}\right)}_e\frac{(\delta x{)}_e}{2}+\frac{1}{2}{\left(\frac{d^2\varphi }{d{x}^2}\right)}_e{\left(\frac{(\delta x{)}_e}{2}\right)}^2+...\\ & {\varphi }_C={\varphi }_w+{\left(\frac{d\varphi }{dx}\right)}_w\frac{(\delta x{)}_w}{2}+\frac{1}{2}{\left(\frac{d^2\varphi }{d{x}^2}\right)}_w{\left(\frac{(\delta x{)}_w}{2}\right)}^2+...\end{array}$$
截去二阶项和高阶项，将其代入到对流项中，离散方程的左端项变为
$$\begin{array}{rl}& (\rho u\Delta y{)}_e{\varphi }_E-(\rho u\Delta y{)}_w{\varphi }_C-\left[{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }\frac{d\varphi }{dx}\Delta y\right)}_w\right]=\\ & (\rho u\Delta y{)}_e{\varphi }_e-(\rho u\Delta y{)}_w{\varphi }_w-\left[{\left({\Gamma }^{\varphi }-\rho u\frac{\delta x}{2}\right)}_e{\left(\frac{d\varphi }{dx}\Delta y\right)}_e-{\left({\Gamma }^{\varphi }-\rho u\frac{\delta x}{2}\right)}_w{\left(\frac{d\varphi }{dx}\Delta y\right)}_w\right]\end{array}$$
那么对于背风格式，其推出的离散方程中的数值扩散是负的，等于
$${\Gamma }_{truncation}^{\varphi }=-\rho u\frac{\delta x}{2}$$
这作用是减小扩散系数，这会导致反扩散错误。使用背风格式的结果会造成clipping of the advected profiles（没搞懂啥意思）。实际上使用背风格式得到的一维对流扩散方程的解比中心差分格式的振荡还要剧烈。

3.3 中心差分格式

有限体积方法中，对中心差分格式的截断误差分析略嫌麻烦，因为梯度的计算需要用到单元面上的插值，而并非直接用到节点处的已有值。假设各处的速度是已知的，且网格是均匀的，网格尺度为$\Delta x$，那么$({\varphi }_e-{\varphi }_w)$为
$$\underset{Interpolated}{\underbrace{({\varphi }_e-{\varphi }_w)}}=\frac{1}{2}({\varphi }_E+{\varphi }_C)-\frac{1}{2}({\varphi }_C+{\varphi }_W)=\underset{Exact}{\underbrace{({\varphi }_e-{\varphi }_w)}}+TE$$
式中的$TE$代表截断误差，为了算出截断误差，做Taylor展开，基于精确的${\varphi }_e$和${\varphi }_w$做展开，得
$$\begin{array}{rl}& {\varphi }_W={\varphi }_w-\frac{\Delta x}{2}{\varphi }_w^{\prime }+\frac{\Delta x^2}{8}{\varphi }_w^{\prime \prime }-\frac{\Delta x^3}{48}{\varphi }_w^{\prime \prime \prime }+\frac{\Delta x^4}{384}{\varphi }_w^{iv}-\frac{\Delta x^5}{3840}{\varphi }_w^v+...\\ & {\varphi }_C={\varphi }_w+\frac{\Delta x}{2}{\varphi }_w^{\prime }+\frac{\Delta x^2}{8}{\varphi }_w^{\prime \prime }+\frac{\Delta x^3}{48}{\varphi }_w^{\prime \prime \prime }+\frac{\Delta x^4}{384}{\varphi }_w^{iv}+\frac{\Delta x^5}{3840}{\varphi }_w^v+...\\ & {\varphi }_C={\varphi }_e-\frac{\Delta x}{2}{\varphi }_e^{\prime }+\frac{\Delta x^2}{8}{\varphi }_e^{\prime \prime }-\frac{\Delta x^3}{48}{\varphi }_e^{\prime \prime \prime }+\frac{\Delta x^4}{384}{\varphi }_e^{iv}-\frac{\Delta x^5}{3840}{\varphi }_e^v+...\\ & {\varphi }_E={\varphi }_e+\frac{\Delta x}{2}{\varphi }_e^{\prime }+\frac{\Delta x^2}{8}{\varphi }_e^{\prime \prime }+\frac{\Delta x^3}{48}{\varphi }_e^{\prime \prime \prime }+\frac{\Delta x^4}{384}{\varphi }_e^{iv}+\frac{\Delta x^5}{3840}{\varphi }_e^v+...\end{array}$$
使用这些展开形式，面上的平均值为
$$\begin{array}{rl}\frac{1}{2}({\varphi }_E+{\varphi }_C)& ={\varphi }_e+\frac{\Delta x^2}{8}{\varphi }_e^{\prime \prime }+\frac{\Delta x^4}{384}{\varphi }_e^{iv}+...\\ \Rightarrow \frac{\Delta x^2}{4}{\varphi }_e^{\prime \prime }& =({\varphi }_C-2{\varphi }_e+{\varphi }_E)-\frac{\Delta x^4}{192}{\varphi }_e^{iv}+...\\ \frac{1}{2}({\varphi }_C+{\varphi }_W)& ={\varphi }_w+\frac{\Delta x^2}{8}{\varphi }_w^{\prime \prime }+\frac{\Delta x^4}{384}{\varphi }_w^{iv}+...\\ \Rightarrow \frac{\Delta x^2}{4}{\varphi }_w^{\prime \prime }& =({\varphi }_W-2{\varphi }_w+{\varphi }_C)-\frac{\Delta x^4}{192}{\varphi }_w^{iv}+...\end{array}$$
上两式相减，得到
$$\frac{1}{2}({\varphi }_E+{\varphi }_C)-\frac{1}{2}({\varphi }_C+{\varphi }_W)=({\varphi }_e-{\varphi }_w)+\frac{\Delta x^2}{8}({\varphi }_e^{\prime \prime }-{\varphi }_w^{\prime \prime })+\frac{\Delta x^4}{384}({\varphi }_e^{iv}-{\varphi }_w^{iv})+...$$
那么，如果把${\varphi }_e$和${\varphi }_w$的精确值用${\varphi }_C$展开，为
$$\begin{array}{rl}& {\varphi }_e={\varphi }_C+\frac{\Delta x}{2}{\varphi }_C^{\prime }+\frac{\Delta x^2}{8}{\varphi }_C^{\prime \prime }+\frac{\Delta x^3}{48}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^4}{384}{\varphi }_C^{iv}+\frac{\Delta x^5}{3840}{\varphi }_C^v+...\\ & {\varphi }_w={\varphi }_C-\frac{\Delta x}{2}{\varphi }_C^{\prime }+\frac{\Delta x^2}{8}{\varphi }_C^{\prime \prime }-\frac{\Delta x^3}{48}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^4}{384}{\varphi }_C^{iv}-\frac{\Delta x^5}{3840}{\varphi }_C^v+...\end{array}$$
此外，${\varphi }_E$和${\varphi }_W$也可以用${\varphi }_C$展开为
$$\begin{array}{rl}& {\varphi }_E={\varphi }_C+\Delta x{\varphi }_C^{\prime }+\frac{\Delta x^2}{2}{\varphi }_C^{\prime \prime }+\frac{\Delta x^3}{6}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^4}{24}{\varphi }_C^{iv}+\frac{\Delta x^5}{120}{\varphi }_C^v+...\\ & {\varphi }_W={\varphi }_C-\Delta x{\varphi }_C^{\prime }+\frac{\Delta x^2}{2}{\varphi }_C^{\prime \prime }-\frac{\Delta x^3}{6}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^4}{24}{\varphi }_C^{iv}-\frac{\Delta x^5}{120}{\varphi }_C^v+...\end{array}$$
从而得到
$$\begin{array}{rl}{\varphi }_C-2{\varphi }_e+{\varphi }_E=& {\varphi }_C-2{\varphi }_C-\Delta x{\varphi }_C^{\prime }-\frac{\Delta x^2}{4}{\varphi }_C^{\prime \prime }-\frac{\Delta x^3}{24}{\varphi }_C^{\prime \prime \prime }-\frac{\Delta x^4}{192}{\varphi }_C^{iv}-\frac{\Delta x^5}{1920}{\varphi }_C^v+...\\ & +{\varphi }_C+\Delta x{\varphi }_C^{\prime }+\frac{\Delta x^2}{2}{\varphi }_C^{\prime \prime }+\frac{\Delta x^3}{6}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^4}{24}{\varphi }_C^{iv}+\frac{\Delta x^5}{120}{\varphi }_C^v+...\\ =& \frac{\Delta x^2}{4}{\varphi }_C^{\prime \prime }+\frac{\Delta x^3}{8}{\varphi }_C^{\prime \prime \prime }+\frac{7\Delta x^4}{192}{\varphi }_C^{iv}+\frac{15\Delta x^5}{1920}{\varphi }_C^v+...\\ {\varphi }_W-2{\varphi }_w+{\varphi }_C=& \frac{\Delta x^2}{4}{\varphi }_C^{\prime \prime }-\frac{\Delta x^3}{8}{\varphi }_C^{\prime \prime \prime }+\frac{7\Delta x^4}{192}{\varphi }_C^{iv}-\frac{15\Delta x^5}{1920}{\varphi }_C^v+...\end{array}$$
得到
$$\begin{array}{rl}\frac{\Delta x^2}{8}({\varphi }_e^{\prime \prime }-{\varphi }_w^{\prime \prime })& =\frac{1}{2}[({\varphi }_C-2{\varphi }_e+{\varphi }_E)-({\varphi }_W-2{\varphi }_w+{\varphi }_C)]+...\\ & =\frac{\Delta x^3}{8}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^5}{128}{\varphi }_C^v+...\end{array}$$
代入到前式$\frac{1}{2}({\varphi }_E+{\varphi }_C)-\frac{1}{2}({\varphi }_C+{\varphi }_W)=({\varphi }_e-{\varphi }_w)+\frac{\Delta x^2}{8}({\varphi }_e^{\prime \prime }-{\varphi }_w^{\prime \prime })+\frac{\Delta x^4}{384}({\varphi }_e^{iv}-{\varphi }_w^{iv})+...$，得
$$\frac{1}{2}({\varphi }_E+{\varphi }_C)-\frac{1}{2}({\varphi }_C+{\varphi }_W)=({\varphi }_e-{\varphi }_w)+\frac{\Delta x^3}{8}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^5}{128}{\varphi }_C^v+...$$
上式两端同除以$\Delta x$，得到计算梯度时的截断误差为
$$TE=\frac{\Delta x^2}{8}{\varphi }_C^{\prime \prime \prime }+\frac{\Delta x^4}{128}{\varphi }_C^v+...$$
这表明该格式具有二阶精度。

4 数值稳定性

在前面的截断误差分析中，一阶精度格式能得到符合物理意义的解，二阶精度反而会得到毫无意义的解，这着实很让人迷惑不解。中心差分格式对于扩散项是如此精确，而对于对流项应该也差不多精确才是啊。然而，前面咱们已经讲了，中心差分格式针对对流项会产生毫无物理意义的解。产生这种现象的原因就在于，当用于推导奇数阶（odd order）比如对流项时，中心差分格式并没有继承对流稳定性（convective stability）。

为了解释中心差分格式应用于对流占优问题时所出现的数值解的振荡问题，Leonard提出了对流稳定性的概念。Leonard使用通有的一维非定常对流扩散方程，且速度为常数，来描述该问题
$$\frac{\partial (\rho \varphi )}{\partial t}=-\frac{\partial (\rho u\varphi )}{\partial x}+\frac{\partial }{\partial x}\left({\Gamma }^{\varphi }\frac{\partial \varphi }{\partial x}\right)+Q^{\varphi }$$

当把这个方程应用于上图所示的单元$C$时，其左端项（LHS）代表了在控制单元上单位时间内${\varphi }_C$的变化率，其右端项（RHS）代表着通过单元面的净流入通量和单元中的源项，它们影响着${\varphi }_C$的值。如果RHS有数值误差，那么算出来的${\varphi }_C$要么增加要么减小，这得视离散过程采用的离散格式而定。在不稳定的格式中，对正确值${\varphi }_C$的微小偏差将会造成RHS所代表的净流入通量的增加或减小。当使用迭代算法来求解时，净流入通量的增加或减小将通过每个迭代步进一步增加或减小${\varphi }_C$的值。而在稳定格式中，${\varphi }_C$的变化导致RHS的误差将会负反馈到RHS，其具有自我修正的特性。显然，对于这种类型的数值稳定性，RHS应该满足
$$\frac{\partial (RHS)}{\partial {\varphi }_C}<0$$
这表明，方程右端项的离散格式应能保证其对${\varphi }_C$的敏感性是负的，这样才能让${\varphi }_C$的增加或减小造成净流入通量的减小或增加，反过来使得${\varphi }_C$减小或增加来趋近其正确值。然而，稳定性不能和有界性、精确性混淆起来。稳定的数值格式也有可能是无界限的，产生过冲和下冲（over/under shoots），产生振荡和起伏，产生过度扩散，也可能给出很低精度的解。这里的稳定性指的是把数值误差控制在有界范围内不让它无穷增大，就跟使用中心差分格式时出现的${\varphi }_C$的无量纲值随$Pe$从$-\infty$到$+\infty$的变化而从$-\infty$变化到$+\infty$那样，实际上${\varphi }_C$的正确值应该是从$1$变化到$0$。当然，迎风格式既拥有有界性也拥有稳定性。

基于这种思想，给出RHS的通有离散形式
$$\begin{array}{rl}RHS=& -(\rho u\Delta y{)}_e{\varphi }_e+(\rho u\Delta y{)}_w{\varphi }_w\\ & +\left[{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]+Q_C^{\varphi }{V}_C\end{array}$$
使用中心差分格式，上式变为
$$\begin{array}{rl}RH{S}_{CD}=& -(\rho u\Delta y{)}_e\frac{{\varphi }_E+{\varphi }_C}{2}+(\rho u\Delta y{)}_w\frac{{\varphi }_C+{\varphi }_W}{2}\\ & +\left[{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]+Q_C^{\varphi }{V}_C\end{array}$$
使用前述判据来分析中心差分格式，发现扩散项的敏感性为
$$\frac{\partial (RH{S}_{CD}^{Diff})}{\partial {\varphi }_C}=-2{\Gamma }^{\varphi }\frac{\Delta y}{\delta x}$$
因为${\Gamma }^{\varphi }$是正值，所以上式值为负的，表明是稳定格式。然而，对于对流项，其敏感性方程为
$$\frac{\partial (RH{S}_{CD}^{Conv})}{\partial {\varphi }_C}=-\frac{1}{2}({\dot{m}}_e+{\dot{m}}_w)$$
对于稳态问题，上式为0，但是对于非稳态问题则未必是0，对于非稳态情形，减速流动时其值是正的。在一般的流动中，这些区域作为起伏的区域，在Peclet数过大的时候很容易就造成了数值崩溃。即便是对于稳态流动，由于其值为0，所以无法很好地负向反馈给方程来做自修正处理。上式同时表明，对于稳态问题，CD格式算得的净对流通量是和${\varphi }_C$值无关的。因此，即便是不同的${\varphi }_C$值，算得的单元$C$的对流通量也是一样的。

接下来看迎风格式，其推得的RHS为
$$\begin{array}{rl}RH{S}_{Upwind}=& -||{\dot{m}}_e,0||{\varphi }_C+||-{\dot{m}}_e,0||{\varphi }_E-||{\dot{m}}_w,0||{\varphi }_C+||-{\dot{m}}_w,0||{\varphi }_W\\ & +\left[{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]+Q_C^{\varphi }{V}_C\end{array}$$
跟预想的一样，其敏感性是负的，因为这个迎风格式对于所有的Peclet数都是稳定的，即
$$\frac{\partial (RH{S}_{Upwind}^{Conv})}{\partial {\varphi }_C}=-||{\dot{m}}_e,0||-||{\dot{m}}_w,0||$$
其为负值或0值（对于所有流动情况），其等于0的情况是，两个质量流量都是负的，这种情况在一维定截面问题中不会出现。再考虑到一阶近似所带来的假扩散，上式表明迎风格式是非常稳定的。然而，其稳定性的实现是以精确性为代价的。

对于背风格式，RHS为
$$\begin{array}{rl}RH{S}_{Downwind}=& -||{\dot{m}}_e,0||{\varphi }_E+||-{\dot{m}}_e,0||{\varphi }_C-||{\dot{m}}_w,0||{\varphi }_W+||-{\dot{m}}_w,0||{\varphi }_C\\ & +\left[{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{\Delta y}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]+Q_C^{\varphi }{V}_C\end{array}$$
其敏感性为
$$\frac{\partial (RH{S}_{Upwind}^{Conv})}{\partial {\varphi }_C}=||-{\dot{m}}_e,0||+||-{\dot{m}}_w,0||$$
其总是正值或者0值（对于所有流动情况），再考虑到反扩散效应，背风格式是高度的不稳定格式！

5 高阶迎风格式

前面小节中讲到，迎风格式和中心差分格式都有严重的缺陷，前者由于数值扩散使得其精度较差，后者的稳定性差存在数值dispersion（色散）误差。这些缺点促使学者们，基于高阶偏迎风插值廓线，提出精确和稳定的对流项格式。这些高阶格式旨在发展至少具有二阶精度的解，同时保持无条件的稳定性。

为了强调单元面通量的守恒特性，将用$D、C、U$表示某特定面的背风Downwind、迎风Upwind、远迎风Far Upwind节点，用$DD$和$UU$表示$D$的下游downstream节点和$U$的上游upstream节点，对应的值用${\varphi }_{DD},{\varphi }_D,{\varphi }_C,{\varphi }_U,{\varphi }_{UU}$来表示。如下图所示，图中给出了面上的速度为正（$\to$）和负（$\leftarrow$）两种情况。

5.1 二阶迎风格式（Second Order Upwind Scheme SOU）

二阶格式使用的仍旧是线性廓线，跟在中心差分格式中的廓线一样，然而，与对称廓线不同，其使用的是偏迎风架构的廓线（即其架构点在迎风方向比背风方向要更多一些）。如下图所示，通过在节点$C$和$U$处的$\varphi$值来构造线性廓线，因此，面上的值实际上是通过外插而非内插来算得的。

5.2 插值廓线

从线性廓线开始
$$\varphi (x)=k_0+k_1(x-x_C)$$
将其用于$x_C$和$x_U$节点的$\varphi$值${\varphi }_C$和${\varphi }_U$，廓线变为
$$\varphi (x)={\varphi }_C+\frac{{\varphi }_C-{\varphi }_U}{x_C-x_U}(x-x_C)$$
使用上式，面$f$处的$\varphi$值，为
$${\varphi }_f=\varphi (x_f)={\varphi }_C+\frac{{\varphi }_C-{\varphi }_U}{x_C-x_U}(x_f-x_C)$$
对于均匀网格，为
$${\varphi }_f=\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_U$$

5.3 离散方程

将这种近似界面值廓线的方式应用于离散一维对流扩散方程，可获得面通量为
$$\begin{array}{rl}{\dot{m}}_e{\varphi }_e& =\left(\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_W\right)||{\dot{m}}_e,0||-\left(\frac{3}{2}{\varphi }_E-\frac{1}{2}{\varphi }_{EE}\right)||-{\dot{m}}_e,0||\\ {\dot{m}}_w{\varphi }_w& =\left(\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_E\right)||{\dot{m}}_w,0||-\left(\frac{3}{2}{\varphi }_W-\frac{1}{2}{\varphi }_{WW}\right)||-{\dot{m}}_w,0||\end{array}$$
将上式代入到半离散通有形式中，得
$$\begin{array}{rl}& \left(\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_W\right)||{\dot{m}}_e,0||-\left(\frac{3}{2}{\varphi }_E-\frac{1}{2}{\varphi }_{EE}\right)||-{\dot{m}}_e,0||+\left(\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_E\right)||{\dot{m}}_w,0||\\ & -\left(\frac{3}{2}{\varphi }_W-\frac{1}{2}{\varphi }_{WW}\right)||-{\dot{m}}_w,0||-\left[{\left({\Gamma }^{\varphi }\frac{S}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{S}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]=0\end{array}$$
可写成如下形式
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W+a_{EE}{\varphi }_{EE}+a_{WW}{\varphi }_{WW}=0$$
其中
$$\begin{array}{rl}{a}_E& =Flux{F}_e=-{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}-\frac{3}{2}||-{\dot{m}}_e,0||-\frac{1}{2}||{\dot{m}}_w,0||\\ a_W& =Flux{F}_w=-{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}-\frac{3}{2}||-{\dot{m}}_w,0||-\frac{1}{2}||{\dot{m}}_e,0||\\ a_{EE}& =Flux{F}_{ee}=\frac{1}{2}||-{\dot{m}}_e,0||\\ a_{WW}& =Flux{F}_{ww}=\frac{1}{2}||-{\dot{m}}_w,0||\\ a_C& =\sum _{f\sim nb(C)}Flux{C}_f\\ & ={\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}+{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}+\frac{3}{2}||{\dot{m}}_e,0||+\frac{3}{2}||{\dot{m}}_w,0||\\ & =-(a_E+a_W+a_{EE}+a_{WW})+({\dot{m}}_e+{\dot{m}}_w)\end{array}$$

5.4 截断误差

跟前面中心差分格式的分析过程类似，可推得二阶迎风格式的截断误差为
$$TE=-\frac{3}{8}\Delta x^2{\varphi }_C^{\prime \prime \prime }-\frac{1}{4}\Delta x^3{\varphi }_C^{iv}+...$$
表明其为二阶精度。

5.5 稳定性分析

跟前面的稳定性分析过程类似，推得离散对流项为
$$\begin{array}{rl}RH{S}_{convection}=& -\left(\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_W\right)||{\dot{m}}_e,0||+\left(\frac{3}{2}{\varphi }_E-\frac{1}{2}{\varphi }_{EE}\right)||-{\dot{m}}_e,0||\\ & -\left(\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_E\right)||{\dot{m}}_w,0||+\left(\frac{3}{2}{\varphi }_W-\frac{1}{2}{\varphi }_{WW}\right)||-{\dot{m}}_w,0||\end{array}$$
该项对${\varphi }_C$的变化率为
$$\frac{\partial (RH{S}_{convection})}{\partial {\varphi }_C}=-\frac{3}{2}||{\dot{m}}_e,0||-\frac{3}{2}||-{\dot{m}}_w,0||$$
该值总是负的，表明这是稳定的格式，即，解的误差将在可控范围内。需要注意的是，这只是针对推导过程中假设的条件才得到的结论，并不是针对一般的情况，比如速度不为常数的时候就未必如此了。

5.6 QUICK格式

对流运动学的平方迎风插值（Quadratic Upstream Interpolation for Convective Kinematics QUICK）格式是由Leonard提出的，该方法的思想是，每个面上值在插值计算时，用偏向迎风方向的平方多项式近似，如下图所示。

5.7 插值廓线

对于一维对流扩散方程，$\varphi$值用下式计算
$$\varphi (x)=k_0+k_{1}x+k_2{x}^2$$
有
$$\varphi =\{\begin{array}{llc}{\varphi }_U& at& x=x_U\\ {\varphi }_C& at& x=x_C\\ {\varphi }_D& at& x=x_D\end{array}$$
推得
$$\varphi ={\varphi }_U+\frac{(x-x_U)(x-x_C)}{(x_D-x_U)(x_D-x_C)}({\varphi }_D-{\varphi }_U)+\frac{(x-x_U)(x-x_D)}{(x_C-x_U)(x_C-x_D)}({\varphi }_C-{\varphi }_U)$$
如果是均匀网格，那么单元面$f$上的$\varphi$值为
$${\varphi }_f=\frac{{\varphi }_C+{\varphi }_D}{2}-\frac{{\varphi }_D-2{\varphi }_C+{\varphi }_U}{8}=\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_U+\frac{3}{8}{\varphi }_D$$
这样，单元面上的对流通量为
$$\begin{array}{rl}{\dot{m}}_e{\varphi }_e& =\left(\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_W+\frac{3}{8}{\varphi }_E\right)||{\dot{m}}_e,0||-\left(\frac{3}{4}{\varphi }_E-\frac{1}{8}{\varphi }_{EE}+\frac{3}{8}{\varphi }_C\right)||-{\dot{m}}_e,0||\\ {\dot{m}}_w{\varphi }_w& =\left(\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_E+\frac{3}{8}{\varphi }_W\right)||{\dot{m}}_w,0||-\left(\frac{3}{4}{\varphi }_W-\frac{1}{8}{\varphi }_{WW}+\frac{3}{8}{\varphi }_C\right)||-{\dot{m}}_w,0||\end{array}$$
代回到一维对流扩散方程中，其离散形式为
$$\begin{array}{rl}& \left(\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_W+\frac{3}{8}{\varphi }_E\right)||{\dot{m}}_e,0||-\left(\frac{3}{4}{\varphi }_E-\frac{1}{8}{\varphi }_{EE}+\frac{3}{8}{\varphi }_C\right)||-{\dot{m}}_e,0||\\ & +\left(\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_E+\frac{3}{8}{\varphi }_W\right)||{\dot{m}}_w,0||-\left(\frac{3}{4}{\varphi }_W-\frac{1}{8}{\varphi }_{WW}+\frac{3}{8}{\varphi }_C\right)||-{\dot{m}}_w,0||\\ & -\left[{\left({\Gamma }^{\varphi }\frac{S}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{S}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]=0\end{array}$$
可写成如下形式
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W+a_{EE}{\varphi }_{EE}+a_{WW}{\varphi }_{WW}=0$$
其中
$$\begin{array}{rl}{a}_E& =Flux{F}_e=-{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}-\frac{3}{4}||-{\dot{m}}_e,0||+\frac{3}{8}||{\dot{m}}_e,0||-\frac{1}{8}||{\dot{m}}_w,0||\\ a_W& =Flux{F}_w=-{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}-\frac{3}{4}||-{\dot{m}}_w,0||+\frac{3}{8}||{\dot{m}}_w,0||-\frac{1}{8}||{\dot{m}}_e,0||\\ a_{EE}& =Flux{F}_{ee}=\frac{1}{8}||-{\dot{m}}_e,0||\\ a_{WW}& =Flux{F}_{ww}=\frac{1}{8}||-{\dot{m}}_w,0||\\ a_C& =\sum _{f\sim nb(C)}Flux{C}_f\\ & ={\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}+{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}+\frac{3}{4}||{\dot{m}}_e,0||-\frac{3}{8}||-{\dot{m}}_e,0||+\frac{3}{4}||{\dot{m}}_w,0||-\frac{3}{8}||-{\dot{m}}_w,0||\\ & =-(a_E+a_W+a_{EE}+a_{WW})+({\dot{m}}_e+{\dot{m}}_w)\end{array}$$

5.8 截断误差

跟前面一样，推得
$$TE=\frac{1}{16}\Delta x^3{\varphi }_C^{iv}-\frac{3}{128}\Delta x^4{\varphi }_C^v+...$$
显然，有三阶精度。

5.9 稳定性分析

跟前面一样，推得
$$\begin{array}{rl}RH{S}_{convection}=& -\left(\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_W+\frac{3}{8}{\varphi }_E\right)||{\dot{m}}_e,0||+\left(\frac{3}{4}{\varphi }_E-\frac{1}{8}{\varphi }_{EE}+\frac{3}{8}{\varphi }_C\right)||-{\dot{m}}_e,0||\\ & -\left(\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_E+\frac{3}{8}{\varphi }_W\right)||{\dot{m}}_w,0||+\left(\frac{3}{4}{\varphi }_W-\frac{1}{8}{\varphi }_{WW}+\frac{3}{8}{\varphi }_C\right)||-{\dot{m}}_w,0||\end{array}$$
该项对${\varphi }_C$的变化率为
$$\begin{array}{rl}\frac{\partial (RH{S}_{convection})}{\partial {\varphi }_C}& =-\frac{3}{4}||{\dot{m}}_e,0||+\frac{3}{8}||-{\dot{m}}_e,0||-\frac{3}{4}||{\dot{m}}_w,0||+\frac{3}{8}||-{\dot{m}}_w,0||\\ & =-\frac{3}{8}||{\dot{m}}_e,0||-\frac{3}{8}\underset{={\dot{m}}_e}{\underbrace{(||{\dot{m}}_e,0||-||-{\dot{m}}_e,0||)}}-\frac{3}{8}||{\dot{m}}_w,0||-\frac{3}{8}\underset{={\dot{m}}_w}{\underbrace{(||{\dot{m}}_w,0||-||-{\dot{m}}_w,0||)}}\\ & \\ & =-\frac{3}{8}||{\dot{m}}_e,0||-\frac{3}{8}||{\dot{m}}_w,0||-\frac{3}{8}({\dot{m}}_e+{\dot{m}}_w)\end{array}$$
对于均匀速度场来说，其总是一个负值，表明这是稳定的格式。然而在一般的非均匀速度情况下，该格式就未必能保证有界性了。

5.10 FROMM格式

FROMM格式是用界面两侧的远迎风节点（U）和下游节点（D）做线性廓线，如下图所示，与对称廓线不同，其用的是偏迎风的框架来计算面$f$上的值$\varphi$，其假设${\varphi }_U$、${\varphi }_C$和${\varphi }_D$是共线的。

5.11 插值廓线

线性廓线
$$\varphi (x)=k_0+k_1(x-x_C)$$
用在$x_D$和$x_U$处的$\varphi$值${\varphi }_D$和${\varphi }_U$来拟合该廓线，得
$$\varphi (x)={\varphi }_U+\frac{{\varphi }_D-{\varphi }_U}{x_D-x_U}(x-x_U)$$
使用上式，可得迎风节点$C$处的$\varphi$值
$${\varphi }_C={\varphi }_U+\frac{{\varphi }_D-{\varphi }_U}{x_D-x_U}(x_C-x_U)$$
均匀网格上，上式变为
$${\varphi }_C=\frac{{\varphi }_D+{\varphi }_U}{2}$$
面$f$上的值为
$${\varphi }_f=\varphi (f)={\varphi }_U+\frac{{\varphi }_D-{\varphi }_U}{x_D-x_U}(x_f-x_U)={\varphi }_C+\frac{x_f-x_C}{x_D-x_U}({\varphi }_D-{\varphi }_U)$$
均匀网格上，上式变为
$${\varphi }_f={\varphi }_C+\frac{{\varphi }_D-{\varphi }_U}{4}$$

5.12 离散方程

将这种近似界面值廓线的方式应用于离散一维对流扩散方程，可获得面通量为
$$\begin{array}{rl}{\dot{m}}_e{\varphi }_e& =\left({\varphi }_C-\frac{1}{4}{\varphi }_W+\frac{1}{4}{\varphi }_E\right)||{\dot{m}}_e,0||-\left({\varphi }_E-\frac{1}{4}{\varphi }_{EE}+\frac{1}{4}{\varphi }_C\right)||-{\dot{m}}_e,0||\\ {\dot{m}}_w{\varphi }_w& =\left({\varphi }_C-\frac{1}{4}{\varphi }_E+\frac{1}{4}{\varphi }_W\right)||{\dot{m}}_w,0||-\left({\varphi }_W-\frac{1}{4}{\varphi }_{WW}+\frac{1}{4}{\varphi }_C\right)||-{\dot{m}}_w,0||\end{array}$$
将上式代入到半离散通有形式中，得
$$\begin{array}{rl}& \left({\varphi }_C-\frac{1}{4}{\varphi }_W+\frac{1}{4}{\varphi }_E\right)||{\dot{m}}_e,0||-\left({\varphi }_E-\frac{1}{4}{\varphi }_{EE}+\frac{1}{4}{\varphi }_C\right)||-{\dot{m}}_e,0||\\ & +\left({\varphi }_C-\frac{1}{4}{\varphi }_E+\frac{1}{4}{\varphi }_W\right)||{\dot{m}}_w,0||-\left({\varphi }_W-\frac{1}{4}{\varphi }_{WW}+\frac{1}{4}{\varphi }_C\right)||-{\dot{m}}_w,0||\\ & -\left[{\left({\Gamma }^{\varphi }\frac{S}{\delta x}\right)}_e({\varphi }_E-{\varphi }_C)-{\left({\Gamma }^{\varphi }\frac{S}{\delta x}\right)}_w({\varphi }_C-{\varphi }_W)\right]=0\end{array}$$
可写成如下形式
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W+a_{EE}{\varphi }_{EE}+a_{WW}{\varphi }_{WW}=0$$
其中
$$\begin{array}{rl}{a}_E& =Flux{F}_e=-{\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}+\frac{1}{4}||{\dot{m}}_e,0||-||-{\dot{m}}_e,0||-\frac{1}{4}||{\dot{m}}_w,0||\\ a_W& =Flux{F}_w=-{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}+\frac{1}{4}||{\dot{m}}_w,0||-||-{\dot{m}}_w,0||-\frac{1}{4}||{\dot{m}}_e,0||\\ a_{EE}& =Flux{F}_{ee}=\frac{1}{4}||-{\dot{m}}_e,0||\\ a_{WW}& =Flux{F}_{ww}=\frac{1}{4}||-{\dot{m}}_w,0||\\ a_C& =\sum _{f\sim nb(C)}Flux{C}_f\\ & ={\Gamma }_e^{\varphi }\frac{S_e}{\delta x_e}+{\Gamma }_w^{\varphi }\frac{S_w}{\delta x_w}+||{\dot{m}}_e,0||-\frac{1}{4}||-{\dot{m}}_e,0||+||{\dot{m}}_w,0||-\frac{1}{4}||-{\dot{m}}_w,0||\\ & =-(a_E+a_W+a_{EE}+a_{WW})+({\dot{m}}_e+{\dot{m}}_w)\end{array}$$

5.13 截断误差

$$TE=O(\Delta x^2)$$
表明为二阶精度。

5.14 稳定性分析

$$\begin{array}{rl}RH{S}_{convection}=& -\left({\varphi }_C-\frac{1}{4}{\varphi }_W+\frac{1}{4}{\varphi }_E\right)||{\dot{m}}_e,0||+\left({\varphi }_E-\frac{1}{4}{\varphi }_{EE}+\frac{1}{4}{\varphi }_C\right)||-{\dot{m}}_e,0||\\ & -\left({\varphi }_C-\frac{1}{4}{\varphi }_E+\frac{1}{4}{\varphi }_W\right)||{\dot{m}}_w,0||+\left({\varphi }_W-\frac{1}{4}{\varphi }_{WW}+\frac{1}{4}{\varphi }_C\right)||-{\dot{m}}_w,0||\end{array}$$
该项对${\varphi }_C$的变化率为
$$\begin{array}{rl}\frac{\partial (RH{S}_{convection})}{\partial {\varphi }_C}& =-\frac{3}{4}||{\dot{m}}_e,0||-\frac{3}{4}||{\dot{m}}_w,0||-\frac{1}{4}({\dot{m}}_e+{\dot{m}}_w)\end{array}$$
对于定值速度场，其值总是负的，表明这是稳定的格式，但是对于速度并非定值的一般情况则未必如此。

5.15 不同格式的对比

将前面几种格式的稳定性做个对比，明显，负值最大的是二阶迎风SOU格式（系数是-3/2），然后是迎风格式（系数是-1），接下来是FROMM格式（系数为-3/4），再下来是QUICK格式（系数-3/8），最后是中心差分CD格式（系数是0，即，对于${\varphi }_C$的变化是中性稳定的）。前面讲的自修正特性既包含了对流项也包含了扩散项，那么迎风格式的假扩散特性增加了其稳定性，即便其系数只有-1，其也是最稳定的格式。下图展示了不同格式求解对流扩散问题的结果，算例的计算域长度$L=1$，Dirichlet边界条件$x=0$处$\varphi =1$，$x=1$处$\varphi =0$。在较低的Peclet数（Pe=1，如下图（a）所示）下所有的解都是稳定的，在较高的Peclet数（Pe=10，如下图（b）所示）下，CD、FROMM、QUICK格式的解出现了波动。

在低Pe下，CD和QUICK格式的精度相当且与精确解非常接近，而一阶迎风格式的解精度最低，FROMM格式的解比SOU格式的解更为精确，SOU又比迎风格式的精度略高，FROMM和SOU两者的精度稍逊于QUICK格式。即
$$ExactSolution\approx CD\approx QUICK>FROMM>SOU>UpwindatlowPe$$

在高Pe下，各种格式的特性发生了变化。仅仅是迎风格式和SOU格式获得了稳定的解，两者的曲线几乎完全相同，表明SOU格式是高度扩散的。CD格式在整个区域上都是波动的，而FROMM和QUICK格式，虽然有波动却幅值很小。之所以出现这些波动，是因为给区域的出口施加了边界条件。由于高Pe的现象是对流主导的，所以其主要受上游值的影响。在区域的出口，解需要满足一个强制的边界值，就会产生较大的梯度导致发生过冲和下冲。扩散性较高的迎风和SOU格式是仅仅基于上游节点构造的，所以较为光滑且对于所有Pe值并不会表现出过冲和下冲现象，因为它们的解是和施加的出口值无关的。然而SOU格式在较高的梯度区域比如像激波这种情况存在时，依旧会产生振荡（过冲和下冲）。

5.16 均匀和非均匀网格的函数关系

前面展示的不同插值廓线是在一维笛卡尔网格上推导的，如果是沿着曲线坐标轴的话，可以采用同样的函数关系，只需要把笛卡尔坐标轴$x$替换为曲线轴$\zeta$就好了，如上图所示。在均匀网格上，函数关系是完全一样的，不管是笛卡尔网格还是曲线网格系统。在非均匀网格上，函数关系中的独立变量$x$应该用更加通有的$\zeta$来替代，$\zeta$代表沿着坐标轴的距离。如果$O$为坐标原点，那么${\zeta }_U$可用下式计算
$$\begin{array}{rl}{\zeta }_U& ={\zeta }_{UU}+({\zeta }_U-{\zeta }_{UU})\\ {\zeta }_{UU}& =\sqrt{(x_{UU}-x_O{)}^2+(y_{UU}-y_O{)}^2+(z_{UU}-z_O{)}^2}\\ {\zeta }_U-{\zeta }_{UU}& =\sqrt{(x_U-x_{UU}{)}^2+(y_U-y_{UU}{)}^2+(z_U-z_{UU}{)}^2}\end{array}$$
因此，距离的计算可以通过把曲线划分成很多小段，然后每段用直线来代替，通常可用下式计算
$${\zeta }_1-{\zeta }_2=\sqrt{(x_1-x_2{)}^2+(y_1-y_2{)}^2+(z_1-z_2{)}^2}$$
在更加通有的坐标系统上，针对均匀和非均匀网格的不同格式汇总在下表中。

格式	均匀网格	非均匀网格
迎风Upwind	${\varphi }_f={\varphi }_C$	${\varphi }_f={\varphi }_C$
背风Downwind	${\varphi }_f={\varphi }_D$	${\varphi }_f={\varphi }_D$
中心差分CD	${\varphi }_f=\frac{1}{2}({\varphi }_C+{\varphi }_D)$	${\varphi }_f={\varphi }_C+\frac{{\varphi }_D-{\varphi }_C}{{\zeta }_D-{\zeta }_C}({\zeta }_f-{\zeta }_C)$
二阶迎风SOU	${\varphi }_f=\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_U$	${\varphi }_f={\varphi }_C+\frac{{\varphi }_C-{\varphi }_U}{{\zeta }_C-{\zeta }_U}({\zeta }_f-{\zeta }_C)$
QUICK	${\varphi }_f=\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_U+\frac{3}{8}{\varphi }_D$	$\begin{array}{rl}{\varphi }_f={\varphi }_U& +\frac{({\zeta }_f-{\zeta }_U)({\zeta }_f-{\zeta }_C)}{({\zeta }_D-{\zeta }_U)({\zeta }_D-{\zeta }_C)}({\varphi }_D-{\varphi }_U)\\ & +\frac{({\zeta }_f-{\zeta }_U)({\zeta }_f-{\zeta }_D)}{({\zeta }_C-{\zeta }_U)({\zeta }_C-{\zeta }_D)}({\varphi }_C-{\varphi }_U)\end{array}$
FROMM	${\varphi }_f={\varphi }_C+\frac{1}{4}({\varphi }_D-{\varphi }_U)$	${\varphi }_f={\varphi }_C+\frac{{\zeta }_f-{\zeta }_C}{{\zeta }_D-{\zeta }_U}({\varphi }_D-{\varphi }_U)$

6 稳态二维对流

稳态二维对流方程为
$$\nabla \cdot (\rho \mathbf{v}\varphi )=0$$
对如图所示的二维单元$C$上针对体积$V_C$做体积分，使用散度定理，将体积分转化为面积分，并将面积分用单元面的形式表示，得

$$\sum _{f\sim nb(C)}\left({\int }_f{\mathbf{J}}^{\varphi ,C}\cdot d\mathbf{S}\right)=0$$
对于面积分使用单点Gauss积分，上式RHS变为
$$\sum _{f\sim nb(C)}\left({\int }_f{\mathbf{J}}^{\varphi ,C}\cdot d\mathbf{S}\right)=\sum _{f\sim nb(C)}\left({\mathbf{J}}_f^{\varphi ,C}\cdot {\mathbf{S}}_f\right)=\sum _{f\sim nb(C)}(\rho \mathbf{v}\varphi {)}_f\cdot {\mathbf{S}}_f$$
即
$$\sum _{f\sim nb(C)}(\rho \mathbf{v}\varphi {)}_f\cdot {\mathbf{S}}_f=0$$
其在笛卡尔网格上的完全离散形式为
$$\underset{{\dot{m}}_e}{\underbrace{(\rho u\Delta y\varphi {)}_e}}-\underset{-{\dot{m}}_w}{\underbrace{(\rho u\Delta y\varphi {)}_w}}+\underset{{\dot{m}}_n}{\underbrace{(\rho u\Delta x\varphi {)}_n}}-\underset{-{\dot{m}}_s}{\underbrace{(\rho u\Delta x\varphi {)}_s}}=0$$
对于每个坐标方向把流动当成局部一维的情况来处理，采用迎风格式离散，上式变为
$$a_C{\varphi }_C+a_E{\varphi }_E+a_W{\varphi }_W+a_N{\varphi }_N+a_S{\varphi }_S=0$$
系数为
$$\begin{array}{rl}{a}_E& =Flux{F}_e=-||-{\dot{m}}_e,0||\\ a_W& =Flux{F}_w=-||-{\dot{m}}_w,0||\\ a_N& =Flux{F}_n=-||-{\dot{m}}_n,0||\\ a_S& =Flux{F}_s=-||-{\dot{m}}_s,0||\\ a_C& =\sum _{f\sim nb(C)}Flux{C}_f\\ & =||{\dot{m}}_e,0||+||{\dot{m}}_w,0||+||{\dot{m}}_n,0||+||{\dot{m}}_s,0||\\ & =-\underset{{\sum }_{F\sim NB(C)}{a}_F}{\underbrace{(a_E+a_W+a_N+a_S)}}+{\dot{m}}_e+{\dot{m}}_w+{\dot{m}}_n+{\dot{m}}_s\end{array}$$
如果使用QUICK格式，离散方程变为
$$\begin{array}{rl}{a}_C{\varphi }_C& +a_E{\varphi }_E+a_W{\varphi }_W+a_{EE}{\varphi }_{EE}+a_{WW}{\varphi }_{WW}\\ & +a_N{\varphi }_N+a_S{\varphi }_S+a_{NN}{\varphi }_{NN}+a_{SS}{\varphi }_{SS}=0\end{array}$$
其系数为
$$\begin{array}{rl}{a}_E& =Flux{F}_e=-\frac{3}{4}||-{\dot{m}}_e,0||+\frac{3}{8}||{\dot{m}}_e,0||-\frac{1}{8}||{\dot{m}}_w,0||\\ a_W& =Flux{F}_w=-\frac{3}{4}||-{\dot{m}}_w,0||+\frac{3}{8}||{\dot{m}}_w,0||-\frac{1}{8}||{\dot{m}}_e,0||\\ a_{EE}& =Flux{F}_{ee}=\frac{1}{8}||-{\dot{m}}_e,0||\\ a_{WW}& =Flux{F}_{ww}=\frac{1}{8}||-{\dot{m}}_w,0||\\ a_N& =Flux{F}_n=-\frac{3}{4}||-{\dot{m}}_n,0||+\frac{3}{8}||{\dot{m}}_n,0||-\frac{1}{8}||{\dot{m}}_s,0||\\ a_S& =Flux{F}_s=-\frac{3}{4}||-{\dot{m}}_s,0||+\frac{3}{8}||{\dot{m}}_s,0||-\frac{1}{8}||{\dot{m}}_n,0||\\ a_{NN}& =Flux{F}_n=\frac{1}{8}||-{\dot{m}}_n,0||\\ a_{SS}& =Flux{F}_s=\frac{1}{8}||-{\dot{m}}_s,0||\\ a_C& =\sum _{f\sim nb(C)}Flux{C}_f=-\sum _{F\sim NB(C)}{a}_F+{\dot{m}}_e+{\dot{m}}_w+{\dot{m}}_n+{\dot{m}}_s\end{array}$$
采用这两个离散方程来求解在倾斜速度场中含台阶廓线的纯对流问题，如下图所示，矩形计算域内，在速度场$\mathbf{v}(1,1)$的驱动下，变量$\varphi$的纯对流问题，计算域左侧$\varphi$为1，底部$\varphi$为0。在毫无扩散的情况下，如图（a）所示，精确解是对角线上为$\varphi =1$，对角线下为$\varphi =0$。图（b）给出了$x=0.5$处的精确廓线，以及使用迎风格式和QUICK格式得到的数值解。

与精确解相比，迎风格式所得廓线较为模糊且非常不准确但是非常光滑。该不精确性是由一种新形式的误差——交叉扩散（cross-stream diffusion）误差引起的，其是由使用的一维插值廓线导致的，即，将流动视为局部一维的这一假设引起的。起初交叉扩散被Patankar和Stubley定义为多维现象，当流动方向与网格不平齐的时候就产生了这一现象。针对二维问题，Vahl Davis和Mallinson给出了交叉扩散的近似表达式
$${\Gamma }_{false}^{\varphi }=\frac{\rho |\mathbf{v}|\Delta x\Delta y\sin (2\theta )}{4[\Delta y{\sin }^3(\theta )+\Delta x{\cos }^3(\theta )]}$$
其中$|\mathbf{v}|$是速度幅值，而$\theta$是速度矢量与$x$坐标轴的夹角。通过使用高阶格式，可减小该误差，比如QUICK格式。QUICK格式比迎风格式给出的廓线更加陡峭也更加准确，然而在尖锐梯度附近其表现出了过冲和下冲。前面提到过，这种误差叫做色散误差（dispersion error），其生成了求解域内的最大/最小值，正是高阶格式的特征。

6.1 误差来源

基于前面章节的讨论，对流通量离散过程的数值误差来源可分为数值扩散（numerical diffusion）和数值色散（numerical dispersion）。

数值扩散会把尖锐的梯度抹平，其可分为流向和交叉流向的数值扩散。可以通过提高插值廓线的阶数来减小流向数值扩散，这会使廓线更加陡峭，但是会在梯度较大处产生上冲和下冲。交叉流向的数值扩散，是由假设廓线为一维的特性引起的，可以通过在流动方向上做插值或者使用一维高阶插值廓线来减少交叉数值扩散。

数值色散误差表现为在梯度较大处的振荡，使得解变得无界。除了迎风格式，数值色散在所有插值格式中都存在，其是由所假设的插值廓线的非物理特性决定的（这句话的意思是说，在梯度较大的地方，这些格式中假设的变量分布规律跟实际情况是不符合的）。

下图（a）展示的为数值扩散，下图（b）展示的为数值色散。

为评估该误差，针对非定常一维对流扩散方程，忽略其扩散项和源项，并假设密度和速度都是定值，则
$$\frac{\partial \varphi }{\partial t}+u\frac{\partial \varphi }{\partial x}=0$$
假设精确解的形式为
$$\varphi (x,t)=\varphi (t)e^{jkx}$$
其中$j$为虚数（$j^2=-1$），由该精确解算得的梯度为
$$\frac{\partial }{\partial x}\varphi (x,t)=jk\varphi (t)e^{jkx}=jk\varphi (x,t)$$
当使用插值廓线的时候，梯度的数值近似是用在位置$-M\Delta x,(-M+1)\Delta x,...,\Delta x,2\Delta x,...,N\Delta x$处的$\varphi$值来获取的，即
$$\begin{array}{rl}\frac{\partial }{\partial x}\varphi (x,t)& \approx \frac{1}{\Delta x}\sum _{n=-M}^N{a}_n\varphi (x+n\Delta x,t)\\ & =\frac{1}{\Delta x}\sum _{n=-M}^N{a}_n\varphi (t)e^{jk(x+n\Delta x)}\\ & =\frac{1}{\Delta x}\sum _{n=-M}^N{a}_n\varphi (t)e^{jkx}{e}^{jkn\Delta x}\\ & =\frac{1}{\Delta x}\sum _{n=-M}^N{a}_n{e}^{jkn\Delta x}\varphi (x,t)\end{array}$$
与精确解对比，得
$$jk\varphi (x,t)\approx \frac{1}{\Delta x}\sum _{n=-M}^N{a}_n{e}^{jkn\Delta x}\varphi (x,t)\Rightarrow k\approx -\frac{j}{\Delta x}\sum _{n=-M}^N{a}_n{e}^{jkn\Delta x}$$
一般$k$是个虚数，将其写成实部和虚部的形式
$$k=\text{Re}(k)+j\text{Im}(k)$$
其中的Re和Im分别代表实部和虚部，将上式代入精确解$\varphi (x,t)=\varphi (t)e^{jkx}$，得
$$\begin{array}{rl}\varphi (x,t)& =\varphi (t)e^{jkx}\\ & \approx \varphi (t)e^{j[\text{Re}(k)+j\text{Im}(k)]x}\\ & =\varphi (t)\underset{\begin{array}{rl}& Phase\\ & Dispersion\end{array}}{\underbrace{e^{j\text{Re}(k)x}}}\underset{\begin{array}{rl}& Amplitude\\ & Dissipation\end{array}}{\underbrace{e^{-\text{Im}(k)x}}}\end{array}$$
因此，数值解包含了扩散diffusion（或耗散dissipative）和色散dispersion误差。如果$k$为实数，则只有色散误差。如果$k$是复数，则两种误差都会出现。基于此分析，检查下迎风格式和中心差分格式的$k$值，对于迎风格式，其梯度算出来是
$$\begin{array}{rl}\frac{\partial \varphi }{\partial x}& \approx \frac{{\varphi }_e-{\varphi }_w}{\Delta x}=\frac{{\varphi }_C-{\varphi }_W}{\Delta x}=\frac{e^{jkx}-e^{jk(x-\delta x)}}{\Delta x}\varphi (x,t)\\ & =\frac{1-\cos (k\delta x)+j\sin (k\delta x)}{\Delta x}\varphi (x,t)\\ & =jk\varphi (x,t)\\ \Rightarrow k& =\frac{\sin (k\delta x)}{\Delta x}-j\frac{1-\cos (k\delta x)}{\Delta x}\end{array}$$
显然，迎风格式既有扩散误差又有色散误差。对于中心差分格式，其梯度项为
$$\begin{array}{rl}\frac{\partial \varphi }{\partial x}& \approx \frac{{\varphi }_e-{\varphi }_w}{\Delta x}=\frac{{\varphi }_E-{\varphi }_W}{\Delta x}=\frac{e^{jk(x+\delta x)}-e^{jk(x-\delta x)}}{\Delta x}\varphi (x,t)\\ & =\frac{j\sin (k\delta x)}{\Delta x}\varphi (x,t)\\ & =jk\varphi (x,t)\\ \Rightarrow k& =\frac{\sin (k\delta x)}{\Delta x}\end{array}$$
因为$k$是实数，所以只有色散误差，色散误差导致解中出现震荡、上冲和下冲。

对数值色散有了更加深入的理解后，就需要针对对流项发展高阶精度的无震荡格式，这可一直是研究的热点，直到人们理解了对流通量的界定bounding。下一章将会详细讲解这类格式的发展。

7 非结构网格上的高阶格式

在结构网格上，高阶精度的函数关系是用远迎风$U$、迎风$C$、背风$D$节点的值来定义的，然而在非结构网格上，界面的$C$和$D$节点是确定的，而$U$节点则无法知晓在哪，如图。针对该难题，一种简单的处理方式是把高阶格式重新定义成节点$C$和$D$处梯度的形式，另一种方法是重构一个伪节点$U$，这个在下一章中再讲，其主要用于发展高分辨率格式。

7.1 将高阶格式写成梯度项的形式

前面讲过的QUICK格式为例，其函数关系为
$$\begin{array}{rl}{\varphi }_f& =\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_U+\frac{3}{8}{\varphi }_D\\ & ={\varphi }_C+\frac{1}{4}\left(\frac{{\varphi }_D-{\varphi }_U}{2}\right)+\frac{1}{4}({\varphi }_D-{\varphi }_C)\end{array}$$

沿着${\mathbf{d}}_{CF}$或$\zeta$方向，位于节点$C$和面$f$处的梯度为
$$\frac{\partial {\varphi }_C}{\partial \zeta }=\frac{{\varphi }_D-{\varphi }_U}{2\Delta \zeta }\frac{\partial {\varphi }_f}{\partial \zeta }=\frac{{\varphi }_D-{\varphi }_C}{\Delta \zeta }$$
跟QUICK格式对比，发现
$$\begin{array}{rl}{\varphi }_f& ={\varphi }_C+\frac{1}{4}\left(\frac{{\varphi }_D-{\varphi }_U}{2}\right)+\frac{1}{4}({\varphi }_D-{\varphi }_C)\\ & ={\varphi }_C+\frac{1}{2}\left(\frac{{\varphi }_D-{\varphi }_U}{2\Delta \zeta }\right)\frac{\Delta \zeta }{2}+\frac{1}{2}\left(\frac{{\varphi }_D-{\varphi }_C}{\Delta \zeta }\right)\frac{\Delta \zeta }{2}\end{array}$$
即
$${\varphi }_f={\varphi }_C+\frac{1}{2}\frac{\partial {\varphi }_C}{\partial \zeta }\frac{\Delta \zeta }{2}+\frac{1}{2}\frac{\partial {\varphi }_f}{\partial \zeta }\frac{\Delta \zeta }{2}$$
把$C$和$f$之间的矢量定义成${\mathbf{d}}_{Cf}$，上式变为
$${\varphi }_f={\varphi }_C+\frac{1}{2}\nabla {\varphi }_C\cdot {\mathbf{d}}_{Cf}+\frac{1}{2}\nabla {\varphi }_f\cdot {\mathbf{d}}_{Cf}$$
在非结构网格上使用这种形式是非常适合的，因为它只需要提供$C$和$f$处的梯度，只要对于这些梯度的计算是二阶精度的就好了，梯度的计算前面章节已经详细讲过了。与其最初的格式相比，这种格式也更加灵活，尤其适合在非结构网格上的计算。同样的，上式表明，基于3点的格式可以写成下面这种梯度的组合形式
$${\varphi }_f=a{\varphi }_C+b\nabla {\varphi }_C\cdot {\mathbf{d}}_{Cf}+c\nabla {\varphi }_f\cdot {\mathbf{d}}_{Cf}$$
其中的$a,b,c$是系数，由在非结构网格上确定${\varphi }_f$时所采用的假设廓线分布所决定。这个通有离散格式可以用于其他格式的推导，先把它写成梯度值近似表达的形式，即
$$\begin{array}{rl}{\varphi }_f& =a{\varphi }_C+b\nabla {\varphi }_C\cdot {\mathbf{d}}_{Cf}+c\nabla {\varphi }_f\cdot {\mathbf{d}}_{Cf}\\ & =a{\varphi }_C+b\frac{{\varphi }_D-{\varphi }_U}{2\Delta \zeta }\frac{\Delta \zeta }{2}+c\frac{{\varphi }_D-{\varphi }_C}{\Delta \zeta }\frac{\Delta \zeta }{2}\end{array}$$
然后可得最终形式为
$${\varphi }_f=\left(a-\frac{c}{2}\right){\varphi }_C+\left(\frac{b}{4}+\frac{c}{2}\right){\varphi }_D-\frac{b}{4}{\varphi }_U$$
使用上式，很容易就能算出$a,b,c$系数来。以二阶迎风格式SOD为例
$$\begin{array}{rl}& {\varphi }_f=\left(a-\frac{c}{2}\right){\varphi }_C+\left(\frac{b}{4}+\frac{c}{2}\right){\varphi }_D-\frac{b}{4}{\varphi }_U\\ & {\varphi }_f=\frac{3}{2}{\varphi }_C-\frac{1}{2}{\varphi }_U\end{array}\}\Rightarrow \{\begin{array}{rl}-\frac{b}{4}=-\frac{1}{2}& \Rightarrow b=2\\ \frac{b}{4}+\frac{c}{2}=0& \Rightarrow c=-1\\ a-\frac{c}{2}=\frac{3}{2}& \Rightarrow a=1\end{array}$$
推得二阶迎风格式的等效形式
$${\varphi }_f={\varphi }_C+\left(2\nabla {\varphi }_C-\nabla {\varphi }_f\right)\cdot {\mathbf{d}}_{Cf}$$
按照同样的方法，可以推得前面所讲格式的等效形式，汇总如下

格式	非结构网格等效形式
迎风Upwind	${\varphi }_f={\varphi }_C$
中心差分CD	${\varphi }_f={\varphi }_C+\nabla {\varphi }_f\cdot {\mathbf{d}}_{Cf}$
二阶迎风SOU	${\varphi }_f={\varphi }_C+\left(2\nabla {\varphi }_C-\nabla {\varphi }_f\right)\cdot {\mathbf{d}}_{Cf}$
FROMM	${\varphi }_f={\varphi }_C+\nabla {\varphi }_C\cdot {\mathbf{d}}_{Cf}$
QUICK	${\varphi }_f={\varphi }_C+\frac{1}{2}\left(\nabla {\varphi }_C+\nabla {\varphi }_f\right)\cdot {\mathbf{d}}_{Cf}$
背风Downwind	${\varphi }_f={\varphi }_C+2\nabla {\varphi }_f\cdot {\mathbf{d}}_{Cf}$

关于$\nabla {\varphi }_C$和$\nabla {\varphi }_f$的算法，前面第9章已经详细讲过了，如果忘掉的话，去看下吧。

8 延迟修正方法（The Deferred Correction Approach）

Khosla和Rubin提出的延迟修正（DC）过程，是一种紧凑的技术，其允许最初是针对低阶算法编写的代码，可以来使用高阶格式，而不会影响稳定性。该方法既可适用于结构网格也可适用于非结构网格，其思想是把使用高阶格式HO的单元面$f$的对流通量写成
$${\dot{m}}_f{\varphi }_f^{HO}=\underset{implicit}{\underbrace{{\dot{m}}_f{\varphi }_f^U}}+\underset{explicit}{\underbrace{{\dot{m}}_f({\varphi }_f^{HO}-{\varphi }_f^U)}}$$
其中上标$U$和$HO$分别表示迎风和高阶格式，通过把对流项写成这种格式，RHS的第1项是隐式计算的，其用的是节点值，而RHS的第2项是显示计算的，其用的是最近的$\varphi$值，即，在迭代过程中上一次的迭代值。上式可展开为节点值的形式
$$\begin{array}{rl}{\dot{m}}_f{\varphi }_f^{HO}& =||{\dot{m}}_f,0||{\varphi }_C-||-{\dot{m}}_f,0||{\varphi }_F+({\dot{m}}_f{\varphi }_f^{HO}-||{\dot{m}}_f,0||{\varphi }_C+||-{\dot{m}}_f,0||{\varphi }_F)\\ & =\underset{implicit}{\underbrace{Flux{C}_f{\varphi }_C+Flux{F}_f{\varphi }_F}}+\underset{explicit}{\underbrace{Flux{V}_f}}\end{array}$$
其中
$$\begin{array}{rl}& Flux{C}_f=||{\dot{m}}_f,0||\\ & Flux{F}_f=-||-{\dot{m}}_f,0||\\ & Flux{V}_f={\dot{m}}_f{\varphi }_f^{HO}-Flux{C}_f{\varphi }_C-Flux{F}_f{\varphi }_F\end{array}$$
可以推得最终的代数方程为
$$a_C{\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F=b_C$$
其中
$$\begin{array}{rl}{a}_F& =Flux{F}_f=-||-{\dot{m}}_f,0||\\ a_C& =\sum _{f\sim nb(C)}Flux{C}_f=\sum _{f\sim nb(C)}||{\dot{m}}_f,0||=\sum _{f\sim nb(C)}({\dot{m}}_f+||-{\dot{m}}_f,0||)\\ & =-\sum _{F\sim NB(C)}{a}_F+\sum _{f\sim nb(C)}{\dot{m}}_f\\ b_C& =Q_C^{\varphi }{V}_C-\underset{b_C^{DC}}{\underbrace{\sum _{f\sim nb(C)}Flux{V}_f}}=Q_C^{\varphi }{V}_C-\underset{b_C^{DC}}{\underbrace{\sum _{f\sim nb(C)}{\dot{m}}_f({\varphi }_f^{HO}-{\varphi }_f^U)}}\end{array}$$
其中的$b_C^{DC}$代表了延迟修正（DC）过程所引入的额外源项，此外，DC技术生成的方程中系数矩阵总是对角占优的，因为其使用了迎风格式。

---

例2 采用QUICK格式，推导延迟修正的系数。

解
守恒方程将离散成如下的代数方程形式
$$a_C{\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F=b_C$$
QUICK给出的单元面的值为（下式中的C、U、D表示迎风、远迎风、背风节点）
$${\varphi }_f=\frac{3}{4}{\varphi }_C-\frac{1}{8}{\varphi }_U+\frac{3}{8}{\varphi }_D$$
延迟修正方法是基于迎风格式得到代数方程的系数的，这些系数为
$$\begin{array}{rl}{a}_F& =Flux{F}_f=-||-{\dot{m}}_f,0||\\ a_C& =-\sum _{F\sim NB(C)}{a}_F+\sum _{f\sim nb(C)}{\dot{m}}_f\end{array}$$
迎风和QUICK格式的差别在于其源项中的显示计算过程不同，即
$$\begin{array}{rl}{b}_C& =S_C^{\varphi }{V}_C-\sum _{f\sim nb(C)}{\dot{m}}_f\left({\varphi }_f^{HO}-{\varphi }_f^{Upwind}\right)\\ & =S_C^{\varphi }{V}_C-\sum _{f\sim nb(C)}{\dot{m}}_f\left(-\frac{1}{4}{\varphi }_C-\frac{1}{8}{\varphi }_U+\frac{3}{8}{\varphi }_D\right)\end{array}$$

---

该紧凑流程非常容易实现，然而随着迎风格式和高阶格式算出来的单元面值的差别的增大，其收敛速度也变慢了。

9 代码讲解

9.1 uFVM

在uFVM中，对流项的组装是在cfdAssembleConvectionTerm中进行的，采用迎风格式先计算内部面的（cfdAssembleConvectionTermInterior），然后再对不同边界片Patch来计算边界面的（cfdAssembleConvectionTermInletBC，cfdAssembleConvectionTermOutletBC等）。

内部面的组装过程如下，首先找到内部面的owner和neighbour标识，然后使用pos来定义迎风标识（=1表示迎风），最后计算迎风格式的相关系数。

theMdotName = ['Mdot' theFluidTag];
mdotField = cfdGetMeshField(theMdotName,'Faces');	% 质量流量场

mdot_f = mdotField.phi(iFaces);						% 面的流出质量流量

iOwners = [theMesh.faces(iFaces).iOwner];			% 面所属单元标识
iNeighbours = [theMesh.faces(iFaces).iNeighbour];	% 面邻近单元标识
pos = zeros(size(mdot_f));
pos((mdot_f>0))=1;									% 依据mdot_f的正负判断是否为迎风方向节点

theFluxes.FLUXC1f(iFaces,1) = mdot_f.*pos;			% 给迎风格式的相关系数赋值
theFluxes.FLUXC2f(iFaces,1) = mdot_f.*(1-pos);
theFluxes.FLUXVf(iFaces,1) = 0;

对于高阶格式的添加是在迎风格式组装之后使用延迟修正方法来实现的，QUICK格式的组装在cfdAssembleConvectionTermDCQUICK中实现，如下

iUpwind = pos.*iOwners + (1-pos).*iNeighbours;	% 计算迎风侧单元标识
%get the upwind gradient at the interior faces
phiGradCf = phiGrad(iUpwind,:,iComponent); 		% 提取界面迎风侧单元的梯度
%interpolated gradient to interior faces
iOwners = [theMesh.faces(iFaces).iOwner]';
iNeighbours = [theMesh.faces(iFaces).iNeighbour]';
pos = zeros(size(mdot_f));
pos(mdot_f>0)=1;
% 插值计算面上的梯度值
phiGradf = ...
cfdInterpolateGradientsFromElementsToInteriorFaces('Average:Corrected'
,phiGrad,phi);
rc = [theMesh.elements(iUpwind).centroid]';
rf = [theMesh.faces(iFaces).centroid]';
rCf = rf-rc;	% 单元形心到面形心的距离矢量
% 计算延迟修正中的修正项
corr = mdot_f .* dot(phiGradCf'+phiGradf',rCf')'*0.5;
% 把修正项加入到源项中去
theFluxes.FLUXTf(iFaces) = theFluxes.FLUXTf(iFaces) + corr;

上述代码中对内部面的修正项corr的计算，即
$$\begin{array}{rl}{\dot{m}}_f\left({\varphi }_f^{HO}-{\varphi }_f^{Upwind}\right)& ={\dot{m}}_f\left[\underset{QUICK}{\underbrace{{\varphi }_C+\frac{1}{2}\left(\nabla {\varphi }_C+\nabla {\varphi }_f\right)\cdot {\mathbf{d}}_{Cf}}}-\underset{Upwind}{\underbrace{{\varphi }_C}}\right]\\ & ={\dot{m}}_f\frac{1}{2}\left(\nabla {\varphi }_C+\nabla {\varphi }_f\right)\cdot {\mathbf{d}}_{Cf}\end{array}$$

9.2 OpenFOAM

to be continued…