啞鈴型氧間隙缺陷在UO2材料內擴散的分子動力學研究

張華陽，王慶宇，李忠宇

(哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001)

二氧化鈾(UO2)是核反應堆中最常見的燃料，反應堆運行期間核燃料會受到裂變碎片、氣體擴散和熱應力應變等效應的影響，造成微觀結構的損傷并引起宏觀性能的變化，進而影響核反應堆的安全運行。所以，UO2材料的輻照損傷效應研究一直是核燃料研究的重要課題之一。用于研究輻照損傷效應的方法主要有兩種，一種是利用離子束輻照實驗，能研究材料輻照一定劑量后的微觀結構和宏觀性能變化，另一種是利用多尺度模擬計算，能研究材料微觀結構的連續變化過程。分子動力學方法是多尺度模擬計算中的一種重要方法，能模擬計算nm空間尺度和ps時間尺度內材料微觀結構的連續變化，揭示輻照損傷效應的機理。對于UO2材料，其自身的氧原子由于輻照效應，有概率離開原來的平衡位置，在其他位置形成氧間隙原子，氧間隙原子的遷移將影響UO2材料的微觀結構以及宏觀性能。UO2材料內間隙原子的擴散現象一直是核燃料輻照效應的研究熱點。賀新福等[1]曾研究過氦間隙原子的擴散現象，關于氧間隙原子的研究[2-4]主要是在沒有應變的條件下研究擴散現象。本文利用分子動力學方法研究應變條件下啞鈴型氧間隙缺陷的擴散現象，為實驗研究提供理論依據。

1 方法

本文的模擬計算均利用分子動力學軟件LAMMPS[5]進行，并結合軟件OVITO[6]對模擬計算結果進行分析。LAMMPS是大規模原子分子并行模擬器，能模擬計算百萬級的原子或分子系統的微觀結構演化，OVITO是一款可視化分析軟件，通過讀取LAMMPS的計算結果能顯示每個原子或分子的運動情況，并能對材料的微觀結構進行相關分析。對于UO2材料體系，本文利用Yakub等[7]開發的分數離子模型勢函數描述體系中各原子之間的相互作用力，此勢函數能精確描述原子間的相互作用以及每個原子的運動情況。

模擬計算需設定UO2材料體系的晶格常數、弛豫溫度、拉伸方向、拉伸應變以及時間步長和步數等參數。首先由LAMMPS軟件建立UO2材料體系，接著在正則系綜(NVT)下進行充分弛豫達到熱力學平衡，弛豫溫度分別為1 000、1 300、1 500、1 800 K。圖1為1 300 K弛豫溫度下的UO2體系模型。其中，圖1a中z軸方向對應UO2材料體系的〈100〉晶向，尺寸為4.42 nm×4.42 nm×13.27 nm，包含6 144個U原子和12 288個O原子；圖1b中z軸方向對應UO2材料體系的〈110〉晶向，尺寸為

圖1 達到熱力學平衡時UO2材料體系初始模型Fig.1 Initial model of UO2 system under thermodynamic equilibrium

4.42 nm×3.91 nm×9.39 nm，包含3 840個U原子和7 680個O原子；圖1c中z軸方向對應UO2材料體系的〈111〉晶向，尺寸為4.06 nm×3.91 nm×11.49 nm，包含4 320個U原子和8 640個O原子。然后在z軸方向上施加拉伸應變，應變比例為1%、2%和3%，UO2泊松比為0.325[8]。最后利用Bai等[9]的方法在體系中間位置添加1個氧間隙原子，在微正則系綜(NVE)下，模擬計算整個體系在6 ns內的微觀結構變化情況，同時利用Voronoi cell[10]方法判斷體系內的間隙原子并跟蹤啞鈴型氧間隙缺陷質心的運動情況。

結合LAMMPS程序計算得到的啞鈴型氧間隙缺陷的質心坐標和愛因斯坦方程[11](式(1))，計算不同溫度下啞鈴型氧間隙缺陷的擴散系數。 計算時，將模擬計算時間內質心坐標隨時間的變化分為K段，每段時間為τK(合理取值范圍為20～250 ps[12])。

(1)

其中：D(T)為T溫度下的擴散系數；Ri(T)為第i段時間內啞鈴型氧間隙缺陷的質心位移；n為啞鈴型氧間隙缺陷的擴散維度。計算過程中，在τK的每種取值情況下，選取不同的10個數據起點用于計算D(T)后取平均值，作為該τK對應的擴散系數。繪制τK與D(T)的關系曲線，擬合D(T)值，作為該溫度下啞鈴型氧間隙缺陷的擴散系數。

2 結果與討論

2.1 間隙原子排列形態和擴散現象

