Cahn–Hilliard 方程式で記述する金属組織シミュレーション

Cahn–Hilliard（CH）法は「保存量である濃度（組成）場」の時間発展を、自由エネルギー低下と拡散の整合で記述する枠組みである。相分離（スピノーダル分解）や析出・粗大化（オストワルド成長）を、界面エネルギーを含む連続体場として同じ数式骨格で扱えるのが特徴である。

参考ドキュメント

1. J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, Journal of Chemical Physics 28 (1958) 258. https://pubs.aip.org/aip/jcp/article/28/2/258/74794/Free-Energy-of-a-Nonuniform-System-I-Interfacial
2. Phase Field Method Recommended Practices Guide（Numerical Implementation / Problem Set-Up）, NIST https://pages.nist.gov/pf-recommended-practices/bp-guide-gh/ch2-numerical-implementation.html
3. 小山敏幸, Phase-field法に関する最近の進展と今後の展望, まてりあ 42 (2003) 397. https://www.jstage.jst.go.jp/article/materia1994/42/5/42_5_397/_pdf

1. 何を Cahn–Hilliard で解くのか

CH法が直接扱う主変数は、保存則に従う場である。

• 二元合金の濃度：$c(\mathbf{x},t)$（例：溶質原子分率）
• 多元合金の濃度ベクトル：$\mathbf{c}=(c_1,c_2,\dots )$
• 点欠陥（空孔など）の濃度：$c_v(\mathbf{x},t)$（拡散と相互作用のモデル化が可能）
• 規則度や相指標などの非保存量は、通常 Allen–Cahn（AC）に担わせ、CHは濃度の保存発展に集中させる

整理すると次の通りである。

項目	Cahn–Hilliard（CH）	Allen–Cahn（AC）
主変数	濃度など保存量 $c$	秩序変数など非保存量 $\varphi$
基本構造	連続の式（フラックス発散）	勾配流（緩和）
典型現象	相分離、析出・粗大化、濃度拡散	粒界移動、規則化、相の切替

金属組織で多い構成は「相の指標（AC）＋濃度（CH）」の連成であり、溶質分配と相移動を同居させる。

2. 自由エネルギー汎関数

CH法は自由エネルギー汎関数 $F[c]$ を起点とする。最も基本的な形は

$$F[c]={\int }_{\mathrm{\Omega }}(f_{\mathrm{chem}}(c,T)+\frac{\kappa }{2}|\mathrm{\nabla }c{|}^2)dV$$

である。

• $f_{\mathrm{chem}}(c,T)$：局所（化学）自由エネルギー密度
  • 正則溶体（regular solution）近似、あるいはCALPHADの自由エネルギーを用いる設計が多い
• $\kappa$：勾配エネルギー係数
  • 濃度勾配に罰則を与え、結果として有限幅の界面と界面エネルギーを生む
  • 解析的には「界面を持つ非一様系の自由エネルギー」の理論に対応する

界面厚み（拡散界面幅）$\delta$ と界面エネルギー $\sigma$ は、$f_{\mathrm{chem}}$ と $\kappa$ の組により決まり、ここがメッシュ設計や定量性に直結する。

3. Cahn–Hilliard 方程式の導出：保存則＋化学ポテンシャル

本質は「質量保存」と「拡散フラックスが化学ポテンシャル勾配で駆動される」という2点である。

• 保存則（連続の式）

$$\frac{\partial c}{\partial t}=-\mathrm{\nabla }\cdot \mathbf{J}$$

• 拡散フラックス（オンザガー型）

$$\mathbf{J}=-M(c,T)\mathrm{\nabla }\mu$$

• 化学ポテンシャル（汎関数微分）

$$\mu =\frac{\delta F}{\delta c}=\frac{\partial f_{\mathrm{chem}}}{\partial c}-\kappa {\mathrm{\nabla }}^{2}c$$

以上より

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }[\frac{\partial f_{\mathrm{chem}}}{\partial c}-\kappa {\mathrm{\nabla }}^{2}c]\right)$$

が CH 方程式である。$M$ と $\kappa$ が定数なら

$$\frac{\partial c}{\partial t}=M{\mathrm{\nabla }}^2\left(\frac{\partial f_{\mathrm{chem}}}{\partial c}\right)-M\kappa {\mathrm{\nabla }}^{4}c$$

となり、4階空間微分を含むのが数値的な特徴である。

3.1 物理的に重要な2つの保存・散逸

適切な境界条件（周期境界、あるいは無流束 $\mathbf{n}\cdot \mathbf{J}=0$）の下で次が成り立つ。

• 質量保存：

$$\frac{d}{dt}{\int }_{\mathrm{\Omega }}cdV=0$$

