© H. Belanger et al., Published by EDP Sciences, 2026

1. Introduction

Despite the growing prevalence of the use of Monte Carlo methods in nuclear reactor analysis, deterministic methods, particularly lattice physics codes, remain the dominant work horse for most industry LWR design calculations. Examples of such codes include APOLLO3 [1], STREAM [2], GALAXY [3], CASMO-5 [4], and DRAGON [5], just to name a few. Most lattice physics codes in use today are closed source¹. This can be a challenging limitation in a teaching and research context, as it limits students’ access to appropriate tools to perform realistic reactor calculations, and impedes graduate students from learning the details about how these methods are implemented and applied “under the hood”. Over the past decade, OpenMC has become a widely used Monte Carlo transport tool that is open source, permitting easy access to all students [6]. However, using only the Monte Carlo method, it is impractical, if not impossible, for students to perform calculations beyond the scope of a single fuel assembly using their laptops. Full core PWR studies related to depletion, fuel loading patterns, transients, and more, are therefore out of reach when one only has access to a Monte Carlo code and a laptop.

To put reactor physics tools into the hands of students, a new suite of deterministic neutron transport calculation tools has been developed, called Scarabée. This library was developed from the ground up, with students and ease of use in mind. It has an assortment of various transport and diffusion solvers which can be accessed directly. Additionally, a lattice physics tool chain is provided for PWR fuel assembly analysis. The objective of this work is to present the current features and development status of the Scarabée library, highlighting the architecture of the code, the methods employed, and some example calculations.

2. Software and architecture

Scarabée has been designed to be natively multi-platform, and supports Linux, Windows, and macOS. A majority of the code base is written in C++20, using an object-oriented approach. Fundamental objects such as cross sections and transport solvers are all written as C++ classes to ensure optimal performance. General users of the library, however, are not expected to know or use any C++. Instead, all publicly accessible objects which are intended for user interaction are made accessible via a Python interface using the Pybind11 library [7]. As such, to use the library, users write an “input” file for Scarabée using standard Python syntax and the exposed Python interface for the library. This design choice was inspired by the approach used in the TRIPOLI-5 Monte Carlo code [8], and is very reminiscent of the current approach used by OpenMC [9]. In fact, some of the syntax in Scarabée, such as for building materials, has been made to be nearly identical to the OpenMC Python interface, facilitating easy translation of models.

While low level objects which represent cross sections and various transport solvers are all written in C++, a portion of the library is also written in pure Python. This includes classes and functionality related to automated PWR fuel assembly lattice calculations. These portions of the library can be written in Python as they are largely based on coupling the outputs and inputs of several different transport solvers in the correct order. Such operations are not typically bottle necks on performance, and writing them in Python allows them to be easily examined and modified without needing to recompile the library. Figure 1 provides a depiction for how the Python portion of the library is layered on top of the C++ portion of the library. Multidimensional arrays are used quite extensively in the library, on both the C++ and the Python side. In Python, the de facto standard for multidimensional arrays is Numpy [10]. To allow for seamless integration between the arrays on both the Python and C++ sides of the code, the xtensor library is used, which has native support for Pybind11 [11, 12]. Using xtensor, Numpy arrays from Python can be viewed and modified on the C++ side without even needing to make a copy of the underlying data, allowing for in-memory modifications to data arrays.

Fig. 1. Depiction of how Pybind11 is used to interface low-level C++ classes in Scarabée with the higher level features written in Python.

Other external libraries which are used in Scarabée include Eigen [13] and HDF5 [14]. Eigen is a linear algebra library which is used for solving all linear systems on the C++ side. It also provides a sparse matrix representation which keeps memory usage at reasonable levels when using certain transport solvers. As will be outlined in Section 3.1, microscopic nuclear data for Scarabée are stored in an HDF5 file, so the HDF5 library is required to read this information. Because Python on Windows now requires the user to give explicit permission to load dynamic libraries, Scarabée compiles and statically links all dependencies to prevent errors and avoid problems with incompatible versions which might already be on a given system.

The project is hosted on GitHub where the community may download and review the source code [15]. The cibuildwheel tool is used to automatically compile Python wheels of Scarabée for all supported platforms and Python versions as an automatic GitHub Action when new versions are released. As such, users do not even need to have a compiler present on their system, and can simply install it with pip, the standard Python package manager. Currently Python versions 3.10 through 3.14 are supported; only living versions of Python are supported by the project. External users are welcome to contribute, as the goal is for Scarabée to become a community developed and maintained tool for nuclear engineering education.

3. Methods and theory

3.1. Nuclear data and self-shielding

3.1.1. Nuclear data library

The cross section library in Scarabée uses a custom formatted HDF5 file, partially inspired by XPZLIB [16]. Each nuclide in the library is contained in a single group² within the HDF5 file. Information such as the nuclide’s atomic weight ratio, potential cross section, and energy released per fission are stored as attributes on the group. All of the microscopic cross sections are stored in various datasets contained within the group. All nuclides contain datasets for the absorption cross section, transport correction cross section, the P0 scattering matrix, and the P1 scattering matrix, at infinite dilution and for all tabulated temperatures. If a nuclide has reactions such as fission, (n, γ), (n, 2n), (n, 3n), (n, p), or (n, α), these may also have datasets present with these cross sections at infinite dilution tabulated as a function of temperature. The naming convention for these datasets is the prefix “inf” followed by the reaction name (i.e. inf-absorption, inf-scatter, or inf-(n,gamma)). Higher scattering moment matrices may be provided in the library. If the P3 scattering matrix is provided, then there would be an inf-p3-scatter dataset provided in the nuclide’s group. Scattering matrices are stored in a compressed format, similar to the MATXS format [17]. Certain reactions, such as absorption, scattering, fission, and (n, γ), may exhibit resonant behavior for certain nuclides (i.e. ²³⁸U). When this is the case, an additional dataset is provided for that reaction with the “res” prefix, where the cross sections are stored as a function of temperature and dilution. The reactions (n, 2n), (n, 3n), (n, p), and (n, α) are never assumed to be resonant. If a resonant table is provided for a given nuclide and reaction, it only provides data for the energy groups which are in the resonance region of the library, and is dependent on the chosen group structure. Figure 2 provides a depiction of what the HDF5 file structure looks like for a single nuclide. The Scarabée source repository contains several scripts which can automate the generation of a minimal nuclear data library for ENDF/B-VIII.0 in various group structures [18]. These scripts are not distributed in the Python wheel for Scarabée, as they require other dependencies such as ENDFtk [19], PapillonNDL [20], and FRENDY [21, 22]. Currently, nuclear data libraries can only be generated for Scarabée using the FRENDY nuclear data processing code. This choice was made due to FRENDY’s ease of use and advanced multi-group processing capabilities [22, 23]. While the use of another data processing code is theoretically possible, no infrastructure currently exists to support other codes. Therefore, most users should use FRENDY to generate their cross section libraries.

Fig. 2. Depiction of the HDF5 file structure for the microscopic nuclear data used by Scarabée.

3.1.2. Transport correction

The transport correction is a common approximation that is used in many neutronics solvers to treat the effects of anisotropic scattering. In general, the total cross section and the scattering matrix are modified by subtracting a correction term:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{tr,g}& ={\mathrm{\Sigma }}_{t,g}-\mathrm{\Delta }{\mathrm{\Sigma }}_{tr,g},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{str,g^{\prime }\to g}& ={\mathrm{\Sigma }}_{s,0,g^{\prime }\to g}-{\delta }_{g^{\prime },g}\mathrm{\Delta }{\mathrm{\Sigma }}_{tr,g},\hfill \end{array}$$(2)

where Σ_t, g is the total cross section in group g, Σ_tr, g is the transport cross section, Σ_{str, g′→g} is an element of the transport corrected scattering matrix, Σ_{s, 0, g′→g} is the P0 moment of the scattering matrix, and ΔΣ_tr, g is the transport correction cross section. After application of this correction, calculations are performed assuming isotropic scattering [24]. Traditionally, the transport correction cross section in many codes is defined using the “outflow” approximation where

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{tr,g}=\sum _{g^{\prime }=0}^G{\mathrm{\Sigma }}_{s,1,g\to g^{\prime }}=\overline{\mu }{\mathrm{\Sigma }}_{s,g},\end{array}$$(3)

Σ_{s, 1, g → g′} being the P1 scattering moment when scattering from group g to group g′, and $\overline{\mu }$ is the average cosine of the scattering angle for all scatter events in group g [25]. Previous work by Choi and Smith have indicated that the use of the outflow approximation is quite inadequate for LWR analysis, particularly for modeling neutron scattering off of ¹H bound in water [26, 27]. As an alternative, they recommend using the “inflow” transport correction, defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{tr,g}=\frac{{\displaystyle \sum _{g^{\prime }=0}^G{\mathrm{\Sigma }}_{s,1,g^{\prime }\to g}{\varphi }_{1,g^{\prime }}}}{{\varphi }_{1,g}},\end{array}$$(4)

where ϕ_1, g is the P1 flux moment in group g.

As the outflow and inflow transport corrections are quite different from one another, Scarabée has a unique feature where the transport correction cross section is tabulated for each nuclide, with both an infinite dilution table and a resonant table (if the nuclide is resonant). In the library preparation scripts provided with Scarabée, the outflow transport correction is used for all nuclides, with the exception of ¹H in water and ²H in heavy water. For these two nuclides, the approach taken by Herman is employed, where the homogeneous P1 equations are solved using the fission spectrum of ²³⁵U as a source [28]. The computational strategy outlined by Stamm’ler and Abbate is used to solve the homogeneous P1 equations, allowing the inflow transport correction cross section from equation (4) to be obtained [25].

3.1.3. Intermediate resonance parameters

In the field of neutron slowing down and resonance self-shielding, two approximations are commonly employed: the narrow resonance (NR) approximation and the wide resonance (WR) approximation [29, 30]. Both permit one to obtain an analytical solution for the neutron slowing down equation in an infinite homogeneous medium. In the NR approximation, it is assumed that all non-absorbing nuclei (the moderators) are very light, which means that the average width of a resonance is much narrower than the average amount of energy a neutron will lose in a scatter. In the WR approximation, the absorbing (resonant) nuclei are considered to be infinitely massive, so a neutron loses no energy during a scattering event so the resonances appear to be very wide. Practically speaking, whether a nuclide appears narrow or wide depends on the width of the individual resonance, the mass of the nuclide with which the neutron is scattering, and the incident energy of the neutron. This effect of nuclides being neither perfectly “narrow” nor “wide” can be accounted for by scaling its potential scattering cross section by a factor λ ∈ [0, 1], known as the intermediate resonance parameter. A value of λ = 0 corresponds to the WR approximation (the nuclide does not contribute to slowing down) and λ = 1 corresponds to the NR approximation. Only ¹H has λ = 1, and this is due to the fact that a neutron can lose all of its energy in an elastic scatter with this isotope.

