© G. Billiet et al., Published by EDP Sciences, 2024

1. Introduction

The French energy and electricity mix is undergoing significant transformations. The phase-out of fossil fuels, the electrification, and the closure of historical Pressurised Water Reactors (PWRs), built in the seventies, pave the way for deploying decarbonized energies such as renewables and new nuclear technologies. Consequently, the structure of the nuclear fleet, and its constraints regarding load following are changing. France is unique as its strategy is to recycle all plutonium from spent fuel into MOX fuels, to stabilize the UOX spent fuel stockpile. The industrial implementation of the downstream nuclear fuel cycle decreases the consumption of natural resources and enhances energy independence. It also prepares for the future deployment of closed fuel cycles and plutonium multi-recycling in fast neutron spectrum reactors such as Sodium Fast Reactors for instance.

Presently, the French nuclear fuel cycle allows for the introduction of three fuel types in the reactors. Irradiated uranium from Enriched Natural Uranium (ENU) fuel, is recycled in Enriched Reprocessed Uranium (ERU) fuel. Plutonium obtained from the irradiation of ENU fuel is also recycled using a blend of depleted uranium oxides and plutonium (MOX fuel).

The study of the whole fleet in evolution is possible thanks to dynamical fuel cycle simulations or scenario studies. They quantify the impacts of any evolution of the fuel reprocessing strategy by simulating all the units in the whole fleet, such as reactors but also fuel fabrication units, reprocessing units, cooling pools, and geological waste disposals. Such simulations provide a way to model an entire nuclear fleet accurately enough to track isotopic inventories’ evolution. The difficulties in such studies lie in reactor modeling as soon as fissile material reprocessing is involved. Indeed, in such cases, the fresh fuel compositions are not known a priori and so is the composition at each reactor unloading. Similarly, the achievable burn-up of different fresh fuels depends also on its composition. The CLASS code [1], developed at the Centre National de la Recherche Scientifique (CNRS), is such a nuclear reactor fleet simulator, that is based on precise reactor models that allow, among other, plutonium reprocessing studies. As this code does not solve neutronics equations (Boltzmann nor Bateman’s équations) directly, it may be linked with depletion codes, such as, for example, DRAGON and DONJON [2, 3], developed at Polytechnique Montréal. Direct couplings of DRAGON and DONJON with CLASS were previously conceived for PWRs, leading to a fresh fuel fabrication model, and a material irradiation model. The first model aims to calculate the plutonium content in fresh MOX fuels, according to the targeted burn-up and the isotopic composition of the plutonium available in different stocks. The second one links the fresh fuel composition and the reactor burn-up, to the isotopy at reactor discharge.

DRAGON and DONJON are themselves linked in order to perform 3D full-core simulations. DRAGON solves the neutronics transport equation on a 2D infinite lattice of similar fuel assemblies. It then exports lattice averaged macroscopic cross-section and diffusion coefficient (i.e. diffusion data) on an object called Multicompo (MCPO), a library of neutronics data tabulated as a function of different parameters. DONJON uses these MCPO, interpolating at the right parameters for each spatial node (in each fuel assembly), to solve the diffusion equation on a 3D full-core reactor model.

Initial CLASS reactor models were based on infinite assembly depletion calculations, for instance, made by DRAGON 2D lattice code, without modeling the full-scale reactor [4]. Then, DONJON 3D core simulations were developed and linked to CLASS, for homogeneous PWRs (900 MWe), containing only one fuel type [4, 5]. Heterogeneous DONJON simulations were then implemented, allowing for correct modeling of UOX/MOX reactors[6–8]. While considerably improving the accuracy of fleet simulations, those coupling methods come with a high computation time.

In the meantime, neutronics thermal-hydraulics couplings for PWRs were added to DRAGON and DONJON codes [9], but remained unused for CLASS applications, where an isothermal model with a uniform fuel temperature map (equal to 900 K) was used in previous works. Such a hypothesis could create significant biases, due to the Doppler effect, which affects neutronics data. Moreover, the radial power distribution is incorrectly simulated in isothermal simulations, and a uniform fuel temperature may overestimate the power emitted by normally hot assemblies, and underestimate the power produced by colder assemblies, resulting, for example, in an overestimation of the power factor¹.

This work presents a new coupling scheme between CLASS, DONJON, and DRAGON that would allow neutronics thermal-hydraulics coupled 3D core calculations in order to quantify numerical biases induced by the isothermal hypothesis in dynamical fuel cycle simulations. To do so, a new interpolation method for diffusion data estimation for the 3D neutron flux distribution at the core level has been developed in order to greatly accelerate core calculation. Without this improvement, the core calculation with DONJON would have been too costly to be included in a fuel cycle simulation where hundreds of depletion simulations are needed. This article focuses only on ENU and MOX fuels, but the methodology can be easily extended to other current or future types of fuels.

This paper first presents this new interpolation method, and its improved performances regarding numerical cost and accuracy, when applied to DRAGON-DONJON simulations used by CLASS. Then, thanks to this improved efficiency, the paper introduces multi-physics couplings and thermohydraulics equations in DRAGON and DONJON models used by CLASS. Biases of the isothermal hypothesis on spent fuel isotopic composition are quantified, and a discussion is proposed on the relevance of neutronics thermal-hydraulics coupled 3D simulations for nuclear scenario studies. Finally, those biases are propagated on a simple academic scenario and on a more realistic one, inspired by the French institutional scenario, in order to quantify errors induced by reactor modeling in dynamical fuel cycle simulations.

2. Lattice calculation, recombination and interpolation

DONJON’s high computation times are mainly caused by the interpolation of neutronics data, produced by DRAGON for a sampling batch, and stored in a lattice-core interface object, named MCPO [10]. When modeling nuclear reactor fleets with CLASS, this MCPO shall contain neutronics data describing the whole fleet, whereas in reactor neutronics codes, it usually only describes fuel assembly behavior for one reactor cycle. Thus, the MCPOs used by CLASS are heavier, because they contain new data dimensions, such as the initial isotopic compositions, sampled in order to describe the variety of fuel assemblies that may be loaded at the beginning of each reactor cycle.

