フェーズフィールド法の基礎と数値計算

自由エネルギー汎関数に基づく連続体モデリングと数値計算の体系

本書は、界面を拡散化した連続場によって組織形成やドメイン形成を記述するフェーズフィールド法を、熱力学・変分法・数値解法の順に体系化する構成である。材料組織、破壊、電気化学、強誘電・磁性などの多様な対象に対して、同一の数理骨格からモデルを立ち上げ、検証可能な計算へ落とし込むための視点を与えるものである。

目次

• 第1章 フェーズフィールド法の位置づけ
• 第2章 熱力学と自由エネルギー汎関数
• 第3章 変分法と支配方程式
• 第4章 界面、勾配エネルギー、異方性
• 第5章 多相・多成分・化学ポテンシャルの取り扱い
• 第6章 弾性・塑性・欠陥場との連成
• 第7章 熱・流体・凝固と成長
• 第8章 粒成長・再結晶・結晶方位場
• 第9章 破壊・損傷・界面剥離のフェーズフィールド
• 第10章 電気化学・電池・腐食のフェーズフィールド
• 第11章 強誘電・磁性・秩序変数場
• 第12章 数値解法と離散化設計
• 第13章 検証、妥当性、再現性
• 第14章 計算基盤と高速化
• おわりに

第1章 フェーズフィールド法の位置づけ

1.1 拡散界面という考え方

相の境界は、教科書的には鋭い幾何学的境界として描かれることが多い。しかし実際の材料では、原子配列や組成は界面近傍で連続的に変化し、有限幅の遷移層を持つ場合が多い。フェーズフィールド法は、この遷移層を明示的に取り込み、界面を有限幅の連続場として表す方法である。

界面を表す代表的な変数は秩序変数と呼ばれ、例えば固相を1、液相を0とするように、相の状態を連続的な値で表す。ここで重要なのは、秩序変数は観測装置で直接測る物理量というより、相の状態を表現するために導入する連続体の内部変数である点である。したがって、秩序変数の定義と自由エネルギーの設計が、そのまま界面の性質や界面移動の法則を決める。

拡散界面の利点は、界面の位置を追跡する別アルゴリズムが不要になることである。界面の生成・消滅、分岐、合体が起きても、秩序変数の空間分布が自動的にそれを表すため、トポロジー変化を自然に扱える。これは、粒界ネットワーク、デンドライト分岐、亀裂分岐など、幾何学的な追跡が難しい現象で特に効く。

一方で、界面が有限幅であることは、数値計算ではメッシュで界面幅を解像する必要があることを意味する。界面幅を極端に薄くすると鋭い界面へ近づくが、メッシュが極端に細かく必要となり計算コストが増える。界面幅を厚くすると計算は楽になるが、界面の曲率効果や界面反応の表現に影響が出る。したがって、モデルで表現したい物理と計算規模の間で、界面幅とパラメータの整合を取ることが出発点になる。

1.2 保存場と非保存場

フェーズフィールド法で扱う場は、大きく保存場と非保存場に分けられる。保存場とは、領域全体で積分した量が時間発展で変わらない（あるいは境界からの流入出でのみ変化する）場である。代表例は組成場であり、濃度 $c(\mathbf{r},t)$ が拡散で移動しても、閉じた系の総量 ${\int }_{\mathrm{\Omega }}cd\mathbf{r}$ は一定になる。これは物質保存則に対応する。

一方で、非保存場とは、局所的に値が増減し得る場である。秩序変数 $\eta (\mathbf{r},t)$ はその代表であり、相の状態が変化すること自体が局所的な生成・消滅として表れる。非保存場は、例えば固液変態で固相分率が局所的に増えるとき、別の場所で同量減る必要はないという直感に対応する。

この区別は、支配方程式の形を根本から分ける。非保存場は、自由エネルギーを下げる方向に局所的に緩和する形になりやすく、時間発展は勾配流として

$$\frac{\partial \eta }{\partial t}\propto -\frac{\delta F}{\delta \eta }$$

の形を取る。保存場は、保存則から連続の式

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

がまず立ち、フラックス $\mathbf{J}$ を化学ポテンシャル勾配で与えることで

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\mu \right),\mu =\frac{\delta F}{\delta c}$$

へ至る。つまり、非保存場では局所緩和、保存場では輸送（拡散）が本質である。

保存か非保存かの選択は、材料の現象論に対応する。例えば、組成による相分離は保存場の拡散で進み、秩序・無秩序転移は秩序変数の緩和で進む。実際の問題では、両者が同時に起きる場合が多く、秩序変数と組成を同時に持つ連成モデルが基本となる。

1.3 適用スケールと他手法との関係

フェーズフィールド法は、原子を個別に追う方法ではなく、連続場でメゾスケールの形態進化を扱う方法である。典型的な空間スケールは、界面幅より十分大きく、組織スケール（粒径、析出物間隔、ドメイン幅）を直接表せる範囲に置かれる。時間スケールも、原子振動や衝突の時間ではなく、拡散や界面移動、緩和の時間を主に扱う。

他手法との関係は、入力と検証の両面で整理できる。第一原理計算は、相安定性、格子ミスマッチ、弾性定数、界面エネルギーの推定など、自由エネルギーや連成項の材料定数に関与する。分子動力学は、界面モビリティ、拡散係数、核生成に関する局所過程の情報を与え得る。CALPHADは、多成分系の熱力学データベースとして $f_{\mathrm{loc}}(c,T)$ の実装に結びつく。これらはフェーズフィールドの外側にありつつ、自由エネルギー設計を通じて内部へ入ってくる。

鋭い界面を前提にした古典的な界面追跡法（ステファン問題など）との関係も重要である。フェーズフィールドが鋭い界面の極限で同じ界面条件へ一致するためには、界面幅を小さくするだけでなく、パラメータの取り方が条件を満たす必要がある。どの極限でどの界面条件を再現したいのかを明確にしないと、同じ見かけの形態でも速度や選択則が変わり得る。

以下は、手法の位置づけを整理するための比較である。

手法	記述の単位	得意な現象	制約の中心
第一原理計算	電子状態	相安定性、界面・欠陥の基礎量	系サイズが小さい
分子動力学	原子	界面の局所過程、拡散の素過程	時間が短い
フェーズフィールド	連続場	組織・ドメイン・界面の形態進化	材料定数の同定が要る
鋭い界面モデル	境界面	界面条件が既知な成長・移動	トポロジー変化が難しい

この表は優劣を示すものではなく、スケールと記述の単位の違いを示すための整理である。フェーズフィールドは、複雑形態の時間発展を扱える代わりに、自由エネルギーと運動係数という入力設計が決定的に重要となる。

1.4 研究史とコミュニティ

フェーズフィールド法の理論的基盤は、界面の自由エネルギーを勾配項として書き下し、変分原理と散逸の仮定から時間発展を構成する点にある。この枠組みは材料の凝固、相分離、析出、粒成長、ドメイン形成などへ拡張され、レビュー論文で統一的に整理されてきた。英語圏では、モデルの設計思想と適用範囲をまとめた総説が参照される。

近年は、同じ現象を別実装で計算した際に結果が一致するかを比較しやすくするため、ベンチマークや推奨事項が整備されている。ここで扱われるのは、単なるプログラム上の取り決めではなく、無次元化、界面幅の解像、時間刻み、初期条件など、物理と離散化が交差する要点である。計算が現象の本質を写しているのか、離散化の癖を見ているのかを切り分けるための共通言語として機能する。

国内でも、凝固や組織形成を中心に講義資料・解説記事・入門書が整備されてきた。国内資料は、材料学の文脈（相図、凝固組織、結晶粒）と数理の橋渡しを意識したものが多く、初学者がモデルの意味づけを掴む助けになる。フェーズフィールド法は抽象度が高いが、国内外の整理を参照しつつ、自分の対象に合わせて秩序変数と自由エネルギーを設計する姿勢が重要である。

1.5 本記事としての到達点

本記事で到達したいのは、与えられた方程式を使うだけでなく、自由エネルギーから支配方程式を自分で導出できる状態である。導出は、(i) どの場が保存でどの場が非保存かを判断し、(ii) 自由エネルギー汎関数 $F$ を書き、(iii) 機能微分 $\delta F/\delta \varphi$ を計算し、(iv) 散逸則と保存則を組み合わせて時間発展式を作る、という順序で進む。この順序を崩すと、式の各項の意味が曖昧になりやすい。

また、数式を導いた後に必ず確認すべきことがある。第一に、単位整合である。$f_{\mathrm{loc}}$、勾配係数、運動係数が整合していないと、計算時間や界面速度が物理と対応しなくなる。第二に、エネルギーが時間とともに減少する構造になっているかである。多くのフェーズフィールド方程式は、適切な係数を用いれば $dF/dt\le 0$ を満たすように作れるため、その確認はモデルの健全性の診断になる。

さらに、対象とする現象がどのスケールの情報を要求しているかを明確にし、必要な材料定数がどこから来るかを見通しておく必要がある。フェーズフィールド法は、形態の時間発展を与えるが、材料定数の供給源は外部にあることが多い。どの定数を既知とし、どれを同定対象とするかを最初に整理することで、研究としての見通しが立ちやすくなる。

第2章 熱力学と自由エネルギー汎関数

2.1 局所自由エネルギー密度

フェーズフィールド法の出発点は、状態を特徴づける自由エネルギー密度を定めることである。最も基本的には、局所自由エネルギー密度 $f_{\mathrm{loc}}(\{{\eta }_i\},c,T)$ を導入し、秩序変数と組成と温度に依存する関数として与える。ここで $\{{\eta }_i\}$ は複数の相や結晶方位を表す秩序変数の集合であり、$c$ は保存される組成場である。

局所自由エネルギー密度は、相の安定性を表現するために、多井戸ポテンシャルの形を取ることが多い。例えば二相系で秩序変数 $\eta$ を用いる場合、$\eta =0$ と $\eta =1$ に安定点（極小）を持つ形を作ると、相の共存が表現できる。加えて、組成依存性を入れると、組成に応じてどちらの相が安定か、どの程度の駆動力があるかが表現される。

多成分合金や温度依存性を正確に扱う場合、$f_{\mathrm{loc}}$ を相図に整合する形で与えることが重要になる。CALPHAD由来のギブズ自由エネルギー $G(T,c)$ を局所項として用いると、熱力学データベースと接続できる。初学者の段階では、まず多井戸ポテンシャルで相安定性を設計し、その後、相図整合のための $G(T,c)$ へ拡張する流れが理解しやすい。

2.2 自由エネルギー汎関数の基本形

局所自由エネルギー密度だけでは、界面のエネルギーが表現できない。界面は空間的に変化する領域であり、その変化の急峻さにエネルギーがかかると考えると自然である。これを表すため、勾配エネルギー項を加えた自由エネルギー汎関数を定義する。

基本形は、領域 $\mathrm{\Omega }$ に対して

$$F[\{{\eta }_i\},c]={\int }_{\mathrm{\Omega }}(f_{\mathrm{loc}}(\{{\eta }_i\},c,T)+\sum _i\frac{{\kappa }_i}{2}|\mathrm{\nabla }{\eta }_i{|}^2+\frac{{\kappa }_c}{2}|\mathrm{\nabla }c{|}^2+f_{\mathrm{el}}(\{{\eta }_i\},c,\mathbf{u})+f_{\mathrm{ext}})d\mathbf{r}$$

の形で書かれる。${\kappa }_i$ は秩序変数の勾配係数であり、${\kappa }_c$ は組成勾配係数である。$f_{\mathrm{el}}$ は弾性エネルギー（変位 $\mathbf{u}$ を含む）であり、$f_{\mathrm{ext}}$ は外場（電場、磁場、化学外場など）との結合を表す。

この汎関数の重要点は、どの項もエネルギー密度であり、積分するとエネルギーになることである。したがって、各項の単位が一致していなければならない。例えば $f_{\mathrm{loc}}$ が J/m${}^3$ の単位を持つなら、勾配項も同じ単位になるように ${\kappa }_i$ の単位が決まる。弾性項や外場項も同じ基準で整合させる必要がある。

また、勾配項は界面の有限幅を生む。局所項だけなら相は鋭く切り替わってよいが、勾配項があると空間変化が急峻になるほどエネルギーが増えるため、界面はある幅に広がる。この界面幅が数値計算では解像すべき長さスケールとなり、モデル設計と離散化設計を結びつける。

2.3 散逸とエネルギー減少

フェーズフィールド方程式は、自由エネルギーを下げる方向に時間発展するように作ることが多い。この性質は、物理的には不可逆過程の散逸に対応し、数学的には勾配流の構造として表現できる。

非保存場 $\eta$ の場合、Allen–Cahn型の時間発展を

$$\frac{\partial \eta }{\partial t}=-L\frac{\delta F}{\delta \eta }$$

と置く。$L>0$ は運動係数である。このとき、汎関数の時間微分は

$$\frac{dF}{dt}={\int }_{\mathrm{\Omega }}\frac{\delta F}{\delta \eta }\frac{\partial \eta }{\partial t}d\mathbf{r}=-{\int }_{\mathrm{\Omega }}L{\left(\frac{\delta F}{\delta \eta }\right)}^{2}d\mathbf{r}\le 0$$

となり、自由エネルギーが単調に減少する。ここでは境界項が消えるような境界条件（例えば勾配が境界で0など）を想定している。重要なのは、運動係数が正であれば、方程式の構造だけで $dF/dt\le 0$ が保証される点である。

保存場 $c$ の場合、連続の式から

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

であり、フラックスを

$$\mathbf{J}=-M\mathrm{\nabla }\mu ,\mu =\frac{\delta F}{\delta c}$$

と仮定する。すると

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\mu \right)$$

となる。エネルギー減少は

$$\frac{dF}{dt}={\int }_{\mathrm{\Omega }}\frac{\delta F}{\delta c}\frac{\partial c}{\partial t}d\mathbf{r}={\int }_{\mathrm{\Omega }}\mu \mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )d\mathbf{r}=-{\int }_{\mathrm{\Omega }}M|\mathrm{\nabla }\mu {|}^{2}d\mathbf{r}\le 0$$

となり、$M>0$ であれば減少が示される。途中で部分積分を用い、境界項が消える条件（例えば $\mathbf{n}\cdot \mathrm{\nabla }\mu =0$）を仮定した。この導出により、Cahn–Hilliard方程式が自由エネルギーを減らす拡散として理解できる。

このように、散逸の仮定は、方程式を恣意的に選ぶのではなく、熱力学的整合（エネルギー減少）を満たす形に固定する役割を持つ。連成が増えるほど、どの変数がどの形で散逸するかを明示することが、モデルの理解と拡張の両方に効く。

2.4 無次元化と単位整合

フェーズフィールド法では、界面幅や拡散係数など複数の長さ・時間スケールが登場し、式の係数が桁違いになりやすい。無次元化は、代表スケールを選び、変数と式をスケール化して数値的に扱いやすい形へ変換する操作である。無次元化は単なる書き換えではなく、支配的な項と無視できる項を見通す道具にもなる。

例えば長さスケールを ${\ell }_0$、エネルギー密度スケールを $f_0$、時間スケールを $t_0$ と置き、

$$\tilde{\mathbf{r}}=\mathbf{r}/{\ell }_0,\tilde{t}=t/t_0,\tilde{F}=F/(f_0{\ell }_0^d)$$

のように無次元量を定義する（$d$ は次元である）。このとき勾配は $\mathrm{\nabla }=(1/{\ell }_0)\tilde{\mathrm{\nabla }}$ となるため、勾配項の係数は $\kappa /(f_0{\ell }_0^2)$ のような無次元比として現れる。ここで界面幅を ${\ell }_0$ と同程度に取ると、界面近傍の変化が無次元空間で $O(1)$ に揃い、離散化設計が考えやすくなる。

単位整合は、無次元化の前提条件である。$f_{\mathrm{loc}}$ が J/m${}^3$ なら、$\kappa |\mathrm{\nabla }\eta {|}^2$ も J/m${}^3$ でなければならず、$\kappa$ の単位は J/m となる。さらに時間発展式で ${\partial }_t\eta$ に現れる係数 $L$ は、$\eta$ が無次元なら $L$ の単位が (m${}^3$)/(J·s) のように決まる。保存場では $M$ の単位が拡散係数と化学ポテンシャルの単位の組で決まる。ここを曖昧にすると、計算で得られた時間が現実の秒へ対応しなくなる。

連成項（弾性、電場、磁場など）を加えると、エネルギー密度の単位整合がさらに重要になる。弾性エネルギー密度は Pa = J/m${}^3$ と一致するため整合しやすいが、外場項は定義次第で係数の単位が変わる。したがって、自由エネルギーをすべて同じ単位系へ統一し、どの項がどのスケールを支配するかを無次元比として整理することが、安定したモデル構築につながる。

2.5 乱雑さの導入

核生成や熱揺らぎの影響を表現したい場合、確率的な項を導入することがある。例えば非保存場では

$$\frac{\partial \eta }{\partial t}=-L\frac{\delta F}{\delta \eta }+\xi (\mathbf{r},t)$$

のように雑音項 $\xi$ を加える。保存場ではフラックスに雑音を加える形で保存則を保つ構成が用いられることが多い。どちらも、場の揺らぎが成長して相分離や核生成を誘起する状況を表すための操作である。

ここで注意すべき点は、雑音を入れることは物理仮定の追加であり、同時に数値仮定の追加でもある点である。連続体で白色雑音を仮定しても、離散化すると格子幅と時間刻みが揺らぎ強度に関係する。したがって、雑音強度を選ぶときには、どのスケールの揺らぎを表現したいのかを明確にし、離散化の変更で揺らぎが不自然に変わらないように定義を与える必要がある。

初学者の段階では、まず雑音なしの決定論的モデルで、自由エネルギーと駆動力が作る形態進化を理解することが重要である。その上で、核生成を扱いたい、微小な初期ゆらぎからの発展を扱いたい、といった目的が出たときに、雑音を熱力学（揺動散逸の考え方）と整合する形で追加するのが自然である。

第3章 変分法と支配方程式

3.1 機能微分とオイラー・ラグランジュ方程式

フェーズフィールド法の支配方程式は、自由エネルギー汎関数 $F[\varphi ]$ の機能微分 $\delta F/\delta \varphi$ を計算することから得られる。機能微分とは、関数 $\varphi (\mathbf{r})$ を微小に変化させたときの $F$ の変化を、線形項として書き下したときに現れる密度である。通常の微分が変数の微小変化に対する関数値の変化を与えるのと同様に、機能微分は場の微小変化に対する汎関数の変化を与える。

定義として、$\varphi \to \varphi +ϵ\psi$ と変化させたとき

$${\frac{d}{dϵ}F[\varphi +ϵ\psi ]|}_{ϵ=0}={\int }_{\mathrm{\Omega }}\frac{\delta F}{\delta \varphi }\psi d\mathbf{r}$$

を満たす量が $\delta F/\delta \varphi$ である。ここで $\psi$ は任意の試験関数である。この形に落とすためには、部分積分が必須になる。特に勾配項は $\mathrm{\nabla }\varphi$ を含むため、$\psi$ にかかる微分を $\varphi$ 側に移す操作を行い、境界項を境界条件で消す。

基本例として

$$F[\varphi ]={\int }_{\mathrm{\Omega }}(f(\varphi )+\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2)d\mathbf{r}$$

を考える。変分を取ると

$$\delta F={\int }_{\mathrm{\Omega }}(\frac{\partial f}{\partial \varphi }\delta \varphi +\kappa \mathrm{\nabla }\varphi \cdot \mathrm{\nabla }(\delta \varphi ))d\mathbf{r}$$

となる。第2項を部分積分すると

$${\int }_{\mathrm{\Omega }}\kappa \mathrm{\nabla }\varphi \cdot \mathrm{\nabla }(\delta \varphi )d\mathbf{r}=-{\int }_{\mathrm{\Omega }}\kappa ({\mathrm{\nabla }}^2\varphi )\delta \varphi d\mathbf{r}+{\int }_{\partial \mathrm{\Omega }}\kappa (\mathbf{n}\cdot \mathrm{\nabla }\varphi )\delta \varphi dS$$

である。境界で $\mathbf{n}\cdot \mathrm{\nabla }\varphi =0$ あるいは $\delta \varphi =0$ を課すと境界項が消え、結局

$$\frac{\delta F}{\delta \varphi }=\frac{\partial f}{\partial \varphi }-\kappa {\mathrm{\nabla }}^2\varphi$$

を得る。フェーズフィールド方程式に現れるラプラシアンは、この変分計算から自然に出るものである。

