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

京都工芸繊維大学 [院] 高木 知弘
神戸大学 [院]冨田 佳宏

Construction of Deformation Simulation Strategy of Material
by Phase Field Crystal Method
Tomoyuki HIROUCHI, Tomohiro TAKAKI and Yoshihiro TOMITA

1緒言

∆x

∆ x+dBefore deformation

After deformation

工業用材料の物理化学的な特性は, 材料の微視組織内に存在する欠陥に強く影響されることが知られている. そのため, 材料の更なる高機能化および新機能の発現のためには, 材料の微視組織に依存する巨視的Fig. 1な材料特性評価を可能とする数理モデルの構築が必要不可欠である. 最近, Elder らによって, 拡散時間0.3スケールで原子の挙動を解析可能とする Phase Field

0.2

The technique of tensile deformation in one
dimension.
a : The distance between atoms

0.00150.

n : Number of difference latticesが提案された. 一般的に, 塑性変形過程における金属材料の体積は変化しないが, これまでの PFC 法を用いた多結晶体の変形シミュレーション(1)(2)

では, 境界条件の設定が困難とされ, 体積一定条件下での塑性変形を再現できていない.

0.1

Phase field

Crystal(PFC) 法(1)

Dimensionless stress σn0-0.10.001

08101214161820220.0005
0.

-0.2

Analysis stress-0.3050

そこで, 本研究では, 新たな変形シミュレーション

10000.005

Distance

Strain(a)

(b)

0.01

手法を提案する. 次いで, それを用いて差分格子サイズの応力評価におよぼす影響を検討するため, 一次元Fig. 2変形シミュレーションを行う. 更に, 二次元問題へ拡(a) Phase ﬁeld proﬁle at the initial condition and (b) stress-strain curves.

する. 提案した手法は, 差分格子の大きさを制御する∂ϕ= Γ ∇2 [{α∆T + λ(q0 2 + ∇2 )2 }ϕ + uϕ3 ] (2)∂t
ここで, α∆T は駆動力, u, λ, q0 は平衡状態時の物質に方法で, 二次元問題では, 体積一定の塑性変形を再現依存するパラメータ, Γ は ϕ のモビリティーである.

張し, ナノ結晶体の弾性変形挙動および塑性変形挙動を評価し, 構築した手法の有用性および妥当性を検討

3できるところに特徴を有する.

2Phase Field Crystal 法
PFC 法は, 単結晶中の原子配列のような周期状態を表現するために, 秩序変数が周期的な値をとることが大きな特徴であり, 従来の Phase ﬁeld 法と大きく異一次元変形シミュレーション

本章では，提案した手法を用いて一次元変形シミュレーションを行い，差分格子サイズの応力評価におよぼす影響を検討する．

3·1

差分格子の制御

図 1 に変形解析におけるなる点である. 本手法にて用いる秩序変数 ϕ は, 時間差分格子の制御方法を模式的に示す. 各時間ステップ

平均化された原子の数密度と定義され, 保存場である.

∆t 毎に差分格子の x 軸方向の長さを d = ε̇∆x∆t だ

固相領域においては原子の周期構造を表現するため,け増加させて, 引張り変形を表現している. ここで, ε̇ϕ は周期的な値をとり, ϕ の最大点を原子の位置とする. 液相領域においては時間に依存するランダムな原は無次元ひずみ速度である.

子配置を時間的に平均化して一定値 ϕ0 としている.

分布を示す. 原子間距離を a とし, 16a の解析領域を

3·2解析モデル

図 2(a) に phase ﬁeld の初期

ここで，この秩序変数を用いた結晶の周期構造を表用いる. 差分格子間隔 ∆x を ∆x = a/n とし, 差分格

