角度相关共轭函数的蒙卡自动减方差算法研究

Research on Monte Carlo Variance Reduction Algorithm for the Angle-Dependent Adjoint Function

摘要: 一致性共轭驱动重要抽样(CADIS)方法是目前广泛采用的自动减方差方法，该方法采用共轭标通量密度来一致性生成减方差参数，但缺乏角度偏倚信息，在处理角通量密度各向异性强的屏蔽问题时加速效果受到限制。本文基于CADIS-Ω方法，以多维离散纵标输运计算程序ARES为载体，开发角度相关减方差参数计算模块，构造正向加权重要性函数，自动生成角度相关源偏倚参数与权窗参数，弥补现有减方差参数在角度信息上的缺失，为蒙卡求解强各向异性深穿透屏蔽问题提供匹配度更高的自动减方差参数。采用IRI-TUB国际基准题进行验证。数值结果表明，CADIS-Ω方法能够在保证结果无偏的情况下获得较CADIS方法更优的加速效果。初步验证可说明CADIS-Ω方法对于复杂的强各向异性深穿透问题具有较好的应用前景。

Abstract: Consistent Adjoint Driven Importance Sampling (CADIS) is a widely utilized automatic variance reduction method. This method employs adjoint scalar flux to consistently generate variance reduction parameters. However, it lacks angle biasing information and its acceleration effect is limited when dealing with deep-penetration shielding problems characterized by strong anisotropy. Building upon the CADIS-Ω method and utilizing the multi-dimensional discrete ordinates transport calculation code ARES, this paper develops a calculation module for angle-dependent variance reduction parameters. It constructs a forward-weighted importance function and automatically generates angle-dependent source biasing parameters and weight window parameters to compensate for the absence of angle information in existing variance reduction parameters. This paper introduces a variance reduction parameter with a higher matching degree for solving the strong anisotropic deep penetration shielding problem by MCNP. The international benchmark problem IRI-TUB was used for verification purposes. Numerical results demonstrate that the CADIS-Ω method can achieve better acceleration effects compared to the CADIS method without biasing. Preliminary verification indicates that the CADIS-Ω method holds promising potential for addressing complex strong anisotropic deep penetration problems.

1. 引言

蒙特卡罗方法因误差来源单一，几何描述性强等优势被广泛用于辐射屏蔽计算，然而该方法的计算效率受到模型复杂度、模型尺寸与屏蔽层厚度等因素的限制。耦合式减方差方法结合离散纵标法(S_N)与蒙特卡罗方法(MCNP)的优势，成为降低蒙卡抽样误差，提高计算效率的重要手段。S_N共轭函数又称重要性函数，可用于衡量中子对目标物理量的相对贡献。利用S_N共轭计算获得重要性分布来指导蒙特卡罗抽样，能够有效提高蒙特卡罗抽样效率。国外早在1963年对重要性函数进行了研究，Kalos将共轭通量密度的近似解作为重要性函数，提出重要性抽样方法，在尺寸2400 cm的一维模型中取得了很好的加速效果[1]。随后数十年里，许多学者均将重要性抽样方法应用至不同的实际问题中，例如乏燃料桶的剂量计算[2]、测井计算[3] [4] [5]与真空管道问题的计算[6]。1997年，Haghighat与Wanger提出了一致性共轭驱动重要性抽样方法(CADIS) [7]，该方法使用确定论方法计算得到的共轭标通量密度，同时设置源偏倚参数和权窗参数，以保证源偏倚与输运偏倚之间的一致性。2007年，Peplow等人在Cooper方法[8]与CADIS方法基础之上，提出了FW-CADIS方法[9]，用于求解全局问题。2016年，Munk提出了一种针对强各向异性深穿透辐射屏蔽问题的减方差方法，并将其命名为CADIS-Ω方法。该方法利用正向加权的共轭标通量密度来生成减方差参数[10]。国内对共轭函数自动生成减方差参数的研究起步较晚，且大多数研究均基于CADIS方法。例如2016年，郑华庆[11]等人基于CADIS方法，提出了基于自动网格划分与权窗平滑的自适应减方差方法。2018年，刘聪[12]等人基于CADIS方法，自主开发了ARES共轭计算模块。2021年，郑琪[13]等人将CADIS方法与自动粗化网格方法结合，开发了MC耦合S_N程序NECP-MCX程序，2023年，西安交通大学大学的舒翰林等人将非结构网格方法与CADIS方法结合(UM-based CADIS) [14]，自动生成减方差参数加速蒙特卡罗计算。

