概要

　フェーズフィールド法を用いた亀裂進展解析は、XFEMに比べ幾何学的な処理がなく、ペリダイナミックスに比べ境界部分の処理が容易であることが知られています。前回の記事で紹介しました。

　破壊には大きく分けて２つあり、脆性破壊と延性破壊がある。延性破壊の場合、塑性変形が起こり破壊に至ります。今回の記事では、延性破壊におけるフェーズフィールド法について紹介したいと思います。

以下の文献を参考にしています。

最初に

　論文ではフェーズフィールド変数を

\begin{align}
s(x) =  \begin{cases}
    0 & \text{完全破壊} \\
    1 & \text{非破壊} \\
  \end{cases}
\end{align}

と定義しているが、今回の記事では

\begin{align}
\phi(x) =  \begin{cases}
    1 & \text{完全破壊} \\
    0 & \text{非破壊} \\
  \end{cases}
\end{align}

とする。論文のパラメータ$l$は、

\begin{align}
l_0 = \frac{l}{2}
\end{align}

と対応する。

　記事では以下の関数を用いる。

\begin{align}
 \left\langle x \right\rangle _{\pm} =\frac{1}{2}\left(x\pm|x|\right)
\end{align}

\begin{align}
H(x) = 
\begin{cases} 
0 & \text{if } x < 0 \\ 
1 & \text{if } x \geq 0 
\end{cases}
\end{align}

また、$\lambda$と$\mu$はラメ定数で

\begin{align}
K = \lambda + \frac{2}{3} \mu
\end{align}

は体積弾性率である。

Von Mises モデルに基づくフェーズフィールドモデル

　全エネルギー$\Psi$は以下となる。

\begin{align}
\Psi = \int_{\Omega} \mathrm{d}\Omega \{g(\phi)\psi_e^+ + \psi_e^- + \psi_p\}+ \int_{\Gamma} \mathrm{d}S G_c 
    - \int_{\Omega} \mathrm{d}\Omega {\bf b} \cdot {\bf u} - \int_{\partial\Omega} \mathrm{d}S  {\bf f} \cdot {\bf u} 
\end{align}

各項の意味は

第1項：弾性エネルギー・塑性エネルギー
第2項：破壊エネルギー
第3項：物体力による仕事
第4項：表面力による仕事

である。弾性エネルギーの分割は、Amor et al. モデル

\begin{align*}
\psi_e^{+} &= \frac{1}{2}K \left\langle \mathrm{Tr}(\varepsilon) \right\rangle _{+} ^2 + \mu \varepsilon_{div} : \varepsilon_{div}
\\
\psi_e^{-} &= \frac{1}{2}K \left\langle \mathrm{Tr}(\varepsilon) \right\rangle_{-} ^2 
\end{align*}

塑性エネルギーは、Von Mises モデルを用いる。

\begin{align*}
\psi_p &= \sigma_Y \alpha + \frac{1}{2}h\alpha^2
\end{align*}

$\sigma_Y $は降伏応力、$h$は硬化係数、$\alpha $は硬化変数である。破壊エネルギーは以下である。

\begin{align*}
\int_{\Gamma} \mathrm{d}S G_c = \int_{\Omega} \mathrm{d}\Omega \frac{G_c}{2l_0} \left\{ \phi^2 + l_0^2 \nabla\phi \cdot  \nabla\phi \right\}
\end{align*}

　延性破壊は、塑性変形を伴いながら亀裂が生じる必要がある。そのため、適切な変位とフェーズフィールド変数の結合関数$g(\phi)$を設定する必要がある。この論文では、結合関数を

\begin{align*}
g(\phi) = (1-\phi)^{2 q} + k
\end{align*}

と設定する。$k$は数値的不安定を回避するためのパラメータである。$q$は累積塑性ひずみ$\varepsilon_{ep}^p $

\begin{align*}
\varepsilon_{ep}^p = \sqrt{\frac{2}{3}} \int_0^t \mathrm{d}\tau \sqrt{\dot{{\bf \varepsilon}}^p :\dot{{\bf \varepsilon}}^p }
\end{align*}

を用いて定義される変数である。

\begin{align*}
q = \frac{\varepsilon_{ep}^p}{\varepsilon_{ep,crit}^p}
\end{align*}

$\varepsilon_{ep,crit}^p $はしきい値である。結合関数を上記のように設定することで、塑性の蓄積により亀裂が生じる。

