基于MCNP 的反应堆建模方法

王武，夏虹，李伟，巴少华

1.哈尔滨工程大学 核安全与先进核能技术重点实验室，黑龙江 哈尔滨 150001

2.哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001

由于蒙特卡罗(Monte Carlo,MC)方法相比于确定性方法(如Boltzmann 输运方程的数值解)，能更准确地表示几何结构和核数据，例如：确定性方法由于使用数值求解技术，因此需要合理地简化研究对象的几何形状，并且使用多群截面数据来近似连续能量截面数据，而MC 方法既可以处理简单的几何和多群截面数据，也可以处理复杂的几何形状和连续的截面数据。因此对于许多研究对象而言，由于其内部结构十分复杂且材料不均匀，MC 是一种更好的模拟方法。

MCNP 作为核能领域的MC 计算程序，发展最为成熟、最具代表性。近年来，国内外学者利用MCNP 在反应堆有效增殖因子[1-3]、反应性[4-5]、中子通量和功率分布[3,6-8]以及中子能谱[9-10]计算等方面做了诸多研究。虽然前人利用MCNP 模拟反应堆临界做了一系列工作，但是，在建模过程中，多数学者只关心反应堆堆芯部件的几何结构、材料以及数据库选取，没有关注影响模型准确性的其他因素，如不可分辨共振区概率表、自由气体热散射的温度修正以及热中子S(α,β) 模型等。因此，本文以某压水反应堆为研究对象，综合分析了上述因素对MCNP 模拟的反应堆有效增殖因子、堆芯能谱以及功率分布的影响，从而提出一种合理的MCNP 反应堆建模方法。

1 反应堆模型描述

1.1 堆芯几何描述

本文以某压水反应堆堆芯为研究对象，选取0 个有效满功率运行天(0 effective full power operating days，0EFPD)燃耗时刻堆芯热态满功率运行工况进行建模。堆芯组件采用7×7 排布方式，组件按控制棒和可燃毒物棒的装载方式不同可分为7 类，分别用A、B、C、D、E、F、G 表示。每个组件中都留有16 根导向管，用于容纳控制棒和可燃毒物棒，每根导向管栅格占据4 个燃料棒栅格的位置。A、B、C、D、E 为含有控制棒的组件，其中A、B、C 为中心棒组，各插入12 根控制棒和4 根可燃毒物棒，D、E 为外围棒组，各插入16 根控制棒；F、G 为不含控制棒的组件，其中F 类组件中插入16 根可燃毒物棒，G 类组件中插入12 根可燃毒物棒。各类组件在堆芯中呈中心对称布置如图1 所示。堆芯内部结构均按照实际尺寸进行精细化建模，燃料棒气隙、可燃毒物棒气隙和组件间水隙等微小结构也考虑在内。

图1 组件在堆芯中的排布

使用MCNP 中的lattice(lat)卡、fill 卡以及universe(u)卡对燃料组件几何进行分级描述。以E 类组件为例，组件结构分三级描述，第一级为组件，第二级为填充组件的栅格，第三级为填充栅格的燃料棒、控制棒以及水。首先，描述燃料棒以及周围的水并分配世界体编号u=2。然后，描述控制棒，由于一个控制棒栅格占据4 个燃料棒栅格的位置，而使用lat 卡和fill 卡进行矩阵式填充时栅格的大小必须一致，所以需要将控制棒栅格划分成4 个大小相同的栅格并分别分配世界体编号u=3、u=4、u=5 和u=6，每个栅格中都包含有1/4 控制棒、1/4导向管以及除控制棒和导向管以外的水。接着，定义参考栅格，并使用lat 卡和fill 卡将各个世界体填充的栅格按矩阵形式进行重复，其中参考栅格作为矩阵的元素。为所有栅格分配世界体编号u=1，最后使用fill=1 将栅格填入组件中。图2 为MCNP5 构建的E 类组件模型的径向截面图，其中不同的颜色代表不同的材料。其他E 类组件只需要使用like…but…卡进行复制和坐标变换。除燃料组件以外其他堆芯结构(如围板、吊篮、压力容器等)的几何描述相对比较简单，因此不再详述。图3 为MCNP 构建的堆芯整体模型的径向截面图。

图2 E 类组件径向截面

图3 反应堆MCNP 径向截面

1.2 堆芯核素份额计算