Fig. 1. 2-level calculation scheme principle for 3D core simulation.

In this work, there is one MCPO per fuel type (ENU and MOX, ERU is not implemented). Fresh ENU isotopic composition, containing ²³⁵U, and ²³⁸U, can simply be described by using uranium enrichment (the total mass being determined by the reactor power). But as soon as reprocessed fissile material (ERU and MOX) is loaded, the fresh isotopic composition becomes more complex due to the presence of multiple new isotopes (plutonium, other uranium isotopes, etc). Because neutronics data produced by DRAGON depletion calculations strongly depend on those fresh isotopies, the initial idea in previous CLASS and DONJON coupling was to add one dimension per isotope (e.g. %²³⁸Pu…%²⁴²Pu for MOX fuel) to the MCPO. As a result, the parameterized space, and its sampling size, are greatly increased leading to huge MCPOs and high computation times during the DONJON interpolation process (more than 80% of the total DONJON computational time). It should be pointed out that the uranium of MOX fuel is depleted curation, hence, the proportion of ²³⁵U is kept constant at 0.25%. While heavy, those calculations were still possible, but the thermal-hydraulics coupling would further increase the number of parameters. Consequently, in order to implement such coupled calculations, the computation time of the core calculations and the structure of MCPOs have to be improved.

2.1. Diffusion cross sections interpolation principles

3D modeling of power reactors, such as 900 MW_e PWR, often requires a 2-level calculation scheme. In this paper, the first level (DRAGON) aims to solve the precise neutron transport equation on a simplified 2D geometry, and the Bateman evolution equation, in order to produce 2 groups of burn-up dependant diffusion data along reactor cycles. The second level (DONJON) is dedicated to the core level where the diffusion equation is solved at each time step to evaluate the neutron flux and the power distributions. For each time step, and for each spatial mesh, the core solver, DONJON in this case, has to compute all the diffusion data, regarding the local characterization, by interpolating the MCPO. Different interpolation modules in DONJON are implementing a multi-dimensional linear or cubic interpolation scheme between neighbor tabulated points. The principle of such a 2-level scheme is presented in Figure 1.

Consequently, parameters that may vary at the core level, and which have significant neutronics effects, should be considered at the lattice level. Boron concentration in the core, mean burn-up over an assembly, and effective fuel temperature (added for thermal-hydraulics purposes) are dimensions that do not describe the fuel isotopy but which are implemented, in order to simulate the variety of conditions under which the fuel is irradiated. They are often (or even always) computed in all-purpose neutronics simulations. Here, we use the burn-up as a time variable; the boron concentration is adjusted at each time step in order to maintain neutronics criticity; fuel temperatures are added in order to take into account the Doppler effect (see Sect. 3). Additionally, for fuel cycle simulations, initial loaded isotopies also have to be considered as parameters, as they are solved step-by-step during the scenario simulations (but fixed during one reactor simulation). For UOX fuels, one isotopic dimension needs to be added to the MCPO, representing the ²³⁵U enrichment. In the case of MOX fuels, six isotopic dimensions should be considered for the MCPO: the total loaded plutonium content %_Pu in MOX, and the loaded plutonium relative isotopic vector $\overrightarrow{\mathrm{Pu}}=({\%}_{{}^{238}\mathrm{Pu}},{\%}_{{}^{240}\mathrm{Pu}},{\%}_{{}^{241}\mathrm{Pu}},{\%}_{{}^{242}\mathrm{Pu}},{\%}_{{}^{241}\mathrm{Am}})$, in % of total Pu-Am contents. ²³⁹Pu contents in MOX is used as a buffer and can be determined using $\overrightarrow{\mathrm{Pu}}$ and %_Pu.

The definition of the MCPO parameters, their range, and the way they are sampled are of prime importance for DRAGON lattice calculations. A simple but heavy way, however quite precise, would be to perform one specific depletion simulation for each possible set of parameters, as it has historically been done for the DONJON-CLASS coupling [4]. In this paper, such MCPOs are called regular MCPOs. In that case, starting from a reference sampling point, all cross terms, describing simultaneous variations of different isotope initial quantities, are taken into account in the interpolation process. The main consequence is obviously the high computational costs this approach requires. Indeed, the introduction of plutonium isotopic content as new parameters for lattice simulations leads to some MCPOs nearly 1000 times heavier, with only 3 sampling points for each of the six additional dimensions that describe the fresh fuel composition. Supposing that DONJON uses, for each sampling segment, a multidimensional linear interpolation method of the data contained in the MCPO², and supposing that there are only 3 parameters (named hereafter x, y, and z) considered for lattice simulations, the standard interpolation process is proposed in equation 1 where interpolation is replaced by linear Taylor expansions around a tabulated point (presented only as a mathematical illustration):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Sigma }(x,y,z)=& \mathrm{\Sigma }(x_{\mathrm{ref}},y_{\mathrm{ref}},z_{\mathrm{ref}})\hfill \\ \hfill & +(x-x_{\mathrm{ref}})\frac{\partial \mathrm{\Sigma }}{\partial x}+(y-y_{\mathrm{ref}})\frac{\partial \mathrm{\Sigma }}{\partial y}\hfill \\ \hfill & +(z-z_{\mathrm{ref}})\frac{\partial \mathrm{\Sigma }}{\partial z}\hfill \\ \hfill & +(x-x_{\mathrm{ref}})(y-y_{\mathrm{ref}})\frac{{\partial }^2\mathrm{\Sigma }}{\partial x\partial y}\hfill \\ \hfill & +(x-x_{\mathrm{ref}})(z-z_{\mathrm{ref}})\frac{{\partial }^2\mathrm{\Sigma }}{\partial x\partial z}\hfill \\ \hfill & +(z-z_{\mathrm{ref}})(y-y_{\mathrm{ref}})\frac{{\partial }^2\mathrm{\Sigma }}{\partial z\partial y}\hfill \\ \hfill & +(x-x_{\mathrm{ref}})(y-y_{\mathrm{ref}})(z-z_{\mathrm{ref}})·\frac{{\partial }^3\mathrm{\Sigma }}{\partial x\partial y\partial z}\hfill \end{array}\end{array}$$(1)

