有限体积法（12）——SIMPLE算法

承接上篇《交错网格》，本文介绍流动方程在交错网格上的离散以及SIMPLE算法。

方程离散

二维稳态的动量方程和连续性方程如下：
$$\begin{array}{rl}\frac{\partial }{\partial x}(\rho uu)+\frac{\partial }{\partial y}(\rho vu)& =\frac{\partial }{\partial x}\left(\mu \frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial u}{\partial y}\right)-\frac{\partial p}{\partial x}+S_u\end{array}\text{\hspace{1em}(1-a)}$$

$$\begin{array}{rl}\frac{\partial }{\partial x}(\rho uv)+\frac{\partial }{\partial y}(\rho vv)& =\frac{\partial }{\partial x}\left(\mu \frac{\partial v}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial v}{\partial y}\right)-\frac{\partial p}{\partial y}+S_v\end{array}\text{\hspace{1em}(1-b)}$$

$$\begin{array}{rl}\frac{\partial }{\partial x}(\rho u)+\frac{\partial }{\partial y}(\rho v)& =0\end{array}\text{\hspace{1em}(1-c)}$$
交错网格的编号如下图：

图1 交错网格

$u$动量方程离散

$u$ 动量方程在 $u$ 网格的离散方程可以写成如下形式：
$$\begin{array}{r}{a}_{i,J}{u}_{i,J}=\Sigma a_{nb}{u}_{nb}-\frac{p_{I,J}-p_{I-1,J}}{\delta x_u}\Delta V_u+\bar{S}\Delta V_u\end{array}\text{\hspace{1em}(2-a)}$$

或
$$\begin{array}{r}{a}_{i,J}{u}_{i,J}=\Sigma a_{nb}{u}_{nb}-(p_{I,J}-p_{I-1,J})A_{i,J}+b_{i,J}\end{array}\text{\hspace{1em}(2-b)}$$

式中，$\Delta V_u$是$u$网格单元的体积，$b_{i,J}=\bar{S}\Delta V_u$是动量方程的源项，$A_{i,J}$是$u$网格单元的边界面积（$w$边界或$e$边界）。压力梯度项使用了中心差分格式来近似，并使用了$u$网格的边界压力值，即压力网格的节点值。
对于二维流场，方程的第二项$\Sigma a_{nb}{b}_{nb}$包含了$E、W、N$和$S$四个方向上相邻单元的速度值，其网格下标分别为$(i+1,J)、(i-1,J)、(i,J+1)$和$(i,J-1)$。如下图所示，

图2 u控制体单元及相邻节点速度分量

公式中的系数$a_{i,J}$和$a_{nb}$的计算对应不同的差分格式有不同的公式，事实上无论选用何种差分格式，方程系数$a_{i,J}$和$a_{nb}$都是控制体单元边界处的单位面积对流量$F(=\rho u)$和单位面积扩散量$D(=\Gamma /\delta x)$的组合。在交错网格（二维均匀网格）的编号系统下，$u$网格单元$(i,J)$的$e、w、n$和$s$边界处的$F$值和$D$值的计算公式为
$$\begin{array}{rl}{F}_w& =(\rho u{)}_w=\frac{F_{i,J}+F_{i-1,J}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I-1,J}}{2}\right)u_{i,J}+\left(\frac{{\rho }_{I-1,J}+{\rho }_{I-2,J}}{2}\right)u_{i-1,J}\right]\end{array}\text{\hspace{1em}(3-a)}$$

$$\begin{array}{rl}{F}_e& =(\rho u{)}_e=\frac{F_{i+1,J}+F_{i,J}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I+1,J}+{\rho }_{I,J}}{2}\right)u_{i+1,J}+\left(\frac{{\rho }_{I,J}+{\rho }_{I-1,J}}{2}\right)u_{i,J}\right]\end{array}\text{\hspace{1em}(3-b)}$$

$$\begin{array}{rl}{F}_s& =(\rho v{)}_s=\frac{F_{I,j}+F_{I-1,j}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I,J-1}}{2}\right)v_{I,j}+\left(\frac{{\rho }_{I-1,J}+{\rho }_{I-1,J-1}}{2}\right)v_{I-1,j}\right]\end{array}\text{\hspace{1em}(3-c)}$$

