CAP1400数值反应堆系统关键技术研究及示范应用

曹良志，邓 力，杨 波，刘宙宇，刘 鹏4,，汤春桃，史敦福4,，陈荣华，田文喜，彭良辉，万承辉，张旻婉，毕光文，费敬然，许晓北，李 帆

(1.西安交通大学 核科学与技术学院，陕西 西安 710049；2.北京应用物理与计算数学研究所，北京 100088；3.上海核工程研究设计院有限公司，上海 200233；4.中物院高性能数值模拟软件中心，北京 100088)

数值反应堆技术通过高精度、高分辨率、高置信度的高保真数值模拟方法，基于高性能计算和多物理耦合实现对核反应堆堆芯内部多种物理现象的刻画，从机理上揭示核反应堆内各种物理现象的内在规律，准确预测核反应堆在服役过程中的关键安全参数，对于核反应堆的研发、设计、安全分析、运行支持和退役的全生命周期都具有重要作用。随着高性能计算技术的快速发展和公众对核能安全性与经济性要求的不断提高，数值反应堆技术在世界核发达国家均受到了广泛重视，如美国的I-NERI[1]、CASL[2]和NEAMS[3]等项目，欧盟的NURE系列计划[4-6]、CORTEX[7]和HPMC/McSAFE计划[8]等项目。国内也通过科技部等相关部门设立了数值反应堆项目研究[9-10]。

在“大型先进压水堆和高温气冷堆核电站”国家科技重大专项的支持下，由北京应用物理与计算数学研究所牵头，上海核工程研究设计院有限公司和西安交通大学参与的课题“CAP1400数值反应堆关键技术研究”，重点针对我国自主三代核电“国和一号”(CAP1400)的设计验证和关键参数预测，开展了基于全堆一步法输运计算的堆芯物理确定论程序NECP-X、基于蒙特卡罗方法的高保真输运计算程序JMCT、基于pin-by-pin先进子通道模型的热工水力程序SACOS、基于精细网格下多物理耦合的先进数值燃料性能分析程序YUAN、基于直接数值扰动和统计抽样方法的不确定度量化分析程序UNICORN的研发，并将所研发的软件进行集成，利用国际基准题VERA、我国AP1000核电厂实际运行测量值对数值反应堆系统进行验证与确认。在此基础上，针对CAP1400堆芯设计方案，对首循环堆芯启动物理试验、功率运行跟踪、负荷跟踪及主蒸汽管道破裂事故典型工况开展高保真模拟分析，为CAP1400的调试运行提供重要的参考依据。

1 CAP1400数值反应堆系统关键技术研究

CAP1400数值反应堆系统的核心体系包括核反应堆物理、反应堆热工水力、燃料性能和不确定性分析等4个领域，并在此基础上实现CAP1400数值反应堆的多物理耦合集成开发、不确定性分析和示范应用。

1.1 基于确定论方法的高保真堆芯物理计算关键技术研究与程序开发

课题攻克了基于确定论方法的高保真堆芯物理关键技术，研制了基于确定论方法的高保真堆芯物理计算程序NECP-X[11-12]。基于ENDF/B-Ⅶ.0评价核数据库制作了堆用连续能量核数据库和多群核数据库，并针对钨控制棒的截面进行了全面验证。在全局-局部耦合共振计算技术的基础上，针对CAP1400的机械补偿运行模式，研究了针对控制棒的超栅元共振计算技术，并攻克了2D/1D耦合输运计算处理控制棒连续移动的难题。最终开发了NECP-X程序稳态计算、反应性计算、循环燃耗计算功能，实现了针对CAP1400堆芯特点的确定论全堆芯高保真堆芯物理计算能力。

