概要
土木の分野において、粒子法の研究が盛んにおこなわれています。地盤工学などでよく使われる MPM は、メッシュ（グリッド）と粒子を組み合わせた計算方法であり、FEM（有限要素法）で生じるメッシュのねじれがなく、計算の途中でメッシュを使うので、SPH や MPS にくらべ計算が安定していると言われています。
MPMの欠点として、cell crossing error と呼ばれる、粒子がセルを横切った際に符号が入れ替わり数値振動が起こることが知られています。Total Lagrangian Material Point Method (TLMPM) と呼ばれる方法は、グリッド上で支配方程式を解くとき、初期に配置されていた粒子の情報を使い計算するので、粒子がセルをまたぐことがなく、cell crossing error が起こらない利点があります。前回の記事で、少し紹介しました。
通常の MPM における接触力の計算は、グリッド上で行います。しかし、TLMPM は、初期に配置されていた粒子の情報を使い計算するので、ある時刻にどこのセルに粒子がいるかは、必要な情報でないため、その情報がなくグリッド上で接触力の計算ができません。
今回の記事では、個別要素法（DEM）のような方法、直接粒子同士の接触力を計算する方法を紹介します。論文は、以下を引用しています。
物体同士の接触力
2つの粒子 $p,q$ が衝突することを考える。粒子 $p$ にける体積を $V_p$ とすれば、粒子の半径を $R_p=1/2 V_p^{1/3} $ と書くことができる。以下、下付き添え字 $p,q,\ldots$ は、粒子の番号とする。
粒子 $p,q$ が衝突すること、または、貫通することは、2つの粒子の距離 $x_{pq}$ より2つの粒子の半径の和の方が大きくなる場合で
\begin{align}
\delta &\equiv R_p+R_q - \|x_{pq}\| \geq 0 \\
x_{pq} &=x_q-x_p
\end{align}
の条件を満たす。粒子 $p$ における粒子 $q$ の接触力を $F_{pq}$ 、粒子 $q$ における粒子 $p$ の接触力を $F_{qp}$ とすると、粒子 $p,q$ における粒子の位置は、
\begin{align}
x^{t+\Delta t}_p &= x^{t}_p + \Delta t^2 \frac{F_{pq}}{m_p} \\
x^{t+\Delta t}_q &= x^{t}_q + \Delta t^2 \frac{F_{qp}}{m_q}
\end{align}
と書くことができる。$t$ は時間であり、$\Delta t$ は時間増分である。
物体は貫通しない（非貫通条件）とすれば、作用反作用の関係 $F_{qp} = -F_{pq}$ を使い、
\begin{align}
& \delta^{t+\Delta t} = 0 \\
\therefore & \ \ R_p+R_q - \|x^{t+\Delta t}_q-x^{t+\Delta t}_p\| = 0 \\
\therefore & \ \ R_p+R_q - \left\|x^{t}_q-x^{t}_p - \Delta t^2 F_{pq}\left( \frac{1}{m_p}+\frac{1}{m_q} \right)\right\| = 0 \\
\end{align}
以上より、
\begin{align}
F_{pq} &= F_{pq}^n = \frac{1}{\Delta t^2} \frac{m_p m_q}{m_p+m_q}\left(\|x_{pq}^t\|-(R_p+R_q) \right) \frac{x_{pq}^t}{\|x_{pq}^t\|}\\
\end{align}
である。$F_{pq}^n $ は、法線方向の接触力である。
次に摩擦がある場合を考える。接線ベクトル $m$ は、粒子 $p,q$ の接線方向の相対速度 $v_{pq}^{t,T}$ を使い
\begin{align}
m_{pq} &= \frac{v_{pq}^{t,T} }{\|v_{pq}^{t,T} \| } \\
v_{pq}^{t,T} &= v_{pq}^{t} - \frac{x_{pq}^t\cdot v_{pq}^t}{\|x_{pq}^t\|^2} x_{pq}^t \ \ , \ \ v_{pq}^{t} = v_q^{t}-v_p^{t}
\end{align}
とかける。接線方向の剛性係数（通常のばね係数なら次元があわない？）を $K_t$ とすると、接線方向の接触力は、
\begin{align}
F_{pq}^t = K_t v_{pq}^{t,T}
\end{align}
と書くことができる。
摩擦すべりが生じるとき、つまり、$\mu |F_{pq}^n| \leq |F_{pq}^t| $ のときは
\begin{align}
F_{pq}^t = \mu \|F_{pq}^n\| m_{pq}
\end{align}
と書くことができる。
粒子 $p$ における接触力は、粒子 $p$ に接触したすべての粒子に対し和をとり
\begin{align}
F_{p} &= \sum_{q} F_{pq} \\
F_{pq} &= F_{pq}^n + F_{pq}^t
\end{align}
と求まる。
接触力は、粒子の情報からグリッド上の力を計算する際に、外力として
\begin{align}
f_I^{ext} &= \sum_{p} N_{I}(X_p)m_p b_p^{} +\frac{1}{M_I}\sum_{p} N_{I}(X_p)m_p F_{p} \\
\end{align}
を計算する。
上記の式 $f_I^{ext}$ は、前回の記事において、
Particle to node mapping を計算する際に現れる。
物体と壁の接触力（参考にした論文と少し異なるので注意）
壁との接触も上記と同様の計算を行う。接触した壁の座標を $x_w$ とすると、接触条件（または貫通条件）は、
\begin{align}
\delta &\equiv R_p - \|x_w-x_p \| \geq 0 \\
\end{align}
である。壁は衝突によって動かないとする。壁との接触によって
\begin{align}
x^{t+\Delta t}_p &= x^{t}_p + \Delta t^2 \frac{F_{pw}}{m_p} \\
\end{align}
粒子 $p$ の位置が更新されたとする。すると、非貫通条件より
\begin{align}
& \delta^{t+\Delta t} = 0 \\
\therefore & \ \ R_p- \|x_w-x^{t+\Delta t}_p\| = 0 \\
\therefore & \ \ R_p - \left\|x_w-x^{t}_p - \Delta t^2 \frac{F_{pw}}{m_p} \right\| = 0 \\
\end{align}
したがって、壁との接触力 $F_{pw}$ は、
\begin{align}
F_{pw} &= F_{pw}^n = \frac{m_p }{\Delta t^2} \left(\|x_{pw}^t\| - R_p \right) \frac{x_{pw}^t}{\|x_{pw}^t\|}\\
x^t_{pw} &= x_w-x^t_p
\end{align}
と求まる。
次に接線方向の接触力を求める。壁に衝突した際の接線方向の相対速度は
\begin{align}
m_{pw} &= \frac{v_{pw}^{t,T} }{\|v_{pw}^{t,T} \| } \\
v_{pw}^{t,T} &= v_{pw}^{t} - \frac{x_{pw}^t\cdot v_{pw}^t}{\|x_{pw}^t\|^2} x_{pw}^t \ \ , \ \ v_{pw}^{t} = v_w-v_p^{t}
\end{align}
とかける。接線方向の剛性係数を $K_t$ とすると、接線方向の接触力は、
\begin{align}
F_{pw}^t = K_t v_{pw}^{t,T}
\end{align}
と書くことができる。
摩擦すべりが生じるとき、つまり、$\mu |F_{pw}^n| \leq |F_{pw}^t| $ のときは
\begin{align}
F_{pw}^t = \mu \|F_{pw}^n\| m_{pw}
\end{align}
と書くことができる。
最後に、粒子 $p$ における接触力に、壁との接触から得られた力を追加する。
\begin{align}
F_{p} \leftarrow F_{p} + F_{pw}^n + F_{pw}^t
\end{align}
壁との接触力を求めるには、壁の座標 $x_w$ を求める必要がある。 壁を面の方程式で表現するときは、前回の記事で紹介した方法が使える。
2つの圧縮性 Neo Hookean rings の衝突
2つの圧縮性 Neo Hookean rings の衝突問題を考える。Neo-Hookean モデルにおける、Piola-Kirchhoff 応力テンソルは
\begin{align}
P_{ij} = \mu(F - F^{-T})_{ij} + \lambda \ln(J)(F^{-T})_{ij}
\end{align}
である。$F$ は変形勾配、$J$ は変形勾配の行列式、$\mu,\lambda$ はラメ定数である。
以下の図のように、内半径 $0.03m$ 、外半径 $0.04m$ のリングを用意し、お互いに、速度 $30m/s$ で衝突させる。
物性値は、bulk modulus $K=121.7 MPa$、shear modulus $G=26.1 MPa$、密度 $1010 kg/m^3$ として計算を行った。粒子の直径は $0.00125m$ であり、1セルに4つの粒子が入るようにセルサイズを設定した。摩擦はゼロとした。
以下の動画は、シミュレーション結果である。正確に比較はできないが、論文とほぼ同じ挙動であった。
動画のヒートマップは、Von mises 応力 $q$
\begin{align}
q &= \sqrt{3 J_2(s)}\\
J_2(s)&=\frac{1}{2}s\colon s \ \ , \ \ P = \frac{1}{3}\mathrm{Tr}\sigma \ \ , \ \ s = \sigma - P I
\end{align}
を示す。
論文では、MPM と TLMPM の比較を行っている。MPM では、衝突の際に2つの物体が接触せず、ギャップが生じてしまう問題点を指摘している。
アルミニウムに高速でディスクを衝突
直径 $0.00953m$ のディスクを、速度 $1160m/s$ でアルミニウムに衝突させる。ディスクとアルミニウムの距離は、$0.01m$ 離す。以下の図の object1 がアルミニウムに対応しており、object2 がディスクに対応している。
（論文の図13 では、ディスクの重心から$0.01m$ 離しているように見えるが、そうすると図17 の結果がおかしくなる。つまり、ディスクの重心から $0.01m$ だと $0.01/1160 \approx 4\mu s$ でアルミニウムに衝突するが、論文では 約 $9\mu s$ で衝突している。）
ディスクの物性値は、ヤング率 $200 GPa$、ポアソン比 $0.3$、密度 $7850 kg/m^3$ とする。構成則は、線形弾性体を用いる。
アルミニウムの物性値は、ヤング率 $78.2 GPa$、ポアソン比 $0.3$、密度 $2700 kg/m^3$ とする。構成則は、弾塑性モデルである Von Mises モデルを用いた。
\begin{align}
\Phi(\sigma,c)&=\sqrt{3J_2(s)} -\sigma_y \\
J_2(s)&=\frac{1}{2}s\colon s \ \ , \ \ P = \frac{1}{3}\mathrm{Tr}\sigma \ \ , \ \ s = \sigma - P I
\end{align}
降伏応力 $\sigma_y = 300 MPa$ とした。
粒子の直径は $0.002m$ であり（論文は$0.001m$ ）、1セルに4つの粒子が入るようにセルサイズを設定し、摩擦はゼロとした。
以下の動画は、シミュレーション結果である。
以下の図は、ディスクがアルミニウムに食い込んだ深さを示している。正確に論文と比較はできないが、ほぼ同じ挙動をしている。
論文では、MPM と TLMPM の比較を行っている。MPM では、衝突してアルミニウムの内部に食い込む際に、アルミニウム表面に存在していた粒子がディスクに引きずられて、表面の粒子とその表面と隣同士であった粒子が完全に分かれてしまう問題点を指摘している。
最後に
固体と液体の接触が可能か試してみたい。