FVM in CFD 學習筆記_第13章_時域離散：瞬態項

學習自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 13 Temporal Discretization: The Transient Term

前面章節的討論中均假設是穩態條件，因此無需對瞬態項進行離散。若是考慮瞬態現象的話，則至關於給問題增長了一個新的維度。然而，因爲瞬態項的性質是拋物型的，故而不須要再額外定義時間域上的場，即，不像在空間域上須要用 $\phi_i,i=1,2,...,N$那樣子來離散節點值。通常來講，瞬態項的離散只須要額外存儲一到兩個時間層上的變量場便可，由所選擇格式的數值精度而定。與穩態問題的另外一個不一樣點是，瞬態系統是用時間步長推動方式來模擬的，即，以 $t=t_0$時刻的初始條件開始，求解算法向前推動並找到 $t_1=t_0+\Delta t_1$時刻的解，這個找到的解再做爲初始條件去找尋 $t_2=t_1+\Delta t_2$時刻的解，該過程不斷重複直到達到所需的時刻爲止。本章的重點是瞬態項離散的技術，將展現兩種發展瞬態格式的方法。其一是使用Taylor展開來把瞬態項展開成節點值的形式，這在有限差分方法中很是奏效；其二是有限體積方法中經常使用的僞時間單元方法，和在對流項中的僞節點很是相似。將展現一些瞬態格式，並討論它們的特性。

1 引言

對於瞬態模擬而言，控制方程的離散是在空間域和時間域上進行的。空間域上的空間離散和穩態問題同樣，時間離散則要設置一個時間座標軸，沿着該座標軸來計算瞬態項的導數（有限差分方法）或積分（有限體積方法）值。

通常而言，某個變量 $\phi$的瞬態特性的表達式，或其時間演化，是由下述形式的方程控制的
$\frac{\partial (\rho \phi)}{\partial t}+L(\phi)=0$
其中 $L(\phi)$爲空間算子，其包含了全部的非瞬態項（擴散項、對流項、源項，等）， ${\partial (\rho \phi)}/{\partial t}$爲瞬態算子，二者以下圖所示。

在單元 $C$上對上式作積分，得
$\int_{V_C}\frac{\partial (\rho \phi)}{\partial t}dV+\int_{V_C}L(\phi)dV=0$
空間離散後，變爲
$\frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0$
其中 $V_C$爲離散單元的體積， $L(\phi_C^t)$爲空間算子在某參考時刻 $t$的離散形式，可寫成以下代數方程
$L(\phi_C^t)=a_C\phi_C^t+\sum_{F\sim NB(C)}a_F\phi_F^t-b_C$
在上上式中，當 $t\rightarrow\infin$時其變爲穩態離散方程。經過時間推動也能達到穩態，即， $\phi_C^{t+\Delta t}=\phi_C^t$。這保證了瞬態問題在到達穩態時所得到的解，和直接由穩態問題求得的解，是同樣的。

在瞬態項的離散中，最經典的法子是用有限差分法的思想，把 ${\partial (\rho \phi)}/{\partial t}$作Taylor級數展開，獲得用離散點的值所表示的導數項形式。在本章中，將展現另外一種與有限體積方法更爲契合的方法，對 ${\partial (\rho \phi)}/{\partial t}$在時間單元上作積分，並轉化爲面通量的形式，這和對流格式的處理方法很像，只是如今離散是沿着時間軸進行的。

2 有限差分方法

如圖，在瞬態空間中的網格是結構化的，因此對瞬態項使用有限差分方法是很是稀鬆日常的操做。在該方法中，空間算子 $L(\phi)$是在時刻 $t$離散的，而瞬態偏導數則可用時刻 $t$的Taylor展開而推導出不一樣的瞬態格式，其中的一些格式呈現以下。

2.1 前向Euler格式（Forward Euler Scheme）

爲了估算瞬態項，所推導量的Taylor展開須要指定時間方向。這裏，展開是針對時間向前進行的。對某函數 $T$，其在 $t+\Delta t$時刻的值，可用其在時刻 $t$的值和導數值作Taylor級數展開，以下
$T(t+\Delta t)=T(t)+\frac{\partial T(t)}{\partial t}\Delta t+\frac{\partial^2 T(t)}{\partial t^2}\frac{\Delta t^2}{2!}+...$
截去二階以上項，獲得
$\frac{\partial T(t)}{\partial t}=\frac{T(t+\Delta t)-T(t)}{\Delta t}+O(\Delta t)$
只有一階精度，把上式中的 $T$替換成 $(\rho\phi)$，並代入到 $\frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0$，離散方程變爲
$\frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^t}{\Delta t}V_C+L(\phi_C^t)=0$

