Underestimation of nuclear fuel burnup – theory, demonstration and solution in numerical models

Monte Carlo methodology provides reference statistical solution of neutron transport criticality problems of nuclear systems. Estimated reaction rates can be applied as an input to Bateman equations that govern isotopic evolution of reactor materials. Because statistical solution of Boltzmann equation is computationally expensive, it is in practice applied to time steps of limited length. In this paper we show that simple staircase step model leads to underprediction of numerical fuel burnup (Fissions per Initial Metal Atom FIMA). Theoretical considerations indicates that this error is inversely proportional to the length of the time step and origins from the variation of heating per source neutron. The bias can be diminished by application of predictor-corrector step model. A set of burnup simulations with various step length and coupling schemes has been performed. SERPENT code version 1.17 has been applied to the model of a typical fuel assembly from Pressurized Water Reactor. In reference case FIMA reaches 6.24% that is equivalent to about 60 GWD/tHM of industrial burnup. The discrepancies up to 1% have been observed depending on time step model and theoretical predictions are consistent with numerical results. Conclusions presented in this paper are important for research and development concerning nuclear fuel cycle also in the context of Gen4 systems.


Introduction
Continuous energy Monte Carlo burnup codes such as MCB [1] and SERPENT [2] are precise tools for depletion calculations in nuclear critical and subcritical systems.The biggest advantages of this approach is exact treatment of neutron transport physics in arbitrary geometries.This leads to reliable estimation of reaction rates that are applied in transmutation procedure.Solutions obtained by Monte Carlo burnup codes are considered as reference ones and the methodology is widely applied in fuel cycle studies.
In this paper we present a numerical problem that may exhibit as underestimation of numerical burnup produced by the code.The bias can be observed in all kinds of burnup codes, also deterministic ones.As will be shown, the observable underestimation of nuclear fuel FIMA (Fissions per Initial Metal Atom) can be present when long time steps are applied and the simplest coupling scheme is in use.The problem is independent from the nuclear data applied in neutron transport and numerical treatment of the Bateman equations.Problem is relevant for comparative studies of fuel cycles, benchmark calculations and experimental validations, in which numerical FIMA may be underestimated in comparison to the burnup in GWd/tHM imposed by code user or measured by the operator of industrial facility.
In Section 2 the methodology of burnup calculation is presented.Section 3 shows origins of numerical bias.Model applied in numerical tests is described in Section 4 and the results are analyzed in Section 5. Section 6 summarizes the research and provides conclusions for the good practices in fuel cycle studies.

Methodology
The approach applied in large variety of currently existing Monte Carlo burnup codes is explained in detail in the paper [3].The key aspect of this methodology is a statistical solution of neutron transport Boltzmann equation obtained at each point of declared time mesh.Monte Carlo code run provides physical quantities (so called tallies) such as reaction rates of nuclear reactions (ℝ) and heating per source neutron (h) that are corresponding to the critical configuration of nuclear system.Nuclear fuel isotopic composition ℕ at given point ̅ and time t is driven by Bateman Equations presented below: Transmutation matrix appearing in Bateman equations is a product of reaction rates and neutron source intensity plus decay matrix .The formula is presented below: It can be deduced that neutron source intensity describes dynamics of fuel depletion and stands for a normalization factor in transmutation problems.can be obtained using power normalization equation, which can be formulated in the following way: where : P is total fission power constrained in the system model.
The set of Bateman equations can be solved using a variety of numerical methods optimized for this purpose and several of them are discussed in the work by [4].For example MCB applies Transmutation Trajectory Analysis (TTA) and Chebyshev Rational Approximation Method (CRAM) is implemented in SERPENT.Formal mathematical solution of the Bateman equations is called matrix exponent presented below: where : ℕ is isotopic vector at the beginning of step.
It should be noted that coefficients of transmutation matrix are assumed constant and on the other hand ℝ and can in general vary significantly over time in realistic nuclear systems.The compromise commonly applied in burnup codes is a limited length of burnup steps, during which the variation of reactor model properties is neglected.The most basic step model is called Euler predictor (aka.staircase), however predictor-corrector methodology becomes more commonly applied in currently developed codes [5].In the next section it will be shown that if variation of transmutation coefficients is not taken into account, the results suffer from inherent numerical bias.

