修士論文

原子弾性剛性係数に基づく局所不安定性解析：fcc ニッケル中のき裂先端の力学特性評価

指導教員:田中克志

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

147T327T
片山寛

Master thesis

Local lattice instability analysis
based on atomic elastic stiﬀness:study on crack tip in fcc nickel.

Hiroshi KATAYAMA

February 2016
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan

要約
本論文では，モード I 型き裂を有する fcc-Ni を対象とし，分子動力学法を用いて引張シαα
ミュレーションを行い，き裂進展挙動と原子弾性剛性係数 Bijの行列式 detBij，ならびにその固有値 η α の関係について検討を行った．まず [001](010)，[1̄10](001)，[1̄1̄2](111)の三種類の結晶方位で平板状の周期セル中央に貫通き裂を有する系に対し引張シミュレーα
ションを行った．いずれのモデルも無負荷平衡状態では detBij< 0 の原子は存在してα
おらず，き裂近傍に detBij< 0 となる原子が生じた後，き裂先端から欠陥が発生し応

力-ひずみ線図に非線形性が現れた．また固有値の変化を調べた結果，構造変化が起こる
点 (欠陥が発生する点) では η α(1) の最小値がほぼ一定の負の値 ([001](010)，[1̄1̄2](111)は-40GPa，[1̄10](001) は-20GPa）となっていた．次に周期セルの中央と四隅に貫通き裂αを配置したダブルき裂モデルでの解析を行った．欠陥生成時の detBij< 0 原子の数はいずれも単一き裂モデルのほぼ二倍であり，き裂先端の不安定領域は一定の寸法である可能性が示唆された．次に周期セルの厚さ方向を 1.5 倍とし，厚さ方向自由度を変えたモデル
で解析を行った．厚さ方向の自由度を増やしたことによって変形モードがわずかに変化しα
たが，detBij< 0 の原子が初めて現れるひずみにはどのモデルも変化がなく，その原子α数も 1.5 倍 ([001](010) のみ 1.4 倍) となっていることから，き裂先端に detBij< 0 の原子が発生するときの力学条件は厚さ方向の自由度の影響はない．全ての解析を総合して，
α< 0 が出現するときの局所力学状
結晶方位による差はあるが，き裂先端に初めて detBijα
態は常に同じ (detBij< 0 の原子列＝不安定領域が一定) であり，またその後き裂先端かαら転位射出しき裂が鈍化・開口するときの Bijの第一固有値が一定の値を示すという結果を得た．

Summary
In order to discuss the local stability at a crack tip, we performed molecular dynamics (MD) simulations on mode I through crack in fcc nickel and evaluate theαdeterminant of atomic elastic stiﬀness Bij, and its eigenvalue η α . First, we per-formed tensile simulations on single crack models with three crystal orientation of theα[001](010)，[1̄10](001)，and [1̄1̄2](111), respectively. There is no detBij< 0 atoms
αin the no load equilibrium, and detBij< 0 atoms appear at the crack tip at thestrain 0.035，0.042，and 0.026.Then inelastic deformation takes place at the tip andthe stress-strain curves show nonlinearity at the point. The minimum 1st η α(1) ofαBijbecomes negative and shows almost constant value,-40GPa for [001](010), and-20GPa for [1̄10](001), during crack propagation. Next, we simulate with the double
crack model of center and corner crack in the periodic cell in order to discuss the eﬀect
αof diﬀerent crack array. The number of detBij< 0 atoms is simply doubled at thedefect nucleation, so that we conclude that the size of the unstable region is almost
constant, despite of the diﬀerent force field. Finally, we analyses thicker model by 1.5
times to change the flexibility of the thickness direction. Deformation mode is slightly
changed only in the [1̄1̄2](111) crack; however, the strain at the first emerging point ofαα
detBij< 0 atoms is identical to the precious ones. The number of detBij< 0 atoms

is just 1.5 times (the [001](010) model 1.4 times) larger than that of the base models.
So the flexibility of thickness direction doesn’t aﬀect the mechianical condition for theαemergence of detBij< 0 atoms at the crack tip. As a summary, we observed the local
αcriteria for the emergence of detBij< 0 atoms at the crack tip, the 1st eigenvalue ofαdetBijfor the subsequent defect nucleation which leads the crack propagation.

目次
1

緒論

1

2

解析手法の基礎42.1

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

42.2

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

52.3

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

7

2.4

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

82.5

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

9

モード I 型き裂の引張シミュレーション

15

3.1

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

15

3.2

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

17

ダブルき裂モデルでのシミュレーション

29

4.1

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

29

4.2

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

30

厚さ方向自由度の影響

38

5.1

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

38

5.2

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

38

結論

46

3

4

5

6
参考文献

48

謝辞

50

i

1 緒論

分子動力学シミュレーションは材料の特性を個々の原子の運動にまで遡って解析する手法であり，材料科学だけでなく，高分子化学や生体化学など様々な分野で活用されている．計算機が飛躍的に性能の向上，およびナノレベルでの材料の変形・破壊実験が可能となった今日，ナノ〜サブミクロン領域の現象を疑似的に観察する手法として分子動力学シ
ミュレーションは広く用いられるようになってきた．しかしながら，分子動力学シミュ
レーションの中で見出した原子レベルでの変形メカニズムについて，どのような指標を用いて普遍的・定量的な議論に結び付けるかということは依然として問題である．実在の材料の挙動をできるだけ精密に再現するポテンシャル関数の開発に膨大な努力がなされているが，「似ているがそのものではない」モデル材料の挙動であるため，局所エネルギーや応力などのしきい値を求めても意味がない．そのため応力拡大係数などの連続体理論を指標とした議論に頼りがちだが，き裂先端などの原子構造に起因したわずかな違いなどは議論できず，モデル材料における連続体理論の適用限界を示すのみとなる．
著者らのグループでは，局所変形抵抗を表す物理量として原子弾性剛性係数 (Atomic

Elastic Stiﬀness, AES) に着目し，その正値性を基準とすることで，用いたモデル材料やポテンシャルエネルギーによらない普遍的な支配力学法則の解明を目的とした議論を行ってきた 1) ．AES による検討は，fcc の Ni や bcc 構造の Fe，アモルファス金属やダイヤモンド構造を持つ Si などに広く展開してきた．1)–4) 本研究で対象とする fcc 構造の Ni については，単結晶ナノワイヤの引張シミュレーションより，表面から部分転位が発生・移動する際には AES が負の原子が転位前縁に現れることを明らかにしている 1) ．アモルファ
ス構造への適用では，AES の正値性によって局所クラスタの崩壊及び再構成が理解でき1

