Ti2AlNb 合金的嵌入原子勢*

劉杰 劉艷俠

(遼寧大學(xué)物理學(xué)院,沈陽 110036)

分子動力學(xué)模擬是一種行之有效的計算機模擬方法;然而,由于缺少合適的多元合金原子間勢,因而限制了分子動力學(xué)模擬的應(yīng)用.多元合金原子間勢的開發(fā)一直具有挑戰(zhàn)性.本文在嵌入原子勢模型的框架下,提出一種適用于三元有序合金的原子間勢構(gòu)建方法,并開發(fā)了適用于原子尺度力學(xué)行為模擬的Ti2AlNb 合金新型原子間勢.該勢能夠很好地再現(xiàn)B2-Ti2AlNb 的彈性常數(shù)、未弛豫的空位形成能、置換原子形成能、換位原子形成能、表面能和三種有序構(gòu)型(B2 相、D019 相、O 相)在不同體積下的內(nèi)聚能.為了進一步檢驗勢函數(shù),計算了B2 相的E-V 曲線,結(jié)果與Rose 曲線符合得很好;利用分子動力學(xué)模擬研究了B2 相的熔化轉(zhuǎn)變過程,結(jié)果大致反映了實驗情況.本文的工作一方面為開發(fā)多元合金原子間勢提供一種途徑,另一方面為模擬計算Ti2AlNb 合金的工作者提供一種選擇.

1 引言

目前開發(fā)原子間勢主要有兩類方法: 一類是近幾年發(fā)展起來的機器學(xué)習(xí)法[1];另一類是傳統(tǒng)的擬合原子間勢參數(shù)的方法.

隨著人工智能的快速發(fā)展,機器學(xué)習(xí)法開發(fā)原子間勢應(yīng)用而生,已經(jīng)有許多成功應(yīng)用的案例[2,3],Smith等[4]在2020 年用機器學(xué)習(xí)的方法構(gòu)建了鋁元素的原子間勢,他們應(yīng)用了多種神經(jīng)網(wǎng)絡(luò)來訓(xùn)練大量數(shù)據(jù),比如由金屬的不同結(jié)構(gòu)的多種性質(zhì),最后得到了計算精度非常高的原子間勢.機器學(xué)習(xí)法構(gòu)建原子間勢的優(yōu)勢在于減少了人為干預(yù)的成分,規(guī)避了傳統(tǒng)方法建勢的參數(shù)化問題,且可以達到與第一原理方法相媲美的精度[5].但是機器學(xué)習(xí)勢的非物理形式阻礙了對模型參數(shù)的直接解釋,并阻止了對參考數(shù)據(jù)集沒有很好描述的構(gòu)型空間區(qū)域的外推,同時機器學(xué)習(xí)勢需要準備大量準確的訓(xùn)練數(shù)據(jù)集.

傳統(tǒng)方法構(gòu)建原子間勢也同樣具有一定的優(yōu)勢,首先構(gòu)建的勢函數(shù)形式簡單,與分子動力學(xué)接口友好,需要準備的擬合數(shù)據(jù)數(shù)量遠少于機器學(xué)習(xí)勢,加之人工特有的跨學(xué)科優(yōu)勢,擬合過程中可根據(jù)物理理論進行人為干預(yù),可以更加靈活地調(diào)整參數(shù),且擬合過程中的每一步都具有清晰的物理意義.

Ti2AlNb 合金具有很多優(yōu)異性能,如高強度、高比模量以及較高的抗氧化性等[6-9],使其成為航空航天發(fā)動機的首選材料之一.關(guān)于Ti2AlNb 合金的實驗研究,可采用摻雜、多種強化機制及熱處理等工藝來改善合金的使用性能[10].理論研究主要采用第一原理計算的方法.如Pathak 和Singh[11]用第一原理方法計算了三種有序構(gòu)型(B2 相、D019相、O 相)結(jié)構(gòu)的穩(wěn)定性,并給出了各自的晶體結(jié)構(gòu)、空間群、原子坐標等信息.

