钠冷快堆先进建模与高性能粒子输运数值模拟进展

周培德，薛小刚，吴明宇,*，王事喜，吴宗芸，霍兴凯，高鑫钊，张 强，付元光，周生诚，刘 晨，马建明

(1.中国原子能科学研究院，北京 102413；2.中国工程物理研究院 高性能数值模拟软件中心，北京 100088)

钠冷快堆是第四代先进核能系统的候选堆型之一。快堆和压水堆在系统布置、冷却剂、燃料组件形式、中子能谱、运行工况等方面存在差别。钠冷快堆堆芯组件采用更为紧凑的燃料元件棒三角形排列，组件外部为六角形外套管形式。燃料元件棒通过绕丝定位，按三角形栅格布置，燃料元件的精细功率分布是燃耗分析、热工流体设计等重要输入，且相互耦合。燃料管理、换料优化、次锕系核素嬗变、堆芯燃料增殖，以及乏燃料贮存、运输和后处理等计算分析对燃耗及同位素成分提出越来越高的精度要求，需基于精细化、高精度的三维输运-燃耗耦合计算。池式结构快堆堆本体几何尺寸大、结构复杂，对中子学计算建模与输运计算也提出挑战，全范围三维精细建模才能给出可靠、准确的计算结果。

基于先进建模和大规模并行计算技术的数值模拟是国内外反应堆技术研究前沿热点。美国制定并执行了CASL计划[1-2]，借助高性能数值模拟，解决轻水堆核电站运行功率提升、寿期延长、燃耗加深等问题。后续开展的NEAMS计划[3]则在开发面向对象的多物理耦合集成模拟平台基础上，正在从燃料产品线和反应堆系统级产品线开展物理-热工-结构-材料及燃料的高精细耦合模拟平台开发。欧洲先进数值模拟研发从反应堆安全性出发，自2004年起分NURESIM、NURESP、NURNEXT[4-5]3个阶段实施，通过集成最先进的物理、热工水力、燃料程序构建具有多尺度、多物理场的仿真平台。俄罗斯在突破计划“Breakthrough”[6]中提出了针对钠冷、铅冷、铅铋冷的快堆谱系的新一代设计与安全分析软件开发计划，突出多尺度建模、精细化模拟、集成计算与减小不确定度。

近几年，我国在多物理、多尺度、先进数值模拟方面也取得进展。吴宜灿等[7]开发了多物理耦合分析自动建模软件SuperMC；邓力等[8]、刘晓晶等[9]均开展了数值反应堆多物理耦合关键技术分析；杨文等[10]开展数值反应堆原型系统开发及示范应用研究，开发基于超算的精细化模拟软件；余红星等[11]对数字反应堆技术进行了调研、分析和论证；北京应用物理与计算数学研究联合上海核工程研究设计院、西安交通大学共同开展CAP1400数值反应堆关键技术研究，突破堆芯高保真模拟与不确定度分析技术；吴宏春等[12]对液态金属快堆专业软件进行了总结并展望了发展趋势；曹良志等[13]研发了高保真数值反应堆程序NECP-X，开展核设计不确定度分析。

反应堆数值计算发展趋势正从单一物理过程向多物理系统集成、从简化的等效模型向精细化先进建模与高性能全范围模拟、从保守的安全裕量包络性计算向降低不确定度与高可靠性方向等发展。

本文综述国内外反应堆高性能数值模拟进展，介绍钠冷快堆高性能粒子输运系统的可视化先进建模技术和计算功能，以及在快堆堆芯精细化功率分布、输运-燃耗耦合计算、堆本体屏蔽计算、堆外探测器深穿透计算、堆芯核-热耦合计算分析等的应用进展。

1 快堆先进建模与数值模拟系统

基于大型池式快堆物理、屏蔽三维精细化计算需求，中国原子能科学研究院与中国工程物理研究院高性能数值模拟软件中心合作，开发了一套快堆专用(钠冷快堆、铅铋冷快堆)、国内自主的高性能粒子输运计算系统JPTS-CFR，如图1所示，包括前处理引擎、支持结构网格[14]、组合几何体[15]、并行中间件框架以及后处理可视化展示，可达百万亿次计算规模的数值模拟。系统包括高性能蒙特卡罗程序JMCT[16-17]，快堆高性能三维SN程序JSNT-CFR[18]，并在此基础上扩展燃耗耦合计算功能和快堆热工子通道耦合计算功能，开展堆芯与堆本体计算。

