不同溫度下α-Fe中級聯(lián)碰撞分子動力學(xué)模擬研究

曹 晗，賀新福，王東杰，吳 石，賈麗霞，豆艷坤，楊 文

(中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

BCC鐵基材料(α-Fe)是反應(yīng)堆中的重要結(jié)構(gòu)材料，低溫輻照脆化是這類材料服役面臨的一個主要問題[1]。輻照導(dǎo)致材料中形成的基體損傷是引起輻照脆化的關(guān)鍵因素之一。在輻照材料中，基體損傷起源于入射粒子與晶格點陣原子相互作用引發(fā)的級聯(lián)碰撞[2]，此過程會產(chǎn)生大量的缺陷。由于級聯(lián)碰撞過程發(fā)生在原子尺度，并在約10 ps時間內(nèi)完成，實驗上無法直接觀察，在此時間和空間尺度下，分子動力學(xué)方法是最有效、應(yīng)用最廣泛的一種模擬級聯(lián)碰撞過程中缺陷產(chǎn)生與演化的方法。與此同時，級聯(lián)碰撞的分子動力學(xué)模擬對于輻照機(jī)理認(rèn)識，如缺陷退火機(jī)制、位錯環(huán)形核長大機(jī)制，具有重要的指導(dǎo)意義，并能為動力學(xué)蒙特卡羅(kinetic Monte Carlo, KMC)和速率理論(rate theory, RT)等其他微觀結(jié)構(gòu)演化模擬方法提供如缺陷的類型、空間分布、存活率和成團(tuán)率等必要的參數(shù)，其中，存活率為缺陷存活數(shù)量NFP占離位原子數(shù)量Nd的比例，成團(tuán)率PC為間隙原子/空位團(tuán)簇中含有的缺陷數(shù)量NC占缺陷存活數(shù)量的比例。

自20世紀(jì)90年代以來，利用分子動力學(xué)模擬方法研究α-Fe的級聯(lián)碰撞一直是輻照損傷模擬的熱點問題：王寧等[3]通過對α-Fe的級聯(lián)碰撞模擬，得到了自間隙原子的取向與遷移規(guī)律；楊莉等[4]研究了He參與下α-Fe在低溫時的級聯(lián)碰撞過程，分析得到了He泡的形成機(jī)制；Peng等[5]研究了100 keV以上初級離位原子(PKA)引發(fā)的級聯(lián)碰撞，并發(fā)現(xiàn)了位錯環(huán)的生成，提出了〈100〉間隙型位錯環(huán)形成的punch-out機(jī)制。然而，截至目前，在α-Fe的級聯(lián)碰撞模擬研究中，依然存在著亟待解決的問題：業(yè)界一直嘗試在級聯(lián)碰撞產(chǎn)生的缺陷存活數(shù)量與PKA能量之間建立準(zhǔn)確的數(shù)量關(guān)系，1955年，Kinchin和Pease建立Kinchin-Pease模型[6]，首次預(yù)測了離位原子數(shù)量與PKA能量的數(shù)量關(guān)系，經(jīng)Norgett等[7]于1975年進(jìn)行修正后，建立了NRT模型；對于缺陷存活數(shù)量，國內(nèi)外多項研究顯示[8-10]，級聯(lián)碰撞產(chǎn)生的穩(wěn)定缺陷數(shù)量與PKA能量并不是NRT模型的正比例關(guān)系，而是冪次關(guān)系；然而，有研究顯示，由于PKA能量達(dá)20 keV以上時將會發(fā)生子級聯(lián)[11]，導(dǎo)致不能用同一個擬合公式描述全部能量范圍內(nèi)存活缺陷數(shù)量與PKA能量的數(shù)值關(guān)系，不同PKA能量范圍的模型存在區(qū)別的機(jī)理尚不明確，這種區(qū)別也會極大地影響模型的準(zhǔn)確性和簡潔性。

本文將利用分子動力學(xué)方法進(jìn)行一系列不同溫度、PKA能量條件的α-Fe級聯(lián)碰撞模擬研究，研究級聯(lián)碰撞后產(chǎn)生的缺陷構(gòu)型和分布，擬為團(tuán)簇、位錯環(huán)等缺陷在微觀尺度上演化的RT和KMC等模擬方法提供必要的參數(shù)，以確定溫度與不同PKA能量范圍下NFP-EPKA擬合曲線存在區(qū)別的相關(guān)性，為修正該模型奠定基礎(chǔ)，進(jìn)而完善反應(yīng)堆用BCC鐵基材料的多尺度模擬。

1 計算方法

本工作采用Marinca等在2007年建立的M07純鐵勢函數(shù)[12-13]，該勢函數(shù)與常用的M03和A04兩種勢函數(shù)相比，模擬結(jié)果能更準(zhǔn)確地與第一性原理以及實驗結(jié)果符合[12]。利用分子動力學(xué)軟件Large-scale Atomic/Molecular Massively Parallel Simulator(LAMMPS)[14]研究了500、600和700 K 3個溫度條件下BCC純鐵材料的級聯(lián)碰撞過程，模擬盒子尺寸根據(jù)PKA能量大小分為80a0×80a0×80a0與150a0×150a0×150a0兩種，其中a0表示晶格常數(shù)。計算過程中采用三維周期邊界條件，在級聯(lián)碰撞開始前將系統(tǒng)升溫至模擬輻照溫度并保溫20～50 ps以確保系統(tǒng)達(dá)到熱穩(wěn)定狀態(tài)。級聯(lián)碰撞過程須在模擬初始狀態(tài)下引入PKA，PKA的能量與模擬的盒子內(nèi)原子數(shù)量和模擬時間的具體對應(yīng)關(guān)系列于表1。為避免產(chǎn)生溝道效應(yīng)，并盡量使入射方向隨機(jī)分布，根據(jù)前人研究經(jīng)驗[11]，本工作針對每種PKA能量分別沿〈122〉、〈135〉和〈235〉 3個方向引入PKA。為使模擬結(jié)果達(dá)到統(tǒng)計效果，每種模擬條件計算3次。引入PKA后，每隔0.05 ps輸出1次所有原子的位置信息，0.5 ps后，每隔0.25 ps輸出1次，5.5 ps后，每隔0.5 ps輸出1次，用以分析級聯(lián)碰撞過程中缺陷的演化過程。為保證級聯(lián)碰撞充分退火，根據(jù)PKA能量的不同，模擬總時長也有區(qū)別，PKA能量低于40 keV時，模擬時間為15 ps，PKA能量為40 keV和50 keV時，模擬時間分別為30 ps和40 ps。

表1 不同PKA能量的級聯(lián)碰撞模擬盒子內(nèi)原子數(shù)量Table 1 Number of atoms in displacement cascadessimulation boxes with different PKA energy

級聯(lián)過程模擬完成后，采用Open Visualization Tool(OVITO)軟件中的Wigner-Seitz晶胞方法判斷間隙原子和空位[15]，缺陷類型根據(jù)晶胞內(nèi)原子的數(shù)量判定。對級聯(lián)碰撞模擬結(jié)果進(jìn)行缺陷團(tuán)簇分析，采用的間隙團(tuán)簇截斷半徑為第二近鄰距離，空位團(tuán)簇截斷半徑為第三近鄰距離，進(jìn)而利用式(1)計算缺陷的成團(tuán)率。

(1)

利用OVITO軟件中的dislocation analysis(DXA)方法對模擬完成后產(chǎn)生的團(tuán)簇進(jìn)行位錯環(huán)分析。

2 結(jié)果與討論

2.1 級聯(lián)碰撞模擬過程

一次完整的由中子引發(fā)的級聯(lián)碰撞包括4個過程，分別為缺陷的產(chǎn)生、離位峰的出現(xiàn)、缺陷的湮滅以及退火形成穩(wěn)定缺陷[2]。圖1為600 K時由30 keV的PKA在含有1 024 000個原子的α-Fe中導(dǎo)致的級聯(lián)碰撞缺陷演化過程可視化截圖。從圖中可看出，30 keV的PKA形成后，在0.1 ps的時間內(nèi)便使PKA附近產(chǎn)生大量點缺陷；此算例在1 ps前后，F(xiàn)renkel缺陷數(shù)量達(dá)到峰值；之后，經(jīng)過14 ps的缺陷 復(fù)合，體系內(nèi)最終產(chǎn)生了3個1/2〈111〉間隙型位錯環(huán)、1個〈100〉空位型位錯環(huán)、不同尺寸的間隙和空位團(tuán)簇以及多個分布于基體中的點缺陷。在此過程中，由于空位的遷移能較間隙原子的遷移能高[16]，導(dǎo)致空位型缺陷不易遷移，所以空位型的位錯環(huán)、團(tuán)簇和單空位缺陷較間隙型缺陷更靠近PKA區(qū)域，從而呈現(xiàn)出間隙原子包圍著空位的現(xiàn)象。在此算例中，產(chǎn)生的間隙原子快速擴(kuò)散逃離了級聯(lián)區(qū)域，而空位由于擴(kuò)散能力較差，大部分仍留在級聯(lián)區(qū)域，從而導(dǎo)致該區(qū)域形成空位型位錯環(huán)[17]。

圖1 600 K時30 keV PKA在α-Fe中導(dǎo)致的級聯(lián)碰撞過程中產(chǎn)生的缺陷演化過程可視化截圖Fig.1 Visual screenshot of defect evolution during displacement cascades induced by 30 keV PKA in BCC iron at 600 K

2.2 缺陷構(gòu)型

對于級聯(lián)碰撞后產(chǎn)生的穩(wěn)定缺陷(圖2)，主要分為3類。1) 彌散分布于基體中的點缺陷，包括間隙原子和空位兩種，間隙原子以啞鈴狀的自間隙原子(dumbbell)存在，如圖2b所示，對dumbbell中的兩原子位置坐標(biāo)進(jìn)行分析，考慮到熱振動會導(dǎo)致原子的位置發(fā)生微小偏移，可判斷出dumbbell為〈110〉取向，點缺陷的遷移能相對較低，在長時間(s量級以上)的演化中，易相互結(jié)合或擴(kuò)散至晶界、位錯、表面等界面處而被湮滅，實際上，中子的不斷入射和點缺陷的湮滅會使基體中的點缺陷濃度維持在較穩(wěn)定的水平。2) 由少量間隙原子或空位組成的相應(yīng)類型的小型團(tuán)簇，間隙型小團(tuán)簇中最近鄰原子的排列方向與啞鈴狀自間隙原子一致，如圖2c所示；而在空位型小團(tuán)簇中兩個最近鄰空位的排列方向有〈100〉和〈111〉兩種，如圖2d所示。3) 由大量間隙原子或空位組成的位錯環(huán)，在所有的算例中，主要形成了兩種類型的位錯環(huán)，分別為1/2〈111〉間隙型位錯環(huán)和〈100〉空位型位錯環(huán)，在間隙型位錯環(huán)中，如圖2e統(tǒng)計顯示，兩個最近鄰原子的相對方向大部分為〈111〉，少量為〈112〉和其他方向，這是因為自間隙原子大量聚集時，〈111〉方向的dumbbell組成的團(tuán)簇結(jié)合能更高[18]；在空位型位錯環(huán)中，兩個最近鄰空位的相對方向則大多為〈100〉方向。

