修士論文

原子モデルによる鋼の酸化物分散強化機構の検討：
分散形態ならびに転位芯との相互作用

指導教員：屋代 如月

山口 明宏

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

Atomistic Study on Oxide Dispersion
Strengthening Mechanism in Steel:Morphology and Core Interaction with Dislocations

February 2012
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan
Akihiro YAMAGUCHI要約ODS 鋼の強化機構解明の一助とすべく，Y と O，Ti と O を陽に区別せずに平均化した「YO」「TiO2 」粒子として扱う近似を用いて，モンテカルロ法，分子動力学法，NEB法など様々な計算力学手法により分散形態ならびに力学変形特性等を検討した．さらに原子弾性剛性係数 (AES) という，局所の変形抵抗の喪失 (ソフト化) の観点からも変形メカニズムを議論した．まず，モンテカルロ法で酸化物分散形態を調べた．バルク
では，YO 粒子が組成，温度によらず微細クラスターを形成するのに対し，TiO2 粒子
は粗大析出物を形成し，YO 粒子の組成および温度で寸法が変化した．また，YO 粒子
が TiO2 析出物の粗大化を抑制すること，Al は YO 粒子を凝集しにくくすること，原子空孔は YO 粒子の周囲に集中すること，TiO2 粒子の粒界への偏析などを明らかにした．分子動力学法による様々な組成での単結晶バルクのせん断シミュレーションでは，酸化物が微細分散した系が，粗大析出物を有する系に比べて高い弾性限界を示す結果が得られた．後者は，析出物界面で AES が負となる Fe 原子が集団的に発生し，そこが双晶変形の起点となっていたが，前者では，TiO2 の界面だけでなく YO 粒子周辺からも AES が負となる Fe 原子が様々な方向，場所で発生，成長することを明らかにした．ねじれ粒界を導入した系の引張シミュレーションでは，すべりに対する最初の応力ピークは組成により変わらないのに対し，その後の破壊過程は組成により大きく異なり，応力-ひずみ曲線も異なる様相を呈した．いずれも，最初の応力ピーク前後に発生したすべりが系全体に伝播した後，追加発生したすべりによる変形が飽和したとき，粒界のへき開または粒内き裂が発生しており，組成によっては粒界での TiO2 粒子の凝集が粒内き裂発生を促した可能性などが示唆された．らせん転位を導入して析出物
構造中を通過させるシミュレーションでは，粗大析出物周囲で交差すべりによる転位
ループの形成や TiO2 界面での欠陥形成，微細析出物間での転位ループ形成メカニズムなどを明らかにした．さらに，NEB 法によりらせん転位が粗大析出物をカッティング
するときのエネルギー障壁を求め，らせん転位が TiO2 の中央に位置するよりわずかにずれた位置にピークを持つこと，したがって中央に位置するときは極小となるがピークとの差は小さく，らせん転位は TiO2 にトラップされないことなどを明らかにした．

Summary2

For a new insight on the strengthening mechanism in ODS steels, we have performed
Monte Carlo (MC), molecular dynamics (MD) simulations and NEB calculations withthe simple pair potential functions for monatomic “pseudo” atoms of Y2 O3 and TiO2 ,that is, each Y, Ti and O atom isn’t distinguished in these simulations but their parameters are fitted to correctly represent the lattice length and bulk modulus of real
Y2 O3 and TiO2 according to the ab initio DFT results. We also discussed the deformation mechanism with the positive definiteness of the atomic elastic stiﬀness (AES),
which indicates the lose of local resistance against deformation or “local softeningatom”. First, MC simulations are performed to discuss the eﬀect of atom composition, vacancy, temperature, etc., on the oxide dispersion morphology. In general, YO
particles form small clusters with the size of a few atoms, and exist homogeneously inthe simulation cell regardless of the oxide composition and temperature. On the other
hand, TiO2 particles tend to form large clusters of a few nanometer size, changing
its size by temperature and compositions. We also found the following facts; (1) YO
particles prohibit coarsening of TiO2 precipitates, (2) Al atoms also prevent the aggregation of YO particles, (3) vacancies accumulate around YO clusters, and (4) TiO2
particles segregate on grain boundary. Then, we have performed shear simulations
by MD on four oxide compositions obtained by MC. One alloy shows higher elastic
limit than the others, and it has the finest TiO2 oxide. From the direct observationof the deformation by central symmetry parameter, it is revealed that all the alloys
cause the stress drop by twining initiated from the interface of TiO2 oxide, while theinitiation and growth of twining are suppressed in the alloy with finer TiO2 clusters.
AES observation also suggests that large domains of negative AES macleate and growat the coarse TiO2 interface, while the negative AES domains are rather small andhave diﬀerent orientations, not only at the TiO2 interface but also at the YO clusters
in the alloy with the finest TiO2 . Tensile simulations are also performed on the bulk
laminate structure with the Σ5 twist grain boundaries, to discuss the eﬀect of grainboundary segregation. The first stress drop is caused by slip in the grains, so that
there is little diﬀerence in pure Fe and four diﬀerent oxide compositions. On the other
hand, the fracture process and the stress-strain response after the first stress peak
are remarkably changed. From the change in the potential energy of each atoms, we
can find shear band like deformation after the first stress peak; then the pure Fe andalloys show the stress increase again and the high energy bands disappear. It means
the grains reach its deformation limit. Then alloys show inter/intra grain cracking depending on its composition. Finally, we have set a screw dislocation in two alloys andglide it by applying shear strain, resulting in the following observations, e.g., formationof dislocation loop by cross slip around coarse precipitates, generation of defects on
the TiO2 surface without cross slip, and so on. We also evaluated the energy barrier
for a screw dislocation cutting a coarse TiO2 precipitate by NEB analysis. The peak ofenergy barrier are evaluated as 1.189[eV/nm] and 1.222[eV/nm], not just at the centerof TiO2 but slightly deviated positions. Thus the energy barrier shows local minimumof 1.144[eV/nm] at the center position; however, the energy diﬀerence is so small that
we can conclude the screw dislocation would not be trapped by the TiO2 cluster.

目次
1

2

緒論

1

1.1

ODS 鋼の開発状況 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

ODS 鋼の数値シミュレーション . . . . . . . . . . . . . . . . . . . . . .

3

1.3

研究目的と各章の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

解析手法

6

2.1

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

6

2.2

モンテカルロ法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2.1

マルコフ過程と Metropolis-Hasting 法 . . . . . . . . . . . . . . .

8

2.2.2

NPT アンサンブル . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.3

擬似乱数発生器 . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.3

高速化手法

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

11

2.4

原子弾性剛性 (Atomic Elastic Stiﬀness: AES) . . . . . . . . . . . . . .

13

2.4.1

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

13

2.4.2

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

13

2.5

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

16

2.6

ポテンシャルフィッティング . . . . . . . . . . . . . . . . . . . . . . . .

17

i

目次

3

モンテカルロ法による酸化物分散形態の検討

3.1

3.2

3.3
4

22

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

22

3.1.1

バルクにおける検討 . . . . . . . . . . . . . . . . . . . . . . . .

22

3.1.2

結晶粒界における検討 . . . . . . . . . . . . . . . . . . . . . . .

24

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

25

3.2.1

MC シミュレーション中のエネルギー変化 . . . . . . . . . . . .

25

3.2.2

Fe-YO-TiO2 の系における組成の影響 . . . . . . . . . . . . . . .

27

3.2.3

Al 添加の効果 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.2.4

原子空孔の影響 . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

3.2.5

結晶粒界近傍の分散形態 . . . . . . . . . . . . . . . . . . . . . .

33

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

35

単結晶のせん断変形シミュレーション

37

4.1

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

37

4.2

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

38

4.2.1

応力-ひずみ応答と欠陥発生の様相 . . . . . . . . . . . . . . . . .

38

4.2.2

系の弾性剛性係数 . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.2.3

AES の評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

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

45

4.3
5

ii

ねじれ粒界のへき開シミュレーション

47

5.1

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

47

5.2

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

49

5.3

5.2.1

応力-ひずみ応答. . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.2.2

変形破壊の様相 . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

5.2.3

系の弾性剛性係数と AES の評価 . . . . . . . . . . . . . . . . . .

55

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

59

目次

6

らせん転位と析出物の相互作用

6.1

6.2

6.3

iii
62

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

62

6.1.1

分子動力学シミュレーション . . . . . . . . . . . . . . . . . . .

62

6.1.2

NEB 法によるエネルギーバリア評価 . . . . . . . . . . . . . . .

64

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

65

6.2.1

応力-ひずみ応答および弾性剛性係数 . . . . . . . . . . . . . . .

65

6.2.2

転位の挙動 (粗大析出物) . . . . . . . . . . . . . . . . . . . . . .

66

6.2.3

転位の挙動 (微細析出物) . . . . . . . . . . . . . . . . . . . . . .

70

6.2.4

NEB 法によるカッティングバリア評価 . . . . . . . . . . . . . .

72

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

73

結論

75

参考文献

79

7

A

関連学術論文・講演

謝辞

84
117

第1章 緒 論

1.1

ODS 鋼の開発状況

燃料被覆管等に代表される原子炉炉心材料は，高温クリープ強度，耐中性子照射性，ボイドスウェリングに対する形状安定性，冷却剤・減速剤に対する耐食性，製管に対する加工性と経済性が同時に求められる．現行の軽水炉には主にジルカロイが使用されており，通常の金属とは異なり中性子照射に対する脆化を起こさないことから現行では十分な特性を備えていると言える (1) ．しかし，軽水炉の更なる高燃焼度化や高速増殖炉等の次世代炉においては高温クリープ強度等の抜本的な向上が求められている．
このような背景の中，次世代原子炉燃料被覆管用材料として酸化物分散型強化 (Oxide

Dispersion Strengthened: ODS) 鋼の開発が進められてきた．ODS 鋼はメカニカルア
ロイング (Mechanical Alloying: MA) によって製造され，酸化イットリウム (Y2 O3 ) 等
を Fe 母相中にナノメートルオーダーで微細に分散させた鋼である．MA により製造された ODS 鋼は優れた高温クリープ特性を示し (2)−(7) ，中性子照射に対して脆化を起こしにくく (2)(6)−(9) ，ボイドスウェリングに対する形状安定性に優れる (6)(7)(10) ことが報告されている．これら重要な特性は，酸化物ナノクラスターによるものであり，ナノクラスターは可動転位に対して強い抵抗となることが示唆されている (2)(5)(11)−(14) ．
さらに，ナノクラスターは 1000◦ C 前後の熱間押出温度において α-γ 変態の抵抗となり (5)(15) ，ナノクラスターにより高転位密度となっている硬い残留フェライトを生じる (4)−(7)(11)(14)−(20) ．ゆえに，ODS 鋼は硬い残留フェライトと緩和された柔らかいマ
ルテンサイトから成る複合材料のような特性を持つ (14) ．この残留フェライトが優れた高温クリープ強度をもたらしており，その主要な変形モードは比較的弱い結晶粒界や焼きならし時のマルテンサイト変態によって生じるパケット境界でのすべりである

1

と考えられている (14)(20) ．そのため，転位に対する障壁と結晶粒界における応力集中との関係で残留フェライトには最適量が存在し，その量は強力なフェライト系元素である Ti の量に左右される (20) ．また，Ti はナノクラスターの微細化を担う元素であり(11)(17)(24)(25)，この点でも粒界における応力集中とトレードオフの関係を持つ．中性子や重イオンの照射に対してナノクラスターが組織の粗雑化と崩壊に抵抗する (2)(21) とともにボイドをトラップする (10)(22)(23) ことが示唆されており，このことが中性子照射に対する優れた特性をもたらすと考えられる．
これら特性と強化メカニズムは全てナノクラスターが起点となっているため，その構造が大変重要となる．原子プローブ断層撮影等の結果からナノクラスターは中心部
分で Y リッチなコア-シェル構造であるとの報告がある (26) ．また最も細かいナノクラ
スターの組成は Y2 Ti2 O7 であるとする実験観察結果がある (24)(27) ．Fe 母相と Y2 O3 ク
ラスターとの結晶方位に強い相関があることも示唆されている (28) ．また，球面収差補正装置を搭載した走査型透過電子顕微鏡による詳細な高散乱角環状暗視野像によりナノクラスターが岩塩型構造であることが示唆され，多くの欠陥を含んでいるとの報告もある (29) ．
利用に向けた研究もなされており，一般的なステンレス鋼に比べて水への耐食性に優れ (30) ，Cr 量の最適化と Al の添加により耐食性が大きく向上する．Al の添加は耐食性向上に欠かせないが，同時に粒界析出物の粗大化とそれに付随する高温クリープ強度の低下を招く．しかし，Al より酸化物を作りやすい Zr や Hf の添加により粒界析出
物の粗大化を抑え，高温クリープ強度を維持したスーパー ODS 鋼がすでに開発されている (7) ．また，熱間押出の際に生じる被覆管の長手方向に細長い残留フェライトの量を熱処理により最適化することで構造異方性を緩和し，周方向の延性と製管に対する加工性を大きく改善できることが分かっている (5)(19) ．

2

1.2

ODS 鋼の数値シミュレーション

このように ODS 鋼の開発は強化メカニズムの考察を伴って順調な進展を見せている．以上はすべて実験観察によるものであり，強化メカニズムの考察は一般的な金属における知見を援用したものである．これら考察はスーパー ODS 鋼の開発につながったことからも妥当かつ極めて有用であるが，ナノクラスターの詳細な構造や転位のダ
イナミクス，粒界のすべりの形態のようなナノ・メゾスケールの知見において未知の部分を多く残している．このような詳細な検討は実験観察と理論のみの考察では限界に近づいており，転位・原子の詳細なダイナミクスやエネルギーバリアを扱うことのできるシミュレーション (数値実験) による検討が期待されている．
密度汎関数理論に基づき電子状態からエネルギーや応力等を計算する第一原理計算
は何ら実験結果を必要とせずに高い精度で検討できるが，膨大な計算量から数十から数百原子程度の取扱いに限られる．そのため ODS 鋼関連では主として Y,Ti,O などの元素の固溶形態の検討に用いられてきた (31)−(34) ．その結果として酸素原子が添加原子間の強力な結合を与えるとともに Y2 O3 ナノクラスターが成長するのに必要となる空孔の集中を促す等の知見を得ている (32) ．我々の研究グループでも，bcc-Fe 中に Y,Ti,O を添加しての安定性を検討し (35) ，現在はそれら構造において格子不安定性を用いた力学的評価を行っている．また，すでに述べたように転位とナノクラスターの相互作用は高温クリープ強度等に影響する重要な要素であり，離散転位動力学 (Discrete Dislocation

Dynamics: DDD) を用いた検討も行われている (36)−(41) ．DDD は任意の曲線形状を有する転位線を多数の直線セグメントに離散化し，転位のダイナミクスを追跡する手法である．DDD を用いた検討から臨界分解せん断応力は Y2 O3 クラスターと母相との格子ミスマッチに強く依存する (39) 等の知見が得られている．また，照射の影響も議論され，ナノクラスターに転位が堆積し内部応力が上昇することで交差すべりした転位が照射欠陥を越えるとの考察もある (40) ．DDD による検討は熱揺動の影響を考慮でき，多数の転位を取扱うことができるという点で大変有用であるが，基本的に転位論と連続体力学によって定式化され，ナノクラスターと転位芯の相互作用もモデル化を通し3て考慮する必要があり，原子レベルでの新たな知見を得ることは難しい．
第一原理計算と DDD の間にあたるスケールでの検討は原子論的モンテカルロ (Monte

Carlo: MC) 法や分子動力学 (Molecular Daynamics: MD) 法による．MC 法は主として構造探索に用いられ，MD 法は転位のダイナミクスやナノクラスター周辺の逐次変化の追跡に用いられてきた．これらは第一原理計算とは違い，原子間相互作用のモデル化を行う必要があるが，数十万から数千万程度の原子の取扱いが可能である．また，

MD 法は DDD と違って熱揺動を考慮できるほどの十分長い過程を現実の計算時間で追跡することはできないが，原子を直接扱うことができ，連続体力学の枠組から離れた知見を得る可能性を持つ．Alinger らは格子モンテカルロ (Lattice Monte Carlo: LMC)法を用いて Fe-Y-Ti-O 系におけるナノクラスターの構造と周辺組成について調べている (42) ．また Hin らは MC ステップを物理的素過程に対応させることで MC 法でも逐次変化を追跡できる動的モンテカルロ (Kinetic Monte Carlo: KMC) 法を用いて単純格子モデルにおけるナノクラスター形成シミュレーションを行った (43) ．一方，ODS 鋼を対象とした MD 法による検討はほとんど見当たらない．これは ODS 鋼中のナノクラスターを MD 法で取扱う際，ポテンシャル関数でイオン結合，配位結合，金属結合を同時に表現しなければならないという困難が背景にある．そこで先行研究では，第一原理計算で求めた Y2 O3 のエネルギー曲線を，Y と O を陽に区別せず平均化した仮想的な「YO 粒子」によってモデル化し，転位と YO ナノクラスターとの相互作用のMD シミュレーションを行ってきた (44)(45) ．

1.3

研究目的と各章の構成

本研究では，ナノクラスターの分散形態および分散形態が力学特性に与える影響を，先行研究と同様に平均粒子を用いて MC 法と MD 法により検討する．ただし，本研究ではさらに Ti の影響も検討する．本研究の特徴は第一原理計算，MC 法，MD 法を駆使することで何ら実験結果を使わず，欠陥の生成過程の検討や転位のダイナミクス評価にまで昇華させる点にある．その際，原子弾性剛性係数 (Atomic Elastic Stiﬀness: AES)

