一种基于蒙特卡罗自适应降低方差的辐射屏蔽计算方法

摘要

本发明公开了一种基于蒙特卡罗自适应降低方差的辐射屏蔽计算方法，在传统的基于权窗减方差的蒙特卡罗方法的基础上，根据其确定论方法伴随通量计算结果的具体情况，进而计算得到重要性参数，将这些参数写入基于蒙特卡罗的辐射屏蔽计算程序的输入文件中，再进行辐射屏蔽计算。在基于蒙特卡罗的辐射屏蔽计算过程中，在根据前一步得到的重要性参数对粒子进行偏倚输运的同时，对粒子径迹进行保存，即给出每一个粒子经过的栅元，帮助用户在模拟过程中自动干预调整几何区域重要性参数的设置，并在后续计算中使用新的栅元重要性参数，实现自动偏倚。本发明避免了传统的基于经验式权窗减方差的蒙特卡罗方法的使用，达到精确快速进行辐射屏蔽计算的效果。

说明书

技术领域

本发明涉及一种基于蒙特卡罗自适应降低方差的辐射屏蔽计算方法，属于核物理、核技 术应用等领域中的辐射屏蔽计算分析方向。

背景技术

辐射屏蔽分析是核装置物理设计与分析的关键环节。蒙特卡罗方法具备处理复杂几何和 模拟精细的物理过程的能力，因而在屏蔽分析方面应用越来越广泛。对于厚屏蔽与结构复杂 问题，由于模拟的粒子难以穿透介质继续输运，难以统计到有效的粒子，不仅计算效率难以 接受且难以得到收敛的解。减方差方法可以加快蒙特卡罗方法的收敛速度，因此在应用蒙特 卡罗方法解决屏蔽计算问题的过程中，有效地使用减方差方法是必不可少的。但是传统减方 差方法需要使用者具备丰富的经验和反复的试算，才能得到比较好的效果。如何避免经验式 的减方差方法使用，如何更高效的发挥减方差方法的使用，是基于蒙特卡罗辐射屏蔽计算方 法的关键技术之一。

减方差不仅仅是蒙特卡罗粒子输运模拟问题加速的一种技巧，并且是得到答案所完全必 须的条件。比如深穿透问题和管探测器问题，若不使用足够的减方差技巧，将会因为许多因 素而运行的很慢。因此，这就要求用户必须很熟练的使用蒙特卡罗粒子输运程序中的减方差 技巧。

现有的蒙特卡罗粒子输运程序中有四类减方差技术：截断技术(能量截断和时间截断)、 总体控制方法(几何分裂与轮盘赌，能量分裂与轮盘赌，时间分裂与轮盘赌，权截断和权窗)、 修正抽样方法(指数变换，隐俘获，强迫碰撞，源偏倚)、部分确定性方法(点探测器，DXTRAN 和相关抽样)。其中权窗是应用最为广泛的技术之一。

权窗是相空间的分裂和俄罗斯轮盘赌技巧，相空间可以是空间-能量、空间-时间或者空 间。作为能量、几何位置或者两者的一部分，低权重粒子通过俄罗斯轮盘赌淘汰，高权重粒 子分裂。这项技术可以使得权重分布保持在一个合理的界限内。

传统的权窗减方差方法在使用的过程中，虽然实践证明它的有效性和可用性，但它决不 是解决所有重要抽样问题的万金油。多数情况下，用户还要判断计算出的权窗是否合理。用 户需牢记，产生的参数是统计结果，有可能带很大误差，除非在设置的相空间中已充分抽样。 为此，用户应在使用权窗时预先给相空间粗略设置比较合理的重要性，经过计算-设置-再计 算的多次过程才能得到较好的权窗；要对计算结果加以判断，对那些值得怀疑的参数重新设 置。对于使用经验不足的用户，反复调试权窗参数极为耗时。

发明内容