a——600 K時30 keV PKA在α-Fe中引發(fā)的級聯(lián)碰撞產(chǎn)生的缺陷整體形貌；b——dumbbell中兩原子取向矢量(由最近鄰兩原子位置坐標(biāo)相減得到，余同)；c——間隙型小團(tuán)簇中最近鄰原子取向矢量；d——空位型小團(tuán)簇中最近鄰空位取向矢量；e——1/2〈111〉間隙型位錯環(huán)中最近鄰原子取向分布柱形圖；f——空位型位錯環(huán)構(gòu)型圖2 級聯(lián)碰撞后產(chǎn)生的缺陷Fig.2 Defects produced by displacement cascades

2.3 缺陷存活數(shù)量

圖3定量給出3種溫度下不同PKA能量級聯(lián)碰撞產(chǎn)生的缺陷數(shù)隨時間的演化曲線。經(jīng)分析得知，在誤差允許范圍內(nèi)，隨著PKA能量或溫度的升高，出現(xiàn)缺陷數(shù)量峰值以及最終達(dá)到穩(wěn)定狀態(tài)所需的時間均會變長。圖4為3種溫度下不同PKA能量引起的級聯(lián)碰撞導(dǎo)致的缺陷存活數(shù)量圖。

圖3 不同溫度下級聯(lián)碰撞缺陷數(shù)隨時間的演化Fig.3 Evolution of number of displacement cascades’ defects with time at different temperatures

