この記事はtraP Advent Calendar2017の11/17の記事だよ❤️
というわけでｷｬｯﾙﾙｰﾝ♪そばやだよ～☆

今回はｱﾄﾞﾍﾞﾝﾄｶﾚﾝﾀﾞｰということでぇ　流体シミュやってみよっかな～って思いましたぁ

というノリを期待されてたような気もするのですが、疲れるので普通に書きます。

以前から興味のあった流体をちょいと本読んでやってみました。

できあがったものはこちら。

ソースはこちら。

以下は解説です。

流体シミュレーションって？

流体シミュレーションとは、難しい方程式を解かないとわからないような流体の動きをプログラムで近似するものです。
剛体シミュレーションで運動方程式

$$F=ma$$

を用いるところを、流体シミュレーションではナビエ・ストークス方程式(以下NS方程式)というものを使います。その昔、ナビエさんとストークスさんが発明したらしいです。
wikipedia曰く、その全貌は以下のようになっています。

$$\frac{Dv}{Dt}=-\frac{1}{\rho }grad(p)+\frac{\mu }{\rho }\Delta v+\frac{\lambda +\mu }{\rho }grad(\Theta )+\frac{\Theta }{\rho }grad(\lambda +\mu )+\frac{1}{\rho }grad(v\cdot grad(\mu ))+\frac{1}{\rho }rot(v\times grad(\mu ))-\frac{1}{\rho }v\Delta \mu +g$$

わわわからぁぁぁぁぁん!!!!
何言ってんだこいつは！ラグランジュ微分ってなんだ！ふざけるな！！
ということで限定条件下で式を簡略化します。
今回使う条件は「粘性率が一定の非圧縮性流れ」です。
粘性率とは、まぁ粘り気です。
非圧縮性とは、どんなに力を加えても流体の体積が一定に保たれるというものです。(現実にはありえません)
すると使う式は

$$\frac{Dv}{Dt}=-\frac{1}{\rho }grad(p)+\nu \Delta v+g$$

となります。
これならなんとか食らいつけそうです。1項ずつ見ていきましょう。

$\frac{Dv}{Dt}$
$\frac{Dv}{Dt}$はラグランジュ微分なのですが、まぁとりあえず普通に時間微分だと思って良いでしょう。$v$は速度なので、左辺はいわゆる加速度です。

$-\frac{1}{\rho }grad(p)$
右辺第一項は圧力勾配項と呼ばれ、圧力の大きいところから小さいところにむけて加速度を作り出す作用のある項です。$\rho$は密度ですが、実質質量ですね。力を質量で割って加速度を得るようなかんじです。

$\nu \Delta v$
右辺第二項は粘性項と呼ばれ、文字通り粘り気のある動きをさせるための項です。
より具体的には、周囲の速度に自分の速度を近づけるような作用をします。$\nu$は粘性係数と呼ばれる定数です。

$g$
右辺第三項は外力項です。
私はそばやであって物理屋ではないので、詳しい内容はわからないのですが、噛み砕くと重力です。

と、以上の計算を以て加速度が割り出されます。
結局流体でも「なんらかの計算によって加速度が定まる」という方式は通用するのです。
なんかできるような気がしてくるでしょ？

ところで、NS方程式って何にたいして適用されるのでしょう？

剛体の場合には「これが1つの物体です！！」と明確に言えるものがありますよね。
ですからその物体の質量だとか速度だとかを議論するのは自然な発想でした。
ですが流体の場合には、「これこそが流体の1単位です！！」といえるようなものがあるでしょうか？ここには議論の余地があって、これによって流体シミュレーションには大きく分けて2つの手法が存在しています。

流体シミュレーションの種類

流体シミュレーションの手法は大きく分けて2種類。

1つは格子法です。
格子法とは、計算する空間を適当な数の格子で分割し、その小領域それぞれに対し速度などのパラメータをあてていきます。
つまり、「1つの格子が流体の1単位です！」という主張なわけですね。
実際、NS方程式は速度場に対する加速度などを考えるものなので、これは結構自然な発想だと思います。

もう1つは粒子法です。今回テーマにしたいのはこちらですね。
粒子法は、流体を粒子の集合として表現する方法です。
粒子1つ1つに速度などが存在します。
これもなんとなく納得できそうなかんじです。