この計算手順を理解すると、多変数（${\eta }_i,c$）や連成項（弾性、外場）でも、どの項がどの微分演算子を生むかを体系的に追えるようになる。初学者が最初につまずきやすいのは、部分積分と境界項の扱いであるが、境界条件は単なる数学的条件ではなく、系が外部とどう接しているかを表す物理条件である点を意識すると理解が進む。

3.2 Allen–Cahn 方程式

Allen–Cahn方程式は非保存場の時間発展を表す最も基本的な式である。構成の考え方は、秩序変数が自由エネルギーを下げる方向へ緩和すると仮定し、その速度が熱力学的駆動力 $\delta F/\delta \eta$ に比例すると置くことである。すなわち

$$\frac{\partial \eta }{\partial t}=-L\frac{\delta F}{\delta \eta }$$

であり、$L>0$ が緩和の速さを決める。

この式が自然である理由は、1次元の最急降下法と対応している点にある。有限次元でエネルギー $E(x)$ を最小化したいとき、$\dot{x}=-\lambda dE/dx$ とすれば $dE/dt=-\lambda (dE/dx{)}^2\le 0$ が成り立つ。Allen–Cahn方程式は、その場 $\eta (\mathbf{r})$ 版であり、勾配が機能微分に置き換わったものと見なせる。

具体的に $F[\eta ]=\int (f(\eta )+\kappa |\mathrm{\nabla }\eta {|}^2/2)d\mathbf{r}$ とすると、前節の結果から

$$\frac{\delta F}{\delta \eta }=\frac{\partial f}{\partial \eta }-\kappa {\mathrm{\nabla }}^2\eta$$

である。したがって

$$\frac{\partial \eta }{\partial t}=-L(\frac{\partial f}{\partial \eta }-\kappa {\mathrm{\nabla }}^2\eta )=-L\frac{\partial f}{\partial \eta }+L\kappa {\mathrm{\nabla }}^2\eta$$

となり、局所ポテンシャルが秩序変数を相の安定点へ引き寄せ、ラプラシアンが空間変化を滑らかにする役割を持つことが読み取れる。界面は局所ポテンシャルが二つの相を選ばせようとする力と、勾配項が急変を嫌う力の釣り合いで現れる。

Allen–Cahn方程式は、相変態のフロント移動、ドメイン壁移動、粒界移動など、界面が動く問題に広く用いられる。ただし、秩序変数が非保存であるため、秩序変数の領域平均は一般に保存されない。固相分率を保存したいような問題では、体積拘束項やラグランジュ未定乗数を加える設計が必要になる場合がある。何を保存すべきかは物理に依存するため、秩序変数の意味づけとともに最初に整理しておく必要がある。

3.3 Cahn–Hilliard 方程式

Cahn–Hilliard方程式は保存場の時間発展を表す基本式である。導出は、保存則と拡散の構成則を組み合わせる手順で行う。まず保存則として連続の式

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

を立てる。次に、拡散の駆動力は化学ポテンシャル $\mu$ の勾配であると仮定し、線形応答として

$$\mathbf{J}=-M\mathrm{\nabla }\mu$$

を置く。$M>0$ はモビリティである。ここで化学ポテンシャルは自由エネルギー汎関数の機能微分で

$$\mu =\frac{\delta F}{\delta c}$$

と定義する。これらを代入して

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\frac{\delta F}{\delta c}\right)$$

を得る。これがCahn–Hilliard方程式である。

具体例として

$$F[c]={\int }_{\mathrm{\Omega }}(f(c)+\frac{{\kappa }_c}{2}|\mathrm{\nabla }c{|}^2)d\mathbf{r}$$

を取ると

$$\mu =\frac{\partial f}{\partial c}-{\kappa }_c{\mathrm{\nabla }}^{2}c$$

である。したがって

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }(\frac{\partial f}{\partial c}-{\kappa }_c{\mathrm{\nabla }}^{2}c)\right)$$

となり、右辺に ${\mathrm{\nabla }}^{4}c$ が現れる。これがCahn–Hilliard方程式が4階の空間微分を含むと言われる理由である。物理的には、界面エネルギー（勾配項）が拡散に影響し、界面近傍の濃度分布を滑らかにしつつ、相分離の駆動力を与える。

Cahn–Hilliard方程式は、スピノーダル分解のように微小ゆらぎが増幅して相分離する現象を自然に表せる。初期の濃度がほぼ一様でも、自由エネルギー曲率が負の領域にあるとき、濃度の空間変動が成長し、二相へ分かれていく。これは、$f(c)$ の形（特に ${\partial }^{2}f/\partial c^2$ の符号）が本質であり、熱力学と数理の結びつきが直感的に見える題材である。

3.4 弱形式と有限要素への接続

偏微分方程式を数値的に解くには、差分法のように強形式をそのまま離散化する方法と、弱形式へ変換してから離散化する方法がある。有限要素法は弱形式を基本とする。弱形式は、支配方程式に試験関数を掛けて領域で積分し、部分積分により微分次数を下げることで得る。

例としてAllen–Cahn方程式

$$\frac{\partial \eta }{\partial t}=-L(\frac{\partial f}{\partial \eta }-\kappa {\mathrm{\nabla }}^2\eta )$$

を考える。試験関数 $v$ を掛けて積分すると

$${\int }_{\mathrm{\Omega }}v\frac{\partial \eta }{\partial t}d\mathbf{r}=-L{\int }_{\mathrm{\Omega }}v\frac{\partial f}{\partial \eta }d\mathbf{r}+L\kappa {\int }_{\mathrm{\Omega }}v{\mathrm{\nabla }}^2\eta d\mathbf{r}$$

である。最後の項を部分積分すると

$${\int }_{\mathrm{\Omega }}v{\mathrm{\nabla }}^2\eta d\mathbf{r}=-{\int }_{\mathrm{\Omega }}\mathrm{\nabla }v\cdot \mathrm{\nabla }\eta d\mathbf{r}+{\int }_{\partial \mathrm{\Omega }}v(\mathbf{n}\cdot \mathrm{\nabla }\eta )dS$$

となる。境界項が消える条件を仮定すれば、弱形式は

$${\int }_{\mathrm{\Omega }}v\frac{\partial \eta }{\partial t}d\mathbf{r}=-L{\int }_{\mathrm{\Omega }}v\frac{\partial f}{\partial \eta }d\mathbf{r}-L\kappa {\int }_{\mathrm{\Omega }}\mathrm{\nabla }v\cdot \mathrm{\nabla }\eta d\mathbf{r}$$

となる。ここで微分が $\eta$ の2階から1階へ下がっており、有限要素の近似空間（1階微分が意味を持つ関数空間）と整合する。

Cahn–Hilliard方程式は4階微分を含むため、弱形式をそのまま作ると2回の部分積分が必要となり、より高い滑らかさを要求する近似空間が必要になる。これを避けるために、化学ポテンシャル $\mu$ を未知量として導入し、

$$\mu =\frac{\partial f}{\partial c}-{\kappa }_c{\mathrm{\nabla }}^{2}c,\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )$$

という2階の連立系へ分割する方法が広く用いられる。この分割により、両方の式が2階となり、標準的な有限要素空間で扱いやすくなる。どの形に変形して解くかは数値手法の選択に属するが、変形の動機は微分次数と近似空間の整合にあると理解すると見通しが良い。

3.5 方程式の比較表

第3章までで現れる2つの基本方程式を、導出の出発点と構造の違いが見える形で整理する。

観点	Allen–Cahn（非保存）	Cahn–Hilliard（保存）
出発点	散逸により局所緩和する仮定	保存則とフラックスの構成
基本式	${\partial }_t\eta =-L\delta F/\delta \eta$	${\partial }_{t}c=\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu ),\mu =\delta F/\delta c$
エネルギー減少	$dF/dt=-\int L(\delta F/\delta \eta {)}^2\le 0$	$dF/dt=-\int M
微分次数	2階（勾配項がある場合）	4階（分割で2階系へ）
物理像	相状態が局所的に切り替わる	物質が流れて濃度が再配分される

この表は暗記のためではなく、どの仮定が式を決めるのかを見える化するための整理である。今後、弾性や外場を加えても、機能微分と保存則・散逸則の組み合わせという骨格は変わらない。したがって、第3章の導出を自分で再現できることが、より複雑なモデルへ進むための最短経路である。

第1〜3章では、フェーズフィールド法が拡散界面の連続場として界面問題を表し、自由エネルギー汎関数と機能微分から支配方程式を導く方法であることを示した。次段階では、界面エネルギーと異方性、多相・多成分の熱力学、弾性や外場などの連成を、同じ導出手順の延長として体系化し、計算結果が物理と対応するための条件をより具体的に詰めることになる。

小まとめ

第1〜3章では、フェーズフィールド法が拡散界面の連続場として界面問題を表し、自由エネルギー汎関数と機能微分から支配方程式を導く方法であることを示した。次段階では、界面エネルギーと異方性、多相・多成分の熱力学、弾性や外場などの連成を、同じ導出手順の延長として体系化し、計算結果が物理と対応するための条件をより具体的に詰めることになる。

参考文献

[1] L.-Q. Chen, Phase-Field Models for Microstructure Evolution, Annual Review of Materials Research 32, 2002.
[2] N. Moelans, B. Blanpain, P. Wollants, An introduction to phase-field modeling of microstructure evolution, Calphad 32, 2008.
[3] PFHub community, Phase Field Method Recommended Practices, NIST.
[4] 日本機械学会 計算力学部門, フェーズフィールド法（MED Wiki）, 2022.
[5] 小山敏幸, フェーズフィールド法（解説記事）, 2025.

第4章 界面、勾配エネルギー、異方性

本章では、フェーズフィールド（拡散界面）モデルが「有限幅の界面」と「界面エネルギー」をどのように生み出すかを、自由エネルギー汎関数の形から順に導く。さらに、界面の曲率が移動の駆動力になる理由、結晶方位に依存する異方性界面エネルギーの入れ方、三重点（多相接点）で接触角が決まる条件、粒界と方位場の表現、応力・電場・磁場など外場が界面応答へ与える効果を、同じ枠組み（自由エネルギー設計）として統一的に整理する。

4.0 拡散界面モデルで何が「物理量」になるか

フェーズフィールドでは、界面は「鋭い幾何学的境界」ではなく、秩序変数（相を表す連続場）$\varphi (\mathbf{r})$ が滑らかに遷移する有限幅の層として表現される。これにより

• 界面を追跡する特別なアルゴリズムなしに、界面の生成・消滅・合体・分岐が自然に起こる
• 曲率駆動・濡れ（接触角）・異方性による形態選択などが、エネルギー最小化と散逸ダイナミクスから一貫して導かれる

という利点が得られる。一方で「界面幅をどれだけ取るか」「その界面幅をメッシュが解像できるか」という数値設計が、物理忠実度と計算コストを左右する。本章ではまず、界面幅と界面エネルギーが自由エネルギーのどの項から決まるかを式で確かめる。

4.1 勾配エネルギー項と平衡界面の導出

最も基本的な二相モデルとして、秩序変数 $\varphi (\mathbf{r})$ を用い、自由エネルギー汎関数を

$$F[\varphi ]={\int }_{\mathrm{\Omega }}(f(\varphi )+\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2)d\mathbf{r}$$

と置く。ここで $f(\varphi )$ は局所自由エネルギー密度（相の安定性を決める項）、$\kappa >0$ は勾配係数である。$|\mathrm{\nabla }\varphi {|}^2$ は $\varphi$ が空間変化するときにエネルギーが増えることを表し、これが有限幅界面を生み出す。

平衡状態は $F$ が極小となる条件、すなわち

$$\frac{\delta F}{\delta \varphi }=0$$

で与えられる。変分計算により

$$\frac{\delta F}{\delta \varphi }=\frac{\partial f}{\partial \varphi }-\kappa {\mathrm{\nabla }}^2\varphi$$

であるから、平衡条件は

$$\kappa {\mathrm{\nabla }}^2\varphi =\frac{\partial f}{\partial \varphi }$$

となる。界面が一次元で法線方向を $x$ とし $\varphi =\varphi (x)$ とすると

$$\kappa \frac{d^2\varphi }{d{x}^2}=\frac{df}{d\varphi }$$

を得る。

局所自由エネルギーに二井戸ポテンシャルを採用する：

$$f(\varphi )=\frac{A}{4}(1-{\varphi }^2{)}^2$$

ここで $\varphi =\pm 1$ が二相の安定点、$A$ が障壁高さ（井戸の深さのスケール）である。このとき

$$\kappa \frac{d^2\varphi }{d{x}^2}=A\varphi ({\varphi }^2-1)$$

を得る。両辺に $d\varphi /dx$ を掛けて積分すると

$$\kappa \frac{d^2\varphi }{d{x}^2}\frac{d\varphi }{dx}=\frac{df}{d\varphi }\frac{d\varphi }{dx}=\frac{df}{dx}$$

より

$$\frac{d}{dx}\left(\frac{\kappa }{2}{\left(\frac{d\varphi }{dx}\right)}^2\right)=\frac{df}{dx}$$

従って

$$\frac{\kappa }{2}{\left(\frac{d\varphi }{dx}\right)}^2=f(\varphi )+C$$

界面から十分離れたバルクでは $\varphi =\pm 1$、$f(\pm 1)=0$、$d\varphi /dx\to 0$ なので $C=0$ である。ゆえに

$$\frac{\kappa }{2}{\left(\frac{d\varphi }{dx}\right)}^2=f(\varphi )(\text{平衡界面の「等分配」関係})$$

が成立する。これを用いると、平衡解は

$$\varphi (x)=\tanh \left(\frac{x}{\delta }\right),\delta =\sqrt{\frac{2\kappa }{A}}$$

となる。$\delta$ が界面厚さ（界面幅）の代表スケールである。

4.2 界面幅と界面エネルギー

4.2.1 界面幅（厚さ）の定義と数値解像条件

界面幅は厳密には定義の揺れがある。例えば

• $-0.9$ から $+0.9$ へ変化する距離
• 勾配が最大となる領域の幅
• $\varphi =0$ を中心にした特定の閾値幅

などが用いられるが、いずれもスケールとしては $\delta =O(\sqrt{\kappa /A})$ に比例する。

数値計算では、この界面幅をメッシュで「見える」ようにする必要がある。典型的な目安として、格子間隔（要素サイズ）を $\mathrm{\Delta }x$ とすると

• 最低でも $\delta /\mathrm{\Delta }x\gtrsim 4$（粗いが破綻しにくい下限）
• できれば $\delta /\mathrm{\Delta }x\gtrsim 6\sim 10$（曲率や異方性を扱うなら推奨）

程度を確保するのが経験的に安全である。界面を薄くしすぎる（$\delta$ を小さくする）と、物理的には鋭い界面極限へ近づくが、数値的には界面が解像できず、界面エネルギーや移動速度が誤って表現される。

4.2.2 界面エネルギー $\gamma$ の導出

界面エネルギー $\gamma$ は、平衡界面を持つ系の自由エネルギーからバルク寄与を引いた「界面に局在する余剰」として

$$\gamma ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}(\frac{\kappa }{2}{\left(\frac{d\varphi }{dx}\right)}^2+f(\varphi ))dx$$

で定義できる。平衡では等分配関係 $\frac{\kappa }{2}({\varphi }^{\prime }{)}^2=f(\varphi )$ があるので

