一种基于多群虚截面的连续能点蒙特卡罗粒子输运方法

摘要

本发明涉及一种基于多群虚截面的连续能点蒙特卡罗粒子输运方法，建立多群统一能量网格，并预计算材料与核素对应于多群统一能量网格的截面数据，使得每步输运中只需要对多群统一能量网格进行一次搜索查找即可得到输运计算中需要的各种截面数据，从而极大的地减少计算时间；同时，正常核素连续能量网格具有上万或者更多的能点，常用的统一能量网格方法中网格点更是高达几十万甚至上百万，多群统一能量网格中网格点非常少，因此在大大减少内存的同时，由于单次搜索时间地减少，也节约了计算时间，保证计算精度不受到损失。

说明书

技术领域

本发明涉及一种基于多群虚截面的连续能点蒙特卡罗粒子输运方法，属于原子能科学技术领域。

背景技术

蒙特卡罗(MC)方法是目前粒子输运模拟中较常用的方法。MC粒子输运模拟从连续能量点截面的核数据中获取粒子与物质发生各种反应的相关截面数据用于指导输运过程中的各种抽样过程。对于不同核素连续能量点的能量网格划分各不相同，在MC计算中，需要频繁的根据入射能量对每个核素的能量网格进行搜索找到对应的网格索引，以得到相应的截面。这些频繁的网格搜索导致核素多时计算很慢，在输运燃耗耦合计算中尤其明显。

蒙特卡罗软件Serpent与OpenMC等提出了统一能量网格的思想，减少粒子输运计算中网格搜索的次数。其基本思想是：构造所用核素适用的能量网格，计算时，只需对该能量网格进行一次搜索，即可得到当前能量下所有核素的截面数据，避免了因各核素能量网格不同而进行的重复的能量网格搜索。但常用的统一能量网格方法中(如专利《一种蒙特卡罗粒子输运模拟中核截面数据处理优化方法》与文献《The OpenMC Monte Carlo particle transport code》)，通常将问题涉及的所有核素的能量网格进行合并，去除重复能点后直接作为统一能量网格使用，这种处理方法实现简单，能完整保留原始核数据的精度，但是由于单核素的能量网格点数一般为几千甚至几万、几十万，这种方法产生的统一能量网格点将达到几百万甚至几千万的级别，导致对此网格进行搜索比较耗时，更为严重的是，某些方法中为了进一步提高计算速度，针对统一能量网格给每个核素建立新的截面数组或者指针数组，这会非常大的增加内存开销。

文献《Two practical methods for unionized energy grid construction in continuous-energy Monte Carlo neutron transport calculation》中发展了一种优化的统一能量网格方法，对每个核素的能量网格进行分析，只将截面数据改变趋势发生改变的能点与一些反应(非弹性散射、裂变等)的阈值能量挑选出来，去重后，形成统一能量网格。这种方法中统一能量网格中能点会减少很多(50％以上)，但是仍然比较多，而且处理复杂，很容易误操作导致核数据精度的丢失。

发明内容

本发明要解决的技术问题为：本发明的目的在于改善蒙特卡罗粒子输运计算方法在核截面数据处理上的不足，提供了一种基于多群虚截面的连续能点蒙特卡罗粒子输运方法，采用多群统一能量网格的方法提高计算速度降低内存消耗，并采用多群虚截面拒绝抽样的方式保证计算精度。

本发明的技术方案如下：一种基于多群虚截面的连续能点蒙特卡罗粒子输运方法，包括以下步骤：

步骤(1)对蒙特卡罗粒子输运模拟中从核数据库中读取的连续能点核数据进行预处理，包括以下内容：

(11)多群统一能量网格数组采用用户给定的或者程序内置的多群统一能量网格数组，记为MultiGErgArray[N_MG]，其中N_MG表示多群统一能量网格数组中的能量网格点个数(N_MG建议为50-400)；

(12)对比每个核素i的能量网格数组EryArray_i[N_i]，N_i表示核素i的能量网格点个数，与步骤(11)中得到的多群统一能量网格数组MultiGErgArray[N_MG]比较，找出多群统一能量网格数组中的每个能量网格点MultiGErgArray[j]在核素i的能量网格数组EryArray_i[N_i]中的位置p，p满足：

