修士論文

原子弾性剛性係数に基づく局所格子不安定性解析：
fcc 金属の表面・界面力学特性評価

指導教員：屋代 如月

池宮 一繁

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

Master thesis

Local lattice instability analysis
based on atomic elastic stiﬀness:
Mechanical properties of surface
and interface in fcc metals

Kazushige IKEMIYA
Feburary 2013
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan要約
fcc 金属の Al ならびに Ni において，表面や粒界という初期不均一構造が引張時の系の不安定挙動=応力急減におよぼす影響を明らかにすることを目的として，原子弾性剛性係数 detBijα の正値性という独自の観点から検討した．まず，表面による構造不均一性を有する系の不安定挙動を議論した．(100),(110),(111) 表面の Al，Ni 薄膜およびナノワイヤについて検討を行い，Al の (111) 薄膜，Ni の薄膜全てと (110) 細線は，応力ひずみのピーク点付近で detBijα の平均が負となる「系の不安定条件」に達して応力急減=引張不安定を生じた．それ以外のナノワイヤでも Ni は全体的に detBijα の平均は 0に近い値をとり，系の不安定に近づいたとき変形を生じることがわかった．Al では，
応力ピーク前に detBijα の平均値は大きな正の値をとることが多く，Ni とは異なり系の不安定挙動の予測が難しい．次に様々な角度の対称ねじれ粒界ならびに対称傾角粒界を作成し，ミスオリエンテーション角度 θ に対する粒界エネルギーならびに，detBijαの平均や最小値，detBijα ＜0 の数を無負荷平衡状態で調べた．Al，Ni ともに，Σ3 ねじ
れ粒界，(111)Σ3，(311)Σ11 傾角粒界においては低いエネルギーカスプを示し，detBijα＜0 の原子が存在しないことを明らかにした．それ以外の角度では detBijα の最小値やdetBijα ＜0 の数は必ずしも粒界エネルギーと対応しなかった．さらに Σ3 の粒界構造について，粒界に垂直方向の引張を分子動力学法により行い，粒界固有の強度について

detBijα の視点から検討した．Σ3 ねじれ粒界の引張では Ni では初期に粒界にキンクが発生した．Al，Ni とも応力ピーク時には detBijα の平均が 0 となる「系の不安定条件」に達しており，粒界から転位が発生した．(111)Σ3 傾角粒界の引張では，Al はやはり系の不安定に達していたが Ni では detBijα の平均はわずかに正値であった．Al，Ni ともに応力ピーク後に粒内に欠陥原子が集団的に発生し，相変態的な挙動を示していた．

Summary
In order to clarify the inﬂuence of inhomogeneity by surface or interface on the stressdecrease on the system instability under tensilon, we investigated the unstable behavior
or “elastic limit”, from the unique view point of the atomic elastic stiﬀness (AES) in fcc
metals of Al and Ni．First，we discussed about the eﬀect of structural inhomogeneity
by surface. We performed molecular dynamics simulations on Al and Ni thin plateand nanowire models with (100),(110) and (111) surface．Al (111) thin plate, all of
Ni thin plate and Ni (110) nanowire showed unstable stress drop, when the average of
detBijα became negative，or reached to the “system instability”. Other Ni nanowires
also showed near-zero value in the average of detBijα at the stress-strain peak, while itis diﬃcult to evaluate the instability behavior for other Al models, since they showed
high positive value of the average at the stress-strain peak. Next we made twist andtilt grain boundaries against various misorientation angle, and evaluated the grainboundary energy, the average and minimum of Bijα , the number of detBijα < 0 under
the no-load equilibrium. Both Al and Ni showed low energy cusp for the Σ3 twist,(111)Σ3 and (311)Σ11 tilt grain boundaries and there is no detBijα < 0 atoms inthese grain boundaries. The tendency of the minimum of detBijα and the number of
detBijα < 0 doesn’t coincide with that of the grain boundary energy for the othermisorientation angles. We have then performed molecular dynamics simulations oftension normal to the Σ3 twist and tilt grain boundaries to investigate the inherent
strength of grain boundaries from the view point of detBijα . In the case of twist grainboundaries of Al and Ni, the average of detBijα become negative at the stress-strain
peak, then dislocations nucleated from the grain boundaries. Kink formation wasobserved at the twist grain boundary of Ni in the early stage of tension. In the caseof tilt grain boundaries, Al also showed negative values in the average of detBijα at theelastic limit, while Ni did slightly positive but near zero value. After the stress peak,
we found many “detBijα < 0” atoms in the grains and the stress drop caused by thephase transformation in the grains.

目次
第1章 緒 論1

第 2 章 解析手法の基礎

42.1

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

4

2.2

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

5

EAM ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2.1
2.3

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

72.4

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

8

2.4.1

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

8

2.4.2

原子系への拡張 . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.4.3

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

11

第 3 章 薄膜・細線モデルによる解析

12

3.1

解析モデル及び解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . .

12

3.2

Al の解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.2.1

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

15

3.2.2

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

17

3.2.3

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

19

3.2.4

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

22

3.2.5

(110)-(1̄10) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

25

3.2.6

(111)-(1̄10) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

28

Ni の解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.3.1

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

30

3.3.2

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

32

3.3.3

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

34

3.3

i

3.4

目次

ii

3.3.4

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

36

3.3.5

(110)-(1̄10) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

38

3.3.6

(111)-(1̄10) 細線モデル . . . . . . . . . . . . . . . . . . . . . . .

40

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

第 4 章 粒界の原子弾性剛性係数評価

43

無負荷平衡時の粒界エネルギーと原子弾性剛性係数 . . . . . . . . . . .

43

4.1.1

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

43

4.1.2

解析結果および考察 . . . . . . . . . . . . . . . . . . . . . . . .

45

4.2

Σ3 対称ねじれ粒界の引張シミュレーション . . . . . . . . . . . . . . .

51

4.3

Σ3 対称傾角粒界における引張シミュレーション . . . . . . . . . . . . .

61

4.4

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

4.1

第 5 章 結論

70

参考文献

72

学術論文・学術講演

75

謝辞

81

第1章
緒論
バルクと異なる構造・力学特性を示す表面は，原子の不整合性に起因する表面エネ
ルギーが存在する．巨視的な物体では，体積は表面積に対して非常に大きいことから，物体の表面積効果は大きくない．しかしながら，ナノ材料等のように系が小さくなればなるほど，体積に対する表面積の割合が増加するため，表面（界面）効果を考慮する必要がでてくる．このため，ナノスケールでの材料加工が可能となった現代においては，表面や粒界の原子構造をも考慮した新しい表面・界面物理を模索していく必要がある．
表面がナノスケールで及ぼす影響を明らかにすることを目的として，薄膜や半導体，
純金属など幅広い材料に対してナノインデンテーション試験が行われている [1], [2] ．そこでは荷重-変位関係における不安定現象として，弾性変形から押し込み荷重一定のまま変位が急増する「変位バースト」現象が報告されている [3]–[5] ．実験だけではなくナ
ノ薄膜及びナノワイヤの分子動力学法シミュレーションも広く行われてきた．北村ら
は Ni の微小細線・薄膜の引張シミュレーションを行い，拘束条件による局所的へき開
破壊の発生の可能性を示している [6] ．Y.Gan らは，Cu のナノ薄膜の引張シミュレー
ションを行い，小さいほど柔らかいが，破断強度は高くなることを報告している [7] ．

Y.Tong らは，Bi2 Te3 ナノ薄膜の引張シミュレーションを行い，緩和後にバルクでは
圧縮応力，表面では引張応力が働いていること，厚さ及び温度によるヤング率の変化を報告している [8] ．ナノワイヤの解析では，H.A.Wu らは，Cu ナノワイヤの単軸引張シミュレーションを行い，温度効果，寸法効果による弾性係数変化及び応力分布変化

1を報告している [9] ．Y.H. Wen らは，Ni ナノワイヤの単軸引張シミュレーションを行
い，ひずみ速度の大きさによって構造変化に大きな差が生じることを報告している [10] ．

Y.Jing らはシリコンの単軸引張シミュレーションを行い，限界荷重が系の直径やひず
み速度，温度に大きく影響していること，融解温度が直径の減少に伴って低下することを報告している [11] ．また，表面力学の分野において，表面欠陥がもたらす影響を分子動力学により検討する研究も行われており，M.Grujicic らは，Fe-Ni 合金中のマルテ
ンサイト変態について MD シミュレーションを用いて表面の持つ結晶欠陥がマルテン
サイトの構成に大きく影響していることを報告している [12] ．
また表面だけでなく，結晶材料中に存在する結晶粒界も強度特性に大きく影響し，それが原因で本来延性を示す材料がぜい性的に破壊することが知られている [13] ．そのような粒界に起因した破壊メカニズムの解明のため，高温でのアルミニウムの粒界拡散の解析 [14], [15] や粒界溝の変形解析など原子シュミレーションが広く行われてきた [16]–[18] ．
また粒界構造や粒界エネルギーおよび粒界破壊強度に関するシミュレーションも銅 [19]やモリブデン [20]–[22] ，α-鉄 [23] ，マグネシウム [24], [25] ，およびパラジウム [26], [27] と幅広く行われている．
著者らのグループでは，個々の原子位置における原子弾性剛性係数 (Atomic Elastic

Stiﬀness Coeﬃcient，AES と称している) の正値性に基づいて，局所変形の開始ならびに構造安定性を評価することを提案している．AES による検討は，fcc，アモルファ
ス金属，そしてダイヤモンド構造のシリコンについて進められてきた [28]–[31] ．初期 Ni
ナノワイヤでの検討では [28] ，ナノワイヤのくびれ発生 → 部分転位の発生メカニズム
が AES の正値性で評価できることを報告している．アモルファス構造 (Ni，Al そして

NiAl 合金系) での検討では，AES の正値性によって局所クラスタの崩壊及び再構成が評価できることを報告している．シリコンの検討では，表面 (内部不均一を有する系)の解析で，弾性限界までは引張によって系全体で均一化していることを AES の標準偏差で明らかにしている．一方，FS ポテンシャルによる bcc-鉄の薄膜の引張では，局所変形を生じる前に多くのケースで全ての原子の AES が負となり系全体が不安定になっていた．このことは，bcc 鉄では表面から構造変化=すべりを発生するよりも相変態な
どの全体的な変形が発生しやすいと考えることができる [32] ．このように AES による

2

局所格子評価の試みは必ずしも「負の AES 発生=局所不安定」と単純なルールではな
いものの，様々な変形破壊現象の理解のために着実に進められてきた．
本研究では，fcc 金属の Al と Ni をモデル材料とし，表面ならびに粒界による初期不均一構造が引張時の系の不安定挙動=応力急減におよぼす影響を AES の観点明らかにすることを目的とする．(100)，(110)，(111) などの表面について，無限平板モデルの引張シミュレーションならびに，表面と表面が会合したエッジ部の効果について議論するため，細線モデルの引張シミュレーションを分子動力学法により行う．粒界については，対称ねじれ粒界および傾角粒界について，ミスオリエンテーション角度 θ に対する粒界エネルギー評価に対応させた AES 評価を行うとともに，Σ3 ねじれ・傾角
粒界について引張シミュレーションを行う．
以下に各章の概略を示す．第 2 章では分子動力学法の概要ならびに粒子間相互作用の評価に用いたポテンシャルについて述べる．その後 AES による力学評価手法についても説明する．第 3 章では，fcc 金属の Al と Ni について，引張時の系の不安定挙動=
応力急減におよぼす表面不均一性の影響を調べるため，(100)，(110)，(111) 表面を有する薄膜や細線の引張シミュレーションを行って応力-ひずみ曲線を求めるとともに，原子弾性剛性係数の行列式（detBijα ）の正値性の観点から検討する．第 4 章では様々な角度の対称ねじれ粒界ならびに対称傾角粒界を作成し，無負荷平衡状態における回転角度 θ に対する粒界エネルギーならびに，detBijα の平均や最小値，detBijα ＜0 となる
原子数を調べる．その後，Σ3 ねじれ・傾角粒界について，粒界に垂直方向の引張を分子動力学法により行い，粒界固有の強度について detBijα の視点から検討する．最後に第 5 章で本研究の総括を述べる．

