© G. K. Delipei et al., published by EDP Sciences, 2020

1 Introduction

The improvement of nuclear reactors computational modeling together with the increasing requirements in safety analyses lead to the development of Best Estimate Plus Uncertainty (BEPU) approaches [1]. BEPU is a systematic approach where Best Estimate (BE) codes are used to calculate safety criteria with an estimation of their respective uncertainties. BEPU is a step forward from conservative approaches allowing for better quantification of the safety margins. BE codes are considered the codes that take into account the most important underlying physical phenomena, governing the different studied scenarios, with reasonable approximations. Different BE codes usually model different physics. The correct modeling of transient scenarios can require the coupling between different BE codes due to the multi-physics interactions. This creates many challenges for the uncertainty quantification (UQ). With the term UQ we mean the estimation of the statistical distributions of the outputs of interest (in particular, their means and standard deviations) and the global sensitivity analysis of the outputs (in particular, their Shapley indices). Many code evaluations are required for the estimation of the statistical quantities and of the sensitivity indices rendering thus the computational expensive modeling the main limitation. A large scope of undergoing research, in both academia and industry, is the multifidelity High to Low modeling approaches that can reduce in an efficient way the computational cost while retaining a good predictive modeling.

The multifidelity approaches combine different sources of information coming from high fidelity models, low fidelity models and experiments. The fundamentals work of [2] and [3] created the theoretical framework that lead to papers such as [4] in structural dynamics, where a cheap surrogate model based on Polynomial Chaos was considered as a low fidelity model and corrected by a high fidelity code in order to accelerate the Bayesian calibration process using experimental data. In [5] a methodology for combining information from multifidelity codes and experiments in a Bayesian hierarchical model was developed using Gaussian processes. In [6] a novel machine learning based multifidelity approach was discussed in the context of Bayesian optimization. A comprehensive and detailed survey of multifidelity methods in uncertainty quantification, statistical inference and optimization can be found in [7]. In nuclear and more particularly in thermal-hydraulics the result of the PREMIUM project [8] indicated a strong user effect in the input uncertainty quantification. For this purpose in [9] a set of good practices is discussed regarding the input uncertainty quantification.

In this work we focus on the study of the Rod Ejection Accident (REA) in Pressurized Water Reactors (PWR) with a multi-physics coupling in an UQ framework. During the REA strong multi-physics coupling effects occur between neutronics, fuel-thermomechanics and thermal-hydraulics. This necessitates a multi-physics modeling to take into account the interdisciplinary interactions increasing significantly the computational cost. The uncertaintly quantification for REA becomes thus very challenging. In this coupling the fuel-thermomechanics modeling is the most expensive since it solves at the fuel rod scale the thermal conduction equations coupled with the thermo-mechanical ones at a fine discretization. A High to Low modeling is necessary to introduce and use a simplified model for fuel-thermomechanics reducing thecomputational cost. The high fidelity modeling manages to model complex thermo-mechanical and physico-chemical phenomena occurring during a REA. Some examples are the high burn-up structure, fission gas releases in the pellet-cladding gap and the pellet cladding mechanical interaction (PCMI). All these phenomena impact the fuel rod behavior through the pellet cladding gap heat transfer (H_gap). For this reason the low fidelity fuel-thermomechanics modeling usually does not model explicitly these phenomena but directly use a predefined H_gap value. There is, therefore, a strong interest in developing High to Low models for H_gap in REA. In this work, we develop and propose a methodology to develop such High to Low models for H_gap. The high fidelity fuel-thermomechanics code ALCYONE1 [10] is used to calibrate a simplified H_gap model basedon fuel thermal expansion. Compared to previous works such as [5], the low fidelity model is not a fitted surrogate but a physics based model that is calibrated in a preliminary way by minimizing the discrepancy with ALCYONE1 predictions. The model is then introduced into the CEA developed CORPUS [11] Best Estimate (BE) multi-physics tool, where core neutronics code APOLLO3^® [12] is coupled with core thermal-hydraulics code FLICA4 [13]. This creates an Improved Best Estimate (IBE) modeling of REA that is used for UQ. To the knowledge of the authors this is one of the first attempts to address H_gap High to Low modeling in transients. Some considerations have been discussed in the context of the UAM (Uncertainty Analysis in Modelling) benchmark [14]. The benchmark consists of three phases starting from stand-alone neutronics uncertainty analysis for steady state and reaching up to full coupling at the system level for transient scenarios. In this benchmark the H_gap is treated through lookup tables for the steady state and evolution calculations. However, there is not a dedicated treatment for the transient scenarios and this work and its prospects could be potentially integrated in the benchmark. Additionally in [15] the pellet cladding modeling was optimized including experimental data. The model’s uncertainty was estimated as well.

In Section 2 we detail the two available coupling frameworks at CEA for the modeling of REA. The computationally expensive Best Effort coupling and the cheaper BE coupling. In Section 3 the specifications of the case study are presented. The geometry of the studied PWR core and the selected nodalization are provided together with the REA characteristics and the inputs and outputs of interest. The High to Low H_gap modeling methodology is based on the results of an UQ analysis performed with the BE coupling (BE UQ). In Section 4 we present these results. In Section 5 we present the calibration of the High to Low Hgap model. The simplified H_gap model is detailed and the calibration results are presented. As mentioned above, the obtained H_gap model is introduced into the previous coupling creating the IBE. In Section 6 an UQ analysis is performed with the IBE coupling (IBE UQ) and the results are compared with the BE results of Section 4. Finally, we end this article with the main conclusions and perspectives of this work.

2 REA coupling framework

