<M&M2008 材料力学カンファレンス・2008 年 9 月 16 日～18 日>
Copyright © 2008 社団法人 日本機械学会

OS0410Phase-Field-Crystal 法の基本特性評価
Phase-Field-Crystal Method and Its Basic Performance○高木 知弘・京工繊大院

廣内 智之・神戸大院

冨田 佳宏・神戸大院

Tomohiro Takaki, Graduate School of Scienceand Technology, Kyoto Institute of Technology

Tomoyuki Hirouchi, Graduate School ofEngineering, Kobe University

Yoshihiro Tomita, Graduate School ofEngineering, Kobe University

Key Words: Phase-field-crystal method, Deformation, Phase diagram, one-dimensional problem
Phase-Field-Crystal (PFC) 法を材料変形シミュレーション手法として工学的なレベルで用いる
新たなモデルを構築するための基礎データを蓄積するために，最も基本的な PFC 法の導出および基本的な特性評価を行った．また，1 次元変形シミュレーションを行い数値シミュレーション精度について検討を行った．

1. はじめにPhase-field (PF) 法は複雑な材料組織発展を精度良く評価することを可能とする強力なツールとして発展している(1-4)．
一般的な PF 法では，相もしくは粒内では一定値を取り，界面において値が変化するような phase field と呼ばれる秩序変数を導入し，phase field の時間発展方程式を解くことで，組織変化を予測する．この場合，固体の原子構造に依存して発現する界面異方性や方位差依存性は，PF 方程式中の係数を方位や方位差の関数とすることで表現している．
古典的な PF 法の応用モデルとして，Phase-field-crystal(PFC) 法が Elder ら(5) によって提案されている．PFC 法は，固体内の規則正しい原子配列を周期系(periodic system)として捕らえ，原子一つ一つを秩序変数 phase field によって区別するモデルである．そのため，原子構造に起因する結晶方位，弾性・塑性特性，転位の生成・消滅などを自然に含んだモデルとなっている．さらに，発展方程式が拡散方程式であるため，原子スケールの構造を対象としながら非常にゆっくりした挙動の評価が可能である．このように，PFC 法は新しいマルチスケール手法として期待されている．しかしながら，これまで報告されている PFC 法を用いた研究では，仮定したパラメータを用いた仮想的な材料を対象としており，PFC 法を実用的な手法として用いるためには，PFC パラメータと物性値との関連付けが不可欠である．
本研究では最も基本的な PFC モデルの導出および基本的な特性を考察することで，工学的に実用性の高い PFC モデル構築のための基礎データの蓄積を行うことを目的とする．
2. Phase-Field-Crystal モデル
一般的な PF 法では，以下の手順により解析を進める(2)．
① 秩序変数の決定② 自由エネルギー汎関数の構築③ 時間発展方程式の導出④ パラメータと物性値の関係付け⑤ 離散化および数値シミュレーション
基本的に PFC 法も同様の手順を踏み，以下に本手順に沿って PFC モデルの導出を概説する．
PFC 法で用いる phase field 変数φ は，図 1(b)のように液相もしくは気相内で一定値を取り，固相内で格子サイズ a を波長とする周期構造を形成する．物理的には原子の存在確立もしくは原子の数密度を表すため，phase field φ は保存量である必要がある．

Fig. 1 Order parameter (a) atomic distribution around interface,(b) phase field for PFC method and (c) phase field for traditional
PF method
以上のように定義された phase field φ を用いて自由エネルギー汎関数を構築する．固相内におけるφ の周期構造によってエネルギーを最小化するように，つまりφ の勾配が存在するとエネルギーを減らすように−|∇φ |2 を導入する．しかしながら，この項だけでは勾配が無限大になってしまうため，1 つ高い次数を有する|∇2φ |2 を用いる必要がある．さらに，1次相転位を表現するために double well 関数を導入する．以上
より，PFC 法の自由エネルギー汎関数は以下のように表すことができる．
22F = ∫ ⎛⎜ f d (φ ) − a3 ∇φ + a 4 ∇ 2φ ⎞⎟dVV⎝⎠()= ∫ {f (φ ) + φ (a ∇ + a ∇ )φ }dV= ∫ f d (φ ) + a3φ∇ 2φ + a 4φ∇ 4φ dV(1)

V2Vd344

ここで，a3 と a4 は正の定数，fd(φ)は double well 関数であり，式の変形過程において部分積分とガウスの発散定理を用いている．式(1)の具体形として，最も基本的な形である次のSwift-Hohenberg モデルを用いる．

{)}

