鋯自輻射損傷的分子動力學(xué)模擬與缺陷判定研究

鐘 睿, 付寶勤, 崔節(jié)超, 侯 氫

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

1 引 言

由于中子和鋯靶原子的彈性碰撞截面很小,發(fā)生相鄰兩次彈性碰撞事件距離在厘米量級,而級聯(lián)碰撞影響區(qū)域在納米量級， 所以載能中子致材料輻射損傷的模擬一般分為兩個獨立的模擬過程[1]. 首先是中子在材料中輸運的模擬計算,以獲得反沖原子能譜及其空間分布. 就如崔振國等[1]在其研究中指出的,1 MeV的中子在鋯材料的靶材中產(chǎn)生的大部分PKA的能量都分布在1～15 keV之間,不同PKA的空間分布大多在6～8 mm之間. 其次再根據(jù)反沖原子能譜抽樣計算各反沖原子能量下的級聯(lián)碰撞過程. 本文涉及的是后一個過程,即級聯(lián)碰撞過程的分子動力學(xué)模擬計算.

目前模擬原子級聯(lián)碰撞過程最廣泛采用的方法是以二體碰撞理論為基礎(chǔ)的蒙特卡洛模擬方法[2-4]. 相應(yīng)的代表性程序包有廣為人知的SRIM等. 不過這種基于兩體碰撞近似的模擬不能考慮級聯(lián)碰撞過程中空位對離位原子的再捕獲,因此對級聯(lián)碰撞過程更精確仔細的模擬需要采用分子動力學(xué)模擬. 雖然隨著計算技術(shù)的發(fā)展分子動力學(xué)已被越來越多地應(yīng)用于材料損傷過程的研究[5-17],但如何用分子動力學(xué)模擬足夠樣本量的級聯(lián)碰撞過程,并從中自動提取需要的統(tǒng)計物理量仍是值得仔細研究的問題.

與BCC等結(jié)構(gòu)的金屬相比較,關(guān)于鋯級聯(lián)碰撞的分子動力學(xué)模擬報道相對較少[18-22],其中比較有代表意義的包含,1998年Wooding等[19]采用Finnis-Sinclair勢函數(shù)研究了PKA能量En達到20 keV時在Zr中產(chǎn)生的損傷,并指出Frenkel-pairs的數(shù)量和En呈一種冪相關(guān)性. 2001年Gao等[20]采用Ackland勢函數(shù)研究了溫度在100～600 K范圍內(nèi)Frenkel-pairs數(shù)量以及自間隙原子團簇隨溫度的變化情況. Voskoboinikov等[22]還研究了Zr金屬Frenkel-pairs的數(shù)量和團簇形式與溫度和PKA能量的相關(guān)性. 文獻[1]報道了鋯中級聯(lián)碰撞的分子動力學(xué)模擬結(jié)果,但是PKA能量限于10 keV以下能區(qū). 另外,文獻[1]在計算離位原子數(shù)方面,靶原子相對其初始位置的位移量超過一個閾值則計入離位原子,由于沒有考慮這些靶原子可以被其他離位原子留下的空位捕獲,因此會高估離位原子數(shù).

本文采用分子動力學(xué)模擬了10～50 keV能區(qū)的初級反沖原子在HCP結(jié)構(gòu)Zr中的原子級聯(lián)碰撞,并從方法學(xué)層面著重于考察了三種方法對點缺陷判定的影響. 結(jié)果表明,一種以Wigner-Seitz(WS)原胞分割為基礎(chǔ)的判定方法更適合用于Zr在初級輻射損傷中產(chǎn)生的缺陷數(shù)統(tǒng)計. 該結(jié)論對于如何精確評估中子在Zr中產(chǎn)生的輻照損傷有指導(dǎo)意義. 此外,本文還研究了根據(jù)WS分割方法確定的自間隙原子點缺陷的遷移特征.

2 物理模型與分析方法

2.1 物理模型

本文采用分子動力學(xué)方法模擬Zr原子自輻射產(chǎn)生的級聯(lián)碰撞,并通過統(tǒng)計級碰撞產(chǎn)生的點缺陷數(shù)來研究Zr受到輻射后的損傷程度. 其中,PKA能量范圍設(shè)置為10～50 keV. 針對不同PKA能量,模擬體系采用不同原子數(shù)的盒子以提高計算效率. 盒子大小的選擇依據(jù)是,級聯(lián)碰撞主要影響區(qū)域不能因周期性邊界條件出現(xiàn)重疊(擊穿盒子),否則會造成晶體缺陷數(shù)減少[12]. 研究中,10～30 keV條件下模擬盒子包含原子數(shù)為414 720、 810 000和1 399 680;40和50 keV條件下包含原子數(shù)為6 480 000.

