数理解析研究所講究録
第 1791 巻 2012 年 182-188

182
き裂進展フェーズフィールドモデルとその応用について

広島国際学院大学・情報デザイン学部 高石 武史 (Takeshi Takaishi )Faculty of Information Design,
Hiroshima Kokusai Gakuin University

1

き裂進展を記述するフェーズフィールドモデル
筆者と木村は板状物質において板面の垂直方向の変位 (モード III) で進展するき裂を，Francfor-Marigo[5] のエネルギー表式を用い，き裂の有無を表すフェーズフィールドの時間発展を板の変位の満たす方程式と組み合わせた，2 変数の反応拡散方程式で表現できることを示した [11,12]. また，この方程式を解くことで，特殊な数値解法を必要とせずに亀裂進展現象の数値シミュレーションができることを実際に示した．

(1)$\{\begin{array}{l}\alpha_{1^{\frac{\partial u}{\partial t}}} = \mu div((1-z)^{2}\nabla u)+f x\in\Omega, t>0\alpha_{2}\frac{\partial z}{\partial t} = (\epsilon div(\gamma(x)z)-\frac{\gamma(x)}{\epsilon}z+\mu|\nabla u|^{2}(1 \text{一} z))_{+} x\in\Omega, t>0u(x, t) = g(x,t) x\in\Gamma_{D}, t>0\frac{\partial u}{\partial n}=0 x\in\Gamma_{N}, t>0, \frac{\partial z}{\partial n}=0 x\in\Gamma, t>0u(x,0)=u_{0}(x), z(x,0)=z_{0}(x) x\in\Omega\end{array}$

ここで， を有界な 2 次元領域で，その境界 は区分的に滑らかであるとする．また， は正$t\geq 0$ は時刻，の長さを持つ有限個の連結成分からなるの閉部分集合で，とおく．
$\Omega$$\Gamma$$\Gamma_{D}$$\Gamma$$\Gamma_{N}:=\Gamma\backslash \Gamma_{D}$$u(x,t)$ は $x\in\Omega$ における物質の板状面に対して垂直な$f(x,t)$ は板状面に対して垂直方向の外力，$g(x,t)$ は境界方向の微小変位であり，において与えられた垂直方向の変位である．変数 $z(x,t)$はき裂の有無を表すフェーズフィールドである．また， は正規化のための微小パラメータ，は破壊しづらさを表す破壊靭性値， は Lam\’e 定数の一つである．
$\Gamma_{D}$$\epsilon$$\gamma(x)$$\mu$

このモデルでは純粋な板のき裂を扱っているが，実際には他の要因との相互作用でき裂が進んでいく場合が多い．(1) は比較的シンプルな方程式であるため，そこへ現実的ないろいろな効果を付与していくことでモデルを拡張することが可能である．ここでは，2 においてき裂進展モデルの中で，このモデルがどのように位置づけられるかについて延べ，3 では，応用として，き裂から破断に至る前段階としての微小き裂のモデル化について検討し，様々なき裂進展現象へのこのモデルの摘要の容易さについて述べる．

2モデルの位置づけ

Griffith は，仮想的にき裂面の形成を考えた際に，き裂面の形成に必要とされるエネルギー量が，き裂面の形成による周囲の物体から解放される弾性エネルギー量とが釣り合った場合に，実際にきlE-mail: t.takaishi@hkg.ac.jp

183
裂として実現されると考えた [6]. また，Griffith は他にも化学エネルギー，熱エネルギー，運動エネルギー等がき裂進展に伴いバランスしていると考えた．近代のき裂進展に関する研究はこの考え方をベースにいろいろ発展してきた (エネルギー変分を元にしたき裂の形状に関する研究は，例えば [10] 等).

Francfort と Marigo は，Griffith の考え方を基に，対象とする物質全体に対して，き裂面形成により増加する系のエネルギーと，バルクの (つまりき裂面以外での) 弾性エネルギーとの和に注目した．簡単のために平板のモード III き裂進展に話を絞ることにする．系全体を , き裂面をとし，次のエネルギー表式を提案した [5].
$\Omega$$\{\begin{array}{ll}E(\Sigma) = E_{1}(u, \Sigma)+E_{2}(\Sigma), E_{1}(u, \Sigma) ;= \min E_{1}(v, \Sigma) (u\in V(g, \Omega\backslash \Sigma)),E_{2}(\Sigma) ;= \int_{\Sigma}\gamma(x)dsv\in V(g’\Omega\backslash \Sigma).

$\Sigma$(2)\end{array}$

彼らは，この合計エネルギーが最小になるようにき裂が形成されると考え，エネルギー変分から準静的にき裂が進展する様子を表現できると提案した．
さてこのモデルから実際にき裂進展の計算をするのは難しいが，Ambrosio と Tortorelli のアイデアによって利用に向かって大きく前進した．彼らは Mumford-Shah functional の弱形式を考え，正則化パラメータを導入することで image segmentation(画像の抽出処理) を行えることを示した[1]. 彼らの考えたモデルは，正則化パラメータをすることで元のモデルに収束することがわかっており，彼らの提案したエネルギーの最小状態を調べることで元のモデルにおけるエネルギーの最小状態を見つけることができる．
$\epsilonarrow 0$

Bourdin, Francfort. Marigo は，物体の微小変位$u$と補助変数$\Gamma$$z$を用いて，Francfort-Marigoのエネルギー表式を Ambrosion-Tortorelli の正則化の手法で次のようなエネルギー表式に書き直し，き裂形状を実際に数値計算できることを示した [3,4].
(3)$E_{\epsilon}(u, z):= \frac{\mu}{2}\int_{\Omega}(1-z)^{2}|\nabla u|^{2}dx-\int_{\Omega}fudx+\frac{1}{2}\int_{\Omega}\gamma(x)(\epsilon|\nabla z|^{2}+\frac{1}{\epsilon}z^{2})dx$

尚，Bourdin らはという表式を用いているが，ここでは筆者らのモデルに合わせて書き直している．ここで，スカラー変数 が Kachanov の damage parameter[7] と似た役割を果たしているが，ここではき裂面で $z\sim 1$ , それ以外では $z\sim 0$ をとるものとして考えている．彼らは，初期状態と境界条件を与えこのエネルギーのグローバルな最小状態を求めることで，き裂の形状を計算している．このモデルをベースにした数値計算は，Bourdin[2] のみならず，いろいろ行われている ([9] 等).
さて，ここまでのモデルでは準静的な脆性き裂を考えており，各離散時間においてき裂は瞬間的に (エネルギー最小状態をとることで) 有限長伸びることも可能だとしている．しかし，現実的にはそうであろうか．木村と私は，き裂は時間とともに連続的に伸びていくもので，その形状は現時点でのき裂形状に依存し，局所的にエネルギー最小となるパスを探しながら進展するのではないかと考えた．そこで，時間発展の Ginzburg-Landau の理論を適用し，エネルギー (3) を系の状態を記述するフリーエネルギーとし，このエネルギーが各時刻において準平衡状態にあると考え，その上でエネルギーの勾配に従い系の状態が時間発展するものとした．また，き裂の有無を表しているは物体の状態を表すフェーズフィールドと考えることができ，奇しくもフェーズフィールド法で時$\iota\ovalbox{\tt\small REJECT}=1-z$

$z$

$z$

間発展方程式導出の際に用いられる時間発展の Ginzburg-Landau の理論を適用することとなった(例えば [8] 等). 外力 $f,g$ は時間変化が非常に緩やかで，微小変位$u$とき裂の有無を表す$z$は比184
Griffith

gradient-flow
図 1: それぞれのモデルの位置づけ

較的短い時間で準平衡状態に到達できると仮定する．さらに，き裂が修復しない条件
加$)$$(z$は単調増と言う条件を加味すると，系の平衡状態近くでのダイナミクスは次のように表すことができる．

$\{\begin{array}{l}\alpha_{1}\frac{\partial u}{\partial l} = -\frac{\delta E_{\epsilon}(u,z)}{\delta u}\alpha_{2}\frac{\partial z}{\partial t} = [-\frac{\delta E_{\epsilon}(u,z)}{\delta z}]_{+}\end{array}$(4)

ここで，は適当な時定数である．その結果 (1) が導出され，モード III き裂進展を 2変数の反応拡散方程式で表現していることから，数値計算の容易さのみならず．き裂形状や進展速度などの理論的な解析が容易であることも期待される．またこのモデルから得られるき裂形状は，$\alpha_{1}\geq 0,$$\alpha_{2}>0$

系全体の最低エネルギー状態の実現ではなく，き裂進展の途中経過に依存することがわかっている[12].
Bourdin らの用いたモデルは$\{\begin{array}{l}(u_{i},z_{i}) = argmin E_{\epsilon}(u, z)u=g(t_{:}),z=1on\Sigma_{i-1}E_{\epsilon}(u, z) ;= \frac{\mu}{2}\int_{\Omega}(1-z)^{2}|\nabla u|^{2}dx+\frac{\gamma}{2}\int_{\Omega}(\epsilon|\nabla z|^{2}+\frac{z^{2}}{\epsilon})d_{X}\Sigma_{i} ;= \{x|z_{i}\leq 1-\eta\}\end{array}$(5)となるのに対し，筆者らのモデルは$(u_{i}, z_{i})= u=g(t\dot{.}),z\geq z:-1\arg\min(E_{\epsilon}(u,z)+\frac{1}{2}\frac{\alpha_{1}||u-u_{i}||_{L^{2}(\Omega)}^{2}+\alpha_{2}||z-z_{i}||_{L^{2}(\Omega)}^{2}}{t_{i}-t_{i-1}})$(6)と書き表せる．Bourdin らが様々な手法を用いて global minimizer を数値的に探しているのに比べ，local minimizer を探す筆者らのモデルは数値計算が圧倒的に簡単になっている．
一般に，き裂進展はき裂を開く方向への変位 (モード I, 開ロモード) に大きく依存することが知られている．ここまでモード III のき裂進展について見てきたが，同じ手法でモード I, II も含む
き裂進展の時間発展方程が導出できる．

185$t=0$

$t=1$

$t=2$

図 2: モード I,II のき裂進展におけるフェーズフィールドの空間プロファイル．

$\{\begin{array}{ll}\alpha_{i}\frac{\partial u_{i}}{\partial t} = \mu(div((1-z)^{2}\nabla u_{i})+\frac{\partial}{\partial x_{j}}((1-z)^{2}\frac{\partial u_{j}}{\partial x_{i}})+\mu\frac{\partial}{\partial x_{i}}((1-z)^{2}divu))+f_{i} x\in\Omega, t>0i=1,2,3 \alpha_{z}\frac{\partial z}{\partial t} = (\epsilon div(\gamma(x)z)-\frac{\gamma(x)}{\epsilon}z+(1-z)(\lambda(divu)^{2}+\mu(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}})^{2}))_{+} x\in\Omega, t>0with B.C and I.C. \end{array}$

(7)

ここで，

$u=(u_{1}, u_{2}, u_{3})$3$f=(f_{1}, f_{2}, f_{3})$ は板に作用する外力である．
は板の変位，微小き裂を含んだ材質のモデル化
ここまで， はき裂の有無を表現するフェーズフィールドとして考えてきたが，単調増加関数と$z$して扱っているため，数値シミュレーション結果から$z$

はき裂近傍の損傷も表していることがわかった．そこで， を damage parameter として再定義し，モデルを微小き裂を含んだ物体での経年変化へ応用していくことにする．
$z$

超微細なき裂 (マイクロクラック) は固い物質内によく見られる現象で，シリコンウェハーの製品不良や，歯のエナメル質の脆化等，大きなき裂へ進展することによる弊害を含む場合が多い．材質にき裂が生じる前段階としてまず微小き裂が発生し，それらが集まり成長することで貫通き裂となり，さらにそのき裂が進展することで材質の破断に至ると考えることができる．

筆者らのモデルにおいては，フェーズフィールドからき裂断面積を見積り，結晶内の原子が表面にさらされることによるエネルギー上昇を勘案している．この際に，微小き裂においては i) 非常に微細な貫通き裂としての見方と，ii) エネルギー上昇が中間的な値をとる状態としての見方が考えられる．前者では，このモデルで初期状態の設定として無数の欠陥を与えて計算することで実現が可能である．後者は，本来ほとんどの領域で $z\sim O$ か $z\sim 1$ という値をとるべきフェーズフィール
ドに中間的な値を設定することで，転移のようにき裂に至らないエネルギー上昇を表現することができる．そのように考えることで，貫通き裂が成長する過程から，破断に至るまでの一貫した解析が可能になることが期待される．数値シミュレーションの結果，初めに微小なき裂が連結していったのち，進展して貫通き裂に成長し，その後でき裂の進展により破断に至る過程が再現された．
Griffith の理論ではき裂面の形成によって系のエネルギーが上昇し，一方，き裂の進展による材質内部の弾性エネルギーが減少することから，このエネルギーのバランスでき裂の進展が決定すると考えられている．Rancfort-Marigo のエネルギーをベースとする筆者らのフェーズフィールド

186

図 3: き裂がある場合 (実線) と，(b) 微小き裂がある場合 (破線) におけるフェーズフィールドの空間プロフアイル．

モデルにおいては，結晶内部のエネルギー的に安定な原子に対して，き裂面では表面原子となるために結合の腕が余ることからエネルギー状態が上がることを勘案し，き裂面の面積 (2 次元問題では長さ) に比例してエネルギーが上昇すると仮定している．ここまで，フェーズフィールドはなめらかな関数ではあるが，材質中のほとんどの領域で z(き裂あり) のO(き裂なし), または状態であるとした．ここで，フェーズフィールドの定義づけを拡張し，転移，結晶構造の歪みなどによるエネルギーの微小な増加を扱うために，エネルギー的に中間的な状態 $0<z<1$ を取りうるものとして考えることにする．一方で，微小き裂とは非常に小さいき裂として考えることも可能$\approx$$z\approx 1$である．

(a)

(b)

(c)

図 4: 初期微小き裂 $z(x, 0)=0.5 \max(\cos(10x)\cos(20y), 0)$ の空間プロファイル (a) と破断の様子( $u(b),$$z(c)$. 左側面中央の初期き裂が進展).

$t=0$ において微小なき裂 $z(x, 0)=0.5 \max(\cos(10x)\cos(20y), 0)$ が存在したとして，$\alpha_{1}=$$0,$$\alpha_{2}=0.01,$ $\epsilon=0.01,$ $\gamma=0.5,g(x)=10,\Gamma_{D}=\{(x_{1},x_{2})|x_{1}\in(-1,1), x_{2}=\pm 1\}$て数値シミュレーションを行った．ALBERTA toolbox を使用することで，時間方向は陰解法のアダプティブメッシュ有限要素法で計算した．初期値において結晶構造の歪み (不完全な穴)$\tilde{z}_{0}=0.5$$($

図 $- 4 (b))$

が存在すると，き裂が進展し破断に至ったところ，断面に凸凹が生じていることが分かる．
このように，フェーズフィールドの意味づけを拡張することでき裂の前段階であるマイクロク

ラックを扱うことができれば，材質の経年変化の予測など広範囲に適用することが期待される．また，化学反応によるエネルギーの変化や破壊靭性値の変化を考慮することで，金属材料におけるさ187
びや化学物質との反応による脆化がき裂の進展にどのような影響を及ぼすかといった研究への適用も同様に期待される．筆者らのモデルは，平衡に近い状態にある系がエネルギーの勾配流によって時間発展するという単純な原理で構築されているため，様々な応用が可能になっている．

4

まとめ

き裂進展現象を数値計算するためには，き裂面に起因する特異性や境界面の変更や，陽的なき裂進展方向の特定方法が無いことから，多くの計算方法が提案されてきた．筆者らのフェーズフィー
ルドを用いたモデルでは，このモデルはエネルギー勾配に従う時間発展の方程式であることから種々の効果の考慮が容易であり，また，数値シミュレーションに適しているため，大規模な数値シミュレーションや迅速な予測等での応用が期待できる．

まだまだモデルの検証の途上ではあるが，数値シミュレーションに好適なモデルとして今後が期待される．

参考文献[1] L. Ambrosio and V. M. Tortorelli, On the approximation of free discontinuity problems,Boll. Un. Mat. Ital. (7) $6-B$ , (1992) 105-123.
[2] B. Bourdin, Numerical implementation of the variational formulation of brittle fracture,
Interfaces Free Bound. 9 (2007), 411-430.

[3] B. Bourdin, G. A. Rancfort and J.-J. Marigo. Numerical experiments in revisited brittlefracture, J. Mech. Phys. Solids 48 (2000), 797-826.

[4] B. Bourdin, G. A. Rancfort and J.-J. Marigo, The variational approach to fracture. (2009),
Springer.
[5] G. A. Francfort and J.-J. Marigo, Revisiting brittle fracture as an energy minimization
problem, J. Mech. Phys. Solids 46 (1998), 1319-1342.

[6] A. A. Griffith, The phenomenon of rupture and flow in solids, Phil. TYans. Royal Soc. London
A221 (1920), 163-198.
[7] L. M. Kachanov, On the life-time under creep conditions, Izv. Akad. Nauk SSSR 8 (1958),

26-31.
[8] R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Physica
63 (1993), 410-423.

$D$[9] G. Launcioni and G. Royer-Carfagni, The variational approach to fracture mechanics, J.
Elast. 95 (2009), 1-30.
[10] Y. Sumi, S. Nemat-Nasser and L. M. Keer, On crack path stability in a finite body, Eng.
Frac. Mech. 22 (1983), 759-771.

188[11] T. Takaishi and M. Kimura, Phase field model for mode III cmck growth. Kybernetika 45(2009), pp605-614.
[12] 高石武史，モード III 亀裂進展のフェーズフィールドモデルとその数値計算．日本応用数理学会論文誌 19 (2009). 351-369.