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

前面幾章講解的關於變量 $\phi$的通常輸運方程的離散和求解流程，均是創建在事先已知速度場的前提下。可是通常狀況下，速度場是未知的，且須要經過求解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$動量方程開始，即
$u_C+H_C[u]=B_C^u-D_C^u\left(\frac{\partial p}{\partial x}\right)_C \\ u_F+H_F[u]=B_F^u-D_F^u\left(\frac{\partial p}{\partial x}\right)_F$
$u_f$速度方程與上式相似，只是壓力梯度是與局部鄰近壓力值相關的（ $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{aligned} H_f[u] &= \frac{1}{2}(H_C[u]+H_F[u])=\overline{H_f}[u] \\ B_f^u &= \frac{1}{2}(B_C^u+B_F^u)=\overline{B_f^u} \\ D_f^u &= \frac{1}{2}(D_C^u+D_F^u)=\overline{D_f^u} \end{aligned}$
使用這些插值算得的係數，單元面處的僞動量方程變爲
$u_f+\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{aligned} \overline{H_f}[u] &= \frac{1}{2}(H_C[u]+H_F[u]) \\ &= \frac{1}{2}\left[ -u_C+B_C^u-D_C^u\left(\frac{\partial p}{\partial x}\right)_C -u_F+B_F^u-D_F^u\left(\frac{\partial p}{\partial x}\right)_F \right] \\ &= -\overline{u_f}-\overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f}+\overline{B_f^u} \end{aligned}$
其中用到了係數的近似
$\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{aligned} \overline{D_f^u\left(\frac{\partial p}{\partial x}\right)_f}- \overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f} = & \frac{1}{2}\left[ D_C^u\left(\frac{\partial p}{\partial x}\right)_C+ D_F^u\left(\frac{\partial p}{\partial x}\right)_F \right] \\ &-\frac{1}{2}(D_C^u+D_F^u)\times\frac{1}{2} \left[\left(\frac{\partial p}{\partial x}\right)_C+\left(\frac{\partial p}{\partial x}\right)_F\right] \\ =& \frac{1}{4}D_C^u \left[\left(\frac{\partial p}{\partial x}\right)_C-\left(\frac{\partial p}{\partial x}\right)_F\right] +\frac{1}{4}D_F^u \left[\left(\frac{\partial p}{\partial x}\right)_F-\left(\frac{\partial p}{\partial x}\right)_C\right] \\ \approx & O(\Delta x^2) \end{aligned}$
把 $\overline{H_f}[u]$的表達式代入到單元面 $f$處的僞動量方程中，可求得以Rhie-Chow插值方法表示的單元面處的速度 $u_f$爲
$\begin{aligned} u_f &= -\overline{H_f}[u]+\overline{B_f^u} -\overline{D_f^u}\left(\frac{\partial p}{\partial x}\right)_f \\ &= -\left[ -\overline{u_f}-\overline{D_f^u}\overline{\left(\frac{\partial p}{\partial x}\right)_f}+\overline{B_f^u} \right] +\overline{B_f^u} -\overline{D_f^u}\left(\frac{\partial p}{\partial x}\right)_f \\ &= \underbrace{\overline{u_f}}_{average~velocity}- \underbrace{\overline{D_f^u}\left[ \left(\frac{\partial p}{\partial x}\right)_f-\overline{\left(\frac{\partial p}{\partial x}\right)_f} \right]}_{correction~term} \end{aligned}$