where (x_ref, y_ref, z_ref) is the set of reference parameters, usually at the center of the sampling, and Σ a macroscopic cross-section (or diffusion coefficient) is required for the DONJON simulation. In DONJON, each partial derivative function is calculated thanks to the set of pre-calculated DRAGON simulations.

For all-purpose neutronics simulation, the standard way to proceed [11] is different: first, a specific depletion calculation is performed for a reference situation (fuel temperature, boron concentration, …). Then, at each burn-up step, several successive transport calculations are performed to cover the full list of perturbed parameters and generate perturbation MCPOs that can be used later for interpolation in DONJON. This option has to be implemented carefully, as it leads to biases due to the impact of parameters on the depletion simulation. The regular sampling of equation 1 can thus be lightened, preferably for the new isotopic dimensions, regarding computation times. Instead of using a regular grid sampling, containing all combinations of points for each dimension, a perturbative grid sampling can be created. A reference isotopic vector is defined, which represents the first sampling state of the MCPO. Then, this reference state is perturbed, and isotopic dimensions are sampled one by one: only one of those dimensions can have a sampling value different from the reference.

The resulting MCPO is called in the following, perturbed MCPO. With such a perturbative sampling, all the cross terms of equation 1 are impossible to compute and, hence neglected. With all dimensions considered as perturbed in the sampling process, this equation becomes expression 2, where $\overrightarrow{x}$ represents a vector of n parameters (instead of the previous three).

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Sigma }\left(\overrightarrow{x_t}\right)& =\mathrm{\Sigma }\left(\overrightarrow{x_{\mathrm{ref}}}\right)+\sum _{i=1}^n(x_t^i-x_{\mathrm{ref}}^i)\frac{\partial \mathrm{\Sigma }}{\partial x^i}{|}_{\begin{array}{c}\overrightarrow{x_{\mathrm{ref}}}\end{array}}\hfill \end{array}\end{array}$$(2)

where $\mathrm{\Sigma }(\overrightarrow{x})$ is the cross-section data contained in the MCPO for an isotopic vector $\overrightarrow{x}$ describing perturbative dimensions, where $\overrightarrow{x_{\mathrm{ref}}}=(x_{\mathrm{ref}}^1,...,x_{\mathrm{ref}}^n)$ is the isotopic reference for perturbative sampling, and where $\overrightarrow{x_t}=(x_t^1,...,x_t^n)$ is the target isotopic vector, i.e. the isotopic content loaded into the reactor.

Perturbative interpolation is faster but less accurate than regular interpolation. It is possible to consider an intermediate approach combining perturbed interpolation and regular interpolation. First, some non-isotopic dimensions such as burn-up and boron concentration always remain regular. For the UOX fuel, the isotopy being simply described by the uranium enrichment, the UOX MCPO remains fully regular. Then, for the MOX fuel, the %_Pu dimension has a strong influence on the neutronics behavior of fuel assembly. This dimension thus remains regular. Other MOX isotopic dimensions, i.e. $\overrightarrow{\mathrm{Pu}}$ dimensions, are perturbed. In that case, the perturbative grid sampling, corresponding to the $\overrightarrow{\mathrm{Pu}}$ sampling in the MCPO, exists for each value of %_Pu (and of the burn-up and boron concentration). Cross derivative terms between %_Pu and each Pu_i are considered in the interpolation process but not the cross terms between different elements of the $\overrightarrow{\mathrm{Pu}}$ vector.

For the purposes of this work, new MCPOs shall include variables representing the effective fuel temperature, averaged over the fuel assembly³. The overall purpose of the current section is optimizing computation times, this new dimension should be integrated inside the MOX and UOX MCPOs as a perturbative parameter. For the MOX MCPO, fuel temperature sampling thus exists for a single plutonium vector $\overrightarrow{{\mathrm{Pu}}_{\mathrm{ref}}}$, an isotopic reference vector for temperature interpolations⁴. This implementation of the fuel temperature means that the Doppler effect, previously based on a homogeneous fuel temperature map, will now simultaneously vary over fuel temperature, and over other parameters such as burn-ups, boron concentrations, %_Pu (and DONJON spatial mesh), but not over the initial plutonium isotopy $\overrightarrow{\mathrm{Pu}}$ at fixed total plutonium rate.

MCPO structures and dimension sampling types are summed in Table 1. For the UOX MCPO, the fuel temperature dimension is added as a regular parameter because there is no other perturbative dimension in the UOX MCPO: a dimension cannot be perturbed alone.

Implementing this sampling method for MOX fuel greatly reduces the number of isotopic sampling points. Our chosen sampling grid contains 7 points for %_Pu and 3 points for each $\overrightarrow{\mathrm{Pu}}$ dimension: thus, the total sampling size due to those isotopic dimensions is 7 × (1 + 2 × 5)=77 states, whereas the same resolution would have given a sampling size of 7 × 3⁵ = 1701 states with an all-regular sampling method. The computational time when using those optimized MCPOs is reduced by a factor of 22, even if the number of dimensions remains identical in both cases. This direct correlation between computation time and MCPO isotopic sampling size also proves that interpolation time accounts for the majority of DONJON computation time, with such parameterization. Small but significant biases have been observed, mainly directed by the distance between the reference isotopy $\overrightarrow{x^{\mathrm{ref}}}$ and the loaded isotopy $\overrightarrow{x}$. Those biases are displayed in Figure 2 and studied in Section 2.2, along with biases associated to a new recombination method, which eliminates most of the bias created.

