Nb基合金中氧原子擴散及團簇化的第一性原理

祝令剛，胡青苗，孫志梅，*

1.北京航空航天大學(xué) 材料科學(xué)與工程學(xué)院，北京 100191

2.中國科學(xué)院金屬研究所，遼寧 沈陽 110016

隨著對發(fā)動機推力的要求越來越高，提高燃燒氣體溫度和燃料利用率勢在必行。當(dāng)前實際應(yīng)用的渦輪葉片為Ni基高溫合金，其使用溫度達到1 150℃，低于其熔點200℃。對于Ni基材料的葉片，如果燃氣溫度進一步提高，需要依賴更好的隔熱涂層和冷卻方法防止葉片過熱甚至融化，然而冷卻葉片意味著能量白白的損耗。另一個選擇是尋找能夠替代Ni基合金的材料。受苛刻服役環(huán)境的限制，潛在的替代材料首先要滿足高熔點，以及高的抗蠕變性能。陶瓷材料、金屬間化合物和難熔金屬是有可能的替代材料。它們各自有優(yōu)缺點：陶瓷和金屬間化合物有比強度高以及抗氧化的優(yōu)點，但其低溫脆性使其難以加工，且較低溫度區(qū)間工作時容易失效；難熔金屬在加工方面有其優(yōu)勢，但最大的問題是自身的抗氧化性能差，遠遠不能滿足高溫合金的服役要求。在難熔金屬中備受關(guān)注的是Nb，熔點為2 469℃［1-2］。

對于純金屬Nb，氧化時環(huán)境中的氧穿過氧化膜向內(nèi)擴散導(dǎo)致持續(xù)氧化。在氧化的初始階段，形成保護性的氧化膜Nb2O5，隨著氧化膜的生長，氧化物在氧化膜-金屬界面生成，在氧化膜中產(chǎn)生內(nèi)應(yīng)力，應(yīng)力的積累可導(dǎo)致氧化膜的開裂，并導(dǎo)致“失穩(wěn)”的線性氧化規(guī)律。例如在600℃ 氧化時，初始階段，氧化膜的生長遵從拋物線規(guī)律，利用聲發(fā)射技術(shù)原位測量時，未發(fā)現(xiàn)氧化膜開裂現(xiàn)象，但15 min后，聲發(fā)射探測到氧化膜發(fā)生了明顯的開裂，氧化動力學(xué)迅速變成了線性規(guī)律，因此純Nb顯然不能在高溫下長期暴露［2-3］。

添加合金元素是提高基體抗氧化性的重要途徑，一般情況下Al、Cr、Si為常用元素，因為它們的氧化物Al2O3、Cr2O3、SiO2具有良好的穩(wěn)定性和對基體的保護性。每種氧化物發(fā)揮保護性的溫度區(qū)間和有效工作環(huán)境是有限定的：對于Al2O3而言，在溫度低于1 300℃ 時生長速度較慢；溫度更高時SiO2的保護性更加有優(yōu)勢，但SiO2在低溫低氧壓下可與Si反應(yīng)生成極具揮發(fā)性的SiO；對于Cr2O3而言，高溫高氧壓下與氧進一步反應(yīng)生成揮發(fā)性的Cr2O3，使氧化膜不斷減薄，除此之外，Cr2O3在恒溫下容易發(fā)生起皺和破裂，主要是由生長應(yīng)力引起的，1 200℃左右金屬Cr表面的Cr2O3即喪失保護性。雖然這3種保護性的氧化物都有缺陷，但是目前仍是設(shè)計耐熱氧化合金時的主要選擇［2］。

