© C. Takasugi et al., Published by EDP Sciences, 2023

1. Introduction

Griffin is a Multiphysics Object-Oriented Simulation Environment (MOOSE) based reactor-multiphysics code jointly developed by Idaho National Laboratory and Argonne National Laboratory that is applicable for both steady-state and transient analyses. It leverages the transport solvers from its predecessors Rattlesnake [1] and Proteus [2]. Griffin also supports coarse-mesh diffusion calculations with equivalence techniques, such as superhomogenization (SPH) or discontinuity factors (DF), which are routinely used in production reactor physics codes. Griffin has the flexibility to model any geometry and, when employed with the proper set of cross sections, can model virtually any type of reactor spectrum, making it a tool of choice for modeling advanced reactor designs. The Monte Carlo code Serpent [3] has been employed at Idaho National Laboratory for the generation of multigroup cross sections and the kinetics parameters, which are used in the downstream Griffin calculations.

The homogenized few group cross sections and kinetic parameters used in Griffin are generated from the full 3D geometry. The SPH equivalence technique relies on non-fundamental mode conditions [4] and is implemented as an intermediate step between the Monte Carlo calculation and low-order core calculations. A 3D macro-geometry is used to generate reference multigroup fluxes in Serpent, and SPH factors are computed via a Griffin SPH Preconditioned Jacobian-Free Newton Krylov (PJFNK) procedure for each statepoint [5]. One key specificity of this approach is that the SPH factors can optionally be computed for only a portion of the macro-geometry (at minimum, the fueled region) and still result in a preservation of the k-eigenvalue and reaction rates [5]. However, early applications of the Serpent–Griffin computational scheme observed large sensitivities of kinetics parameters (i.e., effective delayed neutron fraction β_eff and mean generation time Λ) to the computational scheme options selected for homogenization and the generation of group constants. Specifically, the energy grid used for the cross-section generation step in Serpent and thus in the Griffin diffusion calculations, as well as the regional application of SPH equivalence factors to preserve reaction rates could significantly impact the preservation of these kinetic parameters.

Previous studies have investigated the impact of bilinear, adjoint-weighted quantities through homogenization and energy group condensation processes for cross section generation. For example, results in [6] indicate no significant difference between linear and bilinear weighting methods for steady-state calculations, however demonstrated the use of bilinear weighting in determining β_eff to improve transient power predictions. In this previous work, the option to use either a leakage adjoint spectrum or an infinite medium adjoint spectrum were assessed, however this consideration is not necessary within the Serpent–Griffin scheme as the cross section generation is done on the full 3D geometry. Further, the impact of equivalence procedures such as SPH can result in modified flux values in the homogeneous solution used for the bilinear weighting of effective kinetics parameters, and has not been studied before. Hence, this work aims at evaluating the impact of homogenization and condensation with non-fundamental SPH equivalence on adjoint-weighted quantities.

For this study, the IPEN/MB-01 reactor kinetics benchmark, available through the International Reactor Physics Experiment Evaluation Project handbook and described in [7], was selected for evaluation by the Serpent–Griffin computation scheme. This kinetics benchmark has been evaluated in full or in part by several codes with various nuclear data libraries, including MCNP (ENDF/B-VII.0, ENDF/B-VII.1, JEFF-3.1.1, JENDL-3.3, JENDL-4.0) [7], CASMO5 (ENDF/B-VII.1, Tuttle, JENDL-4.0) [8], MORET 5.D.1 (ENDF/B-VII.1) [9], MPACT (ENDF/B-VII.1, Keepin et al., Tuttle, JENDL-4.0, JEFF-3.3, Santamarina et al. [10]) [11], TRIPOLI-4 (JEFF-3.1.1) and APOLLO2.8 (JEFF-3.1.1, ENDF/B-VII.0) [10]. In addition to the benchmark model used in these studies, which is simplified from the as-built data, the IPEN/MB-01 reactor has been used for kinetics experiments of varying configurations. As a result, some code evaluations have been validated against different benchmark kinetics parameters (e.g., MCNP [JEFF-3.1 and ENDF/B-VII.0] [12, 13] and Serpent [JEFF-3.1.1] [3]) and may not be directly comparable due to model differences.

The MCNP study [7], which made comparisons across five nuclear data libraries, concluded that the β_eff and Λ evaluations were satisfactory across all libraries; however, reactivity estimates were not. The reactivity estimates based on the Inhour equation varied significantly between nuclear data libraries, including up to ≈20% relative difference from the benchmark values using the ENDF/B-VII.0 library and up to ≈50% relative difference with ENDF/B-VII.1. These differences were attributed in large part to the difference in decay constants between the libraries, demonstrated by reducing the relative difference in ENDF/B-VII.1 reactivity estimates to a maximum of ≈10% by fixing the experimental first delayed group decay constant λ₁. The CASMO5 evaluation [8] further supports the observations made in the MCNP reactivity estimates adding that, while the discrepancies in ENDF/B-VII.1 results are primarily due to λ₁, λ₂ and λ₃ also play a significant role.

In addition to the comparison to the benchmark as validation of the Serpent–Griffin scheme, variant test cases are modeled based on the benchmark to compare different Serpent kinetics parameter estimation methods and to study the impact of modeling assumptions on the generation of effective, adjoint-weighted kinetics parameters in Griffin.

Section 2 provides an overview of the Serpent–Griffin computational scheme, describing the typical procedure for a transient calculation and relevant information on kinetics parameter generation. Section 3 describes the IPEN/MB-01 reactor kinetics benchmark used for analysis of modeling options and compares the results of test cases against benchmark values. Section 4 discusses the modeling options and the introduced differences within the context of test case results. Section 5 describes supplementary test cases based on models of the Godiva and Flattop critical experiments. Section 6 provides a summary of the conclusions and potential future extensions of this work.

2. Serpent–Griffin computation scheme

This study focuses on the assessment of kinetics parameter preservation between Serpent and Griffin in a two-step computation scheme that has been originally developed for modeling transient testings in the Transient Reactor Test (TREAT) facility at Idaho National Laboratory. The numerical methodology employed was first described in [14–16] and is summarized below. In this scheme, Serpent is used to generate multigroup cross sections and kinetics parameters as well as provide equivalence data for downstream Griffin calculations. These group constants are used in a series of Griffin calculation steps leading to Griffin transient calculations. The following outlines the sequence of calculations to prepare a standard Griffin transient calculation:

