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 Flows

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

由于本章内容繁杂，篇幅较长，故分成了四部分来讲解，各部分主要内容分别为：交错网格、同位网格、边界条件、SIMPLE家族算法。

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

3 交错网格的缺陷

交错网格的使用对于SIMPLE算法的发展来说至关重要，然而采用交错网格架构也有其缺点。如前所述，在二维和三维情况下，分别需要三套和四套交错网格系统，来存放速度的两到三个不同分量和压力值。

这样一来，就需要对每个速度分量存储一套网格系统，对压力和其他变量再存储一套网格系统，内存开销是蛮大的。除此之外，对于非笛卡尔网格，交错网格系统会产生很多问题，更有甚者，对于非结构网格而言，交错网格系统基本上是无法构造出来的。

在曲线网格中，使用笛卡尔速度分量会导致的问题是，当一个或多个面和交错速度分量的方向平齐而非垂直的时候，比如上图中的最上端的两个单元中，由于每个面上的交错速度与面平齐，算得的质量流量都几乎是零，这样造成其错误地自动满足了连续方程，显然这是不符合物理意义的，会导致SIMPLE计算过程的失效。

因此，这种情况下较好的替代方法是使用协变的（covariant）或抗变（contra-variant）速度分量，如上图a和b所示。

使用抗变速度分量的交错网格例子如下图。

不幸的是，在曲线坐标系上离散动量方程时，复杂性大为提升，表现在对扩散项的处理上，因为方程中有非守恒项。

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

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

结果是咱们还得回到最初的单元中心的同位网格系统上，如下图所示，所有的变量都存在相同的位置处（单元形心），使用这种同位网格系统才是最佳方法（兜兜转转饶了一大圈又转回来了哈）。值得注意的是，虽然速度分量是存储在单元形心处的，与压力和其他变量存储位置一样，然而质量流量（是个标量）在同位网格中是存储在面单元上的。质量流量实际上可以视作抗变（contra-variant）分量，只是此时质量流量的计算通常是用离散动量方程来插值计算的（不再像之前那样子用界面速度，而界面速度又用两侧单元值插值，从而出现了连续方程中用到的是交替单元的速度的情况，无法感知棋盘型的非均匀流场），这便是著名的Rhie-Chow插值，其正是下面章节的主题。

4 Rhie-Chow插值

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

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

（说到这，连我都被绕晕了，没关系，接下来看看如何处理就会豁然开朗了）

由单元$C$和$F$的离散$x$动量方程开始，即
$$\begin{gathered} 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 \end{gathered}$$
$u_f$速度方程与上式类似，只是压力梯度是与局部邻近压力值相关的（$f$处压力梯度是用$C$和$F$压力值计算的，并非间隔了一个单元的交替值），如上图所示，$u_f$方程的形式为
$$u_f+H_f[u]=B_f^u-D_f^u{\left(\frac{\partial p}{\partial x}\right)}_f$$
既然在同位网格上，上式中的系数无法直接算得，那么就用邻近节点的系数插值来近似计算这些系数好了，即
$$\begin{array}{rl}{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{array}$$
使用这些插值算得的系数，单元面处的伪动量方程变为
$$u_f+\overline{H_f}[u]=\overline{B_f^u}-\overline{D_f^u}{\left(\frac{\partial p}{\partial x}\right)}_f$$
这实际上是使用同位网格上的动量方程系数来构建能感知“交错”网格的动量方程。

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

再来看系数$\overline{H_f}[u]$，其是速度的函数，想办法把它用其它量来表示，即
$$\begin{array}{rl}\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{array}$$
其中用到了系数的近似
$$\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}$$
该系数近似有两阶精度，即
$$\begin{array}{rl}\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{array}$$
把$\overline{H_f}[u]$的表达式代入到单元面$f$处的伪动量方程中，可求得以Rhie-Chow插值方法表示的单元面处的速度$u_f$为
$$\begin{array}{rl}{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\\ & =\underset{averagevelocity}{\underbrace{\overline{u_f}}}-\underset{correctionterm}{\underbrace{\overline{D_f^u}\left[{\left(\frac{\partial p}{\partial x}\right)}_f-\overline{{\left(\frac{\partial p}{\partial x}\right)}_f}\right]}}\end{array}$$

注意，对于一维均匀网格而言，以$e$界面为例，上式中
$$\overline{u_e}=0.5(u_C+u_E)$$
这并没有什么特殊的，而
$$\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]$$
用的依旧是交替单元节点值，然而
$${\left(\frac{\partial p}{\partial x}\right)}_e=\frac{p_E-p_C}{\Delta x}$$
则用到了面两侧的邻近连续单元节点的压力值，从而有效避免了棋盘压力场的出现！