$$\begin{array}{rl}{F}_n& =(\rho v{)}_n=\frac{F_{I,j+1}+F_{I-1,j+1}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I,J+1}}{2}\right)v_{I,j+1}+\left(\frac{{\rho }_{I-1,J}+{\rho }_{I-1,J+1}}{2}\right)v_{I-1,j+1}\right]\end{array}\text{\hspace{1em}(3-d)}$$
$$\begin{array}{r}{D}_w=\frac{{\Gamma }_{I-1,J}}{x_i-x_{i-1}}\end{array}\text{\hspace{1em}(3-e)}$$

$$\begin{array}{r}{D}_e=\frac{{\Gamma }_{I,J}}{x_{i+1}-x_i}\end{array}\text{\hspace{1em}(3-f)}$$

$$\begin{array}{r}{D}_s=\frac{{\Gamma }_{i,j}}{y_J-y_{J-1}}=\frac{{\Gamma }_{I-1,J}+{\Gamma }_{I-1,J-1}+{\Gamma }_{I,J-1}+{\Gamma }_{I,J}}{4(y_J-y_{J-1})}\end{array}\text{\hspace{1em}(3-g)}$$
$$\begin{array}{r}{D}_n=\frac{{\Gamma }_{i,j+1}}{y_{J+1}-y_J}=\frac{{\Gamma }_{I-1,J+1}+{\Gamma }_{I-1,J}+{\Gamma }_{I,J}+{\Gamma }_{I,J+1}}{4(y_{J+1}-y_J)}\end{array}\text{\hspace{1em}(3-h)}$$
在交错网格中，密度$\rho$和扩散系数$\Gamma$和压力一样都是标量，因此这些标量值都保存在主网格节点（图中的黑色圆点）中，公式$(4)-(7)$中的边界密度值和公式$(10)$和$(11)$中的扩散系数都周围节点值的平均值来近似。

$v$动量方程离散

同理，$v$动量方程在$v$网格上的离散方程可以写为如下形式：
$$\begin{array}{r}{a}_{I,j}{v}_{I,j}=\Sigma a_{nb}{v}_{nb}+(p_{I,J-1}-p_{I,J})A_{I,j}+b_{I,j}\end{array}\text{\hspace{1em}(4)}$$
相邻单元的网格编号如下图

图3 v控制体单元及相邻节点速度分量

同样，$v$离散方程中的$a_{I,j}$和$a_{nb}$也是$F$和$D$的组合，其在$v$网格上的计算公式如下
$$\begin{array}{rl}{F}_w& =(\rho u{)}_w=\frac{F_{i,J}+F_{i,J-1}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I-1,J}}{2}\right)u_{i,J}+\left(\frac{{\rho }_{I-1,J-1}+{\rho }_{I,J-1}}{2}\right)u_{i,J-1}\right]\end{array}\text{\hspace{1em}(5-a)}$$

$$\begin{array}{rl}{F}_e& =(\rho u{)}_e=\frac{F_{i+1,J}+F_{i+1,J-1}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I+1,J}}{2}\right)u_{i+1,J}+\left(\frac{{\rho }_{I,J-1}+{\rho }_{I+1,J-1}}{2}\right)u_{i+1,J-1}\right]\end{array}\text{\hspace{1em}(5-b)}$$

$$\begin{array}{rl}{F}_s& =(\rho v{)}_s=\frac{F_{I,j}+F_{I,j-1}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I,J-1}}{2}\right)v_{I,j}+\left(\frac{{\rho }_{I,J-1}+{\rho }_{I,J-2}}{2}\right)v_{I,j-1}\right]\end{array}\text{\hspace{1em}(5-c)}$$

$$\begin{array}{rl}{F}_n& =(\rho v{)}_n=\frac{F_{I,j+1}+F_{I,j}}{2}\\ & =\frac{1}{2}\left[\left(\frac{{\rho }_{I,J}+{\rho }_{I,J+1}}{2}\right)v_{I,j+1}+\left(\frac{{\rho }_{I,J-1}+{\rho }_{I,J}}{2}\right)v_{I,j}\right]\end{array}\text{\hspace{1em}(5-d)}$$
$$\begin{array}{r}{D}_w=\frac{{\Gamma }_{I-1,J}+{\Gamma }_{I-1,J-1}+{\Gamma }_{I,J-1}+{\Gamma }_{I,J}}{4(x_I-x_{I-1})}\end{array}\text{\hspace{1em}(5-e)}$$

