【CFD/粒子法】半陰解法による粒子移動の計算【粒子法入門 #6 MPS法編】

粒子法を用いた流体解析（CFD）の解説を行っています。第６回は、半陰解法を用いた粒子移動の計算について解説します。

前回は粒子数密度と勾配・発散・ラプラシアンについて解説しました。まだ見てないかたは下記リンクからどうぞ。ざっくり説明すると、MPS法では粒子数密度という値を用いて、微分演算子を計算します。計算結果が少し特殊なので、計算式よりはイメージで覚えておくと良いですよーという話をしました。

粒子数密度と勾配・発散・ラプラシアン　【粒子法入門　#5　MPS法編】

陽解法と陰解法

今回出てくる重要な単語が半陰解法です。これは陽解法と陰解法を両方知っておかないと理解できないので、まずは陽解法と陰解法についてそれぞれ簡単に説明します。

陽解法

陽解法とは、現在の結果を使って未来を予想することです。何も記述されていない式はほとんど陽解法で解きます。まずはナビエ・ストークス方程式を見てください。

$$\frac {D v}{D t}=\frac {1}{\rho} \nabla p + { \nu {\nabla}^2 v}+{g}$$

このままだとPC上で計算できないので、式の左辺を時間で離散化します。

$$\frac {D v}{D t} = \frac{v^{n+1} – v^n}{\Delta t}$$

これは、速度に対して時間差分を取っています。分子は速度vの未来と現在の値の差分で、分母は時間の差分である時間ステップを使用しています。nは現在、n+1は未来の値を表しています。

このように離散化することで、ナビエ・ストークス方程式を用いて未来の値を計算することができるようになります。ここで、左辺の項は全て現在の時間の値を使用することとします。

$$\frac{v^{n+1} – v^n}{\Delta t} =\frac {1}{\rho} {\nabla p}^n + { \nu {{\nabla}^2 v}^n}+{g} $$

左辺にはn+1、右辺にはnの項が集まるように整理します。

$$v^{n+1} =v^n + \Delta t (\frac {1}{\rho} {\nabla p}^n + { \nu {{\nabla}^2 v}^n}+{g}) $$

これで、未来の値である$v^{n+1}$を計算できる式ができました。圧力項や粘性項の計算に工夫は必要ですが、次の時間ステップの速度が計算できるイメージは持てるでしょう。

陰解法

陽解法では、現在の値から未来の値を予測しました。一方で陰解法では、未来の値から未来の値を予測します。そんなことできるのか疑わしいと思いますが、実は数値計算で工夫をすれば結構良い結果が出ます。CFDではむしろ陰解法のほうがメジャーだったりします。

先程のナビエ・ストークス方程式に対して、陰解法を適用します。右辺の文字がnからn+1に変わっただけです。

$$\frac{v^{n+1} – v^n}{\Delta t} =\frac {1}{\rho} {\nabla p}^{n+1} + { \nu {{\nabla}^2 v}^{n+1}}+{g} $$

左辺にn+1、右辺にnで整理すると下記のようになります。

$$v^{n+1} + \Delta t (\frac {1}{\rho} {\nabla p}^{n+1} + { \nu {{\nabla}^2 v}^{n+1}) = v^n + \Delta t {g} $$

これだと左辺にたくさん項があり、$v^{n+1}$をまともに計算できません。よって、陰解法では行列計算を用いた収束計算で対応します。収束計算のイメージとしては、両辺を少しずつ修正していって等号が成り立つ条件を探します。地道なので大量の計算が必要ですが、収束性が高い手法もたくさん確立されているので結構良い結果が出ます。

今回は陰解法については細かい話はしませんので、これから話す半陰解法からイメージを膨らましてください。

半陰解法

半陰解法では、粘性項は陽的、圧力項は陰的に解きます。よって、ナビエ・ストークス方程式は下記のようになります。

$$\frac{v^{n+1} – v^n}{\Delta t} =\frac {1}{\rho} {\nabla p}^{n+1} + { \nu {{\nabla}^2 v}^{n}}+{g} $$

MPS法では、圧力とその他の項（粘性・重力項）を分けて計算します。このような解き方を分離解法、フラクショナルステップ法、圧力分離型ソルバーと言ったりします。色々呼び方がありますが、とにかく分けて計算するというニュアンスが伝わればOKです。

$$v^{*} =v^n + \Delta t ({ \nu {{\nabla}^2 v}^n}+{g}) $$

ここで、粘性項と重力項によって更新した速度を仮の速度$u^{*}$とします。これを圧力計算に使用します。${{\nabla}^2 v}$は前回紹介したラプラシアンの式を使用して計算します。

$$ {\nabla}^2 v = \frac{2d}{\lambda n^0} \sum_{j \neq i} (v^n_j-v^n_i)w(|r_j-r_i|) $$

$$ \lambda = \frac{1}{n^0} \sum_{j \neq i} |r_j-r_i|^2 w(|r_j-r_i|) $$

詳細は前回の記事を確認してください。

仮の速度が求まったので、次に仮の位置を求めます。

$$r^{*} = r^n + \Delta t v^{*}$$

仮の速度$v^{*}$で時間ステップ$\Delta t$分だけ進んだのが仮の位置になります。

次は圧力を求めます。更新後の粒子位置から粒子密度を計算でき、粒子数密度から圧力が得られます。下記のポアソン方程式から圧力を求めます。

$${\nabla}^2 p^{n+1} = -\rho \frac{1}{{\Delta t}^2} \frac{n^* – n^0}{n^0}$$

ポアソン方程式は連続の式$\nabla \cdot v = 0$とナビエ・ストークス方程式から得られます。詳細は書籍を参照ください。

右辺にある$n^{*}$が粒子の位置情報から粒子数密度です。注意として、この式は右辺を計算するだけではありません。最終的にはPを求める必要があるので、左辺にラプラシアンモデル（第５回参照）を適用して、行列の収束計算をする必要があります。

最後に圧力の計算結果をナビエ・ストークス方程式に代入します。

$$v^{n+1} =v^* + \Delta t (\frac {1}{\rho} {\nabla p}^{n+1})$$

次のステップの速度が求まったので、位置も更新しておきましょう。

$$r^{n+1} = r^* + \Delta t (v^{n+1} – v^*)$$

ここまでが半陰解法の一連の流れです。お疲れさまでした。

計算手順まとめ

１．ナビエ・ストークス方程式を使って、粘性項と重力項から仮の速度を求める
２．仮の速度で仮の位置を求める
３．仮の位置から粒子数密度を求める（数式は第５回参照）
４．ポアソン方程式を使って、粒子数密度から圧力を求める
５．ナビエ・ストークス方程式を使って、仮の速度と圧力項から次のステップの速度を求める
６．次のステップの速度から次のステップの位置を求める
７．１に戻る

上記のループを回していくことで、粒子の移動を計算することができます。特に時間がかかるのは圧力のポアソン方程式であり、この部分でいかに工夫するかが計算屋の腕の見せ所になります。