A Nodal-iterative Technique for Criticality Calculations in Multigroup Neutron Diffusion Models

ABSTRACT

In this work, a nodal and iterative technique to evaluate the effective multiplication factor as well as the neutron flux, in multigroup diffusion problems, is presented. An iterative scheme, similar to the source iteration method, is implemented to decouple the system of differential equations which is the fundamental mathematical model. Then, analytical solutions are derived for the one-dimensional transverse integrated equations, of each energy group, resulting from a nodal approach. Constant approximations are assumed for the unknown transverse leakage terms in the contours of the nodes. In addition, constant and linear representations are investigated to express the fluxes in the source term to be updated in the iterative process. Numerical results for the effective multiplication factor were obtained for a series of twodimensional multigroup problems with upscattering and downscattering. The procedure is simple, fast, the analysis of the results indicated a satisfactory agreement with results available in the literature and the use of different approximations to the source term seems to be a good alternative, instead of using higher-order approximations on the contour of the nodes, to improve accuracy.

Keywords:
multigroup neutron diffusion equation; source iteration method; nodal technique; effective multiplication factor

1 INTRODUCTION

A fundamental issue in nuclear research is to estimate the neutron population in nuclear systems. It is a great challenge from the physical and mathematical point of view and a crucial problem in the study and analysis of the nuclear reactors ^{4), (14}. In this context, although the theory of neutron diffusion for global calculations in reactor physics has limited validity, it is widely used because it provides satisfactory results for many applications. One of them is the analysis of the criticality of nuclear reactors ^{1), (3), (7), (15}, where the main issue is to establish the ratio between the numbers of neutrons generated in successive fission reactions. In this way, the criticality is evaluated by the dominant eigenvalue of the stationary Neutron Diffusion Equation and its corresponding eigenvector is associated with the neutron scalar flux. Accurate calculations of criticality and neutron flux are essential for the licensing and secure operations of nuclear power plants.

Due to the need to precisely determining the neutron flux, many methods were developed with this purpose. Such methods take into consideration the geometry and composition of the reactor core. Among those methods, we can mention the calculation of the neutron flux from the multigroup diffusion equation through the finite difference method (FDM) ^{16), (17} and finite volume method (FVM) ¹. According to Ref. ¹², these methods are simple and of easy implementation, however, they require refined meshes to obtain the desired precision. For this reason, the efficiency of these codes may be limited for application in large-scale coupled problems ¹⁹. In order to avoid the great computational effort of those methods, other methodologies have been proposed, such as nodal methods.

The main idea of the nodal methods is the integration of the diffusion equation over a region (node) of the domain, defining average fluxes in each node and reducing the complexity of the model. However, unknown variables are introduced in this process, the leakage terms, and so, a fundamental problem in nodal methods is to accurately determine the spatial leakage coefficients. Nodal methods have been widely used in reactor physics analysis for many years. According to Refs. ^{5), (8} the first simulator based on this method is the code FLARE developed in 1964. Still, according to the literature ¹⁹ this approach is more efficient and accurate in thick meshes if compared to FVM and FDM, so that the discrete mesh number, computer storage, and computational cost can be greatly reduced for a desired precision. Because of the advantages of transverse-integrated nodal methods, some of them have been applied to the more diverse models in reactor physics such as multi-region diffusion problems in two and three dimensions. In this context, a first analytical approach was proposed, the so-called analytical nodal method (ANM) ^{11), (13}. In the ANM, the equations are derived for each spatial variable separately, such that a set of ordinary differential equations (ODE’s) is obtained. The ODE’s are analytically solved and the only approximation is the transverse leakage term in each direction ¹³. The application of the boundary conditions imposes the relation between the directions and the coupling of the system. However, for problems with more than two groups, the calculation of the nodal coupling matrices turned out to be a complex task. More recently ⁶ presented a simplified formulation of ANM, which is more easily applicable to an arbitrary number of energy groups. In this case, linear transformation techniques are used to allow the diagonalization of a relevant matrix, when possible. Then the ANM procedure is applied to a decoupled problem. After that, the bissection method is used to find the effective multiplication factor K _eff as the root of the characteristic equation derived from a global matrix. Results for the scalar fluxes, however, are not provided.

In this context, the solution of the eigenvalue problem is a permanent topic of investigation, as in ², where block iterative methods were introduced: the block inverse-free preconditioned Arnoldi method and the modified block Newton method. All of these iterative solvers are initialized using a block multilevel technique. The authors also proposed a hybrid multilevel method, based on the combination of the preconditioned block inverse-free Krylov method and the generalized modified block Newton method. In solving the same class of problems, ³, have used a reduced order method, in particular reduced basis finite element method. According the authors, in such approach, efficiency and accuracy are ensured by a full decoupling of the finite element scheme and the reduced order model (the so-called offline-online procedure).

In 2018, ¹⁸ proposed a methodology to solve the multigroup multilayer stationary neutron diffusion equation in two-dimensional Cartesian geometry through the source iteration method. The Fourier Transform was used to develop an analytical solution in each step of the source iteration, instead of using traditional numerical methods. Once that problem was solved the inverse Fourier Transform was required to reconstruct the neutron flux in the original variables. Nevertheless, the usual update requested in the source term for every step leads to additional terms to be inverted and the process may become too laborious. To overcome this problem, the neutron flux was approximated by polynomial interpolation. Still, in order to reduce computational time, the original domain was subdivided into small fictitious regions. In this way, the interpolating polynomials in each region could be chosen to be of a low order, reducing the dimensions of the matrices involved and, consequently, the computational time. This idea, however, was not effective for dealing with two-dimensional problems in heterogeneous media defined by asymmetric boundaries.

In this paper, we present a procedure that is able to determine not only the effective multiplication factor but also to estimate the neutron average scalar flux. The technique is applied in multigroup diffusion problems in heterogeneous media that may be defined by asymmetric boundaries and involve upscattering and downscattering. We derive a nodal formulation from the original model in which the solution is explicitly written in terms of the spatial variables. Firstly, we decouple the system of differential equations assuming as known the source terms associated with fission and scattering, in the right-hand side (source iteration). It allows us to solve the diffusion equation analytically for each group independently, starting from group one. In this way, the unknowns of the ODE’s solutions can be determined for each group separately, through smaller order linear systems. These systems are established when we impose that these ODE’s solutions must satisfy the boundary conditions and the flux and current density continuity conditions at the node interfaces. Once the neutron average scalar flux is obtained for each group, we update the value of K _eff and we continue the iterative process evaluating the flux and K _eff up to a certain stop criterion. The procedure may be considered simple and has shown to be successfully applied to a series of problems. To implement this iterative procedure two approximations are required. Firstly, the unknown leakage terms on the contours, as usual in nodal techniques, are here assumed to be constant. Although higher-order approximations have been proposed along with nodal methods, here, differently, we choose to explore alternative representations to the fluxes updates on the source term. In particular, we investigate constant and linear approximations.