該瞬態格式如上圖所示，計算 $t+\Delta t$時刻的 $(\rho_C \phi_C)$是不須要求解方程組系統的。由於全部的空間項都是在舊時刻 $t$上衡量的，因此 $t+\Delta t$時刻的 $\phi_C$可直接由上個時間步的值顯式算出。這樣的格式屬於顯式瞬態格式。顯式瞬態格式的主要特色是經過在時域推動即可得到解，而無需在每一個時間層求解方程組系統。這樣一來，求解效率很是高，且對計算網格的並行處理十分方便。然而鮮有商業軟件採用該方法，由於其重大缺陷就是時間步長 $\Delta t$極爲受限，這將在下一小節詳細展開。

把空間算子的離散代數方程代入到上式，可得完整的代數方程爲
$a_C^{t+\Delta t} \phi_C^{t+\Delta t} + a_C^t \phi_C^t = b_C-\left(a_C\phi_C^t+\sum_{F\sim NB(C)}a_F\phi_F^t\right)$
其中
$\begin{aligned} a_C^{t+\Delta t} &= \frac{\rho_C ^{t+\Delta t}V_C}{\Delta t} \\ a_C^t &= -\frac{\rho_C^tV_C}{\Delta t} \end{aligned}$
在上述方程中， $a_C^{t+\Delta t}$和 $a_C^t$爲對角係數，源自於瞬態項的離散， $\phi_C^{t+\Delta t}$和 $\phi_C^t$爲時間層 $t+\Delta t$和 $t$上的值，而 $a_C$、 $a_F$、 $b_C$則是空間離散後的係數。

爲了簡化標記，在本章中，前一個時間步的變量值用上標 $^\circ$標記，前兩個時間步的變量值用上標 $^{\circ\circ}$標記，如果沒有上標，則表示當前時間步的變量值，除了在非定常項中乘到 $\phi_C$上的係數是用上標 $^\bullet$來標記的。基於這些標記，上式寫做
$a_C^{\bullet} \phi_C + a_C^\circ \phi_C^\circ = b_C-\left(a_C\phi_C^\circ+\sum_{F\sim NB(C)}a_F\phi_F^\circ\right)$
其中
$\begin{aligned} a_C^{\bullet} &= \frac{\rho_C V_C}{\Delta t} \\ a_C^\circ &= -\frac{\rho_C^\circ V_C}{\Delta t} \end{aligned}$
改寫爲以下形式
$\phi_C = \frac{b_C-\left[(a_C+a_C^\circ)\phi_C^\circ+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F^\circ\right]}{a_C^{\bullet}}$
顯然，當前時間步的 $\phi$值是由顯式關係算出來的，不須要求解方程組系統。

2.2 前向Euler格式的穩定性

數值格式的收斂性和穩定性最初是由Courant、Friedrichs、Lewy所探究的，他們代表，爲了讓差分方程的解收斂到原偏微分方程的解，數值格式中必須用到影響解的初始數據所包含的全部信息。該要求隨後就變成了著名的CFL條件。

實際上，CFL條件能夠被簡單解釋爲，把係數所要知足的基本規則之一，即，反符號準則，擴展成把瞬態係數也包含在內。這樣，正如 $\phi_F$是 $\phi_C$的「空間」鄰居， $\phi_C^\circ$也是 $\phi_C$的「時間」鄰居，那麼反符號準則對他倆應該平等適用。注意對角線係數如今是 $a_C^\bullet$，其「時間」鄰居的係數是 $(a_C+a_C^\circ)$，反符號準則變爲
$a_C+a_C^\circ\le0$

2.2.1 瞬態-對流狀況的穩定性

