基于MOOSE平台的棒状燃料元件性能分析程序开发与验证

邓超群，向烽瑞，贺亚男，牛钰航，巫英伟，田文喜，秋穗正，苏光辉

(西安交通大学 动力工程多相流国家重点实验室，陕西 西安 710049)

福岛核电站事故后，事故容错燃料(ATF)成为研究热点[1]，涂层Zr包壳[2]、复合SiC包壳[3]和U3Si2芯块[4]等多种ATF设计被提出。传统燃料性能分析程序FRAPCON[5]、FRAPTRAN[6]和TRANSURANUS[7]等难以实现ATF复杂结构(如多层包壳结构)的建模和分析。因此，亟需开发新的燃料性能分析程序对ATF进行性能分析及设计改进。

当前，美国爱达荷国家实验室(INL)开发了多物理场燃料性能分析程序BISON[8]，并针对UO2-Zr燃料开展了大量验证，但由于技术保密等原因，我国难以获取其使用权限。Deng和He等[9-10]分别基于FRAPCON4.0、FRAPTRAN2.0程序添加了多层SiC稳态、瞬态分析模型，实现了多层SiC包壳燃料棒全堆芯稳态及典型棒瞬态分析，但程序采用差分离散算法，计算结果并不精细。

为实现对ATF行为的模拟，西安交通大学核反应堆热工水力研究室(NuTHeL)基于开源的多物理场有限元平台MOOSE[11]开发棒状燃料元件性能分析程序BEEs[4]，可实现二维柱坐标下的稳态棒状燃料性能分析。由于国际上ATF实验数据较少，目前主要针对BEEs程序中添加的UO2-Zr燃料物性与行为模型开展验证工作，为后续基于BEEs程序实现ATF性能分析奠定基础。

1 控制方程

BEEs程序中燃料热传导方程如下：

(1)

变形计算中假定芯块、包壳均处于静态平衡状态：

(2)

其中：σ为柯西应力张量，Pa；f为体积力，N·kg-1。

2 数学物理模型

当前BEEs程序中添加的UO2-Zr燃料物性及行为模型主要来自MATPRO报告[12]、FRAPCON[5]和BISON理论手册[13]。

2.1 燃耗计算

燃耗计算采用TUBRNP模型[14]，局部平均核子密度为：

(3)

2.2 间隙换热

燃料包壳和芯块的间隙换热模型[15]如下：

hgap=hg+hs+hr

(4)

其中：hgap为间隙换热系数，W·m-2·K-1；hg为气体换热系数，W·m-2·K-1；hs为芯块包壳接触换热系数，W·m-2·K-1；hr为芯块外表面、包壳内表面辐射换热系数，W·m-2·K-1。

2.3 包壳氧化腐蚀

燃料包壳氧化层厚度计算采用以下模型[16]：

(5)

其中：i和i+1分别为上一时间步和当前时间步；s为包壳氧化层厚度，m；Δw为氧化物增重，g·cm-2；t为时间，d；A为氧化系数，A=6.3×109μm3·d；R为气体常数；stran为模型转换厚度，stran=2.0×10-6m；Ti为氧化层交界面温度，Ti=To+q″s/λ，K；To为冷却剂温度，K；q″为包壳外表面热通量，W·m-2；λ为ZrO2导热系数，W·m-1·K-1；γ为氧化层增重与厚度的转换因子，γ=0.678 9 cm-3·g；k0=11 863+3.5×104×(1.91×10-15φ)0.24g·cm-2·d，φ为快中子通量密度；Q1和Q2为不同阶段包壳氧化激活能，cal·mol-1。对于M5包壳，Q1=27 446 cal·mol-1，Q2=29 816 cal·mol-1；对于ZIRLO包壳，Q1=27 446 cal·mol-1，Q2=27 354 cal·mol-1；对于Zr-2和Zr-4包壳，Q1=32 289 cal·mol-1，Q2=27 354 cal·mol-1。

2.4 密实化、重定位与辐照肿胀

