修 士 論 文(題 目)fcc, bcc, hcp 単結晶のせん断変形
シミュレーションと双晶変形発生機構の
原子弾性剛性係数による検討
令和 3 年度

岐阜大学大学院
自然科学技術研究科
物質・ものづくり工学専攻

氏名 土屋周平目次
第 1 章 緒言1

第 2 章 解析手法の基礎

4

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

4

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

5

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

6

2.4 高速化手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

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

8

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

9

2.7 原子応力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.8 原子弾性係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

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

12

2.10 原子弾性剛性係数の固有値と固有ベクトル . . . . . . . . . . . . . . .

12

第 3 章 bcc 単結晶の [111](1̄1̄2) せん断解析

14

3.1 シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . .

14i

3.2 静力学解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.3 応力-ひずみ曲線および固有値変化 . . . . . . . . . . . . . . . . . . . .

20

3.4 bcc の局所変形発生の詳細 . . . . . . . . . . . . . . . . . . . . . . . .

22

3.5 固有ベクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

第 4 章 fcc 単結晶の [1̄1̄2](111) せん断解析

31

4.1 シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.2 静力学解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.3 応力-ひずみ曲線および固有値変化 . . . . . . . . . . . . . . . . . . . .

37

4.4 Ag の局所変形発生の詳細 . . . . . . . . . . . . . . . . . . . . . . . .

41

4.4.1

欠陥発生と η a(1) <0 の分布 . . . . . . . . . . . . . . . . . . . . 41

4.4.2

η a(2) <0 の分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.4.3

固有ベクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.5 Cu の局所変形発生の詳細 . . . . . . . . . . . . . . . . . . . . . . . .

47

4.5.1

欠陥発生と η a(1) <0 の分布 . . . . . . . . . . . . . . . . . . . . 47

4.5.2

η a(1) と η a(2) の分布 . . . . . . . . . . . . . . . . . . . . . . . . 49

4.5.3

原子移動による変形モード評価 . . . . . . . . . . . . . . . . . 51

4.6 Pt の局所変形発生の詳細 . . . . . . . . . . . . . . . . . . . . . . . . .

53

第 5 章 hcp 単結晶の [1̄012](1̄011) せん断解析

57

5.1 シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . .

57

ii

5.2 静力学解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

5.3 応力-ひずみ曲線および固有値変化 . . . . . . . . . . . . . . . . . . . .

63

5.4 Mg の局所変形発生の詳細 . . . . . . . . . . . . . . . . . . . . . . . .

66

第 6 章 結言

72

参考文献

75

iii

第 1 章 緒言
分子動力学 (MD) シミュレーションは微視的メカニズムを観察することが出来る
有力なツールではあるが，ポテンシャルの精度の問題から得られる臨界ひずみや応力などはモデル材料のものであり，定量的な議論に結びつけることは難しい．我々のグループでは，原子弾性剛性係数 (Atomic Elastic Stiffiness, AES)Bija = ∆σia /∆εjの固有値と固有ベクトルに基づいて局所変形開始を統一的に議論する検討を行ってきた (1) ．ここで ∆σia は原子応力，∆εj はひずみで指標 i, j は i, j = 1 ∼ 6 =xx, yy, zz, yz, zx, xy を表す．固有方程式 Bija ∆εj =η a ∆εi の第一固有値 η a(1) が負となった原子を変形方向と逆に応力が変化する不安定状態とみなし，対応する固有ベ
クトル ∆εi ={∆ε1 , · · · , ∆ε6 }={∆εxx , · · · , ∆γxy } から変形モードを議論する．変形時の臨界応力やエネルギーなどのしきい値でなく，「正か負か」という指標で議論することで普遍的なメカニズムの抽出を目指している．
前報では双晶変形の発生を AES で議論することを目的として，bcc Fe 単結晶のせん断シミュレーションを行った (1) ．その結果，応力ピーク直前に各原子面毎の η a(1)にゆらぎが発生し，一定間隔で η a(1) <0 となった周期的な不安定構造を生じること，
その後実際に η a(1) <0 の位置に双晶変形を生じること，周期構造の間隔は周期方向にセルサイズを変えても不変である事などを明らかにしている．
双晶は同種の単結晶がある一定の角度で対称となるように接合されたものであり，双晶変形は結晶が力を受けて双晶を形成するような変形である (Fig. 1. 1)．双晶変形
はすべりが起こり難い結晶または低温・高速変形の条件下で起こりやすく，双晶変形1を生じる際の原子移動の基本単位が格子ベクトルと一致しないことなどから，すべりとマルテンサイト変態による変形様式の中間に位置する変形とみなせる (2) ．すべ
りにすべり面やすべり方向があるように双晶にも双晶面や双晶方向があり，双晶の有無で強度が変化する事でも知られている．双晶変形に関する MD シミュレーショ
ン研究として，Sainath と Choundhary は bcc Fe のナノワイヤーの変形挙動に対する向きによる影響を MD により調べ，引張・圧縮の変形方向の違いで双晶が発生する結晶方位が異なることを示している (3) ．Wu らは bcc Fe の異なる温度およびひずみ速度での引張 MD シミュレーションを行い，温度が 400K 未満になると双晶の成長に必要な応力が劇的に低下することを報告している (4) ．Zhang らは高圧力による双晶発生の促進を議論するため，衝撃荷重下での bcc 結晶の双晶変形のシミュレー
ションを行っている (5) ．Yamakov らは 30∼100nm の粒径の柱状ナノ結晶 Al の変形プロセスを MD により解析し，双晶がひずみ硬化をもたらす可能性があるという実験結果と同様の傾向を報告している (6) ．Wu と Zhu は Ni ナノ結晶中の双晶変形のMD シミュレーションを行い，実験で報告されている逆粒径依存性とシミュレーショ
ンの差を報告している (7) ．Qi らは約 2nm の Cu ナノ銅ワイヤーの引張 MD を行い，変形が双晶とすべりによることを報告している (8) ．Shao らは hcp Mg のせん断 MD
シミュレーションを行い，複数層の双晶の発生・成長には高い応力が必要である可能性を示している (9) ．Jian らは高圧衝撃荷重下での Mg の可塑性と相転移に関して

MD シミュレーションを行っている (10) ．Minrong らは hcp Ti の引張・圧縮の MD
シミュレーションを行い，双晶変形が発生結晶方位と荷重方向に強く依存していることを示している (11) ．
本研究では前報と同様のせん断シミュレーションを，同じ GEAM ポテンシャル(12)−(14)で bcc 4 元素 (既報の Fe 含む)，fcc 8 元素，hcp 4 元素について行い，原子弾性剛性係数 Bija = ∆σia /∆εj の固有値・固有ベクトルの観点から議論する．ただし，前報のようにせん断面に平行な双晶変形が発生しないケースが多く，発生した局所

2

変形 (せん断方向とは異なる) の発生メカニズムについて議論する．

Fig. 1.1

Schematic of twin deformation.

3

第 2 章 解析手法の基礎
分子動力学法2.1

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

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

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

α

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

dr α= v α (t)dt

(2.1.5)4

であり，式 (2.1.1) と式 (2.1.5) を式 (2.1.3) と式 (2.1.4) に代入すると，)((∆t)2 F α (t)
3
+
O
(∆t)
2
mα
)
(
(∆t)2 F α (t)
3
+
O
(∆t)
r α (t − ∆t) = r α (t) − ∆tv α (t) +2mαr α (t + ∆t) = r α (t) + ∆tv α (t) +(2.1.6)(2.1.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.1.8)(2.1.9)が得られる．∆t3 以上の高次項は無視できるとすると，時刻 t + ∆t での位置ベクトルと t での速度は

r α (t + ∆t) = 2r α (t) − r α (t − ∆t) + (∆t)2v α (t) =1{r α (t + ∆t) − r α (t − ∆t)}2∆t

F α (t)mα(2.1.10)(2.1.11)と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．初期の計算 (t = 0) では，t = ∆t での座標 r α (∆t) は式 (2.1.6) と初速度から得ることができる．

2.2

原子間ポテンシャル

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

5

EAM ポテンシャル

2.3

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

Etot =∑F α (ρα ) +α1 ∑ ∑ αβ αβφ (r )
2 α β(̸=α)(2.3.1)

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

で評価する．
本研究では，多数の元素に対してパラメータが提案されている GEAM(Generalized

Embedded Atom Method) ポテンシャルを用いた (12)–(14) ．GEAM ポテンシャルでは，
Etot =∑αEα =∑{ 1 ∑α2 β(̸=α)φαβ (rαβ ) + F α (ρα )}(2.3.3)

で系のエネルギーを評価する．E α は原子 α のエネルギーを表す．

2.4

高速化手法

原子数 N の系において粒子間の全相互作用を評価すると，1 ステップに N ×(N − 1)回の計算が必要となり，N が大きくなると極めて膨大な計算量となる．実際には，

6

一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ
半径 rc ) 内の粒子からの寄与を効率よく計算することにより高速化できる．本研究で用いたブロック分割法は，シミュレートする系をカットオフ距離程度の格子状に分割し，各ブロックに属する粒子をメモリーに記憶する．着目している粒子に作用する力を評価する際には，Fig. 2.4.1 に示すように，その粒子が属するブロックおよび隣接するブロックから相互作用する粒子を探索して行う．粒子が属するブロック
は，粒子の位置座標をブロックの辺長 bx，by で除した際の整数により判断できるので，ブロック登録時の計算負荷は粒子数 N のオーダーとなる．