如上圖，對於一維純對流問題，流動方向向右，使用迎風格式把變量插值到單元面上去，則單元 $C$的離散方程中的係數 $a_C$和 $a_C^\circ$爲
$\begin{aligned} a_C&=\dot m_e^\circ=\rho_C^\circ u_C^\circ \Delta y_C \\ a_C^\circ&=-\frac{\rho_C^\circ V_C}{\Delta t}=-\frac{\rho_C^\circ \Delta x_C \Delta y_C}{\Delta t} \end{aligned}$
所以，CFL條件爲
$a_C+a_C^\circ\le0 \Rightarrow \rho_C^\circ u_C^\circ \Delta y_C-\frac{\rho_C^\circ \Delta x_C \Delta y_C}{\Delta t}\le0$
即
$\Delta t\le\frac{\Delta x_C}{u_C^\circ}$
對於對流控制的流動，定義CFL數爲
$CFL^{conv}=\frac{|\bold v_C^\circ|\Delta t}{\Delta x_C}$
這意味着爲了數值穩定性，CFL數應該知足
$CFL^{conv}\le1$

2.2.2 瞬態-擴散狀況的穩定性

對於純擴散問題，CFL數的表達式是不一樣的。所以，分析如上圖所示的一維純擴散問題。

使用線性插值廓線，單元 $C$的離散方程中的係數 $a_C$和 $a_C^\circ$爲
$\begin{aligned} a_C&=\frac{\Gamma_e^\phi \Delta y_C}{\delta x_e} + \frac{\Gamma_w^\phi \Delta y_C}{\delta x_w} \\ a_C^\circ&=-\frac{\rho_C^\circ V_C}{\Delta t}=-\frac{\rho_C^\circ \Delta x_C \Delta y_C}{\Delta t} \end{aligned}$
所以，CFL條件須要
$a_C+a_C^\circ\le0 \Rightarrow \frac{\Gamma_e^\phi \Delta y_C}{\delta x_e} + \frac{\Gamma_w^\phi \Delta y_C}{\delta x_w}-\frac{\rho_C^\circ \Delta x_C \Delta y_C}{\Delta t}\le0$
即
$\Delta t\le\frac{\rho_C^\circ\Delta x_C}{\dfrac{\Gamma_e^\phi}{\delta x_e} + \dfrac{\Gamma_w^\phi}{\delta x_w}}$
若網格均勻，且擴散係數爲常數，則上式變爲
$\Delta t\le\frac{\rho_C^\circ(\Delta x_C)^2}{2\Gamma_C^\phi}$
對於擴散控制的問題，定義CFL數爲
$CFL^{diff}=\frac{\Gamma_C^\phi\Delta t}{\rho_C^\circ(\Delta x_C)^2}$
這意味着爲了數值穩定性，CFL數應該知足
$CFL^{diff}\le \frac{1}{2}$

2.2.3 瞬態-對流-擴散狀況的穩定性

對於多維非定常對流擴散問題，基於上一章（第12章）的推導，係數 $a_C$和 $a_C^\circ$爲
$\begin{aligned} a_C&=\sum_{f\sim nb(C)}\left( \Gamma_f^\phi\frac{E_f}{d_{CF}} + ||\dot m_f^\circ,0|| \right) \\ a_C^\circ&=-\frac{\rho_C^\circ V_C}{\Delta t} \end{aligned}$
CFL條件變爲
$a_C+a_C^\circ\le0 \Rightarrow \sum_{f\sim nb(C)}\left( \Gamma_f^\phi\frac{E_f}{d_{CF}} + ||\dot m_f^\circ,0|| \right) -\frac{\rho_C^\circ V_C}{\Delta t} \le0$
即以下對時間步長的限制關係
$\Delta t \le \frac{\rho_C^\circ V_C}{\displaystyle\sum_{f\sim nb(C)}\left( \Gamma_f^\phi\frac{E_f}{d_{CF}} + ||\dot m_f^\circ,0|| \right) }$
上式一般是顯式瞬態格式的穩定性要求，實際上，以前對一維區域的純對流和純擴散的條件，能夠視做是上式的特例。好比一維擴散問題，有均勻網格單元尺度 $\Delta x$，常密度 $\rho$，均勻擴散係數 $\Gamma^\phi$，上式變爲
$\Delta t \le \frac{\rho_C^\circ V_C} {\displaystyle\sum_{f\sim nb(C)}\left( \Gamma_f^\phi \frac{E_f}{d_{CF}} + ||\dot m_f^\circ,0|| \right) }= \frac{\rho_C^\circ \Delta x_C \Delta y_C} {(\Gamma_e^\phi+\Gamma_w^\phi) \dfrac{\Delta y_C}{\Delta x_C} }= \frac{\rho_C^\circ(\Delta x_C)^2}{2\Gamma_C^\phi}$
而對於純對流的一維問題，使用迎風格式離散對流項，流體從左向右流動，則其變爲了
$\Delta t \le \frac{\rho_C^\circ V_C} {\displaystyle\sum_{f\sim nb(C)}\left( \Gamma_f^\phi \frac{E_f}{d_{CF}} + ||\dot m_f^\circ,0|| \right) }= \frac{\rho_C^\circ \Delta x_C \Delta y_C}{\rho_C^\circ u_C^\circ \Delta y_C }=\frac{\Delta x_C}{u_C^\circ}$
穩定性限制條件是嚴格的，且限制性很強，它迫使在求解瞬態問題時採用很是小的時間步長。這樣作的後果是，儘管每一個時間步的計算消耗很小（相比於每一個時間步要求解方程組系統而言），然而CFL條件使得必須求解不少時間步才能推動到終了時刻。這樣一來，每一個時間步減小的計算量的優點就毫無心義了，由於反而須要更多的時間步來求解了。此外，CFL條件代表，若是經過減少網格尺度來提升空間精度的話，那麼會更進一步地減少爲保證穩定性所能使用的最大時間步長。

