反应堆核-热-燃耗多物理耦合框架研究与应用

吴明宇，朱 迎，卢 旭，苗 雪，吴宗芸，李 龙，胡 赟，赵民富

(1.中国原子能科学研究院 反应堆工程技术研究所，北京 102413；2.北京科技大学 计算机与通信工程学院，北京 100083)

反应堆中子学、热学、流体力学是相互耦合的复杂系统，每个物理过程由偏微分方程或方程组描述。例如：中子物理过程由中子输运方程描述，燃耗过程由燃耗方程组描述，瞬态过程由中子动力学方程描述，冷却剂的换热流动过程由传热方程和流体力学方程组描述。上述方程中求解的物理量往往是互相影响、相互耦合的。例如，中子输运计算得到反应堆内的中子通量分布情况，是进行燃耗计算的必需数据，也为热工水力计算提供了相对功率分布；由热工水力输运计算得到的堆内冷却剂的密度和燃料温度分布，可通过冷却剂核素密度的变化及多普勒效应改变材料的宏观截面，从而对中子的输运过程产生影响；燃耗计算得到的新的燃料核素及核素密度，不但会通过改变燃料总截面影响中子输运过程，也给热工水力计算提供了总功率及各燃耗区的实际功率。实际反应堆的运行过程是中子输运、燃耗、热工水力学等多物理过程相互作用、相互影响的结果。

随着计算机硬件计算能力的不断提升，可视化先进建模、多物理、多尺度、精细化模拟获得了快速发展。美国制定并执行了CASL计划，构建了反应堆应用虚拟环境(VERA)，主要是充分利用现有软件建立松耦合环境，产生一个相比传统工具具有更高分辨率的耦合计算系统，对反应堆整个寿期运行历史进行预测。进一步的NEAMS计划是基于美国爱达荷国家实验室开发的面向对象多物理模拟(MOOSE)框架[1]，开发多物理耦合集成平台，实现系统-物理-热工-结构-材料及燃料的高精细耦合模拟软件开发。MOOSE框架采用了模块化的方法，允许科学家和工程师创建新的完全耦合的多物理应用程序。法国电力公司的数值模拟软件平台Salome能够用于数值模拟的通用预处理、后处理和工作流管理操作，能便捷地进行软件间的相互耦合运算和处理。近些年，国内核能领域也逐渐开始重视数值堆相关研究开发工作，并取得了不同程度的进展。中国工程物理研究院高性能数值模拟软件中心分析了数值反应堆多物理耦合关键技术[2]，提出采用中间件把许多共性的计算方法、领域并行编程的复杂问题，通过区域分解、数据分解、多级并行得到化解。

本文根据反应堆堆芯多物理耦合分析需求，研究多物理耦合算法，并完成耦合程序开发，实现中子物理、燃耗、热工子通道的多物理耦合计算。

1 多物理耦合求解方法

堆芯多物理耦合示意图如图1所示。耦合系统求解方法分为两种，如图2所示，一种是采用等效近似来进行物理量的解耦，然后通过参数传递和迭代计算来得到满意的结果；另外一种方法就是构建微分方程组直接联立求解，所有物理量同时更新，如近几年研究的JFNK算法[3-4]。

图1 堆芯多物理耦合示意图Fig.1 Schematic diagram of reactor core multi-physical coupling

图2 多物理耦合计算方法示意图Fig.2 Schematic diagram of multi-physical coupling calculation method

根据求解方法的不同，多物理系统设计中采用的耦合方式主要有以下3种：1) 松耦合，通常是在现有计算分析程序的基础上开发接口程序，程序间的数据传递基于双方的输入输出文件展开，这种耦合方式对程序的改动小，能充分利用已有程序，但程序间的相容性不好，容易存在收敛性的问题；2) 紧耦合，将耦合系统中描述不同物理过程的偏微分方程或方程组进行联立，在同一时间步长内，同时求解耦合系统中所有的非线性微分方程并进行变量更新，在计算过程中全部或部分使用内存进行数据的传递，这类程序涉及大量源代码的编写、重构与修改，需要解决相容性和通用性，是目前研发的热点；3) 耦合框架，搭建统一的I/O平台进行数据的处理和传递，各专业程序分别作为独立求解器内嵌于平台中，数据可通过文件或内存传递，是前两种耦合方式的集成，代表程序有美国爱达荷国家实验室开发的MOOSE系统[5-6]，ANASYS和COMSOL等大型商用软件也具备耦合框架平台的特点。

