金屬鎢級聯碰撞中勢函數的影響

劉麗霞，陳陽春，邱榮陽，胡望宇,*，鄧輝球

(1.湖南大學 材料科學與工程學院，湖南 長沙 410082；2.湖南大學 物理與微電子科學學院，湖南 長沙 410082)

金屬鎢(W)材料因其具有高熔點、高熱導性和良好的抗濺射性等優良特性，已被認為是用于高溫和輻照應用最具前景的第一壁候選材料[1-4]。在聚變反應堆中，14.1 MeV能量的中子輻照是面向等離子體材料所要承受的主要考驗。當高能中子轟擊到第一壁材料表面，與材料發生劇烈碰撞，會在一定程度上改變材料內部的微觀結構，產生一系列輻照缺陷，如點缺陷、缺陷團簇、位錯環、孔洞等，從而影響材料的宏觀性能，如輻照硬化、脆化、腫脹、蠕變等[5-7]。

在研究輻照模擬時，由于高能強流中子難以直接實現輻照模擬，通常以具有一定動能的初級碰撞原子(PKA)，即被中子撞擊后獲得一部分能量的晶格原子來撞擊材料的點陣原子模擬級聯碰撞過程，通過分析碰撞級聯之后材料發生的變化來分析其輻照損傷情況及其抗輻照性能。由于實驗裝置的限制及科技的發展，計算模擬在研究材料的抗輻照性能方面發揮著越來越重要的作用。目前被廣泛應用的計算模擬方法通常有第一性原理(FP)、分子動力學(MD)、動力學蒙特卡羅(KMC)、團簇動力學(CD)和有限元方法(FEM)等，其中分子動力學方法因其可模擬整個碰撞級聯過程中每個原子的動態演化過程，且所模擬的時間和空間尺度與中子輻照初期級聯碰撞的尺度相一致，而成為模擬級聯碰撞的首選方法。近年來，在缺乏高能中子裝置進行輻照實驗的情況下，盡管國內外在聚變反應堆材料模擬的中子輻照方面進行了許多工作，但仍缺乏系統且完善的模擬平臺來實現聚變材料基于中子輻照下從微觀到宏觀尺度的變化。因此，實施分子動力學研究來建立完善的級聯碰撞數據庫，為后續的模擬提供輸入參考，從而實現不同模擬尺度間信息的有效傳遞是非常有必要的，將會促進核聚變反應堆中鎢材料的發展，具有重要的潛在應用價值。

Nordlund[8]回顧了近年來材料輻照效應計算機模擬的歷史，詳細介紹了輻照模擬的研究方法，充分闡明了分子動力學模擬的重要性。Fikar等[9-10]使用分子動力學方法研究了不同鎢勢函數的級聯模擬，發現盡管不同勢函數的離位閾能不同，各勢函數級聯碰撞之后會產生相似的結果，即穩定狀態下的總缺陷數目相差不大，但缺陷的空間分布有一定的差異，這可能與間隙遷移能有關。Setyawan等[11]使用LAMMPS[12]代碼通過分子動力學方法研究了鎢中高能入射PKA在不同溫度梯度下(300～2 050 K)的級聯碰撞模擬，結果表明穩定狀態幸存的缺陷數目對PKA能量有較強的依賴性，而對溫度的依賴性很??；Yang等[13]研究了金屬鎢的級聯碰撞模擬，結果表明在穩定狀態幸存的缺陷數目雖然隨輻照溫度的升高而略有減少，但對溫度的依賴效果并不顯著；Warrier等[14]研究了面心立方銅和體心立方鎢兩種材料中的級聯碰撞模擬，發現在1～5 keV的PKA低能區間內，當進行使離位原子數目穩定在平均值附近的模擬次數統計時，在鎢體系中比在銅體系中所需統計的次數要少，即鎢體系中更快達到平衡。Fellman等[15]利用分子動力學方法模擬了鎢中含空位團簇和孔洞的級聯碰撞模擬，結果表明當級聯發生區域與空位團簇所在區域重疊時，最后產生的穩定狀態下的缺陷數目會減少，這與Byggm?star等[16]研究的鎢中含間隙缺陷團簇的級聯碰撞模擬結果相似，Byggm?star等發現當級聯發生區域與間隙缺陷團簇所在區域重疊時，級聯產生的最終穩定狀態下的缺陷數目有所下降，且當缺陷團簇的尺寸與級聯區域體積的大小相當時，級聯產生的新增缺陷對數量幾乎為0。Fu等[17]使用一個新開發的WRe合金勢函數研究了鎢及鎢錸合金中高能PKA的級聯模擬，發現缺陷團簇的尺寸和數目均隨PKA能量的增加而增加，且純鎢中間隙團簇和位錯環的移動性比鎢錸合金中要高。Zhang等[18]使用分子動力學方法進行了鎢中在晶界附近不同溫度下的級聯碰撞模擬，結果表明空位缺陷的數目與溫度無關。

