修 士 論 文(題 目)

原子弾性剛性係数の固有値・固有ベク
トルによるき裂進展挙動解析：Fe,Mg
中のモード I き裂と介在物の相互作用

令和 4 年度

岐阜大学大学院
自然科学技術研究科 博士前期課程
物質・ものづくり工学専攻

氏名 浅井瞭

目次
1

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

12

解析手法の基礎 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1

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

4

2.2

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

5

2.3

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

6

2.4

高速化手法

6

2.5

速度スケーリング法

2.6

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

92.7

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

92.8

原子弾性係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.9

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

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

8

2.10 固有値と固有ベクトル，変形モード . . . . . . . . . . . . . . . . . . . 13
Fe 中のき裂と Ni/Co/Zr/Mo 介在物の相互作用 . . . . . . . . . . .

3

15

3.1

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

3.2

応力ひずみ線図と最終破断形態 . . . . . . . . . . . . . . . . . . . . . 17

3.3

き裂進展経路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.4

き裂進展時の原子応力と AES 固有値 . . . . . . . . . . . . . . . . . . 23

3.5

き裂−介在物間の不安定モード . . . . . . . . . . . . . . . . . . . . . 29

Mg 中のき裂と Al/Zr/Pb 介在物の相互作用 . . . . . . . . . . . . . .

4

35

4.1

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

4.2

応力ひずみ曲線とき裂進展経路 . . . . . . . . . . . . . . . . . . . . . 36

4.3

き裂進展時の原子応力 . . . . . . . . . . . . . . . . . . . . . . . . . . 43

i

4.4

き裂進展時の AES 固有値

4.5

き裂−介在物間の不安定モード . . . . . . . . . . . . . . . . . . . . . 49

. . . . . . . . . . . . . . . . . . . . . . . . 47

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

54

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

57

謝辞 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

5

ii

1

緒言
進展するき裂の先端では「材料の分離」が生じ，「物質点は常に連続している」という連続体の定義が破綻する．したがって古くから離散系の分子動力学 (MD) シミュ
レーションによる検討がなされてきた (1) ．その検討は単純なき裂にとどまらず，き裂と粒界や析出物の相互作用など，より複雑な条件での検討もなされている．Nishimuraと Miyazaki はα-Fe のき裂と粒界の MD シミュレーションを 2001 年に行っている(2)．Noronha らは MD と転位動力学 (DD) のマルチスケール的なアプローチにより

Al および α-Fe について，転位運動を阻害するブロックを設けた時の挙動を検討している (3) ．Rafii-Taber らは fcc 材料中のモード I き裂の前方にき裂に対して垂直の棒状または円柱状の原子クラスターを配置した系の MD シミュレーションを行っている (4) ．Yang らは Cu/SiC 界面上に配置したき裂の進展挙動の MD シミュレーショ
ンを行っている (5) ．Yan らはき裂と α/γ-鉄相境界の相互作用について相境界が転位やすべりを妨げることを報告している (6) ．Yu らは Ni/Ni3 Al 界面について，先のCu/SiC のような界面はく離のき裂と，Ni 相中のき裂が Ni3 Al 界面に垂直に衝突す
る MD シミュレーションを行っている (7) ．最近では Gi らは Cu/Ag バイメタル界面に沿った混合モードき裂について破壊靭性値を調べている (8) ．Yin らは Cu 析出物を含む α-Fe の [100](001) き裂の MD シミュレーションを行っている (9) ．
バイメタルやトリメタルの分野ではその強度や界面のき裂挙動の実験的な検討もなされている．Liming らはガスアーク溶接によって成形した軟鋼/シリコン青銅のバイメタル界面の引張試験を行い，界面の特性を評価している (10) ．Wengang らは電子ビーム溶解で製造した Ti48Al2Cr2Nb ブレードと Ti6Al4V のバイメタル材料の1

製造と引張試験による評価を行っている (11) ．Liang らは L415/N08825 バイメタル
複合管溶接継手の三相帯に予き裂を導入し破壊挙動を観察している (12) ．

MD 法は各原子に働く力を計算し，ニュートンの運動方程式 F = ma を，コンピュータを用いて多数の原子について解き続けることで原子運動をシミュレー
ションする方法である．MD シミュレーションは微視的メカニズムを観察するための有力なツールではあるが，原子間ポテンシャルの精度から実在の材料と全く同じではないモデル材料とならざるを得ないため，定量的な議論は不得意である．そこで我々のグループでは，局所変形の開始を定量的に議論する指標として原子弾性剛性係数 (Atomic Elastic Stiﬀness，以下 AES と称する)B αij =∆σiα /∆εj の固有値
ならびに固有ベクトルに着目した検討を行ってきた (13) ．ここで ∆σiα は原子応力，∆εj はひずみで指標 i, j は i, j = 1 ∼ 6 = xx, yy, zz, yz, zx, xy を表す．固有方程式
B αij ∆εj =η α ∆εi の固有値 η α が負となった原子を不安定とみなし，対応する固有ベク
トル ∆εi ={∆ε1 ,・・・,∆ε6 }T ={∆εxx ,・・・,∆γxy }T から変形モードを議論する．進展するき裂先端の原子の不安定モードについては Fe(14) ，Si(13) ，Mg(15) ，SiC(16) そして異種金属界面 (17) のき裂について検討してきた．異種金属界面のき裂では，表面エネルギーの差から界面が分離するようにき裂が進展することはなく，表面エネル
ギーが低い=弾性係数が小さい (18) 金属側を進展することを示し，AES の第一固有値 η α(1) がそれぞれの金属相の相対的な「硬さ」（弾性限界への達しにくさ）を表していることを報告している．
本研究では，き裂と析出物との相互作用について AES の固有値ならびに固有ベ
クトルの視点から議論するために，3 章では Yin らのシミュレーションと同様に Fe

2

中のモード I き裂の前方に原子クラスター (介在物) を配置した系の引張シミュレー
ションを 4 種の元素 (Ni,Co,Zr,Mo) について行った．4 章では Fe を Mg に変えて 3 種の介在物 (Al,Pb,Zr) について同様の引張シミュレーションを行った．なお，Yin らの研究と異なり，実在する材料を想定したものではなく，前報 (17) と同様に弾性係数の大小関係から介在物元素を選んでいる．

3

2

解析手法の基礎2.1

分子動力学法

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

(2.1)

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

Fα = −

∂Etot
∂r α

(2.2)

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

α

(2.3)
(2.4)

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

dr α= v α (t)
dt

(2.5)