圖4 不同溫度下級聯(lián)碰撞缺陷存活數(shù)量隨PKA能量變化Fig.4 Number of survival defects at different temperatures vs energy of PKA

1975年，Norgett等對Kinchin-Pease模型進(jìn)行了修正，建立了式(2)所示的NRT模型，定量描述了在級聯(lián)碰撞中產(chǎn)生的離位原子數(shù)量與PKA能量的數(shù)量關(guān)系。

(2)

其中，Ed為靶材料的離位閾能。經(jīng)多年的模擬研究，在NRT模型的基礎(chǔ)上，業(yè)界普遍認(rèn)為級聯(lián)碰撞產(chǎn)生的穩(wěn)定缺陷數(shù)量與PKA能量并不是NRT模型的正比例關(guān)系，而是冪次關(guān)系，即：

(3)

其中，a、b為系數(shù)，b取0.8左右，a和b均由級聯(lián)碰撞模擬后得到的缺陷數(shù)量和PKA能量擬合公式得出。

圖5 EPKA=30 keV前后擬合曲線中的bFig.5 b in fitting curve before and after EPKA=30 keV

2.4 團(tuán)簇分析

圖6為3種溫度下不同PKA能量引起的級聯(lián)碰撞導(dǎo)致的間隙和空位團(tuán)簇尺寸、數(shù)量三維柱形圖。從圖可看出，PKA能量越高，出現(xiàn)大尺寸團(tuán)簇的概率越大，團(tuán)簇總數(shù)量越多；溫度的升高則會使低能PKA導(dǎo)致的級聯(lián)碰撞發(fā)生間隙團(tuán)簇的形核長大行為變得容易，空位團(tuán)簇則無明顯區(qū)別。圖7展示了3種溫度下不同PKA能量引起的級聯(lián)碰撞導(dǎo)致的間隙原子和空位的成團(tuán)率?？煽闯?，在誤差允許的范圍內(nèi)，PKA能量的升高會使間隙原子和空位的成團(tuán)率增大，溫度的升高會使這種增大的趨勢變緩，即在700 K時缺陷的成團(tuán)率基本不隨PKA能量變化，分別穩(wěn)定于50%和40%左右，因此，在PKA能量變化時，高溫時b的變化量與低溫時的相比較小。