粒子输运计算广泛应用于反应堆物理、屏蔽计算、临界安全、辐射输运计算等领域，要开展高分辨率数值模拟工作，基于CAD模型的可视建模与网格自动剖分成为一种理想选择。JPTS-CFR系统支持CAD建模、网格剖分前处理、网格-材料自动映射、蒙特卡罗并行计算、三维SN并行计算、领域编程框架、可视化后处理功能。前处理功能具备可视化交互、重复几何阵列操作、自动网格划分、几何体相关性自动转换的功能，提供层次化几何树的输出文件，保证计算效率，降低建模出错率，无缝对接蒙特卡罗与三维SN求解算法。通过与高性能计算技术结合，可完成大型复杂结构体的精细化建模与数值模拟计算。

图1 高性能粒子输运计算系统JPTS-CFRFig.1 High-performance particle transport computing system JPTS-CFR

池式结构快堆堆本体主要由堆芯、堆内屏蔽、堆顶固定防护平台、一回路主泵、中间热交换器、独立热交换器、主容器、生物屏蔽等构成。对于1 500～2 500 MW热功率池式结构快堆，其反应堆主容器直径可达约15～18 m，主容器底部到堆顶固定防护平台的高度可达17～20 m，堆本体结构的复杂程度和几何规模对建模带来较大难度。

采用JPTS-CFR系统对池式结构钠冷快堆堆本体开展精细建模和高性能数值模拟计算流程如图2所示，以实体结构的三维堆本体模型作为前处理输入。前处理建模采用无穷嵌套父子映射关系的几何树结构，支持栅元快速检索、隐式重复建模、高效重复根节点下的所有实体，保证相同建模复杂度的情况下提高可视化效率，前处理建模几何树结构如图3所示。

图2 JPTS-CFR计算流程示意图Fig.2 Schematic diagram of JPTS-CFR calculation flow

图3 前处理建模几何树结构Fig.3 Geometric tree structure of preprocessing modeling

示例中堆本体全范围(含堆坑及其生物屏蔽)构建了16 000余个不规则几何实体，形成大规模实体阵列，通过大量布尔运算和实体的任意旋转平移，自动转换生成了近6万行的几何映射模型文本文件，实现了手动建模无法达到的复杂规模。示例最后生成的可视化模型如图4所示。

图4 池式结构快堆堆本体粒子输运计算可视化模型Fig.4 Calculation model of pool type sodium-cooled fast reactor block

采用JPTS-CFR系统对快堆堆芯结构进行全堆芯范围精细化建模。快堆燃料组件由六角形外套管和燃料元件组成，堆芯区域燃料元件数目达几万根。根据几何结构和材料成分定制燃料元件-组件-堆芯多层级实体模型，通过重复几何阵列操作快速实现堆芯pin-by-pin精细化建模，并提供层次化几何树的输出文件。

2 计算应用

2.1 快堆堆本体屏蔽计算

池式结构快堆屏蔽计算面临的问题主要有两方面：池式结构快堆结构复杂、尺寸大、粒子深穿透问题。快堆中子通量高、中子能谱较硬，使得靠近堆芯的部件辐照损伤严重，堆芯围桶、堆芯支承等不可更换构件所能承受的辐照损伤是决定反应堆寿命的关键因素之一；池式快堆二回路钠需要流经布置在主容器内的中间热交换器，为降低中间热交换器中二次钠的活化，保证反应堆在运行过程中二回路系统的可接近性，是池式结构快堆屏蔽设计的另一个关键问题[19-20]。

SN方法是解决深穿透屏蔽问题国际上公认的较好方法，大规模并行SN计算程序依托高性能计算，在深穿透和获取物理量的全局分布上具有明显的优势。JSNT-CFR采用群平衡、空间再平衡、上散射再平衡加速等内迭代与外迭代加速算法、区域分解技术实现SN离散方程的大规模并行计算，提高精细化能力和计算效率。根据堆本体径向、轴向、中子引出通道的材料布置，通过共振自屏修正和一维并群计算生成多群截面库[21-22]，然后根据堆本体网格剖分模型开展屏蔽计算。