為了優(yōu)化氧化膜結(jié)構(gòu)并進一步提高Nb的抗氧化性能，基體中要添加更多的合金元素。20世紀(jì)50年代，Sims等［1］研究了各種二元合金Nb X（X=Ti，V，Zr，Al，Cr，Co，F(xiàn)e，Mn等）的氧化行為，發(fā)現(xiàn)Ti對氧化性能的提高最為明顯，而且當(dāng)Ti的含量達到25at% 時，生成的氧化物也只有Nb2O5；1 300℃ 下，Nb可以溶解8at% 的Al，氧化產(chǎn)物仍為Nb2O5。如果要在Nb合金中生成保護性的氧化膜Al2O3，Al的原子分數(shù)需要達到40%，此時體系中會生成大量的金屬間化合物Nb Al3和Nb Al2等，這些金屬間化合物可以為Al2O3的生長提供Al元素。Murakami等［4］發(fā)現(xiàn)Nb-47Si-20Al合金（主要為 Nb3Si5Al2相）在1 023 K氧化時形成了保護性的Al2O3，大大提高了合金的抗氧化性能。同樣對于Nb-Si合金，體系中的組織多為兩相，一相為Nb基固溶體，另外一相為化合物Nb5Si3，且后者占有相當(dāng)高的體積分數(shù)，通用電氣的Bewlay等［5］在Si原子分數(shù)為16%的多元素Nb合金中發(fā)現(xiàn)Nb5Si3的體積分數(shù)高達46%。對于Nb-Si基合金，若要使其在氧化中生成保護性SiO2，Si的含量要求很高。Cheng等［6］在 Nb-（34at%）Si-（10at%）B 的氧化物中發(fā)現(xiàn)了晶態(tài)但仍然不連續(xù)的SiO2，產(chǎn)物主要還是Nb2O5?？梢妼b合金抗氧化性的提高仍然還有很多問題需要解決［2］。

Nb基合金除了可以作為高溫結(jié)構(gòu)材料，在Snoek型阻尼材料應(yīng)用方面也有些報道［7-9］。在外界沖擊等作用力下，阻尼材料中的間隙原子如氧原子等，發(fā)生跳躍遷移，從而消耗外力對基體的影響。實際中可通過添加合金元素，改變間隙原子在基體中的擴散激活能，從而達到改善阻尼合金性能的目的。

從原子尺度上理解上述兩種應(yīng)用背景下Nb合金的材料科學(xué)問題是十分復(fù)雜的，尤其是氧化問題，涉及氧原子的擴散和聚集、氧化物的形核長大、多層氧化物以及氧化物與基體之間的結(jié)合等過程?；诿芏确汉碚摰牡谝恍栽碛嬎汶m然在時間和空間尺度上均不能完整地描述合金的氧化或阻尼過程，但對復(fù)雜過程中涉及的細節(jié)問題是可以獨立出來精確處理的，比如氧原子在金屬基體中的擴散問題［10］、氧化初期基體表面氧化物前驅(qū)體的結(jié)構(gòu)［11-13］、基體與氧化物界面的結(jié)構(gòu)和強度等［14-15］。這些問題的第一性原理計算研究，除了為實驗現(xiàn)象提供物理圖像外，對更接近實際過程的大尺度計算模擬也提供著基本數(shù)據(jù)的支持，同樣對計算模擬輔助材料設(shè)計有重要的意義。

基于上述考慮，本文利用第一性原理計算首先研究了氧原子在基體Nb中的擴散問題，作為潛在的發(fā)動機葉片材料，Nb基合金服役環(huán)境中不斷受到應(yīng)力作用，考慮了單軸應(yīng)變下氧原子擴散行為的變化；研究了氧原子在Nb中形成原子對的傾向性，以及合金元素Ti和Zr對其的影響。

1 計算參數(shù)設(shè)置

采用奧地利維也納大學(xué)開發(fā)的第一性原理VASP（Vienna Ab-initio Simulation Package）軟件包［16-17］進行計算。電子間的交換關(guān)聯(lián)泛函采用廣義梯度近似PW91-GGA［18］。電子與原子核的相互作用通過PAW贗勢［19］描述。為了保證計算的精確性，這里對于涉及的過渡金屬原子Nb、Ti和Zr，其內(nèi)層p電子也被處理成價電子。平面波截斷能設(shè)置為300 eV，k點網(wǎng)格為4×4×4。計算過程中費米能級附近展寬采用Methfessel-Paxton方法，寬度為0.2 eV。對于純Nb，計算得到的晶格常數(shù)為3.321?，與實驗數(shù)據(jù)的3.307 ?［20］符合較好。模擬合金用的超晶胞為3×3×3，包含54個原子，考慮合金結(jié)構(gòu)時，其中1個Nb原子被置換，這樣所模擬的合金的濃度為1.85%。

