数值反应堆原型系统开发及示范应用研究进展

杨 文，胡长军，刘天才，吴明宇，王昭顺，姜 强

(1.中国原子能科学研究院，北京 102413；2.北京科技大学，北京 100083)

更安全高效的核反应堆对快速精准的堆芯多物理耦合理论计算、燃料和材料服役老化行为及剩余寿命预测评估等提出了迫切需求。基于先进耦合建模和大规模并行计算技术的数值反应堆(简称数值堆)已成为国内外核反应堆工程与技术领域前沿热点。数值堆不仅有望使上述需求成为可能，更可能为先进反应堆的设计优化、不同工况运行模拟优化、严重事故序列演示预测及燃料和材料研发提供经济、高效的数值试验验证平台。美欧等核能先进国家高度重视数值堆技术研发与应用推广[1]。

近些年，在CASL[2-3]、NEAMS[4-5]、NURESIM[6-7]等系列大型项目持续支持下，美欧在数值堆技术研发方面取得了明显的研发进展及初步工业应用。美国于2010年部署立项了CASL、NEAMS等大型多部门联合研发项目。CASL项目致力于发展先进模拟与仿真技术，以改善在役轻水堆电厂的性能和效率，包括延长寿命、提升功率和增加燃耗，着重于压力容器内部的物理现象模拟。最终，CASL推出了一套反应堆分析程序包VERA，并于2020年取得商业应用许可证[8]。NEAMS项目则主要面向四代快堆、小型模块堆和其他新型反应堆，开发反应堆和核燃料循环系统先进数值模拟技术，重点开发了燃料产品线FPL、反应堆产品线RPL、集成产品线IPL 3个产品系列[9]。近期，美国又立项支持了第3个数值堆国家级研发项目CESAR[10]，其核心目标是在CASL和NEAMS基础上，升级开发可在正在开发的E级超算上部署运行的确定性中子输运(UNIC)、Monte Carlo中子输运(OpenMC)、热工CFD(Nek5000) 3个软件系统，最终形成耦合的下一代反应堆仿真工具包(TRIDENT)，能在E级平台高效执行，用于快堆设计、许可和安全分析等。欧洲数值堆的研发主要体现在NUR系列项目上。NUR旨在建立一个供欧洲核反应堆仿真通用的参考平台，提供高精度的物理模型、软件工具和评价结果，包括NURESIM(2005—2008)、NURISP(2009—2012)[11]、NURESAFE(2013—2016)[7]及NURE-NEXT等阶段，是一整套适用于现有轻水堆及未来堆型反应堆设计、安全分析及运行应用程序的开发和验证。

我国三大核电集团(中国核工业集团有限公司、中国广核集团有限公司、国家电力投资集团公司)及部分高校等紧跟国际前沿，近几年积极开展数值堆技术研发，并取得一定进展。在科技部“十三五”国家重点研发计划“数值反应堆原型系统开发及示范应用”重点专项支持下，中国原子能科学研究院联合北京科技大学、中国科学院计算机网络信息中心等积极开展了数值堆关键技术研发。目前初步完成了基于大规模并行的数值反应堆原型系统CVR1.0主要模块的开发，正在开展示范应用推广工作。本文介绍数值反应堆原型系统、核心软件开发与验证、反应堆物理-热工-结构多物理耦合计算，以及数值反应堆原型系统初步示范应用等研究进展。

1 数值反应堆原型系统CVR1.0

CVR1.0是核反应堆堆芯典型稳态运行工况数值模拟计算系统，可实现稳态运行工况下反应堆中子物理、热工水力、结构力学以及反应堆燃料和材料等典型物理过程和失效行为的多物理耦合高精细模拟计算和分析预测。其核心是一套基于E级超算的涉及反应堆中子物理、热工水力、结构力学、燃料和材料5个专业的高性能模拟计算软件，采用多物理耦合与软件测试架构(表1)以及“反应堆物理、热工、力学多尺度精细化耦合计算技术”“材料辐照脆化和辐照肿胀多尺度模拟计算技术”“基于E级计算机的超大规模并行求解技术”“多源数据、多模型、多物理装置结合的模拟验证技术”4项关键技术。

表1 数值核反应堆原型系统CVR1.0Table 1 Virtual reactor system CVR1.0

2 CVR1.0核心专业软件开发与验证

2.1 中子物理计算软件研发

1) 中子输运-特征线计算软件ANT-MOC

