二维柱坐标系中子输运问题的边界型算法

A Boundary Type Algorithm for Neutron Transport Problems in Two-Dimensional Cylindrical Coordinate Systems

摘要: 传统的确定论方法在研究高维度的中子输运方程时需要求解大型矩阵的逆矩阵，该过程将增加计算量、降低计算效率，因此本文采用了一种新型的边界型方法——半边界法(the Half Boundary Method)来解决此类问题。半边界方法利用已有的中子输运方程推导得出相邻两点之间的中子通量密度关系通式，然后迭代地应用此通式求出模型内部任一点与边界点之间的中子通量密度关系，即可得到整个模型的中子通量密度分布。半边界方法在已知一半边界条件处的中子通量密度即可求得模型内部任一点的中子通量密度，提高了计算效率。本文将半边界法在真空边界条件模型中得出的结果与NECP-MCX软件得到的数值结果进行比较，通过分析发现不同空间变量和角度变量离散数量对半边界法数值结果的影响，当空间变量M = 20，N = 40和角度变量H = 50，K = 32时，半边界法的平均误差为0.76%，具有较高精确度。

Abstract: Traditional deterministic methods require solving the inverse matrix of large matrices when studying high-dimensional neutron transport equations, which increases computational complexity and reduces computational efficiency. Therefore, this paper adopts a new type of boundary type method—the Half Boundary Method to solve such problems. The half boundary method uses existing neutron transport equations to derive a general formula for the neutron flux density relationship between adjacent points, and then iteratively applies this formula to obtain the neutron flux density relationship between any point inside the model and the boundary point, thereby obtaining the neutron flux density distribution of the entire model. The half boundary method can obtain the neutron flux density at any point inside the model by knowing the neutron flux density at half of the boundary conditions, which improves computational efficiency. This article compares the results obtained by the half boundary method in the vacuum boundary condition model with the numerical results obtained by NECP-MCX software. Through analysis, it is found that the influence of different spatial and angular variables on the numerical results of the half-boundary method is significant. When the spatial variables M = 20, N = 40, and angular variables H = 50, K=32, the average error of the half boundary method is 0.76%, indicating high accuracy.

1. 引言

中子输运理论的数值解法研究是反应堆物理分析、安全控制等相关领域的关键组成部分，在科研基础研究和工程实际应用中发挥重要作用。目前，研究中子输运问题的方法主要分为非确定论方法和确定论方法。非确定论方法主要是基于统计理论的蒙特卡罗方法[1]，该方法被广泛应用于核工程领域，如探测器场、核反应堆等。然而，为了获得足够准确的结果，蒙特卡罗方法需要大量的计算样本和计算时间来消除其存在的统计涨落问题以提高计算精度。

离散纵坐标法(Sn) [2]和球谐函数法(Pn) [3]是现阶段常用的两种确定论方法，其思想都是通过对不同变量进行离散得到中子输运方程的数值解。Sn法[4]和Pn法相较于非确定方法中的蒙特卡罗方法具有更高的计算效率和更稳定的计算误差。然而这两种确定论方法在求解更精确解时各有局限，Sn法的数值结果大小受射线效应[5]和不同求积组的选取影响，Pn法在计算过程中需要求解复杂的角度耦合方程。近年来，人们提出了一种新的确定论方法，即半边界方法(HBM)。该方法已被广泛运用于导热[6]、熔化[7]、对流扩散[8]和中子输运问题中，其中包括核反应堆压力容器[9]和燃料棒上的导热熔化[10]。此外黄美团队等人将半边界法运用到了具有不同边界条件的二维稳态对流扩散方程中[11]。但现阶段半边界法在中子输运方程中的应用还很少，仅仅停留在简单的一维中子输运问题[12]，对于高维度的中子输运问题研究甚少。

因此，本文将从柱坐标系下的二维中子输运方程[13]出发，将半边界法用于求解真空边界条件下的中子输运问题。该方法主要是将基础的中子输运方程进行相应地空间和角度离散后经过一些近似处理得到模型相邻两点之间的递推通式，然后由递推关系式迭代推导出模型中内部点和边界点的关系。该方法将原本求解大型矩阵的逆矩阵过程转换为了求解边界点与内部点的关系过程，显著提高了计算效率和计算精度。为了验证半边界法在中子输运方程中的适用性和准确性，本文将半边界法得到的结果同NECP-MCX软件[14]得到的结果进行了比较分析，该分析表明，在选定合理的空间离散数目M、N和角度离散数目H、K后，半边界法具有极高的计算精度和计算效率。

2. 半边界法在圆柱坐标下的基本方程

2.1. 柱坐标系下的中子输运方程

非稳态情况下的中子输运方程是中子输运理论和反应堆物理分析的基础，其主要表达式为：

$\frac{1}{v}\frac{\partial \varphi }{\partial t}+\Omega \cdot \nabla \varphi +{\Sigma }_t\left(r,E\right)\varphi =S_0\left(r,E,\Omega ,t\right)+Q_s\left(r,E,\Omega ,t\right)+Q_f\left(r,E,\Omega ,t\right)$ (1)

