修 士 論 文(題 目)

き裂の進展・発生に関する
局所格子不安定性解析

平成 29 年度

岐阜大学大学院
工学研究科 博士前期課程
機械システム工学専攻

氏名 西川涼一郎

目次

第1章

緒言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

第2章

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

4

2.1

分子動力学法の概要. . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1.1

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

5

2.1.2

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

6

2.1.3

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

8

2.1.4

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

8

2.1.5

原子弾性係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1.6

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

2.1.7

局所格子不安定性における原子の変形モード . . . . . . . . . . 11

第3章

相互作用遮蔽シートを用いたき裂進展に関する検討 . . . . . . .

3.1

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

3.2

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

13

3.2.1

応力ひずみ曲線 . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.2.2

原子応力の分布 . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2.3

不安定解析と AES の固有値に関する検討 . . . . . . . . . . . . 21

第4章

ねじり粒界モデルに関する検討 . . . . . . . . . . . . . . . . . .

4.1

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

4.2

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

i30

4.2.1

Σ5 ねじり粒界の結果 . . . . . . . . . . . . . . . . . . . . . . . 32

4.2.2

Σ13 ねじり粒界の結果 . . . . . . . . . . . . . . . . . . . . . . 38

4.2.3

Σ17 ねじり粒界の結果 . . . . . . . . . . . . . . . . . . . . . . 44

第5章

円孔の応力集中部からのき裂発生に関する検討 . . . . . . . . . .

5.1

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

5.2

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

第6章

49

5.2.1

応力ひずみ線図および最小固有値の変化 . . . . . . . . . . . . 51

5.2.2

固有ベクトルによる検討 . . . . . . . . . . . . . . . . . . . . . 55

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

58

参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

1

第 1 章 緒言
車体の振動の解析に用いられる有限要素法や車体周りの気流の流れを解析する

CFD，新材料の開発には第一原理計算など，工業に使われる計算科学技術の発展は著しい．この科学技術計算を支えるのが CPU や計算高速化のための GPGPU に用いられる GPU などのハードウェアであるがその高速化・高機能化はナノレベルの微細加工技術によって実現されている．最小配線の寸法は数ナノレベルで微細化されており，有限要素法等に代表される連続体近似に基づくシミュレーションモデル
の適用ができない領域に達している．
ナノレベルの材料挙動をシミュレーションする手法には粒子法 (MPS) や分子動力学法 (MD) が挙げられる．粒子法は連続体を有限の粒子として離散化する方法であるが，物質の挙動を詳しくに理解するには物質の最小構成単位である原子を考える
必要がある．分子動力学法により材料の変形挙動を原子レベルから検討しようとする試みは古くから行われており，1982 年には弾性体とリンクさせた解析が行われている (13) ．また 2003 年には松本ら (3) が，アモルファス金属のき裂の進展挙動について連続体モデルと MD 法による原子モデルと比較した検討を行っている．こういったマルチスケールな解析だけでなく，連続体近似を用いずに MD のみでナノレベル
での破壊へ適応するプローチも行われており，須賀ら (1) は界面き裂に MD 法を適用した研究を 1991 年に発表している．破壊力学が適用困難な界面破壊を想定して，負荷モードの一致などから，MD 法が界面き裂の挙動解析に有効であることを示している．1997 年には中谷ら (2) が，bcc 結晶のき裂先端からの転位挙動について検討している．Po ら (25) はニッケル単結晶の片側き裂進展に対して MD シミュレーションを行い，Ni における (111) すべり挙動の観察や，き裂進展方向の温度により変化す1

ることや臨界応力の結晶方位依存性について報告している．また都留ら (29) は Al に対して 1 億原子の大規模なせん断 MD シミュレーションを行い，塑性変形の開始が結晶粒内の転位源の張り出しによるものであり，転位理論と MD で得られた臨界応力が近い値をとることを報告している．こういった MD は Fe, Si, Al, Ni(30) などの汎用的な金属以外にも行われており，Qing ら (26) は超硬材料 WC 対してナノインデ
ンテーションを行い，(0001)WC が押し込みにより {11̄01} 面に転位が射出され転位ループすることで押し込み反力が低下することを報告している．久保田ら (30) は純鉄の多結晶体に対して MD を適用し，変形速度が遅い場合はき裂の進展に対する抵抗が粒径が小さくなるほど下がることを報告している．
以上のように，変形・破壊現象を原子レベルから検討する試みは多数なされてきたが，問題となるのは観察した事象に対する定量的な考察である．原子間ポテンシャ
ルを用いた動的なシミュレーションは，多かれ少なかれ実際の材料の特性から外れた「モデル材料」にすぎない．原子間力を電子状態から精密に評価する第一原理計
算分子動力学法では扱える原子数はスーパーコンピュータを用いても 1 万個程度であり，さらに熱揺動の効果は時間の問題から考慮できない．これまでの研究の多く
は，原子シミュレーションの結果を応力拡大係数などの連続体理論と比較した議論を行っているが，連続体近似の有効性を追認するだけになりがちである．したがって連続体近似が適用できない事象を観察した時，定量的な議論をするより基準がな
い．
我々のグループでは，原子シミュレーションにおける定量的な評価基準として局所
格子不安定性に着目した検討を行ってきた (7)–(11) ．非線形弾性における応力ーひずみ関係を表す弾性剛性係数 Bij = (∆σi /∆εj ) を基に，原子弾性剛性係数 (Atomic Elastic

Stiﬃness;AES)Bijα = (∆σiα /∆εj ) を定義し，その正値性から局所の変形開始を議論す
る．局所不安定については，北村, 梅野ら (4), (5) が，原子構造体の全自由度 (N 原子な
ら全体の並進・回転自由度を覗いた M = 3N −6) について系のエネルギーの 2 次導関数（ヘッセ行列）を直接求め，M × M の固有値および固有ベクトルから局所不安定を議論する手法を提案しているが，マトリックスが膨大になるため，10 万原子の系を直2接解析することは難しい．一方，Bijα は 6×6 マトリックスなので，それを全ての原子について計算するのは MD における力計算とほぼ変わらず 100 万原子でも実行可能である．これまで AES によるき裂への適用は Si, Fe, Mg について行われてきた (7)–(11) ．
固有方程式 Bijα ∆εj = η α ∆εi の解 η α(k) (k = 1 ∼ 6) の第一固有値 η α(1) が負の原子の固有ベクトル {εi } から変形モードをひずみ空間 {∆εxx , ∆εyy , ∆εzz , ∆εyz , ∆εzx , ∆εxy }で示し，さらにはこのひずみ成分 ∆εij からなるひずみテンソルの主軸から変形モードを明らかにすることなどを試みている．本研究では，特にき裂の発生に着目した検討を bcc-Fe について行う．3 章では原子間力を人為的に遮へいすることでき裂開口させる理想き裂のシミュレーションを行った．4 章ではき裂発生源として CSL ね
じり粒界の引張りシミュレーションを行った．5 章では，応力集中部からのき裂発生シミュレーションとして円孔の引張りシミュレーションを行った．

3

第 2 章 分子動力学法
2.1

分子動力学法の概要

分子動力学法 (molecular dynamics method，略して MD 法) は，系を構成する各粒子についてニュートンの運動方程式mαd2 r α= Fα2dt(2.1.1)をたて，これを数値積分することにより粒子の軌跡を求める方法である (14) ．ここで，mα ，r α はそれぞれ原子 α の質量および位置ベクトルである．原子 α に作用す
る力 F α は，系のポテンシャルエネルギー Etot の各位置における空間勾配として次式により求められる．

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

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

(2.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)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.1.6)(2.1.7)となる．両式の和と差をとると，
2 F

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

α