As previously mentioned, λ will depend on the nuclide with which the scattering is occurring, the energy of the incident neutron, and the width of the particular resonance in question. We could never tabulate λ for all three of these parameters. Instead, we follow the suggestions of Knott, Yamamoto, and Gibson, only tabulating them for the nuclide, i, and the incident neutron energy group g [29, 30]. Since ²³⁸U is typically the most important resonant isotope in an LWR, we will calculate λ_i, g using the ²³⁸U cross section in group g. The method used by Scarabée to determine λ_i, g is known as the hydrogen equivalence procedure, which is well outlined by Xu, Gao, and Downar [31]. First, a dilution table is generated for the (n, γ) cross section of ²³⁸U, where the only other add-mixed moderator is ¹H. When creating this dilution table, the slowing-down equation is numerically solved on a log-spaced 500,000 point energy grid over the interval of [10 μeV, 20 MeV]. We now consider a reference dilution, σ_0, ref, taken to be 50 b [29]. At this reference dilution, the ratio of the ¹H number density to that of ²³⁸U would be

$$\begin{array}{c}\hfill R_{\mathrm{ref}}=\frac{{\sigma }_{0,\mathrm{ref}}}{{\sigma }_{p,H1}}·\end{array}$$(5)

Subsequently, one additional slowing down calculation is performed where both ¹H and nuclide i are add-mixed with the ²³⁸U, assuming that σ_t(E)=σ_p for each of these add-mixed isotopes. The fractions attributed to ¹H and i are determined from R_ref, where the ratio of ¹H to ²³⁸U is taken to be R_H1 = 0.95R_ref, and the ratio of nuclide i to ²³⁸U is taken to be R_i = 0.05R_ref. Once the new grouped cross sections are obtained from this flux spectrum, the capture cross section σ_c, g is used to search in the previously generated ¹H-only dilution table to find what equivalent hydrogen dilution, σ_0, equiv(σ_c, g), is needed to obtain σ_c, g. Once known, the intermediate resonance parameter may be calculated using

$$\begin{array}{c}\hfill {\lambda }_{i,g}=\frac{{\sigma }_{0,\mathrm{equiv}}({\sigma }_{c,g})-R_{H1}{\sigma }_{p,H1}}{R_i{\sigma }_{p,i}}·\end{array}$$(6)

The astute reader might notice that equation (6) does not necessarily enforce λ_r, g ∈ [0, 1]. Sanchez has pointed out that while one might like this condition to be met (to correspond with the NR and WR approximations as the extrema), there is no mathematical reason why this must be the case [32]. If one were to consider the intermediate resonance parameter as a “glorified fudge factor” which allows one to retrieve the correct resonance integral, then it could indeed be any real quantity. At the moment, Scarabée enforces λ_r, g ∈ [0, 1], and sets λ_r, g = 1 for groups where this was not the case. If, after this process, there are resonance groups with λ_r, g = 1, then this procedure is repeated with ²³⁵U, as it is assumed that the reason for an invalid value (resulting in λ_r, g = 1) was potentially due to a lack of ²³⁸U resonances within group g. Any groups where λ_r, g = 1 for ²³⁸U are replaced with the λ_r, g value obtained using ²³⁵U (which might also be unity). In the future, the current choice of enforcing λ_r, g ∈ [0, 1] will be reexamined.

3.1.4. Self-shielding

Like many other lattice physics codes, Scarabée uses equivalence in dilution for determining appropriate multi-group cross sections. By considering the slowing down equation for neutrons in a cylindrical fuel rod surrounded by moderator, one is able to show that after applying the Narrow Resonance (NR) approximation, the slowing down equation becomes

$$\begin{array}{cc}\hfill V_f{\mathrm{\Sigma }}_{t,f}(E){\varphi }_f(E)& =\frac{1}{E}[(1-P_{f\to m}(E))V_f{\mathrm{\Sigma }}_{p,f}\hfill \\ \hfill & +P_{m\to f}(E)V_m{\mathrm{\Sigma }}_{p,m}],\hfill \end{array}$$(7)

where V_f and V_m are the volume (area) of the fuel and moderator, Σ_p, f and Σ_p, m are the potential scattering cross sections for the fuel and moderator, P_f → m(E) is the probability of a neutron to fly from the fuel to the moderator and have a collision, and P_m → f(E) is the probability of the neutron to fly from the moderator to the fuel and have a collision [2, 3, 29, 30, 33]. After applying the reciprocity relation from collision probabilities (provided in Eq. (16)), we can express P_m → f(E) in terms of P_f → m(E), arriving at an expression for the flux spectrum in the fuel:

$$\begin{array}{c}\hfill {\varphi }_f(E)=\frac{1}{E}[(1-P_{f\to m}(E))\frac{{\mathrm{\Sigma }}_{p,f}}{{\mathrm{\Sigma }}_{t,f}(E)}+P_{f\to m}(E)].\end{array}$$(8)

This then implies that if the proper expression for P_f → m(E) can be found, it would be possible to obtain an approximation for the flux spectrum in the fuel.

To approximate the probability of a neutron to escape the fuel into the moderator, a classic 2-term rational approximation is used such that

$$\begin{array}{c}\hfill P_{f\to m}(E)=\sum _{n=1}^2{b}_n\frac{a_n{\mathrm{\Sigma }}_e}{{\mathrm{\Sigma }}_{t,f}(E)+a_n{\mathrm{\Sigma }}_e},\end{array}$$(9)

where b_n and a_n are geometry specific quantities, and Σ_e is known as the “escape cross section”, calculated as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_e=\frac{S}{4V},\end{array}$$(10)

where S is the surface area (perimeter) of the fuel region and V is the volume (area). Equation (9) is also subject to the constraint that b₁ + b₂ = 1 [29]. Substituting equation (9) into equation (8), one finds that

$$\begin{array}{c}\hfill {\varphi }_f(E)=\frac{1}{E}\sum _{n=1}^2{b}_n\frac{{\mathrm{\Sigma }}_{p,f}+a_n{\mathrm{\Sigma }}_e}{{\mathrm{\Sigma }}_{t,f}(E)+a_n{\mathrm{\Sigma }}_e}·\end{array}$$(11)

At this point, we will consider only one resonant nuclide at a time, denoted with the r subscript, and we will assume that the cross section for all other nuclides in the fuel can be approximated as λ_i σ_p, i, where λ_i is nuclide i’s intermediate resonance parameter and σ_p, i is its potential scattering cross section. This allows us to rewrite equation (11) as

$$\begin{array}{c}\hfill {\varphi }_r(E)=\frac{1}{E}\sum _{n=1}^2{b}_n\frac{{\lambda }_r{\sigma }_{p,r}+{\sigma }_{0,n}}{{\sigma }_{a,r}(E)+{\lambda }_r{\sigma }_{s,r}(E)+{\sigma }_{0,n}},\end{array}$$(12)

where the term-dependent dilution is computed as

$$\begin{array}{c}\hfill {\sigma }_{0,n}=\frac{{\displaystyle \sum _{i\ne r}{N}_i{\lambda }_i{\sigma }_{p,i}+a_n{\mathrm{\Sigma }}_e}}{N_r}·\end{array}$$(13)

From equation (12), we can compute the multi-group cross section for any given reaction x in group g as

$$\begin{array}{c}\hfill {\sigma }_{r,x,g}=\frac{{\displaystyle {\int }_{E_g}^{E_{g-1}}{\sigma }_{x,r}(E){\varphi }_r(E)\mathrm{d}E}}{{\displaystyle {\int }_{E_g}^{E_{g-1}}{\varphi }_r(E)\mathrm{d}E}}=\frac{{\displaystyle \sum _{n=1}^2{b}_n{\sigma }_{r,x,g,n}{\varphi }_{r,g,n}}}{{\displaystyle \sum _{n=1}^2{b}_n{\varphi }_{r,g,n}}},\end{array}$$(14)

where the flux ϕ_r, g, n can be retrieved as

$$\begin{array}{c}\hfill {\varphi }_{r,g,n}=\frac{{\lambda }_{r,g}{\sigma }_{p,r}+{\sigma }_{0,n}}{{\sigma }_{r,a,g,n}+{\lambda }_{r,g}{\sigma }_{r,s,g,n}+{\lambda }_{r,g}{\sigma }_{p,r}+{\sigma }_{0,n}},\end{array}$$(15)

and σ_{r, x, g, n} is the effective cross section for resonant nuclide r for reaction x in group g, which was interpolated from the pre-generated dilution tables at a dilution of σ_0, n (computed using the group-dependent intermediate resonant parameters for all nuclides) [30, 33].

The only parameters which have not yet been defined (or obtained from the previously described dilution tables from Sect. 3.1) are the coefficients a_n and b_n for the multi-term rational approximation of P_f → m(E), provided in equation (9). Since Scarabée is written with LWR reactor analysis in mind, all fuel regions are considered to be cylindrical. This allows us to use Carlvik’s approximation for a cylindrical fuel region where b₁ = 2, b₂ = −1, a₁ = 2, and a₂ = 3 [29, 30]. While often neglected in many codes, the isotopes in the fuel cladding do exhibit resonant behavior, and more accurate results can be obtained if the cladding cross sections are also self-shielded. For these regions, Scarabée approximates the cladding as an infinite slab, in the same manner taken by the STREAM lattice physics code [2]. For this geometry, Román’s approximation for slab geometry may be used, where b₁ = 1.1, b₂ = −0.1, a₁ = 1.4, and a₂ = 5.4 [29, 30]. Both the Carlvik and Román approximations, however, assume that the resonant region is the only one in the vicinity, and that there are no other resonant regions from which neutrons could arrive. In the tightly packed lattices of a nuclear reactor, this is not the case, as a neutron could easily fly from one fuel pin to another, or from one cladding region to another. If these lattice effects are not taken into account, the predicted flux spectrum will not be valid, and the multi-group cross sections will not accurately predict the behavior of the system. To account for these effects, the Dancoff correction, C, is obtained for every single fuel and cladding region in the geometry using the neutron current method [34]. This ensures that local spectral effects unique to each pin (such as the presence of an empty guide tube or inter-assembly gap) are accounted for in the self-shielded cross sections of each fuel pin. The Dancoff correction is then used to obtain a new set of parameters to be used in equation (9), the formulas for which can be found in Gibson’s PhD thesis [30].

Throughout this approach to self-shielding, it has been assumed that one set of multi-group cross sections is being obtained for the entire fuel pin. In general, for standard UO₂ or MOX fuel, it is perfectly fine to consider the cross sections as being homogeneous in the fuel pellet. However, many modern LWR designs use Gadolinium or Erbium in fuel pins to act as a poison which is burned out over the course of the first fuel cycle. Both of these materials have very large neutron absorption cross sections, which leads to very strong spatial self-shielding effects across the fuel pellet. To account for the strong spatial self-shielding within the fuel pellet, the pin is discretized into a user provided number of equal volume rings, and each ring has its own set of self-shielded cross sections, obtained using the Stoker-Weiss method [29, 35]. The interested reader is encouraged to read the paper by Ferrer and Hykes which gives a detailed outline of the method [33].

3.2. Cylindrical collision probabilities solver

For annular 1D geometries, the method of collision probabilities is fast and effective, as the problem has no azimuthal dependence, so the integration process is very efficient. This makes it an ideal choice for studying single pin-cell problems. The implementation in Scarabée closely follows the general approach initially developed by Carlvik and well outlined by Stamm’ler and Abbate [25]. If ${\tilde{p}}_{j,i}$ is the probability that a neutron born uniformly in ring j undergoes its first collision in ring i, then the reciprocity relation states that