式中v、t分别为中子运动速度和运动时间，$\varphi =\varphi \left(r,E,\Omega ,t\right)$ 表示中子角通量密度，${\Sigma }_t$ 为宏观截面，$S_0$ 为

固定中子源项，$Q_s={\displaystyle {\int }_0^{\infty }\text{d}{E}^{\prime }}{\displaystyle {\int }_{{\Omega }^{\prime }}{\Sigma }_s\left(r,E^{\prime },{\Omega }^{\prime }\to E,\Omega \right)\varphi \left(r,E^{\prime },{\Omega }^{\prime }\right)\text{d}{\Omega }^{\prime }}$ 为中子散射源项，$Q_f=\frac{\chi \left(E\right)}{4\pi }{\displaystyle {\int }_0^{\infty }\text{d}{E}^{\prime }}{\displaystyle {\int }_{4\pi }\nu {\Sigma }_f\left(r,E^{\prime }\right)\varphi \left(r,E^{\prime },{\Omega }^{\prime }\right)\text{d}{\Omega }^{\prime }}$ 为中子裂变源。

柱坐标系下的中子运动描述如图1所示，图中空间点的位置坐标由$r\left(r,\phi ,z\right)$ 确定，Ω为中子运动方向向量，$\theta$ 为Ω和$e_z$ 的夹角(其余弦$\epsilon =\cos \theta$ )，$\omega$ 为向量Ω在$e_{\theta }$ 和$e_r$ 所形成的平面的投影与向量$e_r$ 的夹角，$\eta$ 为向量Ω和向量$e_{\theta }$ 的余弦($\eta =\sin \omega \sqrt{1-{\xi }^2}$ )，$\mu$ 为向量Ω和向量$e_r$ 的余弦($\mu =\cos \omega \sqrt{1-{\xi }^2}$ )。因此柱坐标系下的一般中子输运方程可以表示为：

$\frac{1}{\nu }\frac{\partial \varphi }{\partial t}+\mu \frac{\partial \varphi }{\partial r}+\frac{\eta }{r}\frac{\partial \varphi }{\partial \phi }-\frac{\eta }{r}\frac{\partial \varphi }{\partial \omega }+\xi \frac{\partial \varphi }{\partial z}+{\Sigma }_t\left(r,E\right)\varphi =Q_s\left(r,E,\Omega ,t\right)+Q_f\left(r,E,\Omega ,t\right)+S_0\left(r,E,\Omega ,t\right)$ (2)

Figure 1. Description of neutron motion direction in cylindrical coordinate system

图1. 柱坐标系下中子运动方向的描述

多维中子输运方程表达式(1)中含有时间、空间等多个变量，求解起来耗时长、难度大。出于对低维中子输运方程在实际中运用多的考虑，以及为了更好地了解半边界法的精髓和适用性，本文将在稳态g能群和固定中子源的情况下，运用半边界法对二维圆柱对称模型进行中子输运方程的求解。对式(2)稍作修改即可，因此柱坐标下稳态、单能的二维中子输运问题可以写为：

$\left|\sin \theta \right|\cos \omega \frac{\partial {\phi }_g}{\partial r}-\frac{\left|\sin \theta \right|\sin \omega }{r}\frac{\partial {\phi }_g}{\partial \omega }+\cos \theta \frac{\partial {\phi }_g}{\partial z}+{\Sigma }_{g,t}\left(r,z\right){\phi }_g=S_g\left(r,z,\Omega \right)$ (3)

式(3)中，$\theta$ 和$\omega$ 为角度方向变量，r为空间变量，g为能群编号，本文将基于式(3)利用半边界法解决柱坐标下的中子输运问题。

2.2. 半边界法求解二维中子输运方程的推导柱坐标系下的离散模型

中子输运问题具有众多的变量，即使经过简化依然具有一定的求解难度。因此在求解中子输运问题之前，需要对二维圆柱模型进行在空间和角度上的变量离散。

Figure 2. Discrete schematic diagram of the model

图2. 模型的离散示意图

首先对空间变量r、z和角度变量$\omega$ 进行离散，二维圆柱的三个变量离散模型如图2所示，沿二维柱模型的半径方向以圆柱轴心为起点，柱面上的点为终点，离散为I个部分，共I + 1个点，第i点和第(i + 1)点间的距离为$\Delta r_i=r_{i+1}-r_i$ 。沿二维圆柱z轴方向，以下底面圆心为起始点，对空间变量z进行等距离离散，离散为J个部分，共J + 1个点，第j点和第(j + 1)点间的距离为$\Delta z_j=z_{j+1}-z_j$ 。在$\omega$ 角度方向由于圆柱几何的对称性，因此只需考虑离散一半角度，即$0^{\circ }<\omega <{180}^{\circ }$ 这部分。将$\omega$ 在$\left(0^{\circ },{180}^{\circ }\right)$ 内离散为M部分，离散步长$\Delta \omega ={\omega }_{m+1}-{\omega }_m$ 。