ることを報告している 2) ．Si の検討では，表面や界面による内部不均一性に着目した検討を行い，これらの欠陥部分の AES は無負荷平衡状態で負または 0 に近い値を示し，欠陥から離れたバルクの部分と大きな差があるものの，外力を与えるとバルク部分の AES がそれに近づく「系の均一化」を生じていることを報告している 3) ．また，bcc 鉄中のらせん転位に関する検討では，転位芯原子の AES の行列式の値は無負荷平衡状態では負とならず，外力を受けて転位が移動する際に初めて負となることなどを明らかにしている 4) ．
分子動力学シミュレーションによるき裂の解析は古くから行われてきた．計算機能力が十分でなかった 1970 年代に既に，き裂の開口時の原子結合切断を簡略化してモデル化することで，格子力学的な視点からき裂の進展を議論した検討がなされている 5) ．Gumbsch
らはき裂の進展を局所格子の破断と捉えたモデルを提案し，第一原理計算を用いた精密な評価を行っている 6) ．また，Shimada らは系の全自由度に対してヘッシアンを計算するという厳密な手法を用いて，Si き裂についてき裂進展の解析を行い，き裂先端の一格子程度の領域のみが不安定となっていることを明らかにした 7) ．しかし系のヘッシアンを直接扱う方法では計算量が膨大となるため，100 万原子程度の大きな系では直接解くことがで
きず，き裂近傍の一部をとりだして弾性体中の境界条件を新たに与えて解いている．一方
著者らのグループで提案している AES 解析は、後で述べるように 6 × 6 のマトリックス
で表されるので，全原子について固有値解析を行うことが容易である．そのため大規模な系での解析を行うことが可能であり，また熱搖動の影響やき裂の鈍化についても考慮することができる．き裂を対象とした AES 解析はこれまでに bcc-鉄, シリコンについて進めてきた 8), 9) ．bcc 鉄での解析では，き裂が進展する際に塑性変形がほとんど発生しない結晶方位 ([001](110),[1̄1̄2](111)) では，き裂の大きさや数，配置の仕方に関わらずき裂先端に現れる負の AES 原子数が一定であることを明らかにしている 9) ．また Si の解析では，
モード I 負荷条件下でき裂先端に非常に大きな値の正の AES 原子が現れており，これが2

転位の射出を妨げ Si のぜい性的な挙動をもたらしている可能性や，き裂進展直前に大きな負の第一固有値を生じる前駆挙動があること，またその固有ベクトルは [001](010) き裂
はモード I，[001](110) き裂はモード II，[112](111) き裂ではモード III の変形モードであることなどを明らかにしている 8) ．
本研究では，fcc 金属の Ni をモデル材料として同様のモード I 型き裂の引張シミュレー
ションを行い，き裂進展/鈍化時の AES を明らかにすることを目的とする．以下に各章の概略を示す．第 2 章では MD 法の基礎理論を述べ，用いた原子間ポテンシャルにおける
原子弾性剛性係数の定式化を行う．第 3 章ではモード I 型き裂の引張シミュレーションを行い，その変形挙動を AES 解析を用いて検討を行う．第 4 章ではき裂を複数導入することでき裂配置の影響について検討を行う．第 5 章ではシミュレーションセルの厚さ方向の原子数を変えたモデルを用いて解析を行い，厚さ方向の自由度の影響について検討する．
第 6 章では得られた結果を総括する．

3

2 解析手法の基礎

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

(2.1)

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

Fα = −

∂Etot
∂r α

(2.2)

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

(2.6)
(2.7)

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

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

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

2 F

α

(t)

(2.8)
(2.9)

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

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

2 F

α

(t)

mα

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

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

3
1 α α αm vi vi = kB T
2
2
5

(2.12)

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

(

3kB T0mα

viα0 =)0.5(2.13)

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

(
viα =

3kB T
mα

)0.5
(2.14)

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

viα0
=
viα

(

T0
T

)0.5
(2.15)

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

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

Fα
i (t)
(2.16)mα

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

6

2.3 原子埋め込み法ポテンシャル
分子シミュレーションにおいて，系の振る舞いを精度良く再現するためには分子に作用するポテンシャルエネルギーが重要となる．量子力学に基づき，電子や原子核のハミルトニアンから系のポテンシャルエネルギーを精密に求める第一原理計算も試みられているが，計算量が極めて膨大となるため扱える粒子数が限られてしまう．そこで多数の原子の挙動を求める場合には原子間相互作用を簡略評価する原子間ポテンシャルが通常用いられる．これまでに数々のポテンシャル関数が提案され用いられてきたが，解析の目的に対して適切なポテンシャルの設定を誤ると，いかに計算を重ねたところで無意味な検討となってしまうため，ポテンシャル関数の選択は慎重に行わなければならない．
原子埋め込み法 (Embedded atom method)12)13) は Daw と Baskes によって Stot とZaremba の擬原子法 14) の考えを用いて提案されたポテンシャル関数で，金属中の多体効果を精度良く再現することから広く用いられている．EAM では密度汎関数理論に基づき，まず金属材料における系全体のポテンシャルエネルギー Etot は原子を価電子雲中に埋め込むエネルギーと原子間の２体間相互作用の和で与えられるとする．さらに，埋め込みエネルギーは埋め込む位置の電子密度にのみ依存すると仮定することによって，系全体のエネルギーは次式のように表される．

Etot =∑1 ∑
Vij (rij ) +Fi (ρi )2
i(2.17)

i,j(̸=i)

ここで，ρi は，原子 i の位置における電子密度，Fi(ρi ) は電子密度 ρi の位置に原子を埋め込むエネルギー，Vij (rij ) は距離 rij 離れた原子 i と j のクーロン相互作用である．

7

電子密度 ρα は周囲の原子 β からの寄与 ϕ(rαβ ) の重ね合わせで与えられると仮定しαρ =neighbor∑ϕ(rαβ )

(2.18)

β(̸=α)

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

2.4 高速化手法
原子論に立脚したシミュレーションでは原子間相互作用の計算に大部分の負荷が集中する．しかし実際の結晶中では近接原子の遮蔽効果により，離れた原子からの相互作用はほとんど及ばない．そのため，相互作用打ち切り (カットオフ) 半径 rc を導入し，その半径内の原子からの寄与のみを考慮すればよい．例えば 2 体ポテンシャルを用いる場合，系の全粒子数を N とすると系の内部だけでも組み合わせは N (N − 1) 通り存在し，粒子数が増加することで飛躍的に計算コストが増えることがわかる．この計算コストを減らす方法として領域分割法を用いた．領域分割法では系を間隔が rc 程度である格子状に分割し，対象とする原子の存在する領域 (図 2.1 着色部) および隣接する領域 (図 2.1 斜線部) の原子の内，距離が rc 以下であるものを探索する.

