平成 30 年度

文部科学省 国家課題対応型研究開発推進事業
英知を結集した原子力科学技術・人材育成推進事業

構造健全性評価の信頼性向上に向けた計算科学基盤の構築と破壊挙動の解明

成果報告書

平成 31 年 3 月
学校法人 東京理科大学

本報告書は、文部科学省の「英知を結集した原子力科学技術・人材育成推進事業」による委託業務として、学校法人

東京理科大学が実施した平成 28－30 年度「構造健全性評価の信頼性向上に向けた計算科学基盤の構築と破壊挙動の解明」の成果を取りまとめたもの
です。

目次
頁
概略 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ ⅲ
1. はじめに ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 1-1
2. 業務計画
2.1 全体計画

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 2.1-1

3. 業務の実施内容及び成果
3.1 転位の運動特性の分子動力学モデリング・・・・・・・・・・・・・・・・ 3.1-1

3.1.1 任意の混合転位の分子動力学モデリング・・・・・・・・・・・・・・ 3.1-1

3.1.2 任意の混合転位の運動特性の分子動力学モデリング（再委託先：電力中央研究所）

・・・・・・・・・・・・・・・・・・ 3.1-7

3.1.3 合金や不純物元素を含む金属中での転位挙動解析の準備（再委託先：電力中央研究所）

・・・・・・・・・・・・・・・・・・ 3.1-7

3.2 き裂先端での転位の放出・運動・弾性場を考慮した破壊靭性値の評価のための計算科学的基盤技術の構築・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3.2-1

3.2.1 転位の弾性場を考慮した破壊靭性値の評価法の確立・・・・・・・・ 3.2-1

3.2.2 照射や熱時効による破壊靭性値への影響評価（再委託先：電力中央研究所）

・・・・・・・・・・・・・・・・・・ 3.2-8

3.3 き裂先端における力学的状態と破壊靭性値の関係の評価

・・・・・・・・・ 3.3-1

3.3.1 大規模計算結果の可視化（再委託先：東洋大学）

・・・・・・・・・・ 3.3-1

3.3.2 超大規模解析に対して適している計算アルゴリズムの検討とコード整備（再委託先：東洋大学） ・・・・・・・・・・・・・・・・・・・・・・ 3.3-5
3.4 研究推進 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3.5-1
4. 結言
4.1 研究の到達点 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 4-1
4.2 今後の展望・次のステップに向けて・・・・・・・・・・・・・・・・・・ 4-2

i

表一覧
頁
表 3.1.2-1 Θi 値

・・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-15

表 3.1.3-1 代表的な Fe 単元系 EAM ポテンシャルのリスト

・・・・・・・・・・3.1-20

表 3.1.3-2 代表的な Fe 単元系 EAM ポテンシャルの基礎的な物性・・・・・・・3.1-20

表 3.1.3-3 Fe-Cu 系ポテンシャルの Cu-空孔のバインディングエネルギ
表 3.1.3-4 Fe-Cu 系ポテンシャルの空孔ジャンプエネルギ障壁・・・・3.1-24

・・・・・・・・3.1-25

表 3.1.3-5 Fe-Ni 系ポテンシャルの FeNi 規則合金の生成エネルギ

・・・・・・3.1-25

表 3.1.3-6 Fe-Ni 系ポテンシャルの BCC-Fe 中の欠陥生成エネルギ

・・・・・・3.1-26

表 3.1.3-7 Fe-Ni 系ポテンシャルの FCC-Ni 中の欠陥生成エネルギ

・・・・・・3.1-26

表 3.1.3-8 Fe-Cu-Ni 系ポテンシャルによる欠陥のバインディングエネルギ

・・3.1-27

図一覧頁
図 2.1-1 全体計画図 ・・・・・・・・・・・・・・・・・・・・・・・・・・・2.1-1
図 2.1-2 平成 28 年度月別全体計画図

・・・・・・・・・・・・・・・・・・ 2.1-1

図 2.1-3 平成 29 年度月別全体計画図

・・・・・・・・・・・・・・・・・・ 2.1-2

図 2.1-4 平成 30 年度月別全体計画図

・・・・・・・・・・・・・・・・・・ 2.1-2

図 3.1.1-1 {110}における原子配列と単位格子の抽出法

・・・・・・・・・・・3.1-1

図 3.1.1-2 混合転位の転位芯構造 ・・・・・・・・・・・・・・・・・・・・・3.1-3
図 3.1.1-3 分子動力学法により得られた刃状転位を含む結晶の応力―ひずみ線図
・・・・・・・・・・・・・・・・・・・3.1-4
図 3.1.1-4 {110}上の各転位のパイエルス応力の分子動力学法による計算結果

・3.1-5

図 3.1.1-5 {112}上の各転位のパイエルス応力の分子動力学法による計算結果

・3.1-6

図 3.1.1-6 {123}上の各転位のパイエルス応力の分子動力学法による計算結果

・3.1-6

図 3.1.1-7 らせん転位上のキンク対形成過程におけるエネルギ変化 ・・・・・・3.1-8
図 3.1.1-8 せん断応力を負荷した場合のらせん転位におけるキンク対の活性化エネルギ
・・・・・・・・・・・・・・・・・・・3.1-8
図 3.1.1-9 せん断応力を負荷した場合の混合転位におけるキンク付の活性化エネルギ
・・・・・・・・・・・・・・・・・・・3.1-8
図 3.1.2-1 モデル結晶の概略 ・・・・・・・・・・・・・・・・・・・・・・・3.1-11
図 3.1.2-2 変形中の転位芯位置と時間の関係 ・・・・・・・・・・・・・・・・3.1-11
図 3.1.2-3 負荷応力と転位速度の関係 ・・・・・・・・・・・・・・・・・・・3.1-11
図 3.1.2-4 易動度および摩擦応力の例 ・・・・・・・・・・・・・・・・・・・3.1-12
図 3.1.2-5 摩擦応力と転位線とバーガースベクトルのなす角度との関係（112 面、温度 100K） ・・・・・・・・・・・・・・・・・・・・・・3.1-12
図 3.1.2-6 摩擦応力と転位線とバーガースベクトルのなす角度との関係（112 面、温度 300K） ・・・・・・・・・・・・・・・・・・・・・・3.1-13ii

図 3.1.2-7 摩擦応力と転位線とバーガースベクトルのなす角度との関係（123 面、温度 100K） ・・・・・・・・・・・・・・・・・・・・・・3.1-13
図 3.1.2-8 摩擦応力と転位線とバーガースベクトルのなす角度との関係（123 面、温度 300K） ・・・・・・・・・・・・・・・・・・・・・・3.1-13
図 3.1.2-9 転位の易動度と転位線とバーガースベクトルのなす角度との関係（112 面、温度 100K） ・・・・・・・・・・・・・・・・・・・・・・3.1-14
図 3.1.2-10 転位の易動度と転位線とバーガースベクトルのなす角度との関係（112 面、温度 300K） ・・・・・・・・・・・・・・・・・・・・・・3.1-14
図 3.1.2-11 転位の易動度と転位線とバーガースベクトルのなす角度との関係（123 面、温度 100K） ・・・・・・・・・・・・・・・・・・・・・・3.1-14
図 3.1.2-12 転位の易動度と転位線とバーガースベクトルのなす角度との関係（123 面、温度 300K） ・・・・・・・・・・・・・・・・・・・・・・3.1-15
図 3.1.2-13 易動度のモデリング結果（112 面、100K）・・・・・・・・・・・・3.1-16

図 3.1.2-14 易動度のモデリング結果（112 面、300K）

・・・・・・・・・・・・3.1-16

図 3.1.2-15 易動度のモデリング結果（123 面、100K）・・・・・・・・・・・・3.1-16

図 3.1.2-16 易動度のモデリング結果（123 面、300K）

・・・・・・・・・・・・3.1-17

図 3.1.3-1 原子間ポテンシャルの分類 ・・・・・・・・・・・・・・・・・・・3.1-18
図 3.1.3-2 代表的な Fe 単元系の原子間ポテンシャル論文引用数 ・・・・・・・・3.1-21
図 3.1.3-3 らせん転位の第一原理計算と MendelvII ポテンシャルによる
Differential Displacement map

・・・・・・・・・・・・・・・・3.1-21

図 3.1.3-4 MendelevII、p10 ポテンシャルの熱膨張率・・・・・・・・・・・・3.1-21

図 3.1.3-5 らせん転位移動の最小エネルギ経路 ・・・・・・・・・・・・・・・3.1-22
図 3.1.3-6 Fe-Cu 系ポテンシャルの Cu の Fe への固溶限

・・・・・・・・・・・3.1-24

図 3.1.3-7 Fe-Ni 系ポテンシャルによる FeNi 合金の生成エネルギと規則合金の生成エネルギとの比較

・・・・・・・・・・・・・・・・3.1-27

図 3.1.3-8 Fe-Ni 系ポテンシャルによる Fe-Ni 相図・・・・・・・・・・・・・3.1-27

図 3.1.3-9 Fe-Cu-Ni 系ポテンシャルによる Cu-Ni 合金生成エネルギ
図 3.1.3-10 Fe-Cu-Ni 系ポテンシャルによる Cu-Ni 合金相図

・・・・・・・・・3.1-28

図 3.1.3-11 Bonny らの BCC 系 Fe 合金ポテンシャルの引用数の推移
図 3.1.3-12 Fe-Cu 合金の格子定数の温度依存性

・・・・・・3.1-28

・・・・・・・・・・・・・・3.1-31

図 3.1.3-13 Fe-Cu 合金の格子定数の Cu 濃度依存性
図 3.1.3-14 Fe-Cu 合金の弾性定数・・・・・3.1-28

・・・・・・・・・・・・・3.1-31

・・・・・・・・・・・・・・・・・・・・3.1-32

図 3.1.3-15 Fe-Cu 合金の{110}上の積層欠陥エネルギ

・・・・・・・・・・・・3.1-32

図 3.1.3-16 Fe-Cu 合金の{112}上の積層欠陥エネルギ

・・・・・・・・・・・・3.1-32

図 3.1.3-17 Fe-Cu 合金の{123}上の積層欠陥エネルギ

・・・・・・・・・・・・3.1-33

図 3.1.3-18 Fe-Cu 合金中の転位の摩擦応力と転位線とバーガースベクトルのなす角との関係（110 面、温度 100K）・・・・・・・・・・・・・・・・・・・・・3.1-33

図 3.1.3-19 Fe-Cu 合金中の転位の摩擦応力と転位線とバーガースベクトルのなす角との関係（110 面、温度 300K）

・・・・・・・・・・・・・・・・・・・・・3.1-33
iii

図 3.1.3-20 Fe-Cu 合金中の転位の易動度と転位線とバーガースベクトルのなす角との関係（110 面、温度 100K）・・・・・・・・・・・・・・・・・・・・・3.1-34

図 3.1.3-20 Fe-Cu 合金中の転位の易動度と転位線とバーガースベクトルのなす角との関係（110 面、温度 300K）

・・・・・・・・・・・・・・・・・・・・・3.1-34

図 3.2.1-1 貫通き裂問題の計算モデル

・・・・・・・・・・・・・・・・・・・3.2-2

図 3.2.1-2 貫通き裂の数値解析におけるバーガースベクトルと転位の数の関係

・3.2-2

図 3.2.1-3 異なる大きさのバーガースベクトルを用いた貫通き裂の応力拡大係数の計算結果

・・・・・・・・・・・・・・・・・・・・・・・・・・・3.2-3

図 3.2.1-4 バーガースベクトルを一定として計算した応力速度によるき裂解析結果
・・・・・・・・・・・・・・・・・・・3.2-4
図 3.2.1-5 バーガースベクトルの大きさを可変とした時の応力拡大係数の計算結果
・・・・・・・・・・・・・・・・・・・3.2-5
図 3.2.1-6 転位の数の制御数を変化させた応力拡大係数の計算結果
図 3.2.1-7 転位を用いたき裂解析の計算時間・・・・・・3.2-6

・・・・・・・・・・・・・・・・3.2-6

図 3.2.1-8 き裂の転位を用いた数値解析における不安定現象
図 3.2.1-9 応力レベルと計算の不安定現象の発生・・・・・・・・・3.2-7

・・・・・・・・・・・・・・3.2-8

図 3.2.1-10 き裂先端を含むすべり面上へのき裂先端からの転位の放出モデル
図 3.2.1-11 き裂先端から放出された転位同士の合体

・・3.2-9・・・・・・・・・・・・・3.2-10

図 3.2.1-12 転位の遮蔽効果の解析モデル

・・・・・・・・・・・・・・・・・・3.2-10

図 3.2.1-13 き裂先端から放出された転位・・・・・・・・・・・・・・・・・・3.2-11

図 3.2.1-14 転位によるき裂先端の遮蔽効果の計算結果
図 3.2.1-15 転位の易動度と破壊靭性値の関係

・・・・・・・・・・・・3.2-12

・・・・・・・・・・・・・・・・3.2-12

図 3.2.1-16 き裂先端から任意のすべり面への転位の放出モデル

・・・・・・・・3.2-13

図 3.2.1-17 き裂面から 45°傾いたすべり面への転位の放出を伴う数値解析の結果・・・・・・・・・・・・・・・・・・・3.2-14
図 3.2.1-18 き裂先端から放出された転位の挙動・・・・・・・・・・・・・・・3.2-15

図 3.2.1-19 温度条件を変化させた転位によるき裂先端の遮蔽効果の計算結果
図 3.2.1-20 破壊靭性値の計算結果と Tanaka らによる実験結果
図 3.2.1-21 斜めのすべり面上の転位による遮蔽の効果
図 3.2.2-1 Fe-Cu 合金中の格子定数の温度依存性

・・3.2-16・・・・・・・・3.2-16

・・・・・・・・・・・・3.2-17

・・・・・・・・・・・・・・3.2-19

図 3.2.2-2 Fe-Cu 合金中の格子定数の Cu 濃度依存性・・・・・・・・・・・・3.2-20

図 3.2.2-3 Fe-Cu 合金中の弾性定数の Cu 濃度依存性・・・・・・・・・・・・3.2-20

図 3.2.2-4 Fe-Cu 合金中の一般化積層欠陥エネルギ（110 面）

・・・・・・・・3.2-20

図 3.2.2-5 Fe-Cu 合金中の一般化積層欠陥エネルギ（112 面）

・・・・・・・・3.2-21

図 3.2.2-6 Fe-Cu 合金中の一般化積層欠陥エネルギ（123 面）

・・・・・・・・3.2-21

図 3.2.2-7 転位動力学解析モデルの概要

・・・・・・・・・・・・・・・・・・3.2-22

図 3.2.2-8 き裂・転位相互作用の転位動力学解析結果の例

・・・・・・・・・・3.2-23

図 3.2.2-9 負荷応力から計算される応力拡大係数とき裂先端の応力拡大係数の比較（クラスター間隔の影響）・・・・・・・・・・・・・・・・・・・3.2-23iv

図 3.2.2-10 負荷応力から計算される応力拡大係数とき裂先端の応力拡大係数の比較（クラスターサイズの影響）・・・・・・・・・・・・・・・・・・3.2-24

図 3.3.1-1 応力分布：xx 応力成分の等値面表示
図 3.3.1-2 応力分布：相当応力の断面コンター表示
図 3.3.1-3 転位線の可視化

・・・・・・・・・・・・・・3.3-2
・・・・・・・・・・・・3.3-2

・・・・・・・・・・・・・・・・・・・・・・・3.3-3

図 3.3.1-4 応力分布と対応する転位線の同時表示

・・・・・・・・・・・・・3.3-3

図 3.3.2-1 き裂入り試験片のマクロ解析とき裂先端でのミクロ解析
図 3.3.2-2 DD-FEM 連成解析のデータフローの様子
図 3.3.2-3 連成解析によるき裂の様子（体心立方）

CNA : Common Neighbor Analysis （共通近接分析）
DD: Dislocation Dynamics

・・・・・・・・・・・・・3.3-5

・・・・・・・・・・・・・・・・・・3.3-6

略語一覧
BCC : Body Centered Cubic

・・・・・3.3-4(転位動力学)EAM : Embedded Atom Method(埋め込み原子法)

FEM : Finite Element Method(有限要素法)NEB : Nudged Elastic Band

（弾性バンド）v

概略
供用期間中の炉心部における中性子の照射や熱時効による影響によって低合金鋼（原子炉圧力容器鋼）、ステンレス鋳鋼やオーステナイトステンレス鋼は脆化することが知られている。すな
わちこれらの金属材料は、延性脆性遷移温度の上昇や、上部棚エネルギが低下するため、機器の健全性評価を行う上で考慮すべき重要な因子となっている。今後、国のエネルギ政策に基づき原子力発電による電力供給割合を 20～22％確保していくためには、安全性が確保された軽水炉を、40 年を超えて利用してゆくことが必要となるが、40 年を越えた軽水炉の継続運転の判断の時点のみならず、その後も継続して金属材料の劣化に関する理解の改善に努め、最新の知見に基づいて健全性を保証するとともに、広く社会全般への説明性を向上させることが必要である。このためには、原子炉構造材料の脆化や延性—脆性遷移のメカニズムを詳細に把握し、脆化の程度を精度良く予測・評価することを可能にする新たな基盤技術の開発への不断の取組みが不可欠である。
金属材料の脆化や延性—脆性遷移のメカニズムに関する研究として、①実験的研究を用いた材料の組織観察及び機械的特性の測定や、②計算科学を応用したマルチスケール材料モデリングによるアプローチがこれまでにも行われてきた。現状では、①の実験的研究が大きく先行しており、ミクロ組織の変化や機械的特性の変化に関する知見は非常に充実してきている。しかし、実験的研究に基づくミクロ組織と機械的特性を繋ぐメカニズムの理解は必ずしも十分ではなく、経験的なアプローチが用いられることが多い。この問題の解決には現在の実験技術は十分でなく、②のマルチスケール材料モデリングの活用が不可欠である。マルチスケール材料モデリングによるアプローチでは、個々の長さ・時間スケールにおける現象をそのスケールに適した計算科学的手法を用いてモデル化し、そのモデルを次のスケールで用いることによってボトムアップ式にミクロな（原子論的な）特性からマクロな機械的特性を導くことが検討されてきている。すなわち、マルチスケール材料モデリングにより金属材料の脆化のメカニズムを知ることができれば、ミクロからマクロにわたる脆化のメカニズムの全貌を知ることが可能となる。しかし、これまでのマル
チスケール材料モデリングによる研究では、マルチスケール材料モデリングの全体像を描く一方で、個々のスケールでの現象の計算機シミュレーションや実験的測定に留まり、スケール間の相互作用やミクロからマクロまでを横断した脆化メカニズムの理解には至っていない。
本研究では、原子炉構造材料の延性—脆性遷移に焦点を絞り、延性—脆性遷移の基本的なメカニズムを解明するためのマルチスケール材料モデリングを基礎とする計算科学的基盤技術の構築を行う。さらに、本研究で構築した計算科学的基盤技術によって得られる延性—脆性遷移のメカニ
ズムに関連する情報と、これまでに蓄積されてきた実験的観察結果を統合することにより、工学的に必要なパラメータ（例えば破壊靭性値）を導き出すことを目指す。具体的には、延性—脆性
遷移のメカニズム解明のための計算科学的アプローチを描き、それを実行可能とする分子動力学、転位動力学及び有限要素法のコードを整備する。構築した計算科学的基盤技術を用いて、フェラ
イト鋼に見られる延性—脆性遷移のメカニズムを解明する。さらに、延性—脆性遷移とき裂先端の力学場との関係を明らかにする。これにより、供用中の材料の脆化の程度を、計算科学的アプローチを用いて予測・評価することを可能にする。
本研究の成果として以下のことを得ることができた。混合転位の運動の性質の温度依存性の解明とモデル化においては、平成 29 年度までに、パイエルス応力と転位芯構造の関係を明らかにし、特にいくつかの混合転位においてパイエルス応力が高くなることが明らかとなった。また、vi