注意，對於一維均勻網格而言，以 $e$界面爲例，上式中
$\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{aligned} v_f &= \overline{v_f}-\overline{D_f^v}\left[ \left(\frac{\partial p}{\partial y}\right)_f-\overline{\left(\frac{\partial p}{\partial y}\right)_f} \right] \\ w_f &= \overline{w_f}-\overline{D_f^w}\left[ \left(\frac{\partial p}{\partial z}\right)_f-\overline{\left(\frac{\partial p}{\partial z}\right)_f} \right] \end{aligned}$
上面仨式子可寫成矢量形式，更適用於多維壓力修正方程的推導，即
$\bold v_f = \overline{\bold v_f}-\overline{\bold D_f^{\bold v}}( \nabla p_f-\overline{\nabla p_f} )$
其中
$\overline{\bold D_f^{\bold v}}= \begin{bmatrix} \overline{D_f^u} & 0 & 0 \\ 0 & \overline{D_f^v} & 0 \\ 0 & 0 & \overline{D_f^w} \end{bmatrix}$
其中 $\nabla p_f$的計算，見第9章第4節，即
$\nabla p_f=\overline{\nabla p_f}+\underbrace{\left[\frac{p_F-p_C}{d_{CF}}-(\overline{\nabla p_f}\cdot\bold e_{CF})\right]\bold e_{CF}}_{Correction~to~interpolated~face~gradient}$
其中
$\overline{\nabla p_f}=g_C\nabla p_C+g_F\nabla p_F$
這樣，可推導出在 $CF$方向上的架構僅由鄰近單元值 $p_F$和 $p_C$給出
$\begin{aligned} \nabla p_f \cdot \bold e_{CF}&=\overline{\nabla p_f} \cdot \bold e_{CF}+ \left[\frac{p_F-p_C}{d_{CF}}-(\overline{\nabla p_f}\cdot\bold e_{CF})\right]\bold e_{CF} \cdot \bold e_{CF}\\ &=\frac{p_F-p_C}{d_{CF}} \end{aligned}$
因爲面速度是與其鄰近連續單元（而非交替單元）的壓力密切相關的，因此棋盤型變量場將會被算法很好地感知，即使在同位網格上也不會出現不符合物理意義的棋盤場了。

5 通常推導

在繼續講解多維同位網格上的壓力修正方程以前，首先講解多維動量方程的離散。

5.1 動量方程離散

將動量方程
$\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$
略做修改，寫爲下述形式
$\underline{\frac{\partial}{\partial t}(\rho \bold v)}+\underline{\nabla\cdot(\rho\bold v\bold v)}= -\nabla p+\underline{\nabla\cdot (\mu \nabla\bold v)}+\nabla\cdot [\mu(\nabla\bold v)^\text{T}] + \bold f_b$
上式的離散形式將在時間區間 $[t-\Delta t/2, t+\Delta t/2]$和單元 $C$上尋找，以下圖所示。