α

(

α

(t)

mα

r α (t + ∆t) − r α (t − ∆t) = 2∆tv α (t) + O (∆t)3)(+ O (∆t)4)(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.1.1

EAM ポテンシャル

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

5

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

Etot =N∑

F (ρ̄α ) +

α

N ∑
N
1∑
ϕ(rαβ )
2 α β(̸=α)(2.1.12)

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

(2.1.13)β(̸=α)

で評価する．
本研究では，鉄のシミュレーションで広く用いられている FS(Finnis Sinclar) ポテ
ンシャル (17) を用いた．FS ポテンシャルでは，

Etot =

∑

α

E =

α

∑{ 1 ∑α2 β̸=α}ϕ(r ) − AF (ρ̄ )αβα(2.1.14)

で系のエネルギーを評価する．ここで

ρ̄α =

∑

ρ(rαβ )()

ρ(rαβ ) = (r − d)2 + β(r − d)3 /d(2.1.15)β̸=αであり，d, β, A はパラメータである．また埋め込みエネルギー F (ρ) は F (ρ) =√ρ

と単純な関数形であるのが特徴である．d, β, A の値を表 2.1.1 に示す．

2.1.2

速度スケーリング法

分子動力学における温度制御には時間をスケーリングする自由度を導入し連成させて解く Nose-Hoover 法 (20) と，ここで説明する速度スケーリング法が用いられる．

6

Table 2.1.1

Potential parameters of FS potential.
A[eV]

d [nm]

β1.828905 0.3569745

1.8

速度スケーリング法は，統計熱力学より導かれる次の原子の運動エネルギーと温度の関係を用いて，以下のように制御する．

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

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

viα0 =

3kB T0mα)0.5(2.1.17)

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

viα =

3kB T
mα

)0.5

(2.1.18)式 (2.1.17) と式 (2.1.18) より以下の式が得られる．
(

viα0
T0
=
α
vi
T

)0.5

(2.1.19)

つまり，系の温度を T から T0 にするには，式 (2.1.19) の右辺を現在の速度に掛けてやればよい．ただ，これだけでは原子配置に反映されないので，Verlet 法における座標更新の式 (2.1.10) を以下のように置き換える必要がある．
√
ααr (t + ∆t) = r (t) +(F α (t)T0r α (t) − r α (t − ∆t) + (∆t)2T
mα)(2.1.20)

平衡状態では，Nose-Hoover 法 (20) によって得られるカノニカルアンサンブルに一致することが示されている．

7

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

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

1
=
V

σij

(

∂F
∂ηij

)

1
, Cijkl =
V

(

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

Bijkl ≡

(

∆σij
∆εkl

)

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

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

2.1.4

原子応力

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

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

σij =1 ∂EtotV ∂ηij(2.1.23)8

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

r α = J r̄ α(2.1.24)と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子α と 原子 β の間の距離 rαβ には(rαβ )2 = r̄iαβ Gij r̄jαβ(2.1.25)

なる関係が成立する．ただし，Gij = Jki Jkj である．仮想変形の Lagrange ひずみテ
ンソル ηij は]1[
Gij − δij2

ηij =(2.1.26)

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

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

∂rαβ
∂ηij

=

r̄iαβ r̄jαβ
rαβ(2.1.28)

以上の関係を用いると FS ポテンシャルにおける応力は次式で評価される
(

σij

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

) αβ αβ
N ∑
N (
r r
1 ∑
1 ′ αβ′ α ′ αβ=ϕ (r ) − AF (ρ̄ )ρ (r ) i αβjV α β(̸=α) 2r

ここで，Ω = V /N とすれば，各原子位置における原子応力を

σijα =) αβ αβN (r r
1 ∑1 ′ αβ′ α ′ αβϕ (r ) − AF (ρ̄ )ρ (r ) i αβjΩ β̸=α 2r(2.1.29)と定義すると，系の応力は

σij

N
1 ∑
=
σijα
N α

(2.1.30)となる．

9

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

1 ∂ 2 Etot=
V ∂ηij ∂ηkl

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

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

{

′

α

αβ

riαβ rjαβ }{ ∑ ′ αβ rkαβ rlαβ }−Aρ (r ) αβρ (r ) αβF (ρ̄ )
r
r
α
β̸=α
β̸=α
N
∑

′′

α

{∑

′′

′

]

αβ

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

α
Cijkl

} αβ αβ αβ αβN {1 1∑ϕ′ (rαβ ) ri rj rk rl′′ αβ=
ϕ (r ) −Ω 2 β̸=αrαβ(rαβ )2} αβ αβ αβ αβρ′ (rαβ ) ri rj rk rl−A
F (ρ̄ ) ρ (r ) −rαβ(rαβ )2
β̸=α
N
∑

′

{

α

αβ

]

rαβ rαβ }{ ∑ ′ αβ rkαβ rlαβ }ρ (r ) i αβjρ (r ) αβ− AF (ρ̄ )
r
r
β̸=α
β̸=α
′′

α

