消息传递并行燃耗程序MCBMPI的栅元验证

杨万奎，刘耀光，马纪敏，杨 鑫，，王冠博，佘 顶

（1.中国工程物理研究院 核物理与化学研究所，四川 绵阳 621900；2.清华大学 工程物理系，北京 100084）

反应堆燃耗计算关系到反应堆核燃料管理、乏燃料后处理、燃耗信任制等，是反应堆领域的一重要分支。实现该目的的基本方法为中子输运计算和核素燃耗计算耦合。中子输运计算部分，基于确定论方法的一、二维几何处理能力的输运计算方式已不能满足反应堆燃耗计算的需求，随着计算机硬件的发展，基于蒙特卡罗方法的三维输运计算程序MCNP 得到了广泛应用［1］；核素的燃耗计算部分，基于矩阵指数法的单群点燃耗计算程序ORIGEN2占统治地位［2］。

因此，国内外均开展了基于MCNP 和ORIGEN 耦合系统的开发。一种为外耦合方式，将两个程序的计算结果相互提取作为自身的输入文件，如美国爱达荷国家实验室的MOCUP［3］、阿 拉 莫 斯 国 家 实 验 室 的Monteburns［4］、清 华 大 学 的MCBurn［5］等；另 一 种 直接在接口程序中加入核素转换、衰变计算功能，取代ORIGEN2程序计算，系统仅由MCNP和接口程序组成，减少了读取文件的时间，直接通过变量传递，如美国Sandia国家实验室开发的BURNCAL程 序［6］、清 华 大 学 的RMC［7］。本文采用第2种耦合方法，直接在程序中加入核素转换和衰变的计算功能，采用MPI并行方式（国内做过相关研究对比［8］），研制并行化蒙特卡罗燃耗程序MCBMPI，并以压水堆栅元燃耗基准题对程序进行验证。

1 程序原理

燃耗方程的求解采用3种方法结合，即截断泰勒展开的矩阵指数法、TTA 线性子链解析法和高斯-赛德尔迭代法。核素的分类参考ORIGEN2，分为锕系核素及其子体、裂变产物和活化产物3类，且初始数据库（衰变数据库、截面文件库）均采用该点燃耗程序配套库，后续燃耗步实时计算生成点截面库。

截断泰勒展开的矩阵指数法的适用对象为长寿命核素。指数函数的级数表示矩阵指数，并将级数展开至足够多项以达到指定的精度要求。其中，最关键的级数项的计算可用递归算法方便实现。TTA 线性子链解析法的适用对象为无长寿命先驱核的短寿命核素。此类核素导致分离本征值，不适于矩阵指数法求解，且在1个时间步长内，可认为此类核素已达平衡，可用解析解直接表示。高斯-赛德尔迭代法的适用对象为具有长寿命先驱核的短寿命核素，通过迭代法直接求解方程组。

该燃耗计算方法的另一特点是非空系数矩阵的存储技术。该方程组的系数矩阵维数可达1 700 多，且为大型稀疏矩阵。理论上可通过存储所有的方程组系数，求解方程组得到所有核素的浓度。但随核素的增多，系数矩阵的存储量将急剧增大（若以双精度浮点数8字节存储1个系数，则1 700×1 700的系数矩阵存储容量为22.049 MB），对系数索引极不利。因此采用系数索引方法，只存储非零元素。从数据库文件中读取各核素的半衰期、衰变分支比、截面、裂变产额等计算相关数据后，确定各核转换类型及其产物，用1个数组存放转换率及产物标识；通过查询该数组以获取每个子代核素Xi的父代核素Xj标识，将Xi的所有父代核素Xj的转换率存于一维浮点数组中；将父代核素标识存放于一维整型数组中；用计数器确定索引位置，实现非零系数矩阵的存储。

MCBMPI程序采用模块化思想，其中蒙特卡罗输运采用MCNP 程序的并行化版本MCNP5MPI，自编并行子程序MPIORI2实现核素截面替换、核素转变和衰变计算、核素信息的提取等功能，并行策略为对多燃耗区采用区域分解的MPI消息传递并行。程序整体框架如图1所示。