Fig. 2.4.1

Schematic of domain decomposition method.

7

2.5

速度スケーリング法

分子動力学解析における温度制御には一般的には速度スケーリング法が用いられる．統計熱力学より系のエネルギーと温度には以下の関係が成立する．

1 α α α 3
m v i v i = kB T
22

(2.5.1)

ここで mα は粒子 α の質量，viα は 温度 T での粒子 α の速度ベクトルの i 方向成分，

kB は Boltzmann 定数で kB =1.38 × 10−23 [J/K] である．目標の温度 T0 における原子α の速度成分を viα0 とすると
viα0 =√

3kB T0mα(2.5.2)であり，温度 T の時の原子 α の速度成分は

viα =

√

3kB T
mα

(2.5.3)となる．したがって

viα0 =

√

T0 α
v
T i

(2.5.4)となる．現在の速度から求めた T を用いて，各原子の速度を√

T0 /T スケールすれば目的の温度 T0 となる．ただ，これでは次の原子配置に反映されないので，Verlet
法における式 (2.1.10) を以下のように変更する．

r α (t + ∆t) = r α (t) +√T0 α
F α (t){r (t) − r α (t − ∆t) + i α ∆t2 }Tm(2.5.5)平衡状態では，能勢の方法 (18) など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている．

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

2.6

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

1
=
V

σij

(

∂F
∂ηij

)

1
, Cijkl =
V

(

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

Bijkl ≡

(

∆σij
∆εkl

)

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

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

2.7

原子応力

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

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

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

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

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

ηij =(2.7.4)であり，その微小量dηij =

1dGij
2(2.7.5)と式 (2.7.3) の関係から次の関係が得られる．

∂rαβ
∂ηij

r̄iαβ r̄jαβ=
rαβ(2.7.6)

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

σij

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

) αβ αβ
N ∑
N (
r r
1 ∑
1 ′ αβ αβ′α α′ β αβ=φ (r ) + F (ρ )f (r ) i αβjV α β(̸=α) 2r

10ここで，原子体積を Ω = V /N とすれば，各原子のエネルギー寄与

Etot =∑F α (ρα ) +α1 ∑ ∑ αβ αβφ (r )
2 α β(̸=α)の微分から次のように原子応力を定義できる．

σijα =

2.8) αβ αβN (r r
1 ∑1 ′ αβ αβφ (r ) + F ′ α (ρα )f ′ β (rαβ ) i αβjΩ β(̸=α) 2r

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

1 ∂ 2 Etot
V ∂ηij ∂ηkl

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

Cijkl

(N ∑
} αβ αβ )
N
N {
N ∑
r r
1 ∑
∂rαβ ∂ ∑1 ′ αβ αβ′α α′ β αβ=φ (r ) + F (ρ )f (r ) i αβjαβ
2V α β(̸=α) ∂ηkl ∂rr
α β(̸=α) 2} αβ αβ αβ αβ
N {N ∑φ′ αβ (rαβ ) ri rj rk rl
1 1∑′′ αβ αβφ (r ) −=V 2 α β(̸=α)
rαβ(rαβ )2
[

+

N
N ∑
∑

{

F ′ α (ρα ) f ′′ β (rαβ ) −α β(̸=α)

+

N
∑

F

′′ α

α

} αβ αβ αβ αβ
f ′ β (rαβ ) ri rj rk rl
rαβ(rαβ )2rαβ rαβ }{ ∑ ′ β αβ rkαβ rlαβ }f (r ) αβ(ρ )f (r ) i αβjrrβ(̸=α)
β(̸=α)
α

{ ∑

′β

]

αβ

(2.8.2)応力と同様に，各原子位置における原子弾性係数を以下のように定義する．
αCijkl} αβ αβ αβ αβ
N {
1 1 ∑
φ′ αβ (rαβ ) ri rj rk rl=
φ′′ αβ (rαβ ) −Ω 2 β(̸=α)
rαβ(rαβ )2
[

+

N
∑

{

F ′ α (ρα ) f ′′ β (rαβ ) −β(̸=α)+F′′ α} αβ αβ αβ αβ
f ′ β (rαβ ) ri rj rk rl
rαβ(rαβ )2

riαβ rjαβ }{ ∑ ′ β αβ rkαβ rlαβ }f (r ) αβf (r ) αβ(ρ )rrβ(̸=α)
β(̸=α)
α

{ ∑

′β

αβ

11

]

(2.8.3)

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

原子弾性剛性係数

2.9以上で定義した原子応力，弾性係数から，Voigt 対称性を持たせた原子弾性剛性係数は以下で評価できる．
α
α
α
α
Bijkl
= Cijkl+ (σilα δjk + σjlα δik + σik
δjl + σjk
δil − 2σijα δkl )/2(2.9.1)α
ここで i〜l は直交座標の指標 (1,2,3) もしくは (x，y，z) であり，Bijkl
は 4 階のテンソルであるが，Voigt 対称性を持たせた場合，弾性係数テンソルと同じく独立な成分
は 21 個となり，xx，yy，zz，yz，zx，xy を 1-6 とする Voigt 表記を用いれば 6 × 6 のマトリックス Bijα として次のように表すことができる．

[




] 

Bijα = Cijα +



] [

2.10

2σ1α

−(σ1α + σ2α )2σ2α−(σ1α + σ3α )−(σ2α + σ3α )
2σ3α

sym.

−σ4α
σ4α
σ4α
α
σ2 + σ3α

σ5α
−σ5α
σ5α
σ6α
α
σ3 + σ1α


σ6α

σ6α

α

−σ6
/2
α

σ5

α

σ4
α
α
σ 1 + σ2(2.9.2)原子弾性剛性係数の固有値と固有ベクトル

Bijα の固有方程式
Bijα ∆εj = η α ∆εi(2.10.1)を解くことは，6×6 マトリクス [Bijα ] を次のように対角化することと等しい．

12

α
B11

α
B12
α
B22

α
B13
α
B23
α
B33

α
B14
α
B24
α
B34
α
B44

α
B15
α
B25
α
B35
α
B45
α
B55

sym.

α
η α(1)
B16
α
0
B26
α
B36
0
→
α
B46
0
α
B56
0
α
B66
0

0
η α(2)000000η α(3)000000η α(4)000000η α(5)000000η α(6)

6 つの固有値 η a(1) , η a(2) ,. . . , η a(6) に対応する固有ベクトル ∆εi ={∆ε1 , · · · , ∆ε6 }={∆εxx , · · · , ∆γxy }のせん断ひずみ成分を 1/2 してテンソルひずみに変換し，各固有ベクトルに対して
∆εxx ∆εxy ∆εxz
∆εyy ∆εyz 


sym.
∆εzz

(2.10.2)の 3 × 3 テンソルを求める．さらにこれの主軸，すなわちこのテンソルの固有ベク
トル x1 < x2 < x3 を求め，x1 + x3 を最大せん断方向と定義して評価する．

13

第3章

3.1

bcc 単結晶の [111](1̄1̄2) せん
断解析

シミュレーション条件

双晶変形が発生する双晶面は bcc は (112) 面とされているので，双晶面を xz 面とする．Fig. 3.1.1 は x=[111]，y=[1̄1̄2]，z=[11̄0] としたときの bcc の直方体単位格子である．これを 23×46×16 並べたセルで解析を行った．総原子数は 101 568 個であり, 全方
向周期境界を適用している．MD 計算では熱による影響を排除するために温度を 0.1Kとした．10 000step の緩和計算を行った後に xy 方向にひずみ速度 ∆γxy = 1.0×10−6 /fs
で周期セルを斜方体に変形させてせん断変形させた (Fig. 3.1.2)．静力学解析では MD計算を行わずに全ての原子を格子点に固定してせん断させた．静力学解析では全ての原子は同じ状態である．なお Fe の MD 結果は前報 (1) のセル 25×100×4 の結果を流用し，今回は Mo, Ta, W について新たに行った．

14

B

B

A

A

B

B
A

y

y

A
B

B
A

A

z

x
Fig. 3.1.1

bcc unit lattice.

=

y'
y
x

Fig. 3.1.2

xSimple shear under oblique periodic boundary condition.

15

3.2

静力学解析

Fig. 3.2.1 に bcc 原子 4 種類 (Fe, Mo, Ta, W) の静力学計算によって得られた応力ひずみ曲線と第一固有値 η a(1) の変化を示す．静力学なので応力ピーク後も局所変形による応力急減を生じることはなく，正弦波状の滑らかな曲線を描き，全ての原子でひずみ 0.339∼0.340 で応力が 0 となる．Fig. 3.2.2 に Mo の γxy =0.34 におけるス
ナップショットを示す．この構造は x[100]，y[010]，z[001] の bcc 構造となっているが，固有値が負であることから安定な格子長さではない．その後はこの点について
点対称に圧縮側に応力ピークを示す．ピークひずみは W, Mo と Fe, Ta がそれぞれほぼ同じで，各ピーク点は η a(1) が負になる点と一致する．すなわち系の不安定点であることが示される．Fe と Ta は応力が 0 となるひずみを頂点として 2 次関数の様に変化し，Mo と Ta は同ひずみで左右対称となる 4 次関数の様に変化した．Fig. 3.2.3
は各元素の他の固有値の変化も示した図である．固有値 η a(1) ∼ η a(6) の 6 つのうち小さい順に第四固有値まで示している．固有値が負となったのは 4 原子種とも η a(1)のみである．Fe, Mo, Ta は γxy = 0.1∼1.5 付近で第一固有値 η a(1) と第二固有値 η a(2)が交差している．Fig. 3.2.1 下部の η a(1) のグラフにてこの点でグラフが折れ曲がって見えたのはこのためである．それぞれの連続した曲線は同じ固有ベクトルすなわち
変形モードを表している．したがって，無ひずみのとき 2 番目に小さい固有値の変形モードが，せん断下で η a(1) <0 となった固有値の固有ベクトルとなっている．な
お W は初期ひずみ部を拡大しても η a(1) と η a(2) は入れ替わっていない．Table 3.2.1に η a(1) が負になった時の固有ベクトルの 6 つの成分 (∆εxx ∼∆γxy ) を示す．固有ベ
クトルは大きさが 1 となるように規格化されているため，数値そのものに物理的な意味はない．いずれも垂直方向 (∆ε) は僅かでせん断方向である ∆γxy に対する不安定モードである．

16

Fig. 3.2.1

Stress-strain curves and change in the 1st eigenvalue under [111](1̄1̄2)
static shear (bcc).

17

Fig. 3.2.2

Stress free configuration at γxy =0.34 under [111](1̄1̄2) shear (bcc Mo).
The orientaion is equivalent to x[100]，y[010], z[001] with elongatedlattice length.

Table 3.2.1

Components of eigenvector at the first point of η a(1) <0 in bcc under[111](1̄1̄2) shear.
∆εxx ∆εyy ∆εzz ∆γyz ∆γzx ∆γxy
Fe -0.028 0.144 -0.121
0
0
0.982
Mo 0.078 0.242 -0.082
0
0
0.964
Ta -0.163 0.059 0.007
0
0
0.998
W 0.0789 0.246 -0.095
0
0
0.961

18

Fig. 3.2.3

Change in the eigenvalues (η a(1) ∼η a(4) ) under static shear (bcc(Fe, Mo,
Ta, W)).

19

3.3

応力-ひずみ曲線および固有値変化

Fig. 3.3.1 に MD 計算で得られた応力-ひずみ, ならびに全ての原子の第一固有値
a(1)

a(1)
η α(1) の平均値 ηavg
，最小値 ηmin の変化を示す．左目盛が応力，右目盛りが固有値である．応力ピークひずみを縦の破線で，固有値が 0 となる水平位置を横の破線で示している．セル中に不均一変形を生じるまでは，各原子の固有値は 0.1K の熱揺らa(1)ぎ下では差がないので平均値 ηavgと最小値はほぼ同一線上にある．せん断によって
a(1)

a(1)
ηavg
および ηmin は減少する．Fe, Mo, Ta には η α(1) の折れ曲がりが見られるが，これは静力学計算で示したように第一固有値の変形モードが他のモードとスイッチしたことによる．応力ピーク付近を拡大した図を図中に丸で囲って示した．Fe 以外は
a(1)

a(1)
応力ピーク前に ηavg
<0 となっており，系全体の不安定状態に達している．ηmin に着目するといずれも応力ピーク直前に急落する前駆的な挙動を示している．その後，後述のようにセル内に局所不均一変形を生じると個々の原子の力学状態は安定化さ
a(1)れ，ηavgが大きな正値を示す．

20

Stress
10

Stress(1)
Average of eigenvalue, ηαavg(1)
Minimum of eigenvalue, ηαmin

14

4
-50

τ① (static)

12

50
Stress, τ , GPa

Stress, τ , GPa

0
6

✂
✁

Eigenvalue, ηα (1), GPaτ① (static)
8

100

✂✁

10
0

8
6

-50

Eigenvalue, ηα (1), GPa(1)
Average of eigenvalue, ηαavg(1)
Minimum of eigenvalue, ηαmin
50

4
2

-100
2
0.1
Strain, γ①

0.2

-100

0

0.05

Fe

Stress, τ , GPa

8

✂✁

-150
0.25

(1)
Average of eigenvalue, ηαavg(1)
Minimum of eigenvalue, ηαmin

Stress

50

0
6
-50
4

τxy(static)

150

15
100
Stress, τxy, GPaτ① (static)

Eigenvalue, ηα (1), GPa

10

0.2

Mo

(1)
Average of eigenvalue, ηαavg(1)
Minimum of eigenvalue, ηαmin

Stress0.1
0.15
Strain, γ①

50
10
0
-50

-100

2

5
-100
-150

0

Fig. 3.3.1

0.1

0.2

-150
0.3

Strain, γ①00.1Strain, γxy

Ta

W

0.2

x

Change in the stress, minimum and average of eigenvalue η a(1) under[111](1̄1̄2) shear (bcc).

21

Eigenvalue, ηα (1), GPa

0

3.4

bcc の局所変形発生の詳細

Fig. 3.4.1 に応力ピーク前後でセル内に双晶変形が発生し始めるときの Fe のセルの一部のスナップショット (右)，同じ位置を Central Symmetry Parameter（CSP)(22)より格子欠陥と判定された原子のみ表示したもの (中央)，η a(1) の正負で色分けしたもの (緑色が正，赤色が負)(左) を再掲する (1) ．前報で示したように，Fe は η a(1) <0の原子が帯状に出現し，最初に η a(1) <0 原子が出現した帯の位置に双晶欠陥が発生した．Fig. 3.4.2 は Mo についての図である．Fe と異なり，前節で述べたようにまずほぼ全ての原子が η a(1) <0 となっている (図 (a))，その後，η a(1) が回復して正負となる部分が出現し，やはり帯状の不安定バンドが形成される．(図 (b))，双晶はやはりこの不安定バンドに発生し (図 (c))，双晶が形成されると η a(1) が正値に回復する(図 (d))．Fig. 3.4.3 の W は Mo と同様の変化を示した．Fig. 3.4.4 の Ta は Mo, W とa(1)
同様に系の平均 ηavg
が負になることからほぼ全ての原子が η a(1) <0 を示す (図 (a))．

