学习自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 Termhtml

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

1 引言

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

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

2 稳态一维对流和扩散

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

2.1 解析解

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

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

2.2 数值解

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

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

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

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

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

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

$a_E=-\Gamma^\phi_e \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^\phi} > 2$
定义Peclet数为
$Pe=\frac{\rho u\delta x}{\Gamma^\phi}$
这样，对于均匀网格来讲， $Pe=Pe_L/2$，那么 $a_E>0$的条件变为
$Pe>2$
所以，单元的Peclet数（Pe）大于2的话，离散过程变得不相符，由于此时邻居值的增长会致使 $C$处值的下降，这又会反过来进一步增长邻居值，致使错误的放大。

固然，经过减少网格尺度让单元Peclet数小于2，就能够避免这种状况。然而，对于不少实际状况，这样作的话将会显著增长存储量和计算量，以致超出计算机负荷。并且，对于纯对流的问题（即不考虑粘性扩散效应，例如，Euler流动），这么作的话也是不可行的，因此须要发展修正方法。

2.4 迎风格式（Upwind Scheme）

如上图，与对流过程更加契合的格式为迎风格式，迎风格式基本上是模仿了对流的物理特性，即，单元面上的值是依赖于迎风节点上的值的，换言之，依赖于流动方向的。图中所示的单元面上的值用下式给出
$\phi_e=\begin{cases} \phi_C & if & \dot m_e > 0 \\ \phi_E & if & \dot m_e < 0 \end{cases} \quad\quad\quad \phi_w=\begin{cases} \phi_C & if & \dot m_w > 0 \\ \phi_W & if & \dot m_w < 0 \end{cases}$
（注意，这里的 $\phi_w$没有写反，由于质量流量 $\dot m$跟速度 $u$是不一样的，看下面说明）

其中 $\dot m_e$和 $\dot m_w$是在面 $e$和 $w$的质量流量（面积矢量的正方向是朝外的，而非指向单元中心 $C$的）
$\begin{aligned} & \dot m_e=(\rho \bold v \cdot \bold S)_e=(\rho u S)_e=(\rho u \Delta y)_e \\ & \dot m_w=(\rho \bold v \cdot \bold S)_w=-(\rho u S)_w=-(\rho u \Delta y)_w \end{aligned}$
如此，面 $e$的对流通量可写成
$\begin{aligned} \dot m_e \phi_e &= ||\dot m_e,0||\phi_C-||-\dot m_e,0||\phi_E \\ &= FluxC_e^{Conv}\phi_C+FluxF_e^{Conv}\phi_E+FluxV_e^{Conv} \end{aligned}$
其中
$\begin{aligned} &FluxC_e^{Conv}= ||\dot m_e,0|| \\ &FluxF_e^{Conv}=-||-\dot m_e,0|| \\ &FluxV_e^{Conv}=0 \end{aligned}$
其中的 $||a,b||=\max(a,b)$表明了 $a$和 $b$中的最大值。一样手法，也可推导面 $w$的对流通量
$\begin{aligned} \dot m_w \phi_w &= ||\dot m_w,0||\phi_C-||-\dot m_w,0||\phi_W \\ &= FluxC_w^{Conv}\phi_C+FluxF_w^{Conv}\phi_W+FluxV_w^{Conv} \end{aligned}$
其中
$\begin{aligned} &FluxC_w^{Conv}= ||\dot m_w,0|| \\ &FluxF_w^{Conv}=-||-\dot m_w,0|| \\ &FluxV_w^{Conv}=0 \end{aligned}$
把扩散项用上标 ${Diff}$表示，推得一维稳态对流扩散方程的离散格式为
$\begin{aligned} (FluxC_e^{Conv}&+FluxC_e^{Diff}+FluxC_w^{Conv}+FluxC_w^{Diff})\phi_C\\ &+(FluxF_e^{Conv}+FluxF_e^{Diff})\phi_E+(FluxF_w^{Conv}+FluxF_w^{Diff})\phi_W=0 \end{aligned}$
可写成下面形式
$a_C\phi_C+a_E\phi_E+a_W\phi_W=b_C$
其中
$\begin{aligned} a_E&= FluxF_e^{Conv}+FluxF_e^{Diff} \\ &= -||-\dot m_e,0|| - \Gamma^\phi_e \frac{S_e}{\delta x_e} \\ a_W&= FluxF_w^{Conv}+FluxF_w^{Diff} \\ &=-||-\dot m_w,0|| -\Gamma^\phi_w \frac{S_w}{\delta x_w} \\ a_C&= \sum_f \left( FluxC_f^{Conv} + FluxC_f^{Diff} \right)\\ &=||\dot m_e,0||+||\dot m_w,0||+\Gamma^\phi_e \frac{S_e}{\delta x_e}+\Gamma^\phi_w \frac{S_w}{\delta x_w} \\ &=-(a_E+a_W)+\underbrace{(\dot m_e + \dot m_w)}_{=0} \\ b_C&=-\sum_f \left( FluxV_f^{Conv} + FluxV_f^{Diff} \right)\\ &=0 \end{aligned}$
显而易见，迎风格式推出了负的邻居系数，并且在知足连续方程的状况下（即， $\dot m_e + \dot m_w=0$），主节点上的系数为
$a_C=-(a_W+a_E)$
这确保了有界性。

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

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

