用于厚屏蔽小探测器的蒙特卡罗模拟减方差方法研究

杲申申，李君利，武 祯，马锐垚，王 鑫，邱 睿，李春艳，张 辉

(1.清华大学 工程物理系，北京 100084；2.粒子技术与辐射成像教育部重点实验室，北京 100084；3.同方威视技术股份有限公司，北京 100084；4.北京应用物理与计算数学研究所，北京 100094)

蒙特卡罗方法是屏蔽计算的重要工具，它可精确描述三维复杂几何与物理过程，较为容易地给出计算误差。但蒙特卡罗方法收敛速度较慢，在厚屏蔽问题及小探测器问题中统计误差较大。厚屏蔽问题的特点是源和统计区域间的屏蔽体足够厚，以至于粒子穿透屏蔽体到达统计区域的概率极小，使用常规蒙特卡罗输运统计误差较大。在小探测器问题中，探测器体积小导致粒子输运到探测器中的概率小，同样会增大蒙特卡罗统计误差。在反应堆屏蔽计算中，厚屏蔽和小探测器问题经常同时存在，常规蒙特卡罗方法难以有效解决。近年来，随着计算机性能的提升和高性能计算科学的发展，大规模并行计算可模拟海量的粒子数，在一定程度上缓解了深穿透问题带来的计算压力。同时，国内外学者在减方差技巧方面也做了大量的研究工作。在实际应用中，常用的减方差技巧有几何重要性、网格权窗、指数变换及针对小探测器的点探测器、DXTRAN方法等[1-2]，其中，确定论伴随-蒙特卡罗耦合计算方式是目前屏蔽计算减方差的主要方法之一。该方法的流程是先由确定论方法求解原问题的伴随方程，生成蒙特卡罗计算的网格权窗参数，再用蒙特卡罗程序进行计算。但该方法需额外重建几何并进行两次输运计算，且确定论程序对复杂几何的描述也不够精确[3-6]。

2004年，文献[7]针对厚屏蔽、深穿透问题提出了自动重要抽样(AIS)方法；文献[8-9]运用该方法在内照射剂量计算、探测器校正因子计算等方面获得了较好的效果；文献[10-14]从效率和实用性等方面对该方法进行了改进。AIS方法基于重要抽样、统计估计法和分层输运的思想，无需设置重要性或权窗参数，正确性和计算效率在一系列基准题上得到了验证，效率相比直接模拟可提高1～2个数量级。对厚屏蔽体外小体积探测器问题，AIS方法虽可生成虚粒子，使较多的粒子到达屏蔽体外，但因为探测器体积较小，实际抵达探测器的粒子仍然不足。而基于指向概率的DXTRAN方法[15]对小探测器的估计更为准确。因此本文基于AIS方法和指向概率方法，提出并实现了小探测器自动重要抽样(SDAIS)方法，同时针对小探测器问题，对AIS方法原有的虚粒子赌分裂算法进行改进，提出基于探测器位置的虚粒子赌分裂算法，使用NUREG/CR-6115 PWR基准题算例，对该方法应用于厚屏蔽小探测器问题的能力进行验证。

1 方法描述

1.1 AIS方法

AIS方法自2004年提出以来，经多次改进，在反应堆屏蔽、人体体素模型计算等问题的应用中取得了较好的效果。该方法基于重要性抽样和统计估计法，引入虚拟面将空间划分为多层子空间，在虚拟面上产生虚粒子向下一层子空间输运，并进行自动的粒子权重调整和数量控制。图1为AIS方法流程图，其流程为：1) 引入K个虚拟面，输运空间划分为K+1个子空间；2) 当源项抽样或粒子发生碰撞时，在虚拟面上生成虚粒子，并维持虚拟面上的虚粒子数等于源粒子数；3) 粒子若在输运中到达当前虚拟面，则将其杀死；4) 当前虚拟面上的虚粒子作为下一子空间的源粒子继续输运。

虚粒子位置为粒子沿输运方向与虚拟面的交点。权重w为碰撞粒子权重乘以碰撞粒子沿输运方向不经历碰撞直接抵达虚拟面的概率：

(1)

其中：w0为输运粒子权重；d为源或碰撞点到虚粒子位置的距离；Σt为总宏观反应截面。

图1 AIS方法流程图

1.2 MCShield程序

