2.1 Even parity discrete ordinates method for neutron transport equation

To provide a detailed description of the angular flux nonlinear finite difference equation, we begin with the multigroup neutron transport equation in the Cartesian coordinate system, as follows:

where vector r is the spatial position relative to the reference point, and vector Ω is the unit direction of flying with respect to the reference. Subscript g denotes the energy group index using multigroup approximation. This equation can represent the steady-state or transient neutron transport equation. When it represents a transient equation, the time variable is eliminated for brevity, as follows:

where S_s, S_f, S_d, and S_e are the scattering, fission, delayed, and external sources, respectively. The external source is often disregarded, except when it must be considered. The fission and delayed neutron sources are isotropic. For the steady state, the time derivative term is zero, and the delayed neutron source is included in the fission source. If we discretize the time derivative term $\frac{\partial }{\partial t}{\psi }_g\left(r,\Omega \right)$ with fluxes at the current and last time steps, the transient equation will become a transient fixed source equation, which can be solved in a manner similar to the steady-state equation. Therefore, we will only discuss the steady-state equation, and Eq. (2a) can be expressed as

The scattering and fission source can be expressed in terms of the angular flux as follows:

where ψ, φ, Σ, χ, μ₀, and k_eff are the angular flux, scalar flux, macro cross sections, fission spectrum, cosine of scattering angle in the center of mass system, and eigenvalue with a standard definition in nuclear reactor physics, respectively. Using the discrete ordinates method, we can represent the angular space with selected angles as follows:

where ${\Omega }_i$ and ${\omega }_i$ are the selected angles and weights, respectively, and are known as the S_N quadrature set; M is the total number of discrete directions. Here, we use subscript i to represent quantities at angle ${\Omega }_i$ and subscript -i to represent quantities at angle $-{\Omega }_i$. The equations for angular fluxes in opposite directions are as follows:

The even and odd parity fluxes are expressed as follows:

The even and odd parity fluxes exhibit the following properties:

Using Eqs. (4) and (5), we can obtain the even and odd parity transport equations as follows:

where $S_{g,i}^{+}\left(r\right)$ and $S_{g,i}^{-}\left(r\right)$ are the even and odd parity sources, respectively, expressed as

As the fission sources are isotropic, we have

where $S_{\text{s},g,i}^{+}\left(r\right)$ and $S_{\text{s},g,i}^{-}\left(r\right)$ are the even and odd parity scattering sources, expressed as

Here, $S_{\text{s},g,i}\left(r\right)$ is the angle-dependent scattering source, expressed as

μ₀ is the cosine of the scattering angle in the center of mass system, where ${\mu }_0=\Omega \cdot {\Omega }^{\prime }$. With a P_n expansion of the scattering cross sections, we obtain

where P_k are the Legendre polynomials and ${\Sigma }_{s,k}$ are P_k scattering moments. By expanding the even and odd parity scattering sources with a P₀ expansion, we have

Combined with Eqs. (9) and (13), the even and odd parity sources can be expressed as follows:

Integrate the even parity transport equation in Eq. (7a) over a region with constant a cross section as follows:

where subscript r is the region index; $V_r$ is the volume of region r; $\Gamma$ is the surface of region r; $\widehat{n}$ is the unit normal direction vector of the surface; ${\overline{\overline{\psi }}}_{g,i,r}^{+}$ and ${\overline{\overline{S}}}_{g,i,r}^{+}$ are the node average even parity flux and node average even parity source weighted by the region volume, respectively. ${\Sigma }_{\text{t},g,i,r}$ is the angle-dependent total macro cross section, defined as follows:

It is noteworthy that the divergence theorem was used for deriving Eq. (15).

The surface of region r is divided into n_r segments, where n_r is the number of neighbors. Each segment is the interface between region r and one of its neighbors. Define the average odd parity flux on each segment as follows:

where n is the neighbor node of node r, and ${\Gamma }_{rn}$ is the interface between r and n. The average odd parity flux on the surface can represent the net number of neutrons passing through the interface in a certain direction, and it exhibits the following property:

By rearranging Eq. (15), we can obtain the even parity flux equation as follows:

The even parity flux equation exactly preserves the neutron balance in each direction and in each coarse mesh. If the relationship between the average odd parity flux ${\overline{\psi }}_{g,i,rn}^{-}$ and the node average even parity flux ${\overline{\overline{\psi }}}_{g,i,r}^{+}$ is obtained, then Eq. (20) can be solved further. This relationship can be obtained from the odd parity transport equation.

Integrating the odd parity transport equation (Eq. (7b)) over each segment of the surface, we obtain

The left hand of Eq. (21) is the average derivative of the even parity flux on interface ${\Gamma }_{rn}$ in the ${\Omega }_i$ direction and can be defined as follows:

Substituting Eq. (22) into Eq. (21), we have

Through a finite difference approximation for ${\overline{\psi }}_{g,i,rn}^{+}{}^{\prime }$, we obtain the following odd parity flux equation:

where ${\overline{\psi }}_{g,i,rn}^{+}$ is the average even parity flux on the surface between r and n, and $h_{i,rn}$ is the thickness of the node center to the surface center in the direction of angle ${\Omega }_i$.

The average even parity flux on the surface exhibits following priority:

2.2 Definition and application of discontinuity factor

To render all ${\overline{\overline{\psi }}}_{g,i,r}^{+}$, ${\overline{\psi }}_{g,i,rn}^{+}$ and ${\overline{\psi }}_{g,i,rn}^{-}$ in the odd parity flux equation (Eq. (24)) satisfy the MOC transport calculation, we introduced a factor known as the even parity discontinuity factor, $f_{g,i,rn}^{+}$, as follows:

$f_{g,i,rn}^{+}$ will be greater than zero when the denominator of Eq. (27) is positive. However, the denominator can less than or equal to zero. In this case, $f_{g,i,rn}^{+}$ will be infinite or negative. To avoid this, we introduced another discontinuity factor to Eq. (27), i.e., the odd parity discontinuity factor, $f_{g,i,rn}^{-}$.

An algorithm was designed to determine the even and odd parity discontinuity factors. This algorithm begins by evaluating the following conditional:

If the condition in Eq. (29) is true, then the even and odd parity discontinuity factors can be calculated from the MOC solutions as follows:

If the condition is false, then different equations are designed to calculate the even and odd parity discontinuity factors as follows:

Using this approach, the even and odd parity discontinuity factors will be positive at all times. Rearranging Eq. (28) yields

If neighbor node n is within the computational domain, then a similar equation is obtained for the other side of the interface ${\Gamma }_{rn}$, as follows:

Combining Eqs. (32) and (33) and the continuity condition ${\overline{\psi }}_{g,i,rn}^{+}={\overline{\psi }}_{g,i,nr}^{+}$, ${\overline{\psi }}_{g,i,rn}^{-}=-{\overline{\psi }}_{g,i,nr}^{-}$, we obtain

As for the vacuum boundary surface, the incident angular flux is zero. Therefore, we have the equation ${\overline{\psi }}_{g,i,rn}^{-}={\overline{\psi }}_{g,i,rn}^{+}$, which can be derived from Eqs. (18a) and (25a). Combining this equation with Eq. (32), we obtain

For the reflective boundary surface, the angular fluxes in the reflect directions are equal. Therefore, we have the continuity condition ${\overline{\psi }}_{g,i,rn}^{+}={\overline{\psi }}_{g,i*,rn}^{+}$ and ${\overline{\psi }}_{g,i,rn}^{-}=-{\overline{\psi }}_{g,i*,rn}^{-}$, where subscript i* and i represent a pair of angles in the incident and reflective directions, respectively. Fig. 1 shows an illustration of the reflective boundary condition.

Fig. 1.Full-screen View

Illustration of reflective boundary condition

Combined with Eq. (32) for both the angle i and i*, we obtain

