1から始める流体シミュレーション~弐~
はじめに
こんにちは～
こんにちは、tardigradeです。本記事が流体シミュレーションについての2本目の記事となります。
前回の記事はこちら↓
高校生が隙間時間を使ってだらっと書いている記事ですので、あまり硬くならずに読んでいただければ幸いです。また、書いている記事の内容の正確性については保証しません。気になる間違いがございましたら、コメントかtwitterでご指摘いただければ修正いたします。
粒子法について語るよ
前回までの記事の最後で、粒子法にも実はいくつかの種類があるんだよ、ということを言いました。本記事では、その辺について詳しく、また実装の方法なんかについてもより具体的な話ができたらと思います。
実装についての具体的なお話
どの方法を使うのか
前述したとおり粒子法にもいくつか種類があるわけですが、それぞれで「出来ること・出来ないこと」、「得意なこと・不得意なこと」が違ってきます。グーグル先生と相談した結果、今回作成する流体シミュレーションは、粒子法の一種である「MPS法」を用いて実装することにしました。MPS法を選んだ理由は、「比較的実装が楽そうだったから」、「参考にできそうな文献がいくつか見つかったから」、「水のシミュレーション専用の手法だから」等です。主観的な理由が多いですが、ぶっちゃけどの方法使っても完璧なシミュレーションは無理なので、最後は気分で決めちゃいました。ちなみに、他の手法で有名どころだと「SPH法」なんかが挙げられます。MPS法と違ってSPH法は気体のシミュレーションなんかもできる万能な手法です。MPS法の実装がひと段落したら手を出してみたい…|ω・)
可視化はどうするの？
計算自体はMPS法でガチャガチャやってしまいますが、最終的には結果が可視化がされないと困ります。可視化までが目標と前回定めたわけですし、そもそも大量の数値を見せられて「これが結果だよ☆」とか言われても僕はコンピューターではないので理解不能です。そこで今回は「Paraview」というアプリを使って結果を可視化しようと考えています。ほかにも方法はいろいろあると思うのですが、
- 可視化に計算量をあまり取られない
- 実装があまり手間じゃない
- 学習コストが安い
の三点で絞った結果、これが一番楽かなぁという感じになりました。特に最後の学習コストは結構大事で、僕はそこまでプログラミングつよつよ勢ではないので、初心者でもできそうな方法を選びました。
Paraviewの具体的な使い方は実際にコーディングしていくときに紹介していきたいと思います(まだよくわかってないとか言えない…)。
MPS法の気持ちを考える~前提知識part~
ここからは、具体的にMPS法がどういった考え方なのかを紹介していきたいと思います。僕は流体力学の専門家ではないので、理解に間違いがあるかもしれませんが、あらかじめご了承ください。
離散化という考え方
MPS法では、流体を細かい粒子の集合体とみなして計算をします。より正確には、「連続していて、力を加えると変形してしまう流体」を、「不連続で、力を加えても変形しない細かい粒子の集合体」に変換します。なぜわざわざそんなことをするのかというと、コンピューターは連続しているものを計算するのが苦手だからです。
高校物理を学習中の身からすると、「流体は力を加えても変形しちゃうし、どこまでで一個かもわからないからf=maが使えない( ﾉД`)ｼｸｼｸ…」とずっと思っていたので、この考え方は割と納得です。
このように、連続したものを不連続なものの集合で近似して、コンピューターで計算しやすい形に直す手法を「 離散化 」と呼びます。身近な例だと、アニメーションや動画は画像を非常に短い間隔で切り替えて滑らかな映像に見せていますが、あれも時間を離散化しているというように言えると思います。MPS法に限らず、流体シミュレーションでは「計算対象となる物体」「空間」「時間」「支配方程式」「関数」等々を片っ端から離散化していくのが普通です。
流体の支配方程式
さて、計算対象を離散化してコンピューターで計算しやすくしたわけですが、結局、細かく分けた粒子一つ一つの物理量(速度、圧力等)が求められないとシミュレーションができません。ここで登場するのが「ナビエストークス方程式」です。こいつを使って各粒子の位置や速度を計算していきます。平たく言えば流体専用の運動方程式ですね。ピースサインが可愛いことで有名です。
高校で習う運動方程式と比較するとそのすさまじさがわかりますね…。この方程式の導出から理解したい方は、ヨビノリ先生がYouTubeに解説動画を上げているのでそちらをご覧ください。1)ナビエストークス方程式の動画を見る前に、「偏微分」「全微分」「テイラー展開」「ベクトル解析」を理解していると話が早いと思います(これもヨビノリ先生が動画を上げてくださっています)。
ちなみに、僕は2週間近くかけて一通り動画を見ました。「完全に理解した」とは到底言えないですが、見てよかったと思っています。知らなくても実装ができないわけではないですが、モチベ維持にもつながりますし、メリットはそれなりにある気がします。
離散化しなきゃ…
上でナビエストークス方程式について紹介しましたが、実際にこの方程式を計算する際には、方程式を離散化をする必要があります。より具体的には、方程式に含まれているベクトル微分演算子(三角形・逆三角形の記号)を離散化が必要になります。
「離散化が必要？つまりはなんか連続してるってこと？」、はい、その通りです。このベクトル微分演算子がどういうものなのか、わからない人向けに、誤解を恐れずに言うならば、「周囲の空間の状態からある点の状態を求める」ような演算子です。あんまりピンと来ない説明かもしれませんが、今はそういうもんなんだと思ってください。数学的、物理的にはあまり正確な言い方ではないと自分でも思います…。
話を戻しますが、このベクトル微分演算子をコンピューターでそのまま計算させることはできません。なぜでしょうか？ポイントは「周囲の空間」です。空間とは連続しているものなので、コンピューターでは上手く扱えず、よってそのままの形ではベクトル微分演算子を計算することができないのです。
これで、ナビエストークス方程式を離散化する必要性を理解いただけたでしょうか…？
ベクトル微分演算子についてより深く知ろう
ベクトル微分演算子を離散化したいわけですが、その方法を理解するにはベクトル微分演算子の性質をより深く知らなければなりません。これもめちゃめちゃざっくりとした言い方になってしまいますが、ベクトル微分演算子がある点の状態を求めるとき、近くの空間の状態は色濃く、逆に遠くの空間の状態ほど薄く反映します。
例えば、密閉された室内のある点における温度を求めることを考えてみましょう。ある点の周りの空間の温度がわかっている場合、そのある点の温度も大体予想がつくことは理解できるかと思います。
このとき、周囲の空間の温度を一部変化させることを考えます。求めたい点に近い位置の温度をがくんと下げた場合、点の温度もがくんと下がってしまいますよね。しかし、求めたい点からかなり離れた位置の温度を下げた場合はどうでしょうか？どう考えても近い位置の温度を下げるより、影響は少ないはずです。
ベクトル微分演算子もこんな感じで、近くの空間からの影響は強く受けますが、遠くの空間からの影響になるほど弱く、無限に近い距離の空間からの影響はほぼ無視できるほどになるのです。
重み付き平均と影響半径
ベクトル微分演算子について詳しく知れたところで、いよいよ離散化の方法に移ります。
ベクトル微分演算子の離散化に登場するのが「 重み付き平均 」という考え方です。「重み付き平均」でググるとほかにわかりやすいサイトが見つかると思いますが、簡単に定義を説明すると、「各値の重み(価値や重要さとも言い換えられる)を考慮したうえでとった平均」です。
実際に複数の粒子を使って例えると、「粒子 の物理量は、 の周りに存在する粒子(以下、近傍粒子) の物理量の平均値に近くなるようにしよーよ。そんとき近くにいるに粒子の値を重くとってね」という感じです。なんとなくベクトル微分演算子と似ているのがわかるでしょうか？
ちなみに、どんなふうに重みを取るのかを表した関数を「重み関数」と言います。
また、重み付き平均を取る際には「影響半径」というものを定義します。
上のイラストのピンクの円の半径が影響半径です。一般に影響半径内にある粒子を、粒子 の「近傍粒子」と言います。
だんだん話がこんがらがってきましたね…。僕も書いててパンクしそうです…w
重み付き平均は近傍粒子からとるので、影響半径を定義することは、重み付き平均を取る際に考える粒子の数を制限することに等しいです。同時に、影響半径外の粒子からの影響はないものとみなすということでもあります。
実はこれらの考え方をうまく用いることで…ベクトル微分演算子の離散化ができるのです…。
ナビエストークス方程式の離散化まとめ
さて、話が長くなったのでまとめたいと思います。
- 流体を複数の粒子に分割する
- 各粒子の物理量はナビエストークス方程式から導出する
- ナビエストークス方程式は離散化しなければ計算できない
- より具体的にはベクトル微分演算子の離散化が必要
- ベクトル微分演算子は周囲の状態から自身の状態を評価するもの
- 重み付き平均という考え方を使って離散化する
- 重み関数と影響半径というものを定義する
- 7を使うことでうまい感じに離散化できる
ざっとこんな感じでしょうか。最後の部分、不満がある方が多いと思いますので補足説明をすると、もともと離散化する前のベクトル微分演算子は「空間(0～∞までの距離に存在するあらゆる点)の状態」から自身の状態を求めていたのに対し、重み付き平均を使った場合だと「影響半径内に存在する近傍粒子の状態」から自身の状態を求めています。離散化とは要するに不連続なもので近似するということでしたので、なんとなくこれで近似できるような気がしませんか…？すみませんこれが僕にできる精いっぱいの説明です_(._.)_
次回予告
次回は重み関数やベクトル微分演算子を離散化した、具体的な式と、大まかなMPS法の計算ステップについて説明出来たらいいなぁと思います。
わかりにくい部分が多かったかと思いますが、最後まで読んでいただきありがとうございました。
引用・参考文献
1)
2)
https://core.ac.uk/download/pdf/147692524.pdf
1から始める流体シミュレーション~壱~
はじめに
はじめまして
こんにちは、tardigradeです。本記事が初投稿となります。高校生が隙間時間を使ってだらっと書いている記事ですので、あまり硬くならずに読んでいただければ幸いです。また、書いている記事の内容の正確性については保証しません。気になる間違いがございましたら、コメントかtwitterでご指摘いただければ修正いたします。
流体シミュレーションを実装したい
ある日突然、どうしても流体シミュレーションを実装してみたい気分になりました。CG用じゃなくで解析用で使うようなガチガチのやつです。どうにもこの衝動が抑えられそうになかったので、試しに「流体シミュレーション 方法」でググってみました。
……「粒子法」、「格子ボルツマン法」、「ラグランジュ微分」、どのサイトにも難しそうな用語が並んでいて一筋縄ではいかなそうです。まぁ前提として、流体力学や大学数学を学んでいないのに、流体シミュレーションを理解して実装しようというのは無理があるのですが。とはいえこのまま引き下がることなどできません。僕はどうしても流体シミュレーションを実装したいのです。
ブログを書くか
根気よくいろいろなサイトを読み進めていった結果、なんとなく実装の雰囲気だけはつかめた気がしました。あくまでも「気がした」だけです…。まぁ細かいことは置いといて一回は実装してみるべきだなと思い、パソコンと向かいあったのですが、ふと「これ自分の中だけで完結させるのもったいなぁ」という考えが頭をよぎりました。そこで、将来僕と同じような衝動に襲われるかもしれない中高生のために、ブログにその実装方法や考え方を詳細に記録していくことに決めました。僕としても考えを整理できるし、一石二鳥です。ちなみにこれを書いている時点で、まだコードは一行も書いていないので、この先きちんと実装しきれるかどうかは未定です…。あしからずご了承ください<(_ _)>
実装の方針について考える
どんなシミュレーションをしたい?
何をするにもゴールを先に決めることは重要です。「流体シミュレーション」と聞いてぱっと思いつくのが、水柱(某有名作品とは関係ないです)が崩壊するやつだったので、今回はそれを実装してみます。数値だけの計算ではなく、計算結果を可視化するところまでが目標です。いまいちイメージがつかめない方は「流体シミュレーション」で検索すると下の図に似た画像が見つかると思います。
「流体」というと水以外にもケチャップやスライムや、なんなら気体も流体ですが、それらのシミュレーションはまた別の機会にということで。
水のシミュレーションの手法
色々調べてみたところ、大きく分けて「格子法」と「粒子法」の2種類があるみたいです。具体的にそれらがどういう方法なのかはとりあえず置いといて、両者のメリット・デメリットについて比較してみました。1)
格子法
- シンプルで堅牢性があり高速
- リアルタイム実装向きである
- 物理的な正確さに欠ける
粒子法
- 大変形に強い
- 剛体との連成解析が比較的楽
- 計算量が比較的多い
ざっくり言うと、格子法は精度が微妙だけど計算は軽い、粒子法は精度はいい感じだけど計算が重い、という感じでしょうか。ちなみに、最近だと両者を混ぜた手法なんかも考案されているようですが、参考文献があまり見つからなかったのでここでは割愛します。
さて、どちらを採用するのかという話ですが、個人的に格子法の「物理的な正確さに欠ける」という点が気になったので、今回は粒子法で実装することにしました。時間があれば、格子法も実装してみて両者の比較をしてみたいですね。
次回予告
ところで、実は粒子法にもいくつか種類があります。次回のブログでは、具体的にどの方法を使うのか、そしてそれはどのような考え方でどのように実装するのかを書いていきたいと思います。
