LoginSignup
0
0

More than 1 year has passed since last update.

Two phase double point 形式の MPM (material point method)

Last updated at Posted at 2022-04-07

概要

 土木の分野において、粒子法の研究が盛んにおこなわれています。例えば、離散要素法(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) 粒子の物理量を更新し、グリッドをリセットする。

image.png

変形とひずみ

 変形とひずみに関して、簡単に説明する。
 時刻 $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$は、以下の通りである。

image.png

計算する際は、$(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$しかとらないので、粒子がグリッドを横切った際に符号が入れ替わり、数値振動が起こる。

 現時点では、調べきれてないので、機会があれば紹介したい。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0