• 1.

  group constant generation with Serpent. This step consists of a Serpent calculation on the full 3D geometry to provide multigroup cross sections and kinetic parameters, spatially homogenized and condensed to a few energy groups. Reference fluxes are also computed for the equivalence procedure. This calculation is performed for each statepoint, defined as a full combination of all the parameters deemed to have an important impact on cross sections and/or equivalence factors (burnup, fuel temperature, moderator temperature, etc.), and the resulting dataset is stored in cross section and equivalence libraries generated by the ISOXML module in Griffin. For this study, selected benchmarks are focused on kinetics parameters rather than actual transient evaluations, and therefore only a single statepoint is considered at benchmark specified conditions.

• 2.

  Intermediate equivalence procedure with Griffin. SPH factors are generated by the Griffin SPH Preconditioned Jacobian-free Newton Krylov (PJFNK) procedure for each statepoint, using the multiplication factor (k_eff) calculated by Serpent and reference fluxes stored in the equivalence library. Different options are possible in terms of where the SPH correction is applied, which can range from all regions in the core to individually selected regions, such as fuel assemblies or specific fuel-containing elements. As long as all fuel-containing regions are SPH corrected, the eigenvalue, power distribution, and total leakage out of SPH regions can be reproduced [5]. In this study, a few different options are explored for regional application of SPH equivalence procedures to assess the effect on kinetics parameters determined by Griffin.

• 3.

  k-eigenvalue calculation with Griffin. In most cases, this calculation step relies on a coupled multiphysics simulation to obtain the right thermal-hydraulic or thermo-mechanical feedback (accounted for via changes in the cross sections) on the Griffin neutronics calculation. The benchmarks selected for this study are taken at prescribed conditions, and therefore no consideration of multiphysics feedbacks is necessary.

• 4.

  k-eigenvalue adjoint calculation with Griffin. The adjoint form of the k-eigenvalue transport equation is solved and results in an identical eigenvalue to that found in the previous step. This step is only needed when point kinetics parameters are required in Griffin, as it provides the adjoint flux to be used in the bilinear weighting of the kinetics parameters. Note that additional simplifications specific to the Griffin Improved Quasi-Static (IQS) implementation might be applied here, notably the adjoint flux is not updated at every time step. These updates to the adjoint are irrelevant for this study since there is no perturbation to the core conditions, however the adjoint calculation with Griffin is necessary to evaluate kinetics parameters.

• 5.

  Null transient calculation with Griffin. The Griffin calculation restarts from the previous one, except the time-dependent diffusion (or transport) equation is solved instead of the k-eigenvalue one. The fission operator is divided by the k_eff obtained so as to preserve the neutron balance in the diffusion operator. Since no perturbation is applied, the solution should be time insensitive. This step is intended as a verification of the numerical stability and is important for catching issues or discrepancies between the k-eigenvalue model and the time-dependent one.

• 6.

  Transient calculation with Griffin. The time-dependent Griffin calculation restarts from the k-eigenvalue one. The fission operator is divided by the k_eff obtained so as to preserve the neutron balance in the diffusion operator. A perturbation (movement of a control rod, thermal-hydraulic condition change, etc.) is introduced and Griffin simulates the neutronics response of the system as a result of the perturbation. This step is not carried out within the present study.

• 7.

  An alternative step is to perform an IQS transient simulation. First, this calculation employs the forward and adjoint flux solution obtained in steps 3 and 4 to collapse the core-wise kinetic parameters using a bilinear product. The core-wise kinetic parameters (Λ – mean generation time, β – delayed neutron fraction, and ρ – reactivity) can be injected into the Inhour equation to obtain the corresponding period (T) of the reactor. Since Griffin can output the point kinetics parameters, this step is useful when performing comparisons to transient experiments for which the measured reactor period is available and to determine a corresponding calculated reactivity that will lead to the same reactor period in the Griffin model. These core-wise kinetic parameters can also be employed in the point kinetics solver available in Griffin or internally by the IQS solver, which relies on a factorization of the flux into a 3D shape and a 0D time-dependent amplitude, obtained by solving the point kinetics equations [17]. For this study, this step is used only to collapse the kinetics parameters for comparison to benchmark studies.

Step 4 is only necessary when core-wise point kinetics parameters are involved, which is the case in many reactor physics applications such as the modeling of small experimental reactors like TREAT. Another application of bilinear-weighted kinetic parameters is when an IQS transient is performed, which provides access to a 3D time-dependent flux solution while being less computationally intensive than a regular space-time kinetics calculation by decoupling the flux into a product of a three-dimensional shape and a time-dependent amplitude solved via a regular point kinetics model. Thus, it is crucial to demonstrate that Griffin can preserve the kinetic parameters obtained from Serpent, in order to not introduce any additional error during the resolution of the amplitude function compared to a point kinetics calculation utilizing the straight Serpent kinetics parameters. The current study assesses the preservation of the kinetics parameters employed within the IQS transient option of the described computation scheme, carrying out only the portions of the typical procedure necessary to obtain estimates from Griffin. Note also that no perturbations are introduced which would result in transient updates to the neutron flux. For the critical experiments modeled in this work (IPEN/MB-01, Godiva, Flattop-23), the same computational approach as for TREAT is used, except that the process stops at Step 5 where Griffin computed kinetic parameters can be compared against experimentally measured ones.

In the described workflow, equivalence techniques, such as SPH, are available to guarantee that the Serpent and Griffin 3D calculations match for the homogenized, few group fluxes, and reaction rates (steps 1, 2, and 3). Note that these homogenized group constants (cross sections, diffusion coefficients) are all generated via a direct flux-volume weighting:

(1)

and

(2)

(4)

(6)

where and are the multigroup cross sections and flux from Serpent, respectively (thus denoted “reference”), and where the integral is performed over region m corresponding to the spatial domain r ∈ V = V_m being homogenized and an energy domain E ∈ [E_g, E_g − 1] corresponding to each energy group. Equivalence procedures, such as the one implemented in Griffin, enforce the preservation of the reaction rates for individual SPH regions m and energy group g via the introduction of an SPH factor μ_m, g:

(3)

where ϕ_m, g is the integral of the flux obtained from the Griffin diffusion calculation over the homogenized SPH region V_m.

Kinetic parameters are typically collapsed via a bilinear forward/adjoint flux weighting. In Serpent, effective kinetics parameters are derived using a Monte Carlo estimate of the following mathematical definitions [3], which can be written using notations from [18]:

(5)

and

where:

• is the effective delayed neutron fraction in homogenized region m and precursor group j.

• is the effective mean generation time.

• v is the neutron speed.

• ν is the average total number of neutrons emitted per fission.

• is the energy spectrum of delayed neutrons from precursor group j from isotope i.

• χ is the energy spectrum of neutrons (both prompt and delayed).

• β_j ⁱ is the delayed neutron fraction in precursor group j from fission in isotope i.

• Σ_f ⁱ is the macroscopic fission cross section for isotope i.

• N is the number of isotopes.

• N_d is the number of delayed neutron precursor groups.

• ϕ^* is the adjoint neutron flux.

• ϕ is the neutron flux.

For equation (4), corresponding to the process used by Serpent for computing the delayed neutron fractions, the integral over V can be either the full domain, leading to a core-averaged value , or to a value per fueled homogenized region V_m, which is denoted . Serpent produces both values, the former from Meulekamp’s method leading to one set of delayed neutron fractions per fueled homogenized region m and the latter from Iterated Fission Probability (IFP) and Nauchi methods which provide only core-wise kinetics parameters. Griffin space-time kinetics calculations can then directly use these on a spatially-varying basis (if taken from Meulekamp’s method), or on a core-wise basis if taken from the other approaches. For cases in this work, only one fuel region with fresh fuel is homogenized into a single region, so all the Serpent methods for estimating the effective delayed neutron fractions produce only one set of . For the mean generation time in equation (6) the integration is performed over the full core.

The generation of effective kinetics parameters in continuous-energy Monte Carlo codes such as Serpent does not involve formally a bilinear product, which would require determining the adjoint flux, but rather proceed via specific Monte Carlo estimates of the importance of the neutrons in the fission chain, as discussed in more detail in Section 3.1. The core-wise effective kinetics parameters from Serpent and can then be directly used in point-kinetics calculations. As discussed above, Griffin can also employ a bilinear weighting to compute core-wise effective kinetics parameters for use in IQS method or for comparison against measured values, using:

(7)

and

(8)

where:

• β_eff, j is the effective delayed neutron fraction for precursor group j computed by Griffin.

• Λ_eff is the mean generation time computed by Griffin.

• is the effective delayed neutron fraction for precursor group j for homogenized region m defined in equation (4) and computed by Serpent.

• is the energy spectrum of delayed neutrons from precursor group j in energy group g homogenized on spatial domain m from Serpent. It is averaged by Serpent over all fissile isotopes.

• χ_m, g is the total energy spectrum in group g homogenized on spatial domain m from Serpent. It is averaged by Serpent over all fissile isotopes.

• μ_m, g is the SPH-factor computed during the equivalence procedure in Griffin for homogenized region m.

• νΣ_{f, m, g′}^ref the nu-fission macroscopic cross-section on homogenized region m, from Serpent.

• is the neutron velocity in group g on homogenized region m, from Serpent.

• N_G is the number of neutron energy groups.

• M is the number of homogenized regions.

• ϕ_g^* is the Griffin multigroup adjoint flux.

• ϕ_g is the Griffin multigroup flux.

The spatial integration in equations (7) and (8) is carried in two steps, first by integrating over each spatial domain for which homogenized group constants are available (note the spatial dependency of the Griffin adjoint and direct fluxes), and second by summing over these integrals over all the different homogenized zone. If considering an IQS transient with perturbations, the kinetics parameters will be recalculated with updated fluxes based on cross sections interpolated between multiple statepoints according to the perturbed conditions.

As opposed to the reaction rates which are preserved via SPH equivalence theory, there is no guarantee that the kinetic parameters from Serpent are actually preserved in the Griffin calculations. In other words, applying a bilinear weighting using the forward and adjoint multigroup flux computed by Griffin to collapse the kinetic parameters β_eff, j and Λ_eff do not necessarily guarantee preservation of the reference values obtained from Serpent. SPH equivalence theory provides a way to preserve the reaction rates via a modification of cross sections, but the theoretical aspects involving the homogenization of adjoint-weighted quantities, such as kinetic parameters, are lacking. Indeed, preserving in Griffin the straight kinetics parameters (non-adjoint weighted) from the Serpent calculation is supported by equivalence theory since these quantities are ratios of reaction rates (themselves preserved by the SPH equivalence process). However, there is no guarantee that the adjoint flux between both codes will be the same, thereby inducing potential errors between both codes for bilinear-weighted quantities. Considering the physical interpretation of the adjoint flux as a measure of importance of a single neutron to the continuation of the fission chain reaction, there may be differences between both adjoint flux solutions even though the fission reaction rates are identical. Other effects such as energy group condensation may also contribute to differences between the Griffin and Serpent kinetics parameters. Another complexity stems from the usage of SPH-corrected multigroup cross sections to calculating the adjoint flux in Griffin. In this case, the adjoint flux in equations (7) and (8) becomes ϕ_m, g^*/μ_m, g as a consequence of using SPH-corrected cross sections in the steady-state adjoint calculation, which can be verified by writing down the adjoint operator with SPH-corrected cross-sections [19]. Thus, one of the objectives of this work is to assess the impact of the SPH equivalence process on kinetics parameters calculated by Griffin.

One additional caveat to be aware of is that, since the Griffin transient calculation starts from the steady-state fission source normalized by k_eff to preserve the neutron balance at the start of the transient, the mean generation time from Griffin needs to be divided by k_eff. This can be deduced by observing equation (8), where the fission source term appears at the denominator. For β_eff, the impact of this normalization cancels out.

3. IPEN/MB-01 reactor kinetics benchmark

To investigate the validity of the described computational scheme and the effect of model settings, the IPEN/MB-01 research reactor facility was modeled according to the simplified benchmark model described in the International Reactor Physics Experiment Evaluation Project handbook and in a separate document by the same authors [7]. The benchmark includes kinetics parameters as shown in Table 1 as well as reactivity values corresponding to various reactor periods listed in Table 2.

