学习自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

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

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

这里是第一部分，主要讲解不可压缩流动问题的求解在常规网格上所碰到的问题，以及，交错网格上的SIMPLE算法是如何解决该问题的。

1 主要困难

前面几章所处理的通用守恒方程，可以改写成与连续方程或动量方程相似的形式。然而，目前为止所呈现的数值技术，并不足以求解Navier-Stokes方程，因为求解一般流动的算法要能够处理压力速度的耦合问题。为理解该问题，先把连续方程和动量方程写出来
∂ρ∂t+∇⋅(ρv)=0∂∂t(ρv)+∇⋅(ρvv)=−∇p+∇⋅{μ[∇v+(∇v)T]}+fb\begin{aligned} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bold v)&=0 \\ \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 \end{aligned} ∂t∂ρ+∇⋅(ρv)∂t∂(ρv)+∇⋅(ρvv)=0=−∇p+∇⋅{μ[∇v+(∇v)T]}+fb
这俩方程是非线性的，但也并不是不可克服的困难，因为这种问题通常可采用迭代方法来应对。此外，第二式是矢量方程，若把其展开成分量形式的话，则是个标量方程系统，需要顺序求解。再者，应力张量可改写成类似扩散项的形式，并可隐式处理，其第2项（即，速度梯度的转置(∇v)T(\nabla\bold v)^\text{T}(∇v)T）则基于前次迭代值做显式估算并添加到源项中去。针对上述连续方程和动量方程，从通有标量方程的数值方法（前面几章讲解的方法）中无法汲取灵感来解决的主要难题在于，出现在动量方程中的压力场无法获得其显式计算方程。

从这俩方程中可以看出，纵然速度场可以直接使用动量方程来计算，然而出现在动量方程中的压力场无法直接从连续方程中计算。该耦合关系（压力和速度的耦合）是如此强烈却又如此隐含，为便于理解，可把方程组写成如下矩阵形式
Au=[FBTB0][vp]=[fb0]\bold A \bold u=\begin{bmatrix} \bold F & \bold B^\text T \\ \bold B & \bold 0 \end{bmatrix} \begin{bmatrix} \bold v \\ \bold p \end{bmatrix}= \begin{bmatrix} \bold f_b \\ \bold 0 \end{bmatrix} Au=[FBBT0][vp]=[fb0]
在该形式中，上式中出现了一个零对角块，这是鞍点（saddle point）问题的特性，表明其无法通过任何迭代手法来求解压力和速度场。因此，需要推导出关于压力的方程才好。

一种处理方法是，通过把矩阵A\bold AA分解成下三角L\bold LL和上三角U\bold UU矩阵，将动量和连续方程系统重新改写成下述形式
A=[FBTB0]=[F0B−BF−1BT][IF−1BT0I]=LU\bold A = \begin{bmatrix} \bold F & \bold B^\text T \\ \bold B & \bold 0 \end{bmatrix}= \begin{bmatrix} \bold F & \bold 0 \\ \bold B & -\bold B\bold F^{-1}\bold B^\text T \end{bmatrix} \begin{bmatrix} \bold I & \bold F^{-1}\bold B^\text T \\ \bold 0 & \bold I \end{bmatrix}= \bold L \bold U A=[FBBT0]=[FB0−BF−1BT][I0F−1BTI]=LU
其中−BF−1BT-\bold B\bold F^{-1}\bold B^\text T−BF−1BT为Schur补矩阵（Schur complement matrix）。

这正是下面要讲的迭代求解Navier-Stokes方程的方法的本质，Patankar和Spalding将该技术嵌入在经典的分离式SIMPLE（Semi Implicit Method for Pressure Linked Equations）算法中。

求解流程是，将Navier-Stokes方程重新整理成动量方程和压力方程的形式，然后做离散和顺序求解。压力方程是通过联合半离散动量方程和连续方程（Schur补矩阵）来创建的。

算法是由Picard类型的迭代流程驱动的，其率先求解动量方程，使用的是前次迭代的压力场。求出的速度场是满足动量方程的但是未必满足质量方程。接下来用该速度场来构建压力方程，并用所得解来修正压力和速度场以便满足质量守恒。然后开始新的迭代过程，并重复上述流程，直到速度和压力场既满足质量守恒又满足动量守恒为止。

该算法可用如下的矩阵形式来描述
[ID−1BT0I][vp]=[v∗p∗]\begin{bmatrix} \bold I & \bold D^{-1}\bold B^\text T \\ \bold 0 & \bold I \end{bmatrix}\begin{bmatrix} \bold v \\ \bold p \end{bmatrix}= \begin{bmatrix} \bold v^* \\ \bold p^* \end{bmatrix} [I0D−1BTI][vp]=[v∗p∗]
接下来使用下式来更新速度场
[F0B−BD−1BT][v∗p∗]=[fb0]\begin{bmatrix} \bold F & \bold 0 \\ \bold B & -\bold B\bold D^{-1}\bold B^\text T \end{bmatrix} \begin{bmatrix} \bold v^* \\ \bold p^* \end{bmatrix}= \begin{bmatrix} \bold f_b \\ \bold 0 \end{bmatrix} [FB0−BD−1BT][v∗p∗]=[fb0]
上式和上上式中的F−1\bold F^{-1}F−1用其逆对角阵D−1\bold D^{-1}D−1来近似，上标∗^*∗表示当前迭代的中间值。流程汇总如下：

求解：Fv∗=fb\bold F \bold v^*=\bold f_bFv∗=fb
求解：−BD−1BTp∗=−Bv∗-\bold B \bold D^{-1} \bold B^\text T \bold p^*=-\bold B \bold v^*−BD−1BTp∗=−Bv∗
更新：v=v∗−D−1BTp∗\bold v=\bold v^*-\bold D^{-1}\bold B^\text T\bold p^*v=v∗−D−1BTp∗
更新：p=p∗\bold p=\bold p^*p=p∗

这种类型的分裂与SIMPLE系列算法中所用的分裂是相同的，这正是本章的主题。

2 初步推导

针对不可压缩流动问题发展求解算法的困难之处，可通过对上图所示的一维空间均匀网格问题的离散处理来加深理解。为简便起见，假设流动是定常的。简化的连续和动量方程（写成守恒形式）为
∂(ρu)∂x=0∂(ρuu)∂x=∂∂x(μ∂u∂x)−∂p∂x\begin{aligned} \frac{\partial (\rho u)}{\partial x}&=0 \\ ~ \\ \frac{\partial(\rho u u)}{\partial x}&=\frac{\partial}{\partial x} \left(\mu \frac{\partial u}{\partial x} \right) -\frac{\partial p}{\partial x} \end{aligned} ∂x∂(ρu) ∂x∂(ρuu)=0=∂x∂(μ∂x∂u)−∂x∂p

2.1 动量方程离散

动量方程的离散，首先要对其在单元CCC上做积分，得
∫VC∂(ρuu)∂xdV=∫VC∂∂x(μ∂u∂x)dV−∫VC∂p∂xdV\int_{V_C}\frac{\partial(\rho u u)}{\partial x}dV= \int_{V_C}\frac{\partial}{\partial x} \left(\mu \frac{\partial u}{\partial x} \right)dV- \int_{V_C}\frac{\partial p}{\partial x}dV ∫VC∂x∂(ρuu)dV=∫VC∂x∂(μ∂x∂u)dV−∫VC∂x∂pdV
使用散度定理，将对流项和扩散项的体积分转化成面积分
∫∂VC(ρuu)dy=∫∂VCμ∂u∂xdy−∫VC∂p∂xdV\int_{\partial V_C}(\rho u u)dy=\int_{\partial V_C}\mu \frac{\partial u}{\partial x}dy - \int_{V_C}\frac{\partial p}{\partial x}dV ∫∂VC(ρuu)dy=∫∂VCμ∂x∂udy−∫VC∂x∂pdV
把面积分用流过单元面的通量加和来替换，并使用单Gauss点来做面积分，得到半离散形式为
(ρuΔy)e⏟m˙eue+−(ρuΔy)w⏟m˙wuw=(μ∂u∂xΔy)e−(μ∂u∂xΔy)w−∫VC∂p∂xdV\underbrace{(\rho u \Delta y)_e}_{\dot m_e} u_e+ \underbrace{-(\rho u \Delta y)_w}_{\dot m_w} u_w= \left( \mu \frac{\partial u}{\partial x}\Delta y \right)_e- \left( \mu \frac{\partial u}{\partial x}\Delta y \right)_w- \int_{V_C}\frac{\partial p}{\partial x}dV m˙e(ρuΔy)eue+m˙w−(ρuΔy)wuw=(μ∂x∂uΔy)e−(μ∂x∂uΔy)w−∫VC∂x∂pdV
其可重写为
m˙eue+m˙wuw⏟Convection−[(μ∂u∂xΔy)e−(μ∂u∂xΔy)w]⏟Diffusion=−∫VC∂p∂xdV\underbrace{\dot m_e u_e+\dot m_w u_w}_{Convection}- \underbrace{\left[\left( \mu \frac{\partial u}{\partial x}\Delta y \right)_e- \left( \mu \frac{\partial u}{\partial x}\Delta y \right)_w\right]}_{Diffusion}=- \int_{V_C}\frac{\partial p}{\partial x}dV Convectionm˙eue+m˙wuw−Diffusion[(μ∂x∂uΔy)e−(μ∂x∂uΔy)w]=−∫VC∂x∂pdV
对流和扩散项可用之前章节中所讲的方法来离散化处理，得到如下形式的代数方程
aCuuC+∑F∼NB(C)(aFuuF)=bCu−∫VC∂p∂xdVa_C^u u_C+\sum_{F\sim NB(C)}(a_F^u u_F)=b_C^u-\int_{V_C}\frac{\partial p}{\partial x}dV aCuuC+F∼NB(C)∑(aFuuF)=bCu−∫VC∂x∂pdV
压力项的离散待咱们讲完连续方程的离散后再说。

