FVM in CFD 学习笔记_第13章_时域离散:瞬态项

2020年06月10日 阅读数:236
这篇文章主要向大家介绍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 Termhtml

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

1 引言

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

通常而言,某个变量 ϕ \phi 的瞬态特性的表达式,或其时间演化,是由下述形式的方程控制的
( ρ ϕ ) t + L ( ϕ ) = 0 \frac{\partial (\rho \phi)}{\partial t}+L(\phi)=0
其中 L ( ϕ ) L(\phi) 为空间算子,其包含了全部的非瞬态项(扩散项、对流项、源项,等), ( ρ ϕ ) / t {\partial (\rho \phi)}/{\partial t} 为瞬态算子,二者以下图所示。架构

在这里插入图片描述

在这里插入图片描述

在单元 C C 上对上式作积分,得
V C ( ρ ϕ ) t d V + V C L ( ϕ ) d V = 0 \int_{V_C}\frac{\partial (\rho \phi)}{\partial t}dV+\int_{V_C}L(\phi)dV=0
空间离散后,变为
( ρ C ϕ C ) t V C + L ( ϕ C t ) = 0 \frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0
其中 V C V_C 为离散单元的体积, L ( ϕ C t ) L(\phi_C^t) 为空间算子在某参考时刻 t t 的离散形式,可写成以下代数方程
L ( ϕ C t ) = a C ϕ C t + F N B ( C ) a F ϕ F t b C L(\phi_C^t)=a_C\phi_C^t+\sum_{F\sim NB(C)}a_F\phi_F^t-b_C
在上上式中,当 t t\rightarrow\infin 时其变为稳态离散方程。经过时间推动也能达到稳态,即, ϕ C t + Δ t = ϕ C t \phi_C^{t+\Delta t}=\phi_C^t 。这保证了瞬态问题在到达稳态时所得到的解,和直接由稳态问题求得的解,是同样的。app

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

2 有限差分方法

在这里插入图片描述

如图,在瞬态空间中的网格是结构化的,因此对瞬态项使用有限差分方法是很是稀松日常的操做。在该方法中,空间算子 L ( ϕ ) L(\phi) 是在时刻 t t 离散的,而瞬态偏导数则可用时刻 t t 的Taylor展开而推导出不一样的瞬态格式,其中的一些格式呈现以下。svg

2.1 前向Euler格式(Forward Euler Scheme)