综上，国内外众多学者对蒙特卡罗减方差方法进行了研究，大量研究工作均基于CADIS方法展开，该方法采用共轭标通量密度求解蒙卡减方差参数，然而针对强各向异性问题，角度依赖性非常强，采用仅与空间变量和能量变量的共轭标通量密度计算出的减方差参数并非最优参数，加速效果将受到限制。CADIS-Ω方法在共轭函数中引入角度信息，适用于求解强各向异性问题，然而以往的研究未考虑实际工程应用，且需熟知确定论共轭理论知识，手动设置共轭计算所需参数。

本文基于CADIS-Ω方法，对ARES共轭计算模块进行优化，开发角度相关减方差参数生成模块，构造正向加权重要性函数，自动生成角度相关减方差参数，弥补现有减方差参数在角度信息上的缺失，利用角度相关共轭函数求解角度相关减方差参数应用至IRI-TUB实际基准实验计算，对比加速效果。此外，对ARES共轭计算模块的优化可形成S_N正向–共轭计算一体化体系，减少用户对共轭理论知识的依赖，提高工作效率。

2. 理论基础

2.1. CADIS-Ω方法

在计算平均自由程大于10的深穿透屏蔽问题的减方差参数时，利用空间、角度与能量相关的共轭角通量密度需要消耗大量的计算机内存[10]。同时S_N方法在数值离散过程中对角度变量仅采用了有限个离散方向，无法得到中子通量密度在角度上的连续分布[12]。在CADIS方法中，常采用与空间和能量变量相关的共轭函数进行计算：

${\varphi }^{†}\left(\overrightarrow{r},E\right)={\displaystyle {\int }_{\Omega }{\psi }^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\text{d}\overrightarrow{\Omega }}$(1)

式中，${\psi }^{†}\left(\overrightarrow{r},E,\widehat{\Omega }\right)$ 为共轭角通量密度，${\varphi }^{†}\left(\overrightarrow{r},E\right)$ 为共轭标通量密度。然而在通量密度的角度依赖性非常

显著的深穿透问题中，尤其是涉及到强中子流、孔道问题时，仅与空间能量相关的共轭函数产生的加速效果并不显著。因此，CADIS-Ω方法引入角度信息，采用角度相关共轭函数来一致性的偏倚蒙特卡罗输运计算，角度相关共轭函数通过以下式子计算：

${\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)=\frac{{\displaystyle {\int }_{\Omega }{\psi }^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\text{d}\overrightarrow{\Omega }}}{{\displaystyle {\int }_{\Omega }\psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\text{d}\overrightarrow{\Omega }}}$(2)

式中，${\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)$ 为角度相关共轭标通量，${\psi }^{†}\left(\overrightarrow{r},E,\widehat{\Omega }\right)$ 为共轭角通量密度，$\psi \left(\overrightarrow{r},E,\widehat{\Omega }\right)$ 为正向角通量密度，共轭角通量密度与正向角通量密度的乘积称为贡献通量密度。因此角度相关共轭函数也称正向加权的共轭通量密度，此时对应的角度相关源偏倚参数为：

${\widehat{q}}_{\Omega }=\frac{{\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)q\left(\overrightarrow{r},E\right)}{{\displaystyle \int {\displaystyle \int {\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)q\left(\overrightarrow{r},E\right)\text{d}E\text{d}\overrightarrow{r}}}}=\frac{{\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)q\left(\overrightarrow{r},E\right)}{R_{\Omega }}$(3)

粒子的偏倚权重也需按：

$W_{\Omega }\left(p\right)=\frac{R_{\Omega }}{{\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)}$ (4)

进行修正，保证结果的无偏性。

2.1.1. 源偏倚

源偏倚是蒙特卡罗计算中常用的减方差技巧，是针对源粒子发射方向、能量和位置的偏倚[15]。为提高抽样效率，可以通过增加对更关注的方向、能量段或位置的抽样来实现。源方向偏倚主要通过调整粒子发射方向的概率密度函数，使粒子朝着更加关注的方向输运；源能量偏倚通过改变粒子在不同能量段的抽样概率，增加粒子在更关注的能量区域的抽样，从而提高抽样概率；而源位置偏倚通过多次抽样更关注的位置来提高抽样效率。