2.2 连续方程的离散

连续方程的离散，也是要对其在单元CCC上做积分，得
∫VC∂(ρu)∂xdV=0\int_{V_C}\frac{\partial(\rho u)}{\partial x}dV=0 ∫VC∂x∂(ρu)dV=0
再次使用散度定理将体积分转化成面积分，然后转化成单元面通量的加和形式，连续方程的离散形式为
∑f∼nb(C)(ρuΔy)f=(ρuΔy)e−(ρuΔy)w=0\sum_{f\sim nb(C)}(\rho u \Delta y)_f=(\rho u \Delta y)_e - (\rho u \Delta y)_w=0 f∼nb(C)∑(ρuΔy)f=(ρuΔy)e−(ρuΔy)w=0
或
∑f∼nb(C)m˙f=m˙e+m˙w=0\sum_{f\sim nb(C)}\dot m_f=\dot m_e + \dot m_w=0 f∼nb(C)∑m˙f=m˙e+m˙w=0

2.3 棋盘问题（Checkerboard Problem）

压力项离散可由下面两种方法中的任意一个来实现。第一种方法，通过单点Gauss积分来计算体积分
∫VC∂p∂xdV=(∂p∂x)CVC\int_{V_C}\frac{\partial p}{\partial x}dV=\left(\frac{\partial p}{\partial x}\right)_CV_C ∫VC∂x∂pdV=(∂x∂p)CVC
使用中心差分格式，上式的离散形式为
∫VC∂p∂xdV=pE−pW2ΔxVC\int_{V_C}\frac{\partial p}{\partial x}dV=\frac{p_E-p_W}{2\Delta x} V_C ∫VC∂x∂pdV=2ΔxpE−pWVC
第二种方法，把压力梯度项的体积分转化成面积分
∫VC∂p∂xdV=∫∂VCpdy\int_{V_C}\frac{\partial p}{\partial x}dV=\int_{\partial V_C}p dy ∫VC∂x∂pdV=∫∂VCpdy
把面积分写成单元面通量加和的形式，有
∫VC∂p∂xdV=∫∂VCpdy=peΔye−pwΔyw=(pe−pw)Δy=(pe−pw)VCΔx\int_{V_C}\frac{\partial p}{\partial x}dV=\int_{\partial V_C}p dy= p_e\Delta y_e - p_w \Delta y_w=(p_e-p_w)\Delta y= (p_e-p_w)\frac{V_C}{\Delta x} ∫VC∂x∂pdV=∫∂VCpdy=peΔye−pwΔyw=(pe−pw)Δy=(pe−pw)ΔxVC
压力的变化选用线性插值廓线，把压力梯度项重新写成主节点上压力值的函数形式
∫VC∂p∂xdV=[12(pE+pC)−12(pC+pW)]VCΔx=pE−pW2ΔxVC\int_{V_C}\frac{\partial p}{\partial x}dV= \left[\frac12(p_E+p_C)-\frac12(p_C+p_W)\right]\frac{V_C}{\Delta x}= \frac{p_E-p_W}{2\Delta x} V_C ∫VC∂x∂pdV=[21(pE+pC)−21(pC+pW)]ΔxVC=2ΔxpE−pWVC
两个方法导出了同样的表达式，即，包含交替节点EEE和WWW之间压力差值的表达式（这里的交替指的是，EEE和WWW节点之间间隔了一个节点CCC）。

采用同样的方法，使用线性插值廓线，并注意到密度是常数，且有(Δy)e=(Δy)w=(Δy)C(\Delta y)_e=(\Delta y)_w=(\Delta y)_C(Δy)e=(Δy)w=(Δy)C，则连续方程为
uE−uW=0u_E-u_W=0 uE−uW=0
用的也是两个交替网格节点的速度值。

在上上式中，单元CCC处的压力梯度项是由两个交替的（而非连续的）横跨单元两侧的网格节点上的压力值所决定的。连续方程亦是如此，这样一来就只是对交替的速度单元施加了守恒特性。这意味着不符合物理意义的之字形（或棋盘型，国际象棋中的黑白格子间隔排列）的压力场或速度场（如下图所示的那样），在该数值方法中将会被感知成均匀场。换句话说，像这种…-A-B-A-B-…排列的场，原本并非是均匀场，但是由于数值算法中的压力梯度和连续方程均考虑的是间隔单元间的特性，所以反倒是把其错误滴认定为是均匀场了。

针对上图所示的压力和速度值，在点WWW、CCC、EEE处的压力梯度为
∫VW∂p∂xdV=(pC−pWW)VW2ΔxW=(10−10)VW2ΔxW=0∫VC∂p∂xdV=(pE−pW)VC2ΔxC=(−100+100)VC2ΔxC=0∫VE∂p∂xdV=(pEE−pC)VE2ΔxE=(10−10)VE2ΔxE=0\begin{aligned} &\int_{V_W}\frac{\partial p}{\partial x}dV=(p_C-p_{WW})\frac{V_W}{2\Delta x_W}=(10-10)\frac{V_W}{2\Delta x_W}=0 \\ &\int_{V_C}\frac{\partial p}{\partial x}dV=(p_E-p_{W})\frac{V_C}{2\Delta x_C}=(-100+100)\frac{V_C}{2\Delta x_C}=0 \\ &\int_{V_E}\frac{\partial p}{\partial x}dV=(p_{EE}-p_{C})\frac{V_E}{2\Delta x_E}=(10-10)\frac{V_E}{2\Delta x_E}=0 \\ \end{aligned} ∫VW∂x∂pdV=(pC−pWW)2ΔxWVW=(10−10)2ΔxWVW=0∫VC∂x∂pdV=(pE−pW)2ΔxCVC=(−100+100)2ΔxCVC=0∫VE∂x∂pdV=(pEE−pC)2ΔxEVE=(10−10)2ΔxEVE=0
而且连续方程看起来在每个单元上都是满足的，因为
∫VW∂u∂xdV=(uC−uWW)VW2ΔxW=(1−1)VW2ΔxW=0∫VC∂u∂xdV=(uE−uW)VC2ΔxC=(10−10)VC2ΔxC=0∫VE∂u∂xdV=(uEE−uC)VE2ΔxE=(1−1)VE2ΔxE=0\begin{aligned} & \int_{V_W}\frac{\partial u}{\partial x}dV=(u_C-u_{WW})\frac{V_W}{2\Delta x_W}=(1-1)\frac{V_W}{2\Delta x_W}=0 \\ & \int_{V_C}\frac{\partial u}{\partial x}dV=(u_E-u_{W})\frac{V_C}{2\Delta x_C}=(10-10)\frac{V_C}{2\Delta x_C}=0 \\ & \int_{V_E}\frac{\partial u}{\partial x}dV=(u_{EE}-u_{C})\frac{V_E}{2\Delta x_E}=(1-1)\frac{V_E}{2\Delta x_E}=0 \\ \end{aligned} ∫VW∂x∂udV=(uC−uWW)2ΔxWVW=(1−1)2ΔxWVW=0∫VC∂x∂udV=(uE−uW)2ΔxCVC=(10−10)2ΔxCVC=0∫VE∂x∂udV=(uEE−uC)2ΔxEVE=(1−1)2ΔxEVE=0
在多维情况下，也会产生相似的非物理特性，即使很难具象化展示（其实如果是二维结构网格的话，这种情况跟国际象棋的黑白格子就很像了，即，黑格子是一个值，而白格子是另一个值，但是按照上面这种数值方法处理起来，整场各个单元反倒是满足了连续方程，而且压力梯度也都是零了）。这为下面展示一种处理该问题的方法做好了铺垫。

2.4 交错网格

