吴宏春，曹良志，郑友琦，李云召，祖铁军，刘宙宇

(西安交通大学 核科学与技术学院，陕西 西安 710049)

近年来，随着核反应堆用途的不断拓展，人们对核能安全性和经济性的要求不断提高，核反应堆的设计越来越趋于复杂化和精细化，对传统的核反应堆设计技术提出了诸多新的挑战，特别是对核反应堆物理计算方法和数值工具提出了更高的要求。在上述需求的驱动下，国际核强国在过去20年的时间中对其物理计算方法和工具进行了深入、广泛的改进，甚至重新研发了新一代的反应堆物理计算程序系统以应对新堆型研发和设计要求提高对物理计算所提出的挑战，具体包括以下两方面。

1) 物理计算方法和程序的通用化。目前反应堆物理程序开发的进展表明，物理计算程序具有更宽的能谱适用范围，满足不同反应堆的物理计算需求已是大势所趋。国际知名的核设计程序如美国西屋公司开发的APA、法国CEA开发的APOLLO-3等均已发展为能同时用于热堆/快堆的物理计算。

2) 物理计算建模的精细化、精确化。随着数值反应堆概念的提出并被广泛接受，基于高性能计算条件的“一步法”反应堆计算方法成为另一重要的研究方向，为实现对核反应堆堆芯的高分辨率建模和高精度计算，产生了一批国际知名的“高保真”物理计算程序，如美国阿贡国家实验室的PROTEUS及密歇根大学的MPACT等。

国内在核反应堆设计方法及相应的计算程序开发上一直受到国外的技术封锁，为打破封锁、实现核电软件的完全自主化，我国在最近10年间开展了大量针对国内在建、在研压水堆的自主化物理程序研发，推出了TORCH、PCM、cosNU等自主化程序[1-2]。然而，上述程序在通用化、精细化或精确化等方面仍有待进一步加强。

西安交通大学反应堆物理团队在核数据及共振理论、中子输运理论及数值方法等方面开展了大量的基础研究，发展了子群方法[3]、超细群方法[4]，提出了小波展开方法[5]、小波-子群耦合方法[6]、全局-局部耦合方法[7]等共振计算方法；发展了模块化特征线方法[8]、2D/1D特征线方法[9]，提出了非结构几何SN节块方法[10]等中子输运数值方法，并最终形成了一整套具有强通用性和高精度的物理计算程序系统，称之为NECP软件包。NECP软件包主体包含NECP-Atlas、NECP-Bamboo、NECP-X和NECP-SARAX(简称为Atlas、Bamboo、X和SARAX)等自主化程序。Atlas程序采用改进的截面处理方法制作高精度的堆用多群核数据；Bamboo程序基于组件均匀化或栅元均匀化，采用改进的子群方法、特征线方法和变分节块方法进行压水堆燃料管理计算；X程序基于全堆一步计算，采用全局-局部耦合共振计算方法、2D/1D特征线方法进行轻水堆高保真计算；SARAX程序基于快谱均匀化思想，采用超细群方法、输运节块方法进行快谱反应堆计算。NECP软件包可实现对各类压水堆和快堆的高精度核设计和燃料管理计算，并具备了进行高保真物理计算的能力。

本文从NECP软件包的核心算法模型出发，对该软件包内主体程序的特点和验证与确认情况进行介绍，并对软件包在实际工程计算中的应用情况进行代表性的阐述。

1 核数据处理程序Atlas

1.1 程序简介

Atlas程序[11]是为满足物理计算程序的堆用核数据需求而开发的核数据处理程序，可将评价核数据库中的数据加工为Bamboo、X、SARAX所需的核截面。该程序可处理ENDF/B-Ⅷ.0、ENDF/B-Ⅶ.1、ENDF/B-Ⅶ.0、ENDF/B-Ⅵ.8、CENDL-3.1、JEFF-3.2及JENDL-4.0等评价核数据，输出WIMS-D、MATXS和ACE格式数据库。

1.2 理论介绍

Atlas程序可将评价核数据库中0 K下的数据转化到所需温度下的连续截面或多群截面，通过对传统核数据处理程序中使用的近似模型进行改进优化，提高了堆用核截面数据的精度。