堆芯各结构部件的材料为：燃料棒和可燃毒物棒包壳、导向管为Zr-4 合金；围板、吊篮为OCr18Ni10Ti；压力容器焊层为奥氏体不锈钢；压力容器、上下管座为508-Ⅲ钢。除堆芯结构材料外，准确描述燃料棒、可燃毒物棒以及控制棒中核素含量对于临界计算至关重要。由于控制棒由天然铪(Hf-nat)组成，因此无需计算核素份额。对于燃料棒以及可燃毒物棒，设计资料只给出0EFPD燃耗时刻堆芯主要同位素235U、238U、10B 的含量分别为100 300、2 848 000 和151 g。因此，其他核素含量需要根据给出的核素含量进行推算。

0EFPD 燃耗时刻燃料棒中所含的核素有235U、238U 和16O，可燃毒物棒中所含核素有10B、11B、12C和Zr-nat，忽略Zr-2 合金中的微量元素。现根据设计参数计算燃料棒中以及可燃毒物棒中各同位素质量份额。

堆芯中UO2总含量为

式中：na为堆芯组件数，nf为组件中燃料棒数，Vpin为每根燃料棒的体积，ρf为燃料密度设计值。

UO2燃料中235U 的富集度为

式中C5为U 中235U 的核子数份额。

UO2燃料中235U 的质量份额为

由此可以计算出堆芯中235U 的总含量为

由于根据堆芯物理参数计算出的235U 的总含量与设计资料给出的值之间的相对误差为9.78%，为了保证堆芯235U 总量与设计值相符，需要对燃料密度进行如下修正：

式中：ρf′为修正后的燃料密度，M5为设计资料中给出的堆芯中235U 的总含量。

通过式(2)可以计算出修正后的燃料密度为9.361 1 g/cm3。

UO2燃料中238U 的质量份额为

UO2燃料中16O 的质量份额为

根据式(1)、式(3)和式(4)计算出的235U、238U、16O这3 种同位素的质量份额分别为0.029 97、0.851 47和0.118 56。

对于可燃毒物，首先计算单位体积可燃毒物中B 的含量：

式中：np为堆芯内可燃毒物棒总数；Vp为每根可燃毒物棒芯体的体积；MB为堆芯中B 的总量，MB=M10/ε10，其中，M10为堆芯中10B 的总量，ε10为B-nat中10B 的丰度。

然后便可以计算出单位体积可燃毒物中10B、11B、12C 的含量分别为

式中：CC和CB分别为MCNP5 手册中给出的核级B4C 中的C 和B 的核子数密度，分别取0.217 4 和0.782 6；AC和AB分别为12C、B-nat 的原子质量数，分别取12 和10.8。

最后，单位体积可燃毒物中Zr 的含量为

式中 ρp为可燃毒物芯体的密度。

由此便可以计算出可燃毒物中10B、11B、12C、Zr-nat 这4 种元素的质量份额分别为0.000 900、0.003 957、0.001 509 和0.993 634。

将上述材料份额分别填写到MCNP 中对应栅元的材料卡Mn 中完成材料描述。由于MCNP 使用的是NJOY 处理的ACE 格式截面库，NJOY 处理的截面库已经包括在截面库评价温度下的弹性、俘获、裂变和其他低阈值吸收截面(小于1 eV)的多普勒展宽，因此为了考虑多普勒效应，本研究在选择堆芯不同构件中相关核素的截面库时，选择评价温度与该构件温度最相近的截面库，对于235U、238U 还要考虑截面库是否包含缓发中子截面数据。由于临界计算基本不考虑粒子的抽样问题，因此，设定堆芯及其周围结构的栅元重要性IMP:N=1；设定安全壳以外区域中的栅元重要性IMP:N=0，某个栅元的IMP:N=0 表示中子一旦进入该栅元域便被“杀死”，中子历史终止。

1.3 不可分辨共振区概率表

在不可分辨的共振区内(例如，ENDF/B-VI中，235U 为2.25～25 keV、238U 为10～149.03 keV、239Pu 为2.5～30 keV)，连续能量中子截面显示为能量的平滑函数。这是因为共振峰的距离太近以至于无法分辨，而并不是共振峰的缺失。平滑变化的截面未能考虑能量自屏蔽效应，但它对于谱峰在或接近不可分辨共振区系统可能是重要的。