$$\gamma ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\kappa {\left(\frac{d\varphi }{dx}\right)}^{2}dx=2{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(\varphi )dx$$

と書ける。$\varphi =\tanh (x/\delta )$ を用いると、よく知られた結果として

$$\gamma =\frac{2\sqrt{2}}{3}\sqrt{\kappa A}$$

が得られる。したがって

$$\delta \propto \sqrt{\frac{\kappa }{A}},\gamma \propto \sqrt{\kappa A}$$

であり、同じ $(\kappa ,A)$ から界面幅と界面エネルギーが同時に決まる。これは重要で、単純な二井戸＋勾配項モデルでは「界面幅」と「界面エネルギー」を独立に指定できないという制約がある。後の定量モデル（薄界面極限、KKSやグランドポテンシャル、補正項など）の多くは、この制約を緩めて「物理的な $\gamma$ を保ちながら数値的に扱いやすい $\delta$ を取る」ための工夫として理解できる。

4.3 曲率効果と界面移動の力学

非保存秩序変数（相場）$\varphi$ の典型的な時間発展は Allen–Cahn 型である：

$$\frac{\partial \varphi }{\partial t}=-L\frac{\delta F}{\delta \varphi }=-L(\frac{\partial f}{\partial \varphi }-\kappa {\mathrm{\nabla }}^2\varphi )(L>0)$$

これは「自由エネルギーを減少させる最急降下（勾配流）」であり、実際

$$\frac{dF}{dt}=\int \frac{\delta F}{\delta \varphi }\frac{\partial \varphi }{\partial t}d\mathbf{r}=-\int L{\left(\frac{\delta F}{\delta \varphi }\right)}^{2}d\mathbf{r}\le 0$$

となって、時間とともに $F$ が単調減少する（散逸する）構造を持つ。

曲率効果の直感は「曲がった界面は面積が大きく、界面エネルギー $\gamma$ を下げるために縮もうとする」というものである。鋭い界面極限（界面幅が十分小さい極限）で解析すると、界面の法線速度 $V_n$ は平均曲率 ${\kappa }_m$ に比例して

$$V_n=-M_{\gamma }\gamma {\kappa }_m$$

の形になる（$M_{\gamma }$ はモデルパラメータ $L$ と界面内部構造から決まる有効モビリティである）。符号は、凸（正曲率）の界面が縮む方向へ進むことを示す。

この結果は粒界移動、析出物の粗大化（オストワルド成長）などの基礎である。小さな粒子ほど曲率が大きいので消滅しやすく、大きな粒子が成長しやすい、という経験則がフェーズフィールドでは界面追跡なしで自然に再現される。

4.4 異方性界面エネルギーと方位依存

結晶材料では、界面エネルギーは界面の向き（界面法線）に依存し得る。界面法線を

$$\mathbf{n}=\frac{\mathrm{\nabla }\varphi }{|\mathrm{\nabla }\varphi |}$$

とおくと、等方界面では $\gamma (\mathbf{n})={\gamma }_0$ だが、異方性があると

$$\gamma (\mathbf{n})={\gamma }_{0}a(\mathbf{n})$$

のように書ける。2次元なら角度 $\theta$ を用いて

$$a(\theta )=1+{ϵ}_a\cos (m\theta )$$

などが典型である（$m$ は対称性、${ϵ}_a$ は異方性強度）。

4.4.1 勾配項を方位依存にする

等方の勾配項

$$\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2$$

を

$$\frac{\kappa }{2}a(\mathbf{n}{)}^2|\mathrm{\nabla }\varphi {|}^2$$

へ置き換えるのが基本である。ここで注意すべき点は、$a(\mathbf{n})$ が $\mathrm{\nabla }\varphi$ を通じて $\varphi$ に依存するため、機能微分 $\delta F/\delta \varphi$ には単純な ${\mathrm{\nabla }}^2\varphi$ だけでなく「異方性由来の追加項」が現れることである。初学者は次の順序で理解するとよい：

1. まず $a$ を定数に近いとみなして「異方性があるとどの向きが好まれるか」を把握する
2. 次に厳密な変分を追い、追加項が「界面法線の変化に対するエネルギー勾配」を表すことを理解する

4.4.2 強い異方性と界面剛性

異方性が強いと、界面はファセット（平坦面）を持つ形へ近づくが、そのとき「界面剛性（stiffness）」と呼ばれる量が負になり得る。2次元では概念的に

$$\tilde{\gamma }(\theta )=\gamma (\theta )+{\gamma }^{″}(\theta )$$

が正であることが安定性に関わる（厳密な形はモデルにより異なるが、直感として「向きの微小変化がエネルギーを下げてしまうと界面が細かく波打つ」）。$\tilde{\gamma }<0$ に入ると、計算上ギザギザが発生しやすい。

これを避ける代表的な方法は

• 異方性関数 $a(\theta )$ を「剛性が負にならない」範囲に制限する
• 高次勾配項（正則化）を加える（例：$|{\mathrm{\nabla }}^2\varphi {|}^2$ 型の項）
• 数値的には界面幅を十分解像し、時間積分の安定化を行う

である。重要なのは、異方性は形態選択の主要因である一方、強い異方性は連続体モデルの安定性条件とも密接に結びつく、という点である。

4.5 異方性が形態を選ぶ

平衡形状は Wulff 構築で与えられる。概念的には「中心から法線方向 $\mathbf{n}$ へ距離 $\gamma (\mathbf{n})$ の面を引き、その包絡線が平衡形状になる」という作図である。等方なら球（2Dなら円）だが、異方性があると多面体状（2Dなら多角形状）へ近づく。

フェーズフィールドでは界面が有限幅であるため、理想的な鋭いファセットは丸みを帯びる。しかし

• 異方性を強くする
• 界面幅を薄くする（ただし解像条件を満たす）

ことでファセット的形状が現れやすくなる。凝固のデンドライトでは、主枝方向や枝分かれ角度が異方性で決まることが多いため、どの対称性（$m$）を入れるか、${ϵ}_a$ をどの範囲で使うかは、早い段階で設計方針を固める必要がある。

4.6 多相界面と三重点

三重点（3つの界面が交わる点）では、界面エネルギーの釣り合いが接触角として現れる。鋭い界面理論では、界面張力ベクトルの釣り合い

$${\gamma }_{12}{\mathbf{t}}_{12}+{\gamma }_{23}{\mathbf{t}}_{23}+{\gamma }_{31}{\mathbf{t}}_{31}=\mathbf{0}$$

が成立する。等方で ${\gamma }_{12}={\gamma }_{23}={\gamma }_{31}$ なら、接触角は ${120}^{\circ }$ ずつになる。

固体–液体–気体の濡れでよく出る Young の式は、壁面上の力の釣り合いとして

$${\gamma }_{SV}={\gamma }_{SL}+{\gamma }_{LV}\cos \theta$$

の形で書ける（$S$ 固体、$L$ 液体、$V$ 気相）。フェーズフィールドで濡れを扱うときも、本質は「自由エネルギーに壁–相互作用を入れ、その最小化が接触角を決める」という点にある。

4.7 多秩序変数モデル

多相を表す代表的な方法は、多秩序変数 $\{{\eta }_i{\}}_{i=1}^N$ を導入し

$$\sum _{i=1}^N{\eta }_i=1,0\le {\eta }_i\le 1$$

という拘束を課すことである。自由エネルギーの一例として

$$F[\{{\eta }_i\}]=\int [W\sum _i{\eta }_i^2(1-{\eta }_i{)}^2+W_{\mathrm{mix}}\sum _{i<j}{\eta }_i^2{\eta }_j^2+\sum _i\frac{{\kappa }_i}{2}|\mathrm{\nabla }{\eta }_i{|}^2]d\mathbf{r}$$

などが用いられる。ここで

• 第1項：各相が純粋（${\eta }_i=0$ または 1）になるよう促す
• 第2項：相の混合状態を罰し、二相界面や三重点の局在構造を作る
• 第3項：界面に有限幅を与える

という役割を持つ。

三重点の接触角は、モデルパラメータから決まる界面エネルギー ${\gamma }_{ij}$ の組で決まる。したがって「狙った濡れ挙動・接触角」を再現するには、${\gamma }_{ij}$ が意図した値になるように $(W,\kappa )$ や混合罰則の設計を行う必要がある。拘束条件（和が1など）と界面エネルギー設計が不整合だと、物理的でない濡れ（例えば勝手な相が三重点で広がる、角度がメッシュ依存で変わるなど）が生じ得る。

4.8 粒界と結晶方位の表現

多結晶では粒界は「同じ相の中で結晶方位が異なる領域の境界」であり、相境界よりも結晶学的自由度が重要になる。

4.8.1 多秩序変数で粒を表す

各結晶粒を ${\eta }_i$ で表すと、粒界は ${\eta }_i$ と ${\eta }_j$ の遷移層として表現でき、粒界ネットワークや三重点は扱いやすい。一方で粒数が増えると変数が増え、計算量やメモリが増大する。

4.8.2 方位場で連続表現する

結晶方位を連続変数（例えば2Dなら角度場 $\theta (\mathbf{r})$）として表すと、粒界エネルギーを方位差 $\mathrm{\Delta }\theta$ の関数 ${\gamma }_{gb}(\mathrm{\Delta }\theta )$ としてモデル化しやすい。例えば概念的に

$$F_{gb}=\int \frac{K(\varphi )}{2}|\mathrm{\nabla }\theta {|}^{2}d\mathbf{r}$$

のような項は「方位の空間変化（粒界）にエネルギーを課す」ことを意味する。ただし、方位は周期性（$\theta$ と $\theta +2\pi$ の同一性）を持つため、数値的拘束や特異点の扱いが難しくなる。

初学者の段階では、まず多秩序変数モデルで粒界ネットワークの幾何と移動（曲率駆動）を再現し、その後に方位差依存を段階的に組み込む流れが理解しやすい。

4.9 外場による界面応答

応力・電場・磁場などの外場は自由エネルギーに外場項 $f_{\mathrm{ext}}$ として加えるのが基本である：

$$F\to F+\int f_{\mathrm{ext}}(\varphi ,\text{外場},\text{他場})d\mathbf{r}$$

重要なのは「外場そのもの」だけでなく「外場の境界条件」を明確にすることである。代表例を挙げる。

4.9.1 弾性（応力）場との連成

相により固有ひずみ（格子ミスマッチ）があるとき、弾性エネルギー密度を

$$f_{\mathrm{el}}=\frac{1}{2}(\epsilon -{\epsilon }^{\ast }(\varphi )):\mathbf{C}(\varphi ):(\epsilon -{\epsilon }^{\ast }(\varphi ))$$

のように入れる（$\mathbf{C}$ は弾性定数テンソル、$\epsilon$ はひずみ、${\epsilon }^{\ast }$ は固有ひずみ）。このとき境界条件は「変位固定」「外力（応力）固定」「周期境界」などで意味が大きく変わる。

4.9.2 電場との連成

誘電率が相で異なる、あるいは電荷が偏る場合、電場エネルギーを

$$f_{\mathrm{elec}}=\frac{1}{2}ϵ(\varphi )|\mathbf{E}{|}^2(\mathbf{E}=-\mathrm{\nabla }V)$$

と書ける。境界条件としては「電位固定（Dirichlet）」「電束固定（Neumann）」が代表的であり、どちらを課すかで界面の動きが変わることがある。

4.9.3 磁場との連成

磁性が関わる場合、相による磁化や磁気異方性が異なり得る。磁場エネルギーは扱い方が複数あるが、概念的には「磁気エネルギー（Zeeman + 異方性 + 交換 + 静磁）」が $\varphi$ に依存し、界面移動に寄与する。境界条件としては外部磁場の与え方、開境界（無限遠）近似の扱いが鍵となる。

外場項は、界面移動速度・安定形状・ドメイン安定性を変えるため、外場の物理モデルだけでなく境界条件・幾何（電極配置、固定端条件など）を含めた定式化が必須である。

4.10 実装設計の指針

• 界面幅 $\delta$ を決める：物理忠実度と計算負担の折衷である
• メッシュ解像：$\delta /\mathrm{\Delta }x\gtrsim 6$ を一つの目標にする
• 異方性：強くするほど不安定化しやすいので、剛性条件と正則化を意識する
• 多相接触角：拘束条件と界面エネルギー設計の整合を確認する
• 外場：外場項と同じくらい境界条件が重要である

第5章 多相・多成分・化学ポテンシャルの取り扱い

本章では、合金や複数相が関与する系で、組成場（保存量）と秩序変数（非保存または準保存）をどのように整合させるかを、化学ポテンシャルの導出を軸に整理する。特に、界面での局所平衡条件をどう実装するかが、物理性（相図・分配）と数値安定性（振動・スティフネス）の両方に直結する。さらに、反応項（生成消滅）を含む場合のエネルギー整合、多成分拡散で必須となる交差項（Onsager 行列）と正定値性、CALPHAD との接続性を丁寧に扱う。

5.0 保存場と非保存場、化学ポテンシャル

相変態や拡散を扱うとき、変数には大きく2種類ある。

• 保存場：組成 $c(\mathbf{r},t)$ のように全量が保存される（源・吸い込みがなければ）
• 非保存場：秩序変数 $\varphi (\mathbf{r},t)$ のように局所的に増減できる（相の変換を表す）

保存場の駆動力は化学ポテンシャル

$$\mu =\frac{\delta F}{\delta c}$$

であり、フラックスは（最も基本には）

$$\mathbf{J}=-M\mathrm{\nabla }\mu$$

として書ける。これを保存則

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

に代入すると Cahn–Hilliard 型

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\mu \right)$$

が得られる。つまり「何を自由エネルギーに入れるか」が「化学ポテンシャルの形」を決め、拡散の駆動力を決める。

5.1 二相二元系の基本設計：補間型自由エネルギー

二相 $\alpha ,\beta$ を秩序変数 $\varphi$ で表し、組成を $c$ とする。最も単純には

$$f(\varphi ,c)=h(\varphi )f^{\alpha }(c)+(1-h(\varphi ))f^{\beta }(c)+Ag(\varphi )$$

と設計する。ここで

• $h(\varphi )$：補間関数（$h(0)=0,h(1)=1$、滑らか）
• $g(\varphi )$：二井戸（$\varphi =0,1$ を安定化）
• $A$：相境界の障壁（界面構造に効く）

である。勾配項を加え

$$F[\varphi ,c]=\int (f(\varphi ,c)+\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2)d\mathbf{r}$$

とすると、相の熱力学（$f^{\alpha },f^{\beta }$）と界面（$A,\kappa$）が同一汎関数で統合される。

このとき化学ポテンシャルは

$$\mu =\frac{\delta F}{\delta c}=h(\varphi )\frac{d{f}^{\alpha }}{dc}+(1-h(\varphi ))\frac{d{f}^{\beta }}{dc}$$

である。ここでの落とし穴は、界面（$0<\varphi <1$）で $\mu$ が相ごとの化学ポテンシャルの「単純混合」になり、局所平衡（相ごとの $\mu$ が一致する）を自動では満たさない場合がある点である。これが界面での不自然な濃度プロファイルや余分な界面自由エネルギーの原因となり得る。

5.2 化学ポテンシャルと勾配項：Cahn–Hilliard を導く

組成場にも勾配エネルギーを入れることがある（組成界面やスピノーダル分解の代表形）：

$$F[c]=\int (f(c)+\frac{{\kappa }_c}{2}|\mathrm{\nabla }c{|}^2)d\mathbf{r}$$

このとき

$$\mu =\frac{\delta F}{\delta c}=\frac{df}{dc}-{\kappa }_c{\mathrm{\nabla }}^{2}c$$

であり、時間発展は

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\mu \right)=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }(\frac{df}{dc}-{\kappa }_c{\mathrm{\nabla }}^{2}c)\right)$$

となる。右辺に ${\mathrm{\nabla }}^{4}c$ が現れるため、数値的には時間刻み制約が厳しくなりやすい。後の数値解法章で扱う「分割（$\mu$ を導入して2階化）」が、ここで必然的に重要になる。

5.3 単純補間の限界：界面での局所平衡が崩れる

単純補間型は式が短く理解しやすいが、次の問題が起こり得る。

• 相図整合：界面での平衡組成が相図のタイラインとずれる
• 余剰自由エネルギー：界面に不要な化学自由エネルギーが乗り、実効界面エネルギーが変わる
• 数値依存：界面幅や補間関数を変えると分配係数や成長速度が変わりやすい

これらは「界面が有限幅である」ことと「相の自由エネルギーを混ぜる」ことの組合せから来る。よって、単純補間を使う場合は、界面幅・補間関数・メッシュを変えた感度試験を必ず行い、物理的な出力（分配、相分率、成長速度）がどの程度変化するか確認する必要がある。

5.4 KKS 型モデル：局所平衡条件を導く

KKS（Kim–Kim–Suzuki）では、界面での局所平衡を満たすために相ごとの組成を補助変数として導入する。相組成を $c^{\alpha },c^{\beta }$ とし、全組成は

$$c=h(\varphi )c^{\alpha }+(1-h(\varphi ))c^{\beta }$$

で与える。自由エネルギー密度を

$$f_{\mathrm{KKS}}(\varphi ,c^{\alpha },c^{\beta })=h(\varphi )f^{\alpha }(c^{\alpha })+(1-h(\varphi ))f^{\beta }(c^{\beta })+Ag(\varphi )$$

と置く。局所平衡条件は「固定した $(\varphi ,c)$ のもとで $(c^{\alpha },c^{\beta })$ が自由エネルギーを最小化する」ことで導ける。

ラグランジュ未定乗数 $\mu$ を導入し

$$\mathcal{L}=h{f}^{\alpha }(c^{\alpha })+(1-h)f^{\beta }(c^{\beta })+\mu (c-h{c}^{\alpha }-(1-h)c^{\beta })$$

停留条件より（界面では $0<h<1$）

$$\frac{d{f}^{\alpha }}{d{c}^{\alpha }}=\mu ,\frac{d{f}^{\beta }}{d{c}^{\beta }}=\mu$$

すなわち

$${\mu }^{\alpha }={\mu }^{\beta }=\mu$$

が得られる。これは界面での局所平衡（化学ポテンシャル一致）が、拘束付き最小化から自然に現れることを示す。

5.5 KKS の数値実装

KKS では、各格子点・各時刻で与えられた $(\varphi ,c)$ に対し、未知 $(c^{\alpha },c^{\beta },\mu )$ を

$$\{\begin{array}{l}\frac{d{f}^{\alpha }}{d{c}^{\alpha }}=\mu \\ \frac{d{f}^{\beta }}{d{c}^{\beta }}=\mu \\ c=h{c}^{\alpha }+(1-h)c^{\beta }\end{array}$$

で同時に満たすように解く。一般に非線形なのでニュートン法等を用いる。

ここで重要なのは「方程式を増やして複雑にした」のではなく、局所平衡を満たすために必要な自由度を導入し、その代わり拘束式で閉じているという構造である。界面での分配が重要な問題では、この構造が物理性を大きく改善する。

5.6 多成分系への拡張：拡散ポテンシャルと独立性

成分数が $N$ のとき、モル分率 $\mathbf{c}=(c_1,\dots ,c_N)$ は拘束

$$\sum _{i=1}^N{c}_i=1$$

を持つため独立変数は $N-1$ 個である。したがって駆動力も $N-1$ 個の独立な「拡散ポテンシャル」として扱うのが自然である。例えば成分 $N$ を基準にして

$${\tilde{\mu }}_i={\mu }_i-{\mu }_N(i=1,\dots ,N-1)$$

を用いると拘束と整合しやすい。

保存則は

$$\frac{\partial c_i}{\partial t}+\mathrm{\nabla }\cdot {\mathbf{J}}_i=0,\sum _i{\mathbf{J}}_i=\mathbf{0}$$

が必要である。線形非平衡熱力学（Onsager）に従う一般形は

$${\mathbf{J}}_i=-\sum _{j=1}^{N-1}{L}_{ij}\mathrm{\nabla }{\tilde{\mu }}_j$$

であり、$L_{ij}$ は Onsager 行列である。

5.7 多成分KKS

多成分KKSでは相ごとに ${\mathbf{c}}^{\alpha },{\mathbf{c}}^{\beta }$ を導入し

$$\mathbf{c}=h(\varphi ){\mathbf{c}}^{\alpha }+(1-h(\varphi )){\mathbf{c}}^{\beta }$$

局所平衡条件は独立成分 $j=1,\dots ,N-1$ について

$${\tilde{\mu }}_j^{\alpha }={\tilde{\mu }}_j^{\beta }$$

となる。未知数は増えるが、「相図（タイライン）と整合する界面分配」を得るにはこの構造が本質的である。

5.8 グランドポテンシャル形式

多成分化・相図接続を整理しやすい枠組みとして、グランドポテンシャル（大ポテンシャル）形式がある。二元系の要点を示す。

自由エネルギー密度 $f(c)$ に対し

$$\omega (\mu )=f(c)-\mu c$$

を定義する。ただし $\mu =df/dc$ から $c=c(\mu )$ を解き、$\omega$ を $\mu$ の関数とみなす。相ごとに ${\omega }^{\alpha }(\mu ),{\omega }^{\beta }(\mu )$ を持つと、相境界で $\mu$ が連続であることが自然に入る。

汎関数を

$$\mathrm{\Omega }[\varphi ,\mu ]=\int (h(\varphi ){\omega }^{\alpha }(\mu )+(1-h(\varphi )){\omega }^{\beta }(\mu )+Ag(\varphi )+\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2)d\mathbf{r}$$

とし、組成は従属的に

$$c(\varphi ,\mu )=-\frac{\partial }{\partial \mu }(h{\omega }^{\alpha }+(1-h){\omega }^{\beta })=h{c}^{\alpha }(\mu )+(1-h)c^{\beta }(\mu )$$

で求める（$c=-\partial \omega /\partial \mu$ を使用）。時間発展は

$$\frac{\partial \varphi }{\partial t}=-L\frac{\delta \mathrm{\Omega }}{\delta \varphi },\frac{\partial c(\varphi ,\mu )}{\partial t}=\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )$$

で与える。左辺は連鎖律で

$$\frac{\partial c}{\partial t}=\frac{\partial c}{\partial \varphi }\frac{\partial \varphi }{\partial t}+\frac{\partial c}{\partial \mu }\frac{\partial \mu }{\partial t}$$

となり、$\partial c/\partial \mu$ が「容量（感受率）」として $\mu$ の時間変化に効く。初心者はまず「$\mu$ を主変数にすると局所平衡が構造として入り、組成は後から決まる」という流れを掴むとよい。

5.9 CALPHAD との接続

実合金の相図を反映するには CALPHAD の $f^{\alpha }(\mathbf{c}),f^{\beta }(\mathbf{c})$ を使うのが強力である。しかし自由エネルギーの形が複雑になるため

• KKS では各点で解く非線形方程式が重くなる
• グランドポテンシャルでは $\omega (\mathit{\mu })$ の評価と $c(\mathit{\mu })$ の逆変換が重要になる

という実装上の負担が増す。研究目的に応じて「相図忠実度」と「計算効率」のバランスを設計する必要がある。

5.10 反応項と生成消滅

析出反応、相変態、電気化学反応などでは、保存則に反応項（生成消滅）$R$ が加わる：

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

このとき重要な原則は、$R$ が熱力学的駆動力（自由エネルギー差）と整合する形で与えられることである。整合しないと、系が自由エネルギーを増加させる方向へ進み、物理的不安定や数値発散を招き得る。

最も基本的な「エネルギー整合」設計として、反応が局所的な自由エネルギーを減らすよう

$$R\propto -\frac{\delta F}{\delta c}=-\mu$$

のような形を採用する発想がある（実際には化学量論、制約、活性化障壁などを入れるため単純ではないが、設計の核は「$R$ が散逸を生む」ことにある）。

5.11 反応と界面：相変態（Allen–Cahn）と拡散（Cahn–Hilliard）

実際の合金相変態では

• 相の生成消滅（非保存）：$\varphi$ の Allen–Cahn
• 溶質再分配（保存）：$c$ の Cahn–Hilliard（あるいは拡散方程式）

が同時に起こる。典型形は

$$\frac{\partial \varphi }{\partial t}=-L\frac{\delta F}{\delta \varphi },\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )+R(\varphi ,c,\mu )$$

のように書ける。ここで $R$ を入れるなら、どの量が保存され、どの量が反応で変わるか（質量保存・電荷保存など）を先に整理し、整合する反応速度式を置く必要がある。

5.12 多成分拡散と交差項：$-D\mathrm{\nabla }c$ では足りない理由

多成分拡散では、フラックスは単純な

$${\mathbf{J}}_i=-D_i\mathrm{\nabla }{c}_i$$

では一般に不十分であり、基本は「化学ポテンシャル勾配が駆動する」：

$${\mathbf{J}}_i=-\sum _{j=1}^{N-1}{L}_{ij}\mathrm{\nabla }{\tilde{\mu }}_j$$

交差項（$i\ne j$）は「成分 $j$ の勾配が成分 $i$ の流れを誘起する」ことを表し、実合金では普通に起こる。

5.13 Onsager 行列の条件：対称性と正定値性

$L_{ij}$ は次を満たすことが望ましい。

• Onsager 対称性：$L_{ij}=L_{ji}$（適用条件の範囲で）
• 正定値性：任意のベクトル $\mathbf{x}\ne \mathbf{0}$ に対し ${\mathbf{x}}^{T}L\mathbf{x}>0$

正定値性は散逸（自由エネルギーの単調減少）と対応し、数値的にも安定性に直結する。モデル設計段階で「$L$ は正定値である」「どの近似でそうなるか」を明記しておくことが重要である。

5.14 移動度（モビリティ）の設計

移動度 $M$（または $L_{ij}$）は拡散速度・成長速度を決める「運動学の材料定数」である。相ごとに拡散が違うなら