对于多维问题，也可对$y$和$z$速度分量推导出相似的插值公式，即
$$\begin{array}{rl}{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{array}$$
上面仨式子可写成矢量形式，更适用于多维压力修正方程的推导，即
$${\mathbf{v}}_f=\overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})$$
其中
$$\overline{{\mathbf{D}}_f^{\mathbf{v}}}=\left[\begin{array}{ccc}\overline{D_f^u}& 0& 0\\ 0& \overline{D_f^v}& 0\\ 0& 0& \overline{D_f^w}\end{array}\right]$$
其中$\nabla p_f$的计算，见第9章第4节，即
$$\nabla p_f=\overline{\nabla p_f}+\underset{Correctiontointerpolatedfacegradient}{\underbrace{\left[\frac{p_F-p_C}{d_{CF}}-(\overline{\nabla p_f}\cdot {\mathbf{e}}_{CF})\right]{\mathbf{e}}_{CF}}}$$
其中
$$\overline{\nabla p_f}=g_C\nabla p_C+g_F\nabla p_F$$
这样，可推导出在$CF$方向上的架构仅由邻近单元值$p_F$和$p_C$给出
$$\begin{array}{rl}\nabla p_f\cdot {\mathbf{e}}_{CF}& =\overline{\nabla p_f}\cdot {\mathbf{e}}_{CF}+\left[\frac{p_F-p_C}{d_{CF}}-(\overline{\nabla p_f}\cdot {\mathbf{e}}_{CF})\right]{\mathbf{e}}_{CF}\cdot {\mathbf{e}}_{CF}\\ & =\frac{p_F-p_C}{d_{CF}}\end{array}$$
由于面速度是与其邻近连续单元（而非交替单元）的压力密切相关的，所以棋盘型变量场将会被算法很好地感知，即便在同位网格上也不会出现不符合物理意义的棋盘场了。

5 一般推导

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

5.1 动量方程离散

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

上式中，三个用下划线标记的项，从左到右分别是非定常项、对流项和扩散项。这些项的离散在前面章节已经讲过了。剩下的那些项将显式估算并作为源项来对待。对于剪切应力项的第二项的体积分将使用散度定理转化为面积分，并写成面通量加和的形式，即
$${\int }_{V_C}\nabla \cdot [\mu (\nabla \mathbf{v}{)}^{\text{T}}]dV={\int }_{\partial V_C}[\mu (\nabla \mathbf{v}{)}^{\text{T}}]\cdot d\mathbf{S}=\sum _{f\sim nb(C)}\mu (\nabla \mathbf{v}{)}^{\text{T}}\cdot {\mathbf{S}}_f$$
其中$(\nabla \mathbf{v}{)}^{\text{T}}\cdot {\mathbf{S}}_f$在三维坐标系统中的展开形式为
$$(\nabla \mathbf{v}{)}^{\text{T}}\cdot {\mathbf{S}}_f=\left[\begin{array}{c}\frac{\partial u}{\partial x}{S}_f^x+\frac{\partial u}{\partial y}{S}_f^y+\frac{\partial u}{\partial z}{S}_f^z\\ \frac{\partial v}{\partial x}{S}_f^x+\frac{\partial v}{\partial y}{S}_f^y+\frac{\partial v}{\partial z}{S}_f^z\\ \frac{\partial w}{\partial x}{S}_f^x+\frac{\partial w}{\partial y}{S}_f^y+\frac{\partial w}{\partial z}{S}_f^z\end{array}\right]$$
压力梯度的体积分也同样作为源项处理并显式估算，即
$${\int }_{V_C}\nabla pdV=(\nabla p{)}_C{V}_C$$
或者转化成面积分来处理，即
$${\int }_{V_C}\nabla pdV={\int }_{\partial V_C}pd\mathbf{S}=\sum _{f\sim nb(C)}{p}_f{\mathbf{S}}_f$$
体积力项直接在控制体上积分，得
$${\int }_{V_C}{\mathbf{f}}_{b}dV=({\mathbf{f}}_b{)}_C{V}_C$$
使用一阶Euler格式来离散非定常项，使用高分辨率格式来离散对流项并通过延迟修正方法来将其实现，至于扩散通量，则将其分解为与网格平齐的隐式部分和显式的交叉扩散部分，这样的离散动量方程写成矢量形式为
$$a_C^{\mathbf{v}}{\mathbf{v}}_C+\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F={\mathbf{b}}_C^{\mathbf{v}}$$
其中的系数为
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& =Flux{C}_C+\sum _{f\sim nb(C)}Flux{C}_f\\ a_F^{\mathbf{v}}& =Flux{F}_f\\ {\mathbf{b}}_C^{\mathbf{v}}& =-Flux{V}_C-\sum _{f\sim nb(C)}Flux{V}_f\\ & +\underset{\begin{array}{c}fromthe2ndpartof\\ theshearstressterm\end{array}}{\underbrace{\sum _{f\sim nb(C)}{\mu }_f(\nabla \mathbf{v}{)}_f^{\text{T}}\cdot {\mathbf{S}}_f}}-\underset{\begin{array}{c}fromthe\\ pressuregradient\end{array}}{\underbrace{(\nabla p{)}_C{V}_C}}\end{array}$$
其中的面通量用下式计算
$$\begin{array}{rl}Flux{C}_f& =\underset{\begin{array}{c}convection\\ contribution\end{array}}{\underbrace{||{\dot{m}}_f,0||}}+\underset{\begin{array}{c}diffusion\\ contribution\end{array}}{\underbrace{{\mu }_f\frac{E_f}{d_{CF}}}}\\ Flux{F}_f& =-\underset{\begin{array}{c}convection\\ contribution\end{array}}{\underbrace{||-{\dot{m}}_f,0||}}-\underset{\begin{array}{c}diffusion\\ contribution\end{array}}{\underbrace{{\mu }_f\frac{E_f}{d_{CF}}}}\\ Flux{V}_f& =\underset{\begin{array}{c}crossdiffusion\\ contribution\end{array}}{\underbrace{-{\mu }_f(\nabla \mathbf{v}{)}_f\cdot {\mathbf{T}}_f}}+\underset{\begin{array}{c}highresolution\\ convectionscheme\\ contribution\end{array}}{\underbrace{{\dot{m}}_f({\mathbf{v}}_f^{HR}-{\mathbf{v}}_f^U)}}\end{array}$$
单元通量则用下式计算
$$\begin{array}{rl}Flux{C}_C& =\frac{{\rho }_C{V}_C}{\Delta t}\\ Flux{V}_C& =\underset{\begin{array}{c}transient\\ contribution\end{array}}{\underbrace{-\frac{{\rho }_C^{\circ }{V}_C}{\Delta t}{\mathbf{v}}_C^{\circ }}}-\underset{\begin{array}{c}sourceterm\\ contribution\end{array}}{\underbrace{({\mathbf{f}}_b{)}_C{V}_C}}\end{array}$$
尽管动量方程的离散代数方程是线性的，然而其系数是依赖于速度和压力场的。该非线性是要通过迭代过程来处理的，即，在每次迭代开始时，系数是基于前次迭代所获得的因变量的值而计算的。这些系数值的变化会引起$\mathbf{v}$的巨大变化，从而影响收敛速率，甚至会导致发散。为了减慢变化，可应用欠松弛方法，尤其是所用的瞬态时间步长较大时。将欠松弛因子定义为${\lambda }^{\mathbf{v}}$，并采用Patankar隐式欠松弛手法，则欠松弛动量方程可写为
$$\frac{a_C^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}{\mathbf{v}}_C+\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F={\mathbf{b}}_C^{\mathbf{v}}+\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}{a}_C^{\mathbf{v}}{\mathbf{v}}_C^{(n)}$$
通过重新定义$a_C^{\mathbf{v}}$和$b_C^{\mathbf{v}}$，即
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& \leftarrow \frac{a_C^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\\ {\mathbf{b}}_C^{\mathbf{v}}& \leftarrow {\mathbf{b}}_C^{\mathbf{v}}+\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}{a}_C^{\mathbf{v}}{\mathbf{v}}_C^{(n)}\end{array}$$
欠松弛动量方程可重新写为其原来的样子
$$a_C^{\mathbf{v}}{\mathbf{v}}_C+\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F={\mathbf{b}}_C^{\mathbf{v}}$$
为了推导出同位网格上的压力修正方程，将压力梯度从源项${\mathbf{b}}_C^{\mathbf{v}}$中剥离出来，并显式计算，即
$${\mathbf{b}}_C^{\mathbf{v}}=-V_C(\nabla p{)}_C+{\hat{\mathbf{b}}}_C^{\mathbf{v}}$$
于是，动量方程最终变成了下式
$${\mathbf{v}}_C+\sum _{F\sim NB(C)}\frac{a_F^{\mathbf{v}}}{a_C^{\mathbf{v}}}{\mathbf{v}}_F=-\frac{V_C}{a_C^{\mathbf{v}}}(\nabla p{)}_C+\frac{{\hat{\mathbf{b}}}_C^{\mathbf{v}}}{a_C^{\mathbf{v}}}$$
定义如下矢量算子
$$\begin{array}{rl}{\mathbf{H}}_C[\mathbf{v}]& =\sum _{F\sim NB(C)}\frac{a_F^{\mathbf{v}}}{a_C^{\mathbf{v}}}{\mathbf{v}}_F\\ {\mathbf{B}}_C^{\mathbf{v}}& =\frac{{\hat{\mathbf{b}}}_C^{\mathbf{v}}}{a_C^{\mathbf{v}}}\\ {\mathbf{D}}_C^{\mathbf{v}}& =\frac{V_C}{a_C^{\mathbf{v}}}\end{array}$$
动量方程整理为
$${\mathbf{v}}_C+{\mathbf{H}}_C[\mathbf{v}]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$$
该形式将用于后续的推导中。

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

