大学院生のための数値流体力学入門④：偏微分方程式の解法

Navier-Stokes方程式の両辺に\(\nabla\)を内積する

\begin{align}
&\frac{\partial \mathbf{u}}{\partial t}+\nabla\cdot\left(\mathbf{uu}\right)
=-\frac{\nabla p}{\rho}
+\nu_{_{E}}\Delta\mathbf{u}
+\mathbf{f} \\
\\ \rightarrow \quad
&\frac{\partial \left(\nabla\cdot\mathbf{u}\right)}{\partial t}
+\nabla\cdot\nabla\cdot\left(\mathbf{uu}\right)
=-\frac{\nabla\cdot\nabla p}{\rho}+\nu_{_{E}}\Delta\left(\nabla\cdot\mathbf{u}\right)
+\nabla\cdot\mathbf{f} \\
\end{align}

連続の式 \(\nabla\cdot\mathbf{u}=0\) より

\begin{align}
&\nabla\cdot\nabla\cdot\left(\mathbf{uu}\right)
=-\frac{\nabla^{2} p}{\rho}
+\nabla\cdot\mathbf{f} \\
\\ \rightarrow \quad
&\frac{\nabla^{2} p}{\rho}
=-\nabla\cdot\nabla\cdot\left(\mathbf{uu}\right)
+\nabla\cdot\mathbf{f} \\
\end{align}

簡単のため，次の図のような1次元流れの有限体積法による離散化を考える

\(P\)が注目セル，\(W\)と\(E\)は隣接セル（West，East），\(w\)と\(e\)は注目セルと隣接セルの境界面である

一般に，物理量\(\phi\)の輸送方程式は時間項，対流項，粘性項，ソース項（生成項＋散逸項）から構成され，次のようにあらわされる

\begin{equation}
\frac{\partial}{\partial t}(\rho\phi)
+\frac{\partial}{\partial x}(\rho u\phi)
-\frac{\partial}{\partial x}\left(\Gamma\frac{\partial \phi}{\partial x}\right)
=
S
\tag{i}
\end{equation}

ここで，

• \(\phi=u\)，\(\Gamma=\mu+\mu_{t}\)，\(S=-\frac{\partial p}{\partial x}\)とすればNavier-Stokes方程式
• \(\phi=k\)，\(\Gamma=\mu+\frac{\mu_{t}}{\sigma_{k}}\)，\(S=P_{k}-\varepsilon\)とすれば標準k-εモデルの\(k\)の輸送方程式，
• \(\phi=\varepsilon\)，\(\Gamma=\mu+\frac{\mu_{t}}{\sigma_{\varepsilon}}\)，\(S=\left(C_{\varepsilon_{1}}P_{k}-C_{\varepsilon_{2}}\varepsilon\right)\frac{\varepsilon}{k}\)とすれば標準k-εモデルの\(\varepsilon\)の輸送方程式

である

有限体積法は保存則をセルの境界からの流入/流出およびセル内部での発生/消滅によって表すので，(i)をセル\(P\)について積分する

\begin{equation}
\int{\frac{\partial}{\partial t}(\rho\phi)}dx
+\int{\frac{\partial}{\partial x}(\rho u\phi)}dx
-\int{\frac{\partial}{\partial x}\left(\Gamma\frac{\partial \phi}{\partial x}\right)}dx
=
\int{S}dx
\tag{ii}
\end{equation}

(ii)の各項について離散化を行う

時間項

陰解法である後退オイラー法を用いて次のように離散化する

\begin{align}
\int{\frac{\partial}{\partial t}(\rho\phi)}dx
&\simeq \frac{(\rho\phi)_{P}-(\rho\phi)_{P}^{(n-1)}}{\Delta t}\Delta x \\
&=\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}-\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}^{(n-1)}
\end{align}

ただし，上付き添え字\(^{(n-1)}\)は1つ前の時間ステップを表している

対流項

2次精度中心差分と2次精度中心補間を用いて次のように離散化する