8

y

by
0

x

bx

Fig.2.1 Schematic of domain decomposition.

2.5 原子弾性剛性係数
2.5.1

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

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

1
σij =
V

(

∂F
∂ηij

)

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

(

)
∂σij
Bijkl ≡
∂ε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 の最小固有値が負であれば，対応するひずみに対して応力が負になる自発的変形モードが存在することになり，不安定と判別される．

10

2.5.2原子系への拡張αα
局所の安定性を評価するための原子弾性剛性係数 Bijklの算出に必要な原子応力 σij
なα
らびに原子弾性係数 Cijkl
は，各原子周りの微小ひずみに対するポテンシャルエネルギーの 1 次，2 次変化量として導出される．

2.5.3

原子応力

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

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

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

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

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

ηij =

]
1[
Gij − δij
2(2.24)

1dGij
2(2.25)であり，その微小量dηij =と式 (2.23) の関係から次の関係が得られる．

r̄iαβ r̄jαβ
∂rαβ
=
∂ηij
rαβ

(2.26)

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

( ∑
)
N ∑
N
∂rαβ ∂Etot
12 α∂ηij ∂rαβ

1 ∂Etot1σij ==V ∂ηij
V
1
=
V

N
∑

β(̸=α)

N (
∑

′

α

′

F (ρ )ρ (rαβ

α β(̸=α)) αβ αβ1 ′ αβ ri rj) + ϕ (r )2rαβ

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

α
σij
=

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

(2.27)

β(̸=α)と定義すると，系の応力は

1 ∑ ασ
N α ij
N

σij =

となる．

12

(2.28)

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

1 ∂ 2 Etot
V ∂ηij ∂ηkl

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

(∑
αγ αγ )
NN
N ∑N
1 ∑ ∑ ∂rαβ ∂1 ′ αγ ri rj′ α ′ αγCijkl ={F (ρ )ρ (r ) + ϕ (r )} αγ2V α∂ηkl ∂rαβ α2
rβ(̸=α)

=

1
V

[∑
N ∑
N
α β(̸=α)γ(̸=α)() αβ αβ αβ αβρ′ (rαβ ) ri rj rk rl
′ α
′′ αβ
F (ρ ) ρ (r ) −rαβ(rαβ )2

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

α(2.30)

β(̸=α)

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

α
Cijkl
=

1
V /N

[ ∑
N
β(̸=α)() αβ αβ αβ αβρ′ (rαβ ) ri rj rk rl′′ αβ
F (ρ ) ρ (r ) −rαβ(rαβ )2
′

α

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

β(̸=α)

これより，系の弾性係数は

1 ∑ α
Cijkl =C
N α ijkl
N13(2.32)

のように原子弾性係数の平均となる．

2.5.5

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

14

3 モード I 型き裂の引張シミュレーション

3.1 シミュレーション条件
図 3.1 に示すように，平板状の周期セル中央に貫通き裂を有する系を解析対象とする．
き裂幅 2c はセル寸法 Lx の 1/5 とした．図の 3 つの結晶方位を対象とし，[001](010)，[1̄10](001)，[1̄1̄2](111) のき裂の解析を行った．き裂変位 h は Voter の EAM ポテンシャ
ル 15) での Ni のカットオフ半径 0.479nm より大きく，かつき裂先端半径 ρ がグリフィス
き裂の条件 8a/π（a は Ni の格子長さ 0.352nm）を満たすよう，ρ=0.4nm，h=2ρ=0.8nmとする．き裂部の原子を取り除いた各モデルの原子総数はそれぞれ 178200, 179074,

173864 である．
各モデルについて，全方向周期境界の下で垂直応力が 0 となるようセル辺長を調整しな
がら 100000fs の緩和計算を行った後，y 方向にひずみ増分 ∆εyy = 1.0 × 10−6 [/fs] を毎
ステップ与えることで引張シミュレーションを行う．数値積分の時間ステップは 1fs であるので，ひずみ速度に換算すると 1.0 × 109 [/s] である．温度は熱搖動の影響を排除する
ため 0.1K とし，速度スケーリングにより制御している．応力-ひずみ応答，ならびに，高エネルギーの原子位置から決定したき裂長さ変化を求めると共に，1000fs 毎 (き裂進展近αならびに原傍では 100fs 毎) に記録した原子座標データから，各原子位置の原子応力 σijαααα
子弾性係数 Cijから Bijを評価し，detBij
および固有値 η α(k) を求めた．なお detBijの値は 0[K] 完全結晶における値 8.8796 × 1012 [GPa] で無化元化している．

15

[001](010) crack

[010]

[001]

[-1-12](111) crack[-110](001) crack

[100]

[001]

[-110][110]

Fig.3.1 Simulation models for mode I through crack.

16

[111]
[-1-12][-110]

3.2 シミュレーション結果および考察
応力-ひずみ線図と負の AES 原子数変化

3.2.1

図 3.2 に応力ひずみ線図及び負の AES 原子数の変化を示す．いずれのモデルにおいて
αも，き裂が存在するにも関わらず無負荷平衡状態では detBij< 0 の原子は存在しない．
α
ひずみ εyy = 0.03 近傍から detBij
が負となる原子が現れ始める．初めて負の原子が現れた点を応力-ひずみ曲線上に×印で示した．また図中の○印は，可視化ソフト atomeye18)の central symmetry parameter によりき裂先端からの欠陥生成（表面原子以外に生じた欠陥）が確認された点を示している．応力-ひずみ線図を見ると，全てのモデルにおいて○印の後から非線形性が現れており，非弾性変形の開始に対応していることがわかる．
[001](010)model[-110](001)model[-1-12](111)model[001](010)model[-110](001)model[-1-12](111)model

Number of negative detBαij atoms

Stress, σyy, [GPa]

15

10

5

0

0

0.1

1000

0

0

0.1

strain, εyystrain, εyy(i) stress-strain curves.

2000

α < 0 atoms.
(ii) Change in the number of detBijαFig.3.2 Stress-strain curves and change in the number of negative detBij
atoms.

17図 3.3 にシミュレーション中の各モデルの変形の様子を atomeye18) の Mises ひずみで色分けして示す．Mises ひずみは 0〜0.2 の範囲で色分けしており，ひずみの値が大きい場
所ほど明るい色で表示されている．[001](010) モデルでは，応力-ひずみ線図で×で示した
点 (図 (a)(ii)) をみると，き裂の先端に高い Mises ひずみが確認される．その後，き裂先端からの転位射出と共にき裂はぜい性的に進展した (図 (a)(iii)，(iv))．ひずみ εyy = 0.057
でき裂進展開始し，その後一定応力で破壊が進行する．[1̄10](001) モデル，[1̄1̄2](111) モ
デルでも同様に非線形性が現れる点 (図 (ii)) ではき裂先端に高い Mises ひずみが観察されるが，その後は転位射出によってき裂が鈍化する (図 (iii)).[1̄10](001) モデルではき裂
から転位射出することで応力が鈍化し，εyy = 0.059 でき裂進展開始する．[1̄1̄2](111) モ
デルではき裂の左側からは多数の転位が発生し，右側はき裂上端のみから転位が射出された．

