Ti, Nb, Al及其二元合金狀態方程的第一原理計算及其應用

程 超， 馬云莉， 曹超銘， 孫嘉興， 劉艷俠

(遼寧大學物理學院，沈陽 110036)

1 引 言

如果金屬狀態方程可被求解，那么金屬的許多能量機制就可以被了解，金屬狀態方程可用于求解金屬性質的諸多領域，其中一個重要應用用于擬合或驗證原子間相互作用勢. 原子間相互作用勢是有關原子尺度計算機模擬的基礎，它的精確與否將直接的影響到模擬的準確性. 如此重要的同時，原子間相互作用勢的構建也是方法眾多，并且涵蓋材料物理性質也不盡相同. 但如此眾多的不同中存在著相同之處，就是都要滿足體系能量隨原子間距離的變化關系 (E-r關系).

然而，金屬狀態方程一直以來是難以確定的，直到1984年，Rose等人[1]針對金屬狀態方程進行研究，并提出了金屬狀態方程的普適形式. 時至今日，Rose的狀態方程在原子間相互作用勢等計算物理領域有重要的應用. 本文將Rose狀態方程可視化，畫出E-r曲線作為參考，對比第一原理計算結果，探究計算E-r曲線計算的最優方法. 之所以采用第一原理計算是由于，早年研究[2-5]表明，在擬合原子間相互作用勢時， 如果結合第一原理數據可以顯著提高勢函數的精確性. 例如，Baskes等人[4]通過擬合勢函數和第一原理研究了Al原子間相互作用力的范圍. 他們發現，在擬合過程中，使用第一原理數據擬合的勢比使用實驗數據擬合更準確. 但使用第一原理探究金屬及其合金E-r關系卻鮮有報道.

本文通過使用基于密度泛函理論的VASP軟件包進行第一原理計算[6, 7]，計算并探究了航空航天常用材料Ti及合金化元素Nb和Al在使用不同贗勢計算E-r曲線的差異性，研究了由三種金屬分別構成的二元合金β相的Ti-25Nb(at.%)合金， L10結構的TiAl合金以及DO22結構的NbAl3合金的第一原理計算精確度. 并且我們以Ti-Nb系合金為例，分析了合金成分不同時，其E-r關系滿足的一般規律. 根據材料的E-r關系，使用最小二乘法擬合了材料的體彈性模量. 根據合金曲線表現出的規律，提出合金E-r關系的新的計算方法，并且發現使用該方法得到的合金的E-r關系數據擬合體彈性模量的精確度非常高. 本文工作旨在為原子間相互作用勢的構建及驗證提供理論數據支持.

2 計算方法和模型

2.1 計算方法

本文第一原理計算采用基于密度泛函理論[8](DFT)的VASP軟件包[9, 10](Vienna Ab initio Simulation Package)進行的. 選擇廣義梯度近似GGA[11, 12](Generalized gradient approximation)以及局域密度近似LDA[13, 14](Local-density approximation)來處理電子的交換關聯能. 自洽精度取1×10-5eV/atom，通過采用Monkhorst-Pack方法[15]選取布里淵區k空間網格點來進行體系總能量的布里淵區積分計算，且計算考慮了自旋極化. 自洽迭代循環中，總能的收斂標準選取 1×10-4， 原子弛豫收斂標準取 0.01，原子弛豫最大步數設置為 100，并采用共軛梯度算法來優化原子的位置. 對于2×2×2超胞的Ti-Nb合金，K-points 取 6×6×6，對于1×2×1超胞的TiAl合金以及NbAl3合金K-points 取 6×4×6.

對所選體系進行結構優化，再取結構優化后的模型，縮放其晶胞尺寸計算能量. 對于作為對比分析的Rose曲線而言，我們根據Rose的普適性狀態方程，并用數學計算軟件Mathematica計算數據，計算選取的對應金屬實驗值在本文表2中列出.

2.2 計算模型

如圖1所示，本文計算模型結構均為常溫下穩定結構. 純金屬為單胞的HCP-Ti，BCC-Nb和FCC-Al，合金為2×2×2超胞的β相Ti-25Nb(at.%)合金，1×2×1超胞L10結構的TiAl合金以及1×2×1超胞DO22結構的NbAl3合金. 圖1給出的模型結構(010)面向上，右側為(100)面.

圖1 合金結構 (a) Ti-25Nb(at.%); (b) TiAl; (c) NbAl3Fig.1 Structures of alloy (a) Ti-25Nb (at.%), (b) TiAl, (c) NbAl3

3 結果與討論

3.1 Ti,Nb和Al的能量以及 E-r曲線分析

體系結合能計算公式為：

(1)