觀察組：男、女占比各為26：19；年齡段在49～83歲之間，經(jīng)計算后中位年齡為（66.21±1.34）歲。

分子動力學(xué)模擬采用Mendelev提出的Zr-Zr勢函數(shù)[23],并采用了經(jīng)過GPU算法優(yōu)化的MD程序包MDSCU[24]以提高模擬計算速度. 每種PKA能量并行模擬5次. 模擬盒子溫度設(shè)置為300 K,熱量耗散采用聲電耦合(E-P coupling， EPC)方法[25]. PKA隨機選擇入射角度. MD模擬采用自適應(yīng)時間步長,保證運動最快的原子在一個時間步內(nèi)的位移不能超過0.01a0(a0為鋯晶格常數(shù)). 同時,每1 000步記錄并輸出一次分子動力學(xué)模擬結(jié)果.

此外,自間隙原子擴散研究方面,采用6 481個原子的盒子,并行模擬1 000次,模擬時間30 ps,記錄時間間隔0.5 fs.

2.2 分析方法

根據(jù)分子動力學(xué)模擬產(chǎn)生的原子運動軌跡,本文比較了3種判定和計算間隙原子及空位的方法.

方法一,離位原子判定方法. 離位原子判定方法在蒙特卡洛模擬中使用較多,用于評估不同能量PKA造成影響的程度. 具體如圖1a所示. 圖1a中標(biāo)記“1”和“2”分別表示一個原子的初始位置和當(dāng)前位置,R表示位移距離. 當(dāng)R值超過設(shè)定閾值,則該原子就被判定為離位原子,而并不考慮空位對間隙原子的捕獲情況. 本研究中設(shè)定R=1a0,a0為Zr晶體晶格常數(shù).

圖1 點缺陷判定方法Fig.1 Point defects identification methods

圖方向可視化結(jié)果(a=0.1 ps; b=1 ps; c=3 ps; d=10 ps)Fig.2 Visual result of cascade, direction (a=0.1 ps; b=1 ps; c=3 ps; d=10 ps)

方法二,傳統(tǒng)間隙原子及空位判定方法. 該方法以計算原子和晶格格點之間的距離為基礎(chǔ),具體實現(xiàn)如圖1b所示. 圖1b中小圓環(huán)為晶體格點位置,實心圓點為原子位置. 以每個格點位置為球心設(shè)置半徑為R的范圍,若原子落入此范圍內(nèi),被判定為普通原子,如標(biāo)識“1”所示. 若原子不落入任何一個格點的R范圍內(nèi),則被判定為間隙原子,如標(biāo)識“2”所示. 若一個格點的R范圍內(nèi)無任何原子,則該格點位置被判定為空位,如標(biāo)識“3”所示. 根據(jù)算法描述,通過該方法獲得的間隙原子和空位數(shù)不一定匹配,即不一定都能構(gòu)成Frenkel-pairs. 此外,采用該方法還需要考慮閾值R設(shè)定對判定結(jié)果的影響.

方法三,Wigner-Seitz (W-S)原胞分割方法. W-S原胞在空間上可以認為是以格點位置為中心點的Voronoi分割. 完備晶格中每個原胞中包含一個格點,若某原子在此原胞中,則該原子距離所在原胞格點的距離比其它原胞都短. 以此為依據(jù),該判定方法如圖1c所示. 圖中小圓環(huán)為格點位置,實心圓點為原子位置. 判定時首先計算每個原子和所有格點的距離,設(shè)定每個原子“歸屬”距離其最近的格點,并判定如下:(1)若只有一個原子“歸屬”某個格點,判定為正常情況. 如圖中標(biāo)識“c”所示原子和“L3”所示格點位置關(guān)系. (2)若多于一個原子“歸屬”某個格點,則該格點存在間隙原子. 本文定義此格點位置為LSIA(location of self-interstitial atoms). 如圖中標(biāo)識“a”,“b”所示原子和“L1”所示格點位置關(guān)系,“L1”即LSIA. (3)若沒有任何原子“歸屬”某個格點,則該格點位置被判定為空位. 如圖中“L2”所示格點位置.

3 模擬結(jié)果