UO2燃料芯块密实化、重定位和辐照肿胀采用ESCORE[17]和MATPRO[12]经验模型：

(6)

(7)

exp(-0.016 2(2 800-T))·

exp(-0.017 8ρBu)

(8)

(9)

(10)

(11)

其中：T为燃料温度，℃；q′为平均线功率，kW·ft-1。

2.5 芯块和包壳蠕变

UO2芯块蠕变采用MATPRO FCREEP[12]模型，包壳热蠕变和辐照蠕变计算分别采用Power-law模型[18]和Hoppe模型[19]：

(12)

(13)

(14)

2.6 裂变气体释放

BEEs采用Forsberg-Massih模型[20]模拟裂变气体释放过程，裂变气体在晶粒内的扩散方程如下：

(15)

其中：C为晶粒内裂变气体浓度，m-3；Deff为裂变气体有效扩散系数，m2·s-1；β为裂变气体产生率，m-3·s-1。

边界条件为：

(16)

其中：r为距晶粒中心的距离，m；a为晶粒半径，m；b为重溶率，%；λ为重溶半径，m；Nb为晶界气体浓度，m-2。

当晶界气体浓度达到饱和值时就会释放。晶界气体饱和浓度Ns为：

(17)

其中：γ为表面张力，J·m-2；θ为二面角半角，(°)；F为气泡体积转换因子；Vc为晶界临界覆盖率，%；KB为玻尔兹曼常量；pext为气泡外压，Pa。

3 程序框架

MOOSE是基于有限元方法、JFNK算法的开源平台。当前MOOSE平台中已开发了多种物理场和丰富的扩展接口，支持多场耦合和全耦合、全隐性求解计算。

BEEs程序总体框架如图1所示。通过继承MOOSE平台C++类模板，本文基于MOOSE平台的张量力学模块、热传导模块和接触模块，添加了压水堆燃料相关行为模型和非线性燃料元件热膨胀系数、导热系数等物性模型，实现了棒状燃料元件性能分析程序的开发。求解计算时，采用有限元库Libmesh完成网格读取与有限元方程离散，由PETSc进行串行或并行求解。

图1 BEEs程序框架Fig.1 Diagram of BEEs code

4 程序模块验证与校核

在前期开发工作中，通过对比BISON和BAQUS程序，验证了BEEs程序针对UO2-Zr燃料温度、空腔内压、间隙宽度变化的计算功能[4]。本文通过对比FRAPCON程序和实验数据，进一步验证BEEs针对长期稳态工况下燃料燃耗、氧化层厚度和燃料温度模拟的准确性。

4.1 燃耗计算校核

随燃耗的加深，径向功率分布形状由于边缘效应将发生显著变化，而燃料性能参数和裂变产物肿胀等行为模拟都与局部燃耗密切相关。FRAPCON程序[5]的燃耗计算模块经大量验证，计算结果可靠。本文针对同一基准题进行模拟，计算UO2燃料径向功率分布，计算参数列于表1。图2示出径向功率分布因子对比。由图2可见，在燃耗初始阶段及较高燃耗阶段，BEEs程序计算结果与FRAPCON的均符合较好，校核了程序燃耗模块的正确性。

表1 燃耗基准题计算参数Table 1 Calculation parameter of burnup benchmark problem

4.2 氧化层计算验证

Oconee Rod 15309实验[21]在Oconee PWR中完成了5个辐照周期，并对燃耗末期的氧化层厚度进行了测量，实验中燃料棒设计参数列于表2，燃料类型为UO2，包壳类型为Zr-4。图3示出氧化层厚度轴向分布对比。由图3可知，BEEs程序计算结果与FRAPCON程序计算结果和实验值均符合良好，验证了BEEs程序计算氧化层厚度结果的合理性。

a——零燃耗；b——平均燃耗为33 MW·d/kgU图2 径向功率分布因子对比Fig.2 Comparison of radial power factor

表2 Oconee Rod 15309实验设计参数Table 2 Design parameter of Oconee Rod 15309 test