Figure 3. Discrete schematic diagram of angle $\theta$

图3. 角度$\theta$ 的离散示意图

$\theta$ 角度方向的离散范围为$0^{\circ }<\omega <{180}^{\circ }$ ，$\theta$ 方向离散模型如图3所示，将$\theta$ 在$\left(0^{\circ },{180}^{\circ }\right)$ 内离散为H部分，任意${\theta }_h$ 上的中子通量密度相互独立互不影响，离散步长$\Delta \theta ={\theta }_{h+1}-{\theta }_h$ ，这样某一$\left(r_i,z_j\right)$ 处任意$\Omega$ 方向都可由$\theta$ 和$\omega$ 描述。

2.3. 半边界法求解二维中子输运方程的推导

在求解本文中子输运问题式(3)时，为方便计算引入一个新的变量$V_g=\frac{\partial {\varphi }_g}{\partial w}$ 。当角度变量离散后，考虑某一$\left({\theta }_h,{\omega }_m\right)$ 方向，上式可写为：

$\{\begin{array}{l}\left|\sin {\theta }_h\right|\cos {\omega }_m\frac{\partial {\varphi }_{r,z}^{h,m}}{\partial r}-\frac{\left|\sin {\theta }_h\right|\sin {\omega }_m}{r}{V}_{r,z}^{h,m}+\cos {\theta }_h\frac{\partial {\phi }_{r,z}^{h,m}}{\partial z}+{\Sigma }_{g,t}\left(r,z\right){\varphi }_{r,z}^{h,m}=S_g^{h,m}\left(r,z,\Omega \right)\\ V_{r,z}^{h,m}=\frac{{\varphi }_{r,z}^{h,m+1/2}-{\varphi }_{r,z}^{h,m-1/2}}{\Delta w_m}\end{array}$ (4)

式(4)中，${\varphi }_{r,z}^{h,m}$ 表示在$\left(r,z\right)$ 处$\left({\theta }_h,{\omega }_m\right)$ 方向上的中子角通量密度；${\varphi }_{r,z}^{h,m\pm 1/2}$ 表示$\left(r,z\right)$ 点处$\left({\theta }_h,{\omega }_{m\pm 1/2}\right)$ 方向上的中子角通量密度，其中${\omega }_{m\pm 1/2}={\omega }_m\pm \frac{\Delta {\omega }_m}{2}$ ；而${\theta }_h$ 和${\omega }_m$ 为不同离散变量且各自独立互不影响，式(4)为圆柱坐标下描述本文所有中子运动方向的通式。

按上节2.2方法离散完二维圆柱中子输运模型后，在$\Delta r\times \Delta z$ 区域对式(4)中的r和z进行求积分，求解后得到以下方程：

$\{\begin{array}{l}\left|\sin {\theta }_h\right|\cos {\omega }_m\Delta z{r}_{i+0.5}\frac{\left({\varphi }_{i+1,j+1}^{h,m}+{\varphi }_{i+1,j}^{h,m}\right)-\left({\varphi }_{i,j+1}^{h,m}+{\varphi }_{i,j}^{h,m}\right)}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\Delta r\Delta z\frac{V_{i+1,j+1}^{h,m}+V_{i+1,j}^{h,m}+V_{i,j+1}^{h,m}+V_{i,j}^{h,m}}{4}\\ +\cos {\theta }_h{r}_{i+0.5}\Delta r\frac{\left({\varphi }_{i+1,j+1}^{h,m}+{\varphi }_{i,j+1}^{h,m}\right)-\left({\varphi }_{i+1,j}^{h,m}+{\varphi }_{i,j}^{h,m}\right)}{2}+r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z\frac{{\varphi }_{i+1,j+1}^{h,m}+{\varphi }_{i+1,j}^{h,m}+{\varphi }_{i,j+1}^{h,m}+{\varphi }_{i,j}^{h,m}}{4}\\ =r_{i+0.5}\Delta r\Delta z\overline{S_{i+0.5,j+0.5}^{h,m}}\\ V_{i/i+1,j/j+1}^{h,m}=\frac{{\varphi }_{i/i+1,j/j+1}^{h,m+1/2}-{\varphi }_{i/i+1,j/j+1}^{h,m-1/2}}{\Delta w_m}\end{array}$ (5)

其中$V_{i/i+1,j/j+1}^{h,m}$ 为$\left({\theta }_h,{\omega }_m\right)$ 方向上$\left(r_i,z_j\right)$ 处、$\left(r_{i+1},z_j\right)$ 处、$\left(r_i,z_{j+1}\right)$ 处或$\left(r_{i+1},z_{j+1}\right)$ 处变量V的值；${\Sigma }_{i+0.5,j+0.5}^{g,t}$ 为在节点$\left(i,j\right)$ 和$\left(i+1,j+1\right)$ 点间的平均宏观截面；$\overline{S_{i+0.5,j+0.5}^{h,m}}$ 为$\left({\theta }_h,{\omega }_m\right)$ 方向上$\left(i,j\right)$ 和$\left(i+1,j+1\right)$ 点间的平均固定中子源。

为进一步推导二维圆柱内部任一点与和边界点的中子通量密度通用关系式，引入“菱形”格式，对${\omega }_m$ 方向上的中子通量密度进行近似，使中子通量密度在$\Delta \omega$ 内呈线性变化：

${\varphi }_{r,z}^{h,m}\approx \frac{1}{2}\left({\phi }_{r,z}^{h,m+1/2}+{\varphi }_{r,z}^{h,m-1/2}\right)$ (6)

将式(6)带入式(5)第二个方程中消去${\varphi }_{r,z}^{h,m-1/2}$ 项，同时将式(5)第二个方程代入式(5)第一个方程，整理后得到$\left({\theta }_h,{\omega }_m\right)$ 方向上$\left(r_{i+1},z_{j+1}\right)$ 处和$\left(r_{i+1},z_j\right)$ 、$\left(r_i,z_{j+1}\right)$ 、$\left(r_i,z_j\right)$ 处的中子通量密度关系：

$\begin{array}{l}\left(\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}+\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}\right){\Phi }_{i+1,j+1}^{h,m}\\ =\left(-\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}-\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}\right){\Phi }_{i+1,j}^{h,m}\\ +\left(\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}-\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}-\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}\right){\Phi }_{i,j+1}^{h,m}\end{array}$