下面會看到，這樣的限制條件對於隱式格式來講並不存在，由於其瞬態項的符號老是合適的。

2.3 後向Euler格式（Backward Euler Scheme）

爲了推出後向Euler格式，把在時刻 $t-\Delta t$的函數 $T$值用在時刻 $t$的函數 $T$值作Taylor級數展開，有
$T(t-\Delta t)=T(t)-\frac{\partial T(t)}{\partial t}\Delta t+\frac{\partial^2 T(t)}{\partial t^2}\frac{\Delta t^2}{2!}+...$
整理後，有
$\frac{\partial T(t)}{\partial t}=\frac{T(t)-T(t-\Delta t)}{\Delta t}+\frac{\partial^2 T(t)}{\partial t^2}\frac{\Delta t}{2}+...$
把上式中的 $T$替換成 $(\rho\phi)$，並代入到 $\frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0$，離散方程變爲
$\frac{(\rho_C \phi_C)^{t}-(\rho_C \phi_C)^{t-\Delta t}}{\Delta t}V_C+L(\phi_C^t)=0$
接下來，引入空間離散項的代數方程，並根據以前設置的上標，瞬態標量方程的完整代數形式爲
$(a_C^{\bullet}+a_C) \phi_C +\sum_{F\sim NB(C)}a_F\phi_F =b_C-a_C^\circ\phi_C^\circ$
其中
$\begin{aligned} a_C^{\bullet} &= \frac{\rho_C V_C}{\Delta t} \\ a_C^\circ &= -\frac{\rho_C^\circ V_C}{\Delta t} \end{aligned}$
其架構以下圖所示，顯然，空間算子是在和新時間係數相同的時間層上計算的，因此，求解新時刻的 $\phi$場是須要求解方程組系統的。這種須要求解方程組系統的解法即是隱式格式。

從上式中能夠看出， $a_C$和 $a_C^\circ$的符號是相反的（知足反號準則，即對角係數和鄰居係數符號是相反的），這樣即可保證 $\phi_C$是由其 當前時間步的空間鄰近點的值 和 前一時間步 $t-\Delta t$的時間鄰近點的值 所限定的。這也意味着該格式老是穩定的，而無論用的是何種時間步長，從而能夠用很大的時間步長向前快速推動。縱然如此，其並不是是理想格式，由於其精度較低，這種格式得到的解精確性差，除非用很小的時間步長才好，這讓其使用起來頗爲尷尬。若採用大時間步長來提升計算效率，則得到的解精度不好。若使用小時間步長來得到精確解，則計算效率會很是低下。

2.4 Crank-Nicolson格式

