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

(1)

修士論文

原子モデルによる鋼の酸化物分散強化機構の検討：

分散形態ならびに転位芯との相互作用

指導教員：屋代如月

山口明宏

2012

年

2 月

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

(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

(3)

要約

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にトラップされないことなどを明らかにした．

(4)

2

Summary

For a new insight on the strengthening mechanism in ODS steels, we have performed Monte Carlo (MC), molecular dynamics (MD) simulations and NEB calculations with the simple pair potential functions for monatomic “pseudo” atoms of Y2O3 and TiO2, that is, each Y, Ti and O atom isn’t distinguished in these simulations but their pa-rameters are fitted to correctly represent the lattice length and bulk modulus of real Y2O3 and TiO2 according to the ab initio DFT results. We also discussed the defor-mation mechanism with the positive definiteness of the atomic elastic stiffness (AES), which indicates the lose of local resistance against deformation or “local softening atom”. First, MC simulations are performed to discuss the effect of atom composi-tion, 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 in the 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 ag-gregation 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 observation of 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 the initiation 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 grow at the coarse TiO2 interface, while the negative AES domains are rather small and have different 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 effect of grain boundary segregation. The first stress drop is caused by slip in the grains, so that there is little difference in pure Fe and four different 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 and alloys 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 de-pending on its composition. Finally, we have set a screw dislocation in two alloys and glide it by applying shear strain, resulting in the following observations, e.g., formation of 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 of energy barrier are evaluated as 1.189[eV/nm] and 1.222[eV/nm], not just at the center of TiO2 but slightly deviated positions. Thus the energy barrier shows local minimum of 1.144[eV/nm] at the center position; however, the energy difference is so small that we can conclude the screw dislocation would not be trapped by the TiO2 cluster.

(5)

1 緒論 1 1.1 ODS鋼の開発状況 . . . . 1 1.2 ODS鋼の数値シミュレーション . . . . 3 1.3 研究目的と各章の構成 . . . . 4 2 解析手法 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

(6)

目次 ii 3 モンテカルロ法による酸化物分散形態の検討 22 3.1 解析条件 . . . . 22 3.1.1 バルクにおける検討 . . . . 22 3.1.2 結晶粒界における検討 . . . . 24 3.2 解析結果 . . . . 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 3.3 結言 . . . . 35 4 単結晶のせん断変形シミュレーション 37 4.1 解析条件 . . . . 37 4.2 解析結果 . . . . 38 4.2.1 応力-ひずみ応答と欠陥発生の様相 . . . . 38 4.2.2 系の弾性剛性係数 . . . . 41 4.2.3 AESの評価 . . . . 42 4.3 結言 . . . . 45 5 ねじれ粒界のへき開シミュレーション 47 5.1 解析条件 . . . . 47 5.2 解析結果 . . . . 49 5.2.1 応力-ひずみ応答 . . . . 49 5.2.2 変形破壊の様相 . . . . 50 5.2.3 系の弾性剛性係数と AES の評価 . . . . 55 5.3 結言 . . . . 59

(7)

目次 iii 6 らせん転位と析出物の相互作用 62 6.1 解析条件 . . . . 62 6.1.1 分子動力学シミュレーション . . . . 62 6.1.2 NEB法によるエネルギーバリア評価 . . . . 64 6.2 解析結果 . . . . 65 6.2.1 応力-ひずみ応答および弾性剛性係数 . . . . 65 6.2.2 転位の挙動 (粗大析出物) . . . . 66 6.2.3 転位の挙動 (微細析出物) . . . . 70 6.2.4 NEB法によるカッティングバリア評価 . . . . 72 6.3 結言 . . . . 73 7 結論 75 参考文献 79 A 関連学術論文・講演 84 謝辞 117

(8)

第

1

_{章緒論}

1.1 ODS

鋼の開発状況

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

Dispersion Strengthened: ODS)鋼の開発が進められてきた．ODS 鋼はメカニカルア

ロイング (Mechanical Alloying: MA) によって製造され，酸化イットリウム (Y2O3)等を 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)．この残留フェライトが優れた高温クリープ強度をもたらしており，その主要な変形モードは比較的弱い結晶粒界や焼きならし時のマルテンサイト変態によって生じるパケット境界でのすべりである