$\begin{array}{l}+\left(\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}-\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}\right){\Phi }_{i,j}^{h,m}\\ +r_{i+0.5}\Delta r\Delta z\overline{S_{i+0.5,j+0.5}^{h,m}}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}\left({\phi }_{i+1,j+1}^{h,m+1/2}+{\phi }_{i+1,j}^{h,m+1/2}+{\phi }_{i,j+1}^{h,m+1/2}+{\phi }_{i,j}^{h,m+1/2}\right)\end{array}$ (7)

将各个位置处的中子角通量密度和系数进行整理得到：

$\begin{array}{c}{\varphi }_{i+1,j+1,m}^h=X_{1ijm}{\varphi }_{i+1,j,m}^h+X_{2ijm}{\varphi }_{i,j+1,m}^h+X_{3ijm}{\varphi }_{i,j,m}^h+X_{4ijm}({\varphi }_{i+1,j+1,m+1/2}^h+{\varphi }_{i+1,j,m+1/2}^h\\ +{\varphi }_{i,j+1,m+1/2}^h+{\varphi }_{i,j,m+1/2}^h)+X_{5ijm},\{\begin{array}{l}i=0,1,2,\cdots ,I-1\\ j=0,1,2,\cdots ,J-1\\ m=1,2,3,\cdots ,K\end{array}\end{array}$ (8)

式中系数分别为：

$X_{1ijm}=\frac{-\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}-\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}+\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}$

$X_{2ijm}=\frac{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}-\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}-\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}+\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}$

$X_{3ijm}=\frac{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}-\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}-\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}+\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}$

$X_{4ijm}=\frac{\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}}{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}+\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}$

$X_{5ijm}=\frac{r_{i+0.5}\Delta r\Delta z\overline{S_{i+0.5,j+0.5}^{h,m}}}{\frac{\left|\sin {\theta }_h\right|\cos {\omega }_m{r}_{i+0.5}\Delta z}{2}+\left|\sin {\theta }_h\right|\sin {\omega }_m\frac{\Delta r\Delta z}{2\Delta {\omega }_m}+\frac{r_{i+0.5}\cos {\theta }_h\Delta r}{2}+\frac{r_{i+0.5}{\Sigma }_{i+0.5,j+0.5}^{g,t}\Delta r\Delta z}{4}}$

公式(8)为$\left({\theta }_h,{\omega }_m\right)$ 方向上任一节点$\left(r_{i+1},z_{j+1}\right)$ 与节点$\left(r_{i+1},z_j\right)$ 、$\left(r_i,z_{j+1}\right)$ 、$\left(r_i,z_j\right)$ 和$\left({\theta }_h,{\omega }_{m+1/2}\right)$ 方向上节点$\left(r_{i+1},z_{j+1}\right)$ 、$\left(r_{i+1},z_j\right)$ 、$\left(r_i,z_{j+1}\right)$ 、$\left(r_i,z_j\right)$ 之间的关系式。因此如果已知$\left({\theta }_h,{\omega }_{m+1/2}\right)$ 方向上4个节点和$\left({\theta }_h,{\omega }_m\right)$ 方向上任意3个节点，则可以通过以上通式推导出节点$\left(r_{i+1},z_{j+1}\right)$ 在$\left({\theta }_h,{\omega }_m\right)$ 方向上的中子通量。利用公式(8)重复推导不同节点之间的关系，则可以得到模型图2内部任一点与边界面BCE、边界面BDE、边界面ACE、边界面ADE、边界面BCF、边界面BDF、边界面ACF、边界面ADF内节点之间的关系，具体关系如下：

