はじめに
前回はMAC法の発展形であるSMAC法について記載しました。今回はこれをさらに発展させたHSMAC法（Highly Simplified Marker And Cell method、別名SOLA法）についてまとめていきます。
SMAC法のおさらい
SMAC法は中間速度を定義することで、ナビエ・ストークス方程式を以下のように二段階に分けた。
\frac{1}{\Delta t}
\left [
\left(
\begin{array}{c}
u^{*} \\
v^{*}\\
\end{array}
\right)
–
\left(
\begin{array}{c}
u^{n} \\
v^{n}\\
\end{array}
\right)
\right]
=
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
p^{n}
+
ν
\left(
\begin{array}{c}
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}\\
\frac{\partial ^2 }{\partial x^2} +\frac{\partial ^2 }{\partial y^2}\\
\end{array}
\right)
\left(
\begin{array}{c}
u^{n} \\
v^{n}\\
\end{array}
\right)
–
\left[
\left(
\begin{array}{c}
u^{n} \\
v^{n}\\
\end{array}
\right)
\cdot
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\right]
\left(
\begin{array}{c}
u^{n} \\
v^{n}\\
\end{array}
\right)\tag{1}
\end{align}
\frac{1}{\Delta t}
\left [
\left(
\begin{array}{c}
u^{n+1} \\
v^{n+1}\\
\end{array}
\right)
–
\left(
\begin{array}{c}
u^{*}\\
v^{*}\\
\end{array}
\right)
\right]
=
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\left(p^{n+1}-p^{n}\right)
\tag{2}
\end{align}
これにより圧力のポアソン方程式をシンプルに変更し収束計算を簡略化することができた。しかしながら、依然として圧力ポアソン方程式を解くため、圧力の境界条件や行列計算が必要である。
またSMAC法（およびMAC法）では圧力のポアソン方程式において、\(D^{n+1}=0\)と置くことで、次ステップで連続の式を満たすような圧力を求めた。
これはこれで正しいが、連続の式の誤差である\(D^{n+1}\)をダイレクトに評価しているわけではないので求めた流速\(u^{n+1},v^{n+1}\)が\(D^{n+1}=0\)をどの程度満たしているのかは分からなかった。仮にポアソン方程式を解く際に十分に反復回数を確保し、収束判定を厳しく設定したとしても、わかるのはポアソン方程式がどれだけ正確に解けたかということだけで、間接的な評価にしかならない。
これらの課題を解決するのが、今回解説するHSMAC法である。
HSMAC法
SMAC法の時と同様に式\((2)\)の両辺で発散を取れば以下のようになる。（詳細は前回記事参照）
\frac{1}{\Delta t}
&\left [
\left(
\frac{\partial u^{n+1}}{\partial x}+\frac{\partial v^{n+1}}{\partial y}
\right)
–
\left(
\frac{\partial u^{*}}{\partial x}+\frac{\partial v^{*}}{\partial y}
\right)
\right]
=-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\left(p^{n+1}-p^{n}\right)
\\\
&→
\frac{D^{n+1}-D^{*}}{\Delta t}
=-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\left(p^{n+1}-p^{n}\right)\\\
&→
D^{n+1}
=D^{*}-\frac{\Delta t}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\left(p^{n+1}-p^{n}\right)
\tag{3}
\end{align}
この式は時間に関してのみ離散化されているが、空間に関しても離散化していくことを考える。空間の離散化は以下に示したようなスタガード格子を用いることとする。
これをもとに式\((3)\)に関して空間の離散化を考えると以下のようになる。
D^{n+1}_{i,j}
=D^{*}_{i,j}-\frac{\Delta t}{ρ}
\left[\frac{p^{n+1}_{i+1,j}-2p^{n+1}_{i,j}+p^{n+1}_{i-1,j}}{(\Delta x)^2}\right]
-\frac{\Delta t}{ρ}
\left[\frac{p^{n+1}_{i,j+1}-2p^{n+1}_{i,j}+p^{n+1}_{i,j-1}}{(\Delta y)^2}\right]
+\frac{\Delta t}{ρ}
\left[\frac{p^{n}_{i+1,j}-2p^{n}_{i,j}+p^{n}_{i-1,j}}{(\Delta x)^2}\right]
+\frac{\Delta t}{ρ}
\left[\frac{p^{n}_{i,j+1}-2p^{n}_{i,j}+p^{n}_{i,j-1}}{(\Delta y)^2}\right]
\tag{4}
\end{align}
ここで圧力の二階微分の離散化には中心差分を用いている。また\(D^{n+1}_{i,j},
D^{*}_{i,j}\)の中身は共通して以下のように書ける。
D_{i,j}^{n+1}
=\frac{u_{i,j}^{n+1}-u_{i-1,j}^{n+1}}{\Delta x}+\frac{v_{i,j}^{n+1}-v_{i-1,j}^{n+1}}{\Delta y}
\tag{5}
\end{align}
D_{i,j}^{*}
=\frac{u_{i,j}^{*}-u_{i-1,j}^{*}}{\Delta x}+\frac{v_{i,j}^{*}-v_{i-1,j}^{*}}{\Delta y}
\tag{6}
\end{align}
ついでに後の議論のため、式\((2)\)に関しても圧力格子と流速格子の位置関係に注意しながら空間離散化しておくと
\left(
\begin{array}{c}
u^{n+1}_{i,j}\\
v^{n+1}_{i,j}\\
\end{array}
\right)
=
\left(
\begin{array}{c}
u^{*}_{i,j}\\
v^{*}_{i,j}\\
\end{array}
\right)
-\frac{\Delta t}{ρ}
\left(
\begin{array}{c}
\frac{p^{n+1}_{i+1,j}-p^{n+1}_{i,j}}{\Delta x} \\
\frac{p^{n+1}_{i,j+1}-p^{n+1}_{i,j} }{\Delta y}\\
\end{array}
\right)
-\frac{\Delta t}{ρ}
\left(
\begin{array}{c}
\frac{p^{n}_{i+1,j}-p^{n}_{i,j}}{\Delta x} \\
\frac{p^{n}_{i,j+1}-p^{n}_{i,j}}{\Delta y}\\
\end{array}
\right)
\tag{7}\end{align}
一つ格子をずらせば、
\left(
\begin{array}{c}
u^{n+1}_{i-1,j}\\
v^{n+1}_{i,j-1}\\
\end{array}
\right)
=
\left(
\begin{array}{c}
u^{*}_{i-1,j}\\
v^{*}_{i,j-1}\\
\end{array}
\right)
-\frac{\Delta t}{ρ}
\left(
\begin{array}{c}
\frac{p^{n+1}_{i,j}-p^{n+1}_{i-1,j}}{\Delta x} \\
\frac{p^{n+1}_{i,j}-p^{n+1}_{i,j-1} }{\Delta y}\\
\end{array}
\right)
-\frac{\Delta t}{ρ}
\left(
\begin{array}{c}
\frac{p^{n}_{i,j}-p^{n}_{i-1,j}}{\Delta x} \\
\frac{p^{n}_{i,j}-p^{n}_{i,j-1}}{\Delta y}\\
\end{array}
\right)
\tag{8}
\end{align}
ここまでで後の議論に必要な式の準備が整った。ここからがHSMAC法のポイントである。SMAC法では式\((4)\)に関して\(D^{n+1}_{i,j}=0\)と置くことで、圧力のポアソン方程式を解きに行った。
しかしながら、HSMAC法では式\((4)\)に関して\(D^{n+1}_{i,j}(p^{n+1}_{i,j})\)とみなし、ニュートン法により\(D^{n+1}_{i,j}(p^{n+1}_{i,j})＝0\)を満たす\(p^{n+1}_{i,j}\)を直接的に見つけに行くことで、圧力のポアソン方程式を解くことを回避する。
ニュートン法とは以下のようなものだった。
\(f(x)=0\)の解は以下の式の反復計算により求まる。ここで\(m\)は反復回数を示す。
^{m+1}x={}^{m}x -\frac{f({}^{m}x)}{f^{\prime}({}^{m}x)}={}^{m}x+\delta x
, {}^{m}\delta x = -\frac{f({}^{m}x)}{f^{\prime}({}^{m}x)}
\end{align}
これを\(D^{n+1}_{i,j}(p^{n+1}_{i,j})=0\)に当てはめると、
{}^{m+1}p^{n+1}_{i,j}={}^{m}p^{n+1}_{i,j} – \frac{D^{n+1}({}^{m}u^{n+1},{}^{m}v^{n+1} )}{\left(\frac{\partial D^{n+1}(p^{n+1}_{i,j})}{\partial p^{n+1}_{i,j}}\right)_{p^{n+1}_{i,j}={}^{m}p^{n+1}_{i,j}}}
\tag{9}
\end{align}
この右辺の分子は式\((5)\)から計算することができるので、流速の関数として表示した。右辺の分母については以下のように考える。式\((5)\)を代入すると
\frac{\partial D^{n+1}(p^{n+1}_{i,j})}{\partial p^{n+1}_{i,j}}
&=
\frac{\partial }{\partial p^{n+1}_{i,j}}\left(
\frac{u^{n+1}_{i,j}-u^{n+1}_{i-1,j}}{\Delta x}+\frac{\ v^{n+1}_{i,j}-v^{n+1}_{i,j-1}}{\Delta y}
\right)\\\
&=
\frac{1}{\Delta x}
\left(
\frac{\partial u^{n+1}_{i,j}}{\partial p^{n+1}_{i,j}}-\frac{\partial u^{n+1}_{i-1,j}}{\partial p^{n+1}_{i,j}}\right)
+\frac{1}{\Delta y}
\left(
\frac{\partial v^{n+1}_{i,j}}{\partial p^{n+1}_{i,j}}-\frac{\partial v^{n+1}_{i,j-1}}{\partial p^{n+1}_{i,j}}\right)\tag{10}
\end{align}
式\((7),(8)\)を代入すると\(p^{n+1}_{i,j}\)に関連する項のみ残るので、
&=
\frac{1}{\Delta x}
\left(
\frac{\Delta t}{\rho \Delta x}+\frac{\Delta t}{\rho \Delta x}\right)
+\frac{1}{\Delta y}
\left(
\frac{\Delta t}{\rho \Delta y}+\frac{\Delta t}{\rho \Delta y}\right)
\\\
&=
\frac{2\Delta t}{\rho}\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)\tag{11}
\end{align}
以上より式\((9)\)は以下のように書ける。
&{}^{m+1}p^{n+1}_{i,j}={}^{m}p^{n+1}_{i,j} – \frac{D_{i,j}^{n+1}({}^{m}u^{n+1},{}^{m}v^{n+1} )}{\frac{2\Delta t}{\rho}\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)}\\\
&→
{}^{m}\delta p_{i,j} = {}^{m+1}p^{n+1}_{i,j}-{}^{m}p^{n+1}_{i,j}=- \frac{D_{i,j}^{n+1}({}^{m}u^{n+1},{}^{m}v^{n+1} )}{\frac{2\Delta t}{\rho}\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)}
\tag{12}
\end{align}
さらに実際には収束を加速するために加速係数\(\omega\)を加えた以下の式が用いられている。
{}^{m}\delta p_{i,j} =-\omega \frac{D_{i,j}^{n+1}({}^{m}u^{n+1},{}^{m}v^{n+1} )}{\frac{2\Delta t}{\rho}\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)}
, 1\leq \omega \leq2 \tag{13}
\end{align}
この反復計算を繰り返すことで、次ステップで\(D^{n+1}_{i,j}=0\)となる圧力の修正量\({}^{m}\delta p_{i,j}\)を求められそうである。
しかしここで注意したいのは、圧力を更新する際にはそれに応じて流速も更新しなければならない点である。具体的には式\((5)\)より右辺分子\(D^{n+1}_{i,j}\)は\(u^{n+1}_{i,j},u^{n+1}_{i-1,j},v^{n+1}_{i,j},v^{n+1}_{i,j-1}\)の4つの流速からなるので、これらも同時に更新する必要がある。これらの更新方法は以下のように考える。
流速\({}^{m+1}u^{n+1}_{i,j}\)を圧力\({}^{m+1}p^{n+1}_{i,j}\)でテイラー展開すると
{}^{m+1}u^{n+1}_{i,j}={}^{m}u^{n+1}_{i,j} + \frac{\partial {}^{m}u^{n+1}_{i,j}}{\partial {}^{m}p^{n+1}_{i,j}}{}^{m}\delta p_{i,j}\tag{14}
\end{align}
式\((7)\)を右辺に用いると
{}^{m+1}u^{n+1}_{i,j}={}^{m}u^{n+1}_{i,j} + \frac{\Delta t}{\rho \Delta x}{}^{m}\delta p_{i,j}
\tag{15}
\end{align}
この式を用いれば圧力の修正量\({}^{m}\delta p\)から流速を更新することができる。式\((7),(8)\)を用いて、\(u^{n+1}_{i-1,j},v^{n+1}_{i,j},v^{n+1}_{i,j-1}\)についても同様に考えれば、
{}^{m+1}u^{n+1}_{i-1,j}={}^{m}u^{n+1}_{i-1,j} – \frac{\Delta t}{\rho \Delta x}{}^{m}\delta p_{i,j}
\tag{16}\\\
{}^{m+1}v^{n+1}_{i,j}={}^{m}v^{n+1}_{i,j} + \frac{\Delta t}{\rho \Delta y}{}^{m}\delta p_{i,j}
\tag{17}\\\
{}^{m+1}v^{n+1}_{i,j-1}={}^{m}v^{n+1}_{i,j-1} – \frac{\Delta t}{\rho \Delta y}{}^{m}\delta p_{i,j}
\tag{18}\\\
\end{align}
以上、式\((15)\sim(18)\)から圧力から流速を更新する方法が分かった。これにより、ついに式\((13)\)の反復計算から圧力の修正量\({}^{m}\delta p\)を求めることができるようになった。
なおこの時の収束判定には式\((5)\)の誤差が一定値以下になることを用いればよいだろう。式\((5)\)を反復回数\(m\)を含んだ形で書けば以下のようになる。これは連続の式をダイレクトに評価していることにほかならない。
D_{i,j}^{n+1}
=\frac{{}^{m}u_{i,j}^{n+1}-{}^{m}u_{i-1,j}^{n+1}}{\Delta x}+\frac{{}^{m}v_{i,j}^{n+1}-{}^{m}v_{i-1,j}^{n+1}}{\Delta y}
\tag{19}
\end{align}
以上がHSMAC法の定式化になるが、ここまでの式変形をより直感的に解釈してみる。式\((13)\)において、\(D_{i,j}^{n+1}>0\)の場合、\(\delta p_{i,j}<0\)となる。さらに式\((15)\sim(18)\)において\(\delta p_{i,j}<0\)の場合、符号を考えるとそのセルの中心に向かって各流速が生じていることが分かる。
これは\((i,j)\)番目のセルで質量の湧き出しがある場合、そのセルの圧力を下げることで質量を吸い込む方向に調整されることを意味している。同様に符号が逆の場合は圧力を上げることで質量が湧き出す方向に調整される。このように、連続の式を満たすように圧力を調整しに行くのが、HSMAC法の本質的な部分なのだろう。
最後にHSMAC法のアルゴリズムをまとめると以下のようになる。
おわりに
今回はHSMAC法についてまとめました。HSMAC法はニュートン法を用いることでポアソン方程式の求解を回避している点が優れていると思います。また今回でMAC系解法は一通り記事にしたことになります。陰解法であるSIMPLE系解法についてもいつか記事にしたいですが実装のほうにトライしたいので、今のところはこのぐらいにしておきます。ではでは。