冉琴 王歡 鐘睿 伍建春 鄒宇 汪俊

(四川大學原子核科學技術研究所,輻射物理及技術教育部重點實驗室,成都 610064)

1 引 言

鎢因其具有高熔點、低濺射產(chǎn)額、高熱導率等特性而被視為聚變堆中非常具有應用前景的面向等離子體候選材料[1].這些材料處于嚴苛的聚變環(huán)境,將受到中子和氦氫同位素的轟擊,產(chǎn)生空位、自間隙原子 (self-interstitial atoms,SIAs)以及位錯等缺陷.這些缺陷的擴散、聚集及演化將導致材料宏觀性能的改變(如: 腫脹和變形),最終會影響到裝置的安全性問題[2?6].不同于空位,SIAs 及其團簇具有快速的移動性,其擴散行為被認為是決定材料微觀性能的關鍵因素[7?10].其擴散系數(shù)及擴散機制等方面的精確知識是我們能正確理解和預測材料微觀演化行為的必要條件.因此,鎢中SIAs的動力學行為研究至關重要.

目前,鎢中SIAs及其團簇的相關特性已經(jīng)得到了廣泛研究[11?17].分子動力學模擬和第一性原理計算表明鎢中單自間隙原子(single selfinterstitial atom,1-SIA)最穩(wěn)定的構型是結構[11,12,18].Tsong和Casanova[13]通過場粒子顯微鏡觀察到鎢中1-SIA與雙自間隙原子(di-interstitial atoms,2-SIAs)團簇沿著方向的擴散.Amino等[14]利用透射電子顯微鏡觀察了納米級SIAs團簇的形成過程,并用動力學蒙特卡羅研究了團簇的擴散行為.他們得出低溫下鎢中SIAs擴散是一維運動,這些實驗結果都與第一性原理[1]和分子動力學[4,15]研究結果一致.值得一提的是,同為體心立方 (body center cubic,BCC)晶體結構的Fe中SIAs的擴散卻是三維運動,這表明在BCC金屬中存在不同的SIAs擴散機制.Zhou等[5,15]利用分子動力學方法研究了鎢中SIAs團簇的形成和運動特性,結果表明2-SIAs形成的團簇在較低溫度下非常穩(wěn)定,即使在900 K也未發(fā)生轉向,保持一維運動,SIAs團簇的遷移能與團簇尺寸大小無關.最近,Swinburne 等[19]用密度泛函理論結合量子分子動力學方法研究了金屬鎢中SIAs在低溫下的擴散行為,研究表明SIAs在10 K左右的溫度下也具有快速的遷移行為.盡管有大量文獻研究金屬鎢中SIAs,但對鎢中具有不同構型SIAs團簇的擴散行為鮮見報道,對其擴散機理仍有待研究.分子動力學不僅能在原子尺度模擬缺陷的微觀演化行為,揭示擴散機理,而且在多尺度模擬中起著重要的作用.它能為動力學蒙特卡羅模擬提供輸入?yún)?shù)(如: 擴散系數(shù)和擴散機制),從而進行長時間、大空間尺度的模擬以獲得實驗上可觀察的結果,進而解釋和預測實驗.2-SIAs作為一種重要的間隙團簇生長核,是大尺寸SIAs團簇的生長基礎,研究其不同構型無論是對理論模擬還是實驗研究都有一定的指導意義.

本文基于新發(fā)展的鎢原子間相互作用勢,采用分子動力學方法研究了鎢中SIAs的擴散行為.主要考察了具有不同構型的2-SIAs隨溫度變化的擴散行為.通過將微動彈性帶 (nudged elastic band,NEB)算法獲得的SIAs遷移能與阿倫尼烏斯(Arrhenius)擬合的結果進行對比,闡明了鎢中1-SIA與2-SIAs的擴散系數(shù)隨溫度的變化規(guī)律.

2 模型與方法

分子動力學模擬的核心問題是原子間作用勢的選取,它直接關系到能否準確描述原子間的相互作用以及真實可靠地反映所模擬的物理過程.Bonny等[20]根據(jù)實驗和密度泛函理論計算結果,詳細評估了19種鎢勢的基本物理屬性,特別是點缺陷和位錯屬性.結果表明,Marinica 等[21]發(fā)展的EAM2勢總體性能最佳.因此,本文也采用EAM2勢,它是根據(jù)鎢晶體的晶格常數(shù)、彈性常數(shù)、內(nèi)聚能、點缺陷形成能等物理參數(shù)和來自于第一性原理的原子間作用力數(shù)據(jù)擬合而成,形式如下:

式中 ? 為兩體作用函數(shù); ρ 為有效電子密度;F 為嵌入函數(shù); θ(x) 是 Heaviside 階躍函數(shù),當時,其他情況下上述參數(shù)的具體取值見文獻[21].

