1から始める流体シミュレーション~弐~

はじめに

こんにちは~

こんにちは、tardigradeです。本記事が流体シミュレーションについての2本目の記事となります。

前回の記事はこちら↓

tardigradecfd.hatenablog.com

高校生が隙間時間を使ってだらっと書いている記事ですので、あまり硬くならずに読んでいただければ幸いです。また、書いている記事の内容の正確性については保証しません。気になる間違いがございましたら、コメントでご指摘いただければ修正いたします。

 

粒子法について語るよ

前回までの記事の最後で、粒子法にも実はいくつかの種類があるんだよ、ということを言いました。本記事では、その辺について詳しく、また実装の方法なんかについてもより具体的な話ができたらと思います。

実装についての具体的なお話

どの方法を使うのか

前述したとおり粒子法にもいくつか種類があるわけですが、それぞれで「出来ること・出来ないこと」、「得意なこと・不得意なこと」が違ってきます。グーグル先生と相談した結果、今回作成する流体シミュレーションは、粒子法の一種である「MPS法」を用いて実装することにしました。MPS法を選んだ理由は、「比較的実装が楽そうだったから」、「参考にできそうな文献がいくつか見つかったから」、「水のシミュレーション専用の手法だから」等です。主観的な理由が多いですが、ぶっちゃけどの方法使っても完璧なシミュレーションは無理なので、最後は気分で決めちゃいました。ちなみに、他の手法で有名どころだと「SPH法」なんかが挙げられます。MPS法と違ってSPH法は気体のシミュレーションなんかもできる万能な手法です。MPS法の実装がひと段落したら手を出してみたい…|ω・)

可視化はどうするの?

計算自体はMPS法でガチャガチャやってしまいますが、最終的には結果が可視化がされないと困ります。可視化までが目標と前回定めたわけですし、そもそも大量の数値を見せられて「これが結果だよ☆」とか言われても僕はコンピューターではないので理解不能です。そこで今回は「Paraview」というアプリを使って結果を可視化しようと考えています。ほかにも方法はいろいろあると思うのですが、

  • 可視化に計算量をあまり取られない
  • 実装があまり手間じゃない
  • 学習コストが安い

の三点で絞った結果、これが一番楽かなぁという感じになりました。特に最後の学習コストは結構大事で、僕はそこまでプログラミングつよつよ勢ではないので、初心者でもできそうな方法を選びました。

Paraviewの具体的な使い方は実際にコーディングしていくときに紹介していきたいと思います(まだよくわかってないとか言えない…)。

MPS法の気持ちを考える~前提知識part~

ここからは、具体的にMPS法がどういった考え方なのかを紹介していきたいと思います。僕は流体力学の専門家ではないので、理解に間違いがあるかもしれませんが、あらかじめご了承ください。

離散化という考え方

 MPS法では、流体を細かい粒子の集合体とみなして計算をします。より正確には、「連続していて、力を加えると変形してしまう流体」を、「不連続で、力を加えても変形しない細かい粒子の集合体」に変換します。なぜわざわざそんなことをするのかというと、コンピューターは連続しているものを計算するのが苦手だからです。

f:id:tardigrade:20210110002529p:plain

数学風に言えば「 ∫ 」を「 Σ 」に変換しているようなイメージ

 