为了估算瞬态项,所推导量的Taylor展开须要指定时间方向。这里,展开是针对时间向前进行的。对某函数 T T ,其在 t + Δ t t+\Delta t 时刻的值,可用其在时刻 t t 的值和导数值作Taylor级数展开,以下
T ( t + Δ t ) = T ( t ) + T ( t ) t Δ t + 2 T ( t ) t 2 Δ 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!}+...
截去二阶以上项,获得
T ( t ) t = T ( t + Δ t ) T ( t ) Δ t + O ( Δ t ) \frac{\partial T(t)}{\partial t}=\frac{T(t+\Delta t)-T(t)}{\Delta t}+O(\Delta t)
只有一阶精度,把上式中的 T T 替换成 ( ρ ϕ ) (\rho\phi) ,并代入到 ( ρ C ϕ C ) t V C + L ( ϕ C t ) = 0 \frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0 ,离散方程变为
( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t V C + L ( ϕ 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 + Δ t t+\Delta t 时刻的 ( ρ C ϕ C ) (\rho_C \phi_C) 是不须要求解方程组系统的。由于全部的空间项都是在旧时刻 t t 上衡量的,因此 t + Δ t t+\Delta t 时刻的 ϕ C \phi_C 可直接由上个时间步的值显式算出。这样的格式属于显式瞬态格式。显式瞬态格式的主要特色是经过在时域推动即可得到解,而无需在每一个时间层求解方程组系统。这样一来,求解效率很是高,且对计算网格的并行处理十分方便。然而鲜有商业软件采用该方法,由于其重大缺陷就是时间步长 Δ t \Delta t 极为受限,这将在下一小节详细展开。学习

把空间算子的离散代数方程代入到上式,可得完整的代数方程为
a C t + Δ t ϕ C t + Δ t + a C t ϕ C t = b C ( a C ϕ C t + F N B ( C ) a F ϕ F 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)
其中
a C t + Δ t = ρ C t + Δ t V C Δ t a C t = ρ C t V C Δ t \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 + Δ t a_C^{t+\Delta t} a C t a_C^t 为对角系数,源自于瞬态项的离散, ϕ C t + Δ t \phi_C^{t+\Delta t} ϕ C t \phi_C^t 为时间层 t + Δ t t+\Delta t t t 上的值,而 a C a_C a F a_F b C b_C 则是空间离散后的系数。ui

为了简化标记,在本章中,前一个时间步的变量值用上标 ^\circ 标记,前两个时间步的变量值用上标 ^{\circ\circ} 标记,如果没有上标,则表示当前时间步的变量值,除了在非定常项中乘到 ϕ C \phi_C 上的系数是用上标 ^\bullet 来标记的。基于这些标记,上式写做
a C ϕ C + a C ϕ C = b C ( a C ϕ C + F N B ( C ) a F ϕ F ) 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)
其中
a C = ρ C V C Δ t a C = ρ C V C Δ t \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}
改写为以下形式
ϕ C = b C [ ( a C + a C ) ϕ C + F N B ( C ) a F ϕ F ] a C \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条件能够被简单解释为,把系数所要知足的基本规则之一,即,反符号准则,扩展成把瞬态系数也包含在内。这样,正如 ϕ F \phi_F ϕ C \phi_C 的“空间”邻居, ϕ C \phi_C^\circ 也是 ϕ C \phi_C 的“时间”邻居,那么反符号准则对他俩应该平等适用。注意对角线系数如今是 a C a_C^\bullet ,其“时间”邻居的系数是 ( a C + a C ) (a_C+a_C^\circ) ,反符号准则变为
a C + a C 0 a_C+a_C^\circ\le0

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

在这里插入图片描述

如上图,对于一维纯对流问题,流动方向向右,使用迎风格式把变量插值到单元面上去,则单元 C C 的离散方程中的系数 a C a_C a C a_C^\circ
a C = m ˙ e = ρ C u C Δ y C a C = ρ C V C Δ t = ρ C Δ x C Δ y C Δ t \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 0 ρ C u C Δ y C ρ C Δ x C Δ y C Δ t 0 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

Δ t Δ x C u C \Delta t\le\frac{\Delta x_C}{u_C^\circ}
对于对流控制的流动,定义CFL数为
C F L c o n v = v C Δ t Δ x C CFL^{conv}=\frac{|\bold v_C^\circ|\Delta t}{\Delta x_C}
这意味着为了数值稳定性,CFL数应该知足
C F L c o n v 1 CFL^{conv}\le1

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

在这里插入图片描述

对于纯扩散问题,CFL数的表达式是不一样的。所以,分析如上图所示的一维纯扩散问题。

使用线性插值廓线,单元 C C 的离散方程中的系数 a C a_C a C a_C^\circ
a C = Γ e ϕ Δ y C δ x e + Γ w ϕ Δ y C δ x w a C = ρ C V C Δ t = ρ C Δ x C Δ y C Δ t \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 0 Γ e ϕ Δ y C δ x e + Γ w ϕ Δ y C δ x w ρ C Δ x C Δ y C Δ t 0 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

Δ t ρ C Δ x C Γ e ϕ δ x e + Γ w ϕ δ x w \Delta t\le\frac{\rho_C^\circ\Delta x_C}{\dfrac{\Gamma_e^\phi}{\delta x_e} + \dfrac{\Gamma_w^\phi}{\delta x_w}}
若网格均匀,且扩散系数为常数,则上式变为
Δ t ρ C ( Δ x C ) 2 2 Γ C ϕ \Delta t\le\frac{\rho_C^\circ(\Delta x_C)^2}{2\Gamma_C^\phi}
对于扩散控制的问题,定义CFL数为
C F L d i f f = Γ C ϕ Δ t ρ C ( Δ x C ) 2 CFL^{diff}=\frac{\Gamma_C^\phi\Delta t}{\rho_C^\circ(\Delta x_C)^2}
这意味着为了数值稳定性,CFL数应该知足
C F L d i f f 1 2 CFL^{diff}\le \frac{1}{2}

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