18

(i) εyy=0(ii) εyy=0.048(i) εyy=0(i) εyy=0(iii) εyy=0.057(iv) εyy=0.072(ii) εyy=0.041(iii) εyy=0.059(b)[110](001)crack(iv) εyy=0.086(ii) εyy=0.030(iii) εyy=0.051-(c)[112](111)crack(iv) εyy=0.072(a)[001](010)crack

Fig.3.3 Snapshots of atoms colored with Mises strain in AtomEye．

19

3.2.2

き裂先端の負の AES 原子の分布α
最初に detBij
が負の原子が観察されたステップにおけるき裂近傍のスナップショットαを図 3.4 に示す．図では detBij< 0 の原子を赤で，それ以外を緑で着色して示している．
α
いずれのモデルでもき裂先端に detBij
が負の原子が観察されており，奥行き方向 (周期方α< 0 となっている．[001](010)，[1̄1̄2](111) モデルでは表面原子が向) の原子列が detBij

負となっているのに対して，[1̄10](001) モデルでは表面から一原子内側の原子が負となっている．また図 3.5 に，Atomeye の Central Symmetry Parameter(CSP) によって，き裂表面以外は欠陥と判定される原子が存在しないステップでの負の AES 原子の分布 (左)と，2000〜4000fs 後で欠陥と判定された原子のスナップショット (右) を示す．ぜい性的なき裂進展を示した [001](010) モデルでは，き裂先端に負の AES 原子が多く密集している．一方，転位射出によりき裂が鈍化した [1̄10](001),[1̄1̄2](111) き裂では負の AES 原子を起点として転位が射出されている．

(i)[001](010) model(ii)[110](001) model-(iii)[112](111) model

Fig.3.4 Snapshots of atoms at emerging point of negative AES atoms.

図 3.6 は応力-ひずみに非線形性が現れてしばらく後のスナップショットを，AES と

CSP で表したものである．先の図 3.5 と違い，同一ステップでの原子配置を AES と CSPで色分けしている点に注意されたい．[001](010) き裂では，き裂先端の原子配置が乱れ，

20

その前縁に負の AES 原子が多くみられる．[1̄10](001) き裂では転位射出後に積層欠陥は消滅するので，先の図 3.5 のすべり面に沿った欠陥はなくなり，き裂底に生じたステップ
部分にのみ CSP で欠陥が認められる．これに対応して，負の AES 原子は図 3.5 のようにき裂の底部だけではなくステップ中央にも生じている．[1̄1̄2](111) き裂では，き裂右側から射出された完全転位 (右上) が不安定となっているが，左側ではき裂先端から発生した
部分転位のループ (積層欠陥部) が途中で止まっている．この転位の前縁には負の AES 原子がみられないため成長が止まった可能性がある．

21

(i)εyy=0.048(i)εyy=0.043(i)εyy=0.029(a) [001](010) crack(b) [110](010) crack--(c) [112](111) crack(ii)εyy=0.052(ii)εyy=0.045(ii)εyy=0.031

Fig.3.5 Snapshots of negative AES atoms(left) and defect atoms visualized bycentral cymmetry parameter(right).

22

(a) ε=0.058([001](010) crack)(b) ε=0.06([110](001) crack)-(c) ε=0.42([112](111)crack)Fig.3.6 Snapshots of negative AES atoms(left) and defect atoms visualized bycentral cymmetry parameter(right).

23

3.2.3固有値の変化αBijεj = η α εi の固有方程式の解 η α(k)(k = 1〜6) を各原子について求め，その変化を調べた．図 3.7 に引張シミュレーション中の第一固有値が負の原子数の変化を示す．初α
めて第一固有値が負の原子が現れる点は detBij< 0 の原子が現れる点と完全に一致して
α< 0 の原子数変化のグラフとほぼ同じ概形となっている．また図 3.2(ii) に示した detBij

