Int. J. Microgravity Sci. Appl. Vol. 30 No. 1 2013 (24–29)

IIIII 特集：理論化・モデリング（解説）その 3 IIIII

凝固・結晶成長のためのフェーズフィールド法大野宗一

Phase-Field Models for Solidification and Crystal Growth
Munekazu OHNOAbstract
The phase-field models developed in the field of solidification and crystal growth are explained with a focus on the models
for dendritic growth in pure substances and alloys. The diffuse interface description, free energy functional, and variational
derivation of phase-field equation with inclusion of anisotropy of solid-liquid interfacial energy are described. In addition,
the essential points of the Kim-Kim-Suzuki’s model and the quantitative phase-field model for alloy solidification are briefly
discussed.
Keyword(s): phase-field model, simulation, solidification microstructure, dendrite

1. はじめに本稿で紹介するフェーズフィールド法は，非平衡過程における多種多様なパターン形成を記述する汎用的かつ強力な手法として発展しており，現在，この手法が応用される分野は多岐に亘ります．この手法は，固相―固相変態の分野においては Cahn-Hilliard 方程式 1)と AllenCahn 方程式 2)に基礎を置いた物理モデルとして発展していますが，凝固・結晶成長の分野に関して言えば，diffuse interface に基づいて界面移動問題を効率的に解くための数値計算手法として，現在までの発展を遂げています．従いまして，そのモデリングの philosophy も必然的に異なります．凝固・結晶成長の分野におけるフェーズフィールド・モデリングの主眼は，組織レベルのスケールにおける固液界面の移動を如何に効率的，そして適切に記述するのかという点にあるといえます．
凝固中の固液界面の移動には，(i) Gibbs-Thomson 則，(ii)固液界面の移動に伴う溶質濃度分配とエネルギー変化に対する保存則(Stefan 条件)，(iii)固相と液相中の溶質・熱拡散（バルク中の拡散），の三種類の法則が関与します
ので，これら三種類の法則に基づいて界面の移動を記述する必要があります．これらの法則を同時に成立させながら界面の移動を記述する問題とは，つまり自由境界問
題です．凝固・結晶成長のフェーズフィールド法は，自由境界問題を効率的に解く手法として発展してきました．
組成的過冷却理論，Mullins-Sekerka の摂動論 3)，KGTモデル 4)や LKT モデル 5)など，凝固工学・科学で発展してきた理論の多くが，(i)から(iii)の法則（の一部）に基礎を置いていることを考えますと，(i)から(iii)に立脚して複雑形態の組織を再現するフェーズフィールド法は，この分野において極めて重要な役割を担う手法であるといえます．
フェーズフィールド法に関する解説は既に多数報告されており 6-9)，いくつかの書籍も出版されています．本稿
では，凝固・結晶成長において発展してきたフェーズフィールド法について概説します．現在までに考案されたフェーズフィールド法は複数あり，それぞれ改良が試みられてきました．本稿では，それらの発展の経緯や応用範囲を俯瞰するようなことはせずに，それらのモデルの中で汎用性・計算効率・計算精度の観点から最も優れていると考えられるモデルを中心に，その計算に必要な方程式を簡潔に説明します．まず純物質の凝固を対象として，この手法が立脚する diffuse interface について述べ，自由エネルギー汎関数，時間発展方程式や界面エネルギーの異方性の取り扱いについて説明したいと思います．
その後，合金の凝固過程を記述するモデルについて紹介します．

北海道大学 大学院工学研究院 材料科学部門 〒060-8628 札幌市北区北 13 条西 8 丁目
Division of Materials Science and Engineering, Faculty of Engineering, Hokkaido University,Sapporo, Hokkaido 060-8628, Japan(E-mail: mohno@eng.hokudai.ac.jp)− 24 −

Kita 13, Nishi 8, Kita-ku,

24

大野宗一

2. 純物質のデンドライト成長

Dendrite structure

