Preliminary uncertainty and sensitivity analysis of the Molten Salt Fast Reactor steady-state using a Polynomial Chaos

Preliminary uncertainty and sensitivity analysis of the Molten Salt Fast Reactor steady-state using a Polynomial Chaos Expansion method / Santanoceto, Mario; Tiberga, Marco; Perkó, Zoltán; Dulla, Sandra; Lathouwers, Danny. In: ANNALS OF NUCLEAR ENERGY. ISSN 0306-4549. 159(2021), p. 108311. [10.1016/j.anucene.2021.108311] Original Preliminary uncertainty and sensitivity analysis of the Molten Salt Fast Reactor steady-state using a Polynomial Chaos Expansion method


Introduction
Thanks to its low-carbon energy production, the nuclear sector can play a vital role in accommodating the world's rising energy demand while fighting climate change (Brook et al., 2014). However, over the last few decades concerns of accidents, proliferation, and waste management have led to a general public and political opposition to the technology and the construction of new power plants, and subsequently the share of nuclear power in the global electricity supply has decreased (International Energy Agency, 2019). In order to address these issues, in 2001 the Generation IV International Forum (GIF) proposed a roadmap aimed at designing Generation IV (Gen-IV) nuclear reactors, which excel in sustainability, safety, reliability, and proliferation resistance, and are economically attractive (Generation IV International Forum, 2002).
Among these, the Molten Salt Fast Reactor (MSFR) (Allibert et al., 2016;Gerardin et al., 2017) offers many advantages thanks to the adoption of a liquid fuel (LeBlanc, 2010;Serp and Allibert, 2014). The high operating temperatures and the excellent salt thermal properties guarantee high thermodynamic performances, while the low operating pressure and the strongly negative temperature feedback coefficient improve stability and safety.
However, due to the lack of know-how and operational experience with the technology, the development of the MSFR design needs extensive research from a number of viewpoints. To aid these efforts, Sensitivity Analysis (SA) and Uncertainty Quantification (UQ) play a vital role by evaluating how uncertain physical/chemical properties or unknown design specifications impact key operational and safety parameters. The need for accurate SA and UQ is further emphasized by the fact that design and licensing procedures for nuclear power plants increasingly rely on Best Estimate Plus Uncertainty approaches, applying the most precise and complete computational system models accompanied by thorough evaluation of the computed results' variability and sensitivity (Wilson, 2013).
Assuming an uncertain set of data, models and tools, UQ estimates the statistical information of system responses. This work focuses only on aleatory variables related to data, assuming that parameters have an accurate, but unknown value. The propagation of aleatory input data uncertainties to the outputs (i.e., the system responses) can be performed in multiple ways. Historically, UQ has been carried out using Monte Carlo (MC) methods or Latin Hypercube Sampling, which consists of generating a sufficiently large number of input parameter realizations according to their known (or assumed) statistical properties, followed by evaluating the corresponding system responses using the computational model of the system (Le Maitre and Knio, 2010). Finally, the statistical information of the outputs, such as means, standard deviations, or the complete probability density functions (pdfs), is obtained using appropriate estimators based on the large number of calculated responses. However, since the design of novel nuclear reactors requires accurate high-fidelity simulations that are typically computationally expensive, traditional sampling based methodologies are infeasible in these situations (Sudret, 2008). As a result, novel methodologies have been developed to mitigate the computational burden, such as stochastic collocation (Sankaran and Marsden, 2011), response surface approximation, or spectral projection (Cavallo and Sinha, 2007;Shih and Williams, 2009;Roy and Oberkampf, 2011).
The Polynomial Chaos Expansion (PCE) method adopted in this work belongs to the latter class. It consists of constructing a polynomial approximation of the responses as function of the stochastic inputs (Le Maitre and Knio, 2010), and it is potentially capable of strongly reducing the number of model evaluations with respect to MC analyses. Moreover, we opted for a non-intrusive approach, which treats the reference model as a ''black-box", thus requiring no modifications of the underlying equations. These characteristics make the adopted PCE approach particularly attractive for nuclear applications (see, for example, Fichtl and Prinja, 2011;Gilli et al., 2012Gilli et al., , 2013Perkó et al., 2014b;Turati et al., 2018;Peng et al., 2019;Skarbeli and Álvarez-Velarde, 2020), where the reference models are typically coupled multi-physics simulations often implemented in difficult to manipulate legacy codes and -in general -are computationally expensive.
The published literature on UQ/SA approaches applied to molten salt reactor systems is limited. Monte-Carlo based UQ analyses were performed for solid-fueled fluoride-salt cooled reactor designs, but the studies were limited to thermal hydraulics inputs, without taking into account neutronics uncertainties (Wang et al., 2017;Yang et al., 2017;Romatoski and Hu, 2019). The propagation of nuclear data uncertainty to the effective multiplication factor of the MSFR was preliminarily studied by Abrate et al. (2019) using an extended generalized perturbation theory approach implemented in the MC neutron transport code Serpent 2 ; however, the effects of temperature feedback and fuel motion were neglected. Coupling effects between neutronics and thermal-hydraulics were taken into account in the adjoint-based sensitivity analysis of the MSFR steady-state by Jeong et al. (2020); however, strongly simplifying assumptions were adopted for the reactor model (e.g, 1D geometry, constant mass flow rate, single-group diffusion). Finally, a non-intrusive reduced order model approach was recently proposed and applied for UQ and SA of both a simplified liquid-fueled reactor system (Alsayyari et al., 2020) and the full MSFR.
In this paper, we present the use of PCE to perform uncertainty quantification and sensitivity analyses of complex, multi-physics molten salt reactor models, assuming that both neutronics and thermal-hydraulic input data are stochastic. The approach -whose validity was proven in an earlier work of the authors  where a simplified liquid-fuel nuclear reactor was considered -is here extended to analyze the MSFR under steady-state conditions. The considered MSFR design is that investigated within the H2020 SAMOFAR project (http://samofar.eu/), while its reference model is constituted by a three-dimensional, multi-physics solver recently developed at Delft University of Technology (Tiberga et al., 2020b). Given the high computational cost of the model, the use of PCE techniques -able to minimize the number of simulations -is crucial to make UQ analysis feasible. We present preliminary statistical information on key neutronics and thermalhydraulics design aspects and rank the stochastic input parameters according to their contribution to the responses' variance, thus providing useful indications for future research efforts.
The reminder of the paper is outlined as follows. 2 presents the adopted PCE method, providing details on the mathematical formulation of the problem and the post processing step necessary to derive the sensitivity indices for the system responses. Section 3 is devoted to the description of the MSFR design under investigation along with the neutronics and thermal-hydraulics simulation tool used as reference model. Section 4 lists the considered system responses and the sources of uncertainties. Section 5 describes the procedures adopted to reduce the computational cost of the analysis. 6 reports and discusses the results of our work. Finally, 7 draws the main conclusions and gives suggestions for future research.

Polynomial Chaos Expansion method
Polynomial Chaos Expansion (PCE) is a technique based on the approximation of a response with a set of polynomials up to an appropriate order. After having determined each expansion coefficient by performing sufficient model evaluations, the obtained expression contains the response's analytical dependence on the stochastic inputs. Consequently, the complete statistical information of the response (such as mean, standard deviation, full pdf) can be obtained by random sampling of the cheap polynomial model rather than the full model itself, requiring only seconds for even millions of samples. PCE was originally developed by Wiener (Wiener, 1938) with the application of Hermite polynomials suitable for normally distributed parameters alone, but later was generalized as generalized Polynomial Chaos (gPC) to many other distributions, using different kinds of polynomial (Xiu and Karniadakis, 2002). PCE offers the possibility to greatly reduce the number of model evaluations compared to MC analyses. However, widening the input parameter space increases the number of simulations rapidly, which can become comparable with conventional analyses. This issue is known as ''curse of dimensionality".
In the following, we present a short theoretical background of the PCE method adopted in this work, which is implemented in the OpenGPC code developed at Delft University of Technology (Perkó et al., 2014a,b;Perkó, 2020). This tool is able to build PCE approximation for a large number of system responses and uses Smolyak sparse grids to alleviate the ''curse of dimensionality".

Problem definition
We define the probability space ðH; R; PÞ, where H is the sample space, containing the set of h 2 H, random events, R is the ralgebra (collection of subsets of H), and P is the probability measure. Stochastic inputs parameterizing the random events are nðhÞ 2 !, where ! is the support of the joint probability density function p n ðnÞ. The present work focuses only on independent random variables, so the joint pdf is where N is the number of stochastic inputs and n ¼ ðn 1 ; n 2 . . . n N Þ.
The system responses (outputs) are defined as functionals R n h ð Þ ð Þ: H ! R and are assumed to belong to the L 2 ðH; PÞ space of second order random variables with finite variance, equipped with the usual inner product Following the Riesz representation theorem, any response can then be expressed in the form of where a i are the expansion coefficients and W i n ð Þ are the orthogonal basis functions. For practical reasons, the infinite sum has to be truncated, and so the approximation contains only the first P þ 1 terms.

Basis set definition
In PCE, the W i n ð Þ basis functions are defined as multidimensional polynomials. The i th basis is built by tensorization of one dimensional polynomials, according to where c i ¼ ðc i;1 ; c i;2 ; . . . ; c i;N Þ is the multi-index in which each component is the polynomial order of the j th random variable n j . w j;c i;j is the polynomial type corresponding to the random variable n j of order c i;j and it can belong to different families (e.g., Hermite, Legendre, Laguerre) depending on the probability density function of the random variable according to the Wiener-Askey scheme (Le Maitre and Knio, 2010;Xiu and Karniadakis, 2002). As customary, we considered basis functions up to a maximum polynomial order O. Hence, their set is defined by Consequently, the number of polynomial terms is which is also equal to the number of a i coefficients to be calculated for each response.

Non-Intrusive Spectral Projections
If the response approximation is obtained by a projection on the basis ðW 0 ; W 1 ; . . . ; W P Þ, the approach is called Non-Intrusive Spectral Projections (NISP). Since the one dimensional polynomials w j;c i;j n j À Á are orthogonal for each input variable, the obtained multidimensional basis functions W i n ð Þ are also orthogonal with respect to joint pdf of the inputs, thus holds, where f 2 i is a (positive) constant and d i;j is the Kronecker delta. The advantage of using orthogonal polynomials is that the coefficients can be computed via the ratio The problem is therefore shifted to the efficient calculation of the integrals in Eq. (8). We performed these calculations through a deterministic approach using Gauss quadrature. The advantage of Gauss formulae is that they ensure exact integral evaluation for high polynomial orders, up to order 2n lev À 1, where n lev is the number of quadrature points associated with quadrature level lev. The values of weights and abscissas depend on the kind of polynomial used for the variable representation. A disadvantage of Gauss quadrature is that they are only merely nested. Since the calculation of expansion coefficients now only requires the evaluation of responses corresponding to prescribed realization of the inputs, this method can be applied without any alteration of the original computational model, thus is ''non-intrusive".

Smolyak sparse grids
In this work, the ''curse of dimensionality" is tackled by combining Gauss formulae with Smolyak sparse grids. Let Q ð1Þ lev be the quadrature rule of level lev. The superscript ''ð1Þ" indicates that the quadrature rule is applied to a single variable n j as Q ð1Þ lev f ¼ where n k j and w k are the abscissas and weights according to the Gauss quadrature rule. We define the following special notation for sparse grids: Thus, D ð1Þ lev f is a quadrature rule itself, having the nodes of both Q ð1Þ lev and Q ð1Þ levÀ1 , with weights equal to the difference between the weights of grid lev and lev À 1. From (9), it follows that We then introduce the level multi-index l ¼ l 1 ; . . . ; l N ð Þto extend the quadrature definition to N variables, signaling the level of quadrature used to integrate along the different dimensions for a given N dimensional cubature. Adopting Smolyak sparse grids and assuming that the maximum integration level lev is the same in each direction, the set of level multi-indices included in the approximation of the integrals of 8 is given by Thus, the Smolyak sparse grid of level lev can be expressed as where we used tensorization to construct the multi-dimensional cubatures approximating the multi-dimensional integrals of 8 from the one dimensional quadratures. In this way, a considerable number of function evaluations f n ð Þ (i.e., model runs) can be eliminated without significant degradation of accuracy.

Post Processing: Sensitivity analysis
Due to orthogonality of the basis vectors, the response mean value can be expressed as Independently from the polynomial type, W 0 ðnÞ ¼ 1 and thus the integral in (13) is the projection of the response on the 0th order polynomial, equal to a 0 . Similarly, the variance V can be computed as It is important to note that, unlike in MC estimators, the mean value and the variance are not explicitly related to the number of evaluations, since they are calculated directly from expansion coefficients.
To study the sensitivity of different responses to the uncertain input parameters, we adopted a variance-based global sensitivity analysis performed with Sobol sensitivity indices. This approach is model-independent and it can quantify inter-dependencies among inputs (Saltelli et al., 2008). Sobol indices split the variance of the response into contributions from each input (and combinations of them). This approach is therefore extremely useful to identify which uncertain parameters are the most critical and which have negligible effects on the considered system responses.
Assuming a set of independent and uniformly distributed inputs Þin the interval n i 2 ½0; 1 8i ¼ 1; . . . ; N, they combine to give the responses R n ð Þ through the general function f, such that R n ð Þ ¼ f ðnÞ. Following Saltelli et al. (2010), the response R n ð Þ can be written as and so on. E n $i indicates the conditional expected value computed considering each input except for n i . E n $ij and the other terms follow analogous definitions. Subsequently, the variance contributions are defined as and so on. Thus, V i is the variance relative to the input n i , after having averaged the contributions of other inputs. The variance of R can be then decomposed as The first order sensitivity index relative to the input n i is finally computed by normalization with respect to the total variance as It is therefore clear that S i provides a measure of the influence on the variance of each input singularly. However, the response variance can have non negligible higher-order contributions related to the combination of variables. We define a total sensitivity index for each input, which contains all the contributions to the variance related to the input n i including the interactions with any other variable as well. It is defined mathematically as The sum over the inputs of the total sensitivity indices is always P 1, since higher order contributions are accounted for more than once. If P d i¼1 S tot;i ¼ 1, the model is purely additive and the contribution of inter-dependencies is zero.
With a PCE meta-model, all Sobol sensitivity indices can easily be evaluated using analytic expressions of the PCE expansion coefficients (Saltelli et al., 2008;Perkó et al., 2014b).

MSFR design and reference model
In this section, we briefly describe the MSFR design and the main characteristics of the multi-physics computational model used in this work. The current design is a 3000MW th fast reactor operating with a liquid fluorides mixture (of lithium, thorium, 13% molar enriched uranium, and other fissile nuclides). The salt, which is both fuel and coolant, freely flows upwards within a toroidal core cavity, without any moderator or control rod. Then, it exits towards sixteen identical sectors where a pump and a heat exchanger are placed. Treatment units are also present to separate the helium bubbles dispersed in the salt, used to remove gaseous and metallic fission products and to control reactivity (Delpech et al., 2009). The core is surrounded by a fertile blanket and by reflectors to improve neutron economy. The main nominal design parameters considered in this study are summarized in Table 1.

Reactor design
The power plant is completed by an intermediate circuit where an inert salt circulates transferring the thermal energy from the fuel circuit to an energy conversion system consisting of a Joule-Brayton cycle.

Multi-physics reactor model
To model the MSFR steady-state behavior, we employed a highfidelity multi-physics simulation tool recently developed at Delft University of Technology. Fig. 2 schematically shows its structure.
A Reynolds-Averaged Navier-Stokes and k À equations solver (DGFlows) is coupled with PHANTOM-S N , a multi-group neutron transport solver. PHANTOM-S N is a discrete ordinates solver based on the discontinuous Galerkin finite element method and is capable of handling steady state and time dependent problems on three-dimensional unstructured meshes. It is also capable of per- forming criticality and time-eigenvalue problems. More details can be found in (Koph'azi and Lathouwers, 2012). The DGFlows CFD code is a Navier-Stokes solver for low-Mach number flows including RANS turbulence models. The code uses a pressure correction method for momentum-pressure coupling. The viscous terms are discretized with the Symmetric Interior Penalty method and convection with Lax-Friedrichs fluxes. Its geometric capabilities are the same as of PHANTOM-S N . Details can be found in (Tiberga et al., 2020a;Hennink et al., 2021). Explicit coupling is realized via exchange of data between the two tools: the fission power density (P dens ) is transferred to DGFlows as it is a source term in the energy equation. The inverse route is followed by the velocity and eddy viscosity fields (u and m t ) that influence the distribution of delayed neutron precursors. Finally, the average temperature in each element (T E ) is used to correct the element-wise macroscopic cross sections (R), taken from a reference library, in order to properly model density and Doppler feedback. The exchange of data is effected by the Galerkin projection between the meshes used in the neutron transport and flow codes. As our meshes are chosen to be hierarchic (hence both the fluid mesh and the neutron transport mesh result from refinement of a common master mesh), such projections are straightforward. Note that both codes support local refinement to be able to focus on relevant physics, hence in some locations the radiation transport mesh may be finer whereas in other regions the reverse may be true. In the present paper, local refinement was not used though.
The steady-state solution is found by iterating the two solvers until convergence. For each iteration, DGFlows seeks the new steady solution through a pseudo-transient, while PHANTOM-S N solves a criticality eigenvalue problem using the desired total reactor power as normalization criterion. As only the steady-state behavior of the MSFR was investigated in this work, decay heat was not modeled. Table 2 reports the composition of the fuel salt mixture along with the properties considered for the nominal, non-stochastic case. Neutronics data were condensed into six-energy groups (whose structure is reported in Table 3, together with the average scalar fluxes in the 6 groups under nominal, non-stochastic conditions) and evaluated at temperature T 0 ¼ 900K with Serpent (Leppanen et al., 2015) starting from the JEFF-3.1.1 nuclear data library (Santamarina et al., 2009). Eight families of delayed neutron precursors were considered.
In this work, all calculations were performed considering only a single recirculation loop as computational domain, thus exploiting the symmetry of the problem. The geometry used is illustrated in Fig. 3. The complete domain was meshed into 52869 tetrahedra for neutronics calculations. As heat transfer within blanket and reflectors was not modeled, these regions were removed from the CFD mesh, which therefore consists of fewer elements (46793). Fig. 4 shows the two meshes employed. A second-order polynomial discretization was used for the velocity field, while all other unknowns were discretized with first-order polynomials.
The pump was modeled as a momentum source, and buoyancy was taken into account adopting the Boussinesq approximation ( Table 2 reports the salt density and thermal expansion coefficient Table 1 MSFR design parameters considered for the nominal, non-stochastic case (Allibert et al., 2016;Gerardin et al., 2017  evaluated at the reference temperature of 973.15 K). A friction factor was imposed in the heat exchanger region to reproduce the resistance to flow, while salt cooling was modeled by a volumetric heat sink term equal to cðT int À TÞ where c and T int are a volumetric heat transfer coefficient and the average temperature of the intermediate salt whose nominal values are reported in Table 1. The interested reader is referred to (Tiberga et al., 2020b) for a more detailed description of the multi-physics code and the MSFR modeling approach.

Sources of uncertainty and considered responses
In this work, the PCE method described in Section 2 was employed to analyze key neutronics and thermal-hydraulics design aspects of the MSFR in steady-state conditions. Regarding the thermal-hydraulics part, emphasis was put on the thermal performance, considering the maximum, minimum and average salt temperatures (T max ; T min ; T avg ) along with the complete temperature field. The reactor's neutronics response was studied in terms of the effective multiplication factor (k eff ). The reference responses, obtained with the nominal, non-stochastic input data, are reported in Table 4 and Fig. 5.
The considered set of stochastic input parameters included several neutronics and thermal-hydraulics data, which are described in Sections 4.1 and 4.2 and are summarized in Table 5.

Neutronics uncertainties
We focused the analysis on the 6-group fission cross sections (R f ;g , where g is the group index) and the 8-family precursor decay constants and delayed neutron fractions (b i ; k i , where i indicates the family). This choice was made since R f ;g are expected to have a large effect on the k eff , being directly related to the production of neutrons, whereas the contribution of b i and k i may be relevant due to the precursors' drift in the complex reactor geometry and turbulent flow, which can affect k eff in an unpredictable way. Since capture and scatter are less affected by the precursors' drift, they were not considered in this preliminary uncertainty analysis study in order to keep the number of stochastic parameters limited.
These parameters are material properties and are generally measured with a sufficiently high number of experiments. Therefore, their statistics can be approximated with a normal distribution (Oberkampf et al., 2000;Roy and Oberkampf, 2011). We assumed the mean value to be equal to the nominal value adopted in non-stochastic simulations, while the relative standard deviation (RSD) was set to 5% (with respect to the mean value) in absence of more precise data. Fig. 4. Mesh adopted for the MSFR model. The neutronics mesh (left) consists of 52869 tetrahedra, while the CFD mesh (right).

Table 4
Values of the reactor responses of interest in nominal conditions for non-stochastic case.

Thermal-hydraulics uncertainties
Generally, the salt thermodynamic properties and the reactor operational conditions have the biggest effect on the salt temperature distribution. For this reason, we chose the salt thermal conductivity (k salt ), the specific heat capacity (c p ), the heat transfer coefficient describing the heat exchanger (c) and the reactor power (P) as uncertain inputs.
Analogously to the neutronics data, we considered the material properties k salt and c p as normally distributed, with mean value equal to the nominal value of non-stochastic simulations. The statistical information of the former were derived from (Benes et al., 2013), where k salt was found to be 1:7 þ À0:4 Wm À1 K À1 . The uncertainty of 0:4 Wm À1 K À1 can be interpreted as a standard deviation (i.e., the RSD is about 23%). The RSD of c p was set to 5% in absence of more precise data. The controllable parameters c and P were approximated with uniform distributions with mean value equal to the nominal value and a half width of variation of 20% with respect to the mean value, and this higher variation was chosen to be more representative of combined data, design and operational condition uncertainties.

Preliminary calculations
Given the 26-dimensional input parameter phase space considered and the use of a high fidelity, complex and computationally expensive multi-physics MSFR model, it was not feasible to perform our study with a single PCE analysis considering all 26 inputs together. Instead, to tackle the ''curse of dimensionality", first preliminary calculations were made to reduce the parameter space without degrading accuracy. Moreover, it was investigated whether the MSFR modeling can be simplified to reduce the com-putational cost of each model evaluation. This section describes these two aspects in detail.

Reducing parameter space by single-physics neutronics simulations
At first, single-physics neutronics simulations were carried out to reduce the problem dimensionality. To evaluate the contribution of each stochastic input, or group of inputs to the reactor k eff , we performed decoupled simulations, perturbing only neutronics data (fission cross sections, delayed neutron fractions or decay constants), while keeping thermal-hydraulics parameters fixed (i.e., considering the salt temperature and flow fields fixed at nominal state). Therefore, in these preliminary calculations we neglected thermal feedback effects and assumed that the overall neutronics parameter ranking is the same regardless of the kind of model evaluation (single-physics or fully-coupled).
To further reduce the number of simulations, we performed separate PCE analyses on small groups of parameters, called classes in the following, assuming no interaction among them. According to the discussion in Section 4, three parameter classes were identified: b (8-members), k (8-members) and R f (6-members). Adopting level-3 quadrature rules and Smolyak sparse grids (lev ¼ 3 in 11) proved sufficient, requiring 161 and 97 model evaluations for the 8-member and the 6-member classes, respectively. Table 6 reports a comparison of the variance of the response k eff relative to each parameter class. The variance due to fission cross sections is almost five orders of magnitude higher than that caused by the other 2 classes. Considering that our study focused on analyzing steady-state conditions, this is not surprising, because the delayed neutron data primarily affects transients. Therefore, the b and k classes were excluded from the set of stochastic input parameters for the coupled calculations.
Secondly, the importance of different energy groups in the R f class was computed in order to potentially eliminate groups with negligible impact on the response variance. The first order and total Sobol sensitivity indices reported in Table 7 show that groups 2, 3, 4, and 5 are responsible for more than 90% of the class variance and have similar sensitivity indices. This is well in line with the reactor spectrum (see the average flux values in each energy group reported in Table 3), since most of the fission reactions take place in these energy groups. Moreover, no significant differences can be seen between the first order and total Sobol indices, indicating minimal higher interactions between these parameters. Consequently, the set of stochastic input parameters (Table 5) was reduced to 4 macroscopic fission cross sections (groups 2 to 5), along with all the thermal-hydraulic parameters previously described (i.e., reactor power, heat exchanger heat transfer coefficient, salt specific heat capacity, and salt thermal conductivity), totaling 8 parameters.

Reducing the computational cost of coupled simulations
The complex MSFR model described in Section 3 is computationally expensive due to the number of degrees of freedom necessary to obtain a sufficiently accurate solution. Moreover, the time step used in DGFlows for the pseudo transient towards the steady-state solution must be limited to avoid numerical instabilities. Solving for the full set of equations in the CFD model at each time step, and for each iteration with the coupled neutronics solver, resulted in an unacceptable computational time needed to obtain a new steady-state solution for each set of input parameters required by the PCE analysis. Therefore, the flow field was fixed assuming it does not change from the nominal steady-state, which significantly reduced the computational cost of each model evaluation. This approximation is reasonable, given that none of the considered stochastic input parameters has a direct impact on flow field (salt viscosity, turbulence parameters, and pump specifications do not vary), and that the contribution of natural circulation to the total flow rate is negligible. Therefore, secondary effects on the flow field induced by the salt temperature variation, due to perturbations in the reactor power, the heat exchanger heat transfer coefficient, and the salt specific heat capacity or thermal conductivity, were considered negligible too.
To validate this computational cost saving approach, we calculated a new steady-state MSFR solution both with and without fixing the flow field. This validation was done for an extreme scenario, considering simultaneous large perturbations to all input parameters at the same time in a manner that could supposedly lead to the highest salt temperature and the minimum effective multiplication factor. These perturbed inputs are summarized in Table 8, and represent a combination of extreme quadrature points along the individual parameter axes for all parameters. For the uniformly distributed P and c the highest and lowest quadrature points were chosen from a level 2 grid, respectively, while for most normally distributed parameters the lowest point from a level 3 grid were taken, in order to yield % 15% variation for all parameters. For the fission cross section of groups 2 and 3 only the second most extreme (level 3) quadrature point was taken. This choice is a compromise between the more extreme scenario (with 5% lower fission cross sections) and computational speed during the initial phase of our study, since using the lower cross sections (i.e., 86% instead of 91%) required excessively long times for convergence with the first, unnecessarily fine computational grids.
Two simulations were performed for the cost saving validation: one in which the flow field was considered fixed at the nominal state (indicated with NoFlow in the following); and the second where the full set of equations was solved for (WithFlow). Fig. 6 shows the evolution of the three thermal-hydraulics responses during the pseudo-transient necessary to reach the new CFD steady-state solution for both cases. These results confirm that considering the flow field fixed at the nominal state has no significant influence on the considered salt temperatures. Moreover, they indicate that the adopted approach is conservative, since the salt maximum temperature is slightly overestimated in the NoFlow case, with a maximum error of 0.11% with respect to the WithFlow simulation. Both simulations the resulting reactor k eff was 0.94025.
A posteriori, to ensure that our approach was conservative, we have also verified that the validation calculation was indeed more extreme than all 161 simulations used to build the PCE model, yielding lower effective multiplication factor (0.94025) and higher maximum temperature ð1140 KÞ than the minimum (0.9815) and maximum ð1123 KÞ across the 161 cases, respectively.

Results of the multi-physics PCE analysis
After reducing the problem dimensionality and the computational cost of each reference model evaluation, we constructed the reactor PCE meta-model by sampling the responses described in Section 4 performing full multi-physics steady-state simulations of the MSFR. The total number of model evaluations required to compute the polynomial expansion coefficients (according to Eq. 8), adopting level-3 quadrature rules and Smolyak sparse grids, for the set of the eight relevant input parameters, was 161. Each multi-physics simulation needed approximately 20 h on average on a high performance computing unit to be completed. The full dataset containing the 161 quadrature points and the corresponding temperature field and neutronic data is openly available .
The PCE approximation was then used to derive uncertainty estimates, in form of pdfs, for each reactor response. In this work, every pdf was obtained with 10 5 samples of the meta-model, which needed only a few seconds in total. To check the accuracy of the PCE approximation, we compared the results obtained with polynomial approximations of order 3 and 4, since the adopted Table 7 Single-physics preliminary neutronics calculations: Sensitivity indices for the R f class reported for each energy group. Only groups 2 to 5 give relevant effects on the effective multiplication factor.

Effective multiplication factor uncertainty
Minimal differences can be noticed between the two polynomial approximations, indicating that the adopted PCE approximation is sufficiently accurate. The k eff mean value is 1.00977, close to the nominal value reported in Table 4, with a variance of 2:070 Â 10 À4 (corresponding to a standard deviation of % 1440 pcm). One can notice that this variance is lower than the one found in the preliminary decoupled neutronics calculations (2:220 Â 10 À4 , see Section 5.1). Density and Doppler thermal feedback have therefore a stabilizing effect on the system from the neutronics point of view, reducing the k eff variance. This is reasonable: wherever fission reactions increase, following a perturbation, the salt temperature increases too, thus introducing negative reactivity which eventually reduces the probability of having fissions.
The normal probability plot (Fig. 7b) shows that the k eff pdf is a nearly perfect normal distribution, indicating that the system has an almost linear behavior from the neutronics point of view, with negligible thermal feedback effects on the response. However, deviations from normality are detected at the tails, probably due to small contributions of the higher-order feedback effects. Table 9 reports the k eff sensitivity indices. The reactor power, the heat transfer coefficient, and the salt thermal properties, which all affect the salt temperature distribution, have indices almost three orders of magnitude (at least) lower than those of the fission cross sections, responsible for most of the total variance. This confirms that thermal feedback effects have very little influence on the Table 8 Perturbed inputs for the computational time reduction test case. They were adopted in both the Withflow and NoFlow cases. Values were chosen to maximize the maximum temperature and minimize the effective multiplication factor, modeling an approximately worst case, highly unlikely scenario. The exact values correspond to the most extreme quadrature points along the single dimensions yielding % 15% perturbations, except for the group 2 and 3 fission cross section, for which the second most extreme point was chosen to speed up convergence.  6. Evolution of the salt maximum (Tmax), minimum (T min ), and average (Tave) temperatures during the pseudo-transient necessary to reach the new CFD steady-state solution after the input parameters perturbation summarized in Table 8. Considering the flow fixed at the nominal state (NoFlow case) has no significant impact. system neutronics response. Finally, as P i S tot;i ¼ 1:000, we deduce that second or higher order interactions between the parameters give negligible contribution to the total k eff variance. Higher order interactions are only relevant for the salt thermal conductivity, for which the total sensitivity coefficient is significantly higher than the first order coefficient due to the interaction with power (S k salt ;P ¼ 2:32 Â 10 À8 ).
6.2. Maximum, minimum, and average temperature uncertainties 6.2.1. Maximum temperature Fig. 8a compares the pdfs of the maximum temperature obtained with 3 rd and 4 th order PCE approximations. Only minor differences can be noticed, in the range T 2 ½1050 K; 1070 K. Analyzing the PCE coefficients, fourth order terms derive from interaction among thermal-hydraulic and neutronics parameters and, thus, they can be associated to thermal feedback effects. Taking fourth order polynomials, the T max mean value is 1083 K (very close to the nominal value reported in Table 4) and has a variance of 478 K 2 . The pdf is far from a normal distribution, as evident from the normality plot reported in Fig. 8b, not only at the tails but also in the central region, especially in the interval T 2 ½1050 K; 1100 K, where the pdf resembles a uniform distribution.
The T max sensitivity indices are reported in Table 10. The salt specific heat capacity has a significant impact, since for a given reactor power and flow rate it influences the temperature difference across the heat exchanger, and therefore also the maximum salt temperature in the system. Unsurprisingly, the reactor power P is the most relevant parameter, contributing to almost 85% of the total variance. In fact, the maximum salt temperature is directly influenced by the power-to-flow ratio. On the other hand, the contribution of c is limited to less than 5%, despite the parameter having an obvious impact on the salt temperature distribution. This is explained by the fact that the maximum salt temperature is located close to the core vertical symmetry axis in proximity of the upper reflector (Fig. 5), whereas the heat exchanger is located in the outer-core region, far from this location.
Fission cross sections have negligible impact on the salt maximum temperature. This is expected, because they influence only mildly the power density distribution, but not the total reactor power that is considered an independent input parameter. The significant difference between the first order and total Sobol indices is due to the interaction with power for all fission cross sections. The contribution of the salt thermal conductivity is also minimal, since the salt flow is highly turbulent, so any perturbation of the molecular thermal conductivity has negligible influence on the total effective one. Finally, it can be seen that P i S tot;i ¼ 1:002 % 1, so overall second and higher-order interactions between parameters are not very relevant.

Minimum temperature
The minimum salt temperature is a key parameter for the safety and reliability of a molten salt reactor, due to the risk of salt solidification. However, in this work the average temperature of the intermediate salt at the secondary side of the heat exchanger was considered fixed at 908:15 K, which is higher than the salt solidification temperature. As a consequence, the fuel salt minimum temperature is limited from below and the presence of stochastic data can only influence the heat exchanger performances. Fig. 9 shows the T min pdfs obtained with 3 rd and 4 th order PCE approximations, again showing almost perfect overlap. The mean value of T min is 924 K (equal to the nominal value reported in Table 4), while the variance is 26 K 2 . Given the lower bound imposed on the minimum salt temperature by the boundary conditions, the pdf is highly skewed. In fact, for T < 920 K the pdf sharply goes to zero (as heat transfer to the secondary side is hindered by the smaller DT). The normality plot reported in Fig. 9b shows that the pdf shape is nearly normal in the central region but not at the tails.
With the given set of inputs, the probability of having a minimum temperature above 923 K, which is one of the design constraints considered in the nominal case to have a sufficient safety margin from the salt solidification temperature, is limited to only 53.2%, which may be concerning. Table 9 Multi-physics PCE results: Sensitivity indices of k eff . Fission cross sections are responsible for most of the response variance, whereas other thermal fluid-dynamics parameters have negligible influence.
The T min sensitivity indices are reported in Table 11. The minimum temperature is found at the exit of the heat exchanger, far from the central core region where the fission power density is the highest. Thus, it is reasonable to see that the c total sensitivity index is by far the highest (around 0.8), whereas the reactor power, despite being still relevant, has only a secondary impact, along with the c p .
Similarly to T max , all other parameters have negligible first order effects, with their total indices of the order of 10 À4 . The sum of the total indices is P i S tot;i ¼ 1:012, which indicates some nonnegligible higher-order effects, mainly coming from the interaction between the power and c (S P;c ¼ 0:01). Fig. 10a compares the average temperature pdfs obtained with 3 rd and 4 th order PCE approximations. Differences are minimal and limited to the interval T 2 ½950 K; 970 K. The T ave mean value is 966 K, while its variance is 73 K 2 . The pdf deviates substantially from a normal distribution at the tails, as illustrated in Fig. 10b. The T ave sensitivity indices are reported in Table 12. As for the previously analyzed temperature responses, fission cross sections and thermal conductivity have negligible effects, and just like before, the difference between total and first order indices are due to interaction with power. The reactor power and the exchanger heat transfer coefficient together contribute to more than 99% of the total variance, with the former being almost twice more important than the latter. This is reasonable, as any perturbation in the heat transfer coefficient has an impact only on the average salt temperature in the heat exchanger, and therefore acts only indirectly on the global average salt temperature. The salt c p contribution is very limited, as it has an impact mainly on the temperature difference across the heat exchanger but not on the average salt temperature. Finally, P i S tot;i ¼ 1:005, which indicates that, as for T min , higher-order effects due to the interaction between parameters are small, but non-negligible.

Complete temperature field uncertainty
An additional study was carried out on the complete salt temperature distribution, outputted by the reference model in form of volume-averaged temperatures per each mesh finite element (T E ). A PCE approximation was constructed for each averaged value, so 46793 separate responses were considered.
Knowing the temperature field is particularly important to ensure salt chemical stability, satisfactory structural material properties, and effectiveness of heat exchange between primary and secondary loop. Since the whole temperature field stochastic, our probabilistic analysis aims to discover possible variations from the nominal design condition that pose potential risks. First, it may be relevant to quantify the probability that the temperature is lower or higher than a certain threshold in each element. This is the case for example for the upper and lower thresholds of 1023 K and 923 K, which are the design targets for the outlet and inlet temperatures in the heat exchanger (Mathieu et al., 2009;Allibert et al., 2016). Second, it may be important to identify the parts of the reactor where the temperature is limited below or   above a threshold with a given confidence level (i.e. with a probability higher than a requested value). For the present study, we considered a probability of 0:95, and again 923 K or 1023 K as threshold temperatures. Fig. 11 shows the probability in each element of having a temperature below 1023K. Despite the design criteria that in the domain T E cannot be higher than the threshold, there are large regions, near the upper reflector in which the probability is far below 95% or even null. This is reasonable because, even in the nominal non-stochastic case, the temperature in that region is higher than 1023 K (Fig. 5). Fig. 11b clearly highlights that in the entire upper region of the core there is a probability lower than 95% that T E < 1023 K . As known from the nominal case, the pump and the upper reflector are easily subjected to temperatures higher Fig. 10. Multi-physics PCE results: (a) Tave pdfs, obtained with third and fourth order polynomials approximations and 10 5 samples. Only minor differences can be noticed in the interval T 2 ½950 K; 970 K; and (b) Normal probability plot relative to the fourth order polynomial PCE approximation. The pdf follows a normal distribution in the central region whereas discrepancies appear at the tails.

Table 12
Multi-physics PCE results: Sensitivity indices Tave. P and c have the greatest impact on the response, followed by cp, while other thermodynamic properties and neutronics parameters first order interactions can be neglected. However, total sensitivity indices show importance of second or higher order interactions.
S i 6:55 Â 10 À1 3:29 Â 10 À1 1:18 Â 10 À2 4:46 Â 10 À9 1:90 Â 10 À8 5:39 Â 10 À8 3:52 Â 10 À8 4:02 Â 10 À9 Stot 6:59 Â 10 À1 3:33 Â 10 À1 1:21 Â 10 À2 4:23 Â 10 À5 4:23 Â 10 À5 4:23 Â 10 À5 4:23 Â 10 À5 4:23 Â 10 À5 Fig. 11. Multi-physics PCE results: (a) Probability in each element of finding a temperature lower than 1023 K (mid-plane cut). In most of the domain it is impossible that TE is higher than this threshold, but there are wide regions, close to the upper reflector, where the probability is very low or even null; and (b) Distinction between regions in which the probability is higher (red) and lower (blue) than 95%. than 1023 K. Analyzing the inlet of the heat exchanger is more interesting. Fig. 12 displays multiple cross sections of the heat exchanger at different heights (measured from the outlet). It can be noticed that at inlet, and still 5 cm downstream, there is a wide region in which the probability that T E < 1023 K is below 95%. Only at 10 cm from the inlet T E is always below 1023 K with the given confidence level. This highlights that at least the first 10 cm are subject to T E P 1023 K with high probability, which was not predicted by the non-stochastic nominal simulations. Therefore, stronger high-temperature protection is likely required for this component. The T min pdf (Fig. 9a) showed that there is a non-negligible probability to have a minimum temperature below 923 K. For this reason, following the same approach as above, we investigated the probability to have a temperature higher than 923 K in each element. Ensuring a sufficient safety margin from the salt solidification temperature is in fact paramount for the safety of reactor operations. Fig. 13 reports the results, which are again influenced by the fact that the salt temperature is limited from below by the boundary condition on the average intermediate salt temperature imposed in the reference model. There is a quite extensive area between the heat exchanger and the core inlet in which elements can have a T E < 923 K with a high probability (up to almost 30%).
Focusing on the heat exchanger, the cross sections reported in Fig. 14 show that the portion between the outlet and 14 cm upstream can be subjected to temperatures lower than the desired threshold.

Conclusions
In this paper, we have successfully applied a previously validated Polynomial Chaos Expansion (PCE) method to perform a preliminary uncertainty and sensitivity analysis of the Molten Salt Fast Reactor (MSFR) behavior at steady-state. We have adopted a high-fidelity multi-physics simulation tool developed at Delft University of Technology as reference model. PCE with Non-Intrusive Spectral Projections and sparse grids has been able to greatly reduce the number of model evaluations to retrieve the statistical information on some selected reactor responses, namely the effective multiplication factor, the maximum, minimum and average salt temperatures, together with the complete salt temperature field.
Starting from a set of 26 stochastic input parameters, both neutronics and thermal-hydraulics related, we have isolated the ones that gave relevant contributions on the responses' variance. Using these preliminary results to reduce the parameters space, we built Fig. 12. Multi-physics PCE results: Multiple cross sections of the heat exchanger, at different heights (measured from the outlet), showing the region in which TE < 1023 K with a probability higher (red) or lower (blue) than 95%. As energy is progressively transferred to the secondary loop, the salt cools down increasing the probability to have a TE below the threshold. However, the first 10 cm are likely subjected to temperatures higher than 1023 K. Fig. 13. Multi-physics PCE results: (a) Probability in each element of finding a temperature higher than 923 K (mid-plane cut). A quite large portion between the heat exchanger outlet and the core inlet might be subjected to lower temperatures, and (b) Distinction between regions in which the probability is higher (red) and lower (blue) than 95%.
an accurate PCE approximation of the full multi-physics model with only 161 evaluations.
Results have shown that the quite large standard deviation of the effective multiplication factor, around 1440 pcm, is mostly affected by the uncertainties of the macroscopic cross sections. In particular, group cross sections in the energy range between 0.7 keV and 0.5 MeV have the largest contribution to the total variance. Therefore, research should focus on reducing as much as possible the measurement uncertainty of these parameters, in order to narrow the multiplication factor pdf. This is further underlined by the fact that our study did not consider capture and scatter cross section uncertainties, hence most probably represents only a lower limit on the k-effective uncertainty.
It has been found that the maximum, minimum, and average salt temperatures are sensitive mainly to the reactor power, the heat transfer coefficient and only partially to the specific heat capacity (especially the maximum salt temperature). Unlike the effective multiplication factor, reducing the variance of these responses is more related to the better definition of the MSFR operational conditions and the components design, such as the heat exchanger. The salt thermal conductivity has proven to have almost no effect on the considered responses. This indicates that further improvements on the statistical information of this parameter might not be necessary.
Finally, the analysis on the complete salt temperature distribution has highlighted that the inlet portion of the heat exchanger can be subjected, with high-probability, to temperatures higher than the 1023 K (750°C) currently considered as design parameter. Therefore, this component should be equipped with stronger hightemperature protections. At the same time, close to the heat exchanger outlet, the fuel salt can likely reach temperatures lower than 923 K (650°C). In the nominal case, this value is the lower limit imposed on the salt temperature to guarantee a sufficient safety margin from salt solidification, whose effects are particularly detrimental for the operations of the heat exchanger (solidified salt might clog channels, drastically increasing pressure drops and compromising heat transfer capabilities). Therefore, to improve the safety and reliability of the MSFR operations, future efforts should focus on increasing the likelihood that a sufficient safety margin from salt solidification is guaranteed.

Disclaimer
The content of this paper does not reflect the official opinion of the European Union. Responsibility for the information and/or views expressed therein lies entirely with the authors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.