The REA is initiated by a control rod ejection due to mechanical malfunction inserting positive reactivity in the core. As a consequence power increases adiabatically until the fuel temperature starts increasing as well introducing a Doppler negative feedback that creates a power peak. The power then continues to decrease and when the heat flux reaches the moderator its density will start decreasing adding another negative feedback due to the negative moderator temperature effect. Depending on the core state at the moment of the ejection the power evolution and its damage on the first containment barrier (cladding) can vary significantly. During the transient strong multi-physics interaction effects occur between neutronics, fuel-thermomechanics and thermal-hydraulics necessitating a multi-physics coupling framework.

At CEA two different coupling schemes are established [16] for the modeling of REA based on CORPUS multi-physics tool as can be seen in Figure 1a. APOLLO3^® code is used for core neutronics, FLICA4 code for core thermal-hydraulics and ALCYONE1 code for pin fuel-thermomechanics. The Best Effort scheme involves the full coupling between the codes. At each time step the codes exchange 3D fields. APOLLO3^® uses the fuel Doppler temperature T_f from ALCYONE1 and moderator temperature T_m from FLICA4 to compute the power generated in the fuel P_f and the fluid P_l. FLICA4 uses P_l from APOLLO3^® and the cladding wall heat flux Φ_w from ALCYONE1 to compute theT_m and the external cladding wall temperature $T_c^{ext}$. ALCYONE1 uses the P_f from APOLLO3^® and the $T_c^{ext}$ from FLICA4 to compute the T_f and Φ_w. The Best Effort scheme models each discipline accurately but it is computationally very expensive for a full core REA and thus its use for the uncertainty propagation and global sensitivity analysis in the UQ framework is prohibited. For this reason the BE scheme is used where the coupling is reduced between APOLLO3^® and FLICA4. The complex fuel-thermomechanics phenomena during REA are modeled by the fuel-thermal module in FLICA4. In this module only the thermal aspects are modeled and all the other phenomena are taken into account through a constant H_gap during the transient. This simplification is very limiting for the UQ because it leads to the use of large H_gap uncertainties in order to penalize for the poor thermomechanics modeling. In this work we develop and propose an offline High to Low H_gap modeling methodology that alleviates the poor H_gap modeling issue during REA. This creates the IBE coupling scheme illustrated in Figure 1b.

Fig. 1 CEA coupling framework used for REA (a) and the proposed improvement in the BE scheme (b).

3 Case study

3.1 PWR geometry and modeling

The REA is studied in a large scale PWR core. The geometry has an 1∕8 symmetry illustrated in Figure 2. It consists of 193 fuel assemblies with UO₂ and UO₂ − GdO₃ fuel compositions. Two different types of control rods are inserted at different depths. The black rods (B) with high neutrons absorption that are typically used for the shutdown of the reactor and the grey rods with less neutrons absorption that are used in the day to day reactivity control. The core is at Hot Zero Power (HZP) conditions at the end of the cycle. Around the fuel assemblies there is one ring of water reflector assemblies. The total height of the core is 468.72 cm with a bottom and top reflector of 21 cm leading to a fuel active height of 426.72 cm. Each assembly is a 17 × 17 lattice of fuel pins with pitch 21.504 cm. The control rod that will be ejected initiating the REA is located on the periphery as highlighted in Figure 2. It is inserted 97 cm from the top into the fuel active region. Due to the extraction of the control rod there is 1∕2 symmetry for the REA.

The REA in the PWR core geometry is modeled using the BE and IBE coupling scheme discussed in Section 2. For core neutronics, APOLLO3^® code is usedwith a two group Diffusion approximation and void boundary conditions on the neutron current. The two group macroscopic cross-sections are parameterized in burn-up, boron concentration, moderator density and fuel temperature. The radial discretization is at the level of the quarter of assembly. For the axial discretization 34 meshes are used of which 30 are for the fuel active height and the rest for the top and bottom reflectors. For core thermal-hydraulics, FLICA4 is used with a 4 equations porous modeling and a multi-1D axial flow approximation. The system of 4 equations consists of: mixture mass balance, vapor mass balance, mixture momentum balance and mixture energy balance. The boundary conditions are the inlet mass flow and enthalpy and the outlet pressure. For the radial discretization one thermal-hydraulic channel is used for each assembly. For the axial discretization only the fuel active part is modeled using 30 meshes in accordance to APOLLO3^® modeling. One average fuel pin per thermal-hydraulic channel with 1D radial modeling is used in FLICA4 with a discretization of 25 radial regions for the fuel and 3 for the cladding. For the time discretization of the REA an adaptive end time is adopted based on the integral power evolution. For each transient in the UQ when the power surpasses half its nominal value then a SCRAM signal is sent. It is considered that from this time on 0.6 s are needed in order for the SCRAM to take place and end the modeling of the transient. The incremental time step is constant and equal to 0.001 s. The control rod is ejected in 0.1 s.

Fig. 2 PWR 1∕8 core geometry and characteristic dimensions. B indicates assemblies with black control rods, G with grey control rods and N with no control rods. The ejected control rod location is highlighted with red borders.

3.2 Initial state and reference transient

The PWR core is critical at the end of the cycle with HZP condition, meaning that the temperature is around 290 deg C and the power negligible (3.8W). Since the core is at the end of the cycle the boron concentration is quite low at 95.5 ppm. Other characteristic conditions of the core are provided in Table 1. There is a burn-up distribution resulting from core evolution calculation performed in a previous internal CEA study. The obtained radial burn-up distribution at the end of cycle is illustrated in Figure 3. The burn-up radially averaged at the level of the assembly ranges from 10 to 52 GWd∕t. The location of the ejected control rod is close to the periphery. In the core there is a radial and axial Xenon distribution that increases the control rod worth leading to an increased reactivity insertion and thus a violent prompt driven transient.