{∑

′′

′

αβ

(2.1.33)

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

2.1.6

原子弾性剛性係数以上で定義した原子応力，弾性係数から，原子弾性剛性係数 (Atomic Elastic Stiﬀ-

ness; AES) は以下で評価できる．
α
α
α
α
α
δil − 2σijα δkl )/2
δjl + σjk
δik + σik
+ (σilα δjk + σjl= Cijkl
Bijkl

10(2.1.35)

α symααWang らによる提案 (22), (23) に従い，Voigt 対称性をもたせた Bijkl≡ (Bijkl+Blkji)/2を用いて安定性評価を行う．なお，ここでの i〜l は直交座標の指標 (1,2,3) もしくは
α(x，y，z) であり，Bijkl
は 4 階のテンソルであるが，Voigt 対称性を持たせた場合，弾性係数テンソルと同じく独立な成分は 21 個となり，xx，yy，zz，yz，zx，xy を 1-6とする Voigt 表記を用いれば 6 × 6 のマトリックスとして Bijα と表すことができる．
その場合の式 (2.1.35) の対応を表 2.1.2 に示す．

Table 2.1.2
α
C11
+ σ1α

σ1α + σ2α
2
α
+ σ2α
C22

α
C12
−

Atomic elastic stiﬃnes matrix in Voigt notation.
σ1α + σ3α2
σ α + σ3α
α
− 2
C23
2
α
+ σ3α
C33
α
C13
−

σ4α
2
σα
α
+ 4
C24
2
σα
α
+ 4
C34
2
σ α + σ3α
α
+ 2
C44
2
α
C14
−

σ5α
2
σα
α
− 5
C25
2
σα
α
+ 5
C35
2
σα
α
+ 6
C45
2
σ α + σ1α
α
+ 3
C55
2
α
C15
+

sym.

2.1.7

σ6α
2
σα
α
+ 6
C26
2
σα
α
− 6
C36
2
σα
α
+ 5
C46
2
σα
α
+ 4
C56
2
σ α + σ2α
α
+ 1
C66
2
α
C16
+

局所格子不安定性における原子の変形モード

AES の各固有値に対する固有ベクトル {∆εi } は








∆ε1
∆εxx




 ∆ε2  ∆εyy 












∆ε
∆ε
3
zz
α
α
α(n)
α(n) 
α(n) 
∆σi = Bij ∆εj == η∆εi = η=η ∆ε4  ∆εyz  ∆ε  ∆ε 
5 
xz 


∆ε6
∆εxy
と書ける．上記のように固有ベクトル {∆εi } は 6 次元のひずみ空間におけるベクトルであり，変形モードを想定するのは困難である．そのため，{∆εi } の成分 ∆εxx , ∆εyy ,∆εzz , ∆εyz , ∆εzx , ∆εxy を成分とする以下のひずみテンソル [∆εij ] の主軸を議論する．

11





∆εxx ∆εxy ∆εxz

[∆εij ] = ∆εyy ∆εyz 

sym.
∆εzz
η α(1) < 0 の原子の固有ベクトルの x1 , x2 , x3 について，[∆εij ] の最も大きい固有値に対応した x3 を原子の主変形モードとして議論した．Bijα , [∆εij ] の固有値・固有ベク
トル解析には数値計算ライブラリ LAPACK を使用した．

12

第 3 章 相互作用遮蔽シートを用いたき裂進展に関する検討
き裂のメカニズムを検討していくため，まず bcc-Fe 完全結晶の中に初期欠陥としてき裂を導入し，その挙動を検討した．ここでは bcc-Fe 完全結晶内に相互作用を遮蔽するシートを挿入し，曲率半径 ρ が無限小となる理想き裂を導入して，まずはき裂進展の様子を観察するとともに局所格子不安定性の観点からき裂進展時の様子を
考察した．

3.1

シミュレーション条件

図 3.1.1 左に模式的に示すように，平板状の周期セル中央に貫通き裂を配置した系を解析対象とする．き裂前縁方向とき裂面の表記で [001](110)，[112̄](111) の系を対象とし，それぞれの結晶方位は x, y, z=[110], [1̄10], [001] および x, y, z=[110],[111], [112̄] である．表 3.1.1 にセル長さ Lx , Ly , Lz および総原子数を示した．き裂上下の原子は互いに力を及ぼさないように制御することでシート状の理想き裂を導入した．き裂幅 2c は 0.3Lx (x 方向セル長さ) とした．熱揺動による影響を排除するた
め，温度は極低温 (0.1K) とし，速度スケーリング法により制御した．10 ps の緩和計算を行った後，x, z 方向のセル長さを固定して y 方向に引張ひずみを毎ステップ∆εyy = 1.0 × 10−7 増加させて引張シミュレーションを行った．数値積分の時間ス
テップは 1 fs であるのでひずみ速度に換算すると ε̇yy = 1.0 × 108 /s である．初期緩和を除いて，x，z 方向の応力制御は行っていない．

13

Table 3.1.1

Cell size and number of atoms.
(110)

［111］

Lx [nm]24.3
Ly [nm]32.4
Lz [nm]2.0
Atoms 127,337

24.3
29.8
2.1
122,157

Ly

model

2c=0.3Lx

y

Lx

Lz

x
z

Fig. 3.1.1

Simulation model.

14

3.2

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

3.2.1

応力ひずみ曲線

図 3.2.1 に [001](110) モデルの応力ひずみ曲線を示す．応力ひずみ曲線は引張り開始から線形に上昇し，εyy = 0.0415 で応力ピークに達して急激に応力が低下し破断した．

1st peak
εyy=0.0415

Stress, σyy [GPa]

10

5

0

0.02

0.04

0.06

Applied strain, εyyFig. 3.2.1

Stress-strain curve of [001](110) crack.

15図 3.2.2 に [112̄](111) モデルの応力ひずみ線図を示す．引張り開始とともに引張応力は増加していき εyy = 0.0403 で応力ピークとなった．その後は (110) き裂の応力
ひずみ曲線とは異なり，応力は急激はせずに徐々に減少していく形となった．

1st peak

εyy=0.0403

Stress, σyy [GPa]

10

5

0

0.02

0.04

0.06

Applied strain, εyyFig. 3.2.2

Stress-strain curves of [112̄](111) crack.

16

3.2.2原子応力の分布図 3.2.3 から図 3.2.6 は，それぞれ [001](110) と [112̄](111) の応力ピーク時におけαる引張方向の原子応力 σyyに関する図である．図 3.2.3 と図 3.2.5 は (a) が各原子をασyyで色付けしたスナップショットであり，色付けの閾値は-10GPa〜80GPa としている．(b) は格子欠陥を可視化するため可視化ソフト AtomEye の Central Symmetry

Parameter(CSP) にて欠陥と判定された原子を示している．(110) き裂ではき裂の上αα
下斜め方向に高い σyyが分布している．また (b) を見ると高い σyyをもつ領域は bcc

構造とは異なった格子欠陥であることが確認できる．(111) き裂は，き裂先端の応力は (110) のそれよりも低い．CSP による欠陥も殆どはき裂表面近傍にのみ存在す
る．図 3.2.4 と図 3.2.6 はき裂面上 (0.5Ly 地点) を y = 0 として，y = 0 から 2d(dαは FS ポテンシャルのカットオフ値) の範囲内に存在する原子の σyyをプロットしαた図であり，(a) がき裂先端にとった座標を，(b) が σyyをプロットしたグラフを示している．なお，図 (a) の赤い原子は後述する η α(1) < 0 の原子である．(b) は横軸αが原子の x 方向の座標であり，縦軸が原子応力 σyyの値を示している．同図の赤線は線形破壊力学による無限体中のき裂の応力分布に，周期条件によって隣接するき√][ √裂の応力場を重ね合わせた式 (11) σyy (r) = KIC 1/ 2πr + 1/ 2π(d − r) である．rはき裂先端からの x 方向の距離を示しており，d はき裂先端から隣接周期セルのき裂先端までの x 方向の距離を示している．破壊じん性値 KIC は等方弾性体を仮定
して KIC =√2Eγs /(1 − ν 2 ) から求めた．ヤング率 E とポアソン比 ν はそれぞれE = (C11 − C12 )(C11 + 2C12 )/(C11 + C12 )，ν = C12 /(C11 + C12 ) より算出した．(110)と (111) の表面エネルギー γs および得られた破壊じん性値 KIC を表 3.2.1 に示す．

(111) き裂では応力ピーク点での応力分布は破壊力学による予測によく一致する．一方，(110) き裂では線形破壊力学の式よりも高い応力を示しており，Griﬃth の理論よりも高い応力になってからき裂進展した事がわかる．(110) き裂では応力ピーク前αの分布も大きく変わっていることがにき裂から欠陥が生じており，き裂先端の σyy

この差の一因と考える．

17α(a) σyysnapshotFig. 3.2.3

Snapshots of atoms([001](110)crack).

(b) CSP view
andyy

defectsatthestress-strain

peak

80
[001](110)Stress, σyy [GPa]Crack tip

60

40

20

2d

0.5Ly

r

0
150

200

250

Position of atoms, r [Å]αdistribution(b) σyyα
distribution(a) model of σyy

Fig. 3.2.4

αon y = 0 crack plane at the stressDistribution of atomic stress σyy
strain peak ([001](110)crack).

18α(a) σyysnapshotFig. 3.2.5

Snapshots of atoms([112̄](111)crack).

(b) CSP view
andyy

defectsatthestress-strain

peak

80

Stress, σyy [GPa][112](111)

Crack tip

60

40

20

2d

0.5Ly

r

0
150

200

250

Position of atoms, r [Å]αdistribution(b) σyyα
distribution(a) model of σyy

Fig. 3.2.6

αon y = 0 crack plane at the stressDistribution of atomic stress σyy
strain peak ([112̄](111)crack).

19

Table 3.2.1

Young’s modulus, Poisson ratio, surface energies and calculated fracture toughness.
model[001](110) [112̄](111)
E [GPa]νγs [J/m2 ]√KIC [MPa m]

143.039
0.362
1.703
2.139
0.749
0.839

20

3.2.3不安定解析と AES の固有値に関する検討
α(1)

図 3.2.7 は (110) き裂における最小固有値 ηmin の変化であり，破線は η α(1) = 0 の境α(1)界である．ここで，ηmin の値は各ひずみにおいて全原子の中で最も小さかったものであり，ある原子の η α(1) をモニターし続けたものではないことに注意されたい．図中に示した (a)〜(j) の時点におけるスナップショットを η α(1) < 0 の原子を赤く着色して図3.2.8 に示した．[001](110) モデルでは引張り初期は系内に η α(1) < 0 となる原子は存在していない．引張り開始からすぐに最小固有値の値は低下して負の原子が現れる．
図 3.2.8(a) に示すように (a)εyy = 0.0031 で相互作用を遮蔽したセル中央部に η α(1) < 0となる原子が発生し，引張りに伴って (b)εyy = 0.0033,(c)εyy = 0.0055 の図のように相互作用を遮蔽した上下部で不安定原子が出現と消失を繰り返す．(d)εyy = 0.0074
以降，き裂先端にのみ η α(1) < 0 の原子が存在するように変化し，(e)εyy = 0.0167
α(1)で ηmin が最大の負値を示した．その後，(e) → (f) で η α(1) < 0 の原子は消失するが，(f)εyy = 0.0197 で再びき裂先端に η α(1) < 0 の原子が現れる．(g)εyy = 0.0219 になる
とき裂先端から円弧を描くように結晶内部に η α(1) < 0 の原子が発生し，欠陥の発生によってバタフライ状に η α(1) < 0 の原子が分布している (図 3.2.8(i))．(110) き裂では (i)εyy = 0.0415 の応力ピーク点の後，き裂が (110) 面を不安定に進展して破断する(図 3.2.8(j))．進展時の最小固有値の値はほぼ一定 (図 3.2.7，-200< η α(1) <-150GPa)
で，バタフライ形状も保たれている．完全に破断した後には η α(1) の値は正値に回復する ((110) 表面原子の値)．

21Minimum of eigen value, η a(1), [GPa]100η α(1)
(a)

0
(b)

-100

(c) (d)

(g)

(h)

-200
(j)
(i) stress peak

-300
(e)

-400

Fig. 3.2.7

(f)

0

0.01

0.02

0.03

0.04

Applied strain, εyy

0.05

α(1)

0.06

Change in the minimum 1st eigen value ηmin ([001](110)crack).

22

Fig. 3.2.8(a) εyy = 0.0031(b) εyy = 0.0033(c) εyy = 0.0055(d) εyy = 0.0074(e) εyy = 0.0167(f) εyy = 0.0197(g) εyy = 0.0219(h) εyy = 0.0267(i) εyy = 0.0415(j) εyy = 0.0465

Snapshots of atoms during tension in the [001](110)crack. Red circles
indicate η α(1) < 0 atoms.

23図 3.2.9 は (111) き裂における最小固有値の変化であり，図 3.2.10 はセル全体の変形の様子（ただし (a)〜(i) は図 3.2.9 とは対応していない）である．引張り開始からき裂が開口するが，η α(1) < 0 の原子が現れたのは εyy = 0.0206 であった．その後，応力ピークの点の図 3.2.10(c) でも η α(1) < 0 の原子は少なく，(c) → (d) でき裂右端は水平にわずかに進展し，左端も左下方向へわずかに進展する．(110) き裂とは異なり，き裂先端から転位が射出され，周期境界の下で交差した (図 3.2.10(e) および (f)の黒丸で囲った原子)．その後き裂先端が鈍化し (図 (g))，き裂の上下には先の (110)
き裂と同様のバタフライ状の η α(1) < 0 の原子が現れている．その後，き裂右端がジα(1)グザグに進展して，図 3.2.10(i) で x 周期境界まで達した．この間，図 3.2.9 の ηminは約-150GPa で一定である．

24

応力ピーク前後のき裂部分を拡大して表示したものが図 3.2.11 である．図の (a)〜α(1)(f) は，ηmin の図 3.2.9 に示した点に対応する．図では central symmetry parameterによる欠陥もあわせて示している．(110) と同じように引張り初期は不安定原子は系内に存在せず，εyy = 0.0206 でき裂両端に不安定原子が分布する．図 3.2.9 の応力α(1)ピーク点で ηmin が急減しているが，図 3.2.11(b) の CSP を見るとき裂先端の角部に欠陥が成長しはじめていることがわかる．図 (b) の η α(1) < 0 の分布を詳細に見ると，左の先端には η α(1) < 0 の原子が存在し，右はき裂面に η α(1) < 0 の原子が見られな
い．CSP の形状に違いは見られないが，次の図 (c) では左は斜め下に開口するとともに積層欠陥が上に発生している．その後き裂先端では積層欠陥が拡大し，右端からはまず右下方向に転位が射出される ((e)εyy = 0.0412)．さらに右端からは図 (f) のように下方にも転位が射出された．これらは (112) および (001) のすべり面である．
図 3.2.12 はき裂左端を拡大したものである．(112) 面に開口した際に先端から積層欠陥を生じている．右端と左端でこのような違いを生じたのは力の遮蔽シートが原子配置に対して完全に対称となっていなかったためと考える．

25

Minimum of eigen value, η a(1), [GPa]100η α(1)
0
(a)

(e)
(b)

-100

(c)
(d)

-200

(f)

-300

-400

Fig. 3.2.9

ε yy=0.0403
stress peak

0

0.02

0.04

0.06

Applied strain, εyyα(1)Change in the minimum 1st eigen value ηmin ([112̄](111)crack).

26

(a) εyy = 0.0000(b) εyy = 0.0206(c) εyy = 0.0403(d) εyy = 0.0406(e) εyy = 0.0415(f) εyy = 0.0420(g) εyy = 0.0555(h) εyy = 0.0621(i) εyy = 0.0669

Fig. 3.2.10

Deformation process during tension in [112̄](111)crack. Red circles
indicate η α(1) < 0 atoms.

27

Fig. 3.2.11(a) εyy = 0.0206(b) εyy = 0.0403(c) εyy = 0.0405(d) εyy = 0.0406(e) εyy = 0.0412(f) εyy = 0.0415

Deformation process before and after the stress peak in (111)crack.
Top: Red circles indicate η α(1) < 0 atoms. Bottom: CSP view.
28

(a) εyy = 0.0403(b) εyy = 0.0404(c) εyy = 0.0405(d) εyy = 0.0413
Fig. 3.2.12

Snapshots of crack propagation and blunting during tension after
stress peak in (111). Left: Red circles indicate η α(1) < 0 atoms. Right:
CSP view.

29

第 4 章 ねじり粒界モデルに関する検討
前章では，予め欠陥としてき裂を導入したモデルに対してシミュレーションを行っ
たが，き裂が発生する状態を議論するためには，き裂がない状態のモデルを検討する必要がある．今回はき裂の発生源としてねじり粒界を想定し，3 つの CSL(Consident

Site Lattice) ねじり粒界モデルを用いて粒界き裂の発生について議論した．

4.1

シミュレーション条件

図 4.1.1 にシミュレーションモデルを示す．図 4.1.1(a) に示すように，セルの上半分を z 軸 [001] 方向を軸として回転させて x，y 方向で周期性が満たされるサイズで切り出し，セル中央面にねじり粒界を作成した．回転角度は θ = 36.9◦ , θ = 22.6◦ , θ = 28.1◦の 3 つで，それぞれ Σ5, Σ13, Σ17 ねじり粒界となる．その後，図 4.1.1(b) に示すよ
うに上下周期境界上の粒界がセル内に入るように周期境界の下で 0.25Lz だけ z 方向に原子をシフトさせた．表 4.1.1 に系の寸法および総原子数を示す．全方向周期境界の下で x, y, z 方向の各応力が 0 となるように応力制御しながら 5, 000 [fs] の緩和計算を行った後，z 方向に引張ひずみを毎ステップ ∆εzz = 1.0 × 10−7 増加させて引張シ
ミュレーションを行った．数値積分の時間ステップは 1 [fs] であるのでひずみ速度に換算すると ε̇zz = 1.0 × 108 [sec−1 ] である．引張り中の x, y 方向のセル辺長は固定している．熱揺動による影響を排除するため，温度は極低温 (0.1[K]) とし速度スケー
リング法を用いて制御した．

30

1/4 Lz

Lz

Ly

Twin grain boundaryTwistz [001]0Lx

x [100](a)(b)

Fig. 4.1.1

Simulation model for twist grain boundaries.

Table 4.1.1

Size, number of atoms and rotation angle.

model

Σ5

Σ13

Σ17

Lx [nm]
8.6
11.2
9.7
Ly [nm]
8.6
11.2
9.7
Lz [nm]
25.8
33.8
29.2
Atoms 162,000 358,956 235,824
θ[deg] 36.8667 22.6167 28.0667

31

4.2

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

4.2.1

Σ5 ねじり粒界の結果
α(1)

図 4.2.1 に Σ5 モデルの応力ひずみ曲線，第 1 固有値の最小値 ηmin の変化と，緩和
終了後および主要な点におけるスナップショットを示す．スナップショット中の赤い原子は η α(1) < 0 の原子である．応力ひずみ曲線は εzz = 0.0973 と εzz = 0.1415 でα(1)ピークを示し，εzz = 0.1470 でほぼ破断した．ηmin の変化では，応力ひずみ曲線では大きな変化がない εzz < 0.0973 の領域で著しい変化を示している．緩和計算終了後の εzz = 0.0 において，右図に見られるように粒界部分に η α(1) < 0 の原子が存在する
ため，引張り前から負値 (−200GPa 程度) を示している．第一ピーク前の η α(1) < 0の変化については後述するが，右図のように最初のピークでは粒界近傍に η α(1) < 0の原子が多く現れ，第二ピークの図では粒内に縞状に η α(1) < 0 の原子が分布している．これは (112) の双晶境界で後で詳細に述べる．破断は粒界部分から生じた．
α(1)ηmin の変化の図に示した (a)〜(k) の点および破断後の (l)（範囲外）に対応する
η α(1) < 0 の原子の分布を図 4.2.2 に示す．図では η α(1) < 0 原子の分布を議論するために斜めから見た図や，η α(1) > 0 の原子も表示させた図などを混在させている．図のセル境界の赤，緑，紫色の直線はそれぞれ x, y, z 軸に対応する．緩和終了後の (a)では上下の粒界部分に η α(1) < 0 となっている不安定原子が層状に存在している．図α(1)4.2.1 の ηmin の変化において，まず (a) → (b) は η α(1) の値がほぼ一定で振動している．この時，図 (b) のように η α(1) < 0 の原子層が上下に広がる．(b) → (c) の変化では粒界中央部の原子の η α(1) が正に回復する．この時，粒界部は乱れているが bcc 構
造のままであった．これは CSL 粒界から数層分の原子が引張りを駆動力として原子α(1)の再配置が行われたためと考える．粒界緩和が完了した点が (d) で，ηmin が上昇している．(e) は η α(1) > 0 の原子も表示したものである．(d) → (f) では一度回復した
粒界部分に η α(1) < 0 の原子を生じ，再び (f) → (g) で粒界面の原子の η α(1) は回復し，粒界から離れた部分が二層の η α(1) < 0 となる構造になっている．その後は (g) → (i)α(1)と ηmin が減少するが，η α(1) < 0 の層の厚みが増大している．応力ピーク直前の (i)

32

α(1)

点で中央と下部結晶の内部に η α(1) < 0 の原子が現れ，図 4.2.1 の ηmin の変化でも不安定に大きく低下している．その後第一応力ピークで粒内に多数の変形双晶が発生
し，応力減少を生じる．図 4.2.3 は初期配置で {112} 面を抽出したものと，第 2 応力
ピーク時の yz 平面のスナップショットを示したものである．粒内の縞は 90◦ 回転させると一致する．調べた結果，(121) 面であるので第一ピーク後の欠陥は双晶変形によるものである．
図 4.2.4 はΣ 5 における第 2 応力ピーク後から破断までのスナップショットである．
上の全体図の，中央部の結晶粒と上の粒界近傍を見るとすべり面に沿った結晶内のη α(1) < 0 の縞が別のすべり面のそれと合体している．この後，下に拡大して示したように εzz = 0.1432 から上部粒界の η α(1) < 0 のすべり面と粒界が交差する部分からき裂が発生した．その後は，y の正方向と奥行き方向にき裂が進展し，εzz = 0.1455のように破断する．

33

εzz=0.0

εzz=0.0973
1st stress peak

εzz=0.1415
ε =0.1500
2nd stress peak zz
Fig. 4.2.1α(1)Stress-strain curve, change in the minimum 1st eigen value ηmin andsnapshots of fracture process in Σ5 twist grain boundaries.

34

(a) εzz = 0.0(b) εzz = 0.0270(c) εzz = 0.0383(d) εzz = 0.0411(e) εzz = 0.0411(f) εzz = 0.0668(g) εzz = 0.0682(h) εzz = 0.0815(i) εzz = 0.0960(j) εzz = 0.1043(k) εzz = 0.1415(l) εzz = 0.1470

Fig. 4.2.2

Snapshots of fracture process in Σ5 twist grain boundaries. Red circles
indicate η α(1) < 0 atoms.

35

θ=79.695
Twin plane(121)

θ

Tw
in p

lan
e

(11
2)

+90[deg]zy(a)Fig. 4.2.3(b)Snapshots of slip plane in Σ5 twist grain boundaries. (a) (112) slip
plane in initial structure, (b) η α(1) < 0 distribution after the 1st stress
peak.

36

εzz =0.1415
2nd peakεzz =0.1431

εzz =0.1431

εzz =0.1432

εzz =0.1433

εzz =0.1434

εzz =0.1435

εzz =0.1455zyFig. 4.2.4

Crack propagation after 2nd stress peakin (Σ5 twist grain boundary).

37

4.2.2

Σ13 ねじり粒界の結果
α(1)

図 4.2.5 に Σ13 モデルの応力ひずみ曲線，最小固有値 ηmin の変化とスナップショッ
トを示す．応力ひずみ曲線は Σ5 と同じく 2 つのピーク点を持っているが，最小固α(1)有値 ηmin の変化は Σ5 と比較して著しい変化はなく，振動しながら低下している．
スナップショットでは粒界部分に η α(1) < 0 原子が分布 (εzz = 0.0) し，応力ひずみ曲線の第 1 ピークで広がっているのは Σ5 と同じである．第一ピーク後に粒内に縞状の η α(1) < 0 の分布を生じるのも Σ5 と同じであるが，Σ5 と異なり Σ13 では縞は上下の粒内で平行となっている．第二ピーク後，εzz = 0.1530 で粒界から破断した．
α(1)

図 4.2.6 は図 4.2.5 の ηmin の (a)〜(l) に対応させて η α(1) < 0 の分布を見たものである．Σ5 と異なり，粒界部分の回復などはなく粒界近傍の η α(1) < 0 の層が増加する．
応力ピーク直前の η α(1) < 0 の分布 (e) と (f) の CSP 表示を見ると欠陥だけでなく上下の原子層も η α(1) < 0 となっている．応力ピーク後の (g) → (h) では下部粒界から

η α(1) < 0 の集団が面状に成長している．Σ5 モデルでは，結晶内に η α(1) < 0 の原子が数個現れた後に粒内に一気に双晶が形成されたが，Σ13 では転位ループのように粒界から伝ぱしている．しかしながら，応力ひずみ曲線の谷に相当する εzz = 0.09994の (i),(j) 図では，縞状に不安定原子が分布している (j) に対して (k) の CSP で多くの原子が欠陥原子と判定されており，1 つの転位がすべり面を横断したような変形ではない．
図 4.2.7 は応力第 2 ピークにおける yz 平面のスナップショットを初期配置における {112} すべり面と比較して示したものである．η α(1) < 0 の縞の向きは Σ5 と異なり，{110} すべり面に一致する．上下で模様の差が少ないのはねじり角 θ が小さいためと考えられる．
図 4.2.8 にΣ 13 におけるき裂進展の様子を示す．左の全体スナップショットを見ると第 2 応力ピーク時 (εzz = 0.1435) に中央結晶内に発生した縞模様が εzz = 0.1495
では広い範囲で消失している．Σ5 と同様，粒界からの割れを生じることで，縞状に分布している (110) 面の η α(1) < 0 が回復したものと考えられる．Σ5 のそれより回38復範囲が大きいことと，先の図 4.2.6(k) の CSP から面間隔が大きくひずんだ部分がη α(1) < 0 の縞状の構造と推測する．下の拡大図を見ると斜めの η α(1) < 0 の縞によって粒界がジグザグに波打った状態となっている．εzz = 0.1496 で上部粒界の周期境
界上の上下の縞状の η α(1) < 0 が交差する粒界近傍からき裂が発生した．

39

Fig. 4.2.5

εzz=0.0

εzz=0.0909
1st stress peak

εzz=0.1435
2nd stress peak

εzz=0.1530
α(1)Stress-strain curve, change in the minimum 1st eigen value ηmin andsnapshots of fracture process in Σ13 twist grain boundaries.

40

(a) εzz = 0.0(b) εzz = 0.0316(c) εzz = 0.0525(d) εzz = 0.0690(e) εzz = 0.0909(f) εzz = 0.909(CSP)(g) εzz = 0.0910(h) εzz = 0.0913(i) εzz = 0.0925(j) εzz = 0.0994(k) εzz = 0.994(CSP)(l) εzz = 0.1530

Fig. 4.2.6

Snapshots of fracture process in Σ13 twist grain boundaries. Red circles
indicate η α(1) < 0 atoms and snapshots by CSP.
41θ'

θ' = 69.964

θ

np

twi

θ = 80.049

lan
2}