The Serpent–Griffin computational scheme was used to generate kinetics parameters and a series of reactivity values for comparison with the benchmark reference values. Specifically, the Serpent-evaluated kinetics parameters were compared to benchmark values as a validation step. This step provides the relevant information for determining which Serpent method should be used to generate delayed neutron data for downstream Griffin calculations. A second, independent step was to verify the degree to which IQS transient-predicted Griffin kinetics parameters preserve the Serpent-generated ones. The process was repeated with model variations in order to observe the changes in kinetics parameters as a result of:

• Different SPH equivalence schemes

  • No SPH equivalence

  • SPH factors computed on the core only

  • SPH factors computed on the core and part of the reflector region

• Energy group count and energy grid selection

  • CAS2 predefined two-group energy grid

  • CAS4 predefined four-group energy grid

  • CAS8 predefined eight-group energy grid

  • C4G custom four-group energy grid based on CAS2 and refinement of thermal energy group present in the CAS7 predefined energy grid

Through this series of tests, a better understanding of modeling assumptions may be achieved and provide recommendations for future use of the Serpent–Griffin computational scheme.

3.1. Serpent kinetics parameters

Serpent provides different estimates of effective kinetics parameters, β_eff and Λ, and results for each are printed to the Serpent output file [3]. The three Serpent methods considered in this study are Meulekamp, Nauchi, and IFP:

• the Meulekamp method approximates β_eff as a fraction of the fissions caused by delayed neutrons. This is done by labeling neutrons as either prompt or delayed at birth, tracking the number of fissions caused by delayed neutrons and comparing it to the total number of fissions [20]. It is the simplest technique for approximating β_eff in a Monte Carlo calculation, and it does not rely on an adjoint flux estimator, contrary to the two other techniques listed below.

• The Next Fission Probability (NFP) method estimates the expected number of fissions of a single neutron, dependent on position, energy, and direction, which is then utilized as an importance weighting function in place of the adjoint flux in kinetics parameter estimations. This technique is termed the Nauchi method in Serpent [21] and will be referred as such in the remainder of the paper.

• The IFP method is an extension of the NFP method and calculates the iterated fission probability, or expected number of descendent fission neutrons in each generation, as a product of successive k_eff estimates in each cycle of a Monte Carlo calculation. This probability is proportional to the adjoint flux and the method utilizes it as an importance weighting function in kinetics parameter calculations in place of the adjoint flux. Even though both NFP and IFP methods tend to produce similar results, the proportionality of these estimates to the adjoint flux is only mathematically proven for IFP [22].

3.1.1. Model description

The IPEN/MB-01 reactor description and configurations used in the benchmark experiments are found in the IRPhE handbook; specific IRPhE report numbers for the core layout and materials can be found in [7]. The Serpent model was developed based on the simplified benchmark model. The core follows a 28 × 26 configuration, with control rod banks inserted to a critical state. The core includes a lower baseplate and is located within a large cylindrical reflector. The benchmark model is simplified from the as-built conditions, removing thermocouples, instrumentation tubes, etc., which are accounted for in the final benchmark model and specifications. The final geometry used for the Serpent model is shown in Figure 1.

Fig. 1. Serpent model geometry. (a) Top view of core region; (b) side view of core region.

The Serpent model was developed with a stacked lattice for the core and a portion of the nearby reflector in axial and radial directions, as indicated by the grid lines in Figure 1. This allowed for the simple integration of a series of lattice-type detectors for generating the reference multigroup fluxes for use in downstream equivalence procedures in Griffin. A modification of the model was made to simplify this process, resizing the baseplate to align with the extended detector lattice; this modification was verified to have minimal impact on kinetics parameter results. Nine cross-section sets were generating corresponding to different compositions of the lattice-defined homogenized cells: base plate, reflector, alumina, spacers, safety rods, control rods, empty control rod guide tubes, control rod end plugs and fuel rods. Ten million particles were simulated with 1,000 active and 100 inactive cycles. For the IFP method, the Serpent default of 15 latent generations was utilized.

The ENDF/B-VII.1 cross section library was selected for the Serpent-to-benchmark comparison. Redundant Serpent runs were performed from a fixed seed for each separate energy grid to generate group constants for Griffin; details of these energy grids are presented in the Griffin model description. k_eff calculated by Serpent was 1.00271 ± 0.00001.

3.1.2. Validation of Serpent kinetics parameters

Presented in Table 3 are the Meulekamp, Nauchi, and IFP method results and statistical errors, noting that the Meulekamp method does not predict a Λ value. The relative difference in these results are compared against the benchmark values in Table 4. In comparison to the benchmark values, the Meulekamp method provides the best prediction of β_eff, and Nauchi provides a better estimate of Λ than IFP. Group-wise delayed neutron fractions and decay constants from the various Serpent methods are shown in Table 5.

Reactivities corresponding to the benchmark reactor periods were evaluated using the Inhour equation:

(9)

where T_p represents the corresponding reactor period taken from the benchmark, λ_j is the decay constant of delayed neutron group j, and the prompt neutron generation time l^* is the mean generation time Λ. Nauchi and IFP reactivity values corresponding to benchmark reactor periods were evaluated by the Inhour equation and the resulting errors compared against the benchmark values are given in Table 6; Meulekamp is excluded since no Λ value is predicted by Serpent’s Meulekamp method.

Each of the Serpent evaluations of β_eff had less than a 1% relative error against the benchmark values, while Λ was estimated well by the Nauchi method (less than a 1% relative error) and ≈5% relative error using the IFP method. Reactivities predicted by the Inhour equation contain a significant relative error but align well with MCNP and CASMO5 evaluations of the IPEN/MB-01 benchmark utilizing the ENDF/B-VII.1 library [7,8]. The MCNP evaluation found that significant improvement in reactivity estimates may be obtained by utilizing experimental λ₁ value, rather than the one provided in the ENDF/B-VII.1 library.

3.2. Griffin kinetics parameters

Griffin kinetics parameters were compared to the same results from Serpent to observe the degree to which kinetics parameters were preserved through the later parts of the computational scheme. Griffin kinetics parameters that utilize delayed neutron data generated by Serpent’s Nauchi method were compared to the Serpent Nauchi method results. This comparison was repeated for several Griffin variants to observe the effect of model options.

