Development of a Reduced Order Model for Fuel Burnup Analysis

Abstract

Fuel burnup analysis requires a high computational cost for full core calculations, due to the amount of the information processed for the total reaction rates in many burnup regions. Indeed, they reach the order of millions or more by a subdivision into radial and axial regions in a pin-by-pin description. In addition, if multi-physics approaches are adopted to consider the effects of temperature and density fields on fuel consumption, the computational load grows further. In this way, the need to find a compromise between computational cost and solution accuracy is a crucial issue in burnup analysis. To overcome this problem, the present work aims to develop a methodological approach to implement a Reduced Order Model (ROM), based on Proper Orthogonal Decomposition (POD), in fuel burnup analysis. We verify the approach on 4 years of burnup of the TMI-1 unit cell benchmark, by reconstructing fuel materials and burnup matrices over time with different levels of approximation. The results show that the modeling approach is able to reproduce reactivity and nuclide densities over time, where the accuracy increases with the number of basis functions employed.

1. Introduction

In nuclear reactor analysis, the time evolution of the fuel composition and the reactivity during the in-core reactor irradiation is relevant for the fuel management and for determining the amount of long-lived radionuclides in spent nuclear fuel [1,2]. Furthermore, the neutronics characterization of a nuclear system depends on the isotopic composition of the fuel material. It follows that, in nuclear reactors, burnup calculations are needed to evaluate the concentrations of the various nuclides over time.

Over the years, dedicated tools have been developed to handle burnup analysis [3,4,5,6,7], by obtaining reaction rates from neutronics and solving burnup equations through different numerical methods and algorithms. On the one hand, neutronics can be simulated by deterministic codes that solve approximations of the transport equation and implement a multi-stage computational procedure from meso-scales (pin cell, fuel assembly) to macro-scales (reactor core) [8]. These codes are widely used for industrial applications, since they provide results in a fast and reasonable computational times. On the other hand, in the last years, the increase of the computational power has led a growing interest in the adoption of Monte Carlo codes, where a probabilistic approach allows one to solve the transport problem by simulating the life histories of individual neutrons of the system. Indeed, Monte Carlo codes have the advantages to be flexible for different reactor geometries, where all scales for neutronics can be resolved once. In addition, they produce accurate results; for this reason, they are also used as reference for benchmarking results of deterministic codes.

Furthermore, in burnup analysis, Monte Carlo codes provide a detailed pin-by-pin characterization of the fuel consumption in full core calculations. Recently, their capabilities have been extended to perform local burnup calculations, subdividing the fuel pins into radial and axial regions, in order to reach a high level of accuracy in the study of the fuel cycle. It follows that the problem of the computational cost is enhanced because the information processed is huge for the total number of burnup regions, in the order of millions or more [9,10]. Moreover, if multi-physics approaches are adopted to consider the local effects of temperature and density fields on fuel consumption, the computational load grows further [11,12,13]. So the need to find a compromise between computational cost and solution accuracy remains a critical issue in reactor analysis, as it does for deterministic codes, and it is not only limited to burnup, but for each problem where it is important to obtain a reliable solution in a fast running time.

To this aim, Reduced Order Methods (ROMs) [14,15,16] have been recently employed in different fields of engineering. The philosophy behind ROMs is that the behavior of a system can be well described by a small number of dominant modes, which allow approximating a physical problem, described by a PDE (Partial Differential Equation) or a set of ODEs (Ordinary Differential Equations), into a system with fewer degrees of freedom. In this way, calculations become faster, so in some cases the capability of a personal computer is enough for simulations, instead of employing HPC technologies.

ROMs usually adopt an offline/online strategy. In the offline stage, the most computationally expensive, the basis functions are generated from solutions (i.e., snapshots) of the Full Order Model (FOM). In the online stage, the simulation of the ROM (we refer to "ROM" as Reduced Order Method or Reduced Order Model, depending on the context) can be run as many times as required.

In the ROM framework, the Proper Orthogonal Decomposition (POD) with the snapshots technique [17,18,19] is one of the most valuable methods. Indeed, POD has the powerful capability to select the most energetic modes that represent the main significant features of the system [20]. After that, a low-dimensional approximation of the original system can be obtained by its Galerkin projection into the space spanned by the POD basis functions.

In the nuclear field, POD has been successfully applied to neutronics [21,22,23], thermal-hydraulics [24], sensitivity and uncertainty quantification [25,26]. Notwithstanding, there is a lack of literature in the employment in burnup calculations, nevertheless it would be necessary to solve the issue of computational burden previously discussed.

Furthermore, in this work, we are motivated to adopt POD by the following consideration. In burnup calculations, a first linear differential problem is usually expressed in matrix notation for each time step, where the transmutation terms of the burnup matrix correspond to the multiplication between the cross sections and the neutron flux. If we consider the neutron flux expressed as a linear combination of basis functions, from the linearity of the transport equation, we can obtain a linear combination of burnup matrix basis functions with the same linear coefficients of the flux basis functions.

Thanks to this consideration, the proposed methodology aims to reconstruct the burnup matrix from the linear combination of few POD matrix basis functions, by calculating the linear coefficients from the projection of fluxes at each time step. This may decrease the computational burden with respect to the standard procedure, where all entries of the burnup matrix (dimensions ~ 1500 × 1500) are computed by the neutron transport.

Furthermore, it would be interesting to find POD basis functions for nuclide vectors as well, in order to decrease the total number of the materials studied in burnup analysis. For this reason, we also investigate the adoption of POD to obtain few representative basis functions to reconstruct nuclide vectors.

Finally, it is important to provide a physical interpretation of the different POD modes of interest, since they carry out physical information of the system.

We test the approach on 4 years of burnup of the 2D TMI-1 unit cell benchmark, by simulating neutron transport with the Serpent Monte Carlo code [3]. In order to provide the axial power profile of the 3D fuel pin in 2D modeling, we reproduce different power levels. Furthermore, we obtain the radial power by dividing the fuel into radial regions. The reference case is modeled with uniform temperature and density fields, since the present work, focusing only in burnup analysis, does not consider a multi-physics characterization, as reported in

Table 1. Neutronics and thermal-hydraulics modeling in the reference case.

Physics Area	Model Assumption
Neutronics	Monte Carlo
Thermal-Hydraulics	Uniform temperature/density fields

.