模型内部中子通量密度和边界面BCE中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{f=i}^{I-1}{a}_{i,j,m,f}{\varphi }_{f0k}^h}}+{\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{l=1}^j{b}_{i,j,m,l}{\varphi }_{Ilk}^h}}+{\displaystyle \sum _{l=1}^j{\displaystyle \sum _{f=i}^{I-1}{c}_{i,j,m,k}{\varphi }_{i,j,K}^h}}+q_{i,j,m}^h,\left(0\le i<I,0<j\le J,0\le m<K\right)$ (9)

模型内部中子通量密度和边界面BDE中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{f=1}^i{a}_{i,j,m,f}{\varphi }_{f0k}^h}}+{\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{l=1}^j{b}_{i,j,m,l}{\varphi }_{0lk}^h}}+{\displaystyle \sum _{l=1}^j{\displaystyle \sum _{f=1}^i{c}_{i,j,m,k}{\varphi }_{i,j,K}^h}}+q_{i,j,m}^h,\left(0<i\le I,0<j\le J,0\le m<K\right)$ (10)

模型内部中子通量密度和边界面ACE中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=0}^m{\displaystyle \sum _{f=i}^{I-1}{a}_{i,j,m,f}{\varphi }_{fJk}^h}}+{\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{l=j}^{J-1}{b}_{i,j,m,l}{\varphi }_{Ilk}^h}}+{\displaystyle \sum _{l=j}^{J-1}{\displaystyle \sum _{f=i}^{I-1}{c}_{i,j,m,k}{\varphi }_{i,j,K}^h}}+q_{i,j,m}^h,\left(0\le i<I,0\le j<J,0\le m<K\right)$ (11)

模型内部中子通量密度和边界面ADE中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{f=1}^i{a}_{i,j,m,f}{\varphi }_{fJk}^h}}+{\displaystyle \sum _{k=m}^{K-1}{\displaystyle \sum _{l=j}^{J-1}{b}_{i,j,m,l}{\varphi }_{0lk}^h}}+{\displaystyle \sum _{l=j}^{J-1}{\displaystyle \sum _{f=1}^i{c}_{i,j,m,k}{\varphi }_{i,j,K}^h}}+q_{i,j,m}^h,\left(0<i\le I,0\le j<J,0\le m<K\right)$ (12)

模型内部中子通量密度和边界面BCF中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=1}^m{\displaystyle \sum _{f=i}^{I-1}{a}_{i,j,m,f}{\varphi }_{f0k}^h}}+{\displaystyle \sum _{k=1}^m{\displaystyle \sum _{l=1}^j{b}_{i,j,m,l}{\varphi }_{Ilk}^h}}+{\displaystyle \sum _{l=1}^j{\displaystyle \sum _{f=i}^{I-1}{c}_{i,j,m,k}{\varphi }_{i,j,0}^h}}+q_{i,j,m}^h,\left(0\le i<I,0<j\le J,0<m\le K\right)$ (13)

模型内部中子通量密度和边界面BDF中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=1}^m{\displaystyle \sum _{f=1}^i{a}_{i,j,m,f}{\varphi }_{f0k}^h}}+{\displaystyle \sum _{k=1}^m{\displaystyle \sum _{l=1}^j{b}_{i,j,m,l}{\varphi }_{0lk}^h}}+{\displaystyle \sum _{l=1}^j{\displaystyle \sum _{f=1}^i{c}_{i,j,m,k}{\varphi }_{i,j,0}^h}}+q_{i,j,m}^h,\left(0<i\le I,0<j\le J,0<m\le K\right)$ (14)

模型内部中子通量密度和边界面ACF中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=1}^m{\displaystyle \sum _{f=i}^{I-1}{a}_{i,j,m,f}{\varphi }_{fJk}^h}}+{\displaystyle \sum _{k=1}^m{\displaystyle \sum _{l=j}^{J-1}{b}_{i,j,m,l}{\varphi }_{Ilk}^h}}+{\displaystyle \sum _{l=j}^{J-1}{\displaystyle \sum _{f=i}^{I-1}{c}_{i,j,m,k}{\varphi }_{i,j,0}^h}}+q_{i,j,m}^h,\left(0\le i<I,0\le j<J,0<m\le K\right)$ (15)

模型内部中子通量密度和边界面ADF中子通量密度的关系：

${\varphi }_{i,j,m}^h={\displaystyle \sum _{k=1}^m{\displaystyle \sum _{f=1}^i{a}_{i,j,m,f}{\varphi }_{fJk}^h}}+{\displaystyle \sum _{k=1}^m{\displaystyle \sum _{l=j}^{J-1}{b}_{i,j,m,l}{\varphi }_{0lk}^h}}+{\displaystyle \sum _{l=j}^{J-1}{\displaystyle \sum _{f=1}^i{c}_{i,j,m,k}{\varphi }_{i,j,0}^h}}+q_{i,j,m}^h,\left(0<i\le I,0\le j<J,0<m\le K\right)$ (16)