對于氧原子在Nb中的擴散，采用NEB方法（Nudged Elastic Band Method）［21］得到擴散路徑和勢壘。使用NEB方法時要先得到擴散前后對應(yīng)的2個結(jié)構(gòu)，即初始態(tài)和終了態(tài)；在2個態(tài)之間線性插值得到一些中間態(tài)（Image），并建立起初始的擴散路徑；在路徑優(yōu)化過程中，平行于路徑方向上人為施加彈簧力（Spring Force），以保證中間態(tài)不會弛豫到初始或終了結(jié)構(gòu)；優(yōu)化完成后要使各中間態(tài)在垂直于路徑方向的分力為0。本文計算過程中選取了4個中間態(tài)，NEB計算收斂標(biāo)準(zhǔn)為各中間態(tài)原子上的力小于0.05 e V/?。這里考慮了在［001］方向施加單向應(yīng)變時氧原子擴散行為的變化，為得到施加應(yīng)變后超晶胞的結(jié)構(gòu)，計算過程中固定單方向的晶格常數(shù)，其他兩個方向自由弛豫。

2 計算結(jié)果與討論

2.1 單向應(yīng)變下氧原子的擴散

本部分計算了氧原子在施加單向應(yīng)變的體心立方Nb中的擴散。如圖1所示，應(yīng)變?yōu)?時（S=0），八面體間隙1、2和3是等同的，此時氧原子在任意相鄰八面體間隙之間的擴散均是等價的，勢壘大小為0.92 eV，與實驗值1.11～1.14 e V［22］符合較好；而且擴散路徑為直線，過渡態(tài)結(jié)構(gòu)為氧原子占據(jù)兩八面體間隙之間的四面體間隙位置，即擴散勢壘的大小可以理解為氧原子占據(jù)八面體間隙和四面體間隙時對應(yīng)體系的能量差［23］。施加的單向應(yīng)變大小為0.02、0.04、0.06、0.08和0.10，方向如圖1中箭頭所示。施加單向應(yīng)變后，八面體間隙均會發(fā)生畸變，但間隙1和間隙2仍是等價的，計算發(fā)現(xiàn)氧原子在這2個間隙之間擴散的勢壘隨著應(yīng)變的增大而不斷減小，如圖2（a）所示，也就是說施加應(yīng)變后氧原子占據(jù)八面體間隙和四面體間隙時對應(yīng)體系的能量差不斷減小。反應(yīng)坐標(biāo)（橫坐標(biāo)）顯示出2個間隙之間的距離隨著應(yīng)變的施加而逐漸被拉開。

圖1 單向應(yīng)變下氧原子在Nb中的擴散路徑Fig.1 Pathways for oxygen atom diffusion in uniaxialstrained Nb

圖2 不同單向拉伸應(yīng)變量下氧原子擴散的勢壘曲線Fig.2 Barrier curve of oxygen atom diffusion under different uniaxial tensile strains

間隙2和間隙3在應(yīng)變作用下不再等價。對于間隙2，扁八面體的2個頂點距離進一步被壓縮，而間隙3相應(yīng)的2個頂點間距不斷拉長?？梢韵胂箝g隙3更適合氧原子的占據(jù)，這從圖2（b）中的擴散勢壘曲線可以直接反映出，擴散前后體系的能量變化隨著應(yīng)變的增大而變大，氧原子更傾向占據(jù)間隙3。間隙2到間隙3的擴散勢壘隨著應(yīng)變的增大而變小，即從動力學(xué)擴散的角度講，氧原子也傾向于向間隙3偏聚。