転位芯に関する分析から、転位が一度に移動するための応力、またはキンクが転位に沿って移動するための応力であるかについて、パイエルス応力の意味が異なることが明らかになった。NEB法を用いたキンク対の形成に必要な活性化エネルギを計算し、活性化エネルギのせん断応力依存性を調べた。さらに転位の運動特性について転位の易動度を詳細に調べ、その結果を、サイン関数を用いた式でフィッティングすることによって、様々な転位の易動度が表現可能であることを示した。鉄中に銅原子が固溶した場合の、弾性定数、格子定数、積層欠陥エネルギを、銅濃度を変化せて調べた。さらに、転位の易動度と摩擦応力について調べた結果、銅の含有量は、転位の易動度に対する影響は小さく、摩擦応力に対して大きく影響することがわかった。
転位運動モデルを用いた、き裂先端での転位の発生・運動による遮蔽の効果を転位動力学に基づいて考慮することによる破壊靭性値の温度依存性の計算技術の開発においては、平成 29 年度
までに、転位を用いたき裂の数値解析法を貫通き裂に拡張し、応力速度を境界条件とする数値解析において高精度に応力拡大係数を算出する方法を開発することに成功した。き裂の転位動力学
コードにき裂先端からの転位の放出モデルを実装することで、き裂先端からの転位の放出およびき裂先端の遮蔽効果の計算を可能にした。平成 30 年度では、分子動力学法による転位の運動特性のモデリングの結果として得られた転位の易動度の温度および方位依存性モデルを実装し、転位の遮蔽に対する温度の影響をシミュレートすることを可能にした。このモデルを用いて貫通き裂における転位の遮蔽効果のシミュレーションを実施した結果、脆性―延性遷移温度が 140K から 160K の間であることが計算された。実験結果では脆性―延性遷移温度はおよそ 150K となっており、本計算結果は実験結果と定量的に一致した。すなわち、本研究によって開発された計算科
学基盤を用いて計算することによって、実験結果を再現する結果を得ることに成功した。三次元の転位の運動に基づく計算機シミュレーションによってこのような結果を得たのは、世界初の成果である。さらに転位のすべり面やすべり方向によっては、逆遮蔽を起こす場合があることがわかった。
照射・熱によって金属中に形成される照射損傷を転位動力学計算で考慮することで、照射・熱による破壊靭性値の変化を計算する技術の開発においては、平成 29 年度までに鉄と銅の積層欠
陥エネルギを計算し、銅濃度が高くなるにつれ、積層欠陥エネルギが低下することがわかった。
平成 30 年度では、これらのパラメータを用いた転位と銅析出物の相互作用モデルを実装し、転位の遮蔽効果への銅析出物の影響について計算を実施し結果、析出物が大きくなるにつれ、さらに析出物間間隔が狭くなるにつれ転位による遮蔽の効果は弱まり破壊靭性値が小さくなる傾向（脆化する傾向）を得た。これは、析出物の分布や大きさが与える破壊靭性値への影響を計算機
シミュレーションを用いて捉えた世界初の結果である。
き裂先端での転位動力学が有限要素法による変形解析を融合することで種々の力学的状態における破壊靭性値の変化を計算し、き裂先端での拘束条件の影響を評価するために、まず、計算結果の可視化環境の構築を行なった。平成 29 年度までに、転位が作る応力の可視化のためのフレームワークを構築し、転位が作る複雑な応力場の可視化を可能にした。平成 30 年度では、転位が作る応力分布の可視化、および転位線との重ね合わせが行えるようなシステムの試作および実装を行い、転位線セグメント群の表示と応力の断面表示および等値面表示機能などの各種機能を実装した。一方、き裂先端での拘束条件の影響の評価するために、平成 29 年度までに、転位動力学法と有限要素法を重ね合わせの原理により融合させることに成功した。平成 30 年度では、
vii

これまで行ってきた単純引っ張り応力の境界条件に加え、せん断引っ張り、曲げ引っ張り応力な
どの異なる応力条件によるシミュレーションにも成功し、破壊靭性値との関係を整理することができた。
本研究によって達成された破壊靭性値に関するマルチスケール材料モデリングは、原子力分野に限らず計算科学、計算力学の分野においても大きな意義があると言える。これまでの計算機シミュレーションを用いた原子力材料の脆化に関する研究は、照射による微視的組織の形成や、微視的組織と転位の相互作用、材料の降伏応力の上昇に関するものに留まり、延性―脆性遷移の本
質的な（直接的な）物性値である破壊靭性値を調べることが可能になったことは、今後の材料研究において多大なインパクトを与えるものである。
本研究課題では、マルチスケール材料モデリングの考え方を用いて、BCC 鉄の転位運動のモデ
ルを転位動力学法に実装することに成功している。この転位動力学法を用いることで、これまで
計算機シミュレーションを用いて数値的に解析することが困難であった BCC 鉄の塑性変形挙動および降伏応力のシミュレーションが可能になる。このような転位動力学法による降伏応力の数値解析と、本研究課題で実現した破壊靭性値の数値解析を並行して行うことにより、降伏応力の変化と破壊靭性値の変化の関係や、これらに寄与する微視的組織の実態について研究することが可能になる。さらに、本研究で開発した転位動力学法をこれに応用すれば、BCC 鉄の硬さの変化についての信頼性の高い研究を展開することが期待できる。以上のように、本研究課題で開発した転位動力学法に基礎をおく計算科学基盤は、破壊靭性値の数値解析法としての役割のみならず、
降伏応力の変化や硬さの変化に関する数値解析的アプローチの信頼性向上にも資するものである。
これらの異なる物性値の研究を統一的に実施することで、巨視的な物性値と微視的組織の関係のより詳細な関係を導き出し、さらに新しい試験法に向けた知見を与える発展性を秘めている。

viii

1． はじめに供用期間中の炉心部における中性子の照射や熱時効による影響によって低合金鋼（原子炉圧力容器鋼）、ステンレス鋳鋼やオーステナイトステンレス鋼は脆化することが知られている。すな
わちこれらの金属材料は、延性脆性遷移温度の上昇や、上部棚エネルギが低下するため、機器の健全性評価を行う上で考慮すべき重要な因子となっている。今後、国のエネルギ政策に基づき原子力発電による電力供給割合を 20～22％確保していくためには、安全性が確保された軽水炉を、40 年を超えて利用してゆくことが必要となるが、40 年を越えた軽水炉の継続運転の判断の時点のみならず、その後も継続して金属材料の劣化に関する理解の改善に努め、最新の知見に基づいて健全性を保証するとともに、広く社会全般への説明性を向上させることが必要である。このためには、原子炉構造材料の脆化や延性—脆性遷移のメカニズムを詳細に把握し、脆化の程度を精度良く予測・評価することを可能にする新たな基盤技術の開発への不断の取組みが不可欠である。
金属材料の脆化や延性—脆性遷移のメカニズムに関する研究として、①実験的研究を用いた材料の組織観察及び機械的特性の測定や、②計算科学を応用したマルチスケール材料モデリングによるアプローチがこれまでにも行われてきた。現状では、①の実験的研究が大きく先行しており、ミクロ組織の変化や機械的特性の変化に関する知見は非常に充実してきている。しかし、実験的研究に基づくミクロ組織と機械的特性を繋ぐメカニズムの理解は必ずしも十分ではなく、経験的なアプローチが用いられることが多い。この問題の解決には現在の実験技術は十分でなく、②のマルチスケール材料モデリングの活用が不可欠である。マルチスケール材料モデリングによるアプローチでは、個々の長さ・時間スケールにおける現象をそのスケールに適した計算科学的手法を用いてモデル化し、そのモデルを次のスケールで用いることによってボトムアップ式にミクロな（原子論的な）特性からマクロな機械的特性を導くことが検討されてきている。すなわち、マルチスケール材料モデリングにより金属材料の脆化のメカニズムを知ることができれば、ミクロからマクロにわたる脆化のメカニズムの全貌を知ることが可能となる。しかし、これまでのマル
チスケール材料モデリングによる研究では、マルチスケール材料モデリングの全体像を描く一方で、個々のスケールでの現象の計算機シミュレーションや実験的測定に留まり、スケール間の相互作用やミクロからマクロまでを横断した脆化メカニズムの理解には至っていない。
本研究では、原子炉構造材料の延性—脆性遷移に焦点を絞り、延性—脆性遷移の基本的なメカニズムを解明するためのマルチスケール材料モデリングを基礎とする計算科学的基盤技術の構築を行う。さらに、本研究で構築した計算科学的基盤技術によって得られる延性—脆性遷移のメカニ
ズムに関連する情報と、これまでに蓄積されてきた実験的観察結果を統合することにより、工学的に必要なパラメータ（例えば破壊靭性値）を導き出すことを目指す。本研究で構築する延性—
脆性遷移のメカニズムの解明を目的とした計算科学的基盤技術には以下の項目があり、それぞれの項目についてその前後の項目との相互作用を考える。
・ 混合転位の分子動力学解析によって、刃状かららせん転位にわたる全ての混合転位の運動の性質の温度依存性を解明し、モデル化する。
・ 転位運動モデルを用い、き裂先端での転位の発生・運動による遮蔽の効果を転位動力学に基づいて考慮することで、破壊靭性値の温度依存性計算する技術を開発する。
・ 照射・熱によって金属中に形成される照射損傷を転位動力学計算で考慮することで、照射・熱による破壊靱性値の変化（脆化）を計算する技術を開発する。
1−1

・ き裂先端での転位動力学や有限要素法による変形解析を融合することで種々の力学的状態における破壊靱性値の変化を計算し、き裂先端の拘束条件の影響を定量的に評価する。
金属材料の延性—脆性遷移のメカニズムを調べるための、計算科学的基盤技術を構築する。
具体的には、延性—脆性遷移のメカニズム解明のための計算科学的アプローチを描き、それを実行可能とする分子動力学、転位動力学及び有限要素法のコードを整備する。構築した計算科学的基盤技術を用いて、フェライト鋼に見られる延性—脆性遷移のメカニズムを解明する。さらに、
延性—脆性遷移とき裂先端の力学場との関係を明らかにする。これにより、供用中の材料の脆化の程度を計算科学的アプローチを用いて予測・評価することを可能にする。

1−2

2．業務計画
2.1 全体計画
本業務の全体計画図を図 2.1-1 に示す。平成 28 年度から平成 30 年度までの月別全体計画図を
図 2.1-2（平成 28 年度）、図 2.1-3（平成 29 年度）、図 2.1-4（平成 30 年度）にそれぞれ示す。

図 2.1-1 全体計画図

図 2.1-2 平成 28 年度月別全体計画図

2.1-1

図 2.1-3 平成 29 年度月別全体計画図

図 2.1-4 平成 30 年度月別全体計画図

2.1-2

3．業務の実施内容及び成果
3.1 転位の運動特性の分子動力学モデリング （H28～H29）
3.1.1 任意の混合転位の分子動力学モデリング(1) 任意の混合転位を含む原子シミュレーションモデルの作成方法
鉄を対象として、体心立方格子中のすべり系に注目し、任意の混合転位を含む原子系を作成する方向を確立した。具体的には、体心立方格子（Body Centered Cubic：BCC）のす
べり面である{110}、{112}、{123}の原子配列を分析することによって、任意の方向について周期性が最小となる単位格子を見つける方法を考案し、任意の混合転位の原子モデルを作成する方法を確立にした。単位格子を x、y、z 方向にそれぞれ 3 次元的に配列することによって、混合転位を含む任意の大きさの原子モデルを作成することが可能である。その一例として{110}における原子配列と、任意の方向についての単位格子の作成の考え方について図 3.1.1-1 に示す。

図 3.1.1-1 {110}における原子配列と単位格子の抽出法
規則的に原子が配列している結晶構造中で単位格子が周期性を満足する十分条件は「格子が直方体である」「格子の頂点に原子がある」の 2 点である。すべり面法線方向（Y 軸）を中心に単位格子を回転させるので XZ 平面に注目すると、この 2 点を満足する単位格子では 図 3.1.1-1 に示すような三角形 A と B は相似の関係にある。この相似より任意の角度φにあるベクトル x を決定することでベクトル z を導くことができ、混合転位の単位格子(黄枠)を作成することができる。その時の z = (z1, z2) は x = (x1, x2) と A と B の相似比 s を用いて式 (3.1.1-1) のようになる。
z1 =

sx2

z2 = sx1(3.1.1-1)p三角形の相似比 s は XZ 平面の最小の単位格子の縦横比より s = 2 2と求まる。しかし、図3.1.1-1 の Z’,Z”が頂点となることも考えるとこの単位格子は最小の単位格子ではないので、ユークリッドの互除法に剰余算を使って最小の単位格子を求めた。以上によって、特定のすべり面を XZ 平面として、その他の 2 方向を任意の結晶方位とする原子モデルの作成が可能になった。同様に、三角形の相似比について、{112}、{123}について求めた結果、
ppそれぞれ s = 2 6、 s = 2 14となった。単位格子を用いて作成した原子モデルの妥当性を確認するために、混合転位を作成していない原子モデルのポテンシャルエネルギを計算し、最も簡単な X、Y、Z 方向を[100]、[010]、[001]とした完全結晶の原子モデルを用いて計算した結果と比較を行った。計算に用いたポテンシャルは鉄に対して提案された埋め込み原子法（Embedded Atom Method:EAM）に基づく原子間ポテンシャルを用いた(2)。具体的には、
3.1−1

原子モデルを作成し、全方向に周期境界条件を与えた。その状態におけるポテンシャルエ
ネルギを、原子間ポテンシャルを用いて計算し、原子系全体のエネルギから原子数を除することによって、凝集エネルギを求めた。その結果、簡単な原子系を用いた場合の凝集エ
ネルギは-4.013eV であった。一方、本研究で開発した単位格子作成法を用いて作成した原子系を用いた計算結果は-4.013eV であり、簡単な原子系を用いた計算結果と一致する数値を得ることができた。これにより、単位格子が正しく作成されていることを確認することができた。次に、作成した任意の結晶包囲を有する原子モデルに転位を導入する。混合転位は変位に刃状成分とらせん成分を含んだ転位である。すなわち、混合転位の作成は、刃状転位の作成方法とらせん転位の作成方法の重ね合わせで実現することができる。刃状成分の導入は、すべり面より上方にある原子モデルに対し、以下の変位を与えることで行う。
dLEdge = b sin ✓（3.1.1-2）

一方、らせん成分についても同様に、すべり面より上方にある原子モデルに対し、以下の変位を与えることで導入する。
dLScrew = b cos ✓（3.1.1-3）以上の変形を原子系に与えることによって、Z 方向に平行な直線形状の混合転位を原子モデル中に作成することができる。また、このモデルに対し、Z 方向に周期境界条件を与えることで、無限に長い転位を仮定することができ、さらに X 方向にも周期境界条件を与えることで、無限に長い転位が周期的に配列したモデルを仮定することができる。Y 方向の上下端面付近の原子に対し、Y 方向の変位を固定し、混合転位のバーガースベクトル方
向に変位を与えることでせん断応力を与え、転位の運動を調べる分子動力学シミュレーシ
ョンを行うことができる。このような原子系のモデルは Periodic Array of Dislocations(PAD)モデルと呼ばれている(1)。
(2) 混合転位の転位芯構造の分析
任意の混合転位の PAD モデルを用いて、混合転位の芯（転位芯）構造の分析を行った。
計算に用いた原子系は、Y 方向の大きさが、104.99Å、X 方向と Z 方向はそれぞれ 100、20Å以上の大きさとした。X、Z 方向の大きさは、混合転位の種類（結晶方位）によって、単位格子の大きさが異なるため、一定の値を用いて様々な混合転位のシミュレーションを実施することはできない。転位芯の構造解析においては、釣り合い状態にある原子系における転位芯構造を調べる必要がある。そこで本研究では、原子間ポテンシャルを用いた力の計算と、数値急冷法と共役勾配法を用いた構造緩和計算を組み合わせた分子動力学法を用いて、原子系の釣り合い状態を求めた。その結果を図 3.1.1-2 に示す。

3.1−2

(a)θ＝0°(b)θ＝54.7356°（d）θ＝5.5706°(e)θ＝8.0495°(c)θ＝109.4712°(f)θ＝14.4201°

図 3.1.1-2 混合転位の転位芯構造
図に示す球は、Common Neighbor Analysis（CNA）(3)を用いた局所原子構造解析の結果として抽出した、転位芯付近の原子である。太線は、Dislocation Extraction Algorithm（DXA）(4)

によって決定した転位線である。図 3.1.1-2(a)〜(c)を見ると、DXA によって求めた転位線は直線となっており、また転位芯付近にある原子の配列が、転位線を対称線として、左右対称の構造を取っていることが分かる。一方、図 3.1.1-2（d）〜（f）を見ると、転位
線に段差があることがわかった。この段差はキンクと呼ばれるものである。以上から、直線的かつ左右対称の転位芯構造を持った転位と、潜在的にキンクを持つ転位があることがわかった。
(3) 混合転位のパイエルス応力の計算
次に、混合転位のパイエルス応力の計算を行った。用いた原子系モデルは、転位芯構造を求める計算に用いたものと同様である。境界条件として、Y 方向上下 10.6Åの領域に含まれる原子に対して、バーガースベクトルと同じ方向に強制変位を与えた。なお、上下の領域の原子に対して、それぞれ反対向きに変位を与え、せん断変形を加えた。パイエルス応力とは、熱による原子の振動の影響などを含めず、原子配列（格子構造）由来の転位の移動に対する障壁であるため、本計算では、8.24×10-4Åを強制変位として与え、数値急冷法と共役勾配法を用いて構造緩和を行い、各変形レベルにおいて力学的な釣り合い状態を追従するシミュレーションを行った。原子系に与えたせん断応力は、強制変位を与えた原子に対して、その他の内部に位置する原子が与える力を計算し、その総和をすべり面の面積で除することによって計算した。せん断ひずみは、与えた強制変位を原子モデルの高さで除することで求めた。図 3.1.1-3 に刃状転位を含む原子モデルを用いた時の応力—ひず
3.1−3み関係の出力結果を示す。初期の段階において、応力とひずみの間に線形の関係が見られる。これは原子モデルが弾性的に（塑性変形を含まずに）変形している挙動を示している。
その後、せん断応力がおよそ 80MPa に達すると、急激に応力が低下している挙動が見られる。これは、せん断応力がパイエルス応力に達したことを意味しており、転位が現在の安定位置から隣の安定位置まで移動した結果である。隣の安定位置に移動すると、転位の移動による塑性変形によって、原子モデルの弾性変形が緩和され、その結果応力が低下した
ため、転位は継続して移動することができず、再びせん断応力がパイエルス応力に達するまで、そこの位置に止まっている。同様の計算を刃状転位以外の混合転位に対して行い、
パイエルス応力についてまとめた。その結果を図 3.1.1-4 に示す。図を見ると、多くの混合転位のパイエルス応力は 50MPa 前後であるのに対し、らせん転位やいくつかの混合転位について、非常に大きなパイエルス応力が計算された。らせん転位のパイエルス応力が高
いことはこれまでの計算により明らかにされてきているが、その他の混合転位においても高いパイエルス応力を示すことがあることは新たな発見である。このようなパイエルス応
力の違いの原因（メカニズム）を理解するために、図 3.1.1-2 に示した転位芯構造とパイ
エルス応力の大きさの関係について調べた。図 3.1.1-2(a)〜(c)に示した直線的な転位芯構造を持つ転位について、計算されたパイエルス応力の大きさを調べると、これらの転位が非常に高いパイエルス応力を持つ転位であることが分かる。一方、図 3.1.1-2(d)〜(f)に示した転位芯構造をもつ転位のパイエルス応力を見ると、低いパイエルス応力であることがわかった。また、パイエルス応力の計算過程における転位の移動の様子を確認すると、直線的な転位芯構造を持つ転位は、全ての転位が一度に隣の安定位置に移動（ジャンプ）
していることがわかった。一方、キンクを含む転位芯構造を持つ転位の場合は、キンクが転位線に沿って移動することで、転位が徐々に前方に進む様子を確認することができた。
すなわち、直線的な転位芯構造を持つ転位のパイエルス応力は、転位が一度に隣の安定位
置にジャンプするために必要なせん断応力であり、一方キンクを持つ転位のパイエルス応
力は、キンクが転位線に沿って移動するために必要なせん断応力であると理解できる。以上のように、パイエルス応力の大きさは、転位芯構造に強く依存し、パイエルス応力の意
味も転位芯構造によって異なってくることがわかった。

図 3.1.1-3 分子動力学法により得られた刃状転位を含む結晶の応力—ひずみ線図
3.1−4

図 3.1.1-4 {110}上の各転位のパイエルス応力の分子動力学法による計算結果
同様に、{112}上の混合転位についてもパイエルス応力の計算を行った。原子モデルや境界条件は、{110}上の混合転位のパイエルス応力の計算のために用いたものと同様である。転位芯構造について調べた結果、{110}上の混合転位で見られたような直線的な転位芯構造を持つ転位と、潜在的にキンクを持つ転位の２種類に分類できることがわかった。
次に、Y 方向上下 10.6Åの領域に含まれる原子に対して、バーガースベクトルと同じ方向に強制変位を与え、数値急冷法と共役勾配法を用いることによって、パイエルス応力の計算を行なった。{112}をすべり面とするせん断変形においては、同一のせん断方向であっても、せん断する向きによって変形の性質が異なることが知られている。これは双晶、反双晶変形と呼ばれ、らせん転位の反双晶変形において双晶変形に比べてせん断変形の抵抗が大きいことが知られている。本研究では、混合転位のパイエルス応力においても、らせん転位に対して知られている双晶と反双晶変形の影響が見られるかについても注目した。
すなわち{112}上の混合転位のパイエルス応力の計算においては、強制変位を同一の方向に対して、その向きを変えた計算も行った。図 3.1.1-5 に計算によって求めたパイエルス応力を示す。双晶変形を与えた場合（図 3.1.1-5(a)）において、{110}と同様に、いくつかの角度を有する混合転位のパイエルス応力について、その他の混合転位のパイエルス応
力よりも大きい値を示していることが分かる。転位芯構造を確認すると、このような高い
パイエルス応力を有する混合転位においては、直線的な形状の転位芯構造を有していた。
一方、低いパイエルス応力を有する混合転位においては、その転位芯構造において、潜在的にキンクが存在し、{110}上の混合転位の転位芯構造とパイエルス応力の関係と定性的に一致した傾向であることがわかった。一方、双晶（図 3.1.1-5(a)）と反双晶変形（図3.1.1-5(b)）におけるパイエルス応力を比較すると、らせん転位のパイエルス応力について、およそ 2 倍程度双晶変形の場合よりも反双晶変形の場合の方が大きいことが分かる。
しかしその一方で、混合転位についてはパイエルス応力の大きな違いは見られない。以上
3.1−5

