有限体积法（12）—

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

方程离散

二维稳态的动量方程和连续性方程如下：
∂∂x(ρuu)+∂∂y(ρvu)=∂∂x(μ∂u∂x)+∂∂y(μ∂u∂y)−∂p∂x+Su(1-a)\begin{aligned} \frac{\partial}{\partial x} (\rho u u) + \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{aligned} \tag{1-a}∂x∂(ρuu)+∂y∂(ρvu)=∂x∂(μ∂x∂u)+∂y∂(μ∂y∂u)−∂x∂p+Su(1-a)

∂∂x(ρuv)+∂∂y(ρvv)=∂∂x(μ∂v∂x)+∂∂y(μ∂v∂y)−∂p∂y+Sv(1-b)\begin{aligned} \frac{\partial}{\partial x} (\rho u v) + \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{aligned} \tag{1-b}∂x∂(ρuv)+∂y∂(ρvv)=∂x∂(μ∂x∂v)+∂y∂(μ∂y∂v)−∂y∂p+Sv(1-b)

∂∂x(ρu)+∂∂y(ρv)=0(1-c)\begin{aligned} \frac{\partial}{\partial x}(\rho u) + \frac{\partial}{\partial y}(\rho v)&=0 \qquad \qquad \qquad \qquad\qquad\qquad\qquad\quad \end{aligned} \tag{1-c}∂x∂(ρu)+∂y∂(ρv)=0(1-c)
交错网格的编号如下图：

图1 交错网格

uuu动量方程离散

uuu 动量方程在 uuu 网格的离散方程可以写成如下形式：
ai,Jui,J=Σanbunb−pI,J−pI−1,JδxuΔVu+SˉΔVu(2-a)\begin{aligned} 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{aligned} \tag{2-a}ai,Jui,J=Σanbunb−δxupI,J−pI−1,JΔVu+SˉΔVu(2-a)

或
ai,Jui,J=Σanbunb−(pI,J−pI−1,J)Ai,J+bi,J(2-b)\begin{aligned} 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{aligned} \tag{2-b}ai,Jui,J=Σanbunb−(pI,J−pI−1,J)Ai,J+bi,J(2-b)

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

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

公式中的系数ai,Ja_{i,J}ai,J和anba_{nb}anb的计算对应不同的差分格式有不同的公式，事实上无论选用何种差分格式，方程系数ai,Ja_{i,J}ai,J和anba_{nb}anb都是控制体单元边界处的单位面积对流量F(=ρu)F(=\rho u)F(=ρu)和单位面积扩散量D(=Γ/δx)D(=\Gamma/\delta x)D(=Γ/δx)的组合。在交错网格（二维均匀网格）的编号系统下，uuu网格单元(i,J)(i,J)(i,J)的e、w、ne、w、ne、w、n和sss边界处的FFF值和DDD值的计算公式为
Fw=(ρu)w=Fi,J+Fi−1,J2=12[(ρI,J+ρI−1,J2)ui,J+(ρI−1,J+ρI−2,J2)ui−1,J](3-a)\begin{aligned} 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{aligned} \tag{3-a}Fw=(ρu)w=2Fi,J+Fi−1,J=21[(2ρI,J+ρI−1,J)ui,J+(2ρI−1,J+ρI−2,J)ui−1,J](3-a)

Fe=(ρu)e=Fi+1,J+Fi,J2=12[(ρI+1,J+ρI,J2)ui+1,J+(ρI,J+ρI−1,J2)ui,J](3-b)\begin{aligned} 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{aligned} \tag{3-b}Fe=(ρu)e=2Fi+1,J+Fi,J=21[(2ρI+1,J+ρI,J)ui+1,J+(2ρI,J+ρI−1,J)ui,J](3-b)