$$\begin{array}{r}{D}_e=\frac{{\Gamma }_{I+1,J}+{\Gamma }_{I+1,J-1}+{\Gamma }_{I,J-1}+{\Gamma }_{I,J}}{4(x_{I+1}-x_I)}\end{array}\text{\hspace{1em}(5-f)}$$

$$\begin{array}{r}{D}_s=\frac{{\Gamma }_{I,J-1}}{y_j-y_{j-1}}\end{array}\text{\hspace{1em}(5-g)}$$

$$\begin{array}{r}{D}_n=\frac{{\Gamma }_{I,J}}{y_{j+1}-y_j}\end{array}\text{\hspace{1em}(5-h)}$$

连续性方程离散

连续性方程是在主网格上进行离散的，即在标量$p、\rho 、\Gamma$的网格上，网格单元及相邻节点如图1的$P$单元所示。连续性方程的离散过程与动量方程类似，其离散后的方程为
$$\begin{array}{r}[(\rho uA{)}_{i+1,J}-(\rho uA{)}_{i,J}]+[(\rho vA{)}_{I,j+1}-(\rho vA{)}_{I,j}]=0\end{array}\text{\hspace{1em}(6-a)}$$

或者写成
$$\begin{array}{r}[F_e{A}_{i+1,J}-F_w{A}_{i,J}]+[F_n{A}_{I,j+1}-F_s{A}_{I,j}]=0\end{array}\text{\hspace{1em}(6-b)}$$
边界处$F$值计算公式如下
$$\begin{array}{r}{F}_e=(\rho u{)}_e=\left(\frac{{\rho }_{I,J}+{\rho }_{I+1,J}}{2}\right)u_{i+1,J}\end{array}\text{\hspace{1em}(7-a)}$$

$$\begin{array}{r}{F}_w=(\rho u{)}_w=\left(\frac{{\rho }_{I,J}+{\rho }_{I-1,J}}{2}\right)u_{i,J}\end{array}\text{\hspace{1em}(7-b)}$$

$$\begin{array}{r}{F}_n=(\rho v{)}_n=\left(\frac{{\rho }_{I,J}+{\rho }_{I,J+1}}{2}\right)v_{I,j+1}\end{array}\text{\hspace{1em}(7-c)}$$

$$\begin{array}{r}{F}_s=(\rho v{)}_s=\left(\frac{{\rho }_{I,J}+{\rho }_{I,J-1}}{2}\right)v_{I,j}\end{array}\text{\hspace{1em}(7-d)}$$

