修士

論文

原子弾性剛性係数による

Al の粒界およびき裂の局所格子不安定性解析

指導教員： 田中

津川

克志

悠太

2015 年 2 月

神戸大学大学院 工学研究科 博士課程前期課程

機械工学科専攻

要約
本研究では原子弾性剛性係数 (AES) によりナノレベルの変形・破壊現象を統一的に議論する研究の一環として，対称ねじれ粒界および，貫通き裂を有する Al の引張シ
ミュレーションを行い，系の応力-ひずみ挙動と原子弾性剛性係数 Bijα の行列式 detBijαの関係について検討した．まず，様々な Σ 値の対称ねじれ粒界の引張シミュレーショ
ンを行った．Σ3,19,17 では無負荷での粒界エネルギーの大小によって弾性限界を議論できたが，Σ9 と Σ11 では逆の結果となった．Σ3,Σ11, Σ17 では引張ひずみが増加す
ると粒界近傍にて負の AES 原子が増加していき，粒界から転位が発生して応力が急減
し，Σ9 および Σ19 粒界は応力-ピーク直前にて粒内の原子がほとんどすべて負になった後に応力急減した．次に，貫通き裂を有する結晶方位の異なる 3 つの系の引張シミュ
レーションを行い，AES による局所不安定性解析をモード I 型き裂先端に適用した．き裂があるにもかかわらず応力- ひずみ関係に非線形性が現れ始めるまで負の AES 原子は出現しなかった．[1̄1̄2](111) き裂では先端の対称性が崩れた後に先端から転位が射出された．[001](010) き裂ではき裂先端が丸みを帯びた直方体形状に鈍化した後に左右の対角線上角部から，[1̄10](001) き裂ではき裂先端の真の底部に負の AES 原子が出現
し，その後転位射出が確認された．

Summary
In order to clarify the relationship between the unstabe behavior and the atomic
elastic stiﬀness(AES), we performed molecular dynamics simulations on twist grainboundaries(G.B.) and mode I crack tips in fcc Al. First, we performed tensile simulation on twist G.B.s of several Σ value. Σ17 shows the smallest strain for the
stress-strain peak, and has the highest G.B. energy in the simulated grain boundaries.
Σ3 and Σ11 also reach the elastic limit in the order of the G.B. energy. The grainboundary of Σ19 and Σ9 show much higher strain for the stress-strain peak than thethree Σ value above mentioned. Although the Σ19 has the lowest G.B. energy, however the Σ9 has the higher G.B. energy than Σ11. That is, the magnitude relationsof the elastic limit of Σ11 and Σ19 contradict to the G.B. energy. On the change in
the AES, negative AES atoms gradually increase around the Σ3, Σ11 and Σ17 G.B.s
and dislocation emission is observed before the stress-strain peak. On the other hand,negative AES atoms hardly increase just before the stress peak in the Σ9 and Σ19
G.B.. Then negative AES atoms built in the bulk part of the grains. After negative
AES atoms increased rapidly within a grain, Σ9 and Σ19 grain boundary reached theelastic limit.
Next, we performed molecular dynamics simulation on the [1̄1̄2](111), [001](010) and[1̄10](001) through cracks under mode I loading. There is no negative AES atoms in the
no load equilibrium. The negative AES atoms energies at the point where the stressstrain curve begins to show nonlinearity. In the [1̄1̄2](111) crack case, the symmetryof the crack tip is broken and the negative AES atoms appear at the shoulders of theopened crack tip, then dislocations are emitted. The crack tip of the [001](010) crack issome what square, and the negative AES atoms and dislocation emission occur in the
opposite direction of the diagonal line. The [1̄10](001) crack also shows the emergenceof negative AES atoms and following dislocation emission, only in the one-side tip.

目次
第1章 緒 論1

第 2 章 解析手法の基礎

42.1

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

4

2.2

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

5

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

6

2.2.1
2.3

高速化手法

2.4

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

92.5

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

10

2.5.1

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

10

2.5.2

原子系への拡張 . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.5.3

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

11

2.5.4

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

12

2.5.5

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

13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 3 章 ねじれ粒界の引張シミュレーションへの適用

7

14

3.1

解析モデル及び解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3.2

解析結果および考察

16. . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 4 章 モード I 型き裂先端の AES 評価

23

4.1

解析モデル及び解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . .

23

4.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25第5章 結 論

38

参考文献

40

i

目次

ii

学術論文・学術講演

43

謝辞

46