Fig. 2. Biases between DONJON simulations using different sampling for MCPOs, on plutonium isotope inventories at reactor discharges. Biases are displayed, for each isotope, over the absolute value of the relative distance between loaded and reference isotopic value, for this same isotope. Biases are displayed for perturbative MCPO (in orange) and for recombined perturbative MCPO (in green), relative to a full MCPO with regular grid sampling.

2.2. Use of recombination for scenario simulations

The accuracy of core simulations using MCPOs with perturbative sampling depends greatly on the definition of the reference set of parameters and more specifically on the distance between the reference and the considered set of parameters. For this work, the initial plutonium isotopic vector has been considered as a perturbative dimension, but may vary during a fuel cycle simulation. In order to reduce this kind of bias, due to the choice of the reference isotopic vector, a code overlay is added, before each call to DONJON by CLASS. This overlay, as described in Figure 3, allows to create a local MCPO, from a large regular grid of MCPOs with better sampling resolution. This large grid contains respectively 6, 4, 4, 5, 5 and 3 points for %_Pu, %_²³⁸Pu, %_²⁴⁰Pu, %_²⁴¹Pu, %_²⁴²Pu and %_²⁴¹Am, for a total loaded isotopy size of 7200 states. For each set of those parameters, burn-up and boron concentration are considered as regular parameters.

Fig. 3. DONJON-CLASS chaining method, with recombination code overlay, used to produce efficient irradiation models.

The regular sampling of this large database allows us to choose the state whose isotopic values are closest to the loaded isotopy, as the reference isotopic vector. The same semi-perturbative sampling (with %_Pu as a regular dimension) is then constructed around this reference state, in order to produce a new dedicated MCPO for DONJON. This recombined MCPO contains only 3 points per fuel isotopic dimension, including %_Pu, for a total of 33 isotopic states, and is thus lighter than the previous perturbative MCPO due to the removal of unnecessary %_Pu sampling points. Moreover, the distance between the reference and the set of parameters of interest is always kept minimum. Thanks to the removal of the %_Pu dimension, and despite the initial call to a recombination DRAGON5-Python3 algorithm, the computation time is again improved by a factor of 2. The fuel temperature dimension is dealt with by simply adding the corresponding reference and sampled elementary MPCOs to the local database.

Biases on isotopic compositions at the end of irradiation calculated by DONJON, resulting from the use of those new MCPOs, are presented in Figure 2. PWR simulations loaded with 28, 28, and 10 different plutonium isotopic vectors are launched, using respectively a regular MCPO⁵, a perturbative MCPO, and a recombined perturbative MCPO. The regular sampling is considered a reference simulation as the interpolation process is the most complete (because it takes into account all cross-terms). Isotopic inventories at reactor discharge are compared to the reference calculations for the two other built MCPOs, and the figure represents the deviation for each plutonium isotope (defined as $\frac{({\mathrm{Pu}}_i-{\mathrm{Pu}}_i^{\mathrm{ref}})}{{\mathrm{Pu}}_i^{\mathrm{ref}}}$). Biases are displayed over the distance between the loaded plutonium isotopic concentration and the reference ones. End-of-irradiation biases are studied for isotopes that compose the initial isotopy. Biases are displayed, for each of these isotopes, over the absolute value of the relative distance between loaded and reference isotopic value, for this same isotope (for instance, bias on unloaded ²³⁸Pu is plotted over the deviation between loaded ²³⁸Pu and reference value for ²³⁸Pu). Observed biases using the perturbative MCPO, with regard to the regular one, stress the importance of choosing an accurate isotopic reference. Indeed, those biases increase when the loaded plutonium isotopy moves away from the reference plutonium isotopy.

Using a recombined perturbative sampling reduces significantly those isotopic biases to less than 0.2%, falling beyond the targeted accuracy of (0.1–1)% for reactor modeling in scenario simulations, and thus validating this method.

2.3. Core simulation upstream interpolation

During DONJON simulations, cross-section data interpolations are required before each flux evaluation, at each time step. An iterative algorithm is used to compute the 2-groups flux distribution to get the bore-critical state or to compute the Doppler effect (see Sect. 3). Even if MCPO dimensions corresponding to initial isotopic compositions are necessary to describe the whole fleet of reactors, once the fuel has been loaded by CLASS into the reactor in order to simulate irradiation by DONJON, those dimensions are always interpolated on the same target isotopy, which is the loaded one. Therefore, there is a redundancy for some of the interpolated dimensions. This redundancy can be treated by pre-interpolating those dimensions, using another code overlay, before CLASS starts calling DONJON, after the recombination code overlay (cf. Sect. 2.2).

Pre-interpolation can only be performed on dimensions that are fixed during irradiation, such as the initial composition for instance, but not the boron concentration, the burn-up, or the fuel temperatures.

Consequently, there are four types of dimensions to be treated by this code overlay: regular pre-interpolated dimensions (%_Pu for MOX, ²³⁵U enrichment for UOX), perturbative pre-interpolated dimensions ($\overrightarrow{\mathrm{Pu}}$ forMOX), other regular dimensions (burn-up, boron concentration), and other perturbative dimensions (fuel temperature). Cubic spline interpolation is performed on pre-interpolated dimensions, a method which slightly differs from DONJON segment cubic interpolation with Ceschino derivative estimation [12]. Using this interpolation, the database is flattened along regular pre-interpolated dimensions, and perturbative cross sections pre-interpolations, expressed in the right term of equation 2, are added as constants to all remaining data values.

The resulting database produces DONJON unloading isotopic biases of less than 10⁻²% and a 30% shorter computation time, with regard to a recombined but not pre-interpolated MCPO, which thus validates this second upstream code overlay.

3. Fuel depletion with fuel temperature distribution