In this way, this paper is organized such that in Section 2 we present the diffusion eigenvalue problem of interest, as well as the development of the nodal and iterative scheme for the solution to the proposed problem. In Section 3, we solve four test problems, for the cases of two and four energy groups. Finally, in Section 4, we discuss concluding aspects and ongoing projects.

2 MATHEMATICAL FORMULATION

We begin with the two-dimensional multigroup steady-state neutron diffusion equation in a rectangular domain Ω, written as ⁴

$$\begin{gathered} \frac{\partial }{\partial x}\left(-D_g(x,y)\frac{\partial }{\partial x}{\varphi }_g(x,y)\right)+\frac{\partial }{\partial y}\left(-D_g(x,y)\frac{\partial }{\partial y}{\varphi }_g(x,y)\right)+{\Sigma }_{Rg}(x,y){\varphi }_g(x,y)= \\ \frac{{\chi }_g}{K_{eff}}\sum _{g\text{'}=1}^G{\nu }_{g\text{'}}{\Sigma }_{fg\text{'}}(x,y){\varphi }_{g\text{'}}(x,y)+\sum _{g\text{'}=1 \\ g\text{'}\ne g}^G{\Sigma }_{sg\text{'}g}(x,y){\varphi }_{g\text{'}}(x,y), \end{gathered}$$ (2.1)

where g = 1, . . . , G stands for neutron energy groups. For each energy group g: ϕ _g (x, y) is the neutron scalar flux; D _g (x, y) is the diffusion coefficient; Σ_Rg (x, y) is the removal cross section; χ _g is the integrated fission spectrum; ν _g is the average number of neutrons emitted by fission; Σ_fg (x, y) is the fission cross section; Σ_sg'g (x, y) is the scattering cross section from energy group g' to g and K _eff is the effective multiplication factor. It is worth mentioning that, the cross section represents the probability of the occurrence of a specific interaction between neutrons and target nuclei. For example, the fission cross section is the probability that a neutron and a nucleus interact to form a new nucleus that undergoes fission ¹⁴.

We complete the formulation of our problem by imposing boundary conditions written in a general form as

${\alpha }_g{\varphi }_g+{\beta }_g\frac{\partial {\varphi }_g}{\partial n}=0,$ (2.2)

where n may be x or y depending on the boundary under consideration, α _g and β _g are constants and $\left|{\alpha }_g\right|+\left|{\beta }_g\right|>0$. In addition, when solving multi-region problems we also consider continuity conditions at the interfaces.

A well known approach to deal with this type of problem is the Source Iteration method (SI), which is an iterative procedure ¹⁰. In each step, we suppose the right-hand side of Eq. (2.1) is known and we then need to find a solution to the resulting equation. In what follows, we detail our nodal approach for solving the referred resulting equation.

2.1 A nodal approach

We divide the domain Ω in rectangular nodes (i, j), where $x\in [x_{i-1},x_i],y\in [y_{j-1},y_j]$ and the nuclear parameters are assumed to be constant over each node. Then, we rewrite Eq. (2.1), for the group g in the node (i, j), as

$$\begin{gathered} \frac{\partial }{\partial x}\left(-D_g^{(i,j)}\frac{\partial }{\partial x}{\varphi }_g^{(i,j)}(x,y)\right)+\frac{\partial }{\partial y}\left(-D_g^{(i,j)}\frac{\partial }{\partial y}{\varphi }_g^{(i,j)}(x,y)\right)+{\Sigma }_{Rg}^{(i,j)}{\varphi }_g^{(i,j)}(x,y)= \\ \frac{{\chi }_g}{K_{eff}}\sum _{g\text{'}=1}^G{\nu }_{g\text{'}}{\Sigma }_{fg\text{'}}^{(i,j)}{\varphi }_{g\text{'}}^{(i,j)}(x,y)+\sum _{g\text{'}=1 \\ g\text{'}\ne g}^G{\Sigma }_{sg\text{'}g}^{(i,j)}{\varphi }_{g\text{'}}^{(i,j)}(x,y). \end{gathered}$$ (2.3)

As we mentioned, we propose to use the SI method in Eq. (2.3), such that K _eff and ${\varphi }_g^{(i,j)}(x,y)$, on the right side of the equation, are assumed to be known, in each iteration. Besides that, we express the unknown fluxes on the right-hand side, at each node (i, j), as linear polynomials

${\varphi }_{g\text{'}}^{(i,j)}(x,y)=a_{0g\text{'}}^{(i,j)}+a_{1g\text{'}}^{(i,j)}x+a_{2g\text{'}}^{(i,j)}y.$ (2.4)

In the first iteration, initial guesses are given to the constants $a_{0g\text{'}}^{(i,j)},a_{1g\text{'}}^{(i,j)},a_{2g\text{'}}^{(i,j)}$ and K _eff . From the second iteration, the constants $a_{0g\text{'}}^{(i,j)},a_{1g\text{'}}^{(i,j)}$ and $a_{2g\text{'}}^{(i,j)}$ will be determined from known values of the solution of the problem, as we will discuss later on in this text. We substitute Eq. (2.4) on the right-hand side of the Eq. (2.3) to obtain

$$\begin{gathered} \frac{\partial }{\partial x}\left(-D_g^{(i,j)}\frac{\partial }{\partial x}{\varphi }_g^{(i,j)}(x,y)\right)+\frac{\partial }{\partial y}\left(-D_g^{(i,j)}\frac{\partial }{\partial y}{\varphi }_g^{(i,j)}(x,y)\right)+{\Sigma }_{Rg}^{(i,j)}{\varphi }_g^{(i,j)}(x,y)= \\ {c}_{g0}^{(i,j)}+c_{g1}^{(i,j)}x+c_{g2}^{(i,j)}y, \end{gathered}$$ (2.5)

where