Fe, Mo, W とは異なり，斜め方向に不安定バンドが形成された (図 (b))．この方向は (1̄1̄1) 面である．その後，欠陥も (1̄1̄1) 面に沿って発生している (図 (c) 中)．欠陥が形成されると 1 原子層を残し η a(1) <0 から正値に回復する (図 (d) 左)．このとき，右の図でも結晶のずれは (1̄1̄1) 面に沿っていることが分かるが，さらにせん断が進
むと (1̄1̄2) 面の双晶欠陥となった (図 (e))．なお，前報 (1) で Fe でも x 軸方向のセル
長さを変えると斜方向に不安定バンドおよび欠陥が形成されることを報告している．

Ta についても，x 軸方向のセル長さを短くすると Mo, W と同様に zx 面に平行に不安定バンドが発生した．

22

η a(1) <0 atoms

Defect atoms by CSP

Snapshot of atoms(a)γxy =0.190(b)γxy =0.191(c)γxy =0.192

Fig. 3.4.1

Distribution of η a(1) <0 atoms and defect atoms judged by CSP (bcc
Fe).

23

η a(1) <0 atoms

Defect atoms by CSP

Snapshot of atoms(a)γxy =0.150(b)γxy =0.153(c)γxy =0.155(c)γxy =0.159

Fig. 3.4.2

Distribution of η a(1) <0 atoms and defect atoms judged by CSP (bccMo).
24

η a(1) <0 atoms

Defect atoms by CSP

Snapshot of atoms(a)γxy =0.149(b)γxy =0.151(c)γxy =0.152(d)γxy =0.153

Fig. 3.4.3

Distribution of η a(1) <0 atoms and defect atoms judged by CSP (bccW).
25

η a(1) <0 atoms

Defect atoms by CSP

Snapshot of atoms(a)γxy =0.204(b)γxy =0.206(c)γxy =0.208(d)γxy =0.209(e)γxy =0.270
Fig. 3.4.4

Distribution of η a(1) <0 atoms and defect atoms judged by CSP (bcc
Ta).

26

3.5

