MCNP源子程序在停机剂量率计算中的应用

吴明昌，吴亮亮，郝丽娟，郑华庆，宋 婧，FDS团队

（1.中国科学技术大学，安徽 合肥 230027；2.中国科学院核能安全技术研究所，安徽 合肥230031）

托卡马克核聚变装置在进行等离子体放电实验时，将产生聚变中子。装置停止运行后，结构材料受到这些聚变中子照射产生衰变光子，从而引起核聚变装置的停机辐射剂量。近年来随着核聚变装置的告诉发展、实验功率水平的不断提高，D－T或D－D等离子体实验产生大量聚变中子，装置结构材料的活化情况严重，导致装置周围的停机剂量难以在较短的时间内衰减到可以忽略的程度。装置运行停止后，实验人员有时需要进入装置内部进行实验测量工作，较高的辐射剂量会对实验人员的人身安全产生影响。因此，需要准确地评估实验停止以后装置周围的停机剂量率水平，为装置辐射屏蔽系统的设计、实验方案的制订以及实验人员允许进入装置内部的时间等提供重要的参考[1］。

随着对停机剂量率分析研究工作的深入，文献[2］提出了一种基于网格计数的停机剂量率计算方法——严格二步法，用于精确计算几何结构复杂的托卡马克装置的停堆剂量率。根据该方法，FDS团队开发了一套精确的停机剂量率计算程序，并将其集成进FDS团队开发的大型集成中子学计算分析系统 VisualBUS[3－4］中。该程序通过耦合三维蒙特卡罗程序MCNP[5］和欧洲活化计算程序 FISPACT[6］，经过一系列中间计算，最终获得装置停机后周围空间的剂量率分布。

由于聚变装置几何和材料的高度复杂性，其基于网格计数的停机剂量率计算中涉及的光子源的栅元数目超过了1 000个，难以使用MCNP的标准源卡描述。为了解决SDEF卡对源栅元数目的限制，本文使用了MCNP自定义源子程序方法对衰变光子源进行直接抽样。

1 停机剂量率的计算

聚变堆的停机剂量率是由聚变中子活化的结构材料所发射出的光子造成的。根据直接二步法，停机剂量率的计算过程如下：

（1）中子输运计算，计算运行时中子注量率的空间分布；

（2）活化计算，使用第一步得到的中子注量率分布，求出衰变光子在时间、空间及能量上的分布；

（3）光子输运计算，根据第二步计算得到的衰变光子分布，求出装置停止运行后，其系统在不同时间段的空间剂量率分布。

2 停机剂量率源子程序的设计实现

2.1 总体设计

中子注量率分布与光子分布的结果的精确性直接影响着剂量率计算结果的精确性。为了精确计算剂量率，停机剂量率计算利用MCNP网格计算能力，将停机剂量计算所涉及的几何划分成若干网格，光子计算时网格被视为光子源栅元。当网格划分的越多越精细时，中子注量率分布以及光子分布也就计算得越精确，从而剂量率结果便就越精确。当网格数目（也即光子计算时源栅元数目）过多时，MCNP的SDEF卡无法描述，系统需要通过MNCP源子程序对光子的分布直接进行抽样。

每个网格中的材料因受到的中子辐射程度不同，被活化的程度也不同，光子源的分布特性也不相同。因此，在使用源子程序对光子源进行抽样时，应先根据每个网格的光子强度抽样决定光子所在的网格（每个网格的光子强度由第二步计算得到，存储在sourceShape与source文件中，其计算方法与格式见文献[7］）。每个网格内的材料是均匀的，这些材料被活化后出射的衰变光子也是在网格内均匀分布，因此光子的位置在网格内均匀抽样。大部分衰变核素都具有各向同性特点，光子的方向根据该特点按照各向同性处理。然后根据第二步计算得到的光子能谱抽样光子的能量，最后根据光子的位置计算光子所在的栅元。源子程序的流程如图1所示。

2.2 随机抽样的基本方法

源子程序方法是根据源的分布规律抽样而得到源粒子的位置、方向、能量、权重、寿命的方法。

由已知分布的随机抽样，是借助于随机数进行的，MCNP提供了RANG（）函数产生随机数。该函数可以产生[0，1）均匀分布的伪随机数。用户在编写源子程序依据源的分布规律对源的状态变量进行随机抽样时，可以调用该函数。

图1 源子程序流程Fig.1 Procedure of Source subroutine

对于任意离散型分布

其中x1，x2，…为离散型分布函数的跳跃点，P1，P2，…为相应的概率。F（x）的直接抽样方法为对于连续型分布，如果分布函数F（x）的反函数F－1（x）存在。则抽样方法是

其中ξ是[0，1）均匀分布的随机数，用户可以在源子程序中调用RANG（）产生该随机数。

2.3 源子程序变量的抽样

停机剂量率计算中，衰变光子的状态取决于光子出射的网格，所以源子程序所涉及第一个随机抽样过程是确定源粒子所处的网格。这是一个离散型随机分布，参考抽样公式（2），其抽样方法为

当确定了源粒子所在的网格后，就根据该网格的位置和尺寸在网格均匀抽样以确定源粒子的具体位置（xxx，yyy，zzz）。

下一步即是确定光子的方向，由于放射性核素的衰变各向同性特点，故而光子的方向按照各项同性处理，即方位角φ根据U[0，2π］分布，极角的余弦值μ遵循U[－1，1］分布，根据公式（3）可以得到方向矢量（uuu，vvv，www）为

然后利用网格的光子谱和光子强度信息，计算出光子的能量erg。最后对粒子赋初始权重值（通常为1.0）和寿命（通常为0.0，表示粒子刚刚产生）。至此粒子的位置（xxx，yyy，zzz）、方向（uuu，vvv，www）、能量（erg）、权重（wgt）以及寿命（tme）都已经抽样得出。

MCNP要求在任何时刻粒子所在的栅元都必须是明确的，并由icl保存该栅元的编号。因此，用户编写源子程序时，还需要确定粒子所在的几何体的编号。为确定粒子所在的栅元，源子程序会根据粒子的位置（xxx，yyy，zzz）逐个检查所有栅元，最终确定粒子所在的栅元。MCNP提供了CHECEL子程序逐个检查栅元，用户可以调用该子程序实现几何检测。该子程序的用法参考文献[5］。

3 基准测试

本文采用ITER（国际热核聚变实验堆，International Thermonuclear Experimental Reactor）组织发布的停机剂量率基准例题[9］以及ITER－T426停机剂量率实验例题[10］对源子程序进行校验。

本文使用核与辐射输运计算自动建模软件MCAM[11－14］对ITER 停机剂量率基准例题的模型的几何进行划分，共划分为1 210个网格，采用源子程序对该模型的停堆剂量率进行计算，结果如图2所示。本程序的计算结果（FDS）与英国卡拉姆聚变中心（CCFE）以及荷兰核研究与咨询研究集团（NRG）的结果相近。

图2 ITER基准测试例题模型的网格划分图Fig.2 Mesh of ITER benchmark

ITER－T426例题的几何被划分为5 000多个网格。计算程序对该例题的计算结果如图3所示，程序计算结果与实验测量结果趋势相吻合。

图3 不同程序对ITER基准例题停机剂量计算结果比较[7］Fig.3 Shutdown dose rate calculated by different codes of ITER benchmark[7］

从该图中可以看出，当冷却时间比较短时，程序计算停堆剂量率结果比实验结果最大处约低估9.4%，随着冷却时间的增加，误差稳定在3%左右。根据前文可知，计算的精确性主要依赖于中子通量和光子分布的计算结果，而其计算又依赖于网格的划分精细程度，所以源子程序计算结果的误差与网格划分的精细度程序有关。

在冷却阶段的早期，停机剂量主要由放射性核素55Mn的衰变光子产生，其半衰期为2.58h。随后58Co成为光子的主要来源，其半衰期为70.86d。理论上，剂量率的在早期阶段衰减得比较迅速，随后渐渐变缓。实验结果与理论计算都符合预期。

图4 ITER－T426网格划分Fig.4 Mesh of ITER－T426

图5 ITER－T426剂量率随时间变化图[7］Fig.5 Curves of shutdown dos rate of ITER－T426in time[7］

4 结论

针对基于网格计数的停机剂量率计算时MCNP通用源卡栅元个数受限制问题，本文采用MCNP自定义源子程序技术克服这一问题。文章系统地介绍了源子程序的编写方法，并将该方法应用于网格计数的停机剂量率计算程序中。通过国际基准的停堆剂量率实验例题对源子程序的可用性与正确性进行了验证。

[1]吴亮亮，应栋川，邱岳峰，等.三维停机剂量率计算程序研发及其在EAST上的应用[J］.核科学与工程，2011，21（1）：80－82.

[2]陈义学，吴宜灿，U.Fischer.聚变装置停机剂量率计算的严格两步（R2S）法[J］.核技术，2003，26（10），763－764.

[3]吴宜灿，李静惊，李莹，等.大型集成多功能中子学计算与分析系统VisualBUS的研究与发展[J］.核科学与工程，2007，27（4）：365－368.

[4]曾勤，邱岳峰，王国忠，等.大型集成中子学程序系统VisualBUS研发进展[A］.第四届全国反应堆物理与核材料学术研讨会论文集（C），2009.

[5]X－5Monte Carlo Team. MCNP5－A General Monte Carlo N－Particle Transport Code，Volume II：Users Guide[R］.New Mexico：Los Alamos National Library，2003.

[6]Forrest R A.FISPACT 2007：User Manual [R］.UKAEA Fusion，Report UKAEA FUS 534，March 2007.

[7]吴亮亮.基于网格计数的停机剂量率计算程序设计研究[D］.中国科学技术大学硕士论文，2012.

[8]许淑艳.蒙特卡罗方法在实验核物理中的应用[M］.北京：原子能出版社，2006.

[9]应栋川.ITER偏虑器中子学性能及其更换维修辐射剂量场分析研究[D］.中国科学院等离子体物理研究所，2011.

[10]P.Batistoni. Experimental validation of shutdown dose rates calculations inside ITER cryostat[J］.Fusion Engineering and Design，2011，58－59：613－616.

[11]吴宜灿，李莹，卢磊，等.蒙特卡罗粒子输运计算自动建模程序系统的研究与发展[J］.核科学与工程，2006，26（1）：20－27.

[12]曾勤，卢磊，李莹，等.蒙特卡罗粒子输运计算自动建模程序MCAM在ITER核分析建模中的应用[J］.原子核物理评论，2006，23（2）：138－141.

[13]王国忠，党同强，熊健，等.MCAM4.8在ITER建筑大厅中子学建模中的应用[J］.核科学与工程，2011，31（4）：352－353.

[14]吴宜灿，胡丽琴，龙鹏程，等.先进核能系统设计分析软件与数据库研发进展[J］.核科学与工程，2010，30（1）：55－64.