程序实现了三维全堆芯精细几何建模和万核以上大规模并行计算的能力，1/4堆芯全堆芯高保真共振计算和输运计算总计约200核时，天河二号超算的万核并行效率约为70%，在国内首先实现了基于确定论方法的全堆芯一步法非均匀输运精细计算。程序经过VERA基准题、BEAVRS基准题等大量的国际知名基准题验证，数值结果表明程序具备很高的计算精度与效率，并成功应用于第3代先进压水堆AP1000和华龙一号、第二代改进型M310压水堆、西安脉冲堆、JRR-3M、高通量研究堆等堆型的模拟计算，获得了反应堆在服役过程中的有效增殖因数、功率、温度等关键安全参数随时间的演变情况，M310堆芯第1循环的临界硼计算结果误差在燃耗寿期内均在30 ppm以内，图1为不同燃耗下的精细功率分布，与测量值相比，最大组件功率相对偏差在5%以内。大量的验证结果表明，NECP-X可为商用大型压水堆、研究堆等各类堆型的堆芯中子学设计提供可靠的工具。

a——0 EFPD；b——304.1 EFPD

1.2 基于蒙特卡罗方法的高保真堆芯物理计算关键技术研究与程序开发

基于蒙特卡罗模拟软件JMCT研制了蒙特卡罗方法的CAP1400数值反应堆物理计算程序。JMCT基于并行组合几何编程框架JCOGIN进行研发，依靠JCOGIN框架高性能的数据结构和高效成熟的通信算法，使程序具备粒子数超千亿、计数超千亿的大规模并行计算的能力。JCOGIN框架采用MPI+OpenMP的混合并行模式，支持数十万核规模的并行。JMCT利用BEAVRS模型在天河二号服务器进行了测试，12万核并行效率可达80%[13]。

JMCT以ENDF/B-Ⅶ为基础制作ACE格式的多温度连续能量核参数库。针对CAP1400反应性计算及核热耦合计算对截面的要求，JMCT建立了两套多温度策略以支持多温度问题的模拟，一是离散温度截面选择(discrete temperature crosssection select, DTCS)策略，二是在线多普勒展宽(on the fly Doppler, OTF)策略。DTCS策略是根据应用需求预先制作多个温度下核素的参数截面，并在JMCT运行过程中选择与真实温度最近的截面加以使用。由于真实模型的温度是连续的，DTCS一般会存在一些误差，但其计算速度快、内存开销小，主要用于截面对温度变化不太敏感的核素，如Fe、Zr等结构材料。OTF策略是在输运过程中根据真实模型温度获得对应的截面，因此对温度的处理是连续的，模拟精度很高，针对核反应堆计算时，只对温度变化比较敏感的核素进行OTF处理。对于热散射截面，在300～1 000 K范围内设置多个温度点，每50 K为温度间隔，程序自动寻找最近温度点的截面数据。

为了支持CAP1400的数值模拟，建立了硼调节和控制棒调节两种反应性调节模式。临界硼浓度搜索采用了基于中子平衡的搜索算法，仅需要一次临界计算，并且避免了统计方差对插值造成的影响，精度较高[14]。对于控制棒调节功能，创新地建立了“延迟材料替换”算法，在粒子输运径迹计算过程中实现棒位的虚拟位移，在不影响几何间的拓扑关系的条件下，解决了核反应堆的控制棒移动问题。

JMCT的堆芯计算功能通过BEAVRS热态零功率模型[15]、NURISP项目的迷你堆模型和VERA基准题进行了较为充分的验证。并且在堆芯的设计和计算中获得了初步应用，如宁德核电厂1号机组的换料周期计算[16]，如图2所示。

图2 宁德核电厂1号机组堆芯首循环模拟结果

1.3 pin-by-pin先进子通道热工水力计算关键技术研究与程序开发

针对数值反应堆系统高保真计算的需求，开发了基于两流体模型的全堆芯先进子通道分析程序SACOS。pin-by-pin全堆芯计算的关键难题是求解全局耦合的压力矩阵，如数百个组件的压水堆堆芯，考虑轴向和径向划分，全局压力矩阵系数的维度将达到百万级，计算时间和计算内存要求非常高。项目克服了大型压力矩阵的大规模并行求解难题，实现了全堆芯pin-by-pin的高效计算。

