修士論文

Bcc 鉄の界面力学特性評価 ：原子弾性剛性係数に基づく検討

指導教員：屋代如月

平成 26 年 2 月
神戸大学大学院工学研究科
機械工学専攻124T380T
山根 尭

Master Thesis
Mechanical properties of interface in bcc Fe :An examination based on atomic elastic stiﬀness

Takashi YAMANE
Feburary 2014
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan

要約
本論文では，表面薄膜やナノワイヤ，結晶粒界等の不均一構造を有した bcc 構造の鉄に対して，分子動力学法を用いて変形シミュレーション行い，不安定挙動と原子弾性剛性係数 Bijα の行列式 detBijα の関係について検討した．まず，単結晶バルク，表面薄膜，
ナノワイヤの引張シミュレーションを行った．単結晶バルクのひずみ制御による [001]引張では，応力ピーク後に局所変形による欠陥生成が生じす，なだらかな応力ひずみ曲線を示す．detBijα の平均値が応力ピーク時に負となることから，bcc 鉄では局所変形より前に系全体の不安定に達することが明らかとなった．つまり，応力制御で引張シミュレーションを行えば，detBijα が負となる点で応力急減が生じるものと考えられる．(110) 表面薄膜，(111) 薄膜および (111)-(1̄10) ナノワイヤも同様に，応力ひずみ応
答はなだらかな曲線を示し，応力ピーク時に detBijα の平均が負となる．これらの表面
ならびにエッジ部は転位生成のような局所変形ではなく，相変態などの集団的挙動による変形が支配的となることが明らかとなった．一方，(100) 表面薄膜や (100)-(010)，(110)-(1̄10) ナノワイヤでは，系全体の不安定よりも前に局所変形による応力急減が見られた．次に，Σ 値が 1〜99 の対称傾角粒界ならびに対称ねじれ粒界構造を作成した．
回転角度に対する粒界エネルギー変化と無負荷平衡状態における detBijα のバルクからの差や detBijα < 0 の数との対応を議論する．さらに，(111)Σ3 や (332)Σ11, (112)Σ3
傾角粒界の粒界面に垂直方向に引張シミュレーションを行った．(111)Σ3 粒界では欠陥原子が 1 原子層程粒界から拡張した後に粒界面からへき開した．また，引張前では

detBijα < 0 の原子が粒界近傍にのみ存在した．一方，(332)Σ11 粒界では 1 つ目の応力ピークから結晶構造の乱れの前縁が反対側に到達するまで引張は続いた．そして，粒界移動による構造変化が生じたことからひずみ硬化し最終的に粒界でへき開した．また，その前縁の原子の detBijα が負となることが明らかになった．(112)Σ3 粒界では負の

detBijα 原子が同様の挙動を示したが，Central symmetry parameter では双晶面を判別できなかった．最後に，(111)Σ3 や (332)Σ11, (112)Σ3 傾角粒界に対して粒界面に平行にせん断変形を行った．(332)Σ11 傾角粒界では先述したとおり，界面移動による粒界拡張が示された応力ピーク以降に粒界面が上下に移動し，欠陥部分が拡大するような変化を生じた．(111)Σ3 粒界は粒界を起点とした欠陥生成 (転位) を生じたが，(112)Σ3
粒界では粒内に欠陥が生じた．両者の傾角粒界では，負の detBijα 原子が前駆的に欠陥原子が示された範囲に生じた．

Summary
In order to discuss about the ralationships between unstable stress drop and determinant of atomic elastic stiﬀness Bijα , we performed molecular dynamics simulations
on thin plate, nanowire of bcc iron, as well as grain boundaries. First, we performed
tensile simulation on bulk, thin plate and nanowire. Bulk model shows smooth arch ofstress-strain curve against the [001] tension under the strain control, without generating defects by local deformation until later point after the peak. It is revealed that the
average of detBijα becomes negative at the stress strain peak so that we can conclude
that the bcc-Fe reaches global instability before the local deformation. That is, thestress-strain curve will show rapid elongation at the point if we implemented the simulation by stress control. (110), (111) surface thin plates, (111)-(1̄10) nanowire also
show smooth arch of stress-strain and the average of detBijα becomes negative at thepeak. Thus we can conclude that these surface and edge structure are not preferable
for local deformation such as dislocation nucleation so that the system would deform
rather by collective changes, e.g., phase transformation. On the other hand, (001) surface thin plate, (100)-(010) and (110)-(1̄10) nanowires show rapid stress drop by the
emergence of local deformation, before the global instability. Then we made laminate
structure of tilt and twist grain boundaries for various Sigma value ranging from 1 to
99. The change in the grain boundary energy against misorientation angle is discussedwith the diﬀerence of detBijα from bulk state and the number of negative detBijα atoms
under the no-load equilibrium. Furthermore, we performed tensile simulations on the(111)Σ3, (332)Σ11 and (112)Σ3 tilt grain boundaries, in the direction normal to them.
The (111)Σ3 shows cleavage cracking at the boundary, after showing the extension ofboundary disorder only 1 atomic layer just above and beneath the grain boundary. Before the extension, the atoms show negative detBijα on the extended adjacent layer. On
the other hand, the boundary disorder continues to extend after the first stress-strain
peak until the border of crystal/disorder travels to the opposite side in the laminateof (332)Σ11. The laminate shows strain hardening and finally cleavage cracking at theboundary, since the crystal structure is changed by the travel of the boundary. We canobserve negative detBijα atoms at the leading edge of the traveling border. The (112)Σ3
shows same tendency in the motion of negative detBijα , while the central symmetry
parameter cannot distinguish the twin boundary of (112)Σ3. Finally, we perform theshear simulation to the (111)Σ3, (332)Σ11 and (112)Σ3 tilt grain boundaries, in the
direction parallel to them. (332)Σ11 shows the extension of boundary disorder by thetravel of the border above mentioned. Defect nucleation or dislocation emission fromgrain boundary is observed in (111)Σ3, while disorder emerges in the middle of crystal
far from the boundary in the laminate of (112)Σ3. In the both cases, negative detBijα
atoms are observed in the defect area precursory.

目次
第1章 緒 論1

第 2 章 解析手法の基礎

62.1

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

62.2

原子間ポテンシャル. . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.3

高速化手法

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.4

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

11

2.5

原子弾性剛性係数 (Atomic Elastic Stiﬀness) . . . . . . . . . . . . . . .

12

2.5.1

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

12

2.5.2

原子弾性剛性係数 . . . . . . . . . . . . . . . . . . . . . . . . . .

15

第 3 章 薄膜・ナノワイヤの引張シミュレーションならびに AES に基づく検討

16

3.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.2

単結晶バルクの [001] 方向引張 . . . . . . . . . . . . . . . . . . . . . . .

18

3.3

表面薄膜モデルの引張 . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.3.1

(100) 表面薄膜モデル . . . . . . . . . . . . . . . . . . . . . . . .

19

3.3.2

(110) 表面薄膜モデル . . . . . . . . . . . . . . . . . . . . . . . .

22

3.3.3

(111) 表面薄膜モデル . . . . . . . . . . . . . . . . . . . . . . . .

25

細線モデルの引張 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.4.1

(100)-(010) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

27

3.4.2

(110)-(110) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

30

3.4.3

(111)-(110) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

33

3.4

第 4 章 対称粒界の粒界エネルギーと AES，ならびに引張シミュレーション

4.1

無負荷平衡状態の粒界エネルギーと AES . . . . . . . . . . . . . . . . .

i

36
36

4.2

目次

ii

4.1.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.1.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

対称傾角粒界の引張シミュレーション. . . . . . . . . . . . . . . . . .

41

4.2.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.2.2

(332)Σ11 傾角粒界 . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.2.3

(111)Σ3 傾角粒界 . . . . . . . . . . . . . . . . . . . . . . . . . .

47

4.2.4

(112)Σ3 傾角粒界 . . . . . . . . . . . . . . . . . . . . . . . . . .

49

第 5 章 対称傾角粒界のせん断シミュレーションならびに AES 解析

51

5.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.2.1

(332)Σ11 傾角粒界 . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.2.2

(111)Σ3 傾角粒界 . . . . . . . . . . . . . . . . . . . . . . . . . .

57

5.2.3

(112)Σ3 傾角粒界 . . . . . . . . . . . . . . . . . . . . . . . . . .

61

第 6 章 結論

64

参考文献

66

学術論文・学術講演

69

謝辞

76

第1章
緒論
金属工学は，多種な金属・合金材料の特性を生かし，新しい機能・特性を持った材料の開発を進めてきた．同時に，加工技術の飛躍的な発展に伴い，合金化による特性の向上という観点だけではなく，組織制御のスケールをナノレベルにまで小さくする
ことで新しい特性を発現することが期待されている．ナノ材料の特性を左右する要因の一つとして「表面」が挙げられ，その割合はマクロな材料に比べて極めて大きくなることから，通常の寸法では顕在化しない表面特有の性質が現れることが考えられる．
ナノスケールの材料の例として，ナノピラーと呼ばれる直径ならびに高さ数十〜数百

nm の柱状のナノ構造体がある (図 1.1)[1]–[3] ．ナノピラーは，ナノインプリントと呼ばれる基板上に塗布した樹脂材料に押し付けて形状を転写する微細加工技術を用いて作製され，免疫分析チップとして DNA 解析等に利用されている．ナノピラーは微小な面積の基板上に大きな表面積を作ることができるため，表面での反応を促進させ，また，微細チャンネルを生かした成分分離の機能を有する．ナノピラー間の距離や表面積を変化させることで，表面に吸着させる抗体の密度も制御が可能となる [4] ．
また，「結晶粒界」は表面と同様，二次元格子欠陥の一種であり，多結晶材料の力学特性を大きく左右する重要な要因である．一般的に，Hall-Petch 則により，多結晶体を構成する結晶粒の大きさを微細化することで，強度の向上を示すことが知られている．しかし，これはあくまでマクロスケール (数 µm 以上) で適用されるものであり結晶粒の寸法がナノサイズになると，粒界が占める割合が上昇し，従来材料にない性質を示すようになる．実際に，ECAP[6]–[10] 等の強ひずみ加工により結晶粒の寸法をナノ

1

2μm

Fig.1.1 Image of 36 x 36 silicon nanopillar array.[5]メートルとした材料では，強度が飛躍的に上昇すること，強度と靱性が共に向上する等の現象が報告されている [11] ．
ナノスケールにおける表面の特性を明らかにすることを目的として，ナノ薄膜及びナノワイヤの分子動力学法シミュレーション研究が広く行われている．薄膜について
は，Y.Gan らは，Cu のナノ薄膜の引張シミュレーションを行い，小さいほど柔らか
く，薄膜強度が高いことが報告している [12] ．Y.Tong らは，Bi2 Te3 ナノ薄膜の引張シ
ミュレーションを行い，緩和後にバルクでは圧縮応力，表面では引張応力が働いていること，厚さ及び温度によるヤング率の変化を報告している [13] ．ナノワイヤ解析では，H.A.Wu らは，Cu のナノワイヤの単軸引張シミュレーションを行い，温度効果，
寸法効果による弾性係数変化及び応力分布変化を報告している [14] ．Y.H. Wen らは，

