300＃研究堆安全棒中子注量率计算中的减方差方法对比及应用

唐凤平，刘耀光，杨万奎，杨 鑫，2

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

300＃研究堆的安全棒为碳化硼吸收体，为研究碳化硼在整个反应堆寿期内的辐照性能，需要获取整根碳化硼吸收体沿高度方向的中子注量率分布情况。但安全棒顶端离堆芯较远，且有效吸收体体积较小，此类问题即属于远距离小体积计数，经计算，直接使用径迹长度或面通量计数根本得不到安全棒顶端的计算结果，因此需要考虑减方差方法［1］。类似的小体积远距离输运问题是蒙特卡罗方法应用中的一个经典问题，这是由粒子到达计数区域附近的概率较小所致。清华大学已开展内照射小器官剂量计算中减方差技巧的比较和应用［2］，通过此方法获得了较好的结果。本文针对该安全棒顶端的计数问题开展若干减方差方法的对比研究。

1 堆芯布置及安全棒位置

300＃研究堆额定功率3 MW，是一座游泳池式研究堆，堆芯包括活性区（燃料组件）、反射层（铍组件和石墨组件）以及垂直孔道。堆芯布置如图1所示［3］。堆芯中心为燃料组件，铍组件和石墨组件作为反射层。其中燃料组件22个（16根燃料元件的燃料组件12个，15根燃料元件的燃料组件10个），铍组件30个，石墨组件33个，控制棒12根（位于石墨盒的控制棒导管内和15根燃料元件的元件盒内），垂直空腔孔道1个，垂直水腔孔道4个。

图1 反应堆堆芯布置及安全棒位置示意图Fig.1 Schematic diagram of reactor core loading and safety rod's location

安全棒为碳化硼吸收体，有效吸收尺寸为φ17.5mm×600 mm。1AB 代表1号安全棒，2AB代表2号安全棒，1AB位于中央孔道的左侧，2AB位于中央孔道的右侧，其高度位置示意图如图2所示，即反应堆启动时将两根安全棒提升至反应堆顶部，之后再提升其余棒使反应堆达临界。

图2 安全棒高度位置示意图Fig.2 Schematic diagram of safety rod's location in height

2 减方差原理

2.1 一般原理

按减方差原理，可将若干减方差方法分为4类［1］：

1）截断方法

该类方法是最简单直接的方法，主要包括几何截断、能量截断和时间截断。

2）数量控制方法

该类方法通过粒子分裂和俄罗斯轮盘赌来控制抽样数量。在重要的区域增加抽样，但每个粒子的权重降低；在不重要的区域减少抽样，但每个粒子的权重增大。该方法通过调整权重以保证结果无偏。主要包括几何、能量、时间相关的赌分裂、权重截断以及权窗。

3）修正抽样方法

该类方法通过改变统计抽样以增加每个粒子的计数率。对任意随机事件，可从任意分布进行抽样，而不是严格按照物理概率抽样，前提是需调整粒子权重来进行补偿修正。因此，修正后的抽样方法使得抽样可按希望的方向发射粒子、粒子可进入感兴趣的相空间（时间、能量），以此改变碰撞位置和类型。主要包括指数变换、隐俘获、强迫碰撞、源偏倚和中子引发的光子产生偏倚。

4）部分确定论方法

该方法使用类似确定论方法来避开标准随机游走过程。例如：下一事件估计、控制随机数序列。MCNP 中该类方法包括点探测器、DXTRAN以及相关抽样。

成功实现MCNP 的减方差方法通常是较为困难的，且艺术技巧性高于科学性［1］。每项降方差方法均有其自身的优点，但也有部分通用的原则。该通用原则就是要对物理问题有比较清晰的理解，依据物理本质来选取相应的方差降低方法。当选定了一项方法后，需指定适合的参数，这通常较选择方法更困难。这些参数的首次猜测通常来自于相似问题的计算经验，或同一问题的相似计算。指定好参数后，还需进行一系列短时间的运行以监测这些参数的有效性。

MCNP输出文件包含了有助于了解抽样的许多有用信息，需查看以下信息来检测方差降低方法参数设置的有效性：改善了计数粒子的抽样；各种降方差方法协作运行，不会严重相互影响；FOM 表稳定，该表可检测抽样不足；结果不会出现明显错误。在确定降方差方法对计算问题能起到改善作用后，用户还需进行多次短时间运行，以此逐步改善参数。这其中参数设置有效性的评判标准主要包括相对误差R和品质因子FOM。