$$M(\varphi )=h(\varphi )M^{\alpha }+(1-h(\varphi ))M^{\beta }$$

のように相依存を入れるのが自然である。さらに

• 界面での数値安定化のために、界面領域だけモビリティを調整する
• 低温での拡散凍結を表すために Arrhenius 型の温度依存を入れる

など、研究目的に応じて設計する。ただし、モビリティをいじると時間スケールが変わり、定量比較が難しくなるため、最終的にどの物理量と同定するか（実験拡散係数、MD、第一原理、文献値）を意識する必要がある。

5.15 自由エネルギー形式の比較

形式	主変数	局所平衡の入り方	追加未知量	特徴
単純補間	$(\varphi ,c)$	自動では満たされない場合がある	なし	簡潔・導入に最適だが界面物理が崩れることがある
KKS	$(\varphi ,c,c^{\alpha },c^{\beta })$	${\mu }^{\alpha }={\mu }^{\beta }$ を拘束で満たす	相組成	分配・相図整合に強いが各点で非線形解が必要
グランドポテンシャル	$(\varphi ,\mu )$	$c=-\partial \omega /\partial \mu$ により構造的に入る	なし（組成は従属）	多成分・CALPHAD接続で整理しやすい

「どれが正しいか」ではなく、「どの物理条件をどこで課しているか」という観点で選ぶのが重要である。

小まとめ

第4章では、勾配エネルギー項から平衡界面の形（$\tanh$ プロファイル）、界面幅 $\delta$、界面エネルギー $\gamma$ を導き、それらが $(\kappa ,A)$ により同時に決まることを示した。さらに、異方性界面エネルギーの導入法、強異方性での安定性、三重点の接触角が界面エネルギーの釣り合いで決まること、粒界の表現として多秩序変数と方位場の考え方、外場（応力・電場・磁場）を自由エネルギーへ入れる際に境界条件が本質であることを整理した。

第5章では、二相二元系から出発して化学ポテンシャルと拡散方程式の導出を確認し、界面で局所平衡を満たすための KKS（拘束付き最小化）とグランドポテンシャル形式（$\mu$ 主変数）を導いた。さらに、反応項を入れる場合の熱力学整合、多成分拡散における交差項と Onsager 行列の正定値性が、物理と数値安定性を同時に守る条件であることを強調した。次に進むべき課題は、これらの定式化を薄界面極限・鋭い界面極限で整理し、材料定数（$\gamma$、モビリティ、拡散係数）との対応を「実装可能な同定手順」として明確化することである。

第6章 弾性・塑性・欠陥場との連成

本章では、フェーズフィールド法を材料の熱力学に基づく場の理論として捉えたうえで、力学場・熱場・流体場・結晶方位場を同一の自由エネルギー枠組みで結び付ける方法を述べる。特に、各場の支配方程式がどのように自由エネルギー汎関数から導かれ、数値計算としてどのような連立問題に落ちるかを、線形代数・微積分・基礎物理を前提に丁寧に記述する。

6.1 弾性エネルギーと固有ひずみ

固相内の相変態や析出では、相ごとに格子定数や結晶対称性が異なるため、同一の空間に相が共存すると変形の不整合が生じる。この不整合を、相（秩序変数）に依存する固有ひずみ（変態ひずみ）${\mathit{\epsilon }}^{\ast }(\mathit{\eta })$ として導入し、弾性エネルギーを自由エネルギーに加えるのが基本である。小ひずみ・線形弾性の範囲では、変位場$\mathbf{u}$からひずみが

$$\mathit{\epsilon }(\mathbf{u})=\frac{1}{2}(\mathrm{\nabla }\mathbf{u}+(\mathrm{\nabla }\mathbf{u}{)}^{\mathsf{T}})$$

で与えられ、弾性ひずみは全ひずみから固有ひずみを差し引いた

$${\mathit{\epsilon }}^{\mathrm{e}}=\mathit{\epsilon }(\mathbf{u})-{\mathit{\epsilon }}^{\ast }(\mathit{\eta })$$

として定義できる。

弾性エネルギー密度は、弾性定数テンソル$\mathbb{C}$を用いて

$$f_{\mathrm{el}}=\frac{1}{2}{\mathit{\epsilon }}^{\mathrm{e}}:\mathbb{C}:{\mathit{\epsilon }}^{\mathrm{e}}=\frac{1}{2}(\mathit{\epsilon }-{\mathit{\epsilon }}^{\ast }):\mathbb{C}:(\mathit{\epsilon }-{\mathit{\epsilon }}^{\ast })$$

である。ここで「$: $」\mathrm{は}\mathrm{二}\mathrm{階}\mathrm{テ}\mathrm{ン}\mathrm{ソ}\mathrm{ル}\mathrm{の}\mathrm{二}\mathrm{重}\mathrm{内}\mathrm{積}\mathrm{を}\mathrm{表}\mathrm{し}、\mathrm{成}\mathrm{分}\mathrm{表}\mathrm{示}\mathrm{で}\mathrm{は}$\mathbf{A}:\mathbf{B}=\sum_{i,j}A_{ij}B_{ij}$\mathrm{で}\mathrm{あ}\mathrm{る}。\mathrm{固}\mathrm{有}\mathrm{ひ}\mathrm{ず}\mathrm{み}\mathrm{は}、\mathrm{析}\mathrm{出}\mathrm{相}\mathrm{の}\mathrm{体}\mathrm{積}\mathrm{ミ}\mathrm{ス}\mathrm{マ}\mathrm{ッ}\mathrm{チ}（\mathrm{等}\mathrm{方}\mathrm{的}）\mathrm{で}\mathrm{あ}\mathrm{れ}\mathrm{ば}$\boldsymbol{\varepsilon}^*=\varepsilon_0\mathbf{I}$の形になり、マルテンサイト変態のようにせん断成分を持てば非対角成分を含む形になる。どの成分が非ゼロであるかが、最終的な形態（板状・針状・ラメラなど）や配列秩序に強い影響を与える。

重要なのは、弾性が長距離相互作用である点である。${\mathit{\epsilon }}^{\ast }$が局所的に変化しても、機械平衡を満たすために変位場$\mathbf{u}$が空間全体で調整され、応力が遠方まで伝播する。フェーズフィールドの局所エネルギー（化学自由エネルギー）や勾配エネルギーが主として界面近傍を支配するのに対し、弾性項は組織スケールの選択則そのものを変え得る。したがって、弾性連成を入れると、同じ化学系でも等方的なドメイン分割から方向性のある配列へと組織が変わることがある。

さらに、弾性不均質（相ごとに$\mathbb{C}$が異なる）を扱う場合、弾性エネルギーは単に局所式を足し合わせるだけでは閉じず、実装上は機械平衡を通じた場の解として定まる。これを理論的に整理する枠組みとして、均質媒質に写像して有効な固有ひずみを定義する「マイクロエラスティシティ」の考え方が広く用いられる。結果として、弾性エネルギーは秩序変数場の汎関数となり、フェーズフィールド方程式における変分$\delta F/\delta \eta$に弾性由来の非局所項が現れる。

6.2 機械平衡と境界条件

準静的弾性（慣性を無視）では、機械平衡は

$$\mathrm{\nabla }\cdot \mathit{\sigma }=\mathbf{0}$$

で与えられる。応力$\mathit{\sigma }$は線形弾性で

$$\mathit{\sigma }=\mathbb{C}:{\mathit{\epsilon }}^{\mathrm{e}}=\mathbb{C}:(\mathit{\epsilon }(\mathbf{u})-{\mathit{\epsilon }}^{\ast }(\mathit{\eta }))$$

である。ここまでで、未知量は変位$\mathbf{u}$と秩序変数$\mathit{\eta }$（あるいは組成$c$）であり、弾性連成は「$\mathit{\eta }$が${\mathit{\epsilon }}^{\ast }$を作り、$\mathbf{u}$が平衡を満たすように決まり、その結果の$\mathit{\sigma }$が再び$\mathit{\eta }$の駆動力に入る」という循環構造を持つ。

境界条件は、(1) 変位境界（ディリクレ条件）$\mathbf{u}=\overline{\mathbf{u}}$、(2) 表面力（トラクション）境界（ノイマン条件）$\mathit{\sigma }\mathbf{n}=\overline{\mathbf{t}}$、(3) 周期境界、のいずれか（または組合せ）で与えることが多い。例えば、外部から引張りを与える問題では、平均ひずみを規定する方法（周期境界＋平均ひずみ制約）と、端面にトラクションを与える方法で、内部応力の分布が変わり得る。固有ひずみがある場合、拘束が強いほど内部応力が大きくなり、析出形状や配列がより顕著に変化しやすい。

数値実装に向けては、弱形式を導くと見通しがよい。仮想変位$\delta \mathbf{u}$に対して

$${\int }_{\mathrm{\Omega }}\mathit{\sigma }:\delta \mathit{\epsilon }dV-{\int }_{\partial {\mathrm{\Omega }}_t}\overline{\mathbf{t}}\cdot \delta \mathbf{u}dS=0$$

が成り立つ（$\partial {\mathrm{\Omega }}_t$はトラクション境界）。ここで$\delta \mathit{\epsilon }=\frac{1}{2}(\mathrm{\nabla }\delta \mathbf{u}+(\mathrm{\nabla }\delta \mathbf{u}{)}^{\mathsf{T}})$である。この形にすると、有限要素法では通常の線形弾性問題として組み立てられ、右辺に固有ひずみ由来の「等価体積力」的な寄与が入ると理解できる。スペクトル法（フーリエ法）を用いる場合は、周期境界のもとでフーリエ空間に変換し、グリーン関数により$\mathbf{u}$や$\mathit{\sigma }$を効率的に求める構成がよく用いられる。

境界条件の選択は、単なる数学的選好ではなく物理の一部である。薄膜・基板拘束では面内平均ひずみがほぼ固定され、バルク材料とは異なる析出配列が生じることがある。微小領域の計算でも、周期境界を使うと「無限に繰り返す材料」を暗に仮定するため、孤立析出物の応力緩和を表しにくい場合がある。したがって、対象の実験条件（拘束の有無、自由表面の存在、基板や治具）を力学境界に翻訳する作業が必要である。

6.3 塑性と不可逆性の導入

弾性だけでなく塑性が関与する問題（加工熱処理、変態誘起塑性、応力緩和を伴う析出など）では、ひずみ分解を拡張するのが基本である。小ひずみ近似のもとで

$$\mathit{\epsilon }={\mathit{\epsilon }}^{\mathrm{e}}+{\mathit{\epsilon }}^{\mathrm{p}}+{\mathit{\epsilon }}^{\ast }$$

とし、${\mathit{\epsilon }}^{\mathrm{p}}$を塑性ひずみとして内部変数にする。すると応力は

$$\mathit{\sigma }=\mathbb{C}:{\mathit{\epsilon }}^{\mathrm{e}}=\mathbb{C}:(\mathit{\epsilon }-{\mathit{\epsilon }}^{\mathrm{p}}-{\mathit{\epsilon }}^{\ast })$$

で与えられる。ここで${\mathit{\epsilon }}^{\mathrm{p}}$が時間発展することが、不可逆性（履歴依存）を表現する核になる。

塑性の発展式は、エネルギー散逸と整合するように構成する必要がある。最も基礎的には、散逸$D$を

$$D=\mathit{\sigma }:{\dot{\mathit{\epsilon }}}^{\mathrm{p}}\ge 0$$

となるように定める（点ごとの局所散逸である）。この不等式が破れると、計算上は外力がなくても塑性が勝手に減少してエネルギーが増えたり、非物理的な振動が生じたりし得る。散逸を保証する構成として、(i) 降伏関数$f(\mathit{\sigma },\dots )\le 0$と整合条件（古典塑性）を入れる方法、(ii) 粘塑性（Perzyna型など）として

$${\dot{\mathit{\epsilon }}}^{\mathrm{p}}=\dot{\gamma }\frac{\partial g}{\partial \mathit{\sigma }},\dot{\gamma }={⟨\frac{f}{{\sigma }_0}⟩}^m$$

のように過応力で滑らかに流す方法、が用いられる。粘塑性は「降伏を厳密に満たす」代わりに、時間刻みの中で徐々に降伏面へ近づく表現であり、数値的に扱いやすいことが多いと捉えるとよい。

フェーズフィールドとの結合は主に2通りである。第一は、塑性を「応力緩和機構」として導入し、固有ひずみが作る応力集中を塑性流動で緩和させる結合である。第二は、塑性による転位密度の蓄積を貯蔵エネルギーとして自由エネルギーに入れ、再結晶や粒成長の駆動力に変換する結合である（第8章と接続する）。いずれにせよ、塑性は拡散や界面移動よりも速いことも遅いこともあり得るため、時間スケールの相対関係が計算の難しさを決める。強連成では、塑性更新とフェーズフィールド更新を交互に行うだけでは収束しない場合があり、同時解法や内点反復が必要になることがある。

6.4 欠陥場、空孔、拡散と応力

空孔や溶質拡散が応力と結合する現象は、熱力学的には「応力が化学ポテンシャルを変える」と表現できる。最も基本的な導出は、自由エネルギー汎関数を

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

とし、組成$c$に関する変分で化学ポテンシャル（より正確には拡散を駆動するポテンシャル）を

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

と定義することである。最後の項$\partial f_{\mathrm{el}}/\partial c$が「応力による寄与」を生み、例えば格子膨張係数（化学ひずみ）を${\mathit{\epsilon }}^{\ast }(c)=\beta (c-c_0)\mathbf{I}$と置けば、$f_{\mathrm{el}}$を介して$\mu$に静水圧応力に比例する項が現れる。直感的には、圧縮応力場では原子の入れ替わりが起きにくくなったり、逆に空孔が集まりやすくなったりする方向に駆動力が変化する、ということである。

拡散方程式は保存則として

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

と書けるため、結局

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\mu \right)$$

というCahn–Hilliard型の形になる。ここで$M$はモビリティであり、温度や相によって変えることも多い（固相では小さく、液相では大きいなど）。この式の中に弾性由来の$\mu$が入ることで、応力場が拡散を偏らせ、偏った拡散がさらに固有ひずみを変えて応力場を変えるという双方向結合が成立する。

欠陥場を連続場として扱う場合、空孔濃度$c_{\mathrm{v}}$や転位密度$\rho$を場変数として導入し、自由エネルギーに欠陥の生成エネルギーや相互作用（例えば$\rho$に比例する貯蔵エネルギー）を加える。たとえば転位密度に対して

$$f_{\mathrm{store}}(\rho )=\frac{1}{2}G{b}^2\rho$$

のような形（$G$は剛性率、$b$はバーガースベクトルの大きさ）を入れると、塑性加工で$\rho$が増える領域が高い内部エネルギーを持つことになり、熱処理時の再結晶や粒成長の駆動力として働く。より進んだモデルでは、転位が作る内部応力場を弾性場に組み込むために、転位密度テンソルや滑り系ごとの内部変数を導入し、相変態と転位構造の相互作用を表すこともある。ただし学部レベルの理解としては、欠陥は「自由エネルギーに追加の蓄積項と散逸の道筋を与えるもの」であり、拡散・相変態の駆動力に新しい空間不均一を持ち込む、と把握しておくのが有効である。

6.5 数値実装上の連成の形

弾性・塑性・欠陥場をフェーズフィールドに結合すると、楕円型方程式（機械平衡）と放物型方程式（拡散、秩序変数の緩和）が同居する連立問題になる。機械平衡は時間微分を含まないため、時間積分の各ステップで「その時刻の秩序変数・組成に整合する弾性場を解く」形になりやすい。この解き方は、計算手順としては自然である一方、強連成では反復が必要になり、反復回数が計算量を支配することがある。

大別すると、(i) 分離型（逐次型）と、(ii) 同時型（モノリシック型）に整理できる。分離型では、時刻$n\to n+1$の更新において、まず弾性場$\mathbf{u}$を解いて$f_{\mathrm{el}}$や$\mathit{\sigma }$を得てから、秩序変数や組成を時間発展させる。この形は実装が比較的単純であり、既存の弾性ソルバを利用しやすい。一方、同時型では、未知ベクトルを$\{\mathbf{u},c,\mathit{\eta },{\mathit{\epsilon }}^{\mathrm{p}},\dots \}$としてまとめ、ニュートン法などで同一反復の中で同時に更新する。強く結合した問題では同時型が収束しやすいことがあるが、ヤコビアン（線形化）の組み立てと前処理が難しくなる。

連成の難しさは、物理の強さだけでなく、数式の剛性（stiffness）にも由来する。界面幅を小さくすると勾配エネルギー項が大きな二階微分を含み、時間刻みが制約されやすい。塑性を硬く（降伏を鋭く）すると、応力更新が非線形になり反復が増えやすい。さらに、拡散モビリティが相で大きく変わると、同じ時間刻みでも場所によって実効的な拡散距離が異なり、数値的に不均一性が増す。これらは「どれか一つの場を高精度化すると他の場の更新が難しくなる」という形で現れるため、モデル化（界面幅、モビリティ、粘性化）と数値解法（陰解法、線形ソルバ、前処理）を同時に設計する必要がある。

以下に、連成の取り扱いを整理するための比較表を示す。

方式	主要な未知	長所	短所
分離型（逐次反復）	まず$\mathbf{u}$、次に$c,\mathit{\eta }$	実装が比較的容易で既存ソルバを使いやすい	強連成で反復が増えやすく、収束しない場合がある
同時型（モノリシック）	$\{\mathbf{u},c,\mathit{\eta },{\mathit{\epsilon }}^{\mathrm{p}},\dots \}$	強連成で収束しやすい場合がある	ヤコビアンと前処理の設計が難しい
スペクトル弾性＋PF	周期境界の$\mathbf{u}$をフーリエで解く	大規模3Dで弾性が高速	周期境界が前提になりやすい
FEM弾性＋PF	任意境界の$\mathbf{u}$を要素で解く	複雑形状・境界条件に強い	反復と線形解法が計算量を支配しやすい

第7章 熱・流体・凝固と成長

7.1 熱伝導と潜熱

凝固や固相変態が起きると、潜熱の放出・吸収によって温度場が変化し、温度が変わることで相変態の駆動力が変化する。したがって、温度$T$は外部条件として与えるだけでなく、エネルギー保存則に基づく場として扱うことが多い。最も基本的な形は熱伝導方程式

$$\rho c_p\frac{\partial T}{\partial t}=\mathrm{\nabla }\cdot (k\mathrm{\nabla }T)+Q$$

であり、$\rho$は密度、$c_p$は比熱、$k$は熱伝導率、$Q$は単位体積あたりの発熱項である。

潜熱は、相の進行を表す秩序変数$\varphi$（液相$\varphi =0$、固相$\varphi =1$など）に結び付けて表すのが分かりやすい。例えば、潜熱$L$を用いて

$$Q=-\rho L\frac{\partial \varphi }{\partial t}$$

と置くと、$\varphi$が増える（凝固）ときに$Q>0$となり発熱（潜熱放出）を表現できる。符号の取り方は$\varphi$の定義に依存するため、凝固で温度が上がる向きに整合するように置くのが要点である。さらに厳密には、温度依存の自由エネルギー$f(\varphi ,c,T)$を用意し、$\varphi$の時間発展式（Allen–Cahn型）を

$$\frac{\partial \varphi }{\partial t}=-L_{\varphi }\frac{\delta F}{\delta \varphi }$$

と書けば、$\delta F/\delta \varphi$に温度依存項が入り、温度場が相変態を駆動する道筋が明示される。

実装上の代表的な工夫として、エンタルピー法（有効比熱法）がある。エンタルピー$H(T,\varphi )$を導入し、

$$\frac{\partial H}{\partial t}=\mathrm{\nabla }\cdot (k\mathrm{\nabla }T)$$

の形で解くと、潜熱を$H$の中に吸収させた扱いができる。これにより、潜熱項を陽に扱うときに起きやすい温度の過剰な振動を抑えやすい。学部段階では、熱方程式は「保存則なので形が崩れにくいが、潜熱が入ると非線形結合が増える」と理解するとよい。

7.2 対流と溶質輸送

溶融池や液相が連続的に存在する凝固では、対流が溶質輸送と界面形態に強く影響する。溶質境界層が対流で薄くなると、同じ冷却条件でも偏析の強さやセル／デンドライトの間隔が変わり得る。流体は非圧縮とみなせる場合が多く、ナビエ・ストークス方程式は

$$\rho (\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \mathrm{\nabla }\mathbf{v})=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot (\mu (\mathrm{\nabla }\mathbf{v}+(\mathrm{\nabla }\mathbf{v}{)}^{\mathsf{T}}))+\mathbf{f},\mathrm{\nabla }\cdot \mathbf{v}=0$$