公式中Etot為體系總能量，E1與E2為體系中不同元素的孤立原子能量，n和m為對應元素的原子個數. 表1為Ti,Nb和Al使用不同贗勢計算出的結合能并與實驗值進行對比. 如表1所示，三種金屬計算結果表明，選用GGA方法處理的結果明顯優于LDA方法，而相比較之下同樣是廣義梯度近似，使用GGA-PBE[16]贗勢計算金屬Ti和Al的結果比GGA-PW91[17]贗勢的計算結果更加的接近實驗值. 而對于金屬Nb而言，使用GGA-PW91贗勢計算的結果與實驗值更接近.

表1 不同贗勢計算的結合能

Table 1 The results of cohesive energy calculated by different pseudopotentials

ElementEPW91(eV)EPBE(eV)ELDA(eV)Experimental(eV)Ti5.33634.93975.4844.85[18]Nb7.48136.49147.98317.57[18]Al3.37293.40783.90673.39[18]

金屬的狀態方程對于了解金屬能量學等方面至關重要. Rose就金屬的狀態方程做過研究并且給出了普適的狀態方程，忽略高次項為：

表2 Rose曲線計算所用數據

圖2 不同方法計算的E-r曲線(a) Ti, (b) Nb, (c) Al Fig.2 The comparison of E-r curves calculated by different methods (a) Ti, (b) Nb, (c) Al

如圖2所示，第一原理計算和Rose狀態方程所得到的E-r關系對比分析，發現，第一原理使用不同贗勢計算得出的曲線大體上近似，卻有著不可忽視的差別. 所有的探究都表現出，平衡態右側部分的第一原理計算結果與Rose曲線稍有偏差，而平衡態左側部分，金屬Al和Nb的第一原理計算結果與Rose曲線吻合得很好，但金屬Ti卻相差較大. 如圖2(a)中，對比發現無論在高能量的排斥區還是低能量的吸引區域，GGA-PBE贗勢都要比另外兩種方法計算的結果更加的接近Rose曲線，并且發現在平衡態處計算的能量與Rose曲線最接近. 這都展現出，對于金屬Ti而言PBE贗勢計算的更加準確. 在圖2(b)中，GGA-PW91贗勢的計算結果更接近于Rose曲線，同時結合表1的計算結果也表明，采用GGA-PW91贗勢計算Nb更為精確. 如圖2(c)，可以看出，第一原理計算的結果與Rose曲線吻合度很高，而且PW91計算結果與PBE計算結果幾乎重合，如圖2(c),可以看到PW91贗勢計算結果與PBE贗勢計算結果的細小差別，難以分辨哪個結果與Rose更接近，但是從結合能計算的結果來看，如表1所示，使用PW91贗勢計算結合能為3.3729 eV，而PBE贗勢為3.4078 eV，相比較而言，PBE贗勢計算結果更加的精確.

3.2 合金體系的E-r曲線分析

如圖3，分別為β相Ti-25Nb(at.%)合金，L10結構TiAl合金以及DO22結構NbAl3合金E-r關系的第一原理計算結果. 對于每種合金的合金化元素我們都選取上述工作中最優的計算結果，比如，金屬Ti和Al選用GGA-PBE的計算結果，而金屬Nb則選用GGA-PW91計算結果. 從圖3可以看到，合金的E-r曲線分布處于其合金化元素E-r曲線之間，這與已有的研究結果一致[23]，可見計算結果吻合度不錯.

為了探究和合金化元素成分不同對合金E-r曲線的影響，我們又選取了Ti-37.5Nb(at.%)合金[24]以及Ti-50Nb(at.%)合金[21]進行了對比計算，如圖3(e)所示. 可以看到，合金的成分雖然改變但曲線仍然處在其合金化元素曲線之間，且隨著Nb元素含量的提高，合金曲線越趨近于Nb的曲線，Nb含量較少的Ti-25Nb(at.%)合金曲線更接近于Ti的曲線. 這也就表明同種合金系合金的E-r關系，處于其合金化元素之間，不會超出或低于該范圍；并且合金的曲線會更靠近合金內含量較高的元素的E-r曲線.

圖3 二元合金體系的E-r曲線 (a) β相Ti-25Nb(at.%)合金; (b) L10結構TiAl合金; (c) DO22結構NbAl3合金 (d) Ti-Nb系合金對比Fig. 3 Comparison of E-r curves for the binary alloys: (a) β-phase Ti-25Nb (at.%) alloy, (b) L10 structure TiAl alloy, (c) DO22 structure NbAl3 alloy and (d) Ti-Nb alloy

