DFT計算と熱-力学連成Phase Field法によるLaB₆セラミックスの熱衝撃亀裂解析#

PHASE logo

六ホウ化ランタン（LaB₆）は、低い仕事関数と高い融点から電子顕微鏡や電子線描画装置の熱電子放出源（電子銃カソード）として用いられますが、運転時の急加熱・急冷といった熱衝撃による割れが実用上の課題です。脆性セラミックスの耐熱衝撃性は、巨視的な破壊靭性（割れにくさ）に支配され、その破壊靭性は原子スケールで決まる「表面エネルギー」に直結します。さらに LaB₆ では、結晶成長時の合成環境（化学ポテンシャル）に応じて表面終端（La終端／B₆終端）が変化し、表面エネルギー $\gamma$ が約3倍も異なります。本事例では、第一原理計算ソフトウェア Advance/PHASE で算出した LaB₆(001) 面の表面エネルギー $\gamma$ を、グリフィス則（ $G_c = 2\gamma$ ）を介して、マルチフィジックス有限要素法フレームワーク MOOSE 上の熱伝導・熱応力・Phase Field 破壊の連成モデルへ橋渡しします。これにより、合成環境の違い（La過剰／B過剰）が、同一の熱衝撃に対する割れやすさをどれだけ変えるかを検証したシミュレーションの解析結果を紹介します。

Keywords: 第一原理計算（DFT）, Phase Field法, 熱衝撃 (Thermal Shock), 熱応力, 熱-力学連成, LaB₆, 表面エネルギー, グリフィス理論, 耐熱衝撃性, MOOSE

解析モデルと理論背景#

1. マルチスケール・マルチフィジックス連携のアプローチ#

本解析は「応力場」「温度場（熱伝導）」「Phase Field（亀裂場）」の3つの物理場を連成させます。Advance/PHASE による LaB₆ 表面エネルギー $\gamma$ の第一原理計算 [1] では、(001) 面の $\gamma$ が終端構造と化学ポテンシャルに強く依存し、La過剰（La-rich）で安定なLa終端面と B過剰（B-rich）で安定な B₆終端面とで約3倍異なることが示されています。この $\gamma$ を $G_c = 2\gamma$ で破壊靭性に換算し、Phase Field 破壊モデルへ与えます。

巨視モデルは 1.0 mm $\times$ 1.0 mm の2次元平板（平面ひずみ）で、左端中央に長さ 0.2 mm の初期亀裂（損傷場 $d$ の初期条件として導入）を置きます（図1）。一様な高温（800 K、無応力基準温度）の状態から、左端を急冷（クエンチ）して表面に引張熱応力を発生させます。力学境界条件は剛体運動の除去のみ（下端 $y$ 固定、右端 $x$ 固定）とし、熱収縮を一切阻害しない静定支持とします。連成・破壊計算はオープンソースの有限要素法フレームワーク MOOSE [2] 上で実行しました。

初期亀裂モデルの状態
図1. モデルの初期状態。1.0 mm $\times$ 1.0 mm の平板（一様 800 K）の左端に長さ 0.2 mm の初期亀裂（赤色が損傷変数 $d=1$ ）を設定。この左端を急冷して熱衝撃を与えます。※ $d$ は Miehe らのフェーズフィールド破壊での表記、図中の $c$ は MOOSE での入力変数名であり、両者は同一の損傷変数（0＝健全、1＝破断）です。

2. 支配方程式 — 熱から巨視破壊への連成#

熱応力を介して表面エネルギーを Phase Field 破壊へ連成させます [3-5]。

A. 過渡熱伝導と熱ひずみ（温度 → 応力）#

クエンチによる温度場 $T$ の時間発展を過渡熱伝導方程式で解き、温度変化による熱ひずみ（固有ひずみ） $\boldsymbol{\varepsilon}_\text{th}$ を介して熱応力を生じさせます。

$\rho c_p \frac{\partial T}{\partial t} = \nabla\cdot(k\,\nabla T), \qquad \boldsymbol{\varepsilon}_\text{th} = \alpha\,(T - T_0)\,\mathbf{I}$

ここで $\rho$ は密度、 $c_p$ は比熱、 $k$ は熱伝導率、 $\alpha$ は線膨張係数、 $T_0$ は無応力基準温度です。冷えた表面層は収縮しようとしますが内部（高温）に拘束されるため、表面に平行な引張熱応力が生じ、これが亀裂を駆動します。

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

劈開破壊では、臨界エネルギー解放率 $G_c$ は新生する2つの自由表面のエネルギーに等しく、 $G_c = 2\gamma$ で与えられます。表面エネルギー $\gamma$ を代入することで、各合成環境の破壊靭性が一意に定まります。

$G_c = 2\gamma$

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

亀裂を連続的な損傷変数 $d$ （0で健全、1で完全破壊）で表現し、全エネルギー $\Pi$ を最小化するように進展させます。

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

