鎢中自間隙原子團簇擴散行為的分子動力學模擬

李小椿，潘新東,2，周海山,*，羅廣南,2

(1.中國科學院 等離子體物理研究所，安徽 合肥 230031；2.中國科學技術大學，安徽 合肥 230026)

核聚變能因其對環境友好和不產生放射性廢料的高安全性等優勢成為一種清潔的安全能源，其燃料——氘在海水中幾乎取之不盡、用之不竭，是解決能源危機最有效的方法之一。然而，可控熱核聚變技術現在還處于一個相當不成熟的階段。因為庫侖排斥作用使核聚變反應非常困難，須使反應氣體電離成為等離子體并加熱至極端的高溫。當前一般使用磁約束熱核聚變的方法，托卡馬克裝置是迄今為止最有發展前途的磁約束熱核聚變裝置[1-2]。國際上已啟動國際熱核聚變實驗反應堆(ITER)計劃[3-4]，為將來進一步發展聚變反應堆奠定基礎。中國也正在進行聚變試驗堆(CFETR)的工程設計[5-7]。

聚變能源應用的最終實現除要解決可控熱核聚變這一物理問題外，在很大程度上取決于可控熱核聚變裝置托卡馬克以及未來反應堆中關鍵材料問題的解決。核能界公認聚變堆材料是開發核聚變能的最關鍵技術之一，其中面向等離子體材料(PFM)的選擇尤為關鍵[8]。在聚變堆嚴酷的工況環境下，PFM受三重輻照：等離子體輻照、熱輻照、中子輻照。高能聚變中子輻照導致PFM中產生大量缺陷，這些缺陷的遷移和聚集會引起材料的腫脹和變形，從而嚴重影響材料的力學性能，降低部件的使用壽命，最終必將導致材料在服役條件下的失效。如此極端的服役條件對PFM的開發提出了巨大的挑戰[9]。

鎢(W)材料被視為未來聚變堆中最可能全面使用的PFM，其主要優勢有高熔點、高熱導率、不與氫發生化學刻蝕以及相對低的氫滯留等[10-12]。中國科學院等離子體研究所的EAST托卡馬克將逐步過渡到全W的PFM[13]，同時W也是ITER和CFETR的重要候選材料[5,14]。

在未來聚變堆真實環境下，氘氚聚變產生的14 MeV高能中子輻照將在材料中產生嚴重的原子離位損傷和各種缺陷積累[15]。其中自間隙原子(self-interstitial atom, SIA)及其團簇是中子輻照損傷中最常見的缺陷種類，因此研究SIA團簇聚集及其動力學擴散行為具有重要意義。由于純中子源輻照均存在樣品活化帶來的后續問題，因此國際上一直試圖利用離子束輻照產生缺陷，以避免中子活化帶來的巨大麻煩[16-17]。材料模擬是研究極端條件下材料行為的有效工具，其中分子動力學(MD)模擬是研究材料輻照損傷效應的重要工具[18]，是材料研究領域中對微觀結構和微觀過程研究的重要手段。本文采用MD模擬研究不同尺寸SIA團簇的擴散行為，擬為更大尺度的動力學蒙特卡羅和團簇動力學模擬提供準確的輸入參數，為正確掌握和評價W中子輻照行為提供依據。

1 研究方法

在分子靜力學和MD模擬方面采用開源并行程序LAMMPS[19]，LAMMPS是經典的免費開源MD模擬軟件。本文所涉及的模擬均通過LAMMPS軟件實現。MD模擬的準確性極大依賴于模擬選用的勢函數，為正確掌握SIA團簇的性質，針對不同的性質，本文選用兩套勢函數開展相關模擬工作。

Ef(SIAn)=Etot(SIAn)-NWEc(W)

(1)

其中：Etot(SIAn)為包含n個SIA的體系總能量；NW為總的W原子數；Ec(W)為單個W原子在完美體心立方晶格中的內聚能。單個SIA與SIA團簇的結合能定義為：

Eb(SIA+SIAn-1)=

Ef(SIA)+Ef(SIAn-1)-Ef(SIAn)

(2)

v=v0exp(-Em/kBT)

(3)