Fs=(ρv)s=FI,j+FI−1,j2=12[(ρI,J+ρI,J−12)vI,j+(ρI−1,J+ρI−1,J−12)vI−1,j](3-c)\begin{aligned} 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{aligned} \tag{3-c}Fs=(ρv)s=2FI,j+FI−1,j=21[(2ρI,J+ρI,J−1)vI,j+(2ρI−1,J+ρI−1,J−1)vI−1,j](3-c)

Fn=(ρv)n=FI,j+1+FI−1,j+12=12[(ρI,J+ρI,J+12)vI,j+1+(ρI−1,J+ρI−1,J+12)vI−1,j+1](3-d)\begin{aligned} 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{aligned} \tag{3-d}Fn=(ρv)n=2FI,j+1+FI−1,j+1=21[(2ρI,J+ρI,J+1)vI,j+1+(2ρI−1,J+ρI−1,J+1)vI−1,j+1](3-d)
Dw=ΓI−1,Jxi−xi−1(3-e)\begin{aligned} D_w = \frac{\Gamma_{I-1,J}}{x_i - x_{i-1}} \end{aligned} \tag{3-e}Dw=xi−xi−1ΓI−1,J(3-e)

De=ΓI,Jxi+1−xi(3-f)\begin{aligned} D_e = \frac{\Gamma_{I,J}}{x_{i+1} - x_{i}} \end{aligned} \tag{3-f}De=xi+1−xiΓI,J(3-f)

Ds=Γi,jyJ−yJ−1=ΓI−1,J+ΓI−1,J−1+ΓI,J−1+ΓI,J4(yJ−yJ−1)(3-g)\begin{aligned} 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{aligned} \tag{3-g}Ds=yJ−yJ−1Γi,j=4(yJ−yJ−1)ΓI−1,J+ΓI−1,J−1+ΓI,J−1+ΓI,J(3-g)
Dn=Γi,j+1yJ+1−yJ=ΓI−1,J+1+ΓI−1,J+ΓI,J+ΓI,J+14(yJ+1−yJ)(3-h)\begin{aligned} 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{aligned} \tag{3-h}Dn=yJ+1−yJΓi,j+1=4(yJ+1−yJ)ΓI−1,J+1+ΓI−1,J+ΓI,J+ΓI,J+1(3-h)
在交错网格中，密度ρ\rhoρ和扩散系数Γ\GammaΓ和压力一样都是标量，因此这些标量值都保存在主网格节点（图中的黑色圆点）中，公式(4)−(7)(4) - (7)(4)−(7)中的边界密度值和公式(10)(10)(10)和(11)(11)(11)中的扩散系数都周围节点值的平均值来近似。

vvv动量方程离散

同理，vvv动量方程在vvv网格上的离散方程可以写为如下形式：
aI,jvI,j=Σanbvnb+(pI,J−1−pI,J)AI,j+bI,j(4)\begin{aligned} 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{aligned} \tag{4}aI,jvI,j=Σanbvnb+(pI,J−1−pI,J)AI,j+bI,j(4)
相邻单元的网格编号如下图

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

同样，vvv离散方程中的aI,ja_{I,j}aI,j和anba_{nb}anb也是FFF和DDD的组合，其在vvv网格上的计算公式如下
Fw=(ρu)w=Fi,J+Fi,J−12=12[(ρI,J+ρI−1,J2)ui,J+(ρI−1,J−1+ρI,J−12)ui,J−1](5-a)\begin{aligned} 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{aligned} \tag{5-a}Fw=(ρu)w=2Fi,J+Fi,J−1=21[(2ρI,J+ρI−1,J)ui,J+(2ρI−1,J−1+ρI,J−1)ui,J−1](5-a)

Fe=(ρu)e=Fi+1,J+Fi+1,J−12=12[(ρI,J+ρI+1,J2)ui+1,J+(ρI,J−1+ρI+1,J−12)ui+1,J−1](5-b)\begin{aligned} 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{aligned} \tag{5-b}Fe=(ρu)e=2Fi+1,J+Fi+1,J−1=21[(2ρI,J+ρI+1,J)ui+1,J+(2ρI,J−1+ρI+1,J−1)ui+1,J−1](5-b)