The previous DONJON-CLASS chaining structure was supposed a uniform 900 K fuel temperature map for DONJON calculations. This hypothesis deprecates a well-documented phenomenon, known as the Doppler effect, on depletion simulations, that could also be responsible for significant bias on fleet-scaled codes like CLASS. The Doppler effect, which is defined as the cross-section dependence on temperature, leads to a negative reactivity coefficient in PWRs. Consequently, inside a fuel assembly, an increase in the fuel temperature causes a decrease in the thermal power, which is appreciated for stability reasons. Thus, when using an homogeneous fuel temperature map, the power of hotter assemblies is overestimated, and the power of the colder ones is underestimated, which could lead to an overestimation of the power factor, one of the main indicators used by the fuel fabrication models in CLASS. Moreover, as previously said, the fuel temperature affects the neutron cross-sections and the isothermal hypothesis of the core could also lead to biases of the irradiation models of CLASS. Hence, neutronics and thermal-hydraulics coupled 3D DONJON simulations are necessary for biases quantification at the core level, detailed below, and at the whole fleet level (detailed in Sect. 4).

Fig. 4. DONJON thermal-hydraulics neutronics computation scheme. The reactor is boron-critical, meaning that the boron concentration in the coolant is adjusted to maintain neutronics criticality and a controlled chain reaction (keff = 1). EOC stands for End Of Cycle (or irradiation).

3.1. Coupling principle and simulation issues

Thanks to efficiency improvements made in Section 2 regarding the interpolation process in DONJON simulations, coupled thermal-hydraulic and neutronics models can now be implemented. For this work, void coefficient distribution, due to local coolant density variations is not implemented: a uniform coolant density is used when interpolating cross-section data from the MCPO. The moderator temperature is also considered a constant. This choice has been driven by the need to study the effect of local fuel temperatures alone, for which biases are expected, on irradiation and on fuel fabrication models. Considering a uniform coolant density also helps simplify new MCPO structures and coupling methods. Nevertheless, the moderator temperature coefficient is also known to be systematically negative in PWRs, and taking into account local moderator density should modify the axial power distribution. The moderator temperature evolution in the core should not exceed 20 K, and its radial variation is very limited (smaller to 1 K), whereas the fuel temperature radial distribution may vary more than 100 K. We therefore assume that if the radial temperature distribution does not affect power redistribution, the impact of the moderator temperature distribution is negligible. As in CLASS the models focus on radial power distributions, this refinement does not seem necessary for this work, but local coolant densities should be implemented in future simulations⁶.

A thermal-hydraulic module named THM, calibrated on 900 MWe PWRs, has already been implemented in DONJON5 [9], and is capable of computing fuel temperature map and other thermal-hydraulic variables on the core mesh, based on the power distribution provided by neutronics. To reflect that thermal-hydraulics feedback on neutronics, MCPOs are sampled along thermal-hydraulics dimensions (see Sect. 2). Because local thermal-hydraulics originate from 3D core simulations, the corresponding dimension added to MCPOs is an effective fuel rod temperature, averaged over the lattice assembly. In DONJON, an iterative method, cycling flux calculations providing power maps using MCPOs, and thermal-hydraulics calculations providing fuel temperatures, is implemented at each thermal-hydraulics update.

Fig. 5. Histograms of unloaded isotopic relative differences for a variety of loading isotopic compositions, when adding thermal-hydraulics coupling, with regard to isothermal core models, plutonium isotopes, for UOX and MOX fuels, for PARITE and TIER fuel loading managements. The fuel loading management correspond to 4-rounds (PARITE) and 3-rounds (TIER) MOX and UOX irradiations, as described in [14]. Histograms are stacked with a width of 0.1% relative difference value. The references are the multi-physic simulations.

The new DONJON calculation scheme for depletion simulations is displayed in Figure 4. Thermal-hydraulics variables are updated at each time step, after computing boron critical concentration. The power map provided by the resulting boron-critical updated-thermal-hydraulics state is then used to perform the next burn-up step. Initial loading state stabilization consists of two successive boron-critical and thermal-hydraulics calculations, in order to obtain critical boron concentration from nonuniform thermal-hydraulic maps. More in-depth couplings, using boron-THM successive iterations at each time step and coupled convergence criteria, have not been implemented, due to the very high computation time, and a negligible bias.

The updated computation scheme coupling neutronics with thermal-hydraulics may lead to various oscillations of neutronics and thermal-hydraulics variables, as observed in preliminary simulations and in recent literature [13]. Those oscillations has been resolved by a high fuel temperature convergence criterion, small burn-up steps, and average temperatures and power depositions over two successive iterations in the temperature calculation loop, such as in [13].

3.2. Impact of fuel temperature distribution on depletion simulations

The updated computation scheme produces biases on unloaded isotopic inventories. Added thermal-hydraulics calculus allows for modeling a core at a better overall averaged fuel temperature (which is, actually, not at 900 K), and with a fuel temperature distribution over core assemblies. Those two improvements are estimated to account respectively for 1/3 and 2/3 of the total thermal-hydraulic biases. Figure 5 shows biases of non multi-physics calculations compared to neutronics thermal-hydraulics, on plutonium isotopic composition at reactor discharge, for UOX and MOX fuel assemblies, for two different fuel management: TIER and PARITE [8]. Biases are computed for 10 different random initial isotopic compositions. As illustrated by this figure, potentially significant discrepancies (> 1%) can be observed on unloaded isotopic inventories. Output biases seem to have a strong systematic component along loading isotopy, at fixed fuel type and fixed fuel loading management. Those isotopic biases can also evolve with cycle length or reactor size. PARITE fuel loading management is subject to aggravated biases, probably due to a greater cycle length of 45GWj/t, instead of 35GWj/t for TIER fuel loading management.

It should be noted that unloaded isotopic compositions calculated with multi-physics models are characterized by a systematic degradation of plutonium quality, defined by the ratio of fissile plutonium isotopes (²³⁹Pu, ²⁴¹Pu) over plutonium isotopes. This quantity is a rough observable of spent fuel recycling potential. On UOX and MOX assemblies, for both TIER and PARITE fuel loading managements, fissile plutonium isotopes are negatively impacted by the thermal-hydraulics feedback, whereas other plutonium isotopes are less under-produced, or even over-produced, with regard to previous isothermal core models. This degraded plutonium quality could impact CLASS recycling scenario simulations (see Sect. 4).