漸近解析とフェーズフィールド変数の挙動

　フェーズフィールド方程式は、以下のように変分原理により導かれる。

\begin{align}
- l_0^2 \nabla^2 \phi + \phi - \frac{2q l_0}{G_c}(1-\phi)^{2q-1}\psi_e^{+} =0
\end{align}

ここで、変数を$L,A,E$を用いて

\begin{align}
x_i/L \rightarrow x_i \ , \  u_i/A \rightarrow u_i \ , \   D_{ijkl}/E \rightarrow D_{ijkl} \ , \  
l_0/L \rightarrow l_0 \ , \  
\end{align}

のように変換する。$L$と$A$は長さ、$E$は圧力の次元を持つ。すると

\begin{align}
- l_0^2 \nabla^2 \phi + \phi - \frac{2q l_0 L}{G_c}(1-\phi)^{2q-1}E\left(\frac{A}{L}\right)^2\psi_e^{+} =0
\end{align}

ここで、

\begin{align}
A = \sqrt{\frac{G_c L}{2E}}
\end{align}

とすれば、

\begin{align}
- l_0^2 \nabla^2 \phi + \phi - l_0 q(1-\phi)^{2q-1}\psi_e^{+} =0
\end{align}

となる。ここで$l_0$は、亀裂の幅を表すため $0 \leq l_0 << 1$である。フェーズフィールド変数$\phi$を$l_0$で展開する。

\begin{align}
\phi = \phi_0 + l_0 \phi_1 + l_0^2 \phi_2 + \ldots
\end{align}

これをフェーズフィールド方程式に代入して、恒等式から

\begin{align}
\mbox{0次の項} : \phi_0 &= 0  \\
\mbox{1次の項} : \phi_1 &= q(1-\phi_0)^{2q-1}\psi_e^{+}  \\
\mbox{2次の項} : \phi_2 &= \nabla^2 \phi_0 - q(2q-1)(1-\phi_0)^{2q-1}\phi_1 \psi_e^{+}
\end{align}

が得られる。$\phi$の展開式に代入すると、

\begin{align}
\phi &= l_0 q\psi_e^{+} - (2q-1)(l_0q\psi_e^{+})^2 + \ldots \\
\end{align}

となる。変数をもとに戻すと、

\begin{align}
 l_0 q\psi_e^{+} \rightarrow \frac{l_0}{L} \frac{L^2}{EA^2} l_0 q\psi_e^{+}= \frac{2l_0}{G_c} q\psi_e^{+}
\end{align}

　従って、フェーズフィールド変数は

\begin{align}
\phi \approx \frac{2l_0}{G_c} q\psi_e^{+}
\end{align}

のように振る舞う。式から

$\psi_e^{+}$と$q$が大きければ亀裂が発生しやすい
$G_c$が大きければ亀裂が発生しにくい

が確認できる。$q$は累積塑性ひずみをしきい値で割った値なので、塑性変形が起きることでフェーズフィールド変数が発展することが確認できる。

残差と剛性行列

　変位とフェーズフィールド変数の離散化は、形状関数を使い

\begin{align}
{\bf u} = \sum_I N_I {\bf u}_I \ \ \ , \ \ \  
\phi = \sum_I N_I \phi_I 
\end{align}

と表現する。

　また、スタッガード法を使う。スタッガード法は、変位の方程式を解く際はフェーズフィールド変数を固定し解き、

\begin{align}
K_{IJ}^u|_{\phi} {\bf u}_J = - {\bf R}_I^u |_{\phi}
\end{align}

フェーズフィールド変数の方程式を解く際は変位を固定し解く。

\begin{align}
K_{IJ}^{\phi}|_{u} {\bf \phi}_J = - {\bf R}_I^{\phi} |_{u}
\end{align}

これを繰り返し収束させる方法である。

変分原理

　変分原理により

\begin{align}
\delta\Psi = 0
\end{align}

となる変位$u$とフェーズフィールド変数$\phi$を求めたい。（汎関数）微分を行えば残差は以下となり、

\begin{align*}
{\bf R}_I^u = \frac{\partial}{\partial {\bf u_I}}\Psi &= \int_{\Omega} \mathrm{d}\Omega B_I {\bf \sigma} 
    - \int_{\Omega} \mathrm{d}\Omega N_I {\bf b} - \int_{\partial\Omega} \mathrm{d}S N_I {\bf f} 