对于热谱反应堆，由于重核素的热运动及其弹性散射截面的共振，使得中子散射后有很大的概率获得能量发生上散射。传统核数据处理程序中并未对共振弹性散射进行专门的处理，而是假设超热能区靶核处于静止状态并且其弹性散射截面为常数，这样处理后中子发生散射后总是损失能量。针对该问题，Atlas程序采用不同阶数的散射矩进行多普勒展宽，精确考虑靶核热运动及散射截面的共振对散射中子能量分布的影响。在超热能区可认为靶核不受化学键束缚，将靶核的运动按自由气体模型处理，可得到不同温度T下各阶散射矩的计算式：

(1)

对于快谱或中间能谱反应堆，不可分辨共振区核数据的处理对后续计算精度的影响显著。Atlas程序使用Multi-Level Breit-Wigner公式。该公式可避免采用Single-Level Breit-Wigner公式时出现的负弹性散射截面现象，提高了结果的鲁棒性，且可实现对竞争反应概率表的计算，方便堆物理计算程序实现阈能反应的自屏计算。此外，对于积分统计方法，Atlas程序对波动积分中的孤立项、序列内干涉项和序列间干涉项进行了精确的计算，提高了有效自屏截面的精度。

1.3 程序验证

对Atlas程序进行了广泛的基准题验证，包括国际原子能机构的WLUP基准题、OECD的ICSBEP基准题、IRPhEP基准题及VERA基准题等。反应堆的能谱涉及快谱、中间能谱和热谱。图1示出采用Atlas程序制作的数据库和NJOY2016制作的数据库对ICSBEP基准题的计算结果，计算程序为MCNP。由图1可见，MCNP使用两个数据库的计算结果偏差小于25 pcm。采用BEAVRS基准题验证了Atlas程序中共振弹性散射效应的处理精度，表1列出3种计算结果的对比，包括采用Atlas程序制作的未考虑共振弹性散射的数据库(方法1)、考虑共振弹性散射的数据库(方法2)和直接采用多普勒舍弃修正方法制作的数据库(方法3)，计算程序均采用OpenMC程序。表1中，偏差1表示方法1、2间的偏差(方法2的值-方法1的值)，偏差2表示方法2、3间的偏差(方法3的值-方法2的值)。结果表明，对共振弹性散射的准确处理将使压水堆keff计算误差降低100～200 pcm。对不可分辨共振区截面处理的改进相比于传统方法可得到更准确的多普勒常数，计算结果列于表2。

图1 ICSBEP基准题计算结果Fig.1 Calculation result of ICSBEP benchmark

表1 BEAVRS基准题 keff计算结果Table 1 keff calculation result of BEAVRS benchmark

注：括号内为蒙特卡罗程序OpenMC的统计偏差

表2 钠冷快堆燃料组件无限增殖因子kinf和多普勒常数KD的结果Table 2 kinf and KD results of subassembly in sodium-cooled fast reactor

2 压水堆堆芯燃料管理两步法计算程序Bamboo

2.1 程序简介

Bamboo程序[12-13]由二维栅格计算程序Bamboo-Lattice、三维稳态分析程序Bamboo-Core和三维瞬态分析程序Bamboo-Transient组成。程序可用于压水堆核设计和燃料管理计算，具有临界和目标轴向功率偏移搜索、启动物理试验模拟、基负荷运行模拟、功率历史跟踪模拟及标签换料计算等主要功能。目前，Bamboo程序已开发了两个版本，其中，Bamboo1.0是基于组件均匀化的计算程序，Bamboo2.0是基于更加先进的栅元均匀化的pin-by-pin计算程序。

2.2 理论模型

Bamboo程序基于由ENDF/B-Ⅶ评价核数据库经由Atlas程序加工制作的69群多群截面数据库、核素燃耗数据库和动力学参数库，包含了258种核素，可直接分析硼、银-铟-镉、钨、铪、钆等多种中子吸收体材料。采用全组件一步子群共振计算方法，通过式(2)给出有效自屏截面：

(2)

组件输运计算采用组件模块化特征线方法，显式考虑涂硼可燃毒物(IFBA)和水隙等重要精细结构的影响。采用显式中子俘获释热方法，如式(3)所示，以提高对可燃毒物棒功率计算的精度。同时采用高阶燃耗方法，实现对各种可燃毒物燃耗效应的精细计算。

(3)

其中：P′m，g为第g群中子在第m个区域内核反应发出的功率；m为空间区域索引号；κ′i为第i个核素单次裂变释放的能量；Ni为第i个核素的原子核密度；σf，i，g为第g群第i个核素的微观裂变截面；qi为第i个核素的中子俘获释能；σc，i，g为第g群第i个核素的微观俘获截面。