Bias of numerical burnup
FIMA indicator of nuclear fuel burnup can defined as a ratio between fissioned heavy metal atoms (∆ ) and their initial number ( , ) in fresh fuel as shown below: In practice, it is easier to measure the amount energy released in industrial reactor per initial heavy metal mass ( ) and equivalent formulation of burnup ( ) is often in use: As long as the energy released per fission in system is approximately constant (roughly 200MeV), it can be assumed that both formulations are physically equivalent.
Most often the users of burnup codes declare the burnup of fuel using its second formulation (6) and the staircase function.It appears that numerical solution of Bateman equations does not automatically take into account that heating per source neutron h decreases over time step.This should normally result as an increase of neutron source intensity S, however in the simplest case only the beginning of step value is applied.The problem is schematically shown in Fig. 1 below.Neutron source intensity governs the speed of nuclear reactions in model during each step and it can be assumed that the number of fissions is directly proportional to the area under the curve in Fig. 1.
The error of the fuel burnup at each time step is proportional to the area between the real curve and staircase function (green color in Fig. 1).The considered bias of FIMA indicator over irradiation period T is defined by the ratio of integrals: As will be shown, the variation of S in numerical models can be well represented by linear function, thus the absolute or relative bias can be estimated analytically.The error can be described as the sum of n triangular surfaces representing contributions from each time step: This prediction has prominent implications for the accuracy of the Euler predictor step model.The underestimation of FIMA is dependent on time grid applied in numerical model.The bias is directly proportional to length of time step.Logarithm can be applied to both sides of Eq. ( 8) so as to provide the form more useful in comparative analyses: where: = • ∆ /

∆
. Under such formulation the analyzed error should vary linearly with logarithm of applied time step.This model will be verified in the following sections.More advanced coupling schemes such as predictorcorrector variants (described by [6]) takes into account the variation of neutron source intensity and reaction rates.The simplest solution is application of an average values from the beginning and end of each step as an effective input for transmutation coefficients.In such a way the error originating from variation of S should disappear.Simultaneously this approach requires at least one more Monte Carlo neutron transport run at each time step, which is a computational drawback.
The bias of FIMA described in this section is only one among several sources of systematic errors possibly present in burnup simulations.More examples of numerical and modeling difficulties and corresponding solutions are presented in the papers by [7] and [8].

Test model of PWR assembly
In order to demonstrate the presence of the problem explained in previous section, the model of typical fuel assembly from Pressurized Water Reactor (PWR) has been prepared.It consists of 17x17 lattice of fuel pins among which 25 control rod guide tubes are localized.Uranium dioxide is surrounded by zirconium cladding and light water coolant.The section of this 2-dimensional model is shown in Fig. 2 below.The geometry and compositions applied in this model are presented in detail in the article by [7] and also very similar model has been applied in the study by [6].Both descriptions are sufficient for reproducing our results and the model is representative for commercially applied PWR reactors.The only modification applied in the current research is the enrichment of uranium, which has been increased up to 5%.
The burnup studies have been performed using the SERPENT code version 1.18, which belongs to family of continuous energy Monte Carlo burnup codes.Bateman equation solver has been set to CRAM method and neutron transport precision has been set to 40 pcm of standard deviation on keff.Cross-section libraries applied to neutron transport is JEFF3.1 and the rest of nuclear data (decay constants, fission yields, branching ratios) have been provided by the code author [9].Power in the system has been normalized to 14MWth per assembly and various number of time steps during irradiation period of 2000 days has been analyzed.In each case we compared staircase and predictor-corrector time step models.