由于缺少合適的Ti2AlNb 合金原子間勢,因而限制了使用分子動力學(xué)模擬的方法研究Ti2AlNb合金的性能.Ti2AlNb 基合金一般由α2,O,B2 中的兩相或三相構(gòu)成.β/B2 都是bcc 結(jié)構(gòu),β是無序相,B2 是有序相.Nb是β相的穩(wěn)定元素,由于Ti2AlNb 基合金中Nb 含量比較高,所以β相穩(wěn)定元素一直存在,且B2 相滑移系較多,是塑性相,因此,本文采用傳統(tǒng)方法主要針對B2 相開發(fā)Ti2AlNb合金原子間勢,期望能為分子動力學(xué)模擬工作者提供一種選擇.

2 原子間勢的構(gòu)建

2.1 模型

本文選用嵌入原子 (embedded-atom method,EAM)勢模型[12],在該模型下體系的總能為

其中α,β,γ,?,σ,μ為擬合參數(shù);rc為截斷距離;h為截斷函數(shù)形狀的調(diào)整參數(shù),截斷函數(shù)形式為

對于本文的三元合金原子間勢的構(gòu)建,需要擬合三種交叉對勢,分別為?(rTiAl),?(rTiNb)和?(rNbAl),包含18 個擬合參數(shù),再加上sx和gx(下標x表示Ti,Al,Nb 的任意一種元素),共24 個擬合參數(shù),截斷距離選用相同的值為rc.

在擬合合金的原子間勢時,應(yīng)用了關(guān)于以下規(guī)范變換的不變性[14,15]:

x和y代表Ti,Al,Nb 三種元素中的任意一種,上述變換中sx和gx是任意常數(shù),對于三元有序構(gòu)型比如B2-Ti2AlNb 合金,不是嚴格意義滿足能量對其公式變換的不變性(本質(zhì)上這種變換關(guān)系對于合金需要滿足高度有序的對稱性,規(guī)范變換會改變勢函數(shù)的形狀).對于A-B-C 體系,可以使變換后的多體勢滿足[14],這樣多體勢函數(shù)就可以在平衡狀態(tài)下對應(yīng)原子的電子密度處取得極值(通常情況下為極小值)[16],本文在B2 相的環(huán)境中使其滿足的值由歸一化條件確定,此外本文盡量使和的值接近1.

2.2 擬合數(shù)據(jù)庫

關(guān)于Ti2AlNb 合金性能的實驗值非常有限,因此本文主要使用第一原理計算的數(shù)據(jù)來參數(shù)化原子間勢.使用基于密度泛函理論的VASP 軟件(Viennaab-initiosimulation package),交換關(guān)聯(lián)泛函選擇Perdew-Wang 91 形式,平面波截斷能為500 eV,總能收斂標準為10—5,原子弛豫最大步數(shù)為100,采用共軛梯度算法(CG 算法)優(yōu)化原子的位置,k點選取因晶體結(jié)構(gòu)的差異和計算物理量的不同而略有不同.本文用于擬合的數(shù)據(jù)庫包括:B2 相的晶格常數(shù)如表1 所列;三種純金屬Ti,Al,Nb 的內(nèi)聚能,如表2 所列;B2 相的彈性常數(shù)如表4 所列;B2 相的未弛豫的單空位形成能和雙空位形成能如表5 所列;B2 相的置換原子形成能和換位原子形成能如表6 所列;B2 相的內(nèi)聚能如表8 所列.

表1 B2-Ti2AlNb 的晶格常數(shù)和生成焓(單位: eV)Table 1.Lattice constants and formation enthalpy for B2-Ti2AlNb (in eV).

表2 單質(zhì)Ti,Al,Nb 的內(nèi)聚能(單位: eV)Table 2.Cohesive energies for Ti,Al and Nb (in eV).

表1 中同時列出了實驗值和文獻[11]的計算結(jié)果,可以看出,本文的計算結(jié)果與實驗值和文獻[11]的結(jié)果都非常接近.

2.3 原子間勢函數(shù)參數(shù)化