であり，式 (2.1) と式 (2.5) を式 (2.3) と式 (2.4) に代入すると，((∆t)2 F α (t)3)
r (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α
α

α

α

4

(2.6)
(2.7)

となる．両式の和と差をとると，()
F α (t)+ O (∆t)4αm()
r α (t + ∆t) − r α (t − ∆t) = 2∆tv α (t) + O (∆t)3r α (t + ∆t) + r α (t − ∆t) = 2r α (t) + (∆t)2

(2.8)
(2.9)

が得られる．∆t3 以上の高次項は無視できるとすると，時刻 t + ∆t での位置ベクトルと t での速度は

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

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

2.2

原子間ポテンシャル

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

5

2.3

EAM ポテンシャル

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

Etot =∑F α (ρα ) +α1 ∑ ∑ αβ αβϕ (r )
2 α(2.12)

β(̸=α)

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

β(̸=α)

で評価する．
本研究では，Zhou らが多数の元素に対してパラメータを提案しているポテンシャ
ルを用いた (22)–(24) ．

2.4

高速化手法

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

6

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

Fig. 2.1

Schematic of domain decomposition method.

7

2.5

速度スケーリング法

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

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

ここで mα は粒子 α の質量，viα は 温度 T での粒子 α の速度の i 方向成分，kB はBoltzmann 定数で kB =1.38 × 10−23 [J/K] である．したがって，目標の温度 T0 における原子 α の速度成分を viα0 とおくと viα0 は√
viα0 =

3kB T0mα(2.15)となる．同様に，温度 T の時の原子 α の速度成分は

√
viα =

3kB T
mα

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

viα0
=
viα

√

T0
T

(2.17)

つまり，系の温度を T から T0 にするには，式 (2.17) の右辺を現在の速度に掛けてやればよい．ただ，これだけでは数値積分に反映されないので，Verlet 法における式 (2.10) で t →t + ∆t の変化分である次式を√

T0 /T 倍して数値積分する．

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

α

α

α

α

2F

α

(t)
(2.18)mα

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

8

2.6

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

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

σij

1
=
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

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

(
Bijkl ≡

∆σij
∆εkl

)

= Cijkl + (σil δjk + σjl δik + σik δjl + σjk δil − σij δkl − σkl δij )/2(2.20)
ここで，δij はクロネッカーのデルタである．Wang, Yip らは，ひずみの対称性を考
sym
慮したテンソル Bijkl≡ (Bijkl + Blkji )/2 の正値性によって結晶の安定性を評価することを提案している (27), (28) ．

2.7

原子応力

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

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

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

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

(2.22)と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子α と 原子 β の間の距離 rαβ には(rαβ )2 = r̄iαβ Gij r̄jαβ(2.23)なる関係が成立する．ただし，Gij = Jki Jkj である．ここで i，j は自由指標，k は総和規約に従うダミー指標であることに注意されたい．仮想変形の Lagrange ひずみテンソル ηij は]
1[
Gij − δij2

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

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

∂rαβ
∂ηij

r̄iαβ r̄jαβ=
rαβ(2.26)
10これより EAM ポテンシャルにおける応力は次式で評価される．

σij

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

N ∑
N (
∑
α β(̸=α)

( ∑
)
N ∑
N
1
∂rαβ ∂Etot
2 α∂ηij ∂rαββ(̸=α)) αβ αβri rj
1 ′ αβ αβ′α α′ β αβϕ (r ) + F (ρ )f (r )2rαβ

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

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

1 ∑
Ω
N

σijα =(

β(̸=α)

2.8) αβ αβri rj
1 ′ αβ αβ′α α′ β αβϕ (r ) + F (ρ )f (r )2rαβ

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

1 ∂ 2 Etot
V ∂ηij ∂ηkl

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

Cijkl

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

+

N
∑
α

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

riαβ rjαβ
f ′ β (rαβ ) αβ
r
11

}{ ∑
β(̸=α)

rαβ rαβf ′ β (rαβ ) k αβlr}](2.28)

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

[
} αβ αβ αβ αβ
N {
1 1 ∑
ϕ′ αβ (rαβ ) ri rj rk rl′′ αβ αβ=ϕ (r ) −Ω 2rαβ(rαβ )2β(̸=α){} αβ αβ αβ αβ
f ′ β (rαβ ) ri rj rk rl′′ β αβ+F (ρ ) f (r ) −rαβ(rαβ )2β(̸=α)]{∑αβ αβ }{ ∑αβ αβ }rr
r
r
ij
+ F ′′ α (ρα )f ′ β (rαβ ) αβ
f ′ β (rαβ ) k αβl(2.29)
r
r
N
∑

′α

α

β(̸=α)

β(̸=α)

原子弾性剛性係数

2.9以上で定義した原子応力，弾性係数から，Voigt 対称性をもたせた原子弾性剛性係数を以下で求める．
α
α
α
α
α
Bijkl
= Cijkl+ (σilα δjk + σjl
δik + σik
δjl + σjk
δil − 2σijα δkl )/2(2.30)αBijkl
は 4 階のテンソルであるが，Voigt 対称性を持たせた場合，弾性係数テンソルと同じく独立な成分は 21 個となり，xx，yy，zz，yz，zx，xy を 1∼6 とする Voigt 表
記を用いれば 6 × 6 のマトリックスとして Bijα と表すことができる．その場合の式 (2.30) の成分は以下のようにマトリクスで表現できる．
 α
α
α
α
α
α
α

[

Bijα

]

[
=

Cijα

]


1
+ 
2


2σ1

−σ1 − σ22σ2α

−σ1 − σ3−σ2α − σ3α
2σ3α

sym.

12

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

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



σ6α
σ6α

−σ6α 


σ5α

σ4α
α
α
σ1 + σ2(2.31)

2.10

固有値と固有ベクトル，変形モード

固有方程式





 ∆ε1  ∆ε1 

[ ]
.
.
α
α
.
.
Bij
. =η  . 





∆ε6
∆ε6

(2.32)

を解くことは Bijα を対角化することと等価である．

 α(1)η[ ]
αBij ⇒ 0η α(2)00η α(3)

0
0
0

0
0
0
0

η α(4)00000η α(5)



(2.33)

η α(6)

sym.

せん断成分は互いに独立なので，6×6 の Bijα マトリクスは通常は次のような形となる．

 α
α
α
B11

B12

B13

0

0

α
α
α
0
0
B23
B22
B21
[ ] 
α
α
α

0
0
B
B
B
α
33
32
33
Bij = 
α
 0
0
0
0
B44
 0
α
0
0
0
B55

0

0

0

0

0

0
0 
0 

0 

0
α
B66

(2.34)αα
したがって B44
∼ B66
は固有方程式 B αij ∆εj =η α ∆εi の解で，それぞれの固有ベクト

ルは {0, 0, 0, ∆γyz , 0, 0}，{0, 0, 0, 0, ∆γzx , 0}，{0, 0, 0, 0, 0, ∆γxy } である．垂直ひずみに対する 3×3 の部分マトリクスの固有値は，無ひずみの等方弾性体では特性方程式

C11 − ηC12
C12

(2.35)det 
C11 − ηC12  = 0sym.
C11 − ηの解から η=C11 −C12 (重解)，C11 +2C12 となる (C11 , C12 は弾性係数)．この C11 +2C12に相当する固有値は他の固有値より大きく (第 6 固有値)，負になることはない．等
α
α
との大小関
∼ B66
方でない場合，C11 − C12 に相当する固有値は重解とならず，B44

係で第 1∼5 固有値のいずれかとなる．(2.32) 式の固有ベクトルの ∆ε4 ∼ ∆ε6 は工学
せん断ひずみ ∆γyz ，∆γzx ，∆γxy なので，1/2 して


∆ε
∆ε
∆ε
xx
xy
xz
[
]


∆εij = ∆εyy ∆εyz 
sym.
∆εzz

13

(2.36)

の 3×3 テンソルとし，これの主軸，すなわちこのテンソルの固有ベクトル x1 , x2 , x3を求めることができる (固有値 ∆ε1 < ∆ε2 < ∆ε3 の順)．本論文では最大せん断方向の x1 + x3 を変形モードとして定義する．

14

3

Fe 中のき裂と Ni/Co/Zr/Mo 介在物の相互作用

3.1

シミュレーション条件