3

第2章
解析手法の基礎2.1

分子動力学法

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

(2.1)

をたて，これを数値積分することにより粒子の軌跡を求める方法である．ここで，mα ，r α はそれぞれ原子 α の質量および位置ベクトルである．原子 α に作用する力 F α は，系のポテンシャルエネルギー Etot の各位置における空間勾配として次式により求められる．

Fα = −

∂Etot
∂r α

(2.2)

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

(2.5)

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

α

α

(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) と初速度から得ることができる．

2.2

原子間ポテンシャル

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

5

2.2.1

EAM ポテンシャル

EAM は金属中の多体効果を良好に再現することから広く用いられている．密度汎関数理論に基づき，まず金属材料における系のポテンシャルエネルギー Etot は原子を
価電子雲中に埋め込むエネルギーと原子間の 2 体間相互作用の和で与えられるとする．
さらに，埋め込みエネルギーは埋め込む位置の電子密度にのみ依存すると仮定することによって，系全体のエネルギーは次式のように表わされる．

Etot =N∑F (ρα ) +

α

N ∑
N
1∑
V (rαβ )2 α β(̸=α)(2.12)

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

(2.13)

β(̸=α)

で評価する．本解析では Mishin が Ni，Al に対し弾性定数，格子定数，原子エネルギー
等へのフィッティングを行った EAM ポテンシャルを使用する [35] ．ここでの二体間ポテンシャル V (rαβ )，二体間関数 ρα ，埋め込み関数 F (ρα ) はすべてスプラインデータで提案されており [36] ，数式としては与えられていない．本研究においてはこれらを 3 次のスプライン関数にフィッティングすることにより使用している．

6

2.3

速度スケーリング法

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

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

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

(

viα0 =

3kB T0mα)0.5(2.15)

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

viα =

3kB T
mα

)0.5

(2.16)

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

viα0
T0
=
α
vi
T

)0.5

(2.17)

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

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.18)mα

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

7

2.4

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

2.4.1

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

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

1
σij =
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

∂ 2F
∂ηij ∂ηkl)(2.19)と定義される [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.20)

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

8

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

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

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

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

ここで，V は平衡状態における系の体積であり，下付添字のローマ文字はテンソルのデカルト座標成分を表す．(2.21) 式の微分を求めるため，平衡状態からの仮想的な均一変形を考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によってr α = J r̄ α

(2.22)と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子 αと 原子 β の間の距離 rαβ には(rαβ )2 = r̄iαβ Gij r̄jαβ(2.23)なる関係が成立する．ただし，Gij = Jki Jkj である．仮想変形の Lagrange ひずみテン
ソル ηij は

ηij =

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

1dηij = dGij29(2.25)

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

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

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

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

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

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

