こんにちは。流体の運動を可視化するにはNavir-Stokes方程式を数値的に解く必要がある。それに使用される手法は大きく分けて格子法と粒子法の二種類がある。これから、粒子法の基礎について書いていく。また、粒子法も大きく分けるとSPH法とMPS法に分けられ、今回はSPH法について紹介する。
目次
1.粒子法と格子法の違い
2.流体の支配方程式
3.SPH法の積分補完子とkernel関数
4.微分演算子の積分補間表現
5.WCSPH法～陽解法型非圧縮SPH法～
6.人工粘性
7.人工粘性～Balsaraスイッチ～
1. 粒子法と格子法の違い
数値計算で用いられる手法は一長一短であることが多い。粒子法と格子法もそれぞれ特徴を持っているので、扱う問題に応じて使い分ける必要がある。大まかには、格子法が工業利用、粒子法がCG業界での利用と考えられることが多い。
| 項目 | 格子法 | 粒子法 |
|---|---|---|
| 物理量の定義 | 格子の結節点 | 粒子 |
| 運動の記述法 | Eular法 | Lagrange法 |
| 計算 | 固定 | 移動 |
| 格子形成 | 必要（前処理） | 不要 |
| 近傍粒子探索 | 不要 | 必要 |
| 大変形問題 | 困難 | 容易 |
粒子法はLagrange法的な計算を行う手法である。Lagrange法は特定の粒子を追跡して追っていく手法であるのに対して、Eular法は空間に固定された点に立って場を観測する手法である。Lagrange手法は非線形項を取り扱わずに済むという利点を持っているが、非現実的な観測を行わなければならないという欠点がある。
また、粒子法は粒子が存在する点が計算点になるので、事前に格子の生成法を考えなくて済むという利点があるが、相互に影響し合う粒子を抽出する計算を毎フレーム計算しなければならないという欠点がある。
2. 流体の支配方程式
圧縮性流体の支配方程式を記す。
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \boldsymbol{u}) =0 \tag{1}
\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \cdot \nabla )\boldsymbol{u} = -\frac{1}{\rho} (\nabla p - \nabla \cdot \Pi~) + \boldsymbol{f} \tag{2}
\nu = \frac{\mu}{\rho} \tag{3}
p = \rho RT \tag{4}
R = C_p-C_v~,~~\gamma = \frac{C_p}{C_v} \tag{5}
C_s = \sqrt{\gamma \frac{p}{\rho}} = \sqrt{\gamma R T} \tag{6}
\frac{\partial p}{\partial t} = -\rho C^2_s\nabla\cdot \boldsymbol{u} - \frac{C^2_s}{\gamma} \boldsymbol{u}\cdot \nabla p \tag{7}
\nabla \cdot \Pi = \mu \Big\{ \nabla^2 \boldsymbol{u} + \frac{1}{3}\nabla (\nabla \cdot \boldsymbol{u}) \Big\} \tag{8}
非圧縮流体の場合は、以下の条件を加えることで、(1)と(8)に対して変更が加わる。
\rho = const \tag{9}
\nabla \cdot \boldsymbol{u} = 0 \tag{10}
\nabla \cdot \Pi = \mu \nabla^2 \boldsymbol{u} \tag{11}
(1)：圧縮性流体の連続体の式
(2)：ナビエストークス方程式
(3)：$\nu$は動粘性係数、$\mu$は粘性係数
(4)：熱力学方程式、$R$は気体定数
(5)：マイヤーの関係式、$C_p$は等圧比熱、$C_v$は等積比熱、$\gamma$は比熱比
(6)：$C_s$は音速
(7)：圧力の時間発展方程式
(8)：$\Pi$は粘性応力テンソル
(9)(10)(11)：$\rho$が一定だという条件(非圧縮)を用いると、(1)と(8)から(10)と(11)が導ける。
3. SPH法の積分補完子とkernel関数
粒子法は注目している粒子に対する近傍粒子からの影響量をkernel関数を用いて評価する。kernel関数$w(|\boldsymbol{r}|,h)$は上図のようにガウシアンのような形をした多項式で与える。$h$は影響範囲を表すパラメータ(smoothing length)である。中心粒子の物理量は半径2hの円内にある粒子からkernelに応じた重みをもって寄与しているという考え方から、粒子法では以下のような積分補間の式を与えている。
\phi(\boldsymbol{x}) = \int \int \int_\Omega \phi(\boldsymbol{\xi}) w(|\boldsymbol{r}|,h) dV \tag{12}
ところで、kernel関数として使える関数には以下の7つの条件がある。
➀正規化：影響域内で積分すると1になる。
\int \int \int_\Omega w(|\boldsymbol{r}|,h) dV =1 \tag{13}
➁Compact Support：影響域が有限区間で、区間の端では関数値が0になる。
w(|\boldsymbol{r}|,h) =0 ~~~~~~(|r| > 2h) \tag{14}
➂非負性：kernelは定義域内で非負。
w(|\boldsymbol{r}|,h) \geq 0 ~~~~~~(|r| \leq 2h) \tag{15}
➃単調減少性：近接している粒子ほど強く影響する。
➄デルタ関数への収束性：hを小さくすると、kernelはデルタ関数になる。
\lim_{h \to 0} w(|\boldsymbol{r}|,h) = \delta(|\boldsymbol{r}|) \tag{16}
➅偶関数：kernelは中心に対して点対称。
➆平滑性：kernelとその微分は連続かつ滑らかな関数形状。
以上の条件を満たすkernel関数はこれまで多く紹介されてきた。その中の例として、以下のkernel関数を紹介する。Monaghan and Lattanzio(1985)においては、
q \equiv \frac{r}{h} \tag{17}
3次B-Spline kernel
w_{BS}(q,h) = \alpha_{D_s}
\begin{cases}
1-\frac{3}{2}q^2 + \frac{3}{4}q^3 ~~~~(0 \leq q < 1)\\
\frac{1}{4}(2-q)^3 ~~~~~~~~~~~~~(1 \leq q < 2)\\
0~~~~~~~~~~~~~~~~~~~~~~~~~~~(2 \leq q) \tag{18}
\end{cases}
\alpha^{2D}_{D_s} = \frac{10}{7\pi h^2},~~~~~\alpha^{3D}_{D_s} = \frac{1}{\pi h^3} \tag{19}
3次B-Spline kernelは2階微分が一次関数になり、滑らかな関数ではなくなるので、ラプラシアンを含む項にはより高次のkernel関数が用いられる。
5次Q-Spline kernel
w_{QS}(q,h) = \alpha_{D_s}
\begin{cases}
(3-q)^5 -6(2-q)^5 + 15(1-q)^5 ~~~~(0 \leq q < 1)\\
(3-q)^5 -6(2-q)^5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~(1 \leq q < 2)\\
(3-q)^5~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(2 \leq q < 3) \\
0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (3 \leq q) \tag{18}
\end{cases}
\alpha^{2D}_{D_s} = \frac{7}{478\pi h^2},~~~~~\alpha^{3D}_{D_s} = \frac{1}{120\pi h^3} \tag{19}
4. 微分演算子の積分補間表現
この節では、スカラー量やベクトル量について、その勾配や発散、ラプラシアンを離散的な形で書き表す。
V_j = \frac{1}{\sum_k w(|\boldsymbol{r_k - r_j}|,h)} \tag{20}
として、スカラー量、ベクトル量はいずれも次のように離散化される。
\phi(\boldsymbol{r}_i) = \sum_j \phi(\boldsymbol{r_j}) w(|\boldsymbol{r_j - r_i}|,h) V_j \tag{21}
\boldsymbol{\phi}(\boldsymbol{r}_i) = \sum_j \boldsymbol{\phi}(\boldsymbol{r_j}) w(|\boldsymbol{r_j - r_i}|,h) V_j \tag{22}
スカラー量の勾配と発散は
\nabla_i \phi(\boldsymbol{r}_i) = \sum_{j \neq i} \phi(\boldsymbol{r_j}) \nabla w(|\boldsymbol{r_j - r_i}|,h) V_j \tag{23}
\nabla_i \cdot \boldsymbol{\phi}(\boldsymbol{r}_i) = \sum_{j \neq i} \boldsymbol{\phi}(\boldsymbol{r_j})\cdot \nabla w(|\boldsymbol{r_j - r_i}|,h) V_j \tag{24}
ここからは$w_{ij} = w(|\boldsymbol{r_j - r_i}|,h)$とする。圧力勾配項は、
\bigg \langle \frac{\nabla p}{\rho} \bigg \rangle _i = \sum_j m_j \bigg( \frac{p_j}{\rho^2_j} + \frac{p_i}{\rho^2_i}\bigg) \nabla_i w_{ij} \tag{25}
粘性項のラプラシアンには次の離散化を用いる。
\bigg \langle \nabla \cdot \bigg (\frac{\nabla p}{\rho} \bigg) \bigg \rangle _i = \sum_j \frac{8m_j}{(\rho_i + \rho_j)^2} (\phi_j - \phi_i) \frac{(\boldsymbol{r_j - r_i}) \cdot \nabla_i w_{ij}}{|\boldsymbol{r_j - r_i}|^2} \tag{26}
5. WCSPH法～陽解法型非圧縮SPH法～
現状での非圧縮SPH法は(1)～(8)が方程式系だが、温度に関する時間発展の式が無いので、閉じた問題になっていない。Monaghan(1994)で提案されたWCSPH法は(4)とは別の状態方程式(Tait型)を用いる。
\frac{p + B}{p_0 + B} = \bigg(\frac{\rho}{\rho_0} \bigg)^{\gamma_T} \tag{27}
水の場合、
\rho_0 = 1000(kg/m^3)~,~~C_{s0} = 1466(m/s)~,~~\gamma_T=7~,~~B = \frac{\rho C^2_{s0}}{\gamma_T} \tag{28}
6. 人工粘性
SPH法は本来、圧縮性流体に適応される手法で、衝撃波を産むような高速な流れのシミュレーションを行うために用いられる物である。しかし、衝撃波面で物理量が不連続になり、数値的に不安定となると言われている。そこで、人工粘性を取り入れることで、衝撃波面の厚さを計算格子幅程度にすることで物理量を連続的にし、解の発散を防ぐ。
\Pi_{ij} =
\begin{cases}
\frac{-\alpha_{\Pi} \bar{C}_{sij}\mu_{ij} + \beta_{\Pi} \mu^2_{ij}}{\bar{\rho}_{ij}} ~~~~~~~~ (\boldsymbol{u}_{ij} \cdot \boldsymbol{r}_{ij} < 0) \\
0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(\boldsymbol{u}_{ij} \cdot \boldsymbol{r}_{ij} \geq 0) \tag{29}
\end{cases}
\mu_{ij} = \frac{h \boldsymbol{u}_{ij} \cdot \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij}|^2 + \epsilon h^2} \tag{30}
\bar{C}_{sij} = \frac{C_{si} + C_{sj}}{2}~~,~~~\bar{\rho}_{ij} = \frac{\rho_{i} + \rho_{j}}{2} \tag{31}
(29)の定数に関して、おおよそ$\alpha～1$、$\beta=2\alpha$が用いられる。
7. 人工粘性～Balsaraスイッチ～
人工粘性には体積粘性と剪断粘性の２つが含まれる。ただし、衝撃波面を貫通する粒子を阻止する働きをあるのは体積粘性のみなので、剪断粘性による角運動量の損失は考慮してはならない。そのため、以下の速度勾配と回転
(\nabla \cdot \boldsymbol{u})_i = \frac{1}{\rho_i} \sum_j m_j (\boldsymbol{u}_i - \boldsymbol{u}_j) \cdot \nabla w_{ij} \tag{32}
(\nabla \times \boldsymbol{u})_i = \frac{1}{\rho_i} \sum_j m_j (\boldsymbol{u}_i - \boldsymbol{u}_j) \times \nabla w_{ij} \tag{32}
のうち、どちらが支配的であるかに応じて、人工粘性の影響度を変えなければならない。
\Pi_{ij} \rightarrow \frac{|f_{Bi} + f_{Bj}|}{2} \Pi_{ij} \tag{33}
f_{Bi} = \frac{|\nabla \cdot \boldsymbol{u}|_i}{|\nabla \cdot \boldsymbol{u}|_i + |\nabla \times \boldsymbol{u}|_i + \eta_B C_s/h_i} ~~,~~~\eta_B \approx 10^{-4} \tag{34}
f_{Bi} =
\begin{cases}
1~~~~~~(|\nabla \cdot \boldsymbol{u}|_i \ll |\nabla \times \boldsymbol{u}|_i) \\
0~~~~~~(|\nabla \cdot \boldsymbol{u}|_i \gg |\nabla \times \boldsymbol{u}|_i) \tag{35}
\end{cases}
8. 最後に
SPH法の基本的な立式はこの記事に書いたとおりである。実装にするにあたっては、時間発展や境界条件などについて考える必要があるが、それは今後、実装についての記事を書くときに回すとする。