3.3. Fuel loading models for fuel cycle dynamics calculations

Thermal-hydraulics coupling could also affect fabrication model accuracy. Indeed, local thermal-hydraulics is known to produce a lower estimation of core power factors on PWRs, due to the negative thermal feedback.

The power factor over assemblies is used by CLASS for MOX fuel fabrication model, which is processed thanks to spent UOX plutonium [7]. On the one hand, plutonium isotopic composition $\overrightarrow{\mathrm{Pu}}$ to be loaded in MOX fresh fuels, is defined by stock compositions during the CLASS simulation. On the other hand, plutonium content %_Pu inside the MOX assembly is a variable computed by the fabrication models. This %_Pu rate is adjusted in order to minimize the core radial power factor, as a function of the UOX fresh composition, using neural networks trained on data produced by DONJON isothermal calculations [15]. Figure 6 shows the radial power factor averaged over the whole irradiation cycle, as a function of the plutonium content %_Pu in the MOX, for six different loaded isotopic compositions $\overrightarrow{\mathrm{Pu}}$. As expected, the power factor is systematically lowered when using the thermal-hydraulics coupling simulations, compared to the isothermal ones. %_Pu minimizing the power factor is however not clearly discernible between coupled simulations from corresponding isothermal ones. Associated bias on %_Pu falls beyond the neural network accuracy that is presented in [6]. Isothermal MOX fabrication models thus remain satisfactory for scenario simulations presented in Section 4: neural networks predicting optimized %_Pu does not need to be trained again on a new database of reactor cores simulations containing thermal-hydraulics models.

Fig. 6. Radial power factor averaged over irradiation cycle, over MOX plutonium level. 900 MWe PWR with PARITE fuel loading management, for six different loading isotopies, iso1 to iso6, for isothermal and THM irradiation models. MOX level stands for the plutonium content inside the MOX fresh fuel.

CLASS UOX fabrication models are used to obtain the critical ²³⁵U enrichment for a given target cycle length, for a given loading pattern without pre-determined enrichment [6]. This search for critical enrichment is, again, made thanks to DONJON isothermal core models. This process is costly, because a complete 8-cycle core simulation is needed to bring the core to an accurate loading – and thus unloading – steady state, for each enrichment value which is tested by a dichotomous algorithm. Compared critical enrichment between isothermal core models and thermal-hydraulics core models always lead to relative differences below 0.05%, for cycle lengths between 250FPED and 400FPED⁷. As thermal-hydraulics coupling increases computation time by a factor of 2, UOX fabrication models are, again, kept isothermal, when modeling nuclear reactor fleets in Section 4.

4. Scenario simulation

The new DONJON-CLASS coupling method developed for this work allows efficient neutronics and thermal-hydraulics reactor-coupled simulations that can be used for fleet-scale simulations. Calculations on elementary scenarios, then on a more complete transition scenario, are produced and analyzed here in order to assess the impact of this new model with regards to the previous isothermal hypothesis for reactor models. Comparison between the two types of models allows a quantification of the biases that may be compared to other source of uncertainties such as in [8] or in [16].

4.1. Elementary scenario description

Elementary scenarios model a reduced number of reactors, over a short period of time. This simplified nuclear fleet facilitates the data interpretation, and makes it easier to isolate the impacts attributed to the new reactor models, by reducing variations due to interactions between the plants and factories simulated.

The elementary fleet is composed of 10 900 MWe PWRs loaded with UOX assemblies, fed by an UOX fabrication plant linked to an infinite uranium stock. The fleet additionally contains one heterogeneous 900 MWe PWR, loaded with 30% MOX – 70% UOX fuel assemblies, and following a PARITE fuel loading management [8]. Its fuel comes from the UOX reprocessing plant that ensures the plutonium separation and from a MOX fabrication plant. Spent UOX and MOX fuels are transferred to their cooling pools after irradiation, in which they stay for a few years. Spent fuels being cooled in those pools cannot be used to fabricate recycled fuel. They are eventually transferred to definitive stocks where they can be retrieved for recycling. A schematic view of such a scenario is presented in Figure 7. The MOX fueled reactor starts in year 2000 whereas those fueled with UOX start in year 1978. The 22-year delay allows the build-up of the UOX spent fuel stockpile, to ensure the amount of plutonium for MOX fuel fabrication.

Fig. 7. Schematic view of elementary fuel cycle simulations. The heterogeneous reactor start-up is delayed 22 years after the start-up of UOX reactors.

Fabrication models, for UOX and MOX fuels, are based on isothermal simulations, as justified in Section 3. Irradiation models, for UOX and MOX reactors, can use new DONJON simulations containing thermal-hydraulics improvements. Scenarios with those new irradiation models are compared to equivalent scenarios with isothermal irradiation models.

Several fuel cycle parameters are used as variables in the CLASS simulations to create a batch of different trajectories: the fuel fabrication time (2–4 years), the fuel cooling time (5–10 years), and the fuel inventory management (LIFO or FIFO⁸) when collecting plutonium from UOX spent fuel for MOX fuel fabrications. Variables are taken identical between the UOX fuel assemblies of UOX reactors, the UOX fuel assemblies of the UOX/MOX reactor, and the MOX fuel assemblies of the UOX/MOX reactor. Overall, 8 set of parameters have been considered, defining 8 different fuel cycle scenarios that may be simulated either with isothermal reactor models or with neutronics thermal-hydraulics coupled models. Those simple scenario variations are implemented to avoid falling onto a specific case, and to assess overall thermal-hydraulic biases for fleet simulations.

The scenario ends after the fabrication and cooling times to ensure the same number of MOX irradiation cycles over all the scenario parametrizations, which corresponds approximately to 2030. All reactors evolve at full power from fuel reloading to target discharge burn-up: load following is not studied here.

4.2. Simulation results