本文示例中JSNT-CFR与DOT3.5计算均采用34群中子和22群γ光子的耦合截面库[23]。JSNT-CFR计算示例为堆本体横切面1/2模型，整个模型为长(30 m)×宽(15 m)×高(20 m)的空间区域，按照小于不同材料区域的平均中子吸收自由程进行网格划分为454×233×337，约3 500万结构网格，离散方向采取S16，各项异性取P5阶展开[24]。在并行计算服务器上，首次实现了大尺度堆本体、全范围三维屏蔽计算，最大并行规模为1000核，时间约为7 h。池式结构快堆堆内设备布置复杂且有一定的非对称性，屏蔽设计中需计算堆芯组件、冷却剂、堆芯支撑、换料机、控制棒动导管、含硼石墨屏蔽柱、主泵、中间热交换器等设备的中子辐照损伤或中子活化水平。如采用DOT3.5二维计算程序，需选取主泵、或中间热交换器、或设备间隙、或中子引出通道等不同角度方位的纵剖面进行二维等效计算。在同样收敛准则情况下，采用JSNT-CFR程序，可直接获得堆本体范围三维精细化中子注量率、光子注量率、剂量率等分布，计算结果通过Teravap可视化后处理，全范围展示出堆本体中子注量率三维分布云图，相较于以往的DOT3.5二维计算，显著提升了展示效果，中子引出通道方位二维纵剖面DOT3.5计算结果展示与JSNT-CFR三维计算展示效果对比如图5所示。

图5 堆本体中子注量率分布展示效果对比Fig.5 Display comparison of neutron fluence rate distribution in reactor block

为验证三维计算结果的正确性，在堆内不可更换构件的纵剖面，采用DOT3.5进行局部二维等效计算，与三维计算结果较为接近，最大积分中子注量对比列于表1。

表1 快堆积分中子注量对比Table 1 Integrated neutron flux comparison in fast reactor

通过CAD前处理建模、网格剖分、高性能计算、可视化展示能够更为全面和精细地给出大型池式钠冷快堆堆本体三维屏蔽计算结果，在局部关键位置与已有DOT3.5程序进行对比，验证方法合理可行，为进一步降低设计保守裕度和局部结构优化提供数据支撑。

2.2 堆外探测器中子注量率深穿透计算

池式结构快堆中子注量率监测系统有堆内和堆外两套系统组成。其中，堆内测量系统使用高温电离室探测器，位于钠池中的探测孔道内，用于装料、换料工况的中子注量率监测；堆外监测系统的中子探测器位于反应堆保护容器外的石墨腔内(图6)，测量反应堆启动到满功率运行期间的中子注量率。

图6 堆芯中子引出方向示意图Fig.6 Schematic diagram of core neutron extraction direction

堆外探测器区域等效热中子注量率[25]计算公式如下：

(1)

传统方式采用二维等效保守计算，通过中子引出通道到达堆外电离室石墨空腔的局部近似处理，用DOT3.5进行二维计算，给出堆外石墨腔中子注量率的保守下限值。基于JPTS-CFR开展堆本体三维精细化建模，分别用蒙特卡罗方法和确定论SN方法开展计算，在蒙特卡罗程序JMCT计算中采用了基于源偏倚和权窗方法的减方差技巧，尝试减小深穿透问题带来的计算偏差。

蒙特卡罗方法基于中子输运积分方程构建概率模型，然后进行大量的随机试验并确定各种统计量的估计值，对于大的几何形状系统和小概率问题，例如屏蔽计算中的深穿透问题，会导致统计方差较大，计算结果较真值偏低，几何区域越大，这种偏低现象越明显。裴鹿成[26-27]和王瑞宏等[28]在文献中对上述问题也有所描述，并提出了改进方法。关于深穿透问题国内许多学者也都开展了大量研究工作，并取得了一定的进展。张寅等[29]、李锐等[30]开展屏蔽深穿透问题研究，潘清泉等[31]开展能量偏倚减方差方法研究，崔甜甜等[32]开展多种减方差技术在屏蔽计算中的组合应用，张斌等[33]、郑征等[34]分别提出蒙特卡罗程序与三维SN程序耦合计算以解决复杂结构的深穿透问题。

JPTS-CFR系统在在300万亿次计算中心完成了1 000 MWe快堆堆本体建模、千万网格、千亿粒子数、固定源、4800核并行规模的中子输运计算，中子在堆本体内中子注量率分布云图三维展示如图7所示。

图7 堆本体中子注量率计算结果示意图Fig.7 Schematic diagram of neutron fluence rate of reactor block

JMCT、JMCT减方差、JSNT-CFR的计算结果与DOT3.5二维等效计算结果对比列于表2。

表2 中子引出通道方向中子注量率对比Table 2 Comparison of neutron fluence rate in neutron extraction channel direction

