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

MCNP作为核能领域的MC计算程序，发展最为成熟、最具代表性。近年来，国内外学者利用MCNP在反应堆有效增殖因子^[1-3]、反应性^[4-5]、中子通量和功率分布^[3,6-8]以及中子能谱^[9-10]计算等方面做了诸多研究。虽然前人利用MCNP模拟反应堆临界做了一系列工作，但是，在建模过程中，多数学者只关心反应堆堆芯部件的几何结构、材料以及数据库选取，没有关注影响模型准确性的其他因素，如不可分辨共振区概率表、自由气体热散射的温度修正以及热中子 $S(\alpha ,\beta )$ 模型等。因此，本文以某压水反应堆为研究对象，综合分析了上述因素对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所示。堆芯内部结构均按照实际尺寸进行精细化建模，燃料棒气隙、可燃毒物棒气隙和组件间水隙等微小结构也考虑在内。

使用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卡将各个世界体填充的栅格按矩阵形式进行重复，其中参考栅格作为矩阵的 $\left[ {\begin{array}{*{20}{c}} 0&0&0 \end{array}} \right]$ 元素。为所有栅格分配世界体编号u=1，最后使用fill=1将栅格填入组件中。图2为MCNP5构建的E类组件模型的径向截面图，其中不同的颜色代表不同的材料。其他E类组件只需要使用like…but…卡进行复制和坐标变换。除燃料组件以外其他堆芯结构(如围板、吊篮、压力容器等)的几何描述相对比较简单，因此不再详述。图3为MCNP构建的堆芯整体模型的径向截面图。

1.2 堆芯核素份额计算

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

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

堆芯中UO₂总含量为

${M_{\rm{f}}}{\rm{ = }}{n_{\rm{a}}}{n_{\rm{f}}}{V_{{\rm{pin}}}}{\rho _{\rm{f}}}$

式中： ${n_{\rm{a}}}$ 为堆芯组件数， ${n_{\rm{f}}}$ 为组件中燃料棒数， ${V_{{\rm{pin}}}}$ 为每根燃料棒的体积， ${\rho _{\rm{f}}}$ 为燃料密度设计值。

UO₂燃料中²³⁵U的富集度为

${\varepsilon _5}{\rm{ = }}\frac{{235{C_5}}}{{235{C_5} + 238(1 - {C_5})}}$

式中 ${C_5}$ 为U中²³⁵U的核子数份额。

UO₂燃料中²³⁵U的质量份额为

${\omega _5}{\rm{ = }}\frac{{235{C_5}}}{{235{C_5} + 238(1 - {C_5}) + 32}}$

(1)

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

${M_5}^\prime {\rm{ = }}{\omega _5}{M_{\rm{f}}}$

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

${\rho _{\rm{f}}}^\prime = \frac{{{M_5}}}{{{M_5}^\prime }}\rho _{\rm{f}}^{}$

(2)

式中： $\;{\rho _{\rm{f}}}^\prime $ 为修正后的燃料密度， ${M_5}$ 为设计资料中给出的堆芯中²³⁵U的总含量。

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

UO₂燃料中²³⁸U的质量份额为

${\omega _8}{\rm{ = }}\frac{{238(1 - {C_5})}}{{235{C_5} + 238(1 - {C_5}) + 32}}$

(3)

UO₂燃料中¹⁶O的质量份额为

${\omega _{16}}{\rm{ = }}\frac{{32}}{{235{C_5} + 238(1 - {C_5}) + 32}}$

(4)

根据式(1)、式(3)和式(4)计算出的²³⁵U、²³⁸U、¹⁶O这3种同位素的质量份额分别为0.029 97、0.851 47和0.118 56。

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

${m_{\rm{B}}}{\rm{ = }}\frac{{{M_{\rm{B}}}}}{{{n_{\rm{p}}}{V_{\rm{p}}}}}$

式中： ${n_{\rm{p}}}$ 为堆芯内可燃毒物棒总数； ${V_{\rm{p}}}$ 为每根可燃毒物棒芯体的体积； ${M_{\rm{B}}}$ 为堆芯中B的总量， ${M_{\rm{B}}}{\rm{ = }} $ $ {{{M_{10}}} / {{\varepsilon _{10}}}}$ ，其中， ${M_{10}}$ 为堆芯中¹⁰B的总量， ${\varepsilon _{10}}$ 为B-nat中¹⁰B的丰度。

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

${m_{10}}{\rm{ = }}{m_{\rm{B}}}{\varepsilon _{10}}$

${m_{11}} = {m_{\rm{B}}}(1 - {\varepsilon _{10}})$

${m_{\rm{C}}} = \frac{{{C_{\rm{C}}}{A_{\rm{C}}}}}{{{C_{\rm{B}}}{A_{\rm{B}}}}}{m_{\rm{B}}}$

式中： ${C_{\rm{C}}}$ 和 ${C_{\rm{B}}}$ 分别为MCNP5手册中给出的核级B₄C中的C和B的核子数密度，分别取0.217 4和0.782 6； ${A_{\rm{C}}}$ 和 ${A_{\rm{B}}}$ 分别为¹²C、B-nat的原子质量数，分别取12和10.8。

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

${m_{{\rm{Zr}}}} = {\rho _{\rm{p}}} - {m_{\rm{B}}} - {m_{\rm{C}}}$

式中 ${\rho _{\rm{p}}}$ 为可燃毒物芯体的密度。

由此便可以计算出可燃毒物中¹⁰B、¹¹B、¹²C、Zr-nat这4种元素的质量份额分别为0.000900、0.003957、0.001509和0.993634。

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

1.3 不可分辨共振区概率表

在不可分辨的共振区内(例如，ENDF/B-VI中，²³⁵U为2.25~25 keV、²³⁸U为10~149.03 keV、²³⁹Pu为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(\alpha ,\beta )$ 模型

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

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

2 临界计算

临界计算最重要的是判断裂变源分布与K_eff是否收敛，因此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(\alpha ,\beta )$ 模型等因素以及考虑上述所有因素下的堆芯临界，并使用基于网格计数的*FMESHn卡记录堆芯径向各组件功率，使用基于栅元计数的Fn卡和En卡记录堆芯能谱。

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

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

计算case1—case4的堆芯归一化能谱，并使用全谱绝对误差积分以及各能区相对误差的二范数说明各因素对能谱的影响程度。case2、case3、case4能谱与case1能谱的绝对误差的全谱积分分别为 $4.435\;7 \times {10^{ - 6}}$ 、 $4.943\;4 \times {10^{{\rm{ - }}5}}$ 和 $1.894\;6 \times {10^{{\rm{ - }}4}}$ ，低能区(小于1 eV)、中能区(1 eV~0.1 MeV)和高能区(0.1 MeV~20 MeV)相对误差的均方根值分别为(0.0550，0.0099，0.0619)、(0.0571，0.0121，0.0668)和(0.3100，0.0344，0.0835)。case4的误差明显大于case2和case3的误差，而且尤其在低能区误差最为显著，相对误差的均方根值高达0.31，这也是其K_eff计算结果非保守的原因。所以在计算堆芯能谱时必须考虑自由气体热散射的温度修正以及热中子 $S(\alpha ,\beta )$ 模型，而其他2个因素相对于热中子 $S(\alpha ,\beta )$ 模型可以忽略。

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

3 结论

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

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

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

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