吕炜枫 陈明亮 黄倩倩 宋 誉 洪 浩

（中广核工程有限公司核电安全监控技术与装备国家重点实验室 深圳 518172）

在压水堆核电站功率运行期间，来自堆芯的中子将会对堆外设备和墙体材料进行活化，从而产生活化产物源项，进而对停堆期间和退役期间核电站内工作人员产生辐射影响，而且随着核电站运行年限的增长，中子活化产生的长寿命核素（如60Co）的活度将逐渐增加，因此有必要建立堆外设备和墙体材料内活化产物源项的精细评估手段，从而为后端的屏蔽设计和废物处理提供输入。

美国橡树岭实验室基于矩阵指数法开发了求解燃耗方程的燃料组件中裂变产物核素和活化产物核素源项的计算程序 ORIGEN-S［1］。ORIGEN-S 程序基于点燃耗进行计算，其已广泛应用于国内外核电站设计。

欧洲原子能共同体和英国原子能管理局共同开发了欧洲压水堆活化产物源项计算程序FISPACT-2007［2］。FISPACT-2007基于欧拉指数法开发，可计算高能中子装置的材料活化以及压水堆和聚变堆中材料的活化情况。其计算配套的核数据库为EAF-2007。

华北电力大学韩文静等［3-5］分别基于欧拉指数法和切比雪夫有理近似法（Chebyshev Rational Approximation Method，CRAM）对燃耗方程进行求解，采用EAF-2007数据库开发了活化产物源项计算程序EuACT和ABURN。其中，EuACT程序为基于欧拉指数法开发，ABURN程序为基于CRAM方法开发。其基于ABURN程序计算了聚变堆第一壁活化例题和UO2燃料燃耗例题，并将ABURN程序的计算结果与欧洲活化产物源项计算程序FISPACT-2007进行对比。结果表明：ABURN程序可达到FISPACT-2007程序同等精度，并且由于采用了CRAM方法，程序在燃耗步设置方面具有高度的灵活性，初步验证了ABURN程序的可用性与准确性。

中国核动力研究设计院胡建军等［6］采用二维离散纵标输运程序DORT及ENDF/B-VI库计算压水堆压力容器内中子注量率分布，用自行研制的活化产物源项计算程序计算冷却剂16N和17N源项，并验证堆芯径向与轴向功率分布对计算结果的影响。其自行研制的活化源项计算程序名称为SCACT。

中国科学院核能安全技术研究所张彬航等［7］基于超级蒙特卡罗核模拟软件系统SuperMC，采用CRAM方法发展了基于深度优先搜索的活化链动态构建方法和大规模矩阵自适应降阶方法，进行了活化产物源项计算的初步研究。通过选用IAEA-ACB（InternationalAtomic Energy Agency-Activation Calculation Benchmark）国际基准例题及压水堆燃料包壳材料例题，初步验证了该活化产物源项计算方法的正确性，且测试结果表明发展的大规模矩阵自适应降阶方法能够有效地提高活化系数矩阵的求解效率。

综上所述，国内外各研究机构针对压水堆核电站活化产物源项计算程序进行了大量的研究工作，但因程序与工程设计脱节以及软件认证的问题，除了国外开发的相关程序，尚未有自主研发的活化产物源项计算程序在工程设计中大规模应用。在华龙一号“走出去”的大背景下，有必要研发可应用于工程实际的具有自主知识产权的压水堆活化产物源项计算程序。

1 核电站内中子活化产物源项计算模型

1.1 计算方程组的建立

核电站内受中子辐照的非裂变材料内放射性核素产生途径有以下两种：先驱核素衰变产生和先驱核素中子俘获产生。核电站内受中子辐照的非裂变材料内放射性核素消失途径有以下两种：放射性衰变消失以及与中子反应消失。可分别建立放射性核素i的产生率和消失率计算公式如下：

式中：Si(t)为t时刻放射性核素i的产生率，s-1；Pi(t)为t时刻放射性核素i的消失率，s-1；Φk为k群中子注量率，m-2·s-1；Ni(t)为t时刻中子辐照材料内放射性核素i的积存量；V为中子辐照材料体积，m3；nj为先驱核素 j的核子密度，m-3；σɑ，i，k为 k群中子对先驱核素i的微观吸收截面，m2；σx，j，k为k群中子对先驱核素j的微观x反应截面，x反应指（n，r）、（n，2n）、（n，p）、（n，ɑ）、（n，t）等中子反应，m2；λi为 i核素的衰变常数，s-1；Ki，j为衰变分支比，无量纲；Γx，i，j为先驱核素 j通过与中子的x反应产生核素i的产额，无量纲。

根据以上放射性核素产生率和消失率的计算公式，可建立放射性核素i积存量的微分方程如下：

求解以上微分方程可得中子辐照材料内放射性核素i的积存量。如需计算n种核素，则以上方程可写成矩阵的形式：

其中：