现在我们得到了交错网格下二维压力速度耦合问题的有限体积法离散方程组，即公式$(2)、(4)$和$(6)$的组合，有
$$\{\begin{array}{r}{a}_{i,J}{u}_{i,J}=\Sigma a_{nb}{u}_{nb}-(p_{I,J}-p_{I-1,J})A_{i,J}+b_{i,J}\\ \\ a_{I,j}{v}_{I,j}=\Sigma a_{nb}{v}_{nb}+(p_{I,J-1}-p_{I,J})A_{I,j}+b_{I,j}\\ \\ [(\rho uA{)}_{i+1,J}-(\rho uA{)}_{i,J}]+[(\rho vA{)}_{I,j+1}-(\rho vA{)}_{I,j}]=0\end{array}\text{\hspace{1em}(8)}$$

两个关键问题

为降低问题的复杂度，我们先考虑不可压缩流体的流动情况，即方程组$(8)$中的密度$\rho$是已知量。省去求解密度后，需要求解的变量还有三个，速度$u、v$和$p$，然而要求解这个离散方程组我们还需要解决两个关键问题：

问题1

方程组中，未知量$u、v$和$p$是耦合的，且动量方程是非线性的。
非线性问题的解决方法是通过假设系数中的速度场$u$和$v$是已知的，给定一组初始值然后进行迭代。
虽然方程组是耦合的，但理论上是可以直接求解的，如果采用同时求解$u、v、p$的直接方法，计算过程大概如下：

1. 先假设一组流场的初始值，由流场初始值计算出系数的初始值；
2. 在给定的系数下，方程组联合求解各节点上$u、v$和$p$的值；
3. 根据所得的$u、v$和$p$的新值更新系数，再在新系数下计算各节点的变量值，如此反复迭代，直至计算收敛。

这种在全计算域范围内同时求解$u、v、p$的方法（耦合求解法）对计算资源的要求较高，在解决工程问题中不够经济，因此实际中大多采用所谓的分离式求解法（顺序求解法），即将$u、v、p$分别独立求解：

• 求解$u$速度场时认为$v$和$p$的分布是已知的；
• 求解$v$场时，认为$u$和$p$的分布是已知的；
• 求解$p$场时认为$u$和$v$的分布是已知的。
  分离式求解法的计算过程与耦合法大致类似，不同的是第二步。耦合法是联合求解方程组，同时计算出$u、v、p$，每个方程中除了系数是用的初始值或上一步迭代结果外，其他$u、v、p$都是要求解的量；而分离法是在计算$u$时，系数和所有的$v、p$都认为是已知量（初始值或上一步迭代值），然后求解$u$，计算$v$时也类似。

问题2

虽然分离求解法可以解决流动方程组耦合的问题，但从动量方程可见，速度场的变化会影响到压力场的分布，而压力场的变化反过来也会影响速度场的分布。那么在进行求解过程中，动量方程可以根据速度场和压力场的初始值更新速度场，速度场得到了更新，那下一步就应该是根据更新后的速度场来更新压力场。然而，方程组$(8)$中没有关于描述压力$p$的独立方程，压力是作为动量源项出现在动量方程中的，而压力和速度的耦合关系隐含在连续性方程中。因此，压力场如何单独求解，如何根据当前更新后的速度场来更新压力场就是方程组求解过程中的一个关键问题。解决这个问题的办法就是SIMPLE算法。

SIMPLE算法

SIMPLE的全称是“Semi-Implicit Method for Pressure-Linked Equations”，即求解压力耦合方程的半隐式方法，“Semi-Implicit”翻译过来就是半隐式，其指代含义后面会提到。SIMPLE算法由Patankar和Spalding（1972）提出来，是在交错网格上计算压力的一个guess-and-correct的过程。下面以笛卡尔坐标系下二维稳态不可压流体的层流流动为例来介绍SIMPLE算法的具体推导过程。
SIMPLE算法首先需要假设一个压力场（$p^{\ast }$）和一个流场（$u$和$v$），然后动量方程就可以在初始压力场和初始流场上计算速度分量$u$和$v$：
$$\begin{array}{r}{a}_{i,J}{u}_{i,J}^{\ast }=\Sigma a_{nb}{u}_{nb}^{\ast }+(p_{I-1,J}^{\ast }-p_{I,J}^{\ast })A_{i,J}+b_{i,J}\\ \\ a_{I,j}{v}_{I,j}^{\ast }=\Sigma a_{nb}{v}_{nb}^{\ast }+(p_{I,J-1}^{\ast }-p_{I,J}^{\ast })A_{I,j}+b_{I,j}\end{array}\}\text{\hspace{1em}(9)}$$
上式中，等号右边的$u^{\ast }、v^{\ast }、p^{\ast }$都是假设的初始值，等号左边的$u^{\ast }、v^{\ast }$是计算的初始速度分量，需要指出的是，系数$a_{i,J}、a_{I,j}、a_{nb}$的计算公式中也有$u^{\ast }$和$v^{\ast }$，也是使用假设值计算。然而，一般情况下，这样计算得到的速度分量$u$和$v$是不满足连续性方程的，因为假设的压力和速度分布一般是与真实情况是不符的。因此，为了得到合理的计算结果，就需要对压力场$p^{\ast }$和速度场$u^{\ast }、v^{\ast }$进行修正。定义压力的修正量和修正后的正确性分别为$p^{\prime }$和$p$，速度的修正量和修正后的正确值分别为$u^{\prime }、v^{\prime }$和$u、v$，则其之间的关系为
$$\begin{array}{r}p=p^{\ast }+p^{\prime }\end{array}\text{\hspace{1em}(10-a)}$$

$$\begin{array}{r}u=u^{\ast }+u^{\prime }\end{array}\text{\hspace{1em}(10-b)}$$

