0.   引　言

离散纵标方法（S_N）是确定论算法研究的主流方法之一，在反应堆物理分析和核武器研究中被广泛应用^[1-4]。S_N方法在计算孤立源、大空腔介质中子输运问题时，射线效应尤为明显^[5]，使得计算结果失真，这在一定程度上限制了应用。目前首次碰撞源法被广泛应用于缓解射线效应^[6-10]，该方法将中子注量率拆分为碰撞与未碰撞两部分, 使孤立源转换为全局源，进而缓解射线效应。首次碰撞源方法的核心是如何获得计算网格的未碰撞中子通量矩，目前常用的方法是网格中心法和网格角点平均值法。这2种方法计算网格中子通量矩均存在一定的偏差，使得未碰撞中子输运方程的中子守恒原则无法保证，导致计算结果失真。因此，本文采用全局因子修正未碰撞中子通量矩，确保未碰撞中子数守恒，以提高计算精度。

1.   计算原理

1.1   首次碰撞源方法

首次碰撞源方法的思想是将孤立源问题转化为全局源问题，从而缓解射线效应。以单群固定源问题为例，将中子角注量率分解为2部分，即：$ \psi {\text{(}}{\boldsymbol{r}}{\text{,}}{\boldsymbol{\varOmega }}{\text{)}} = {\psi _{\text{u}}}{\text{(}}{\boldsymbol{r}}{\text{,}}{\boldsymbol{\varOmega }}{\text{) + }}{\psi _{\text{c}}}{\text{(}}{\boldsymbol{r}}{\text{,}}{\boldsymbol{\varOmega }}{\text{)}} $，则单群中子输运方程可写成：

式中，$ {\boldsymbol{\varOmega }} $为中子运动方向；$ {\boldsymbol{r}} $为中子空间位置；n和m分别为散射截面和中子注量率的展开阶数；$ \varSigma _{\text{t}}^{}{\text{(}}{\boldsymbol{r}}{\text{)}} $为r处的宏观总截面，cm⁻¹；$ \varSigma _{{\text{s,}}n}^{}{\text{(}}{\boldsymbol{r}}{\text{)}} $为第n阶勒让德展开宏观散射截面，cm⁻¹；N为勒让德展开总阶数；${\phi _{n,m,c}}{\text{(}}{\boldsymbol{r}}{\text{)}}$为第(n,m)阶碰撞中子通量矩，cm⁻²·s⁻¹；$ {Y_{n{\text{,}}m}}{\text{(}}{\boldsymbol{\varOmega }}{\text{)}} $为第(n,m)阶球谐函数；$ s{\text{(}}{\boldsymbol{r}},{\boldsymbol{\varOmega }}{\text{)}} $为外源项，cm⁻³·s⁻¹；$ {\psi _{\text{u}}}{\text{(}}{\boldsymbol{r}}{\text{,}}{\boldsymbol{\varOmega }}{\text{)}} $为未碰撞中子角注量率， cm⁻²·s⁻¹；$ {\psi _{\text{c}}}{\text{(}}{\boldsymbol{r}}{\text{,}}{\boldsymbol{\varOmega }}{\text{)}} $为碰撞中子角注量率, cm⁻²·s⁻¹；$ {s_{{\text{fcs}}}}{\text{(}}{\boldsymbol{r}}{\text{,}}{\boldsymbol{\varOmega }}{\text{)}} $为首次碰撞源（是指中子在介质输运过程中，在 r 处单位体积内发生首次碰撞的中子数）, cm⁻³·s⁻¹； $ \nabla $为梯度算符。

首次碰撞源的表达式为：

式中，$ {\phi _{n{\text{,}}m{\text{,u}}}} $为第(n,m)阶未碰撞中子通量矩，其计算公式为：

$ {\psi _{\text{u}}}({\boldsymbol{r}},{\boldsymbol{\varOmega }}) $通过式（1）可获得：

式中，${{\boldsymbol{\varOmega}} _{p \to r}}$为位置r_p指向位置r的方向；$\tau ({\boldsymbol{r}},{{\boldsymbol{r}}_p})$为中子在几何空间中从位置 r_p穿行到位置r的光学距离，计算表达式为：

式中，l为中子沿飞行方向穿行的距离，cm。

基于上述分析可知，首次碰撞源方法的基本计算流程为：①利用式（6）计算源中子穿行的光学距离；②基于式（4）和式（5）求解不同位置处的未碰撞中子角注量率和未碰撞中子通量矩；③通过式（3）计算获得首次碰撞源，进而求解式（2）得到碰撞中子角注量率；④最后对碰撞中子角注量率和未碰撞中子角注量率求和计算实际的中子角注量率。

1.2   全局因子修正法

