一种基于蒙特卡罗方法与特征线方法耦合的计算辐射屏蔽的方法

技术领域

本发明涉及一种基于蒙特卡罗方法与特征线方法耦合的计算辐射屏蔽的方法，属于核科 学中辐射输运与数值计算模拟研究领域。

背景技术

辐射屏蔽计算主要是利用粒子输运计算方法来对大型核装置进行粒子输运模拟计算，从 而获得核装置生物屏蔽区域的辐射粒子通量分布的过程，它可以辅助核工程设计人员对核装 置进行有效地辐射屏蔽设计与分析。目前，粒子输运计算方法一般分为蒙特卡罗方法和确定 论方法。蒙特卡罗方法是以统计学理论为基础，通过对大量粒子物理事件和过程的随机模拟 来获取粒子通量密度。确定论方法常见的有特征线方法、离散纵标法、碰撞概率方法等，此 类方法主要是从粒子输运方程本身出发进行数值求解来获得粒子通量密度。

蒙特卡罗方法的优点在于能够比较逼真地描述具有随机性质的事物的特点及过程，并且 对几何和材料的限制小，可以精确模拟复杂材料、复杂几何下的粒子输运问题；它的缺点在 于计算比较耗时，收敛速度慢，且在粒子输运问题中，计算结果与系统大小有关，对于大系 统或小概率事件的计算问题，计算结果往往比真实值偏低，特别针对大型核装置进行屏蔽计 算时，容易产生深穿透效应而无法获得精确的计算结果。确定论方法的最大优点在于计算收 敛速度快；但它的缺点在于存在一定程度的几何限制性，难以处理大型复杂几何模型。

传统的辐射屏蔽计算方法通常基于单一的蒙特卡罗方法或者确定论方法，从而两种方法 本身的缺点严重限制了它们在大型复杂核装置（如ITER模型、IFMIF装置）辐射屏蔽计算 中的应用。蒙特卡罗方法屏蔽计算时对大尺度的屏蔽区域难以获得精确的结果，统计误差大， 且收敛速度慢，计算耗时；而确定论方法存在一定的几何限制性，难以处理大型复杂核装置 的辐射屏蔽计算问题。

发明内容

本发明要解决的技术问题为：本发明的目的在于克服传统的辐射屏蔽计算方法在计算精 度、计算效率和几何适应性上的不足，提供了一种蒙特卡罗-特征线方法耦合的辐射屏蔽计 算方法，不仅可以获得更好的计算精度和计算效率，而且增强了对计算模型的几何适应性。

本发明解决上述技术问题采用的技术方案为：一种基于蒙特卡罗方法与特征线方法耦合 的计算辐射屏蔽的方法，包括以下步骤：

步骤（1）、几何模型划分：

分析待计算的几何模型，根据几何模型各个部分几何结构的复杂度，将几何复杂的辐射 源及其周围支撑结构部分划为蒙卡子模型，将几何简单的远离辐射源的外围屏蔽层区域划为 特征子模型，所述的远离辐射源为≥10个自由程，两个子模型之间通过一系列公共耦合面连 接起来组成整个完整的几何模型；

步骤（2）、获得计算参数，包括以下内容：

对蒙卡子模型进行粒子输运计算，获得蒙卡子模型与特征子模型之间公共耦合面上的面 积分流量矩阵：

J_I,G＝[j_i，g]，(i＝1~I,g＝1~G)

其中，j_i,g是穿过耦合面上第i个面元的第g群的面积分流量，I为耦合面上面元的总数， G为材料宏观截面划分的能群总数；

蒙卡子模型的物理计算模型信息，由用户手工给出或者利用自动建模技术处理获得；

特征子模型中各几何体的重要性信息，由用户给出；

特征子模型的网格栅元划分信息，由用户给出；

耦合面面元参数信息，可根据特征线子模型的网格栅元划分信息获得；

粒子运动方向离散信息，由用户给出；

特征子模型的特征线信息，利用通用软件ACIS软件处理获得；

特征子模型中各材料区的宏观截面信息，可利用截面处理软件对评价核数据库ENDF/B 进行并群处理得到；

步骤（3）、利用如下非均匀修正的特征线方法进行计算：

