UNCERTAINTY QUANTIFICATION IN STEADY STATE SIMULATIONS OF A MOLTEN SALT SYSTEM USING POLYNOMIAL CHAOS EXPANSION ANALYSIS

Uncertainty Quantification (UQ) of numerical simulations is highly relevant in the study and design of complex systems. Among the various approaches available, Polynomial Chaos Expansion (PCE) analysis has recently attracted great interest. It belongs to nonintrusive spectral projection methods and consists of constructing system responses as polynomial functions of the stochastic inputs. The limited number of required model evaluations and the possibility to apply it to codes without any modification make this technique extremely attractive. In this work, we propose the use of PCE to perform UQ of complex, multi-physics models for liquid fueled reactors, addressing key design aspects of neutronics and thermal fluid dynamics. Our PCE approach uses Smolyak sparse grids designed to estimate the PCE coefficients. To test its potential, the PCE method was applied to a 2D problem representative of the Molten Salt Fast Reactor physics. An in-house multi-physics tool constitutes the reference model. The studied responses are the maximum temperature and the effective multiplication factor. Results, validated by comparison with the reference model on 10 Monte-Carlo sampled points, prove the effectiveness of our PCE approach in assessing uncertainties of complex coupled models.


INTRODUCTION
The Molten Salt Fast Reactor (MSFR) is one of the six new concept reactors proposed by the Generation IV International Forum. The liquid nuclear fuel represents the main design feature and guarantees many intrinsic safety advantages [1]. On the other hand, this reactor requires strong efforts in R&D, both from the experimental and the numerical point of view, due to the early stage of development. In this perspective, Uncertainty Quantification (UQ) is a powerful tool to assess the sensitivity of key design aspects to the variation of intrinsically uncertain input parameters.
UQ was developed to infer the statistical information on the response(s) of a system assuming a variability in the set of data, models, and tools. It has gained increasing importance in the last decades thanks to a change in the design procedure of nuclear reactors from a conservative approach to the so-called Best Estimate Plus Uncertainty methodology, which focuses on the evaluation of the variability of the results in the most precise and complete way [2]. UQ on aleatory variables has been traditionally carried out using collocation methods, such as Monte Carlo (MC) or Latin Hypercube Sampling, which generally need a large number of model evaluations to infer statistical information [3]. Accurate simulations of complex systems are commonly computationally expensive, making these methodologies unfeasible in many situations [4]. For this reason, novel techniques have been developed.
Among these, Polynomial Chaos Expansion (PCE) is chosen in the present work. It belongs to nonintrusive spectral projection methods and consists of constructing system responses as polynomial functions of the stochastic inputs. The potential to strongly reduce the number of model evaluations without needing any modification of underlying codes used to sample the system responses makes this technique extremely attractive [3].
We propose the use of PCE to perform UQ of complex, multi-physics models for liquid fueled reactors, addressing key design aspects of neutronics and thermal fluid dynamics in case of high problem dimensionality. In a recent work [5], a Reduced Order Modeling approach was proposed for the same purpose. To the best of our knowledge, this is the first application of PCE to this kind of problems. Knowing the inputs' statistical information, the PCE meta-model is created by choosing a set of polynomial basis vectors, and using Smolyak sparse grid based numerical integration to compute the basis coefficients [6]. This is described in more detail in Section 2. To test its potential, the PCE method was applied to a simplified 2D problem representative of the MSFR physics, as explained in Section 3. An in-house multi-physics tool [7], coupling an S N radiation transport code with an incompressible Navier-Stokes solver, was considered as the reference model of the molten salt reactor system. The resulting PCE meta model was used to provide uncertainty estimates (in the form of pdfs) of the system maximum temperature and the effective multiplication factor, and to assess their sensitivity to the variation of input parameters. The outcomes of this analysis are detailed and discussed in Section 4.

POLYNOMIAL CHAOS EXPANSION METHOD
The Polynomial Chaos Expansion (PCE) method consists of approximating a response with a polynomial of appropriate order, which contains the response's statistical information as function of the stochastic inputs. Then, for applications such as sensitivity analysis, the sampling is performed on the polynomial approximation constituting the meta-model, requiring just a few seconds for millions of evaluations.
The PCE mathematical treatment is here introduced in a schematic way. Define the probability space Γ = (Θ, Σ, P), where Θ is the sample space such that θ ∈ Θ, being θ the random event, Σ is the σ-algebra, and P is the value of probability. Random inputs ξ(θ) belongs to Υ, the support of their probability density function (pdf) p(ξ). The present work focuses only on independent random variables. Thus, any joint pdf is pξ where N is the number of stochastic inputs. The response is defined such that whereξ = (ξ 1 , ξ 2 ...ξ N ), and L 2 (Θ, P) is the space of second order random variables in which the inner product is defined and it is finite. The response can be expressed in the form where a i represents the i th coefficient and Ψ i the i th basis. For practical reasons, the sum has to be truncated to P + 1 terms.
PC involves the definition of multi-dimensional bases. The i th basis function is built by tensorization of one-dimensional polynomials, according to is the multi-index of the i th basis, whose components are the polynomial orders of the j th random variable ξ j . ψ j,γ i,j is the polynomial type referred to the random variable ξ j of order γ i,j . It can belong to different families depending on the pdf of ξ j . In this work, we focused only on uniform and normal distributions (see Section 4), so we chose Legendre polynomials for the former and Hermite polynomials for the latter. See [3] for a more complete definition of the basis set.