$$\begin{array}{r}v=v^{\ast }+v^{\prime }\end{array}\text{\hspace{1em}(10-c)}$$
将压力场的正确值$p$带入到动量方程$(2-b)$和$(4)$就能计算出正确的速度场$(u,v)$，如果用方程$(2-b)$和$(4)$减去公式$(9)$的方程，则得到
$$\begin{array}{r}{a}_{i,J}(u_{i,J}-u_{i,J}^{\ast })=\Sigma a_{nb}(u_{nb}-u_{nb}^{\ast })+[(p_{I-1,J}-p_{I-1,J}^{\ast })-(p_{I,J}-p_{I,J}^{\ast })]A_{i,J}\end{array}\text{\hspace{1em}(11-a)}$$

$$\begin{array}{r}{a}_{I,j}(v_{I,j}-v_{I,j}^{\ast })=\Sigma a_{nb}(v_{nb}-v_{nb}^{\ast })+[(p_{I,J-1}-p_{I,J-1}^{\ast })-(p_{I,J}-p_{I,J}^{\ast })]A_{I,j}\end{array}\text{\hspace{1em}(11-a)}$$

注：这里其实做了假设，即假设离散动量方程是线性的，实际上系数中都是包含速度的。

把公式$(10)$的关系带入到公式$(11)$中，将得到修正值的方程：
$$\begin{array}{r}{a}_{i,J}{u}_{i,J}^{\prime }=\Sigma a_{nb}{u}_{nb}^{\prime }+(p_{I-1,J}^{\prime }-p_{I,J}^{\prime })A_{i,J}\end{array}\text{\hspace{1em}(12-a)}$$

$$\begin{array}{r}{a}_{I,j}{v}_{I,j}^{\prime }=\Sigma a_{nb}{v}_{nb}^{\prime }+(p_{I,J-1}^{\prime }-p_{I,J}^{\prime })A_{I,j}\end{array}\text{\hspace{1em}(12-b)}$$
这里引入一个近似：将$\Sigma a_{nb}{u}_{nb}^{\prime }$和$\Sigma a_{nb}{v}_{nb}^{\prime }$两项省略掉。这两项的省略是SIMPLE算法中主要的近似操作。然后，我们就得到
$$\begin{array}{r}{u}_{i,J}^{\prime }=d_{i,J}(p_{I-1,J}^{\prime }-p_{I,J}^{\prime })\end{array}\text{\hspace{1em}(13-a)}$$

$$\begin{array}{r}{v}_{I,j}^{\prime }=d_{I,j}(p_{I,J-1}^{\prime }-p_{I,J}^{\prime })\end{array}\text{\hspace{1em}(13-b)}$$
其中，$d_{i,J}=\frac{A_{i,J}}{a_{i,J}}$，$d_{I,j}=\frac{A_{I,j}}{a_{I,j}}$.

把公式$(13-a)$和$(13-b)$代入到公式$(10-b)$和$(10-c)$中，有
$$\begin{array}{r}{u}_{i,J}=u_{i,J}^{\ast }+d_{i,J}(p_{I-1,J}^{\prime }-p_{I,J}^{\prime })\end{array}\text{\hspace{1em}(14-a)}$$

$$\begin{array}{r}{v}_{I,j}=v_{I,j}^{\ast }+d_{I,j}(p_{I,J-1}^{\prime }-p_{I,J}^{\prime })\end{array}\text{\hspace{1em}(14-b)}$$
同样对于$u_{i+1,J}$和$v_{I,j+1}$有
$$\begin{array}{r}{u}_{i+1,J}=u_{i+1,J}^{\ast }+d_{i+1,J}(p_{I,J}^{\prime }-p_{I+1,J}^{\prime })\end{array}\text{\hspace{1em}(15-a)}$$

$$\begin{array}{r}{v}_{I,j+1}=v_{I,j+1}^{\ast }+d_{I,j+1}(p_{I,J}^{\prime }-p_{I,J+1}^{\prime })\end{array}\text{\hspace{1em}(15-b)}$$
其中，$d_{i+1,J}=\frac{A_{i+1,J}}{a_{i+1,J}}$，$d_{I,j+1}=\frac{A_{I,j+1}}{a_{I,j+1}}$.

到此为止，我们只考虑了动量方程，如上所述，如果压力修正$p^{\prime }$已知，速度场就可以更新了，那么现在就剩下最后一个问题：如何利用连续性方程来计算压力修正$p^{\prime }$。
把更新的速度值代入到离散的连续性方程
$$[(\rho uA{)}_{i+1,J}-(\rho uA{)}_{i,J}]+[(\rho vA{)}_{I,j+1}-(\rho vA{)}_{I,j}]=0\text{\hspace{1em}(16)}$$