$$\begin{array}{c}\hfill V_j{\mathrm{\Sigma }}_{tr,j}{\tilde{p}}_{j,i}=V_i{\mathrm{\Sigma }}_{tr,i}{\tilde{p}}_{i,j},\end{array}$$(16)

where Σ_tr, j is the transport cross section and V_j is the volume/area, both in ring j. This then indicates that a symmetric collision probability could instead be defined, such that $p_{j,i}\equiv V_j{\mathrm{\Sigma }}_{tr,j}{\tilde{p}}_{j,i}$, where p_j, i = p_i, j. These symmetric collision probabilities can then be calculated as

$$\begin{array}{c}\hfill p_{i,j}={\mathrm{\Sigma }}_{tr,j}{V}_j{\delta }_{i,j}+2(S_{i-1,j-1}-S_{i-1,j}-S_{i,j-1}+S_{i,j})\end{array}$$(17)

given

$$\begin{array}{c}\hfill S_{i,j}={\int }_0^{R_i}{\mathrm{Ki}}_3({\tau }_{i,j}^{+}(y))-{\mathrm{Ki}}_3({\tau }_{i,j}^{-}(y))\mathrm{d}y,\end{array}$$(18)

where Ki₃(x) is the 3rd Bickley–Naylor function, and τ_i, j ^±(y) are the optical depths in mean free paths, as indicated by Figure 3 [25]. The integral in equation (18) is broken up such that each annulus is integrated separately:

$$\begin{array}{cc}\hfill S_{i,j}& ={\int }_0^{R_1}{\mathrm{Ki}}_3({\tau }_{i,j}^{+}(y))-{\mathrm{Ki}}_3({\tau }_{i,j}^{-}(y))\mathrm{d}y\hfill \\ \hfill & +{\int }_{R_1}^{R_2}{\mathrm{Ki}}_3({\tau }_{i,j}^{+}(y))-{\mathrm{Ki}}_3({\tau }_{i,j}^{-}(y))\mathrm{d}y\hfill \\ \hfill & +\cdots +{\int }_{R_{i-1}}^{R_i}{\mathrm{Ki}}_3({\tau }_{i,j}^{+}(y))-{\mathrm{Ki}}_3({\tau }_{i,j}^{-}(y))\mathrm{d}y.\hfill \end{array}$$(19)

Fig. 3. Depiction of how the optical depths, τi, j ±(y), are computed as a function of y in the 1D annular collision probabilities solver.

Each sub-integral over an annulus is performed using a 15 point Gauss-Kronrod quadrature. The Bickley-Naylor functions are computationally quite expensive to compute, and most lattice physics codes approximate its solution to optimize performance. In Scarabée, the range x ∈ [0, 15] has been broken into 8 regions with each one having a unique 25 term Chebyshev approximation, maintaining high accuracy and fast evaluation times. In equation (17), S_i, j = 0 if either i or j are zero. Since the symmetric collision probabilities are used, one only needs to compute p_i, j∀j ≤ i for all energy groups.

The symmetric collision probabilities are only indirectly used to compute the flux within ring i. To begin, the Matrix $\widehat{M}$ is constructed such that its elements are

$$\begin{array}{c}\hfill M_{i,j}={\mathrm{\Sigma }}_{tr,i}{V}_i{\delta }_{i,j}-\frac{{\mathrm{\Sigma }}_{str,i}}{{\mathrm{\Sigma }}_{tr,i}}{p}_{j,i},\end{array}$$(20)

where Σ_str, i ≡ Σ_{str, i, g → g} and g is the energy group being solved. The matrix $\widehat{M}$ is then used to compute the set of X_i ^k and Y_i. X_i ^k is the flux in ring i due to a uniform unit source in ring k, and Y_i is the flux in ring i due to a unit neutron current entering the external boundary of the cell with a cosine angular distribution; both sets of quantities assume a zero albedo boundary condition on the periphery of the cell. The vector of flux responses for a source in ring k can be calculated as

$$\begin{array}{c}\hfill \widehat{M}{\mathit{X}}^{\mathit{k}}=\frac{1}{{\mathrm{\Sigma }}_{tr,k}}{\mathit{p}}^{\mathit{k}},\end{array}$$(21)

where X ^k is the vector of X_i ^k entries, and p ^k is the vector of p_k, i entries. A similar matrix equation is solved for the current response vector:

$$\begin{array}{c}\hfill \widehat{M}\mathit{Y}=\mathit{\gamma },\end{array}$$(22)

where Y is the vector of of Y_i entries and γ is the vector of first-flight collision blackness values, γ_i, which can be calculated as

$$\begin{array}{c}\hfill {\gamma }_i=\frac{4}{S_B}({\mathrm{\Sigma }}_{tr,i}{V}_i-\sum _{j=1}^N{p}_{j,j}),\end{array}$$(23)

and S_B is the circumference of the cell. If there are N rings, then N + 1 systems of equations must be solved for each energy group. The same matrix is used for all systems in a given group, so it can be inverted once and then reused. In Scarabée, the Eigen library is used to solve the system of equations using a column pivoting Householder QR decomposition.

With the set of source and current responses (X_i ^k and Y_i) for black boundary conditions (zero albedo) determined, the multi-collision blackness of the cell can be computed as

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\sum _{i=1}^N{\mathrm{\Sigma }}_{r,i}{V}_i{Y}_i,\end{array}$$(24)

Σ_r, i being the removal cross section for ring i. This term is needed in computing the current and source responses for an arbitrary albedo, a, at the cell boundary:

$$\begin{array}{cc}\hfill Y_i(a)& =\frac{Y_i}{1-a(1-\mathrm{\Gamma })},\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill X_i^k(a)& =X_i^k+a\left(\frac{S_B{V}_i{Y}_i}{4}\right)Y_i(a).\hfill \end{array}$$(26)

At this point, the flux in ring i for a given energy group is calculated using

$$\begin{array}{c}\hfill {\varphi }_i=\sum _{k=1}^N{Q}_k{X}_i^k(a)+J_{\mathrm{ext}}{Y}_i(a),\end{array}$$(27)

where Q_k is the source in ring k and J_ext is the external current on the boundary of the problem [25]. The source term contains the sum of the scattering sources, fission sources, and any possible external sources that may exist in the problem:

$$\begin{array}{c}\hfill Q_{i,g}=\sum _{g^{\prime }\ne g}{\mathrm{\Sigma }}_{str,i,g^{\prime }\to g}{\varphi }_{i,g^{\prime }}+\frac{{\chi }_{i,g}}{k_{\mathrm{eff}}}\sum _{g^{\prime }}\nu {\mathrm{\Sigma }}_{f,i,g^{\prime }}{\varphi }_{i,g^{\prime }}+Q_{\mathrm{ext},i,g}.\end{array}$$(28)

Here, νΣ_f, i, g is the fission production cross section in ring i and group g, χ_i, g is the fission spectrum in ring i group g, and k_eff is the eigenvalue of the problem (if being run in eigenvalue mode). When run in k-eigenvalue mode, the standard power iteration method is used to converge to the fundamental eigenstate [36].

The cylindrical collision probabilities solver in Scarabée is very general, and can support a user provided boundary albedo, external current and sources, and any number of annular regions with any number of energy groups. It can be run in fixed-source mode or k-eigenvalue mode. The most common use case, however, is to run single pin-cell calculations with an albedo of unity and no external current in k-eigenvalue mode.

3.3. Discrete ordinates solver for slab geometry

Scarabée has a specialized discrete ordinates (S_n) solver for 1D slab geometries with a reflective boundary condition on the −x side and a vacuum boundary condition on the +x side. Currently, it can only solve k-eigenvalue problems, and is primarily used to generate few-group nodal diffusion parameters for the baffle / reflector region using a scheme very similar to the Baff-Refl method employed in Framatome’s SCIENCE chain, and well outlined by Machach [37, 38]. The solver can be run with isotropic scattering and transport corrected cross sections, or with explicit modeling of anisotropic scattering.

Both the isotropic and anisotropic sweep algorithms make use of a Gauss-Legendre quadrature. Since there is always a vacuum boundary on the right hand side and a reflective boundary on the left, sweeps start on the right hand boundary with all angles μ <  0 moving into the slab, as the entry angular flux is always known to be zero. Subsequently, all angles μ >  0 are considered, where the sweep starts on the left hand boundary using the known reflected flux from the complementary μ <  0 sweep. We will use the notation that ϕ_{i, n, in/out} is the angular flux along direction μ_n that is entering/leaving region i ³. This solver uses the diamond-difference approximation for the average value of the angular flux in region i:

$$\begin{array}{c}\hfill {\varphi }_{i,n}=\frac{{\varphi }_{i,n,\mathrm{in}}+{\varphi }_{i,n,\mathrm{out}}}{2}·\end{array}$$(29)

Using equation (29), it is possible to derive a relation for the outgoing angular flux in the region from the neutron balance equation, resulting in

$$\begin{array}{c}\hfill {\varphi }_{i,n,\mathrm{out}}=\frac{2\mathrm{\Delta }{x}_i{Q}_{i,n}+(2|{\mu }_n|-\mathrm{\Delta }{x}_i{\mathrm{\Sigma }}_{tr,i}){\varphi }_{i,n,\mathrm{in}}}{\mathrm{\Delta }x{\mathrm{\Sigma }}_{tr,i}+2|{\mu }_n|},\end{array}$$(30)

where Δx_i is the width of region i and Q_i, n is the source term for neutrons in region i that enter direction μ_n [24]. In this presentation, we have neglected the energy dependence, but in general, the flux, source, and cross sections have an implicit index for the energy group g. The source term is determined from the flux obtained in a previous iteration, and in the assumption of isotropic scattering with transport corrected data, is calculated as

$$\begin{array}{c}\hfill Q_{i,g,n}=\frac{1}{2}\sum _{g^{\prime }}[{\mathrm{\Sigma }}_{str,g^{\prime }\to g}+\frac{{\chi }_{i,g}}{k_{\mathrm{eff}}}\nu {\mathrm{\Sigma }}_{f,i,g^{\prime }}]{\varphi }_{i,g^{\prime }}.\end{array}$$(31)

The scalar flux in region i is accumulated during a transport sweep as

$$\begin{array}{c}\hfill {\varphi }_i=\sum _n{w}_n\frac{{\varphi }_{i,n,\mathrm{in}}+{\varphi }_{i,n,\mathrm{out}}}{2},\end{array}$$(32)

where w_n is the weight associated with μ_n from the Gauss-Legendre quadrature. As this solver is primarily used for generating nodal diffusion parameters for reflectors, the net current on the boundaries of all regions is also scored for use in generating discontinuity factors [37]. The current at boundary b of region i is calculated as

$$\begin{array}{c}\hfill J_b=\sum _n{w}_n{\mu }_n{\varphi }_{b,n},\end{array}$$(33)

