卒 業 論 文(題 目)マグネシウム中のモード I き裂先端の局所格子不安定性解析

平成 27 年度

岐阜大学工学部
機械システム工学科

氏名 吉原大志目次
1

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

12

解析手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.1

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

2

2.2

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

3

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

3

2.2.1
2.3

高速化手法

2.4

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

72.5

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

8

2.5.1

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

8

2.5.2

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

8

2.5.3

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

2.5.4

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

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

4

3

シミュレーション方法. . . . . . . . . . . . . . . . . . . . . . . . . .

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

13

5

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

22

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

23

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

24

i

1

緒言
マグネシウムは密度が鉄の約 4 分の 1, アルミニウムの約 3 分の 2 と実用金属の中で最も軽い．また比強度も高く，地球上に豊富に存在し，リサイクル性を持つことから優れた金属材料であると考えられている．そのためマグネシウムは人体へのイ
ンプラント (1) や輸送機器 (2) などの様々な分野での利用法が考えられてきた．特に輸送機器の分野においては，金属資源の枯渇や地球温暖化問題が深刻化しており，輸送機器の軽量化・金属部品のリサイクル性の観点から，マグネシウムが注目を集めている．少量の遷移金属 (Zn) と希土類金属 (Y) を付加することで LPSO(Long-Period

Stacking Ordered) 構造 (3) と呼ばれる積層構造を持った新しいマグネシウム合金も開発され，その構造安定性や変形メカニズムが盛んに研究されている．したがってマグネシウムの特性を把握するために様々な検討がされてきた．原子シミュレーショ
ンを用いた研究では，単結晶中のボイド成長 (4) を検討した例があり，ボイド付近の塑性変形に影響されてボイド成長が試料のサイズと結晶軸方向に強く依存することが明らかになっている．変形双晶のメカニズムの検討 (5) ではすべり面と双晶の関係が報告されている. くり返し変形下のき裂進展挙動 (6) を検討した研究では，き裂伝
播は結晶軸方向に依存することを報告している．また，LPSO 構造の変形挙動 (7) については，fcc 構造層が双晶変形を抑制するために底面すべりを生じることなどが明らかにされている．
我々の研究グループでは局所変形抵抗を表す物理量である原子弾性剛性係数 (Atomic

Elastic Stiﬀness, AES) の正値性に着目し，それに基づいて局所変形の開始を統一的に議論する研究を展開してきた (8) ．本論文では Liu らの EAM ポテンシャル (9) を用
いて hcp マグネシウムのモード型き裂の分子動力学シミュレーションを行うとともに，AES の正値性の観点からき裂挙動を議論する．

1

2

解析手法

2.1

分子動力学法

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

(2.1)

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

Fα = −

∂Etot
∂r α

(2.2)

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

α

(2.3)
(2.4)

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

dr α= v α (t)
dt

(2.5)

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

α

α

2

(2.6)
(2.7)

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

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

(

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

)

(2.8)
(2.9)

が得られる．∆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.10)(2.11)と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．第一ステップの計算 (t = 0) では，t = ∆t での座標 r α (∆t) は式 (2.6) と初速度から得ることができる．

2.2

原子間ポテンシャル

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

2.2.1

EAM ポテンシャル

EAM は金属中の多体効果を良好に再現することから広く用いられている．密度汎関数理論に基づき，まず金属材料における系のポテンシャルエネルギー Etot は原

3

子を価電子雲中に埋め込むエネルギーと原子間の 2 体間相互作用の和で与えられるとする．さらに，埋め込みエネルギーは埋め込む位置の電子密度にのみ依存すると仮定することによって，系全体のエネルギーは次式のように表わされる．

Etot =N∑F (ρα ) +

α

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

ここで，ρα は原子 α の位置における多体効果を考慮する密度を表し，F (ρα ) は密度

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

(2.13)

β(̸=α)

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

2.3

高速化手法

原子数 N の系において粒子間の全相互作用を評価すると，1step に N × (N − 1)回の計算が必要となり，N が大きくなると極めて膨大な計算量となる．実際には，
一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ
半径 rc ) 内の粒子からの寄与を効率よく計算することにより高速化できる．従来よく用いられてきた高速化手法に粒子登録法がある．これは，図 2.1 に示したように，