图3 氧化层厚度轴向分布对比Fig.3 Comparison of axial distribution of oxide thickness

4.3 燃料温度计算验证

Halden IFA432实验[22]是于Halden反应堆进行的针对BWR-6燃料棒的长期稳态辐照实验，实验中由下端伸入的热电偶测得燃料棒中心温度，实验中燃料棒Rod1和Rod3的设计参数列于表3，燃料类型为UO2，包壳类型为Zr-2。图4示出Rod1和Rod3的燃料中心温度的变化。由图4可知，BEEs程序计算结果与实验值总体符合良好，与FRAPCON计算结果基本一致，验证了BEEs程序在长期稳态工况下燃料温度变化预测的准确性。

表3 Halden IFA432实验燃料棒Rod1和Rod3的设计参数Table 3 Design parameter of Rod1 and Rod3 for Halden IFA432 test

图4 Rod1(a)和Rod3(b)燃料棒中心温度对比Fig.4 Fuel rod centerline temperature comparison of Rod1 (a) and Rod3 (b)

5 BR3 Rod实验算例整体分析

BR3 Rod实验[23]是于比利时BR3反应堆开展的高燃耗辐照实验。本文采用BR3 Rod实验进行整体分析并与FRAPCON及实验结果进行对比，进一步验证BEEs程序对压水堆棒状燃料行为模拟的准确性。实验中燃料棒的相关参数列于表4，燃料类型为UO2，包壳类型为Zr-4，平均线功率的变化如图5所示。

表4 BR3 Rod实验燃料棒设计参数Table 4 Design parameter of fuel rod for BR3 Rod test

图5 燃料平均线功率随时间的变化Fig.5 Average linear power of fuel vs. time

模拟时将碟形加倒角燃料芯块层叠形状简化单一燃料棒，燃料棒顶部为平端面。网格划分时在燃料径向划分20个节点，轴向划分100个节点，包壳径向划分4个节点，轴向划分80个节点，网格类型为QUAD4。

图6示出包壳内表面、芯块外表面和芯块中心最高温度随时间的变化。结合图5燃料平均线功率变化可知，在279 d时，由于功率激增，燃料棒整体迅速升温，在热膨胀作用下芯块、包壳间隙减小，直至最终间隙闭合，如图7所示。

图6 包壳内表面、芯块外表面和芯块中心最高温度随时间的变化Fig.6 Maximum temperature at cladding inner surface, fuel outer surface and fuel centerline vs. time

图7 间隙尺寸随时间变化Fig.7 Gap width vs. time

选取燃料棒中心高度处的包壳内表面和芯块中心线温度进行对比，结果分别如图8所示。由图8可知，BEEs程序预测结果合理，其计算结果与FRAPCON的基本一致。

裂变气体释放计算结果对比示于图9。辐照前期，产生的少量裂变气体扩散至晶界沉积。在279 d时，由于功率激增，裂变气体生成量显著增加，同时燃料棒温度提升到较高水平，使得裂变气体扩散率较大，积累在晶界的气体浓度超过饱和值，裂变气体开始释放。BEEs程序预测裂变气体释放份额与FRAPCON程序的符合良好，但总体高于实验测量值，计算结果偏保守。

图8 包壳内表面(a)和芯块(b)在中心高度处温度的对比Fig.8 Comparison of temperature of cladding inner surface (a) and fuel centerline (b) at axial mid-plane

图9 裂变气体释放份额对比Fig.9 Comparison of fission gas release portion

6 总结

本文基于开源的多物理场有限元平台MOOSE开发了棒状燃料元件性能分析程序BEEs。通过与FRAPCON程序和相关实验测量值的对比，初步验证了程序的稳态工况模拟能力。结果表明，BEEs能实现燃料燃耗、氧化层厚度和燃料温度的准确计算及裂变气体释放的合理预测，且能针对典型压水堆燃料元件在长期稳态工况下的整体性能进行合理分析。