目前已有大量使用分子動力學方法模擬鎢基材料級聯碰撞的研究，且所使用的原子間勢函數大多各有優缺點，為了使其更好地用于輻照損傷模擬，鎢勢函數還在不斷地進一步優化改進。隨著勢函數的不斷發展，近幾年來出現了一些新的優化版本的鎢基勢函數。不同原子間勢函數在輻照級聯模擬中的影響有待進一步研究，尤其是不同勢函數在級聯碰撞過程中輻照點缺陷、缺陷團簇和位錯環等分布信息，可為鎢基材料初級輻照損傷的理解及退火過程缺陷長時間的演化模擬提供基礎，為用于輻照模擬鎢勢函數的選擇和優化提供參考。

本文對金屬鎢的勢函數進行比較測試，分析輻照過程中缺陷的產生、缺陷團簇和位錯環分布信息。

1 不同勢函數的中子輻照級聯模擬

折衷考慮到勢函數預測準確性和在大尺度計算時的效率，主要選取的勢函數類型有Embedded-atom方法(EAM)[19]和Finnis-Sinclair(F-S)形式[20]兩種類型。首先選取了兩個目前為止使用得較頻繁的勢函數：一個是由Ackland等[21]開發的F-S形式的鎢勢函數(勢函數AT)，另一個是Juslin等[22]在勢函數AT的基礎上進行了優化和修正的鎢勢函數(勢函數JW)。然后選取了最近發表的幾個各有優勢的鎢基合金勢函數：Marinica等[23]發表了3個EAM勢函數(勢函數MV2、MV3和MV4)，其中的勢函數MV2被Bonny等[24]進行了更新和優化，開發了一個新的WRe合金勢函數，即勢函數MV2-B；勢函數MV4由Setyawan等[25]進行了優化，開發了另一個新的WRe合金勢函數，即勢函數MV4-S；Chen等[26]新發表了一個可準確預測鎢中位錯環形成能的F-S形式WRe合金勢函數(勢函數Chen)。本工作選取上述5個勢函數進行動態級聯碰撞模擬測試。

為更好地模擬在級聯過程中原子間的短程相互作用，首先對各勢函數進行Ziegler-Biersack-Littmark(ZBL)修正，并計算了其離位閾能，具體的連接參數和形式參考文獻[27]。級聯模擬采用三維周期性邊界條件，模擬晶胞的最外層6層原子(3×a0，a0為0.316 52 nm)設置為邊界區域，如圖1所示Region-Ⅰ，模擬均在300 K溫度下進行。在輻照級聯模擬啟動前，首先使用共軛梯度法將系統進行靜態馳豫，然后整個模擬晶胞均通過Nose-Hoover熱浴和控壓方法在等溫等壓(NPT)系綜下進行動態弛豫20 ps，時間步長設置為1 fs。弛豫充分之后，為避免溝道效應，且使得級聯碰撞在盒子中心區域產生，參考文獻[17，28-29]，在模擬晶胞的中心附近隨機選取1個鎢原子作為PKA原子，沿典型高指數〈135〉晶格方向啟動級聯碰撞模擬，一旦超出模擬盒子邊界，則該模擬結果數據無效。在整個級聯碰撞過程中，模擬晶胞外層邊界區域(Region-Ⅰ)通過速度標定法來控制該區域晶格體系溫度保持為300 K，內部區域(Region-Ⅱ)則在微正則(NVE)系綜下自發進行。級聯過程使用變時間步長，且控制在10-3～10-7ps之間，模擬的總時間根據PKA能量(EPKA)大小的不同選擇為20～60 ps。具體的模擬時間、模擬盒子的邊長和模擬事件總數列于表1。由于級聯碰撞過程具有較大的隨機性，為降低誤差，本文采用隨機選擇不同PKA原子的方式對所有勢函數和所模擬PKA能量下的每種情況都進行了15次模擬級聯碰撞過程。本文所使用的分子動力學模擬軟件為LAMMPS[12]，輻照級聯后產生的點缺陷分析方法為Winger-Seitz(W-S)原胞方法[30]，而位錯分析提取為DXA位錯算法[31]。此外，這兩種方法均可通過可視化和數據分析軟件OVITO[32]來實現。