与交错网格的时候一样，最开始用猜测值或者前次迭代值$({\mathbf{v}}^{(n)},{\dot{m}}^{(n)},p^{(n)})$，求解动量方程（上节最后推出的式子）
$${\mathbf{v}}_C+{\mathbf{H}}_C[\mathbf{v}]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$$
来获得满足动量方程的速度场${\mathbf{v}}^{\ast }$，这样所得解满足
$${\mathbf{v}}_C^{\ast }+{\mathbf{H}}_C[{\mathbf{v}}^{\ast }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{(n)}{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$$
最终的解是要满足上上式的，那么这两个方程（上式和上上式）的差别在于，满足上上式的速度场既满足动量方程又满足连续方程，而满足上式的速度场则只满足动量方程而未必满足连续方程，因为上式做了线性化处理，其计算系数和压力梯度项所用的压力和速度是基于前次迭代的值。因此，需要对速度、压力、质量流量场进行修正，以使它们满足质量守恒。将这些修正量用$({\mathbf{v}}^{\prime },p^{\prime },{\dot{m}}^{\prime })$表示，则精确解和计算值之间的关系可写为
$$\begin{array}{rl}\mathbf{v}& ={\mathbf{v}}^{\ast }+{\mathbf{v}}^{\prime }\\ p& =p^{(n)}+p^{\prime }\\ \dot{m}& ={\dot{m}}^{\ast }+{\dot{m}}^{\prime }\end{array}$$
将质量流量的修正方程代入到连续方程中，得
$$\sum _{f\sim nb(C)}{\dot{m}}_f^{\prime }=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }{\dot{m}}_f^{\ast }={\rho }_f{\mathbf{v}}_f^{\ast }\cdot {\mathbf{S}}_f$$
速度${\mathbf{v}}_f^{\ast }$使用Rhie-Chow插值计算，即
$${\mathbf{v}}_f^{\ast }=\overline{{\mathbf{v}}_f^{\ast }}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}})$$
当算得的质量流量场满足守恒时，${\sum }_{f\sim nb(C)}{\dot{m}}_f^{\prime }=-{\sum }_{f\sim nb(C)}{\dot{m}}_f^{\ast }$的右端项为零，表明无需修正流场。另一方面，若速度场不正确的话将会产生不平衡的质量流量，并使得该式右端项非零，表明需要修正流场以使其满足质量守恒条件。

