蒙特卡罗粒子输运权窗减方差方法研究及在ITER屏蔽分析中应用

郑华庆，宋 婧，郝丽娟，孙光耀，胡丽琴，FDS团队

（中国科学院核能安全技术研究所，安徽 合肥230031）

蒙特卡罗方法具备处理复杂几何和模拟精细的物理过程的能力，方法本身和程序结构的简单性，因而在解决粒子输运问题方面应用越来越广泛。屏蔽分析是聚变装置物理设计与分析的关键环节，但蒙特卡罗方法在处理屏蔽问题时，方法本身收敛速度比较慢是一个关键的问题，因此往往需要很长的计算时间，在深穿透情况下，这个问题尤其突出。减方差方法可以加快蒙特卡罗方法的收敛速度，因此在应用蒙特卡罗方法解决屏蔽计算问题的过程中，有效地使用减方差方法是必不可少的环节。

现有的蒙特卡罗粒子输运程序中有四类减方差方法［1］：截断方法（能量截断和时间截断）、总体控制方法（几何分裂与轮盘赌，能量分裂与轮盘赌，时间分裂与轮盘赌，权截断和权窗）、修正抽样方法（指数变换，隐俘获，强迫碰撞，源偏倚）和部分确定性方法（点探测器，DXTRAN和相关抽样）。

Booth和 Hendrick提出的权窗方法［2］，在多年的实际应用中逐渐被证实是目前最为有效和最为通用的减方差方法之一。

超级蒙特卡罗计算软件SuperMC［3］是由FDS团队自主研发，采用蒙特卡罗方法并耦合确定论方法，定位于基于先进计算机技术实现输运（稳态粒子输运、中子动力学）、燃耗与活化（同位素燃耗、材料活化、停机剂量）、多物理耦合（热工水力学、燃料性能、结构力学）过程的高效能精细模拟，集建模、计算、可视化分析于一体，可广泛应用于反应堆物理、辐射屏蔽、医学物理、核探测、高能物理等领域。

本文在中国科学院核能安全技术研究所·FDS团队前期研究工作的基础［4－12］上，将权窗减方差方法在SuperMC中实现，并通过基准例题的校验和国际热核聚变实验堆ITER屏蔽分析的应用，证明了方法和程序的正确性和有效性。

1 权窗减方差方法的设计与实现

权窗是相空间（空间－能量、空间－时间）上的分裂与轮盘赌，本文将权窗减方差方法在SuperMC实现的方法如下：

用户首先给出权窗基本参数CU（权窗上限参数）、CS（赌存活参数）、MXSPIN（最大分裂参数／最大赌存活参数）。一般默认的参数设置是：CU＝5，CS＝3，MXSPIN＝5。

对于每一个相空间栅元有一个权重下限WL（用户给出）和上限 WU（WU＝WL×CU），当粒子的权重WGT小于下限WL时，发生轮盘赌：赌赢则粒子权重变为 min（CS×WL，WGT×MXSPIN），赌输则该粒子被杀死；当粒子的权重WGT大于上限WU时，发生分裂：分裂出npa（npa＝min（int（WGT／WU ＋1），MXSPIN）个权重为 WGT／npa的全同粒子；当粒子权重WGT在WU和WL之间时，粒子不发生任何变化。低权重粒子发生轮盘赌，使得不会在权重无意义的粒子上浪费时间；权重高于权窗，粒子发生分裂能有效避免极高权重的粒子对计数的扰动，其中流程如下图所示：

图1 权窗减方差流程图Fig.1 Flow chart of weight windows variance reduction

在SuperMC的用户图形界面（GUI）中，如图2所示，用户可以通过可视化几何交互的方式进行权窗参数的设置。在处理复杂的几何模型时，根据放射源和计数区的位置，用户可以方便地点击每个栅元，设置相应栅元的权窗参数。其设置的原则是：粒子由放射源所在栅元开始到计数区为止，考虑其可能经过的所有的栅元，在这些栅元上根据粒子通过该栅元的可能性来设置这些栅元的权窗下限，保证粒子通过合理的引导进入计数区。对于其他对计数区贡献可以忽略的栅元，可以将权窗下限设为0，让进入该区域的粒子进行权截断。

图2 SuperMC几何建模功能用户界面Fig.2 Interface of geometry modeling of SuperMC

2 程序校验及在ITER屏蔽分析中应用

2.1 混凝土基准校验例题

2.1.1 模型描述

混凝土是聚变堆典型的屏蔽体材料。计算例题选自参考文献［13］，其基本的信息如下：200cm厚的混凝土屏蔽层，材料密度为ρ＝2.03g／cm3，被分成20层，每层是10cm，如图3所示。放射源是位于（X＝0cm；Y＝0cm；Z＝1.0×10－6cm），单能（14MeV）、单方向（0，0，1）的中子点源。计数栅元是混凝土圆柱最上层的栅元，计数类型为体通量计算。三组权窗的栅元权重下限设置参见表1。

图3 混凝土基准测试例题几何示意图Fig.3 Geometry of concrete benchmark case

表1 不同权窗的栅元权重下限设置Table 1 Weight low limit in each cell of different weight windows

2.1.2 测试结果

SuperMC与MCNP使用相同的数据库HENDL［14］，计算结果如表2所示，在没有添加权窗、添加权窗1、添加权窗2和添加权窗3，这三种计算条件下，SuperMC通量计算结果与参考程序MCNP的计算结果基本吻合，偏差均小于0.5%。

