Thermodynamic Analysis of the Formation Process of Metastable Carbides in Iron–Carbon Martensite

Abstract

In this study, the free energy of iron-carbon binary BCT martensite was calculated using the first-principles calculation and the cluster expansion and variational method. The free-energy curves of BCT martensite show the possibility of promoting the clustering of carbon atoms in the tempering process, because there is two-phase separation associated with the formation of metastable BCT-Fe₂C ordered structure. This BCT-Fe₂C structure was found to have many crystallographic similarities to η-carbide (Fe₂C). Then, the energy barrier required for the transition from the BCT-Fe₂C ordered structure to η-carbide was calculated by means of the G-SSNEB method. The obtained activation energy was sufficiently small, suggesting that η-carbide may be formed through the BCT-Fe₂C ordered structure. According to these findings, it was suggested that η-carbide in the low-temperature annealing process of BCT martensite may precipitate through a two-step process in which BCT-Fe₂C ordered structure is formed by the two-phase separating tendency in BCT martensite, and the ordered structure transitions to η-carbide over the energy barrier.

1. 緒言

鉄-炭素二元系におけるマルテンサイト変態は非常に硬い組織を造り出すことが可能であるため，材料の高強度化の手法として広く利用されている¹⁾。この強化機構には，マルテンサイト変態に伴って結晶格子内に過剰に固溶される炭素による格子ひずみの誘起，ラス状またはレンズ状の微細組織の形成，マルテンサイト変態による多量の転位の導入の三つの効果が寄与していると考えられている。

急冷直後のマルテンサイトは脆く靭性が不足しているため，焼戻し処理を施すことで強度と靭性のバランスを調整する。焼戻し処理における昇温過程では温度領域ごとに様々な準安定炭化物が生成するが，これらの炭化物の析出は，次のような四つの段階を経て進行する²⁾。すなわち第一段階は100°C程度までの低温度領域での時効であり，マルテンサイト中の固溶炭素のクラスタリングが生じ，局所的な高炭素領域が生じるとされる^3,4)。第二段階は100°Cから200°Cの温度領域で進行し，Fe₂Cのη炭化物あるいはFe_2+xCのε炭化物が析出する⁵⁾。第三段階は，200°Cから300°Cの領域での時効であり，マルテンサイト変態を起こさなかった残留オーステナイトがセメンタイト(θ-Fe₃C)とフェライトへ分解する^2,6)。第四段階はそれ以上の温度領域での変化であり，χ-Fe₅C₂炭化物⁷⁾やFe_2n+1C_n組成のθ’セメンタイト構造⁸⁾の析出が確認され，最終的には安定なセメンタイトへと遷移していく。500°C以上の温度域では炭化物はセメンタイトのみとなり，温度の上昇に従ってセメンタイトの粗大化，球状化が確認されている。

このようにさまざまな炭化物が異なる温度領域で生成するが，材料組織の制御にはこれらの析出機構を十分理解する必要がある。しかしながら，この領域に関する研究は析出物と母相との方位関係や析出物の形態，焼戻し中の強度に関する研究に集中し，炭化物生成の熱力学的な考察は十分に行われているとは言えない。そのため，準安定炭化物の生成起源については，いまだに不可解な点が多く残されている。例えば，より高炭素濃度をもつη炭化物やε炭化物などはより安定なセメンタイトよりも優先して生成するが，その熱力学的起源については明らかにされていない。

この興味深い現象に関して，Fangら⁹⁾は第一原理計算に基づいた調和振動子近似による自由エネルギー計算から，60°C以下の温度領域においては，セメンタイトよりもη炭化物が熱力学的に安定であることを示し，η炭化物の焼戻し中における析出について議論している。また，焼戻しの初期段階では炭素のクラスタリング挙動が報告されているが^3,4)，このようなマルテンサイト中のクラスタリングと準安定炭化物生成の関連についてはこれまで明確にはされていない。これは一つにはマルテンサイト変態が準安定状態で起こるために，その熱力学的性質が不明であることに原因があると考えられる。そこでこのようなクラスタリング挙動や準安定炭化物の遷移を議論するためには，まずマルテンサイト中での固溶炭素のエネルギー状態を明らかにすることが重要である。

そこで，本研究では体心正方晶(以下BCTとする)マルテンサイトの自由エネルギーを第一原理計算によって評価し，低温焼戻し中に生成するη炭化物との関係を調べた。さらにその遷移過程のエネルギー障壁を計算することで，構造遷移の可能性について熱力学的観点から考察した。

本論文の構成は以下の通りである。2章では鉄-炭素二元系マルテンサイトの構造モデルの構築と，クラスター展開・変分法およびGeneralized Solid-State Nudged Elastic Band(以下G-SSNEBとする)法の計算手法について説明する。3章では，3・1節において炭素原子の固溶サイトについての計算結果を示し，3・2節でクラスター展開・変分法による自由エネルギーの計算結果を示す。3・3節においてG-SSNEB法から求めた構造遷移にともなうエネルギー障壁について議論する。4章で本研究の結果を総括する。

2. 計算条件

2・1　マルテンサイトの構造モデルの構築

