Phase Field Crystal 法による変形シミュレーション
神戸大学 [院] ○廣内 智之

京都工芸繊維大学 [院] 高木 知弘
冨田 佳宏
神戸大学 [院]Deformation Simulation by Phase Field Crystal Method
Tomoyuki HIROUCHI, Tomohiro TAKAKI and Yoshihiro TOMITA

1緒

言

Initial state ∆ x = ∆ yDeformed state ∆ x' = ∆ y'

工業用材料の機能特性は, 材料内部の微視組織に存在すsAnalysi
domainる欠陥に強く影響されることが知られている. そのため,材料の更なる高機能化および新機能発現のためには, 微視組織に依存する巨視的な材料特性評価を可能とする数理yモデルの構築が必要不可欠である. 現在, 材料の力学的特∆yx

SS∆ x'∆xS = ∆ x∆ y= ∆ x'∆ y'∆ y'

性の評価手法として, 有限要素法と分子動力学法が広く用

Fig. 1: Proposed numerical scheme of tensile deforma-

いられているが, これらが対象としている時空間スケー

tion in two dimensions.

ルには大きな隔たりがある. このため, 原子配列に起因す力学を繋ぐ新たなマルチスケール手法として期待されている. 本手法は, これまでに様々な物理現象に適用されてDistance between atoms

60a (600 lattices)
Periodic boundary conditions in all direction

手法としての可能性を検討する必要がある.
そこで, 本研究では, 多様な結晶構造を対象として変形シミュレーションを行い, PFC 法の変形解析手法としての可能性を検討・拡張し, 実用的な変形解析手法として確立することを目的としている.

2101020304050Position(b)

Fig. 2: (a) Phase ﬁeld proﬁle at the initial condition ofsingle-crystal and (b) atomic displacement distributions.
2·2
体積一定変形の解析手法
これまでの研究成(2)
果において, PFC 法を用いて金属材料の体積一定で進行する塑性変形が再現できないことを示した. その問題点を解決するために, 差分格子の大きさを制御する新たなPhase Field Crystal 法

変形解析手法を提案する. 変形解析における差分格子の2·1
秩序変数と支配方程式
PFC 法は, 単結晶中の原子配列のような周期状態を表現するために, 秩序変数が周期的な値をとることが大きな特徴であり, 従来のPhase ﬁeld 法と大きく異なる点である. 本手法にて用いる秩序変数 phase ﬁeld φ は, 時間平均化された原子の数密度と定義され, 保存場である

0.002
0.004
0.0062(a)きたが, 多結晶構造の変形挙動を詳細に解析した研究報告はないようである. そのため, 本手法の一般的な変形解析

Strain0.42Displacement

Phase Field Crystal(PFC) 法 (1) は, 拡散時間スケールで個々の原子を追跡できることから, 分子動力学法と連続体

φ -0.61

3a (30 lattices)する数値解析手法の構築が急務である. 最近開発された

3Tensile direction

る材料挙動評価を実用的な時空間スケールで評価可能と(1). 固相領域においては,原子が時間に依存せず周期的な配置をとる. その周期状態を表現するために φ は周期的な値をとる. また, φ の最大点を原子の位置とする. 液相領域においては, ランダムな原子配置を時間的に平均化して一定値 φ0 とする.

分格子の x 軸方向の長さを d(= ε̇∆x∆t) だけ増加させて∆x とし, 常に差分格子の面積 S(= ∆x∆y = ∆x ∆y  ) が一定となるように y 軸方向の長さを決定する. ここで, ε̇は無次元ひずみ速度, ∆x, ∆y は x, y 軸方向の初期差分格子間隔, ∆x , ∆y  は変形過程での x, y 軸方向の差分格子間隔である. 本手法は, 全ての原子に一定のひずみ速度を加えるため, 変形状態はアフィン変形状態に相当する.

3変形シミュレーション

本章では, 提案した手法を用いて, 単結晶体, 二結晶体,ここで，この秩序変数を用いた無次元化された自由エネルギー汎関数 (1) を以下に示す．
 
α 2 φ4φ2 2F =φ ++ (1 + ∇ ) φ dV242制御法を模式的に図 1 に示す. 各時間ステップ ∆t 毎に差そして, ナノ多結晶体の変形シミュレーションを行い, 弾性および塑性変形挙動を評価する.

(1)

3·1単結晶体のシミュレーション
図 2(a) に単結晶体の phase ﬁeld の初期分布を示す. 変形前後の原子のここで, α は温度に関係する駆動力項で, 過冷度に比例す位置を容易に比較できるように, 引張方向に対して平行る係数である. 秩序変数を保存場とするため, φ の時間発に原子を配置している. 解析領域の縦横比の関係上, 中央展方程式は次のように表すことができる.

部分は省略している. 実際のひずみ速度でおよそ 100 /s〜∂φ= ∇2 [αφ + φ3 + (1 + ∇2 )2 φ]∂t(2)102 /s に相当する複数の変形速度で変形解析を行う.
図 2(a) に各変形量に対する原子の変位分布を示す. 点