The reference (without uncertainties) REA characteristics obtained with the BE coupling are presented in Table 2 and Figure 4. The control rod worth is ρ_worth = 1.2 indicating aprompt driven transient. In the Figure we can observe the created power pulse of width Γ = 38 ms with a maximum power of ${\text{P}}_{\text{core}}^{\text{max}}=2.54{\text{P}}_{\text{nom}}$ at instant t^max = 292 ms. The nominal power is P_nom = 3800 MWth. Additionally, the F_xyz deformationfactor evolution in time is plotted. It starts from a value of 5 and reaches up to 25 when the control rod is fully ejected (0.1 s).

Fig. 3 PWR core burn-up radial distribution in the core.

Fig. 4 Integral power and Fxyz deformation factor evolution for the reference REA.

3.3 Input and outputs identification

Before presenting the BE UQ results we define the uncertain inputs that will be considered together with the outputs of interest. The inputs and outputs for the three disciplines are defined in Table 3. The statistical distributions of the inputs are presented in [17] and are included in Appendix A. In neutronics the two-group macroscopic cross-sections and the kinetic parameters are considered with a multivariate normal distribution. The rest of the inputs are considered independent of the neutronic inputs and between them. The thermal-hydraulic input distributions are mainly based on CEA experts opinions and are applied as random multiplicative factors with mean 1 on the different models. The UAM recommendations are used for the fuel-thermal inputs, i.e. for the thermal conductivities and specific heat capacities. The H_gap is the quantity of focus in this work. In the BE coupling scheme H_gap is an uncertain constant with a uniform distribution between the value for a complete open gap (2e³ Wm⁻²K⁻¹) and the one for a pellet-cladding contact (5e⁴ Wm⁻²K⁻¹). In the IBE H_gap is calculated by a model and thus the uncertainty will be carried by the model parameters. For the outputs 3 scalars and 1 functional output are considered. In neutronics we consider the maximum in time and space of the local linear power during the REA $P_{lin}^{max}$ and the radial linear power distribution $P_{lin}^{2D}$ at the time instance andaxial plane of $P_{lin}^{max}$. In fuel-thermomechanics we selected the maximum in time and space of the local stored enthalpy $H_f^{max}$ while in thermal-hydraulics the minimum in time of the distance to Departure from Nucleate Boiling (DNB) DNB^min. The latter output is defined as the difference between the Departure from Nucleate Boiling Ratio (DNBR) and the DNB threshold R_crit. The UQ methodology is detailed and tested in [17] and [18]. The methodology consists of an initial screening of the inputs and the use of kriging models [19] to perform uncertainty propagation and global sensitivity analysis using Shapley indices [20, 21]. More details about the methodology are discussed in Appendix B. The obtained results for the uncertainty propagation and global sensitivity analysis will be presented directly for both BE and IBE coupling schemes.

4 Best estimate uncertainty quantification

In this section we present the results of the BE UQ for the REA. The uncertainty propagation is performed by training kriging models between each scalar output and a reduced number of screened inputs. The screening is performed with a method based on HSIC statistical significance tests [22, 23] presented in [17]. The kriging models are trained on a Latin Hypercube Sampling (LHS) of size 250. For the functional output Principal Component Analysis (PCA) [24] is used and kriging models are trained for two principal components that represent 95% of the outputs variances. For all the quantities ofinterest the obtained Q² of the kriging models, evaluated on a separate LHS of size 125, is close to 1 with smallest value 0.98 for DNB^min and 0.945 for the second principal component of $P_{lin}^{2D}$. More detailed results are provided in Appendix B. The kriging models are used for brute force Monte Carlo uncertainty propagation with 10⁵ samples andthe results are presented in Figure 5. In each histogram we plot also the pdf of a normal distribution with the estimated mean and standard deviation (blue line) for a visual comparison of the normality. In Appendix B we present also the estimated relative standard deviations for the scalar quantities directly using the training LHS. It is important to notice that it is not recommended to use surrogate models for accurate estimation of the tails of the output pdf since they consist of rare events that the surrogate model has not been trained for.

We observe a quite large relative standard deviation for $P_{lin}^{max}$ of around 57%. The result for $P_{lin}^{2D}$ shows that the spatial distribution of the relative standard deviation does not vary significantly radially. For $H_f^{max}$ a 20% relative standard deviation is obtained and the distribution is approximately normal. For DNB^min a large mean value is obtained with large relative standard deviation of 57% resulting in a very small probability of reaching boiling crisis.

For the global sensitivity analysis, Shapley indices are estimated for the scalar outputs and aggregate Shapley indices for the functional ones as discussed in [17]. All the Shapley indices are estimated using the kriging models. For a more detailed focus on the Shapley indices estimation with surrogate models we refer to [25]. The R language [26] package “sensitivity” [27] (function “shapleyPermRand”) was used for the estimation of the Shapley indices and their confidence intervals, that quantify the uncertainty of the estimation method. The results are presented in Figure 6. The inputs in this Figure are the selected ones by the screening process of the UQ methodology. For the outputs $P_{lin}^{max}$, $P_{lin}^{2D}$ and $H_f^{max}$ the β_eff is responsible for 50% of the output variance with the other 50% attributed to the cross-sections TD₁, D₁, S_1→2 and IV₁. All these cross-sections are highly correlated to each other and thus it is difficult to distinguish their separate contributions. For the output DNB^min the H_gap is the dominant input responsible for 50% of the output variance while the remaining 50% is attributed to the cross-sections and the β_eff. This is due to the large H_gap uncertainty ranges and thus gives the incentive to improve its modeling.

Fig. 5 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin histograms and $P_{lin}^{2D}$ relative standard deviation distribution for BE UQ analysis.

Fig. 6 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin Shapley indices and $P_{lin}^{2D}$ aggregate Shapley indices for BE UQ analysis.

5 High to low pellet cladding gap heat transfer modeling methodology

5.1 Introduction