Fs=(ρv)s=FI,j+FI,j−12=12[(ρI,J+ρI,J−12)vI,j+(ρI,J−1+ρI,J−22)vI,j−1](5-c)\begin{aligned} 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{aligned} \tag{5-c}Fs=(ρv)s=2FI,j+FI,j−1=21[(2ρI,J+ρI,J−1)vI,j+(2ρI,J−1+ρI,J−2)vI,j−1](5-c)

Fn=(ρv)n=FI,j+1+FI,j2=12[(ρI,J+ρI,J+12)vI,j+1+(ρI,J−1+ρI,J2)vI,j](5-d)\begin{aligned} 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{aligned} \tag{5-d}Fn=(ρv)n=2FI,j+1+FI,j=21[(2ρI,J+ρI,J+1)vI,j+1+(2ρI,J−1+ρI,J)vI,j](5-d)
Dw=ΓI−1,J+ΓI−1,J−1+ΓI,J−1+ΓI,J4(xI−xI−1)(5-e)\begin{aligned} 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{aligned} \tag{5-e}Dw=4(xI−xI−1)ΓI−1,J+ΓI−1,J−1+ΓI,J−1+ΓI,J(5-e)

De=ΓI+1,J+ΓI+1,J−1+ΓI,J−1+ΓI,J4(xI+1−xI)(5-f)\begin{aligned} 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{aligned} \tag{5-f}De=4(xI+1−xI)ΓI+1,J+ΓI+1,J−1+ΓI,J−1+ΓI,J(5-f)

Ds=ΓI,J−1yj−yj−1(5-g)\begin{aligned} D_s = \frac{\Gamma_{I,J-1}}{y_{j} - y_{j-1}} \end{aligned} \tag{5-g}Ds=yj−yj−1ΓI,J−1(5-g)

Dn=ΓI,Jyj+1−yj(5-h)\begin{aligned} D_n = \frac{\Gamma_{I,J}}{y_{j+1} - y_{j}} \end{aligned} \tag{5-h}Dn=yj+1−yjΓI,J(5-h)

连续性方程离散

连续性方程是在主网格上进行离散的，即在标量p、ρ、Γp、\rho、\Gammap、ρ、Γ的网格上，网格单元及相邻节点如图1的PPP单元所示。连续性方程的离散过程与动量方程类似，其离散后的方程为
[(ρuA)i+1,J−(ρuA)i,J]+[(ρvA)I,j+1−(ρvA)I,j]=0(6-a)\begin{aligned} [(\rho u A)_{i+1,J} - (\rho u A)_{i,J}] + [(\rho v A)_{I,j+1} - (\rho v A)_{I,j}] = 0 \end{aligned} \tag{6-a}[(ρuA)i+1,J−(ρuA)i,J]+[(ρvA)I,j+1−(ρvA)I,j]=0(6-a)

或者写成
[FeAi+1,J−FwAi,J]+[FnAI,j+1−FsAI,j]=0(6-b)\begin{aligned} [F_e A_{i+1,J} - F_w A_{i,J}] + [F_n A_{I,j+1} - F_s A_{I,j}] = 0 \end{aligned} \tag{6-b}[FeAi+1,J−FwAi,J]+[FnAI,j+1−FsAI,j]=0(6-b)
边界处FFF值计算公式如下
Fe=(ρu)e=(ρI,J+ρI+1,J2)ui+1,J(7-a)\begin{aligned} F_e = (\rho u)_e = \left( \frac{\rho_{I,J} + \rho_{I+1,J}}{2} \right )u_{i+1,J} \end{aligned} \tag{7-a}Fe=(ρu)e=(2ρI,J+ρI+1,J)ui+1,J(7-a)

