确定论数值反应堆程序NECP-X的开发及应用

刘宙宇，许晓北，温兴坚,2，王 博,2，曹 璐,3，周欣宇，张旻婉,3，陈俊辑,2，王习宁，李 帆，李帅铮，马升泽，邵世豪，曹良志，吴宏春

(1.西安交通大学 核科学与技术学院，陕西 西安 710049；2.中国核动力研究设计院，四川 成都 610213；3.西北核技术研究所，陕西 西安 710024)

核能技术自20世纪40年代问世以来，因其清洁性、稳定性以及能量密度高等优点得到大规模的发展和应用。但传统反应堆模拟技术由于一系列的理论近似，导致其计算分辨率低、计算精度不足和适用范围小等问题，严重制约了现有核反应堆进一步提高经济性和安全性，以及新型核反应堆系统的设计。近年来，计算机技术和数值方法的快速发展以及物理模型的改进推动了数值反应堆技术的发展。数值反应堆基于大规模并行计算平台，利用先进的物理模型和数值模拟算法，采用精细化建模，从而可精确模拟反应堆在正常运行与事故工况中发生的各类物理现象。相比于传统模拟技术，数值反应堆技术具有精度高、分辨率高等优点，对于提高核能系统设计的经济性、安全性具有重要意义。目前数值反应堆研究已在全世界得到广泛开展，如美国和韩国联合开展的I-NERI计划[1-2]、美国的CASL计划[3]、NEAMS计划[4]、欧盟的NURESIM计划[5]等。

本文简要介绍针对数值反应堆的中子学计算方法研究，包括共振计算方法、输运计算方法、输运-燃耗耦合方法等及相应的数值反应堆物理程序，以及在此基础上进行的物理-热工-燃料性能分析的多物理耦合研究，并将其应用于商用大型压水堆、研究堆等各类堆型中，以实现反应堆的高精度、精细几何建模、多物理耦合的大规模并行数值模拟，从而准确预测反应堆在运行过程中的关键安全参数随时间的演变，为反应堆的设计及安全分析提供指导。

1 理论模型

1.1 连续能量和多群堆用核数据库

对基于ENDF/B-Ⅶ.0评价核数据库制作连续能量截面数据库和多群截面数据库的制作方法进行研究，可为共振计算和输运计算提供截面信息。WIMS格式69群能群结构原有的共振能量范围为9.118 keV～4 eV，即15～27群，由于无法考虑Pu在低能量区的共振，对MOX燃料的计算存在较大偏差。因此将共振能量段的下限拓展到0.625 eV，将共振能量段的上限拓展到24.78 keV。拓展之后的共振能群为15～45群。WLUP发布的多群核数据库只有28个核素当作共振核素处理，因此在WIMS库的基础上将共振核素的数目从28增加到66以上。WLUP发布的WIMS格式多群核数据库中的平均裂变谱是根据235U、238U和239Pu的裂变谱按照0.54、0.08和0.38的比例加权得到的，但该裂变谱不是对所有问题都适用，因为各核素的裂变谱和平均裂变谱有所区别，因此在数据库中给出了平均裂变谱及可裂变核素和易裂变核素的裂变谱。连续能量数据库覆盖的温度范围为293～1 800 K，温度点间隔为50 K。多群数据库温度点同样为293～1 800 K，温度点间隔为300 K。

1.2 基于全局-局部耦合的共振计算方法

提出了全局-局部耦合方法[8]进行共振计算。其思想是：将共振中的有关效应分为全局效应与局部效应。全局效应为全局的空间自屏效应、全局的温度分布效应和和全局的共振干涉效应；局部效应为局部的空间自屏效应、共振弹性散射效应、局部温度分布效应、局部共振干涉效应、多群等效效应和边缘效应。采用丹可夫修正因子处理全局效应；采用共振伪核素子群方法处理局部效应。

对于全局效应，采用中子流方法计算丹可夫修正因子Cb：

(1)