從圖2中2個擴散勢壘曲線可以看出，當(dāng)施加單向應(yīng)變后，氧原子會以比無應(yīng)變下更快的擴散速度（勢壘更?。┌l(fā)生偏聚，而且是沿施加應(yīng)變的方向排列。顯然氧原子的這種偏聚會導(dǎo)致應(yīng)力的集中，容易引起裂紋萌生，降低材料的力學(xué)性能。

2.2 Nb合金中氧原子對的形成能力

此部分討論了2個氧原子在Nb中占位較近并形成原子對的可能性，圖3給出了可能形成氧原子對的全部構(gòu)型。圖3中字母代表各種八面體間隙位置，這里只考慮了合金元素X的最近鄰和次近鄰位置，因為研究發(fā)現(xiàn)［9］，第3近鄰位置時合金元素與氧原子的相互作用已經(jīng)可以忽略不計。

圖3 合金元素X周圍第1和第2近鄰八面體間隙位置Fig.3 First-and second-nearest octahedral interstice near alloying element X

本文中氧原子對的構(gòu)型用圖3中的2個字母表示，比如原子對AH表示一個氧原子位于A間隙，另一個氧原子位于H間隙的構(gòu)型。這里計算氧原子形成原子對能力的方法具體為：將氧原子形成原子對的過程分成2個階段，第1階段（Stage 1）對應(yīng)氧原子濃度很低的情況，此時氧原子都有機會在最強的吸引位置（對應(yīng)體系的能量最低），即每個氧原子都可以單獨“享有”一個合金原子，此時對應(yīng)的形成氧原子對的能力可以用能量表示為

式中：EXOO為含合金元素和氧原子對的超胞的能量；EX為僅包含合金原子的超胞的能量；EXO為含合金元素以及單個氧原子的超胞的能量且此氧原子處于最強吸引間隙位置。第2階段（Stage 2），對應(yīng)氧濃度稍微高些的情況，此時假設(shè)每個合金原子周圍已經(jīng)有一個氧原子占據(jù)在能量最優(yōu)的位置，此時氧原子對的能力可用能量表示為

式中：E為純Nb超胞的能量；ENbO為Nb基體中含一個氧原子時超胞的能量，可見ΔE2也表示把氧原子從基體內(nèi)部移動到合金原子周圍與另外的氧形成原子對時體系能量的變化。

除了純Nb體系，這里考慮了2種合金元素Ti和Zr，之所以選擇這2種合金元素是因為前期研究發(fā)現(xiàn)Ti對第1近鄰間隙位置的氧原子作用最強，而Zr對位于其第2近鄰間隙的氧原子作用最強，V與Ti的情況類似，而Hf與Zr的情況類似，即這里選擇的Ti和Zr在合金元素與氧原子相互作用方面是具有代表性的［9］。Stage 1和Stage 2對應(yīng)的相互作用能計算數(shù)值如表1所示，對于Stage 1，能量均為正值，表示氧原子的濃度很低時，更傾向于單獨占據(jù)在合金原子附近的最低能量對應(yīng)的間隙位置，不會形成原子對。對于Stage 2，Ti和Zr都出現(xiàn)了負值能量，表示此時氧原子會從基體移動到合金原子周圍形成原子對。還可看出對于Zr，氧原子對形成的可能性更大，因其負值能量較多且絕對值較大。因此可以判定，對于合金原子為Zr時，氧原子對的形成可能性比Ti對應(yīng)的情況大得多。同樣可以推斷X=Hf時，情況和Zr類似。

可見，加入不同合金元素之后，會對基體中氧原子的分布產(chǎn)生不同的影響，而氧原子的分布對應(yīng)力集中，甚至是氧化物的形核都有重要作用。當(dāng)然這里只計算了2個氧原子形成氧原子對的可能性，要研究對材料氧化等宏觀性能的直接影響，還需進一步考慮更多的氧原子形成團簇的情況。

表1 兩階段中氧形成不同構(gòu)型原子對的能力Table 1 Pairing ability of oxygen atoms into certain configurations at two stages