• 自由エネルギー減少（散逸）：

$$\frac{dF}{dt}=-{\int }_{\mathrm{\Omega }}M|\mathrm{\nabla }\mu {|}^{2}dV\le 0$$

CH法に基づくシミュレーションでは、これらの量をモニタして数値解が「それらしい」だけでなく「条件を満たす」ことを確認するのが基本である。

4. 金属組織への対応

4.1 スピノーダル分解（相分離初期）と線形安定解析

二元系で一様濃度 $c_0$ の近傍を考える。$f_{\mathrm{chem}}^{″}(c_0)<0$ の領域では微小ゆらぎが成長し、相分離が始まる。

線形化すると、フーリエ成分（波数 $k$）の成長率は概念的に

$$R(k)\propto -M{k}^2[f_{\mathrm{chem}}^{″}(c_0)+\kappa k^2]$$

となり、

• 不安定帯域：$0<k<\sqrt{-f^{″}/\kappa }$
• 最も成長しやすい波数：$k_m=\sqrt{-f^{″}/(2\kappa )}$
• 代表長さ：${\lambda }_m=2\pi /k_m$ が得られる。

この議論は「初期パターンの波長」と「$\kappa$ の役割」を定性的に理解する上で有効である。

4.2 析出・粗大化（オストワルド成長）

CH法は、界面エネルギーを含む拡散駆動で析出相の成長・消滅を記述できる。

• 初期：核生成をどう扱うかが設計点
  • 初期条件に微小ノイズを入れる
  • あるいは古典核生成理論（別モデル）で核を配置してCHで成長・粗大化へ移る
• 後期：粗大化（coarsening）
  • 典型的には代表長さが時間に対してべき則で成長することが多い
  • ただし、弾性拘束や拡散係数の濃度依存、体積分率によってスケーリングは変化する

4.3 多元合金：濃度ベクトルと移動度行列

多元系では濃度 $\mathbf{c}$ と化学ポテンシャル $\mathit{\mu }$ をベクトル化し、移動度を行列 $M_{ij}$ として

$$\frac{\partial c_i}{\partial t}=\mathrm{\nabla }\cdot \left(\sum _j{M}_{ij}\mathrm{\nabla }{\mu }_j\right),{\mu }_i=\frac{\delta F}{\delta c_i}$$

と書く。実装上は

• $M_{ij}$ を対称・正定にする（散逸則の維持）
• 独立変数数の削減（$\sum _i{c}_i=1$ の制約） などの工夫が必要になる。

4.4 弾性との連成：コヒーレント析出・格子ミスフィット

コヒーレント析出では弾性エネルギーが形態選択（板状・針状）や配列（ラメラ・整列）を支配しうる。自由エネルギーに

$$F_{\mathrm{el}}={\int }_{\mathrm{\Omega }}{f}_{\mathrm{el}}(\epsilon ,{\epsilon }^{\ast }(c),\mathbf{C})dV$$

を加えると、化学ポテンシャルに $\delta F_{\mathrm{el}}/\delta c$ が加算され、濃度拡散が応力場により偏る。

• ${\epsilon }^{\ast }(c)$：固溶・析出に伴う固有ひずみ（Vegard則近似など）
• $\mathbf{C}$：弾性定数テンソル
• 連成が強いと剛性が増し、陰的解法や反復解法が必要になりやすい

4.5 ノイズを含む拡張：Cahn–Hilliard–Cook（CHC）

核生成や熱揺らぎを扱う目的で、確率フラックスを加えた形式が議論される。

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )+\mathrm{\nabla }\cdot \mathit{\xi }$$

ここで $\mathit{\xi }$ はゆらぎを表す確率過程であり、離散化においては相関構造や保存則との整合が注意点となる。

5. モデルパラメータの設計：$f_{\mathrm{chem}}$・$\kappa$・$M$

5.1 化学自由エネルギー $f_{\mathrm{chem}}$

• 正則溶体モデル（理解用）$$f_{\mathrm{chem}}(c,T)=\mathrm{\Omega }c(1-c)+k_{B}T[c\ln c+(1-c)\ln (1-c)]$$$\mathrm{\Omega }$ が混合エンタルピー、対数項が配置エントロピーである
• CALPHAD接続（定量志向）
  • 実系に近い自由エネルギー曲面を使えるが、相領域の扱い（多相・多成分の切替）をどうするかが設計点となる

5.2 勾配係数 $\kappa$ と界面エネルギー

$\kappa$ は界面幅と界面エネルギーに効く。数値的には