鉄-炭素二元系におけるマルテンサイト変態は，結晶構造が原子の拡散を伴わず，せん断変形によって生じる無拡散変態である。高温相であるFCC-Feを急冷するとBCC-Feに変態するが，FCC-Fe中に炭素原子が固溶している場合は，BCCへの変態に伴って固溶炭素原子による一軸方向の格子ひずみが生じ，BCT格子のマルテンサイトが生成する。一軸的な格子ひずみが誘起する機構は，BCCおよびFCCの炭素固溶サイトの対応関係から説明することができる。BCCにおける炭素が固溶する八面体空隙サイトをFig.1に示した。BCCでは等価な八面体空隙サイト位置が鉄原子一個あたり三種類存在する。Fig.1中に□で示した(1/2, 0, 0)位置をaサイト，△で示した (0, 1/2, 0)位置をbサイト，×で示した(0, 0, 1/2)位置をcサイトとして呼称する。また，(0, 0, 1/2)位置のcサイトを中心とした八面体Fe格子を破線で示した。図からZ軸方向に最隣接原子である二つの鉄原子が存在し，八面体は一軸方向に押しつぶされた形状をしていることがわかる。同様に，aサイトではX軸方向，bサイトではY軸方向に二つの最隣接の鉄原子が存在する。炭素原子がそれらのサイトに固溶すると最隣接の鉄原子が離れる方向に変位するため，Fe格子がその方向に伸張する。それぞれのサイトに同じ分率の炭素原子が固溶した場合は，同時にそれぞれの軸方向に伸張するため，軸比の変化は生じず等方的な体積膨張のみが生じる。しかし，炭素を固溶したFCCからBCCへ変態する場合は，BCC中での固溶サイトに偏りが生じ，それが軸比を発現させる。Fig.2にFCCとBCCのBain対応¹⁰⁾を示した。BCCに対応する格子を太線で示し，Fig.1と対応するa，b，cの三つのサイトを同様に図中に□，△，×で示した。この図の関係から，FCC-FeはBCC構造をZ軸に方向に伸張し，X，Y方向に収縮させた構造に対応することがわかる。そのため，a，bサイトは最隣接鉄原子との距離がさらに小さくなるが，cサイトは逆に最近接鉄原子との距離が大きくなる。このため，3・1節において実際の計算結果が示されるが，FCC-Fe中では×で示したサイトに炭素が固溶した状態が他の固溶サイトに固溶するよりも大きく安定であり，炭素原子はこのサイトのみを占有する。したがって，固溶炭素が存在するFCC-FeがBCCへと変態すると，cサイトのみに炭素が固溶した状態をとる。その結果一軸方向だけが伸張することになり，軸比c/aが1よりも大きいBCT構造となる。本研究ではこの結晶学的特徴から，c/a>1のBCT-Fe格子に対してcサイトの空孔が炭素と置き換わる構造モデルを採用した。

Fig. 1.

Octahedral interstitial solution sites in BCC structure.

Fig. 2.

Bain lattice correspondence of FCC with BCC structures.

BCC構造では等価なa，bサイトにも炭素原子が固溶するが，同様にc/aが1に近いBCT構造においてもa，bサイトに炭素原子は固溶する可能性がある。しかし，詳細は3・1節で述べるが，本研究ではa，bサイトへの炭素の固溶に要するエネルギーがcサイトに固溶する場合に比べて非常に大きかったことから，この二つのサイトへの炭素の固溶は取り扱わず，cサイトだけに固溶する構造モデルを採用した。また，マルテンサイトでは炭素固溶量の増加に伴って，軸比とともに体積も膨張する。そこで，本研究では，体積と軸比の影響を考慮してクラスター展開・変分法から自由エネルギーの計算を行った。

2・2　クラスター展開・変分法

本研究では，BCT-Fe中に固溶する炭素原子間の相互作用をクラスター展開法により評価した。クラスター展開法では，有効クラスター相互作用(Effective Cluster Interactions:以下ECIsとする)を用いることで規則構造の全エネルギーが式(1)のように，ECIsとクラスターの密度に相当する相関関数の積の総和として表される。

  

$$E^{\varphi }\left(V\right)={\displaystyle \sum }_{\alpha }^{{\alpha }_{\max }}{J}_{\alpha }\left(V\right)\cdot {\xi }_{\alpha }$$(1)

式(1)中のφは規則構造を，αはクラスターを区別する添字であり，E^φ(V)は第一原理計算から求められる体積Vにおける規則構造φの全エネルギー，ξ_αはクラスターαの相関関数，J_α(V)は体積VにおけるECIsである。エネルギーE^φは無限大のサイズのクラスターを用いて正確に表すことが可能であるが，実際の計算ではクラスターのサイズを有限サイズで打ち切る必要がある。そこで，最大のサイズとなるクラスターα_maxを定義しておき，α_max中に内包されるクラスターについて式(1)の和をとることで規則構造のエネルギーを近似する。このα_maxを最大クラスターと称し，最大クラスターに内包されるクラスター群のことをサブクラスターと呼ぶ。

一方，クラスター相関関数ξ_αは，結晶格子中の任意のサイトiにおける占有スピン演算子σ(i)を用いることで表現される関数である。鉄-炭素二元系合金において炭素原子と空孔が置換し合うモデルでは，サイトiに炭素原子が存在する場合σ(i)=0，空孔が存在する場合をσ(i)=1として占有スピン演算子を割り当てた。n体のクラスター相関関数ξ_nは結晶格子中に存在する結合の総数N(n)を用いて次式のように記述される。

  

$${\xi }_n=\frac{1}{N\left(n\right)}{\displaystyle \sum _{i_1}^N{\displaystyle \sum _{i_2}^N\dots }{\displaystyle \sum _{i_n}^N\sigma \left(i_1\right)\sigma \left(i_2\right)\dots \sigma \left(i_n\right)}}$$(2)