在以上方程中，非放射性核素也在n种核素之内，由此先驱核素中子俘获产生项可包含放射性核素和非放射性核素中子俘获活化的贡献。

对于式（3）中方程，其解为：

式中：X(t)是核素积存量矩阵；A是包含衰变和中子反应的转换矩阵（活化矩阵）。X(0)是初始核素积存量矩阵。

1.2 矩阵指数法求解

对式（8）中指数函数部分进行Taylor级数展开，有：

式中：I表示单位矩阵。矩阵指数法在理论上可以精确求解燃耗问题，但是在具体的例子中，求矩阵幂的精确值时会产生各种困难。在这种情形下，就需要借助Bateman方程或者使用Gauss-Seidel迭代法来解决。

燃耗矩阵A的规模通常在2000×2000量级，但是在实际的活化源项计算中，往往仅关注于几个主要的核素，极少出现一次性计算2000×2000的情况，考虑到核电厂实际情况，这个规模将大大缩小。

1.3 CRAM方法求解

CRAM方法的核心是在[0，∞)上对e-x进行有理逼近。定义πk，l表示所有如下形式实有理函数rk，l(x)的集合：

式中：pj是不高于j阶的多项式。根据逼近理论可知存在唯一的 r^k，l∈ πk，l满足：

对于给定的k和l，可通过Remez算法或Carathéodory-Fejér方法得到逼近表达式。由式（10）可知，有理函数r^k，k的极点{θ1，θ2，...，θk}是分立极点，因此可以通过部分分式展开得到：

式中：α0是函数在无穷远的极限；αi表示极点θi处的留数。式（12）也可以通过留数定理得到。另外，注意到αi和θi共轭出现，因此，式（12）可以简化成：

对于CRAM方法计算矩阵指数，主要的难点在于对给定的k值求出相应的有理函数的系数。原则上来说，p^k和q^k的系数都可以通过Remez型方法得到，但是计算过程需要相当有技巧的算法和非常高精度的计算过程。文献调研显示14阶系数的计算结果已经具有较高的精度［5，8］。基于此，CPAP程序采用16阶系数以进一步提高精度，采用的CRAM系数见表1。

2 活化产物源项计算程序CPAP的开发

基于以上中子活化产物源项计算模型，采用面向对象的编程思维以及当前主流的图形程序界面，开发了具有良好人机界面的各堆型压水堆核电站通用的中子活化源项计算程序CPAP。其中，程序图形化界面和计算模块均采用C++语言开发。CPAP程序的计算流程图见图1和图2。

CPAP程序有如下特点：

1）各堆型压水堆核电站在计算中子活化产物源项时，一般采用矩阵指数法或者CRAM方法，由此CPAP程序中提供了矩阵指数法和CRAM方法两种选项；

2）CPAP程序中采用的CRAM方法基于16阶系数，计算精度进一步提高；

3）CPAP程序采用的中子反应多群截面数据库基于ENDF/B-VII.0中的中子反应子库（NSUB=10）及EAF-2010中子活化反应数据库加工得出，涉及860种核素；

4）CPAP程序的截面数据库采用47群结构，能群划分参照BUGLE-96多群数据库的能群边界；加工多群截面数据库时，采用传统核截面加工处理软件NJOY中的相关运行模块，最终制成多群数据文件时采用压水堆堆芯外的中子能谱作为权重谱；截面数据库涉及的反应类型共11种，包括总反应、弹性散射反应、非弹性散射反应、辐射俘获反应、（n，2n）、（n，3n）、（n，α）、（n，p）、（n，t）、（n，2p）和裂变反应；

表116 阶CRAM系数Table 1 16th order coefficient of CRAM method

图1 CPAP程序计算流程图1Fig.1 Calculation flowchart 1 of CPAP code

5）具有良好的人机界面，根据用户需求预先嵌入大量常见材料数据库。

CPAP程序界面见图3。

3 活化产物源项计算程序CPAP的测试验证

对于活化产物源项计算程序CPAP的测试包括算法的测试以及数据库的测试两部分。以下分别说明。

3.1 算法的测试

选取了部分国际上通用的基准题以及自主建立的典型例题对CPAP程序的两种算法进行测试。算法测试时CPAP程序所使用的核素丰度、反应类型、截面数据和衰变常数均取自EAF-2007数据库，与FISPACT-2007程序一致。本文以两个典型的活化算例说明程序的算法测试结果。

3.1.1 AISI304活化计算测试

在本测试算例中，AISI304的初始质量为1 kg，其密度为8.04 g·cm-3。元素质量成分见表2。中子通量密度为 3.641×1015cm-2·s-1，中子辐照时间为1 a。中子能谱为100群的GAM-II格式能谱［9］。

AISI304活化例题的计算结果以及其与FISPACT-2007程序计算结果之间的对比如表3所示。由于计算所得核素较多，仅选取典型核素的结果展示。

图2 CPAP程序计算流程图2Fig.2 Calculation flowchart 2 of CPAP code