圖2為通過可視化軟件觀察到的PKA能量為10 keV時級聯(lián)碰撞的一個例子.

圖2a為級聯(lián)碰撞的初始階段,只有少量原子被碰撞出來,形成離位原子;圖2b為隨著級聯(lián)碰撞發(fā)展,更多原子成為離位原子,并逐步達到離位峰階段;圖2c表示部分間隙原子重新被空位捕獲,晶體缺陷數(shù)減少階段;圖2d表示模擬系統(tǒng)恢復(fù)穩(wěn)定,缺陷進入由溫度驅(qū)動的自由演化階段. 整個模擬過程中點缺陷數(shù)的變化呈現(xiàn)迅速升高到峰值、逐步減少、基本穩(wěn)定等三個過程,這與通常對級聯(lián)碰撞過程的認識一致.

上述過程還可以通過模擬系統(tǒng)溫度場的變化得到印證. 本研究截取了級聯(lián)碰撞最活躍的區(qū)域,通過不同顏色標(biāo)記每一個原子的動能范圍,繪制了該區(qū)域溫度場的變化情況,如圖3所示.

圖3 溫度場(a=0.1 ps; b=1 ps; c=10 ps)Fig.3 Temperature field (a=0.1 ps; b=1 ps; c=10 ps)

圖3中,高能原子用紅色表示,低能原子用藍色表示. 級聯(lián)碰撞之初,存在大量高能原子,這些原子和周圍原子發(fā)生碰撞時,新原子獲得的能量易超過離位閾能,從而被激發(fā)出來形成離位原子并使得級聯(lián)碰撞可以持續(xù)下去(如圖3a所示);離位峰出現(xiàn)后,大部分原子的能量已經(jīng)降低,不足以激發(fā)出更多離位原子原子(如圖3b所示);模擬系統(tǒng)恢復(fù)平衡后,原子基本恢復(fù)到模擬盒子本底溫度,只在熱驅(qū)動下演化(如圖3c所示).

本文采用2.2節(jié)中介紹的3種方法分別對MD模擬結(jié)果進行了點缺陷判定. 以PKA能量10 keV為例,點缺陷的變化趨勢對比如圖4所示.

圖4 點缺陷演化比較(a:方法1; b:方法2; c:方法3)Fig.4 Point defects evolution comparison(a: method 1; b: mehod 2; c: method 3)

圖5 統(tǒng)計漲落及魯棒性比較(a:方法1; b:方法2; c:方法3)Fig.5 Fluctuation and robustness comparison (a: method 1; b: mehod 2; c: method 3)

由圖可見,離位原子數(shù)變化呈上升趨勢,盡管模擬系統(tǒng)恢復(fù)穩(wěn)定后變化趨緩,但總體未出現(xiàn)減少過程,這表明方法一不能體現(xiàn)間隙原子被空位重新捕獲的情況. 方法二和方法三中,盡管峰值以及穩(wěn)定后的缺陷數(shù)不同,但變化趨勢總體一致,和圖2結(jié)論自洽,可以完整反映級聯(lián)碰撞各個階段缺陷數(shù)的變化情況.

模擬系統(tǒng)恢復(fù)穩(wěn)定后,缺陷在溫度驅(qū)動下自由演化,此時的缺陷數(shù)可以用于反映晶體的損傷程度. 此階段3種判定方法獲得的陷數(shù)統(tǒng)計漲落以及判定方法魯棒性的對比如圖5所示.

由圖5可見, 5次分子動力學(xué)模擬獲得的離位原子數(shù)統(tǒng)計漲落最大,在約800～1 600;方法二獲得的點缺陷數(shù)統(tǒng)計漲落也較大,在約60～120之間;W-S分割方法獲得的點缺陷數(shù)統(tǒng)計漲落最小,在19～33之間. 統(tǒng)計漲落越小,意味著獲得較準(zhǔn)確平均值所需的模擬次數(shù)越少,可以相應(yīng)節(jié)約計算開銷.

判定方法的魯棒性比較方面,方法一和方法二獲得的缺陷數(shù)都因缺陷在溫度驅(qū)動下的自由演化而出現(xiàn)明顯振蕩,而W-S分割方法獲得的缺陷數(shù)非常穩(wěn)定,很少出現(xiàn)變化,表現(xiàn)出良好的魯棒性. 魯棒性越好,意味著單次模擬所需的時間越短,同樣可以節(jié)約計算開銷.