式中：φf,0为孤立棒问题的中子标通量密度；φf,1为真实问题的中子标通量密度。

对于局部效应，首先需要计算慢化剂区的外半径。设燃料至慢化剂的首次飞行碰撞概率为Pf→M，其为慢化剂外半径的函数，因此丹可夫因子可写成：

(2)

式中：rM为慢化剂区的外半径；Pe,f为中子逃脱概率；Pf→M为燃料至慢化剂的首次飞行碰撞概率。对于给定的丹可夫修正因子，rM可通过二分法得到。

局部效应的计算基于共振伪核素子群方法，其步骤为：

(3)

2)采用最小二乘拟合获得共振伪核素的子群参数，并获得各个组成的共振核素子群截面；

3)求解子群固定源问题，并获得各共振核素的有效自屏截面；

4)采用SPH修正以保证输运计算的反应率守恒。

在此基础上研究并提出了可处理复杂几何的等效几何共振计算方法。该方法的思想为：首先，针对复杂几何燃料的孤立问题，基于燃料的逃脱概率守恒，建立复杂几何燃料模型的等效一维圆柱(或平板)燃料模型；其次，基于燃料到外围结构材料区的碰撞概率守恒，获得燃料外围结构材料的等效尺寸；然后，根据复杂几何燃料的丹可夫因子守恒，建立等效一维圆柱(或平板)燃料外围的慢化剂尺寸；最后，针对等效一维圆柱(或平板)燃料模型，采用伪核素子群方法进行有效自屏截面计算。

1.3 基于2维/1维耦合的输运计算方法

开展了高效稳定的2D/1D耦合输运计算方法[9]研究。角度、能群离散后的三维中子输运方程为：

Qg(x,y,z)

(4)

式中：ε、η、μ分别为x、y、z方向的方向余弦；Σt,g(r)为总截面；ψg,m(x,y,z)为空间(x,y,z)上m方向的角通量，即中子输运方程中的未知量；Qg(x,y,z)为总源项，包括裂变源及散射源。

首先，对式(4)在第L层轴向高度ΔzL区间内进行积分得到径向二维方程：

Σt,g(x,y)ψg,m(x,y)=

(5)

(6)

之后，对原始三维问题在棒p上沿径向在ΔxiΔyj区间内进行积分，得到轴向一维方程：

(7)

ψg,m,i,j-1/2(z)]

(8)

针对2D/1D耦合方法中由于泄漏项导致的负源问题，研究了改进的泄漏项分割技术[11]，通过燃料棒内角度的各项异性重构获得每个平源区的角通量，即：

(9)

为进一步提高计算效率，开发了基于区域分解松耦合的3级粗网有限差分(CMFD)方法[10,12]实现输运计算的加速。该方法将CMFD线性系统在空间和能量上解耦成若干个子线性系统，在并行计算中，每个计算核心独立求解各自的子线性系统，求解完成后再进行通信。该方法可显著提升CMFD并行计算效率[12]。但由于该方法将空间及能群解耦，造成迭代格式退化，导致CMFD外迭代次数增加。为进一步提升CMFD的计算效率，NECP-X采用三重CMFD方法[12]。在该方法中，通过单群组件CMFD计算加速单群栅元CMFD计算，再通过单群栅元CMFD计算加速多群栅元CMFD计算，最终通过多群栅元CMFD计算加速输运计算。对于单群栅元CMFD计算，其系数可基于多群栅元CMFD计算结果进行能群归并，得到：

(10)

对于单群组件CMFD计算，其系数可基于单群栅元CMFD计算结果进行空间归并，得到：

(11)

1.4 基于反应率预测的预估校正方法

在燃耗计算中，通过求解每个核素的Bateman方程得到反应堆中各核素的原子核密度：

(12)

式(12)的解可用矩阵形式表示：

N(t)=etAN(0)

(13)

式中：N(t)为t时刻的原子核密度矩阵；A为燃耗矩阵；N(0)为初始时刻的原子核密度矩阵。式(13)可通过切比雪夫近似方法[13]进行高效求解。