である（$\mathbf{v}$は流速、$p$は圧力、$\mu$は粘性）。熱対流を入れる場合はBoussinesq近似で体積力$\mathbf{f}$に温度依存の浮力項を入れるが、付加製造の溶融池では表面張力勾配（マランゴニ力）が支配的になることもある。

溶質（組成）輸送は、移流拡散方程式

$$\frac{\partial c}{\partial t}+\mathbf{v}\cdot \mathrm{\nabla }c=\mathrm{\nabla }\cdot (D\mathrm{\nabla }c)$$

で表される（$D$は拡散係数）。ここに相場$\varphi$を組み込むには、$D$を相依存にして固相で小さく液相で大きくする、あるいはフェーズフィールドに整合した「有効拡散方程式」を用いる。さらに、固液界面では分配（$k$値）や界面輸送（溶質の取り込み）が現れるため、薄界面極限の議論と整合させた補正（反捕捉流など）が重要になる（第4章・第5章で扱う薄界面極限と接続する）。

固液共存領域に流体を入れる場合、固相の剛体性を表すために、多孔質抵抗に似た減衰項を運動方程式に追加することがある。例えば

$${\mathbf{f}}_{\mathrm{drag}}=-A(\varphi )\mathbf{v}$$

とし、固相で$A(\varphi )$が大きく液相で小さくなるように設計すると、固相内の流速が抑えられる。こうしたモデル化は、現象の本質（対流が効くのは主に液相）を式に反映するための手段である。

7.3 多相凝固と結晶粒競合

多結晶凝固では、同じ液相から異なる結晶方位が成長し、競合によって柱状晶・等軸晶・混合組織が形成される。フェーズフィールドでこれを扱うには、複数の秩序変数$\{{\varphi }_i\}$で各粒（方位）を表す方法（マルチオーダーパラメータ）と、方位を連続変数として表す方法（方位場モデル）がある（第8章で後者を詳述する）。前者では、各点で

$$\sum _{i=1}^N{\varphi }_i=1,0\le {\varphi }_i\le 1$$

の拘束を課し、${\varphi }_i$が優勢な領域を粒$i$とみなす。

自由エネルギーは、界面エネルギーと障壁項を含む形として

$$F[\{{\varphi }_i\}]={\int }_{\mathrm{\Omega }}(\sum _{i<j}\frac{{\sigma }_{ij}}{2}|\mathrm{\nabla }({\varphi }_i-{\varphi }_j){|}^2+W\sum _{i<j}{\varphi }_i^2{\varphi }_j^2)dV+\text{（熱・溶質項）}$$

のように構成できる（${\sigma }_{ij}$は粒界エネルギー、$W$は障壁の強さ）。この形の利点は、三重点や多重点が界面追跡なしに自然に出現し、粒界ネットワークが自律的に更新される点にある。凝固ではさらに、界面の異方性（結晶方位に依存する界面エネルギーや界面動力学）を入れる必要があるため、勾配項や動力学係数を方位に依存させる。

粒競合の理解には、局所の過冷却と界面法線方向が重要である。温度勾配が大きいと柱状晶が前方へ選択成長しやすいが、核生成（等軸晶）が頻発すると等軸晶が優勢になり得る。フェーズフィールドでは、核生成を確率的に入れる方法もあるが、まずは「既に複数粒が存在する初期条件で競合成長させる」構成で、選択成長の幾何学と異方性の役割を理解するのが有効である。

7.4 付加製造条件での非平衡

付加製造では、局所の冷却速度・凝固速度が非常に大きく、固液界面での局所平衡（分配平衡）が破れる条件が現れる。これを材料学では溶質捕捉や非平衡分配として捉え、分配係数$k$が界面速度$V$に依存して変化する、と考えることが多い。フェーズフィールドで非平衡を表すには、(i) 界面動力学（速度に比例した過冷却）を明示する、(ii) 分配や拡散の時間スケールを界面幅と整合するようにモデル化する、という2つの方向がある。

一例として、界面に有限の運動学的抵抗があると、自由境界問題の境界条件に

$$V={\mu }_k\left(\mathrm{\Delta }G\right)$$

のような形（${\mu }_k$は界面動力学係数、$\mathrm{\Delta }G$は自由エネルギー差）を入れることになる。フェーズフィールドでは、Allen–Cahn型の係数や結合関数を調整することで、この界面動力学を再現する。非平衡が強いと、薄界面極限で成立した補正（反捕捉流など）の前提が崩れる場合もあり、計算格子や時間刻みの要求が厳しくなる。これは「物理として速い」だけでなく、「数値として剛い」問題になりやすいことを意味する。

また、付加製造は熱源が移動するため、温度場が強い時空間勾配を持つ。したがって、相場・組成場・温度場・場合によっては流体場を同時に解くことになり、各場の時間刻みの整合が支配的になる。学部段階では、非平衡凝固を「局所平衡を仮定できないため、界面で起きる輸送を式の中に残す必要がある」と整理するとよい。結果として、モデルの係数同定は難しくなるが、実験条件（走査速度、投入熱量、雰囲気）を直接反映できる利点がある。

7.5 連成計算の目的変数

凝固・成長問題では、出力として何を評価したいかを明確にすることで、必要な物理と数値分解能が定まりやすい。例えば、デンドライトアーム間隔$\lambda$を評価したいなら、界面の曲率と溶質境界層が解像できる格子が必要である。一方、偏析量の統計（最大濃度、分布幅）を評価したいなら、長時間の拡散緩和や対流の寄与が重要になり得る。

代表的な量として、(i) 形態のスケール（${\lambda }_1,{\lambda }_2$など）、(ii) 偏析指標（相ごとの平均濃度、最大／最小、分布）、(iii) 粒径・粒数密度、(iv) テクスチャ（方位分布関数、柱状晶比率）、がある。これらはフェーズフィールドの場から後処理で計算できるが、後処理の定義は一意ではないため、実験での測定定義（例えばEBSDの解析条件、EDSの空間分解能）と揃える必要がある。特に付加製造では、同じ材料でも熱履歴が位置で大きく変わるため、評価量を「どの領域で、どのタイミングで」定義するかが結果を左右する。

さらに、設計問題（条件最適化）へ接続する場合、目的量は複数になりやすい。例えば、微細化（小さい$\lambda$）は強度に有利でも、偏析増大は脆性や耐食性に不利になることがある。フェーズフィールドは、これらのトレードオフを場のレベルで可視化できるが、何を第一に評価するかで必要な連成（熱のみ、熱＋溶質、熱＋溶質＋流体、さらに弾性）も変わる。したがって、出力量と必要物理を対応付ける設計が重要である。

第8章 粒成長・再結晶・結晶方位場

8.1 粒界移動の駆動力

粒成長は、粒界エネルギーを下げる方向に粒界が移動する現象である。曲率$\kappa$を持つ粒界では、界面エネルギー$\gamma$とモビリティ$M$により、界面の法線速度$v_n$が

$$v_n=M\gamma \kappa$$

に比例する、という関係が基礎として知られている（符号は曲率の定義に依存する）。フェーズフィールドでは、粒界を秩序変数の連続遷移層として表現するため、曲率駆動が方程式の中に自然に組み込まれる。例えば単一秩序変数$\varphi$のAllen–Cahn型

$$\frac{\partial \varphi }{\partial t}=-L\frac{\delta F}{\delta \varphi }$$

において、自由エネルギーを

$$F[\varphi ]={\int }_{\mathrm{\Omega }}(W{\varphi }^2(1-\varphi {)}^2+\frac{{ϵ}^2}{2}|\mathrm{\nabla }\varphi {|}^2)dV$$

と取ると、界面では${\mathrm{\nabla }}^2\varphi$が曲率と結び付くため、界面が曲率で動くことになる。

学部レベルの理解としては、勾配項$\frac{{ϵ}^2}{2}|\mathrm{\nabla }\varphi {|}^2$が「界面を作るとエネルギーがかかる」ことを表し、障壁項$W{\varphi }^2(1-\varphi {)}^2$が「相が混ざる中間状態を避ける」ことを表す、と捉えるとよい。界面幅はおおよそ$\delta \sim ϵ/\sqrt{W}$で与えられ、界面エネルギー$\gamma$は$ϵ$と$W$の組合せで決まる。したがって、粒界移動の速度や粒成長の時間スケールは、$L$だけでなく$ϵ,W$にも依存し、物性（実材料の$\gamma ,M$）に合わせるにはこれらの同定が必要となる。

8.2 マルチフェーズフィールド（MPF）

多数の粒を扱うには、粒ごとに秩序変数${\varphi }_i$を割り当て、各点で$\sum _i{\varphi }_i=1$の拘束を課すMPFが基本である。MPFの利点は、粒界のネットワークが自動的に更新され、粒界三重点の運動や粒の消滅が自然に表現できる点にある。さらに、粒界エネルギーを${\gamma }_{ij}$として粒対ごとに与えれば、異方的な粒界エネルギー分布を含む粒成長も表現できる。

自由エネルギーの構成は複数あるが、考え方の核は「(i) 勾配項で界面エネルギー、(ii) 障壁項で相の分離、(iii) 拘束で相の総和を保つ」を同時に満たすことである。ラグランジュ未定乗数$\lambda (\mathbf{x})$を用いれば

$$\tilde{F}=F[\{{\varphi }_i\}]+{\int }_{\mathrm{\Omega }}\lambda (\mathbf{x})(\sum _i{\varphi }_i-1)dV$$

と書け、変分により拘束を満たす方程式が得られる。数値的には、拘束を厳密に満たす更新式を設計する方法と、ペナルティ項で近似的に満たす方法がある。粒数が多いほど未知数が増えるため、計算コストと安定性の両面から、拘束処理の設計が中心課題になる。

MPFで注意すべき点は、粒界が2相界面であるにもかかわらず、多相が不必要に混ざることを抑える必要があることである。三重点近傍では複数の${\varphi }_i$が同時に存在し得るが、それ以外の粒界では基本的に2つの相だけが混ざる状態が望ましい。自由エネルギーと運動方程式の設計により、2相界面の安定性と三重点の力学（ヤングの式に相当する力の釣合い）を再現する。これは、MPFを「粒界ネットワークの連続体模型」として信頼できる形にするための条件である。

8.3 再結晶と貯蔵エネルギー

再結晶は、塑性加工で導入された転位やサブグレインが熱処理により整理され、新しい低欠陥の結晶粒が核生成して成長する現象である。駆動力は主として貯蔵エネルギーの低下であり、転位密度$\rho$が大きい領域ほど内部エネルギーが高い。先に述べたように

$$f_{\mathrm{store}}(\rho )=\frac{1}{2}G{b}^2\rho$$

のような形を自由エネルギーに入れると、$\rho$の空間分布がそのまま再結晶の駆動力分布になる。

フェーズフィールドで再結晶を扱う基本は、(i) 既存粒の方位場（あるいは粒ID）を保持しつつ、(ii) 新粒の核生成と成長を秩序変数で表し、(iii) 駆動力として$f_{\mathrm{store}}$を結合することである。核生成の表現には、確率的核生成（過冷却に基づく核生成率）を入れる方法と、加工組織（EBSDや結晶塑性計算）から「核になりやすい領域」を与える方法がある。学部レベルでは、再結晶は「界面移動（粒界移動）に、追加の体積駆動力（貯蔵エネルギー差）が乗る問題」と捉えると理解しやすい。

式としては、粒$i$を表す秩序変数${\varphi }_i$の運動に

$$\frac{\partial {\varphi }_i}{\partial t}=-L_i\left(\frac{\delta F}{\delta {\varphi }_i}\right)$$

を用意し、$F$の中に貯蔵エネルギーの項を${\varphi }_i$で重み付けして入れる。すると、貯蔵エネルギーの高い領域では新粒が成長してエネルギーを下げ、低い領域では成長が止まりやすくなる。さらに、転位密度$\rho$自体も回復（減少）するため、$\rho$の時間発展式を別に持たせて

$$\frac{\partial \rho }{\partial t}=-K_{\mathrm{rec}}\rho +\text{（変形があれば生成項）}$$

のように結合することもある。ここで$K_{\mathrm{rec}}$は回復速度を表す係数である。こうして、回復→核生成→成長という現象連鎖を同じ枠組みで記述できる。

8.4 方位場の連続表現

多数の粒を${\varphi }_i$で表すMPFは直感的である一方、粒数が増えるほど未知変数が増える。これに対して、方位を連続変数として表す方位場モデルは、未知変数を抑えながら多結晶を表現する道を与える。基本的な考え方は、(i) 相（固液や結晶／非晶）を表す秩序変数$\eta$と、(ii) 結晶方位を表す場$\theta$（2次元なら角度、3次元なら回転を表す変数）を導入し、粒界を「$\theta$が急変する領域」として表すことである。

2次元の基礎形として、自由エネルギーに

$$F[\eta ,\theta ]={\int }_{\mathrm{\Omega }}(f(\eta )+\frac{{ϵ}_{\eta }^2}{2}|\mathrm{\nabla }\eta {|}^2+g(\eta )\frac{{ϵ}_{\theta }^2}{2}|\mathrm{\nabla }\theta {|}^2)dV$$

のような項を入れる。ここで$g(\eta )$は結晶内部（$\eta \approx 1$）では$\mathrm{\nabla }\theta$のコストを小さく、界面や粒界近傍で大きくするように設計されることが多い。直感的には、結晶内部で方位が滑らかに変わること（サブグレイン）と、粒界で方位が急に変わること（高角粒界）を同じ式で表すための重みである。運動方程式は$\eta$にAllen–Cahn型、$\theta$に拡散型（あるいは緩和型）を組み合わせることが多く、粒界移動と方位の緩和が同時に進む。

3次元では、方位は回転群（SO(3)）の元であり、単一の角度では表せない。オイラー角、回転ベクトル、クォータニオンなどの表現があり、それぞれに数値的な特性がある。例えばオイラー角は特異点（ジンバルロック）を持ち得るため、広い方位分布を扱うと不安定になりやすい。クォータニオンは特異点を避けやすいが、単位長制約（$|q|=1$）を保つ必要がある。このように、方位場モデルは「粒界を連続場で表す」代わりに「方位の数学的構造を尊重した変数選び」が重要になる。

8.5 実験との接続

EBSDなどで得られる方位マップは、フェーズフィールド計算の初期条件として非常に有用である。実験マップを格子に写像し、各セルに方位（あるいは粒ID）を割り当てれば、熱処理による粒成長・再結晶の進行を直接比較できる。比較の際に中心となる同定対象は、粒界エネルギー$\gamma$と粒界モビリティ$M$である。粒成長は速度論としては$v_n=M(\gamma \kappa +\mathrm{\Delta }f)$（$\mathrm{\Delta }f$は体積駆動力）に支配されるため、$\gamma$と$M$の組合せが時間スケールを決める。

粒界エネルギーは方位差（ミスオリエンテーション）に依存し、低角粒界では転位配列として解釈できるため、Read–Shockley型の関係が用いられることがある。例えば低角域で

$$\gamma (\theta )={\gamma }_0\theta (A-\ln \theta )$$

のような形が採用される（$\theta$は方位差、$A$は定数）と、方位差が小さい粒界ほどエネルギーが低くなる性質を表現できる。一方で高角域では$\gamma$が飽和することが多く、全域を一つの式で表すには補間や実験同定が必要である。モビリティ$M$も温度依存が強く、Arrhenius型

$$M(T)=M_0\exp (-\frac{Q}{RT})$$

で近似することが多い（$Q$は活性化エネルギー）。したがって、温度履歴が与えられる実験との比較では、$M(T)$の同定が結果に直結する。

再結晶を含む場合、加工履歴をどのように初期条件へ落とすかが鍵になる。結晶塑性計算で得た転位密度や貯蔵エネルギー分布をフェーズフィールドへ渡す方法、あるいはEBSDから推定した局所方位勾配（KAMなど）を貯蔵エネルギーに結び付ける方法がある。これらは「欠陥場を介して力学履歴を熱処理の組織進化へ接続する」考え方であり、第6章の塑性・欠陥場の導入と自然に接続する。最終的には、粒径分布、粒界長さ密度、テクスチャ、再結晶分率など、実験で定義された量と同じ定義で計算から評価することが重要である。

小まとめ

第6章では、固有ひずみに基づく弾性連成を出発点として、塑性と欠陥場を内部変数として導入することで、履歴依存の組織進化を自由エネルギーと散逸の両面から記述できることを示した。第7章では、潜熱を含む熱方程式、対流を含む輸送方程式をフェーズフィールドと結合し、凝固・成長が複数の保存則と緩和則の連立問題として定式化されることを示した。第8章では、粒界移動の速度論をフェーズフィールドで表す考え方、MPFや方位場モデルによる多結晶表現、再結晶を貯蔵エネルギーと結び付ける枠組みを整理した。

今後の発展としては、付加製造のような強非平衡条件に対して、界面動力学・分配非平衡・流体自由表面を整合的に取り込む定式化が重要になる。また、EBSDやその場観察と結び付けた同定問題（$\gamma (\theta )$や$M(T)$、貯蔵エネルギー変換則）を体系化すると、フェーズフィールドは再現計算から予測計算へと役割を広げる。さらに、弾性・塑性・欠陥・熱・流体・方位を同時に扱う計算では、モデル化と数値解法の同時設計が不可欠であり、計算科学としての体系（安定性・収束性・計算量の見積り）を教科書レベルで整理する余地が大きい。

参考ドキュメント

1. Y. U. Wang, Y. M. Jin, A. M. Cuitiño, A. G. Khachaturyan, Phase-field microelasticity theory and modeling of elastically and structurally inhomogeneous solids, Journal of Applied Physics, 92, 1351 (2002).
2. I. Steinbach, Phase-field models in materials science, Modelling and Simulation in Materials Science and Engineering, 17, 073001 (2009).
3. M. Ohno, Quantitative Phase-field Modeling and Simulations of Solidification, ISIJ International, 60, 12 (2020).（J-STAGE）
4. 村松眞由 ほか, サブグレインからの核生成と核成長に関する静的再結晶Phase-Fieldモデル, 材料, 59, 11, 853 (2010).（J-STAGE）
5. 廣内智之 ほか, 高次Multi-phase-fieldモデリングによる粒成長過程の予測, 日本機械学会論文集A編, 77, 782, 1723 (2011).（J-STAGE）
6. H. Henry, H. R. M. Rehberg, N. Goldenfeld, Phase-field model of polycrystalline solidification with a singular coupling between order and orientation, Physical Review B, 86, 054117 (2012).
7. S. Xu, Recent progress in the phase-field dislocation dynamics method, Computational Materials Science, 211, 111 (2022).

第11章 強誘電・磁性・秩序変数場

強誘電体や磁性体では、分極や磁化といった秩序変数が空間的に分布し、ドメインやドメイン壁として観測される。フェーズフィールドは、これらを連続場として記述し、自由エネルギー最小化と散逸を整合させて時間発展を与える枠組みである。

11.1 Landau 型自由エネルギー

強誘電体の相転移を最も基本的に表すには、秩序変数として分極 $\mathbf{P}=(P_x,P_y,P_z)$ を導入し、自由エネルギー密度を $P$ の多項式で近似する。これは、転移点近傍では秩序変数が小さく、自由エネルギーが解析的に展開できるという考えに基づくものである。強誘電の結晶対称性に合わせて許される不変量のみを残すことで、対称性を破る相の安定性が決まる。

最も単純な一次元（単一成分）モデルでは、

$$f_{\mathrm{Landau}}(P;T)=\frac{a(T)}{2}{P}^2+\frac{b}{4}{P}^4+\frac{c}{6}{P}^6$$

のように書く。係数 $a(T)$ は温度で符号が変わることが多く、$a(T)=\alpha (T-T_0)$ と近似すれば、$T<T_0$ で自発分極が生じる状況を表せる。平衡状態は $\partial f/\partial P=0$ を満たす $P$ で与えられ、

$$\frac{\partial f}{\partial P}=a(T)P+b{P}^3+c{P}^5=0$$

となるため、$P=0$ に加えて $P\ne 0$ の解が現れる条件が相転移条件になる。$b>0,c=0$ の単純化では $P^2=-a/b$ が $a<0$ のときに成立し、$T$ を下げると自発分極が連続的に立ち上がることが分かる。$b<0$ で $c>0$ を残せば一次転移も表現でき、ヒステリシスの起源も自由エネルギーの多重極小として理解できる。