SACOS由西安交通大学热工水力实验室(XJTU-NuThel)开发,采用成熟的两流体八方程模型[17]，相较于均相流模型(COBRA-EN等)和三流场模型(CTF等)，在保证模型精度的条件下计算效率更高且模型不确定因素更少。同时程序采用带非线性通量限制器(如OSPRE格式)的高阶稳定离散格式，相比于常规的一阶向上迎风格式，在数值求解中提高了计算精度且减少了数值扩散程度。对于大型压力矩阵求解，程序采用基于MPI的空间区域分解策略，通过划分的并行域可将计算负担分到各个CPU上，各个CPU只需产生其负责的区域的系数矩阵并分块存储。形成压力方程的大型系数矩阵后，调用并行矩阵求解器实现快速求解。与国际上先进的并行子通道程序CTF和ESCOT相比，摒弃了与组件严格对应的径向区域分解，采用更灵活可扩展程度更高的并行域划分方式[18]，图3为全堆芯划分87个计算域由87个计算核心参与计算。可实现全堆芯、全通道高效率的热工水力并行模拟，可获得稳态和瞬态流场、压力场、温度场以及空泡份额分布和临界热流密度。基于典型压水堆堆芯布置，利用SACOS在高性能计算集群上实现了总控制体数超过百万的全堆芯子通道稳态工况模拟，如图4所示。随着进程数的增加，计算速度明显加快，252个处理器的加速比可达到180，计算时间可减少到3 min以内，表明SACOS具有高效的并行效率，可以满足数值反应堆的堆芯pin-by-pin精细计算需求。

图3 全堆芯径向并行域划分

图4 全堆芯pin-by-pin热工水力分析结果

1.4 核反应堆先进数值燃料棒关键技术研究与程序开发

数值燃料棒程序系统YUAN是一个燃料棒多物理数值模拟平台，主要功能是在燃耗过程中基于可变形的浮动网格耦合求解中子学、热工水力、机械力学问题，模拟燃料棒在反应堆内极端环境下的热传导、机械变形、中子输运、核素燃耗等行为。数值燃料棒适用于研究燃料行为及其演化机理，从多物理角度全面分析燃料棒性能。

利用YUAN对CAP1400首循环堆芯的典型燃料棒进行了计算分析。燃料棒基于正方形栅元，采用UO2芯块、锆合金包壳，间隙填充氦气，芯块直径、包壳内外径、栅距分别为0.819 2、0.835 7、0.95、1.26 cm，有效高度为426.72 cm。燃料棒径向温度分布、径向等效应力及径向功率分布计算结果如图5所示。

a——径向温度分布；b——径向等效应力分布；c——径向相对功率分布

1.5 数值反应堆的不确定度量化分析关键技术研究

基于直接数值扰动的敏感性和不确定性分析理论方法，在核反应堆数值计算敏感性和不确定性分析平台UNICORN[19]的基础上，针对数值反应堆多物理耦合计算程序，实现了核反应截面、制造参数等的敏感性和不确定性分析功能。针对传统抽样方法在样本容量有限时具有较大的误差、无法准确描述输入参数的不确定度信息的问题，研究了面向协方差矩阵的样本变换(COST)方法，克服了抽样过程中样本数少无法满足标准正态分布引入的误差，有效减少直接数值扰动方法进行不确定度量化的计算量，为中子学、热工水力及物理-热工-燃耗耦合的不确定度量化分析建立了基础。

针对中子学计算的不确定度量化分析，重点研究了核数据和制造参数不确定来源，其中制造参数的不确定性包括燃料内径和外径、包壳厚度、燃料芯块直径、燃料密度、气隙压力和燃料富集度，量化了这些因素对有效增殖因数、功率分布及燃耗过程中的核子密度、燃耗深度分布的不确定度。

热工水力分析程序基于简化的守恒方程，经验公式众多，且公式的有效范围有限。此外，边界条件和初始条件通常不是完全已知，具体包括：1)边界条件的不确定度，如系统压力、冷却剂流量、系统功率、入口冷却剂温度和功率分布；2)几何的不确定度，如燃料棒偏移和燃料棒外径；3)模型的不确定度，如单相混合系数、两相混合系数、传热系数、每种流态的阻力系数、定位格架损失系数和壁面摩擦模型等。本课题在热工水力方面分析了稳态和瞬态过程中燃料温度、包壳温度和压降的不确定度范围，如针对TM1组件问题[20]，瞬态过程最大冷却剂温度、最大包壳温度及其不确定度如图6所示。在初始状态下，包壳最大温度为(338.00±1.27)℃，冷却剂最大温度为(326.00±1.06)℃；在5.0 s时，包壳最大温度为(345.00±0.62)℃，冷却剂最大温度为(335.00±1.14)℃。