在Crank-Nicolson格式中，爲了得到瞬態項更加精確的近似，採用了在 $t-\Delta t$時刻的函數 $T$值和在 $t+\Delta t$時刻的函數 $T$值的展開形式，一樣是用 $t$時刻的量來展開的，即
$\begin{aligned} T(t+\Delta t)=T(t)+\frac{\partial T(t)}{\partial t}\Delta t+\frac{\partial^2 T(t)}{\partial t^2}\frac{\Delta t^2}{2!}+\frac{\partial^3 T(t)}{\partial t^3}\frac{\Delta t^3}{3!}+... \\ T(t-\Delta t)=T(t)-\frac{\partial T(t)}{\partial t}\Delta t+\frac{\partial^2 T(t)}{\partial t^2}\frac{\Delta t^2}{2!}-\frac{\partial^3 T(t)}{\partial t^3}\frac{\Delta t^3}{3!}+... \end{aligned}$
兩式相減，得
$\frac{\partial T(t)}{\partial t}=\frac{T(t+\Delta t)-T(t-\Delta t)}{2\Delta t}+O(\Delta t^2)$
注意，偏導的精度的階數爲 $O(\Delta t^2)$，由於二階偏導被徹底消掉了。

把上式中的 $T$替換成 $(\rho\phi)$，並代入到 $\frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0$，離散方程變爲
$\frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^{t-\Delta t}}{2\Delta t}V_C+L(\phi_C^t)=0$
接下來，引入空間離散項的代數方程，並根據以前設置的上標，瞬態標量方程的完整代數形式爲
$a_C^{\bullet} \phi_C = b_C-\left(a_C\phi_C^\circ+\sum_{F\sim NB(C)}a_F\phi_F^\circ\right) - a_C^{\circ\circ} \phi_C^{\circ\circ}$
其中
$\begin{aligned} a_C^{\bullet} &= \frac{\rho_C V_C}{2\Delta t} \\ a_C^{\circ\circ} &= -\frac{\rho_C^{\circ\circ} V_C}{2\Delta t} \end{aligned}$
其架構以下圖所示，顯然，這個格式是顯示算法，由於 $(\rho\phi)^{t+\Delta t}$的獲取只用到了舊時刻的值。然而如今須要存儲兩個舊時刻的值，但空間離散算子仍是隻須要其中一箇舊時刻的值就行了。

對於CN格式的穩定性分析，能夠把其最初的方程作些許修改來進行，使用以下近似：
$\phi^\circ\approx\frac{\phi+\phi^{\circ\circ}}{2}$
代數方程變爲
$a_C^{\bullet} \phi_C+0.5\left(a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F\right) = b_C-0.5\left((a_C+2a_C^{\circ\circ})\phi_C^{\circ\circ}+\sum_{F\sim NB(C)}a_F\phi_F^{\circ\circ}\right)$
如此一來，穩定性條件就變成了
$a_C+2a_C^{\circ\circ} \le 0$
對於一維瞬態對流問題，有
$\Delta t\le \frac{2\rho_C^{\circ\circ}V_C}{\dot m_e^\circ}=\frac{2\rho_C^{\circ\circ}\Delta x_C\Delta y_C}{\rho_C^\circ u_C^\circ \Delta y_C}\approx\frac{2\Delta x_C}{|\bold v_e^\circ|}$
其中對流項是用迎風格式離散的。使用前面定義的對流CFL數，上式變爲
$CFL^{conv} \le 2$
這個CFL數的最大限制可要比正向Eluer格式的大了，很是愜意，然而精度的提高則是更有意義，由於其讓精確解的獲取無需經過訴諸於更小的時間步長來實現，尤爲是其把二階導數都從偏差中消掉了。隨後章節中將展開更爲詳細的精確性分析。

2.5 實現細節

CN格式也能夠經過把前向和後向Euler格式相加而推出，以下
$\begin{aligned} & \text{Forward~Euler} \rightarrow \frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^t}{\Delta t}V_C+L(\phi_C^t)=0 \\ & \text{Backward~Euler} \rightarrow \frac{(\rho_C \phi_C)^{t}-(\rho_C \phi_C)^{t-\Delta t}}{\Delta t}V_C+L(\phi_C^t)=0 \end{aligned}$
二者相加
$\begin{aligned} &\frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^t}{\Delta t}V_C+ \frac{(\rho_C \phi_C)^{t}-(\rho_C \phi_C)^{t-\Delta t}}{\Delta t}V_C+2L(\phi_C^t)=0 \\ \rightarrow &\frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^{t-\Delta t}}{2\Delta t}V_C+L(\phi_C^t)=0 \end{aligned}$
即Crank-Nicolson格式