Computational results
In this chapter we present computational results concerning the output FIMA in numerical model.Firstly, the fuel burnup at the End-of-Irradiation for considered step models and lengths of step (between 5 and 2000 days) is shown in Fig. 3 below.It is clearly visible that both the applied step length and the coupling scheme have systematic influence on the output burnup of fuel.We analyze here FIMA indicator only, however its discrepancies indicate differences for all the isotopic concentration separately as well as.For staircase model the bias quickly increases with time step length reaching 1% in the extreme case.On the other case, predictor-corrector model gives roughly same FIMA regardless from step length applied.
Next, the most exact result (predictor-corrector, 5days step) is treated as the reference solution and the absolute discrepancies from its FIMA are studied.The observed differences are presented in Fig. 4. The error of results obtained using staircase model varies linearly as predicted in Section 3. In case of predictor-corrector the numerical bias is always over one order of magnitude lower.For smaller lengths of steps we observe saturated behavior of error due to the limited numerical precision of variables in output files.Finally, the evolution of the neutron source intensity and the heating per source neutron (obtained with identical model with MCB code) is exposed in Fig. 5.We deduce from the figure that an assumption about linear evolution of S is reasonable and its variation is nonnegligible in practical applications.Obtained linear trend can be applied to compare theoretical prediction (Eq.8) and the value of error from numerical tests.We performed calculation for the case of 2000-day time step.Analytical approach gives the bias value of -1.7% while the numerical test shows FIMA underestimation of about -1.1%.The values are not exactly equal as the other sources of error exist in the numerical model and can partially compensate each other.In this case we attribute the difference to spectral changes of neutron flux, which influence reaction rates in the system and partly compensate bias caused by the variation of S. Observed bias is meaningful in comparison to the FIMA at the end of cycle (~6.24%).

Conclusions and summary
The problem of numerical underestimation of Fissions per Initial Metal Atoms in burnup codes has been studied theoretically and numerically.Analytical considerations indicate that step model assuming beginning-of-step coefficients gives systematic error in predicted fuel burnup.The bias value is directly proportional to the length of time step of numerical model.The problem originates from the variation of heating per source neutron in a system, which modifies transmutation coefficients present in the Bateman equations.
Numerical tests performed with SERPENT code confirmed presence of the problem for the most basic step model.The second applied coupling scheme -predictorcorrector -provided more reliable estimation of fuel burnup for any step length.We conclude that advanced step models eliminate analyzed bias in first approximation.Such solution can be recommended for all upcoming burnup simulations based on computationally expensive Monte Carlo neutron transport method.Usage of predictor-corrector coupling scheme provides more exact results than simple reduction of time step length while decreasing or maintaining overall computational cost.
Analyzed type of error is not specific to Monte Carlo burnup methodology, but can be present also in results of deterministic codes.It belongs to systematic errors present in modeling of nuclear systems and awareness about it inevitable to provide reliable results.Fortunately, this type of bias can be effectively excluded/diminished in the codes currently available for nuclear reactor engineers.The neglect of considered bias can arise following problems: 1) overestimation of energy extracted from the fuel, 2) underestimation of the heavy metal consumption, 3) propagated errors of burnup-dependent temperature coefficients (in case of safety assessments).
The results and analysis are meaningful basically for all types of calculation performed in research and development of nuclear fuel cycle including: validation of codes using experimental measurements, code-to-code comparison, benchmark problems, comparison of various fuel cycles scenarios, nuclear wastes radioactivity and transmutation assessments.In fact, all types of safety assessments account for burnup feedback.Reliable calculations in this field are fundamental for upcoming concepts of Gen-III+ and Gen-IV reactors.

Figure 1 .
Figure 1.Scheme explaining the source of possible numerical burnup underestimation in depletion simulations.

Figure 2 .
Figure 2. Horizontal section of geometry applied in model.Fuel pins are navy blue, guide tubes are red and coolant/moderator is yellow.

Figure 3 .
Figure 3. Nuclear fuel FIMA at the End-of-Irradiation (2000 days) for various time discretization and coupling schemes.

Figure 4 .
Figure 4. Convergence of the output FIMA values for two considered step models.

Figure 5 .
Figure 5. Variation of heating per source neutron during irradiation and resulting change of neutron source intensity.