為研究啞鈴型氧間隙缺陷的擴散行為，以z軸方向(對應〈100〉晶軸方向)上3%拉伸應變和無應變兩種情況對比，去除正常的原子僅保留間隙原子，觀察氧間隙原子的排列情況，結果列于表1。可見，在向體系內添加1個氧原子的情況下，形成了穩定的氧間隙原子缺陷，其中包含2個氧原子，且2個氧原子沿〈111〉方向排列，呈啞鈴型，即形成了啞鈴型氧間隙缺陷。表1中藍色為原體系中的氧原子，黃色為添加的1個氧間隙原子，氧原子排布方向構型為添加1個氧間隙原子后進行短暫弛豫，使添加原子后的新體系在預定溫度下重新建立熱力學平衡后觀察到的啞鈴型氧間隙缺陷的形態。此時，較低溫度(1 000 K和1 300 K)時，添加的氧間隙原子(黃色)在添加位置附近與最近的氧原子(藍色)形成穩定的氧間隙缺陷；而較高溫度(1 500 K和1 800 K)時，在重新建立熱力學平衡的過程中，通過碰撞運動，添加的氧間隙原子(黃色)變成正常晶格點位的氧原子，而在新的位置由原本正常點位的氧原子(藍色)形成新的氧間隙原子，并在其附近與最近的氧原子(藍色)形成穩定的氧間隙缺陷。

表1 啞鈴型氧間隙缺陷排列方向Table 1 Configuration of oxygen dumbbell interstitial defect

圖2為1 300 K下沿〈100〉軸向無應變和施加3%拉伸應變時啞鈴型氧間隙缺陷的位置。可看出，無應變時啞鈴型氧間隙缺陷的運動集中在初始位置附近，擴散運動現象不明顯；施加3%拉伸應變時，啞鈴型氧間隙缺陷呈現明顯的三維擴散現象。

2.2 溫度對啞鈴型氧間隙缺陷擴散的影響

為研究溫度對啞鈴型氧間隙缺陷運動的影響，分別將體系弛豫到1 000、1 300、1 500、1 800 K溫度，在無應變和〈100〉方向施加3%拉伸應變條件下模擬計算啞鈴型氧間隙缺陷的擴散現象。利用愛因斯坦方程計算擴散系數，并繪制擴散系數D和1/kBT的關系曲線，結果如圖3所示。由圖3可見，模擬計算得到的擴散系數與1/kBT在對數坐標下呈線性關系，與Arrhenius方程描述的一致；隨著溫度的升高，擴散系數增大，表明在高溫情況下由于原子的熱運動增強，增大了啞鈴型氧間隙缺陷的擴散能力；同時，3%拉伸應變條件下的擴散系數大于無應變情況下的擴散系數，表明拉伸應變能提升啞鈴型氧間隙缺陷的擴散能力。

圖2 1 300 K下沿〈100〉軸向無應變和施加3%拉伸應變時啞鈴型氧間隙缺陷的運動情況對比Fig.2 Comparison of oxygen dumbbell interstitial defect diffusion under 0% strain and 3% strain along 〈100〉 axial with 1 300 K

圖3 〈100〉方向拉伸應變下 擴散系數與溫度的關系Fig.3 Diffusion coefficient as a function of temperature under tensile strain along 〈100〉 direction

Arrhenius方程表達式如式(2)所示：

D(T)=D0e-Em/kBT

(2)

其中：D0為頻率因子，由擴散系數和1/kBT的曲線關系獲得；kB為玻爾茲曼常數；Em為擴散激活能。通過計算可得頻率因子D0和擴散激活能Em，結果列于表2。

由表2可知，沿〈100〉方向施加3%應變較無應變的頻率因子D0更大，這與固體中的擴散理論[13]相符。該理論表明，D0正比于晶格常數的平方，由于對體系施加3%應變，實際原子間距離增大，因此具有更大的頻率因子，對啞鈴型氧間隙缺陷的擴散運動具有促進作用。利用Arrhenius方程計算出〈100〉方向施加3%應變條件下啞鈴型氧間隙缺陷的擴散激活能為0.919 2 eV，無應變條件下為0.876 1 eV。

表2 頻率因子D0和擴散激活能Em擬合結果Table 2 Fitting result of frequency factor and migration energy

2.3 應變方向對啞鈴型氧間隙缺陷擴散的影響

晶體中原子的排布存在各向異性，所以不同晶向的晶體性質存在各向異性[14-15]，間隙原子沿不同晶向的運動同樣存在各向異性[16]。為研究啞鈴型氧間隙缺陷沿不同晶向的運動情況，建立了3個體系模型并充分弛豫到1 300 K溫度，z軸分別對應UO2材料的〈100〉、〈110〉和〈111〉晶軸方向，模擬計算過程中拉伸應變方向均為z軸方向。分別模擬了無應變和1%、2%、3%拉伸應變條件下UO2中啞鈴型氧間隙缺陷的擴散過程，計算了每種情況的擴散系數，結果如圖4所示。