where ϕ_b, n is either ϕ_i, n, in or ϕ_{i, n, out}, depending on the direction of μ_n and wether the boundary b is on the positive or negative side of the region i. The standard power iteration method is used to converge onto the fundamental eigenstate and obtain k_eff [36]. Unfortunately, using isotropic scattering with transport corrected data, the discrete ordinates solver can fail to converge at times. This is caused by negative values in some groups for Σ_tr, g and Σ_{str, g → g}, leading to negative flux values and negative source values, becoming more problematic with more energy groups. To permit convergence, the stabilization method of Gunow is employed [39]. Since implementing this technique, no transport corrected simulations have failed to converge.

When explicitly modeling anisotropic scattering, a few modification must be made. Firstly, the appearances of the transport cross section, Σ_tr, i, in equation (30) must be replaced with the total cross section, Σ_t, i. Second, instead of only accumulating the scalar flux during sweeps, we must accumulate all required moments of the angular flux. If P3 anisotropic scattering data is being used, then one must tally the angular flux in each region up to the l = 3 Legendre moment. This can be accomplished by replacing equation (29) with

$$\begin{array}{c}\hfill {\varphi }_{i,l}=\sum _n{w}_n{P}_l({\mu }_n){\varphi }_{i,n},\end{array}$$(34)

P_l(μ) being the lth order Legendre polynomial. The source Q_i, n is now computed from from the Legendre moments of the source:

$$\begin{array}{c}\hfill Q_{i,n}=\sum _n{P}_l({\mu }_n)Q_{i,n,l},\end{array}$$(35)

where the lth Legendre moment is computed as

$$\begin{array}{cc}\hfill Q_{i,n,l,g}& =\frac{2l+1}{2}\sum _{g^{\prime }}{\mathrm{\Sigma }}_{s,i,g^{\prime }\to g,l}{\varphi }_{i,l,g^{\prime }}\hfill \\ \hfill & +{\delta }_{l,0}\frac{{\chi }_{i,g}}{2k_{\mathrm{eff}}}\sum _{g^{\prime }}\nu {\mathrm{\Sigma }}_{f,i,g^{\prime }}{\varphi }_{i,l,g^{\prime }},\hfill \end{array}$$(36)

which assumes that fission is perfectly isotropic and therefore only contributes to the l = 0 source moment. As the transport correction is not employed when modeling anisotropic scattering, there is no need to apply Gunow’s stabilization technique.

3.4. Method of characteristics solver

Scarabée has a transport solver for 2D geometries with rectilinear outer boundaries based on the method of characteristics (MOC) with the flat source region approximation [29]. While some MOC codes are quite generalized and allow the user to define almost any geometry [40], Scarabée only provides the user with a few distinct primitive sub-geometries or “cells” which are applicable to LWRs [29]. These include cell types for PWR pin cells, BWR pin cells in the corner of a channel box, and empty homogeneous cells. This somewhat limits the systems which can be accurately modeled, but is more than sufficient for LWR analysis and greatly reduces the user effort required to define a given geometry. The solver can be applied to both fixed-source and k-eigenvalue simulations with isotropic scattering and the transport correction or with explicit anisotropic scattering.

The user must provide the number of azimuthal directions that will be examined within the x-y plane of the geometry, the maximum track width, and the polar quadrature used for the z-direction. Based on this information, the exact azimuthal angles and track widths are slightly modified such that cyclic tracks can be traced across the geometry [29]. This means that at each point where a track exits the problem geometry, there will always be a corresponding track at that same position which has a complementary angle. For a given azimuthal angle, one can compute the starting positions of all tracks associated with the angle in advance. From this starting position, the “segment” information is generated for each track, which includes the index of the Flat Source Region (FSR) being traversed in addition to the distance traveled by the track in that FSR. These tracks are only generated for azimuthal angles φ ∈ [0, π], as following them in the “backwards” direction then gives you the corresponding φ ∈ [π, 2π] angles. Each azimuthal angle φ_a(t) (where t is the track index) has a corresponding weight w_a(t). A product quadrature is used to accomplish integration in the z-direction, perpendicular to the tracks which are generated in the x-y plane. This product quadrature provides a set of polar angles (θ_p) and weights (w_p) which are used to project the track lengths into 3D space, as the length of a segment in the x-y plane (l_s) projected along the polar angle θ_p becomes l_ssin(θ_p). The width of track t is denoted as δ_t.

At the beginning of the simulation, a uniform scalar flux is assumed in all energy groups of all FSRs, which is used for computing the initial source terms. The angular flux is only stored for both end points of a track, which then has the angular flux for that track’s azimuthal angle, along all polar angles, in all energy groups. During a “sweep”, all tracks are followed along the forwards and backwards directions, computing the change in the angular flux along each segment. Once all segments have been traversed, one has obtained the angular fluxes in all groups for the corresponding cyclically connected re-entry track at the opposing end of the geometry. One difference between the MOC solver in Scarabée and other codes is that Scarabée uses total source iterations (similar to OpenMOC [40]), where the scattering source and fission source (including the eigenvalue) are updated before each sweep. Many other codes keep the fission source constant, perform the necessary number of sweeps to converge the scattering source, and only then update the fission source and eigenvalue. Testing seemed to indicate that while this leads to a fewer number of power iterations (where the fission source and eigenvalue are updated), a shorter wall-clock compute time can be achieved with total source iteration.

In this derivation, we shall denote each FSR with the index i, each track with the index t, each segment along a track with the index s, each unique azimuthal direction with the index a (dependent on t), and each unique polar angle with the index p. Under the assumption of isotropic scattering with the transport correction, if ϕ_i, s, p(0) is the angular flux at the beginning of segment s on track t in FSR i, traveling in the direction specified by azimuthal angle a and polar angle p, then the angular flux leaving that segment after traveling distance l_s, p is ϕ_i, s, p(l_s, p), and can be shown to be computed as [29, 40]

$$\begin{array}{c}\hfill {\varphi }_{i,s,p}(l_{s,p})={\varphi }_{i,s,p}(0)e^{-{\mathrm{\Sigma }}_{tr,i}{l}_{s,p}}+\frac{Q_i}{{\mathrm{\Sigma }}_{tr,i}}(1-e^{-{\mathrm{\Sigma }}_{tr,i}{l}_{s,p}}).\end{array}$$(37)

Therefore, the change in the angular flux along the segment is

$$\begin{array}{c}\hfill \mathrm{\Delta }{\varphi }_{i,s,p}=({\varphi }_{i,s,p}(0)-\frac{Q_i}{{\mathrm{\Sigma }}_{tr,i}})(1-e^{-{\mathrm{\Sigma }}_{tr,i}{l}_{s,p}}).\end{array}$$(38)

These changes in the angular flux can be integrated over all tracks and angles, leading to the following form for the scalar flux in FSR i:

$$\begin{array}{c}\hfill {\varphi }_i=\frac{4\pi }{{\mathrm{\Sigma }}_{tr,i}}[Q_i+\frac{1}{A_i}\sum _{t\in i}\sum _{s\in i}\sum _p{\delta }_t{w}_a(t)w_{p}sin({\theta }_p)\mathrm{\Delta }{\varphi }_{i,s,p}],\end{array}$$(39)

where A_i is the 2D area of FSR i. As before, the group index on all flux, source, and cross section terms has been neglected for conciseness. The source term for a sweep can be computed from the flux obtained by the previous iteration:

$$\begin{array}{c}\hfill Q_{i,g}=\frac{1}{4\pi }\sum _{g^{\prime }}[{\mathrm{\Sigma }}_{str,g^{\prime }\to g}+\frac{{\chi }_{i,g}}{k_{\mathrm{eff}}}\nu {\mathrm{\Sigma }}_{f,i,g^{\prime }}]{\varphi }_{i,g^{\prime }}.\end{array}$$(40)

When using isotropic scattering with the transport correction, Gunow’s stabilization method is employed to avoid convergence problems caused by negative transport cross sections and negative diagonal elements of the transport corrected scattering matrices.

3.4.1. Anisotropic scattering

The explicit anisotropic scattering treatment in the MOC solver uses a spherical harmonics expansion. The main difference between the isotropic and anisotropic treatment is that the scattering source term for each segment crossing each flat source region is now different, and depends on the direction. We can write the source term for segment s crossing FSR i as [41]

$$\begin{array}{cc}\hfill Q_{i,s,p,g}& =\sum _{g^{\prime }}[\sum _{l=0}^N\left({\mathrm{\Sigma }}_{s,l,g^{\prime }\to g}\sum _{j=-l}^l{\varphi }_{i,l,j,g^{\prime }}{Y}_{l,j}({\widehat{\mathbf{\Omega }}}_{s,p})\right)\hfill \\ \hfill & +\frac{1}{\sqrt{4\pi }}\frac{{\chi }_{i,g}}{k_{\mathrm{eff}}}\nu {\mathrm{\Sigma }}_{f,i,g^{\prime }}{\varphi }_{i,0,0,g}].\hfill \end{array}$$(41)

Here, the source term Q_{i, s, p, g} is dependent on the track’s azimuthal angle and polar angle, in addition to the angular flux moments, ϕ_{i, l, j, g}, which describe the shape of the angular flux on the unit sphere, and $Y_{l,j}({\widehat{\mathbf{\Omega }}}_{s,p})$ which is the real spherical harmonic function [24]. Therefore, one must calculate the angular flux moments during transport sweeps for all js and ls, where j varies from −l to l for each l = 0, 1, …, N if PN anisotropic scattering data is being used. One will notice that equation (41) assumed the fission source is perfectly isotropic, and only the l = 0, j = 0 flux moment contributes to fission source component.

The expressions for the angular flux leaving a segment, given by equation (37), and the change in the angular flux along the segment, given by equation (38), can be used in the anisotropic scattering treatment with only two minor changes. First, Q_i will be replaced by Q_i, s, p, which incorporates the anisotropic scattering of neutrons from one direction to another. Secondly, the transport corrected cross section will be replaced by the total cross section, Σ_t. However, since the higher spherical harmonic moments of the angular flux are required when calculating the source term, it is no longer sufficient to only score the scalar flux. We must also score the ϕ_{i, l, j, g} angular flux moments, which can be accomplished with the following relation [41]:

$$\begin{array}{cc}& {\varphi }_{i,l,j,g}=\frac{4\pi }{{\mathrm{\Sigma }}_{t,i,g}{A}_i}[\sum _{t\in i}\sum _{s\in i}\sum _p{\delta }_t{w}_a(t)w_p\hfill \\ \hfill & (sin({\theta }_p)\mathrm{\Delta }{\varphi }_{i,s,p,l,j,g}+\sum _{l=0}^N\sum _{j=-l}^l{l}_{s,p}{Q}_{i,s,p,g})Y_{l,j}({\widehat{\mathbf{\Omega }}}_{s,p})].\hfill \end{array}$$(42)

In the explicit anisotropic scattering treatment, the transport correction is not employed. Therefore, there is no need to apply Gunow’s stabilization technique to fix the flux moments. At this point, it is likely clear to the reader that modeling anisotropic scattering with a PN approximation requires (N + 1)² as much memory as the isotropic scattering simulation. It will also require more time to perform the computation, as equation (42) must now be scored for all angular flux moments. Therefore, anisotropic scattering should only be used when required to accurately model the problem of interest.

3.4.2. Coarse Mesh Finite-Difference (CMFD) acceleration