上式中，三個用下劃線標記的項，從左到右分別是非定常項、對流項和擴散項。這些項的離散在前面章節已經講過了。剩下的那些項將顯式估算並做爲源項來對待。對於剪切應力項的第二項的體積分將使用散度定理轉化爲面積分，並寫成面通量加和的形式，即
$\int_{V_C}\nabla\cdot [\mu(\nabla\bold v)^\text{T}] dV=\int_{\partial V_C}[\mu(\nabla\bold v)^\text{T}]\cdot d\bold S= \sum_{f\sim nb(C)}\mu(\nabla\bold v)^\text{T}\cdot \bold S_f$
其中 $(\nabla\bold v)^\text{T}\cdot \bold S_f$在三維座標系統中的展開形式爲
$(\nabla\bold v)^\text{T}\cdot \bold S_f= \begin{bmatrix} \dfrac{\partial u}{\partial x}S_f^x+\dfrac{\partial u}{\partial y}S_f^y+\dfrac{\partial u}{\partial z}S_f^z \\ \dfrac{\partial v}{\partial x}S_f^x+\dfrac{\partial v}{\partial y}S_f^y+\dfrac{\partial v}{\partial z}S_f^z \\ \dfrac{\partial w}{\partial x}S_f^x+\dfrac{\partial w}{\partial y}S_f^y+\dfrac{\partial w}{\partial z}S_f^z \\ \end{bmatrix}$
壓力梯度的體積分也一樣做爲源項處理並顯式估算，即
$\int_{V_C} \nabla p dV=(\nabla p)_C V_C$
或者轉化成面積分來處理，即
$\int_{V_C} \nabla p dV=\int_{\partial V_C}p d\bold S=\sum_{f\sim nb(C)}p_f\bold S_f$
體積力項直接在控制體上積分，得
$\int_{V_C} \bold f_b dV= (\bold f_b)_C V_C$
使用一階Euler格式來離散非定常項，使用高分辨率格式來離散對流項並經過延遲修正方法來將其實現，至於擴散通量，則將其分解爲與網格平齊的隱式部分和顯式的交叉擴散部分，這樣的離散動量方程寫成矢量形式爲
$a_C^{\bold v}\bold v_C+\sum_{F\sim NB(C)}a_F^{\bold v}\bold v_F=\bold b_C^{\bold v}$
其中的係數爲
$\begin{aligned} a_C^{\bold v} &= FluxC_C+\sum_{f\sim nb(C)}FluxC_f \\ a_F^{\bold v} &= FluxF_f \\ \bold b_C^{\bold v} &= -FluxV_C-\sum_{f\sim nb(C)} FluxV_f \\ &+ \underbrace{\sum_{f\sim nb(C)}\mu_f(\nabla\bold v)_f^\text{T}\cdot \bold S_f} _{\footnotesize{\begin{matrix}from~the~2nd~part~of\\~the~shear~stress~term\end{matrix}}} - \underbrace{(\nabla p)_C V_C}_ {\footnotesize{\begin{matrix}from~the\\pressure~gradient\end{matrix}}} \end{aligned}$
其中的面通量用下式計算
$\begin{aligned} FluxC_f &= \underbrace{||\dot m_f, 0||}_{\footnotesize{\begin{matrix}convection \\ contribution\end{matrix}}} + \underbrace{\mu_f \frac{E_f}{d_{CF}}}_{\footnotesize{\begin{matrix}diffusion\\contribution\end{matrix}}} \\ FluxF_f &= -\underbrace{||-\dot m_f, 0||}_{\footnotesize{\begin{matrix}convection \\ contribution\end{matrix}}} - \underbrace{\mu_f \frac{E_f}{d_{CF}}}_{\footnotesize{\begin{matrix}diffusion\\contribution\end{matrix}}} \\ FluxV_f &= \underbrace{-\mu_f(\nabla\bold v)_f\cdot \bold T_f}_{\footnotesize{\begin{matrix}cross~diffusion\\contribution\end{matrix}}} + \underbrace{\dot m_f(\bold v_f^{HR}-\bold v_f^U)}_{\footnotesize{\begin{matrix}high~resolution\\convection~scheme\\contribution\end{matrix}}} \end{aligned}$
單元通量則用下式計算
$\begin{aligned} FluxC_C &= \frac{\rho_C V_C}{\Delta t} \\ FluxV_C &= \underbrace{-\frac{\rho_C^\circ V_C}{\Delta t}\bold v_C^\circ}_ {\footnotesize{\begin{matrix}transient \\ contribution\end{matrix}}}- \underbrace{(\bold f_b)_C V_C}_{\footnotesize{\begin{matrix}source~term\\contribution\end{matrix}}} \end{aligned}$
儘管動量方程的離散代數方程是線性的，然而其係數是依賴於速度和壓力場的。該非線性是要經過迭代過程來處理的，即，在每次迭代開始時，係數是基於前次迭代所得到的因變量的值而計算的。這些係數值的變化會引發 $\bold v$的巨大變化，從而影響收斂速率，甚至會致使發散。爲了減慢變化，可應用欠鬆弛方法，尤爲是所用的瞬態時間步長較大時。將欠鬆弛因子定義爲 $\lambda^{\bold v}$，並採用Patankar隱式欠鬆弛手法，則欠鬆弛動量方程可寫爲
$\frac{a_C^{\bold v}}{\lambda^{\bold v}}\bold v_C+\sum_{F\sim NB(C)}a_F^{\bold v}\bold v_F=\bold b_C^{\bold v}+ \frac{1-\lambda^{\bold v}}{\lambda^{\bold v}}a_C^{\bold v}\bold v_C^{(n)}$
經過從新定義 $a_C^{\bold v}$和 $b_C^{\bold v}$，即
$\begin{aligned} a_C^{\bold v} &\leftarrow \frac{a_C^{\bold v}}{\lambda^{\bold v}} \\ \bold b_C^{\bold v} &\leftarrow \bold b_C^{\bold v}+ \frac{1-\lambda^{\bold v}}{\lambda^{\bold v}}a_C^{\bold v}\bold v_C^{(n)} \end{aligned}$
欠鬆弛動量方程可從新寫爲其原來的樣子
$a_C^{\bold v}\bold v_C+\sum_{F\sim NB(C)}a_F^{\bold v}\bold v_F=\bold b_C^{\bold v}$
爲了推導出同位網格上的壓力修正方程，將壓力梯度從源項 $\bold b_C^{\bold v}$中剝離出來，並顯式計算，即
$\bold b_C^{\bold v}=-V_C(\nabla p)_C+\hat{\bold b}_C^{\bold v}$
因而，動量方程最終變成了下式
$\bold v_C+\sum_{F\sim NB(C)}\frac{a_F^{\bold v}}{a_C^{\bold v}}\bold v_F= -\frac{V_C}{a_C^{\bold v}}(\nabla p)_C+\frac{\hat{\bold b}_C^{\bold v}}{a_C^{\bold v}}$
定義以下矢量算子
$\begin{aligned} \bold H_C[\bold v] &= \sum_{F\sim NB(C)}\frac{a_F^{\bold v}}{a_C^{\bold v}}\bold v_F \\ \bold B_C^{\bold v} &= \frac{\hat{\bold b}_C^{\bold v}}{a_C^{\bold v}} \\ \bold D_C^{\bold v} &= \frac{V_C}{a_C^{\bold v}} \end{aligned}$
動量方程整理爲
$\bold v_C+\bold H_C[\bold v] = -\bold D_C^{\bold v}(\nabla p)_C + \bold B_C^{\bold v}$
該形式將用於後續的推導中。