図 3.1 に示すように，中央にモード I き裂を有する Fe 完全結晶の薄板セルにおいて，円状領域の Fe 原子を介在物 (Ni,Co,Zr,Mo) に置換したセルに全方向周期境界を適用した系を解析対象とした．図では Fe 原子を青色，介在物原子を黄色，初期き裂の上下の原子をそれぞれ緑色，赤色で示している．結晶方位は x，y，z 方向は [1̄1̄2]，[11̄0]，[111] であり，薄板セルの大きさは 50 nm×40 nm で厚みは 3 nm とした．き裂
上下の原子 (緑と赤) 間には互いの力を無視する仮想遮蔽シートを設けている．遮蔽
シートの横幅は 1/3Lx (x 方向セル長さ) で，y 方向の厚さは y=0.5Ly の原子層を中心にカットオフ半径の 2 倍とした．セルの y 方向長さを基準としたため緑の層は 3 原
子分，赤の層は 2 原子分と上下対称ではないが，予き裂導入部であるため重要ではない．介在物は遮蔽シートの端部より 6.66 nm 離れた位置の，半径 1.5 nm の Fe 原
子を別の原子種に置換して作成した．したがって介在物の結晶構造・格子長さはその元素の安定構造ではないことに注意されたい．総原子数は 506088 である．原子間相互作用には Zhou(22) らの EAM ポテンシャルを用いた．熱による影響を排除する
ため，温度は極めて低温 (0.1 K) とし速度スケーリング法により制御した．垂直応力が 0 になるようにセル寸法を調整しながら 10 ps の緩和ステップ計算を行った後，x，z 方向のセル長さを固定して y 方向のみに引張ひずみを毎ステップ ∆εyy =1.0 × 10−7増加させる引張シミュレーションを行った．

15

2rcutthickness : 3nm

Fe

Ly = 40nm

6.66nmLx / 3

R = 1.5nm
Lx / 3yx

Lx = 50nm

Fig. 3.1

Simulation model.

16

3.2

応力ひずみ線図と最終破断形態

図 3.2 に応力ひずみ曲線を示す．(a) は比較のために行った介在物なしの系の結果
で，(b)〜(e) がそれぞれ Ni，Co，Zr，Mo を介在物に含む系である．また，図 3.3 に原子種と各原子の AES の第一固有値 η α(1) の正負で色分けした最終破断時のスナッ
プショット (片側のき裂先端を拡大) を示す．ただし Zr 介在物の (d) はき裂が停滞し，

εyy =0.1 の範囲で破断しなかったので εyy =0.127 の図を示している．介在物のない(a) では応力ピーク後にき裂が不安定成長して周期境界下でつながり破断し，応力がほぼ 0 となる．介在物が存在する場合，(a) と同様に最初の応力ピーク後にき裂が左右に広がるが，すぐに介在物にぶつかり停滞し，介在物を貫通/回避する挙動に応じて鋸刃状の応答を示しながら (d) の Zr 以外は最終的に破断して応力が 0 となった．

(d) は最初のピークから応力急減がほとんどなく，εyy =0.1 までの引っ張りでは応力が減少せず破断しなかった．最終破断の図 3.3 を見ると (b)，(c) はき裂が介在物の中を貫通して破断し，(e) は介在物を回避して破断している．(d) は介在物がき裂先
端を大きく鈍化させて進展を阻止している．(e) は Ni や Co と比べると最初のピー
ク時のひずみ・応力が他に比べて小さいが，これについては後のき裂先端の応力分
布で述べる．

17

6

6

6

(e)
(f)

4

2

0

0.05

0.1

4

(a)
(d)
(b)

2

0

0.05

(a) No obstacle

4

(b)

(d)

2

0

0.1

(c)
(a)

0.05

Strain, εyyStrain, εyy(b) Fe / Ni(c) Fe / Co

Strain, εyy

0.1

6

6

(e)

(a)

Stress , σ yy , GPa(b)Stress , σ yy , GPa(e)Stress , σ yy , GPa

Stress , σ yy , GPa

Stress , σ yy , GPa

(c)

(d)

4
(c)

2

0

0.05

0.1

(a)
(d)
(c)

4

(f)

(b)

(e)

2

0

0.05

0.1

Strain, εyyStrain, εyy(d) Fe / Zr(e) Fe / Mo

Fig. 3.2(a) No obstacle,= 0.033

Stress-strain curves.

(b) Fe/Ni,= 0.069(c) Fe/Co,= 0.093Fe,
Fe,
Ni, Co, Zr, Mo,
Ni, Co, Zr, Mo,(d) Fe/Zr,

Fig. 3.3= 0.127(e) Fe/Mo,= 0.101

Snapshots of final fracture morphology. Atoms are colored by the atom
type and the sign of the 1st eigenvalue η α(1) .

18

3.3

き裂進展経路図 3.4 は Ni 介在物の応力ひずみ曲線 (図中に再掲) の点 (a)〜(f) における，左側の介在物周辺のスナップショットを原子種と η α(1) の正負で 4 種類に色分けして示したものである．ただし η α(1) については後で議論するためここでは参考程度に示す．なお，基本的には左右の介在物で挙動（回避や貫通）に大きな差はなく，き裂が介在物に到達するのもほぼ同じタイミングであるが，Ni の場合は左側が先に貫通したた
めここでは左側を示している．他の元素についても先に進展した方を示す．(a) が応力ピーク点で，この後き裂が Fe 中を進展し，Ni 介在物の中に進入して応力急減が止まる (図 3.4(b))．(b)→(c) は応力が再び上昇して抵抗を示しているが，この時き裂先
端は介在物の中で鈍化している (図 3.4(c))．(c)→(d) ではわずかに応力低下するが，これはき裂が介在物中を１格子分進展してすぐに停止したためで，介在物中の「安定成長」である．(e)→(f) ではき裂先端が赤い点線で示した方向にすべり変形しながら進展した．(f) では介在物の前方の界面からボイドが発生し，そこから Fe マトリ
クス中をき裂が不安定進展して周期境界下で反対側のき裂と合体して破断した (図3.3(b))．図 3.5 は Co 介在物とき裂の相互作用のスナップショットであり左側の介在
物周辺を拡大している．先の Ni 介在物と似た進展挙動を示すが，界面からき裂を生じることなく介在物を完全に分離するように進展した．応力-ひずみ曲線の (b)→(c)は介在物中でのき裂の安定成長，(c)→(d) は左のマトリクス界面までの不安定成長
で，(e) 点はマトリクス内への不安定成長開始点に対応する．図 3.6 に Zr 介在物のき裂進展過程について右側を拡大して示す．Zr 介在物は初期に配置した bcc 構造から乱れてアモルファスに近い状態となっている．(a)→(b) でき裂がマトリクス中を進19

展する時も応力上昇しており安定成長しながら Zr 介在物に達する．その後一定の応力で介在物中心まで進入し ((b)→(c)), き裂先端が大きく開口し鈍化した ((d)→(e))．
鈍化したき裂先端には Zr 原子がキャップのように存在している．図 3.7 は Mo 介在
物のき裂との相互作用である．最初の応力ピークでき裂が進展するが，介在物から離れた位置で停滞し，先端が 120 °の多角形に鈍化した ((b)→(c)). (c)→(d) はき裂が進まずに開口量が増加している．その後，(e) で 120 °の先端から新しいき裂が発生する形で介在物に到達するとともに，介在物上側界面にボイドが発生した．右図(e)→(f) で応力上昇して最初のピークより高い応力を示したときは，主き裂の開口とともに介在物上部のボイドが成長している．その後，ボイドがき裂先端と合体して介在物を回避するように不安定成長して破断した．介在物は塑性変形することなく初期形状を維持している．また，ボイドは Fe/Mo 界面そのものではなく Fe 側の格子間で生じており，Mo 介在物の周りに Fe が一層だけ残っている．