In general, the relationships between the surface average odd parity flux and node average even parity flux for both the interior surface and boundary surface are established through Eqs. (34), (35), and (36). Moreover, these equations satisfy the high-order MOC solution through the strict definitions of the even and odd parity discontinuity factors.

2.3 Construction of angular flux nonlinear finite difference linear system

Based on Eqs. (34), (35), and (36), the odd parity fluxes on the interior and boundary surfaces can be expressed as a linear function of the neighbor average even parity fluxes using the coupling coefficient C.

For the interior surface, we have

where

For the vacuum boundary surface, we have

where

For the reflective boundary surface, we have

where

Substituting Eqs. (37)–(39) into Eq. (20), we obtain the following angular flux nonlinear difference linear system:

If the homogenized cross sections and coupling coefficients are obtained from the results of the MOC transport calculation, then the ANFD linear system can be solved. As the even parity fluxes in opposite directions are equal, only the ANFD linear system need to be solved on a half angular space.

3.1 Decoupling of energy group and angle

To solve the linear system by group and angle separately, we would like to keep $C_{g,i,rn}=C_{g,i,nr}$ in the interior surfaces such that the linear system will be symmetric. A symmetric linear system is much easier to solve than an asymmetric system. To create a symmetric ANFD linear system, the coupling coefficients for the interior surface are redefined as follows:

where ${\bcancel{\psi }}_{g,i,rn}$ is the asymmetric term of Eq. (42a), which can be determined using even parity fluxes. Substituting Eq. (42) into Eq. (40) and placing the asymmetric term to the right side of the equation, the ANFD linear system is rearranged as follows:

(43)

where N, S, W, E, T, and B represent the north, south, west, east, top, and bottom surfaces of the node, respectively. The superscript k denotes the power iteration number. As the symmetric and diagonal dominant matrix is a symmetric positive definite matrix, Eq. (43) will become an expression for a symmetric positive definite linear system. Furthermore, the preconditioned conjugate gradient (CG) method ^[21] can be used to solve this linear system.

It is noteworthy that decoupling the angle and energy group will increase the iteration number of the ANFD method. However, the ANFD computing costs are low compared with those of MOC calculations, and they can be reduced significantly using the multilevel acceleration method, which has been implemented in several established core simulators ^[22-24]. The main idea of the multilevel acceleration method is to use a few-group CMFD linear system, such as a one- or two-group CMFD, to reduce the computing burden of the MG CMFD linear system. The few-group CMFD linear system is used to provide an effective estimation of the eigenvalue at the beginning of each power iteration, thereby significantly reducing the outer iteration number of the MG CMFD linear system. Therefore, a new approach similar to the multilevel acceleration method can be used to reduce the iteration number of ANFD required for convergence. This will be investigated further in future studies.

Meanwhile, if we wish to solve the linear system with all angles coupled or all energy groups coupled, it would be meaningless to maintain $C_{g,i,rn}=C_{g,i,nr}$. In this case, we will define the interior surface coupling coefficients directly from Eq. (37), and the preconditioned CG method will not be used.

3.2 Preconditioned Krylov subspace method

The direct solution of this linear system is unpractical as the matrix can become highly ill-conditioned, causing problems in the numerical solution. The Krylov subspace method has been proven to be extremely effective for solving large, nonsymmetric, sparse linear systems, particularly when a good preconditioner is applied. In this study, the right preconditioned generalized minimal residual (GMRES) method was used to solve the ANFD linear system. The right preconditioned GMRES algorithm can be summarized as follows:

Algorithm 1:Full-screen View

Right Preconditioned GMRES

Meanwhile, if we decoupling the energy group and angle to construct a symmetric positive definite linear system, then the preconditioned conjugate gradient method can be used to solve the ANFD linear system. The left preconditioned conjugate gradient algorithm can be summarized as follows:

Algorithm 2:Full-screen View

Left Preconditioned CG