2.1 diffuse interface
フェーズフィールド法においてまず特筆すべき点は，界面を diffuse interface として表現することです．この表現法のおかげで，組織が複雑になっても計算が煩雑になるようなことはありません．その diffuse interface について簡単に説明します．
凝固における組織形成を記述・予測するためには，一般に，時々刻々と変化する固液界面の位置を追跡する必要があります．しかし，デンドライト形状のように組織の形態が複雑になりますと，その組織における全ての固液界面の位置を記録し，その時間変化を追跡するのは大変困難になります．フェーズフィールド法の最大の特徴
は，界面の存在を陽に取り扱わないことです．この手法
では，相の種類を秩序変数によって区別します(Fig. 1 参照)．その秩序変数が空間的に均一（値が一定）の領域をバルク相の領域と定義し，その変数が空間的に変化する領域を異相界面とみなします．本稿では，秩序変数を で表記し，=0 のときに液相，=1 のときに固相であると定義します．の値が 0 と 1 の間で変化する領域が固液界面に相当します．ここで重要なことは，固液界面の移動を計算するために界面の空間座標を特定し，その座標の時間変化を記述するする必要がないということです．の空間プロファイルの時間変化を記述することで，界面移動が計算されます．言い換えますと，システムをメッシュに分割し，その各格子点上での値の時間変化を計算することで，組織形成過程が計算されることになります．これがフェーズフィールド法における組織の表現法で，デンドライトのような複雑形状の組織を簡便にシミュレートすることを可能にする重要な点です．
固液界面を対象とした原子レベルのシミュレーションによれば 10)，現実の固液界面の幅は数Å～数 nm 程度です．フェーズフィールド法の界面領域においては，が連続的に変化する必要がありますから，計算に用いる空間グリッドサイズはその界面幅よりも小さくする必要があります．従いまして，ある方向に，例えば 1000 メッシュを割り当てたとしても，その全長はせいぜい 1m 程度です．一方で，三次元の全ての方向に 1000 メッシュを割り当てると 109 メッシュの計算になります．これは現在の最新のスーパーコンピュータを必要とする計算規模です．
つまり，超大規模計算を実行できたとしても，その計算が対象とできるシステムのサイズは非常に小さなものです．従いまして，典型的な凝固組織のスケールである数十～数百m スケールのシステムを対象とするためには，空間グリッドサイズを大きくする必要があります．この手法において，界面は常に有限の幅を持ちますから，空間グリッドサイズを大きくすると必然的に界面幅も大きくなります．つまり，フェーズフィールド法における界面幅とは現実の界面幅を意味するわけではありません．

25Int. J. Microgravity Sci. Appl. Vol. 30 No. 1 2013

Diffuse interface+1 Solid0Liquid

Spatial coordinate, xFig. 1 Diffuse interface description
この手法において界面幅は任意の定数で，通常よく用いられる界面幅は，現実の値よりも数桁大きいものです．
従いまして，界面幅を変化させても界面エネルギーなどの物性を変化させず，計算結果も変化しないようにモデルを構築することが必要となります．この辺の詳細については後ほど議論したいと思います．

2.2 自由エネルギーと発展方程式
フェーズフィールド法においては，を含めた状態変数の時間発展を，システムの自由エネルギーが時間と共に減少するという原理に立脚して計算します．ここで，フェーズフィールド法において対象とする自由エネルギーとは，不均一系の自由エネルギーです．状態変数は空間座標（さらに時間 t）に依存する関数  (r,t)であり，自由エネルギーは関数 (r,t)の関数，つまり汎関数で与えられます．今，状態変数として のみを考えると，不均一系の
自 由 エ ネ ル ギ ー 汎 関 数 ， 言 い 換 え ま す と ， GinzburgLandau 型の自由エネルギーは以下のように与えられます．
1𝐹 = ∫ ( 𝜀 2 |∇ϕ|2 + 𝑓(𝜙)) 𝑑𝑣
2(1)

ここで，積分は体積積分に相当します．被積分項の第一項は勾配エネルギーと呼ばれ，はこの大きさを決定する係数，勾配エネルギー係数です．その形式から明らかなように，が空間的に非一様なときには，この第一項のエネルギーが大きくなります．
式(1)の被積分項の第二項は，系が均一のときの自由エネルギー変化量を表しており，一般に二つの項からなります．
𝑓(𝜙) = 𝑓𝑑𝑤 (𝜙) + 𝑓𝑏𝑢𝑙𝑘 (𝜙)(2)

ここで，fdw()は固相と液相の自由エネルギー障壁を表した項であり，この議論の中では，以下の二重井戸型ポテンシャルで近似します．

− 25 −

凝固・結晶成長のためのフェーズフィールド法𝑓𝑑𝑤 (𝜙) = 𝜔𝜙 2 (1 − 𝜙)2(3)