固有ベクトルFig. 3.5.1 は Mo の γxy =0.151 のセルの一部について，左図の負の原子位置を右図で〇で示し，固有ベクトルの成分から得た最大せん断方向を，固有値の大きさでスケーリングした矢印で表示したものである．Fig. 3.4.2 と同じ箇所であるが見やすさの為 y 軸上下方向を広く表示し，表示断面の不安定モードのみを表示している．図(a) の xy 断面では 2 種類のモードが見られ，(100) 面と (110) 面に沿うモードである．yz 断面の図 (b) では (1̄1̄2) 面，(11̄0) 面，(211) 面に沿った 3 つの方向が見られる．成分は {∆εxx , ∆εyy ,∆εzz , ∆γyz , ∆γzx , ∆γxy }={0.075∼0.114, -0.002∼0.001, ±0,±0.604, ±0.633, -0.495∼0.496} と静力学の ∆γxy 以外のモードが大きくなっていた．
Fig. 3.5.2 は Ta についての図である．図 (b) は図 (a) 中の破線で示した断面を表示している．図 (a) では不安定バンド方向に平行な固有モードが多い．図 (b) で大きい
モードは (1̄01) 面に沿っている．成分は {∆εxx , ∆εyy ,∆εzz , ∆γyz , ∆γzx , ∆γxy }={-

1.73∼-0.139, ±0, ±0.475, ±0.766, ±0.458} であり，こちらも ∆γxy 以外のせん断成
分モードが大きくなっている．Fig. 3.5.3 は Ta の x 軸方向のセル長さを短くしたセルで，zx 面に平行に不安定バンドが発生した時の不安定モードである．自由度が制限されたことにより，図 (a) の xy 断面は先の Fig. 3.5.2 に比べてモードがそろっており，同一 xz 面でのモードは同じである．yz 断面は自由度が変わっていないので図(b) では Fig. 3.5.2(b) とあまり変わらないように見える．ただし，固有ベクトルの成分は {∆εxx , ∆εyy ,∆εzz , ∆γyz , ∆γzx , ∆γxy }={-0.165∼0.147, -0.02∼0.01, -0.02∼0.04,-0.587∼0.586, ±0.791, ±0.466} であり，yz 面でのせん断成分が大きくなっている．

27

y
x
(a)xy planeyz(b)yz plane
Fig. 3.5.1

Deformation mode of η a(1) <0 atoms at γxy =0.151 (bcc Mo).

28

A'

y
A

x

(a)xy plane

A-A' cross sectionyz(b)yz plane
Fig. 3.5.2

Deformation mode of η a(1) <0 atoms at γxy =0.207 (bcc Ta).

29

y
x
(a)xy planeyz(b)yz plane
Fig. 3.5.3

Deformation mode of η a(1) <0 atoms at γxy =0.216 (bcc Ta, shorter cellin x-direction).

30

第4章

4.1

fcc 単結晶の [1̄1̄2](111) せん断
解析

シミュレーション条件

対象は Ag, Al, Au, Cu, Ni, Pb, Pd, Pt の 8 元素とした．用いた fcc の直方体単位格子は Fig. 4.1.1 に示すように x=[1̄1̄2]，y=[111]，z=[11̄0] の座標系で，fcc の双晶面とされる (111) 面を xz 面とした．単位格子を 14×46×26 並べたものをシミュレー
ションセルとし，他は前章の bcc と同条件で静力学解析および MD シミュレーショ
ンを行った．

A
A

B
A

Fig. 4.1.1

A B

B

A B
A B

B

fcc unit lattice for [1̄1̄2](111) simple shear.

31

4.2

静力学解析

Fig. 4.2.1 に静力学条件下での応力-ひずみ曲線と，第一固有値 η a(1) および第二固有値 η a(2) の変化を示す．fcc の積層欠陥形成の変形経路 ( a6 [2̄11] + a6 [1̄21]，a は格子定数) の解析 (23) では負の固有値が発生しないことが報告されているが，[1̄1̄2] 方向へのせん断ではいずれも γxy =0.126∼0.243 で η a(1) <0 となる．bcc と異なり，η a(1) <0 となる点は応力ピークと一致しないので，応力が急減する不安定点ではなくこの変形経路よりもエネルギーが低い経路の存在を表す「分岐点」と考えられる．すなわち，

γxy のせん断に対する抵抗を失った点ではなく，別の方向への変形が生じる点である．fcc では η a(2) も γxy =0.277∼0.430 で η a(2) <0 となった．各原子毎の η a(4) までの固有値の変化を Fig. 4.2.2 および Fig. 4.2.3 に示す．Al, Pt 以外の 6 元素は η a(1) がせん
断とともに減少し，最初に η a(1) <0 となるまでは曲線は互いに交差することはない．
一方，Al, Pt は η a(1) <0 となるまでに η a(1) , η a(2) の曲線が交差している．Table 4.2.1,

Table 4.2.2 に η a(1) , η a(2) が負になった時の固有ベクトルの成分 (∆εxx ∼∆γxy ) をそれ
ぞれ示す．η a(1) の主成分は ∆γyz 方向のせん断モードである．垂直ひずみ成分 (∆ε)
とせん断を加えている方向 (∆γxy ) の成分は 0 だった．一方，η a(2) は垂直ひずみ成分(∆ε) とせん断を加えているモード (∆γxy ) のみであり，せん断方向の変形モードである．

32

Table 4.2.1

Components of eigenvector at the first point of η a(1) <0 in fcc under[1̄1̄2](111) shear.
∆εxx ∆εyy ∆εzz ∆γyz
∆γzx ∆γxyAg0000.974 -0.228
0
Al
0
0
0
-0.978 -0.211
0
Au
0
0
0
0.966 -0.257
0
Cu
0
0
0
0.992 -0.123
0
Ni
0
0
0
0.953 -0.302
0
Pb
0
0
0
-0.941 0.339
0
Pd
0
0
0
0.977 -0.234
0
Pt
0
0
0
-0.974 0.227
0Table 4.2.2

Components of eigenvector at the first point of η a(2) <0 in fcc under[1̄1̄2](111) shear.
∆εxx
∆εyy
∆εzz ∆γyz ∆γzx ∆γxy
Ag 0.210 -0.363 -0.038
00-0.907
Al -0.304 0.307 0.223
0
0
0.874
Au -0.224 0.354
0
0
0
0.908
Cu 0.219 -0.349 -0.083
00-0.907
Ni 0.111 -0.137 -0.036
00-0.941
Pb 0.193 -0.354 0.003
00-0.915
Pd -0.226 0.346 0.043
0
0
0.910
Pt -0.152 0.196 0.039
0
0
0.968

33

Fig. 4.2.1

Stress-strain curves and change in the 1st and 2nd eigenvalues under
static shear (fcc).

34

Fig. 4.2.2

Changes in eigenvalues (η a(1) ∼η a(4) ) under [1̄1̄2] static shear (fcc(Ag,Al, Au, Cu)).

35

Fig. 4.2.3

Changes in eigenvalues (η a(1) ∼η a(4) ) under [1̄1̄2] static shear (fcc(Ni,Pb, Pd, Pt)).

36

4.3

応力-ひずみ曲線および固有値変化

Fig. 4.3.1∼Fig. 4.3.3 に MD 計算で得られた fcc 8 元素の応力-ひずみ, ならびに全て
a(1)

a(2)

の原子の第一固有値の最小値 ηmin および第二固有値の最小値 ηmin の変化を示す．左軸目盛が応力，右軸目盛が固有値である．左のグラフはひずみ 0∼0.25 の全体を表示
し，右側に応力ピーク付近の変化を拡大して示している．応力ピークひずみを縦の破線で，固有値が 0 となる水平位置を横の破線で示している．Fig. 4.3.3 の図 (b) では左軸 (応力) の最小値は 0 でない為，応力が 0 となる水平位置にも横の破線を示している．応力はピークを示し急減するまで静力学解析のそれと一致している．各図
a(1)

a(2)

a(1)

a(2)

の左側の全体図より，Pt のみ応力ピーク前に ηmin と ηmin が交差して固有ベクトルが入れかわっている．Al はひずみ 0.05 近傍で ηmin と ηmin の数値が接近しているが，図中に拡大して示したように交差はしていない．次に右の拡大図について議論する．
a(1)Ag, Au, Ni, Pb, Pd の 5 元素は応力ピーク前に前駆的に ηmin が急減している．一方a(1)Al, Cu, Pt の 3 元素は応力が急減しはじめた後に ηmin が変化 (いずれも応力ピーク
a(1)

a(1)

点の後 ηmin は僅かに上昇し，その後減少) しているので，これらの元素は ηmin <0 の原子の発生がトリガーとなって局所変形を生じたのではなく，先に局所変形を生じa(2)始めて原子の固有値が変化したと考えられる．左図に戻って固有値 ηmin に着目する
a(1)と，第一固有値 ηmin が応力ピーク前に前駆挙動を示した 5 元素と Pt が応力ピーク
a(2)

a(2)

の後の変形過程で ηmin <0 となった．Al と Cu は ηmin <0 を生じていない．次節では，a(1)ηmin が前駆挙動を示した Ag, Au, Ni, Pb, Pd の 5 元素を代表して Ag を，前駆挙動を示さない Al, Cu, Pt のうち Cu と Pt について議論する．

37

Stress, τxy , GPaηαmin(1)

4

10
3
0
2
τxy (MD)
x

-10

1

0

0.05

0.1
0.15
Strain, γxy

0.2

30

x

x

5

20

α (2)
ηmin

min

20

τxy (static)

τxy(MD)

Stress, τxy , GPa

(2)
ηαmin

