王輝，蘇福永，劉文麗，溫治

(北京科技大學能源與環(huán)境工程學院，北京，100083)

碳是鐵基合金中最常見的間隙原子之一，通過形成碳化物[1]、晶界碳偏析[2]、硬化和脆化效應[3]等，可顯著提高鋼的強度和硬度[4]。碳在鐵中的分布和擴散可以影響鋼中許多過程，如碳化物沉淀、馬氏體老化和鐵素體轉變等[5]。從微觀上研究鐵中碳原子的擴散行為，有利于更好地解釋碳原子對鐵合金性能的影響規(guī)律和機理。只有明確碳原子在鐵晶體內(nèi)的擴散路徑，才能通過改進工藝設計和合金成分來減弱或增加碳原子的吸附和擴散能力，改善金屬性能。

為改善鋼的性能，常在鐵基合金添加元素。例如，添加Cr 元素可使鋼材具有更好淬透性，增加其碳化物的穩(wěn)定性，延長鋼材的使用壽命。封輝等[6]發(fā)現(xiàn)隨著Cr元素質量分數(shù)增加，低合金高強耐候鋼的強度上升，塑性下降。Fe-Cr 多層體系還因其巨磁阻效應而廣泛應用于存儲器件中[7?8]。侯瑞東等[9]通過實驗研究發(fā)現(xiàn)，900～950 ℃溫度條件下，隨Cr 質量分數(shù)增加，鋼的抗氧化性增強。然而，Cr 元素的加入給鐵基合金帶來一系列優(yōu)良性質的同時也將與其中的其他原子發(fā)生相互作用，例如它會與碳原子之間產(chǎn)生化學或電子效應[4]，從而影響鋼中碳的行為，如碳溶解度極限[10]，顯微偏析[11]，擴散[12]和碳化物沉淀[13]等。

迄今為止，基于密度泛函理論(DFT)的第一性原理可以準確評估原子相互作用，是研究一些基本原子現(xiàn)象的有力工具。JIANG等[14]采用該方法計算了C 在α-Fe 和γ-Fe 中的溶解和擴散，發(fā)現(xiàn)α-Fe中2 個八面體位置之間的碳擴散的最小能量路徑(MEP)具有0.86 eV 的勢壘，與實驗值0.87 eV 接近。HERSCHBERG 等[15]計算在不同局部化學環(huán)境中1組C擴散的結合能和遷移能壘，采用動力學蒙特卡羅方法模擬Fe-Cr合金中C的擴散特征，發(fā)現(xiàn)在Fe-Cr合金中C擴散的平均遷移勢壘隨著Cr濃度增加而逐漸增加。孫興廣等[16]對Fe-Cr-Ni合金中點缺陷相互作用能的計算表明，相鄰空位表現(xiàn)出相互吸引作用，會促進孔洞形成。VON APPEN 等[17]利用密度泛函理論研究了奧氏體Fe-Mn 合金中的氫間隙，發(fā)現(xiàn)八面體H的間隙位置更傾向于含Mn原子的局部環(huán)境而不是Fe原子。SANDERG等[18]采用第一性原理密度泛函理論計算了體心立方Fe 和Cr中C溶液焓和擴散活化焓，2個體系中計算得到的焓與實驗焓的偏差均小于0.05 eV。SIMONOVIC等[19]計算了間隙碳與取代硅的相互作用以及這種相互作用對體心立方鐵中碳的擴散的影響，并預測計算擴散的活化能和擴散前因子。

目前，雖然學者們研究了α-Fe 中的碳擴散特征，但是都沒有詳細解釋α-Fe 中Cr 原子和C 原子的相互作用和鉻對碳擴散影響。因此，本文首先基于第一性原理方法計算了α-Fe 中Cr 原子和C 原子的相互作用能；其次，通過態(tài)密度計算分析了各原子之間的成鍵情況以及其對結構穩(wěn)定性的影響；最后，用攀升圖像微擾彈性帶(CI-NEB)方法對α-Fe 中C 原子的擴散特征進行了分析，同時還研究了Cr原子對α-Fe中C原子的擴散的影響。

1 計算方法與結構模型

本文中第一性原理計算基于自旋極化密度泛函理論[20?21]，在維也納大學開發(fā)的Vienna Ab initio simulation package(VASP)[22?23]軟件包中計算。利用VASP 求解具有周期性邊界條件和平面波基組的Kohn-Sham 方程。離子和電子的相互作用采用投影綴加波方法(PAW)[24?25]計算，電子交換關聯(lián)泛函采用廣義梯度近似GGA[26]下的PBE 泛函形式。通過收斂測試，平面波截斷能取450 eV，使能量收斂到每個原子0.1 meV 以內(nèi)。布里淵區(qū)積分采用Г點為中心的Monhkorst-Pack方法產(chǎn)生的k網(wǎng)格點方法，在2×2×2(16 個原子)，3×3×3(54 個原子)和4×4×4(128 個原子)的α-Fe 超胞的計算中分別采用選取8×8×8，6×6×6 和4×4×4 的k 網(wǎng)格，使原子所受力小于0.2 eV/nm。

