金屬Zr中O原子擴散的第一性原理研究

趙晨曦 謝耀平 姚美意 胡麗娟

(上海大學材料研究所，上海 200072)

鋯的熱中子吸收截面小，在高溫高壓水中具有良好的耐腐蝕性能和力學性能，因而鋯合金被廣泛用作壓水堆核電站核燃料元件的包殼材料，且在未來20～30年內(nèi)鋯合金仍將是壓水堆核燃料包殼的首選材料[1]。2011年福島核事故的發(fā)生主要是由于核燃料包殼鋯合金在冷卻劑缺失(loss of coolant accident, LOCA)時反應堆堆芯得不到充分的冷卻，堆內(nèi)溫度持續(xù)上升使冷卻劑變成水蒸汽的同時包殼溫度急劇上升，與高溫蒸汽發(fā)生劇烈的氧化反應產(chǎn)生大量熱量和氫氣，導致了反應堆堆芯熔化和氫氣爆炸[2]。冷卻劑缺失事故也即失水事故發(fā)生時，堆芯內(nèi)冷卻條件惡化，積蓄的大量熱量和裂變產(chǎn)物的衰變熱無法導出，甚至可能導致核電站防御體系的屏障(燃料包殼、一回路壓力邊界、安全殼)功能全部喪失。由此可見LOCA是反應堆運行過程中非常嚴重的事故，因此設計反應堆時將LOCA作為設計基礎事故考慮[3]。包殼是確保裂變產(chǎn)物不泄露的第一道屏障，因此有必要深入研究鋯合金的LOCA行為，為研發(fā)性能更優(yōu)異的新型鋯合金提供依據(jù)和指導。

LOCA工況下，鋯合金的殘余塑性是決定包殼安全性的重要因素。LOCA工況下鋯合金包殼先后經(jīng)歷高溫蒸汽氧化的兩個峰值和淬火：反應堆一回路失水，壓力下降，引發(fā)安注系統(tǒng)啟動向堆內(nèi)注水以冷卻迅速升溫的堆芯。但這僅能帶走部分余熱，堆內(nèi)冷卻水仍會很快漏完，燃料元件表面仍產(chǎn)生泡核沸騰和膜態(tài)沸騰，堆芯因冷卻不良致使包殼溫度上升并達到第一個峰值；接著高壓儲水罐和低壓注水泵相繼啟動注水，壓力容器內(nèi)水位開始上升并達到活性區(qū)的下端，這時包殼溫度開始下降，但由于冷卻不良及衰變熱的作用會使包殼溫度進一步升高，并達到第二個峰值，即失水事故中包殼最高溫度(peaking cladding temperature, PCT)；隨后由于冷卻條件的改善使包殼溫度逐漸降低，在LOCA末期包殼因淬火而迅速降溫[4]。高溫蒸汽氧化淬火后鋯合金包殼進一步脆化[5]，嚴重威脅核反應堆的安全，因此高溫蒸汽氧化后的組織需有一定的殘余塑性。

鋯合金包殼經(jīng)過高溫蒸汽氧化后的典型顯微組織有3層結(jié)構(gòu)：外層為氧化物，稱為ZrO2層；中間層為O穩(wěn)定的α相，稱為α- Zr(O)層；內(nèi)層為原α相。α- Zr(O)層的塑性較差，會降低鋯合金的殘余塑性，因此LOCA高溫氧化時鋯合金中α- Zr(O)層的厚度是評價鋯合金淬火后殘余塑性的重要指標[6]。由于在LOCA工況下鋯合金發(fā)生α→β轉(zhuǎn)變的同時O在鋯合金中大量固溶，隨著鋯合金中O含量的上升，α→β轉(zhuǎn)變溫度升高[7]，因此阻礙鋯合金發(fā)生α→β轉(zhuǎn)變而形成O穩(wěn)定的α- Zr(O)層。在給定的氧化溫度和時間條件下，鋯合金中α- Zr(O)層的厚度僅取決于ZrO2/α- Zr界面的氧濃度[8]。鋯合金在高溫蒸汽中氧化形成的各層顯微組織中O的分布和含量與O的擴散行為密切相關[9]，所以O在鋯合金中的擴散是影響α- Zr(O)層厚度的主要因素。為了深入了解O在鋯合金中的擴散行為，本文采用第一性原理計算方法，研究了O原子在金屬Zr中的熱力學和動力學性質(zhì)，以更好地理解鋯合金的LOCA行為。