5

6

min

x

Minimum of eigenvalue, η α, GPa

τxy (static)

4
10
α (1)

3

η min

0
2
-10

1

-20
0.25

0
0.15

-20
0.16

0.155
Strain, γxyxMinimum of eigenvalue, η α, GPa

30

6

x

(a)Ag
4

4
20

η min

10
2

ηαmin(2)10ηαmin(1)

0

3

min

α (1)

0.1
Strain, γxy

0.2

Stress, τxy , GPa

Stress, τxy , GPa

(2)
ηαmin

x

5
(2)
ηαmin

2

ηαmin(1)
0

1

-10

0

0.145

0.15

min

x

3

τxy(static)x

Minimum of eigenvalue, ηα , GPaxτxy (MD)10τxy (MD)

Minimum of eigenvalue, ηα , GPa

τxy(static)-5

Strain, γxy

x

x

(b)Al
6
20

Stress, τxy, GPa

5
4

10

(1)
ηαmin

3
0
2
τxy(MD)

-10

1

0

0.1
Strain, γxy

0.2

(2)
ηαmin

τxy(static)
20

5
τxy(MD)
Stress, τxy, GPa

τxy(static)

Minimum of eigenvalue, ηαmin(1), GPaηαmin(2)

4
10
3
(1)
ηαmin

2
0
1

-20

0

0.15

0.155
Strain, γxy(c)Au
Fig. 4.3.1

a(1)

a(2)

Change in the stress, minimum eigenvalue ηmin and ηmin under[1̄1̄2](111) shear (Ag, Al, Au).
38

-10

(1)
Minimum of eigenvalue, ηαmin
, GPa

6

30

α (2)
ηmin

8

τxy(static)

6
10
4

0
τxy(MD)

2

-10

0

0.1
Strain, γxy

0.2

20

min

min

Stress, τxy, GPa

τxy(MD)
6
(1)
ηαmin

10

4

0

2

-20

Minimum of eigenvalue, ηα , GPa

Minimum of eigenvalue, ηα , GPa

8
20

(1)
ηαmin

Stress, τxy, GPa30

ηαmin(2)τxy(static)

0
0.145

-10
0.155

0.15
Strain, γxy

(a)Cu
15

α (1)

η min

20
0
5
τxy(MD)

-20

20

10
α (1)

η min

0
5

-20

min

10

τxy(MD)min

Stress, τxy, GPa

40

Stress, τxy, GPa

τxy(static)40
τxy(static)

(2)
ηαmin

60
Minimum of eigenvalue, ηα , GPaηαmin(2)Minimum of eigenvalue, ηα , GPa

15

-40
0

0.1
Strain, γxy

0
0.155

0.2

0.16
Strain, γxy

0.165

(b)Ni
10

1.5

τxy(static)

η min

0.5

0
τxy(MD)0

0.1
Strain, γxy1τxy(MD)

1

min

5
α (1)

Stress, τxy, GPa

1

min

η min

τxy(static)

Minimum of eigenvalue, ηα , GPa

Stress, τxy, GPa

α (2)

2

ηαmin(2)

(1)
ηαmin

0

0.5
-1

-5
0.2

0
0.115

0.12
Strain, γxy

-2
0.125

(c)Pb
Fig. 4.3.2

a(1)

a(2)

Change in the stress, minimum eigenvalue ηmin and ηmin under[1̄1̄2](111) shear (Cu, Ni, Pb).
39

Minimum of eigenvalue, ηα , GPa

1.5

40

6

10
0

4

-10

τxy(MD)

2

-20
0

0.1
Strain, γxy

0.2

20

τxy(MD)
610

ηαmin(1)
4

0

2

-30

-10

0

0.145

0.15

min

Stress, τxy, GPa

η min

30
τxy(static)

8

min

20

α (1)

α (2)
ηmin

30
Minimum of eigenvalue, ηα , GPa

8Stress, τxy, GPa

τxy(static)

ηαmin(2)10Minimum of eigenvalue, ηα , GPa

10

-20

Strain, γxy

(a)Pd
60

10

40

8

τxy(static)20

40τxy(MD)-20

τxy(MD)

min

η min
6

2

Stress, τxy, GPa

α (1)

min

Stress, τxy, GPa

τxy(static)

Minimum of eigenvalue, ηα , GPa

η min
8

10

(2)
ηαmin

α (2)

6

8

ηαmin(1)
4

2
6

0

0
0

0.1
Strain, γxy

0.2

-40

0.19

0.195
Strain, γxy(b)PtFig. 4.3.3

a(1)

a(2)

Change in the stress, minimum eigenvalue ηmin and ηmin under[1̄1̄2](111) shear (Pd, Pt).

40

Minimum of eigenvalue, ηα , GPa

10

4.4

Ag の局所変形発生の詳細

4.4.1

欠陥発生と η a(1) <0 の分布

Fig. 4.5.1 に応力ピーク (γxy =0.15432) 前後でセル内に欠陥が発生し始めるときのAg のセルの一部のスナップショット (右)，同じ位置を CSP により格子欠陥と判定された原子のみを表示したもの (中央)，η a(1) の正負で色分したもの (緑色が正，赤色が負)(左) を示す．応力ピーク直前の γxy =0.154 では欠陥はもちろん η a(1) <0 原子も存在しないが，先の Fig. 4.3.1(a) 右では前駆的に η a(1) が急減している．応力ピーク
直後の γxy =0.155 では η a(1) <0 原子がバンド状に現れている (図 (b))．不安定バンドは双晶面である xz 面 ((111) 面) ではなくて，(221) 面に沿って発生している．この方向は積層欠陥の発生経路に対応するが，この時点では欠陥は確認できない．その後η a(1) <0 原子が増え，不安定バンドが拡大するとともに，CSP で欠陥と判定された原子が発生する (図 (c))．その後，図 (d)∼(e) で双晶面である xz 面 (111̄) に沿って双晶欠陥が形成された．欠陥が形成されると対応箇所の η a(1) は回復し正値となる．前章で bcc Ta が xz でない面に欠陥が現れ，双晶面に再配置される結果を得たが，fcc
でも η a(1) が前駆挙動を示した Ag, Au, Ni, Pb, Pd は同様の挙動を示した．なお，この 5 元素はセルを小さくすると次節で示す Cu のように η a(1) <0 が発生せず双晶が発
生した．

41

η a(1) <0 atoms

Defect atoms by CSP

Snapshot of atoms(a)γxy =0.154(b)γxy =0.155(c)γxy =0.156(d)γxy =0.157(e)γxy =0.175
Fig. 4.4.1

Distribution of η a(1) <0 atoms and defect atoms judged by CSP (fcc
Ag).

42

η a(2) <0 の分布

4.4.2

Fig. 4.4.2 に η a(2) <0 原子が生じたひずみにおけるセルの一部の同じ位置を η a(1) ，η a(2) のそれぞれを正負で色分けしたもの (緑色が正，赤色が負)，CSP により格子欠陥と判定された原子のみを表示したものを左から順に示す．η a(2) <0 原子は紙面垂直方向に 1 原子列，線状に発生した (図 (a)(b))．η a(2) <0 原子列の位置を丸で囲った数
字で示した．η a(2) <0 の原子は当然 η a(1) <0 であるが，η a(2) <0 の原子は双晶欠陥が形成されるときの余剰原子面に対応しており「転位」と考える．

η a(1) <0 atoms

η a(2) <0 atoms①①

Defect atoms by CSP①(a)γxy =0.166

②

②

②

(b)γxy =0.170
Fig. 4.4.2

Distribution of η a(1) <0, η a(2) <0 atoms and defect atoms judged by CSP(fcc Ag).

43

4.4.3

固有ベクトルFig. 4.4.3 は γxy =0.155 における不安定バンド (左図) について，赤い η a(1) <0 原子の位置を右図で〇で示すとともに，固有ベクトルの成分から求めた最大せん断方向を，固有値の大きさによってスケーリングした矢印で可視化した図である．yz 面の図は図 (a) 中に破線で示した断面を表示している．また右図でセル内すべての η a(1) <0
原子を重ねて表示すると見にくくなるので，左図に見える原子についてのみ表示している．xy 面での矢印の大きさ，方向にはバラつきがあるが，yz 断面では同一 z 軸上の原子は皆同じモードを示している．Fig. 4.4.3 右図を拡大しつつ，η a(1) が正の原
子位置も〇で示したものが Fig. 4.4.4 である．左の xy 断面には (112̄) 面を実線で示した．yz 断面のモードは同一方向を向いており，長い矢印は 5 原子面上の原子位置を向いている．

44yx(a)xy planeyz(b)yz plane
Fig. 4.4.3

Deformation mode of η a(1) <0 atoms at γxy =0.155 (fcc Ag).

45

A'

(112)

y

y

A

z

x
Fig. 4.4.4

Deformation mode of η a(1) <0 atoms at γxy =0.155 (fcc Ag)(zoom).

46

4.5

Cu の局所変形発生の詳細

4.5.1

欠陥発生と η a(1) <0 の分布