图6 瞬态过程中冷却剂和包壳最大温度及其不确定度(3σ)

2 CAP1400数值反应堆系统开发及验证

2.1 CAP1400数值反应堆系统多物理耦合集成开发

在上述核心模块的基础上，采用开源的版本控制软件、测试软件、软件性能分析软件等，进行核心软件的可信开发和开发过程的有效控制，并建立了并行网格映射、数据传递和高负载平衡的并行任务分配算法进行核心软件的耦合，实现了CAP1400数值反应堆系统多物理耦合的集成开发。为了提高数值反应度系统的使用友好性，进一步建立了数值反应堆系统可视化建模和后处理工具，完成数值反应堆的高效建模和数值结果的高效率可视化处理。其中JLAMT是中国工程物理研究院高性能数值模拟软件中心基于UG研发的可视化建模工具，采用基本体(CSG)的建模方式，可直接生成计算可用的几何文件，适用于蒙特卡罗、MOC等计算方法。针对压水堆和快堆的堆芯模型，JLAMT提供模型定制功能，可实现棒内精细结构建模、任意燃料棒布置和组件布置，可以任意轴向剖分，可快速生成千万量级几何单元的输入文件。基于开源引擎研发，实现了分布体系结构的可视化后处理功能，充分利用大规模并行机的计算能力及本地工作站的图形加速性能，实现全堆芯精细数值结果的可视化。CAP1400数值反应堆系统的结构如图7所示。

图7 CAP1400数值反应堆系统结构

2.2 CAP1400数值反应堆系统集成验证

在数值反应堆系统核心软件独立验证的基础上，课题针对数值反应堆集成系统开展了大量的验证和测试，重点分析了针对CASL计划系列基准题的验证情况和利用先进三代压水堆AP1000实测数据进行的确认情况。

1)VERA系列问题

CASL计划基于Watts Bar核电站，建立了包括单栅元问题、单栅格、单组件、全堆芯、首循环的一系列基准题[21]，能够比较系统地验证核反应堆数值模拟方法和程序在压水堆中的适用性和计算精度。此处给出了覆盖不同温度、不同控制棒、不同可燃毒物、是否包含格架等各种工况范围的典型算例结果，包括蒙特卡罗和确定论方法的数值反应堆程序的特征值、功率分布及与参考结果的比较情况。表1～3为单栅元、栅格和三维单组件的特征值结果，与KENO-Ⅵ程序[22]相比，除含AIC控制棒的栅格问题偏差为-106 pcm外，其余特征值偏差均在±100 pcm以内。表4、5为全堆芯不同控制棒棒位的特征值结果和温度系数的计算结果，棒价值和温度系数均与参考结果吻合很好。对于VERA#5算例2，JMCT在曙光I620-G30计算节点，Intel Xeon Gold 5118(2.3 GHz)核心，利用120核并行，采用200万粒子/代，共1 000代(其中300非活跃代)耗时24 417 s；NECP-X程序计算条件为0.03 cm特征线宽度，192个角度方向，基于天河2号超算(Intel Xeon E5-2692 v2(2.2 GHz)核心)利用2 056核并行，耗时408 s。

表1 VERA#1单栅元基准题特征值计算结果

表2 VERA#2二维单组件基准题特征值计算结果

表3 VERA#3三维单组件基准题特征值计算结果

对于核-热-燃耗耦合计算，图8给出了热态满功率0燃耗的全堆芯的功率分布比较结果，与MC21/CTF耦合程序的结果对比表明，无论是基于蒙特卡罗方法还是基于确定论方法，CAP1400数值反应堆系统的功率分布相对偏差均在2.5%以内，具有很高的计算精度。图9、10给出了VERA#9首循环CAP1400数值反应堆程序与临界硼浓度测量值的偏差和全堆芯pin-by-pin功率分布，最大临界硼浓度偏差小于30 ppm。