ANT-MOC是基于三维特征线法的堆芯中子输运计算软件，软件采用直接三维特征线方法，攻克了“特征线法多堆型高精细建模与预处理技术”“面向E级架构的并行特征线法优化技术”等关键技术，具备稳态中子输运方程求解能力[12]。

ANT-MOC通过了C5G7、Takeda等多个国际基准题验证。表2和图1分别为ANT-MOC求解C5G7基准例题[13]的有效增殖因数keff和中子通量密度分布计算结果，3种配置下keff与Monte Carlo程序计算值的相对偏差分别为0.140%、0.175%、0.270%，符合较好。目前，ANT-MOC已实现了中国实验快堆(CEFR)活性区127组件模拟，正积极推广其他应用。

图1 C5G7基准测试中子通量密度分布Fig.1 Neutron flux density distribution of C5G7 benchmark test

表2 C5G7算例测试结果Table 2 Test result of C5G7

2) 中子输运-Monte Carlo软件ANT-RMC

ANT-RMC是基于Monte Carlo方法的堆芯中子输运计算软件，由项目团队与清华大学合作开发完成。ANT-RMC能处理复杂几何结构、采用连续能量点截面对复杂能谱和材料进行描述，并能根据实际问题的需要对临界问题本征值和本征函数计算、精细核素链燃耗模拟、中子动态参数计算、在线核截面处理、核热耦合等进行计算。支持燃耗区超百万、计数超千亿的大规模并行算法。

目前，ANT-RMC通过了多个国际基准题(VERA#3基准题、VERA#5基准题、VERA#8基准题、VENUS-Ⅱ基准题和深穿透基准题等)的验证。表3为ANT-RMC求解VERA#5全堆基准题算例2[14]的keff计算结果，和MCNP的计算结果符合较好。

表3 VERA#5测试结果Table 3 Test result of VERA#5

2.2 热工水力模拟软件研发

1) 子通道分析软件CVR-PASA

CVR-PASA是基于数据库进行数据管理的高性能并行子通道分析软件，采用MPI+OpenMP混合编程模型实现，使用区域分解方式进行并行任务划分，在划分时采用“适应于多几何堆芯类型的自适应全堆芯子通道任务划分与映射”关键技术划分全堆芯子通道任意求解域个数。CVR-PASA可实现全堆芯、全通道热工水力模拟，可获得稳态流场、压力场、温度场以及空泡份额分布和临界热流密度。

CVR-PASA通过了OECD/NRC标准题等算例、压水堆实堆数据的验证，在天河2号超算上进行压水堆实堆全堆芯157组件(控制体约600万)的稳态模拟，在8 792核时模拟时间为354.66 s，同类型算例CTF需约23 min。表4和图2是CVR-PASA模拟商用压水堆实堆全堆芯的冷却剂压力等和燃料棒包壳外表面温度分布计算结果，与实堆数据基本一致，能反映各参数在堆内的真实分布。

表4 CVR-PASA计算结果与实测数据对比Table 4 Comparison of CVR-PASA results with measured data

图2 燃料棒包壳外表面温度分布Fig.2 Temperature distribution on the outer surface of fuel rod cladding

2) 计算流体力学模拟软件CVR-PACA

CVR-PACA是数值堆单相精细化计算流体力学热工水力模拟软件。软件采用精确大涡模拟模型的高精度谱元方法(最高可达24阶精度，常用商用软件一般不超过5阶精度)，支持异构架构，在“神威·太湖之光”超算上采用主从核异构并行方案。CVR-PACA可对反应堆的堆芯进行精确的数学物理建模，支持模拟的网格数量超百亿，选用合适的湍流模型来闭合平均质量、动量和能量守恒方程，精确求解反应堆容器内部的压力场、流场、温度场[15-16]。CVR-PACA通过了Matis-H等多个基准题验证。从Matis-H基准题[17]撕裂式搅混翼下游不同位置处各时均速度分量与实验值的对比[15]可看出，CVR-PACA的LES模型及URANS SSTk-ω模型均与实验值较接近，能较好预测由于撕裂式搅混翼引起的湍流流动状态。

2.3 结构力学模拟软件研发

HARSA是用于反应堆堆芯结构力学行为数值模拟的软件，实现大规模撕裂有限元求解方法及负载均衡技术、性能分析方法，支持静力学、流致振动、磨损等计算，具备全堆芯组件静力变形、流致振动分析能力[18]。