該方程指明瞭CN格式的簡單實現方法，即，在隱式框架下的兩步法來實現。第一步是後向Euler格式，用於隱式找到 $(\rho\phi)^t$
$(\rho_C \phi_C)^{t}+\frac{\Delta t}{V_C}L(\phi_C^t)=(\rho_C \phi_C)^{t-\Delta t}$
而後在第2步中，顯式算得 $t+\Delta t$時刻的CN值
$\begin{aligned} & \frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^t}{\Delta t}V_C=-L(\phi_C^t)= \frac{(\rho_C \phi_C)^{t}-(\rho_C \phi_C)^{t-\Delta t}}{\Delta t}V_C \\ \Rightarrow & (\rho_C \phi_C)^{t+\Delta t}=2(\rho_C \phi_C)^{t}-(\rho_C \phi_C)^{t-\Delta t} \end{aligned}$
在該推導過程當中，假設瞬態時間步長 $\Delta t$分爲了兩個等長的局部時間步長 $\Delta t_{local}$，且 $\Delta t_{local}$爲設置時間步長 $\Delta t$的一半。

值得注意的是，儘管CN格式是二階精度，可是它仍舊是顯式格式，且受限於CFL條件，如前所示。

2.6 Adams-Moulton格式

二階Adams-Moulton格式的發展須要把 $T$在時刻 $t-\Delta t$和 $t-2\Delta t$的值使用Taylor級數展開成 $t$時刻的形式，即
$\begin{aligned} T(t-2\Delta t)&=T(t)-\frac{\partial T(t)}{\partial t}2\Delta t+\frac{\partial^2 T(t)}{\partial t^2}\frac{4\Delta t^2}{2!}+... \\ T(t-\Delta t)&=T(t)-\frac{\partial T(t)}{\partial t}\Delta t+\frac{\partial^2 T(t)}{\partial t^2}\frac{\Delta t^2}{2!}+... \end{aligned}$
從上面兩個式子出發，想辦法把二階偏導項消掉，得
$\frac{\partial T(t)}{\partial t}=\frac{3T(t)-4T(t-\Delta t)+T(t-2\Delta t)}{2\Delta t}$
把上式中的 $T$替換成 $(\rho\phi)$，並代入到 $\frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0$，離散方程變爲
$\frac{3(\rho_C \phi_C)^t-4(\rho_C \phi_C)^{t-\Delta t}+(\rho_C \phi_C)^{t-2\Delta t}}{2\Delta t}V_C+L(\phi_C^t)=0$
接下來，引入空間離散項的代數方程，並根據以前設置的上標，瞬態標量方程的完整代數形式爲
$(a_C^{\bullet}+a_C) \phi_C +\sum_{F\sim NB(C)}a_F\phi_F = b_C-a_C^{\circ}\phi_C^\circ - a_C^{\circ\circ} \phi_C^{\circ\circ}$
其中
$\begin{aligned} a_C^{\bullet} &= \frac{3\rho_C V_C}{2\Delta t} \\ a_C^{\circ} &= -\frac{2\rho_C^{\circ} V_C}{\Delta t} \\ a_C^{\circ\circ} &= \frac{\rho_C^{\circ\circ} V_C}{2\Delta t} \end{aligned}$
顯然， $a_C^{\circ\circ}$係數是正值，因此 $\phi_C^{\circ\circ}$的增長會致使 $\phi_C$的減少。因爲 $a_C^\circ$係數是負的，因此若 $a_C^\circ$較大，則可略微削減該不利影響。如此看來，儘管該格式是穩定的，然而其是無界的，在某些狀況下會出現非物理振盪。

3 有限體積方法

有限體積方法對瞬態項的離散 和 對流項的離散 很是類似，只是積分是在時域單元而非空間單元上進行的，以下圖所示。

對式 $\frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0$在時間區間 $[t-\Delta t/2, t+\Delta t/2]$作積分，即
$\int_{t-\Delta t/2}^{t+\Delta t/2}\frac{\partial (\rho_C \phi_C)}{\partial t}V_Cdt+\int_{t-\Delta t/2}^{t+\Delta t/2}L(\phi_C^t)dt=0$
把 $V_C$做爲常數處理，則第1項變爲面通量差分形式，第2項用積分中值定理來作體積分，得
$V_C(\rho_C \phi_C)^{t+\Delta t/2}-V_C(\rho_C \phi_C)^{t-\Delta t/2}+L(\phi_C^t)\Delta t=0$
該半離散瞬態方程可寫成更加標準的形式，兩邊都除以 $\Delta t$，得
$\frac{(\rho_C \phi_C)^{t+\Delta t/2}-(\rho_C \phi_C)^{t-\Delta t/2}}{\Delta t}V_C+L(\phi_C^t)=0$
爲了獲得徹底離散方程，須要用插值廓線把在 $t+\Delta t/2$和 $t-\Delta t/2$處的面值用 $t$、 $t-\Delta t$等處的單元值來表示。該廓線的選擇可從對流項離散中汲取靈感。選擇不一樣的廓線，毫無疑問會影響到方法的精確性和穩定性。基於此，須要說一下，空間算子的積分是時域二階精度的，可是空間算子自己的精度（空間精度）則是由其離散時所採用的格式所決定的。