Fig. 4.5.1 は応力ピーク (γxy =0.14936) 前後のスナップショット，CSP による欠陥
原子 (中央)，η a(1) <0 原子の分布である．4.3 節で述べたように Cu は η a(1) <0 原子の発生より局所変形=欠陥生成が先行する．図 (a) → (c) で右のスナップショットではずれは確認できないが，CSP による判定から欠陥生成に向けた不均一なずれが発生
していることがわかる．その後 x 軸 ([1̄1̄2] 方向) に沿った双晶欠陥がスナップショッ
トからでも確認できる (図 (c) → (d))．CSP の欠陥原子の間の一原子面が η a(1) <0 となっており，積層欠陥と考えられる (図 (d))．すなわち，Cu では積層欠陥を生じるまでの局所変形は η a(1) の正負では検出できない．

47η a(1) <0 atoms

Defect atoms by CSP

Snapshot of atoms(a)γxy =0.150(b)γxy =0.151(c)γxy =0.152(d)γxy =0.157

Fig. 4.5.1

Distribution of η a(1) <0 atoms and defect atoms judged by CSP (fcc
Cu).
48

4.5.2

η a(1) と η a(2) の分布

Cu では η a(1) <0 原子の発生より前に局所変形が起こったため，η a(1) の正負ではなく正値の η a(1) , η a(2) の大小関係から議論する．Fig. 4.5.2 は Bija の第一固有値 η a(1) (上)，第二固有値 η a(2) (下) の大小で着色した図である．Fig. 4.5.1 で示した部分を図 (g)(h) に赤枠で示した．γxy = 0.148 ではカラーバーを非常に狭い範囲にしているので原子のゆらぎによるまだら模様を示しているがほぼ均一である．応力ピーク直前の γxy =0.149
で η a(1) , η a(2) ともにモデル上部にわずかに固有値が高い部分が発生する．γxy =0.150
は先の Fig. 4.5.1 の CSP では構造に差が認められないにも関わらず，η a(1) , η a(2) ともはっきりと xz 面毎に差を生じており，力学的にはすでに不均一となっていることがわかる．γxy = 0.150 の時点では第二固有値 η a(2) の図 (g) は η a(1) の大小図 (c) と完全に対応しているが，CSP で欠陥が確認された γxy =0.151 の図では η a(1) が回復した暖色部分は η a(2) が低い寒色となりコントラストが逆転している．

49

(a)γxy =0.148(b)γxy =0.149(c)γxy =0.150(d)γxy =0.151(g)γxy =0.150(h)γxy =0.151

η a(1)(e)γxy =0.148(f)γxy =0.149
η a(2)

Fig. 4.5.2

Snapshot of atoms colored by magnitude of η a(1) and η a(2) (fcc Cu).

50

4.5.3原子移動による変形モード評価

Fig. 4.5.3 は最初に格子欠陥が発生した部分についてさらに拡大して，原子位置を〇印で示し，原子の移動方向を移動量の大きさでスケーリングした矢印で示したものである．図 (a) は応力ピーク前の γxy = 0.148∼149 の間，図 (b) は応力ピークを挟む γxy = 0.149∼150 間の移動量と方向を，xy 面と yz 面の 2 方向から表示してある．
図中には一つの原子を中心にカットオフ半径を示している．図 (a) は局所変形発生前であるため原子の移動方向に差は見られず，xy 面と yz の両者ともベクトルは同一方向を向いている．xy 面の原子移動が大きいのはせん断によるアフィン変形による．一方，図 (b) は xy 面の移動方向は図 (a) と変わらないのに対し，yz 面の移動は(111) 原子層毎に交互に移動するモードが確認できる．図に示したカットオフ半径内に，隣接する上下のモードが含まれるため，原子弾性剛性係数の固有値では検出できなかった可能性がある．

51

rcut

rcut

y

y

z

x

(a)γxy = 0.148 ∼ 0.149

cut
rrcut

rrcutcut

y
y

y

z
z

x

(b)γxy = 0.149 ∼ 0.150
Fig. 4.5.3

Migration direction of atoms just before defect nucleation (fcc Cu).

52

4.6

Pt の局所変形発生の詳細

Fig. 4.6.1 に応力ピーク (γxy =0.19306) 前後の Pt のセル全体の CSP による格子欠陥原子と γxy =0.250 のスナップショット (η a(1) の正負で着色) を示す．図 (c) 中の四角い枠は後で拡大して示す箇所である．Pt は応力ピーク時に η a(1) の最小値が負になっておらず，局所変形の発生過程で η a(1) <0 の原子が現れる．CSP による判定では，応力ピーク後にモデルの中央下部に欠陥と判定された原子が現れる (図 (b))．複数の原子が欠陥と判定されているが，この時点では双晶面である xz 面に沿って分布しているようには見られない．その後，y 軸方向に拡大して欠陥がブロック状に分布した (図 (c))．静力学解析では Al, Pt は固有モードの入れ替わりが生じるが，それによる相変態的な挙動を生じている可能性がある．その後，上部・下部にも新たな欠陥が確認され欠陥が成長したが，セル全体が欠陥判定されることはなかった (図(d))．この間にも η a(1) <0 の原子は発生していない．図 (e) に示すようにその後もせん断を継続するとブロック状から xz 面に沿った面欠陥に変化し，多数の双晶欠陥または積層欠陥となった．η a(1) <0 の原子は図 (e)(f) に矢印で示した積層欠陥のキンク部分に認められた．Fig. 4.6.2 に初期緩和終了時と応力急減直前の原子配置を最大級に拡大したスナップショットを示す．図ではせん断による [110], [111], [112] の角度変化を示している．応力ピーク直前では y 軸方向に対して [110], [112] 方向のなす角度がほぼ等しくなっている．Fig. 4.6.3 は Fig. 4.6.1(c) で示した枠線部分の原子配置について, 応力降下時と CSP の欠陥原子が分布した時の原子移動方向と量を矢印で示したものである．図 (a) に 1 原子に着目したカットオフ距離も重ねて示している．
右図において xz 面の層毎に原子が交互に移動する様子が確認できる．Cu と同様に
カットオフ半径内に隣接する上下のモードが含まれるため原子弾性剛性係数が平均化され，η a(1) <0 とならなかった可能性がある．Cu の Fig. 4.5.3 のように yz 面で互い違いに原子が移動しているが，Cu に比べて y 方向の移動が少なく，欠陥と判定さ53

れた図 (b) では y 軸方向の成分がなくなり xz 面上の移動のみとなる．

(a)γxy =0.193(d)γxy =0.196
Fig. 4.6.1(b)γxy =0.194(c)γxy =0.195(e)γxy =0.250(CSP) (f)γxy =0.250(η a(1) )Distribution of η a(1) <0 atoms and defect atoms judged by CSP (fcc Pt).

54

[111]
[112]

[110]

<112>

<110>

27.70°

28.16°

y35.27°

y19.48°

x19.48°

x(a)γxy =0.00

Fig. 4.6.2

35.27°(b)γxy =0.193

Snapshot of atoms and angle of atomic planes under [1̄1̄2](111) shear(fcc Pt).

55

rcut

rcut

y

y

z

x
(a)γxy =0.193

y

y

z

x
(b)γxy =0.195
Fig. 4.6.3

Deformation mode of atom during defect nucleation (fcc Pt).

56

第5章

5.1

hcp 単結晶の [1̄012](1̄011) せん断解析

シミュレーション条件

Mg, Co, Ti, Zr の 4 元素を対象とした．hcp の双晶面は (1̄011) 面で Fig. 5.1.1 に示すように x=[1̄012]，y=[101̄2]，z=[12̄10] とする斜方座標の単位格子を用いて解析を行った．軸比 c/a および軸比によって変化する x, y 方向の斜方セル長さ l と角度 θ はGEAM ポテンシャルでポテンシャルエネルギーが最小となる値とした (Table 5.1.1)．
総原子数は 100 800 個である．斜方座標で周期境界を全方向に適用し，3,4 章と逆向
きにせん断を与える．その他は bcc, fcc と同条件で静力学解析および MD シミュレー
ションを行った．

57

Fig. 5.1.1

Table 5.1.1

hcp unit lattice for [1̄012] (1̄011) shear.

Axial ratio c/a, cell length l , cell angle θ [deg.].
Mg
Co
Ti
Zr
c/a 1.628
1.621
1.624
1.624
l 3.688a 3.676a 3.680a 3.680a
θ 56.03
56.22
56.15
56.15

58

5.2

静力学解析

Fig. 5.2.1 に hcp の静力学計算における応力-ひずみ曲線と，第一固有値 η a(1) の変化を示す．γxy =0.307∼0.339 で 4 元素とも η a(1) <0 となった．η a(1) <0 となる点は応力ピークより前であり，fcc と同様に低エネルギーの変形経路への分岐点と考えられる．Fig. 5.2.2 は各元素における η a(1) ∼η a(4) までの変化を示した図である．無ひずみでは η a(2) と η a(3) の値がほぼ同じで，Mg はせん断開始直後に入れ替わっている．負となった固有値のモードは，Mg は ∆γxy =0 における η a(2) のモード (固有ベクトル)，