標準的なAT2（Ambrosio-Tortorelli）モデルを採用し、 $\psi_0^{+}$ は（熱ひずみを差し引いた）弾性ひずみエネルギーの引張成分（引張のみが亀裂を駆動）、 $l$ は亀裂幅を規定する長さスケールです。損傷の不可逆性は $d$ の単調増加拘束として課しています。

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

純 LaB₆ を対象とし、Advance/PHASE の第一原理計算値 [1] と文献の熱・弾性物性を組み合わせました。第一原理計算からは破壊靭性を決める「表面エネルギー」のみを抽出し、弾性・熱物性は文献値で補っています。なお、本事例で文献値を用いた弾性定数・線膨張係数・比熱などの物性も、第一原理計算（弾性定数テンソルの直接計算、フォノン計算に基づく熱膨張係数・比熱、状態方程式など）から算出可能であり、必要に応じて当該パラメータを第一原理ベースで与えることもできます。

表1. 熱衝撃亀裂シミュレーションに用いた主なパラメータ

カテゴリー	パラメータ	設定値と出典
第一原理計算 (Advance/PHASE)	表面エネルギー $\gamma$（La終端 / La-rich）	0.814 J/m² [1]
第一原理計算 (Advance/PHASE)	表面エネルギー $\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 = 2\gamma$)	臨界エネルギー解放率 $G_c$（B-rich）	4.876 J/m²（= 4.88 N/m）
熱物性 (文献値)	線膨張係数 ($\alpha$)	2.8 × 10^-6 /K（線膨張係数, 800 K）[6]
	熱伝導率 ($k$)	47 W/(m·K)（代表値）[7,8]
	密度 × 比熱 ($\rho c_p$)	2.27 × 10⁶ J/(m³·K)（$\rho$=4.72 g/cm³ × $c_p$≈480 J/(kg·K)）
弾性・幾何・負荷 (文献値等)	ヤング率 $E$ / ポアソン比 $\nu$	337 GPa / 0.18 [9]
	長さスケール $l$ / 初期亀裂長	0.04 mm / 0.2 mm
	クエンチ温度差 $\Delta T$（左端 800→620 K）	180 K（0.02 s で降下）

長さスケール $l$ と熱物性について： $l$ =0.04 mm は要素寸法 $h\approx0.0083$ mm（1 mm を 120 分割）の約5倍とし、損傷帯を複数要素で解像しています。熱拡散率は $D=k/(\rho c_p)\approx20.7\,\text{mm}^2/\text{s}$ で、1 mm を拡散する時間は $\tau\sim0.05$ s です。なお熱伝導率 $k$ と体積熱容量 $\rho c_p$ は便覧値（ $\rho$ =4.72 g/cm³ は標準値）で、これらは冷却前線の時間スケールを規定するのみで、熱応力の大きさや La/B のコントラストには影響しません。

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

1. 温度場（冷却前線）と熱応力の発生#

図2は、左端クエンチによる温度場 $T$ の推移です（熱物性が同一のため La/B 共通）。急冷された左端から冷却前線が内部へ進行します。前線近傍の冷えた層は収縮しようとして内部に拘束されるため、表面に平行な引張熱応力が生じ、これが初期亀裂の先端付近に集中し、損傷を駆動します。

温度場（冷却前線）の推移
図2. 温度場の推移（左：t=8.6 ms、右：t=20.8 ms）。左端の急冷で生じた冷却前線が内部へ進行する様子。冷却面と内部の温度差（勾配）はクエンチのランプ終了付近（t≈20 ms）で最大となり、このとき最大の引張熱応力が生じます。

2. 合成環境による耐熱衝撃性のコントラスト#

図3は、同一の熱衝撃（ $\Delta T$ =180 K）を与えた La-rich と B-rich の損傷場 $d$ の比較です。結果は明確に分かれました。

La-richとB-richの熱衝撃亀裂コントラスト
図3. 同一熱衝撃（ $\Delta T$ =180 K）に対する損傷場 $d$ （図中のc）の比較。左: La-rich では冷却面（左端）に沿って広い損傷域（中心線で c>0.5 が x≈0.308 mm、さらに冷却面全体に拡がる）を生じる。右: B-rich では損傷が初期亀裂近傍（x≈0.225 mm）にとどまり母材はほぼ無損傷。なお完全亀裂（c≈1）の長さ自体は両者とも初期 0.20 mm 付近で、差は主に準臨界損傷域の広がり（面積）に現れます。