$$\begin{gathered} c_{g0}^{(i,j)}=\frac{{\chi }_g}{K_{eff}}\sum _{g\text{'}=1}^G{\nu }_{g\text{'}}{\Sigma }_{fg\text{'}}^{(i,j)}{a}_{0g\text{'}}^{(i,j)}+\sum _{g\text{'}=1 \\ g\text{'}\ne g}^G{\Sigma }_{sg\text{'}g}^{(i,j)}{a}_{0g\text{'}}^{(i,j)}, \end{gathered}$$ (2.6)

$$\begin{gathered} c_{g1}^{(i,j)}=\frac{{\chi }_g}{K_{eff}}\sum _{g\text{'}=1}^G{\nu }_{g\text{'}}{\Sigma }_{fg\text{'}}^{(i,j)}{a}_{1g\text{'}}^{(i,j)}+\sum _{g\text{'}=1 \\ g\text{'}\ne g}^G{\Sigma }_{sg\text{'}g}^{(i,j)}{a}_{1g\text{'}}^{(i,j)} \end{gathered}$$ (2.7)

and

$$\begin{gathered} c_{g2}^{(i,j)}=\frac{{\chi }_g}{K_{eff}}\sum _{g\text{'}=1}^G{\nu }_{g\text{'}}{\Sigma }_{fg\text{'}}^{(i,j)}{a}_{2g\text{'}}^{(i,j)}+\sum _{g\text{'}=1 \\ g\text{'}\ne g}^G{\Sigma }_{sg\text{'}g}^{(i,j)}{a}_{2g\text{'}}^{(i,j)}. \end{gathered}$$ (2.8)

Finally, we integrate Eq. (2.5) over all y in [y _j−1 , y _j ] and divide by $\Delta y^{(i,j)}=y_j-y_{j-1}$,

$$\begin{gathered} -D_g^{(i,j)}\frac{d}{d{x}^2}{\varphi }_{gy}^{(i,j)}\left(x\right)+\frac{1}{\Delta y^{(i,j)}}\left(-D_g^{(i,j)}\frac{\partial }{\partial y}{\varphi }_g^{(i,j)}(x,y){|}_{y_{j-1}}^{y_j}\right)+{\Sigma }_{Rg}^{(i,j)}{\varphi }_{gy}^{(i,j)}\left(x\right)= \\ {c}_{g0}^{(i,j)}+c_{g1}^{(i,j)}x+\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2\Delta y^{(i,j)}}. \end{gathered}$$ (2.9)

In this previous derivation, we defined the average scalar flux in group g along the y-direction as

${\varphi }_{gy}^{(i,j)}\left(x\right)=\frac{1}{\Delta y^{(i,j)}}{\int }_{y_{j-1}}^{y_j}{\varphi }_g^{(i,j)}(x,y)dy.$ (2.10)

At this point, it is important to note that, due to the integration procedure, additional unknowns terms appear in Eq. (2.9). In particular, the second term on the left-hand side of Eq. (2.9) is, in general, unknown for an arbitrary node (i, j). That is a classical issue of nodal schemes. Here, we write the x-directed transverse leakage term as

$J_{gx}^{(i,j)}\left(y\right)\approx -D_g^{(i,j)}\frac{\partial }{\partial y}{\varphi }_g^{(i,j)}(x,y),$ (2.11)

and we assume it is constant, $J_{gx}^{(i,j)}\left(y_j\right)$ and $J_{gx}^{(i,j)}\left(y_{j-1}\right)$, in the contours y _j and y _j-1 , respectively, of the node (i, j). We proceed and substitute Eq. (2.11) into Eq. (2.9) to obtain

$$\begin{gathered} -D_g^{(i,j)}\frac{d}{d{x}^2}{\varphi }_{gy}^{(i,j)}\left(x\right)+\frac{1}{\Delta y^{(i,j)}}\left(J_{gx}^{(i,j)}\left(y_j\right)-J_{gx}^{(i,j)}\left(y_{j-1}\right)\right)+{\Sigma }_{Rg}^{(i,j)}{\varphi }_{gy}^{(i,j)}\left(x\right)= \\ {c}_{g0}^{(i,j)}+c_{g1}^{(i,j)}x+\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2\Delta y^{(i,j)}}. \end{gathered}$$ (2.12)

We multiply Eq. (2.12) by $1/D_g^{(i,j)}$ to obtain

$$\begin{gathered} \frac{d}{d{x}^2}{\varphi }_{gy}^{(i,j)}\left(x\right)-\left(\frac{J_{gx}^{(i,j)}\left(y_j\right)-J_{gx}^{(i,j)}\left(y_{j-1}\right)}{D_g^{(i,j)}\Delta y^{(i,j)}}\right)-{\varpi }_g^{(i,j)}{\varphi }_{gy}^{(i,j)}\left(x\right)= \\ -\frac{c_{g0}^{(i,j)}}{D_g^{(i,j)}}-\frac{c_{g1}^{(i,j)}}{D_g^{(i,j)}}x-\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2D_g^{(i,j)}\Delta y^{(i,j)}}, \end{gathered}$$ (2.13)

where ${\varpi }_g^{(i,j)}=\frac{{\Sigma }_{Rg}^{(i,j)}}{D_g^{(i,j)}}$, to rewrite Eq. (2.13) as

$$\begin{gathered} \frac{d^2}{d{x}^2}{\varphi }_{gy}^{(i,j)}\left(x\right)-{\varpi }_g^{(i,j)}{\varphi }_{gy}^{(i,j)}\left(x\right)=\left(\frac{J_{gx}^{(i,j)}\left(y_j\right)-J_{gx}^{(i,j)}\left(y_{j-1}\right)}{D_g^{(i,j)}\Delta y^{(i,j)}}\right) \\ -\frac{c_{g0}^{(i,j)}}{D_g^{(i,j)}}-\frac{c_{g1}^{(i,j)}}{D_g^{(i,j)}}x-\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2D_g^{(i,j)}\Delta y^{(i,j)}}. \end{gathered}$$ (2.14)

Eq. (2.14) is an ordinary differential equation (ODE), that has an analytical solution of the form

$$\begin{gathered} {\varphi }_{gy}^{(i,j)}\left(x\right)=A_x^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x-x_{i-1})}+B_x^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x)}-\frac{J_{gx}^{(i,j)}\left(y_j\right)-J_{gx}^{(i,j)}\left(y_{j-1}\right)}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta y^{(i,j)}} \\ +\frac{c_{g0}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}+\frac{c_{g1}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}x+\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta y^{(i,j)}}. \end{gathered}$$ (2.15)