Further comparison of Griffin-evaluated kinetics parameters to the benchmark data was used as a characterization of the net effect of the two-step Serpent–Griffin scheme on kinetics parameters.

3.2.1. Model descriptions

The size of the reflector region for the IPEN/MB-01 reactor kinetics benchmark makes it impractical to perform SPH equivalence procedures on the entirety of the reflector region as low flux in the far-reflector region would result in low reaction rates with large statistical uncertainty in corresponding Serpent detectors, leading to potential poor convergence of the SPH procedure. In addition, the MOOSE native meshing tool, MooseMesh, was determined to be inadequate at the time of this study for development of the complex mesh required for the IPEN/MB-01 reactor. To circumvent these issues, a Cartesian mesh was generated for the core and a near-reflector region utilizing MooseMesh and corresponding SPH equivalence lattice detectors were implemented in the Serpent model. A separate mesh of a cylinder was generated in MooseMesh, as an extruded concentric circular mesh, representing the far-reflector regions. The two meshes were merged and the cylindrical portion was re-meshed with tetrahedral elements as part of the merging process, utilizing Cubit [23]. This results in a final mesh preserving the original Cartesian mesh of the near-core region in agreement with Serpent lattice detectors while maintaining the full, cylindrical reflector region. Materials and equivalence numbers mapping the Serpent-generated cross sections and reference fluxes were assigned to the mesh via a Python script to facilitate SPH within the Cartesian portion of the model. A partial top-view cross section of the resultant final mesh near the core-region is shown in Figure 2.

Fig. 2. Top-view cross section of Griffin mesh.

The Cartesian portion of the mesh includes both the fuel and the near-core reflector regions within which SPH can optionally be performed for all or individual sub-regions such as the fuel alone, the core or portions of the near-core reflector. Indeed, previous studies involving Griffin calculations with Serpent-based cross sections observed a sensitivity of kinetics parameters to the selection of regional SPH application [15]. As demonstrated in [5], preserving the k_eff, power distribution, and leakage out of SPH regions from the 3D Serpent model via the Griffin SPH methodology requires at least the SPH correction of fuel regions. However, no discussion was made as to what is required regarding the preservation of kinetics parameters, especially what could be the impact of restricting the SPH correction to the fuel region only. Thus, to gain some insight into this observation, data was generated both for the SPH equivalence over the entire Cartesian portion of the model and for SPH equivalence applied in only the fuel-containing elements. Finally, the option of a Griffin diffusion calculation without any SPH correction was also assessed.

An additional consideration arose regarding the theoretical applicability of SPH equivalence for an evaluation of the steady-state adjoint flux, since the SPH factors are generated to preserve forward flux reaction rates. The adjoint flux is calculated using the same SPH-corrected cross sections as for the forward one (see step 4 in Section 2), and a further investigation into how SPH equivalence at this specific step may affect the kinetics parameters was also desirable.

The Griffin-evaluated kinetics parameters are significantly dependent upon the group constant energy grid selection. Four energy grids were evaluated in this study: CASMO 2-group (CAS2) predefined energy grid, CASMO 4-group (CAS4) predefined energy grid, CASMO 8-group (CAS8) predefined energy grid, and a custom 4-group energy grid. The CAS2 energy grid was able to reproduce the effective delayed neutron fraction from Serpent, while CAS4 and CAS8 were not. This qualitative difference in error present in Griffin-calculated delayed neutron fractions is a consequence of the refinement of the epithermal and fast neutron energy groups in the CAS4 and CAS8 energy grids. A custom 4-group energy (C4G) grid was included to provide an additional point of comparison to CAS2 for mean generation times, and its development was based on the refinement of thermal energy groups present in the CASMO 7-group (CAS7) energy grid, while neglecting the epithermal and fast energy group refinement. The four energy grids utilized in this study are presented in Table 7.

3.2.2. Verification of kinetics parameters with Serpent

The relative error of Griffin kinetics parameters against corresponding Serpent Nauchi kinetics parameters for all tested methodology variants are shown in Table 8. Separate Serpent evaluations utilizing the same seed were used for group constant generation and SPH equivalence in each energy grid. Note that entries of zero indicate identical values when Griffin results are rounded to the last displayed Serpent digit.

The following observations may be made from the results in Table 8:

• β_eff was identically replicated with the selection of energy grids CAS2 and C4G. Error is introduced in β_eff due to the refinement of the epithermal/fast energy groups in CAS4 and CAS8.

• SPH equivalence procedures were observed to only significantly affect β_eff when SPH is applied to the adjoint. Cases without SPH in the adjoint yielded similar β_eff values to non-SPH cases.

• The relative error in Λ was significantly reduced by applying SPH to the fuel only (rather than the fuel+reflector portion of the model). This effect became less pronounced with the refinement of the energy grid.

• Errors in Λ were reduced by the refinement of the energy grid, with the exception of CAS2 cases with only fuel-region SPH. Improvement in β_eff is observed in moving from CAS4 to CAS8.

The Inhour equation, equation (9), was utilized with Griffin mean generation time results to determine reactivities for the benchmark-specified reactor periods. Because Serpent kinetics parameters were, in general, reasonably preserved by Griffin (within a few percent in most cases) and delayed neutron data is provided directly by Serpent, reactivities were also well preserved. The cases with minimum and maximum relative errors of each energy grid are compared with Serpent reactivity estimates in Table 9. The specific case corresponding to the minimum or maximum relative error for each energy grid is indicated by a letter code.

Results of the reactivities comparison are consistent with the kinetics parameter observations. For most cases, the minimum deviation from Serpent reactivities was found in cases with SPH in the fuel only that did not include SPH in the adjoint calculation step. Further, the energy grids with the refinement of epithermal/fast energy ranges similarly introduced further deviation from Serpent results. In general, however, the reactivities are only slightly changed – in none of the cases was a greater than 0.1% change from the Serpent reactivity observed. This is expected due to the primary dependence on delayed neutron data taken directly from Serpent.

3.2.3. Validation of Serpent–Griffin scheme kinetics parameters

