1. 引言

近年来，随着计算机的发展和反应堆全堆芯pin-by-pin精细建模的需求，蒙卡方法越来越受到重视。针对于核反应堆安全分析、反应堆临界与屏蔽设计、辐射医学物理等研究领域，国际上很多国家和单位研制了蒙特卡罗粒子输运计算软件，如国际上的MCNP、GEANT4、EGS、MCBEND、MONK等程序，国内清华大学工程物理系研制的RMC程序 [1]，中科院合肥等离子体所正在研制开发SuperMC程序 [2]，北京应用物理与计算数学研究所正在开发的JMCT中子光子耦合输运程序 [3]。而在三维核与辐射防护仿真领域，采用蒙特卡罗方法模拟计算各种三维复杂几何下的中子输运、光子输运及中子光子耦合输运问题，快速响应三维辐射剂量动态计算并应用于核设施退役工程仿真，目前国内研究还较少 [4]。

本文针对目前三维核与辐射防护及核设施退役工程需求，为实现三维复杂几何建模，快速实现辐射场剂量计算，基于蒙卡方法，跟踪模拟粒子输运与碰撞过程，对辐射场快速计算方法进行研究，采用分裂、轮盘赌、偏倚抽样等三类减方差技巧加速收敛，采用MPI进行并行加速，设计并验证一套适用于辐射场快速计算的蒙卡计算程序NPTS，并应用于辐射场剂量计算及退役研究。

2. 输运过程模拟

不同于确定论计算方法直接求解粒子输运方程，蒙卡方法求解积分输运方程。对于中子、光子等粒子的输运问题，蒙卡方法对单个粒子的随机输运计算，获得大量粒子输运样本，利用统计方法给出某个随机变量数值特征的估计量，作为问题的解，可对粒子通量、泄漏量、有效增殖因子等物理量进行计算 [5] [6] [7]。

通过模拟每个源粒子的历史径迹，从源粒子产生直到终止条件(包含死亡、泄漏、截断、赌等)，粒子每步输运根据输运截面数据进行随机抽样 [8]。通常粒子历史模拟包括如下步骤：源参数抽样、输运过程模拟、碰撞模拟以及结果统计等。通常有两种计算模式：一是固定源模式，广泛应用于屏蔽计算，二是源迭代模式，应用于反应堆临界计算中 [9] [10]。

2.1. 飞行过程模拟

相比于带电粒子模拟，中子和光子在飞行过程中方向不发生变化，故而输运过程模拟仅需要对其飞行距离进行抽样，下面以中子为例说明，给定一个有固定组分的几何体，中子沿飞行方向飞行距离s不发生碰撞，而在s至 $s+\text{d}s$ 之间发生第一次碰撞的概率 $p\left(s\right)\text{d}s$ 为：

$p\left(s\right)\text{d}s={\text{e}}^{-{\Sigma }_{t}s}{\Sigma }_t\text{d}s$ (1)

其中 ${\Sigma }_t$ 是该几何体内的物质的宏观总截面，也可以理解为该物质中中子每发生一次碰撞的平均飞行距离。对公式(1)进行积分，可以得到中子飞过距离l时，发生碰撞或反应的概率为：

$\xi ={\displaystyle {\int }_0^l{\text{e}}^{-{\Sigma }_{t}s}{\Sigma }_t\text{d}s}=1-{\text{e}}^{-{\Sigma }_{t}l}$ (2)

其中 $\xi \in \left[0,1\right)$，如果我们实现从 $\left[0,1\right)$ 随机抽取一个 $\xi$ 值，则通过公式(2)可以计算出中子的飞行距离为： $\xi \le \frac{{\sigma }_a}{{\sigma }_t}$，因为 $1-\xi$ 和 $\xi$ 均为 $\left[0,1\right)$ 的均匀分布，所以 $l=-\frac{1}{{\Sigma }_t}\ln \xi$。

