FVM in CFD 学习笔记_第15章_流动计算:不可压缩流动_2_同位网格上的SIMPLE算法

2020年06月10日 阅读数:317
这篇文章主要向大家介绍FVM in CFD 学习笔记_第15章_流动计算:不可压缩流动_2_同位网格上的SIMPLE算法,主要内容包括基础应用、实用技巧、原理机制等方面,希望对大家有所帮助。

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 15 Fluid Flow Computation: Incompressible Flowshtml

前面几章讲解的关于变量 ϕ \phi 的通常输运方程的离散和求解流程,均是创建在事先已知速度场的前提下。可是通常状况下,速度场是未知的,且须要经过求解Navier-Stokes方程组来获取。对于不可压缩流动而言,该项工做尤其复杂,由于压力和速度是强烈耦合在一块儿的,并且压力并不以主变量的形式出如今动量或是连续方程中。本章的重点在于展现解决上述两个问题的方法,以及不可压缩流动问题的流场计算方法。先讲解一维交错网格,而后是一维同位网格,最后是同位三维非结构网格。除了阐明SIMPLE、SIMPLEC、PRIME和PISO算法的前因后果,还清晰阐明了Rhie-Chow插值方法,以及将其扩展到瞬态、松弛和体积力项的方法。最后,展示了一些常见边界条件的添加细节。web

因为本章内容繁杂,篇幅较长,故分红了四部分来说解,各部分主要内容分别为:交错网格、同位网格、边界条件、SIMPLE家族算法。算法

这里是第二部分,主要讲解交错网格的缺陷,以及如何不用交错网格,而直接在原来的同位网格(即最初的网格,且是三维非结构的复杂网格)上开展SIMPLE算法。架构

3 交错网格的缺陷

交错网格的使用对于SIMPLE算法的发展来讲相当重要,然而采用交错网格架构也有其缺点。如前所述,在二维和三维状况下,分别须要三套和四套交错网格系统,来存放速度的两到三个不一样份量和压力值。app

这样一来,就须要对每一个速度份量存储一套网格系统,对压力和其余变量再存储一套网格系统,内存开销是蛮大的。除此以外,对于非笛卡尔网格,交错网格系统会产生不少问题,更有甚者,对于非结构网格而言,交错网格系统基本上是没法构造出来的。框架

在这里插入图片描述

在曲线网格中,使用笛卡尔速度份量会致使的问题是,当一个或多个面和交错速度份量的方向平齐而非垂直的时候,好比上图中的最上端的两个单元中,因为每一个面上的交错速度与面平齐,算得的质量流量都几乎是零,这样形成其错误地自动知足了连续方程,显然这是不符合物理意义的,会致使SIMPLE计算过程的失效。
在这里插入图片描述ide

所以,这种状况下较好的替代方法是使用协变的(covariant)或抗变(contra-variant)速度份量,如上图a和b所示。svg

使用抗变速度份量的交错网格例子以下图。函数

不幸的是,在曲线坐标系上离散动量方程时,复杂性大为提高,表如今对扩散项的处理上,由于方程中有非守恒项。
在这里插入图片描述学习

另外一种处理方法是把在全部方向的笛卡尔速度份量所有作交错处理,这样在每一个面上都含有全部的速度份量了,以下图。这样,待求解动量方程的数目将会是两倍(二维问题)或三倍(三维问题)。

在这里插入图片描述

对非结构网格状况来讲,问题变得更加复杂了,由于此时并无一个明显的交错方向,那么应用交错概念的惟一方法是经过改变压力和速度份量单元的大小,或者经过诉诸于在每一个面上交错全部的速度份量,这无疑会急剧增长待求解变量的数目。最终,所存储的几何构型将远大于两倍,由于须要全新的非结构网格来存储速度份量。

结果是我们还得回到最初的单元中心的同位网格系统上,以下图所示,全部的变量都存在相同的位置处(单元形心),使用这种同位网格系统才是最佳方法(兜兜转转饶了一大圈又转回来了哈)。值得注意的是,虽然速度份量是存储在单元形心处的,与压力和其余变量存储位置同样,然而质量流量(是个标量)在同位网格中是存储在面单元上的。质量流量实际上能够视做抗变(contra-variant)份量,只是此时质量流量的计算一般是用离散动量方程来插值计算的(再也不像以前那样子用界面速度,而界面速度又用两侧单元值插值,从而出现了连续方程中用到的是交替单元的速度的状况,没法感知棋盘型的非均匀流场),这即是著名的Rhie-Chow插值,其正是下面章节的主题。
在这里插入图片描述

4 Rhie-Chow插值