AIS方法由范佳锦[7]在MCNP上进行了实现，虚拟面的设置使用RDUM卡，与普通几何设置独立，需额外编写代码判断粒子是否会穿过虚拟面，因此虚拟面的类型受限，仅支持平面、球面、圆柱面等简单几何。MCShield程序[13,16]是清华大学工程物理系辐射防护实验室自主开发的中子、光子、电子耦合输运辐射屏蔽蒙特卡罗计算程序，采用CSG几何架构，支持CAD模型几何导入和计算结果三维可视化，具有强大的前、后处理能力。本文将AIS方法在MCShield程序中进行实现和改进，虚拟面的设置利用Geant4[17]的平行几何思想进行实现。平行几何是独立于真实几何的一套几何系统，只有几何形状信息，没有材料等信息。AIS方法需在粒子每次碰撞时判断与虚拟面的距离，且虚粒子需在虚拟面驻留。在平行几何中，粒子在真实几何中进行几何边界检查、驻留及碰撞，在平行几何中仅做几何边界检查和驻留。因此，平行几何能满足算法实现的需求，又可尽量减少对常规输运的效率影响。

平行几何在创建中可正常使用交、并、补及其他复杂几何结构，可与真实几何重叠，因此理论上支持任意形状的虚拟面，解决了原AIS方法虚拟面类型受限的问题，同时也增强了MCShield程序对屏蔽问题的计算能力。

1.3 SDAIS方法

对于厚屏蔽体外小探测器问题，AIS方法未利用探测器信息，产生的虚粒子并不会向探测器聚集，而是大致在虚拟面上均匀分布。当探测器体积较小时，实际抵达探测器的粒子仍然不足。在这种情况下，DXTRAN方法对小探测器的估计更为准确。因此本文将AIS方法和DXTRAN方法结合，提出SDAIS方法。该方法在源与探测器之间建立K个AIS虚拟面，将空间划分为K+1个子空间，探测器位于第K+1个，即最后1个子空间中。在最后1个子空间中，设置DXTRAN球覆盖探测器。粒子在前K层子空间输运时遵循AIS方法，在最后1个子空间中输运时，每当粒子发生碰撞，就会在DXTRAN球上产生虚粒子，使更多的粒子输运到探测器周围。图2为SDAIS方法的流程图，其具体流程为：1) 当最后1层AIS虚拟面外的粒子发生碰撞时，在DXTRAN球上产生虚粒子；2) 若普通粒子进入DXTRAN球则将其杀死；3) 当普通粒子输运完毕后，再输运DXTRAN球上的虚粒子。若粒子穿出DXTRAN球则进行轮盘赌，存活下来的粒子转化为普通粒子继续输运。

DXTRAN球上的虚粒子通过抽样产生，并被赋予相应的权重，以保证无偏性。图3为DXTRAN球上的虚粒子抽样方法。设粒子在P1发生碰撞，DXTRAN球心为P0，DXTRAN外球半径为RO，内球半径为RI，P1P0长度为L。DXTRAN外球、内球在P1P0方向上极角最大值分别为θO、θI，其余弦值分别为ηO、ηI。虚粒子位置PS在DXTRAN外球上，具体坐标通过抽样以P1P0方向为基准的极角和方位角确定。极角η余弦值按式(2)的概率密度函数p(η)抽样。

图2 SDAIS方法流程图

图3 DXTRAN球上虚粒子的抽样方法

(2)

其中：Q为调整系数，在该算法中取5；U=Q(1-ηI)+ηI-ηO。方位角按各向同性抽样。

虚粒子权重w为：

(3)

其中，p(μ)为粒子碰撞散射角余弦为μ的概率密度。

1.4 基于探测器位置的虚粒子赌分裂算法

1) 对新产生的虚粒子，计算w/r2；

该算法不仅可应用于SDAIS方法中，也可直接应用于AIS方法，作为AIS方法统计小探测器的优化选项。

2 算例验证

本文对美国NRC发布的NUREG/CR-6115 PWR注量率计算基准题进行了求解。NUREG/CR-6115 PWR基准题由美国布鲁克海文国家实验室(BNL)开发。NUREG/CR-6115 PWR基准题的反应堆总热功率为2 527.73 MW，堆芯含有204个燃耗深度不同的燃料组件。整个压水堆结构主要包括堆芯、吊篮、热屏蔽层、压力容器和混凝土生物屏蔽体等，其结构如图4所示。混凝土生物屏蔽体范围为R=335.915～549.275 cm，中子源强度为1.44×1020s-1 [18]。

a——全堆模型俯视图；b——θ=0°处的切面视图

本文选取基准题中的标准核燃料组件布置方案作为计算对象，计算了参考文献[18]中的堆腔中子剂量仪能谱，堆腔中子剂量仪位于堆腔内部R=319.56～320.56 cm、Z=176.27～178.27 cm处(图4)。为了比较计算的正确性与计算效率，分别应用MCShield几何分裂和轮盘赌(IMP-MC)方法、AIS方法及SDAIS方法进行了计算，确定论SN程序DORT的计算结果作为参考解。

