表面エネルギーのDFT計算とPhase Field法の連携によるLaB₆セラミックスの劈開破壊解析

六ホウ化ランタン（$\text{LaB}_6$）は、低い仕事関数と高い融点を併せ持つことから、電子顕微鏡や電子線描画装置の熱電子放出源（電子銃カソード）として広く用いられる重要なデバイス材料です。一方で本材料は硬質な脆性セラミックスであり、熱衝撃や機械的負荷による劈開割れ・クラック進展が実用上の課題となります。脆性破壊の予測が難しい本質的な理由は、巨視的な破壊強度が、原子スケールで決まる「表面エネルギー」に直接支配される点にあります。さらに LaB₆では、結晶成長時の合成環境（化学ポテンシャル）に応じて表面終端（La終端／B₆終端）が変化し、表面エネルギーが大きく異なります。本事例では、第一原理計算ソフトウェア Advance/PHASE で算出した LaB₆(001) 面の表面エネルギー $\gamma$ を出発点とし、グリフィスの脆性破壊理論（$G_c = 2\gamma$）を介して、マルチフィジックス有限要素法フレームワーク MOOSE 上の Phase Field 破壊力学モデルへ橋渡しします。これにより、「合成環境の違い（La過剰／B過剰）が、最終製品の巨視的な破壊強度をどれだけ変えるか」について、シミュレーションによる解析結果を紹介します。

Keywords: 第一原理計算（DFT）, Phase Field法, 劈開破壊 (Cleavage Fracture), LaB₆, 表面エネルギー, グリフィス理論, 化学ポテンシャル, マルチスケール解析, MOOSE, 電子銃陰極

解析モデルと理論背景

1. マルチスケール連携のアプローチ

本解析の核心は、原子スケールの第一原理計算 (DFT)と巨視的な連続体破壊力学を、「表面エネルギー」という物理量で連結する点にあります。Advance/PHASE による LaB₆ 表面エネルギーの第一原理計算 [1] では、スラブモデルを用いて (001) 面の表面エネルギーが終端構造と化学ポテンシャルに強く依存することが示されています。すなわち、La過剰（La-rich）環境で安定なLa終端面と、B過剰（B-rich）環境で安定な B₆終端面とでは、表面エネルギー $\gamma$ が約3倍も異なります。

この $\gamma$ を、グリフィスの劈開破壊条件 $G_c = 2\gamma$ により臨界エネルギー解放率へ換算し、Phase Field 破壊モデルの破壊靭性パラメータとして与えます。巨視モデルは図1に示す 1.0 mm $\times$ 1.0 mm の2次元平板（平面ひずみ）とし、全体を(001)面劈開を仮定した LaB₆ 単結晶と見なして、左端中央に長さ 0.2 mm の初期亀裂を設定しました。初期亀裂は、メッシュを幾何学的に切り離したスリットではなく、損傷場 $d$ の初期条件として導入しています（亀裂線に沿って $d \approx 1$ となる、幅 $\sim l$ の拡散的な初期損傷帯）。損傷の不可逆条件（下限拘束）により、この初期亀裂は解析を通じて保持されます。上端に一定速度の引張変位を与え、下端は $y$ 方向を固定、上端の $x$ 方向は自由としてポアソン収縮を許容しています。連成・破壊計算は、オープンソースの有限要素法フレームワーク MOOSE [2] 上で実行しました。

図1. 2次元平板モデルの初期状態。1.0 mm $\times$ 1.0 mm の領域の左端に長さ 0.2 mm の初期亀裂（赤色部分が損傷変数 $d=1$）を設定し、上端に引張変位を与えます。

2. 支配方程式 — 第一原理から巨視破壊への橋渡し

Phase Field法は、亀裂などの不連続な界面を連続的な損傷変数で表現し、メッシュに依存せず複雑な亀裂進展を再現する手法です。本事例では、第一原理計算で得た表面エネルギーを、以下の2段階で巨視破壊モデルに連携させています [3-5]。

A. グリフィス則による橋渡し（ミクロ → マクロ）

劈開破壊では、亀裂が単位面積だけ進展するのに必要なエネルギー（臨界エネルギー解放率 $G_c$）は、新たに生成される2つの自由表面のエネルギーに等しく、次式で与えられます。第一原理計算による $\gamma$ を代入することで、各合成環境に対応する破壊靭性が一意に定まります。

$$G_c = 2\gamma$$

B. Phase Field 破壊力学モデル（マクロ）

亀裂を不連続な面ではなく、連続的な損傷変数 $d$（0で健全、1で完全破壊）として表現します。系全体の全エネルギー $\Pi$ を最小化するように亀裂が進展します。