2 堆芯多物理耦合系统架构设计

含有多个偏微分方程的耦合系统直接求解的难度很大，即便是精细化求解单一物理过程的偏微分方程也存在诸多困难，因此本文工作首先采用基于参数传递的松耦合方法，基于单一物理过程的蒙特卡罗程序、特征线程序、燃耗程序、子通道程序，采用Master-Slave主从架构、Python语言编程，构建中子物理、燃耗、热工子通道的多物理耦合系统，如图3所示。在耦合架构设计中采用分层与封装策略，按照功能需求分为系统控制层、驱动层、数据层、执行层。同时基于网格类型和模拟参数，建立网格数据和参数的表示规范以及统一接口，节省存储资源、提高访问效率；基于不同空间离散形式的高效定位数据索引与差值算法，实现参数的高效转换与传递，提高耦合效率；基于具体耦合计算模型，建立耦合参数传递函数的映射关系，提高耦合参数传递的准确性；基于单级与多级时间步长控制方法，平衡模拟精度与计算时间。针对不同模拟模型，研究收敛控制方式的选择对程序稳定性和收敛性的影响。

图3 堆芯多物理耦合系统架构设计Fig.3 Architecture design of core multi-physical coupling system

3 耦合系统初步测试

3.1 输运-燃耗耦合计算测试

输运-燃耗耦合计算中，堆芯的中子通量密度、能谱分布及同位素核密度的参数传递与迭代计算原理如图4所示。蒙特卡罗程序进行临界计算，得到中子通量、单群反应截面等数据，传递给点燃耗程序；点燃耗程序进行燃耗计算，得到新的核素密度，传递给蒙特卡罗程序。通过数据的往返传递，从而处理燃耗计算的全过程。本次工作中子输运计算程序采用清华大学反应堆工程计算分析实验室(REAL)研发的三维粒子输运蒙特卡罗程序ANT-RMC[7-8]，该程序包含燃耗计算模块DEPTH[9]，采用程序内耦合模式开发，能够处理含1 500余种核素的精细燃耗链。程序首先通过临界计算(连续能量)模块得到中子通量、单群反应截面等数据，传递给点燃耗模块DEPTH，通过燃耗计算得到新的核素密度，传递给临界计算模块。

图4 输运-燃耗耦合计算流程Fig.4 Coupled calculation process of transport-burnup

ANT-RMC程序通过负载平衡和区域分解能实现大规模输运燃耗耦合计算，本文采用ANT-RMC程序对Hoogenboom-Martin基准题进行验证。Hoogenboom-Martin基准题[10]是一个含有燃耗计算的全堆芯基准计算模型，堆芯共包含241个燃料组件，每个燃料组件包含17×17个栅元，其中包括264个燃料棒，25个控制棒，燃料棒的高度为366 cm，沿轴向划分成24等份，堆芯结构与区域分解如图5所示。

图5 Hoogenboom-Martin基准题堆芯结构与区域分解Fig.5 Core structure and domain decomposition of Hoogenboom-Martin benchmark

对Hoogenboom-Martin基准题模型进行精细化pin-by-pin描述，采用96核并行规模，实现了百万燃耗区全堆三维精细燃耗计算。时间步长为250 d时的径向功率分布与有效增殖因数keff变化如图6所示。

图6 径向功率分布与keff的变化Fig.6 Change of radial power distribution and keff

3.2 核热耦合计算测试