The gap heat transfer (H_gap) modeling in the BE coupling is performed through a constant value during the REA with a uniform distribution over a large interval. This is one of themost important modeling differences between Best Effort and BE coupling. There is a strong interest thus to improve the H_gap modeling and to introduce it into a refined coupling that we will call Improved Best Estimate (IBE). In this subsection we address this challenge by developing a High to Low H_gap modeling methodology. The methodology involves the calibration of a simplified H_gap model that is based on fuel thermal expansion. We consider that this model is adequate for the REA and especially for the gap closing phase. The calibration is performed through decoupled ALCYONE V1.4 calculations with imposed power evolution. In this Section we first detail the model with its calibration parameters. Afterwards, we present the methodology where we discuss issues such as how the power pulses are selected, how many H_gap models will be created for the different fuel assemblies etc. Finally, we present the calibration results for the PWR core.

5.2 Pellet cladding gap heat transfer simplified model

The sharp power increase in the REA leads to a corresponding sharp fuel temperature increase. We assume that the H_gap evolution is driven by the gap closing due to fuel thermal expansion and by the gas conductivity evolution in the gap. This is used to derive the simplified formulation for H_gap defined by equations (1a)–(1c). $$H_{gap}=\frac{{\lambda }_g(T_g,E_f)}{e(T_f),}$$(1a) $${\lambda }_g={\lambda }_{g,init}\left(1+{\theta }_1\frac{T_g-T_{g,init}}{T_{g,init}}+{\theta }_2{E}_f\right)$$(1b) $$e(T_f)=r_{c,init}-r_f(T_f)$$(1c)

where:

• λ_g is the gasconductivity in the gap and ${\lambda }_{g,init}=\frac{H_{gap}^{init}}{e_{init}}$ its initialvalue prior to the REA. The latter is calculated by the initial gap heat transfer $H_{gap}^{init}$ and initial gap width e_init.

• T_g is the gastemperature and T_g,init the initial gas temperature prior to the REA.

• E_f is the energy stored in the fuel during the REA.

• θ₁ and θ₂ are two calibration parameters that have to be estimated.

• e is the pellet-cladding gap width. It is assumed that only the fuel expansion is responsible for the gap evolution.

• r_c,init is the initial internal cladding radius prior to the REA.

• r_f = r_f,initα_f(T_f) is the fuel external radius. The fuel expansion is modeled using the fuel expansion coefficient α_f(T_f) in ALCYONE1 which is a cubic function of the fuel temperature T_f.

The H_gap predicted by the proposed simplified model is based on fuel thermal expansion and depends on the evolution of the λ_g and e. The λ_g is considered a linear function of T_g and E_f. The latter allows to include a historical effect on the H_gap. The two calibration parameters to be determined θ₁ and θ₂ are the coefficients of T_g and E_f respectively. The e evolution is assumed to depend only on the fuel thermal expansion while the cladding radius remains constant. For the modeling of the gas temperature T_g the average between the external fuel temperature and the internal cladding temperature is used.

5.3 Methodology

Having defined the simplified H_gap model the next step is to calibrate it. To this purpose we developed the High to Low methodology illustrated in Figure 7.

The starting point is the application of the UQ methodology in the BE coupling of Section 4. The temperature evolution in the fuel depends on the fuel assembly burn-up and the power seen by this assembly at its position in the core. Since we need to build a H_gap model for every fuel spatial mesh (1080 meshes) a clustering is necessary. In step 1.1, we cluster the assemblies by similar radial average burn-up. This means that we will construct one H_gap model for every identified cluster. While each cluster will have the same model, their different fuel initial conditions and the different power seen during the REA will lead to a 3D H_gap.

The next step 1.2 consists in selecting representative pulses from the BE UQ presented in Section 4. The pulses must cover most of the possible H_gap variations inside the cluster due to both statistical and spatial aspects. For the statistical aspect three pulses are extracted: namely the ones which have produced the mean, the 97.5% upper quantile, and the 2.5% lower quantile of $P_{lin}^{max}$. For the spatial aspects, one selects the assembly that achieves the maximal power for the pulses corresponding to the mean and upper 97.5% quantile and one selects the assembly that achieves the minimal power for the pulse corresponding to the lower 2.5% quantile. For the spatial aspects, when the mean and the upper 97.5% quantile are imposed then the assembly seeing the maximum power is selected. Correspondingly, when the lower quantile is imposed the assembly seeing the lowest power is selected. This creates for each H_gap model different representative axial and temporal profiles of linear power and external cladding temperature.

In step 2, the selected profiles are extracted and imposed in a decoupled ALCYONE1 REA transient calculation. One representative fuel pin is modeled. The resulting temperature and gap heat transfer profiles from ALCYONE1 are extracted and used for the H_gap model calibration. In step 3.1 the calibration is carried out by minimizing an objective function using the BFGS optimization method. The mean square error on the H_gap maximal and final values during the REA for each axial slice is considered as the objective function. Once the parameters that minimize this function are estimated the final step 3.2 is to quantify the calibrated model’s uncertainty. The two main sources of uncertainties are the initial conditions and the calibrated parameters. The former one is quantified as two multiplicative factors on the initial gap width and the initial H_gap with normal distribution with mean one and standard deviation 0.1. This is a result of a previous uncertainty propagation for fuel evolution calculations with ALCYONE1. The results also showed that the initial gap width and H_gap are almost fully negatively correlated. This leads to make the assumption that the two coefficients are fully negatively correlated (ρ = −1) rendering thus one effective uncertain quantity for the initial conditions H_g,i. The latter uncertainty source is the calibrated parameters. For simplifying their uncertainty quantification, they are considered as fully positively correlated (ρ = 1) with uniform distributions. The bounds of the distributions are calculated in order to account for the calibration error. The effective uncertain input representing the calibration uncertainty is H_g,m. Finally, the model of equations (1a)–(1c), including its two effective uncertain parameters, is introduced into the multi-physics coupling creating the IBE modeling of Figure 1b. It can be seen as an intermediate modeling between the BE and Best Effort modelings. Since the H_gap model is a simple analytic function the computational cost of IBE is similar to BE amounting to a significant gain compared to Best Effort coupling. The average Best Effort computational time for one REA is 3 h while the average BE and IBE is 6 min.

