修士論文

ナノ多結晶体の内部不均一性に関する
局所格子不安定性解析

指導教員：冨田 佳宏

067T395N 樋口 昌宏

2008 年 2 月
神戸大学大学院 自然科学研究科 博士課程前期課程 機械工学専攻

Master’s Thesis

Local Lattice Instability Analysis on
Inhomogeneity in Nano-Polycrystals

Masahiro HIGUCHI
February 2008
Department of Mechanical Engineering,
Graduate School of Science and Technology,
Kobe University, Kobe, Japan要約
α
内部不均一を有する系の変形挙動を，原子弾性剛性係数 BIJの正値性の観点から議論する研究の一環として，2 次元状解析モデルで結晶粒形状および寸法の異なるナノ多結晶体を対象とした種々の分子動力学シミュレーションを行った．まず，無負荷平衡状
α態における系の det BIJ および局所の det BIJの平均が，結晶粒形状や寸法でどのように変化するかを調べた．次に，結晶寸法を変えたナノ多結晶体の単軸引張シミュレー
ションを行い，0.2% 耐力には粒形状および寸法の影響が少ないこと，ボロノイ分割した多結晶体では，最大応力が粒径 10nm 近傍でピークとなること，等を明らかにした．
最後に，引張変形下における不安定原子の物理的描像を明らかにするため，結晶寸法α約 10nm のナノ多結晶体について，det BIJの正負および局所構造の対称性から同定した fcc，hcp，その他の欠陥 (other) で原子を分類して，その割合および応力変化を詳細αに検討した．その結果，(1) 無負荷における det BIJ< 0 の不安定原子の多くは other である，(2)fcc，other の不安定原子は無負荷で応力 0 であるのに対し，粒界部の小数の不安定 hcp 原子が高い静水圧引張応力を示す，(3)0.2% 耐力ひずみを超えると，粒内に積層欠陥が発生するため fcc → hcp の間で割合の変化を生じるが，other はほとんど変わらない，(4)0.2% 耐力ひずみを超えると other の不安定原子のみが増加する．これは構造上「非結晶」な原子の割合は変わらないまま，それらの原子が不安定となっていくことを意味する，(5)(2) で述べた不安定 hcp 原子の応力は常に他より高いが，0.2%
耐力ひずみを超えると急減した．これは 3 次元のナノ多結晶体ではみられず，また 2 次
元でも引張方向に垂直な粒界のないナノ多結晶体ではみられない，等を明らかにした．

Summary
As a series of studies which discuss the deformation behavior of inhomogeneous systems
αbased on the positive deﬁniteness of atomic elastic stiﬀness coeﬃcients, BIJ, we haveconducted several molecular dynamics simulations on nano-polycrystals with diﬀerent
grain shapes and sizes, as the two-dimensional model. First, we have discussed theαchange in det BIJ for the whole the crystal and the average of det BIJ
by grain shapeand size before loading. Then, the several nano-polycrystals with diﬀerent grain sizes
are subjected to uniaxal tension, revealing that 0.2% proof stress is not inﬂuenced by
the grain shape and size, and the polycrystals with Voronoi grains show the highest
maximum stress at the grain size of 10nm. The last, in order to clarify the physicalmeaning of unstable atoms under tensile deformation, we have classiﬁed all the atoms
αby the positive/negative of det BIJ, and local crystallographic structures, i.e. fcc,

hcp, and the other defect judged by CNA, in the simulations of polycrystals withαgrain size of 10nm. It reveals the followings; (1) many of the negative det BIJ
atoms

at εyy = 0 are neither fcc nor hcp, but the other disordered defects, (2) the initialresidual stresses on the unstable fcc and unstable other defects are zero at εyy = 0. On
the other hand, the fraction of the unstable hcp atoms is very small but shows very
high positive hydrostatic stress, (3) beyond the strain of 0.2% proof stress, stacking
faults are emitted from grain boundaries leading to the change in the ratio between
fcc and hcp atoms. The ratio of the other defects stays almost constant throughout
tension, (4) the number of unstable atoms increase only in the other defects, after
the strain for 0.2% proof stress. This means that the ratio of the other defects does
not change, but the stable → unstable transition takes place in the defects, (5) the
unstable hcp atoms always show highest stress during tension, while the stress shows
rapid drop at the strain of 0.2% proof stress. This stress drop is not observed neither inthe three-dimensional nano-polycrystals nor in the two-dimensional polycrystals whichhave no grain boundary normal to the tensile direction.

目次
第1章 緒 論1

第 2 章 解析手法

5

2.1 分子動力学法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.2 原子埋め込み法ポテンシャル. . . . . . . . . . . . . . . . . . . . . . .

6

2.3 速度スケーリング法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.4 領域分割による高速化手法 . . . . . . . . . . . . . . . . . . . . . . . . .

92.5 動径分布関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.6 ボロノイ分割 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.7 局所格子不安定性解析 . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.7.1

弾性剛性係数と格子不安定性 . . . . . . . . . . . . . . . . . . .

13

2.7.2

原子応力ならびに原子弾性係数 . . . . . . . . . . . . . . . . . .

13

第 3 章 無負荷平衡状態における局所格子不安定性

17

3.1 シミュレーション方法 . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.1

解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.2

結晶寸法の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.1.3

初期緩和シミュレーション . . . . . . . . . . . . . . . . . . . . .

20

3.1.4

アモルファス構造の作成ならびに構造評価 . . . . . . . . . . . .

20

3.2 シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

21

3.2.1

不安定原子の分布 . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.2.2

不安定原子に生じる応力 . . . . . . . . . . . . . . . . . . . . . .

24

3.2.3

粒界および結晶内部の不安定原子 . . . . . . . . . . . . . . . . .

25

3.3 結 言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

i

目次
第 4 章 引張変形下の力学応答ならびに内部構造変化

ii
29

4.1 シミュレーション方法 . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

4.2 シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

31

4.2.1

応力-ひずみ関係

4.2.2

0.2% 耐力および最大応力の結晶寸法依存性. . . . . . . . . . .

32

4.2.3

引張変形時の原子構造変化 . . . . . . . . . . . . . . . . . . . . .

33

4.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35. . . . . . . . . . . . . . . . . . . . . . . . . .

第 5 章 引張変形下の力学応答と不安定原子の力学状態

3136

5.1 解析モデルおよびシミュレーション方法 . . . . . . . . . . . . . . . . .

36

5.2 シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

38

5.2.1

応力-ひずみ関係. . . . . . . . . . . . . . . . . . . . . . . . . .

38

5.2.2

不安定原子に生じる応力 . . . . . . . . . . . . . . . . . . . . . .

39

5.2.3

結晶構造毎の応力変化 . . . . . . . . . . . . . . . . . . . . . . .

40

5.2.4

結晶構造と不安定原子 . . . . . . . . . . . . . . . . . . . . . . .

41

5.2.5

引張軸に垂直な粒界原子の力学状態 . . . . . . . . . . . . . . . .

42

5.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

第6章 結 論

46

参考文献

48

関連講演論文

52

謝辞

67

第1章
緒論
物質を原子レベルで制御することで，革新的な材料・デバイスの開発を目指すナノテ
クノロジーが注目されている (1)−(3) ．ナノテクノロジーの対象範囲は IT，バイオ，医療など各種産業の広範な分野の基盤に関わるものであり，今後の科学技術の発展に大きく寄与することからその重要性は非常に高い．材料分野においても，新規機能材料の創製に向けナノテクノロジーに大きな期待が寄せられている．炭素原子単体でも多様な物質 (グラファイト，ダイアモンド，フラーレン，カーボンナノチューブ) が存在するように，原子レベルでの構造変化は材料特性に大きな変化をもたらす．材料の微細組織を積極的に制御することで，材料が持つ固有の特性を極限まで引き出そうという試みがなされている．金属材料におけるナノテクノロジーとしては，結晶粒をナノメートルレベルまで微細化したナノ多結晶体が挙げられる．
ナノ多結晶体は 1980 年代の中頃に Gleiter らが 10 ナノメートル以下の鉄の超微粉末から凝縮して作成したのが始まりとされる (4) ．Gleiter らは X 線散乱や Mössbauer 効果，磁気測定等を通して，従来の粗大結晶粒からなる多結晶体とは異なる新しい固体構造を持つことを報告している．実際，ナノ多結晶体は特異な力学特性を示すことが知られている．従来の粗大粒からなる多結晶体では，結晶粒を微細化すると転位運動の障害となる結晶粒界の割合が増加し，低温において降伏応力や変形応力が増加する
硬化現象を生じる．この硬化現象は Hall-Petch の関係として古くから経験的に知られている．しかしながら，ナノサイズまで微細化した多結晶体では，図 1.1 に模式的に示すように，逆に軟化する逆 Hall-Petch の関係が報告されている (5),(6) ．Chokshi らは1

2

Nano Crystals
Submicron Crystals

Macro Crystals

Inverse Hall-Petch
Relation

Amorphous

Yield stressHall-Petch Relation
σ=σ0+kd -1/2

〜 10-20nm〜 1µ mGrain size
Fig.1.1 Change in yield stress by grain size.

ナノ結晶構造を有する Cu および Pd の試験片に対して硬さ試験を行い，6 ∼ 16nm の粒径範囲で微細化に伴う軟化現象を確認した (5) ．また，McFadden らは Ni，Al 合金(1420-Al) および Ni3 Al のナノ多結晶体について超塑性の検討を行い，Ni はこれまで
よりも 470K 低い 0.36Tm (ここで，Tm は融点) で，また Ni3 Al についてはミクロ結晶性領域よりも 450K も低い温度で超塑性を示すことを報告している (7) ．Al 合金についても 1 × 10−1 1/s という高ひずみ速度条件下で超塑性が発現する．辻らは ARB 法による強ひずみ加工を用いて作成したサブミクロンサイズの種々の Al 多結晶体について単軸引張試験を行い，耐力およびひずみ硬化指数に Hall-Petch の関係を確認したが，粒径1µm 以下では延性の低下が生じることを報告している (8) ．これらナノ多結晶体特有の力学特性の発現は，粒界割合の増大による粒間変形の促進によって，従来の粗大粒多結晶体とは異なる独自の変形メカニズムを持つためと考えられている．
試料作製および実験観察の困難さを避け，またモデルを単純化することで支配的なメカニズムを切り分けるには，分子動力学法に代表される計算機シミュレーションが有効である．Schiøtz らは，分子動力学法により Cu ナノ多結晶体の単軸引張解析を行

3い，粒界領域の体積占有率の増加がナノ多結晶体を軟化させることを報告している (9) ．
中谷らは Al ナノ多結晶体について，積層欠陥の生成を伴う部分転位によるすべりが主たる変形を担うことを確認し，粒の回転に伴う変形集合組織の発現について報告している (10) ．Van Swygenhoven らは，大角粒界と小角粒界をそれぞれ多く含む 2 種類の