為進一步理解合金元素周圍氧原子對形成能力不同的物理層次原因，圖4（a）～圖4（c）給出了不同合金成分下的電子態(tài)密度圖，插圖對應(yīng)的能量范圍為-8.5～-6.5 eV。這里以AH和BG 2種構(gòu)型為例，因為此構(gòu)型下氧原子分別在Ti和Zr周圍配對的傾向性最大，即此2種構(gòu)型具有代表性。

圖4 氧原子對構(gòu)型為AH和BG時對應(yīng)超胞的態(tài)密度Fig.4 Density of states for superstructure with oxygen atom pair having configuration of AH and BG

從圖4可以看出氧原子對組成AH和BG 2種構(gòu)型時體系的態(tài)密度圖差別很小，主要因為對應(yīng)超胞的總能量差別較小。相對明顯的差異在能級為-7.5 e V左右的態(tài)密度峰（圖4中的插圖），對于X=Nb和Ti，構(gòu)型BG對應(yīng)的此峰向高能級偏移明顯，處于較高的能量水平，即BG對應(yīng)構(gòu)型的能量也較高，這與表1中列出的能量數(shù)據(jù)吻合。對于X=Zr時，-7.5 e V左右的峰出現(xiàn)了明顯的分裂，BG對應(yīng)峰的左側(cè)部分仍出現(xiàn)了高能級的偏移（相對AH構(gòu)型），但其偏移量較小，而且其右側(cè)的分裂出來的峰出現(xiàn)了向低能級移動的趨勢，可以看出這是導(dǎo)致X=Zr時BG對應(yīng)的構(gòu)型能量較低，而X=Nb和Ti時則不然的一個原因。

3 結(jié) 論

1）對Nb施加單向應(yīng)變，即產(chǎn)生了被拉長的八面體間隙（在應(yīng)變施加的方向被拉長）和被進一步壓扁的八面體間隙，計算發(fā)現(xiàn)氧原子更傾向于占據(jù)被拉長的間隙中。當(dāng)應(yīng)變?yōu)?.1時，占位優(yōu)勢（能量差）達到了1.0 eV，而且向這種類型間隙擴散的速度會加快（對應(yīng)的勢壘減小）。應(yīng)變最終導(dǎo)致氧原子沿應(yīng)變方向偏聚排列，且隨著應(yīng)變的增大這種現(xiàn)象更加明顯。

2）在純Nb中，氧原子濃度不高時，傾向于均勻分布在基體中，無形成原子對的傾向。在元素Ti和Zr摻雜后，由于2種元素對氧原子均有吸引作用，氧原子的分布情況會類似于合金元素的分布，通常也是均勻分布在基體中的；當(dāng)氧原子的濃度升高后，氧原子傾向于在合金元素周圍形成原子對，原子對的形成可使體系的能量降低最多達0.29 e V。計算結(jié)果也表明，Zr比Ti更容易使得氧原子形成原子對。電子結(jié)構(gòu)分析表明，態(tài)密度峰的分裂和偏移是導(dǎo)致形成不同構(gòu)型原子對的重要原因。

參 考 文 獻

［1］ Ghosh G，Olson G B.Integrated design of Nb-based superalloys：Ab initio calculations，computational thermodynamics and kinetics，and experimental results［J］.Acta Materialia，2007，55（10）：3281-3303.

［2］ Zhu L G.Theoretical study on the oxidation of Ti and Nbbased alloys［D］.Shenyang：Insititute of Metal Research，Chinese Academy of Sciences，2014.（in Chinese）祝令剛.Ti和Nb合金中氧化問題的理論研究［D］.沈陽：中國科學(xué)院金屬研究所，2014.

［3］ Birks N，Meier G H，Pettit F S.Introduction to the hightemperature oxidation of metals［M］.Xin L，Wang W，translated.Beijing：High Education Press，2010：84-85.（in Chinese）Birks N，Meier G H，Pettit F S.金屬高溫氧化導(dǎo)論［M］.辛麗，王文，譯.北京：高等教育出版社，2010：84-85.