HARSA通过了IAEA基准算例1、单根组件受限热弯曲试验算例、IAEA基准算例6等算例验证，在天河2号上，模拟网格量达137亿。表5和图3是使用HARSA按照IAEA基准算例1的要求建立单盒组件自由变形模型，开展有限元静力计算的计算结果。与IAEA基准算例1给出的理论解对比，HARSA计算结果的偏差在3.8%以内，计算准确性较好。

表5 组件轴向伸长量计算结果对比Table 5 Comparison of fuel assembly axial elongation deformation displacement

图3 组件轴向伸长形变位移模拟Fig.3 Simulation of fuel assembly axial elongation deformation displacement

2.4 燃料性能分析软件研发

ATHENA是用于对压水堆燃料元件进行单棒和多棒性能分析的软件，软件开发了燃料温度分布、变形计算、裂变气体释放及内压等模型，结合燃料元件热工-力学多物理耦合计算分析耦合方案，基于先进并行计算方法，建立了燃料并行分析能力。支持计算燃料元件温度分布、变形、裂变气体释放、燃料棒内压、包壳腐蚀等参数，具备高效、精细化全堆芯燃料性能分析的能力[19]。

ATHENA通过国际标准基准例题、同类程序计算结果、典型商用压水堆核电站数据进行了验证，获得了燃料温度、燃料变形、燃料内压等参数计算对比结果，分析表明ATHENA软件计算结果合理可靠。图4为ATHENA和对标软件METEOR基于典型商用压水堆核电站燃料数据计算得到的燃料芯块中心温度随时间变化的趋势图，计算结果基本符合。

图4 峰值温度随时间变化Fig.4 Peak temperature change with time

2.5 材料多尺度模拟软件研发

MISA是项目团队自主研发的系列反应堆用结构材料辐照损伤多尺度模拟软件集，包括微观尺度的分子动力学模拟软件MISA-MD、空位跃迁机制的动力学Monte Carlo模拟软件MISA-KMC、空位和间隙双跃迁机制的动力学Monte Carlo模拟软件MISA-AKMC[20]、介观尺度的相场模拟软件MISA-PF、随机团簇动力学模拟软件MISA-SCD。软件集可对核材料辐照损伤进行多尺度模拟，获得用于RPV钢辐照脆化、堆内构件辐照肿胀预测所需的微观信息。

MISA系列软件主要通过国外软件模拟算例及结果进行验证，模拟结果与对比软件较吻合。图5分别为对LAMMPS和MISA-MD模拟结果作统计分析得到的Frenkel缺陷对数量Nfp随时间变化趋势图，结果较吻合。

图5 级联碰撞过程中Frenkel缺陷对数量随时间变化图Fig.5 Evolution of number of displacement cascades’ Frenkel defect with time

3 多物理耦合系统与数值模拟软件测试框架研发

3.1 多物理耦合系统研发

数值堆多物理耦合模拟平台SPIDER是一款能实现核反应堆高精细耦合模拟的集成软件平台，其整体架构如图6所示，支持各物理模块之间的不断交互和统一协作，并针对各模块设计建立数据库，用于促进核反应多物理耦合过程中的数据交互和模拟数据的统一管理，具备数据支撑能力。

图6 SPIDER系统架构Fig.6 System structure of SPIDER

SPIDER平台中包括核热耦合和流固耦合两大核心模块。核热耦合过程十分复杂，存在着强烈的相互反馈，一方面，通过物理计算可得出裂变反应率和中子通量分布，从而计算出功率分布，并进一步求得反应堆的传热特性和热工水力特性；另一方面，热工水力的计算会对物理过程进行反馈，影响中子通量分布。SPIDER平台中核热耦合模块针对自主研发的并行子通道模拟软件CVR-PASA、特征线法中子输运模拟软件ANT-MOC设计实现了基于文件的数据交互耦合流程，为反应堆物理-热工耦合的数值模拟提供了更为准确的堆内参数变化。流固耦合模块针对自主研发的计算流体力学软件CVR-PACA和结构力学软件HARSA同样设计并实现了基于文件的数据交互耦合流程[21]。

SPIDER通过了C5G7、中国实验快堆燃料组件、Hoogenboom-Martin等多个算例验证[22]。图7为采用OECD的C5G7-Rodded B例题[13]进行核热耦合功能验证的计算结果，燃料棒表面温度分布合理，与文献符合较好。采用中国实验快堆燃料组件算例进行流固耦合功能验证的计算，燃料棒束最大位移出现在中心燃料棒束的顶端，为27.83 mm[21]。