In the FOM, we carry out the burnup analysis with the Serpent Monte Carlo code. In the ROM, we use Serpent to solve neutron transport and MATLAB for burnup calculations. We compare ROM and FOM results for reactivity and nuclide concentrations over time. Indeed, FOM against ROM is the standard procedure employed in the evaluation of the quality of the approximation in ROM [21,24]. Nevertheless, it should be noted that comparisons with experimental data and other approaches for burnup analysis, such as those developed for neutronics deterministic solvers, are mandatory for the assessment of the methodology. However, since the present work is focused on the development of a new approach and its application on a simplified test case for demonstration purposes, extensive verifications and validations are demanded of future work.

The paper is organized as follow. In Section 2, we describe the test case and we present in detail the methodological approach adopted in the work. In Section 3, we show and analyze the results of the ROM in comparison with FOM. Section 4 discusses the computational costs of simulations. Section 5 is left to conclusions and final remarks.

2. Methodological Approach

2.1. Motivation of the Methodology

In fuel burnup analysis, the following set of first order linear differential equations, called Bateman equations [27], are solved to calculate the time evolution of nuclide concentrations:

$${\displaystyle \frac{d{N}_j}{dt}}={\displaystyle \sum }_{i\ne j}^k(f_{ij}{\lambda }_i+{\sigma }_{ij}\varphi )N_i-({\lambda }_j+\varphi {\displaystyle \sum }_i^k{\sigma }_{ji})N_j Nj\left(0\right)=N_0 j=1,\dots ,k$$

(1)

where $N_j$ is the concentration of j-th nuclide, ${\lambda }_i$ the decay constant of the i-th nuclides, $\varphi$ the neutron flux, ${\sigma }_{ij}$ the microscopic cross section for the reaction that generates the j-th nuclides from the i-th nuclide, $f_{ij}$ the decay fraction from the i-th nuclide to the j-th nuclide and k the total number of nuclide inventories. In compact form, the equation can be expressed as follows:

$${\displaystyle \frac{d{N}_j}{dt}}={\displaystyle \sum }_{i\ne j}{\gamma }_{ij}{N}_i-{\gamma }_j{N}_{j}Nj\left(0\right)=N_0 j=1,\dots ,k$$

(2)

where ${\gamma }_{ij}$ and ${\gamma }_j$ are the coefficients of production-destruction rates. The calculations performed by neutron transport codes, Monte Carlo or deterministic, divide a burnup history into time steps. For each of them, the coefficients ${\gamma }_{ij}$ and ${\gamma }_j$ are fixed, constant and are computed by the neutron transport. It is important to note that ${\gamma }_{ij}$ and ${\gamma }_j$ are constant for a single solution of the Bateman equations in a single burnup step. However, more solutions are usually calculated for a single burnup step, by employing predictor-corrector methods [28]. Equation (2) can be written in matrix notation:

$$\underline{\dot{n}}=\underline{\underline{A}}·\underline{n} {\underline{n}}_0=\underline{n}\left(0\right)$$

(3)

where $\underline{n}$ is the nuclide vector of dimension k, with entries corresponding to nuclide densities; $\underline{\underline{A}}$ is the k × k burnup matrix, which contains the coefficients ${\gamma }_{ij}$ and ${\gamma }_j$; and ${\underline{n}}_0$ the vector of initial concentrations. The running time of a burnup analysis is mainly affected by the solution of neutron transport, since Equation (3) is solved very quickly by using different algorithms [29]. In terms of memory usage, the solution of the Bateman equations becomes problematic in full core calculations, if the needed RAM memory reaches hundreds of gigabytes to compute reaction rates for thousands of burnup regions.

The methodological approach of this work aims to reduce the memory consumption of burnup calculations, by approximating $\underline{\underline{A}}$ as a linear combination of basis functions:

$$\underline{\underline{A}}\simeq {\displaystyle \sum }_{i=1}^{r_{\mathsf{\Pi }}}{a}_j·\underline{\underline{{\mathsf{\Pi }}_j}}$$

(4)

where $a_j$ is the time dependent coefficient of the $\underline{\underline{{\mathsf{\Pi }}_j}}$ basis function and $r_{\mathsf{\Pi }}$ the number of basis functions employed. In this way, it is possible to restrict the computation only for the coefficients $a_j$, instead of calculating all entries of burnup matrices. For this scope, we consider the burnup matrix without the decay terms in Equation (1), in order to only take into account its dependency on neutron flux. Indeed, if we also express the neutron flux in energy groups as a linear combination of basis functions, from the linearity of transport equation, we obtain:

$$\left\{\{\begin{array}{cc}\hfill & \underline{\underline{A}}=f\left(\underline{\varphi }\right)\hfill \\ \hfill & \underline{\varphi }\simeq {\displaystyle \sum }_{j=1}^{r_{\psi }}{b}_j·\underline{{\psi }_j}\hfill \end{array}\right.\Rightarrow \underline{\underline{A}}\simeq f({\displaystyle \sum }_{j=1}^{r_{\psi }}{b}_j·\underline{{\psi }_j})\simeq {\displaystyle \sum }_{j=1}^{r_{\psi }}{b}_j·f\left(\underline{{\psi }_j}\right)={\displaystyle \sum }_{i=1}^{r_{\psi }}{b}_j·\underline{\underline{{\mathsf{\Pi }}_j}}$$

(5)

where, through a suitable choice of basis functions, the linear coefficients of Equation (4) and Equation (5) are the same; i.e., $a_j=b_j$. In this way, the coefficients can be computed from the projection of fluxes onto their basis functions, since the dimensions of flux vector, corresponding to the number of neutron energy groups, can be minor respect to the dimensions of the burnup matrix. The fluxes can be obtained from a deterministic or a Monte Carlo code, since the previous considerations depend on the physics equations and are generally valid.

As stated in the introduction, we select the basis functions for $\underline{\underline{A}}$ from POD, a technique of the ROM framework, which has the ability to provide the most energetic modes corresponding to the most significant features of a physical system [15]. In particular, since POD is based on a maximum energy criterion, the accuracy of the approximation is usually guaranteed by the first POD modes, with the neglection of higher order modes.

Moreover, it is possible to associate a physical meaning to POD modes. Indeed, POD was introduced to study the coherent structures in experimental turbulent flows. In the nuclear field (as in other fields of engineering), it has been shown that POD basis functions carry out physical information on distribution of scalar and energy group fluxes, fluid velocity distributions in thermal-hydraulics and cross section uncertainties [21,22,23,24]. For this reason, we expect that POD modes for burnup matrices and fluxes provide the main physical information on fuel burnup. This assumption also suggests to investigate the possibility to obtain representive basis functions for nuclide vectors; i.e., fuel materials. In this way, the compositions of fuel materials can be approximated as follows:

$$\underline{n}\simeq {\displaystyle \sum }_{j=1}^{r_{\xi }}{c}_j·\underline{{\xi }_j}$$

(6)

where $\underline{{\xi }_j}$ and $r_{\xi }$ are the first basis functionsm and $c_j$ is the coefficient of the linear combination.

In conclusion, we have to evaluate the level of approximation and the physical information of POD with different number of basis functions for burnup matrices, fluxes and nuclide vectors. In the next section, we present the scheme of the approach adopted to implement POD both for fuel materials and burnup matrix/fluxes.

2.2. Reference Case Description

In order to assess the methodological approach based on ROM, we carry out the burnup analysis on the two-dimensional TMI-1 PWR cell (Figure 1a), described by the Exercise I-1b in the framework of OECD benchmarks for UAM (Uncertainty Analysis in Modelling) of LWR [30]. The design parameters are reported in

Table 2. Design parameters of TMI-1 PWR unit cell.

Parameter	Value
Fuel type	UO2
Fuel density	10.283 g/cm${}^3$
Enrichment in 235U	4.85%
Pin pitch	1.4427 cm
Fuel Pellet Diameter	9.391 mm
Cladding Inner Diameter	9.582 mm
Cladding Outer diameter	10.928 mm
Cladding material	Zircaloy-4

.

We adopt Serpent in the simulations (Figure 1b), a continuous-energy Monte Carlo code developed by VTT Technical Research Centre of Finland, which we proved in the past to be a reliable tool for reactor analysis [31,32,33]. Indeed, this code provides powerful features to perform neutron transport and burnup simulations, in addition to its fast running time and reactor geometry pre-implementation. Moreover, Serpent allows one to easily set the burnup history, to obtain the burnup matrices and concentrations for all nuclide inventory adopted at each time step, in which the burnup history is divided.

Serpent reads the continuous-energy interaction data from the ACE format cross section libraries; in burnup analysis, radiocative decay and fission yield data are imported from the standard ENDF format data libraries. In this work, we employ JEFF-3.1.1 [34] library for the transport and depletion calculations.

Since the model of the fuel cell is two dimensional, we simulate five power levels at 52, 150, 233, 294 and 326 W/cm, in order to reproduce the axial profile of the 3D pin during the offline phase. In particular, they correspond to five different burnup histories. The powers are obtained from the assumption of the typical axial cosine shape for the average fuel pin in TMI-1 reactor.

The time points of the burnup histories are 53 by covering 4 years of depletion for all nuclide inventory, composed of 1555 nuclides, with initial finer step widths (1, 1, 5, 8 and 15 days) and following time steps of 30 days. Fuel material is divided into five radial regions of equal area. For burnup calculations, we adopt the constant extrapolation scheme, in which the Bateman equation is solved with Beginning Of Step (BOS) cross sections and fluxes. In addition, the xenon is set at the equilibrium. From the transport calculations of each burnup region, we detect fluxes with the ECCO1968 group structure [35], which have the finest energy discretizations over all energy pre-defined group structures available in Serpent.

At each burnup step, we simulate 40 million active neutron histories, divided into 2000 source neutrons per 20,000 cycles. The inactive cycles are 200 in number. This statistics assures uncertainties less than 10 pcm on k_eff.

2.3. Scheme of the Methodological Approach

The reduction method of this work is based on an offline/online strategy. In the offline phase, the most computationally expensive POD modes are generated from solutions (i.e., snapshots) of the FOM. In the online phase, completely decoupled from the previous stage, the ROM is run as many times as required.

In general, ROM is obtained by a Galerkin projection of the original system onto the space spanned by snapshots, where a PDE (Partial Differential Equation) or a set of ODEs (Ordinary Differential Equations), is reduced to a system with fewer degrees of freedom. Differently, in this case, we refer to the ROM as the system obtained by the linear combination of basis functions in Equations (4)–(6). In other words, the ROM consists of the reconstruction of suitable approximations of burnup matrix/fluxes and nuclide vectors at different time step.

The next paragraph describes in detail the offline and the online stages of the modeling approach, shown schematically in Figure 2 and Figure 3.

2.3.1. Offline Phase

The offline phase involves the computation of POD modes of $\underline{\underline{A}}$, $\underline{\varphi }$ and $\underline{n}$ from solutions of the FOM; i.e., snapshots. In particular, we use the latter to build a snapshots matrix, in order to perform a Singular Value Decomposition (SVD), usually used in POD for basis functions generation [19].

In this case, for neutron fluxes and burnup matrices, a snapshot is referred to a burnup step between two time points, for a burnup region at a specific power level. In this way, the total number of snapshots is 1300 (52 burnup steps × 5 burnup regions × 5 power levels).

Differently, for nuclide vectors, we consider the concentrations at each time point, instead of burnup step. It follows that we also include the concentrations at the end of the burnup history. In this way, we obtain a total of 1325 snapshots (53 time points × 5 burnup regions × 5 power levels).

For the sake of clarity, we summarize the offline phase, reported in Figure 2, with the following steps:

• Generation of snapshots for $\underline{\underline{A}}$, $\underline{n}$ and $\underline{\varphi }$ from Monte Carlo burnup calculations of the FOM.
• Construction of the matrix of snapshots $\underline{\underline{S_A}}$, $\underline{\underline{S_{\varphi }}}$ and $\underline{\underline{S_n}}$. Depending on the different types of basis functions that we investigate, we get the following snapshots matrices:
  • $\underline{\underline{S_n}}=[{\underline{n}}^1,{\underline{n}}^2,\dots ,{\underline{n}}^{1325}]$, with dimensions 1555 × 1325, for snapshots of nuclide vectors;
  • $\underline{\underline{S_A}}=[{\underline{A}}^1,{\underline{A}}^2,\dots ,{\underline{A}}^{1300}]$, with dimensions 2,418,025 × 1300, obtained by reshaping each burnup matrix $\underline{\underline{A^i}}$ into a single column A${}^i$;
  • $\underline{\underline{S_{\varphi }}}=[{\underline{\varphi }}^1,{\underline{\varphi }}^2,\dots ,{\underline{\varphi }}^{1300}]$, with dimensions 1968 × 1300, for snapshots of fluxes.
• Application of SVD to obtain $\underline{\underline{S_A}}$, $\underline{\underline{S_{\varphi }}}$ and $\underline{\underline{S_n}}$ to obtain:

  $$\underline{\underline{S_A}}=\underline{\underline{U_A}}·\underline{\underline{{\mathsf{\Sigma }}_A}}·{\underline{\underline{V_A}}}^T$$

  (7)

  $$\underline{\underline{S_{\varphi }}}=\underline{\underline{U_{\varphi }}}·\underline{\underline{{\mathsf{\Sigma }}_{\varphi }}}·{\underline{\underline{V_{\varphi }}}}^T$$

  (8)

  $$\underline{\underline{S_n}}=\underline{\underline{U_n}}·\underline{\underline{{\mathsf{\Sigma }}_n}}·{\underline{\underline{V_n}}}^T$$

  (9)

  where the i-th columns of ${\underline{\underline{U}}}_A$, ${\underline{\underline{U}}}_{\varphi }$ and ${\underline{\underline{U}}}_n$ are the POD modes $\underline{u_A^i}$, $\underline{u_{\varphi }^i}$ and $\underline{u_n^i}$, while $\underline{\underline{{\mathsf{\Sigma }}_A}}$, $\underline{\underline{{\mathsf{\Sigma }}_{\varphi }}}$ and $\underline{\underline{{\mathsf{\Sigma }}_n}}$ are the diagonal matrices of the singular values $s_A^i$, $s_{\varphi }^i$ and $s_n^i$. The latter are sorted in descending order of relevance and represent the energy of each mode; i.e., the information of the system carried by the mode itself. The square of the singular value $s_i$ corresponds to the eigenvalue ${\beta }_i$. $\underline{\underline{V}}$ is the matrix containing the coefficients of the linear combination of basis functions.
• Referring to Equation (5), from the projection of $\underline{\varphi }$ onto suitable $\underline{\varphi }$ basis functions, it is possible to obtain the same linear coefficients of $\underline{\underline{A}}$. Firstly, it is necessary to verify if the POD modes of $\underline{\varphi }$ are suitable for the scope. Otherwise, the matrix $\underline{\underline{V_A}}$, which contains the coefficients of $\underline{\underline{A}}$, can be employed to create new flux functions by the following manipulation:

  $$\underline{\underline{B_{\varphi }}}=\underline{\underline{S_{\varphi }}}·\underline{\underline{V_A}}$$

  (10)

  where $\underline{\underline{B_{\mathsf{\Phi }}}}=[{\underline{\psi }}_{dep}^1,{\underline{\psi }}_{dep}^2,\dots ,{\underline{\psi }}_{dep}^{N_s}]$. We call ${\underline{\psi }}_{dep}^j$ the j-th depletion flux functions, not necessarily orthogonal like POD modes, but suitable to be employed in burnup calculation of the online phase.

First of all, it is important in this phase to analyze the composition of $\underline{n}$ basis functions, in order to identify the nuclides with the highest concentrations. We expect that these nuclides are the most important in burnup, expecially if columns of $\underline{\underline{S_n}}$ are weighted over the cross sections of fissile or fissionable nuclides (weighted POD). Moreover, also $\underline{\varphi }$ and $\underline{\underline{A}}$ modes have to be associated to their physical meaning. Secondly, before the manipulation in point 4, it has to be evaluated if the POD basis functions of $\underline{\varphi }$ are compatible with $\underline{\underline{A}}$.

2.3.2. Online Phase

In the online phase, we decouple the reconstruction of $\underline{n}$ from the reconstruction of $\underline{\underline{A}}$ by splitting the procedure in two parts.

In the first part, as shown by the flow chart on the left side of Figure 3, we run criticality calculations with approximations of $\underline{n}$ at different time points. To sum up, the procedure consists of the following steps:

• Perform the reconstruction of $\underline{n}$ for each snapshot by approximating $\underline{\underline{S_n}}$ for the first $r_n$ basis functions. This is done by setting to zero the singular values ${\sigma }_n^i$ with $i>r_n$, in order to obtain the reduced matrix ${\underline{\underline{\mathsf{\Sigma }}}}_n^{r_n}$. So we obtain an optimal rank approximation of ${\underline{\underline{S}}}_n$ by its SVD form ${\underline{\underline{S}}}_n^{r_n}$:

  $${\underline{\underline{S}}}_n^{r_n}=\underline{\underline{U_n}}·{\underline{\underline{\mathsf{\Sigma }}}}_n^{r_n}·{\underline{\underline{V_n}}}^T$$

  (11)
• Run Monte Criticality calculations at different time points, corresponding to nuclide concentrations imported by different snapshots of ${\underline{\underline{S}}}_n^{r_n}$.

From the first part of the online procedure, it is possible to evaluate the possibility to restrict the fuel burnup analysis only for some fuel material basis functions. It follows that, if $r_n$ is lower than the number of all fuel materials of the FOM, we save computational load. This is not the case of this article, since here we aim to develop and verify the methodological approach.

In the second part of the online phase, as reported on the right side of Figure 3, we run burnup analysis by solving the burnup equations with approximations of $\underline{\underline{A}}$. The procedure is summarized by the following loop, related to each fuel burnup region at every time step:

• Run Monte Carlo criticality calculations to obtain $\underline{\varphi }$;
• Perform the projection of $\underline{\varphi }$ onto $B_{\varphi }$, to obtain coefficients {a${}_i{\}}_{i=1,\dots ,r}$;
• Construct $\underline{\underline{A}}$ from a linear combination of the first r basis functions, as in Equation (4);
• Solve the Bateman equation to obtain nuclide vectors for the next burnup step;
• Import new densities into the Monte Carlo model and restart at point 1, for the next time step.

Differently from the FOM, we use the neutron transport to calculate only $\underline{\varphi }$ and not each element of $\underline{\underline{A}}$. In this way, the memory requirements depends on the dimension of vector $\underline{\varphi }$, corresponding to the number of energy groups employed.

We point out that the two parts of the online phase are not to be considered completely separated. In future developments, the burnup analysis of the online phase can be restricted to few basis functions of burnup materials, with the aim to decrease the cost of calculations with respect to those needed for the total number of burnup regions in the FOM.

3. Results

At the beginning, we present the results obtained by the offline and the online phase with $\underline{n}$ reconstruction. After that, we show the results from $\underline{\underline{A}}$ reconstruction.

3.1. ROM of Nuclide Vectors

POD Modes (Offline Phase)

From the SVD applied to $\underline{\underline{S_n}}$ (Equation (9)), we obtain the POD basis functions $\{{\underline{\xi }}_1,{\underline{\xi }}_2,\dots ,{\underline{\xi }}_{N_{\xi }}\}$ of nuclide vectors. Figure 4 shows the basis function #1, in function of isotopes, indicated by ZAI index (ZAI = 10,000 × Z + 10 × A + I, where Z is the atomic number, A the mass number and I the isomeric state). The distribution is typical of a nuclide vector during the burnup [29], where the concentrations of fission products have two peaks around the atomic numbers of 40 and 60 respectively. Additionally, other basis functions have the same distributions on a logarithmic scale, without the possibility to distinguish the most important nuclides. For this reason, we plot in Figure 5 the highest ten concentrations for the first five basis functions.

Nuclide basis function $\#1$ shows ${}^{16}$O as the first nuclide, which comprises UO${}_2$, fresh fuel. Indeed, oxygen is not a neutron absorbent, and its concentration remains high along the entire burnup history; i.e., ${}^{16}$O is also present as the third isotope in the basis functions $\#2$ and $\#3$. Furthermore, fissile and fissionable nuclides provide important contributions for the basis functions $\#1$, $\#2$, $\#3$ and $\#4$.

We also find strong presence of ${}^{144}$Nd, ${}^{136}$Xe and ${}^{135}$Cs in the basis functions #4 and #5. It is interested to note that, on one hand, ${}^{144}$Nd is proportional to fuel consumption and is traditionally used as indicator of radial burnup [36]. On the other hand, ${}^{136}$Xe and ${}^{135}$Cs are sensitive indicators of the thermal neutron flux [37].

For a more rigorous analysis, we report in Figure 6 the normalized eigenvalues ${\beta }_i^{norm}$ obtained by POD. They represent the information, called energy, carried out by each POD mode. We observe the typical exponential decrease, with values under 10${}^{-15}$ after the basis function $\#50$.

We perform the snapshots’ reconstruction with 5, 10 and 25 basis functions, to obtain approximations that do not exceed few tens of basis functions. The truncation errors is defined as:

$$e_{POD}=1-{\displaystyle \sum }_{i=1}^{r_n}{\beta }_i^{norm}$$

(12)

where $r_n$ is the number of basis functions employed, and ${\beta }_i^{norm}$ the normalized eigenvalues. The value of e${}_{POD}$ is lower than 10⁻⁷, as reported in

Table 3. Truncation error e_POD respect to number of basis functions N${}_r$.

N${}_{\mathit{r}}$	ePOD
5	1.26 × 10${}^{-8}$
10	1.30 × 10${}^{-12}$
25	1.44 × 10${}^{-13}$

, and it decreases with the number of basis functions. However, we point out the e${}_{POD}$ is not representative of the error of the Reduced Order Model but has only indicative purposes [24].

3.2. Nuclide Vectors Reconstruction (Online Phase)

In order to obtain approximations of nuclide vectors with 5, 10 and 25 basis functions, we adopt the optimal rank approximations of $\underline{\underline{S}}$, with matrices $\underline{\underline{S_n^5}}$, $\underline{\underline{S_n^{10}}}$ and $\underline{\underline{S_n^{25}}}$. From the latter, we run transport calculations with fuel materials at 1, 2, 3 and 4 years, at each power level, corresponding to a total of 60 criticality simulations. The Monte Carlo statistics are the same of the FOM. Figure 7a shows the results for the absolute variation $\mathsf{\Delta }$k_eff in pcm with respect to the FOM, in function of burnup in MWd/kgU.

For all cases, taking into account other verifications against Monte Carlo codes, as in [38,39,40,41,42], variations of hundreds of pcm represent a good level of approximation. Major differences are shown by the reconstruction with five basis functions, reaching absolute variations of criticality greater than 400 pcm for low (<10 MWd/kgU) and high (>60 MWd/kgU) burnup. Moreover, as expected, the best approximation was obtained with 25 basis functions, where $\mathsf{\Delta }$k_eff is between 0 and 200 pcm. Furthermore, the trend of $\mathsf{\Delta }$k_eff changes over time at different approximations. It oscillates with 5 and 10 basis functions; finally, it becomes more stable with 25 basis functions.

We report in

Table 4. Means and standard deviations of $\mathsf{\Delta }$$k_{eff}$ with different numbers of POD basis functions in the vector nuclide reconstruction. The results are also calculated by excluding the highest and lowest burnup. The last column reports the average value of the Monte Carlo statistical error.

	ALL BU		NO Highest/Lowest BU		Serpent MC
N${}_{\mathbf{r}}$	$\overline{\mathsf{\Delta }{\mathit{k}}_{\mathit{eff}}}$[pcm]	${\mathbf{\sigma }}_{\mathsf{\Delta }{\mathit{k}}_{\mathit{eff}}}$[pcm]	$\overline{\mathsf{\Delta }{\mathit{k}}_{\mathit{eff}}}$[pcm]	${\mathbf{\sigma }}_{\mathsf{\Delta }{\mathit{k}}_{\mathit{eff}}}$[pcm]	${\mathit{\sigma }}_{\mathit{MC}}$[pcm]
Bas. Fun. 5	222	151	195	136	16
Bas. Fun. 10	149	112	130	99	16
Bas. Fun. 25	108	30	107	29	16
Weig. Bas. Fun. 5	160	154	114	77	16

the mean value and the standard deviation of $\mathsf{\Delta }$k_eff, ${\sigma }_{\mathsf{\Delta }{k}_{eff}}$, for each set of simulations. The results are in agreement with the previous considerations, where the approximation improves as the number of basis functions increases.

The accuracy of k_eff depends on the approximation of nuclide concentrations, particularly for the main nuclides that affect the criticality. For this reason, we show in Figure 7b the average relative variation (%) with respect to the full order cases, over densities of ${}^{235}$U, ${}^{238}$U and ${}^{239}$Pu. We select the latter nuclides since they are the most important fissile and fissionable nuclides, in relation to their high absorption cross sections and concentrations over time. In particular, Figure 7b shows that the average relative variation improves with the increase of the number of basis functions, as observed for k_eff. In addition, the values are under 1%, assuring a good level of accuracy, with changes of orders of magnitudes over the different levels of approximation.

3.3. Weighted POD

The previous application of POD does not take into account the physics of the reaction rates, which depend on the cross sections of the nuclides. For this reason, we weigh each nuclide, by choosing as weight the fraction between the absorption cross section of the nuclide of interest and the absorption cross section of ²³⁸U; i.e., $w={\sigma }_{abs}/{\sigma }_{abs}^{238}$. We fix the lower and the upper limits at 10${}^3$ and 10${}^{-3}$ respectively. We adopt ${\sigma }_{abs}^{238}$ as a reference, since ²³⁸U has the highest concentration over fissile/fissionable nuclides, and provides a fundamental contribution to criticality and reaction rates for all burnup history.

Figure 8 reports the w values of the nuclides over the ZAI index. High weights are associated to actinides (ZAI > 90,000), where it is evident the reference weight of ${}^{238}$U corresponds to w = 10${}^0$ = 1. Moreover, important weights are shown for poisons with ZAI between 40,000 and 80,000, such as ${}^{135}$Xe (ZAI = 54,135) and ${}^{149}$Sm (ZAI = 62,149). For the majority of nuclide inventory, the weights are minimal at 10${}^{-3}$.

Here we employ a weighted POD, where the snapshots are composed of the variation of concentrations with respect to the case of having fresh fuel, multiplied by the weight

$$\underline{n_w^i}=({\underline{n}}^i-{\underline{n}}^{fresh})*w \mathrm{with} {\underline{n}}^i\ne {\underline{n}}^{fresh}, {10}^{-3}<w<{10}^3$$

(13)

$$\underline{\underline{S_{n_w}}}=[{\underline{n_w}}^1,{\underline{n_w}}^2,\dots ,{\underline{n_w}}^{1300}]$$

(14)

where $n_w$ is the weighted nuclide vector.

We report in Figure 9 the highest ten nuclides of the first five POD modes, obtained by SVD applied to $\underline{\underline{S_{n_w}}}$. As expected, the contribution of fissile and fissionable isotopes is increased with respect to the non-weighted POD, in which they are present in the first three positions for all basis functions, excepted for ${}^{232}$ Am in the basis function #5. The latter is generated by the decay of ${}^{243}$Pu and accumulates over time [43], and its presence is justified since it acts as a parasitic absorbent. The oxygen is not found in first position; indeed, its presence in non-weighted POD is only due to the high concentration in the fuel, while in weighted POD its weight is very low (w_O$=0.0033$).

In the online phase, we run criticality calculations with materials reconstructed from five weighted POD basis functions, for the same snapshots of the previous analysis. In Figure 10a, we show $\mathsf{\Delta }k$_eff obtained from weighted and non-weighted POD with five nuclide basis functions. For both cases, the approximations of k_eff are worst for the lowest (<10 MWd/kgU) and highest burnup (>70 MWd/kgU). We explain it by considering the first five basis functions mainly representative only of the most important fissile nuclides; i.e., ${}^{235}$U and Pu isotopes. Indeed, at low burnup, the presence of ${}^{239}$Pu, under formation, is negligible; at high burnup, the consumption of ${}^{235}$U provides a minor contribution to fission rate. At medium burnup, the presence of all fissile nuclides is more important, where $\mathsf{\Delta }$k_eff decreases for both sets of basis functions.

In Table 4, by excluding the highest and lowest burnup, $\mathsf{\Delta }$k_eff of weighted POD shows a minor oscillation around the average value, with a standard deviation of 77 pcm against 136 pcm of the non-weighted case. This value is comparable to that at 99 pcm, obtained by non-weighted POD with 10 basis functions. This suggests that weighted POD should be able to provide a better accuracy with minor number of basis functions than non-weighted POD.

Figure 10b reports the average relative variation (%), with respect to the full order cases, over densities of ${}^{235}$U, ${}^{238}$U and ${}^{239}$Pu. As expected, the approximation improves with weighted POD for all burnup points.

3.4. ROM of Burnup Matrices

POD Modes (Offline Phase)

After the presentation of results for $\underline{n}$ reconstruction, we analyze basis functions of $\underline{\underline{A}}$ and $\underline{\varphi }$. Figure 11 shows the sparsity pattern of the first POD mode obtained for burnup matrix, with distributions in logarithmic scale. The non zero elements are concentrated along the diagonal (absorption and decay rates) and fission product distributions on the right columns, as is typical for burnup matrices [29]. We only report the first mode because the other basis functions do not show noticeable differences in the structure. It follows that less energetic modes provide high order corrections on the reconstruction.

On the other hand, the differences over POD modes are evident for neutron fluxes. In this way, we are able to provide a specific physical interpretation of the POD modes, as reported in Figure 12, for the first four neutron flux basis functions, normalized in the lethargy scale. We observe that:

• Flux basis function $\#1$ represents the typical spectrum of a PWR reactor [44] with two distinct regions for thermal and fast neutrons;
• Flux basis function $\#2$ arises the contribution for the thermal part of the spectrum, representing the neutron thermalization;
• $\#3$ shows the resonances between the thermal and epithermal region; in particular, for ${}^{239}$Pu, ${}^{240}$Pu and ${}^{236}$U at 0.3, 1 and 5.46 eV respectively;
• Flux basis function $\#4$ represents the epithermal part of the spectrum, characterized by resonances of all fissile and fissionable nuclides.

For a more rigorous analysis, in Figure 13, we plot the eigenvalues obtained from POD of fluxes and burnup matrices. As expected, each trend shows the typical exponential decrease.

In order to evaluate if the projections of $\underline{\varphi }$ on POD basis functions provide suitable linear coefficients for $\underline{\underline{A}}$ reconstruction, we compare the coefficients of $\underline{\underline{V_{\mathsf{\Phi }}}}$ and $\underline{\underline{V_A}}$. We find average relative difference of 2% for the coefficients of the first basis functions and differences of hundreds of percentage points for less energetic basis functions.

For this reason, from Equation (10), we obtain the depletion $\underline{\varphi }$ basis functions, corresponding to the columns of $\underline{\underline{B_{\mathsf{\Phi }}}}$. The first four depletion basis functions, normalized in the lethargy scale, are shown in Figure 14. In comparison with POD modes, depletion basis function $\#1$ maintains the PWR spectrum shape while other basis functions change. However, they conserve part of their original physical meaning, since we distinguish thermal and fast regions, in addition to energy resonances of plutonium and the epithermal region.

3.5. Cases Description for ROM with Burnup Matrices (Online Phase)

We manage the online stage by a bash script that handles input and output, by following the flow chart on the right side of Figure 3. In each transport calculation, we run Serpent with 10 million active neutron histories. In burnup calculations, we solve the Bateman equations with MATLAB [45] through the solver ode15s [46].

At the beginning, we run ROM for the case of the snapshot generation with the highest power level, at 326 W/cm, by employing 5, 10, 12, 20, 25, 50 and 80 basis functions of $\underline{\underline{A}}$. In this way, we evaluate the different levels of approximation as the number of basis functions increases. We chose this case because it reaches the highest fuel consumption in the FOM, by covering the entire burnup range of the snapshots.

After that, we carry out FOM and ROM with 20 and 50 basis functions for three cases, not included in snapshot generation:

(a)

Change of total power at 200 W/cm (enrichment at 4.85%);

(b)

Change of enrichment at 5.5% (power at 326 W/cm);

(c)

Change of enrichment at 3.5% (power at 326 W/cm).

On one hand, the case (a) is an interpolation with respect to the space of the power parameter, which we only considered to obtain basis functions $\left\{{\mathsf{\Pi }}_i\right\}$. On the other hand, the cases (b) and (c) are an extrapolation of the space of the enrichment parameter. We discuss the results in the following sections.

3.5.1. ROM Case: Power at 326 W/cm, the Case of the Snapshot Generation

In order to evaluate the accuracy of the ROM, we compare criticality results with the FOM for the case of the snapshot generation at 326 W/cm.

Figure 15a shows the multiplication factor k_eff over time, obtained by FOM and ROM simulations for 5, 10, 20, 25, 50 and 80 basis functions. The trend of k_eff, which decreases under 0.9 at the end of the burnup history, is correctly reproduced by the ROM over time for all levels of approximation.

Figure 15b shows the difference $\mathsf{\Delta }$k_eff in pcm between FOM and ROM, to quantify the level of approximations of ROM with different number of basis functions. The simulations with 80 and 50 basis functions provide the best accuracy, with $\mathsf{\Delta }$k_eff between 200 and −200 pcm. Respecting the continuous Monte Carlo approach, variations of 200 pcm are also observed in the accuracy assessment of burnup analysis in [38]. Here, in addition to the validation with experimental data, code-to-code comparisons are carried out between the lattice codes TRITON/NEWT and Polaris against the Monte Carlo code KENO, which are modules contained in the SCALE Code System [6].

For 20 and 25 basis functions, $\mathsf{\Delta }$k_eff decreases over time, with values of −400 and −800 pcm at the end of the burnup. As expected, the case with five basis functions is less stable, with a fluctuation that reaches a maximum around 200 pcm and minimum at around −800 pcm. However, taking into account the code-to-code comparisons in [38,39,40,41,42], we can consider acceptable $\mathsf{\Delta }$k_eff in the order of hundreds of pcm.

Furthermore, the approximations of k_eff are explained for all cases by the variation of nuclide concentrations during the burnup. Figure 16 reports the relative percentage variation, with respect to the FOM, for the nuclide densities of ${}^{235}$U, ${}^{238}$U and ${}^{239}$Pu. The oscillation with five basis functions results, evidently, for all isotopes, where positive and negative variation of fissile increases and decreases the fission rates, leading to positive and negative variation of k_eff respectively. The opposite effect happens for ${}^{238}$U, where its concentration provides higher or lower contributions to radiactive capture.

By comparing the other cases for ${}^{239}$Pu and ${}^{235}$U, relative variations with 20 and 25 basis functions decrease over time, differently than 50 and 80 basis functions, where relative variations increase. Nevertheless, after 600 days, $\mathsf{\Delta }$k_eff remains negative for all cases, since the concentrations of ${}^{238}$U are higher than the ones calculated in the FOM.

For ${}^{235}$U, except for the case with five basis functions, the approximations of concentrations are acceptable. Indeed, relative variations until 1.5% at the end of the burnup history are in agreement with the work in [38,42]. In the latter, comparisons are carried out with Serpent over different burnup schemes, used both for lattice and Monte Carlo codes.

For ${}^{238}$U, the relative variations are acceptable since they are very low, under 0.15%.

For ${}^{239}$Pu, the relative variations are high for the initial steps at all approximations, reaching values between −10 and +5%. Indeed, at this time, the worst approximations are justified by the lower nuclide density of plutonium, which it is under formation. At higher burnup, as in [38,42], relative variations of ²³⁹Pu are under 3% with 50 and 80 basis functions, and slightly higher with 20 and 25 basis functions. Overall, we consider the results acceptable for basis functions not less than 20, since they all remain at the order of few percents.

Respect to the FOM, a higher amount of ${}^{238}$U and a lower amount ${}^{239}$Pu should be noted, from which we interpret both an under estimation of the breeding and an over estimation of the burning of plutonium, verified by checking the reaction rates in the burnup matrices.

For the sake of brevity, we do not report isotopes higher in the build-up chain, for which (in particular ${}^{242}$Pu, Am and Cm) the relative variations are under 6% for all cases, according to [38].

Cases with 10 and 12 Basis Functions

A separate discussion concerns the simulations with 10 and 12 basis functions. As reported in Figure 17, they show unphysical results for k_eff, which decrease exponentially after 300 days. At 730 days, $\mathsf{\Delta }$k_eff reaches negative values of thousands of pcm for both cases.

We point out that these unphysical results are not present with five basis functions. This is further shown by nuclide concentration of ${}^{235}$U, reported in Figure 18 for all cases analyzed. We observe an exponentially decrease until −2 and −4% for ${}^{235}$U with 10 and 12 basis functions, different than the stable trend with five basis functions. The increase of the number of basis functions does not entail a better prediction of the multiplication factor. This comes from the fact that the basis functions are built on an optimal decomposition of the burnup matrix reconstruction. This means that both the optimality and the monothonic decrease of the error can be taken for granted just in the case of burnup matrix reconstruction and not in general (for the k_eff estimation in this case). However, further investigations of this aspect are required in future.

In conclusion, considering the approximations of ROM both on k_eff and nuclide concentrations, the results are acceptable with a number of burnup matrix basis functions not less than 20, with very good agreement for 50 and 80 basis functions.

3.5.2. ROM of Case [a]: Interpolation of Power

After the previous analysis, it is important to provide verifications also for cases not investigated during the offline phase. For this reason, we perform FOM and ROM simulations with 20 and 50 basis functions of the TMI fuel cell at 200 W/cm, which represents an interpolation between the minimum and maximum power employed in the snapshot generation.

The results of criticality are reported in Figure 19. The trend of k_eff is reproduced correctly by the ROM. Moreover, as expected, a better accuracy is reached with 50 basis functions, obtaining absolute values of $\mathsf{\Delta }$k_eff less than 400 pcm at the end of the burnup history.

Considering the previous discussion in Section 3.5.1, variations of hundreds of pcm are acceptable for both cases, suggesting suitable burnup matrix basis functions for the parameter space of powers.

3.5.3. ROM of Case [b] and [c]: Extrapolation of Enrichments

The burnup of a fuel pin depends on many parameters, related to its composition and different operational conditions inside the reactors. It follows that, if we consider all space of parameters for this physical system, snapshot generation may be hugely computationally expensive. For this reason, we aim to verify whether it is possible to obtain an extrapolation with POD basis functions for different enrichments, not investigated during the offline phase. In particular, we simulate FOM and ROM with enrichments at 5.5% and 3.5%, to take into account higher and lower enrichment with respect the reference case at 4.85%. The number of basis functions are 20 and 50 for ROM.

Figure 20 and Figure 21 show the results for criticality values:

• For case (b), the trend of k_eff is well reproduced by 50 basis functions, with variation under 500 pcm for the entire burnup history. Unphysical results, after 400 days, are shown by 20 basis functions with variations until 3500 pcm with respect to the full order case.
• For case (c), $\mathsf{\Delta }$k_eff reaches −2000 pcm with 50 basis functions. We consider this result physical, since the trend of k_eff is correctly reproduced. With 20 basis functions, the k_eff becomes unphysical after 400 days, reaching a $\mathsf{\Delta }$k_eff of −14,000 pcm at the end of burnup.

Regarding case (a), as expected, we obtained the worst approximation to FOM. However, the case (b) shows better results than the case (c). We explain this fact by considering that the lower enrichment of the case (c) with respect to (b), with the same power at 326 W/cm, allows one to reach lower nuclide density of ${}^{235}$U at the end of the burnup history, with concentrations not covered during the offline phase. On the other hand, the case (b) reaches levels of nuclide density never covered by snapshot generation. In conclusion, we evaluate that the trend of k_eff is correctly reproduced with 50 basis functions for both cases.

4. Computational Cost

From the results of the present work, we evaluate the dimension of $\underline{\varphi }$ and the high-resolution energy spectrum ${\underline{\varphi }}_{HR}$, composed of 955,231 energy groups. These correspond to the unionized energy grid points, employed by Serpent for all reaction cross sections in its continuous-energy Monte Carlo neutron transport routine [47]. In burnup analysis, ${\underline{\varphi }}_{HR}$ is computed to average the transmutation cross sections for the calculation of the entries in $\underline{\underline{A}}$. In this way, the neutron transport calculation in ROM computes $\underline{\varphi }$, with 1968 energy groups, instead of ${\underline{\varphi }}_{HR}$, obtaining a gain of a factor of around 500. In this paper, the comparison is always made with the full order model; i.e., the Serpent model based on continuous Monte Carlo energy. Additional comparisons with other approaches as multi-group energy Monte Carlo or deterministic codes are demanded as future works.

Furthermore, the ROM skips the computation of each entry of $\underline{\underline{A}}$, by decreasing the number of calculations required. These aspects may become crucial with a higher number of burnup regions, at the order of hundreds of thousands.

In order to complete the analysis of the computational cost, we show in

Table 5. Hard disk and RAM memory usage for the snapshot matrices.

	n	$\underline{\underline{\mathit{A}}}$	$\mathit{\varphi }$
Number of snapshots	1325	1300	1300
Memory for snapshots matrix	16 MB	25,147 MB	20 MB
Number of total data	2.06$·{10}^6$	3.14$·{10}^9$	2.56$·{10}^6$

the memory requirement for the construction of snapshot matrices. In particular, this aspect becomes problematic for $\underline{\underline{A}}$ matrices, when we load them as not sparsely in MATLAB. For this reason, we carried out data analysis in the cluster of the Department of Energy at Politecnico di Milano, which supports 60 GB of RAM. To overcome the latter issue in the future, we suggest to load snapshots as ${\underline{\underline{A}}}_{sparse}$ matrices; i.e., stored with the nonzero elements and their indices. Table 5 reports that $\underline{\underline{S_A}}$ required around 25 GB of RAM, while $\underline{n}$ and $\underline{\varphi }$ required 16 and 20 MB respectively.

In future works, we shall divide the PWR fuel cell into a number of regions more minor than the basis functions employed in the ROMs; we aim to maintain tens of basis functions for the extension of the method to cases with 10${}^3$–10${}^6$ of burnup regions as well. In this way, it should be possible to save the memory usage for burnup matrices, also by orders of magnitude.

Finally, in order to save computational cost in snapshot generation for complex geometries, the offline stage can be limited to a single fuel cell. Indeed, the basis functions generated can be employed for the burnup of pins in fuel assembly and full core regions.

5. Conclusions

In this work, we developed a methodological approach, based on POD, to implement reduced order modeling in fuel burnup analysis.

The methodology was applied to the burnup of the TMI-1 fuel cell and was based in an offline/online phase, where the online phase is divided into two parts. In the first part, we reconstructed the material compositions by an optimal rank approximation of the snasphot matrix. At different levels of approximation with 5, 10 and 25 basis functions, k_eff was correctly reproduced, with differences of hundreds of pcm with respect to the FOM. In the second part, we carried out a burnup analysis, based on reconstruction of burnup matrix from POD basis functions. In reproducing k_eff and nuclide concentrations, taking into account code-to-code comparisons in literature, the results show good agreement with the FOM for 50 and 80 basis functions, and acceptable results for basis functions not less than 20.

An unphysical decrease of k_eff was observed for number of matrix basis functions between 5 and 20. This comes from the fact that the basis functions are built on an optimal decomposition of the burnup matrix reconstruction. It follows that the optimality and the monothonic decrease of the error can be taken for granted just in the case of burnup matrix reconstruction, and not in general. Further investigations of this aspect are needed in the future.

We also performed an extrapolation with POD basis functions for the enrichment parameter, at 5.5 and 3.5%, reproducing acceptable results for k_eff through 50 basis functions for both cases, with variations of hundreds of pcm.

From the analysis of the computational cost, the calculation of the flux with the ECCO 1968 group structure in ROM gains a factor of 500 in memory saving. In addition, the possibility to only compute the linear coefficients for fluxes in ROM allows one to decrease the number of calculations required.

In conclusion, this work has shown a new methodology with POD for fuel burnup analysis, tested on the TMI-1 unit cell for demonstration purposes. In future works, verifications and validations of the approach with computational and experimental benchmarks should be mandatory. For future applications in fuel assembly and full core design, the offline stage can be limited to a single fuel cell, by simulating the different operational conditions of interest.

Furthermore, since we proved that it is possible to reconstruct nuclide vectors from POD modes, it would be interesting to investigate the possibility to restrict the burnup analysis to tens of basis functions for fuel materials. This may allow us to further decrease the computational cost with respect to that needed for the total number of burnup regions in the full order model.