前面展现的原始同位网格上的缺陷(即没法感知棋盘型压力场和速度场),究其缘由,是因为采用了线性插值来计算单元面上的速度。该插值致使了在单元层次的压力和速度值的解耦(即,用到的是单元 C C 两侧单元 E E W W 的值来计算压力梯度和连续方程中的通量,而并非用的紧邻单元值),从而产生了棋盘型问题(对棋盘型的非均匀场将被错误地感知成均匀场)。1983年,Rhie和Chow发表了一种插值策略,其容许SIMPLE算法在同位网格上实现。以下图所示,在他们的方法中有个耗散项,其表明着两种预估单元面压力梯度的差值(即,后面要讲到的 D f u [ ( p / x ) f ( p / x ) f ] \overline{D_f^u}\left[({\partial p}/{\partial x})_f-\overline{({\partial p}/{\partial x})_f}\right] 项),该耗散项被添加进单元面速度 u f u_f 的线性插值中。注意,这两个压力梯度的预估值 ( p / x ) f \left({\partial p}/{\partial x}\right)_f ( p / x ) f \overline{\left({\partial p}/{\partial x}\right)_f} 是基于不一样的网格框架的,后者是线性插值(用到四个单元形心点),前者则是Rhie-Chow插值(用到面邻近连续俩单元形心点)。

在这里插入图片描述

该流程等效于在单元面处构建一个伪动量方程,其系数是从该面两侧单元形心处的动量方程中的系数线性插值而来,而压力梯度的计算则使用最小网格框架(即 C C F F 点的压力值来计算)。如此一来,Rhie-Chow插值轻易地模仿出了交错网格配置下的压力-速度耦合最小框架。

(说到这,连我都被绕晕了,不要紧,接下来看看如何处理就会豁然开朗了)

由单元 C C F F 的离散 x x 动量方程开始,即
u C + H C [ u ] = B C u D C u ( p x ) C u F + H F [ u ] = B F u D F u ( p x ) F u_C+H_C[u]=B_C^u-D_C^u\left(\frac{\partial p}{\partial x}\right)_C \\ u_F+H_F[u]=B_F^u-D_F^u\left(\frac{\partial p}{\partial x}\right)_F
u f u_f 速度方程与上式相似,只是压力梯度是与局部邻近压力值相关的( f f 处压力梯度是用 C C F F 压力值计算的,并不是间隔了一个单元的交替值),如上图所示, u f u_f 方程的形式为
u f + H f [ u ] = B f u D f u ( p x ) f u_f+H_f[u]=B_f^u-D_f^u\left(\frac{\partial p}{\partial x}\right)_f
既然在同位网格上,上式中的系数没法直接算得,那么就用邻近节点的系数插值来近似计算这些系数好了,即
H f [ u ] = 1 2 ( H C [ u ] + H F [ u ] ) = H f [ u ] B f u = 1 2 ( B C u + B F u ) = B f u D f u = 1 2 ( D C u + D F u ) = D f u \begin{aligned} H_f[u] &= \frac{1}{2}(H_C[u]+H_F[u])=\overline{H_f}[u] \\ B_f^u &= \frac{1}{2}(B_C^u+B_F^u)=\overline{B_f^u} \\ D_f^u &= \frac{1}{2}(D_C^u+D_F^u)=\overline{D_f^u} \end{aligned}
使用这些插值算得的系数,单元面处的伪动量方程变为
u f + H f [ u ] = B f u D f u ( p x ) f u_f+\overline{H_f}[u]=\overline{B_f^u} -\overline{D_f^u}\left(\frac{\partial p}{\partial x}\right)_f
这其实是使用同位网格上的动量方程系数来构建能感知“交错”网格的动量方程。

在上面的方程和后面的方程中,有上划线的值是由 C C F F 的值作线性插值而获得的,即
f = g C C + g F F \overline{\square}_f=g_C\overline{\square}_C+g_F\overline{\square}_F
其中 g C g_C g F g_F 是几何插值系数,与单元面 f f 相对节点 C C F F 的位置相关,这个在前几章已经详细讲过了。