Non-intrusive spectral projections
If the response approximation is obtained by projection on the basis [Ψ 0 , Ψ 1 , ..., Ψ N ], the approach is called Non-Intrusive Spectral Projections (NISP). Then, if the basis is also othogonal, we have where λ 2 i is a (positive) constant, and δ i,j is the Kronecker delta. The advantage of using polynomials with these features reflects on the coefficient calculation. In fact, they can be computed by the ratio The problem is therefore shifted to the efficient calculation of the integral. As the name "Non-Intrusive" suggests, this property allows applying this method without any alteration of the original model. Thus, the PC meta-model can be directly superposed to the already implemented codes for the response evaluation. In this work, the integral is evaluated adopting Gauss quadrature formulae, which allow for exact integral computation for high polynomial order, up to order 2 n lev − 1, where n lev is the number of quadrature points corresponding to level lev. Weights and abscissas depend on the kind of polynomial used for the variable representation (see [3,6] for further information).
Carefully taking the evaluation points, the expansion coefficients in Equation (5) can be computed very efficiently, requiring a number of simulations often order of magnitudes lower than MC methods. Conversely, in case of large number of inputs, as in our case, the number of coefficients drastically increases, requesting a large amount of model evaluations. This issue is known as "curse of dimensionality". In this work, the problem is tackled by combining the Gauss formulae with Smolyak sparse grids (see [6] for more details).

TEST PROBLEM: A SIMPLIFIED MOLTEN SALT REACTOR SYSTEM
In this work, we applied the PCE method to a simplified system representative of the main characteristics of the MSFR: fast spectrum, strong negative temperature feedback, precursors movement, and thus strong coupling between neutronics and thermal-hydraulics. The problem was developed as a benchmark for multi-physics tools dedicated to liquid-fuel fast reactors [8,9]. Figure 1a shows the problem domain, which is a 2 m x 2 m cavity filled with molten salt at initial temperature of 900 K. Its composition is reported in Table 1. The reactor is surrounded by vacuum. All walls are insulated, and a heat sink equal to γ(T ext − T ), where T ext = 900 K and γ is a volumetric heat transfer coefficient (HTC), reproduces the salt cooling. A zero-velocity boundary condition is imposed at all walls, except the top lid, which moves at v lid . The steady-state solution is sought with criticality eigenvalue calculations normalizing the reactor power to P . Fluid properties are constant with temperature and uniform in space. Neutronics data are condensed into 6 energy groups, while delayed neutron precursors are grouped into 8 families. Doppler feedback is ignored, so cross sections are corrected with temperature only via the fuel density feedback. The flow is laminar and buoyancy effects are modeled via the Boussinesq approximation. We refer to [8,9] for a complete description of the problem.
An in-house multi-physics tool based on the Discontinuous Galerkin Finite Element space discretization is used to model the system. Figure 1b displays its structure. It couples a solver for the incompressible Navier-Stokes equations (DGFlows) with a neutronics code solving the multi-group S N Boltzmann equation coupled with the transport equations for the delayed neutron precursors (PHANTOM-S N ). Data are exchanged due to the coupling between the physics characterizing the molten salt reactor: the average temperature on each element (T avg ), based on which cross sections are corrected with the density feedback from a library at 900 K; the velocity field (u), which is an input for the precursors equation; and the fission power density (P f iss ), which is a source for the energy equation. The codes are iterated until convergence to find the steady state solution. More details can be found in [7]. Simulations were performed choosing a 50 by 50 uniform structured mesh, a second-order polynomial discretization for the velocity, and a first-order one for all other quantities. An S 2 discretization was chosen for the angular variable.

RESULTS
PCE was employed to analyze the system response in terms of maximum temperature (T max ) and effective multiplication factor (k ef f ). The initial set of stochastic input parameters were: lid velocity, thermal expansion coefficient, salt viscosity, HTC (γ), reactor power (P ), 6-groups fission cross sections (Σ f,g ), 6-groups scattering cross sections, 6-groups average number of neutrons produced per fission event (ν g ), 6-groups fission prompt spectrum, and 8-families precursors decay constants and fractions. This set includes both material properties and controllable parameters, which have different statistics. The former generally come from a large number of measurements, so they could be represented with normal distribution with mean value equal to the nominal one. In this work, the same variance was associated to each input to understand which has the most relevant impact on the response. In absence of precise data, it was set to 5 % of the mean value. Controllable parameters were statistically represented with uniform distribution, since they can assume any kind of value in an interval of modulation. The mean value coincides with their nominal one and the "standard deviation" (from now on, to be intended as half of the variation interval for uniform pdfs) was taken as 20 % of the mean value in absence of precise data.
A way to select the more important parameters for the responses was needed first, to reduce the problem dimensionality, thus decreasing the number of model evaluations. Therefore, we performed preliminary calculations on decoupled problems (i.e., neutronics with fixed nominal flow and temperature fields, and fluid dynamics with fixed nominal fission power density). Results, reported in [10] only for brevity, showed that only six inputs contribute to the responses variance in a relevant way (i.e., at least to 70 % of the total response variance). We assumed this ranking to be valid also for the fully coupled problem, given the secondary contribution of feedback effects on the responses. The parameters considered are reported in Table 2, together with their statistical properties and mean values, which correspond to the nominal ones. The nominal responses obtained for a coupled calculation with these parameters are T max = 1333 K and k ef f = 0.99525.