在首次碰撞源的计算中，其关键是计算网格的未碰撞中子通量矩。目前常采用网格中心法和网格角点平均值法计算离散网格的未碰撞中子通量矩。下面以二维问题为例对这2种方法进行原理说明，如图1所示，空间r_p位置有一各向同性、强度为S的中子点源。网格中心法计算任意r位置空间网格的未碰撞通量矩：先计算点源r_p与r位置空间网格中心O点之间的光学距离，然后通过式（5）求解中心O点的未碰撞中子角注量率，进而获得中心O点的首次碰撞源，用O点的首次碰撞源近似代替r位置空间网格的首次碰撞源；而网格角点平均值法计算任意r位置空间网格的未碰撞通量矩：先计算点源r_p与r位置空间网格4个角点(a,b,c,d)之间的光学距离，然后通过式（5）求解4个角点的未碰撞中子角注量率，进而获得各角点的未碰撞中子通量矩，用各角点的未碰撞中子通量矩算术平均值近似代替r位置空间网格的未碰撞中子通量矩来求解首次碰撞源。

图  1  二维问题示意图

Figure  1.  Diagram of Two-Dimensional Problem

下载: 全尺寸图片 幻灯片

网格中心法和网格角点平均值法计算得到的未碰撞中子角注量率均存在近似，使得计算获得未碰撞中子角注量率无法满足式（1），即未碰撞中子数不守恒。本文利用全局因子对网格角点平均值法计算获得的未碰撞中子注量率进行修正，使未碰撞中子数守恒，其基本原理如下：

对式（1）关于角度和空间进行积分可得：

式中，$ {R_{{\text{L,u}}}} $是未碰撞中子泄漏率， s⁻¹；$ {R_{{\text{T,u}}}} $是未碰撞中子总的反应率， s⁻¹；Q是外中子源发射率， s⁻¹。

按照物理意义，$ {R_{{\text{T,u}}}} $的表达式为：

式中，K是总的网格数；$ {\varSigma _{{\text{t}},n}} $表示第n个网格材料的宏观总反应截面，cm⁻¹；$ {V_n} $表示n个网格的体积，cm³；$ {\varPhi _{n,{\text{u}}}} $表示第n个网格的未碰撞中子注量率，cm⁻²·s⁻¹。

$ {R_{{\text{L,u}}}} $的计算需要对未碰撞中子流密度关于整个计算区域的外表面积分，本文采用7点Gauss-Legendre积分公式^[11]计算$ {R_{{\text{L,u}}}} $。

$ {R_{{\text{L,u}}}} $是基于非守恒网格角点平均值法中的未碰撞中子注量率计算得到的，由于计算过程中存在近似，式（7）的等号两边不相等。为了保证未碰撞中子数守恒，引入全局修正因子α：

式中，$ R_{{T} ,{\text{u}},{\text{new}}}^{} $是按照式（7）新算出来的未碰撞中子总反应率，即：

$ R_{{T} ,{\text{u}},{\text{new}}}^{} $保证了未碰撞中子数守恒，结合式（8）可得到修正后的未碰撞中子注量率为：

1.3   算法实现

JSNT-S是由中国工程物理研究院高性能数值模拟软件中心自主研发的三维离散纵标中子/光子输运软件^[12]。 JSNT-S利用S_N方法求解中子输运方程，能够获得三维问题内的中子和光子注量率分布。基于JSNT-S实现了全局因子修正首次碰撞源方法（FCS-C）的计算模块，其计算流程如图2所示，主要步骤为：①利用射线跟踪法计算粒子穿行的光学距离；②基于光学距离计算空间网格各角点的未碰撞中子角注量率，进而采用网格角点平均值法计算各空间网格的未碰撞中子注量率和未碰撞中子通量矩；③计算全局修正因子；④通过修正后的未碰撞中子注量率得到每个空间网格的首次碰撞源；⑤基于首次碰撞源计算碰撞中子注量率，进而得到总的中子注量率。

图  2  FCS-C的计算流程

Figure  2.  Calculation Process of FCS-C Based on Global Factor　　　　　

下载: 全尺寸图片 幻灯片

2.   数值分析

2.1   计算模型

采用日本Keisuke Kobayashi 教授提出的辐射输运计算基准题中的模型1进行数据分析^[13]，模型1从里到外依次由孤立源区、真空区和屏蔽区组成，文献[13]提供了详细的几何信息。表1给出了孤立源区的中子源强和各计算区域材料的宏观截面信息。计算过程中采用均匀网格，网格尺寸大小为2 cm。

表  1  截面参数

Table  1.  Section Parameters

区域	源强 /(cm−3·s−1)	宏观总截面/cm−1	宏观散射截面/ cm−1
源区	1	0.1	0.05
真空区	0	0.0001	0.00005
屏蔽区	0	0.1	0.05

下载: 导出CSV 

| 显示表格

2.2   计算结果与分析