${\phi }_{j,g}=\frac{1}{V_j}{\int }_{4\pi }\underset{k=1}{\overset{K}{\Sigma }}{A}_{j,k}\left(\Omega \right){\int }_0^{L_{j,k}\left(\Omega \right)}[{\psi }_{j,k,g}^{\mathrm{in}}\left(\Omega \right)e^{{-\Sigma }_{j,k,g}s\left(\Omega \right)}+\frac{Q_{j,g}\left(\Omega \right)}{{\Sigma }_{j,k,g}}(1-e^{{-\Sigma }_{j,k,g}s\left(\Omega \right)})]{\mathrm{dsd}}^2\Omega$

其中，φ_j,g为网格栅元j里的第g群的粒子标通量密度；V_j为网格栅元j的体积；A_j,k(Ω) 和L_j，k(Ω)分别为网格栅元j里、离散方向为Ω的第k条特征线的宽度和长度；s(Ω)为网格 栅元j里、离散方向为Ω的第k条特征线上某一点距离该特征线起点的距离；∑_j，k,g为网格栅 元j里材料的第g群宏观总截面；Q_j,g(Ω)为网格栅元j里、离散方向为Ω的第g群的辐射源 项；为网格栅元边界上的特征线入射粒子角通量密度，此处即为网格栅元j里、沿 Ω方向上的第k条特征线的第g群的特征线入射粒子角通量密度；

进行特征线方法计算时，根据各几何体的重要性进行网格剖分，对重要性高的几何体， 采用精细网格剖分；对重要性低的几何体，采用粗网格剖分。

其中，所述的几何模型划分，将待计算的几何模型划分为蒙卡子模型和特征子模型两部 分，并采用自动建模技术进行自动建模与模型划分，无需人工干预，可快速、方便、高效地 生成蒙卡子模型和特征子模型的物理计算模型，这样与传统的手工建模相比，自动建模可以 处理更复杂的几何模型，提高了生成物理模型的正确性和处理效率。

其中，所述的特征子模型网格栅元划分，根据组成特征子模型各几何体的重要性进行网 格划分，对重要性高的几何体，采用精细网格剖分，对重要性低的几何体采用粗网格剖分， 这种处理方法可以在保证用户关心区域计算精度的同时，提高计算效率；同时，减少网格数 目亦可以减少内存开销，提高内存使用效率，增强其对大型复杂几何模型的几何适应性。

其中，所述的非均匀修正的特征线方法，在模拟计算过程中，对材料非均匀的网格，需 要先进行材料宏观总截面修正，再利用修正后的宏观总截面进行特征线方法计算，这样与传 统的特征线方法相比，提高了计算精度，同时，也提高了计算效率。

本发明的原理为：

几何模型的划分可利用自动建模技术来处理，获得蒙特卡罗计算需要的物理计算模型信 息文件和特征线法计算时需要的特征线信息文件；

耦合面上的面积分流量矩阵J_I,G可由国际上通用的蒙特卡罗粒子输运计算程序（如 GEANT4）计算获得；

蒙特卡罗方法计算时材料采用的能量点截面数据库采用中国原子能院发布的 CENDL2.0；

特征线方法计算时采用的多群宏观截面数据采用国际标准例题的参考数据或者利用评 价核数据库并群处理获得；

对屏蔽区域进行粒子输运模拟时，基于特征线理论，采用源迭代法进行迭代计算；

本发明公布的辐射屏蔽计算方法，不仅可以得到更高的计算精度，而且能够获得更好的 计算效率，并增强了几何适应性。

附图说明

图1为本发明几何模型划分示意图；

图2为本发明蒙卡子模型与特征子模型之间数据耦合处理示意图；

图3为本发明特征线示意图；

图4为本发明计算流程示意图；

图5为本发明计算实例几何模型示意图。

图中，1为蒙卡子模型，2为中子源，3为公共耦合面，4为特征子模型（屏蔽层区域）， 面元法向向量n，中子运动方向Ω，面元ΔS，5为栅元边界。

具体实施方式

下面结合附图和具体实施例进一步说明本发明。

一种基于蒙特卡罗方法与特征线方法耦合的计算辐射屏蔽的方法，