20

6

(a)

(b)

= 0.025= 0.035(c)

Stress , σ yy , GPa(c)= 0.044

4

(e)
(f)

(a)
(d)
(b)

2

0

0.05

Strain, εyy(d)= 0.047(e)(f)= 0.062

Ni,

Ni,

Fe,

Fe,

Fig. 3.4= 0.056

Interaction between crack and Ni obstacle in Fe matrix. Atoms are colored by the atom type and the sign of the 1st eigenvalue η α(1) .

(a)

(b)

= 0.025= 0.036(c)= 0.055

6

Stress , σ yy , GPa

(e)

4

(c)
(a)

(b)

(d)

2

0

(d)

Fe,

Fig. 3.5= 0.068(e)Fe,= 0.083

Co,

0.05

0.1

Strain, εyy

Co,Interaction between crack and Co obstacle in Fe matrix. Atoms arecolored by the atom type and the sign of the 1st eigenvalue η α(1) .
21

0.1

(a)

(b)

= 0.025= 0.038(c)= 0.065

6(e)Stress , σ yy , GPa

(b)
(a)

(d)

4
(c)

2

0

(d)

= 0.082(e)Fig. 3.6

0.1

Zr,

Zr,

Fe,

Fe,

0.05

Strain, εyy= 0.1Interaction between crack and Zr obstacle in Fe matrix. Atoms are colored by the atom type and the sign of the 1st eigenvalue η α(1) .

6

(a)

(b)

= 0.020= 0.023(c)= 0.040

Stress , σ yy , GPa

(a)
(d)
(c)

4

(f)

(b)

(e)

2

0

0.05

Strain, εyy

(d)

Fe,

Fig. 3.7= 0.052

(e)

Fe,

(f)

= 0.057= 0.077

Mo,

Mo,

Interaction between crack and Mo obstacle in Fe matrix. Atoms arecolored by the atom type and the sign of the 1st eigenvalue η α(1) .
22

0.1

3.4

き裂進展時の原子応力と AES 固有値

表 3.1 に各元素の完全結晶における最近接原子間距離および弾性特性を示す．また，
図 3.8 にき裂進展前の同一ひずみにおける右のき裂先端と介在物のスナップショット
αを，引張方向の原子応力 σyy
で着色して示す．暖色が引張応力を，寒色が圧縮応力を示す．(a) の Fe/Ni は表 3.1 の最近接原子間距離はほぼ同じであるが，bcc に配置された Ni 介在物中には不均一な応力分布が見られる．Ni は bcc 構造を保っている
が，安定でない構造のため局所的にひずんでいると考えられる．(b) の Co 介在物はbcc に配置されても介在物内部の応力はほぼ均一で，周囲の Fe に対してわずかに引張応力を示している．(c) の Zr 介在物は格子長さの差が大きく，緩和計算の時に bcc
格子点から膨張してアモルファス化し，Fe マトリクス内にも転位を生じている．Fe
マトリクス中の青/黄のコントラストが強い部分が転位芯である．(d) の Mo は同じ

bcc だが最近接原子間距離は 0.273 nm で Fe より少し大きいので，介在物内は均一な圧縮状態となっている．また，y 方向の引張ひずみによって介在物の赤道方向周
囲に引張応力を生じているが，き裂先端の応力場と重なり合ってき裂−介在物間に一様な引張応力が見られる．このためき裂進展が開始するひずみが介在物なしや他の介在物に比べて小さい (図 3.2)．図 3.9 は図 3.8 の応力分布の，y=0 のき裂面上のαをプロットしたものである．なお，介在物が無い Fe き裂周りの分
原子応力分布 σyy

布も合わせて示している．横軸はき裂先端からの距離を示しており，図中の実線は√線形破壊力学による無限体中の単一き裂の応力分布 σyy =KI / 2πr である．青線はき裂の表面エネルギー γs とヤング率 E，ポアソン比 ν から算出される Griﬃth き裂の破壊靭性値 KIC =√2Eγs /(1 − ν 2 )=22 MPa・m1/2 をプロットしたものであるが，23

本き裂はそれよりも大きなき裂先端応力を示しており，KIC による予測は大きく外れる．赤線は (a) の応力分布に合う KI (=90 MPa・m1/2 ) をプロットしたもので，介√在物を有する系もき裂近傍の応力分布は赤の 90/ 2πr の分布に沿っているが，(e)の Mo 介在物は介在物界面でも高い引張応力となりき裂-介在物間の応力が全体的に高くなっている．このため他に比べ低ひずみ・低応力で最初のき裂進展が起こったものと考える．

24

Table 3.1

element
structurenearest neighbordistance [nm]Young’s modulus [GPa]1st eigenvalueη α(1) [GPa]Properties in bulk perfect lattice.

Fe
bcc

Ni
fcc

Co
hcp

Zr
hcp

Mo
bcc

0.248

0.249

0.251

0.320

0.273

129.0
94
(C11 − C12 )

136.3
99
(C11 − C12 )

251.0
65
(C44 )

121.7
25(C44 )

367.7
113(C44 )

40 GPa40 GPa-2 GPa-5 GPa(a) Fe/Ni,(c) Fe/Zr,
Fig. 3.8(b) Fe/Co,= 0.022= 0.02250 GPa40 GPa-30 GPa-15 GPa= 0.022(d) Fe/Mo,= 0.022αaround the right crack edge and obstacle.
Distribution of σyy

25

30

30

crack

90/(2πr)1/2

20

10

0

Ni

20

10

10

0

0

22/(2πr)1/2
0

5

10

0

Distance from cracktip, r, nm

crack

5

10

Distance from cracktip, r, nm

Mo

crack

Zr

0

yy, GPa

10

0

Stress,Stress, σ yy, GPa

20

10

-10

0

-10

-20

0

5

10

Distance from cracktip, r, nm(d) Fe / Zr
Fig. 3.90510Distance from cracktip, r, nm(e) Fe / Mo
Stress distribution on y=0 crack plane at εyy =0.022.

26

5

10

Distance from cracktip, r, nm(c) Fe / Co(b) Fe / Ni(a) No obstacle
20

Co

crack
20

Stress, σ yy, GPa

Stress, σ yy, GPa

Stress,yy, GPacrack

図 3.10 は図 3.8 と同じスナップショットを各原子の η α(1) の値で着色して示したものである．また，表 3.2 にシミュレーションセル中の全ての Fe 原子と介在物原子でそれぞれ平均した η α(1) の値を示した．図では η α(1) が小さい値，すなわち不安定側の原子を暖色で示している．無負荷完全結晶での η α(1) の値は先の表 3.1 に示している (括弧内は最小固有値となる弾性係数成分を表す) が，引張ひずみの増加によって全体的に低下し，Fe 母相はいずれも約 48 GPa と半分程度となっている．き裂が介在物に進入する Ni，Co，Zr はいずれも周囲の Fe より η α(1) が小さくなっており，Fe母相に比べて「柔らかい」（塑性変形が生じやすい）状態にある (15) ．(d) のみ介在物の η α(1) が Fe 母相より高い値を示し「硬い」状態にある．

27

60 GPa

60 GPa0 GPa

15 GPa(a) Fe/Ni,= 0.022(b) Fe/Co,= 0.022110 GPa