Fw=(ρu)w=(ρI,J+ρI−1,J2)ui,J(7-b)\begin{aligned} F_w = (\rho u)_w = \left( \frac{\rho_{I,J} + \rho_{I-1,J}}{2} \right )u_{i,J} \end{aligned} \tag{7-b}Fw=(ρu)w=(2ρI,J+ρI−1,J)ui,J(7-b)

Fn=(ρv)n=(ρI,J+ρI,J+12)vI,j+1(7-c)\begin{aligned} F_n = (\rho v)_n = \left( \frac{\rho_{I,J} + \rho_{I,J+1}}{2} \right )v_{I,j+1} \end{aligned} \tag{7-c}Fn=(ρv)n=(2ρI,J+ρI,J+1)vI,j+1(7-c)

Fs=(ρv)s=(ρI,J+ρI,J−12)vI,j(7-d)\begin{aligned} F_s = (\rho v)_s = \left( \frac{\rho_{I,J} + \rho_{I,J-1}}{2} \right )v_{I,j} \end{aligned} \tag{7-d}Fs=(ρv)s=(2ρI,J+ρI,J−1)vI,j(7-d)

现在我们得到了交错网格下二维压力速度耦合问题的有限体积法离散方程组，即公式(2)、(4)(2)、(4)(2)、(4)和(6)(6)(6)的组合，有
{ai,Jui,J=Σanbunb−(pI,J−pI−1,J)Ai,J+bi,JaI,jvI,j=Σanbvnb+(pI,J−1−pI,J)AI,j+bI,j[(ρuA)i+1,J−(ρuA)i,J]+[(ρvA)I,j+1−(ρvA)I,j]=0(8)\left \{ \begin{aligned} 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 u A)_{i+1,J} - (\rho u A)_{i,J}] + [(\rho v A)_{I,j+1} - (\rho v A)_{I,j}] = 0 \end{aligned} \right. \tag{8}⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧ai,Jui,J=Σanbunb−(pI,J−pI−1,J)Ai,J+bi,JaI,jvI,j=Σanbvnb+(pI,J−1−pI,J)AI,j+bI,j[(ρuA)i+1,J−(ρuA)i,J]+[(ρvA)I,j+1−(ρvA)I,j]=0(8)

两个关键问题

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

问题1

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

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

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

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

问题2

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

SIMPLE算法

SIMPLE的全称是“Semi-Implicit Method for Pressure-Linked Equations”，即求解压力耦合方程的半隐式方法，“Semi-Implicit”翻译过来就是半隐式，其指代含义后面会提到。SIMPLE算法由Patankar和Spalding（1972）提出来，是在交错网格上计算压力的一个guess-and-correct的过程。下面以笛卡尔坐标系下二维稳态不可压流体的层流流动为例来介绍SIMPLE算法的具体推导过程。
SIMPLE算法首先需要假设一个压力场（p∗p^*p∗）和一个流场（uuu和vvv），然后动量方程就可以在初始压力场和初始流场上计算速度分量uuu和vvv：
ai,Jui,J∗=Σanbunb∗+(pI−1,J∗−pI,J∗)Ai,J+bi,JaI,jvI,j∗=Σanbvnb∗+(pI,J−1∗−pI,J∗)AI,j+bI,j}(9)\left. \begin{aligned} a_{i,J}u^*_{i,J} = \Sigma a_{nb} u^*_{nb} + (p^*_{I-1,J} - p^*_{I,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} \end{aligned} \right \} \tag{9}ai,Jui,J∗=Σanbunb∗+(pI−1,J∗−pI,J∗)Ai,J+bi,JaI,jvI,j∗=Σanbvnb∗+(pI,J−1∗−pI,J∗)AI,j+bI,j⎭⎪⎬⎪⎫(9)
上式中，等号右边的u∗、v∗、p∗u^*、v^*、p^*u∗、v∗、p∗都是假设的初始值，等号左边的u∗、v∗u^*、v^*u∗、v∗是计算的初始速度分量，需要指出的是，系数ai,J、aI,j、anba_{i,J}、a_{I,j}、a_{nb}ai,J、aI,j、anb的计算公式中也有u∗u^*u∗和v∗v^*v∗，也是使用假设值计算。然而，一般情况下，这样计算得到的速度分量uuu和vvv是不满足连续性方程的，因为假设的压力和速度分布一般是与真实情况是不符的。因此，为了得到合理的计算结果，就需要对压力场p∗p^*p∗和速度场u∗、v∗u^*、v^*u∗、v∗进行修正。定义压力的修正量和修正后的正确性分别为p′p^\primep′和ppp，速度的修正量和修正后的正确值分别为u′、v′u^\prime 、v^\primeu′、v′和u、vu、vu、v，则其之间的关系为
p=p∗+p′(10-a)\begin{aligned} p=p^* + p^\prime \end{aligned} \tag{10-a}p=p∗+p′(10-a)