1.1 区域分解策略

区域分解通过自动读取输运计算中建模中的燃耗区几何栅元，并依次编号实现，之后对编号进行计算核心分配。一种方式为采用常用等间距分配（图2a），按计算核心1～N，接收燃耗区域1～TotalBurnupZone的计算任务。这种方式的区域分解方法实现起来很直接。但每个计算核心得到的结果在汇总时，材料序号的排列不是顺序排列的。虽可再排序，但增加了额外的计算任务。因此，需考虑一种可直接得到材料号顺序排列的区域分解方法：即让每个计算核心所计算的燃耗区域序号是连续的（图2b）。

图1 MCBMPI程序流程图Fig.1 Calculation flow chart of MCBMPI

图2 区域分解方式Fig.2 Type of zone decomposition

按此区域分解方法，程序实现过程中为包含所有的燃耗区，需用到向上取整函数ceiling（）。至此，完成了燃耗区的区域分解，将燃耗计算任务分配给各计算核心，通过获取相关截面信息后可并行调用燃耗方程求解器实现燃耗的并行计算。

1.2 燃耗并行存储策略

通常的燃耗计算并行存储方式为对各燃耗区分别建立对应的文件夹来存储，但在燃耗区较多时，文件夹的数量过多，管理起来较凌乱。出于以上考虑，在保存完整信息的前提下，根据MPI相关的文件I／O 并行处理方法，采取一新的并行存储策略，即根据计算核心的数量来建立相应的文件夹数量存储单个计算核心计算生成的结果。

燃耗并行计算中所需存储的结果主要包括燃耗计算的输出文件及输出文件中提取出的核素信息。这些核素信息可分为3类：下一次输运计算所需的核素成分；下一次核素转变和衰变计算的继承核素；下一次燃耗计算所需的替换截面核素。

通过采取上述的燃耗计算并行存储策略可对所有必需的结果实现并行化处理，在全部燃耗区都完成核素的转变与衰变计算后，可提供输运计算所需的全部参数。

1.3 单群截面的生成

通过MCNP5MPI输运计算可得到当前燃耗深度下系统的有效增殖因数、各燃耗区的中子通量密度以及可通过中子通量密度响应得到的量，如反应率、功率等［9］。另一方面，大部分锕系核素和裂变产物的各反应截面与具体装置对象的能谱相关，在进行1个新的燃耗计算前需对重要核素的截面进行修正或直接替换。一种有效的截面替换方法是通过输运计算实时地产生。

所需替换的截面类型根据不同的核素类型有所不同，但均只需考虑以下反应截面［2］：锕系核素的（n，γ）、（n，2n）、（n，3n）、（n，f）；裂变产物的（n，γ）、（n，2n）、（n，α）、（n，p）。各反应类型在MCNP5输运截面库中的对应关系列于表1。

表1 ORIGEN2反应类型与MCNP5截面库之间的对应关系Table 1 Relationship between ORIGEN2reaction type and MCNP5cross section library

本燃耗计算子程序中采用的核素截面为单群微观截面，相对于多群微观截面，更有利于MCNP5输运计算实时生成相应的替换截面，并在输运计算时的误差估计方面也放宽了要求，即为达到规定的误差要求所需的迭代次数减少了。单群微观截面的生成需保证反应率守恒，通常采用能谱权重得到，如式（1）所示。

输运计算中，只能分别统计得到式（1）中分子和分母。分母为全部能量区间的总通量，MCNP5中只需应用体通量计数即可得到。分子为微观反应率，MCNP5中可通过Fq乘子卡结合体通量计数得到。最终通过对MCNP5结果文件中通量和反应率的提取，得到替换核素的单群微观截面。

2 程序验证

2.1 OECD／NEA 压水堆栅元基准题