质量流量修正项可写成速度修正项的形式，先用精确解的动量方程${\mathbf{v}}_C+{\mathbf{H}}_C[\mathbf{v}]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$减去计算值的动量方程${\mathbf{v}}_C^{\ast }+{\mathbf{H}}_C[{\mathbf{v}}^{\ast }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{(n)}{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$，即可
$${\mathbf{v}}_C^{\prime }+{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C$$
同样也可推导出$F$单元的速度修正值动量方程
$${\mathbf{v}}_F^{\prime }+{\mathbf{H}}_F[{\mathbf{v}}^{\prime }]=-{\mathbf{D}}_F^{\mathbf{v}}(\nabla p^{\prime }{)}_F$$
单元面上的质量流量修正可写为
$${\dot{m}}_f^{\prime }={\rho }_f{\mathbf{v}}_f^{\prime }\cdot {\mathbf{S}}_f$$
其中面速度修正是由${\mathbf{v}}_f=\overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})$减去${\mathbf{v}}_f^{\ast }=\overline{{\mathbf{v}}_f^{\ast }}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}})$而得到的
$${\mathbf{v}}_f^{\prime }=\overline{{\mathbf{v}}_f^{\prime }}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}((\nabla p^{\prime }{)}_f-\overline{\nabla p_f^{\prime }})$$
把上式和上上式代入到连续方程${\sum }_{f\sim nb(C)}{\dot{m}}_f^{\prime }=-{\sum }_{f\sim nb(C)}{\dot{m}}_f^{\ast }$中，推导出如下的压力修正方程
$$\begin{array}{r}\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{v}}_f^{\prime }}\cdot {\mathbf{S}}_f\right)+\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\overline{\nabla p_f^{\prime }}\cdot {\mathbf{S}}_f\right)}-\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f\right)\\ =-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }\end{array}$$
在该方程中，下划线部分代表了邻近速度对于所考虑单元的速度修正的影响。为更清晰地理解该影响，可通过将前面的单元速度修正方程式${\mathbf{v}}_C^{\prime }+{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C$和式${\mathbf{v}}_F^{\prime }+{\mathbf{H}}_F[{\mathbf{v}}^{\prime }]=-{\mathbf{D}}_F^{\mathbf{v}}(\nabla p^{\prime }{)}_F$两者插值到面上去，推出上式下划线部分的如下的等效表达式。
$$\begin{array}{rl}& \overline{{\mathbf{v}}_f^{\prime }}+\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]=-\overline{{\mathbf{D}}_f^{\mathbf{v}}}\overline{(\nabla p^{\prime }{)}_f}\\ \Rightarrow & \overline{{\mathbf{v}}_f^{\prime }}+\overline{{\mathbf{D}}_f^{\mathbf{v}}}\overline{(\nabla p^{\prime }{)}_f}=-\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\end{array}$$
将上式代入到上上式中，得
$$\sum _{f\sim nb(C)}\left(-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f\right)=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}$$
或者写成更加显式的形式
$$\sum _{f\sim nb(C)}\left(-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f\right)=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\left(\overline{\sum _{F\sim NB(C)}\frac{a_F^{\mathbf{v}}}{a_C^{\mathbf{v}}}{\mathbf{v}}_F^{\prime }}\right)\cdot {\mathbf{S}}_f\right)}$$
在上式或上上式中，对于下划线部分的处理至关重要，其决定着方程可解与否。在最初的SIMPLE算法中，是直接把该项忽略掉的，即，认为在某点处的速度修正是只跟压力相关，而跟周围的速度毫不相干。由于这只是个修正方程，所以对其中的某些项作以修改甚至是直接扔掉也不会影响到最终的解，因为在收敛情况下修正会变成零的（可以这么来理解，对于某点速度（压力）的修正，一方面来自于压力，一方面来自于邻近速度，那么如果认为周围速度的影响是零（如果收敛的话的确如此），只考虑压力的影响，相当于只考虑了部分修正，一样可以到达最终的收敛状态（此时压力和邻近速度的影响都是零了））。这么个近似虽然对最终收敛解没有影响，但是它会影响到收敛的速度哈，因为在每步迭代中，忽略的项越大，那么带来的近似误差就越高了。

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