占有スピン演算子の配列から相関関数の組み合わせは任意の規則構造に対して一意に決まるために，クラスターの濃度と相関関数は一対一の関係にある。したがって式(1)中の未知数はECIsであるJのみであるが，構築した複数の規則構造のエネルギーを第一原理計算から計算することで，式(1)の連立方程式が得られるため，これを解くことでECIsを決定することが可能である。しかし，前述のとおりα_maxにより有限サイズのクラスターで計算を打ち切るために，式(1)の左辺のエネルギーを完全に再現するECIsを決定することはできない。そこで求めるべきECIsの数，すなわちサブクラスターの数よりも十分多くの規則構造のエネルギーを計算し，得られた複数の関係式をよく再現するECIsを最小二乗法から決定する。また，α_maxやサブクラスターの選択には任意性があるため，本研究ではleave-one-out交差検証を行い，予測誤差が小さくなるようにα_maxおよびサブクラスターを選択した¹¹⁾。

このようなクラスターの配置エントロピーを考慮した自由エネルギーは式(3)のように表される。

  

$$F\left(T;V\right)={\displaystyle \sum }_{\alpha }^{{\alpha }_{\max }}{J}_{\alpha }\left(V\right)\cdot {\xi }_{\alpha }-T{\displaystyle \sum }_{\alpha }^{{\alpha }_{\max }}{\gamma }_{\alpha }\cdot S_{\alpha }$$(3)

ここでγ_αはKikuchi-Barker係数であり，S_αはクラスターαからのエントロピーの寄与である¹²⁾。この計算式を用いて，任意の温度や組成において最小となるFを与えるV，ξの組み合わせを変分法から求めることで自由エネルギーを計算した。クラスター展開法ではそれぞれの規則構造に対して，体積を変えた10-20点程度のエネルギーを求め，ECIsの体積依存性を決定した。

一方，軸比c/aも体積と同様に自由エネルギーに影響を与える。クラスター展開・変分法では，体積Vと同様の手法でc/aも式(1)，式(3)の有効クラスター相互作用Jの変数として記述することができる。これにより自由エネルギーが最安定となるc/aを変分法により決定することが可能である。しかし，変分法を適用するには，Jがc/aに対して滑らかな関数として定義される必要があり，そのためにはc/aを細かく変えた規則構造のエネルギー計算が必要となる。体積と同様にc/aを細かく変える場合計算量が膨大となるため，本研究ではc/a=1.01，1.05，1.10の三つの条件において計算した。マルテンサイトの厳密な自由エネルギーはVおよびc/aを変数としたエネルギー曲面上の最安定点になるが，本研究ではこの三つのc/aのエネルギー曲線から，マルテンサイトの自由エネルギーの特徴を調べた。

第一原理計算は密度汎関数による計算コードVienna Ab initio Simulation Package(VASP)^13–15)を使用し，交換相関汎関数にはGGA-PBE^16,17)を用いた。平面波のカットオフエネルギーは400 eVとし，k点サンプリングはFe原子サイト1つからなる菱面体単位胞に対し17×17×17のグリッド点を用い，異なるサイズのセルに対しては逆格子空間上のk点数の密度がほぼ等しくなるように設定して計算を行った。磁気構造モデルは強磁性とし，磁気モーメントの初期値として鉄原子を+2.5 μ_B，炭素原子を0 μ_Bとして与えた。また，本研究が対象とする侵入型原子と空孔が置換する構造モデルでは，原子位置の緩和に伴う構造のひずみが大きい規則構造が多数生成される。このような大きくひずんだ構造のエネルギーは，格子点上に原子が存在することを前提として解析するクラスター展開法の精度を著しく低下させるため，本研究では原子位置の局所緩和は行わず，各原子はBCT格子上に固定した状態でエネルギーを求めた。したがって本研究では，設定した軸比および原子位置を格子点に固定した状態で等方的に格子長を変えることで体積を変えたエネルギー計算を行い，規則構造の体積依存性を求めた。

本研究では，最大結合長が第8隣接原子の6体までクラスターを最大クラスターとして，52種類のクラスターに体積依存性を有したECIsを決定した。それぞれの空孔および窒素サイト一つあたりの交差相関誤差はc/a=1.01の構造で28 meV，c/a=1.05の構造で31 meV，c/a=1.10の構造で31 meVとなった。

2・3　Generalized Solid-State Nudged Elastic Band(G-SSNEB)法

物質がある状態から別の状態に遷移する時，その遷移過程は最も低いエネルギー経路Minimum Energy Path(以下MEPとする)を通る。MEP上のエネルギーが最も高い点を鞍点と呼び，鞍点のエネルギーがその遷移過程において乗り越えなければならない活性化エネルギー障壁を与える。このMEPおよび鞍点を遷移過程の始点と終点における結晶構造から求める手法がNudged Elastic Band(以下NEBとする)法である^18–20)。この手法では初めに，始点から終点へ至るN+1個の構造イメージ(R₀，R₁，…，R_N)を作成する。ここでR₀，R_Nは始状態および終状態の構造イメージである。各構造イメージR_iは遷移経路上の位相点と見做し，それらが一連のバネから構成され連結している状態として捉える。位相点R_iに作用する力F_i^NEBを

  

$$F_i^{\text{NEB}}=F_i^{\nabla \perp }+F_i^{s\parallel }$$(4)

と定義する。ここでF_i^∇⊥はポテンシャル勾配からの力の経路に対する垂直成分であり，F_i^s||はバネからの力の遷移経路に対する平行成分である。位相点R_iにおける遷移経路の接線の単位ベクトルを${\widehat{\tau }}_i$として，それぞれの力は式(5)，式(6)で表される。

  

$$F_i^{\nabla \perp }=F_i^{\nabla }-\left(F_i^{\nabla }\cdot {\widehat{\tau }}_i\right){\widehat{\tau }}_i$$(5)

  