在MCNP5计算中，常采用SEDF卡来定义源粒子抽样信息。包括利用SI、SP与SB卡来分别定义相关信息。SI卡定义源粒子抽样区间，包括源区尺寸范围与源能量区间，SP卡定义初始源抽样概率，SB卡定义偏倚抽样概率。在CADIS-Ω方法中，角度相关共轭通量密度${\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)$ 是一个与网格–能量变量相关的函数，涉及到几何网格信息与能量区间信息，因此SI卡与SP卡中填入的信息均需要与S_N计算中的对应。在SB卡中填入偏倚抽样概率时，首先需计算探测器响应值$R_{\Omega }$ ，即：

$R_{\Omega }={\displaystyle {\int }_r{\displaystyle {\int }_E{\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)q\left(\overrightarrow{r},E\right)\text{d}E\text{d}\overrightarrow{r}}}$(5)

随后按照：

$\begin{array}{l}{\widehat{q}}_{\Omega }\left(\overrightarrow{r}\right)={\displaystyle \int q_{\Omega }\left(\overrightarrow{r},E\right){\varphi }_{\Omega }^{†}}\left(\overrightarrow{r},E\right)\text{d}E/R_{\Omega }\\ {\widehat{q}}_{\Omega }\left(E\right)={\displaystyle \int q_{\Omega }\left(\overrightarrow{r},E\right){\varphi }_{\Omega }^{†}}\left(\overrightarrow{r},E\right)\text{d}V/R_{\Omega }\end{array}$(6)

分别计算关于空间和能量上的偏倚分布。

2.1.2. 权窗

权窗是一种有效减方差方法，结合了分裂与轮盘赌技巧，用于控制与空间能量变量相关的粒子权重。其主要目标是将所有抽样粒子的权重限制在特定的权窗上限值与权窗下限值范围内。如图1所表示的权窗原理图，当粒子的权重位于权窗上限值之上时，触发分裂操作，高权重的粒子将分裂为多个权重低权重的粒子；当粒子的权重位于权窗下限值之下时，执行轮盘赌操作，较低权重的粒子将会被杀死或以权窗范围内的权重存活下来[16]。这一过程有助于有效地控制粒子权重的分布，减小方差并提高蒙特卡罗计算效率。基于CADIS-Ω方法，权窗下限值的具体计算公式为：

$W_1\left(\overrightarrow{r},E\right)=\frac{1}{S}\frac{R_{\Omega }}{{\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)}\frac{2}{C_u+1}$ (7)

式中，$C_u$ 为权窗上下限值之比，一般为5，S为S_N正向计算总源强值，按$S={\displaystyle {\int }_0^{\infty }\text{d}E}{\displaystyle {\int }_{v}Q\left(\overrightarrow{r},E\right)\text{d}V}$进行计算，除以总源强的目的在于将权窗下限值归一化，这是由于在使用MCNP5进行输运计算时，源强值设为1.0 n∙s⁻¹。

Figure 1. Diagram of the weight window

图1. 权窗原理示意图

在MCNP5计算中，采用WWP卡定义权窗相关的参数，包括权窗文件读取位置、权窗上下限之比和粒子分裂数等其他相关变量。通过CADIS-Ω方法，ARES程序能够自动生成与网格–能量变量相关的权窗参数，包括网格尺寸信息，能群划分结构与每个能群下每一个网格对应的权窗下限值，这些参数独立储存在名为WWINP的权窗文件中，具体的数据格式可参考文献[17]。通过读取WWP卡中定义的权窗文件位置读取信息，实现网格能量相关的权窗下限值的读取，用于限制模型不同位置处的粒子权重的分布，指导MCNP5抽样。

2.2. 共轭函数

共轭函数又称中子价值函数或重要性函数，可用于测量某一空间位置的粒子对响应值的贡献[18]。为了得到正确反映粒子重要性分布的中子价值函数，需要求解共轭中子输运方程。

在中子输运理论中，针对稳态固定源问题，中子输运方程可表示为：

$\begin{array}{l}\overrightarrow{\Omega }\cdot \nabla \psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)+{\Sigma }_t\left(\overrightarrow{r},E\right)\psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\\ ={\displaystyle \underset{4\pi }{\int }{\displaystyle {\int }_0^{\infty }{\Sigma }_s\left(E\text{'}\to E,\overrightarrow{\Omega }\text{'}\to \overrightarrow{\Omega }\right)\psi \left(\overrightarrow{r},E\text{'},\overrightarrow{\Omega }\text{'}\right)\text{d}E\text{'}\text{d}\overrightarrow{\Omega }\text{'}+q_e\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)}}\end{array}$(8)