ここで，  は定数です．式(3)の形から明らかなように，fdw()は =0，+1 で最小値をとる関数です．つまり，系が固相と液相以外の状態になったとき，この項のエネルギーが大きくなります．
式(2)における fbulk()はバルクの自由エネルギーの寄与です．ここで，熱力学において相平衡を考える際に用いる自由エネルギーを，均一系自由エネルギーと呼び，fs とfl を，それぞれ固相と液相の均一系自由エネルギーとして定義すると，fbulk()は fbulk(=0)=fl，fbulk(=1)=fs となる関数です．これは，g(=0)=0，g(=1)=1 を満足するの単調増加関数 g()を用いて，𝑓𝑏𝑢𝑙𝑘 (𝜙) = 𝑔(𝜙)𝑓𝑠 + (1 − 𝑔(𝜙))𝑓𝑙(4)という形で近似されます．また，に関する g()の一次導関数 g’()が，g’(=0)=g’(=1)=0 の条件を満足するようにg()を定義する必要があります．ここでは，一般的に用いられる関数の一つとして，𝑔(𝜙) = 𝜙 3 (10 − 15𝜙 + 6𝜙 2 )を用いることにします．
既に述べたように，組織の時間発展は，システムの自由エネルギーが時間と共に常に減少するという原理に立脚して計算されます．組織の時間発展とは，状態変数の空間プロファイルの時間発展に相当します．状態変数が保存量，非保存量であるかによって時間発展方程式の形式は異なります．ここで対象とする秩序変数 は非保存量です．その時間発展は以下の方程式で記述されます．
𝜕𝜙𝛿𝐹= −𝑀𝜙𝜕𝑡𝛿𝜙(5)

ここで，Mは定数です．式(5)の右辺は汎関数微分を意味しています．式(1)を式(5)に代入することで， の時間変化が計算されます．
1 𝜕𝜙𝑑𝑓𝑑𝑤 (𝜙) 𝑑𝑔(𝜙)(𝑓𝑠 − 𝑓𝑙 )= 𝜀 2 ∇2 𝜙 −−𝑀𝜙 𝜕𝑡𝑑𝜙𝑑𝜙(7)ここで，()は熱拡散率です．固相と液相中の熱拡散率をs ， l とすると， (=1)=s ， (=0)=l であり，これもfbulk()と同様に の単調増加関数を使って近似します．右辺第二項は潜熱の効果を表しており，L は単位体積当たりInt. J. Microgravity Sci. Appl. Vol. 30 No. 1 2013

2.3

物性値

ここでは，平衡状態における平滑な固液界面を対象とします．従いまして，式(6)において∂𝜙⁄∂𝑡 = 0 ，さらにfs=fl となります．界面と垂直な方向を x 方向とし，x→+∞で=0（液相），x→∞で =1（固相）の境界条件を考えると，式(6)の解として以下の方程式が得られます．
1𝑥𝜙(𝑥) = [1 − tanh ()]2√2𝑊(8)

ここで，W は，W=/1/2 です．W が大きいと，(x)のプロファイルは緩やかになり，小さいと急峻になることが式(8)からわかります．つまり，界面幅は W の大小で決められます．ここでは W を界面幅であるとして議論を進めます．式(8)を式(1)に代入すると，界面が存在することによる過剰なエネルギー，すなわち固液界面エネルギー が以下のように求まります．
𝛾=(6)上式の右辺第三項が凝固の駆動力に相当します．ここで，明記していませんが，fs 及び fl は温度に依存しています．
物質が凝固する際には潜熱が放出されますので，温度場の変化を以下の熱伝導方程式で計算します．
𝜕𝑇
𝐿 𝑑ℎ(𝜙)= ∇(𝛼(ϕ)∇𝑇) +𝜕𝑡
𝐶𝑝 𝑑𝑡の潜熱，Cp は単位体積当たりの比熱です．h()は の単調増加関数であり，本稿では h()=と定義します．なお，式(7)の熱伝導方程式は，エントロピー汎関数から導出することが可能です 11)．また，式(6)と同様に自由エネルギー汎関数から導出する方法も示されています 12)．式(6)における右辺第三項の駆動力 fsfl の温度依存性を介して，温度場と のダイナミックスがカップリングすることになります．
式(6)と(7)を連立して解くことで，純物質の凝固過程が計算されます．その計算には， ， ，M，fs，fl ， s， l，L，そして Cp の値を入力する必要があります．これらの中で，， ，Mはフェーズフィールド法に固有の物理量ですので，一般に文献値は存在しません．さらに，フェーズフィールド法においては界面幅を任意の定数として扱うので，界面幅も入力する必要があります．これらの物理量をどのように与えるのかについて次節で説明しま
す．