［4］ Murakami T，Sasaki S，Ichikawa K，et al.Microstructure，mechanical properties and oxidation behavior of Nb-Si-Al and Nb-Si-N powder compacts prepared by spark plasma sintering［J］.Intermetallics，2001，9（7）：621-627.

［5］ Bewlay B P，Jackson M R，Zhao J C，et al.A review of very-high-temperature Nb-silicide-based composites［J］.Metallurgical and Materials Transactions A：Physical Metallurgy and Materials Science，2003，34 （10）：2043-2052.

［6］ Cheng J C，Yi S，Park J S.Oxidation behaviors of Nb-Si-B ternary alloys at 1100°C under ambient atmosphere［J］.Intermetallics，2012，23：12-19.

［7］ Yin F X，Yu L M，Ping D H.Modeling and control of the high damping behavior in Ti-Nb-O alloys［J］.Materials Science and Engineering：A，2009，521-522：372-375.

［8］ Yu L M，Yin F X，Ping D H.Natural mechanism of the broadened Snoek relaxation profile in ternary body-centered-cubic alloys［J］.Physical Review B，2007，75（17）：174105-174112.

［9］ Yu L M，Yin F X.Internal friction of Niobium-Titanium-Oxygen alloys［J］.Journal of Materials Science，2007，42（18）：7819-7826.

［10］ Zhu L G，Wang H，Hu Q M，et al.Trapping of interstitial defects：filling the gap between the experimental measurements and DFT calculations［J］.Journal of Physics：Condensed Matter，2013，25（43）：435402-435408.

［11］ Wang L，Shang J X，Wang F H，et al.Oxygen adsorption onγ-TiAl surfaces and the related surface phase diagrams：A density-functional theory study［J］.Acta Materialia，2013，61（5）：1726-1738.

［12］ Wang Q G，Shang J X，Yang Z.First-principles study on the initial oxidization of a Nb（100）surface［J］.The Journal of Physical Chemistry C， 2012， 116 （44）：23371-23376.

［13］ Todorova M，Reuter K，Scheffler M.Density-functional theory study of the initial oxygen incorporation in Pd（111）［J］.Physical Review B，2005，71（19）：195403-195410.

［14］ Zhu L G，Hu Q M，Yang R，et al.Binding of an oxide layer to a metal：The case of Ti（10-10）/TiO2（100）［J］.The Journal of Physical Chemistry C，2012，116（6）：4224-4233.

［15］ Zhu L G，Hu Q M，Yang R，et al.Atomic-scale modeling of the dynamics of titanium oxidation［J］.The Journal of Physical Chemistry C，2012，116（45）：24201-24205.

［16］ Kresse G，F(xiàn)urthmuller J.Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set［J］.Computational Materials Science，1996，6（1）：15-50.

［17］ Kresse G，F(xiàn)urthmüller J.Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set［J］.Physical Review B，1996，54（16）：11169-11186.

［18］ Perdew J P，Wang Y.Pair-distribution function and its coupling-constant average for the spin-polarized electron gas［J］.Physical Review B，1992，46（20）：12947-12954.

［19］ Bl?chl P E.Projector augmented-wave method［J］.Physical Review B，1994，50（24）：17953-17979.

［20］ Ikehata H，Nagasako N，F(xiàn)uruta T，et al.First-principles calculations for development of low elastic modulus Ti alloys［J］.Physical Review B，2004，70（17）：174113-174120.

［21］ Henkelman G，Uberuaga B P，Jonsson H.A climbing image nudged elastic band method for finding saddle points and minimum energy paths［J］.The Journal of Chemical Physics，2000，113（22）：9901-9904.

［22］ Vykhodets V B，Kurennykh T E，Lakhtin A S，et al.Activation energy of hydrogen，oxygen，and nitrogen diffusion in metals［J］.Doklady Physical Chemistry，2005，401（4）：56-58.

［23］ Li C X，Luo H B，Hu Q M，et al.Theoretical investigations of interstitial atoms in bcc metals：Local lattice distortion and diffusion barrier［J］.Computational Materials Science，2012，58：67-70.