第1章
緒論
分子動力学法は材料の特性を個々の原子の運動にまで遡って解析する手法であり，材
料科学だけでなく，高分子科学や生体化学など様々な分野で活用されている. 近年の計算機の性能の向上，およびナノレベルでの材料の変形・破壊実験が可能になってきた今日では，ナノ〜サブミクロン領域の現象を対象とする分子動力学シミュレーショ
ンは一般に用いられるようになってきた．しかしながら，分子動力学シミュレーショ
ンの中で見出した局所の変形メカニズムについて，どのような指標を用いて普遍的・定量的な議論に結びつけるかということは依然として問題である．通常はポテンシャ
ルエネルギーや局所構造，局所応力やひずみ等によって議論されるが，原子間ポテン
シャルを用いたシミュレーションは仮想材料の結果であるため，臨界エネルギーや臨界ひずみ等のしきい値を導出することはあまり意味を持たないように思われる．したがって使用する原子間ポテンシャルによらない普遍的な基準を模索する必要があると考える．
そこで著者らのグループでは局所変形抵抗を表す物理量として原子弾性剛性係数(Atomic Elastic Stiﬀness，AES) に着目し，その正値性に基づいて局所変形の開始クライテリアを統一的に議論することを試みてきた．AES による検討は，fcc，アモルファ
ス金属，そしてダイヤモンド構造のシリコンについて進められてきた [1]–[4] ．fcc 構造のNi についての検討では，完全結晶に近いナノワイヤの引張で表面から部分転位が発生する際に AES が負の原子が転位芯の位置に発生すること，アモルファス構造 (Ni-Al)
での検討では，AES の正値性によって局所クラスタの崩壊及び再構築が表されること，1

またシリコンの検討では，表面 (内部不均一を有する系) の解析において弾性限界まで
はひずみ増加による変形状態が系全体で均一化すること，傾角粒界を有した bcc 鉄の検討では負の detBijα 原子が前駆的に欠陥原子が示された範囲に生じることなどを報告している．
結晶粒界は二次元格子欠陥の一種であり，多結晶材料の特性を大きく左右する要因である．一般的に，Hall-Petch 則によって多結晶体を構成する結晶粒の大きさを微細化することで，強度の向上を示すことが知られる．しかしこのことはマクロスケールで適用される話であり，結晶粒がナノサイズになると，粒界が占める割合が上昇し，従来材料にない性質を示すようになり，実際に ECPA[5]–[9] などの強ひずみ加工により結晶粒をナノメートルとした材料では，強度が飛躍的に上昇すること，強度と靱性がともに向上することが報告されている [8] ．一方で，結晶粒界が原因で本来延性を示す材
料が脆性的に破壊することが報告されている [10] ．そのような粒界に起因した破壊メカ
ニズムの解明のため，高温でのアルミニウムの粒界拡散の解析 [11], [12] や粒界溝の変形解析などの原子シミュレーションが広く行われてきた [13]–[15] ．また粒界構造や粒界エ
ネルギーおよび粒界破壊強度に関するシミュレーションも銅やモリブデン，α-鉄，マグネシウム，およびパラジウムと幅広く行われている [16]–[24] ．
また，材料内に存在するき裂は材料の疲労破壊の原因として重要である. 金属が一度では到底破断しないような小さな負荷であっても，その負荷が繰り返しかかることによってき裂が生じる．き裂先端部分は応力が局所的に増大し，き裂の進展が始まる．
そうすると，き裂は急速に成長していき短時間のうちに材料を破壊する．金属材料内
のき裂伝ぱ観察や評価は鉄 [25] ，ステンレス鋼 [26] ，マグネシウム合金 [27] など種々の金属材料で行われてきた．鉄による進展する疲労亀裂の直接観察によると，疲労亀裂進
展速度と成長モードの間には密接な関係があることが報告されている [28] ．モード 型
繰り返し負荷を受けるアルミニウム合金の検討によると，アルミニウム合金のモード
型疲労亀裂進展速度はモード

型進展速度に比べて大きいことを報告している [29] ．

き裂進展実験だけでなく数値シミュレーションによるアプローチも行われており，有限要素法による銅のき裂成長評価 [30] ，分子動力学法によるアモルファス金属のき裂成
長シミュレーション [31] など幅広く行われている．モンテカルロ法によるアルミニウ

2ム合金のき裂伝ぱ過程の検討によると，微小き裂が結晶粒界と相互作用しながら伝ぱする挙動をよく表し，き裂伝ぱの解析に有効であることを示した [32] ．
本研究では fcc 金属の Al をモデル材料とし，対称粒界やき裂を導入した系に対して引張を与え，非弾性変形時の AES について明らかにすることを目的とする．第 2 章では分子動力学法の基礎理論を述べ，使用する原子間ポテンシャル，不安定性解析に利用する原子弾性剛性係数について説明する．第 3 章では先行研究 [33] をもとに，Σ値の異なるいくつかの対称ねじれ粒界について分子動力学法により引張シミュレーション
を行い，粒界固有の強度について AES の視点から検討する．第 4 章ではき裂について

AES による局不安定解析を Al のモード I 型き裂先端に適用した結果を報告する．

3

第2章
解析手法の基礎2.1

分子動力学法

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

(2.1)

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

Fα = −

∂Etot
∂r α

(2.2)

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

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

(2.3)
(2.4)

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

dr α= v α (t)
dt
4

(2.5)

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

α

α

(2.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)

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

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

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

5

2.2.1

EAM ポテンシャル

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

Etot =N∑F (ρα ) +

α

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

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

(2.13)

β(̸=α)

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

6

2.3

高速化手法