$$F_i^{s\parallel }=k\left(\left|R_{i+1}-R_i\right|-\left|R_i-R_{i-1}\right|\right){\widehat{\tau }}_i$$(6)

ここでF_i^∇はR_iがエネルギー曲面から受けるポテンシャル力，kは仮想的に設定したバネのバネ定数である。このF_i^NEBが全ての位相で最小となるように，各位相の構造を最適化することでMEPを探索する。ただし，このNEB法だけの計算ではエネルギーの低い位相点に収束するため，イメージの間隔を十分細かく設定しない限り，正確な鞍点のエネルギーを得ることができない。そこでClimbing-Image NEB(CI-NEB)法が用いられる。これは収束したMEPにおいて最もエネルギーの高い位相点R_iを求めるために，ポテンシャル力とバネに平行な成分の符号を反転させたものの合力

  

$$F_{\text{Cl}}=F_i^{\nabla }-2\left(F_i^{\nabla }\cdot {\widehat{\tau }}_i\right){\widehat{\tau }}_i$$(7)

を計算し，これを最小化することで鞍点を求める。

また，従来のNEB法では単位胞の格子形状は不変の条件下で計算されるため，相変態のような格子変形を自由度として含む遷移過程を考慮することができない。そこで，Sheppardら²¹⁾により格子の変形も考慮したNEB法であるG-SSNEB法が考案された。G-SSNEB法では，各位相点の変位ベクトルΔℝ_iを格子変形ɛ_iと単位胞中の原子座標の変位ΔR_iの二つの変数を用いて表現する。

  

$$\Delta {ℝ}_i=\left\{J{\epsilon }_i,\Delta R_i\right\}$$(8)

式(8)中のɛ_iは式(9)で表されるヤコビアンJを用いてΔR_iとスケールが一致するように補正される。

  

$$J={\Omega }^{1/3}{N}^{1/6}$$(9)

ここで，Ωは単位胞の体積，Nは単位胞内の原子数である。NEBと同様に，位相点にかかる力${\mathbb{F}}_i^{\text{G}-\text{SSNEB}}$はF_i^NEBと格子変形に対する応力テンソルσ_iを用いて式(10)のように表される。

  

$${\mathbb{F}}_i^{\text{G}-\text{SSNEB}}=\left(-\frac{\Omega {\sigma }_i}{J},F_i^{\text{NEB}}\right)$$(10)

G-SSNEB法では各位相の${\mathbb{F}}_i^{\text{G}-\text{SSNEB}}$を最小化することで，MEPを求める。本研究ではVASPコードに組み込まれたG-SSNEB法を用いて相変態のMEPを求めた。

3. 結果と考察

3・1　炭素原子の固溶サイトの検討

BCC構造では，Fig.1に示した三種類の炭素原子の固溶サイトは完全に等価であるが，c/a>1の場合はa，bサイトは狭く，cサイトの空隙は広くなる。したがって，軸比に応じた空隙サイズの変化に伴って，cサイトへの固溶状態のエネルギーとa，bサイトへの固溶状態のエネルギー差は変化する。Table 1に，aあるいはbサイトに炭素原子が固溶した状態とcサイトに炭素原子が固溶した状態のエネルギー差を求めた結果を示した。計算にはBCC構造の単位胞を3×3×3に拡張したスーパーセルを用いて，炭素原子一つを八面体侵入位置に配置した構造を使用した。本研究のクラスター展開法では原子位置の局所緩和を行わない規則構造を用いたが，その条件におけるエネルギー差をTable 1の第2列に示した。c/a=1.01，1.05，1.10のいずれの軸比においてもcサイトに固溶した状態が熱力学的に大きく安定である。最もエネルギー差が小さいc/a=1.01(ΔE=0.680 eV)のエネルギー差から温度T=400 Kにおけるa，bサイトとcサイトの炭素固溶の分配比をボルツマン因子から求めると，cサイトの固溶量に対するa，bサイトの固溶量は2.63×10^-9倍となり，無視できるほど小さい。この理由により，本研究ではa，bサイトへ炭素原子が固溶しない構造モデルを用いた。ただし，このエネルギー差について，本研究では考慮していない局所原子位置の緩和が影響することに留意しておきたい。局所原子位置の緩和を行なった場合，近すぎるFe-C原子間距離は適切な距離まで変位して安定化するために，a，bサイトに固溶した場合とcサイトに固溶した場合のエネルギー差はより小さくなる。Table 1第3列には局所原子位置の緩和を行なった条件で計算した固溶状態のエネルギー差を示した。c/a=1.01におけるエネルギー差ΔEは0.054 eVまで低下するため，この値を用いてT=400 Kのa，bサイト中の炭素固溶量を計算するとcサイトの炭素固溶量に対して0.208倍となる。一方，c=1.05におけるΔEは0.216 eVであり，この値から同様に計算される分配比は1.89×10^-4と依然として小さい値となる。したがって，軸比の小さな領域では，固溶サイト間のエネルギー差が小さいためにa，bサイトにも固溶した状態が起こり得るが，前述のように本研究のモデルではこの効果は取り入れていない。また，3・2節にてマルテンサイトの軸比の炭素濃度依存性について改めて議論するが，実験ではおよそ2.75%以上の炭素濃度では軸比が炭素濃度増加に伴い増加し，それ以下の炭素濃度においては軸比が発現せず立方晶をとることが確認されている²²⁾。このような実験事実からも低炭素濃度領域ではa，b，c三つのサイトに混在して固溶することが示唆される。一方，十分大きな軸比が発現しているような炭素濃度が高いマルテンサイトにおいては，この緩和のエネルギーを考慮してもcサイトに固溶した方が熱力学的に安定であり，本研究で採用した計算モデルのように炭素原子はcサイトだけを占有すると考えられる。