いることから，不安定原子の大部分は第一固有値のみが負の原子を示していることがわかる．図 3.8 は，各ひずみで η α(1) が最も小さい原子を調べ，その値をプロットしたものであαる．図中にはき裂長さの変化を青線で示してある．また図中の×印は初めて detBij<0の原子が現れる点である．いずれも×印の前から減少しており，[001](010)，[1̄1̄2](111)αεj = η α εi = σi より) で一定の値をとった後，き裂の長き裂ではほぼ-40GPa(単位は Bij

さが変化している．[1̄10](001) き裂では一度正値に回復するような複雑な挙動を示してい

Number of negative ηα(1) atoms[001](010)model[-110](001)model[-1-12](111)model

2000

1000

0
0

0.1

strain, εyy

Fig.3.7 Change in the number of negative η α(1) atoms.

24

るが，転位射出によってき裂近傍の不安定性が回復した可能性がある．き裂長さが増加す
る際は η α(1) の最小値は約半分の-20GPa 程度となっている．表 3.1 に×印，非弾性変形
開始点 (CSP により欠陥原子を確認した点) におけるひずみの値と，そのときの負の η α(1)を有する原子数，最小の η α(1) の値をまとめて示した．×印における負の原子数，η α(1)の値に大きな差は無く，局所不安定が発生するときの力学条件はモデルによる違いがないと考えられる．一方，その後欠陥となって実際に構造変化を生じる点では負の原子数 (領域) には大きな違いがあるが，η α(1) の最小値 (変形モードに対応) は-15〜-35GPa の範囲
内にある．

Minimum of eigenvalue
Crack length

20
0

0.4

-20
0.2

-40
-60
0

0.05

20
0

0.4

-20
0.2

-40
-60
0

0.1

0.05

strain, εyystrain, εyy(a)[001](010) crack(b)[110](001) crack

0.1-

Minimum of eigenvalue
Crack length

0.6

20
0

0.4

-20
0.2

-40

Crack length

Minimum of eigenvalue, ηα(1)

0.6

-60
0

0.05

strain, εyy

0.1

--

(c)[112](111) crack

Fig.3.8 Change in the minimum 1st eigenvalue and crack length.
25

Crack length

Minimum of eigenvalue, ηα(1)0.6Crack length

Minimum of eigenvalue, ηα(1)Minimum of eigenvalue
Crack lengthTable 3.1 Strain, number of η α(1) < 0 atoms, and minimum of η α(1) at emerging pointof negative AES and defect.
emergence point of η α(1) < 0

defect nucleation by CSP
εyy

atoms

min. of η α(1)

-5.08

0.050

294

-35.91

56

-2.99

0.043

168

-14.63

48

-3.31

0.030

69

-22.44εyyatoms[001](010) crack model

0.035

40

[1̄10](001) crack model

0.042[1̄1̄2](111) crack model

0.026

min. of ηα(1)

図 3.9 は引張シミュレーション中に第二固有値が負となった原子数変化を示したものである．図中にはき裂長さの変化を青線で示してある．[001](010) き裂では，第二固有値が負の原子数がパルス状に増加した直後にき裂が進展した．[1̄10](001) モデルでは非弾性変形の開始点 (εyy = 0.043) では第二固有値が負の原子は観察されず，εyy = 0.07 近傍で初めて第二固有値が負の原子が現れる．青線のき裂長さはこの点から急激に増加している．[1̄1̄2](111) き裂では他に比べ早くから第二固有値が負の原子がみられるがその数は少ない．転位射出による鈍化のため，き裂長さ変化は他に比べ緩やかであるが第二固有値が負の原子が出現した点からき裂の成長が始まっており，き裂の進展は第二固有値が負の原子によって測定できることが示唆される．

26

40

Number of negative ηα(2) atoms

20
0.2
10

0

0.05

20
0.2
10
0

0.1

0.4

30

0

0.05

(ii)[110](001)model(i)[001](010)model

0.6

50

0.4

30
20

crack length

Number of negative ηα(2) atoms

40

0.2
10
0

0

0.05

0.1

-(iii)[112](111)model

Fig.3.9 Number of negative η α(2) atoms and change in crack length．

27

0.1

crack length

0.4

30

crack length

Number of negative ηα(2) atoms

40

0

0.6

50

0.6

50

図 3.10 に第二固有値が負の原子が出現したステップでのき裂近傍のスナップショッ
トを，η α(1) < 0 の原子を赤，η α(2) < 0 の原子を青，それ以外を緑で示す．第二固有値αが負となった原子は第一固有値も負となっている．したがって，detBijの判定では正と判定される原子であることに注意が必要である．[001](010) き裂ではき裂先端の表面側に，第二固有値が負の原子が二列 (奥行き方向全ての原子)，き裂面上下に発生している．

[1̄10](001) き裂では転位射出により生じた図中右下の入り込み部に一列，また第一固有値が負の原子列が凝集した右上結晶内部に一列，第二固有値が負の原子列が認められる．

[1̄1̄2](111) モデルでは同様に転位射出による右上の入り込み部に第二固有値が負の原子列があり，この部分はわずかにき裂進展した．左下にも第二固有値が負の原子がみられる
が，この部分では奥行き方向に連続しておらずき裂進展は生じていない．

(i)ε=0.048(ii)ε=0.071(iii)ε=0.061

Fig.3.10 Snapshots of atoms at emerging point of η α(2) < 0 atoms.

28

4 ダブルき裂モデルでのシミュレーション

4.1 シミュレーション条件
図 4.1 に示すように，前章で用いたセルの四隅に 4 分の 1 貫通き裂を配置したモデル
で解析を行った．き裂幅 2c は前章と同様にセル寸法 Lx の 1/5 とした．き裂部の原子を取り除いた各モデルの原子総数はそれぞれ 178200, 177828, 173008 である．前章と同様に，全方向周期境界の元で垂直応力が 0 となるようセル辺長を調整しながら温度 0.1K で100000fs の緩和計算を行った後，y 方向にひずみ増分 ∆εyy = 1.0 × 10−6 [/fs] を毎ステッ
プ与えることで引張シミュレーションを行う．その後の引張条件についても前章と同様とした．

[001](010)crack

[010]

[001]

[-1-12](111)crack[-110](001)crack

[100]

[001]

[-110]

[110]

[111]
[-1-12]

Fig.4.1 Simulation models for mode I through crack (2 crack model).

29[-110]

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

応力-ひずみ線図と負の AES 原子数変化

図 4.2 に応力ひずみ線図及び負の AES 原子数の変化を示す．初めて負の原子が現れた
点，ならびに CSP で欠陥が認められた点を応力-ひずみ上にそれぞれ×，○印で示した．

[110](001) き裂，[112](111) き裂は前章と同様，○印以降非線形性が現れた後応力上昇が鈍化し，弾完全塑性体のような挙動を示している．一方，[001](010) き裂では応力-ひずみ曲線はピーク後著しく低下した．図 4.3 に変形の様子を，負の AES 原子を赤，それ以外を緑で着色して示す．図 (i) は各モデルにおける応力-ひずみ線図の×印近傍でのスナップ
ショットを示している．その後 [001](010) き裂では，き裂先端に AES が負の原子が集中[001](010)model[-110](001)model[-1-12](111)model

Number of negative detBαij atoms

Stress, σyy, [GPa][001](010)model[-110](001)model[-1-12](111)model

10

5

0

0

0.1

strain, εyy(i) stress-strain curves.

2000

1000

0

0

0.1

strain, εyy

α < 0 atoms.
(ii) Change in the number of detBijαFig.4.2 Stress-strain curves and change in the number of negative detBij
atoms.

30

して生じた後に，四隅のき裂がぜい性的に成長している (図 (iii),(iv))．変形後期では四隅のき裂が合体して材料がほぼ分離した状態となったため先の応力急減を生じたものと結論づけられる．[1̄10](001)，[1̄1̄2](111) き裂ではき裂先端に負の AES 原子が現れるが，前章と同様に転位射出によるき裂鈍化を示した．(図 (iv))．

(i)ε=0.037(ii)ε=0.051(iii)ε=0.070(iv)ε=0.113(iii)ε=0.077(iv)ε=0.097(iii)ε=0.054(iv)ε=0.086(a)[001](010)crack(i)ε=0.045(ii)ε=0.067(i)ε=0.027(ii)ε=0.040(b)[110](001)crack-(c)[112](111)crackαFig.4.3 Snapshots of atoms colored by positive/negative detBij.

31

4.2.2

負の AES 原子が発生するひずみと分布形態

負の AES 原子が発生したステップでのひずみの値を表 4.4 にまとめた．いずれのモデ
ルもわずかに単一き裂の場合のほうがひずみが小さいが，ほぼ同じ値となっている．複数
き裂の方が応力遮蔽によって局所不安定に達するのが単一き裂より遅くなったものと考えられる．負の AES 原子発生時の中央のき裂と四隅のき裂周囲での負の原子分布を，前章で用いた単一き裂モデルと合わせて図 4.4 に示す．図中では負の原子を赤，それ以外を緑で示している．[001](010) き裂では，単一き裂時にはき裂の両端に奥行き方向 (周期方向)に整列して負の AES 原子が発生したが，複数き裂の場合では中央のき裂の右下と，四隅のき裂の左下にのみ観察され，また奥行き方向にも全ての原子列が負になってはいなかった．ただし，この直後のステップ (1000fs 後) のスナップショットでは複数き裂でも単一き裂と同様の分布となった．[1̄10](001) き裂では，単一き裂では一原子内部の原子のみが負となったのに対して，複数き裂ではどちらのき裂も表面原子と，その一つ内部の原子が負の AES を示し，単一き裂時に負の値を示した原子列と異なる原子列が負となっている．

[1̄1̄2](111) き裂では，単一き裂ではき裂の先端原子とその上下の原子列が負となったが複数き裂では先端原子のみが負となった．ただし，[001](010) き裂の場合と同様，1000fs 後の原子配置では負の原子分布は単一き裂と同じになった．

Table 4.1 Strain at emerging point of negative AES atoms.
single crack model

double crack model[001](010)

0.035

0.037

[1̄10](001)

0.042

0.044

[1̄1̄2](111)

0.026

0.027

32

(a)single crack(b)double crack (left:center right:corner)(i)[001](010)crack-(ii)[110](001)crack--(iii)[112](111)crack

Fig.4.4 Snapshots of atoms at emerging point of negative AES atoms.

33

表 4.2 に各モデルでの非弾性変形開始直前 (CSP での欠陥生成点，非線形性が現れる
点) のひずみと負の AES 原子数を，前章で用いた単一き裂のものと合わせて示す．ただし[1̄10](001) き裂ではき裂の右端と左端で欠陥原子が出るタイミングが異なるため，それぞれ欠陥原子が出る直前のステップでの原子数を求め，足したものを用いている．鉄による先行研究 9) では，ほとんど塑性変形を生じずへき開する結晶方位 ([001](110),[1̄1̄2](111))
では，き裂の大きさや数，配置に関わらずき裂先端に現れる負の AES 原子数が一定であることが明らかにされている．fcc-Ni では [001](010) き裂では非弾性変形時にき裂が進展し，[1̄10](001)，[1̄1̄2](111) き裂では転位射出によるき裂の鈍化が起こった．しかしいずれのモデルにおいても非弾性変形が起こるひずみがほぼ同じで，ダブルき裂モデルは単一き裂モデルのほぼ二倍の負の AES 原子を生じている．鉄の場合と同様に，き裂先端の不安定領域は一定の寸法である可能性がある．