从结果能够看出采用JSNT-CFR计算结果与DOT3.5的计算结果较为接近。中子引出通道方向上的围桶，堆内探测器位置的中子注量率JMCT、JSNT-CFR和DOT3.5均较为接近。随着输运距离的增加，在挤钠器和堆外探测器位置处JSNT-CFR与DOT3.5计算结果仍较为接近，JMCT在未采用减方差技巧的情况下计算结果明显偏低，堆外探测器中子注量率只有105cm-2·s-1的量级，相对确定论计算结果偏低约3个量级。采用权窗和源偏倚计算的情况下计算结果明显改善，略微偏低，快堆复杂结构深穿透计算仍需开展进一步研究。

JPTS-CFR系统通过先进建模与并行计算技术，提高了精细化描述能力，扩展了计算维度，对快堆中子探测器设计提供了支撑。

2.3 堆芯输运-燃耗耦合计算

快堆的燃料循环特性、堆芯装载方案等通常通过输运-燃耗耦合计算来进行模拟分析。关于输运-燃耗耦合，国内外均开展了大量研究工作，如TRIPOLI[35]、Serpent[36]、MCNPX[37]程序系统，国内薛小刚等[38]开发的MCNP和ORIGEN2耦合程序，清华大学REAL团队开发的堆用蒙特卡罗程序RMC[39]，成昱廷等[40]开发的CRAM方法燃耗程序及耦合系统，李静惊等[41]开发的基于确定论的三维燃耗程序等。STEP1.0是自主研制的点燃耗计算程序[42]，采用线性子链(TTA)方法进行燃耗方程求解。

燃耗计算本质上属于一种细粒度、高频调用的计算模式，存在天然的可并行性，具备大规模并行计算的条件。以1 000 MWe钠冷快堆典型设计为示例，堆芯活性区含有541燃料盒组件，每盒组件采用169棒束，含有9万多根燃料棒，若考虑轴向功率差异，需要将燃料棒轴向分成多层，输运-燃耗耦合pin-by-pin模拟，燃耗区会达到数十万，若考虑每个燃耗区内的能谱、反应率计数，则计数量会突破千万。JMCT程序通过大规模并行技术实现了百万几何体和计数区域的pin-by-pin输运模拟。在此基础上，通过对STEP1.0程序进行并行重构工作，将STEP1.0与JMCT进行内耦合，形成一套新的输运-燃耗耦合计算系统JMCT+STEP2.0，具有输运-燃耗耦合大规模并行模拟能力，能够运行在百万亿次大型服务器上，提升数值模拟的分辨率。

示例将每根棒轴向分成5层，燃耗区超过40余万。采用6个燃耗步长迭代，模拟160 d的燃耗过程，6000核并行计算规模，计算时间为27 h。燃耗总反应性变化与快堆核设计程序NAS[43-44]分析结果较为接近(表3)。

表3 燃耗反应性效应计算对比Table 3 Comparison of burnup reactivity effect

示例表明进行全堆芯范围精细化输运-燃耗耦合计算时间可接受，且计算结果的可靠性有保证，能够分析燃料中的超铀元素成分的变化对堆芯局部特性的影响，为快堆堆芯方案设计提供指导。

2.4 快堆核-热子通道耦合计算

钠冷快堆热堆芯热工设计需要根据物理计算的功率分布，考虑堆芯总流量和燃料组件、转换区组件以及堆芯其他组件的流量分配，并根据流量分配结果和堆芯压降确定燃料组件管脚开孔。通过热工计算给出各流量区最热燃料组件的钠出口温度、燃料棒包壳中壁热点温度和最高燃料温度，以及全堆芯燃料平均温度，为安全评价提供依据。

子通道分析程序能够兼顾准确性以及高效性，在堆芯热工设计计算中，通常给出每个流量最热燃料组件及最热燃料棒的子通道计算分析，按照包络性设计满足热工设计限值，往往存在进一步优化的空间。为了更为精确地分析核设计边界条件，反应堆核热耦合计算一直是研究的热点，蒙卡程序MC21[45]、TRIPOLI[46]、Serpent[47]均开发了核热耦合功能，国内有唐琪等[48]、叶辛欧文等[49]、刘鹏等[50]、卢旭等[51]先后开展了核热子通道耦合方法研究及程序开发。子通道分析方法是将流体的质量、能量、动量守恒方程应用于燃料棒束通道的流动区域，通过对质量守恒、能量守恒、动量守恒方程的离散得到形式上为Ax=b的线性代数方程组，如式(2)所示。

Ax=b,A∈Rn×n,x,b∈Rn

(2)

