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

前面章节的讨论中均假设是稳态条件，所以无需对瞬态项进行离散。如果考虑瞬态现象的话，则相当于给问题增加了一个新的维度。然而，由于瞬态项的性质是抛物型的，故而不需要再额外定义时间域上的场，即，不像在空间域上需要用${\varphi }_i,i=1,2,...,N$那样子来离散节点值。一般来说，瞬态项的离散只需要额外存储一到两个时间层上的变量场即可，由所选择格式的数值精度而定。与稳态问题的另一个不同点是，瞬态系统是用时间步长推进方式来模拟的，即，以$t=t_0$时刻的初始条件开始，求解算法向前推进并找到$t_1=t_0+\Delta t_1$时刻的解，这个找到的解再作为初始条件去找寻$t_2=t_1+\Delta t_2$时刻的解，该过程不断重复直到达到所需的时刻为止。本章的重点是瞬态项离散的技术，将展示两种发展瞬态格式的方法。其一是使用Taylor展开来把瞬态项展开成节点值的形式，这在有限差分方法中非常奏效；其二是有限体积方法中常用的伪时间单元方法，和在对流项中的伪节点非常类似。将展示一些瞬态格式，并讨论它们的特性。

1 引言

对于瞬态模拟而言，控制方程的离散是在空间域和时间域上进行的。空间域上的空间离散和稳态问题一样，时间离散则要设置一个时间坐标轴，沿着该坐标轴来计算瞬态项的导数（有限差分方法）或积分（有限体积方法）值。

一般而言，某个变量$\varphi$的瞬态特性的表达式，或其时间演化，是由下述形式的方程控制的
$$\frac{\partial (\rho \varphi )}{\partial t}+L(\varphi )=0$$
其中$L(\varphi )$为空间算子，其包含了所有的非瞬态项（扩散项、对流项、源项，等），$\partial (\rho \varphi )/\partial t$为瞬态算子，两者如下图所示。

在单元$C$上对上式做积分，得
$${\int }_{V_C}\frac{\partial (\rho \varphi )}{\partial t}dV+{\int }_{V_C}L(\varphi )dV=0$$
空间离散后，变为
$$\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_C+L({\varphi }_C^t)=0$$
其中$V_C$为离散单元的体积，$L({\varphi }_C^t)$为空间算子在某参考时刻$t$的离散形式，可写成如下代数方程
$$L({\varphi }_C^t)=a_C{\varphi }_C^t+\sum _{F\sim NB(C)}{a}_F{\varphi }_F^t-b_C$$
在上上式中，当$t\to \infty$时其变为稳态离散方程。通过时间推进也能达到稳态，即，${\varphi }_C^{t+\Delta t}={\varphi }_C^t$。这保证了瞬态问题在到达稳态时所获得的解，和直接由稳态问题求得的解，是一样的。

在瞬态项的离散中，最经典的法子是用有限差分法的思想，把$\partial (\rho \varphi )/\partial t$做Taylor级数展开，得到用离散点的值所表示的导数项形式。在本章中，将展示另一种与有限体积方法更为契合的方法，对$\partial (\rho \varphi )/\partial t$在时间单元上做积分，并转化为面通量的形式，这和对流格式的处理方法很像，只是现在离散是沿着时间轴进行的。

2 有限差分方法

如图，在瞬态空间中的网格是结构化的，所以对瞬态项使用有限差分方法是非常稀松平常的操作。在该方法中，空间算子$L(\varphi )$是在时刻$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 \varphi )$，并代入到$\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_C+L({\varphi }_C^t)=0$，离散方程变为
$$\frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^t}{\Delta t}{V}_C+L({\varphi }_C^t)=0$$

该瞬态格式如上图所示，计算$t+\Delta t$时刻的$({\rho }_C{\varphi }_C)$是不需要求解方程组系统的。因为所有的空间项都是在旧时刻$t$上衡量的，所以$t+\Delta t$时刻的${\varphi }_C$可直接由上个时间步的值显式算出。这样的格式属于显式瞬态格式。显式瞬态格式的主要特点是通过在时域推进便可获得解，而无需在每个时间层求解方程组系统。这样一来，求解效率非常高，且对计算网格的并行处理十分方便。然而鲜有商业软件采用该方法，因为其重大缺陷就是时间步长$\Delta t$极为受限，这将在下一小节详细展开。

把空间算子的离散代数方程代入到上式，可得完整的代数方程为
$$a_C^{t+\Delta t}{\varphi }_C^{t+\Delta t}+a_C^t{\varphi }_C^t=b_C-\left(a_C{\varphi }_C^t+\sum _{F\sim NB(C)}{a}_F{\varphi }_F^t\right)$$
其中
$$\begin{array}{rl}{a}_C^{t+\Delta t}& =\frac{{\rho }_C^{t+\Delta t}{V}_C}{\Delta t}\\ a_C^t& =-\frac{{\rho }_C^t{V}_C}{\Delta t}\end{array}$$
在上述方程中，$a_C^{t+\Delta t}$和$a_C^t$为对角系数，源自于瞬态项的离散，${\varphi }_C^{t+\Delta t}$和${\varphi }_C^t$为时间层$t+\Delta t$和$t$上的值，而$a_C$、$a_F$、$b_C$则是空间离散后的系数。

为了简化标记，在本章中，前一个时间步的变量值用上标${}^{\circ }$标记，前两个时间步的变量值用上标${}^{\circ \circ }$标记，若是没有上标，则表示当前时间步的变量值，除了在非定常项中乘到${\varphi }_C$上的系数是用上标${}^{\bullet }$来标记的。基于这些标记，上式写作
$$a_C^{\bullet }{\varphi }_C+a_C^{\circ }{\varphi }_C^{\circ }=b_C-\left(a_C{\varphi }_C^{\circ }+\sum _{F\sim NB(C)}{a}_F{\varphi }_F^{\circ }\right)$$
其中
$$\begin{array}{rl}{a}_C^{\bullet }& =\frac{{\rho }_C{V}_C}{\Delta t}\\ a_C^{\circ }& =-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\end{array}$$
改写为如下形式
$${\varphi }_C=\frac{b_C-\left[(a_C+a_C^{\circ }){\varphi }_C^{\circ }+\sum _{F\sim NB(C)}{a}_F{\varphi }_F^{\circ }\right]}{a_C^{\bullet }}$$
显然，当前时间步的$\varphi$值是由显式关系算出来的，不需要求解方程组系统。