Ni について，粒径 3∼12nm の範囲で，密度，粒界密度および過剰エンタルピーを比較し，個々の粒界特性が塑性変形下の挙動に大きく影響する可能性を示している (11) ．下川らは，六角形粒を有する粒径 5∼80nm の Al ナノ多結晶体について単軸引張解析を行い，粒内変形および粒間変形に対する変形抵抗の関係から，最高強度を実現する最適粒径は 30nm 近傍であることを報告している (12) ．
著者らのグループでは，局所格子安定性 (13) という独自の視点から，ナノ多結晶体
やアモルファス等の内部不均一を有する系の力学応答を議論してきた (13)−(21) ．局所格
子安定性は，各原子位置におけるエネルギーの空間勾配に相当する物理量 (原子弾性剛αα
性係数 BIJ) の正値性を議論するもので，det BIJ< 0 の正負で単結晶表面からの転位の発生 (13),(16) や微視的へき開の発生 (14) ，ナノインデンテーションにおける圧子下の転位発生 (22) などが評価できることが示されている．3 次元立方体セルをボロノイ分割
または格子状に分割したナノ多結晶体と，アモルファスを対象とした詳細な検討では(18)−(21)α，(1) 無負荷平衡状態において det BIJ< 0 の「不安定な」原子は，ナノ多結晶体では静水圧引張状態に，アモルファスでは静水圧圧縮状態にあること，(2) ナノ多結α晶体およびアモルファスいずれにおいても，det BIJ< 0 の原子は周囲よりも大きく変形し高い応力を示す「弱い」部分であること，等を明らかにしている．
本研究では，様々な結晶粒形状および寸法の Ni ナノ多結晶体について，薄板状セルに周期境界条件を適用した 2 次元状モデルで検討する．3 次元から 2 次元への移行は，本質から逆行するように思われるかもしれないが，3 次元では観察困難な内部の現象が2 次元にすることで観察が容易になること，Hall-Petch から逆 Hall-Petch への遷移粒
径 (10∼20nm) が扱えること，等の利点がある．また 3 次元から 2 次元への単純化により，3 次元との本質的違いが明確になるものと期待される．以下，各章の概略を示す．
第 2 章では解析手法について概説する．分子動力学法の基礎方程式を示し，原子間相互作用の評価に用いた原子埋め込み法ポテンシャルの概要について述べる．そして，

4
局所不安定性の評価を行うための局所格子不安定性解析について概説する．第 3 章では，種々の結晶寸法のナノ多結晶体およびアモルファスの作成と初期緩和シュミレー
ションを行った後，無負荷平衡状態における系ならびに局所の安定性を評価し，結晶
寸法ならびに形状による変化を検討する．第 4 章では，寸法の違うナノ多結晶体の単軸引張シミュレーションを行い，0.2% 耐力および最大応力の粒寸法依存性について検α
討する．第 5 章では，det BIJ< 0 の基準だけでなく局所構造の対称性から結晶，非結晶そして積層欠陥等に分類して，引張変形下のナノ多結晶体中の不安定原子の物理的描像を明らかにする．最後に，第 6 章で本研究の総括を述べる．

第2章
解析手法
本章では，分子動力学法の基礎方程式を示し，本研究で原子間相互作用の評価に用いた原子埋め込み法ポテンシャルの概要について説明する．続いて分子動力学計算における温度制御のために用いた速度スケーリング法，領域分割による高速化手法について述べる．さらに，アモルファス金属の原子構造を解析するために用いた動径分布
関数，また多結晶体作成の際に用いたボロノイ分割について説明する．最後に，局所の安定性を評価するために用いた局所格子不安定性解析について概説する．

2.1

分子動力学法

分子動力学法 (Molecular Dynamics Method;MD)(23),(24) は，系を構成する個々の原子についてニュートンの運動方程式mαd2 r α
= Fα
dt2

(2.1)

を作成し，これを数値積分することによって全原子の運動を追跡する手法である．ここで t は時間，r α ，mα はそれぞれ原子 α の位置ベクトル及び質量である．原子 α に作用する力 F α は，系全体のポテンシャルエネルギー Etot の空間座標についての勾配ベクトルとして次式のように求められる．

Fα = −

∂Etot
∂r α

(2.2)

式 (2.1) の数値積分には，Verlet の方法 (23) が簡便で高精度が得られるため MD 法で5

2.2 原子埋め込み法ポテンシャル6はよく用いられる．時刻 t + ∆t と t − ∆t での原子 α の位置ベクトル r α (t ± ∆t) を

Taylor 展開すると，()
dr α (t) (∆t)2 d2 r α (t)
3
+
+
O
(∆t)
dt
2
dt2
2 2 α
α
(
)
dr (t) (∆t) d r (t)
3
+
O
(∆t)
r α (t − ∆t) = r α (t) − ∆t
+
dt
2
dt2

r α (t + ∆t) = r α (t) + ∆t

(2.3)
(2.4)

となる．ここで，v α を時刻 t における原子 α の速度とすると，

dr α= v α (t)
dt

(2.5)

であり，式 (2.1) と式 (2.5) を，式 (2.3) および式 (2.4) に代入すると，()(∆t)2 F α (t)
3
+
O
(∆t)
2
mα
2
(
)
(∆t) F α (t)3r α (t − ∆t) = r α (t) − ∆tv α (t) +
+
O
(∆t)
2
mα

r α (t + ∆t) = r α (t) + ∆tv α (t) +

(2.6)
(2.7)

となる．両式の和と差をとると，

r α (t + ∆t) + r α (t − ∆t) = 2r α (t) + (∆t)2()
F α (t)
4
+
O
(∆t)
mα

(

r α (t + ∆t) − r α (t − ∆t) = 2∆tv α (t) + O (∆t)3

)

(2.8)
(2.9)

が得られる．これより，時刻 t + ∆t での位置ベクトルと t での速度は

r α (t + ∆t) = 2r α (t) − r α (t − ∆t) + (∆t)2v α (t) =()
F α (t)
4
+
O
(∆t)
mα

(
)
1
{r α (t + ∆t) − r α (t − ∆t)} + O (∆t)22∆t(2.10)(2.11)と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．初期の計算 (t = 0) では t = ∆t での座標 r α (∆t) は式 (2.6) と初速度から求められる．r α (∆t) と r α (0) が既知であれば，式 (2.10) を繰り返し適用することにより各粒子の座標を求められる．

2.2

原子埋め込み法ポテンシャル

式 (2.2) で示したように，原子 α に作用する力 F α は系のエネルギー Etot をポテン
シャルとして決定される．したがって，系のポテンシャルエネルギー Etot をいかに精2.2 原子埋め込み法ポテンシャル

7

度よく評価するかが重要となる．量子力学に基づき，電子や原子核のハミルトニアン
から系のポテンシャルエネルギーを精密に求めて原子の運動を追跡する第一原理分子
動力学法も試みられているが，計算量が極めて膨大になるため，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する原子間ポテンシャルが通常用いられる．
原子埋め込み法 (Embedded Atom Method ; EAM)(25),(26) は，Daw，Baskes らによって提案された原子間ポテンシャルであり，金属中の多体効果を良好に再現することか
ら広く用いられている．EAM では密度汎関数理論に基づき，まず金属材料における系のポテンシャルエネルギー Etot は原子を価電子雲中に埋め込むエネルギーと原子間の2 体間相互作用の和で与えられるとする．さらに，埋め込みエネルギーは埋め込む位置の電子密度にのみ依存すると仮定することによって，系全体のエネルギーは次式のように表す．
N ∑N
1∑ϕαβ (rαβ )Etot =Fα (ρ̄α ) +2 α β(̸=α)α
N∑(2.12)

ここで，ρ̄α は原子 α の位置における電子密度，Fα (ρ̄α ) は電子密度 ρ̄α の位置に原子を埋め込むエネルギー，ϕαβ (rαβ ) は距離 rαβ 離れた原子 α と β のクーロン相互作用である．また，ρ̄α は周囲の原子 β からの寄与 ρβ (rαβ ) の重ね合わせで与えられると仮定しρ̄α =neighbor∑ρβ (rαβ )(2.13)

β(̸=α)

で評価する．

Voter ら (27),(28) は Ni 単結晶に対して昇華エネルギー，空孔形成エネルギー，弾性定数，格子定数等へのフィッティングを行い，ρ(r)，ϕ(r) について以下の関数形を提案している．

ρ(r) = S r6 (e−β r + 29 e−2 β r )(2.14)

ϕ(r) = D {1 − exp[−α (r − R)]}2 − D − 2 g ρ(r)(2.15)

式中のパラメーターの値は表 2.1，2.2 に示す．Fα (ρ̄α ) は原論文 (29) と同様に Rose らの凝集エネルギー関数 (30) を用いて数値的に求め， 3 次のスプライン関数によりフィッ2.3 速度スケーリング法

8

ティングした．フィッティング範囲は 0.0 ≤ ρ̄ ≤ 1.0 で，スプラインノードの間隔は∆ρ̄ = 0.01 とした．
Table 2.1 Potential parameters for ρ (r).
β (Å−1 )3.6408Ni

S
1.0000

Table 2.2 Potential parameters for ϕ (r).
Ni-Ni

2.3

α (Å−1 )1.7728D (eV)
1.5335R (Å)
2.2053g (eV Å3 )6.5145

速度スケーリング法

分子動力学法における温度制御には，もっとも簡単で直接的な方法として速度スケー
リング法 (23),(24) がよく用いられる．熱統計力学より系の運動エネルギー K は次のよ
うに系の温度と関連づけられる．

K=

N
1∑
3
mα v α ·v α = N kB T
2 α=12

(2.16)

ここで，mα は原子 α の質量，v α は原子 α の速度，N は系の全原子数，kB はボルツマ
ン定数，T は系の温度である．式 (2.16) より，系の温度 T は原子速度を用いて，次のように求められる．
∑T =mα v α ·v α
3N kB

(2.17)

設定温度が TC ，式 (2.17) より求めたある時刻の温度が T のとき，速度スケーリング
法では，各原子の速度 v α を√

TC /T 倍し設定温度 TC に近づける．ベルレ法では，∆r α (t + ∆t) = r α (t + ∆t) − r α (t) = r α (t) − r α (t − ∆t) + (∆t)2√をF α (t)(2.18)
m

TC /T ∆r α (t + ∆t) で置き換えることに相当する．平衡状態では，能勢の方法 (23)など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている．

2.4 領域分割による高速化手法

2.49

領域分割による高速化手法

式 (2.12) からわかるように，N 個の原子からなる系では，Etot の評価に N×(N−1) 回の原子対の計算が必要となる．一方，実際の結晶中では近接原子による遮蔽 (screening)
効果により第二近接距離程度より離れた原子はほとんど作用を及ぼさないことが知られている．このため，分子動力学計算では相互作用打ち切り (カットオフ) 半径 rc を導入し (図 2.1)，その半径内の原子からの寄与のみを考慮する．しかしながら，相互作用する原子対の検索に N × (N − 1) 回の試行を要するため，系が大きくなるにつれ計算負荷が飛躍的に増加する．これを避けるために rc よりひとまわり大きい半径 rfc (図 2.1)
内の原子をメモリーに記憶し，rfc 内での原子対の探索とすることによりオーダー N の計算に近づける方法 (粒子登録法 (24) ) がこれまでよく用いられてきた．しかしながら，粒子登録法では rfc 半径より外の原子が r c 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．このため，系がある程度の規模以上に大きくなると，粒子登録による高速化は
登録更新の負荷により打ち消される．

rfc

rc

Fig.2.1 Schematic of bookkeeping method.

2.4 領域分割による高速化手法

10

そこで本解析では計算負荷の増大なく効率よく rfc 内の原子対を探索可能な領域分割法を用いる．領域分割法では，まず図 2.2 に模式的に示すようにシミュレートする系
をカットオフ距離程度の格子状に分割する．ある原子に作用する力を評価する際には，その原子が属する領域（図 2.2 の着色部）と隣接領域内 (図 2.2 の斜線部) の原子からカットオフ距離内の原子を探索する．原子が属する領域は，位置座標を領域ブロックの辺長 bx , by で除した際の整数により判断できるので，領域分割そのものの計算負荷は小さい．領域分割法は，粒子登録法において登録更新の負荷が大きくなるような大規模な系の高速化に適している．

y
by
0

x
bx

Fig.2.2 Schematic of domain decomposition.