rc よりひとまわり大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) 計算負荷が減少される．

4

しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計算負荷により打ち消される．

Fig. 2.1

Schematic of bookkeeping method.

5別の高速化手法としてブロック分割法がある．図 2.2 に示すように，シミュレー
トする系をカットオフ距離程度の格子状に分割し，各ブロックに属する粒子をメモ
リーに記憶する．着目している粒子に作用する力を評価する際には，図 2.2 に示すよ
うに，その粒子が属するブロックおよび隣接するブロックから相互作用する粒子を探索して行う．粒子が属するブロックは，粒子の位置座標をブロックの辺長 bx，byで除した際の整数により判断できるので，ブロック登録時の計算負荷は粒子数 N のオーダーとなる．したがって，粒子登録法では登録更新の負荷が大きくなるような大規模な系でも高速化が可能である．本研究はブロック分割による高速化を採用している．

Fig. 2.2

Schematic of domain decomposition method.

6

2.4

速度スケーリング法

分子動力学解析における温度制御には一般的には速度スケーリング法が用いられる．以下にその説明を記す．温度 T は，統計熱力学より導かれる次式で表される．

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

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

viα0 =

3kB T0mα)0.5(2.15)

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

viα =

3kB T
mα

)0.5

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

viα0
T0
=
α
vi
T

)0.5

(2.17)

系の温度を T から T0 にするには，式 (2.17) の右辺を現在の速度に掛けてやればよい．式 (2.10) を変形させると Verlet 法における式は

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

(2.18)F α (t)∆tmα(2.19)v α (t + ∆t) = v α (t) +とみなせるので，√ααr (t + ∆t) = r (t) +{T0F α (t) 2r α (t) − r α (t − ∆t) +
∆t
T
mα
7

}

(2.20)

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

2.5

原子弾性剛性係数2.5.1

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

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

σij

1
=
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

∂ 2F
∂ηij ∂ηkl)(2.21)と定義される．ここで，F は Helmholtz の自由エネルギー (断熱過程では内部エネル
ギー U )，V は結晶の体積，ηij は平衡状態 (無負荷とは限らない) からの仮想的な微小ひずみである．一方，無負荷平衡状態を基準とするひずみ εij と応力 σij の関係は，

2 つの平衡状態間の変形を考えて導出される次の弾性剛性係数によって表される．
Bijkl ≡

(

∆σij
∆εkl

)

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

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

2.5.2

原子応力

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

8ルエネルギーの単位体積当たりの変化として与えられる．

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

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

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

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

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

∂rαβ
∂ηij

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

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

σij

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

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

9

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

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

σij

N
1 ∑
σα
=
N α ij(2.30)

となる．

2.5.3

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

1 ∂ 2 Etot
V ∂ηij ∂ηkl

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

Cijkl

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

′′

α

