铅基快堆氧化物燃料性能分析程序FUTURE开发与初步验证

王 冠，顾 龙,3，于 锐，王 挺，王 兆，袁 和，恽 迪,*

(1.中国科学院 近代物理研究所，甘肃 兰州 730000；2.中国科学院大学，北京 100049；3.兰州大学，甘肃 兰州 730000；4.西安交通大学，陕西 西安 710049)

21世纪以来，第4代核能系统国际论坛(GIF)面向未来能源需求，从安全性、经济性、可持续性等多个角度出发，选择了6种最具前景的先进反应堆系统[1]。目前，以铅冷快堆、钠冷快堆等为代表的第4代核能技术正逐步走向工程化应用的轨道，这种趋势对快堆的理论认识、计算模拟和安全设计提出了很高的要求。相比于金属燃料，U-Pu混合氧化物(MOX)燃料的稳定性较好，在高燃耗时肿胀更小，并具有较高的熔点，因此MOX燃料是当前快中子增殖堆的首选燃料[2-3]。

在反应堆内部高温、强辐照环境下，燃料内部发生着从微观原子尺度到宏观结构尺度的演化，同时也经历着微秒量级的中子学反馈和数年的服役寿期[4]。相比于热中子堆，快堆内中子能谱更硬、功率密度和燃耗深度更大，造成了MOX燃料元件内部存在显著的燃料重构、元素迁移、裂变气体释放、芯包相互作用以及各种化学效应等[5-6]。这些现象之间强烈的非线性特征和复杂的耦合关系使得燃料元件性能分析只能通过计算机数值模拟实现。因此从20世纪开始，许多国家相继开发了一系列针对液态金属快堆氧化物燃料元件性能分析程序从而实现对燃料元件性能和行为的数值模拟和定量分析，例如欧洲超铀元素研究所(TUI)开发的TRANSURANUS[7]，美国麻省理工学院(MIT)开发的FEAST-OXIDE[8]，法国原子能与可替代能源委员会(CEA)开发的GERMINAL[9]以及日本原子能委员会开发的CEDAR[10]等。

为实现对铅基快堆以及CiADS[11-12]次临界堆氧化物燃料元件性能和行为的模拟，中国科学院近代物理研究所反应堆室和西安交通大学NUSOL团队联合开发了铅基快堆氧化物燃料棒性能分析程序FUTURE。该程序基于串行的半隐式耦合求解方法，可实现对铅基快堆氧化物燃料的热工-力学-辐照-化学等多物理场耦合分析计算。相比于国际上的1.5维快堆燃料性能分析程序，一方面FUTURE程序采用了两步分析法耦合多物理场计算，在实现较高精度求解的同时提升了计算效率，具有后发优势；另一方面，国际上的燃料性能分析程序多针对钠冷快堆进行开发和研究，能够实现铅基快堆燃料性能分析的程序较少，而FUTURE程序初期即面向铅基快堆进行开发，并添加了目前较为先进的液态铅铋(LBE)腐蚀模型。

本文针对FUTURE程序中实现的快谱条件下氧化物燃料特异性行为进行模型的初步验证，为后续基于程序实现铅基快堆氧化物燃料棒性能分析和模拟计算奠定基础。

1 程序介绍

1.1 计算建模

燃料元件是反应堆内具有独立结构的燃料使用单元，也是燃料组件的核心构件。典型的液态金属快堆燃料元件呈棒状结构，外部为细长的金属包壳管(长约为2～3 m，直径约为5～10 mm，厚度约为0.4～0.6 mm[6])，内部包容有燃料活性区、气腔及其他具有定位、测量、屏蔽等功能的部件。FUTURE程序将燃料元件的燃料区和气腔作为主要模拟计算区域，采用1.5维轴向堆叠假设和有限体积法进行简化和区域离散，如图1所示。离散后的控制体为相互嵌套的节点环，从而生成计算网格。为了对快堆MOX燃料棒的复杂行为进行模拟计算，程序基于Fortran95语言进行模块化编程，各计算模块实现相应的计算和分析功能并进行数据传递，如图2所示。

图1 FUTURE程序燃料计算模型Fig.1 Fuel calculation model of FUTURE code

图2 FUTURE程序模块框架Fig.2 Framework of FUTURE code

1.2 计算流程