2.5 動径分布関数

2.5

11

動径分布関数結晶のように規則的な格子配置を取らず，均質で一様とみなせる液体やアモルファスのような構造では，どの点をとってみても，原子の平均密度は同じである．つまり，一体の分布関数は場所によらず数密度 N/V (≡ ρ) に等しい一定値をとる (24) ．
液体やアモルファス中の原子配置の特徴は，2 体相関関数 n(2) (r, r ′ ) や動径分布関
数 g(r, r ′ ) に現れる．n(2) (r, r ′ ) は 2 個の粒子を体積要素 dr と dr ′ 内に見出す確率がn(2) (r, r ′ )drdr ′ に比例するように定義される．これに対し動径分布関数は 1 個の粒子が r に存在するとき，位置 r ′ の体積要素 dr ′ = (dx′ , dy ′ , dz ′ ) 内に存在する平均粒子数を ρg(r, r ′ )dr ′ に等しいとおくことによって定義される．dr ′ 内の平均粒子数 ρdr ′ は体
積要素 dr ′ の大きさが同じであれば，r ′ をどこにとっても変わらない．ところが，位置 r に他の粒子が存在したとき，その影響が r ′ に存在する粒子にどう現れるかを表わすのが g(r, r ′ ) である．n(2) (r ′ , r ′′ ) を体積 V について r ′ と r ′′ で積分した値が∫ ∫n(2) (r ′ , r ′′ )dr ′ dr ′′ =N!
(N − 2)!

(2.19)となるように規格化すると，g(r ′ , r ′′ ) との間には N/V を一定に保って N → ∞ としたときn(2) (r ′ , r ′′ ) = ρ2 g(r, r ′ )(2.20)の関係があることが示される (31) . 均質等方で一様とみなせる場合，相対ベクトル r ′ −r ′′をどのように平行移動したり回転させても n(2) (r ′ , r ′′ ) は変わらないから |r ′ − r ′′ | ≡ rのみの関数となる，つまりg(r ′ , r ′′ ) = g(r)(2.21)

lim g(r) = 0(2.22)

lim g(r) = 1(2.23)

が成立する．またr→0r→∞

2.6 ボロノイ分割

12の関係がある．(2.22) 式はすでに粒子が存在する位置に他の粒子が入り込めなくなること，(2.23) 式は 2 粒子が遠く離れると，互いに存在する影響がなくなってしまうことを示している．

MD 計算により得られた粒子の座標のデータから動径分布関数を求めるには以下のようにする．nk (r, t) を時刻 t に粒子 k を中心とした半径 r − ∆r/2 と r + ∆r/2 の 2 球
面ではさまれた球殻中の粒子数とする．nk (r, t) の平均 ⟨n (r)⟩ を求めると (24) ，⟨n (r)⟩ = ρg(r)4π 2 ∆r

(2.24)であるので，これより g(r) が得られる．

2.6

ボロノイ分割本研究では多結晶モデルの作成にボロノイ分割を用いる．ボロノイ分割とは平面上に配置された無数の点群において母点と呼ばれる幾つかの点を選び，点群内の点をどの母点に最も近いかによって分割する手法である．分割される領域のことをボロノイ領域と呼び，ボロノイ領域の境界線をボロノイ境界, またはボロノイ辺と呼び，その交点をボロノイ点と呼ぶ．ボロノイ分割の一例を図 2.3 に示す．

Fig.2.3 Voronoi division.

2.7 局所格子不安定性解析

2.7

局所格子不安定性解析

2.7.1

弾性剛性係数と格子不安定性13

応力 σij および弾性係数 Cijkl は，等温過程では

1
σij =
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

∂ 2F
∂ηij ∂ηkl)(2.25)と定義される (32) ．ここで，F は Helmholtz の自由エネルギー (断熱過程では内部エネ
ルギー U )，V は結晶の体積，ηij は平衡状態 (無負荷とは限らない) からの仮想的な微小ひずみである．一方，無負荷平衡状態を基準とするひずみ εij と応力 σij の関係は，2つの平衡状態間の変形を考えて導出される次の弾性剛性係数によって表される (32) ．
(

Bijkl ≡

∂σij
∂εkl

)

= Cijkl + (σil δjk + σjl δik+ σik δjl + σjk δil − 2σij δkl )/2(2.26)

ここで，δij はクロネッカーのデルタである．すなわち，Bijkl は非線形弾性における応力-ひずみ関係の勾配を表す．Wang, Yip らは，ひずみの対称性を考慮したテンソ
sym
ル Bijkl≡ (Bijkl + Blkji )/2 の正値性によって格子不安定性を評価することを提案してsymsym
いる (33) ．det Bijkl≥ 0 であれば安定，det Bijkl< 0 であれば不安定である．すなわ
sym
ち，変形抵抗の喪失点が不安定の物理的描像である．なお，Bijkl
は Voigt 対称性 (32)symを有するので，Bijkl
を Voigt 表記した 6 × 6 行列 BIJ を用いれば上述の不安定条件はdet BIJ < 0 と表される．

2.7.2

原子応力ならびに原子弾性係数αの算出に必要な微視的応力な局所の安定性を評価するための原子弾性剛性係数 BIJ

らびに弾性係数は，各原子周りの微小ひずみに対するポテンシャルエネルギーの 1 次，

2 次変化量として導出される．

2.7 局所格子不安定性解析

14

応力
簡単のため，結晶の内部エネルギー U が系全体のポテンシャルエネルギー Etot に等しいとする．このとき，応力は平衡状態からの微小ひずみ η に対するポテンシャル
エネルギーの単位体積当たりの変化として与えられる (32) ．

σij =1 ∂EtotV ∂ηij(2.27)

ここで，V は平衡状態における系の体積であり，下付添字のローマ文字はテンソルのデカルト座標成分を表す．(2.27) 式の微分を求めるため，平衡状態からの仮想的な均一変形を考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によってr α = J r̄ α(2.28)と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子 αと 原子 β の間の距離 rαβ には(rαβ )2 = r̄iαβ Gij r̄jαβ(2.29)なる関係が成立する．ただし，Gij = Jki Jkj である．仮想変形の Lagrange ひずみテン
ソル ηij は