堆芯扩散计算采用变分节块方法，直接在三维空间内将中子通量密度进行正交基函数展开，如式(4)所示。利用变分原理获得相应展开系数矩的离散代数方程组，避免了横向积分，可完全消除控制棒尖齿效应，能直接给出节块内连续的中子通量密度分布。

(4)

其中：φg(r)为在空间位置r处第g群的中子通量密度；φi，g为第g群中子通量密度的第i阶展开矩；fi(r)为正交基函数系的第i阶基函数；f为fi(r)构成的向量函数；φg为φi，g构成的向量。

2.3 验证与应用

Bamboo程序经过了大量模块验证、子程序验证和系统验证，已成功应用于BEAVRS、Wats-Bar、M310、CNP300、华龙一号、AP1000和EPR等多种型号压水堆的堆芯物理计算。以AP1000堆芯为例，表3列出首次临界试验参数的计算结果，表4列出控制棒组积分价值的计算结果。与实际测量值相比，计算结果均满足核安全审评的要求，表明Bamboo达到了与AP1000核设计程序APA系统相当的计算精度。

表3 启动物理试验参数的计算结果Table 3 Calculation result of startup physical test parameter

表4 控制棒组积分价值的计算结果Table 4 Calculation result of control rod integral worth

3 轻水堆全堆一步法高保真物理程序X

3.1 程序简介

X程序[15]是面向轻水堆开发的基于一步法计算的反应堆物理计算程序。程序基于高分辨率精细化建模、全堆共振计算、全堆大规模并行输运计算及pin-by-pin物理热工耦合实现了全堆芯高保真物理计算，给出了反应堆全三维精细的功率分布。

3.2 理论方法

X程序包括基于ENDF/B-Ⅵ.8和ENDF/B-Ⅶ.0评价数据库经由Atlas程序加工制作的多群数据库和连续能量数据库。为能精细处理复杂的反应堆堆芯结构，采用构造实体几何(CSG)建模方法，可直接描述包括格架、堆芯围板、吊篮等堆内构件。

为实现准确高效的全堆芯共振计算，提出了全局-局部耦合共振计算方法，利用中子流方法获得每个燃料棒的丹可夫因子Cb：

(5)

其中：φf，1为真实几何通过特征线方法(MOC)计算的中子通量密度；φf，0为利用碰撞概率方法求解的孤立棒问题中子通量密度。利用丹可夫因子守恒，将全堆芯共振问题等效为一系列一维圆柱问题，通过共振伪核素子群方法求解等效后的一维圆柱问题，获得每个燃料棒的共振自屏截面。

考虑到堆芯的非均匀性主要集中在径向，输运计算采取了基于泄漏项分割技术的2D/1D耦合方法，其中径向和轴向方程如式(6)、(7)所示。利用MOC求解二维方程、SN方法求解一维方程。

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

(6)

(7)

为提高计算效率，采用空间并行、角度并行和特征线并行三重并行技术，实现大规模并行计算。采用无矩阵粗网有限差分(CMFD)加速方法，通过扫描方式计算矩阵向量积，实现CMFD方程的高效求解和有效加速。

X程序拥有精细燃耗数据库(含有1 547个核素)和经过优化压缩的燃耗数据库(包含233个核素)以满足不同应用需求。利用切比雪夫有理近似方法(CRAM)进行点燃耗方程的求解。在此基础上通过预估校正、子步方法实现输运计算和燃耗计算之间的耦合。基于棒内子区尺度实现了精细的物理计算和pin-by-pin子通道热工计算的耦合，结合内部耦合和外部耦合的混合耦合方式实现了耦合数据的高效传递。

3.3 验证与应用

利用CASL计划的VERA系列基准题对X程序进行了详细的对比验证。计算基于全堆精细几何建模，如图2、3所示，直接给出了每根燃料棒的功率。结果表明：棒功率最大相对偏差和均方根偏差为：二维组件分别在1.5%和1.0%以内，二维堆芯分别在2.5%和1.0%以内。表5列出VERA系列基准题中代表性的VERA5三维堆芯算例的详细结果。由表5可见，三维堆芯特征值偏差在100 pcm以内。组件径向归一化功率最大相对偏差为1.38%，均方根偏差为0.69%。

图2 堆芯布置Fig.2 Reactor core layout

图3 精细几何建模Fig.3 High-fidelity geometry modeling

表5 VERA5三维堆芯计算结果比较Table 5 Calculation result comparison of 3D core of VERA5