σijα =) αβ αβ
N (
r r
1
1 ∑
F ′ (ρα )ρ′ (rαβ ) + ϕ′ (rαβ ) i αβj
V /N β(̸=α)
2r(2.27)と定義すると，系の応力は

σij =

N
1 ∑
σα
N α ij(2.28)となる．

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

1 ∂ 2 Etot
V ∂ηij ∂ηkl

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

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

′′

α

′

αβ

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

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

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

α

{ ∑
N

′

αβ

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

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

2.4.3

原子弾性剛性係数以上で定義した原子応力，弾性係数から，原子弾性剛性係数は以下で評価できる．
α
α
α
Bijkl
= Cijkl+ (σilα δjk + σjl
δik
α
α
+ σik
δjl + σjk
δil − 2σijα δkl )/2(2.33)α sym
ααWang らによる提案 [38] に従い，Voigt 対称性をもたせた Bijkl≡ (Bijkl+ Blkji)/2 を用いて安定性評価を行う．なお，ここでの i〜l は直交座標の指標 (1,2,3) もしくは (x，y，z) であるが，以降では xx，yy，zz，yz，zx，xy を 1〜6 とする Voigt 表記を用いて
Bijα と表すことにする．

11

第3章
薄膜・細線モデルによる解析
3.1

解析モデル及び解析条件

立方体セルの y，z 方向に周期境界を適用した薄膜モデル，ならびに z 方向のみ周期境界を適用した細線モデルを解析対象とする (図 3.1)．x，y，z 軸をそれぞれ fcc の結晶
方位 [100]，[010]，[001] とする (100) 表面モデル以外に，結晶方位を回転させて (110)や (111) 表面で構成されるモデルを Al，Ni それぞれについて 6 種類作成した．

(a)(100) 表面及び (100)-(010) ナノワイヤ
x 軸を [100] 方向，y 軸を [010] 方向，z 軸を [001] とする座標系において，単位格子を 20個ずつ並べた総原子数 32000 個のシミュレーションセルを用いる (図 3.2(a))．

(b)(110) 表面及び (110)-(110) ナノワイヤ
z 軸を中心に 45◦ 回転させ，x 軸を [110] 方向，y 軸を [110] 方向とする座標系において，(a) と同程度の大きさの立方体セルを切り出して用いる．総原子数は 31360 個である
(図 3.2(b))．
(c)(111) 表面及び (111)-(110) ナノワイヤ
x 軸を [111] 方向，y 軸を [110] 方向，z 軸を [112] とする座標系において，ほぼ同程度の形状のシミュレーションセルを切り出して用いる．総原子数は 32256 個である (図3.2(c))．

12

Ni，Al それぞれのセルの初期寸法を表 4.1 にまとめて示す．初期応力が 0 となるよ
うにセル寸法を微調整しながら 40000fs の初期緩和シミュレーションを行った．その後，z 方向に毎ステップひずみ増分 ∆εzz =1.0×10−6 を与えて引張シミュレーションを行う．この時，∆εzz =1.0×10−4 毎に原子座標データを記録し，各原子位置の原子応力α
σijα ならびに原子弾性係数 Cijkl
を求め，AES の行列式 detBijα (Voigt 表記) を評価した．

なお，薄膜モデルの際は横方向 (y 軸方向) の応力が 0 となるようにセル寸法をスケーリングしている．

Table 3.1 Analysis condition．
Coordination axis Number ofCell length[nm]x,y,zunit cellAlNi[100],[010],[001]
20×20×20 8.10×8.10×8.10 7.04×7.04×7.04[110],[110],[001]
28×28×20 8.02×8.02×8.10 6.97×6.97×7.04[111],[110],[112]12×28×16 8.42×8.02×7.94 7.32×6.97×6.90

13

Atomic number
32000
31360
32256

[001][010][100](a) (100) surface modelzyx(i) Thin plate model

[001]
[110]
[110]

(b) (110) surface model

[112]
[110]
[111]

(c) (111) surface modelzyx(ii) Nanowire model

Fig.3.2 Atom conﬁguration in the simulaFig.3.1 Schematic of simulation models．
tion cells．
14

3.2

Al の解析結果

3.2.1(100) 表面薄膜モデル(100) 表面薄膜のシミュレーション結果を図 3.3 にまとめて示す．0[K] 完全結晶バル
クの値で無次元化した detBijα の平均，および，最大・最小の変化を，それぞれ黒・赤・青線で左軸に目盛りをとって示している．なお，detBijα の標準偏差の範囲をエラーバー
で示しているので，平均の黒線はほぼその中心にあり色は見えづらくなっている．また図中には右軸に目盛をとって，灰色で応力-ひずみ曲線を示している．さらに下には detBijα ＜0 の原子数の変化をあわせて示している．応力-ひずみ曲線を見ると，応力εzz = 0.1665 で応力ピークを示し，その後応力が急減している．
detBijα の変化を見ると，シリコンや鉄での報告 [30], [32] と同様にエラーバーの幅が引張によって小さくなっており，引張による系の均一化を生じているものと考えられる．
図 3.3(b) は応力急減近傍の挙動を拡大したものである (エラーバーは表示していない)．
応力急減を生じる点でも detBijα の平均は負になっておらず，系全体の Bijα に対する不安定条件には達していない．また，応力ピークの εzz = 0.1665 の点においても detBijα が負の原子は存在しておらず，応力が急減しはじめてすぐの εzz = 0.1678 の時に detBijαが負の原子が発生する．すなわち Al の (100) 薄膜モデルでは，応力-ひずみのピークがdetBijα ＜0 の原子が現れる点より先にあるため，局所不安定の発生がトリガーになって系の応力急減を生じたのではない．応力-ひずみのピーク以降は detBijα の最大・最小が急激に拡大しているので均一化の限界ともとらえることができるが，応力急減す
なわち系の不安定が発生したためにエラーバーが拡大したともみなせ，現段階ではどちらが先か結論づけられない．図 3.4 に応力ピーク近傍の主要な点における原子配置
を，detBijα < 0 となった原子を緑色で着色して示す．先述のように，応力-ひずみの折れ曲がりであるひずみ (a)εzz = 0.1670 でも，detBijα ＜0 の原子は存在していない．折れ曲がり直後のひずみ (b)εzz = 0.1678 において detBijα ＜0 の原子が表面部分に現れ，その点を中心に表面だけではなく内部にも不安定原子が広がっている．

15

0.6

Maximum
Average

1.2

εzz=0.1665

0.8

0.4

Minimum000Number of detBαij<0 atoms

5

0.05

0.1

0.15

2000

εzz=0.1678

0
0

0.05

0.1

0.15

0.2

Strain, εzz

Stress-Strain curve
10

0.4

Maximum

εzz=0.1665

5

0.2

Average
Minimum

0
0.15

0.2

3000

1000

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

1.6

10

Number of detBαij<0 atoms

Stress-Strain curve

Stress, σzz GPa

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

2

0

0.155

0.16

0.165

0.17

0.155

0.16

0.165

0.17

0.175

0.18

3000
2000
1000

εzz=0.1678

0
0.15

0.175

0.18

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

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

(a) εzz=0.1670

[001]
[010]

(b) εzz=0.1678[100](c) εzz=0.1681(d) εzz=0.1683

Fig.3.4 Snapshot of (100) surface thin plate. Green circles indicate detBijα <0atoms（Al）．
16

3.2.2(110) 表面薄膜モデル(110) 表面薄膜のシミュレーション結果を図 3.5 にまとめて示す．εzz = 0.1237 で応力がわずかに低下して折れ曲がりを示した後，εzz = 0.1365 に再びピークを示して応力が急減した. (100) 表面薄膜とは異なり，応力-ひずみピークまでのエラーバーの幅はあまり小さくならない．また，(100) 表面薄膜よりも detBijα の平均が高い値で系の不安定挙動を生じている．図 3.5(b) は応力急減近傍の挙動を拡大したものである．detBijα＜0 となった原子数を見ると，εzz = 0.1237 の応力-ひずみの折れ曲がり直後に発生しているものの，やはり detBijα ＜0 の原子が発生して系の不安定を引き起こすというバ
ルクシリコンで観察されたストーリーとは異なる．図 3.6 に応力ピーク近傍の主要な点における原子配置を，detBijα < 0 となった原子を緑色で着色して示す．応力-ひずみの折れ曲がりであるひずみ (a)εzz = 0.1237 においては，detBijα ＜0 の原子は存在していない．折れ曲がり直後のひずみ (b)εzz = 0.1241 において，detBijα ＜0 の原子が現れ，表面の原子列に沿って拡大している．

17

10

Stress-Strain curve

1.2

Average
5

0.8

0.4

0

Minimum
0.05

0.1

0.15

0.2

3000
2000

εzz=0.1241

1000
0
0

0.05

0.1

0.15

1.5

0.2

Strain, εzz

10Stress-Strain curve1

εzz=0.1237
5

Maximum
0.5

Average00Minimum
0.105

Number of detBαij<0 atoms

00

Number of detBαij<0 atoms

εzz=0.1237

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

Maximum
1.6

Stress, σzz GPa

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

2

0.11

0.115

0.12

0.125

0.11

0.115

0.12

0.125

0.13

0.135

300
200
100

(a) Outline0

εzz=0.1241

0.105

0.13

0.135

Strain, εzz(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（Al）．

[001](a) εzz=0.1237[110](b) εzz=0.1241[110](c) εzz=0.1242

Fig.3.6 Snapshot of (110) surface thin plate. Green circles indicate detBijα <0atoms（Al）．

18

3.2.3(111) 表面薄膜モデル(111) 表面薄膜のシミュレーション結果を図 3.7 にまとめて示す．これまでの 2 つの薄膜と異なり，εzz = 0.1146 の応力ピークより前のひずみで detBijα ＜0 となる原子が現れている．応力急減近傍の挙動を拡大した図 3.7(b) を見ると，青線の detBijα の最小値がひずみ εzz = 0.1050 で前駆的な折れまがりを生じ，εzz = 0.1106 から一定数（1788個）の原子が detBijα ＜0 となる．detBijα ＜0 の原子はその後応力ピーク時に急増するが，このとき系の detBijα の平均も負となっており「系の不安定」条件に達している．
図 3.8 に応力ピーク近傍の主要な点における原子配置を，detBijα < 0 となった原子を緑色で着色して示した．左が detBijα が正の原子を点として表示し 3 次元的に見たもの，右の 2 つは正の原子も灰色の球で表示しセル断面を上から見たものである．前駆的に発生した detBijα ＜0 の原子は表面から 2 原子層目の原子であることが右図からわかる．その後応力ピークのひずみ εzz = 0.1147 において，表面ではなく内部から負の原子が増加し，系の detBijα の平均が負となる εzz = 0.1151 の図 3.8(d) ではセルを半分占めるように原子が detBijα < 0 となっている．
この (111) 表面薄膜では，detBijα ＜0 の原子が系の不安定変形より先に生じている
ため，さらにその不安定条件を明らかにするため 6×6 行列の Bijα の固有値解析を行った．図 3.9 は，図 3.8 で示した原子配置を，固有値 (λ1 ) の値で色付けしたものである．
先の detBijα ＜0 の図に対応し，固有値でも表面 2 原子層目が最初に負の値を示している (図 (a) の水色原子)．その後 (a)→(c) で系の内部の値が減少し全体に負値に近い状態となっている．ここで図 (c) の赤で囲った部分は図 3.8(c) で内部に detBijα ＜0 となった原子が現れた部分である（周期境界下で [110] 方向にシフトして表示)．図からわかるように固有値では赤で囲った部分と周囲の detBijα が正の部分には大きな差がなく，水色と緑色の λ のゆらぎで確率的に発生したものと推測される．なお 2 番目に小さな固有値 (λ2 ) を調べたが負となる原子は存在していないので，固有値 λ1 による評価とdetBijα の正値性による評価は本シミュレーション結果では一致する．

19

0.6
10

1.6 Maximum

Stress-Strain curve

5

0.8

0.4

Stress, σzz GPa

1.2

Average00Minimum
0

0.05

0.1

0.15

20000

εzz=0.1147

10000
0
0.05

0.1

0.15

εzz=0.1146(peak)
10

Stress-Strain curve
0.4

0.2

Strain, εzz(a) Outline

εzz=0.1151
5

0.2

Maximum

Average
0
0.1

0.2

30000

0

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

εzz=0.1146

Number of detBαij<0 atoms

Number of detBαij<0 atoms

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice20Minimum
0.105

0.11

0.115

0.12

0.125

0.13

0.115

0.12

0.125

0.13

30000
20000
10000
0
0.1

εzz=0.1147
0.105

0.11

Strain, εzz(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（Al）．

[112](a) εzz=0.1106

[110]

[112]

[110]

[111]

(c) εzz=0.1147(a’)εzz=0.1106(b) εzz=0.1146(d) εzz=0.1151[111](b’) εzz=0.1146

Fig.3.8 Snapshot of (111) surface thin plate. Green circles indicate detBijα <0atoms（Al）．

20

5.0

eigen value λ1

[112]

[110]
(b) εzz=0.1146(a) εzz=0.1106

[111]

-3.0

(c) εzz=0.1147(d) εzz=0.1151

Fig.3.9 Snapshot of (111) surface thin plate，colored by the minimum eigen
value of λ1 ．

21

3.2.4(100)-(010) 細線モデル(100)-(010) 細線モデルのシミュレーション結果を図 3.10 にまとめて示す．(100) 薄
膜よりも小さいひずみ εzz = 0.1381 で応力ピークを示し，その後原子配置の乱れを生じて応力が急減した．detBijα ＜0 の原子数変化を見ると，εzz = 0.1045 までは引張前の値 80 で一定で，その後から減少し εzz = 0.1225 で一度 0 となる．その後再び増加す
るが，図 3.10(b) の拡大図を見ると εzz = 0.1381 の応力ピークにおいて detBijα の最大値・最小値は急激に変化しているのに対し detBijα < 0 の増加は連続的かつなだらかである．すなわち系の不安定変形時には detBijα < 0 の原子数増加よりも detBijα の固有値の変化が重要になっている可能性がある．また，他の系では応力ピーク後に detBijα の最大・最小値が急変しているのに対し，この (100)-(010) 細線のみピーク前に急変しており，内部の構造変化がピークをもたらしている可能性がある．
図 3.11 に応力ピーク前後における detBijα < 0 の原子配置を示す．先述の，εzz =0.1045 までの一定数の負の原子はエッジ部分に存在する．負の原子は一度消滅し，図(b) で示したように εzz = 0.1226 で再びエッジ部に出現する．(c)→(d) で負の原子がエッジ部に増加し，応力ピークの (d)εzz = 0.1381 では矢印で示したように系内部に伝播している．
図 3.12，図 3.13 は図 3.11 に対応する固有値 λ1 ，λ2 の値の分布である．detBijα の正値条件では一度エッジ部の負の原子は消滅していたが，固有値で見るとエッジ部の λ1は負値を示している（図 (b)）．同様に図 3.13 の λ2 の値の分布を見ると，こちらもエッジ部の λ2 は負値を示している（図 (b)）．このように，detBijα < 0 の条件では見えない
変化が固有値 λ の値で分かることがあるが，多くの場合は detBijα < 0 の条件と λ < 0の条件は一致している．

22

Stress-Strain curve

6

1.2

Average

0.8

4

2

0.4

εzz=0.138100Minimum
0

0.05

0.1

0.15

0.2

3000
2000
1000
0
0

0.05

0.1

0.15

0.2

Strain, εzz

Stress-Strain curve
1.5

1

10

Maximum
5

0.5

0
0.12

Average

εzz=0.1381

Minimum

0.125

0.13

Stress, σzz GPa

8

Number of detBαij<0 atoms

1.6

Stress, σzz GPa

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice
Number of detBαij<0 atoms

Maximum

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

10

2

0
0.135

0.14

0.145

0.15

0.135

0.14

0.145

0.15

1500
1000
500
0
0.12

(a) Outline

εzz=0.1225
0.125

0.13

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

Fig.3.10 Change in average，standard deviation of detBijα ，and stress-strain
curve of (100)-(010) surface nanowire（Al）．

[001](a) εzz=0.1045[010](b) εzz=0.1226[100](c) εzz=0.1330(d) εzz=0.1381

Fig.3.11 Snapshot of the (100)-(010) surface nanowire．Green circles indicateddetBijα < 0 atoms（Al）．

23

eigen value λ1

25

[001]
(a) εzz=0.1045[010](b) εzz=0.1226

[100]

-7

(c) εzz=0.1330(d) εzz=0.1381

Fig.3.12 Snapshot of the (100)-(010) surface nanowire，colored by the minimum eigen value of λ1 ．

eigen value λ2

30

[001]
(a) εzz=0.1045[010](b) εzz=0.1226

[100]

-5

(c) εzz=0.1330(d) εzz=0.1381

Fig.3.13 Snapshot of the (100)-(010) surface nanowire，colored by the 2nd
minimum eigen value of λ2 ．

24

3.2.5(110)-(1̄10) 細線モデル(110)-(110) 細線モデルのシミュレーション結果を図 3.14 にまとめて示す．やはり薄膜の場合より小さなひずみ εzz = 0.0938 で応力が頭打ちになり，この点が系の弾性限界と考えられる．一方，detBijα の最小値変化をみると εzz = 0.0521 とピークひずみよりずっと前のひずみで急減し負の値を示している．一方，detBijα ＜0 となった原子数の変化では，εzz = 0.0691 までは揺らいではいるがほぼ一定数（30 個）となっており,detBijαの最小値の折れ曲がりとは対応していない．ひずみ εzz = 0.0691 から detBijα が負の原子は増加しているが，急激に増えるのではなく，ゆるやかに増加している．図 (b) の拡
大図をみると，応力-ひずみの折れ曲がり点から detBijα < 0 の原子数が急増している．

detBijα の最大・最小が急変するのは変形が発生する破線より後である．
図 3.15 に応力ピーク前後における detBijα < 0 の原子配置を示す．detBijα の最小値が前駆的に急減する点の前後の図 (a)，(b) を見ると負の原子は対角線方向の 2 つのエッジ部に生じており，(b)εzz = 0.0522 では矢印に示したようにエッジ部横の 1 原子も負となっている．応力ピークの (d) では別の対角線方向のエッジ部付近に負の原子が発生しており，最初に負の原子が存在していたエッジではない所から変形が生じていたことがわかる．
固有値 λ1 の分布を図 3.16 に示す．detBijα < 0 による評価と大きな差はない．また，

λ2 が負の原子はほとんど存在しなかった．

25

Maximum1.6

8

Stress-Strain curve

εzz=0.0938

1.2

6

0.8

4

0.4

2

Average
0

0

1.5

10

Maximum

εzz=0.0938
Stress-Strain curve
15
0.5

Average

0

Stress, σzz GPa

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

10

Stress, σzz GPa

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice20Minimum

0.05

0.1

0.15

0.2

Number of detBαij<0 atoms

Number of detBαij<0 atoms

Minimum
0

0.075

0.08

0.085

0.08

0.085

0.09

0.095

0.1

0.105

0.09

0.095

0.1

0.105

1500

3000

1000

2000
1000
0
0

0.05

0.1

0.15

0.2

500
0
0.075

Strain, εzz

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

Fig.3.14 Change in average，standard deviation of detBijα ，and stress-strain
curve of (110)-(110) surface nanowire（Al）．

(a) εzz=0.0521

[001]
[110]

(b) εzz=0.0522[110](c) εzz=0.0691(d) εzz=0.0938

Fig.3.15 Snapshot of the (110)-(110) surface nanowire．Green circles indicateddetBijα < 0 atoms（Al）．

26

eigen value λ130(a) εzz=0.0521

[001]
[110]

(b) εzz=0.0522

[110]

-15

(d) εzz=0.0938(c) εzz=0.0691

Fig.3.16 Snapshot of the (110)-(110) surface nanowire，colored by the minimum eigen value of λ1 ．

27

(111)-(1̄10) 細線モデル

3.2.6(111)-(110) 細線モデルのシミュレーション結果を図 3.17 にまとめて示す．他のナノワ
イヤ同様，detBijα の最小値は初期から負に近い値でほとんど変化なく，また εzz = 0.0919の応力ピークよりも前に負値を示している．detBijα ＜0 となった原子数の変化を見ると，他のナノワイヤと異なり引張前に detBijα ＜0 の原子は存在せず，εzz = 0.0515 で 0→64，

εzz = 0.0596 で 64→128 と段階的に増加している．その後応力ピークの ∆εzz = 0.0515
前の εzz = 0.0912 から detBijα ＜0 の原子が急増している．
図 3.18 にいくつかのひずみにおける detBijα < 0 の原子配置を示す．(a)εzz = 0.0515
における 64 個の負の原子はエッジ部第一層で，(b)εzz = 0.0596 において 2 原子層目も負になり 128 個となる．一方，detBijα ＜0 の原子がさらに増加する (c)εzz = 0.0912 ではエッジ部第三層目に負の原子が発生しはじめ，応力ピークの (d)εzz = 0.0919 においてはエッジ部 3 原子層が全て負となっている．
図 3.19 に，λ1 の分布を示す．(110)-(110) 細線と同様に detBijα < 0 による評価と大きな差はなかっため，以降では固有値による解析は省略する．
10

1.2

6

0.8

4

2

0.4

Average
εzz=0.0919

00Minimum
0

0.05

0.1

0.15

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect latticeStress-Strain curve8

Stress, σzz GPa

1.6

0.2

Number of detBαij<0 atoms

Number of detBαij<0 atoms

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

Maximum

3000
2000

εzz=0.0596

1000 ε =0.0515
zz
0
0

0.05

0.1

0.15

0.2

εzz=0.0919 (peak)
1.5

10Stress-Strain curve
1
5

Maximum
0.5

Stress, σzz GPa

2

Average00.07

0Minimum
0.075

0.08

0.085

0.09

0.095

0.1

0.095

0.1

1500
1000

εzz=0.0912

500
0
0.07

0.075

0.08

0.085

0.09

Strain, εzz

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

Fig.3.17 Change in average，standard deviation of detBijα ，and stress-strain
curve of (111)-(110) surface nanowire（Al）．

28

(a) εzz=0.0515

[112]

[111]

(b) εzz=0.0596[110](d) εzz=0.0919(c) εzz=0.0912

Fig.3.18 Snapshot of the (111)-(110) surface nanowire．Green coloured circles
indicated detBijα < 0 atoms（Al）．

25

eigen value λ1

[112]

[111]

(a) εzz=0.0515(b) εzz=0.0596

[110]

-25

(c) εzz=0.0912(d) εzz=0.0919

Fig.3.19 Snapshot of the (111)-(110) surface nanowire，colored by the minimum eigen value of λ1 ．

29

3.3

Ni の解析結果

3.3.1(100) 表面薄膜モデル

Ni の (100) 表面薄膜のシミュレーション結果を図 3.20 に Al の場合と同様にまとめて示す．εzz = 0.1165 で応力ピークを示すが，Al のように折れ曲がりを示すことなく滑らかな曲線を伴ったまま急減した．応力急減近傍の挙動を拡大した図 3.20(b) を見ると，応力ピークを示す点の直前に detBijα の平均が負に達しており，本 (100) 薄膜においては「系全体の不安定条件」に達して応力急減したと考えることができる．さらに，

detBijα ＜0 となった原子数を見ると，無負荷平衡状態では detBijα が負の原子数は存在せず，εzz = 0.0907 において 1600 個の原子が負となる．また，きわめて短時間の挙動ではあるが，εzz = 0.1160 から detBijα ＜0 の原子数が急増し，応力ピークの直前に一気に増加してすべての原子が負となっている．したがって本 (100) 薄膜ではごく短時間の間に，局所不安定の発生 ⇒ 系の不安定条件 ⇒ 応力ピークという現象を生じているものと捉えることができる．
応力急減時の detBijα < 0 の原子配置を図 3.21 に示す．1600 個の負の原子は表面第一層に存在する（図（a））．εzz = 0.1160 からの変化は (b)εzz = 0.1160，(c)εzz = 0.1161のように表面から内部へ一原子層ずつ広がり，(d)εzz = 0.1163 ではすべての原子がdetBijα ＜0 となっている．

30

1
10

εzz=0.1165
0.5

5

Average00Minimum
0.05

0.1

0.15

0.2

30000
20000

εzz=0.0907

10000
0
0

0.05

0.1

0.15

Stress, σzz GPa

20

15

0

Number of detBαij<0 atoms

Maximum

Stress-Strain curve

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

1.5

0.6

0.2

Strain, εzz

25

εzz=0.1165 (peak)
20
0.4
15

Stress-Strain curve
10

Maximum
0.2

5Stress, σzz GPa

25

Average00Minimum
0.1

Number of detBαij<0 atoms

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

2

0.105

0.11

0.115

0.12

0.125

-5
0.13

0.12

0.125

0.13

30000
20000
10000
0
0.1

(a) Outline

εzz=0.1160
0.105

0.11

0.115

(b)Zoom-up around stress drop

Fig.3.20 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (100) surface（Ni）．

(a) εzz=0.1159

[001]
[010]

(b) εzz=0.1160[100](c) εzz=0.1161(d) εzz=0.1163

Fig.3.21 Snapshot of (100) thin plate surface. Green circles indicate detBαij < 0atoms（Ni）．

31

3.3.2(110) 表面薄膜モデル

Ni の (110) 表面薄膜のシミュレーション結果を図 3.22 にまとめて示す．本モデルでは εzz = 0.1350 の応力ピークは Al 同様急激な折れ曲がりを示している．detBijα の最大値・最小値ともにピークひずみよりずっと前で折れ曲がりを示しており，特に detBijαの最小値は εzz = 0.1042 以降，負値を示している．detBijα ＜0 の原子数変化を見ると，

εzz = 0.1042 までは 0 で，その後 2240 個（εzz = 0.1042）→3360 個（εzz = 0.1216）→
4480 個（εzz = 0.1257）と段階的に増加している．応力ピークより前の εzz = 0.1263 で急激に増加し，このとき detBijα も平均が負となっている（図（b））．detBijα の最大値が急変するのは応力ピーク以降だが，先述のように detBijα の平均はピークよりも前のひずみ εzz = 0.1263 で負になっている．系の不安定に達していることから，変形速度の関係でピークがオーバーシュートした可能性がある．
図 3.23 に応力ピーク前後における detBijα < 0 の原子配置を示す．図 3.23 を見ると，負の原子が発生する (a)εzz = 0.1042 においては表面 1 原子層が負となり，(b)εzz = 0.1044
では 2 原子層と増えていく．負の原子数が急激に増加する前の (c)εzz = 0.1262 では表面 4 原子層分が負となっている．(d)εzz = 0.1263 から detBijα が負の原子数の増加が内部で始まり，応力ピーク (f) の時，表面のみではなく系が全体的に負となっている．

32

0.6

20

1.5

15
1

εzz=0.1350

10Maximum

0.5

5

Average
00Stress, σzz GPa

Stress-Strain curve

20
0.4

0.1

0.15

15

Maximum

10
5

Average00Minimum
0.12

0.2

Number of detBαij<0 atoms

Number of detBαij<0 atoms

0.05

Stress-Strain curve

0.2

Minimum025

εzz=0.1350 (peak)Stress, σzz GPa

25

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

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

2

0.125

0.13

0.135

0.14

0.145

-5
0.15

0.14

0.145

0.15

30000 εzz=0.1257

30000
20000

20000

εzz=0.1216

εzz=0.1262

10000

10000
0
0

0.05

0.1

0.15

0.2

0
0.12

0.125

0.13

0.135

Strain, εzz

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

Fig.3.22 Change in the average，standard deviation of detBijα ，and stressstrain curve of thin plate of (110) surface（Ni）．

[001](a) εzz=0.1042[110](b) εzz=0.1044(c) εzz=0.1262(e) εzz=0.1264(f) εzz=0.1350[110](d) εzz=0.1263

Fig.3.23 Snapshot of (110) thin plate surface. Green circles indicate detBαij <0atoms（Ni）．
33

3.3.3(111) 表面薄膜モデル

Ni の (111) 表面薄膜のシミュレーション結果を図 3.24 および図 3.25 にまとめて示す．Ni では (100) 表面薄膜，(110) 表面薄膜とも応力ピーク前に detBijα の平均が負となっていたが，(111) 表面薄膜では応力ピークの εzz = 0.1191 でも detBijα の平均は，わずかではあるが正値を示している．detBijα ＜0 となった原子数の変化も，破線で示した構造変化後に急増している．図 3.25 の detBijα < 0 の原子配置を見ると，不安定変形を生じる前の εzz = 0.1161 においては Al と同様に表面第 2 層が負値をとっている．
応力ピークの εzz = 0.1191 でも変化はなく，その後の構造変化の発生によって内部に

detBijα < 0 の原子が急増する（図 (c)，(d)）．したがって Ni の他の薄膜と違い (111) 薄膜では局所の detBijα < 0 がトリガーになったり，または系全体が detBijα < 0 の条件に達したのではなく，何が応力の頭打ちをもたらしているのか不明である．

34

10

εzz=0.1191
0.5

5

Average
Minimum

0

0.05

0
0.1

0.15

25
20

0.4

10

εzz=0.1191 (peak)

0.2

Maximum

5

Average
0 Minimum

0.1

0.2

15

Stress-Strain curve

0.105

Stress, σzz GPa

1

Stress, σzz GPa

15

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

Maximum

Number of detBαij<0 atoms

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

20

Stress-Strain curve

1.50

Number of detBαij<0 atoms

0.6

25

2

0

0.11

0.115

0.12

0.125

-5
0.13

30000

30000

20000

20000

εzz=0.1161

10000
0
0

0.05

0.1

εzz=0.1194

10000

0.15

0.2

Strain, εzz

0
0.1

0.105

0.11

0.115

0.12

0.125

0.13

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

Fig.3.24 Change in average，standard deviation of detBijα ，and stress-strain
curve of thin plate of (111) surface（Ni）．

[112]

[110]
(b) εzz=0.1191(a) εzz=0.1161[111](c) εzz=0.1194(d) εzz=0.1195

Fig.3.25 Snapshot of (111) thin plate surface. Green circles indicate detBαij <0atoms（Ni）．

35

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

Ni の (100)-(010) 細線モデルのシミュレーション結果を図 3.26 および図 3.27 にまとめて示す．ただし図 3.27 は表面部の原子により見えにくいため細線断面を上から見た図として示している．detBijα ＜0 の原子数変化を見ると，εzz = 0.0502 から detBijα が負の原子が現れるが，図 3.27(a) に示したようにエッジ部の原子列が負となっている．
その後図 3.27(b)，(c) のようにエッジ部から表面全体を覆うように detBijα < 0 の原子が増加する．ここで (b)εzz = 0.0970 は図 3.27(a) において detBijα の最小値 (青線) が折れ曲がる点に対応する．応力ピークの点 (εzz = 0.0183) ではエッジ部から内部に変形を生じていることが図 3.27(d) から推測されるが，この (a)→(d) における detBijα < 0の変化は図 3.26(b) 下図からわかるように連続で，detBijα < 0 の数による線引きはできない．ただし，detBijα の最大・最小は εzz = 0.0183 の破線よりも前に急激に拡大しているので，わずかにではあるが応力ピークよりも先に構造変化=系の不安定挙動が発生していた可能性がある．

36

1
10

εzz=0.1183
0.5

Stress, σzz GPa

15

5

Average00Minimum
0

Number of detBαij<0 atoms

20

Maximum

0.05

0.1

0.15

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect latticeStress-Strain curve
1.5

1.2

0.2

25

Stress-Strain curve
20
0.8
15

εzz=0.1183

Maximum

10

0.4
5

Average
0

0.1

Stress, σzz GPa

25

Number of detBαij<0 atoms

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice20Minimum
0.105

0.11

0.105

0.11

0.115

0.12

0.125

-5
0.13

0.115

0.12

0.125

0.13

10000

30000
20000
10000
0
0

0.05

0.1

0.15

0.2

5000

0
0.1

Strain, εzz

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

Fig.3.26 Change in average，standard deviation of detBijα ，and stress-strain
curve of (100)-(010) surface nanowire（Ni）．

(a) εzz=0.0502

[010]

[001]

(b) εzz=0.0970[100](c) εzz=0.1130(d) εzz=0.1183

Fig.3.27 Snapshot of the (100)-(010) surface nanowire．Green circles indicate
detBijα < 0 atoms（Ni，top view）．

37

3.3.5(110)-(1̄10) 細線モデル

Ni の (110)-(110) 細線モデルのシミュレーション結果を図 3.28 に示す．応力ピーク
は εzz = 0.1267 であり，detBijα の平均はほぼ 0 であるが負値をとるのは応力ピークのやや後になる．detBijα ＜0 となった原子数は，εzz = 0.0754 まではほぼ 0 である．その後少しづつ増えていき，detBijα の最小値が折れ曲がる点 εzz = 0.0782 においては 36 とそれほど大きな変化はなかった．その後 εzz = 0.1007 までは 200 個の原子が負となってしばらく一定となる．さらに εzz = 0.1220 までに 4100 と階段状に増え，応力ピーク
より前の εzz = 0.1221 からゆるやかに増加し，ピーク後急増する（図 3.28(b) 下）．
図 3.29 に応力ピーク前後における detBijα < 0 の原子配置について断面を上から見たものを示す．detBijα の最小値が折れ曲がる（最小値が負値を示す）点 (a) ではエッジ
線状に原子列が存在するが，エッジのその原子列ではなく 1 つ内側の 2 原子列が負となっている．εzz = 0.1007 から 0.1220 までに 4100 個の原子が負になるときには図 (c)のように表面 2 原子層分が負となっている．応力ピーク点の (d) においてはさらに内側の 1 原子層分が負となっている．

38

25

εzz=0.1267

1.2

20

1.5

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect latticeStress-Strain curve

20

1
10

0.5

5

Average00Maximum

0.15

εzz=0.1008

10000

εzz=0.0782

0
0

0.05

0.1

0.15

0

0

0.11

30000
20000

5

Average

Minimum

0.2

Number of detBαij<0 atoms

Number of detBαij<0 atoms

0.1

10

0.4

Minimum
0.05

15

εzz=0.1267(peak)Stress, σzz GPa

15

0

25

Stress-Strain curve

0.8

Maximum

Stress, σzz GPa

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

2

0.2

0.115

0.12

0.125

0.13

0.135

-5
0.14

0.125

0.13

0.135

0.14

20000

εzz=0.1220

10000

0
0.11

0.115

0.12

Strain, εzz

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

Fig.3.28 Change in average，standard deviation of detBijα ，and stress-strain
curve of (110)-(110) surface nanowire（Ni）．

[110](a) εzz=0.0782(b) εzz=0.1007
[001]

[110]

(c) εzz=0.1220(d) εzz=0.1267

Fig.3.29 Snapshot of the (110)-(110) surface nanowire．Green circles indicate
detBijα < 0 atoms（Ni，top view）．

39

3.3.6(111)-(1̄10) 細線モデル

Ni の (111)-(110) 細線モデルのシミュレーション結果を図 3.30 に示す．他のナノワイ
ヤ同様，detBijα の最小値は引張初期から 0 に近い値でほぼ一定を示した後，応力ピー
ク (εzz = 0.1058) よりも前のひずみ εzz = 0.0885 で負値へと急減している．detBijα の最大値が急上昇する点（構造変化による回復）は応力ピークにほぼ一致するがわずかにピーク直後であり，detBijα の平均は負になっていない．
図 3.31 に応力ピーク前後における detBijα < 0 の分布を断面で示す．detBijα の最小値が負を示す (a) において負の原子はエッジ部に発生している．その後，detBijα ＜0 の原子は図 3.30(b) 下に拡大して示したように凹凸を示しながら増加するが，図 (b)，(c)に示したようにエッジ部付近に拡大している．その後 (110) 表面に負の原子が拡大し，応力ピーク (d) では表面原子が負となっている．また，その後 detBijα ＜0 となった原子数が応力ピーク後に極小を迎える (e) においては，[110] 表面から系内部へと，負の原子が発生しており，変形が内部に伝播していることが示唆される．

40

15

1
10

0.5

5

Average00Minimum
0

0.05

0.1

0.15

0.2

30000
20000
10000

εzz=0.0885

0
0

0.05

0.1

Strain, εzz

0.15

Stress, σzz GPa

Maximum

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

20

1.5

Number of detBαij<0 atoms

1.2

εzz=0.1058

Stress-Strain curve

0.2

25

Stress-Strain curve
20
0.8
15

εzz=0.1058
10
0.4

Maximum
5

Stress, σzz GPa

25

Average00Minimum
0.09

Number of detBαij<0 atoms

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

2

0.095

0.1

0.105

0.11

0.115

-5
0.12

0.105

0.11

0.115

0.12

5000

2500

εzz=0.0924

0
0.09

0.095

0.1

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

Fig.3.30 Change in average，standard deviation of detBijα ，and stress-strain
curve of (111)-(110) surface nanowire（Ni）．

(a) εzz=0.0885(b) εzz=0.0924(c) εzz=0.0925(d) εzz=0.1058(e) εzz=0.1086

[110]

[112]

[111]

Fig.3.31 Snapshot of the (111)-(110) surface nanowire．Green coloured circles
indicate detBijα < 0 atoms（Ni，top view）．
41

3.4

結言

fcc 金属の Al と Ni について，引張時の系の不安定挙動=応力急減におよぼす表面不均一性の影響を調べるため，(100)，(110)，(111) 表面を有する薄膜や細線の引張シミュ
レーションを行い，挙動を調べるとともに原子弾性剛性係数（detBijα ）の正値性の観点から検討した．得られた結果を以下に示す．

(1) Al(111) の表面薄膜，Ni の全薄膜 ((100)，(110)，(111)) と (110) 表面細線は，応力-ひずみのピーク点付近で detBijα の平均が負となっていることから「系の不安定条件」に達したため系の不安定挙動が発生したものと考えられる．

(2) それ以外の系でも Ni は全体的に detBijα の平均は 0 に近い値をとっているが，Alの他の系は detBijα の平均が比較的大きな正の値のとき応力ピークを生じており，系全体の不安定に達する前に応力急減している．

(3) Al の (100) 薄膜，および (110) 薄膜以外では，応力ピークよりかなり前に detBijαの最小値が減少して負値をとり，表面やエッジ部に負の原子が発生・拡大している．

(4) Al の (100) 薄膜，(110) 薄膜は応力ピーク後に負の原子が発生していること，かつ応力ピーク時にも detBijα の平均が大きな正値をとっていることから，現時点では応力ピークをもたらすメカニズムを detBijα の観点から理解することは難しい．

42

第4章
粒界の原子弾性剛性係数評価
前章では，表面による構造不均一性を有する系の不安定挙動を detBijα の正値性に基づいて議論することを検討した．本章では粒界による構造不均一性について，同様に

detBijα に基づいて検討する．

4.1

無負荷平衡時の粒界エネルギーと原子弾性剛性係数

4.1.1

解析条件

ねじれ粒界モデル
図 4.1 に示すように，x 軸を [110] 方向，y 軸を [001] 方向，z 軸を [110] 方向とする座標系において，z 軸中心に ± θ/2 回転させて接合し対称ねじれ粒界を作成する．粒界間の距離が 5[nm] 以上となるようにセルの z 方向寸法をとっている．Σ 値が 1∼80 のものを解析対象とし，全方向周期境界条件下において，初期応力が 0 となるようにセル寸法を微調整しながら 0.1[K] で 100000fs の初期緩和シミュレーションを行った後，粒界エネルギーならびに detBijα を評価した．

yxz

Fig.4.1 Schematic of simulation cell for twist grain boundary．

43

傾角粒界モデル
図 4.2 に模式的に示すように，x 軸を [110] 方向，y 軸を [001] 方向，z 軸を [110] 方向にとり，z 軸を回転軸として鏡面対称となるように，それぞれ角度 θ/2 回転させることで対称傾角粒界を作成する．なお，粒界近傍で原子が完全結晶の第一近接距離の 0.8倍より接近する場合，片方を取り除き，もう片方をその組み合わせの中心（つまり粒界面）に配置して作成している (図 4.3 に一例を示す)．ここでも y 軸方向の周期境界について，粒界間の距離が 5[nm] 以上になるようにセル寸法をとっている．Σ 値が 1∼80のものを解析対象としてねじれ粒界と同じ条件で初期緩和を行い，粒界エネルギーと

detBijα の評価を行った．

x

z

y

Fig.4.2 Schematic of simulation cell for tilt grain boundary．

(a) Initial placement．

(b)Replacement．

(c) After Relaxation.

Fig.4.3 Remove process for overlapping atoms on tilt grain boundary．

44

4.1.2

解析結果および考察

ねじれ粒界の結果
図 4.4 に Al の対称ねじれ粒界の粒界エネルギーを，(a)detBijα を系全体で平均した
値，(b)detBijα の合計から同じ原子数の完全結晶の detBijα の合計を引き，粒界面の面積で割ったもの，(c)detBijα の最小値，(d)detBijα ＜0 の原子数と比較してそれぞれ示した．よく知られているように粒界エネルギーは Σ3 で低いエネルギーカスプを示している．図 (a) の detBijα の平均は，ねじれ角に対する依存性を見出せないが，detBijα の平均はバルク部分の原子が多くなるにつれ粒界部分の影響が小さくなる．そこで (b)で差をとり，かつ単位面積当たりの値で比較すると Σ 値による差は明確になっているが，粒界エネルギーの分布形態とは必ずしも対応していないことがわかる．図 (c) のdetBijα の最小値を見るとほとんどの粒界で負の原子が存在しているが，もっとも粒界エネルギーが小さい Σ3 では最小値が正値であり，detBijα ＜0 の原子が存在しない構造であることがわかる．図 (d) の detBijα ＜0 の原子分布では前述のように Σ3 で 0 となっている．また，図に矢印で示している低 Σ 値の粒界は負の原子数が比較的小さくグラ
フ上で極小点となっており，安定な構造を有していることが予想される．
図 4.5 に Ni の対称ねじれ粒界の解析結果をまとめて示す．Al と同様に粒界エネル
ギーは Σ3 でエネルギーカスプを示し，それ以外の角度では大きなエネルギー極小点は存在せず上に凸の曲線分布を示している．図 (a)，(b) の detBijα の平均ならびに粒界部分の寄与のどちらを見てもほぼ同じ値を示しておりねじれ角による変化は見られな
い．図 (c) の detBijα の最小値は角度によって変化があり，Al と同様に粒界エネルギーが低い Σ3 は 0 に近いが正値である．また Al の場合はねじれ角が |θ| ≦ 22.8◦ の範囲
では大きな負値を示し，他はほぼ 0 近傍の値をとっているのに対し，Ni も同様の傾向は見られるもののその差はあまり大きくない．Ni の場合はむしろ図 (d) の detBijα ＜0の数からわかるように小さなねじれ角では多くの原子が負となっていることがわかる．
低 Σ の粒界で detBijα ＜0 の数が比較的少なくなるのは Al と同じだが，Σ19 については隣の Σ57 がほぼ 0 になった．

45

0.5

0.5

0

0
0

30

60

90

120

150

Σ19

1

0.5

0.5

0

0

0

Σ11 Σ3 Σ17 Σ3

1

1

0.5

0.5

0

0

-0.5

-0.5
180

0

30

60

90

120

150

Σ19
Boundary Energy [J/m2]

Σ19

60

90

120

150

Boundary Energy
Number of detBαij<0 atoms

1.5

1.5

Minimum of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice

Boundary Energy [J/m 2 ]

1.5

30

-0.5
180

Misorientation angle, θ [degrees](b) Boundary energy and detBαij on GB
at Al in Twist model.

Misorientation angle, θ [degrees](a) Boundary energy andaverage of detBαij at Al in Twist model.
Boundary Energy
Minimum of detBαij

Σ11 Σ3 Σ17 Σ3

1

-0.5

180

1.5

Difference ∆AES
(ΣdetBijGB-ΣdetBijBulk)/2A[ /m 2 ]1

1

Boundary Energy [J/m 2 ]Σ11 Σ3 Σ17 Σ3

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

Boundary Energy [J/m 2 ]Σ19Boundary Energy
detBαij on GB

1.5

1.5

300

Σ11 Σ3 Σ17 Σ3

1

200

0.5

100

0

0
0

30

60

90

120

150

180

Misorientation angle, θ [degrees](d) Boundary energy and Number of detBαij<0 atoms
at Al in Twist model.

Misorientation angle, θ [degrees](c) Boundary energy and Minimumof detBαij at Al in Twist model.

Fig.4.4 Simulation results of twist grain boundary（Al）．

46

Number of detBαij<0 atoms

Boundary Energy
average of detBαij

1.5

2

2

1

1

0

0
0

30

60

90

120

150

Σ19

2

1

1

0

0

0

Σ11 Σ3 Σ17 Σ3

2

2

1

1

0

0

-1

0

30

60

90

120

150

Σ19
Boundary Energy [J/m2]

Σ19

60

90

120

150

Boundary Energy
Number of detBαij<0 atoms

3

3

Minimum of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice

Boundary Energy [J/m 2 ]

3

30

-1
180

Misorientation angle, θ [degrees](b) Boundary energy and detBαij on GB
at Ni in Twist model.

Misorientation angle, θ [degrees](a) Boundary energy and average of
detBαij at Ni in Twist model.
Boundary Energy
Minimum of detBαij

Σ11 Σ3 Σ17 Σ3

2

-1

180

3

Difference ∆AES
(ΣdetBijGB-ΣdetBijBulk)/2A[ /m 2 ]Σ11 Σ3 Σ17 Σ3

Boundary Energy [J/m 2]

Σ19

3

3

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

Boundary Energy [J/m 2 ]3

Boundary Energy
detBαij on GB

-1
180

300

Σ11 Σ3 Σ17 Σ3

2

200

1

100

0

0
0

30

60

90

120

150

180

Misorientation angle, θ [degrees](d) Boundary energy and Number of detBαij<0 atoms
at Ni in Twist model.

Misorientation angle, θ [degrees](c) Boundary energy and Minimumof detBαij at Ni in Twist model.

Fig.4.5 Simulation results of twist grain boundary（Ni）．

47

Number of detBαij<0 atoms

Boundary Energy
average of detBαij

傾角粒界の結果
図 4.6 に Al の対称傾角粒界の解析結果をまとめて示す．⟨110⟩ 対称傾角粒界の粒界エ
ネルギーは先の ⟨110⟩ 対称ねじれ粒界よりもエネルギーカスプが強く現れ，(111)Σ3，(311)Σ11 粒界の粒界エネルギーが非常に低い値となっている．傾角粒界モデルでは Σ
値によってセルの大きさが変わり，バルク部分の割合がまちまちなので (a) の detBijαの平均は参考までに示している．図 (b) の粒界近傍の detBijα を抽出した値も，θ に比例して右下がりになっているように見える以外は特筆すべき点は見受けられない．一方，
図 (c) の detBijα の最小値を見ると粒界エネルギーに対応した分布を示している．粒界エネルギーが非常に低い (111)Σ3，(311)Σ11 粒界では detBijα の最小値は大きな正値をとっている．図 (d) の detBijα ＜0 の原子数分布を見ると，傾角粒界では負の原子が 0の粒界が多く存在する．
図 4.7 に Ni の対称傾角粒界の解析結果をまとめて示す．(111)Σ3 の粒界エネルギーが極めて小さく，バルクとの差がほとんどない．図 (b) の粒界の detBijα の平均も (111)Σ3
ではバルクとの差がほとんどない．それでも，図 (c) の detBijα の最小値は，Al と同様に (111)Σ3 で高い値を示しているものの，バルクの半分程度の値となっており力学的には差があることがわかる．図 (d) の detBijα ＜0 の原子数のグラフを見ると，ねじれ角 θ が小さいときに負の原子が多く存在する傾向があり，117.5◦ ≦ |θ| ≦ 163.9◦ では
負の原子は存在しない．

48

(332)(111)(334) (112)(311)Σ19

Σ11 Σ3 Σ17 Σ3 Σ11
1.5

0.6
1
0.3
0.5

0

0
0

30

60

90

120

150

0.6

1

0.3

0.5

0

0

-0.3

180

0

1.5
0.6
1
0.3
0.5

0
0

30

60

90

120

150

90

120

150

-0.5
180

100

(331) (332) (111)(334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ11

Boundary Energy [J/m 2 ]Boundary Energy [J/m2](331) (332) (111) (334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ110

60

Boundary Energy
Number of detBαij<0 atoms

0.92

Minimum of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice

0.9

30

Misorientation angle, θ [degrees](b) Boundary energy anddetBαij on GB
at Al in Tilt model.

Misorientation angle, θ [degrees](a) Boundary energy andaverage of detBαij at Al in Tilt model.

Boundary Energy
Minimum of detBαij

1.5(331) (332) (111)(334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ11

Boundary Energy [J/m 2]

(331)

0.9

Difference ∆AES
(ΣdetBijGB-ΣdetBijBulk)/2A[J/m 2 ]2

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

Boundary Energy [J/m 2 ]0.9

Boundary Energy
detBαij on GB

0.6
50
0.3

0

0
0

180

Number of detBαij<0 atoms

Boundary Energy
average of detBαij

30

60

90

120

150

180

Misorientation angle, θ [degrees]

Misorientation angle, θ [degrees](c) Boundary energy and Minimumof detBαij at Al in Tilt model.

(d) Boundary energy and Number of detBαij<0 atoms
at Al in Tilt model.

Fig.4.6 Simulation results of tilt grain boundary（Al）．

49

3

2
1.5

1

1
0.5
0.5
0

0
0

30

60

90

120

150

1.5
1.5
1

1
0.5

0.5

0

0

-0.5

-0.5

180

0

Misorientation angle, θ [degrees](a) Boundary energy andaverage of detBαij at Ni in Tilt model.
Boundary Energy
Minimum of detBαij

1.5
2
1.5
1
0.5
0.5
0

0
0

30

60

90

120

150

90

120

150

(331) (332) (111)(334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ11

Boundary Energy [J/m 2]

Boundary Energy [J/m 2 ]

2.5

1

60

Boundary Energy
Number of detBαij<0 atoms

23(331) (332) (111)(334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ11

30

180

Misorientation angle, θ [degrees](b) Boundary energy and detBαij on GB
at Ni in Tilt model.

Minimum of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice2

Difference ∆AES
(ΣdetBijGB-ΣdetBijBulk)/2A[J/m 2 ]

2.5

1.5

2
(331) (332) (111)(334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ11

2

Boundary Energy [J/m 2]

Boundary Energy [J/m 2](331) (332) (111)(334) (112) (311)Σ19 Σ11 Σ3 Σ17 Σ3 Σ11

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

Boundary Energy
detBαij on GB

200

1.5

150
1
100
0.5

50

0

0
0

180

250

Number of detBαij<0 atoms

Boundary Energy
average of detBαij

30

60

90

120

150

180

Misorientation angle, θ [degrees]

Misorientation angle, θ [degrees](c) Boundary energy and Minimumof detBαij at Ni in Tilt model.

(d) Boundary energy and Number of detBαij<0 atoms
at Ni in Tilt model.

Fig.4.7 Simulation results of tilt grain boundary（Ni）．

50

4.2

Σ3 対称ねじれ粒界の引張シミュレーション

周期セルにおける粒界面の面積が前章の表面の大きさ程度になるようにセルの大きさを決定し，総原子数約 6 万の系を対象とする（表 4.1）．初期応力が 0 となるようにセル寸法を微調整しながら 100000fs の初期緩和シミュレーションを行った後，粒界と垂直方向（z 方向）に毎ステップひずみ増分 ∆εzz =1.0×10−6 を与えて引張シミュレー
ションを行った．この時，∆εzz =1.0×10−4 毎に原子座標データを記録し，各原子位置ααの原子応力 σij
ならびに原子弾性係数 Cijkl
を求め，AES の行列式 detBijα (Voigt 表記)を評価した．なお横方向 (x,y 軸方向) の応力は 0 となるようにセル寸法をスケーリン
グしている．

Table 4.1 Simulation conditions for Σ3 twist grain boundary．
CellΣ value Misorientation
length[nm]Atomicangle
AlNinumber◦Σ3
70.5
7.93×7.72×16.04 6.90×6.71×13.94 58432

Al の結果図 4.8 に Al のシミュレーション結果を示す．応力ピーク時のひずみ（εzz = 0.0902）
において detBijα の平均が 0 になっており，やはり系全体の不安定に達している．また，

detBijα の最小値はピークひずみより前のひずみ εzz = 0.0699 で急減し始めている．図4.8 左の最小値変化を見ると，引張直後は detBijα ＜0 の原子が発生するがその後 0 となる．これは次節で明確に観察された Ni の粒界キンクの発生と同じ現象がわずかながら起きていると考えられる．その後は detBijα の最小値が再び急減する εzz = 0.0699 まで detBijα ＜0 の原子は存在していない．図 4.8 右に示すように εzz = 0.0774 以降にゆるやかに増加しはじめ，その後 εzz = 0.0850 より前で急増する．
図 4.9 に応力ピーク近傍の主要な点における横から見た原子配置を detBijα 0 の値で色分けし，図 4.10 にその時の detBijα < 0 となった原子のみを緑色で着色して示した．負51の原子が増加し始める (c) では，図 4.10 より粒界面と近傍の原子が負となっていることがわかり，それは図 4.9(c) 赤丸で囲った部分と対応する．その後結晶内部に負の原子が拡大し，応力ピーク (e) においては粒界部分以外の原子が負となっている．
図 4.11 に Atom eye[40] の central symmetry parameter により「欠陥」と判別された原子を示す．応力急減前のひずみ (a)εzz = 0.0900 においては粒界面にしか存在していないが，応力-ひずみ曲線が折れ曲がる応力ピークの図 (b)εzz = 0.0902 において，粒界面から結晶内部に転位が発生し，(c)，(d) と欠陥が内部に成長していることがわかる．このように応力急減時には粒界の detBijα は負ではなく，結晶部分の原子が全体に負となった状態で粒界から転位が発生し結晶内部に成長する．

2

Maximum

1.5
2

εzz=0.0902

1
0.5

0

0

Minimum

Number of detBαij<0 atoms

0

0.05

0.1

0.15

40000
30000
20000
10000
0
0.05

0.1

0.15

0.2

1

4

Maximum2

εzz=0.0902
0.5

0
0.07

0.2

50000

0

Stress, σzz GPa

2

Stress-Strain curve

Stress-Strain curve
1.50Minimum
0.075

0.08

0.075

0.08

Stress, σzz GPa

4

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

2.5

Number of detBαij<0 atoms

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

3

0.085

0.09

0.095

0.1

0.085

0.09

0.095

0.1

50000
40000
30000
20000
10000
0
0.07

Strain, εzz

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

Fig.4.8 Change in average，standard deviation of detBijα ，and stress-strain
curve of Σ3 twist grain boundary（Al）．

52

value of det Bαij

1.0(a) εzz=0.0773(b) εzz=0.0774(c) εzz=0.0815

0.0

z

y

x

-1.0
(d) εzz=0.0884(e) εzz=0.0902

Fig.4.9 Snapshot of Σ3 twist grain boundary，colored by the detBijα（Al，side view）．

53

(a) εzz=0.0773(b) εzz=0.0774(c) εzz=0.0815zy(d) εzz=0.0884x(e) εzz=0.0902

Fig.4.10 Snapshot of Σ3 twist grain boundary．Green circles indicate
detBijα < 0 atoms（Al，side view）．

54

z

y

(a) εzz=0.0900(b) εzz=0.0902(c) εzz=0.0904(d) εzz=0.0906

x

Fig.4.11 Snapshot of Σ3 twist grain boundary．Circles indicate defect atoms（Al，side view）．

55

Ni の結果
図 4.12 に Ni の解析結果を図 4.12 に示す．Al と違い，薄膜と同じ程度のひずみ εzz =0.1119 で応力ピークを示して応力急減するが，ずっと前のひずみ εzz = 0.0206 において一度，応力-ひずみ曲線に小さいピークを生じ，それに対応して detBijα の平均が不連続になっている．detBijα ＜0 となった原子数や，detBijα の最小値にもこの初期変化に対応した挙動がみられる．右の拡大図をみると Ni のねじれ粒界も応力ピーク時にはdetBijα の平均が負になっており，系全体の不安定に達している．
図 4.13 に原子配置を detBijα の値で着色し，図 4.14 ではその時の detBijα ＜0 の物のみ着色して示した．応力ひずみの最初のピークの (c) において，粒界近傍の負または零の原子面が湾曲している．応力低下がとまるひずみ εzz = 0.0206 の図 4.13(d) から，
小さなピークは粒界面のキンク変形によるものであることがわかる．その後応力ピー
ク近傍では全ての原子が寒色系となっており全体に detBijα ＜0 に近づいている．応力急減時には図 4.13，4.14(f) に赤丸で囲ったキンク部分から内部に detBijα ＜0 の原子が発生する．

Central symmetry parameter により欠陥原子を表した図 4.15 を確認すると，(a)⇒(d)で先ほどの粒界部分のキンク変形が顕著に表れている．先ほど負の原子が発生した(e)εzz = 0.1008 ではキンク部からの欠陥原子は見受けられなく，応力ピーク直前の(f)εzz = 0.1015 においてキンク部からの転位発生が確認できる．

56

εzz=0.1119

0.55

εzz=0.0206
0

Average0

Minimum
0.05

0.1

0.15

30000
20000
10000
0
0.05

0.1

0.15

0.2

Stress-Strain curve
1

Strain, εzz(a) Outline10

εzz=0.1119

0.5

5

Maximum
Average00Minimum

0.09

0.2

40000

0

Stress, σzz GPa

10

15

Stress, σzz GPa

Maximum

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice10

Number of detBαij<0 atoms

1.5

15

Stress-Strain curve

Number of detBαij<0 atoms

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

1.5

0.095

0.1

0.095

0.1

0.105

0.11

0.115

0.12

0.105

0.11

0.115

0.12

40000
30000
20000
10000
0
0.09

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

Fig.4.12 Change in average，standard deviation of detBijα ，and stress-strain
curve of Σ3 twist grain boundary（Ni）．

57

value of det Bαij

1.0z(b) εzz=0.0175(a) εzz=0.0171
0.0 y(c) εzz=0.0182x-1.0(d) εzz=0.0206(e) εzz=0.0774(f) εzz=0.1008

Fig.4.13 Snapshot of Σ3 twist grain boundary，colored by the detBijα（Ni，side view）．

58

z

y

(a) εzz=0.0171(b) εzz=0.0175(c) εzz=0.0182(d) εzz=0.0206(e) εzz=0.0774(f) εzz=0.1008

xFig.4.14 Snapshot of Σ3 twist grain boundary．Green circles indicate
detBijα < 0 atoms（Ni，side view）．

59

z

y

(a) εzz=0.0171(b) εzz=0.0175(c) εzz=0.0182(d) εzz=0.0206(e) εzz=0.1008(f) εzz=0.1105

x

Fig.4.15 Snapshot of Σ3 twist grain boundary．Circles indicate defect atoms（Ni，side view）．

60

4.3

Σ3 対称傾角粒界における引張シミュレーション

対称ねじれ粒界と同様に，粒界面の面積を前章の表面の大きさ程度とした，総原子数約 6 万の系を対象とする．表 4.2 に寸法を示す．初期応力が 0 となるようにセル寸法を微調整しながら 100000fs の初期緩和シミュレーションを行った後，粒界と垂直方向（y 方向）に毎ステップひずみ増分 ∆εyy =1.0×10−6 を与えて引張シミュレーションを行った．∆εyy =1.0×10−4 毎の原子座標データを記録して，AES の行列式 detBijα (Voigt
表記) を評価した．また横方向 (x,z 軸方向) の応力が 0 となるようにセル寸法をスケーリングしている．

Table 4.2 Simulation conditions for Σ3 tilt grain boundary．
GBboundary
plain
{111}

Cell
length[nm]Σ value Misorientation
angle
Al
◦
Σ3
70.5
7.93×16.86×8.02

61

AtomicNi
number
6.90×14.64×6.97 64512

Al の結果
Al の傾角粒界の解析結果を図 4.16〜4.18 に示す (detBijα で色分けした原子配置図は省略)．薄膜モデルとそれほど変わらないひずみ εyy = 0.1382 で応力ピークを示す．ピー
ク時に detBijα の平均が 0 になっており，やはり系全体の不安定に達している．また，

detBijα ＜0 の原子はピーク直前まで 0 であり，かつ応力ピーク直前に detBijα の最小値
が折れ曲がって負値を示すことから局所不安定をトリガーとした系全体の不安定ととらえることができる．
図 4.17 の detBijα ＜0 となった原子を見ると，(a) の εyy = 0.1372 から (c) の εyy =0.1374 のごくわずかな時間に detBijα ＜0 の原子が粒内から急増し，応力ピークの εyy =0.1382 ではやはり粒界部以外のすべての部分で負となっている．
図 4.18 の欠陥原子を確認すると，応力ピークを示す (a)εyy = 0.1382 においては結晶
内部に格子欠陥は存在しておらず，その後応力ピーク後の (b)〜(d) で粒内に集団的な欠陥原子が現れており，ねじれ粒界のときとは違い粒界部分からの転位ではなく，粒内で相変態的な変形を生じている．

62

10
1

Maximum50.5

εyy=0.1382
Minimum

0

0.05

50000
40000
30000
20000
10000
0

0
0.1

0.15

εyy=0.1372
0

0.05

0.1

Strain, εyy

0.15

0.2

Stress-Strain curve
100.25(a) Outline

εyy=0.1382
5

Maximum
Average
0
0.12

0.2

15

Stress, σyy GPa

Stress-Strain curve

Stress, σyy GPa
Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice

1.50

Number of detBαij<0 atoms

0.5

15

Number of detBαij<0 atoms

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect lattice20Minimum
0.125

0.13

0.135

0.14

0.145

0.15

50000
40000
30000
20000
10000
0
0.12

εyy=0.1372
0.13

0.14

0.15

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

Fig.4.16 Change in average，standard deviation of detBijα ，and stress-strain
curve of (111) Σ3 tilt grain boundary（Al）．

63

y

z

(a) εyy=0.1372(b) εyy=0.1374(c) εyy=0.1375(d) εyy=0.1382

x

Fig.4.17 Snapshot of (111) Σ3 tilt grain boundary．Green circles indicate
detBijα < 0 atoms（Al，side view）．

64

y

z

(a) εyy=0.1382(b) εyy=0.1390(c) εyy=0.1393(d) εyy=0.1396

x

Fig.4.18 Snapshot of Σ3 tilt grain boundary．Circles indicate defect atoms（Al，side view）．

65

Ni の結果
図 4.19〜4.21 に Ni における解析結果を示す．これまでの系の倍近くのひずみ εyy =0.1872 で応力ピークを示している. detBijα が負の原子数は初期状態では 0 で，応力ピー
クを迎える直前から負の原子が現れ急増する．他の粒界モデルは応力-ひずみのピーク
点付近で detBijα の平均が負となっていることから「系の不安定条件」に達したため系の不安定挙動が発生したものと考えられるが，Ni の (111)Σ3 対称傾角粒界においては0 に近づいてはいるが負値にはならなかった．また，detBijα の最小値もきわめて高い値を示しているのも特徴である．図 4.20 を見ると，こちらも Al の Σ3 傾角粒界と同様に，粒界部からではなく結晶内部から負の原子が発生している．図 4.21 の欠陥原子を確認すると，応力ピークを示す (a)εyy = 0.1382 においては，格子欠陥は存在していな
い．その後 detBijα ＜0 の原子が集団的に発生している部分に欠陥と判定された原子が
現れる．

0.5

2

20

10

Minimum
0

0
0.05

0.1

0.15

α

10000
5000
0
0

0.05

0.1

Strain, εyy

0.15

0.2

(a) Outline

Maximum

20

10

Average00Minimum
0.17

0.2

15000

εyy=0.1872

0.25

Number of detB ij<0 atoms

0.5

0Number of detBαij<0 atoms

εyy=0.1872

Average

1

Stress, σyy GPa

Maximum

30

Stress, σyy GPa

1.5

30

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

Average of 6x6 determinant, ΣdetBαij/N,normalizzed by detBαij at perfect latticeStress-Strain curve
Stress-Strain curve

0.175

0.18

0.175

0.18

0.185

0.19

0.195

0.2

0.19

0.195

0.2

15000
10000
5000
0
0.17

εyy=0.1870
0.185

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

Fig.4.19 Change in average，standard deviation of detBijα ，and stress-strain
curve of (111) Σ3 tilt grain boundary（Ni）．

66

y

z

(a) εyy=0.1870(b) εyy=0.1871(c) εyy=0.1872(d) εyy=0.1873

xFig.4.20 Snapshot of (111) Σ3 tilt grain boundary．Green circles indicate
detBijα < 0 atoms（Ni，side view）．

67

y

z

(a) εyy=0.1872(b) εyy=0.1873(c) εyy=0.1874(d) εyy=0.1875

xFig.4.21 Snapshot of Σ3 tilt grain boundary．Circles indicate defect atoms（Ni，side view）．

68

4.4

結言

様々な角度の対称ねじれ粒界ならびに対称傾角粒界を作成し，無負荷平衡状態における回転角度 θ に対する粒界エネルギーならびに，detBijα の平均や最小値，detBijα ＜0 の数を調べた．その後，Σ3 ねじれ粒界ならびに Σ3 傾角粒界について，粒界に垂直方向の引張を分子動力学法により行い，粒界固有の強度について detBijα の視点から検討した．得られた結果を以下に示す．

(1) 低いエネルギーカスプを示す，Σ3 ねじれ粒界，(111)Σ3，(311)Σ11 傾角粒界では detBijα の最小値が正で detBijα ＜0 の原子は存在しない．

(2) detBijα の最小値や detBijα ＜0 の数は必ずしも粒界エネルギーの大小関係とは対応しないが，Σ3 粒界などの最安定な構造については detBijα の最小値がバルクの半分程度の高い値を示すなどの一定の対応関係も見られた．

(3) Σ3 ねじれ粒界の引張シミュレーションにおいて，Ni は初期に粒界キンクが発生した．Ni，Al とも応力ピーク時には detBijα の平均が 0 となり，「系の不安定条件」に達していた．

(4) (3) の応力ピーク時には粒界以外の原子は全て detBijα ＜0 となっていた．Atomeye の central symmetry parameter により欠陥原子の発生を観察したところ，粒界面を起点として結晶内部に転位が発生していた．

(5) (111)Σ3 傾角粒界の引張では，Al，Ni ともピークひずみが大きくなり，特に Niは他の薄膜や細線などのシミュレーションの倍程度と大きな値を示した．

(6) Al，Ni ともに粒界ではなく粒内にクラスター状に detBijα ＜0 となった原子群が発生し，Al の場合は (4) と同様，粒界以外が負となり系の不安定に達していた．

(7) Central symmetry parameter による欠陥生成を観察した結果，Ni，Al とも (6) のクラスター状の負の原子群部分が欠陥と判定されていた．

69

第5章
結論
本研究では fcc 金属の Al ならびに Ni について，表面や粒界という初期不均一構造が引張時の系の不安定挙動=応力急減におよぼす影響を明らかにすることを目的として，各原子位置における原子弾性剛性係数 detBijα の正値性に基づいて検討した．以下に得られた結果を総括する．

2 章では本研究で用いた解析手法の基礎について述べた．まず，分子動力学法の概要ならびに基礎方程式を示し，本研究で用いた数値積分法について説明した．その後
粒子間相互作用の評価に用いた EAM ポテンシャルについて述べた．また，AES による局所力学状態評価についても説明した．

3 章では引張時の系の不安定挙動=応力急減におよぼす表面不均一性の影響を調べるため，(100)，(110)，(111) 表面を有する薄膜や細線の引張シミュレーションを行い，応力-ひずみ曲線を求めるとともに detBijα の正値性の観点から検討した．Al の (111) 表面
薄膜，Ni の全表面薄膜 ((100)，(110)，(111)) と (110) 表面細線は，応力-ひずみのピー
ク点付近で detBijα の平均が負となっていることから「系の不安定条件」に達して応力急減=引張不安定を生じたと結論づけられた．それ以外の系でも Ni は全体的に detBijαの平均は 0 に近い値をとっているが，Al の他の系は比較的大きな正の値のときに応力ピークを生じており，系全体の不安定に達する前に応力急減している．多くの系では，応力ピークよりかなり前に detBijα の最小値が減少して負値をとり，表面やエッジ部に負の原子が発生・拡大していた．しかしながら Al(100) 薄膜，(110) 薄膜は応力ピーク
後に負の原子が発生していること，かつ応力ピーク時にも detBijα の平均が大きな正値70をとっていることから，系の不安定挙動の予測は困難だった．

4 章では様々な角度の対称ねじれ粒界ならびに対称傾角粒界を作成し，無負荷平衡状態におけるミスオリエンテーション角度 θ に対する粒界エネルギーならびに，detBijαの平均や最小値，detBijα ＜0 の数を調べた．Σ3 ねじれ粒界，(111)Σ3，(311)Σ11 傾角
粒界においては低いエネルギーカスプを示し，detBijα の最小値が正値で detBijα ＜0 の原子は存在しないことを明らかにした．detBijα の最小値や detBijα ＜0 の数は必ずしも
粒界エネルギーと対応しないが，Σ3 粒界などの最安定な構造については detBijα の最小値がバルクの半分程度の高い値を示すなどの一定の対応関係も見られた．Σ3 ねじれ
粒界の引張では Ni は引張初期に粒界キンクの形成が見られた．Ni，Al とも応力ピー
ク=detBijα の平均が 0 となり，「系の不安定条件」に達している．応力ピーク時は粒界以外が全て detBijα ＜0 となっており，Atom eye の central symmetry parameter により欠陥原子を確認すると，粒界から結晶内部に転位が発生していた．一方 (111)Σ3 傾
角粒界の引張では，Al，Ni ともピークのひずみが大きくなり，特に Ni は他の薄膜や細線などのシミュレーションの倍程度と大きな値を示した．応力ピーク時の detBijα の平均は Ni はわずかに正値であったが Al は負となっており系の不安定に達していた．さらに，Ni，Al ともに粒界ではなく粒内にクラスター状に detBijα ＜0 となった原子群が発生し，Al の場合はねじれ粒界と同様に，粒界以外が負となり系の不安定に達した．
欠陥原子を見ると，Ni，Al ともに応力ピーク後に粒内に欠陥原子が集団的に発生しており，相変態的な挙動を示した．

71

参考文献
[1] L.Zhou and Y.Yao，Mater.Sci.Eng.A，460-461，(2007)，95-100．
[2] 吉田・李木・紀・奥宮・恒川，J. Japan Inst. Metals，73-2，(2009)，95-102．
[3] Gerberivh，W.W. et al.，Acta Materialia，44-9，(1996)，3585-3598．
[4] Suresh，S．et．
，Scripta Materialia，41-9，(1999)，951-957．
[5] H.A.Wu，Mechcanics Reserch Communications，33-1，(2006)，9．
[6] 北村，材料，44-2，(1995)，25．
[7] Y.Gan，J.K.chen，Mechanics Reserch Communications，36-7，(2009)，838．
[8] Y.Tong，F.Yi，Lisheng.Liu，Pengcheng.Zhai，Qinging.Zhang，Phisica B: Condensed Matter，405-15，(2010)，3190．
[9] H.A.Wu，Mechcanics Reserch Communications，33-1，(2006)，9．
[10] Y.H.Wen，Zi-Zhong.Zhu，Ru-Zeng.Zhu，Computational Materials Science，414，(2008)，553．
[11] A.J.Devey，R.Grau-Crespo and N.H. de Leeuw，J.Phys. Chem. C 112(29)，(2008)，10960．
[12] M.Grujicic，P.Dang，Matrials Science and Engineering，201-1，(1995)，194．
[13] G.A.Chadwick and D.A.Smith，Chap，10，(1976)，Acamedic Press．
[14] T.Kwork, P.S.Ho and S.Yip，Sci，144，(1984)，44．
[15] 北村, 梅野，大谷，日本機械学会論文集 (A 編)，62, (1996)，2791．
[16] 北村・大谷・山中・屋代，日本機械学会論文集 (A 編)，60, (1994)，2843．
72

[17] 岩崎・佐々木・千葉・阿部，日本機械学会論文集 (A 編)，61, (1995)，1376．
[18] 北川・中谷・尾方・齋藤・前川，日本機械学会論文集 (A 編)，62, (1996)，853．
[19] K. Morita and H. Nakashima，Mater. Sci. Eng.，A234-236, (1997), 1053-1056.
[20] 森田・中島，九州大学大学院総合理工学研究科報告，19, (1997), 1-8.
[21] 森田・上原・連川・英治，日本金属学会誌，61，(1997)，251-260.
[22] 池田・森田・中島・阿部，日本金属学会誌，63，(1999)，179-186.
[23] 中島・竹内，鉄と鋼，86，(2000)，357-362.
[24] 周致霆・黒木康洋・池田賢一・吉田冬樹・中島英治・阿部弘，九州大学大学院総合理工学報告，26，(2004)，333-336.

[25] 周致霆・池田賢一・中島英治，日本金属学会誌，69，(2005)，303-307.
[26] N. Takata，T. Mizuguchi，K. Ikeda，and H. Nakashima，Meter. Trans.，45，(2004)，2099-2105.
[27] 水口・高田・池田・中島英治，日本金属学会誌，69，(2005)，1010-1015.
[28] 屋代・冨田，日本機械学会論文集 (A 編)，67,-656 (2001)，678．
[29] 屋代・西村・冨田，材料，54(10)，(2005)，1053．
[30] 屋代・藤原，材料，60(11)，(2011)，968．
[31] Yashiro, K. and Fujihara, M.,Modell. Simul. Mater. Sci. Eng., 20-4, (2012)．
[32] 屋代他，第 17 回分子動力学シンポジウム，2012[33] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
[34] 上田顯, 分子シミュレーション-古典系から量子系手法まで-, (2003), 裳華房.
[35] Y.Misin，Acta Materialia，52，(2004)，1451-146773

[36] National
ment

InstituteofStandards

Laboratory，Interatomic

andPotentials

Technology，Material
Repository

Measure-Project，(2004)，(http://www.ctcms.nist.gov/potentials/Al-Ni.html)[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] Ju Li，Modelling Simul.Mater.Sci.Eng, 11-173, (2003).

74学術論文・学術講演
学術講演

 Ni 基超合金γ / γ’ 構造の安定性ならびに形態変化に関する
MC シミュレーション
池宮

一繁，山口

明宏，屋代

如月

日本材料学会 第 1 回マルチスケールマテリアルモデリングシンポジウム (第 16回分子動力学シンポジウム), 大阪大学, (2011.5)

 原子弾性剛性係数による Fcc 金属の表面・界面力学特性解析
池宮

一繁，屋代

如月

日本機械学会 第 25 回計算力学講演会 CMD2012, 神戸ポートアイランド南地区,(2012.10)75

論文貼り付け用ページ76論文貼り付け用ページ

77

論文貼り付け用ページ

78論文貼り付け用ページ

79論文貼り付け用ページ80

謝辞
本研究を遂行するに当たり，屋代如月准教授の研究の指針および要所で深い洞察を伴った適切なご助言を下さり，日々の研究を円滑に行うことのできる環境を整えて下さったことは著者が研究を進める上で大きな助けとなりました．浅識な著者が本論文を完遂できましたのも屋代准教授の厚い御指導の賜物であり，心より御礼申し上げます．また，研究報告会などで広い視点から御助言を下さり，著者の見解を広げて下さった田中克志教授と木下博准教授 (現岡山大学) に感謝致します．

本研究室配属後，温かく接していただいた表面・界面工学研究室の諸先輩方には研究室生活全般におよびご助力，ご指導いただきましたことに感謝致します．特に研究内容ついて詳しく教えていただいた山口明宏様（現三菱重工業株式会社）に深く感謝いたします．研究生活を共に過ごしてきた同輩の住谷昂大氏，西村英晃氏，中田伸哉氏，藤原正大氏，にはお互いに励ましあい，支えてくれたことに感謝しております．
特に，同じ研究チームの藤原正大氏ならびに，後輩ながらも互いに切磋琢磨し合った，
山根尭氏には更なる感謝の意を申し上げます．また，今日に至るまで数々のご支援とご協力を頂いた, 表面・界面工学研究室諸氏にも，ここに記することで感謝の意を表します．
最後に, 6 年間の大学生活を支えて頂いた家族や友人に心から最大限の感謝の意を表します．
平成 25 年 2 月
池宮

81

一繁