再来看系数 H f [ u ] \overline{H_f}[u] ,其是速度的函数,想办法把它用其它量来表示,即
H f [ u ] = 1 2 ( H C [ u ] + H F [ u ] ) = 1 2 [ u C + B C u D C u ( p x ) C u F + B F u D F u ( p x ) F ] = u f D f u ( p x ) f + B f u \begin{aligned} \overline{H_f}[u] &= \frac{1}{2}(H_C[u]+H_F[u]) \\ &= \frac{1}{2}\left[ -u_C+B_C^u-D_C^u\left(\frac{\partial p}{\partial x}\right)_C -u_F+B_F^u-D_F^u\left(\frac{\partial p}{\partial x}\right)_F \right] \\ &= -\overline{u_f}-\overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f}+\overline{B_f^u} \end{aligned}
其中用到了系数的近似
1 2 [ D C u ( p x ) C + D F u ( p x ) F ] = D f u ( p x ) f D f u ( p x ) f \frac{1}{2}\left[ D_C^u\left(\frac{\partial p}{\partial x}\right)_C+ D_F^u\left(\frac{\partial p}{\partial x}\right)_F \right] =\overline{D_f^u\left(\frac{\partial p}{\partial x}\right)_f} \approx \overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f}
该系数近似有两阶精度,即
D f u ( p x ) f D f u ( p x ) f = 1 2 [ D C u ( p x ) C + D F u ( p x ) F ] 1 2 ( D C u + D F u ) × 1 2 [ ( p x ) C + ( p x ) F ] = 1 4 D C u [ ( p x ) C ( p x ) F ] + 1 4 D F u [ ( p x ) F ( p x ) C ] O ( Δ x 2 ) \begin{aligned} \overline{D_f^u\left(\frac{\partial p}{\partial x}\right)_f}- \overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f} = & \frac{1}{2}\left[ D_C^u\left(\frac{\partial p}{\partial x}\right)_C+ D_F^u\left(\frac{\partial p}{\partial x}\right)_F \right] \\ &-\frac{1}{2}(D_C^u+D_F^u)\times\frac{1}{2} \left[\left(\frac{\partial p}{\partial x}\right)_C+\left(\frac{\partial p}{\partial x}\right)_F\right] \\ =& \frac{1}{4}D_C^u \left[\left(\frac{\partial p}{\partial x}\right)_C-\left(\frac{\partial p}{\partial x}\right)_F\right] +\frac{1}{4}D_F^u \left[\left(\frac{\partial p}{\partial x}\right)_F-\left(\frac{\partial p}{\partial x}\right)_C\right] \\ \approx & O(\Delta x^2) \end{aligned}
H f [ u ] \overline{H_f}[u] 的表达式代入到单元面 f f 处的伪动量方程中,可求得以Rhie-Chow插值方法表示的单元面处的速度 u f u_f
u f = H f [ u ] + B f u D f u ( p x ) f = [ u f D f u ( p x ) f + B f u ] + B f u D f u ( p x ) f = u f a v e r a g e   v e l o c i t y D f u [ ( p x ) f ( p x ) f ] c o r r e c t i o n   t e r m \begin{aligned} u_f &= -\overline{H_f}[u]+\overline{B_f^u} -\overline{D_f^u}\left(\frac{\partial p}{\partial x}\right)_f \\ &= -\left[ -\overline{u_f}-\overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f}+\overline{B_f^u} \right] +\overline{B_f^u} -\overline{D_f^u}\left(\frac{\partial p}{\partial x}\right)_f \\ &= \underbrace{\overline{u_f}}_{average~velocity}- \underbrace{\overline{D_f^u}\left[ \left(\frac{\partial p}{\partial x}\right)_f-\overline{\left(\frac{\partial p}{\partial x}\right)_f} \right]}_{correction~term} \end{aligned}

注意,对于一维均匀网格而言,以 e e 界面为例,上式中
u e = 0.5 ( u C + u E ) \overline{u_e}=0.5(u_C+u_E)
这并无什么特殊的,而
( p x ) e = 0.5 [ ( p x ) C + ( p x ) E ] = 0.5 [ p E p W 2 Δ x + p E E p C 2 Δ x ] \overline{\left(\frac{\partial p}{\partial x}\right)_e}=0.5\left[\left(\frac{\partial p}{\partial x}\right)_C+\left(\frac{\partial p}{\partial x}\right)_E\right]=0.5\left[\frac{p_E-p_W}{2\Delta x}+\frac{p_{EE}-p_C}{2\Delta x}\right]
用的依旧是交替单元节点值,然而
( p x ) e = p E p C Δ x \left(\frac{\partial p}{\partial x}\right)_e=\frac{p_E-p_C}{\Delta x}
则用到了面两侧的邻近连续单元节点的压力值,从而有效避免了棋盘压力场的出现!