13√2𝜀√𝜔(9)このように，とは，固液界面幅 W と固液界面エネルギーと関連付けられます．従って，の測定値または計算値，そして W の値を決めると，との値はそれに応じて必然的に決定されることになります．
Mも同様に測定可能な物性値と関連付けられます．ただし，この議論のためには，方程式の漸近解析，sharpinterface limit および thin-interface limit13)に関する詳細な説明を要します．紙面の制約のため，それらの詳細を割愛して，ここでは，thin-interface limit から得られる結果のみを以下に記します．

− 26 −26

大野

𝑀𝜙 = 𝐶𝜔𝐶𝑝 𝑇𝑚 𝛼𝜀 2 𝐿2宗一𝛿𝐹 𝑑𝑓𝑑𝑤 (𝜙) 𝑑𝑔(𝜙)(𝑓𝑠 − 𝑓𝑙 ) − 𝜀 2 ∇2 𝜙=+𝛿𝜙𝑑𝜙𝑑𝜙
𝜕2𝜙 𝜕2𝜙− 𝜀ε′ [sin(2𝜃) ( 2 − 2 )𝜕𝑦𝜕𝑥2
𝜕 𝜙+ 2cos(2𝜃)]𝜕𝑥𝜕𝑦1
𝜕2𝜙+ (𝜀′2 + 𝜀𝜀") [2sin(2𝜃)2𝜕𝑥𝜕𝑦2−∇ 𝜙𝜕2𝜙 𝜕2𝜙− cos(2𝜃) ( 2 − 2 )]𝜕𝑦𝜕𝑥(10)C は定数で，その値は式(6)中の fdw()，g()，式(7)中の()，h()の形式に依存します．また，この式は =s=lとして導出したものです 13)．

2.4異方性

既に述べましたように，式(6)と(7)を連立して解くことで，凝固過程を記述することができます．ただし，このモデルは界面エネルギーが等方的であると仮定したものです．デンドライトを再現するためには，界面エネルギーの異方性，つまり方位依存性を考慮する必要があります．例えば，4 回対称の異方性は，以下の式によって近似されます．
𝛾(𝜃) = 𝛾0 (1 + 𝜀4 cos(4𝜃))(11)

ここで，は界面の法線方向とシステムの x 軸方向とのなす角度であり，0 及び4 は定数です．4 が異方性の強さを表します．式(9)から理解できるように， は に比例します．従って，式(11)に示した の異方性を考慮するために，に同様の異方性を導入することになります．
𝜀(𝜃) = 𝜀0 (1 + 𝜀4 cos(4𝜃))𝜃 = tan−1(𝜕𝜙⁄𝜕𝑦)(𝜕𝜙⁄𝜕𝑥)(13)

従って， が プロファイルに依存します．このことを考慮すると，式(5)の汎関数微分は以下のように与えられます．
𝛿𝐹 𝑑𝑓𝑑𝑤 (𝜙) 𝑑𝑔(𝜙)(𝑓𝑠 − 𝑓𝑙 ) − ∇ ∙ (𝜀 2 ∇𝜙)=+𝛿𝜙𝑑𝜙𝑑𝜙𝜕𝜕𝜙𝜕
𝜕𝜙+(𝜀𝜀′ ) −(𝜀𝜀′ )𝜕𝑥𝜕𝑦𝜕𝑦𝜕𝑥

式(15)には， ()の二回微分が表れています．式(11)に示したような 4 回対称の異方性を導入する場合には，の二次導関数は，4 の 16 倍の異方性を示します．従いまして，式(11)において，界面エネルギーが方位によって，高々数%しか違わないとしても，組織形態には決定的な違いが現れます．実際，金属における界面エネルギーの大きさは，方向によって数%程度しか違わないことが多々ありますが，この数%の異方性を如何に正確に見積もるのかが，デンドライト組織の計算の精度にとって重要になります．

3. 合金凝固のデンドライト成長(12)

ここで， 0 は定数です． は を使って，以下のように与えられます．

(14)

ここで， ’=d/d です．式(14)を式(5)に代入して， の時間発展を解くことで，デンドライト組織が再現されます．
その際，式(10)に示したように，Mはに依存しますので，M にも  依存性が考慮されることになります．また，式(14)は，計算上，より便宜的な形式である式(15)の形に変形することができます．

(15)

3.1

合金凝固の駆動力二元系合金を対象として，合金凝固のモデルを紹介します．前節で説明したように，純物質の凝固の場合は，秩序変数と温度 T のみがカップリングすることを考慮してモデリングを行います．合金の凝固を対象とする場合には，これらに加えて，溶質濃度を変数として導入する必要があります．濃度の時間発展方程式は，式(1)に示したような Ginzburg-Landau 型の自由エネルギー汎関数を，濃 度 場 に 依 存 す る 形 式 で 再 定 義 し ， 最 終 的 に CahnHilliard 方程式の形で求めることができます．しかし，凝固の自由境界問題を記述するにあたって，濃度に関する勾配エネルギーを導入する必要はありません．現在最も良く用いられるモデルにおいては，式(4)の fbulk()における固相と液相の均一系自由エネルギーを，それぞれ濃度の関数にすることで濃度とのカップリングを導入します．つまり，二元系合金の溶質濃度を c で表すと，fs=fs(c,T)および fl=fl(c,T)と定義することになります．
ここで，濃度の時間発展方程式を求める前に，の時間発展方程式について考えます．上で述べたような濃度依存性を導入した自由エネルギーを，式(5)に代入すると，式(6)と同じ形式の方程式が得られます．ただし，式(6)の右 辺 第 三 項 の 駆 動 力 の 項 は ， 𝑑𝑔(𝜙)⁄𝑑𝜙 (𝑓𝑠 (𝑐, 𝑇) −𝑓𝑙 (𝑐, 𝑇))に書き直され，濃度に依存することになります．
ここで注意を要することは，純物質においては fs=fl が平衡状態の条件に相当しますが，合金においてこの条件は二相平衡の条件に相当しないということです．合金中の27Int. J. Microgravity Sci. Appl. Vol. 30 No. 1 2013− 27 −

凝固・結晶成長のためのフェーズフィールド法

二相平衡は，両相の化学ポテンシャルが等しいときに成立します．従いまして，合金の二相平衡状態において，式(6)の第三項の寄与は消滅しません．これは言い換えますと，平衡状態，さらには定常状態における のプロファイルに第三項の寄与が現れることを意味します．従って，定常状態の の空間プロファイルは式(8)と異なり，濃度に依存することになります．
純物質のときには，式(8)に示したように，温度場とのカップリングなしに のプロファイルが決定されるため，界面幅 W は予め決定することができます．しかし，合金の凝固の場合，上記の手法を採用すると，の空間プロフ
ァイルが濃度プロファイルに依存するために，W の値を予め決定することはできません．既に述べたように，フェーズフィールド法においては，diffuse interface を表現するため，界面幅よりも空間グリッドサイズを小さくする必要があります．小さな界面幅領域に対しては，それ以上に小さな空間グリッドサイズを割り当てる必要があります．従って，界面幅が濃度プロファイルに依存すると，計算効率が濃度プロファイルに依存することになり，上記の方法では計算効率が極めて悪くなります．この問題を解決するいくつかの方法が提案されていますが，最も汎用性があり，物理的にも見通しの良い方法が Kim,Kim, Suzuki らによって提案されています 14)．そのモデ
ルは KKS モデルと呼ばれており，次節においてそのKKS モデルのエッセンスを説明します．

