概要
土木の分野において、粒子法の研究が盛んにおこなわれています。FEM（有限要素法）で生じるメッシュのねじれがなく、大変形の解析に適している場合があります。
地盤の変形は、基本的には塑性変形です。地盤の塑性変形において、よく使われているモデルは、Mohr Coulomb モデルや Drucker Prager モデルが挙げられます。
今回の記事では、地盤の解析でよく使われる Drucker Prager モデルについて説明したいと思います。
以下の書籍を引用しています。
論文は以下を引用しています。
論文では、Total Lagrangian 形式のSPH（Smoothed Particle Hydrodynamics）を使っていますが、この記事ではMPMを用います。
TLMPM は、前回の記事で紹介しました。
Mohr Coulomb の破壊基準の簡単な説明
Mohr Coulomb の破壊基準は、「塑性現象は、物体の粒子間に生じる摩擦すべりに起因する」という仮定から導かれる。
図のようなある物体に、上方向から $P$、横方向から $Q$ と力がかかっていたとする。通常なら $P$ の符号は負であるが、説明の都合上、正とする。$Q$ の値を増加させ、ある値で滑り出したとする。横方向における力の釣り合い式は、粘着力 （滑りを邪魔する力）$c$ を考慮して
\begin{align}
(Q - c )\cos \phi = P \sin \phi
\end{align}
と書くことができる。$Q$ の単位面積当たりの力はせん断応力 $\tau$ であり、$P$ の単位面積当たりの力は垂直応力 $\sigma$ である。
面積と$\cos\phi$で割り、まとめると Coulomb の式（暗に$\phi\neq\pi/2$として）
\begin{align}
\tau = \sigma \tan \phi + c
\end{align}
が得られる。$\phi$ は摩擦角と呼ばれる。（粘着力$c$ の単位は $\mbox{Pa}$。）
すべりが起こらない場合、つまり
\begin{align}
\tau - \sigma \tan \phi - c \leq 0
\end{align}
の場合、塑性変形は起こっていないと考えることができる。
次に、Mohrの応力円を使い主応力成分の最大値 $\sigma_{\max}$ と最小値 $\sigma_{\min}$ で Coulomb の式を表す。
Coulomb の式 が Mohrの応力円より上にある場合、
\begin{align}
\tau - \sigma \tan \phi - c < 0
\end{align}
であり塑性変形は起こらないので、Coulomb の式 が Mohrの応力円に接する場合を考える。図より、 Mohrの応力円の中心は
\begin{align}
(\sigma_c ,\tau_c) = \left(\frac{\sigma_{\max}+\sigma_{\min}}{2},0 \right)
\end{align}
である。Mohrの応力円の半径は、図から
\begin{align}
r= \frac{\sigma_{\max}+\sigma_{\min}}{2} - \sigma_{\min} = \frac{\sigma_{\max}-\sigma_{\min}}{2}
\end{align}
である。Mohrの応力円の中心から、Coulomb の式 が接する点までの$\tau$方向の距離は
\begin{align}
r \cos\phi = \frac{\sigma_{\max}-\sigma_{\min}}{2} \cos\phi
\end{align}
Mohrの応力円の中心から、Coulomb の式 が接する点までの$\sigma$方向の距離は
\begin{align}
r \sin\phi = \frac{\sigma_{\max}-\sigma_{\min}}{2} \sin\phi
\end{align}
であることがわかり、最終的にCoulomb の式 が Mohrの応力円に接する座標$(\sigma,\tau)$ は
\begin{align}
(\sigma ,\tau) = \left(\frac{\sigma_{\max}+\sigma_{\min}}{2}-\frac{\sigma_{\max}-\sigma_{\min}}{2} \sin\phi ,\frac{\sigma_{\max}-\sigma_{\min}}{2} \cos\phi \right)
\end{align}
である。この座標をCoulomb の式に代入しまとめると
\begin{align}
\frac{\sigma_{\max}-\sigma_{\min}}{2} \cos\phi - \left(\frac{\sigma_{\max}+\sigma_{\min}}{2}-\frac{\sigma_{\max}-\sigma_{\min}}{2} \sin\phi\right) \tan \phi - c = 0 \\
\frac{\sigma_{\max}-\sigma_{\min}}{2} \cos^2\phi - \left(\frac{\sigma_{\max}+\sigma_{\min}}{2}-\frac{\sigma_{\max}-\sigma_{\min}}{2} \sin\phi\right) \sin \phi - c\cos\phi =0 \\
\frac{\sigma_{\max}-\sigma_{\min}}{2}-\frac{\sigma_{\max}+\sigma_{\min}}{2} \sin \phi - c\cos\phi=0 \\
\end{align}
Mohr Coulomb の破壊基準
\begin{align}
\Phi(\sigma,c)=
(\sigma_{\max}-\sigma_{\min})-(\sigma_{\max}+\sigma_{\min}) \sin \phi - 2c\cos\phi \\
\end{align}
が得られる。通常なら $P$ の符号は負であるが、正としたため、$P$の符号を入れ替える必要がある。Coulomb の式より $\phi\to\pi-\phi$ とすればよいので、
\begin{align}
\Phi(\sigma,c)=
(\sigma_{\max}-\sigma_{\min})+(\sigma_{\max}+\sigma_{\min}) \sin \phi - 2c\cos\phi \\
\end{align}
となる。
3次元の場合、主応力成分は3個であり、Mohr Coulomb モデルは6個の関数を使い表される。Mohr Coulomb モデルの降伏関数（弾性状態か塑性状態か判断する関数）は、
\begin{align}
\Phi_1(\sigma,c)&=
(\sigma_{1}-\sigma_{3})+(\sigma_{1}+\sigma_{3}) \sin \phi - 2c\cos\phi \\
\Phi_2(\sigma,c)&=
(\sigma_{2}-\sigma_{3})+(\sigma_{2}+\sigma_{3}) \sin \phi - 2c\cos\phi \\
\Phi_3(\sigma,c)&=
(\sigma_{2}-\sigma_{1})+(\sigma_{2}+\sigma_{1}) \sin \phi - 2c\cos\phi \\
\Phi_4(\sigma,c)&=
(\sigma_{3}-\sigma_{1})+(\sigma_{3}+\sigma_{1}) \sin \phi - 2c\cos\phi \\
\Phi_5(\sigma,c)&=
(\sigma_{3}-\sigma_{2})+(\sigma_{3}+\sigma_{2}) \sin \phi - 2c\cos\phi \\
\Phi_6(\sigma,c)&=
(\sigma_{1}-\sigma_{2})+(\sigma_{1}+\sigma_{2}) \sin \phi - 2c\cos\phi \\
\end{align}
であり、主応力空間でプロットすると6角錐の形になる。
流れ則を計算する際に、塑性ポテンシャルを用いるが、降伏関数を塑性ポテンシャルとして用いる方法を関連流れ則と言い、そうでない場合を非関連流れ則と言う。
体積膨張（ダイレンタシー）現象が過大になることがあり、地盤材料においては、非関連流れ則の方が良い場合がある。それを説明するため、塑性ひずみ速度を考える。
\begin{align}
\dot{\varepsilon}^p = \dot{\gamma}N \ \ , \ \ N = \frac{\partial \Phi}{\partial \sigma}
\end{align}
$\dot{\gamma}$ は塑性乗数である。Mohr Coulomb モデルの場合（上付き添え字はベキ乗ではない）
\begin{align}
\dot{\varepsilon}^p = \sum_{i=1}^6\dot{\gamma}^iN^i \ \ , \ \ N = \frac{\partial \Phi^i}{\partial \sigma}
\end{align}
である。また
\begin{align}
\frac{\partial \sigma_i}{\partial \sigma} = e_i\otimes e_i
\end{align}
と計算でき、塑性ひずみ速度の$ e_1\otimes e_1$ 成分は、$\Phi_1,\Phi_3,\Phi_4,\Phi_6$ から
\begin{align}
\dot{\varepsilon}^p_1 = \dot{\gamma}^1(\sin\phi+1)+\dot{\gamma}^3(\sin\phi-1)+\dot{\gamma}^4(\sin\phi-1)+\dot{\gamma}^6(\sin\phi+1)
\end{align}
、塑性ひずみ速度の$ e_2\otimes e_2$ 成分は、$\Phi_2,\Phi_3,\Phi_5,\Phi_6$ から
\begin{align}
\dot{\varepsilon}^p_2 = \dot{\gamma}^2(\sin\phi+1)+\dot{\gamma}^3(\sin\phi+1)+\dot{\gamma}^5(\sin\phi-1)+\dot{\gamma}^6(\sin\phi-1)
\end{align}
、塑性ひずみ速度の$ e_3\otimes e_3$ 成分は、$\Phi_1,\Phi_2,\Phi_4,\Phi_5$ から
\begin{align}
\dot{\varepsilon}^p_3 = \dot{\gamma}^1(\sin\phi-1)+\dot{\gamma}^2(\sin\phi-1)+\dot{\gamma}^4(\sin\phi+1)+\dot{\gamma}^5(\sin\phi+1)
\end{align}
と計算される。塑性ひずみ速度の体積成分は、 $\dot{\gamma}^i\geq0$ 、$0\leq\phi<\pi/2$ なので、
\begin{align}
\dot{\varepsilon}^p_V &= \dot{\varepsilon}^p_1+\dot{\varepsilon}^p_2+\dot{\varepsilon}^p_3 \\
&= 2\sin\phi (\dot{\gamma}^1+\dot{\gamma}^2+\dot{\gamma}^3+\dot{\gamma}^4+\dot{\gamma}^5+\dot{\gamma}^6)\geq0
\end{align}
となる。よって、塑性ひずみ速度の体積成分は膨張を表す。これは地盤変形において、塑性変形中に体積膨張（ダイレンタシー）現象が起きることを意味する。ただし、数値計算においてダイレンタシー現象が過大になることがある。その場合、摩擦角 $\phi$ より低角度のダイレンタシー角 $\psi$ に置き換えた塑性ポテンシャルを用いて計算される。
Drucker Prager モデル
Drucker Prager モデルの降伏関数は、Mohr Coulomb モデルをなめらかに近似したモデルであり
\begin{align}
\Phi(\sigma,c)&=\sqrt{J_2(s)} + \eta P -\xi c(\bar{\varepsilon}^p) \\
J_2(s)&=\frac{1}{2}s\colon s \ \ , \ \ P = \frac{1}{3}\mathrm{Tr}\sigma \ \ , \ \ s = \sigma - P I
\end{align}
である。$\sigma$ は応力、$P$ は静水圧、$s$ は偏差ひずみである。主応力空間でプロットすると円錐の形になる。
また、累積塑性ひずみは
\begin{align}
\dot{\bar{\varepsilon}^p} \equiv \dot{\gamma}\xi
\end{align}
と定義される。
パラメータ $\eta,\xi$ はMohr Coulomb モデルに応じて決められ、6角錐の外側の稜線を通る円と一致させる場合、
\begin{align}
\eta \equiv \frac{6\sin\phi}{\sqrt{3} (3-\sin\phi)} \ \ , \ \
\xi \equiv \frac{6\cos\phi}{\sqrt{3} (3-\sin\phi)}
\end{align}
6角錐の内側の稜線を通る円と一致させる場合、
\begin{align}
\eta \equiv \frac{6\sin\phi}{\sqrt{3} (3+\sin\phi)} \ \ , \ \
\xi \equiv \frac{6\cos\phi}{\sqrt{3} (3+\sin\phi)}
\end{align}
とする。
塑性ポテンシャルは、
\begin{align}
\Psi(\sigma,c)&=\sqrt{J_2(s)} + \bar{\eta} P
\end{align}
を使い、パラメータ $\bar{\eta}$ は、摩擦角 $\phi$ より低角度のダイレンタシー角 $\psi$ を使い、6角錐の外側の稜線を通る円と一致させる場合、
\begin{align}
\bar{\eta} \equiv \frac{6\sin\psi}{\sqrt{3} (3-\sin\psi)}
\end{align}
6角錐の内側の稜線を通る円と一致させる場合、
\begin{align}
\bar{\eta} \equiv \frac{6\sin\psi}{\sqrt{3} (3+\sin\psi)}
\end{align}
とする。
1. Return Mapping アルゴリズム
以下、応力の計算は、return mapping を用いて計算を行う。return mapping アルゴリズムの詳細は省略する。塑性変形は、降伏関数がゼロを超えれば発生すると仮定している。return mapping アルゴリズムは、図のように降伏関数をゼロに戻すように、つまり、弾性領域（図の赤い領域に囲まれた部分）に応力を修正するアルゴリズムである。$\sqrt{J_2(s)}$ は、ゼロ以上の値を必ず取るが、return mapping を行うと負の値 $(\sqrt{J_2(s)}<0)$ となる場合がある。その場合は、頂点に値を戻す必要がある。
（最初に登場したMohrの円とは、横軸の符号が逆なので図が反転している。）
return mapping を使い応力を計算したいので、流れ則
\begin{align}
\dot{\varepsilon}^p = \dot{\gamma}N \ \ , \ \ N = \frac{\partial \Phi}{\partial \sigma}
\end{align}
を具体的に計算する。離散化した式は
\begin{align}
\Delta \varepsilon^p = \Delta \gamma N
\end{align}
である。
まずは、静水圧と偏差応力を応力で微分すると
\begin{align}
\frac{\partial P}{\partial \sigma_{ij}} &= \frac{1}{3} \frac{\partial }{\partial \sigma_{ij}} \sigma_{mm} = \frac{1}{3} \delta_{ij} \\
\frac{\partial s_{kl}}{\partial \sigma_{ij}} & = \frac{\partial }{\partial \sigma_{ij}} \left(\sigma_{kl} - \frac{1}{3} \delta_{kl} \sigma_{mm} \right) = \delta_{ik}\delta_{jl} - \frac{1}{3} \delta_{kl}\delta_{ij}
\end{align}
であり、偏差応力はトレースレスであるので、塑性ポテンシャルの微分は
\begin{align}
\frac{\partial \Psi}{\partial \sigma_{ij}} &= \frac{\partial }{\partial \sigma_{ij}} \left(\sqrt{J_2(s)} + \bar{\eta} P \right) \\
&=\frac{1}{2\sqrt{J_2(s)}}\left(\delta_{ik}\delta_{jl} - \frac{1}{3} \delta_{kl}\delta_{ij} \right) s_{kl} + \frac{\bar{\eta}}{3}\delta_{ij} \\
&= \frac{1}{2\sqrt{J_2(s)}} s_{ij} + \frac{\bar{\eta}}{3}\delta_{ij}
\end{align}
と求まる。
応力の更新は、return mapping の更新式
\begin{align}
\sigma_{ij} &= \sigma_{ij}^{\ast} - \Delta\gamma D^e_{ijkl} N_{kl} \\
D^e_{ijkl} &= 2G I^{div}_{ijkl} + K\delta_{ij}\delta_{kl} \ \ , \ \ I^{div}_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})-\frac{1}{3}\delta_{ij}\delta_{kl}
\end{align}
を使う。$K$ は体積弾性係数、$G$ はせん断弾性係数である。$\sigma^{\ast} $は試行状態を表し、これは弾性状態として計算した応力である。試行状態の偏差応力と静水圧は、$s^{\ast},P^{\ast}$と書く。
（有限変形理論では、単純にひずみを弾性成分と塑性成分に分解できない。ただし、hencky 理論のように、対数ひずみを使えば、微小変形理論のようひずみを弾性成分と塑性成分に分解される。）
FEM などで準静的な解析を行う場合、ニュートン法で変位を更新するとき応力をひずみで微分した係数、接線係数を求める必要がある。今回の記事では、陽解法で解くので接線係数は不要である。
2. 頂点以外の Return Mapping アルゴリズム
$\sqrt{J_2(s)}$ が、ゼロ以上の値をとる場合を説明する。Drucker Prager モデルをreturn mapping の更新式に代入すると
\begin{align}
\sigma_{ij} &= \sigma_{ij}^{\ast} - \Delta\gamma D^e_{ijkl} \frac{\partial \Psi}{\partial \sigma_{kl}} \\
& = \sigma_{ij}^{\ast} - \Delta\gamma \left(\frac{G}{\sqrt{J_2(s)}} s_{ij} + \bar{\eta} K\delta_{ij} \right)
\end{align}
と求まる。
上式では偏差応力と塑性乗数が未知である。偏差応力を計算してみると
\begin{align}
s_{ij} &= \sigma_{ij} -\frac{\delta_{ij}}{3} \sigma_{mm} \\
&=\sigma_{ij}^{\ast} - \Delta\gamma \left(\frac{G}{\sqrt{J_2(s)}} s_{ij} + \bar{\eta} K\delta_{ij} \right)
- \frac{\delta_{ij}}{3}\left(\sigma_{mm}^{\ast} - 3\Delta\gamma
\bar{\eta}K \right) \\
&= s_{ij}^{\ast} - \frac{\Delta\gamma G}{\sqrt{J_2(s)}} s_{ij}\\
\\
\therefore \ \ & \left(1+ \frac{\Delta\gamma G}{\sqrt{J_2(s)}} \right)s_{ij} = s_{ij}^{\ast}
\end{align}
となる。これは偏差応力が、試行状態の偏差応力に比例することを意味するので
\begin{align}
J_2(s^{\ast}) &=\frac{1}{2} s_{ij}^{\ast}s_{ij}^{\ast} = \frac{1}{2}\left(1+ \frac{\Delta\gamma G}{\sqrt{J_2(s)}} \right)^2s_{ij}s_{ij} \\
&=\left(1+ \frac{\Delta\gamma G}{\sqrt{J_2(s)}} \right)^2J_2(s) \\
\\
\therefore \ \ \frac{s_{ij}^{\ast}}{\sqrt{J_2(s^{\ast})}} &= \frac{\left(1+ \frac{\Delta\gamma G}{\sqrt{J_2(s)}} \right)s_{ij}}{\sqrt{\left(1+ \frac{\Delta\gamma G}{\sqrt{J_2(s)}} \right)^2J_2(s) }} =\frac{s_{ij}}{\sqrt{J_2(s)}}
\end{align}
応力は
\begin{align}
\sigma_{ij} &= \sigma_{ij}^{\ast} - \Delta\gamma \left(\frac{G}{\sqrt{J_2(s^{\ast})}} s_{ij}^{\ast} + \bar{\eta} K\delta_{ij} \right)
\end{align}
と更新される。よって、静水圧と偏差応力は
\begin{align}
P &= P^{\ast} - K\bar{\eta}\Delta \gamma \\
s_{ij} &= \left(1 - \frac{\Delta\gamma G}{\sqrt{J_2(s^{\ast} )}} \right)s_{ij}^{\ast} \\
\sqrt{J_2(s)} &= \left(1 - \frac{\Delta\gamma G}{\sqrt{J_2(s^{\ast} )}} \right)\sqrt{J_2(s^{\ast} )} = \sqrt{J_2(s^{\ast} )} - \Delta\gamma G
\end{align}
と更新する。累積塑性ひずみは
\begin{align}
\dot{\bar{\varepsilon}^p} \equiv \dot{\gamma}\xi
\end{align}
と定義されるので、離散化すると
\begin{align}
\Delta \bar{\varepsilon}^p = \Delta \gamma\xi
\end{align}
と書ける。累積塑性ひずみの更新は
\begin{align}
\bar{\varepsilon}^p \leftarrow \bar{\varepsilon}^p + \Delta \bar{\varepsilon}^p
= \bar{\varepsilon}^p + \Delta \gamma\xi
\end{align}
と更新する。
上記の計算において、塑性乗数が未知である。塑性変形は、降伏関数がゼロを超えれば発生すると仮定している。return mapping アルゴリズムは、降伏関数をゼロに戻すように応力を修正するので
\begin{align}
\Phi(\sigma,c)&=\sqrt{J_2(s)} + \eta P -\xi c(\bar{\varepsilon}^p) =0 \\
\tilde{\Phi}(\Delta\gamma)&=\sqrt{J_2(s^{\ast} )} - \Delta\gamma G + \eta \left( P^{\ast} - K\bar{\eta}\Delta \gamma \right)-\xi c(\bar{\varepsilon}^p)=0
\end{align}
を$\Delta \gamma$ について解けばよい。しかし、粘着力は、一般に、累積塑性ひずみの関数である。よって、ニュートン法などで繰り返し計算を行う必要がある。今回の記事では、粘着力は以下の線形関係、線形硬化モデルを使う。
\begin{align}
c(\bar{\varepsilon}^p) = c+H\bar{\varepsilon}^p(\Delta\gamma)
\end{align}
$H$は硬化係数である。硬化係数の項は、塑性ひずみが蓄積されると粘着力が強くなる（硬化）ことを意味し、単純に塑性変形が起こりにくくなる。
よって、$\Delta \gamma$ は
\begin{align}
&\sqrt{J_2(s^{\ast} )} - \Delta\gamma G + \eta \left( P^{\ast} - K\bar{\eta}\Delta \gamma \right)-\xi \{ c+H \bar{\varepsilon}^p + H\Delta \gamma\xi \}=0 \\
\\
\therefore \ \ & \Delta \gamma = \frac{\sqrt{J_2(s^{\ast} )}+ \eta P^{\ast}-\xi( c+H \bar{\varepsilon}^p) }{G+\eta\bar{\eta}K+\xi^2H}
\end{align}
と求まる。
3. 頂点の Return Mapping アルゴリズム
上記で計算した $\sqrt{J_2(s)}$ の値が負になる場合、つまり、
\begin{align}
\sqrt{J_2(s)} & = \sqrt{J_2(s^{\ast} )} - \Delta\gamma G < 0
\end{align}
となる場合、図のような頂点への return mapping が必要となる。このとき、偏差応力は頂点に引き戻されるので
\begin{align}
s_{ij} &= 0
\end{align}
となる。応力の更新式は
\begin{align}
\sigma_{ij} &= \sigma_{ij}^{\ast} - \Delta\gamma D^e_{ijkl} \frac{\partial \Psi}{\partial \sigma_{kl}} \\
& = \sigma_{ij}^{\ast} - \Delta\gamma \bar{\eta} K\delta_{ij} \\
& \equiv \sigma_{ij}^{\ast} - \Delta \varepsilon_V^p K\delta_{ij} \\
\end{align}
と書ける。$\Delta \varepsilon_V^p $ は、頂点の return mapping アルゴリズムにおける塑性乗数である。
累積塑性モデルは、体積成分しかないので
\begin{align}
\bar{\varepsilon}^p \leftarrow \bar{\varepsilon}^p + \frac{\xi}{\eta} \Delta \varepsilon_V^p
\end{align}
と更新する。降伏関数を、
\begin{align}
\tilde{\Phi}(\Delta \varepsilon_V^p )&= P^{\ast} - K\Delta \varepsilon_V^p -\frac{\xi}{\bar{\eta}} c(\bar{\varepsilon}^p)=0
\end{align}
と変更する。線形硬化モデルを使えば、$\Delta \varepsilon_V^p$ は
\begin{align}
& P^{\ast} - K\Delta \varepsilon_V^p -\frac{\xi}{\bar{\eta}} \left\{ c+H \bar{\varepsilon}^p + H\frac{\xi}{\eta} \Delta \varepsilon_V^p \right\}=0 \\
\\
\therefore & \ \ \Delta \varepsilon_V^p = \frac{ P^{\ast}-\xi( c+H \bar{\varepsilon}^p)\frac{\xi}{\bar{\eta}} }{K+\frac{\xi}{\eta}\frac{\xi}{\bar{\eta}}H}
\end{align}
と求まる。最終的に、静水圧と偏差応力は
\begin{align}
P &= P^{\ast} - K\Delta \varepsilon_V^p \\
s_{ij} &= 0
\end{align}
と更新する。
ひずみ・応力の求め方
上記の return mapping アルゴリズムを適用するためには、試行状態の応力を求める必要がある。変形勾配 $F$ を計算し、変形速度勾配
\begin{align}
L = \dot{F}F^{-1}
\end{align}
を求める。変形速度勾配の対称部分は、ひずみ速度と呼ばれ
\begin{align}
\dot{\varepsilon} = \frac{1}{2}\left(L+L^T\right)
\end{align}
と計算される。反対称成分は、スピンテンソルと呼ばれ
\begin{align}
\omega = \frac{1}{2}\left(L-L^T\right)
\end{align}
と計算される。それぞれ、離散化を行うと
\begin{align}
\Delta \varepsilon = \Delta t \dot{\varepsilon} \ \ , \ \ \Delta \omega = \Delta t \omega
\end{align}
となる。ラメ定数 $\lambda,\mu$ を使い、応力は
\begin{align}
\sigma_{ij} \leftarrow \sigma_{ij} + \Delta \sigma_{ij} = \sigma_{ij} +\lambda \delta_{ij}\Delta \varepsilon_{mm} +2\mu \Delta \varepsilon_{ij}
\end{align}
と更新される。 Jaumman 速度を使う場合、
\begin{align}
\mathrm{\sigma}_{ij}^{\bigtriangledown} =\dot{\mathrm{\sigma}}_{ij} -\omega _{im}\mathrm{\sigma}_{mj} +\mathrm{\sigma}_{im} \omega _{mj}
\end{align}
離散化すると応力は、
\begin{align}
\sigma_{ij} \leftarrow \sigma_{ij} +\lambda \delta_{ij}\Delta \varepsilon_{mm} +2\mu \Delta \varepsilon_{ij}+ \Delta \omega_{im}\mathrm{\sigma}_{mj} -\sigma_{im} \Delta \omega_{mj}
\end{align}
と更新される。
今回の記事では、Total Lagrangian 形式で計算を行うので、上記の return mapping アルゴリズムを適用した後、 Piola-Kirchhoff 応力
\begin{align}
P = J \sigma F^{-T} \ \ , \ \ J = \mathrm{det}F
\end{align}
を計算する。
（引用した論文では、恐らく、上記の方法で Piola-Kirchhoff 応力を計算している。個人的にはこのような方法で有限変形の塑性問題の計算をしてよいか悩む。）
地すべりのシミュレーション
TLMPM を使い、長さ $20\mbox{m}$、高さ $7\mbox{m}$ の 盛土の地すべりシミュレーションを行った。引用した論文では、2次元の解析（平面ひずみ）を行っているが、この記事では、3次元の解析で奥行き方向には周期境界条件を課している。（奥行き方向に粒子が移動しないようにDirichet 境界条件を課した。）
格子間隔は、$0.1\mbox{m}$ とし１つのセルに4つの粒子が入るようにした。 密度は $1850 \mbox{kg}/\mbox{m}^3$、ヤング率 $100 \times 10^{6} \mbox{pa}$、 ポアソン比 $ 0.3$、重力 $9.8 \mbox{m}/\mbox{s}^2$ とした。
Drucker Prager モデルのパラメータは、摩擦角 $\pi/6$ 、ダイレンタシー角 $\pi/6$ 、粘着力 $5.0\times 10^{3} \mbox{pa}$ とした。論文では、安全率 $f_s$ で摩擦角・粘着力を割っているので、$f_s=2.0$ とした。
シミュレーション結果は、以下の通りである。
以下のヒートマップは、$x$ 方向の変位を表している。
以下のヒートマップは、累積塑性ひずみを表している。
論文の結果と異なるが、SPH と MPM の違いもあり、また平面ひずみの問題を3次元で境界条件を課して解いているからだろう。（はっきりとした理由は検討するかもしれない。）
最後に
TLMPM を使い弾塑性解析を行った。次は、土水連成計算を行いたい。