Phase-Field-Crystal 法と分子動力学法の比較
京都工芸繊維大学

高木 知弘

Phase-Field-Crystal Method and Molecular Dynamics Method
Tomohiro TAKAKI
１．はじめにうな phase field 変数を用いる．今日では，PF 法および MPF

Phase-field-crystal (PFC) 法は，金属材料内の原子構造の規則法は，多くの物理・化学・工学的な問題に応用されているが，正しい配列を，周期的に分布する秩序変数（phase field 変数）(1)phase field 変数の考え方は基本的に同じである．

によって表現する新しいモデルである ．PFC 法は，分子動Phase-Field-Crystal (PFC) 法は，上記のような一般的な PF

力学 (MD) 法と同程度の空間スケールを対象とするが，非常法と phase field 変数の定義が異なる．図 1 (b)に PFC 法で用いに大きな時間増分をとることが可能であり，新しいマルチスる phase field φ の界面近傍におけるプロファイルを示していケール手法として期待されている．

る．液相内では通常の PF 法と同様，一定値φ = φ0 をとるが，PFC 法を変形特性評価に用いる場合，外力の与え方など境固相内では原子 1 つに対して phase field 変数の 1 周期が対応界条件の設定が困難となる．著者らは，差分格子間隔を変形する．に応じて変化させ，塑性変形中の体積一定変形を再現する手

２・２ 自由エネルギー汎関数と相図 以上のように定義さ法を構築した ．また，本手法を用いて，様々な条件の下でれたphase field φ を用いて自由エネルギー汎関数を構築する．のナノ多結晶体の変形シミュレーションを行うことで，PFC固相内におけるφ の周期構造によってエネルギーを最小化す

(2)

(3)

法が変形評価手法として有用であることを示した ．

るように，つまりφ の勾配が存在するとエネルギーを減らす

本稿では，PFC 法の基本的なメカニズムを概説し，広い温ように − ∇φ 2 を導入する．しかしながら，この項だけでは勾度範囲での変形シミュレーションを試みる．また，MD 法と配が無限大になってしまうため，1 つ高い次数を有するの比較についても考察を行う．

∇ 2φ を用いる必要がある．さらに，1 次相転位を表現するためにダブルウェル関数 fd (φ)を導入する．以上より，PFC 法

2の自由エネルギー汎関数は以下のように表すことができる．
22F = ∫ ⎛⎜ f d (φ ) − a3 ∇φ + a4 ∇ 2φ ⎞⎟dVV⎝⎠()= ∫ {f (φ ) + φ (a ∇ + a ∇ )φ }dV= ∫ f d (φ ) + a3φ∇ 2φ + a4φ∇ 4φ dV

(1)

V2V

d434

ここで，a3 と a4 は正の定数であり，式の変形過程において部分積分とガウスの発散定理を用いている．式 (1) の具体形として，最も基本的な形である次の Swift-Hohenberg モデルを用いる．
2 ⎤φ1⎡αF = ∫ ⎢ φ 2 + φ 4 + 1 + ∇ 2 φ ⎥ dVV 242⎣⎦()

ここで，α は温度に比例する係数である．
Fig. 1 Phase field variables (a) phase field for normal PF method, (b)phase field for PFC method and (c) atomic distribution aroundsolid-liquid interface
２．PFC 法２・１ Phase field 変数 一般的な PF 法で用いる phase field
変数φ は，図 1 (a)に示すように，ある相内で一定値を取り，界面において滑らかに変化する．例えば，凝固問題の場合，液相内でφ = 0，固相内でφ = 1 などと定義し，固液界面を 0 < φ< 1 の領域と考える．また，多相問題や多結晶粒成長問題では，phase field 変数を複数用いた Multi-phase-field (MPF) 法が多く用いられている．この場合，i 番目の相もしくは粒内においてφi = 1，他の相もしくは他の粒内においてφi = 0，i 番目の相もしくは粒の異相界面もしくは粒界において 0 < φi < 1 となるよ

Fig. 2 Phase diagram(2)式(2)より，φ の平均値φ 0 とα の関係として，図 2 のような相図を描くことができる．グレーの領域は各相の共存領域である．
２・３ 時間発展方程式 Phase field 変数φ は保存量であるので，Chan-Hilliard 方程式より次式の発展方程式を得ることができる．

[) ]2δF∂φ= ∇2= ∇ 2 αφ + φ 3 + 1 + ∇ 2 φ∂tδφ((3)

３．シミュレーション
３・１ α の影響評価 図 2 中に示す 3 種類のα を用いて単結晶の変形シミュレーションを行う．α は温度を表す係数でFig. 5 Stress-strain curve corresponding to the results of Fig. 4あり，負側に大きくなるほど温度は低下する．図 3 左図は用いた解析領域 (150×150 FD lattices)と初期 phase field プロフ
ァイルを示している．図 3 右図は，10% (ε = 0.1) まで x 方向に引張りシミュレーションを行った際の応力‐ひずみ曲線を示している．この場合，非常に小さい完全結晶に対して拘束の強い変形シミュレーションを行っているため，転位などの欠陥は生じず，ほぼ弾性的な応答を示した．しかしながら，

α の影響は非常に顕著に応力－ひずみ関係に現れていることがわかる．

Fig. 6 Atomic arrangement formed from Fig. 4 (a) for MD
simulation
ー分布として表している．また，図 5 は対応する応力－ひずみ曲線である．荷重を負荷すると，始め微視構造に大きな変化は見られず，応力は単調に増加する．ε = 0.06 あたりになると粒界を構成していた転位が上下の粒内にすべり面に沿って移動し，局部的な回転を誘発させることで，2 結晶間の方位Fig. 3 Phase field profile (left) and stress – strain curves (right)差を小さくしようと寄与する．さらに変形が加わると，転位は中央の粒内にも進入し，解析領域全体において単結晶化が３・２ 2 結晶変形シミュレーション 図 4 (a) に示す 2 結晶生じていることがわかる．図 5 の応力－ひずみ関係は，上記もでる(300×300 FD lattices)を用意し，x 方向への引張シミュのような転位の挙動を反映した結果となっている．弾性領域

レーションを行う．ここで，α = - 0.5 としている．図 4(b)(c)(d)のひずみが非常に大きいが，粒界からの転位の生成や，すべは，それぞれε = 0.06, 0.07, 0.08 時の変形組織を自由エネルギり面に沿った転位の移動や消滅など，転位挙動を良く再現できている．
図 6 は図 4(a)の初期 phase field プロファイルから作り出したMD シミュレーションで用いる原子構造である．この場合，phase field プロファイルの局大値に原子を置くことで原子構造を形成している．図 6 のモデルを用いて 2 次元 MD シミュ
レーションを試みたが，非常に早い段階で粒界にボイドが発生し，PFC シミュレーションと比較できるような良好な結果を得ることができなかった．詳細は講演会当日に報告する．
参考文献

1)

K. R. Elder, M. Katakowski, M. Haataja, M. Grant, Phys.
Rev. Lett. 88 (2002) 245701.

2)

T. Hirouchi, T. Takaki, Y. Tomita, Comp. Mater. Sci. 44(2009) 1192.

3)

T. Hirouchi, T. Takaki, Y. Tomita, submitted to Int. J. Mech.
Sci. (2009).

Fig. 4 Microstructures at ε = (a) 0.00, (b) 0.06, (c) 0.07 and (d) 0.08