高校物理を学習中の身からすると、「流体は力を加えても変形しちゃうし、どこまでで一個かもわからないからf=maが使えない( ノД`)シクシク…」とずっと思っていたので、この考え方は割と納得です。
このように、連続したものを不連続なものの集合で近似して、コンピューターで計算しやすい形に直す手法を「 離散化 」と呼びます。身近な例だと、アニメーションや動画は画像を非常に短い間隔で切り替えて滑らかな映像に見せていますが、あれも時間を離散化しているというように言えると思います。MPS法に限らず、流体シミュレーションでは「計算対象となる物体」「空間」「時間」「支配方程式」「関数」等々を片っ端から離散化していくのが普通です。

流体の支配方程式

さて、計算対象を離散化してコンピューターで計算しやすくしたわけですが、結局、細かく分けた粒子一つ一つの物理量(速度、圧力等)が求められないとシミュレーションができません。ここで登場するのが「ナビエストークス方程式」です。こいつを使って各粒子の位置や速度を計算していきます。平たく言えば流体専用の運動方程式ですね。ピースサインが可愛いことで有名です。

f:id:tardigrade:20210111023822p:plain

運動方程式の比較

高校で習う運動方程式と比較するとそのすさまじさがわかりますね…。この方程式の導出から理解したい方は、ヨビノリ先生がYouTubeに解説動画を上げているのでそちらをご覧ください。ナビエストークス方程式の動画を見る前に、「偏微分」「全微分」「テイラー展開」「ベクトル解析」を理解していると話が早いと思います(これもヨビノリ先生が動画を上げてくださっています)。

ちなみに、僕は2週間近くかけて一通り動画を見ました。「完全に理解した」とは到底言えないですが、見てよかったと思っています。知らなくても実装ができないわけではないですが、モチベ維持にもつながりますし、メリットはそれなりにある気がします。

離散化しなきゃ…

上でナビエストークス方程式について紹介しましたが、実際にこの方程式を計算する際には、方程式を離散化をする必要があります。より具体的には、方程式に含まれているベクトル微分演算子(三角形・逆三角形の記号)を離散化が必要になります。

「離散化が必要?つまりはなんか連続してるってこと?」、はい、その通りです。このベクトル微分演算子がどういうものなのか、わからない人向けに、誤解を恐れずに言うならば、「周囲の空間の状態からある点の状態を求める」ような演算子です。あんまりピンと来ない説明かもしれませんが、今はそういうもんなんだと思ってください。数学的、物理的にはあまり正確な言い方ではないと自分でも思います…。

話を戻しますが、このベクトル微分演算子をコンピューターでそのまま計算させることはできません。なぜでしょうか?ポイントは「周囲の空間」です。空間とは連続しているものなので、コンピューターでは上手く扱えず、よってそのままの形ではベクトル微分演算子を計算することができないのです。

これで、ナビエストークス方程式を離散化する必要性を理解いただけたでしょうか…?

ベクトル微分演算子についてより深く知ろう

ベクトル微分演算子を離散化したいわけですが、その方法を理解するにはベクトル微分演算子の性質をより深く知らなければなりません。これもめちゃめちゃざっくりとした言い方になってしまいますが、ベクトル微分演算子がある点の状態を求めるとき、近くの空間の状態は色濃く、逆に遠くの空間の状態ほど薄く反映します。

例えば、密閉された室内のある点における温度を求めることを考えてみましょう。ある点の周りの空間の温度がわかっている場合、そのある点の温度も大体予想がつくことは理解できるかと思います。

f:id:tardigrade:20210111150011p:plain

温度の予想

このとき、周囲の空間の温度を一部変化させることを考えます。求めたい点に近い位置の温度をがくんと下げた場合、点の温度もがくんと下がってしまいますよね。しかし、求めたい点からかなり離れた位置の温度を下げた場合はどうでしょうか?どう考えても近い位置の温度を下げるより、影響は少ないはずです。

ベクトル微分演算子もこんな感じで、近くの空間からの影響は強く受けますが、遠くの空間からの影響になるほど弱く、無限に近い距離の空間からの影響はほぼ無視できるほどになのです。

重み付き平均と影響半径

ベクトル微分演算子について詳しく知れたところで、いよいよ離散化の方法に移ります。

ベクトル微分演算子の離散化に登場するのが「 重み付き平均 」という考え方です。「重み付き平均」でググるとほかにわかりやすいサイトが見つかると思いますが、簡単に定義を説明すると、「各値の重み(価値や重要さとも言い換えられる)を考慮したうえでとった平均」です。

実際に複数の粒子を使って例えると、「粒子i の物理量は、i の周りに存在する粒子(以下、近傍粒子) の物理量の平均値に近くなるようにしよーよ。そんとき近くにいるに粒子の値を重くとってね」という感じです。なんとなくベクトル微分演算子と似ているのがわかるでしょうか?

f:id:tardigrade:20210111021143p:plain

重み付き平均を取ると…

 

ちなみに、どんなふうに重みを取るのかを表した関数を「重み関数」と言います。

また、重み付き平均を取る際には「影響半径」というものを定義します。

f:id:tardigrade:20210111033230p:plain

影響半径

上のイラストのピンクの円の半径が影響半径です。一般に影響半径内にある粒子を、粒子i の「近傍粒子」と言います。

だんだん話がこんがらがってきましたね…。僕も書いててパンクしそうです…w

重み付き平均は近傍粒子からとるので、影響半径を定義することは、重み付き平均を取る際に考える粒子の数を制限することに等しいです。同時に、影響半径外の粒子からの影響はないものとみなすということでもあります。

実はこれで…ベクトル微分演算子の離散化ができています…。

ナビエストークス方程式の離散化まとめ

さて、話が長くなったのでまとめたいと思います。

  1. 流体を粒子に分割
  2. 複数の粒子に分けて考える
  3. 各粒子の物理量はナビエストークス方程式から導出する
  4. ナビエストークス方程式は離散化しなければ計算できない
  5. より具体的にはベクトル微分演算子の離散化
  6. ベクトル微分演算子は周囲の状態から自身の状態を評価するもの
  7. 重み付き平均という考え方を使って離散化をする
  8. 重み関数と影響半径を定義してあれこれする
  9. あら不思議、離散化ができちゃいました

ざっとこんな感じでしょうか。最後の部分、不満がある方が多いと思いますので補足説明をすると、もともと離散化する前のベクトル微分演算子は「空間(0~∞までの距離に存在するあらゆる点)の状態」から自身の状態を求めていたのに対し、重み付き平均を使った場合だと「影響半径内に存在する近傍粒子の状態」から自身の状態を求めています。離散化とは要するに不連続なもので近似するということでしたので、なんとなくこれでできているがしませんか…?すみませんこれが僕にできる精いっぱいの説明です_(._.)_

次回予告

次回は具体的な式と、大まかなMPS法の計算ステップについて説明出来たらいいなぁと思います。

わかりにくい部分が多かったかと思いますが、最後まで読んでいただきありがとうございました。