2.2 前向Euler格式的稳定性

数值格式的收敛性和稳定性最初是由Courant、Friedrichs、Lewy所探究的，他们表明，为了让差分方程的解收敛到原偏微分方程的解，数值格式中必须用到影响解的初始数据所包含的所有信息。该要求随后就变成了著名的CFL条件。

实际上，CFL条件可以被简单解释为，把系数所要满足的基本规则之一，即，反符号准则，扩展成把瞬态系数也包含在内。这样，正如${\varphi }_F$是${\varphi }_C$的“空间”邻居，${\varphi }_C^{\circ }$也是${\varphi }_C$的“时间”邻居，那么反符号准则对他俩应该平等适用。注意对角线系数现在是$a_C^{\bullet }$，其“时间”邻居的系数是$(a_C+a_C^{\circ })$，反符号准则变为
$$a_C+a_C^{\circ }\le 0$$

2.2.1 瞬态-对流情况的稳定性

如上图，对于一维纯对流问题，流动方向向右，使用迎风格式把变量插值到单元面上去，则单元$C$的离散方程中的系数$a_C$和$a_C^{\circ }$为
$$\begin{array}{rl}{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{array}$$
因此，CFL条件为
$$a_C+a_C^{\circ }\le 0\Rightarrow {\rho }_C^{\circ }{u}_C^{\circ }\Delta y_C-\frac{{\rho }_C^{\circ }\Delta x_C\Delta y_C}{\Delta t}\le 0$$
即
$$\Delta t\le \frac{\Delta x_C}{u_C^{\circ }}$$
对于对流控制的流动，定义CFL数为
$$CF{L}^{conv}=\frac{|{\mathbf{v}}_C^{\circ }|\Delta t}{\Delta x_C}$$
这意味着为了数值稳定性，CFL数应该满足
$$CF{L}^{conv}\le 1$$

2.2.2 瞬态-扩散情况的稳定性

对于纯扩散问题，CFL数的表达式是不同的。因此，分析如上图所示的一维纯扩散问题。

使用线性插值廓线，单元$C$的离散方程中的系数$a_C$和$a_C^{\circ }$为
$$\begin{array}{rl}{a}_C& =\frac{{\Gamma }_e^{\varphi }\Delta y_C}{\delta x_e}+\frac{{\Gamma }_w^{\varphi }\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{array}$$
因此，CFL条件需要
$$a_C+a_C^{\circ }\le 0\Rightarrow \frac{{\Gamma }_e^{\varphi }\Delta y_C}{\delta x_e}+\frac{{\Gamma }_w^{\varphi }\Delta y_C}{\delta x_w}-\frac{{\rho }_C^{\circ }\Delta x_C\Delta y_C}{\Delta t}\le 0$$
即
$$\Delta t\le \frac{{\rho }_C^{\circ }\Delta x_C}{\frac{{\Gamma }_e^{\varphi }}{\delta x_e}+\frac{{\Gamma }_w^{\varphi }}{\delta x_w}}$$
若网格均匀，且扩散系数为常数，则上式变为
$$\Delta t\le \frac{{\rho }_C^{\circ }(\Delta x_C{)}^2}{2{\Gamma }_C^{\varphi }}$$
对于扩散控制的问题，定义CFL数为
$$CF{L}^{diff}=\frac{{\Gamma }_C^{\varphi }\Delta t}{{\rho }_C^{\circ }(\Delta x_C{)}^2}$$
这意味着为了数值稳定性，CFL数应该满足
$$CF{L}^{diff}\le \frac{1}{2}$$

2.2.3 瞬态-对流-扩散情况的稳定性

对于多维非定常对流扩散问题，基于上一章（第12章）的推导，系数$a_C$和$a_C^{\circ }$为
$$\begin{array}{rl}{a}_C& =\sum _{f\sim nb(C)}\left({\Gamma }_f^{\varphi }\frac{E_f}{d_{CF}}+||{\dot{m}}_f^{\circ },0||\right)\\ a_C^{\circ }& =-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\end{array}$$
CFL条件变为
$$a_C+a_C^{\circ }\le 0\Rightarrow \sum _{f\sim nb(C)}\left({\Gamma }_f^{\varphi }\frac{E_f}{d_{CF}}+||{\dot{m}}_f^{\circ },0||\right)-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\le 0$$
即如下对时间步长的限制关系
$$\Delta t\le \frac{{\rho }_C^{\circ }{V}_C}{\sum _{f\sim nb(C)}\left({\Gamma }_f^{\varphi }\frac{E_f}{d_{CF}}+||{\dot{m}}_f^{\circ },0||\right)}$$
上式通常是显式瞬态格式的稳定性要求，实际上，之前对一维区域的纯对流和纯扩散的条件，可以视作是上式的特例。比如一维扩散问题，有均匀网格单元尺度$\Delta x$，常密度$\rho$，均匀扩散系数${\Gamma }^{\varphi }$，上式变为
$$\Delta t\le \frac{{\rho }_C^{\circ }{V}_C}{\sum _{f\sim nb(C)}\left({\Gamma }_f^{\varphi }\frac{E_f}{d_{CF}}+||{\dot{m}}_f^{\circ },0||\right)}=\frac{{\rho }_C^{\circ }\Delta x_C\Delta y_C}{({\Gamma }_e^{\varphi }+{\Gamma }_w^{\varphi })\frac{\Delta y_C}{\Delta x_C}}=\frac{{\rho }_C^{\circ }(\Delta x_C{)}^2}{2{\Gamma }_C^{\varphi }}$$
而对于纯对流的一维问题，使用迎风格式离散对流项，流体从左向右流动，则其变为了
$$\Delta t\le \frac{{\rho }_C^{\circ }{V}_C}{\sum _{f\sim nb(C)}\left({\Gamma }_f^{\varphi }\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 \varphi )$，并代入到$\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_C+L({\varphi }_C^t)=0$，离散方程变为
$$\frac{({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{\Delta t}{V}_C+L({\varphi }_C^t)=0$$
接下来，引入空间离散项的代数方程，并根据之前设置的上标，瞬态标量方程的完整代数形式为
$$(a_C^{\bullet }+a_C){\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F=b_C-a_C^{\circ }{\varphi }_C^{\circ }$$
其中
$$\begin{array}{rl}{a}_C^{\bullet }& =\frac{{\rho }_C{V}_C}{\Delta t}\\ a_C^{\circ }& =-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\end{array}$$
其架构如下图所示，显然，空间算子是在和新时间系数相同的时间层上计算的，所以，求解新时刻的$\varphi$场是需要求解方程组系统的。这种需要求解方程组系统的解法便是隐式格式。

从上式中可以看出，$a_C$和$a_C^{\circ }$的符号是相反的（满足反号准则，即对角系数和邻居系数符号是相反的），这样便可保证${\varphi }_C$是由其 当前时间步的空间邻近点的值 和 前一时间步$t-\Delta t$的时间邻近点的值 所限定的。这也意味着该格式总是稳定的，而不管用的是何种时间步长，从而可以用很大的时间步长向前快速推进。纵然如此，其并非是理想格式，因为其精度较低，这种格式获得的解精确性差，除非用很小的时间步长才好，这让其使用起来颇为尴尬。若采用大时间步长来提高计算效率，则获得的解精度很差。若使用小时间步长来获得精确解，则计算效率会非常低下。

2.4 Crank-Nicolson格式

在Crank-Nicolson格式中，为了获得瞬态项更加精确的近似，采用了在$t-\Delta t$时刻的函数$T$值和在$t+\Delta t$时刻的函数$T$值的展开形式，同样是用$t$时刻的量来展开的，即
$$\begin{array}{r}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{array}$$
两式相减，得
$$\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 \varphi )$，并代入到$\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_C+L({\varphi }_C^t)=0$，离散方程变为
$$\frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{2\Delta t}{V}_C+L({\varphi }_C^t)=0$$
接下来，引入空间离散项的代数方程，并根据之前设置的上标，瞬态标量方程的完整代数形式为
$$a_C^{\bullet }{\varphi }_C=b_C-\left(a_C{\varphi }_C^{\circ }+\sum _{F\sim NB(C)}{a}_F{\varphi }_F^{\circ }\right)-a_C^{\circ \circ }{\varphi }_C^{\circ \circ }$$
其中
$$\begin{array}{rl}{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{array}$$
其架构如下图所示，显然，这个格式是显示算法，因为$(\rho \varphi {)}^{t+\Delta t}$的获取只用到了旧时刻的值。然而现在需要存储两个旧时刻的值，但空间离散算子还是只需要其中一个旧时刻的值就好了。

对于CN格式的稳定性分析，可以把其最初的方程做些许修改来进行，使用如下近似：
$${\varphi }^{\circ }\approx \frac{\varphi +{\varphi }^{\circ \circ }}{2}$$
代数方程变为
$$a_C^{\bullet }{\varphi }_C+0.5\left(a_C{\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F\right)=b_C-0.5\left((a_C+2a_C^{\circ \circ }){\varphi }_C^{\circ \circ }+\sum _{F\sim NB(C)}{a}_F{\varphi }_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}{|{\mathbf{v}}_e^{\circ }|}$$
其中对流项是用迎风格式离散的。使用前面定义的对流CFL数，上式变为
$$CF{L}^{conv}\le 2$$
这个CFL数的最大限制可要比正向Eluer格式的大了，非常惬意，然而精度的提升则是更有意义，因为其让精确解的获取无需通过诉诸于更小的时间步长来实现，尤其是其把二阶导数都从误差中消掉了。随后章节中将展开更为详细的精确性分析。

2.5 实现细节

CN格式也可以通过把前向和后向Euler格式相加而推出，如下
$$\begin{array}{rl}& \text{Forward Euler}\to \frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^t}{\Delta t}{V}_C+L({\varphi }_C^t)=0\\ & \text{Backward Euler}\to \frac{({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{\Delta t}{V}_C+L({\varphi }_C^t)=0\end{array}$$
两者相加
$$\begin{array}{rl}& \frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^t}{\Delta t}{V}_C+\frac{({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{\Delta t}{V}_C+2L({\varphi }_C^t)=0\\ \to & \frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{2\Delta t}{V}_C+L({\varphi }_C^t)=0\end{array}$$
即Crank-Nicolson格式

该方程指明了CN格式的简单实现方法，即，在隐式框架下的两步法来实现。第一步是后向Euler格式，用于隐式找到$(\rho \varphi {)}^t$
$$({\rho }_C{\varphi }_C{)}^t+\frac{\Delta t}{V_C}L({\varphi }_C^t)=({\rho }_C{\varphi }_C{)}^{t-\Delta t}$$
然后在第2步中，显式算得$t+\Delta t$时刻的CN值
$$\begin{array}{rl}& \frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^t}{\Delta t}{V}_C=-L({\varphi }_C^t)=\frac{({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{\Delta t}{V}_C\\ \Rightarrow & ({\rho }_C{\varphi }_C{)}^{t+\Delta t}=2({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-\Delta t}\end{array}$$
在该推导过程中，假设瞬态时间步长$\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{array}{rl}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{array}$$
从上面两个式子出发，想办法把二阶偏导项消掉，得
$$\frac{\partial T(t)}{\partial t}=\frac{3T(t)-4T(t-\Delta t)+T(t-2\Delta t)}{2\Delta t}$$
把上式中的$T$替换成$(\rho \varphi )$，并代入到$\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_C+L({\varphi }_C^t)=0$，离散方程变为
$$\frac{3({\rho }_C{\varphi }_C{)}^t-4({\rho }_C{\varphi }_C{)}^{t-\Delta t}+({\rho }_C{\varphi }_C{)}^{t-2\Delta t}}{2\Delta t}{V}_C+L({\varphi }_C^t)=0$$
接下来，引入空间离散项的代数方程，并根据之前设置的上标，瞬态标量方程的完整代数形式为
$$(a_C^{\bullet }+a_C){\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F=b_C-a_C^{\circ }{\varphi }_C^{\circ }-a_C^{\circ \circ }{\varphi }_C^{\circ \circ }$$
其中
$$\begin{array}{rl}{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{array}$$
显然，$a_C^{\circ \circ }$系数是正值，所以${\varphi }_C^{\circ \circ }$的增加会导致${\varphi }_C$的减小。由于$a_C^{\circ }$系数是负的，所以若$a_C^{\circ }$较大，则可略微削减该不利影响。如此看来，尽管该格式是稳定的，然而其是无界的，在某些情况下会出现非物理振荡。

3 有限体积方法

有限体积方法对瞬态项的离散 和 对流项的离散 非常相似，只是积分是在时域单元而非空间单元上进行的，如下图所示。

对式$\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_C+L({\varphi }_C^t)=0$在时间区间$[t-\Delta t/2,t+\Delta t/2]$做积分，即
$${\int }_{t-\Delta t/2}^{t+\Delta t/2}\frac{\partial ({\rho }_C{\varphi }_C)}{\partial t}{V}_{C}dt+{\int }_{t-\Delta t/2}^{t+\Delta t/2}L({\varphi }_C^t)dt=0$$
把$V_C$作为常数处理，则第1项变为面通量差分形式，第2项用积分中值定理来做体积分，得
$$V_C({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}-V_C({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}+L({\varphi }_C^t)\Delta t=0$$
该半离散瞬态方程可写成更加标准的形式，两边都除以$\Delta t$，得
$$\frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}-({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}}{\Delta t}{V}_C+L({\varphi }_C^t)=0$$
为了得到完全离散方程，需要用插值廓线把在$t+\Delta t/2$和$t-\Delta t/2$处的面值用$t$、$t-\Delta t$等处的单元值来表示。该廓线的选择可从对流项离散中汲取灵感。选择不同的廓线，毫无疑问会影响到方法的精确性和稳定性。基于此，需要说一下，空间算子的积分是时域二阶精度的，但是空间算子本身的精度（空间精度）则是由其离散时所采用的格式所决定的。

不管是使用的哪种廓线，通量都可用新值和旧值线性化为
$$FluxT=FluxC{\varphi }_C+Flux{C}^{\circ }{\varphi }_C^{\circ }+FluxV$$
其中上标${}^{\circ }$仍旧代表旧值。完成该线性化处理后，代数方程的相关系数更替为下式，即可
$$\begin{array}{rl}{a}_C& \leftarrow a_C+FluxC\\ b_C& \leftarrow b_C-Flux{C}^{\circ }{\varphi }_C^{\circ }-FluxV\end{array}$$
接下来展示某些廓线形式下的离散过程。

3.1 一阶瞬态格式

一阶隐式和显示Euler格式，可分别通过采用迎风和背风瞬态插值廓线来构造，如下。

3.2 一阶隐式Euler格式

使用一阶迎风插值廓线，即得到瞬态一阶隐式Euler格式。如上图所示，在时间单元面处的$\rho \varphi$值，等于迎风单元形心值，即
$$({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}=({\rho }_C{\varphi }_C{)}^t({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}=({\rho }_C{\varphi }_C{)}^{t-\Delta t}$$
那么半离散方程为
$$\frac{({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{\Delta t}{V}_C+L({\varphi }_C^t)=0$$
即一阶隐式Euler格式，线性化后的相关系数为
$$\begin{array}{rl}FluxC& =\frac{{\rho }_C{V}_C}{\Delta t}\\ Flux{C}^{\circ }& =-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\\ FluxV& =0\end{array}$$

3.2.1 数值扩散

由于这是一阶格式，那么，基于前面对流格式中获得的知识，该格式会产生数值扩散。可通过在时刻$t$的Taylor级数展开，将其恢复到最初控制方程的形式，来确定数值扩散的值，即
$$(\rho \varphi {)}^{t-\Delta t}=(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}+O(\Delta t^3)$$
调整为
$$\frac{(\rho \varphi {)}^t-(\rho \varphi {)}^{t-\Delta t}}{\Delta t}={\frac{\partial (\rho \varphi )}{\partial t}|}_t-\left(\frac{\Delta t}{2}\right){\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t-O(\Delta t^2)$$
把上式代入到离散格式中，得
$${\frac{\partial (\rho \varphi )}{\partial t}|}_t+\frac{1}{V_C}L({\varphi }_C^t)=\underset{\text{Numerical diffusion term}}{\underbrace{\left(\frac{\Delta t}{2}\right){\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t}}+O(\Delta t^2)$$
实际上，给方程添加了个数值扩散项，且其与时间步长成正比，这跟对流格式中的迎风格式是相似的。所以呢，这个格式是无条件稳定的，其可对大时间步长产生稳定解。

3.3 一阶显示Euler格式

使用一阶背风插值廓线，即可得到瞬态一阶显示Euler格式，如上图所示，在时间单元面处的$\rho \varphi$值，等于背风单元形心值，即
$$({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}=({\rho }_C{\varphi }_C{)}^{t+\Delta t}({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}=({\rho }_C{\varphi }_C{)}^t$$
那么半离散方程为
$$\frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^t}{\Delta t}{V}_C+L({\varphi }_C^t)=0$$
即一阶显式Euler格式，线性化后的相关系数为
$$\begin{array}{rl}FluxC& =\frac{{\rho }_C{V}_C}{\Delta t}\\ Flux{C}^{\circ }& =-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\\ FluxV& =0\end{array}$$
注意，现在新时刻是$t+\Delta t$，空间算子则是在旧时刻$t$上衡量的，这样一来，右端项就完全已知了，可以直接用来得到$t+\Delta t$时刻的$\rho \varphi$，而不需要求解线性代数方程组。这是显式格式。

3.3.1 数值反扩散

基于$t$时刻做Taylor展开
$$(\rho \varphi {)}^{t+\Delta t}=(\rho \varphi {)}^t+{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}+O(\Delta t^3)$$
调整为
$$\frac{(\rho \varphi {)}^{t+\Delta t}-(\rho \varphi {)}^t}{\Delta t}={\frac{\partial (\rho \varphi )}{\partial t}|}_t+\left(\frac{\Delta t}{2}\right){\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t+O(\Delta t^2)$$
把上式代入到离散格式中，得
$${\frac{\partial (\rho \varphi )}{\partial t}|}_t+\frac{1}{V_C}L({\varphi }_C^t)=-\underset{\text{Numerical anti-diffusion term}}{\underbrace{\left(\frac{\Delta t}{2}\right){\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t}}+O(\Delta t^2)$$
现在二阶差分项的符号是负的，类似于负扩散或反扩散，廓线上有压缩效应，这与对流格式中的背风格式非常像。反扩散项也是和时间步长呈正比关系的。当其与迎风对流格式联合起来，且Courant数为1的时候，对流项的数值扩散和$CF{L}^{conv}=1$的显式Euler格式的数值反扩散，两者幅值相等且符号相反，因此两者相互抵消，产生了近乎精确的解。尽管如此，确保$CF{L}^{conv}=1$是不切实际的，而且简单的一维网格也并非是实际问题的再现。

与反扩散效应相关的问题是数值不稳定性，其随着$\Delta t$的增加而增加，因此需要对时间步长施加较强的限制，可以通过负邻近系数准则来评估最大稳定时间步。

3.4 二阶瞬态Euler格式

与对流格式类似，二阶瞬态格式可通过线性插值廓线来构造。可选用对称廓线（中心差分）来产生Crank-Nicolson（CN）格式，也可选择迎风（二阶迎风格式）来产生Adams-Moulton格式，即，著名的隐式格式，二阶迎风Euler（Second Order Upwind Euler（SOUE））。

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

采用在迎风和背风节点之间的线性插值来计算$\rho \varphi$值，便可获得上图所示的Crank-Nicolson格式。

对于均匀时间步，其数学表达式为
$$\begin{array}{rl}& ({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}=\frac{1}{2}({\rho }_C{\varphi }_C{)}^{t+\Delta t}+\frac{1}{2}({\rho }_C{\varphi }_C{)}^t\\ & ({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}=\frac{1}{2}({\rho }_C{\varphi }_C{)}^t+\frac{1}{2}({\rho }_C{\varphi }_C{)}^{t-\Delta t}\end{array}$$
那么半离散方程为
$$\frac{({\rho }_C{\varphi }_C{)}^{t+\Delta t}-({\rho }_C{\varphi }_C{)}^{t-\Delta t}}{2\Delta t}{V}_C+L({\varphi }_C^t)=0$$
即一阶显式Euler格式，线性化后的相关系数为
$$\begin{array}{rl}FluxC& =\frac{{\rho }_C{V}_C}{2\Delta t}\\ Flux{C}^{\circ }& =0\\ FluxV& =-\frac{{\rho }_C^{\circ \circ }{V}_C}{2\Delta t}{\varphi }_C^{\circ \circ }\end{array}$$
该格式是由$t$和$t-\Delta t$时刻的值直接显式算得$t+\Delta t$时刻的值，所以是显式格式。其稳定性受CFL限制。

跟在有限差分框架中讲过的那样，上式可改写成两步走的样子，即，一阶隐式Euler步，和，插值形式的修正显式Euler步。

3.5.1 数值精度

基于$t$时刻做Taylor展开
$$\begin{gathered} (\rho \varphi {)}^{t+\Delta t}=(\rho \varphi {)}^t+{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}+{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^3}{6}+O(\Delta t^4) \\ (\rho \varphi {)}^{t-\Delta t}=(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}-{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^3}{6}+O(\Delta t^4) \end{gathered}$$
两式相减，得
$$\frac{(\rho \varphi {)}^{t+\Delta t}-(\rho \varphi {)}^{t-\Delta t}}{2\Delta t}={\frac{\partial (\rho \varphi )}{\partial t}|}_t+{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^2}{6}-O(\Delta t^3)$$
把上式代入到离散格式中，得
$${\frac{\partial (\rho \varphi )}{\partial t}|}_t+\frac{1}{V_C}L({\varphi }_C^t)=-{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^2}{6}+O(\Delta t^3)$$
表明该格式具有二阶精度，第三项是色散项，会导致不稳定性。

3.6 二阶迎风Euler（SOUE）格式

使用上图所示的二阶迎风插值廓线，面$\rho \varphi$值的近似为
$$\begin{array}{rl}& ({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}=\frac{3}{2}({\rho }_C{\varphi }_C{)}^t-\frac{1}{2}({\rho }_C{\varphi }_C{)}^{t-\Delta t}\\ & ({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}=\frac{3}{2}({\rho }_C{\varphi }_C{)}^{t-\Delta t}-\frac{1}{2}({\rho }_C{\varphi }_C{)}^{t-2\Delta t}\end{array}$$
那么半离散方程为
$$\frac{3({\rho }_C{\varphi }_C{)}^t-4({\rho }_C{\varphi }_C{)}^{t-\Delta t}+({\rho }_C{\varphi }_C{)}^{t-2\Delta t}}{2\Delta t}{V}_C+L({\varphi }_C^t)=0$$
这是隐式二阶迎风Euler（SOUE）格式，该格式中需要存储两个旧时刻的$\rho \varphi$值，其线性化系数为
$$\begin{array}{rl}FluxC& =\frac{3{\rho }_C{V}_C}{2\Delta t}\\ Flux{C}^{\circ }& =-\frac{2{\rho }_C^{\circ }{V}_C}{\Delta t}\\ FluxV& =\frac{{\rho }_C^{\circ \circ }{V}_C{\varphi }_C^{\circ \circ }}{2\Delta t}\end{array}$$

3.6.1 数值精度

基于$t$时刻做Taylor展开
$$\begin{array}{rl}(\rho \varphi {)}^{t-\Delta t}& =(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}-{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^3}{6}+O(\Delta t^4)\\ (\rho \varphi {)}^{t-2\Delta t}& =(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_{t}2\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{4\Delta t^2}{2}-{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{8\Delta t^3}{6}+O(\Delta t^4)\end{array}$$
第一个式子乘以4并与第二个式子相减去，得
$$\frac{3({\rho }_C{\varphi }_C{)}^t-4({\rho }_C{\varphi }_C{)}^{t-\Delta t}+({\rho }_C{\varphi }_C{)}^{t-2\Delta t}}{2\Delta t}={\frac{\partial (\rho \varphi )}{\partial t}|}_t-{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^2}{3}-O(\Delta t^3)$$
把上式代入到离散格式中，得
$${\frac{\partial (\rho \varphi )}{\partial t}|}_t+\frac{1}{V_C}L({\varphi }_C^t)={\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^2}{3}+O(\Delta t^3)$$
表明其有二阶精度，虽然有第三项数值色散项，但没有数值扩散效应。

3.7 FV方法的初始条件

有限体积公式的实现是非常直接的，只是初始条件不好处理。如上图所示，第一个时间单元为边界时间单元，所以它没有迎风邻近点。而是把下面的单元面处的值用来表示面上的值（有点拗口，这句话意思是说，本来要用这个下边面的邻近迎风点的值来表征面上值，但是没有这个迎风点，所以只好用面上值来代替了，刚好呢边界面的值是已知的），这就导致梯度实际上只有正确数值的一半。因为计算差分的时候用的是${\varphi }_C^{t_{initial}+\Delta t/2}$和${\varphi }_C^{t_{initial}}$，刚好是相距了$\Delta t/2$，所以如果是把它们除以完整时间步长$\Delta t$的话，是会导致不容忽视的初始误差。

为了说明这个问题，采用一阶隐式Euler格式，写出第一个时间单元的离散方程，有
$$\frac{({\rho }_C{\varphi }_C{)}^{t_{initial}+\Delta t/2}-({\rho }_C{\varphi }_C{)}^{t_{initial}}}{\Delta t}{V}_C+L({\varphi }_C^{t_{initial}+\Delta t/2})=0$$
对于第一个时间单元，迎风插值推出的梯度是由在$t_{initial}+\Delta t/2$和$t_{initial}$时刻的$\rho \varphi$值相减再除以$\Delta t$得到的。然而，对于常规单元，比如第二个单元来说，其梯度实际上是由$t_{initial}+3\Delta t/2$和$t_{initial}+\Delta t/2$时刻的$\rho \varphi$值相减再除以$\Delta t$得到的。两种梯度算法的差别巨大，不管是何种格式，用上式作为梯度的起始计算都会产生较大的初始误差，并严重影响后续时间步的计算。如果采用下图所示的时间单元，则可很好地避免该错误，此时有限差分和有限体积方法都是相同的（相当于错位了半个时间网格）。

采用该方法，跨过时间区间$[t_{initial}+\Delta t/2,t_{initial}+3\Delta t/2]$的第一个单元的迎风值为
$$\begin{array}{rl}({\rho }_C{\varphi }_C{)}^{t_{initial}+3\Delta t/2}& =({\rho }_C{\varphi }_C{)}^{t_{initial}+\Delta t}\\ ({\rho }_C{\varphi }_C{)}^{t_{initial}+\Delta t/2}& =({\rho }_C{\varphi }_C{)}^{t_{initial}}\end{array}$$
代入到半离散格式中，有
$$\frac{({\rho }_C{\varphi }_C{)}^{t_{initial}+\Delta t}-({\rho }_C{\varphi }_C{)}^{t_{initial}}}{\Delta t}{V}_C+L({\varphi }_C^{t_{initial}+\Delta t})=0$$
这与任何内部单元的形式都是一样的。

4 非均匀时间步长

前面讲的都是均匀时间步长，在实际应用中常使用的是变时间步长，即，通过在每个时间步选择不违反CFL条件的最大许用时间步长，来减小计算消耗。

对于一阶格式来说，时间步长变化与否对其离散格式并无影响。对二阶瞬态格式来说就并非如此了，因为在它们的框架中用到了两个时间层次的值。对于Crank-Nicolson瞬态格式的两步法实现而言，只需要把两步中每步的时间步长做更改即可，这当然会影响到离散的精度了，因为其不再是在时间单元的中心了。对其它的二阶方法来说，插值廓线需要加以更改，以便适应不相等的时间步长。下面讲下非均匀瞬态网格上标准CN和SOUE格式对瞬态项的离散方法。尽管有限体积和有限差分方法在均匀网格上的代数方程是一样的，然而在可变时间步长的情况下则未必如此。

4.1 非均匀时间步长的有限差分方法

4.1.1 Crank-Nicolson格式

如上图所示，为了推倒非均匀时间步长时的CN格式，把$t+\Delta t$和$t-\Delta t^{\circ }$时刻的值$\rho \varphi$用$t$时刻的值和偏导做Taylor级数展开，得
$$\begin{array}{rl}(\rho \varphi {)}^{t+\Delta t}& =(\rho \varphi {)}^t+{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}+{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{\Delta t^3}{6}+...\\ (\rho \varphi {)}^{t-\Delta t^{\circ }}& =(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t^{\circ }+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{(\Delta t^{\circ }{)}^2}{2}-{\frac{{\partial }^3(\rho \varphi )}{\partial t^3}|}_t\frac{(\Delta t^{\circ }{)}^3}{6}+...\end{array}$$
把第一个式子乘上$(\Delta t^{\circ }{)}^2$，把第二个式子乘上$\Delta t^2$，并两式相减，得到
$${\frac{\partial (\rho \varphi )}{\partial t}|}_t\approx \frac{(\Delta t^{\circ }{)}^2(\rho \varphi {)}^{t+\Delta t}-[(\Delta t^{\circ }{)}^2-\Delta t^2](\rho \varphi {)}^t-\Delta t^2(\rho \varphi {)}^{t-\Delta t^{\circ }}}{\Delta t(\Delta t^{\circ }{)}^2+\Delta t^{\circ }\Delta t^2}$$
那么半离散方程为
$$\frac{(\Delta t^{\circ }{)}^2(\rho \varphi )-[(\Delta t^{\circ }{)}^2-\Delta t^2](\rho \varphi {)}^{\circ }-\Delta t^2(\rho \varphi {)}^{\circ \circ }}{\Delta t\Delta t^{\circ }(\Delta t+\Delta t^{\circ })}{V}_C+L({\varphi }_C^{\circ })=0$$
空间项展开，最后的代数方程为
$$a_C^{\bullet }{\varphi }_C=b_C-\left(a_C{\varphi }_C^{\circ }+\sum _{F\sim NB(C)}{a}_F{\varphi }_F^{\circ }\right)-a_C^{\circ }{\varphi }_C^{\circ }-a_C^{\circ \circ }{\varphi }_C^{\circ \circ }$$
其中
$$\begin{array}{rl}{a}_C^{\bullet }& =\frac{\Delta t^{\circ }}{\Delta t(\Delta t+\Delta t^{\circ })}{\rho }_C{V}_C\\ a_C^{\circ }& =\frac{\Delta t-\Delta t^{\circ }}{\Delta t\Delta t^{\circ }}{\rho }_C^{\circ }{V}_C\\ a_C^{\circ \circ }& =-\frac{\Delta t}{\Delta t^{\circ }(\Delta t+\Delta t^{\circ })}{\rho }_C^{\circ \circ }{V}_C\end{array}$$
显然，对$\Delta t=\Delta t^{\circ }$，这些系数跟之前推倒的均匀时间步长的系数是吻合的。

4.1.2 Adams-Moulton(或SOUE)格式

如上图所示，Adams-Moulton格式，即SOUE格式，在非均匀时间步长情况下的形式，可把$t-\Delta t$和$t-\Delta t-\Delta t^{\circ }$时刻的值$\rho \varphi$用$t$时刻的值和偏导做Taylor级数展开，得
$$\begin{array}{rl}(\rho \varphi {)}^{t-\Delta t}& =(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_t\Delta t+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{\Delta t^2}{2}+O(\Delta t^3)\\ (\rho \varphi {)}^{t-\Delta t-\Delta t^{\circ }}& =(\rho \varphi {)}^t-{\frac{\partial (\rho \varphi )}{\partial t}|}_t(\Delta t+\Delta t^{\circ })+{\frac{{\partial }^2(\rho \varphi )}{\partial t^2}|}_t\frac{(\Delta t+\Delta t^{\circ }{)}^2}{2}+O(\Delta t^3)\end{array}$$
把第一个式子乘上$(\Delta t+\Delta t^{\circ }{)}^2/\Delta t^2$，并减去第二个式子，得到
$${\frac{\partial (\rho \varphi )}{\partial t}|}_t\approx \frac{1}{\Delta t}\left[\left(1+\frac{\Delta t}{\Delta t+\Delta t^{\circ }}\right)(\rho \varphi {)}^t-\left(1+\frac{\Delta t}{\Delta t^{\circ }}\right)(\rho \varphi {)}^{t-\Delta t}+\left(\frac{\Delta t^2}{\Delta t^{\circ }(\Delta t+\Delta t^{\circ })}\right)(\rho \varphi {)}^{t-\Delta t-\Delta t^{\circ }}\right]$$
那么半离散方程为
$$\left(\frac{1}{\Delta t}+\frac{1}{\Delta t+\Delta t^{\circ }}\right)({\rho }_C{\varphi }_C)V_C-\left(\frac{1}{\Delta t}+\frac{1}{\Delta t^{\circ }}\right)({\rho }_C{\varphi }_C{)}^{\circ }{V}_C+\left(\frac{\Delta t}{\Delta t^{\circ }(\Delta t+\Delta t^{\circ })}\right)({\rho }_C{\varphi }_C{)}^{\circ \circ }{V}_C+L({\varphi }_C)=0$$
空间项展开，最后的代数方程为
$$(a_C^{\bullet }+a_C){\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F=b_C-a_C^{\circ }{\varphi }_C^{\circ }-a_C^{\circ \circ }{\varphi }_C^{\circ \circ }$$
其中
$$\begin{array}{rl}{a}_C^{\bullet }& =\left(\frac{1}{\Delta t}+\frac{1}{\Delta t+\Delta t^{\circ }}\right){\rho }_C{V}_C\\ a_C^{\circ }& =-\left(\frac{1}{\Delta t}+\frac{1}{\Delta t^{\circ }}\right){\rho }_C^{\circ }{V}_C\\ a_C^{\circ \circ }& =\frac{\Delta t}{\Delta t^{\circ }(\Delta t+\Delta t^{\circ })}{\rho }_C^{\circ \circ }{V}_C\end{array}$$
同样地，若$\Delta t=\Delta t^{\circ }$，则这些系数跟之前推导的均匀时间步长的系数是吻合的。

4.2 非均匀时间步长的有限体积方法

根据FVM中的术语，时间单元的长度定义为$\Delta t$，而两个相继时间单元的形心间的距离定义为为$\delta t$。对于均匀网格步长，它们俩是相等的，即$\Delta t=\delta t$，而且有限差分方法中两个相续计算流场间的时间差为$\Delta t$，等于有限体积法中的时间差$\delta t$。对于非均匀时间步长，有限差分方法仍旧是$\Delta t$，而有限体积方法则变成$\delta t=(\Delta t+\Delta t^{\circ })/2$，这导致两者的公式是不同的。

和在有限差分方法中类似，在非均匀时间步长中，当前和旧时刻的值影响着格式的插值廓线和有限体积离散。这跟在结构化的非均匀网格上把对流项写成廓线形式非常相像。将用CN和SOUE格式来阐明该流程，把该思想推广到其它廓线上也是非常方便的。

4.2.1 Crank-Nicolson格式

如上图所示，CN格式可以通过把面上的$\rho \varphi$值用该面两侧的主节点$\rho \varphi$值的平均来获得，即
$$\begin{array}{rl}({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}& =\frac{\Delta t^{\circ }}{\Delta t+\Delta t^{\circ }}({\rho }_C{\varphi }_C{)}^t+\frac{\Delta t}{\Delta t+\Delta t^{\circ }}({\rho }_C{\varphi }_C{)}^{t-(\Delta t^{\circ }+\Delta t)/2}\\ ({\rho }_C{\varphi }_C{)}^{t-\Delta t/2-\Delta t^{\circ }}& =\frac{\Delta t^{\circ \circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}({\rho }_C{\varphi }_C{)}^{t-(\Delta t^{\circ }+\Delta t)/2}+\frac{\Delta t^{\circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}({\rho }_C{\varphi }_C{)}^{t-\Delta t^{\circ }-(\Delta t+\Delta t^{\circ \circ })/2}\end{array}$$
那么半离散方程为
$$\begin{array}{rl}\frac{\Delta t^{\circ }}{\Delta t+\Delta t^{\circ }}\frac{V_C}{\Delta t}({\rho }_C{\varphi }_C)& +\left(\frac{\Delta t}{\Delta t+\Delta t^{\circ }}-\frac{\Delta t^{\circ \circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\right)\frac{V_C}{\Delta t}({\rho }_C{\varphi }_C{)}^{\circ }\\ & -\frac{\Delta t^{\circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\frac{V_C}{\Delta t}({\rho }_C{\varphi }_C{)}^{\circ \circ }+L({\varphi }_C^{\circ })=0\end{array}$$
推得CN格式在非均匀时间步长时的线性化系数为
$$\begin{array}{rl}FluxC& =\frac{\Delta t^{\circ }}{\Delta t+\Delta t^{\circ }}\frac{{\rho }_C{V}_C}{\Delta t}\\ Flux{C}^{\circ }& =\left(\frac{\Delta t}{\Delta t+\Delta t^{\circ }}-\frac{\Delta t^{\circ \circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\right)\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\\ FluxV& =-\frac{\Delta t^{\circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\frac{{\rho }_C^{\circ \circ }{V}_C{\varphi }_C^{\circ \circ }}{\Delta t}\end{array}$$
和在固定时间步长中一样，该显式方法需要存储前两个时间步的值。此外，令$\Delta t=\Delta t^{\circ }=\Delta t^{\circ \circ }$即可将其恢复到均匀时间步长的公式。

4.2.2 Adams-Moulton(或SOUE)格式

如上图所示，采用二阶迎风插值廓线，把在$t+\Delta t/2$和$t-\Delta t/2$面上的$\rho \varphi$值可表示为
$$\begin{array}{rl}({\rho }_C{\varphi }_C{)}^{t+\Delta t/2}& =({\rho }_C{\varphi }_C{)}^t+\frac{\Delta t}{\Delta t+\Delta t^{\circ }}\left[({\rho }_C{\varphi }_C{)}^t-({\rho }_C{\varphi }_C{)}^{t-(\Delta t^{\circ }+\Delta t)/2}\right]\\ & \\ ({\rho }_C{\varphi }_C{)}^{t-\Delta t/2}& =({\rho }_C{\varphi }_C{)}^{t-(\Delta t^{\circ }+\Delta t)/2}+\frac{\Delta t^{\circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\left[({\rho }_C{\varphi }_C{)}^{t-(\Delta t^{\circ }+\Delta t)/2}-({\rho }_C{\varphi }_C{)}^{t-\Delta t^{\circ }-(\Delta t+\Delta t^{\circ \circ })/2}\right]\end{array}$$
那么半离散方程为
$$\begin{array}{rl}\left(1+\frac{\Delta t}{\Delta t+\Delta t^{\circ }}\right)\frac{V_C}{\Delta t}({\rho }_C{\varphi }_C)& -\left(1+\frac{\Delta t}{\Delta t+\Delta t^{\circ }}+\frac{\Delta t^{\circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\right)\frac{V_C}{\Delta t}({\rho }_C{\varphi }_C{)}^{\circ }\\ & +\frac{\Delta t^{\circ }}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\frac{V_C}{\Delta t}({\rho }_C{\varphi }_C{)}^{\circ \circ }+L({\varphi }_C)=0\end{array}$$
推得SOUE格式在非均匀时间步长时的线性化系数为
$$\begin{array}{rl}FluxC& =\left(\frac{1}{\Delta t}+\frac{1}{\Delta t+\Delta t^{\circ }}\right){\rho }_C{V}_C\\ Flux{C}^{\circ }& =-\left(\frac{1}{\Delta t}+\frac{1}{\Delta t+\Delta t^{\circ }}+\frac{\Delta t^{\circ }/\Delta t}{\Delta t^{\circ }+\Delta t^{\circ \circ }}\right){\rho }_C^{\circ }{V}_C\\ FluxV& =\frac{\Delta t^{\circ }/\Delta t}{\Delta t^{\circ }+\Delta t^{\circ \circ }}{\rho }_C^{\circ \circ }{V}_C{\varphi }_C^{\circ \circ }\end{array}$$
与固定时间步情况类似，该方法是隐式的，其需要求解方程组系统来获得在每个时间步上的$\varphi$值。令$\Delta t=\Delta t^{\circ }=\Delta t^{\circ \circ }$即可将其恢复到均匀时间步长的公式。

5 代码讲解

5.1 uFVM

在uFVM中，瞬态项的离散时通过有限体积方法在隐式框架下完成的。用一阶后向Euler格式所得到的的瞬态通量的组装过程如下：

%
theDensityField = cfdGetMeshField(['Density' theFluidTag]);
density = theDensityField.phi(iElements);
density_old = theDensityField.phi_old(iElements);

volumes = [theMesh.elements(iElements).volume]';

theFluxes.FLUXCE(iElements) = volumes .* density / dt;
theFluxes.FLUXCEOLD(iElements) = - volumes .* density_old / dt;
theFluxes.FLUXTE(iElements) = theFluxes.FLUXCE .* phi;
theFluxes.FLUXTEOLD(iElements) = theFluxes.FLUXCEOLD .* phi_old ;

即，一阶隐式Euler格式的线性化后的相关系数的计算和组装工作
$$\begin{array}{rl}FluxC& =\frac{{\rho }_C{V}_C}{\Delta t}\\ Flux{C}^{\circ }& =-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}\\ FluxV& =0\end{array}$$

5.2 OpenFOAM

to be continued…