在快堆高温、大温度梯度和强中子辐照作用下，MOX燃料内部存在显著的孔隙迁移、燃料重构、氧和钚的径向重分布、辐照肿胀、蠕变、裂变气体释放和芯包相互作用等现象。这些现象以及它们之间的相互作用十分复杂，为了对燃料在服役期间的各种行为进行模拟和定量计算，需要通过燃料元件性能分析程序对燃料、包壳和冷却剂的重要特性随辐照时间的演变加以追踪，从而预测燃料棒在各种运行工况下的性能[13-14]。基于串行的半隐式计算架构，程序采用燃料棒全域热力分析和局部行为模拟耦合的两步分析法，能够较好地反映和处理燃料棒在堆内服役过程中的多物理场耦合和多尺度演化。FUTURE程序计算流程如图3所示。

图3 FUTURE程序计算流程Fig.3 Flow chart of FUTURE code

2 数学物理模型

2.1 传热分析

燃料棒服役期间的许多行为都强烈地依赖燃料棒内部的温度分布并且能够反作用于传热过程，为精确计算燃料棒内部的温度分布和演化过程，程序通过耦合力学和燃料行为模型，采用高斯-赛德尔迭代方法对燃料棒传热方程进行数值求解。

芯块和包壳内部的热量传递过程由一维固体热传导方程描述：

(1)

其中：ρ为材料密度，kg/m3；c为材料比热容，J/(kg·K)；T为温度，K；t为时间，s；r为径向位置，m；λ为热导率，W/(m·K)；Q为燃料释热率，W/m3。

选择围绕单根燃料棒的独立冷却剂通道作为流体研究区域，将冷却剂流道类比成内壁面加热、外壁面绝热的环管，采用单通道分析计算冷却剂轴向传热过程，冷却剂轴向能量方程为：

(2)

其中：Al为燃料棒冷却剂流通面积，m2；w为冷却剂轴向流速，m/s；z为轴向位置，m；qf为包壳-冷却剂界面热通量，W/m。

2.2 力学分析

在力学模拟计算中，主要考虑了材料在外部载荷及高温辐照环境下发生的弹性、热膨胀、肿胀、重定位、密实化、蠕变和塑性变形等。为了降低程序力学计算的复杂度，计算引入如下基本假设[15]：1) 连续性假设；2) 轴对称假设；3) 广义平面应变假设；4) 各向同性假设；5) 小变形假设；6) 准静态假设。

力学计算从基本方程出发，采用位移法建立应力、应变与位移的对应关系，并通过全主元高斯-若当消去法进行数值求解，力学控制方程如下。

1) 径向平衡方程

该方程由极坐标下微元体的受力平衡得到，不考虑径向体积力时方程为：

(3)

其中，σr、σθ为径向和环向主应力，Pa。

2) 几何方程

基于小变形假设和广义平面应变假设，柱坐标系下燃料棒模型适用的几何方程为：

(4)

其中：εr、εθ、εz为径向、环向和轴向的主应变；u、w为径向、轴向位移，m；C为常数。

3) 物理本构方程

采用广义胡克定律描述应力-应变的本构关系，在弹性项的基础上添加了热膨胀、肿胀、密实化、重定位以及塑性和蠕变应变，其形式为：

(5)

2.3 燃料重构

在高温和大温度梯度作用下，快堆内运行的氧化物燃料由于初始孔隙的径向迁移和原始晶粒的不断长大，其几何形态和微观结构都会发生很大变化。程序通过求解孔隙率守恒方程[5]计算和模拟燃料内部孔隙的迁移和中心开孔现象，并耦合传热和力学分析，复现快堆柱状芯块向环形芯块的演化过程。孔隙率守恒方程为：

(6)

其中：P为燃料孔隙率；vp为孔隙迁移速率，m/s。

程序采用Sens[16]的孔隙迁移速率公式：

vp=5×10-16×(0.988+6.395×10-6T+

3.543×10-9T2+3×10-12T3)·

(7)

其中：ΔH为蒸发焓，ergs/mol；p0为蒸汽压，dynes/cm2；R为理想气体常数。

2.4 元素重分布

MOX燃料中铀和钚存在许多价态，所以这些元素偏离准确的化学计量是允许的，这种偏离会对燃料芯块的性能产生很大影响。程序基于固体热扩散的迁移动力学理论，计算了MOX燃料中氧和钚的重分布过程，其中用于间隙氧原子和氧空位的迁移动力学方程[17]为：

(8)

其中：cO为间隙氧原子/氧空位的原子分数；DO为扩散系数，m2/s；QO为有效摩尔输运热，J/mol。