采用OECD／NEA 发布的压水堆栅元计算基准题［10］，对该程序进行验证。该基准题的目的在于对不同研究机构提供的针对该简化栅元模型的核素成分的计算结果进行比较。基准题结果发布时，包含了世界范围内16家研究机构的21套计算结果，主要是对系统反应性影响较显著的核素在寿期末的浓度进行对比，以验证燃耗计算程序及数据库。对比的核素包括12种 锕 系 核 素（234U、235U、236U、238U、237Np、238Pu、239Pu、240Pu、241Pu、242Pu、241Am、243Am）和16 种 裂 变 产 物（95Mo、99Tc、101Ru、103Rh、109Ag、133Cs、135Cs、143Nd、145Nd、147Sm、149Sm、150Sm、151Sm、152Sm、153Eu、155Gd）。

该燃料样品的化学分析结果来自于CE对14×14 的组件设计过程。为满足该基准题的目的，使用实际棒尺寸，但为使铀水比与实际的二维排列相当，对燃料芯体的尺寸进行了修改。燃料样品经历4 个完整的燃耗循环，各反应堆循环的辐照时间和衰变时间列于表2。该基准题包含3 种工况分别对应于参考燃料棒3个轴向位置处的燃料样品，各样品的总燃耗深度不同。

表2 栅元计算基准题的运行历史参数Table 2 Operation history data for benchmark problem pin-cell calculation

2.2 验证结果及分析

通过对3个工况的燃料样品应用MCBMPI程序进行计算，得到寿期末3个工况下的样品核素浓度，为评价MCBMPI程序的计算结果，分别给出了MCBMPI计算值与测量值和计算平均值两者的相对误差，将3个工况的相对误差集中示于图3。

根据各核素浓度的量级范围与核素浓度的测量精度，基准题中所采用的相对误差限为10%，即相对误差在±10%以内即认为理论与实验符合得很好。由图3a可知：239Pu、240Pu、150Sm、152Sm、153Eu 5 种 核 素 浓 度 的MCBMPI计算值与测量值的相对误差在10%附近；只有149Sm 的相对误差较大；其余核素的相对误差均在误差限以内，符合较好。

图3b 中，237Np、239Pu、240Pu、150Sm、152Sm 5种核素浓度的计算值与平均值的相对误差在10%附近；243Am、149Sm、155Gd 3种 核 素 浓 度 的相对误差较大，但其相对误差与基准题中给出的计算平均值的标准偏差相当，且所有计算值均在其他计算值的波动范围内；其余核素符合得很好，尤其是U 元素的各同位素和95Mo、99Tc、101Ru、103Rh、109Ag 等 核 素，相 对 误 差 在3%以内。

图3 浓度计算值与测量值（a）和计算平均值（b）的相对误差Fig.3 Relative error between calculation and experiment（a）and averaged calculation（b）concentrations

另外，浓度的误差较大的核素大部分微量，如149Sm 和155Gd，其浓度在给出的28种核素中最低，此处的浓度是1g UO2中目标核素的质量，3 个工况中这两种核素的浓度均在μg·g-1量级，243Am 的则为10-4μg·g-1量级，很小的偏差就能引起较大的相对误差；另一方面，截面及裂变产额的因素也不可忽视，这也是基准 题 中 给 出 的243Am、149Sm、155Gd 3 种 核 素 的计算平均值存在较大标准偏差的原因［10］。导致截面差异的原因除基准库不同外，燃耗步长是重要因素。因为这3 种核素的热中子吸收截面较大，对中子通量和能谱的影响较大，反过来又影响自身核素的燃耗计算，故燃耗步长对其影响较大。程序对该基准题兼顾计算时间与计算精度，设置的燃耗步数为4，计算结果也在基准题给出的计算值范围内。另外，程序中计算参数的选取亦可导致计算结果的偏差，需根据试算结果进行优化，如核素选取时的截断份额（核素继承截断份额设置为-25 次方，核素吸收截断份额设置为-6次方）、核素转换计算时内步长的确定（内步长设置为10d）。

3 加速比与并行效率

以工况A 为例说明采用并行计算在计算时间方面的优势。计算硬件条件为Intel®CoreTM2 Quad CPU Q8300 ＠ 2.50 GHz，3.25G 内存，计算耗时及加速情况列于表3。其中，加速比为多核计算耗时与单核计算耗时之比，并行效率为加速比与计算核数之比。