基于Kobayashi提出的上述基准题进行数值分析，表2给出了S₈、S₁₂、S₁₆不同离散角度条件下JSNT-S直接的计算结果，可以发现随着求积组的离散角度数增加，计算结精度提高。图3给出了JSNT-S利用S₁₆，采用直接计算和本文的FCS-C计算的中子注量率分布，JSNT-S直接计算存在明显的射线效应（锯齿状），而采用FCS-C计算方法可很好地缓解射线效应。

表  2  不同离散角度下的 JSNT-S 直接计算结果

Table  2.  Calculation Results of JSNT-S at Different Discrete Angle

位置坐标/cm	中子注量率/(cm−2·s−1)				误差/%
	参考解	S8	S12	S16	S8	S12	S16
(5, 5, 5)	8.293×100	8.389×100	8.309×100	8.292×100	1.16	0.19	0
(15,15,15)	6.632×10−1	6.614×10−1	6.121×10−1	6.081×10−1	−0.27	−7.7	−8.32
(25,25,25)	2.688×10−1	2.087×10−1	3.169×10−1	2.905×10−1	−22.36	17.88	8.06
(35,35,35)	1.567×10−1	2.052×10−1	2.144×10−1	1.337×10−1	30.94	36.85	−14.66
(45,45,45)	1.044×10−1	2.008×10−1	1.373×10−1	9.749×10−2	92.29	31.5	−6.63
(55,55,55)	3.021×10−2	8.466×10−2	2.500×10−2	3.101×10−2	180.23	−17.24	2.64
(65,65,65)	4.066×10−3	1.515×10−2	2.534×10−3	4.686×10−3	272.65	−37.68	15.26
(75,75,75)	5.861×10−4	2.718×10−3	3.345×10−4	6.726×10−4	363.77	−42.94	14.75
(85,85,85)	8.661×10−5	4.753×10−4	4.935×10−5	9.518×10−5	448.76	−43.03	9.9
(95,95,95)	1.129×10−5	8.038×10−5	5.795×10−6	1.204×10−5	611.99	−48.68	6.66

下载: 导出CSV 

| 显示表格

图  3  中子注量率分布

Figure  3.  Neutron Flux Distribution

下载: 全尺寸图片 幻灯片

表3给出了各程序的中子注量率计算结果及其误差，其中参考解（Reference）为蒙特卡罗洛程序计算结果^[13]； FCS表示采用网格角点平均值法计算中子注量率的结果。由表3计算结果可以看出，FCS可以有效地缓解射线效应，将最大误差从15.26%降到6.15%， FCS-C进一步提高了计算精度，最大误差从6.15%降到3.71%。

表  3  JSWT-S与FCS、FCS-C中子注量率计算结果与误差

Table  3.  Calculation Results and Errors of Neutron Flux Rates for JSNT-S, FCS and FCS-C

位置坐标/cm	中子注量率/(cm−2·s−1)				误差/%
	参考解	JSNT-S	FCS	FCS-C	JSNT-S	FCS	FCS-C
(5, 5, 5)	8.293×100	8.292×100	8.132×100	8.274×100	0	−1.94	0.23
(15,15,15)	6.632×10−1	6.081×10−1	6.429×10−1	6.549×10−1	−8.32	−3.06	1.25
(25,25,25)	2.688×10−1	2.905×10−1	2.690×10−1	2.745×10−1	8.06	0.08	−2.13
(35,35,35)	1.567×10−1	1.337×10−1	1.470×10−1	1.509×10−1	−14.66	−6.15	3.71
(45,45,45)	1.044×10−1	9.749×10−2	1.016×10−1	1.042×10−1	−6.63	−2.7	0.14
(55,55,55)	3.021×10−2	3.101×10−2	2.877×10−2	2.949×10−2	2.64	−4.77	2.38
(65,65,65)	4.066×10−3	4.686×10−3	4.077×10−3	4.159×10−3	15.26	0.29	−2.3
(75,75,75)	5.861×10−4	6.726×10−4	5.872×10−4	5.985×10−4	14.75	0.18	−2.12
(85,85,85)	8.661×10−5	9.518×10−5	8.582×10−5	8.747×10−5	9.9	−0.91	−0.99
(95,95,95)	1.129×10−5	1.204×10−5	1.116×10−5	1.137×10−5	6.66	−1.18	−0.71

下载: 导出CSV 

| 显示表格

3.   结 论

S_N计算孤立源、大空腔问题存在射线效应。首次碰撞源方法能有效缓解射线效应，但在计算过程中，常采用非守恒网格角点平均值法计算未碰撞中子注量率，导致未碰撞中子数不守恒。针对这一问题，本文基于离散纵标程序JSNT-S实现了全局因子修正首次碰撞源方法，对未碰撞中子注量率进行修正，确保未碰撞中子输运方程满足中子守恒原则。利用Kobayashi输运基准题对全局因子修正的首次碰撞源模块进行了数值分析，计算结果表明首次碰撞源方法能有效地缓解射线效应，但局部存在较大偏差，采用全局因子修正方法对首次碰撞源方法进行修正后，计算精度进一步提高。