1、几何模型划分与物理计算模型信息获取：

第一步，定义模型蒙卡子模型和特征子模型的范围。对大型核装置而言（如图1），通常 其芯部结构比较复杂，主要体现出芯部区域几何结构组成非常复杂、中子源非常规源分布、 芯部材料非均匀性强等特点，因此具有这些特点的模型适合用蒙卡方法来处理，故将芯部区 域划为蒙卡子模型。在远离装置芯部、靠近外界环境的屏蔽区域通常几何结构和材料组成相 对比较简单，适合用特征线方法来处理，故将其划为特征子模型。

第二步，如图1所示，将整个模型用公共耦合面分割成蒙卡子模型和特征子模型两部分。

利用蒙卡子模型的物理计算模型信息，对于简单模型可以手工撰写，对于复杂模型，可 利用自动建模技术来获得；特征子模型的特征线信息，可利用AutoCAD、ACIS处理获得， 同时，在进行特征子模型网格栅元划分时，针对重要性高的几何采用精细网格划分，对重要 性低的几何采用粗网格划分，最终处理获得特征子模型的物理计算模型信息。

2、核数据库准备：

第一步，为蒙特卡罗计算准备核数据库。可以直接利用中国原子能院发布的CENDL2.0 数据库。

第二步，为特征线法计算准备宏观截面数据文件。根据美国布鲁克海文国家实验室发布 的评价核数据库ENDB/B，利用截面处理软件对其进行并群处理来获得各材料区的宏观截面 数据；若是国际标准例题，也可以根据相关机构发布的标准例题文档直接查阅宏观截面信息。

3、对蒙卡子模型进行粒子输运计算：

利用第1步获得的蒙特卡罗计算输入文件和第2步准备的CENDL2.0数据库，对几何模 型的蒙卡子模型部分进行粒子输运模拟计算，根据如下的物理公式：

$J(r,\mathrm{\Delta S},\mathrm{\Delta E},\mathrm{\Delta \Omega })=\underset{\mathrm{\Delta S}}{\int }\underset{\mathrm{\Delta E}}{\int }\underset{\mathrm{\Delta \Omega }}{\int }|\Omega ·n|\phi (r,E,\Omega )\mathrm{d\Omega dEdr}$

$=\underset{n=1}{\overset{N}{\Sigma }}{w}_n\underset{\mathrm{\Delta S}}{\int }\mathrm{dr}\underset{\mathrm{\Delta E}}{\int }\mathrm{dE}\underset{\mathrm{\Delta \Omega }}{\int }\mathrm{d\Omega \delta }(\Omega ,{\Omega }_n)\delta (E,E_n)\delta (r,r_n)$

其中：

φ(r,E,Ω)为r处的粒子角通量密度，即单位时间内能量为E的穿过与Ω方向相垂直的 单位面积上的粒子数；

J(r,ΔS,ΔE,ΔΩ)为在r处沿方向Ω穿过面元ΔS且能量在E处ΔE间隔内的粒子数 目；

w_n为穿过面元ΔS的粒子的权重；

r_n、E_n、Ω_n分别为粒子的状态信息，即空间位置、能量和运动方向；

由上面的物理模型可知，只要粒子的运动状态落在相空间（ΔS、ΔE、ΔΩ）内，则该粒 子就会对J(r,ΔS,ΔE,ΔΩ)有贡献，从而最终获得耦合面上各网格面元的粒子面积分流量 J_I,G。

4、蒙卡子模型与特征子模型之间的数据耦合处理：

蒙卡子模型与特征子模型之间是通过公共面上的粒子面积分流量矩阵耦合起来的。特征 线方法进行计算时，耦合面对特征子模型而言是一个有源表面，因此需要知道该源表面上的 边界条件，即表面上的的入射粒子通量密度因此需要利用蒙特卡罗计算获得的 面积分流量矩阵J_I,G进行转换来获得。如图2所示，在耦合面上针对某一面元ΔS，其指向 特征子模型的法线向量为n，将该面元源法向拉伸δ形成体元，该体元的体积可近似为： ΔV＝δΔS

当粒子沿Ω方向穿过面元时，粒子在体元ΔV内轨迹长度为：