前面公式中所呈现的问题是压力和速度场的解耦问题。可以通过把不同的变量存储在交错位置上来加强它们的耦合特性，因为这样一来，就不需要用插值手段来计算动量方程中的压力梯度以及连续方程中的速度场了。这样的交错网格如下图所示，在交错网格中速度场是存储在单元面上的（图a），而压力场则是存储在单元形心上的（图b）。

单元CCC的连续方程离散为
∑f∼nb(C)m˙f=m˙e+m˙w=0\sum_{f\sim nb(C)}\dot m_f=\dot m_e+\dot m_w=0 f∼nb(C)∑m˙f=m˙e+m˙w=0
或
ue−uw=0u_e-u_w=0 ue−uw=0
无需再做插值了，因为速度在界面eee和www处是可以获得的。此外，动量方程的积分是在单元eee上进行的，其产生如下的离散动量方程形式
aeuue+∑f∼NB(e)afuuf=beu−Ve(∇p)e=beu−VepE−pCδxea_e^u u_e+\sum_{f\sim NB(e)}a_f^u u_f=b_e^u-V_e(\nabla p)_e= b_e^u-V_e\frac{p_E-p_C}{\delta x_e} aeuue+f∼NB(e)∑afuuf=beu−Ve(∇p)e=beu−VeδxepE−pC
压力梯度的计算用的是横跨单元面两侧的连续单元形心节点值，故不需要做插值。因此棋盘压力场和速度场的情况将不会出现，因为这些情况将很容易地被数值方法探测并衰减掉。

2.5 压力修正方程

接下来要呈现的推导是基于Patankar和Spalding的工作，他们发展了SIMPLE（Semi Implicit Method for Pressure Linked Equations）算法的最初实用形式。

由上图所示网格上的离散连续方程和动量方程出发
∑f∼nb(C)m˙f=0aeuue+∑f∼NB(e)afuuf=beu−Ve(∂p∂x)e\begin{aligned} \sum_{f\sim nb(C)}\dot m_f &= 0 \\ a_e^u u_e+\sum_{f\sim NB(e)}a_f^u u_f&= b_e^u-V_e\left(\frac{\partial p}{\partial x}\right)_e \end{aligned} f∼nb(C)∑m˙faeuue+f∼NB(e)∑afuuf=0=beu−Ve(∂x∂p)e
求解过程首先需要提供初始的猜测速度和压力场。把初始猜测的或是迭代开始时的解用上标(n)^{(n)}(n)表示，那么这些个暂定的（初始猜测的或是未达到收敛的不确定的）速度和压力场用u(n)u^{(n)}u(n)和p(n)p^{(n)}p(n)表示。在每步迭代中，先求解动量方程来获取速度场，所得解用上标∗^*∗表示，因为其并非当前迭代步的最终解（注意，由于并未考虑连续方程，所以这时候的速度场并不满足连续方程）。这样，动量方程满足
aeuue∗+∑f∼NB(e)afuuf∗=beu−Ve(∂p(n)∂x)ea_e^u u_e^*+\sum_{f\sim NB(e)}a_f^u u_f^*= b_e^u-V_e\left(\frac{\partial p^{(n)}}{\partial x}\right)_e aeuue∗+f∼NB(e)∑afuuf∗=beu−Ve(∂x∂p(n))e
其中压力场仍旧采用的是前次迭代值，那么，算出来的速度场u∗u^*u∗满足动量方程，但是未必满足连续方程，因为压力场并不是精确值。因此，需要寻找修正措施来确保速度（或质量流量）和压力满足连续方程。

将修正场用上标撇来表示，即，（u′u'u′、p′p'p′），那么速度和压力修正为
u=u∗+u′p=p∗+p′u=u^*+u' \\ p=p^*+p' u=u∗+u′p=p∗+p′
注意，单元面处的质量流量也修正为
m˙f=m˙f∗+ρu′Sfx=m˙f∗+mf′\dot m_f=\dot m_f^*+\rho u' S_f^x=\dot m_f^*+m_f' m˙f=m˙f∗+ρu′Sfx=m˙f∗+mf′
这样，精确质量流量应满足连续方程，即
m˙e+m˙w=m˙e∗+m˙e′+m˙w∗+m˙w′=0\dot m_e+\dot m_w=\dot m_e^*+\dot m_e'+\dot m_w^*+\dot m_w'=0 m˙e+m˙w=m˙e∗+m˙e′+m˙w∗+m˙w′=0
可以写为
m˙e′+m˙w′=−m˙e∗−m˙w∗\dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^* m˙e′+m˙w′=−m˙e∗−m˙w∗
这个形式的连续方程非常有趣，其表明，一旦算得的质量流量是满足连续方程的精确解的话，那么RHS（右端项）是零，从而修正场也是零（即，如果满足连续方程话就无需再修正了）。因此，正是当前场的质量守恒误差驱动着修正场，单元面处的质量流量和质量流量修正为
m˙e∗=ρve∗⋅Se=ρue∗Sex=ρue∗Δyem˙w∗=ρvw∗⋅Sw=ρuw∗Swx=−ρuw∗Δyw\begin{aligned} \dot m_e^* &= \rho \bold v_e^* \cdot \bold S_e=\rho u_e^* S_e^x = \rho u_e^* \Delta y_e \\ \dot m_w^* &= \rho \bold v_w^* \cdot \bold S_w=\rho u_w^* S_w^x = -\rho u_w^* \Delta y_w \end{aligned} m˙e∗m˙w∗=ρve∗⋅Se=ρue∗Sex=ρue∗Δye=ρvw∗⋅Sw=ρuw∗Swx=−ρuw∗Δyw
和
m˙e′=ρve′⋅Se=ρue′Sex=ρue′Δyem˙w′=ρvw′⋅Sw=ρuw′Swx=−ρuw′Δyw\begin{aligned} \dot m_e' &= \rho \bold v_e' \cdot \bold S_e=\rho u_e' S_e^x = \rho u_e' \Delta y_e \\ \dot m_w' &= \rho \bold v_w' \cdot \bold S_w=\rho u_w' S_w^x = -\rho u_w' \Delta y_w \end{aligned} m˙e′m˙w′=ρve′⋅Se=ρue′Sex=ρue′Δye=ρvw′⋅Sw=ρuw′Swx=−ρuw′Δyw
上式中用到了Sex=ΔyeS_e^x=\Delta y_eSex=Δye和Swx=−ΔywS_w^x=-\Delta y_wSwx=−Δyw。