EryArray_i[p]＜MultiGErgArray[j]＜EryArray_i[p+1]，

标记p为Position_i[j]，存入核素i的多群指针数组，标记为Position_i[N_MG]，建立起每个核素的核素多群指针数组；

(13)每个核素i对应多群统一能量网格数组中的能量网格点MultiGErgArray[j]的多群虚总截面标记为σ_MG,tot,i[j](MG表示多群，tot表示总截面)，按照如下方式获得：

σ_MG,tot,i[j]＝max(σ_tot,i[k],k＝Position_i[j],Position_i[j+1]+1)，

其中σ_tot，i[k]是单核素i的能量网格数组EryArray_i[N_i]中的能量网格点EryArray_i[k]对应的核素i的总截面值，将σ_MG,tot,i[j]存入核素i的多群虚总截面数组，标记为σ_MG,tot,i[N_MG]，从而建立起每个核素的多群虚总截面数组；

(14)在步骤(13)中建立核素的多群虚总截面数组的基础上，将材料k中所有核素的多群虚总截面值按密度比例求和得到材料的多群虚宏观总截面数组，标记为Σ_MG,k[N_MG]，从而建立起每个材料的多群虚总截面数组；

步骤(2)依据步骤(1)中预处理得到的每个核素的核素多群指针数组、每个核素的多群虚总截面数组和每个材料的多群虚总截面数组，进行基于多群虚截面的连续能点蒙特卡罗粒子输运模拟，包括以下步骤：

(21)使用二分查找法搜索当前粒子能量erg在步骤(11)中产生的多群统一能量网格数组MultiGErgArray[N_MG]中的位置，记为index_MG，(index_MG满足MultiGErgArray[index_MG]≤erg＜MultiGErgArray[index_MG+1])，当前粒子能量下当前材料k的多群虚宏观总截面值Σ_MG，按照如下方式获得：

Σ_MG＝max(Σ_MG,k[index_MG]，Σ_MG,k[index_MG+1])；

当前粒子能量下当前材料k中的每个核素i的多群虚总截面值σ_MG,tot,i，按照如下方式获得：

σ_MG,tot,i＝max(σ_MG,tot,i[index_MG]，σ_MG,tot,i[index_MG+1])；

(22)使用步骤(21)中得到的当前粒子能量下当前材料k的多群虚宏观总截面代替传统的当前材料真实宏观总截面进行输运长度抽样，并根据抽样得到的输运长度、粒子的原始位置和粒子的原始运动方向，更新粒子的位置；

(23)使用步骤(21)中得到的当前粒子能量下当前材料中各核素的多群虚总截面与当前材料k的多群虚宏观总截面的比值，代替传统方法中的材料中各核素的真实总截面与当前材料真实宏观总截面的比值，进行反应核素的抽样，得到反应核素为i_col；

(24)从核数据库中读取的连续能点核数据计算反应核素i_col的真实总截面，除以核素i_col的多群虚总截面σ_MG,tot,i的，得到拒绝比率ρ，抽取一个随机数ξ，如果ξ＞ρ，则反应被拒绝，为虚反应，粒子在此处不发生反应，转至步骤(22)继续抽样输运长度进行下一步输运；如果ξ≤ρ，反应是真反应，粒子在此处真实发生反应，按照传统蒙卡粒子输运方法，根据从核数据库中读取的反应核素的连续能点核数据得到真实弹性散射、非弹性散射、吸收反应截面，进行反应类型抽样、反应后次级粒子产生和粒子状态的抽样，更新粒子的能量、出射角度信息，完成粒子输运的一步，转至步骤(21)进行粒子输运的下一步，输运停止的条件与传统蒙特卡罗粒子输运方法中相同，从而实现基于多群虚截面的连续能点蒙特卡罗粒子输运模拟。

核数据可以从国际原子能机构(IAEA)的官网下载获得；

蒙特卡罗粒子输运程序可以使用发展的SuperMC软件或者开源的Geant4软件。

