压水堆燃料辐照性能多物理耦合并行分析程序开发

韩智杰，何晓军,*，明 春，杨衍康，任 帅，胡长军，杨 文

(1.中国原子能科学研究院 反应堆工程技术研究所，北京 102413；2.北京科技大学，北京 100083)

核燃料元件是反应堆系统中的核心部件，是反应堆运行过程中热量及裂变反应的主要来源，故燃料元件的堆内性能对反应堆经济性和安全性有重要影响。正常运行条件下，燃料元件设计必须保证其在高温、高压、强辐照环境下的完整性。利用计算模拟的方法可以用较低的成本预测燃料元件堆内行为，达到评价元件堆内辐照服役性能、指导燃料元件设计的目的[1]。因此，燃料元件性能分析程序对燃料研发[2]、安全评价具有重要意义[3]。

燃料元件的堆内行为涉及中子物理、热工水力、固体力学、化学腐蚀等多种物理现象[4]，典型商用压水堆中装载上万根燃料元件，传统分析评价方法通常选取典型或极限工况燃料元件进行分析，存在一定的局限性，故有必要建立先进多物理耦合并行化燃料性能分析软件进行多棒并行分析，高效率、高精度评价反应堆燃料性能，为提高反应堆安全性及经济性提供有效依据[5]。

本工作建立燃料元件温度、变形、裂变气体释放及内压等主要分析模型，通过热工-力学-内压多物理耦合，基于先进并行计算方法，开发先进多物理燃料性能分析程序Athena。

1 数理模型

1.1 程序介绍

Athena是中国数值反应堆原型系统(CVR1.0)的重要组成程序，CVR1.0目标是开发面向E级超算的反应堆模拟软件，包括中子物理、热工水力、结构力学、燃料性能和材料性能等核心分析软件及多物理耦合环境[6]。为充分发挥超算优势，适应全堆芯pin-by-pin[7]物理-热工-燃料精细化多物理耦合应用需求[8]，Athena利用先进并行算法，建立了燃料并行分析能力，具有高效、精细化全堆芯燃料性能分析优势。

1.2 温度分析模型

燃料元件温度计算在性能分析中占有重要地位，它是所有分析的基础。温度分析模型的主要功能是获得燃料元件的温度分布，即计算热量在冷却剂、包壳、芯块-包壳间隙、芯块内部的传递。

将燃料元件及冷却剂简化为单通道单相传热条件，冷却剂温度通过求解质量、动量和能量守恒方程得到。冷却剂到包壳传热方式为对流换热，包壳表面温度可通过牛顿冷却定律得到[9]。

包壳及燃料芯块可简化为不同热导率材料组成的轴对称圆柱体，故其内部热量传递可用稳态傅里叶导热方程(式(1))描述，差分求解节点划分示意图如图1所示。图1中：δ为径向节点间距，m；N为径向节点数量。

图1 燃料棒温度差分求解节点示意图Fig.1 Schematic of fuel rod node for temperature model

(1)

式中：k为热导率，W/(m·K)；S为体积元表面积，m2；n为表面标准单位向量；q(r)为热源，W/m3；T为温度，K；V为体积元体积，m3；r为空间径向坐标，m。

假设芯块中心线处的径向温度梯度为零，则径向轴对称条件和表面温度边界条件为：

(2)

式中：rfi和rfo分别为芯块中心线处半径和芯块外表面半径；Tfo为芯块外表面温度。

温度分布计算模型的关键之一是材料热导率的计算[10]。芯块及包壳热导率可通过材料物性分析模型得到。芯块-包壳间隙传热系数考虑由气体导热、辐射传热及接触导热3种组成，即：

hgap=hsolid+hr+hgas

(3)

式中：hgap为间隙传热系数；hsolid为接触传热系数；hr为辐射传热系数；hgas为气体传热系数。

1.3 变形计算模型

燃料应力应变计算是包壳失效分析的基础，且会通过芯块-包壳间隙尺寸直接影响温度分布计算。特定温度、内外压力作用下燃料变形物理过程描述方程包括平衡方程、几何方程以及本构方程[11]。

平衡方程为：

(4)

几何方程为：

(5)

本构方程为：

(6)

(7)

(8)

在外力作用下，材料可能发生屈服，塑性变形材料应变的函数可表示为σ=f(ε)，如图2所示。则对于塑性应变εP，屈服应力可通过下式求得：

图2 应力-应变曲线Fig.2 Stress-strain curve

(9)

可写为：

(10)

该非线性方程可通过牛顿迭代法求解得到：

(11)

1.4 裂变气体释放及内压计算模型

伴随燃料在堆内的裂变反应会产生气态和固态裂变产物。裂变产物会引起燃料肿胀，加重芯块与包壳的机械相互作用(PCMI)。气态裂变产物在燃料内溶解度较低，会从燃料芯块释放到燃料棒气空间，增加燃料棒内压，限制燃料棒的寿命；同时降低棒内气体热导率，从而使燃料温度升高。故裂变气体释放行为一直都是燃料元件性能分析的重要内容之一。

裂变气体累积产量与发生的裂变次数相关：

(12)

式中：GPT为该区域裂变气体累计产量；Bu为该区域累积燃耗；VF为该区域燃料体积；Av为阿伏伽德罗常数；PR为裂变气体生成比率，主要考虑的裂变气体为Kr和Xe。

裂变气体通常认为通过扩散进行释放，利用晶内扩散及晶间行为模拟。假设晶体为球体，裂变气体从中心点向外扩散，通过求解扩散方程可得到晶界位置的裂变气体浓度，基本扩散方程[12]为：

(13)

2 程序开发

2.1 程序结构