对于燃耗方程与输运方程的耦合计算，常采用预估-校正方法。为提升计算效率，NECP-X中研究了基于反应率预测的预估校正方法[14]。通过在前几个燃耗点预测预估步和校正步核反应率的变化规律，减少后续燃耗点校正步的物理-热工耦合计算，从而有效提高物理-热工-燃耗耦合的计算效率。

1.5 基于Picard迭代的物理-热工-燃料性能分析耦合计算方法

核反应堆是一个多物理耦合系统，涉及中子通量、核子密度、温度、密度、应力、应变等物理场，这些物理场相互影响。要准确预测反应堆的状态，必须实现多物理耦合计算分析。NECP-X通过Picard方法实现多个物理场的耦合。其思想为：将相互耦合的多物理场分解为多个独立的物理场，再逐次求解各物理场，在求解某个物理场过程中，其他物理场采用初始值或上一迭代步计算得到的值，各物理场的耦合物理变量相互迭代，直至收敛。如对于反应堆模拟中普遍的核热耦合问题，将耦合问题分解为中子学、热工水力两个子问题。首先求解中子学问题，其所需温度采用上一时间点求解得到的温度，可得到堆芯中子通量及对应的功率；接着基于此功率进行热工水力求解，得到的温度再反馈给中子学问题进行再次求解，依此往复，直至通量、温度场均收敛。Picard方法的优点是可充分利用各物理场现有的成熟计算程序，易实现。

在核反应堆堆芯模拟中，反应堆物理-热工耦合(也称核热耦合)已得到广泛应用[15]。但对燃料元件在堆芯服役时的弹塑性形变、蠕变、裂变气体释放等燃料行为的分析是离线分析，其并未考虑各种燃料行为(如间隙传热、热膨胀、肿胀、密实化、蠕变等)对物理和热工计算的影响，结果并不准确。为更精确地进行核反应堆内多物理耦合现象的模拟，应对其进行物理-热工-燃料性能的紧密耦合计算，因此，数值反应堆中必须开展多物理耦合计算研究[16-18]。基于NECP-X程序，研究了物理-热工-燃料性能耦合计算，其计算流程如图1所示。在物理-热工-燃料性能全耦合计算中，每个物理场在每个燃耗点进行Picard迭代，收敛后进入下一个时间步。耦合计算中所传递的物理量如图2所示。

图1 NECP-X物理-热工-燃料性能全耦合计算的流程

图2 NECP-X物理-热工-燃料性能全耦合中传递的物理量

2 程序开发与验证

基于上述理论开发了数值反应堆物理计算程序NECP-X。该程序基于面向对象思想，采用Fortran2003语言；基于Git与CMake进行版本管理与跨平台编译，并针对每个独立模块建立单元测试以保证代码质量。程序具有稳态计算、瞬态计算、燃耗计算、源项计算、多物理耦合计算等功能，具备三维全堆芯精细几何建模能力与万核以上大规模并行计算能力。目前，程序已经过大量的国际知名基准问题(如VERA基准题、BEAVRS基准题等)验证[14,19]，具备很高的计算精度与计算效率。

3 应用

3.1 商用大型压水堆模拟

1)二代改进型压水堆M310

我国多个电厂采用的是二代改进型M310核电技术，该反应堆一回路拥有3个环路，堆芯共157个燃料组件，使用的燃料组件呈17×17布置，并包含多种类型的燃料组件。堆芯从中心到边缘可分成3个区，不同富集度的燃料组件呈棋盘状布置，同时为补偿后备反应性，部分组件使用不同数量的Pyrex可燃毒物棒。采用NECP-X针对M310堆芯进行高保真精细几何建模，其堆芯径向布置如图3所示。

图3 几何模型堆芯径向布置