u=u∗+u′(10-b)\begin{aligned} u=u^* + u^\prime \end{aligned} \tag{10-b}u=u∗+u′(10-b)

v=v∗+v′(10-c)\begin{aligned} v=v^* + v^\prime \end{aligned} \tag{10-c}v=v∗+v′(10-c)
将压力场的正确值ppp带入到动量方程(2−b)(2-b)(2−b)和(4)(4)(4)就能计算出正确的速度场(u,v)(u,v)(u,v)，如果用方程(2−b)(2-b)(2−b)和(4)(4)(4)减去公式(9)(9)(9)的方程，则得到
ai,J(ui,J−ui,J∗)=Σanb(unb−unb∗)+[(pI−1,J−pI−1,J∗)−(pI,J−pI,J∗)]Ai,J(11-a)\begin{aligned} a_{i,J} (u_{i,J} - u^*_{i,J}) = \Sigma a_{nb}(u_{nb} - u^*_{nb}) + [(p_{I-1,J}-p^*_{I-1,J}) - (p_{I,J}-p^*_{I,J})]A_{i,J} \end{aligned} \tag{11-a}ai,J(ui,J−ui,J∗)=Σanb(unb−unb∗)+[(pI−1,J−pI−1,J∗)−(pI,J−pI,J∗)]Ai,J(11-a)

aI,j(vI,j−vI,j∗)=Σanb(vnb−vnb∗)+[(pI,J−1−pI,J−1∗)−(pI,J−pI,J∗)]AI,j(11-a)\begin{aligned} a_{I,j} (v_{I,j} - v^*_{I,j}) = \Sigma a_{nb}(v_{nb} - v^*_{nb}) + [(p_{I,J-1}-p^*_{I,J-1}) - (p_{I,J}-p^*_{I,J})]A_{I,j} \end{aligned} \tag{11-a}aI,j(vI,j−vI,j∗)=Σanb(vnb−vnb∗)+[(pI,J−1−pI,J−1∗)−(pI,J−pI,J∗)]AI,j(11-a)

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

把公式(10)(10)(10)的关系带入到公式(11)(11)(11)中，将得到修正值的方程：
ai,Jui,J′=Σanbunb′+(pI−1,J′−pI,J′)Ai,J(12-a)\begin{aligned} a_{i,J} u^\prime_{i,J} = \Sigma a_{nb} u^\prime_{nb} + (p^\prime_{I-1,J}-p^\prime_{I,J})A_{i,J} \end{aligned} \tag{12-a}ai,Jui,J′=Σanbunb′+(pI−1,J′−pI,J′)Ai,J(12-a)