2.5 背风格式（Downwind Scheme）

若是把迎风格式反过来用，即，用背风格式的话，会发生什么有趣的事情呢？如上图所示，在这种插值廓线中，是把界面背风方向节点处的值视做界面上的值，这样，面 $e$和面 $w$处的值为
$\phi_e=\begin{cases} \phi_E & if & \dot m_e>0 \\ \phi_C & if & \dot m_e<0 \end{cases} \quad\quad\quad \phi_w=\begin{cases} \phi_W & if & \dot m_w>0 \\ \phi_C & if & \dot m_w<0 \end{cases}$
面上的对流通量能够写成
$\begin{aligned} \dot m_e\phi_e & = - ||-\dot m_e,0||\phi_C+||\dot m_e,0||\phi_E \\ &= FluxC_e^{Conv}\phi_C+FluxF_e^{Conv}\phi_E+FluxV_e^{Conv} \\ \dot m_w\phi_w & = - ||-\dot m_w,0||\phi_C+||\dot m_w,0||\phi_W \\ &= FluxC_w^{Conv}\phi_C+FluxF_w^{Conv}\phi_W+FluxV_w^{Conv} \end{aligned}$
可进一步推得一维稳态对流扩散方程的离散形式为
$\begin{aligned} (FluxC_e^{Conv}&+FluxC_e^{Diff}+FluxC_w^{Conv}+FluxC_w^{Diff})\phi_C\\ &+(FluxF_e^{Conv}+FluxF_e^{Diff})\phi_E+(FluxF_w^{Conv}+FluxF_w^{Diff})\phi_W=0 \end{aligned}$
可写成下面形式
$a_C\phi_C+a_E\phi_E+a_W\phi_W=b_C$
其中
$\begin{aligned} a_E&= FluxF_e^{Conv}+FluxF_e^{Diff} \\ &= ||\dot m_e,0|| - \Gamma^\phi_e \frac{S_e}{\delta x_e} \\ a_W&= FluxF_w^{Conv}+FluxF_w^{Diff} \\ &=||\dot m_w,0|| -\Gamma^\phi_w \frac{S_w}{\delta x_w} \\ a_C&= \sum_f \left( FluxC_f^{Conv} + FluxC_f^{Diff} \right)\\ &=-||-\dot m_e,0||-||-\dot m_w,0||+\Gamma^\phi_e \frac{S_e}{\delta x_e}+\Gamma^\phi_w \frac{S_w}{\delta x_w} \\ &=-(a_E+a_W)+\underbrace{(\dot m_e + \dot m_w)}_{=0} \\ b_C&=-\sum_f \left( FluxV_f^{Conv} + FluxV_f^{Diff} \right)\\ &=0 \end{aligned}$
在连续方程知足的前提下，假设为均匀网格和常扩散系数，用 $\phi_E$和 $\phi_W$所给出的 $\phi_C$值为
$\frac{\phi_C-\phi_W}{\phi_E-\phi_W} =\frac{2-||Pe_L,0||}{4-||-Pe_L,0||-||Pe_L,0||} =\frac{2-||Pe_L,0||}{4-|Pe_L|}$
不用把上式的图形画出来，也能够看出，当 $|Pe_L|\rightarrow4$的时候，解变得彻底无界。背风格式与其它格式联合起来预测突变界面时是有用的，然而这里引入背风格式的目的主要是为了便于在下一节给出更好的稳定性分析。

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

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

3.1 迎风格式

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

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

3.2 背风格式

跟上一小节相似，假设流动是沿着 $x$正方向的，将 $\phi_E$和 $\phi_C$分别写成 $\phi_e$和 $\phi_w$函数的形式，背风格式为
$\phi_e=\phi_E\quad\quad\quad\phi_w=\phi_C$
这样，使用背风格式的一维对流扩散方程离散形式为
$(\rho u \Delta y)_e\phi_E-(\rho u \Delta y)_w\phi_C- \left[\left(\Gamma^\phi \frac{d\phi}{dx}\Delta y\right)_e-\left(\Gamma^\phi \frac{d\phi}{dx}\Delta y\right)_w\right]=0$
在均匀网格上，将 $\phi_E$和 $\phi_C$作一维Taylor展开，分别以单元面 $e$和 $w$的值为基准，得
$\begin{array}{rl}{\displaystyle }& {\displaystyle {\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+\mathrm{.}\mathrm{.}\mathrm{.}}\end{array}$