はじめに
非圧縮Navier-Stokes方程式の数値解析実装に向けて、流体解析の圧力-速度連成解法の1種であるMAC法（Marker And Cell method）について2次元の例で説明します。
基礎方程式
基礎方程式は、下記の2次元非圧縮ナビエ・ストークス方程式
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}=-\frac{1}{ρ} \frac{\partial p}{\partial x} +ν\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)\tag{1}
\end{align}
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}=-\frac{1}{ρ} \frac{\partial p}{\partial y} +ν\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)\tag{2}
\end{align}
及び連続の方程式
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0\tag{3}
\end{align}
の3つとなる。ちなみにナビエ・ストークス方程式\((1),(2)\)に関してベクトル表示して、非定常項のみ左辺に残すと
\frac{\partial }{\partial t}
\left(
\begin{array}{c}
u \\
v
\end{array}
\right)
+
\left[
\left(
\begin{array}{c}
u \\
v
\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 \\
v
\end{array}
\right)
=-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}
\end{array}
\right)
p
+
ν
\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 \\
v
\end{array}
\right) \\\\
→
\frac{\partial }{\partial t}
\left(
\begin{array}{c}
u \\
v
\end{array}
\right)
=
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}
\end{array}
\right)
p
+
ν
\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 \\
v
\end{array}
\right)
–
\left[
\left(
\begin{array}{c}
u \\
v
\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 \\
v
\end{array}
\right)\tag{4}
\end{align}
ここで未知数は\(u,v,p\)の3つであるから、方程式の数から考えても、これらは解くことができるはずである。では具体的にどう求めようか…。これらの式をじっと見てみると、ナビエ・ストークス方程式には流速の時間微分項がある。これを離散化すれば次のタイムステップの流速を求めることができそうだ。空間微分はひとまずほっておいて左辺の時間微分をオイラー前進差分法で離散化してみると、
&\frac{1}{\Delta t}
\left [
\left(
\begin{array}{c}
u^{n+1} \\
v^{n+1}
\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+1}
+
ν
\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)
\\\\
→
&
\left(
\begin{array}{c}
u^{n+1} \\
v^{n+1}
\end{array}
\right)
=
\left(
\begin{array}{c}
u^{n} \\
v^{n}
\end{array}
\right)
+
\Delta t
\left\{
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}
\end{array}
\right)
p^{n+1}
+
ν
\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)
\right\}
\tag{5}
\end{align}
添え字の\(n\)は\(n\)番目のタイムステップを表している。ここで圧力は未来のタイムステップ\(n+1\)で陰的に表現し、流速は現在のタイムステップ\(n\)で陽的に表現していることに注意する。圧力は流速と異なり方程式内に時間微分項はないため、陽的に時間追跡できない（この辺がナビエストークス方程式が難しいポイントの一つなのだろう）。この意味で圧力項はどの解法でも陰解法として解かれる。
また、後述のとおりMAC法は圧力の初期値を不要とする解法である。つまり最初は圧力が未知数でも問題ないのである。仮に現ステップの情報は全て揃っているという視点に立つならば、\(p^{n+1}\)と表現しておいたほうが、未知数感が出てよいだろう。
このように変数毎に採用するタイムステップを変えてしまうことについて、私自身は当初違和感を持った。しかしながら\(\Delta t\)が十分小さければ、これらの操作はなんら問題ない。これについては補足１にて最後に記載しようと思う。
ちなみにここであえて\(p^{n+1}\)を\(p^{n}\)と陽的に表現してMAC法を定式化することも可能である。例えばこちらの記事などは\(p^{n}\)としてMAC法を記載している。内容としては、ポアソン方程式で求まる圧力が\(p^{n+1}\)から\(p^n\)に変わるだけである。個人的にはむしろ流速と圧力を平等に扱ってる分こちらのほうが自然な定式化だと感じる。しかしながら、別記事で説明予定のSMAC法はMAC法とは異なり圧力の初期値（現ステップの圧力値）が必要となる解法である。つまり式内に明示的に\(p^n\)を組み込む必要がある。したがって、初期圧力が不要なMAC法では圧力項は\(p^{n+1}\)と陰的に表現されたほうが後々統一的に表現できてよさそうだ。
式\((5)\)に話を戻す。この式を見る限り次ステップの流速\(u^{n+1},v^{n+1}\)を決めるには、次ステップの圧力\(p^{n+1}\)と現ステップの流速\(u^{n},v^{n}\)が必要である。仮に\(u^{n},v^{n}\)を初期条件として与えたとしても、\(p^{n+1}\)を求める手段（＝圧力を時間発展させる手段）が必要である。そこで提案されているのがSIMPLE法やMAC法といった流体解析における速度-圧力連成解法である。今回はMAC法について記載する。
MAC法
さて、ナビエ・ストークス方程式\((5)\)と連続方程式\((3)\)から実際に\(p^{n+1}\)をもとめられないか考えていく。圧力項が含まれるのはナビエ・ストークス方程式なので、まずナビエ・ストークス方程式に連続の式を組み込めないかと考えるのが自然な発想だろう。具体的には、式\((5)\)の両辺で発散\(div\)を取れば、そのまま連続の方程式を代入できそうである。実際に計算してみる。長くなるので左辺と右辺を分けて考える。
(左辺)&=
\frac{1}{\Delta t}
\left [
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}
\end{array}
\right)
\cdot
\left(
\begin{array}{c}
u^{n+1} \\
v^{n+1}
\end{array}
\right)
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}
\end{array}
\right)
\cdot
\left(
\begin{array}{c}
u^{n} \\
v^{n}
\end{array}
\right)
\right]\\\\
&=
\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^{n}}{\partial x}+\frac{\partial v^{n}}{\partial y}
\right)
\right]\tag{6}
\end{align}
\scriptsize{
(右辺)=
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
p^{n+1}
+
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
ν
\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(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\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)\\\
=-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
\left(
\begin{array}{c}
\frac{\partial }{\partial x}\\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\left(
\begin{array}{c}
u^{n}\\
v^{n}\\
\end{array}
\right)
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x}\\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\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)\\\
=-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
\left(
\frac{\partial u^n}{\partial x} +\frac{\partial v^n}{\partial y}
\right)
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\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{7}
\end{align}
右辺第二項の式変形に注意する。（実際に展開してみるとわかる。） 以上、式\((6),(7)\)より、両辺で発散を取ったナビエ・ストークス方程式は
&\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^{n}}{\partial x}+\frac{\partial v^{n}}{\partial y}
\right)
\right]\\\
&=
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
\left(
\frac{\partial u^n}{\partial x} +\frac{\partial v^n}{\partial y}
\right)
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\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)\\\
&→
\frac{D^{n+1}-D^{n}}{\Delta t}
=
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D^{n}
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\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)\\\
&→
D^{n+1}
=
D^{n}
+\Delta t
\left\{
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D^{n}
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x}\\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\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)
\right\}
\tag{8}
\end{align}
ここで、
D^n=\frac{\partial u^{n}}{\partial x}+\frac{\partial v^{n}}{\partial y}
\tag{9}
\end{align}
式\((8)\)は流速ベクトルの発散\(D^{n}\)（連続の式の左辺部分）の時間発展を記述している。普通に考えれば連続の方程式\((3)\)より、\(D^{n},D^{n+1}\)はゼロになるはずである。連続の方程式は過去から未来に渡って成り立つ流体力学の基本法則なのだから、現実世界では任意のタイムステップで間違いなくこれらはゼロである。（ゼロにしたときに得られる方程式はいわゆる圧力のポアソン方程式と呼ばれている。）
ただし、コンピューター上で行われる数値計算には誤差がある。この為、\(D^{n},D^{n+1}\)を厳密にゼロとすることはできない。詳しくは補足2にて述べるが\(D^{n},D^{n+1}=0\)として計算を進めていくと、次第に連続の式をまったく満たさない意味のない解が生じることとなる。
ではどうするか。このままでは\(p^{n+1}\)を求めることはできなそうだ。そこでMAC法では次のようなテクニックを適用する。
まず、\(D^{n}\)は誤差により完全にゼロではないこと、また現ステップの流速から計算可能な量であることを踏まえ、残しておくこととする。そのうえで、あえて\(D^{n+1}=0\)と置くこととする。すると式\((8)\)は以下のようになる。
D^{n}
+\Delta t
\left\{
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D^{n}
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\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)
\right\}
=0\tag{10}
\end{align}
この式は「例え現ステップで連続の方程式を満たさない\(（D^{n}\neq 0）\)としても、それを考慮したうえで、次ステップでは連続の方程式を満たす\(（D^{n+1}=0）\)ように現ステップの流速\(u^{n},v^{n}\)から次ステップの圧力\(p^{n+1}\)を決めにいく。」ということを意味している。（なぜなら式\((10)\)を式\((8)\)に代入すれば、\(D^{n+1}=0\)になるのだから。）
圧力\(p^{n+1}\)のみが左辺にある形に変形すると
&
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D^{n}
–
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\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)
=-\frac{D^{n}}{\Delta t}\\\
&→
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
=-\frac{D^{n}}{\Delta t}
–
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D^{n}
+
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\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)\\\
&→
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
=
\frac{ρD^{n}}{\Delta t}
+
ρν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D^{n}
–
ρ\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\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{11}
\end{align}
右辺第一項はきわめて小さい\(\Delta t\)で割り算しているため、多くの場合、同様に\(D^{n}\)を有する右辺第二項よりも圧倒的に大きいとされるらしい。従って、右辺第二項は無視してしまうこととする。
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p^{n+1}
=
\frac{ρD^{n}}{\Delta t}
-ρ\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y} \\
\end{array}
\right)
\cdot
\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{12}
\end{align}
現ステップの流速\(u^{n},v^{n}\)が既知だとすると右辺は計算可能なので、この方程式は左辺の圧力\(p^{n+1}\)に関するポアソン方程式になっていることがわかる。つまり、現ステップの流速\(u^{n},v^{n}\)を初期値としてを与えてやれば、次のステップで連続の式を満たすような\(（D^{n+1}=0）\)、次ステップの圧力\(p^{n+1}\)を求めることができる。
そして次ステップの圧力\(p^{n+1}\)が求まれば、式\((5)\)を用いて、次ステップの流速\(u^{n+1},v^{n+1}\)を求めることができる。この時\(u^{n+1},v^{n+1}\)は、連続の式を満たしているはずである。なぜなら、そうなるように式\((12)\)で次ステップの圧力\(p^{n+1}\)を調整しているのだから。（もちろん厳密には\(D^{n+1}=0\)にならないが、極力満たすように維持することができる。）
\(u^{n+1},v^{n+1}\)が求まったならば、再度、式\((12)\)を用いて、\(p^{n+2}\)を求めることができる。このように式\((5),(12)\)を利用することで、ついに流速・圧力の時間発展を計算することができるようになった。
この手法の特筆すべきは点は、式\((8)\)に対して、「\(D^{n}\)は残しておき、\(D^{n+1}\)はゼロとする」とする点だと思う。これにより現タイムステップの連続の方程式の誤差を許容しつつも、次ステップで連続の方程式を満たすように制約を課すことに成功している。
最後にMAC法のアルゴリズムをもう一度まとめることとする。
| （１） | 流速\(u^{n},v^{n}\)を用いて、式\((12)\)より圧力値\(p^{n+1}\)を求める。 |
| （２） | 圧力\(p^{n+1}\)を用いて、式\((5)\)より\(n+1\)ステップにおける流速\(u^{n+1},v^{n+1}\)を求める。 |
| （３） | （１）に戻る。 |
おわりに
今回はMAC法について記載しました。連続の式をあえて満たさないとするところが、すごい発想だと思います。次はSMAC法について記載していきます。
補足１：項ごとに陰的・陽的な取り扱いを変えていいのか？
式\((5)\)を見ると時間に関して離散化する際に、流速は現在のタイムステップ\(n\)で陽的に、圧力は一つ先のタイムステップ\(n+1\)で陰的に表現している。このような操作は許されるのだろうか？運動方程式とは現在に対して成り立つのだから、変数ごとに採用するタイムステップを変えてしまうのはおかしい気がする…と私自身、当初違和感を持っていたが、以下のような考え方で一応納得している。
まず以下のように外力が複数働く運動方程式を考えてみる。
m\frac{dv(t)}{dt} = f_{1}(t)+f_{2}(t)\tag{13}
\end{align}
これを時間に関して離散化することを考えるが、この時\(f_{1}(t)\)は陰的に、\(f_{2}(t)\)は陽的に扱い以下にように変形してみる。
m\frac{v(t+\Delta t)-v(t)}{\Delta t} = f_{1}(t+\Delta t)+f_{2}(t)\\\\
→v(t+\Delta t)＝v(t)+\frac{\Delta t}{m}\left(f_{1}(t+\Delta t)+f_{2}(t)\right)\tag{14}
\end{align}
ここで\(f_{1}(t)\)は式\((5)\)における圧力を、\(f_{2}(t)\)は流速を想定しており、正にこういった離散化が可能なのかが問題である。
ここで\( f_{1}(t+\Delta t)\)に関して変数\((t+\Delta t)\)でテイラー展開してみると
f_{1}(t+\Delta t)＝f_{1}(0)+f_{1}^{\prime}(0)(t+\Delta t)+\frac{1}{2}f_{1}^{\prime\prime}(0)(t+\Delta t)^2+\cdots\tag{15}
\end{align}
式\((15)\)の二次までを式\((14)\)に代入すれば
&v(t+\Delta t)＝v(t)+\frac{\Delta t}{m}\left(f_{1}(0)+f_{1}^{\prime}(0)(t+\Delta t)+\frac{1}{2}f_{1}^{\prime\prime}(0)(t+\Delta t)^2\right)+\frac{\Delta t}{m}f_{2}(t)\\\\
&→
v(t+\Delta t)＝v(t)+\frac{\Delta t}{m}f_{1}(0)+\frac{\Delta t}{m}f_{1}^{\prime}(0)(t+\Delta t)+\frac{\Delta t}{2m}f_{1}^{\prime\prime}(0)(t+\Delta t)^2+\frac{\Delta t}{m}f_{2}(t)\\\\
&→
v(t+\Delta t)＝v(t)+\frac{\Delta t}{m}f_{1}(0)+\frac{\Delta t}{m}f_{1}^{\prime}(0)t+\frac{\Delta t}{m}f_{1}^{\prime}(0)\Delta t+\frac{\Delta t}{2m}f_{1}^{\prime\prime}(0)(t^2+2t\Delta t+t(\Delta t)^2)+\frac{\Delta t}{m}f_{2}(t)
\tag{16}
\end{align}
\((\Delta t)^2\)以上の高次の項は無視すると
v(t+\Delta t)＝v(t)+\frac{\Delta t}{m}f_{1}(0)+\frac{\Delta t}{m}f_{1}^{\prime}(0)t+\frac{\Delta t}{2m}f_{1}^{\prime\prime}(0)t^2+\frac{\Delta t}{m}f_{2}(t)
\tag{17}
\end{align}
右辺２～４項は\(f_1(t)\)を2次までテイラー展開した形に等しいので、これらを\(f_1(t)\)と置けば
v(t+\Delta t)＝v(t)+\frac{\Delta t}{m}f_{1}(t)+\frac{\Delta t}{m}f_{2}(t)\\\\
→
m\frac{v(t+\Delta t)-v(t)}{\Delta t}＝f_{1}(t)+f_{2}(t)\\\\
\tag{18}
\end{align}
これは式\((13)\)において、\(f_1(t),f_2(t)\)ともに陽的に扱った場合に等しい。つまり当初\(f_1(t)\)を陰的に扱ったにもかかわらず、\(\Delta t\)が十分小さいとしてテイラー展開を施せば、陽的に取り扱ったことと等しくなる。
以上より、\(\Delta t\)が十分小さい限り、各項を陰的に取り扱おうが、陽的に取り扱おうが、数式上は何の問題もないものと思われる。
補足２：\(D=0\)として扱ってはいけない理由
式\((8)\)は離散化しない式として書けば以下のように書ける。
\frac{\partial D}{\partial t}
=
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p
+
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D-
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\left[
\left(
\begin{array}{c}
u\\
v\\
\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\\
v\\
\end{array}
\right)
\tag{19}
\end{align}
ここで連続の式が成り立つとして\(D=0\)とすると、
&0
=
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p-
\left(
\begin{array}{c}
\frac{\partial }{\partial x}\\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\left[
\left(
\begin{array}{c}
u\\
v\\
\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\\
v\\
\end{array}
\right)\\\
→
&\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
p
=-
\left(
\begin{array}{c}
\frac{\partial }{\partial x}\\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\cdot
\left[
\left(
\begin{array}{c}
u\\
v\\
\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\\
v\\
\end{array}
\right)
\tag{20}
\end{align}
これはいわゆる圧力のポアソン方程式と呼ばれる式であり、流速を与えることで圧力を決めることができる。では仮にこの式を起点として、初期の流速から圧力を決めたとしよう。その流速と圧力を再度式\((19)\)に代入すると、\(D\)に関する以下の関係式が得られる。
\frac{\partial D}{\partial t}
=
ν
\left(
\frac{\partial ^2 }{\partial x^2} + \frac{\partial ^2 }{\partial y^2}
\right)
D
\tag{21}
\end{align}
これは連続の式\(D=0\)とは異なる。つまり、もともと式\((19)\)に連続の式を課すことで、圧力のポアソン方程式\((20)\)が得られたが、圧力のポアソン方程式が成り立ったからといって連続の式\(D=0\)が成り立つとは限らないのである。
もちろん初期条件として境界も含めて完全に\(D=0\)を満たす流れ場を与えることができれば、式\((21)\)の右辺は０になるから、次のステップでも連続の式は満たされることになる。しかしながら、そのような条件を与えるのは一般的に困難である。逆を言えば\(D\neq0\)の初期条件を与えてしまうと、仮にその値が微小だったとしても式\((21)\)の右辺は値を持つことになるので、時間とともに\(D\)がどんどん増幅していくこととなる。この為、何らかの方法で、\(D\sim0\)を保ち続けるギミックが必要であり、この1つが上述のMAC法によるアプローチとなる。