為了對原子間勢函數(shù)進行參數(shù)化,以及檢驗原子間勢的質(zhì)量,需要使用原子間勢表示各個物理量.本文考慮了長程相互作用,原子間勢的截斷距離設(shè)置在第六和第七近鄰之間,且需要給出Ti,Al,Nb 三種元素之間的對勢及電子云之間的相互作用,因此解析表達式比較復(fù)雜.為此本文主要通過程序設(shè)計來建立物理量與原子間勢的關(guān)系.設(shè)計思想是,先建立原子構(gòu)型模塊,包含原子空間坐標、原子類型和原子編號等信息,然后在此基礎(chǔ)上,根據(jù)相應(yīng)的理論進行各種物理量的計算.這里只介紹本文使用的彈性常數(shù)和空位形成能的計算公式.

根據(jù)彈性理論,可以得到平衡狀態(tài)下彈性常數(shù)與原子間勢之間的關(guān)系為[20]

根據(jù)(1)式,可以得到平衡狀態(tài)下空位形成能與原子間勢的關(guān)系式為[21]

在確立了物理量與原子間勢的關(guān)系和擬合數(shù)據(jù)庫之后,利用能量最小化的方法,同時考慮了物理量的權(quán)重,擬合了原子間勢函數(shù)的參數(shù),如表3 所列.

表3 Ti2AlNb 合金EAM 勢的擬合參數(shù)Table 3.Fitting parameters of EAM potential for Ti2AlNb alloys.

在擬合之前,已將B2 相構(gòu)型中的Ti 原子的背景電子密度歸一化,由此將限制sTi,sAl,sNb三者之間的關(guān)系,為了確保電子密度函數(shù)曲線在x軸上方,在擬合過程中已限制三者的符號均為正.對于Ti 原子,在電子密度歸一點處嵌入能函數(shù)達到最低點,而對于Al 和Nb 原子而言,在擬合過程中盡力使其接近1.其中參數(shù)sTi,sAl,sNb的值分別為0.3333,0.4646,4.1557,參 數(shù)gTi,gAl,gNb的值分別為—7.3047×10—7,—0.0640,—0.6733,由此可以得到EAM 型原子間勢的各函數(shù)如圖1 所示.

3 結(jié)果與討論

作為對原子間勢的測試及應(yīng)用,利用開發(fā)的Ti2AlNb 合金的EAM 勢研究了B2-Ti2AlNb 的彈性常數(shù)、點缺陷、表面能、熔點性能,以及Ti2AlNb合金三種有序相的E-V關(guān)系.

3.1 彈性常數(shù)

B2-Ti2AlNb 晶胞的空間點陣屬于正方晶系,其晶格常數(shù)滿足a=b≠c,α=β=γ=90°,彈性常數(shù)[22]滿足C11=C22,C55=C66,C13=C23,計算比較了B2 相的6 個彈性常數(shù)(C11,C33,C44,C55,C12,C13),利用本文開發(fā)的EAM 勢和第一原理的計算結(jié)果如表4 所列.

由表4 可以看出,本文開發(fā)的EAM 勢的計算結(jié)果與第一原理(DFT)計算結(jié)果基本一致.

表4 B2-Ti2AlNb 的彈性常數(shù)(單位: GPa)Table 4.Elastic constants for B2-Ti2AlNb (in GPa).

3.2 空位形成能

單空位及雙空位形成能[23]的計算方法分別為:

其中,Etot(N)表示理想超胞總能,Etot(N-1)表示去掉一個原子的超胞總能,Etot(N-2)表示去掉兩個原子的超胞總能,和分別表示去掉的單質(zhì)原子A 和原子B 的內(nèi)聚能.本文所有關(guān)于空位形成能的計算結(jié)果均列于表5 中,表5 中第二列表示去掉的原子類型,后面的數(shù)字表示第幾近鄰,例如AlNb-1 表示去掉超胞中最近鄰的Al 和Nb 兩個原子.

表5 B2-Ti2AlNb 單空位形成能和雙空位形成能(單位: eV)Table 5.Mono-vacancy formation energies and divacancy formation energies for B2-Ti2AlNb (in eV).