$T_l=\frac{\delta }{\cos \theta }=\frac{\delta }{n·\Omega }$

则根据粒子通量密度的定义，则体元ΔV内的平均粒子通量密度为：

$\phi (r,E,\Omega )=\underset{\delta \to 0}{\lim }\left(\frac{{\mathrm{WT}}_l}{\mathrm{\Delta V}}\right)=\underset{\delta \to 0}{\lim }\left(\frac{\mathrm{W\delta }}{\mathrm{\delta \Delta S}(\delta ·\Omega )}\right)=\frac{W}{\mathrm{\Delta S}(\delta ·\Omega )}$

其中W为Ω方向穿过面元ΔS的能量为E的粒子数目，该参数即为蒙特卡罗计算时得 到的面积分流量；φ(r,E,Ω)则为特征线法计算时需要的边界条件。

5、对特征子模型进行特征线法计算：

如图3所示，第1步获得的特征线信息文件内容主要包括：

（a）粒子运动方向的离散求积组信息{Ω_m,ω_m}；

（b）各离散方向下、各几何栅元所有特征线的宽度A_j，k(Ω)；

（c）各离散方向下、各几何栅元所有特征线段长度L_j，k(Ω)；

（d）各几何栅元的体积V_j；

利用特征线信息文件和第4步计算获得的耦合面上的源边界条件，根据如下计算模型计 算各几何栅元各能群的粒子标通量密度φ_j,g：

${\phi }_{j,g}=\frac{1}{V_j}{\int }_{4\pi }\underset{k=1}{\overset{K}{\Sigma }}{A}_{j,k}\left(\Omega \right){\int }_0^{L_{j,k}\left(\Omega \right)}[{\psi }_{j,k,g}^{\mathrm{in}}\left(\Omega \right)e^{{-\Sigma }_{j,k,g}s\left(\Omega \right)}+\frac{Q_{j,g}\left(\Omega \right)}{{\Sigma }_{j,k,g}}(1-e^{{-\Sigma }_{j,k,g}s\left(\Omega \right)})]{\mathrm{dsd}}^2\Omega$

$=\frac{1}{V_j}\underset{m=1}{\overset{M}{\Sigma }}{\omega }_m\underset{k=1}{\overset{K}{\Sigma }}{A}_{j,k}\left({\Omega }_m\right)L_{j,k}\left({\Omega }_m\right)(\frac{Q_{j,g}\left({\Omega }_m\right)}{{\Sigma }_{j,k,g}}+\frac{{\psi }_{j,k,g}^{\mathrm{in}}\left(\Omega \right)-{\psi }_{j,k,g}^{\mathrm{out}}\left({\Omega }_m\right)}{{\Sigma }_{j,k,g}{L}_{j,k}\left({\Omega }_m\right)})$

其中为特征线段k在栅元边界处沿Ω方向的出射粒子角通量密度，可由如下 公式计算得到：

${\psi }_{j,k,g}^{\mathrm{out}}\left({\Omega }_m\right)={\psi }_{j,k,g}^{\mathrm{in}}\left({\Omega }_m\right)e^{{-\Sigma }_{j,k,g}{L}_{j,k}\left({\Omega }_m\right)}+\frac{Q_{j,g}\left({\Omega }_m\right)}{{\Sigma }_{j,k,g}}(1-e^{{-\Sigma }_{j,k,g}{L}_{j,k}\left({\Omega }_m\right)})$

其中，对于非均匀网格栅元进行宏观总截面修正的条件为：当时进行修正：

为当前网格栅元的估算体积；V_真实为当前网格 栅元的真实体积；α为非均匀修正系数，通常大于1.0；修正后的宏观总截面为：

即用修正后的宏观总截面代替原始的宏观总截面进行计算。

6、整个计算的流程如图4所示，最终计算输出结果包含以下内容：

特征子模型的粒子通量密度分布文件。

7、计算实例，包含以下内容：

计算实例的几何模型（图5）是一个20×20×140cm的长方体模型，在底部中心有一个 高为6cm，半径为2cm的圆柱中子源，中子源能量为14MeV，方向服从各向同性分布，模 型填充的屏蔽材料为60%Fe和40%H₂O的混合物，需要计算屏蔽层区域的中子通量分布。