The analysis of these scenarios focuses on isotopic inventories in nuclear plants and fuel cycle factories. The total plutonium inventory biases contained in spent UOX and MOX fuel are displayed in Figure 8. The relative plutonium content of the thermal-hydraulics-coupled scenarios is presented, with regard to the isothermal scenarios.

Fig. 8. Relative difference of plutonium inventory in spent fuel stockpile between fuel cycle simulations with isothermal and multi-physics reactor models. The reference simulation is the multi-physic one.

Unsurprisingly, patterns on unloaded plutonium inventories match the results obtained on a reactor scale (see Sect. 3). Plutonium is slightly under-produced in UOX assemblies, for scenarios with thermohydraulics couplings, compared with scenarios simulated with the isothermal models. In MOX spent fuel, the total plutonium inventory does not seem strongly affected, although, as seen in Section 3, odd isotopes are under-produced, and even isotopes are over-produced, also resulting in a loss of plutonium quality.

Thermohydraulics scenario simulations also reveal new interactions between fabrication and irradiation models. Retrieving a fixed quantity of plutonium in both thermohydraulics and isothermal scenarios mathematically increases the relative error between the two models, despite the fact that the absolute error remains unchanged because the reference value of the relative error is lower. But here, additionally, in order to compensate for the loss of plutonium quality in spent UOX and to keep the energetic equivalence with the UOX fuel of MOXed reactors, fabrication models adjust the %_Pu value and increase plutonium content in new MOX assemblies. This aggravates the lack of plutonium in spent UOX stocks. As explained in Section 3, fabrication models are not impacted by thermohydraulics biases when considering a fixed input plutonium isotopy; but here, the deterioration of input plutonium isotopy leads to a significant impact on the output isotopic content of the fabricated MOX fuel.

4.3. Full scale transition scenario definition

For this section, a more complex scenario has been considered. It is an academic translation of an institutional French scenario produced by RTE [17].

The French fleet is composed of numerous PWRs (56 on this day), of various thermal powers and various effective production rates. In this work, this French nuclear fleet is fully modeled. Its schematic view is presented in Figure 9. Some assumptions are made for this work, here are the main one’s. First, the evolution of historically installed nuclear power plant capacity is modeled since the beginning of PWR exploitation. Assumptions are also made to define the future installed nuclear capacity until 2050, based on the work of scenario N03 of RTE [17]⁹, which notably includes the hypothesis that reactors will be shut down at age 60. Second, changes in fuel management are also modeled: the progressive introduction of MOX assemblies on some 900 MWe reactors of the fleet (following a PARITE fuel loading management [8] with 70%UOX/30%MOX assemblies) is taken into account. Hypothesis such as MOXing some 1300 MWe PWRs, in order to compensate for the closure of 900 MWe MOXed PWRs, have also been made. The ERU fuel assemblies, loaded into some of the French reactors, are not modeled, and here replaced by UOX assemblies. The overall scheme linking UOX infinite source, UOX spent stocks, MOX fuels, and other cycle plants remains the same as in the elementary simulations. The fabrication and irradiation models used are also similar to those of the elementary scenarios. Even if the French PWRs fleet contains reactors of 900 MWe, 1300 MWe, 1450 MWe, and soon, 1650 MWe, all irradiation models are based on 900 MWe PWRs and isotopic inventories are up-scaled to the real reactor thermal power. This hypothesis should bring limited uncertainties and allow us to clearly show the impact of multi-physics reactor models on plutonium inventories without adding new parameters such as different loading patterns for each of the different types of PWRs.

Fig. 9. Schematic view of the complete transition scenario.

The scenario is simulated twice: a first time with an isothermal irradiation model and a second one with multi-physics irradiation models, each time for all the reactors. The comparison between plutonium inventories in spent fuel stockpiles allows to quantify properly the biases induced by this isothermal hypothesis at the full fleet level.

4.4. Full-scale scenario simulation results

The results are presented in Figures 10–12 It confirms results obtained on elementary scenarios. Firstly, the plutonium quality of spent UOX fuel is systematically smaller leading to a smaller plutonium inventory with multi-physics reactor models as shown in Figure 10.

Fig. 10. Plutonium quality of UOX spent fuel stockpile between isothermal and multi-physics reactor models in a complete fuel cycle dynamic simulation. The reference simulation is the multi-physic one.

Fig. 11. Relative difference in fabrication plant plutonium inventory between isothermal and multi-physics reactor models. The reference simulation is the multi-physic one.

Fig. 12. Differences in spent fuel stockpile plutonium inventories between isothermal and multi-physics reactor models in a complete fuel cycle dynamic simulation. The reference simulation is the multi-physic’s one.

Figure 11, which shows the difference in plutonium inventory in the fabrication plant, clearly underlines the compensation of a lower plutonium quality by a higher plutonium content in fresh MOX fuel assemblies, as it has been evoked in Section 4.2.

Then, results show that the amplitude of the spent UOX stock thermal-hydraulics biases is related to the ratio of MOX assemblies in the fleet, to the total number of assemblies. In those fleet-scale scenarios, for which the ratio of MOX assemblies stabilizes around in 10% (instead of 2.7% for the elementary scenarios), differences of plutonium quantity in spent UOX stocks induced by thermal-hydraulics modeling reach more than 10%. This figure is far from negligible in front of other assumptions such as the loading pattern in the reactor for instance [8, 18]. Biases of the total plutonium inventory in spent fuel stockpile is smaller and remain within a few % at the end of the scenario.

As it has been stressed, nuclear fleet scenarios with thermohydraulics coupling allow a significant improvement of previous isothermal models. The total plutonium inventory is slightly diminished, but the distribution of plutonium between the stocks is changed, and more plutonium is retrieved to make the MOX fuel when using multi-physics models.

The impact of considering a fuel assembly local temperature, instead of an overall core temperature, reiterates the need for scenario simulations that can take into account the variations of the thermal power of each nuclear reactor in the fleet. Indeed a change of power results in a change of overall and local temperatures, leading to variation in the reactor irradiation models used for the scenario simulation.