图7 燃料棒表面温度分布Fig.7 Temperature distribution of fuel rod surface

3.2 数值模拟软件测试框架研发

数值模拟软件测试框架(NuSTA)是为数值模拟程序提供一种可用的测试框架，可生成算例，进行单元测试、蜕变测试、属性测试、精度测试等。框架的主要功能为：统一测试用例定义规范，支持多种不同输入输出格式的高性能数值模拟软件；支持Verification & Validation流程中的主要活动，实现传统测试、蜕变测试、精度阶分析、回归测试、误差分析等功能接口的集成和封装；完成蜕变测试、误差分析工具(ScDebug)的开发；实现测试数据的管理和展示。框架的技术特点在于：将蜕变测试和数值模拟程序相结合，有效缓解测试用例少且难构造的问题；研发ScDebug精度测试工具，相比目前动态精度分析工具性能提高70%。NuSTA框架的总体结构设计如图8所示。

图8 NuSTA框架总体结构Fig.8 Overall structure of NuSTA framework

4 数值反应堆原型系统示范应用

完成CVR1.0的开发和验证后，目前正开展相关示范应用工作，重点在关键材料失效分析和反应堆工程设计或设计验证。

RPV钢辐照脆化多尺度模拟研究，主要利用自主开发软件针对RPV钢及其模型合金的辐照脆化机理开展研究。项目团队利用MISA-KMC软件研究了合金元素Mn对富Cu团簇形核的影响，分析了Mn含量变化对团簇析出的影响[23]；利用位错动力学方法结合分子动力学和分子静力学计算获得的缺陷钉扎力，研究了FeCu模型合金中Cu析出物导致硬化的机理，该工作将微观结构演化模拟与宏观性能预测建立桥梁，对深入研究RPV钢辐照硬化机理以及预测辐照脆化趋势具有重要意义[24]。

基于MISA-MD与MISA-SCD耦合模拟了Fe0.3at.%Cu合金中微观结构演化及其导致的辐照硬化、脆化行为，图9为Fe0.3at.%Cu合金经中子辐照后微观结构的数密度和平均半径随损伤剂量(dpa)的变化，辐照后产生的微观结构主要为纯Cu团簇、Cu_Vac复合团簇、孔洞和位错环。位错环的直径小于1 nm，接近0.75 nm，可认为位错环对位错运动阻碍作用较弱；孔洞的直径约为2 nm，但孔洞在演化后期才出现，且数密度很低，可忽略其对位错运动阻碍作用；Cu_Vac类型团簇尺寸较大，但辐照前期Cu_Vac数密度较低，可不考虑其对位错运动阻碍作用。当仅考虑纯Cu团簇对辐照硬化、脆化的影响时，基于Orowan硬化模型及RPV钢辐照后硬化与脆化关系，可获得纯Cu团簇引起材料韧脆转变温度增量随损伤剂量变化关系。图10为纯Cu团簇引起材料韧脆转变温度增量(DBTT)与实验结果对比，由图可知，模拟结果与实验结果[25]较吻合，验证了结果的可靠性。

图9 中子辐照下Fe0.3at.%Cu微观结构半径和数密度随损伤剂量演化关系Fig.9 Relationship of radius microstructure and number density of Fe0.3at.%Cu and damage dose under neutron irradiation

图10 韧脆转变温度增量随损伤剂量变化关系Fig.10 Predicted increasement of DBTT as a function of damage dose

5 总结与展望

在国家重点研发计划“数值反应堆原型系统开发及示范应用”重点专项支持下，初步完成了基于大规模并行的数值反应堆原型系统CVR1.0主要模块的开发，并在反应堆压力容器辐照脆化机理认知与辐照脆化趋势预测、快堆设计分析等方面实现了初步示范应用。

CVR1.0核心是基于E级超算的反应堆中子物理、热工水力、结构力学、燃料和材料5个专业高性能模拟计算软件，可实现反应堆堆芯稳态运行工况下中子物理、热工水力、结构力学以及反应堆燃料和材料等典型物理过程和失效行为的多物理耦合高精细模拟计算和分析预测。

数值堆后续研发重点是持续推进CVR1.0应用反馈和升级优化的同时，进一步拓展CVR1.0至全堆系统、全工况运行模拟，并逐步实现反应堆设计研发与运行工程应用。