圖1 模擬盒子示意圖Fig.1 Schematic of simulation box

表1 勢函數級聯碰撞模擬參數Table 1 Collision cascade simulation parameter for potential

2 缺陷的產生

表2列出不同勢函數在PKA能量為10 keV和50 keV級聯過后穩定狀態下弗蘭克爾缺陷對(FPs)的數目。由表2可見，盡管各勢函數的離位閾能有差別[27]，但在PKA能量為10 keV下，各勢函數FPs數目的差別很小，這與之前Fikar等[9-10]獲得的結果一致；在PKA能量為50 keV下，不同勢函數之間的FPs數目略有一些差別，如果增加模擬次數其差異可能會縮小。此外，FPs的數目隨著PKA能量的增加而顯著增加。

圖2示出各勢函數FPs數目在PKA能量為50 keV下隨時間的演化，表3列出各勢函數達到熱峰狀態的時間和演化趨于平衡的時間。結果表明，50 keV的PKA能量下，各勢函數到達級聯熱峰的時間(T熱峰)和趨于平衡狀態的時間(T平衡)均相差不大，均在約1 ps產生最多的缺陷，達到熱峰狀態，隨著級聯的演化，間隙與空位缺陷快速湮滅復合，且在10 ps內趨于平衡，缺陷狀態基本穩定，各勢函數間最后穩定狀態下的FPs數目差別不大。

表2 不同勢函數在穩定狀態下的FPs數目Table 2 Number of FPs in stable state obtained by different potentials

圖2 不同勢函數下FPs數目隨時間的演化Fig.2 Number of FPs as a function of simulation time obtained by different potentials

表3 不同勢函數下級聯到達熱峰和平衡穩定狀態的時間Table 3 Time for cascade simulation with different potentials to reach peak state and equilibrium stable state

紅色和藍色小球分別代表鎢間隙原子和空位a，c——熱峰狀態；b，d——最終穩定狀態圖3 級聯碰撞兩種典型的缺陷分布Fig.3 Two typical defects distributionswith collision cascade

此外，注意到級聯熱峰階段缺陷分布的空間形貌主要有兩種：一種是集中型，其熱峰狀態構型如圖3a所示，另外一種是較為分散的連續型次級聯結構，如圖3c所示。相對而言，熱峰階段為集中型時，隨著級聯的后續演化，間隙和空位的湮滅復合后，穩定狀態級聯中心可能會形成略大的空位核區域，由較大的間隙團簇和單間隙原子包圍，如圖3b所示，且最后穩定狀態下的FPs數目會略多于分散型形貌，缺陷團簇的尺寸也會更大，總體上可能會使其擁有更高的團簇分數。熱峰階段為較分散的連續型次級聯結構時，最后穩定狀態的缺陷多以單間隙或單空位以及小尺寸的團簇形式存在，很少能觀察到較大的缺陷團簇，如圖3d所示。這與Fu等[17]在鎢及WRe合金高能中子輻照中觀察到的結果一致，且此類現象在高能PKA級聯時會更顯著。

圖4示出各勢函數在PKA能量為50 keV下的級聯過后穩定狀態的缺陷構型。由圖4可見，間隙缺陷多以單間隙〈111〉啞鈴(dumbbell)形式存在，缺陷團簇數目不多，總體上間隙型缺陷團簇比空位型多。且PKA能量為50 keV下觀察到有位錯環和位錯線的出現，位錯環大多是由間隙型團簇形成。而在PKA能量為10 keV情況下，5種勢函數在穩定狀態下的間隙缺陷構型大多數為單間隙〈111〉啞鈴，團簇數目較少，可忽略。本工作中，定義包含2個及2個以上的凈缺陷數(間隙原子/空位)為缺陷團簇。在已有文獻[17，33-35]中，間隙為第3近鄰、空位為第2近鄰的缺陷團簇截斷距離判據已被廣泛使用。本文以勢函數Chen為例，計算了其穩定狀態下間隙與空位缺陷的徑向分布函數g(r)，如圖5所示。由圖5可見，級聯過后對于間隙原子，主峰在第3近鄰(NN3)處，而對于空位，主峰在第2近鄰(NN2)處?；谝陨辖Y果，本文中缺陷團簇截斷距離判據選取間隙為第3近鄰、空位為第2近鄰。