可以把左端项修改成梯度点积形式
$$\begin{array}{rl}\left(\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p^{\prime }{)}_f\right)\cdot {\mathbf{S}}_f& =\left((\nabla p^{\prime }{)}_f{\overline{{\mathbf{D}}_f^{\mathbf{v}}}}^T\right)\cdot {\mathbf{S}}_f\\ & =(\nabla p^{\prime }{)}_f\cdot \left({\overline{{\mathbf{D}}_f^{\mathbf{v}}}}^T{\mathbf{S}}_f\right)\\ & =(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f^{\prime }\end{array}$$
其中${\mathbf{S}}_f^{\prime }$的展开形式为
$${\mathbf{S}}_f^{\prime }={\overline{{\mathbf{D}}_f^{\mathbf{v}}}}^T{\mathbf{S}}_f=\left[\begin{array}{ccc}\overline{D_f^u}& 0& 0\\ 0& \overline{D_f^v}& 0\\ 0& 0& \overline{D_f^w}\end{array}\right]\left[\begin{array}{c}{S}_f^x\\ S_f^y\\ S_f^z\end{array}\right]=\left[\begin{array}{c}\overline{D_f^u}{S}_f^x\\ \overline{D_f^v}{S}_f^y\\ \overline{D_f^w}{S}_f^z\end{array}\right]$$
使用${\mathbf{S}}_f^{\prime }$，压力修正梯度项的离散过程和之前一样，得到
$$\begin{array}{rl}(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f^{\prime }& =(\nabla p^{\prime }{)}_f\cdot {\mathbf{E}}_f+(\nabla p^{\prime }{)}_f\cdot {\mathbf{T}}_f\\ & =\frac{E_f}{d_{CF}}(p_F^{\prime }-p_C^{\prime })+\underline{(\nabla p^{\prime }{)}_f\cdot {\mathbf{T}}_f}\end{array}$$
其中用到了对${\mathbf{S}}_f^{\prime }$如下分解
$${\mathbf{S}}_f^{\prime }={\mathbf{E}}_f+{\mathbf{T}}_f$$
分解的类型在第8章中已经讲过（三种类型：最小修正、正交修正、过松弛修正方法）。下划线项是由网格非正交特性所引起的，可以选择忽略或保留。如果忽略，其不会影响最终解，因为这只是修正项。如果保留，其将会在内循环中显式处理（OpenFOAM中的非正交循环）。由于每次迭代的求解都是从零压力修正场开始的，在求解方程时该项需要迭代更新。

忽略掉非正交项的作用，压力修正方程的线性化形式为
$$\begin{array}{rl}(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f^{\prime }& =\frac{E_f}{d_{CF}}(p_F^{\prime }-p_C^{\prime })\\ & ={\mathcal{D}}_f(p_F^{\prime }-p_C^{\prime })\end{array}$$
代回到压力修正方程$\sum _{f\sim nb(C)}\left(-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p^{\prime }{)}_f\cdot {\mathbf{S}}_f\right)=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}$中，可得其离散代数形式
$$a_C^{p^{\prime }}{p}_C^{\prime }+\sum _{F\sim NB(C)}{a}_F^{p^{\prime }}{p}_F^{\prime }=b_C^{p^{\prime }}$$
其中的系数为
$$\begin{array}{rl}{a}_F^{p^{\prime }}& =Flux{F}_f=-{\rho }_f{\mathcal{D}}_f\\ a_C^{p^{\prime }}& =-\sum _{f\sim nb(c)}Flux{F}_f=-\sum _{F\sim NB(C)}{a}_F^{p^{\prime }}\\ b_C^{p^{\prime }}& =-\sum _{f\sim nb(C)}Flux{V}_f+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}\\ & =-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}\end{array}$$
注意，对下划线项的不同近似将产生不同的算法，在最初的SIMPLE算法中是把这些项直接给刨掉了。

最后，上式中的质量流量${\dot{m}}_f^{\ast }$，是在求解完动量方程之后，基于最新的速度场（动量方程求得的速度场），采用惯常的Rhie-Chow插值技术算出的，即
$${\dot{m}}_f^{\ast }={\rho }_f{\mathbf{v}}_f^{\ast }\cdot {\mathbf{S}}_f={\rho }_f\overline{{\mathbf{v}}_f^{\ast }}\cdot {\mathbf{S}}_f-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\left(\nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}}\right)\cdot {\mathbf{S}}_f$$
基于算得的压力和速度修正场，单元形心处的压力和速度 以及 单元面处的质量流量都可以修正了。如前所述，在上上式中，下划线项在SIMPLE算法中的忽略处理 会 产生较大的压力修正值，其会减慢收敛速度或是造成发散。为了增加健壮性和提高收敛性，上上式所获得的压力修正值可做显式松弛处理。而更新速度场和质量流量场则无需欠松弛处理，因为压力修正是要来保证速度场与质量流量场的质量守恒的（如果再做欠松弛的话，就无法保证守恒特性了）。把欠松弛因子定义为${\lambda }^p$，使用如下修正方程
$$\begin{array}{rlrl}{\mathbf{v}}_C^{\ast \ast }& ={\mathbf{v}}_C^{\ast }+{\mathbf{v}}_C^{\prime }& {\mathbf{v}}_C^{\prime }& =-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C\\ {\dot{m}}_f^{\ast \ast }& ={\dot{m}}_f^{\ast }+{\dot{m}}_f^{\prime }& {\dot{m}}_f^{\prime }& =-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\nabla p_f^{\prime }\cdot {\mathbf{S}}_f\\ p_C^{\ast }& =p_C^{(n)}+{\lambda }^p{p}_C^{\prime }& & \end{array}$$

5.3 ${\mathcal{D}}_f$项的计算

上一小节讲到${\mathbf{S}}_f^{\prime }={\mathbf{E}}_f+{\mathbf{T}}_f$的分解方法可以是第8章中讲到的任何一个，那么不同的分解方法会产生不同的${\mathcal{D}}_f$表达式，如下。

5.3.1 最小修正方法 Minimum Correction Approach