The Scarabée MOC solver includes optional non-linear diffusion acceleration in the form of a Coarse Mesh Finite-Difference (CMFD) solver [29, 42]. CMFD is a common method of accelerating transport problems, particularly in large geometries, and has been implemented in several other codes [4, 40, 42]. In CMFD, a coarse Cartesian mesh is superimposed over the fine FSR mesh used in the MOC solver⁴, with each coarse CMFD cell typically containing multiple FSRs. In order to remain consistent with the transport solution from the MOC solver, a non-linear diffusion correction is added to the surface diffusion coefficients used in the finite-difference diffusion solver, based on the currents tallied during a MOC sweep. During a MOC sweep, the net current is tallied on each surface, S, of each CMFD cell with

$$\begin{array}{c}\hfill J_{G,S}=\frac{1}{l_S}\sum _{g\in G}\sum _{t\in {\mathcal{T}}_S}\sum _p{\omega }_a(t){\omega }_{p}sin{\theta }_p{\delta }_t{\psi }_{t,g,p}\mathrm{sign}(cos{\phi }_a(t)),\end{array}$$(43)

where 𝒯_S is the set of all tracks crossing surface S and l_S is the length of the surface. When using CMFD, it is possible to specify a condensed group structure that is used to solve the diffusion problem and update the MOC scalar flux. In equation (43), g is the MOC energy group corresponding to the condensed CMFD energy group G. Performing CMFD on a condensed energy grid is generally beneficial for run times.

In the CMFD formalism, the net current crossing the surface S is defined as

$$\begin{array}{c}\hfill J_{G,S}=-{\widehat{D}}_S({\varphi }_{G,S+}-{\varphi }_{G,S-})-{\tilde{D}}_S({\varphi }_{G,S+}+{\varphi }_{G,S-}),\end{array}$$(44)

where S+ represents the CMFD cell on the positive side of the surface, and S− represents the CMFD cell on the negative surface. ${\widehat{D}}_{G,S}$ is the standard diffusion coefficient for the surface, while ${\tilde{D}}_{G,S}$ is the non-linear correction which is computed to ensure that the current tallied in the MOC solver is preserved in the CMFD formalism [43]. If the MOC problem has a reflective boundary at surface S, it is assumed that J_G, S = 0 and the tallied current is not utilized; a similar assumption is made for the case of a vacuum boundary conditions⁵ [24]. The choice of how diffusion coefficients are calculated can have important consequences for the stability and convergence of the CMFD solver, particularly for cells with a large optical thickness. Therefore, Scarabée includes two optional methods to deal with this very problem, one being Larsen’s effective diffusion coefficient [44] and the other being to adjust the diffusion coefficient to an optimized value based on the cell optical thickness [45]. By default, the latter method is used in Scarabée, but this can be disabled or replaced with Larsen’s method should the user desire.

Once the surface diffusion coefficients are calculated, Scarabée solves the finite-difference diffusion problem in the typical manner for either k-eigenvalue or fixed-source problems. CMFD is applied at every MOC iteration, and the solution is used to update the MOC scalar flux in each FSR before starting the next MOC iteration:

$$\begin{array}{c}\hfill {\varphi }_{g,i}^{n+1}={\varphi }_{g,i}^n\frac{{\varphi }_{G,I}^{\mathrm{CMFD}}}{{\varphi }_{G,I}^{\mathrm{MOC}}}\forall (g,i)\in (G,I),\end{array}$$(45)

where the FSR i is contained within the CMFD cell I [43]. The CMFD method greatly accelerates the convergence of the flux in the MOC problem, particularly in problems with large reflectors. An example of this acceleration is provided in Section 4.2.

3.5. Diffusion solvers

While it is ideal to perform neutronics calculations using transport theory whenever possible, this is still a cumbersome endeavor for 3D full core problems on laptops or small workstations. To facilitate the analysis of full core problems, Scarabée provides two different solvers based on diffusion theory for Cartesian geometries.

3.5.1. Finite-difference solver

The first of these is a finite-difference solver based on the method of cell-centered finite-differences [24]. It can solve 1D, 2D, or 3D problems with generalized albedo boundary conditions, and supports the use of assembly discontinuity factors [46]. The theory behind such solvers has been covered in many different texts, though we encourage the reader to examine the excellent explanation provided by Hébert [24].

3.5.2. Nodal expansion method solver

For large problems, the finite-difference methodology can still be somewhat slow, particularly if one must perform many thousands of calculations or examine reactor transients. As an alternative, users can instead use the nodal diffusion solver based on the Nodal Expansion Method (NEM) [47]. In nodal diffusion methods, the problem domain is discretized into relatively large (sometimes assembly size) “nodes” which have homogeneous material properties. Within node i, we want to describe the neutron flux as a function of position, ϕ_i(x, y, z), which is described by the diffusion equation which we can write as

$$\begin{array}{c}\hfill \mathbf{\nabla }·{\mathit{J}}_i(x,y,z)+{\mathrm{\Sigma }}_{r,i}{\varphi }_i(x,y,z)=Q_i(x,y,z),\end{array}$$(46)

where J _i(x, y, z)= − D_i ∇ ϕ_i(x, y, z) is the neutron current, Σ_r, i is the removal cross section and Q_i(x, y, z) is the position dependent source within the node. To simplify the problem, only one axis (x, y, or z) is solved for at a time. If solving for the flux along the x direction, equation (46) is integrated along the y and z axes, applying the operator

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}(·)\mathrm{d}z\mathrm{d}y\end{array}$$(47)

to every term, where Δy is the width of the node along the y direction and similarly for Δz, with the node being centered at the origin. Following this approach, the average value of the flux in the node as a function of the x position is

$$\begin{array}{c}\hfill {\varphi }_{x,i}(x)\equiv \frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}{\varphi }_i(x,y,z)\mathrm{d}z\mathrm{d}y,\end{array}$$(48)

and the average value of the source in the node as a function of the x position is

$$\begin{array}{c}\hfill Q_{x,i}(x)\equiv \frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}{Q}_i(x,y,z)\mathrm{d}z\mathrm{d}y.\end{array}$$(49)

For each of the three terms from the divergence of the current in equation (46) we have a different resulting expression. For the current along the x direction, we obtain

$$\begin{array}{cc}& \frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}\frac{\partial }{\partial x}{J}_{x,i}(x,y,z)\mathrm{d}z\mathrm{d}y\hfill \\ \hfill & =\frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}-D_i\frac{{\partial }^2{\varphi }_i(x,y,z)}{\partial x^2}\mathrm{d}z\mathrm{d}y\hfill \\ \hfill & =-D_i\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}{\varphi }_{x,i}(x).\hfill \end{array}$$(50)

For the currents in the two transverse directions, one obtains

$$\begin{array}{cc}& \frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}\frac{\partial }{\partial y}{J}_{y,i}(x,y,z)\mathrm{d}z\mathrm{d}y\hfill \\ \hfill & =\frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}[J_{y,i}(x,\mathrm{\Delta }y/2,z)-J_{y,i}(x,-\mathrm{\Delta }y/2,z)]\mathrm{d}z\hfill \\ \hfill & \equiv \frac{1}{\mathrm{\Delta }y}{L}_{y,i}(x)\hfill \end{array}$$(51)

for the transverse leakage in the y direction and

$$\begin{array}{cc}& \frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}{\int }_{-\mathrm{\Delta }z/2}^{\mathrm{\Delta }z/2}\frac{\partial }{\partial z}{J}_{z,i}(x,y,z)\mathrm{d}z\mathrm{d}y\hfill \\ \hfill & =\frac{1}{\mathrm{\Delta }y\mathrm{\Delta }z}{\int }_{-\mathrm{\Delta }y/2}^{\mathrm{\Delta }y/2}[J_{z,i}(x,y,\mathrm{\Delta }z/2)-J_{z,i}(x,y,-\mathrm{\Delta }z/2)]\mathrm{d}y\hfill \\ \hfill & \equiv \frac{1}{\mathrm{\Delta }z}{L}_{z,i}(x)\hfill \end{array}$$(52)

for the transverse leakage in the z direction. Combining these definitions we arrive at the second order ordinary differential equation for the average value of the flux in the x direction:

$$\begin{array}{cc}& -D_i\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}{\varphi }_{x,i}(x)+{\mathrm{\Sigma }}_{r,i}{\varphi }_{x,i}(x)\hfill \\ \hfill & =Q_{x,i}(x)-\frac{1}{\mathrm{\Delta }y}{L}_{y,i}(x)-\frac{1}{\mathrm{\Delta }z}{L}_{z,i}(x).\hfill \end{array}$$(53)

Naturally, there are corresponding equations for ϕ_y, i(y) and ϕ_z, i(z) that must be solved.

In the NEM approximation, we solve equation (53) by approximating the flux as a 4th order polynomial which takes the form

$$\begin{array}{c}\hfill {\varphi }_{x,i}(x)={\overline{\varphi }}_i+\sum _{n=1}^4{a}_{x,i,n}{f}_n(x),\end{array}$$(54)

where the basis functions are

$$\begin{array}{cc}\hfill f_1(x)& =\frac{x}{\mathrm{\Delta }x}\hfill \\ \hfill f_2(x)& =3{\left(\frac{x}{\mathrm{\Delta }x}\right)}^2-\frac{1}{4}\hfill \\ \hfill f_3(x)& =\frac{x}{\mathrm{\Delta }x}(\frac{x}{\mathrm{\Delta }x}-\frac{1}{2})(\frac{x}{\mathrm{\Delta }x}+\frac{1}{2})\hfill \\ \hfill f_4(x)& =({\left(\frac{x}{\mathrm{\Delta }x}\right)}^2-\frac{1}{20})(\frac{x}{\mathrm{\Delta }x}-\frac{1}{2})(\frac{x}{\mathrm{\Delta }x}+\frac{1}{2}).\hfill \end{array}$$(55)

As outlined by Lawrence, using a 4th order polynomial for the flux requires defining two supplementary moment equations, taking the inner product of equation (53) with f₁(x) and f₂(x) over the node [47]. Additionally, one must define the functional forms of the transverse leakage terms L_y, i(x) and L_z, i(x). This is approximated as a quadratic polynomial such that the integral of L_y, i(x) over the domain of the node in question (i) and the two adjacent nodes (i − 1 and i + 1) separately recovers the average transverse leakage of the node in question. While this approximation has been shown to yield excellent results, it is not evident how to construct such a quadratic polynomial when the node i lies next to a boundary condition, and the node i ± 1 does not exist. For these nodes, Scarabée assumes that L_y, i(x) is a constant equal to the value of the transverse leakage in node i, as Zerkle indicated that this was a relatively accurate and simple approach for these particular nodes [48]. To obtain the solution to a problem, the values of ${\overline{\varphi }}_i$ and a_x, i, n are iteratively solved for across all nodes and energy groups, making sure that the entry current of one node equals the exit current of the adjacent node, taking the derivative of equation (54). In Scarabée, this is accomplished by using the interface current response matrix formalism which has been previously outline by Lawrence [47].