Table 1. Energy difference between the cases in which carbon atom is located in a or b site and that in c site. The second column shows the results without considering the relaxation of the ionic positions, while the third one shows the results when the ionic positions are relaxed. The unit is eV.

c/a	Without relaxation	With relaxation
1.01	0.680	0.054
1.05	3.823	0.214
1.10	4.147	0.324

3・2　CE-CVMの計算結果

鉄-炭素二元系マルテンサイトの自由エネルギーをT=400 Kの温度条件で計算した結果をFig.3に示した。Fig.3(a)における横軸は原子分率で表した炭素濃度であり，終端組成は全てのcサイトに炭素が固溶した状態のFeCとなる。鉄および炭素のエネルギーの基準状態はBCC-Feおよびグラファイトとした。またFig.3(b)は低炭素濃度領域の自由エネルギーの拡大図であるが，ここには各軸比の自由エネルギーの交点を白抜きの丸で示した。炭素の組成分率をx_Cとすると，およそx_C=0.012以下の組成範囲ではc/a=1.01が，x_C=0.012-0.029の範囲ではc/a=1.05，x_C=0.029以上の組成範囲ではc/a=1.10がエネルギー最小となることがわかる。このように炭素濃度の増加に伴って，軸比が高い状態が安定化する計算結果が得られた。この結果は，定性的にはマルテンサイトのc/aが炭素濃度の増加とともに増加する実験事実²²⁾と整合する。しかし，本研究による炭素濃度に対する軸比の増加量は実験によるc/a=1+x_Cの関係よりも傾きが大きく，軸比の変化量がより過剰に評価されている。この原因は，本研究で考慮しなかった原子位置の局所緩和の影響が考えられる。Table 2に局所原子位置緩和なしとありの条件で行なった平衡格子定数および単位原子あたりの全エネルギーの比較を示す。構造モデルは3・1節で用いたものと同様，BCC構造の単位胞を3×3×3に拡張したスーパーセルによって作成したFe₅₄Cを用いた。Table 2の格子定数に注目すると，a軸とc軸のいずれもが局所緩和を行なった計算結果の方が小さくなることがわかる。また全エネルギーは局所緩和ありの計算では軸比1.01が計算した中で最もエネルギーが低い結果となったが，局所緩和なしの計算では軸比の増加とともにエネルギーが低下する。この計算結果は，以下のように解釈される。すなわち2・1節で述べたようにcサイトに存在する炭素原子は上下の鉄原子との原子間距離が最も近く，八面体FeはZ軸方向に押しつぶされた形状をしている。そのためBCCに炭素が固溶すると，近すぎる最隣接Fe-Cの原子間距離が適切な距離まで離れることでエネルギーが低下する。原子位置を格子点上に固定した本研究における計算モデルの場合は，c/aを大きくすることでしかこの最隣接のFe-C原子間距離を大きくできないが，原子位置を緩和する場合は八面体上下の鉄原子位置を局所的に変位させることで，炭素と鉄原子間の距離を変えることが可能である。そのため，局所原子位置を緩和した場合，小さな軸比においても結晶格子中の鉄原子が変位することにより炭素と鉄原子間距離が適切な距離をとることが可能になると考えられる。したがって最安定の軸比は格子点上に固定して計算したものより小さくなる。また，局所緩和ありとなしの全エネルギーの差から緩和に伴うエネルギーの低下量を計算すると，c/a=1.01では-0.087 eV，c/a=1.05では-0.069 eV，c/a=1.10では-0.051 eVと軸比が小さい構造ほどエネルギーの低下量が大きいことがわかる。このことから，自由エネルギーは局所緩和の効果を取り入れることで軸比が小さいものほど大きく安定化すると予想され，Fig.3(b)における自由エネルギーの交点は高炭素濃度側にシフトすると考えられる。したがって，局所緩和の効果を導入することで実験結果との乖離を解消する傾向にあると期待される。これを検証する目的で，スーパーセル法によって炭素濃度を変えた構造モデルを用いて，最安定な軸比と炭素濃度の関係を局所原子位置の緩和有無に対して比較した結果をFig.4に示す。計算ではBCC構造を4×4×4に拡張したスーパーセル中に炭素原子を1つ配置したFe₁₂₈Cおよび3×3×3に拡張したスーパーセル中に炭素原子を1つあるいは2つ配置したFe₅₄C，Fe₅₄C₂の3種類の構造を用いた。黒丸で示したものが局所原子位置を固定して計算したもので，白丸で示したものは局所原子位置を緩和した計算結果である。局所原子位置を固定した計算結果に比べて，局所原子位置の緩和を行うと安定な軸比が低下し，点線で示したc/a=1+x_Cの経験式をよく再現する。また，このような構造緩和を行った構造モデルで定量的にも実験結果を再現する同様の結果は先行研究においても報告されている^23,24)。したがって本研究においても，CE-CVM計算に局所原子位置緩和を取り入れることで，c/aの炭素濃度依存性の関係を再現可能になると推定される。

Fig. 3.

(a) Free energy curves of BCT martensitic structure in the Fe-C system at T=400 K by means of CE-CVM. Calculation was performed for c/a=1.01, 1.05, and 1.10. Figure (b) shows the enlarged portion of lower concentration region of figure (a).

Table 2. Comparison of lattice parameters and total energies of BCT-Fe₅₄C structure before and after relaxing the ionic positions.

	Without relaxation			With relaxation