As for precondition technology, a self-developed RSILU preconditioner has been used to improve parallel computing capability and minimize computational burden. The RSILU preconditioner does not require a multicolor ordering strategy for parallel computing and can provide convergence rates comparable to the standard incomplete lower–upper preconditioner. More details are available in the literature ^[25].

3.3 Iteration scheme of ANFD method

Once the node and surface angular flux obtained from the MOC calculation become available, the ANFD homogenized cross sections as well as the even and odd parity discontinuity factors can be easily obtained using Eq. (16) and (30), respectively. The strategy of decoupling the energy group and angle was used to construct the ANFD linear system, and the coupling coefficients were calculated using Eqs. (37)–(39). Hence, an equivalent ANFD linear system can be established. It is noteworthy that the asymmetric term arising from the construction of a symmetric positive definite matrix should be updated with the even party angular flux when solving the ANFD linear system. The right preconditioned GMRES solver was used to solve the ANFD linear system in this study. The exit from the GMRES solver to the ANFD cycle occurs when the required reduction in the relative L2-norm of the residual vector is achieved.

After solving the ANFD linear system, new even parity flux distributions and eigenvalues were obtained. To update the 2D MOC parameters, the average scalar flux on each node and the average angular flux on each surface were reconstructed as follows, indicated as ${\overline{\varphi }}_{g,r}$ and ${\overline{\psi }}_{g,i,rn}$, respectively:

The former was used to update the scalar flux in the flat source regions, whereas the latter was used to update the incident angular flux on the boundary surface. The M in Eq. (44a) denotes the total number of discrete directions, which is consistent with the number of MOC azimuthal angles. After updating the average node scalar flux ${\overline{\varphi }}_{g,r}$ and average surface angular flux ${\overline{\psi }}_{g,i,rn}$, 2D MOC calculations were performed, and the cell-averaged cross sections as well as the even and odd parity discontinuity factors were updated. Subsequently, ANFD was performed using the updated cross sections as well as the even and odd parity discontinuity factors. This process was repeated until all convergence criteria were satisfied. The iteration scheme of the ANFD method is depicted in Fig. 2.

Fig. 2.Full-screen View

Iteration scheme of ANFD method

4.1 2D pin cell problem

As indicated in Fig. 3(a), the side length of the cell was 1.26 cm, and the radius of the fuel pins was 0.54 cm for the 2D pin cell problem. The cross sections data for the UO2 fuel and moderator were obtained from the 2D C5G7 benchmark report. A vacuum boundary was applied to the east and south boundary surfaces, and the reflected boundary was applied to the north and west boundary surfaces. Three moderator rings and five fuel rings with 16 azimuthal divisions were selected for 128 flat-source regions.

Fig. 3Full-screen View

(Color online) Layout of 2D pin and 2D assembly problems

For the 2D pin cell problem, Fig. 5(a) shows the angular distribution of the even parity flux in the first energy group, and Fig. 5(b) shows the energy group distribution of the even parity flux in the first azimuth angle. It is clear that the even parity flux distribution of the ANFD was consistent with the MOC solution in both the angle and energy groups. Table 1 shows the eigenvalue difference and overall even parity flux relative error. The difference in the eigenvalue was 0 pcm, and the flux maximum relative error was 3.683E-06.

Fig. 5.Full-screen View

Normalized distribution of even parity flux in azimuth angle and energy group for 2D pin cell problem

4.2 2D assembly problem

As indicated in Fig. 3(b), the 2D assembly problem comprised a 17 × 17 lattice of square pin cells containing a UO2 fuel pin cell, guide tube, and fission chamber. The cross sections data for these materials were obtained from the 2D C5G7 benchmark problem. The division of the flat-source region was the same as that of the 2D pin cell problem. As for the boundary conditions, vacuum boundary conditions were applied to the east and south boundary surfaces, whereas the reflected boundary conditions were applied to the north and west boundary surfaces.