(9)

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

(10)

1.2 ODS

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

このように ODS 鋼の開発は強化メカニズムの考察を伴って順調な進展を見せている．以上はすべて実験観察によるものであり，強化メカニズムの考察は一般的な金属における知見を援用したものである．これら考察はスーパー ODS 鋼の開発につながったことからも妥当かつ極めて有用であるが，ナノクラスターの詳細な構造や転位のダイナミクス，粒界のすべりの形態のようなナノ・メゾスケールの知見において未知の部分を多く残している．このような詳細な検討は実験観察と理論のみの考察では限界に近づいており，転位・原子の詳細なダイナミクスやエネルギーバリアを扱うことのできるシミュレーション (数値実験) による検討が期待されている．密度汎関数理論に基づき電子状態からエネルギーや応力等を計算する第一原理計算は何ら実験結果を必要とせずに高い精度で検討できるが，膨大な計算量から数十から数百原子程度の取扱いに限られる．そのため ODS 鋼関連では主として Y,Ti,O などの元素の固溶形態の検討に用いられてきた(31)−(34)_{．その結果として酸素原子が添加原子間} の強力な結合を与えるとともに Y2O3ナノクラスターが成長するのに必要となる空孔の集中を促す等の知見を得ている(32)．我々の研究グループでも，bcc-Fe 中に Y,Ti,O を添加しての安定性を検討し(35)，現在はそれら構造において格子不安定性を用いた力学的評価を行っている．また，すでに述べたように転位とナノクラスターの相互作用は高温クリープ強度等に影響する重要な要素であり，離散転位動力学 (Discrete Dislocation Dynamics: DDD)を用いた検討も行われている(36)−(41)_{．DDD は任意の曲線形状を有} する転位線を多数の直線セグメントに離散化し，転位のダイナミクスを追跡する手法である．DDD を用いた検討から臨界分解せん断応力は Y2O3クラスターと母相との格子ミスマッチに強く依存する(39)_{等の知見が得られている．また，照射の影響も議論} され，ナノクラスターに転位が堆積し内部応力が上昇することで交差すべりした転位が照射欠陥を越えるとの考察もある(40)．DDD による検討は熱揺動の影響を考慮でき，多数の転位を取扱うことができるという点で大変有用であるが，基本的に転位論と連続体力学によって定式化され，ナノクラスターと転位芯の相互作用もモデル化を通し

(11)

て考慮する必要があり，原子レベルでの新たな知見を得ることは難しい．

第一原理計算と 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 法で取扱う際，ポテンシャル関数でイオン結合，配位結合，金属結合を同時に表現しなければならないという困難が背景にある．そこで先行研究では，第一原理計算で求めた Y2O3のエネルギー曲線を，Y と O を陽に区別せず平均化した仮想的な「YO 粒子」によってモデル化し，転位と YO ナノクラスターとの相互作用の MDシミュレーションを行ってきた(44)(45)．

1.3 研究目的と各章の構成

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

(12)

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

(13)

第

2

_{章解析手法}

2.1 分子動力学法