本文采用HENKELMAN等[27?28]開發(fā)的過渡態(tài)包(VTST)中的攀升圖像微擾彈性帶(climbing image nudged elastic band，CI-NEB)方法計算C 原子擴散過程中的能量壁壘，這種方法也提供了一種尋找最小能量路徑(MEP)的方法。計算過程中圖像弛豫的收斂標準為作用于原子上的力小于0.1 eV/nm。

C原子的原子半徑較小，通常認為在α-Fe中占據(jù)間隙位置。為了確定C原子占據(jù)的是四面體位置還是八面體位置，本文計算了C 原子在α-Fe 中的溶解焓ΔHs，其表達式為

式中：E(FenC)為包含n個Fe原子和1個C原子的α-Fe 晶胞的能量，eV；nE(Fe)為包含n個Fe 原子的α-Fe 晶胞的能量，eV；E(C)為單個石墨原子的能量，eV。為了能更好消除誤差，E(FenC)和nE(Fe)這2 項能量采用相同的參數(shù)(截斷能及K 點等)進行計算。石墨分子采用8×8×4的網(wǎng)格，450 eV的截斷能進行計算。

當α-Fe 中溶質原子X和Y原子之間相距Rs時，有效相互作用能表達式為[19]

式中：s 為特定相鄰位置的殼層；E(FenXY,Rs)為溶質原子X和Y之間相距Rs時體系的能量，eV；因為計算所采用的晶胞尺寸有限且計算時周期性的邊界條件，所以某些相對位置的計算中溶質原子X和Y之間會發(fā)生多次相互作用，本文用m來表示這種情況。ΔE(X)和ΔE(Y)分別為單個溶質原子X和Y的能量，eV。其表達式為[19]：

式中：E(FenX)為包含n個鐵原子和一個溶質原子X的α-Fe晶胞的能量。

2 計算結果與討論

2.1 α-Fe（C）-Cr體系結構分析

基于密度泛函理論第一性原理，分別采用超軟贗勢(OTFG ultrasoft)和投影綴加波(PAW)方法對3×3×3的α-Fe超胞進行了優(yōu)化和能量計算，計算得到的純鐵磁性α-Fe 的晶格邊長均為0.283 nm，每個Fe 原子的磁矩分別為2.18μB和2.21μB(μB=9.272 6×10?24A·m2)，如表1所示。這與之前學者計算得到的結果相近，但投影綴加波方法計算得到的結果和實驗結果符合更好，故本文的計算均采用投影綴加波方法。

表1 UUSP和PAW計算結果與文獻結果的對比Table 1 Comparison of UUSP and PAW calculation results with literature results

由于Cr原子的半徑與Fe原子的半徑相似，當Cr原子加入到α-Fe中時，會直接替代Fe原子的位置。α-Fe 中Fe 的占位有2 種，分別為體心位置和頂角位置。所以，為了確定Cr 原子在α-Fe 中的占位，本文分別計算了Cr 原子在2×2×2 的α-Fe 晶胞中占據(jù)2種不同位置時晶胞的體積和總能量，然后計算了Cr 原子替代Fe 原子后體系的體積變化率，體積變化率ε計算公式如下：

式中：V為Fe15Cr 晶胞的體積；V0為Fe16晶胞的體積。計算結果如表2所示。

表2 Cr原子占位分析Table 2 Cr atomic occupancy analysis

由表2可見，Cr 原子加入到α-Fe 中引起了晶格畸變，Cr 原子占據(jù)體心位置時的體積變化率略小于占據(jù)頂角位置時的體積變化率，但差別不大，說明Cr 原子在α-Fe 中可能是隨機占位。此外，Cr原子占據(jù)2種位置時體系總能量沒有差別，這也可以驗證上述觀點。由于Cr 原子占據(jù)體心位置時的體積變化率略小于占據(jù)頂角位置時的體積變化率，本文后續(xù)的工作都按Cr 原子占據(jù)體心位置來進行計算。通過對Fe(Cr)體系的計算，Cr 原子在α-Fe中表現(xiàn)為反鐵磁性，后續(xù)計算中Fe 原子按照鐵磁性，Cr原子按照反鐵磁性設置。

