132 PhaseFieldCrystal 変形シミュレーションによる温度・粒径依存性の特性評価
Evaluation of Temperature and Grain-size Dependences in Phase Field Crystal Deformation Simulation
 廣内 智之 · 神戸大院

高木 知弘 · 京工繊大院

冨田 佳宏 · 神戸大院

Tomoyuki HIROUCHI,
Graduate School of Engineering,
Kobe University

Tomohiro TAKAKI,
Graduate School of Science and Technology,
Kyoto Institute of Technology

Yoshihiro TOMITA,
Graduate School of Engineering,
Kobe University

Key Words : Phase Field Crystal Method, Multi Scale Model, Nanopolycrystalline Metal, Plastic Deformation

2 Phase Field Crystal 法
2.1 秩序変数と支配方程式
PFC 法は, 単結晶中の原子配列のような周期状態を表現するために, 秩序変数が周期的な値をとることを大きな特徴としており, 従来の Phase ﬁeld 法と大きく異なる点である.
本手法にて用いる秩序変数 φ は, 時間平均化された原子の数密度と定義され, 保存場である (1) . 固相領域においては原子の周期構造を表現するため, φ は周期的な値をとる. また, φの最大点を原子の位置とする. 液相領域においては原子のランダムな配置を時間的に平均化して一定値 φ0 としている.
ここで，この秩序変数を用いた無次元化された自由エネルギー汎関数 (1) を以下に示す．
α 2 φ4φ2 2
F = dx(1)φ ++ (1 + ∇ ) φ242
ここで, α は温度に関係する駆動力項で, 過冷度に比例する係数である. 秩序変数を保存場とするため, φ の時間発展方程式は次のように表すことができる.
∂φ(2)= ∇2 [αφ + φ3 + (1 + ∇2 )2 φ]∂t

2.2 体積一定変形の解析手法
変形解析における差分格子の制御法を模式的に図 1 に示す. 各時間ステップ ∆t 毎に差分格子の x 軸方向の長さをd(= ε̇∆x∆t) だけ増加させて ∆x とし, 常に差分格子の面積 S(= ∆x∆y = ∆x ∆y  ) が一定となるように y 軸方向の長さを決定する. ここで, ε̇ は無次元ひずみ速度, ∆x, ∆y はx, y 軸方向の初期差分格子間隔, ∆x , ∆y  は変形過程における x, y 軸方向の差分格子間隔である. 本手法は, 全ての原

Initial state ∆ x = ∆ yDeformed state ∆ x' = ∆ y's
Analysidomainy∆yx

Fig. 1S∆xS∆ x'∆ y'S = ∆ x∆ y
= ∆ x'∆ y'

Numerical scheme of tensile deformation simulation in two dimensions.
Tensile direction

75a (750 lattices)
Tensile direction

y58.1

35.2

41.9
Crystal orientation
4.69
0.32

35.2

21.2

18.6

21.2

x

Fig. 244.8

41.7

14atoms 42atoms

Crystal orientation

39.0

39.0
59.9

21.0

21.0

132a (1320 lattices)

17.4

15.9

19.1

58.1

1.64

57.4

1.64

37.9

41.7

44.8

Periodic condition in all directions

44a (440 lattices)

50a (500 lattices)
Periodic conditions in all direction

1 緒言
実材料の機能特性は, 材料内部の微視組織に存在する欠陥に強く影響を受ける. そのため, 材料の高機能化や新規素材の開発により材質向上を図るためには, 材料の微視組織に起因する巨視的な材料特性を評価できる数理モデルを構築することが急務とされている. 現在, 原子レベルの空間スケー
ルと, 拡散時間スケールにおいて物理現象を再現可能とする
Phase Field Crystal(PFC) 法 (1) が新たなマルチスケール手法として注目されている. 本手法は, 様々な物理現象に適用されているが, これまでの PFC 法を用いた多結晶体の変形シミュレーション (1)(2) では, 境界条件の設定が困難とされ,体積一定状態で進行する金属材料の塑性変形を再現できていない.
そこで, 本研究では, 体積一定の塑性変形を再現可能とする差分格子を制御する変形手法を提案する. 提案した手法を用いて, 温度と粒径を変化させたナノ多結晶体の変形シミュレーションを行い, 変形挙動に対する温度および粒径依存性を評価する.

37.9(a) Model A(b) Model B

Desired nanopolycrystalline structures with regular hexagonal grains.

子に一定のひずみ速度を加えるため, 変形状態はアフィン変形に相当する.

3 変形シミュレーション
本章では, 提案した手法を用いて, 温度と粒径を変化させたナノ多結晶体の変形シミュレーションを行い, 塑性変形挙動に対する温度・粒径依存性の影響を評価する. 初期構造は凝固シミュレーションにより生成する.
3.1 解析モデル
図 2(a)(b) に示す 2 種類の粒径を有する正六角形定型粒のナノ多結晶体モデル Model A と Model B を用いる. 図 2(a)の Model A では粒数を 12, 図 2(b) の Model B では, 計算負荷を軽減するため粒数を 6 とする. 各々の粒の一辺に, それぞれおよそ 14, 42 個の原子が含まれており, 粒径比は 1:3 である. また, 解析領域サイズはそれぞれ 44a×50a, 132a×75aであり, 440 × 500, 1320 × 750 の差分格子にて分割している.
ここで, a は最近接原子間距離である. 図 2 に示す初期構造を形成させるため, 正六角形の中心に図示の結晶方位を有する核を配置し, 凝固シミュレーションを行う. 高温状態を想定し (3) , 凝固過程を安定に解析するため, パラメータとして(α, φ0 ) = (−0.25, −0.285) を用いる. 生成後の初期構造は,
図 3(a) および図 4(a) に示され, 所定の多結晶構造が得られていることが確認できる. また, 格子欠陥の位置を容易に把〔№08-33〕日本機械学会第21回計算力学講演会ＣＤ-ＲＯＭ論文集〔2008-11.1〜3・沖縄県西原町〕