4

という基準を用いて個々の原子の力学状態を考慮した検討を行うとともに，Nudged

Elastic Band(NEB) 法を用いたエネルギーバリアの定量評価にまで発展させる．
以下に, 各章の概略を示す. 第 2 章では MC 法，MD 法における基礎理論について述
べ，それらに必要な原子間相互作用のモデル化（ポテンシャルパラメータの決定）について説明する．第 3 章では，作成したポテンシャル関数を用いた MC 法による酸化物分散形態の検討について述べる．第 4 章では，MC 法で探索した構造の基礎的な力
学的評価としてせん断変形を議論する．第 5 章では，強度に大きく影響する粒界への酸化物の影響を検討するため Σ5 ねじれ粒界のへき開シミュレーションについて述べる．第 6 章では，MC 法で作成した構造にらせん転位を導入し，転位芯とナノクラス
ターの相互作用について検討する．第 7 章では本研究で得られた結果を総括する．

5

第 2 章 解析手法

2.1

分子動力学法

分子動力学法 (Molecular Dynamics Method: MD) は，系を構成する個々の原子についてニュートンの運動方程式
2 α
αd r
= Fα
m
2
dt

(2.1)

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

Fα = −

∂Etot
∂r α

(2.2)

式 (2.1) の数値積分には，Verlet の方法が簡便でエネルギー等の保存性が良いため分子動力学法ではよく用いられる (46)(47) ．時刻 t±∆t での原子 α の座標 r α (t ± ∆t) をテー
ラー展開すると

r α (t ± ∆t) = r α (t) ± ∆t

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

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

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

(2.4)

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

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

(2.5)

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

(2.6)

設定温度が TC ，式 (2.6) より求めたある時刻の温度が T のとき，速度スケーリング
法では，各原子の速度 v α を√

TC /T 倍し設定温度 TC に近づける．Verlet 法では，∆r α (t + ∆t) = r α (t + ∆t) − r α (t) = r α (t) − r α (t − ∆t) + (∆t)2√をF α (t)m(2.7)

TC /T ∆r α (t + ∆t) で置き換えることに相当する．平衡状態では，能勢の方法など

外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られる
カノニカルアンサンブルに一致することが示されている (48) ．

2.2

モンテカルロ法

統計力学の集団理論に基づいて確率過程の生成により原子・分子の空間配置を作り出す方法は統計物理学においてモンテカルロ法と呼ばれる．単にモンテカルロ法という場合，多くはマルコフ過程の確率理論から定式化される Metropolis-Hasting 法を指す．分子動力学法において長時間平均によって求められる力学量とモンテカルロ法において統計平均（確率過程での平均）によって求められる力学量は一致することが分かっている．しかし，モンテカルロ法は分子動力学法とは異なり確率過程の素過程を選ぶ余地があり，分子動力学法では現実の計算時間では到達不可能な平衡状態の探索に力を発揮する．

7

2.2.1

マルコフ過程と Metropolis-Hasting 法

1 つ手前の状態から現在の状態が確率的に決定される確率過程をマルコフ過程という．対象とするマルコフ過程を時間に依存しない定常であるとすると，状態 i から状(n)

態 j に n 回の試行で遷移する確率 Pij は
(n)

Pij =

∑

(n−1)

Pik

Pkj

(2.8)

k
(1)

で表される．ここで Pkj = Pkj である. 式 (2.8) によって生成される状態の連なりのことをマルコフ連鎖と呼ぶ．マルコフ過程の任意の状態から任意の状態への遷移が起こり得るとして，このときマルコフ過程は既約であるという．また，状態 j から出発し(l)て l の試行回数で初めて状態 j に再帰するときの確率を fj で表す．
∞
∑

(l)

fj = 1(2.9)

l=1が成立するとき状態 j は再帰確実な状態であるという．平均再帰時間を

µj =

∞
∑

(l)

lfj

(2.10)

l=1

で定義する．また，同じ試行回数毎に再帰するときその状態を周期的な状態という. マ
ルコフ過程の全状態が再帰確実で非周期的かつ平均再帰時間が有限である状態のとき
マルコフ過程はエルゴード的と呼ばれる．既約でエルゴード的なマルコフ過程において任意の状態 j で(n)lim Pij = wj > 0n→∞(2.11)となる定常分布 {wj = 1/µj } が存在し，

wj =

∑

wi Pij(2.12)

iが成立する (49) ．式 (2.12) を書き改めて∑

wi Pij =

i̸=j

∑
i̸=j

8

wj Pji(2.13)

を得る．この十分条件としてwi Pij = wj Pji(2.14)