Co, Ti, Zr は η a(3) のモードである．Table 5.2.1∼Table 5.2.3 に無負荷での第三固有値までの固有ベクトル成分を示す．第一固有値は同じ z 方向のせん断に関するモードである．第二固有値は Mg のみ ∆γxy のせん断＋垂直方向を含むモードで，他は第一固有値のモードの ∆γyz , ∆γzx 成分を入れ替えたモードである．逆に第三固有値はMg のみ ∆γyz と ∆γzx のせん断で他の元素は Mg の第二固有値のモードである．せん
断ひずみ下で最初に負となったはいずれもこの ∆γxy ＋垂直ひずみのモードであり，垂直ひずみに他の経路への変形が加わることで，せん断方向への不安定を生じる．

59

30

hcp

Stress, τ , GPa

Co
20
Ti

①

Zr
10
Mg

0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

Strain, γ✁✂80

1st eigenvalue, ηα , GPa

60
Co
40
Ti
Zr

20
Mg
0
-20

hcp
-40

Fig. 5.2.1

0

0.1

0.2
0.3
Strain, γ✄☎Stress-strain curves and change in the 1st eigenvalue under static[1̄012](1̄011) shear (hcp).

60

200

50

Mg

Co

40

150

30
η

Eigenvalue, ηα , GPa

Eigenvalue, ηα , GPaηα (4)α (3)20ηα (2)10ηα (3)

0

0

0.1

0.2
0.3
Strain, γxyηα (2)50
ηα (1)

0.4

-50

0.5

0

0.1

0.2
0.3
Strain, γxy

0.4

0.5

0.4

0.5

100

Ti

Zr

ηα (4)

ηα (4)8080
ηα (3)

60
40

Eigenvalue, ηα , GPa

Eigenvalue, ηα , GPa

100ηα (1)100

ηα (2)ηα (1)20
0-20

ηα (3)0ηα (2)-10ηα (4)

ηα (3)6040

ηα (2)ηα (1)

20
0

0

0.1

Fig. 5.2.2

0.2
0.3
Strain, γxy

0.4

-20

0.5

0

0.1

0.2
0.3
Strain, γxy

Changes in eigenvalues (η a(1) ∼η a(4) ) under static [1̄012](1̄011) shear.

61

Table 5.2.1

Components of eigenvector for η a(1) of perfect hcp lattice in the
x=[1̄012]，y ′ =[101̄2]，z=[12̄10] coordination.
∆εxx ∆εyy ∆εzz ∆γyz
∆γzx ∆γxyMg
0000.470 0.883
0
Co
0
0
0
0.471 0.882
0
Ti
0
0
0
-0.471 -0.882
0
Zr
0
0
0
0.471 0.882
0Table 5.2.2

Components of eigenvector for η a(2) of perfect hcp lattice in the
x=[1̄012]，y ′ =[101̄2]，z=[12̄10] coordination.
∆εxx ∆εyy ∆εzz ∆γyz ∆γzx ∆γxy
Mg -0.266 0.185 0.067
0
0
0.944
Co
0
0
0
0.882 -0.471
0
Ti
0
0
0
0.882 -0.471
0
Zr
0
0
0
0.882 -0.471
0Table 5.2.3

Components of eigenvector for η a(3) of perfect hcp lattice in the
x=[1̄012]，y ′ =[101̄2]，z=[12̄10] coordination.
∆εxx ∆εyy ∆εzz ∆γyz ∆γzx ∆γxyMg
0000.883 -0.470
0Co -0.253 0.127 0.092
0
0
0.955
Ti -0.258 0.140 0.880
0
0
0.952
Zr -0.272 0.161 0.087
0
0
0.945

62

5.3

応力-ひずみ曲線および固有値変化

MD での応力-ひずみ曲線および第一固有値の最小値 ηmin の変化を Fig. 5.3.1，Fig. 5.3.2a(1)に示す．bcc, fcc に比べて MD の応力が静力学より小さい値を示して静力学の破線からずれており，特に Zr ではその差が著しいが内部構造緩和が原因と考える．応力ピーク近傍の拡大図 (右) を見ると，いずれも応力ピーク直前に ηmin が急落する前駆a(1)挙動が確認できるが，負になるのはピークを過ぎてから，すなわち局所変形が生じてからである．変形過程を観察した結果，いずれも同様の応答であったため，次節では代表して Mg について議論する．

63

20

4

20
τ (static)τ (static)

Stress, τ , GPa

2

0

1

-10

τ (MD)x0

0.1
Strain, γ

0.2

3

min

10

τ (MD)

min

η min

3

x

x

Stress, τ , GPa

α (1)Minimum of eigenvalue, ηα , GPax102ηαmin(1)
0
1

0
0.155

-20

Minimum of eigenvalue, ηα , GPa

4

-10
0.165

0.16
Strain, γ

x

x

(a)Mg

15

60

ηαmin(1)1030τxy(static)
x

0
5
-30

τxy (MD)x

min

Stress, τxy , GPa

Stress, τxy , GPa

min

60

10
30
5
τxy (MD)
x

0
-60
0

0.1
Strain, γxyτxy (static)

0
0.14

0.2

x

0.145
Strain, γ xy
x

x

(b)Co
Fig. 5.3.1a(1)Change in the stress, minimum eigenvalue ηmin under [1̄012] (1̄011)
shear (Mg, Co).

64

Minimum of eigenvalue, ηα , GPaη minMinimum of eigenvalue, ηα , GPa

15
α (1)

0.15

10

min

40

20

x

4
0
τxy(MD)2

0

0.1
Strain, γxy-20

ηαmin(1)20xτxy (MD)
6

-40

0.2

τxy (static)

8

Stress, τxy , GPa

6Stress, τxy, GPa

40

min

ηmin

Minimum of eigenvalue, ηα , GPa

τxy (static)α (1)Minimum of eigenvalue, ηα , GPa

8

0

4

0.145

0.15

-20

Strain, γxy

(a)Ti

40

40

τ① (static)✂✁

02
τ① (MD)0

0.1
Strain, γ①

✂✁

0
2

-20

0
0.145

0.2

0.15
Strain, γ①(b)ZrFig. 5.3.2a(1)Change in the stress, minimum eigenvalue ηmin under [1̄012] (1̄011)
shear (Ti, Zr).

65

min

min

-20

0

20

(1)
ηαmin

4

Minimum of eigenvalue, ηα , GPa

4Stress, τ , GPa

Stress, τ , GPa

20

6

Minimum of eigenvalue, ηα , GPaηαmin(1)6τ① (MD)5.4

Mg の局所変形発生の詳細

Mg の応力ピーク (γxy =0.15959) 前後の変化について，CSP では hcp 中の欠陥判定が不明瞭であるため，原子のエネルギー E a で着色して Fig. 5.4.1 に示す．原子エネルギーは高い方が暖色である．図 (a) は η a(1) の最小値が前駆的に急落するひずみの図であるが，次の図 (b) にて帯状にエネルギーの高い原子が斜方向の (202̄1) 面に沿って現れる．図 (a) 中央をよく見るとわずかに同じ位置に明るい部分を確認できる．図 (c) で η a(1) <0 の原子が高エネルギー帯に発生し，E a の差も明瞭になる．右の灰色のスナップショットでもずれが認められるようになる．したがってこの時点
では η a(1) <0 の原子は欠陥原子である．図 (d) で右のスナップショットのずれが明瞭になり，その部分の E a は著しく高い赤を示すが，エネルギーが低い青色と交互になっており，左の η a(1) も負から回復している．

66

η a(1) <0 atoms

Distribution of energy

Snapshot of atomseVeV(a)γxy =0.160
eVeV(b)γxy =0.161
eVeV(c)γxy =0.162
eVeV(d)γxy =0.163
Fig. 5.4.1

Distribution of η a(1) <0 atoms and atom energy E a (hcp Mg).

67

Fig. 5.4.1 では η a(1) <0 の原子が最初に発現したセルの一部を示したが，Fig. 5.4.2にモデル全体の変化を E a で着色して示す．図 (b)(c) の下側の赤い帯が Fig. 5.4.1 で説明した変化で発生した欠陥である．モデル全てで帯状の欠陥が発生したわけではなく，帯状の欠陥が発生した範囲を領域 I，欠陥が確認できなかった部分を領域 II とする (図 (b))．(b) → (c) で領域 II の上下に高エネルギーの欠陥が残り，それ以外の
斜方向に帯状に分布した欠陥が消滅した (図 (d)(e))．Fig. 5.4.2 図 (e) 中の黒枠で示
す領域 I, II の境界面の拡大図を Fig. 5.4.3 に示す．図中の実線はそれぞれ ⟨12̄10⟩ であり，結晶回転により，傾角粒界が形成された．

68

(a)

0.160 (b)Fig. 5.4.2

0.170

(c)

0.180

(d)

0.190

(e)

0.200

Distribution of E a in whole cell at γxy =0.160∼0.200 (hcp Mg).

I

65.74°
58.86°

II

I
Fig. 5.4.3

Snapshot of boundary surface at γxy =0.200 (hcp Mg).

69