对于多维问题,也可对 y y z z 速度份量推导出类似的插值公式,即
v f = v f D f v [ ( p y ) f ( p y ) f ] w f = w f D f w [ ( p z ) f ( p z ) f ] \begin{aligned} v_f &= \overline{v_f}-\overline{D_f^v}\left[ \left(\frac{\partial p}{\partial y}\right)_f-\overline{\left(\frac{\partial p}{\partial y}\right)_f} \right] \\ w_f &= \overline{w_f}-\overline{D_f^w}\left[ \left(\frac{\partial p}{\partial z}\right)_f-\overline{\left(\frac{\partial p}{\partial z}\right)_f} \right] \end{aligned}
上面仨式子可写成矢量形式,更适用于多维压力修正方程的推导,即
v f = v f D f v ( p f p f ) \bold v_f = \overline{\bold v_f}-\overline{\bold D_f^{\bold v}}( \nabla p_f-\overline{\nabla p_f} )
其中
D f v = [ D f u 0 0 0 D f v 0 0 0 D f w ] \overline{\bold D_f^{\bold v}}= \begin{bmatrix} \overline{D_f^u} & 0 & 0 \\ 0 & \overline{D_f^v} & 0 \\ 0 & 0 & \overline{D_f^w} \end{bmatrix}
其中 p f \nabla p_f 的计算,见第9章第4节,即
p f = p f + [ p F p C d C F ( p f e C F ) ] e C F C o r r e c t i o n   t o   i n t e r p o l a t e d   f a c e   g r a d i e n t \nabla p_f=\overline{\nabla p_f}+\underbrace{\left[\frac{p_F-p_C}{d_{CF}}-(\overline{\nabla p_f}\cdot\bold e_{CF})\right]\bold e_{CF}}_{Correction~to~interpolated~face~gradient}
其中
p f = g C p C + g F p F \overline{\nabla p_f}=g_C\nabla p_C+g_F\nabla p_F
这样,可推导出在 C F CF 方向上的架构仅由邻近单元值 p F p_F p C p_C 给出
p f e C F = p f e C F + [ p F p C d C F ( p f e C F ) ] e C F e C F = p F p C d C F \begin{aligned} \nabla p_f \cdot \bold e_{CF}&=\overline{\nabla p_f} \cdot \bold e_{CF}+ \left[\frac{p_F-p_C}{d_{CF}}-(\overline{\nabla p_f}\cdot\bold e_{CF})\right]\bold e_{CF} \cdot \bold e_{CF}\\ &=\frac{p_F-p_C}{d_{CF}} \end{aligned}
因为面速度是与其邻近连续单元(而非交替单元)的压力密切相关的,因此棋盘型变量场将会被算法很好地感知,即使在同位网格上也不会出现不符合物理意义的棋盘场了。

5 通常推导

在继续讲解多维同位网格上的压力修正方程以前,首先讲解多维动量方程的离散。

5.1 动量方程离散

将动量方程
t ( ρ v ) + ( ρ v v ) = p + { μ [ v + ( v ) T ] } + f b \frac{\partial}{\partial t}(\rho \bold v)+\nabla\cdot(\rho\bold v\bold v)= -\nabla p+\nabla\cdot\left\{ \mu\left[ \nabla\bold v+(\nabla\bold v)^\text{T} \right] \right\} + \bold f_b
略做修改,写为下述形式
t ( ρ v ) + ( ρ v v ) = p + ( μ v ) + [ μ ( v ) T ] + f b \underline{\frac{\partial}{\partial t}(\rho \bold v)}+\underline{\nabla\cdot(\rho\bold v\bold v)}= -\nabla p+\underline{\nabla\cdot (\mu \nabla\bold v)}+\nabla\cdot [\mu(\nabla\bold v)^\text{T}] + \bold f_b
上式的离散形式将在时间区间 [ t Δ t / 2 , t + Δ t / 2 ] [t-\Delta t/2, t+\Delta t/2] 和单元 C C 上寻找,以下图所示。

在这里插入图片描述