相对误差R 是MCNP 计算结果可信度的评判标准，其定义为估计值的标准差与估计值的比值，即：

品质因子FOM 为MCNP 计算效率的评判标准，其定义为：

其中，T 为计算耗时。实际上，可从计数波动表（TFC）中监测相对误差R 和品质因子FOM 随粒子数的变化趋势。

2.2 3种减方差方法

为获取安全棒顶端的中子注量率情况，根据目前的建模方式，即整根棒为1个栅元，整个上方水池为1个栅元，选取3种减方差方法：点探测器、点探测器＋DXTRAN、点探测器＋强迫碰撞。

1）探测器贡献［1］

由于粒子精确输运到某点的概率是非常小的，因此使用伪粒子来定向到该点。每产生1个粒子（可以是源粒子，也可以是碰撞后的粒子），就要求有1个伪粒子在指定的空间位置点被计数。对小体积计数而言，很可能由于粒子到达该体积的概率极低而不能得到径迹长度估计或穿过表面的通量估计，而点探测器仍可得到较好的通量估计，这就是其一大优点。正是由于每产生1个粒子就要在点探测器位置处进行计数操作，故将影响整体计算效率。

2）DXTRAN［1］

DXTRAN 即 确 定 性 输 运（deterministic transport），它同时偏倚散射方向和源的方向。由于粒子几乎散射不到一些小区域内，因此会造成小区域抽样不充分。为改善这个情形，可在输入文件中指定1 个DXTRAN 球，它需包含该小区域。一旦粒子在球外发生碰撞，或球外有源，DXTRAN 技巧就会生成1 个特殊的DXTRAN 粒子，且会确定性地散射到这个DXTRAN 球内。

3）强迫碰撞［1］

强迫碰撞方法将粒子分裂为碰撞部分和非碰撞部分。碰撞部分在当前栅元内发生强迫碰撞，非碰撞部分离开该栅元且不发生碰撞，但该部分是暂时存储起来，等到该粒子的径迹在栅元边界上时才继续追踪。非碰撞部分继续游走时的权为W0e-Σtd，其中W0为粒子初始权，Σt为宏观总截面，d 为距离。即非碰撞部分的权等于当前粒子权乘以无碰撞地离开该栅元的概率。相应地，强迫碰撞部分的权为W0（1-e-Σtd），即总的权是守恒的。

3 计算结果分析

3.1 初步对比分析

为获取安全棒最顶端的中子注量率情况，首先用较少的粒子数计算，检验3种方法的优劣。主要考察3种方法在相同计算条件下的相对误差R 和品质因子FOM 随粒子数nps的变化趋势。计算方式为基于MPICH 的消息传递并行计算［4］，计算条件为：1）相同的临界计算控制卡kcode 3000 1.0 10 40；2）相同的初始源，均来自同一个srctp源分布文件；3）相同的随机数序列和随机数种子；4）相同的栅元重要性，均设为1。

计算结果如图3和4所示。图中，F5代表点探测器计数；F5＋DXTRAN 代表点探测器加DXTRAN 减方差方法；F5＋FCL 代表点探测器加强迫碰撞减方差方法。

图3 3种减方差方法的计算相对误差Fig.3 Relative error of three variance reduction methods

图4 3种减方差方法的FOMFig.4 FOM of three variance reduction methods

由图3和4的计算结果可知，在相同的计算条件下，点探测器加强迫碰撞的相对误差R较小，品质因子FOM 较高，其减方差效果较好。故采用该减方差组合进行增加粒子数计算，计算结果的相对误差及品质因子如图5所示。

图5 增加粒子数运行结果Fig.5 Run result of increasing particle histories

由图5可知，即便是在增加粒子数后，相对误差R 也未明显减小，呈现波动现象，同时，品质因子FOM 逐渐下降，即计算效率越来越低。且在粒子数增加到一定程度后，相对误差发生阶跃升高，计算获得的结果数据发散。

由此可知，针对本实际问题，安全棒最顶端位置采用点探测器＋强迫碰撞的方式并不能达到预期的减方差效果，故需要采取其他方法，如赌分裂，将粒子引至感兴趣的区域，但由于当前的几何模型中安全棒是由1个栅元定义，且堆芯上部水层也是由1个栅元定义，直接指定栅元重要性并不能实现减方差效果，故考虑将其分层建模。