3.2

合金凝固の駆動力と KKS モデル

KKS モデルにおいては，固相領域，液相領域に関わらず，システム全体にわたって，固相濃度 cs，液相濃度 clを定義します．そして，合金の濃度 c は
𝑐 = ℎ(𝜙)𝑐𝑠 + (1 − ℎ(𝜙))𝑐𝑙∂𝑓𝑠 ∂𝑓𝑙=∂𝑐𝑠 ∂𝑐𝑙(17)つまり，固相と液相の化学ポテンシャルが各空間点で等しいと仮定します．この条件によって，cs と cl は独立ではなくなり，濃度に関する変数は一つになります．希薄固溶体を対象とすると，式(17)は cs=kcl の条件と等価になります．
式(16)と(17)を考慮して，式(5)の自由エネルギーの汎関数微分を行うと，式(6)のの時間発展方程式は，1 𝜕𝜙𝑑𝑓𝑑𝑤 (𝜙)= 𝜀 2 ∇2 𝜙 −𝑀𝜙 𝜕𝑡𝑑𝜙(18)𝑑𝑔(𝜙)−(𝑓𝑠 − 𝑓𝑙 − 𝜇𝑐 (𝑐𝑠 − 𝑐𝑙 ))𝑑𝜙と書き表されます．第三項の駆動力が修正されており，式(18)の駆動力は二相平衡において 0 になります．この修正によって，定常状態における の空間プロファイルは，式(8)で与えられることになり，界面幅に濃度依存性が現れません．以上が KKS モデルにおける重要な特徴です．
濃度の時間発展方程式を導出する際にも式(16)と(17)を考慮して，以下の式が導かれます．
∂𝑐= ∇𝐷(𝜙)∇𝑐 + ∇[𝐷(𝜙)(𝑐𝑙 − 𝑐𝑠 )∇𝜙]∂𝑡(19)