－420－

G

BA

(a)

(b)

Fig. 3(a)(c)

Free energy 0.021(b)Free energy -0.016

0.041

Deformation process for Model A of α = −0.25
at ε = (a) 0.00, (b) 0.05 and (c) 0.10.

Fig. 5(c)0.012

Deformation process for Model A of α = −0.50
at ε = (a) 0.00, (b) 0.05 and (c) 0.10.

FKLECJHI

D(a)
Free energy 0.021

Fig. 4(b)(a)(c)(b)

Free energy -0.016

0.041

Deformation process for Model B of α = −0.25
at ε = (a) 0.00, (b) 0.05 and (c) 0.10.

握できるようにするため, 自由エネルギーを用いて表示している. これらの初期構造に対して, α = −0.50 として緩和シミュレーションを行うことにより, 前者よりも低温領域を想定した多結晶構造を形成させる (図 5(a), 図 6(a)). また, 無次元ひずみ速度 ε̇ = 1.84 × 10−6 の下で x 軸方向への引張り変形シミュレーションを行う. これは, 実際のひずみ速度では, 100 s−1 〜101 s−1 に相当する. また, 無次元時間ステップを ∆t = 3.0 × 10−3 とし, 境界条件は周期境界条件を適用する.

3.2 解析結果
図 3 に, Model A(α = −0.25) の変形組織を示す. 変形開始から粒が回転し, 図 3(a)(b) に示す矢印の方向に転位が移動することが分かる. 特に転位 A が粒界で吸収され, 粒が融合している (図 3(b)). このため, 粒界上の三重点は消失し,融合粒周辺の粒界は不安定な状態になる. 更に, 図 3(b) の楕円 B で示す粒界は, 曲率を駆動力として移動し, 融合粒の粗大化過程が確認できる. このようにして, 粒回転および粒界移動による粒成長が支配的な変形メカニズムとなることが確認できる.
図 4 に, Model B(α = −0.25) の変形過程を示す. 図 3 と同様, 粒界上での転位の移動が観察できるが, 粒径が大きいため粒回転は生じ難く, 変形量が大きいときも粒の融合は確認できない. 一方, 粒界の形状が変化することが分かる. 例
えば, 図 4(b) の楕円 C で示すように, 粒界が曲率を有する形状に変化することが確認できる. これは, 粒界 C の粒界エネルギーよりも応力による弾性エネルギーが大きくなるために生じると考えられる. 更に, 粒界 C のエネルギーの均衡を保つため, 図 4(b) の粒界 D が矢印の方向に移動し, 粒界C が直線形状になることが確認できる (図 4(b)(c)). 対照的に, 図 4(b) の楕円 E, F で示す粒界では形状が大きく変化しFig. 6(c)0.012Deformation process for Model B of α = −0.50
at ε = (a) 0.00, (b) 0.05 and (c) 0.10.

ないことが確認できる. これは, 粒界 E, F が対称性の強い結晶粒から構成され, 更に, 粒界と荷重の方向が一致するため, 粒界における構造変化が生じ難くなると考えられる.
図 5 に, Model A(α = −0.50) の変形挙動を示す. 例えば,
図 3 の転位 A と図 5 の転位 G を比較すると, 転位 G の移動量が小さいことが分かる. これは, 温度の低下によって転位の移動に対する抵抗力が高くなることと一致している. また, 主に粒回転や粒界移動が生じて構造変化に対応しているが, 粒の融合および粗大化は確認できないことから, 図 3 に比べてこれらの変形機構が生じ難くなることが分かる. 一方, 図 5(b) に示す粒界上の転位が, 粒内を矢印の方向に移動して粒界で吸収される. このような転位の移動は, 図 3 では確認できないことから, 構造変化に対応する粒回転や粒界移動などの変形機構を補填する役割を果たすと考えられる.
図 6 に, Model B(α = −0.50) の変形過程を示す. 図 4 と同様で, 粒径が大きいため粒回転が生じ難いことが分かる.
一方, 図 6(a) の粒界 H, I のように移動, もしくは, 粒界 J,K, L のように形状が変わることが確認できる. これは, 粒界を構成する結晶粒の対称性の強さに起因すると考えられる.
更に変形が進むと, 図 6(b) の粒界 H, I 上の転位は, 図 5 で観察されるように移動して粒界で吸収される. 更に, 図 6(b)の粒界 K, L 上からも転位が移動することが確認できる. これは, 粒界形状が変化して非直線性が強くなるため, 応力集中を受けて転位が移動すると考えられる. また, Model A に比べて粒界長が長いため, 粒界形状が容易に変化し, 多数の転位が粒内に移動すると考えられる (図 6(c)).

－421－

参考文献
(1) K. R. Elder et al, Phys. Rev. E 70, 051605 (2004).
(2) T. Hirouchi et al, Proc. of Third APCOM (CD-ROM), MS316-3, 1-8., (2007). MS31-6-3, 1-8., (2007).
(3) J. Berry et al, Phys. Rev. E 73, 031609 (2006).