Once the nodal diffusion problem has been solved, one only has values of the node average flux, and the partial currents for each surface of each node. Therefore, we are not directly given any information regarding the intranodal shape of the flux, and therefore the intranodal shape of the power distribution. Without taking the intranodal power distribution into account, one will have very poor estimates for relative pin powers. Therefore, a “reconstruction” technique must be applied where the node average flux and partial currents are used to approximate the intranodal flux shape in each group from which the intranodal power distribution can also be obtained. Many other works have considered this problem, and most make assumptions regarding the number of energy groups being used or the flux behavior in specific groups [49, 50]. While such approximations may have been necessary in the past to obtain reasonable performance, this is no-longer the case. Scarabée applies the Analysis of Variance - High Dimensional Model Representation (ANOVA-HDMR) method to each energy group [51]. Once the reconstruction step has been completed, the user can interrogate the value of the flux or power at any position in the problem domain. All values are normalized such that the average power of all nodes in the problem is unity.

3.6. Lattice calculations

Scarabée also has lattice physics capabilities which allow the user to perform 2D simulations of single PWR fuel assemblies. This can be used to generate few-group diffusion cross sections and discontinuity factors for use in the nodal diffusion solver, and also includes fuel depletion capabilities. Geometries which include normal fuel pins, guide tubes, gadolinium fuel pins, wet burnable absorber rods, pyrex rods, and fuel pin spacer grids with grid sleeves can all be simulated. Currently, there is no modeling of control rods, though this feature is planned for the near future.

Unlike other lattice physics codes, the assembly calculations are only carried out in the group structure of the nuclear data library using the MOC solver; no energy condensation step is performed. This certainly leads to longer computation times, but also more accurate results and eliminates the need to come up with ad-hoc methods for determining the flux spectrum which should be used when condensing cross sections for the grid sleeve and assembly gap. By default, CMFD is used to accelerate flux convergence, helping to reduce run times, particularly when performing depletion calculations. Unique Dancoff corrections are computed for each fuel pin and cladding (including guide tubes) using the neutron current method [29, 34]. These Dancoff corrections are then used in the self-shielding step, along with the geometry and material composition for each fuel pin, to interpolate unique macroscopic cross sections for each fuel pin. When simulating assemblies with gadolinium fuel pins, the user can specify how many rings should be used for spatial self-shielding and depletion, and the Stoker-Weiss method is used to obtain unique cross sections for each fuel ring. While isotropic scattering with the transport correction is used by default, users may choose to explicitly model anisotropic scattering, up to the P3 moment, if the nuclear data library provides the necessary scattering matrices. Single assembly calculations can not directly model the spectral effects caused by neutrons “leaking” between assemblies to achieve criticality, as would occur in an operating reactor. To approximate this effect on the flux spectrum, either a P1, B1, or Fundamental Mode leakage correction can be applied to the MOC results [24, 25, 29]. By default, a P1 leakage model is used, based on the suggestion and reasoning provided by Smith [27]. To generate the few-group diffusion cross sections, standard flux-volume weighting is used across the entire assembly [29]. As suggested by Smith, the diffusion coefficient is directly condensed in energy, instead of first condensing the transport cross section and then converting to the diffusion coefficient [27].

3.6.1. ADF calculation from MOC results

In addition to the few-group cross sections, assembly discontinuity factors are required for the generalized equivalence theory proposed by Smith [52], which greatly improves the results of subsequent full core diffusion simulations. The definition of the discontinuity factor is typically given as

$$\begin{array}{c}\hfill f_{x+}=\frac{{\varphi }_{\mathrm{het}}(x_{+})}{{\varphi }_{\mathrm{hom}}(x_{+})},\end{array}$$(56)

where ϕ_het(x₊) is the average value of the flux along the x+ surface of the assembly from the detailed MOC calculation, and ϕ_hom(x + ) is the average value of the flux along the x+ surface from the homogeneous (i.e. nodal diffusion) solution to the equivalent problem. Since assembly calculations are often performed with reflective boundary conditions, ϕ_hom(x₊) is typically just the average value of the assembly flux, and therefore independent of the surface being considered. To obtain ϕ_het(x₊), the average flux along the FSRs on the boundary of the desired surface is computed using length weighting [53]. A similar process is used when computing the corner discontinuity factors which are required for the flux-reconstruction processing in the nodal diffusion solver.

3.6.2. ADF calculation from CMFD results

When Smith initially developed and proposed the use of ADFs, assembly calculations were often being performed with homogenized pin cell calculations. Smith has since remarked that the naïve approach for calculating the ADF based on equation (56) is incorrect for modern lattice physics methods where the true heterogeneous geometry is modeled, such as in the MOC⁶ [27]. The reasoning for this is that ADFs are meant to account for heterogeneity (meaning different types of cells) within a given assembly. It does not make sense to consider the heterogeneity of the assembly at the sub pin cell level, as the downstream nodal diffusion simulation ignores such fine detail, and is looking for the smooth homogeneous solution for the reactor flux. With this in mind, if one were to simulate an assembly with only one type of fuel pin, no guide tubes, and no gap on the exterior of the assembly, the ADFs (and CDFs) in each group should be exactly unity, as there would be no heterogeneity within the assembly. Clearly, computing the ADFs (and CDFs) using the method outlined in Section 3.6.1 will not result in unity, as the flux in the moderator will be used to compute ϕ_het along the edge of an assembly, and this is not representative of the behavior of the homogenized pin cell flux on the edge of the assembly. To obtain valid ADFs from MOC lattice calculations, Smith instead proposes tallying the currents for the interface of each pin cell along the exterior of the assembly along with homogenizing diffusion cross sections for those cells. Then, by performing a fixed-source nodal calculation on each pin cell with the reference currents, one can recover the homogenized surface flux within that pin cell which can be integrated along the length of the assembly to compute ϕ_het [27].

In Scarabée, a method to compute ADFs based on this concept has been implemented, leveraging the currents which are tallied when CMFD is used to accelerate the calculation. The CMFD mesh in Scarabée corresponds to the pin cells of the assembly, with edge pin cells also having the extra flat source regions which may be present for the moderator gap between assemblies and for the grid sleeve that may be present. The CMFD current is generally tallied in an intermediate group structure, so the currents for the four sides of a pin cell are first condensed down to the few-group structure used for the nodal diffusion parameters. Subsequently, a set of few-group diffusion cross sections are homogenized and condensed for the pin cell. Few-group average fluxes within the pin cell are also computed. All these parameters can then be used to perform a 1D fixed-source calculation for the pin cell. We shall assume that the average flux in group G for the pin cell is ${\overline{\varphi }}_G$, and the net currents on the +x and −x sides of the cell are J_G ^+x and J_G ^−x , respectively. Machach has shown that the four coefficients, a_x, n, can be solved for using the following system of equations⁷ [37, 38]:

$$\begin{array}{c}\hfill -\frac{D_G}{2\mathrm{\Delta }{x}^2}{a}_{x,3}+\frac{{\mathrm{\Sigma }}_{r,G}}{12}(a_{x,1}-\frac{a_{x,3}}{10})\\ \hfill -\sum _{H\ne G}(\frac{{\mathrm{\Sigma }}_{s,H\to G}}{12}+\frac{{\chi }_G}{k_{\mathrm{eff}}}\frac{\nu {\mathrm{\Sigma }}_{f,H}}{12})(a_{x,1}-\frac{a_{x,3}}{10})& =0\hfill \end{array}$$(57)

$$\begin{array}{c}\hfill -\frac{D_G}{5\mathrm{\Delta }{x}^2}{a}_{x,4}+\frac{{\mathrm{\Sigma }}_{r,G}}{20}(a_{x,2}-\frac{a_{x,4}}{35})\\ \hfill -\sum _{H\ne G}(\frac{{\mathrm{\Sigma }}_{s,H\to G}}{20}+\frac{{\chi }_G}{k_{\mathrm{eff}}}\frac{\nu {\mathrm{\Sigma }}_{f,H}}{20})(a_{x,2}-\frac{a_{x,4}}{35})& =0\hfill \end{array}$$(58)

$$\begin{array}{cc}\hfill \frac{D_G}{\mathrm{\Delta }x}(-a_{x,1}-3a_{x,2}-\frac{a_{x,3}}{2}-\frac{a_{x,4}}{5})& =J_G^{+x}\hfill \end{array}$$(59)

$$\begin{array}{cc}\hfill \frac{D_G}{\mathrm{\Delta }x}(-a_{x,1}+3a_{x,2}-\frac{a_{x,3}}{2}+\frac{a_{x,4}}{5})& =J_G^{-x}\hfill \end{array}$$(60)

where Δx is the width of the cell in the x direction. With all a_x, n determined for each group, one can use equation (54) to demonstrate that the average value of the homogenized pin cell flux at the +x boundary of cell c in group G is

$$\begin{array}{c}\hfill {\varphi }_{c,G,+x}={\overline{\varphi }}_G+\frac{a_{x,1,G}}{2}+\frac{a_{x,2,G}}{2}·\end{array}$$(61)

Similarly, for the average value of the homogenized pin cell flux at the −x boundary we find that

$$\begin{array}{c}\hfill {\varphi }_{c,G,-x}={\overline{\varphi }}_G-\frac{a_{x,1,G}}{2}+\frac{a_{x,2,G}}{2}·\end{array}$$(62)

Using these average homogenized pin cell fluxes on the +x boundary, we can compute ϕ_het, G(+x) as

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{het},G}(+x)=\frac{{\displaystyle \sum _{c\in +x}{l}_{c,+x}{\varphi }_{c,G,+x}}}{{\displaystyle \sum _{c\in +x}{l}_{c,+x}}},\end{array}$$(63)

where the sum is carried out for all cells c which touch the +x boundary of the assembly and l_c, +x is the length of cell c along the +x assembly boundary.

The homogeneous flux used in the denominator of equation (56) is calculated in a similar manner, but on the scale of the entire assembly. The average net current on each face of the assembly is condensed to the few-group structure from the CMFD results, and few-group diffusion cross sections are also generated. Then, equation (57) through (60) are solved for the assembly, and the average homogenized assembly surface flux can be computed using equation (61) or equation (62). CDFs are computed in a very similar procedure, where the heterogeneous corner flux is computed from a 2D nodal calculation for just the corner pin cell, and the homogeneous corner flux is computed from a 2D nodal calculation for the entire assembly. The value of the corner flux from the 2D nodal calculation is based on the ANOVA-HDRM reconstruction without considering any x-y cross terms [51]. Using this procedure on perfectly homogeneous assemblies with no gap, ADFs and CDFs well within 0.1% of unity can be obtained. When CMFD is used to accelerate assembly problems in Scarabée, this improved technique is used to generate ADFs and CDFs.

Fig. 4. Diagram of the simplified 1D reflector model used to generate few-group diffusion cross sections and discontinuity factors for reflector nodes. This scheme is inspired by the Baff-Ref method [37, 38].

3.6.3. Reflector calculations