60 GPa-10 GPa-50 GPa(c) Fe/Zr,= 0.022

Fig. 3.10(d) Fe/Mo,= 0.022

Distribution of η α(1) around the right obstacle.

Table 3.2 Average of η α(1) in Fe matrix and obstacle at εyy =0.022.
Element
Fe
Ni
Fe
Co
Fe
Zr
Fe
Mo
α(1)
Average η[GPa] 48.4 10.7 48.5 10.8 47.5 27.1 48.1 100.4

28

3.5

き裂−介在物間の不安定モード

図 3.11 はき裂先端が Ni 介在物に到達した瞬間と，介在物内に進入した状態のスナッ
プショットを，η α(1) の値で着色したもの ((a),(d))，原子種と η α(1) の正負で色分けしたもの ((b),(e)) で，さらに η α(1) < 0 の原子について，固有ベクトル {∆εxx ,・・・,∆γxy }Tから求めた最大せん断方向 (15) を η α(1) の大きさによってスケーリングした矢印で xy
面に投影した図 ((c),(f)) である．最大せん断方向の矢印の始点は (b)，(e) スナップ
ショットの赤または水色の原子位置であり，スナップショットの色を薄めてその上に矢印を描いている．前報 (14) で示したように進展するき裂先端には η α(1) < 0 の原子が現れるが，図 (b) でも介在物に達したき裂先端に赤色の原子が認められる．(c),(f)の図では Ni 介在物の η α(1) < 0 原子の不安定モードが大きく様々な方向を向いているが，(b),(c) のき裂先端の赤原子，(e),(f) のき裂先端の水色原子はき裂開口のモードを示している (赤矢印)．固有値の正負ではなく η α(1) の値 ((a),(d)) を見ると，介在物内でき裂進展に伴うすべり変形が生じると介在物内部の η α(1) が回復し，(a) に比べ (d) では Ni 内部が青色に安定化している．図 3.12 は Co 介在物についての図である．図 (b) においてき裂先端から Co 介在物中に，x 軸に対して ±60 °方向に負の固有値を持つ水色の原子列が認められる．η α(1) < 0 原子の変形モード (図 (c)) は，Fe原子の赤矢印が大きく，また Co 介在物の η α(1) < 0 原子のモードは 60 °の原子列に対して垂直に近い．変形モードがすべり面 (原子列) に平行でないのは，連続体ではなく離散系のため，Co の中ですべりを生じる際に図 (c) 上の模式図に示すような原子の凹凸 (パイエルスポテンシャル) を波打つように移動する必要があり，それが反映したモードと考える．介在物内にき裂が進入すると介在物中の η α(1) の値が回復す29

る (図 (d))．図 (c),(f) ではき裂先端から少し離れた ±60 °方向に大きな変形モード(矢印) が見られ，この時点ではき裂の開口よりも介在物中の塑性変形が優位である
ことを示唆している．図 3.13 に Zr 介在物の図を示す．図 (b) のき裂先端に負の固有値を持つ赤い原子が存在するが，そのモードは介在物や界面に存在する η α(1) < 0 原子に比べると小さく (c) 図では矢印として見えない. 介在物と母相の境界は互いに混ざりあっており，Fe の η α(1) < 0 原子が多く見られる．図 (e) では，き裂表面に沿うように黄色の Fe 原子が介在物原子の中に存在していることから，「割れ」により進入したのではなく介在物表面を「押し込む」形でき裂が進展したものと推測される．
アモルファス化した介在物の Zr 原子の η α(1) は負に近い暖色のままで，Ni，Co のよ
うに回復しない．図 3.14 は Mo 介在物の図で，先端が 120 °に鈍化した時と，その後頂点から新たに開口して介在物に達するとともに上の界面からボイドが発生している時について示した．(a) で介在物の界面と周囲に η α(1) が低い原子が見られ，介在物界面から転位を生じている．図 (b) で η α(1) が負なのはき裂前方の Mo 原子 (水色) と上の界面の Fe 原子でこれらの変形モードは (c) に示したとおりである．後者の Mo 原子のモードはき裂表面へ引き寄せられる方向であり，圧縮状態にある Mo
介在物がき裂の接近によってつながろうとするモードを生じたものと考える．その後，実際に 120 °のき裂先端から新たなき裂が発生して介在物に達する．一方，上の界面 Fe のモードは，下の図 (f) に破線¬で示した Fe と介在物の境界面のすべりを誘発し，ボイドが発生した．この界面すべりのモードが図 (f) のき裂と介在物とすべりの三交点に認められる．

30

Fig. 3.11

Snapshots of crack tip just reaching and penetrating Ni obstacle. Left:colored by η α(1) . Middle: colored by atom type and the sign of η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

31

Fig. 3.12

Snapshots of crack tip just reaching and penetrating Co obstacle. Left:colored by η α(1) . Middle: colored by atom type and the sign of η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

32

Fig. 3.13

Snapshots of crack tip just reaching and penetrating Zr obstacle. Left:colored by η α(1) . Middle: colored by atom type and the sign of η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

33

Fig. 3.14

Snapshots of crack tip just reaching and bypassing Mo obstacle. Left:colored by η α(1) . Middle: colored by atom type and the sign of η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

34

4

Mg 中のき裂と Al/Zr/Pb 介在物の相互作用

4.1

シミュレーション条件
2rcut

thickness : 3nm

Ly = 40nm

Mg
6.66nmLx / 3

R = 1.5nm
Lx / 3yx

Lx = 50nm

Fig. 4.1

Simulation model.

前章と同様に，中央にモード I き裂を有する Mg 完全結晶の薄板セルにおいて，円状領域の Mg 原子を介在物 (Al,Zr,Pb) に置換したセルに全方向周期境界を適用した系を解析対象とした (図 4.1)．Mg の結晶方位は x [1̄1̄20]，y [0001]，z [101̄0] であり，薄板セルの大きさは 50 nm×40 nm で厚みは 3 nm とした．総原子数は 240240 である．温度，緩和条件，y 方向へのひずみ速度 ∆εyy =1.0 × 10−7 は前章と同じにした．

35

4.2

応力ひずみ曲線とき裂進展経路

図 4.2 に応力ひずみ曲線を示す．(a) は比較のために行った介在物なしの系の結果
で，(b)〜(d) がそれぞれ Al，Zr，Pb を介在物に含む系である．(b) と (d) の応力ピー
クは介在物なしのものと比べてひずみと応力が大きくなっており，介在物によってき
裂進展が遅延したものと考える．介在物がない場合は周期境界下でき裂がつながり応力が 0 となるが，介在物がある系はいずれも一度き裂進展が止まり 0.5 GPa 程度の応力で低下が止まる．その後 (b) と (d) は εyy =0.1 まで続いており，ほぼ一定の応力で変形が進行するが，(c) では εyy =0.045 付近を境に再び緩やかに減少して εyy =0.08
を超えたあたりで応力はほぼ 0 となる．図中に示した (a)〜(d) 点のスナップショッ
トを図 4.3∼ 図 4.5 に示す．

36

1.5

1.5

Stress,σ ,GPaStress,σ ,GPa

(a)

1

1
(c)

0.5

0.5

(b)
(d)
0

0.02

0.04

0.06

Strain,ε

0.08

0

0.1

0.02

0.06

0.08

0.1

(b) Mg / Al(a) No obstacle
1.5

1.5
(a)

Stress,σ ,GPa(a)Stress,σ ,GPa

0.04

1
(c)
(d)
0.5