上式中,三个用下划线标记的项,从左到右分别是非定常项、对流项和扩散项。这些项的离散在前面章节已经讲过了。剩下的那些项将显式估算并做为源项来对待。对于剪切应力项的第二项的体积分将使用散度定理转化为面积分,并写成面通量加和的形式,即
V C [ μ ( v ) T ] d V = V C [ μ ( v ) T ] d S = f n b ( C ) μ ( v ) T S f \int_{V_C}\nabla\cdot [\mu(\nabla\bold v)^\text{T}] dV=\int_{\partial V_C}[\mu(\nabla\bold v)^\text{T}]\cdot d\bold S= \sum_{f\sim nb(C)}\mu(\nabla\bold v)^\text{T}\cdot \bold S_f
其中 ( v ) T S f (\nabla\bold v)^\text{T}\cdot \bold S_f 在三维坐标系统中的展开形式为
( v ) T S f = [ u x S f x + u y S f y + u z S f z v x S f x + v y S f y + v z S f z w x S f x + w y S f y + w z S f z ] (\nabla\bold v)^\text{T}\cdot \bold S_f= \begin{bmatrix} \dfrac{\partial u}{\partial x}S_f^x+\dfrac{\partial u}{\partial y}S_f^y+\dfrac{\partial u}{\partial z}S_f^z \\ \dfrac{\partial v}{\partial x}S_f^x+\dfrac{\partial v}{\partial y}S_f^y+\dfrac{\partial v}{\partial z}S_f^z \\ \dfrac{\partial w}{\partial x}S_f^x+\dfrac{\partial w}{\partial y}S_f^y+\dfrac{\partial w}{\partial z}S_f^z \\ \end{bmatrix}
压力梯度的体积分也一样做为源项处理并显式估算,即
V C p d V = ( p ) C V C \int_{V_C} \nabla p dV=(\nabla p)_C V_C
或者转化成面积分来处理,即
V C p d V = V C p d S = f n b ( C ) p f S f \int_{V_C} \nabla p dV=\int_{\partial V_C}p d\bold S=\sum_{f\sim nb(C)}p_f\bold S_f
体积力项直接在控制体上积分,得
V C f b d V = ( f b ) C V C \int_{V_C} \bold f_b dV= (\bold f_b)_C V_C
使用一阶Euler格式来离散非定常项,使用高分辨率格式来离散对流项并经过延迟修正方法来将其实现,至于扩散通量,则将其分解为与网格平齐的隐式部分和显式的交叉扩散部分,这样的离散动量方程写成矢量形式为
a C v v C + F N B ( C ) a F v v F = b C v a_C^{\bold v}\bold v_C+\sum_{F\sim NB(C)}a_F^{\bold v}\bold v_F=\bold b_C^{\bold v}
其中的系数为
a C v = F l u x C C + f n b ( C ) F l u x C f a F v = F l u x F f b C v = F l u x V C f n b ( C ) F l u x V f + f n b ( C ) μ f ( v ) f T S f f r o m   t h e   2 n d   p a r t   o f   t h e   s h e a r   s t r e s s   t e r m ( p ) C V C f r o m   t h e p r e s s u r e   g r a d i e n t \begin{aligned} a_C^{\bold v} &= FluxC_C+\sum_{f\sim nb(C)}FluxC_f \\ a_F^{\bold v} &= FluxF_f \\ \bold b_C^{\bold v} &= -FluxV_C-\sum_{f\sim nb(C)} FluxV_f \\ &+ \underbrace{\sum_{f\sim nb(C)}\mu_f(\nabla\bold v)_f^\text{T}\cdot \bold S_f} _{\footnotesize{\begin{matrix}from~the~2nd~part~of\\~the~shear~stress~term\end{matrix}}} - \underbrace{(\nabla p)_C V_C}_ {\footnotesize{\begin{matrix}from~the\\pressure~gradient\end{matrix}}} \end{aligned}
其中的面通量用下式计算
F l u x C f = m ˙ f , 0 c o n v e c t i o n c o n t r i b u t i o n + μ f E f d C F d i f f u s i o n c o n t r i b u t i o n F l u x F f = m ˙ f , 0 c o n v e c t i o n c o n t r i b u t i o n μ f E f d C F d i f f u s i o n c o n t r i b u t i o n F l u x V f = μ f ( v ) f T f c r o s s   d i f f u s i o n c o n t r i b u t i o n + m ˙ f ( v f H R v f U ) h i g h   r e s o l u t i o n c o n v e c t i o n   s c h e m e c o n t r i b u t i o n \begin{aligned} FluxC_f &= \underbrace{||\dot m_f, 0||}_{\footnotesize{\begin{matrix}convection \\ contribution\end{matrix}}} + \underbrace{\mu_f \frac{E_f}{d_{CF}}}_{\footnotesize{\begin{matrix}diffusion\\contribution\end{matrix}}} \\ FluxF_f &= -\underbrace{||-\dot m_f, 0||}_{\footnotesize{\begin{matrix}convection \\ contribution\end{matrix}}} - \underbrace{\mu_f \frac{E_f}{d_{CF}}}_{\footnotesize{\begin{matrix}diffusion\\contribution\end{matrix}}} \\ FluxV_f &= \underbrace{-\mu_f(\nabla\bold v)_f\cdot \bold T_f}_{\footnotesize{\begin{matrix}cross~diffusion\\contribution\end{matrix}}} + \underbrace{\dot m_f(\bold v_f^{HR}-\bold v_f^U)}_{\footnotesize{\begin{matrix}high~resolution\\convection~scheme\\contribution\end{matrix}}} \end{aligned}
单元通量则用下式计算
F l u x C C = ρ C V C Δ t F l u x V C = ρ C V C Δ t v C t r a n s i e n t c o n t r i b u t i o n ( f b ) C V C s o u r c e   t e r m c o n t r i b u t i o n \begin{aligned} FluxC_C &= \frac{\rho_C V_C}{\Delta t} \\ FluxV_C &= \underbrace{-\frac{\rho_C^\circ V_C}{\Delta t}\bold v_C^\circ}_ {\footnotesize{\begin{matrix}transient \\ contribution\end{matrix}}}- \underbrace{(\bold f_b)_C V_C}_{\footnotesize{\begin{matrix}source~term\\contribution\end{matrix}}} \end{aligned}
尽管动量方程的离散代数方程是线性的,然而其系数是依赖于速度和压力场的。该非线性是要经过迭代过程来处理的,即,在每次迭代开始时,系数是基于前次迭代所得到的因变量的值而计算的。这些系数值的变化会引发 v \bold v 的巨大变化,从而影响收敛速率,甚至会致使发散。为了减慢变化,可应用欠松弛方法,尤为是所用的瞬态时间步长较大时。将欠松弛因子定义为 λ v \lambda^{\bold v} ,并采用Patankar隐式欠松弛手法,则欠松弛动量方程可写为
a C v λ v v C + F N B ( C ) a F v v F = b C v + 1 λ v λ v a C v v C ( n ) \frac{a_C^{\bold v}}{\lambda^{\bold v}}\bold v_C+\sum_{F\sim NB(C)}a_F^{\bold v}\bold v_F=\bold b_C^{\bold v}+ \frac{1-\lambda^{\bold v}}{\lambda^{\bold v}}a_C^{\bold v}\bold v_C^{(n)}
经过从新定义 a C v a_C^{\bold v} b C v b_C^{\bold v} ,即
a C v a C v λ v b C v b C v + 1 λ v λ v a C v v C ( n ) \begin{aligned} a_C^{\bold v} &\leftarrow \frac{a_C^{\bold v}}{\lambda^{\bold v}} \\ \bold b_C^{\bold v} &\leftarrow \bold b_C^{\bold v}+ \frac{1-\lambda^{\bold v}}{\lambda^{\bold v}}a_C^{\bold v}\bold v_C^{(n)} \end{aligned}
欠松弛动量方程可从新写为其原来的样子
a C v v C + F N B ( C ) a F v v F = b C v a_C^{\bold v}\bold v_C+\sum_{F\sim NB(C)}a_F^{\bold v}\bold v_F=\bold b_C^{\bold v}
为了推导出同位网格上的压力修正方程,将压力梯度从源项 b C v \bold b_C^{\bold v} 中剥离出来,并显式计算,即
b C v = V C ( p ) C + b ^ C v \bold b_C^{\bold v}=-V_C(\nabla p)_C+\hat{\bold b}_C^{\bold v}
因而,动量方程最终变成了下式
v C + F N B ( C ) a F v a C v v F = V C a C v ( p ) C + b ^ C v a C v \bold v_C+\sum_{F\sim NB(C)}\frac{a_F^{\bold v}}{a_C^{\bold v}}\bold v_F= -\frac{V_C}{a_C^{\bold v}}(\nabla p)_C+\frac{\hat{\bold b}_C^{\bold v}}{a_C^{\bold v}}
定义以下矢量算子
H C [ v ] = F N B ( C ) a F v a C v v F B C v = b ^ C v a C v D C v = V C a C v \begin{aligned} \bold H_C[\bold v] &= \sum_{F\sim NB(C)}\frac{a_F^{\bold v}}{a_C^{\bold v}}\bold v_F \\ \bold B_C^{\bold v} &= \frac{\hat{\bold b}_C^{\bold v}}{a_C^{\bold v}} \\ \bold D_C^{\bold v} &= \frac{V_C}{a_C^{\bold v}} \end{aligned}
动量方程整理为
v C + H C [ v ] = D C v ( p ) C + B C v \bold v_C+\bold H_C[\bold v] = -\bold D_C^{\bold v}(\nabla p)_C + \bold B_C^{\bold v}
该形式将用于后续的推导中。