To generate few-group diffusion cross sections with discontinuity factors for reflector regions, an approach similar to the improved Baff-Ref technique outlined by Machach is used [37, 38]. A model fuel assembly from the periphery of the core is chosen, and then simulated using a standard k-eigenvalue calculation. This assembly is homogenized to generate cross sections in the group structure of the library. Subsequently, a simplified 1D geometry representing a reflector is created. It consists of a reflective boundary condition on the −x boundary, followed by one assembly width of the homogenized assembly material, a moderator region the width of the gap between the assemblies and the baffle, a baffle region, and then the remaining moderator out to the vacuum boundary condition on the +x boundary. A depiction of this geometry can be found in Figure 4. This geometry is discretized and then a k-eigenvalue problem is solved using the 1D discrete ordinates solver outlined in Section 3.3. Once this simulation is finished, the resulting scalar flux in the reflector regions is used to flux-volume average the cross sections in the reflector before performing energy condensation down to the few-group structure. These homogenized few-group cross sections, along with the few-group net currents which were tallied at the assembly-reflector interface and the vacuum boundary are used to perform a 1D nodal fixed-source calculation, from which the shape of the homogeneous flux can be obtained. This homogeneous flux, along with the heterogeneous flux from the discrete ordinates solve, are used to determine the discontinuity factor at the assembly-reflector interface using equation (56) and the general approach outlined in Section 3.6.2. The discontinuity factor on the reflector side is then normalized by the discontinuity factor on the fuel assembly side; in the full core diffusion model, the reflector discontinuity factor is then multiplied by the adjacent fuel assembly discontinuity factor, as suggested by Smith [27]. It has been observed that the homogeneous flux from the fixed source nodal calculation can be negative at the vacuum boundary, which is of course nonphysical. When this is observed, the previously computed diffusion coefficients from all groups are iteratively increased by a small fraction until the homogeneous flux is positive at the vacuum boundary. Arbitrarily changing the diffusion coefficients in this manner is licit, as the discontinuity factors ensure that the required leakage is preserved [52]. These resulting diffusion coefficients and discontinuity factors are then edited for use in the full core simulation.

3.6.4. Depletion

One of the primary uses of lattice physics codes is to perform depletion calculations so that reactor designers can predict the life time of a fuel cycle and the fuel composition at various time points. This is accomplished by performing a series of k-eigenvalue simulations, where the resulting reaction rates and user provided power normalization are used to perform time integration of the Bateman equation, predicting the isotopic composition of materials at future time points. For a given material region, the Bateman equation which models the time rate of change in isotope i is

$$\begin{array}{cc}& \frac{\mathrm{d}{N}_i(t)}{\mathrm{d}t}\hfill \\ \hfill & =\sum _j{N}_j(t)[\sum _r{\int }_0^{\infty }{\nu }_{r,j\to i}(E){\sigma }_{j,r}(E)\varphi (E,t)\mathrm{d}E+{\lambda }_{j\to i}]\hfill \\ \hfill & -N_i(t)[{\int }_0^{\infty }{\sigma }_j(E)\varphi (E,t)\mathrm{d}E+\sum _j{\lambda }_{i\to j}].\hfill \end{array}$$(64)

In this form, N_i(t) is the concentration of isotope i at time t, σ_i is the total cross section of isotope i, σ_i, r is isotope i’s cross section for reaction r, ν_{r, i → j} is the yield for reaction r on isotope i producing isotope j, λ_i → j is the decay constant for isotope i to decay to isotope j, and ϕ(E, t) is the scalar neutron flux. Often, this system of coupled differential equations is written as a matrix equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathit{N}(t)}{\mathrm{d}t}=\widehat{\mathit{A}}(t)\mathit{N}(t),\end{array}$$(65)

where N is a vector which contains the concentration of all isotopes in the material, and $\widehat{\mathit{A}}$ is the “depletion matrix” which represents the rate at which isotope i is transforming into isotope j. While in reality there are thousands of isotopes which might appear in a single material region of a nuclear reactor, many lattice physics codes might only track approximately 200 of these [29].

The main difficulty with performing the integration of equation (65) is that the depletion matrix is time dependent. This is due to the fact that the flux in equation (64) is time dependent, as changes to N will result in changes to ϕ(E). Many time integration schemes have been devised so as to adequately model this time dependence in the burn-up matrix [54]. In Scarabée, a predictor-corrector scheme known as LE/QI (Linear Extrapolation on the predictor and Quadratic Interpolation on the corrector) is used with a 4th order commutator free Magnus expansion integrator [54]. To initialize the depletion calculation, a single CE/LI (Constant Extrapolation predictor and Linear Interpolation corrector) time step is performed. The choice of LE/QI was primarily driven by prior observations from Isotalo et al. that other lower order depletion schemes can suffer from problems related to the power normalization not remaining constant over the depletion step [55].

In practice, the depletion matrix is quite sparse. To reduce memory requirements, a specialized data structure based on the Eigen sparse matrix has been implemented on the C++ side of Scarabée, with corresponding Python bindings. The time integration procedure requires that several matrix exponentials be computed [6, 54]. To accomplish this task, Scarabée employs the 16th order Chebyshev Rational Approximation Method (CRAM) to compute the product of a matrix exponential with a vector [56, 57]. To reduce the amount of nuclear data which needs to be stored in the nuclear data library, currently only (n, fission), (n, γ), (n, α), (n, p), (n, 2n), and (n, 3n) reactions are tracked for depletion. Any arbitrary number of isotopes could be modeled (baring memory constraints), and a user can provide an OpenMC formatted depletion chain when creating a nuclear data library [58]. When running a PWR assembly problem with depletion, each fuel pin is automatically depleted independently. If the pin was discretized into rings, each ring is treated as a separate material and depleted independently. In addition to fuel materials, the poison regions of burnable absorber rods are also depleted. No other materials in the simulation are considered to undergo transmutation or radioactive decay. To build the depletion matrix, $\widehat{\mathit{A}}$, the user must specify the known linear power rating of the entire assembly (without considering symmetry) in units of kW cm⁻¹. Currently, Scarabée only computes the assembly power based on the average energy released per fission for the fissile isotopes. Heating due to other reactions or the effects of photon transport on the heating distribution are not currently considered.

4. Example usage and results

4.1. 1D annular pin cell

This solver is primarily used for performing independent pin cell flux spectrum calculations in the fine-group structure of the nuclear data library. As an example, a single pin-cell calculation has been performed with both Scarabée and OpenMC for the 2.4% enriched pin of the BEAVRS benchmark, with all materials at a temperature of 575 K [59]. For the Scarabée calculation, Dancoff factors were calculated with the 1D CP solver as opposed to the MOC solver, as the OpenMC calculation was also performed in annular geometry with white boundary conditions on the outer surface; additionally, only one annular ring was used per material region (i.e. the fuel pellet was not discretized into several rings) and the convergence criteria was a maximum relative difference in the flux of 10⁻¹⁰. The OpenMC calculation used 10⁵ particles per batch with 2000 active generations and 100 inactive generations. Using ENDF/B-VIII.0, OpenMC estimated k_eff = 1.14534 ± 0.00006, while Scarabée obtained k_eff = 1.14547, a difference of only 13 pcm. Figure 5 presents the flux spectrum in the fuel region calculated by both OpenMC and Scarabée, along with the relative difference between the two. With the exception of just a few energy groups, the flux spectrum predicted by Scarabée is largely within 5% relative difference of the Monte Carlo results. The regions where the relative difference exceedes 5% are areas where the flux has very small values (such as in the 21 eV resonance of ²³⁸U or above 10 MeV), so the relative difference becomes a poor metric.

Fig. 5. Comparison of the flux spectrum from Scarabée and OpenMC for a 1D annular pin cell problem using Scarabée’s 1D annular collision probabilities solver.

4.2. Method of characteristics solver

To demonstrate the MOC solver, we have chosen to examine the C5G7 benchmark, as many models exist for various transport codes, thus facilitating comparisons [60]. In 2D, this benchmark has four 17 × 17 PWR fuel assemblies in a 1/4 mini-core arrangement, and is defined with 7 energy groups. Figure 6 provides an example of the interactive geometry plotter provided in Scarabée for visualizing MOC geometries. Users can choose wether to color by FSR, or material, and can mark boundaries and change colors if desired. This tools is very useful for verifying that a core geometry as been correctly defined. The plotter renders the image using the same geometry interface used in tracing the characteristics, so it is guaranteed to be representative of the geometry “seen” in the simulation.

Fig. 6. A portion of the C5G7 geometry in the graphical plotter tool for examining MOC simulation geometries.

When simulating C5G7, 64 azimuthal angles over [0, 2π] were used, along with a maximum track spacing of 0.05 cm. The Yamamoto-Tabuchi 6 point polar quadrature (3 points over [0, π/2]) was employed to form the product quadrature, which is an excellent choice for problems with isotropic scattering [61]. In addition to these basic simulation settings, CMFD was also employed to accelerate the convergence of the problem. CMFD tiles were chosen to be the same size as a pin cell (1.26 cm side length) over both the assemblies and the reflector. Only 3 energy groups were used for the CMFD, with the two fast groups, the three epithermal groups, and the two thermal groups all being condensed separately. As a point of reference, a Monte Carlo solution obtained with OpenMC resulted in an eigenvalue of k_eff = 1.18648 ± 0.00002 using 2 × 10⁴ active generations (500 inactive), each with 10⁵ particles. The Scarabée MOC solver resulted in an eigenvalue of k_eff = 1.18657, and only took approximately 51 seconds on a Lenovo Yoga Pro 9 16IMH9 laptop to converge to a maximum relative difference in the flux of 10⁻⁵, corresponding to 34 iterations. Without CMFD, achieving a similar level of convergence requires 897 iterations and approximately 903 seconds⁸. The resulting thermal flux (group 7) is provided in Figure 7. Figure 8 gives the convergence rate of the problem when CMFD is or is not used for acceleration. From this, we can see that the CMFD implementation in Scarabée is highly effective at accelerating the flux convergence. Additionally, the eigenvalues obtained with and without CMFD are exactly the same, to within floating point precision, assuming the flux was converged to near floating point differences. This indicates that the CMFD implementation is correct, and is correctly treating nuances, such as corner crossings [62].

Fig. 7. Thermal flux from the C5G7 benchmark predicted by Scarabée.

Fig. 8. Convergence of the C5G7 benchmark problem with and without CMFD to accelerate flux convergence.

Fig. 9. Comparison of depletion results from Scarabée with OpenMC.

Fig. 10. Comparison of the pin power distributions obtained by OpenMC and Scarabée for the 2D variant of the BEAVRS benchmark.

4.3. Assembly depletion calculation