表5 中標注“This work”這一列的值是利用本文開發(fā)的EAM 勢編程計算得到的,標注“Thislmp”的這一列表示利用本文開發(fā)的EAM 勢,使用Lammps 軟件進行模擬計算的結(jié)果.標注“EAM-fa”這一列表示本文利用Farkas[24]開發(fā)的EAM 勢的計算結(jié)果,計算工具使用的依然是Lammps 軟件.由圖2 所示的空位形成能的對比圖可以看出,本文開發(fā)的EAM 勢,用兩種方法的計算結(jié)果幾乎完全一致,相比于用Farkas[24]的勢整體上更接近第一原理計算結(jié)果.

圖2 B2-Ti2AlNb 單空位形成能與雙空位形成能Fig.2.Mono-vacancy formation energies and di-vacancy formation energies for B2-Ti2AlNb.

3.3 置換和換位原子形成能

利用本文開發(fā)的EAM 勢計算了B2 相的置換原子形成能,同時計算了截斷距離范圍內(nèi)交換原子位置的形成能,即換位原子形成能.在圖3 中已經(jīng)清晰地表示出原子換位及其表征情況,其中紅色、綠色、藍色的原子分別表示Ti,Al,Nb.置換原子形成能為

圖3 B2-Ti2AlNb 原子位置換位示意圖Fig.3.Atom position transposition diagram for B2-Ti2AlNb.

EN-x-y-z+k表示有缺陷的超胞總能,EN表示理想超胞的總能,表示單質(zhì)Ti,Al,Nb,K原子的內(nèi)聚能,K表示用于置換的原子的類型.x,y,z的值可以是0 或±1,分別表示被置換原子Ti,Al,Nb 的數(shù)目,k表示用于置換的原子的數(shù)目.例如滿足x=z=0,y=1,K=Ti 的組合表示在Al 原子的點陣上置換一個Ti 原子,用(T iAl)表征,其他置換原子的形成能也采用這種表征方式.對于換位原子的形成能可表示為

式中表示換位后超胞的總能,表示理想超胞的總能.

從表6 可以看出,本文開發(fā)的EAM 勢計算的多種置換原子和換位原子的形成能更接近第一原理計算結(jié)果,從圖4 和圖5 可以更清晰地看到這一結(jié)果,與第一原理相比,各種形成能的誤差最大也沒有超過0.233 eV,表明該EAM 勢能很好地描述Ti2AlNb 合金的點缺陷性質(zhì).

圖4 B2-Ti2AlNb 置換原子形成能的誤差Fig.4.Deviations of substitutional atom formation energies for B2-Ti2AlNb.

圖5 B2-Ti2AlNb 換位原子形成能的誤差Fig.5.Deviations of transposition atom formation energies for B2-Ti2AlNb.

表6 B2-Ti2AlNb 置換原子和換位原子形成能(單位: eV)Table 6.Substitutional atom and transposition atom formation energies for B2-Ti2AlNb (in eV).

3.4 表面能

平衡態(tài)下,理想晶胞中的原子受周圍原子作用的合力為零,處于能量最低的狀態(tài),但位于表面上的原子所受的周圍原子的合力不為零,因而能量較高,減去沒有表面時體系的能量再除以該表面的面積,即為該表面的表面能,表面能的計算公式為

式中為表面能;ES表示具有表面時體系的總能;Et表示沒有表面時體系的總能;Ssur為體系中形成表面的面積.本文計算了B2 相中(100)面和(110)面的表面能,形成表面的體系在平行于表面的方向上利用了周期性邊界條件,在垂直于表面的方向上設(shè)置了擴胞和真空層,具體的計算結(jié)果如表7 所示,表7 中標記一撇“'”的表面表示對該表面進行了能量最小化.從表7 中可以發(fā)現(xiàn)本文開發(fā)的EAM 勢關(guān)于(100)和(110)兩表面能的計算結(jié)果與第一原理計算結(jié)果趨勢一致,都表明(110)面的表面能更小.