aI,jvI,j′=Σanbvnb′+(pI,J−1′−pI,J′)AI,j(12-b)\begin{aligned} a_{I,j} v^\prime_{I,j} = \Sigma a_{nb} v^\prime_{nb} + (p^\prime_{I,J-1}-p^\prime_{I,J})A_{I,j} \end{aligned} \tag{12-b}aI,jvI,j′=Σanbvnb′+(pI,J−1′−pI,J′)AI,j(12-b)
这里引入一个近似：将Σanbunb′\Sigma a_{nb}u^\prime_{nb}Σanbunb′和Σanbvnb′\Sigma a_{nb}v^\prime_{nb}Σanbvnb′两项省略掉。这两项的省略是SIMPLE算法中主要的近似操作。然后，我们就得到
ui,J′=di,J(pI−1,J′−pI,J′)(13-a)\begin{aligned} u^\prime_{i,J} = d_{i,J} (p^\prime_{I-1,J} - p^\prime_{I,J}) \end{aligned} \tag{13-a}ui,J′=di,J(pI−1,J′−pI,J′)(13-a)

vI,j′=dI,j(pI,J−1′−pI,J′)(13-b)\begin{aligned} v^\prime_{I,j} = d_{I,j} (p^\prime_{I,J-1} - p^\prime_{I,J}) \end{aligned} \tag{13-b}vI,j′=dI,j(pI,J−1′−pI,J′)(13-b)
其中，di,J=Ai,Jai,Jd_{i,J}=\frac{A_{i,J}}{a_{i,J}}di,J=ai,JAi,J，dI,j=AI,jaI,jd_{I,j}=\frac{A_{I,j}}{a_{I,j}}dI,j=aI,jAI,j.

把公式(13−a)(13-a)(13−a)和(13−b)(13-b)(13−b)代入到公式(10−b)(10-b)(10−b)和(10−c)(10-c)(10−c)中，有
ui,J=ui,J∗+di,J(pI−1,J′−pI,J′)(14-a)\begin{aligned} u_{i,J}=u^*_{i,J} + d_{i,J} (p^\prime_{I-1,J} - p^\prime_{I,J}) \end{aligned} \tag{14-a}ui,J=ui,J∗+di,J(pI−1,J′−pI,J′)(14-a)

vI,j=vI,j∗+dI,j(pI,J−1′−pI,J′)(14-b)\begin{aligned} v_{I,j} = v^*_{I,j} + d_{I,j} (p^\prime_{I,J-1} - p^\prime_{I,J}) \end{aligned} \tag{14-b}vI,j=vI,j∗+dI,j(pI,J−1′−pI,J′)(14-b)
同样对于ui+1,Ju_{i+1,J}ui+1,J和vI,j+1v_{I,j+1}vI,j+1有
ui+1,J=ui+1,J∗+di+1,J(pI,J′−pI+1,J′)(15-a)\begin{aligned} u_{i+1,J}=u^*_{i+1,J} + d_{i+1,J} (p^\prime_{I,J} - p^\prime_{I+1,J}) \end{aligned} \tag{15-a}ui+1,J=ui+1,J∗+di+1,J(pI,J′−pI+1,J′)(15-a)

vI,j+1=vI,j+1∗+dI,j+1(pI,J′−pI,J+1′)(15-b)\begin{aligned} v_{I,j+1} = v^*_{I,j+1} + d_{I,j+1} (p^\prime_{I,J} - p^\prime_{I,J+1}) \end{aligned} \tag{15-b}vI,j+1=vI,j+1∗+dI,j+1(pI,J′−pI,J+1′)(15-b)
其中，di+1,J=Ai+1,Jai+1,Jd_{i+1,J}=\frac{A_{i+1,J}}{a_{i+1,J}}di+1,J=ai+1,JAi+1,J，dI,j+1=AI,j+1aI,j+1d_{I,j+1}=\frac{A_{I,j+1}}{a_{I,j+1}}dI,j+1=aI,j+1AI,j+1.

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