5. Conclusion

This paper quantifies precisely reactor model biases in fuel cycle dynamic simulations induced by the isothermal hypothesis in core depletion simulations. Fuel cycle simulations, or scenario simulations for the future evolution of any national fleet, are based on complex modeling of fresh fuel fabrication and fuel irradiation. On one hand, these models identify the composition of the fuels to be loaded as a function of the desired cycle time and the isotopic composition in the stocks of materials to be recycled, and, on the other hand, calculate the discharged composition as a function of the loaded one. Previous couplings between CLASS and DONJON have shown the importance of 3D full core modeling in uncertainty reduction induced by reactor models; nevertheless, the reactor simulations were done without 3D thermal-hydraulics considerations. The development of a multi-physic depletion scheme raised the issues of model biases when omitting the local fuel temperature distributions in the core.

Our results show that reactivity evolution is not affected by this local temperature distribution. The same initial fuel composition reaches the same global burn-up whereas the core simulation takes into account thermal-hydraulics coupling or not. However, fuel irradiation, and specifically plutonium isotopic composition, is affected by this hypothesis. As a result, omitting local fuel temperature distribution leads to non-negligible biases of a few percent on plutonium composition at reactor discharge. More importantly, the plutonium quality of spent UOX and MOX fuel assemblies, defined as the fissile proportion in the plutonium vector, is systematically degraded with the new advanced multi-physics calculation scheme. Those observations are not dependent on the plutonium initial isotopic vector.

To propagate those reactor model biases in fuel cycle simulation, the 3D DONJON computation costs had to be reduced by improving the interpolation process of the diffusion data. To do so, DONJON was given perturbed MCPOs instead of full MCPOs, optimized with two smart upstream overlays, allowing for recombination of local MCPOs and for pre-interpolation of redundant dimensions inside a single reactor cycle, in order to get a representative light MCPO, containing very few evolving parameters such as burn-up, boron concentration and effective fuel temperature. This new interpolation methodology does not create significant bias, and allows for a much more efficient coupling between CLASS and DONJON.

It is then possible to propagate reactor biases on a scenario simulation by simulating each fuel cycle trajectory with and without the advanced reactor simulation scheme developed for this work. Two type of scenario were considered here: an academic scenario without dynamic transition with only one heterogeneous reactor loaded with UOX and MOX spent fuel and a much more complex scenario inspired by an institutional French scenario. Results show that the systematic bias on the plutonium quality of spent UOX has to be compensated by a higher plutonium content in fresh MOX. In the full-scale scenario, this bias leads to a discrepancy of several percent in the UOX spent fuel inventory in 2050. This work shows that the isothermal hypothesis in reactor modeling for scenario simulations induces biases of the same order of magnitude as other reactor parameters like the reactor loading pattern. Further work should be dedicated to propagating those biases in plutonium multi-recycling scenarios.

Funding

The author would like to thank the NEEDS French program and the Natural Sciences and Engineering Research of Council of (NSERC) of Canada for partly funding this work.

Conflicts of interest

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

Data availability statement

This article has no generated nor analyzed data associated with this article.

Author contribution statement

All the authors have read and approved the final manuscript, they all have participated in the simulations analysis that have been performed by G.B.

References

All Tables

All Figures

	Fig. 1. 2-level calculation scheme principle for 3D core simulation.
In the text

Fig. 2. Biases between DONJON simulations using different sampling for MCPOs, on plutonium isotope inventories at reactor discharges. Biases are displayed, for each isotope, over the absolute value of the relative distance between loaded and reference isotopic value, for this same isotope. Biases are displayed for perturbative MCPO (in orange) and for recombined perturbative MCPO (in green), relative to a full MCPO with regular grid sampling.

In the text

	Fig. 3. DONJON-CLASS chaining method, with recombination code overlay, used to produce efficient irradiation models.
In the text

	Fig. 4. DONJON thermal-hydraulics neutronics computation scheme. The reactor is boron-critical, meaning that the boron concentration in the coolant is adjusted to maintain neutronics criticality and a controlled chain reaction (keff = 1). EOC stands for End Of Cycle (or irradiation).
In the text

Fig. 5. Histograms of unloaded isotopic relative differences for a variety of loading isotopic compositions, when adding thermal-hydraulics coupling, with regard to isothermal core models, plutonium isotopes, for UOX and MOX fuels, for PARITE and TIER fuel loading managements. The fuel loading management correspond to 4-rounds (PARITE) and 3-rounds (TIER) MOX and UOX irradiations, as described in [14]. Histograms are stacked with a width of 0.1% relative difference value. The references are the multi-physic simulations.

In the text

	Fig. 6. Radial power factor averaged over irradiation cycle, over MOX plutonium level. 900 MWe PWR with PARITE fuel loading management, for six different loading isotopies, iso1 to iso6, for isothermal and THM irradiation models. MOX level stands for the plutonium content inside the MOX fresh fuel.
In the text

	Fig. 7. Schematic view of elementary fuel cycle simulations. The heterogeneous reactor start-up is delayed 22 years after the start-up of UOX reactors.
In the text

	Fig. 8. Relative difference of plutonium inventory in spent fuel stockpile between fuel cycle simulations with isothermal and multi-physics reactor models. The reference simulation is the multi-physic one.
In the text

	Fig. 9. Schematic view of the complete transition scenario.
In the text

	Fig. 10. Plutonium quality of UOX spent fuel stockpile between isothermal and multi-physics reactor models in a complete fuel cycle dynamic simulation. The reference simulation is the multi-physic one.
In the text

	Fig. 11. Relative difference in fabrication plant plutonium inventory between isothermal and multi-physics reactor models. The reference simulation is the multi-physic one.
In the text

	Fig. 12. Differences in spent fuel stockpile plutonium inventories between isothermal and multi-physics reactor models in a complete fuel cycle dynamic simulation. The reference simulation is the multi-physic’s one.
In the text