修士論文

Phase Field 法・NEB 法によるL12 構造中の超転位ダイナミクスの解明

指導教員：屋代 如月

蟹川淳

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

Study on Dynamics of Super Dislocationin L12 structure based on
Phase Field Simulation and NEB Calculation

February 2011
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan
Jun KANIKAWA要約
本研究は，Ni 基単結晶超合金中の超転位をモデル化することを目的として，分子静/動力学による逆位相境界 (APB) エネルギーの評価，Phase Field 法によるシミュレー
ション，NEB 法による APB 形成および γ/γ 0 界面を転位が貫通するときのエネルギー
バリアの評価等を行った．まず，熱膨張を想定し，平均格子間隔を広げて静力学的に解析した時の APB エネルギーと，実際に熱揺動の効果を考慮したときのそれの違いを議論するために，EAM ポテンシャルによる静力学解析ならびに動力学解析を行った．熱膨張を想定して体積 (格子長さ) を増加させると，(111) 面での APB エネルギーは 0K の格子長さにおける値 180.4mJ/m2 からなだらかに単調減少した．一方，10〜

1000K の温度条件下で分子動力学解析を行ったところ，APB エネルギーはいずれも145.6±2.2mJ/m2 となり，静力学解析では APB エネルギーを過大に評価してしまうことがわかった．次に，γ 0 相の逆位相境界（APB）のエネルギーを考慮した超転位のモデル化を行った．γ 0 相 (L12 構造) 単相中の 2 本の転位のシミュレーションでは，与えた APB エネルギーの値に応じて，転位論から導かれる超転位の幅を再現できること，

γ 相中に 2 本の転位を配置し，異なるせん断応力を与えて前方の γ 0 相をカッティング
させるシミュレーションでは，τ =300MPa では転位は γ 0 相に侵入できず界面で停滞するのに対し，τ =400MPa では超転位として γ 0 相に侵入し一定の幅 (約 30b) を保ち運動すること，などを示した．最後に，NEB 法による APB 形成および γ/γ 0 界面を転位が貫通するときのエネルギーバリアの評価を行った．APB を形成する途中の最小エネルギー経路での各位相点におけるすべり面上下の原子配置をみたところ，fcc におけるショックレーの部分転位と同じ原子移動であること，そのときの最初と後続の部分転位に相当するエネルギーバリアはそれぞれ 371.8mJ/m2 ，479.7mJ/m2 であること，Ni 相と Ni3 Al 相の (010) 整合界面を貫通するときの最小エネルギー経路は，単一ピークをもつエネルギーバリアとなり，その値は 0.80 × 10−3 eV/Å であること，などがわかった．

2
The main objective of this study is to propose a computational model for super
dislocation in Ni-based superalloys. We have conducted molecular statics/dynamics
simulations for evaluation of the anti-phase boundary (APB) energy, phase field simulations for dislocation motions in γ 0 phase and at γ/γ 0 interface, and NEB analysis forenergy barriers of APB formation and dislocation cutting through the γ/γ 0 interface.
First, in order to discuss whether we can evaluate the APB energy at high temperatureby static calculation only with changing the lattice spacing or not, we have compared
the molecular statics and dynamics calculation using the same EAM potential. In thestatic calculation, the APB energy has decreased monotonically from 180.4mJ/m2 at
0K, if we increase the lattice length mimicking the thermal expansion. In contrast, theAPB energy is always evaluated as 145.6±2.2mJ/m2 by molecular dynamics simulation, despite of the diﬀerent temperature of 10K, 100K, 300K and 1000K. We have thenproposed a phase field model of super dislocation in the γ 0 phase (L12 structure) byadding the APB energy term as similar as the Peierls energy term for Shockley Partials
in fcc. We have first demonstrated that this model can correctly reproduce the widthof a super dislocation in γ 0 phase according to the given APB energy. Then we haveapplied this PF model to more complicated situation, i.e. dislocation cutting throughthe γ/γ 0 interface. Dislocations cannot penetrate into the γ 0 phase, if we set two edge
dislocations in γ phase and proceed to the interface with shear stress of 300MPa. On
the other hand, if we increase the shear stress to 400MPa, the two dislocations cut into
the γ 0 phase and form a super dislocation with constant width of 30b (Burgers vector).
Finally, the energy barriers against APB formation and dislocation cutting throughthe γ/γ 0 interface are evaluated by NEB analysis. We have obtained the minimum
energy path for APB formation and shown that the atom migration is same as the
Shockley partials in fcc and their energy barriers are 371.8mJ/m2 and 479.7mJ/m2 forleading and trailing partials, respectively. On the other hand, the minimum energy
path for an edge dislocation has a single peak to cut into the γ 0 phase, and it’s value
is evaluated as 0.80 × 10−3 eV/Å.

目次
第 1 章 緒論

1

第 2 章 解析手法

5

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

5

2.1.1

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

5

2.1.2

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

6

2.1.3

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

6

2.1.4

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

8

2.1.5

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

92.2 転位の Phase Field モデル . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.1

秩序変数と固有ひずみ . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.2

自由エネルギー汎関数 . . . . . . . . . . . . . . . . . . . . . . .

13

2.2.3

TDGL 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.2.4

TDGL 方程式の数値解法 . . . . . . . . . . . . . . . . . . . . . .

18

NEB(Nudged Elastic Band) 法 . . . . . . . . . . . . . . . . . . . . . . .

23

2.3

第 3 章 EAM ポテンシャルによる

APB エネルギー評価

25

3.1 逆位相境界・超転位 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

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

27

3.2.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.2.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.3 動力学解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.3.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.3.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.4 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

第 4 章 APB エネルギーを考慮した

Phase Field モデル
4.1

APB エネルギー項の導入

35. . . . . . . . . . . . . . . . . . . . . . . . .

i

35

目次

ii

4.2 γ 0 単相中の超転位のシミュレーション . . . . . . . . . . . . . . . . . .

37

4.2.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.2.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.3 γ 0 相をカッティングする転位のシミュレーション . . . . . . . . . . . .

43

4.3.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.3.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

44

4.4 複数の転位を導入した
γ 0 相カッティングシミュレーション . . . . . . . . . . . . . . . . . . . .

48

4.4.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

4.4.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

48

4.5 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

第 5 章 NEB 法による反応経路解析

5.1

APB 形成の反応経路解析

52. . . . . . . . . . . . . . . . . . . . . . . . .

52

5.1.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5.1.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.2 γ/γ 0 界面を貫通する転位の反応経路解析 . . . . . . . . . . . . . . . . .

56

5.2.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

5.2.2

解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . .

59

5.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

第 6 章 結論

62

参考文献

64

第 A 章 講演論文

68

謝辞

74

第1章
緒論
省エネルギー化と地球環境への負荷低減の観点から，ガスタービン燃焼ガス温度の高温化による発電効率の向上が求められている．セラミック材料はきわめて高温まで耐えらえるものの，ぜい性破壊するため実機での使用は甚大な被害をもたらす危険性がある．このため耐熱性は劣るが，やはり金属材料にかわるものはなく，耐熱性を向上させた超合金の開発が進められてきた．現在，最も有力な耐熱超合金として Ni 基単結晶超合金があげられる (1) ．

Ni 基単結晶超合金は，高温でのクリープ破壊の起点となる結晶粒界を排除すること
で耐熱性を向上させてきた．本材は，fcc 構造を持つ Ni を主成分とする γ 相の中に，

L12 型金属間化合物の Ni3 Al を主成分とする γ 0 相を 0.5µm 程度の立方体形状で格子状に析出させた特徴的な構造を有する（図 1.1）．また高温強度，耐腐食性，耐酸化性などを向上させるため Co，Cr，Mo，W など様々な合金元素が添加されており，第 3 世代と呼ばれる Ni 基超合金では Re を 5〜6mass ％添加し，高温クリープ強度と耐高温
腐食性の両方を効果的に向上させている (3) ．Ru など白金属元素を添加した第 4，5 世

Fig.1.1 Typical microstructure of nickel-base superalloys(2) ．

1

代 Ni 基単結晶超合金も開発されており (4) ，その組成改良が現在も盛んに行われている．これらの合金元素は，その添加割合を変化させることにより γ ，γ 0 相間の格子ミ
スフィットを制御することができる．γ/γ 0 界面における格子ミスフィットを持つ合金の方がより良好なクリープ特性・寿命を持つこと (5) ，Mo を添加し負の格子ミスフィッ
トを大きくすると後述の，ラフト化した γ 0 相表面に観察される界面転位網を緻密化する効果があること (6) など多数の実験的検討が行われている．この γ/γ 0 構造の Ni 基超合金の変形特性に与える影響については，γ 相の大きさや形状の違いによる降伏応力の変化 (7) や，強度に及ぼす格子ミスフィットや熱膨張係数の差の影響など (8),(9) が報告されている．
処女材におけるこの微細な析出形態は，クリープ変形を受けると γ 0 相が応力軸に対して平行または垂直方向に板状に粗大化するラフト化挙動を示す (10) ．ラフト化が完了した後はクリープ変形が加速し最終破断にいたる．ラフト現象には Al 原子の拡散が主たる因子と思われるが，ラフト化した γ/γ 0 構造中には多数の転位が確認されており，ラフト化におけるその役割が議論の対象となっている (11) ．

γ/γ 0 界面における転位挙動を明らかにするためには，転位芯すなわち原子配置まで考慮する必要がある．そのため，原子レベルでの変形・破壊現象を動的・連続的に観察可能な分子動力学法 (Molecular Dynamics ; MD) による研究も多くなされている．ミスフィット転位芯の構造 (12) ，圧子下で生じるプリズマティック転位ループとの相互作
用 (13) や刃状およびらせん転位との相互作用 (14) などの検討例があり，従来の転位論にはない新たな知見を与えている．また，原子構造を直接扱わず多数の転位の挙動を解
析する一階層上位スケールの解析手法として離散転位動力学法 (Discrete Dislocation

Dynamics ; DDD) がある．離散転位動力学法は，任意の曲線形状を有する転位線を多数の直線セグメントに離散化し，各セグメントと他の転位との相互作用を考えること
により，転位の運動を逐次追跡する手法である．DDD では，厳密には原子構造を考慮しなければならない転位芯レベルの相互作用 (対消滅，交差すべり，ジョグ・ジャンクション形成など) を臨界距離や臨界力などのローカルルールを導入することで再現可能にし，計算負荷を軽減している．また有限要素解析と連成解析することで，介在物および異種材料界面における転位挙動について検討した例 (15) や，境界要素法との連成解析で，ナノインデンテーションで観察されるナノスケールの塑性現象についての検討をした例 (16) もある．
しかしながら，原子を直接扱うことにより原子の熱揺動を表現できるまで時間スケー
ルを小さくとらねばならず，結果として MD では有限温度下で生じる確率的な現象（拡散や転位の上昇等）を扱うことはできない．一方，DDD は転位を線として扱うこと2

で，MD では扱うことが困難な多数の転位や熱活性化過程における現象を再現できるが，転位構造が複雑になると離散化したセグメントの多体問題が MD のそれ同様厳し
くなり，すぐ計算機能力の限界に達する．そこで，離散化した個々の物質座標の運動を追跡するのではなく，空間座標の「場の」時間発展を追う Phase Field 法が注目されている．Phase Field 法はもともと固液界面の問題（ステファン問題）において，界面を直接扱わずに，場の秩序変数の発展として扱うことで多大な成功を収めてきた．すな
わち，界面を直接離散化してその時間発展を追跡するのではなく，界面を含む全ての空間を分割し，空間に固定された各点の秩序変数の値の違いで界面を定義し，秩序変数の変化から界面の運動を追跡する．凝固の分野では，相変態，組織形成および形態形成など非常に多くのシミュレーションが行われている (17)∼(19) ．転位を結晶がすべった領域と未すべり領域の界面と考えれば，転位の発展の問題は同じように秩序変数で扱うことができる．このような考えから Phase Field 法を転位のシミュレーションに適用する試みも盛んになされてきた．Ortiz ら (20) によって 1 本の転位と障害物の切り合
いシミュレーションがなされた後，Ghoniem(21) ，Wang(22) のグループにより相次いで