実際の強誘電（例：ペロブスカイト）では $P_x,P_y,P_z$ が同等に扱われ、異方性や結晶対称性に応じて交差項が現れる。例えば立方晶の高温相から正方晶相へ移るような状況では、$P_x^2{P}_y^2$ のような項が相対安定性を分ける。したがって Landau 多項式は単なる近似ではなく、どの方向に分極が揃いやすいか、どの相が競合するかを「係数の組」で表す装置である。

空間構造（ドメイン壁）を表すには勾配エネルギーを加える。強誘電では一般に

$$f_{\mathrm{grad}}=\frac{1}{2}\sum _{i,j}{g}_{ij}{\partial }_{i}P{\partial }_{j}P(\text{多成分なら }{g}_{ijkl}{\partial }_i{P}_j{\partial }_k{P}_l)$$

と書き、秩序変数の急峻な変化に罰則を与える。ここからドメイン壁の厚さのスケールが出る。一次元で $g$ を定数とみなし、二重点ポテンシャル（$P=\pm P_s$ が安定）を想定すると、ドメイン壁の幅は概ね

$${\ell }_{\mathrm{DW}}\sim \sqrt{\frac{g}{|a|}}$$

で見積もられ、$g$ が大きいほど壁が広がり、$|a|$ が大きいほど壁が薄くなることが理解できる。

磁性体では秩序変数として磁化 $\mathbf{M}$ を置けるが、多くの強磁性では $|\mathbf{M}|=M_s$ がほぼ一定で方向が変化するため、単位ベクトル $\mathbf{m}=\mathbf{M}/M_s$ を用いる記述が便利である。自由エネルギーは交換・磁気異方性・外場（ゼーマン）・静磁エネルギーの和として与えられ、例えば

$$F[\mathbf{m}]={\int }_{\mathrm{\Omega }}(A|\mathrm{\nabla }\mathbf{m}{|}^2+f_{\mathrm{ani}}(\mathbf{m})-{\mu }_0{M}_s\mathbf{H}\cdot \mathbf{m}+f_{\mathrm{ms}}(\mathbf{m}))dV$$

のように書ける。ここで $A$ は交換剛性であり、強誘電の勾配項と同じく「向きの急変」を抑えて磁壁幅を決める役割を持つ。したがって、強誘電と磁性はいずれも、局所の多項式（相の選択）と勾配項（壁のコスト）の競合としてドメインが生まれるという意味で、数学構造がよく似ている。

外場結合は、強誘電では $-\mathbf{E}\cdot \mathbf{P}$、磁性では $-{\mu }_0\mathbf{H}\cdot \mathbf{M}$ である。外場はドメインの体積分率や壁の移動方向を決め、準静的には自由エネルギーの傾きを与える駆動力とみなせる。特に非線形多項式を含む場合、外場はポテンシャルの非対称化を通じて一方の極小を深くし、ドメインの消滅・生成の障壁を変える。

時間発展式の導入では、強誘電の分極はしばしば緩和型（勾配流）として

$$\frac{\partial P_i}{\partial t}=-L_{ij}\frac{\delta F}{\delta P_j}$$

と置かれる。ここで $L_{ij}$ は運動係数であり、$L_{ij}$ が正定値ならば $dF/dt\le 0$ が保証され、計算はエネルギー散逸と整合する。磁性は保存的回転（歳差運動）を含むため、一般に Landau–Lifshitz–Gilbert (LLG) 方程式

$$\frac{\partial \mathbf{M}}{\partial t}=-\gamma \mathbf{M}\times {\mathbf{H}}_{\mathrm{eff}}+\frac{\alpha }{M_s}\mathbf{M}\times \frac{\partial \mathbf{M}}{\partial t}$$

で与えられる。${\mathbf{H}}_{\mathrm{eff}}$ は有効磁場で、${\mathbf{H}}_{\mathrm{eff}}=-(1/{\mu }_0)\delta F/\delta \mathbf{M}$ として自由エネルギー汎関数から導かれるため、ここでも変分構造が中心にあることが分かる。

11.2 静電場・磁場の長距離相互作用

強誘電ドメインの形態を決める最重要因子の一つは、静電場が生む長距離相互作用である。分極が空間変化すると束縛電荷 ${\rho }_b=-\mathrm{\nabla }\cdot \mathbf{P}$ が生じ、これが電場を作る。電場は再び分極の自由エネルギーに寄与するため、局所的な Landau 項だけでは決まらない「遠隔結合」が現れる。

基本方程式はガウスの法則であり、電束密度 $\mathbf{D}$ に対して

$$\mathrm{\nabla }\cdot \mathbf{D}={\rho }_{\mathrm{free}}$$

が成立する。線形誘電体として $\mathbf{D}={\epsilon }_0{\epsilon }_r\mathbf{E}+\mathbf{P}$ を置き、静電ポテンシャル $\varphi$ を $\mathbf{E}=-\mathrm{\nabla }\varphi$ と定義すれば、

$$\mathrm{\nabla }\cdot ({\epsilon }_0{\epsilon }_r\mathrm{\nabla }\varphi )=\mathrm{\nabla }\cdot \mathbf{P}-{\rho }_{\mathrm{free}}$$

となる。右辺の $\mathrm{\nabla }\cdot \mathbf{P}$ が分極起源のソースであり、電極・自由電荷・界面の誘電率差があれば境界条件として反映される。例えば短絡（電極が完全に遮蔽）なら電極上で $\varphi$ が固定され、開放（電荷供給なし）なら法線方向の $\mathbf{D}\cdot \mathbf{n}$ が制約される。この違いだけで、薄膜における分極の向きが面内に倒れたり、渦状ドメインが現れたりと、形態が大きく変わり得る。

静電エネルギーの書き方は複数あるが、よく用いられる表現の一つは

$$F_{\mathrm{es}}=\frac{1}{2}{\int }_{\mathrm{\Omega }}\mathbf{E}\cdot \mathbf{D}dV$$

である。$\mathbf{D}$ に $\mathbf{P}$ が含まれるため、$F_{\mathrm{es}}$ は分極場とポテンシャル場の連成エネルギーであり、$\varphi$ を消去した非局所汎関数としても見なせる。数値計算では、$P$ を更新するたびにポアソン型方程式を解き直し、電場を自己無撞着に得るという手順になる。

磁性でも、静磁場（反磁界、ストレイフィールド）は長距離であり、ドメイン形成の本質である。準静的で自由電流がないとき、マクスウェル方程式は

$$\mathrm{\nabla }\cdot \mathbf{B}=0,\mathrm{\nabla }\times \mathbf{H}=0$$

である。$\mathbf{B}={\mu }_0(\mathbf{H}+\mathbf{M})$ を用いると

$$\mathrm{\nabla }\cdot \mathbf{H}=-\mathrm{\nabla }\cdot \mathbf{M}$$

が得られ、さらに $\mathrm{\nabla }\times \mathbf{H}=0$ より磁気スカラーポテンシャル $U$ を $\mathbf{H}=-\mathrm{\nabla }U$ と置ける。すると

$${\mathrm{\nabla }}^{2}U=\mathrm{\nabla }\cdot \mathbf{M}$$

というポアソン方程式が現れ、強誘電と同様に「秩序変数の発散」がソースとなって長距離場を生む構造が分かる。反磁界エネルギーは

$$F_{\mathrm{ms}}=\frac{{\mu }_0}{2}{\int }_{{\mathbb{R}}^3}|\mathrm{\nabla }U{|}^{2}dV=\frac{{\mu }_0}{2}{\int }_{\mathrm{\Omega }}\mathbf{M}\cdot {\mathbf{H}}_{d}dV$$

のように書け、境界を含む全空間問題としての性質を持つ。有限領域で計算する場合は、無限遠での減衰条件、開放境界の近似、境界要素法、吸収層、あるいは周期境界に対する Ewald 型和など、境界条件の扱いが計算結果に直結する。

長距離相互作用は計算コストも支配する。周期境界を仮定できる場合、畳み込みを FFT で評価し、$O(N\log N)$ で反磁界や静電場を計算する方法が広く使われる。非周期・複雑形状では有限要素法や境界要素法が自然であるが、連立の規模が大きくなるため、多重格子や前処理付き反復法を用いた解法設計が重要になる。したがって、強誘電・磁性のフェーズフィールドでは、秩序変数の更新則そのものよりも、ポアソン型ソルバの精度と速度が実用的制約になる場面が多い。

11.3 電気機械・磁気機械の連成

強誘電では分極が格子ひずみと強く結合し、圧電応答やドメインスイッチングのしやすさを左右する。結晶学的には、分極の発生に伴って単位胞が歪み、分極方向が変わると歪みも変わるため、力学境界条件（基板拘束、外力、残留応力）が相の安定性を直接変える。

連成の基本は、弾性ひずみ $\mathit{\epsilon }$ と固有ひずみ ${\mathit{\epsilon }}^{\ast }$ を用いて

$$f_{\mathrm{el}}=\frac{1}{2}(\mathit{\epsilon }-{\mathit{\epsilon }}^{\ast }):\mathbb{C}:(\mathit{\epsilon }-{\mathit{\epsilon }}^{\ast })$$

と書くことである。強誘電の場合、固有ひずみは主に電歪（electrostriction）として

$${\epsilon }_{ij}^{\ast }=Q_{ijkl}{P}_k{P}_l$$

のように $P$ の二次で与えられる。ここで $Q_{ijkl}$ は電歪係数であり、分極の向きが変わると ${\epsilon }^{\ast }$ の主軸も変わるため、ドメイン壁の近傍で応力集中が生じやすい。薄膜では面内ひずみが基板で拘束されることが多く、同じ材料でもバルクとは異なるドメインが安定化する。

機械平衡は準静的に

$$\mathrm{\nabla }\cdot \mathit{\sigma }=\mathbf{0},\mathit{\sigma }=\mathbb{C}:(\mathit{\epsilon }-{\mathit{\epsilon }}^{\ast })$$

で与えられ、変位場 $\mathbf{u}$ を用いれば $\mathit{\epsilon }=(\mathrm{\nabla }\mathbf{u}+(\mathrm{\nabla }\mathbf{u}{)}^T)/2$ である。したがって、強誘電フェーズフィールドは $(\mathbf{P},\varphi ,\mathbf{u})$ の三者連成になり、ポアソン方程式と弾性平衡方程式を同時に満たす必要がある。数値的には、時間ステップごとに $\mathbf{P}$ を仮更新し、$\varphi$ と $\mathbf{u}$ を解き、得られた $\mathbf{E}$ と $\mathit{\sigma }$ を使って $\mathbf{P}$ を補正するという自己無撞着手順で整合を取る。

磁性でも磁歪がひずみと磁化方向を結び付ける。固有ひずみを磁化の関数として与える代表的な書き方として、等方近似では

$${\epsilon }_{ij}^{\ast }=\frac{3}{2}{\lambda }_s(m_i{m}_j-\frac{1}{3}{\delta }_{ij})$$

が用いられることがある。より一般には立方晶の磁気弾性エネルギーを

$$f_{\mathrm{me}}=b_1({\epsilon }_{xx}{m}_x^2+{\epsilon }_{yy}{m}_y^2+{\epsilon }_{zz}{m}_z^2)+2b_2({\epsilon }_{xy}{m}_x{m}_y+{\epsilon }_{yz}{m}_y{m}_z+{\epsilon }_{zx}{m}_z{m}_x)$$

のように表し、$b_1,b_2$ が磁歪起源の結合強度となる。ここから分かるのは、ひずみは単に磁壁を押す外力ではなく、異方性項そのものを作り替え、ドメインの容易軸・困難軸を変えるという点である。逆に、磁区が変われば固有ひずみが変わるので、弾性エネルギーも変化し、磁化とひずみは双方向に影響する。

電気機械と磁気機械の共通点は、連成項がドメイン壁のエネルギー地形を変えることである。強誘電では電場と応力を同時に与えると、単独では起きないスイッチング経路が現れることがある。磁性でも応力（あるいは磁歪応力）が磁壁のピン止め条件を変え、磁化反転の臨界場や動力学を変える。したがって連成を入れることは、現象を増やすためではなく、実験条件に対応した自由エネルギーを正しく構成するために不可欠である。

11.4 ナノ構造と境界条件の重要性

強誘電・磁性のドメインは、バルクよりも薄膜・多層膜・ナノパターンで顕著に変わる。理由は、長距離場が境界条件に敏感であり、表面や界面での電荷遮蔽・反磁界・機械拘束が、体積の Landau 項に匹敵する寄与を持つからである。フェーズフィールドでは境界条件が方程式の一部であり、境界条件の選択自体がモデルの物理仮定を表している。

強誘電薄膜で最も重要な概念の一つは脱分極場である。面外分極があると表面に束縛電荷が現れ、これが遮蔽されない場合、分極に逆向きの電場が生じて面外分極を不安定化する。短絡電極なら遮蔽が強く、面外一様分極が安定になり得るが、開放表面や不完全遮蔽では、分極が面内に回る、あるいは渦状・閉曲線状の配置で束縛電荷を減らす方向にドメインが組み替わる。このとき、静電方程式の境界条件（電位固定か、$\mathbf{D}\cdot \mathbf{n}$ 制約か、界面誘電率や死層の有無）が、直接ドメイン形態へ写る。

機械境界条件も同等に重要である。基板拘束により面内ひずみが固定されると、電歪で生じる固有ひずみと衝突し、ある分極方向が不利になる。その結果、バルクで安定な相が薄膜では不安定になったり、分極回転が起きたりする。さらに、ナノパターンや探針印加では応力分布が非一様となり、局所応力が核生成サイトを作る。したがって、電気境界条件と機械境界条件は独立ではなく、同時に与えたときの整合が観測量を決める。

磁性でも形状が決定的である。反磁界は試料外形に強く依存し、細長い形状では長軸方向が容易軸になる形状異方性が現れる。薄膜では面外成分が反磁界で不利になりやすく、面内磁化が優勢になる。多層膜では層間の反磁界結合や界面異方性が加わり、スキルミオンや迷路状磁区などが安定化する条件が変わる。これらはすべて、ポアソン方程式で求める静磁場が境界を通じて変わることに由来する。

境界における秩序変数の条件も整理しておく必要がある。勾配項を含む自由エネルギーから自然境界条件を導くと、例えば強誘電では

$$\mathbf{n}\cdot \frac{\partial f_{\mathrm{grad}}}{\partial (\mathrm{\nabla }{P}_i)}=0$$

のような形になり、単純化すれば $\mathbf{n}\cdot \mathrm{\nabla }{P}_i=0$（法線方向勾配ゼロ）に相当する。これは「表面で分極が自由に調整できる」仮定であり、表面吸着や化学状態で分極が固定される状況とは異なる。磁性でも交換エネルギーに由来する境界条件として $\mathbf{n}\cdot \mathrm{\nabla }\mathbf{m}=0$ が現れるが、実際には表面異方性や粗さで境界条件が変わる。したがって、ナノ構造の計算では境界のモデル化が支配的になり得る。

また、境界条件を時間依存で与える設計は、実験再現にとどまらず、材料設計問題に直結する。例えば探針電圧の波形や走査速度を入力として与え、得られるドメイン配置を目的に沿って誘導するという考え方は、境界条件が制御入力であることを意味する。フェーズフィールドはそのまま「制御入力に対する空間場の応答系」と見なせるため、最適化や逆問題へ展開しやすい。

11.5 機械学習による代替モデル

強誘電や磁性の三次元計算は、長距離場（静電・静磁）と弾性場の連成により、計算負荷が急激に増える。特にナノ構造で高解像度を要求すると、ポアソン型方程式を多数回解く必要があり、時間発展を長く追うことが難しくなる。そこで、数値計算の出力データを学習し、場の更新を近似する代替モデルが導入されることがある。

最も直感的な代替は、状態 $\mathbf{P}(\mathbf{x},t)$ や $\mathbf{m}(\mathbf{x},t)$ から次時刻の状態を直接予測する写像

$${\mathcal{G}}_{\theta }:\text{state}(t)\mapsto \text{state}(t+\mathrm{\Delta }t)$$

を学習する方法である。畳み込みニューラルネットワーク（U-Net など）を用いれば、局所パターン（ドメイン壁、核、欠陥近傍）を捉えた更新が可能になる。一方で、長距離相互作用を含む現象は「局所演算の積み重ね」だけでは表現が難しいため、フーリエ領域で演算する operator learning（Fourier Neural Operator など）や、ポアソン方程式の解作用素そのものを近似するモデルが検討されている。これにより、電場・反磁界の計算を明示的に解かずに、効果を学習器に内包させる方針が立つ。

しかし、代替モデルが物理に反しないためには制約が必要である。強誘電の緩和型発展では $dF/dt\le 0$ が基本であるため、学習器がエネルギーを増大させる更新を繰り返すと破綻する。そこで、エネルギー汎関数を同時に学習して勾配流を模倣する、あるいは損失関数に散逸条件を組み込むなど、変分構造を守る工夫が導入される。磁性では $|\mathbf{m}|=1$ 制約が本質であり、予測後に正規化するだけでは歳差運動の位相が壊れ得るため、球面上の力学を保つパラメータ化（例えば指数写像）や、LLG の構造を組み込んだネットワーク設計が望ましい。

代替モデルの価値は、単に高速化にとどまらない。計算が速くなれば、境界条件や材料係数の探索回数が増やせ、データ同化や逆問題に踏み込みやすくなる。例えば探針印加の条件から望みのドメインを得る設計では、反復探索が多くなるため、代替モデルが研究戦略を変える可能性がある。一方で、学習データの分布から外れた条件（電極形状の変更、材料定数の大幅変更、欠陥導入）では誤差が急増し得るため、適用範囲を数理的に点検し、必要ならば高忠実度計算へ戻る運用が必要である。

強誘電に関しては、三次元の探針誘起スイッチングを対象に、機械学習でフェーズフィールド計算を加速する研究が報告されている。加えて、強誘電フェーズフィールドを対象にした学習器の設計・比較が進み、精度・頑健性・汎化の観点で議論が積み上がっている。磁性でも、反磁界計算の高速化や、スピンダイナミクスの時系列予測を狙った代替が検討されており、長距離相互作用と制約条件をどう組み込むかが共通課題である。

小まとめ

強誘電・磁性のフェーズフィールドは、局所の Landau（あるいは異方性）項、勾配項、長距離場（静電・静磁）という三層の競合がドメインを作るという点で統一的に理解できる。さらに弾性連成が入ると、境界条件を通じて相の安定性と壁運動の経路が変わり、薄膜・ナノ構造ではその効果が支配的になり得る。

今後は、長距離場ソルバと連成の数値安定性を押さえた上で、境界条件を制御入力として扱う設計問題へ展開する方向が有望である。機械学習の代替モデルは探索回数を増やすための手段になり得るが、変分構造や制約条件を守る形で統合し、適用範囲の点検を体系化することが、強誘電・磁性のフェーズフィールドを設計へ接続する鍵になる。

第12章 数値解法と離散化設計

フェーズフィールド法は偏微分方程式の連立として実装されるため、空間と時間の離散化の選び方が、得られる形態・速度・保存則の再現性を直接左右する。ここでは、代表的な離散化と解法を、式の形と数値的な意味づけが追える形で整理する。

12.1 空間離散（差分・有限要素・スペクトル）

フェーズフィールドで扱う場は、秩序変数 $\varphi (\mathbf{x},t)$ や組成 $c(\mathbf{x},t)$ のように連続場であり、勾配 $\mathrm{\nabla }$ やラプラシアン ${\mathrm{\nabla }}^2$ が頻繁に現れる。まず連続式を離散化して計算機上の有限自由度に落とす必要があるが、そのときに何を保持し、何を近似誤差として受け入れるかが設計である。特に界面が有限幅で表現されるため、界面の勾配を正しく表す格子設計が最重要となる。

差分法（有限差分・有限体積を含む）は、規則格子上で導関数を差分で置き換える方法である。例えば 1 次元で格子幅 $h$ のとき、中心差分は

$${\frac{\partial \varphi }{\partial x}|}_i\approx \frac{{\varphi }_{i+1}-{\varphi }_{i-1}}{2h},{\frac{{\partial }^2\varphi }{\partial x^2}|}_i\approx \frac{{\varphi }_{i+1}-2{\varphi }_i+{\varphi }_{i-1}}{h^2}$$

で与えられる。2 次元・3 次元ではラプラシアンの近似に 5 点・7 点・27 点など複数のステンシルがあり、等方性誤差（格子方向に依存した見かけの異方性）が形態に混入し得る。界面張力や異方性を議論する問題では、この離散化由来の方向性が、物理的異方性と区別しづらくなるため注意が要る。

有限要素法（FEM）は、複雑形状や多様な境界条件を自然に扱える方法である。FEM では、支配方程式を弱形式（変分形式）に直してから離散化する。例えば Allen–Cahn 型の勾配流

