はじめに
今回は前回の記事の続編として、流体解析の圧力-速度連成解法の一種であるSMAC法(Simplified Marker And Cell method)について解説します。
MAC法のおさらい
MAC法は連続の式の誤差を考慮した圧力のポアソン方程式
\begin{align}
\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{1}
\end{align}
\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{1}
\end{align}
を用いて、現ステップの流速\(u^n,v^n\)から次ステップの圧力\(p^{n+1}\)を求め、さらに離散化されたナビエ・ストークス方程式
\begin{align}
\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)\tag{2}
\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^{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)\tag{2}
\end{align}
を用いて、現ステップの流速\(u^n,v^n\)および次ステップの圧力\(p^{n+1}\)から次ステップの流速\(u^{n+1},v^{n+1}\)を求めることができた。
しかし圧力のポアソン方程式\((1)\)は解くのに行列演算が必要である。また、右辺に移流項を含んでいるため煩雑な形をしている。SMAC法はこの部分を少し簡略化する手法である。
SMAC法
式\((2)\)を以下のように二つに分割することを考える。
\begin{align}
\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{3}
\end{align}
\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{3}
\end{align}
\begin{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)
=
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\delta p
\tag{4}
\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)
=
-\frac{1}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\delta p
\tag{4}
\end{align}
ここで\(u^{*},v^{*}\)は移流項を切り離すために定義される中間速度である。これらの式を足し合わせると式\((2)\)になっていることが分かる。またMAC法にはなかった現在ステップの圧力\(p^{n}\)が新たに出現していることにも注意する。つまりSMAC法は圧力の初期値が必要な計算手法なのである。
ここで速度補正量\( \delta u,\delta v\)を以下のように定義する。
\begin{align}
\left(
\begin{array}{c}
\delta u\\
\delta v\\
\end{array}
\right)
=
\left(
\begin{array}{c}
u^{n+1} \\
v^{n+1}\\
\end{array}
\right)
–
\left(
\begin{array}{c}
u^{*}\\
v^{*}\\
\end{array}
\right)
\tag{5}
\end{align}
\left(
\begin{array}{c}
\delta u\\
\delta v\\
\end{array}
\right)
=
\left(
\begin{array}{c}
u^{n+1} \\
v^{n+1}\\
\end{array}
\right)
–
\left(
\begin{array}{c}
u^{*}\\
v^{*}\\
\end{array}
\right)
\tag{5}
\end{align}
これ用いると式\((4)\)は以下のようにも書ける。
\begin{align}
\left(
\begin{array}{c}
\delta u\\
\delta v\\
\end{array}
\right)
=
-\frac{\Delta t}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\delta p
\tag{6}
\end{align}
\left(
\begin{array}{c}
\delta u\\
\delta v\\
\end{array}
\right)
=
-\frac{\Delta t}{ρ}
\left(
\begin{array}{c}
\frac{\partial }{\partial x} \\
\frac{\partial }{\partial y}\\
\end{array}
\right)
\delta p
\tag{6}
\end{align}
つまり式\((4)\)は圧力の補正量\(\delta p\)から流速の補正量\( \delta u,\delta v\)が生じることを意味している。
MAC法の時と同様に式\((4)\)の両辺で発散を取ると
\begin{align}
（左辺）
&=
\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^{*}\\
v^{*}\\
\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^{*}}{\partial x}+\frac{\partial v^{*}}{\partial y}
\right)
\right]
\tag{7}
\end{align}
（左辺）
&=
\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^{*}\\
v^{*}\\
\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^{*}}{\partial x}+\frac{\partial v^{*}}{\partial y}
\right)
\right]
\tag{7}
\end{align}
\begin{align}
（右辺）=
-\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)
\delta p
=
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\delta p
\tag{8}
\end{align}
よって
（右辺）=
-\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)
\delta p
=
-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\delta p
\tag{8}
\end{align}
\begin{align}
\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)
\delta p
\\\
&→
\frac{D^{n+1}-D^{*}}{\Delta t}
=-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\delta p
\tag{9}
\end{align}
ここで\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)
\delta p
\\\
&→
\frac{D^{n+1}-D^{*}}{\Delta t}
=-\frac{1}{ρ}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\delta p
\tag{9}
\end{align}
\begin{align}
D^{*}
=\frac{\partial u^{*}}{\partial x}+\frac{\partial v^{*}}{\partial y}
\tag{10}
\end{align}
D^{*}
=\frac{\partial u^{*}}{\partial x}+\frac{\partial v^{*}}{\partial y}
\tag{10}
\end{align}
と置いている。
MAC法の時と同じ要領で\(D^{*}\)は残し、\(D^{n+1}=0\)とすると式\((9)\)は
\begin{align}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\delta p
=
\frac{ρD^{*}}{\Delta t}
\tag{11}
\end{align}
\left(
\frac{\partial ^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\right)
\delta p
=
\frac{ρD^{*}}{\Delta t}
\tag{11}
\end{align}
この式はMAC法における式\((1)\)に相当する。式\((1)\)と比較すると、右辺の移流項が消えており、よりシンプルな形となっている。これを解くことで\(\delta p\)を求めることができる。
以上よりSMAC法のアルゴリズムは下記のとおりである。
SMAC法はMAC法に比べて、ポアソン方程式式\((11)\)の移流項が省かれ単純になっているのが利点である。しかしながら、行列を解かなければいけないという点はMAC法と同様である。
おわりに
今回はSMAC法について解説しました。中間速度の導入によりポアソン方程式がシンプルになり、収束計算が簡略化される点が優れていると思います。次はHSMAC法について記載していきます。