本发明的目的在于：克服传统权窗减方差方法在操作中的不足(权窗减方差方法在操作 中的不足是不行的，没有与技术领域相结合，提供一种自适应降低方差的方法，使得用于基 于蒙特卡罗的辐射屏蔽计算的粒子自动从源区域向计数区域进行偏倚，避免了传统的基于经 验式权窗减方差的蒙特卡罗方法的使用，达到精确快速进行辐射屏蔽计算的效果。

本发明的技术方案如下：

本发明根据确定论方法伴随通量计算结果的具体情况，进而计算得到蒙特卡罗正向模拟 的重要性参数(源偏倚和权窗)；在正向粒子输运模拟过程中，在根据上一步得到的重要性 参数对粒子进行偏倚输运的同时，对粒子径迹进行保存，即给出每一个粒子经过的栅元(包 括目标计算栅元，判断粒子是否对目标计数有贡献)，帮助用户在模拟过程中自动干预调整 几何区域重要性参数的设置，并在后续计算中使用新的栅元重要性参数，在粒子输运物理过 程实现自动偏倚。

本发明所采用的技术方案实现如下：一种基于蒙特卡罗自适应降低方差的辐射屏蔽计算 方法，包含如下步骤：

(1)基于蒙特卡罗的辐射屏蔽计算方法中初始减方差参数获取：

a)根据需计算的基于蒙特卡罗的辐射屏蔽计算问题，建立其确定论方法伴随计算的计 算模型；$\frac{1}{\upsilon }\frac{\partial \Phi }{\partial t}=\int \int \Sigma (\overrightarrow{r},E^{\prime })f(r,E^{\prime }\to E,{\overrightarrow{\Omega }}^{\prime }\to \overrightarrow{\Omega })\Phi (\overrightarrow{r},E^{\prime },{\overrightarrow{\Omega }}^{\prime },t)d{\overrightarrow{\Omega }}^{\prime }d{E}^{\prime }+Q+S-\overrightarrow{\Omega }·▿\Phi -\Sigma (\overrightarrow{r},E)\Phi$

其中：

Φ表示为中子角通量；

为位置处能量为E的中子与核发生反应的宏观截面；

为位置处能量为E′的中子与核发生反应的宏观截面；

为处发生某种类型的核反应后，中子由能量为E′和方向为 变为能量为E方向为在能量为dE′范围内方向为范围内的几率；

Q表示为单位时间、单位体积和单位立体角内在处t时刻，一个能量 为E方向为的中子出现的几率；

S表示为独立中子源。

b)根据确定论方法粒子输运伴随计算的结果，生成基于蒙特卡罗的辐射屏蔽计算的减 方差参数，包括：源偏倚参数和权窗参数；

源偏倚参数：

其中，S*(p)为正向计算时源在位置p处的源偏倚参数，S(p)为正向计算时源在位 置p处的初始重要性，为位置p处的伴随通量，R为伴随计算的响应值，

权窗参数：

其中，E是模拟粒子的能量，r是模拟粒子的位置，W_L(r,E)是模拟粒子所在的能量E， 位置r处的权窗下限，S为正向计算源强，C_u是权窗上限与权窗下限的比值，一般为5，W 是模拟粒子的权重，为粒子所在的能量E，位置r处的的伴随通量，R为伴随计算 的响应值；

(2)基于蒙特卡罗的辐射屏蔽计算中粒子径迹跟踪及减方差参数调节：

a)分批进行基于蒙特卡罗的辐射屏蔽计算，首批用于辐射屏蔽计算的粒子模拟使用步 骤(1)获得的减方差参数；

$F=\frac{{\mathrm{WT}}_{l}E}{V}$

其中，F是模拟粒子在计数栅元的通量，W是模拟粒子在计数栅元的权重，T_l是模拟粒 子在计数栅元的径迹长度，E是模拟粒子在计数栅元的能量，V是模拟粒子在计数栅元的体 积；

b)在粒子模拟的过程中，对每一个粒子的径迹进行跟踪，保留对目标计数有贡献的粒 子径迹信息，包括：粒子穿过栅元的编号，粒子在目标计数栅元中的权重；

c)在该批粒子模拟结束后，根据公式计算每个栅元的重要性，并按照重要性的倒数，重新定义各栅元的权窗参数，用于下批粒子 模拟；

重复步骤(2)中的b、c部分，直至所有的分批粒子模拟结束，最后统计得到辐射屏蔽 计算的结果。