5.2 同位网格上的压力修正方程

与交错网格的时候同样,最开始用猜想值或者前次迭代值 ( v ( n ) , m ˙ ( n ) , p ( n ) ) (\bold v^{(n)},\dot m^{(n)},p^{(n)}) ,求解动量方程(上节最后推出的式子)
v C + H C [ v ] = D C v ( p ) C + B C v \bold v_C+\bold H_C[\bold v] = -\bold D_C^{\bold v}(\nabla p)_C + \bold B_C^{\bold v}
来得到知足动量方程的速度场 v \bold v^* ,这样所得解知足
v C + H C [ v ] = D C v ( p ( n ) ) C + B C v \bold v_C^*+\bold H_C[\bold v^*] = -\bold D_C^{\bold v}(\nabla p^{(n)})_C + \bold B_C^{\bold v}
最终的解是要知足上上式的,那么这两个方程(上式和上上式)的差异在于,知足上上式的速度场既知足动量方程又知足连续方程,而知足上式的速度场则只知足动量方程而未必知足连续方程,由于上式作了线性化处理,其计算系数和压力梯度项所用的压力和速度是基于前次迭代的值。所以,须要对速度、压力、质量流量场进行修正,以使它们知足质量守恒。将这些修正量用 ( v , p , m ˙ ) (\bold v',p',\dot m') 表示,则精确解和计算值之间的关系可写为
v = v + v p = p ( n ) + p m ˙ = m ˙ + m ˙ \begin{aligned} \bold v &= \bold v^* + \bold v' \\ p &= p^{(n)} + p' \\ \dot m &= \dot m^* + \dot m' \end{aligned}
将质量流量的修正方程代入到连续方程中,得
f n b ( C ) m ˙ f = f n b ( C ) m ˙ f m ˙ f = ρ f v f S f \sum_{f\sim nb(C)}\dot m_f' = -\sum_{f\sim nb(C)} \dot m_f^* \quad\quad \dot m_f^* = \rho_f \bold v_f^* \cdot \bold S_f
速度 v f \bold v_f^* 使用Rhie-Chow插值计算,即
v f = v f D f v ( p f ( n ) p f ( n ) ) \bold v_f^* = \overline{\bold v_f^*}-\overline{\bold D_f^{\bold v}}( \nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}} )
当算得的质量流量场知足守恒时, f n b ( C ) m ˙ f = f n b ( C ) m ˙ f \sum_{f\sim nb(C)}\dot m_f' = -\sum_{f\sim nb(C)} \dot m_f^* 的右端项为零,代表无需修正流场。另外一方面,若速度场不正确的话将会产生不平衡的质量流量,并使得该式右端项非零,代表须要修正流场以使其知足质量守恒条件。