Ni のナノワイヤの単軸引張シミュレーションを行い，ひずみ速度の大きさによって構造変化に大きな差が生じることを報告している [15] ．Y.Jing らはシリコンの単軸引張シ
ミュレーションを行い，限界荷重が系の直径やひずみ速度，温度に大きく影響していること，融解温度が直径の減少に伴って低下することを報告している [16] ．また，表面力学の分野において，表面欠陥がもたらす影響を解析する研究も数多く行われており，

2

M.Grujicic らは，Fe-Ni 合金中のマルテンサイト変態について MD シミュレーションを用いて表面の持つ結晶欠陥がマルテンサイトの構成に大きく影響していることを報告している [17] ．
粒界構造がナノ材料に及ぼす影響についても，粒界エネルギーと原子構造の関係性や，第一原理計算による電子状態の解析など，原子種に依らず，広く研究が行われている．J.Chou らは Mg[11̄00] 対称傾角粒界について [18] ，中島らは α-Fe⟨110⟩ 対称傾角
粒界について [19] ，池田らは，Mo⟨001⟩ 対称傾角粒界と Cu⟨110⟩ 対称傾角粒界について[20]，それぞれ分子動力学法を用いて粒界エネルギーとその原子構造の依存性について明らかにし，基本ユニットとした構造ユニットモデルから対称傾角粒界を構成可能であることを報告している．水口らは，Pd について ⟨110⟩ 対称傾角粒界の粒界エネルギーと原子構造を，分子動力学法を用いて系統的に明らかにし，Hartree-Fock-Slater 法を用いて，完全結晶と粒界に水素原子を配置したときの電子状態を計算することで，粒界の電子状態が粒界構造に依存することを明らかにしている．[21] ．R.Wang らは，Alと Cu の傾角，ねじれ粒界における，界面結合のエネルギーを第一原理計算を用いて解析し，原子による積層欠陥エネルギーの差や，界面結合の挙動の違いを明らかにしている [22] ．北村らは，Al について Σ5 傾角粒界と ⟨100⟩ 表面の会合部の原子構造及び形成エネルギーを分子動力学法を用いて明らかにし，会合部の原子構造に乱れが生じて，表面に凹凸が形成されることを報告している [23] ．また，Douglas らは，Cu を対象として，⟨100⟩ 軸，⟨110⟩ 軸を回転軸とした様々な Σ 値の傾角粒界について分子動力学法を用いて引張強度の比較を行い，転位の生成過程を明らかにしている [24] ．また，

Tschopp らは，Fe を対象として，分子動力学法を用いて，空孔形成エネルギーと自己格子間原子エネルギーについて，粒界からの距離，Σ 値，回転角度，粒界エネルギーとの依存性を粒界構造を交えながら明らかにしている [25] ．また，S.P.Chen らは，Ni3 Alの ⟨001⟩ 対称傾角粒界を対象として，分子動力学法を用いて粒界エネルギーを計算し，純金属と同様の粒界構造を有するがゆがみや多様性を持つこと，Ni の粒界エネルギー
ならびに結合エネルギーとほぼ同等の値を持つことを明らかにしている [26] ．
対称粒界だけではなく，非対称の粒界に対しても同様の検討が数多く行われている．

Tschopp らは Al と Cu の Σ3 非対称傾角粒界 [27] ，Cu⟨110⟩ 非対称傾角粒界について [28] ，3

分子動力学法を用いて粒界からの転位生成を議論している．また，Cu と Al の低 Σ 値の粒界に対して，その構造と粒界エネルギーも報告している [29] ．また，J.A.Brown らは，CuΣ11⟨110⟩ 非対称傾角粒界について，分子動力学法を用いて粒界エネルギーとその構造を計算し，高い粒界エネルギーを持つ粒界は，ナノサイズの範囲で，異なる Σ
値，かつ CSL にも属さない粒界から構成されること等を明らかにしている [30] ．また，
シミュレーションの分野だけではなく，実験の観点から石田らは粒界すべりの機構や粒界孔の生成機構などを明らかにしている [31] ．
一方，著者らのグループでは，連続体近似が成立しないナノスケールでの変形・破壊を統一的に議論する基準として，個々の原子位置における原子弾性剛性係数 (AES)の正値性に基づいた検討を行ってきた [32] ．AES による検討は，fcc，アモルファス金属，そしてダイヤモンド構造のシリコンについて広く進めてきた [32]–[34] ．fcc 構造の Niについての検討では，転位の発生が AES の正値性によって評価可能であることを報告している [32] ．アモルファス構造 (Ni-Al) への検討では，AES の正値性によって局所クラスタの崩壊及び再構成が表されることを報告している [33] ．また，シリコンの検討では，表面 (内部不均一を有する系) の解析で，弾性限界まではひずみ増加による変形状
態が系全体で均一化することを報告している [34] ．
本研究では，bcc 金属について，Fe をモデル材料として初期不均一構造をもつ表面
薄膜，ナノワイヤ，対称粒界を導入した系に対して引張やせん断変形を与え，非弾性
変形開始時の AES について明らかにすることを目的とする．以下に，各章の概略を示
す．第 2 章では，MD 法の基礎理論を述べ，使用する原子間ポテンシャル，不安定性を判別する際に用いる原子弾性剛性係数について説明する．第 3 章では，表面薄膜ならびにナノワイヤの引張シミュレーションを行い，その変形挙動を AES 解析を用いて検討を行う．低指数面である (100)，(110), (111) の表面について薄膜・ナノワイヤの傾向について不安定挙動の比較を行う．第 4 章では，対称ねじれ粒界ならびに対称傾角粒界の粒界エネルギーと AES の角度依存性について検討を行うとともに，引張シミュ
レーションを行い，AES と不安定挙動 (=応力急減) の関係性について議論する．第 5
章では，第 4 章で導入した対称傾角粒界に対して，異なる変形モードであるせん断変
形下の不安定挙動を観察し，AES との関係性について検討を行う．最後に，第 6 章で

4

結論を述べる．

5

第2章
解析手法の基礎2.1

分子動力学法

分子動力学法 (molecular dynamics method，略して MD 法) は，系を構成する各粒子についてニュートンの運動方程式

mid2 r i
= Fi
dt2

(2.1)

をたて，これを数値積分することにより粒子の軌跡を求める方法である．ここで，mi ，

r i はそれぞれ粒子 i の質量および位置ベクトルである．粒子 i に作用する力 F i は，系のポテンシャルエネルギー Etot の各位置における空間勾配として次式により求められる．

Fi = −

∂Etot
∂r i

(2.2)

式 (2.1) の数値積分には，Verlet の方法，予測子–修正子法等がよく用いられる．本
研究では，以下に示す Verlet の方法を用いた．
時刻 t + ∆t と t − ∆t での粒子 i の位置ベクトル r i (t ± ∆t) を Taylor 展開すると，()
dr i (t) (∆t)2 d2 r i (t)
3
+
+
O
(∆t)
dt
2
dt2
2 2
(
)
dr i (t) (∆t) d r i (t)3r i (t − ∆t) = r i (t) − ∆t
+
+
O
(∆t)
dt
2
dt2

r i (t + ∆t) = r i (t) + ∆t

(2.3)
(2.4)

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

dr i= v i (t)
dt
6

(2.5)

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

(2.6)
(2.7)

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

r i (t + ∆t) + r i (t − ∆t) = 2r i (t) + (∆t)2()
F i (t)+ O (∆t)4mi(

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

)

(2.8)
(2.9)

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

r i (t + ∆t) = 2r i (t) − r i (t − ∆t) + (∆t)2v i (t) =()
F i (t)+ O (∆t)4
mi()
1{r i (t + ∆t) − r i (t − ∆t)} + O (∆t)22∆t(2.10)(2.11)と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．初期の計算 (t = 0) では，t = ∆t での座標 r i (∆t) は式 (2.6) と初速度から得ることができる．

2.2

原子間ポテンシャル

式 (2.2) で示したように，原子 α に作用する力 F α は系のエネルギー Etot をポテン
シャルとして決定される．したがって，系のポテンシャルエネルギー Etot をいかに精度よく評価するかが重要となる．量子力学に基づき，電子や原子核のハミルトニアン
から系のポテンシャルエネルギーを精密に求めて原子の運動を追跡する第一原理分子
動力学法も試みられているが，計算量が極めて膨大になるため，ごく少数の原子しか扱うことができず，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する原子間ポテンシャルが通常用いられる．

7

2.2.0.1

FS ポテンシャル

本研究では，金属系の MD シミュレーションで広く用いられている EAM(原子埋め込み法) ポテンシャルの一つである FS(Finnis Sinclair) ポテンシャルを使用した [35] ．FS
ポテンシャルは，鉄のシミュレーションによく用いられている．FS ポテンシャルでは系全体のエネルギーを次式で表す．
N
∑

N ∑
N
1∑
Etot = −AF (ρ ) +V (rαβ )2 α β(̸=α)αα(2.12)

ここで，ρα は原子 α の位置における多体効果を考慮する密度を表し，F (ρα ) は密度 ραの位置に原子を埋め込むエネルギー，V (rαβ ) は距離 rαβ 離れた原子 α と β のクーロン相互作用であり，A はポテンシャルパラメータである．密度 ρα は周囲の原子 β からの寄与 ϕ(rαβ ) の重ね合わせで与えられると仮定しρα =neighbor∑ϕ(rαβ )

(2.13)

β(̸=α)

で評価する．ϕ(rαβ )，V (rαβ ) の具体系は次式である．
()

V rαβ =()

ϕ rαβ =

)2 (
)
 (
αβ
αβ
αβ 2

r
−
c
c
+
c
r
+
c
(r
)
(c < rαβ )

0
1
2


(c ≫ rαβ )

0



)2
(
)3
 (
αβ
αβ

r
−
d
+
β
r
−
d
/d



0




(d < rαβ )(d ≪ rαβ )(2.14)(2.15)

c, c0 , c1 , c2 , d, β は，FS ポテンシャルのフィッティングパラメータであり，表 2.1 に bcc
鉄の値を示す．

Table 2.1 Potential parameters for FS potential.
A[eV]1.828905c0
c11.2371147 −0.3592185

c2−0.0385607

8

c [nm]
0.340

d [nm]
0.3569745

β
1.8

2.3

高速化手法

原子数 N の系において粒子間の全相互作用を評価すると，1step に N × (N − 1) 回の計算が必要となり，N が大きくなると極めて膨大な計算量となる．実際には，一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ半径 rc )内の粒子からの寄与を効率よく計算することにより高速化できる．従来よく用いられてきた高速化手法に粒子登録法がある．これは，図 2.1 に示したように，rc よりひとま
わり大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) に計算負荷が減少される．しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計算負荷により打ち消される．