(b)

1
(d)

0.5
(c)

(b)
0

0.02

0.04

0.06

Strain,ε

0.08

0.1

0.02

0.04

0.06

Strain,ε(d) Mg / Pb(c) Mg / Zr
Fig. 4.2

0Stress-strain curves.

37

0.08

0.1

図 4.3 は Al 介在物を含む系のスナップショットである．上の 4 枚は引張方向の原α
子応力 σyyで着色したものであり，上部に示しているカラーバーの閾値で暖色になるほど高い引張応力を示している．下の 4 枚は原子種と各原子の AES の第一固有値

η α(1) の正負で色分けしている．前章同様 AES の正負はここでは参考程度に示す．また，前章では右介在物付近を拡大して示したが，本章では介在物周辺以外での転位や双晶変形が重要となるためセル全体を示した．(a) はき裂が進展し始める時のス
ナップショットで，き裂が介在物に到達するとすぐに (b) の上図に黄色線で強調した
ように介在物とき裂先端から pyramidal 面に双晶変形の帯が形成される．左側の双晶変形は介在物の上側から発生して介在物上部で屈折してから左下方向に成長して，周期境界下で右からの双晶とつながる．双晶変形によってき裂先端の形は台形となる．右側のき裂先端は介在物を押し込むようにしてき裂先端が鈍化しており，介在
物はつぶれるように変形している (下図 (c))．さらに引っ張ると双晶変形が形成された下側に介在物を押しつぶしながら回避して進展した．応力はき裂先端が大きく鈍化したためか，あまり減少していない．また，き裂が大きく開口したことによって双晶変形の帯の幅が広がっている (図 4.3(c)→(d) 黄線)．図 4.4 に Zr 介在物を含む系のスナップショットを示す．応力ピークでき裂進展が開始してから (b)，(c) のように下側に双晶変形の帯が形成され，右側のき裂先端が双晶変形を生じた下側に回避した．(a)(b) の下図では青い部分が介在物を示すが全く変形しておらず初期形状のままでき裂先端のみが変化していることが確認できる．図 4.2(c) で応力は (b)→(c) で波を打ちながら徐々に増加しているが，き裂が介在物を回避した後は緩やかに減少した．図 4.5 に Zr 介在物を含む系のスナップショットを示す．(a) における応力ピー

38

クでは上図 (a) の拡大図に示すようにき裂先端が原子 1 個分き裂が開口して鈍化しており，き裂進展が遅れている．この開口は図 4.2(d) の応力ひずみ曲線の応力ピー
ク前の赤矢印で示した段差の時に生じている．さらにき裂先端は原子 2〜3 個分開口して鈍化する．(b) で右側のき裂先端が介在物に達するが，左側の介在物から周期境
界を通して右側の介在物に線状欠陥が発生している．(a)→(b) の応力変化はこれらのき裂鈍化や欠陥発生による．この欠陥は下図の η α(1) で見ると η α(1) が負 (赤原子)となっている．(c) では線状欠陥は双晶変形へと変化して，右側のき裂先端が介在物の中に進入している (下図 (c))．この (b)→(c) の過程で応力が急減した．(c)→(d) はき裂が介在物の中を安定成長している．εyy =0.1 以降は前章の Fe/Zr のようにき裂が開口するとともにき裂先端を覆うように介在物が大きく変形した．

39

twin deformation
-3GPa

(a)

(b)

= 0.023Mg,(a)= 0.023

Fig. 4.3

(b)

13GPa

= 0.037

(c)

Mg,

Al,

= 0.037(c)= 0.054

(d)

= 0.1

(d)

= 0.1

Al,

= 0.054

α
Top:Snapshots colored by σyycorresponding to Points (a)∼(d) in
Fig.4.2(b). Bottom:Snapshots colored by the atom type and the signof the 1st eigenvalue η α(1) .

40

-3GPa(a)(b)= 0.024Mg,(a)= 0.024

Fig. 4.4

(b)

15GPa

= 0.031

(c)

Mg,

Zr,

= 0.031(c)= 0.045(d)= 0.055(d)= 0.055

Zr,= 0.045α
Top:Snapshots colored by σyycorresponding to Points (a)∼(d) in
Fig.4.2(c). Bottom:Snapshots colored by the atom type and the signof the 1st eigenvalue η α(1) .

41

-4GPa(a)(b)= 0.028Mg,(a)= 0.028

Fig. 4.5

(b)

7GPa

= 0.039

(c)

Mg,

Pb,

= 0.039(c)= 0.049

(d)

= 0.1

(d)

= 0.1

Pb,

= 0.049

α
Top:Snapshots colored by σyycorresponding to Points (a)∼(d) in
Fig.4.2(d). Bottom:Snapshots colored by the atom type and the signof the 1st eigenvalue η α(1) .

42

4.3

き裂進展時の原子応力表 4.1 に各元素の完全結晶における最近接原子間距離および弾性特性を示す．また，図 4.6 にき裂進展が開始する直前の応力ピークにおける右のき裂先端と介在物
αのスナップショットを，引張方向の原子応力 σyy
で着色して示す．暖色が引張応力を，寒色が圧縮応力を示す．(b) の Mg/Al は介在物との界面において介在物側に高い引張応力，Mg 母相に圧縮応力がかかっている．Mg の最近接原子間距離は表 4.1
より 0.32 nm で Al が 0.29 nm なので，Al 介在物はやや引張となる．周囲の Mg 母相は介在物の子午線方向に引張，赤道方向に圧縮となっている．介在物からき裂先端にかけて介在物高さで原子構造が変化しているように見えるが，後述のように介在物から上下にミスフィット転位が拡張しているためで，構造は hcp のままである.

(c) の Mg/Zr は界面に高い応力が局在しているが，介在物中は Mg 母相と応力差がない．Zr は hcp 構造で最近接原子間距離も Mg と同じためである．(d) は Mg よりもPb の方が最近接原子間距離が大きいため，Mg の格子長さに設置した Pb 介在物が膨張して中はほぼ均一な値で圧縮状態になっている．ここでもミスフィット転位が界面上下に拡張している．

43

Table 4.1
element
structurenearest neighbordistance [nm]Young’s modulus [GPa]1st eigenvalueη α(1) [GPa]Properties in bulk perfect lattice.
Mg
Al
Zr
hcp
fcc
hcp

Pb
fcc

0.32

0.29

0.32

0.35

53.5
13.9
(C44 )

63.2
46.4(C11 − C12 )

121.7
25
(C44 )

11.1
7.6(C11 − C12 )12 GPa12 GPa-3 GPa(a) No obstacle,= 0.023-3 GPa(b) Mg/Al,= 0.023

15 GPa12 GPa-3 GPa(c) Mg/Zr,
Fig. 4.6-4 GPa(d) Mg/Pb,= 0.024= 0.028αaround the right crack tip and obstacle.
Distribution of σyy44α
図 4.7 は図 4.6 の応力分布の，y=0 のき裂面上の原子応力分布 σyyをプロットしたものである．図中の赤線は Mg の表面エネルギー γs とヤング率 E，ポアソン比 ν か√ら算出される Griﬃth き裂の破壊靭性値 KIC = 2Eγs /(1 − ν 2 )=24 MPa・m1/2 でプ√ロットした線形破壊力学での応力分布 σyy =KI / 2πr である．Fe と異なり，介在物なしの Mg ではき裂進展時の応力分布は KIC におおよそ一致した．介在物を含む系もほぼ一致しているが Al を含む (b) のモデルではき裂-介在物間の応力は (a) のそれよりも低い．ただしき裂先端の応力は (a) のそれとほぼ同じである．(c),(d) は介在物の界面の応力に変化があるだけで，前章の Fe 中の介在物に比べて母相と介在物の応力差は大きくない．