式中，$\overrightarrow{r}$ 、E与$\overrightarrow{\Omega }$ 分别为空间变量、能量变量与角度变量，$\psi$ 为角通量密度，${\Sigma }_t$ 为宏观总截面，${\Sigma }_s$ 为宏观散射截面，$q_e$ 为外源。方程等号右边的产生项仅包括散射源项与外中子源，不包括核裂变引起的产生率。

定义算子L与$L{}^{†}$ ，设函数及定义于域$G\left(P\right)$ 内，其中P为相空间上$\overrightarrow{r}\times E\times \overrightarrow{\Omega }$ 的点，将函数$\psi \left(P\right)$ 与${\psi }^{†}\left(P\right)$ 做内积：

$\left(\psi ,{\psi }^{†}\right)={\displaystyle {\int }_G\psi {\psi }^{†}\text{d}P}$ (9)

若$L\psi$ 及$L^{†}{\psi }^{†}$ 满足：

$\left({\psi }^{†},L\psi \right)=\left(\psi ,L^{†}{\psi }^{†}\right)$ (10)

则称$L{}^{†}$ 为L的共轭算子。

根据价值理论可推到出共轭中子输运方程，具体推导过程见文献[19]。稳态共轭中子输运方程可表示为：

$\begin{array}{l}-\overrightarrow{\Omega }\cdot \nabla {\psi }^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)+{\Sigma }_t\left(\overrightarrow{r},E\right){\psi }^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\\ ={\displaystyle \underset{4\pi }{\int }{\displaystyle {\int }_0^{\infty }{\Sigma }_s\left(E\to E\text{'},\overrightarrow{\Omega }\to \overrightarrow{\Omega }\text{'}\right){\psi }^{†}\left(\overrightarrow{r},E\text{'},\overrightarrow{\Omega }\text{'}\right)\text{d}E\text{'}\text{d}\overrightarrow{\Omega }\text{'}+q_e^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)}}\end{array}$(11)

式中，${\psi }^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)$ 为共轭角通量密度，$q_e^{†}$ 为共轭源。由共轭中子输运方程可知，粒子是从能量E向能量$E\text{'}$ 输运，角度从$\overrightarrow{\Omega }$ 至$\overrightarrow{\Omega }\text{'}$ ，即能量发生上散射、角度变化与正向计算相反，固定源也从$q_e$ 变为$q_e^{†}$ 。

在求解式(11)时，需将空间变量，角度变量与能量变量分别进行离散，确定边界条件以及初始值。共轭源项是求解共轭输运方程时需正确定义的初始条件，其数值大小将直接影响计算出的重要性函数，从而导致减方差参数的变化，影响减方差方法的加速效果。在大多数蒙特卡罗粒子输运问题中，目标是计算特定位置的响应值(例如反应率、通量密度或剂量率)。当目标函数为反应率时，可用以下式子表示：

$R_{tally}=\langle \psi ,{\Sigma }_d\rangle ={\displaystyle {\int }_v{\displaystyle {\int }_E{\displaystyle {\int }_{\Omega }\psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right){\Sigma }_d\left(\overrightarrow{r},E\right)\text{d}\overrightarrow{\Omega }\text{d}E\text{d}V}}}$(12)

式中，V为探测器体积，E表示能量，$\overrightarrow{\Omega }$ 表示角度，$R_{tally}$ 表示探测器响应值，$\psi$ 为中子角通量密度，${\Sigma }_d$ 为探测器响应截面。由等式$\langle \psi ,q^{†}\rangle =\langle {\psi }^{†},q\rangle$ 可得，共轭源项此时应设为探测器响应截面${\Sigma }_d$ ，探测器响应值也可由共轭角通量密度与正向源分布积分获得：

$R_{tally}={\displaystyle {\int }_v{\displaystyle {\int }_E{\displaystyle {\int }_{\Omega }{\psi }^{†}\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)q\left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\text{d}\overrightarrow{\Omega }\text{d}E\text{d}V}}}$(13)

由此可得，探测器响应值$R_{tally}$ 就是目标函数，${\psi }^{†}$ 表示粒子对探测器响应(即目标函数)的贡献。