5.2 同位網格上的壓力修正方程

與交錯網格的時候同樣，最開始用猜想值或者前次迭代值 $(\bold v^{(n)},\dot m^{(n)},p^{(n)})$，求解動量方程（上節最後推出的式子）
$\bold v_C+\bold H_C[\bold v] = -\bold D_C^{\bold v}(\nabla p)_C + \bold B_C^{\bold v}$
來得到知足動量方程的速度場 $\bold v^*$，這樣所得解知足
$\bold v_C^*+\bold H_C[\bold v^*] = -\bold D_C^{\bold v}(\nabla p^{(n)})_C + \bold B_C^{\bold v}$
最終的解是要知足上上式的，那麼這兩個方程（上式和上上式）的差異在於，知足上上式的速度場既知足動量方程又知足連續方程，而知足上式的速度場則只知足動量方程而未必知足連續方程，由於上式作了線性化處理，其計算係數和壓力梯度項所用的壓力和速度是基於前次迭代的值。所以，須要對速度、壓力、質量流量場進行修正，以使它們知足質量守恆。將這些修正量用 $(\bold v',p',\dot m')$表示，則精確解和計算值之間的關係可寫爲
$\begin{aligned} \bold v &= \bold v^* + \bold v' \\ p &= p^{(n)} + p' \\ \dot m &= \dot m^* + \dot m' \end{aligned}$
將質量流量的修正方程代入到連續方程中，得
$\sum_{f\sim nb(C)}\dot m_f' = -\sum_{f\sim nb(C)} \dot m_f^* \quad\quad \dot m_f^* = \rho_f \bold v_f^* \cdot \bold S_f$
速度 $\bold v_f^*$使用Rhie-Chow插值計算，即
$\bold v_f^* = \overline{\bold v_f^*}-\overline{\bold D_f^{\bold v}}( \nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}} )$
當算得的質量流量場知足守恆時， $\sum_{f\sim nb(C)}\dot m_f' = -\sum_{f\sim nb(C)} \dot m_f^*$的右端項爲零，代表無需修正流場。另外一方面，若速度場不正確的話將會產生不平衡的質量流量，並使得該式右端項非零，代表須要修正流場以使其知足質量守恆條件。