3.2 分层建模

通过重新建模，将堆芯上部水反射层分层建模，并通过重要性设置，希望进入顶层水的中子径迹数增多［5］。这种方法的优点是，当划分的区域和确定的重要性合适时，可使更多的粒子被引到感兴趣区域的附近，从而可缓解常规随机游动的粒子很少能到达计数区域附近的问题。为便于查看减方差效果，将安全棒顶端的小块按样品研究所需的大小（φ17.5 mm×30mm）单独建模，称该栅元为目标栅元。水反射层分层建模如图6所示。

栅元重要性设置、中子在各水层和目标栅元中的径迹及碰撞情况列于表1。第1 次计算时，将远离堆芯的水层栅元重要性依次加倍，并将目标栅元与其所在水层的重要性均设置为最高的131 072。计算结果显示各水层中的中子径迹与碰撞情况都得到了提升，但目标栅元中的中子径迹却较少。为实现目标栅元中的中子径迹与其同高度的水层相当，将目标栅元重要性调高为149 000／70×131 072＝2.79×108。从新的重要性分配下的第2次计算结果可知，目标栅元中的中子径迹及碰撞数均达到较高值，起到了将中子引至目标栅元的效果。

再增加粒子数计算，以临界控制条件kcode 3000 1.0 10 6010进行计算，对目标栅元采用径迹长度体通量计数（F4），得到的计算结果相对误差及品质因子如图7所示。该体通量计数结果通过全部10个判断标准，并从图7可知：相对误差随运行粒子数的增加逐渐减小，且减小趋势符合；品质因子FOM 虽整体较小，计算效率相对较低，但在收敛后基本保持稳定，也说明了该计算的稳定性。总体来讲，依靠分层建模，进行栅元重要性分配，将粒子引向目标栅元，实现了减方差的目的。

图6 水反射层分层建模Fig.6 Multiple layers modeling of water reflector

表1 栅元中的中子径迹情况Table 1 Neutron track activity in cells

图7 分层建模增加粒子数计算结果Fig.7 Run result of multiple layers modeling after increasing particle histories

4 结论

根据本研究堆几何模型特点，首先，选取了点探测器、点探测器＋DXTRAN、点探测器＋强迫碰撞3种减方差方法进行对比计算，并对其中效果较好的点探测器＋强迫碰撞方法进行了进一步的增加粒子数计算，计算结果显示：随着粒子数的增加，计算发散，未起到减方差的作用。之后，考虑使用栅元重要性分配，将中子引向目标栅元。为实现此目的，将原有几何模型重新分层建模，并分配适当的栅元重要性。计数结果表明：计数相对误差在5%以内，品质因子保持稳定，实现了减方差的目的。

［1］ JUDITH F B.MCNP：A general Monte Carlo N-particle transport code version 4C，LA-13709-M［R］.USA：Los Alamos National Laboratory，2000.

［2］ 张崭，李君利，武祯，等.内照射小器官剂量计算中减方差技巧的比较和应用［J］.清华大学学报：自然科学版，2007，47（S1）：1 051-1 056.ZHANG Zhan，LI Junli，WU Zhen，et al.Comparison and application of variance-reduction techniques used in internal radiation dose calculations for small organs［J］.Journal of Tsinghua University：Science and Technology，2007，47（S1）：1 051-1 056（in Chinese）.

［3］ 窦海峰，代君龙.SPRR-300 反应堆辐照孔道中子注量率的MCNP 程序计算［J］.核动力工程，2006，27（1）：51-54.DOU Haifeng，DAI Junlong.Calculation of neutron flux in SPRR-300reactor irradiation channel with MCNP code［J］.Nuclear Power Engineering，2006，27（1）：51-54（in Chinese）.

［4］ 杨万奎，曾和荣，冷军，等.300＃研究堆首炉中央孔道中子通量密度计算［J］.强激光与粒子束，2012，24（12）：3 000-3 005.YANG Wankui，ZENG Herong，LENG Jun，et al.Neutron flux calculation for central channel in first cycle of SPRR-300［J］.High Power Laser and Particle Beams，2012，24（12）：3 000-3 005（in Chinese）.

［5］ THOMAS E B.A sample problem for variance reduction in MCNP，LA-10363-MS［R］.USA：Los Alamos National Laboratory，1985.