As a demonstration of the lattice physics calculation capabilities of Scarabée, we have chosen to demonstrate a depletion simulation of a PWR assembly based on the initial EPR core specification [63]. The chosen assembly has a 3 w/o enrichment, and 20 Gadolinia poison pins which were 8 w/o Gadolinia, and the remaining 92 w/o was UO₂ with an enrichment of 0.25 w/o. This specific assembly was chosen as the depletion of Gadolinia pins is particularly challenging, due to the very strong spatial self-shielding effects and the fast burn-up rate. It therefore serves as an excellent test and verification of the Stoker-Weiss self-shielding implementation, as well as the depletion capabilities of Scarabée. Listing 1 provides the first half of the Python script for Scarabée where the materials are defined. For the reader who is already accustomed to the OpenMC Python interface [9], this input should feel very familiar. Listing 2 contains the second portion of the Python script, where the fuel pins, guide tubes, and assembly geometry are defined. Additional parameters such as the moderator conditions, the linear power rating, and depletion time steps are also specified. While a P1 critical leakage model is used by default in Scarabée, it has been disabled in this case to facilitate comparisons with OpenMC which is not able to incorporate such a leakage model. An equivalent depletion simulation was performed using OpenMC, with the same time integration and power normalization scheme, to obtain a reference solution. The OpenMC solution used 10⁵ particles per generation over a full assembly (as opposed to a quarter assembly in Scarabée), and 100 active generations per transport step, with an additional 10 generations to permit for source convergence. An LE/QI predictor-corrector time integration scheme was chosen, matching Scarabée’s only provided depletion method. The depletion simulation was carried out to a burn-up of 30 MW d kg⁻¹ using five steps of 0.01 MW d kg⁻¹, one step of 0.15 MW d kg⁻¹, one step of 0.8 MW d kg⁻¹, and twenty-nine steps of 1 MW d kg⁻¹, which can be seen at the bottom of Listing 2. Power normalization was performed using only the average energy released from fission (obtained from the intermediate ACE files when the Scarabée nuclear data library is produced), and identical values were used by both codes to ensure an equivalent comparison. A linear power ratting of 42 kW cm⁻¹ was used for the full assembly. The simplified depletion chain proposed by the CASL project was used for the simulation, albeit with a few modifications [64]. Nuclides with a half-life of less than 24 h were removed from the chain and replaced with their decay products. The two exceptions to this were ¹³⁵I and ¹³⁵Xe, which are kept due to the importance of ¹³⁵Xe poisoning in reactor cores. Additionally, nuclides for which no evaluation exists in ENDF/B-VIII.0 were also removed and replaced with their immediate decay products; this would include ¹¹⁵Cd, ¹⁰²Rh, ^102mRh, ¹²⁷Sb, and ⁸²Br. The resulting chain contains 185 nuclides in total, and fuel which starts as enriched UO₂ will evolve with 165 tracked nuclides.

Listing 1. Portion of the Python script to define all materials used in the simulation.

Listing 2. Portion of the Python script to define all fuel pins, the assembly, and perform the depletion simulation.

The depletion results were obtained using the ENDF/B-VIII.0 nuclear data library in a SHEM-281 group structure, and are presented in Figure 9. In Figures 9(a) and 9(b), a comparison of the k_eff value over the burnup of the assembly is presented. For the duration of the simulation, Scarabée is well within 200 pcm of the reference OpenMC values, and generally within 100 pcm, even up to 30 MW d kg⁻¹. Figures 9(c) and 9(d) show the assembly averaged densities of ²³⁵U and ²³⁹Pu, respectively. For both nuclides, the relative difference with OpenMC is well within 0.5% for the duration of the simulation. Similarly, the densities of important poisons, such as ¹³⁵Xe and ¹⁵⁵Gd are provided in Figures 9(e) and 9(f). Apart from a discrepancy at the beginning of the burnup sequence, when the ¹³⁵Xe density was still very low before achieving equilibrium, the relative error is well within 1%. ¹⁵⁵Gd is the primary absorber which is consumed in the Gadolinia fuel pins, and the relative difference with the OpenMC atom density is below 2%, until the ¹⁵⁵Gd is almost completely consumed at the time of maximum reactivity, where the relative error reaches its maximum at just shy of 8%. These results demonstrate that Scarabée is capable of achieving very accurate results for PWR assembly depletion problems, and at a far lower computational cost than using a Monte Carlo code such as OpenMC. As a point of comparison, the OpenMC depletion simulation required 6387 minutes (4 days, 10 hours, 27 minutes) to complete on a PC with an i9-13900K CPU (32 threads) and 32 GB of RAM, while the Scarabée simulation only took 30 minutes on a Lenovo Yoga Pro 9 16IMH9 laptop, both using shared memory parallelism.

4.4. 2D BEAVRS benchmark

As a final example, a 2D variation of the BEAVRS PWR benchmark was modeled [59]. This simulation assumed Hot Zero Power conditions, with a simplified reflector model that only considered the baffle and moderator; the core barrel and neutron shields were neglected for the purposes of this demonstration. All assemblies were modeled with the ENDF/B-VIII.0 library in a SHEM-281 group structure and had a P1 homogeneous leakage spectrum correction. The full core nodal diffusion calculation was performed using 4 nodes per assembly (2 × 2) in a 2-group structure. Following the advice of Smith, the cross sections for the two assemblies with asymmetric burnable poison configurations were generated independently for each quadrant of those assemblies [27]; each independent quadrant therefore had unique ADFs on each side, as well as unique form factors. While Smith specifically signaled the asymmetric 15 burnable poison pin assemblies, the 6 poison pin assemblies on the core periphery, with the poisons toward the core center, should also receive this treatment. It was initially found that while this treatment yields excellent results for the assembly powers, it produces less accurate results for the reconstructed pin powers in the asymmetric assemblies. The reason for this is that if independent quarter assemblies are not used for all assemblies, then the flux reconstruction system does not have access to the most accurate CDFs for the symmetric assemblies at their mid-points along an edge. Generating “independent” quadrant data for a symmetric assembly comes at no additional computational cost (a single transport calculation on a 1/4 assembly is needed in either case), but was found to greatly improve the errors in pin powers. We have also observed that this treatment is not necessary if a core only has symmetric assembly configurations.

The reference OpenMC solution, using 10⁶ particles with 2000 active generations and 200 inactive generations, obtained a multiplication factor of k_eff = 1.00476 ± 0.00001, while the Scarabée nodal diffusion solver reported k_eff = 1.00406, a difference of only 70 pcm. While the OpenMC results initially seemed well converged regarding uncertainty and entropy, there was still a clear asymmetry in the predicted assembly and pin powers across the full core. Therefore, for all comparisons made in this paper for both the assembly and pin powers, the OpenMC results were averaged across all eight octants to enforce symmetry and to eliminate unphysical “tilts” in the relative differences not caused by the Scarabée tool chain itself. The relative assembly powers of a single octant predicted by Scarabée are presented in Table 1, along with the relative % difference with respect to the OpenMC reference. The root mean squared difference in assembly powers was only 0.37% compared to the OpenMC reference, with a minimum error of -0.68% and a maximum error of 0.86%. The assembly with the largest relative error is the highly asymmetric 15 poison pin assembly, located at H13 in Table 1. With such irregular poison configurations, it is likely that a rehomogenization model would be needed to achieve better results [65–67].

Regarding the pin power distributions, OpenMC predicted a maximum pin peaking factor of 1.545, and a minimum pin peaking factor of 0.207. Using Scarabée, the maximum pin peaking factor was 1.547, with a minimum of 0.206. Figure 10 provides a detailed comparison of the pin power distributions predicted by the two codes. There was a 0.82 % root mean squared difference between Scarabée and OpenMC, with the minimum relative difference being –3.74% and the maximum being 5.75%. As can be observed in Figure 10(b), the largest pin power errors occur at the core-reflector interface, and in the asymmetric assemblies.

5. Conclusions

Scarabée is a new deterministic neutron transport library, providing users with a variety of transport and diffusion solvers. It has tool chains to operate as a full lattice physics code with depletion capabilities, but users can also easily write their own computational schemes. It takes the form of a Python library, giving users access to the full power of the Python scripting language, as well as many other useful open source projects. While indented for use from Python, it is mostly written in C++, ensuring high performance and reasonable run times. It currently targets the traditional two-step computational scheme for LWR analysis, and therefore mostly contains 1D and 2D transport solvers, and 3D simulations (generally assumed to be at the core level) can be performed in diffusion theory. A notable difference from several other open source deterministic tool chains is that Scarabée can use a custom mutli-group nuclear data library, and perform self-shielding calculations to obtain appropriate multi-group cross sections for many different LWR problem geometries; there is no need to run a Monte Carlo simulation before hand to generate these cross sections.

Intended to be an easily accessible and usable tool for teaching nuclear engineering students, it is free and open source software, licensed under the GNU LGPLv3 license. While the very first version was released under the GNU GPLv3 [68], the authors determined that it would likely be more beneficial to the community to use the LGPL, as it would facilitate the library’s use in other open source projects that might use a permissive license. Pre-compiled Python wheels are made available on the Python Package Index, PyPI, and can easily be installed with the pip command on Windows, macOS, and Linux, making adoption in the classroom very easy. The source code for the library can be found at https://github.com/scarabee-dev/scarabee, and the documentation at https://scarabee.readthedocs.io. Many future developments are planned, such as the generation of case-matrices and reactor databases, the handling of IFBA layers on fuel pellets, as well as the creation of a full fledged “core simulator”. Long term objectives include support for other reactor types and geometries such as BWR assemblies, hexagonal VVER / fast reactor assemblies, and CANDU fuel clusters. Other developers are more than welcome to contribute to the code, as the hope is that Scarabée can become a thriving community driven project facilitating reactor physics training and education.

Funding

This research received no external funding.

Conflicts of interest

The authors have no conflicts of interest to disclose.

Data availability statement

No datasets are associated with this article.

Author contribution statement

H. Belanger: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Writing – Original Draft Preparation, Writing – Review & Editing, Visualization, Supervision; P. Singh: Software, Validation, Writing – Original Draft Preparation; J. Gaiardelli: Software, Validation, Writing – Original Draft Preparation;

Acknowledgments

We would like to thank Prof. Kord Smith for his valuable advice and insightful suggestions on modeling the 2D BEAVRS benchmark, which enabled the excellent level of agreement achieved with the Monte Carlo reference calculations.

References

All Tables

All Figures

	Fig. 1. Depiction of how Pybind11 is used to interface low-level C++ classes in Scarabée with the higher level features written in Python.
In the text

	Fig. 2. Depiction of the HDF5 file structure for the microscopic nuclear data used by Scarabée.
In the text

	Fig. 3. Depiction of how the optical depths, τi, j ±(y), are computed as a function of y in the 1D annular collision probabilities solver.
In the text

	Fig. 4. Diagram of the simplified 1D reflector model used to generate few-group diffusion cross sections and discontinuity factors for reflector nodes. This scheme is inspired by the Baff-Ref method [37, 38].
In the text

	Fig. 5. Comparison of the flux spectrum from Scarabée and OpenMC for a 1D annular pin cell problem using Scarabée’s 1D annular collision probabilities solver.
In the text

	Fig. 6. A portion of the C5G7 geometry in the graphical plotter tool for examining MOC simulation geometries.
In the text

	Fig. 7. Thermal flux from the C5G7 benchmark predicted by Scarabée.
In the text

	Fig. 8. Convergence of the C5G7 benchmark problem with and without CMFD to accelerate flux convergence.
In the text

	Fig. 9. Comparison of depletion results from Scarabée with OpenMC.
In the text

	Fig. 10. Comparison of the pin power distributions obtained by OpenMC and Scarabée for the 2D variant of the BEAVRS benchmark.
In the text

	Listing 1. Portion of the Python script to define all materials used in the simulation.
In the text

	Listing 2. Portion of the Python script to define all fuel pins, the assembly, and perform the depletion simulation.
In the text