CN104346533A - 一种蒙特卡罗粒子输运模拟中核截面数据处理优化方法 - Google Patents

Abstract

本发明公布了一种蒙特卡罗粒子输运模拟中核截面数据处理优化方法，在传统MC粒子输运模拟连续点截面处理与应用的方法基础上，归并问题涉及的所用核素的能量网格点形成统一能量网格，使得每步输运只需对统一能量网格进行一次搜索查找，由预处理得到的核素指针数组直接找到当前能量在每个核素能量网格中的具体位置，替代传统方法对材料中每个核素的不同能量网格的重复搜索查找，极大地减少了蒙特卡罗粒子输运计算中对能量网格的搜索查找次数，在不损失计算精度的情况下提高计算速度；并根据点状核数据能量点分布的特点，引入能量分段思想形成能量网格的双层查找模式，减少在分布密度较高区域的无效搜索次数，提高单次网格搜索的效率。

Description

一种蒙特卡罗粒子输运模拟中核截面数据处理优化方法

技术领域

本发明涉及蒙特卡罗粒子输运计算中的一种核截面数据处理优化方法，属于原子能科学技术领域。

背景技术

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

统一能量网格是Serpent等软件使用的一种加速MC计算的方法。其基本思想是：构造所用核素适用的能量网格，计算时，只需对该能量网格进行一次搜索，即可得到当前能量下所有核素的截面数据，避免了因各核素能量网格不同而进行的重复的能量网格搜索。但是，已有的统一能量网格方法，考虑到计算速度、计算精度与内存开销，通常有两种处理方法：一种为了保证原始核数据的精度，通常将所有核素的能量网格合并，并对所有能点预计算所有核素的各种截面和每种材料的各种宏观截面，存储在内存中供输运计算使用，因此会占用非常大的内存，甚至因为需要内存太大而不能进行计算；另外一种为了避免大内存开销，在建立统一能量网格时，不是合并所有核素的能量网格，而是自定义一套能量网格或者对所有核素的能量网格进行精简，这样的处理导致了原始核数据精度的部分丢失。如何在保持计算精度且不显著增加内存使用的情况下，提高统一能量网格方法对MC计算的加速效果是需要解决的关键问题。

发明内容

本发明要解决的技术问题为：本发明的目的在于改善传统蒙特卡罗粒子输运计算方法在核截面数据处理与使用效率上的不足，提供了一种蒙特卡罗核截面数据处理优化方法，在保持计算精度的前提下，提高计算速度。

本发明的技术方案如下：一种蒙特卡罗程序中核截面数据处理优化方法，其特征在于包括以下步骤：

(1)对蒙特卡罗程序从数据库中读取的核数据进行预处理，包括以下内容：

a.统一能量网格的建立：在读取核素数据的时候，采用归并法合并问题涉及的所用核素的能量网格形成统一能量网格数组UnionErgArray[N_tot]，其中N_tot表示统一能量网格中的网格数；

b.核素指针数组的建立：对比每个单核素i的能量网格EryArray_i[N_i]，N_i表示网格数，与统一能量网格数组UnionErgArray[N_tot]，找出统一能量网格数组中的每个能点UnionErgArray[j]在单核素i能量网格中的位置，标记为Position_i[j]，存入数组，形成核素i的核素指针数组，标记为Position_i[N_tot]；

c.核素总截面数组与材料宏观总截面数组的建立：在核素指针数组建立的同时，依据单核素i的能量网格值计算线性插值因子：

$\mathrm{factor}=\frac{\mathrm{UnionEryArray}\left[j\right]-\mathrm{EryArra}{y}_i\left[{\mathrm{Position}}_i\right[j\left]\right]}{{\mathrm{EryArray}}_i\left[{\mathrm{Position}}_i\right[j]+1]-{\mathrm{EryArray}}_i\left[{\mathrm{Position}}_i\right[j\left]\right]},$

插值得到核素i的总截面：

σ＝σ[Position_i[j]]+factor*(σ[Position_i[j]+1]-σ[Position_i[j]])，

其中σ[Position_i[j]]表示单核素能量网格对应的总截面数组中的第Position_i[j]个值，将总截面值存入数组，形成单核素i的总截面数组，标记为σ_i[N_tot]；将材料k中所有核素的总截面值按密度比例求和得到材料的宏观总截面数组，标记为Σ_k[N_tot]；

d.分段能点指针数组的建立：根据用户给定的或者程序内置的分段能点数组PointwiseErgArray[N_pw]，其中N_pw表示分段能点个数，遍历统一能量网格，找出分段能点数组中每个能点在统一能量网格中的位置，存入数组，形成分段能点指针数组，标记为Position_pw[N_pw]；

(2)在蒙特卡罗粒子输运模拟中需要使用截面数据时，对各类截面进行计算，包括以下步骤：