無論是使用的哪一種廓線，通量均可用新值和舊值線性化爲
$FluxT=FluxC\phi_C+FluxC^\circ\phi_C^\circ+FluxV$
其中上標 $^\circ$仍舊錶明舊值。完成該線性化處理後，代數方程的相關係數更替爲下式，便可
$\begin{aligned} a_C & \leftarrow a_C+FluxC \\ b_C & \leftarrow b_C-FluxC^\circ\phi_C^\circ-FluxV \end{aligned}$
接下來展現某些廓線形式下的離散過程。

3.1 一階瞬態格式

一階隱式和顯示Euler格式，可分別經過採用迎風和背風瞬態插值廓線來構造，以下。

3.2 一階隱式Euler格式

使用一階迎風插值廓線，即獲得瞬態一階隱式Euler格式。如上圖所示，在時間單元面處的 $\rho\phi$值，等於迎風單元形心值，即
$(\rho_C \phi_C)^{t+\Delta t/2}=(\rho_C \phi_C)^{t}\quad\quad (\rho_C \phi_C)^{t-\Delta t/2}=(\rho_C \phi_C)^{t-\Delta t}$
那麼半離散方程爲
$\frac{(\rho_C \phi_C)^{t}-(\rho_C \phi_C)^{t-\Delta t}}{\Delta t}V_C+L(\phi_C^t)=0$
即一階隱式Euler格式，線性化後的相關係數爲
$\begin{aligned} FluxC &= \frac{\rho_C V_C}{\Delta t} \\ FluxC^\circ &= -\frac{\rho_C^\circ V_C}{\Delta t} \\ FluxV &= 0 \end{aligned}$

3.2.1 數值擴散

因爲這是一階格式，那麼，基於前面對流格式中得到的知識，該格式會產生數值擴散。可經過在時刻 $t$的Taylor級數展開，將其恢復到最初控制方程的形式，來肯定數值擴散的值，即
$(\rho \phi)^{t-\Delta t}=(\rho \phi)^t-\left.\frac{\partial (\rho \phi)}{\partial t}\right|_t\Delta t +\left.\frac{\partial^2 (\rho \phi)}{\partial t^2}\right|_t\frac{\Delta t^2}{2}+O(\Delta t^3)$
調整爲
$\frac{(\rho \phi)^{t}-(\rho \phi)^{t-\Delta t}}{\Delta t}= \left.\frac{\partial (\rho \phi)}{\partial t}\right|_t -\left(\frac{\Delta t}{2}\right)\left.\frac{\partial^2 (\rho \phi)}{\partial t^2}\right|_t-O(\Delta t^2)$
把上式代入到離散格式中，得
$\left.\frac{\partial (\rho \phi)}{\partial t}\right|_t+\frac{1}{V_C}L(\phi_C^t)= \underbrace{\left(\frac{\Delta t}{2}\right)\left.\frac{\partial^2 (\rho \phi)}{\partial t^2}\right|_t}_\text{Numerical~diffusion~term}+O(\Delta t^2)$
實際上，給方程添加了個數值擴散項，且其與時間步長成正比，這跟對流格式中的迎風格式是類似的。因此呢，這個格式是無條件穩定的，其可對大時間步長產生穩定解。

3.3 一階顯示Euler格式