图3 CPAP程序界面展示Fig.3 Interface display of CPAP code

表2 AISI304的元素成分Table 2 Elemental content ofAISI304

表3中，Nf表示使用FISPACT-2007程序计算出的核素的原子数量；Nc表示使用CPAP程序中的CRAM求解器计算出的核素的原子数量；Nm表示使用CPAP程序中的矩阵指数求解器计算出的核素的 原 子 数 量 ；delta_c=（Nc-Nf）/Nf，delta_m=（Nm-Nf）/Nf。

3.1.2 混凝土活化计算测试

在本测试算例中，混凝土的初始质量为1 kg，其密度为2.35 g·cm-3。元素质量成分见表4。中子通量密度为3.641×1015cm-2·s-1，中子辐照时间为1 a。中子能谱为100群的GAM-II格式能谱［9］。

混凝土活化例题的计算结果以及其与FISPACT-2007程序计算结果之间的对比如表5所示。由于计算所得核素较多，仅选取典型核素的结果进行展示。

3.1.3 测试结论

经测试，对于测试算例，CPAP程序的两种算法在活化计算时大部分核素的计算结果与FISPACT-2007程序的计算结果之间的误差范围均在1%以下。但在H和He的计算时CPAP程序的计算结果与FISPACT-2007程序有较大误差，这主要是因为FISPACT-2007程序在计算H和He时实际上是合并计算了（n，p）和（n，α）反应产生的质子和α粒子，而CPAP程序不把（n，p）和（n，α）反应产生的质子和α粒子计算为H和He元素，所以FISPACT-2007程序计算结果高于CPAP程序的计算结果。

表3 算法测试时AISI304活化例题的计算结果Table 3 Calculation result ofAISI304 activation during algorithm test

表4 混凝土的元素成分Table 4 Elemental content of concrete

3.2 截面库和衰变数据库的测试

3.2.1 测试方法

CPAP程序的截面库和衰变数据库的测试采用与ORIGEN-S程序对比分析的方式。具体方法为首先使用CPAP程序的截面库和衰变数据库替换ORIGEN-S程序自带的截面库和衰变数据库，然后基于替换截面库和衰变数据库后的ORIGEN-S程序构建典型算例，并与原版ORIGEN-S程序的计算结果进行对比。

由于CPAP程序的截面数据库为47群的多群截面库，而ORIGEN-S程序最终读取和使用的中子反应截面为单群截面，由此需要对CPAP程序的截面数据库进行并群处理。并群所基于的中子能谱为使用蒙特卡罗输运计算程序对典型压水堆的堆芯进行输运计算后按照BUGLE-96库中能群边界统计得出的分为47群的中子能谱。

表5 算法测试时混凝土活化例题的计算结果Table 5 Calculation result of concrete activation during algorithm test

3.2.2 测试算例

对于CPAP程序的截面库和衰变数据库的测试算例见表6。其中混凝土和AISI304的元素成分分别见表2和表4。

选取各算例中典型核素的计算结果与ORIGEN-S程序原版计算结果进行对比，对比结果分别见表7和表8。由对比结果可知，部分核素存在小于10%的偏差，其可能原因为截面数据库替换时并群处理考虑的中子谱的差异以及截面加工所基于的评价库版本的差异。

表6 CPAP数据库的测试算例Table 6 List of test cases for database of CPAPcode

表7和表8中，Ns表示使用ORIGEN-S程序原数据库计算出的核素的原子密度，cm-3；Nc表示使用CPAP程序中的数据库替ORIGEN-S程序数据库后计算出的核素的原子数量，cm-3。delta_c=（Nc-Ns）/Ns。

表7 数据库测试时AISI304活化例题的计算结果Table 7 Calculation result ofAISI304 activation during database test

3.2.3 测试结论

由表7和表8中测试结果可知，对于所测试的算例，使用CPAP程序的截面库和衰变数据库所得绝大部分核素的计算结果与ORIGEN-S程序自带数据库所得计算结果之间的偏差在1%以内，少量核素存在不超过10%的偏差。考虑到测试时截面数据库替换过程中并群处理考虑的中子谱的差异以及截面加工所基于的评价库版本的差异，以上偏差可认为合理。

4 结语

本文基于堆外设备和墙体材料中活化产物的产生和迁移的机理建立了堆外设备和墙体材料中活化产物源项的计算模型，并编制了计算程序CPAP1.0版。基于典型材料的活化例题，采用国外同类型软件对比计算的方式对CPAP程序的算法和数据库分别进行了测试，测试结果显示：CPAP程序计算精度可满足工程设计的需要。CPAP程序具有人机界面友好以及求解器可选的优点，可广泛应用于压水堆核电站的设计、运行和退役阶段，解决堆外设备和墙体材料中活化源项计算的相关问题。

表8 数据库测试时混凝土活化例题的计算结果Table 8 Calculation result of concrete activation during database test