$$\frac{\partial \varphi }{\partial t}=-M\mu ,\mu =\frac{\delta F}{\delta \varphi }$$

を考え、$\mu$ にラプラシアンが含まれる典型形として

$$\mu =\frac{\partial f(\varphi )}{\partial \varphi }-\kappa {\mathrm{\nabla }}^2\varphi$$

を置く。試験関数 $v$ を掛けて領域 $\mathrm{\Omega }$ で積分すると

$${\int }_{\mathrm{\Omega }}v\frac{\partial \varphi }{\partial t}d\mathrm{\Omega }=-M{\int }_{\mathrm{\Omega }}v(\frac{\partial f}{\partial \varphi }-\kappa {\mathrm{\nabla }}^2\varphi )d\mathrm{\Omega }$$

となり、部分積分で

$${\int }_{\mathrm{\Omega }}v{\mathrm{\nabla }}^2\varphi d\mathrm{\Omega }=-{\int }_{\mathrm{\Omega }}\mathrm{\nabla }v\cdot \mathrm{\nabla }\varphi d\mathrm{\Omega }+{\int }_{\partial \mathrm{\Omega }}v(\mathrm{\nabla }\varphi \cdot \mathbf{n})dS$$

を用いると境界条件が自然に現れる。Neumann（フラックス）条件は境界積分項で制御でき、Dirichlet（値固定）条件は関数空間の拘束として与える。フェーズフィールドでは $\mathrm{\nabla }\varphi \cdot \mathbf{n}=0$ のような勾配ゼロ条件（界面の法線方向勾配が境界で消える条件）を採用することが多く、弱形式との整合が取りやすい。

スペクトル法（主にフーリエスペクトル）は、周期境界条件のもとで強力である。$\varphi (\mathbf{x})$ をフーリエ展開 $\hat{\varphi }(\mathbf{k})$ に変換すると、導関数が代数演算に置き換わる。例えば

$$\widehat{\mathrm{\nabla }\varphi }=i\mathbf{k}\hat{\varphi },\widehat{{\mathrm{\nabla }}^2\varphi }=-|\mathbf{k}{|}^2\hat{\varphi }$$

である。Cahn–Hilliard のように高階微分が現れても、$|\mathbf{k}{|}^4$ の乗算に変わるため実装が簡潔で、線形部分を厳密に扱いやすい。一方で複雑形状・非周期境界では使いにくく、また FFT に伴う全体通信が並列性能を支配する場合がある。

以上より、規則格子かつ周期境界であれば差分法やスペクトル法が高効率であり、複雑形状・多様な境界・局所細分を要する場合は FEM が優位である。どの方式でも、界面幅 $\ell$ を格子幅 $h$ で十分解像する条件（例えば界面に 6〜10 点程度を置く設計）が必要であり、これは後述の時間刻み制約とも結びつく。

12.2 時間積分（陽・陰・半陰）

時間方向の離散化は、安定性と計算コストの交換である。フェーズフィールドでは拡散項（${\mathrm{\nabla }}^2$）や高階拡散（${\mathrm{\nabla }}^4$）が支配的になりやすく、方程式が剛性を持つ。剛性とは、空間分解能を上げるほど、安定に積分するために必要な時間刻みが急激に小さくなる性質である。

まず Allen–Cahn（非保存）と Cahn–Hilliard（保存）の基本形を確認する。自由エネルギー汎関数を

$$F[\varphi ]={\int }_{\mathrm{\Omega }}(f(\varphi )+\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2)d\mathrm{\Omega }$$

とすると、変分微分は

$$\mu =\frac{\delta F}{\delta \varphi }=\frac{\partial f}{\partial \varphi }-\kappa {\mathrm{\nabla }}^2\varphi$$

である。Allen–Cahn は

$$\frac{\partial \varphi }{\partial t}=-M\mu$$

であり、Cahn–Hilliard は保存則を組み込んで

$$\frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot \left(M\mathrm{\nabla }\mu \right),\mu =\frac{\delta F}{\delta c}$$

となる。Cahn–Hilliard は $\mu$ に ${\mathrm{\nabla }}^{2}c$ が含まれ、さらに $\mathrm{\nabla }\cdot (M\mathrm{\nabla }\mu )$ によって ${\mathrm{\nabla }}^{4}c$ が現れるため、時間積分がより厳しくなる。

陽解法（明示法）は、例えば前進オイラーで

$${\varphi }^{n+1}={\varphi }^n+\mathrm{\Delta }t\mathcal{L}({\varphi }^n)$$

の形で進める。実装は簡潔であるが、拡散項は一般に $\mathrm{\Delta }t\le C{h}^2$（Cahn–Hilliard ではさらに厳しい $\mathrm{\Delta }t\le C{h}^4$ に近い制約）を要求し、3 次元や界面幅を小さくした設定では現実的でない刻みになる。これは、短波長（高波数）モードが強く減衰・増幅するために起こる現象であり、線形安定性解析で確認できる。

陰解法（暗示法）は

$${\varphi }^{n+1}={\varphi }^n+\mathrm{\Delta }t\mathcal{L}({\varphi }^{n+1})$$

のように未知量で右辺を評価するため、大きな $\mathrm{\Delta }t$ を許容し得る。ただし毎ステップで非線形連立を解く必要があり、線形ソルバ・前処理が性能を決める。Crank–Nicolson のような 2 次精度法は精度面で有利であるが、非線形性との組合せで反復が難しくなる場合がある。

半陰解法（semi-implicit）は、剛性の主因となる線形部分を陰に、非線形部分を陽に分ける設計である。例えば Allen–Cahn を

$$\frac{\partial \varphi }{\partial t}=-M(f^{\prime }(\varphi )-\kappa {\mathrm{\nabla }}^2\varphi )$$

とし、$f^{\prime }(\varphi )$ を陽、${\mathrm{\nabla }}^2\varphi$ を陰として

$$\frac{{\varphi }^{n+1}-{\varphi }^n}{\mathrm{\Delta }t}=-M(f^{\prime }({\varphi }^n)-\kappa {\mathrm{\nabla }}^2{\varphi }^{n+1})$$

とすれば、未知量は線形楕円型方程式として解ける。Cahn–Hilliard でも同様に、線形の高階項を陰に取り込み、非線形の化学ポテンシャル部分を陽に置くことで、安定性と計算量の両立を図る。半陰解法は、設計の自由度が高い反面、エネルギーが単調減少する性質（勾配流としての物理性）を離散化で破らない工夫が必要になる。

時間積分の選択は、界面が移動する時間スケール、拡散の時間スケール、外場や連成場の時間スケール（弾性は準静的、流体は慣性を含むなど）を分けて考えると整理しやすい。特に凝固や電気化学では、熱・電位・輸送の時定数が異なるため、最も厳しい制約がどの式から来るかを特定し、その式に合わせて半陰化や分割を設計するのが筋道である。

12.3 非線形解法とヤコビアン

非線形項は、二重井戸ポテンシャルの $f^{\prime }(\varphi )$、濃度依存モビリティ $M(c)$、弾性エネルギーの $\varphi$ 依存、電気化学反応の指数則など、多くの場所に現れる。非線形連立は一般に、離散化後の未知ベクトル $\mathbf{u}$ に対して

$$\mathbf{R}(\mathbf{u})=\mathbf{0}$$

という残差方程式として書ける。ニュートン法は、反復 $k$ において

$$\mathbf{J}({\mathbf{u}}^k)\mathrm{\Delta }{\mathbf{u}}^k=-\mathbf{R}({\mathbf{u}}^k),{\mathbf{u}}^{k+1}={\mathbf{u}}^k+\alpha \mathrm{\Delta }{\mathbf{u}}^k$$

を解く方法であり、$\mathbf{J}=\partial \mathbf{R}/\partial \mathbf{u}$ がヤコビアンである。$\alpha \in (0,1]$ は減衰係数であり、強非線形や大きい時間刻みでは $\alpha <1$ が収束に必要な場合がある。

ヤコビアンの構築は、収束回数と総計算時間を左右する。厳密ヤコビアンは反復回数を減らすが、構築コストが大きい。近似ヤコビアンは 1 回の反復が軽いが、反復回数が増えやすい。実際には、(i) 自動微分や記号微分を用いて厳密性を確保する、(ii) ブロック構造を保った近似（物理ごとの近似）を作る、(iii) 行列を作らず作用だけを与える行列フリーで Krylov 法を回す、などの方針が取られる。

多物理連成では、未知量が $\mathbf{u}=(\varphi ,c,{\mathbf{u}}_{\mathrm{disp}},T,\mathrm{\Phi },\dots )$ のように増え、ヤコビアンはブロック行列になる。例えば $(c,\mu )$ の Cahn–Hilliard 分割では

$$\left(\begin{array}{cc}I/\mathrm{\Delta }t& -\mathrm{\nabla }\cdot (M\mathrm{\nabla }\cdot )\\ -\frac{{\partial }^{2}f}{\partial c^2}-\kappa {\mathrm{\nabla }}^2& I\end{array}\right)\left(\begin{array}{c}\mathrm{\Delta }c\\ \mathrm{\Delta }\mu \end{array}\right)=-\left(\begin{array}{c}{R}_c\\ R_{\mu }\end{array}\right)$$

のような形が現れ、前処理はこのブロック構造に合わせるのが定石である。ブロック対角近似、Schur 補完近似、AMG（代数マルチグリッド）をブロックごとに当てる設計などが、計算規模が大きいほど有効になる。

線形ソルバの収束は、物理パラメータにも依存する。例えば界面幅を小さくすると勾配項の係数が効き、剛性が増して条件数が悪化しやすい。ヤコビアンのスケーリング（未知量の無次元化や対角スケーリング）、境界条件の扱い、時間刻みの増加戦略（小さい $\mathrm{\Delta }t$ から徐々に増やす）を組み合わせることで、ニュートン反復の挙動が安定する。

12.4 Cahn–Hilliard の分割と安定化

Cahn–Hilliard は 4 階微分を含むため、そのまま 1 変数で離散化すると、C^1 連続要素や特殊差分が必要になる。そこで化学ポテンシャル $\mu$ を導入して 2 階系に分割するのが広く用いられる。自由エネルギーを

$$F[c]={\int }_{\mathrm{\Omega }}(f(c)+\frac{\kappa }{2}|\mathrm{\nabla }c{|}^2)d\mathrm{\Omega }$$

とすると

$$\mu =\frac{\delta F}{\delta c}=f^{\prime }(c)-\kappa {\mathrm{\nabla }}^{2}c$$

であり、支配式は

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

である。これにより未知量は $(c,\mu )$ の 2 変数となるが、空間微分は 2 階までで済む。FEM では $c,\mu \in H^1(\mathrm{\Omega })$ の標準的な関数空間に乗るため、実装とソルバ設計が容易になる。

ただし分割すると、数値安定性の中心は、非線形な $f^{\prime }(c)$ をどう時間離散化するかに移る。特にスピノーダル分解では、初期の微小ゆらぎが指数成長する波数帯を持つため、離散化がエネルギーの時間変化を不自然にすることがある。勾配流としての性質を離散化でも保ちたい場合、エネルギー安定な時間離散化を採用するのが有効である。

凸分割（convex splitting）は、自由エネルギー密度 $f(c)$ を凸部 $f_{\mathrm{cvx}}$ と凹部 $f_{\mathrm{ccv}}$ に分け、

$$f(c)=f_{\mathrm{cvx}}(c)+f_{\mathrm{ccv}}(c)$$

として、凸部を陰に、凹部を陽に扱う設計である。例えば

$${\mu }^{n+1}=f_{\mathrm{cvx}}^{\prime }(c^{n+1})+f_{\mathrm{ccv}}^{\prime }(c^n)-\kappa {\mathrm{\nabla }}^2{c}^{n+1}$$

のように置くと、適切な条件のもとで離散エネルギーが単調減少する性質を示せる。これにより、時間刻みを大きくしてもエネルギー発散を起こしにくくなり、長時間計算での信頼性が増す。

安定化の別系統として、線形安定化項を加えた半陰法や、SAV/IEQ と呼ばれる補助変数を導入してエネルギー安定性を確保する方法もある。重要なのは、どの方式でも、保存則（総量 $\int cd\mathrm{\Omega }$ の保存）とエネルギー散逸（$dF/dt\le 0$）が、離散化でどの程度保たれるかを確認できる形に式を整えることである。確認は、時間ステップごとの $\int cd\mathrm{\Omega }$ と $F$ の履歴を出力し、理論と整合する単調性・保存性が得られているかを見ることで実行できる。

12.5 適応メッシュと誤差制御

界面幅が有限であるという特徴は、界面近傍だけ細かい分解能を置けば十分であることも意味する。例えば界面幅を $\ell$、格子幅を $h$ とすると、界面に必要な点数 $N_{\ell }$ を

$$N_{\ell }\approx \frac{\ell }{h}$$

として設計でき、$N_{\ell }$ を一定以上に保つことが、界面張力・曲率駆動・界面異方性の再現性に直結する。3 次元で一様格子にすると計算量が爆発しやすいので、界面近傍に自由度を集中する適応メッシュは大きな効果を持つ。

適応の鍵は、どこを細かくするかを定める誤差指標である。フェーズフィールドでは、最も単純には $|\mathrm{\nabla }\varphi |$ や $|\mathrm{\nabla }c|$ を用い、勾配が大きい領域（界面）を細分する設計が取れる。より理論的には、離散残差 $\mathbf{R}(\mathbf{u})$ の局所大きさや、回復勾配（数値解の勾配を滑らかな空間に射影して差を取る）を用いた事後誤差評価を用いる。破壊の損傷場や、弾性・静電のような長距離場が連成する場合は、界面近傍だけでなく応力集中や電位勾配が支配的になる領域にも自由度が必要であり、単一指標では不足し得る。

適応メッシュは、時間積分とも結び付く。局所的に $h$ を小さくすると陽解法の安定条件が厳しくなるため、半陰化や陰解法の採用が実質的に必要になる。FEM では適応に伴って線形系の性質が変わるため、AMG などの前処理を安定に機能させる工夫（スムーザ選択、レベル構成、係数スケーリング）が性能を左右する。

誤差制御は、数値誤差とモデル誤差を分けて扱う視点を含む。例えば界面幅 $\ell$ を物理量として小さくするほどモデルは鋭い界面に近づくが、同時に離散化誤差が増えやすい。したがって、(i) $\ell$ を固定して格子収束を見る試験、(ii) 格子比 $h/\ell$ を固定して $\ell$ を変える試験、の双方を行うと、数値誤差とモデル側の近似（拡散界面近似）の影響を分離して議論できる。

12.6 手法比較表

観点	差分法	有限要素法	スペクトル法
形状適応	弱い（規則格子向き）	強い（複雑形状・局所細分）	周期境界で強い
実装量	少ない	中〜多	中
高階項	工夫が必要（分割や広いステンシル）	分割により標準要素で対応	フーリエで自然に扱える
長距離場連成	工夫が必要	既存基盤と親和	フーリエで高速化可能だが境界制約あり
大規模並列	容易	容易（前処理が鍵）	FFT 通信が鍵
適応自由度	実装に工夫	標準的に扱える	原理的に相性が悪い

この表は、方法そのものの優劣ではなく、境界条件・形状・連成場・必要な解像度が選択を決めることを示す整理である。例えば凝固のデンドライトで周期境界・大規模統計を狙う場合は差分法やスペクトル法が効率的であり、電池の電極粒子形状や複雑境界を含む場合は FEM の利点が表に出やすい。最終的には、数値誤差が目的変数の不確かさに与える影響を見積もり、必要十分な設計を選ぶのが筋道である。

第13章 検証、妥当性、再現性

フェーズフィールド計算は、モデルの仮定と数値解法の仮定が重なって結果を作るため、何が正しさの根拠であるかを明示しなければならない。本章では、解析的整合、コード同士の比較、パラメータ同定、実験比較、記録方法を、式と手順の形で整理する。

13.1 解析解・収束試験

最初に行うべき検証は、解析的に性質が分かる問題で、離散化の誤りをあぶり出すことである。例えば 1 次元の静的界面は、二重井戸 $f(\varphi )=\frac{W}{4}({\varphi }^2-1{)}^2$ を用いた Allen–Cahn の平衡条件 $\mu =0$、すなわち

$$f^{\prime }(\varphi )-\kappa \frac{d^2\varphi }{d{x}^2}=0$$

から導かれる。解は

$$\varphi (x)=\tanh \left(\frac{x}{\sqrt{2}\xi }\right),\xi =\sqrt{\frac{\kappa }{W}}$$

の形になり、界面幅が $\xi$ でスケールすることが分かる。数値実装では、この解を初期条件にして時間発展させ、形が保たれること、$\mu$ が小さいこと、エネルギーが時間で変化しないことを確認できる。これにより勾配項、境界条件、時間積分の符号ミスなどが早期に検出できる。

収束試験は、格子幅 $h$ と時間刻み $\mathrm{\Delta }t$ を変えたときに、数値解が一定の次数で収束するかを確認する試験である。例えば解 ${\varphi }_h$ と参照解 ${\varphi }_{\mathrm{ref}}$ の差の $L^2$ ノルム

$$\parallel {\varphi }_h-{\varphi }_{\mathrm{ref}}{\parallel }_{L^2}={({\int }_{\mathrm{\Omega }}({\varphi }_h-{\varphi }_{\mathrm{ref}}{)}^{2}d\mathrm{\Omega })}^{1/2}$$

を計算し、$h$ を半分にしたときに誤差が $p$ 次で減る（$\propto h^p$）かを見る。時間方向も同様に $\mathrm{\Delta }t$ を変えて収束次数を見る。フェーズフィールドでは界面幅が有限であるため、$h$ を小さくしても、モデル側の近似（鋭い界面に対する拡散界面近似）が残る。したがって、$h/\ell$ を固定して $h$ を変える試験と、$\ell$ を変える試験を分けることで、離散化誤差とモデル近似の影響を分離できる。

解析解が無い場合には、製造解（manufactured solution）を用いる。任意の滑らかな関数 ${\varphi }^{\star }(\mathbf{x},t)$ を用意し、それを支配方程式に代入して外力項（ソース）$s(\mathbf{x},t)$ を逆算し、

$$\frac{\partial \varphi }{\partial t}=\mathcal{L}(\varphi )+s$$

という形にして $s$ を実装する。数値解が ${\varphi }^{\star }$ に収束すれば、空間・時間離散化、境界条件実装、ヤコビアン実装の整合が確認できる。連成系でも同様に、各場に製造解を与えて整合を取ることで、連成項の符号や係数の取り違えを検出しやすい。

13.2 ベンチマーク問題の活用

解析解と収束試験は局所的な正しさを保証するが、複雑な非線形発展に対しては、コミュニティで合意されたベンチマークが有効である。ベンチマークは、初期条件、境界条件、無次元化、評価指標（エネルギー履歴、相分率、構造因子など）を固定し、コード間で比較可能にしたものである。フェーズフィールドでは、拡散や界面項の実装差、時間積分の違い、ヤコビアン近似の違いが、どの出力にどの程度現れるかを定量的に観察できる。

ベンチマークを使うときは、(i) 仕様の無次元量をそのまま実装する、(ii) 自分のコードの内部表現が有次元である場合は、変換式と単位系を明示する、(iii) 仕様が要求する出力（例えば $F(t)$、平均組成、構造因子 $S(k)$）を同じ定義で計算する、という順序で整える。ここで定義の違いが入ると、見かけの差が物理差なのか後処理差なのか判別できない。特にスペクトル法と FEM の比較では、構造因子の離散フーリエ変換の定義や、メッシュ上の補間が差を生みやすい。

ベンチマーク結果の比較は、単に一致・不一致を判定するだけでなく、差の原因を分解する材料になる。例えば、エネルギーが単調減少しないなら時間積分設計に問題がある可能性が高い。保存量がずれるなら、Cahn–Hilliard の離散化や境界フラックス実装を疑うのが筋道である。ベンチマークは、研究対象が違っても数値基盤の品質保証に使えるため、研究テーマごとに個別に検証を作る労力を減らす効果がある。

13.3 パラメータ同定

フェーズフィールドのパラメータは、自由エネルギー密度 $f$、勾配係数 $\kappa$、モビリティ $M$（あるいは緩和係数）などである。これらは、実験や第一原理計算、分子動力学から与えることが多いが、与え方の根拠が曖昧だと、結果の信頼性が議論できない。したがって、パラメータがどの観測量を再現するために選ばれたかを式で結び付ける必要がある。