質量流量修正項可寫成速度修正項的形式，先用精確解的動量方程 $\bold v_C+\bold H_C[\bold v] = -\bold D_C^{\bold v}(\nabla p)_C + \bold B_C^{\bold v}$減去計算值的動量方程 $\bold v_C^*+\bold H_C[\bold v^*] = -\bold D_C^{\bold v}(\nabla p^{(n)})_C + \bold B_C^{\bold v}$，便可
$\bold v_C'+\bold H_C[\bold v'] = -\bold D_C^{\bold v}(\nabla p')_C$
一樣也可推導出 $F$單元的速度修正值動量方程
$\bold v_F'+\bold H_F[\bold v'] = -\bold D_F^{\bold v}(\nabla p')_F$
單元面上的質量流量修正可寫爲
$\dot m_f'=\rho_f \bold v_f' \cdot \bold S_f$
其中面速度修正是由 $\bold v_f = \overline{\bold v_f}-\overline{\bold D_f^{\bold v}}( \nabla p_f-\overline{\nabla p_f} )$減去 $\bold v_f^*=\overline{\bold v_f^*}-\overline{\bold D_f^{\bold v}}( \nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}})$而獲得的
$\bold v_f' = \overline{\bold v_f'}-\overline{\bold D_f^{\bold v}}( (\nabla p')_f-\overline{\nabla p_f'})$
把上式和上上式代入到連續方程 $\sum_{f\sim nb(C)}\dot m_f' = -\sum_{f\sim nb(C)} \dot m_f^*$中，推導出以下的壓力修正方程
$\begin{aligned} \underline{ \sum_{f\sim nb(C)}\left(\rho_f \overline{\bold v_f'} \cdot \bold S_f \right) +\sum_{f\sim nb(C)}\left(\rho_f \overline{\bold D_f^{\bold v}}~\overline{\nabla p_f'} \cdot \bold S_f\right) } -\sum_{f\sim nb(C)}\left(\rho_f \overline{\bold D_f^{\bold v}} (\nabla p')_f \cdot \bold S_f\right) \\ = -\sum_{f\sim nb(C)} \dot m_f^* \end{aligned}$
在該方程中，下劃線部分表明了鄰近速度對於所考慮單元的速度修正的影響。爲更清晰地理解該影響，可經過將前面的單元速度修正方程式 $\bold v_C'+\bold H_C[\bold v'] = -\bold D_C^{\bold v}(\nabla p')_C$和式 $\bold v_F'+\bold H_F[\bold v'] = -\bold D_F^{\bold v}(\nabla p')_F$二者插值到面上去，推出上式下劃線部分的以下的等效表達式。
$\begin{aligned} & \overline{\bold v_f'}+\overline{\bold H_f}[\bold v'] = -\overline{\bold D_f^{\bold v}}~\overline{({\nabla p'})_f} \\ \Rightarrow & \overline{\bold v_f'}+\overline{\bold D_f^{\bold v}}~\overline{({\nabla p'})_f}=-\overline{\bold H_f}[\bold v'] \end{aligned}$
將上式代入到上上式中，得
$\sum_{f\sim nb(C)}\left(-\rho_f \overline{\bold D_f^{\bold v}} (\nabla p')_f \cdot \bold S_f\right) = -\sum_{f\sim nb(C)} \dot m_f^*+ \underline{\sum_{f\sim nb(C)}\left(\rho_f \overline{\bold H_f}[\bold v'] \cdot \bold S_f\right) }$
或者寫成更加顯式的形式
$\sum_{f\sim nb(C)}\left(-\rho_f \overline{\bold D_f^{\bold v}} (\nabla p')_f \cdot \bold S_f\right) = -\sum_{f\sim nb(C)} \dot m_f^*+ \underline{\sum_{f\sim nb(C)}\left(\rho_f \left(\overline{\sum_{F\sim NB(C)}\frac{a_F^{\bold v}}{a_C^{\bold v}}\bold v'_F} \right) \cdot \bold S_f\right) }\\$
在上式或上上式中，對於下劃線部分的處理相當重要，其決定着方程可解與否。在最初的SIMPLE算法中，是直接把該項忽略掉的，即，認爲在某點處的速度修正是隻跟壓力相關，而跟周圍的速度絕不相干。因爲這只是個修正方程，因此對其中的某些項做以修改甚至是直接扔掉也不會影響到最終的解，由於在收斂狀況下修正會變成零的（能夠這麼來理解，對於某點速度（壓力）的修正，一方面來自於壓力，一方面來自於鄰近速度，那麼若是認爲周圍速度的影響是零（若是收斂的話的確如此），只考慮壓力的影響，至關於只考慮了部分修正，同樣能夠到達最終的收斂狀態（此時壓力和鄰近速度的影響都是零了））。這麼個近似雖然對最終收斂解沒有影響，可是它會影響到收斂的速度哈，由於在每步迭代中，忽略的項越大，那麼帶來的近似偏差就越高了。

上式和上上式中，剩餘項的處理就很簡單了，壓力修正方程的係數，是按照第8章中所講的擴散項的離散方法來得到的，具體來講，是非各向同性（各向異性）擴散的處理方法。

能夠把左端項修改爲梯度點積形式
$\begin{aligned} \left( \overline{\bold D_f^{\bold v}} (\nabla p')_f\right) \cdot \bold S_f &= \left( (\nabla p')_f \overline{\bold D_f^{\bold v}}^T \right) \cdot \bold S_f \\ &= (\nabla p')_f \cdot \left( \overline{\bold D_f^{\bold v}}^T \bold S_f \right)\\ &= (\nabla p')_f \cdot \bold S'_f \end{aligned}$
其中 $\bold S'_f$的展開形式爲
$\bold S'_f=\overline{\bold D_f^{\bold v}}^T \bold S_f= \begin{bmatrix} \overline{D_f^u} & 0 & 0 \\ 0 & \overline{D_f^v} & 0 \\ 0 & 0 & \overline{D_f^w} \end{bmatrix} \begin{bmatrix} S_f^x \\ S_f^y \\ S_f^z \end{bmatrix}= \begin{bmatrix} \overline{D_f^u}S_f^x \\ \overline{D_f^v}S_f^y \\ \overline{D_f^w}S_f^z \end{bmatrix}$
使用 $\bold S'_f$，壓力修正梯度項的離散過程和以前同樣，獲得
$\begin{aligned} (\nabla p')_f \cdot \bold S'_f &= (\nabla p')_f \cdot \bold E_f + (\nabla p')_f \cdot \bold T_f \\ &= \frac{E_f}{d_{CF}}(p'_F-p'_C)+\underline{ (\nabla p')_f\cdot\bold T_f} \end{aligned}$
其中用到了對 $\bold S'_f$以下分解
$\bold S'_f=\bold E_f+\bold T_f$
分解的類型在第8章中已經講過（三種類型：最小修正、正交修正、過鬆弛修正方法）。下劃線項是由網格非正交特性所引發的，能夠選擇忽略或保留。若是忽略，其不會影響最終解，由於這只是修正項。若是保留，其將會在內循環中顯式處理（OpenFOAM中的非正交循環）。因爲每次迭代的求解都是從零壓力修正場開始的，在求解方程時該項須要迭代更新。

忽略掉非正交項的做用，壓力修正方程的線性化形式爲
$\begin{array}{rl}{\displaystyle (\mathrm{\nabla }{p}^{\mathrm{\prime }}{)}_f\cdot {\mathbf{S}}_f^{\mathrm{\prime }}}& {\displaystyle =\frac{E_f}{}}\end{array}$