45

Stress,

Distance from cracktip, r, nm(b) Mg / Al

Fig. 4.7

Stress distribution on y=0 crack plane at sress peak point.

46

4.4

き裂進展時の AES 固有値

図 4.8 は図 4.6 と同じスナップショットを各原子の η α(1) の値で着色して示したものである．図の上部にはシミュレーションセル中の全ての Mg 原子と介在物原子でそれぞれ平均した η α(1) の値を示している．図では η α(1) が小さい値，すなわち不安定側の原子を暖色で示している．Al と Pb は介在物の上下から真横の basal 面に η α(1)の値が下がっている部分が見られる．これは緩和計算時に発生したもので Mg 母相と Al，Pb 介在物の格子ミスフィットのため，hcp 構造の主すべり系である basal 面に拡張した転位である．(c) の Zr 介在物ではミスフィット転位は生じていない．Feの時と比べ Mg は完全結晶での η α(1) の値が低く，引張による低下は小さい．一方，

Al は完全結晶の値 46.4 GPa に対して介在物の状態では引張前から-0.4 GPa と平均で不安定となっている．しかし，Al 内の分布を見ると η α(1) が著しく低い「赤」の原子が界面に点在し，内部は後の Pb に比べて高い．このため Al 介在物は初期に変形するものの，双晶境界と界面の会合部に沿った割れの方が支配的になってき裂は介在物を回避したものと考えられる．Mg 母相より η α(1) が高い (c) の Zr は母相より硬く，変形することなくき裂は回避した．バルクでは Al よりもヤング率が低い Pb であるが，介在物では平均の η α(1) は 1.0 GPa と Al よりは高い．しかし介在物内部が一様に大きな負の値を示しているのは Pb で，前章の Fe 中の Zr のようにキャップのようにき裂先端を覆ってき裂進展を阻害した．

47

Average

Mg : 11.1 GPa
Al : -0.4 GPa

AverageMg : 11.2 GPa

16 GPa16 GPa-10 GPa-10 GPa(a) No obstacle,

Average(b) Mg/Al,= 0.023

Mg : 11.1 GPa
Zr : 27.1 GPa

Average= 0.023

Mg : 10.8 GPa
Pb : 1.0 GPa
16 GPa30 GPa-10 GPa-10 GPa(c) Mg/Zr,

Fig. 4.8(d) Mg/Pb,= 0.024= 0.028

Distribution of η α(1) around the right obstacle.

48

4.5

き裂−介在物間の不安定モード

図 4.9 はき裂先端が Al 介在物に到達して双晶変形が発生する瞬間と，介在物を回避する瞬間のスナップショットである．原子種と η α(1) の正負で色分けしたもの((a),(d))，η α(1) の値で着色したもの ((b),(e))，η α(1) < 0 の原子の η α(1) の固有ベクト
ル {∆εxx ,・・・,∆γxy }T から求めた最大せん断方向を xy 面に投影した図 ((c),(f)) を示している．矢印は (a) と (d) の色を薄めたスナップショットの上に重ねて表示した．

(a) でき裂先端から下向きの pyramidal 面に双晶変形の元となる転位が発生し，そこには η α(1) < 0 の赤色の原子群が確認できる．(b) の η α(1) 分布を見ると，先述の拡張したミスフィット転位の低い η α(1) を示し，また介在物上下の Mg の η α(1) の値が下がっている．介在物の中ではき裂先端に近い部分は濃い青の η α(1) > 0 の原子が多い(図 (a))．不安定モードの図 (c) はやはり介在物中のモードは様々な方向を向いているが，界面部に大きな矢印が見られる．εyy =0.045 では双晶変形が広がって先端形状が台形となり，き裂先端と双晶境界と界面の交点から割れを生じ始める (図 (d))．
この交点部分には η α(1) < 0 の赤，水色原子が多い．この「新しい」き裂先端は，この後 η α(1) が低い方向へ流れるように介在物界面 → 双晶境界部へと進んだ．図 (f) では割れを生じる部分の η α(1) < 0 の介在物原子が 10 時方向に並んだ矢印を示しており，割れに伴う集団的な移動方向を反映している可能性がある．図 4.10 に Zr 介在物に達したき裂先端から割れが発生し，母相-介在物の界面ですべりを生じてボイドが発生した時のスナップショットを示す．き裂先端 → 双晶境界の割れは先の Al の図

4.9 と同じである．Zr 介在物には η α(1) < 0 の原子がないため，割れ前方の η α(1) < 0
原子のモードが明白である．図 (b) より，双晶境界の η α(1) は大幅に下がっていて局49

所変形が生じやすいことが分かる．εyy =0.047 までに双晶境界に沿ったすべりが生じて，介在物下部にボイドが形成されている．すべりの起点 (き裂先端) の η α(1) < 0
原子 (図 (d) 赤原子) に大きな矢印が確認できる．図では 11∼12 時の方向を向いているが最大せん断の ± は等価であるので，逆向き方向に原子が動いてすべりを生じると考えられる．図 4.11 は Pb 介在物にき裂先端が進入した時とその後介在物の中でき裂が開口した時のスナップショットである．図 (a) を見ると Pb 介在物は外側の原子が η α(1) > 0 で内側が軟化している．(b) からも介在物の中心ほど暖色となっており η α(1) が下がっている．図 (c) の介在物の η α(1) < 0 の変形モードは放射状になっていて何らかの特徴を見出すことは難しい．なお，この時点で介在物は格子長さの違いのため膨張しているが内部の構造は hcp 構造のままであることを確認している．

εyy =0.038 になると介在物の内側の η α(1) が一部回復している．この部分は構造が変化しており fcc 構造もしくはそれに近い構造になっている．図 (f) は (c) と同じく介在物中のモードは複雑であるが，全体的にその大きさは (c) よりも大きくなっている．また，図 (d) 右下の赤い原子群は図 4.5 で示した線状欠陥であり，図 (f) で線状
欠陥の端部には 2∼3 時方向の大きな変形モードが認められる．

50

Mg,
Mg,
Al,
Al,

16 GPa-10 GPa= 0.032(a) atom type and sign of(b)(c) eigenvector(d) atom type and sign of(e)(f) eigenvector= 0.045

Fig. 4.9

Snapshots of crack tip just reaching and penetrating Al obstacle. Left:colored by atom type and the sign of η α(1) . Middle: colored by η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

51

Mg,
Mg,
Zr,
Zr,

-5 GPa32 GPa= 0.043(a) atom type and sign of(b)(c) eigenvector(d) atom type and sign of(e)(f) eigenvector= 0.047

Fig. 4.10

Snapshots of crack tip just reaching and penetrating Zr obstacle. Left:colored by atom type and the sign of η α(1) . Middle: colored by η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

52

Mg,
Mg,
Pb,
Pb,

10 GPa-15 GPa= 0.032(a) atom type and sign of(b)(c) eigenvector(d) atom type and sign of(e)(f) eigenvector= 0.038

Fig. 4.11

Snapshots of crack tip just reaching and penetrating Pb obstacle. Left:colored by atom type and the sign of η α(1) . Middle: colored by η α(1) .
Right: maximum shear direction of η α(1) < 0 atoms calculated from theeigenvector of η α(1) .

