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

鐘 睿, 付寶勤, 崔節超, 侯 氫

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

1 引 言

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

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

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

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

2 物理模型與分析方法

2.1 物理模型

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

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

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

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

2.2 分析方法

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

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

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

圖方向可視化結果(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)

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

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

3 模擬結果

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

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

上述過程還可以通過模擬系統溫度場的變化得到印證. 本研究截取了級聯碰撞最活躍的區域,通過不同顏色標記每一個原子的動能范圍,繪制了該區域溫度場的變化情況,如圖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中,高能原子用紅色表示,低能原子用藍色表示. 級聯碰撞之初,存在大量高能原子,這些原子和周圍原子發生碰撞時,新原子獲得的能量易超過離位閾能,從而被激發出來形成離位原子并使得級聯碰撞可以持續下去(如圖3a所示);離位峰出現后,大部分原子的能量已經降低,不足以激發出更多離位原子原子(如圖3b所示);模擬系統恢復平衡后,原子基本恢復到模擬盒子本底溫度,只在熱驅動下演化(如圖3c所示).

本文采用2.2節中介紹的3種方法分別對MD模擬結果進行了點缺陷判定. 以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 統計漲落及魯棒性比較(a:方法1; b:方法2; c:方法3)Fig.5 Fluctuation and robustness comparison (a: method 1; b: mehod 2; c: method 3)

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

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

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

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

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

圖6 閾值R影響

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

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

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

圖7 不同PKA能量產生的缺陷數Fig.7 Defects number in different PKA energies

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

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

表1 ILJ,OLJ,OPJ比例

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

4 結 論

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

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