どっちがいいの？

2種類あると比較したくなるのが人情ってもんです。

まず格子法ですが、こちらは数学的な「場」の概念に則った自然なモデルです。
ですから、(こういうとおかしいかもしれませんが)物理法則が通用します。
計算もとても安定しています。
欠点としては

• 完全に流体の存在する領域を決定しないと使えない
• まともな見た目のものを作ると、計算負荷が高い
• 水面を作りづらい

などが考えられます。

一方粒子法です。
粒子法について詳細の説明はこれからとなりますが、私の感想としては「粒子法は格子法に比べて、だいぶ無理がある」といったかんじです。
若干こじつけじみたところがあり、そのせいで数値的安定性を欠く傾向があります。
もちろんそのぶん得られるメリットもあります。
メリットとしては

• 存在する場所を予め決めておかなくても良い
• 計算負荷の調整がある程度容易
• 水面や飛沫が表現できる

といったところです。
水面が動く部分や水飛沫などは見た目に面白いので、粒子法のほうがエンターテインメント性に富んでいると思います。

いいから早くやれ

はいでは粒子法やっていきます。
で、さらに分岐があるんですね。
一口に粒子法といってもいくつか種類があるようで、今回はMoving Particle Semi-implicit(以下MPS法)という手法を使います。

Moving 動くんですね。

Particle 粒子が。まぁ粒子法ですしね。

Semi-Implicit ココが肝なんですよね。あとでご説明します。

事前準備

長いですね～。
まぁ準備は必要なんでしますね。
まず最初に考えなきゃいけないのは、「物理的なパラメータをどうやって粒子に紐付けるか」です。例えばなんですけど、NS方程式の中に$grad$とか$\Delta$とかありましたよね？$grad$ってのは変数を$dx$とか$dy$とかで微分したものなので、その粒子の位置におけるパラメータそのものではなく、その粒子の付近でのパラメータの上下関係が知りたいわけです。
そうすると本来なら「自分のパラメータと隣のパラメータの差分」をとるのが良いわけじゃないですか。実際、格子法ではまさにそうやって$grad$とかを計算します。
ですが、粒子法の場合には違います。粒子法の場合は隣に粒子がいないかもしれません。それに、格子と違って粒子は乱雑に並んでいるので、そもそも「隣」ってなんだよって話です。
そこで、粒子法では(あいまいですが)「近傍」の粒子を全部みてそこから勾配とか発散とかを計算します。
ただし、いくら近傍とはいえ距離が遠いものの寄与が少ないほうがよいので、なんかいいかんじの重み関数を使って平均化します。
重み関数は距離が遠ければ遠いほど小さくなると良いので、MPSでは以下を使います。

