「その１」で流体の挙動を最低限シミュレートできたはずである。
次は少しステップアップして、流体の可視化と、流体の中に存在する「動く障害物」に対してどう処理するかというところを見ていこう。
【可視化】
これはかなり独学なところもあったが、基本は「渦度」を画面に出力すればらしい映像が得られる、という認識で間違いないと思う。
百聞は一見に如かず
画像左半分が128*128格子の渦度の可視化画面
ちなみに右は圧力の可視化。
渦度を求める方法はいろいろすっ飛ばした結果以下のとおり。
渦度=(A+B+C+D)/4
AとCは、ｘ速度の定義点におけるｙ速度、BとDはｙ速度の定義点におけるｘ速度を表している。
それぞれの速度は、定義点から外れているため、近傍４点から平均を取ってくるなどの方法をとってうまく計算する。
渦度の画面出力は、マイナス値を赤、プラス値を緑などとし、値の大きさと色をフィーリングで比例させればだいたい↑のようなレンダリング結果が得られるだろう。
あと計算高速化のアイディアとして、細かいことをきにしなくていいなら、平均処理をしないで渦度を求めることだってできる。
渦度＝周囲４点の時計（又は反時計） 回りへ の速度を平均したもの
とラフに考えれば、格子の「角」における渦度はどうなるか。
これならｘ速度、ｙ速度の値をそのまま使えるため計算コストが低く抑えられる。
次に、渦度だけでは画面が寂しいというのなら、幾万ものドット粒子を浮遊させるというのもよくある良い可視化手段である。
粒子の移動速度は当然流体の速度と同じで、近傍の速度定義点から求めてやればいい。
ちなみに３次補間と比べ精度の低い線形補間（バイリニア補間）で求めても粒子の動きは見た目問題ないほど非常になめらかであった。
【流体の中で動く障害物】
流体の計算の醍醐味は、障害物の存在下で流体がどのような挙動を取るのかを見て楽しむことにある。
このように、移動する障害物を流体力学の計算にどのように組み込んだらいいのか。
ここでもスタガード格子ということで考えてみる。
まず格子ごとに、この格子が障害物か流体かを識別するための配列変数「kabe」、を作る。
可動障害物でも静止障害物でも、ある時刻ｔにおけるある格子が障害物であるときkabeは0の値を取り、流体であるとき1の値を取る。
次に、流体内に、移動しない障害物があるときの処理をおさらいする。
このとき、障害物の境界で強制的に流速を0にすれば良かったはずだ。
一方、障害物が移動することを考えると。
ここのサイトによると、このときの境界の条件は、流速=障害物の移動速度、となる。
最初に、一番単純化して格子の長さが1m×1mとして、1格子からなる最小単位の障害物がｘ方向に1m/sで動くプログラムを考える。
またΔt=0.2sとする。
最初のt=0.0sの時点では当然こうなる。
濃い灰色と薄い灰色の速度は固定である。そのため移流や圧力の影響を受けない。
また変数kabeは1と0だけで定義できる状態だ。
さて次にt=0.2sではどうなるか。
理想ではこうしたいところではあるが、障害物の側面が速度定義点から外れる上、１格子内に壁成分と流体成分が混じっているため、 変数kabeは1と0だけで定義できない状態となる。
つまり正解はこうなる。
これはt=0.0sの状況と全く同じである。
これでいいのかと思うかもしれないが、これで良い！
そしてt=1.0sまでタイムステップを進めたとき、↓のようにすれば良い。
ということはt=0.8sの時点で一番誤差が大きくなり、1m^2の面積の格子に0.8m^2分の誤差流体が入っていることになる。
こんな言い方をしたくないが、こればかりは仕方がない。
そしてつぎに斜め移動の場合を考える。例えば斜め右下30°に1.0m/sの速度を持つ壁としたら、境界条件はｘ速度=0.87m/s 、ｙ速度=0.5m/s
に固定すれば良い。
そこら辺はさっきのリンク ここ の下の方の説明画像を見ればわかると思う。
ところで話はかわり、タコマ橋という有名な橋があるのだが、これは完成後すぐに横風によるカルマン渦が発生して上下に振動し、ついには崩壊してしまったアメリカの橋だ。
このように、障害物（壁）が流体の力を受けて移動する場合の計算は、ポアソン方程式で圧力計算が終わったあとの障害物周囲の圧力を障害物に働く力として計算すればよい。
これで可動障害物の話は終わり。
最後に。
iPhoneアプリで面白い物を見つけた。
「wind tunnel lite」という格子法で計算されているCFDアプリだ。
指で流体を動かせる機能も付いている。
当然フルにスペックを使って計算されているはずなので電池消耗アプリとしていいかもしれない・・・（そんな馬鹿な）
そして風上差分の大欠点である数値拡散を抑えるため、CIP法の勉強を最近始めたので
次の記事にはCIP法についてわかったことを書いていこうと思う。
今回お世話になった参考文献。
http://news.mynavi.jp/articles/2007/10/09/cedec01/006.html
（やはり西川善司さんの講座は参考になる）
次回：流体力学のプログラムを作りたい！その３
次は少しステップアップして、流体の可視化と、流体の中に存在する「動く障害物」に対してどう処理するかというところを見ていこう。
【可視化】
これはかなり独学なところもあったが、基本は「渦度」を画面に出力すればらしい映像が得られる、という認識で間違いないと思う。
百聞は一見に如かず
画像左半分が128*128格子の渦度の可視化画面
ちなみに右は圧力の可視化。
渦度を求める方法はいろいろすっ飛ばした結果以下のとおり。
渦度=(A+B+C+D)/4
AとCは、ｘ速度の定義点におけるｙ速度、BとDはｙ速度の定義点におけるｘ速度を表している。
それぞれの速度は、定義点から外れているため、近傍４点から平均を取ってくるなどの方法をとってうまく計算する。
渦度の画面出力は、マイナス値を赤、プラス値を緑などとし、値の大きさと色をフィーリングで比例させればだいたい↑のようなレンダリング結果が得られるだろう。
あと計算高速化のアイディアとして、細かいことをきにしなくていいなら、平均処理をしないで渦度を求めることだってできる。
渦度＝周囲４点の時計（又は反時計） 回りへ の速度を平均したもの
とラフに考えれば、格子の「角」における渦度はどうなるか。
これならｘ速度、ｙ速度の値をそのまま使えるため計算コストが低く抑えられる。
次に、渦度だけでは画面が寂しいというのなら、幾万ものドット粒子を浮遊させるというのもよくある良い可視化手段である。
粒子の移動速度は当然流体の速度と同じで、近傍の速度定義点から求めてやればいい。
ちなみに３次補間と比べ精度の低い線形補間（バイリニア補間）で求めても粒子の動きは見た目問題ないほど非常になめらかであった。
【流体の中で動く障害物】
流体の計算の醍醐味は、障害物の存在下で流体がどのような挙動を取るのかを見て楽しむことにある。
このように、移動する障害物を流体力学の計算にどのように組み込んだらいいのか。
ここでもスタガード格子ということで考えてみる。
まず格子ごとに、この格子が障害物か流体かを識別するための配列変数「kabe」、を作る。
可動障害物でも静止障害物でも、ある時刻ｔにおけるある格子が障害物であるときkabeは0の値を取り、流体であるとき1の値を取る。
次に、流体内に、移動しない障害物があるときの処理をおさらいする。
このとき、障害物の境界で強制的に流速を0にすれば良かったはずだ。
一方、障害物が移動することを考えると。
ここのサイトによると、このときの境界の条件は、流速=障害物の移動速度、となる。
最初に、一番単純化して格子の長さが1m×1mとして、1格子からなる最小単位の障害物がｘ方向に1m/sで動くプログラムを考える。
またΔt=0.2sとする。
最初のt=0.0sの時点では当然こうなる。
濃い灰色と薄い灰色の速度は固定である。そのため移流や圧力の影響を受けない。
また変数kabeは1と0だけで定義できる状態だ。
さて次にt=0.2sではどうなるか。
理想ではこうしたいところではあるが、障害物の側面が速度定義点から外れる上、１格子内に壁成分と流体成分が混じっているため、 変数kabeは1と0だけで定義できない状態となる。
つまり正解はこうなる。
これはt=0.0sの状況と全く同じである。
これでいいのかと思うかもしれないが、これで良い！
そしてt=1.0sまでタイムステップを進めたとき、↓のようにすれば良い。
ということはt=0.8sの時点で一番誤差が大きくなり、1m^2の面積の格子に0.8m^2分の誤差流体が入っていることになる。
こんな言い方をしたくないが、こればかりは仕方がない。
そしてつぎに斜め移動の場合を考える。例えば斜め右下30°に1.0m/sの速度を持つ壁としたら、境界条件はｘ速度=0.87m/s 、ｙ速度=0.5m/s
に固定すれば良い。
そこら辺はさっきのリンク ここ の下の方の説明画像を見ればわかると思う。
ところで話はかわり、タコマ橋という有名な橋があるのだが、これは完成後すぐに横風によるカルマン渦が発生して上下に振動し、ついには崩壊してしまったアメリカの橋だ。
このように、障害物（壁）が流体の力を受けて移動する場合の計算は、ポアソン方程式で圧力計算が終わったあとの障害物周囲の圧力を障害物に働く力として計算すればよい。
これで可動障害物の話は終わり。
最後に。
iPhoneアプリで面白い物を見つけた。
「wind tunnel lite」という格子法で計算されているCFDアプリだ。
指で流体を動かせる機能も付いている。
当然フルにスペックを使って計算されているはずなので電池消耗アプリとしていいかもしれない・・・（そんな馬鹿な）
そして風上差分の大欠点である数値拡散を抑えるため、CIP法の勉強を最近始めたので
次の記事にはCIP法についてわかったことを書いていこうと思う。
今回お世話になった参考文献。
http://news.mynavi.jp/articles/2007/10/09/cedec01/006.html
（やはり西川善司さんの講座は参考になる）
次回：流体力学のプログラムを作りたい！その３