表7 B2-Ti2AlNb 的表面能(單位: eV)Table 7.Surface energies for B2-Ti2AlNb (in eV).

3.5 E-V 曲線

為了檢驗原子間勢偏離平衡態(tài)時的表現(xiàn),利用本文開發(fā)的EAM 勢推導(dǎo)了Ti2AlNb 合金三種有序相的E-V關(guān)系.考慮到合金中Ti,Al,Nb 三種原子的周圍環(huán)境不同,因此計算了分別以Ti,Al,Nb 三種原子為參考原子的E-V關(guān)系,合金中以Ti 原子為參考原子的E-V關(guān)系為

式中,r1表示B2 相合金中最近鄰距離,r2表示次近鄰距離,以此類推,rn表示第n近鄰距離.電子密度函數(shù)fTi(rn) 表示為第n近鄰的Ti 原子貢獻的電子密度,前面的數(shù)字表示該原子的數(shù)目.

同理,可以得到分別以合金中Al 原子和Nb 原子為參考原子的E-V關(guān)系為

對上述三種E-V關(guān)系進行加權(quán)平均,便可得到Ti2AlNb 合金中平均每個原子的E-V關(guān)系,如(20)式所示:

本文計算的B2相E-V關(guān)系曲線如圖6 所示,同時為了比較,根據(jù)Rose 的經(jīng)驗公式[25]計算了E-V關(guān)系,同時畫在圖6 中.

由圖6 可以看出,在距離小于約3.3 ?的范圍內(nèi)擬合曲線與Rose 曲線吻合非常好,同時表明曲線在體系處于平衡狀態(tài)時的曲率相近,彈性模量一致.超過這個范圍兩曲線雖然略有差距,但整體變化平穩(wěn),沒有出現(xiàn)過多的振蕩.

圖6 B2-Ti2AlNb的E-V 曲線與Rose 曲線Fig.6.E-V curve and Rose curve for B2-Ti2AlNb.

為了驗證本文開發(fā)的EAM 勢對于Ti2AlNb合金不同相的描述能力,使用第一原理方法計算了平衡狀態(tài)下B2 相、D019相、O 相在不同體積下的內(nèi)聚能,計算結(jié)果見表8,三種有序相的E-V曲線如圖7 所示.從圖7 可以發(fā)現(xiàn): 對于B2 相,該EAM 勢的計算結(jié)果非常接近第一原理的計算結(jié)果,Farkas[24]的結(jié)果在體積膨脹超過20% 時變化非常緩慢,不太適合這個范圍;對于D019相,本文的EAM 勢的計算結(jié)果相比Farkas[24]的計算結(jié)果更接近第一原理計算結(jié)果;對于O 相,本文的結(jié)果不僅與第一原理的結(jié)果的趨勢完全一致,偏差也非常小,Farkas[24]的結(jié)果與第一原理的計算偏差也不大,而且曲線的左側(cè),即當晶體被壓縮時,Farkas的結(jié)果與第一原理的結(jié)果趨勢也完全一致,但是曲線右側(cè)的變化趨勢與第一原理計算結(jié)果不太一致,也就是當晶體膨脹時,Farkas 的結(jié)果比第一原理的結(jié)果變化緩慢.

圖7 Ti2AlNb 合金的E-V曲線 (a) B2 相;(b) D019 相;(c) O相Fig.7.E-V curve of Ti2AlNb alloys: (a) B2 phase;(b) D019 phase;(c) O phase.

3.6 熔點

利用本文開發(fā)的EAM 勢通過分子動力學(xué)模擬研究了B2 相的熔化轉(zhuǎn)變過程,模擬盒子的大小為3456 個原子的超胞,在初始B2 相構(gòu)型下保持初始溫度,并在NVT 系綜下進行弛豫,升溫的過程在NPT 系綜下進行,期間原子體積在0 至1400 K 附近的溫度范圍內(nèi)呈線性增加,如圖8 所示.圖8中可以看出B2-Ti2AlNb 合金在1400 K附近發(fā)生一級相變,單位體積發(fā)生突變,正如圖9 所示,體系在此溫度附近前后從有序構(gòu)型轉(zhuǎn)變?yōu)闊o序狀態(tài).圖9 中紅色、綠色、藍色的原子分別表示Ti,Al,Nb.