原子数 N の系において粒子間の全相互作用を評価すると，1step に N × (N − 1) 回の計算が必要となり，N が大きくなると極めて膨大な計算量となる．実際には，一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ半径 rc )内の粒子からの寄与を効率よく計算することにより高速化できる．従来よく用いられてきた高速化手法に粒子登録法がある．これは，図 2.1 に示したように，rc よりひとま
わり大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) に計算負荷が減少される．しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計算負荷により打ち消される．

rc

rfc

Fig.2.1 Schematic of bookkeeping method.

7

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

y

by
0

x

bx

Fig.2.2 Schematic of domain decomposition method.

8

2.4

速度スケーリング法

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

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

mα : 粒子αの質量
viα : 温度 T での粒子αの速度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) の右辺を現在の速度に掛けてやればよい．ただ，これだけでは原子配置に反映されないので，Verlet 法における ∆r αi (t+∆t)√（式 2.18）を

T0 /T ∆r αi (t + ∆t) と置き換える必要がある．

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

F αi (t)(2.18)mα

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

9

2.5

原子弾性剛性係数2.5.1

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

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

1
σij =
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

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

Bijkl ≡

∂σij
∂εkl

)

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

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

10

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

ギーの 1 次，2 次変化量として導出される．

2.5.3

原子応力

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

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

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

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