按上图的网格标号，离散的连续性方程则成为
$$\begin{array}{rl}[{\rho }_{i+1,J}& A_{i+1,J}(u_{i+1,J}^{\ast }+d_{i+1,J}(p_{I,J}^{\prime }-p_{I+1,J}^{\prime }))-{\rho }_{i,J}{A}_{i,J}(u_{i,J}^{\ast }+\\ & \\ & d_{i,J}(p_{I-1,J}^{\prime }-p_{I,J}^{\prime }))]+[{\rho }_{I,j+1}{A}_{I,j+1}(v_{I,j+1}^{\ast }+d_{I,j+1}(p_{I,J}^{\prime }-p_{I,J+1}^{\prime }))\\ & \\ & -{\rho }_{I,j}{A}_{I,j}(v_{I,j}^{\ast }+d_{I,j}(p_{I,J-1}^{\prime }-p_{I,J}^{\prime }))]=0\end{array}\text{\hspace{1em}(17)}$$
把$p^{\prime }$看作未知变量，重新整理一下，就变成
$$\begin{array}{rl}[(\rho dA{)}_{i+1,J}& +(\rho dA{)}_{i,J}+(\rho dA{)}_{I,j+1}+(\rho dA{)}_{I,J}]p_{I,J}^{\prime }=\\ & \\ & (\rho dA{)}_{i+1,J}{p}_{I+1,J}^{\prime }+(\rho dA{)}_{i,J}{p}_{I-1,J}^{\prime }+(\rho dA{)}_{I,J+1}{p}_{I,J+1}^{\prime }+(\rho dA{)}_{I,j}{p}_{I,J-1}^{\prime }\\ & \\ & +[(\rho u^{\ast }A{)}_{i,J}-(\rho u^{\ast }A{)}_{i+1,J}+(\rho v^{\ast }A{)}_{I,j}-(\rho v^{\ast }A{)}_{I,j+1}]\end{array}\text{\hspace{1em}(18)}$$
写出熟悉的样子就是
$$\begin{array}{r}{a}_{I,J}{p}_{I,J}^{\prime }=a_{I+1,J}{p}_{I+1,J}^{\prime }+a_{I-1,J}{p}_{I-1,J}^{\prime }+a_{I,J+1}{p}_{I,J+1}^{\prime }+a_{I,J-1}{p}_{I,J-1}^{\prime }+b_{I,J}^{\prime }\end{array}\text{\hspace{1em}(19)}$$
其中系数为

$a_{I+1,J}$	$(\rho dA{)}_{i+1,J}$
$a_{I-1,J}$	$(\rho dA{)}_{i,J}$
$a_{I,J+1}$	$(\rho dA{)}_{I,j+1}$
$a_{I,J-1}$	$(\rho dA{)}_{I,j}$
$b_{I,J}$	$(\rho u^{\ast }A{)}_{i,J}-(\rho u^{\ast }A{)}_{i+1,J}+(\rho v^{\ast }A{)}_{I,j}-(\rho v^{\ast }A{)}_{I,j+1}$
$a_{I,J}$	$a_{I+1,J}+a_{I-1,J}+a_{I,J+1}+a_{I,J-1}$

方程式$(19)$是用压力修正$p^{\prime }$描述的离散连续性方程，其中源项$b^{\prime }$代表的是由于不准确的$u^{\ast }、v^{\ast }$导致的不守恒流量，多次迭代修正后$b^{\prime }$会趋于零，因此$b^{\prime }$可以作为判断迭代收敛过程是否满足要求的判据。
根据假设的初始速度场，通过求解压力修正方程$(19)$可以得出计算域所有网格节点的压力修正$p^{\prime }$，一旦压力修正值知道了，那么压力场和速度场就可以根据公式$(14)$进一步更新了，然后按照这个过程不断迭代计算，压力修正值和速度修正值会越来越小，最终速度场和压力场会收敛于真实值。最后，SIMPLE算法的计算流程如下图所示。

关于SIMPLE算法的几点说明