$$w(r)=\{\begin{array}{l}{r}_e/r-1(r<r_e)\\ 0(r\ge r_e)\end{array}$$

だいたい距離$r$に反比例するような式ですね。これを用いて、物理量$\varphi$の勾配は以下のように定義されます。

$${⟨\nabla \varphi ⟩}_i=\frac{d}{n^0}\sum _{j\ne i}\left[\frac{{\varphi }_j-{\varphi }_i}{|{\mathbf{r}}_j-{\mathbf{r}}_i{|}^2}({\mathbf{r}}_j-{\mathbf{r}}_i)w(|{\mathbf{r}}_j-{\mathbf{r}}_i|)\right]$$

同様に、発散は以下

$${⟨\nabla \cdot \mathbf{u}⟩}_i=\frac{2d}{n^0}\sum _{j\ne i}\frac{{\mathbf{u}}_{ij}\cdot ({\mathbf{r}}_j-{\mathbf{r}}_i)}{|{\mathbf{r}}_j-{\mathbf{r}}_i{|}^2}w(|{\mathbf{r}}_j-{\mathbf{r}}_i|)$$

同様に、ラプラシアンは以下

$${⟨{\nabla }^2\varphi ⟩}_i=\frac{2d}{\lambda n^0}\sum _{j\ne i}\left[({\varphi }_j-{\varphi }_i)w(|{\mathbf{r}}_j-{\mathbf{r}}_i|)\right]$$

ただし、$\lambda ,n^0$は初期粒子配置で十分内部にある粒子において以下のように計算したものを用います。

$$\lambda =\frac{{\sum }_{j\ne i}|{\mathbf{r}}_j-{\mathbf{r}}_i{|}^{2}w(|{\mathbf{r}}_j-{\mathbf{r}}_i|)}{{\sum }_{j\ne i}w(|{\mathbf{r}}_j-{\mathbf{r}}_i|)}$$

$$n^0=\sum _{j\ne i}w(|{\mathbf{r}}_i-{\mathbf{r}}_j|)$$

それと、$d$は次元です。今回は2次元のシミュレーションをするので、ずっと$d=2$です。

途端になんかややこしい式が出てきましたが、基本的には

1. 周囲の粒子との物理量の差分をとる
2. 重み関数で平均化する
3. 定数倍する

みたいなステップで計算されます。

「初期配置状態で十分内部の」とかいうあいまいな文言が出てきますが、これは「安定状態のときの粒子の状態」を保存するような気持ちがあるようです。
粒子の動きが安定化してきたらおおよそ初期配置状態と同じような密度になるんじゃね？ということらしいです。

ともあれ、これで準備完了です。

MPS法のタイムステップ

非圧縮性流れにおける支配方程式は2つです。

• $\frac{D\rho }{Dt}=0$
• $\frac{D\mathbf{u}}{Dt}=-\frac{1}{\rho }\nabla P+\nu {\nabla }^2\mathbf{u}+\mathbf{g}$
  都合により、速度は$\mathbf{u}$と書かせていただきます。
  下の式はすでにご説明したNS方程式ですが、上の式はなんでしょう？
  密度の時間微分が$0$、これは非圧縮性を示しています。
  圧縮したら一部の密度だけ増えますもんね。

これらから、各粒子に対する加速度を計算し、速度と位置を更新していきます。
ここでポイントなのが、MPS法ではNS方程式の右辺を半陰的(Semi-Implicit)に解きます。
半陰的ってなんだよ？ってなりますよね。
数値計算では、いま自分が$k$フレーム目にいるとして$k+1$フレーム目の状態を計算するときに、$k$フレーム目の値を使うことを陽的、$k+1$フレーム目の値を使うことを陰的といいます。
MPS法ではこれらを混ぜて使います。
すなわち、NS方程式第2,3項は陽的に、NS方程式第1項と非圧縮の式は陰的に解きます。

$${\left[\frac{D\rho }{Dt}\right]}^{k+1}=0$$

$$\frac{D\mathbf{u}}{Dt}={\left[-\frac{1}{\rho }\nabla P\right]}^{k+1}+{\left[\nu {\nabla }^2\mathbf{u}\right]}^k+{\left[\mathbf{g}\right]}^k$$

陰的に解くと陽的に解くのに比べ、計算誤差により値が発散するのを防ぐ効果があります。ですが、陰的に解くほうが大変です。
われわれは$k$ステップ目の値は知っていますが、$k+1$ステップ目の値はまだ計算していないからです。
とりあえず$\frac{D\mathbf{u}}{Dt}$を時間で離散化して、$\frac{\Delta \mathbf{u}}{\Delta t}$としましょう。
さらにそれを2つに分けて、

$$\frac{\Delta {\mathbf{u}}^{\ast }}{\Delta t}={\left[\nu {\nabla }^2\mathbf{u}\right]}^k+{\left[\mathbf{g}\right]}^k$$

$$\frac{\Delta {\mathbf{u}}^{\prime }}{\Delta t}={\left[-\frac{1}{\rho }\nabla P\right]}^{k+1}$$

とします。
まず、上の式は未知数が$\Delta {\mathbf{u}}^{\ast }$のみなので普通に計算できます。
計算しましょう。
ついでに位置も動かしときましょう。

$${\mathbf{r}}^{\ast }=\mathbf{r}+\Delta {\mathbf{u}}^{\ast }\Delta t$$

下の式ですが、未知数が$\Delta {\mathbf{u}}^{\prime }$と${\left[\nabla P\right]}^{k+1}$なので計算できません。
先程のラプラシアンの式を使ったり非圧縮の式を使ったりごにょごにょすると、以下のようになります。

$${\nabla }^{2}P=\frac{2d}{\lambda n^0}\sum _{j\ne i}\left[(P_j^{k+1}-P_i^{k+1})w(|{\mathbf{r}}_j^{\ast }-{\mathbf{r}}_j^{\ast }|)\right]=-\frac{{\rho }_0}{\Delta t^2}\frac{n^{\ast }-n^0}{n^0}$$

$n^{\ast }$は粒子数密度というやつで、

$$n^{\ast }=\sum _{j\ne i}w(|{\mathbf{r}}_i^{\ast }-{\mathbf{r}}_j^{\ast }|)$$

です。重みの和です。粒子が集中している場所では大きくなる値です。
ともあれ、式をよく見てみると、未知数は$P_i^{k+1}$のみです。しかも線形です。
$P_i^{k+1}$は粒子$i$のところでの$k+1$ステップ目の圧力なので、結局これは連立一次方程式ってことになります。未知数の数と式の数は合っているのでいけそうな気がします。

って思ったのですが、実はこれ、解けないです。（これは本に書いてなかった)
というのも、左辺をすべての粒子に対して足すと、常に$0$になるのですが、右辺は一般に$0$になりません。したがって、この連立方程式は一般に解を持ちません。
ゆえに、これを解くのを掃き出し法のような決定的なアルゴリズムでやろうとすると破綻します。
この連立方程式を解くにはGauss-Seidel法のような繰り返しで解に収束していくタイプの解き方じゃないといけません。繰り返し系の場合には、解が存在しなくても解に近い値に収束してくれます。(たまに死にます)
まぁとにかく解いたら圧力が求まります。
なんやかんやあって、

