はじめに
以前の記事
『最速で流体シミュレーションを実装するチュートリアル』で
MAC法について触れていましたが，あまりにも説明が適当過ぎたので補足として解説していきます．
まずは上記の記事を一度ご一読することをお勧めします．
前回，圧力場を更新するために以下のような式を用いるという話をしました．
\[
\nabla^2 p = - \left\{ \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} \right\}
\]
上記の式は，ナビエ・ストークス方程式を無次元化，外力をゼロとした場合の式です．無次元化をしないで，外力も考慮すると以下のような式になります．
\[
\frac{1}{\rho} \nabla^2 p = - \left\{ \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} \right\} + \nabla \cdot \vec{g}
\]
今回はこの式がどのような過程で導出されているのかについて記していきたいと思います．
ナビエ・ストークス方程式の発散を計算する
はじめに，ナビエ・ストークス方程式はこんな感じの式でした．
\[
\frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla) \vec{v} = - \frac{1}{\rho} \nabla p + \nu \nabla^2 \vec{v} + \vec{g}
\]
そして流体が非圧縮という条件を加えます．流速\( \vec{v} \)の発散が0なので，ある点から急に流体が湧き出したり，または吸い込まれて消えることはないということを示している式でした．
\[
\nabla \cdot \vec{v} = 0
\]
ここで，これら2つの式をうまく組み合わせていきたいのです．連立方程式を代入法で解くような感覚で式を1本にまとめていければ楽なのですが，そのままでは何かを代入しても綺麗にまとまらないことがわかると思います．
\( \nabla \cdot \vec{v} = 0 \)という条件があるので，ナビエ・ストークス方程式の両辺に対して発散（\( \mathrm{div} □ = \nabla \cdot □ \)）を加えると何となくうまくいきそうな気がします．
ということで，ナビエ・ストークス方程式の発散を計算してみます．元の式は全ての項がベクトルでしたが，以下のように発散を加えると全ての項がスカラーになります．
\[
\nabla \cdot \frac{\partial \vec{v}}{\partial t} + \nabla \cdot \{ (\vec{v} \cdot \nabla) \vec{v} \} = - \frac{1}{\rho} \nabla \cdot (\nabla p) + \nu \nabla \cdot ( \nabla^2 \vec{v} ) + \nabla \cdot \vec{g}
\]
左辺の第1項に着目してみます．もし，この項を以下のように変形出来たら式の中に\( \nabla \cdot \vec{v} \)を作ることができます．実際，\( \nabla \)を展開して計算すると以下のような変形が可能であることがわかると思います．
\( \nabla \)を展開するときは必ず，①括弧の中から計算する，②左から計算するということに注意してください．ベクトル解析の計算なので非可換性には十分気を付ける必要があります．
\[
\nabla \cdot \frac{\partial \vec{v}}{\partial t} = \frac{\partial ( \nabla \cdot \vec{v} )}{\partial t}
\]
次に右辺の第1項を見ていきます．一度，\( \nabla \)を展開するとラプラシアン\( \nabla^2 \)で括れることがわかると思います．
\[
- \frac{1}{\rho} \nabla \cdot (\nabla p) = - \frac{1}{\rho} \nabla^2 p
\]
続いて，右辺の第2項を見ていきます．これについても以下のようにちょっとややこしい変形をすることで\( \nabla \cdot \vec{v} \)を作り出せます．
\[
\begin{align}
& \nu \nabla \cdot ( \nabla^2 \vec{v} ) \\
&= \nu \nabla \cdot (\nabla \cdot \nabla \vec{v}) \\
&= \nu \left( \frac{\partial^3 v_x}{\partial x^3} + \frac{\partial^3 v_x}{\partial x \partial y^2} + \frac{\partial^3 v_y}{\partial x^2 \partial y} + \frac{\partial^3 v_y}{\partial y^3} \right) \\
&= \nu \left\{ \frac{\partial^2}{\partial x^2} \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) + \frac{\partial^2}{\partial y^2} \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) \right\} \\
&= \nu \left\{ \frac{\partial^2}{\partial x^2} \left( \nabla \cdot \vec{v} \right) + \frac{\partial^2}{\partial y^2} \left( \nabla \cdot \vec{v} \right) \right\} \\
&= \nu \nabla^2 ( \nabla \cdot \vec{v} )
\end{align}
\]
一旦、ここまでの式を整理します．
\[
\frac{\partial ( \nabla \cdot \vec{v} )}{\partial t} + \nabla \cdot \{ (\vec{v} \cdot \nabla) \vec{v} \} = - \frac{1}{\rho} \nabla^2 p + \nu \nabla^2 ( \nabla \cdot \vec{v} ) + \nabla \cdot \vec{g}
\]
\( \nabla \cdot \vec{v} = 0 \)なので，
\[
\nabla \cdot \{ (\vec{v} \cdot \nabla) \vec{v} \} = - \frac{1}{\rho} \nabla^2 p + \nabla \cdot \vec{g}
\]
移項してさらに整理します．
\[
\frac{1}{\rho} \nabla^2 p = - \nabla \cdot \{ (\vec{v} \cdot \nabla) \vec{v} \} + \nabla \cdot \vec{g}
\]
かなり形になってきました．あとは右辺第1項のごちゃごちゃしたところを何とかすれば良さそうです．
∇・{ (v・∇)v } を計算する
ここからも基本的な流れは同じです．ごちゃごちゃしていますが\( \nabla \)を展開して整理していくだけです．計算ルールは同じです
①括弧の中から計算する，②左から計算するです．
ほとんど説明することが無いので途中式をひたすら記しておきます．
\[
\begin{align}
& \nabla \cdot \{ (\vec{v} \cdot \nabla) \vec{v} \} \\
&= \nabla \cdot \left( v_x \frac{\partial \vec{v}}{\partial x} + v_y \frac{\partial \vec{v}}{\partial y} \right) \\
&= \frac{\partial}{\partial x} \left( v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} \right) + \frac{\partial}{\partial y} \left( v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} \right)
\end{align}
\]
この先を計算したいのですが，積で連なっている関数を微分しようとしているので，
積の微分公式を用いて計算していきます．\( (fg)^{\prime} = f^{\prime}g + fg^{\prime} \)ってやつです．ちゃんと覚えていましたか？
\[
\begin{align}
&= \left( \frac{\partial v_x}{\partial x} \right)^2 + v_x \frac{\partial^2 v_x}{\partial x^2} + \frac{\partial v_y}{\partial x} \frac{\partial v_x}{\partial y} + v_y \frac{\partial^2 v_x}{\partial x \partial y} \\
&+ \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} + v_x \frac{\partial^2 v_y}{\partial x \partial y} + \left( \frac{\partial v_y}{\partial y} \right)^2 + v_y \frac{\partial^2 v_y}{\partial y^2} \\
&= \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} \\
&+ v_x \frac{\partial}{\partial x} \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) + v_y \frac{\partial}{\partial y} \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) \\
&= \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} \\
&+ v_x \frac{\partial}{\partial x} \left( \nabla \cdot \vec{v} \right) + v_y \frac{\partial}{\partial y} \left( \nabla \cdot \vec{v} \right)
\end{align}
\]
\( \nabla \cdot \vec{v} = 0 \)なので，
\[
\begin{align}
&= \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x}
\end{align}
\]
\( \nabla \cdot \{ (\vec{v} \cdot \nabla) \vec{v} \} \)の計算ができたので，前項で求めた式に代入すると一番最初に示した式の形になりました（やったぜ）．
\[
\frac{1}{\rho} \nabla^2 p = - \left\{ \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} \right\} + \nabla \cdot \vec{g}
\]
この式，\( \nabla \)で括り切れてないんですよね．今回は2次元で計算してみましたが，空間の次元が増えると式の形も変わってきます．
3次元空間の場合を考えてみる
折角なので3次元も考えてみたいと思います．これについても2次元と同じ考え方で導出できます．しかし，ここに途中式を書くと大変なことになるので最終的な結果だけ載せておきます．（数式詰め込み過ぎるとページがどんどん重くなっていくんですよ）
\[
\begin{align}
\small{
\frac{1}{\rho} \nabla^2 p = - \left\{ \left( \frac{\partial v_x}{\partial x} \right)^2 + \left( \frac{\partial v_y}{\partial y} \right)^2 + \left( \frac{\partial v_z}{\partial z} \right)^2
+ 2 \frac{\partial v_x}{\partial y} \frac{\partial v_y}{\partial x} + 2 \frac{\partial v_y}{\partial z} \frac{\partial v_z}{\partial y} + 2 \frac{\partial v_z}{\partial x} \frac{\partial v_x}{\partial z} \right\}
+ \nabla \cdot \vec{g}
}
\end{align}
\]
最後に
これ合ってるんですかねえ，有識者の意見を聞いてみたい．（素直にそういう専門書を買って勉強するのがセオリーなんでしょうが，本が苦手なのです）
MAC法は何回も反復して計算させるので計算量が多くなってしまうのが難点ですよね．これを解消するために高速な
SMAC法なんて方法があるらしいですよ，なんでも圧力場を更新する計算が陽的にできるそうなのです．
つまり\( p^{(n + 1)} = ○○○ \)といった感じに1回の計算で圧力場が更新できるという方法らしいので，時間と心に余裕が出来たらその辺も勉強してみたいですねえ．
更新履歴
2024/6/7：粘性項の消し方について，嘘の説明が含まれていたので内容を修正しました．