a.查找当前能量在统一能量网格数组中的位置，首先，在分段能点数组PointwiseErgArray[N_pw]中搜索当前能量erg，得到erg在分段能点数组中的位置为N_pw1与N_pw1+1之间；其次，根据分段能点指针数组，得到erg在统一能量网格数组中的位置为p1＝Position_pw[N_pw1]与p2＝Position_pw[N_pw1+1]之间；再次，使用二分查找法在统一能量网格数组UnionErgArray[N_tot]的p1到p2区间上进行搜索得到erg在统一能量网格中的精确位置Mainindex，并得到插值因子：

$\mathrm{fa}\min =\frac{\mathrm{erg}-\mathrm{UnionEryArray}\left[\mathrm{Mainindex}\right]}{\mathrm{UnionEryArray}[\mathrm{Mainindex}+1]-\mathrm{UnionEryArray}\left[\mathrm{Mainindex}\right]};$

b.在抽样粒子输运长度时需要计算材料的宏观截面，直接依据当前材料的宏观总截面数组Σ_k[N_tot]进行线性插值得到材料的宏观截面Σ：

Σ＝Σ_k[Mainindex]+fmain×(Σ_k[Mainindex]-Σ_k[Mainindex])；

c.在抽样反应核素时需要使用当前材料中各核素的总截面，根据步骤(2)a中得到的当前能量erg在统一能量网格中的精确位置Mainindex和插值因子fmain，同步骤(2)b中，直接依据每个核素的总截面数组σ_i[N_tot]进行线性插值得到每个核素的总截面σ_i：

σ_i＝σ_i[Mainindex]+fmain×(σ_i[Mainindex]-σ_i[Mainindex])；

d.在抽样反应类型或物理参数计算时需要计算当前反应核素i的多种反应截面，根据当前能量erg在统一能量网格中的精确位置Mainindex，对应当前核素i的核素指针数组Position_i[Mainindex]得到erg在核素i的能量网格中的位置p，线性插值因子f由下式得到：

f＝(erg-EryArray_i[p])/(EryArray_i[p+1]-EryArray_i[p])，

线性插值得到各种反应截面。

所述的蒙特卡罗核截面数据处理优化方法，其特征在于使用归并法建立问题涉及的所用核素的统一能量网格，每步输运只对统一能量网格进行一次搜索查找，由步骤(1)b中得到的核素指针数组直接找到当前能量在每个核素能量网格中的具体位置，替代了对材料中每个核素的不同能量网格的搜索查找，极大地减少了蒙特卡罗粒子输运计算中能量网格的搜索查找次数，在不损失计算精度的情况下提高计算速度；在计算材料宏观截面的时候，直接使用预计算得到的宏观截面数组进行插值，而避免了对材料中每个核素截面插值再求和的计算过程；并根据点状核数据能量点分布的特点，引入能量分段思想形成能量网格的双层查找模式，减少在分布密度较高区域的无效搜索次数，提高单次网格搜索的效率。

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

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

本发明与现有技术相比的优点在于：本发明的蒙特卡罗核截面数据处理的优化方法，可以在不降低计算精度且不显著增大内存开销的情况下，得到更好的计算速度。

附图说明

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

具体实施方式

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

1.对蒙特卡罗程序从数据库中读取的核数据进行预处理，通过如下步骤实现：

第一步，建立统一能量网格：在读取核素数据的时候，采用归并法合并问题涉及的所用核素的能量网格形成统一能量网格数组UnionErgArray[N_tot]；

第二步，核素指针数组的建立：对比单核素的能量网格EryArray_i[N_i]与统一能量网格数组UnionErgArray[N_tot]，对统一能量网格中的每个能点UnionErgArray[j]标注其在单核素能量网格中的位置Position_i[j]，存入数组，形成每个核素的核素指针数组Position_i[N_tot]；

第三步，核素与材料的总截面数组建立：在核素指针数组建立的同时，依据单核素能量网格值计算线性插值因子：

$\mathrm{factor}=\frac{\mathrm{UnionEryArray}\left[j\right]-\mathrm{EryArra}{y}_i\left[{\mathrm{Position}}_i\right[j\left]\right]}{{\mathrm{EryArray}}_i\left[{\mathrm{Position}}_i\right[j]+1]-{\mathrm{EryArray}}_i\left[{\mathrm{Position}}_i\right[j\left]\right]},$

插值得到核素的总截面：

σ＝σ[Position_i[j]]+factor*(σ[Position_i[j]+1]-σ[Position_i[j]])，

存入数组，形成每个核素的总截面数组σ_i[N_tot]；对每个材料k中涉及的核素按能量密度比例求和得到每个材料的宏观总截面数组Σ_k[N_tot]；

