フェーズフィールド法を用いた亀裂進展解析

Last updated at 2025-09-09Posted at 2025-09-05

概要

　亀裂進展解析では、代表的な方法として以下の３点が挙げられます。

リメッシュの方法
XFEM（拡張有限要素法）
ペリダイナミックス

　リメッシュ法とXFEMは、亀裂進展時の幾何学的な処理が複雑であることが知られています。ペリダイナミックスは粒子法の一つで、支配方程式が微分方程式ではなく積分方程式であるため、積分ができない領域（境界部分）の取り扱いが困難とされています。

　今回の記事では、XFEMに比べ幾何学的な処理がなく、ペリダイナミックスに比べ境界部分の処理が容易である（FEMと一緒である）フェーズフィールド法を用いた解析の紹介をしたいと思います。

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

フェーズフィールド法における亀裂の簡単な説明

　$x=0$ の点で亀裂が生じているとする。このとき、亀裂を表現する関数は以下のように表せる。

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

　フェーズフィールド法の基本的な考えは、この不連続な亀裂を幅$l_0$とした指数関数で近似する。

\begin{align}
\phi(x) = \exp\left(-\frac{|x|}{l_0} \right)
\end{align}

　この指数関数は、以下の境界条件

$x=0$で亀裂が生じていること
無限遠ではゼロであること

\begin{align}
\phi(0) = 1 \ ,\ \phi(\pm\infty) = 0
\end{align}

を満たすオイラー方程式

\begin{align}
\phi(x) - l_0^2 \frac{\mathrm{d}^2\phi(x)}{\mathrm{d} x^2} = 0
\end{align}

の解であることがわかる。これを３次元に拡張して、破壊エネルギーは以下で表せる。

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

$\Gamma$は亀裂面を表す。$G_c$は臨界エネルギー解放率であり、$l_0$は亀裂の幅に相当する。FEMの解析の際は、メッシュサイズより大きな値にする必要がある。

全エネルギーと支配方程式

　有限要素方程式を導出するために、全エネルギーを記述する。全エネルギー$\Psi$は以下となる。

\begin{align}
\Psi = \int_{\Omega} \mathrm{d}\Omega \psi_e + \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項：表面力による仕事

　ひずみエネルギーは、破壊に寄与するエネルギー$\psi_e^+$と残留エネルギー$\psi_e^-$ に分離し、フェーズフィールド変数と相互作用させ

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

とする。$k$は数値的不安定を回避するためのパラメータである。

　変分原理により

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

となる変位$u$とフェーズフィールド変数$\phi$を求めれば良い。これより支配方程式は

\begin{align}
&\nabla \cdot {\bf \sigma} + {\bf b} =  {\bf 0}\\
&\left[2(1-k)H + \frac{G_c}{l_0} \right]\phi - G_c l_0 \nabla^2 \phi = 2(1-k)H 
\end{align}

と求まる。応力は

\begin{align}
{\bf \sigma} = \frac{\partial \psi_e}{\partial \varepsilon}
\end{align}

である。$H$は純粋な変分では$ \psi_e^+$となるが、亀裂がもとに戻ることはないことを考慮する必要がある。従って、以下のKuhu-Tucker 条件を満たす必要がある。

\begin{align}
\psi^{+}_{e} - H \leq 0, \quad \dot{H} \geq 0, \quad \dot{H} (\psi^{+}_{e} - H) = 0
\end{align}

第1項：不可逆性条件
第2項：単調性条件
第3項：相補性条件

プログラム実装では、以下とすれば良い。

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

$s$は荷重ステップである。Kuhu-Tucker 条件の1項目と２項目は満たす事がわかる。3項目は、現在のステップ$t$でエネルギーが最大なら

\begin{align}
H = \psi^+_e(t)
\end{align}

であり、過去のエネルギーのほうが大きいならば$H$は更新されないので

\begin{align}
\dot{H} = 0
\end{align}

となり満たされる。
　他の方法としては、ペナルティ法がある。

有限要素方程式とスタッガード法

　次に離散化を行う。変位とフェーズフィールド変数を形状関数を使い

\begin{align}
{\bf u} = \sum_I N_I {\bf u}_I \ \ \ , \ \ \  
\phi = \sum_I N_I \phi_I 
\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\} - 2\int_{\Omega} \mathrm{d}\Omega (1-k)H N_I (1-\phi)
\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) + \left(\frac{G_c}{l_0} +2(1-k)H \right)N_I N_J \right\}
\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}
 \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}

は体積弾性率である。

Amor et al. model

　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}

すると応力は、