の結果から、{112}上の混合転位においては、{110}上の混合転位と同様の転位芯構造とパイエルス応力の関係が見られ、双晶と反双晶変形については、らせん転位以外の混合転位
においては、パイエルス応力に大きな違いが見られないことがわかった。

(a)双晶変形(b)反双晶変形

図 3.1.1-5 {112}上の各転位のパイエルス応力の分子動力学法による計算結果
図 3.1.1-6 に{123}上の混合転位のパイエルス応力の計算結果を示す。双晶変形と反双晶
変形ともに角度が 0°であるらせん転位の場合に極端に大きなパイエルス応力となっていることがわかる。これは、鉄中で多くのらせん転位が実験観察されている結果と整合性のある結果である。さらに、角度が 60°および 120°の時にピークがみられ、それ以外のい
くつかの角度においてもパイエルス応力のピークが見られる。この角度は、転位線の方向が体心立方格子の近接方向と一致する。このことは、{110}と{112}の転位においても見られた同様の傾向である。また、双晶変形と反双晶変形を比較すると、定量的にはやや違いが見られるものの、全体としての定性的な傾向に大きな違いは見られなかった。

(a)双晶変形(b)反双晶変形

図 3.1.1-6 {123}上の各転位のパイエルス応力の分子動力学法による計算結果(4) 転位の熱的運動におけるキンク対形成
高いパイエルス応力を持つ転位は、転位線上の一部が張り出し、キンクの対が形成され、そのキンクが転位線に沿って移動することで、転位全体の移動が起こることが知られている。このようなメカニズムのことをキンク対メカニズムと呼び、熱活性化過程で発生する
3.1−6

現象であることが知られている。したがって、パイエルス応力の高い転位の運動を知るためには、キンク対の発生を知る必要があり、そのためにはキンク対の活性化エネルギを知る必要がある。本研究では、Nudged Elastic Band (NEB)法(5)と呼ばれる最小エネルギ経路を探索する方法を用いて、キンク対形成における最小エネルギ経路を求めることにより、
キンク対形成に必要な活性化エネルギを算出する。計算モデルは、パイエルス応力の計算に用いたモデルと同様である。NEB 法では、初期状態を始点、最終状態を終点とし、その間を複数のイメージによって表現する。始点、終点、イメージは隣接するイメージとバネで接続されており、ポテンシャルの微分から計算されるエネルギを最小化する方向である
力の成分から、イメージ間方向を取り除いた力とイメージ間を繋ぐバネから受ける力によって、全てのイメージの原子の位置を更新する。これにより、始点から終点までを繋いだイメージの計算が可能であり、各イメージのエネルギをプロットすると、始点から終点に至るために必要なエネルギを知ることが可能である。らせん転位についてキンク対の形成過程におけるエネルギ変化の計算結果を図 3.1.1-7 に示す。キンク対の形成段階において、エネルギの上昇が見られる。その後、キンク対が転位に沿って移動し対消滅することによってエネルギが低下している。また、らせん転位においてキンク対の活性化エネルギはお
よそ 0.7eV であることがわかった。次に、せん断応力を負荷した場合のキンク対の活性化
エネルギについて計算を行った。本計算では、NEB 法による計算時に与えるせん断応力に相当する力を、Y 方向に変位を固定している原子に対して与えることによって行った。計
算結果を図 3.1.1-8 に示す。図中には、Aono らによって行われた純鉄による実験結果(6)も示す。計算結果と実験結果を比較すると、定量的に良い一致が見られる。すなわち、NEB法を用いたキンク対の活性化エネルギの計算結果の信頼性および精度を確認することができる。さらに、せん断応力が大きくなるにつれて、キンク対を形成するために必要な活性化エネルギが低くなることがわかる。図 3.1.1-9 にらせん転位およびパイエルス応力の計算において高いパイエルス応力であった混合転位におけるキンク対を形成するために必要な活性化エネルギの計算結果を示す。この計算においても与えるせん断応力を変化させ、その依存性を確認した。いずれの混合転位においても、せん断応力が大きくなるにつれて、キンク対の形成に必要な活性化エネルギが小さくなる傾向を見ることができる。さらに、
キンク対の形成に必要な活性化エネルギは、らせん転位において最大で、その他の混合転
位では、転位線の方向の近接番号が大きくなるにつれて（原子間の近接距離が大きくなる
につれて）、活性化エネルギが低くなることがわかった。

3.1−7

図 3.1.1-7 らせん転位上のキンク対形成過程におけるエネルギ変化

図 3.1.1-8 せん断応力を負荷した場合のらせん転位におけるキンク対の活性化エネルギ

図 3.1.1-9 せん断応力を負荷した場合の混合転位におけるキンク対の活性化エネルギ

3.1−8

(5) まとめ鉄を対象として、体心立方格子中のすべり系に注目し、任意の混合転位を含む原子系を作成する方法を確立した。原子モデルの作成は、結晶の周期性の単位を構成する単位格子を周期的に配置することによって行うため、任意の結晶方位に対する単位格子を決定する方法を考案した。作成した原子モデルを用いて、混合転位の転位芯構造を計算し、作成した原子モデルの妥当性を評価するとともに混合転位の転位芯構造に関する分析を行った。
また、転位芯の構造と転位の運動に注目し、その関係について考察した。開発した混合転位の原子モデルを用いて、混合転位のパイエルス応力を求めた。さらに、らせん転位のキンク対の形成に必要な活性化エネルギを計算し、その負荷せん断応力依存性を明らかにした。これらの情報は、転位動力学法に転位の運動の温度依存性を実装するための基本的な情報であるため、これらを転位動力学法に実装することで、転位の運動の温度依存性を考慮したシミュレーションが可能になる。
参考文献(1) Y.N. Osetsky, D.J. Bacon, “An atomic-level model for studying the dynamics
of edge dislocations in metals”, Modelling and Simulation in Materials
Science and Engineering, Vol. 11, pp. 427-446 (2003).
(2) G.J. Ackland, M.I. Mendelev, D.J. Srolovitz, S. Han, A.V. Barashev,“Development of an interatomic potential for phosphorus impurities in α-iron”,
Journal of Physics: Condensed Matter, Vol. 16, S2629 (2004).
(3)

A. Stukowski, V.V. Bulatov, A. Arsenlis , “Automated identification andindexing of dislocations in crystal interfaces”, Modellin and Simulation inMaterials Engineering and Science, Vol. 20, p. 085007 (2012).

(4) A. Stukowski, K. Albe, “Extracting dislocations and non-dislocation crystal
defects from atomistic simulation data”, Modelling and Simulation inMaterials Engineering and Science, Vol. 18, p. 085001 (2010).
(5)

G. Henkelman, H. Jonsson, “Improved tangent estimate in the nudged elastic
band method for finding minimum energy paths and saddle points”, Hournal ofChemical Physics, Vol. 113, pp. 9978-9985 (2000).

(6)

Y. Aono, E. Kuramoto, K. Kitajima, Dislocation Dynamics and Plasticity,
Springer, Berlin (1991).

3.1−9

3.1.2 任意の混合転位の運動特性の分子動力学モデリング（再委託先：電力中央研究所）(1) 純 Fe 中の混合転位のモデリング
ここでは、転位動力学法で用いるための BCC 鉄中の混合転位のすべり運動における摩擦
応力と易動度について分子動力学計算を用いて算出した。BCC 鉄中における活動すべり系には 1 111 {110}、 1 111 {112}、 1 111 {123}が挙げられるが、ここでは、 1 111 {112}転位お2
2
22
よび 1 111 {123}転位の計算を行った。図 3.1.2-1 に作成した結晶モデルの概略を示す。モ2デル結晶の大きさは約 10×10×10nm3 であり、転位の進行方向を x 方向、すべり面法線方
向を y 方向とし、y 方向中央位置にすべり面を設定した。混合転位は、東京理科大学高橋研究室により開発されたコードを用い、すべり面上の x 方向中央位置に直線状の転位を 1本導入した。混合転位は転位線とバーガースベクトルのなす角度αにより定義し、0°≦α＜180°の範囲で、 1 111 {112}、 1 111 {123}についてそれぞれ 39、 38 種の角度につい22て計算を行った。このときa=0°が純らせん転位、a=90°が純刃状転位を示す。なお、純
らせん転位については、特に有限温度ではキンク対形成を伴う運動が想定されるが、現状では{110}上においても分子動力学法によるすべり運動の再現が難しい状況であることを考慮し、ここでは計算の対象外とした。また、BCC 金属の{112}上でのすべり変形では、双晶を形成する双晶方向と双晶を形成しない反双晶方向が存在し、両者で変形のしやすさが異なることが知られており(1)、また{123}上のすべり変形においても同様の変形方向の依存性が存在し得る。本研究においても双晶方向と反双晶方向の両方を計算の対象とした。
分子動力学計算は電力中央研究所が所有する大型計算機 SGI ICE X を用いて行った。計
算コードには東京理科大学高橋研究室で開発された分子動力学計算コードに、電力中央研
究所の大型計算機での動作のための修正を加えたコードを用いた。分子動力学計算では、初めに原子数・体積・温度一定の条件（NVT 一定条件）にて結晶を所定の温度に設定した後、温度制御を外し、結晶上下の xz 面それぞれ端部 10Åに対してバーガースベクトル方向（111 方向）に一定せん断応力を負荷して結晶をせん断変形させ、転位を運動させた。z方向には周期境界条件を、x 方向には転位のらせん成分を考慮した周期境界条件を適用した。温度は 100K および 300K とし、計算時間ステップは 2fs でそれぞれ 1.2ns、 0.6ns まで計算を行った。計算により得られた原子座標から Common Neighbor Analysis (CNA)法(2)により転位を抽出し、転位の進行方向に対する転位芯位置を求めた。原子間ポテンシャ

ルには Ackland らによる Embedded Atom Method （EAM）ポテンシャル(3)を用いた。

3.1−10y [112] or [123]

混合転位

x
z

図 3.1.2-1 モデル結晶の概略
図 3.1.2-2 に変形中の混合転位の転位芯位置の変化の例を示した。この条件下においては、転位芯位置は時間に対して概ね直線的に変化していることがわかる。この直線の傾きから各負荷
応力下における転位の速度を求めた。図 3.1.2-3 に負荷応力と転位速度の関係を示す。

350
10MPa
30MPa

300

50MPa
70MPa

250

転位芯位置 [nm]112面双晶方向
温度=100K
角度=15.2°

90MPa

200
150
100
50
0
0.0E+00

2.0E-10

4.0E-10

6.0E-108.0E-101.0E-09

1.2E-09

時間 [s]

図 3.1.2-2 変形中の転位芯位置と時間の関係
300

112面双晶方向
温度=100K
角度=15.2°

転位速度 [m/s]

250

200

150

100

50

0
0

20

40

60

80

100

負荷応力 [MPa]図 3.1.2-3 負荷応力と転位速度の関係

負荷応力と転位速度の関係から、以下の転位速度と応力の関係を用いて転位の易動度と摩擦応力を求めた。

3.1−11v = Mb (t - t f )ここで v は転位速度、M は易動度、b はバーガースベクトルの大きさ、τ はせん断応力、τf は摩擦応力である。摩擦応力はある条件下で転位を運動させるための臨界せん断応力と考えることができ、負荷応力と転位速度の関係において v=0 の時の応力軸との切片が摩擦応力として得られる（図 3.1.2-4）。なお、今回の計算条件においては、転位速度が概ね 20m/s 以上の変形において転位芯位置と時間の間に直線的な関係が認められたため易動度、摩擦応力の計算には転位速度 20m/s 以上の結果を用いた。
160

112面双晶方向
温度=100K
角度=39.2°140

転位速度 [m/s]

120
100

Mb
M:易動度
b:バーガースベクトル

80
60
40

摩擦応力

20
0
0

20

40

60

80

100

負荷応力 [MPa]図 3.1.2-4 易動度および摩擦応力の例
図 3.1.2-5～図 3.1.2-8 に摩擦応力を転位線とバーガースベクトルがなす角度についてまとめた。{112}上で温度 100K での転位運動では、角度 39.2°、90°、140.8°において高い摩擦応力が認められ、パイエルス応力の影響が残存していることが示された。300K においては総じて摩擦応力は低くなった。また{123}上の転位運動においても 100K で摩擦応力が高い混合転位が認められるが、これもパイエルス応力が高い角度にあたり、{123}上の転位運動においても100K ではパイエルス応力の影響が見られることがわかる。一方 300K では摩擦応力は総じて低
く、転位線とバーガースベクトルがなす角度依存性は認められなかった。
50
双晶方向

112面
温度= 100K

45

反双晶方向

40

摩擦応力[MPa]

35
30
25
20
15
10
5
0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-5 摩擦応力と転位線とバーガースベクトルのなす角度との関係（112 面、温度 100K）

3.1−12

50
双晶方向

112面
温度= 300K

45

反双晶方向

40

摩擦応力[MPa]

35
30
25
20
15
10
5
0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-6 摩擦応力と転位線とバーガースベクトルのなす角度との関係（112 面、温度 300K）
60
双晶方向

123面
温度= 100K

50

反双晶方向

摩擦応力[MPa]

40
30
20
10
0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-7 摩擦応力と転位線とバーガースベクトルのなす角度との関係（123 面、温度 100K）
60
双晶報告

123面
温度= 300K

50

反双晶報告

摩擦応力[MPa]

40
30
20
10
0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-8 摩擦応力と転位線とバーガースベクトルのなす角度との関係（123 面、温度 300K）
図 3.1.2-9～図 3.1.2-12 には転位の易動度を転位線とバーガースベクトルがなす角度についてまとめた。{112}上、{123}上ともに 100K においては、角度 20°付近、160°付近において易動

3.1−13度が高い。角度 60°、120°付近では転位が動かないケースが見られ、グラフ中に値は示されないがこれらの混合転位の易動度は非常に小さいことが推察される。一方、300K においては、
90°付近で易動度が高いほぼ単一ピークを持つ角度依存性が認められた。
16000
14000

易動度[m2/Ns]

12000
10000
8000
6000
4000
双晶方向

112面
温度= 100K

2000

反双晶方向

0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-9 転位の易動度と転位線とバーガースベクトルのなす角度との関係（112 面、温度 100K）
16000

112面
温度= 300K

14000

易動度[m2/Ns]

12000
10000
8000
6000
4000
双晶方向

2000

反双晶方向
0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-10 転位の易動度と転位線とバーガースベクトルのなす角度との関係（112 面、温度 300K）
16000
14000

易動度[m2/Ns]

12000
10000
8000
6000
4000
双晶方向

123面
温度= 100K

2000

反双晶方向

0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-11 転位の易動度と転位線とバーガースベクトルのなす角度との関係（123 面、温度 100K）

3.1−14

12000易動度[m2/Ns]

10000
8000
6000
4000
2000

双晶方向

123面
温度= 300K

反双晶方向

0
0

20

40

60

80

100

120

140

160

180

角度 [degree]

図 3.1.2-12 転位の易動度と転位線とバーガースベクトルのなす角度との関係（123 面、温度 300K）

以上の結果をもとに、易動度のモデリングを行った。ここでは転位の易動度は下記の式を用いて記述することとした。
n{}M (q ) = M 0 + å M i (1 - sin q - Qi )i =1

ここでθは転位線とバーガースベクトルのなす角度、n は曲線上の極大値と極小値の個数であり、 Q i は各すべり系で 100K において易動度が極大あるいは極小の近傍となる値を設定した。ここではすべり系、温度によらず n=9 を用いた。これを分子動力学法で得られた値と式で得られた値の二乗誤差が最小となるようにフィッティングを行った。表 3.1.2-1に用いた Q i の値を示す。図 3.1.2-13～図 3.1.2-16 にモデリング結果を示す。いずれの場
合も易動度の角度依存性を概ね再現できていることがわかる。
表 3.1.2-1

Qi 値

i

{112}

{123}

1

0

0

2

15.2252

25.06583

58.5178

61.8745477.7748

76.8377

5

90

90

6

105.2252

115.06587

121.4822

123.7498

164.7748

156.2411

9

180

180

3.1−15

16000

反双晶方向
MD
モデル

14000

モデル

12000

転位の易動度[m^2/Ns]

転位の易動度[m^2/Ns]

14000

16000

双晶方向

MD

10000
8000
6000
4000

12000
10000
8000
6000
4000

2000

2000

0

0

0

10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

混合転位の角度[degree]

混合転位の角度[degree]

図 3.1.2-13 易動度のモデリング結果（112 面、100K）

12000

12000

双晶方向

10000

転位の易動度[m^2/Ns]

転位の易動度[m^2/Ns]

10000

8000
6000
4000
2000

反双晶方向

MD
モデル

8000
6000
4000

MD
モデル

2000
0

0
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

混合転位の角度[degree]

混合転位角度[degree]

図 3.1.2-14 易動度のモデリング結果（112 面、300K）

16000

16000

双晶方向

14000

モデル

12000

転位の易動度[m^2/Ns]

転位の易動度[m^2/Ns]

14000

MD

10000
8000
6000
4000

MD

反双晶方向

モデル

12000
10000
8000
6000
4000
2000

2000

0

0

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

混合転位の角度[degree]

混合転位の角度[degree]

図 3.1.2-15 易動度のモデリング結果（123 面、100K）

3.1−16

12000

12000

双晶方向

反双晶方向
10000

転位の易動度[m^2/Ns]

転位の易動度[m^2/Ns]

10000
8000
6000
4000

8000
6000
4000
2000

2000

MD
モデル

0

0

MD
モデル

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

混合転位の角度[degree]

混合転位角度[degree]

図 3.1.2-16 易動度のモデリング結果（123 面、300K）(2) まとめ{112}、{123}上の混合転位の各温度における易動度と摩擦応力を計算し、易動度の近似式を用いて結果のフィッティングを行った結果、高い再現性を得ることができた。以上により、{112}、{123}上の任意の混合転位の易動度を計算することが可能になった。さらに、近似式を転位動力学法に実装することによって、各温度における混合転位の運動をより定量的にシミュレートすることが可能になる。
参考文献(1) M. Yamaguchi, V. Vitek, “Core structures of non screw 1/2 (111)Temperature
Dependence of Mixed Dislocation Mobility in BCC Iron dislocations on (112)planes in bcc crystals. II. Peierls stresses and the effects of an external
shear stress on the cores,” Journal of Physics F: Metal Physics, Vol. 5, p.
11 (1975).
(2) D. J. Honeycutt, H. C. Andersen, “Molecular dynamics study of melting andfreezing of small Lennard-Jones clusters,” The Journal of Physical Chemistry,
Vol. 91, p. 4950–4963 (1987).
(3) G. J. Ackland, M. I. Mendelev, D. J. Srolovitz, S. H. a. A. V. Barashev,“Development of an interatomic potential for phosphorus impurities in αiron,” Journal of Physics: Condensed Matter, Vol. 16, pp. S2629-S2642 (2004).

3.1−17

3.1.3 合金や不純物元素を含む金属中での転位挙動解析の準備（再委託先：電力中央研究所）(1) 原子間ポテンシャルの検討
原子間ポテンシャルは多原子系のエネルギを計算する関数（ポテンシャル関数）であり、古典分子動力学法では、その微分値を原子間力として用い、 原子位置を更新して系の時
間発展を行う。結合の特性ごとに様々なポテンシャル関数があり、4 種類への分類の仕方がある (1) （図 3.1.3-1）。本研究では、このうち、多元系金属で用いられることの多い
EAM ポテンシャルを中心に調査を進めた。

図 3.1.3-1 原子間ポテンシャルの分類(1)

EAM ポテンシャルの内部エネルギΦは下記のように表現することができる。
1
Φ = % ϕ'( )r'( + + % U' (ρ1' )
2
',(

'

