概要
土木の分野において、粒子法の研究が盛んにおこなわれています。例えば、離散要素法（DEM）、Smoothed Particle Hydrodynamics（SPH） 、Material Point Method（MPM）などがあります。
この記事では、土粒子と水粒子を分けた double point 形式の MPM について、簡単に説明したいと思います。
土と水を混合した粒子（2相混合）で、つまり、1つの粒子で2相を計算するMPMは、Two phase single point 形式と呼ばれる。
以下の論文の内容を参考にしています。
MPM の大まかな説明
MPM の特徴として、メッシュ（グリッド）と粒子を組み合わせた計算方法である。利点としては、FEM（有限要素法）で生じるメッシュのねじれがなく、計算の途中でメッシュを使うので、SPH や MPS にくらべ計算が安定していると言われている。
計算途中にメッシュを用いる、つまりEuler 座標を用いるが、粒子の更新は Lagrange 座標を用いるので移流項の計算は不要である。（Euler 座標とLagrange 座標を交互に使う。）
大まかな計算手順は以下の4手順である。
(a) 粒子の情報を使いグリッド点での物理量を計算する。
(b) 支配方程式を用いてグリッド点での物理量を更新する。
(c) グリッドの情報を用いて粒子の物理量を計算する。
(d) 粒子の物理量を更新し、グリッドをリセットする。
変形とひずみ
変形とひずみに関して、簡単に説明する。
時刻 $t=0$ における初期配位 $X$ 、現時刻における配位を $x$ とする。初期配位 $X$ から配位 $x$ への写像 $\phi$ とすると
\begin{align*}
&\phi : \Omega^0 \rightarrow \Omega^t \ \ \ \ \ \ \Omega^0 , \Omega^t \in \mathbb{R}^d ,X\in\Omega^0 \\
\end{align*}
Lagrange 表示では
\begin{align*}
x=x(X,t)=\phi(X,t)
\end{align*}
Euler 表示では
\begin{align*}
X=X(x,t)=\phi^{-1}(x,t)
\end{align*}
と書ける。
変形を特徴づける量、変形勾配テンソルは、
\begin{align*}
\mathbf{F}(X,t)&=\frac{\partial \phi(X,t)}{\partial \mathbf{X}}=\frac{\partial\bf{x}}{\partial \mathbf{X}} \ \ ,
F_{ij}=\frac{\partial x_i}{\partial X_j}
\end{align*}
と書くことができる。変形勾配テンソルの行列式は、体積の増分を表す。$x$ から$X$ の座標変換を考えると、体積は
\begin{align*}
\mbox{Vol} = \int \mathrm{d}x=\int \mathrm{d}X \mathrm{det}\left|\frac{\partial \bf{x}}{\partial \mathbf{X}}\right|=\int \mathrm{d}X \mathrm{det}\left|\mathbf{F} \right|=\int \mathrm{d}X \mathrm{J}
\end{align*}
と書くことができる。よって、変形勾配テンソルの行列式 $J>1$ のとき体積が増加し、$J<1$ のとき体積が減少する。
次に Lagrange 的、および Euler 的な速度と加速度を定義する。速度と加速度は、
\begin{align*}
\mathrm{V}(X,t)&=\frac{\partial \mathrm{\phi} (X,t )}{\partial t} \\ \mathrm{A}(X,t)&=\frac{\partial^2 \mathrm{\phi} (X,t )}{\partial t^2}=\frac{\partial \mathrm{V} (X,t )}{\partial t}
\end{align*}
とする。そして、Lagrange 的な速度と加速度は、
\begin{align*}
\mathrm{V}(X,t)&=\mathrm{v}\left(\phi(X,t),t \right) \\
\mathrm{A}(X,t)&=\mathrm{a}\left(\phi(X,t),t \right)
\end{align*}
Euler 的な速度と加速度は
\begin{align*}
\mathrm{v}(x,t)&=\mathrm{V}\left(\phi^{-1}(x,t),t \right) \\
\mathrm{a}(x,t)&=\mathrm{A}\left(\phi^{-1}(x,t),t \right)
\end{align*}
と定義する。Lagrange 的な速度と加速度は、速度と加速度の定義であるが、Euler 的な速度は
\begin{align*}
v_j(x,t)=V_j(\phi^{-1}(x,t),t)=\frac{\partial \phi_j (\phi^{-1}(x,t),t )}{\partial t}
\end{align*}
であり、それを微分した加速度は、
\begin{align*}
a_i(x,t)&=A_i(\phi^{-1}(x,t),t)=\frac{\partial}{\partial t} V_i (\phi^{-1}(x,t),t ) \\
&=\frac{\partial v_i(x,t)}{\partial t} +\frac{\partial v_i(x,t )}{\partial \phi_j}\frac{\partial \phi_j(\phi^{-1}(x,t),t)}{\partial t} \\
&=\frac{\partial v_i(x,t)}{\partial t} +\frac{\partial v_i(x,t )}{\partial x_j} v_j(x,t)
\end{align*}
となる。最後の行の2項目は移流項であり、グリッド上で解くEuler 的な計算方法では移流項を計算する必要があり、移流項が数値的不安定性を引き起こす場合がある。
最後にひずみについて、説明する。変形勾配テンソルを時間微分すると
\begin{align*}
\frac{\partial }{\partial t} F_{ij}(X,t )=\frac{\partial }{\partial X_j} \frac{\partial \phi_i(X,t)}{\partial t}=\frac{\partial \phi_k}{\partial X_j} \frac{\partial v_i(X,t)}{\partial \phi_k} =F_{kj}\frac{\partial v_i(X,t)}{\partial x_k}
\end{align*}
と書くことができ、速度勾配 $L$ を
\begin{align*}
\mathrm{L}=\frac{\partial \mathrm{v}}{\partial \mathrm{x}} = \frac{\partial \mathrm{F}}{\partial t} \mathrm{F}^{-1} =\dot{\mathrm{F}}\mathrm{F}^{-1}
\end{align*}
と書く。速度勾配 $L$ を対称・反対称部分に分ける。
\begin{align*}
\mathrm{L}&=\mathrm{D}+\mathrm{\Omega}\\
\mathrm{D}&=\frac{1}{2}\left(\mathrm{L}+\mathrm{L}^T \right) \ \ , \ \ \mathrm{\Omega}=\frac{1}{2}\left(\mathrm{L}-\mathrm{L}^T \right)
\end{align*}
対称部分 $D$ をひずみ速度、反対称部分 $\Omega$ をスピンテンソルと呼ぶ。ひずみテンソルとひずみ速度は、（時間微分を $\dot{}$ で表す場合もある）
\begin{align*}
\dot{\mathrm{\varepsilon}} = \mathrm{D}
\end{align*}
の関係がある。応力の時間微分は、弾性係数テンソル $C_{ijkl} $を使えば
\begin{align*}
\dot{\sigma}_{ij} = C_{ijkl} \dot{\varepsilon}_{kl}
\end{align*}
となり、積分すれば応力が求まる。しかし、この応力の微分は客観性がない。客観性とは、観測者を並進移動や回転させても、運動が変わらないことを意味する。
写像 $\phi$ を、並進変換・回転変換を行うと
\begin{align*}
\mathrm{\phi}(X,t) &\rightarrow \mathrm{\phi}^*(X,t) =\mathrm{y}(t) +\mathrm{Q}(t) \left(\mathrm{\phi}(X,t) -\mathrm{x}_0 \right)
\end{align*}
となり、変形勾配テンソルは、
\begin{align*}
\mathrm{F}(X,t) &\rightarrow \frac{\partial }{\partial \mathrm{X}} \mathrm{\phi}^*(X,t) =\mathrm{Q}(t) F(X,t)
\end{align*}
と変換される。
テンソル・ベクトル・スカラーは、並進変換・回転変換を行い
\begin{align*}
\sigma_{ij}(t) &\rightarrow Q_{ik}(t) \sigma_{kl}(t) Q^T_{lj}(t) \\
q_{i}(t) &\rightarrow Q_{ik}(t) q_{k}(t) \\
g(t) &\rightarrow g(t) \\
\end{align*}
のように変換されれば、客観性をもつ。しかし、上記の応力の微分は、
\begin{align*}
\dot{\sigma}(t) \rightarrow \dot{Q}(t) \sigma(t) Q^T(t) + Q(t) \dot{\sigma}(t) Q^T(t) + Q(t) \sigma(t) \dot{Q}^T(t)
\end{align*}
となり、客観性をもたない。
客観性をもつ応力の時間微分として、 Jaumann 速度がある。
\begin{align*}
\mathrm{\sigma}^{\bigtriangledown} =\dot{\mathrm{\sigma}} -\mathrm{\Omega }\mathrm{\sigma} +\mathrm{\sigma} \mathrm{\Omega }
\end{align*}
速度勾配を、並進変換・回転変換を行うと
\begin{align*}
\dot{\mathrm{F}}\mathrm{F}^{-1} \rightarrow \dot{\mathrm{Q}}\mathrm{Q}^T + \mathrm{Q}\dot{\mathrm{F}}\mathrm{F}^{-1} \mathrm{Q}^T
\end{align*}
さらに回転テンソルは、
\begin{align*}
\mathrm{Q}\mathrm{Q}^T = \mathrm{I} \Rightarrow \frac{\mathrm{d}} {\mathrm{d} t}(\mathrm{Q}\mathrm{Q}^T) =0 \Rightarrow \dot{\mathrm{Q}}\mathrm{Q}^T =-\mathrm{Q}\dot{\mathrm{Q}}^T
\end{align*}
となるので、Jaumann 速度を変換すると
\begin{align*}
\mathrm{\sigma}^ {\bigtriangledown} &=\dot{\mathrm{\sigma}} -\mathrm{\Omega }\mathrm{\sigma} +\mathrm{\sigma} \mathrm{\Omega } \\
& = \dot{\mathrm{\sigma}} - \frac{1}{2}\left(\mathrm{L}-\mathrm{L}^T \right) \mathrm{\sigma} + \frac{1}{2}\mathrm{\sigma}\left(\mathrm{L}-\mathrm{L}^T \right)\\
\rightarrow & \dot{Q}\sigma Q^T+ Q\dot{\sigma}Q^T + Q\sigma\dot{Q}^T \\
-& \frac{1}{2}\left(2\dot{\mathrm{Q}}\mathrm{Q}^T + \mathrm{Q}\dot{\mathrm{F}}\mathrm{F}^{-1} \mathrm{Q}^T- \mathrm{Q}\left(\dot{\mathrm{F}}\mathrm{F}^{-1}\right)^T \mathrm{Q}^T \right) Q\mathrm{\sigma}Q^T \\
+& \frac{1}{2}Q\mathrm{\sigma}Q^T \left(2\mathrm{Q}\dot{\mathrm{Q}}^T + \mathrm{Q}\dot{\mathrm{F}}\mathrm{F}^{-1} \mathrm{Q}^T- \mathrm{Q}\left(\dot{\mathrm{F}}\mathrm{F}^{-1}\right)^T \mathrm{Q}^T \right)\\
&= Q\dot{\sigma}Q^T - \frac{1}{2}\mathrm{Q}\left(L -L^T\right)\mathrm{\sigma}Q^T + \frac{1}{2}\mathrm{Q}\mathrm{\sigma}\left(\mathrm{L}-\mathrm{L}^T \right)\mathrm{Q}^T\\
&=\mathrm{Q}\mathrm{\sigma}^ {\bigtriangledown} \mathrm{Q}^T
\end{align*}
となり、客観性がある。
密度の増分に関する補足
密度の増分を考えたいので、変形勾配テンソルの増分を考える。変形勾配テンソルの増分は、微小変化$\Delta t$ を考慮して
\begin{align*}
\Delta F_{ij}=\frac{\partial \Delta x_i}{\partial X_j}= \frac{\partial}{\partial X_j}(X_i+v_i\Delta t)\approx\delta_{ij} + \frac{\partial v_i}{\partial x_j}\Delta t
\end{align*}
となる。ひずみの増分は
\begin{align*}
\Delta \varepsilon_{ij} =& \frac{\Delta t}{2} \left(\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right)
\end{align*}
であり、体積の増分は、変形勾配テンソルの行列式なので
\begin{align*}
\mathrm{det}\left|\Delta \mathbf{F} \right|\approx 1+\frac{\partial v_i}{\partial x_i}\Delta t = 1+\Delta\varepsilon_{ii}
\end{align*}
とかける。
質量一定とする。密度は体積に反比例（$\rho = m/V$）するので、時間ステップを$k$ とすると、$k+1$ステップの密度は、
\begin{align*}
\rho^{k+1} =\frac{\rho^{k}}{1+\Delta\varepsilon_{ii}} \approx \rho^{k}(1-\Delta\varepsilon_{ii})
\end{align*}
と書け、密度の増分$\Delta \rho$は
\begin{align*}
\Delta \rho = \rho^{k+1} -\rho^{k} = -\Delta\varepsilon_{ii}\rho^k
\end{align*}
と書ける。
または、一般的に質量保存の式は、
\begin{align}
\frac{\mathrm{d} \rho}{\mathrm{d} t} + \rho \frac{\partial }{\partial x_i}v_i = 0
\end{align}
と書けるので、離散化すると
\begin{align}
\Delta\rho = -\rho \Delta t \frac{\partial v_i}{\partial x_i} = -\Delta\varepsilon_{ii}\rho
\end{align}
となる。
形状関数に関する補足
以降の説明では、物理量を離散化するとき形状関数を用いる。空間座標の添え字は小文字$i$、形状関数のノード座標は大文字$I$を用いる。
形状関数 $N_{I}(x)$ を用いると
\begin{align*}
x_i(X,t)&=\sum_I N_I(X)x_{iI}(t) \\
X_i(t)&=\sum_I N_I(X)X_{iI}(t) \\
u_i(X,t)&=x_i(X,t) -X_i(t)\equiv\sum_I N_I(X)u_{iI}(t)
\end{align*}
のように離散化される。微分は
\begin{align*}
\ddot{u}_i(X,t)&= \sum_I N_I(X)\ddot{u}_{iI}(t) \\
\frac{\partial}{\partial X_j} u_{i}(X,t)&= \sum_I\frac{\partial N_{I}(X)}{\partial X_j}u_{iI}(t) \\
\end{align*}
となる。形状関数 $N_{I}(x)$の具体的な形は、6 面体要素を使い
\begin{align*}
N_{I}(\xi,\eta,\zeta) =\frac{1}{8}(1+\xi_I\xi)(1+\eta_I\eta)(1+\zeta_I\zeta)
\end{align*}
グリッド上の値 $\xi_I,\eta_I,\zeta_I$は、以下の通りである。
計算する際は、$(x,y,z)$座標を $(\xi,\eta,\zeta)$ 座標に変換する。
MPMでは、Euler 座標とLagrange 座標を交互に使う。つまり、空間方向の離散化はEuler 座標を使い、時間方向の離散化は、Lagrange 座標を使うので、以降では、初期配位 $X$ と現時刻における配位を $x$ を区別して書かない。
水粒子の支配方程式と離散化
水粒子は圧縮流体を仮定する。質量保存の式は、
\begin{align*}
\frac{\partial P_{ij}}{\partial t} = \frac{k_w}{\bar{\rho}_w}\left( \frac{\partial m_{w}}{\partial t} - \frac{\bar{\rho}_w}{n} \frac{\partial \varepsilon_V}{\partial t} \right)I_{ij}
\end{align*}
とかける。$\varepsilon_{V}$ はひずみテンソルの対角和であり、$P$ は間隙水圧で、平均化された水の密度は
\begin{align*}
\bar{\rho}_w = n\rho_w
\end{align*}
と書ける。$n$ は間隙率であり、$m_{w}$は土の間隙における水の質量である。流速は、ダルシー則を使い
\begin{align*}
{V^{\ast}_w}_i &= {V_w}_i -{V_s}_i \\
&= -\frac{K_{ij}}{\bar{\rho}_w g} \left[\frac{\partial }{\partial x_k}P_{kj} + \frac{\bar{\rho}_w}{n} \left(g_j - {a_s}_j \right) \right]
\end{align*}
と書ける。$g =g_z=g_3$ であり、${V^{\ast}_w}_i$ は相対流速、 $a_s$ は土粒子の加速度であり土粒子の計算の際に求める。
最初に、ダルシー則（相対流速）の離散化について説明する。テスト関数 $\omega_i(x)$ を使い全領域で積分すると、
\begin{align*}
\int_{\Omega} dV \bar{\rho}_w \omega_i {V^{\ast}_w}_i
&= -\frac{K_{ij}}{g}\int_{\Omega} dV \omega_i\left[\frac{\partial }{\partial x_k}P_{kj} + \frac{\bar{\rho}_w}{n} \left(g_j - {a_s}_j \right) \right] \\
&= -\frac{K_{ij}}{g}\int_{\Omega} dV \frac{\partial }{\partial x_k}\left(\omega_i P_{kj} \right) + \frac{K_{ij}}{g}\int_{\Omega} dV \frac{\partial \omega_i}{\partial x_k}P_{kj} - \frac{K_{ij}}{g}\int_{\Omega} dV \omega_i\frac{\bar{\rho}_w}{n} \left(g_j - {a_s}_j \right) \\
&= -\frac{K_{ij}}{g}\int_{\partial \Omega} dS n_k \omega_i P_{kj} + \frac{K_{ij}}{g}\int_{\Omega} dV \frac{\partial \omega_i}{\partial x_k}P_{kj} - \frac{K_{ij}}{g}\int_{\Omega} dV \omega_i\frac{\bar{\rho}_w}{n} \left(g_j - {a_s}_j \right) \\
\end{align*}
と書くことができる。
水粒子の密度は、水粒子の質量$m_{wp}$を用いて
\begin{align*}
\bar{\rho}_w(x)=\sum_{wp}^{N_{wp}} m_{wp} \delta(x-x_{wp})
\end{align*}
とかける。形状関数を使い離散化すると、
\begin{align*}
u_{i}(x_p)=\sum_{I}N_I(x_p)u_{iI}
\end{align*}
と書けるので、これらを使えば、
\begin{align*}
\int_{\Omega} dV \bar{\rho}_w \omega_i {V^{\ast}_w}_i &=
\int_{\Omega} dV \sum_{wp} m_{wp} \delta(x-x_{wp})\sum_{I,J}N_I(x)N_J(x) \omega_{iI} {V^{\ast}_w}_{iJ} \\
& = \sum_{I,J}\sum_{wp} m_{wp}N_I(x_{wp})N_J(x_{wp}) \omega_{iI} {V^{\ast}_w}_{iJ}\\
&\approx\sum_{I}\sum_{wp} m_{wp}N_I(x_{wp}) \omega_{iI} {V^{\ast}_w}_{iI}\\
\end{align*}
\begin{align*}
\frac{K_{ij}}{g}\int_{\partial \Omega} dS n_k \omega_i P_{kj} &=
h^{-1}\frac{K_{ij}}{g} \int_{\Omega} dV \frac{\bar{\rho}_w }{\bar{\rho}_w}n_k \omega_i P_{kj}(x) \\
& = h^{-1}\frac{K_{ij}}{g} \int_{\Omega} dV \sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w}\delta(x-x_{wp})\sum_{I}N_I(x)\omega_{iI} P_{jk}(x)n_k \\
& = h^{-1}\frac{K_{ij}}{g}\sum_{I} \sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w}N_I(x_{wp})\omega_{iI} P_{jk}(x_{wp})n_k
\end{align*}
\begin{align*}
\frac{K_{ij}}{g}\int_{\Omega} dV \frac{\partial \omega_i}{\partial x_k}P_{kj} &= \frac{K_{ij}}{g} \int_{\Omega} dV \sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w}\delta(x-x_{wp}) \sum_{I}\frac{\partial N_I(x)}{\partial x_k}\omega_{iI}P_{kj}(x) \\
&= \frac{K_{ij}}{g}\sum_{I}\sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w} \frac{\partial N_I(x_{wp})}{\partial x_k}\omega_{iI}P_{kj}(x_{wp}) \\
\end{align*}
\begin{align*}
\frac{K_{ij}}{g}\int_{\Omega} dV \omega_i\frac{\bar{\rho}_w}{n} \left(g_j - {a_s}_j \right)
&=\frac{K_{ij}}{g}\int_{\Omega} dV \sum_{wp}\frac{ m_{wp} }{n(x)}\delta(x-x_{wp})\sum_{I,J}N_I(x)N_J(x) \omega_{iI} \left(g_j - {a_s}_{jJ} \right) \\
&= \frac{K_{ij}}{g}\sum_{I,J}\sum_{wp}\frac{m_{wp}}{n(x_{wp})}N_I(x_{wp})N_J(x_{wp})\omega_{iI} \left(g_j - {a_s}_{jJ} \right) \\
&\approx \frac{K_{ij}}{g}\sum_{I}\sum_{wp}\frac{m_{wp}}{n(x_{wp})}N_I(x_{wp})\omega_{iI} \left(g_j - {a_s}_{jI} \right)
\end{align*}
と計算できる。$\approx $ は、lamp massを使ったこと、つまり、
\begin{align*}
\sum_{I,J}m_{wp}N_I(x_{wp})N_J(x_{wp})\approx \sum_{I}m_{wp}N_I(x_{wp})N_I(x_{wp}) = \sum_{I}m_{wp}N_I(x_{wp})
\end{align*}
として、対角成分のみで計算を行う。また$h$ は、境界の幅である。
グリッド$I$の相対流速の運動量を$P^*_{iI}$ とすると
\begin{align*}
P^*_{iI} = {f^{ext}_w}_{iI} + {f^{int}_w}_{iI}
\end{align*}
\begin{align*}
P^*_{iI} &= \sum_{wp} m_{wp}N_I(x_{wp}) {V^{\ast}_w}_{iI} \\
{f^{ext}_w}_{iI} &= - h^{-1}\frac{K_{ij}}{g} \sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w}N_I(x_{wp}) {t_w}_j - \frac{K_{ij}}{g}\sum_{wp}\frac{m_{wp}}{n(x_{wp})}N_I(x_{wp}) \left(g_j - {a_s}_{jI} \right)\\
{f^{int}_w}_{iI} &= \frac{K_{ij}}{g}\sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w} \frac{\partial N_I(x_{wp})}{\partial x_k} P_{kj}(x_{wp})
\end{align*}
と書ける。また、
\begin{align*}
{t_w}_j(x_{wp})=P_{jk}(x_{wp})n_k
\end{align*}
と置いた。グリッド$I$の質量を$M^*_{I}$ とすると
\begin{align*}
M^*_{I} &= \sum_{wp} m_{wp}N_I(x_{wp}) \\
\end{align*}
グリッド$I$の相対流速は、
\begin{align*}
{V^{\ast}_w}_{iI} = \frac{P^*_{iI}}{M^*_{I}}
\end{align*}
と書ける。
次に、質量保存の式の離散化について説明する。粒子$wp$ における、 時間間隔$\Delta t$ における間隙水圧の増分は、質量保存の式より、
\begin{align*}
\Delta P_{ij}(x_{wp}) = \frac{k_w}{\bar{\rho}_w}\left( \Delta \bar{\rho}^{\ast}_w(x_{wp}) + \frac{\bar{\rho}_w}{n(x_{wp})} \Delta \varepsilon_{sV}(x_{wp}) \right) I_{ij}
\end{align*}
と書ける。添え字$s$は、グリッド上の計算では、土の速度を使うことを意味する。
相対流速における、水粒子のひずみテンソルの増分は、
\begin{align*}
\Delta {\varepsilon_{w}}_{ij}(x_{wp}) =& \frac{\Delta t}{2} \left(\frac{\partial {V^{\ast}_w}_i}{\partial x_j} + \frac{\partial {V^{\ast}_w}_j}{\partial x_i} \right) \\
= & \frac{\Delta t}{2} \sum_{I} \left(\frac{\partial N_I(x_{wp}) }{\partial x_j}{V^{\ast}_w}_{iI} + \frac{\partial N_I(x_{wp})}{\partial x_i} {V^{\ast}_w}_{jI}\right)
\end{align*}
であり、密度の増分は
\begin{align*}
\Delta \bar{\rho}^{\ast}_w(x_{wp}) = -\Delta \varepsilon_{wV}(x_{wp})\bar{\rho}_w = -\Delta {\varepsilon_{w}}_{ii}(x_{wp})\bar{\rho}_w
\end{align*}
となる。
土粒子のひずみテンソルの増分は、
\begin{align*}
\Delta {\varepsilon_{s}}_{ij}(x_{wp}) =& \frac{\Delta t}{2} \left(\frac{\partial {V_s}_i}{\partial x_j} + \frac{\partial {V_s}_j}{\partial x_i} \right) \\
= & \frac{\Delta t}{2} \sum_{I} \left(\frac{\partial N_I(x_{wp}) }{\partial x_j}{V_s}_{iI} + \frac{\partial N_I(x_{wp})}{\partial x_i} {V_s}_{jI}\right)
\end{align*}
であり、密度の増分は
\begin{align*}
\Delta \bar{\rho}_s(x_{wp}) = -\Delta \varepsilon_{sV}(x_{wp})\bar{\rho}_w = -\Delta {\varepsilon_{s}}_{ii}(x_{wp})\bar{\rho}_w
\end{align*}
となる。
よって、間隙水圧の増分は、密度の増分を使い
\begin{align*}
\Delta P_{ij}(x_{wp}) = \frac{k_w}{\bar{\rho}_w}\left( \Delta \bar{\rho}^{\ast}_w(x_{wp}) - \frac{\Delta \bar{\rho}_s(x_{wp})}{n(x_{wp})} \right) I_{ij}
\end{align*}
と書ける。時間ステップを$k$ とすると、$k+1$ステップの間隙水圧は、
\begin{align*}
P_{ij}(x^{k+1}_{wp}) = P_{ij}(x^{k}_{wp}) + \Delta P_{ij}(x^{k}_{wp})
\end{align*}
と更新される。
土粒子の支配方程式と離散化
土粒子における支配方程式は、
\begin{align*}
\rho\frac{\mathrm{d} {V_s}_{i} }{\mathrm{d} t} = \frac{\partial }{\partial x_j}\sigma_{ij} + \rho b_i
\end{align*}
である。応力は、有効応力 $\sigma^{\prime}$と間隙水圧を使い
\begin{align*}
\sigma_{ij} = \sigma^{\prime}_{ij} +P_{ij}
\end{align*}
と書ける。密度に関しては、
\begin{align*}
\rho &= \bar{\rho}_s + \bar{\rho}_w \\
\bar{\rho}_s &= (1-n) \rho_s
\end{align*}
である。土粒子の密度は、土粒子の質量$m_{sp}$を用いて
\begin{align*}
\bar{\rho}_s(x)=\sum_{sp}^{N_{sp}} m_{sp} \delta(x-x_{sp})
\end{align*}
と書ける。
支配方程式の離散化を行う。水粒子のときと同様に、テスト関数 $\omega_i(x)$ を使い全領域で積分すると、
\begin{align*}
\int_{\Omega} dV \rho \omega_i {a_s}_{i} &=\int_{\Omega} dV \omega_i \frac{\partial \sigma_{ij}}{\partial x_j} + \int_{\Omega} dV \rho \omega_i b_i \\
& = \int_{\partial \Omega} dS \omega_i \sigma_{ij} n_j - \int_{\Omega} dV \bar{\rho}_s \frac{\partial \omega_i}{\partial x_j} \frac{\sigma_{ij}^{\prime}}{\bar{\rho}_s}
- \int_{\Omega} dV \bar{\rho}_w \frac{\partial \omega_i}{\partial x_j} \frac{P_{ij}}{\bar{\rho}_w} \\
& +\int_{\Omega} dV \bar{\rho}_w \omega_i b_i+\int_{\Omega} dV \bar{\rho}_s \omega_i b_i
\end{align*}
形状関数および密度を代入すれば、
\begin{align*}
\int_{\Omega} dV \rho \omega_i {a_s}_{i} \approx \sum_{I} \left(\sum_{wp}m_{wp}N_I(x_{wp})+\sum_{sp}m_{sp}N_I(x_{sp}) \right)\omega_{iI}{a_s}_{iI}
\end{align*}
\begin{align*}
\int_{\partial \Omega} dS \omega_i \sigma_{ij} n_j = h^{-1}\sum_{I} \sum_{sp}\frac{ m_{sp} }{\bar{\rho}_s}N_I(x_{sp})\omega_{iI} {t_s}_i(x_{sp}) + h^{-1}\sum_{I} \sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w}N_I(x_{wp})\omega_{iI} {t_w}_i(x_{wp})
\end{align*}
\begin{align*}
\int_{\Omega} dV \bar{\rho}_s \frac{\partial \omega_i}{\partial x_j} \frac{\sigma_{ij}^{\prime}}{\bar{\rho}_s} = \sum_{I}\sum_{sp}\frac{ m_{sp} }{\bar{\rho}_s} \frac{\partial N_I(x_{sp})}{\partial x_j}\omega_{iI}\sigma_{ij}^{\prime}(x_{sp}) \\
\end{align*}
\begin{align*}
\int_{\Omega} dV \bar{\rho}_w \frac{\partial \omega_i}{\partial x_j} \frac{P_{ij}}{\bar{\rho}_w} = \sum_{I}\sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w} \frac{\partial N_I(x_{wp})}{\partial x_j}\omega_{iI}P_{ij}(x_{wp}) \\
\end{align*}
\begin{align*}
\int_{\Omega} dV \bar{\rho}_w \omega_i b_i = \sum_{I}\sum_{wp} m_{wp} N_I(x_{wp}) \omega_{iI}b_i
\end{align*}
\begin{align*}
\int_{\Omega} dV \bar{\rho}_s \omega_i b_i = \sum_{I}\sum_{sp} m_{sp} N_I(x_{sp}) \omega_{iI}b_i
\end{align*}
となる。また、
\begin{align*}
{t_s}_j(x_{sp})=\sigma^{\prime}_{jk}(x_{sp})n_k
\end{align*}
と置いた。
グリッド$I$の土粒子の運動量を${P_s}_{iI}$ とすると
\begin{align*}
\dot{{P_s}}_{iI} = {f^{ext}_s}_{iI} +{f^{ext}_w}_{iI} + f^{int}_{iI}
\end{align*}
\begin{align*}
\dot{{P_s}}_{iI} &= \left(\sum_{wp}m_{wp}N_I(x_{wp})+\sum_{sp}m_{sp}N_I(x_{sp}) \right){a_s}_{iI} \\
{f^{ext}_s}_{iI} &=h^{-1} \sum_{sp}\frac{ m_{sp} }{\bar{\rho}_s}N_I(x_{sp}) {t_s}_i(x_{sp}) + \sum_{sp} m_{sp} N_I(x_{sp}) b_i \\
{f^{ext}_w}_{iI} &= h^{-1}\sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w}N_I(x_{wp}) {t_w}_i(x_{wp}) + \sum_{wp} m_{wp} N_I(x_{wp}) b_i\\
f^{int}_{iI} &= -\sum_{sp}\frac{ m_{sp} }{\bar{\rho}_s} \frac{\partial N_I(x_{sp})}{\partial x_j}\sigma_{ij}^{\prime}(x_{sp}) - \sum_{wp}\frac{ m_{wp} }{\bar{\rho}_w} \frac{\partial N_I(x_{wp})}{\partial x_j}P_{ij}(x_{wp})
\end{align*}
と書ける。
グリッド$I$の土の質量を${M_s}_{I}$ とすると
\begin{align*}
M_{I} &=\left(\sum_{wp}m_{wp}N_I(x_{wp})+\sum_{sp}m_{sp}N_I(x_{sp}) \right) \\
\end{align*}
グリッド$I$の加速度は、
\begin{align*}
{a_s}_{iI} &=\frac{\dot{{P_s}}_{iI}}{M_{I}} \\
\end{align*}
となる。
粒子の位置とグリッドの速度の更新
MPM の特徴として、Euler 座標（グリッド）とLagrange 座標（粒子）を交互に使うことである。最初にEuler 座標で支配方程式の計算を行う。時間ステップを$k$ とすると、粒子の情報を使い、各グリッドノードごとに、一時的な速度$v^L$を計算する。
\begin{align*}
{V^L_s}_{iI} &= {V^k_s}_{iI} + {a^k_s}_{iI} \Delta t \\
{V^L_w}_{iI} &= \alpha_{ij}{V^{\ast}_w}^k_{jI} + {V^L_s}_{iI}
\end{align*}
$\alpha_{ij}$は、capillary 効果（毛細管現象）を表す定数である。
グリッド上の速度を使い、次に（Lagrange 座標で）粒子の更新を行う。粒子の更新は、
\begin{align*}
{v_{sp}}_i^{k+1} &= {v_{sp}}_i^{k} + \Delta t\sum_{I} {a^k_s}_{iI}N_{I}(x^{k}_{sp})\\
{x_{sp}}_i^{k+1} &= {x_{sp}}_i^{k} + \Delta t\sum_{I} {V^L_s}_{iI}N_{I}(x^{k}_{sp})\\
{x_{wp}}_i^{k+1} &= {x_{wp}}_i^{k} + \Delta t\sum_{I} {V^L_w}_{iI}N_{I}(x^{k}_{wp})
\end{align*}
と計算される。
さらに、粒子の情報を使い、各グリッドノードごとに速度の更新を行う。
\begin{align*}
{M_s}_I &= \sum_{sp} m_{sp} N_{I}(x^{k}_{sp}) \\
{V^{k+1}_s}_{iI} &= \frac{\sum_{sp} m_{sp} {V_s}_i^{k+1} N_{I}(x^{k}_{sp})}{{M_s}_I} \\
{V^{k+1}_w}_{iI} &= \alpha_{ij}{V^{\ast}_w}^k_{jI} + {V^{k+1}_s}_{iI}
\end{align*}
形状関数の値は、$k$ステップ時の値を使う
以上のように、グリッドノードの速度の更新を行い、次に粒子の更新を行う。さらに、再度グリッドノードの速度の更新を行い、最後に応力・ひずみ・密度の更新を行う。このような計算方法を、Modified Update Stress Last (MUSL)と呼ばれる。
応力・ひずみ・密度の更新
応力・ひずみ・密度の更新について説明する。
土・水粒子のひずみの増分は、
\begin{align*}
\Delta {\varepsilon_{s}}^{k+1}_{ij}(x^{k}_{sp}) = & \frac{\Delta t}{2} \sum_{I} \left(\frac{\partial N_I(x^{k}_{sp}) }{\partial x_j}{V_s}^{k+1}_{iI} + \frac{\partial N_I(x^{k}_{sp})}{\partial x_i} {V_s}^{k+1}_{jI}\right)
\end{align*}
\begin{align*}
\Delta {\varepsilon_{w}}^{k+1}_{ij}(x^{k}_{wp}) = & \frac{\Delta t}{2} \sum_{I} \left(\frac{\partial N_I(x^{k}_{wp}) }{\partial x_j}{V_w}^{k+1}_{iI} + \frac{\partial N_I(x^{k}_{wp})}{\partial x_i} {V_w}^{k+1}_{jI}\right)
\end{align*}
と計算され、ひずみは
\begin{align*}
{\varepsilon_{s}}_{ij}(x^{k+1}_{sp}) &= {\varepsilon_{s}}_{ij}(x^{k}_{sp}) + \Delta {\varepsilon_{s}}^{k+1}_{ij}(x^{k}_{sp}) \\
{\varepsilon_{w}}_{ij}(x^{k+1}_{wp}) &= {\varepsilon_{w}}_{ij}(x^{k}_{wp}) + \Delta {\varepsilon_{w}}^{k+1}_{ij}(x^{k}_{wp}) \\
\end{align*}
と更新される。形状関数の値は、$k$ステップ時の値を使う。
有効応力は、弾性係数テンソル $C_{ijkl} $とすれば、
\begin{align*}
\sigma^{\prime}_{ij}(x^{k+1}_{sp}) &= \sigma^{\prime}_{ij}(x^{k}_{sp}) + \Delta \sigma^{\prime}_{ij}(x^{k+1}_{sp}) \\
\Delta {\sigma}^{\prime}_{ij}(x^{k+1}_{sp}) &= C_{ijkl}\Delta {\varepsilon_{s}}^{k+1}_{ij}(x^{k}_{sp}) \\
\end{align*}
と更新される。Jaumann 速度を使う場合、
\begin{align*}
\sigma^{\prime}_{ij}(x^{k+1}_{sp}) &= \sigma^{\prime}_{ij}(x^{k}_{sp}) +\Delta t
\left(\sigma^{\prime}_{il}(x^{k}_{sp}) W_{lj}(x^{k}_{sp}) -W_{il}(x^{k}_{sp}) \sigma^{\prime}_{lj}(x^{k}_{sp})\right)
+ \Delta \sigma^{\prime}_{ij}(x^{k+1}_{sp})\\
W_{ij}(x^{k}_{sp}) &= \frac{1}{2} \sum_{I} \left(\frac{\partial N_I(x^{k}_{sp}) }{\partial x_j}{V_s}^{k+1}_{iI} - \frac{\partial N_I(x^{k}_{sp})}{\partial x_i} {V_s}^{k+1}_{jI}\right)
\end{align*}
を計算する。
間隙水圧も同様に、
\begin{align*}
P_{ij}(x^{k+1}_{wp}) &= P_{ij}(x^{k}_{wp}) + \Delta P^{k+1}_{ij}(x^{k}_{wp}) \\
\Delta P^{k+1}_{ij}(x^{k}_{wp}) &= \frac{k_w}{\bar{\rho}^{k}_w}\left( \Delta \bar{\rho^{\ast}}^{k}_w(x^{k}_{wp}) - \frac{\Delta \bar{\rho}^{k+1}_s(x^{k}_{wp})}{n^{k}(x_{wp})} \right) I_{ij}
\end{align*}
と更新する。$k+1$ の上付き添え字は、グリッドノードの速度は $k+1$ ステップの値を使うことを意味する。
密度は、ひずみテンソルを使い、土・水それぞれ
\begin{align*}
\bar{\rho}^{k+1}_s &=\frac{\bar{\rho}^{k}_s}{1+\Delta {\varepsilon_{s}}^{k+1}_{ij}(x^{k}_{sp}) } \\
\bar{\rho}^{k+1}_w &=\frac{\bar{\rho}^{k}_w}{1+\Delta {\varepsilon_{w}}^{k+1}_{ij}(x^{k}_{wp}) } \\
\end{align*}
と更新する。
最後に間隙率の更新について説明する。土は非圧縮性を仮定するので、土の初期密度を$\rho_s$とすれば
\begin{align*}
\bar{\rho}^{k+1}_s = (1-n^{k+1}(x_{sp}))\rho_s
\end{align*}
が成り立ち、土の間隙率は
\begin{align*}
n^{k+1}(x_{sp}) = 1- \frac{\bar{\rho}^{k+1}_s}{\rho_s}
\end{align*}
と更新する。水の間隙率は、計算したい粒子があるセル（メッシュ）にある土の間隙率を平均として
\begin{align*}
n^{k+1}(x_{wp}) = \frac{1}{N_s} \sum_{sp=1}^{N_s}n^{k+1}(x_{sp})
\end{align*}
と計算される。$N_s$ は計算したい粒子があるセルにいる土粒子の数である。
最後に
以上の説明は、土粒子と水粒子を別々に計算する double point 形式のMPMについて説明した。ただし、以上の説明で改善点がいくつかある。
1つ目は、重力の問題である。上記で説明した方法でシミュレーションしてしまうと、初期の応力がゼロのため、重力に対する反作用が働かず、一瞬つぶれたり振動したりする。
2つ目は、形状関数に関する問題である。例えば、形状関数の $\xi$ 微分は、
\begin{align*}
\frac{\partial N_{I}(\xi,\eta,\zeta)}{\partial \xi } =\frac{1}{8}\xi_I(1+\eta_I\eta)(1+\zeta_I\zeta)
\end{align*}
となる。$\xi_I$ は、$\pm 1$しかとらないので、粒子がグリッドを横切った際に符号が入れ替わり、数値振動が起こる。
現時点では、調べきれてないので、機会があれば紹介したい。