ηij =

]
1[
Gij − δij
2(2.30)であり，その微小量

1dηij = dGij2(2.31)と式 (2.29) の関係から次の関係が得られる．

r̄αβ r̄αβ∂rαβ= i αβj∂ηijr(2.32)

これより EAM ポテンシャルにおける応力は次式で評価される
(

σij =

N ∑
N
1 ∂Etot
1 1∑∂rαβ ∂Etot=
V ∂ηijV 2 α β(̸=α) ∂ηij ∂rαβ)

} αβ αβ
N ∑
N {
r r
1 ∑
′ α
′ β
′ αβ
′ αβ
=
[F (ρ ) + F (ρ )]ρ̃ (r ) + ϕ (r ) i αβj2V α β(̸=α)
r} αβ αβN ∑N {1 ∑1 ′ αβ ri rj′ α ′ αβ=F (ρ )ρ̃ (r ) + ϕ (r )V α β(̸=α)
2
rαβ(2.33)

2.7 局所格子不安定性解析

15

ここで，各原子位置における原子応力を

σijα =} αβ αβN {1 ′ αβ ri rj
1 ∑′ α ′ αβF (ρ )ρ̃ (r ) + ϕ (r )V /N β(̸=α)
2
rαβ(2.34)と定義すると，系の応力は

σij =

N
1 ∑
σα
N α ij(2.35)となる．
弾性係数
弾性係数も応力と同様に U ≈ Etot の場合にはCijkl =

1 ∂ 2 Etot
V ∂ηij ∂ηkl

(2.36)であるので，平衡状態からの仮想均一変形を考えると EAM ポテンシャルにおける弾性係数は以下のようになる．
(N ∑} αγ αγ )N {N ∑N
1 ′ αγ ri rj
1 ∑∂rαβ ∂ ∑′ α ′ αγF (ρ )ρ̃ (r ) + ϕ (r )Cijkl =
2V α β(̸=α) ∂ηkl ∂rαβ α γ(̸=α)
2
rαγ
[N ∑
{
} αβ αβ αβ αβN1 ∑ρ̃′ (rαβ ) ri rj rk rl
′ α
′′ αβ
=
F (ρ ) ρ̃ (r ) −V α β(̸=α)
rαβ(rαβ )2

+

N
∑

}
N
rαβ rαβ }{ ∑rαγ rαγρ̃′ (rαγ ) k αγlρ̃′ (rαβ ) i αβjrrγ(̸=α)
β(̸=α){ ∑N

F ′′ (ρα )α} αβ αβ αβ αβN ∑N {ϕ′ (rαβ ) ri rj rk rl
1∑ϕ′′ (rαβ ) −+2 α β(̸=α)
rαβ(rαβ )2](2.37)

応力と同様に，各原子位置における原子弾性係数を以下のように定義する．
α
=
Cijkl

[ ∑
} αβ αβ αβ αβ{N1ρ̃′ (rαβ ) ri rj rk rl
′ α
′′ αβ
F (ρ ) ρ̃ (r ) −V /N β(̸=α)
rαβ(rαβ )2αγ αγ }
Nriαβ rjαβ }{ ∑′ αγ rk rl+ F (ρ )ρ̃ (r ) αβρ̃ (r ) αγrrβ(̸=α)
γ(̸=α)
′′

α

{ ∑
N

′

αβ

} αβ αβ αβ αβN {1 ∑ϕ′ (rαβ ) ri rj rk rl+ϕ′′ (rαβ ) −2 β(̸=α)rαβ(rαβ )2](2.38)

これより，系の弾性係数はCijkl =N
1 ∑CαN α ijka(2.39)

2.7 局所格子不安定性解析

16のように原子弾性係数の平均となる．
以上で定義した原子応力，弾性係数から，原子弾性剛性係数は以下で評価できる．
α
α
α
Bijkl
= Cijkl+ (σilα δjk + σjl
δik
α
α
+ σik
δjl + σjk
δil − 2σijα δkl )/2(2.40)α sym
ααWang らによる提案 (33) に従い，Voigt 対称性をもたせた Bijkl≡ (Bijkl+ Blkji
)/2 を
α sym
α
用いて安定性評価を行う．以降では Bijkl
を Voigt 表記した BIJを用いて原子弾性剛性係数を表す．

第3章
無負荷平衡状態における局所格子不安定性α
本章では，原子弾性剛性係数 BIJが負となった「不安定な」原子の分布および系全体としての剛性との関係を明らかにするため，結晶寸法および形状を様々に変えたナノ多結晶体および，アモルファスの無負荷平衡状態について検討する．不安定原子の粒寸法による分布および割合の変化と不安定原子の応力状態，局所密度等について議
論する．

3.1

シミュレーション方法

3.1.1

解析モデル

図 3.1 に模式的に示すように，大きさ 53.89nm × 54.87nm × 1.74nm の薄板状セルを解析領域とする．この薄板状セルを，周期境界条件の下で 2 次元ボロノイ分割および規則分割して多結晶体形状を作成する．その後，薄板状セルの厚さ方向は fcc の [011]方向に固定して，結晶方位を平面内でランダムに回転させた Ni 結晶格子を粒内に充填した．このとき，セル内の結晶粒数を変えることで，平均結晶寸法の違う種々の多結晶体を作成した．図 3.2 に作成したナノ多結晶体の例を示す．原子数は系によりばらつきはあるが約 50 万である．

17

3.1 シミュレーション方法

z direction

1.7

53.89nmy4nm54.87nmWith PBC in
x,y,

x
z

Fig.3.1 Simulation cell.

(a-i) 6 grains(a-ii) 10 grains(a) Voronoi cells(a-iii) 30 grains(b-i) 3×3 division(b-ii) 5×5 division(b) Regular cells(b-iii) 12×12 divisionyx

Fig.3.2 Polycrystals with Voronoi and regular grains.

18

3.1 シミュレーション方法

3.1.2

19

結晶寸法の定義

ボロノイ分割されたナノ多結晶体の代表寸法を以下のように定義した．図 3.3 に模式的に示すように系の各結晶粒を，面積 (体積) の等しい円形結晶粒に置き換えること
で寸法 (直径) を算出し，セル内の全結晶粒の寸法を平均したものをその解析セルの代
表結晶寸法とする．fcc 結晶における原子容 Ω は

Ω=

la3
4

(3.1)

である．ここで，la は初期格子長さである．粒内の原子数を N とすると，粒の体積 V
はV = NΩ(3.2)となる．これを薄板の厚さ t で除せば表面積 S が得られる．したがって，これに等しい面積の円形結晶粒はπ r2 = S =

V
t

(3.3)

となる．これより
√

r=

NΩ
πt

(3.4)

を求め，代表結晶寸法 d = 2r とした．

S

tt

r

V
d

Fig.3.3 Evaluation of representative grain diameter for polycrystals.

3.1 シミュレーション方法

3.1.3

20

初期緩和シミュレーション

作成した各多結晶体について，全方向の垂直応力が零になるように各セル辺長をス
ケーリングしながら 5ps の初期緩和シミュレーションを行った．温度は 10K とし，速度スケーリングにより制御した．(2.1) 式の数値積分は Verlet 法（式 (2.10)）により行
い，積分の時間ステップ ∆t は 1 fs = 10−15 s とした．得られた無負荷平衡状態においαて，原子弾性剛性係数の行列式 det BIJ
を全ての原子について評価するとともに，系の弾性剛性係数の行列式 det BIJ を求めた．

3.1.4

アモルファス構造の作成ならびに構造評価比較のために，同一寸法の薄板セルにアモルファス構造を作成した．まず Ni 原子を

fcc 格子点上に配置した後，全方向周期境界条件を適用して 3000K で 10ps の分子動力学計算を行って溶融させた．その後 10K まで冷却速度 −1 × 1014 K/s で急冷却し，さらに 10K で 5ps の緩和シミュレーションを行った．図 3.4 に 3000K で溶融した状態における動径分布関数，および急冷・緩和シミュレーション後の平衡状態における動径分
布関数を示す．急冷・緩和後の動径分布関数には，アモルファス金属特有の第 2 ピー
クの分岐が現れている．また，この二つのサブピークの位置 r2 , r3 は，第 1 ピークの
位置 r1 に対して，r2 /r1 = 1.75, r3 /r1 = 1.98 となり，サブピークの大小関係と共に Fe
蒸着アモルファス膜の実験値 (34) とほぼ一致する．
Radial distribution function, g(r), 1/nm

80
After melt-quench simulation
Liquid
60

40

20

0

0.2

0.4

0.6

0.8

Distance, r, nm

Fig.3.4 Radial distribution functions.

3.2 シミュレーション結果および考察

3.2

シミュレーション結果および考察

3.2.1

不安定原子の分布21α
図 3.5 に原子弾性剛性係数の行列式 det BIJが負となった不安定原子の分布の例を示す．ボロノイ分割ならびに規則分割の各多結晶体それぞれについて，結晶寸法が同程度のものを寸法が大きい系，小さい系，その中間から選び示している．図中に濃く着α色した原子が det BIJ< 0 となった不安定原子であり，その多くは粒界部分に点在している．

(a) Voronoi polycrystal 1
Grain size = 2.10 nm(b) Voronoi polycrystal 2
Grain size = 10.84 nm(c) Voronoi polycrystal 3
Grain size = 30.60 nm(d) Regular polycrystal 1
Grain size = 2.11 nm(e) Regular polycrystal 2
Grain size = 10.23 nm(f) Regular polycrystal 3
Grain size = 30.59 nm

Unstable atoms : detBαIJ < 0
Stable atoms : detBαIJ > 0

Fig.3.5 Snapshots of unstable(dark shade) and stable(light shade) atoms.

3.2 シミュレーション結果および考察

22α
図 3.6 に各分割の多結晶体における det BIJ< 0 の原子の割合を示す．縦軸に系全体の原子数に対する不安定原子の割合をとり，横軸には 3.1.2 節で示した代表結晶寸法 dをとっている．なお，アモルファスは欠陥構造の究極として結晶寸法 d = 0nm の位置にプロットしている．特に d = 5nm より小さな寸法では微細化するにつれ不安定原子の割合が急激に増加している．Voronoi と規則分割で結晶粒の形状による差異はない．

d → 0 の極限にアモルファスにおける不安定原子の割合があるようにも見えるが，□印の最小の d(=0.55nm) より小さな粒径は，単位格子と同程度となるため，厚さ方向

Ratio of unstable (detBIJα <0) atoms, N unstable / N all

を [011] 方向に固定する条件では作成不能な寸法である．

0.125

Voronoi
Regular
Amorphous

0.1

0.075

0.05

0.025

0

5

10

15

20

25

30

Grain size, d, nm
α< 0) atoms by grain size
Fig.3.6 Change in the ratio of unstable(det BIJ(before loading).

3.2 シミュレーション結果および考察

23α
図 3.7 に原子弾性剛性係数の行列式 det BIJの平均と，系の弾性剛性係数の行列式

det BIJ の結晶寸法による変化を示す．系の弾性剛性係数 BIJ は個々の原子弾性剛性係
α
α
α
数 BIJ
の平均であるが，一般に det BIJ ̸= Σ det BIJ/N である．det BIJの平均は先のα
det BIJ< 0 の原子割合を反映して結晶粒が小さくなるにつれ減少し，粒形状による差異 (Vorononi と規則格子の違い) も見られない．一方，系全体の det BIJ は系全体としての変形に対する剛性 BIJ の変化を反映し，ボロノイ分割した系では内部不均一によるばらつきを生じている．ただし，全体的にはいずれも d =10nm 近傍で最大値をとり，それより大きな d では det BIJ はわずかに減少している．Ni 系合金のナノ多結晶において，Hall-Petch から逆 Hall-Petch へと変化する結晶寸法は 10〜20nm であることが報告されている (35) ．ただし，本解析結果は引張る前の無負荷平衡状態での値であαるため，現時点で直接結びつけられるものではない．図 3.6 と同様，det BIJ ，det BIJの平均いずれも結晶寸法が小さくなるにつれアモルファスのそれに近づいている．ただし，先述のように規則分割可能な最小の d = 0.55nm ではアモルファスのそれとは一致していない．3 次元の解析 (14)−(21) では結晶粒を細分化していくと最終的にはアモル
ファスのそれに一致することが報告されているが (18) ，本解析では 2 次元状の薄板セル
で板厚方向を [011] に固定しているためこのような差を生じたものと考える．
Magnitude of 6×6 determinant, detBIJ, ×1012GPa6

15
Voronoi : global stability

12

Regular : global stability
Amorphous : global stability

9
Voronoi : local stability

6

Regular : local stability

30

Amorphous : local stability

5

10

15

20

25

Polycrystals
VoronoiαΣdetBIJ /N
detBIJ
RegularαΣdetBIJ /N
detBIJ
AmorphousαΣdetBIJ /N
detBIJ

30

Grain size, d, nm

Fig.3.7 Change in global and local stabilities by grain size (before loading).

3.2 シミュレーション結果および考察

24

不安定原子に生じる応力

3.2.2

α
3 次元の解析において，det BIJ< 0 の正値性に基づき「安定な」原子と「不安定な」原子に生じる応力を分けて評価した結果，無負荷平衡状態ではナノ多結晶体中の不安定原子には静水圧引張応力が，アモルファスのそれには静水圧圧縮応力が生じている
ことが報告されている (20) ．本解析で対象とした薄板状セルの，無負荷平衡状態においαて det BIJ
< 0 の原子に働く各方向垂直応力の平均値を，結晶寸法を横軸にとって示したものが図 3.8 である．3.1.3 節で述べたように，系全体の平均応力はいずれも 0 に制御されている．アモルファスならびに，d = 5nm までのナノ多結晶体は，x, y と z 方向の周期境界の違いにも関わらず，応力の異方性は強くなく，3 次元のときと同様にナαノ多結晶体の det BIJ< 0 の原子には引張応力が，アモルファスには圧縮応力が生じてαいる．一方，図 3.6 において det BIJ< 0 の原子割合が収束しはじめる d = 5nm より

粒寸法が大きいナノ多結晶体では異方性が大きくなり，また平均粒径の違いによるばらつきも顕著となっている．このばらつきは，結晶寸法が大きくなるとセル内の結晶粒数が少なくなり，その特定の結晶方位ならびに粒界構造の寄与が強くなったためとα考えられる．また，規則分割した系の det BIJ< 0 の原子に生じる応力はボロノイ分割のそれより常に大きい．これは規則分割した系では x, y 方向に粒界が連続しているためと考えられる．

Stress in unstable (detBα
IJ <0) atoms, Σσyy /N, GPa

Stress in unstable (detBα
IJ <0) atoms, Σσxx /N, GPa

Voronoi
Regular
Amorphous

8

6

4

2

0
0

5

10

15

20

25

30

10
Stress in unstable (detBα
IJ <0) atoms, Σσzz /N, GPa

10

10

Voronoi
Regular
Amorphous

8

6

4

2

0
0

5

10

15

20

25

30

Voronoi
Regular
Amorphous

8

6

4

2

0
0

5

10

15

20

25

Grain size, d, nm

Grain size, d, nm

Grain size, d, nm(a) Stress in x-direction(b) Stress in y-direction(c) Stress in z-directionα< 0) atoms (before loading).
Fig.3.8 Normal stress on unstable(det BIJ

30

3.2 シミュレーション結果および考察

25

粒界および結晶内部の不安定原子

3.2.3

先述のように多結晶体における不安定原子の多くは粒界部分に存在するが，図 3.5 では粒径が小さくなると，結晶粒内にも少なからず不安定原子が存在している．粒界と結晶内部におけるそれぞれの不安定原子について議論するため，CNA(Common Neighbor

Analysis)(36) により粒界原子と結晶粒内原子に分けて検討を行う．CNA により判断した粒界原子数および粒内原子数の粒寸法による変化を図 3.9(a) に示す．また図 3.9(b)には，粒界と結晶内部における不安定原子数をそれぞれの原子数で除して求めた不安
定原子の割合を示す．同程度の寸法を持つ Voronoi と Regular では，粒界および粒内の原子数に大きな差異は見られない．寸法が小さくなると，粒内原子は粒界原子の増加に伴い減少するが，不安定原子の割合は増加する．d = 5nm よりも小さな寸法では
特に高い割合を示しており，d = 2nm の系における結晶粒内に存在する不安定原子の分布を図 3.10 に示す．緑色で示した粒内不安定原子は，桃色の粒界近傍に存在していることが分かる．結晶寸法が大きな系においてもその数は少ないが，同様に結晶内部で不安定となる原子は粒界近傍に位置している．以上より，結晶粒内でも粒界近傍の原子は不安定となる可能性があり，寸法が小さくなると粒界近傍の粒内原子が増える
ため，結晶部分における不安定原子の割合が高まったと考えられる．
0.075

400000Voronoi
300000

NGB : Grain boundary
NIN : In crystal
Regular
NGB : Grain boundary
NIN : In crystal

200000100000

0

5

10

15

20

25

Grain size, d, nm(a) Number of atoms at grain boundary/in crystal

30

unstable
Ratio of unstable (detBIJα <0) atoms, NGB/IN/ NGB/IN

Number of atoms in grain boudary, NGB/IN

500000Voronoi
unstable
NGB/NGB
unstable
NIN
/NIN

0.05

Regularunstable
NGB/NGB
unstable
NIN
/NIN

0.025

0
0

5

10

15

20

25

30

Grain size, d, nm(b) Ratio of detBIJα <0 atoms at grain boundary/in crystal

Fig.3.9 Change in the number of atoms(grain boundary or crystal) and theα< 0) atoms by grain size.
ratio of unstable(det BIJ

3.2 シミュレーション結果および考察

Unstable atoms in crystal
Stable atoms in crystal
Grain boundary atomsyx(a) Voronoi polycrystals(b) Regular polycrystals

Grain size = 2.09 nm
No. of grains : 800

Grain size = 2.11 nm
No. of grains : 841

α< 0) atoms in crystal.
Fig.3.10 Snapshots of unstable(det BIJ

26

3.2 シミュレーション結果および考察

27α
図 3.5 に示したように，すべての粒界原子が det BIJ< 0 となっているわけではない．

「欠陥の中の欠陥」ともいうべき粒界中の不安定原子の局所構造についてさらに議論するため，「安定な」粒界原子と「不安定な」それの局所密度を評価した (図 3.11(a))．また結晶粒内についても同様に評価し図 3.11(b) に示す．ここで，局所密度は各原子からカットオフ半径内に存在する原子体積から求めたものを平均化して示している．安定原子の局所密度は結晶部と粒界の構造の違いから差異が生じているが，結晶粒内，粒界いずれにおいても，不安定原子は同程度の局所密度を示して安定原子のそれよりも低い．これは結晶内部の不安定原子が粒界近傍に存在することから，粒界またはその近傍には局所構造の違いから特に疎な部分が存在し，そこに位置する原子が不安定原子として評価されていることを意味する．不安定原子として明らかになった粒界およびその近傍の特に疎な部分は，負荷時の系の力学応答に何らかの影響を及ぼすことが予想され，これについては第 5 章で検討する．

Stable (detBIJα >0) grain boundary atoms
Unstable (detBIJα <0) grain boundary atoms

Stable (detBIJα >0) atoms in crystal
Unstable (detBIJα <0) atoms in crystal

8.8

8.8
Voronoi

8.7
Local density, Σρ α /NGB, g /cm3

Local density, Σρ α /NGB, g /cm3

8.7
8.6
8.5
8.4
8.3

Regular

8.6

Voronoi

Regular

8.5
8.4
8.3
8.2

8.2
0

5

10

15

20

Grain size, d, nm(a) At grain boundary

25

30

0

5

10

15

20

Grain size, d, nm(b) In crystal

Fig.3.11 Change in the local density by grain size.

25

30

3.3 結 言

3.3

28

結言

ナノ多結晶体の変形メカニズムを原子レベルから明らかにするため，分子動力学法を用いて初期緩和シミュレーションを行い，無負荷平衡状態における局所格子不安定性について検討を行った．得られた結果を要約して以下に示す．
α) が負と評価される不安定原子は，ナノ多結(1) 原子弾性剛性係数の行列式 (det BIJ

晶体中では粒界部分に存在し，微細化による粒界増加に伴いその数は増加する．
α(2) ナノ多結晶体の det BIJの平均 (局所の安定性) は粒界割合に応じて変化し，分割による差異なく微細化に伴い単調に減少する．系の弾性剛性係数の行列式 det BIJ (系の安定性) には分割よる差異が見られ，また各分割とも結晶寸法 10nm 近傍で最α
大値となった．ナノ多結晶体の det BIJ ，det BIJの平均とも微細化によりアモル

ファスに近づくが，いずれも高い値となる．

(3) 無負荷平衡状態においても多結晶体中の不安定原子には正の各方向垂直応力が，アモルファス中の不安定原子には負のそれが生じる．薄板状セルのため，寸法が大きくなると各方向応力に異方性が生じる．また特定の結晶粒および粒界構造の寄与から，系によるばらつきが大きく表れる．寸法が大きな系では規則分割の系
がボロノイよりも高い応力を示した．

(4) 微細化に伴い粒界では不安定原子の割合が高まり，粒界割合の増加に伴って粒界近傍の結晶粒内原子が増えることで，結晶部でも不安定となる原子が増加する．

(5) 無負荷平行状態における不安定原子は，粒界および粒界近傍の結晶部に存在し，他の安定な粒界原子に比べ周囲の構造が疎であることが示された．

第4章
引張変形下の力学応答
ならびに内部構造変化
本章では，寸法の異なるナノ多結晶体およびアモルファスについて単軸引張シミュ
レーションを行い，その力学応答と内部構造変化について検討する。0.2% 耐力により評価した弾性限界，ならびに引張変形中の最大応力について，結晶粒の寸法および形状による違いを議論する．

4.1

シミュレーション方法

前章で得られた無負荷平衡状態のナノ多結晶体とアモルファスについて，y 軸方向に引張るシミュレーションを行った．引張はひずみ制御で行い，全原子間の y 方向距離を均等に拡げることでひずみを与えた．毎ステップ与える微小ひずみ増分は ∆εyy = 1.0 × 10−6であり，ひずみ速度に換算すると 1.0 × 109 s−1 となる．引張中は Poisson 収縮を考慮し，横方向の応力が零 (σxx = σzz = 0) となるようにセル寸法を制御した．計算時間は 30psとし，他の条件については前章と同じとした．計算時間の問題から，すべてのナノ多結晶体への引張を行ったのではなく，表 4.1 に示す同程度の結晶寸法を持つ 5 組の系について行った．無負荷平衡状態における各系の原子配置を CNA により結晶部分 (fcc)，
積層欠陥 (hcp)，それ以外 (other) と分けて着色したものを図 4.1 に示す．

29

4.1 シミュレーション方法

30

Table 4.1 Simulation groups for tensile simulation.
Group Name

GS30

GS20

GS10

GS4

GS2

Model type & No.
No. of grains
Grain size [nm]No. of atoms
Model type & No.
No. of grains
Grain size [nm]No. of atoms

Voronoi 1
430.59
472 752
Regular 1
430.59
470 008

Voronoi 2
920.24
472 787
Regular 2
920.45
472 619

Voronoi 3
30
10.84
472 640
Regular 3
36
10.23
473 270

Voronoi 4
200
4.30
472 528
Regular 4
196
4.38
472 878

Voronoi 5
800
2.09
447 300
Regular 5
841
2.11
471 695(a-i) Voronoi 1(b-i) Voronoi 2(c-i) Voronoi 3(d-i) Voronoi 4(e-i) Voronoi 5(a-ii) Regular 1(b-ii) Regular 2(c-ii) Regular 3(d-ii) Regular 4(e-ii) Regular 5(a) GS30(b) GS20(c) GS10(d) GS4(e) GS2

Hcp structure
Other defective structure ( grain boundary )Fcc structure

Fig.4.1 Snapshots of atoms in the simulation groups with diﬀerent grain size (GS).

4.2 シミュレーション結果および考察

4.2

シミュレーション結果および考察

4.2.1

応力-ひずみ関係

31

引張シミュレーションにより得られた応力-ひずみ関係を，表 4.1 のグループ (GS) 毎に図 4.2 に示した．図中にはアモルファスの応力-ひずみ曲線を細線で示している．図
中の直線は，最小二乗近似により求めた εyy = 0 ∼ 0.01 の範囲の平均勾配である．この勾配，すなわち初期弾性応答は，結晶粒寸法が大きくなるにつれ急になる傾向がある．ただし，Voronoi 分割した系では GS10 の勾配がもっとも急で，GS20，30 と粒が大きくなるにつれ再び勾配がゆるやかになっている．応力-ひずみ曲線は，この初期応答の直線から外れはじめた後は，系によって様々に複雑な挙動を示し，Voronoi および Regular の結晶粒形状，および寸法などに統一的な傾向は見出せない．アモルファ
スは最も小さな (緩やかな) 初期応答を示し，定常流動変形時もいずれの多結晶体よりも低い応力を示した．
Gradient at ε yy = 0~0.01

9Gradient at ε yy = 0~0.01

9

d = 20 nm

RegularVoronoi

6
Voronoi

Regular

3

0.1

0.15

0

0.2

0.05

Voronoi

Regular
3
Amorphous

0.1

0.15

0.2

0

0.05

(a) GS30

Strain, ε yy(b) GS20
Gradient at ε yy = 0~0.01

9(c) GS10
Gradient at εyy = 0~0.01

9

d = 2 nmd = 4 nmVoronoi

6

Regular3

6
RegularVoronoi3

Amorphous

Amorphous

00.1

Strain, εyyStrain, εyyStress, σyy , GPa

0.05

6

Amorphous

Amorphous0

d = 10 nm
Stress, σyy , GPa

Stress, σyy , GPa

6Stress, σyy , GPa

Stress, σyy , GPad = 30 nm

3Gradient at εyy = 0~0.01

9

0.05

0.1

0.15

0.2

0

0.05

0.1

Strain, ε yy

Strain, εyy(d) GS4(e) GS2

Fig.4.2 Stress-strain curves.

0.15

0.2

0.15

0.2

4.2 シミュレーション結果および考察

4.2.2

32

0.2% 耐力および最大応力の結晶寸法依存性

非線形変形の開始点を，通常の引張試験における 0.2% 耐力に従って評価した．図 4.2の初期応答直線を εp = 0.002 シフトさせ，その直線と応力-ひずみ曲線の交点から 0.2%
耐力を求めた．また，εyy = 0 ∼ 0.2 の範囲内で示した最大応力を評価し，0.2% 耐力とともに粒寸法を横軸にとって図 4.3 に示した．d = 30nm の系 (GS30) を除き 0.2% 耐力には形状に違いによる差異は見られず，また d = 5 ∼ 20nm の範囲では寸法依存性も小さい．また Voronoi では d = 10nm で最大応力がピークを示し，この点で Hall-Petchから逆 Hall-Petch へと移行しているようにみえるが，Regular では GS20 の最大応力が最も低くなりそのような傾向は見られない．特に，d = 30nm の Regular の系は 0.2%
耐力，最大応力ともに Regular 中で最大となっており，結晶粒の少ない系におけるばらつき (18) によるものと考えられる．d = 5nm より小さな寸法では 0.2% 耐力，最大応
力に Voronoi と Regular で差異はなく，ともに単調減少してアモルファスに近づいているが，アモルファスのそれよりも高い値をとっている．

Stress, σyy , GPa

9

6Voronoi
3

Regular

Amorphous
0

σ0.2%
σmax

Regular

5

σ0.2%
σmax
Amorphous
σ0.2%
σmax

Voronoi

10
15
20
25
Grain size, d, nm

30

Fig.4.3 Change in the 0.2% proof stress and maximum stress by grain size.

4.2 シミュレーション結果および考察

4.2.3

33

引張変形時の原子構造変化

図 4.2 に示したように，応力-ひずみ曲線は複雑に変動している．応力変動と内部構造変化を明らかにするため，d = 30nm の Regular の系について詳細に検討した．図4.4 にその応力-ひずみ曲線を再掲し，(a)∼(i) 各点における原子配置を，CNA により判断した原子構造 (hcp，fcc，その他の欠陥構造) に従って着色し図 4.5 に示した．ここ
で，図上部に示したように，系全体よりも一回り小さな四角 (図中の最上部) の範囲を拡大して示している．なお図 4.5 の (a)∼(i) は図 4.4 のひずみに対応している．図 4.5(a)より，粒界部分は青の hcp 原子とピンクの欠陥原子からなる．0.2% 耐力のひずみ (図4.5(b)) において，粒界を起点として積層欠陥 (hcp) が結晶内部に発生している．また，粒内の積層欠陥の端点はピンク色に着色されており，部分転位の芯である．(b) から(f) までの遷移領域では粒内に多数の積層欠陥を生じている．ここで，図 4.4(c)∼(d) の間で応力停滞→最上昇しているが，図 4.5 の (c) をみると右下の粒内に積層欠陥を生じておらず，(d) ではみられる．このことから，結晶方位の関係でもっとも変形抵抗の高かった右下の結晶粒に転位を発生する応力が図 4.4 の (d) のピークであったことが分かる．最大応力を示す (f) まで，他の結晶粒では 2 つのすべり系が駆動しているが，右下の粒は 1 すべり系のみである．また (f) までは初期の粒界構造を保っている．一方，最大応力を示した (f) 以降は，図中 (f)∼(i) に白抜き矢印で示したように粒界が変形に伴い移動している．

(a) (b)(c) (d)(e) (f)

(g)

(h)

(i)

0.15

0.2

Stress, σyy , GPa

9

6

3

0

0.05

ε0.2%

0.1
Strain, εyy

Fig.4.4 Stress-strain curves on regular polycrystal of GS30.

4.2 シミュレーション結果および考察

Hcp structure
Other defective structure
Fcc structureyx(a) ε yy = 0(b) ε yy = 0.017(c) εεyy = 0.028

0.2% proof stress(d) εεyy = 0.042(e) εεyy = 0.064(f) εεyy = 0.083
Maximum stress(g) εεyy = 0.13434(h) εεyy = 0.156(i) εεyy = 0.2

Fig.4.5 Snapshots of hcp, fcc and other defective atoms.

4.3 結言

4.3

35

結言

引張変形下の応力-ひずみ応答の寸法依存性，ならびに内部構造変化を明らかにするために，ナノ多結晶体およびアモルファスについて単軸引張シミュレーションを行った．得られた結果を要約して以下に示す．

(1) εyy = 0 ∼ 0.01 の初期弾性応答は，規則分割したナノ多結晶体は結晶寸法が小さ
くなるにつれ変形抵抗が減少 (勾配が緩やかになる) 傾向が見られた．

(2) Voronoi 分割した系では，d ≤ 10nm では (1) と同じ傾向であったが，10nm よりも大きな寸法では結晶寸法が大きくなっても初期勾配が減少した．

(3) アモルファスは最も小さな初期勾配を示した．
(4) 初期勾配の直線から逸脱し始めた後の応力-ひずみ応答は，系によって様々に複雑な応答を示し，結晶寸法等による統一的な傾向は認められなかった．

(5) 0.2% 耐力および最大応力と結晶粒径の関係を調べた結果，0.2% 耐力には粒形状(Voronoi と Regular) の影響が少なく，d = 30nm の系を除きほぼ一致した．また最大応力は，ボロノイ分割の多結晶体では寸法 10nm でピークを示し，Hall-Petchから逆 Hall-Petch への移行が見られたが，Regular 分割ではそのような傾向は見られなかった．

(6) CNA 解析により引張変形下の原子構造変化を観察し，0.2% 耐力近傍から積層欠陥が発生していること，最大応力を示した点以降では粒界構造自体の変化を生じることを示した．

第5章
引張変形下の力学応答と不安定原子の力学状態
本章では，結晶寸法約 10nm のナノ多結晶体について，引張変形下の力学応答と不安定原子の割合や分布の変化について検討する．不安定原子に生じる応力および割合の変化を局所構造ごとに評価し，系の応力-ひずみ曲線に表れる変化との関係を議論する．なおボロノイ分割および規則分割した多結晶体に加え，新たに正六角形分割した多結晶体についても検討を行う．

5.1

解析モデルおよびシミュレーション方法

解析モデルには 4.1 節で示した GS10 の多結晶体 (Voroni，Regular) を用いた．また両者の中間の結晶形状として，大きさ 51.96nm × 54.0nm × 1.74nm のセルを粒寸法約10nm の正六角形分割した系も対象とした．各多結晶体の初期構造を表 5.1 に示す．正六角形分割の多結晶体は 3.1.1 節と同様の手順で作成し，3.1.3 節に示した条件で初期緩和を行った．無負荷平衡状態の各多結晶体を初期状態とし，前章と同様の条件で単軸引張シミュレーションを行った．

36

5.1 解析モデルおよびシミュレーション方法

37

Table 5.1 The tensile simulation models of polycrystals.
Type of Polycrystal

Voronoi

RegularHexagon

Simulation Cell

Cell Size

53.89nm×54.87nm×1.74nm

51.96nm×54.0nm×1.74nm

Grain Size

10.84nm10.23nm10.91nmNo. of Atoms

472 640473 270448 385

No. of Grains30

36 (6×6)30

5.2 シミュレーション結果および考察

5.2

シミュレーション結果および考察

5.2.1

応力-ひずみ関係

38

図 5.1 に得られた応力-ひずみ曲線を示す．図中には初期弾性応答として，最小二乗近似により求めた εyy = 0 ∼ 0.01 の範囲の平均勾配を破線で示す．Voronoi および Regularについては図 4.2(c) と同じものである．初期の線形応答には，結晶粒形状による差は小さい．その後，ゆらぎながらの応力上昇は，ボロノイ分割の系が最も高い応力を示し，規則分割および六角形分割の多結晶体では，εyy = 0.025 近傍で生じる小さなピー
クにわずかに差が見られるが，εyy = 0.06 近傍までの応力変動に大きな差はない．

Gradient at εyy = 0~0.01

9Stress, σyy , GPad = 10 nm

6VoronoiHexagon

3

0

0.05

0.1

Regular

0.15

Strain, εyy

Fig.5.1 Stress-strain curves.

0.2

5.2 シミュレーション結果および考察

39

不安定原子に生じる応力

5.2.2

α
α
3 次元の解析では，det BIJ< 0 の不安定原子に作用する応力は，det BIJ> 0 の安定原子のそれよりも高く「大きく変形」していることが報告された (20) ．個々の不安定原
α子に生じる引張方向の原子応力 σyy
を平均して評価し，系の応力-ひずみ曲線と併せて

図 5.2 に示す．安定原子に生じる原子応力を平均したものは，系の応力-ひずみ曲線とほぼ等しいので省略している．前章で述べたように不安定原子の多くは粒界に存在するので，結晶内部に比べて周囲の原子が「疎」な状態にあり，εyy = 0 で高い静水圧引
張応力を示す (20) ．ここで，3 次元のナノ多結晶体 (d = 3.5 ∼ 5.3nm) では不安定原子の初期応力は 2∼3GPa であったのに対し，2 次元の本解析では 4.5∼6GPa と倍程度大きくなっている．引張変形時の応力も 3 次元と全く異なる挙動を示している．安定原子より常に高い応力を示しているのは同じであるが，3 次元では系の応力-ひずみ関係を上にシフトさせたような応答であったのに対し，2 次元では 0.2% 耐力を示すひずみ近傍で不安定原子の応力は最大となり，その後大きく減少した．前章の図 4.5 において，0.2% 耐力のひずみから結晶内部に部分転位を生じていることからわかるように，不安
定原子の応力ピークは粒界部分での変形が最大となった「弱い部分の変形限界」である．その後，粒界に積層欠陥が導入されると，厚さ方向に粒界が貫通している 2 次元周期セルでは応力の緩和が容易となり，3 次元では見られない顕著な応力低下を生じたものと考える．

6

3

9Stress, σyy , GPa

9Stress, σyy , GPa

Stress, σyy , GPa

9

6

3

Unstable (detBIJα <0) atoms
System
0

εp
ε0.2%

0.05

0.1
Strain, εyy(a) Voronoi

0.15

6

3

Unstable (detBIJα <0) atoms
System
0.2

0

εp
ε0.2%

0.05

0.1
Strain, εyy(b) Regular

0.15

Unstable (detBIJα <0) atoms
System
0.2

0

εp
ε0.2%

0.05

0.1
Strain, εyy

0.15(c) Hexagonα< 0) atoms.
Fig.5.2 Stress-strain curves on system and unstable (det BIJ

0.2

5.2 シミュレーション結果および考察

40

結晶構造毎の応力変化

5.2.3

CNA で fcc，hcp，その他の欠陥構造 (other) と判定した原子の割合変化を図 5.3 に示す．図中には 0.2% 耐力を示すひずみを破線で示している．いずれの多結晶体も，0.2%
耐力を示すひずみ近傍から積層欠陥 (hcp) が増加し，結晶内部の fcc 構造原子が減少する．ひずみ 0.1 近傍で積層欠陥の増加が飽和し，fcc 構造原子の減少も収束しているが，これは図 5.1 の応力-ひずみ曲線で応力上昇が飽和し，定常流動変形する領域に対応する．
αdet BIJの正負ではなく，これらの局所構造別に平均応力を評価しその変化を示したものが図 5.4 である．いずれも εyy = 0 の初期状態では hcp 構造原子が高い応力を示している．図 4.5(a) で示したように，粒界構造は hcp およびその他の欠陥構造の原子から構造されているが，粒界においても高い引張応力を生じている部分とほぼ零の部分があることが分かる．引張ひずみを与えると，0.2% 耐力を示すひずみにおいて，hcpとその他の欠陥構造の原子の応力-ひずみ応答に折れ曲がりを生じる．前節で述べたよ
うに，初期の粒界構造を保ったまま変形を吸収する時の，粒界部分の変形限界に対応するものと考える．その後は hcp と fcc 原子の担う応力との差が縮まり，その他の欠陥原子は常に低い応力を示す．図 5.3 に示したように，0.2% 耐力を示すひずみ以降は

0.8
fcc
0.6

0.4

hcp

0.2

other
0

εp

0.05

0.1
Strain, εyy(a) Voronoi

0.15

0.2

1

Ratio of the atoms (hcp, fcc, other), NA/Nall

1

Ratio of the atoms (hcp, fcc, other), NA/Nall

Ratio of the atoms (hcp, fcc, other), NA/Nall

hcp 原子が増加するが，これは粒内に生じた積層欠陥によるものである．したがって

0.8

fcc
0.6

0.4
hcp

0.2

other
0

εp

0.05

0.1
Strain, εyy(b) Regular

0.15

0.2

1

0.8
fcc
0.6

0.4
hcp

0.2

other
0

εp

0.05

0.1
Strain, εyy(c) Hexagon

Fig.5.3 Change in the ratio of the fcc, hcp, other atoms judged by CNA.

0.15

0.2

5.2 シミュレーション結果および考察

9

6
hcp
fcc

other

3

Stress, σyy , GPa

9Stress, σyy , GPa

Stress, σyy , GPa

9

41

6
hcp
fcc
3

6
hcp

fcc
other

3

other

0

εp

0.05

0.1

0.15

0.2

0

εp

0.05

0.1

0.15

0.2

0

εp

0.05

0.1

Strain, εyyStrain, εyyStrain, εyy(a) Voronoi(b) Regular(c) Hexagon

0.15

Fig.5.4 Stress-strain curves on the fcc, hcp, other atoms judged by CNA.
hcp 原子は粒界に属するものと粒内の積層欠陥の 2 つが存在する．後者の応力状態は周囲の fcc 原子と大きな差はない．このため，粒内に積層欠陥が生じはじめる 0.2% 耐力で hcp の応力が折れ曲がり，hcp 原子の増加とともに fcc の応力に近づく．その他の欠陥原子 (other) についても，粒内に多数の積層欠陥を生じはじめると，粒界部分の応力が緩和するため応力の折れ曲がりを生じたものと考えられる．

5.2.4

結晶構造と不安定原子

図 5.3 では，fcc，hcp，その他の欠陥構造 (other) の原子数の割合変化を示したが，αdet BIJが負となった原子の割合変化を，fcc，hcp，other で分けて示したものが図 5.5である．図中には系の応力-ひずみ曲線も併記している．εyy = 0 における割合をみてわα
かるように，det BIJ< 0 の原子のほとんどは fcc でも hcp でもない乱れた構造原子であることがわかる．0.2% 耐力のひずみ近傍まで不安定原子の割合はいずれもほとんど変化がなく，その後 other の不安定原子が顕著に増加していることから，fcc でも hcpαα< 0 の原< 0 となっていることがわかる．fcc，hcp の det BIJ
でもない原子が det BIJ

子割合は，変形後期に増加しているが other に比べると小さい．以上総合すると，0.2%
α< 0 の原子増加は発生した部分転位の芯に対応するもの
耐力のひずみ以降の det BIJと考えられる．しかしながら，ここまでの議論でいくつかの矛盾を生じる．それはα< 0 の原子に生じる応力が大きくて，図 5.4 では other の原子に生
図 5.2 では det BIJ

0.2

42

9
0.01

σyy(system) : right axis

0.0075

other
6

other

0.005
other

σyy(system) : right axis

σyy(system) : right axis

30.0025
fcc

fcc

0

εp

0.05

0.1
Strain, εyy

hcp
0.15

(a) Voronoi

0.2 0

εp

0.05

0.1
Strain, εyy

hcp
0.15

(b) Regular

fcc

0.2 0

εp

0.05

0.1
Strain, εyy

hcp
0.15

(c) Hexagon

αFig.5.5 Change in the ratio of unstable (det BIJ< 0) atoms in the fcc, hcp,
other atoms judged by CNA.

じる応力が小さいこと，α
図 5.5 の other の det BIJ< 0 原子の増加により転位芯の増加を示唆したが，図 5.3 の other は増えていないこと，である．これらについて次節で詳細に検討する．

5.2.5

引張軸に垂直な粒界原子の力学状態

安定・不安定原子，fcc・hcp・その他の欠陥原子 (other) の各力学状態をさらに解明
αすべく，図 5.4 で示した各構造原子の応力を，さらに det BIJの正負で分けて評価したものを図 5.6 に示す．図中には系の応力-ひずみ曲線ならびに 0.2% 耐力のひずみを併α記している．図 5.5 で示したように，fcc，hcp で det BIJ< 0 の原子は少ないため，応α
力のゆらぎが大きい．fcc 構造でかつ det BIJ> 0 の原子の応力-ひずみ曲線は，系のそα< 0 の原子でれとほぼ一致している．最も高い応力を示しているのは hcp で det BIJ

ある．hcp の中で安定な原子の応力は，0.2% 耐力のひずみ以降は系の応力と一致しており，先の議論のように結晶粒内の安定な積層欠陥に対応するものと考えられる．このように，図 5.4 で示した hcp 原子の高い応力は，少数であるが極めて高い応力の不安定 hcp 原子によるものであり，図 5.2 の不安定原子の応力とも矛盾しない．ただし，

Hexagon については 0.2% 耐力のひずみを超えても不安定 hcp 原子の応力は減少しておらず，図 5.2 の傾向をもたらしたのは other の不安定原子である．ここで，Voronoiの不安定 hcp 原子の応力も，Regular に比べると応力減少が比較的小さい．Hexagon は0
0.2

Stress, σyy , GPa

Ratio of unstable (detBIJα <0) atoms,

NAinst / N all

5.2 シミュレーション結果および考察

5.2 シミュレーション結果および考察

43

引張方向に垂直な粒界が存在せず，また Voronoi も Regular に比べると少ない．すなαわち，hcp で det BIJ< 0 の原子が引張方向に垂直な粒界で高い引張応力を示す局所構

造である可能性がある．
detBIJα <0 atoms

5

System & detBIJα >0 atoms

0

15

10

10

detBIJα >0 atoms
5

System

0
-5

-5
-10

15

Stress, σyy , GPa

detBIJα <0 atoms

10

Stress, σyy , GPa

Stress, σyy , GPa

15

0 εp

0.05

0.1

0.15

0.2

-10

detBIJα <0 atoms

5
detBIJα >0 atoms
0

System

-5

0 εp

0.05

Strain, εyy(a-i) Fcc

0.1

0.15

0.2

-10

0 εp

0.05

0.1

0.15

Strain, εyyStrain, εyy(a-ii) Hcp(a-iii) Other

0.2(a) Voronoi
detBIJα <0 atoms

15

15

15

5
System & detBIJα >0 atoms
0
-5

-10

10
5
detBIJα >0 atoms
0
System

-5

0 εp

0.05

0.1

0.15

0.2

Stress, σyy , GPa

detBIJ <0 atoms

10

Stress, σyy , GPa

Stress, σyy , GPa

detBIJα <0 atoms
α

-10

0 εp

(b-i) Fcc

System5
detBIJα >0 atoms
0
-5

0.05

Strain, εyy

10

0.1

0.15

0.2

-10

0 εp

0.05

0.1

0.15

Strain, εyyStrain, εyy(b-ii) Hcp(b-iii) Other

0.2(b) Regular15

detBIJα <0 atoms

15

15

5
System &detBIJα >0 atoms
010Stress, σyy , GPa

10Stress, σyy , GPa

Stress, σyy , GPa

detBIJα <0 atoms

detBIJα >0 atoms

50

detBIJα <0 atoms

10
5

detBIJα >0 atoms
0
System-5-10

System

-5

0 εp

0.05

0.1

Strain, εyy(c-i) Fcc

0.15

0.2

-10

0 εp

0.05

-5

0.1

0.15

0.2

-10

0 εp

0.05

0.1

0.15

Strain, εyyStrain, εyy(c-ii) Hcp(c-iii) Other(c) Hexagonα< 0) atoms in
Fig.5.6 Stress-strain curves on system and unstable (det BIJ
the fcc, hcp, other atoms judged by CNA.

0.2

5.3 結言

5.3

44

結言

引張変形下における不安定原子の役割ならびに系の応答との関係を明らかにすることを目的として，結晶寸法約 10nm のナノ多結晶体について単軸引張シミュレーショ
α
ンを行い，det BIJの正値性ならびに CNA による結晶構造の観点から詳細に検討した．

得られた結果を要約して以下に示す．
αα(1) det BIJの正負で応力を分けて評価した結果，3 次元の場合と同様，det BIJ<0の原子の平均応力は静水圧引張状態にあることが示された．一方，引張下では0.2% 耐力で不安定原子の応力が最大となり，その後減少するという 3 次元とは異なる挙動を示した．

(2) CNA により fcc，hcp，その他の欠陥 (other) と判定された原子の割合変化を調べたところ，0.2% 耐力のひずみ近傍から hcp が増加し，fcc が減少した．その他の欠陥原子は引張中ほぼ一定であった．

(3) (2) で分類した原子毎に応力変化を調べた結果，εyy = 0 で静水圧引張応力を示したのは hcp 原子であった．fcc，other は初期応力は 0 である．0.2% 耐力のひずみ
で hcp と other が同時に折れ曲がりを示し，その後 hcp の応力は fcc のそれに一致した．other はそれらより常に低い応力を示した．
α(4) det BIJ< 0 の原子の割合変化を，(2) の分類で分けて調べたところ，εyy = 0 でα< 0 の原子はほとんど 0 であった．引は other が最も多く，fcc，hcp で det BIJ

張時，0.2% 耐力ひずみまではいずれも変化せず，それ以降は other の不安定原子が顕著に増加した．
α< 0 の正負で分け(5) (1)∼(4) は矛盾するようであるが，(2) の分類をさらに det BIJてそれぞれの原子に生じる応力を分けて評価 (fcc，hcp，other× 安定・不安定の6 通り) した所，少ない不安定 hcp 原子に最も高い応力を生じていることが明らかになり，矛盾なく説明できることが示された．

5.3 結言45(6) 引張方向に垂直な粒界をもたないナノ多結晶体では 0.2% 耐力ひずみ時の不安定
hcp 原子の応力減少がみられなかった．このことは，不安定 hcp 原子が引張方向に垂直な粒界に特異な力学状態である可能性を示す．

第6章
結論
本研究では，ナノ多結晶体の変形メカニズムについて，格子不安定性という独自の視点から解明することを目的として，分子動力学法により結晶寸法や粒形状を変えたナノ多結晶の引張シミュレーションを行った．以下に，得られた結果を総括する．
第 2 章では解析手法について説明した．始めに分子動力学法の基礎方程式を示し，原子間相互作用の評価に用いた原子埋め込み法ポテンシャルの概要について説明した．
さらに，原子弾性剛性係数を用いて局所の安定性を評価する局所格子不安定性解析について概説した．
第 3 章では，結晶寸法および粒形状の違う種々のナノ多結晶体を作成して初期緩和α
シミュレーションを行い，原子弾性剛性係数 BIJが負となった「不安定な」原子の力

学状態および局所構造の変化と，系の剛性との関係について検討した．また比較のため，溶融-急冷シミュレーションによりアモルファスを作成し，同様の検討を行った．
α
その結果，ナノ多結晶体の det BIJの平均 (局所の安定性) は粒界割合に応じて変化し，粒形状による差異なく微細化に伴い単調に減少すること，系の弾性剛性係数の行列式

det BIJ (系の安定性) には粒形状による差異が見られるが，ボロノイ分割および規則分αの割ともに結晶寸法 10nm 近傍で最大値となること，ナノ多結晶体の det BIJ ，det BIJ

平均とも微細化によりアモルファスに近づくが，いずれも高い値となること，などを明らかにした．また，無負荷平衡状態においても不安定原子には，多結晶体中では正の各方向垂直応力が，アモルファス中では負のそれが生じることを示し，寸法が大きなナノ多結晶体に生じる各方向応力の異方性および系によるばらつきを明らかにした．
さらに，微細化に伴って粒界および結晶内の不安定原子の割合が高まること，無負荷46

47
平行状態における不安定原子は，粒界および粒界近傍の結晶部に存在し，他の安定な粒界原子に比べ周囲の構造が疎であることを示した．
第 4 章では，引張変形下の応力-ひずみ応答の寸法依存性，ならびに内部構造変化を明らかにするために，寸法の異なるナノ多結晶体およびアモルファスについて単軸引張シミュレーションを行った．初期弾性応答の応力-ひずみ勾配は，規則分割したナノ多結晶体では結晶寸法が小さくなるにつれ減少傾向を示すこと，ボロノイ分割した系
では寸法 10nm で最大となって 10nm よりも大きな寸法では初期勾配が減少すること，
アモルファスは最も小さな初期勾配を示すことを明らかにした．0.2% 耐力および最大応力と結晶粒径の関係を調べた結果，0.2% 耐力には粒形状の影響が少なく，最大応力にはボロノイ分割の多結晶体で，寸法 10nm 近傍での Hall-Petch から逆 Hall-Petch への移行が見られた．さらに CNA 解析により引張変形下の原子構造変化を観察し，0.2%
耐力近傍から積層欠陥が粒内に発生していること，最大応力を示した点以降では粒界構造自体の変化を生じることを示した．
第 5 章では，引張変形下における不安定原子の役割ならびに系の力学応答との関係を明らかにするため，結晶寸法約 10nm のナノ多結晶体について単軸引張シミュレーα
ションを行い，det BIJの正値性ならびに CNA による結晶構造の観点から詳細に検討
αした．det BIJの正負だけでなく，fcc，hcp，その他の欠陥 (other) で分類すると，(1)αα
無負荷では det BIJ< 0 の原子の多くは other であり，fcc，hcp で det BIJ< 0 の原子は少ない，(2)fcc，other の不安定原子は無負荷で応力 0 であり，3 次元の解析で報告した静水圧引張応力を示すのは，小数ではあるが粒界の不安定 hcp 原子である，(3)0.2%
耐力ひずみを超えると不安定で other な原子割合が増える，(4)other 原子の割合自体は引張変形中ほぼ一定で，0.2% 耐力ひずみを超えると fcc の割合が減り，それに対して

hcp が増える．これは粒内の積層欠陥増加に対応する，(5)(3) および (4) より，0.2% 耐力ひずみ以降に不安定となるのは主として「非結晶」な部分である，(6) 不安定 hcp 原子に生じる応力は常に高く，かつ 0.2% 耐力ひずみでピークを示し急減した．これは 3次元では見られない現象であった，(7) 引張方向に垂直な粒界をもたないナノ多結晶体
では (6) の不安定 hcp 原子の応力低下がみられず，引張方向に垂直な粒界を特徴づけるものである可能性がある，等を明らかにした．

参考文献(1) 川合知二, ナノテクノロジー入門, (2002), オーム社.
(2) 平尾一之, 基礎から学ぶナノテクノロジー, (2003), 東京化学同人.
(3) 平尾一之, ナノテクへの挑戦 人と技術で究極のものづくり, (2004), 東京化学
同人.

(4) Birringer, R., Gleiter, H., Klein, H. P. and Marquardt, P., Nanocrystalline Materials, An Approach to a Novel Solid Structure with Gas-Like Disorder?, Phys.
Lett., 102A-8(1984), 365–369.
(5) A. H. Chokshi, A. Rosen, J. Karch and H. Gleiter, On the Validity of theHall-Petch Relationship in Nanocrystalline Materials, Scripta Metall Mater, 23–
10(1989), 1679–1683(6) G. E. Fougere, J. R. Weertman, R. W. Siegel and S. Kim, Grain-Size Dependent
Hardening and Softening of Nanocrystalline Cu and Pd, Scripta Metall Mater,
26–12(1989), 1879–1883(7) S. X. McFadden, R. S. Mishra, R. Valiev, A. P. Zhilyaev and A. K. Mukherjee,
Low-Temperature Superplasticity in Nanostructured Nickel and Metal Alloys,
Nature, 398–6729(1999), 684–686.
(8) 伊藤良規, 辻 伸泰, 斎藤好弘, 宇都宮裕, 左海哲夫, 超微細結晶粒アルミニウムの焼鈍にともなう組織と機械的性質の変化, 日本金属学会誌, 64–6(2000), 429–437(9) Schiøtz, J., Tolla,F.D.D. and Jacobsen, K. W., Softening of Nanocrystalline
Metals at Very Small Grain Size, Nature, 391(1998), 561–563.

48

参考文献

49

(10) 中谷彰宏, 北川浩, 下川智嗣, 分子動力学法によるナノ多結晶体の力学特性と引張荷重下における非弾性変形機構の検討, 日本機械学会論文集, A 編, 66–643(2000), 435–441.
(11) H.Van Swygenhoven, M. Spaczér and Jacobsen, K. W., Characterisation of theMicrostructure of Nanophase Ni : A Molecular Dynamics Simulation Study,
Nanostruct Mater, 12(1999), 631–636.
(12) T. Shimokawa, A. Nakatani and H. Kitagawa, Grain-size Dependence of the Relationship between Intergranular and Intragranular Deformation of Nanocrystalline Al by Molecular Dynamics Simulations, Phys. Rev. B, 71–224110 (2005),
1–8.
(13) 屋代如月, 北村隆行, 大谷隆一, 野田茂之, 面心立方晶の転位発生の分子動力学シ
ミュレーションおよび前駆局所ひずみ集中の格子不安定条件による特性解析, 日本機械学会論文集, A 編, 64–626 (1998), 2583–2589.

(14) 屋代如月, 北村隆行, 大谷隆一, 引張軸に対して直角方向に拘束を受けるニッケル
結晶のへき開破壊に関する原子レベル解析, 日本機械学会論文集, A 編, 63–615(1997), 2380–2386.
(15) 屋代如月, 金属結晶の変形と破壊の分子動力学シミュレーション並びに局所格子
不安定性解析, 学位論文, 京都大学, (1998).

(16) 屋代如月, 冨田佳宏, 転位の発生および運動時の局所格子不安定性に関する分子動力学研究, 日本機械学会論文集, A 編, 67–656 (2001), 678–684.

(17) 屋代如月, 西村正臣, 冨田佳宏, 引張変形下のアモルファス金属における局所格
子不安定性：分子動力学による検討, 材料, 54–10 (2005), 1053–1059.

(18) 屋代如月, 西村正臣, 冨田佳宏, 不均一系の大域的・局所的安定性：ナノ多結晶体とアモルファス金属の局所格子不安定性解析, 日本機械学会論文集,A 編, 73–726(2007), 216–223.

参考文献

50

(19) 屋代如月, 西村正臣 久馬雅彦, 冨田佳宏, 不均一系の大域的・局所的安定性：外部境界条件に関する局所格子不安定性解析, 日本機械学会論文集,A 編, 73–725(2007), 66–72.
(20) 屋代如月, 西村正臣, 冨田佳宏, 単軸引張下のナノ多結晶体およびアモルファス
金属の不安定挙動：局所格子不安定解析, 材料, 57–2 (2008), 掲載予定.

(21) 屋代如月, 西村正臣, 冨田佳宏, アモルファス金属の局所格子不安定解析：Ｉ．単
元系アモルファスにおける不安定原子の役割, 材料, (2008), 掲載予定.

(22) Li, J., Van Vliet, K. J., Zhu, T., Yip, S. and Suresh, S., Atomistic Mechanisms
Governing Elastic Limit and Incipient Plasticity in Crystals, Nature, 418 (2002),
307–310.
(23) 北川浩, 渋谷陽二, 北村隆行, 中谷彰宏, 初心者のための分子動力学法, (1997), 養賢堂.

(24) 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
(25) Daw, M. S. and Baskes, M. I., Semiempirical, Quantum Mechanical Calculationof Hydrogen Embrittlement in Metals, Phys. Rev. Lett., 50–17 (1983), 1285–
1288.
(26) Daw, M. S. and Baskes, M. I., Embedded–Atom Method : Derivation and Application to Impurities, Surfaces and Other Defects in Metals, Phys. Rev. B,
29–12 (1984), 6443–6453.
(27) Voter, A. F. and Chen, S. P., Accurate Interatomic Potentials for Ni, Al andNi3 Al, Mat. Res. Soc. Symp. Proc., 82 (1987), 175–180.
(28) Voter, A. F., Chapter 4 The Embedded Atom Method, John Wiley & Sons,(1994), 77–90.

参考文献

51

(29) Foiles, S. M., Calculation of the Surface Segregation of Ni-Cu Alloys with theUse of the Embedded-Atom Method, Phys. Rev. B, 32–12 (1985), 7685–7693.
(30) Rose, J. H., Smith, J. R., Guinea, F. and Ferrante, J., Universal Features of theEquation of State of Metals, Phys. Rev. B, 29–6 (1984), 2963–2969.
(31) 戸田盛和, 松田博嗣, 樋渡保秋, 和達三樹, 液体の構造と物性 (1976), 岩波書店.
(32) Wallace, P. C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
(33) Wang, J., Li, J., Yip, S., Phillpot, S. and Wolf, D., Mechanical Instabilities of
Homogeneous Crystals, Phys. Rev. B, 52–17 (1995), 12627–12635.
(34) 増本健, アモルファス金属の基礎, (1982), オーム社.
(35) 渡邊忠雄, 結晶粒微細化への新アプローチ, 日本金属学会編, (2000).
(36) J. D. Honeycutt, and H. C. Andersen, Molecular Dynamics Study of Melting
and Freezing of Small Lennard-Jones Clusters, Journal of physical Chemistry,82 (1987), 4950–4963.

関連講演論文
 樋口 昌宏・西村 正臣・屋代 如月・冨田 佳宏, 不均一系の大域的・局所的安定性：
ボロノイ・規則分割したナノ多結晶体の局所格子不安定性解析, 日本機械学会第19 回計算力学講演会講演論文集, (2006), 637-638.
 屋代 如月・西村 正臣・樋口 昌宏・冨田 佳宏, ナノ多結晶体およびアモルファス金属の局所不安定解析, 第 12 回分子動力学シンポジウム講演論文集, (2007), 56-61.

 屋代 如月・西村 正臣・樋口 昌宏・冨田 佳宏, せん断変形を受けるナノ多結晶体
およびアモルファス金属の局所格子不安定性解析, M&M2007 材料力学カンファ
レンス講演論文集, (2007), 500-501.

 樋口 昌宏・西村 正臣・屋代 如月・冨田 佳宏, ナノ多結晶体の局所格子不安定性解
析：結晶粒形状の影響, 日本機械学会第 20 回計算力学講演会講演論文集, (2007),

255-256.
 屋代 如月・西村 正臣・樋口 昌宏・冨田 佳宏, アモルファス金属における不安定原子の役割：局所格子不安定性解析, 第 51 回日本学術会議材料工学連合講演会講
演論文集, (2007), 222-223.

 Kisaragi YASHIRO・Masaomi NISHIMURA・Masahiro HIGUCHI・Yoshihiro TOMITA,
Local Lattice Instability Analysis on Nano-polycrystalline and Amorphous Metals: Role of Unstable Atoms under Tension, Compression and Shear, The 2nd
International Conference on Heterogeneous Materials Mechanics, (2008), 講演
予定.

52

謝辞
本研究を遂行するにあたり，終始懇切丁寧な御指導賜りました冨田佳宏教授に心より感謝致します．研究生活を通じて，健全な見識と料簡について数多く学ぶことができ，今後の貴重な指針を享受することができました．また，研究活動を直接御指導く
ださり，迷妄な著者を座視することなく，一貫して熱心な御教示を頂きました屋代如
月准教授に深く感謝致します．多くの的確かつ明敏な考察と，透徹した論理を導出する諸過程に触れることができました．この体験から大いに刺激を受けると共に，改めて自己鍛錬に励む必要性を感じました．浅学非才の著者が本修士論文を完成することができたのも，両先生のお陰であり，深甚に感謝をささげる次第であります．
また，研究内容について有益な議論と御助言を頂き，研究活動促進に貢献下さいました長谷部忠司准教授，高木知弘准教授 (現 京都工芸繊維大学) に心から感謝致します．そして研究環境づくりに御尽力下さり，日々の小さな相談にも親切に耳を傾けてくださった古宇田由夫技術職員に心から御礼申し上げます．
同じ研究グループの先輩として，著者が大学四年生時から懇切克明に御助言・御支
援頂きました西村正臣氏に万謝致します．研究室生活を送るにあたり，日頃から大変お世話になった神戸大学大学院自然科学研究科固体力学研究室の先輩，後輩諸氏に御礼申し上げます．そして三年間の研究室生活を充実と共に送ることができたのも，同輩である朝西誠君，久馬雅彦君，古賀裕君，田中克昌君の存在による処が大きく，感謝の念に堪えません．
最後に，6 年間に及ぶ大学・大学院生活は，家族の温情ある理解と惜しみない支援のお陰で，様々な物事に触れて多くの事を学ぶことができた貴重な時間となりました．
父 樋口宏と，母 容子をはじめとする家族の皆に，最大限の敬意と感謝の意を表してここに記します．ありがとうございました．
平成 20 年 2 月
樋口67

昌宏