图9 VERA#9问题第1循环的临界硼结果

表5 VERA#5控制棒价值和反应性计算结果

a——基于确定论方法；b——基于蒙特卡罗方法

CASL计划系列基准题的数值结果表明，基于确定论和蒙特卡罗方法的CAP1400数值反应堆系统均具有可靠的计算精度。

a——0 EFPD；b——194.3 EFPD；c——392.3 EFPD

利用CAP1400数值反应堆物理-热工-燃料性能耦合计算功能，对VERA#9问题进行了多物理耦合计算，获得了随燃耗变化的燃料和包壳应力、应变，燃料和包壳之间的间隙宽度。图11给出了堆芯每个燃料组件内最大包壳米塞斯应力随燃耗的变化，数值结果与理论预期吻合较好。

a——0 EFPD；b——194.3 EFPD；c——392.3 EFPD

2)AP1000压水堆

AP1000压水堆是由美国西屋公司研发的非能动安全压水堆，首循环堆芯采用低泄漏设计，使用从天然铀至富集度为4.8%的5个燃料区。采用MSHIM控制策略，即在正常运行和变功率阶段，主要通过调节钨控制棒组和银-铟-铬控制棒组，再配合极少的堆芯冷却剂硼浓度变化实现堆芯反应性控制和功率分布调节。同时AP1000采用了两类可燃毒物棒，即将薄的ZrB2涂层沉积到燃料表面的整体式可燃毒物棒IFBA及湿环形离散式可燃毒物棒WABA，对于离散式WABA而言，根据毒物长度的不同分为长、中和短WABA。AP1000堆芯的布置如图12所示[23]，为了增加经济性，实现低泄漏首循环，同时考虑到展平功率分布和反应性控制，堆芯径向采用了五区装料方案。按照燃料富集度由低至高分为A、B、C、D及E区，其富集度依次为0.73%、1.58%、3.20%、3.78%和4.38%。总体而言，相对现役的二代和二代改进型压水堆，其堆芯设计更加精细，给数值模拟带来了更大的挑战。

图12 AP1000首循环1/4堆芯布置

使用数值反应堆程序对AP1000堆芯进行精细的几何建模，如图13所示，并计算了AP1000启动物理试验参数，包括临界硼搜索、硼微分价值及各类温度系数，如表6、7所列。可以看出CAP1400数值反应堆系统的临界硼搜索计算结果、硼微分价值及各类温度系数计算结果与测量值符合得很好。控制棒价值与测量值的结果比较如图14所示，无论是灰棒价值还是黑棒价值均与测量值符合很好。以机组A为例，以实测值为参考，工业级程序、确定论方法和蒙特卡罗方法的数值反应堆程序偏差如下：1)灰棒价值最大偏差分别为-10.3 pcm、7.3 pcm和5.7 pcm；2)黑棒价值最大相对偏差分别为3.2%、-1.7%和0.8%；3)ITC偏差分别为-0.597 pcm/℃、0.280 pcm/℃和1.570 pcm/℃。数值反应堆程序与测量值的结果更加接近。热态零功率的功率分布偏差统计情况如图15所示，确定论方法和蒙特卡罗方法的数值反应堆程序功率分布偏差很小，相比而言工业级程序与蒙特卡罗方法的功率分布偏差更大，体现了数值反应堆程序具有高的计算精度。

图13 堆芯几何建模图

表6 启动物理试验临界硼浓度计算结果

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

图14 机组A(a)和机组B(b)的棒价值偏差

a——确定论数值反应堆程序与蒙特卡罗数值反应堆程序；b——工业级程序与蒙特卡罗数值反应堆程序

3 CAP1400数值反应堆系统示范应用