\\
R_I^{\phi} = \frac{\partial}{\partial \phi_I}\Psi &= \int_{\Omega} \mathrm{d}\Omega G_c\left\{l_0(\nabla N_I)\cdot\nabla \phi + \frac{1}{l_0} N_I \phi \right\} \\
&- \int_{\Omega} \mathrm{d}\Omega 2q N_I (1-\phi)^{2q-1} H^+ 
\end{align*}

剛性行列は以下となる。

\begin{align*}
K_{IJ}^u &= \frac{\partial^2}{\partial {\bf u_I} \partial {\bf u_J}}\Psi = \int_{\Omega} \mathrm{d}\Omega B_I D B_J 
\\
K_{IJ}^{\phi} &= \frac{\partial^2}{\partial \phi_I \partial \phi_J}\Psi \\
&= \int_{\Omega} \mathrm{d}\Omega \left\{G_cl_0(\nabla N_I)\cdot(\nabla N_J) + \frac{G_c}{l_0} N_I N_J \right\} \\
& + \int_{\Omega} \mathrm{d}\Omega 2q(2q-1) (1-\phi)^{2q-2} N_I N_J H^+ 
\end{align*}

Kuhu-Tucker 条件を満たすように計算しなければならないので、破壊に寄与するエネルギーは以下とする。

\begin{align*}
H^+ = \max_{s\in [0,t]} \psi^+_e(s)
\end{align*}

Return Mapping アルゴリズム

　Von Mises モデルの降伏関数は、

\begin{align*}
f(\sigma,\alpha) = \sqrt{3J_2}-(\sigma_Y + h\alpha)
\end{align*}

である。偏差ひずみを

\begin{align*}
e_{ij} = \varepsilon_{ij}-\frac{1}{3}\mathrm{Tr}(\varepsilon)\delta_{ij}
\end{align*}

とすれば、第2不変量$J2$は

\begin{align*}
J_2 = 2(\mu g)^2 e_{kl}e_{kl}
\end{align*}

である。

　発展則は塑性乗数$\gamma$を使い

\begin{align*}
{\dot \varepsilon^{p}}_{ij} &= {\dot \gamma}\frac{\partial f}{\partial \sigma_{ij}} \\
{\dot \alpha} &= {\dot \gamma} 
\end{align*}

である。離散化を実施すれば、

\begin{align*}
\Delta\varepsilon^{p}_{ij} &= \sqrt{\frac{3}{2}} \Delta\gamma n_{ij} \\
\Delta\alpha &= \Delta\gamma
\end{align*}

とかける。また、

\begin{align*}
n_{ij} = \frac{e_{ij}} {\sqrt{e_{kl}e_{kl}}}
\end{align*}

とおいた。従って、塑性ひずみと硬化変数は

\begin{align*}
\varepsilon^{p}_{ij} &= \varepsilon^{p}_{nij} + \sqrt{\frac{3}{2}} \Delta\gamma n_{ij} \\
\alpha &= \alpha_n +\Delta\alpha 
\end{align*}

と更新する。下付き添字$n$は更新前（Return Mapping を適用前）の値を表す。この式を降伏関数に代入すると塑性乗数

\begin{align*}
\Delta\gamma = \frac{\sqrt{3{J_2}_n}-(\sigma_Y+h\alpha_n ) }{3\mu g + h}
\end{align*}

が求まる。${J_2}_n$は更新前（Return Mapping を適用前）の$J_2$の値である。

　累積塑性ひずみは

\begin{align*}
\varepsilon_{ep}^p(t+\Delta t) &= \sqrt{\frac{2}{3}} \int_0^{t+\Delta t} \mathrm{d}\tau \sqrt{\dot{{\bf \varepsilon}}^p :\dot{{\bf \varepsilon}}^p } \\
&= \varepsilon_{ep}^p(t) + \sqrt{\frac{2}{3}} \int_t^{t+\Delta t} \mathrm{d}\tau \sqrt{\dot{{\bf \varepsilon}}^p :\dot{{\bf \varepsilon}}^p }
\\
&\approx \varepsilon_{ep}^p(t) + \sqrt{\frac{2}{3}}\Delta t \sqrt{\Delta\varepsilon^{p}_{ij}\Delta\varepsilon^{p}_{ij} }
\\
&= \varepsilon_{ep}^p(t)+ \Delta\gamma
\end{align*}

と計算できる。

　上記をまとめると

$f \leq 0$ の場合