a——T=500 K，間隙團(tuán)簇；b——T=600 K，間隙團(tuán)簇；c——T=700 K，間隙團(tuán)簇；d——T=500 K，空位團(tuán)簇；e——T=600 K，空位團(tuán)簇；f——T=700 K，空位團(tuán)簇圖6 級聯(lián)碰撞間隙和空位團(tuán)簇尺寸和數(shù)量三維柱形圖Fig.6 3D cylindrical map of size and quantity of interstitial and vacancy clusters of displacement cascades

圖7 不同溫度下間隙原子和空位成團(tuán)率Fig.7 Cluster fraction of interstitials and vacancies for different temperatures

3 結(jié)論

通過對α-Fe在500、600和700 K下的級聯(lián)碰撞進(jìn)行一系列的分子動力學(xué)模擬研究，得到以下結(jié)論。

1) 級聯(lián)碰撞會使體系內(nèi)產(chǎn)生大量輻照缺陷，包括點缺陷、間隙和空位團(tuán)簇等，缺陷的尺寸決定了缺陷內(nèi)最近鄰原子的排列方向：間隙型小團(tuán)簇中最近鄰原子的排列方向與啞鈴狀自間隙原子一致，均為〈110〉方向；間隙型位錯環(huán)中兩個最近鄰原子的相對方向大部分為〈111〉。

2) 高能PKA使得產(chǎn)生的位錯環(huán)級別的大尺寸團(tuán)簇的概率升高，從而有更多數(shù)量的缺陷存活，導(dǎo)致同一溫度下不同PKA能量階段的NFP-EPKA擬合曲線的b存在差別。

3) 在分段NFP-EPKA擬合曲線中，700 K時，EPKA=30 keV前后兩階段b相差最小，原因在于：溫度的升高會使低能PKA導(dǎo)致的級聯(lián)碰撞發(fā)生間隙團(tuán)簇的形核長大行為變得容易，使得700 K時缺陷的成團(tuán)率基本不隨PKA能量變化，從而導(dǎo)致不同PKA能量引發(fā)的級聯(lián)碰撞產(chǎn)生的缺陷數(shù)量增長趨勢未發(fā)生明顯變化。