采用NECP-X模拟计算某M310堆芯第一循环，其临界硼浓度的计算偏差示于图4，其中FPC为传统预估校正方法，IFPC为基于反应率预测的预估校正方法。采用IFPC方法时，程序于第4～6燃耗时间点获得多群中子能谱偏差随燃耗时间的二阶规律，并从第7燃耗时间点才开始多群中子能谱偏差的插值计算以预测更新反应率，使IFPC方法真正发挥作用，因此与FPC的不同。由图4可知，总体上NECP-X计算值与测量值符合良好；基于反应率预测的预估校正方法的临界硼浓度计算精度更高。采用NECP-X程序计算的第一循环精细棒功率分布示于图5。由图5可更直观地看出核反应堆棒功率分布的展平现象。此外还采用NECP-X模拟了第一循环不同燃耗时间点的燃料有效温度(在耦合问题中共振自屏计算所采用的燃料栅元代表温度，通过对热工计算得到的燃料栅元内部温度分布进行处理得到，本文采用的有效温度计算公式为：有效温度=0.3×燃料中心温度+0.7×燃料表面温度)分布及冷却剂温度分布，结果示于图6。由图6可见，堆芯燃料有效温度分布随燃耗的进行逐渐变得平坦；第一循环燃料有效温度分布在576.0～1 027.5 K之间，且燃料有效温度最大值出现在循环初，并随燃耗的进行不断降低，同时燃料有效温度最大值所在轴向位置不断上移。第一循环不同燃耗时间点三维非均匀冷却剂温度分布示于图7。由图7可见，冷却剂温度分布变化与燃料温度变化趋势保持一致；其最大值小于613 K，且与功率水平呈正相关，其中第一循环末期功率下降时，冷却剂出口温度下降显著。

图4 第一循环临界硼浓度NECP-X模拟计算偏差

图5 NECP-X模拟计算的第一循环精细棒功率分布

图6 NECP-X模拟计算的第一循环有效燃料温度分布

图7 NECP-X模拟计算的第一循环慢化剂温度分布

2)第三代压水堆AP1000

AP1000压水堆是由美国西屋公司研发的非能动安全压水堆。其首循环堆芯是一循环期为18个月的低泄漏设计，使用从天然铀至富集度为4.38%的5个燃料区。AP1000采用了两类可燃毒物棒，即整体式可燃毒物棒IFBA及湿环形离散式可燃毒物棒WABA，对于离散式WABA，根据毒物长度的不同分为长、中和短WABA。AP1000堆芯的布置如图8所示，为增加经济性，实现18个月低泄漏循环，同时考虑到展平功率分布和反应性控制，径向AP1000采用5区装料方案。按照燃料富集度由低至高分为A、B、C、D和E区，其富集度依次为0.74%、1.58%、3.20%、3.78%和4.38%。

图8 AP1000首循环四分之一堆芯布置

基于上述几何模型，采用NECP-X进行AP1000启动物理试验计算，并与蒙特卡罗程序KENO计算得到的参考解进行对比。其中，计算参数包括硼微分价值及各类温度系数，对于临界硼浓度，NECP-X与蒙特卡罗程序的偏差为5 ppm，与测量值的偏差为24 ppm；硼微分价值(DBW)、等温温度系数(ITC)、燃料温度系数(DTC)和慢化剂温度系数(MTC)列于表1。由表1可看出，硼微分价值及各类温度系数的NECP-X计算结果与蒙特卡罗程序计算结果符合得很好。将热态零功率有效增殖因数和组件归一化功率的NECP-X计算结果与蒙特卡罗程序JMCT的计算结果进行对比，JMCT基于径迹长度统计功率。蒙特卡罗程序JMCT和NECP-X程序计算的热态零功率工况下的有效增殖因数keff分别为1.000 44和0.999 85，二者的偏差为59 pcm。两程序的热态零功率组件归一化功率计算结果示于图9。从图9可见，除靠近围板处实际功率较小导致相对偏差为2%外，其余位置的偏差均小于1.5%，可见，NECP-X有很高的计算精度。

表1 启动物理试验中硼微分价值及各类温度系数计算结果

图9 热态零功率组件归一化功率计算结果