Table 4.2 Number of negative AES atoms at the onset of inelastic deformation.
single crack model

double crack modelstrainatoms

strainatoms[001](010)

0.050

294

0.053

602

[1̄10](001)

0.044

168

0.045

266

[1̄1̄2](111)

0.030

69

0.030

137

34

4.2.3固有値による検討

図 4.5 に最小固有値 (η α(1) のうち最も小さいもの) の変化を示す．図中にはき裂の長さの変化を赤線と青線で示してある．[001](010) き裂では εyy = 0.025 近傍から最小固有値の値が前駆的に減少し始める．その後すぐにき裂長さにステップ上の変化があらわれているが，このときはまだ η α(1) は負ではない．ひずみ εyy = 0.025 以降 η α(1) は負値となる
αが，η α(1) が-40GPa(Bijεj = σiα = η α εi より) 程度になったとき，四隅のき裂が急激に成長している．その値を保った状態となるが，この点からき裂の長さが急激に伸び始めている．．[1̄10](001) き裂ではやはり×印の少し前から最小値が減少し始める．この結晶方位ではき裂開口せず転位射出するが，負値に達した後のピーク状の変化が転位の射出に対応するものと考えられる．[1̄1̄2](111) き裂では η α(1) の減少が最も急激で，き裂成長時-40〜-60GPa の一定の値を示す．図 4.6 に第二固有値が負の原子数の変化を，き裂長さの変化と共に示す．前章と同様，[001](010) き裂，[1̄10](001) き裂では負の原子数がパルス状の増加を示したのちにき裂長さが変化している．[1̄1̄2](111) き裂も他より早い段階で負の原子が観察されるとともに，穏やかにき裂長さが変化している．ダブルき裂モデルにおいてもき裂長さが変化し始める点と第二固有値が負の原子が発生する点が対応している．

35

minimum of eigenvalue
cracklength(center crack)cracklength(corner crack)

20
0

0.4

-20
0.2

-40
-60
0.05

20
0

0.4

-20
0.2

-40
-60

0.1

strain, εyy

0

0.05

strain, εyy-(a)[001](010)crack(b)[110](001)crack
minimum of eigenvalue
cracklength(center crack)cracklength(corner crack)

0.6

20
0

0.4

-20
0.2

-40

Crack length

Minimum of eigenvalue, ηα(1)

0

-60
0

0.05

strain, εyy

0.1

--

(c)[112](111)crack

Fig.4.5 Change in the minimum 1st eigenvalue and crack length.

36

0.6

0.1

Crack length

Minimum of eigenvalue, ηα(1)0.6Crack length

Minimum of eigenvalue, ηα(1)minimum of eigenvalue
cracklength(center crack)cracklength(corner crack)minimum of eigenvalue
cracklength(center crack)cracklength(corner crack)

0.6

0.2

0.05

0.1

strain, εyy

0.4

50

25

0

0.2

0

0.05

strain, εyy-(b) [110](001) crack(a) [001](010) crack
minimum of eigenvalue
cracklength(center crack)cracklength(corner crack)
75

0.4

50

25

0

0.6

0.2

0

0.05

Crack length

Number of negative ηα(2) atoms

0

Number of negative ηα(2) atoms

25

0.1

strain, εyy--(c) [112](111) crack

Fig.4.6 Change in the minimum 2nd eigenvalue and crack length.

37

0.6

0.1

Crack length

0.4

50

75

Crack length

Number of negative ηα(2) atoms

75

0minimum of eigenvalue
cracklength(center crack)cracklength(corner crack)5 厚さ方向自由度の影響

5.1 シミュレーション条件
本章では第 3 章で用いたモデルを厚さ方向に 1.5 倍して，き裂前縁方向の自由度を増やした系を解析対象とする．平板状の周期セル中央に貫通き裂を導入し，き裂寸法，結晶方位は第 3 章のものと同様の条件である．き裂部の原子を取り除いた原子総数はそれぞれ249480, 268611, 260796 である．各モデルについて，全方向周期境界の下で垂直応力が 0となるようセル辺長を調整しながら 100000fs の緩和計算を行った後，y 方向にひずみ増分 ∆εyy = 1.0 × 10−6 [/fs] を毎ステップ与えることで引張シミュレーションを行う．

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