The validation of Serpent kinetics parameters against benchmark values and verification of Griffin kinetics parameters against Serpent ones provide useful comparisons for the analysis and assessment of modeling options. Griffin, however, was also compared against the benchmark values to provide a validation of the two-step Serpent–Griffin scheme. Relative differences for kinetics parameters of this comparison are provided in Table 10.

The data presented shows the impact of error compensation between the two steps; where Griffin and Serpent results have opposing error values, the net result is less extreme than in either case. This is evident in many of the Λ values, which were over-predicted by Serpent’s Nauchi method, but Griffin’s results were lower than Serpent ones. Because Serpent Nauchi results were generally close to the reference benchmark values, the relative difference in the kinetics parameters originates primarily from the Griffin evaluation, with the exception of β_eff in CAS2 and C4G where it was perfectly preserved.

Comparisons to experimental reactivities are not presented for Griffin due to observations made in the analysis of individual steps of the calculation scheme. Specifically, near all relative difference from the benchmark reactivity values is attributable to Serpent and the source library delayed neutron data that remains unchanged between Griffin and Serpent evaluations (i.e., β_j and λ_j differences from the benchmark values). The evaluation of core homogenized kinetics parameters in Griffin introduces only a small additional error as shown in Table 9, and results for Griffin closely resemble the Serpent results presented in Table 6.

4. Discussion

Kinetics parameters evaluated by various Serpent methods were compared against the IPEN/MB-01 reactor kinetics benchmark values and found to be within a few percent relative error in most cases. Reactivities were estimated by the Inhour equation, based on kinetics parameters and generated delayed neutron data. These reactivities had significant error compared with benchmark values, particularly for large or negative reactor periods, but were consistent with MCNP evaluations of the benchmark utilizing one of the same nuclear data libraries.

Serpent reactivities (as evaluated by the Inhour equation) were well preserved by Griffin in all cases. This may be due in part to the equation being strongly dependent on delayed neutron data, which comes from Serpent, with only a small dependence on the introduced difference in Griffin’s generation of kinetics parameters. Though only a small relative error was observed, much of the qualitative behavior from kinetics parameter comparisons remained observable. Despite significant differences in kinetics parameters between Serpent and Griffin, the reactivities corresponding to various periods, and therefore expected transient behavior, between the two codes remain well-aligned due to the shared delayed neutron data.

4.1. Energy grid sensitivity

The significant sensitivity of kinetics parameters to the number of energy groups was similarly observed in accordance with previous studies and was observed to be mitigated by specific energy grid selections. An observation of specific interest is that a qualitative shift occurred from effectively no difference in β_eff in the two group and custom four group thermal-energy refined cases to large differences in those with refinement of epithermal/fast energies. The model for collapsing Griffin’s effective delayed neutron fraction explains this qualitative behavior.

Griffin calculates core-averaged β_eff by summing β_eff, j terms calculated via a ratio of a bilinear weighting of the delayed fission rates over the total fission rates, utilizing both the forward and adjoint flux as shown in equation (10).

(10)

The definition of the average fission spectrum used in the denominator, χ_m, g, is as follows (the index m is dropped for clarity):

(11)

with

(12)

Since Serpent provides a single delayed spectrum for all the delayed precursor groups,

(13)

and equation (11) becomes:

(14)

In cases where , that is when the prompt and delayed spectrum are identical per energy group equation (14) simplifies further. For example, the four-group and eight-group energy grids used in this study each have distinct prompt and delayed spectra and no simplifications can be made. However, in the two-group and custom four-group energy grids both prompt and delayed components exclusively produce neutrons in the highest energy group, allowing the following simplification:

(15)

and in this situation, equation (10) simplifies to:

(16)

For this reason, β_eff was perfectly preserved in the CAS2 and C4G energy grid portions of the Griffin verification step. However, whenever the prompt and delayed spectrum are not identical per group, this simplification is not possible and β_eff in Griffin is not necessarily identical to that of Serpent.

Since the qualitative shift caused by the mismatch between the prompt and delayed spectra is a major source of the relative difference observed in the overall Serpent–Griffin scheme, it may be desirable to mitigate this issue by optimally selecting the energy grid. These spectra are fuel-dependent, and it should be possible to assess energy grid structures for a degree of conservation, even with simplified models. It may be possible to develop an infinite homogeneous model containing fuel representative of the complex model and carry out a test of various candidate energy grids similar to this study. The results of these homogeneous tests may compare the Griffin kinetics parameters against the Serpent ones to down-select or optimize energy grids for the preservation of kinetics parameters in the desired complex model transient.

4.2. SPH equivalence domain sensitivity

The sensitivity to regional application of SPH equivalence procedures found consistent improvement in Griffin’s preservation of mean generation time when applied only within the fuel. The cases that had a relative difference in effective delayed neutron fraction resulted in a slight increase in difference when SPH was only performed in the fuel region.

An application of SPH equivalence in the steady-state adjoint flux calculation step was found to increase the relative difference in mean generation time for all cases. In contrast, the difference in effective delayed neutron fraction was slightly reduced by the application of SPH in the adjoint (in cases where Griffin did not perfectly reproduce it).

Observations suggest that application of SPH while neglecting equivalence procedures in the steady-state adjoint step may provide significant improvement to the kinetics parameter preservation over the alternative approach. Additionally, the sensitivity of kinetics parameter results to regional application on either the full Cartesian region or fuel region only was reduced with the refinement of the energy grid. This further elevates the importance of gaining a better understanding and mitigating of the energy grid sensitivity.

While cases with SPH in either the full Cartesian region or fuel were capable of reproducing Serpent k_eff values, cases with no SPH equivalence procedure obtained significantly different k_eff values as listed in Table 11. These k_eff values are similar in CAS2 and C4G energy grids but different from the CAS4 and CAS8 energy grids, which are similar to one another. This difference is qualitatively similar to that observed in the kinetics parameters assessment and may arise from the same distortion of reaction rates. Note also that in all the Griffin calculations performed in this work, k_eff are identical between forward and adjoint eigenvalue calculations.

