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}$