以上即为蒙卡方法中模拟粒子飞行的基本原理，根据指数分布的“无记忆性”特点，屏蔽计算中采用Ray-tracking方法对飞行过程进行模拟，抽样飞行距离。Ray-tracking方法的主要思路为逐层考虑粒子是否在该层介质内发生碰撞的方法来确定飞行距离和碰撞点的位置 [11] [12]。假设中子在系统中处于状态 $\left(\stackrel{\to }{r},E,\stackrel{\to }{\Omega },t,w\right)$，那么中子在 $\stackrel{\to }{r}$ 处沿着方向 $\stackrel{\to }{\Omega }$ 继续飞行，在飞行中将与介质发生碰撞，那么其飞行距离的概率密度函数为：

$f\left(l\right)={\Sigma }_t\left(\stackrel{\to }{r}+l\stackrel{\to }{\Omega },E\right)\exp \left\{-{\displaystyle \underset{0}{\overset{l}{\int }}{\Sigma }_t\left(\stackrel{\to }{r}+l^{\prime }\stackrel{\to }{\Omega },E\right)\text{d}{l}^{\prime }}\right\}$ (3)

蒙卡模拟中，虽然不需要对空间位置、飞行角度以及能量等进行离散，但是对几何体中材料还是以均匀化的方式进行处理，即，认为在特定的几何体中，材料参数是一致的。这样在式(3)中的积分就可以变为求和的形式，如式(4)所示：

${\displaystyle \underset{0}{\overset{l}{\int }}{\Sigma }_t\left(\stackrel{\to }{r}+l^{\prime }\stackrel{\to }{\Omega },E\right)\text{d}{l}^{\prime }}\to {\displaystyle \underset{i=1}{\overset{k-1}{\sum }}{\Sigma }_{t,i}{l}_i}+{\Sigma }_{t,k}\left(l-{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}{l}_i}\right)$ (4)

其中l满足关系式： ${\displaystyle \underset{i=1}{\overset{k-1}{\sum }}{l}_i}<l\le {\displaystyle \underset{i=1}{\overset{k}{\sum }}{l}_i}$，其中 $l_i$ 为飞行方向上每层介质的厚度。

图1为粒子飞行时穿过多层介质的示意图。利用ray-tracking方法进行计算时，首先考虑第一层介质，选取随机数 $\xi$ ( $\xi$ 为 $\left[0,1\right)$ 上均匀分布的随机数)，在 ${\Sigma }_{t,1}{\text{e}}^{{\Sigma }_{t,1}l}$ 的分布中抽样得到 $l=-\ln \xi /{\Sigma }_{t,1}$，若 $l\le l_1$ 则认为中子在第一层介质内发生碰撞并且飞行距离为l；若 $l>l_1$，则认为中子在第一层介质内不发生碰撞，将模拟点移动到飞行路径和1、2介质的交界面处，再从第二层介质开始抽样，直至最终确定碰撞位置，得到相应的飞行距离。

2.2. 碰撞过程模拟

如果中子穿过介质时与几何体内材料发生了碰撞，接下来将按照以下顺序依次进行抽样模拟：

a) 确定发生碰撞的核素；

b) 如果中子的能量足够低，且存在合适的 $s\left(\alpha ,\beta \right)$ 模型可以使用，则采用 $s\left(\alpha ,\beta \right)$ 模型；否则，对于低能的中子需要抽样靶核的速度，使用自由气体模型进行处理；

c) 根据模型确定是否产生光子，如果为中子–光子耦合输运则产生光子并放入堆栈，进行后续的处理；

d) 处理中子的吸收反应；

e) 如果经b)判断需要使用 $s\left(\alpha ,\beta \right)$ 模型，则使用该模型处理中子和靶核的反应，并计算得到中子的出射方向和能量；如果不使用 $s\left(\alpha ,\beta \right)$ 模型，则判断是发生弹性散射还是非弹性散射等反应(包括裂变)，进一步计算出射中子的方向和能量。

蒙特卡罗方法解粒子输运问题的程序，一般都可分为：源抽样，空间输运过程，碰撞过程，记录过程和结果的处理与输出等部分 [13]，如图2所示。

3. 计算功能实现

3.1. 固定源计算功能

基于MC方法的NPTS程序(Neutron-Photon Transport Simulation program)，在处理固定源计算过程中，在输入粒子类型、能量、数量等计算所需参数后进行计算，计算完成后统计结果并输出，固定源计算主要流程如图3所示。