本发明与现有技术相比的优点在于：本发明的一种基于多群虚截面的连续能点蒙特卡罗粒子输运方法，与传统蒙特卡罗粒子输运方法相比，通过建立多群统一能量网格，每步输运计算中只对多群统一能量网格进行一次搜索查找，替代了传统方法对材料中每个核素的不同能量网格的重复搜索查找，极大地减少了蒙特卡罗粒子输运计算中能量网格的搜索查找次数，极大地提高计算速度；与现有统一能量网格方法相比，多群统一能量网格的网格点数非常少，单次搜索查找的时间显著减少，提高了截面数据的计算速度，同时，多群统一能量网格的内存开销远远低于现有统一能量网格方法，而且采用多群虚截面拒绝抽样的输运模式避免了多群统一能量网格导致的计算精度损失的问题，实现了在保证计算精度的情况下，消耗更少的内存，得到更好的蒙特卡罗粒子输运计算速度。

附图说明

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

具体实施方式

如图1所示，本发明的实现过程，包括以下内容：

1、对蒙特卡罗粒子输运模拟中从核数据库中读取的连续能点核数据进行预处理，通过如下步骤实现：

第一步，多群统一能量网格数组采用用户给定的或者程序内置的多群统一能量网格数组，记为MultiGErgArray[N_MG]，其中N_MG表示多群统一能量网格数组中的网格数(N_MG建议为50-400)；

第二步，多群核素指针数组的建立：对比每个单核素i的能量网格数组EryArray_i[N_i]，N_i表示网格数，与多群统一能量网格数组MultiGErgArray[N_MG]，找出多群统一能量网格数组中的每个能点MultiGErgArray[j]在单核素i能量网格EryArray_i[N_i]中的位置p，p满足：

EryArray_i[p]＜MultiGErgArray[j]＜EryArray_i[p+1]，

标记p为Position_i[j]，存入核素i的多群指针数组，标记为Position_i[N_MG]；

第三步，核素的多群虚总截面数组的建立：每个核素i对应多群统一能量网格数组中的能点MultiGErgArray[j]的多群虚总截面标记为σ_MG,tot,i[j](MG表示多群，tot表示总截面)，按照如下方式获得：

σ_MG,tot,i[j]＝max(σ_tot,i[k],k＝Position_i[j],Position_i[j+1]+1)，

其中σ_tot，i[k]是单核素i的能量网格EryArray_i[N_i]中的能点EryArray_i[k]对应的核素i的总截面值，将σ_MG,tot,i[j]存入单核素i的多群虚总截面数组，标记为σ_MG,tot,i[N_MG]；

第四步，材料的多群虚宏观总截面数组的建立：基于核素的多群虚总截面数组，将材料k中所有核素的多群虚总截面值按密度比例求和得到材料的多群虚宏观总截面数组，标记为Σ_MG,k[N_MG]；

2、基于多群虚截面的蒙特卡罗粒子输运模拟，通过如下步骤实现：

第一步，当前粒子能量下当前材料k的多群虚宏观总截面以及材料中各核素的多群虚总截面计算：使用二分查找法搜索当前粒子能量erg在多群能量网格数组MultiGErgArray[N_MG]中的位置，记为index_MG，(index_MG满足MultiGErgArray[index_MG]≤erg＜MultiGErgArray[index_MG+1])，当前粒子能量下当前材料k的多群虚宏观总截面值Σ_MG，按照如下方式获得：

Σ_MG＝max(Σ_MG,k[index_MG]，Σ_MG,k[index_MG+1])；

当前粒子能量下核素i的多群虚总截面值σ_MG,tot,i，按照如下方式获得：

σ_MG,tot,i＝max(σ_MG,tot,i[index_MG]，σ_MG,tot,i[index_MG+1])；

第二步，粒子输运长度抽样：使用当前粒子能量下当前材料k的多群虚宏观总截面代替传统的当前材料真实宏观总截面进行输运长度抽样，并根据输运长度、粒子的原始位置和粒子的运动方向，更新粒子的位置；