總體而言，盡管各勢函數的離位閾能略有差別[27]，但最后穩定狀態下的FPs數目的差別很小，且各勢函數FPs數目隨PKA能量的增加均顯著增加。

3 缺陷團簇分布

紅色和藍色小球分別代表鎢間隙原子和空位，間隙以啞鈴形式顯示圖4 勢函數PKA能量為50 keV下穩定狀態的缺陷構型Fig.4 Defect configuration in stable state obtained by different potentials at PKA energy of 50 keV

圖5 穩定狀態下間隙和空位的徑向分布函數Fig.5 Radial distribution function of interstitial and vacancy in stable state

表4 不同勢函數下缺陷團簇分布情況Table 4 Distribution of defect cluster obtained by different potentials

圖6示出各勢函數PKA能量為50 keV下的團簇分數。由圖6可見：所有勢函數的間隙團簇分數均大于空位團簇分數，這表明間隙更容易形成團簇；團簇分數較大的勢函數AT和MV2-B，在15次級聯模擬中熱峰階段缺陷主要呈現集中型形貌(所占比例分別為60%和80%)，與之相比，較為分散的次級聯型形貌出現概率略?。黄渌麆莺瘮?Chen、JW、MV4-S)熱峰階段呈現集中型形貌的概率為30%～40%，致使這些勢函數擁有略小的團簇分數，且在勢函數MV4-S中，間隙和空位的團簇分數相差不大。此現象可能與各勢函數的間隙遷移能有關[9-10]，勢函數AT、JW、MV2-B、MV4-S和Chen所計算的〈111〉啞鈴遷移能分別為0.03、0.01、0.04、0.03、0.13 eV?！?11〉單間隙擁有較小的遷移能會使得間隙缺陷具有更強的運動性，級聯后它們能快速逃離級聯中心區域，形成更廣的缺陷分布范圍，降低了與空位復合的概率和形成團簇(尤其是較大尺寸團簇)的可能性，最終使得缺陷團簇分數變小。勢函數JW的間隙遷移能比其他勢函數略小，其擁有最小的間隙團簇分數；而盡管勢函數Chen所計算的〈111〉啞鈴遷移能略大于其他勢函數，但它在本文所模擬的結果中卻未擁有最高的團簇分數，這與在目前總模擬次數中勢函數Chen在熱峰階段呈現集中型形貌的概率約為40%有關(低于勢函數AT和MV2-B)，繼續增加模擬次數可能會有所變化。

圖6 勢函數PKA能量為50 keV下的缺陷團簇分數Fig.6 Defect cluster fraction obtained by different potentials at PKA energy of 50 keV

圖7示出各勢函數PKA能量為50 keV下的團簇信息分布。圖7中5種勢函數的間隙和空位團簇均以小尺寸(包含2～20個間隙/空位)團簇為主。由圖7a可見，勢函數AT和MV2-B出現稍大尺寸間隙團簇的概率比勢函數Chen、JW和MV4-S的要大；圖7b表明，所有勢函數均以小尺寸的空位團簇為主，略大一些的空位團簇主要在勢函數AT和JW及MV4-S中觀察到。

從缺陷團簇分布方面來看，各勢函數的缺陷團簇均以小尺寸(包含2～20個間隙/空位)的團簇為主；各勢函數的間隙遷移能對其級聯過程中缺陷的空間分布范圍有影響，而缺陷的空間分布與團簇分數有關，熱峰時缺陷為集中型形貌時，最終穩定狀態下產生缺陷團簇的數目會更多且尺寸會更大，這種情況下獲得一個較大團簇分數的可能性更高。

圖7 勢函數PKA能量為50 keV下的缺陷團簇分布Fig.7 Defect cluster distribution obtained by different potentials at PKA energy of 50 keV

4 位錯信息分布