Analysis of Responses: Maximum Temperature
The T max pdf is reported in Figure 2. To prove the correctness of the PCE approximation, we first compared the polynomial approximations of order 3, with 10 5 samples, and of order 4, with 10 6 samples. From Figure 2a, it is clear that no relevant discrepancies appear in the two pdfs, which are almost superposed. The main difference between the two cases lies in the number  Figure 2a with those obtained with PCE. The agreement is excellent. The PCE mean value (1337 K) and variance (3793 K 2 ) obtained with the 3 rd order polynomial approximation differ from the MC estimates by 0.17 % and 1.00 % respectively. Moreover, using the 3 rd -order PCE model, we evaluated the responses for the same MC samples, finding a maximum relative error with respect to the MC reference of 0.2 %. Figure 2b reports the normality test performed on the response set obtained with 10 5 samples of the validated 3 rd order PCE model. Despite the curve is almost Gaussian in the range [1270, 1470] K, it shows substantial deviation at the tails. The reason of the pdf asymmetry is related to a small non-linearity of the response, in particular with respect to the HTC.
Finally, we employed the PCE model to perform global sensitivity analysis. The total Sobol indices, S tot (see [4] for their rigorous definition), reported in Table 3 show that γ and P are the only important parameters, being responsible for the 45 % and 55 % of the T max variance, respectively, whereas neutronics-related parameters have a negligible effect. This is reasonable, because γ and P directly control the amount of energy present in the system. Finally, we notice that second or higher order interactions between parameters are negligible, as i S tot,i ≈ 1.

Analysis of Responses: Effective Multiplication Factor
A similar analysis was performed for the other system response, k ef f . Figure 3a compares the pdfs obtained with PCE approximations of order 3 (10 5 samples) and order 4 (10 6 samples). As they are almost perfectly superposed, order 3 polynomial was again chosen as best compromise. The pdf obtained with the 10 3 MC samples of the reference model is also reported in Figure 3a. The agreement is again excellent. The 3 rd -order PCE mean value is 0.99253 and the variance 5.28 × 10 −4 , with a relative error with respect to the MC estimates of 0.04 % and 0.86 % respectively. Moreover, the maximum error of the PCE responses for the same 10 3 samples is 0.03 %. Figure 3b shows a nearly exact normal distribution, due to the irrelevance of thermal fluid dynamics parameters and a linearity of the system, with respect to the neutronics parameters. In fact, the power and the HTC play a secondary role on k ef f , since they affect only indirectly the neutron population, as highlighted by the total sensitivity indices reported in Table 4. The HTC and the power indices are two orders of magnitude lower than the ones of the neutronics inputs. In addition, also in this case, second or higher order interactions between parameters are negligible.

CONCLUSIONS
This paper has presented an application of a PCE method to complex, multi-physics models for liquid fueled reactors, addressing key design aspects of neutronics and thermal fluid dynamics. The PCE method is based on a Non-Intrusive Spectral Projection approach and Smolyak sparse grids designed to estimate the PCE coefficients. It was chosen to reduce as much as possible the number of model evaluations, while retrieving accurate statistical information on responses like the system maximum temperature and effective multiplication factor.       To test its potential, the PCE technique was applied to a simplified 2D problem representative of the MSFR physics, using an in-house multi-physics code as reference model. Having reduced the set of inputs to only six relevant parameters with preliminary, decoupled calculations, the PCE-model required less than 100 simulations to be constructed. The resulting k ef f pdf shows an almost perfect normal distribution, while the one for T max deviates from a Gaussian behavior at the tails. The maximum temperature is mainly sensitive to the reactor power and heat removal coefficient, while the k ef f is affected mostly by neutronics parameters. PCE results were validated by comparison with the reference model sampled on 10 3 Monte Carlo generated points. The maximum relative error was 0.2 % on the responses, lower than 0.2 % on the pdfs means, and 1.0 % on the pdfs variances. It is worth highlighting that the reference model needed around 1500 h to generate the 10 3 responses, while sampling the PCE meta-model required a few seconds for 10 6 evaluations. We conclude that this PCE technique can be effectively used to perform UQ on multi-physics models as those describing the MSFR, thus helping progressing the development of this new reactor concept.