53

5

結言
単結晶の薄板状周期セル中のモード I き裂について，進展方向の原子を他原子に置換して円柱状の原子クラスター (介在物) を配置した系の分子動力学シミュレーショ
ンを行うとともに，き裂進展挙動を原子弾性剛性係数 B αij =∆σiα /∆εj の第一固有値

η α(1) ，ならびに固有ベクトルの視点から議論した.
3 章では母相を Fe とし，介在物を 4 種の元素 (Ni,Co,Zr,Mo) として解析を行った．
その結果を以下に示す．

1. き裂は Ni,Co 介在物を貫通した．Zr 介在物は Fe との格子長さの違いからアモルファスに近い状態で存在し，き裂先端を包み込むようにして大きく鈍化して進展を著しく阻害した．Mo 介在物は接近したき裂を数原子離れた位置に停滞させた後，Fe/Mo 界面ですべりを生じて介在物上部にボイドを発生させ，これがき裂先端と合体してき裂は介在物を回避して最終破断した．

2. き裂が進入した Ni，Co，Zr 介在物のいずれも周囲の Fe より第一固有値 η α(1)が小さく，Fe 母相に比べて「柔らかい」状態であった．また，回避する挙動を示した Mo 介在物は η α(1) が Fe 母相より高い値を示し「硬い」状態であった．

3. き裂が介在物に達する前後について，各原子の η α(1) の分布，η α(1) < 0 原子の位置，そして η α(1) < 0 原子の固有ベクトル {∆εxx ,・・・,∆γxy }T から算出した
最大せん断方向 (不安定モード) を可視化した．

4. Ni,Co 介在物では介在物内にき裂が進入すると介在物原子の η α(1) が回復し安定化した．また，き裂開口やすべり変形部に η α(1) < 0 原子の対応するモード

54

が確認された．

5. アモルファス化した Mo 介在物は Fe 原子と境界で互いに混ざりあっており，介在物内部に進入して鈍化したき裂表面には Fe 原子が存在していた．4 のような介在物中の η α(1) の回復は見られなかった．

6. Mo 介在物では，先端が 120 °に鈍化して停滞していたき裂とつながろうとするような不安定モードが介在物界面原子に確認できた．また，鈍化したき裂先
端から新たに発生したき裂と介在物表面に沿うすべり面の三交点に，ボイドを発生させたすべり変形の不安定モードが認められた．

4 章では Mg を母相とし，3 種の介在物 (Al,Pb,Zr) について検討した．その結果を以下に示す．

1. き裂は Al,Zr 介在物を回避し，Pb 介在物は貫通した．Al と Zr は介在物とき裂
先端から pyramidal 面に双晶変形の帯が形成され，き裂は双晶の方向に回避した．Pb 介在物ではき裂が介在物に達する前に介在物から線状欠陥の発生が見られた．

2. Al と Pb 介在物は引張前にミスフィット転位を生じ，hcp 構造の主すべり系である basal 面に拡張した．

3. き裂が回避した Zr 介在物は周囲の Mg より η α(1) が大きく，き裂が進入した Pb
は周囲の Mg より η α(1) が小さいのは前章と同じである．しかしながら，Al 介
在物は Mg より η α(1) が低かったが，双晶境界と介在物界面の会合部に沿った割れが支配的になってき裂は回避挙動を示した．

55

4. き裂が回避した Al,Zr 介在物はいずれもき裂先端と介在物の交点から pyramidal
面に沿った双晶境界に η α(1) が低い領域が存在し，そこを起点として回避挙動を生じていることを確認した．

5. 介在物中の η α(1) < 0 原子のモードは複雑かつ大きな負の固有値を持つため特徴的な挙動を見出すことは難しいが，介在物が η α(1) > 0 の Zr の場合は 4 の挙動における η α(1) < 0 のモードを明らかにすることができた．

56

参考文献
(1) F. F. Abraham, et al., Physical Review Letters, Vol.73, No.2, pp.272-276,(1994).
(2) K. Nishimura, et al., Computer Modeling in Engineering and Sciences Vol.2,
No.2, pp.143-154, (2001).
(3) S. J. Noronha and D. Farks, Materials Science and Engineering, Vol.365, Issue
1-2, pp.156-165, (2004).
(4) H. Rafii-Taber, et al., Mechanics of Materials., Vol.38, pp.243-252, (2006).
(5) Z. Yang, et al., Computational material science, Vol.112, pp.17-25, (2014).
(6) Z. Yan, et al., Computational Materials Science, Vol.50, pp. 1754-1762, (2011).

(7) J. Yu, et al., RSC Advances, Vol.4, pp.32749-32754, (2014).
(8) H. Gi, et al., Computational Materials Science, Vol.183, 109806, (2020).
(9) J. Yin, et al., International Journal of Pressure Vessels and Piping, Vol.194,
104519, (2021).
(10) L. Liming, et al., The International Journal of Advanced Manufacturing Technolog, Vol.69, pp.2131-2137, (2013).

57

(11) D. Wenjiang, et al., Compasites Part B: Engineering, Vol.68, pp.8-17, (2015).

(12) W. Liang, et al., Engineering Fracture Mechanics, Vol.276, 108898, (2022).
(13) K. Yashiro, Computational material science, Vol.112, pp.120-127, (2016).
(14) K. Yashiro, et al., The Japan Society of Mechanical Engineers, Vol.81, No.829,(2015).
(15) K. Yashiro, Computational material science, Vol.131, pp.220-229, (2017).
(16) K. Yashiro, Computational material science, Vol.147, pp.72-80, (2018).
(17) K. Yashiro, Philosophical Transactions of the Royal Society A, 379, 20200124(2021).
(18) K. Yashiro, et al., Journal of the Society of Material Science, Japan, Vol.71,
No.2, pp.127-134, (2022).
(19) 上田顯, コンピュータシミュレーション, 朝倉書店（1990）．
(20) 洲之内治男, サイエンスライブラリ–理工系の数学=15, 数値計算, サイエンス社,(1978)．

(21) 香山正憲, 固体材料の電子状態の基礎, 材料学会勉強会資料, (1993)．
(22) X. W. Zhou, et al., Physical Review B, Vol.69, 144113, (2004).
(23) X. W. Zhou, et al., Physical Review B, Vol.69, 035402, (2004).
58

(24) X. W. Zhou, et al., Acta mater. Vol.49, No.19, pp.4005-4015, (2001).
(25) S. Nosé, The Journal of Chemical Physics, Vol.81, No.1, pp.511-519, (1984)．
(26) D. C. Wallace, Thermodynamics of Crystals, Wiley, Newyork, (1972)．
(27) J. Wang, et al., Physical Review Letters, Vol.71, No.25, pp.4182-4185, (1993)．

(28) J. Wang, et al., Physical Review Letters, Vol.52, No.17, pp.12627-12635,(1995)．

59

謝辞
本研究を遂行するにあたり，屋代如月教授には浅学非才な著者に対し懇切丁寧に指導していただきました．ここに心より御礼申し上げます．本論文を完成させるにあたり，広い視野から研究全般に対して多くのご助言をいただきました内藤圭史助
教にも心より感謝いたします．ともに切磋琢磨し合った加藤典子氏を始めとする研究室メンバー，面倒を見てくださった棚橋直哉氏を始めとする先輩方にも御礼申し上げます．最後に，編入学してから 4 年間の学生生活を暖かく見守り精神的にも経済的にも支えて頂いた家族に心より感謝いたします．ありがとうございました．

60