如图5所示，在距离底部10cm的位置作一个耦合面，将整个模型分成上、下的特征子 模型和蒙卡子模型两部分，并将特征子模型划分成2×2×2cm的网格模型。同时，蒙卡子模 型的物理计算模型信息文件采用自动建模技术获得，特征线信息文件利用ACIS处理获得。

材料的数据库采用评价数据库CENDL2.0和利用截面处理软件并群得到的宏观截面数 据。

首先根据蒙卡子模型的物理计算模型信息文件和CENDL2.0数据库，利用蒙特卡罗粒子 输运计算程序GEANT4对蒙卡子模型进行中子输运模拟，获得耦合面上的面积分流量矩阵 J_I,G：

$J_{I,G}\left|\left(\begin{array}{cccc}{j}_{\mathrm{1,1}}& j_{\mathrm{1,2}}& ······& j_{1,G}\\ j_{\mathrm{2,1}}& j_{\mathrm{2,2}}& ······& j_{2,G}\\ ····& ····& ······& ····\\ j_{I,1}& j_{I,2}& ······& j_{I,G}\end{array}\right)\right|(i=1~I,g=1~G)$

其中j_i,g为耦合面上面元i的、第g群的、各离散方向下的中子面积分流量，可表示为： j_i,g＝|j_i,g(Ω₁)j_i,g(Ω₂)j_i，g(Ω₃)......j_i,g(Ω_M)|

特征子模型耦合面上的源边界条件由J_I,G转换得到（以几何栅元j为例，且该栅元在耦 合面上，假设栅元j与面元i对应）：

${\psi }_{j,k,g}^{\mathrm{in}}\left(\Omega \right)=\frac{j_{i,g}\left(\Omega \right)}{{\mathrm{\Delta S}}_i(\delta ·\Omega )}$

得到特征子模型在耦合面上的源边界条件以后，即可根据特征线模型：

${\phi }_{j,g}=\frac{1}{V_j}{\int }_{4\pi }\underset{k=1}{\overset{K}{\Sigma }}{A}_{j,k}\left(\Omega \right){\int }_0^{L_{j,k}\left(\Omega \right)}[{\psi }_{j,k,g}^{\mathrm{in}}\left(\Omega \right)e^{{-\Sigma }_{j,k,g}s\left(\Omega \right)}+\frac{Q_{j,g}\left(\Omega \right)}{{\Sigma }_{j,k,g}}(1-e^{{-\Sigma }_{j,k,g}s\left(\Omega \right)})]{\mathrm{dsd}}^2\Omega$

利用特征线信息文件和源边界条件，得到几何栅元j的第g群的中子标通量密度为：

${\phi }_{j,g}=\frac{1}{V_j}\underset{m=1}{\overset{M}{\Sigma }}{\omega }_m\underset{k=1}{\overset{K}{\Sigma }}{A}_{j,k}\left({\Omega }_m\right)L_{j,k}\left({\Omega }_m\right)(\frac{Q_{j,g}\left({\Omega }_m\right)}{{\Sigma }_{j,k,g}}+\frac{{\psi }_{j,k,g}^{\mathrm{in}}\left(\Omega \right)-{\psi }_{j,k,g}^{\mathrm{out}}\left({\Omega }_m\right)}{{\Sigma }_{j,k,g}{L}_{j,k}\left({\Omega }_m\right)})$

基于特征线方法，采用源迭代方法对每次计算得到的各栅元的中子通量密度φ_j,g进行数 值迭代，直到相邻两次迭代的和满足如下收敛条件：

$\left|\frac{{\phi }_{j,g}^{(k+1)}-{\phi }_{j,g}^{\left(k\right)}}{{\phi }_{j,g}^{(k+1)}}\right|<ϵ$

其中ε为收敛准则，通常ε＝10^-3~10^-4。最终当特征区域的所有栅元均满足如上收敛 条件后，得到的收敛结果即为屏蔽区域的中子通量分布结果。

本发明未详细阐述的部分属于本领域公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解 本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来 讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而 易见的，一切利用本发明构思的发明创造均在保护之列。