$F=\frac{{\mathrm{WT}}_{l}E}{V}$

其中，F是模拟粒子在计数栅元的通量，W是模拟粒子在计数栅元的权重，T_l是模拟粒 子在计数栅元的径迹长度，E是模拟粒子在计数栅元的能量，V是模拟粒子在计数栅元的体 积；

步骤(2)中的c)部分计算栅元重要性时，根据公式 计算时，当计算得到的栅元重要性小于0.01 时，则定义栅元的权窗参数为0；当计算得到的栅元重要性小于0.0001时，则定义栅元的权 窗参数为-1，减小对目标计数栅元贡献较小的粒子输运的时间。

本发明与现有技术具有如下优点：

(1)根据伴随计算结果自适应减小蒙特卡罗模拟的计算方差，无需用户反复调试减方 差的参数；避免了传统的基于经验式权窗减方差的蒙特卡罗方法的使用，达到精确快速进行 辐射屏蔽计算的效果。

(2)在蒙特卡罗模拟计算中根据粒子标识，对减方差参数进行二次优化，进一步提高 计算的效率；

(3)当计算得到的栅元重要性小于0.01时，则定义栅元的权窗参数为0；当计算得到 的栅元重要性小于0.0001时，则定义栅元的权窗参数为-1，减小对目标计数栅元贡献较小的 粒子输运的时间，达到进一步精确快速进行辐射屏蔽计算的效果。

附图说明

图1为本发明的实现流程图。

具体实施方式

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

(1)蒙特卡罗粒子输运正向计算初始减方差参数获取：

a)根据需计算ITER-T426例题的蒙特卡罗粒子程序(例如：MCNP等)输运正向计算 问题，建立确定论方法粒子输运伴随计算的计算模型(例如：三维离散纵标程序TORT的输 入文件)，补充TORT输入文件中材料信息、放射源信息、计数信息，调用TORT进行伴随 计算。需要说明的是，在补充材料信息之前，需要使用GIP(橡树岭国家实验室开发的多群 截面处理程序)生成TORT伴随计算的截面数据。另外，伴随计算的放射源的空间分布对应 正向计算时的计数区域，伴随计算的能谱为响应能谱。而伴随计算的计数区的空间分布对应 正向计算时的放射源区域，能谱对应正向计算时的放射源能群；

b)根据ITER-T426例题确定论方法粒子输运伴随计算的结果根据TORT伴随 计算结果，得到不同网格、不同能群的伴随通量，将与该位置、该能群的正向计算实际源强 的乘积求和，就得到了响应量R。根据以下公式计算MCNP正向计算时源偏倚参数和权窗参 数：

其中，S为正向计算源强，为伴随通量，C_u一般为5。E是模拟粒子的能量，r是模拟 粒子的位置，W是模拟粒子的权重，W_L是模拟粒子所在的能量E，位置r处的权窗下限。

(2)使用蒙特卡罗粒子输运模拟计算程序(例如：MCNP等)进行蒙特卡罗粒子输运 正向计算，在计算的过程中对粒子径迹跟踪及减方差参数调节：

a)分批进行蒙特卡罗粒子输运正向计算，首批粒子模拟使用步骤(1)获得的减方差参 数；

b)在粒子模拟的过程中，对每一个粒子NPS的径迹进行跟踪，保留对目标计数有贡献 的粒子径迹信息，包括：粒子穿过栅元的编号CELLnum，粒子在目标计数栅元中的权重W_NPS；

b)在该批粒子模拟结束后，根据公式计算每个栅元的重要性，并按照重要性的倒数，重新定义各栅元的权窗参数W_L，用于下批 粒子模拟；

d)重复步骤(2)中的b、c部分，直至所有的分批粒子模拟结束。

(3)通过比较含有自适应减方差方法的蒙特卡罗正向模拟计算结果和直接蒙特卡罗正 向模拟的计算结果，自适应减方差方法达到加快计算收敛，避免了传统的基于经验式权窗减 方差的蒙特卡罗方法的使用，达到精确快速进行辐射屏蔽计算的效果。