3 次元の転位の Phase Field モデルが提案されている．完全転位を対象としたこれらの定式化は，Shen ら (23) によって一般化積層欠陥エネルギー (GSF) を考慮した Shockleyの部分転位へ拡張された．Ni 基超合金の転位に関しては，Hu と Chen(24) が，整合析出
構造の転位による拡散プロセスをモデル化するため，整合ひずみならびに転位のひずみ場を固有ひずみで表し，γ 0 相の成長への寄与を Phase Field 法により検討している．
著者らのグループでは，分子動力学法，転位動力学法および第一原理計算を援用して γ/γ 0 界面・構造中における転位挙動を検討してきた (25)∼(27) ． しかしながら，ラフ
ト化プロセスにおいては，転位と γ 0 相の相互作用だけでなく，Al 原子の拡散・γ 0 相の形態変化なども同時に扱わなければならない．元素の拡散と形態変化を同時に扱うことのできる最も有力な手法は先に述べた Phase Field 法であり，最終的には Al 原子の拡散 (γ 0 相の形態変化) および転位の両者を同時に扱った Phase Field シミュレーショ
ンの構築が可能であると考える．ここで，γ 0 相は Ni3 Al を主成分とする L12 構造であるため，γ 0 相中の転位挙動を扱うためにはまず逆位相境界 (APB) を考慮した超転位のモデル化を行わなければならない．転位動力学による検討では，これまでに第一原
理計算により APB エネルギーを求め，それをバックフォースモデルというローカル
ルールを適用することで γ/γ 0 構造中の複雑な転位挙動を再現している (26) ．しかしながら，第一原理計算では，原子の熱揺動を考慮できるような長時間の計算は不可能であるため，APB エネルギーはあくまでも「絶対零度」のものであった．そこで本研究では，APB エネルギーの温度依存性を評価するため分子静/動力学シミュレーション

3

による評価を行い，その結果を元に Phase Field シミュレーションにおける超転位のモデル化を行う．さらに，APB エネルギーを考慮するだけでなく，APB を形成する途
中のエネルギー障壁の評価が重要である，という視点に立ち，遷移状態理論に基づく

NEB(Nudged Elastic Band) 法による解析を行う．NEB(Nudged Elastic Band) 法は，反応における最小エネルギー経路から活性化エネルギーおよび活性化体積を求めることができ，金属界面からの転位の生成メカニズムの解明 (28) や転位と界面の相互作用(29)などに適用されている．以下に各章の概略を示す．

第 2 章では，本研究で用いた分子静/動力学，Phase Field 法，NEB 法など様々な計算力学手法の基礎について説明する．第 3 章では，熱膨張を想定し，平均格子間隔を広げて静力学的に解析した時の APB エネルギーと，実際に熱揺動の効果を考慮したときのそれの違いを議論するために，EAM ポテンシャルによる静力学解析ならびに動力学解析を行う．第 4 章では，これまでに提案されてきた転位の Phase Field モデ
ルを参考に γ 0 相の APB を考慮したシミュレーションモデルを提案する．第 5 章では，APB を形成する途中のエネルギー障壁，および，Ni/Ni3 Al の (010) 整合界面に停止した fcc の刃状転位が，Ni3 Al 相に侵入して leading partial となるときのエネルギー障壁
を NEB 法により評価する．最後に，第 6 章で本研究の総括を述べる．

4

第2章
解析手法
2.1

分子動力学法2.1.1

分子動力学法の概要

分子動力学法 (Molecular dynamics method; MD) は，系を構成する個々の原子についてニュートンの運動方程式

mid2 r i
= Fi
dt2

(2.1)

を作成し，これを数値積分することによって全原子の運動を追跡する手法である．ここで t は時間，r i ，mi はそれぞれ原子 i の位置ベクトル及び質量である．原子 i に作用する力 F i は系全体のポテンシャルエネルギー Etot の空間座標についての勾配ベクトルから次式のように求められる．

Fi = −

∂Etot
∂r i

(2.2)

式 (2.1) の数値積分には，Verlet の方法が簡便で高精度が得られるため MD 法ではよく用いられる．時刻 t±∆t での原子 i の座標 r i (t ± ∆t) を Taylor 展開すると

r i (t ± ∆t) = r i (t) ± ∆t

dr i (∆t)2 d2 r i (∆t)3 d3 r i+±+ O((∆t)4 )dt
2 dt2
3! dt3(2.3)となる．両式の和をとり式 (2.1) を代入すると

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

F i (t)+ O((∆t)4 )
m

(2.4)

を得る．これより，時刻 t − ∆t と t における全原子の位置を既知として，時刻 t + ∆t
における任意の原子 i の位置を求めることができる．

5

2.1.2

原子間ポテンシャル

系のポテンシャルエネルギー Etot の評価は，以下の 3 つに大別される．

(1) 経験的ポテンシャル(2) 半経験的ポテンシャル(3) 非経験的手法 (第一原理計算)
経験的ポテンシャルは，量子力学の厳密な理論に基づいて決定されるのではなく，ポテンシャルを微分可能な未定係数を含む簡単な関数形で仮定し，従来の実験的事実に合致するようにその未定係数が決められる．半経験的ポテンシャルは，密度汎関数論
より導出される形で定義されるが，そのポテンシャルパラメータは平衡状態でのマク
ロな物性値や，あるいは ab-initio な計算により求められた値に対してフィッティング
される．非経験的手法とは，従来の特性値などを一切用いず，原子核の位置ならびに種類のみを必要情報とし，各時刻における電子状態を量子力学に基づいて解くことで，遂次原子に働く力を精密に評価する手法である．

MD 法において，原子 i に作用する力 F i は系のエネルギー Etot の空間微分によって求めるため (式 (2.2))，系のポテンシャルエネルギー Etot をいかに精度よく評価するかが重要となる．(3) の第一原理分子動力学法は，計算量が極めて膨大になるため，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する (1)，(2) の原子間ポテンシャルが通常用いられる．本解
析の原子間ポテンシャルでは，Ni 基超合金には Daw，Baskes らによって提案された原子埋め込み法 (Embedded atom method; EAM)(30),(31) を用いている．

2.1.3

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

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

Etot =N∑

F (ρ̄i ) +

i

N ∑
N
1∑
φ(rij )
2 i j(6=i)(2.5)

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

た，ρ̄i は周囲の原子 j からの寄与 ρ(rij ) の重ね合わせで与えられると仮定しρ̄i =neighbor∑ρ(rij )

(2.6)

j(6=i)

で評価する．

2 種類以上の原子を含む系では，系のエネルギーは次式で表される．
Etot =N∑

Fti (ρ̄i ) +

i

ρ̄i =

N ∑
N
1∑
φt t (rij )2 i j(6=i) i j

neighbor∑ρtj (rij )

(2.7)

(2.8)

j(6=i)

ここで Fti (ρ̄) は原子種 ti の埋め込みエネルギー関数，φtitj (r) は原子種 ti と tj 間の 2 体間相互作用，ρtj (r) は原子種 tj の電子密度関数である．

Ni, Al の 2 原子系への適用を考える場合，FNi (ρ̄)，FAl (ρ̄)，ρNi (r)，ρAl (r)，φNiNi (r)，φAlAl (r)，φNiAl (r)，の関数形が必要となるが，
0
φNiNi(r) = φNiNi (r) − 2 gNi ρNi (r)(2.9)

FNi0 (ρ̄) = FNi (ρ̄) + gNi ρ̄Ni

(2.10)

なるパラメーター gNi (または gAl ) を導入して関数形の変換を行っても，式 (2.7) で Ni,

Al それぞれの純金属系におけるエネルギーは変化しない．さらに0ρAl(r) = sAl ρAl (r)(2.11)

FAl0 (ρ̄) = FAl (ρ̄/sAl )(2.12)のように Al の価電子密度をスケーリングするパラメーター sAl を導入して Ni 単原子系の電子密度の範囲に合わせる変換を施しても不変である．したがって，Ni, Al それぞれ単原子系で FNi (ρ̄)，FAl (ρ̄)，ρNi (r)，ρAl (r)，φNiNi (r)，φAlAl (r) を決定した後に，Ni–Al
合金系のエネルギーに対し φNiAl (r)，gNi ，gAl ，sAl を最適化することにより，Ni–Al 合
金系のエネルギーを正確に表すポテンシャル関数が決定できる．
本解析で対象とする Ni ならびに Ni3 Al については，Voter らが Ni，Al，Ni3 Al それ
ぞれ単結晶に対して昇華エネルギー，空孔形成エネルギー，弾性定数，格子定数等へのフィッティングを行い，ρ(r)，φ(r) の関数形として，

ρ(r) = s r6 (e−β r + 29 e−2 β r )7(2.13)

Table 2.1 Potential parameters for ρ(r)β (Å−1 )3.64083.3232Ni
Als
1.00000.6172

Table 2.2 Potential parameters for φ(r)
Ni-Ni
Al-Al
Ni-Al

α (Å−1 )1.77281.48591.6277

D (eV)
1.55353.77603.0322R (Å)
2.2053
2.11762.0896g (eVÅ3 )6.5145-0.2205
0

φ(r) = D {1 − exp[−α (r − R)]}2 − D − 2 g ρ(r)(2.14)

を提案している (32)∼(34) ．式中のパラメーターの値を表 2.1，2.2 に示す．
埋め込みエネルギー関数 FNi (ρ̄)，FAl (ρ̄) については，原論文中でその具体的な関数形は示されていないが，Foiles(36) が提案している Rose らの純金属系の凝集エネルギー関数 (35) を用いる方法で数値的に求めることができる．本研究では，この方法により埋め込みエネルギー関数を 3 次のスプライン関数でフィッティングした．フィッティ
ング範囲は 0.0 ≤ ρ̄ ≤ 1.0 で，スプラインノードの間隔は ∆ρ̄ = 0.01 とした．

2.1.4

速度スケーリング法

分子動力学法で温度制御する場合，もっとも簡単で直接的な方法として速度スケー
リング法がよく用いられる．熱統計力学より系の運動エネルギー K は次のように表さ
れる．

K=

N
1∑
3
mi (v i ·v i ) = N kB T
2 i=1
2

(2.15)

ここで，mi は原子 i の質量，v i は原子 i の速度，N は系の全原子数，kB はボルツマン
定数，T は系の温度である．式 (2.15) より，系の温度 T は原子速度を用いて，次のよ
うに求められる．

∑T =mi (v i ·v i )3N kB(2.16)

設定温度が TC ，式 (2.16) より求めたある時刻の温度が T のとき，速度スケーリング√

法では，各原子の速度 v i を

TC /T 倍し設定温度 TC に近づける．ベルレ法では，∆r i (t + ∆t) = r i (t + ∆t) − r i (t) = r i (t) − r i (t − ∆t) + (∆t)2
8

F i (t)m(2.17)

√をTC /T ∆r i (t + ∆t) で置き換えることに相当する．平衡状態では，能勢の方法 (37)など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている．

2.1.5

高速化手法

領域分割による高速化

N 個の原子からなる系では，Etot の評価に N × (N − 1) 回の原子対の計算が必要となる．一方，実際の結晶中では近接原子による遮蔽 (screening) 効果により第二近接距離程度より離れた原子はほとんど作用を及ぼさないことが知られている．このため，分子動力学計算では相互作用打ち切り (カットオフ) 半径 rc を導入し (図 2.1)，その半径内の原子からの寄与のみを考慮する．
しかしながら，相互作用する原子対の検索に N × (N − 1) 回の試行を要するため，系が大きくなるにつれ計算負荷が飛躍的に増加する．これを避けるために rc よりひとまわり大きい半径 rfc (図 2.1) 内の原子をメモリーに記憶し，rfc 内での原子対の探索とすることによりオーダー N の計算に近づける方法 (粒子登録法 (37) ) がこれまでよく用いられてきた．しかしながら，粒子登録法では rfc 半径より外の原子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1)回の探査) を行わなければならない．このため，系がある程度の規模以上に大きくなると，粒子登録による高速化は登録更新の負荷により打ち消される．
領域分割法では，まず図 2.2 に模式的に示すようにシミュレートする系をカットオ
フ距離程度の格子状に分割する．ある原子に作用する力を評価する際には，その原子が属する領域（図 2.2 の着色部）と隣接領域内 (図 2.2 の斜線部) の原子からカットオフ

rc

rfc

Fig.2.1 Schematic of bookkeeping method

9

y
by
0

x
bx

Fig.2.2 Schematic of domain decomposition
距離内の原子を探索する．原子が属する領域は，位置座標を領域ブロックの辺長 bx , by
で除した際の整数により判断できるので，領域分割そのものの計算負荷は小さい．領域分割法は，粒子登録法において登録更新の負荷が大きくなるような大規模な系の高速化に適している．

