金属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原子在八面体间隙位之间的扩散与在四面体间隙位之间的扩散难易程度相当。