で表される詳細釣合い条件式を得る．
式 (2.14) を用いて任意の確率分布を実現する方法は Metropolis-Hasting 法 (50) と呼ばれる．扱う統計集団が既約でエルゴード的，すなわち任意の 2 つの状態は互いに平均再帰時間が有限で到達可能であり周期的でないとする．詳細釣合い条件式 (2.14) から状態 i から状態 j への遷移確率を
{

Pij =

wj /wi (wj <= wi )
1(wj > wi )(2.15)とすると，確率過程の極限分布として指定した分布 {wj } を十分条件として実現できる．

2.2.2

NPT アンサンブル

Lill らは系の変形によるエネルギーを明示的にポテンシャルエネルギーから分離することで系全体の変形を考慮する圧力一定のモンテカルロ法を定式化した (51) ．Lill らの方法に沿った粒子数 N ，圧力 P ，温度 T 一定の系のモンテカルロ法の概要を以下に示す．
シミュレーションセルを平行六面体として 1 つの頂点を直交座標系の原点とする．平行六面体の基本ベクトルを基底としたセル内座標系 si と直交座標系 ri との変換行列を−1
aij = ri t s−1
j ，変形の基準となる変換行列 aij を特に bij ，変形勾配テンソルを Jij = aik bkj ，Green のひずみテンソルを ηij = 21 (Jki Jkj − δij )，指定する真応力テンソルを σij ，第−12piola-kirchhoﬀ 応力テンソルを τij = |J|(Jikσkl Jjl−1 ) とする．また，粒子 n の位置を
(n)

(n)

ri ，速度を qi

(1)

(N )

で表し drN = dri . . . dri
(n)

(1)

(N )

, dqN = dqi . . . dqi，シミュレーション(n)セルのハミルトニアンを H(ri , qi ) とすると，N, σij , T 一定のアンサンブルの分配関数は1
ΞN σT =
N !h3N
p

∫∫∫



(n)

(n)



H(ri , qi ) − |b|τij ηij  N N
exp −
dq dr dJij
kB T

9(2.16)

で表される．ただし，hp はプランク定数，kB はボルツマン定数である．ポテンシャル(n)エネルギーを U (ri ) とすると式 (2.16) より微小体積要素 drN dJij に系を見出す確率は




(n)





(n)

U (si ) − |b|τij ηij 
U (ri ) − |b|τij ηij  N
exp −
dr dJij = exp −|aij |N dsN dJij (2.17)
kB T
kB T
に比例する．ここで式 (2.17) の右辺は左辺を変数変換して得られるものである．W =U − |b|τij ηij − N kB T ln |aij | を定義し, 式 (2.17) を式 (2.15) に当てはめると，十分条件
として状態 l から状態 m への遷移確率
{

Plm =

1
exp {−(Wm − Wl )/(kB T )}(Wm <= Wl )(Wm > Wl )(2.18)

を得る．

2.2.3

擬似乱数発生器

モンテカルロシミュレーションにおいて使用する一様乱数は周期が長く，高次元に対して一様均等分布することが望ましい．そのような乱数を得る手段として熱振動や原子核分裂などの本質的に等確率で起こるされている自然現象を利用する物理乱数が挙げられる．しかし，物理乱数は高価であったり速度が遅かったりする問題がある．さらに，シミュレーションの再現性を確保するには使用した乱数をすべて記録しておく必要がある．それに対して，一様性と均等分布性を合わせ持つ数列を真の乱数の代わりに使うことがモンテカルロシミュレーションでは多い．このような数列を擬似乱数といい，発生させるアルゴリズムのことを擬似乱数発生器という．
現在，最も普及している擬似乱数発生器の 1 つは Lehmer による線形合同法であり，次の線形漸化式により表される．

xn+1 = axn + b (mod M )(2.19)

この擬似乱数は，fortran やＣ言語の標準乱数として組込み関数で与えられている場合が多く扱いやすい．しかし，法を M とするため発生する擬似乱数の周期は高々M であるという問題点を持つ．また，高次元空間でのプロットに対する欠陥も報告されている (52) ．

10

そこで，本研究では Mersenne twister 法 (53) を用いる．今，x，a が要素が 0 または

1 の w 次元の行ベクトルすると漸化式は次式で与えられる．
xk+n = xk+m ⊕ (xuk |xlk+1 )A,(k = 0, 1, . . .)(2.20)

ここで，⊕ は排他的論理和を表し，n は繰り返しを表す整数，m は 0 <= n となる=m<整数である．(xuk |xlk+1 ) は xk の上位 (w − r)bit と xk+1 の下位 rbit を並べて得られる行ベクトルで 0 <= w − 1 を満たす．また，行列 A は次式で与えられる．
=r<





A=







1
1
..

.








1 


(2.21)

aw−1 aw−2 · · · · · · a0

本研究では n = 624，w = 32，r = 31 の場合である Mt19937(53) を用いる．この場合，

2nw−r − 1 = 219937 − 1 がメルセンヌ素数であることから 219937 − 1 という長大な周期を持つ．さらに，623 次元の均等分布も保証されており，プログラムに実装した際のワーキングスペースが小さいという利点を持つ．(53)(54) また，CPU の SIMD 命令や
並列パイプライン処理を生かした，Mersenne twister 法より生成速度が速く周期が長い SIMD-oriented Fast Mersenne Twister(SFMT) 法 (55) も開発されているが，ソース
コードが Mt19937 より複雑で管理が難しいこと，Mt19937 で十分な性能が得られることを踏まえて今回は導入を見送った．

2.3

高速化手法

原子論に立脚したシミュレーションでは原子間相互作用の計算に大部分の負荷が集中する．2 体ポテンシャルを用いる場合，系の全粒子数を N とすれば系の内部だけでも組み合わせは N (N − 1) 通りで粒子数が増加すると飛躍的に計算コストがかかる．しかし，実際の結晶中では近接原子の遮蔽効果により離れた原子からの相互作用はほとんど及ばない．そのため，相互作用打ち切り (カットオフ) 半径 rc を導入し，その半径11内の原子からの寄与のみを考慮すればよい．そこで本研究では領域分割法を用いている．領域分割法では系を間隔が rc 程度である格子状に分割し，対象とする原子の存在する領域 (図 2.1 着色部) および隣接する領域 (図 2.1 斜線部) の原子の内，距離が rc 以
下であるものを探索する．
また，本研究では分子動力学法においてメモリ共有環境を対象とした並列化を行っている．Message Passing Interface(MPI) を用いたメモリ分散環境を対象とする並列化
は，マシンに搭載する CPU 数に捕らわれないが，分散したメモリ間での通信が必要となり，アルゴリズム構築に手間がかかる．一方，OpenMP を用いたメモリ共有環境を対象とする並列化は通信の必要がなく，逐次プログラムで並列化させる部分を指定するのみで済む．領域分割法により原子間相互作用を計算する際，領域毎に独立に計算することができるため，その部分を並列化した．

y
by
0

x
bx

Fig.2.1 Schematic of domain decomposition.

12

2.4

原子弾性剛性 (Atomic Elastic Stiﬀness: AES)

2.4.1

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

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

1
σij =
V

(

∂F
∂ηij

)

1
, Cijkl =
V

(

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

Bijkl ≡

∂σij
∂εkl

)

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

ここで，δij はクロネッカーのデルタである．すなわち，Bijkl は非線形弾性における応力-ひずみ関係の勾配を表す．Wang, Yip らは，ひずみの対称性を考慮したテンソル
sym
Bijkl
≡ (Bijkl + Blkji )/2 の正値性によって結晶の安定性を評価することを提案していsym
る (57) ．すなわち，変形抵抗の喪失点を不安定の物理的描像としている．Bijklの行列式の正負によって格子不安定性が評価される．

2.4.2

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

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

13

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

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

ここで，V は平衡状態における系の体積であり，下付添字のローマ文字はテンソルのデカルト座標成分を表す．(2.24) 式の微分を求めるため，平衡状態からの仮想的な均一変形を考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によってr α = J r̄ α(2.25)と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子 αと 原子 β の間の距離 rαβ には(rαβ )2 = r̄iαβ Gij r̄jαβ(2.26)なる関係が成立する．ただし，Gij = Jki Jkj である．仮想変形の Lagrange ひずみテン
ソル ηij は

ηij =

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

1dηij = dGij2(2.28)と式 (2.26) の関係から次の関係が得られる．

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

これより 2 体ポテンシャル ϕ(rαβ ) における応力は次式で評価される(

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

)

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

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

σijα =αβ αβN
1 ∑1 ′ αβ ri rj
ϕ (r ) αβ
V /N β(̸=α) 2r

14(2.30)

と定義すると，系の応力は
N
1 ∑
σij =
σijα
N α

(2.31)となる．
原子弾性係数弾性係数も応力と同様に U ≈ Etot の場合にはCijkl =

1 ∂ 2 Etot
V ∂ηij ∂ηkl

(2.32)であるので，平衡状態からの仮想均一変形を考えると 2 体ポテンシャルにおける弾性係数は以下のようになる．
(N ∑αγ αγ )N ∑NN
1 ∑∂rαβ ∂ ∑1 ′ αγ ri rj
Cijkl =ϕ (r ) αγ
2V α β(̸=α) ∂ηkl ∂rαβ α γ(̸=α) 2r[ N ∑} αβ αβ αβ αβN {1 1∑ϕ′ (rαβ ) ri rj rk rl′′ αβ=
ϕ (r ) −V 2 α β(̸=α)
rαβ(rαβ )2](2.33)

応力と同様に，各原子位置における原子弾性係数を以下のように定義する．
[} αβ αβ αβ αβ
N {
1 1 ∑
ϕ′ (rαβ ) ri rj rk rlα′′ αβCijkl =ϕ (r ) −V /N 2 β(̸=α)
rαβ(rαβ )2](2.34)

これより，系の弾性係数はCijkl =N
1 ∑CαN α ijkl(2.35)のように原子弾性係数の平均となる．
原子弾性剛性係数以上で定義した原子応力，弾性係数から，原子弾性剛性係数 (AES) は以下で評価で
きる．
α
α
α
α
δil − 2σijα δkl )/2
δjl + σjk
+ (σilα δjk + σjlα δik + σik= Cijkl
Bijkl(2.36)

α sym
α
α
)/2 を
+ BlkjiWang らによる提案 (57) に従い，Voigt 対称性をもたせた Bijkl≡ (Bijklα symαを原子弾性剛性係数
用いて安定性評価を行う．以降では Bijkl
を Voigt 表記した BIJとして用いる．

15

2.5

Nudged Elastic Band (NEB) 法

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

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

End

Energy stateEnd

Energy state

Start

Start

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

16

F⃗iN EB = F⃗i0 + f (ϕi )(F⃗is − (F⃗is・τ̂∥ )τ̂∥ )(2.38)1f (ϕi ) = (1 + cos(π(cos ϕi )))2(2.39)

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

2.6

ポテンシャルフィッティング

本研究では Y や Ti が酸化物を構成している影響を Y と O，Ti と O を陽に区別せず「YO」「TiO2 」粒子として扱うことで簡便に評価する．
まずはモンテカルロ法における原子間相互作用として次式の Morse Potential を採用
し，そのパラメータを決定する．

ϕij = De [exp {−2A(rij − r0 )} − 2 exp {−A(rij − r0 )}](2.40)

ここで，rij を粒子間距離，De , A, r0 はそれぞれエネルギー，長さの逆数，長さの次元を持つパラメータである．
パラメータ De , A, r0 の決定は，Kresse らにより開発された平面波基底ウルトラソ
フト擬ポテンシャル法に基づく第一原理バンド計算コード VASP(58) (Vienna Ab-initio

Simulation Package) により行った．表 2.1 に第一原理計算での条件を示す．以下フィッティングの手順と詳細を示す．まず，図 2.3(a) に示す hcp 単位格子の 2 つのサイトに

Y および O を配置した単位格子と，図 2.3(b) に示す高温で安定な α-TiO2 (ルチル) の単位格子に対し，第一原理計算により系の全自由エネルギーを格子定数 la を変化させながら求めた．次に，Y，Ti，O がそれぞれ原子 1 個が孤立した状態で持つ全自由エネ

17

ルギーを求め，各原子の個数分を単位格子の持つ全自由エネルギーから差し引くこと
で単位格子あたりの結合エネルギーを求めた (図 2.4，図 2.5 実線)．このエネルギー曲
線を．同じ配置であるが Y と O，Ti と O の区別をなくして YO，TiO2 平均粒子からなる「単元系」として，最小二乗法によって Morse potential にフィッティングし，パラメータ De , A, r0 を決定した．表 2.2 に得られたパラメータを示す．平均粒子からなる系の単位格子あたりのポテンシャルエネルギー (図 2.4 破線) を計算するにあたって，

YO 平均粒子の場合，第 7 近接まで考慮して Morse potential を計算している．一方で，TiO2 平均粒子の場合，その構造の複雑さから考慮する近接原子の数によって変動する．
そこで図 2.5 に示すように第 8 近接と第 9 近接の中間まで考慮したとき第一原理計算の結果がその間に入るようにフィッティングした．また，Al 同士の相互作用も Morsepotential で評価するものとし，ポテンシャルパラメータは岩崎らの文献 (59) を参照した．表 2.3 に引用したパラメータを示す．ただし，エネルギーの単位は原論文から変換して表示している．

Fe 同士の相互作用は次式の Johnson potential(60) を採用した．
ϕij = −C1 (rij − C2 )3 + C3 rij − C4 [eV](2.41)

パラメータ C1 ∼ C4 の値を表 2.4 に示す．このポテンシャルは，bcc 構造をとる α-Feのために作られた経験的 2 体ポテンシャルである．
また，異種間粒子の場合は，それぞれの粒子が同種の組み合わせのときにもつ結合あたりのポテンシャルエネルギーを斥力項と引力項にわけ，それぞれを相乗平均をとることによって異種間の組み合わせの結合がもつポテンシャルエネルギーとした.

Table 2.1 DFT calculation conditons.
YO TiO2
Y
Ti
Number of ions2
61
1
Number of electrons
16
32
3
4
Number of bands
12
28
4
4
Cutoﬀ energy [eV] 296.9 495.0 148.9 226.5
Number of k points
153
63
73
73

18

O
1
6
5
495.0
73

Y

Ti

1.633laO

0.644la
O
la

la

la

la

(b) unit lattice for TiO2 .

(a) hcp unit lattice for YO.

Fig.2.3 Unit cells for calculations of total free energy.
-7.0

5

VASP
Morse potential

VASP
n=8
n=9

0

-7.2
-7.4

Energy [eV]

Energy [eV]

-5

-7.6
-7.8

-10
-15
-20
-25
-30

-8.0

-35
-8.2
0.23

0.24

0.25

0.26

0.27

0.28

0.29

0.30

Lattice parameter la [nm]Fig.2.4 YO potential fitting.

-40

0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56
Lattice parameter la [nm]Fig.2.5 TiO2 potential fitting.

Table 2.2 Morse potential parameters for YO and TiO2 .
De [eV] A [nm−1 ] r0 [nm]YO-YO
0.232029 7.16842 0.364215
TiO2 -TiO2 0.676051 13.2322 0.253578
Table 2.3 Morse potential parameters for Al.
Al-Al

De [eV]
0.1193058

A [nm−1 ]23.53643

r0 [nm]0.2864944

Table 2.4 Johnson potential parameters for Fe-Fe.
rij range [Å ]
C1
C2
C3
C4
1.9-2.4
2.195976 3.097910 2.704060 7.436448
2.4-3.0
0.639230 3.115829 0.477871 1.581570
3.0-3.44
1.115035 3.066403 0.466892 1.547967

19

More potential はパラメータが少なく，合わせ込む際の手間も比較的少ないが，その分，関数の自由度が小さい．そこで，分子動力学法向けに平均粒子においても Fe と同じ Johnson potential のパラメータを決定した．Johnson potential は第 1 次微分までの連続性を考慮すると独立なパラメータの数は 9 個であり，Morse potential より多い．そのため，より詳細に体積弾性率を追従できる可能性を持つ．合わせ込みの手順と VASP により計算したエネルギー曲線は Morse potential の場合と全く同様である．
また，異種粒子間では同種粒子間の２体ポテンシャルエネルギー値の相加平均をとり，最小二乗法で合わせ込むことで Johnson potential パラメータをあらかじめ決定しておいた．得られたパラメータを表 2.5∼2.9 に示す．

Table 2.5 Johnson potential parameters for YO-YO.
rij range [Å ]-2.80
2.80-3.68
3.68-4.29

C1
0.154700.418820.20371

C2
5.212133.73824
2.90755

C3
2.40528
0.81108
1.16981C4
9.54316
3.254104.48026

Table 2.6 Johnson potential parameters for TiO2 -TiO2 .
rij range [Å ]-2.52
2.52-2.95
2.95-3.66

C1
3.24926
2.65600
0.34211C2
2.91738
2.86962
1.71521C3
2.934652.36968
3.88863C4
9.391887.87776
11.7160

Table 2.7 Johnson potential parameters for Fe-YO.
rij range [Å ]-2.38
2.38-3.46
3.46-4.26

C1
1.15840
0.566610.25457

C2
3.25778
3.29399
3.55386

C3
1.687310.427130.38527C4
5.04371
1.69365
1.55162

Table 2.8 Johnson potential parameters for Fe-TiO2 .
rij range [Å ]-2.43
2.43-3.11
3.11-3.64

C1
2.737091.72268
0.12507C2
2.984102.914210.40886

C3
2.69997
1.37733
3.91689C4
8.076294.5922210.0382

Table 2.9 Johnson potential parameters for YO-TiO2 .
rij range [Å ]-2.94
2.94-3.63
3.63-4.30

C1
1.89672
0.47657
0.04912C22.93498
2.92525
1.7241920C31.14839
1.148620.97924

C4
4.15864
4.15929
3.37124

以下，モンテカルロシミュレーションには Morse potential，分子動力学シミュレー
ションには Johnson potential を使用する．参考に，モンテカルロ法，分子動力学法における同種粒子間の 2 体ポテンシャル曲線をそれぞれ図 2.6，図 2.7 に示す．

6

Fe-Fe
YO-YO
TiO2-TiO2
Al-Al

Pair potential energy [eV]

5
4
3
2
1
0
-1
-2
0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

Distance [nm]Fig.2.6 Pair potential energy curves for Monte Carlo simulations.

6

Fe-Fe
YO-YO
TiO2-TiO2

Pair potential energy [eV]

5
4
3
2
1
0
-1
-2
0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

Distance [nm]Fig.2.7 Pair potential energy for molecular dynamics simulations.

21

第 3 章 モンテカルロ法による酸化物
分散形態の検討
3.1

解析条件

3.1.1

バルクにおける検討

定性的な分散形態の検討にあたり，粒子数・体積・温度一定のカノニカルアンサン
ブルを採用する．すなわち，式 (2.18) の遷移確率を用いるが，系全体の変形に関わる探索は行わない．シミュレーションは以下の手順で行う．まず，大局的な添加粒子の安定サイトを探るため，原子の移動を bcc 格子点間での入れ換えに限定する格子モン
テカルロ法を行う．次に，格子モンテカルロシミュレーションに続いて格子点以外への探索を行うモンテカルロシミュレーションを行い，酸化物ならびに空孔の分散形態を検討する．
図 3.1 に示すように立方体シミュレーションセルを解析対象とする．セルの寸法はbcc の単位格子 50 × 50 × 50 と 42 × 42 × 42 の 2 つを用いた．bcc 構造の Fe 完全結晶から，乱数を用いて Fe を YO，TiO2 ，Al に置換することで初期配置を作成する．Fe-

YO-TiO2 の系では 50 × 50 × 50 のセルを用い，YO，TiO2 の組成をそれぞれ表 3.1 に示すように合計 3 × 2 = 6 通りとした．また，ODS 鋼では耐食性の向上のため Al が添加されるため，42 × 42 × 42 のセルに Fe-0.4wt%YO-0.2wt% TiO2 -1.0wt%Al の組成で

YO，TiO2 ，Al を分散させた系でも検討した (表 3.2)．温度はいずれのシミュレーショ
ンでも 1000[K]，1400[K] の 2 通りを考慮し，全方向周期境界条件を課した．
また，Fe-YO-TiO2 の系においては原子空孔を導入した場合も検討した．原子空孔を取り扱うために，いずれの粒子との２体ポテンシャルエネルギーも 0[eV] となる「空孔22

粒子」を新たに定義し，LMC シミュレーションに限定して解析を行った．解析条件を
表 3.3 に示す．bcc の単位格子 50 × 50 × 50 のセルを用い，Fe-0.2wt%YO-0.2wt% TiO2の組成とし，全粒子数 (250,000 個) の 0.5%(1250 個) が原子空孔であるとした．温度は1000[K] としている．

L

L
L

Fig.3.1 Simlation cell for MC simulation in bulk.
Table 3.1 Condition of Monte Carlo simulations for bulk Fe-YO-TiO2 alloys.
Cell structure
50 × 50 × 50 bcc unit lattices
Cell length L [nm]14.35
Number of atoms
250000Temperature [K]1000,1400
YO composition [wt%]0.2, 0.4, 0.6
TiO2 composition [wt%]0.2, 0.4
Boundary conditon
Periodic boundary condition
Table 3.2 Condition of Monte Carlo simulations for bulk Fe-0.4wt%YO0.2wt% TiO2 -1.0wt%Al alloy.
Cell structure
42 × 42 × 42 bcc unit lattices
Cell length L [nm]12.054
Number of atoms
148176
Temperature [K]1000,1400
Boundary conditon Periodic boundary condition
Table 3.3 Condition of Monte Carlo simulations for bulk Fe-0.4wt%YO0.2wt% TiO2 -1.0wt%-vacancy system.
Cell structure
50 × 50 × 50 bcc unit lattices
Cell length L [nm]14.35
Number of atoms
250000Temperature [K]1000
Boundary conditon Periodic boundary condition23

3.1.2

結晶粒界における検討

結晶粒界が酸化物の分散形態にどのような影響を与えるかを検討するため傾角対称
粒界を導入した LMC シミュレーションを行った．シミュレーション条件を表 3.4 に示
す．シミュレーションセルは図 3.2(a) に示す直方体セルとした．bcc-Fe 単位格子を並べた後，図 3.2(b) のように y 方向中央の面の上下で [001] 方向を軸に 18.435◦ ずつ図中黄色の線に対称に回転させ，方位角 36.87◦ を持つ Σ5 の傾角対称粒界を導入した．YO 粒子の組成を 0.2wt%,0.6wt%，TiO2 粒子の組成を 0.2wt%,0.4wt% とし，合計 2 × 2 = 4
通りの場合を解析した．バルクと同様，一様乱数を用いて Fe 原子を添加粒子に置換することで初期配置とする．全粒子数は 200,000 であり，原子空孔も全粒子の 1.0% 導入
した．3.1.1 節と同様にカノニカルアンサンブルを仮定し，格子点サイトにおいて粒子種の交換を行う NVT-LMC 法を全方向周期境界条件，温度 1000[K] の条件下で行った．

Table 3.4 Condition of LMC simulations for Fe-YO-TiO2 alloy with G.B.
Cell Length [nm]9.08 × 45.38 × 5.74
Sigma number
5 (misorientation angle: 36.87 deg.)
Number of atoms
200000Temperature [K]1000
Boundary conditon
Periodic boundary condition[010]
Grain Boundary

[100]
22.69

Grain Boundary
misorientation angle
36.87

45.38

5.74

[100]

y

9.08

y

x
z

[010]

unit of length : [nm]x(a) Simulation cell(b) Tilt grain boundary

Fig.3.2 Simulation cell for Σ5 tilt grain boundary.

24

3.2

解析結果

3.2.1

MC シミュレーション中のエネルギー変化

図 3.3 に 1400[K] における Fe-0.6wt%YO-0.4wt%TiO2 の系のシミュレーション中のポテンシャルエネルギー変化を例として示す．LMC シミュレーションにおける MCstep
を LMCstep と呼び，格子点以外への探索を行うシミュレーションにおける MCstep と区別する．LMC シミュレーションおよび格子点以外への探索を行うシミュレーショ
ンに共通して，初期の数百ステップで急激にポテンシャルエネルギーが低下し，その後，ステップが進むにつれて変化の割合を小さくしながら徐々にポテンシャルエネル
ギーが低下していく．図 3.4 に各 LMC，MC ステップにおけるスナップショットを示す．図 3.4(a),(b),(c) に示すように，106 LMCstep 以降，析出物の全体的な位置関係は
変わっておらず，一部粒子がわずかに変化しているだけである．2 × 106 LMCstep の状態を LMC 法での平衡点として格子点以外への探索 (MC) を行うと，図 3.3(b) に示したようにポテンシャルエネルギーはさらに大きく低下し，2.3 × 105 MCstep 計算してもまだエネルギーがわずかに低下している．しかしながら，図 3.4(d),(e),(f) を見ると，5 × 104 MCstep 以降は粒子が大きく移動しておらず，1.5 × 105 MCstep のときと2.3 × 105 MCstep のときの違いを図 3.5 に示すが，白線で囲ったようなごくわずかな違
いしかない．他のシミュレーションでも同様であり，計算時間との兼ね合いもあるため本研究では 1.5 × 105 MCstep 後の粒子配置を平衡状態とみなし，以降で酸化物の分
散形態および析出物内部の構造を検討する．
-450.0-455.0-455.0Potential energy [keV]

Potential energy [keV]-454.5-455.5-456.0-456.5-457.0-457.5-458.0-460.0-465.0-470.0-475.0-480.0-485.0-490.0-495.0-500.0-458.5
0

5

10

15

-505.0

20

0

5

5

10

15

20

23

4

10 LMCstep

10 MCstep(a) LMC simulation.

(b) Random MC simulation.

Fig.3.3 Change in potential energy of Fe-0.6wt%YO-0.4wt%TiO2 alloy under
1400[K].

25

[010]

[010]

[100]

[010]

[100]

[010]

[100]

(b) At 1.6 × 106 LMCstep.

(a) At 106 LMCstep.

(c) At 2 × 106 LMCstep(0MCstep).

[010]

[100]

(d) At 5 × 104 MCstep.

[010]

[100]

[100]

(e) At 1.5 × 105 MCstep.

(f) At 2.3 × 105 MCstep.

Fig.3.4 Snapshots during LMC and MC trial from initial random structure(Fe-0.6wt%YO-0.4wt%TiO2 alloy under 1400[K]).

[010]

[010]

[100]

[100]

(a) At 1.5 × 105 MCstep.

(b) At 2.3 × 105 MCstep.

Fig.3.5 Slight diﬀerences in the atom morphology between 1.5 × 105 MCstep
and 2.3 × 105 MCstep.

26

3.2.2Fe-YO-TiO2 の系における組成の影響

図 3.6 に温度 1000[K]，図 3.7 に温度 1400[K] での各組成におけるシミュレーション
セル全体の YO 粒子と TiO2 粒子の分布を示す．温度，組成によらず共通する点として，YO 粒子 (図中黄色の粒子) 数個からなるクラスターがシミュレーションセル全体に均等分布している．また，TiO2 粒子 (図中青色の粒子) は YO 粒子を内包する比較的大きい析出物を形成した．組成の違いによる影響としては，同じ温度で TiO2 粒子の組成が同じ場合を比較すると，TiO2 粒子の組成が YO 粒子の組成より小さくなるほ
ど TiO2 析出物が小さく複数に分かれる傾向が見られる．温度の違いによる影響としては，図 3.6(a) と図 3.7(a)，図 3.6(b) と図 3.7(b)，図 3.6(d) と図 3.7(d) をそれぞれ比較す
ると，1000[K] において TiO2 析出物がある程度の大きさで複数個形成されたのに対して，1400[K] においては 1 つの大きな析出物に集約されていた．しかし，Fe-0.6wt%YO-

0.2wt%TiO2 の組成のときでは，比較的大きな TiO2 析出物は 1000[K](図 3.6(e)) に比べて 1400[K](図 3.7(e)) で少なくなっていた．
図 3.8 に温度 1000[K]，図 3.9 に温度 1400[K] におけるシミュレーションセル内の TiO2
析出物近傍を抜き出して示す．いずれの組成，温度の場合も，析出物内部に注目すると，TiO2 粒子が単純立方格子の位置に存在し，Fe 原子または少量の YO 粒子が TiO2粒子の間に存在した．析出物内の YO 粒子は，Fe 母相中に分散している YO 粒子と同様に YO 粒子クラスターを形成していた．
本研究では，異種間のポテンシャルエネルギーを引力項，斥力項に分けて相乗平均で評価している．したがって，図 2.6 から分かるように，TiO2 同種粒子間の引力相互作用は他の全ての組み合わせに比べて大きく，析出物が大きくなる傾向を持つと考えられる．また，YO 同種粒子間の 2 体引力相互作用は他の組み合わせに比べて遠距離まで働く．よって，YO 粒子の組成が TiO2 粒子の組成に比べて大きい場合，TiO2 粒子が大きい析出物を形成する前に，YO 粒子付近に TiO2 粒子が集まったと考えられる．
また，Metropolis 法では温度が高いほどポテンシャルエネルギーの増加する遷移の確率が大きくなる．ゆえに，1000[K] では越えられないエネルギー障壁を 1400[K] では越えることで，1000[K] のときに比べて 1400[K] では TiO2 粒子が大きな析出物としてま27とまったと考えられる．したがって，遠距離の引力相互作用を及ぼす「YO 粒子によるトラップ」と「高温によるエネルギー障壁の乗り越え」のバランスによって，TiO2 が
析出物を形成するか YO 粒子の周囲にトラップするかが分かれることが示唆される．

[001]

[001]

[010]

[010]

[100]

[100]

[001]
[010]
[100]

(a) Fe-0.2wt%YO-0.2wt%TiO2 .

(b) Fe-0.2wt%YO-0.4wt%TiO2 .

(c) Fe-0.4wt%YO-0.2wt%TiO2 .

[001]

[001]

[001]

[010]

[010]

[100]

[100]

(d) Fe-0.4wt%YO-0.4wt%TiO2 .

[010]
[100]

(e) Fe-0.6wt%YO-0.2wt%TiO2 .

(f) Fe-0.6wt%YO-0.4wt%TiO2 .

Fig.3.6 Dispersion morphology of YO and TiO2 particles under 1000[K].

[001]

[001]
[010]

[010]

[100]

[100]

[001]
[010]
[100]

(a) Fe-0.2wt%YO-0.2wt%TiO2 .

(b) Fe-0.2wt%YO-0.4wt%TiO2 .

(c) Fe-0.4wt%YO-0.2wt%TiO2 .

[001]

[001]

[001]

[010]
[100]

(d) Fe-0.4wt%YO-0.4wt%TiO2 .

[010]

[010]

[100]

[100]

(e) Fe-0.6wt%YO-0.2wt%TiO2 .

(f) Fe-0.6wt%YO-0.4wt%TiO2 .

Fig.3.7 Dispersion morphology of YO and TiO2 particles under 1400[K].
28

TiO 2
YO

[010]

Fe

[010]

[100]

[010]

[100]

[100]

(a) Fe-0.2wt%YO-0.2wt%TiO2 .

(b) Fe-0.2wt%YO-0.4wt%TiO2 .

(c) Fe-0.4wt%YO-0.2wt%TiO2 .

[010]

[010]

[010]

[100]

[100]

(d) Fe-0.4wt%YO-0.4wt%TiO2 .

[100](e) Fe-0.6wt%YO-0.2wt%TiO2 .

(f) Fe-0.6wt%YO-0.4wt%TiO2 .

Fig.3.8 Magnified view around a precipitate of YO and TiO2 in the simulation
cell under 1000[K].

TiO 2

YO

Fe

[010]

[010]

[100]

[010]

[100]

[100]

(a) Fe-0.2wt%YO-0.2wt%TiO2 .

(b) Fe-0.2wt%YO-0.4wt%TiO2 .

(c) Fe-0.4wt%YO-0.2wt%TiO2 .

[010]

[010]

[010]

[100]

(d) Fe-0.4wt%YO-0.4wt%TiO2 .

[100]

[100]

(e) Fe-0.6wt%YO-0.2wt%TiO2 .

(f) Fe-0.6wt%YO-0.4wt%TiO2 .

Fig.3.9 Magnified view around a precipitate of YO and TiO2 in the simulation
cell under 1400[K].

29

3.2.3

Al 添加の効果

図 3.10 に 1000[K] における Al 原子を添加した場合のシミュレーション結果を示す．

Al 原子を添加しなかったときと同様に，YO 粒子は原子数個のクラスターを形成，TiO2粒子は比較的大きい析出物を形成している．詳細に観察すると Al 原子を添加した場合は YO 粒子間の距離が広がり，また YO 粒子クラスターから外れて孤立した YO 粒子も多く存在しており，Al 原子の添加により YO 粒子クラスターの形成が妨げられる傾向が確認できる．Al 原子の分布には粗密があるもののシミュレーションセルに全体的に分布し，Al 原子の密な部分には YO 粒子が存在している．同一温度，Al の存在以外
同一組成の図 3.6(c) と図 3.10(a) を比較すると，Al 原子を添加しなかった場合に見られた YO 粒子周辺にトラップされた小さな TiO2 クラスターは Al 原子を添加することでほぼ見られなくなった．また，Al 原子は図 3.10(b) に拡大して示すように，YO 粒子や Fe と同様に TiO2 析出物の内部にも存在している．
図 3.11 に 1400[K] におけるシミュレーション結果を示す．TiO2 粒子の大きさに顕著な変化はないが，1000[K] のときに比べて黄色の YO 粒子が集まる傾向があり YO 粒子クラスターがわずかに粗大化している．また，1000[K] のときに比べて Al 原子も同様の傾向を示している．図 3.11(b) に拡大して示したように，1000[K] のときと同様に

Fe，YO，Al が TiO2 クラスター内部に侵入している．Al 原子が析出物や YO 粒子ク
ラスター，Fe 原子間サイトに分布しているのは，図 2.6 に見られるように Al 原子の結合エネルギー (ポテンシャルエネルギーの極小値) が他の組み合わせに比べて小さいためと考えられる．

30

Fe

YO

TiO2

Al

[001]
[010]

[010]

[100]

[100]
(a) Dispersion forms of YO, TiO2 and Al.

(b) Magnified view.

Fig.3.10 Final morphology of Fe-0.4wt%YO-0.2wt% TiO2 -1.0wt%Al alloy under 1000[K].

Fe

YO

TiO2

Al

[001]
[010]

[010]

[100]

[100]

(a) Dispersion forms of YO, TiO2 and Al.

(b) Magnified view.

Fig.3.11 Final morphology of Fe-0.4wt%YO-0.2wt% TiO2 -1.0wt%Al alloy under 1400[K].

31

3.2.4

原子空孔の影響

図 3.12(a) に原子空孔を導入した場合の LMC シミュレーション結果を示す．図では
空孔サイトを「緑の粒子」として表示している．原子空孔は数十個程度毎に集積し，ある程度の大きさを持ったボイドとなってシミュレーションセルに一様に分布していた．
詳細に観察すると，拡大図 (図 3.12(b)) に示すように原子空孔 (図中緑色の粒子) は YO
粒子 (図中黄色の粒子) が数個程度集まってできたクラスターの周辺に集中している．
より詳細に原子空孔の YO 粒子近傍への集中を示すため，図 3.16 に YO 粒子から n 番
目の格子点サイトにおける各原子の占有率を示す．原子空孔はほぼ全て YO 粒子の最近接サイトに存在しており，第 2 近接以降に原子空孔はほとんど見られなかった．最近接サイトの距離 (0.249[nm]) では，図 2.6 より Fe-Fe および TiO2 -TiO2 はエネルギーが低いのに対し，YO-YO 結合では YO ポテンシャルエネルギーが正となり強い斥力相
互作用を受ける．したがって，YO 粒子の最近接サイトに原子空孔が存在するとポテ
ンシャルエネルギーがより低位になり安定化する．
なお，原子空孔を含まない同じ組成の場合 (図 3.6(a)) と比較すると，YO 粒子の分散形態および TiO2 粒子の析出形態は大きく変化しておらず，この結果からは原子空孔
は酸化物の分散形態には大きく関与しないことが示唆される．
YO

TiO2

Atom Vacancy

[001]

[001]

[010]

[010]

[100]

[100]

(a) Dispersion morphology of YO and TiO2(b) Magnified view.

particles.

Fig.3.12 Result in Fe-0.2wt%YO-0.2wt%TiO2 alloy with 0.5% atom vacancies
under 1000[K].

32

Vacancy

100

TiO2

TiO2

TiO2

Ratio [%]

80

60

YO

40
Fe
20

0

1

2

3

4

5

6

7

8

9

10

the "n"th nearest sites from YO particles

Fig.3.13 Site occupancy at the “n”th nearest sites from YO particles.

3.2.5

結晶粒界近傍の分散形態

図 3.14 に傾角対称粒界に格子モンテカルロ法を施した結果を示す．図では Fe 原子
は表示していない．いずれの組成の場合においても結晶粒界面に沿って薄く均等に集積する青色の TiO2 粒子が目立つ．また，図 3.14(b) と (d) を比較すると，YO 粒子の増加とともに粒界から離れた位置にある TiO2 析出物は微細化する．これは，バルクでの検討で述べたように，比較的 YO 粒子が遠距離まで引力相互作用を及ぼし TiO2 粒子をトラップするためである．

0.2wt%YO-0.2wt%TiO2 系における結晶粒界近傍の様子を拡大して図 3.15 に示す．
図 3.15(a) より TiO2 粒子同士あるいは TiO2 粒子と原子空孔は y 方向（粒界に垂直方向）に 2 原子の組になっている．また，図 3.15(b) を見ると，TiO2 粒子と原子空孔は粒界面から第 1 近接距離程度に限定されていることがわかる．TiO2 粒子および原子空孔がどの程度粒界面に偏析しているかを定量的に示すため，シミュレーションセルを y
方向に 0.01[nm] 幅のブロックに分け，それぞれの種類の粒子の総数に対する各ブロッ
ク内の粒子の数の割合を図 3.16 にプロットした．y 方向中央の上下 2 つのピークおよび端面の 2 つのピークを合わせて簡便に偏析度合として比較すると，YO 粒子が約 8%，33

空孔約 14% の偏析に対して TiO2 粒子 は約 52% と半数に上る．
YO

TiO2

Atom Vacancy

y

y
x

x

z

z

(a) In Fe-0.2wt%YO-0.2wt%TiO2 alloy.

(b) In Fe-0.2wt%YO-0.4wt%TiO2 alloy.

y

y
x

x

z

z

(c) In Fe-0.6wt%YO-0.2wt%TiO2 alloy.

(d) In Fe-0.6wt%YO-0.4wt%TiO2 alloy.

Fig.3.14 Dispersion morphology of YO and TiO2 particles.

Grain Boundary

x
y
x

y

z

(a) Bird view.

(b) x-y plane.

Fig.3.15 Dispersion morphology of YO and TiO2 particles near the grainboundary in Fe-0.2wt%YO-0.2wt%TiO2 alloy.

34

G.B.

G.B.

G.B.

14

YO
TiO2
Vacancy

12

Ratio [%]

10
8
6
4
2
0
0

5

10

15

20

25

30

35

40

45

y coordinate [nm]Fig.3.16 The ratio of the distribution of YO, TiO2 and vacancy along they coordinate (normal to the grain boundaries) in Fe-0.2wt%YO0.2wt%TiO2 alloy.

3.3

結言

高温クリープ強度など ODS 鋼における優れた特性は，微細に分散した酸化物ナノクラスターによると考えられている．酸化物分散形態を原子レベルから検討するため，
本章ではモンテカルロシミュレーションにより ODS 鋼中の酸化物の分散形態を検討した．以下に得られた結果を示す．

(1) バルクの Fe-YO-TiO2 モンテカルロシミュレーションを様々な組成，温度条件で行った．YO 粒子は組成，温度に関わらず粒子数個程度のクラスターを形成し，シミュレーションセル内に均等に分布した．TiO2 粒子は比較的大きなクラスターを形成し，その寸法は YO 粒子の組成および温度条件によって変化した．

(2) 温度 1000[K] と 1400[K] での TiO2 粒子の析出物寸法を比較すると，やはり 1400[K]の方が粗大化する傾向が認められた．ただし，YO 粒子の組成が比較的高い Fe-

0.6wt%YO-0.2wt%TiO2 の組成では，1400[K] でも TiO2 粒子の粗大化が抑制さ35れていた．

(3) Fe-YO-TiO2 -Al 系のシミュレーションにより，Al 原子の添加はもともと小さく分散する YO 粒子をより凝集しにくくすることがわかった．また，YO 粒子は Al原子近傍に位置していた．

(4) 原子空孔を導入した Fe-YO-TiO2 系の LMC シミュレーションにより，原子空孔
は YO 粒子の周囲に集中することがわかった．その理由は，Fe の bcc 格子において，最近接距離での YO-YO 結合の強い反力を原子空孔が緩和するためである．

(5) Σ5 傾角対称粒界近傍についての検討では，結晶粒界から離れた領域ではバルクでの検討と同じ分散形態が得られる一方，結晶粒界面から第 1 近接面に結晶粒界を挟んで対を成すように TiO2 粒子と原子空孔の偏析がみられた．TiO2 粒子が特に顕著に粒界に偏析するが，YO 粒子にはそのような傾向は認められなかった．

36

第 4 章 単結晶のせん断変形シミュ
レーション
4.1

解析条件

図 4.1 にシミュレーション手順を示す．まず，bcc 構造を対象とした圧力 0[Pa] のNPT アンサンブル格子モンテカルロ法 (NPT-LMC 法) により，1500[K] から 600[K] まで 20[K] 刻みで温度を下げる焼鈍し法 (Simulated Annealing: SA) で初期構造を作成した．図 4.2 に得られた構造の例 (粒子数 241,920) を示す．(1̄1̄2) 面をすべり面とし，[111]
方向へのせん断を行うために，x,y,z 軸をそれぞれ bcc の [111]，[1̄10]，[1̄1̄2] としている．NPT-LMC 法で作成した分散形態は，前章の NVT アンサンブルで作成した物とほぼ同じである．その後，分子動力学法により 600[K] で 20[ps] の緩和計算を行った．
さらに，10[K] で 20[ps]，応力が 0 となるようにセル寸法を制御しながら緩和計算を行い，せん断シミュレーションの初期状態とした．ここで，10[K] の極低温としたのは転位に関する先行研究 (44)(45) と同じく温度の影響と外力のそれを切り分けるためである．

zx 成分のせん断は，毎ステップひずみ増分 ∆εzx = 5.0 × 10−7 を与えることで行った．温度は 10[K] とし，せん断シミュレーション中はせん断方向 (zx 方向) を除いた応力成分 (σxx , σyy , σzz ,σxy ,σyz ) を 0[GPa] に制御した．YO 粒子の組成は 0.2, 0.6wt%，TiO2粒子の組成を 0.2, 0.4wt% の 4 通りとし，酸化物を添加しない Fe 完全結晶の場合も解析した．いずれも，全方向周期境界条件でバルクとしている．

37

Structual Searchby Simulated Annealing
with NPT-LMC
from 1500[K] to 600[K]

13.0[nm]

Structual Relaxation 1
by Molecular Dynamics
at 600[K] for 20[ps]
19.7
[nm]

Structual Relaxation 2
by Molecular Dynamics
at 10[K] for 20[ps]-- z [112]y [110]x [111]

Shear Deformation
by Molecular Dynamics
at 10[K]strain rate ε zx 5.0x10 -7[1/fs]

11.2[nm]Fig.4.2 Example of simulation cell forshear deformation.

Fig.4.1 Simulation procedure for shear deformation.

4.2

解析結果

4.2.1

応力-ひずみ応答と欠陥発生の様相

図 4.3 に応力-ひずみ線図を示す．酸化物粒子の添加割合が非常に小さいため，応力ピークのひずみは異なるものの，いずれの組成においてもひずみ 0.03 まではほぼ同じ曲線となり，鉄のみの場合と弾性率は変わらない．一方，応力ピークひずみを比較す
ると，酸化物を添加した場合はすべての組成で Fe 単体より小さい．これは，Fe 単体が欠陥を一切含まない「完全結晶」であるのに対し，酸化物を含有する系は内部不均一性を有していることに起因する．また，Fe 単体では応力ピークを迎えてもすぐに応力が急減しないが，酸化物を添加した場合はピーク直後で急減する．酸化物を添加した組成同士を比較すると，Fe- 0.2wt%YO -0.2wt% TiO2 ，Fe- 0.2wt%YO -0.4wt%TiO2 ，

Fe- 0.6wt%YO -0.4wt%TiO2 の組成では応力ピークひずみに大きな違いは見られないのに対し，Fe-0.6wt%YO-0.2wt%TiO2 の組成では他に比べ高いピークひずみ・応力を示した．図 4.4 にせん断前の各組成での分散形態を示す．3.2.2 節で示したように，Fe-38

0.6wt%YO-0.2wt%TiO2 の組成が他の場合に比べて TiO2 酸化物が細かく分散している．応力ピークが最も低ひずみで生じた Fe-0.2wt%YO-0.4wt%TiO2 の内部欠陥・構造変化を図 4.5 に，最も高ひずみを示した Fe-0.6wt%YO-0.2wt%TiO2 のそれを図 4.6 に示す．図には可視化ソフト AtomEye(61) の Central Symmetry 値によって判別した，結晶構造からずれた「欠陥」を示している．粗大化した TiO2 の界面や YO 周りの乱れた部分が欠陥として表示されているが，図 4.5，図 4.6 中の矢印で示した析出物と母相の界面付近の乱れた部分は後述するように双晶であり，互いに合体し成長する．図 4.5 を見ると，大きい析出物の界面では比較的低ひずみで発生している．一方，TiO2 酸化物
クラスターが小さい 0.6wt%YO-0.2wt%TiO2 (図 4.6) は，界面からの双晶発生が図 4.5に比べて抑えられている．
Fe
0.2wt%YO-0.4wt%TiO2
0.6wt%YO-0.4wt%TiO2
6

0.2wt%YO-0.2wt%TiO2
0.6wt%YO-0.2wt%TiO2
Peak stress point

Shear stress σzx [GPa]5Peak stress5.56[GPa]Fe

4( εzx =0.0519 )4.20[GPa]

30.2wt%YO0.2wt%TiO2( εzx =0.0313 )20.2wt%YO0.4wt%TiO2

4.14[GPa]( εzx =0.0310 )1

0.6wt%YO0.2wt%TiO2

0

0.6wt%YO0.4wt%TiO2

4.91[GPa]( εzx =0.0392 )4.33[GPa]( εzx =0.0325 )

-1

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Shear strain εzx

Fig.4.3 Stress-Strain curves under shear deformation.

-- z [112]y [110]-- z [112]y [110]x [111]-- z [112]y [110]x [111](a) Fe-0.2wt%YO-0.2wt%TiO2.

-- z [112]y [110]x [111](b) Fe-0.2wt%YO-0.4wt%TiO2.

(c) Fe-0.6wt%YO-0.2wt%TiO2.

x [111](d) Fe-0.6wt%YO-0.4wt%TiO2.

Fig.4.4 Dispersion morphology before shear loading.
39

-z [112]
x [111]y [110](a) εzx = 0.030(b) εzx = 0.031(c) εzx = 0.032(d) εzx = 0.033

Fig.4.5 Central symmetry parameter visualization under shear deformation(Fe-0.2wt%YO-0.4wt%TiO2 ).

-z [112]x [111]y [110](a) εzx = 0.038(b) εzx = 0.039(c) εzx = 0.0395(d) εzx = 0.040

Fig.4.6 Central symmetry parameter visualization under shear deformation(Fe-0.6wt%YO-0.2wt%TiO2 ).

Fe-0.6wt%YO-0.2wt%TiO2 における欠陥発生後 (εzx = 0.05) のポテンシャルエネル
ギー分布を図 4.7 に示す．局所的なポテンシャルエネルギーが高くなり，(1̄1̄2) 面に欠陥が入っている様子が認められるが，欠陥周辺の原子配列は大きな乱れがなく，むしろ規則的に配列している．欠陥を境に結晶の向きが変化しているのみで，z 方向に関して同じ結晶方向を持つ領域が交互に現れることから双晶変形であることが確認でき
る．bcc 金属においては双晶面が {112} で双晶方向が ⟨111⟩ であり (62)(63) ，低温かつ高速に変形させる本シミュレーション条件を考慮すれば双晶が生じやすかったためと考
えられる．

40

[eV]
-1.43

-z [112]-1.57

x [111]y [110]

Fig.4.7 Potential energy distribution in Fe-0.6wt%YO-0.2wt%TiO2 alloy at
εzx = 0.05.

4.2.2

系の弾性剛性係数

原子弾性剛性係数 (AES) の全原子平均で求めた系の弾性剛性係数 BIJ の行列式変化(0)を図 4.8 に示す．図中 (b) には各組成の初期値 detBIJ で無次元化した変化も合わせて示している．図中 (a) の εzx = 0 において，酸化物を添加した組成では純鉄より detBIJが高い値を示している．また，TiO2 粒子の組成が高い程 detBIJ は高くなるが，YO 粒子の量を増やすとわずかに detBIJ が低下するものの TiO2 ほどの変化はない．図 (b)の無次元化した量で変形による変化を見ると，ひずみの増加とともにいずれの組成の場合も detBIJ が低下する．純鉄の場合は負の値，すなわち系の不安定 (57) に達した後に非弾性変形 = 内部でのすべり変形が発生し，それにより系の detBIJ も回復する．一方，初期不均一性を有する系では系全体が不安定となる前に非弾性変形が発生して系の detBIJ が回復する．その条件 (detBIJ がどの程度になれば非弾性変形が発生するのか) についてはマクロの降伏と同じく，未だ解明されていない問題である．

41

(0)

6

Normalized 6x6 determinant, det BIJ / det BIJ

6x6 determinant of elastic stiffness coefficient, det BIJ ,[ 1012 GPa6 ]0.2wt%YO-0.2wt%TiO2
0.6wt%YO-0.2wt%TiO2

Fe
0.2wt%YO-0.4wt%TiO2
0.6wt%YO-0.4wt%TiO2

5
4
3
2
1
Fe
0.2wt%YO-0.2wt%TiO2
0 0.2wt%YO-0.4wt%TiO2
0.6wt%YO-0.2wt%TiO2
0.6wt%YO-0.4wt%TiO2-1
00.01 0.02 0.03

0.04

0.05

0.06

0.07

Shear strain εzx

1.4
1.2
1.0
0.8
Peak stress point

0.6

1 εzx =0.0519
2 εzx =0.0392

0.4

3 εzx =0.0324

0.2
4 εzx =0.0315

0-0.2

5 εzx =0.0314

0

0.01

0.02

3
4
2

5
0.03

1
0.04

0.05

0.06

0.07

Shear strain εzx(a) Magnitude of detBIJ .

(b) Normalized value.

Fig.4.8 Change in the 6 × 6 determinant of system elastic stiﬀness coeﬃcient.

4.2.3

AES の評価

全原子の原子弾性剛性行列 (AES) を調べ，AES が負の原子の割合変化を図 4.9 に示した．完全結晶の Fe は前述のように系全体が不安定となり，ほとんど全ての原子のAES が負となるため，図 4.9(a) の縦軸のスケールが他とは異なることに注意されたい．
一方，酸化物を含有する不均一な系では無負荷平衡状態においても AES が負の原子が存在する．図 4.9(f) にせん断前の負の AES 原子の割合を示す．ほとんど YO 粒子であり，YO の添加割合に従って初期の負の AES 原子が増加することが分かる．ひずみ

0.01 を越えたあたりから，負の AES 原子が増加し始めるが，そのほとんどは Fe 原子
であり，YO 粒子はその割合をほとんど変えない．また，TiO2 粒子は AES が負となる原子は皆無であった．

42

Ratio of negative AES atoms, Nneg/Nall [%]Ratio of negative AES atoms, Nneg/Nall [%]

100

80

60

40
Peak stress point
εzx =0.0519
20

0

0.01

0.02

0.03

0.04

0.05

0.06

5

4

3

2

Peak stress point
εzx =0.0313

1

0

0.07

Fe
YO
TiO2

0.02

0.01

Ratio of negative AES atoms, Nneg/Nall [%]Ratio of negative AES atoms, Nneg/Nall [%]

5

Fe
YO
TiO2

4

3

Peak stress point
εzx =0.0310

1

0.01

0.02

0.03

0.04

0.05

0.06

0.07

5

3

2

Peak stress point
εzx =0.0390

1

0

0.01

0.02

Ratio of negative AES atoms, Nneg/Nall [%]Ratio of negative AES atoms, Nneg/Nall [%]3

Peak stress point
εzx =0.0325

1

0.04

0.05

0.06

0.07

0.30

4

0.03

0.04

(d) Fe-0.6wt%YO-0.2wt%TiO2

Fe
YO
TiO2

0.02

0.03

Shear strain εzx

5

0.01

0.07

4

(c) Fe-0.2wt%YO-0.4wt%TiO2

0

0.06

Fe
YO
TiO2

Shear strain εzx

2

0.05

(b) Fe-0.2wt%YO-0.2wt%TiO2(a) Fe00.04

Shear strain εzx

Shear strain εzx

2

0.03

0.05

0.06

0.25

Fe
YO
TiO2 : All 0%

0.20

0.15

Fe-0.2wt%YO-0.2wt%TiO2

0.10

0.05
0%
0

0.07

Fe only

Shear strain εzx(e) Fe-0.6wt%YO-0.4wt%TiO2

Fe-0.6wt%YO-0.2wt%TiO2

Fe-0.2wt%YO-0.4wt%TiO2

Fe-0.6wt%YO-0.4wt%TiO2(f) At initial.

Fig.4.9 Change in the ratio of negative AES atoms under shear deformation.

43

図 4.10，図 4.11 にそれぞれ Fe-0.2wt%YO-0.4wt%TiO2 ，Fe-0.6wt%YO-0.2wt%TiO2の負の AES 原子の分布を示す．AES が負となった Fe 原子を緑で，酸化物粒子は AESの正負で色分けして表示している．Fe-0.2wt%YO-0.4wt%TiO2 では図 4.10(c)〜(e) 中の矢印で示したように，ひずみの増加と共に大きな TiO2 析出物の界面に AES 行列式が負となる Fe 原子が集団的に発生，拡大している．双晶変形が発生するのはこれらの Fe の負の AES 原子が合体したときである．一方，大きな TiO2 析出物を持たない

Fe-0.6wt%YO-0.2wt%TiO2 の場合 (図 4.11)，TiO2 界面に負の Fe 原子が現れるが，図4.11(b)〜(f) 中の矢印で示したように，YO 粒子 (正負の区別なく) の周辺からも AESが負となる Fe 原子が発生している．この系では YO と TiO2 が複合した微細析出物を形成しているため，負の Fe 原子は図 4.10 のように決まった結晶方位を持つ集団で拡大するのではなく，様々な方向，場所で発生する．そのため，それぞれの不安定変形
モードが重なり，比較的高ひずみまで合体成長せずに増加していたと考えられる．

(a) εzx = 0(b) εzx = 0.010(c) εzx = 0.020(d) εzx = 0.025

det BIJ >0 det BIJ <0Fe

HiddenYO
TiO2-z [112]x [111]y [110](e) εzx = 0.030(f) εzx = 0.032(g) εzx = 0.0335

Fig.4.10 Distribution of negative AES atoms in Fe-0.2wt%YO-0.4wt%TiO2 alloy.

44

(a) εzx = 0(b) εzx = 0.020(c) εzx = 0.030(d) εzx = 0.035

det BIJ >0 det BIJ <0Fe

HiddenYO
TiO2-z [112]x [111]y [110](e) εzx = 0.0375(f) εzx = 0.0385(g) εzx = 0.0395

Fig.4.11 Distribution of negative AES atoms in Fe-0.6wt%YO-0.2wt%TiO2 alloy.

4.3

結言

酸化物分散形態による力学応答変化を調べるために，単結晶バルクを対象に (1̄1̄2)
面をすべり面とし，[111] 方向へ単純せん断させるシミュレーションを分子動力学法により行った．以下に得られた結果を示す．

(1) 酸化物を添加した場合は，完全結晶でないため，純鉄より応力ピークひずみが低下した．完全結晶の Fe のピークひずみは 0.05 であり，YO，TiO2 の組成を変えた 4 つの系のうち 3 つは応力ピークひずみ 0.03 程度を示したが，酸化物が細か
く均等に分散した構造を有する Fe-0.6wt%YO-0.2wt%TiO2 のみ 0.04 と他に比べて高い弾性限界を示した．

(2) 応力ピーク前後の原子構造を観察した結果，析出物界面近傍から欠陥が発生し，ひずみの増加と共に互いに合体し成長することで応力急減していた．

45

(3) 結晶方位と原子エネルギー分布から，(2) の欠陥発生→合体で生じる変形は双晶変形であることを示した．

(4) 高い弾性限界を示した Fe-0.6wt%YO-0.2wt%TiO2 の組成では，TiO2 クラスターが小さく界面からの欠陥生成が抑制されていた．

(5) 系全体の弾性剛性係数の行列式の変化 (系の安定性) を調べた結果，完全結晶 Fe
は応力ピークが系の不安定点 (行列式が負となる点) と一致した．酸化物を添加した系では内部不均一を有するため，系全体の不安定に達する前に応力急減していた．

(6) 行列式の値は TiO2 の組成によって大きく変わるが，YO による変化は小さかった．
(7) 各原子の原子弾性剛性係数の行列式の値 (AES) を調べた結果，無負荷平衡状態でも負の AES 原子が存在し，そのほとんどは YO 粒子であった．TiO2 粒子は負の AES を有するものは皆無であった．

(8) せん断変形による負の AES 原子の割合変化を調べた結果，ひずみ 0.01 近傍から
負の Fe 原子が急増した．負の AES を有する YO 粒子の割合は初期のそれからほとんど変わらなかった．TiO2 は応力ピーク後の変形下においてもほとんど負の原子は認められなかった．

(9) AES の正負で内部構造変化を観察した結果，比較的粗大な TiO2 析出物を形成し，弾性限界が 0.03 程度であった系では，TiO2 析出物-母相界面で AES が負となる

Fe 原子が集団的に発生し，そこが双晶変形の起点となっていた．一方，酸化物が細かく均等に分散し，高い弾性限界を示した Fe-0.6wt%YO-0.2wt%TiO2 の組成では，TiO2 の界面だけでなく YO 粒子周辺からも AES が負となる Fe 原子が様々な方向，場所で発生，成長していた．

(10) (9) より，Fe-0.6wt%YO-0.2wt%TiO2 では，負の AES となった Fe 原子の不安定
モード，方向が異なるため，界面からの双晶発生が抑制されたものと結論づけられる．

46

第 5 章 ねじれ粒界のへき開シミュ
レーション
5.1

解析条件

図 5.1 に模式的に示すように，bcc 格子の [001] 方向を z 軸として上半分，下半分をそれぞれ 18.45◦ ，−18.45◦ 回転させることで [001]Σ5 ねじれ粒界を導入する (粒子数240,000)．次に，前節と同様に，NPT-LMC 法により 1500[K] から 600[K] まで 20[K] 刻
みで温度を下げながら酸化物分散形態を作成する．図 5.2，図 5.3 に得られた分散形態を示す．3.2.5 項で示した [001]Σ5 傾角対称粒界の場合と同様に粒界面の TiO2 粒子の偏析が認められる．また，いずれの組成も粒界から離れたバルク部分での TiO2 粒子の粗大析出物は見られず，特に TiO2 が 0.2wt% の図 5.2(a)(c) では TiO2 粒子がすべて粒界近傍に位置している．すなわち，3 章の傾角粒界よりもねじれ粒界の方が TiO2 の偏析が大きい．
この構造を，分子動力学法により 600[K] で 20[ps]，さらに 10[K] で応力制御しながら 20[ps] の緩和計算を行うことで引張シミュレーションの初期状態とした．得られたシミュレーションセルに対して，z 方向に毎ステップひずみ増分 ∆εzz = 1.0 × 10−6 を与えることで粒界のへき開シミュレーションを行った．YO 粒子の組成を 0.2, 0.6wt%，

TiO2 粒子の組成を 0.2, 0.4wt% の 4 通りとし，純鉄の場合も解析した．温度は 10[K]，全方向周期境界条件とし，引張時の x, y 方向のセル辺長は不変 (ε′ = 0) とした．

47

[001]

[nm]
[010]

17.2

[100]

z[001]

[001]

17.2

[010]

x

[100]

y
9.03
9.03

Fig.5.1 Schematic of simulation cell for Σ5 twist grain boundary.

z [001]z [001]yyxx(a) Fe-0.2wt%YO-0.2wt%TiO2 .

(b) Fe-0.2wt%YO-0.4wt%TiO2 .

z [001]z [001]

y

y

x

x

(c) Fe-0.6wt%YO-0.2wt%TiO2 .

(d) Fe-0.6wt%YO-0.4wt%TiO2 .

Fig.5.2 Dispersion morphology of YO and TiO2 particles in bulk laminate
structure of Σ5 twist grain boundary.

48

y

y
x

x

z [001]z [001](a) Fe-0.2wt%YO-0.2wt%TiO2 .

(b) Fe-0.2wt%YO-0.4wt%TiO2 .

y

y
x

x

z [001]z [001](c) Fe-0.6wt%YO-0.2wt%TiO2 .

(d) Fe-0.6wt%YO-0.4wt%TiO2 .

Fig.5.3 Dispersion morphology of YO and TiO2 on Σ5 grain boundary.

5.2

解析結果

5.2.1

応力-ひずみ応答

図 5.4 に応力-ひずみ応答を示す．いずれの場合も εzz = 0.06〜0.07 近傍で，最初の応力ピークを示しており，酸化物を添加した場合は総じて純鉄よりわずかにピークひずみが小さくなるぐらいでほとんど同じ応答を示している．応力急減後，各組成とも

εzz = 0.08 付近で応力 σzz が下げ止まり，再び上昇に転じるのもほぼ同じである．一方，ひずみ εzz = 0.11 近傍からの応答は，後述するように内部の変形・破壊形態によって様々である．純鉄は εzz =0.125 近傍で応力がわずかに減少するが，εzz =0.15 を越えてから再び大きな応力ピークを示して粒界でのへき開発生により応力急減する．

Fe-0.6wt%YO-0.4wt%TiO2 は純鉄と似た応答を示し，εzz =0.193 近傍で第 1 ピーク応力より高い応力 (9.86[GPa]) を示して粒界でへき開した．他の系は第 2 ピークが不明確であり，特に，Fe-0.2wt%YO-0.4wt%TiO2 の組成では，第 1 ピーク後 εzz = 0.139 で最49

大 8.10[GPa] の応力を示した後，粒界でへき開した系に比べると段階的な応力減少を示している．
Pure Iron(a) Fe-0.2wt%YO-0.2wt%TiO2(c) Fe-0.6wt%YO-0.2wt%TiO2(b) Fe-0.2wt%YO-0.4wt%TiO2(d) Fe-0.6wt%YO-0.4wt%TiO2

10
9

Tensile stress σ zz [GPa]8

1st peak stressfracture point

7Fe

9.38[GPa]( εzz =0.0639 )

8.33[GPa]( εzz =0.177 )60.2wt%YO0.2wt%TiO2

9.27[GPa]( εzz =0.0609 )

8.49[GPa]( εzz =0.179 )40.2wt%YO0.4wt%TiO2

9.24[GPa]( εzz =0.0602 )8.10[GPa]( εzz =0.139 )30.6wt%YO0.2wt%TiO2

7.98[GPa]( εzz =0.168 )2

0.6wt%YO0.4wt%TiO2

9.28[GPa]( εzz =0.0620 )9.07[GPa]( εzz =0.0595 )

5

1
0
0

0.05

0.10

0.15

0.20

9.86[GPa]( εzz =0.193 )0.22

Tensile strain εzz
Fig.5.4 Stress-strain curves under tension (bulk laminate structure of Σ5
twist G.B.).

5.2.2

変形破壊の様相図 5.5〜図 5.9 に各組成の原子配置変化を，各原子のポテンシャルエネルギーで色分けして示す．純鉄の場合を見ると，図 5.5(c) と (d) の間 (第 1 応力ピーク前後) に発生した帯状の分布が図 5.5(d) と (e) の間 (応力の下げ止まる点) にシミュレーションセル
全体に伝播している．帯状の分布の境目で相対的なすべり変形を生じており，特に図5.5(e) では比較的高エネルギーを示す領域に原子配置のずれが見られる．酸化物を添加した他の系の図 5.5〜図 5.9 の (a)〜(g) を見ても，せん断帯の方向に違いがあるが，εzz =0.01 までの変化は大きく変わらない．しかし，その後の変形の進行と破壊の形態
は各組成で大きく異なるため，以下それぞれの系について記述する．

Pure-Fe
図 5.5(f) を見ると帯状の分布の一部境界領域に高エネルギー箇所が生じており，その後図 5.5(g)，さらに (h) にかけて系全体に高エネルギー箇所が現れ，帯状分布は不明

50瞭となる．図 5.5(f) と (h) を参照すると，高エネルギーとなっている部分は (f) の帯状
領域の境目ではなく，新たなすべりが生じていることが確認できる．新たなすべりが系全体に伝播した図 5.5(g) と (h) の間から応力が急増し，その変形が飽和した後，最終的に図 5.5(i) に示すように粒界でのへき開を生じた．

Fe-0.2wt%YO-0.2wt%TiO2
図 5.6 で同じ位置に点在する濃い青の原子はエネルギーが低い酸化物粒子である．図5.2(a) と比べると YO 粒子であることが分かる．純鉄の場合と同様に，結晶内部での変形が図 5.6(h) で飽和した後，すべりと粒界が交差する高いエネルギー部分からへき開を生じて破断した．

Fe-0.2wt%YO-0.4wt%TiO2
応力の第 2 ピークが最も低ひずみで生じ，その後の応力減少も段階的であった系で，結晶での変形飽和後の図 5.7(h) を見ると，粒内き裂による高エネルギー領域 (赤い帯状の領域) が下の結晶をほぼ縦断している．その形態は，粒界と結晶内部のすべり交差部を接続したようなものとなっている．その後も ε′ = 0 の条件のため，このき裂が開口するように変形が生じている (図 5.7(i))．

Fe-0.6wt%YO-0.2wt%TiO2
この系も応力第 2 ピーク後，段階的に応力が減少していたが，やはり粒内き裂を生じている (図 5.8(h))．ただし，先述の Fe-0.2wt%YO-0.4wt%TiO2 とは異なり，粒界を起点とせず，完全に粒内から発生している．図 5.2(c) を見ると多数の YO クラスター
が結晶粒内に存在しており，図 5.8 でも青い粒子として確認できるが，図 5.8(g) の上部で斜めの高エネルギー領域に複数の YO 粒子が認められる．図 5.8(h) のき裂を見るとその YO を起点としていることが確認できる．

Fe-0.6wt%YO-0.4wt%TiO2
最も高い応力を示した系である．εzz =0.2 の図 5.9(i) で粒界近傍でのへき開が発生し始めているが，粒内き裂も認められる．また，へき開は TiO2 が偏析している粒界

51

そのものではなく，その周囲の Fe 原子部分で生じており，この系の界面の結合が強いことが分かる．図 5.3 を再度確認すると，この系では TiO2 が粒界面に均等に分布していることが分かる．一方，同様に TiO2 が粒界に多い Fe-0.2wt%YO-0.4wt%TiO2
は第 2 ピークが最も低ひずみで粒内き裂と生じているが，図 5.3(b) を見ると TiO2 の分布にむらがあり，図 5.3(d) よりも TiO2 同士が密接しているように見られる．これ
が粒界-粒内交差すべり部へのき裂発生を促した可能性がある．粒界面上部の高エネル
ギー領域が起点となってへき開が発生しているが，同時に Fe-0.2wt%YO-0.4wt%TiO2 ，Fe-0.6wt%YO-0.4wt%TiO2 と同様なき裂と思われる高エネルギー領域も粒界の近くで発生している．しかし，前項で示したように応力-ひずみ応答は粒界でへき開を起こした純鉄と同様の挙動となっていることから粒界でのへき開が支配的な破壊モードと考えられる．

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

[eV]
-1.37

-1.56
z [001]x(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.5 Snapshots colored by potential energy (pure iron).
52y(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

[eV]
-1.37

-1.56
z [001]yx(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.6 Snapshots colored by potential energy (Fe-0.2wt%YO-0.2wt%TiO2 ).

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

[eV]
-1.37

-1.56
z [001]x(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175y(i) εzz = 0.200

Fig.5.7 Snapshots colored by potential energy (Fe-0.2wt%YO-0.4wt%TiO2 ).
53

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

[eV]
-1.37

-1.56
z [001]yx(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.8 Snapshots colored by potential energy (Fe-0.6wt%YO-0.2wt%TiO2 ).

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

[eV]
-1.37

-1.56
z [001]x(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175y(i) εzz = 0.200

Fig.5.9 Snapshots colored by potential energy (Fe-0.6wt%YO-0.4wt%TiO2 ).
54

系の弾性剛性係数と AES の評価

5.2.3

系の弾性剛性係数
原子弾性剛性係数 (AES) の全原子平均で求めた系の弾性剛性係数 BIJ の行列式変化(0)を初期値 detBIJ で無次元化して図 5.10 に示す．前章とは異なり，純鉄でも粒界による内部不均一性を有するため，系全体の行列式は負に達することなく最初の応力ピー
クを迎える．また，図 5.4 の応力変化と同様に，組成毎の違いが少ない．

Normalized 6x6 determinant, det BIJ / det BIJ(0)Pure Iron(a) Fe-0.2wt%YO-0.2wt%TiO2(c) Fe-0.6wt%YO-0.2wt%TiO2(b) Fe-0.2wt%YO-0.4wt%TiO2(d) Fe-0.6wt%YO-0.4wt%TiO2

1.0

0.50

0.8
0.45

First peak point

3 εzz =0.0609

1 εzz =0.0639

4 εzz =0.0602

2 εzz =0.0620

5 εzz =0.0595

0.6
0.40
0.4
0.35

2

0.2

3
4

1

5

0.30
0

0.05

0.10

0.15

0.20

0.22

0.05

0.055

0.06

0.065

0.07

0.075

Tensile strain εzz

Fig.5.10 Change in the 6 × 6 determinant of system elastic stiﬀness coeﬃcient.

AES の評価
図 5.11〜図 5.15 に AES が負となった原子の分布を示す．純鉄を見ると，図 5.11(a)の引張前に負の AES 原子が粒界に認められる一方，図 5.11(b) では粒界面の負の AESが減少し，応力ピークに近い εzz =0.05 の図 5.11(c) で再び増加している．これは引張による力学状態の均一化⇒界面原子の安定化と，その後の均一限界後の不安定原子の増加 (64) によるものと考えられる．その後，図 5.11(d)〜(f) で粒界から結晶粒内にすべりが発生するため，粒内部分に多数の負の AES 原子が認められる．一方，図 5.11(h)は応力-ひずみの第 2 ピークにほぼ一致するが，(g) → (h) で応力上昇する際に再び内部

55

では「均一化」を生じるため，AES が負となる原子が大幅に減少している．また，へき開直前の図 5.11(h) からへき開を生じる粒界に負の AES 原子が多く存在することが分かる．粒界でのへき開を起こした Fe-0.2wt%YO-0.2wt%TiO2 の場合 (図 5.12) もほぼ同様のメカニズムで，図 5.12(h) が応力の第 2 ピークに相当するが，やはりへき開す
る粒界に多数の負の Fe 原子が認められる．その後の図 5.12(i) から分かるように，粒界面に沿って偏析した TiO2 の間でへき開しているが，TiO2 そのものは負となっていない．
粒内き裂を生じた Fe-0.2wt%YO-0.4wt%TiO2 の場合，表面原子は必ずしも AES が負ではないためそれほど明確ではないが，図 5.7(h) と図 5.13(h) を比較すると，き裂
位置に AES が負となる Fe 原子が帯状に連なって見える．その帯には YO 粒子およびTiO2 粒子が含まれており，これらを起点としてき裂が発生したことが分かる．粒界を起点とせず，完全に粒内からき裂が発生した Fe-0.6wt%YO-0.2wt%TiO2 では，引張後期 (εzz = 0.15) でも下部の結晶に比べ上部の結晶部分に負の AES 原子が多く残っていたことが分かる．図 5.8(h) と図 5.14(h) を比較するとやはりき裂近傍に負の AES 原子が存在し，多くの YO 粒子を含んでいるがき裂面や方向などはこの図では議論できな
い．第 2 ピークが第 1 ピークより大きく，へき開ひずみが最大であった Fe-0.6wt%YO-

0.4wt%TiO2 の場合，最初の応力ピーク後に発生したすべり面に沿った負の AES 原子が，変形後期 (図 5.15(g)) でも解消されず，図 5.15(g) → (h) では 2 つの負の AES 原子
面が 1 つになるような変化をしている．これは最初のピーク後に導入された双晶欠陥間で結晶の回転または相変態を生じたものと推測される．なぜこの系の下の結晶だけにこれが生じたのか，また YO や TiO2 の分布形態との関係などは現時点では不明である．

56

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

det BIJ >0 det BIJ <0Fe

HiddenYOTiO2z [001]yx(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.11 Distribution of negative AES atoms (pure iron).

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

det BIJ >0 det BIJ <0Fe

HiddenYOTiO2z [001]x(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.12 Distribution of negative AES atoms (Fe-0.2wt%YO-0.2wt%TiO2 ).
57y(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

det BIJ >0 det BIJ <0Fe

HiddenYOTiO2z [001]yx(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.13 Distribution of negative AES atoms (Fe-0.2wt%YO-0.4wt%TiO2 ).

(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

det BIJ >0 det BIJ <0Fe

HiddenYOTiO2z [001]x(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175(i) εzz = 0.200

Fig.5.14 Distribution of negative AES atoms (Fe-0.6wt%YO-0.2wt%TiO2 ).
58y(a) εzz = 0(b) εzz = 0.025(c) εzz = 0.050(d) εzz = 0.075(e) εzz = 0.100

det BIJ >0 det BIJ <0Fe

HiddenYOTiO2z [001]x(f) εzz = 0.125(g) εzz = 0.150(h) εzz = 0.175y(i) εzz = 0.200

Fig.5.15 Distribution of negative AES atoms (Fe-0.6wt%YO-0.4wt%TiO2 ).

5.3

結言

酸化物の組成の違いによる粒界の力学特性変化を調べるため，分子動力学法により

Σ5 ねじれ粒界を導入した系の [001] 方向への引張シミュレーションを行った．以下に得られた結果を示す．

(1) 焼鈍し法と分子動力学法により得られた [001]Σ5 ねじれ粒界は，3.2.5 項で示した [001]Σ5 傾角対称粒界の場合と同様に粒界面の TiO2 粒子の偏析が認められた．
しかし，いずれの組成もバルク部分での TiO2 粒子の粗大析出物は見られず，傾角粒界よりもねじれ粒界の方が TiO2 の偏析が大きい．

(2) 純鉄および酸化物を添加した 4 つの組成いずれも εzz = 0.06〜0.07 近傍で第 1 応力
59

ピーク (9.07〜9.38[GPa]) に達した後，εzz = 0.08 近傍で下げ止まり，εzz = 0.011
近傍まで同様の応答を示したが，以降は組成により応力-ひずみ曲線が大きく異なった．

(3) 純鉄の場合，εzz = 0.177 で第 2 応力ピーク (8.33G[Pa]) を示し粒界でへき開した．Fe-0.6wt%YO-0.4wt%TiO2 は純鉄と同様の挙動を示し，εzz = 0.193 で第 1
応力ピークより高い第 2 応力ピーク (9.86[GPa]) を示し粒界でへき開した．その他の組成では第 2 ピークが不明確であり，特に Fe-0.2wt%YO-0.4wt%TiO2 では

εzz = 0.139 で第 1 ピーク後の最大応力を示したあと，へき開した系に比較して段階的に応力が減少した．

(4) 原子毎のポテンシャルエネルギー分布により内部構造の変化を調べた．いずれの組成も，εzz = 0.1 近傍までで第 1 応力ピーク前後に発生した帯状の分布が各結晶領域全体に伝播した．その後，組成で応力-ひずみ応答変化に違いが生じ始め
る εzz = 0.1〜0.15 の間で, いずれの組成でも追加発生したすべりによる高エネル
ギー領域が発生した．

(5) (4) での新たなすべり面による変形が飽和したとき，純鉄では粒界のへき開が起こった．Fe-0.2wt%YO-0.2wt%TiO2 も同様ですべりと粒界が交差する高エネル
ギー領域からへき開が生じた．Fe-0.2wt%YO-0.4wt%TiO2 では粒界と結晶内部のすべり交差部を接続するような粒内き裂が発生した．Fe-0.6wt%YO-0.2wt%TiO2
では完全に粒内からき裂が発生していた．Fe-0.6wt%YO-0.4wt%TiO2 では粒界に偏析した TiO2 とその周囲の Fe 原子の間でへき開を起こすと同時に粒内き裂も認められた．

(6) Fe-0.2wt%YO-0.4wt%TiO2 と Fe-0.6wt%YO-0.4wt%TiO2 の粒界での TiO2 粒子の分布の違いから，むらのある，一部密集した構造が粒界-粒内交差すべりへのき裂発生を促した可能性が示唆された．

(7) 前章と異なり，粒界による内部不均一性により純鉄でも系の弾性剛性係数行列式が負となることはなかった．また，組成による変化の違いも少ない．

60

(8) 負の AES 原子分布により内部構造変化を調べた．純鉄において引張初期で粒界での負の AES 原子が，減少した後，増加に転じた．これは引張による力学状態の均一化⇒界面原子の安定化と，その後の均一限界によるものと考えられる．

(9) へき開を生じる粒界には AES が負となる Fe 原子が多く存在した．Fe-0.2wt%YO0.2wt%TiO2 の場合，偏析した TiO2 間でへき開しているが TiO2 そのものは負となっていない．

(10) Fe-0.6wt%YO-0.2wt%TiO2 の場合，引張後期 (εzz = 0.15) でも下部の結晶に比べ上部の結晶に負の AES 原子が多く残っていた．き裂を生じた後，き裂近傍に負の AES 原子が存在していた．

(11) Fe-0.6wt%YO-0.4wt%TiO2 の場合，引張後期 (εzz = 0.15) でも (10) と同じく最初の応力ピークで生じた AES 負原子が解消されず，εzz = 0.15 → 0.175 で 2 つの負の AES 原子面が 1 つになるように変化した．これは最初の応力ピーク後に導入された双晶欠陥間で結晶の回転または相変態を生じたものと推測される．

61

第 6 章 らせん転位と析出物の相互
作用

6.1

解析条件

6.1.1

分子動力学シミュレーション

NPT-LMC 法による焼鈍と分子動力学法による構造緩和によって作成した析出構造(図 4.4，粒子数 241,920) を基にシミュレーションセルを作成した．まず，周期境界を考慮してセルを x[111]-y[1̄10] 平面に 3×3 並べ拡大する (図 6.1，粒子数 2,177,280)．次に，
すべり面法線方向の z[1̄1̄2] 方向にセルを 100 分割し，TiO2 がもっとも多く存在する領域が z 方向の中央すべり面となるように周期境界下でシフトさせる．y[1̄10] 方向についても同様の操作を行い，y 方向セル辺長の 1/3 および 2/3(周期を 3 倍にしたため) の距離で TiO2 の濃度が最大となるようにした．その後，y,z 方向の周期境界を外し，図6.2 に示すように y 方向端面から 1/5 までの領域を上下でそれぞれ ⃗b/2[111]，−⃗b/2[111]ずつ x 方向にずらすことで，セル辺長の 1/5 の位置にらせん転位を導入した．

y,z 方向は端面 3 原子層を固定，x 方向は周期境界条件下で 10[K]，20[ps] 間，セルをスケーリングしながらの分子動力学計算を行い，垂直応力 (σxx , σyy , σzz ) を 0[GPa] にしながら緩和させた．組成は，TiO2 が粗大な Fe-0.2wt%YO-0.4wt%TiO2 と比較的微細な Fe-0.6wt%YO-0.2wt%TiO2 の 2 通りとした．転位を導入した中央すべり面の構造をポテンシャルエネルギーで着色して図 6.3 に示す．エネルギーの低い青色の粒子は酸化物粒子，赤色の粒子が転位芯付近のエネルギーの高い原子である．ただし，端面の表面原子も赤く着色している．このセルに，毎ステップひずみ増分 ∆εzx = 2 × 10−7 与えることで zx 方向にせん断ひずみを与え，転位を [1̄10] 方向に進行させるシミュレー

62

ションを行った．温度は 10[K] とし，垂直応力 (σxx , σyy , σzz ) は 0[GPa] に制御した．

-z [112]y [110]19.7[nm]x [111]39.0[nm]33.6[nm]Fig.6.1 Dimensions of simulation cell for interaction between screw dislocations and oxides.
dislocation line-z [112]y [110]x [111]

b/2
b/2

Fig.6.2 Setting of a screw dislocation.

[eV]

-1.44

-z [112]y [110]x [111]-1.59(a) Fe-0.2wt%YO-0.4wt%TiO2 .

(b) Fe-0.6wt%YO-0.2wt%TiO2 .

Fig.6.3 Snapshots on center slip plane colored by potential energy after initial
relaxation.

63

6.1.2

NEB 法によるエネルギーバリア評価

NEB 法によるエネルギー経路探索は計算量が増えるため，拡大しない元のセル (図4.4，粒子数 241,920) を用いた．分子動力学計算と同じ様に y 方向，z 方向それぞれの中央面で TiO2 粒子の濃度が最大となるようにシフトさせた後，y 方向端面からセル
辺長の 1/4 にらせん転位を導入した (図 6.4)．図では z 方向中央の転位芯だけでなく端面においても転位を導入するためにずらした部分の「乱れ」が表示されている．析
出物の x 方向直径は 2.89[nm]，析出物間距離は 11.1[nm] である．その後，NEB 解析を行う際に必要な各位相点 (原子配置) を作成するため，分子動力学法によるせん断シ
ミュレーションを行った．具体的には，0.1[K] で 20[ps] 初期緩和させた後，毎ステップ∆εzx = 2.5 × 10−7 のせん断を与えて転位を動かした．境界条件は前節の分子動力学解析と同じである．
可視化ソフト AtomEye で転位と析出物が相互作用する過程を観察し，時系列順に

21 点の原子配置を抽出し，与えたせん断ひずみ分戻して (アフィン変換)NEB 解析を開始する際の初期位相点とした．そして，始めと終わりの位相点を温度 10−4 [K] の分子動力学計算によりポテンシャルエネルギーが Local Minimum となるようにし，NEB
法における始めと終わりの状態とした．以上の手順で作成した位相点を用いて，最急降下法による収束計算により NEB 計算を行った．始状態の位相点 (3N 個の成分

x1 , y1 , z1 , x2 , y2 , z2 , · · · , xN , yN , zN ) から終状態への経路を直線 (終点から始点を引いたesxe1 − xs1 , y1e − y1s , · · · , zN
− zN
) としてそれを分割したものを NEB 開始時の中間点とせずに，分子動力学計算の結果を初期点と用いることで発散することを防いでいる．
2.89[nm](Diameter on x direction )-z [112]y [110]19.6[nm]

x [111]

11.1[nm](Distance between precipitates )12.9[nm]Fig.6.4 Setting of a screw dislocation for NEB calculation.

64

6.2

解析結果

6.2.1

応力-ひずみ応答および弾性剛性係数

図 6.5 に応力-ひずみ線図を示す．Fe-0.2wt%YO-0.4wt%TiO2 はひずみ εzx = 0.010
近傍，Fe-0.6wt%YO-0.2wt%TiO2 はひずみ εzx = 0.012 近傍から転位が動き出しているが，εzx = 0.017 付近まで応力-ひずみ曲線に現れる変化は極めて微小であり，差がないように見える．その後，応力上昇が鈍化するが，Fe-0.2wt%YO-0.4wt%TiO2 では

εzx = 0.0230 近傍，Fe-0.2wt%YO-0.4wt%TiO2 では εzx = 0.0216 近傍で再び線形的に上昇する．内部を観察すると転位がセル壁に達していた．
図 6.6 に系の弾性剛性係数の行列式変化を示す．εzx = 0.016 近傍まで組成で大きな違いは見られない．応力-ひずみに差を生じた点に対応して，εzx = 0.017 近傍から，

Fe-0.6wt%YO-0.2wt%TiO2 の減少速度がまず鈍化し，応力-ひずみ曲線を上下反転させたような挙動を示している．

Normalized 6x6 determinant, det BIJ / det BIJ

(0)

2.5

Shear stress σzx [GPa]

2.0

1.5

1.0

0.5
Fe-0.2wt%YO-0.4wt%TiO2
Fe-0.6wt%YO-0.2wt%TiO2
0

0.004

0.008

0.012

0.016

0.020

0.024

Shear strain εzx

1.005
1.000
0.995
0.990
0.985
0.980
0.975
0.970
0.965
Fe-0.2wt%YO-0.4wt%TiO2
Fe-0.6wt%YO-0.2wt%TiO2

0.960
0.955

0

0.004

0.008

0.012

0.016

0.020

0.024

Shear strain εzx

Fig.6.5 Stress-strain diagram under shear Fig.6.6 Change in the 6×6 determinant of
deformation (system with a screw
system elastic stiﬀness coeﬃcient.
dislocation).

65

6.2.2

転位の挙動 (粗大析出物)

図 6.7 に Fe-0.2wt%YO-0.4wt%TiO2 の系でらせん転位がセル壁に抜けるまでの挙動を示す．すべり面付近のみを抽出して，central symmetry で判別した画像である．図6.7(c)〜(e) にかけて，転位は図 6.7(c) で緑色の矢印で示した析出物は大きな抵抗なく通過するが，赤色の矢印で示した中央の 2 つの析出物の位置ではピンニングされている．図 6.7(f) では析出物の前縁 (転位が最初に析出物にぶつかる箇所) を転位線とつな
ぐような欠陥が生成されており，転位通過後の図 6.7(g) でも赤色の矢印で示した部分に欠陥を残している．図 6.8 に AtomEye の coordination number を用いてバルク以外の原子を抜き出し，上と正面から見たものを示す．転位と析出物がぶつかる部分で生じた欠陥がその後，析出物-母相界面を沿うように移動する様子が見られた (6.8 矢印)．
ただし，析出物を挟んで転位がすべり面を変えた様子はない．したがって，図 6.9 のようにバルクと析出物での misfit を埋め合わせるように欠陥が生じたと考えるがその詳細は不明である．
図 6.7 に戻ると，(g)〜(l) では析出物との切り合いによって転位ループが形成されている．このループ形成時には，図 6.10 の front view から分かるように明確に交差すべりを起こしている．したがって，図 6.11 に示すように交差部を起点として格子不整合を埋め合わせる欠陥が生成し，すべり面を変えて運動していた転位が析出物界面で再び同じすべり面に戻ることで転位ループを形成したと考えられる．
このような欠陥生成が析出物と転位の切り合いの過程で多数認められ，様々な欠陥を析出物間に残しながら転位は通過した．
図 6.12 に転位が動き出す際の負の AES 原子の分布を示す．せん断シミュレーションの比較的初期において，図 6.12(a)〜(c) を見ると，転位芯付近の原子の AES が正から負となり，わずかに転位が動いた後，再び正となり転位が動かない状態への推移が見られた．その後，断続的に転位芯の一部で AES が負となる原子が生じ，図 6.12(d) →(e) に見られるように転位芯全体に広がり，その数が急増した後，転位は動き続けるよ
うになった．

66

(a) εzx = 0.011(b) εzx = 0.012(c) εzx = 0.014(d) εzx = 0.0144(e) εzx = 0.0152(f) εzx = 0.0160(g) εzx = 0.0168(h) εzx = 0.0176(i) εzx = 0.0184(j) εzx = 0.0192(k) εzx = 0.0200(l) εzx = 0.0208-z [112]x [111](m) εzx = 0.0216(n) εzx = 0.0230

Fig.6.7 Dislocation behavior. (Fe-0.2wt%YO-0.4wt%TiO2 )67

y [110]

Top
view

Front
view

εzx = 0.0142

εzx = 0.0146

εzx = 0.0150

Top
view

-z [112]

Front
view

y [110]εzx = 0.0154

εzx = 0.0158

x [111]

Fig.6.8 Top and front view of defect formation in Figs. 6.7(c)〜(f) by coordination number. (Fe-0.2wt%YO-0.4wt%TiO2 )Fig.6.9 Schematic of defect formation in Figs. 6.7(c)〜(f). (Fe-0.2wt%YO0.4wt%TiO2 )

68

Top
view

-z [112]

Front
view

y [110]x [111]εzx = 0.0166

εzx = 0.0184

Cross slip and change slip plane.

Fig.6.10 Top and front view of defect formation during Figs. 6.7(g)〜(l) bycoordination number. (Fe-0.2wt%YO-0.4wt%TiO2 )Fig.6.11 Schematic of formation mechanism of dislocation loop during Figs.
6.7(g)〜(l).

69

(a) εzx = 0.0016(b) εzx = 0.0018(c) εzx = 0.0020-z [112]det BIJ >0 det BIJ <0y [110]Fe

HiddenYOx [111]TiO2(d) εzx = 0.0108(e) εzx = 0.0110

Fig.6.12 Negative AES atoms in Fe-0.2wt%YO-0.4wt%TiO2 alloy.

6.2.3

転位の挙動 (微細析出物)

図 6.13(e) に Fe-0.2wt%YO-0.4wt%TiO2 の系に残された欠陥の様子を centarl sym-

metry で判別して示す．Fe-0.6wt%YO-0.2wt%TiO2 で見られたのと同様の「しっぽ状」の欠陥が見られた．転位がセル壁に抜けた後，さらにひずみを与えてこの欠陥がループ状であることを確認している．図 6.13(a)〜(c) に見られるように転位が複合酸化物とぶつかってできた欠陥が x 方向 ([111] 方向) に動いて，別の析出物にぶつかり欠陥がy 方向 ([1̄10] 方向) に成長している．
図 6.14 に欠陥が生じるときに，AES が負となった原子の分布を示す．転位とぶつかったナノクラスター周辺に負の Fe 原子が生じており，その集団が x 方向に移動した後，y 方向に成長している．その際，YO 粒子は負の AES となるものが生じるが，TiO2 には生じていない．

70

(a) εzx = 0.0136(b) εzx = 0.0140(c) εzx = 0.0144-z [112]y [110]x [111](d) εzx = 0.0148(e) εzx = 0.0220

Fig.6.13 Defect nucleation by cutting of screw dislocation in Fe-0.6wt%YO0.2wt%TiO2 alloy.

(a) εzx = 0.0136(b) εzx = 0.0140-z [112](c) εzx = 0.0144

det BIJ >0 det BIJ <0y [110]Fe

HiddenYOx [111]TiO2(d) εzx = 0.0148

Fig.6.14 Negative AES atoms in Fe-0.6wt%YO-0.2wt%TiO2 alloy.

71

6.2.4

NEB 法によるカッティングバリア評価

図 6.15 に NEB 法により求めたエネルギーバリアを，図 6.16 に各位相点に対応する
原子配置のスナップショットを示す．横軸の Reaction Coordinate(RC) は始状態と終状態を結ぶ位相空間上の点の数であって実際の空間の位置などを表すものでないことに注意されたい．図 6.15 では，RC0 での緩和が不十分で，RC1(図 6.16(a)，転位が中央の TiO2 に接触している状態) でエネルギー極小値をとったため，その点を基準に表示している．RC9 で 1.189[eV/nm] の極大値を示すが，図 6.16 を見ると転位が中央のTiO2 をカッティングし中心近傍に存在する状態である．RC11 でわずかにエネルギー
バリアは低下するが，1 つ目のピークと同程度のピーク RC13(1.222[eV/nm]) をとっている．図 6.16(c)〜(d) を見ると，この 2 つのピークはらせん転位が TiO2 析出物の中心付近を通過するときであり，RC11 の極小点はまさに析出物中心に転位が存在するときに対応する．この谷が深ければ転位は析出物にトラップされるものと考えるが，本結果からはほとんどトラップしないものと考える．

1.4
1.2219[eV/nm]

1.1893[eV/nm]

Energy barrier [eV/nm]

1.2
1.0

1.1444[eV/nm]

0.8
0.6
0.4
0.2
0

0

5

10

15

20

Reaction coordinate

Fig.6.15 Energy barrier for screw dislocation cutting a TiO2 precipitate.

72

(a) RC1(b) RC5(c) RC9(d) RC11(e) RC13(f) RC19

Fig.6.16 Snapshots of atoms at each reaction coordinate (central symmetry
parameter).

6.3

結言

本章では，第 4 章のせん断シミュレーションで作成した析出構造を拡大してらせん転位を導入し，分子動力学法によりせん断ひずみを与え [1̄10] 方向に転位を運動させる
シミュレーションを行った．対象は，比較的大きな TiO2 析出物を持つ Fe-0.2wt%YO-

0.4wt%TiO2 および酸化物が微細に分散した Fe-0.6wt%YO-0.2wt%TiO2 とした．また，
Fe-0.2wt%YO-0.2wt%TiO2 については，第 4 章で扱ったセルにらせん転位を導入し，NEB 法によるカッティングバリア評価も行った．以下に得られた結果を示す．
(1) Fe-0.2wt%YO-0.4wt%TiO2 はせん断ひずみ約 0.010，Fe-0.6wt%YO-0.2wt%TiO2
は約 0.012 で転位が動き出したが，応力-ひずみ曲線は εzx = 0.017 までほとんど同じであった．εzx = 0.017 からの差は，微細析出した Fe-0.6wt%YO-0.2wt%TiO2
が粗大析出の Fe-0.2wt%YO-0.4wt%TiO2 より応力上昇の鈍化が大きくなったた73

めである．弾性剛性係数の行列式の変化も同様で，εzx = 0.016 までほとんど同じで，そこから Fe-0.6wt%YO-0.2wt%TiO2 の行列式低下が停滞したために，先の応力-ひずみ線図を上下反転させたような挙動を示した．

(2) 比較的大きな TiO2 析出物を持つ Fe-0.2wt%YO-0.4wt%TiO2 の系では，大きなTiO2 析出物と転位の切り合いにより多数の欠陥を残した．その際，[111] 方向において析出物を挟んですべり面が明確に変わる場合と変わらない場合があり，変わる場合は交差部を起点とした転位ループの形成が確認された．すべり面を変えないで析出物界面に欠陥を残したのはミスフィットによるものと推測されるが，その詳細は現時点では不明である．

(3) 負の AES 原子の分布を調べた結果，転位芯付近の原子はせん断シミュレーショ
ン初期には AES が負とならなかったが，転位が動き続けるときに転位芯の Fe 原子 AES がになっていることが分かった．

(4) 酸化物が細かく分散する Fe-0.6wt%YO-0.2wt%TiO2 の組成でもまた，(2) と同様の転位ループ形成が見られた．その際，転位と析出物の切り合いから欠陥が生じ，[111] 方向へ動いた後，[1̄10] 方向へ欠陥が成長する過程が観察された．AES
負原子の分布を調べると，上記欠陥生成時に, 析出物-母相界面から負の Fe 原子が生じており，負の YO 粒子も観察された．TiO2 には AES が負の原子は生じな
かった．

(5) NEB 法によりらせん転位が TiO2 析出物 (直径:2.89[nm]，析出物間距離:11.1[nm])をカッティングするときのエネルギー障壁を求めた結果，析出物中心よりわずかに前後にずれた位置でピークをとり，それぞれ 1.189[eV/nm]，1.222[eV/nm] であった．転位が析出物中心にあるときは 1.144[eV/nm] とわずかに低い値をとり，エネルギー極小値をとることが分かったが，その谷が小さいためらせん転位は析
出物にトラップされないものと結論づけられる．

74

第7章

結論

本研究では，第一原理計算により実際の Y2 O3 ,TiO2 の自由エネルギー曲線を求めた
後，Y と O，Ti と O を陽に区別しない平均粒子と近似したポテンシャル関数にフィッティングして，ODS 鋼中の酸化物分散形態のモンテカルロシミュレーション，分子動力学法による変形シミュレーション等を行った．また，原子弾性剛性係数 (AES) という基準を用いて変形中の個々の原子の力学状態について検討した．さらに，らせん転位と析出物の相互作用について，分子動力学シミュレーションだけではなく，NEB 法によるエネルギーバリアの定量評価も行った．以下に本論文を総括する.
第 2 章では本研究で用いた分子動力学法，モンテカルロ法，AES による力学評価手法，NEB 法について説明した．また，平均粒子のポテンシャルパラメータの決定方法についても述べた．
第 3 章では，モンテカルロ法により酸化物分散形態の検討を行った．バルクでの検討により YO 粒子は組成，温度に関わらず粒子数個程度のクラスターを形成し，シミュレーションセル内に均等に分布するのに対し，TiO2 粒子は比較的大きなクラスターを形成し，その寸法は YO 粒子の組成および温度条件によって変化することを明らかにした．温度 1000[K] と 1400[K] での TiO2 粒子の析出物寸法を比較すると，やはり 1400[K]の方が粗大化する傾向が認められたが，YO 粒子の組成が比較的高い Fe-0.6wt%YO-

0.2wt%TiO2 の組成では，1400[K] でも TiO2 粒子の粗大化が抑制されていた．また，その他の原子や空孔の影響も検討し，Al 原子の添加はもともと小さく分散する YO 粒子
をより凝集しにくくすること，原子空孔は YO 粒子の周囲に集中すること等を明らかにした．Σ5 傾角対称粒界近傍についての検討では，結晶粒界近傍に TiO2 粒子が特に顕著に粒界に偏析するが，YO 粒子にはそのような傾向は認められないことを明らか
にした．

75

第 4 章では，酸化物分散形態による力学応答変化を調べるために，モンテカルロ法
により作成した様々な組成の単結晶バルクを対象に，(1̄1̄2) 面をすべり面とし，[111]方向へ単純せん断させるシミュレーションを行った．酸化物が細かく均等に分散した

Fe-0.6wt%YO-0.2wt%TiO2 は，粗大 TiO2 析出物を有する他の組成に比べて高い弾性限界を示した．応力ピーク前後の原子構造を観察した結果，析出物界面近傍から双晶が発生し，ひずみの増加と共に互いに合体し成長することで応力急減しており，高い弾性限界を示した Fe-0.6wt%YO-0.2wt%TiO2 の組成では，TiO2 クラスターが小さく界面からの欠陥生成が抑制されていた．また，AES の正負で内部構造変化を観察した
結果，比較的粗大な TiO2 析出物を形成した系では，TiO2 析出物-母相界面で AES が負となる Fe 原子が集団的に発生し，そこが双晶変形の起点となっていた．一方，酸化物が細かく均等に分散し，高い弾性限界を示した Fe-0.6wt%YO-0.2wt%TiO2 の組成では，TiO2 の界面だけでなく YO 粒子周辺からも AES が負となる Fe 原子が様々な方向，
場所で発生，成長しており，負の AES となった Fe 原子の不安定モード，方向が異なるため，界面からの双晶発生が抑制されることを明らかにした．
酸化物が粒界での力学特性に与える影響を調べるため，第 5 章では，Σ5 ねじれ粒界を導入した系の [001] 方向への引張シミュレーションを行った．純鉄も含め，いずれの組成も εzz = 0.06〜0.07 近傍で第 1 応力ピークに達した後，εzz = 0.08 近傍で下げ止まり，εzz = 0.011 近傍まで同様の応答を示したが，以降は組成により応答が大きく異なった．純鉄および Fe-0.6wt%YO-0.4wt%TiO2 は第 2 応力ピークを示し粒界でへき開したが，その他の組成では第 2 ピークが不明確であり，特に Fe-0.2wt%YO-0.4wt%TiO2 では第 1 ピーク後の最大応力を示したあと，へき開した系に比較して段階的に応力が減少した．原子毎のポテンシャルエネルギー分布により構造の変化を調べた結果，いずれの組成も，第 1 応力ピーク前後に発生した帯状の分布が各結晶領域全体に伝播した後，追加発生したすべりによる変形が飽和したとき粒界のへき開または粒内き裂による破壊が起こった．Fe-0.2wt%YO-0.2wt%TiO2 ではすべりと粒界が交差する高エネルギー領域から純鉄と同様へき開が生じた．Fe-0.2wt%YO-0.4wt%TiO2 では粒界と結晶内部のすべり交差部を接続するような粒内き裂が発生した．Fe-0.6wt%YO-0.2wt%TiO2 で76

は完全に粒内からき裂が発生していた．Fe-0.6wt%YO-0.4wt%TiO2 では粒界に偏析した TiO2 とその周囲の Fe 原子の間でへき開を起こすと同時に粒内き裂も認められた．

Fe-0.2wt%YO-0.4wt%TiO2 と Fe-0.6wt%YO-0.4wt%TiO2 の粒界での TiO2 粒子の分布の違いからむらのある，一部密集した構造が粒界-粒内交差すべりへのき裂発生を促した可能性が示唆された．また，AES 負原子の分布により内部構造変化を調べた結果，

Fe-0.6wt%YO-0.4wt%TiO2 の組成において，引張後期でも最初の応力ピークで生じた
AES 負原子が解消されず，2 つの負の AES 原子面が 1 つになるように変化し，最初の応力ピーク後に導入された双晶欠陥間で結晶の回転または相変態を生じた可能性が示された．
第 6 章では，第 4 章のせん断シミュレーションで用いたセルを拡大して (1̄1̄2) 面をすべり面とするらせん転位を導入し，せん断ひずみを与えて転位を運動させるシミュ
レーションを行った．比較的大きな TiO2 析出物を持つ Fe-0.2wt%YO-0.4wt%TiO2 の系では，TiO2 析出物と転位の切り合いにより多数の欠陥を残した．その際，[111] 方
向において析出物を挟んですべり面が明確に変わる場合と変わらない場合があり，変わる場合は交差部を起点とした転位ループの形成が確認された．一方，すべり面を変えずに TiO2 析出物界面に欠陥を残したケースも観察された．これは界面の凹凸やミスフィットに起因するものと思われるがその詳細は不明である．酸化物が細かく分散する Fe-0.6wt%YO-0.2wt%TiO2 の組成でも，転位ループ形成が見られた．その際，転位と小さな析出物の切り合いから欠陥が生じ，析出物間で [111] 方向へ動いた後，[1̄10]方向へ欠陥が成長する過程が観察された．NEB 法によりらせん転位が TiO2 析出物(直径:2.89[nm]，析出物間距離:11.1[nm]) をカッティングするときのエネルギー障壁を求めた結果，析出物中心よりわずかに前後にずれた位置でピークをとり，それぞれ

1.189[eV/nm]，1.222[eV/nm] であった．転位が析出物中心にあるときは 1.144[eV/nm]とわずかに低い値をとり，エネルギー極小値をとることが分かったが，その谷が小さいためらせん転位は析出物にトラップされないものと結論づけられる．

ODS 鋼における原子シミュレーションでは，イオン結合，配位結合，金属結合を同時に表現しなければならないという困難がある．本研究ではそれを「平均粒子」とい77

う簡便なモデル化により回避した．このような大胆な近似に関わらず，以上の結果は実験観察で得られた知見に符合するばかりでなく，原子論を起点としたより詳細な検討を行うことに成功している．今後，高精度な原子間相互作用による検討や新たな実験観察結果から，本論文で得られた知見が検証・比較されることで新たな知見が得られることを期待する．

78

参考文献
(1) 木村晃彦，軽水炉高燃焼度化のための炉心材料に関する技術開発，革新的実用
原子力技術開発平成 14 年度成果報告書（2003）(2) Ukai,S. and Fujiwara,M., Journal of Nuclear Materials, 307 (2002), 749(3) Schaeublin, R. Leguey,T., Spatig,P., Baluc,N. and Victoria,M., Journal of Nuclear Materials, 307 (2002), 778(4) Hayashi,T., Sarosi,P.M., Schneibel,J.H. and Mills,M.J., Acta Materialia, 56(2008), 1407-1416(5) 鵜飼重治，ODS フェライト−マルテンサイト鋼被覆管の組織制御技術開発, サイクル機構技報No.7(2000)(6) 核燃料サイクル開発機構，株式会社コベルコ科研，残留α粒を有する高温強度に優れたマルテンサイト系酸化物分散強化型鋼の製造方法，特開 2005-76087，

2005-3-24(7) 京都大学，北海道大学，日本原子力研究開発機構，株式会社コベルコ科研，スーパー ODS 鋼，特開 2010-65302，2010-3-25(8) McClintock,D.A., Sokolov,M.A., Hoelzer,D.T. and Nanstad, R.K., Journal of
Nuclear Materials, 392 (2009), 353(9) Jitsukawa,S., Suzuki,K., Okubo,N., Ando,M. and Shiba,K., Nucl. Fusion, 49(2009), 115006(10) Zhang,C.H., Jang,J., Cho,H.D. and Yang,Y.T., Journal of Nuclear Materials,
386 (2009), 457(11) Kishimoto,H., Alinger,M.J., Odette,G.R. and Yamamoto,T, Journal of Nuclear
Materials, 329 (2004), 369(12) Miller,M.K., Kenik,E.A., Russell,K.F., Heatherly,L., Hoelzer,D.T. and Maziasz,P.J., Materials Science and Engineering A, 353 (2003), 140-14579

(13) Wasilkowaska,A., Bartsch,M., Messerschmidt,U., Herzog,R. and CzyrskaFilemonowicz,A., Journal of Materials Procesing Technology, 133 (2003), 218224(14) Ukaia,S., Ohtsuka,S., Kaito,T., Sakasegawa,H., Chikata,N., Hayashi,S. andOhnuki,S., Materials Science and Engineering A, 510-511 (2009), 115-120(15) Yamamoto,M., Ukai,S., Hayashi,S., Kaito,T. and Ohtsuka,S., Materials Scienceand Engineering A, 527 (2010), 4418-4423(16) Sakasegawa,H., Ukai,S., Tamura,M., Ohtsuka,S., Tanigawa,H., Ogiwara,H., Kohyama,A. and Fujiwara,M., Journal of Nuclear Materials, 373 (2008), 82-89(17) Ohtsuka,S., Ukai,S., Fujiwara,M., Kaito,T. and Narita,T., Journal of Nuclear
Materials, 329 (2004), 372(18) Shinozuka,K., Tamura,M., Esaka,H., Tanigawa,H. and Shiba,K., Journal of theSociety of Materials Science, Japan, 59-8 (2010), 637-644(19) 成田健，鵜飼重治，皆藤威二，大塚智史，小林十思美，再結晶組織を有する ODS
フェライト鋼被覆管の二段軟化熱処理技術開発，サイクル機構技報

No.18(2003)(20) Ohtsuka, S., Ukai,S., Sakasegawa,H., Fujiwara,M., Kaito,T. and Narita,T., Journal of Nuclear Materials, 367-370 (2007), 160-165(21) Pareige,P., Miller,M.K., Stoller,R.E., Hoelzer,D.T., Cadel,E. and Radiguet,B.,
Journal of Nuclear Materials, 360 (2007), 136(22) Kim, I.-S., Hunn,J.D., Hashimoto,N., Larson,D.L., Maziasz,P.J., Miyahara,K.
and Lee,E.H., Journal of Nuclear Materials, 280 (2000), 264(23) Kimura,A., Sugano,R., Matsushita,Y. and Ukai,S., Journal of Physics andChemistry of Solids, 66 (2005), 504-508(24) Ohnuma,M., Suzuki,J., Ohtsuka,S., Kim,S.-W., Kaito,T., Inoue,M. and Kitazawa,H., Acta Materialia, 57 (2009), 5571
(25) Alinger,M.J., Odette,G.R. and Hoelzer,D.T., Acta Materialia, 57 (2009), 392(26) Marquis,E.A., APPLIED PHYSICS LETTERS, 93 (2008), 181904
80(27) Klimiankou,M., Lindau,R. and Moslang,A., Micron, 36 (2005), 1-8(28) Klimiankou,M., Lindau,R. and Moslang,A., Journal of Crystal Growth, 249(2003), 381-387(29) Hirata,A., Fujita,T., Wen,Y.R., Schneibel,J.H., Liu,C.T. and Chen,M.W., Nature Materials, 10 (2011), 922-926(30) Narita,T., Ukai,S., Kaito,T., Ohtsuka,S. and Matsuda,Y., Journal of Nuclear
Science and Technology, 45-2 (2008), 99-102(31) Gopejenko,A., Zhukovskii,Y.F., Vladimirov,P.V., Kotomin,E.A. and Moslang,A,Journal of Nuclear Materials, 406 (2010), 345-350(32) Gopejenko,A., Zhukovskii,Y.F., Vladimirov,P.V., Kotomin,E.A. and Moslang,A,Journal of Nuclear Materials, 416 (2011), 40-44(33) Jiang,Y., Smith,J.R. and Odette,G.R., PHYSICAL REVIEW B, 79 (2009),
064103(34) Barnard,L., Odette,G.R., Szlufarska,I. and Morgan,D., Acta Materialia, 60(2012), 935-947(35) 山本智，神戸大学修士論文（2010）(36) Bako,B., Weygand,D., Samaras,M., Chen,J., Pouchon,M.A., Gumbsch,P. and
Hoﬀelner,W., Philosophical Magazine, 87-24 (2007), 3645-3656(37) Bako,B., Weygand,D., Samaras,M., Hoﬀelner,W. and Zaiser,M., Physical Review B, 78 (2008), 144104(38) Bako,B., Zaiser,M., Weygand,D., Samaras,M. and Hoﬀelner,W., Journal of Nuclear Materials, 385 (2009), 284-287(39) Takahashi,A., Chen,Z., Ghoniem,N. and Kioussis,N., Journal of Nuclear Materials, 417 (2011), 1098-1101(40) Robertson,C. and Gururaj,K., Journal of Nuclear Materials, 415 (2011), 167178

81

(41) Gururaj,K. and Robertson,C., Energy Procedia, 7 (2011), 279-285(42) Alinger,M.J., Wirth,B.D., Lee,H.J. and Odette,G.R., Journal of Nuclear Materials, 367 (2007), 153.
(43) Hin,C., Wirth,B.D. and Neaton3,J.B., Physical Review B, 80 (2009), 134118(44) 睦門賢憲，神戸大学修士論文（2011）(45) Yashiro,K., Yamaguchi,A., Tanaka,M., Okuda,T., Koga,K. and Segi,T., ”Molecular Dynamics Simulation on Interaction between Screw Dislocation and Pseudo
Yttrium Oxide in Bcc-Fe”, Materials Transactions, 53-2 (2011), in press(46) 上田顯, コンピューターシミュレーション, (1990), 朝倉書店
(47) 上田顯, 分子シミュレーション-古典系から量子系手法まで-, (2003), 裳華房(48) 北川浩, 渋谷陽二, 北村隆行, 中谷彰宏, 初心者のための分子動力学法, (1997), 養賢堂(49) 寺元英, マルコフ過程と力学過程, (1961), 槇書店(50) Hastings,W.K., Biometrika, 57 (1970), 97-109.
(51) Lill,J.V. and Broughton,J.Q., Physical Review B, 49(17) (1994), 11619.
(52) 松本眞, 擬似乱数ユーザーの方へ-Mersenne Twister 法開発者より, 日本統計学
会誌, 35(2) (2006), 165(53) Matsumoto,M. and Nishimura,T., ACM Transactions on Modeling and Computer Simulation, 8(1) (1998), 3(54) 二宮祥一, 松本眞さんの疑似乱数の発生法について, 数学通信, 10(2) (2005), 59(55) Saito,M. and Matsumoto,M., SIMD-oriented Fast Mersenne Twister: a 128-bit
Pseudorandom Number Generator,Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, (2008), 607(56) Wallace,D.C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
(57) Wang, J., Li, J., Yip, S., Phillpot, S. R. and Wolf, D., Mechanical Instabilities
of Homogeneous Crystals Physical Review, B, 52(17), pp.12627–12635, (1995).
82

(58) Kresse,G. and Furthmuller,J., Physical Review B, 54 (1996), 11169.
(59) 岩崎富生, 佐々木直哉, 保川彰夫, 千葉矩正, 日本機械学会論文集 (A 編), 62(597)(1996), 1242.
(60) Johnson,R.A., Physical Review, 134 (1964), A1329.
(61) Li,J., Modelling Simul. Mater. Sci. Eng., 11 (2003), 173(62) 日本金属学会, 転位論―その金属学への応用―, (1971), 丸善(63) 日本金属学会, 材料強度の原子論, (1985), 丸善(64) 屋代如月, 藤原正大, 材料, 59(8) (2010), 631-636.

83

付録 A

関連学術論文・講演

関連学術論文 K. Yashiro, A. Yamaguchi, M. Tanaka, T. Okuda, K. Koga, T. Segi, Molecular Dynamics Simulation on Interaction between Screw Dislocation and Pseudo
Yttrium Oxide in Bcc-Fe, Materials Transactions, 53(2)(2011), in press.

関連学術講演 山口 明宏，田仲 稔，睦門 賢憲，屋代 如月，ＯＤＳ鋼中の酸化物分散形態
予測：モンテカルロシミュレーション，日本材料学会第 15 回分子動力学シンポジウム，札幌コンベンションセンター (2010.5)．

 山口 明宏，田仲 稔，奥田 隆成，古賀 健治，世木 隆，屋代 如月，ＯＤ
Ｓ鋼中における結晶粒界と酸化物分散形態：モンテカルロシミュレーション，日本機械学会第 23 回計算力学講演会，北見工業大学 (2010.9)．

 山口 明宏，田仲 稔，奥田 隆成，古賀 健治，世木 隆，屋代 如月，ＯＤ
Ｓ鋼中の格子欠陥が酸化物分散形態に及ぼす影響：モンテカルロシミュレーショ
ン，第 1 回マルチスケールマテリアルモデリングシンポジウム，大阪大学コンベ
ンションセンター (2011.5)．

 K. Yashiro, T. Mutsukado, M. Tanaka, A. Yamaguchi, K. Koga, T. Segi andT. Okuda, Molecular Dynamics Simulations on Interaction between Dislocation
and Y2O3 Nanocluster in Fe, XI International Conference on Computational
Plasticity (COMPLAS XI), Barcelona, Spain,(2011)．

84 山口 明宏，屋代 如月，古賀 健治，世木 隆，奥田 隆成，NEB 法による
ODS 鋼中の転位反応経路解析，日本機械学会第 24 回計算力学講演会，岡山大学
津島キャンパス (2011.10)．

 K. Yashiro, A. Yamaguchi, K. Koga, T. Segi and T. Okuda, Molecular Dynamics
Simulation, Local Lattice Instability Analysis and Nudged Elastic Band Calculations on Interaction between Screw Dislocation and Pseudo-Y2O3 Nanoparticle
in Bcc-Fe, Plasticity 2012, San Juan, PR, USA,(2012)．

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

謝辞
本論文は，著者が学部 4 年生で研究室に所属して以来，博士前期課程を修了するまでに得られた研究成果をまとめたものです．研究全般に渡り，終始懇切丁寧にご指導いただきました表面・界面工学研究室 屋代如月准教授に深謝いたします．浅学な著者が本論文の執筆を完遂するに至ったのは，ひとえに，先生が私の自主性を尊重しつつ，真摯な研究姿勢により導いてくださったお陰であり，先生の下で研究を行えたことは
著者のこの上ない幸運でありました．また，広い知識と鋭い洞察により研究に関して非常に有益なご助言をいただくとともに，日々の研究活動が円滑となるようご支援い
ただいた同研究室の田中克志教授に心より感謝いたします．
本研究を遂行するにあたり，研究環境の整備および第一原理計算や実験観察の分野で株式会社コベルコ科研の古賀健治氏，世木隆氏，奥田隆成氏には多大なるご支援ご
協力を承りました．厚く御礼申し上げます．また，本論文の完成にあたって，転位論
や金属材料に関してご助言いただきました固体力学研究室 長谷部忠司准教授に深く感謝いたします．
研究を進めるにあたり，ご支援ご協力いただいた表面・界面工学研究室および固体力学研究室の諸氏に心より感謝の意を表します．特に，同じ研究グループにおいて，睦門賢憲氏 (株式会社神戸製鋼所) には原子レベルでの転位の取扱いで，田仲稔氏には第一原理計算で多岐に渡ってご協力いただきました．また，表面・界面工学研究室の谷口正典氏には分子動力学法で多くのご助言をいただきました．ありがとうございます．
最後に，長きに渡る学生生活において深い理解を以って支えてくれた両親に最大限の感謝の意を表します．

平成 24 年 2 月

117

山口

明宏