第四步，分段能点指针数组的建立：根据用户给定的或者程序内置的分段能点数组PointwiseErgArray[N_pw]，遍历统一能量网格，找出其在统一能量网格中的位置，存入数组，形成分段能点指针数组Position_pw[N_pw]；

2.在蒙特卡罗粒子输运程序中需要使用截面数据时，对各类截面进行计算，通过如下步骤实现：

第一步，查找当前能量在统一能量网格数组中的位置，首先，在分段能点数组PointwiseErgArray[N_pw]中搜索当前能量erg，得到erg在分段能点数组中的位置为N_pw1与N_pw1+1之间；其次，根据分段能点指针数组，得到erg在统一能量网格数组中的位置为p1＝Position_pw[N_pw1]与p2＝Position_pw[N_pw1+1]之间；再次，使用二分查找法在统一能量网格数组UnionErgArray[N_tot]的p1到p2区间上进行搜索得到erg在统一能量网格中的精确位置Mainindex，并得到插值因子：

$\mathrm{fa}\min =\frac{\mathrm{erg}-\mathrm{UnionEryArray}\left[\mathrm{Mainindex}\right]}{\mathrm{UnionEryArray}[\mathrm{Mainindex}+1]-\mathrm{UnionEryArray}\left[\mathrm{Mainindex}\right]};$

第二步，在抽样粒子输运长度时需要计算材料的宏观截面，直接依据当前材料的宏观总截面数组Σ_k[N_tot]进行线性插值得到材料的宏观截面Σ：

Σ＝Σ_k[Mainindex]+fmain×(Σ_k[Mainindex]-Σ_k[Mainindex])；

第三步，在抽样反应核素时需要使用当前材料中各核素的总截面，根据上一步中得到的当前能量erg在统一能量网格中的精确位置Mainindex和插值因子fmain，直接依据每个核素的总截面数组σ_i[N_tot]进行插值可以得到每个核素的总截面σ_i；

第四步，在抽样反应类型时需要计算当前反应核素i的多种反应截面，根据当前能量erg在统一能量网格中的精确位置Mainindex，对应当前核素i的核素指针数组Position_i[Mainindex]得到erg在核素i的能量网格中的位置p，线性插值因子：

f＝(erg-EryArray_i[p])/(EryArray_i[p+1]-EryArray_i[p])，

线性插值得到各种反应截面；

第四步，在临界计算或者计算反应率等物理量时，需要使用当前材料中每个核素的吸收、散射等反应截面，计算方法同上一步。

3.输出结果在MC粒子输运程序中直接使用，包含以下内容：

材料在当前能量时的宏观截面、材料中的各核素在当前能量时总截面、吸收截面、弹散截面、裂变截面、裂变中子数、热参数。

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

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

从数据库中读取据核素U-233、U-234、U-235的能量网格，分别为25036、27637、77165个数据点，使用归并法合并三个核素的能量网格得到统一能量网格，含118991个网格点。

对比U-233的能量网格与统一能量网格，得到统一能量网格中的每一个点在U-233的能量网格中的位置，即核素指针数组，并用插值法计算统一能量网格中的每一个点上U-233的总截面值；同样对U-234、U-235建立核素指针数组并计算总截面。

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

对比分段能量网格与统一能量网格，得到分段能量网格中的每个一点在统一能量网格中的位置，即分段能点指针数组。

MC计算中，从源中抽样源粒子，能量为14MeV。

粒子输运长度计算：

首先查找14MeV在分段网格中的位置，得到位置p＝3。

读取分段能点指针数组的第p＝3个值(110326)，与第p+1＝4个值(118991)，得到14MeV在统一能量网格中的位置为第110326个点到118991个点之间。

在统一能量网格中的第110326个点到118991个点之间使用二分查找法得到14MeV在第118556(E＝14MeV)与118557(E＝14.5MeV)点之间，即index＝118556，使用线性插值得到插值因子为0。

读取材料的宏观截面数组中的第index与index+1个值，使用线性插值得到14MeV时的材料宏观截面，抽样得到粒子输运长度L＝-ln/。

反应核素抽样：

读取核素U-233、U-234、U-235的总截面数组中的第index与index+1个值，使用线性插值得到14MeV时的核素反应总截面，使用离散抽样法得到反应核素(U-233)。

反应类型抽样：

根据U-233的核素指针数组中的第index值得到14MeV在核素U-233的能量网格中的位置p＝24833，用线性插值法得到插值因子frac。取出核素U-233数据库中的第p和p+1个吸收截面、弹散截面，线性插值得到14MeV时U-233的吸收截面与弹散截面值。非弹截面＝总截面-吸收截面-弹散截面值。使用离散抽样方法得到反应类型。

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

Claims (1)

1.一种蒙特卡罗粒子输运模拟中核截面数据处理优化方法，其特征在于包括以下步骤：

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