1 計算方法和模型

1.1 計算方法

本文采用基于密度泛函理論(density functional theory, DFT)[10- 11]的第一性原理方法，通過VASP(Vienna ab- initio simulation package)[12- 14]軟件進行理論計算，交換關聯(lián)泛函采用廣義梯度近似(general gradient approximation, GGA)[15]下的Perdew- Burke- Ernzerhof泛函[16]。波函數(shù)采用平面波基組展開，計算中所有截斷能均為400 eV。為了研究O原子在金屬Zr中的熱力學性質(zhì)，需計算O原子在超胞中不同位置的形成能，計算公式為：

Ef int=E(ZrnO)-E(Zrn)-E(O)

(1)

Ef sut=E(Zrn-1O)-E(Zrn)/n×

(n-1)-E(O)

(2)

式中：n為完美超胞中Zr原子的數(shù)目，Ef int和Ef sub分別為O原子位于間隙位置和置換位置時體系的形成能；E(ZrnO)和E(Zrn-1O)分別對應O原子位于超胞的間隙位置和置換位置時的體系總能；E(Zrn)和E(O)分別為完美超胞的體系總能和一個O原子的能量。

1.2 模型構(gòu)建

首先，計算密排六方結(jié)構(gòu)(hcp)α- Zr的晶格常數(shù)為：a=3.227 ?,c=5.175 ?，與試驗結(jié)果(a=3.23 ?,c=5.15 ?)[17]及前人的理論計算結(jié)果(a=3.230 ?,c=5.176 ?)[18]均較吻合，說明本文的計算是可靠的。接著，構(gòu)建如圖1(a)所示的54個原子hcp超胞，在hcp- Zr54超胞中計算O原子位于不同的高對稱位置時的形成能、空位(Vacancy, Va)與O的相互作用以及O原子的擴散勢壘。在計算O原子的擴散勢壘時，主要采用微動彈性帶(nudged elastic band, NEB)[19]方法。計算時所設幀數(shù)值為16，將初態(tài)到末態(tài)的變化分為多個過程，從而計算出從初態(tài)到末態(tài)的具體變化路徑和能量變化曲線，最后得到從初態(tài)擴散到末態(tài)所需克服的最大能量。本文在計算體系總能時，布里淵區(qū)K點取樣網(wǎng)格均為2×2×1，計算擴散能壘時，K點取樣網(wǎng)格為1×1×1。

2 結(jié)果與討論

2.1 hcp- Zr中O的穩(wěn)定性研究

hcp- Zr中有兩個典型的高對稱性間隙位置，即八面體間隙位置(octahedral interstice position，oct)和四面體間隙位置(tetrahedral interstice position，tet)，如圖1(b,c)所示。結(jié)合如圖1(d)所示的置換位置(substitutional position，sub)，本文研究了O在這3種占位時晶體的基本性質(zhì)。采用式(1)和式(2)分別計算得到O在上述3種占位時的基本性質(zhì)，并將結(jié)果列于表1。

圖1 hcp- Zr54超胞及O原子在hcp- Zr54超胞中的3種位置

表1 O原子在hcp- Zr54超胞中3種位置(oct、tet和sub)對應的計算結(jié)果