Fig. 5.4.4 は η a(1) <0 の原子が発生した γxy =0.162 の図 (Fig. 5.4.1) について，η a(1) <0の原子の固有ベクトルの最大せん断方向を示したものである．矢印の大きさは負の η a(1) の固有値の大きさを示す．静力学での不安定モードは ∆γxy に対するものであったが，不均一変形を生じている途中の MD では {∆εxx , ∆εyy ,∆εzz , ∆γyz , ∆γzx ,∆γxy }={-0.265∼0.257, -0.200∼0.278, -0.248∼0.250, 0.171∼0.481, -0.706∼0.778, 0.434∼0.500}
であり，∆γzx が最も大きい．また，右図で I, II で示した原子の矢印を左図で倍の長さに拡大して重ねて表示した．図 (a) の xy 断面に示した細い実線はすべり面である(0001) 面で，不安定モードはこの面に沿ったものと考える．

70

A'(0001)
Ⅰ

Ⅰ

y' y
x

A
(a)xy plane

A-A' cross section

Ⅱ

Ⅱ

y'
z

(b)yz plane
Fig. 5.4.4

Deformation mode of η a(1) <0 atoms at γxy =0.162 (hcp Mg).

71

第 6 章 結言
双晶面に平行に単純せん断を与える分子動力学シミュレーションを Zhou の GEAM
ポテンシャルで bcc 4 元素，fcc 8 元素，hcp 4 元素について行った．均一なせん断から双晶や積層欠陥など不均一変形を生じて応力が低下する点について，原子弾性
剛性係数 (AES)Bija = ∆σia /∆εj の固有値 η a(1) と固有ベクトルの観点から検討した．

3 章では bcc 4 元素 (Fe, Mo, Ta, W) について (1̄1̄2) 面 [111] 方向のせん断解析を行った．その結果，

1. 全ての原子を格子点に固定した静力学解析では，第一固有値 η a(1) のみが負とな
りそのひずみ点は応力ピークすなわちせん断方向への変形抵抗喪失点であった．
対応する固有ベクトルの成分もせん断方向である ∆γxy に対するものである．

2. MD シミュレーションでは，Fe は応力が急減する直前に η a(1) <0 の原子が現れてそこに双晶が形成されたが，他の Mo, Ta, W は応力ピーク前にほぼ全ての原子が η a(1) <0 となっていた．

3. Mo, W は 2. のあと η a(1) が一部回復して，Fe と同様に zx 面 (1̄1̄2) に沿った
η a(1) <0 の不安定バンドが形成され，そこに双晶変形を生じた．
4. Ta は 3. の不安定バンドが (1̄1̄2) 面ではなく斜め方向 ((1̄1̄1) 面) に形成された．
欠陥も初め (1̄1̄1) 面に沿って発生したが，せん断が進むと (1̄1̄2) 面の双晶欠陥となった．なお，x 軸方向のセル長さを短くすると 3. の Mo, W と同様に zx 面に平行に不安定バンドが発生した．

72

5. Mo, Ta について η a(1) <0 となった原子の固有ベクトル {∆εxx , ∆εyy ,∆εzz , ∆γyz ,∆γzx , ∆γxy } の主軸から最大せん断方向を算出し，不安定モードとして可視化した．

4 章では fcc 8 元素 (Ag, Al, Au, Cu, Ni, Pb, Pd, Pd, Pt) の (111) 面 [1̄1̄2] 方向のせん断解析を行った．その結果，

6. 静力学解析では η a(1) と η a(2) が負値を示した．負となる点の固有ベクトルは
η a(1) は ∆γyz 方向を主とするのモードで，η a(2) は垂直ひずみ成分 (∆ε) とせん
断を加えている方向の変形モード (∆γxy ) であった．

7. bcc と異なり η a(1) <0 となる点は応力ピークより前にあり，変形方向の抵抗喪
失点ではなく低いエネルギー経路の存在を表す「分岐点」と考えられる．

8. MD シミュレーションを行うと Ag, Au, Ni, Pb, Pd の 5 元素は応力ピーク直前a(1)に前駆的に η a(1) の最小値 ηmin が急減したが，Al, Cu, Pt の 3 元素は応力が急a(1)減しはじめた後に ηmin が変化した．

9. Ag, Cu, Pt について，局所変形を生じて応力が急減するときの内部構造変化と η a(1) <0 原子の分布，固有ベクトルや原子移動方向を議論した．

10. Ag は 4. の Ta のように斜め方向に η a(1) <0 の不安定バンドが発生し，その後
xz 面 (111) に沿った双晶が形成された．η a(2) <0 の原子列も，双晶欠陥が形成される途中の転位芯に出現した．

11. Cu では η a(1) <0 原子の発生より欠陥生成が先行した．欠陥は x 軸 ([1̄1̄2] 方向)に沿った積層欠陥である．欠陥生成直前の原子移動を確認した結果，yz 面で(111) 原子層毎に互い違いの方向に移動するモードが確認できた．このモード
は原子エネルギーのカットオフ半径内で相殺されてしまい，原子弾性剛性係数の固有値では検出できないと考える．

73

12. Pt も η a(1) <0 原子が出ずに，セル半分が欠陥と判定される相変態的な挙動を生じていた．最終的には (111) 面に沿った複数の積層欠陥・双晶欠陥となったが，原子移動を見ると 11. の Cu と同様に yz 面で交互の変形モードを生じていた．

5 章では hcp 4 元素 (Mg, Co, Ti, Zr) の (1̄011) 面 [1̄012] 方向のせん断解析を行った．その結果13. 静力学解析では 7. と同じく η a(1) <0 となる点は応力ピークより前の分岐点であった．元素によって η a(1) ∼η a(3) の 固有モードが入れ替わっているが，負となった η a(1) の固有ベクトルはいずれも ∆γxy ＋垂直ひずみのモードである．

14. MD シミュレーションを行うと，η a(1) <0 原子はやはり欠陥生成してから発生
a(1)したが，ηmin は応力急減前に前駆的に減少していた．

15. Mg の局所変形発生を原子のエネルギー E a で可視化した．まず斜め方向の(202̄1) 面に欠陥が発生し，その後，結晶回転を生じて，zx 面に平行に傾角粒界が形成された．

16. η a(1) <0 の原子は (202̄1) 面方向の欠陥発生時に出現した．対応する固有ベクト
ルの最大せん断方向は (0001) 面に沿っていた．

74

参考文献
(1) Yashiro, K., et al., Computation Material Science, Vol.183, 109804, (2020).
(2) Narita, N., 日本機械学会会報, Vol.24(12), pp.984-992, (1985).
(3) Sainath, G. and Choundhary, B. K., Computation Material Science, Vol.111, pp.406415, (2016).
(4) Wu, Q., et al., Journal of Engineering Materials and Technology, Vol.143, 031007, (2021).
(5) Zhang, R. F., et al, Philosophical Magazine Letters, Vol.91(12), pp.731740, (2011).
(6) Yamakov, V., Acta Materialia, Vol.151(14), pp.4135-4147, (2003).
(7) Wu, X. L. and Zhu, Y. T., Physical Review Letters, Vol.101, 025503, (2008).
(8) Qi, Y., et al., MRS Online Proceedings Library, Vol.554, 367, (1998).
(9) Shao-Lin, M., et al., Scripta Materialia, Vol.180(15), pp.40-44, (2020).
(10) Jian, Z., et al., Journal of Physics: Condensed Matter, Vol.34, 115401, (2022).
(11) Minrong, A., et al., Journal of Physics: Condensed Matter, Vol.29, 046201, (2020).
(12) Zhou, X. W., et al., Physical Review B, Vol.69, 144113, (2004).
(13) Zhou, X. W., et al., Physical Review B, Vol.69, 035402, (2004).
75

(14) Zhou, X. W., et al., Acta Materialia, Vol.49(19), pp.4005-4015, (2001).
(15) 上田顯, コンピュータシミュレーション, 朝倉書店, (1990).
(16) 洲之内治男, サイエンスライブラリ–理工系の数学=15 数値計算 サイエンス社, (1978).

(17) 香山正憲, 固体材料の電子状態の基礎, 材料学会勉強会資料, (1993).
(18) Nosé, S., The Journal of Chemical Physics, Vol.81(1), pp.511-519 (1984).
(19) Wallace, D. C., Thermodynamics of Crystals, Wiley, Newyork, (1972).
(20) Wang, J., et al., Physical Review Letters, Vol.71(25), pp.4182-4185, (1993).
(21) Wang, J., et al, Physical Review B, (17).
(22) Li, J., Modelling and Simulation in Materials Science and Engineering, Vol.11(2), pp.173178, (2003).
(23) Yashiro, K., AIP Advance, Vol.9(6), 065014, (2019).

76

謝辞
本研究を進めるにあたり，無知無学な著者に対し多くのご指導を頂いた屋代如月
教授に心より御礼申し上げます．研究に関して様々な視点からご助言を頂いた内藤圭史助教にも心から感謝致します．また共に研究を進め，切磋琢磨した小川晃氏，郷卓磨氏，広瀬領祐氏，藤村侑樹氏，村井裕也氏，これまで数多くのアドバイスをいただいた研究室の先輩方，後輩方に御礼申し上げます．ありがとうございました．

77