在進行分子動力學模擬之前,首先構建一個BCC 結構的鎢晶體,基底共 20 層,每層 10 × 10個原子排列在(001)面.根據(jù)研究的具體需要,引入SIAs.本文所涉及的SIAs結構包括1-SIA (圖1)和多種構型的 2-SIAs (圖 2),2-SIAs結構包括最近鄰(2-SIAs-1st)、次近鄰(2-SIAs-2nd)、三近鄰結構(2-SIAs-3rd)以及非平行結構.為了將體系擴展到無限大情形,分別沿著 [100],[010],[001]方向加上周期性邊界條件.值得一提的是,當SIAs運動出基底盒子時,由于周期性邊界條件的作用,SIAs將會從另外一面進入基底,這將會丟失掉真實的擴散信息.為了獲得SIAs在無限大體系中的運動軌跡,在SIAs擴散出基底時修正了其位置信息,具體的修正方式見文獻 [22,23].本文所有的模擬都在NVT (固定粒子數(shù)、體積、溫度)系綜下進行,為了將系統(tǒng)升至所需溫度,將系統(tǒng)中的每個原子按Maxwell隨機分配速度,然后弛豫.在演化過程中利用 Wigner-Seitz cell方法判斷 SIAs的位置,并按一定的時間間隔保存位置信息,從而得到SIAs 運動的軌跡.演化的時間步長設為 10–15s (1 fs),每一種結構配置演化時間為10 ns.

圖1 鎢中 1-SIA 的結構圖 (紫色球為 SIAs結構,藍色球為格點原子)Fig.1.1-SIA configuration in W.The purple sphere represents the SIA;the blue one represents the lattice atom.

圖2 2-SIAs的不同構型圖 (a),(b),(c),(d)分別代表最近鄰、次近鄰、三近鄰以及非平行結構的結構示意圖;右上方的插圖分別代表這幾種結構方向的視圖;紫色球為SIAs,藍色球為格點原子Fig.2.Schematic illustrations of the 2-SIAs with different configurations: (a),(b),(c),(d) Represent the configuration of the 2-SIAs-1st,2-SIAs-2nd,2-SIAs-3rd and the non-parallel SIAs,respectively.Insets represent the views corresponding to their orientations;the purple sphere stands for the SIA and the blue one stands for the lattice atom.

3 結果與討論

3.1 形成能與束縛能

為了考察各種SIAs結構的相對穩(wěn)定性,首先從能量角度考察了其形成能和束縛能.重點研究了2-SIAs的平行結構和非平行結構,平行結構包括圖2(a)—(c)所示的最近鄰(2-SIAs-1st)、次近鄰 (2-SIAs-2nd)、三近鄰結構 (2-SIAs-3rd);非平行結構如圖2(d)所示.SIAs的形成能為

2-SIAs的束縛能公式:

其中 EB為束縛能,E1為 SIAs 之間無相互作用的系統(tǒng)能量,E2代表兩者結合形成穩(wěn)定結構時系統(tǒng)的最低能量.形成能和束縛能的計算結果,分別見表1和表2.通過計算得到單個SIAs的形成能為 10.406 eV,這與文獻 [20]結果吻合.可以看出,對于2-SIAs,最近鄰結構具有最低的形成能(18.574 eV),同時具有最高的束縛能 (2.245 eV),是最穩(wěn)定的結構.有趣的是,非平行結構的穩(wěn)定性僅次于最近鄰結構,而高于次近鄰和三近鄰結構,最不穩(wěn)定的結構是三近鄰結構.這些結論與本文隨后的動力學模擬結果一致.

表1 不同缺陷結構的形成能Table 1.Formation energies of self-interstitials with different configurations.

表2 不同缺陷結構的束縛能Table 2.Binding energies of self-interstitials with different configurations.

3.2 單個自間隙原子的擴散

在BCC金屬中,基于不同的材料有兩種最穩(wěn)定的自間隙結構,一種是結構,另一種是結構.金屬鐵中的最穩(wěn)定結構是前一種[5],而鎢中是后一種[11,12].兩種不同的穩(wěn)定結構導致它們呈現(xiàn)出不同的擴散行為.為了了解鎢中1-SIAs的擴散行為,考察了鎢中結構的1-SIA在100—1500 K溫度范圍內(nèi)的擴散,并采用Wigner-Seitz cell方法得到其擴散徑跡(圖3).結果表明:1-SIA 在 600 K 以下,沿著方向表現(xiàn)為快速的一維運動.這與分子動力學[4,5]和第一性原理計算[1]低溫下的擴散結果一致.當溫度達到700 K時,SIA的運動方向會隨機轉向其他三個等價的方向,從而呈現(xiàn)三維運動形式.如圖3所示,隨溫度升高,SIA 轉向越頻繁.需要注意的是,對于不同的勢函數(shù),鎢中SIAs的旋轉溫度不同[22].