该方法中，${\mathbf{E}}_f$与${\mathbf{T}}_f$正交，从而使得修正量最小，即
$${\mathbf{E}}_f=({\mathbf{e}}_{CF}\cdot {\mathbf{S}}_f^{\prime }){\mathbf{e}}_{CF}$$
其中${\mathbf{e}}_{CF}$是$CF$方向的单位矢量，结合上节${\mathbf{S}}_f^{\prime }$的表达式${\mathbf{S}}_f^{\prime }={\left[\begin{array}{ccc}\overline{D_f^u}{S}_f^x& \overline{D_f^v}{S}_f^y& \overline{D_f^w}{S}_f^z\end{array}\right]}^T$，算得${\mathbf{E}}_f$为
$${\mathbf{E}}_f=\frac{d_{CF}^x\overline{D_f^u}{S}_f^x+d_{CF}^y\overline{D_f^v}{S}_f^y+d_{CF}^z\overline{D_f^w}{S}_f^z}{d_{CF}^2}{\mathbf{d}}_{CF}$$
使用上式，推得${\mathcal{D}}_f$的表达式为
$${\mathcal{D}}_f=\frac{E_f}{d_{CF}}=\frac{d_{CF}^x\overline{D_f^u}{S}_f^x+d_{CF}^y\overline{D_f^v}{S}_f^y+d_{CF}^z\overline{D_f^w}{S}_f^z}{(d_{CF}^x{)}^2+(d_{CF}^y{)}^2+(d_{CF}^z{)}^2}$$

5.3.2 正交修正方法 Orthogonal Correction Approach

这个修正方法是让${\mathbf{E}}_f$和${\mathbf{S}}_f^{\prime }$的幅值相同，即
$${\mathbf{E}}_f={\mathbf{S}}_f^{\prime }{\mathbf{e}}_{CF}$$
这样一来，${\mathcal{D}}_f$的表达式为
$${\mathcal{D}}_f=\frac{E_f}{d_{CF}}=\sqrt{\frac{(\overline{D_f^u}{S}_f^x{)}^2+(\overline{D_f^v}{S}_f^y{)}^2+(\overline{D_f^w}{S}_f^z{)}^2}{(d_{CF}^x{)}^2+(d_{CF}^y{)}^2+(d_{CF}^z{)}^2}}$$

5.3.3 过松弛方法 Over-Relaxed Approach

这个方法是让${\mathbf{T}}_f$和和${\mathbf{S}}_f^{\prime }$垂直，那么${\mathbf{E}}_f$为
$${\mathbf{E}}_f=\frac{{\mathbf{S}}_f^{\prime }\cdot {\mathbf{S}}_f^{\prime }}{{\mathbf{d}}_{CF}\cdot {\mathbf{S}}_f^{\prime }}{\mathbf{d}}_{CF}$$
推得${\mathcal{D}}_f$的表达式为
$${\mathcal{D}}_f=\frac{E_f}{d_{CF}}=\frac{(\overline{D_f^u}{S}_f^x{)}^2+(\overline{D_f^v}{S}_f^y{)}^2+(\overline{D_f^w}{S}_f^z{)}^2}{d_{CF}^x\overline{D_f^u}{S}_f^x+d_{CF}^y\overline{D_f^v}{S}_f^y+d_{CF}^z\overline{D_f^w}{S}_f^z}$$
上述方法中的任何一个皆可用于计算${\mathcal{D}}_f$的值。

5.4 同位（网格）SIMPLE算法

将同位网格上SIMPLE算法的流程汇总如下：

1. 为计算时刻$t+\Delta t$的解，用时刻$t$的压力、速度和质量流量场$p_C^t$、${\mathbf{v}}_C^t$、${\dot{m}}_f^t$作为初始猜测值$p_C^{(n)}$、${\mathbf{v}}_C^{(n)}$、${\dot{m}}_f^{(n)}$；
2. 组装并求解动量方程以获得新速度场${\mathbf{v}}_C^{\ast }$；
   $${\mathbf{v}}_C^{\ast }+{\mathbf{H}}_C[{\mathbf{v}}^{\ast }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{(n)}{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$$
3. 基于新速度场${\mathbf{v}}_C^{\ast }$，使用Rhie-Chow插值来更新单元面处的质量流量，算得满足动量方程的质量流量场${\dot{m}}_f^{\ast }$；
   $${\dot{m}}_f^{\ast }={\rho }_f{\mathbf{v}}_f^{\ast }\cdot {\mathbf{S}}_f={\rho }_f\overline{{\mathbf{v}}_f^{\ast }}\cdot {\mathbf{S}}_f-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\left(\nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}}\right)\cdot {\mathbf{S}}_f$$
4. 使用新的质量流量${\dot{m}}_f^{\ast }$组装压力修正方程，并求解该方程以获得压力修正场$p_C^{\prime }$；
   $$a_C^{p^{\prime }}{p}_C^{\prime }+\sum _{F\sim NB(C)}{a}_F^{p^{\prime }}{p}_F^{\prime }=b_C^{p^{\prime }}$$
   其中
   $$\begin{array}{rl}{a}_F^{p^{\prime }}& =Flux{F}_f=-{\rho }_f{\mathcal{D}}_f\\ a_C^{p^{\prime }}& =-\sum _{f\sim nb(c)}Flux{F}_f=-\sum _{F\sim NB(C)}{a}_F^{p^{\prime }}\\ b_C^{p^{\prime }}& =-\sum _{f\sim nb(C)}Flux{V}_f+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}\\ & =-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}\end{array}$$