　Amor et al. モデルの計算と同じで、弾性行列は以下として計算する。

\begin{align}
D_{ijkl}=\frac{\partial \sigma_{ij}}{\partial \varepsilon_{kl}} 
&= K\left\{gH(\mathrm{Tr}(\varepsilon))+H(-\mathrm{Tr}(\varepsilon)) \right\}\delta_{ij} \delta_{kl} \\
&  + 2\mu g\left(\delta_{ik} \delta_{jl}  - \frac{1}{3} \delta_{ij} \delta_{kl}\right)
\end{align}

応力は、計算されたひずみから以下を計算する。

\begin{align*}
\sigma_{ij} = g(\phi) \{ K\left\langle \mathrm{Tr}(\varepsilon) \right\rangle_{+}\delta_{ij} + 2 \mu e_{ij} \} +K\left\langle\mathrm{Tr}(\varepsilon) \right\rangle_{-}\delta_{ij}
\end{align*}

$f > 0$ の場合

　降伏関数を使い塑性乗数を求める。

\begin{align*}
\Delta\gamma = \frac{f}{3\mu g + h}
\end{align*}

弾性行列は以下として計算する。

\begin{align}
D_{ijkl}=\frac{\partial \sigma_{ij}}{\partial \varepsilon_{kl}} 
&= K\left\{gH(\mathrm{Tr}(\varepsilon))+H(-\mathrm{Tr}(\varepsilon)) \right\}\delta_{ij} \delta_{kl} \\
& +2\mu g\left(1-\frac{3\mu g \Delta\gamma}{\sqrt{3J_2}} \right)\left(\delta_{ik} \delta_{jl}  - \frac{1}{3} \delta_{ij} \delta_{kl}\right)
\\
& + 6(\mu g)^2 \left(\frac{\Delta\gamma}{\sqrt{3J_2}} - \frac{1}{3\mu g + h} \right) n_{ij} n_{kl}
\end{align}

塑性ひずみと硬化変数、累積塑性ひずみを更新する。

\begin{align*}
\varepsilon^{p}_{ij} & \leftarrow \varepsilon^{p}_{ij} + \sqrt{\frac{3}{2}} \Delta\gamma n_{ij} \\
\alpha & \leftarrow  \alpha +\Delta\alpha 
\\
\varepsilon_{ep}^p &\leftarrow \varepsilon_{ep}^p + \Delta\gamma
\end{align*}

ひずみの弾性成分を計算し

\begin{align*}
\varepsilon_{ij} \leftarrow \varepsilon_{ij} - \sqrt{\frac{3}{2}} \Delta\gamma n_{ij} \\
\end{align*}

応力は塑性成分を取り除いたひずみから計算する。

\begin{align*}
\sigma_{ij} = g(\phi) \{ K\left\langle \mathrm{Tr}(\varepsilon) \right\rangle_{+}\delta_{ij} + 2 \mu e_{ij} \} +K\left\langle\mathrm{Tr}(\varepsilon) \right\rangle_{-}\delta_{ij}
\end{align*}

計算結果

　パラメータは以下として、シミュレーションを行った。連立方程式は共役勾配法（CG法）を用いている。

パラメータ	値
せん断弾性係数	$27.28\rm{Gpa}$
体積弾性係数	$71.66\rm{Gpa}$
降伏応力	$0.345\rm{Gpa}$
硬化係数	$0.250\rm{Gpa}$
$G_c$	$9310 \rm{J/m^2}$
$l_0$	$1.6\times 10^{-1} \rm{mm}$
累積塑性ひずみのしきい値	$7.5\times 10^{-3} \rm{mm}$
1メッシュの節点数	$4$
積分点	$4$
ステップ数	$20000$
変位増分	$5.0\times 10^{-5} \rm{mm}$

　解析モデルは以下である。

　以下は累積塑性ひずみの分布であり、左から14000ステップ、16000ステップ、17000ステップである。プロットは積分点の値である。

　以下はフェーズフィールド変数の分布である。

累積塑性ひずみが高い領域で亀裂が生じていることが確認できる。

計算結果

　今回の記事では、延性破壊におけるフェーズフィールド法について紹介した。結合関数$g(\phi)$に累積塑性ひずみの情報を与えることで、塑性域から亀裂が生じる。

　延性破壊におけるフェーズフィールドにおいては、単純に塑性エネルギーを$H^+$に加える方法がある。次回はその方法を紹介したい。