\begin{align}
\int{\frac{\partial}{\partial x}(\rho u\phi)}dx &\simeq \frac{(\rho u \phi)_{e}-(\rho u \phi)_{w}}{\Delta x}\Delta x =(\rho u)_{e} \phi_{e}-(\rho u)_{w} \phi_{w} \\
&=c_{e}\phi_{e}-c_{w}\phi_{w} =c_{e}\frac{\phi_{_{E}}+\phi_{{_P}}}{2}-c_{w}\frac{\phi_{{_P}}+\phi_{_{W}}}{2} \\ \\
\rightarrow \quad \int{\frac{\partial}{\partial x}(\rho u\phi)}dx &=\left(\frac{1}{2}c_{e}-\frac{1}{2}c_{w}\right)\phi_{_{P}}-\frac{1}{2}c_{w}\phi_{_{W}}+\frac{1}{2}c_{e}\phi_{_{E}}
\end{align}

\begin{align}
c_{e}=(\rho u)_{e}=\rho\frac{u_{_{E}}+u_{_{P}}}{2}, \quad
c_{w}=(\rho u)_{w}=\rho\frac{u_{_{P}}+u_{_{W}}}{2}
\end{align}

拡散項

2次精度中心差分を用いて次のように離散化する

\begin{align}
\int{\frac{\partial}{\partial x}\left(\Gamma\frac{\partial \phi}{\partial x}\right)}dx
&\simeq \frac{\left(\Gamma\frac{\partial \phi}{\partial x}\right)_{e}-\left(\Gamma\frac{\partial \phi}{\partial x}\right)_{w}}{\Delta x}\Delta x
=\left(\Gamma\frac{\partial \phi}{\partial x}\right)_{e}-\left(\Gamma\frac{\partial \phi}{\partial x}\right)_{w} \\
&=\frac{\Gamma\left(\phi_{_{E}}-\phi_{_{P}}\right)}{\Delta x}-\frac{\Gamma\left(\phi_{_{P}}-\phi_{_{W}}\right)}{\Delta x} \\
&=\left(\frac{\Gamma}{\Delta x}\right)_{e}\left(\phi_{_{E}}-\phi_{_{P}}\right)
-\left(\frac{\Gamma}{\Delta x}\right)_{w}\left(\phi_{_{P}}-\phi_{_{W}}\right)\\ \\
&=d_{e}\left(\phi_{_{E}}-\phi_{_{P}}\right)-d_{w}\left(\phi_{_{P}}-\phi_{_{W}}\right)\\ \\

\rightarrow \quad
\int{\frac{\partial}{\partial x}\left(\Gamma\frac{\partial \phi}{\partial x}\right)}dx
&=\left(-d_{e}-d_{w}\right)\phi_{_{P}}+d_{w}\phi_{_{W}}+d_{e}\phi_{_{E}}
\end{align}

\begin{align}
d_{e}=\left(\frac{\Gamma}{\Delta x}\right)_{e}, \quad
d_{w}=\left(\frac{\Gamma}{\Delta x}\right)_{w}
\end{align}

ソース項

\begin{align}
\int{S}dx
&\simeq S\Delta x
\end{align}

すべての項を(ii)に代入すると

\begin{align}
&\int{\frac{\partial}{\partial t}(\rho\phi)}dx
+\int{\frac{\partial}{\partial x}(\rho u\phi)}dx
-\int{\frac{\partial}{\partial x}\left(\Gamma\frac{\partial \phi}{\partial x}\right)}dx
=
\int{S}dx \\ \\
\rightarrow \quad
&\left\{\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}-\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}^{(n-1)}\right\}
+\left\{\left(\frac{1}{2}c_{e}-\frac{1}{2}c_{w}\right)\phi_{_{P}}-\frac{1}{2}c_{w}\phi_{_{W}}+\frac{1}{2}c_{e}\phi_{_{E}}\right\} \\ \\
&\quad -\left\{\left(-d_{e}-d_{w}\right)\phi_{_{P}}+d_{w}\phi_{_{W}}+d_{e}\phi_{_{E}}\right\}
=
S\Delta x \\ \\
\rightarrow \quad
& \left(\frac{\rho\Delta x}{\Delta t}+\frac{1}{2}c_{e}-\frac{1}{2}c_{w}+d_{e}+d_{w}\right)\phi_{_{P}}
+\left(-\frac{1}{2}c_{w}-d_{w}\right)\phi_{_{W}}
+\left(\frac{1}{2}c_{e}-d_{e}\right)\phi_{_{E}} \\ \\
&\quad =
S\Delta x
+\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}^{(n-1)}
\end{align}