5. 基于压力修正场，更新单元形心处的压力和速度场，以及单元面处的质量流量场，以获得满足连续方程的场，这些场分别用$p_C^{\ast }$、${\mathbf{v}}_C^{\ast \ast }$、${\dot{m}}_f^{\ast \ast }$表示；
   $$\begin{array}{rlrl}{\mathbf{v}}_C^{\ast \ast }& ={\mathbf{v}}_C^{\ast }+{\mathbf{v}}_C^{\prime }& {\mathbf{v}}_C^{\prime }& =-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C\\ {\dot{m}}_f^{\ast \ast }& ={\dot{m}}_f^{\ast }+{\dot{m}}_f^{\prime }& {\dot{m}}_f^{\prime }& =-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\nabla p_f^{\prime }\cdot {\mathbf{S}}_f\\ p_C^{\ast }& =p_C^{(n)}+{\lambda }^p{p}_C^{\prime }& & \end{array}$$
6. 将所得解作为新的初始值，返回第2步，继续该流程直到收敛；
   $$p_C^{(n)}=p_C^{\ast },{\mathbf{v}}_C^{(n)}={\mathbf{v}}_C^{\ast \ast },{\dot{m}}_f^{(n)}={\dot{m}}_f^{\ast \ast }$$
7. 将所得收敛解作为$t+\Delta t$时刻的解；
   $$p_C^{t+\Delta t}=p_C^{(n)}=p_C^{\ast },{\mathbf{v}}_C^{t+\Delta t}={\mathbf{v}}_C^{(n)}={\mathbf{v}}_C^{\ast \ast },{\dot{m}}_C^{t+\Delta t}={\dot{m}}_f^{(n)}={\dot{m}}_f^{\ast \ast }$$
8. 推进到下一个时间步；
   $$t=t+\Delta t$$
9. 返回第1步并重复该流程直至到达最终时间步。

其流程图为：

此时需要一个小例子来让算法更容易理解。

---

例3

如图，二维问题如下所示，给定$p_W=100$、$p_N=20$、$p_E=50$，面$s$处的进口条件为$v_s=20$和零压力梯度。

流动是稳态的，且密度是均匀的，密度值为1。$u_C$和$v_C$的动量方程为（这里只是个简单的小例子，所以忽略了瞬态项、扩散项、对流项，只保留了压力梯度项）
$$\begin{array}{rl}{u}_C& =-d_x(p_e-p_w)\\ v_C& =-d_y(p_n-p_s)\end{array}$$
其中常数$d_x=1$和$d_y=0.25$，单元有$\Delta x=\Delta y=1$。使用同位SIMPLE算法，推导压力修正方程并求解方程以找到单元$C$的压力值。以$p_C^{(n)}=70$作为$C$处压力的初始猜测值。

解
进口的零压力梯度条件可用于计算
$$p_s=P_C=P_C^{(n)}=70$$
计算$u_C^{(n)}$和$v_C^{(n)}$
$$\begin{array}{rl}{u}_C^{\ast }& =-d_x(p_e^{(n)}-p_w^{(n)})=-1(60-85)=25\\ v_C^{\ast }& =-d_y(p_n^{(n)}-p_s^{(n)})=-0.25(45-70)=6.25\end{array}$$
计算中用到了$p_e^{(n)}=0.5(p_E^{(n)}+p_C^{(n)})$，$p_n^{(n)}=0.5(p_C^{(n)}+p_N^{(n)})$，$p_w^{(n)}=0.5(p_C^{(n)}+p_W^{(n)})$