Fig. 7 High to Low Hgap modeling methodology scheme.

5.4 Results

The first step 1.1 consists in clustering assemblies with similar burn-ups. For each cluster one H_gap model will be considered. In the PWR core there is a 3D burn-up distribution. This leads to a total of 193 × 4 × 30∕2 = 11580 meshes (due to symmetry) with different H_gap evolutions due to different burn-ups and power histories. In order to avoid constructing one H_gap model for each mesh, the clustering of the assemblies of different burn-ups is carried out. At first we consider only radial burn-ups, by averaging the axial variations. Secondly, we observe that the burn-ups have radially small variations around three main values 15 GWd∕t, 30 GWd∕t and 45 GWd∕t due to the PWR fuel loading pattern. We, therefore, choose to cluster the assemblies based on these three values and we add one cluster for the minimum 10 GWd∕t and one cluster for the maximum burn-up 52 GWd∕tt. This is done in order to have models covering all the burn-up variations. Additionally, it could be potentially used in the future for an application of a full 3D burn-up dependent H_gap model by constructing models that interpolate the calibration parameters.

To summarize we consider a total of 5 fuel assembly clusters and for each cluster one H_gap model will be constructed. We have to select for each cluster representative boundary conditions that vary both randomly and spatially since each cluster includes different assemblies. This is performed in step 1.2. The selected boundary conditions are presented in Figure 8. We know that the REA is a local phenomenon located in the upper part of the core around the ejected control rod position as seen in the radial cross-section of the Figure. We thus expect large variations of H_gap on the upper part and low to negligible variations in the lower part. For the random aspects we use the results from Section 4. More specifically, we consider that $P_{lin}^{max}$ gives a good indicator of H_gap variations in REA. Based on this we select samples corresponding to the mean, the upper 97.5% and lower 2.5% quantiles of $P_{lin}^{max}$. From these samples we extract the linear power and cladding wall temperature axial and temporal evolutions. We combineboth random and spatial aspects by selecting representative assemblies in the upper part for the mean and upper quantiles while we select their mirror assemblies from the lower part. The selected assemblies are presented in Figure 8 where each assembly has the burn-up value of its cluster. The green circles correspond to the selection for the mean and upper quantile while the yellow circles for the lower quantile. For the 10 GWd∕t and 52 GWd∕t clusters there is only one possible assembly for each cluster. For the other clusters, from the many available options we prefer the assemblies close to the ejected control rod location. For 15 GWd∕t and 30 GWd∕t clusters we select two different assemblies while for the 45 GWd∕t we select three, since we consider that this cluster will have the largest H_gap variations due to its high burn-up.

In step 2 the extracted boundary conditions are imposed in ALCYONE1 and the REA stand-alone calculations are performed for each cluster. A representative fuel pin for each selected assembly in each sampling is modeled with the same axial discretization as the BE modeling. The H_gap and radial temperatureaxial and temporal evolutions during the REA are computed. These results together with the stored fuel energy are extracted and used for the H_gap model calibration of each cluster. Only the values corresponding to the time of the maximum and last value of H_gap during the REA for all the axial slices are kept for the calibration. The created dataset size varies for the different clusters depending on the number of representative assemblies. The results will be presented for 15GWd∕t and 30 GWd∕t and 45GWd∕t clusters, because they are the ones around the control rod ejection location. For 15 GWd∕t and 30 GWd∕t cluster the dataset size is: 30 (axial slices) × 3 (quantiles and mean) × 2 (assembly locations) × 2 (H_gap values during REA) = 360. For 45 GWd∕t cluster the size is: 30 (axial slices) × 3 (quantiles and mean) × 3 (assembly locations) × 2 (H_gap values during REA) = 540. The calibrationparameters of the H_gap models are estimated by minimizing the mean square error on these datasets. The resulting calibration errors for the three different models are presented in Figure 9 and the estimated values of the calibration parameters in Table 4. We observe that the calibration errors increase in general with burn-up. The maximum error is of the order of 8% for the 45 GWd∕t cluster.

For the uncertainty quantification of the calibration parameters θ₁, θ₂ we make the assumption that they are positively fully correlated, which makes it possible to simplify significantly their uncertainty quantification. Additionally, as for all the other inputs they are also fully correlated spatially. This means that the calibration parameters of each cluster vary homogeneously. We assign uniform distribution to each calibration parameter with ranges that cover the calibration errors (the uniform distribution maximizes the entropy amongst distributions with prescribed ranges). The results for the estimated ranges for each parameter are shown in Table 5.

The bounds of the uniform distributions are used for the prediction of the H_gap evolution during REA by the calibrated models. The results are compared to the ALCYONE1 calculations and are illustrated in Figures 10, 11 and 12. For the comparison the axial slice with the maximum H_gap value for the three main clusters is presented. The plotted H_gap predictions are also the ones with the largest errors and we can see that for all the predictions the ALCYONE1 calculation is within the uncertainty bounds created by the calibration parameters.

Finally, there are two different sources of uncertainties for the calibrated models. The first one is due to the calibration error and is quantified by the distributions of the calibration parameters. The parameters are assumed fully positively correlated. The second one as discussed previously is due to the initial conditions. At the end, two new uncertain parameters replace the constant H_gap of the BE modeling. The first one is related to the calibration error of the models H_g,m and the other is related to the initialconditions H_g,i.

Fig. 8 Selected assemblies on the symmetric 1∕2 PWR geometry for each cluster. The green circles indicate the selection of the mean and upper quantiles and the yellow circles indicates the selection for the lower quantile.