基于分层抽样技术，MCNP5 生成了不可分辨共振区截面的概率表。只要单次碰撞的平均能量损失比共振的平均宽度大得多，也就是说，如果窄共振近似有效，那么从这些概率表中对随机游动截面抽样就是一种有效的物理近似。

从概率表中抽样截面很简单。在每一个入射能量处，都有一个累积概率表，以及与这些概率相对应的近总、弹性、裂变和辐射俘获截面值以及热沉积数。这些数据补充了通常的连续数据。如果概率表关闭(在MCNP5 输入文件数据卡中添加PHYS:N J J 1 J 卡)，则使用平滑截面；但是如果概率表打开(MCNP5 默认)，并且如果碰撞能量在不可分辨共振区内，则从概率表中抽样截面。

1.4 自由气体热散射模型

中子和原子之间的碰撞受到原子热运动的影响，在大多数情况下，碰撞也受到附近其他原子的存在的影响。MCNP 采用基于自由气体近似的热散射来解释热运动。由于栅元温度会影响弹性散射截面，而自由气体热散射的默认值为室温，因此，只要栅元处于其他温度，就必须在TMP 卡上给出这些值。

MCNP 程序在加载截面时，如果数据库中的总截面和弹性截面的温度与模型中TMP 卡指定的栅元温度不同，则MCNP 将根据TMP 卡指定的栅元温度对数据库中的总截面和弹性截面进行热调整。如果不同的栅元有不同的温度，MCNP 首先将截面调整到0K 下的截面，然后在运输过程中将其再次调整到栅元温度对应的截面。由于反应堆堆芯温度高于室温，调整后的弹性散射截面会增大。本研究在TMP 卡上指定堆芯各栅元的平均温度来模拟温度对自由气体热散射模型的影响。

1.5 热中子S(α,β)模型

S(α,β)模型考虑了化学(分子)结合和晶体效应，当中子的能量足够低时(通常针对能量小于4 ev 的中子)，这些效应就变得很重要。S(α,β)模型允许2 个过程：1)非弹性散射，截面 σin和相互耦合的能量、散射角均源于ENDFS(α,β)散射定律；2)弹性散射，截面σel和散射角由固体晶格参数导出。

对于材料卡Mn 中指定的核素，可以使用自由气体模型直到低于某一能量，此时若S(α,β)数据存在，S(α,β)模型将自动取代自由气体模型。通常，自由气体模型用于材料的每一种同位素，直到中子能量低于几个电子伏特，然后对MTn 卡上指定物质的同位素进行S(α,β)处理。本研究将堆芯中冷却剂轻水设置为S(α,β)模型，选择轻水的S(α,β)表lwtr.03t，lwtr.03t 为500 K 时轻水的热中子S(α,β)表，该表只提供轻水中1H 的热中子非弹性散射数据，16O 仍然是默认的自由气体热散射模型。

2 临界计算

临界计算最重要的是判断裂变源分布与Keff是否收敛，因此MCNP5 对于临界计算要求如下[11-12]：1)对所有可裂变物质进行充分抽样；2)确保在有效循环前获得基本特征值；3)最终的统计标准差应该小于0.001。MCNP5 具有检查功能以确保满足前2 个要求，如果不满足MCNP5 将在输出文件中给出警告信息。

为了对可裂变物质进行充分地抽样，并尽量减少模拟粒子数，本研究在KSRC 卡上给出每个燃料组件的中心燃料棒的中心点作为初始源位置。为确保在有效循环前获得基本特征值，本研究在KCODE 卡中设置跳过的无效循环次数为100。计算结果显示，相关设置均满足要求。

对于统计标准差，其大小主要由计算量决定。用MCNP5 进行临界计算的计算量由KCODE卡中的总循环次数KCT 和每次循环的源粒子数NSRCK 决定。中心极限定理指出，在某置信水平下，统计标准差与模拟粒子数的平方成反比，因此欲使统计标准差减半，需将模拟粒子数增大为原来的4 倍，计算时间也会增大为原来的4 倍。本研究综合考虑计算时间与计算精度，选择KCT为500，NSRCK 为100 000，计算结果显示标准差约为0.000 1，远小于0.001，可以认为足够精确。

在设置好临界源后，分别模拟不考虑不可分辨共振区概率表、自由气体热散射的温度修正、热中子S(α,β)模型等因素以及考虑上述所有因素下的堆芯临界，并使用基于网格计数的*FMESHn卡记录堆芯径向各组件功率，使用基于栅元计数的Fn 卡和En 卡记录堆芯能谱。