其中：v0為躍遷頻率指前因子；Em為躍遷激活能；kB為玻爾茲曼常數；T為溫度。對式(3)兩邊分別取對數，可得到：

lnv=lnv0-Em/kBT

(4)

可見躍遷頻率的對數與溫度的倒數滿足線性關系，線性擬合后就可得到躍遷頻率指前因子和躍遷激活能。

2 結果與討論

2.1 SIA團簇的形成能與結合能

首先采用分子靜力學模擬獲得了數量為1～30個的1/2〈111〉和〈100〉 SIA團簇的穩定結構、形成能和結合能，這部分模擬采用Chen勢完成。圖1示出了1/2〈111〉和〈100〉 SIA團簇的部分穩定構型。表1列出了這些團簇的形成能和結合能，作為參考也列出了第一原理計算的結果[24]。由表1可看出，1/2〈111〉 SIA團簇的形成能要低于〈100〉 SIA團簇，說明1/2〈111〉 SIA團簇較〈100〉 SIA團簇穩定。SIA團簇的形成能隨團簇尺寸的增長而增加，且SIA與SIA團簇的結合能非常高，說明SIA團簇聚集后會穩定存在。

圖1 1/2〈111〉和〈100〉 SIA團簇的最穩定結構Fig.1 The most stable structure of 1/2〈111〉 and 〈100〉 SIA clusters

表1 1/2〈111〉和〈100〉 SIA團簇的形成能與結合能Table 1 Formation energy and binding energy of 1/2〈111〉 and 〈100〉 SIA clusters

為給后續大尺度模擬提供更準確的輸入，研究了1/2〈111〉和〈100〉 SIA團簇尺寸大于30時的形成能，結果如圖2所示。受計算量限制，無法給出每個尺寸的SIA團簇的形成能，可采用以下經驗公式擬合[24]：

(5)

采用團簇尺寸大于30之后的形成能擬合，得到了1/2〈111〉 SIA團簇的擬合參數為：a0=6.801 17、a1=-11.021 26、a2=60.138，而〈100〉 SIA團簇的擬合參數為：a0=8.845 46、a1=-20.636 8、a2=93.665 1。這些參數可為后續的大尺度模擬提供每個尺寸的SIA團簇的形成能。相應地，各團簇之間的結合能也可根據形成能輕易獲得。據此，可為后續的大尺度模擬提供W中最常見的兩種SIA團簇的形成能和結合能的完備數據庫。

圖2 大尺寸SIA團簇的形成能與經驗擬合Fig.2 Formation energy of large SIA clusters and empirical fitting

2.2 SIA團簇的擴散行為

采用Ackland勢模擬研究不同尺寸的1/2〈111〉 SIA團簇在不同溫度下的擴散行為。SIA團簇的尺寸變化范圍為1～109個，溫度范圍為300～900 K。對于每個SIA團簇在每個溫度下的模擬，均通過隨機產生初始速度的方法運行20次，因此有共20 ns的模擬數據，隨后采用OVITO軟件分析SIA團簇的運動軌跡。

通過對SIA運動軌跡的分析發現，在300～600 K溫度條件下，單個SIA在MD模擬時間尺度內未觀察到轉向，表現為一維方向的快速運動，如圖3a所示，SIA只在x方向(〈111〉方向)移動，而在另外兩個方向幾乎不動。只有當溫度大于700 K時，單個SIA才可轉向，如圖3b所示，SIA在0.37 ns發生轉向后，易在0.7 ns時再次轉回穩定的〈111〉方向。隨著溫度的升高，單個SIA越易轉向。在700 K的模擬下，20個模擬里只有2個模擬觀察到了SIA的轉向。而在800 K的模擬下，有13個模擬觀察到SIA的轉向。

圖3 單個SIA在600 K(a)和800 K(b)下的擴散坐標與時間的關系Fig.3 Time dependence of diffusion coordinates of single SIA at 600 K (a) and 800 K (b)