根據上述表現出來合金E-r曲線與合金中元素含量的關系，我們猜想能否使用純金屬的E-r關系，通過合金中對應的比例轉化數據直接得到合金的E-r關系. 以Ti-25Nb(at.%)合金為例，利用純金屬Ti和Nb的E-r關系，按照合金中對應元素的含量進行轉化，得到的結果如圖4所示. 從圖中可以看到，使用第一原理計算得到的曲線(虛線)與使用純金屬E-r關系轉化得到的曲線(實線)吻合的非常好，特別是在平衡態附近曲線幾乎重合. 這表明這種方法是可行的，這也給為計算合金體系E-r關系提供了一種簡單的方法，即利用第一原理計算純金屬的E-r關系，再通過合金中元素的含量，將純金屬E-r關系比例轉化為合金的E-r關系.

圖4 使用不同方法得到的Ti-25Nb(at.%)合金的E-r曲線Fig. 4 The E-r curves of Ti-25Nb (at.%) Alloys obtained using different methods

3.3 態方程的應用——體彈性模量的計算

彈性模量是描寫固體材料彈性形變表征力學性質的重要物理量，是選定機械構件材料的重要依據[25, 26]，也是用于擬合原子間相互作用勢的重要參數. 采用本文中第一原理計算結果，通過最小二乘法擬 合曲線的方式計算體彈性模量. 選取的擬合函數為[27]：

(3)

式中Etot，E0，V，V0以及α分別為總能，平衡態能量，原子體積，平衡態原子體積和擬合參數. 體彈性模量由以下式子給出：

(4)

式中結合能E可由總能Etot與平衡態能量E0差得到，所以平衡態下：

(5)

從(5)式可以看到，擬合得到參數α即可得到體彈性模量. 上述方程使用的是能量體積關系，所以我們將上述E-r關系數據換算為E-V關系，使用數學軟件包Mathematica編寫程序擬合參數. 擬合過程中V0和E0選取實驗值列于表2和表4中，擬合得到結果列于表3中.

表3 體彈性模量B0的計算結果

從體模量的擬合結果與實驗值對比來看，我們發現擬合的結果與實驗值相差10 Gpa左右，誤差較小，擬合的結果很好. 長期以來實驗數據一直是各類研究的參考依據，但是多種因素也影響著測定的準確性，所以理論計算材料物理性質是十分必要的. 從表3的擬合體模量結果的精確度來看，一方面檢驗了本工作中E-r關系計算的準確，另一方面，也表明使用這種方法計算材料體模量的可行性.

在文章中3.2節中提到，可以使用純金屬E-r關系按加權平均轉換為合金的E-r關系. 為了進一步驗證該方法的可行性，我們使用轉化后的合金數據，最小二乘法擬合合金體模量，擬合所用到的實驗值與擬合結果列于表4中. 對比表4與表3中合金體模量擬合結果可以看出，表4中的結果更接近實驗值，擬合結果更好. 這說明使用純金屬E-r關系，按加權平均轉換為合金的E-r關系的方法是可行的，并且使用該方法得到的E-r數據擬合體模量結果更佳.

表4 合金實驗數據及其體模量擬合結果

Table 4 The experimental data and the fitting results about bulk modulus of alloy

Alloya(?)c(?)V0(?3)E0(eV)B(GPa)Ti-25Nb3.273[28]—17.531055.64[32]115.503TiAl3.977[33]4.0565[33]16.04014.51[34]107.733NbAl33.80[35]8.75[35]15.793754.97[36]123.574

4 結 論

本文采用第一原理計算了純金屬Ti,Nb和Al及其二元合金的E-r曲線，計算并討論了使用不同贗勢對于不同的金屬適用程度，及合金化元素的含量對該合金E-r關系的影響. 研究結果表明：

1)對于HCP結構的金屬Ti和FCC結構的Al而言，使用GGA-PBE贗勢計算其平衡態結合能的結果較為優異，比其他贗勢計算結果更接近實驗值. 而對于BCC結構的金屬Nb而言，使用GGA-PW91贗勢計算的結果與實驗值更加接近.

2)純金屬E-r關系計算結果與結合能計算結果相吻合，相比之下使用GGA-PBE贗勢計算Ti和Al的E-r關系與Rose的結果更加吻合，而對于金屬Nb則使用GGA-PW91贗勢計算得到的結果更佳.

3)合金E-r關系計算結果較好，曲線均處于其合金化元素曲線之間，并且，曲線會更靠近于含量較多的合金化元素的E-r曲線.

4)使用計算得到的E-r數據采用最小二乘法擬合得到的純金屬和二元合金的體彈性模量與實驗值相比誤差均較小，結果都很好.

5)可以使用純金屬的E-r關系，根據加權平均的方法計算合金的E-r關系，使用該方法得到的合金的E-r關系計算的體彈性模量與實驗值非常接近.