圖3 1-SIAs在不同溫度下演化 10 ns的擴散徑跡圖 (a) T=100 K;(b) T=700 K;(c) T=1000 KFig.3.Diffusive trajectories of 1-SIA for temperatures of (a) 100 K,(b) 700 K and (c) 1000 K.

3.3 雙自間隙原子擴散

圖4給出了鎢中2-SIAs的最近鄰結構在不同溫度下的擴散徑跡.在達到轉向溫度1400 K之前,最近鄰結構的擴散一直保持一維運動形式.當溫度達到1400 K時,擴散由方向轉向其他三個等價的方向,從而運動形式由一維轉變?yōu)槿S.隨著溫度的進一步升高,轉向越頻繁,但整個過程中一直保持最近鄰狀態(tài),即使在高達2000 K的溫度下也未解離.表1和表2的數(shù)據(jù)也表明了該結構極為穩(wěn)定,其束縛能為 2.245 eV.

圖4 最近鄰結構在不同溫度下演化 10 ns的擴散徑跡圖 (a) T=100 K;(b) T=1400 K;(c) T=2000 KFig.4.Diffusive trajectories of 2-SIAs-1st for temperatures of (a) 100 K,(b) 1400 K and (c) 2000 K.

圖5 T=500 K 時 sessile 結構在不同時間的結構圖 (a) t=0.5 ns;(b) t=2 ns;(c) t=5 nsFig.5.Views of the sessile cluster obtained by molecular dynamics simulation at different time when T=500 K: (a) t=0.5 ns;(b) t=2 ns;(c) t=5 ns.

關于2-SIAs非平行的情況,主要研究了圖2(d)所示的結構在不同溫度下的動力學行為.在實際模擬時,將2-SIAs沿各自方向移動一定距離以便觀察演化過程及結構變化.模擬表明,非平行的2-SIAs,沿各自方向擴散,一旦靠近,就會形成固著結構,具有相當?shù)姆€(wěn)定性,演化徑跡表明了這一過程.通過對其結構的分析,發(fā)現(xiàn)2-SIAs形成了如圖5所示的sessile結構.sessile結構在一定溫度范圍內(nèi)能長時間(ns)穩(wěn)定存在.當溫度達到1000 K時,sessile結構可以轉化形成更為穩(wěn)定的最近鄰結構,隨后按最近鄰擴散方式進行擴散運動.前面束縛能的計算結果也表明最近鄰結構較sessile結構更穩(wěn)定.

SIAs團簇一直被認為具有高遷移率,在空間尺度上可以遷移相當大距離.直到Bacon等[24?27]的發(fā)現(xiàn)改變了人們的認識,即在某些情況下鋯、鈦、鐵金屬中SIAs團簇不再具有這種移動性,而是形成固著性的團簇.在材料受到輻照之后,在級聯(lián)碰撞初期或熱峰過程中SIAs相互作用都可能會形成sessile結構[28].但鎢中sessile結構鮮見文獻報道.本文研究結果表明,鎢中也存在sessile結構,這有助于我們更全面地理解鎢中的SIAs行為.如前所述,SIAs的動力學行為是決定材料微觀演化過程的重要因素,這方面的精確知識是我們能正確理解和預測材料微觀演化行為的必要條件.sessile結構在一定溫度范圍內(nèi)不移動,與其他構型具有完全不同的行為屬性.在動力學蒙特卡羅模擬中,如果單純地考慮一種構型的SIAs,尤其是忽略掉sessile結構,很可能會導致模擬結果與實際情況出現(xiàn)大的偏差.該結構對材料微觀演化行為的影響程度也是今后理論模擬和實驗研究值得進一步探索的課題.

3.4 擴散系數(shù)

材料中缺陷的擴散行為會直接影響到缺陷的分布、生長演化等過程,從而對材料性能產(chǎn)生影響.擴散系數(shù)是衡量物質擴散能力的重要參數(shù).為了進一步了解鎢中SIAs不同構型的擴散特性,我們計算了它們在不同溫度下的擴散系數(shù).根據(jù)愛因斯坦公式[29],粒子的擴散系數(shù)可以表示為

圖6 不同結構的擴散系數(shù) (圖中實線是根據(jù) Arrhenius關系擬合的結果) (a) 1-SIA;(b) 2-SIAs-1st;(c) 2-SIAs-2ndFig.6.Arrhenius plots of diffusion coefficients of single SIA and di-interstitial atoms in tungsten,which is determined using MD simulations and plotted as a function of the absolute temperature T: (a) 1-SIA;(b) 2-SIAs-1st;(c) 2-SIAs-2nd.