ここで i、j は原子のインデックス、r'( は i 番原子と j 番原子の間の原子間距離、ρ1' は i 番原子における電子密度、φij は 2 体ポテンシャル関数、U' は埋め込み関数である。
ρ1' = % ρ(',(

ここで、ρ( は j 番原子からの電子密度の寄与である。φや U の関数形により様々な EAM
ポテンシャルが存在する。初期の EAM ポテンシャルでは内部エネルギを Rose 関数等で近似し、そこから 2 体ポテンシャルの寄与分を差し引いた値を U とする関数形であったが、φと U を任意の関数で表現し、複雑な形状へとすることにより、徐々に精度が高まっていった。多元系 EAM ポテンシャルでは、基本的には埋め込み関数に異元素間パラメータがないことから、異元素ペア間のφij のポテンシャルパラメータのみを考えればよいことになるため、開発が比較的容易である。ただし、異元素間の電子密度比が考慮される場合もあるが、それも異元素ペア間のみであり、異元素のトリプレットについてはポテンシャル
パラメータを決定する必要はない。U(ρ) = A4ρである Finnis-Sinclair(FS)ポテンシャル(2)も EAM ポテンシャルの一種と考えることができる。
本研究で主な対象としている圧力容器鋼は BCC 構造の Fe を母相とするため、まず、基本となる Fe 単元系のポテンシャルについてまずは調査を行った。Fe 単元系 EAM ポテン
シャルによるき裂の破壊現象に関するレビュー論文(3)に示された代表的な Fe 単元系 EAM
3.1−18

ポテンシャルのリストを表 3.1.3-1 に示す。FS ポテンシャル(2)、 Simonelli ポテンシャル(4)(5)、

Chamati ポテンシャル(6)、MendelevII ポテンシャル(7)、Gordon ポテンシャル(8)、Marinica07 ポテンシャル(9)、Marinica11 ポテンシャル(10)、Chiesa ポテンシャル(11)が比較対象に挙げられている。これに加えて、Mendelev p10(12) (または MendelevIIext)と呼ばれる原子間ポテンシャルが広く使われている。また、Dudarev らによるポテンシャル(13)が言及されることもある。2 band model（2BM）ポテンシャル(14)(15)もあるが、近年ではあまり使われていない。このうち、近年広く使われるようになった並列化古典分子動力学法ソフトウェア LAMMPS コード(16)に搭載されているのは MendelevII、 Mendelev p10 のみである。
表 3.1.3-2 に示されるように、どの Fe 単元系の EAM ポテンシャルも実験値または第一
原理計算値をある程度再現できている。しかしながら、原子間ポテンシャルを古典分子動力学法で使用した場合の有効性を正しく評価することは難しく、現在でも評価方法は確立していない。このため、原子間ポテンシャルの引用数で原子間ポテンシャルの評価がなされることがある。図 3.1.3-2 に 2000 年以降に発表された Fe 単元系の原子間ポテンシャルの論文引用数の推移を示すが、近年では MendelevII、 Mendelev p10 ポテンシャルの引用数が高いことがわかる。このことは、LAMMPS の普及の影響も大きいと考えられる。
Mendelev 系ポテンシャルの特徴として、らせん転位の転位芯構造を再現できることがあり、これが広く使われていることが要因の一つである。図 3.1.3-3 にらせん転位のDifferential Displacement map を示す(17)が、第一原理計算の結果をよく再現できている
ことがわかる。Mendelev p10 と MendelevII の違いは熱膨張挙動である。図 3.1.3-4 に示
すように MendelevII ポテンシャルは低温で負の熱膨張率を示すが、Mendelev p10 ポテン
シャルはこの点が改良され、常に正の熱膨張率を示す(18)。以上のように、Mendelev II、
p10 ポテンシャルはらせん転位への適用性があり、実績も高い。一方、Mendelev ポテンシ
ャルは転位の移動時の Minimum Energy Path が第一原理計算と異なり、分裂することが知られており（Chamati ポテンシャルは分裂しない）（図 3.1.3-5）(19)、この点には留意が必要である。

3.1−19

表 3.1.3-1 代表的な Fe 単元系 EAM ポテンシャルのリスト(3)表 3.1.3-2 代表的な Fe 単元系 EAM ポテンシャルの基礎的な物性(3)3.1−20図 3.1.3-2 代表的な Fe 単元系の原子間ポテンシャル論文の引用数(2000 年以降発表の論文、Scopus(20)による調査)

図 3.1.3-3 らせん転位の第一原理計算と MendelvII ポテンシャルによる Differential
Displacment map(17)図 3.1.3-4 MendelevII、 p10 ポテンシャルの熱膨張率(18)3.1−21図 3.1.3-5 らせん転位移動の最小エネルギ経路(19)次に Fe 合金系のポテンシャルについて記す。Fe の合金系ポテンシャルは主として、任意の元素の組み合わせに用いることのできる原子間ポテンシャルと、 ある元素の組み合
わせの元素について特化した原子間ポテンシャルの 2 種類に分けることができる。これらについてそれぞれ調べた。
（a） Fe を含む任意の元素の組み合わせに用いることのできる原子間ポテンシャル
EAM ポテンシャルを多元系に用いるためには、 異元素の組み合わせの 2 体ポテンシャ
ルを決める必要がある。任意の組み合わせに用いる異元素間の 2 体ポテンシャルは、理論に基づき決定される。このため、ポテンシャル全体としての多元系への適用性は実験値等への合わせ込みによって保証されることはなく、用いた理論の正しさにのみ依存することになる。また、元素間で統一的なポテンシャル関数を用いるので、複雑なポテンシャル関数の導入は今のところできていない。これらで代表的なのが FS ポテンシャル (2) とgeneralized EAM(GEAM)ポテンシャル(21)(22)である。FS ポテンシャルでは、単元素の平均律、 GEAM は Johnson の合金 EAM 理論(23)に基づく。また、MEAM ポテンシャル(24)は厳密には EAM ポテンシャルの枠組みには入らないが、任意の元素の組み合わせに用いることができる。
現状で、これらの Fe を含む原子間ポテンシャルで計算できる組み合わせは下記の通りである。
FS：V、 Nb、 Ta、 Cr、 Mo、 W、 Fe
GEAM：Cu、 Ag、 Au、 Ni、 Pd、 Pt、 Al、 Pb、 Fe、 Mo、 Ta、 W、 Mg、Co、 Ti、 Zr
MEAM：Cu、Ag、Au、Ni、Pd、Pt、Al、Pb、Rh、Ir、Li、Na、K、V、Nb、Ta、Cr、
Mo、W、Fe、C、Si、Ge、H、N、O（b）特定の元素の組み合わせの元素について特化した Fe 系原子間ポテンシャル

3.1−22

特定の元素の組み合わせの Fe 系の原子間ポテンシャルについて調べたところ下記の系について、原子間ポテンシャルが開発されていることがわかった。（ただし、精度が低いと考えられる原子間ポテンシャルはあらかじめ除いてある。）
・Fe-Cr 系：Caro ら* (2005)(25)、 Olson ら✝ (2007)(14)、 Bonny ら*(2011)(15)
・Fe-Ni 系：Bonny ら*(2009)(26)
・Fe-Cu 系：Pasianot ら*(2007)(27)
・Fe-V 系：Mendelev ら✝(2007)(28)
・Fe-P 系：Ackland ら✝(2004)(12)・Fe-Ni-Cr 系：Bonny ら(2011)(29)、 Bonny ら(2013)(30)
・Fe-Cu-Ni 系：Bonny ら*(2009)(31)
・Fe-Cu-Ni-Mn 系：Bonny ら、 Terentyev ら*(2013、2015)(32)(33)（*は Fe-Fe 原子間ポテンシャルが Mendelev II ポテンシャル、 ✝は Mendelev p10 ポテ
ンシャル）
これらのポテンシャルにおいて、Fe-Fe 原子間のポテンシャルには Mendelev II、
Mendelev p10 ポテンシャルが使われていることが多かった。圧力容器鋼あるいは圧力容器鋼を念頭においた Fe モデル合金では、調査した EAM ポテンシャルの中で利用できるのは Fe-Cu、Fe-Cu-Ni、 Fe-Cu-Ni-Mn 系のポテンシャルである。Fe-Cu 系原子間ポテンシ
ャルの開発では、Pasianot、Malerba(27)は、Fe-Cu 原子間ポテンシャル開発時に BCC FeCu bcc 相の混合エネルギと BCC/FCC 2 相共存平衡状態における Cu の Fe への固溶限（図3.1.3-6）に加えて BCC-Fe 中の銅-空孔のバインディングエネルギ（表 3.1.3-3）、空孔ジャンプエネルギ障壁（表 3.1.3-4）が再現されるようにパラメータを決定している。また、Ackland ら(34)、Ludwigy ら(35)の原子間ポテンシャルと比較して、熱力学特性が優れているとしている。Fe-Ni 系原子間ポテンシャルの開発では Bonny らが規則合金の物性（表3.1.3-5）や欠陥のエネルギ（表 3.1.3-6、表 3.1.3-7）を合わせ込んでいる(26)。また、物性の再現性や不規則合金の不安定性から Mayer らの EAM ポテンシャル(36)や Mishin らの
ADP ポテンシャル(37)より優れ、相図も再現できるとしている（図 3.1.3-7、図 3.1.3-8）。
3 元系以上では、従来の原子間ポテンシャルを組み合わせた原子間ポテンシャルもあり、
Fe-Cu-Ni 系原子間ポテンシャルは既存の Fe-Fe(37)、 Fe-Cu(27)、 Fe-Ni(26)、 Cu-Cu(38)、
Ni-Ni(39)系ポテンシャルを組み合わせ、Cu-Ni 系ポテンシャルを追加した原子間ポテンシ
ャルであり、Cu-Ni 系の欠陥のバインディングエネルギ(表 3.1.3-9)と CuNi 合金の生成エ
ネルギ(図 3.1.3-9)、相図(図 3.1.3-10)に合わせ込んでいる。Fe-Cu-Ni-Mn 系原子間ポテン
シャル(32)(33)も Fe-Cu-Ni 系原子間ポテンシャル(30)に Fe-Mn、 Ni-Mn、 Cu-Mn、 Mn-Mn、系原子間ポテンシャルを追加した原子間ポテンシャルである。一方 Fe-Ni-Cr 系原子間ポテ
ンシャル(29)(30)はすべてのポテンシャルパラメータを新たに最適化している。開発者自身による Fe-Cu を含む系における BCC 中の転位に対する適用例がある(40)(41)(42)。また、必ずしも原子間ポテンシャルへの引用だけではないが、引用数もある程度あり、開発者以外による使用例があることもわかる(図 3.1.3-11)。
3.1−23図 3.1.3-6 Fe-Cu 系ポテンシャルの Cu の Fe への固溶限(27)（AB が Ackland らによるポテンシャ
ル(34)、LF が Ludwig らによるポテンシャル(35)、CO5.20 が Pasianot らの原子間ポテンシャル(27)）

表 3.1.3-3 Fe-Cu 系ポテンシャルの Cu-空孔のバインディングエネルギ(27)3.1−24

表 3.1.3-4 Fe-Cu 系ポテンシャルの空孔ジャンプエネルギ障壁(27)表 3.1.3-5 Fe-Ni 系ポテンシャルの FeNi 規則合金の生成エネルギ(31)(MEY は Meyer らの EAM ポ
テンシャル(36)、MISH は Mishin らの ADP ポテンシャル(37))

3.1−25

表 3.1.3-6 Fe-Ni 系ポテンシャルの BCC－Fe 中の欠陥生成エネルギ(26)表 3.1.3-7 Fe-Ni 系ポテンシャルの FCC－Ni 中の欠陥生成エネルギ(26)3.1−26

図 3.1.3-7 Fe-Ni 系ポテンシャルによる FeNi 合金の生成エネルギと規則合金の生成エネルギとの比較(26)

図 3.1.3-8 Fe-Ni 系ポテンシャルによる Fe-Ni 相図(26)表 3.1.3-8 Fe-Cu-Ni 系ポテンシャルによる欠陥のバインディングエネルギ(31)3.1−27

図 3.1.3-9 Fe-Cu-Ni 系ポテンシャルによる Cu-Ni 合金生成エネルギ(31)図 3.1.3-10 Fe-Cu-Ni 系ポテンシャルによる Cu-Ni 合金相図(26)

図 3.1.3-11 Bonny らの BCC 系 Fe 合金ポテンシャルの引用数の推移（Scopus(20)による調査）

3.1−28

以上の調査結果は次のようにまとめることができる。
・Fe 単元系ポテンシャルについては、らせん転位の転位芯構造を再現できることから多く用いられている Ackland ら、 Mendelev らにより開発されたポテンシャルが利用できると考えられる。
・BCC 鉄の合金を対象にした多元系のポテンシャルは多くないが、本研究で主な成果反映
対象として想定している軽水炉圧力容器鋼の照射脆化では、Fe 中の Cu が照射脆化に強く寄与し、加えて Ni、 Si、 Mn といった元素が寄与することが知られており、これに鑑みると、まず、合金系として Fe-Cu2 元系の計算から着手し、将来的には Fe-Cu-Ni３
元系等に対象を広げることが想定される。従って、当研究においては Fe-Cu 系だけでな
く、Fe-Cu-Ni 系への拡張が可能な Bonny らのポテンシャルを用いることが望ましいと考えられる。
(2) Fe 合金中の混合転位のモデリング
本研究では圧力容器鋼の脆化に強く寄与することが知られている Cu を含む 2 元系合金
系から検討を開始することを想定し、Fe-Cu 合金を対象とした。まず、格子定数を求める
ために、原子 16000 個からなる BCC 構造の Fe 結晶中に所定の濃度に相当する数の Cu 原子をランダムに配置したモデル結晶を作成した。このモデル結晶を対象に、Berendesen
thermostat, barostat による原子数、圧力、温度一定条件（NPT 一定条件）下で、時間刻
み 5fs にて、 10000 ステップの分子動力学計算を行った。原子間ポテンシャルには前節で選定した Bonny らによる Fe-Cu-Ni 三元系 EAM ポテンシャル(43)を用いた。温度は 0K〜
1000K とし、温度制御は毎ステップ行った。圧力の制御は 10 ステップ毎に x、 y、 z 軸
方向の平均圧力sxx、 syy、 szz を算出し、平均圧力がゼロとなるように x、 y、 z 軸方向の格子定数を変化させて制御した。平衡格子定数は最終 500 ステップにおける x、 y、 z
方 向 の 格 子 定 数 の 平 均 値 と し た 。 一 般 化 積 層 欠 陥 エ ネ ル ギ は < 111̄ > {110} 、< 111̄ > {112}、 < 111̄ > {123}のすべり系を対象に計算した。Cu 濃度は 0%（純 Fe）、

0.25、 0.5at% とした。基本セルの x、y、z 軸はそれぞれのすべり系に対して < 11̄1 >< 110 >< 1̄12 >、 < 111̄ >< 112 >< 11̄0 >、 < 111̄ >< 123 >< 54̄1 >とした。基本セル内の原子数はそれぞれ、 23328 原子、 23328 原子、 20412 原子であり、 計算対象の組成となるようランダムに Cu 原子を配置した。モデル結晶のサイズは各辺 60~70Åである。モ
デル結晶の y 方向上半分を<111>方向に b/128（b はバーガースベクトルの大きさ）ずつ
128 ステップに分けてずらした。原子の x、z 方向は固定し、 y 方向にだけ移動できるようにし、 ｘ、z 方向には周期境界条件を適用した。原子位置は共役勾配法を行って緩和し、緩和後のエネルギを求めた。弾性定数の計算は、上で算出した 0K における格子定数を用いて原子 16000 個からなる BCC 構造の Fe 結晶中に所定の濃度に相当する数の Cu 原子をランダムに配置して作成したモデル結晶を対象に、 NPT 一定条件で、時間刻み 5fs、 ス
テップ数 10000step の構造緩和を行った。温度は 100K と 300K とし、原子間ポテンシャ
ルには Bonny らによる Fe-Cu-Ni 三元系 EAM ポテンシャル(1)を用いた。この緩和構造にひずみを与える deformation method(44)により弾性定数を計算した. まず、ひずみexx、 eyy、

ezz、 gxy、 gyy、 gzx に相当する変形をセルに与え、 セルにかかる圧力を計算した. ひずみ
3.1−29の値はそれぞれ-0.2%、 -0.1%、 0.1%、 0.2%とした。得られた圧力に、平衡構造（ひずみ 0.0%）の圧力を加えると、 それぞれのひずみ成分に対し、3 種類(ひずみ-0.2%、 0.1%、 0.1%、 0.2%)に対する圧力が得られることになる。このときのexx に対するsxx、

syy、 szz の傾きをそれぞれ C11、 C12、 C13、 eyy に対するsxx、 syy、 szz の傾きをそれ
ぞれ C21、 C22、C23、 ezz に対するsxx、 syy、 szz の傾きをそれぞれ C31、 C32、 C33 とし、 gyz に対するsyz の傾きを C44、 gzx に対するszx の傾きを C55、 gxy に対するsxy の傾
きを C66 とした. 傾きは最小二乗法で求めた。得られた弾性定数から次の式を用いて剛性率の Voigt 平均 GV および体積弾性率の Voigt 平均 Bv を計算した。
GV =1[C11 + C22 + C33
15

BV =

1
[C11 + C22 + C33 + 2 (C12 + C13 + C23 )]
9

C12

C13

C23 + 3 (C44 + C55 + C66 )]ここで、 Cij は弾性定数の ij 成分を示す.
混合転位の運動の計算は主要なすべり系である 1 111 {110}すべりを対象に行った。計算2

方法は 3.1 の BCC Fe を対象にしたものとほぼ同様である。サイズ約 10×10×10nm3 のモデル結晶を対象に、転位の進行方向を x 方向、すべり面法線方向（すなわち 110 方向）をy 方向として、y 方向中央位置にすべり面を設定した。モデル結晶中のすべり面上の x 方向
中心位置に直線状の混合転位を 1 本導入し、転位線とバーガースベクトルのなす角度αが0°≦α＜180°の範囲で 53 種の角度について計算を行った。BCC Fe の時と同様に純らせん転位は計算対象外とした。分子動力学計算は電力中央研究所所有の SGI ICE X を用い、計算コードには大型計算機での動作のための修正を加えたコードを用いた。分子動力学計
算は、NVT アンサンブルにて結晶を所定の温度に設定した後に NVE アンサンブルに切り替え、結晶上下の xz 面それぞれ端部 10Åに対してバーガースベクトル方向（111 方向）に
一定せん断応力を負荷して結晶をせん断変形させ、転位を運動させた。z 方向には周期境界条件を、x 方向にはらせん成分を考慮した周期境界条件を適用した。温度は 100K および 300K とし、計算時間はそれぞれ 1.2ns、 0.6ns、 計算時間ステップは 2fs とした。原子間ポテンシャルは格子定数を求めたときと同じ物を用いた。10ps 間隔で出力した原子座
標から Common Neighbor Analysis （CNA）法により転位を抽出し、転位の進行方向に対する転位芯位置を求めた。
図 3.1.3-12 に各 Cu 濃度の Fe-Cu 合金における格子定数の温度依存性を示す。使用した
原子間ポテンシャルの Fe 部分の元となった原子間ポテンシャル(3)の特性上、格子定数は40K 付近で極小値を示すが、Cu 濃度が 5.0%程度となってもその影響が残ることがわかる。
極小値を示した後は温度が高いほど格子定数が高くなる傾向を示した。図 3.1.3-13 には100K および 300K における格子定数の Cu 濃度依存性を示した。いずれの温度においてもCu 濃度が高いほど格子定数が大きくなり、この濃度範囲においては濃度に対して直線的に格子定数が増加する。図 3.1.3-14 には各 Cu 濃度における弾性定数を示す。表中には体積弾性率 K、剛性率 G を併せて示した。Cu 濃度が 5%以下では、弾性定数の Cu 濃度依存性はほとんど見られないことがわかる。
3.1−30図 3.1.3-15～図 3.1.3-17 には各すべり系における一般化積層欠陥エネルギを示す。
Mendelev らの原子間ポテンシャル(7)を用いて積層欠陥エネルギを計算すると第一原理計算と比較して低い値を示すことが報告されているが(3)、本計算においても純 Fe において同じ結果が得られた。計算領域の小さい第一原理計算では Fe-Cu ランダム固溶体の計算が難し
いため Cu の影響を調べるために本ポテンシャルを用いて Cu の影響を調べたが、Cu 濃度が高いと若干積層欠陥エネルギが低下する傾向が見られるがその差はわずかであり、
1.0%Cu 添加は積層欠陥エネルギにほとんど影響を及ぼさないことが示唆された。
図 3.1.3-18、図 3.1.3-19 にはそれぞれ 100K と 300K における Fe-Cu 合金中の混合転位の摩擦応力を示す。温度 100K では、純 Fe において摩擦応力の高い向きにおいて同じように摩擦応力が高く、またその高さは Cu 濃度が高いほど高くなり、Cu の添加による固溶強化が明確に確認できる。一方、 300K になるとやはり Cu 濃度が高い方が全体的に摩擦応力が高くなる傾向が見られるが、100K と比較するとその程度は低い、図 3.1.3-20、図3.1.3-21 はそれぞれ 100K と 300K における混合転位の易動度を示す。混合転位の易動度は 100K で Cu 濃度が高い方が易動度が低い傾向も見られるが、300K では Cu 濃度によらずほとんど一定である。易動度に対しては 1%程度の Cu 添加であればその影響はほとんど見られないと言える。

2.89

Cu濃度
0.0%
0.25%
0.5%
1.0%
3.0%
5.0%

格子定数（A）

2.88

2.87

2.86

2.85
0

200

400

600

800

1000

温度（K）

図 3.1.3-12 Fe-Cu 合金の格子定数の温度依存性
2.862
100K
2.861

300K

格子定数（A）

2.860
2.859
2.858
2.857
2.856
2.855
0

1

2

3

4

5

6

Cu濃度（at%）

図 3.1.3-13 Fe-Cu 合金の格子定数の Cu 濃度依存性

3.1−31300

250

弾性率[GPa]

200

150

100

50

K[GPa]
C11[GPa]
C44[GPa]00

G[GPa]
C12[GPa]

1

2

3

4

5

Cu濃度

図 3.1.3-14 Fe-Cu 合金の弾性定数

0.7
0.6
0.5

純Fe
Cu=0.25%
Cu=0.5%
Cu=1.0%0.665

エネルギー[J/m2 ]エネルギー[J/m2 ]

0.670

純Fe
Cu=0.25%
Cu=0.5%
Cu=1.0%

0.4
0.3
0.2

0.660

0.655

0.1
0.0

0.650
0.0

0.2

0.4

0.6

0.8

1.0

0.4

0.5

変位量 x/b

0.5

0.6

0.6

0.6

0.6

変位量 x/b(a) 全体(b) ビーク付近の拡大

図 3.1.3-15 Fe-Cu 合金の{110}上の積層欠陥エネルギ
0.9

0.8

0.80

Cu=0.25%

0.79

Cu=0.5%0.7

0.78

Cu=1.0%エネルギー[J/m2 ]エネルギー[J/m2 ]

純Fe

0.6
0.5
0.4
0.3

0.77

0.76
0.75
0.74
0.73

0.2

0.72

0.1

0.71

純Fe
Cu=0.25%
Cu=0.5%
Cu=1.0%

0.70

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.4

0.5

0.5

変位量 x/b

変位量 x/b(a) 全体(b) ビーク付近の拡大
図 3.1.3-16 Fe-Cu 合金の{112}上の積層欠陥エネルギ

3.1−32

0.9

0.79
0.78

エネルギー[J/m2 ]0.7

エネルギー[J/m2 ]

0.80

純Fe
Cu=0.25%
Cu=0.5%
Cu=1.0%

0.8

0.6
0.5
0.4
0.3

0.77

0.76
0.75
0.74
0.73

0.2

0.72

0.1

0.71

0.0

純Fe
Cu=0.25%
Cu=0.5%
Cu=1.0%

0.70
0.0

0.2

0.4

0.6

0.8

1.0

0.4

0.5

0.5

変位量 x/b0.6

変位量 x/b(a) 全体(b) ビーク付近の拡大

図 3.1.3-17 Fe-Cu 合金の{123}上の積層欠陥エネルギ

80

T=100K

70

0%
0.25%
0.5%
1.0%

摩擦応力(MPa)

60

50
40
30

20
10
0
0

30

60

90

120

150

180

角度(degree)

図 3.1.3-18 Fe-Cu 合金中の転位の摩擦応力と転位線とバーガースベクトルのなす角との関係（110 面、温度 100K）

80

T=300K

70

0%
0.25%
0.5%
1.0%

摩擦応力(MPa)

60

50
40
30

20
10
0
0

30

60

90

120

150

180

角度(degree)

図 3.1.3-19 Fe-Cu 合金中の転位の摩擦応力と転位線とバーガースベクトルのなす角との関係（110 面、温度 300K）

3.1−330.6

40000

0%
0.25%
0.5%
1.0%

T=100K

35000

易動度（m2 /Ns）

30000
25000

20000
15000
10000
5000
0
0

30

60

90

120

150

180

角度(degree)

図 3.1.3-20 Fe-Cu 合金中の転位の易動度と転位線とバーガースベクトルのなす角との関係（110 面、温度 100K）

40000

0%
0.25%
0.5%
1.0%

T=300K

35000

易動度（m2 /Ns）

30000
25000

20000
15000
10000
5000
0
0

30

60

90

120

150

180

角度(degree)

図 3.1.3-21 Fe-Cu 合金中の転位の易動度と転位線とバーガースベクトルのなす角との関係（110 面、温度 300K）(3) まとめFe 合金系の原子間ポテンシャルを調査した。まず、基本となる Fe 単元系ポテンシャルについて調査した結果、らせん転位の転位芯構造を再現できることから多く用いられている Ackland ら、Mendelev らにより開発されたポテンシャルが利用できると考えられる。
BCC 鉄の合金を対象にした多元系のポテンシャルは多くないが、例えば３元系では Bonnyらにより Fe-Cu-Ni 型の EAM ポテンシャルが Mendelev らのポテンシャルを元に開発されている。本研究で主な成果反映対象としている軽水炉圧力容器鋼は Cu、Ni、Si、Mn を含む多元系の合金であるが、圧力容器鋼の照射脆化は銅を含むクラスターの形成が主な要因であることから、本研究においては、合金系として Fe-Cu2 元系の計算に着手し、将来的には Fe-Cu-Ni３元系等に拡張することが想定される。従って、当研究においては Fe-Cu
系だけでなく、Fe-Cu-Ni 系への拡張が可能な Bonny らのポテンシャルを用いることが望
ましいと考えられる。
Fe-Cu 合金を対象として、格子定数の温度依存性および Cu 濃度依存性を調べた。格子定
3.1−34数は、40K 付近で最小となることが分かり、Cu 濃度の増加とともに直線的に増加することが分かった。弾性定数の Cu 濃度依存性を調べた結果、Cu 濃度依存性が非常に小さいことが分かった。さらに、{110}、{112}、{123}上における積層欠陥エネルギを計算した。積層欠陥エネルギの情報は、転位と析出物の相互作用モデルの入力情報として重要である。
参考文献(1) A. E. Carlsson, “Beyond Pair Potentials in Elemental Transition Metals andSemiconductors,” Solid State Physics, Vol. 43, pp. 1-91 (1990).
(2) M. W. Finnis, J. E. Sinclair, “A simple empirical N-body potential fortransition metals,” Philosophical Magazine A, Vol. 50, pp. 45-55 (1984).
(3) J. J. Möller,E. Bitzek, “Comparative study of embedded atom potentials foratomistic simulations of fracture in α-iron,” Modelling and Simulation inMaterials Science and Engineering, Vol. 22, p. 045002 (2015).
(4) G. Simonelli, R. Pasianot, E. J. Savino, “Point-defect computer simulation
including angular forces in bcc iron,” Physical Review B, Vol. 50, pp. 727 738 (1994).
(5) G. Simonelli, R. Pasianot, E. J. Savino, “Embedded-atom-method interatomic
potentials for bcc-iron,” Materials Research Society Symposium Proceedings,
Vol. 291, pp. 567-572 (1993).
(6) H. Chamati, N. I. Papanicolaou, Y. Mishin, D. A. Papaconstantopoulos,“Embedded-atom potential for Fe and its application to self-diffusion on
Fe(100),” Surface Science, Vol. 600, p. 1793 (2006).
(7) M. I. Mendelev, S. Han, D. J. Srolovitz, G. J. Ackland, D. Y. Sun, M. Asta,“Development of new interatomic potentials appropriate for crystalline andliquid iron,” Philosophical Magazine, Vol. 83, pp. 3977-3994 (2003).
(8) P. A. Gordon, T. Neeraj, M. I. Mendelev, “Screw dislocation mobility in BCC
Metals: a refined potential description for α-Fe,” Philosophical Magazine,
Vol. 91, pp. 3931-3945 (2011).
(9) L. Malerba, M. C. Marinica, N. Anento, C. Bjorkas, H. Nguyen, C. Domain, F.
Djurabekova, P. Olsson, K. Nordlund, A. Serra, D. Terentyev, F. Willaime , C.
S. Becquart, “Comparison of empirical interatomic potentials for iron applied
to radiation damage studies,” Journal of Nuclear Materials, Vol. 406, pp. 1938 (2010).
(10)
L. Proville, D. Rodney, M. C. Marinica, “Quantum effect on thermally
activated glide of dislocations,” Nature Materials, Vol. 11, pp. 845-849(2012).
(11)
S. Chiesa, P. M. Derlet, S. L. Dudarev, H. V. Swygenhoven, “Optimizationof the magnetic potential for α-Fe,” Journal of Physics: Condensed Matter,
Vol. 23, p. 206001 (2011).
(12)
G. J. Ackland, M. I. Mendelev, D. J. Srolovitz, S. H. a. A. V. Barashev,“Development of an interatomic potential for phosphorus impurities in αiron,” Journal of Physics: Condensed Matter, Vol. 16, pp. S2629-S2642 (2004).
(13)
S. L. Dudarev, P. M. Derlet, “A 'magnetic' interatomic potential formolecular dynamics simulations,” Journal of Physics: Condensed Matter, Vol.
17, pp. 7097-7118 (2005).
(14)
P. Olsson, J. Wallenius, C. Domain, K. Nordlund, L. Malerba, “Two-bandmodeling of α-prime phase formation in Fe-Cr,” Physical Review B, Vol. 72,
p. 214119 (2005).

3.1−35(15)G. Bonny, R. C. Pasianot, D. Terentyev, L. Malerba, “Iron chromium
potential to model high-chromium ferritic alloys,” Philosophical Magazine,
Vol. 91, pp. 1724-1746 (2011).
(16)
S. J. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular
Dynamics,” Journal of Computational Physics, Vol. 117, pp. 1-19 (1995).
(17)
G. Monnet, C. Domain, S. Queyreau, S. Naamane, B. Devincre, “Atomic
and dislocation dynamics simulations of plastic deformation in reactor pressure
vessel steel,” Journal of Nuclear Materials, Vol. 394, pp. 174-181 (2009).
(18)
S. Chiesa, P. M. Derlet, S. L. Dudarev, “Free energy of a <110> dumbbell
interstitial defect in bcc Fe: Harmonic and anharmonic contributions,” Physcal
Review B, Vol. 79, p. 214109 (2009).
(19)
木下惠介, 下川智嗣, 喜成年泰, 澤田英明, 川上和人 , 潮田浩作, “非すべり応力がらせん転位のパイエルスエネルギに与える影響,” 日本機械学会論文集, Vol.
80, p. CM0018 (2014).
(20)“Scopus,” Elsevier, https://www.scopus.com/.
(21)
X. W. Zhou, H. N. G. Wadley, R. A. Johnson, D. J. Larson, N. Tabat, A.
Cerezo, A. K. Petford-Long, G. D. W. Smith, P. H. Clifton, R. L. Martens, T.
F. Kelly, “Atomic scale structure of sputtered metal multilayers,” Acta
Materialia, Vol. 49, pp. 4005-4015 (2001).
(22)
X. W. Zhou, R. A. Johnson, H. N. G. Wadley, “Misfit-energy-increasing
dislocations in vapor-deposited CoFe/NiFe multilayers,” Phisical Review B,
Vol. 69, p. 144113 (2004).
(23)
R. A. Johnson, “Alloy models with the embedded-atom method,” PhysicalReview B, Vol. 39, pp. 12554-12559 (1989).
(24)
M. I. Baskes, “Modified embedded-atom potentials for cubic materials
and impurities,” Physical Review B, Vol. 46, pp. 2727-2742 (1992).
(25)
Caro, D. A. Crowson, M. Caro, “Classical Many-Body Potential forConcentrated Alloys and the Inversion of Order in Iron-Chromium Alloys,”
Physical Review Letters, Vol. 95, p. 075702 (2005).
(26)
G. Bonny, R. C. Pasianot, L. Malerba, “Fe-Ni many-body potential for
metallurgical applications,” Modelling and Simulation in Materials Scienceand Engineering, Vol. 17, p. 025010 (2009).
(27)
R. C. Pasianot, L. Malerba, “Interatomic potentials consistent withthermodynamics:The Fe-Cu system,” Journal of Nuclear Materials, Vol. 360, pp.
118-127 (2007).
(28)
M. I. Mendelev, S. Han, W. Son, G. J. Ackland, D. J. Srolovitz,“Simulation of the interaction between Fe impurities and point defects in V,”
Physical Review B, Vol. 76, p. 214105 (2007).
(29)
G. Bonny, D. Terentyev, R. C. Pasianot, S. P. a. A. Bakaev, “Interatomic
potential to study plasticity in stainless steels: the FeNiCr model alloy,”
Modelling and Simulation in Materials Science and Engineering, Vol. 19, p.
085008 (2011).
(30)
G. Bonny, N. Castin, D. Terentyev, “Interatomic potential for studying
ageing under irradiation in stainless steels: the FeNiCr model alloy,”
Modelling and Simulation in Materials Science and Engineering, Vol.21, p.
085004 (2013).
(31)
G. Bonny, R. C. Pasianot, N. Castin, L. Malerba, “Ternary Fe-Cu-Ni
many-body potential to model reactor pressure vessel steels: First validation
by simulated thermal annealing,” Philosophical Magazine, Vol. 89, pp. 35313546 (2009).

3.1−36(32)G. Bonny, D. Terentyev, A. Bakaev, E. E. Zhurkin, M. Hou, D. V. Neck,
L. Malerba, “On the thermal stability of late blooming phases in reactor
pressure vessel steels: An atomistic study,” Journal of Nuclear Materials,
Vol. 442, pp. 282-291 (2013).
(33)
D. Terentyev, X. He, G. Bonny, A. Bakaev, E. Zhurkin, L. Malerba,“Hardening due to dislocation loop damage in RPV model alloys: Role of Mn
segregation,” Journal of Nuclear Materials, Vol. 457, pp. 173-181 (2015).
(34)
G. J. Ackland, D. J. Bacon, A. F. Calder, T. Harry, “Computer simulationof point defect properties in dilute Fe-Cu alloy using a many-body interatomic
potential,” Philosophical Magazine A, Vol. 75, pp. 713-732 (1997).
(35)
M. Ludwigy, D. Farkasz, D. Pedrazax, S. Schmaudery, “Embedded atom
potential for Fe-Cu interactions and simulations of precipitate-matrix
interfaces,” Modelling and Simulation in Materials Science and Engineering,Vol. 6, p. 19–28 (1998).
(36)
R. Meyer, P. Entel, “Martensite-austenite transition and phonon
dispersion curves of Fe12xNix studied by molecular-dynamics simulations,”
Physical Review B, Vol. 57, p. 5140 (1998).
(37)
Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, “Phase stability in
the Fe3Ni system: Investigation by first-principles calculations and atomistic
simulations,” Acta Materialia, Vol. 53, pp. 4029-4041 (2005).
(38)
Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, A. F. Voter, J. D.
Kress, “Structural stability and lattice defects in copper: Ab initio; tightbinding; and embedded-atom calculations,” Physical Review B, Vol. 63, p.
224106 (2001).
(39)
A. F. Chen, S. P. Voter, “Accurate Interatomic Potentials for Ni; Al
and Ni3Al,” Materials Research Society Symposium Proceedings, Vol. 82, p. 175(1987).
(40)
G. Bonny, D. Terentyev, E. E. Zhurkin, L. Malerba, “Monte Carlo studyof decorated dislocation loops in FeNiMnCu model Alloys,” Journal of Nuclear
Materials, Vol. 452, pp. 486-492 (2014).
(41)
D. Terentyev, L. Malerba, “Interaction of a screw dislocation with Cuprecipitates, nanovoids and Cu–vacancy clusters in BCC iron,” Journal of
Nuclear Materials, Vol. 421, pp. 32-38 (2012).
(42)
D. Terentyev, L. Malerba, G. Bonny, A. Al-Motasem , M. Posselt,“Interaction of an edge dislocation with Cu–Ni-vacancy clusters in bcc iron,”
Journal of Nuclear Materials, 419, p. 134–139, 2011.
(43)
G. Bonny, R. C. Pasianot, N. Castin, L. Malerba, “Ternary Fe-Cu-Ni
many-body potential to model reactor pressure vessel steels: First validation
by simulated thermal annealing,” Philosophical Magazine, Vol. 89, pp. 35313546 (2009).
(44)
L. A. Girifalco, V. G. Weizer, “Application of the Morse Potential
Function to Cubic Metals,” Physical Review, Vol. 114, pp. 687-690 (1959).

3.1−37

3.2 き裂先端での転位の放出・運動・弾性場を考慮した破壊靭性値の評価のための計算科学的基盤技術の構築（H28〜H29）
3.2.1 転位の弾性場を考慮した破壊靭性値の評価法の確立(1) 転位を用いたき裂の計算方法の高精度化（H28）
き裂を表現する方法の一つに、き裂面上に転位密度を分布させ、き裂面位置において伝達応力が０になるように転位密度分布を求める方法がある。研究代表者らはこれまでに、転位密度分布を離散転位に置き換え、転位動力学法を用いてき裂面位置において伝達応力が０になる離散転位の配置を求める方法を開発してきた(1)。貫通き裂を対象とするため、離散転位を用いてき裂を表現する方法を拡張し、無限体中の貫通き裂の計算を可能にした。
離散転位を用いてき裂を表現する方法では、き裂の形状（外形）を表現する転位を配置し、固定する。その内側に転位（ループ）を置き、転位動力学法のアルゴリズムによって、一様な外部応力と転位間の相互作用応力が釣り合う位置まで転位を移動させる。この操作を必要な回数繰り返し、転位の開口や転位周りの応力分布を得ることができる方法である。
すなわち、本計算方法は、き裂を表現する特殊な転位の配置に対して転位動力学法を直接適用する方法である。転位動力学法の核となるのは、転位に作用する力を転位に作用する
局所的な応力からピーチケーラーの式を用いて求めることであり、計算コストのほとんどがその力の計算に費やされる。破壊力学解析において必要となる応力拡大係数の計算は、き裂の形状を表現するために配置し固定した転位に作用する力から計算することが可能である。き裂の形状を表現する転位以外の内部の転位は、釣り合い位置に移動しているため、転位に作用する力（合力）は０となっている。一方、き裂の形状を表現する転位は、き裂の形状となるように固定されているため、内部の転位が作る応力と、一様に分布する応力による力を受けている。この力は転位という特異線に作用する力であり、J 積分と同義の力である。そのため、線形弾性の範囲においては、この力を応力拡大係数に変換することが可能であり、本研究においても、このような力の変換によって応力拡大係数を計算した。
転位に作用する力の計算は転位動力学法に含まれる計算であるため、固定した転位に作用する力の計算のために新たな機能を追加する必要はなかった。これまでに開発してきた離散転位を用いたき裂の計算方法は、材料内部の楕円形き裂に対してのみ適用してきた。本
研究では、まず貫通き裂を計算対象とするため、き裂を表現する転位の取り扱いや境界条
件について修正する必要がある。具体的には、無限体中の貫通き裂を表現するため、き裂が伸びる方向に周期境界条件を適用した。また、楕円き裂を取り扱う場合には、楕円形状の転位をき裂を表現する転位として用いたが、周期境界条件によって仮定する無限に長い
直線形状の転位を取り扱うように修正を行った。き裂内部の転位の発生については、き裂の中央に、転位源とする点を配置し、その位置での応力が閾値を超えた場合に、１対の直線転位を配置することによって実現した。
開発した計算方法の精度評価のため、転位に設定するバーガースベクトルの大きさを変化させ、最終的に計算される応力拡大係数を理論解と比較を行った。図 3.2.1-1 に計算に用いたモデルを示す。図 3.2.1-1(a)に示すように、幅 6μm の貫通き裂を仮定し、その先端に固定した転位を配置した。き裂面に対して垂直方向に 100MPa の引っ張り応力を与え、き裂先端での応力拡大係数を計算した。せん断弾性定数は 26.5GPa、ポアソン比は 0.347 と3.2−1

した。き裂の離散転位を用いた数値解析では、転位に与えるバーガースベクトルは、き裂の開口量と関係するため、結晶の塑性変形解析の転位動力学シミュレーションと異なり任意の大きさを与えることができる。図 3.2.1-2 に計算に用いたバーガースベクトルの大きさに対する必要となった転位の数を示す。

(a)計算モデル(b) き裂形状

図 3.2.1-1 貫通き裂問題の計算モデル

図 3.2.1-2 貫通き裂の数値解析におけるバーガースベクトルと転位の数の関係
図中の発生は、き裂の離散転位を用いた計算を行うために、き裂の中央に転位源を配置し、転位の生成のためのせん断応力の閾値を用いて計算した結果であり、対消滅は、計算初期
3.2−2

に多くの転位をあらかじめ配置し、転位動力学計算の過程で、き裂の中央で左右のき裂を表す転位が互いに対消滅により消える方法で計算した結果である。図に示すように、バーガースベクトルの大きさが小さいほど、必要となる転位の数が増えていることが分かる。
また、図中の点線は、本貫通き裂問題の開口量 Dmax から各バーガースベクトルの大きさの場合において必要となる転位の数 Nth を式（3.2.1-1）のように理論的に求めた結果を示す。
Nth =

Dmax
2(1 ⌫) a=bµb(3.2.1-1)

ここで、νはポアソン比、μはせん断弾性定数、σは与えた引っ張り応力、a はき裂の半幅、b はバーガースベクトルの大きさである。図より、数値計算によって得られた計算に必要な転位の数は、式（3.2.1-1）で理論的に計算される値とほぼ一致していることが分かる。しかし、いずれの計算結果も、転位の数の変化は階段状であり、理論的に求めた転位の数のように滑らかな変化になっていない。これは、転位が持つバーガースベクトルが一定であるため、段階的に必要となる転位の数が変化することと、実際の計算における転位の数は整数であるため、転位の数は滑らかに変化することができないことが原因である。
図 3.2.1-3 に計算した応力拡大係数を式（3.2.1-2）で与えられる応力拡大係数の理論解

Kth で除した結果を示す。
Kth =

p
⇡a

(3.2.1-2)

図 3.2.1-3 異なる大きさのバーガースベクトルを用いた貫通き裂の応力拡大係数の計算
結果
図中の縦軸の値が 1 のとき、理論解と一致した解であることを示している。計算結果を見3.2−3

ると、おおよそ 1 に近い値を取っていることが分かる。全体的な傾向として、バーガース
ベクトルの大きさを小さくするほど、解の精度が向上する傾向にあることが分かる。本き裂の離散転位を用いた計算方法は、転位の位置でき裂面での面力の条件（伝達応力がゼロ）を満たす選点法であるため、転位の数が増えるほど条件を満たす場所が増えることになる。
すなわち転位の数が増えることによって、き裂面の面力の条件をより厳密に満たすことが可能となるため、バーガースベクトルの大きさを小さくし、必要となる転位の数を増やす
ことによって、応力拡大係数の解の精度が向上したと理解することができる。しかし、バーガースベクトルの大きさを小さくし転位の数が増えると計算量が増し、その結果長い計算時間が必要となる。したがって、今後の大規模計算の実施に向けては、計算コストとのバランスを考慮し、計算に用いるバーガースベクトルの大きさを設定する必要がある。
破壊靭性値の評価のための計算においては、ひずみ速度（応力速度）を与えた計算を行う必要がある。そこで、転位を用いてき裂を表現する方法の応力速度の境界条件を与えた場合の精度について評価を行なった。計算モデルは、バーガースベクトルを変化させて応力拡大係数の精度について調べるために用いた計算モデルと同様である。応力速度は、0.1MPa/ns と 0.5MPa/ns の２通りについて計算を行い、応力速度の違いによる応力各係数の計算精度への影響についても検討を行った。計算結果を図 3.2.1-4 に示す。

(a)応力速度: 0.5MPa/ns(b)応力速度：0.1MPa/ns

図 3.2.1-4 バーガースベクトルを一定として計算した応力速度によるき裂解析結果
図中の直線は、式（3.2.1-2）を用いて解析的に計算した応力拡大係数の理論解である。
応力速度が 0.5MPa/ns の場合を見ると、計算初期に転位を用いたき裂の数値解析結果が理論解を下回る結果を与えている。その後応力が大きくなると、理論解に近づくが階段状に
応力拡大係数が増加していく傾向が見られる。一方、応力速度が 0.1MPa/ns の場合は、さらに階段形状の応力拡大係数の計算結果の変化の傾向が強調され完全に理論解周りを振動する結果となった。これは、計算において転位に与えているバーガースベクトルの大きさが一定値であるため、転位の発生とともにステップ状に応力拡大係数が増加することが原因であると考えられる。そこで、本計算の精度向上を目的として、必要となる転位の数が各応力レベルにおいて一定となるような計算方法を構築した。具体的には、転位に与える
バーガースベクトルを与えられている応力と、き裂の開口量の見積もりから変化させ、必要となる転位の数を一定とする方法である。具体的には、式（3.2.1-1）によって見積も3.2−4

られる計算に必要な転位数が一定となるように、式中のバーガースベクトルの大きさを計算し、各転位に与えるバーガースベクトルの大きさを変化させる方法である。図 3.2.1-5に、応力拡大係数の計算結果を示す。図中の A は、これまでに用いてきたバーガースベク
トルを一定の大きさとして、転位の数を変化させながら計算した結果である。一方、B-8は、き裂の片側の転位の数を 8 本（全部で 16 本）に制御して計算した結果である。

図 3.2.1-5 バーガースベクトルの大きさを可変とした時の応力拡大係数の計算結果
図より、バーガースベクトルを可変として、転位の数を制御した計算結果は、ほぼ直線的に応力の上昇に伴って応力拡大計算が増加している。さらに解析結果とよく一致した結果を得ることができた。バーガースベクトルを一定として計算した場合と比較しても、計算精度の安定性を飛躍的に向上させることに成功した。図 3.2.1-6 にバーガースベクトルの大きさを可変として、制御する転位の数を変化させた結果を示す。B-X は、転位の数をき裂の片側に X に制御した結果を示している。

3.2−5

図 3.2.1-6 転位の数の制御数を変化させた応力拡大係数の計算結果
図を見ると、転位の数が 12 の時に最も理論解から離れた応力拡大計算を計算していることが見られる。これは、転位の数が多くなったことにより、緩和のために必要となる時間が増したため、応力の増加に転位の緩和状態が追いついていないためであると考えられる。
しかし、いずれの転位の数においても、理論解と比べて良い精度で計算できているため、本計算方法の有効性を確認することができた。また、計算速度について計測した結果を図3.2.1-7 に示す。

図 3.2.1-7 転位を用いたき裂解析の計算時間

3.2−6バーガースベクトルを一定とした場合は、応力が大きくなるにつれ必要となる転位の数が増えるため、計算量が増し、その結果計算時間が他のモデルに比べて長くなっていることが分かる。一方、転位の数を制御した場合は、計算時間が短縮されていることが分かる。
B-12 の場合については、モデル A よりも長い時間となっている。したがって、制御する転位の数によっては、従来法に比べて計算時間が長くなることを示しているが、図 3.2.1-6を見ても分かるように、モデル B-8 でも十分な精度の計算が可能となっているため、バーガースベクトルの大きさを可変とすることで、計算精度のみではなく、計算速度も向上させることができることが確認できた。
引っ張り応力が大きくなると、図 3.2.1-8 に示すような計算の破綻（不安定）が発生するケースが見られた。これは、応力が大きくなるにつれき裂先端付近の転位間の距離が非常に短くなり、その結果、転位の運動における微小な数値的な振動によって、転位の位置が入れ子になり、き裂内部の転位がき裂の形状を表す固定された転位を超えて漏れ出てしまう現象である。この計算の不安定現象に対するバーガースベクトルを可変とする方法の効果について確認を行った。その結果を図 3.2.1-9 に示す。

図 3.2.1-8 き裂の転位を用いた数値解析における不安定現象

3.2−7

図 3.2.1-9 応力レベルと計算の不安定現象の発生
図中で、応力拡大係数が 0 に落ちる挙動が計算の不安定現象を示している。バーガースベ
クトルを一定とした場合とモデル B-12 においては、応力がおよそ 60MPa で計算が不安定となり、計算が停止していることが分かる。一方、転位の制御数を小さくすると、より大きな応力まで計算を継続できていることが分かる。これは、転位の数が少ない場合においては、バーガースベクトルの大きさが応力の上昇に伴い大きくなり、き裂先端付近の転位間の距離を転位の数が多い場合に比べて広げる効果があるためである。すなわち、少ない転位の数で精度の高い計算を実現することは、計算の安定性にも良い影響を与えることがわかった。
(2) 既存研究の調査（H29）
転位の遮蔽効果を考慮した延性―脆性遷移のメカニズム研究に関する文献を調査した。
Tanaka らは、鉄単結晶で作られたき裂付き試験片の 4 点曲げ試験を行い、温度毎の破壊靭性値を測定した(2)。さらに、き裂先端からの転位の放出、移動、き裂先端の遮蔽効果を考慮した二次元モデルを構築し、理論解析も実施した。その結果、破壊靭性値が急激に上昇する臨界温度について実験結果と解析結果はよく一致していることから、転位の放出、移動、き裂先端の遮蔽効果によって延性―脆性遷移を表現可能であることがわかった。本研究においても、き裂先端からの転位の放出、移動、き裂先端の遮蔽効果を三次元的に取り扱っているため、Tanaka らの研究と同様に、延性―脆性遷移を調べる方法として妥当であることがわかる。さらに本研究では、き裂と転位を三次元的に取り扱っているため、転位と析出物や他の欠陥との相互作用を陽にシミュレーションに取り込むことができる。したがって、本研究のモデルによって、照射や熱時効による破壊靭性値への影響の評価を行う
ことが可能である。

3.2−8

(3) き裂先端からの転位の放出モデルの開発と、転位による遮蔽効果の計算（H29）
き裂先端を含むすべり面上へのき裂先端からの転位の放出モデルの実装を行なった。図3.2.1-10 に転位の放出モデルの概略を示す。き裂先端の前方のすべり面上に転位源（図3.2.1-10 中で緑の点）を仮定した点を配置する。転位動力学法の各時間ステップにおいて、転位源の位置における応力を計算する。その応力から放出を仮定する転位に対する分解せん断応力を計算し、その分解せん断応力があらかじめ設定した臨界値を超えていた場合には、き裂先端から転位が放出したことを仮定し、新たに転位をすべり面上に配置する。新たに配置する転位は、両端の節点がき裂先端位置に配置し、転位源を囲むような半円形状とした。放出された転位は、転位動力学法によって、転位に作用する力によって移動する。
この時、半円の両端の節点は、き裂先端に沿って移動するように変位を拘束した。本研究では、き裂先端からの転位の放出を仮定する分解せん断応力の臨界値として、新たに配置する半円形状の転位の半径から計算されるフランク・リード源の活性化の臨界せん断応力の値を用いた。シミュレーションにおいては、き裂先端に沿って転位源を複数配置し、様々な位置からの転位の放出を取り扱う。複数の転位が放出され移動すると、き裂先端に沿って移動する転位端部の節点が、他の転位の端部の節点と衝突する。この場合は、図3.2.1-11 に示すように、転位端部の節点を削除し、２つの転位を繋ぐことによって転位の合体を仮定した。放出された転位が作る応力を転位動力学法で計算することによって、放出された転位間、き裂を構成する転位と放出された転位間の応力による弾性相互作用を計算することが可能である。すなわち、このき裂を構成する転位と放出された転位間の弾性相互作用が、転位によるき裂先端の遮蔽効果であり、転位動力学法を用いた本計算手法によって考慮することが可能である。

図 3.2.1-10 き裂先端を含むすべり面上へのき裂先端からの転位の放出モデル

3.2−9

図 3.2.1-11 き裂先端から放出された転位同士の合体開発した転位の放出モデルを実装したき裂の転位動力学コードを用いて、転位によるき裂先端の遮蔽効果のシミュレーションを実施した。図 3.2.1-12 に示すように、中央部に平坦なき裂を有する弾性体の計算を行なった。図に示すようにき裂は転位を用いて表現する。
モデルの y 方向には周期境界条件を与え、無限に長いき裂を仮定し、その他の方向には、
境界条件は与えず、開放系とすることによって無限体を仮定した。弾性体中に一様な引っ張り応力を応力速度 0.1MPa/ns によって作用させた。き裂の幅は 300μm とした。また、転位の運動に対する温度の影響を擬似的に模擬するため転位の易動度を 100 から 10000/Pasに変化させた。一般的に、転位の運動に対する温度依存性が強いらせん転位においては、温度が上昇することによって、易動度が大きくなることが知られている。

図 3.2.1-12 転位の遮蔽効果の解析モデル
図 3.2.1-13 に転位の放出と放出された転位の形状を示す。き裂先端から放出した半円形状の転位は、隣接する転位と合体することによって、最終的に一つの無限に長い直線形状の転位になり、き裂先端から遠方に向けて移動をする。このような転位の放出が連続的に発生し、き裂先端前方に直線形状の転位が並びながら移動していることがわかる。図3.2.1-14 に、横軸に与えた外力から計算されるき裂先端の応力拡大係数の理論解（Kapp）
を応力拡大係数の臨界値（KIC）で無次元化した値、縦軸にシミュレーションによって計算されたき裂先端での応力拡大係数（Ktip）を KIC で無次元化した値をとったグラフを示す。
3.2−10ここで、破壊の開始を表す KIC は、き裂面の表面エネルギを考慮した Griffith の破壊靭性値を用いた。図を見ると、応力が小さい場合には、どの易動度を用いた場合も同じ結果となっている。これは、き裂先端からの転位の放出に十分な応力に達していないため、転位が放出されず、き裂の引っ張りのみの計算となっているためである。Kapp/KIC がおよそ 0.5になると、き裂先端から転位が放出され、き裂先端が放出された転位が作る応力によって遮蔽されることによって、き裂先端での応力拡大係数の計算値である Ktip の値に変化がみられる。いずれの易動度を用いた計算結果においても、転位が放出される以前に比べて、転位が放出された後は、Ktip の変化の勾配は小さくなっていることから、転位によるき裂
先端の遮蔽効果によって、き裂先端の応力拡大係数が理論解に比べて小さくなっていることがわかる。また、設定した転位の易動度が大きくなるにつれて、転位による遮蔽効果が強くなり、Ktip が小さくなっていくことがわかる。本研究では、計算されたき裂先端での応力拡大係数 Ktip が応力拡大係数の臨界値 KIC と等しくなる時、すなわち Ktip/KIC が 1 の時の Kapp の値を、転位によるき裂先端の遮蔽効果を考慮した破壊靭性値 Kfrac とする。図3.2.1-14 から求められる各易動度に対する Kfrac をまとめたのが、図 3.2.1-15 である。図より易動度が大きくなるにつれて Kfrac の値が大きくなっていることがわかる。既述の通り、
易動度が大きくなることは、温度が高くなることを模擬しているので、図 3.2.1-15 は温度が高くなると Kfrac が大きくなることを意味している。通常鉄鋼材料は、低温で脆性を示し、破壊靭性値が小さくなることから、本研究の結果は、鉄鋼材料の破壊靭性値の温度依存性
を定性的に表現している。

図 3.2.1-13 き裂先端から放出された転位

3.2−11図 3.2.1-14 転位によるき裂先端の遮蔽効果の計算結果

図 3.2.1-15 転位の易動度と破壊靭性値の関係
3.2−12

(4) 転位の放出モデルの拡張（H29）
ここまでで実装したき裂先端からの転位の放出モデルは、放出される転位のすべり面がき裂先端を含む場合にのみ適用可能である。しかし、現実にはそのようなすべり面が存在しない場合や、それ以外のすべり面に転位が放出されることが想定される。そこで、本研
究では、図 3.2.1-10 に示した、き裂先端からの転位の放出モデルを応用し、き裂先端から任意のすべり面への転位の放出を取り扱うモデルを開発した。図 3.2.1-16 に示すように、任意のすべり面を取り扱う場合、すべり面はき裂先端とは点で交差し、き裂面とは線で交差する。そのため、分解せん断応力を計算する転位源は、すべり面とき裂先端が交差する点からき裂前方かつすべり面上に配置する。また、放出された転位は、形状を円形に近い扇型とし、すべり面とき裂面が交差する線上に転位の両端の節点を配置する。両端の節点の移動は、き裂面から出ている転位であることを維持するために、すべり面とき裂面が交差する線に沿って移動させる。

図 3.2.1-16 き裂先端から任意のすべり面への転位の放出モデル
き裂先端から任意のすべり面への転位の放出モデルの試計算を実施した。計算条件は、
図 3.2.1-12 のモデルを用いた転位の遮蔽効果の数値解析と同様であるが、き裂先端から放出される転位のすべり面として、き裂面に対し 45°傾いたすべり面を設定した。図 3.2.117 に計算結果を示す。き裂先端を含むすべり面への転位の放出では、放出された転位は、隣接する転位と合体し、直線になりながら進む様子が見られたが、すべり面とき裂先端と点で交差する本計算結果では、放出された転位が、き裂先端を取り囲み、すべり面とき裂
先端の交差する点を中心に広がっていく様子が見られる。すなわち、すべり面の方位によって、放出される転位の運動が異なっていることがわかった。

3.2−13図 3.2.1-17 き裂面から 45°傾いたすべり面への転位の放出を伴う数値解析の結果(5) 転位の運動の温度依存性の実装と破壊靭性値への影響の評価（H30）
平成 30 年度では、「3.1 転位の運動特性の分子動力学モデリング」において、分子動力学法の結果から経験的に導いた転位の易動度の式をき裂問題の転位動力学法に実装し、転位の動きの温度依存性および転位の方位の依存性を考慮することを可能にした。また、本研究で開発したき裂先端からの転位の放出モデルにおいては、き裂の前方に配置した転位源の位置におけるせん断応力が閾値を超えた場合に転位を放出させるとしたが、き裂先端での応力拡大係数および温度の依存性を考慮可能な Gordon らによって提案されたき裂先端からの転位の放出条件(3)を実装した。図 3.2.1-12 に示した、モード I き裂の問題を、温度の設定を 50〜300K に変化させて計算した。計算によって得られたき裂先端での転位の放出と、放出された転位の動きを示すスナップショットを図 3.2.1-18 に示す。温度が 50K の場
合の転位の形状は、き裂先端から放出後、台形形状を有していることがわかる。一方、温度が上昇するについれて、台形形状から曲率を有する形状に遷移することがわかる。温度が 300K の場合においても、直線を維持している転位の方向があることがわかる。この方向の転位はらせん転位であり、らせん転位はその他の転位に比べて非常に易動度が低いため直線形状を維持している。低温の場合は、らせん転位に加えて、その他二方向においても直線形状を有する傾向が見られる。これは「3.1 転位の運動特性の分子動力学モデリング」
において、明らかになったように、らせん転位以外の転位においても、易動度が低い方向があり、その方向の転位が直線形状を維持する傾向を示している。

3.2−14図 3.2.1-18 き裂先端から放出された転位の挙動。(a) 50K (b) 100K (c) 200K (d) 300K
図 3.2.1-19 に図 3.2.1-14 と同様に、き裂先端での応力拡大係数の変化、すなわち転位による遮蔽効果を示す。Kapp/KIc がおよそ 0.8 以下の領域では、線形的に応力拡大係数が上昇しているが、およそ 0.8 以降から曲がり始め、勾配が徐々に小さくなる傾向が見られる。
これは、Kapp/KIc がおよそ 0.8 付近からき裂先端からの転位の放出が開始し、転位による遮蔽の影響の結果、Ktip/KIc が低下している。また、温度が高いほど、Ktip/KIc の低下が著しい。
図中水平に引かれた黒線は、Ktip/KIc が 1 の時、すなわち Ktip が KIc と等しい時を示している。本研究では、KIc は Griffith の条件から計算した。破壊靭性値は、Ktip/KIc が 1.0 の時の Kapp の値として計算される。図 3.2.1-20 に各温度における破壊靭性値の計算結果を示す。
比較のため、Tanaka らによって実施された実験による結果(2)も示した。本研究においては、
温度が 140K までにおいては、Ktip/KIc が 1.0 を超えたため、破壊靭性値を計算することが可能であったが、160K 以上の温度においては、転位による遮蔽の影響が強く、Ktip/KIc が 1.0に至ることなくさらに低下していった。すなわち、160K 以上の温度においては、脆性的ではない破壊が起こることを意味している。よって、本計算結果では、延性―脆性遷移温度は 140K から 160K の間であることが得られた。図中に示した実験結果においても、低温側で脆性破壊が観察され、温度が上昇すると高ひずみを伴う破壊に遷移している。実験結果
では、延性―脆性遷移温度をおよそ点線で示した 150K 程度としており、本計算結果と定量的にも一致する温度となっている。すなわち、本研究によって開発した転位動力学法に基づく破壊靭性値評価法は、定量的に実験結果と一致する結果を得ることができることが示された。
3.2−15

図 3.2.1-19 温度条件を変化させた転位によるき裂先端の遮蔽効果の計算結果

図 3.2.1-20 破壊靭性値の計算結果と Tanaka らによる実験結果(2)放出される転位のすべり面の影響を明らかにするために、平成 29 年度に開発した斜め方向のすべり面への転位の放出モデルを使い、一つの転位が与えるき裂先端への遮蔽の効果を調べた。図 3.2.1-21 に結果を示す。転位が囲むき裂先端部（x/l=0.5）においては、
応力拡大係数が低下していることがわかる。これは、転位による遮蔽の効果があることを意味している。一方、その両側において、応力拡大係数が上昇していることがわかった。
これは、転位が存在することでき裂先端の応力拡大係数が上昇することがあることを意味している。この結果から、これまでの転位による遮蔽効果の逆の効果（逆遮蔽）が存在す
ることがわかった。

3.2−16

図 3.2.1-21 斜めのすべり面上の転位による遮蔽の効果(6) まとめ
平成 29 年度までに、転位を用いたき裂の数値解析法の貫通き裂問題への拡張を実現することに成功した。さらに計算精度について検討することができた。その結果、計算精度を向上させる方法を開発し、その精度について検討することができた。転位の放出モデルの開発と実装をおこなった。まず、き裂先端を含むすべり面に転位を放出するモデルを開発した。開発した転位の放出モデルを実装した転位動力学コードを用いて、き裂先端からの転位の放出、運動、弾性場を考慮したき裂の数値解析を実施し、転位によるき裂先端への遮蔽効果を計算可能であることを示した。また、転位の易動度を変化させることによって、温度の変化の影響を擬似的に考慮したシミュレーションを行った結果、低温側で破壊靭性値が低く計算されることを確認した。平成 30 年度では、き裂先端から任意のすべり面への転位の放出モデルを開発および転位動力学コードに実装することによって、任意のすべり面への転位の放出の取り扱いが可能になった。本モデルを実装した転位動力学コード
を用いた数値計算によって、すべり面が異なることによって、き裂先端から放出された転位の形状や挙動が異なることを確認した。さらに発生する転位の方位によっては、逆遮蔽が起こることがわかった。
参考文献(1) A. Takahashi, N.M. Ghoniem, “Fracture mechanics of propagating 3-D fatigue
cracks with paratmetric dislocations”, Philosophical Magazine A, Vol. 93, pp.
2662-2679 (2013).
(2) M. Tanaka, E. Tarleton, S.G. Roberts, “The brittle-ductile transition insingle-crystal iron”, Acta Materialia, Vol. 56, pp. 5123-5129 (2008).
(3) P. A. Gordon, T. Neeraj and M. J. Luton, “Atomistic simulation of dislocation
nucleation barriers from cracktips in α-Fe”, Modelling and Simulation ofMaterials Science and Engineering, Vol. 16, p. 045006 (2008).

3.2−17

3.2.2 照射や熱時効による破壊靭性値への影響評価（再委託先：電力中央研究所）(1) 欠陥クラスターのモデリング（H29）
原子炉圧力容器鋼における主要な脆化要因のひとつに照射により溶質原子が集積したナノメートルスケールの溶質原子クラスターが形成されることが挙げられ、溶質原子クラス
ターは転位運動の障害物となり硬化型の脆化を引き起こすと考えられている。特に Cu を多く含む鋼材では Cu が濃化したクラスターが形成されることが知られている。そこで本研究では圧力容器鋼における溶質原子クラスターの形成による脆化を念頭に、まずは FeCu 二元系における Cu クラスター（析出物）と転位の相互作用について検討を行うこととした。
転位動力学法で用いられる転位と析出物の相互の計算では、析出物の転位運動に対する
障害物強度をモデル化する必要がある。障害物強度を記述するモデルとしては、転位の張出角度に基づくモデルがよく知られている(1)。このモデルでは障害物強度 K は点欠陥との相互作用により張り出した転位の線張力との釣り合いにより計算され、転位が離脱する際の張出角度φをパラメータとして表される。
K=Gb2cos(φ/2)
ここでは G は剛性率、b はバーガースベクトルの大きさである。本モデルはシンプルである一方、張り出した転位同士の相互作用が生じるため障害物近傍での転位線張力の見積もりが難しく、また障害物近傍の転位の張り出し角度を正確に定めることが難しいこと、障害物強度の上限が張出角度 0°となることといった問題がある。また、BCC Fe 中の Cu ク
ラスターの障害物強度を記述するモデルとしては Russel らにより提案されたモデル（Russel-Brown モデル）がよく知られている(2)。Russel-Brown モデルでは、Cu クラス
ターの障害物強度はマトリックスの BCC Fe とクラスターである BCC Cu の剛性率の比として表される。従って、Russel-Brown モデルにおいてはマトリックスとクラスターの剛性率が得られればよい。Fe-Cu 2 元系合金、BCC Cu の剛性率については Liu et al.(3)により第一原理計算を用いて求められており、また、熊谷ら(4)は実機圧力容器鋼に見られる Cuの他に Ni や Si や Mn が含まれるクラスターを念頭に、Fe-Cu-X（X=Ni、 Si、 Mn）3
元系合金の剛性率を第一原理計算により求めている。しかしながら、Russel-Brown モデ
ルでは析出物の剛性率がマトリックスよりも低い場合に成り立つモデルであるものの、熊谷らの計算によればクラスターの化学組成によってはクラスターの剛性率がマトリックス
を上回ることがあり、Fe-Cu 合金以外の利用に適さない可能性がある。Monnet らは障害物の強度を原子シミュレーションから直接評価する手法を提案している(5)。Cu クラスターに対しては、過去の原子シミュレーションの結果からその障害物強度の値が求められている。また、Cu クラスターとの相互作用のモデル化には、転位芯に積層欠陥の存在を考慮した Peierls-Nabaro モデルを拡張した小転位モデル(6)を用いる方法もある。小転位モデルでは一本の転位をいくつかの小転位に分割して転位芯構造を記述し、マトリックスと Cu ク
ラスター中の転位芯構造の違いにより転位運動に対するクラスターの影響を算出することが可能である。このとき小転位間における力の釣り合いは、弾性エネルギによるものと積層欠陥エネルギによる力により記述される。したがって小転位モデルでは、積層欠陥エネ
ルギを求めることにより障害物強度を評価可能となる。Fe-Cu 合金においては、これまで
3.2−18

に{110}上での一般化積層欠陥エネルギが第一原理計算により計算されているが、他のすべり面に関する結果は見当たらない。そこで、以降では Cu クラスターの一般化積層欠陥エ
ネルギについての分子動力学計算を行うこととした。
積層欠陥エネルギおよび格子定数と弾性定数を計算した。Cu 濃度 50at%、 75at%の FeCu 合金について計算を行った。弾性定数、積層欠陥エネルギについては BCC 構造の純 Cuについても計算を行った。
図 3.2.2-1 に Cu 濃度が 50at%および 75at%の Fe-Cu 合金における格子定数の温度依存性を示す。Cu 濃度が低い場合（図中に併せて示す）で見られた 40K 付近の格子定数の極小値はみられなくなり、格子定数は温度に伴い単調増加する。図 3.2.2-2 には 100K および 300K における格子定数の Cu 濃度依存性を示した。いずれの温度においても格子定数
は Cu 濃度に対して単調に増加することがわかる。図 3.2.2-3 には弾性定数の Cu 濃度依存性を示す。Cu 濃度が増加するともにいずれの定数も低下することがわかった。図 3.2.2-4～図 3.2.2-6 には Cu 濃度が 50 at%、 75 at%、 100 at%（BCC 純 Cu）の{110}、 {112}、{123}上の一般化積層欠陥エネルギを示す。参考に純 Fe の結果も示した。いずれの面上でも一般化積層欠陥エネルギは 0.5b の位置に置いて単一ピークを示し、ピーク値すなわち不安定積層欠陥エネルギは BCC 純 Cu で 3.58×10-1 J/m2 であった。これは大沼らによる第一原理計算(7)により求められた BCC 純 Cu の不安定積層欠陥エネルギである 3.31×10-1
J/m2 と近い値であり、本分子動力学計算においても概ね妥当な値が求められているものと思われる。また Cu 濃度 50at%、 75%では純 Fe と純 Cu の間の値を示し、Cu 濃度が高くなるほど一般化積層欠陥エネルギが低下し、純 Fe との差が大きくなることが示された。

2.94

Cu濃度

2.93

0.0%
0.25%
0.5%
1.0%
3.0%
5.0%
50.0%
75.0%

2.92

格子定数（A）

2.91
2.9
2.89
2.88
2.87

2.86
2.85
0

200

400

600

800

1000

温度（K）

図 3.2.2-1 Fe-Cu 合金中の格子定数の温度依存性

3.2−19

2.890
100K

2.885

300K

2.880

格子定数（A）

2.875
2.870
2.865
2.860
2.855
2.850
0

20

40

60

80

Cu濃度（at%）

図 3.2.2-2 Fe-Cu 合金中の格子定数の Cu 濃度依存性

300

K[GPa]
C11[GPa]
C44[GPa]250

G[GPa]
C12[GPa]

弾性率[GPa]

200

150

100

50

0
0

20

40

60
Cu濃度[%]

80100

図 3.2.2-3 Fe-Cu 合金中の弾性定数の Cu 濃度依存性

0.7

純Fe
Cu=50%

エネルギー[J/m2]

0.6

Cu=75%
純Cu

0.5
0.4
0.3
0.2
0.1

0.0
0.0

0.2

0.4

0.6

0.8

1.0

変位量 x/b

図 3.2.2-4 Fe-Cu 合金の一般化積層欠陥エネルギ (110 面)

3.2−20

0.9

純Fe

0.8

Cu=50%
Cu=75%

エネルギー[J/m2]

0.7

純Cu

0.6
0.5
0.4
0.3
0.2
0.1
0.0
0.0

0.2

0.4

0.6

0.8

1.0

変位量 x/b

図 3.2.2-5 Fe-Cu 合金の一般化積層欠陥エネルギ (112 面)

0.9

純Fe

0.8

Cu=50%
Cu=75%

エネルギー[J/m2]

0.7

純Cu

0.6
0.5
0.4
0.3
0.2
0.1
0.0
0.0

0.2

0.4

0.6

0.8

1.0

変位量 x/b

図 3.2.2-6 Fe-Cu 合金の一般化積層欠陥エネルギ (123 面)(2) 照射や熱時効による破壊靭性値への影響評価（H30）
中性子照射された圧力容器鋼では、ナノメートルスケールの溶質原子クラスターや転位
ループ等の照射欠陥が最大で 1023m-3 のオーダーの数密度（数十 nm の欠陥クラスターの間隔に相当）で形成されることが知られている。また圧力容器鋼を模した Fe モデル合金の熱時効材においても、同様にナノメートルサイズの Cu クラスターが形成される。したがって、転位動力学法によるき裂・転位の相互作用をこれらの材料への適用する場合には、上記の空間スケールにおける微細な照射欠陥クラスターを転位動力学解析の中で取り扱える必要がある。本研究ではこれまでに、マイクロメートルスケールのモデル結晶を対象に純鉄におけるき裂・転位の相互作用の計算手法を開発すると共に、数マイクロメートルの欠陥クラスターと転位の相互作用のモデル化を進めてきた。これにナノメートルオーダーの欠陥を導入する場合、以下の点に留意する必要がある。
・結晶中の欠陥クラスター間隔が十ナノメートルのオーダーになると、これまでの転位動
力学解析で設定される転位セグメントと同等の長さとなり正確な計算ができない。
・転位と欠陥クラスター（Cu クラスター）の相互作用を考慮する際に必要な、転位セグメ
ント上における格子復元力の評価点を適切に設定できない。
3.2−21

これらを回避するには、計算における転位セグメントの最小長さを数 nm まで小さくする必要があるが、これにより計算セグメント数が増大し、計算コストが非常に大きくなってしまう。そこで、今回は、転位セグメント長を 5nm と小さく設定すると共に、結晶サイズも同時に小さくして、計算コストとのバランスをとることとした。これに伴い、き裂先端で設定される転位の初期長さと関連の計算パラメータを併せて再設定した。これにより、
実際の圧力容器鋼中性子照射材や鉄モデル合金熱時効材中に形成される欠陥クラスターと同等のスケールの計算が可能となった。
上記のモデルを用いて、純鉄中の 110 すべり面上に 3nm～9nm の Cu クラスターを、17nm～50nm の間隔で正方格子状に配置した系モデルを作成した。モデルの概要を図 3.2.2-7 に示す。計算ボックスのサイズは l=1100nm, t=50nm, h=200nm, 2a=1000nm とし、き裂面は 001
面上に配置した。き裂・転位の相互作用を行い、き裂先端における応力拡大係数を見積もった。図 3.2.2-8 はき裂・転位相互作用の転位動力学解析結果の例で、き裂先端から放出される転位がすべり面上の Cu クラスターと相互作用して運動している様子が認められる。図3.2.2-9 は直径 3nm の Cu クラスターをクラスターの間隔 17nm、25nm、50nm として解析した結果、図 3.2.2-10 はクラスター間隔を 17nm に固定し、クラスターサイズ 1.5nm、3nm,
4.5nm として解析した結果である。横軸は転位動力学計算に基づき負荷応力から計算される応力拡大係数、縦軸はき裂先端における転位の応力遮蔽効果を考慮して算出された、き裂先端の応力拡大係数である。Cu クラスターの間隔が狭いほど、またクラスターサイズが大きいほど、負荷応力から計算される応力拡大係数が同じでも、き裂先端における応力拡
大係数が相対的に大きくなっていることがわかる。応力拡大係数の絶対値について論じるにはさらなる検討が必要なものの、この結果は、圧力容器鋼が照射欠陥クラスターの形成より脆化することと定性的に一致しており、金属組織変化と脆化を結ぶメカニズムの一端が、き裂先端における転位による応力遮蔽効果によって説明できる可能性を示している。

Cuクラスター
き裂面
h
2a
すべり面
t
l

図 3.2.2-7 転位動力学解析モデルの概要

3.2−22

Cu クラスター

き裂転位

格子転位

すべり面
き裂先端

き裂面

10n m

図 3.2.2-8 き裂・転位相互作用の転位動力学解析結果の例

き裂先端における応力拡大係数, Ktip (MPa√m)

0.8
d = 17n m

0.7

d = 25n m
d = 50n m

0.6

Cu クラスター無し

0.5
0.4
0.3
0.2
0.1
Cu クラスターサイズ 3n m

0
0

0.1

0.2

0.3

0.4

0.5

0.6

負荷応力から計算される応力拡大係数, Kapp (MPa√m)

図 3.2.2-9 負荷応力から計算される応力拡大係数とき裂先端の応力拡大係数の比較（クラスター間隔の影響）

3.2−23

き裂先端における応力拡大係数, Ktip (MPa√m)

0.8
d = 4.5n m

0.7

d = 3.0n m
d = 1.5n m

0.6

Cu クラスター無し

0.5
0.4
0.3
0.2
0.1
Cu クラスター間隔 17n m

0
0

0.1

0.2

0.3

0.4

0.5

0.6

負荷応力から計算される応力拡大係数, Kapp (MPa√m)

図 3.2.2-10 負荷応力から計算される応力拡大係数とき裂先端の応力拡大係数の比較（クラスターサイズの影響）(3) 専門家ヒアリング（H29、H30）
本事業の成果を実機の脆性破壊挙動の予測に適用する場合にどのような課題があるかを検討・明確化するために、脆性破壊挙動の実験面からの評価を行っている専門家である
一般財団法人電力中央研究所材料科学研究所の三浦直樹博士との意見交換を 2017 年 10 月
12 日に、同研究所所長の曽根田直樹博士との意見交換を 2019 年 1 月 25 日に行った。両氏に対し、本事業の概要とこれまでの研究成果をご説明した後に質疑応答を行った。三浦氏からは「将来が楽しみな研究でブレークスルーになり得る」との評価を受けた。また、曽根田氏からは、技術的に詳細な議論とともに、「非常に良い成果が上がっているので、論文投稿を推奨したい。」「本技術により破壊靭性と引張特性や硬さ等、他の機械的特性との関係を理論的に明らかにすることで、本技術の実用性が高まる」との助言を受けた。
質疑応答：

三浦直樹博士との質疑応答（Q：質問、A：回答）
Q「き裂先端における真の破壊靭性値 Ktip の臨界値をどう決めるのか？」A「グリフィスのモデルから計算している。」
Q「転位のつり合い状態は求まるものなのか？」A「説明資料にて示した手順で求めている。」
Q「Ktip と見かけの破壊靭性値 Kd の割合はどの程度か？」A「転位のモビリティによって異なる。」
3.2−24

Q「き裂先端は小規模降伏しているわけだが、これは考慮に入ったモデルなのか？」A「転位が動くことが塑性変形なので、考慮に入っている。」
Q「遷移温度域の破壊靱性はばらつくが、これを再現することは可能か？」A「き裂面の結晶方位が関係している可能性があるため、様々な結晶方位のき裂面を試す。」
Q「計算でモビリティを変化させているが、古典分子動力学法の計算によれば、モビリテ
ィの温度依存性はないのではないか？」A「らせん転位のみ温度依存性が見られる。」

曽根田直樹博士との質疑応答（Q：質問、A：回答、C：コメント）
Q「き裂近傍応力場だけでなく、どのような応力場でも転位で表現できるのか？」A「応力がゼロとなるよう転位のつり合い位置が決まるので、基本的には表現できると思われる。」
Q「Kapp（負荷応力から計算される応力拡大係数）と Ktip（き裂先端における応力拡大係数）のグラフにおいて、横軸 0.7 近傍では、温度ごとの変化が見られないとのことだが、モビリティを変化させる計算でもその近傍ではそれほど違いがないので、横軸の値がもっと高ければ、温度ごとの変化がみられるのではないか。」
A「転位がき裂先端から発生すると違いが生じるが、この計算の場合では横軸 0.7 でも転位が発生しており、温度ごとの違いが出るはずである。」
Q「どのような破壊形態を対象としているのか。脆性であれば、多数の発生点から転位が発生するモデルとは一致していないと思われる。」
A「いわゆるき裂先端の小規模降伏に対応しているのではないかと想像している。微視的なブランティング等の影響（の有無）については将来の検討課題。」
Q「DD による破壊靭性計算において、析出物サイズが大きいのではないか。」
A「指摘のとおり、圧力容器鋼で想定されるよりも大きな析出物を導入している。計算の安定性等を考慮したためである。き裂のサイズと析出物サイズの比が問題になる。現在は、
小さな析出物の効果を導入するために、まずき裂の無い系で析出物によるモビリティ変化
を DD 計算し、それをき裂の DD 計算に持ち込むことを検討している。」
Q「き裂と析出物のサイズ比が問題ならば、き裂形状を半無限き裂にしてはどうか。」
A「半無限き裂にすると負荷した破壊靭性を計算することができなくなる。」
C「定量的なデータを出すことということに対しては、計算パラメータ（モビリティ等）
で記述された横軸の相当する物理量（温度など）に変換すれば、実験と比較できるデータになる。」
C「破壊靭性値と硬さなどの別の物性の理論的な関係を明らかにすると実用性が高まる。
今後のインフラの劣化評価等では、実機から微小サンプルを採取した評価も考えられることから、微小サンプルの物性から破壊靭性を得るための根拠を与える研究に展開できれば非常に有用である。」
C「非常に良い成果が上がっているので、論文投稿を推奨したい。」
Q「混合転位はらせん転位と刃状転位の組み合わせで表現できないのか。」
3.2−25

A「弾性場は表現可能であるが、転位芯や転位線の形状やモビリティは難しいことが本研究により分かった。」
Q「古典 MD で転位芯構造を表現することができるのか」A「どの程度詳細にかつ詳しく表現できるかは不明。転位線のおおまかな形状を観察する程度であれば十分ではないかと考えている。」
Q「純鉄の破壊靭性値の実験値とは比較しているのか？」A「今のところ比較は行っていない。実験値がないのではないか。」
Q「結晶粒界の影響はどう考慮するのか。1 結晶粒だけを考えてもよい。」
A「本計算では要に考慮していないが、パイルアップの影響が強く効く可能性がある」(4) まとめ
平成 29 年度では、転位運動の障害物として Cu クラスターを想定し、転位との相互作用計算に必要なパラメータについて調べた。ここでは小転位モデルによる相互作用計算に必要な Cu クラスターの一般化積層欠陥エネルギを算出し、Cu クラスター中の Cu 濃度が増大するほど一般化積層欠陥エネルギが低下し、純 Fe との差が大きくなることがわかった。
また、平成 30 年度では、転位動力学解析モデルに圧力容器鋼中性子照射材や Fe モデル合
金熱時効材中に形成される欠陥クラスターを想定した、Cu クラスターを数十ナノメートル
間隔で導入した。導入にあたっては、ナノメートルスケールのサイズ、間隔の Cu クラス
ターの影響を適切に計算できるよう、計算ボックスサイズ、最小転位セグメントサイズ等の計算パラメータの再調整を行った。この解析モデルを用いて、純鉄中のすべり面上に Cu
クラスターを配置した系でのき裂・転位相互作用の計算を行った結果、負荷応力から計算されるみかけの応力拡大係数が同じでも、Cu クラスターの間隔が狭いほど、またクラスタ
ーサイズが大きいほど、すなわち欠陥クラスターにより転位が動きにくくなるほど、き裂先端における応力拡大係数が相対的に大きくなるとの結果が得られた。応力拡大係数の絶対値について論じるにはさらなる検討が必要なものの、この結果は、圧力容器鋼が照射欠
陥クラスターの形成より脆化することと定性的に一致しており、金属組織変化と脆化を結ぶメカニズムの一端が、き裂先端における転位による応力遮蔽効果によって説明できる可能性を示している。
さらに、破壊靭性に関する実験研究に携わる一般財団法人電力中央研究所の三浦直樹博士、曽根田直樹博士のヒアリングを行い、技術的議論を行うと共に「将来が楽しみな研究でブレークスルーになり得る」「本技術により破壊靭性と引張特性や硬さ等、他の機械的特性との関係を理論的に明らかにすることで、本技術の実用性が高まる」とのコメントを得た。
参考文献(1) 例えば，Bacon and Hull, Introduction to dislocations, Butterworth-Heinemann,(2011)(2) K.C.Russell, L. M. Brown, "A dispersion strengthening model based on differing
elastic moduli applied to the iron-copper system", Acta Metallurgica. Vol. 20,
pp.969-974 (1972).
3.2−26

(3) J. Z. Liu, A. Van De Walle, G. Ghosh, M. Asta, Structure, energetics, andmechanical stability of Fe-Cu bcc alloys from first-principles calculations.
Physical Review B, Vol. 72, 144109 (2005).
(4) 熊谷, 野本, クラスターの硬化に及ぼす化学組成の影響－第一原理計算による 3 元系
Fe-Cu-X(X=Si/Ni/Mn) 合金の剛性率－電力中央研究所研究報告 Q16009(2017)(5) G. Monnet, in: I. Ionescu, S. Bouvier (Eds.), Plasticity of Crystalline
Materials: From Dislocations to Continuum, ISTE Ltd., pp. 3–36 (2011).
(6) Banerjee, N. Ghoniem, G. Lu, N. Kioussis, Non-singular descriptions of
dislocation cores: a hybrid ab initio continuum approach, Philosophical
Magazine. Vol. 87, pp.4131-4150 (2007).
(7) 大沼，中島，別役，曽根田，"一般化積層欠陥エネルギを用いた溶質元素の転位芯構造の計算 -体心立方構造の鉄，銅，ニッケル，マンガン，Fe3Si の転位芯構造-"，電力中央研究所研究報告 Q10035 (2011)

3.2−27

3.3 き裂先端における力学的状態と破壊靭性値の関係の評価（H28〜H30）
3.3.1 大規模計算結果の可視化（再委託先：東洋大学）
き裂先端における力学的状態と破壊靭性値の関係を評価するためには、き裂先端における転位の遮蔽効果を視覚的に理解し分析する必要がある。ここで、転位動力学と有限要素法とを結合させ、大規模な破壊靭性値の計算を実施した場合においては、データ量の多くは有限要素法のデータであり、これまでに有限要素法のスパコン上での解析で培われてきた財産を応用することが可能である。
転位動力学法の可視化に対する要求として、転位の形状を表すだけではなく、転位が作る応力分布を可視化する必要がある。転位動力学法では、計算過程にいて転位線上の応力のみが計算されるため、空間に広がる応力分布の情報は通常計算されない。このような転位が作る応力分布の可視化について、2 つの方法を検討した。
１つ目の方法は、可視化する応力分布の面や空間を固定し、面や空間を格子に分割する。
転位動力学法の計算過程において、その面や空間の格子点の応力を通常の転位動力学法における応力の計算方法で計算し、出力する方法である。可視化ソフトは、出力された応力を表示するのみとなる。
一方、２つ目の方法は、転位動力学法から出力された転位の情報を可視化ソフトが読み込み、可視化ソフト内で可視化の対象となる面や空間における応力を計算し可視化する方法である。この場合、可視化ソフトにおける計算量が問題となり、ディスプレイに応力を可視化するまでの時間が長くなってしまう。この計算時間の対策として、解析対象全体を二次元あるいは三次元直交格子に分割し、転位線を格子の稜線で近似することで計算の高速化を行う方法を検討した。さらに可視化処理を計算の時間方向に（タイムステップで）分割し、各タイムステップの情報を並列計算環境下で分散させることにより、各タイムス
テップの応力の計算を並列に計算することで高速化を行う方法の検討も行なった。
また、スパコン上での大規模計算に対し、必要な部分の可視化処理をスパコン上で行い、ローカルホストへ結果を転送し、ローカルホスト上で可視化を実現するシステムについて、ミクロからマクロまでの解析を含んだマルチスケールへの対応の準備を行った。マルチス
ケールでは、可視化に要求されるデータ量、可視化範囲のスケールも異なることから、効率的なデータの選択および圧縮を行った後に転送を行わないと、ローカルホストでの可視化が実用的な処理時間での実現が極めて困難となることから、解析規模、データ転送量、可視化精度の関係を見積もり、求められる精度に対する、現実的な時間での可視化処理を選択可能にするシステムの構築を検討した。
平成 29 年度において、上記の検討をもとに、スパコン上での大規模計算に対し、必要な部分の可視化処理をスパコン上で行い。ローカルホストへ結果を転送し、ローカルホスト
上で可視化を実現するシステムについて、ミクロからマクロまでの解析を含んだマルチス
ケールへの対応を行った。マルチスケールでは、可視化に要求されるデータ量、可視化範
囲のスケールも異なることから、解析規模、データ転送量、可視化精度の関係を見積もる
ことで、求められる精度に対する、現実的な時間での可視化処理を選択可能にするシステ
ムの構築を行った。
また、転位動力学法において転位が作る応力分布の可視化手法について検討を行った。
3.3−1

可視化する応力分布の面や空間を固定し、面や空間を格子に分割する。転位動力学法の計算過程において、その面や空間の格子点の応力を通常の転位動力学法における応力の計算
方法で計算し、出力する方法である。より具体的には、転位動力学解析結果データおよびFEM 解析結果データの双方について、これを別途３－D 直交格子フォーマットに変換を行ってマージする。直交格子データの可視化についてはすでにマクロな連続体力学シミュレー
ションの分野、特に流体解析における断面表示、等値面、ボリュームレンダリングなどで実績があり、既存の可視化技術や可視化ソフトが利用可能である。
可視化の例として、Frank Read ソースの転位動力学シミュレーション結果に対する応力の等値面表示を図 3.3.1-1 に、断面コンター表示を図 3.3.1-2 に示す。

図 3.3.1-1 応力分布：ｘｘ応力成分の等値面表示

図 3.3.1-2 応力分布：相当応力の断面コンター表示
3.3−2

またこれに続き平成 30 年度では、転位線群およびそれらが作る応力分布との重ね合わせができるようなシステムを構築した。転位動力学法の可視化では、転位の形状を表すことに加え、これに転位が作る応力分布をオーバーラップ形式で重ね合わせ表示する必要がある。転位動力学法では、計算過程において転位線上の応力のみが計算されるため、空間に広がる応力分布の情報は通常計算されない。このような転位が作る応力分布の可視化、および転位線との重ね合わせが行えるようなシステムの試作および実装を行い、転位線セグ
メント群の表示と応力の断面表示および等値面表示機能などの各種機能を実装した。これにより、転位動力学シミュレーション結果の分析をより詳細に行うことが可能となった。
図 3.3.1-3 と図 3.3.1-4 は、転位動力学シミュレーション向けの可視化システムの動作画面であり、それぞれ転位線の可視化と応力分布・転位線の同時可視化の結果を示している。

図 3.3.1-3 転位線の可視化

図 3.3.1-4 応力分布と対応する転位線の同時表示

3.3−3

3.3.2 超大規模解析に対して適している計算アルゴリズムの検討とコード整備（再委託先：
東洋大学）
東洋大学は平成 23 年から 29 年にかけて、JST CREST「ポストペタスケール高性能計算に資するシステムソフトウェア技術の創出」領域プロジェクトの「ポストペタスケールシミ
ュレーションのための階層分割型数値解法ライブラリ開発」チームにおける中核研究センターとして、最先端スーパーコンピュータ上での大規模構造解析シミュレーションの豊富な実績を有し、この技術をベースとして、平成 28 年度には転位動力学法(DD)と有限要素法(FEM)を融合させた有限体中のき裂と転位の数値解析の超大規模化のためのコード整備を行った。具体的には、「き裂」「転位動力学」「破壊靭性値」をキーワードとして、先行研究の調査を進め、超大規模解析に最適な解析アルゴリズムについて検討し、超大規模線形弾性解析コードと転位動力学法のハイブリッドの実装を行うための準備を進めた。
また、外部ライブラリ利用のための、統一データフォーマットの導入を行った。スパコ
ン上で高い計算性能を示している開発ライブラリでは、入出力フォーマットとして ADV/IO
をベースとしたシステムを導入しており、本研究で採用した。これまでに開発されてきた様々なコードについて、データフォーマットを統一することにより、外部ライブラリの利用および、コード間の連携を容易にし、より高速かつ効率的な計算システムを構築する準備を行った。これにより例えば、マクロ―ミクロ連成のための初歩的検討を行うため、図3.3.2-1 に示されるようなき裂入り試験片のマクロ解析結果から、き裂先端におけるミク
ロ解析へ接続することを可能とした。

図 3.3.2-1 き裂入り試験片のマクロ解析とき裂先端でのミクロ解析
平成 29 年度には超大規模解析に適した計算アルゴリズムの検討とコード整備を行い、転位動力学法と有限要素法による連成解析システムを開発した。図 3.3.2-2 に開発システム
における連成解析時のデータフローの様子を示す。
開発システムの計算手順は以下の通りである。
1. 転位動力学法(DD)による計算領域と同じサイズの直方体 3 次元 4 面体要素メッシュを用
3.3−4

意
2. 表面三角形の重心点座標を抽出
3. 重心点座標を DD に入力し、応力値を計算
4. 重心点座標における応力を表面力[F]に変換
5. 表面力[F]を境界条件、直方体底面を完全固定の FEM 用入力ファイルを作成
6. FEM による弾性解析を実行
本システムを用いて、2,400 万要素数、3,200 万節点数、9,600 万自由度数の有限要素モ
デルを用いた解析に成功し、DD-FEM 連成解析による様々な転位の放出に対する応力拡大係数の計算が可能となった。

図 3.3.2-2 DD-FEM 連成解析のデータフローの様子
平成 30 年度には重ね合わせの原理を用いた有限体を取り扱うことのできる転位動力学を用いた数値解析を実施し、破壊靭性値に対する与えた応力条件の影響を調べた。転位動力学法と有限要素法を融合させた有限体中のき裂と転位の数値解析の超大規模化のためのコード整備を進め、超大規模線形弾性解析コードと転位動力学法の連成解析の実装を行い、有限体中の転位を含むき裂解析を、き裂モデル転位動力学・有限要素法連成解析により実現した。これまで行ってきた単純引っ張り応力の境界条件に加え、せん断引っ張り、曲げ
引っ張り応力などの異なる応力条件によるシミュレーションにも成功し、破壊靭性値との関係を整理することができた。図 3.3.2-3 にき裂の様子を示す．また、従来から開発を行ってきた手法の適用により、超大規模解析による性能低下は生じなかったため、これに対するチューニング作業は発生しなかった。

3.3−5

図 3.3.2-3 連成解析によるき裂の様子

3.3−6

3.4 研究推進（H28～H30）
研究代表者の下で各研究項目間における連携を密にして研究を進めるとともに、研究実
施計画等を協議するために、研究委員会を開催した。
平成 28 年度は、研究委員会を 3 回開催した。第１回目の委員会（平成 28 年 12 月 7 日開催）では、キックオフミーティングとして、研究内容の再確認と研究の進め方及びスケジ
ュールの確認が行われた。第２回（平成 29 年 1 月 25 日開催）、第３回（平成 29 年 3 月 21日開催）の委員会では、研究内容及び進捗に関して活発な意見が交換され、研究を遂行す
るにあたり有意義な議論の場となった。また、研究の方向性について参加者全体で確認した。
平成 29 年度は、研究委員会を 6 回開催した。第１回目の委員会（平成 29 年 5 月 17 日開催）では、年度のキックオフミーティングとして、研究内容の再確認と研究の進め方及びスケジュールの確認が行われた。第２回（平成 29 年 7 月 12 日開催）、第３回（平成 29 年
9 月 15 日開催）、第４回（平成 29 年 11 月 24 日）、第５回（平成 30 年 1 月 26 日）、第６回（平成 30 年 3 月 14 日）の委員会では、研究内容及び進捗に関して活発な意見が交換され、研究を遂行するにあたり有意義な議論の場となった。また、研究の方向性について参加者全体で確認した。中でも、第３回は近畿大学、第４回は諏訪東京理科大学、第５回は九州大学で開催し、本研究に携わっていない研究者と議論することができた。
平成 30 年度は、研究委員会を 6 回開催した。第 1 回目の委員会（平成 30 年 5 月 16 日開催）では、過年度の研究成果の確認と、最終年度に向けた研究課題の確認を行った。第 2回目（平成 30 年 7 月 11 日開催）、第３回目（平成 30 年 9 月 14 日開催）、第４回目（平成 30 年 11 月 22 日開催）、第５回目（平成 31 年 1 月 25 日開催）、第６回目（平成 31 年
3 月 12 日開催）の委員会では、研究内容及び進捗に関して活発な意見が交換され、研究を遂行するにあたり有意義な議論の場となった。特に最終回となった第６回では、本研究課
題の今後の展開について議論を行った。また、第３回は八戸工業大学、第４回は徳島大学、
第５回は琉球大学で実施し、関連分野の研究者との議論を通じ、本研究課題の方向性についての再認識を行った。
さらに、国際会議での研究発表 9 件（The 9th International Conference on Multiscale
Materials Modeling での 2 件を含む）、国内会議での研究発表 4 件を行い、国内会議での研究発表を 1 件予定している。また、国際学術論文誌に論文が 2 編掲載され、2 編が現在投稿中である。

3.4-1

4．結言
4.1 研究の到達点
混合転位の運動の性質の温度依存性の解明とモデル化においては、パイエルス応力と転位
芯構造の関係を明らかにし、特にいくつかの混合転位においてパイエルス応力が高くなることが明らかとなった。また、転位芯に関する分析から、転位が一度に移動するための応力、またはキンクが転位に沿って移動するための応力であるかについて、パイエルス応力の意味が異なることが明らかになった。らせん転位のみではなく、いくつかの混合転位においてもパイエルス応力が高くなることが明らかになった。NEB 法を用いたキンク対の形成に必要な活性化エネルギを計算した結果、実験結果と定量的に成功した活性化エネルギを計算することが可能であることを示し、活性化エネルギのせん断応力依存性を調べた。その結果、せん断
応力を大きくすると、活性化エネルギが低くなることがわかった。らせん転位以外のパイエ
ルス応力が高い混合転位についても活性化エネルギを調べた所、転位線の方向の近接距離が長くなるほど、活性化エネルギが低くなることがわかった。さらに転位の運動特性について{112}と{123}面上の転位の易動度を詳細に調べ、その結果を、サイン関数を用いた式でフィッティングすることによって、様々な転位の易動度が表現可能であることを示した。この結果より、フィッティングした式を転位動力学コードに実装すれば、転位のすべり面および結晶方位依存性を転位動力学シミュレーションによって考慮することが可能になる。鉄中に銅原子が固溶した場合の、弾性定数、格子定数、積層欠陥エネルギを、銅濃度を変化せて調べた。さらに、転位の易動度と摩擦応力について調べた結果、銅の含有量は、転位の易動度に対する影響は小さく、摩擦応力に対して大きく影響することがわかった。また、鉄を母材とした合金に関して、既存の研究の調査を行なった結果、本研究の将来的な対象として考え得る Fe-Cu や Fe-Cu-Ni について信頼性があるポテンシャルが開発されていることを確認することができた。
転位運動モデルを用いた、き裂先端での転位の発生・運動による遮蔽の効果を転位動力学に基づいて考慮することによる破壊靭性値の温度依存性の計算技術の開発においては、貫通
き裂に拡張し、応力速度を境界条件とする数値解析において高精度に応力拡大係数を算出する方法を開発することに成功した。き裂の転位動力学コードにき裂先端からの転位の放出モデルを実装することで、き裂先端からの転位の放出およびき裂先端の遮蔽効果の計算を可能にした。さらに、分子動力学法による転位の運動特性のモデリングの結果として得られた転位の易動度の温度および方位依存性モデルを実装し、転位の遮蔽に対する温度の影響をシミ
ュレートすることを可能にした。このモデルを用いて貫通き裂における転位の遮蔽効果のシ
ミュレーションを実施した結果、脆性―延性遷移温度が 140K から 160K の間であることが計算された。実験結果では脆性―延性遷移温度はおよそ 150K となっており、本計算結果は実験結
果と定量的に一致している。すなわち、本研究によって開発された計算科学基盤を用いて計算することによって、実験結果を再現する結果を得ることに成功した。三次元の転位の運動に基づく計算機シミュレーションによってこのような結果を得たのは、世界初の成果である。
さらに転位のすべり面やすべり方向によっては、逆遮蔽を起こす場合があることがわかった。
照射・熱によって金属中に形成される照射損傷を転位動力学計算で考慮することで、照射・熱による破壊靭性値の変化を計算する技術の開発においては、鉄と銅の積層欠陥エネル
4-1

ギを計算し、銅濃度が高くなるにつれ、積層欠陥エネルギが低下することがわかった。これらのパラメータを用いた転位と銅析出物の相互作用モデルを実装し、転位の遮蔽効果への銅析出物の影響について計算を実施した。その結果、析出物が大きくなるにつれ、さらに析出
物間間隔が狭くなるにつれ転位による遮蔽の効果は弱まり破壊靭性値が小さくなる傾向（脆化する傾向）を得た。これは、これまでの実験結果と定性的に整合する結果であり、析出物の分布や大きさが与える破壊靭性値への影響を計算機シミュレーションを用いて捉えた世界初の結果である。今後さらに詳細な相互作用モデルや析出物の配置を入力することによって、本研究によって開発した方法が、照射による脆化のメカニズム解明のための有力なツールとして応用可能であることを意味している。
き裂先端での転位動力学が有限要素法による変形解析を融合することで種々の力学的状態における破壊靭性値の変化を計算し、き裂先端での拘束条件の影響を評価するために、計算結果の可視化では、転位が作る応力の可視化のためのフレームワークを構築し、転位が作る
複雑な応力場の可視化を可能にした。大規模数値解析に向けた検討では、可視化におけるデータの種類や量について検討し、既存の財産を応用可能であることが確認できた。さらに、転位動力学特有の可視化について検討し、いくつかの可能性のあるアルゴリズムについて検討を行った。転位が作る応力分布の可視化、および転位線との重ね合わせが行えるようなシステムの試作および実装を行い、転位線セグメント群の表示と応力の断面表示および等値面
表示機能などの各種機能を実装した。これにより、転位動力学シミュレーション結果の分析をより詳細に行うことが可能となった。また、転位動力学法と有限要素法を重ね合わせの原理により融合させた、任意の荷重条件下での破壊靭性値の計算を目指した研究では、大規模計算において実績を有する有限要素コードと本研究で開発した転位動力学コードを融合させ、本システムを用いて、2,400 万要素数、3,200 万節点数、9,600 万自由度数の有限要素モデル
を用いた解析に成功した。これまで行ってきた単純引っ張り応力の境界条件に加え、せん断
引っ張り、曲げ引っ張り応力などの異なる応力条件によるシミュレーションにも成功し、破壊靭性値との関係を整理することができた。
4.2 今後の展望・次のステップに向けて本研究では、マルチスケール材料モデリングに立脚した原子シミュレーションを用いた転位の運動のモデリングから始まり、モデル化された転位の運動を転位動力学法に実装し、巨視的な量である破壊靭性値を計算可能にする計算科学基盤の構築に成功した。マルチスケー
ル材料モデリングの考え方は、古典的である一方、これまでにスケールを跨いだシミュレー
ションの実現はほとんど実現されておらず、本研究によって達成された破壊靭性値に関する
マルチスケール材料モデリングは、原子力分野に限らず計算科学、計算力学の分野においても大きな意義があると言える。さらに、これまでの計算機シミュレーションを用いた原子力
材料の脆化に関する研究は、照射による微視的組織の形成や、微視的組織と転位の相互作用、材料の降伏応力の上昇に関するものに留まり、延性―脆性遷移の本質的な（直接的な）物性値である破壊靭性値を調べることが可能になったことは、今後の材料研究において多大なインパクトを与えるものである。しかし、本研究課題で達成された計算科学基盤においては、
１）き裂先端からの転位の放出モデルの詳細化、２）多結晶金属中の（結晶粒界の影響を考4-2

慮した）破壊挙動、３）母材中の潜在転位の影響など、さらなる定量化に向けては多くの課題が存在する。しかし、これらの問題は、マルチスケール材料モデリングの考え方を応用することで、その素過程のモデリングと、そのモデルを実装したシミュレーションによって、一つ一つを解決することが可能である。すなわち、本研究の今後の展望の一つの方針として、
マルチスケール材料モデリングを用いたシミュレーションモデルのさらなる詳細化・緻密化が挙げられる。この展望については、今後継続的に研究を行い、本研究課題で開発した計算科学基盤をさらに発展させ、実現していく予定である。
一方、本研究課題の大きなターゲットとして、破壊靭性値の計算を上げてきたが、これに並行して、照射や熱事項による材料の硬さの変化や降伏応力の変化、さらにこれらの物性値間の関係を明らかにする必要性が考えられる。例えば、照射された材料の劣化状態の評価においては、監視試験片等を用いた試験が行われてきた。さらに近年では、ミニチュア試験片の開発およびその妥当性の検証によって、より少ない材料を用いて劣化状態の評価を行うことが進められている。現在では、破壊靭性値の変化を知るためには、破壊靭性試験を実施す
る必要があるため、破壊靭性試験のための試験片の小型化を検討する必要がある。しかし、破壊靭性値と硬さの関係、または破壊靭性値と降伏応力の関係を知ることができれば、より多角的に様々な試験方法を用いて破壊靭性値の評価を行うことが可能になる。従来、転位動
力学法は、転位の挙動に基づいた塑性変形挙動の数値解析手法として開発されてきた。すな
わち、塑性変形の開始に必要な応力である降伏応力の計算を行うことに適した数値解析手法であると言える。本研究課題では、マルチスケール材料モデリングの考え方を用いて、BCC 鉄の転位運動のモデルを転位動力学法に実装することに成功している。この転位動力学法を用いることで、これまで計算機シミュレーションを用いて数値的に解析することが困難であった BCC 鉄の塑性変形挙動および降伏応力のシミュレーションが可能になる。このような転位動
力学法による降伏応力の数値解析と、本研究課題で実現した破壊靭性値の数値解析を並行して行うことにより、降伏応力の変化と破壊靭性値の変化の関係や、これらに寄与する微視的組織の実態について研究することが可能になる。材料の硬さの変化について、これまでに研究例は見ることができないが、ナノインデンテーションのような硬さを調べる実験についての転位動力学シミュレーション研究は行われており、本研究で開発した転位動力学法をこれに応用すれば、BCC 鉄の硬さの変化についての信頼性の高い研究を展開することが期待できる。
以上のように、本研究課題で開発した転位動力学法に基礎をおく計算科学基盤は、破壊靭性値の数値解析法としての役割のみならず、降伏応力の変化や硬さの変化に関する数値解析的アプローチの信頼性向上にも資するものである。これらの異なる物性値の研究を統一的に実施することで、巨視的な物性値と微視的組織の関係のより詳細な関係を導き出し、さらに新しい試験法に向けた知見を与える発展性を秘めている。

4-3