ここで，D()は拡散係数です．
式(18)に式(14)や(15)のように異方性を導入し，式(7)と式(19)と連立して解くことで，合金の凝固におけるデンドライト組織が計算されます．

3.3(16)で与えられるとします．h()は の単調増加関数であり，ここでは h()= とします．従って，固相領域(=1)ではc=cs，液相領域( =0)では c=cl です．ここで，濃度の変数が一つ増えていることに注意して下さい．簡単のために，温度 T は一定であるとすると（等温凝固），溶質濃度 c と秩序変数 が決まれば，システムの熱力学状態は決定されるべきです．しかし，このモデルにおいては，cs と cl が新たに導入されています．式(16)の関係を考慮すると，システムの状態は秩序変数 と二つの濃度によって記述されます．従って，新たな条件を追加して，濃度に関する変数を一つにする必要があります．例えば，平衡分配係数 kを用いて，各空間点で cs = kcl という条件が成立することを仮定すると，濃度に関する変数（独立変数）は c, cs, clの中の一つとなります．実際，このように分配関数を用いたモデルも提案されていますが，KKS モデルでは以下の関係が成立すると仮定します．

Int. J. Microgravity Sci. Appl. Vol. 30 No. 1 2013𝜇𝑐 =定量的フェーズフィールド法KKS モデルは現在までに広く用いられ，その有効性が示されてきました．しかし，KKS モデルを含めた既存の合金凝固のモデルには，その計算結果の定量性に欠陥があることが知られています．詳細は割愛しますが，式(7)，(18) ， (19) の 方 程 式 の 解 の 挙 動 を 解 析 す る と ， GibbsThomson 則と Stefan 条件が正確に成立しないことが示されます．既に少し触れましたが，これは thin interfacelimit に基づいてモデリングが行われていないためであり，既存のモデルの計算結果は任意の定数である界面幅 W の値に依存します．言い換えますと，界面エネルギーやバルクの自由エネルギーなど，計算に入力する物性値を同じ値に設定しても，W の値によって組織のサイズ，形態，成長速度，濃度プロファイルなどが異なってしまいます．
この問題を解決した手法が定量的フェーズフィールド法であり，近年急速に発展・応用されています．
簡単のため，希薄固溶体を対象とすると，式(19)の拡散方程式は，以下のように書き表されます．

− 28 −28大野∂𝑐= ∇[𝐷𝑎 (𝜙)∇𝑐𝑙 ]∂𝑡(20)

ここで，Da()は拡散係数に関連する量です．固相と液相の拡散係数をそれぞれ Ds と Dl と表して，固相濃度と液相 濃 度 の 分 配 係 数 を k と 定 義 す る と ， Da(=1)=kDs ，Da(=0)=Dl を満足します．定量的フェーズフィールド法においては，上記の問題を解決するために，式(20)に新たな流束項を加えます．
∂𝑐1 ε∂𝜙 ∇𝜙(𝑐𝑙 − 𝑐𝑠 )= ∇ [𝐷𝑎 (𝜙)∇𝑐𝑙 +]∂𝑡∂𝑡 |∇𝜙|2√2 √ω

(21)

