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 \phi)}{d x} - \frac{d}{dx}\left( \Gamma^{\phi}\frac{d \phi}{d x} \right)=0$
萬分幸運的是，這個問題存在着解析解，因此它一般用做基準算例來跟各類各樣的數值格式作對照。

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$的變化也逐漸變爲了更加陡峭的廓線形式。

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$須要用鄰近的節點值來獲得，肯定這些面上值的方法就是文獻中所謂的「對流格式」。

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)}^{}}\end{array}$