当目标函数为中子通量密度时，即：

$R_{tally}=\varphi \left(\overrightarrow{r},E\right)={\displaystyle {\int }_{\Omega }\psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\text{d}\overrightarrow{\Omega }}$(14)

则对应的共轭源项为$q^{†}=总能群数$ 。

当目标函数为剂量当量率时，则目标函数为：

$R_{tally}=\varphi \left(\overrightarrow{r},E\right)={\displaystyle {\int }_{\Omega }\psi \left(\overrightarrow{r},E,\overrightarrow{\Omega }\right)\text{d}\overrightarrow{\Omega }}$(15)

式中，$C_g$ 为通量–剂量转换系数。则对应的共轭源项为$q^{†}=C_g$ 。

3. 自动减方差参数计算流程

角度相关减方差参数计算模块的开发和CADIS-Ω方法实现均基于ARES程序。ARES程序[20] [21]是由华北电力大学IRPS团队基于离散纵标法自主开发的程序系统，具备与蒙特卡罗程序接续计算和耦合加速功能，可进行中子/光子耦合问题的输运求解，能够求解任意阶各向异性散射的固定源问题，具备复杂核装置屏蔽计算分析能力。角度相关减方差参数的具体计算过程如图2所示。第一步需要进行确定论正向计算来存储正向角通量密度。第二步进行ARES共轭输运计算，获得共轭角通量密度。根据式(2)得到正向加权的共轭标通量密度${\varphi }_{\Omega }^{†}\left(\overrightarrow{r},E\right)$ 。最后，基于CADIS-Ω方法，通过角度相关减方差参数计算模块，自动生成源偏倚参数与权窗参数，指导MCNP5的计算。

Figure 2. Implementation of the angle-dependent variance reduction parameters

图2. 自动减方差参数计算流程

4. 数值验证

针对不同的屏蔽问题，不同的减方差技巧可能带来不同的加速效果，使用不合适的减方差方法甚至会降低蒙特卡罗计算效率。为衡量某种减方差技巧对MC计算效率的影响，常采用品质因子(Figure Of Merit, FOM)来评价[22]，即：

$\text{FOM}=\frac{1}{R^{2}T}$ (16)

式中，R为相对误差，T为计算时间。由上述式子可知，FOM值的大小与R²T成反比，在较短的时间内获得较低的相对误差，FOM值就越大，MC的计算效率就越高。

4.1. IRI-TUB直孔道实验

IRI-TUB研究反应堆位于布达佩斯技术大学核技术系，功率为100千瓦，常用于研究粒子在直孔道或弯曲孔道内的输运行为并验证粒子输运数值计算方法[23] [24]。IRI-TUB模型为典型的强各向异性相关问题，其管道内存在强各向异性通量密度并且石墨材料中¹²C核素存在强各向异性散射现象。IRI-TUB装置X-Y平面的几何模型如图3所示，反应堆模型由一个堆芯和一个大型辐射通道组成。堆芯活性区域由24根尺寸为7.20 cm × 7.20 cm燃料棒组成，区域周围为反射层，由石墨层和水组成。堆芯外的屏蔽层由金属铝组成。在距离活性区域外25.00 cm处有一个用混凝土填满的大型辐射通道，内有一个半径为11.80 cm、长度为187.00 cm的圆柱形孔道。管道外侧包裹着一层4.50 mm厚的不锈钢。探测器沿中轴线设置，探测点距离管道入口分别为0.00、67.00、121.00、148.00和175.00 cm。选取管道出口处的探测点测量¹¹⁵In(n, n′)^115mIn反应率，具体探测点如图3所示。此外，测量该点在1.10 × 10⁻¹⁰ MeV~1.96 × 10¹ MeV能量范围内的中子能谱，用于评估CADIS-Ω方法的无偏性。

Figure 3. Configuration of the IRI-TUB straight problem at Z = 0 cm

图3. IRI-TUB直孔道X-Y面几何模型