质量流量修正项可写成速度修正项的形式,先用精确解的动量方程 v C + H C [ v ] = D C v ( p ) C + B C v \bold v_C+\bold H_C[\bold v] = -\bold D_C^{\bold v}(\nabla p)_C + \bold B_C^{\bold v} 减去计算值的动量方程 v C + H C [ v ] = D C v ( p ( n ) ) C + B C v \bold v_C^*+\bold H_C[\bold v^*] = -\bold D_C^{\bold v}(\nabla p^{(n)})_C + \bold B_C^{\bold v} ,便可
v C + H C [ v ] = D C v ( p ) C \bold v_C'+\bold H_C[\bold v'] = -\bold D_C^{\bold v}(\nabla p')_C
一样也可推导出 F F 单元的速度修正值动量方程
v F + H F [ v ] = D F v ( p ) F \bold v_F'+\bold H_F[\bold v'] = -\bold D_F^{\bold v}(\nabla p')_F
单元面上的质量流量修正可写为
m ˙ f = ρ f v f S f \dot m_f'=\rho_f \bold v_f' \cdot \bold S_f
其中面速度修正是由 v f = v f D f v ( p f p f ) \bold v_f = \overline{\bold v_f}-\overline{\bold D_f^{\bold v}}( \nabla p_f-\overline{\nabla p_f} ) 减去 v f = v f D f v ( p f ( n ) p f ( n ) ) \bold v_f^*=\overline{\bold v_f^*}-\overline{\bold D_f^{\bold v}}( \nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}}) 而获得的
v f = v f D f v ( ( p ) f p f ) \bold v_f' = \overline{\bold v_f'}-\overline{\bold D_f^{\bold v}}( (\nabla p')_f-\overline{\nabla p_f'})
把上式和上上式代入到连续方程 f n b ( C ) m ˙ f = f n b ( C ) m ˙ f \sum_{f\sim nb(C)}\dot m_f' = -\sum_{f\sim nb(C)} \dot m_f^* 中,推导出以下的压力修正方程
f n b ( C ) ( ρ f v f S f ) + f n b ( C ) ( ρ f D f v   p f S f ) f n b ( C ) ( ρ f D f v ( p ) f S f ) = f n b ( C ) m ˙ f \begin{aligned} \underline{ \sum_{f\sim nb(C)}\left(\rho_f \overline{\bold v_f'} \cdot \bold S_f \right) +\sum_{f\sim nb(C)}\left(\rho_f \overline{\bold D_f^{\bold v}}~\overline{\nabla p_f'} \cdot \bold S_f\right) } -\sum_{f\sim nb(C)}\left(\rho_f \overline{\bold D_f^{\bold v}} (\nabla p')_f \cdot \bold S_f\right) \\ = -\sum_{f\sim nb(C)} \dot m_f^* \end{aligned}
在该方程中,下划线部分表明了邻近速度对于所考虑单元的速度修正的影响。为更清晰地理解该影响,可经过将前面的单元速度修正方程式 v C + H C [ v ] = D C v ( p ) C \bold v_C'+\bold H_C[\bold v'] = -\bold D_C^{\bold v}(\nabla p')_C 和式 v F + H F [ v ] = D F v ( p ) F \bold v_F'+\bold H_F[\bold v'] = -\bold D_F^{\bold v}(\nabla p')_F 二者插值到面上去,推出上式下划线部分的以下的等效表达式。
v f + H f [ v ] = D f v   ( p ) f v f + D f v   ( p ) f = H f [ v ] \begin{aligned} & \overline{\bold v_f'}+\overline{\bold H_f}[\bold v'] = -\overline{\bold D_f^{\bold v}}~\overline{({\nabla p'})_f} \\ \Rightarrow & \overline{\bold v_f'}+\overline{\bold D_f^{\bold v}}~\overline{({\nabla p'})_f}=-\overline{\bold H_f}[\bold v'] \end{aligned}
将上式代入到上上式中,得
f n b ( C ) ( ρ f D f v ( p ) f S f ) = f n b ( C ) m ˙ f + f n b ( C ) ( ρ f H f [ v ] S f ) \sum_{f\sim nb(C)}\left(-\rho_f \overline{\bold D_f^{\bold v}} (\nabla p')_f \cdot \bold S_f\right) = -\sum_{f\sim nb(C)} \dot m_f^*+ \underline{\sum_{f\sim nb(C)}\left(\rho_f \overline{\bold H_f}[\bold v'] \cdot \bold S_f\right) }
或者写成更加显式的形式
f n b ( C ) ( ρ f D f v ( p ) f S f ) = f n b ( C ) m ˙ f + f n b ( C ) ( ρ f ( F N B ( C ) a F v a C v v F ) S f ) \sum_{f\sim nb(C)}\left(-\rho_f \overline{\bold D_f^{\bold v}} (\nabla p')_f \cdot \bold S_f\right) = -\sum_{f\sim nb(C)} \dot m_f^*+ \underline{\sum_{f\sim nb(C)}\left(\rho_f \left(\overline{\sum_{F\sim NB(C)}\frac{a_F^{\bold v}}{a_C^{\bold v}}\bold v'_F} \right) \cdot \bold S_f\right) }\\
在上式或上上式中,对于下划线部分的处理相当重要,其决定着方程可解与否。在最初的SIMPLE算法中,是直接把该项忽略掉的,即,认为在某点处的速度修正是只跟压力相关,而跟周围的速度绝不相干。因为这只是个修正方程,因此对其中的某些项做以修改甚至是直接扔掉也不会影响到最终的解,由于在收敛状况下修正会变成零的(能够这么来理解,对于某点速度(压力)的修正,一方面来自于压力,一方面来自于邻近速度,那么若是认为周围速度的影响是零(若是收敛的话的确如此),只考虑压力的影响,至关于只考虑了部分修正,同样能够到达最终的收敛状态(此时压力和邻近速度的影响都是零了))。这么个近似虽然对最终收敛解没有影响,可是它会影响到收敛的速度哈,由于在每步迭代中,忽略的项越大,那么带来的近似偏差就越高了。

上式和上上式中,剩余项的处理就很简单了,压力修正方程的系数,是按照第8章中所讲的扩散项的离散方法来得到的,具体来讲,是非各向同性(各向异性)扩散的处理方法。

能够把左端项修改为梯度点积形式
( D f v ( p ) f ) S f = ( ( p ) f D f v T ) S f = ( p ) f ( D f v T S f ) = ( p ) f S f \begin{aligned} \left( \overline{\bold D_f^{\bold v}} (\nabla p')_f\right) \cdot \bold S_f &= \left( (\nabla p')_f \overline{\bold D_f^{\bold v}}^T \right) \cdot \bold S_f \\ &= (\nabla p')_f \cdot \left( \overline{\bold D_f^{\bold v}}^T \bold S_f \right)\\ &= (\nabla p')_f \cdot \bold S'_f \end{aligned}
其中 S f \bold S'_f 的展开形式为
S f = D f v T S f = [ D f u 0 0 0 D f v 0 0 0 D f w ] [ S f x S f y S f z ] = [ D f u S f x D f v S f y D f w S f z ] \bold S'_f=\overline{\bold D_f^{\bold v}}^T \bold S_f= \begin{bmatrix} \overline{D_f^u} & 0 & 0 \\ 0 & \overline{D_f^v} & 0 \\ 0 & 0 & \overline{D_f^w} \end{bmatrix} \begin{bmatrix} S_f^x \\ S_f^y \\ S_f^z \end{bmatrix}= \begin{bmatrix} \overline{D_f^u}S_f^x \\ \overline{D_f^v}S_f^y \\ \overline{D_f^w}S_f^z \end{bmatrix}
使用 S f \bold S'_f ,压力修正梯度项的离散过程和以前同样,获得
( p ) f S f = ( p ) f E f + ( p ) f T f = E f d C F ( p F p C ) + ( p ) f T f \begin{aligned} (\nabla p')_f \cdot \bold S'_f &= (\nabla p')_f \cdot \bold E_f + (\nabla p')_f \cdot \bold T_f \\ &= \frac{E_f}{d_{CF}}(p'_F-p'_C)+\underline{ (\nabla p')_f\cdot\bold T_f} \end{aligned}
其中用到了对 S f \bold S'_f 以下分解
S f = E f + T f \bold S'_f=\bold E_f+\bold T_f
分解的类型在第8章中已经讲过(三种类型:最小修正、正交修正、过松弛修正方法)。下划线项是由网格非正交特性所引发的,能够选择忽略或保留。若是忽略,其不会影响最终解,由于这只是修正项。若是保留,其将会在内循环中显式处理(OpenFOAM中的非正交循环)。因为每次迭代的求解都是从零压力修正场开始的,在求解方程时该项须要迭代更新。

忽略掉非正交项的做用,压力修正方程的线性化形式为
( p ) f S f