【1】 公式$(12)$到$(13)$的过程引入了近似处理，即省略掉了$\Sigma a_{nb}{u}_{nb}^{\prime }$和$\Sigma a_{nb}{v}_{nb}^{\prime }$。这一近似处理并不会影响最终的计算结果，因为压力修正和速度修正最终都趋向于零，即最终$p^{\ast }=p,u^{\ast }=u,v^{\ast }=v$。

【2】 从物理意义上讲，$(p_{I-1,J}^{\prime }-p_{I,J}^{\prime })$代表了压力修正对$u_{i,J}^{\prime }$的直接影响，而$\Sigma a_{nb}{u}_{nb}^{\prime }$则反映了压力修正对$u_{i,J}^{\prime }$的间接的或隐含的影响（压力修正通过影响周围节点的$u_{nb}^{\prime }$进而影响中心节点的$u_{i,J}^{\prime }$）。去掉了$\Sigma a_{nb}{u}_{nb}^{\prime }$这一项就称为半隐（Semi-Implicit），若保留这一部分，$u^{\prime }$方程就是一个全隐的代数方程，即网格点上的$u^{\prime }$就是耦合的，必须同时计算出，不像SIMPLE算法中那样可以在每个节点上独立计算。

【3】 一般情况下，如果按公式$(10)$直接用压力修正值来更新压力场往往会导致迭代发散。使用亚松弛迭代方法可以解决这一问题，压力修正方程就成为
$$\begin{array}{r}{p}^{new}=p^{\ast }+{\alpha }_p{p}^{\prime }\end{array}\text{\hspace{1em}(20)}$$

式中，${\alpha }_p$是亚松弛因子。因为修正量随着迭代计算会区域零，所以亚松弛因子的大小并不会影响最终计算的结果。
当${\alpha }_p=1$时，实际上就是直接用压力修正值$p^{\prime }$来修正压力初始值$p^{\ast }$。然而，当压力初始值与真实值相差较大时，这样计算出来的压力修正值$p^{\prime }$会是很大的数，如果直接用$p^{\prime }$来修正$p^{\ast }$，即使$p^{\prime }$的修正方向是正确的，那也会导致矫枉过正，造成计算的不稳定。
当${\alpha }_p$取0~1之间的值时，即人为控制压力修正的幅度。每次迭代压缩了修正的幅度，控制了压力变动的幅度，从而保证计算的稳定性。然而，${\alpha }_p$的值也并不是越小越好，过小的松弛因子虽然可以保证收敛稳定，但要迭代到最后的收敛解就需要更多次的迭代计算，这样会导致收敛速度很慢。计算中所希望的${\alpha }_p$值是尽可能大（以加快收敛）但又不引起计算发散。
速度的修正也需要采用亚松弛迭代，
$$\begin{array}{r}{u}^{new}=u^{(n-1)}+{\alpha }_u{u}^{\prime }=u^{(n-1)}+{\alpha }_u(u^{(n)}-u^{(n-1)})={\alpha }_u{u}^{(n)}+(1-{\alpha }_u)u^{(n-1)}\end{array}\text{\hspace{1em}(21-a)}$$
$$\begin{array}{r}{v}^{new}=v^{(n-1)}+{\alpha }_v{v}^{\prime }=v^{(n-1)}+{\alpha }_v(v^{(n)}-v^{(n-1)})={\alpha }_v{v}^{(n)}+(1-{\alpha }_v)v^{(n-1)}\end{array}\text{\hspace{1em}(21-b)}$$
其中，${\alpha }_u$和${\alpha }_v$是速度$u$和$v$修正的亚松弛因子，取值也是在$0$~$1$之间，$u^{(n-1)}$和$v^{(n-1)}$是上一次迭代计算的结果，$u^{(n)}$和$v^{(n)}$是本次迭代以$p^{\prime }$直接修正后的值。

【4】 最后需要指出的是，虽然亚松弛因子存在一个最优值，使得计算能够保证收敛且收敛最快，但并没有一个通用的最优亚松弛因子，不同计算的最优值是不同的。此外，在计算之前也没有很有效的方法来确定松弛因子的最优值，只能是在计算中试验性地取值，反复调整最终确定松弛因子的最优值。

参考资料

1. Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
2. 陶文铨. 数值传热学(第2版)[M]. 西安交通大学出版社, 2001.
3. 李人宪. 有限体积法基础-第2版[M]. 国防工业出版社, 2008.