在这里插入图片描述

对于多维非定常对流扩散问题,基于上一章(第12章)的推导,系数 a C a_C a C a_C^\circ
a C = f n b ( C ) ( Γ f ϕ E f d C F + m ˙ f , 0 ) a C = ρ C V C Δ t \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 0 f n b ( C ) ( Γ f ϕ E f d C F + m ˙ f , 0 ) ρ C V C Δ t 0 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
即以下对时间步长的限制关系
Δ t ρ C V C f n b ( C ) ( Γ f ϕ E f d C F + m ˙ f , 0 ) \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) }
上式一般是显式瞬态格式的稳定性要求,实际上,以前对一维区域的纯对流和纯扩散的条件,能够视做是上式的特例。好比一维扩散问题,有均匀网格单元尺度 Δ x \Delta x ,常密度 ρ \rho ,均匀扩散系数 Γ ϕ \Gamma^\phi ,上式变为
Δ t ρ C V C f n b ( C ) ( Γ f ϕ E f d C F + m ˙ f , 0 ) = ρ C Δ x C Δ y C ( Γ e ϕ + Γ w ϕ ) Δ y C Δ x C = ρ C ( Δ x C ) 2 2 Γ C ϕ \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}
而对于纯对流的一维问题,使用迎风格式离散对流项,流体从左向右流动,则其变为了
Δ t ρ C V C f n b ( C ) ( Γ f ϕ E f d C F + m ˙ f , 0 ) = ρ C Δ x C Δ y C ρ C u C Δ y C = Δ x C u C \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 Δ t t-\Delta t 的函数 T T 值用在时刻 t t 的函数 T T 值作Taylor级数展开,有
T ( t Δ t ) = T ( t ) T ( t ) t Δ t + 2 T ( t ) t 2 Δ 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!}+...
整理后,有
T ( t ) t = T ( t ) T ( t Δ t ) Δ t + 2 T ( t ) t 2 Δ t 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 T 替换成 ( ρ ϕ ) (\rho\phi) ,并代入到 ( ρ C ϕ C ) t V C + L ( ϕ C t ) = 0 \frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0 ,离散方程变为
( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t Δ t V C + L ( ϕ 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 + a C ) ϕ C + F N B ( C ) a F ϕ F = b C a C ϕ C (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 = ρ C V C Δ t a C = ρ C V C Δ t \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 a C a_C^\circ 的符号是相反的(知足反号准则,即对角系数和邻居系数符号是相反的),这样即可保证 ϕ C \phi_C 是由其 当前时间步的空间邻近点的值 和 前一时间步 t Δ t t-\Delta t 的时间邻近点的值 所限定的。这也意味着该格式老是稳定的,而无论用的是何种时间步长,从而能够用很大的时间步长向前快速推动。纵然如此,其并不是是理想格式,由于其精度较低,这种格式得到的解精确性差,除非用很小的时间步长才好,这让其使用起来颇为尴尬。若采用大时间步长来提升计算效率,则得到的解精度不好。若使用小时间步长来得到精确解,则计算效率会很是低下。

2.4 Crank-Nicolson格式

在Crank-Nicolson格式中,为了得到瞬态项更加精确的近似,采用了在 t Δ t t-\Delta t 时刻的函数 T T 值和在 t + Δ t t+\Delta t 时刻的函数 T T 值的展开形式,一样是用 t t 时刻的量来展开的,即
T ( t + Δ t ) = T ( t ) + T ( t ) t Δ t + 2 T ( t ) t 2 Δ t 2 2 ! + 3 T ( t ) t 3 Δ t 3 3 ! + . . . T ( t Δ t ) = T ( t ) T ( t ) t Δ t + 2 T ( t ) t 2 Δ t 2 2 ! 3 T ( t ) t 3 Δ t 3 3 ! + . . . \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}
两式相减,得
T ( t ) t = T ( t + Δ t ) T ( t Δ t ) 2 Δ t + O ( Δ t 2 ) \frac{\partial T(t)}{\partial t}=\frac{T(t+\Delta t)-T(t-\Delta t)}{2\Delta t}+O(\Delta t^2)
注意,偏导的精度的阶数为 O ( Δ t 2 ) O(\Delta t^2) ,由于二阶偏导被彻底消掉了。

把上式中的 T T 替换成 ( ρ ϕ ) (\rho\phi) ,并代入到 ( ρ C ϕ C ) t V C + L ( ϕ C t ) = 0 \frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0 ,离散方程变为
( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t 2 Δ t V C + L ( ϕ 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 ϕ C = b C ( a C ϕ C + F N B ( C ) a F ϕ F ) a C ϕ C 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}
其中
a C = ρ C V C 2 Δ t a C = ρ C V C 2 Δ t \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}
其架构以下图所示,显然,这个格式是显示算法,由于 ( ρ ϕ ) t + Δ t (\rho\phi)^{t+\Delta t} 的获取只用到了旧时刻的值。然而如今须要存储两个旧时刻的值,但空间离散算子仍是只须要其中一个旧时刻的值就行了。
在这里插入图片描述

对于CN格式的稳定性分析,能够把其最初的方程作些许修改来进行,使用以下近似:
ϕ ϕ + ϕ 2 \phi^\circ\approx\frac{\phi+\phi^{\circ\circ}}{2}
代数方程变为
a C ϕ C + 0.5 ( a C ϕ C + F N B ( C ) a F ϕ F ) = b C 0.5 ( ( a C + 2 a C ) ϕ C + F N B ( C ) a F ϕ F ) 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 + 2 a C 0 a_C+2a_C^{\circ\circ} \le 0
对于一维瞬态对流问题,有
Δ t 2 ρ C V C m ˙ e = 2 ρ C Δ x C Δ y C ρ C u C Δ y C 2 Δ x C v e \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数,上式变为
C F L c o n v 2 CFL^{conv} \le 2
这个CFL数的最大限制可要比正向Eluer格式的大了,很是惬意,然而精度的提高则是更有意义,由于其让精确解的获取无需经过诉诸于更小的时间步长来实现,尤为是其把二阶导数都从偏差中消掉了。随后章节中将展开更为详细的精确性分析。

2.5 实现细节

CN格式也能够经过把前向和后向Euler格式相加而推出,以下
Forward Euler ( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t V C + L ( ϕ C t ) = 0 Backward Euler ( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t Δ t V C + L ( ϕ C t ) = 0 \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}
二者相加
( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t V C + ( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t Δ t V C + 2 L ( ϕ C t ) = 0 ( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t 2 Δ t V C + L ( ϕ C t ) = 0 \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格式,用于隐式找到 ( ρ ϕ ) t (\rho\phi)^t
( ρ C ϕ C ) t + Δ t V C L ( ϕ C t ) = ( ρ C ϕ C ) t Δ t (\rho_C \phi_C)^{t}+\frac{\Delta t}{V_C}L(\phi_C^t)=(\rho_C \phi_C)^{t-\Delta t}
而后在第2步中,显式算得 t + Δ t t+\Delta t 时刻的CN值
( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t V C = L ( ϕ C t ) = ( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t Δ t V C ( ρ C ϕ C ) t + Δ t = 2 ( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t \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}
在该推导过程当中,假设瞬态时间步长 Δ t \Delta t 分为了两个等长的局部时间步长 Δ t l o c a l \Delta t_{local} ,且 Δ t l o c a l \Delta t_{local} 为设置时间步长 Δ t \Delta t 的一半。

值得注意的是,尽管CN格式是二阶精度,可是它仍旧是显式格式,且受限于CFL条件,如前所示。

2.6 Adams-Moulton格式

二阶Adams-Moulton格式的发展须要把 T T 在时刻 t Δ t t-\Delta t t 2 Δ t t-2\Delta t 的值使用Taylor级数展开成 t t 时刻的形式,即
T ( t 2 Δ t ) = T ( t ) T ( t ) t 2 Δ t + 2 T ( t ) t 2 4 Δ t 2 2 ! + . . . T ( t Δ t ) = T ( t ) T ( t ) t Δ t + 2 T ( t ) t 2 Δ t 2 2 ! + . . . \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}
从上面两个式子出发,想办法把二阶偏导项消掉,得
T ( t ) t = 3 T ( t ) 4 T ( t Δ t ) + T ( t 2 Δ t ) 2 Δ t \frac{\partial T(t)}{\partial t}=\frac{3T(t)-4T(t-\Delta t)+T(t-2\Delta t)}{2\Delta t}
把上式中的 T T 替换成 ( ρ ϕ ) (\rho\phi) ,并代入到 ( ρ C ϕ C ) t V C + L ( ϕ C t ) = 0 \frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0 ,离散方程变为
3 ( ρ C ϕ C ) t 4 ( ρ C ϕ C ) t Δ t + ( ρ C ϕ C ) t 2 Δ t 2 Δ t V C + L ( ϕ 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 + a C ) ϕ C + F N B ( C ) a F ϕ F = b C a C ϕ C a C ϕ C (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}
其中
a C = 3 ρ C V C 2 Δ t a C = 2 ρ C V C Δ t a C = ρ C V C 2 Δ t \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 a_C^{\circ\circ} 系数是正值,因此 ϕ C \phi_C^{\circ\circ} 的增长会致使 ϕ C \phi_C 的减少。因为 a C a_C^\circ 系数是负的,因此若 a C a_C^\circ 较大,则可略微削减该不利影响。如此看来,尽管该格式是稳定的,然而其是无界的,在某些状况下会出现非物理振荡。

3 有限体积方法

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

在这里插入图片描述

对式 ( ρ C ϕ C ) t V C + L ( ϕ C t ) = 0 \frac{\partial (\rho_C \phi_C)}{\partial t}V_C+L(\phi_C^t)=0 在时间区间 [ t Δ t / 2 , t + Δ t / 2 ] [t-\Delta t/2, t+\Delta t/2] 作积分,即
t Δ t / 2 t + Δ t / 2 ( ρ C ϕ C ) t V C d t + t Δ t / 2 t + Δ t / 2 L ( ϕ C t ) d t = 0 \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 V_C 做为常数处理,则第1项变为面通量差分形式,第2项用积分中值定理来作体积分,得
V C ( ρ C ϕ C ) t + Δ t / 2 V C ( ρ C ϕ C ) t Δ t / 2 + L ( ϕ C t ) Δ t = 0 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
该半离散瞬态方程可写成更加标准的形式,两边都除以 Δ t \Delta t ,得
( ρ C ϕ C ) t + Δ t / 2 ( ρ C ϕ C ) t Δ t / 2 Δ t V C + L ( ϕ C t ) = 0 \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 + Δ t / 2 t+\Delta t/2 t Δ t / 2 t-\Delta t/2 处的面值用 t t t Δ t t-\Delta t 等处的单元值来表示。该廓线的选择可从对流项离散中汲取灵感。选择不一样的廓线,毫无疑问会影响到方法的精确性和稳定性。基于此,须要说一下,空间算子的积分是时域二阶精度的,可是空间算子自己的精度(空间精度)则是由其离散时所采用的格式所决定的。

无论是使用的哪一种廓线,通量均可用新值和旧值线性化为
F l u x T = F l u x C ϕ C + F l u x C ϕ C + F l u x V FluxT=FluxC\phi_C+FluxC^\circ\phi_C^\circ+FluxV
其中上标 ^\circ 仍旧表明旧值。完成该线性化处理后,代数方程的相关系数更替为下式,便可
a C a C + F l u x C b C b C F l u x C ϕ C F l u x V \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 值,等于迎风单元形心值,即
( ρ C ϕ C ) t + Δ t / 2 = ( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t / 2 = ( ρ C ϕ C ) t Δ t (\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}
那么半离散方程为
( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t Δ t V C + L ( ϕ 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
即一阶隐式Euler格式,线性化后的相关系数为
F l u x C = ρ C V C Δ t F l u x C = ρ C V C Δ t F l u x V = 0 \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 t 的Taylor级数展开,将其恢复到最初控制方程的形式,来肯定数值扩散的值,即
( ρ ϕ ) t Δ t = ( ρ ϕ ) t ( ρ ϕ ) t t Δ t + 2 ( ρ ϕ ) t 2 t Δ t 2 2 + O ( Δ t 3 ) (\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)
调整为
( ρ ϕ ) t ( ρ ϕ ) t Δ t Δ t = ( ρ ϕ ) t t ( Δ t 2 ) 2 ( ρ ϕ ) t 2 t O ( Δ t 2 ) \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)
把上式代入到离散格式中,得
( ρ ϕ ) t t + 1 V C L ( ϕ C t ) = ( Δ t 2 ) 2 ( ρ ϕ ) t 2 t Numerical diffusion term + O ( Δ 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 值,等于背风单元形心值,即
( ρ C ϕ C ) t + Δ t / 2 = ( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t / 2 = ( ρ C ϕ C ) t (\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}
那么半离散方程为
( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t V C + L ( ϕ 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
即一阶显式Euler格式,线性化后的相关系数为
F l u x C = ρ C V C Δ t F l u x C = ρ C V C Δ t F l u x V = 0 \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 + Δ t t+\Delta t ,空间算子则是在旧时刻 t t 上衡量的,这样一来,右端项就彻底已知了,能够直接用来获得 t + Δ t t+\Delta t 时刻的 ρ ϕ \rho\phi ,而不须要求解线性代数方程组。这是显式格式。

3.3.1 数值反扩散

基于 t t 时刻作Taylor展开
( ρ ϕ ) t + Δ t = ( ρ ϕ ) t + ( ρ ϕ ) t t Δ t + 2 ( ρ ϕ ) t 2 t Δ t 2 2 + O ( Δ t 3 ) (\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)
调整为
( ρ ϕ ) t + Δ t ( ρ ϕ ) t Δ t = ( ρ ϕ ) t t + ( Δ t 2 ) 2 ( ρ ϕ ) t 2 t + O ( Δ t 2 ) \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)
把上式代入到离散格式中,得
( ρ ϕ ) t t + 1 V C L ( ϕ C t ) = ( Δ t 2 ) 2 ( ρ ϕ ) t 2 t Numerical anti-diffusion term + O ( Δ 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的时候,对流项的数值扩散和 C F L c o n v = 1 CFL^{conv}=1 的显式Euler格式的数值反扩散,二者幅值相等且符号相反,所以二者相互抵消,产生了近乎精确的解。尽管如此,确保 C F L c o n v = 1 CFL^{conv}=1 是不切实际的,并且简单的一维网格也并不是是实际问题的再现。

与反扩散效应相关的问题是数值不稳定性,其随着 Δ t \Delta t 的增长而增长,所以须要对时间步长施加较强的限制,能够经过负邻近系数准则来评估最大稳定时间步。

3.4 二阶瞬态Euler格式

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

3.5 Crank-Nicolson(中心差分廓线)

在这里插入图片描述

采用在迎风和背风节点之间的线性插值来计算 ρ ϕ \rho\phi 值,即可得到上图所示的Crank-Nicolson格式。

对于均匀时间步,其数学表达式为
( ρ C ϕ C ) t + Δ t / 2 = 1 2 ( ρ C ϕ C ) t + Δ t + 1 2 ( ρ C ϕ C ) t ( ρ C ϕ C ) t Δ t / 2 = 1 2 ( ρ C ϕ C ) t + 1 2 ( ρ C ϕ C ) t Δ t \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}
那么半离散方程为
( ρ C ϕ C ) t + Δ t ( ρ C ϕ C ) t Δ t 2 Δ t V C + L ( ϕ 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
即一阶显式Euler格式,线性化后的相关系数为
F l u x C = ρ C V C 2 Δ t F l u x C = 0 F l u x V = ρ C V C 2 Δ t ϕ C \begin{aligned} FluxC &= \frac{\rho_C V_C}{2\Delta t} \\ FluxC^\circ &= 0 \\ FluxV &= -\frac{\rho_C^{\circ\circ} V_C}{2\Delta t}\phi_C^{\ci