例えば二相の界面エネルギー $\gamma$ と界面幅のスケールは、1 次元界面解から関係づけられる。平衡界面で

$$\gamma ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left(\kappa {\left(\frac{d\varphi }{dx}\right)}^2\right)dx$$

の形が得られ、$\varphi =\tanh (x/(\sqrt{2}\xi ))$ を代入すると $\gamma$ が $\kappa$ と $W$（二重井戸の高さ）で決まることが分かる。これにより、実験や第一原理から得た $\gamma$ をターゲットにして $(\kappa ,W)$ を設定する筋道ができる。多成分・多相では界面ごとに ${\gamma }_{ij}$ が異なり得るため、モデルでその自由度を持たせるか、あるいは平均化した値を使うかを選ぶ必要がある。

動力学パラメータ（モビリティ）$M$ は、界面速度や拡散係数と結び付く。Cahn–Hilliard の線形化では、自由エネルギーの曲率 $f^{″}(c_0)$ により成長率が決まり、初期のスピノーダル分解の時間スケールから $M$ の見積りが可能である。凝固では、鋭い界面理論の毛管長や界面運動学的係数に合わせて、フェーズフィールドパラメータを調整する薄界面解析が用いられる。ここでは、界面厚さを実際より厚くしても、適切な補正を入れることで鋭い界面極限を再現する設計が行われる。

同定には不確かさが付随する。界面エネルギーや拡散係数の誤差が、枝間隔や粒径分布にどう影響するかは非線形であるため、感度解析（パラメータを振って出力の変化を見る）を明示するのが望ましい。複数のパラメータが同じ出力に影響する場合は、同定が非一意になり得るため、複数の観測量（例えば界面エネルギーと共存組成、拡散係数と界面速度）を同時に合わせる設計が必要になる。

13.4 実験データとの比較戦略

フェーズフィールドの出力は場そのものであり、実験で観測される量（粒径、偏析、応力、電圧など）とは表現が異なる。したがって、計算変数を観測量に写像する式を決めてから比較するのが筋道である。例えば組織の比較では、相分率、2 点相関関数、界面曲率分布、粒径分布など、統計量に落として比較する方が頑健である。EBSD で得られる方位マップと粒界ネットワークを比較するなら、粒界角度分布や粒径分布の一致を指標にできる。

物性比較では、追加のモデルが要る場合がある。例えば電気抵抗は、相分布だけでは決まらず、相ごとの抵抗率と連結性が必要である。そこで相場 $\varphi$ から局所導電率 $\sigma (\varphi )$ を定義し、準静的な電位場

$$\mathrm{\nabla }\cdot (\sigma (\varphi )\mathrm{\nabla }V)=0$$

を解いて有効抵抗を算出する、といった連成が必要になる。電池の電圧は、電気化学ポテンシャル差や過電圧の定義と結び付けて計算し、実験の電圧曲線と同じ定義で比較する必要がある。

時間と空間の観測窓の違いも重要である。実験画像は有限の視野と時間分解能を持ち、フェーズフィールドは連続場を任意にサンプリングできる。比較の際は、実験の分解能に合わせて計算結果を畳み込む（空間フィルタを掛ける、閾値で相を二値化する、時間平均を取るなど）ことで、観測条件の差を減らせる。逆に、計算側だけが高分解能で議論すると、実験で観測不能な微細構造の差に議論が寄ってしまう危険があるため、観測条件を揃えた上で議論するのが整合的である。

13.5 記録と公開の作法

再現性は、同じ入力を与えれば同じ出力が得られるという意味に留まらず、第三者が同じ条件を再構成できるという意味を含む。そのためには、入力パラメータだけでなく、単位系、無次元化、メッシュ生成条件、境界条件の実装、ソルバ設定、収束判定、初期条件生成の方法（乱数を含む場合は乱数種）を機械可読に残す必要がある。図の生成条件（閾値、平滑化、カラーバー範囲）も比較に影響するため、後処理の設定も記録対象に含めるのが望ましい。

コードの版管理は、実装が変わるほど重要になる。Git のコミット ID、依存ライブラリの版（MPI、線形代数、ソルバ）、コンパイラと最適化オプションを保存することで、計算結果の差が環境由来かモデル由来かを判別できる。環境を固定する方法としては、コンテナ（Docker/Singularity/Apptainer）や、パッケージマネージャのロックファイルを用いるのが有効である。これにより共同研究者が同じ環境を構築しやすくなる。

公開に際しては、入力デッキと最小限の出力（再計算に必要な初期条件や中間状態）を、永続的識別子の付くリポジトリに置くのが望ましい。データ量が大きい場合は、全時系列ではなく、検証・比較に必要な統計量や代表スナップショットを選ぶ設計も取り得る。重要なのは、第三者が計算の前提（何を固定し、何を自由にしたか）を読み取れる形にまとめることである。

第14章 計算基盤と高速化

フェーズフィールドは局所演算の比率が高く並列化に適する一方、長距離場連成や FFT、非線形ソルバが性能の律速になることが多い。本章では、計算量とメモリを式で見積もり、どこに投資すべきかを判断できるように整理する。

14.1 大規模並列とメモリ設計

3 次元格子で $N_x\times N_y\times N_z$ 点、変数数を $N_{\mathrm{var}}$ とすると、単純な配列格納でもメモリは $O(N_x{N}_y{N}_z{N}_{\mathrm{var}})$ で増える。例えば倍精度 8 byte を用い、${1024}^3$ 点で $N_{\mathrm{var}}=5$ なら、単一配列だけで約 ${1024}^3\times 5\times 8\approx 40$ GB 程度になる。さらに時間積分で過去値や中間ベクトル、ソルバの作業領域が必要であり、実際にはこの数倍が要る。したがって、計算を成立させるためには、変数数の管理（秩序変数を増やし過ぎない、拘束を別変数にせず消去する等）と、ソルバが要求する補助ベクトル数を見積もる必要がある。

領域分割（domain decomposition）は、多くの実装で基本となる。差分法では、各 MPI ランクが部分領域を持ち、境界にゴーストセルを置いて近傍通信を行う。フェーズフィールドの局所演算は近傍参照で済むため、通信は主に表面に比例し、計算は体積に比例する。したがって、十分大きい局所領域を持てば並列効率が上がるが、局所領域が小さくなるほど通信が支配的になる。高階ステンシルや ${\mathrm{\nabla }}^4$ を直接差分化する場合は参照範囲が広がり、ゴースト層が厚くなるため、通信量が増える点も見積もりに入れる。

スペクトル法では FFT が中心となり、全体通信（all-to-all）が不可避になる。高並列では、スラブ分割よりペンシル分割が必要になり、通信パターンが複雑化する。FFT の性能が律速であれば、局所計算を最適化しても全体は速くならないため、問題の設計段階で、スペクトル法の利点（高精度、線形項処理の簡潔さ）と通信コストを比較する必要がある。

FEM では自由度が不規則に分布し、行列・ベクトルの分散格納と前処理が重要になる。AMG は不均一係数や複雑メッシュに強いが、粗視化レベルの構築コストとメモリが大きい場合がある。したがって、メッシュを細分し過ぎる前に、界面幅と必要解像度を固定して、どの程度の自由度で目的変数が収束するかを見積もることが、計算基盤設計の中心になる。

14.2 オープンソース計算枠組み

研究で繰り返し現れる作業は、方程式の弱形式（あるいは差分式）の実装、境界条件の付与、非線形解法、線形ソルバと前処理、入出力、並列化である。これらを自作すると、研究対象の変更に追随できる利点がある一方、ソルバや並列 I/O を含む基盤部分の品質保証に大きな労力が要る。そこで、方程式記述とソルバ基盤を分離した枠組みを用い、研究者はモデルとパラメータに集中する設計が広く採用されている。

FEM 系では、フェーズフィールド方程式を部品として提供する枠組みが存在し、Allen–Cahn、Cahn–Hilliard、弾性連成、反応拡散などを組み合わせやすい。重要なのは、入力が抽象化されているほど、実装ミスの余地が減る反面、内部で採用されている離散化や時間積分の詳細を理解しないと、モデル差と数値差の切り分けが難しくなる点である。したがって、枠組みを使う場合でも、13 章の検証（解析解、収束、ベンチマーク）を自分の設定で必ず通し、基盤の仕様が自分の問題に対して十分であることを確認する必要がある。

差分法系では、単純な規則格子コードが最も軽量であるが、並列化や I/O を含む基盤を整えた枠組みも存在する。規則格子はデンドライトやスピノーダル分解など、周期境界で統計性を取りたい問題と相性が良い。一方で複雑境界が重要な電池電極や多孔質媒体では、FEM の方が境界条件を忠実に実装しやすい。枠組みの選択は、モデルの複雑さよりも、境界条件と幾何の要求によって決まる側面が大きい。

枠組みを使う最大の利点は、ソルバと前処理が整備されていることである。非線形連成の反復は、基盤が提供する Newton–Krylov や AMG によって現実的なコストに抑えられることが多い。逆に、単純な陽解法で済む問題では、枠組みの抽象化が過剰になり、最小実装の方が速い場合もある。したがって、目的変数の要求精度、必要な連成、計算規模を見積もり、基盤を重くする価値がどこにあるかを判断するのが整合的である。

14.3 行列フリーと高性能実装

大規模問題では、疎行列を陽に組み立てて保持することが、メモリと帯域の面で不利になることがある。行列フリー（matrix-free）は、行列 $\mathbf{A}$ を作らず、行列ベクトル積 $\mathbf{y}=\mathbf{A}\mathbf{x}$ を、弱形式や局所演算に基づいて直接計算する設計である。これにより、メモリ使用量を減らし、演算をキャッシュに乗りやすくして高速化できる場合がある。高次要素（高次数 FEM）では、和因子化により 1 自由度あたりの演算量が下がりやすく、行列フリーの利点が出やすい。

ただし Krylov 法の反復回数は前処理に依存し、前処理が弱いと反復回数が増えて総時間が伸びる。したがって、行列フリーは、(i) 作用は速いが反復が多い、という状況になり得る。これを避けるために、幾何マルチグリッドや p-マルチグリッド、近似 Schur 補完など、作用の形に合わせた前処理設計が重要になる。フェーズフィールドでは、拡散型の楕円演算子が中心であるためマルチグリッドと相性が良い一方、強連成（弾性、電位、反応）ではブロック構造を保つ前処理が必要になる。

行列フリーは、ヤコビアンの与え方とも結び付く。ニュートン法では $\mathbf{J}$ の作用が必要であり、厳密ヤコビアンを行列として持たずに作用だけ与える設計が可能である。さらに、$\mathbf{J}$ を厳密に作らず、有限差分近似で作用を評価する Jacobian-free Newton–Krylov もある。これらは実装の自由度を増やすが、収束性と前処理の設計がより重要になる。特に Cahn–Hilliard のような保存系では、保存量を壊さない作用の設計が必要である。

GPU では、局所演算が多い差分法は相性が良く、FEM でも行列フリーと高次要素により GPU 効率を上げる方向が取られる。GPU 実装では、メモリアクセスと通信が性能を支配するため、データ配置、カーネル融合、通信隠蔽が中心課題になる。フェーズフィールドは時間ステップ数が多いことが多いため、単ステップの数 % の改善が総計算時間に大きく効く。

14.4 国内計算基盤と応用例

国内では、凝固・粒成長の大規模フェーズフィールド計算が、GPU スパコンや共用計算資源を用いて報告されている。これらの研究では、デンドライトの競合成長や粒成長統計など、十分大きい系でないと評価できない量を対象にし、計算規模を上げる動機が明確である。例えば粒成長では、粒数が少ないと統計ゆらぎが支配し、理論的スケーリング則との比較が難しくなるため、超大規模計算により統計性を確保する方向が取られる。

凝固では、多結晶の競合成長やミリメートルスケールの指向性凝固を扱うために、領域サイズと時間ステップ数が大きくなる。GPU を用いた並列化は、局所演算が主体であるフェーズフィールドの特性を活かしやすく、長時間計算を可能にしている。さらに、付加製造のように温度勾配・走査速度・溶融池流れを含む場合は、熱・流体・溶質の連成が必要になり、単純な単独フェーズフィールドよりも計算基盤の整備が重要になる。

共用計算資源では、課題設定の段階で、必要な格子数、変数数、時間ステップ数を見積もり、必要なメモリと計算量を提示することが求められる。これは研究者にとっても有益であり、なぜその規模が必要か、どの物理を入れるとどれだけ重くなるかを定量的に理解できる。結果として、モデルの拡張（弾性や電位の導入）が、計算基盤の制約と整合する形で設計されやすくなる。

国内の研究報告を参照する際には、方程式そのものより、並列化の戦略と数値検証の提示に注目すると学べる点が多い。例えば、格子収束やエネルギー散逸の確認が示されていれば、計算規模が大きくても数値の信頼性を議論できる。逆に、規模だけが強調されて検証が薄い場合は、結果の物理的解釈に慎重であるべきである。

14.5 機械学習による時間発展の代替

時間発展をそのまま解く代替として、機械学習で時間発展写像を近似する研究が進んでいる。代表的には、(i) 組織画像を低次元表現に圧縮する（オートエンコーダ等）、(ii) 圧縮空間で時間発展を学習する（RNN、Transformer、Neural ODE 等）、(iii) 必要なら高次元へ復元する、という設計が取られる。これにより、同一条件での長時間予測や、多数条件の走査が可能になる場合がある。

一方で、フェーズフィールドは自由エネルギーに基づく勾配流であり、保存則やエネルギー散逸という構造を持つ。代替モデルがこの構造を壊すと、見かけの形態は似ていても、相分率保存やエネルギーの時間挙動が不自然になる危険がある。したがって、学習の損失関数に保存則の誤差やエネルギー散逸の違反を入れる、あるいは支配方程式残差を損失に入れる物理情報付き学習が有効である。Neural operator や PINO のように、方程式構造を取り込む方向もこの文脈で理解できる。

代替モデルの有効性は、誤差の定義に依存する。画素ごとの差が小さくても、界面位置がわずかにずれると、曲率駆動の速度や粒径統計に大きく影響することがある。したがって、比較指標として、相分率、界面長（面積）、粒径分布、相関関数、構造因子などを用い、物理量として意味のある距離で評価する必要がある。さらに、学習データの範囲外（温度勾配や外場が違う条件）に外挿すると破綻しやすいため、適用範囲を明示し、必要なら不確かさ推定も含めて扱うのが整合的である。

14.6 推奨事項に基づく設計整合

数値設計を体系化するためには、共通の言語で条件を揃える必要がある。その中核は、無次元化とスケールの固定である。界面幅 $\ell$、格子幅 $h$、時間刻み $\mathrm{\Delta }t$、モビリティ $M$、勾配係数 $\kappa$ などを無次元群として整理し、どの群が支配しているかを示すと、実装や論文間の比較が容易になる。例えば $h/\ell$ を固定した比較は、離散化差を減らし、モデル差を見やすくする。

また、ベンチマーク仕様が公開されている場合、それを参照して数値設計を整えると、比較可能性が大きく上がる。エネルギー履歴、保存量、収束基準、出力の定義が揃うことで、議論が数値差から物理差へ移りやすくなる。さらに、入力・出力・後処理の定義を明文化することは、共同研究でのデータ交換を円滑にし、同一条件の再計算を可能にする。

設計整合の実行は、具体的には以下のような形に落ちる。第一に、解析解や製造解を用いた単体試験を用意し、実装の根本を確認する。第二に、保存量とエネルギー散逸を常に出力し、時間積分の性質を監視する。第三に、格子収束と時間収束を行い、目的変数が許容範囲で安定する領域を見つける。これらを最初に整えることで、モデル拡張（弾性、電位、反応、流体、機械学習代替）を行っても、数値の根拠が崩れにくくなる。

小まとめ

本章群では、離散化（空間・時間）と非線形解法の設計が、フェーズフィールド計算の物理的整合（保存則とエネルギー散逸）に直結することを、式の形で確認した。今後は、ベンチマークや公開基盤を用いた比較可能性の確保と、行列フリー・マルチグリッド・GPU を含む高性能実装の発展が、より複雑な連成問題（電気化学、破壊、強誘電・磁性）を現実的な計算量で扱う鍵になると期待される。

おわりに

フェーズフィールド法は、自由エネルギー汎関数を起点として支配方程式を導出するという一貫した数理構造を持ち、界面を含む複雑な材料現象を連続体の枠組みで統一的に記述できる点に本質的な強みがある。本書で扱ってきたように、材料組織形成、破壊、電気化学、強誘電・磁性といった一見異なる現象群も、秩序変数とエネルギー汎関数の設計という共通の視点からモデル化できることは、フェーズフィールド法の汎用性と拡張性を端的に示している。

一方で、モデルの自由度が高いことは、物理的妥当性や数値結果の信頼性を常に問い直す必要があることも意味する。本書では、熱力学的一貫性、変分原理に基づく定式化、数値離散化の選択、検証と再現性といった要素を段階的に整理し、単に「計算が回る」モデルではなく、「検証可能な計算」を構築するための判断軸を示してきた。今後、ベンチマーク問題や推奨事項の体系化が進むことで、手法の成熟度はさらに高まっていくと考えられる。

将来展望としては、(i) 弾性・塑性・電磁場・化学場を含む多物理場連成の高忠実化、(ii) 実験・第一原理計算と接続したパラメータ同定および不確かさの定量化、(iii) 行列フリー法や適応メッシュを含む高性能数値計算技術の導入、(iv) 機械学習を用いた近似モデルや逆問題への展開といった方向が、相互に影響しながら発展していくと予想される。特に、設計問題への応用においては、物理モデルに基づくフェーズフィールド法とデータ駆動的手法をどのように接続するかが重要な課題となる。

本書を通じて読者が身につけるべき最も重要な点は、特定の現象や実装技術そのものではなく、「自由エネルギーから出発して、どのようにモデルを立ち上げ、どの段階で仮定や近似を導入しているのか」を自ら説明できる視点である。その視点を持つことで、新しい対象や未経験の問題に対しても、フェーズフィールド法を単なる計算手法ではなく、思考の枠組みとして柔軟に適用できるようになるだろう。本書が、そのための確かな基盤となることを願っている。

参考ドキュメント

• PFHub: The Phase Field Community Hub（Benchmark Specifications を含む） https://pages.nist.gov/pfhub/
• DeWitt, S. et al., PRISMS-PF: A general framework for phase-field modeling, npj Computational Materials 6, 29 (2020) https://www.nature.com/articles/s41524-020-0298-5
• 大野宗一, 凝固・結晶成長のためのフェーズフィールド法, まてりあ関連の解説（J-STAGE 掲載） https://www.jstage.jst.go.jp/article/ijmsa/30/1/30_24/_pdf
• PFHub Benchmarks Overview
  https://pages.nist.gov/pfhub/benchmarks/
• Jokisaari, A. M. et al., Benchmark problems for numerical implementations of phase field models, Computational Materials Science 126 (2017)
• Eyre の凸分割に基づくエネルギー安定化に関する文献（convex splitting / unconditionally energy stable schemes）
• PETSc documentation（MATSHELL / MatCreateShell など行列フリー関連）
• deal.II documentation（matrix-free infrastructure）
• Miyoshi, E. et al., Ultra-large-scale phase-field simulation study of ideal grain growth, npj Computational Materials 3, 12 (2017)
• Takaki, T. et al., GPU を用いた大規模凝固・結晶成長のフェーズフィールド計算に関する報告（関連論文・講演資料）
• N. Moelans, B. Blanpain, P. Wollants, An introduction to phase-field modeling of microstructure evolution, Calphad, 2008.
• PFHub, Phase Field Method Recommended Practices.
• MOOSE Framework, Phase Field Module documentation.
• S. DeWitt et al., PRISMS-PF: A general framework for phase-field modeling with a matrix-free finite element method, npj Computational Materials, 2020.
• 瀬川正仁ほか, Multi-Phase-Field Simulation of Non-Equilibrium Solidification Microstructure Formation in Additive Manufacturing, 2024.
• Y. Kamikawa et al., Chemo-electro-mechanical phase-field simulation of Li electrodeposition, Communications Materials, 2024.
• K. Alhada–Lahbabi et al., Machine learning surrogate for 3D phase-field modeling of ferroelectric tip-induced electrical switching, npj Computational Materials, 2024.
• HPCI成果報告, フェーズフィールド法と格子ボルツマン法の大規模計算によるデンドライト樹間液相流れの透過率評価, 2021.
• P. Munch et al., Efficient Finite-Element Solver for Phase-Field Simulations, 2024.
• T. Gong et al., Accelerating phase-field simulation of multi-component systems via machine learning, 2025.