3.2. 临界计算功能

在源迭代计算模式中，NPTS程序需要模拟中子在反应堆内的增殖问题，源迭代计算主要流程如图4所示。

临界问题以求解 $k_{eff}$ 为特征，即求解中子输运方程的特征值。计算 $k_{eff}$ 包括估计在一代中每个中子所能产生的裂变中子数的平均值，其中一代指的是中子历史从裂变产生，并通过泄漏、俘获或者裂变吸收反应结束。

3.3. 减方差技巧

在蒙卡计算中使用减方差技巧能够有效缩减计算时间提高计算效率 [14]。NPTS程序通过如下形式估计计数器结果：

$\langle T\rangle ={\displaystyle \int \text{d}\stackrel{\to }{r}}{\displaystyle \int \text{d}\stackrel{\to }{v}}{\displaystyle \int \text{d}t}N\left(\stackrel{\to }{r},\stackrel{\to }{v},t\right)T\left(\stackrel{\to }{r},\stackrel{\to }{v},t\right)$ (6)

其中，N为粒子数密度，T为计数量。

NPTS程序计算采用分裂、轮盘赌、偏倚抽样等三类减方差技巧 [15]。几何分裂可以通过中子分裂而使更多的中子输运到计数区，而使更多中子输运到计数区几何分裂将粒子由重要性小的区域进入重要性大的区域时，增大粒子数，同时降低粒子权重，从而保持粒子总权重不变，保证计数结果无偏，该方法通过设定imp (重要)卡而达到减方差目的；通过轮盘赌可减少非计数区的统计，粒子由重要性大的区域进入重要性小的区域时，粒子以一定的概率存活，存活粒子权重增大，从而保持粒子总权重不变，保证计数结果无偏。通过源的方向偏倚使更多的源中子偏向感兴趣的方向，通过源的能量偏倚使源中子具有相对高的能量，使得较高能量的中子更易穿透到屏蔽层深处，偏倚抽样通过增大特定的对结果有重要贡献的粒子数的抽样，降低粒子权重的方式，保证计数结果无偏。

3.4. 并行计算实现

目前常用的并行编程接口的标准有MPI和OpenMP，考虑OpenMP主要针对单主机上多核/多CPU并行计算，而MPI同时协调多台主机间的并行计算，在规模上的可伸缩性较强，考虑NPTS程序并行进程需求，选取进程级MPI并行模式，并行过程采用主从模式(master-slave mode)，master进程负责接受、散发和整理数据，同时也参与计算，slave进程只负责计算输运过程。具体实现上，反应堆蒙卡分析软件采用了单程序多数据形式的并行算法设计，属于粗粒度的并行算法，充分利用多核系统中各个CPU核的独立性，使一个计算进程包含大的计算任务 [16] [17]。

计算开始时，所有进程同时读输入文件及输入文件中涉及到的中子反应截面数据；然后进入代循环，master (主)进程分配任务并散发任务，slave (子)进程接受任务，这个过程包含了为每一个中子产生一个初始随机数的工作；接下来就是对每一个中子进行模拟；各个slave进程将收集到的中子信息及统计结果发送到master进程，master进程收集这些信息，对其处理后，为下一个代循环做准备。

4. 计算验证及分析

4.1. 临界计算验证

为考察NPTS程序在临界计算方面的问题，选取15个MCNP5自带的临界基准题进行计算验证。临界基准题在材料组分方面涵盖了通常反应堆内的裂变材料以及反射层材料，其基准题活性区材料包括U、Pu金属及其溶液。通过临界基准题计算，验证程序针对快谱、热化谱系统的临界模拟能力，同时也考察了程序对于热化库的计算功能，其详细结果见表1。

从表中计算结果可以看出，对于15个基准题临界问题计算，程序计算结果与MCNP5的计算值吻合很好，其临界模型K_eff相对偏差均小于5‰，验证了该程序处理临界计算功能的正确性。

4.2. 并行效率计算分析