3.2 研究堆模拟

1)脉冲反应堆

西安脉冲堆是我国自行设计建造的国内第一座实用性脉冲反应堆，是在中国核动力研究设计院设计建造的首座原型脉冲堆的基础上，根据用户对该堆的应用要求进行设计、建造的。西安脉冲堆工程属小型研究堆，堆芯以正六边形环形布置，共9圈，采用粗棒型燃料元件，燃料棒直径为3.6 cm，栅距为4.3 cm，两者约为一般压水堆的4倍，燃料采用UZrH1.6芯体，235U富集度为19.75%，不锈钢作包壳材料，芯体中插有锆合金棒，燃料元件101根，稳态控制棒5根，材料为B4C，下部设置有燃料跟随体，脉冲控制棒1根(堆芯右端控制棒)，材料也为B4C，下部是一密闭空气腔，石墨元件86根，充当反射层，不锈钢吸收体元件2根，堆芯有稳态和脉冲两种堆芯布置方式，本文针对其稳态工况进行模拟，其径向布置如图10所示。

图10 脉冲堆径向布置

基于以上模型，使用NECP-X和蒙特卡罗程序MCX分别针对脉冲堆稳态工况进行模拟，其中MCX基于径迹长度针对功率进行统计。MCX和NECP-X计算得到的堆芯有效增殖因数分别为1.015 94±0.000 01和1.015 89，二者的偏差仅为5 pcm，符合得很好。功率分布计算结果示于图11，其中蒙特卡罗程序计算的功率相对统计偏差小于0.01%。由图11可见，NECP-X计算得到的功率偏差在-1.64%～5.2%之间，值得注意的是，功率偏差较大的5.2%和4.98%两根棒是处于全插位置的控制棒，其归一化功率参考解分别为0.155和0.157，远小于正常燃料棒，故功率相对偏差较大。控制棒全提和全插工况下的堆芯有效增殖因数的计算结果列于表2。从表2可见，蒙特卡罗程序和NECP-X在控制棒全插工况下的有效增殖因数的偏差为78 pcm；控制棒全提工况下的偏差为84 pcm。两程序的结果符合得很好。

图11 脉冲堆稳态工况径向归一化功率计算结果

表2 脉冲堆控制棒全提与全插工况有效增殖因数计算结果

2)板状燃料反应堆

JRR-3M(日本第3号改进型研究反应堆)是日本建造的一座高性能多用途研究反应堆，于1990年3月首次达到临界状态，并于当年11月开始运行，最大输出功率为20 MW。JRR-3M是一种池式研究反应堆，使用低浓缩铝化物板状燃料，采用板型组件和方形堆芯排列，第一循环燃料组件装载图如图12所示。堆芯内布置标准燃料组件、跟随燃料组件、控制棒组件和固定式控制棒组件和辐照孔组件，共37盒，控制棒组件中的控制棒有6根(调节棒R1和R2、补偿棒S1和S2、安全棒SA1和SA2)，除这6根控制棒外，第一循环还在3个辐照孔位置处布置有固定式控制棒(N、T1和T2)。

图12 JRR-3M第一循环径向布置

活性区高度为75 cm，上下反射层各30 cm。

采用NECP-X程序对JRR-3M二维全堆芯控制棒全提和全插两种工况进行模拟，计算结果与蒙特卡罗程序MCX计算得到的参考解进行对比，其中MCX基于径迹长度进行功率统计。有效增殖因数的计算结果列于表3，在控制棒全提工况下两程序的偏差为94 pcm，控制棒全插工况下两程序的偏差为97 pcm。可见NECP-X具有较高的计算精度。两种工况下的功率分布及与蒙特卡罗程序计算结果的偏差示于图13，图中每个组件方格的第1行数据为NECP-X计算结果，第2行数据为偏差。采用蒙特卡罗程序计算得到的功率相对统计偏差均在0.01%以内。从图13可见，对于控制棒全提工况，所有计算功率偏差均小于1.2%；对于控制棒全插工况，最大功率偏差为1.8%，大部分计算功率偏差小于1%。表明NECP-X的计算结果与蒙特卡罗程序符合得很好。