按上图的网格标号，离散的连续性方程则成为
[ρi+1,JAi+1,J(ui+1,J∗+di+1,J(pI,J′−pI+1,J′))−ρi,JAi,J(ui,J∗+di,J(pI−1,J′−pI,J′))]+[ρI,j+1AI,j+1(vI,j+1∗+dI,j+1(pI,J′−pI,J+1′))−ρI,jAI,j(vI,j∗+dI,j(pI,J−1′−pI,J′))]=0(17)\begin{aligned} [\rho_{i+1,J} & A_{i+1,J}(u^*_{i+1,J}+d_{i+1,J}(p^\prime_{I,J}-p^\prime_{I+1,J})) - \rho_{i,J}A_{i,J}(u^*_{i,J} + \\\\ &d_{i,J}(p^\prime_{I-1,J}-p^\prime_{I,J})) ] + [\rho_{I,j+1}A_{I,j+1}(v^*_{I,j+1} + d_{I,j+1}(p^\prime_{I,J} - p^\prime_{I,J+1})) \\\\ &-\rho_{I,j}A_{I,j}(v^*_{I,j} + d_{I,j}(p^\prime_{I,J-1} - p^\prime_{I,J}))] =0 \end{aligned} \tag{17}[ρi+1,JAi+1,J(ui+1,J∗+di+1,J(pI,J′−pI+1,J′))−ρi,JAi,J(ui,J∗+di,J(pI−1,J′−pI,J′))]+[ρI,j+1AI,j+1(vI,j+1∗+dI,j+1(pI,J′−pI,J+1′))−ρI,jAI,j(vI,j∗+dI,j(pI,J−1′−pI,J′))]=0(17)
把p′p^\primep′看作未知变量，重新整理一下，就变成
[(ρdA)i+1,J+(ρdA)i,J+(ρdA)I,j+1+(ρdA)I,J]pI,J′=(ρdA)i+1,JpI+1,J′+(ρdA)i,JpI−1,J′+(ρdA)I,J+1pI,J+1′+(ρdA)I,jpI,J−1′+[(ρu∗A)i,J−(ρu∗A)i+1,J+(ρv∗A)I,j−(ρv∗A)I,j+1](18)\begin{aligned} [(\rho d A)_{i+1,J} &+ (\rho d A)_{i,J} + (\rho d A)_{I,j+1} + (\rho d A)_{I,J} ] p^\prime_{I,J} = \\\\&(\rho d A)_{i+1,J} p^\prime_{I+1,J} + (\rho d A)_{i,J} p^\prime_{I-1,J} + (\rho d A)_{I,J+1} p^\prime_{I,J+1} + (\rho d A)_{I,j} p^\prime_{I,J-1} \\\\ &+[(\rho u^*A)_{i,J} - (\rho u^*A)_{i+1,J} + (\rho v^*A)_{I,j} - (\rho v^*A)_{I,j+1}] \end{aligned} \tag{18}[(ρdA)i+1,J+(ρdA)i,J+(ρdA)I,j+1+(ρdA)I,J]pI,J′=(ρdA)i+1,JpI+1,J′+(ρdA)i,JpI−1,J′+(ρdA)I,J+1pI,J+1′+(ρdA)I,jpI,J−1′+[(ρu∗A)i,J−(ρu∗A)i+1,J+(ρv∗A)I,j−(ρv∗A)I,j+1](18)
写出熟悉的样子就是
aI,JpI,J′=aI+1,JpI+1,J′+aI−1,JpI−1,J′+aI,J+1pI,J+1′+aI,J−1pI,J−1′+bI,J′(19)\begin{aligned} a_{I,J} p^\prime_{I,J} = a_{I+1,J} p^\prime_{I+1,J} + a_{I-1,J} p^\prime_{I-1,J} +a_{I,J+1} p^\prime_{I,J+1} + a_{I,J-1} p^\prime_{I,J-1} + b^\prime_{I,J} \end{aligned} \tag{19}aI,JpI,J′=aI+1,JpI+1,J′+aI−1,JpI−1,J′+aI,J+1pI,J+1′+aI,J−1pI,J−1′+bI,J′(19)
其中系数为

