張華陽，王慶宇，李忠宇

(哈爾濱工程大學(xué) 核安全與仿真技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

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

1 方法

本文的模擬計(jì)算均利用分子動(dòng)力學(xué)軟件LAMMPS[5]進(jìn)行，并結(jié)合軟件OVITO[6]對模擬計(jì)算結(jié)果進(jìn)行分析。LAMMPS是大規(guī)模原子分子并行模擬器，能模擬計(jì)算百萬級的原子或分子系統(tǒng)的微觀結(jié)構(gòu)演化，OVITO是一款可視化分析軟件，通過讀取LAMMPS的計(jì)算結(jié)果能顯示每個(gè)原子或分子的運(yùn)動(dòng)情況，并能對材料的微觀結(jié)構(gòu)進(jìn)行相關(guān)分析。對于UO2材料體系，本文利用Yakub等[7]開發(fā)的分?jǐn)?shù)離子模型勢函數(shù)描述體系中各原子之間的相互作用力，此勢函數(shù)能精確描述原子間的相互作用以及每個(gè)原子的運(yùn)動(dòng)情況。

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

圖1 達(dá)到熱力學(xué)平衡時(shí)UO2材料體系初始模型Fig.1 Initial model of UO2 system under thermodynamic equilibrium

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

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

(1)

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

2 結(jié)果與討論

2.1 間隙原子排列形態(tài)和擴(kuò)散現(xiàn)象

為研究啞鈴型氧間隙缺陷的擴(kuò)散行為，以z軸方向(對應(yīng)〈100〉晶軸方向)上3%拉伸應(yīng)變和無應(yīng)變兩種情況對比，去除正常的原子僅保留間隙原子，觀察氧間隙原子的排列情況，結(jié)果列于表1?？梢姡谙蝮w系內(nèi)添加1個(gè)氧原子的情況下，形成了穩(wěn)定的氧間隙原子缺陷，其中包含2個(gè)氧原子，且2個(gè)氧原子沿〈111〉方向排列，呈啞鈴型，即形成了啞鈴型氧間隙缺陷。表1中藍(lán)色為原體系中的氧原子，黃色為添加的1個(gè)氧間隙原子，氧原子排布方向構(gòu)型為添加1個(gè)氧間隙原子后進(jìn)行短暫弛豫，使添加原子后的新體系在預(yù)定溫度下重新建立熱力學(xué)平衡后觀察到的啞鈴型氧間隙缺陷的形態(tài)。此時(shí)，較低溫度(1 000 K和1 300 K)時(shí)，添加的氧間隙原子(黃色)在添加位置附近與最近的氧原子(藍(lán)色)形成穩(wěn)定的氧間隙缺陷；而較高溫度(1 500 K和1 800 K)時(shí)，在重新建立熱力學(xué)平衡的過程中，通過碰撞運(yùn)動(dòng)，添加的氧間隙原子(黃色)變成正常晶格點(diǎn)位的氧原子，而在新的位置由原本正常點(diǎn)位的氧原子(藍(lán)色)形成新的氧間隙原子，并在其附近與最近的氧原子(藍(lán)色)形成穩(wěn)定的氧間隙缺陷。

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

圖2為1 300 K下沿〈100〉軸向無應(yīng)變和施加3%拉伸應(yīng)變時(shí)啞鈴型氧間隙缺陷的位置?？煽闯?，無應(yīng)變時(shí)啞鈴型氧間隙缺陷的運(yùn)動(dòng)集中在初始位置附近，擴(kuò)散運(yùn)動(dòng)現(xiàn)象不明顯；施加3%拉伸應(yīng)變時(shí)，啞鈴型氧間隙缺陷呈現(xiàn)明顯的三維擴(kuò)散現(xiàn)象。

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

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

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

圖3 〈100〉方向拉伸應(yīng)變下 擴(kuò)散系數(shù)與溫度的關(guān)系Fig.3 Diffusion coefficient as a function of temperature under tensile strain along 〈100〉 direction

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

D(T)=D0e-Em/kBT

(2)

其中：D0為頻率因子，由擴(kuò)散系數(shù)和1/kBT的曲線關(guān)系獲得；kB為玻爾茲曼常數(shù)；Em為擴(kuò)散激活能。通過計(jì)算可得頻率因子D0和擴(kuò)散激活能Em，結(jié)果列于表2。

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

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

2.3 應(yīng)變方向?qū)♀徯脱蹰g隙缺陷擴(kuò)散的影響

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

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

由圖4可見，對應(yīng)3種拉伸應(yīng)變方向，擴(kuò)散系數(shù)均隨應(yīng)變的增大而增加，但變化率不同。擴(kuò)散系數(shù)的對數(shù)與應(yīng)變呈線性變化關(guān)系，因此可按式(3)對結(jié)果進(jìn)行擬合分析。

D=AeBε

(3)

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

表3 待定系數(shù)A、B擬合結(jié)果Table 3 Fitting result of coefficient A and B

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

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

2.4 拉伸方向擴(kuò)散系數(shù)的比例

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

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

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

3 結(jié)論

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