圖8 B2 相的每個原子的體積隨溫度的變化Fig.8.Volume variation of each atom for B2 phase as functions of temperature.

圖9 B2 相隨溫度變化的截面快照 (a) 296 K;(b) 754 K;(c) 1428 K;(d) 1390 K;(e) 1449 K;(f) 1728 KFig.9.Cross-section snapshots for B2 phase as functions of temperature: (a) 296 K;(b) 754 K;(c) 1428 K;(d) 1390 K;(e) 1449 K;(f) 1728 K.

從Ti-Al 平衡相圖[26]中看出,Ti3Al 合金的熔點大于1900 K,向Ti3Al 基合金中摻雜Nb 可以提升合金的綜合性能,而Ti-22 Al-26 Nb 合金B(yǎng)2 相的熔點大約為1873 K[27],以此為參考,利用本文開發(fā)的EAM 勢計算的Ti2AlNb 合金的熔點比此溫度低了約25%.由于本文主要針對Ti2AlNb 合金的力學(xué)性能開發(fā)EAM 勢,擬合數(shù)據(jù)庫中未包含與溫度相關(guān)的數(shù)據(jù).一般來說,同一種勢函數(shù)形式很難同時描述低溫和高溫的性能,考慮有兩種改進方案,一是在擬合數(shù)據(jù)庫中加入相變等與溫度相關(guān)的數(shù)據(jù),在保證不影響對力學(xué)性能描述的基礎(chǔ)上,改善勢函數(shù)描述高溫性能的能力;另一種是針對高溫性能開發(fā)新的原子間勢,只用于預(yù)測高溫性能.

表8 Ti2AlNb 合金的內(nèi)聚能(單位: eV)Table 8.Cohesive energies for Ti2AlNb alloys (in eV).

4 結(jié)論

本文使用傳統(tǒng)方法開發(fā)了Ti2AlNb 合金的EAM 勢.主要采用靜態(tài)計算和分子動力學(xué)模擬兩種手段考察了本文開發(fā)的EAM 勢的質(zhì)量.

1) 靜態(tài)計算表明,該EAM 勢很好地再現(xiàn)了B2相Ti2AlNb 合金的彈性常數(shù)、內(nèi)聚能和多種缺陷形成能,能很好地描述Ti2AlNb 合金B(yǎng)2 相、D019相、O 相的E-V系.其中B2 相的E-V關(guān)系在近距離約3.3 ?以內(nèi)與Rose 曲線符合的非常好.

2) 利用Lammps 軟件模擬計算了B2-Ti2AlNb合金的各種空位形成能,計算結(jié)果均與第一原理計算結(jié)果符合得很好.置換原子形成能和換位原子形成能計算結(jié)果與第一原理結(jié)果相比,其大部分缺陷形成能的誤差均在0.100 eV 以下,對于Al 原子置換Ti 原子的置換原子形成能,誤差在0.165 eV 之內(nèi),Nb 原子置換Al 原子的置換原子形成能,誤差在0.233 eV 以內(nèi).利用Lammps 軟件模擬計算了B2 相Ti2AlNb 的表面能,計算結(jié)果表明B2 相的(110)面比(100)面更穩(wěn)定,表面能更低,與第一原理的計算結(jié)果一致.

3) 利用Lammps 軟件模擬研究了Ti2AlNb 合金B(yǎng)2 相的熔化轉(zhuǎn)變過程,熔點的模擬結(jié)果大約為1400 K,低于實驗結(jié)果.主要由于開發(fā)勢函數(shù)的擬合庫中未包含與溫度、相變相關(guān)的物理量.

本文開發(fā)的Ti2AlNb 合金的EAM 勢在描述B2 相的彈性性能,點缺陷性能、表面能及E-V關(guān)系性能方面的能力較好,預(yù)測可用于與這些性能相關(guān)的力學(xué)性質(zhì)方面的研究.