方法一和方法二中需要設(shè)置判定閾值,本文還研究了閾值對分析結(jié)果的影響. 以方法二為例,圖6顯示了R=0.3a0、R=0.35a0和R=0.4a0時獲得的點缺陷數(shù).

圖6 閾值R影響

圖6a為10 ps以前模擬結(jié)果,圖6b為10 ps以后模擬結(jié)果. 整個過程中,采用不同R值進行判定獲得的點缺陷數(shù)差距很大. 這表明以距離閾值為判定依據(jù)的缺陷數(shù)分析方法對閾值的設(shè)定非常敏感.

綜上所述,方法一單純計算原子和初始位置的距離,未考慮原子和其它格點位置的關(guān)系,因此無法識別間隙原子被其它空位再捕獲的情況;方法二本質(zhì)是一種范圍判定,它可以識別間隙原子被其它空位再捕獲的情況,但是如果由于R設(shè)置的敏感性造成部分間隙原子可能落在R范圍內(nèi),這些間隙原子將被識別為普通原子,從而造成間隙原子和空位數(shù)不匹配;方法三采用了相對距離判定,無需設(shè)置閾值參數(shù),而且可以保證間隙原子和空位數(shù)一定相等. 因此,從3種判定方法獲得的晶體缺陷數(shù)變化趨勢、統(tǒng)計漲落;方法的魯棒性、對閾值的敏感性等多方面評價,采用W-S原胞分割判定方法獲得的點缺陷數(shù)更適合準(zhǔn)確反映晶體受到輻射后的損傷程度.

本文還對10、20、30、40和50 keV PKA能量的分子動力學(xué)模擬結(jié)果進行了分析. 采用W-S方法獲得的點缺陷數(shù)變化如圖7所示.

圖7 不同PKA能量產(chǎn)生的缺陷數(shù)Fig.7 Defects number in different PKA energies

不同PKA能量條件下缺陷數(shù)演化過程一致. 缺陷數(shù)和PKA能量基本成正比,缺陷數(shù)峰值出現(xiàn)在1 ps附近,模擬體系恢復(fù)平衡在10 ps附近,以上兩個時刻基本不隨PKA能量改變而發(fā)生變化.

圖8 300 K時LSIA擴散軌跡Fig.8 Typical trajectories of LSIAs for temperature T=300 K

表1 ILJ,OLJ,OPJ比例

由表可見,OLJ類跳躍是主要的跳躍方式. 因此LSIA更傾向于2D擴散.

4 結(jié) 論

本文采用分子動力學(xué)方法模擬了5種PKA能量引起的Zr原子的級聯(lián)碰撞,同時采用了3種不同的點缺陷判定方法對MD模擬結(jié)果進行了分析. 其中,W-S原胞分割方法具有良好的魯棒性,不需要設(shè)置任何閾值參數(shù),多次模擬所獲的點缺陷數(shù)的統(tǒng)計漲落很小,更適合表征晶體的輻射損傷程度. 此外,分子動力學(xué)模擬結(jié)果表明,Zr級聯(lián)碰撞后的點缺陷數(shù)正比于PKA能量,離位峰出現(xiàn)在級聯(lián)碰撞開始后約1 ps左右,系統(tǒng)恢復(fù)穩(wěn)定出現(xiàn)在約10 ps左右. 以上兩個時刻基本不隨PKA能量變化而發(fā)生改變. 此外,SIA的擴散更傾向于是一種2D擴散.

事實上,正如我們熟知的,中子輻射損傷基本可以分為兩個階段,第一個階段是由原子級聯(lián)碰撞造成初級輻射損傷,這個階段發(fā)生在皮秒(ps)量級. 本文主要針對這一階段的物理過程以及點缺陷的判定方法進行研究. 第二個階段是由初級輻射損傷造成的缺陷在溫度驅(qū)動下的演化過程,包含缺陷的擴散、團簇的捕獲與解離等過程. 第二個階段的時間尺度遠遠大于第一個階段,可以采用分子動力學(xué)、KMC或者速率理論等方法進行. 需要指出的是,對第二階段模擬需要缺陷的初始狀態(tài),包括點缺陷數(shù)和分布等,則依賴于由第一個階段的模擬結(jié)果提供. 因此,本文關(guān)注缺陷判定方法的魯棒性等的目的都在于確保第一階段獲得的參數(shù)更為準(zhǔn)確,對下一步多尺度模擬的參數(shù)選擇具有一定的借鑒意義.