采用ARES程序进行共轭输运计算时，探测点所在网格定义为共轭源区，探测器反应截面设为共轭源强值，其值大小为4.89 n∙cm⁻³∙s⁻¹。求解过程中能群采用199群并47群，在满足计算精度的条件下节省相应的计算时间，各向异性散射阶数为P₃阶展开，求积组采用定向θ权重差分格式(EDW)，网格大小约为3厘米。确定论三维几何模型如图4所示。由于计算模型中包含187 cm长的通道，采用首次碰撞源方法减缓S_N计算时产生的射线效应，从而提高计算精度。图5显示了47群能群结构下的归一化正向初始源强分布和偏向源强分布。从归一化的源强分布可以看出，MCNP将增加在0.30 MeV~20.00 MeV能量范围内的源粒子抽样。原因在于¹¹⁵In(n, n′)^115mIn为阈值探测器，高能中子对探测器的响应产生的影响更大，加大对该区域的偏倚抽样能够适当的提高MCNP5的计算效率。

Figure 4. Geometry of the IRI-TUB benchmark problem with straight duct

图4. IRI-TUB直孔道三维可视化几何图

Figure 5. Normalized distribution of unbiased and biased sources of the IRI-TUB benchmark problem with straight duct

图5. IRI-TUB直孔道问题归一化的初始源分布与偏倚源分布

表1列出了孔道出口点(43.2, 248.5, 0.5)的¹¹⁵In(n, n′)^115mIn反应率计算结果，包括MCNP直接模拟、CADIS方法加速模拟和CADIS-Ω方法加速模拟三种抽样模拟方式。基准实验在该点的测量值为4.61E−18s⁻¹，表中C/E是指计算结果与基准值结果的比值，文献[25]中提供了直孔道模型中关键点位置的反应率基准值和采用蒙特卡罗程序MORSE-SGC/S接续二维离散纵标程序DOT3.5计算出的反应率计算结果。从表中可以看出，为确保获得可靠的统计结果，在不采用任何减方差方法的情况下，MCNP直接模拟了10天以上，而采用CADIS方法和CADIS-Ω方法后，在不到7个小时的时间便获得了较为可靠的结果。根据FOM值，与直接模拟相比，CADIS-Ω方法将蒙特卡罗计算效率提高了4个数量级。同时，对比CADIS与CADIS-Ω两种方法计算出的FOM值，结果可知，在相同计算时间内，CADIS-Ω方法产生的FOM值更高，可说明在相同计算时间内，CADIS-Ω方法能够产生更低的统计误差R，从而获得更高的FOM值。对比CADIS方法，CADIS-Ω产生的FOM值约提高了38%。因此，可说明CADIS-Ω方法更适合处理强各向异性深穿透屏蔽问题。

Table 1. ¹¹⁵In(n, n′)^115mIn reaction results and FOMs for the IRI-TUB straight duct problem

表1. 点(43.2, 248.5, 0.5)处¹¹⁵In(n, n′)^115mIn反应率计算结果

计算方式	计算结果(s^−1)	FOM	C/E	计算时间
MCNP	3.06E−18	0.71	0.66	10.0 d
CADIS	3.09E−18	5232.00	0.67	6.5 h
CADIS-Ω	3.09E−18	7122.00	0.67	6.5 h

为了验证CADIS-Ω方法是否能与CADIS方法一样保证统计结果无偏，求解该点(43.2, 248.5, 0.5)在1.0 × 10⁻⁸ MeV~1.0 × 10¹ MeV范围内的中子能谱。图6(a)比较了文献参考值、MCNP直接模拟和CADIS-Ω加速模拟的能谱计算结果。图中可看出，MCNP直接模拟与加速模拟在低能区的结果偏大，整体趋势一致，可以表明三者具有良好的一致性。图6(b)显示了各能量段的统计误差，其中MCNP直接模拟计算时间为225小时(约9天)，CADIS-Ω加速模拟时间为29小时。图中可看出，在MCNP直接模拟计算中，尽管模拟了数天，不同能量段的相对统计误差仍然存在明显差异，范围在0.01~0.08，波动较大。然而，利用CADIS-Ω方法后可将所有能量段的统计误差保持在0.02以下，显著减少不同能量段统计误差之间的差异。通过对比，验证了CADIS-Ω方法能确保无偏的统计结果，并表现出很高的可靠性。

(a) (b)

Figure 6. Neutron spectrum and statistical tally errors for IRI-TUB straight ductproblem: (a) Neutron spectrum calculation results; (b) Statistical tally errors for each energy segment

图6. IRI-TUB直孔道问题点(43.2, 248.5, 0.5)的中子能谱结果与统计误差：(a) 中子能谱计算结果；(b) 各能量段的统计误差

