概要

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

リメッシュの方法
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

2次元の正方形プレートの引張

　2次元の正方形プレートの引張シミュレーションの境界条件は以下の図である。（論文から引用）

　変位の増分は、$1.0\times 10^{-5} \rm{mm}$ とした。

以下は、500ステップの左からAmor et al.モデル、Miehe et al.モデル、ハイブリッドモデルである。

以下は、600ステップ

以下は、700ステップ

である。亀裂が進展するタイミングは異なるが、挙動はどのモデルでも同じである。

2次元の正方形プレートのせん断

　2次元の正方形プレートのせん断シミュレーションの境界条件は以下である。（論文から引用）

　変位の増分は、$2.5\times 10^{-5} \rm{mm}$ とした。

　このシミュレーションにおいて、Amor et al.モデルの挙動は他のモデルと異なることが知られている。下図のような挙動となる。

ただし、ステップ数を増やせば問題なく、Amor et al.モデルのステップ数を他のモデルと比べ10倍にした。

以下は、450ステップ、Amor et al.モデルは4500ステップ。

以下は、550ステップ、Amor et al.モデルは5500ステップ。

以下は、700ステップ、Amor et al.モデルは7000ステップ。

以下は、1000ステップ、Amor et al.モデルは途中で収束性が悪くなったため7400ステップを表示（7400ステップくらいで終了）。

まとめ

　フェーズフィールド法を用いた亀裂進展解析を行った。3つモデルを説明し、2次元の正方形プレートの引張・せん断シミュレーションを行った。

　Miehe et al.モデルが亀裂の現象を良く表現していると考えられる。Miehe et al.モデルを取り扱いやすくしたハイブリッドモデルも少し挙動は異なるが、亀裂の現象を良く表現している。

　Amor et al.モデルは、弾性エネルギーの分割において解釈しやすい。
（引張体積ひずみと偏差ひずみは損傷を受け、圧縮体積ひずみは損傷をうけない）

今回の記事では弾性体を扱ったが、延性破壊では弾塑性を考慮する必要がある。次回は、延性破壊におけるフェーズフィールド法の亀裂進展を紹介したい。