10

2.2

転位の Phase Field モデル

2.2.1

秩序変数と固有ひずみ

転位を定義する一つの方法は，結晶内のすべり面上において，すべった領域とすべっていない領域の境界を考えることである．本研究ではこの考えに基づき，すべりが生じた領域 (図 2.3 中の斜線) に以下の固有ひずみを与える．

ε0ij (α, mα ) =1[bi (α, mα )nj (α) + bj (α, mα )ni (α)]
2d(2.18)

ここで，α はすべり面番号，mα はすべり面 α 上のすべり方向番号である．bi (α, mα ), nj (α), dはそれぞれ，バーガースベクトル，すべり面法線ベクトル，すべり面間距離である．

Phase Field 法では，界面をシャープなものとして取り扱わず，秩序変数 (phase field)の変化領域として考える．本研究では，図 2.4 に模式的に示すように，すべり面上で転位がすべった領域で phase field パラメータ η を 1，それ以外で η = 0 とするとし，その変化する領域を転位芯と考える．この秩序変数を用いると，任意の転位形態を表す

b
n
b

d

Fig.2.3 Schematic of slip plane and Burgers vector.

11phase field

1

slipped planeslipped region
Interface

0.5

unslipped region
0
0

100

Fig.2.4 Phase field profile on slip plane.

固有ひずみの空間分布は次の式で表される．

ε0ij (r) =

p
q
∑
∑

ε0ij (α, mα )η(α, mα , r)(2.19)

α=1 mα =1

ここで，p, q はそれぞれすべり面の数，1 つのすべり面上のすべり方向の数である．

12

2.2.2

自由エネルギー汎関数

転位が存在することにより生ずる付加的なエネルギーを，結晶学的エネルギー E cryst ，勾配エネルギー E grad ，そして弾性ひずみエネルギー E elast の和として以下のように表す．

E = E cryst + E grad + E elast(2.20)

以下では，各エネルギー項について説明する．図 2.5 に模式的に示すように，単結晶中を転位が一格子分運動するとき，転位は原子の並びと同じ周期を持つポテンシャル
場を運動する．この周期的ポテンシャルはパイエルスポテンシャルと呼ばれ，図 2.5(d)のように変化する．このエネルギーを結晶学的エネルギー E cryst と考え，次式により
表す．
∫

E

cryst

ecryst dV

=
V

∫

=

p
q
∑
∑

V α=1 m =1

A sin2 (πη(α, mα , r)) dV(2.21)ατxy(b)x=b/2(c)x=bPeierls potential(a)x=0

0

b/2

b

x

(d)

Fig.2.5 Schematic of dislocation motion and corresponding Peierls potential.

13ecryst は結晶学的エネルギー密度である．係数 A は µ(εo )2 /2π 2 であり，界面領域におけるエネルギー障壁の高さを表している．µ はせん断係数である．秩序変数 η の変化による式 (2.21) の変化を図 2.6 に示す．図 2.6 より η が整数値のときエネルギーはゼロとなりその中間で最大値を取ることから，秩序変数 η が変化する転位芯部のみにおいてエネルギーが生ずることがわかる．
結晶学的エネルギーは局所的であるため近くの境の点との長距離相互作用がなく，秩序変数 η に連続性を与えない．Phase Field 法では，連続性は秩序変数 η の勾配項によって与えられる．また，すべり面垂直方向の秩序変数 η の変化領域においては，表面エネルギーが生じてはいけない (図 2.7)．以上の点を考え，次の勾配エネルギーを用
いる．
∫

E

grad

egrad dV

=
V

[

]

∂bi (α1 , r) ∂bk (α2 , r)1∫ ∑βijkl (α1 , α2 )=dV
2 V α1 ,α2
∂xj
∂xl
1
=
2

∫

p
p
∑
∑

q
∑

V α =1 α =1 m
1

2

[

q
∑

βijkl (α1 , α2 )bi (α1 , mα1 )bk (α2 , mα2 )α1 =1 mα2 =1]∂ηi (α1 , mα1 , r) ∂ηk (α2 , mα2 , r)
dV
∂xj
∂xl

(2.22)

ここで，egrad は勾配エネルギー密度である．βijkl は四次元テンソルであり，すべり面に平行なすべった領域の表面エネルギーを消すように与えられる．また，fcc 結晶におsin2(πη)

1

0.5

0
-2

-1

0

1

η
Fig.2.6 Profile of sin2 (πη).

142change in this direction: ∇η// n∴n×∇η=0

normal vector(slip plane)∇ηn(same to dislocation line)
n×∇η

η=1
∇η
n×∇η

η=0

∇η
gradient vector(normal to dislocation line)

Fig.2.7 Schematic of gradient term.

いて勾配エネルギーは以下のように近似できる．

E grad =β∫Φ(α, mα , r) dV
2 V(2.23)

ここで β は正の定数である．Φ(r) は以下のように与えられる．
([1̄10])[1̄10] 2([01̄1])[01̄1] 2([101̄])[101̄] 2([1̄1̄0])[1̄1̄0] 2

(

[101]

)
[101] 2([01̄1̄])[01̄1̄] 2

Φ(r)= n(111) × ∇η(111) + n(1̄1̄1) × ∇η(1̄1̄1)+ n(111) × ∇η(111) + n(11̄1̄) × ∇η(11̄1̄)+ n(111) × ∇η(111) + n(1̄11̄) × ∇η(1̄11̄)+ n(11̄1̄) × ∇η(11̄1̄) + n(1̄11̄) × ∇η(1̄11̄)+ n(11̄1̄) × ∇η(11̄1̄) + n(1̄1̄1) × ∇η(1̄1̄1)+ n(1̄11̄) × ∇η(1̄11̄) + n(1̄1̄1) × ∇η(1̄1̄1)(2.24)

弾性ひずみエネルギーは次式のように表される．
∫

E

elast

eelast dV

=
V

1∫
=
Dijkl (εij − ε0ij (α, r))(εkl − ε0kl (α, r)) dV
2 V(2.25)

ここで，eelast は弾性ひずみエネルギー密度，Dijkl は弾性係数である．
式 (2.20) に (2.21)(2.23)(2.25) を代入することによって FCC 結晶における自由エネ
ルギー汎関数を得ることができる．
格子サイズ dx がすべり面間距離 d より大きい場合，各エネルギーについてすべり面
間距離と格子サイズの違いを考慮しなければならない．E elast は，固有ひずみを格子サ
イズに修正した 0 = b/dx を用い，その他各項は，dx/d を除したものを使う (22)(23) ．

15

TDGL 方程式

2.2.3

秩序変数 η は非保存量であるため，時間に依存したギンツブルグ−ランダウ (TDGL)
方程式は次のように表される．

∂η
δE
= −Mη
∂t
δη(r)

(2.26)

ここで，Mη , t, はそれぞれ転位の易動度に関係するカイネティック係数と時間である．
自由エネルギー汎関数の微分値は次式のように表される．
(

δE∂ecryst ∂eelast∂ ∂egrad∂ ∂egrad∂ ∂egrad
=
+
−
+
+
δη(r)
∂η(r)
∂η(r)
∂x ∂ηx (r) ∂y ∂ηy (r) ∂z ∂ηz (r))(2.27)

ここで，ηx (r) などの下付き x, y, z はそれぞれ x, y, z による微分を表す．式 (2.27) の各
項は以下のように表される．