\begin{align}
\sigma_{ij} = \{ (1-k)(1-\phi)^2+k \}\{ K\left\langle \mathrm{Tr}(\varepsilon) \right\rangle_{+}\delta_{ij} + 2 \mu \varepsilon^{div}_{ij} \} +K\left\langle\mathrm{Tr}(\varepsilon) \right\rangle_{-}\delta_{ij}
\end{align}

そして弾性行列は

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

となる。

Miehe et al. model

　Miehe et al.モデルは、ひずみを固有値分解を行い

\begin{align}
{\bf \varepsilon} _{\pm} = \sum_{a=1}^d  \left\langle \varepsilon_{a} \right\rangle _{\pm} {\bf n}_a \otimes {\bf n}_a 
\end{align}

弾性エネルギーを以下のように分ける。

\begin{align}
\psi_e^{\pm} = \frac{\lambda}{2} \left\langle \mathrm{Tr}(\varepsilon) \right\rangle _{\pm} ^2 + \mu \mathrm{Tr}(\varepsilon^2_{\pm} )
\end{align}

応力は

\begin{align}
{\bf \sigma} =  \{ (1-k)(1-\phi)^2+k \} &\left[\lambda\left\langle \mathrm{Tr}(\varepsilon) \right\rangle _{+} {\bf I}+ 2\mu \varepsilon_{+}\right] \\
&+ \lambda\left\langle \mathrm{Tr}(\varepsilon) \right\rangle _{-} {\bf I}+ 2\mu \varepsilon_{-}
\end{align}

そして弾性行列は

\begin{align}
D_{ijkl} &= \bar{D}_{ijkl} + \tilde{D}_{ijkl} \\
\bar{D}_{ijkl} &= \lambda \left[ \left\{ (1 - k)(1 - \phi)^2 + k \right\} H_{\epsilon}(\text{tr}(\epsilon)) + H_{\epsilon}(-\text{tr}(\epsilon))\right] \delta_{ij} \delta_{kl}  \\
\tilde{D}_{ijkl} &= 2\mu \left[\left\{(1 - k)(1 - \phi)^2 + k \right\}P^+_{ijkl} + P^-_{ijkl}\right]
\end{align}

\begin{align}
P_{ijkl}^{\pm} = \sum_{a=1}^{3} \sum_{b=1}^{3} & H_{\epsilon}(\pm\epsilon_a) \delta_{ab} n_{ai} n_{aj} n_{bk} n_{bl} \\
&+ \frac{1}{2} \sum_{a=1}^{3} \sum_{b \neq a}^{3}
\frac{\langle \epsilon_a \rangle_{\pm} - \langle \epsilon_b \rangle_{\pm}} {\epsilon_a - \epsilon_b} n_{ai} n_{bj} (n_{ak} n_{bl} + n_{bk} n_{al})
\end{align}

である。固有値が重複する場合は発散してしまうので

\begin{align}
\begin{cases}
\epsilon_1 = \epsilon_1(1 + 10^{-9}) & \text{if } \epsilon_1 = \epsilon_2 \\
\epsilon_3 = \epsilon_3(1 - 10^{-9}) & \text{if } \epsilon_2 = \epsilon_3
\end{cases}
\end{align}

の処理を行う。

ハイブリッド法

　Miehe et al.モデルは、等方モデルと比べて非線形で取り扱いが難しい。ハイブリッド法は、変位を計算するときは$\psi_e^+$と$\psi_e^-$に分離せず計算する。つまり、等方モデルを用いる。応力は

\begin{align}
\sigma_{ij} = \{ (1-k)(1-\phi)^2+k \}\left\{ K \mathrm{Tr}(\varepsilon) \delta_{ij} + 2 \mu \varepsilon^{div}_{ij} \right\}
\end{align}

弾性行列は

\begin{align}
D_{ijkl}=\frac{\partial \sigma_{ij}}{\partial \varepsilon_{kl}} 
&= \{ (1-k)(1-\phi)^2+k \}\left\{K\delta_{ij} \delta_{kl} + 2\mu \left(\delta_{ik} \delta_{jl}  - \frac{1}{3} \delta_{ij} \delta_{kl}\right)\right\}
\end{align}

で計算する。
　フェーズフィールド変数を更新する際は、Miehe et al.モデルの$\psi_e^+$を使い計算する。ハイブリッド法の変位を求める際の剛性行列は線形であり、スタッガード法の収束性が良くなる。

計算結果

パラメータは以下として、シミュレーションを行った。

パラメータ	値
ヤング率	210 $\rm{Gpa}$
ポアソン比	0.3
$G_c$	2700$\rm{J/m^2}$
$l_0$	$1.5\times 10^{-2} \rm{mm}$
メッシュサイズ	$7.5\times 10^{-3} \rm{mm}$
1メッシュの節点数	4
積分点	4
ステップ数	1000