用于钚热扩散的迁移动力学方程[18]为：

(9)

其中：cPu为MOX燃料分子中钚元素的原子分数；DU-Pu为铀钚互扩散系数，m2/s；QU-Pu为有效摩尔输运热，J/mol。

2.5 LBE腐蚀

为评估LBE对燃料元件包壳的影响，FUTURE程序基于Zhang[19-20]开发的长周期氧化-腐蚀动力学模型，计算和预测不锈钢包壳在服役期间的腐蚀和氧化情况。

在外部Fe3O4磁铁层未完全去除之前，尖晶石层并不与LBE接触，所以此时尖晶石层的生长符合抛物线定律：

(10)

其中：δSp为尖晶石层厚度，m；kp为氧化常数，m2/s。

外部Fe3O4磁铁层的生长动力学方程为：

(11)

其中：δFe3O4为磁铁层厚度，m；ρFe3O4、ρSt、ρSp、ρLBE为Fe3O4磁铁层、不锈钢基体、尖晶石层、LBE冷却剂密度，kg/m3；FFe,Fe3O4、FFe，St、FFe，Sp为Fe3O4磁铁层、不锈钢基体、尖晶石层中Fe的质量分数；Rm为Fe的质量传输率，m/s。

当LBE流动造成的Fe去除速率小于通过尖晶石层Fe的扩散速率时，LBE对尖晶石层无腐蚀作用，此时尖晶石层遵循线性生长律：

(12)

其中，t0为外层Fe3O4磁铁层去除时间，s。

随着尖晶石层按照线性律不断生长，其厚度最终会趋于一个临界厚度，在此厚度下LBE流动造成的Fe去除速率等于通过尖晶石层Fe的扩散速率，此时LBE对尖晶石层开始存在腐蚀作用，因此尖晶石的生长开始受到Tedmon模型的控制：

(13)

其中，RSp为尖晶石层腐蚀速率，m/s。

3 初步验证

铅铋快堆氧化物燃料元件性能分析程序FUTURE包含有传热、力学、燃料重构、元素重分布、LBE腐蚀等主要计算模块。由于各模块内部计算过程复杂，模块间耦合接口众多，为了测试各计算模块的稳定性和计算精度，保证程序全耦合条件下的稳定运行和正确计算，分别通过与文献算例、基准公式、现有程序进行对比，开展了程序各主要计算模块的初步验证。

3.1 传热验证

传热模块的验证采用与FEAST-OXIDE温度验证相同的基准题[21]，该基准题针对UO2+He2+HT9燃料元件传热材料体系，基于简单的工况和物性假设，模拟计算特定几何结构的燃料元件内部温度分布。该基准题的温度分布满足下列解析表达式：

(14)

其中：T(r)为温度分布，K；q‴为体功率密度，W/m3；ro为燃料芯块外径，m；k为燃料热导率，W/(m·K)；Ts为燃料芯块外表面温度，K。

FUTURE程序传热分析对比验证如图4所示。图4结果表明，在快堆定常工况条件下，燃料内部温度呈现抛物线分布，当线功率密度较高时，燃料中心温度逼近熔点。基于温度基准题线功率66 kW/m和70 kW/m两个工况的计算结果分析，程序的传热计算比较准确，能较好地模拟燃料元件内部的传热过程和温度分布。

图4 FUTURE程序传热分析对比验证Fig.4 Comparison of FUTURE code heat analysis

3.2 力学验证

3.2.1芯块非均匀热弹性验证 通过FUTURE程序和商用的多物理场耦合平台COMSOL对预设温度场的柱状芯块和环形芯块非均匀热弹性算例进行对比，验证程序基于广义平面应变假设和1.5维轴向堆叠模型力学核心算法的可靠度和准确性。力学验证算例的输入参数列于表1，材料物性热膨胀系数α、杨氏模量E、泊松比ν和密度ρ如下式所列。

α=1.0×10-5K-1

(15)

E=2.334×1011×(1-1.091 5×

10-4T)×(1.0-2.752×

(1-0.925))×(1+0.15×0.2)

(16)

ν=(1-0.2)×(0.316+(T-300.0)/

2 800.0×(0.5-0.316))+0.2×

(0.276+(T-300.0)/2 800.0×(0.5-0.276))

(17)

(18)

表1 力学验证算例输入参数Table 1 Input parameter of mechanical analysis case