右辺括弧内の第二項が新たな流束項で，anti-trappingcurrent と呼ばれます 15,16)．この項を導入することで，上記の定量性に関する問題が解決されます．特に，界面幅W を大きくしても，高精度の計算が可能です．従って，大きな空間グリッドサイズを用いて，正しい計算ができるため，計算時間も著しく短縮できます．このモデルは凝固組織のシミュレーション法の主流になりつつあります．ただし，このモデルはいまだ発展途上にあり，扱えるシステム・現象は限定されています．式(21)は，固相中の拡散が存在せず(Ds=0)，希薄固溶体近似が成立する二元系合金にのみ適用可能です．著者らは，近年，この定量的フェーズフィールド法を固相拡散が存在する単相凝固17)，そして多相凝固 18)の計算に拡張することを試みました．Fig. 2 にその結果の例を示します．詳細は割愛します
が，二元系合金における一方向凝固過程を計算した結果であり，形状・濃度プロファイルなどが極めて高精度に算出されます．なお，デンドライトの二次枝発生のためには，界面形状・濃度・温度の揺らぎを考慮する必要があります．式(7)，(18)，(21)のそれぞれに揺動散逸定理を満たすノイズ項を導入することで揺らぎを記述します．
定量的モデルの最近の進展については，原著論文 15-18)を参考にしてください．

4. おわりにt’=16のカップリングや，多元系合金の取り扱い，多相凝固へ拡張など，この稿で説明しきれなかった内容がいくつかあります．これらの詳細については既に出版されている書籍・解説 6-9)などを参考にして下さい．
フェーズフィールド法の基礎方程式は古くから研究されてきましたが，組織形成のシミュレーション手法としては，ここ 20 年あまりで急速に発展してきた手法であり，現在においてもその応用範囲は拡張されつつあります．
その中で，凝固のフェーズフィールド法は，比較的歴史が古く，この稿の中でも説明したように，様々な改良が施されてきました．しかし，未だに解決されていないモデル上の欠陥が残されており，この先，新たな欠陥が見つかる可能性も十分にあります．この手法における発展方程式の数学的構造は単純ではなく，その数学的表現が物理現象を正しく捉えているのかについては常に注意を要します．しかし，フェーズフィールド法の概念は比較的理解しやすく，一から計算コードを書きあげるのも決して困難ではありません．本稿がフェーズフィールド法を始めるきっかけの一つになれば幸甚です．

参考文献
1)
2)
3)
4)
5)
6)
7)
8)
9)
10)
11)

本稿では，純物質と合金の凝固過程を記述するフェー
ズフィールド法について概説しました．流体の速度場とt’=14

宗一

t’=18

t’=20

12)
13)
14)
15)
16)
17)
18)

Fig. 2 Directional dendrite growth in a binary alloy

29Int. J. Microgravity Sci. Appl. Vol. 30 No. 1 2013

J.W. Cahn and J.E. Hilliard: J. Chem. Phys., 28 (1958)
258.
S.M. Allen and J.W. Cahn: Acta Metall., 27 (1979) 1085.
W.W. Mullins and R.F. Sekerka: J. App. Phys., 35 (1964)
444.
W. Kurz, B. Giovanola and R. Trivedi.: Acta Metall., 34(1986) 823.
J. Lipton, W. Kurz and R. Trivedi: Acta Metall., 35 (1987)
957.
T. Suzuki, S.G. Kim and W.T. Kim: JASMA, 19 (2002) 2(in Japanese).
T. Koyama: J. Japan Inst. Metals, 73 (2009) 891 (inJapanese).
R. Kobayashi, T. Takaki, T. Koyama, M. Ohno, K.
Matsuura and A. Takezawa: J. Japan Soc. Comp. Eng.
Sci., 15 (2010) 2287 (in Japanese).
T. Takaki and A. Yamanaka: Science of machine,60(2009)-63(2010) (in Japanese).
J.J. Hoyt, M. Asta and A. Karma: Mater. Sci. Eng. R, 41(2003) 121.
S-L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R.
Coriell, R.J. Braun and G.B. McFadden: Physica D, 69(1993) 189.
M. Plapp: Phys. Rev. E, 84 (2011) 031601.
A. Karma and W.-J. Rappel: Phys. Rev. E, 57 (1998) 4323.
S.G. Kim, W.T. Kim and T. Suzuki: Phys. Rev. E, 60(1999) 7186.
A. Karma: Phys. Rev. Lett., 87 (2001) 115701.
B. Echebarria, R. Folch, A. Karma and M. Plapp: Phys.
Rev. E, 70 (2004) 061604.
M. Ohno and K. Matsuura: Phys. Rev. E, 79 (2009)
031603.
M. Ohno and K. Matsuura: Acta Mater., 58 (2010) 5479.

(2012 年 10 月 31 日受理，2012 年 12 月 10 日採録)

− 29 −