$$\Delta {\mathbf{u}}^{\prime }=-\frac{\Delta t}{\rho }\nabla P^{k+1}$$

となっているので、$P$さえわかっていれば$\Delta {\mathbf{u}}^{\prime }$もわかり、無事にタイムステップが進められます。

～完～

うまくいかねぇ

実はここまでのことを素直に実装してもうまく動きません。

というのも、基本的に流体は圧力の低い方向に流れていくのですが、流体の存在しない場所は圧力が低くなるため、粒子は互いに遠ざかるだけとなってしまいます。

なので、ある粒子が表面にいるかどうかを判定し、その部分では圧力を定数にしてやることで解決を図ります。
私はそばやであって物理屋ではないので詳しくはわかりませんが、表面張力とか気圧とかいう解釈を与えれば納得できるんですかね？

表面化どうかの判定には、「粒子数密度が一定値以下なら表面にいる」とします。
表面にいる粒子は周りにいる粒子が少なくなるのでこれで判定可能です。

これで粒子は落ち着きを取り戻します。

壁

壁の表現は基本的に粒子と同じです。
粒子と同じパラメータを持たせ、配置します。
唯一の違いは「壁は動かない」という点です。
ただし、壁の粒子は一層だと粒子が壁を抜けてしまうことがあります。
これは壁粒子の部分の粒子数密度が低いため、粒子を跳ね返すのに十分な密度を保てないことに起因します。
ですからそれを解消するために、計算にほとんど入ってこない壁を別途用意します。
これらの役目は壁付近の粒子数密度を保つということだけです。
これで壁ができました。

レンダリング

粒子をメタボールみたいなかんじで描画してやります。
法線もつけて屈折とかさせてやるとそれっぽさが滲み出てくるかもしれません。

できた

なんとかできました。
しかし安定しませんねぇ...
どうも一次方程式が解を持たないことがしばしばあるのでそこで値が爆発するようですね。
ここら辺をどう解消するのかは結局わからないままでした。。。
とりあえずラプラシアンの適用に緩和係数つけたり速度とか圧力に最大値を設けたりしてそれっぽく見せてます。
いったいどうすりゃいいんでしょうねぇ。
流体に詳しい人は是非ともご一報を。
ご指摘くださったs_cyan様に深く感謝を申し上げます。

参考文献

日本計算工学会 編 越塚誠一 著「計算力学レクチャーシリーズ5 粒子法」

明日はDekaさんの「あっ！！！PCをつけたら… 」です！お楽しみに！