The large biases in k_eff without SPH are due to homogenization errors, as well as remaining discretization errors, which are inherently captured by the SPH equivalence procedure in Griffin. Omitting this step leads to strong differences in reaction rates, including fission. Despite the large differences in eigenvalue obtained by test cases without SPH, kinetics parameter results are similar to those obtained with SPH. It can be intuitively observed that these quantities are by definition bilinear ratios, where the weighting functions (forward and adjoint flux) appear at both the numerator and denominator, and thereby, errors therein tend to cancel out and result in similar kinetics parameter estimates even without SPH.

5. Supplementary test cases

Due to the complexity of the Griffin computational scheme, especially when the SPH equivalence procedure is applied on a subset of the whole geometry, simpler critical experiment benchmarks were developed. These tests were motivated by the desire to test SPH application over the full geometry as well as carry out a limited testing of SPH together with the S_N transport solver available in Grffin, which were impractical to perform on the IPEN/MB-01 model. The Godiva and Flattop-23 benchmarks available through the International Criticality Safety Benchmark Evaluation Project (ICSBEP) were modeled using both Serpent and Griffin. Both Godiva and Flattop-23 have a fast neutron spectrum, so different energy grids from the ones used for IPEN/MB-01 are utilized for multigroup cross section generation. The energy grid boundaries used are described in Table 12 and derived by collapsing the 33-group European Cell Code (ECCO) energy grid. This collapsing is necessary when generating group constants with a Monte Carlo code, as low flux tallies induce significant errors on the multigroup cross sections [24].

Utilizing the same Serpent−Griffin computational scheme as for IPEN/MB-01 test cases, diffusion provided poor estimates of kinetics parameters, particularly for the mean generation time Λ. In this instance, the diffusion calculation with an SPH equivalence fully reproduces the reference k_eff from Serpent and β_eff to within a few percent, yet fails with the mean generation time for which there were relative errors up to 15%. As a result, additional tests were performed utilizing Griffin’s S_N transport solver based on the self-adjoint angular flux formulation with a continuous finite element method (SAAF-CFEM-SN) [1] in the otherwise identical computation scheme. The S_N tests utilize a Gauss–Chebyshev angular quadrature with eight polar angles and are both angularly and spatially converged. Although the SPH equivalence procedure is usually targeted for a coarse-mesh diffusion calculation where homogenization, and in the case of Griffin where the cross sections are generated using a 3D model and spatial/angular/energy discretization errors remain, it is theoretically possible to apply an SPH correction with transport as well. Thus, S_N transport calculations were also employed in conjunction with an SPH equivalence, with an additional limitation being that the PJFNK SPH convergence can be much more difficult with S_N than with diffusion, especially with a large number of energy groups. Further investigation is warranted on this topic; however, the Griffin S_N transport calculation should converge towards the reference Serpent results with increased spatial/angular/energy refinement, regardless of the application of SPH procedure.

5.1. Godiva

The Godiva model is developed according to the bare Godiva sphere model specifications found in the ICSBEP handbook report HEU-MET-FAST-001 [25]. The reference measured delayed neutron fraction and mean generation time for the Godiva experiment are shown in Table 13 [26].

Serpent was run with one million particles in 1,000 active and 100 inactive cycles to evaluate kinetics parameters and to generate multigroup cross sections and equivalence data for downstream Griffin calculations. Results from Serpent Meulekamp, Nauchi, and IFP using the ENDF/B-VII.1 neutron cross-section libraries are given in Table 14 and compared to the benchmark values in Table 15. For these tests, the Nauchi results stand out for providing the closest estimate of the benchmark kinetics parameters, and corresponding delayed neutron data was utilized for Griffin calculations.

The Griffin model was developed as a one-dimensional, r-spherical geometry with 20 nodes of equal size in the radial direction. SPH equivalence was based on a single Serpent flux detector over the full sphere, and options of applying equivalence to both direct and adjoint flux evaluations, applying to adjoint flux only, and no equivalence were tested. Further evaluations were made to compare between the cross sections generated from Serpent using the ENDF/B-VII.1 library, between multigroup energy grids ECCO5 and ECCO24 and between Griffin diffusion and Griffin S_N methods. Griffin-evaluated kinetics parameters are compared to corresponding Serpent Nauchi results in Table 16 and to benchmark kinetics parameters in Table 17. Note that the SPH equivalence calculation did not converge using the S_N method with the ECCO24 energy grid, so results for these cases are not included.

Only a small change is observed between Griffin diffusion results using energy grids ECCO5 and ECCO24. While S_N results are limited for ECCO24, a significant reduction in relative difference in Λ is observed between the energy grids for the cases without an SPHequivalence.

5.2. Flattop-23

The Flattop-23 model is developed according to the specifications found in the ICSBEP handbook report U233-MET-FAST-006 [27]. The reference measured delayed neutron fraction and mean generation time for the Flattop-23 experiment are shown in Table 18 [26].

Serpent was run with one million particles in 1,000 active and 100 inactive cycles to evaluate kinetics parameters and to generate multigroup cross sections and equivalence data for downstream Griffin calculations. Results from Serpent Meulekamp, Nauchi, and IFP using the ENDF/B-VII.1 neutron cross-section library are given in Table 19 and compared to the benchmark values in Table 20. Delayed neutron data from Nauchi results were extracted for use in Griffin.

The Griffin model was developed as a one-dimensional, r-spherical geometry with 16 nodes of equal size in the radial direction representing the inner uranium-233 sphere and 80 nodes of equal radial size representing the outer uranium reflector shell. SPH equivalence was based on two Serpent detectors, one for each material region. Options of applying equivalence to both direct and adjoint flux evaluations, applying to adjoint flux only, and no equivalence were tested for both multigroup energy grids ECCO5 and ECCO24, and for Griffin diffusion and Griffin S_N methods. Griffin-evaluated kinetics parameters are compared to corresponding Serpent Nauchi results in Table 21 and to benchmark kinetics parameters in Table 22. Note that the SPH equivalence calculation did not converge using the S_N method with the ECCO24 energy grid, so results for these cases are not included.

Griffin diffusion results for the Flattop-23 model contain significant inconsistencies between SPH application strategies and energy grids utilized. Fortunately, as with the Godiva evaluations, utilizing a S_N transport method provides a significantly closer replication of Serpent results and shows significant improvement in the estimate of Λ with ECCO24 compared to ECCO5 for the case without SPH equivalence. These two critical benchmarks demonstrate that when transport calculations are performed with Griffin, the agreement between the Serpent and Griffin mean generation time improve with a finer energy group structure. Interestingly, diffusion with SPH equivalence fails at reproducing the mean generation time, even with the same energy group structure ECCO24, which showed satisfactory agreement with the S_N calculation. Thus, analysts need to be cognisant of these limitations, especially when performing transient calculations with diffusion.