式(9)至式(16)中，i和j均为整数，$a_{i,j,m,f}$ ，$b_{i,j,m,l}$ ，$c_{i,j,m,k}$ ，$q_{i,j,m}^h$ 为方程系数。为了求解这些系数，本文以式(9)为例，将式(9)中各个位置的中子通量密度均写为式(8)的形式，而后代入式(8)即可得到系数表达式：

$\begin{array}{c}{a}_{i,j,m,f}=X_{1pqn}{a}_{i,j-1,m,f}+X_{2pqn}{a}_{i-1,j,m,f}+X_{3pqn}{a}_{i-1,j-1,m,f}+X_{4pqn}{a}_{i,j,m+1/2,f}+X_{5pqn}{a}_{i,j-1,m+1/2,f}\\ +X_{6pqn}{a}_{i-1,j,m+1/2,f}+X_{7pqn}{a}_{i-1,j-1,m+1/2,f},\left(0<i\le I,0<j\le J,0<m\le K\right)\end{array}$ (17)

$\begin{array}{c}{b}_{i,j,m,l}=X_{1pqn}{b}_{i,j-1,m,l}+X_{2pqn}{b}_{i-1,j,m,l}+X_{3pqn}{b}_{i-1,j-1,m,l}+X_{4pqn}{b}_{i,j,m+1/2,l}+X_{5pqn}{b}_{i,j-1,m+1/2,l}\\ +X_{6pqn}{b}_{i-1,j,m+1/2,l}+X_{7pqn}{b}_{i-1,j-1,m+1/2,l},\left(0<i\le I,0<j\le J,0<m\le K\right)\end{array}$ (18)

$\begin{array}{c}{c}_{i,j,m,k}=X_{1pqn}{c}_{i,j-1,m,k}+X_{2pqn}{c}_{i-1,j,m,k}+X_{3pqn}{c}_{i-1,j-1,m,k}+X_{4pqn}{c}_{i,j,m+1/2,k}+X_{5pqn}{c}_{i,j-1,m+1/2,k}\\ +X_{6pqn}{c}_{i-1,j,m+1/2,k}+X_{7pqn}{c}_{i-1,j-1,m+1/2,k},\left(0<i\le I,0<j\le J,0<m\le K\right)\end{array}$ (19)

$\begin{array}{c}{q}_{i,j,m}=X_{1pqn}{q}_{i,j-1,m}+X_{2pqn}{q}_{i-1,j,m}+X_{3pqn}{q}_{i-1,j-1,m}+X_{4pqn}{q}_{i,j,m+1/2}+X_{5pqn}{q}_{i,j-1,m+1/2}\\ +X_{6pqn}{q}_{i-1,j,m+1/2}+X_{7pqn}{q}_{i-1,j-1,m+1/2}+X_{8pqn},\left(0<i\le I,0<j\le J,0<m\le K\right)\end{array}$ (20)

式(17)到式(20)中，$p=i-1$ ，$q=j-1$ ，$n=m-1$ ；当且仅当$j=m=1$ ，$f=i$ 时$a_{i,j,m,f}=1$ ，否则$a_{i,j,m,f}=0$ ；当且仅当$i=m=1$ ，$l=j$ 时$b_{i,j,m,l}=1$ ，否则$b_{i,j,m,l}=0$ ；当且仅当$i=j=1$ ，$k=m$ 时$c_{i,j,m,f}=1$ ，否则$c_{i,j,m,f}=0$ ，同样地利用相同的方法即可以求解出其他式子中的系数。当在得到式(9)至式(16)的中子通量密度关系式后，根据不同方向下的条件选用相对应的中子通量密度关系式即可得到模型中任意位置相应方向下的中子通量密度。特别地，当${\omega }_m=\pi$ 和0时，式(8)与${\varphi }_{i/i+1,j/j+1,m+1/2}^h$ 项无关，因而当模型中为真空边界条件时，直接利用(8)式便可以得到$\left({\theta }_h,{\omega }_m\right)$ 方向上的中子通量密度分布情况。最后改变${\theta }_h$ 的取值，根据式(9)至式(16)便能得到任意$\left({\theta }_h,{\omega }_m\right)$ 方向上的中子通量密度，而后将所有方向下的中子通量密度值进行累加即可得到模型内任意位置的中子通量密度及模型内所有中子通量密度的分布情况。

3. 半边界法在二维圆柱坐标下的验证

对半边界法在二维圆柱坐标系下进行简单介绍和推导后，本文对真空边界条件的数值例子进行了测试，验证了半边界法在二维中子输运问题中的适用性。所有数值计算均在MATLAB R2021b软件中完成，运算过程中数值精度为双精度。对比分析过程以NECP-MCX软件计算结果作为基准解。

求解圆柱坐标系(r-z)下整个模型的二维稳态中子输运问题，如图所示的圆柱模型半径R = 1 cm，Z = 2 cm，模型内部具有固定的宏观截面(${\Sigma }_t=1{\text{cm}}^{-1}$ )和固定的均匀中子源($S=1\text{n}\cdot {\text{cm}}^{-2}\cdot {\text{s}}^{-1}$ )，模型整个外界面为真空边界条件(图4)。