分子動力学法 (Molecular Dynamics Method: MD) は，系を構成する個々の原子についてニュートンの運動方程式 mαd 2_rα dt2 = F α (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)± ∆tdr α dt + (∆t)2 2 d2_rα dt2 ± (∆t)3 3! d3_rα dt3 + O((∆t) 4₎ _(2.3) となる．両式の和をとり式 (2.1) を代入すると rα(t + ∆t) = 2rα(t)− rα(t− ∆t) + (∆t)2F α (t) m + O((∆t) 4₎ _(2.4) を得る．これより，時刻 t− ∆t と t における全原子の位置を既知として，時刻 t + ∆t における任意の原子 α の位置を求めることができる．

(14)

分子動力学法で温度制御する場合，最も簡単で直接的な方法として速度スケーリング 法がよく用いられる．熱統計力学より系の運動エネルギー K は次のように表される． K = 1 2 N ∑ α=1 mαvα·vα = 3 2N kBT (2.5) ここで，mα_{は原子 α の質量，v}α_{は原子 α の速度，N は系の全原子数，k} Bはボルツマ ン定数，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)2F α_(t) m (2.7) を √ TC/T ∆rα(t + ∆t)で置き換えることに相当する．平衡状態では，能勢の方法など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている(48)_．

2.2 モンテカルロ法

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

(15)

2.2.1 マルコフ過程と

Metropolis-Hasting

法

1つ手前の状態から現在の状態が確率的に決定される確率過程をマルコフ過程とい う．対象とするマルコフ過程を時間に依存しない定常であるとすると，状態 i から状 態 j に n 回の試行で遷移する確率 P_ij(n)は P_ij(n)=∑ k P_ik(n−1)Pkj (2.8) で表される．ここで Pkj = P (1) kj である. 式 (2.8) によって生成される状態の連なりのことをマルコフ連鎖と呼ぶ．マルコフ過程の任意の状態から任意の状態への遷移が起こ り得るとして，このときマルコフ過程は既約であるという．また，状態 j から出発し て l の試行回数で初めて状態 j に再帰するときの確率を fj(l) で表す． ∞ ∑ l=1 f_j(l)= 1 (2.9) が成立するとき状態 j は再帰確実な状態であるという．平均再帰時間を µj = ∞ ∑ l=1 lf_j(l) (2.10) で定義する．また，同じ試行回数毎に再帰するときその状態を周期的な状態という. マルコフ過程の全状態が再帰確実で非周期的かつ平均再帰時間が有限である状態のときマルコフ過程はエルゴード的と呼ばれる．既約でエルゴード的なマルコフ過程におい て任意の状態 j で lim n→∞P (n) ij = wj > 0 (2.11) となる定常分布{wj = 1/µj} が存在し， wj = ∑ i wiPij (2.12) が成立する(49)_{．式 (2.12) を書き改めて} ∑ i̸=j wiPij = ∑ i̸=j wjPji (2.13)

(16)

を得る．この十分条件として wiPij = wjPji (2.14) で表される詳細釣合い条件式を得る．式 (2.14) を用いて任意の確率分布を実現する方法は Metropolis-Hasting 法(50)_と呼ばれる．扱う統計集団が既約でエルゴード的，すなわち任意の 2 つの状態は互いに平均再帰時間が有限で到達可能であり周期的でないとする．詳細釣合い条件式 (2.14) か ら状態 i から状態 j への遷移確率を Pij = { wj/wi (wj <_{= w}i) 1 (wj > wi) (2.15) とすると，確率過程の極限分布として指定した分布{wj} を十分条件として実現できる．

2.2.2 NPT

アンサンブル

Lillらは系の変形によるエネルギーを明示的にポテンシャルエネルギーから分離することで系全体の変形を考慮する圧力一定のモンテカルロ法を定式化した(51)．Lill ら の方法に沿った粒子数 N ，圧力 P ，温度 T 一定の系のモンテカルロ法の概要を以下に 示す．シミュレーションセルを平行六面体として 1 つの頂点を直交座標系の原点とする．平 行六面体の基本ベクトルを基底としたセル内座標系 siと直交座標系 riとの変換行列を aij = rits−1j ，変形の基準となる変換行列 aijを特に bij，変形勾配テンソルを Jij = aikb−1kj， Greenのひずみテンソルを ηij = 1₂(JkiJkj − δij)，指定する真応力テンソルを σij，第 2piola-kirchhoﬀ応力テンソルを τij = |J|(Jik−1σklJjl−1)とする．また，粒子 n の位置を r(n)_i ，速度を q_i(n)で表し drN _{= dr}(1) i . . . dr (N ) i , dqN = dq (1) i . . . dq (N ) i ，シミュレーション セルのハミルトニアンを H(r(n)_i , q_i(n))とすると，N, σij, T 一定のアンサンブルの分配関数は ΞN σT = 1 N !h3N p ∫∫∫ exp  ₋H(r(n)i , q (n) i )− |b|τijηij kBT  _dqN_drN_dJ ij (2.16)

(17)

で表される．ただし，hpはプランク定数，kBはボルツマン定数である．ポテンシャル エネルギーを U (ri(n))とすると式 (2.16) より微小体積要素 drNdJijに系を見出す確率は exp  ₋U (r (n) i )− |b|τijηij kBT  _drN dJij = exp  ₋U (s (n) i )− |b|τijηij kBT  _|a_ij_|N dsNdJij (2.17) に比例する．ここで式 (2.17) の右辺は左辺を変数変換して得られるものである．W = U − |b|τijηij − NkBT ln|aij| を定義し, 式 (2.17) を式 (2.15) に当てはめると，十分条件 として状態 l から状態 m への遷移確率 Plm = { 1 (Wm <= Wl) exp{−(Wm− Wl)/(kBT )} (Wm > Wl) (2.18) を得る．

2.2.3 擬似乱数発生器

モンテカルロシミュレーションにおいて使用する一様乱数は周期が長く，高次元に対して一様均等分布することが望ましい．そのような乱数を得る手段として熱振動や原子核分裂などの本質的に等確率で起こるされている自然現象を利用する物理乱数が挙げられる．しかし，物理乱数は高価であったり速度が遅かったりする問題がある．さらに，シミュレーションの再現性を確保するには使用した乱数をすべて記録しておく必要がある．それに対して，一様性と均等分布性を合わせ持つ数列を真の乱数の代わりに使うことがモンテカルロシミュレーションでは多い．このような数列を擬似乱数といい，発生させるアルゴリズムのことを擬似乱数発生器という．現在，最も普及している擬似乱数発生器の 1 つは Lehmer による線形合同法であり，次の線形漸化式により表される． xn+1 = axn+ b (mod M ) (2.19) この擬似乱数は，fortran やＣ言語の標準乱数として組込み関数で与えられている場合 が多く扱いやすい．しかし，法を M とするため発生する擬似乱数の周期は高々M で あるという問題点を持つ．また，高次元空間でのプロットに対する欠陥も報告されている(52)_．

(18)

そこで，本研究では Mersenne twister 法(53)_{を用いる．今，x，a が要素が 0 または} 1の w 次元の行ベクトルすると漸化式は次式で与えられる． xk+n = xk+m⊕ (xuk|xlk+1)A, (k = 0, 1, . . .) (2.20) ここで，⊕ は排他的論理和を表し，n は繰り返しを表す整数，m は 0 <_{= m <}_{= n となる} 整数である．(xu k|xlk+1)は xkの上位 (w− r)bit と xk+1の下位 rbit を並べて得られる行 ベクトルで 0 <_{= r <}_{= w − 1 を満たす．また，行列 A は次式で与えられる．} A =           1 1 . .. 1 aw−1 aw−2 · · · a0           (2.21) 本研究では 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を導入し，その半径

(19)

内の原子からの寄与のみを考慮すればよい．そこで本研究では領域分割法を用いてい

る．領域分割法では系を間隔が rc程度である格子状に分割し，対象とする原子の存在

する領域 (図 2.1 着色部) および隣接する領域 (図 2.1 斜線部) の原子の内，距離が rc以

下であるものを探索する．

また，本研究では分子動力学法においてメモリ共有環境を対象とした並列化を行っている．Message Passing Interface(MPI) を用いたメモリ分散環境を対象とする並列化は，マシンに搭載する CPU 数に捕らわれないが，分散したメモリ間での通信が必要となり，アルゴリズム構築に手間がかかる．一方，OpenMP を用いたメモリ共有環境を対象とする並列化は通信の必要がなく，逐次プログラムで並列化させる部分を指定するのみで済む．領域分割法により原子間相互作用を計算する際，領域毎に独立に計算することができるため，その部分を並列化した．

x

y

0 bx

by

(20)

2.4 原子弾性剛性

(Atomic Elastic Stiﬀness: AES)

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

応力 σij および弾性係数 Cijkl は，等温過程では σij = 1 V ( ∂F ∂ηij ) , Cijkl = 1 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 らは，ひずみの対称性を考慮したテンソル

B_ijklsym ≡ (Bijkl+ Blkji)/2の正値性によって結晶の安定性を評価することを提案してい

る(57)_{．すなわち，変形抵抗の喪失点を不安定の物理的描像としている．B}sym ijkl の行列式の正負によって格子不安定性が評価される．

2.4.2 原子系への拡張

局所の安定性を評価するための原子弾性剛性係数 Bijklα の算出に必要な原子応力 σαij ならびに原子弾性係数 Cijklα は，各原子周りの微小ひずみに対するポテンシャルエネルギーの 1 次，2 次変化量として導出される．原子応力 簡単のため，結晶の内部エネルギー U が系全体のポテンシャルエネルギー Etot に 等しいとする．このとき，応力は平衡状態からの微小ひずみ η に対するポテンシャル

(21)

エネルギーの単位体積当たりの変化として与えられる (56)_． σij = 1 V ∂Etot ∂ηij (2.24) ここで，V は平衡状態における系の体積であり，下付添字のローマ文字はテンソルの デカルト座標成分を表す．(2.24) 式の微分を求めるため，平衡状態からの仮想的な均 一変形を考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によって rα= J ¯rα (2.25) と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子 α と原子 β の間の距離 rαβ _には (rαβ)2 = ¯r_iαβGijr¯αβj (2.26) なる関係が成立する．ただし，Gij = JkiJkj である．仮想変形の Lagrange ひずみテン ソル ηij は ηij = 1 2 [ Gij − δij ] (2.27) であり，その微小量 dηij = 1 2dGij (2.28) と式 (2.26) の関係から次の関係が得られる． ∂rαβ ∂ηij = r¯ αβ i r¯ αβ j rαβ (2.29) これより 2 体ポテンシャル ϕ(rαβ₎_{における応力は次式で評価される} σij = 1 V ∂Etot ∂ηij = 1 V (₁ 2 N ∑ α N ∑ β(̸=α) ∂rαβ ∂ηij ∂Etot ∂rαβ ) = 1 2V N ∑ α N ∑ β(̸=α) ϕ′(rαβ)r αβ i r αβ j rαβ ここで，各原子位置における原子応力を σ_ijα = 1 V /N N ∑ β(̸=α) 1 2ϕ ′_(rαβ )r αβ i r αβ j rαβ (2.30)

(22)

と定義すると，系の応力は σij = 1 N N ∑ α σ_ijα (2.31) となる．原子弾性係数 弾性係数も応力と同様に U ≈ Etot の場合には Cijkl= 1 V ∂2Etot ∂ηij∂ηkl (2.32) であるので，平衡状態からの仮想均一変形を考えると 2 体ポテンシャルにおける弾性係数は以下のようになる． Cijkl = 1 2V N ∑ α N ∑ β(̸=α) ∂rαβ ∂ηkl ∂ ∂rαβ (_∑N α N ∑ γ(̸=α) 1 2ϕ ′_(rαγ )r αγ i r αγ j rαγ ) = 1 V [ 1 2 N ∑ α N ∑ β(̸=α) { ϕ′′(rαβ)− ϕ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 ] (2.33) 応力と同様に，各原子位置における原子弾性係数を以下のように定義する． C_ijklα = 1 V /N [₁ 2 N ∑ β(̸=α) { ϕ′′(rαβ)−ϕ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 ] (2.34) これより，系の弾性係数は Cijkl= 1 N N ∑ α C_ijklα (2.35) のように原子弾性係数の平均となる．原子弾性剛性係数以上で定義した原子応力，弾性係数から，原子弾性剛性係数 (AES) は以下で評価できる． Bα_ijkl= C_ijklα + (σα_ilδjk+ σjlαδik+ σikαδjl+ σjkαδil− 2σijαδkl)/2 (2.36)

Wangらによる提案(57)_{に従い，Voigt 対称性をもたせた B}α sym

ijkl ≡ (Bijklα + Blkjiα )/2 を

用いて安定性評価を行う．以降では B_ijklα sym を Voigt 表記した Bα

IJ を原子弾性剛性係数

(23)

2.5 Nudged Elastic Band (NEB)

法

Nudged Elastic Band(NEB)法は，初期状態と最終状態の 2 つの状態を結ぶ反応経

路の中で，エネルギーが最も低くなる経路を探索する手法である．N 個の原子から成 る系の状態は，すべての原子の位置座標を軸とする 3N 次元の位相空間において，位 置ベクトル ⃗R= (x1, y1, z1,· · ·) で表される．位相空間における初期状態 ⃗Rsおよび最終 状態 ⃗Reの 2 点を結ぶ経路を複数の位相点に分割し，各位相点を仮想的なばねでつないだ状態で，位相点にかかる力に基づいて位相空間内における低いエネルギー経路を探 索する (図 2.2)．各位相点に作用する力 ⃗F0 i は，経路に垂直なポテンシャルエネルギーの勾配と経路に平行なバネ力の和として次式で表す． ⃗ F_i0 =−⃗▽V ( ⃗Ri)|⊥+ ( ⃗Fis・ˆτ∥)ˆτ∥ (2.37) ここで，V はポテンシャルエネルギー， ⃗Riは i 番目の位相点の位置ベクトル， ⃗Fisは i 番目の位相点に作用するバネ力，ˆτ_∥は経路の接線方向単位ベクトルである．−⃗▽V ( ⃗Ri) は状態 ⃗Riの原子配置において，各原子に働く力の 3N 個の成分そのものである．位 相空間内で，ばねに平行な方向でエネルギーが急激に変化する場合，式 (2.37) では仮 想的なばねの力のみが支配的になる．そこで隣接する角度によって変化する関数 f (ϕi) を導入し，次式で各位相点に働く力を表現することが提案されている． Energy state Start

End _{Energy state}

Start

End

(a) Without spring (b) With spring

(24)

⃗ F_iN EB = ⃗F_i0 + f (ϕi)( ⃗Fis− ( ⃗Fis・ˆτ∥)ˆτ∥) (2.38) f (ϕi) = 1 2(1 + cos(π(cos ϕi))) (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を変化させ

(25)

ルギーを求め，各原子の個数分を単位格子の持つ全自由エネルギーから差し引くことで単位格子あたりの結合エネルギーを求めた (図 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 同士の相互作用も Morse potentialで評価するものとし，ポテンシャルパラメータは岩崎らの文献(59) を参照した．表 2.3 に引用したパラメータを示す．ただし，エネルギーの単位は原論文から変換して表示している． Fe同士の相互作用は次式の Johnson potential(60) _{を採用した．} ϕij =−C1(rij − C2)3+ C3rij− C4 [eV] (2.41) パラメータ C1 ∼ C4の値を表 2.4 に示す．このポテンシャルは，bcc 構造をとる α-Fe のために作られた経験的 2 体ポテンシャルである．また，異種間粒子の場合は，それぞれの粒子が同種の組み合わせのときにもつ結合あたりのポテンシャルエネルギーを斥力項と引力項にわけ，それぞれを相乗平均をとることによって異種間の組み合わせの結合がもつポテンシャルエネルギーとした.

Table 2.1 DFT calculation conditons.

YO TiO2 Y Ti O

Number of ions 2 6 1 1 1

Number of electrons 16 32 3 4 6

Number of bands 12 28 4 4 5

Cutoﬀ energy [eV] 296.9 495.0 148.9 226.5 495.0 Number of k points 153 ₆3 ₇3 ₇3 ₇3

(26)

Ti Y O O l_a la la 0.644la l_a 1.633la

(a) hcp unit lattice for YO. (b) unit lattice for TiO .2

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

-8.2 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30 Energy [eV] Lattice parameter [nm]la VASP Morse potential

Fig.2.4 YO potential fitting.

-40 -35 -30 -25 -20 -15 -10 -5 0 5 VASP n=8 n=9 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 Energy [eV] Lattice parameter [nm]la

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.

De [eV] A [nm−1] r0 [nm]

Al-Al 0.1193058 23.53643 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

(27)

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.

-2.80 0.15470 5.21213 2.40528 9.54316

2.80-3.68 0.41882 3.73824 0.81108 3.25410

3.68-4.29 0.20371 2.90755 1.16981 4.48026

Table 2.6 Johnson potential parameters for TiO2-TiO2.

-2.52 3.24926 2.91738 2.93465 9.39188

2.52-2.95 2.65600 2.86962 2.36968 7.87776

2.95-3.66 0.34211 1.71521 3.88863 11.7160

Table 2.7 Johnson potential parameters for Fe-YO.

-2.38 1.15840 3.25778 1.68731 5.04371

2.38-3.46 0.56661 3.29399 0.42713 1.69365

3.46-4.26 0.25457 3.55386 0.38527 1.55162

Table 2.8 Johnson potential parameters for Fe-TiO2.

-2.43 2.73709 2.98410 2.69997 8.07629

2.43-3.11 1.72268 2.91421 1.37733 4.59222

3.11-3.64 0.12507 0.40886 3.91689 10.0382

Table 2.9 Johnson potential parameters for YO-TiO2.

-2.94 1.89672 2.93498 1.14839 4.15864

2.94-3.63 0.47657 2.92525 1.14862 4.15929

(28)

以下，モンテカルロシミュレーションには Morse potential，分子動力学シミュレーションには Johnson potential を使用する．参考に，モンテカルロ法，分子動力学法における同種粒子間の 2 体ポテンシャル曲線をそれぞれ図 2.6，図 2.7 に示す． -2 -1 0 1 2 3 4 5 6 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50

Pair potential energy [eV]

Distance [nm]

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

Fig.2.6 Pair potential energy curves for Monte Carlo simulations.

-2 -1 0 1 2 3 4 5 6 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50

Pair potential energy [eV]

Distance [nm]

Fe-Fe YO-YO TiO2-TiO2

(29)

第

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の系においては原子空孔を導入した場合も検討した．原子空孔を

(30)

粒子」を新たに定義し，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 250000

Temperature [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%YO-0.2wt%TiO2-1.0wt%Al alloy.

Cell structure 42× 42 × 42 bcc unit lattices

Temperature [K] 1000,1400

Table 3.3 Condition of Monte Carlo simulations for bulk Fe-0.4wt%YO-0.2wt%TiO2-1.0wt%-vacancy system.

Cell structure 50× 50 × 50 bcc unit lattices

Temperature [K] 1000

(31)

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

Temperature [K] 1000

5.74 9.08 unit of length : [nm] 45.38 22.69 Grain Boundary x y z

(a) Simulation cell

Grain Boundary [100] [100] [010] [010] misorientation angle x y 36.87

(b) Tilt grain boundary

(32)

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× 105MCstep 計算してもまだエネルギーがわずかに低下している．しかしながら，図 3.4(d),(e),(f) を見ると，5× 104MCstep以降は粒子が大きく移動しておらず，1.5× 105MCstepのときと 2.3× 105_MCstep_{のときの違いを図 3.5 に示すが，白線で囲ったようなごくわずかな違} いしかない．他のシミュレーションでも同様であり，計算時間との兼ね合いもあるた め本研究では 1.5× 105_MCstep_{後の粒子配置を平衡状態とみなし，以降で酸化物の分} 散形態および析出物内部の構造を検討する． -458.5 -458.0 -457.5 -457.0 -456.5 -456.0 -455.5 -455.0 -454.5 0 5 10 15 20

Potential energy [keV]

10 LMCstep 5 (a) LMC simulation. 10 MCstep 4 -505.0 -500.0 -495.0 -490.0 -485.0 -480.0 -475.0 -470.0 -465.0 -460.0 -455.0 -450.0 0 5 10 15 20 23

Potential energy [keV]

(b) Random MC simulation.

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

(33)

[010] [100] (a) At 106_LMCstep. [010] [100] (b) At 1.6× 106_LMCstep. [010] [100] (c) At 2× 106_{LMCstep(0MCstep).} [010] [100] (d) At 5× 104_MCstep. [010] [100] (e) At 1.5× 105_MCstep. [010] [100] (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] [100] (a) At 1.5× 105_MCstep. [010] [100] (b) At 2.3× 105_MCstep.

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

(34)

3.2.2 Fe-YO-TiO

2

の系における組成の影響

図 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] では越

(35)

とまったと考えられる．したがって，遠距離の引力相互作用を及ぼす「YO 粒子によるトラップ」と「高温によるエネルギー障壁の乗り越え」のバランスによって，TiO2が析出物を形成するか YO 粒子の周囲にトラップするかが分かれることが示唆される． [001] [010] [100] (a) Fe-0.2wt%YO-0.2wt%TiO2. [001] [010] [100] (b) Fe-0.2wt%YO-0.4wt%TiO2. [001] [010] [100] (c) Fe-0.4wt%YO-0.2wt%TiO2. [001] [010] [100] (d) Fe-0.4wt%YO-0.4wt%TiO2. [001] [010] [100] (e) Fe-0.6wt%YO-0.2wt%TiO2. [001] [010] [100] (f) Fe-0.6wt%YO-0.4wt%TiO2.

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

[001] [010] [100] (a) Fe-0.2wt%YO-0.2wt%TiO2. [001] [010] [100] (b) Fe-0.2wt%YO-0.4wt%TiO2. [001] [010] [100] (c) Fe-0.4wt%YO-0.2wt%TiO2. [001] [010] [100] (d) Fe-0.4wt%YO-0.4wt%TiO2. [001] [010] [100] (e) Fe-0.6wt%YO-0.2wt%TiO2. [001] [010] [100] (f) Fe-0.6wt%YO-0.4wt%TiO2.

(36)

TiO2 YO Fe [010] [100] (a) Fe-0.2wt%YO-0.2wt%TiO2. [010] [100] (b) Fe-0.2wt%YO-0.4wt%TiO2. [010] [100] (c) Fe-0.4wt%YO-0.2wt%TiO2. [010] [100] (d) Fe-0.4wt%YO-0.4wt%TiO2. [010] [100] (e) Fe-0.6wt%YO-0.2wt%TiO2. [010] [100] (f) Fe-0.6wt%YO-0.4wt%TiO2.

Fig.3.8 Magnified view around a precipitate of YO and TiO2in the simulation cell under 1000[K]. TiO2 YO _Fe [010] [100] (a) Fe-0.2wt%YO-0.2wt%TiO2. [010] [100] (b) Fe-0.2wt%YO-0.4wt%TiO2. [010] [100] (c) Fe-0.4wt%YO-0.2wt%TiO2. [010] [100] (d) Fe-0.4wt%YO-0.4wt%TiO2. [010] [100] (e) Fe-0.6wt%YO-0.2wt%TiO2. [010] [100] (f) Fe-0.6wt%YO-0.4wt%TiO2.

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

(37)

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 原子の結合エネルギー (ポテンシャルエネルギーの極小値) が他の組み合わせに比べて小さいためと考えられる．

(38)

YO

Fe TiO2 Al

[001] [010]

[100]

(a) Dispersion forms of YO, TiO2and Al.

[010] [100]

(b) Magnified view.

Fig.3.10 Final morphology of Fe-0.4wt%YO-0.2wt%TiO2-1.0wt%Al alloy un-der 1000[K]. YO Fe TiO2 Al [001] [010] [100]

(a) Dispersion forms of YO, TiO2and Al.

[010]

[100]

(b) Magnified view.

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

(39)

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]

[010]

[100]

(a) Dispersion morphology of YO and TiO2

particles.

[001]

[010]

[100]

(b) Magnified view.

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

(40)

0 20 40 60 80 100 2 4 6 8 10 1 3 5 7 9 Ratio [%]

the "n"th nearest sites from YO particles Fe

YO

TiO2 TiO2 TiO2

Vacancy

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%，

(41)

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

YO TiO2 Atom Vacancy

x

z y

(a) In Fe-0.2wt%YO-0.2wt%TiO2alloy.

x z y (b) In Fe-0.2wt%YO-0.4wt%TiO2 alloy. x z y (c) In Fe-0.6wt%YO-0.2wt%TiO2alloy. x z y (d) In Fe-0.6wt%YO-0.4wt%TiO2 alloy.

Fig.3.14 Dispersion morphology of YO and TiO2 particles.

z x y

(a) Bird view.

x

y

Grain Boundary

(b) x-y plane.

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

(42)

0 2 4 6 8 10 12 14 0 5 10 15 20 25 30 35 40 45 y coordinate [nm] Ratio [%] YO TiO2 Vacancy G.B. G.B. G.B.

Fig.3.16 The ratio of the distribution of YO, TiO2 and vacancy along the y coordinate (normal to the grain boundaries) in

Fe-0.2wt%YO-0.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粒子の粗大化が抑制さ

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

修士論文