c/a	a (Å)	c (Å)	Energy (eV)	a (Å)	c (Å)	Energy (eV)
1.01	2.843	2.872	–8.227	2.821	2.849	–8.314
1.05	2.805	2.945	–8.242	2.785	2.924	–8.311
1.10	2.764	3.041	–8.249	2.750	3.025	–8.300

Fig. 4.

Variation of c/a ratio with carbon concentration in BCT Fe-C. The filled and open circles represent the calculated results with and without relaxation of the ionic positions, respectively. The dashed line denotes the empirical relationship given by c/a = 1 + x_c.

Fig.3(a)に示した自由エネルギーは，グラファイト基準では炭素濃度増加に伴い増加するが，エネルギーの増加率は単調ではなく，低炭素濃度領域の自由エネルギーは上に凸の炭素濃度依存性を示している。またいずれの軸比においても約30%の炭素濃度領域の自由エネルギーが下にたわんだ形状をしており，この組成と低炭素濃度の間に二相分離傾向が存在する。

この二相分離の原因を調べるために，クラスター展開・変分法に用いたBCT基規則構造のエネルギーをFig.5にプロットした。エネルギーの基準は，各軸比を有するBCTの純Feと全てのcサイトに炭素が固溶したBCT構造のFeCにとってある。また各規則構造のエネルギーは，それぞれの規則構造の平衡体積位置でのエネルギーである。図中の横軸はcサイト中での炭素のサイト占有率で表記しており，終端組成はFeCに対応する。いずれの軸比においてもサイト占有率が0.5となるFeC_0.5の組成，すなわち図中に白丸で示したFe₂Cの組成比で安定な規則構造が存在することがわかる。なお，このBCT格子中での最安定規則構造を以後BCT-Fe₂Cと表記することにする。これらの図において，FeとBCT-Fe₂Cの間の組成領域の規則構造のエネルギーは，FeとBCT-Fe₂Cのエネルギーを結んだ直線よりも高い位置を占めている。したがって，自由エネルギーはこのBCT-Fe₂C規則構造への規則化に由来して，約30%の組成領域でエネルギーが低下していると考えられ，この規則化に伴ってFeとBCT-Fe₂C間に二相分離が生じているものと考えられる。

Fig. 5.

Formation energies of BCT-based ordered structures varying with site fraction of carbon for: (a) c/a = 1.01, (b) c/a = 1.05, and (c) c/a = 1.10.

CE-CVM計算から得られた自由エネルギーにみられる二相分離領域の内側にはスピノーダル分解領域が存在すると考えられる。このような領域は自由エネルギーの組成に対する二階微分の正負から判別できるが，本研究で評価した自由エネルギーは軸比を断続的に変えて計算したものであり，この評価方法を適用することができない。一方，Renら⁴⁾はおよそ8.5 at.%Cの試料を室温時効することで，鉄-炭素二元系マルテンサイト中にスピノーダル分解による炭素原子のクラスターの生成を報告している。本研究からは，正確なスピノーダル領域の判定は行えないが，Renらの報告によるスピノーダル分解は，本研究で得られた二相分離傾向に由来している可能性がある。以上のように，低温時効によって確認されている炭素原子のクラスタリングは，マルテンサイトの持つ二相分離傾向に起因している可能性があり，この二相分離はBCT-Fe₂Cへの規則化に駆動されて進行するものと考えられる。

3・3　NEB法によるBCT-Fe₂C規則構造からη炭化物への遷移エネルギー

前節では，安定なBCT-Fe₂C規則構造の存在に由来して，鉄-炭素二元系BCT固溶体中にFeとFe₂Cの二相分離傾向が存在する結果を示した。この規則構造は実験で報告されているη炭化物と組成比が一致することから，本研究ではBCT-Fe₂C規則構造とη炭化物との関連性に着目した。そこで，まず両者の結晶構造の観点から考察した。

Fig.6に，Fig.5に示したBCT基規則構造のエネルギーの中で最も安定な規則構造であるBCT-Fe₂C規則構造，およびη炭化物の結晶構造を示した。Fig.6に示したBCT-Fe₂Cは，η炭化物との比較が容易になるように格子軸の方位を取り直した。図中の新たな方位軸a，b，cは，BCT構造の格子ベクトルa'，b'，c'と以下の関係式で示される。

  

$$\left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& 0\\ -1& 0& 1\end{array}\right]\left[\begin{array}{c}{a}^{\prime }\\ b^{\prime }\\ c^{\prime }\end{array}\right]$$(10)

Fig. 6.

Comparison of crystal structures of (a) BCT-Fe₂C and (b) η-Fe₂C. The black and gray circles represent Fe and C atoms, respectively.

また，図中のBCT-Fe₂C構造はBCTのc/aが1.05の構造のものを示した。また，この方位軸で表記したBCT-Fe₂Cおよびη-Fe₂Cの格子定数，軸角，原子座標をTable 3に示した。η-Fe₂Cの格子定数はBCT-Fe₂Cに比べa，b軸が短く，c軸が長くなっており斜方晶となる。一方，分率表記で表した原子座標をみると，a軸方向のFe原子位置が両者でわずかに異なり，さらにBCT-Fe₂CではC面心位置に存在する炭素原子が，η-Fe₂Cでは体心位置に存在する。炭素原子位置を除いて，両者の結晶構造におけるこれらの違いは僅かであるため，面心位置から体心位置への炭素原子の移動と鉄原子位置の緩和によって，BCT-Fe₂C構造からη炭化物へ遷移する可能性がある。この遷移過程が起こり得るかを判断するには，遷移過程におけるエネルギー障壁の計算が必要である。そこで，この遷移障壁をG-SSNEB法により計算した。