在堆芯核热耦合计算中堆芯热工计算采用子通道分析方法。子通道分析方法将堆芯内复杂的流通面积划分子通道，通过子通道间质量守恒、能量守恒、动量守恒方程求解堆芯组件内的流场和温度场。本次计算采用的子通道计算程序为自主研发的子通道模拟程序CVR-PASA，CVR-PASA程序具备压水堆子通道计算功能[11]，同时也开发了面向六角形组件的液态金属快堆子通道计算功能，采用多层级索引实现从堆芯-组件-棒束-子通道多层级映射描述[12]，对钠冷快堆进行精细建模。快堆组件子通道映射示意图如图7所示。

图7 快堆组件子通道映射示意图Fig.7 Mapping diagram of fast reactor assembly subchannel

中子输运计算程序采用三维蒙特卡罗程序ANT-RMC，拥有中子光子输运、全能区在线截面处理、平衡氙、燃耗计算和衰变计算、全堆倒换料、MPI与OpenMP混合并行、接续计算等功能，具备在线截面处理功能，并能够在物理热工燃耗耦合中综合应用。ANT-RMC支持包括外耦合(由第三方脚本控制)、混合耦合(ANT-RMC系统调用控制热工计算)、内耦合(两者融为一个程序)在内的多种核热耦合方式，适用于不同的反应堆分析场景和不同计算平台的要求，具有分析实际反应堆多循环、多物理场(输运、燃耗、热工水力)耦合问题的能力。

核热耦合系统联立公式如图8所示，直接求解带热工反馈的输运方程是非常困难的，目前较常用的方式是Picard迭代方法：首先设定一个初始的温度和密度场的分布，然后开展输运计算，功率分布根据临界计算结果不断地予以更新，新的功率分布作为热工水力计算的条件，通过流动、传热及燃料棒导热计算，得到新的冷却剂温度、密度及燃料棒温度参数，存储在中子物理几何模块内的对应栅元下，迭代过程直到满足收敛条件为止。

图8 核热子通道耦合计算示意图Fig.8 Schematic diagram of nuclear thermal coupling calculation

核热耦合计算模型为商用示范快堆堆芯与组件[13]，采用ANT-RMC开展堆芯pin-by-pin精细化计算，燃料活性区轴向等分为20段，通过输运计算给出组件内部径向与轴向功率分布，然后采用CVR-PASA开展子通道计算，在稳态工况下，经过迭代反馈，燃料组件出口温度的计算结果与文献[14]中的组件子通道出口温度分布相近。耦合计算的温度与冷却剂密度分布如图9所示。通过核热耦合计算能够更为精确地描述轴向与径向功率分布对冷却剂的影响。核热耦合计算后续完善可用于钠冷快堆堆芯热工水力计算及流量优化分析。

图9 快堆组件的温度分布(a)和冷却剂密度分布(b)Fig.9 Temperature distribution (a) coolant density distribution (b) of fast reactor assembly

4 结论与展望

单一的物理计算难以全面模拟堆芯物理过程，多物理耦合计算成为研究多工况下堆芯行为演化规律的重要手段。基于先进耦合建模和大规模并行计算技术的数值反应堆已成为国内外领域前沿热点。数值反应堆不仅使上述需求成为可能，更为先进反应堆的设计优化、不同工况运行模拟优化、严重事故序列演示预测及燃料和材料研发提供一个经济高效的试验平台。本文构建了基于中子输运、燃耗、热工子通道的堆芯多物理耦合系统，通过输运-燃耗耦合计算测试和核热耦合计算测试初步验证了耦合系统的功能。

在后续工作中，将继续拓展多物理耦合系统，形成较为完备的全堆芯中子学、热工分析、燃料性能分析多物理耦合计算系统，同时开展瞬态工况下的紧耦合计算方法研究。通过燃料元件-组件-全堆芯的多尺度建模和物理-热工-燃料性能多物理场耦合计算，给出全堆芯精细化中子场分布、流场分布、温度场分布以及燃料组件性能在整个寿期内的精细化、量化分析结果，可以改变传统的以局部最大功率或最热组件为限值的保守性包络计算方法，为建立堆芯行为演化规律提供重要支撑。