柱状芯块和环形芯块非均匀热弹性算例预设温度分布如图5所示，力学计算对比如图6所示。由图6可见：在高温和大温度梯度条件下，FUTURE程序力学计算与COMSOL的计算结果符合较好，其中反映燃料变形的位移计算基本无偏差；von Mises应力计算趋势相同，结果相近，FUTURE程序柱状芯块和环形芯块应力计算与COMSOL计算结果的最大相对偏差分别为8.27%和6.48%。基于上述燃料非均匀热弹性算例的计算结果分析，FUTURE程序针对不同结构的静力学计算具有较高的准确度和可靠性，能实现快堆MOX燃料几何结构从柱状向环形演化过程的力学分析和模拟。

a——柱状芯块；b——环形芯块

a——柱状芯块；b——环形芯块

3.2.2芯包机械相互作用验证 在快堆高温、高燃耗条件下，燃料元件内部存在着比较严重的芯包机械相互作用(PCMI)。通过FUTURE程序和COMSOL对环形芯块芯包接触工况算例进行对照，验证程序执行PCMI计算的准确度和可靠性。环形芯块PCMI工况算例温度分布如图7所示，力学计算对比如图8所示。

图7 环形芯块PCMI工况温度分布Fig.7 Given temperature distribution of annular pellet PCMI case

由图8可见，在接触工况条件下环形燃料的最大应力出现在芯块的内壁面和外壁面处，而接触压作为芯块和包壳外部载荷影响了内部的应力分布和结构变形。对于芯块而言，环向应力和轴向应力在内壁面呈压应力，外壁面处呈拉应力状态。这种应力状态会导致芯块外缘出程序计算的芯包内部各应力分量和位移与COMSOL结果符合较好，其中径向应力计算基本无偏差，环向应力和轴向应力最大相对偏差为8.19%和11.11%。基于接触工况的对比验证，表明FUTURE程序具备处理和模拟芯包机械相互作用的能力。

图8 环形芯块PCMI工况力学计算对比Fig.8 Comparison of annular pellet PCMI analysis

图9、10为FUTURE程序加载典型快堆功率谱并开启热力耦合(热弹性)功能计算得到的UO2燃料中心线温度和接触压与COMSOL模拟的结果对比。由图9可见，FUTURE程序计算得到的燃料中心线温度随功率的演化与COMSOL计算结果趋势完全相同，最大相对偏差为3.9%。由图10可见，FUTURE程序计算的接触压与COMSOL的相比变化趋势相同，数值相近，但仍存在一定偏差，最大相对偏差在20%左右。但FUTURE程序计算的接触后芯包边界的贯穿深度比COMSOL计算的更好，其贯穿深度更小。从方法和应用上来看，COMSOL是基于有限元法开发的三维多物理场耦合平台，在力学计算上具有更高的精度和灵活性，而FUTURE程序采用了简化的有限差分方法处理力学问题，受方法和网格量限制，精度无法比及COMSOL，但在现有框架下已能较好地模拟快堆燃料棒的热力耦合问题。

图9 快堆工况下UO2燃料中心线温度演化对比Fig. 9 Comparison of UO2 fuel central temperature for fast reactor

图10 快堆工况下UO2燃料FCMI演化对比Fig.10 Comparison of UO2 fuel FCMI simulation for fast reactor

3.3 燃料重构验证

FUTURE程序燃料重构计算选用Sens[16]开发的氧化物燃料孔隙迁移速率公式，并与文献[16]模拟结果进行对比验证，结果如图11所示。

由图11可见，快堆氧化物燃料内部初始分布均匀的孔隙会发生沿径向的定向迁移，这种迁移过程不仅使孔隙率重新分布，影响燃料局部的热力学特性，也会使燃料中心在短时间内富集较多的孔隙导致中心孔的出现。图11a、b结果表明，程序计算的孔隙率迁移与文献结果符合较好，能较好地反映孔隙率的演化过程。图11c显示在孔隙率演化逐渐达到稳态的情况下，FUTURE程序计算的由于孔隙率迁移导致的燃料靠近中心处密实化更大，分析其原因可能为：1) 由于文献[16]中计算模拟的输入条件不明确，程序的孔隙迁移计算采用弱对照的方法，拟合并匹配输入条件，可能会导致计算出现一定的偏差；2) 程序孔隙迁移采用有限体积法进行计算，守恒性好，而文献[16]中采用的数值方法和离散格式未知，可能导致计算结果精度有所差异。

a——30 min；b——1 h；c——2 h

3.4 元素重分布验证