Fig. 9 Hgap model calibration errors for the different assembly clusters where index is the label of the data used for the calculations.

Fig. 10 15 GWd∕t Hgap model prediction with its uncertainty bounds in green.

Fig. 11 30 GWd∕t Hgap model prediction with its uncertainty bounds in green.

Fig. 12 45 GWd∕t Hgap model prediction with its uncertainty bounds in green.

6 Improved best estimate uncertainty quantification

The calibrated H_gap models from Section 5 are introduced into the BE coupling to create an intermediate IBE multi-physics coupling. The complete UQ methodology detailed in [17] is applied on this improved modeling. It is important to mention that the computational cost does not increase. As in the BE UQ, kriging models were trained on a LHS of size 250 and used in brute force Monte Carlo for uncertainty propagation with 10⁵ samples for each output. The results are presented in Figure 13. The validation of the kriging model’s training show a Q² close to 1 with the smallest being 0.987 for DNB^min and 0.945 for the second principal component of $P_{lin}^{2D}$. More detailed results are provided in Appendix B.

The obtained histogram for $H_f^{max}$ is approximately normal as for the BE UQ. We observe small increase of 2% in the mean value and a reduced relative standard deviation from 20% to 16%. For $P_{lin}^{max}$ the mean value significantly decreases by 10% with similar relative standard deviation. The relative standard deviation of $P_{lin}^{2D}$ is not affected by the improved H_gap model. The output quantity that is the most impacted is the DNB^min with an increase of 14% of the mean value and a decrease of the relative standard deviation from 57% to 42%. This means that there is significantly smaller probability to reach boiling crisis. In Appendix B we present also the estimated relative standard deviations of the scalar quantities directly using the training LHS. Although there are some discrepancies due to the limited number of samples, the same trends are observed as well.

The impact on the different mean values is attributed to the more realistic modeling of the H_gap evolution during the REA. In the BE modeling the mean constant value of H_gap is 2.4e⁴W/m^2K, much larger than the one predicted by the calibrated models, and it is applied during the whole duration of the transient. This leads to more heat transferred from the fuel to the coolant. The fuel temperatures are lower with a corresponding weaker Doppler feedback and thus a higher maximum linear power. The increased heat extracted from the coolant in the BE UQ explains also the smaller minimum distance to boiling crisis compared to the IBE UQ. The lower fuel temperatures induce also the observed lower stored enthalpy.

The global sensitivity analysis is performed by estimating Shapley indices for all the outputs. The results are presented in Figure 14. For $P_{lin}^{max}$, $P_{lin}^{2D}$ and $H_f^{max}$ as in the BE modeling the β_eff is responsible for around 50% of the outputs variances and the TD₁ and D₁ are responsible for the remaining 50%. A significant difference in observed for the DNB^min sensitivities. The gap heat transfer is not any more the dominant input, instead as for the other outputs the β_eff and the TD₁ and D₁ account for most of the DNB^min variations. This explains also the significant reduction of the DNB^min relative standard deviation.

Finally, we studied the 3D H_gap ($H_{gap}^{3D}$) as an output of interest and included it in the uncertainty quantification methodology. For $H_{gap}^{3D}$ two principal components are needed to represent 95% of its variance. The estimated means and standard deviations for the radial and axial cross-sections at the location and instant of the local maximum are presented in Figures 15, 16. As expected by the calibrated H_gap models, the assemblies with higher burn-up have also higher H_gap. The maximum value is obtained at the assemblies on the right and left of the assembly with the ejected control rod. This is due to the important power seen by these assemblies in combination with their high burn-up. The relative standard deviation distribution exhibits strong variations from 10% up to 32%. The 10% lower bound on the H_gap is related to the initial conditions uncertainties.

The aggregate Shapley indices for $H_{gap}^{3D}$ are estimated and presented in Figure 17. In this case the H_g,i is the dominant input responsible for 80% of the outputs variances. The remaining 20% is mainly explained by β_eff. This result is not surprising since the initial conditions determine the H_gap evolution.

Fig. 13 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin histograms and $P_{lin}^{2D}$ relative standard deviation distribution for IBE UQ analysis.

Fig. 14 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin Shapley indices and $P_{lin}^{2D}$ aggregate Shapley indices for IBE UQ analysis.

Fig. 15 $H_{gap}^{3D}$ estimated mean and relative standard deviation in the axial cross-section for IBE UQ analysis.

Fig. 16 $H_{gap}^{3D}$ estimated mean and relative standard deviation in the radial cross-section for IBE UQ analysis.

Fig. 17 $H_{gap}^{3D}$ Shapley indices for IBE UQ analysis.

7 Conclusion and perspectives

In this work we have presented a High to Low H_gap model methodology for REA uncertainty quantification in a PWR core. A simplified H_gap model based on fuel thermal expansion has been considered. The CORPUS multi-physics coupling framework is used for thetransient modeling. The fuel-thermomechanics code ALCYONE1 is used as the high fidelity code. The available Best Estimate (BE) coupling scheme involves a coupling between core neutronics code APOLLO3^® and core thermal-hydraulics code FLICA4. For APOLLO3^® a two-group diffusion modeling is used while for FLICA4 an axial multi-channel 1D modeling. The internal fuel thermal modeling of FLICA4 is updated to include the H_gap predicted by the calibrated model. This creates an Improved Best Estimate (IBE) coupling scheme.