C原子的原子半徑較小，通常認為在α-Fe中占據(jù)間隙位置。為了確定C原子的具體占位情況，本文分別在八面體和四面體間隙位置插入了C原子，并計算了C原子在這2種位置的溶解焓，結果如表3所示。由表3可見，C 原子占據(jù)八面體位置時體系的總能量和溶解焓相比于占據(jù)四面體位置時的低，因此，C在α-Fe中優(yōu)先占據(jù)八面體位置，這與DUNLOP等[13]的計算結果一致。

表3 C原子占位分析Table 3 C atomic occupancy analysis

Cr 元素給鐵基合金帶來一系列優(yōu)良性質可能與基體的鍵合作用有關。建立了Fe15CrC的α-Fe-Cr合金體系，并分析了其電子結構，計算了α-Fe-Cr合金體中各原子的分態(tài)密度，分析鍵合作用，結果如圖1所示。

圖1 α-Fe(C)-Cr合金體系的態(tài)密度圖(DOS)和分態(tài)密度圖(PDOS)Fig.1 Density of states(DOS)and partial density of states(PDOS)of α-Fe(C)-Cr alloy systems

從圖1可見：α-Fe-Cr 合金體系成鍵電子主要分布在?75.0～?71.9，?45.0～?42.0，?12.5～?12.0 和?8.0～15.0 eV 等4 個區(qū)間，而從圖1中費米能級處有很明顯的成鍵峰的是Fe 3d 軌道和Cr 3d 軌道，這2 條軌道的雜化導致DOS 圖費米能級處出現(xiàn)了尖銳波峰，而C在此處的態(tài)密度幾乎為0，因此Fe和Cr 在此處強烈成鍵，C 原子雖然在?12.5～?12.0 eV 區(qū)間內(nèi)提供成鍵電子，但Fe 和Cr 在這2個區(qū)間的態(tài)密度卻很小，因此，在α-Fe(C)-Cr合金中，3 種原子都提供成鍵電子，但Fe 和Cr 是成鍵的主要因素。

態(tài)密度能判斷原子間成鍵情況，也能夠具體分析成鍵的軌道，但不能對成鍵進行量化描述，因此，為更直觀地揭示合金元素Cr 與Fe，C 原子間的鍵合強度，從電子結構方面進一步量化成鍵強度，計算了α-Fe(C)-Cr 的Bader 電荷布居，結果如表4所示。

表4 α-Fe(C)-Cr的barder電荷布居Table 4 Barder charge population of α-Fe(C)-Cr

從表4可見：在α-Fe(C)-Cr晶體中Cr原子的電荷數(shù)為11.495 9，小于純鉻晶體中電荷數(shù)12，說明Cr原子加入到α-Fe(C)體系中后失電子，F(xiàn)e原子一部分得到電子，一部分失去電子，但總體來說表現(xiàn)為失電子，只有C原子的電荷數(shù)增加了1.146 9，這說明Cr原子加入α-Fe基體內(nèi)后，Cr原子和Fe原子失去電子，C得到電子，三者之間形成離子鍵的相互作用，導致鐵素體晶胞得到強化。

2.2 α-Fe中Cr與C的相互作用

為了進一步分析Cr原子加入α-Fe(C)體系后對C原子的影響，本節(jié)在4×4×4的α-Fe晶胞中建立了Cr 原子和不同近鄰位置C 原子的結構示意圖，如圖2所示。圖2給出了2 個5 近鄰位置的位置，分別用5 和5′表示，圖2中Cr-C5 之間上沒有其他原子，而Cr-C5'之間有1個Fe原子。

圖2 α-Fe中Cr原子和不同近鄰位置C原子的結構示意圖Fig.2 Schematic diagram of Cr atom structure in α-Fe crystal and C atom in its different neighboring positions

對含有Cr 原子且在其第一近鄰八面體位置加入C原子的α-Fe晶體進行結構優(yōu)化。圖3所示為結構優(yōu)化前后α-Fe 晶體中的Cr 原子和C 原子所在的<010>方向的切面上的結構圖，其中，藍色的圓球代表Fe原子。從圖3可見，與結構優(yōu)化之前相比，結構優(yōu)化后Cr 原子沿<010>方向向左邊移動，Cr原子和C 原子之間的距離由0.142 nm 增加到0.182 nm，這說明Cr 原子與在其第一近鄰八面體位置的C原子存在排斥力。

圖3 α-Fe-Cr-C晶胞結構優(yōu)化前后結構變化圖Fig.3 Structure change diagram of α-Fe-Cr-C cell structure before and after optimization

本文分別在Cr 原子的不同近鄰位置加入C 原子進行結構優(yōu)化，分析結構優(yōu)化前后Cr 原子和C原子之間的距離變化情況，利用式(2)計算它們之間的相互作用能，結果如表5所示。