第三步，反应核素抽样：使用当前粒子能量下材料中各核素的多群虚总截面与当前材料的多群虚宏观总截面的比值，代替传统的材料中各核素的真实总截面与当前材料真实宏观总截面的比值，进行反应核素的抽样，得到反应核素为i_col；

第四步，基于拒绝抽样的虚反应判断：计算反应核素i_col的真实总截面，以及与多群虚总截面σ_MG,tot,i的比值ρ，抽取一个随机数ξ，如果ξ＞ρ，则反应被拒绝，为虚反应，粒子在此处不发生反应，转至步骤2的第二步继续抽样输运长度进行下一步输运；如果ξ≤ρ，反应是真反应，粒子在此处真实发生反应，根据反应核素的真实弹性散射、非弹性散射、吸收反应截面，进行反应类型抽样，以及反应后次级粒子产生和粒子状态的抽样。更新粒子的能量、出射角度信息，完成粒子输运的一步，转至步骤2中的第一步进行粒子输运的下一步，输运停止的条件与传统蒙特卡罗粒子输运方法中相同，从而实现基于多群虚截面的连续能点蒙特卡罗粒子输运模拟。

3.计算实例，包含以下内容：

计算实例是一个半径为20cm的球体，在球心处有一个点源，能量为14MeV，方向服从各向同性缝补，球体中填充的材料为97％的U-233、2％的U-234与1％的U-235的混合物，需要计算球面的面积分流量：

从数据库中读取据核素U-233、U-234、U-235的能量网格数组，分别为25036、27637、77165个数据点，多群统一能量网格数组使用程序内置的175群能量网格。

对比U-233的能量网格数组与多群统一能量网格数组，得到多群统一能量网格数组中的每一个点在U-233的能量网格数组中的位置，即核素指针数组。U-233对应多群统一能量网格数组中每一个能点的虚总截面取原U-233在此能量区域中总截面的最大值；同样对U-234、U-235建立核素指针数组与多群虚总截面数组。

对多群统一能量网格数组中的每一个能量点将U-233、U-234、U-235的虚总截面求和，得到材料的多群虚宏观总截面。

蒙特卡罗粒子输运计算中，从源中抽样源粒子，能量为14MeV。

粒子输运长度计算：

首先查找14MeV在多群统一能量网格数组中的位置，得到位置p＝3。当前能量下材料的多群虚宏观总截面值Σ_MG为多群虚宏观总截面数组中的第3个与第4个值中的最大值，U-233、U-234、U-235三个核素各自的虚总截面σ_MG,U-233，σ_MG,U-234，σ_MG,U-235取各自的多群虚总截面数组的第3个与第4个值中的最大值。抽样一个随机数ξ，粒子输运长度L＝-ln(ξ/Σ_MG)。

反应核素抽样：

根据U-233、U-234、U-235三个核素各自的虚总截面σ_MG,U-233，σ_MG,U-234，σ_MG,U-235与材料的多群虚宏观总截面值Σ_MG的比值进行离散抽样，得到反应核素为U-233。

拒绝抽样：

根据U-233的多群核素指针数组中的第3个值与第4个值，得到14MeV能点在U-233的能量网格中的区域，在此区域进行二分查找得到14MeV在核素U-233的能量网格中的位置p＝24833，用线性插值法得到插值因子frac。用核素U-233数据库中的第p和p+1个总截面，线性插值得到14MeV时U-233的总截面σ_U-233。使用σ_U-233与σ_MG,U-233进行拒绝抽样，抽取随机数ξ₁，如果ξ₁≤σ_U-233/σ_MG,U-233则反应被接收，进行反应类型抽样；如果ξ₁＞σ_U-233/σ_MG,U-233，则反应被拒绝，重新进行粒子输运长度抽样。

反应类型抽样：

取出核素U-233数据库中的第p和p+1个吸收截面、弹散截面，线性插值得到14MeV时U-233的吸收截面与弹散截面值。非弹截面＝总截面-吸收截面-弹散截面值。使用离散抽样方法得到反应类型。

MC计算中，得到反应类型后调用具体反应类型处理函数进行计算，最后对计算结果进行统计得到面积分流量。由于本发明计算实例中的这部分计算与统计与传统处理方法一致，因此不再详细阐述

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。