由表1可知，當O處于間隙位置oct和tet時，O形成能為負值，表明O在進入Zr間隙位置時為放熱過程，因此O在Zr間隙位置時具有熱穩(wěn)定性。而且oct位置比tet位置更穩(wěn)定。當O處于sub位置時，O形成能為正值，表明O進入Zr后取代Zr格點的原子為吸熱過程，因此O在Zr中的替位式占據(jù)形式并不穩(wěn)定。晶格常數(shù)的計算發(fā)現(xiàn)，當O進入Zr后會引起晶格畸變，而且O在不同的位置引起的晶格畸變程度不同。對于純Zr，c和a分別為15.38和9.78 ?，軸比c/a為1.58。當O在oct位置時c/a增大，而另兩種占位均導致c/a減小。

2.2 hcp- Zr超胞中Va對O穩(wěn)定性的影響

如圖2所示，通過去掉超胞中的一個Zr原子，構(gòu)建一個原子Va，并在其周圍標出近鄰的oct位置。利用該模型，通過式(1)可以計算原子Va與其第一到三近鄰間隙(the first to the third near neighbor，1 to 3 nn)位置O原子形成缺陷對時的體系總能，計算結(jié)果如表2所示。從表2可見，間隙O與Va之間的形成能均為負值，表明間隙O與Va結(jié)合過程中會釋放熱量，O與Va形成的復合缺陷是穩(wěn)定的，因此O在Zr中易與Va結(jié)合。從表2還可以看出，當O處于Va的第三近鄰時，其結(jié)構(gòu)在三者中最穩(wěn)定，說明O和Va雖然相互吸引，但若靠得太近，距離小于第三近鄰時，將有輕微的排斥效應。綜上所述，Va的存在使O原子在Zr中的穩(wěn)定性提高。相比純Zr，hcp- Zr54超胞中oct、tet和sub 3種位置加入O原子后，軸比分別增加了1.3%、1.4%和1.2%。

表2 Va- Va近鄰的O (oct)在hcp- Zr54超胞中的計算結(jié)果

圖2 hcp- Zr54超胞中的Va- Va近鄰oct位置

2.3 hcp- Zr中O的擴散性質(zhì)

通過NEB分別計算了O在Zr中oct位置之間和tet位置之間的兩種擴散，其擴散路徑如圖3所示，能量變化如圖4所示。體系狀態(tài)包括：初態(tài)(initial state, IS)、過渡態(tài)(transition state, TS)、中間態(tài)(middle state, MS)和末態(tài)(final state, FS)。由圖4可知，O在Zr中兩個近鄰oct位置之間擴散時需要連續(xù)越過兩個勢壘，兩個勢壘中間有一個亞穩(wěn)態(tài)。O在兩個近鄰tet間隙位置的擴散也有類似的特點。擴散勢壘方面，oct- oct與tet- tet兩條路徑的擴散勢壘分別為1.78和0.89 eV。因此，從oct出發(fā)的O原子需要越過更大的勢壘進行擴散，而從tet出發(fā)的原子只需克服較小的勢能。但是，由于原子擴散難易程度也受絕對勢壘的影響，oct- oct與tet- tet擴散路徑的絕對勢壘分別為-468.40和-468.26 eV，oct- oct擴散路徑的絕對勢壘相對較低。如果以oct位置為參考點，oct- oct和tet- tet的相對勢壘分別為1.78和1.92 eV。因此O原子在兩條路徑擴散需要的勢壘相近，導致擴散難易程度相當。

圖3 O原子的擴散結(jié)構(gòu)圖

圖4 O原子擴散的能量變化曲線(體系狀態(tài)包括初態(tài)、過渡態(tài)1、中間態(tài)、過渡態(tài)2和末態(tài))

3 結(jié)論

采用密度泛函理論研究了O原子在hcp- Zr超胞中的穩(wěn)定占位以及擴散狀況。從計算結(jié)果可知：O在八面體間隙位置時形成能最低，表明O原子在hcp- Zr超胞中穩(wěn)定的占位為八面體間隙位置；在hcp- Zr超胞空位與間隙O原子的形成能為負值，說明空位提高了O原子的穩(wěn)定性；利用NEB方法計算了hcp- Zr中O擴散的路徑，O原子擴散能壘的計算結(jié)果表明O原子在八面體間隙位之間的擴散與在四面體間隙位之間的擴散難易程度相當。