フェーズフィールド法を用いた延性破壊における亀裂進展解析 2/2

Posted at 2025-09-06

概要

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

　破壊には大きく分けて２つあり、脆性破壊と延性破壊がある。延性破壊の場合、塑性変形が起こり破壊に至ります。前回の記事では、結合関数に累積塑性ひずみの情報を与えることで塑性域から亀裂が生じるモデルについて紹介しました。

　今回は塑性エネルギーを破壊に寄与するエネルギーに直接組み込むモデルを紹介したいと思います。以下の文献を参考にしています。

最初に

　$\lambda$と$\mu$はラメ定数で

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

は体積弾性率である。偏差ひずみは

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

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

\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}

これを繰り返し収束させる。

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

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

\begin{align}
\Psi = \int_{\Omega} \mathrm{d}\Omega g(\phi)\{\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}

弾性エネルギーは分割せず、

\begin{align*}
\psi_e = \frac{1}{2}K \mathrm{Tr}(\varepsilon)^2 + \mu e : e
\end{align*}

であり等方モデルを用いる。塑性エネルギーは

\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*}

結合関数

　結合関数は以下の3次関数を用いる。

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

$k$は数値的不安定を回避するためのパラメータである。また$0\leq \beta \leq 2$ である。（$g$は減少関数の必要があるから）

　これまでのqiitaの記事では$\beta=2$に対応する2次関数を用いてきたが、$\beta=2$では劣化するのが早く、塑性が起こる前から破壊が起きてしまう。そのため、$\beta$は0に近い値をとる。

残差と剛性行列

　Kuhu-Tucker 条件を満たすように計算しなければならない。さらに、塑性域から亀裂が生じてほしいので、破壊に寄与するエネルギーは以下とする。

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

残差

\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 g^{\prime}(\phi) N_I 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 g^{\prime\prime}(\phi) N_I N_J H^+ 
\end{align*}

Return Mapping アルゴリズム

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

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

第2不変量$J2$は

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

である。
　発展則は塑性乗数$\gamma$を使い、離散化を実施すれば、

\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*}

とおいた。

　まとめると

$f \leq 0$ の場合

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

\begin{align}
D_{ijkl} = Kg(\phi) \delta_{ij} \delta_{kl} + 2\mu g(\phi) \left(\delta_{ik} \delta_{jl}  - \frac{1}{3} \delta_{ij} \delta_{kl}\right)
\end{align}

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

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

$f > 0$ の場合

　降伏関数を使い塑性乗数を求める。降伏関数に出現する結合関数$g$は、塑性乗数を計算する際は打ち消し合う。

\begin{align*}
\Delta\gamma = \frac{f}{g(\phi)(3\mu + h)}
= \frac{\sqrt{6\mu^2 e_{kl}e_{kl}}-(\sigma_Y + h\alpha)}{3\mu + h}
\end{align*}

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

\begin{align}
D_{ijkl}&= Kg(\phi)\delta_{ij} \delta_{kl}
+2\mu g(\phi)\left\{1-\sqrt{\frac{3}{2}}\frac{\Delta\gamma}{\sqrt{e_{mn}e_{mn}}} \right\} \left\{\delta_{ik} \delta_{jl}-\frac{1}{3}\delta_{ij} \delta_{kl} \right\} \\
& + \sqrt{6}g(\phi)\mu \left(\frac{\Delta\gamma}{\sqrt{e_{mn}e_{mn}}} - \frac{\sqrt{6}\mu}{3\mu + 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 
\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) \left\{ K\mathrm{Tr}(\varepsilon)\delta_{ij} + 2 \mu e_{ij} \right\} 
\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}$
$\beta$	$1.0\times 10^{-6} $
1メッシュの節点数	$4$
積分点	$4$
ステップ数	$16500$
変位増分	$5.0\times 10^{-5} \rm{mm}$