式中：A为与子通道规模相关的n×n矩阵；x为需要求解的n维列向量；b为边界条件n维列向量；Rn为向量空间。

通过求解方程组以得到冷却剂的压力、焓、流量的分布。堆芯中某个子通道仅与附近的几个子通道通过间隙相连，因此在式(2)中，矩阵A为稀疏矩阵。并且矩阵的维数与子通道的数目相同。以采用六边形组件结构的中国实验快堆为例，全堆芯内有712盒组件，包括燃料组件、不锈钢组件、控制棒组件多种类型，燃料组件内包含61根燃料棒[52]，全堆芯模拟往往需要计算上万子通道。对于1个大功率快堆堆芯，对其采用栅元级别的精细化子通道划分时，总的子通道数目会超过十万规模，此时，矩阵A将成为一个大型稀疏矩阵。针对上述需求，中国原子能科学研究院自主研发了通用子通道分析程序VERSA[53]，采用非结构化网格稀疏矩阵存储格式和Krylov子空间迭代法，具备快堆六角形组件全堆精细化子通道计算功能，计算效率满足快堆全堆芯大规模子通道分析。

在快堆全堆核-热耦合方案中，首先通过对堆芯物理计算结果中的功率分布、中子注量率信息的自动识别、提取、映射与转换，然后根据堆芯-组件-棒束-子通道多层级网格映射关系，为堆芯中所有的组件、棒束、子通道和间隙建立一个唯一的标识对象，标识包括编号和类型，供后续模拟时计算机进行逐个遍历。最后，根据标识对象的空间关联关系和位置信息，将物理计算的三维pin-by-pin功率分布信息进行映射和传递，完成全堆物理计算到子通道计算的单向耦合计算，流程如图8所示。

图8 核-热子通道耦合流程示意图Fig.8 Schematic diagram of neutronic-thermal sub-channel coupling process

在1 000 MWe级钠冷快堆堆芯设计中，燃料组件的子通道划分如图9所示。整个组件划分为342个子通道。在轴向上等距划分20个节点。整个反应堆活性区燃料组件为541盒，全堆芯总共包括18万个子通道，VERSA全堆子通道计算约为46 min，计算结果得到的反应堆堆芯出口温度分布如图10所示。

图9 全堆芯燃料组件子通道划分示意图Fig.9 Schematic diagram of sub-channel division of full core fuel assembly

快堆组件及燃料棒束采用三角形排列，组件盒之间存在盒间流，快堆全堆芯高精细的子通道模拟时，棒束/子通道的数量、关联关系的复杂性相对压水堆均成倍增加。通过快堆网格映射，初步实现了快堆全堆芯核-热子通道耦合分析计算。组件出口温度和温升计算结果与SUPERENERGY程序[54-55]单组件子通道模型计算分布相近。对比结果列于表4，表明计算结果可靠，采用的技术路线是可行的。通过全堆芯范围的核-热子通道耦合计算，可以细致分析功率分布导致的出口温度差异，为进一步优化堆芯流量分配、降低设计保守裕量提供可能。

表4 堆芯燃料组件热工计算结果Table 4 Thermal calculation result of fuel assembly

3 总结

基于先进建模、高性能计算、多物理多尺度耦合的数字化工程设计是反应堆专业与工业数值模拟方法、数字信息技术及计算机硬件资源发展到一定阶段的必然趋势，是典型的反应堆专业、计算机科学、软件工程多学科交叉的新型技术。本文介绍了高性能计算在大型池式快堆中堆本体屏蔽、堆外探测器、堆芯输运-燃耗耦合、堆芯核-热耦合方面的初步应用，通过建模和可视化功能，提高了模型的精细化描述水平和计算结果的展示度。

在后续工作中，将快堆燃料元件性能分析程序与物理计算程序、热工计算程序进行耦合，获取全堆芯燃料元件功率、燃耗、冷却剂温度随整个寿期的变化，实现对整个寿期的燃料元件热学、力学、辐照现象的性能分析，预测燃料元件寿命，形成物理-屏蔽-燃耗-热工-燃料性能快堆多物理耦合计算系统。利用高性能计算与先进建模技术，提高大型池式快堆堆本体三维粒子输运屏蔽计算、稳态快堆堆芯中子学计算、全堆芯输运-燃耗耦合计算、核-热子通道耦合计算、核-热-燃料性能耦合计算分析能力，既能够完成单一物理过程的精细化计算，也能够完整的模拟示范快堆中子物理、热工水力及燃耗的多物理耦合过程，对于提高计算精度、保证设计质量，具有重要意义。