表3 MCBMPI的并行耗时及计算加速Table 3 Collapsed time and parallel speedup of program MCBMPI

由于该程序系统是由两个部分组成的，整体加速效果是两部分的综合作用，因此需分别讨论其加速情况。对MCNP5MPI而言，可从单次MCNP5MPI的计算耗时得到该程序对于该基准题的并行加速情况，其结果列于表4。

表4 MCNP5MPI的并行耗时及计算加速Table 4 Collapsed time and parallel speedup of program MCNP5MPI

从表3、4的加速比及并行效率可知，该并行燃耗计算程序系统的并行效率低于MCNP5MPI的并行效率。这是由于程序系统融入了输运计算和燃耗计算两部分，MCNP5MPI只包含输运计算，而燃耗计算采用区域分解并行。燃耗区小于计算核数时，并行效果不能完全体现出来，但若应用到多燃耗区的燃耗计算则能完全体现并行加速，这将在后续的实际应用中得以体现。

4 小结

本文开发的并行蒙特卡罗燃耗计算程序MCBMPI与压水堆栅元燃耗基准题的对比结果表明：各核素的计算值是可信的，且并行计算结果与串行结果完全一致，该程序加速了蒙特卡罗燃耗计算，可用于多燃耗区的燃耗计算。尤其是U 元素的各同位素的计算值与实验值与基准计算值均符合很好，表明该程序在燃料管理与换料分析方面都具有应用前景。

［1］ 邓力，李刚.粒子输运蒙特卡罗模拟现状概述［J］.计算物理，2010，27（6）：791-798.DENG Li，LI Gang.A summarization on Monte Carlo simulation［J］.Chinese Journal of Computational Physics，2010，27（6）：791-798（in Chinese）.

［2］ CROFF A G.A user’s manual for ORIGEN2 computer code，ORNL／TM-7175［R］.US：Oak Ridge National Laboratory，1980.

［3］ MOORE R L，SCHNITZLER B G，WEMPLE C A，et al.MOCUP：MCNP-ORIGEN2coupled utility program，INEL-95／0523US［R］.US：Idaho National Engineering Laboratory，1995.

［4］ TRELLUE H R.Development of Monteburns：A code that links MCNP and ORIGEN2in an automated fashion for burnup calculations［R］.US：Los Alamos National Laboratory，1999.

［5］ 余纲林，王侃，王煜宏.MCBurn：MCNP和ORIGEN耦合程序系统［J］.原子能科学技术，2003，37（3）：250-254.YU Ganglin，WANG Kan，WANG Yuhong.MCBurn：A coupling package of program MCNP and ORIGEN［J］.Atomic Energy Science and Technology，2003，37（3）：250-254（in Chinese）.

［6］ PARMA E J.BURNCAL：A nuclear reactor burnup code using MCNP tallies，SAND2002-3868［R］. US：Sandia National Laboratory，2002.

［7］ 佘顶，王侃，余纲林.堆用蒙卡程序燃耗计算功能开发［J］.核动力工程，2012，33（3）：1-5.SHE Ding，WANG Kan，YU Ganglin.Development of burnup calculation function in reactor Monte Carlo code RMC［J］.Nuclear Power Engineering，2012，33（3）：1-5（in Chinese）.

［8］ 邓力，张文勇，徐涵，等.蒙特卡罗程序MCNP-Ⅱ与MCNP-5并行效率比较［J］.计算机工程与科学，2009，31（A1）：185-187.DENG Li，ZHANG Wenyong，XU Han，et al.The comparison of parallel efficiency between Monte Carlo code MCNP-Ⅱand MCNP-5［J］.Computer Engineering ＆ Science，2009，31（A1）：185-187（in Chinese）.

［9］ X-5 Monte Carlo Team. MCNP：A general Monte Carlo N-particle transport code，Version 5，LA-UR-03-1987［R］.US：Los Alamos National Laboratory，2003.

［10］de HART M D，BRADY M C，PARKS C V.OECD／NEA burnup credit calculational criticality benchmark phase Ⅰ-B results，ORNL-6901［R］.US：Oak Ridge National Laboratory，1996.