$$\Pi = \int_{\Omega} \left[ (1-d)^2\, \psi_0^{+}(\boldsymbol{\varepsilon}) + G_c \left( \frac{d^2}{2l} + \frac{l}{2} |\nabla d|^2 \right) \right] dV$$

ここで $\Omega$ は解析領域です。亀裂表面積を近似する汎関数として標準的なAT2（Ambrosio-Tortorelli）モデルを採用しており、$\psi_0^{+}(\boldsymbol{\varepsilon})$ は弾性ひずみエネルギーの引張成分（スペクトル分解により引張側のみが亀裂を駆動し、圧縮下では割れない）、$l$ は亀裂の幅を規定する長さスケールパラメータです。損傷の不可逆性（自己修復の禁止）は、損傷変数 $d$ の単調増加拘束として課しています。

3. 使用した材料パラメータ

純 $\text{LaB}_6$ を対象とし、Advance/PHASE による第一原理計算値 [1] と、実績のある文献値 [6]を組み合わせて解析を行いました。第一原理計算からは、破壊特性を決定づける最重要パラメータである「表面エネルギー」のみを抽出し、弾性定数等は文献値で補うことで、計算コストを最小化しています。

表1. 劈開破壊シミュレーションに用いた主なパラメータ

カテゴリー	パラメータ	設定値と出典
第一原理計算 (Advance/PHASE)	表面エネルギー \(\gamma\)（La終端 / La-rich）	0.814 J/m² [1]
	表面エネルギー \(\gamma\)（\(\text{B}_6\)終端 / B-rich）	2.438 J/m² [1]
グリフィス換算 (\(G_c = 2\gamma\))	臨界エネルギー解放率 \(G_c\)（La-rich）	1.628 J/m²（= 1.63 N/m）
	臨界エネルギー解放率 \(G_c\)（B-rich）	4.876 J/m²（= 4.88 N/m）
弾性・幾何 (文献値[6]等)	ヤング率 (\(E\))	337 GPa
	ポアソン比 (\(\nu\))	0.18
	長さスケール (\(l\))	0.04 mm
	モデル寸法 / 初期亀裂長さ	1.0 mm × 1.0 mm / 0.2 mm

（注）長さスケール $l$ の設定： $l$ は正則化された亀裂帯（破壊プロセスゾーン）の幅を表すモデルパラメータです。本解析では一様メッシュの要素寸法 $h \approx 0.0083\,\text{mm}$（$1\,\text{mm}$ を 120 分割）に対して $l \approx 5h$ となるよう $l = 0.04\,\text{mm}$ とし、幅 $\sim l$ の損傷帯を複数要素で解像しています（一般に $l \gtrsim 2h$ が目安）。なお $l$ は見かけの破壊強度（$\sigma_c \propto \sqrt{E G_c/l}$）にも影響するため、強度の絶対値はこの $l$ の選択に依存します。

シミュレーション結果と考察

得られた解析結果は、大規模データの処理に優れたオープンソースの可視化ソフトウェア ParaView [7, 8] を用いて可視化・評価を行いました。

1. 荷重–変位応答と合成環境による強度コントラスト

図2は、La-rich・B-rich 両環境の公称応力–変位応答です。いずれも弾性的な立ち上がりののち、明瞭なピークを経て急峻に軟化する、典型的な脆性破壊曲線が得られました。軟化域は、亀裂が進展して荷重負担能力が低下していく過程に対応します。

図2. 合成環境ごとの公称応力 $\langle\sigma_{yy}\rangle$–変位応答。ピーク応力は La-rich で約 29.5 MPa、B-rich で約 49.3 MPa。表面エネルギー（合成環境）の違いが、巨視的な破壊強度の差として明瞭に現れています。

ピーク応力は La-rich で約 29.5 MPa、B-rich で約 49.3 MPa となりました。その比 $49.3/29.5 \approx 1.67$ は、グリフィス則から予想される $\sqrt{G_c^{\,\text{B}}/G_c^{\,\text{La}}} = \sqrt{3} \approx 1.73$ とよく一致します。すなわち、表面エネルギーが約3倍高い B-rich 環境の $\text{LaB}_6$ は、約 1.7 倍高い荷重まで耐えてから破壊します。なお両ピーク応力は、それぞれの均質体強度の目安（La 約 38 MPa、B 約 66 MPa）を下回っており、初期亀裂（切欠き）による応力集中を正しく反映した、物理的に妥当な結果であることを確認しています。