而SIA團簇尺寸超過2個后，在300～900 K時，這些SIA團簇均非常穩定，表現為一維方向的運動。需說明的是，雖然模擬的20 ns的時間內未觀察到2個以上SIA團簇的轉向，并不代表這些尺寸的SIA團簇不會發生轉向，受限于MD的局限，可能需更長的時間才能觀察到。將模擬溫度升高到1 200 K可發現，2個SIA團簇在1 000～1 200 K、3個SIA團簇在1 200 K時可觀察到SIA團簇的轉向，而4個以上的SIA團簇在1 200 K下也未發現轉向。

2.3 SIA團簇的躍遷頻率

為掌握SIA團簇的擴散行為，進一步分析了SIA團簇在不同溫度下的躍遷頻率，結果如圖4所示。SIA團簇尺寸越大，躍遷頻率越低，擴散越緩慢，而溫度越高，SIA團簇擴散越快。根據式(4)對SIA團簇在不同溫度下的躍遷頻率的對數與溫度的倒數進行線性擬合，可得到躍遷頻率指前因子和躍遷激活能，擬合曲線如圖5所示，在300～900 K下，SIA團簇的躍遷頻率vn的對數與溫度的倒數有著良好的線性關系，說明可采用式(3)定義躍遷頻率與溫度的關系。根據擬合結果，在表2列出了各種尺寸SIA團簇的躍遷能壘Em和躍遷頻率指前因子v0。躍遷頻率指前因子隨SIA團簇尺寸的增加而遞減，而躍遷能壘卻幾乎與SIA團簇尺寸無關。通過分析SIA團簇的運動軌跡可發現，SIA團簇里的每個SIA均在一條獨立的〈111〉方向里移動，因此整個團簇的躍遷能壘與單個SIA的躍遷能壘類似。

圖4 SIA團簇在不同溫度下的躍遷頻率Fig.4 Vibration frequencies of SIA clusters at different temperatures

圖5 SIA團簇的躍遷頻率與溫度倒數的關系Fig.5 Relationship between vibration frequency and reciprocal temperature of SIA clusters

表2 SIA團簇的躍遷頻率指前因子和躍遷能壘Table 2 Pre-exponential factor of vibration frequency and diffusion barrier of SIA clusters

在真實的環境中，SIA團簇有各種不同的尺寸，因此需根據經驗公式，外推出各種尺寸SIA團簇的躍遷頻率指前因子和躍遷能壘。SIA團簇的躍遷頻率與團簇尺寸滿足以下關系：

vn=v0n-Sexp(-〈Em〉/kBT)

(6)

其中：S為可擬合的參數；〈Em〉為平均躍遷激活能。對表2中間隙團簇尺寸大于7的數據進行擬合，可得到1/2〈111〉 SIA團簇的相關參數：v0=2.226×1012Hz、S=0.498和〈Em〉=0.021 7 eV。

圖6為SIA團簇的躍遷頻率指前因子與團簇尺寸的關系。從圖6可看出，擬合參數可很好描述1～109個SIA團簇的躍遷頻率指前因子的關系。據此，就可為后續的大尺度模擬提供每個SIA團簇的躍遷頻率和躍遷能壘等參數。

圖6 SIA團簇的躍遷頻率指前因子與團簇尺寸的關系Fig.6 Relationship between pre-exponential factor of vibration frequency and SIA cluster size

3 結論

采用MD模擬研究了W中1/2〈111〉和〈100〉 SIA團簇的穩定結構，發現W中SIA團簇最穩定結構是1/2〈111〉 SIA團簇結構，且SIA與SIA團簇的結合能非常高，說明SIA團簇聚集后會穩定存在。給出了1/2〈111〉和〈100〉 SIA團簇的形成能和結合能，并給出了計算大尺寸SIA團簇形成能的經驗參數。隨后研究了1/2〈111〉 SIA團簇的動力學擴散行為，單個SIA在溫度高于700 K時易擴散和轉向。而兩個以上的SIA團簇在300～900 K時主要表現為一維方向的運動，未觀察到轉向。根據不同溫度下的躍遷頻率，得到各種尺寸SIA團簇的躍遷能壘和躍遷頻率指前因子。為準確描述各種尺寸SIA團簇的動力學行為，給出了一套計算SIA團簇躍遷頻率的經驗參數。目前的模擬結果將有助于理解W中SIA的團簇行為，從微觀角度理解和解釋實驗，并為大尺度的模擬提供準確和完備的輸入參數。