為了更直觀描述Cr 原子C 原子之間相互作用能隨距離的變化，計算所得的Cr 原子與不同近鄰位置C原子的相互作用能，如圖4所示。由表5及圖4可知，Cr 原子和C 原子之間一直存在排斥力。隨著距離增加，Cr 原子和C 原子之間的相互作用能逐漸降低，且在距離超過1.75aα-Fe時逐漸趨于0，其中，aα-Fe為α-Fe 晶胞的晶格常數(shù)。對于5 近鄰位置和5′近鄰位置構型，兩者相互作用能的不同是因為圖中Cr-C5之間上沒有其他原子，而Cr-C5′之間有1個Fe原子。

圖4 α-Fe中間隙Cr原子與八面體位置C原子間的相互作用能隨距離的變化Fig.4 Interaction energy with distance between interstitial Cr atoms and octahedron position C atoms in α-Fe

表5 結構優(yōu)化前后Cr原子與C原子之間間距的變化及相互作用能Table 5 Interaction and distance change between Cr and C atoms before and after structural optimization

2.3 Cr對C在α-Fe中擴散的影響

α-Fe 中Cr 原子與C 原子的相互作用會影響C原子的擴散，但具體的影響程度還不太清楚。

首先，建立了3×3×3 的Fe54C 模型，利用CI-NEB 方法計算了C 原子在2 個相鄰的八面體位置間的擴散。擴散過程能量變化如圖5所示。從圖5可知：C原子位于八面體位置時，體系能量最低。C原子擴散到四面體位置時，體系能量最高，說明四面體位置是C原子在2個相鄰八面體位置之間擴散過程中的過渡態(tài)位置。計算得到擴散激活能為0.89 eV，與DOMAIN 等[31?32]DFT 模擬結果和實驗結果符合很好，驗證了本文計算擴散激活能方法的準確性。

圖5 3×3×3晶胞計算得到的最小能量路徑Fig.5 Minimum energy path calculated from 3×3×3 unit cells

然后，分別計算C 原子在Fe54C 和Fe53CrC 這2種模型中向其第一、第二近鄰位置擴散的路徑和擴散激活能，結果如圖6所示。從圖6可見：α-Fe(C)體系和α-Fe(C)-Cr體系中C原子在向第一近鄰位置擴散的方式為O-T-O(O代表八面體位置，T代表四面體位置)，向第二近鄰位置擴散方式為O-T-O1-T-O2。2 種體系中C 原子在向第一及第二近鄰位置的擴散方式?jīng)]有發(fā)生變化，說明Cr 原子對C原子的擴散方式?jīng)]有影響。

圖6 C在α-Fe和α-Fe(Cr)中不同擴散路徑的最小能量路徑對比Fig.6 Comparison of minimum energy paths of different diffusion paths in α-Fe and α-Fe(Cr)

最后，研究了擴散激活能隨C-Cr 距離的變化情況，建立4×4×4的α-Fe模型，計算了2種體系中C原子的擴散激活能，α-Fe(C)體系中C原子的擴散激活能為0.93 eV。

α-Fe(C)-Cr中，由圖2可知C原子從Cr的第1～2，2～3，3～4，4～5，5～8，8～9 和9～11 近鄰位置的擴散距離都為0.5aα-Fe，所以計算了這些近鄰位置之間C原子的擴散激活能，計算結果如圖7所示。由圖7可見：隨著距離增加，擴散激活能逐漸變大，然后逐漸趨于穩(wěn)定。通過計算得到，C原子從1～2近鄰位置的擴散激活能為0.88 eV，從9～11近鄰位置的擴散激活能為0.92 eV。C原子離Cr原子越遠，體系的能量越低，這表明當Cr 原子和C 原子離得越遠，結構越穩(wěn)定。

圖7 α-Fe中Cr原子周圍不同近鄰位置間C原子的擴散激活能變化Fig.7 Variation of diffusion activation energy of C atom around Cr atom in α-Fe

3 結論

1)C 原子在α-Fe 內(nèi)占據(jù)八面體位置，Cr 原子在α-Fe中隨機替代Fe原子，Cr原子加入α-Fe(C)中后與Fe 原子和C 原子形成了離子鍵，有較好鍵合作用，增強了晶胞穩(wěn)定性，其中Fe和Cr的成鍵占主要作用。

2)Cr 原子和C 原子在α-Fe 中相互排斥，隨著距離增加其相互作用能越來越弱，體系能量越低，結構越穩(wěn)定，在距離超過約1.75aα-Fe時，相互作用能趨于0。

3)Cr 原子加入α-Fe(C)后，當C 原子朝著遠離Cr 原子的方向擴散時，Cr 原子與C 原子之間的相互作用降低了其擴散激活能，且隨著C-Cr 間距離增加，擴散激活能逐漸趨于α-Fe(C)體系中C 原子的擴散激活能。