6. Future work and conclusions

The validation of the Serpent model against the IPEN/MB-01 reactor kinetics benchmark produced kinetics parameter results similar to other evaluations with the ENDF/B-VII.1 cross-section library. By looking at the results obtained on the three reactor benchmarks studied in this work (IPEN/MB-01, Godiva, and Flattop-23), the Serpent result obtained with the IFP method provided the closest effective delayed neutron fraction β_eff against the measured one on average, and the Nauchi method provided a slightly better estimate of the mean generation time Λ than the IFP method. Regarding the reactivity measurements, similar trends are observed with Serpent as with the MCNP results reported in the benchmark. The largest uncertainties obtained for negative periods were identified in the MCNP evaluation of the benchmark and can be mitigated by utilizing an experimental decay constant for the first precursor group λ₁ rather than the one provided in the ENDF/B libraries. This conclusion is similarly applicable in this case.

Verification of Griffin’s preservation of Serpent kinetics parameters highlighted several considerations that may be useful in the development of future studies utilizing the Serpent–Griffin scheme. A comparison of SPH application options suggests that the SPH factors obtained by preservation of the forward Serpent fluxes may not be directly applicable to the adjoint k-eigenvalue calculation, which provides the adjoint flux used for the bilinear collapsing of kinetics parameters. The observation made in previous studies that the kinetics parameters are preserved better when SPH is applied only in the fuel region is further supported by this study. This study also identified the sensitivity of kinetics parameters to the selection of the multigroup energy grid as an area for potential further analysis. However, the additional validation cases consisting of the Godiva and Flattop experiments demonstrate that, when a transport calculation is performed with Griffin, the agreement in kinetics parameters between Serpent and Griffin increases when the energy group structure is refined, while the gain is only marginal for SPH diffusion. The same conclusion can also be reached observing the trend of the diffusion results without SPH equivalence for IPEN/MB-01.

While the observed qualitative shift between Griffin and Serpent kinetics parameters may be attributed to the misalignment of energy spectra in the model implemented in Griffin for the collapsing of effective delayed neutron fraction, the scale of the introduced differences may necessitate the development of a method for energy grid optimization. Further work might also be relevant within the Serpent homogenization approach in order to obtain delayed fission spectra for all the delayed precursor groups, instead of assuming them to be equal as it is likely a poor assumption, especially when many energy groups are considered. If individual spectra were made available in Serpent results, the preservation of kinetics parameters within Griffin may be further assessed within this scope.

The validation of the combined Serpent–Griffin calculation with the IPEN/MB-01 reactor kinetics benchmark provided similar conclusions as the individual step comparisons. However, the potential competing over- and under-prediction of kinetics parameters in the two steps may reduce the observed difference between Griffin and benchmark kinetics parameters. This error compensation is acknowledged in this study by separation of the analysis into two steps.

The following points summarize the key findings pertaining to the use of the Serpent–Griffin computational scheme and the capability of Griffin to accurately reproduce Serpent kinetics parameters:

• Griffin is capable of reproducing the Serpent effective delayed neutron fraction where the prompt and delayed fission spectra originating from the Serpent evaluation are identical. This was the case for CAS2 and C4G energy grids in the IPEN/MB-01 benchmark evaluation, and discussed in detail in Section 4.1. In some cases such as the CAS4 and CAS8 energy grid IPEN/MB-01 evaluations, the fission spectra will differ and may result in inconsistency between Griffin and Serpent kinetics parameters.

• Application of SPH equivalence procedures provided significant improvement to kinetics parameter estimates in the tested cases. For the IPEN/MB-01 benchmark evaluations, Griffin most closely reproduced Serpent kinetics parameters in cases where SPH was applied in the fuel region only, and excluded from the steady-state adjoint flux calculation.

• Test cases utilizing diffusion, with or without SPH equivalence, could lead to poor estimates of the mean generation time, particularly in the supplementary Godiva and Flattop test cases where the diffusion approximation is inappropriate due to the high anisotropy of the neutron flux. Utilizing Griffin’s S_N transport method was shown to rectify this issue for the supplementary test cases, however utilizing S_N for large models or highly refined energy grids may come at significant computational expense and difficulty in converging the SPH equivalence calculation.

• Refinement of the energy grid for the S_N method yielded significant improvement in the evaluation of the mean generation time, as demonstrated by the supplementary Godiva and Flattop test cases.

Finally, future validation work should include comparisons with space-time kinetics transients, such as the SPERT experiments or similar where a coupled multiphysics model is required. For cores with different isotopic compositions per fuel assembly zone, due to burnup distributions for instance, the impact of utilizing one set of kinetics parameters for the whole core versus one per homogenized fuel region should be investigated, as well as introducing in Serpent a bilinear weighting of the neutron velocities.

Conflict of interests

The authors declare that they have no competing interests to report.

Acknowledgments

This research made use of the resources of the High Performance Computing Center at Idaho National Laboratory, which is supported by the Office of Nuclear Energy of the U.S. Department of Energy and the Nuclear Science User Facilities under Contract No. DE-AC07-05ID14517.

Funding

This material is based upon work supported under an Integrated University Program Graduate Fellowship and has been co-authored by Battelle Energy Alliance, LLC under Contract No. DE-AC07-05ID14517 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.

Data availability statement

The input files for Griffin and Serpent developed as part of this work can be obtained by directly contacting the authors.

Author contribution statement

Cole Takasugi: formal analysis, methodology, data curation, writing – original draft. Nicolas Martin: conceptualization, methodology, writing – original draft. Vincent Labouré: writing – review and editing. Javier Ortensi: writing – review and editing. Kostadin Ivanov: writing – review and editing, funding acquisition, Maria Avramova: writing – review and editing, funding acquisition.

References

All Tables

All Figures

	Fig. 1. Serpent model geometry. (a) Top view of core region; (b) side view of core region.
In the text

	Fig. 2. Top-view cross section of Griffin mesh.
In the text