Figure 4. Columnar symmetric model

图4. 柱对称模型

由模型外部为真空边界条件可知：

$\{\begin{array}{ll}{\varphi }_{r,Z,m}^h=0,0\le {\theta }_h<\frac{\pi }{2}\hfill & (\text{a})\hfill \\ {\varphi }_{r,0,m}^h=0,\frac{\pi }{2}\le {\theta }_h\le \pi \hfill & (\text{b})\hfill \\ {\varphi }_{R,z,m}^h=0,\frac{\pi }{2}\le {\omega }_m\le \pi \hfill & (\text{c})\hfill \end{array}$

特别地，当${\omega }_m=\pi$ 时，式(8)与${\varphi }_{i/i+1,j/j+1,m+1/2}^h$ 项无关，由已知的初始条件(b)和(c)可得到$\left({\theta }_h,{\omega }_m\right)$ 方向(其

中${\omega }_m=\pi$ )上任一点的中子通量，再将此结果带入式(9)便可以得到在某一$\left({\theta }_h,{\omega }_m\right)$ 方向(其中$\frac{\pi }{2}<{\omega }_m<\pi$ )

上任意$\left(r_i,z_j\right)$ 点(i从0到I，j从0到J)的中子角通量密度。之后通过圆柱轴心处的对称条件：

$\left(0,j,{\omega }_m\right)=\left(0,j,\pi -{\omega }_m\right),\left(0<{\omega }_m\le \frac{\pi }{2}\right)$

得到$\left(0,z_j\right)$ 处$\left({\theta }_h,{\omega }_m\right)$ 方向(其中$0\le {\omega }_m<\frac{\pi }{2}$ )上的中子角通量密度。在$\left({\theta }_h,{\omega }_m\right)$ 方向(其中${\omega }_m=\frac{\pi }{2}$ )上时，若模型中包含真空边界条件则利用“菱形”格式对${\omega }_m=\frac{\pi }{2}$ 方向上的中子角通量密度进行取值(即$2{\varphi }_{r,z,m}^h={\varphi }_{r,z,m+1/2}^h+{\varphi }_{r,z,m-1/2}^h$ ，这样便得到了$\left(0,z_j\right)$ 处$\left({\theta }_h,{\omega }_m\right)$ 方向(其中$0\le {\omega }_m<\frac{\pi }{2}$ )上的中子角通量密度。然后重复利用式(8)便可以得到所有$\left({\theta }_h,{\omega }_m\right)$ 方向上任意一点的中子通量密度分布，最后将$\left(r_i,z_j\right)$ 位置处所有${\omega }_m$ 方向(m为整数)的中子角通量密度进行累加便得到整个模型其中一个${\theta }_h$ 方向上中子通量密度。接着改变${\theta }_h$ ($0\le {\theta }_h<\frac{\pi }{2}$ )的值，并将$\left(r_i,z_j\right)$ 点上所有${\theta }_h$ ($0\le {\theta }_h<\frac{\pi }{2}$ )方向下的中子通量密度累加便得到模型一半的中子通量分布情况。

当${\theta }_h$ 取值在($\frac{\pi }{2}\le {\theta }_h\le \pi$ )时，初始条件则为(a)和(c)，可得到$\left({\theta }_h,{\omega }_m\right)$ 方向(其中${\omega }_m=\pi$ )上任一点的中子通量，再将此结果带入式(11)，接着利用上述同样的计算方法可得到$\left(r_i,z_j\right)$ 点上所有${\theta }_h$ ($\frac{\pi }{2}\le {\theta }_h\le \pi$ )方向下的中子通量密度，依次累加${\theta }_h$ ($\frac{\pi }{2}\le {\theta }_h\le \pi$ )方向上的中子通量密度，则可以得到模型另一半的中子通量密度，最后将所有${\theta }_h$ 方向上的中子通量密度加和，可得到整个二维圆柱模型的中子通量分布情况。

离散变量的敏感性分析

由2.2节可知，在使用半边界方法研究二维柱模型中的中子分布情况时，需要对中子输运方程中的空间变量R、Z及角度变量$\theta$ 、$\omega$ 进行数值离散，因此以上四个变量的离散数目对半边界方法的计算结果有重大影响。本小节将对半径R、高度Z和角度$\theta$ 、$\omega$ 的离散数目进行敏感性分析，并筛选出其中最佳的离散数目值。

先对角度变量$\theta$ 的离散数目H进行敏感性分析，设置空间变量R的离散份数M = 10，空间变量Z的离散分数N = 10，角度变量$\omega$ 的离散份数K = 50，例题中半边界法的结果如图5所示，不同H下的平均误差如图6所示。

Figure 5. Neutron flux distribution under different H conditions

图5. 不同H下的中子通量分布

Figure 6. Average error under different H values

图6. 不同H下的平均误差

通过图5可知，当H离散值超过50时，半边界方法计算结果和MCX计算结果十分贴合，此现象说明在真空条件下，半边界方法和MCX方法具有很高的准确性；由图6可知，在H由2增加至4时，半边界法的结果与MCX方法的结果的平均误差减少了69.24%，由4增加至20时，半边界法的结果与MCX方法的结果的平均误差减少了51.11%，继续增大H的值，半边界法的结果更加贴合MCX方法的结果，同时误差变化将非常微小，在H由50增加至60时平均误差只减少了0.56%。根据上述数据可知，在H = 50时误差最小，平均误差为1.09%。

综上所述，在M、N和K合理确定的情况下，若H ≥ 50半边界法的结果与MCX结果的误差将十分微小，因而半边界法的结果具有极高的准确性。

对角度变量$\omega$ 的离散份数K的敏感性分析时，设置半径R的离散数目M=10，高度Z的离散数目N = 10，角度变量$\theta$ 的离散份数H = 50，例题中半边界法的结果如图7所示，不同K下的平均误差如图8所示。由图7可知，在K离散数值越大，半边界法的结果与MCX方法的结果越贴合，但无明显的变化。

由图8可知，K的值由20增加至32的过程中半边界法的结果与MCX结果相对更加拟合，两种方法对比的平均误差减少了38.62%；当K ≥ 32后，随K值逐渐增大，MCX方法和半边界方法的结果的平均误差差别十分微小，此误差精度在0.00001内。

综上所述，在M、N和H的数值确定的情况下，增大K的值能明显提高半边界法的准确性，但当离散值超过适当范围，即使增加计算时间，也无法提高计算精度，因而在实际运用过程中需要根据具体的精度要求对K进行合理的取值。

对空间变量R的离散份数M进行敏感性分析时，设置高度Z的离散数目N = 10，角度变量$\theta$ 的离散份数H = 50，角度变量$\omega$ 的离散份数K = 32，计算结果如图9所示，不同M下的平均误差如图10所示。

由图9可知，半边界方法与MCX方法的计算结果十分相近；由图10可知，改变M的离散值对结果影响不大，且当M = 20时，半边界方法与MCX方法的平均误差最小，小的离散数值可有效提高计算时间。

综上所述，在N、H和K的值确定的情况下，增大M值对半边界法的计算结果影响不大，但在实际运用过程中仍需要根据具体的条件要求对M进行合理的取值。

对空间变量Z的离散份数N的敏感性进行分析时，设置半径R的离散数目M = 20，角度变量$\theta$ 的离散份数H = 50，角度变量$\omega$ 的离散份数K = 32，例题中半边界法的结果如图11所示，不同N下的平均误差如图12所示。

Figure 7. Neutron flux distribution under different K conditions

图7. 不同K下的中子通量分布图

Figure 8. Average error under different K values

图8. 不同K下的平均误差

由图10可知，在N < 40时半边界法的结果对N的取值敏感，在N ≥ 40时半边界法和蒙卡方法的结果十分吻合，若继续增大N的值，并不能明显提高半边界法的精度。由图12可知，N由20增加至40时半边界法的结果和MCX结果的误差相对减少最多，平均误差减少了8.38%。

由此可见，在M、H和K的值确定的情况下，选择合适的N值，HBM和MCX的结果误差将会很小且半边界法具有更快的计算效率，所以半边界方法具有拓展到带有反射边界条件的二维中子输运问题中去的能力和前景。

Figure 9. Neutron flux distribution under different M conditions

图9. 不同M下的中子通量分布图

Figure 10. Average error under different M values

图10. 不同M下的平均误差

4. 结论

本文基于半边界法，推导了中子输运方程在二维柱坐标系空间和角度方向上的离散过程，并把半边界法(HBM)应用于所提出的离散模型中，将其数值结果与NECP-MCX软件方法得到的结果进行了比较。与传统的离散纵坐标法不同，HBM得到的中子通量密度分布只与边界值有关，与其他空间位置的通量

Figure 11. Neutron flux distribution under different N conditions

图11. 不同N下的中子通量分布图

Figure 12. Average error under different N values

图12. 不同N下的平均误差

密度无关，显著减少了研究过程中的计算量。此外，在真空边界条件下的例子中，当M = 20、N = 40、H = 50和K = 32时HBM结果具有较高精度，平均误差为0.76%，从而验证了HBM在中子输运问题中的适用性。

基金项目

本研究得到了国家重点研发项目(2022YFB1902700)、国家教育部装备预研联合基金(8091B042203)、国家自然科学基金(11875129)、国家强脉冲辐射模拟与效应重点实验室基金(SKLIPR1810)、辐射应用创新中心基金(KFZC2020020402)、北京大学核物理与技术国家重点实验室基金(NPT2023KFY06)、中国铀业有限责任公司与华东理工大学核资源与环境国家重点实验室联合创新基金(2022NRE-LH-02)、中央高校基础研究基金(2023JG001)等支持。

NOTES

^*通讯作者。