现在要构建压力修正方程，通过把质量通量和其修正代入到质量守恒方程来实现。使用Rhie-Chow插值，算得面质量通量为
$$\begin{array}{rl}{\dot{m}}_e^{\ast }=& \rho u_e^{\ast }\Delta y=u_e^{\ast }\\ =& \overline{u_e^{\ast }}-d_x\underset{\begin{array}{c}pressuredifference\\ acrossface\end{array}}{\underbrace{[\left(p_E^{(n)}-p_C^{(n)}\right)}}-\underset{\begin{array}{c}averagepressuredifference\\ acrossface\end{array}}{\underbrace{0.5\left(p_{ee}^{(n)}-p_e^{(n)}+p_e^{(n)}-p_w^{(n)}\right)]}}\\ =& 0.5(u_C^{\ast }+u_E^{\ast })-d_x\left[\left(p_E^{(n)}-p_C^{(n)}\right)-0.5\left(p_{ee}^{(n)}-p_e^{(n)}+p_e^{(n)}-p_w^{(n)}\right)\right]\\ =& -0.5\left[d_x\left(p_e^{(n)}-p_w^{(n)}\right)+d_x\left(p_{ee}^{(n)}-p_e^{(n)}\right)\right]\\ & -d_x\left[\left(p_E^{(n)}-p_C^{(n)}\right)-0.5\left(p_{ee}^{(n)}-p_e^{(n)}+p_e^{(n)}-p_w^{(n)}\right)\right]\\ =& -d_x\left(p_E^{(n)}-p_C^{(n)}\right)\\ =& -1(50-70)\\ =& 20\end{array}$$
同样可以算得$n$界面的质量通量
$$\begin{array}{rl}{\dot{m}}_n^{\ast }=& \rho v_n^{\ast }\Delta x=v_n^{\ast }\\ =& \overline{v_n^{\ast }}-d_y\left[\left(p_N^{(n)}-p_C^{(n)}\right)-0.5\left(p_{nn}^{(n)}-p_n^{(n)}+p_n^{(n)}-p_s^{(n)}\right)\right]\\ =& 0.5(v_C^{\ast }+v_N^{\ast })-d_y\left[\left(p_N^{(n)}-p_C^{(n)}\right)-0.5\left(p_{nn}^{(n)}-p_n^{(n)}+p_n^{(n)}-p_s^{(n)}\right)\right]\\ =& -0.5\left[d_y\left(p_n^{(n)}-p_s^{(n)}\right)+d_y\left(p_{nn}^{(n)}-p_n^{(n)}\right)\right]\\ & -d_y\left[\left(p_N^{(n)}-p_C^{(n)}\right)-0.5\left(p_{nn}^{(n)}-p_n^{(n)}+p_n^{(n)}-p_s^{(n)}\right)\right]\\ =& -d_y\left(p_N^{(n)}-p_C^{(n)}\right)\\ =& -0.25(20-70)\\ =& -12.5\end{array}$$
还有$w$界面的质量通量
$$\begin{array}{rl}{\dot{m}}_w^{\ast }=& -\rho u_w^{\ast }\Delta y=-u_w^{\ast }\\ =& -\overline{u_w^{\ast }}+d_x\left[\left(p_C^{(n)}-p_W^{(n)}\right)-0.5\left(p_w^{(n)}-p_{ww}^{(n)}+p_e^{(n)}-p_w^{(n)}\right)\right]\\ =& -0.5(u_C^{\ast }+u_W^{\ast })+d_x\left[\left(p_C^{(n)}-p_W^{(n)}\right)-0.5\left(p_w^{(n)}-p_{ww}^{(n)}+p_e^{(n)}-p_w^{(n)}\right)\right]\\ =& -0.5\left[-d_x\left(p_e^{(n)}-p_w^{(n)}\right)-d_x\left(p_w^{(n)}-p_{ww}^{(n)}\right)\right]\\ & +d_x\left[\left(p_C^{(n)}-p_W^{(n)}\right)-0.5\left(p_w^{(n)}-p_{ww}^{(n)}+p_e^{(n)}-p_w^{(n)}\right)\right]\\ =& d_x\left(p_C^{(n)}-p_W^{(n)}\right)\\ =& 1(70-100)\\ =& -30\end{array}$$
在速度进口面$s$处的质量流量为
$${\dot{m}}_s=-\rho v_s\Delta x=-1\ast 20\ast 1=-20$$
压力修正方程将从连续方程中构建，即
$$({\dot{m}}_e^{\ast }+{\dot{m}}_e^{\prime })+({\dot{m}}_n^{\ast }+{\dot{m}}_n^{\prime })+({\dot{m}}_w^{\ast }+{\dot{m}}_w^{\prime })+{\dot{m}}_s=0$$
有
$$\begin{array}{rl}{\dot{m}}_e^{\prime }+{\dot{m}}_n^{\prime }+{\dot{m}}_w^{\prime }=& -({\dot{m}}_e^{\ast }+{\dot{m}}_n^{\ast }+{\dot{m}}_w^{\ast }+{\dot{m}}_s)\\ =& -(20+12.5-30-20)\\ =& 17.5\end{array}$$
根据前面推导，有
$$\begin{array}{rl}{\dot{m}}_e^{\prime }& =-d_x(p_E^{\prime }-p_C^{\prime })\\ {\dot{m}}_n^{\prime }& =-d_y(p_N^{\prime }-p_C^{\prime })\\ {\dot{m}}_w^{\prime }& =d_x(p_C^{\prime }-p_W^{\prime })\end{array}$$
于是有
$$-(p_E^{\prime }-p_C^{\prime })-0.25(p_N^{\prime }-p_C^{\prime })+(p_C^{\prime }-p_W^{\prime })=17.5$$
即
$$2.25p_C^{\prime }-p_E^{\prime }-0.25p_N^{\prime }-p_W^{\prime }=17.5$$
由于$E、W、N$的值都为精确解，所以其修正量为零，故
$$2.25p_C^{\prime }=17.5\Rightarrow p_C^{\prime }=\frac{70}{9}\approx 7.78$$
基于该修正量，接下来修正单元$C$形心处的压力和速度分量
$$\begin{array}{rl}{p}_C^{\ast }& =p_C^{(n)}+p_C^{\prime }=70+7.78=77.78\\ u_C^{\prime }& =-d_x(p_e^{\prime }-p_w^{\prime })=-d_x[0.5(p_E^{\prime }+p_C^{\prime })-0.5(p_W^{\prime }+p_C^{\prime })]=0\\ u_C^{\ast \ast }& =u_C^{\ast }+u_C^{\prime }=u_C^{\ast }=25\\ v_C^{\prime }& =-d_y(p_n^{\prime }-p_s^{\prime })=-d_y[0.5(p_N^{\prime }+p_C^{\prime })-p_C^{\prime }]=-d_y[-0.5p_C^{\prime }]\\ & =-0.25\ast (-0.5)\ast 7.78=0.9725\\ v_C^{\ast \ast }& =v_C^{\ast }+v_C^{\prime }=6.25+0.9725=7.2225\end{array}$$
由于这个例子较为简单，动量方程中并不包含非线性项，所以迭代一步即可获得收敛解，无需再进行后续迭代了。

---