ηij =

]
1[
Gij − δij
2(2.24)であり，その微小量

1dηij = dGij211(2.25)

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

r̄αβ r̄αβ∂rαβ= i αβj∂ηijr(2.26)

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

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

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

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

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

σij =

N
1 ∑
σα
N α ij(2.28)

となる．

2.5.4

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

1 ∂ 2 Etot
V ∂ηij ∂ηkl

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

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

′′

α

′

αβ

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

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

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

α

{ ∑
N

′

αβ

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

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

2.5.5

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

13

第3章
ねじれ粒界の引張シミュレーション
への適用
3.1

解析モデル及び解析条件

図 3.1 に示すように，x，y，z 軸をそれぞれ fcc の [1̄10]，[001]，[110] 方向とする座標軸系において z 軸中心に上下のセルを ±θ/2 ずつ回転させて回転角度 θ の対称ねじれ粒界を作成した．回転角度 θ は表 3.1 に示す 5 つの Σ 値に対応する角度とした．各モ
デルに対して, 全方向周期境界条件下で初期応力が 0 となるようにセル寸法を微調整しながら制御温度 0.1[K] で 50000fs の初期緩和シミュレーションを行った．その後 z 方向に毎ステップひずみ増分 ∆εzz =1.0×10−6 を与える引張シュミレーションを行った．この時，∆εzz =2.0×10−4 毎に原子座標データを記録し，各原子位置の原子応力 σiα ならびに原子弾性係数 Cijα から原子弾性剛性係数 Bijα (i，j は 1〜6 の Voigt 表記) を求め，その行列式 detBijα を評価した．detBijα の値は 0[K] 完全結晶における値で無次元化している．以下では無次元化した detBijα を AES と称している．なお横方向（x，y 軸方向）の応力は 0 となるようにセル寸法をスケーリングしている．

14

lz

z[110]

x[110]y[001]

lylx

Fig.3.1 Schematic of simulation cell of twist grain boundary．

Table 3.1 Simulation conditions for twist grain boundary
Sigma value Misorientation angle
Σ370.53 °
Σ938.94 °
Σ11
50.48 °
Σ17
86.63 °
Σ19
26.53 °

Cell length[nm] Number of atoms
9.8 × 9.9 × 18.3
107520
9.7 × 9.4 × 18.3
101376
9.4 × 9.4 × 18.9101640
9.4 × 10.0 × 18.9
10771210.5 × 9.9 × 18.9
12038415

3.2

解析結果および考察

表 3.2 に緩和計算後の初期平衡状態における AES の平均，最大および最小値，負のAES を持つ原子数，そして粒界エネルギーを示す．粒界から 3 原子層程度離れるとバルクと変わらなくなるため，セル全体で平均した AES には差がなくなる．一方，粒界に存在する負の AES 原子数はΣ値によって差を生じている．いずれの粒界でも完全結晶の値である１より大きな原子が存在しているが，粒界で局所的に密度が高くなっている分布に対応する．特に Σ3 ねじれ粒界は AES の最大値と最小値の幅が大きい．

Table 3.2 After relaxation.
Sigma value Average of
detBijα
Σ3
0.91
Σ9
0.92
Σ11
0.93
Σ17
0.92
Σ19
0.93

Maximum Minimum
detBijαdetBijα
3.25
-0.28
1.68
-0.05
1.41
-0.23
1.81
-0.03
1.89
-0.06

Number of
detBijα < 0 atoms
1496
949
1815
674
566

GB energy[J/nm2 ]
1.63
1.61
1.57
1.70
1.48

図 3.2 に緩和直後の各モデルの粒界面のスナップショットを，AES の正負で着色して示す．一方，図 3.3 は回転して接合した上部を取り除いて表示し，図 3.2 における負の AES 原子のみを 3 次元的に鳥瞰したものである．Σ3 および Σ11 の負の AES 原子の分布はやや不規則だが，表 3.2 の負の AES 原子数からもわかるように多数存在し粒界全面を覆っている．一方，他のモデルでは負の AES 原子の分布は規則的な幾何学形状
で分布している．

16

detBijα <0 atom.

detBijα >0 atom.

(a) ∑3.

(b) ∑9.

(c) ∑11.

yx(d) ∑17.

(e) ∑19.

Fig.3.2 Atomic conﬁgration on grain boundary plane after relaxation．

(a) ∑3.

(b) ∑9.

(d) ∑17.

(c) ∑11.

(e) ∑19.

Fig.3.3 3D bird view of negative AES atoms in Fig.3.2．

17図 3.4 に引張時の応力ひずみ曲線および detBijα の平均値の変化，図 3.5 に負の AES
原子数変化を示す．横方向の結晶方位に差があるものの，引張方向が [110] であるためすべての応力-ひずみ曲線は弾性限界=系の急激な応力減少までは大きな差がない．弾性限界が最も低いひずみとなったのは緑色で示した Σ17 であり，これは表 3.2 で粒界エネルギーが 1.70J/nm2 と最も高かった粒界である．次に黒色の Σ3，青色の Σ11 が弾性限界に達するが，これも粒界エネルギーの順（1.63J/nm2 ，1.57J/nm2 ）となっている．一方，赤色の Σ9，橙色の Σ19 は高い弾性限界を示しており，応力-ひずみの勾配が0 に近い，すなわち系全体の不安定に近いところまで耐えている．これは図 3.5 の (b)の AES の平均が 0 に近い値をとっていることからもわかる．Σ19 は粒界エネルギーが最も低いが Σ9 は 1.61J/nm2 と Σ11 より高く，Σ11 と Σ9 の弾性限界は粒界エネルギーの大小関係と矛盾する．図 3.5 の負の AES 原子数変化に示したように，弾性限界が低ひずみ側の 3 つの粒界は矢印で示した応力ピークの前に負の AES が増加しており (図3.5(a)), ピークひずみより前に局所的な構造変化を生じている可能性が示唆される．一方，高ひずみを示した Σ9 と Σ19(図 3.5(b)) では，先の 3 つの系に比べて負の AES 原
子数はピーク直前までほとんど変化せず，その後急激に増加してピーク応力を示し急減するためパルス状の変化を示している．

18

1

Σ3
Σ9
Σ11
Σ17
Σ19

2

Average detBαij

Stress, σzz, Gpa

2.5

1.5

Σ3
Σ9
Σ11
Σ17
Σ19

0.8
0.6
0.4
0.2

1

0
0.02

0.04

0.06

Strain, εzz

0.08

0.02

0.04

0.06

Strain, εzz

0.08(b) Change in average of detBijα .

(a) Stress-strain curves.

Number of detBαij<0 atoms

Number of detBαij<0 atoms

Fig.3.4 Stress-strain curves and change in the average of detBijα .

Σ3
Σ11
Σ17

5000
4000
3000
2000
1000

0.02

0.04

0.06

0.08

60000

40000

Σ9
Σ19

20000

0
0.02

0.04

0.06

Strain, εzz

Strain, εzz(b) Σ9 and Σ19Bijα twist GBs.

(a) Σ3,Σ11 and Σ17 twist GBs.

Fig.3.5 Change in number of detBijα < 0 atoms.

19

0.08

Σ3 の応力ピーク時 (εzz =0.0606) とその直前の各ひずみのスナップショットを図 3.6に示す．図では左側が atomeye[42] の central symmetry parameter により「欠陥」と判別された原子を，右側には detBijα < 0 となった原子のみを示している．Σ3 の図を見ると応力ピーク直前のひずみ εzz =0.0598 近傍から，粒界面から転位が粒内へと成長していることが atomeye の可視化でわかるが，それに先立って負の AES 原子が伝ぱしていることがわかる．Σ9 粒界における欠陥原子と負の AES 原子の分布を，負の AES が急増し始めたひずみ近傍 (εzz =0.0712〜0.0718) と応力ピーク (εzz =0.0728) 前後について示したものが図 3.7 である．Σ9 では粒内に負の AES 原子が広がっており，そして系全体がほぼ負となった時に粒界面から欠陥原子が発生し始めている．Σ11 と Σ17 粒界は

Σ3 粒界と同じように粒内には負の AES 原子は発生せずに応力ピークを迎え，その後粒界から転位を発生した．Σ19 粒界は Σ9 粒界と同じように粒内に負の AES 原子が急増した後に弾性限界に達した．

20

(a) εzz =0.0596(b) εzz =0.0598(c) εzz =0.0600(d) εzz =0.0602(e) εzz =0.0602(f) εzz =0.0606

zx
Fig.3.6 Distribution of defect and negative AES atoms on Σ3 twist grainboundary. Left ; defect atoms visualized by central symmetry parameter. Right ; negative AES atoms．

21

(a) εzz =0.0712(b) εzz =0.0718(c) εzz =0.0726(d) εzz =0.0728

zx

Fig.3.7 Distribution of defect and negative AES atoms on Σ9 twist grainboundary. Left ; defect atoms visualized by central symmetry parameter. Right ; negative AES atoms．

22

第4章
モード I 型き裂先端の AES 評価
4.1

解析モデル及び解析条件

図 4.1 左側に示すように，平板状の周期セル中央に貫通き裂を有する系を解析対象とする．き裂幅 2c はセル寸法の Lx の 1/5 とした．図 4.1(a)〜(c) の 3 つの結晶方位を対象とし，[1̄1̄2](111)，[001](010)，[1̄10](001) のき裂の解析を行った．き裂変位 h は Voterの EAM ポテンシャルのカットオフ半径 0.555nm より大きく，かつき裂先端半径 ρ がグリフィスき裂の条件 8a/π（a は Al の格子長さ 0.405nm）を満たすように ρ=0.5nm，

h=2ρ=1.0nm としている．各モデルの原子総数は 102992，1056655，104283 である．
各モデルに対して，全方向周期境界条件下において初期応力が 0 となるようにセル寸法を微調整しながら，制御温度 0.1[K] で 100000fs の初期緩和シミュレーションを行った．
その後，き裂垂直方向に毎ステップひずみ増分 ∆εyy =1.0×10−6 を与える引張シュミ
レーションを行った．温度は緩和計算と同じく 0.1[K] の極低温とした．∆εyy =2.0×10−4毎に原子座標データを記録し，各原子位置の原子応力 σiα ならびに原子弾性係数 Cijα から原子弾性剛性係数 Bijα (i，j は 1〜6 の Voigt 表記) を求め，その行列式 detBijα を評価した．detBijα の値は 0[K] 完全結晶における値で無次元化している．

23

y[111]z[112]
Lyh(a)[112](111)crack

2c=0.2Lx

Lx

x[110]

y[010]Lz
z[001](b)[001](010)crack

y[001]z[110]

x[100](c)[110](001)crack

Fig.4.1 Simulation models for mode I through crack．

24

x[110]

4.2

解析結果

表 4.1 に緩和計算後の初期平衡状態における系全体の AES の平均，最大および最小値，負の AES を持つ原子数を示す．き裂形状に違いはなく結晶方位の差のみのため，無負荷では AES の値に大きな差はない．いずれのモデルにおいても完全結晶の値である 1 よりも大きな値を示す原子が存在する．また，無負荷平衡状態では負の AES を持
つ原子は存在しない．

Table 4.1 After relaxation.
crack model

Average of detBijαMax． detBijαMin． detBijα[1̄1̄2](111)[001](010)[1̄10](001)

0.96
0.97
0.97

1.50
1.59
1.55

0.08
0.08
0.11

Number of
detBijα < 0 atoms
000

図 4.2 に，無負荷平衡状態における各き裂周辺の原子を detBijα の値によって色分けしたスナップショットを示す．カラーバーのスケールは最小値の青色が 0.070, 最高値の赤色が 1.606 としている．いずれのき裂モデルもき裂表面の原子はバルクより低い

detBijα の値を示し，0 に近い値をとっている．一方，き裂上下面では 1 原子層内部の表面原子から内部原子は表面から内部に収縮するため，バルクよりも高い detBijα の値を示している．き裂先端部では表面原子は 0 に近い値をとるが，結晶方位と表面構造の違いから一層内側の原子の detBijα はやや高くなっているもののバルクよりは高くなっ
ていない．

25

α

detBij
1.606

0.070(a)[112](111)crack(b)[001](010)crack(c)[110](001)crack

Fig.4.2 Distribution of detBijα after relaxation．

26

図 4.3 に引張時の応力-ひずみ曲線 (a) および負の AES 原子数変化 (b) を示す．いずれのモデルもき裂があるにもかかわらず，引張を開始してしばらくは負の AES 原子が現れなかった．ひずみ εyy =0.03 近傍から負の AES を持つ原子が現れ始める．初めて負の AES 原子が現れた点を応力-ひずみ曲線中に実線矢印で示した．いずれも非線形性が現れる点に対応している．一方，破線矢印は応力ピークとみられるひずみを示している．負の AES 原子は単調に増加し，応力ピーク以降も増加している．ピークひずみが最も小さくなったのは，引張方向が fcc のすべり面に垂直となる [1̄1̄2](111) き裂で

Number of detBαij<0 atoms

あった．[001](010)，[1̄10](001) き裂はピークひずみがほぼ同じで差がない．

Stress, σyy, Gpa

32[112](111)[001](010)[110](001)

1

0

0.05

0.1

0.15

Strain, εyy

0.2[112](111)[001](010)[110](001)

2000

1000

0
0

0.05

0.1

0.15

Strain, εyy(b) Change in number of detBijα < 0 atoms.

(a) Stress-strain curves.

Fig.4.3 Stress-strain curves and change in number of detBijα < 0 atoms.

27

0.2

図 4.4 に [1̄1̄2](111) き裂の負の AES 原子が出現するひずみ εyy =0.0270 までのき裂先端のスナップショットを示す．εyy =0.0270 まではき裂上下表面が拡がり，両端の「<」「>」の角度が開口するが原子配置の組み換えはない．そして先端の対称性が崩れ右上のすべり面方向に負の AES 原子が現れる．図 4.5 に [1̄1̄2](111) き裂の εyy =0.0270 以降の転α位射出の様子を示す．左側には σyy
によって色分けした図を，真ん中には atomeye[42] のcentral symmetry parameter により「欠陥」と判別された原子を，右側には detBijα < 0となった原子を示している．右側の図ではき裂の空隙をグレースケールで塗りつぶして位置を示している．応力-ひずみ曲線に非線形性が現れる εyy =0.0270 以降では，図(b)(c) に示すように central symmetry parameter で欠陥と判別された原子が現れ，その後図 (d)→(e) のように転位が射出される．図 4.6 は射出された転位を central symmetry

parameter で追跡したもので，下方の 2 つは周期境界下で下から下から運動してきた
最初に射出された転位である．この転位は (251) および (21̄5) 面上を運動している．

εyy =0.054 の図 4.6(c) では最初の転位が周期境界下で互いに合体したところであり，この後応力が急減する．すなわち図 4.3 の破線矢印で示したピークひずみは周期境界とすべり面交線で決定づけられるため物理的重要性はない．これらの転位の移動の際には，転位前縁の原子の AES が負になることを動画で確認している．その後もき裂先端から転位射出が繰り返し生じるため，図 4.3(b) に示したように系全体の負の AES 原子数は単調に増加している．

28

(a) εyy =0.000.

(b) εyy =0.004.

(c) εyy =0.008.

(d) εyy =0.012.

(e) εyy =0.016.

(f) εyy =0.020.

(g) εyy =0.024.

(h) εyy =0.026.

(i) εyy =0.027.

Fig.4.4 Snapshots of [1̄1̄2](111)crack tip before the emergence of negative
AES atoms．

29

(a)εyy =0.0270.

(b)εyy =0.0272.

(c)εyy =0.0284.

(d)εyy =0.0292.

y[111]
x[110](e)εyy =0.0300.
α. CenterFig.4.5 Snapshots of [1̄1̄2](111)crack．Left : stress distribution of σyy: defect atoms visualized by central symmetry parameter. Right :negative AES atoms．
30

y[111]x[110](a)εyy=0.046(b)εyy=0.050(c)εyy=0.054

Fig.4.6 Snapshots of defect atoms in [1̄1̄2](111)crack at the strain εyy =0.046〜0.054, visualized by central symmetry parameter．

31

図 4.7 に [001](010) き裂における転位射出の様子を示す．左側には atomeye[42] の cen-tral symmetry parameter による「欠陥」，右側には detBijα < 0 となった原子を緑色に着色して示している．このき裂ではひずみ εyy =0.0392 において初めて負の AES 原子
がき裂先端に出現した．その後，ひずみ εyy =0.0420 においてき裂先端からの転位射出が観察される．射出された転位が移動する際に，転位に先行して detBijα < 0 となる原子が出現する点は先の [1̄1̄2](111) き裂と同様である．[001](010) き裂も引張ひずみを増加し続けるとき裂先端から転位が繰り返し射出されるが，図 4.7(f) のように対角線上の角度から互いに反対方向に射出された転位が周期境界条件下で衝突して対消滅していた．このためピーク後に急激な応力低下を生じなかったものと考えられる．図 4.8 は[001](010) き裂のひずみ εyy =0.0420 までのき裂先端の拡大図である．原子配置に変化がないまま，き裂上下表面が拡がっている．負の AES 原子はエッジ部の原子配置にすき間が多くなった部分に現れている．

32

(a) εyy =0.0420.

(b) εyy =0.0428.

(c) εyy =0.0436.

(d) εyy =0.0444.

(e) εyy =0.0452.

(f) εyy =0.0460.

y[010]
x[100]

Fig.4.7 Snapshots of [001](010)crack．Left : defect atoms visualized by central symmetry parameter. Right : distribution of negative AES
atoms．

33

(a) εyy =0.000.

(b) εyy =0.006.

(c) εyy =0.012.

(d) εyy =0.018.

(e) εyy =0.024.

(f) εyy =0.030.

(g) εyy =0.036.

(h) εyy =0.039.

(i) εyy =0.042.

y[010]
x[100]

Fig.4.8 Snapshots of [001](010)crack tip

34

under tension．

[1̄10](001) き裂では，ひずみ εyy =0.0324 に最初にき裂両端に負の AES 原子が出現
したが，ひずみ εyy =0.037 まで同じ原子位置において負の AES 原子が消滅と出現を繰り返した．図 4.9 は [1̄10](001) き裂における転位射出の様子を示す．これまで同様，
左側は atomeye[42] の central symmetry parameter，右側が負の AES 原子位置である．

[1̄1̄2](111) き裂，[001](110) き裂ではき裂の両端から対称的に転位が射出されたのに対し，[1̄10](001) き裂では片方のき裂先端から転位が射出された．図 4.10 は転位が射出された側のき裂先端の拡大図であり，detBijα < 0 となった原子を緑色に着色して示している．図 (d) ではき裂先端の前縁で負の AES 原子が集団的に現れ，その後転位が射出されその前縁に負の AES 原子が確認される (図 (e)〜(h))．さらにひずみが増加すると再びき裂先端から転位が射出される (図 (i)〜(l))．

35

(a)εyy =0.036.

(b)εyy =0.038.

(c)εyy =0.040.

(d)εyy =0.042.

y[001]
x[110](e)εyy =0.044.
Fig.4.9 Snapshots of [1̄10](001)crack．Left : defect atoms visualized by central symmetry parameter. Right : negative AES atoms．
36

(a) εyy =0.036.

(b) εyy =0.037.

(c) εyy =0.038.

(d) εyy =0.039.

(e) εyy =0.040.

(f) εyy =0.041.

(g) εyy =0.042.

(h) εyy =0.043.

(i) εyy =0.044.

(j) εyy =0.045.

(k) εyy =0.046.

(l) εyy =0.047.

y[001]
x[110]

Fig.4.10 Snapshots of [1̄10](001)crack tip under tension．

37

第5章
結論
本研究では，原子弾性剛性係数 (AES) によりナノレベルの変形・破壊現象を統一的に議論する研究の一環として，fcc 金属の Al について対称ねじれ粒界，および，モード I 型貫通き裂の分子動力学シミュレーションを行い，系の応力-ひずみ挙動と原子弾性剛性係数 Bijα の行列式 detBijα の関係について検討した．以下に本論文で得られた結果を総括する．
第 2 章では解析手法の概要を述べた．はじめに分子動力学法の基礎方程式を示し，原子間相互作用の評価に用いた EAM ポテンシャルの概要を説明した．さらに各原子の安定性評価の指標となる原子弾性剛性係数を概説した．
第 3 章では様々な Σ 値の対称ねじれ粒界を作成し，無負荷平衡状態および引張変形
時の AES と粒界エネルギー等について議論した．無負荷平衡状態において Σ3 および

Σ11 は負の AES 原子が多く，粒界全面を不規則に覆っていた．一方，Σ9，Σ17 および Σ19 では負の AES 原子は比較的少なく，粒界面では規則的な幾何学形状に分布した．これらの各ねじれ粒界モデルに対して粒界面に垂直方向に引張シミュレーションを行った結果，弾性限界が最も低かったのは粒界エネルギーが最も高かった Σ17 であり，ついで Σ3，Σ11 の順に弾性限界に達したが，これも粒界エネルギーの順であった．
一方，Σ9 と Σ19 は高い弾性限界を示し，系全体の不安定条件，すなわちすべての原子の AES が負に近い状態になるまで応力急減を生じなかった．Σ19 は粒界エネルギーが最も低い粒界であるためこの結果は予測できるが，Σ9 のそれは Σ11 より高く，Σ9の弾性限界は無負荷における粒界エネルギーの大小関係で説明できない結果となった．

38

低い弾性限界を示した Σ3，Σ11，Σ17 では，引張ひずみが増加すると負の AES 原子が粒界近傍で増加していき，粒界から転位が発生して応力が急減した．一方，Σ9 および

Σ19 粒界は引張後期まで負の AES 原子はほとんど増加せず，応力-ピーク直前にて粒内の原子が全体的かつ瞬間的に負になった後に応力急減した．
第 4 章では平板状の周期セル中央に貫通き裂を配置したモデルにより，[1̄1̄2](111)，[1̄10](001)，[001](010) き裂の引張シミュレーションを行ってき裂成長と AES の関係について調べた．無負荷平衡状態ではいずれのき裂においても負の AES を持つ原子は存在しなかった．初めて負の AES が出現したのは応力-ひずみ関係に非線形性が現れ始める点であった．いずれのき裂も，き裂先端に負の AES 原子が出現するまでは両端の「<」「>」部分の角度が開口し，き裂上下面が拡がった．[1̄1̄2](111) き裂では先端の対称性が崩れた後に先端から転位が射出され，射出された転位はすべり面上を運動し，周期境界下で反対の先端から射出された転位と合体しそして応力急減した．[001](010)
き裂ではき裂先端が丸みを帯びた直方体形状に鈍化した後に左右の対角線上角部から，[1̄10](001) き裂ではき裂先端の真の底部に負の AES 原子が出現し，その後転位射出が確認された．

39

参考文献[1] 屋代・冨田，日本機械学会論文集（A 編），67，− 656（2001），678．
[2] 屋代・西村・冨田，材料，54(10)，(2005)，1053.
[3] 屋代・藤原，材料，60(11)，(2011)，968.
[4] Yashiro，K．and Fujihara，M．
，Modell．Simul．Mater．Sci．Eng．，20 − 4，（2012）．

[5] V.M.Segal, V.I.Reznikov, A.E.Drobyshevskiy, and V.I.Kopylov,Russian Metallurgy(English trans.), 1,(1981),115[6] V.M.Segal, V.I.Reznikov, V.I.Kopylov, D.A.Pavlik, V.F.Malyshev,Processes of
Plastic Structure Formation of Metalls, (1994),232[7] V.M.Segal, Mater. Sci. Eng, 17,(1995),157[8] V.M.Segal, Mater. Sci. Eng, A345,(2003),36[9] V.M.Segal, Mater. Sci. Eng, A271,(1991),322[10] G.A.Chadwick and D.A.Smith，Chap，10，（1976），Acamedic Press．
[11] T.Kwork，P.S.Ho and S.Yip，Sci，144，（1984），44．
[12] 北村，梅野，大谷，日本機械学会論文集（A 編），62，（1996），2791．
（1994），2843．
[13] 北村・大谷・山中・屋代，日本機械学会論文集（A 編），60，[14] 岩崎・佐々木・千葉・阿部，日本機械学会論文集（A 編），61，（1995），1376．
[15] 北川・中谷・尾方・磨藤・前川，日本機械学会論文集（A 編），62，（1996），853．

40

[16] K・Morita and H．Nakashima，Mater．Sci．Eng．，A234 − 2367,（1997），
1053 − 1056．
[17] 森田・中島，九州大学大学院総合理工学研究科報告，19，（1997），1 − 8．
[18] 森田・上原・連川・英治，日本金属学会誌，61，（1997），251 − 260．
[19] 池田・森田・中島・阿部，日本金属学会誌，63，（1999），179 − 186．
[20] 中島・竹内，鉄と鋼，86，（2000），357 − 362．
[21] 周致震・黒木康洋・池田賢一・吉田冬樹・中島英治・阿部弘，九州大学大学院総合理工学報告，26，（2004），333 − 336．
（2005），303 − 307．
[22] 周致窪・池田賢一・中島英治，日本金属学会誌，69，[23] N・Takata，T．Mizuguchi，K．Ikeda，and H．Nakashima，Meter．Trans．，45，(2004），2099 − 2105.
[24] 水口・高田・池田・中島英治，日本金属学会誌，69，（2005），1010 − 1015．
[25] 粟谷，片桐，白石，材料，25(269), 145-150, 1976
[26] 小寺沢，森，材料，26(289), 948-954,(1977).
[27] 武藤，松元，白岩，日本金属学会誌，78(10),381-387,(2014).
[28] 城野，宗，後藤，日本機械学會論文集 (A 編)， 51(465), 1453-1456, 1985[29] 大塚，森，大島，津山，材料，29(325), 1042-1048, 1980
[30] 樋口，川原，近藤，栗原，日本造船学会論文集，139, 221-226, 1976
[31] 中谷，中谷，北川，杉山，材料，49(3), 275-281, 2000
[32] 板垣，石塚，黄，日本造船学会誌 ,800, (1996) , 114-118[33] 屋代・池宮，神戸大学大学院修士論文，43-69，(2013).
41

[34] 屋代他，第 17 回分子動力学シンポジウム，2012.
[35] 平，林，日本機械学会論文集，33(245),1-10,1967[36] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
[37] 上田顯, 分子シミュレーション-古典系から量子系手法まで-, (2003), 裳華房.
[38] Y. Mishin and D. Farkas，M. J. Mehl and D. A. Papaconstantopoulos,Phys. Rev.
B ,, 59, (1999)，3393-3407[39] National
ment

InstituteofStandards

Laboratory，Interatomic

andPotentials

Technology，Material
Repository

Measure-Project，(2004)，(http://www.ctcms.nist.gov/potentials/Al-Ni.html)[40] Wallace, D. C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
[41] Wang, J., S.Tip, S.R.Phillpot and Wolf, D., Phys.Rev.lett., 71-25, (1993), 4182.
[42] Ju Li， Atomeye： an eﬃcient atomistic conﬁguration viewer．
，Modeling Simul．
Mater． Sci． Eng．，(2003)，11-173[43] A．F． Voter and S．P． Chen，Mat. Res. Soc. Symp. Proc. 82, (1987) 175.

42

学術論文・学術講演
学術講演

 fcc 金属の刃状・らせん転位の運動と原子弾性剛性係数
津川

悠太，池宮

一繁，屋代

如月

日本材料学会第 18 回分子動力学シンポジウム，東京工業大学 (2013.5)43論文貼り付け用ページ

44

論文貼り付け用ページ45

謝辞
本研究を遂行するにあたり，仔細に渡る懇切丁寧な御指導を賜りました屋代如月教
授 (現 岐阜大学) に，心より深く感謝の意を表します．研究を進めるにあたり，日ごろから貴重な御助言を賜りました田中克志教授に心より感謝いたします．また，日々の研究活動が円滑に行えるよう便宜を図っていただきました片山雷太技術職員に心から御礼申し上げます．
本研究室配属後，若輩な私に沢山のアドバイスを下さり，基礎的な事まで丁寧に教えて下さった表面・界面工学研究室の諸先輩方に厚く御礼申し上げます．特に，池宮一繁様 (現 オークマ株式会社)，藤原正大様 (現 株式会社神戸製鋼所) ，山根尭様（現コマツ），寺田駿様（現マツダ株式会社）には更なる感謝の意を申し上げます．

また，互いに切磋琢磨しながら研究に取り組み，苦楽を共にした同輩の伊藤貴彦君，
小出卓矢君，菱谷卓也君，藤麻成貴君，本馬剛徳君，毛利友宙君，森太一君に御礼申し上げます．
最後に．6 年間の大学生活をあらゆる面から支えて頂いた家族には，深く感謝致します．ありがとうございました．
平成 2 ７年 2 月
津川

46

悠太