The calibration methodology involves different steps. The starting point is a previous BE REA uncertainty propagation that creates the set of available results from which representative ones will be selected for the calibration. Then the first step is to cluster the assemblies and use one H_gap model for each cluster. The assemblies are clustered based on their burn-up, which creates 3 main clusters: 15, 30 and 45 GWd∕t. The next step is to select the fuel-thermomechanics boundary conditions on which ALCYONE1 is run for each cluster. REA power pulses are selected from the available BE results in such a way to cover the extreme cases and in locations where the higher and lower power variations are observed. The ALCYONE1 calculations are performed for each cluster power pulses and the H_gap and temperatures evolutions are extracted. The H_gap model calibration is performed for all the axial slices at two time instances: the instance of the maximum H_gap during the REA and the last value of H_gap. The calibration parameters are estimated for these datasets and the calibration error for all the time steps is computed. The assumption has been made that the calibration parameters are fully positively correlated. Their uncertainty bounds are selected in order to create H_gap intervals that cover the ALCYONE1 results. The results indicate that the uncertainty of the parameters increases with burn-up, which shows the limitation of such a simplified model for high burn-up fuel rods. Additionally, an uncertainty of 10% on the initial rod conditions is considered as a result of ALCYONE1 evolution calculations. Finally, the last step is to include the calibrated H_gap model into the BE coupling and to create the IBE.

An uncertainty quantification study is performed using the IBE and compared to the BE results. The most significant impact is on the distance to boiling crisis with an increase of 14% in its mean value reducing the probability of boiling crisis. The relative standard deviation of DNB^min is reduced from 57% to 42%. The reason is the decreased sensitivity on the H_gap due to its better uncertainty quantification and modeling. This is observed in the estimated Shapley indices, where the H_gap sensitivity reduces from 50% to almost 0. These results, although in an academic proof of concept study, show the potential gain of High to Low approaches for H_gap modeling in REA uncertainty quantification.

Many future prospects open following this work. On the one hand a more elaborated H_gap model could be sought with incorporation of the burn-up parameter to avoid creating one model for each burn-up. Besides that, a more direct High to Low modeling could be investigated in order to inform specific low level code models, such as simplified dynamical H_gap models, from the high level code. On the other hand, different transient scenarios and different reactor types could be studied in order to cover a larger spectrum of applications since they all share the common need for High to Low H_gap modeling.

Acknowledgements

APOLLO3^® is registered trademark of CEA. We gratefully acknowledge CEA, Framatome, and EDF for their long-term partnership and support. We would like also to thank the APOLLO3^® and FLICA4 development teams for their support and successful exchanges. This work would not be possible without the continuous help of A. Boulore, E. Royer, J.-M. Martinez and C. Patricot.

Author contribution statement

G. Delipei has performed the calculations and written the manuscript. J. Garnier, J.-C. Le Pallec, and B. Normand have performed scientific advisory roles and provided support through expert judgement, critical verification, and proofreading of the manuscript.

Appendix A Input uncertainty quantification

The input uncertainty quantification is performed by using the UAM recommendations when they are available and by using expert opinions for the rest. The resulting pdf for the inputs are presented in Figure A.1.

Fig. A.1 Input uncertainty quantification, where U(a,b) is a uniform distribution over (a,b) and N(a,b) is a normal distribution with mean a and standard deviation b. Pr is the sum of two independent random variables.

A multivariate normal distribution is used for the neutronic inputs. The mean vector of the pdf is the reference cross-sections produced at CEA Σ_CEA and the covariance matrix C_UAM is empirically estimated by the results of UAM. The UAM provides a dataset of two-group macroscopic cross-sections obtained by neutronic lattice uncertainty propagation simulations. This dataset of 100 realizations was adapted to the cross-sections used by APOLLO3^® by assuming:

• 1.

  Negligible uncertainties on the up-scattering cross-section S_2→1.

• 2.

  Negligible n − 2n, n − 3n… cross-sections uncertainties.

The adapted macroscopic cross-sections dataset is used to estimate the correlation matrix and the relative standard deviations of the cross-sections. Finally, C_UAM is calculatedusing the correlation matrix, the relative standard deviations and the reference CEA cross-sections. The UAM correlation matrix is illustrated in Figure A.1, where we can observe large positive and negative correlations. The TD₁ is strongly positively correlated with S_1→2 and IV₁ and negatively correlated with D₁. The β_eff is also strongly positively correlated with λ_eff.

The rest of the inputs are considered independent of the neutronic inputs and between them. The thermal-hydraulic inputs are mainly based on CEA experts opinions and are applied as multiplicative coefficients on the different parameters. The R_crit usually is penalized to 1.3 and in this work we consider this value as the 95% upper quantile of a normal distribution with mean value 1 and standard deviation 0.15. For H_DNB a uniform distribution is used reflecting the limited current knowledge about this phenomenon. In fuel-thermomechanics inputs and more specifically for the thermal conductivities and specific heat capacities the UAM recommendations are used. For the Rowlands temperature a uniform distribution is considered on the weight fraction of the fuel centerline temperature. In the reference situation the Rowlands temperature has a 4∕9 weight on the fuel centerline temperature and 5∕9 on the fuel external surface temperature. By using an uncertain multiplicative factor with $\mathcal{U}(0,1)$ distribution on the centerline temperature weight we consider that it can only decrease uniformly between 4∕9 and zero with a corresponding increase in the external surface temperature weight. The H_gap is a particular input quantity since in the BE coupling it is an uncertain constant in FLICA4 with uniform pdf bounded from below by its value for a complete open gap (2000 Wm⁻²K⁻¹) and from above by its value for a pellet-cladding contact (50000 Wm⁻²K⁻¹). In the IBE coupling H_gap is a result and thus is not considered, instead the calibration parameter’s uncertainty and the initial conditions of the model are used. Finally, the power radial profile uncertainty is modeled by a multiplicative factor on the fuel external surface power with a pdf resulting from the convolution of a normal and a uniform distribution. The power radial profile is peaked towards the periphery for the high burn-up fuel pins. The deformation is increasing with burn-up and to model it an explicit function of burn-up is used. The uncertainty is a convolution of two effects: the uncertainty of the function used $\mathcal{N}(0,0.0175)$ and the uncertainty due to the presence or not of a guide tube near the fuel pin $\mathcal{U}(0,0.04)$. The latter effect is a result of a previous CEA study.