表面エネルギーが低い La-rich（ $G_c$ =1.63 N/m）では、冷却面（左端）に沿って広い損傷域（ $d$ >0.5 の損傷ノード数は B-rich の約8.5倍、冷却面の平均損傷 $d\approx0.65$ ）を生じました。一方、表面エネルギーが約3倍高い B-rich（ $G_c$ =4.88 N/m）では、損傷が初期亀裂近傍（中心線で約 0.22 mm）にとどまり、母材はほぼ無損傷のまま熱衝撃に耐えました。完全亀裂（ $d$ ≈1）の長さ自体は両者とも初期 0.20 mm 付近ですが、AT2 モデルは強度しきい値を持たないため、低 Gc の La-rich では準臨界の微小損傷が広範に蓄積します。すなわち、グリフィス則が示す $\sqrt{3}$ 倍の破壊靭性差が、「同じ熱衝撃で広範に損傷するか・耐えるか」という巨視的な耐熱衝撃性の差として明瞭に現れています。

【数値解法と注意点】
温度（$T\sim10^2\,\text{K}$）・変位（$\sim10^{-3}\,\text{mm}$）・損傷（$\sim1$）の3場はスケールが大きく異なるため、(i) 自動スケーリングで各方程式の残差を釣り合わせ、(ii) 損傷の不可逆性と $d\le1$ を束縛条件付き（変分不等式）ソルバーで課し、(iii) 適応タイムステップ制御で熱過渡と亀裂進展の両時間スケールを追従し、(iv) クエンチをランプ状に与えて初期の応力スパイクを回避、しています。なお採用した AT2 モデルは強度しきい値を持たないため、強く引張を受ける La-rich の冷却面では巨視き裂に加えて準臨界の微小損傷が広がります（現実の表面微小割れに対応）。B-rich 側では母材は健全に保たれます。

まとめ#

本事例では、第一原理計算ソフトウェア Advance/PHASE による表面エネルギー $\gamma$ の DFT 計算値を、グリフィス則 $G_c=2\gamma$ を介して、熱伝導・熱応力・Phase Field 破壊を連成させた巨視モデルへ橋渡ししました。「破壊靭性を決める表面エネルギーだけを第一原理計算から抽出し、熱・弾性物性は文献値で補う」という低コストなワークフローにより、以下を定量的に実証しています：(1) 左端クエンチによる冷却前線と引張熱応力の発生・分布を再現。(2) 同一の熱衝撃（ $\Delta T$ =180 K）に対し、合成環境の違いが「La-rich は広範に損傷／B-rich は耐える」という耐熱衝撃性の差として現れること（損傷領域で約8.5倍の差）。このアプローチは、「製造プロセス条件（合成雰囲気・組成）が、最終製品の耐熱衝撃性・寿命をどう左右するか」を予測し、より割れにくいカソード材料・プロセス設計に直結するシミュレーションの枠組みとなります。

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

お問い合わせ

参考文献#

電子銃陰極材料LaB6の表面エネルギー計算：表面終端と化学ポテンシャル
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).
A. A. Griffith, "The phenomena of rupture and flow in solids", Phil. Trans. R. Soc. Lond. A 221, 163 (1921).
B. Bourdin, G. A. Francfort, and J.-J. Marigo, "Numerical experiments in revisited brittle fracture", J. Mech. Phys. Solids 48, 797 (2000).
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).
E. K. Akdoğan and E. A. Payzant, "Thermal expansion of LaB₆ from 298 to 998 K", J. Appl. Crystallogr. 58, 1635 (2025).
P. A. Popov, V. V. Novikov, A. A. Sidorov and E. V. Maksimenko, "Thermal conductivity of LaB₆ and SmB₆ in the range 6-300 K", Inorg. Mater. 43, 1187 (2007).
T. Tanaka, "The thermal and electrical conductivities of LaB₆ at high temperatures", J. Phys. C: Solid State Phys. 7 (1974).
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).

カテゴリー	パラメータ	設定値と出典
第一原理計算 (Advance/PHASE)	表面エネルギー \(\gamma\)（La終端 / La-rich）	0.814 J/m² [1]
第一原理計算 (Advance/PHASE)	表面エネルギー \(\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 = 2\gamma\))	臨界エネルギー解放率 \(G_c\)（B-rich）	4.876 J/m²（= 4.88 N/m）
熱物性 (文献値)	線膨張係数 (\(\alpha\))	2.8 × 10^-6 /K（線膨張係数, 800 K）[6]
	熱伝導率 (\(k\))	47 W/(m·K)（代表値）[7,8]
	密度 × 比熱 (\(\rho c_p\))	2.27 × 10⁶ J/(m³·K)（\(\rho\)=4.72 g/cm³ × \(c_p\)≈480 J/(kg·K)）
弾性・幾何・負荷 (文献値等)	ヤング率 \(E\) / ポアソン比 \(\nu\)	337 GPa / 0.18 [9]
	長さスケール \(l\) / 初期亀裂長	0.04 mm / 0.2 mm
	クエンチ温度差 \(\Delta T\)（左端 800→620 K）	180 K（0.02 s で降下）

DFT計算と熱-力学連成Phase Field法によるLaB6セラミックスの熱衝撃亀裂解析#