压力场并没有出现在修正形式的连续方程m˙e′+m˙w′=−m˙e∗−m˙w∗\dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^*m˙e′+m˙w′=−m˙e∗−m˙w∗中，为了将其引入到式中，需要使用动量方程的离散形式。该过程先把离散动量方程aeuue+∑f∼NB(e)afuuf=beu−Ve(∂p∂x)ea_e^u u_e+\displaystyle\sum_{f\sim NB(e)}a_f^u u_f=b_e^u-V_e\left(\dfrac{\partial p}{\partial x}\right)_eaeuue+f∼NB(e)∑afuuf=beu−Ve(∂x∂p)e写成更加紧凑的形式
ue+He(ue)=Beu−Deu(∂p∂x)eu_e+H_e(u_e) = B_e^u - D_e^u\left(\frac{\partial p}{\partial x}\right)_e ue+He(ue)=Beu−Deu(∂x∂p)e
其中
He(ue)=∑f∼NB(e)afuaeuufBeu=beuaeuDeu=VeaeuH_e(u_e) =\sum_{f\sim NB(e)}\frac{a_f^u}{a_e^u} u_f \quad\quad B_e^u = \frac{b_e^u}{a_e^u} \quad\quad D_e^u = \frac{V_e}{a_e^u} He(ue)=f∼NB(e)∑aeuafuufBeu=aeubeuDeu=aeuVe
对于计算的速度场ue∗u_e^*ue∗，上式变为
ue∗+He(ue∗)=Beu−Deu(∂p(n)∂x)eu_e^*+H_e(u_e^*) = B_e^u - D_e^u\left(\frac{\partial p^{(n)}}{\partial x}\right)_e ue∗+He(ue∗)=Beu−Deu(∂x∂p(n))e
用精确速度场表示的动量方程ue+He(ue)=Beu−Deu(∂p∂x)eu_e+H_e(u_e) = B_e^u - D_e^u\left(\frac{\partial p}{\partial x}\right)_eue+He(ue)=Beu−Deu(∂x∂p)e，减去，算得速度场所表示的动量方程ue∗+He(ue∗)=Beu−Deu(∂p(n)∂x)eu_e^*+H_e(u_e^*) = B_e^u - D_e^u\left(\frac{\partial p^{(n)}}{\partial x}\right)_eue∗+He(ue∗)=Beu−Deu(∂x∂p(n))e，可得到修正场的方程为
ue′+He(u′)‾=−Deu(∂p′∂x)eu_e'+\underline{H_e(u')} = - D_e^u\left(\frac{\partial p'}{\partial x}\right)_e ue′+He(u′)=−Deu(∂x∂p′)e
同样也可以推导出www面的修正方程
uw′+Hw(u′)‾=−Dwu(∂p′∂x)wu_w'+\underline{H_w(u')} = - D_w^u\left(\frac{\partial p'}{\partial x}\right)_w uw′+Hw(u′)=−Dwu(∂x∂p′)w
将m˙e′=ρue′Δye\dot m_e' = \rho u_e' \Delta y_em˙e′=ρue′Δye和m˙w′=−ρuw′Δyw\dot m_w' = -\rho u_w' \Delta y_wm˙w′=−ρuw′Δyw代入到连续方程m˙e′+m˙w′=−m˙e∗−m˙w∗\dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^*m˙e′+m˙w′=−m˙e∗−m˙w∗，得
ρeue′Δye+(−ρwuw′Δyw)=−(m˙e∗+m˙w∗)\rho_e u_e' \Delta y_e+(-\rho_w u_w' \Delta y_w)=-(\dot m_e^*+\dot m_w^*) ρeue′Δye+(−ρwuw′Δyw)=−(m˙e∗+m˙w∗)
然后，把上上式和上上上式的ue′u_e'ue′和uw′u_w'uw′的离散形式代入上式，可得到压力修正项的方程
ρe[−He(u′)‾−Deu(∂p′∂x)e]Δye−ρw[−Hw(u′)‾−Dwu(∂p′∂x)w]Δyw=−(m˙e∗+m˙w∗)\begin{aligned} &\rho_e \left[ -\underline{H_e(u')} - D_e^u\left(\frac{\partial p'}{\partial x}\right)_e \right] \Delta y_e \\ &-\rho_w \left[ -\underline{H_w(u')} - D_w^u\left(\frac{\partial p'}{\partial x}\right)_w \right] \Delta y_w=-(\dot m_e^*+\dot m_w^*) \end{aligned} ρe[−He(u′)−Deu(∂x∂p′)e]Δye−ρw[−Hw(u′)−Dwu(∂x∂p′)w]Δyw=−(m˙e∗+m˙w∗)
在该方程中，压力场以扩散项的形式出现，离散后变为
ρe[−He(u′)‾−Deu(pE′−pC′Δx)e]Δye+ρw[−Hw(u′)‾−Dwu(pC′−pW′Δx)w](−Δyw)=−(m˙e∗+m˙w∗)\begin{aligned} &\rho_e \left[ -\underline{H_e(u')} - D_e^u\left(\frac{p_E'-p_C'}{\Delta x}\right)_e \right] \Delta y_e \\ &+\rho_w \left[ -\underline{H_w(u')} - D_w^u\left(\frac{p_C'-p_W'}{\Delta x}\right)_w \right] (-\Delta y_w)=-(\dot m_e^*+\dot m_w^*) \end{aligned} ρe[−He(u′)−Deu(ΔxpE′−pC′)e]Δye+ρw[−Hw(u′)−Dwu(ΔxpC′−pW′)w](−Δyw)=−(m˙e∗+m˙w∗)
或
−ρeDeu(ΔyeΔxe)(pE′−pC′)−ρwDwu(−ΔywΔxw)(pC′−pW′)=−(m˙e∗+m˙w∗)+ρeHe(u′)‾Δye+ρwHw(u′)‾(−Δyw)\begin{aligned} & -\rho_e D_e^u\left(\frac{\Delta y_e}{\Delta x_e}\right)(p_E'-p_C') -\rho_w D_w^u\left(-\frac{\Delta y_w}{\Delta x_w}\right) (p_C'-p_W')\\ &=-(\dot m_e^*+\dot m_w^*)+\rho_e \underline{H_e(u')}\Delta y_e+\rho_w \underline{H_w(u')}(-\Delta y_w) \end{aligned} −ρeDeu(ΔxeΔye)(pE′−pC′)−ρwDwu(−ΔxwΔyw)(pC′−pW′)=−(m˙e∗+m˙w∗)+ρeHe(u′)Δye+ρwHw(u′)(−Δyw)
整理后的压力修正方程形式为
aCp′pC′+aEp′pE′+aWp′pW′=bCp′a_C^{p'}p_C'+a_E^{p'}p_E'+a_W^{p'}p_W'=b_C^{p'} aCp′pC′+aEp′pE′+aWp′pW′=bCp′
其中
aEp′=−ρeDeuΔyeΔxeaWp′=−ρwDwuΔywΔxwaCp′=−(aEp′+aWp′)bCp′=−(m˙e∗+m˙w∗)+[ρeΔyeHe(u′)−ρwΔywHw(u′)]‾\begin{aligned} a_E^{p'} &= -\frac{\rho_e D_e^u\Delta y_e}{\Delta x_e} \\ a_W^{p'} &= -\frac{\rho_w D_w^u\Delta y_w}{\Delta x_w}\\ a_C^{p'} &= -(a_E^{p'}+a_W^{p'}) \\ b_C^{p'} &= -(\dot m_e^*+\dot m_w^*)+ \underline{[\rho_e \Delta y_e H_e(u')-\rho_w \Delta y_w H_w(u')]} \end{aligned} aEp′aWp′aCp′bCp′=−ΔxeρeDeuΔye=−ΔxwρwDwuΔyw=−(aEp′+aWp′)=−(m˙e∗+m˙w∗)+[ρeΔyeHe(u′)−ρwΔywHw(u′)]
注意，对于上式下划线所标记的项，在达到稳态收敛的情况下，其值将变为零。因此，它们对于最终的解将没有影响。对于这些项的不同近似方法将会衍生不同的算法，这将在下面详细展开。在最初的SIMPLE算法中，这些项是简单地直接忽略掉了。此外，对于一维问题，且面积Δy\Delta yΔy为定值的情况，可直接把它设成111并从方程中刨掉即可。

2.6 交错网格上的SIMPLE算法

使用动量方程和压力修正方程，可以获得流动问题的解。在SIMPLE算法中的求解过程是，在每步迭代中，交替地生成压力和速度场，并先后满足动量方程和连续方程，然后再逐步逼近最终解（连续方程和动量方程这两个方程都满足的最终解）。该流程并非是同步进行的，这种求解方程的法子常被称为分离式方法（segregated approach）（实际上还有一种把压力和速度直接耦合起来同步求解的法子呢，不过比较耗内存了）。该分离式的SIMPLE算法的计算流程汇总如下：

从某个猜测的压力场p(n)p^{(n)}p(n)和速度场u(n)u^{(n)}u(n)开始；
求解动量方程来获得新的速度场uf∗u_f^*uf∗；
aeuue∗+∑f∼NB(e)afuuf∗=beu−Ve(∂p(n)∂x)eawuuw∗+∑f∼NB(e)afuuf∗=bwu−Vw(∂p(n)∂x)wa_e^u u_e^*+\sum_{f\sim NB(e)}a_f^u u_f^*= b_e^u-V_e\left(\frac{\partial p^{(n)}}{\partial x}\right)_e \\ a_w^u u_w^*+\sum_{f\sim NB(e)}a_f^u u_f^*= b_w^u-V_w\left(\frac{\partial p^{(n)}}{\partial x}\right)_w aeuue∗+f∼NB(e)∑afuuf∗=beu−Ve(∂x∂p(n))eawuuw∗+f∼NB(e)∑afuuf∗=bwu−Vw(∂x∂p(n))w
使用满足动量方程的速度场uf∗u_f^*uf∗来更新质量流量m˙f∗\dot m_f^*m˙f∗；
m˙f∗=ρvf∗⋅Sf=ρfuf∗Sfx\dot m_f^* = \rho \bold v_f^* \cdot \bold S_f=\rho_f u_f^* S_f^x m˙f∗=ρvf∗⋅Sf=ρfuf∗Sfx
使用新的质量流量m˙f∗\dot m_f^*m˙f∗求解压力修正方程来获取压力修正场p′p'p′；
aCp′pC′+aEp′pE′+aWp′pW′=bCp′a_C^{p'}p_C'+a_E^{p'}p_E'+a_W^{p'}p_W'=b_C^{p'} aCp′pC′+aEp′pE′+aWp′pW′=bCp′
更新压力场和速度场以获取满足连续方程的场
uf∗∗=uf∗+uf′uf′=−Dfu(∂p′∂x)fpC∗=pC(n)+pC′m˙f∗∗=m˙f∗+m˙f′m˙f′=−ρfDfuΔyf(∂p′∂x)f\begin{aligned} u_f^{**} &= u_f^* + u_f' \quad\quad u_f'=-D_f^u\left(\frac{\partial p'}{\partial x}\right)_f \\ p_C^* &=p_C^{(n)} + p_C' \\ \dot m_f^{**} &= \dot m_f^* + \dot m_f' \quad\quad \dot m_f'=-\rho_f D_f^u \Delta y_f \left(\frac{\partial p'}{\partial x}\right)_f \end{aligned} uf∗∗pC∗m˙f∗∗=uf∗+uf′uf′=−Dfu(∂x∂p′)f=pC(n)+pC′=m˙f∗+m˙f′m˙f′=−ρfDfuΔyf(∂x∂p′)f
令u(n)=u∗∗u^{(n)}=u^{**}u(n)=u∗∗，p(n)=p∗p^{(n)}=p^{*}p(n)=p∗；
返回第2步，并重复该流程直到收敛。

理论是晦涩的，而实例则是生动的，用个简单的小例子来展示下SIMPLE的强大是再好不过的了。

例1 管道网络（管网）中的流动

下图展示的是水力管网系统中的一部分，管道流动中的动量方程可以写成
m˙=ρuA=−D∇P\dot m=\rho u A = -D\nabla P m˙=ρuA=−D∇P
其中，DA=0.5D_A=0.5DA=0.5，DB=DF=0.4D_B=D_F=0.4DB=DF=0.4，DC=DE=0.3D_C=D_E=0.3DC=DE=0.3，DD=0.19D_D=0.19DD=0.19，DG=0.1875D_G=0.1875DG=0.1875，DH=0.35D_H=0.35DH=0.35。使用SIMPLE算法，计算系统中的未知质量流量和压力。

解

在该系统中，使用质量流量场作为变量，来替代速度场。这并不会产生什么问题，因为该例中的动量方程已经忽略掉了对流项和扩散项，仅保留了压力梯度项，即认为对流和扩散效应，与压差驱动力相比非常小，直接忽略了。

使用SIMPLE算法的求解要先用假设的压力场来求解一个初始速度场，然后再根据这个计算的速度场预测一个压力场，以满足连续方程。

该流程汇总如下：

由一个猜测的压力场开始；
使用给定的动量方程计算质量流量；
构造压力修正方程来满足连续性（质量守恒），并用其修正压力和速度场。

该例中无需迭代，因为其不含对流项所带来的非线性效应（如果是用SIMPLE来求解不可压缩流动，那是要迭代多次才能收敛的）。

step 1
给定待求解位置处的压力猜测值，如下：
p3(n)=300p6(n)=200p8(n)=120p_3^{(n)}=300 \quad\quad p_6^{(n)}=200 \quad\quad p_8^{(n)}=120 p3(n)=300p6(n)=200p8(n)=120
当然也可以使用其它值。

step 2
基于假设的压力值，使用动量方程计算不同的质量流量：
m˙A∗=DA(p1−p3(n))=0.5∗(400−300)=50m˙B∗=DB(p3(n)−p2)=0.4∗(300−350)=−20m˙C∗=DC(p4−p3(n))=0.3∗(50−300)=−75m˙D∗=DD(p3(n)−p6(n))=0.19∗(300−200)=19m˙E∗=DE(p6(n)−p5)=0.3∗(200−300)=−30m˙F∗=DF(p7−p6(n))=0.4∗(80−200)=−48m˙G∗=DG(p6(n)−p8(n))=0.1875∗(200−120)=15m˙H∗=DH(p9−p8(n))=0.35∗(200−120)=28\begin{aligned} & \dot m_A^*=D_A(p_1-p_3^{(n)})=0.5*(400-300)=50 \\ & \dot m_B^*=D_B(p_3^{(n)}-p_2)=0.4*(300-350)=-20 \\ & \dot m_C^*=D_C(p_4-p_3^{(n)})=0.3*(50-300)=-75 \\ & \dot m_D^*=D_D(p_3^{(n)}-p_6^{(n)})=0.19*(300-200)=19 \\ & \dot m_E^*=D_E(p_6^{(n)}-p_5)=0.3*(200-300)=-30 \\ & \dot m_F^*=D_F(p_7-p_6^{(n)})=0.4*(80-200)=-48 \\ & \dot m_G^*=D_G(p_6^{(n)}-p_8^{(n)})=0.1875*(200-120)=15 \\ & \dot m_H^*=D_H(p_9-p_8^{(n)})=0.35*(200-120)=28 \end{aligned} m˙A∗=DA(p1−p3(n))=0.5∗(400−300)=50m˙B∗=DB(p3(n)−p2)=0.4∗(300−350)=−20m˙C∗=DC(p4−p3(n))=0.3∗(50−300)=−75m˙D∗=DD(p3(n)−p6(n))=0.19∗(300−200)=19m˙E∗=DE(p6(n)−p5)=0.3∗(200−300)=−30m˙F∗=DF(p7−p6(n))=0.4∗(80−200)=−48m˙G∗=DG(p6(n)−p8(n))=0.1875∗(200−120)=15m˙H∗=DH(p9−p8(n))=0.35∗(200−120)=28

step 3
检查质量流量是否满足连续性，对所有内部节点计算∑∼km˙k∗\displaystyle\sum_{\sim k}\dot m_k^*∼k∑m˙k∗，得
Node 3:∑∼k(m˙k∗)=m˙B∗+m˙D∗−m˙A∗−m˙C∗=−20+19−50+75=24Node 6:∑∼k(m˙k∗)=m˙E∗+m˙G∗−m˙D∗−m˙F∗=−30+15−19+48=14Node 8:∑∼k(m˙k∗)=m˙I∗−m˙G∗−m˙H∗=50−15−28=7\begin{aligned} & \text{Node 3:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*=-20+19-50+75=24 \\ & \text{Node 6:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*=-30+15-19+48=14 \\ & \text{Node 8:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_I^*-\dot m_G^*-\dot m_H^*=50-15-28=7 \end{aligned} Node 3:∼k∑(m˙k∗)=m˙B∗+m˙D∗−m˙A∗−m˙C∗=−20+19−50+75=24Node 6:∼k∑(m˙k∗)=m˙E∗+m˙G∗−m˙D∗−m˙F∗=−30+15−19+48=14Node 8:∼k∑(m˙k∗)=m˙I∗−m˙G∗−m˙H∗=50−15−28=7
由于质量守恒并不满足，需要修正处理，压力修正方程推导如下：
∑∼k(m˙k∗+m˙k′)=0⇒∑∼k(m˙k′)=−∑∼k(m˙k∗)\sum_{\sim k}(\dot m_k^*+\dot m_k')=0 \Rightarrow \sum_{\sim k}(\dot m_k')=-\sum_{\sim k}(\dot m_k^*) ∼k∑(m˙k∗+m˙k′)=0⇒∼k∑(m˙k′)=−∼k∑(m˙k∗)
基于动量方程，把质量流量的修正用压力修正表示为
m˙A′=DA(−p3′)m˙B′=DB(p3′)m˙C′=DC(−p3′)m˙D′=DD(p3′−p6′)m˙E′=DE(p6′)m˙F′=DF(−p6′)m˙G′=DG(p6′−p8′)m˙H′=DH(−p8′)\begin{aligned} & \dot m'_A = D_A(-p'_3) \\ & \dot m_B'=D_B(p_3') \\ & \dot m_C'=D_C(-p_3') \\ & \dot m_D'=D_D(p_3'-p_6') \\ & \dot m_E'=D_E(p_6') \\ & \dot m_F'=D_F(-p_6') \\ & \dot m_G'=D_G(p_6'-p_8') \\ & \dot m_H'=D_H(-p_8') \end{aligned} m˙A′=DA(−p3′)m˙B′=DB(p3′)m˙C′=DC(−p3′)m˙D′=DD(p3′−p6′)m˙E′=DE(p6′)m˙F′=DF(−p6′)m˙G′=DG(p6′−p8′)m˙H′=DH(−p8′)
注意p1,2,4,5,7′p'_{1,2,4,5,7}p1,2,4,5,7′为0，因为这些压力值是已知的，本身就是精确值，所以无需修正。

节点3、6、83、6、83、6、8处的质量守恒方程（连续方程）为
m˙B′+m˙D′−m˙A′−m˙C′=−(m˙B∗+m˙D∗−m˙A∗−m˙C∗)m˙E′+m˙G′−m˙D′−m˙F′=−(m˙E∗+m˙G∗−m˙D∗−m˙F∗)−m˙G′−m˙H′=−(m˙I∗−m˙G∗−m˙H∗)\begin{aligned} \dot m_B'+\dot m_D'-\dot m_A'-\dot m_C' &= -(\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*) \\ \dot m_E'+\dot m_G'-\dot m_D'-\dot m_F' &= -(\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*)\\ -\dot m_G'-\dot m_H' &= -(\dot m_I^*-\dot m_G^*-\dot m_H^*) \end{aligned} m˙B′+m˙D′−m˙A′−m˙C′m˙E′+m˙G′−m˙D′−m˙F′−m˙G′−m˙H′=−(m˙B∗+m˙D∗−m˙A∗−m˙C∗)=−(m˙E∗+m˙G∗−m˙D∗−m˙F∗)=−(m˙I∗−m˙G∗−m˙H∗)
把压力修正表示的质量流量修正代入上式，得到压力修正方程
DB(p3′)+DD(p3′−p6′)−DA(−p3′)−DC(−p3′)=−(m˙B∗+m˙D∗−m˙A∗−m˙C∗)DE(p6′)+DG(p6′−p8′)−DD(p3′−p6′)−DF(−p6′)=−(m˙E∗+m˙G∗−m˙D∗−m˙F∗)−DG(p6′−p8′)−DH(−p8′)=−(m˙I∗−m˙G∗−m˙H∗)\begin{aligned} D_B(p_3')+D_D(p_3'-p_6')-D_A(-p'_3)-D_C(-p_3') &= -(\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*) \\ D_E(p_6')+D_G(p_6'-p_8')-D_D(p_3'-p_6')-D_F(-p_6') &= -(\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*)\\ -D_G(p_6'-p_8')-D_H(-p_8') &= -(\dot m_I^*-\dot m_G^*-\dot m_H^*) \end{aligned} DB(p3′)+DD(p3′−p6′)−DA(−p3′)−DC(−p3′)DE(p6′)+DG(p6′−p8′)−DD(p3′−p6′)−DF(−p6′)−DG(p6′−p8′)−DH(−p8′)=−(m˙B∗+m˙D∗−m˙A∗−m˙C∗)=−(m˙E∗+m˙G∗−m˙D∗−m˙F∗)=−(m˙I∗−m˙G∗−m˙H∗)
整理后，压力修正方程变为
1.39(p3′)−0.19(p6′)=−(m˙B∗+m˙D∗−m˙A∗−m˙C∗)1.0775(p6′)−0.1875(p8′)−0.19(p3′)=−(m˙E∗+m˙G∗−m˙D∗−m˙F∗)−0.1875(p6′)+0.5375(p8′)=−(m˙I∗−m˙G∗−m˙H∗)\begin{aligned} 1.39(p_3')-0.19(p_6') &= -(\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*) \\ 1.0775(p_6')-0.1875(p_8')-0.19(p_3') &= -(\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*)\\ -0.1875(p_6')+0.5375(p_8') &= -(\dot m_I^*-\dot m_G^*-\dot m_H^*) \end{aligned} 1.39(p3′)−0.19(p6′)1.0775(p6′)−0.1875(p8′)−0.19(p3′)−0.1875(p6′)+0.5375(p8′)=−(m˙B∗+m˙D∗−m˙A∗−m˙C∗)=−(m˙E∗+m˙G∗−m˙D∗−m˙F∗)=−(m˙I∗−m˙G∗−m˙H∗)
把算得的初始质量流量代入进去，得到修正方程为
1.39(p3′)−0.19(p6′)=−241.0775(p6′)−0.1875(p8′)−0.19(p3′)=−14−0.1875(p6′)+0.5375(p8′)=−7\begin{aligned} 1.39(p_3')-0.19(p_6') &= -24 \\ 1.0775(p_6')-0.1875(p_8')-0.19(p_3') &= -14 \\ -0.1875(p_6')+0.5375(p_8') &= -7 \end{aligned} 1.39(p3′)−0.19(p6′)1.0775(p6′)−0.1875(p8′)−0.19(p3′)−0.1875(p6′)+0.5375(p8′)=−24=−14=−7
求解压力修正方程，得
{p3′=−20p6′=−20p8′=−20\begin{cases} p'_3=-20 \\ p'_6=-20 \\ p'_8=-20 \end{cases} ⎩⎪⎨⎪⎧p3′=−20p6′=−20p8′=−20
使用算得的修正压力，更新速度场和压力场，以产生满足质量守恒的速度场。修正后的质量流量为：
m˙A∗∗=m˙A∗+m˙A′=m˙A∗+DA(−p3′)=50+0.5(−(−20))=60m˙B∗∗=m˙B∗+m˙B′=m˙B∗+DB(p3′)=−20+0.4(−20)=−28m˙C∗∗=m˙C∗+m˙C′=m˙C∗+DC(−p3′)=−75+0.3(−(−20))=−69m˙D∗∗=m˙D∗+m˙D′=m˙D∗+DD(p3′−p6′)=19+0.19(−20−(−20))=19m˙E∗∗=m˙E∗+m˙E′=m˙E∗+DE(p6′)=−30+0.3(−20)=−36m˙F∗∗=m˙F∗+m˙F′=m˙F∗+DF(−p6′)=−48+0.4(−(−20))=−40m˙G∗∗=m˙G∗+m˙G′=m˙G∗+DG(p6′−p8′)=15+0.1875(−20−(−20))=15m˙H∗∗=m˙H∗+m˙H′=m˙H∗+DH(−p8′)=28+0.35(−(−20))=35\begin{aligned} & \dot m_A^{**}=\dot m_A^*+\dot m'_A = \dot m_A^*+D_A(-p'_3)=50+0.5(-(-20))=60 \\ & \dot m_B^{**}=\dot m_B^*+\dot m_B'=\dot m_B^*+D_B(p_3')=-20+0.4(-20)=-28 \\ & \dot m_C^{**}=\dot m_C^*+\dot m_C'=\dot m_C^*+D_C(-p_3')=-75+0.3(-(-20))=-69 \\ & \dot m_D^{**}=\dot m_D^*+\dot m_D'=\dot m_D^*+D_D(p_3'-p_6')=19+0.19(-20-(-20))=19 \\ & \dot m_E^{**}=\dot m_E^*+\dot m_E'=\dot m_E^*+D_E(p_6')=-30+0.3(-20)=-36 \\ & \dot m_F^{**}=\dot m_F^*+\dot m_F'=\dot m_F^*+D_F(-p_6')=-48+0.4(-(-20))=-40 \\ & \dot m_G^{**}=\dot m_G^*+\dot m_G'=\dot m_G^*+D_G(p_6'-p_8')=15+0.1875(-20-(-20))=15 \\ & \dot m_H^{**}=\dot m_H^*+\dot m_H'=\dot m_H^*+D_H(-p_8')=28+0.35(-(-20))=35 \end{aligned} m˙A∗∗=m˙A∗+m˙A′=m˙A∗+DA(−p3′)=50+0.5(−(−20))=60m˙B∗∗=m˙B∗+m˙B′=m˙B∗+DB(p3′)=−20+0.4(−20)=−28m˙C∗∗=m˙C∗+m˙C′=m˙C∗+DC(−p3′)=−75+0.3(−(−20))=−69m˙D∗∗=m˙D∗+m˙D′=m˙D∗+DD(p3′−p6′)=19+0.19(−20−(−20))=19m˙E∗∗=m˙E∗+m˙E′=m˙E∗+DE(p6′)=−30+0.3(−20)=−36m˙F∗∗=m˙F∗+m˙F′=m˙F∗+DF(−p6′)=−48+0.4(−(−20))=−40m˙G∗∗=m˙G∗+m˙G′=m˙G∗+DG(p6′−p8′)=15+0.1875(−20−(−20))=15m˙H∗∗=m˙H∗+m˙H′=m˙H∗+DH(−p8′)=28+0.35(−(−20))=35
压力修正为
p3∗=p3(n)+p3′=300+(−20)=280p6∗=p6(n)+p6′=200+(−20)=180p8∗=p8(n)+p8′=120+(−20)=100\begin{aligned} & p_3^*=p_3^{(n)}+p_3'=300+(-20)=280 \\ & p_6^*=p_6^{(n)}+p_6'=200+(-20)=180 \\ & p_8^*=p_8^{(n)}+p_8'=120+(-20)=100 \\ \end{aligned} p3∗=p3(n)+p3′=300+(−20)=280p6∗=p6(n)+p6′=200+(−20)=180p8∗=p8(n)+p8′=120+(−20)=100
把修正后的值作为新的猜测值，即
p3(n)=p3∗=280p6(n)=p6∗=180p8(n)=p8∗=100\begin{aligned} & p_3^{(n)}=p_3^*=280 \\ & p_6^{(n)}=p_6^*=180 \\ & p_8^{(n)}=p_8^*=100 \end{aligned} p3(n)=p3∗=280p6(n)=p6∗=180p8(n)=p8∗=100
继续上述流程，使用动量方程计算不同的质量流量为
m˙A∗=DA(p1−p3(n))=0.5∗(400−280)=60m˙B∗=DB(p3(n)−p2)=0.4∗(280−350)=−28m˙C∗=DC(p4−p3(n))=0.3∗(50−280)=−69m˙D∗=DD(p3(n)−p6(n))=0.19∗(280−180)=19m˙E∗=DE(p6(n)−p5)=0.3∗(180−300)=−36m˙F∗=DF(p7−p6(n))=0.4∗(80−180)=−40m˙G∗=DG(p6(n)−p8(n))=0.1875∗(180−100)=15m˙H∗=DH(p9−p8(n))=0.35(200−100)=35\begin{aligned} & \dot m_A^*=D_A(p_1-p_3^{(n)})=0.5*(400-280)=60 \\ & \dot m_B^*=D_B(p_3^{(n)}-p_2)=0.4*(280-350)=-28 \\ & \dot m_C^*=D_C(p_4-p_3^{(n)})=0.3*(50-280)=-69 \\ & \dot m_D^*=D_D(p_3^{(n)}-p_6^{(n)})=0.19*(280-180)=19 \\ & \dot m_E^*=D_E(p_6^{(n)}-p_5)=0.3*(180-300)=-36 \\ & \dot m_F^*=D_F(p_7-p_6^{(n)})=0.4*(80-180)=-40 \\ & \dot m_G^*=D_G(p_6^{(n)}-p_8^{(n)})=0.1875*(180-100)=15 \\ & \dot m_H^*=D_H(p_9-p_8^{(n)})=0.35(200-100)=35 \end{aligned} m˙A∗=DA(p1−p3(n))=0.5∗(400−280)=60m˙B∗=DB(p3(n)−p2)=0.4∗(280−350)=−28m˙C∗=DC(p4−p3(n))=0.3∗(50−280)=−69m˙D∗=DD(p3(n)−p6(n))=0.19∗(280−180)=19m˙E∗=DE(p6(n)−p5)=0.3∗(180−300)=−36m˙F∗=DF(p7−p6(n))=0.4∗(80−180)=−40m˙G∗=DG(p6(n)−p8(n))=0.1875∗(180−100)=15m˙H∗=DH(p9−p8(n))=0.35(200−100)=35
节点3、6、83、6、83、6、8处的不平衡质量流量为
Node 3:∑∼k(m˙k∗)=m˙B∗+m˙D∗−m˙A∗−m˙C∗=−28+19−60+69=0Node 6:∑∼k(m˙k∗)=m˙E∗+m˙G∗−m˙D∗−m˙F∗=−36+15−19+40=0Node 8:∑∼k(m˙k∗)=m˙I∗−m˙G∗−m˙H∗=50−15−35=0\begin{aligned} & \text{Node 3:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*=-28+19-60+69=0 \\ & \text{Node 6:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*=-36+15-19+40=0 \\ & \text{Node 8:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_I^*-\dot m_G^*-\dot m_H^*=50-15-35=0 \end{aligned} Node 3:∼k∑(m˙k∗)=m˙B∗+m˙D∗−m˙A∗−m˙C∗=−28+19−60+69=0Node 6:∼k∑(m˙k∗)=m˙E∗+m˙G∗−m˙D∗−m˙F∗=−36+15−19+40=0Node 8:∼k∑(m˙k∗)=m˙I∗−m˙G∗−m˙H∗=50−15−35=0
压力修正方程为
1.39(p3′)−0.19(p6′)=01.0775(p6′)−0.1875(p8′)−0.19(p3′)=0−0.1875(p6′)+0.5375(p8′)=0\begin{aligned} 1.39(p_3')-0.19(p_6') &= 0 \\ 1.0775(p_6')-0.1875(p_8')-0.19(p_3') &= 0 \\ -0.1875(p_6')+0.5375(p_8') &= 0 \end{aligned} 1.39(p3′)−0.19(p6′)1.0775(p6′)−0.1875(p8′)−0.19(p3′)−0.1875(p6′)+0.5375(p8′)=0=0=0
解得压力修正值为
{p3′=0p6′=0p8′=0\begin{cases} p'_3=0 \\ p'_6=0 \\ p'_8=0 \end{cases} ⎩⎪⎨⎪⎧p3′=0p6′=0p8′=0
可见，只用一步迭代就获得了收敛解（其实第2次迭代的时候算出来各节点的不平衡质量流量为零就无需再往后计算压力修正值了），这是因为该例子较为简单，不含非线性效应的缘故。实际若是包含对流项、扩散项和压力梯度项的动量方程的话，SIMPLE算法的求解时需要迭代很多次才能收敛的。

2.7 二维交错笛卡尔网格上的压力修正方程

在二维笛卡尔网格中，将使用三套网格系统来构建交错网格体系。第一套是存储速度分量uuu的，第二套是存储速度分量vvv的，第三套是存储压力和其它变量的，如下图所示。

前面展示的一维区域内压力修正方程的推导方法，可以很容易地扩展到多维情况。对于下图所示的单元CCC，其压力修正方程为

aCp′pC′+aEp′pE′+aWp′pW′+aNp′pN′+aSp′pS′=bCp′a_C^{p'}p'_C+a_E^{p'}p'_E+a_W^{p'}p'_W+a_N^{p'}p'_N+a_S^{p'}p'_S=b_C^{p'} aCp′pC′+aEp′pE′+aWp′pW′+aNp′pN′+aSp′pS′=bCp′
其中
aEp′=−ρeDeuΔyCδxeaWp′=−ρwDwuΔyCδxwaNp′=−ρnDnvΔxCδynaSp′=−ρsDsvΔxCδysaCp′=−(aEp′+aWp′+aNp′+aSp′)bCp′=−(m˙e∗+m˙w∗+m˙n∗+m˙s∗)\begin{aligned} a_E^{p'} &= -\frac{\rho_e D_e^u \Delta y_C}{\delta x_e} \quad\quad a_W^{p'} = -\frac{\rho_w D_w^u \Delta y_C}{\delta x_w} \\ a_N^{p'} &= -\frac{\rho_n D_n^v \Delta x_C}{\delta y_n} \quad\quad a_S^{p'} = -\frac{\rho_s D_s^v \Delta x_C}{\delta y_s} \\ a_C^{p'} &= -(a_E^{p'}+a_W^{p'}+a_N^{p'}+a_S^{p'})\\ b_C^{p'} &= -(\dot m^*_e+\dot m^*_w+\dot m^*_n+\dot m^*_s) \end{aligned} aEp′aNp′aCp′bCp′=−δxeρeDeuΔyCaWp′=−δxwρwDwuΔyC=−δynρnDnvΔxCaSp′=−δysρsDsvΔxC=−(aEp′+aWp′+aNp′+aSp′)=−(m˙e∗+m˙w∗+m˙n∗+m˙s∗)

例2

在下图所示的二维问题中，变量uw=50u_w=50uw=50，vs=20v_s=20vs=20，pN=0p_N=0pN=0，pE=10p_E=10pE=10。

流动是定常的，密度为常数且等于1，则ueu_eue和vnv_nvn的动量方程为（同样不考虑对流项和扩散项，只考虑压力梯度项）
ue=−de(pE−pC)vn=−dn(pN−pC)\begin{aligned} u_e &= -d_e(p_E-p_C) \\ v_n &= -d_n(p_N-p_C) \end{aligned} uevn=−de(pE−pC)=−dn(pN−pC)
其中常量系数de=1d_e=1de=1、dn=0.25d_n=0.25dn=0.25。单元的Δx=Δy=1\Delta x=\Delta y=1Δx=Δy=1。

使用SIMPLE算法计算ueu_eue、vnv_nvn和pCp_CpC的值。

解

对于要寻找的CCC处的压力赋予个猜测值，假设
pC(n)=100p_C^{(n)}=100 pC(n)=100
其他猜测值也都可以的。

基于假设的压力值，使用动量方程计算不同的速度
ue∗=−de(pE−pC(n))=−1∗(10−100)=90vn∗=−dn(pN−pC(n))=−0.25∗(0−100)=25\begin{aligned} u_e^* &= -d_e(p_E-p_C^{(n)})=-1*(10-100) =90 \\ v_n^* &= -d_n(p_N-p_C^{(n)})=-0.25*(0-100)=25 \end{aligned} ue∗vn∗=−de(pE−pC(n))=−1∗(10−100)=90=−dn(pN−pC(n))=−0.25∗(0−100)=25
由于密度是均匀的且等于111，而且有Δx=Δy=1\Delta x=\Delta y=1Δx=Δy=1，所以质量流量
m˙e∗=ue∗=90m˙n∗=vn∗=25\begin{aligned} \dot m_e^* &= u_e^* =90 \\ \dot m_n^* &= v_n^*=25 \end{aligned} m˙e∗m˙n∗=ue∗=90=vn∗=25
检查在单元CCC上的质量流量是否满足连续性
∑∼f(m˙f∗)=m˙e∗−m˙w∗+m˙n∗−m˙s∗=90−50+25−20=45\sum_{\sim f}(\dot m_f^*)=\dot m_e^*-\dot m_w^*+\dot m_n^*-\dot m_s^*=90-50+25-20=45 ∼f∑(m˙f∗)=m˙e∗−m˙w∗+m˙n∗−m˙s∗=90−50+25−20=45
在上面的质量守恒方程中，m˙w∗\dot m_w^*m˙w∗和m˙s∗\dot m_s^*m˙s∗是负的，而且是直接显式使用的，因为它们的值是确定的。因为质量守恒不满足，所以需要修正流场，要推导压力修正方程，其推导过程如下：
∑∼f(m˙f∗+m˙f′)=0⇒∑∼f(m˙f′)=−∑∼f(m˙f∗)\sum_{\sim f}(\dot m_f^*+\dot m_f')=0 \Rightarrow \sum_{\sim f}(\dot m_f')=-\sum_{\sim f}(\dot m_f^*) ∼f∑(m˙f∗+m˙f′)=0⇒∼f∑(m˙f′)=−∼f∑(m˙f∗)
再把质量流量的修正值用压力修正值来表示（从动量方程中推出），即
m˙e′=ue′=depC′=pC′m˙n′=vn′=dnpC′=0.25pC′\begin{aligned} \dot m_e' &= u_e'=d_e p_C'=p_C' \\ \dot m_n' &= v_n'=d_n p_C'=0.25p_C' \end{aligned} m˙e′m˙n′=ue′=depC′=pC′=vn′=dnpC′=0.25pC′
这样，得到压力修正方程为
m˙e′+m˙n′=−∑∼f(m˙f∗)⇒1.25pC′=−45⇒pC′=−36\dot m_e' +\dot m_n' =-\sum_{\sim f}(\dot m_f^*) \Rightarrow 1.25p_C'=-45 \Rightarrow p_C'=-36 m˙e′+m˙n′=−∼f∑(m˙f∗)⇒1.25pC′=−45⇒pC′=−36
获得了压力修正值后，可用其来修正质量流量和压力场，使得流场满足守恒条件，即
m˙e∗∗=m˙e∗+m˙e′=m˙e∗+pC′=90−36=54m˙n∗∗=m˙n∗+m˙n′=m˙n∗+0.25pC′=25−0.25∗36=16pC∗∗=pC(n)+pC′=100−36=64\begin{aligned} \dot m_e^{**} &=\dot m_e^*+\dot m_e'=\dot m_e^*+p_C'= 90-36=54 \\ \dot m_n^{**} &=\dot m_n^*+\dot m_n'=\dot m_n^*+0.25p_C'= 25-0.25*36=16 \\ p_C^{**} &=p_C^{(n)}+p_C'=100-36=64 \end{aligned} m˙e∗∗m˙n∗∗pC∗∗=m˙e∗+m˙e′=m˙e∗+pC′=90−36=54=m˙n∗+m˙n′=m˙n∗+0.25pC′=25−0.25∗36=16=pC(n)+pC′=100−36=64
用修正后的压力场和流量场作为新的猜测值，继续上述过程，即
m˙e∗=m˙e∗∗=54m˙n∗=m˙n∗∗=16pC(n)=pC∗∗=64\begin{aligned} \dot m_e^{*} &=\dot m_e^{**}=54 \\ \dot m_n^{*} &=\dot m_n^{**}=16 \\ p_C^{(n)} &=p_C^{**}=64 \end{aligned} m˙e∗m˙n∗pC(n)=m˙e∗∗=54=m˙n∗∗=16=pC∗∗=64
再次检查在单元CCC上的质量流量是否满足连续性
∑∼f(m˙f∗)=m˙e∗−m˙w∗+m˙n∗−m˙s∗=54−50+16−20=0\sum_{\sim f}(\dot m_f^*)=\dot m_e^*-\dot m_w^*+\dot m_n^*-\dot m_s^*=54-50+16-20=0 ∼f∑(m˙f∗)=m˙e∗−m˙w∗+m˙n∗−m˙s∗=54−50+16−20=0
再次计算压力修正方程
1.25pC′=−∑∼f(m˙f∗)=01.25p_C'=-\sum_{\sim f}(\dot m_f^*)=0 1.25pC′=−∼f∑(m˙f∗)=0
算得压力修正值为
pC′=0p_C'=0 pC′=0
即，无需修正，也就是说，一个迭代步就获得了收敛解
ue=m˙e∗=54vn=m˙n∗=16pC=pC(n)=64\begin{aligned} u_e&=\dot m_e^{*} =54 \\ v_n&=\dot m_n^{*} =16 \\ p_C&=p_C^{(n)} =64 \end{aligned} uevnpC=m˙e∗=54=m˙n∗=16=pC(n)=64

2.8 三维交错笛卡尔网格上的压力修正方程

不再做详细推导和完整展示，直接给出如上图所示的三维交错笛卡尔网格上的压力修正方程，其中uuu、vvv、www速度分量分别存储在(e,w)(e,w)(e,w)、(n,s)(n,s)(n,s)、(t,b)(t,b)(t,b)单元面上，压力修正方程为
aCp′pC′+aEp′pE′+aWp′pW′+aNp′pN′+aSp′pS′+aTp′pT′+aBp′pB′=bCp′a_C^{p'}p'_C+a_E^{p'}p'_E+a_W^{p'}p'_W+a_N^{p'}p'_N+a_S^{p'}p'_S+a_T^{p'}p'_T+a_B^{p'}p'_B=b_C^{p'} aCp′pC′+aEp′pE′+aWp′pW′+aNp′pN′+aSp′pS′+aTp′pT′+aBp′pB′=bCp′
其中
aEp′=−ρeDeuΔyeΔzeδxeaWp′=−ρwDwuΔywΔzwδxwaNp′=−ρnDnvΔxnΔznδynaSp′=−ρsDsvΔxsΔzsδysaTp′=−ρtDtwΔxtΔytδztaBp′=−ρbDbwΔxbΔybδzbaCp′=−(aEp′+aWp′+aNp′+aSp′+aTp′+aBp′)bCp′=−(m˙e∗+m˙w∗+m˙n∗+m˙s∗+m˙t∗+m˙b∗)\begin{aligned} a_E^{p'} &= -\frac{\rho_e D_e^u \Delta y_e \Delta z_e}{\delta x_e} \quad\quad a_W^{p'} = -\frac{\rho_w D_w^u \Delta y_w \Delta z_w}{\delta x_w} \\ a_N^{p'} &= -\frac{\rho_n D_n^v \Delta x_n \Delta z_n}{\delta y_n} \quad\quad a_S^{p'} = -\frac{\rho_s D_s^v \Delta x_s \Delta z_s}{\delta y_s} \\ a_T^{p'} &= -\frac{\rho_t D_t^w \Delta x_t \Delta y_t}{\delta z_t} \quad\quad a_B^{p'} = -\frac{\rho_b D_b^w \Delta x_b \Delta y_b}{\delta z_b} \\ a_C^{p'} &= -(a_E^{p'}+a_W^{p'}+a_N^{p'}+a_S^{p'}+a_T^{p'}+a_B^{p'})\\ b_C^{p'} &= -(\dot m^*_e+\dot m^*_w+\dot m^*_n+\dot m^*_s+\dot m^*_t+\dot m^*_b) \end{aligned} aEp′aNp′aTp′aCp′bCp′=−δxeρeDeuΔyeΔzeaWp′=−δxwρwDwuΔywΔzw=−δynρnDnvΔxnΔznaSp′=−δysρsDsvΔxsΔzs=−δztρtDtwΔxtΔytaBp′=−δzbρbDbwΔxbΔyb=−(aEp′+aWp′+aNp′+aSp′+aTp′+aBp′)=−(m˙e∗+m˙w∗+m˙n∗+m˙s∗+m˙t∗+m˙b∗)

FVM in CFD 学习笔记_第15章_流动计算：不可压缩流动_1_交错网格上的SIMPLE算法相关推荐

FVM in CFD 学习笔记_第9章_梯度计算
学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - ...
一文弄懂元学习 (Meta Learing)（附代码实战）《繁凡的深度学习笔记》第 15 章元学习详解（上）万字中文综述
<繁凡的深度学习笔记>第 15 章元学习详解 (上)万字中文综述(DL笔记整理系列) 3043331995@qq.com https://fanfansann.blog.csdn.net ...
FVM in CFD 学习笔记_第11章_对流项离散
学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - ...
FVM in CFD 学习笔记_第6章_有限体积网格
学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - ...
FVM in CFD 学习笔记_第7章_OpenFOAM和uFVM中的有限体积网格
学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - ...
FVM in CFD 学习笔记_第12章_高分辨率格式
学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - ...
读书笔记_代码大全_第14章_组织直线型代码_第15章_使用条件语句
组织直线型代码 + 使用条件语句希望我的读书笔试能带你翻过18页的书 http://www.cnblogs.com/jerry19880126/ <代码大全>第14章和第15章的内容比较 ...
学习笔记：Java 并发编程①_基础知识入门
若文章内容或图片失效,请留言反馈. 部分素材来自网络,若不小心影响到您的利益,请联系博主删除. 视频链接:https://www.bilibili.com/video/av81461839 视频下载: ...
学习笔记：Java 并发编程②_管程
若文章内容或图片失效,请留言反馈. 部分素材来自网络,若不小心影响到您的利益,请联系博主删除. 视频链接:https://www.bilibili.com/video/av81461839 配套资料: ...

FVM in CFD 学习笔记_第15章_流动计算：不可压缩流动_1_交错网格上的SIMPLE算法