Appendix B Uncertainty quantification methodology additional results

In Section 4 we apply the UQ methodology developed in a previous work [17] on the BE coupling. The methodology consists of an initial screening of the inputs based on HSIC statistical significance tests, a training of a surrogate kriging model between each output and the reduced inputs and the use of the kriging models to perform uncertainty propagation and global sensitivity analysis. The computational cost of this modeling is 3 h in average. In this Appendix we provide additional results regarding the identified reduced input space and the validation of the kriging models. Concerning the screening process, a random sampling of size 125 is used as a Design of Experiments (DOE). The result for the identified important input parameters for each output of interest are gathered in Table B.1. They and can be grouped into two subspaces: I₁ = (TD₁, NF₂, D₁,S_1→2,IV₁,IV₂, β_eff,Cp_f H_gap, T_R) and I₂ = (TD₁, NF₂, D₁,S_1→2,IV₁, β_eff, H_gap, H_c).

Concerning the training of the kriging models a LHS of size 250 optimized in both the complete input space and the two important subspaces is used as DOE. The kriging models are trained on the identified input subspaces. The resulting Q² includes thus the dimension reduction error. The R², the leave-one-out cross validation $R_{loo}^2$ and the Q² are presented in Table B.2. The Q² is estimated on an independent LHS of size 125. For all the outputs the Q² is close to 1 with smallest values for DNB^min (0.98%) and the second principal component of $P_{lin}^{2D}$ (0.945%).

The same methodology is applied afterwards in the IBE coupling, where the calibrated H_gap models are introduced into the BE coupling. It is important to mention that the computational cost does not increase. The uncertain inputs and outputs of the previous study are used with the replacement of the constant H_gap uncertain input by the H_gap model uncertain parameters. The two new uncertain inputs are H_g,m and H_g,i related to the H_gap model calibration error and initial conditions. For the screening, as before, a random sampling of size 125 is used as DOE. The result for the identified important input parameters are gathered in Table B.3. They can be grouped into two subspaces: I₁ = (TD₁, NF₂, D₁, β_eff, Cp_f,H_g,i, T_R) and I₂ = (TD₁, NF₂, D₁, β_eff, H_g,i, H_c, R_crit, T_R).

Compared to the BE study we observe the inclusion of the parameter H_g,i in both subspaces while the parameter H_g,m is rejected. This means that for the outputs of interest the initial conditions are more important than the calibration parameters uncertainties. For the kriging models a training LHS of size 250 with optimized subspaces I₁ and I₂ is constructed. The result for the R², $R_{loo}^2$ and Q² are presented in Table B.4. The Q² is estimated on an independent LHS of size 125. For all the outputs the Q² is close to 1 with smallest values for DNB^min (0.987) and the second principal component of $P_{lin}^{2D}$ (0.945).

Finally, we present in Table B.5 the estimated relative standard deviations using directly the training LHS code evaluations. We observe some differences compared to the kriging model’s brute force Monte Carlo due to the limited number of samples. However, similar trends are observed between the BE and IBE study. The most noticeable is the decrease of the DNB^min relative standard deviation from 50.32% to 36.63%.

References

All Tables

All Figures

	Fig. 1 CEA coupling framework used for REA (a) and the proposed improvement in the BE scheme (b).
In the text

	Fig. 2 PWR 1∕8 core geometry and characteristic dimensions. B indicates assemblies with black control rods, G with grey control rods and N with no control rods. The ejected control rod location is highlighted with red borders.
In the text

	Fig. 3 PWR core burn-up radial distribution in the core.
In the text

	Fig. 4 Integral power and Fxyz deformation factor evolution for the reference REA.
In the text

	Fig. 5 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin histograms and $P_{lin}^{2D}$ relative standard deviation distribution for BE UQ analysis.
In the text

	Fig. 6 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin Shapley indices and $P_{lin}^{2D}$ aggregate Shapley indices for BE UQ analysis.
In the text

	Fig. 7 High to Low Hgap modeling methodology scheme.
In the text

	Fig. 8 Selected assemblies on the symmetric 1∕2 PWR geometry for each cluster. The green circles indicate the selection of the mean and upper quantiles and the yellow circles indicates the selection for the lower quantile.
In the text

	Fig. 9 Hgap model calibration errors for the different assembly clusters where index is the label of the data used for the calculations.
In the text

	Fig. 10 15 GWd∕t Hgap model prediction with its uncertainty bounds in green.
In the text

	Fig. 11 30 GWd∕t Hgap model prediction with its uncertainty bounds in green.
In the text

	Fig. 12 45 GWd∕t Hgap model prediction with its uncertainty bounds in green.
In the text

	Fig. 13 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin histograms and $P_{lin}^{2D}$ relative standard deviation distribution for IBE UQ analysis.
In the text

	Fig. 14 $P_{lin}^{max}$, $H_f^{max}$ and DNBmin Shapley indices and $P_{lin}^{2D}$ aggregate Shapley indices for IBE UQ analysis.
In the text

	Fig. 15 $H_{gap}^{3D}$ estimated mean and relative standard deviation in the axial cross-section for IBE UQ analysis.
In the text

	Fig. 16 $H_{gap}^{3D}$ estimated mean and relative standard deviation in the radial cross-section for IBE UQ analysis.
In the text

	Fig. 17 $H_{gap}^{3D}$ Shapley indices for IBE UQ analysis.
In the text

	Fig. A.1 Input uncertainty quantification, where U(a,b) is a uniform distribution over (a,b) and N(a,b) is a normal distribution with mean a and standard deviation b. Pr is the sum of two independent random variables.
In the text