• 解像できる界面幅を先に決め、そこから $\kappa$ を調整する設計
• 実験的界面エネルギー $\sigma$ を目標にし、$\kappa$ と $f_{\mathrm{chem}}$ の組で合わせる設計 がある。

5.3 移動度 $M(c,T)$ と時間スケール

二元系では拡散係数 $D$ を用い

$$M(c,T)\sim \frac{D(c,T)}{k_{B}T}c(1-c)$$

のような形がよく現れる（因子はモデル化の選び方による）。金属では拡散係数が温度依存（アレニウス）かつ相依存であり、相場（AC）と結合する場合は「相ごとに $D$ や $M$ を切り替える」設計も出てくる。

6. 数値解法：4階PDEをどう解くか

CH方程式は4階であり、安定性の観点から陽解法は刻み制約が厳しいことが多い。代表的な方針は次である。

6.1 分割（mixed / split）で 2階×2 に落とす

化学ポテンシャル $\mu$ を未知数として導入し、

$$\mu -\frac{\partial f_{\mathrm{chem}}}{\partial c}+\kappa {\mathrm{\nabla }}^{2}c=0,\frac{\partial c}{\partial t}-\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )=0$$

の連立（どちらも2階）として解く。FEMでは標準的な $C^0$ 要素で扱える利点が大きい。

6.2 半陰解法（IMEX）とエネルギー安定化

• 拡散（剛な線形項）を陰的に、非線形項を陽的に扱う半陰は定番である
• 自由エネルギー単調減少を離散でも守る設計（凸分割、安定化付き半陰、補助変数法など）は長時間発展で効きやすい

6.3 スペクトル法（FFT）

周期境界かつ矩形領域では、フーリエ空間で線形部を代数的に処理でき、高速である。相分離ベンチマーク（スピノーダル分解）では、FFT＋半陰がよく利用される。

6.4 検証

CHコードの健全性確認として、少なくとも次を追う。

• $\int cdV$ の保存（時間方向にほぼ一定）
• $F(t)$ の非増加（あるいは時間刻みに応じた一貫した減少）
• メッシュ・時間刻みの収束試験（パターン統計量が変わらないか）

7. 典型的な問題設定

7.1 初期条件

• 一様濃度 $c_0$ ＋微小ノイズ：スピノーダル分解の典型
• 既に核（析出相）を埋め込む：粗大化・相互作用評価を狙う場合

7.2 境界条件

• 周期境界：統計量評価やFFTに適する
• 無流束：$\mathbf{n}\cdot \mathbf{J}=0$（閉じた系で質量保存を貫く）
• 表面偏析や外部フラックスを扱う場合：境界での化学ポテンシャル条件や表面自由エネルギー項を追加する設計となる

7.3 観測量（後処理）

• 体積分率、相の連結性、界面面積密度
• 濃度ヒストグラム $p(c)$（二峰性の発達が相分離の進行を示す）
• 構造因子 $S(\mathbf{k})$（代表波長や粗大化の追跡に有効）
• 代表長さ $L(t)$（等濃度面からの曲率統計、または $S(k)$ のピークから推定）

8. 公開実装・ベンチマーク

• PFHub Benchmark Problem 1（Spinodal Decomposition） 相分離の初期〜粗大化を比較可能な形で規定しており、CH実装の検証に使いやすい
• MOOSE Phase Field の CahnHilliard カーネル 自由エネルギーの微分と移動度を与えてCHのFEM実装を組み立てる設計である
• PRISMS-PF 金属の相分離・析出・粒成長などのマイクロ組織進化を大規模並列で扱うための枠組みとして整備されている

（上記は概念整理のための入口であり、モデルの選択は自由エネルギーと移動度の具体形に強く依存する。）

9. よくある落とし穴

• 界面幅を変えると結果が変わる 定量志向では「どの量を不変に保つか（σ、粗大化速度、分配係数など）」を定めて再同定が必要である
• $M(c,T)$ の選び方で時間スケールが大きく変わる 相依存・濃度依存が落ちると、形態が似ていても速度論が外れることがある
• 4階PDEを強引に扱うと数値不安定になる 分割形（2階×2）や陰的・半陰、前処理により成立性が上がる
• 連成（弾性、ACとの結合）で剛性が増す 収束条件や時間刻みを系統的に評価しないと、見かけ上の組織が数値誤差に引きずられる

まとめ

Cahn–Hilliard 法は、濃度保存・界面エネルギー・拡散駆動を一体化して金属組織の相分離や析出・粗大化を記述する基礎式である。自由エネルギー汎関数の設計（化学自由エネルギー、勾配係数、弾性項）と、保存則・散逸則を壊しにくい数値解法（分割形、半陰・陰、検証）が、結果の信頼性を左右する。