Similarly, if we integrate Eq. (2.5) for all x in [x _i−1 , x _i ], divide by $\Delta x^{(i,j)}=x_i-x_{i-1}$ and approximate by a constant the y-directed transverse leakage term in the contours of the node (i, j) as $J_{gy}^{(i,j)}\left(x_i\right)$ and $J_{gy}^{(i,j)}\left(x_{i-1}\right)$, we obtain an ODE in y for the average scalar flux in the x-direction of group g

$$\begin{gathered} \frac{d^2}{d{y}^2}{\varphi }_{gx}^{(i,j)}\left(y\right)-{\varpi }_g^{(i,j)}{\varphi }_{gx}^{(i,j)}\left(y\right)=\left(\frac{J_{gy}^{(i,j)}\left(x_i\right)-J_{gy}^{(i,j)}\left(x_{i-1}\right)}{D_g^{(i,j)}\Delta x^{(i,j)}}\right) \\ -\frac{c_{g1}^{(i,j)}(x_i^2-x_{i-1}^2)}{2D_g^{(i,j)}\Delta x^{(i,j)}}-\frac{c_{g2}^{(i,j)}}{D_g^{(i,j)}}y. \end{gathered}$$ (2.16)

Again, we write the analytical solution of Eq. (2.16) as

$$\begin{gathered} {\varphi }_{gx}^{(i,j)}\left(y\right)=A_y^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y-y_{j-1})}+B_y^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y)}-\frac{J_{gy}^{(i,j)}\left(x_i\right)-J_{gy}^{(i,j)}\left(x_{i-1}\right)}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta x^{(i,j)}} \\ +\frac{c_{g0}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}+\frac{c_{g1}^{(i,j)}(x_i^2-x_{i-1}^2)}{2{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta x^{(i,j)}}+\frac{c_{g2}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}y. \end{gathered}$$ (2.17)

The solutions obtained in Eqs. (2.15) and (2.17) allow us to propose updated approximations for the transversal leakage terms. Thus, when we derive the Eq. (2.15) with respect to x

$$\begin{gathered} \frac{d{\varphi }_{gy}^{(i,j)}\left(x\right)}{dx}=-A_x^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x-x_{i-1})}+B_x^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x)} \\ +\frac{c_{g1}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}, \end{gathered}$$ (2.18)

such that we can now express the y-direction transverse leakage term as

$$\begin{gathered} J_{gy}^{(i,j)}\left(x\right)\approx -D_g^{(i,j)}(-A_x^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x-x_{i-1})} \\ +B_x^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x)}+\frac{c_{g1}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}). \end{gathered}$$ (2.19)

Analogously, when we derive the Eq. (2.17) with respect to y

$$\begin{gathered} \frac{d{\varphi }_{gx}^{(i,j)}\left(y\right)}{dy}=-A_y^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y-y_{j-1})} \\ +B_y^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y)}+\frac{c_{g2}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}, \end{gathered}$$ (2.20)

such that we can now express the x-direction transverse leakage term as

$$\begin{gathered} J_{gx}^{(i,j)}\left(y\right)\approx -D_g^{(i,j)}(-A_y^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y-y_{j-1})} \\ +B_y^{(i,j)}\sqrt{{\varpi }_g^{(i,j)}}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y)}+\frac{c_{g2}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}). \end{gathered}$$ (2.21)

Finally, we substitute Eq. (2.21), evaluated in y _j and y _j−1 , in Eq. (2.15), to rewrite the average scalar flux of group g in the x-direction, for each node (i, j), as

$$\begin{gathered} {\varphi }_{gy}^{(i,j)}\left(x\right)=A_x^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x-x_{i-1})}+B_x^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x)} \\ -A_y^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y_{j-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta y^{(i,j)}}-B_y^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y_{j-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta y^{(i,j)}} \\ +\frac{c_{g0}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}+\frac{c_{g1}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}x+\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta y^{(i,j)}}. \end{gathered}$$ (2.22)

Similarly, if we substitute Eq. (2.19), evaluated in x _i and x _i−1 , in Eq. (2.17), we rewrite the average scalar flux in the y-direction, for each node (i, j), as