本研究中假设考虑所有因素的模型为case1，关闭不可分辨共振区概率表的模型为case2，去掉TMP 卡(忽略栅元温度对自由气体热散射模型的影响)的模型为case3，去掉MTn 卡的模型为case4。将case2、case3、case4 的计算结果分别与case1 的计算结果进行比较，分析各因素的影响。

表1 给出了case1—case4 的Keff计算结果以及case2、case3、case4 的Keff计算结果分别与case1的Keff计算结果的相对误差。可以看出：case3 和case4 的相对误差均大于MCNP5 给出的统计误差，而case1 的相对误差小于MCNP5 给出的统计误差。这表明不可分辨共振区概率表对Keff的计算结果的影响可以忽略，自由气体热散射的温度修正、热中子S(α,β)模型对Keff的计算结果的影响不可忽略。其中，case4 的相对误差最大，而且得出了非保守的Keff计算结果。所以对于反应堆临界安全分析必须考虑热中子S(α,β)模型。

表1 Keff 计算结果

计算case1—case4 的堆芯归一化能谱，并使用全谱绝对误差积分以及各能区相对误差的二范数说明各因素对能谱的影响程度。case2、case3、case4 能谱与case1 能谱的绝对误差的全谱积分分别为4.435 7×10-6、4.943 4×10-5和1.894 6×10-4，低能区(小于1 eV)、中能区(1 eV～0.1 MeV)和高能区(0.1 MeV～20 MeV)相对误差的均方根值分别为(0.055 0，0.009 9，0.061 9)、(0.057 1，0.012 1，0.066 8)和(0.310 0，0.034 4，0.083 5)。case4的误差明显大于case2 和case3 的误差，而且尤其在低能区误差最为显著，相对误差的均方根值高达0.31，这也是其Keff计算结果非保守的原因。所以在计算堆芯能谱时必须考虑自由气体热散射的温度修正以及热中子S(α,β)模型，而其他2 个因素相对于热中子S(α,β)模型可以忽略。

表2 给出了case1—case4 的1/4 堆芯各组件归一化相对功率分布计算结果，以及case2、case3、case4 的各组件功率分布计算结果分别与case1 的各组件功率分布计算结果的相对误差。计算得到case2、case3 和case4 堆芯功率分布相对误差的二范数和正向无穷范数分别为(0.029 4，0.012 59)、(0.031 3，0.015 89)和(0.041 0，0.018 72)，统计误差的二范数为0.004 5。可以看出绝大部分相对误差都大于甚至远大于统计误差，并且相对误差的二范数远大于统计误差的二范数。这表明不可分辨共振区概率表、自由气体热散射的温度修正、热中子S(α,β)模型对堆芯功率分布均有着较大的影响，在计算堆芯功率时不可忽略。并且，case4 相对误差的二范数和无穷范数最大，case3 次之，case2 最小，表明热中子S(α,β)模型对堆芯功率分布影响最大，自由气体热散射的温度修正次之，不可分辨共振区概率表对堆芯功率分布影响最小。

表2 1/4 堆芯各组件相对功率分布值

3 结论

本文利用MCNP5 对某压水反应堆堆芯进行精细的几何描述，并计算了0EFPD 燃耗时刻堆芯燃料棒及可燃毒物棒中主要核素的质量份额，完成了材料描述。在此基础上，分别计算不考虑不可分辨共振区概率表、自由气体热散射的温度修正、热中子S(α,β)模型等因素下的反应堆有效增值因子、堆芯能谱及功率分布，并以考虑上述所有因素下的计算结果为基准进行比较，得出以下结论。

1)不可分辨共振区概率表对Keff值计算结果的影响可以忽略，自由气体热散射的温度修正、热中子S(α,β)模型对Keff值计算结果的影响不可忽略，其中热中子S(α,β)模型对临界值的影响最为显著，若忽略热中子S(α,β)模型将得到非保守的Keff值。

2)S(α,β)模型对能谱有着显著影响，尤其是在低能区影响最大，其他2 个因素对能谱的影响相对于热中子S(α,β)模型可以忽略。

3)3 个因素对于堆芯功率分布均有较大影响，并且热中子S(α,β)模型对堆芯功率分布影响最大，自由气体热散射的温度修正次之，不可分辨共振区概率表对堆芯功率分布影响最小。