表3 Arrhenius擬合所得各個間隙結構的擴散遷移能與前因子Table 3.Migration energy Em (in eV) and prefactor D0 (in cm2/s) for W clusters diffusion obtained by Arrhenius fitting.

圖7 通過 NEB 方法所得不同結構的遷移能壘 (a) 1-SIA;(b) 2-SIAs-1st;(c) 2-SIAs-2ndFig.7.Migration barriers for SIAs with different structures studied by NEB method: (a) 1-SIA;(b) 2-SIAs-1st;(c) 2-SIAs-2nd.

式中D0是擴散前因子,Em是遷移能,kB是玻爾茲曼常數(shù),T是溫度.隨后將各種間隙結構的擴散系數(shù)進行Arrhenius擬合,獲得了各個間隙結構的擴散遷移能.如表3所列,1-SIA的遷移能Em為0.0274 eV,D0為2.52×10–7m2·s–1;最近鄰結構Em為 0.0386 eV,D0為 2.08×10–7m2·s–1;次近鄰結構Em為 0.0248 eV,Do為1.02×10–7m2·s–1.

然而,最新的研究表明,對于單個鎢 SIAs,其擴散規(guī)律并不滿足Arrhenius關系,而是線性關系[30].這主要是因為,在所考察的溫度范圍下,其能量接近或高于鎢SIAs的擴散遷移能.NEB方法通過搜索初末態(tài)之間的最小能量路徑來確定遷移能壘,是最可靠的方法,能夠準確地獲得缺陷的擴散遷移能[31].為了確定上述幾種結構的擴散規(guī)律能否用Arrhenius關系描述,采用NEB方法計算了這些間隙結構的遷移能.然后對比這兩種方法獲得的結果,若兩者接近,則說明在目前的勢函數(shù)下,擴散滿足Arrhenius關系,否則就不滿足.

圖7給出了由NEB所獲得的反應坐標與能量的關系.1-SIA,2-SIAs-1st和 2-SIAs-2nd 的遷移能 分 別 為 0.018,0.0476 和 0.0174 eV,這 與 由Arrhenius關系獲得的結果(分別為0.0274,0.0386和 0.0248 eV)有明顯差異,對比結果見表 3.這說明1-SIA和2-SIAs的擴散規(guī)律都不滿足Arrhenius關系.于是對這些擴散數(shù)據(jù)進行線性擬合,擬合公式如下:

圖8 不同缺陷的擴散系數(shù)Fig.8.Diffusion coefficient for self-interstitials of different configuration in tungsten determined by molecular dynamics simulations and plotted as a function of the absolute temperature T (the solid lines are linear fits).

式中,kB是玻爾茲曼常數(shù),kB=1.38 × 10–23 J·K–1;T是溫度,單位是K;B和A為擬合常數(shù).如圖8所示,這些間隙結構的擴散系數(shù)與溫度形成良好的線性關系.從圖8可以看出2-SIAs結構的擴散能力明顯低于1-SIA,在一定溫度范圍內(nèi)2-SIAs-1st和2-SIAs-2nd的擴散能力接近.

4 結 論

由于SIAs及其團簇對材料微觀性能的重要影響,本文運用分子動力學方法重點研究了鎢中2-SIAs的不同構型在不同溫度下的擴散行為,考察了它們的熱穩(wěn)定性和運動特性.結果表明: 2-SIAs不同構型的穩(wěn)定性差異較大,擴散行為也有明顯不同.彼此互為最近鄰的雙自隙原子結構最穩(wěn)定,隨著溫度的升高,從一維擴散演變成三維擴散,在2000 K的溫度下都沒有出現(xiàn)解離;次近鄰2-SIAs呈現(xiàn)出一維的擴散行為,在溫度高于600 K將解離成兩個獨立運動的SIAs;而三近鄰結構熱穩(wěn)定性最差,在溫度高于 300 K 就將解離.鎢中2-SIAs按非平行方式擴散會形成具有固著性的sessile結構,這一結構具有相當?shù)姆€(wěn)定性,在一定溫度范圍內(nèi)幾乎不移動,隨著溫度升高將轉變成最近鄰結構,然后按最近鄰結構的方式擴散運動.因此,在討論2-SIAs團簇的動力學特性時不能同一而論,需考慮構型的影響.

最后,通過將NEB方法獲得的SIAs遷移能與Arrhenius擬合的結果進行對比,表明了鎢中1-SIAs和2-SIAs的擴散系數(shù)隨溫度變化的規(guī)律不滿足Arrhenius關系,而是滿足線性關系.這些結果給我們提供了新的視點,加深了我們對鎢中SIAs的擴散行為的認識.