使用一階背風插值廓線，便可獲得瞬態一階顯示Euler格式，如上圖所示，在時間單元面處的 $\rho\phi$值，等於背風單元形心值，即
$(\rho_C \phi_C)^{t+\Delta t/2}=(\rho_C \phi_C)^{t+\Delta t}\quad\quad (\rho_C \phi_C)^{t-\Delta t/2}=(\rho_C \phi_C)^{t}$
那麼半離散方程爲
$\frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^{t}}{\Delta t}V_C+L(\phi_C^t)=0$
即一階顯式Euler格式，線性化後的相關係數爲
$\begin{aligned} FluxC &= \frac{\rho_C V_C}{\Delta t} \\ FluxC^\circ &= -\frac{\rho_C^\circ V_C}{\Delta t} \\ FluxV &= 0 \end{aligned}$
注意，如今新時刻是 $t+\Delta t$，空間算子則是在舊時刻 $t$上衡量的，這樣一來，右端項就徹底已知了，能夠直接用來獲得 $t+\Delta t$時刻的 $\rho\phi$，而不須要求解線性代數方程組。這是顯式格式。

3.3.1 數值反擴散

基於 $t$時刻作Taylor展開
$(\rho \phi)^{t+\Delta t}=(\rho \phi)^t+\left.\frac{\partial (\rho \phi)}{\partial t}\right|_t\Delta t +\left.\frac{\partial^2 (\rho \phi)}{\partial t^2}\right|_t\frac{\Delta t^2}{2}+O(\Delta t^3)$
調整爲
$\frac{(\rho \phi)^{t+\Delta t}-(\rho \phi)^{t}}{\Delta t}= \left.\frac{\partial (\rho \phi)}{\partial t}\right|_t +\left(\frac{\Delta t}{2}\right)\left.\frac{\partial^2 (\rho \phi)}{\partial t^2}\right|_t +O(\Delta t^2)$
把上式代入到離散格式中，得
$\left.\frac{\partial (\rho \phi)}{\partial t}\right|_t+\frac{1}{V_C}L(\phi_C^t)= -\underbrace{\left(\frac{\Delta t}{2}\right)\left.\frac{\partial^2 (\rho \phi)}{\partial t^2}\right|_t}_\text{Numerical~anti-diffusion~term}+O(\Delta t^2)$
如今二階差分項的符號是負的，相似於負擴散或反擴散，廓線上有壓縮效應，這與對流格式中的背風格式很是像。反擴散項也是和時間步長呈正比關係的。當其與迎風對流格式聯合起來，且Courant數爲1的時候，對流項的數值擴散和 $CFL^{conv}=1$的顯式Euler格式的數值反擴散，二者幅值相等且符號相反，所以二者相互抵消，產生了近乎精確的解。儘管如此，確保 $CFL^{conv}=1$是不切實際的，並且簡單的一維網格也並不是是實際問題的再現。

與反擴散效應相關的問題是數值不穩定性，其隨着 $\Delta t$的增長而增長，所以須要對時間步長施加較強的限制，能夠經過負鄰近係數準則來評估最大穩定時間步。

3.4 二階瞬態Euler格式

與對流格式相似，二階瞬態格式可經過線性插值廓線來構造。可選用對稱廓線（中心差分）來產生Crank-Nicolson（CN）格式，也可選擇迎風（二階迎風格式）來產生Adams-Moulton格式，即，著名的隱式格式，二階迎風Euler（Second Order Upwind Euler（SOUE））。

3.5 Crank-Nicolson（中心差分廓線）

採用在迎風和背風節點之間的線性插值來計算 $\rho\phi$值，即可得到上圖所示的Crank-Nicolson格式。

對於均勻時間步，其數學表達式爲
$\begin{aligned} & (\rho_C \phi_C)^{t+\Delta t/2}=\frac12(\rho_C \phi_C)^{t+\Delta t}+\frac12(\rho_C \phi_C)^{t} \\ & (\rho_C \phi_C)^{t-\Delta t/2}=\frac12(\rho_C \phi_C)^{t}+\frac12(\rho_C \phi_C)^{t-\Delta t} \end{aligned}$
那麼半離散方程爲
$\frac{(\rho_C \phi_C)^{t+\Delta t}-(\rho_C \phi_C)^{t-\Delta t}}{2\Delta t}V_C+L(\phi_C^t)=0$
即一階顯式Euler格式，線性化後的相關係數爲
$\begin{array}{rl}{\displaystyle FluxC}& {\displaystyle =\frac{{\rho }_C{V}_C}{2\mathrm{\Delta }t}}\\ {\displaystyle Flux{C}^{\circ }}& {\displaystyle =0}\\ {\displaystyle FluxV}& {\displaystyle =-\frac{{\rho }_C^{\circ \circ }{V}_C}{2\mathrm{\Delta }t}{\varphi }_{}^{}}\end{array}$