Table 3. Comparison of lattice parameters and atomic coordinates of BCT-Fe₂C and η-Fe₂C. Differences in atomic coordinates are shown in bold text.

Structure	Lattice parameter	Element	Fractional coordinate
			x	y	z
BCT-Fe2C	a = 4.38 Å b = 3.02 Å c = 4.38 Å α=β=γ=90°	Fe 1	0.75	0	0.75
		Fe 2	0.25	0	0.25
		Fe 3	0.25	0.5	0.75
		Fe 4	0.75	0.5	0.25
		C 1	0	0	0
		C 2	0.5	0.5	0
η-Fe2C	a = 4.71 Å b = 2.83 Å c = 4.28 Å α=β=γ=90°	Fe 1	0.65	0	0.75
		Fe 2	0.34	0	0.25
		Fe 3	0.15	0.5	0.75
		Fe 4	0.84	0.5	0.25
		C 1	0	0	0
		C 2	0.5	0.5	0.5

G-SSNEB法により求めた各軸比のBCT-Fe₂C規則構造からη炭化物への遷移経路のエネルギーをFig.7に示した。軸比の異なる三つの構造を初期構造として用いて計算を行ない，c/a=1.01の結果を○で，c/a=1.05の結果を□で，c/a=1.10の結果を△で示した。また縦軸のエネルギーはBCT-Fe₂C規則構造のエネルギーを基準としており，図中の記号ⓐ~ⓓはそれぞれⓐ:BCT-Fe₂C規則構造，ⓑ:遷移過程における最安定構造，ⓒ:遷移過程の鞍点，ⓓ:η炭化物を指している。横軸は遷移に伴う構造イメージの番号である。初期構造に用いた構造の軸比によって，遷移過程のエネルギー状態に若干の違いが見られるが，いずれの軸比においても類似したエネルギー変化を示した。Fig.8にFig.7におけるⓐ~ⓓ点のb軸方向から見た結晶構造を示した。またFig.8は初期状態としてc/a=1.05の構造を用いた計算結果を示したが，他の軸比を初期構造とした場合においても同様の構造変化の過程が確認された。η炭化物への遷移では，初めにⓐ点に示したBCT-Fe₂C規則構造が格子の形状と原子位置を緩和し局所的な安定構造であるⓑ点の構造を経由した後，丸で囲った面心位置の炭素が体心位置へと移動していく。丸で記した炭素と格子内の鉄原子の距離が最も近くなるⓒ点がエネルギー曲面における鞍点であり，BCT-Fe₂C規則構造からη炭化物へのエネルギー障壁となる。ⓑ点とⓒ点のエネルギー差から求めたエネルギー障壁の高さは，用いた初期構造の軸比に応じて若干の違いがみられるが，計算に用いた単位胞(Fe₄C₂)あたり0.9-1.10 eVとなる。面心位置から体心位置に大きく変位する炭素原子は単位胞あたり一つ存在するので，このエネルギーはこの炭素原子が移動に伴う障壁とみなせる。フェライト中の炭素の拡散障壁は0.86 eVであると報告されており²³⁾，この値とほぼ同じ大きさである。したがって，炭素の拡散が十分起こりうる温度領域では，BCT格子からη炭化物への遷移障壁も十分に乗り越えることが可能であると考えられる。

Fig. 7.

Energy change with the structure transition from BCT-Fe₂C to η-Fe₂C calculated using G-SSNEB method. The unit of vertical axis is eV per formula unit (f.u.) of Fe₄C₂.

Fig. 8.

Crystal structures corresponding to the points ⓐ through ⓓ in Fig.7 observed from b-axis direction. The black and gray circles represent Fe and C atoms, respectively. The initial state of axial ratio c/a was fixed to 1.05.

以上の結果より，BCTマルテンサイトの低温焼戻し過程におけるη炭化物は二段階の過程を経て析出する可能性が示された。すなわち，(1)BCTマルテンサイトが示す二相分離傾向によりBCT-Fe₂C規則構造が形成され，(2)その状態からエネルギー障壁を超えてη炭化物に遷移するという過程である。

このような解釈は低温焼戻し過程において，より安定な炭化物であるセメンタイトよりも先にη炭化物が析出する理由をよく説明できる。すなわちBCTマルテンサイト中の炭素原子は，BCTマルテンサイトの二相分離によってFe₂Cの組成まで濃化するため，セメンタイトの炭素濃度を超えた状態が熱力学的に実現しやすいと考えられる。そのため，炭素濃度が近いη炭化物が優先的に生成すると理解される。このようにBCTマルテンサイトのもつ熱力学的性質が，相対的に不安定である高濃度のη炭化物を生成するという一見不可解な振る舞いの一因になっていると推察される。一方，本研究ではBCTのFe-C固溶体からセメンタイト構造への遷移状態の評価を行なっておらず，変態過程も含めたη炭化物およびセメンタイトの生成過程の比較に至っていない。今後，Fe₃Cセメンタイト構造に関しても，固溶体中の構造を経由して生成する過程を計算し比較することで，より詳細な議論が可能になると期待される。

本研究ではη炭化物を一つの例として解析を進めたが，このような解析や考察は他の炭化物や窒化物についても可能であると考えられ，本手法のような研究を展開することによって，鋼中の炭化物や窒化物の生成過程に関するより幅広い議論が可能になることが期待される。一方，本研究では原子の局所緩和を取り入れていないために，今後の検討が必要な項目もある。例えば，低軸比/低炭素濃度領域ではa，bサイトへの炭素原子固溶が生じると考えられるが，本研究のモデルではa，bサイトへの固溶が過剰に不安定的に評価される傾向にあり，この効果を取り入れることができなかった。またこれに関連して，炭素濃度-軸比の関係は実験で報告されている傾きよりも大きく評価されている。これらの課題は，原子位置の緩和を考慮したモデルを適用することによって改善できると期待される。