{ ∑
N

′

αβ

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

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

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

α

{ ∑
N

′

αβ

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

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

2.5.4

原子弾性剛性係数以上で定義した原子応力，弾性係数から，原子弾性剛性係数は以下で評価できる．
α
α
α
α
α
Bijkl
= Cijkl+ (σilα δjk + σjl
δik + σik
δjl + σjk
δil − 2σijα δkl )/2(2.35)α sym
ααWang らによる提案に従い，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.35) の対応を表 2.1 に示す．

Table 2.1
C11 + σ1

2
C12 − σ1 +σ2
C22 + σ2

6 × 6 matrix of Bij

3
C13 − σ1 +σ2
3
C23 − σ2 +σ2
C33 + σ3

C14 − σ24
C24 + σ24
C34 + σ243
C44 + σ2 +σ
2

sym.

11

C15 + σ25
C25 − σ25
C35 + σ25
C45 + σ26
1
C55 + σ3 +σ2

C16 + σ26
C26 + σ26
C36 − σ26
C46 + σ25
C56 + σ242
C66 + σ1 +σ2

3

シミュレーション方法

図 3.1 に示すように，平板状の周期セル中央に貫通き裂を配置した系を対象とす√
る．系の大きさは 132a × 82 3a × 4.869a (a はマグネシウムの格子長さ 3.206 nm，

c/a = 1.623)，原子数は 15 万である．き裂幅は 0.3lx (lx は方向セル長さ)，き裂先端の曲率半径はグリフィスき裂の条件を満たすようにしている．温度 0.1K で垂直応力が 0 となるように 10000fs の緩和計算を行った後，き裂垂直方向に引張ひずみを毎ステップ 1.0 × 10−6 増加させて引張シミュレーションを行った．き裂貫通方向，進展方向のセル寸法は固定している．応力‐ひずみ曲線を求めるとともに， 1000fs 毎に記録した原子配置データから各原子位置における Bijα および Bijα εj = η α εi の 6 つ

Lyの固有値 η α(k) を求めた．

2c=0.3Lx

y[0001]Lx

Lz

x[1000]z[0110]Fig. 3.1

Simulation cell

12

4

シミュレーション結果および考察
引張シミュレーションで得られた応力‐ひずみ線図を図 4.1 に示す．また，ひずみ εyy = 0.08 までの変化を拡大したものを図 4.2 に示す．ひずみ εyy = 0.04 近傍まではほぼ線形に上昇しているが，その後応力上昇が鈍化し，ひずみ εyy = 0.22 近傍
で最大応力 σ = 3.5GPa を示した後急減した．図 4.2 をみると，ひずみ εyy = 0.14
近傍でごくわずかではあるが下向きに凸となっており，構造変化を生じていることが示唆される．ひずみ εyy = 0.6 までの変形の様子を図 4.3, 図 4.4 に示す．変形は(a) → (b) → (c) → (d) → (e) → (f ) → (g) → (h) の順に進んでいく．図 4.3 の (b)のき裂形状の変化から塑性変形が始まっていることが確認できる．すなわち図 4.2 の鈍化点 (応力‐ひずみに非線形性が現れる点 ) が塑性変形の開始点である．(c) → (d)ではき裂形状が大きく変わっている．応力‐ひずみ曲線が下に凸となった点が (d) に対応しており，後述のように相変態な構造変化が周期き裂をつなぐように (X 状に)生じたものと考える．応力‐ひずみのピーク点 (図 (e)) では相変態の交差部からボイドが発生し，それが成長して破断した．

13

Fig. 4.1

Fig. 4.2

Stress-strain curve

Magnification of Fig.4.1

14

Fig. 4.3(a) yy=0.013(b)yy＝0.04(c)yy＝0.051(d)yy＝0.14

Snapshots of atoms colored by potential energy (εyy = 0.013〜0.107)

15

(e)yy＝0.22(f)yy＝0.35(g)yy＝0.47(h)yy＝0.6

Fig. 4.4

Snapshots of atoms colored by potential energy (εyy = 0.22〜0.6)

16

detBijα および Bijα の固有値 η α(1) ，η α(2) が負になった原子数の変化を図 4.5 に示す．
また，ひずみ εyy = 0.08 までの変化を拡大して図 4.6 に示した．ひずみ εyy = 0.04 近
傍から detBijα ，η α(1) が負の原子が急増し，ひずみ εyy = 0.051 で最大となっている．
後述するようにこの不安定原子の急増は相変態の臨界点に対応するものと考えられ
る．η α(2) が負の原子はひずみ εyy = 0.05 以降に見られるが，これも相変態による局所構造変化時に生じたものと考える．η α(3) 以上は負になった原子はほとんどみられなかった．detBijα = η α(1) η α(2) η α(3) η α(4) η α(5) η α(6) であり，η α(3) 以上はほとんどが常に正であるので detBijα < 0 と η α(1) < 0 の変化に差が生じるのは η α(1) , η α(2) がともに負のときだけである．

Fig. 4.5

Number of atoms with negative determinant/eigenvalue

17

Fig. 4.6

Magnification of Fig.4.3

18

図 4.7 は，η α(1) < 0 となる原子が多く出現し始めたひずみ εyy = 0.013 におけるき裂近傍のスナップショットを，η α(1) < 0 の原子を赤く着色して示したものである．

η α(1) < 0 の原子はき裂先端にのみ存在している．図 4.8 に，応力上昇が鈍化しはじめるひずみ εyy = 0.04 におけるセルの全体のスナップショットを，同様に η α(1) < 0の原子を赤く着色して示す．右のき裂先端から右下方向に離れた位置に η α(1)<0 の原子が確認でき，転位射出によって応力鈍化したものと結論付けられる．図 4.9 にη α(1) < 0 の原子数が最大となったひずみ εyy = 0.051 におけるスナップショットを示す．η α(1) < 0 の原子の分布はすべり方向が定まっておらず，相変態もしくはアモ
ルファス化のような様相を呈している．図 4.10 に応力‐ひずみのピーク点における

detBijα < 0 の原子の分布を示す．detBijα < 0 の原子は転位や積層欠陥などの欠陥部分に見られるが，多結晶体のような分布となっており，周期き裂の中間点にすでにボイドを生じていたことがわかる．

Fig. 4.7

Distribution of η α(1) < 0 atoms (red circles)， εyy = 0.013

19

Fig. 4.8

Distribution of η α(1) < 0 atoms (red circles)，εyy = 0.04

Fig. 4.9

Distribution of η α(1) < 0 atoms (red circles)，εyy = 0.051

20

Fig. 4.10

Distribution of

detBijα < 0 atoms

21(red circles)，εyy = 0.22

5結言
原子弾性剛性係数 (Bijα ) の正値性に基づいてマグネシウム中のモード I 貫通き裂の進展挙動を議論するため，薄板状周期セル中央に貫通き裂を有する単一き裂モデルの引張シミュレーションを行い，Bijα の負の固有値等について検討した．得られた結果をまとめて以下に記す．

1. き裂先端からの転位射出の後，周期き裂間を連結するように相変態的な変形を生じてき裂形状が変化した．

2. 応力‐ひずみのピーク点では， 1 の相変態交差部からボイドが発生して，それが成長・合体して破断した．

3. 応力‐ひずみは， 1 の転位射出点で非線形性が現れており，1 の相変態点でわずかに下に凸の曲線を描いている．

4. Bijα の負の第一固有値 η α(1) を持つ原子は，1 の相変態の前に急増している．
5. 転位射出や相変態の考察は，η α(1) < 0 の原子の分布を確認したことで導かれた．

22

参考文献
(1) E. Zhang，A volume in Woodhead Publishing Series in Biomaterials，pp. 2357，(2015)(2) B.L.Mordike，et al，Materials Science and Engineering，A302，pp. 3745，(2001)(3) M.Tane，Acta Materialia，61，pp.6338-6351，(2013)(4) T. Tang，et al.，Acta Materialia. 58，4742-4759，(2010)(5) R. Aghababaei，et al.，Acta Materialia，69，pp. 326-342，(2014)(6) T. Tang，et al.，Computational Materials Science 48，pp. 426-439，(2010)(7) R.Matsumoto，et al，Materials Transactions，Vol. 54，pp. 686-692，(2013)(8) 屋代如月 他，日本機械学会論文集，81(829)，15-00271，(2015)(9) X.Y.Liu，et al，Modelling and Simulation in Materials Science and Engineering，
pp.293–303，(1996)23謝辞
本研究を行うにあたり，指導教官の屋代如月教授から丁寧な指導を頂きました．著
者にとって新たな分野であるコンピューターシミュレーションの基礎から研究者としてのあり方を示していただき，著者の成長を促して頂きました．ここに感謝の意を表します．また同研究室の内藤圭史助教には，研究環境の準備や他愛もないことで気さくに話しかけていただき，充実した生活を送ることが出来ました．また，精神面でも支えていただき本当にありがとうございました．
そして喜怒哀楽を共有し，切磋琢磨し合い共に過ごした，堤貴文，寺田稜，西川涼一郎，堀広志，山田泰成の 5 名にも同様に感謝の意を表します．
最後に 4 年間の学生生活を送るにあたり，遠方から経済面・精神面で支えて頂いた両親に感謝いたします．

24