応力-ひずみ線図と変形の様子

図 5.1 に応力ひずみ線図及び負の AES 原子数の変化を示す．図中下には比較のため第3 章の図を再掲する．いずれのモデルにおいても無負荷平衡状態では負の原子は存在してαおらず，ε = 0.03 近傍から detBij
が負の原子が現れ始める．初めて負の原子が現れた点

を応力-ひずみ曲線上に×印で，CSP で欠陥を確認したひずみを○で示す．○印から応力ひずみの非線形性が現れるという点は第 3 章，第 4 章の結果と同様である．[1̄1̄2](111) き裂の応答が，第 3 章のそれよりわずかに変化しているが，それ以外はほぼ同じである．

38

[001](010)model[-110](001)model[-1-12](111)model

Number of negative detBαij atoms

3000

Stress, σyy, [GPa]

15

10

5

0

0

[001](010)model[-110](001)model[-1-12](111)model

2000

1000

0

0.1

strain, εyy(a) stress-strain curves.

00.1

strain, εyy

α < 0 atoms.
(b) Change in the number of detBij

thick model[001](010)model[-110](001)model[-1-12](111)model[001](010)model[-110](001)model[-1-12](111)model

Number of negative detBαij atoms

Stress, σyy, [GPa]

15

10

5

0

0

0.1

2000

1000

0

0.1

strain, εyystrain, εyy(c) stress-strain curves.

0

α < 0 atoms.
(d) Change in the number of detBijbase modelα
Fig.5.1 Stress-strain curves and change in the number of negative detBij
atoms.

39

図 5.2 に引張シミュレーション中の変形の様子を，負の AES 原子を赤，それ以外を緑で着色して示す．[001](010) き裂のみぜい性的にき裂が進展し，[1̄10](001)，[1̄1̄2](111) き裂ではき裂先端に負の原子が現れた後転位が射出され転位が鈍化するという点は第 3 章と同じである．図 5.2(i) はそれぞれのモデルで最初に負の AES 原子が現れたステップでのスナップショットを示しているが，その分布は図 3.4 と一致している．しかし [1̄1̄2](111)
き裂では，第 3 章ではき裂左端の負の AES 原子は転位ループが途中で止まっていたが，厚さ方向が大きい今回の結果ではき裂左端からも完全転位として負の AES 原子が移動する様子が観察された．図 5.3 に [1̄1̄2](111) き裂での，図 5.1 の○印の前のステップ (表面原子以外欠陥がない) における負の AES 原子の分布 (左) と，その後欠陥と判定された原子のスナップショット (右) を示す．[1̄1̄2](111) き裂ではき裂前縁長さ方向の自由度を増やしたことによって，両方のき裂先端から完全転位が射出されている．

40

(i) εyy=0.035(ii) εyy=0.048(iii) εyy=0.062(a)[001](010)crack(i) εyy=0.043(i) εyy=0.026(ii) εyy=0.060(iii) εyy=0.070(ii) εyy=0.035(iii) εyy=0.049(b)[110](001)crack-(c)[112](111)crackαFig.5.2 Snapshot colored by positive/negative detBij.

41

(a)ε=0.03(b)ε=0.034

Fig.5.3 Snapshots of [1̄1̄2](111) crack negative AES atoms(left) and defect atoms
visualized by central cymmetry parameter(right).

5.2.2

負の AES 原子数と固有値

図 5.1 の×印におけるひずみと負の原子数を，第 3 章における値と共に表 5.1 に示す．いずれのモデルも負の原子が初めて現れるひずみはほぼ等しく，またその数も 1.5 倍α< 0 の原子が発生するときの([001](010) き裂のみ 1.4 倍) となっていることから，detBij

力学条件には厚さ方向自由度の影響はない．図 5.4 に最小の η α(1) の変化を，第 3 章のものと併せて示す．両者はほとんど同じで，[001](010) き裂，[1̄1̄2](111) き裂では-40GPa
で一定の値をとり，[1̄10](001) き裂では-20GPa 程度の値を示して変形が進行する．また

η α(2) が負になった原子数の変化を，第 3 章のものと図 5.4 併せて示す．[001](010) き裂
ではパルス状の増加を示した際の原子数が増えているが，それ以外はほぼ同じである．

[1̄10](001) き裂は η α(2) < 0 の原子が出現するひずみがわずかに早くなっており，また 2
番目のパルスが非常に大きくなっている．[1̄1̄2](111) き裂では，第 3 章では εyy = 0.03
近傍から η α(2) < 0 の原子が現れていたが，それが現れておらず変形後期にパルス状に現れるようになった点が異なる．このことからも，[1̄1̄2](111) き裂ではき裂前縁方向の自由42[001](010)model[-110](001)model[-1-12](111)model

Minimum of eigenvalue, ηα(1)Minimum of eigenvalue, ηα(1)[001](010)model[-110](001)model[-1-12](111)model

20

0

-20

-40
0

0.02

0.04

0.06

0.08

20

0

-20

-40
0

0.02

0.04

strain, εyystrain, εyy(a) base model(b) thick model

0.06

0.08

Fig. 5.4 Change in the minimum 1st eigenvalue.

度によって変形モードが異なる可能性がある．各モデルでの η α(2) < 0 の原子数の変化を，き裂長さ変化と併せて図 5.6 に示す．[001](010)，[1̄10](001) き裂は η α(2) < 0 の出
現とき裂長さが急増する点が対応しているが，[1̄1̄2](111) き裂はき裂長さ変化が緩慢であるため η α(2) < 0 の出現と対応していると結論づけることは難しい．

Table 5.1 Strain at emerging point of negative AES atoms.
single crack model

thick modelstrainatoms

strainatoms[001](010)

0.035

40

0.035

56

[1̄10](001)

0.042

56

0.043

84

[1̄1̄2](111)

0.026

48

0.026

72

43

[001](010)model[-110](001)model[-1-12](111)model

Number of negative ηα(2) atoms

Number of negative ηα(2) atoms[001](010)model[-110](001)model[-1-12](111)model
80

60

40

20

0

0

0.02

0.04

strain, εyy

0.06

0.08

(a) base model

80

60

40

20

0

0

0.02

0.04

strain, εyy

0.06(b) thick model

Fig.5.5 Change in the minimum 2nd eigenvalue.

44

0.08

Number of negative ηα(2) atoms
Crack length

80

60

0.4

40
0.2

20

0

0.05

strain, εyy

80

60

0.4

40
0.2

20

0

0.1

0

0.05

strain, εyy-(b)[110](001)crack(a)[001](010)crackα(2)Number of negative η
Crack length

atoms
0.6

80

