対象の読者
対象は個人ゲーム製作やCG関係で、流体の精度よりまずは動くプログラムが欲しい、学びたいという方。なので今回は一番単純で簡単な動く流体力学プログラムをつくろうと言うのが目標だ。最終的に
こんな感じのプログラムが動くところまでやる。
導入
まず2Dでも3Dでも流体力学で避けて通れないのはナビエストークスの式（NS式）([1]wikipedia より)
(画像はwikipediaより)
左辺 - 第1項 : 時間微分項、第2項 : 移流項（対流項）
右辺 - 第1項 : 圧力項、第2項 : 粘性項（拡散項）、第3項 : 外力項
ちなみにこれでも簡単なほうらしい。これだけだと未知数が1つおおいままなので、∇・v = 0 という連続の式を連立して解けということだ（ｖは速度ベクトル場）
(画像はwikipediaより)
三角形の逆みたいな記号があるけどなんだこれは？内積の掛け算で使う「・」があるけど、ベクトルの計算式とかあるのか？と、私は最初見た時いろいろわからないとこだらけで、結局理解するのに２ヶ月くらいかかった。
実際のとこ∇をよく理解せずプログラムを組んで動いてしまったので、この記事を書いていた当時は「∇とかわからなくても四則演算でプログラムは作れる」「今はとりあえず動くプログラムを作ることが最優先」とか思っていた。(やっぱそれじゃダメだろと思ったので補足記事書きました)
粒子法か格子法か
今回の記事では格子法についてのみ解説を行う。流体計算は大きく粒子法と格子法の2つがあり、前者は流体を粒子として計算、後者は流体の入った空間を何百という格子で分割して各々の格子の圧力や速度を計算するというもの。
どうも格子法の方が簡単だということがわかった。(粒子法専門の方には大変失礼ですが・・・)
ちなみに粒子法もたまにゲームやCGに組み込まれているのを見る。Position Based Fluidsというものらしい。
格子法で流体力学の完成形
例
[13]http://www.youtube.com/watch?v=4Ov6tguRr8E
格子法では、2Dの場合、各々の格子の中の圧力、x速度、y速度を時間ステップごとに更新していくものだ。
格子の種類
さて格子法でやると決めた後はどの種類の格子にするかを決める必要がある。大きく分けて通常格子、コロケート格子とスタガード格子がある。
通常格子
コロケート格子
スタガード格子
スタガード格子の方が少し複雑で難しい。が、プログラムを作っている途中でわかったが最終的にこっちの方が楽になる。
というわけでスタガード格子でやることに決定！
これは東工大のサイトが参考になった。
[2.1]http://www.eto.titech.ac.jp/contents/sub03/chapter08-1.html
で、次にスタガード格子の速度と圧力をNS式でどのように離散化して時間ステップを進めるかということを考える。
全体像と変数
流体格子の分割は10*10格子とする。
壁を含めると全計算領域は12*12格子である。
必要な変数はまずは圧力変数 「p」 、これは12*12格子
次にｘ速度変数 「vx」 、これは13*12
逆にｙ速度変数 「vy」 は12*13
また後々、ダイバージェンス(∇・v=0)の計算で使う各格子での「発散」値を格納するデータ 「s」 も必要。これは圧力と同じく12*12
「p」「vx」「vy」の一時的な記憶変数として「p_after」「vx_after」「vy_after」もそれぞれ同じ配列要素で用意。
そして、以降話を分り易くするためにΔxやΔyは1.0として考える （少し衝撃的な決断だが）なので以降の式にΔxΔyが省略されているがそこはご了承を。
まともにプログラムが動いてからΔxやΔyをいじってみればいい。またΔtは1.0にすると計算が発散するので、今は0.1～0.2くらいだと思ってもらえればいい。
フローチャート
フローチャートはこうなる
上から順にみていこう
①移流粘性項
最初に移流と粘性を分けて考える。
移流項
移流とは、流れる前の速度変数 vx vy をつかって一定時間流れた後の速度を求める、ということだ。(と思ってましたがこれはあくまで陽解法の説明です。陰解法では...略)流れ後の速度変数をそれぞれvx_after vy_afterとすると
（i,jは配列の要素、整数）
u=vx[i][j]
v=( vy[i-1][j] + vy[i][j] + vy[i-1][j+1] + vy[i][j+1] )/4
if (u>=0.0)&(v>=0.0){
vx_after[i][j] = vx[i][j] - u*(vx[i][j] - vx[i-1][j])*Δt - v*(vx[i][j]-vx[i][j-1])*Δt
}
if (u<0 .0)&(v>=0.0){
vx_after[i][j] = vx[i][j] - u*(vx[i+1][j] - vx[i][j])*Δt - v*(vx[i][j]-vx[i][j-1])*Δt
}
if (u>=0.0)&(v<0 .0){
vx_after[i][j] = vx[i][j] - u*(vx[i][j] - vx[i-1][j])*Δt - v*(vx[i][j+1]-vx[i][j])*Δt
}
if (u<0.0)&(v<0.0){
vx_after[i][j] = vx[i][j] - u*(vx[i+1][j] - vx[i][j])*Δt - v*(vx[i][j+1]-vx[i][j])*Δt
}
ここでuとvは、代入が行なわれる変数vx_after[i][j]の位置にあるｘ速度、ｙ速度を意味する。よってｕ=vx[i][j]でいいがｖは、上のスタガード格子の性質上、ｘ速度の定義点とｙ速度の定義点がずれているため
v=( vy[i-1][j] + vy[i][j] + vy[i-1][j+1] + vy[i][j+1] )/4
と、４つの近くの格子点の情報を足して４で割る必要がある。
ここまででvxの移流の話が終わり、次はvyの移流である。
上の計算でvxとvyを置き換えるだけでいいがｕとｖは
u=( vx[i][j-1] + vx[i+1][j-1] + vx[i+1][j] + vx[i][j] )/4
v=vy[i][j]
とする。
したがって
u=( vx[i][j-1] + vx[i+1][j-1] + vx[i+1][j] + vx[i][j] )/4
v=vy[i][j]
if (u>=0.0)&(v>=0.0){
vy_after[i][j] = vy[i][j] - u*(vy[i][j] - vy[i-1][j])*Δt - v*(vy[i][j]-vy[i][j-1])*Δt
}
if (u<0 .0)&(v>=0.0){
vy_after[i][j] = vy[i][j] - u*(vy[i+1][j] - vy[i][j])*Δt - v*(vy[i][j]-vy[i][j-1])*Δt
}
if (u>=0.0)&(v<0 .0){
vy_after[i][j] = vy[i][j] - u*(vy[i][j] - vy[i-1][j])*Δt - v*(vy[i][j+1]-vy[i][j])*Δt
}
if (u<0.0)&(v<0.0){
vy_after[i][j] = vy[i][j] - u*(vy[i+1][j] - vy[i][j])*Δt - v*(vy[i][j+1]-vy[i][j])*Δt
}
これを計算領域内すべてのi jで計算する。
全計算が終わったならvx vy にvx_after vy_afterを代入する。計算領域はイメージ図１で言う緑の部分と接するか重なるところだけ。
つまり最外層の１層だけは参照されるのみで代入は行なわれない。
またまた東工大のサイトが参考になった。
[2.2]http://www.eto.titech.ac.jp/contents/sub03/chapter03.html
移流のイメージ
移流のイメージとしては例えば以下の画像の黄緑の部分の速度が（1.0,1.0）で白の部分が（0.0,0.0）、話を分り易くするためにΔt = 1.0だとしたらこうなる。
しかし実際はΔtや速度の値が 0.148 とか0.009とか不規則な数であり、これにより移流計算するごとに速度の情報が崩れまくる。
特に一次風上差分法はとてもその形が崩れやすい。黄緑色の領域がじわじわ薄くなって広がってしまうイメージだ。
数値流体力学ではこのプロファイルの崩壊をいかに防いで移流計算するかということに、そのシミュレーション精度がかかっている。
精度を上げるためには、Cubic セミ・ラグランジアン法やCIP法を用いると良いようだ。
今回は一番簡単な一次風上差分を採用した。
このサイトの座標のところがすごく理解の助けになった。(リンク切れ)
[3]http://www.eco.zaq.jp/env_univ/euler.html
粘性項
次に粘性であるがこれはただ単に、１つの点の速度が周囲４つ点に拡散するという意味である。どのくらいの割合で拡散するかはRe（レイノルズ数）で決まる。
vx_after[i][j]=vx[i][j]-1/Re*(vx[i+1][j]+vx[i][j+1]+vx[i-1][j]+vx[i][j-1])*Δt vy_after[i][j]=vy[i][j]-1/Re*(vy[i+1][j]+vy[i][j+1]+vy[i-1][j]+vy[i][j-1])*Δt
全計算が終わったならvx vy にvx_after vy_afterを代入する。
②外力項
マウス操作等で流体に速度をつけたいのなら、ここで速度に加算する。
また、壁と接する速度定義点は本来 速度=0.0でないといけないはずだが、前行程の移流粘性項で近隣の速度定義点から影響を受け0でなくなってしまった定義点があった場合に、ここで0に戻す。
そして流出境界などもここで設定する。
またまた東工大のサイトが参考になった！
[2.3]http://www.eto.titech.ac.jp/contents/sub03/chapter08-4.html
③ダイバージェンス計算
次に速度からダイバージェンス＝発散量を求める。
「s」はダイバージェンスを格納する配列変数。
s[i][j]=( -vx[i][j] -vy[i][j] +vx[i+1][j] +vy[i][j+1] )/Δt
ここでは1格子の空間から出ていく(or入ってくる)水の量をすべてのi,jで計算している。
スタガード格子を採用したので、この矢印を足し引きすれば発散は簡単に求まる。
この発散量を求めていく意味なのだが、記事の上の方で∇・v = 0 を連立するという話につながる。
∇・v = 0 つまり、どのi,jでも s[i][j] が0.0になる必要があると言っている。
移流や粘性、外力項などで速度情報が何度も書き換えられたため、今のままでは∇・v = 0 を満たしていないと考えられる。これでは非圧縮の条件が満たされない。
これをいっきに解決するのが次の④である。
④圧力のポアソン方程式
ここが山場だった。
上で求めた各格子のダイバージェンスから圧力を算出したい。ダイバージェンスから圧力を求めるには、以下の反復法のSOR法を使うといい。
SOR法
p[i][j]=(1.0-ω)*p[i][j] + ω/4*(p[i-1][j] + p[i+1][j] + p[i][j-1] + p[i][j+1] - s[i][j])
ωは加速係数といって自分の経験則では1.7～1.9で一番収束が早くなる。格子数が少なめなら1.7、多めなら1.9が良い(気がする)。
ΔxΔyを変えたい場合は
[2.5]http://www.eto.titech.ac.jp/contents/sub03/chapter07.html
の(7-4)式を参照のこと
これを計算領域内のi,jで、何回も反復することで収束していく。だいたい100回とか1000回とか？
当然初期p値が解から遠ければ、収束までの反復回数は増える。ここが一番計算コストのかかる部分なので、ここをいかに効率良く解くかで、全体の処理速度が決まってくるといってもいい。
反復回数は多分計算領域の大きさや許容できる誤差によるところもある。
さて、ここで２つ問題がある。
1.障害物の扱いについて
最外層以外にも任意の場所に障害物を置きたいことがある。ある流体格子の圧力を求めるとき隣の格子が壁の場合、SORの計算の前にその壁に自分の圧力をコピーしておく
薄茶色が壁、白が流体の存在するところである(さっきと色が逆ですみません)。
このとき圧力定義点の白丸はSOR法で代入も参照も行なわれる。
黄緑色の点は参照のみされる。
黒丸は参照すら行なわれない。
Bの圧力をSOR法で計算しようとする時、Bの上下左右の圧力が参照されるため、Cの圧力が参照されることになるが、Cは壁内にあるためここにはBの圧力を使う。よってBの圧力計算をする前に、CにBの圧力を代入しておく。次にAの圧力をSOR法で計算しようとする時、同じようにCにAの圧力を代入してから計算する。
結果的に黄緑色の部分の本来の圧力というものは無視される。まぁ壁の中に圧力なんてないわけなんだが。
以上の作業はSORの「ループ内」で行う必要がある。「最初に１回だけ壁に圧力をコピーしておいてあとはループ計算を流す」ではうまくいかないはず。
ここでも東工大のサイトが一番参考になった。
[2.6]http://www.eto.titech.ac.jp/contents/sub03/chapter08-2.html
この「壁」を考慮したやり方は他のサイトや書籍など調べてもあまり書いてないものの、発散から圧力を求めるときの連立一次方程式を考えると、このやり方が一番合理的だと思うので自分はこの方式を支持する。
さらにこの東工大の記事によると、格子の壁情報を格納する変数(trueなら壁 falseなら流体とか)を作ると便利だということだ。確かにそうすることで、今自分が壁なら圧力計算をskipという処理が書きやすい。
また余談になるが、SOR法は参照する４つの圧力情報のうち、計算後の情報と計算前の情報の２種類が混じり合っているように見えるがそれで良い。(ガウス・ザイデル法のアレと同じ話)
2.SOR法を何回反復するか
「p_after」に現在計算中の圧力情報を、「p」に1ループ前の圧力情報を格納しているとしてp[i][j] - p_after[i][j] の絶対値が、すべてのi,jで一定値以下なら「収束した」と判断する。ただこんな作業を毎ループやっていると、計算量2倍になり非効率なので、10ループに1回にしたり、i,jを1つ飛ばしでやるなど工夫する。
また残差計算が面倒ならループ数を固定してしまっても良いと思われる。計算の最初のフレームだけ誤差が大きそうだから最初だけ1000回ループとして、それ以降はレンダリングやゲーム等のフレームレート維持のため10ループ固定、などという案も現実的だ。
⑤修正項
圧力により速度が修正される。ここで、修正されるというのは、∇・v = 0 になるように速度が調節されるということである。
スタガード格子の利点はここにも。圧力定義点がx速度点、y速度点を挟み込むように配置されているので修正はとても簡単だ。
vx[i][j]-=( p[i][j]-p[i-1][j] )*Δt vy[i][j]-=( p[i][j]-p[i][j-1] )*Δt
⑤はこれだけ
これでやっと全ての処理が終わり。
vx vy p の値全てが次のタイムステップに進んだことになる。修正項のあとは、可視化ルーチンに回すなり、ランダムに配置した粒子の移動ルーチンに回すなり好きにすればいい。
そして処理は①へ
まとめ
こうして①～⑤を更新してSTEPを進めていくのだが、一発でプログラムが予想した挙動を示してくれるとは限らない。最後に注意すべき点がいくつかある。まずコーディングミス。実装量的にそこそこな量になるし、一つでもミスると結果がアベコベな値になる。かなり慎重に書いていく必要がある。特にxとyの取違い
次に移流項。この記事で紹介した方法ではΔtを0.1～0.2と適度に細かくしておけば発散しない。（詳しくはCFL条件で検索）
次に圧力のポアソン方程式の項。圧力pを毎STEPごとに初期化せず、前STEPのを使いまわすほうが結果的に早く収束することが分かっている(warm startingとかいう手法だった気がする)。ただしその場合、何STEPか計算していくと数値誤差でpの値が正か負に発散していくことがある。pの値は前STEPのpにのみ依存するため誤差の蓄積が原因となる。pは結局その「xy方向差分」だけを計算に使うのでpが発散するまで気づかない。なのでpを使いまわすのであれば「ガス抜き」のような処理がどこかに必要。平均値が0.0になるような処理を。
もう少し突っ込んで話すと、x,yで周期境界条件を作ったとき、圧力の逃げ場がないので数値誤差が積み重なりpが正か負に発散する。境界条件にノイマン条件がどこか1か所でもあればそこが圧力の逃げ場になり発散はしない(と自分では理解している)。
ちなみに、この流体力学のサンプルソースを配布している。
コッチ(Python版)からどうぞ。
おしまい
とりあえずNS式の意味を深く考えず、ここをこうすればいいと計算の手順の部分だけ書いてきたが、一応これで動くはずである！正直、最初に流体力学をやろうと思ったときは「∇」とか「ベクトル場」とか「ラプラシアン」とか全く分からなかったが、見よう見まねでプログラムを動かそうとしているうちにだんだんと式や単語の意味が分かってきた。私は、まずは動くプログラムを作ってから式や単語の意味を考えるのが一番効率がいいと思う。なのでとにかく動かすことを最優先で書いてきた。
2012年当時、ネットで参考になる資料はいくつかあったものの自分の理解できるレベルのものは少なく、東工大のサイトに出会わなかったらずっとCFDプログラムが作れずにいただろう。本当に感謝
参考文献
参考にしたサイトやらなにらや[1]http://ja.wikipedia.org/wiki/%E3%83%8A%E3%83%93%E3%82%A8-%E3%82%B9%E3%83%88%E3%83%BC%E3%82%AF%E3%82%B9%E6%96%B9%E7%A8%8B%E5%BC%8F
[2]http://www.eto.titech.ac.jp/contents/
[3]http://www.eco.zaq.jp/env_univ/euler.html
[4]http://d.hatena.ne.jp/etopirika5/20110425/1303740299
[5]http://www.bee-www.com/log.htm
[6]http://journal.mycom.co.jp/articles/2007/10/09/cedec01/001.html
[7]http://www.4gamer.net/games/000/G000000/20070926041/screenshot.html?num=012
[8]http://fluid.me.seikei.ac.jp/lecture/cfd/cfd-08-NS-algorithm.pdf
[9]http://staff.aist.go.jp/kogaki.t/Dissertation/PDF_files/chap3.pdf
[10]http://homepage.mac.com/catincat/java/index.html
[11]http://oshiro.bpe.es.osaka-u.ac.jp/people/staff/imura/lecture/2007cg/fire071220.pdf
[12]http://www.youtube.com/watch?v=6rPL-QkUFf8&feature=related
[13]http://www.youtube.com/watch?v=4Ov6tguRr8E
他、個人ブログの方が書いたゲームetcに格子法流体を実装するのに参考になりそうなサイト(2020/6更新)
http://el-ement.com/blog/2015/03/16/gpu-fluid/
https://edom18.hateblo.jp/entry/2020/02/13/104744
2018/10追記
東工大の流体のページがいつからか閉鎖させておりアクセスできない状態になっていた。ただWARPという国内の教育機関などに特化したwebアーカイブサイトに保存されていた！！
見たい方はコチラで「コンピュータシミュレーションによる流体解析（発展編）」と検索
アーカイブ自体は私も保存しているのでもしほしい方はコメントを(教育目的のサイトだし著作権的に個人所有は問題ないはず・・)
次回：流体力学のプログラムを作りたい！その２
はじめまして、ナビエストークスの技術資料を探しもとめてこのサイトの情報に出会えました。
アメブロ読んでたんですがね、全くわかりませんでしたが。
こっちにも顔出すようにします。