∂ecryst= 2πA sin(πη(r)) cos(πη(r))
∂η(r)
[
]
∂eelast1∂
0
0=
Dijkl (εij − εij (r))(εkl − εkl (r))∂η(r)∂η(r) 2][∂ε0ij (r)1∂ε0kl (r)00
= − Dijkl (εkl − εkl (r))+ (εij − εij (r))
2
∂η
∂η
[
]
1
= − Dijkl (εij − ε0ij (r))ε0kl (α, mα ) + (εij − ε0ij (r))ε0kl (α, mα )2
= −Dijkl (εkl − ε0kl (r))ε0ij (α, mα )
grad

grad

grad

∂ ∂e
∂ ∂e
∂ ∂e
β
+
+
=
∂x ∂ηx (r)∂y ∂ηy (r) ∂z ∂ηz (r)2((2.28)(2.29))∂ ∂Φ(r)∂ ∂Φ(r)∂ ∂Φ(r)++∂x ∂ηx (r) ∂y ∂ηy (r) ∂z ∂ηz (r)(2.30)

これより，式 (2,28)〜(2.30) を式 (2.27) に代入すると，TDGL 方程式が次のように表
される．
[

∂η(r)
β
= Mη
∂t
2

{

}

∂ ∂Φ(r)∂ ∂Φ(r)∂ ∂Φ(r)++− 2πA sin(πη(r)) cos(πη(r))∂x ∂ηx (r) ∂y ∂ηy (r) ∂z ∂ηz (r)]+Dijkl (εkl − ε0kl (r))ε0ij (α, mα )(2.31)TDGL 方程式 (2.31) において右辺第 1 項は拡散項を表し，界面の η の分布を滑らかにしようとする．残りの項は駆動力項となり η の分布を急峻にしようとする．これらがバランスする事によって η は図 2.4 のようなプロファイルを保ったまま界面の移動を
表現する．

16

結晶座標 (r, θ) は x-y 座標で次式で表される．
√x2 + y 2yθ = tan−1xr=(2.32)(2.33)式 (2.28)〜式 (2.30) を 2 次元 1 すべり系問題に対して表すと以下のようになる．

∂ecryst= 2πA sin(πη) cos(πη)(2.34)∂η∂eelast= −Dijkl (εkl − ε0kl )ε0ij(2.35)∂η()
22∂ ∂egrad∂ ∂egrad∂ ∂egrad
∂2η
2∂ η
2∂ η
+
+
= β n2 2 − 2n1 n2+ n1 2(2.36)∂x ∂ηx∂y ∂ηy∂z ∂ηz
∂x
∂x∂y
∂y
[ (
)
]
2
2
∂η
∂ 2η
2∂ η
2∂ η
0
0
= Mη β n2 2 − 2n1 n2+ n1 2 − 2πA sin(πη) cos(πη) + Dijkl (εkl − εkl )εij
∂t
∂x
∂x∂y
∂y
(2.37)

17

2.2.4

TDGL 方程式の数値解法

TDGL 方程式の数値解析を行うため時間に関しては前進差分法，空間に関しては二階中心差分法を用いて離散化する．図 2.8 に示す差分格子点 (i, j) において，式 (2.37)の時間および空間に関する項はそれぞれ以下のように差分表示される．
t+dttη(i,j)− η(i,j)∂η=∂t∆t
η(i+1,j) − 2η(i,j) + η(i−1,j)
∂ 2η
=
∂x2
(∆x)2
η(i,j+1) − 2η(i,j) + η(i,j−1)
∂ 2η
=
∂y 2
(∆y)2
η(i+1,j+1) − η(i−1,j+1) − η(i+1,j−1) + η(i−1,j−1)
∂ 2η
=
∂x∂y
4∆x∆y

(2.38)(2.39)(2.40)(2.41)

ここで，上付き添え字の t, dt はそれぞれ現在時刻および時間増分を表す．
応力場は有限要素法を用いて離散化する．内部に不均一に固有ひずみが分布した系の応力場を求めるために図 2.9 に示すような 2 次元線形弾性体を考え，以下に示す仮
想仕事の原理式，ひずみ−変位関係式，応力−ひずみ関係式，幾何学的境界条件からなる境界値問題を設定する．
仮想仕事の原理式∫∫SS∫(σx δεx + σy δεy + τxy δγxy )dS − (Px δu + Py δv)dl = 0(2.42)

δ{ε}T {σ}dS −(2.43)(i -1, j +1)

l

∫

l

δ{u}T {P }dl = 0(i , j +1)(i +1, j +1)

Computational grid(i -1, j )(i , j )(i +1, j )dyy(i -1, j -1)(i , j -1)(i +1, j -1)
dx

Fig.2.8 Finite diﬀerence lattice.

18

x

Pi
Sp

y

Su
0

ui

x

Fig.2.9 Two-dimensional linear elastic system.

ひずみ−変位関係式

1
εij = (ui,j + uj,i )2(2.44)

σij = Dijkl (εij − ε0ij )(2.45)

ui = U(2.46)

応力−ひずみ関係式

幾何学的境界条件(on Su )図 2.10(a) に示すような要素内の任意の点における変位 u, v を形状関数 Ni を用いて，
節点 ui , vj の線形結合によって表示する．

u = Ni ui(2.47)

v = Ni vi(2.48)

また，本研究では図 2.10(b) に示す 4 節点アイソパラメトリック要素を採用する．このとき，形状関数 Ni は局所座標 (ξ, η) の関数として以下のように表される．

1
N1 = (1 + ξ)(1 + η)
4
1
N2 = (1 − ξ)(1 + η)
41
N3 = (1 − ξ)(1 − η)
41
N4 = (1 + ξ)(1 − η)
4
19

(2.49)(2.50)(2.51)(2.52)

η4(x4,y4)( ξ 4, η4) = (-1,1)( ξ 3, η3) = (1,1)

4

3

3 (x3,y3)(x,y)( ξ ,η )ξoy

2(x2,y2)1(x1,y1)1

2( ξ 1, η1) = (-1,-1)( ξ 2, η2) = (1,-1)x(a) Global system(b)Local system

Fig.2.10 Four nodes isoparametric element.

式 (2.44) に式 (2.47)，式 (2.48) を代入すると，

∂u
∂N1
∂N2
∂N3
∂N4
=
u1 +
u2 +
u3 +
u4
∂x
∂x
∂x
∂x
∂x
∂N1
∂N2
∂N3
∂N4
∂v
=
v1 +
v2 +
v3 +
v4
εy =
∂y
∂y
∂y
∂y
∂y
∂v ∂u
∂N1
∂N1
∂N2
∂N2
γxy =
+
=
v1 +
u1 +
v2 +
u2
∂x ∂y
∂x
∂y
∂x
∂y
∂N3
∂N4
∂N4
∂N3
v3 +
u3 +
v4 +
u4
+
∂x
∂y
∂x
∂y

εx =

(2.53)(2.54)(2.55)を得る．これをマトリックス表示すると次式を得る．




 εx

ε

y


 γ

xy








N1,x


=
 0



N

1,y

0 N2,x
0 N3,x
0 N4,xN1,y0 N2,y0 N3,y0
N1,x N2,y N2,x N3,y N3,x N4,y u1 











v
1










u


2




0





v
2

(2.56)
N4,y  
u3 






N4,x 

 v3 









u

4 









v4

ここで




N1,x
0 N2,x
0 N3,x
0 N4,x0[B] = N1,y
0 N2,y0 N3,y0 N4,y  0
N1,y N1,x N2,y N2,x N3,y N3,x N4,y N4,x20

(2.57)

と定義すると{ε} = [B]{ue }(2.58)となり，要素内のひずみは節点の変位によって表される．
またフックの法則より，要素内の応力は以下に表される．

{σ} = [De ]({ε} − {εo })(2.59)= [De ][B]{ue } − [De ]{εo }(2.60)

ここで [De ] は弾性係数マトリクスである．
式 (2.47)，式 (2.48)，式 (2.58) を式 (2.43) に代入すると∫
S

δ{ue }T [B]T {σ}dS −∫δ{ue }T [N ]T {P }dl = 0l(2.61)

ここで，任意の δ{ue } に対して成立するため，∫∫[B] {σ}dS − [N ]T {P }dl = 0
TS(2.62)lを得ることができる．また，この式に式 (2.60) を代入すると，∫[B] [D][B]dS{u }dS =
T

S

∫

e

∫

[B] [D]{ε }dS + [N ]T {P }dlTSo(2.63)l[k e ] {ue } = {foe } + {fpe }(2.64)

を得ることが出来る．この式を整理すると要素剛性方程式を得ることができる．
∫[k e ] =[B]T [D][B]dS(2.65)[B]T [D]{εo }dS(2.66)[N ]T {P }dl(2.67)

S{foe } ={fpe } =∫∫Sl

ここで [k e ], {foe }, {fpe } はそれぞれ要素剛性マトリックス，固有ひずみの影響を等価な節点力として表したもの，外力に等価な節点力である．
要素剛性マトリックスに対して，次式に示す全体座標系 (x, y) と局所座標系 (ξ, η) における関数 F (x, y) の要素全域の積分の関係式∫∫∫ 1 ∫ 1F (x, y)dxdy =−1−1

F (ξ, η)|J|dξdη

(2.68)

を用いると，式 (2.65) は以下のように表示できる．
∫∫[k e ] =[B]T [D][B]dxdy∫ 1 ∫ 1=−1

−1

[B]T [D][B]|J|dξdη
21(2.69)ηf (ξ )

Gauss integration points

w1

w2

1
1

-1

√
3
‑1
√
3

o

1

√
3

1

ξ

‑1
‑1
√
3

1
√
3

√
3

-1

1

ξ

-1

Gauss integration points

Fig.2.11 Two points Gaussian integration.

さらに，次式に示す 2 点ガウス積分を行う．
∫ 1 ∫ 1
−1

−1

2 ∑
2
∑

f (ξ, η)dξdη =f (ξi , ηj )wi wj(2.70)

i=1 j=1

ここで，wi , wj は図 2.11 に示す重み，添え字 i, j は積分するガウス点を示している．これより，要素剛性方程式は次のように与えられる．

[k e ]{ue } = {foe } + {fpe }[k e ] =2 ∑2∑(2.71)[B]Tij [D]ij [B]ij |J|ij wi wj(2.72)

i=1 j=1

要素剛性マトリックスを全要素について重ね合わせることで全解析領域に対する全体剛性マトリックスと有限要素方程式を得る．

[K] {U } = {Fo } + {Fp }
[K] =

∑

(2.73)[k e ](2.74){foe }(2.75)

element{Fo } =∑

element{Fp } =∑

{fpe }

element

22(2.76)

2.3

NEB(Nudged Elastic Band) 法

NEB(Nudged Elastic Band) 法は，初期状態と最終状態の 2 つの状態を結ぶ反応経路の中で，エネルギーが最も低くなる経路を探索する手法である．N 個の原子からなる系の状態は，すべての原子の位置座標を軸とする 3N 次元の位相空間において，位~ (x1 , y1 , z1 , ・ ・ ・) で表される．位相空間における初期状態 R~ s およ置ベクトル R=~ e の 2 点を結ぶ直線経路 R~e − R~ s を複数の位相点に分割し，各位相点をび最終状態 R
仮想的なばねでつないだ状態で，位相点にかかる力に基づいて位相空間内における低~ 0 は，経路に垂直ないエネルギー経路を探索する (図 2.12)．各位相点に作用する力 Fi

ポテンシャルエネルギーの勾配と経路に平行なバネ力の和として次式で表す．

~ (R~ i )|⊥ + (F~ s・τ̂k )τ̂k
F~i0 = −5Vi(2.77)~ i は i 番目の位相点の位置ベクトル，F~ s は i
ここで，V はポテンシャルエネルギー，Ri~~ i)番目の位相点に作用するバネ力，τ̂k は経路の接線方向単位ベクトルである．−5V (R~ i の原子配置において，各原子に働く力の 3N 個の成分そのものである．位相は状態 R空間内で，ばねに平行な方向でエネルギーが急激に変化する場合，式 (2) では仮想的なばねの力のみが支配的になる．そこで隣接する角度によって変化する関数 f (φi ) を導入し，次式で各位相点に働く力を表現することが提案されている．

F~iN EB = F~i0 + f (φi )(F~is − (F~is・τ̂k )τ̂k )EndEnergy state(2.78)End

Energy state

Start

Start

(b) With spring(a) Without springFig.2.12 NEB method.

23

1
f (φi ) = (1 + cos(π(cos φi )))2(2.79)

経路に垂直方向の勾配は位相空間における「谷」を表し，経路に平行方向の勾配は谷底を流れる「川」に沿った方向の高低に対応する．ばねの力によって，位相点が「川」に沿ってすべて低いところに流れてしまうのを防ぎ，各位相点における「川」の高低を求めることができる．

24

第3章

EAM ポテンシャルによる
APB エネルギー評価
第一原理計算によれば，Ni3 Al の APB エネルギーは 126mJ/m2 と見積もられている(27)．一方，第一原理計算では計算量が膨大となるため，原子の熱揺動を考慮できるよ

うな動的な計算は現時点では不可能である．そこで本章では，熱膨張を想定して平均格子間隔を広げて静力学的に解析した時の APB エネルギーと，実際に熱揺動の効果を考慮したときのそれの違いを議論するために，EAM ポテンシャルによる静力学解析な
らびに動力学解析を行う．

3.1

逆位相境界・超転位

Ni3 Al の L12 構造は，fcc の Ni 単位格子において頂点原子を Al に置換した構造をとる (図 3.1)．したがって，fcc のすべり方向 [110] を見ると，並進対称性を満足する最Ni atom
Al atom

APB
fcc

L12 structure

Fig.3.1 Schematic of anti-phase boundary in L12 structure.

25短距離は，基礎となる fcc 構造におけるそれの 2 倍となる．すなわち，fcc の完全転位のバーガースベクトル b = a/2h110i では元の結晶とは一致せず，図 3.1 右に模式的に示すようにすべり面上下で位相のずれを生じる．その境界を逆位相境界 (Anti Phase

Boundary;APB) と称する．
積層欠陥と同様に，面欠陥である APB が導入されると系のエネルギーは高くなるので，転位と同様熱力学的には安定に存在しえないものである．一方，1 本目の転位が経過した後に生じる APB は，同じバーガースベクトルを持つ 2 本目の転位が通過することで解消される．fcc における拡張転位と同様に，これら 2 本の転位の間に挟まれた

APB の幅は転位間の反力と APB エネルギーのバランスにより決まる．この 2 本の転位と APB を一組のものとして超転位と呼ぶ (図 3.2)．また，Shockley の部分転位にならい 1 本目の転位を super leading partial，2 本目を super trailing partial と称する．

a/2 [110]a/6 [121] a/6 [211] Stacking fault

Dissociation
Slip plane

Perfect dislocation

Partial dislocation(a) extended dislocation in fcc.
a [110]a/2 [110]a/2 [110]

APB
Dissociation
Slip plane

Super dislocation

Dislocation(b) superdislocation in L12 structure.
Fig.3.2 Schematic of extended dislocation and superdislocation.

26

3.2

静力学解析

3.2.1

解析条件

Ni3 Al 相の (111) 面を xy 平面とし，z 軸方向に積み重ねた直方体セル (5.04nm×5.22
nm×39.36nm) を解析セルとして用いた．この直方体セルを下段，中段，上段と分け，中段の原子を [-110] 方向に 1 原子分ずらすことで APB を 2 つ導入する (図 3.3)．原子数は 92,160 である．APB を導入する前の系 (完全単結晶) と，導入した系それぞれについて，熱膨張を想定して格子長さを変化させ，原子を動かすことなく全方向に周期境界条件の下でエネルギーを計算し，APB エネルギーを次式で算出した．

ΓAPB = (EAPB − Ebulk )/2AAPB(3.1)

ここで，ΓAPB ，Ebulk ，EAPB は APB の持つ単位面積当りのエネルギー，APB を含まない系全体のエネルギー，APB を含む系全体のエネルギーである．また AAPB は APBの面積であり，図 3.3 のセルは APB を 2 面含むため 2AAPB となる．
2nm

5.0

9.84nm

19.68nm

9.84nm

Ni
Al

[111]
z

[112] y

[110]
x

5.21n

m

Fig.3.3 Simulation model for static analysis.

27

3.2.2

解析結果及び考察

図 3.4 に APB を有する系と，単結晶 Ni3 Al の格子長さ変化に対するエネルギー変化を示す．いずれも，体積（格子長さ）が増大するとエネルギーが増大する．当然ながら，APB を含む構造の方が常に系のエネルギーが高い．これらの差から式 (3.1) により算出した APB エネルギーを図 3.5 に示す．体積 (格子長さ) の増加に従い APB エネ
ルギーは減少している．これが静力学解析のみで熱膨張を想定したときの APB エネ
ルギー変化であり，図 3.5 の左端が 0K の APB エネルギーに相当し，右にいくに従い「高温側」となる．

(×10‑16)-2670

Potential energy [J]

Bulk
APB
-2675

-2680

-2685

-2690
1.075

1.100

Volume [ A3]

1.125(×106)

Fig.3.4 Potential energy by static analysis.

28

APB energy [mJ/m2]

200

180

160

140

1.075

1.100
3

Volume [ A ]1.125(×106)

Fig.3.5 APB energy by static analysis.

29

3.3

動力学解析

3.3.1

解析条件

静力学解析で用いたシミュレーションセルを x，y 軸方向にそれぞれ倍としたシミュ
レーションセル (10.08nm×10.43nm×39.36nm) を用いた (図 3.6)．APB を導入する前の系と，導入した系それぞれについて，全方向に周期境界条件の下で 15000step の緩和計算を行い，そのエネルギー差から APB エネルギーを算出した．温度は速度スケー
リング法により制御し，10K，100K，300K，1000K の 4 通りについて検討した．また，垂直応力が 0 となるようにセル辺長を制御している．エネルギー評価は，緩和によるエネルギー減少が収束する 10000〜15000step の時間平均により行った．
8nm

0
10.

9.84nm

19.68nm

9.84nm

Ni
Al

10.43

nm

[111]
z

[112] y[110]xFig.3.6 Simulation model for dynamic simulation.

3.3.2

解析結果及び考察

動力学解析により得られた単結晶 Ni3 Al のエネルギー，APB を含有する系のエネル
ギー，そして式 (3.1) により算出した APB エネルギーを表 3.1 に示す．APB エネルギーの値は温度により大きく変わることなく，145.6±2.2mJ/m2 であった．実験的な検討では，Ni3 Al の L12 の（111）APB エネルギーは，Karnthler らは 175±15mJ/m2(39) ，Sunらは 144±20mJ/m2(40) と見積もっており，特に後者には比較的近い値となっている．

30

Table 3.1 Antiphase boundary energy by EAM potential.
10K
100K
300K
1000K

EBulk (J)×10−16 EAPB (J)×10−16 αAPB (mJ/m2 )-2689.10063-2688.78991
145.3-2682.17078-2681.86014
144.4-2666.42356-2666.11130
144.9-2600.49401-2600.16611
147.8

図 3.7 は静力学解析から計算した APB エネルギーの値 (図 3.5) に, 動力学解析で求めた各温度における体積と APB エネルギーをプロットしたものである．ただし，静力
学解析で用いたシミュレーションセルは動力学解析のそれの 1/4 なので，ここでは 4 倍
して同じ大きさで示している．静力学解析での APB エネルギーは動力学解析より常に高く，低温側 (左側) でその差が大きくなっている．そこで，10K において静力学解析と動力学解析の APB 近傍の原子のエネルギーを調べた (図 3.8)．図 3.8 の縦軸のエネ
ルギーは，図 3.9 に示すように，APB を境として各 (111) 面に属する原子のエネルギーを平均し，APB を導入していない単結晶の平均エネルギーからの差をとったものである．静力学解析では APB から上下 2 原子層のエネルギーが高く，3 層目以上は影響がない．一方，動力学の場合は 3 層目にもエネルギーが高い所が生じているものの，1，2 層目のエネルギーが大きく緩和されていることがわかる．また，10K と 1000K における layer 毎のエネルギーをそれぞれ静力学解析，動力学解析毎に示した (図 3.10(a)，(b))．静力学解析 (図 (a)) では，APB 近傍のエネルギー分布のプロファイルは変わらないまま，各層のエネルギーが減少している．一方，動力学では 1 層目のエネルギーはほとんど変わらず，温度が上昇すると 2 層目以降のエネルギーが高くなり，かつ影
響範囲が 4 層目まで及んでいる．この APB 近傍のエネルギー変化だけみると，動力学解析では温度が上昇すれば APB エネルギーは増加することになり，表 3.1 に示した結果と矛盾する．しかしながら，温度上昇とともに体積が増加し，式 (3.1) の分母が大きくなる．すなわち，動力学解析で APB エネルギーがさほど変化しないのは，APB 近
傍でのエネルギー上昇と体積増加が相殺されたためである．

31

APB energy [mJ/m2]

200

10K
100K
300K
1000K
by MD

180

160

140
4.3

4.4

4.5
(×106)

Volume [ A3]

Fig.3.7 APB energy by statics analysis.

static analysis(10K)
dynamic analysis

Energy [eV]

APB

0.02

0.01

0
-4

-2

2

4

Layer No.

Fig.3.8 Energy profile of each atomic layer near the APB.

32

Layer No.
・
・
・
3
2
1
-1 APB
-2
-3
・
・
・

[111]
z
[110]
x

Fig.3.9 Schematic of layer number.

0.03

0.03

APB

10K
1000K

0.02

Energy [eV]

Energy [eV]

APB

0.01

0

10K
1000K

0.02

0.01

0
-4

-2

2

4

-4

Layer No.

-2

2

4

Layer No.
(b) Dynamic analysis(a) Static analysis

Fig.3.10 Energy each atomic layers in (a) static and (b) dynamic analysis.

33

3.4

結言

本章では，熱膨張を想定するものとして静力学的に平均格子間隔を変化させて算出した時の Ni3 Al の APB エネルギーと，実際に熱揺動の効果を考慮した場合のそれの違いを検討するため，EAM ポテンシャルによる静力学解析ならびに動力学解析を行った．得られた結果を総括して以下に示す．

1. (111) 面の APB を 2 つ含む直方体周期セルを用いて，構造緩和を行わず静力学
的に APB エネルギーを求めた．熱膨張を想定して体積 (格子長さ) を増加させると，APB エネルギーは 0K の格子長さにおける値 180.4mJ/m2 からなだらかに
単調減少した．

2. 同様のシミュレーションセルを用いて，10〜1000K の温度条件下で分子動力学解析を行ったところ，APB エネルギーはいずれも 145.6±2.2mJ/m2 となり，温度による大きな変化は見られなかった．この値は実験的に求められている値は175±15mJ/m2(39) ，144±20mJ/m2(40) と比較的近い．
3. 構造緩和しない静力学解析と動力学解析の APB エネルギーの差は，10K で 34.5mJ/m2 ，1000K で 9.1 mJ/m2 と極めて大きく，特に低温側で顕著となった．
4. APB の上下の原子層毎のエネルギーを見ると，静力学解析では 1，2 原子層のみ単結晶のエネルギーより高い値をとっており，格子長さを広げると分布形状は変わらず各層のエネルギーが少しだけ減少する．

5. 動力学解析では，APB の上下 1 層目のエネルギーは 10K，1000K ともほぼ同じ値を示し，2 層目以降は 1000K の方が 10K よりもエネルギーが高く，かつ広範
囲 (4 層目以上) に広がっていた．

6. 5 の結果だけみると高温ほど APB エネルギーが高いことになる．しかしながら，
2 で示したように APB エネルギーは温度上昇に対して変化しないので，APB 近
傍でのエネルギー上昇は，熱膨張による体積増加によって相殺されるものと考え
られる．

34

第4章

APB エネルギーを考慮した
Phase Field モデル
4.1

APB エネルギー項の導入

Ni 基単結晶超合金の転位挙動をモデル化するためには，γ 0 相中の APB を考慮しなければならない．2.2.2 節で説明したように，完全転位の Phase Field モデルにはパイ
エルスポテンシャルに相当する E cryst 項がある．これを参考に，APB によって発生するエネルギーを考慮した転位の Phase Field モデルを提案する．すなわち，APB エネ
ルギー E APB を考慮した場合の転位を含む系のポテンシャルエネルギーを次式で表す．

E = E cryst + E grad + E elast + E APB(4.1)

γ 0 相に侵入する一本目の転位（転位 1）は γ 0 相で APB を形成させ，転位が侵入するために APB エネルギーだけの過剰エネルギーが必要となる．一方，転位 1 が通過した
後の γ 0 相に侵入する二本目の転位（転位 2）は APB を解消するため，転位 2 は APB
エネルギーだけ減少する．したがって，E APB は，転位 1,2 が γ 相もしくは γ 0 相に存在する場合によって，次のような関数を用いて導入する．

E APB = 0∫

E APB =

p
∑

q
∑

αAP B
π
sin2 ( η(α, mα , r)) dVd2
V α=1 m =1α

(in γ phase)(4.2)(in γ 0 phase)(4.3)

ここで，αAPB は単位面積当りの APB エネルギーである．γ 相にある転位の E APB は
もちろんゼロである．γ 0 相内の差分格子点では，図 4.1 に示すように秩序変数の値によって E APB が変化する．1 本目の転位が通過した点では，phase field パラメータ η は1 となる．fcc 結晶に対するパイエルスポテンシャル式 (2.21) では E cryst は 0 にもどる
35

が，γ 0 相内では E APB は最大となる．逆に 2 本目の転位が同一すべり面を通過すると

η が加算され 2 となり，E APB が 0 となる．
APB エネルギー項を考慮した自由エネルギー汎関数の微分値は次式のように表せる．
()

δE∂ecryst ∂eelast∂ ∂egrad∂ ∂egrad∂ ∂egrad
∂eAPB
=
+
−
+
+
+
(4.4)
δη(r)
∂η(r)
∂η(r)
∂x ∂ηx (r) ∂y ∂ηy (r) ∂z ∂ηz (r)∂η(r)
ここで，egrad は APB エネルギー密度である．そして，上式の APB エネルギー項は以下のように表せる．


0




∂eAPB
=

∂η(r)



 π

(in γ phase)
αAP B
π
π
sin( η(r)) cos( η(r))
d
2
2

(4.5)
0

(in γ phase)

式 (2.26) に式 (4.4) を代入すると，TDGL 方程式は次のように表される．
(

(

)

)

∂η(r)
∂ ∂egrad∂eAPB∂ecryst ∂eelast∂ ∂egrad∂ ∂egrad
+
= −Mη
+
−
+
+
∂t
∂η(r)
∂η(r)
∂x ∂ηx (r) ∂y ∂ηy (r) ∂z ∂ηz (r)
∂η(r)
(4.6)

sin2(πη /2)

1

0.5

0
0

1

2

Phase field ηFig.4.1 Phase field profile of sin2 (πη/2) for superdislocation.

36

4.2

γ 0 単相中の超転位のシミュレーション

4.2.1

解析条件

2 次元 1 すべり系 (fcc) とした γ 0 相単相の解析モデルを図 4.2 に示す．ここで転位は紙面垂直方向に無限に長い直線転位を想定している．解析領域は 117.6nm × 42.8nm(140× 51 差分格子) とし，平面ひずみ問題を仮定する．差分格子の間隔は dx = 0.84nm である．x，y ，z 方向の結晶方位をそれぞれ [1̄10]，[111]，[1̄1̄2] とし，FEM における境界条件は側面と上面が自由境界，下面は図 4.2 のように下端を固定する．また，秩序変数 η については全方向零ノイマン条件とした．転位の初期配置 (η の初期プロファイル) として，中央の差分格子の値を左から 40 までは η = 2，41 から 100 までは η = 1，
その他の差分格子の値は η = 0 とした．縦弾性係数 E は 205GPa, ポアソン比は ν=0.3とした．格子定数 a の値からすべり面間距離 d，バーガースベクトルの大きさ b を導出
し，それぞれの値は a=0.35nm，d=0.42nm，b=0.24nm とする．その他の無次元化したパラメータについては，勾配係数 β ∗ =0.1，パイエルスポテンシャル A∗ = 0.05，易動度 Mη∗ =1.0 とする．また，単位面積当りの APB エネルギー αAPB は前章の動力学解析で算出した 145.6mJ/m2 とした．時間ステップ 4t∗ = 0.2 を 1 ステップとして η の発展を解き，APB を考慮した γ 0 相単相での転位の挙動を調べた．

yzγ´ phase

x

[112]

[110]

Fig.4.2 Simulation model for dislocations in γ 0 phase.

37

Shear free boundary condition[111]

42.8nm (51 lattices)

117.6nm (140 lattices)

4.2.2

解析結果及び考察

図 4.3 に秩序変数 η の各時間ステップでの分布を示す. 秩序変数 η が非整数の値を持
つ遷移領域 (0 ＜η ＜ 1，1 ＜η ＜ 2 の領域) が転位芯領域である．10000step までは，前方 (右側) の転位 1 (0 ＜η ＜ 1) は左に移動し，転位 2 (1 ＜η ＜ 2) が右に移動して両者の間隔が狭まっている．これは，転位 1 は，図 4.1 の 0 ＜η ＜ 1 における正のエネルギー勾配，一方転位 2 は，図 4.1 の η >1 における負のエネルギー勾配からわかるように，転
位 1，転位 2 が運動して APB を解消する方がエネルギー的に安定となるので両者の間隔が狭まっている．転位 1 と転位 2 がある程度接近すると弾性相互作用 (反力) を生じるため，10000step 以降の図では転位間の距離を一定に保ち, 超転位として一組となっている．
図 4.4 に各時間ステップにおけるせん断応力 τxy の分布を示す．正負対の高いせん断応力を示している領域が転位芯の位置である．10000step において APB をはさんで超転位になると転位間の応力分布 (転位 1 の左側と転位 2 の右側) は，互いに打ち消しあい，両転位の外側 (転位 1 の右側と転位 2 の左側) は強くなっている．
図 4.5(a) は，転位 1 と転位 2 の位置の時間変化を，図 4.5(b) は転位間の距離の時間変化を示したものである．5000step 近傍までは先述のように転位 1 と転位 2 は一定速度で近づいている．5000〜7000step において，前述のように転位同士が接近すること
による弾性相互作用を生じ，転位 1 と転位 2 は接近する速度を緩め，一定間隔で停止する．
本モデルの妥当性を検証するために，APB エネルギーが 120.0，145.6，170.0mJ/m2の時の転位全体のエネルギーから導出される刃状転位 2 本の平衡間隔の理論値 d と解析結果を比較する．転位論より理論値 d を以下に示す．

d=

µb2
1
1
・
・
2π αAP B (1 − ν)(4.7)

ここで，µ は横弾性係数である．図 4.6 に，理論値と解析結果を示す．各印は APB エ
ネルギーが 120.0，145.6，170.0mJ/m2 の解析結果に対応している．両者の値はほぼ一致しており，本モデルが超転位を良好に再現できている．

38

500step

5000step

10000step

15000step
0.0

0.5

1.0

1.5

2.0

Phase field ηFig.4.3 Distribution of phase field in γ 0 phase.

39

500step

5000step

10000step

15000step
-2320

-1160

0

1160

2320

Shear stress τ xy [MPa]Fig.4.4 Distribution of shear stress τxy in γ 0 phase.
40

Position of dislocations [b]350dislocation1
dislocation2

300

250

200

150
0

5000

10000

15000

Time steps Δt*

Distance between dislocations [b]

(a)

200

150

100

50

0

5000

10000

15000

Time steps Δt*(b)
Fig.4.5 (a) Change in position of dislocation and (b) distance between dislocation1 and dislocation2.

41

Distance between dislocations d, [b]: αAPB=120.0 mJ/m2: αAPB=145.6 mJ/m2: αAPB=170.0 mJ/m2(simulation result): Dislocation Theory

40

30

20
100

120

140

160

180

200

APB energy [mJ/m2]
Fig.4.6 Change in distance between dislocations according to APB energy.
Solid line shows the result of dislocation theory and open symbols
indicate the value evaluated in the phase field simulation.

42

4.3

γ 0 相をカッティングする転位のシミュレーション

4.3.1

解析条件

前節と同様，2 次元 1 すべり系 (fcc) として半分を γ 相としたモデル (図 4.7) を対象とする．解析領域は転位進行方向を少し大きくとり，252.0nm × 42.8nm(300 × 51 差分格子) とした．結晶方位や上面以外の境界条件等は前節と同じである．上面にはせん断応
力 τ を与える．中央のすべり面に置く転位は，左側面からそれぞれ 70, 140 差分格子の位置に配置した．γ 相と γ 0 相の縦弾性係数 E の差は小さいのでいずれも E = 205GPaとし，ポアソン比も同様に ν=0.3 とした．その他のパラメータは，前節と同じとした．
与えるせん断応力は 300MPa と 400MPa の 2 種類を考慮し，界面における挙動について調べた．

126.0nm (150 lattices)126.0nm (150 lattices)

y

z
[112]

γ phaseγ´ phaseγ/γ´ interface

x

[110]

Fig.4.7 Simulation model for dislocations at γ/γ 0 interface.

43

Shear free boundary condition[111]

42.8nm (51 lattices)

Shear stress4.3.2

解析結果及び考察

図 4.8，図 4.9 はそれぞれ τ = 300MPa および 400MPa を与えたときの phase field ηの変化である．γ と γ 0 の固相の違いをさらに別の phase field パラメータで表せば界面の形態変化も追跡可能であるが，今回のシミュレーションでは導入していないので界面の位置を破線で示した．τ = 300MPa の図 4.8 において，γ 相内にある転位 1, 2 は右に移動するが，転位 1 が界面に到達するとそこで停止する (図 (a))．その後も転位 2 は右に運動しつづけ転位 1 に接近するが，35000〜45000step では転位の幅が一定になったまま静止した. 一方，τ =400MPa を与えた図 4.9 では，図 4.8 と同様，転位 1 は界面で停滞するが，転位 2 は τ =300MPa のときより転位 1 に接近し，その後転位 1 が押し出される形で界面から γ 0 相に侵入して運動しつづけた (図 (d)〜(e))．
せん断応力 τ =300, 400MPa を与えたときの各時間ステップにおける転位 1, 2 の位置を図 4.10 に示す．図中には破線で界面の位置を示している．τ =300MPa の図 (a) では，
転位 1 が界面で静止している．また転位 2 は転位 1 との弾性相互作用によって反力を受け界面からの位置で停滞している．なお，このときの値は約 75b であるが，fcc の γ

phase での 2 本の刃状転位間の距離であり，前節の APB を間に含む γ 0 相中の超転位の間隔約 30b とは全く別の現象であることに注意されたい．一方，せん断応力 τ =400MPaを与えた図 (b) では，20000step 近傍で接近した転位 2 からの力を受け，転位 1 が界面を越えて γ 0 相に侵入している．γ 0 相内では転位 1，2 は一定の距離 (約 30b) を保ち超転位として運動している．

44

γ phaseγ' phase(a)5000step(b)15000step(c)25000step(d)35000step(e)45000step
0.0

0.5

1.0

1.5

2.0

Phase field ηFig.4.8 Distribution of phase field under shear stress of τxy = 300MPa.
45γ phaseγ' phase(a)5000step(b)15000step(c)25000step(d)35000step(e)45000step
0.0

0.5

1.0

1.5

2.0

Phase field ηFig.4.9 Distribution of phase field under shear stress of τxy = 400MPa.

46

Position of dislocations [b]800
dislocation1
dislocation2

700

interface
600

500

400

300

0

10000

20000

30000

40000

Time steps Δt*(a) τ =300 [MPa]Position of dislocations [b]800
dislocation1
dislocation2
700

interface
600

500

400

300

0

10000

20000

30000

40000

Time steps Δt*(a) τ =400 [MPa]Fig.4.10 Change in position of dislocation with shear stress of (a) τ = 300MPa
and (b) τ = 400MPa.
47

4.4

複数の転位を導入した

γ 0 相カッティングシミュレーション
4.4.1

解析条件

APB エネルギーを考慮した離散転位動力学シミュレーションにおいて，接近した後続の転位からの力を受け，低い外部応力下で γ 0 相へのカッティングを生じることが報告されている (41) ．そこで，節 4.3 でのせん断応力 τ =300MPa のモデルにおいて，同
一すべり面上にもう 1 本転位を配置したシミュレーションを行った．中央のすべり面に置く転位は，左側面からそれぞれ 60，120，140 差分格子の位置に配置した．

4.4.2

解析結果及び考察

図 4.11 は刃状転位 3 本を同一すべり面上に配置し，せん断応力 τ =300MPa を与えたときの phase field η の変化である．また，図 4.12(a) は転位 1，転位 2，転位 3(2 < η < 3)の位置の時間変化を，図 4.12(b) は転位 1 と転位 2 の転位間の距離の時間変化を示したものである．先ほどの結果では転位 1，転位 2 とも界面で停滞していた．しかし，今回3 本目の転位を導入し，その転位が界面に近づくことによって弾性相互作用が生じる
ため，10000step 付近で，転位 1 は γ 0 相に侵入し，転位 2 は界面に接近している．その時の転位 1，転位 2 の転位間は一定 (約 58b) である．転位 2 が界面に到達すると容易に γ 0 相に侵入し，その後，両転位間一定の距離 (約 30b) を保ち超転位として運動している．転位 3 は，界面に到達すると APB のエネルギーバリアによって停止している．
バックフォースモデルを用いた離散転位動力学でも同様の結果を示している (41) ．

48

γ phaseγ' phase(a)5000step(b)10000step(c)20000step(d)30000step(e)40000step
0.00

0.75

1.50

2.25

3.00

Phase field ηFig.4.11 Change in position of three dislocations with shear stress τxy =300MPa.
49

dislocation1
dislocation2
dislocation3

Position of dislocations [b]

1000

800

interface
600

400

200
0

20000

40000

Time steps Δt*

Distance between dislocations [b]

(a)

80

60

40

0

20000

40000

Time steps Δt*(b)
Fig.4.12 (a) Change in position of three dislocations and (b) distance between
dislocation1 and dislocatin2.

50

4.5

結言

Phase Field 法による転位のシミュレーションにおいて，γ 0 相の逆位相境界（APB）のエネルギーを考慮した超転位のモデル化を行った．転位の Phase Field 法で用いられている結晶学的エネルギー項 E cryst のパイエルスポテンシャルを参考に，1 本目の転位には前章の分子動力学解析で求めた APB エネルギー分の障壁を生じ，2 本目はそれを解消するような E APB 項を導入したシミュレーションを行い，以下の結果を得た．

1. まず 2 次元平面ひずみで，L12 (fcc) の 1 すべり系を仮定した無限直線状の刃状転位について，γ 0 単相中の超転位のシミュレーションを行った．約 200b 程度の間隔で 2 本の転位を配置したところ，APB の幅を狭めるように互いに接近し，APBの幅が約 30b となったとき，両者は転位間の距離を一定に保った．この幅は転位論から導出される理論解と一致しており，本モデルで APB エネルギーを正しく反映した超転位を良好に再現できることが示された．

2. γ 相中に 2 本の転位を配置し，せん断応力 τ =300MPa, または 400MPa 与えて前方の γ 0 相をカッティングさせるシミュレーションを行った．τ =300MPa では，1本目の転位 (転位 1) は界面で停止し，2 本目の転位 (転位 2) もその反力を受け静止し，γ 0 相へ侵入することはなかった．一方，τ =400MPa では，転位 1 は界面で停止するものの，転位 2 が約 72b まで接近すると γ 0 相に侵入した．転位 2 は界
面で反力を受けることなく侵入し，γ 0 相内で転位 1，2 は一定の距離 (約 30b) を保ち超転位として運動した．

3. 2 のせん断応力 τ =300MPa のモデルにおいて，同一すべり面上にもう 1 本転位を配置したシミュレーションを行った．3 本目の転位 (転位 3) が界面で停滞している転位 1，転位 2 に接近すると，転位 1，2 間の距離が縮められ距離が約 58b になったとき γ 0 相へ進入した．その後は 2 と同様に転位 1，2 は，両転位間一定の距離 (約 30b) を保ち超転位として運動した．転位 3 は，界面に到達すると APBのエネルギーバリアによって停止した．

51

第5章

NEB 法による反応経路解析
転位の Phase Field 法は転位論 (静止した転位に対する) に基づいてモデル化されているため，phase field が変化し転位が運動しているときの正当性は現時点では証明できない．すなわち，外力を受けると結晶学的エネルギー，弾性ひずみエネルギー，そして勾配エネルギーのバランスによって phase field が変化し転位が移動するが，移動中の「非平衡状態」のエネルギーを真に考慮したものではない．APB を例に挙げると，

L12 構造中で「積層欠陥」に相当する APB は，1 本の転位が通過する前に比べるとエネルギーが高いが，実際に原子面を [110] 方向に移動させる途中にはより高いエネル
ギーバリアが存在する．したがって，転位の移動については「終状態」の APB エネ
ルギーではなく，途中の経路のエネルギー障壁を求めなければならない．転位が γ/γ 0
界面を貫通する場合も同様である．そこで本章では，NEB 法により APB 形成ならび

γ/γ 0 界面を貫通する転位のエネルギーバリアを求める「反応経路解析」を行う．

5.1

APB 形成の反応経路解析

5.1.1

解析条件

Ni3 Al 相の (111) 面を xy 平面とし，z 軸方向に積み重ねた直方体セル (2.02nm×1.74nm
×4.92nm) をシミュレーションセルとして用いた（図 5.1）．前章のように z 方向に周期境界を適用すると，2 つの APB が同時に生成するような場合となり NEB 解析が困難となる．そこで x，y 方向は周期境界条件，z 方向は自由境界条件とした．この直方体セルを下段と上段に分け，上段の原子を [1̄10] 方向に 1 原子分ずらすことで APB を導入する．この構造を最急降下法を用いて緩和させて最終状態とした． また，ずらす
前の構造についても同じセルで最急降下法で構造緩和を行い初期状態とした．原子数

52

nm

2.02
2.46nmNi
Al

2.46nmAPB
view point

1.74

nm

(a) Model(b) View point

Fig.5.1 Simulation cell for NEB analysis.

は 1,536 である．この初期状態と最終状態を基に NEB 法で計算を行った．NEB 法におけるイメージ（位相）点数は 20 とした．

53

5.1.2

解析結果及び考察

図 5.2 は初期状態と最終状態を結ぶ最小エネルギー経路に沿ったエネルギーバリアのプロファイルである．エネルギーの値は初期状態のエネルギーを基準とし, そこからの差として表している．横軸は反応「座標」と表示しているが，初期状態と最終状態の間の位相点の番号を表しており，実空間における座標とは無関係である．図で 1 〜

3 で示した点ならびに初期，最終状態について，APB を導入するすべり面上下の原子配置を図 5.3 に示す．黒色で着色した原子は上面にある原子，水色で着色した原子は下面にある原子であり，すべり面法線方向から見ている．分子動力学計算ではないの
で，各スナップショットは時系列に沿った連続したものではなく，NEB 計算によって最も低いエネルギー経路を通る「途中の」構造として得られたものである．fcc では直接 [110] への原子移動を生じるのではなく，図 5.4 に模式的に示すようなショックレー
の部分転位として生じることがよく知られているが，途中の原子配置はそれに対応している．特にエネルギーバリアの 1 ， 3 では，下側の原子間に一直線上に位置するときであり，途中エネルギーの谷となる 2 は fcc の積層欠陥位置に対応する．左側の極大点では 371.8mJ/m2 ，右側の極大点では 479.7mJ/m2 であり，初期状態に近い側より最終状態に近い側のエネルギーが高い値を示した．

500③

Planar fault energy [mJ/m2]

①
400
300

②

200

end
APB energy (0K)

100
0

start
0

5

10

15

20

Reaction coordinate

Fig.5.2 Reaction path for APB formation.

54

[112]
y[110]
x① Reaction coordinate : 8Start② Reaction coordinate : 9 ③ Reaction coordinate : 13

End

Fig.5.3 Snapshots at each reaction coordinate.

BBa/6[211]

C

A

B

B

A

B

A

A

A

A

B

C

a/6[121]AA

a/2[110]

B

A

B

A

C
A

Fig.5.4 Schematic of Shockley’s partial dislocation.

55

5.2

γ/γ 0 界面を貫通する転位の反応経路解析

5.2.1

解析条件

図 5.5 に示すように，Ni 相と Ni3 Al 相が (010) 整合界面を有するようなシミュレーショ
ンセルに刃状転位を導入して解析を行った．セルの寸法は 34.02nm×12.30nm×1.74nm
であり，x，y ，z 方向の結晶方位はそれぞれ [1̄10]，[111]，[1̄1̄2] である．刃状転位は，
図 5.6 に模式的に示したように，端から 2.52nm の幅 3.53nm の領域においてすべり面
上下の原子列の数を変えることで導入する．原子数は 63,600 個である．y 方向 (転位線方向) は周期境界条件を，x，z 方向は最外周部からカットオフ距離にある原子を表面の法線方向には固定する条件で，分子動力学シミュレーションによる 10000fs の初期緩和を行った．このとき．y 方向の応力が 0 となるようにセル寸法を制御している．熱揺動の影響を排除するため温度は 1K とし，速度スケーリングにより制御した．初期緩和後の系に，さらに最急降下法を用いて NEB 解析の初期状態 (図 5.7(a)) とする．図はCommon Neighbor Analysis(CNA)(42) により fcc 以外と判断された原子 (転位芯) のみ表示しており，黄色で着色した部分は leading partial と trailing partial，灰色で着色した部分は局所的には hcp と判定される積層欠陥である (図下の拡大図)．ただし，上の全体図では x，z 方向の端面原子も黄色で表示している．先の図 5.6 のように上下の原子面数を変えたことで，初期緩和シミュレーション中に転位は右に移動し，界面に到33.26nm12.30nm

Ni
Al

[111]
z
[110]
x

[112] y[110]x1.74nminterfaceγ'phaseγ phase16.63nmFig.5.5 Simulation cell.

56

16.63nm達して停止している．一方，NEB 法の最終状態として，界面に止められているこの転位が γ 0 相に侵入した状態を求めるために，初期緩和後の系全体にせん断ひずみ zx を与える分子動力学シミュレーションを温度 1K で行った．せん断は，系の上半分の原子には [110] 方向に，下半分には [1̄10] 方向に毎ステップ dux だけの微小変位を与えることで行った．ここでdux = ζ {tan(γ + dγ) − tan(γ)}(5.1)

である (図 5.8)．dux は変位量，ζ は中央のすべり面からの座標 (距離)，γ はせん断ひずみである．ひずみ速度 dγ は一定とし，dγ = 1.0 × 10−5 1 / fs = 1.0 × 108 1 / s とした．せん断ひずみ γ=0.04 まで与え刃状転位が γ 0 相に侵入した後，せん断ひずみ γ を除した状態で 10000fs の緩和を行い，さらに最急降下法を用いて最終状態 (図 5.7(b))
を作成した．初期状態と最終状態の間のイメージ（位相）点数は 11 として NEB 解析を行った．

2.52nm 3.53nm

Ni
Al

A'

A

[111]
z
[110]
x

la2×14
la2

A'

A

la1
la1×15

Fig.5.6 Introduction of edge dislocation.

57

γ'phase

γphase

[111]
z
[110]
x

Interface

Trailing partial
[111]
z
[110]
[112] yx

Leading partial(a) Initial state(b) Final stateFig.5.7 (a) Intial state and (b) Final state for NEB.

dγ

dux

ζ

γ

[111]
z
[110]
x
Fig.5.8 Variable displacement for constant shear.

58

5.2.2

解析結果及び考察

図 5.9 に初期状態と最終状態を結んだ最小エネルギー経路を示す．前節同様，縦軸のエネルギーは初期状態を基準として表示している．エネルギーは反応座標 10 でピー
クを示し，そのときのエネルギーバリアは 0.80 × 10−3 eV/Å である．先の完全結晶から APB を生じるときには大きなエネルギーバリアが存在し，APB を形成するとエネ
ルギーが大きく低下したが，完全結晶ではすべり面上の全原子がいっせいに移動しなければならないためである．一方，γ/γ 0 界面をカッティングする際のエネルギーバリ
アのピークと，最終状態とのエネルギー差は非常に小さく，0.041 × 10−3 eV/Å 程度である．転位が界面を通過する前後では，界面に段差を生じるものの，全体的な構造緩和は生じないためピークの後のエネルギー低下は小さいものと考える．図 5.9 の各位相点に対応する転位構造を図 5.10 に示す．反応座標 4 において leading partial が γ 相内に侵入し，反応座標 10 の前で trailing partial が侵入し始める．なお，初期状態での転位の幅は 3.1nm であるが，反応座標 10 までの変化で転位の幅 (4.2nm) が広がっている．それ以降急激に転位の幅は狭まり，最終状態での転位の幅は 1.7nm となる．エネ
ルギーバリアのピーク後のわずかなエネルギー減少は転位の幅の変化に対応する．

Energy barrier [eV/A]( 10-3 )
0.8

end

0.6

0.4

0.2

start

0
0

5

10

Reaction coordinate
Fig.5.9 Reaction energy of unit length.

59

γ'phase

InterfaceγphaseInitial state

Reaction coordinate : 4

Reaction coordinate : 8

Reaction coordinate : 10

Leading partial

Trailing partial

Reaction coordinate : 12

Final stateFig.5.10 Reaction images of edge dislocation cutting through the γ/γ 0 interface.
60

5.3

結言

APB エネルギーそのものでなく，APB を形成する途中のエネルギーバリアを NEB法により評価した．また Ni/Ni3 Al の (010) 整合界面に停止した fcc の刃状転位が，Ni3 Al相に侵入して leading partial となるときのエネルギーバリアも評価した．得られた結果を要約して以下に示す．

1. Ni3 Al の完全単結晶を初期状態，APB を導入した Ni3 Al 単結晶を最終状態として，イメージ (位相) 点数 20 で NEB 解析を行った．各位相点におけるすべり面上下の原子配置をみたところ，最小エネルギー経路は fcc の部分転位の経路を通っており，最初の部分転位のバーガースベクトルに相当する原子移動のエネルギー
バリアは 371.8mJ/m2 ，後続の部分転位のそれは 479.7mJ/m2 と評価された．これは最終的に APB を含んだ状態より非常に高い．

2. Ni 相と Ni3 Al 相の (010) 整合界面を作成し，刃状転位を導入して界面に接近・貫通した状態を分子動力学シミュレーションにより求めた．それらを初期状態，最終状態として NEB 解析を行ったところ，leading partial と trailing partial がほぼ Ni3 Al 相に侵入した点でエネルギーの極大値 0.80 × 10−3 eV/Å をとった．

3. 2 の極大点から最終状態までのエネルギー減少はわずかであり，かつそれが Ni3 Al
相侵入時に広がった拡張転位の幅が狭まることによるものであることを明らかに
した．

61

第6章
結論
本研究は，Ni 基単結晶超合金中の超転位をモデル化することを目的として，分子静/動力学による逆位相境界 (APB) エネルギーの評価，Phase Field 法によるシミュレー
ション，NEB 法による APB 形成および γ/γ 0 界面を転位が貫通するときのエネルギー
バリアの評価等を行った．以下に，得られた結果を総括する．
第 2 章では，本研究で用いた解析手法の基礎について概説した．はじめに，分子動力学法の基礎方程式や，原子間相互作用の評価に用いた原子埋め込み法ポテンシャルの概要を説明し，分子動力学計算における温度制御法と大規模分子動力学計算に用いられる領域分割による高速化について述べた．次に，転位の Phase Field 法で用いる
秩序変数および固有ひずみの説明をし，TDGL 方程式で使用する各自由エネルギー汎関数について述べた．また，応力やひずみを求めるため有限要素法の概要も説明した．
最後に NEB 法における位相空間内での最小エネルギー経路探索の概念を述べた．
第 3 章では，熱膨張を想定し，平均格子間隔を広げて静力学的に解析した時の APB
エネルギーと，実際に熱揺動の効果を考慮したときのそれの違いを議論するために，

EAM ポテンシャルによる静力学解析ならびに動力学解析を行った．熱膨張を想定して体積 (格子長さ) を増加させると，(111) 面での APB エネルギーは 0K の格子長さにおける値 180.4mJ/m2 からなだらかに単調減少した．一方，10〜1000K の温度条件下
で分子動力学解析を行ったところ，APB エネルギーはいずれも 145.6±2.2mJ/m2 となり，静力学解析では APB エネルギーを過大に評価してしまうことがわかった．その差を調べるため，APB 近傍の原子層毎にエネルギーを調べた．静力学解析では APB の
上下 1，2 原子層のエネルギーが非常に高く，格子長さを広げてもこの分布形状は変わらずにわずかにその値が減少するのに対し，動力学解析では APB の上下 1 層目のエネ
ルギーは 10K と 1000K でほぼ等しく，2〜4 層目のエネルギー分布が高温ほど広がりかつ高くなる．しかしながら，温度上昇とともに系の体積が増加するので，この APB 近62

傍でのエネルギー上昇は相殺され，APB エネルギーは温度によらずほぼ一定となる．
第 4 章では，これまでに提案されてきた転位の Phase Field モデルを参考に，γ 0 相の

APB を考慮した超転位のシミュレーションモデルを提案した．提案したモデルで，まず γ 0 相 (L12 構造) 単相中の 2 本の転位のシミュレーションを行い，与えた APB エネル
ギーの値に応じて，転位論から導かれる超転位の幅を再現できることを示した．次に，

γ 相 (fcc) 中に 2 本の転位を配置し，せん断応力 τ =300MPa または 400MPa を与えて前方の γ 0 相をカッティングさせるシミュレーションを行った．せん断応力 τ =300MPa では，転位は γ/γ 0 相界面に停滞し γ 0 相に侵入しなかった．一方，せん断応力 τ =400MPa
では，界面に停止した 1 本目の転位が後続の転位からの反力を受けて γ 0 相に侵入した．
この場合，後続の転位は界面で反力を受けることなく侵入し，両者は γ 0 相内で一定の距離 (約 30b) を保ち超転位として運動した．さらに，せん断応力 τ =300MPa のシミュ
レーションにおいて，同一すべり面上にもう 1 本転位を配置した場合も検討し，離散
転位動力学シミュレーションでも報告されているように 3 本目の転位からの反力を受けて低いせん断応力で超転位を形成することを確認した．
第 5 章では，APB を形成する途中のエネルギー障壁，および，Ni/Ni3 Al の (010) 整
合界面に停止した fcc の刃状転位が，Ni3 Al 相に侵入して超転位の leading partial となるときのエネルギー障壁を NEB 法により評価した．まず，Ni3 Al の完全単結晶を初期状態，APB を導入した Ni3 Al 単結晶を最終状態とし NEB 解析を行った．最小エネル
ギー経路での各位相点におけるすべり面上下の原子配置をみたところ，fcc におけるショックレーの部分転位と同じ原子移動であること，そのときの最初と後続の部分転位に相当するエネルギーバリアはそれぞれ 371.8mJ/m2 ，479.7mJ/m2 であること，などが示された．次に，Ni 相と Ni3 Al 相の (010) 整合界面を有するシミュレーションセ
ルに刃状転位を導入し，Ni/Ni3 Al 界面 (γ/γ 0 界面) に停滞した転位が γ 0 相に侵入する際のエネルギーバリアを，極低温 (1K) の分子動力学シミュレーション＋ NEB 解析により求めた．最小エネルギー経路は単一ピークをもつエネルギーバリアとなり，その値は 0.80 × 10−3 eV/Å であった．先の完全結晶からの APB 形成と異なり，この途中のエネルギーバリアと転位が γ 0 相に侵入した最終状態のエネルギーにはほとんど差がな
かった．
以上，分子静/動力学，Phase Field 法，NEB 法など様々な計算力学手法を用いて Ni
基単結晶超合金中の超転位に関する検討を行った．本研究で得た知見は，転位貫通による γ 0 相の界面形態変化などのより複雑な現象のモデル化に反映させることで，ラフト化における転位の役割解明の一助となると期待される．

63

参考文献
(1) 原田広史, Ni 基超合金開発の現状と新材料の可能性, 日本ガスタービン学会誌, 284(2000), 278-284.
(2) Yokokawa, T., Osawa, M., Nishida, K., Kobayashi, T., Koizumi, Y. and Harada,
H., Partitioning behavior of Platinum group metals on the γ and γ 0 phases of
Ni-based superalloys at high temperatures, Scr. Mater.，49 (2003)，1041–1046．
(3) Walston, W. S., O’Hara, K. S., Ross, E. W., Pollock, T. M. and Murphy, W. H.,
Third Generation Single Crystal Superalloy, Superalloys.， (1996), 27-34.
(4) Zhang, J. K., Murakumo, T., Koizumi, Y., Kobayashi, T., Hrada, H. and Masaki,Jr., S., Interfacial Dislocation Networks Strengthening a Fourth-Generation SingleCrystal TMS-138 Superalloy , Metall.Mater.T rans A., 33 (2002), 3741-3746.
(5) 横川忠晴，原田広史，大野勝美，高温 X 線回折による Ni 基超合金の γ/γ 0 格子定数ミスフィットの測定技術，溶接学会誌，67-2 (1998)，34-37．

(6) 小泉裕，張建新，小林敏治，横川忠晴，原田広史，青木祥宏，荒井幹也，白金族元素ルテニウムを含む次世代 Ni 基単結晶超合金の開発，日本金属学会誌，67-9(2003)，468-471．
(7) Tien, J. K. and Copley, S. M., The eﬀect of orientation and sense of applied
uniaxial stress on the morphology of coherent gamma prime precipitates in stressannealed nickel-base superalloy crystals, Phys. Stat. Sol., A, Metall. Trans., 2-2(1971), 543-553.
(8) Grose, D. A. and Ansell, G. S., The influence of coherency strain on the elevatedtemperature tensile behavior of Ni-15Cr-Al-Ti-Mo alloys, Metall. Trans., A12(1981), 1631-1645.
(9) Li, J. and Wahi, R. P., Investigaton of γ/γ 0 lattice misfit in the polycrystalline
nickel–base superalloy IN738LC: influence of heat treatment and creep deformation, Acta Metall. Mater., 43 (1995), 507–517.

64

(10) 石橋浩一, 近藤義宏, 行方二郎, 大井成人, 服部博, Ni 基単結晶超合金 CMSX-2 における γ 0 相のラフト化の引張方位依存性, 耐熱金属材料第 123 委員会研究報告,34–2(1993), 165–172.
(11) N. Zhou, C. Shen, M. J. Mills, Y. Wang, Acta Materialia 55 (2007) 5369-5381(12) Yashiro, K., Pangestu, J. R., Tomita, Y., Structure and Motion of Misfit Dislocation at Ni/Ni3 Al Interface : Molecular Dynamics Study, Journal of T he Societyof M aterials Science, Japan (2006), submitted.
(13) Yashiro, K., Pangestu, J. R., Tomita, Y., Molecular Dynamics Study of Interfacial Dislocation Network at Gamma/Gamma-Prime Interface in Ni-Based Superalloys, T hird International Conf erence on M ultiscale M aterials M odelingM M M 2006 F reiburg, (2006), 582–584.
(14) 鈴木雄風, 屋代如月，冨田佳宏，刃状およびらせん転位とミスフィット転位の転位
芯相互作用：分子動力学による解析 ，日本機械学会論文集，A 編，73–735 (2007),

1217–1224.
(15) Zbib, H. M. and Rubia, T. D., A multiscale model if plastucity, Int. J. Mech.
Sci.，18 (2002)，1133–1163．
(16) 都留智仁, 渋谷陽二, 三次元離散転位力学法と境界要素法のマルチスケール解析とインデンテーション問題への応用日本機械学会論文集, A 編, 74–743 (2008),

933–938.
(17) 小林亮 , 齋藤幸夫, 2001 日本結晶成長学会誌 28 (4) 265-266(18) B. Nestler, A. Adam, Wheeler Computer Physics Communications 147 (2002)
230-233(19) 上原拓也 , 辻野貴洋 , 井上達雄, (2003) 材料力学部門講演会講演論文集 497-498(20) M. Ortiz, (1999) Journal of Applied Mechanics, Transactions ASME 66 (2) 289298(21) N. Ghoniem, S. Tong, L. Sun, (2000) Physical Review B 61 (2) 913-927(22) Y. Wang, Y. Jin, A. Cuitiño, A. Khachaturyan, Acta Materialia, 49 (2001),18471857.
65

(23) C. Shen, Y. Wang, Acta Materialia 52 (2004) 683-691(24) S. Y. Hu, Y. L. Li, Y. X. Zheng, L. Q. Chen, International Journal of Plasticity
20 (2004) 403-425(25) 屋代如月，内藤正登，冨田佳宏，Ni 基超合金格子状析出物に堆積する刃状転位の分子動力学シミュレーション，日本機械学会論文集，A 編，69–677 (2003), .210-217(26) K. Yashiro, F. Kurose, Y. Nakashima, K. Kubo, Y. Tomita, H.M. Zbib, International Journal of Plasticity 22 (2006) 713-723(27) 久保圭佑, 神戸大学修士論文, Ni/Ni3 Al 界面の第一原理分子動力学法による力学特性評価, (2006)(28) T. Zhu, J. Li, A. Samanta, H. G. Kim and S. Suresh, Proceedings of the National
Academy of Sciences, 104, (2007), 3031-3036(29) T. Zhu, J. Li, A. Samanta, A. Leach and K. Gall, Physical Review Letters, 100,
(2008), 025502(30) Daw, M. S. and Baskes, M. I., Semiempirical, Quantum Mechanical Calculationof Hydrogen Embrittlement in Metals, Phys. Rev. Lett., 50–17 (1983), 1285–1288.
(31) Daw, M. S. and Baskes, M. I., Embedded–Atom Method : Derivation and Application to Impurities, Surfaces and Other Defects in Metals, Phys. Rev. B, 29–12(1984), 6443–6453.
(32) Voter, A. F. and Chen, S. P., Accurate Interatomic Potentials for Ni, Al andNi3 Al, Mat. Res. Soc. Symp. Proc., 82 (1987), 175–180.
(33) Chen, S. P., Voter, A. F., Albers, R. C., Boring, A. M. and Hay, P. J., Investigation of the Eﬀects of Boron on Ni3Al Grain Boundaries by Atomistic Simulations,J. Mater. Res., 5–5 (1990), 955–970.
(34) Voter, A. F., Chapter 4 The Embedded Atom Method, John Wiley & Sons,(1994), 77–90.
(35) Rose, J. H., Smith, J. R., Guinea, F. and Ferrante, J., Universal Features of theEquation of State of Metals, Phys. Rev. B, 29–6 (1984), 2963–2969.

66

(36) Foiles, S. M., Calculation of the Surface Segregation of Ni-Cu Alloys with theUse of the Embedded-Atom Method, Phys. Rev. B, 32–12 (1985), 7685–7693.
(37) 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
(38) Zhang,J.X., Murakumo,T., Koizumi,Y., Kobayashi., Harada., Acta Materialia,
51 (2003),3741(39) Karnthaler,H.P., Muhlbacher,E.TH. and Rentenberger,C., Acta Materialia, 44(1996),547.
(40) Sun,J., Lee,C.S., Lai,J.K.L. and Wu,J.S., Intermetallics, 7 (1999),1329(41) K. Yashiro, M. Naito, Y. Tomita, Int. J. Mech. Sci. 44, (2002), 1845-1860(42) Honeycutt, J. D. and Andersen, H. C., Molecular Dynamics Study of Melting andFreezing of Small Lennard–Jones Clusters, J. Phys. Chem., 91 (1987), 4950–4963.
.

67

付録 A
講演論文
 蟹川 淳 ・高木 知弘・屋代 如月，Ni 基単結晶超合金における超転位の Phase
Field シミュレーション，第 14 回分子動力学シンポジウム，愛媛県県民文化会館(2009.5)．
 蟹川 淳 ・高木 知弘・屋代 如月，EAM ポテンシャルによる Ni3Al 中の APB
エネルギー評価と PhaseField シミュレーションによる超転位のモデリング，日本
機械学会関西支部第 85 期定時総会講演会，神戸大学六甲台キャンパス (2010.3)．

 蟹川 淳・屋代如月，原子シミュレーションによるＡＰＢエネルギー評価とＰｈａｓｅ Ｆｉｅｌｄシミュレーションによる超転位のモデル化，日本機械学会第 23 回計算力学講演会，北見工業大学 (2010.9)．

68

69

70

71

72

73

謝辞
本研究を遂行するにあたり，懇切丁寧なご指導賜りました田中克志教授に深甚な感
謝の意を表します．また，阪上隆英教授の有益な御教示や福井工業大学工学部機械工学科冨田佳宏教授の広い御見識と暖かな御人柄に導かれたことは，著者にとって無類の幸福でありました．
浅学非才の著者が本論文を完成させるまで研究全般に渡り直接御指導を頂いたのみ
ならず，日頃の研究生活に対して幅広い視点から御助言・御協力頂きました屋代如月
准教授に深く感謝致します．研究を進めるにあたり，丁寧な御指導と有意義な助言を頂きました京都工芸繊維大学大学院工芸科学研究科高木知弘准教授に心より深く感謝
申し上げます．また有益な議論を賜りました長谷部忠司准教授ならびに木之下博助教
授に深く感謝致します．研究活動が円滑に行えるように研究環境の整備に御尽力して頂ました古宇田由夫技官に厚くお礼を申し上げます．本研究を遂行する上で，数々の御支援と御協力を頂きました山中晃徳先輩 (現 東京工業大学 大学院理工学研究科 助教) を初めとする諸先輩方に深く感謝致します．また，日々の生活において，励まし合
い，時には刺激し合った同期の相賀裕太郎氏，芦田雅樹氏，金谷敬輔氏，坂田了介氏，
岸本和也氏，常見祐介氏，初田祐貴氏，原田将伍氏，福田晃司氏，森田泰博氏，睦門
賢憲氏，村上智宣氏，望月利紀氏に感謝致します．
最後に，六年間の大学生活に対して精神的にも経済的にも多大なる支援と理解頂いた父母に感謝と敬意の意を表しましてここに記します．

平成 23 年 2 月

74

蟹川淳