(11)统一能量网格的建立：在读取核素数据时，采用归并法合并问题涉及的所用核素的能量网格，形成统一能量网格数组，标记为UnionErgArray[N_tot]，其中N_tot表示统一能量网格中的网格数；

(12)核素指针数组的建立：对比每个单核素i的能量网格EryArray_i[N_i]，N_i表示网格数，与步骤(11)中产生的统一能量网格数组UnionErgArray[N_tot]，找出统一能量网格数组中的每个能点UnionErgArray[j]在单核素i能量网格EryArray_i[N_i]中的位置，标记为Position_i[j]，存入数组，形成核素i的核素指针数组，标记为Position_i[N_tot]；

(13)核素总截面数组与材料宏观总截面数组的建立：在步骤(12)中建立核素指针数组的同时，依据单核素i的能量网格值计算线性插值因子：

$\mathrm{factor}=\frac{\mathrm{UnionEryArray}\left[j\right]-{\mathrm{EryArray}}_i\left[{\mathrm{Position}}_i\right[j\left]\right]}{{\mathrm{EryArray}}_i\left[{\mathrm{Posotoon}}_i\right[j]+1]-{\mathrm{EryArray}}_i\left[{\mathrm{Position}}_i\right[j\left]\right]},$

插值得到核素i的总截面：

σ＝σ[Position_i[j]]+factor*(σ[Position_i[j]+1]-σ[Position_i[j]])，

其中σ[Position_i[j]]表示单核素i能量网格对应的总截面数组中的第Position_i[j]个值，将总截面值存入数组，形成单核素i的总截面数组，标记为σ_i[N_tot]，将材料k中所有核素的总截面值按密度比例求和得到材料的宏观总截面数组，标记为Σ_k[N_tot]；

(14)分段能点指针数组的建立：根据用户给定的或者程序内置的分段能点数组PointwiseErgArray[N_pw]，其中N_pw表示分段能点个数，遍历步骤(11)中产生的统一能量网格数组UnionErgArray[N_tot]，找出分段能点数组中每个能点在统一能量网格中的位置，存入数组，形成分段能点指针数组，标记为Position_pw[N_pw]；

步骤(2)在蒙特卡罗粒子输运模拟中需要使用截面数据时，对各类截面进行计算，包括以下步骤：

(21)查找当前能量在步骤(11)中产生的统一能量网格数组UnionErgArray[N_tot]中的位置，首先，在分段能点数组PointwiseErgArray[N_pw]中搜索当前能量erg，得到erg在分段能点数组中的位置为N_pw1与N_pw1+1之间；其次，根据步骤(14)中得到的分段能点指针数组PointwiseErgArray[N_pw]，得到erg在统一能量网格数组中的位置为p1＝Position_pw[N_pw1]与p2＝Position_pw[N_pw1+1]之间；再次，使用二分查找法在统一能量网格数组UnionErgArray[N_tot]的p1到p2区间上进行搜索得到erg在统一能量网格中的精确位置Mainindex，并得到插值因子：

$\mathrm{fmain}=\frac{\mathrm{erg}-\mathrm{UnionEryArray}\left[\mathrm{Mainindex}\right]}{\mathrm{UnionEryArray}[\mathrm{Mainindex}+1]-\mathrm{UnionEryArray}\left[\mathrm{Mainindex}\right]};$

(22)在抽样粒子输运长度时需要计算材料的宏观截面，依据步骤(21)中得到的当前能量erg在统一能量网格中的精确位置Mainindex和插值因子fmain，对步骤(13)中得到的当前材料的宏观总截面数组Σ_k[N_tot]进行线性插值得到材料的宏观截面Σ：

Σ＝Σ_k[Mainindex]+fmain×(Σ_k[Mainindex]-Σ_k[Mainindex])；

(23)在抽样反应核素时需要使用当前材料中各核素的总截面，根据步骤(21)中得到的当前能量erg在统一能量网格中的精确位置Mainindex和插值因子fmain，同步骤(22)中，直接依据步骤(13)中得到的每个核素的总截面数组σ_i[N_tot]进行线性插值得到每个核素的总截面σ_i：

σ_i＝σ_i[Mainindex]+fmain×(σ_i[Mainindex]-σ_i[Mainindex])；

(24)在抽样反应类型或物理参数计算时需要计算当前反应核素i的多种反应截面，根据步骤(21)中得到的当前能量erg在统一能量网格中的精确位置Mainindex，对应步骤(12)中得到的当前核素i的核素指针数组Position_i[Mainindex]得到erg在核素i的能量网格中的位置p，线性插值因子f由下式得到：

f＝(erg-EryArray_i[p])/(EryArray_i[p+1]-EryArray_i[p])，

线性插值得到各种反应截面。