まとめると次のようになる

\begin{align}
&a_{_{P}}\phi_{_{P}}+a_{_{W}}\phi_{_{W}}+a_{_{E}}\phi_{_{E}} = b \\
\\ \rightarrow \quad
&a_{_{P}}\phi_{_{P}}+\sum{a_{nb}\phi_{nb}} = b
\end{align}

ただし，それぞれの係数は以下の通り

\begin{align}
&a_{_{P}}=-a_{_{W}}-a_{_{E}}+c_{e}-c_{w}+\frac{\rho\Delta x}{\Delta t} \qquad&
&c_{e}=(\rho u)_{e}=\rho\frac{u_{_{E}}+u_{_{P}}}{2} \\
&a_{_{W}}=-\frac{1}{2}c_{w}-d_{w} &
&c_{w}=(\rho u)_{w}=\rho\frac{u_{_{P}}+u_{_{W}}}{2}\\
&a_{_{E}}=\frac{1}{2}c_{e}-d_{e} &
&d_{e}=\left(\frac{\Gamma}{\Delta x}\right)_{e} \\
&b=S\Delta x+\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}^{(n-1)} &
&d_{w}=\left(\frac{\Gamma}{\Delta x}\right)_{w}\\
\end{align}

最初に述べたように，\(\phi=u\)，\(\Gamma=\mu+\mu_{t}\)，\(S=-\frac{\partial p}{\partial x}\)，\(s_{_{P}}=\frac{\rho\Delta x}{\Delta t}\phi_{_{P}}^{(n-1)}\)とすれば，離散化されたNavier-Stokes方程式が得られる

\begin{equation}
a_{_{P}}u_{_{P}}+\sum{a_{nb}u_{nb}}
=-\frac{\partial p}{\partial x}+s_{_{P}}
\end{equation}

今回はオイラー後退法と2次精度中心差分を用いたが，それ以外の離散化スキームを用いても同様の式を導出できる

離散化されたNavier-Stokes方程式と連続の式に \(\mathbf{u}_{_{P}}=\mathbf{u}^{*}_{_{P}}+\mathbf{u'}_{_{P}}\) と \(p=p^{*}+p'\)を代入することで，SIMPLE法に必要な式を導出する

（連続の式）

\begin{equation}
\nabla\cdot\mathbf{u}_{_{P}}=0
\end{equation}

（離散化されたNavier-Stokes方程式）

\begin{equation}
\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}}
=-\nabla p+\mathbf{s}_{_{P}}
\end{equation}

基本式

