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

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

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

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

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

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

1 計算方法和模型

1.1 計算方法

本文采用基于密度泛函理論(density functional theory, DFT)[10- 11]的第一性原理方法，通過VASP(Vienna ab- initio simulation package)[12- 14]軟件進行理論計算，交換關聯泛函采用廣義梯度近似(general gradient approximation, GGA)[15]下的Perdew- Burke- Ernzerhof泛函[16]。波函數采用平面波基組展開，計算中所有截斷能均為400 eV。為了研究O原子在金屬Zr中的熱力學性質，需計算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原子的數目，Ef int和Ef sub分別為O原子位于間隙位置和置換位置時體系的形成能；E(ZrnO)和E(Zrn-1O)分別對應O原子位于超胞的間隙位置和置換位置時的體系總能；E(Zrn)和E(O)分別為完美超胞的體系總能和一個O原子的能量。

1.2 模型構建

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

2 結果與討論

2.1 hcp- Zr中O的穩定性研究

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

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

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

由表1可知，當O處于間隙位置oct和tet時，O形成能為負值，表明O在進入Zr間隙位置時為放熱過程，因此O在Zr間隙位置時具有熱穩定性。而且oct位置比tet位置更穩定。當O處于sub位置時，O形成能為正值，表明O進入Zr后取代Zr格點的原子為吸熱過程，因此O在Zr中的替位式占據形式并不穩定。晶格常數的計算發現，當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穩定性的影響

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

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

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

2.3 hcp- Zr中O的擴散性質

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

圖3 O原子的擴散結構圖

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

3 結論

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