aI+1,Ja_{I+1,J}aI+1,J	(ρdA)i+1,J(\rho d A)_{i+1,J}(ρdA)i+1,J
aI−1,Ja_{I-1,J}aI−1,J	(ρdA)i,J(\rho d A)_{i,J}(ρdA)i,J
aI,J+1a_{I,J+1}aI,J+1	(ρdA)I,j+1(\rho d A)_{I,j+1}(ρdA)I,j+1
aI,J−1a_{I,J-1}aI,J−1	(ρdA)I,j(\rho d A)_{I,j}(ρdA)I,j
bI,Jb_{I,J}bI,J	(ρu∗A)i,J−(ρu∗A)i+1,J+(ρv∗A)I,j−(ρv∗A)I,j+1(\rho u^* A)_{i,J}-(\rho u^* A)_{i+1,J} +(\rho v^* A)_{I,j} - (\rho v^* A)_{I,j+1}(ρu∗A)i,J−(ρu∗A)i+1,J+(ρv∗A)I,j−(ρv∗A)I,j+1
aI,Ja_{I,J}aI,J	aI+1,J+aI−1,J+aI,J+1+aI,J−1a_{I+1,J} + a_{I-1,J} + a_{I,J+1} + a_{I,J-1}aI+1,J+aI−1,J+aI,J+1+aI,J−1

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

关于SIMPLE算法的几点说明

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

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

【3】 一般情况下，如果按公式(10)(10)(10)直接用压力修正值来更新压力场往往会导致迭代发散。使用亚松弛迭代方法可以解决这一问题，压力修正方程就成为
pnew=p∗+αpp′(20)\begin{aligned} p^{new} = p^* + \alpha_p p^\prime \end{aligned} \tag{20}pnew=p∗+αpp′(20)

式中，αp\alpha_pαp是亚松弛因子。因为修正量随着迭代计算会区域零，所以亚松弛因子的大小并不会影响最终计算的结果。
当αp=1\alpha_p=1αp=1时，实际上就是直接用压力修正值p′p^\primep′来修正压力初始值p∗p^*p∗。然而，当压力初始值与真实值相差较大时，这样计算出来的压力修正值p′p^\primep′会是很大的数，如果直接用p′p^\primep′来修正p∗p^*p∗，即使p′p^\primep′的修正方向是正确的，那也会导致矫枉过正，造成计算的不稳定。
当αp\alpha_pαp取0~1之间的值时，即人为控制压力修正的幅度。每次迭代压缩了修正的幅度，控制了压力变动的幅度，从而保证计算的稳定性。然而，αp\alpha_pαp的值也并不是越小越好，过小的松弛因子虽然可以保证收敛稳定，但要迭代到最后的收敛解就需要更多次的迭代计算，这样会导致收敛速度很慢。计算中所希望的αp\alpha_pαp值是尽可能大（以加快收敛）但又不引起计算发散。
速度的修正也需要采用亚松弛迭代，
unew=u(n−1)+αuu′=u(n−1)+αu(u(n)−u(n−1))=αuu(n)+(1−αu)u(n−1)(21-a)\begin{aligned} 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{aligned} \tag{21-a}unew=u(n−1)+αuu′=u(n−1)+αu(u(n)−u(n−1))=αuu(n)+(1−αu)u(n−1)(21-a)
vnew=v(n−1)+αvv′=v(n−1)+αv(v(n)−v(n−1))=αvv(n)+(1−αv)v(n−1)(21-b)\begin{aligned} 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{aligned} \tag{21-b}vnew=v(n−1)+αvv′=v(n−1)+αv(v(n)−v(n−1))=αvv(n)+(1−αv)v(n−1)(21-b)
其中，αu\alpha_uαu和αv\alpha_vαv是速度uuu和vvv修正的亚松弛因子，取值也是在000~111之间，u(n−1)u^{(n-1)}u(n−1)和v(n−1)v^{(n-1)}v(n−1)是上一次迭代计算的结果，u(n)u^{(n)}u(n)和v(n)v^{(n)}v(n)是本次迭代以p′p^\primep′直接修正后的值。

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

参考资料

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