【現在の設定の制約：貫通後の残留応力】
図2では、亀裂貫通後も応力が完全には 0 に戻らず、数 MPa の残留応力が残っています。これは、亀裂進展時（$d \to 1$）の数値的な収束性を確保するため、劣化関数に小さな残留剛性 $\eta = 10^{-3}$ を残している（$g(d) = (1-d)^2(1-\eta) + \eta$）ことに起因する、現在の設定の制約です。$\eta$ を小さくすれば残留応力は 0 に近づきますが、完全分離（最終リガメント破断）近傍での収束が難しくなるトレードオフがあります。本事例の主目的であるピーク強度と靭性コントラストの評価には影響しません。

2. 亀裂進展の可視化

図3・図4は、損傷変数 $d$ の分布の推移（初期 → 進展途中 → 貫通直後）を示しています。応力集中により初期亀裂先端で $d$ が完全破壊状態（$d=1$、赤色表示）に達し、鋭利な亀裂となって右端へ走り抜ける様子が、メッシュに依存せず再現されています。

図3. La-rich（La終端、$\gamma$=0.814 J/m²）の亀裂進展。左：初期、中：進展途中、右：貫通直後。低い表面エネルギーのため、比較的低い荷重で亀裂が貫通します。

図4. B-rich（B₆終端、$\gamma$=2.438 J/m²）の亀裂進展。左：初期、中：進展途中、右：貫通直後。表面エネルギーが高いため、より大きな荷重・変位まで耐えてから同様に貫通します。

両者を比較すると、亀裂進展の様相は定性的に同じ（鋭利な単一亀裂の貫通）でありながら、B-rich の方がより高い荷重・変位に達するまで進展を踏みとどまります。これは、合成環境（化学ポテンシャル）が決める表面エネルギーが、巨視的な耐割れ性を直接支配していることを視覚的に示しています。

【不安定亀裂進展の数値的安定化】
変位制御下での脆性亀裂の進展は、強い非線形性（スナップ的不安定）を伴うため、安定に解くには数値解法上の工夫が必要です。本解析では、(i) 損傷の不可逆性と \(d \le 1\) を課す束縛条件付き（変分不等式）ソルバー、(ii) 破断部の極端な剛性コントラストに対応する直接法線形ソルバー、(iii) 力学場と損傷場の残差スケールを揃える自動スケーリング、(iv) 亀裂進展時に時間刻みを自動細分化する適応タイムステップ制御、を組み合わせています。これにより、ピーク強度や破壊靭性比が時間刻みに依存しない、収束した解を得ています。

まとめ

本事例では、Advance/PHASE による第一原理計算のミクロな物理量（表面エネルギー $\gamma$）を、グリフィス則 $G_c = 2\gamma$ を介して、巨視的な連続体力学（Phase Field法）の破壊パラメータへシームレスに連携させた解析を実施しました。「破壊特性を決める最重要パラメータを1つだけ第一原理計算から抽出し、残りを文献値で補う」という実用的かつ低コストなワークフローにより、以下を定量的に実証しています：(1) 合成環境（化学ポテンシャル → 表面終端 → 表面エネルギー）の違いが、巨視的な劈開破壊強度を約 1.7 倍変化させること。(2) その強度コントラストが、グリフィスの脆性破壊理論（$\sqrt{3}$）と定量的に整合すること。このアプローチは、「製造プロセス条件（組成や合成雰囲気）が、最終製品の信頼性（破壊強度・寿命）をどう左右するか」を予測する材料シミュレーションの枠組みとなります。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献

1. 電子銃陰極材料LaB6の表面エネルギー計算：表面終端と化学ポテンシャル
2. D. Gaston, C. Newman, G. Hansen, and D. Lebrun-Grandie, "MOOSE: A parallel computational framework for coupled systems of nonlinear equations", Nucl. Eng. Des. 239, 1768 (2009).
3. A. A. Griffith, "The phenomena of rupture and flow in solids", Phil. Trans. R. Soc. Lond. A 221, 163 (1921).
4. B. Bourdin, G. A. Francfort, and J.-J. Marigo, "Numerical experiments in revisited brittle fracture", J. Mech. Phys. Solids 48, 797 (2000).
5. C. Miehe, M. Hofacker, and F. Welschinger, "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits", Comput. Methods Appl. Mech. Eng. 199, 2765 (2010).
6. X. Zeng, Y. Ye, S. Zou, Q. Gou, and P. Ou, "First-Principles Study of the Nonlinear Elasticity of Rare-Earth Hexaborides REB₆ (RE = La, Ce)", Crystals 7, 320 (2017).
7. J. Ahrens, B. Geveci, and C. Law, "ParaView: An End-User Tool for Large Data Visualization", The Visualization Handbook, p.717-731, Elsevier (2005).
8. 可視化ソフトウェア ParaView 利用サポートサービス

関連ページ

• 第一原理計算ソフトウェア Advance/PHASE
• 解析分野：ナノ・バイオ
• 産業分野：材料・化学