3.4.1氧的重分布验证 FUTURE程序中氧的重分布计算选用Sari等[17]开发的氧重布迁移动力学模型和计算关系式，并与文献[17]模拟结果进行对比，用以验证程序氧分布计算的准确度和可靠性。

快堆氧化物燃料中发生的氧的重分布现象影响了O/M比(氧元素与重金属元素的含量之比)的变化，而O/M比的变化将会直接改变燃料的各种特性并引发一系列现象。快堆氧化物燃料O/M比演化计算对比如图12所示。图12中，R为外半径。由图12可见：在快堆工况下，对于超化学计量的UO2或MOX燃料，氧的定向迁移会使燃料热端的含氧量增加，同时冷端含氧量相应减少；在欠化学计量燃料中，氧的重布则以反向发生，使冷端浓集氧，热端失去氧。通过FUTURE程序计算的不同初始O/M比条件下氧的重分布结果与文献[17]结果对比发现，程序中超化学计量氧的重分布计算与文献[17]符合较好，但在欠化学计量燃料(O/M<2.0)的计算中，推测文献中可能采用固定O/M=2.0为外边界条件，而FUTURE程序采用第二类边界条件，导致外侧区域O/M比的计算有所偏离。

红色线为FUTURE程序计算结果，黑色线为文献[17]结果

3.4.2钚的重分布验证 为模拟和分析快堆MOX燃料中钚元素的重分布现象，FUTURE程序中钚的重分布计算选用Bober等[18]开发的基于固体热扩散机理的动力学模型和经验公式，并与TRANSURANUS程序的模拟结果[22]进行对比，结果如图13所示。

图13 快堆MOX燃料钚元素迁移计算对比Fig.13 Comparison of MOX fuel Pu migration calculation for fast reactor

在快堆MOX燃料运行过程中，由于高温和大温度梯度的作用会在燃料内部发生铀、钚成分的离析。图13表明，在快堆启堆的几十天内，燃料中心处钚的浓度大幅攀升，这种钚的重分布过程会造成局部燃料成分的变化，不仅影响燃料热物理性能，也会由于钚在中心处的富集导致燃料中心温度大幅上升。基于Bober等[18]的热扩散动力学模型，FUTURE程序模拟计算结果与TRANSURANUS程序的符合较好，表明能较好处理和模拟快堆MOX燃料钚的重分布过程。

3.5 LBE腐蚀验证

图14为 FUTURE程序计算的高流速(v=5 m/s)条件下长周期LBE腐蚀情况与文献[20]模拟结果的对比。由图14可见，程序腐蚀计算判定准确，能够较好地切换不同计算模型，从而模拟和追踪长周期HT-9和T-91不锈钢包壳双层氧化层的生长与腐蚀情况。由于文献[20]中部分输入条件不明确，计算采用的数值方法未知，导致程序模拟计算结果在较长周期的腐蚀演化计算中不断出现累积偏差，10 a腐蚀总量的相对偏差为8.82%(HT-9)和9.09%(T-91)。考虑到快堆燃料元件服役寿期远小于10 a，且程序计算具有一定的保守性，因此认为程序LBE腐蚀计算能对不锈钢包壳的腐蚀情况进行较好地模拟和评估。

图14 LBE腐蚀计算对比Fig.14 Comparison of LBE corrosion calculation

4 结论

本文基于串行的半隐式耦合求解方法开发了铅基快堆氧化物燃料性能分析程序FUTURE。FUTURE程序采用两步分析法，能够实现铅基快堆氧化物燃料棒全域热力分析与局部行为模型的多物理场耦合计算。通过与文献算例、基准公式和现有程序的对照，初步验证了程序对快堆氧化物燃料各分离效应的模拟方法，以及对各计算模块求解的准确度、稳定性和误差来源进行了分析。结果表明，FUTURE程序能模拟快堆稳态工况条件下燃料元件内部的温度演化、结构变形、应力分布和相互作用，并实现燃耗模拟、燃料重构、组分迁移、裂变气体释放和LBE腐蚀等计算内容。其中热力分析、燃料重构、氧和钚元素的迁移等计算具有较高的计算精度，燃料服役期内的LBE腐蚀计算也具有较好的计算结果。除此以外，FUTURE程序在实现快堆燃料热力耦合计算中能够保证一定的计算精度和稳定性，为接下来针对程序的全耦合验证和工程应用提供可靠的模拟计算内容，同时也为快堆稳态工况条件下氧化物燃料的行为演化和性能分析奠定基础。