まず，離散化されたNavier-Stokes方程式に \(\mathbf{u}_{_{P}}=\mathbf{u}^{*}_{_{P}}+\mathbf{u'}_{_{P}}\) と \(p=p^{*}+p'\)を代入すると次のようになる

\begin{align}
&\mathbf{A}_{_{P}}\left(\mathbf{u}^{*}_{_{P}}+\mathbf{u'}_{_{P}}\right)
+\sum{\mathbf{A}_{nb}\left(\mathbf{u}^{*}_{nb}+\mathbf{u'}_{nb}\right)}
=-\nabla \left(p^{*}+p'\right)+\mathbf{s}_{_{P}} \\
\\ \rightarrow \quad
&\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{*}+\mathbf{s}_{_{P}}
-\left(\mathbf{A}_{_{P}}\mathbf{u'}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}
+\nabla p'
\right)
\tag{i}
\end{align}

十分に時間進行を行うと，(i)の右辺第3項は0になるので

\begin{align}
&\mathbf{A}_{_{P}}\mathbf{u'}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}
+\nabla p' =0 \\
\\ \rightarrow \quad
&\mathbf{u'}_{_{P}}
=
-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'
\tag{ii}
\end{align}

よって，(i)は次のようになる

\begin{equation}
\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{*}+\mathbf{s}_{_{P}}
\tag{iii}
\end{equation}

次に，連続の式に \(\mathbf{u}_{_{P}}=\mathbf{u}^{*}_{_{P}}+\mathbf{u'}_{_{P}}\) を代入する

\begin{align}
\nabla\cdot\mathbf{u}^{*}_{_{P}}+\nabla\cdot\mathbf{u'}_{_{P}}=0
\end{align}

(ii)を代入すると

\begin{align}
&\nabla\cdot\mathbf{u}^{*}_{_{P}}+\nabla\cdot\left(-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)=0 \\
\\ \rightarrow \quad
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}}
-\nabla\cdot\left(\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}\right)
\tag{iv}
\end{align}

(ii)と(iv)より，速度と圧力の修正量に関する連立方程式が得られた

\begin{align}
&\mathbf{u'}_{_{P}}
=-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'
-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}} \\ \\
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}}
-\nabla\cdot\left(\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}\right)
\end{align}

上式は\(\mathbf{u'}\)と\(p'\)についての完全陰解法であり，\(\mathbf{u'}\)と\(p'\)を求めるには大規模な連立一次方程式を解く必要があるため，多大な計算コストがかかってしまう

そこで，「\(\mathbf{u'}\)と\(p'\)は解が収束したらどうせ0になるのだから，それまでの過程で多少不正確な値であっても解の信頼性には関係ない」という考え方に基づき，\(\mathbf{u'}_{nb}\)に関する項を無視することにする

\begin{align}
&\mathbf{u'}_{_{P}}
=-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'
\tag{5.4}\\ \\
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}}
\tag{5.3}
\end{align}

さらに，(iii)における\(p^{*}\)については1つ前のステップにおける圧力\(p^{(n-1)}\)を用いることにすれば

\begin{align}
&p^{(n)}=p^{(n-1)}+p'
\tag{5.5'} \\ \\
&\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{(n-1)}+\mathbf{s}_{_{P}}
\tag{5.2'}
\end{align}

これにより，\(\mathbf{u}^{*}_{_{P}}\) → \(p'\) → \(\mathbf{u'}_{_{P}}\) → \(p^{(n)}(=p^{*})\) → \(\mathbf{u}^{(n)}_{_{P}}\) という順番で求めれば，すべての変数を計算することができるようになった

本来は完全陰解法で求めるはずの\(\mathbf{u'}_{_{P}}\)と\(p'\)を，\(\mathbf{u'}_{nb}\)に関する項を無視することで陽的に求めるという手法が，SIMPLE法が Semi-Implict Method（半陰解法）と呼ばれる所以である

不足緩和

SIMPLE法では，速度と圧力の修正量として\(\mathbf{u'}_{nb}\)に関する項を無視した不正確な値を使用していることが原因で，時間進行の序盤（修正量が大きいとき）に計算が発散しやすくなってしまう

上記の計算の発散を抑制するため，速度と圧力の更新の際に次のような緩和係数を導入する（これを不足緩和という）

\begin{align}
&\mathbf{u}^{*}_{_{P}}
=\mathbf{u}^{n-1}_{_{P}}+\alpha_{u}\left(\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}-\nabla p^{(n-1)}}{\mathbf{A}_{_{P}}}-\mathbf{u}^{(n-1)}_{_{P}}\right)
\tag{v} \\ \\
&p^{(n)}=p^{(n-1)}+\alpha_{p}p'
\tag{5.5} \\ \\
\end{align}

圧力に対する不足緩和は修正量\(p'\)を小さくするように陽的に行われているが，速度に対する不足緩和は推定値\(\mathbf{u}^{*}_{_{P}}\)の変化を小さくする形で陰的に導入されている

(v)を整理すると，(5.2)が得られる

\begin{align}
&\mathbf{u}^{*}_{_{P}}
=\mathbf{u}^{n-1}_{_{P}}+\alpha_{u}\left(\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}-\nabla p^{(n-1)}}{\mathbf{A}_{_{P}}}-\mathbf{u}^{(n-1)}_{_{P}}\right) \\
\\ \rightarrow \quad
&\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
=\mathbf{A}_{_{P}}\mathbf{u}^{n-1}_{_{P}}
+\alpha_{u}\left(\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}-\nabla p^{(n-1)}-\mathbf{A}_{_{P}}\mathbf{u}^{(n-1)}_{_{P}}\right) \\
\\ \rightarrow \quad
&\frac{1}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
=\frac{1}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}^{n-1}_{_{P}}
+\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}-\nabla p^{(n-1)}
-\mathbf{A}_{_{P}}\mathbf{u}^{(n-1)}_{_{P}} \\
\\ \rightarrow \quad
&\frac{1}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{(n-1)}
+\mathbf{s}_{_{P}}
+\frac{1-\alpha_{u}}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}^{(n-1)}_{_{P}}
\tag{5.2}
\end{align}

参考
≫非圧縮性流体計算の圧力-速度連成手法 (PDF)

先述したSIMPLE法の導出より，次の関係がわかっている

\begin{align}
&\mathbf{u}^{(n)}=\mathbf{u}^{*}+\mathbf{u'} \tag{i} \\
&p^{(n)}=p^{(n-1)}+p' \tag{ii} \\ \\
&\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{(n-1)}+\mathbf{s}_{_{P}}
\tag{iii} \\ \\
&\mathbf{u'}_{_{P}}
=-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'
\tag{iv} \\
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}}
\tag{v}
\end{align}

これらの式を用いて，修正量（\(\mathbf{u'}\)，\(p'\)）を用いないSIMPLE法を導出する

基本形

(iii)を次のように変形する

\begin{align}
&\mathbf{u}^{*}_{_{P}}
=\mathbf{\tilde{u}}^{*}_{_{P}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n-1)}
\tag{vi} \\ \\
&\mathbf{\tilde{u}}^{*}_{_{P}}
=\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}}{\mathbf{A}_{_{P}}}
\tag{6.2} \\ \\
\end{align}

(v)に(vi)を代入し，\(p^{(n)}=p^{(n-1)}+p'\) を適用する

\begin{align}
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\left(\mathbf{\tilde{u}}^{*}_{_{P}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n-1)}\right)\\
\\ \rightarrow \quad
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n-1)}+\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\mathbf{\tilde{u}}^{*}_{_{P}} \\
\\ \rightarrow \quad
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}\right)
=\nabla\cdot\mathbf{\tilde{u}}^{*}_{_{P}}
\tag{6.3} \\
\end{align}

同様に，\(\mathbf{u}^{(n)}=\mathbf{u}^{*}+\mathbf{u'}\) に(vi)と(iv)を代入し，\(p^{(n)}=p^{(n-1)}+p'\) を適用する

\begin{align}
&\mathbf{u}^{(n)}=\mathbf{\tilde{u}}^{*}_{_{P}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n-1)}
+\left(-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)\\
\\ \rightarrow \quad
&\mathbf{u}^{(n)}=\mathbf{\tilde{u}}^{*}_{_{P}}
-\frac{1}{\mathbf{A}_{_{P}}}\left(\nabla p^{(n-1)}
+\nabla p'\right)\\
\\ \rightarrow \quad
&\mathbf{u}^{(n)}=\mathbf{\tilde{u}}^{*}_{_{P}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}
\tag{6.4} \\
\end{align}

これにより，\(\mathbf{u}^{*}_{_{P}}\) → \(\mathbf{\tilde{u}}^{*}_{_{P}}\) → \(p^{(n)}\) → \(\mathbf{u}^{(n)}_{_{P}}\)という順番で求めれば，すべての変数を計算することができるようになった

不足緩和

修正量を用いないSIMPLE法では，速度と圧力の不足緩和を次のように行う

\begin{align}
&\mathbf{u}^{*}_{_{P}}
=\mathbf{u}^{n-1}_{_{P}}+\alpha_{u}\left(\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}-\nabla p^{(n-1)}}{\mathbf{A}_{_{P}}}-\mathbf{u}^{(n-1)}_{_{P}}\right)
\tag{v} \\ \\
&p^{(n)}=p^{(n-1)}+\alpha_{p}\left(p^{(n)}-p^{(n-1)}\right)
\tag{5.5} \\ \\
\end{align}

圧力に対する不足緩和の方法が若干変更されている

参考
≫非圧縮性流体計算の圧力-速度連成手法 (PDF)

ある時刻\(t_{n}\)における速度 \(\mathbf{u}^{(n)}\) と圧力 \(p^{(n)}\) の精度を高めるために，速度と圧力の修正を1次修正量（\(\mathbf{u'}\)，\(p'\)）と2次修正量（\(\mathbf{u^{\prime\prime}}\)，\(p^{\prime\prime}\)）の2回にわけて行い，1次推定値と1次修正量の和を2次推定値（\(\mathbf{u}^{**}\)，\(p^{**}\)）として定義する

\begin{align}
&\mathbf{u}^{(n)}=\mathbf{u}^{*}+\left(\mathbf{u'}+\mathbf{u^{\prime\prime}}\right), \qquad &
&p^{(n)}=p^{*}+\left(p'+p^{\prime\prime}\right) \\
&\mathbf{u}^{**}=\mathbf{u}^{*}+\mathbf{u'}, \qquad &
&p^{**}=p^{*}+p' \tag{i}\\
\\ \rightarrow \quad
&\mathbf{u}^{(n)}=\mathbf{u}^{**}+\mathbf{u^{\prime\prime}}, \qquad &
&p^{(n)}=p^{**}+p^{\prime\prime} \tag{ii}\\
\end{align}

これら2つの修正量を用いて各時刻の速度と圧力を計算する手法がPISO法である

基本形

まず，離散化されたNavier-Stokes方程式に \(\mathbf{u}_{_{P}}=\mathbf{u}^{*}+\mathbf{u'}+\mathbf{u^{\prime\prime}}\) と \(p=p^{*}+p'+p^{\prime\prime}\)を代入すると次のようになる

\begin{align}
\mathbf{A}_{_{P}}\left(\mathbf{u}^{*}_{_{P}}+\mathbf{u'}_{_{P}}+\mathbf{u^{\prime\prime}}_{_{P}}\right)
+\sum{\mathbf{A}_{nb}\left(\mathbf{u}^{*}_{nb}+\mathbf{u'}_{nb}+\mathbf{u^{\prime\prime}}_{nb}\right)}
=-\nabla \left(p^{*}+p'+p^{\prime\prime}\right)+\mathbf{s}_{_{P}}
\end{align}

これを次の2つの段階に分けて考える

\begin{align}
&\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{*}+\mathbf{s}_{_{P}}
-\left(\mathbf{A}_{_{P}}\mathbf{u'}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}
+\nabla p'
\right) \\ \\
&\mathbf{A}_{_{P}}\mathbf{u}^{**}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{**}_{nb}}
=-\nabla p^{**}+\mathbf{s}_{_{P}}
-\left(\mathbf{A}_{_{P}}\mathbf{u^{\prime\prime}}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u^{\prime\prime}}_{nb}}
+\nabla p^{\prime\prime}
\right) \\ \\
\end{align}

各時刻において真の解と推定値が一致するためには，上式の右辺第3項は0にならなければならないので

\begin{align}
&\mathbf{u'}_{_{P}}
=-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p' \tag{A} \\ \\
&\mathbf{u^{\prime\prime}}_{_{P}}
= -\frac{\sum{\mathbf{A}_{nb}\mathbf{u^{\prime\prime}}_{nb}}}{\mathbf{A}_{_{P}}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{\prime\prime} \\ \\
\end{align}

\begin{align}
&\mathbf{A}_{_{P}}\mathbf{u}^{*}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{*}_{nb}}
=-\nabla p^{*}+\mathbf{s}_{_{P}} \tag{iii} \\ \\
&\mathbf{A}_{_{P}}\mathbf{u}^{**}_{_{P}}
+\sum{\mathbf{A}_{nb}\mathbf{u}^{**}_{nb}}
=-\nabla p^{**}+\mathbf{s}_{_{P}}
\end{align}

次に，連続の式に \(\mathbf{u}_{_{P}}=\mathbf{u}^{*}_{_{P}}+\mathbf{u'}_{_{P}}\) と \(\mathbf{u}_{_{P}}=\mathbf{u}^{**}_{_{P}}+\mathbf{u^{\prime\prime}}_{_{P}}\) をそれぞれ代入し，SIMPLE法と同様に\(\mathbf{u}_{nb}\)に関する項を無視すると次のようになる

\begin{align}
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}} \tag{iv} \\
&\mathbf{u'}_{_{P}}
=-\frac{1}{\mathbf{A}_{_{P}}}\nabla p' \tag{v} \\ \\
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{\prime\prime}\right)
=\nabla\cdot\mathbf{u}^{**}_{_{P}} \tag{vi} \\
&\mathbf{u^{\prime\prime}}_{_{P}}
=-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{\prime\prime} \tag{vii} \\ \\
\end{align}

上記の考えに基づいて，PISO法では次に示す①～⑨の手順を繰り返すことによって非定常な速度場と圧力場を計算している

2回の修正は何を意味しているのか

速度の1次推定値\(\mathbf{u}_{_{P}}^{*}\)は離散化された偏微分方程式を解くことで求めているが，速度の2次推定値は1次推定値と1次修正量の和 \(\mathbf{u}_{_{P}}^{**}=\mathbf{u}_{_{P}}^{*}+\mathbf{u'}_{_{P}}\) で求めている

手順⑥の(vi)に \(\mathbf{u}_{_{P}}^{**}=\mathbf{u}_{_{P}}^{*}+\mathbf{u'}_{_{P}}\)を代入すると

\begin{align}
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{\prime\prime}\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}}+\nabla\cdot\mathbf{u'}_{_{P}}
\end{align}

(A)より\(\mathbf{u'}_{_{P}}\)を代入すると次のようになる

\begin{align}
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{\prime\prime}\right)
=\nabla\cdot\mathbf{u}^{*}_{_{P}}
+\nabla\cdot\left(-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'\right)\\
\\ \rightarrow \quad
&\nabla\cdot\left\{\frac{1}{\mathbf{A}_{_{P}}}\nabla \left(p'+p^{\prime\prime}\right)\right\}
=\nabla\cdot\mathbf{u}^{*}_{_{P}}
-\nabla\cdot\left(\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}\right) \\
\end{align}

同様に，\(\mathbf{u'}+\mathbf{u^{\prime\prime}}\)に(A)と(vii)を代入すると次のようになる

\begin{align}
&\mathbf{u'}+\mathbf{u^{\prime\prime}}
=-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}}
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p'
-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{\prime\prime} \\
\\ \rightarrow \quad
&\mathbf{u'}+\mathbf{u^{\prime\prime}}
=-\frac{1}{\mathbf{A}_{_{P}}}\left(\nabla p'+\nabla p^{\prime\prime}\right)
-\frac{\sum{\mathbf{A}_{nb}\mathbf{u'}_{nb}}}{\mathbf{A}_{_{P}}} \\
\end{align}

上式から，PISO法では速度と圧力の修正を行う2回行うことで，速度と圧力の合計修正量に速度の1次修正量\(\mathbf{u'}_{nb}\)の項を考慮していることが分かる

依然として速度の2次修正量による\(\mathbf{u^{\prime\prime}}_{nb}\)の項は無視しているが，速度の修正量による項を完全に無視していたSIMPLE法と比べると，時々刻々における速度と圧力の精度は向上している

PISO法の簡略化

基本形の手順をよく見てみると，③～⑤と⑥～⑧は変数の見た目が違うだけで，全く同じ計算をしていることがわかる

全く同じ計算をしているならプログラムにおいてメモリを使いまわすことができるし，メモリを使いまわすことができるなら修正回数を2回に限定する必要はない

よって，上記に示したPISO法の計算手順は次のように簡略化できる

修正量を用いないPISO法

PISO法もSIMPLE法と同様に修正量を用いない形に書き換えることができ，修正量を用いない簡略化したPISO法が1番最初に示した手順である

不足緩和

PISO法では，速度に対する不足緩和はSIMPLE法と同様に行うが，圧力に対する不足緩和は必要ないとされている

参考
≫非圧縮性流体計算の圧力-速度連成手法 (PDF)