p104模型为自制精细出壳伽马能谱(∆E = 0.01 MeV)对比模型。采用程序NPTS与MCNP5对该模型进行能谱计算，结果如图5所示，可以看出两程序计算结果基本一致，在计算时间上NPTS程序计算时间1.38 min，MCNP5计算时间为2.22 min，初步验证采用减方差技巧的有效性，其中 $\text{NPTS-ERR}=\text{abs}\left(\frac{x_1-x_2}{x_1{\sigma }_1}\right)$， $x_1$ 、 $x_2$ 与 ${\sigma }_1$ 为两次随机模拟结果与相对标准差。

为进一步考核程序的并行效率问题，基于上述P104模型，采用多核开展了NPTS程序并行计算效率测试，得到不同并行核数下的计算时间与并行效率，计算结果如表2所示，从并行计算效率可以看出，随着并行核数的增加，并行效率均高于95%，提高了计算的效率。

注1：节省计算时间，1核和2核计算样本数比其它模型低，表2中计算时间已折算至同等样本数；注2：为降低并行平台载荷带来的计算时间波动，128核计算样本数比其它模型高5倍，表2中计算时间已折算至同等样本数。

4.3. 固定源计算验证分析

Be08i模型为Be材料中子慢化模型，该模型几何较简单，为两个同心球体在x轴方向被圆柱所截断，如图6所示给出XY截面示意图，带能量和时间截断，源抽样复杂，包含了空间、角度和能量分布。

整个模型用一个半径为1000 cm的球包围，统计最外层包围球面上的能谱分布，两个程序中子能谱结果对比如图7所示，两个程序计算结果基本一致，如图中红线和黑线。NPTS程序与MCNP5程序两者计算相对偏差最大绝对值与涨落比值偏差基本都在0~3倍之间。

C29i模型为C材料中子慢化模型，该模型源分布与Be08i一致，几何和材料有所不同，其XY截面如图8所示。

两个程序结果对比如图9所示。两个程序计算结果基本一致，NPTS程序与MCNP5程序计算结果相对偏差绝对值与涨落比值偏差基本小于2倍，两者相对偏差最大绝对值与涨落比值偏差均小于4倍。

FE09i模型为铁中子慢化模型，该模型源分布与Be08i一致，几何模型与C29i模型类似。两个程序结果对比如图10所示。两个程序计算结果趋势基本一致，NPTS程序与MCNP5程序计算结果相对偏差绝对值与涨落比值偏差基本在3倍偏差以内，两者相对偏差最大绝对值与涨落比值偏差均小于4倍。

SKYINP2模型为纯伽马输运模型，考虑了不同栅元重要性和伽马韧致辐射。该模型由多个圆锥和同心球所围成，如图11所示，利用NPTS程序计算伽马能谱，其计算结果如图12所示。

由图可以看出，NPTS程序与MCNP5程序计算结果趋势基本一致，NPTS程序与MCNP5程序计算结果相对偏差绝对值与涨落比值偏差基本在2倍附近，两者相对偏差最大绝对值与涨落比值偏差均小于4倍。

5. 结论

本文基于蒙卡方法，跟踪模拟粒子输运与碰撞过程，对辐射场快速计算方法进行研究，设计并验证一套适用于辐射场快速计算的蒙卡计算程序NPTS，经验证分析，结果表明：

1) 考察NPTS程序在临界计算方面的问题时，选取15个临界基准题进行计算验证在临界计算，计算结果表明临界模型K_eff相对偏差均小于5‰，初步验证了临界模型处理的正确性；

2) 选取4个固定源基准题模型验证程序的固定源处理问题，计算结果表明，模型程序计算能谱与蒙卡程序计算结果趋势基本一致，NPTS程序与MCNP5程序计算结果相对偏差绝对值与涨落比值偏差基本小于2倍，两者相对偏差最大绝对值与涨落比值偏差均小于5倍，验证了粒子输运的精准性；

3) 在程序计算速度和并行效率上，NPTS程序采用加速方法，提高选取模型的计算速度，针对本文选用的计算模型百核并行效率高于95%，缩短了程序的计算时间，但在深穿透问题屏蔽计算方面，后续有待进一步深入研究。以上为建立核与辐射防护剂量快速计算体系及退役研究提供参考和依据。

NOTES

^*第一作者。

参考文献