计算效率的比较如表3所示，在没有添加权窗的情况下，计算品质因子FOM（Figure of Merit）仅为0.593，添加权窗1的情况FOM为4.109，计算效率提升了近6倍；添加权窗2的情况FOM为66.077，计算效率提升了110倍；添加权窗3的情况FOM为90.240，计算效率提升了150倍。

表2 混凝土测试例题通量计算结果Table 2 Flux results of concrete benchmark case

表3 混凝土测试例题计算效率比较Table 3 Efficiency comparison of FOM of concrete benchmark case

2.2 ITER－T426屏蔽分析应用例题

2.2.1 模型描述

国际热核聚变实验堆ITER是目前全球规模最大的国际科研合作项目之一，ITER装置是一个能产生大规模核聚变反应的超导托卡马克。本文选用ITER－T426例题［15］作为应用例题。该例题基于 FNG（Frascati neutron generator）装置，设计一个由不锈钢和有机玻璃组成的屏蔽结构，模拟ITER真空室。本文选用ITER－T426例题，如图4所示。源是位于底层 栅 元 （X＝0cm；Y＝ 0.001cm；Z＝0cm）处，单能（14MeV）、单方向（0，1，0）的中子点源。计数区域是中心点位于（X＝0cm；Y＝74.93cm；Z＝0cm），半径为1.5cm，高度为5cm的不锈钢圆柱，计数类型为体通量计算。

图4 ITER－T426例题在SuperMC的显示图Fig.4 Geometry view of ITER－T426in SuperMC

2.2.2 测试结果

SuperMC与MCNP使用相同的核数据库HENDL，计算结果如表4所示，在没有添加权窗、添加权窗两种计算条件下，SuperMC通量计算结果与参考程序MCNP的结果基本吻合，偏差均小于0.2%。

计算效率的比较如表5所示，在没有添加权窗的情况下，计算品质因子FOM仅为59.707，添加权窗的情况FOM 为905.199，计算效率提升了14倍。

表4 ITER－T426例题计算结果Table 4 Flux results of ITER－T426case

表5 ITER－T426例题计算效率比较Table 5 Efficiency comparison of FOM of ITER－T426case

3 总结

本文对蒙特卡罗粒子输运模拟中权窗减方差方法进行了研究：给出了关键参数设置的原则和程序实现的思路，并基于超级蒙特卡罗计算软件SuperMC进行实现。通过混凝土基准例题的校验和ITER屏蔽分析的应用，证明了方法和程序的正确性和有效性，通过合理设置权窗参数可以提高计算效率。

［1］ X－5Monte Carlo Team，MCNP－A General Monte Carlo N－Particle Transport Code，Version 5［R］.Los Alamos National Laboratory Report LA－UR－03－1987，April 24，2003.

［2］ Booth T E，Hendricks J S.Importance Estimation in Forward Monte Carlo Calculations ［J］. Nuclear Technology Fusion，1984，5：90－100.

［3］ 孙光耀，宋婧，郑华庆，等 .超级蒙特卡罗软件SuperMC2.0中子输运计算校验［J］.原子能科学技术，2013，47（s2）：520－525.

［4］ 吴宜灿，胡丽琴，龙鹏程，等 .先进核能软件发展与核信息学实践［J］.中国科研信息化蓝皮书2013，2013.12：232－244.

［5］ Li Y，Lu L，Ding A P，et al.Benchmarking of MCAM 4.0with the ITER 3DModel［J］.Fusion Engineering and Design，2007，82（15－24）：2861－2866.

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

［7］ Wu Y C，FDS Team.CAD－based Interface Programs for Fusion Neutron Transport Simulation ［J］.Fusion Engineering and Design，2009，84（7－11）：1987－1992.

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

［9］ Wu Y C，Xie Z S，Fischer U.A Discrete Ordinates Nodal Method for One－dimensional Neutron Transport Calculation in Curvilinear Geometries［J］.Nuclear Science and Engineering，1999，133（3）：350－357.

［10］ Wu Y C，FDS Team.Conceptual Design Activities of FDS Series Fusion Power Plants in China［J］.Fusion Engineering and Design，2006，81（23－24）：2713－2718.

［11］ Wu Y C，FDS Team.Design Analysis of the China Dual－functional Lithium Lead （DFLL）Test Blanket Module in ITER［J］.Fusion Engineering and Design，2007，82（15－24）：1893－1903.

［12］ Wu Y C，FDS Team.Overview of Liquid Lithium Lead Breeder Blanket Program in China ［J］. Fusion Engineering and Design，2011，86（9－11）：2343－2346.

［13］ Booth T E，MCNP Variance Reduction Examples［R］.Los Alamos National Laboratory，Diagnostic Applications Group X－5，Mail Stop F663，December 23，2004.

［14］ Zou J，He Z Z，Zeng Q，et al.Development and Testing of Multigroup Library with Correction of Self－shielding Effects in Fusion－Fission Hybrid Reactor［J］.Fusion Engineering and Design，2010，85（7－9）：1587－1590.

［15］ Batistoni P.Experimental Validation of Shutdown Dose Rates Calculations Inside ITER Cryostat［J］.Fusion Engineering and Design，2001，58－59：613－616.