同时为进一步验证SDAIS方法在厚屏蔽情况下的计算效率，本文在距堆芯更远的混凝土内部R=370、400、430、460、490、520 cm，Z=191.06 cm处各设置了直径为1 cm和5 cm的球形探测器，共12个。使用AIS方法、SDAIS方法及点探测器统计的AIS方法(F5AIS)进行了计算。同时为验证基于探测器位置的虚粒子赌分裂算法的通用性，也使用基于该算法的AIS方法和点探测器统计的AIS方法(分别记为AIS*和F5AIS*)进行了计算。

采用FOM因子代表其计算效率，则：

(4)

其中：D为相对统计误差；t为计算时间。FOM因子越大，代表计算效率越高。

3 结果分析

3.1 堆腔中子剂量仪通量(E>0.1 MeV)及能谱

IMP-MC方法的样本数为4×108，计算时间为13 736 min。AIS方法引入了4个圆柱形中子虚拟面，半径分别为188、215、230、300 cm，样本数为5×106，计算时间为206 min。SDAIS方法引入了3个半径分别为188、215、230 cm的圆柱形中子虚拟面和1个DXTRAN球。DXTRAN球的中心为剂量仪中心(317.04 cm，43.85 cm，177.27 cm)，外球半径7 cm，内球半径6.5 cm，使用了基于探测器位置的虚粒子赌分裂算法，样本数为5×106，计算时间为251 min。

统计能量0.1 MeV以上的中子通量，结果列于表1。由表1可见，以DORT程序的结果作为参考基准，几种方法结果均符合得较好，与DORT程序的相对偏差均在5%以内。AIS方法的计算效率是IMP-MC的8倍，而本文提出的SDAIS方法的计算效率约是IMP-MC方法的56倍，相比AIS方法提高了7倍。

采用与基准题中DORT方法一致的47个能量间隔，统计堆腔中子剂量仪的中子能谱。IMP-MC方法、AIS方法和SDAIS方法的中子能谱计算结果及FOM因子如图5、6所示。IMP-MC方法的中子能谱中有统计值的平均蒙特卡罗相对统计误差为26.0%，平均FOM因子为0.001 08 min-1；AIS方法的中子能谱中有统计值的平均蒙特卡罗相对统计误差是32.9%，平均FOM因子为0.044 8 min-1；SDAIS方法的中子能谱中有统计值的平均蒙特卡罗相对统计误差为7.66%，平均FOM因子为0.679 min-1。除AIS方法因统计误差较大导致与其他计算结果偏离较大，其他计算结果符合得较好。SDAIS方法的平均计算效率比AIS方法高约1个量级，比IMP-MC方法高约2个两级。

表1 堆腔中子剂量仪通量

图5 堆腔剂量仪中子能谱

图6 堆腔剂量仪中子能谱FOM因子

3.2 混凝土生物屏蔽体内球形探测器通量

在12个探测器、5种计算方法的通量结果中，除直径1 cm探测器的AIS方法和AIS*方法的统计误差较大外，其他结果的统计误差均在10%以下。直径5 cm的6个探测器的归一化通量结果如图7所示，各方法的计算结果符合得较好，各方法间的平均偏差在10%以内。

图7 生物屏蔽体内5 cm球形探测器中子通量

图8 生物屏蔽体内球形探测器FOM因子比较

对不同位置探测器的FOM因子进行平均，结果如图8所示。对于探测器直径5 cm的情形，AIS*方法效率最高，SDAIS方法次之；对于探测器直径1 cm的情形，SDAIS方法效率最高，传统AIS方法最差。这是因为探测器的直径与直接进入探测器区域的粒子数呈正相关，而后者直接决定了统计误差。当探测器直径较小时，直接进入探测器区域的粒子数很少，此时需使用指向概率的方法进行统计。此外可看到，使用基于探测器位置的虚粒子赌分裂算法的计算效率比传统方法高约1个量级。另外，SDAIS方法和F5AIS*方法均使用了指向概率方法，但SDAIS方法比F5AIS*方法的计算效率高些，这是因为SDAIS方法仅在最后1个虚拟面外才进行指向概率的计算，计算量较小。

4 结论

本文在AIS方法的基础上，针对小探测器问题，优化了虚粒子赌分裂算法，将AIS方法和指向概率方法相结合，提出了SDAIS方法，在自主开发的蒙特卡罗程序MCShield中进行了实现。使用NUREG/CR-6115 PWR基准题对该方法的正确性及计算效率进行了验证。计算结果显示，当探测器较小时，SDAIS方法的计算效率比传统AIS方法高约1个量级，比常规蒙特卡罗方法高约2个量级。基于探测器位置的虚粒子赌分裂算法可较为有效地提升传统AIS方法在小探测器问题中的计算效率。