4. 結論

本研究では，鉄-炭素二元系BCTマルテンサイトの自由エネルギーを第一原理計算とクラスター展開・変分法を用いて計算した。さらにG-SSNEB法を用いて規則化したBCT-Fe₂C構造から，η-Fe₂C炭化物への構造遷移の活性化エネルギーを求めた。本研究により得られた知見は以下の通りである。

(1)BCTマルテンサイトの自由エネルギー曲線にはBCT-Fe₂C規則構造の形成に伴う二相分離が存在し，これが低温焼戻し過程で炭素のクラスタリングを促進する可能性が示された。

(2)BCTマルテンサイト中で準安定的に生成するBCT-Fe₂C規則構造は，η炭化物と結晶学的な共通点の多い類似した構造であることが明らかになった。

(3)BCT-Fe₂C規則構造からη炭化物への遷移に伴う鉄原子および炭素原子の移動に要するエネルギー障壁をG-SSNEB法から求めた結果，活性化エネルギーは単位胞(Fe₄C₂)あたり0.9-1.10 eVであり，BCT-Fe₂C規則構造を介してη炭化物が生成する可能性が示唆された。

(4)以上の結果より，BCTマルテンサイトの低温焼戻し過程におけるη炭化物は，BCTマルテンサイトが示す二相分離傾向によりBCT-Fe₂C規則構造が形成され，その状態からエネルギー障壁を超えてη炭化物に遷移するという二段階の過程を経て析出する可能性が示された。

謝辞

この成果は，国立研究開発法人科学技術振興機構 産学共創基礎基盤研究プログラムの「革新的構造用金属材料創製を目指したヘテロ構造制御に基づく新指導原理の構築」研究課題，JSPS科研費16H02378と17H01330，および公益財団法人JFE21世紀財団2019年度技術研究助成の助成を受けて得られたものです。ここに謝意を表します。

文献

• 1)  Y. Imai and S. Sawada: Bull. Jpn. Inst. Met., 8(1969), 820 (in Japanese).
• 2)  S. Nagakura, Y. Hirotsu, M. Kusunoki, T. Suzuki and Y. Nakamura: Metall. Trans. A, 14(1983), 1025.
• 3)  K.A. Taylor, L. Chang, G.B. Olson, G.D.W. Smith, M. Cohen and J.B. Vander Sande: Metall. Trans. A, 20(1989), 2717.
• 4)  S.B. Ren, T. Tadaki, K. Shimizu and X.T. Wang: Metall. Mater. Trans. A, 26(1995), 2001.
• 5)  Y. Hirotsu and S. Nagakura: Acta Metall., 20(1972), 645.
• 6)  Phase Transformations in Steels, Vol.1, ed. by E.Pereloma and D.V.Edmonds, Woodhead Publishing, Cambridge, U.K., (2012).
• 7)  G. Hägg: Z. Krist., 89(1934), 92.
• 8)  S. Nagakura, S. Suzuki and M. Kusunoki: Trans. Jpn. Inst. Met., 22(1981), 699.
• 9)  C.M. Fang, M.H.F. Sluiter, M.A. Van Huis, C.K. Ande and H.W. Zandbergen: Phys. Rev. Lett., 105(2010), 1.
• 10)  E.C. Bain: Trans. Metall. Soc. AIME, 70(1924), 25.
• 11)  M.H.F. Sluiter, C. Colinet and A. Pasturel: Phys. Rev. B, 73(2006), 174204.
• 12)  R. Kikuchi: Phys. Rev., 81(1951), 988.
• 13)  G. Kresse and J. Furthmüller: Phys. Rev. B, 54(1996), 11169.
• 14)  G. Kresse and D. Joubert: Phys. Rev. B, 59(1999), 1758.
• 15)  P.E. Blöchl: Phys. Rev. B, 50(1994), 17953.
• 16)  J.P. Perdew, A. Ruzsinsky, G.I. Csonka, O.A. Vydrov, G.E. Scuseria, L.A. Constantin, X. Zhou and K. Burke: Phys. Rev. Lett., 100(2008), 136406.; Erratum, Phys. Rev. Lett., 102(2009), 039902.
• 17)  J.P. Perdew, K. Burke and M. Ernzerhof: Phys. Rev. Lett., 77(1996), 3865.
• 18)  Classical and Quantum Dynamics in Condensed Phase Simulations, ed. by B.J Berne, World Scientific, Singapore, (1998).
• 19)  D. Sheppard, R. Terrell and G. Henkelman: J. Chem. Phys., 128(2008), 134106.
• 20)  G. Henkelman and H. Jónsson: J. Chem. Phys., 113(2000), 9978.
• 21)  D. Sheppard, P. Xiao, W. Chemelewski, D.D. Johnson and G. Henkelman: J. Chem. Phys., 136(2012), 074103.
• 22)  O.D. Sherby, J. Wadsworth, D.R. Lesuer and C.K. Syn: Mater. Trans., 49(2008), 2016.
• 23)  D.E. Jiang and E.A. Carter: Phys. Rev. B, 67(2003), 214103.
• 24)  H. Ohtsuka, V.A. Dinh, T. Ohno, K. Tsuzaki, K. Tsuchiya, R. Sahara, H. Kitazawa and T. Nakamura: Tetsu-to-Hagané, 100(2014), 1329 (in Japanese).