す最も基本的な自由エネルギー汎関数を式 (1) に, ϕ子数 n を変化させて解析精度を評価する. 無次元ひずを保存場とするので時間発展方程式は式 (2) に示す.
][∫uϕ{α∆T + λ(q0 2 + ∇2 )2 }ϕ + ϕ4 (1)F = dV24

み速度 ε̇ を ε̇ = 9.55 × 10−6 とし, 各種パラメータを,(α∆T, ϕ0 , u, λ, q0 , Γ ) = (−0.05, 0, 1, 1, 1, 1) とする.
図 2(b) に応力-ひずみ線図を示
3·3
解析結果す. 応力 σ は, 前時間ステップと現ステップのひず∆x

φ -0.61yS∆ x'x60a 600 lattices)
Periodic boundary condition in all direction∆ x+dBefore deformation

Fig. 3Tensile direction

3a (30 lattices)∆xS0.42After deformation

Fig. 4

The technique of tensile deformation in two

The initial conditions of the single crystal.

33Strain

みの変化量 ∆ε に対する自由エネルギーの変化量 ∆Fの割合 σ = ∆F/∆ε により算出する．図 2(b) の点線は, 解析解である. これより, n が大きくなるにつれDisplacement

dimensions.

0.002
0.004
0.006

2211て, 数値解が解析解に近づくことが確認できる. また,n = 8, 10, 14, 22 と変化させた時の誤差は 29, 15, 6, 101020

3040Fig. 5二次元変形シミュレーション

ミュレーションを行い, 弾性変形挙動および塑性変形

4050(b)Atomic displacement distributions at strain

A20deg.

40deg.

20deg.

B挙動を評価する.
差分格子の制御30

Periodic boundary conditions in all directions

本章では, 単結晶体および二結晶体の二次元変形シ4·120rate with (a) 5.17×10−5 and (b) 5.17×10−6 .

60a (600 lattices)410(a)

分格子を用いると, 二次元, 三次元問題では計算容量が膨大となるため, 次章では n = 10 を用いる.

0Position％程度であり, 1 ％程度の誤差で評価するためには，

n = 22 以上の差分格子数が必要であることが分かる.
しかしながら, 一つの原子を表現するために 22 個の差50Grain boundary

図 3 に差分格子の制御方法を模式的に示す. 各時間ステップ ∆t 毎に差分格子60a (600 lattices)(a)Free energy(b)0.021(c)0.041の x 軸方向の長さを d = ε̇∆x∆t だけ増加させ, 常に差分格子の面積が一定となるように y 軸方向の長さを∆x′ = (∆x)2 /(∆x + d) とし, 体積一定変形を表現することができる.
4·2単結晶体の変形シミュレーション
図4に ϕ の初期状態を示す.

変形前後の原子の位置を容易に比較できるように, 変形方向に対して平行にFig. 6(a) The initial conﬁgurations of grainboundaries and atoms in the bicrystal and,(b) and (c) the deformated states withstrain ε = 0.061 and 0.068 .

形成し, 前節と同様の無次元ひずみ速度 ε̇ = 5.17×10−6および各種パラメータを用いて解析を行う.

原子を配置し, 無次元ひずみ速度 ε̇ = 5.17 × 10−5

図 6(b), (c) に二結晶体の変形過程を示す. これよおよび 5.17 × 10−6 について解析を行う. 実際のひり, 粒界から転位が放出され, A に示す 2 つの転位がずみ速度では, およそ 101 /s〜102 /s に相当する. ま相互作用により消滅し, 更に, 転位 B が粒界に吸収さた, 各種パラメータを, (α∆T, ϕ0 , u, λ, q0 , Γ , ∆t) =れることが分かる. 粒界が転位の発生源および吸収源(−0.25, −0.285, 1, 1, 1, 1, 3.2 × 10−3 ) とする.となることは, MD シミュレーションでも報告されて

図 5 に各変形量における原子の変位分布を示す. また, 点線で各変形量の解析解を表示する. これより, いいる.(3) これより, 本手法を用いて, 二結晶体の塑性変形挙動を再現できていると考えられる.

ずれの変形速度領域においても, 変位の線形分布を確認できる. これより, 本手法が線形弾性応答を精度良く表現可能であることを確認できる.

4·3二結晶体の変形シミュレーション

図 6(a)に二結晶体の初期構造の自由エネルギー分布を示す. 自由エネルギーの表示により, 格子欠陥領域が高いエネルギーを示すため, 格子欠陥の位置を容易に把握できる.
各々の粒方位が θ = 20◦ , 40◦ となるように初期構造を

参 考 文 献(1) K. R. Elder and M. Grant, Phys. Rev. E 70,051605 (2004).
(2) T. Hirouchi et al, Proc. of Third APCOM (CDROM), MS31-6-3, 1-8., (2007).
(3) L. Capolungo et al, J. Mech. Phys. Solids 552300-2327 (2007).