圖4 1 300 K下啞鈴型氧間隙缺陷的擴散系數Fig.4 Diffusion coefficient of oxygen dumbbell interstitial defect with 1 300 K

由圖4可見，對應3種拉伸應變方向，擴散系數均隨應變的增大而增加，但變化率不同。擴散系數的對數與應變呈線性變化關系，因此可按式(3)對結果進行擬合分析。

D=AeBε

(3)

其中：ε為應變；A、B為待定系數，A表示無應變條件下啞鈴型氧間隙缺陷擴散系數，B表示擴散系數隨應變變化的顯著程度。表3為待定系數A、B的擬合結果。

表3 待定系數A、B擬合結果Table 3 Fitting result of coefficient A and B

由圖4和待定系數A可看出，無應變條件下，啞鈴型氧間隙缺陷在〈110〉方向的擴散系數最大，在〈111〉方向最小；而在較大應變條件下(3%拉伸應變)，啞鈴型氧間隙缺陷在〈111〉方向的擴散系數最大，在〈110〉方向最小。

由圖4和待定系數B可看出，拉伸應變條件下，啞鈴型氧間隙缺陷在〈111〉方向擴散系數變化最大，在〈110〉方向最小；隨著〈111〉方向拉伸應變的相對增加，擴散系數變化得更加明顯；隨著〈110〉方向拉伸應變的增加，擴散系數變化較為緩慢。所以，〈111〉方向的拉伸應變對啞鈴型氧間隙缺陷的擴散影響最大，反之〈110〉方向的拉伸應變對擴散影響最小，不同晶向的拉伸應變對啞鈴型氧間隙缺陷的擴散系數影響效果不同，存在各向異性。

2.4 拉伸方向擴散系數的比例

由于啞鈴型氧間隙缺陷在UO2體系內的擴散運動呈現三維現象，為研究啞鈴型氧間隙缺陷沿拉伸方向的擴散運動趨勢，引入參數Dz/D來表示應變方向的擴散系數占總擴散系數的比例，其中Dz表示啞鈴型氧間隙缺陷在z軸(應變方向)的分量，利用z方向位移替換愛因斯坦方程中三維均方位移，即可得到Dz。

圖5 1 300 K下拉伸方向的Dz/DFig.5 Dz/D along tensile strain direction with 1 300 K

圖5為沿〈100〉、〈110〉和〈111〉 3種晶軸方向在無應變和1%、2%、3%拉伸應變條件下UO2中啞鈴型氧間隙缺陷沿應變方向的Dz/D。如果三維擴散運動是各向同性的，那么Dz/D將會趨近于0.33。從圖5可看出，無論有無拉伸應變，〈100〉方向的Dz/D約為0.5～0.6，大于0.33，這表明啞鈴型氧間隙缺陷有明顯的沿〈100〉方向擴散的趨勢；隨著拉伸應變的增加，啞鈴型氧間隙缺陷在應變方向的擴散系數比例并沒有明顯改變。由此可知，雖然增加拉伸應變促進了啞鈴型氧間隙缺陷的擴散運動，表現為擴散系數增加，但并未改變拉伸方向的擴散運動比例，且無論有無拉伸應變，〈100〉方向始終是啞鈴型氧間隙缺陷的主要擴散方向。

3 結論

本文利用分子動力學方法，研究了沿〈100〉、〈110〉和〈111〉 3種晶軸方向分別施加無應變和1%、2%、3%拉伸應變的條件下，啞鈴型氧間隙缺陷在UO2材料體系內的擴散現象。觀察到啞鈴型氧間隙缺陷的擴散運動是三維的，在擴散過程中始終呈現〈111〉方向排列，且無論有無拉伸應變，〈100〉方向始終是啞鈴型氧間隙缺陷的主要擴散方向；研究結果驗證了Arrhenius方程對溫度和擴散系數之間關系的描述，并發現了應變和擴散系數的指數變化關系。研究發現，溫度和應變的增加均會促進啞鈴型氧間隙缺陷的擴散運動，且擴散系數隨應變變化的顯著程度受應變方向的影響，其中〈111〉方向的拉伸應變對啞鈴型氧間隙缺陷擴散運動的影響最大；拉伸應變提升了啞鈴型氧間隙缺陷的擴散能力，但并未改變其沿拉伸方向擴散運動的比例。關于應變條件下氧原子擴散現象的研究，還需在氧間隙原子的數量和應變方式等方面進行改進補充，將在以后的工作中繼續深入研究。