表3 JRR-3M有效增殖因数计算结果

图13 JRR-3M在控制棒全提和全插工况下的功率分布及偏差

3.3 多物理耦合模拟

Watts Bar反应堆一号机组(WBN1)位于美国田纳西州，修建于20世纪80—90年代，并于1996年开始组网发电。该反应堆由美国西屋公司设计，其一回路有4个环路，堆芯使用的燃料组件呈17×17布置，并包含多种类型的燃料组件。本文针对Watts Bar反应堆进行物理-热工-燃料性能混合耦合模拟，即中子学程序与热工水力程序首先进行核热耦合计算，然后基于核热耦合的求解结果搜索得到每个组件内的最热棒，并针对这些棒进行燃料性能分析计算。该反应堆第一循环的堆芯装载方案及控制棒所在的组件示于图14，WBN1第一循环的总体参数列于表4。堆芯额定热功率为3 411 MW，堆芯的额定冷却剂质量流量为16.59 t/s，堆芯运行压力为15.5 MPa。

图14 WBN1第一循环堆芯布置

表4 WBN1第一循环参数

计算得到的燃料温度分布、包壳米塞斯应力分布、燃料与包壳间隙宽度分布分别示于图15～17，需要注意的是图中的一格表示堆芯中每一个组件的最热棒的一层。从图15～17可看见，径向和轴向有明显的由功率展平带来的温度展平现象。整个燃耗期间燃料最大温度为1 709 K；包壳米赛斯应力随时间与空间的变化很小，最大米塞斯应力为118.5 MPa；对于燃料与包壳间隙，可看到所有燃料棒的间隙宽度轴向分布均呈两端大、中间小的趋势。这是因为燃料棒中央段功率密度高，燃料热膨胀和肿胀更明显，造成间隙宽度减小。从燃耗初期到燃耗中期，间隙宽度整体增大，这是因为燃料棒的密实化效应导致燃料芯块收缩，进而导致间隙宽度上升；之后，从燃耗中期到燃耗末期，间隙宽度略减小，这是由燃料芯块的肿胀和包壳的蠕变效应导致的。

图15 WBN1第一循环燃料温度分布

图16 WBN1第一循环包壳米塞斯应力分布

图17 WBN1第一循环燃料与包壳间隙宽度分布

4 总结

本文介绍了针对数值反应堆的堆芯模拟计算方法，包括共振计算方法、输运计算方法、输运-燃耗耦合计算方法、多物理耦合计算方法等，以及基于上述计算方法开发的数值反应堆程序NECP-X，并针对程序进行了验证与应用。NECP-X程序采用了基于ENDF/B-Ⅶ.0评价核数据库产生的连续能量堆用核数据库及多群核数据库，研究了全局-局部耦合的共振计算方法；研究了改进泄漏项分割的2D/1D耦合输运计算方法和区域分解松耦合并行的三重CMFD加速方法；基于切比雪夫有理近似的燃耗计算方法，并通过反应率预测的预估-校正方法提高燃耗计算效率；实现了基于Picard迭代的物理-热工-燃料性能全耦合方法。程序采用Fortran2003语言开发，具有稳态计算、瞬态计算、燃耗计算、源项计算、多物理耦合计算等功能，具备三维全堆芯精细几何建模能力与万核以上大规模并行计算能力。

本文将程序应用于第二代反应堆改进型M310、第三代反应堆AP1000、西安脉冲堆、JRR-3M研究堆、Watts Bar反应堆等堆型的模拟计算中，数值结果表明，程序计算结果合理可信，可真实地预测反应堆在服役过程中各类重要物理量(如有效增殖因数、功率、温度、应力、间隙宽度等)随时间的演变情况，可为商用大型压水堆、研究堆等各类堆型的设计及安全分析提供可靠的工具。