将X程序应用到大型商用压水堆Watts Bar热态满功率的计算分析中，获得了pin-by-pin分辨率的燃料棒温度分布、冷却剂温度分布，如图4所示。与基于蒙特卡罗的耦合程序系统MC21/CTF相比，组件径向归一化功率的最大相对偏差和均方根偏差分别为1.7%和0.8%。

图4 燃料棒和冷却剂温度的三维分布Fig.4 3D temperature distribution of fuel rod and coolant

4 快谱反应堆物理计算程序SARAX

4.1 程序简介

SARAX程序[16-17]是针对快谱反应堆和局部装载固体慢化剂的中间能谱反应堆研发的物理计算程序，由快谱/中间能谱专用核数据库、截面产生程序TULIP、三维稳态堆芯程序LAVENDER组成。

4.2 理论介绍

SARAX程序基于ENDF/B-Ⅶ.0评价数据库经由Atlas程序加工制作的点截面数据库和超细群数据库。为处理快堆复杂的共振效应，采用基于点截面的方法计算有效自屏截面。超细群采用1 968群的能群结构，在每个能群中采用插值的方式产生200个能量点的截面数据。鉴于在中高能量区间，中子能谱很好地满足窄共振近似，因此超细群的有效自屏截面可由式(8)获得：

(8)

为处理组件的非均匀效应，基于Wigner单项有理近似，假设在1 968群各能群内的逃脱截面为常数，则非均匀系统下的有效自屏截面通过式(9)计算：

(9)

为准确计算快堆堆芯强烈的输运效应，本文采用非结构几何的SN节块方法求解中子输运方程以获得堆芯中子通量密度分布：

Σtψm(x,y,z)=Qm(x,y,z)

(10)

4.3 验证与应用

为验证SARAX程序的正确性及对于钠冷快堆计算的适用性，分别采用数值基准题和反应堆物理实验基准算例的测量数据对程序进行了验证与确认。在验证方面，主要基于IAEA及OECD发布的各类快堆基准算例。数值结果表明，对于各基准题算例，与蒙特卡罗程序计算的参考解相比，keff计算偏差在100～300 pcm范围内，功率相对偏差在3%范围内，反应性系数相对偏差不超过7%。在确认方面，针对日本长阳实验堆(JOYO)、美国零功率物理试验堆(ZPPR)及俄罗斯大科学装置(BFS)等物理试验进行了建模计算。计算结果表明，各临界试验的计算结果标准差为170 pcm，各类反应性计算结果计算的标准差小于5%，燃耗效应计算相对误差小于10%。

在此基础上，将SARAX程序应用于法国凤凰快堆寿期末物理试验与中国实验快堆(CEFR)启动物理试验的模拟。表6、7分别列出凤凰快堆提棒临界试验和CEFR首次临界试验的模拟结果。数值结果表明，SARAX程序对凤凰快堆的堆芯临界计算误差小于200 pcm，达到甚至优于法、俄专用核设计程序的精度。对CEFR的临界计算误差小于100 pcm。

表6 凤凰快堆控制棒提棒临界的模拟结果 Table 6 Simulation result of critical state at control rod withdrawal of Phoenix reactor

表7 CEFR净堆临界的模拟结果Table 7 Simulation result of critical state of CEFR

5 结论

本文介绍了西安交通大学反应堆物理团队开发的、具有完全自主知识产权的核反应堆物理确定论计算程序系统NECP软件包的核心理论模型和主要验证及应用情况。该软件包是在西安交通大学20多年来堆用核截面处理、共振理论、中子输运计算方法等方面研究的基础上，面向现有压水堆技术改造和新型反应堆研发需求开发的通用化确定论程序系统。软件包采用先进的理论模型和现代软件工程开发架构，形成了标准化、模块化的计算体系。所开展的大量验证和确认结果表明，该软件包功能齐全、计算精度高，可满足以压水堆为代表的热谱堆和以钠冷快堆为代表的快谱堆的计算要求，同时对于压水堆具备了高保真物理计算的能力。

NECP软件包已应用于国内大型压水堆项目、中国示范快堆项目等，国内用户基本包括了各大核电公司和研究院所，国外用户已达9个。已成功应用于AP1000启动物理试验的计算，预测的临界硼浓度与测量值相比误差仅为30 ppm；应用于CEFR启动物理试验的计算，预测的临界堆芯装载与实际测量值相比小于30 pcm。应用情况表明，NECP软件包的计算精度达到甚至优于国外先进的核设计专用程序。