1
e{1

slip plane{110}

z
y

(a)
Fig. 4.2.7(b)Snapshots of slip plane in Σ13 twist grain boundaries. (a)slip planemodel, (b) simulation model with yz view in 2nd stress peak.

42

εzz =0.1435
2nd peakεzz =0.1495

εzz =0.1495

εzz =0.1496

εzz =0.1497

εzz =0.1498

εzz =0.1499

εzz =0.1518

zyFig. 4.2.8

Crack propagation in Σ13 twist grain boundaries after 2nd stress peak.

43

4.2.3

Σ17 ねじり粒界の結果
α(1)

図 4.2.9 に Σ17 モデルの応力ひずみ曲線および最小固有値 ηmin の変化とスナップ
ショットを示す．やはり応力ひずみ曲線は大きな 2 つのピークを示した後破断したα(1)が，最小固有値 ηmin は Σ5 と同じく複雑な変化を第 1 ピーク前に示している．何より特徴的なのは第 1 応力ピークでは粒界部以外の原子が η α(1) < 0 となっていることである (右上スナップショット)．このように，結晶全体が η α(1) < 0 となるときは，転位のような局所的な変化よりも相変態のような集団的な変形を生じやすいことをα(1)

報告している (8) ．その後，第 1 ピーク直後に ηmin は急減し，粒内に Σ13 と同様の縞模様を残して回復する．
α(1)

図 4.2.10 に ηmin の変化の (a)〜(i) に対応させた図を示す．εzz = 0.0 では系内に

η α(1) < 0 の原子は存在しない．(b) → (c)，(d) → (e) などの変化は Σ5 の変化と似ており，粒界近傍の構造緩和と考えられる．(f) の時点ではまだ η α(1) < 0 の原子は粒界近傍のみだが，(g) の第 1 ピークでは図 4.2.9 に示したように粒内すべてが η α(1) < 0 となる．この時を CSP による可視化をしたのが (g) で粒内に欠陥原子は認められない．
その後の応力急減をもたらす変形によって粒内の多くの原子は回復して η α(1) > 0 とα(1)なるが，図 4.2.9 の ηmin の (f) → (g) → (h) の変化からわかるように縞の部分の η α(1)はより大きな負値となる．
図 4.2.11 は応力第 2 ピークにおける yz 平面のスナップショットである．同図 (a),(b) の青およびピンク色の線が {112} すべり面を yz 平面に投影したものである．同図 (b) より，ねじっていない中央結晶の縞は青線に一致せず，{110} すべり面を投影したに一致する．しかし，ねじった上下の結晶内は {112} すべり面を投影したピンク色の線に一致しており，上下の粒内で別の変形を生じたことが推測される．
図 4.2.12 は破断した下の粒界を側面から見たものである．Σ17 では図 4.2.9 の応力ーひずみ曲線のように小さな第 2 ピークの後，第 3 ピークで粒界をはさんで上下の η α(1) < 0 の縞が交差する部分からへき開を生じて破断した．

44

εzz = 0.0

εzz=0.1037
1st stress peak

εzz=0.1535
3rd stress peakFig. 4.2.9

εzz=0.1595α(1)Stress-strain curve, change in the minimum 1st eigen value ηmin andsnapshots of fracture process in Σ17 twist grain boundaries.

45

Fig. 4.2.10(a) εzz = 0.0(b) εzz = 0.0400(c) εzz = 0.0647(d) εzz = 0.0649(e) εzz = 0.0870(f) εzz = 0.0965(g) εzz = 0.1037(CSP)(h) εzz = 0.1295(i) εzz = 0.1580

Snapshots of fracture process in Σ17 twist grain boundaries. Red
circles indicate η α(1) < 0 atoms and snapshots by CSP.

46

θ'
}
2
11
{
ne
a
l
p
θ' = 71.8351
lipslip plane{112}sθ+90[deg]

slip

pla

ne{

112
}

θ

slip plane{110}θ = 78.9174

zx
Fig. 4.2.11

Snapshots of slip plane in Σ17 twist grain boundaries. (a)slip planemodel, (b) simulation model with yz view in 2nd stress peak.

47

εzz =0.1540
3rd stress peakεzz =0.1539

εzz =0.1540

εzz =0.1541

εzz =0.1542

εzz =0.1543

εzz =0.1575zyFig. 4.2.12

Crack propagation in Σ17 twist grain boundaries after 2nd stress peak.

48

第 5 章 円孔の応力集中部からのき裂発生に関する検討
前章では粒界き裂の発生を検討したが，どのモデルもき裂発生前にまず粒界から結晶内部への塑性変形を生じ，粒内に縞状の不安定原子が分布した．粒界近傍の不安定原子の割合などを調べたが，き裂発生箇所を特徴づけるような有意な結果を見出すことは難しかった．そこで今回は奥行き方向のセル長さを小さくすることで平面的な固有ベクトルの観察を実現するとともに，き裂の発生源として応力集中部（円孔）を導入し，応力集中部からのき裂発生について検討した．

5.1

シミュレーション条件

図 5.1.1 に示すような bcc-Fe の完全結晶薄板セルの中央に，0.3Lx の大きさとなる円孔を設けることで応力集中部を形成した．x, y, z を bcc 構造の [110]，[1̄10]，[001] とする [110] モデルと，[110]，[111]，[112̄] とする [111] モデルの 2 つを対象とした．系の寸法および原子数を表 5.1.1 に示す．全方向周期境界条件下で x, y, z 方向の各応力が0 となるように応力制御しながら 10, 000 [fs] の緩和計算を行った後，x, z 方向のセル辺長は固定する条件下で y 方向に引張ひずみを毎ステップ ∆εyy = 1.0×10−7 増加させて引張シミュレーションを行った．数値積分の時間ステップは 1 [fs] = 1.0 × 10−15 [sec]であるのでひずみ速度に換算すると ε̇yy = 1.0 × 108 [sec−1 ] である．熱揺動による影響を排除するため温度は極低温 (0.1[K]) とし，速度スケーリング法を用いて制御した．

49

Table 5.1.1

Cell size and number of atoms.
［110］

［111］

Lx [nm]Ly [nm]Lz [nm]
Atoms

24.3
32.4
2.0
127,337

24.3
29.8
2.1
122,157

Ly

model

y

Lz

Lx
x

Fig. 5.1.1

Simulation model.

50

5.2

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

5.2.1

応力ひずみ線図および最小固有値の変化
α(1)

図 5.2.1 および図 5.2.2 に応力ひずみ曲線，および，全ての原子の中で最小の ηminα(1)の値をプロットしたものを示す．左目盛りに応力，右目盛りに ηmin の値をとっている．水平破線が固有値 0 の境界を示している．応力はいずれもほぼ線形に上昇した
後，ピークを示し応力急減しているが，詳細に見ると [110] モデルでは矢印で示した

εyy = 0.0375 で小さな応力減少を示しており，その後の応力ひずみ曲線にも小さな振α(1)動が認められる．εyy = 0.0375 の直前に AES の最小固有値 ηmin の値が急減し，矢印で示した点以降で大きく揺らいでいる．εyy = 0.0562 の応力ピーク点では，後述する

15
0Stress, σyy GPa[110] model10σyy
η α(1)

Minimum of eigen value, η a(1), GPa

ようにき裂が発生・進展した． [111] モデルの場合，最小固有値が負になるのは応力

-100
5

-200

0

Fig. 5.2.1

0.02

0.04

strain, εyy

0.06

Stress-strain curve and change in the minimum 1st eigen value η α(1)under [110] tension.

51

0Stress, σyy GPa[111] model
10

Minimum of eigen value, η a(1), GPa

15

-100

σyy
η α(1)

-200
5

-300

-400
0

Fig. 5.2.2

0.02

0.04

strain, εyy

0.06

Stress-strain curve and change in the minimum 1st eigen value η α(1)under [111] tension.

52

ピークの直前であり，その前の応力ひずみ曲線には [110] モデルのような変化は認めにくい．応力ピーク点では転位射出が複数回生じ，それに対応して応力が段階的に減少しており，き裂は発生していない．図 5.2.3 および図 5.2.4 は [110]，[111] の応力
ピーク前後のスナップショットである．AES の第 1 固有値 η α(1) が負の原子を赤く着色している．[110] モデルの図 5.2.3(a) において，εyy = 0.0562 で円孔を中心に斜方に不安定原子が分布している部分は AtomEye の Central Symmetry Parameter(CSP)
でも積層欠陥と判定されており，先の応力ひずみ曲線で変化が現れた時に欠陥がバ
タフライ状に成長したことを確認している．ピーク点以降，円孔赤道部からき裂が発生し図 5.2.3(b) の図のように水平方向へ不安定に進展して破断した．図 5.2.4 の [111]モデルでは εyy = 0.0511 の応力ピーク時に [110] モデルと同じく円孔斜方に η α(1) < 0の原子がわずかに存在している．しかし，[110] と異なり円孔赤道部には η α(1) < 0 の原子は発生していない．εyy = 0.0514 近傍で円孔斜方の不安定原子群の根本から赤道部へ向かって不安定原子（転位）が射出された．(111) 面での角度からこの転位は(112) 面と考えられる．図 (c) では η α(1) < 0 の原子が交差し 1 つに見えるが，CSP により欠陥を可視化すると 2 つの転位芯が確認できる．その後 εyy = 0.0527 で x 方向
周期境界に達して反対側の円孔表面から射出された転位と交差し停滞した．このよ
うな切り欠きからの転位射出や他の転位との交差など複雑な応答を繰り返して塑性
変形し，εyy = 0.0800〜0.0900 の変形後期に図 (e) のような小さなき裂を生じるものの，[110] と異なりき裂の不安定成長は生じずに鈍化と進展を繰り返した．

53

cracking(a) εyy = 0.0562
Fig. 5.2.3(b) εyy = 0.0570

Snapshots before and after the stress-strain peak ([110] tension).

CSP view

emission of dislocationcon
cracking
Fig. 5.2.4

Snapshots before and after the stress-strain peak ([111] tension).

54

5.2.2

固有ベクトルによる検討

η α(1) < 0 の原子について，対応する固有ベクトル {∆εi } の成分 ∆εxx , ∆εyy , ∆εzz ,∆εyz , ∆εzx , ∆εxy からひずみテンソル [∆εij ] を求め，その主軸である固有ベクトルx を求めた．固有値 ε1 < ε2 < ε3 に対応する固有ベクトル x1 , x2 , x3 のうち x3 について，ε3 の大きさに応じてスケーリングし xy 平面に投影したものを図 5.2.5 および
図 5.2.6 に示す．詳細観察のため原子座標の出力間隔を細かくしてとりなおしたので図 5.2.3 および図 5.2.4 で示したスナップショットのひずみとはわずかに異なる．
図 5.2.5(a) は応力ひずみ線図のピーク時で，き裂開口する直前である．き裂が発生する図 (b) では左のグラフ上に (010) 面を破線で示した．き裂開口のする先端では

η α(1) < 0 の原子の変形モードは (010) 面に対して垂直であることが分かる．き裂進
展中の図 (c) でも不安定モードは (010) 面に垂直または平行であり，(110) 面で開口しているが開口した側の結晶方位は (010) 面に沿う方向に傾いている．図 (c) にこの
変形モードの方向を矢印で示した．
図 5.2.6(a) は [111] 引張りで，第 1 応力ピーク後に η α(1) がパルス状に大きな負値を示したときの図である．同図左には {112} すべり面を破線で示した．(a) → (b) にかけて円孔斜方から不安定原子の射出が見られる．射出された不安定原子の x3 には規則性は見られないが，転位射出元である円孔斜方の不安定原子は (112) すべり面に対して，垂直もしくは水平な x3 のどちらかを持っている（青丸印）．図 (c) では
円孔斜方から不安定原子が射出され転位となるが，射出された転位は図 (b) と同じ
く {112} すべり面に対して垂直か水平なベクトルが存在している．転位はすべり面上を動いており，不安定原子は転位によって結晶格子の形状を乱されるため，垂直方向のベクトルはパイエルスポテンシャルに相当する原子面間方向のモードと考え
られる．

55

[Å]

(a) εyy = 0.05620

x3
[Å]

[Å]

(b) εyy = 0.05624

x3

[Å]

(c) εyy = 0.05670
Fig. 5.2.5

Principal axis x3 of unstable deformation mode of [∆εij ] under the [110]tension．

56

(a) εyy = 0.05120

[A]

x3

{112} slip plane[A](b) εyy = 0.05130

[A]

x3

{112} slip plane[A](c) εyy = 0.05140
Fig. 5.2.6

Principal axis x3 of unstable deformation mode of [∆εij ] under the [111]tension．

57

第 6 章 結言
原子弾性剛性係数 Bijα ≡ ∆σiα /∆εj の正値性によって変形・破壊を統一的に議論する研究の一環として，bcc-Fe のき裂を対象とした様々な分子動力学シミュレーショ
ンを行い，Bijα の固有値 η α εi の変化や分布，固有ベクトルなどを議論した．

2 章では，分子動力学法および EAM ポテンシャルについて概説し，弾性剛性係数と格子不安定性について説明した．また EAM における原子応力 σiα や原子弾性剛性係数 Bijα の導出とそれから得られる Bijα について示した．

3 章では，原子間相互作用を遮蔽するシートによるき裂の引張りシミュレーションを行った．[001](110) き裂では応力ピーク後にき裂が水平に進展し破断した．[112̄](111)
き裂では，わずかに進展したものの，その直後に転位射出して鈍化した．き裂面の原子応力と線形破壊力学による応力分布を比較した結果，応力ピーク前に欠陥を生じなかった [112̄](111) き裂は線形破壊力学の式に一致したが，[001](110) き裂では応力ピークより前にき裂から欠陥原子が発生し，破壊力学による予測よりも原子応力の方が高い値をとった．全ての原子の AES の第 1 固有値の η α(1) のうち最小のもの
α(1)(ηmin ) の変化を調べたところ，どちらのき裂でも原子構造の組み替え，すなわち欠
α(1)

α(1)

陥を発生する直前に ηmin の最小値が急減し，-200GPa< ηmin <-150GPa で変形が進行した．

4 章ではき裂の発生を議論するため，き裂発生箇所としてねじり粒界を想定し，
Σ5, Σ13, Σ17 の 3 つの CSL(Consident Site Lattice) ねじり粒界を周期境界によってラミネート状に有するバルクの引張りシミュレーションを行った．いずれも応力ひずみ曲線は 2 つ大きなピーク (Σ17 のみその間に階段状の小さなピークを有する) をα(1)示し，後のピーク以降粒界からき裂を生じて破断した．ηmin の変化および η α(1) < 0

58

α(1)

原子の分布を観察し，(1) 引張りによる粒界近傍の原子構造緩和で ηmin の値が複雑に変化すること，(2) 第 1 応力ピーク前には粒界近傍の η α(1) < 0 の層が増加すること，(3)Σ5 および Σ13 は第 1 ピーク後粒内に (112) 面の双晶変形や (110) 面のすべり変形が生じて η α(1) < 0 の縞が発生したが，Σ17 は第 1 ピークで粒界以外が全て

η α(1) < 0 になった後，やはり応力急減後には粒内に (110) 面と (112) 面に沿った縞状の η α(1) < 0 の原子を残して回復すること，などが明らかになった．破断は粒界とこれらの縞状組織の交差部より生じているように思われるが，上述のように非常に複雑な変形の後，破断したため η α(1) < 0 の原子による評価は難しい．

5 章では円孔による応力集中を設けた薄板モデルの引張りシミュレーションを行った．[110] 方向の引張りでは応力集中部からき裂が発生し，そのまま不安定に破断した．[111] 方向の引張りでは応力ピーク後に転位射出や積層欠陥生成など多くの塑性
変形を伴ってジグザグに時間をかけて成長した．最後に，η α(1) が負となった原子について，対応する固有ベクトル [∆εi ] を成分とするひずみテンソル [∆εij ] の主軸 x3から変形モードを議論した．[110] 引張りでは (010) 面に垂直な x3 が表れており，これは [110] 面で開口する際の，最も端の原子移動に対応する．[111] 引張りでは (112)
面の転位射出を生じる方向に平行または垂直な x3 が現れた．平行なすべりだけでな
く垂直方向にも現れるのは転位芯が切り替わるときの原子面方向のパイエルスバリ
アに相当するものと考える．

59

参考文献

(1) 須賀，船橋，高橋，機論，57A（535），pp.678 - 685，（1991）．
(2) 中谷，北川，中谷，材料，46（6），pp.625 - 630，（1997）．
(3) 松本，北川，中谷，機論 69A（683），pp.1074 - 1081，（2003）．
(4) 北村，梅野，辻，中谷，機論 69A（681），pp.945 - 951，（2003）．
(5) 梅野，嶋田，北村，機論，75A（757），pp. 1247 - 1254，（2009）．
(6) 嶋田，阪本，北村，材料，61（2），pp.151 - 161，（2012）．
(7) 屋代，冨田，機論，67A（656），pp.678 - 684，（2001）．
(8) K. Yashiro and T. Yamane，Proc. of ICCM2014，ISSN 2374-3948 （online），
Scientech Publisher，（2014）．
(9) K. Yashiro，R. Nishikawa and M. Yoshihara，Transaction of the JSME （inJapanese），DOI:10.1299/transjsme.15-00271，（2015）．
(10) K. Yashiro，Com. Mater. Sci.，131，pp.220 - 229，（2017）．
(11) K. Yashiro，Com. Mater. Sci.，147，pp.72 - 78，（2017）．
(12) M. W. Finnis and J. E. Sinclair，Phil. Mag. A，50（1），pp.45 - 55，（1984）．
(13) M. Mullins，Scr. Metall. ，16（6），pp.663 - 666，（1982）．
(14) 上田，コンピュータシミュレーション，朝倉書店，（1990）．
60

(15) 洲之内，サイエンスライブラリ–理工系の数学=15，数値計算，サイエンス社，（1978）．

(16) 香山，固体材料の電子状態の基礎，材料学会勉強会資料，（1993）．
(17) M. W. Finnis，and J. E. Sinclair，Phil. Mag.A，50(1)，pp.45 - 55，（1984）．
(18) M. W. Finnis and J. E. Sinclair，Phil. Mag.A，53（1），pp.161，（1986）．
(19) S. Nosé，J. Chem. Phys. ，81（1），pp.511 - 519，（1984）．
(20) Hoover, William G., Phys. Rev. A, 31(3), pp.1695 - 1697
(21) D. C. Wallace，Thermodynamics of Crystals，Wiley，Newyork，（1972）．
(22) J. Wang，et al. ，Phys. Rev. Lett. ，71（25），pp.4182 - 4185，（1993）．
(23) J. Wang，et al. ，Phys. Rev. B，52（17），pp.12627 - 12635，（1995）．
(24) 泉，増田，機械・材料設計に生かす実践分子動力学シミュレーション，森北出
版，（2013）．

(25) P. Sung, T. Chen，Com. Mater. Sci.，102，pp.151 - 158，（2015）．
(26) Q. Feng，et al. ， M aterials&Design，120，pp.193 - 203，（2017）．
(27) Yuan Qi, et al, RSCAdvances, 81, pp.65651 - 66501, (2015)(28) C.B.Cui, et al, Materials Science and Engineering:A, 609, pp.102 - 109, (2014)(29) 都留，青柳，加治，下川，日本金属学会誌，78，第 1 号， pp.45 - 51, (2014)(30) 久保田，松本，宮崎，日本機械学会論文集 (A 編),73，773, pp.989 - 996，(2007)

61

謝辞

本研究を遂行するにあたり，屋代如月教授には浅学非才な著者に対し懇切丁寧に指導していただきました．ここに心よりお礼申し上げます．本論文を完成させるにあたり，広い視野から研究全般に対して多くのご助言をいただきました内藤圭史助
教にも心より感謝いたします．ともに切磋琢磨し合った同輩の寺田稜君や修士 1 年の坪井伶以君には研究に関して多大なアドバイスを頂きました．また他の研究室メンバーや機械システム工学専攻の同期たち，2 年前に卒業した機械システム工学科の方々にも御礼申し上げます．色々な方々の支えがあって今の自分があります．最後に，6 年間の学生生活を暖かく見守り精神的にも経済的にも支えて頂いた母や親族に心より感謝いたします．ありがとうございました．

62