rc

rfc

Fig.2.1 Schematic of bookkeeping method.

9別の高速化手法としてブロック分割法がある．図 2.2 に示すように，シミュレートする系をカットオフ距離程度の格子状に分割し，各ブロックに属する粒子をメモリーに記憶する．着目している粒子に作用する力を評価する際には，その粒子が属するブロックおよび隣接するブロックから相互作用する粒子を探索して行う．粒子が属する
ブロックは，粒子の位置座標をブロックの辺長 bx，by で除した際の整数により判断できるので，ブロック登録時の計算負荷は粒子数 N のオーダーとなる．したがって，粒子登録法では登録更新の負荷が大きくなるような大規模な系でも高速化が可能である．

y

by
0

x

bx

Fig.2.2 Schematic of domain decomposition method.

10

2.4

速度スケーリング法

分子動力学解析における温度制御には一般的には速度スケーリング法が用いられる．
この方法は，統計熱力学より導かれる式 (2.16) を用いて，以下のように制御する．

1 α α α 3m vi vi = k B T
22(2.16)

mα : 粒子αの質量
viα : 温度 T での粒子αの速度kB : Boltzmann 定数 = 1.38 × 10−23 [J/K]
目標の温度 T0 における原子 α の速度を viα0 とおくと viα0 は式 (2.17) のように表される．

(

viα0 =

3kB T0mα)0.5(2.17)

同様に，温度 T の時の原子 α の速度は式 (2.18) のように表される．
(

viα =

3kB T
mα

)0.5

(2.18)

よって，式 (2.17) と式 (2.18) より以下の式が得られる．
(

viα0
T0
=
α
vi
T

)0.5

(2.19)

つまり，系の温度を T から T0 にするには，式 (2.19) の右辺を現在の速度に掛けてやればよい．ただ，これだけでは原子配置に反映されないので，Verlet 法における ∆r αi (t+∆t)√（式 2.20）を

T0 /T ∆r αi (t + ∆t) と置き換える必要がある．

∆r αi (t + ∆t) = r αi (t + ∆t) − r αi (t) = r αi (t) − r αi (t − ∆t) + (∆t)2

F αi (t)(2.20)mα

平衡状態では，能勢の方法 [36] など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されて
いる．

11

2.5

原子弾性剛性係数 (Atomic Elastic Stiﬀness)

2.5.1

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

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

1
σij =
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

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

Bijkl ≡

∂σij
∂εkl

)

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

ここで，δij はクロネッカーのデルタである．Bijkl は非線形弾性における応力-ひずみ関係の勾配を表すことから，Wang, Yip らは，ひずみの対称性を考慮したテンソル
sym
Bijkl
≡ (Bijkl + Blkji )/2 の正値性によって結晶の安定性を評価することを提案している．Bijkl の最小固有値が負であれば，対応するひずみに対して応力が負になる自発的変形モードが存在することになり，不安定と判別される．

2.5.1.1

原子系への拡張ααの算出に必要な原子応力 σij
局所の安定性を評価するための原子弾性剛性係数 Bijklαは，各原子周りの微小ひずみに対するポテンシャルエネル
ならびに原子弾性係数 Cijkl

ギーの 1 次，2 次変化量として導出される．

原子応力
簡単のため，結晶の内部エネルギー U が系全体のポテンシャルエネルギー Etot に等しいとする．このとき，応力は平衡状態からの微小ひずみ η に対するポテンシャル

12エネルギーの単位体積当たりの変化として与えられる [37] ．

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

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

ηij =

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

1dηij = dGij2