圖8示出各勢函數PKA能量為50 keV下的位錯信息分布。在模擬過程中，5種勢函數均觀察到有位錯線及位錯環的產生。根據文獻[26]所計算的各勢函數位錯環的形成能來看，勢函數Chen與第一性原理(DFT)中1/2〈111〉位錯環的形成能低于〈100〉位錯環形成能的趨勢一致，且結果最相符。勢函數AT和JW也與DFT結果有相似的趨勢，但這兩個勢函數1/2〈111〉和〈100〉位錯環之間的形成能差值相對較小，而勢函數MV2-B和MV4-S所計算的位錯環形成能的結果與DFT結果不一致。勢函數Chen、AT、JW和MV4-S 4種勢函數級聯后所產生的位錯環以1/2〈111〉間隙型位錯環為主，與實驗上觀察到的結果一致[36]。而勢函數MV2-B以〈100〉位錯環出現較多，1/2〈111〉位錯環略少，這與勢函數MV2-B所預測的位錯環形成能為〈100〉位錯環低于1/2〈111〉位錯環有關。此外，勢函數AT和JW分出現了1次〈100〉間隙型位錯環和空位環，而勢函數Chen和MV4-S目前未觀察到〈100〉位錯環，這可能與模擬次數的限制有關，也可能與位錯環形成能有關。勢函數Chen所預測的〈100〉位錯環的形成能高于〈111〉位錯環，所以勢函數Chen的模擬結果觀察到的位錯環以〈111〉類型為主。在本文PKA能量為50 keV下模擬結果所觀察到的位錯環以小尺寸的為主，而MV4-S勢函數在小尺寸下(原子數<20)亦是1/2〈111〉位錯環的形成能低于〈100〉位錯環，隨位錯環尺寸的增大，則為1/2〈111〉位錯環的形成能高于〈100〉位錯環，所以勢函數MV4-S在小尺寸位錯環下可能會以〈111〉類型為主，而大尺寸位錯環更可能為〈100〉類型。

圖9示出各勢函數PKA能量為50 keV下位錯環的尺寸和數量分布。由圖9a可見，在勢函數Chen、AT和JW中占主導地位的為1/2〈111〉小尺寸位錯環，略大尺寸的1/2〈111〉位錯環主要在勢函數AT中出現。由圖9b可見：〈100〉位錯環主要在勢函數MV2-B中出現，且以間隙型位錯環為主；勢函數AT中觀察到1次小尺寸的間隙型〈100〉環，勢函數JW中觀察到1次大尺寸的空位〈100〉環。

從位錯信息分布來看，勢函數Chen、AT、JW和MV4-S 4種勢函數級聯后所產生的位錯以伯格斯矢量為1/2〈111〉間隙型位錯環為主，與實驗結果一致；而勢函數MV2-B以〈100〉位錯環出現較多，1/2〈111〉位錯環略少，與實驗結果略有出入。這與各勢函數預測的不同類型的位錯環形成能大小有關，與在第一性原理計算中1/2〈111〉位錯環的形成能低于〈100〉位錯環的結果相比，在這5個勢函數中，勢函數Chen的計算結果與第一性原理計算結果最相符。

結果為15次重復模擬統計的總和a——1/2〈111〉位錯環；b——〈100〉位錯環圖9 勢函數PKA能量為50 keV下位錯環的尺寸和數量分布Fig.9 Size and number distributions of dislocation loop obtained by different potentials at PKA energy of 50 keV

5 小結

本文對5個典型的鎢勢函數進行了PKA能量在10 keV和50 keV下的中子輻照級聯碰撞模擬，系統分析和討論了級聯碰撞過程中缺陷的產生與分布、缺陷團簇和位錯環的數目與結構等信息，獲得的結果對于鎢基材料初級輻照損傷的理解以及退火過程缺陷長時間的演化模擬提供了基礎，為用于輻照模擬鎢勢函數的選擇和優化提供了參考。得到的主要結論如下：

1) 對于級聯碰撞到達熱峰狀態的時間及平衡穩定狀態下FPs數目，不同勢函數的模擬結果沒有明顯差別。

2) 在所模擬的PKA能量下，不同勢函數產生的缺陷團簇均以小尺寸(凈缺陷數<20)為主，對于勢函數AT和MV2-B，缺陷團簇分數比其他勢函數較高，且較易出現大尺寸缺陷團簇，這與其熱峰階段易呈現集中型缺陷形貌有關。

3) 從位錯環的分布來看，勢函數Chen、AT、JW和MV4-S級聯碰撞模擬后所產生的位錯環以1/2〈111〉間隙型位錯環為主，與實驗結果一致，MV2-B勢函數出現的〈100〉間隙位錯環比1/2〈111〉間隙位錯環略多，與實驗結果不符。