Athena程序系统设计如图3所示。为开展全堆芯反应堆物理、热工水力、燃料性能耦合分析，设计接口包括输入输出接口、物理耦合接口、热工耦合接口。输入输出接口用于向燃料软件提供初始输入及计算结果输出。物理耦合接口与反应堆物理计算软件耦合，为燃料提供功率分布。热工耦合接口与热工水力计算软件耦合，为温度模块提供边界温度分布。程序耦合接口可与CVR 1.0耦合，更加真实地描述反应堆中燃料棒工况，为全堆芯并行化分析提供基础。

图3 燃料元件性能分析程序系统Fig.3 Fuel performance analysis code system

2.2 多物理耦合方案

Athena程序多物理耦合计算流程如图4所示。通过中子物理模块获得燃料功率分布，经温度模块计算出燃料棒温度分布，利用力学模块得到燃料应力-应变更新燃料包壳间隙尺寸，通过间隙传热系数进行热工-力学耦合；随后计算裂变气体释放及燃料气腔内压，更新内压并反馈到温度和力学模块计算中，当气压达到收敛时完成耦合计算。

图4 程序耦合计算流程图Fig.4 Coupling flow chart of code

2.3 程序并行化

程序并行方式如图5所示。利用程序接口定义详细燃料功率、冷却剂边界条件后燃料棒之间相互独立[13]，Athena程序通过消息传递接口(MPI)技术中的单程序多数据(SPMD)方法实现多棒并行分析，使用主进程将读取的输入文件广播给所有子进程，以便所有子进程都能获取到输入文件中的参数信息，并根据读取的参数信息进行对应燃料棒温度、应力应变、内压、包壳腐蚀等行为的数值计算，进程执行完成后，得到所有棒的性能分析结果，实现反应堆燃料元件性能并行化处理。

图5 燃料元件多棒并行方式Fig.5 Multi fuel rod parallel mode

3 程序计算结果

采用典型商用压水堆核电站燃料数据和同类软件计算结果对程序进行初步验证。对标同类软件为法国原子能委员会(CEA)、法国电力公司(EDF)和法玛通公司联合开发的轻水堆燃料元件性能分析程序Meteor。

3.1 温度场分析

燃料棒平均燃耗和功率随时间的变化如图6所示。燃料棒所在的组件在反应堆内经历两个辐照循环，第1个循环长度为394 d，第2个循环长度为526 d，两个循环燃料组件在堆芯内的位置不同，平均功率略有不同。燃料棒燃耗随辐照时间加长而逐渐加深。燃料芯块峰值温度随时间的变化如图7所示。反应堆启动后，由于芯块“重定位”作用，芯块直径会瞬间增大，间隙减小，间隙热导随之增加，芯块温度降低。随着反应堆的运行，芯块和包壳间隙逐渐减小甚至接触，即芯块-包壳间隙闭合。间隙闭合后，间隙传热系数相对稳定，燃料芯块中心温度变化与功率变化相似，寿期末，随着燃耗加深，燃料芯块热导率降低，虽然功率降低但芯块峰值温度稍有升高。

图6 燃料棒平均燃耗及平均线功率随时间的变化Fig.6 Fuel rod average burnup and average power vs time

图7 芯块峰值温度随时间的变化Fig.7 Fuel pellet peak temperature vs time

3.2 变形计算

寿期末燃料芯块及包壳直径沿轴向分布如图8所示。燃料包壳受到内、外压作用，由于冷却剂外压大于内压，在长期蠕变作用下包壳寿期末直径小于制造直径。由于受到燃料元件两端端塞影响，燃料包壳两端蠕变变形较小，而分析程序未考虑端塞的支撑作用，同时两端芯块肿胀、热膨胀等变形较小，造成包壳两端向内蠕变变形较大。对于燃料元件中间区段，两程序计算结果与测量值符合较好。

图8 燃料元件寿期末直径沿轴向分布Fig.8 Axial distribution of fuel cladding diameter at end of life

燃料元件包壳轴向伸长随时间的变化如图9所示。由图9可看出，在寿期初包壳由于蠕变等作用长度缩短，随燃耗逐渐增加，辐照生长变形明显增加，包壳轴向变长[14]。平均直径、燃料元件包壳伸长数据列于表1。

表1 寿期末燃料元件尺寸结果Table 1 Fuel element dimension result at end of life

图9 包壳轴向伸长随时间的变化Fig.9 Axial elongation of cladding vs time

3.3 裂变气体释放

裂变气体释放随时间的变化如图10所示，裂变气体释放份额随燃耗的增加逐渐增加，当气态裂变产物在晶界内累积且燃料温度高时，裂变气体释放速率明显增加。

图10 裂变气体释放随时间的变化Fig.10 Fission gas release vs time

裂变气体释放数据列于表2。试验单独测量了Kr和Xe两种裂变气体元素释放，其中Xe释放量大于Kr。Athena相较Meteor程序，燃料棒内腔体积及内压计算结果与测量结果吻合更好，总体裂变气体释放数据预测合理。

表2 裂变气体释放数据Table 2 Data of fission gas release

4 结论

本工作开发了数值反应堆原型系统CVR1.0燃料性能分析程序Athena，能够实现燃料元件热工-力学行为高性能模拟，具备多棒并行计算能力。利用典型商用压水堆核电站数据及同类程序计算结果对Athena程序进行了初步验证，对比分析了燃料温度、变形、裂变气体释放及内压等计算结果，结果表明Athena程序计算结果可靠。

开发高精细、强耦合、高效并行运行的软件系统是准确模拟核反应堆燃料性能的基础[15]。CVR1.0旨在充分利用国产超算的优势解决核反应堆复杂的工程问题，最终将结合多尺度、多物理、先进并行技术推动核反应堆分析软件创新发展。Athena程序充分利用并行技术在精细建模中的优势，为CVR1.0全堆芯多物理耦合提供了重要支持。