(2.27)と式 (2.25) の関係から次の関係が得られる．

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

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

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

) αβ αβN ∑N (1 ∑1 ′ αβ ri rj′ α ′ αβ=F (ρ )ρ (r ) + ϕ (r )V α β(̸=α)
2rαβ

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

σijα =) αβ αβN (1 ∑1 ′ αβ ri rj′ α ′ αβF (ρ )ρ (r ) + ϕ (r )V /N β(̸=α)
2
rαβ

13

(2.29)

と定義すると，系の応力は
N
1 ∑
σij =
σijα
N α

(2.30)となる．

原子弾性係数弾性係数も応力と同様に U ≈ Etot の場合にはCijkl =

1 ∂ 2 Etot
V ∂ηij ∂ηkl

(2.31)であるので，平衡状態からの仮想均一変形を考えると FS ポテンシャルにおける弾性係数は以下のようになる．

Cijkl ==(N ∑αγ αγ )N ∑N
N∂rαβ ∂ ∑1 ′ αγ ri rj
1 ∑′ α ′ αγ
{F
(ρ
)ρ
(r
)
+
ϕ
(r
)}
2V α β(̸=α) ∂ηkl ∂rαβ α γ(̸=α)
2
rαγ
(
[N ∑
) αβ αβ αβ αβNρ′ (rαβ ) ri rj rk rl
1 ∑F ′ (ρα ) ρ′′ (rαβ ) −V α β(̸=α)
rαβ(rαβ )2αγ αγ }
Nriαβ rjαβ }{ ∑′ αγ rk rl
ρ (r ) αγ+
ρ (r ) αβF (ρ )rr
αγ(̸=α)
β(̸=α)
N
∑

′′

α

{ ∑
N

′

αβ

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

応力と同様に，各原子位置における原子弾性係数を以下のように定義する．
α
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.33)

これより，系の弾性係数はCijkl =N
1 ∑CαN α ijklのように原子弾性係数の平均となる．

14(2.34)

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

15

第3章
薄膜・ナノワイヤの引張シミュレー
ションならびに AES に基づく検討
3.1

解析条件

図 3.1(a)，(b) に模式的に示すように，立方体セルの y, z 方向に周期境界を適用した薄膜モデル，ならびに z 方向にのみ周期境界を適用した細線モデルを解析対象とす
る．立方体セルの原子構造は，図 3.1(i)〜(iii) に示すように x，y ，z の結晶方位を回転させてほぼ同じ大きさに切り出した系とする．各セルの総原子数は (a)16000，(b)15680，(c)15456 である. また，完全結晶との比較を行う為，図 3.1(i) の解析モデルに全方向周期境界を適用したバルクのシミュレーションも行う．
原子間相互作用の評価には FS ポテンシャル [35] を用いた．まず，温度 0.1[K] で構造緩和計算を 10000[fs] 行い，その後 z 方向に毎ステップひずみ増分 ∆εzz =1.0×10−6 を与えて引張りひずみを増加させた．この時，∆εzz =1.0×10−4 毎に原子座標データを記録
α [32]
αならびに原子弾性係数 Cijkl
を求め，AES の行列式
し，各原子位置の原子応力 σij

detBijα (Voigt 表記) を評価した．

16

5.73nm

nm

5.73nm

[001]
[010]
[100]

z

(i) (100) surface(100)-(010) nanowire(a) Thin plate model

5.73nm

y
x

5.73

[001]
[110]
[110]

5.68nm5.68

nm

5.62nm(ii) (110) surface(110)-(110) nanowire

[112]

z

[110]
[111]

y
x

5.71nmnm(iii) (111) surface(111)-(110) nanowire(b) Nanowire model
Fig.3.1 Simulation models．

17

5.68

3.2

単結晶バルクの [001] 方向引張

図 3.2 に全方向周期境界条件下のバルクに引張りひずみを与えたシミュレーション
結果を示す．図中黒の実線は応力ひずみ関係であり，その目盛を右軸に，赤い実線はdetBijα の平均値であり，標準偏差の範囲をエラーバーで示している．detBijα の値は 0[K]完全結晶での値で無次元化しており，目盛は左に付している．下の図は detBijα が負の原子数の変化を示す．応力ひずみ応答を見ると ϵzz = 0.0982 で応力ピークを示すが，この時点では原子配置に大きな変化を生じないままわずかに低下し，その後 ϵzz = 0.1252
で結晶構造に乱れを発生して応力が急減した. 応力ピークの ϵzz = 0.0982 では全原子の AES が負値を示し，系全体の不安定となっていることが分かるが，ひずみ制御の引張では先述のように局所変形は生じず応力は急減しない. 応力制御で引張りを行っていればこの点で系は応力上昇できなくなるため不安定となると考えられる．すなわち，

bcc 鉄の [001] 引張では，内部の局所変形発生よりも全体的な不安定 (dilatancy) が先に
15

1.5

Stress strain curve
ε zz=0.1252

ε zz=0.0982

1

10

5
0.5
0
0

Stress, σzz, GPa

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

生じる．

-5

Average of detBαij

Number of detBαij <0 atoms

-0.5
0

0.05

0

0.05

0.1

-10
0.15

0.1

0.15

15000
10000
5000
0

Strain, εzz

Fig.3.2 Change in the average，standard deviation of detBijα ，and stressstrain curve of bulk model.

18

3.3

表面薄膜モデルの引張3.3.1(100) 表面薄膜モデル

図 3.3 は (100) 薄膜モデルの応力ひずみ関係，detBijα の平均値の変化及び負の AES
原子数変化である. また，前節のバルクでは，detBijα の標準偏差を示していたが，表面と結晶内部の差異を明確にする為，表面原子層のみで平均化した detBijα の値と，AESが負の表面原子数の変化を橙色で示した．図 (b) は応力急減近傍を拡大したグラフである. ϵzz = 0.0947 で応力ピークはあるが，この点で，表面原子層の detBijα の平均値が0 になっている．先のバルクの結果と異なり，系全体の平均は負に達しておらず系全体は不安定となっていない．すなわち「局所不安定」が「全体不安定」より先に現れたケースである．右の拡大図を見ると，応力が急減する点 (ϵzz = 0.0965) 以降，detBijαの平均値 (表面原子層) が回復している．右下の負の AES 原子数変化に着目すると，応力急減する点までの負の原子は全て表面原子層のものであるが，この点以降，表面原子層の負の AES 原子数が減少している一方で，全体の負の AES 原子は急増しており，内部に負の原子が増加していることが分かる．応力ピーク点前後の原子配置を図 3.4に，負の AES 原子を赤く着色して示した．応力ピーク点では表面原子の detBijα の平均が負になっているものの，表面原子がすべて負にはなっていない．応力急減直前のϵzz = 0.0955 では表面原子が全て負となり，その後結晶構造にずれを生じて結晶内部に負の AES 原子が現れる．

19

10

0.5

5

0

0
Average of detBijα(All Atoms)

Average of detBijα(Surface Atoms)

-1
6000

Number of detBijα <0 atoms

All Atoms

2000
Surface Atoms
0
0

0.05

0.1

0.15

0.2

Strain, εzz

Stress strain curve
10

0.5

5

0

0

-0.5

-5

-10

4000

1

15

εzz=0.0947
εzz=0.0965

Average of detBijα
Average of detBijα (Surface Atoms)(All Atoms)

-1
6000

Stress, σzz, GPa

Average of 6x6 determinant, ΣdetBijα /N,normalized by detBijα at perfect lattice

εzz=0.0947
1-0.5

Number of detBijα <0 atoms

1.5

15
Stress strain curve

Stress, σzz, GPa

Average of 6x6 determinant, ΣdetBijα /N,normalized by detBijα at perfect lattice

1.5

-5

-10

All Atoms

4000
2000

Surface Atoms
0
0.09

0.095

0.1

0.105

0.11

0.115

0.12

Strain, εzz(a) Outline(b) Zoom-up around stress drop

Fig.3.3 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (100) surface.

20

z
y

x

ε zz =0.0965

ε zz =0.0920

ε zz =0.0947(peak)ε zz =0.0970

ε zz =0.0980

ε zz =0.0955

Fig.3.4 Snapshots of (100) surface thin plate. Red coloured circles indicate
detBijα < 0 atoms.

21

3.3.2(110) 表面薄膜モデル

図 3.5 は (110) 薄膜モデルのシミュレーション結果を示す．先の (100) 薄膜モデルと同じ引張方向 ([001])，同程度の原子数であるにも関わらず，(100) 薄膜とは異なりバルクと似た傾向を示した．ϵzz = 0.1024 で応力ピークを迎えるが，すぐに応力急減することなくなだらかに低下し，ϵzz = 0.1130 で応力急減した．detBijα の変化を見ると，応力ピーク時に系全体の平均値は 0 となり，表面原子層の平均値は応力ピークより前に負値となっている．下の負の AES 原子数の変化を見ると，応力ピーク前に表面の負のAES 原子が 1105 個発生し，ピーク直後にすべての原子の AES が負となっている．応力ピーク付近の原子配置を図 3.6 に示す．ひずみ ϵzz = 0.0913 近傍で表面原子がほぼす
べて負になり，応力ピークの次のステップで全原子の AES が負となっている．その後，表面原子から構造緩和が生じ，表面原子の detBijα が正値となっているが，バルクの時と同様，応力制御の引張であれば，ピーク点以降はひずみ延伸に対して系が応力上昇できず，系がさらに伸ばされて不安定となる．表面が存在するにも関わらず系の局所不安定よりも系全体の不安定が先に現れたことは，[001] 方向の引張に対して (110) 表
面からの局所変形の発生しにくさを示唆している．

22

10

0.5

5

0

0
Average of detBijα(Surface Atoms)-0.5

Stress strain curve -5αAverage of detBij(All Atoms)

Number of detBijα <0 atoms

All Atoms
10000
Surface Atoms
5000
0
0

0.05

0.1

0.15

0.2

Strain, εzz

εzz=0.1130

εzz=0.1024

1

10

0.5

5

0

0

-0.5

Average of detBijα(Surface Atoms)

Average of detBijα(All Atoms)

-1

-10

15000

15
Stress strain curve

15000

Stress, σzz, GPa

εzz=0.1130

Average of 6x6 determinant, ΣdetBijα /N,normalized by detBijα at perfect lattice

εzz=0.1024

1-1
Number of detBijα <0 atoms

1.5

15

Stress, σzz, GPa

Average of 6x6 determinant, ΣdetBijα /N,normalized by detBijα at perfect lattice

1.5

-5

-10
All Atoms

10000

Surface Atoms

5000
0
0.09

0.095

0.1

0.105

0.11

0.115

0.12

Strain, εzz(a) Outline(b) Zoom-up around stress drop

Fig.3.5 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (110) surface.

23z

ε zz =0.0913

ε zz =0.1032yx

ε zz =0.1024(peak)ε zz =0.1130

ε zz =0.1031

ε zz =0.1140

Fig.3.6 Snapshots of (110) surface thin plate. Red coloured circles indicate
detBijα < 0 atoms.

24

(111) 表面薄膜モデル

3.3.3

図 3.7 は (111) 薄膜モデルのシミュレーション結果である．本モデルでも応力ピーク
後になだらかに応力低下した後，局所変形による応力急減を生じているため，応力制御のシミュレーションであればピーク点で系全体の不安定を生じるものと考えられる．

(111) 表面原子の detBijα は引張前はバルクのそれより低いが，引張により大きく低下
するのは内部の detBijα であり，応力ピークの直後は表面と内部の大小関係が逆転する．
応力ピーク付近の原子配置を図 3.8 に示す．応力ピーク時 (ϵzz = 0.0772) には表面から3 原子列目のみに負の AES 原子が生じ，応力急減以前のひずみ ϵzz = 0.0891 で，表面1 原子列を除いて全ての原子が負となる．その後，結晶内部で構造にずれが生じ，バ
ルク部分の不安定原子が回復し，表面原子列に負の AES 原子が伝播している．

10

0.5

5

0

0
Average of detBijα(All Atoms)-0.5

Average of detBijα(Surface Atoms)-5

Number of detBijα <0 atoms

All Atoms

10000
Surface Atoms
5000
0
0

0.05

0.1

0.15

0.2

Strain, εzz

εzz=0.0772
εzz=0.0925

1

10

Stress strain curve
0.5

5

0

0

-0.5 Average of detBijα(All Atoms)

Average of detBijα(Surface Atoms)

-1

-10

15000

15

-5

-10

15000
10000

All Atoms

5000

Surface Atoms

0
0.07

0.075

0.08

0.085

0.09

0.095

0.1

Strain, εzz(a) Outline(b) Zoom-up around stress drop

Fig.3.7 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (111) surface.

25

Stress, σzz, GPa

Stress strain curve
1-1
Number of detBijα <0 atoms

1.5

15

εzz=0.0925

Average of 6x6 determinant, ΣdetBijα /N,normalized by detBijα at perfect lattice

εzz=0.0772

Stress, σzz, GPa

Average of 6x6 determinant, ΣdetBijα /N,normalized by detBijα at perfect lattice

1.5

z
y

ε zz =0.0820

ε zz =0.0772(peak)x

ε zz =0.0891

ε zz =0.0793

ε zz =0.0930

ε zz =0.0810

Fig.3.8 Snapshots of (111) surface thin plate. Red coloured circles indicate
detBijα < 0 atoms.

26

3.4

細線モデルの引張

3.4.1(100)-(010) 細線モデル

図 3.9 は (100)-(010) 細線モデルのシミュレーション結果である．前節では，表面原子層で detBijα を平均化して評価したが，本節では表面と表面の会合部 (エッジ部) で平均した detBijα の変化及び負の AES 原子数の変化を橙色で示した．また，応力急減近傍の AES 変化を図 3.10 に示す．本モデルでは，ϵzz = 0.0822 において ⟨110⟩ 方向のエッジ部が負となっている．ここで，⟨1̄10⟩ 方向が負となっていないのは，図 3.11 に模式的に示すようにバルク bcc から切り出した時にエッジ部分に原子列の有無による違いがあるためである．前節の (100) 表面薄膜モデルと同様，系全体の detBijα の値が正値の時に局所変形発生による応力急減を生じており，系全体より局所不安定が先に現れるケースである．エッジ部の detBijα は正値ではあるがほぼ 0 となっている．図 3.9(b)の右図をみると，応力ピークのわずか前にエッジ部に負の AES 原子が生じ，応力ピー
ク時を境に減少して結晶内部で増加している．

27

10

εzz =0.0871

Stress strain curve

1.5

6
1
4
0.5

2

0

0

Average of detBij(All Atoms)-0.5

Number of detB ijα <0 atoms

8

Stress, σzz , GPa

Average of 6x6 determinant, ΣdetB ijα /N ,normalized by detB ijα at perfect lattice

2

αAverage of det B ij
(Edge Atoms)
α

-2
-4

-1
5000
4000

All Atoms

3000
2000

Edge Atoms

1000
0
0

0.05

0.1

0.15

0.2

Strain, εzz(a) Outline
10

Number of detB ijα <0 atoms

0.4

6
1
Stress strain curve 4
0.5

2

0

0
Average of det Bij(Edge Atoms)

α6Stress strain curve 4

0.2

2
0

0

-2

-2

Average of det Bij(Edge Atoms)

α

-0.2

-4

-1
5000

-4

400

4000

300

All Atoms

3000
2000

8

All Atoms

200
100

Edge Atoms

Edge Atoms

1000
0
0.07

0
0.075

0.08

0.085

0.09

0.095

0.07

0.1

0.075

0.08

0.085

0.09

0.095

0.1

Strain, εzz

Strain, εzz(b) Zoom-up around stress drop
Fig.3.9 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (100)-(010) nanowire.

28

Stress, σzz , GPa

8

1.5-0.5 Average of det Bαij
(All Atoms)

10εzz =0.0871

εzz =0.0871

Stress, σzz , GPa

Average of 6x6 determinant, ΣdetB ijα /N ,normalized by detB ijα at perfect lattice

2

Average of det Bαij
(All Atoms)

z
y

x

ε zz =0.0785

ε zz =0.0871(peak)ε zz =0.0822

ε zz =0.0877

ε zz =0.0870

ε zz =0.0900

Fig.3.10 Snapshots of the (100)-(010) surface nanowire. Red coloured circles
indicate detBijα < 0 atoms.

edge

edge

y

z
y

x

x
(a) Bcc latticecut(b) Top viewFig.3.11 Explanation of the diﬀerence of edge structure in (100)-(010) nanowire.

29

3.4.2(110)-(110) 細線モデル

図 3.12 に，(110)-(110) 細線モデルの引張シミュレーション結果を示す. 薄膜の場合と異なり，(110) 表面でも細線では系全体の不安定に達する前に局所変形発生による応力急減を生じている．一方，負の AES 原子の変化を見ると，応力ピーク前にエッジ部の detBijα が負になるのは先の (100)-(010) 細線と同じであるが，応力急減より前にエッジ部以外にも負の原子が増加している．
図 3.13 の原子配置を見ると，エッジ部に負の AES 原子が生じ (ϵzz = 0.0901)，さらに応力ピーク直前 (ϵzz = 0.0993) では (110) 及び (110) 表面全体が負となっている. その後構造変化により負の AES 原子は結晶内部に現れ，表面やエッジ原子の AES は正値に回復する．

30

Stress strain curve
8

1.5

6
1
4
0.5

2

0

0

-0.5

Number of detB ijα <0 atoms

10εzz =0.0998

Stress, σzz , GPa

Average of 6x6 determinant, ΣdetB ijα /N ,normalized by detB ijα at perfect lattice

2-2

Average of det Bij(Edge Atoms)

αAverage of det Bij(All Atoms)

α

-4

-1
All Atoms

3000
2000
Edge Atoms

1000
0
0

0.05

0.1

0.15

0.2

Strain, εzz(a) Outline
Average of det Bij(All Atoms)

αStress strain curveαAverage of det Bij(All Atoms)

εzz =0.0998

2

0

8

εzz =0.0998

6
4

0.5

0

4
2
0

-2

Average of det Bij(Edge Atoms)

α

6

0.1

0
-2

Average of det Bij(Edge Atoms)

α

-4

-1

-4

-0.1
400

3000

300

All Atoms

2000

All Atoms

200

Edge Atoms

1000
0
0.08

10

Edge Atoms

100
0

0.085

0.09

0.095

0.1

0.105

0.11

0.08

0.085

Strain, εzz

0.09

0.095

0.1

0.105

0.11

Strain, εzz(b) Zoom-up around stress drop
Fig.3.12 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (110)-(110) nanowire.

31

Stress, σzz , GPa

1

0.2

8

1.5

-0.5

Number of detB ijα <0 atoms

10

Stress, σzz , GPa

Average of 6x6 determinant, ΣdetB ijα /N ,normalized by detB ijα at perfect lattice2

Stress strain curve

z
y

x

ε zz =0.0993

ε zz =0.0901

ε zz =0.0998(peak)ε zz =0.0929

ε zz =0.1015

ε zz =0.0950

Fig.3.13 Snapshot of the (110)-(110) surface nanowire. Red coloured circles
indicated detBijα < 0 atoms.

32

3.4.3(111)-(110) 細線モデル

図 3.14 および図 3.15 は (111)-(110) 細線モデルのシミュレーション結果である. これまでの応力ひずみ曲線に比べ，のこ歯状の応答を示しているが，これは引張方向の違いによる応力波の影響と考えられる. ϵzz = 0.0826 で応力はすでに減少しはじめているが系全体の平均 detBijα は減少し続け，ϵzz = 0.0857 で負値となる．図 3.14(b) 右の拡大図で見ると，系全体の detBijα が回復し始めるのは ϵzz = 0.0950 で，このひずみの後に応力減少が急になっているので，この点で構造変化を生じたものと考えらえる．すなわち，系全体と局所変形の発生がバルクと比べ比較的近く分かりにくいが，本細線も応力制御でシミュレーションを行えば先の ϵzz = 0.0857 で系全体の不安定モードが優先的に発生するものと考えられる．従って，bcc 鉄の (111) 表面は薄膜・細線いずれにおいても局所変形の発生箇所にはなりにくいと結論づけられる．

33

15

εzz =0.0826
Stress strain curve

1.5

10
1
5
0.5
0

0

-0.5

Number of detB ijα <0 atoms

Stress, σzz , GPa

Average of 6x6 determinant, ΣdetB ijα /N ,normalized by detB ijα at perfect lattice

2

Average of det Bij(Edge Atoms)

αAverage of det Bij(All Atoms)

α

-5

-1

15000

All Atoms

10000
Edge Atoms

5000
0
0

0.05

0.1

0.15

0.2

Strain, εzz(a) Outline
εzz =0.0826

Stress strain curve

0.4

10

εzz =0.0950

εzz =0.0857

5

0

0

10

0.2

εzz =0.0950
Average of detB ijα(Edge Atoms)

5

0

0

εzz =0.0857-0.5 Average of detB ijα(All Atoms)

Average of detB ijα(Edge Atoms)

-5
-0.2

Average of detB ijα(All Atoms)-5

Number of detB ijα <0 atoms

-1
400

15000

All Atoms

All Atoms

300

10000
200

5000
0
0.07

100

Edge Atoms

Edge Atoms

0

0.075

0.08

0.085

0.09

0.095

0.1

0.07

0.075

Strain, εzz

0.08

0.085

0.09

0.095

0.1

Strain, εzz(b) Zoom-up around stress drop
Fig.3.14 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (111)-(110) nanowire.

34

Stress, σzz , GPaεzz =0.0826

1
0.5

15

Stress, σzz , GPa

Average of 6x6 determinant, ΣdetB ijα /N ,normalized by detB ijα at perfect lattice

1.5

Stress strain curve

15

2

z
y

x

ε zz =0.0776

ε zz =0.0885

ε zz =0.0826(peak)ε zz =0.0900

ε zz =0.0853

ε zz =0.0920

Fig.3.15 Snapshot of the (111)-(110) surface nanowire.Red coloured circles
indicated detBijα < 0 atoms.

35

第4章
対称粒界の粒界エネルギーと AES，
ならびに引張シミュレーション
前章では，バルクと異なる不均一構造として，表面および表面会合部の detBijα の分
布や，引張時の系の不安定挙動に与える影響などを検討した．本章では，結晶粒界による不均一構造について同様の検討を行う．

4.1

無負荷平衡状態の粒界エネルギーと AES

4.1.1

解析条件

対称傾角粒界図 4.1 に示すように，[1̄10],[001],[110] をそれぞれ x, y, z 軸とした座標系において，z軸を回転軸として鏡面対象となるように左右の bcc-Fe の単結晶を ±θ/2 回転させて対称傾角粒界を作成する．作成時，粒界面近傍において原子が最近接距離よりも近づきすぎる場合，重なり合う 2 つの原子のうち 1 つを取り除き，残りは粒界面中心に配置
した (図 4.2)．周期境界条件下で粒界同士が相互作用しないようにするため，粒界間の距離は 5[nm] 以上になるように y 軸方向のセル長を設定した．また，Σ 値が 1〜99 の傾
角粒界を解析対象とした．全方向周期境界条件を適用し，初期応力が 0 となるようにセル寸法を調整しながら初期緩和シミュレーションを 50000[fs] 行った．温度は 0.1[K]とした．その後，粒界エネルギーならびに detBijα を評価した．

36

x z
y
Fig.4.1 Schematic of simulation cell for tilt grain boundary.

G.B.

G.B.

x<113>(a) Overlap at grain boundary(b) Removed structurey<332>
z<110>

Fig.4.2 Detail of atom configuration at (332)Σ11 symmetric tilt grain boundary.

対称ねじれ粒界
傾角粒界と同様の座標系において，図 4.3 に示すように z 軸を回転軸とし，±θ/2 回転させて重ね合わせることで対称ねじれ粒界を作成する．ねじれ粒界では粒界面の原子の接近はないので，初期配置において原子の取り除き等は行わない．傾角粒界と同様に，粒界間の相互作用を除くため，粒界間の距離が 5[nm] 以上になるように z 軸方向のセル長を設定した．Σ 値も 1〜99 を解析対象とした．温度の条件などは傾角粒界

37と同じである．

xyz

Fig.4.3 Schematic of simulation cell for twist grain boundary.

38

4.1.2解析結果

対称傾角粒界図 4.4(a)〜(b) の黒線は粒界エネルギーを傾角 θ に対してプロットしたものである．
赤線は，(a) は detBijα の合計値から同じ原子数の完全結晶の detBijα の合計値を引き，粒界面積で割った値，(b) は負の AES 原子数を単位面積当たりの値に換算したものであり，目盛りはそれぞれ右軸に示す．また，低 Σ 値の粒界を矢印で示している．粒界エネ
ルギーは回転角度に依存し，図のように特定の Σ 値 (ここでは，(332)Σ11 と (112)Σ3)
でエネルギーカスプが生じることが報告されている [19] ．特に，(112)Σ3 傾角粒界は整合双晶界面であり，粒界近傍の原子配列にずれが少ないことから粒界エネルギーが最も低くなっている．図 (a) を見ると，エネルギーカスプが現れた粒界モデルに対応して，AES の値も他のモデルと比較して極大点 (バルクの値 1.0 に近い側) を示している．また，図 (b) を見ると，単位面積当たりの負の AES 原子が 0 またはほぼ 0 の粒界は θ=31.5,34.8,70.5((111)Σ3),86.6,109.4((112)Σ3) の 5 つである．逆に (332)Σ11 粒界はエネルギーカスプを示したにもかかわらず負の AES 原子が多く，そのため図 (a) でもAES が極大となっていない．

1
0.5
0

0
-0.5
-1

-0.5
-1.5

0

30

60

90

120

150

1.5

Grain boundary energy [J/m2]

0.5

Defference of 6x6 determinant,(Σ(detBαij )Bulk-Σ(detBαij )GB laminate)/2A[/m2]

Grain boundary energy [J/m2]

1.5(221) (332) (111) (112)(113) (114)
Σ9 Σ11
Σ3 Σ3 Σ11 Σ9

2

50
1
0.5
0

0
-0.5
-1

Number of negative AES atoms
per unit area [/nm2].

(221) (332) (111) (112)(113) (114)
Σ9 Σ11
Σ3 Σ3 Σ11 Σ9

2

-50
-1.5

180

Misorientation angle [degrees]

0

30

60

90

120

150

180

Misorientation angle [degrees]

Grain boudnary energy
AES

Grain boudnary energy
AES

(a)

(b)

Fig.4.4 Relationships between grain boundary energy and misorientation angle , diﬀerence of AES and number of negative AES atoms per unit
area of tilt grain boudnary.

39対称ねじれ粒界
図 4.5 に対称ねじれ粒界の解析結果をまとめて示す．粒界エネルギーを見ると，Σ3の粒界モデルでエネルギーカスプが現れ，その点に対応して AES に極大点が生じている．先の傾角粒界では，粒界エネルギー・AES ともに回転角によって複雑に変化していたが，ねじれ粒界では Σ3 の粒界を除いてはエネルギーに大きな差はなく，AES もほ
ぼ一定となった．(b) を見ると，Σ3 粒界のみ負値の原子が存在していない．また図 (a)のようにバルクからの差として平均した場合はわからないが，図 (b) では θ=50.4(Σ11)
や 90.6，129.5 は負の原子数の割合が大きいことが分かる．
Σ3 Σ11Σ9

0.5
1
0.25

0

0

-0.5

-0.25

-1
-0.5
-1.5

0

30

60

90

120

150

180

1.5
50
1
0.5
0

0

-0.5
-1
-50
-1.5

0

30

60

90

120

150

Misorientation angle [degrees]

Misorientation angle [degrees]

Grain boudnary energy
AES

Grain boudnary energy
AES

180

(b)

(a)

Fig.4.5 Relationships between grain boundary energy and misorientation angle, diﬀerence of AES and number of negative AES atoms per unit
area of twist grain boudnary.

40

Number of negative AES atoms
per unit area [/nm2].

1.5

Grain boundary energy [J/m2]

Grain boundary energy [J/m2]

0.75

0.5

Σ3 Σ11Σ9

Σ9Σ11 Σ3
2

Defference of 6x6 determinant,(Σ(detBαij )Bulk-Σ(detBαij )GB laminate)/2A[/m2]Σ9Σ11 Σ32

4.2

対称傾角粒界の引張シミュレーション

4.2.1

解析条件

引張シミュレーションを行った傾角粒界モデルの諸条件を表 4.1 に示す．いずれも，
全方向周期境界条件，温度 0.1[K] の下で初期構造緩和を 50000[fs] 行った後，粒界面に垂直な方向 (y 方向) に毎ステップひずみ増分 ∆εyy =1.0×10−6 を与えて引張シミュレー
ションを行った．この時，∆εyy =1.0×10−4 毎に原子座標データを記録し，各原子位置ααの原子応力 σkl
ならびに原子弾性係数 Cklmn
を求め，AES の行列式 detBijα を評価した．

緩和計算後の各原子のポテンシャルエネルギーの分布 (図 4.6(a)) に示すように，粒界から数原子層程度が結晶内部よりも高いエネルギーを示す．以降の解析結果では，図4.6(b) に赤で示したようにこれら初期緩和後の高エネルギー原子および，1 原子層分をGBatoms として登録し，粒界とバルク部に分けた議論を行う．
Table 4.1 Simulation condition for symmetric tilt grain boundary under tension．
G.B. patturnΣ11(332) Σ3(111) Σ3(112)
Misorientation [deg]
50.4
70.5
109.4
Cell length:Lx [Å]
113.88105.35109.21
Ly
214.30217.06224.87Lz
105.57105.43105.36
Number of atoms
218400204360 219648

41

G.B.

x
y
(a) energy

G.B.

(b) G.B. atom distinguished by potential energy at initial equilibrium
Fig.4.6 Snapshot coloured by potential energy and grain boundary region ofsymmetric tilt grain boundary.

42

4.2.2(332)Σ11 傾角粒界

(332)Σ11 対称傾角粒界の引張シミュレーション結果を図 4.7 に示す．前章同様，応力ひずみ，detBijα の平均，負の AES 原子数変化を示している． (332)Σ11 粒界は粒界エネ
ルギーはカスプを示すが，負の AES 原子の割合が高い粒界である．応力ひずみ応答を見ると，ϵyy = 0.0611 付近で応力上昇が一度頭打ちし，減少することなく ϵyy = 0.080 付近
まで横ばいとなる．その後再び応力が上昇し，ϵyy = 0.200 で応力ピークを迎え，直後に応力が急減する．AES 変化に着目すると，応力が頭打ちになる ϵyy = 0.0611 で系全体のdetBijα の減少が底打ちし，応力と同様にほぼ横ばいとなる．粒界近傍の detBijα の平均値
はこの応力折れ曲がり点より前の ϵyy = 0.0474 で回復している．その後，ϵyy = 0.0650
以降は減少し続け，応力が再び上昇しひずみ硬化を示す ϵyy = 0.1530 付近でバルクと粒界の AES に差が無くなる．図 4.8 に応力が頭打ちする点近傍の原子座標を示す．(i) は負の AES 原子を赤く着色して示したもの，(ii) は可視化ソフト Atomeye[40] の Central

Symmetry 値によって「欠陥」と判別された原子を示す．ϵyy = 0.0611 以降は粒界から欠陥と判定された原子が結晶内部に伝ぱしている．負の AES 原子はこの伝ぱする欠陥の前縁に観察される．ϵyy = 0.151 近傍から応力が再び上昇 (ひずみ硬化) するのは構造変化が系全体に伝播した点であり，ϵyy = 0.200 のピーク後は図 4.9 に示すように欠陥原子が系全体に伝播し接合した付近からへき開が生じて破断する．

43

20

εyy=0.0474
0αAverage of detBij(GB Atoms)
0

0.05

0.1

0.15

6000

0.2

-10

All Atoms

4000
2000

GB Atoms

0
0

0.05

0.1

0.15

0.2

Stress strain curve30

εyy=0.0474 εyy=0.0611

20

10

0

0.25

Strain, εyy0αAverage of detBij(GB Atoms)

-0.5
0.045

0.25

Number of detBαij <0 atoms

-0.5

1

0.5

10

0

Average of detBαij(All Atoms)

Stress, σyy, GPa

εyy=0.0611

0.5

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

1

40

1.5

40

Average of detBαij
Stress strain curve(All Atoms)
εyy=0.200 30
Stress, σyy, GPa

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

1.5

0.05

0.055

0.06

-10
0.065

0.07

All Atoms
6000
4000

GB Atoms

2000
0
0.045

0.05

0.055

0.06

0.065

0.07

Strain, εyy(a) Outline(b) Zoom-up around stress drop

Fig.4.7 Change in the average，standard deviation of detBijα ，and stressstrain curve of (332)Σ11 tilt grain boundary.

44

G.B.

G.B.

x

yy =0.0474

yyy =0.0611

yy =0.0650(i)Negative atoms(ii)Central symmetry parameterFig.4.8 Snapshots of (332)Σ11 symmetric tilt grain boundary around the 1ststress peak.

45

x
y

yy =0.133

connect

yy =0.151

cleavage(i)Negative atoms

yy =0.190

yy =0.220(ii)Central symmetry parameterFig.4.9 Snapshots of (332)Σ11 symmetric tilt grain boundary around strain
hardening point and cleavage fracture .

46

(111)Σ3 傾角粒界

4.2.3(111)Σ3 対称傾角粒界の引張シミュレーション結果を図 4.10 に示す．また，主要なひずみでの AES，欠陥原子配置を図 4.11 に示す．(111)Σ3 粒界は粒界エネルギーが高いが，負の AES 原子が無負荷平衡状態にはほぼ存在しない粒界である (図 4.4)．応力
ひずみ応答を見ると，応力ピークを示す ϵyy = 0.099 までに ϵyy = 0.020 付近でわずかな折れ曲がりが見られるが，その点では粒界近傍の detBijα の平均値が回復している．
この時，粒界近傍に負の AES 原子が図 4.11(a) に示すように瞬間的に現れている．同様の瞬間的急増は応力ピーク直前の ϵyy = 0.090 近傍で見られ，この時の負の AES を
図 4.11(b) に示すように粒界面近傍である．図 4.10(b) の応力ピーク近傍の拡大図を見ると，ϵyy = 0.0958 以降負の AES 原子が出現しているが粒界は 1000 程度，全体では2000 程度で頭打ちとなっている．図 4.11(c)〜(f) でこの時の変化を見ると図 (d) のよ
うに粒界からキンク上に欠陥原子が発生し，粒界から数原子層ほど欠陥原子が伝播す

30

εyy=0.0990

εyy=0.020

Stress, σyy, GPa

20

-0.5

0

Average of detBαij(GB Atoms)
0

0.05

0.1

0.15

-10
0.2

All Atoms
2000
GB Atoms
1000
0
0

0.05

0.1

0.15

0.2

Strain, εyy

30

20

Average of detBαij(All Atoms)

10

0

0

Average of detBαij(GB Atoms)

-0.5
0.09

3000

εyy=0.0990

1

0.5

10

0

εyy=0.0958

Stress, σyy, GPa

Average of detBij(All Atoms)
1

Stress strain curve

1.5

α

0.5

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect latticeStress strain curve

1.5

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect latticeる．しかし，内部に欠陥が伝播することなく (図 (e))，その後粒界面からへき開した．

0.095

0.1

0.105

0.11

0.115

-10
0.12

3000
All Atoms
2000
GB Atoms

1000
0
0.09

0.1

0.11

0.12

Strain, εyy(a) Outline(b) Zoom-up around stress drop

Fig.4.10 Change in the average，standard deviation of detBijα ，and stressstrain curve of (111)Σ3 tilt grain boundary.

47

G.B.

x
y

(a) yy =0.020(b) yy =0.090

G.B.

G.B.

(c) yy =0.0958(d) yy =0.0961

cleavage(e) yy =0.0990(peak)(f) yy =0.120(i) Negative atoms(ii) Central symmetry parameterFig.4.11 Snapshots colored by distribution of negative AES atoms and defeat
atoms distinguished by central symmetry parameter around the 1ststress peak ((111)Σ3 tilt grain boundary).
48

(112)Σ3 傾角粒界

4.2.4(112)Σ3 対称傾角粒界の引張シミュレーション結果を図 4.12 及び図 4.13 に示す．応
力ひずみ応答を見ると，ϵyy = 0.0701 付近で一度応力ピークを迎えた後に，応力が急減することなく横ばいになり，先の (332)Σ11 傾角粒界と似た傾向を示す．しかし，ひず
み増加に伴って応力は上昇し続け，ϵyy = 0.3 まで引っ張っても応力急減挙動は示さなかった．(112)Σ3 対称傾角粒界は粒界面が整合双晶であるため，欠陥原子はほぼ表示されず，表示されている原子に関しても原子配置のずれはごくわずかである．(332)Σ11
粒界と同様，ϵyy = 0.0725 以降，粒界近傍から内部へと双晶面が移動していることが示唆される．(332)Σ11 粒界と異なる点として，移動した界面が周期境界条件によって

ϵyy = 0.221 付近で接合した後，引っ張りによって全ての原子の AES が正となり，へき

1.5

Average of detBij(All Atoms)30

εyy=0.0190 εyy=0.0701
Stress, σyy, GPa

1

20

0
α

Average of detBij(GB Atoms)

-0.5

0

0.05

0.1

0.15

1

εyy=0.0673

Stress strain curve 30

εyy=0.0701
20

εyy=0.0763

0.5

10

0

Average of detBαij(All Atoms)

Stress, σyy, GPa

α0.5

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect latticeStress strain curve

1.5

-10
0.2

10

0

0
α

-0.5
0.06

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

開は生じなかった．

Average of detBij(GB Atoms)
-10
0.065

0.07

0.075

0.08

0.085

0.09

15000

15000
All Atoms

All Atoms

10000

10000
5000

GB Atoms

0
0

0.05

0.1

0.15

0.2

5000
0
0.06

GB Atoms

0.065

0.07

0.075

0.08

0.085

0.09

Strain, εyyStrain, εyy(b) Zoom-up around stress drop(a) Outline

Fig.4.12 Change in the average，standard deviation of detBijα ，and stressstrain curve of (112)Σ3 tilt grain boundary.

49

G.B.

G.B.

x

yy =0.0673

yyy =0.0701

yy =0.0725

connect

yy =0.0763

yy =0.221(i)Negative atoms

yy =0.300(ii)Central symmetry parameterFig.4.13 Snapshots colored by distribution of negative AES atoms in (112)Σ3
tilt grain boundary. Red coloured circles indicate detBijα < 0 atoms.

50

第5章
対称傾角粒界のせん断シミュレー
ションならびに AES 解析
3 章，4 章では，表面や粒界等の不均一構造を対象として，引張変形下の不安定挙
動 (応力急減) について検討を行った．本章では前章に導入した対称傾角粒界に対して，異なる変形モードであるせん断変形の際の不安定挙動を観察し，そのときの各原子のdetBijα について検討した．

5.1

解析条件

前章で作成したシミュレーションセルの諸条件を表 5.1 ならびに図 5.1 に再掲する．
いずれも x, z 方向周期境界，y 方向自由境界の条件下で 50000[fs] の初期緩和シミュレー
ションを行った後に，上下 3 原子層をつかみ部として固定しながら yz 方向にせん断ひずみを毎ステップ増加させてせん断シミュレーションを行った．せん断ひずみ増分は∆εyz =1.0×10−6 /[fs] とし，制御温度は 0.1[K] である．この時，∆εyz =1.0×10−4 毎に原ααを求ならびに原子弾性係数 Cklmn
子座標データを記録し，各原子位置の原子応力 σkl

め，AES の行列式 detBijα を評価した．

51

Table 5.1 Simulation condition for symmetric tilt grain boundary under shear．
G.B. patturnΣ11(332) Σ3(111) Σ3(112)
Misorientation [deg]
50.4
70.5
109.4
Cell length:Lx [Å]
113.88105.35109.21
Ly
214.30217.06224.87Lz
105.57105.43105.36
Number of atoms
218400204360 219648

chunk

G.B.
y
x
z

G.B.

[332]

[113]

[112]

[111]

[110]

[112]
(a)Σ11(332)

G.B.

G.B.

[111]

[110]
(b)Σ3(111)

Fig.5.1 Schematic of shear simulation.

52

[110]
(c)Σ3(112)

5.2

解析結果

5.2.1

(332)Σ11 傾角粒界

図 5.2 に (332)Σ11 傾角粒界の解析結果を示す．前章と同様，粒界付近の高エネルギー
原子を GB Atoms として分けて評価した．応力ひずみ応答を見ると，ϵyz = 0.1381 でピークを迎えるが，その後応力は急減することなくなだらかに減少する．応力ピーク付近の detBijα 変化に着目すると，系全体ならびに粒界の平均値は共に正値であるが，粒界近傍の平均値が応力ピークより前の ϵyz = 0.1173 付近で急回復し，ピーク直前のϵyz = 0.1199 では系全体の平均値と逆転している．この変化によって，初期に存在した粒界近傍の負の AES 原子数が急減し，ϵyz = 0.1225 では 0 となっている．
図 5.3 及び 5.4 に (332)Σ11 傾角粒界の原子配置変化 (AES 及び欠陥原子) を示す．粒界の detBijα 平均値が回復する以前のひずみ ϵyz = 0.1173 では，粒界近傍に負の AES 原子が多数存在し，粒界部分の原子層が欠陥と判別されている．ϵyz = 0.1255 の図では初期の粒界部分に負の AES 原子は存在せず，その上下に負の AES 原子が現れている．
欠陥原子はその後拡大し応力ピークでは幅 9 原子層程度となっている．ϵyz = 0.250 までせん断を行ったが，欠陥原子はさらに上下に拡大し結晶内部のほとんどで構造変化が生じている．

53

15
10
5

-0.5

Average of detBαij(GB Atoms)
0

0.05

0.1

0

Stress strain curve
-5
0.15

0.2

GB Atoms

5000

0
0

0.05

1

εyz=0.1199

0.1
0.15
Strain, εyz

0.2

0.25

(a) Outline

εyz=0.1255

εyz=0.1381

20
15

0.5

10

0

Average of detBαij 5(All Atoms)
0

10000

All Atoms

10000

Stress strain curve 25

-0.5
0.115

-10
0.25

Number of detBαij <0 atoms

030

Stress, σyz, GPa

20

0.5

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice25

εyz=0.1381

1

1.5

30

Average of detBαij(All Atoms)

Stress, σyz, GPa

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

1.5

Average of detBαij(GB Atoms)
0.12

0.125

0.13

-5
0.135

0.14

-10
0.145

All Atoms

5000

GB Atoms

0
0.115

0.12

0.125

0.13 0.135
Strain, εyz

0.14

0.145

(b) Zoom-up around stress drop

Fig.5.2 Change in the average，standard deviation of detBijα ，and stressstrain curve of (332)Σ11 tilt grain boundary.

54

G.B.

y

x

z

yz =0.1173

yz =0.1199

yz =0.1255

G.B.

yz =0.1381(peak)yz =0.1400

yz =0.2500

Fig.5.3 Snapshots of (332)Σ11 symmetric tilt grain boundary under shear.
Red coloured circles indicate detBijα < 0 atoms.

55

G.B.

y

x

z

yz =0.1173

yz =0.1199

yz =0.1255

G.B.

yz =0.1381(peak)yz =0.1400

yz =0.2500

Fig.5.4 Snapshots of (332)Σ11 symmetric tilt grain boundary under shear(central symmetry parameter).

56

(111)Σ3 傾角粒界

5.2.2

図 5.5 に (111)Σ3 傾角粒界の解析結果を示す．応力ひずみ応答を見ると，先の粒界と異なり ϵyz = 0.1300 で急激な折れ曲がりを示し急減する．また，ϵyz = 0.1074 付近でわずかな折れ曲がり点が存在し，直後に粒界近傍の負の AES 原子数が急減している．
応力ピーク直後では，系全体及び粒界付近で負の AES 原子数が急増している．図 5.6に最小のわずかな折れ曲がり (ϵyz = 0.1074) 近傍の変化を AES ならびに欠陥原子に分けて示す．折れ曲がり点前では以前では粒界近傍に負の AES 原子が存在し，粒界部分が欠陥原子と判定されている．折れ曲がり後，粒界近傍の AES が回復し，1 原子層結
晶粒側に欠陥原子が拡大している (ϵyz = 0.1091)．大きな応力低下する前後の AES 分
布，欠陥原子分布を図 5.7，5.8 に示す．応力ピーク時には粒界から 3 原子層程度に負のAES 原子が存在しているが，応力急減に伴い粒界を起点として結晶内部へと欠陥 (転εyz=0.130

15

1
10

0.5
5

Average of detBαij(GB Atoms)
0

0.05

0.1

-5
0.15

15

1

Stress strain curve
εyz=0.130
Average of detBαij
10(All Atoms)
5

0

0.2

0

Average of detBαij(GB Atoms)

-5

-0.5
0.1

Number of detBαij <0 atoms020

εyz=0.1074

0.5

0Stress strain curve-0.5

Number of detBαij <0 atoms

1.5

20

Stress, σyz, GPa

Average of detBαij(All Atoms)

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

εyz=0.1074

1.5

Stress, σyz, GPa

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice位) を生じている．

0.11

0.12

0.13

0.14

0.13

0.14

10000

10000

5000

All Atoms

GB Atoms

0
0

0.05

0.1
Strain, εyz

0.15

0.2

All Atoms
5000

0
0.1

GB Atoms
0.11

0.12
Strain, εyz(b) Zoom-up around stress drop(a) Outline

Fig.5.5 Change in the average，standard deviation of detBijα ，and stressstrain curve of (111)Σ3 tilt grain boundary.

57

yz =0.1074

yz =0.1081

yz =0.1091

G.B.

y

x

z

(i) AES

G.B.

(ii) Central symmetry parameter

Fig.5.6 Snapshots colored by distribution of negative AES atoms and defeat
atoms distinguished by central symmetry parameter around ϵyz =0.1074 ((111)Σ3 tilt grain boundary).

58

G.B.

y

x

z

yz =0.1300(peak)yz =0.1305

yz =0.1310

G.B.

yz =0.1350

yz =0.1400

yz =0.1500

Fig.5.7 Snapshots colored by distribution of negative AES atoms aroundstress drop at ϵyz = 0.13 ((111)Σ3 tilt grain boundary).

59

G.B.

y

x

z

yz =0.1300(peak)yz =0.1305

yz =0.1310

G.B.

yz =0.1350

yz =0.1400

yz =0.1500

Fig.5.8 Snapshots of (111)Σ3 symmetric tilt grain boundary under shear
around stress drop at ϵyz = 0.13 (central symmetry parameter).

60

(112)Σ3 傾角粒界

5.2.3

図 5.9 に (112)Σ3 傾角粒界の解析結果を示す．応力ひずみ応答を見ると，ϵyz = 0.2201
で応力ピークを示すがさらに応力急減が顕著となっている．負の AES 原子が現れるのは，ϵyz = 0.178 であるが，応力ピーク直前の ϵyz = 0.192 でさらに増加して粒界から 1
原子層がすべて負の AES 原子となる．図 5.10，5.11 に (112)Σ3 傾角粒界の応力ピーク前後の原子配置をまとめて示す．先述のように，ϵyz = 0.178〜0.192 の変化で粒界から1 原子層が負となっているのが図 5.10(a) で確認できる．先の図 5.9(b) の右下拡大図において，ピーク後粒界以外の原子が負となっていることがわかるが，図 5.10(b)(c) お
よび 5.11(b)(c) からわかるように，粒界ではなく結晶内部から欠陥生成に対して前駆

30
20

0
α

-0.5

0

0.05

60000

0.1

0.15

0.2

-10
0.25

10
0
α

Average of detBij(GB Atoms)
0.215

60000

All Atoms

20

Average of detBαij(All Atoms)

0

-0.5
0.21

Number of detBαij <0 atoms

Average of detBij(GB Atoms)

40
30

0.5

10

Stress strain curve

1

Stress strain curve

εyz=0.2201

Stress, σyz, GPa

40

1

0

50

1.5

50

Average of detBαij
εyz=0.2201(All Atoms)

0.5

Number of detBαij <0 atoms

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

1.5

Stress, σyz, GPa

Average of 6x6 determinant, ΣdetBαij /N,normalized by detBαij at perfect lattice

的に負の AES 原子が発生していることが確認できる．

0.22

0.225

0.23

-10

0.235

0.24

0.235

0.24

All Atoms

45000

45000
30000

30000

GB Atoms

GB Atoms

15000

15000
0
0

0.05

0.1
0.15
Strain, εyz

0.2

0.25

0
0.21

0.215

0.22

0.225 0.23
Strain, εyz(b) Zoom-up around stress drop(a) Outline

Fig.5.9 Change in the average，standard deviation of detBijα ，and stressstrain curve of (112)Σ3 tilt grain boundary.

61

G.B.

y

x

z

(a) yz =0.1920(b) yz =0.2200(c) yz =0.2201(peak)G.B.

(d) yz =0.2207(e) yz =0.2235(f)yz =0.2260

Fig.5.10 Snapshots colored by distribution of negative AES atoms aroundstress drop at ϵyz = 0.2201 ((112)Σ3 tilt grain boundary).

62

G.B.

y

x

z

(a) yz =0.1920(b) yz =0.2200(c) yz =0.2201(peak)G.B.

(d) yz =0.2207(e) yz =0.2235(f)yz =0.2260

Fig.5.11 Snapshots of (112)Σ3 symmetric tilt grain boundary under shear
around stress drop at ϵyz = 0.2201 (central symmetry parameter).

63

第6章
結論
本研究では，表面薄膜やナノワイヤ，結晶粒界等の不均一構造を有した bcc 構造の鉄に対して，分子動力学法を用いて引張・せん断変形シミュレーション行い，不安定
挙動と原子弾性剛性係数 Bijα の行列式 detBijα の関係について検討した．以下に，本論文で得られた結果を総括する．
第２章では解析手法の概要を述べた．はじめに分子動力学法の基礎方程式を示し，原子間相互作用の評価に用いた FS ポテンシャルの概要を説明した．さらに，各原子の安定性評価の指標となる原子弾性剛性係数について概説した．
第３章では表面薄膜ならびにナノワイヤ構造を有したモデルについて分子動力学法
による引張シミュレーションを行い，局所変形と系全体の変形についての評価を検討した．不均一構造の解析前に，単結晶バルクの [001] 方向引張シミュレーションを行う
ことで，bcc 鉄の [001] 引張では局所変形よりも系全体の不安定に達することを明らかにした．(100),(110),(111) の 3 つの表面薄膜モデルについて引張を行った結果，(100)
表面薄膜は，系全体で平均した detBijα (系全体の安定性) が正のときに局所変形を生じて応力急減した．(100)-(010) ナノワイヤや (110)-(1̄10) ナノワイヤでも同様に，系の不安定より先に表面からの局所変形によって応力急減することが示された．一方，(110)
表面薄膜，(111) 薄膜および (111)-(1̄10) ナノワイヤは表面からの局所変形発生より前に系全体の不安定に達することが示された．
第４章では様々な Σ 値の対称傾角粒界ならびに対称ねじれ粒界を作成し，無負荷平衡状態における粒界エネルギーと，detBijα の平均や detBijα < 0 の割合の回転角度依64

存性について議論した．エネルギーカスプを示した粒界 ((112)Σ3 傾角粒界，(110)Σ3
ねじれ粒界) では detBijα の平均値に極大点が見られ，バルクとの差が小さかった．しかし，エネルギーカスプを示した (332)Σ11 傾角粒界では負の AES 原子が多数存在し，detBijα の平均値で見るとエネルギーカスプと対応は見られない．次に，3 つの傾角粒界((332)Σ11,(111)Σ3,(112)Σ3) について粒界面に垂直に引張シミュレーションを行った．
ピラミッド型の応力ひずみ曲線を示して粒界でへき開した (111)Σ3 傾角粒界では，粒界近傍で負の detBijα 原子が出現し欠陥原子が 1 原子層程粒界から拡張しただけであった．(332)Σ11 粒界は応力ひずみに折れ曲がりを示した後ひずみ硬化し最終的に粒界でへき開したが，折れ曲がり点以降欠陥が粒界から粒内を伝ぱすること，その前縁の原子の detBijα が負となっていることが明らかになった．(112)Σ3 粒界では同様のひずみ硬化挙動を示し，ϵyy = 0.25 まで引っ張ってもへき開しなかった．Central symmetry

parameter による欠陥判別では双晶面は認識できないが，負の detBijα 原子変化を見ると先の (332)Σ11 粒界と同様に折れ曲がり点以降粒界面から内部へ伝ぱしていることが確認された．
第 5 章では先の 3 つの対称傾角粒界に対して異なる変形モードであるせん断変形を行い，不安定挙動について検討した．(332)Σ11 傾角粒界は引張変形と同様，応力ピー
ク以降に粒界面が上下に移動し，欠陥部分が拡大するような変化を生じた．(111)Σ3 粒界は粒界を起点とした欠陥生成 (転位) を生じたが，最も急激に応力低下した (112)Σ3
粒界では粒内に負の detBijα 原子が発生し，その箇所に欠陥が生じた．

65

参考文献
[1] A. Miyauchi, K. Kuwabara, and M. Hasegawa presented at 1st Int. Conf. on
Nanoimprint and Nanoprint Technology , (2002)[2] K. Kuwabara and A. Miyauchi，J. Vac. Sci. Technol. B，26 ,(2008)，582[3] K. Kuwabara, A. Miyauchi, H, Sugimura，Jpn. J. Appl. Phys，49, (2010), 106505[4] M.Tabuchi，Y.Baba，Anal.Chem，77，(7090)，2005[5] N.Chekurov，K Grigoras，A Peltonen，Nanotechnology，20，(2009)，065307[6] V. M. Segal, V. I. Reznikov, A. E.Drobyshevskiy, and V. I. Kopylov, RussianMetallurgy (English trans.), 1, (1981), 115[7] V. M. Segal, V. I. Reznikov, V. I. Kopylov, D.A. Pavlik, V. F. Malyshev, Processesof Plastic Structure Formation of Metalls , (1994), 232.
[8] V. M. Segal, Mater. Sci. Eng. , A17, (1995), 157[9] V. M. Segal, Mater. Sci. Eng. , A345, (2003), 36[10] V. M. Segal, Mater. Sci. Eng. , A271, (1999), 322[11] L. Lu, M. L. Sui, and K. Lu, Science, 287, (2000), 1463
[12] Y.Gan，J.K.chen，Mechanics Reserch Communications，36-7，(2009)，838[13] Y.Tong，F.Yi，Lisheng.Liu，Pengcheng.Zhai，Qinging.Zhang，Phisica B: Condensed Matter，405-15，(2010)，3190[14] H.A.Wu，Mechcanics Reserch Communications，33-1，(2006)，966

[15] Y.H.Wen，Zi-Zhong.Zhu，Ru-Zeng.Zhu，Computational Materials Science，414，(2008)，553[16] Y.Jing，Qingyuan.Meng，Wei.Zhao，Physica E:Low-dimensional Systems andNanostructures，41-4，(2009)，685[17] M.Grujicic，P.Dang，Matrials Science and Engineering，201-1，(1995)，194[18] Jyh-Tyng，K.Ikeda，J.Japan Inst.Metals，69-3，(2005)，303[19] 中島・竹内，鉄と銅，86-5，(2000)，73[20] 池田，九州大学情報基盤センター広報 全国共同利用版，6-3，(2007)[21] T.Mizuguchi，N.Tanaka，J.Japan Inst.Metals，69-11，(2009)，1010[22] Ru-Zhi Wang，M.Kohyama，J.Japan Inst.Metals，73-8，(2009)，559[23] T.Kitamura，Y.Imeno，J..Soc.Mat.Sci.Japan，48-9，(1999)，992[24] Douglas E.Spearot，Mark A. Tschopp，Acta Materialia，55，(2007)，705[25] M.A.TSchopp，K.N.Solanki，PHYSICAL REVIEW B，85，(2012)，064108[26] S.P. Chen,D.J. Srolovitz,A.F. Voter，J. Mater. Res.，4-1，(1989)，62[27] M.A.TSchopp，D.L. Mcdowell，International Journal of Plasticity，24，(2008)，
191[28] M.A.TSchopp，G.J. Tucker, D.L. Mcdowell，Computational Materials Science，
44，(2008)，351[29] M.A.TSchopp，D.L. Mcdowell，Philosophical Magazine，87-25，(2007)，3871[30] J.A. Brown，Y. Mishin，PHYSICAL REVIEW B，76，(2007)，134118[31] 石田，UDC，4，(1967)，539

67

[32] 屋代・冨田，日本機械学会論文集 (A 編)，67,-656 (2001)，678[33] 屋代・西村・冨田，材料，54(10)，(2005)，1053[34] 屋代・藤原，材料，60(11)，(2011)，968[35] Finnis，M.W.andSinclair,J.E,Phil.Mag, 50(1), (1984)，45[36] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
[37] Wallace, D. C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
[38] Wang, J., S.Tip, S.R.Phillpot and Wolf, D., Phys.Rev.lett., 71-25, (1993), 4182
[39] Wang, J., Li, S.Tip, A.R.Phillpot and Wolf, D., Phys. Rev., 52-17B(1995), 12627[40] Li,J.，Modelling Simul. Mater. Sci. Eng., 11, (2003), 17368学術論文・学術講演
学術講演

 bcc 鉄の表面・エッジ部の弾性限界：原子弾性剛性係数による検討
山根

尭，屋代

如月

第 17 回分子動力学シンポジウム，東京大学生産技術研究所 (2012.6)

 原子弾性剛性係数 (AES) による bcc-Fe 鉄の結晶粒界の特性評価
山根

尭，屋代

如月

第 25 回計算力学講演会，ポートアイランド南地区 (2012.10)69

論文貼り付け用ページ70

論文貼り付け用ページ

71

論文貼り付け用ページ72

論文貼り付け用ページ

73

論文貼り付け用ページ

74

論文貼り付け用ページ75

謝辞
本研究を遂行するにあたり，終始，仔細に渡る懇切丁寧な御指導を賜りました屋代
如月准教授に, 心より深く感謝の意を表します．また，本研究に対して貴重な御助言を頂いたのみならず，日頃の研究生活に多くの御指導・御助言を頂きました田中克志教
授に深く感謝致します．
本研究室配属後，沢山のアドバイスを下さり，丁寧に教えて下さった表面・界面工学研究室の諸先輩方に厚く御礼申し上げます．特に，浅識な著者にプログラム，材料等に関して幅広く教えて下さった山口明宏様 (現 三菱重工業株式会社)，藤原正大様(現 株式会社神戸製鋼所)，池宮一繁様 (現 オークマ株式会社) には更なる感謝の意
を申し上げます．
また，互いに切磋琢磨しながら研究に取り組み，苦楽を共にした同輩の上野彰宣君，
永田祥一郎君，河野良祐君，田中健君，寺田駿君，平田和幸君に御礼申し上げます．
最後に．6 年間の大学生活をあらゆる面から支えて頂いた家族には，深く感謝致します．ありがとうございました．
平成 26 年 2 月
山根76

尭