2u ⎤⎡φF = ∫ ⎢ r + λ q02 + ∇ 2 φ + φ 4 ⎥ dV
V 24 ⎦⎣((2)

次に，phase fieldφ の時間発展方程式を導出する．φ は保存料であるため，時間発展方程式は次式のように表される．

∂φδF= Mφ∇2∂tδφ[{(= M φ ∇ 2 r + λ q02 + ∇ 2) }φ + uφ ]2(3)3

ここで，Mφ はφ のモビリティーである．この式を差分法などによって離散化して，数値シミュレーションを行う．
ここで，課題であるパラメータ r，λ，q0，u，Mφと物性値の関係付けを行う必要がある．その基礎的な考察を行うために，1 次元問題を考え図 2 に示す相図を用いる(5)．ここで，λ= q0 = u = 1 としている．

ここで，φ0 = -0.2 を基準とし，±10%までの引張りと圧縮ひずみを作用させた場合のエネルギーF/a の変化を r = -0.2, -0.3,-0.4 に対して示している．図 3 の左の図は横軸がひずみε，右の図は横軸が φ0 である．ひずみ ε と φ0 の間には， φ0 =1/(1- ε)φ00 の関係がある．ここで，φ00 は基準状態でありφ00 =-0.2 である．図 3 の状態は図 2 の赤，青，緑の線分に対応している．図 3 より，引張りが生じるとエネルギーが単調に増加するが，圧縮が生じるとエネルギーが減少することが分かる．特に，図 3 左の図は原子間ポテンシャルに対応するため，ε = 0 においてエネルギーが最小となり，ひずみが増加しても減少してもエネルギーが増加するような傾向が好ましい．しかしながら，図 3 よりφ0 が減少するとエネルギーは単調に増加しており，引張り側しか適切に表現できないことが分かる．
次に，棒の引張りシミュレーションを行う．図 4 のようなphase field プロファイルを初期条件とし，差分格子間隔を徐々に増加させることで引張り変形をシミュレートする．ここで，原子間距離 a を 8, 10, 12, 14, 16 の差分格子で分割して精度を評価する．なお，r = -0.05，φ0 = 0 を用いている．図 5は得られた応力－ひずみ関係である．点線は理論解である．
図 5 より，精度の高い評価を行うためには，かなり細かい差分格子を必要とし，さらに変形が大きくなるほど差分格子間隔が広くなるため，精度が低下することが分かる．
より詳細な考察は講演当日に行う．

Fig. 2 Phase diagram for one dimension
図 2 の相図は以下のようにして導出する．液相もしくは気相の一定相φ = φ0 が安定な場合の式(2)のエネルギーは，
1⎡1⎤(4)
F = Funif = ⎢ r + λq04 φ02 + uφ04 ⎥ a4⎦⎣2とできる．また，固相の周期構造が安定な場合，φ の平衡プ
ロファイルをφ = Aeq sin(2π/a x) + φ0 とすると，同様に式(2)は⎡ r2 1⎤5(5)F = Fperi = ⎢−+ − r + λq04 φ02 − uφ04 ⎥ au624⎣⎦()()とできる．ここで，a は図 1 に示す格子サイズである．式(4)の Funif と式(5)の Fperi の大小関係により安定相が決定される．
図 2 中，水色の領域が Funif <Fperi で一定相が安定，黄色の領域が Funif >Fperi で周期相が安定な領域である．ここで，両相が共存可能な不安定平衡状態が存在するため，Maxwell の等面積原理よりその領域を導出したのが図 2 の赤色で示した領域である．
3. 考察
本一連の研究では，PFC 法を材料の変形解析手法として用いることを目標としており，図 2 の黄色い領域で示す周期相を考える．ここで，棒の引張り圧縮問題を考える．PFC モデ
ルの 1 次元問題で考えると，phase field φ の積分量は変形前後で保存されなければならないので，棒が引張りもしくは圧縮され棒の長さが変わるとφ0 は変化しなければならない．つまり，温度一定の状態での棒の引張り圧縮は，図 2 の相図上で r 一定でφ0 が増減することに対応する．この場合のエネル
ギー変化を図示すると，図 3 のようになる．

Fig. 3 Variations of energy F/a with changing φ0Fig. 4 Initial phase field profile for tensile simulation

Fig. 5 Stress – strain curves obtained by tensile simulations
参考文献[1] 小山敏幸, 材料組織設計工学, (社)日本鉄鋼協会フォーラ
ム資料, (2005).
[2] 高木知弘, 計算力学ハンドブック －14.1 フェーズフィールド法による相変化解析－, (2007), 376-389.
[3] K. G. G. Janssens, D. Raabe, E. Kozeschnik, M. A.
Miodownik : Computational Materials Engineering,
ELSEVIER (2007).
[4] N. Moelans, B. Blanpain, P. Wollants, CALPHAD, 32 (2008),
268.
[5] K. R. Elder, M. Grant, Phys. Rev. E, 70 (2004), 051605.