在工程应用方面，数值反应堆的模拟分辨率优于工业级程序，数值反应堆在核设计程序的验证、设计结果的校核及新型反应堆设计的分析等方面均具有重要应用价值。完成数值反应堆程序的开发和验证后，针对CAP1400开展了示范应用，对CAP1400首循环堆芯启动物理试验、功率运行跟踪、负荷跟踪及主蒸汽管道破裂事故典型工况开展了高保真模拟分析。

3.1 启动物理参数计算分析

启动物理试验是反应堆运行调试的重要环节，试验过程中，若设计程序预测值与电厂实测值偏差过大，将直接影响电厂正常启动、升功率等过程，进而造成较大的经济损失。为防止上述问题，可使用数值反应堆对启动物理试验进行高保真模拟分析。使用CAP1400数值反应堆系统对CAP1400反应堆首循环堆芯启动物理试验参数进行了高保真模拟。图16给出控制棒组价值偏差统计结果。以蒙特卡罗方法数值反应堆程序计算结果为参考解，确定论数值反应堆程序对于热态零功率(HZP)控制棒全提(ARO)工况下的临界硼浓度的绝对偏差为-23.6 ppm；对于灰棒组(MA～MD)价值的最大绝对偏差为-5.2 pcm，对于黑棒组(M1、M2、AO、SD1～SD6)价值的最大相对偏差为2.68%；对于黑棒总价值的相对偏差为-0.08%。如图17所示，以蒙特卡罗方法数值反应堆程序为参考，对14种插棒工况(ARO和各组棒单独插入堆芯)下，NECP-X程序组件相对功率分布的计算偏差进行统计，其中相对功率偏差的最正值为2.24%，最负值为-2.47%，平均值为-0.06%，标准偏差为1.02%。

图16 CAP1400首循环堆芯棒价值偏差

图17 组件相对功率分布计算偏差统计

3.2 功率运行跟踪计算分析

反应堆功率运行工况是典型的中子学、热工水力学多物理耦合过程，使用数值反应堆系统对CAP1400首循环堆芯功率运行工况进行高保真模拟，以进一步检验设计结果的准确性。图18中，以工业级程序为参考，给出NECP-X临界硼浓度计算偏差随燃耗的变化，临界硼浓度的最大绝对偏差为27.6 ppm。图19中，以工业级程序为参考，给出NECP-X组件相对功率分布相对计算偏差的最正和最负值随燃耗的变化，各燃耗点组件相对功率分布相对偏差的最正值为3.14%，最负值为-2.93%。

图18 临界硼浓度计算偏差

图19 组件相对功率分布计算偏差

使用数值反应堆程序对CAP1400首循环堆芯进行了高保真模拟分析，给出了CAP1400首循环堆芯启动物理试验参数、功率运行参数的高保真模拟结果，进一步检验了设计结果的准确性，可有效支撑国和一号示范工程反应堆首循环物理启动。

4 总结及展望

CAP1400数值反应堆关键技术研究课题对核反应堆先进数值分析技术开展了深入的研究，突破了基于一步法全堆芯输运计算的确定论数值反应堆物理计算方法、具有控制棒调节功能的蒙特卡罗数值反应堆物理计算方法、具有pin-by-pin和横向搅混功能的子通道分析方法、基于棒内多物理耦合的数值燃料棒计算方法以及堆芯计算不确定性传递原理和量化方法等多项关键技术。依托高性能硬件平台研制了基于确定论方法和蒙特卡罗方法的高保真堆芯物理-先进子通道-燃料性能多物理耦合的CAP1400数值反应堆系统。研制的数值反应堆系统通过了国际基准题VERA的启动物理试验和循环模拟验证，并首次利用国内三代压水堆AP1000电厂实测数据实现了数值反应堆程序的确认，在此基础上，将该系统应用于CAP1400反应堆堆芯精细建模及模拟计算，完成了CAP1400低功率物理启动试验的模拟和部分的功率能力分析计算，为CAP1400设计验证和物理启动提供了重要参考数据，为确保三代非能动核电的安全性、提升经济性提供了重要的分析工具。研究团队将以该课题研发的数值反应堆系统为基础，进一步推动高保真数值模拟技术在核电、核动力等多种先进核反应堆中的验证和应用，为实现“碳达峰、碳中和”目标做出更大的贡献。