Free energy 0.021

Periodic boundary conditions in all direction

0.041

60a (600 lattices)

44.8
1.64

B

41.7
57.4

19.120deg.

40deg.

Approximately
14atoms

20deg.

AA

35.2

41.9
4.69
18.6
44.8(a) ε = 0.00(b) ε = 0.02

bicrystal with ’stepped’ grain boundaries.

60a (600 lattices)

Periodic boundary conditions in all direction

Free energy 0.021

21.2
41.7

50a (500 lattices)(c) ε = 0.04

Fig. 3: Initial structure and deformation behaviors of

35.2

Crystal orientation

0.32
21.2

60a (600 lattices)

1.64
15.9

44a (440 lattices)

Periodic boundary conditions in all direction

Fig. 5: Desired nanopolycrystalline structure with regular hexagonal grains.

0.041

DC
20deg.

40deg.

20deg.

(b) ε = 0.05(a) ε = 0.00

Free energy 0.021

60a (600 lattices)(a) ε = 0.000(b) ε = 0.080(c) ε = 0.085

Fig. 4: Initial structure and deformation behaviors of
bicrystal with ’planar’ grain boundaries.
線で各変形量の線形分布の解析解を表示する. これより,系が線形弾性応答を示すことが確認できる. 更に, いずれの変形速度でも同様の結果が得られることから, ひぞみ速度に依存しないことが確認できる. 以上より, 本手法を用いて, 適切に弾性変形挙動を再現できると考えられる.

3·2

二結晶体のシミュレーション

図 3(a) と図4(a) に二結晶体の初期構造を自由エネルギーによって表示している. 自由エネルギーで表示することで, φ の周期性が乱れた格子欠陥領域ではエネルギーが高く表示されるため, 欠陥の位置を容易に把握できる利点を有する. 図に示すように粒界形状が直線・非直線形状である二種類の二結晶体を対象として変形解析を行う.
図 3(b)(c) に非直線形状の粒界を有する二結晶体の変形過程を示す. 変形が進むにつれて, 粒界から転位が放出される. これらの転位の内, 転位 A は粒内を移動して反対側の粒界に吸収されることが分かる. これより, 粒界が転(c) ε = 0.080.041

Fig. 6: Initial structure and deformated states ofnanopolycrystalline structure.
正六角形定型粒のナノ多結晶体を初期構造として用いる.
このような初期構造を形成させるため, 正六角形の中心に図示の結晶方位を有する核を配置し, 凝固シミュレーションを行う. 生成後の初期構造は図 6(a) に示され, 所定の多結晶構造が得られていることが確認できる.
図 6(b)(c) にナノ多結晶体の変形過程を示す. 変形が進行するにつれて, 粒回転が生じ, 粒界上の転位が移動する.
例えば, 転位 C は移動して粒界で吸収され, 粒が融合する(図 6(b)). これより粒界三重点が消失するため, 融合粒近傍の粒界は不安定な状態になる. 特に, 粒界 D は曲率を駆動力として移動し, 融合粒が粗大化していくことが分かる (図 6(c)). このようにして, 粒回転による粒の融合および粒界移動による粒成長などの粒間変形が支配的な変形機構となることが確認できる. これらの結果は MD シ
ミュレーション (4) でも確認されており, 本手法を用いても, ナノ多結晶体の変形挙動を適切に再現できていると考えられる.

4位の発生源および吸収源となることが確認できる. 更に,今後の展望

二つの転位 B は相互作用によって消滅する. これより, 典

今後, 本研究では, 異なる温度領域における異粒径の多型的な金属材料の塑性変形挙動である転位の生成・移動・結晶体を対象とする変形シミュレーションを行い, 本手法消滅を確認できる.の変形解析手法としての汎用性を示す. 更に, より実用的

図 4(b)(c) に直線形状の粒界を有する二結晶体の変形な変形解析手法として確立するため, 原子間ポテンシャル過程を示す. 図 3(b)(c) と同様の塑性変形挙動を確認できを評価可能とする三次元 PFC モデルを構築し, 実材料へる. 更に, 図 3(b) と図 4(b) を比較すると, 非直線形状のの導入を図る.

粒界構造を有する二結晶体の方が, 転位を放出する傾向にあることが確認でき, 粒界構造に依存する変形挙動を再現できていることが分かる. これらの結果は MD シミュ
レーション (3) でも確認されており, 本手法を用いて二結晶体の塑性変形挙動を適切に再現できると考えられる.
ナノ多結晶体のシミュレーション
図 5 に示3·3
すような各粒の一辺におよそ 14 個程度の原子が含まれる

参 考 文 献(1) K. R. Elder et al, Phys. Rev. E 70, 051605 (2004).
(2) 廣内・高木・冨田, 機論, A(2008), 掲載予定.
(3) L. Capolungo et al, J. Mech. Phys. Solids 55 2300
(2007).
(4) A.J. Haslam et al, Acta Mater. 52 1971 (2004).