$$\begin{gathered} {\varphi }_{gx}^{(i,j)}\left(y\right)=A_y^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y-y_{j-1})}+B_y^{(i,j)}{e}^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y)} \\ -A_x^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x_{i-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta x^{(i,j)}}-B_x^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x_{i-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta x^{(i,j)}} \\ +\frac{c_{g0}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}+\frac{c_{g1}^{(i,j)}(x_i^2-x_{i-1}^2)}{2{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta x^{(i,j)}}+\frac{c_{g2}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}y. \end{gathered}$$ (2.23)

Here the constants A _x , A _y , B _x and B _y are still unknown. To establish the general solutions, we impose that the Eqs. (2.22) and (2.23) must satisfy the boundary conditions along with flux and current density continuity conditions at the interfaces of the nodes. We generate a 4R linear system by energy group, where R is the total number of nodes. We solve a linear system for each group to obtain A _x , A _y , B _x and B _y such that the general solutions, Eqs. (2.22) and (2.23), are completely established.

In order to obtain the average scalar flux, ${\overline{\varphi }}_g^{(i,j)}$, in each node, we integrate either Eq. (2.22) over all x in [x _i−1 , x _i ] or Eq. (2.23) over all y in [y _j−1 , y _j ], then we divide the resulting equations respectively by ∆x ^{(i, j)} or ∆y ^{(i, j)} , to find the same final expression

$$\begin{gathered} {\overline{\varphi }}_g^{(i,j)}=-A_x^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x_{i-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta x^{(i,j)}}-B_x^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(x_i-x_{i-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta x^{(i,j)}} \\ -A_y^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y_{j-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta y^{(i,j)}}-B_y^{(i,j)}\frac{e^{-\sqrt{{\varpi }_g^{(i,j)}}(y_j-y_{j-1})}-1}{\sqrt{{\varpi }_g^{(i,j)}}\Delta y^{(i,j)}} \\ +\frac{c_{g0}^{(i,j)}}{{\varpi }_g^{(i,j)}{D}_g^{(i,j)}}+\frac{c_{g1}^{(i,j)}(x_i^2-x_{i-1}^2)}{2{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta x^{(i,j)}}+\frac{c_{g2}^{(i,j)}(y_j^2-y_{j-1}^2)}{2{\varpi }_g^{(i,j)}{D}_g^{(i,j)}\Delta y^{(i,j)}}. \end{gathered}$$ (2.24)

Once we found the averaged scalar flux in each node, we are ready to update the coefficients of Eq. (2.4), as we mentioned previously. Such coefficients are obtained through an interpolation based on known values of the averaged flux at the node (i, j) and at the neighbords nodes $(i+1,j)$ and $(i,j+1)$

$a_{1g}^{(i,j)}=\frac{2\left({\overline{\varphi }}_g^{(i+1,j)}-{\overline{\varphi }}_g^{(i,j)}\right)}{\Delta x^{(i,j)}+\Delta x^{(i+1,j)}},$ (2.25)

$a_{2g}^{(i,j)}=\frac{2\left({\overline{\varphi }}_g^{(i,j+1)}-{\overline{\varphi }}_g^{(i,j)}\right)}{\Delta y^{(i,j)}+\Delta y^{(i,j+1)}}$ (2.26)

and

$a_{0g}^{(i,j)}={\overline{\varphi }}_g^{(i,j)}-a_{1g}^{(i,j)}{x}_c^{(i,j)}-a_{2g}^{(i,j)}{y}_c^{(i,j)},$ (2.27)

where $x_c^{(i,j)}$ means x-coordinate in the center of the node (i, j) and $y_c^{(i,j)}$ means y-coordinate in the center of the node (i, j). The definitions given in Eqs. (2.25) and (2.26) do not apply to boundary nodes, since they depend on averaged fluxes of the neighbours nodes. In these cases, we use Eq. (2.22) or Eq. (2.23). To be clear, the coefficients for the right boundary nodes are obtained by interpolating the average fluxes ${\overline{\varphi }}_g^{(i,j)}$ and ${\overline{\varphi }}_g^{(i,j+1)}$ and the fluxes evaluated on the right contour, ${\varphi }_{gy}^{(i,j)}\left(x_i\right)$. So, we get the coefficients

$a_{1g}^{(i,j)}=\frac{2\left({\varphi }_{gy}^{(i,j)}\left(x_i\right)-{\overline{\varphi }}_g^{(i,j)}\right)}{\Delta x^{(i,j)}}.$ (2.28)

The coefficients $a_{0g}^{(i,j)}$ and $a_{2g}^{(i,j)}$ for these nodes are given in the same way, as in Eqs. (2.27) and (2.26), respectively.

Similarly, the coefficients for the top boundary nodes are obtained by interpolating the average fluxes ${\overline{\varphi }}_g^{(i,j)}$ and ${\overline{\varphi }}_g^{(i+1,j)}$ and the fluxes of the top contour ${\varphi }_{gx}^{(i,j)}\left(y_j\right)$. So, we get the coefficients

$a_{2g}^{(i,j)}=\frac{2\left({\varphi }_{gx}^{(i,j)}\left(y_j\right)-{\overline{\varphi }}_g^{(i,j)}\right)}{\Delta y^{(i,j)}}.$ (2.29)

The coefficients $a_{0g}^{(i,j)}$ and $a_{1g}^{(i,j)}$ for these nodes are given in the same way, as in Eqs. (2.27) and (2.25), respectively. In the case of the boundary nodes with both contours (rigth and top), the coefficients are those of Eqs. (2.27), (2.28) and (2.29). Note that if the boundary conditions are of zero flux, we do not need to evaluate the fluxes in the contours by Eqs. (2.22) and Eq. (2.23), since ${\varphi }_{gy}^{(i,j)}\left(x_i\right)=0$ and ${\varphi }_{gx}^{(i,j)}\left(y_j\right)=0$.

In order to analyze the influence of the order approximation on the source term, we also tested, in this work, constant functions, in addition to the linear expansions. In the constant case, the fluxes of each node are updated with the average fluxes, Eq. (2.24), obtained in the previous iteration, therefore, we do not need to interpolate. Once we updated the fluxes on the right-hand side, to complete the SI procedure, we update the K _eff , as we describe in the text below.

2.2 An updated value to K _eff

We now introduce the $GR\times 1$ average scalar flux vector Φ, where G is the total number of energy groups and R is the total number of nodes, with components given by Eq. (2.24). In addition, we define the block diagonal matrix F, $GR\times GR$ order, where each one of the R blocks associated with nodes (i, j) is defined by a $G\times G$ square matrix F ^{(i, j)} in the form

${\mathbf{F}}^{(i,j)}=\left[\begin{array}{ccc}{\chi }_1{\nu }_1{\Sigma }_{f1}^{(i,j)}& \cdots & {\chi }_1{\nu }_G{\Sigma }_{fG}^{(i,j)}\\ ⋮& \ddots & ⋮\\ {\chi }_G{\nu }_1{\Sigma }_{f1}^{(i,j)}& \cdots & {\chi }_G{\nu }_G{\Sigma }_{fG}^{(i,j)}\end{array}\right].$ (2.30)

We follow ⁴ and define an iterative scheme for determining K _eff such that in each iteration ”n”

(2.31)

where the symbol ⟨u , v⟩ means the inner product. In each iteration, after obtaining the new values for the average scalar fluxes and K _eff , we update the polynomials in Eq. (2.4) according to Eqs. (2.25), (2.26), (2.27), (2.28) and (2.29). The iterative process follows up to a certain stop criterion as

$\left|\frac{K_{eff}^{\left[n\right]}-K_{eff}^{[n-1]}|}{K_{eff}^{\left[n\right]}}\right|<{ϵ}_1\text{and}\frac{{\displaystyle {||{\overline{\mathrm{\Phi }}}^{\left[n\right]}-{\overline{\mathrm{\Phi }}}^{[n-1]}||}_2}}{{\displaystyle {||{\overline{\mathrm{\Phi }}}^{\left[n\right]}||}_2}}<{ϵ}_2.$ (2.32)

3 RESULTS AND DISCUSSION

We solved four test problems, to which results are available in the literature ^{6), (9), (13), (18}, for the sake of comparisons. In our calculations, the fluxes in the source term are updated by constant and linear approximations, in order to analyze the influence of the approximation orders. Thus, in the tables below we refer to Constant Source and Linear Source. All numerical results were generated by programs implemented in Fortran95 and ran on a computer with an Intel Core i5 − 8250U , 1.60GHz processor and 8GB of RAM. In test problems, we present the graphs of the relative errors of K _{ef f} as a function of the ln(R), where R is the total number of nodes. The stop criterions adopted for the four test problems are ε ₁ = ε ₂ = 10⁻⁵.

3.1 Problem 1

This first problem we consider is a two groups heterogeneous medium problem proposed in Ref. ¹⁸. The geometry is defined as a rectangular domain, [0, 36] × [0, 12] (cm), divided into three 12 cm × 12 cm regions (see Fig. 1). The first and third regions are composed by the same reflective material and the central region is a fissile material. Zero flux is assumed as boundary conditions for all boundaries.

Figure 1:
Geometry of Problem 1.

Due to the symmetry conditions of the problem, we redefine the domain of the problem (dotted line) as [0, 18] ×[0, 6](cm), in order to reduce the computational effort. In this new configuration, zero flux boundary conditions are kept on the left and bottom boundaries while reflection is imposed on the two other boundaries (right and top) for each group g. With respect to the parameters in Eq. (2.1), we follow ¹⁸ and the cross sections for each group are given in Table 1.

Table 1:
Nuclear parameters of Problem 1.

Material	Group	Dg (cm)	ΣRg (cm −1)	νΣ fg (cm −1 )	Σs1g (cm −1)	χ g
1	1	1.6497	0.02309	0.005008	0.00000	1.0
1	2	0.4754	0.07886	0.097130	0.01423	0.0
2	1	0.6702	0.09013	0.000000	0.00000	1.0
2	2	0.1509	0.03277	0.000000	0.09084	0.0

In Table 2 we present results for the effective multiplication factor for several mesh configurations. The behavior of the numerical results in Ref. ¹⁸ is similar to our results when using linear approximations. However, that methodology can not be extended, for instance, to deal with heterogeneous media problems where the domain has to be divided into rectangular regions (subdividing the x- and y-axis) instead of layers. In fact, the results listed in Table 2, presented in Ref. ¹⁸, were generated just considering subdivision (layers) for the interval of definition in the variable x, in order to allow the computational evaluation of the solutions. In Table 2, we also see a better agreement of the linear approximation compared to the constant, furthermore, the computational cost does not increase significantly.

Table 2:
Comparison of effective multiplication factor Problem 1.

∆x (cm)	Ref. 18	Constant Source∗	t (s)	Linear Source∗	t (s)
1	0.13808	0.13561	0.010	0.13795	0.013
0.5	0.13805	0.13742	0.147	0.13802	0.161
0.25	0.13804	0.13788	5.454	0.13803	5.513

• *∆y = ∆x

In Table 3, we present results for the fast neutrons fluxes (g = 1). The results presented in the columns with superscript ^∗ refer to Ref. ¹⁸ and our results are generated by the proposal that approximates the fluxes of the source term by linear polynomials (Linear Source). The results presented in Table 3 show a good agreement between the two methodologies, where the relative error is less than 0.135%. The errors of thermal fluxes (g = 2) are similar to fast fluxes.

Table 3:
Comparison of fast neutrons fluxes Problem 1.

x	ϕ1(x, 2)	ϕ1∗(x, 2)	Error(%)	ϕ1(x, 6)	ϕ1∗(x, 6)	Error(%)
2	0.003647	0.003646	0.027	0.007294	0.007291	0.041
6	0.026370	0.026353	0.064	0.052740	0.052707	0.063
10	0.160665	0.160498	0.104	0.321330	0.320996	0.104
14	0.476699	0.476816	0.024	0.953402	0.953631	0.024
18	0.500000	0.500000	0.000	1.000000	1.000000	0.000

3.2 Problem 2

The second problem we consider is the TWIGL benchmark described by ¹³. It is a problem with two energy groups and three types of material, which geometry (in cm) is presented in Figure 2a.

Figure 2:
Geometry of Problem 2.

In Table 4, we list the nuclear parameters for each type of material and energy group according to Ref. ¹³. Furthermore, we emphasize that, in the case of the stationary problem, materials 1 and 2 are identical.

Table 4:
Nuclear parameters of Problem 2.

Material	Group	Dg (cm)	ΣRg (cm −1)	νΣ fg (cm −1 )	Σs1g (cm −1)	χ g
1 and 2	1	1.4	0.010	0.007	0.00	1.0
1 and 2	2	0.4	0.150	0.200	0.01	0.0
3	1	1.3	0.008	0.003	0.00	1.0
3	2	0.5	0.050	0.060	0.01	0.0

In Figure 2b, we redefine the material regions of the original problem in nine regions: Region I: [0, 24] × [0, 24]; Region II: [24, 56] × [0, 24]; Region III: [56, 80] × [0, 24]; Region IV: [0, 24] × [24, 56]; Region V: [24, 56] ×[24, 56]; Region VI: [56, 80] ×[24, 56]; Region VII: [0, 24] ×[56, 80]; Region VIII: [24, 56] × [56, 80] and Region IX: [56, 80] × [56, 80]. To evaluate the problem we define 4 mesh configurations on these subintervals: Mesh 1: 12 × 12 in Region I, 16 × 12 in Region II, 12 × 12 in Region III, 12 × 16 in Region IV, 16 × 16 in Region V, 12 × 16 in Region VI, 12 × 12 in Region VII, 16 × 12 in Region VIII and 12 × 12 in Region IX, and other meshes are equally spaced; Mesh 2: 8 × 8; Mesh 3: 4 × 4 and Mesh 4: 2 × 2. To complete the formulation of the problem, zero flux is assumed at the top and right boundaries while reflection is considered on the left and bottom boundaries.

The results obtained with our code are compared with Smith’s results generated by the code QUANDRY ¹³. The code QUANDRY implements a nodal methodology where it is assumed quadratic transverse leakage approximations, while we use constant leakage approximations. In Table 5 we present numerical results for the effective multiplication factor obtained for several mesh configurations. We can note that the linear source approximations provide better agreement with the reference results. Results for more refined meshes are not available in Ref. ¹³ and it seems a complicated task the extension of that methodology for problems with more than two groups. We also observe that our results are generated with low computational cost, even in more refined meshes.

Table 5:
Comparison of effective multiplication factor Problem 2.

Mesh	Smith (1979) Ref. 13	Constant Source	t (s)	Linear Source	t (s)
1	0.91323	0.90070	0.007	0.91111	0.009
2	0.91321	0.90859	0.023	0.91293	0.027
3	-	0.91205	0.269	0.91318	0.269
4	-	0.91292	5.445	0.91321	5.552

In Figure 3, we present the graph of the relative error of K _eff produced by the constant and linear approximations of the source, if we consider as a reference value K _eff = 0.91321. In this figure, we can see that linear approximations generate more accurate results than constant approximations, especially in thicker meshes. Moreover, we observed that as we increase the order of approximation of the source term it does not result in a significant increase in the computational cost.

Figure 3:
Illustration of the relative error of Problem 2.

3.3 Problem 3

According to Ref. ⁶ this problem is an ENE-6103 academic benchmark. It is a two groups problem where the domain is a 350 cm × 350 cm square region composed of two types of materials, as we see in Figure 4a.

Figure 4:
Geometry of Problem 3.

Zero flux is assumed at the top and right boundaries while reflection is considered on the left and bottom boundaries. In Table 6 we present the relevant parameters for defining the problem ⁶.

Table 6:
Nuclear parameters of Problem 3.

Material	Group	Dg (cm)	ΣRg (cm −1 )	νΣ fg (cm −1 )	Σs1g (cm −1 )	χ g
1	1	1.2640	8.154 × 10−3	0.0000	0.0000	1.0
1	2	0.9328	4.100 × 10−3	4.562 × 10−3	7.368 × 10−3	0.0
2	1	3100	1.018 × 10−2	0.0000	0.0000	1.0
2	2	0.8695	2.117 × 10−4	0.0000	1.018 × 10−2	0.0

In Figure 4b, we redefine the material regions of the original problem in four regions: Region I:[0, 310] × [0, 310]; Region II: [310, 350] × [0, 310]; Region III: [0, 310] × [310, 350]; Region IV: [310, 350] × [310, 350]. To evaluate the effective multiplication factor we define 3 mesh configurations on these subintervals: Mesh 1: 310 × 310 in Region I, 40 × 310 in Region II, 310 × 40 in Region III and 40 × 40 in Region IV; Mesh 2: 62 × 62 in Region I, 40 × 62 in Region II, 62 × 40 in Region III and 40 × 40 in Region IV and Mesh 3: 31 × 31 in Region I, 20 × 31 in Region II, 31 × 20 in Region III and 20 × 20 in Region IV.

The results obtained with our code are compared with Hébert’s results and results generated by the code QUANDRY both available in Ref. ⁶, (see Table 7). Hérbert developed a nodal methodology where it is assumed either flat or quadratic transverse leakage approximations, while we use only constant leakage approximations. In Table 7 we can see results provided by the QUANDRY code that, for the mesh configurations presented, seem to have up to three correct digits. Those results agree with the ones provided by Herbert using quadratic approximations. The same agreement is found in our results generated by the linear source approximation. We note that in the Ref. ⁶ a diagonalization process is required to decouple the problem. Besides, in our methodology, we can determine not only the effective multiplication factor but also to estimate the neutron average scalar flux for each node simultaneously. In Ref. ⁶ just the K _eff was evaluated.

Table 7:
Comparison of effective multiplication factor Problem 3.

Mesh	Flat	Hérbert (2008) Quadratic	6 QUANDRY	Constant Source	t (s)	Linear Source	t (s)
1	0.991733	0.990257	0.990410	0.9516	0.001	0.9927	0.001
2	0.990186	0.990106	0.990106	0.9759	0.013	0.9907	0.020
3	-	-	-	0.9854	0.101	0.9903	0.132

In Figure 5, we can also see the relative error of the results considering as reference K _eff = 0.990106 generated by QUANDRY code. These errors show that linear approximations for the fluxes of the source term generate results much more accuracy than the constant approximations, although in both cases the approximations on the unknown fluxes on the contours of the nodes are constant.

Figure 5:
Illustration of the relative error of Problem 3.

3.4 Problem 4

In order to test the potential of the proposed methodology, we solve a representative problem as described by Ref. ⁹, of the beginning of life core of Unit 1 of the KOEBERG nuclear plant¹. It is a four groups problem with upscattering and downscattering. The realistic cross sections and other parameters are given in Ref. ⁹. The core is composed of 221 homogeneous square assemblies with 21.608 cm of width and seven different material types with zero flux conditions at the boundaries. Due to the symmetry conditions of the problem, we redefine the problem for a quarter of the core, as we can see in Figure 6. Thus, zero flux is assumed at the top and right boundaries while reflection is considered on the left and bottom boundaries.

Figure 6:
Geometry of Problem 4 with four groups.

Muller and Weiss (1991) ⁹ solved the two-dimensional diffusion problem by the high-order response matrix method, where the currents on nodal interfaces are approximated by expansions in terms of Legendre polynomials. In each node, the diffusion equation is solved via the Fourier series method, generating the nodal response matrices. In Table 8 we present the numerical results for the effective multiplication factor generated by our method and by the method of Muller and Weiss ⁹. We generate our results using ∆x = 10.804/m and ∆y = 10.804/m while in Ref. ⁹ they work with expansions in Legendre Polynomials of different orders, as listed in the Table 8.

Table 8:
Comparison of effective multiplication factor Problem 4.

Legendre Order	Ref. 9	m	Constant Source	t (s)	Linear Source	t (s)
0	1.001457	1	1.0004	0.404	1.0040	0.619
2	1.007959	2	1.0052	5.544	1.0064	6.202
4	1.007954	3	1.0065	31.24	1.0071	32.91
6	1.007954	4	1.0071	143.78	1.0074	145.65

In Table 8, we reported the effective multiplication factor for various mesh configurations and we can observe that as the mesh is refined the K _eff agree up to 4 digits with Ref. ⁹. It is worth mentioning that, differently of Ref. ⁹, we use only constant approximations for the currents (transverse leakages) on nodal interfaces in Eq. (2.11). The results were generated with a low computational cost to obtain up to 4 digits with the reference adopted. In addition, until this moment no special attention has been given to optimizing the computational procedure, for instance, exploring the sparsity of the matrices when solving the linear systems. In Figure 7, we present the graph of the relative error of K _eff produced by the constant and linear approximations of the source, we consider as reference K _eff = 1.007954. In this figure, we can see that the linear approximations generate more accurate results than the constant approximations, mainly in thicker meshes.

Figure 7:
Illustration of the relative error of Problem 4.

4 CONCLUDING REMARKS

The effective multiplication factor and the neutron flux in multigroup diffusion problems were evaluated through a scheme based on a source iteration method and a nodal technique. The decoupling of the system, provided by the source iteration method, and the assumptions of constant transverse leakage approximations and linear representation for the fluxes in the source term, allowed the derivation of analytical expressions for the averaged fluxes in each direction.

In addition to being simple, the procedure is robust since it was successfully applied to solve heterogeneous media multigroup problems with up and downscattering for different geometric configurations in Cartesian geometry. This procedure may be extended for multidimensional problems since the transversal integration always leads to one-dimensional equations that may be solved analytically. It certainly enhanced the previously proposed methodology ¹⁸ based on the use of Fourier Transforms.

With low computational cost, several reference problems available in the literature were solved in this work. The linear representation of the fluxes in the source term produced better results than constant approximations, mainly in coarser meshes. And we observed that the increase in the order of approximation of the term source does not generate significant growth in computational cost. Thus, this strategy proposed in this work appears to be an alternative rather than exploring higher-order approximations to express the leakage terms on the boundary of the nodes, as usual in the literature associated with nodal methods.

Surely, either higher-order expansions or alternative forms of approximations may be explored to represent the transverse leakage on the edges of the nodes and the fluxes of the source term, which may lead to better results in coarser meshes. It is a natural extension for future works. However, it possibly will imply either more extensive analytical derivations or computational cost. The main aspect we found relevant to emphasize reporting the present results was the generality and performance obtained with such an approach. The fact of keeping some analyticity in the formulation is also something to be noted. In this way, we obtained the preliminary endorsement to associate, improve and extend this formalism as part of the solution of time-dependent problems, which is our main goal.

Acknowledgement

This study was financed in part by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) Brasil, the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) Brasil (Finance Code 001) and the Fundação de Amparo à Pesquisa do Estado do Rio Grande do Sul (FAPERGS): (Proc: 19/2551-0001766-9).

REFERENCES

• ¹
  A. Bernal, J.E. Roman, R. Miró & G. Verdú. Calculation of multiple eigenvalues of the neutron diffusion equation discretized with a parallelized finite volume method. Progress in Nuclear Energy, 105 (2018), 271-278. doi:https://doi.org/10.1016/j.pnucene.2018.02.006.
  » https://doi.org/https://doi.org/10.1016/j.pnucene.2018.02.006
• ²
  A. Carreño, A. Vidal-Ferràndiz, D. Ginestar & G. Verdú . Block hybrid multilevel method to compute the dominant λ-modes of the neutron diffusion equation. Annals of Nuclear Energy, 121 (2018), 513-524.
• ³
  Z. Chunyu & C. Gong. Fast solution of neutron diffusion problem by reduced basis finite element method. Annals of Nuclear Energy , 111 (2018), 702-708. doi:https://doi.org/10.1016/j.anucene.2017. 09.044.
  » https://doi.org/https://doi.org/10.1016/j.anucene.2017. 09.044
• ⁴
  J.J. Duderstadt & L.J. Hamilton. “Nuclear reactor analysis”. John Wiley, New York (1976).
• ⁵
  N. Gupta. Nodal methods for three-dimensional simulators. Progress in Nuclear Energy , 7(3) (1981), 127-149.
• ⁶
  A. Hébert. A simplified presentation of the multigroup analytic nodal method in 2-D Cartesian geometry. Annals of Nuclear Energy , 35(11) (2008), 2142-2149.
• ⁷
  S.A. Hosseini. Development of Galerkin finite element method three-dimensional computational code for the multigroup neutron diffusion equation with unstructured tetrahedron elements. Nuclear Engineering and Technology, 48(1) (2016), 43-54.
• ⁸
  R. Lawrence. Progress in nodal methods for the solution of the neutron diffusion and transport equations. Progress in Nuclear Energy , 17(3) (1986), 271-301.
• ⁹
  E. Muller & Z. Weiss. Benchmarking with the multigroup diffusion high-order response matrix method. Annals of Nuclear Energy , 18(9) (1991), 535-544.
• ¹⁰
  H. Sekimoto. “Nuclear reactor theory”. Tokyo Institute of Technology Press, Tokyo (2007).
• ¹¹
  R. Shober & A. Henry. “Nonlinear methods for solving the diffusion equation”. MIT Report NE-196 (1976).
• ¹²
  A.C. Silva, A.S. Martinez & A.C. Gonçalves. Reconstruction of the neutron flux in a slab reactor. World Journal of Nuclear Science and Technology, 2 (2012), 181-186.
• ¹³
  K.S. Smith. “An analytic nodal method for solving the two-group, multidimensional, static and transient neutron diffusion equations”. Ph.D. thesis, Massachusetts Institute of Technology, Department of Nuclear Engineering, Cambridge, MA (1979).
• ¹⁴
  W.M. Stacey. “Nuclear reactor physics”. John Wiley & Sons, New York (2001).
• ¹⁵
  M. Vagheian, D.R. Ochbelagh & M. Gharib. A new moving-mesh finite volume method for the efficient solution of two-dimensional neutron diffusion equation using gradient variations of reactor power. Nuclear Engineering and Technology , 51(5) (2019), 1181-1194.
• ¹⁶
  M. Vagheian , N. Vosoughi & M. Gharib . Enhanced finite difference scheme for the neutron diffusion equation using the importance function. Annals of Nuclear Energy , 96 (2016), 412-421.
• ¹⁷
  W. Wu, Y. Yu, Q. Luo, D. Yao, Q. Li & X. Chai. Calculation of higher eigen-modes of the forward and adjoint neutron diffusion equations using IRAM algorithm based on domain decomposition. Annals of Nuclear Energy , 143 (2020), 107463.
• ¹⁸
  R. Zanette, C.Z. Petersen, M. Schramm & J.R. Zabadal. A modified power method for the multilayer multigroup two-dimensional neutron diffusion equation. Annals of Nuclear Energy , 111 (2018), 136-140.
• ¹⁹
  X. Zhou & F. Li. General nodal expansion method for multi-dimensional neutronics/thermal-hydraulic coupled problems in pebble-bed core systems. Annals of Nuclear Energy , 116 (2018), 10-19.

• 1
  PWR reactor situated on the west coast of South Africa.

Publication Dates

• Publication in this collection
  27 June 2022
• Date of issue
  Apr-Jun 2022

History

• Received
  04 Dec 2020
• Accepted
  23 Nov 2021