Fig. 6 shows the space relative error of the even parity flux in the first energy group and first azimuth angle for the 2D assembly problem. As shown, the even parity flux distribution of the ANFD indicated a relatively good agreement in space compared with the MOC solution. Table 2 shows the eigenvalue difference and overall even parity flux relative error for the 2D assembly problem. The difference in the eigenvalue was 0 pcm, whereas the flux maximum relative error was 1.162E-05.

Fig. 6.Full-screen View

Even parity flux distribution relative error in first energy group and first azimuth angle for 2D assembly problem

4.3 2D C5G7 problem

The 2D C5G7 benchmark problem was originally a pin-cell level heterogeneous benchmark for testing the accuracy of modern deterministic transport methods and codes without spatial homogenization. The 2D C5G7 problem configuration is shown in Fig. 4. The overall geometry of this problem was 64.26 cm × 64.26 cm, and each assembly measured 21.42 cm × 21.42 cm. Each fuel assembly comprised a 17 × 17 lattice of square pin cells, including a moderator, guide tube, fission chamber, MOX fuel pin cell with different enrichment, and UO2 fuel pin cell. For the reflector pin cases, the six additional pin cells of each reflector assembly used 1 × 1 sub mesh, whereas the 11 pin cells near the fuel used a 6 × 6 Cartesian sub mesh. For the other pin cells of the fuel and control rods, 128 flat-source regions comprising five fuel rings and three moderator rings with 16 azimuthal divisions were used. An accurate Monte-Carlo reference solution was provided in the 2D C5G7 benchmark report, which contains a precise eigenvalue solution and normalized pin power (NPP) that can be used as a reference solution. The NPP ($p_r$) in each fuel pin is defined as follows:

Fig. 4Full-screen View

(Color online) Layout of 2D C5G7 benchmark problem

where $R_{\text{f},r}$ is the fission rate in each fuel pin, ${\overline{R}}_{\text{f}}$ the average fission rate of all fuel pins, and $N$ the number of all fuel pins.

Figure 7 indicates the percent error of the NPP for the 2D C5G7 benchmark problem; as shown, the main error appeared in the fuel pin cell near the reflector.

Fig. 7Full-screen View

(Color online) Percent error of pin power for 2D C5G7 benchmark problem

Table 3 shows the eigenvalue difference and NPP error from our results and the MCNP reference solutions. Compared with the reference results, the error of the eigenvalue was 3 pcm, which was less than the stochastic uncertainty of the MCNP reference results. In general, the eigenvalue, NPP, etc. agreed well with the reference results.

4.4 Void region problem

The ANFD method was developed from transport theory and does not require the calculation of the diffusion coefficient. Therefore, the ANFD method may exhibit a better acceleration than the CMFD method in certain strongly anisotropic cases, such as those involving void regions and anisotropic scattering. In this study, a lattice problem involving void regions was constructed and tested; 5 × 5 square pin cells were set to the void regions, whereas UO₂ fuel cells were used in the other regions. The geometry and cross sections of the UO₂ fuel cell were the same as those provided in the 2D C5G7 benchmark report. The reflected boundary conditions were applied to all boundary surfaces. The layout of this void region problem is shown in Fig. 8.

Fig. 8Full-screen View

(Color online) Layout of lattice case involving void regions

The convergence history of the MOC, ANFD, and CMFD are shown in Table 4 and Fig. 9. As expected, the CMFD failed to converge in this problem, whereas the ANFD converged and remained effective.

Fig. 9.Full-screen View

Convergence histories of MOC, ANFD, and CMFD

Void regions are typical in real reactors, such as TREAT ^[27] and HTR-10 ^[28]. In a previous study [29], the quasi-diffusion method was used to address the void region problem in TREAT core analysis. In this study, the quasi-diffusion method was implemented using the nodal expansion method. Meanwhile, the nodal wise average cross sections, Eddington factors, and discontinuity factors were generated from 3D whole-core Serpent Monte-Carlo simulations. By comparison, the ANFD is more concise and only requires even and odd parity discontinuity factors generated from MOC calculations.