60

0.4

40
0.2

20

0

0

0.05

strain, εyy

Crack length

Number of negative ηα(2) atoms

0

0.1

--

(c)[112](111)crack

Fig.5.6 Change in the minimum 2nd eigenvalue.

45

0.6

0.1

Crack length

Number of negative ηα(2) atoms

0.6

Crack length

Number of negative ηα(2) atoms

Number of negative ηα(2) atoms
Crack length

6 結論

本研究では，Ni の薄板周期セルにモード I 型き裂を導入したモデルに，分子動力学法に
ααより引張シミュレーションを行い，き裂進展挙動を原子弾性剛性係数 Bijの行列式 detBijと固有値 η α の観点から議論した．以下に本論文で得られた結果を総括する．
第 2 章では解析手法の概要を述べた．はじめに分子動力学法の基礎方程式を示し，原子間相互作用の評価に用いた EAM ポテンシャルの概要を説明した．さらに，各原子の安定性評価の指標となる原子弾性剛性係数について概説した．
第 3 章では薄板状の周期セル中央に貫通き裂を有する系を解析対象とし，[001](010)，[1̄10](001)，[1̄1̄2](111) の三種類の結晶方位について検討した．いずれのモデルにおいて
αα<0と< 0 の原子は存在しておらず，き裂近傍に detBij
も無負荷平衡状態では detBij

なる原子が生じた後，き裂先端から欠陥が発生し応力-ひずみ線図に非線形性が現れた．
αdetBij< 0 の原子が初めて出現するときのひずみはそれぞれ 0.035，0.042，0.026 で，きα裂先端近傍の原子列がそれぞれ 2，2，3 本，detBij< 0 となっている．固有値の変化を調べた結果，構造変化が起こる点 (欠陥が発生する点) では η α(1) の最小値がほぼ一定の負の値 ([001](010)，[1̄1̄2](111) は-40GPa，[1̄10](001) は-20GPa）となっていた．
第 4 章では，先の 3 つの周期セルの四隅に 4 分の 1 き裂を配置したダブルき裂モデル
αのシミュレーションを行い，周期き裂の配置を変えた影響を調べた．初めて detBij<0の原子が出現するひずみは先の単一き裂に比べてわずかに高ひずみ側 (0.002 程度) となっ
α
たが，これはき裂配置による応力遮蔽効果と考えられる．欠陥生成時の detBij< 0 原子の数はいずれも単一き裂モデルのほぼ二倍であり，き裂先端の不安定領域は一定の寸法である可能性が示唆された．固有値による検討でも欠陥が生成するときの最小の η α (1) は，

46

ほぼ先述の値 (-40GPa と-20GPa) と一致していた．
第 5 章では第 3 章のシミュレーションセルを厚さ方向に 1.5 倍とし，き裂前縁方向の自由度の影響について検討した．[1̄1̄2](111) き裂のみ，厚さ方向の自由度を増やしたことによって変形モードがわずかに変化することが，応力-ひずみ，原子配置変化，第二固有値のα< 0 の原子が初めて現れるひずみにはどのモデル
変化などから示された．ただし detBij

も変化がなく，その原子数も 1.5 倍 ([001](010) のみ 1.4 倍) となっていることから，き裂
α先端に detBij< 0 の原子が発生するときの力学条件は厚さ方向の自由度の影響はない．

欠陥生成してき裂が開口するときの第一固有値の値も先述の値と同じである．以上を総合
α
すると，結晶方位による差はあるが，き裂先端に初めて detBij< 0 が出現するときの局α
所力学状態は常に同じ (detBij< 0 の原子列＝不安定領域が一定) であり，またその後き
α裂先端から転位射出しき裂が鈍化・開口するときの Bijの第一固有値が一定の値を示すことが本研究で得られた事実である．

47

参考文献
[1] 屋代・冨田，日本機械学会論文集 (A 編),67，656,(2001)，678[2] 屋代・西村・冨田，材料，54,10,(2005)，1053[3] 屋代・藤原，材料，60,11,(2011),968[4] Yashiro, K. and Yamane, T., Proceedings of ICCM2014, ISSN2374- 3948(online),
Sientech Publisher (2014).
[5] Chang R, Int. Journ. of Fracture Mech., 6,(1970),pp111-125[6] Gumbsch P, Mater. Sci. Eng.:A,Vol.319-321,(2001),pp.1-7[7] Shimada, T., Ouchi, K., Chihara, Y. and Kitamura, T., Scientific Reports,
No.5(2015),p.8596.
[8] Yashiro, K., Tsugawa Y., and Katayama H., Advanced Structured Materials,
From Creep Damage Mechanics to Homogenization Methods, Springer, pp.557570(2015)[9] 屋代・片山・西川・吉原，日本機械学会論文集, Vol.81, No.829, 15-00271(2015)[10] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
[11] 上田顯, 分子シミュレーション-古典系から量子系手法まで-, (2003), 裳華房.
[12] Daw, M.S., and Baskes, M.I., Physical Review Letters, 50,(1983),1285[13] Daw,M.S., and Baskes,M.I., Physical Review B, 29,(1984),6443[14] Parisi, G., and Wu, Y., The Science of Sin, 24, 483, (1981)[15] Voter, A.F., Chen, S.P., Siegel, R.W., Weertman, J.R., Sinclair, R. (Eds.), Characterization of Defects in Materials, 82,(1987),175[16] Wallace, D.C., Thermodynamics of Crystals, (1972), Wiley, Newyork.

48

[17] Wang, J., Tip, S., Phillpot, S.R., and Wolf, D.,Phys.Rev.lett., 71-25, (1993),
4182.
[18] Li, Ju, Modeling Simul.Mater.Sci.Eng.，(2003)，11-173

49

謝辞
本研究を遂行するにあたり，仔細に渡る懇切丁寧な御指導を賜りました屋代如月教授(現 岐阜大学) に，心より深く感謝の意を表します．また，本研究に対して，貴重な御助言を賜りました田中克志教授，長谷部忠司准教授に，日々の研究活動において数々の励ましの言葉を頂いた片山雷太技術職員に深く感謝致します．
本研究室配属後，若輩な私に沢山のアドバイスを下さり，基礎的な事まで懇切丁寧に教えてくださった表面・界面工学研究室の諸先輩方に厚く御礼申し上げます．
また，互いに切磋琢磨しながら研究に取り組み，苦楽を共にした同輩の穴瀬正典君，池田翔伍君，岩中拓夢君，大橋正博君，佐藤博紀君，向井良多君に御礼申し上げます．
最後に，6 年間の大学生活をあらゆる面から支えて頂いた家族には，深く感謝致します．
ありがとうございました．
平成 28 年 2 月
片山 寛

50