4.2. IRI-TUB弯曲孔道实验

IRI-TUB弯曲孔道模型如图7所示，孔道前端的堆芯区域、模型中的材料信息均与直管道问题中的保持一致，明显的不同之处在于堆芯后部的大辐照坑内存在一个90度的直角弯曲孔道。探测点也相应地发生了变化，沿中轴线分布在距离管道入口0.00、67.00、134.00、161.00和188.00 cm处。选取第4个关键点，测量⁵⁵Mn(n, γ)⁵⁶Mn(无Cd包层)反应率，以验证CADIS-Ω方法的加速效果。

ARES共轭计算中网格总数为26 × 79 × 27，采用定向θ权重差分格式(EDW)求积组，迭代收敛准则为10⁻⁴。问题采用P₃阶各向异性散射和47能群结构进行求解。共轭源区对应探测点(43.2, 234.5, 27.5)所在的几何区域，源强大小为78.70 n∙cm⁻³∙s⁻¹，为47群下⁵⁵Mn(n, γ)⁵⁶Mn反应截面的总和。图8展示了ARES输运计算时的三维可视化图。此外，图9显示了47能群结构下归一化的初始源强分布和偏倚源强分布。与¹¹⁵In(n, n′)^115mIn探测器相比，源粒子偏倚的趋势并不是特别明显，原因在于⁵⁵Mn(n, γ)⁵⁶Mn作为全能区探测器，其响应截面在所有能群上都存在数值，所有能量范围内的中子对响应值均存在一定的贡献。

Figure 7. Configuration of the IRI-TUB bend duct problem at X = 43.2 cm

图7. IRI-TUB弯曲孔道Y-Z面几何模型

Figure 8. Geometry of the IRI-TUB benchmark problem with bend duct

图8. IRI-TUB弯曲孔道三维可视化几何图

Figure 9. Normalized distribution of unbiased and biased sources of the IRI-TUB benchmark problem with bend duct

图9. IRI-TUB弯曲孔道归一化的初始源分布与偏倚源分布

表2中给出了点(43.2, 234.5, 27.5)处的⁵⁵Mn(n, γ)⁵⁶Mn (无Cd包层)反应率计算结果，包括直接模拟、CADIS加速模拟与CADIS-Ω加速模拟三种模拟方式的计算结果，参考文献中在该点的测量值为1.59 × 10⁻¹⁷ s⁻¹。从表中可以看出，使用MCNP5直接模拟需要24天以上才能获得较为可信的统计结果。然而，在采用减方差方法后，统计误差可在24小时内减小至1%以下。表中的FOM值表明，在确保结果无偏的情况下，CADIS-Ω方法可将MCNP5计算效率提高三个数量级，但加速效果与CADIS方法一致。

Table 2. ⁵⁵Mn(n, γ)⁵⁶Mn reaction results and FOMs for the IRI-TUB bend duct problem

表2. 点(43.2, 234.5, 27.5)处⁵⁵Mn(n, γ)⁵⁶Mn反应率计算结果

计算方式	计算结果(s^−1)	FOM	C/E	计算时间
MCNP	2.32E−17	0.01	1.46	576.00 h
CADIS	2.41E−17	12.00	1.51	19.50 h
CADIS-Ω	2.39E−17	12.00	1.50	19.50 h

5. 结论

基于CADIS-Ω方法，开发角度相关减方差参数生成模块，对ARES共轭计算模块进行优化。利用正向加权的共轭函数自动生成蒙特卡罗减方差参数，提高蒙特卡罗计算效率。利用SINBAD数据库中的IRI-TUB基准实验的直孔道模型与90度弯曲孔道模型，验证了CADIS-Ω方法的正确性与适用性。对于IRI-TUB直孔道问题，反应率计算结果表明，与直接模拟相比，CADIS-Ω方法可将计算效率提高4个量级，且加速效果对比CADIS方法提高了约38%。对于IRI-TUB弯曲模型，CADIS-Ω方法可将蒙特卡罗计算效率提高3个量级，但其加速效果与CADIS方法一致。虽然本研究初步证明了CADIS-Ω方法更为优越，但对于弯曲模型其生成的加速结果与CADIS一致的原因还需要进一步研究。同时，模型的差异性以及引起强各向异性的原因多种多样，未来的重点应对强各向异性问题分类，详细评估CADIS-Ω方法的加速效果与适用性。

NOTES

^*通讯作者。