Uncertainty Quantification Approach on Numerical Simulation for Supersonic Jets Performance

One of the main issues addressed in any engineering design problem is to predict the performance of the component or system as accurately and realistically as possible, taking into account the variability of operating conditions or the uncertainty on input data (boundary conditions or geometry tolerance). In this paper, the propagation of uncertainty on boundary conditions through a numerical model of supersonic nozzle is investigated. The evaluation of the statistics of the problem response functions is performed following ‘Surrogate-Based Uncertainty Quantification’. The approach involves: (a) the generation of a response surface starting from a DoE in order to approximate the convergent–divergent ‘physical’ model (expensive to simulate), (b) the application of the UQ technique based on the LHS to the meta-model. Probability Density Functions are introduced for the inlet boundary conditions in order to quantify their effects on the output nozzle performance. The physical problem considered is very relevant for the experimental tests on the UQ approach because of its high non-linearity. A small perturbation to the input data can drive the solution to a completely different output condition. The CFD simulations and the Uncertainty Quantification were performed by coupling the open source Dakota platform with the ANSYS Fluent® CFD commercial software: the process is automated through scripting. The procedure adopted in this work demonstrate the applicability of advanced simulation techniques (such as UQ analysis) to industrial technical problems. Moreover, the analysis highlights the practical use of the uncertainty quantification techniques in predicting the performance of a nozzle design affected by off-design conditions with fluid-dynamic complexity due to strong nonlinearity.


Introduction
In recent years, an important part of the research in numerical simulations, in many engineering sectors, has been dedicated to the performance prediction of systems and components in off-design conditions. The behavior of a system or a component in conditions different from those for which it was designed is not only related to the deterministic variation of the input parameters, but also to the aleatory uncertainty which can characterize the input data and the geometrical tolerances. Hence, in order to improve the accuracy and the reliability of the numerical predictions, it is necessary to understand how the uncertainties can affect the results of the problem under investigation. This is one of the main targets of Uncertainty Quantification (UQ) analysis with direct positive fall-out on engineering problems. Computational Fluid Dynamics (CFD) is one of the disciplines where UQ is increasingly applied within a simulation environment. The UQ analysis is becoming an effective approach for industrial use thanks to the concurrent development of both soft-computing methods and computer performance. The use of CFD coupled to optimization algorithms for the automatic design optimization of industrial components is nowadays a mature technology [1,2]. More recently, the development of Response Surface Methodology (RSM) has significantly enhanced the effectiveness of the design optimization based on simulation [3,4], as an industrial standard or for data mining and diagnostic applications [5,6]. The same RSM approach can be used to develop very efficient frameworks for UQ analysis [7]. A contribution to this field is given in the paper with reference to one of the most interesting examples of compressible flow thanks to its high non-linearity: the adiabatic flow in a variable section duct. It is a relevant case for its physical and mathematical background, and it has also a wide range of engineering applications. The best known application is for spacecraft propulsion in aircraft engines or in supersonic wind tunnel [8]. Convergent-divergent nozzles are also applied in gas burners [9,10] and in supersonic water separators for natural gas purification [11,12]. The preliminary design of a supersonic nozzle (De Laval nozzle) is based on the basic gas-dynamic equations in order to have a reference, ideal, operative condition of a fully expanded jet [8]. However, the theory is based on some assumptions, such as isentropic flow conditions, that can affect the actual nozzle performance. Moreover, the real operating conditions (pressure, temperature, gas mixture composition, etc.) can cause a significant change in the flow structure and in the shock waves generation with direct outcomes on nozzle performance.
The main goal of this paper is to answer to the following questions: a small variation of some input parameters can remarkably affect the performance of a designed convergent-divergent nozzle? Can UQ methods help to quantify the above effect and give more insight into the nozzle behavior?
The uncertainties propagation through the CFD model of the nozzle was performed using an automated procedure developed by the authors within the Dakota open source platform. The UQ approach is a 'surrogate-based' approach: a response surface was generated from a DoE, and then, the Latin Hypercube Sampling (LHS) was applied to the meta-model to perform the UQ analysis. The use of a surrogate model is highly effective in reducing calculation times when a sampling-based UQ method is used: in this case, in fact, a high number of response function evaluations is required to generate converging statistics.
The LHS was performed both for the DoE generation and for the UQ analysis. In LHS, the design space (with dimensions equal to the problem variables) is subdivided into an orthogonal grid with N elements per parameter. Within the grid, N sub-volumes are located so that along each row and column of the grid, only one sub-volume is chosen. Inside each sub-volume, a sample is chosen randomly [13]. The LHS method is widely used because it ensures a better coverage of the design space [14][15][16][17] and a faster convergence with respect to the basic Monte Carlo method.
The surrogate model was generated through the Gaussian Process (also known as Kriging) approach which uses a Gaussian correlation function with parameters that are selected by Maximum Likelihood Estimation (MLE); this correlation function results in a response surface that is C ∞ continuous [14].
In a first application, an optimized nozzle geometry with air flowing at a supersonic Mach number was considered. The commercial software ANSYS Fluent ® was used to generate a series of CFD simulations to identify the influence of the nozzle discharge pressure on the convergent-divergent operating condition, flow structure and performance. A perturbation range of the exit static pressure was chosen with respect to the critical gas-dynamic condition, where the normal shock wave is positioned at the nozzle exit section (referred to the r2 condition in classical literature). The choice of this particular condition is related to the potential performance drop that will occur if the shock wave moves into the divergent part of the nozzle. After a preliminary sensitivity analysis, a full UQ approach was developed to quantify the effects of the input variable distribution (ambient static pressure) on the nozzle performance.
A second application considers the flow of natural gas in a convergent-divergent nozzle to understand the effects of the gas chemical composition on the expansion nozzle performance. The gas composition variation can occur on a daily or seasonal base, according to the geographical area and pipeline used [18], and it is a serious problem for several industrial applications, where large quantities Algorithms 2020, 13, 130 3 of 17 of natural gas are treated and energy consumption is an issue [19,20]. The same investigation approach was applied as for the first case and the effectiveness of the UQ analysis is confirmed.

CFD Model
The fluid domain consists of a 2D axisymmetric De Laval nozzle with an "expansion volume" (Figure 1) at the exit, which represents the discharge environment at constant pressure. The axisymmetric RANS equations were solved using an 'axis' boundary condition at the nozzle axis, as common practice [21]. A structured grid was used for the spatial discretization. As shown in Figure 2, a grid refinement was performed near the nozzle wall, in order to ensure a y+ value close to one. The Reynolds Averaged Navier Stokes Equations (RANS) in the steady form are numerically solved with the addition equation for the turbulence closure. The two equation k-ω SST turbulence model is selected based on previous experience on the application case; it allows a good accuracy both in proximity and far from the walls. A second order upwind scheme was selected for the spatial discretization, because if the flow is not aligned with the structured mesh (e.g., oblique shock waves) the first-order convective discretization is not appropriate.

CFD Model
The fluid domain consists of a 2D axisymmetric De Laval nozzle with an "expansion volume" (Figure 1) at the exit, which represents the discharge environment at constant pressure. The axisymmetric RANS equations were solved using an 'axis' boundary condition at the nozzle axis, as common practice [21]. A structured grid was used for the spatial discretization. As shown in Figure  2, a grid refinement was performed near the nozzle wall, in order to ensure a y+ value close to one. The Reynolds Averaged Navier Stokes Equations (RANS) in the steady form are numerically solved with the addition equation for the turbulence closure. The two equation k-ω SST turbulence model is selected based on previous experience on the application case; it allows a good accuracy both in proximity and far from the walls. A second order upwind scheme was selected for the spatial discretization, because if the flow is not aligned with the structured mesh (e.g., oblique shock waves) the first-order convective discretization is not appropriate.  -total pressure and total temperature at the nozzle inlet; -static pressure on expansion volume (outlet domain) side walls.

CFD Model
The fluid domain consists of a 2D axisymmetric De Laval nozzle with an "expansion volume" (Figure 1) at the exit, which represents the discharge environment at constant pressure. The axisymmetric RANS equations were solved using an 'axis' boundary condition at the nozzle axis, as common practice [21]. A structured grid was used for the spatial discretization. As shown in Figure  2, a grid refinement was performed near the nozzle wall, in order to ensure a y+ value close to one. The Reynolds Averaged Navier Stokes Equations (RANS) in the steady form are numerically solved with the addition equation for the turbulence closure. The two equation k-ω SST turbulence model is selected based on previous experience on the application case; it allows a good accuracy both in proximity and far from the walls. A second order upwind scheme was selected for the spatial discretization, because if the flow is not aligned with the structured mesh (e.g., oblique shock waves) the first-order convective discretization is not appropriate.  -total pressure and total temperature at the nozzle inlet; -static pressure on expansion volume (outlet domain) side walls. total pressure and total temperature at the nozzle inlet; -static pressure on expansion volume (outlet domain) side walls. In the first part of the work, the operating fluid is air as an ideal compressible gas with constant Cp. Sutherland's law with three coefficients (Equation (1)) is adopted to model the variation of viscosity with temperature.
In the second part of the article, the working fluid is natural gas treated as an ideal and compressible gas and the viscosity is computed through Sutherland's law with suitable coefficients for methane.

Sensitivity Analysis on Discharge Pressure
The value of the discharge pressure p s has a remarkable impact on the nozzle flow structure, especially in determining the position and strength of the shock wave: this is evident referring to the 'r2' condition where the shock wave is positioned at the nozzle exit section [8]. As a proof of the reliability of the CFD model, similar results have been obtained in similar conditions in the experimental case by Zapryagaev et al. [22]. Figure 3 shows Mach number contours for five different p s values (80 kPa ÷ 120 kPa) around the 'r2' condition. The flow structure and the shock wave location change significantly even with a small ∆ of 10 kPa (about 10% of the pressure range mean value).  In the first part of the work, the operating fluid is air as an ideal compressible gas with constant Cp. Sutherland's law with three coefficients (Equation (1)) is adopted to model the variation of viscosity with temperature.
In the second part of the article, the working fluid is natural gas treated as an ideal and compressible gas and the viscosity is computed through Sutherland's law with suitable coefficients for methane.

Sensitivity Analysis on Discharge Pressure
The value of the discharge pressure ps has a remarkable impact on the nozzle flow structure, especially in determining the position and strength of the shock wave: this is evident referring to the 'r2' condition where the shock wave is positioned at the nozzle exit section [8]. As a proof of the reliability of the CFD model, similar results have been obtained in similar conditions in the experimental case by Zapryagaev et al. [22]. Figure 3 shows Mach number contours for five different ps values (80 kPa÷120 kPa) around the 'r2' condition. The flow structure and the shock wave location change significantly even with a small Δ of 10 kPa (about 10% of the pressure range mean value). The performance variation of the nozzle within the above range is further confirmed by the charts in Figure 4, which represent the Mach and the Pressure Ratio (PR) along the nozzle axis (Equation (2)). The performance variation of the nozzle within the above range is further confirmed by the charts in Figure 4, which represent the Mach and the Pressure Ratio (PR) along the nozzle axis (Equation (2)).

/p (2)
It is evident that the position of the shock wave close to the nozzle exit is affected by the small change in the environment discharge static pressure. Another important parameter, which confirms the performance variation, is the Δpt (Equation (3)) i.e., the total pressure drop through the nozzle that quantifies the mechanical losses due to viscous or compressible flow effects. In Table 2, the large variation of this quantity was reported: losses tend to decrease with outlet static pressure in a nonlinear way (see the difference between cases 1 and 2 compared to 2 and 3); this is related to the shock wave location (inside/outside the convergent-divergent).
This off-design behavior makes the UQ analysis useful to study the component response to a non-deterministic input variable.

Uncertainty Quantification Analysis
A UQ analysis was carried out considering the discharge pressure as the variable parameter. The Quantities of Interest (QoI) considered as output performance parameters are the Mach number at the nozzle exit (M exit) and the total pressure drop (Δpt). With a fully automated procedure implemented in the Dakota environment, a DoE for the two QoI based on LHS method was performed for the pressure range 80-120 [kPa]. The outlet static pressure values provided by Dakota, according to LHS scheme, were used as the input variable in Fluent CFD software, which allowed us to obtain the corresponding response functions values (see Figure 5). In order to approximate the 'physical' (numerical) model, 121 samples were chosen inside the design space. Equations (4) to (6) show the response functions in the normalized form:  It is evident that the position of the shock wave close to the nozzle exit is affected by the small change in the environment discharge static pressure. Another important parameter, which confirms the performance variation, is the ∆p t (Equation (3)) i.e., the total pressure drop through the nozzle that quantifies the mechanical losses due to viscous or compressible flow effects. In Table 2, the large variation of this quantity was reported: losses tend to decrease with outlet static pressure in a non-linear way (see the difference between cases 1 and 2 compared to 2 and 3); this is related to the shock wave location (inside/outside the convergent-divergent). This off-design behavior makes the UQ analysis useful to study the component response to a non-deterministic input variable.

Uncertainty Quantification Analysis
A UQ analysis was carried out considering the discharge pressure as the variable parameter. The Quantities of Interest (QoI) considered as output performance parameters are the Mach number at the nozzle exit (M exit ) and the total pressure drop (∆p t ). With a fully automated procedure implemented in the Dakota environment, a DoE for the two QoI based on LHS method was performed for the pressure range 80-120 [kPa]. The outlet static pressure values provided by Dakota, according to LHS scheme, were used as the input variable in Fluent CFD software, which allowed us to obtain the corresponding response functions values (see Figure 5). In order to approximate the 'physical' (numerical) model, 121 samples were chosen inside the design space. Equations (4)-(6) show the response functions in the normalized form: Algorithms 2020, 13, 130 6 of 17  Kriging or 'Gaussian Process' Response Surface Methodology is considered and the response surfaces for the above functions were generated from the training points given by the DoE. The UQ analysis was then performed through Latin Hypercube sampling method, considering 1000 samples taken on the generated response surfaces.
In order to ensure the meta-model reliability, a cross-validation analysis was performed. In Dakota, the type of cross-validation which can be exploited is the k-fold cross-validation: at first, the DoE dataset is divided into k partitions and k meta-models are generated, each excluding the k-th partition of training data. Each surrogate is tested at the points that were excluded in its generation and the user-specified diagnostic metrics are computed with respect to the held out data [14]. In this work, a particular type of k-fold cross-validation was performed, i.e., the Leave-one-out crossvalidation or Prediction Error Sum of Squares (PRESS). In this special case, the number of partitions is equal to the number of data points. The results obtained from this analysis, reported in Table 3, include the Root Mean Squared, the Mean Absolute Value and the Maximum Absolute Value of the prediction error (calculated between the observed value and the surrogate model prediction for the training data points). Moreover, three random operating points were selected from the UQ database of 1000 samples, and then, simulated with the CFD solver. The corresponding fluid dynamics results for the QoI were given in input to Dakota as a challenge set of data. The metrics calculated with respect to this dataset are reported in Table 4.
The outcomes of the previous analyses were supplemented by the calculation of the relative percentage error (Table 5) between the CFD and the surrogate model response function values. The obtained results certify the reliability of the meta-model.  Kriging or 'Gaussian Process' Response Surface Methodology is considered and the response surfaces for the above functions were generated from the training points given by the DoE. The UQ analysis was then performed through Latin Hypercube sampling method, considering 1000 samples taken on the generated response surfaces.
In order to ensure the meta-model reliability, a cross-validation analysis was performed. In Dakota, the type of cross-validation which can be exploited is the k-fold cross-validation: at first, the DoE dataset is divided into k partitions and k meta-models are generated, each excluding the k-th partition of training data. Each surrogate is tested at the points that were excluded in its generation and the user-specified diagnostic metrics are computed with respect to the held out data [14]. In this work, a particular type of k-fold cross-validation was performed, i.e., the Leave-one-out cross-validation or Prediction Error Sum of Squares (PRESS). In this special case, the number of partitions is equal to the number of data points. The results obtained from this analysis, reported in Table 3, include the Root Mean Squared, the Mean Absolute Value and the Maximum Absolute Value of the prediction error (calculated between the observed value and the surrogate model prediction for the training data points). Moreover, three random operating points were selected from the UQ database of 1000 samples, and then, simulated with the CFD solver. The corresponding fluid dynamics results for the QoI were given in input to Dakota as a challenge set of data. The metrics calculated with respect to this dataset are reported in Table 4. The outcomes of the previous analyses were supplemented by the calculation of the relative percentage error (Table 5) between the CFD and the surrogate model response function values. The obtained results certify the reliability of the meta-model. At first, a uniform distribution of the p s variable were considered as the input probability density function (pdf): in this case, any pressure values within the established range have the same probability. Results of the uncertainty propagation through the surrogate model are the output functions (M exit and ∆p t ) discretized pdfs, reported in Figure 6. It is evident that the pdfs of the QoI are very far from the uniform distribution of the input variable (p s ). The physical explanation of these non-uniform pdfs can be found, considering that the 'r2' state matches with the case p s = 100 [kPa], i.e., the mean of the uncertain input variable distribution. Around this value, two kinds of solutions are highly probable, as evidenced by the exit Mach and ∆p t pdfs ( Figure 6): • high exit Mach and low losses → overexpanded jet with shock waves outside the nozzle; • low exit Mach and high losses → shock waves inside the nozzle.  At first, a uniform distribution of the ps variable were considered as the input probability density function (pdf): in this case, any pressure values within the established range have the same probability. Results of the uncertainty propagation through the surrogate model are the output functions (M exit and Δpt) discretized pdfs, reported in Figure 6. It is evident that the pdfs of the QoI are very far from the uniform distribution of the input variable (ps). The physical explanation of these non-uniform pdfs can be found, considering that the 'r2' state matches with the case ps = 100 [kPa], i.e., the mean of the uncertain input variable distribution. Around this value, two kinds of solutions are highly probable, as evidenced by the exit Mach and Δpt pdfs ( Figure 6): • high exit Mach and low losses → overexpanded jet with shock waves outside the nozzle; • low exit Mach and high losses → shock waves inside the nozzle.
This means that even with a large variation of the input variable, the most frequent flow structures can be restricted to two main ranges of M exit and Δpt. In order to provide more reliable results for the UQ analysis, a second case with an increased number of samples (10,000, one order of magnitude greater than the previous one) was considered. The comparison among the histograms obtained for the QoI is reported in Figure 7: differences between the two distributions can be detected, but they are so small that the 1000 samples UQ analysis can be considered satisfactory. This means that even with a large variation of the input variable, the most frequent flow structures can be restricted to two main ranges of M exit and ∆p t .
In order to provide more reliable results for the UQ analysis, a second case with an increased number of samples (10,000, one order of magnitude greater than the previous one) was considered. The comparison among the histograms obtained for the QoI is reported in Figure 7: differences between the two distributions can be detected, but they are so small that the 1000 samples UQ analysis can be considered satisfactory. Another case was performed, assigning a normal distribution for the input pdf with a mean value 'm' (100 [kPa]) centered on the previously considered variation range and a standard deviation σ chosen in order to have a small perturbation of the pressure: e.g., σ ≈ 1% of the mean value. The resulting pdf for the input variable ps is shown in Figure 8. Using again the LHS method for the uncertainty quantification propagation (1000 samples as before), the distributions reported in Figure 9 are obtained for the response functions. It is evident, especially for the losses, that even in this case, the output pdfs are different from the input one (normal distribution). Table 6 summarizes the data of the statistical moments calculated for these distributions. The red curves over the histograms of Figure 9 represent the reference Gaussian continuous distributions associated with the mean and standard deviation values of each discretized distribution obtained from Dakota. Comparing to the normal pdf (in red), the M exit distribution is more flat (negative value of Kurtosis) and slightly asymmetrical. The distribution of the total pressure drop is shifted to the right (negative Skewness) and sharper (positive Kurtosis) with respect to the Gaussian pdf. Another case was performed, assigning a normal distribution for the input pdf with a mean value 'm' (100 [kPa]) centered on the previously considered variation range and a standard deviation σ chosen in order to have a small perturbation of the pressure: e.g., σ ≈ 1% of the mean value. The resulting pdf for the input variable ps is shown in Figure 8. Another case was performed, assigning a normal distribution for the input pdf with a mean value 'm' (100 [kPa]) centered on the previously considered variation range and a standard deviation σ chosen in order to have a small perturbation of the pressure: e.g., σ ≈ 1% of the mean value. The resulting pdf for the input variable ps is shown in Figure 8. Using again the LHS method for the uncertainty quantification propagation (1000 samples as before), the distributions reported in Figure 9 are obtained for the response functions. It is evident, especially for the losses, that even in this case, the output pdfs are different from the input one (normal distribution). Table 6 summarizes the data of the statistical moments calculated for these distributions. The red curves over the histograms of Figure 9 represent the reference Gaussian continuous distributions associated with the mean and standard deviation values of each discretized distribution obtained from Dakota. Comparing to the normal pdf (in red), the M exit distribution is more flat (negative value of Kurtosis) and slightly asymmetrical. The distribution of the total pressure drop is shifted to the right (negative Skewness) and sharper (positive Kurtosis) with respect to the Gaussian pdf. Using again the LHS method for the uncertainty quantification propagation (1000 samples as before), the distributions reported in Figure 9 are obtained for the response functions. It is evident, especially for the losses, that even in this case, the output pdfs are different from the input one (normal distribution). Table 6 summarizes the data of the statistical moments calculated for these distributions. The red curves over the histograms of Figure 9 represent the reference Gaussian continuous distributions associated with the mean and standard deviation values of each discretized distribution obtained from Dakota. Comparing to the normal pdf (in red), the M exit distribution is more flat (negative value of Kurtosis) and slightly asymmetrical. The distribution of the total pressure drop is shifted to the right (negative Skewness) and sharper (positive Kurtosis) with respect to the Gaussian pdf.
The results analysis suggests that in this case, the values of the QoI are spread over a wider range of values. Hence, a small variation in the input discharge pressure (with respect to r2 discharge pressure which is approximately 100 [kPa]) can change drastically the operative condition, and consequently, the nozzle performance.   The results analysis suggests that in this case, the values of the QoI are spread over a wider range of values. Hence, a small variation in the input discharge pressure (with respect to r2 discharge pressure which is approximately 100 [kPa]) can change drastically the operative condition, and consequently, the nozzle performance.
The comparison between the pdf histograms for the QoI in both cases of 1000 and 10,000 samples is reported in Figure 10. Even in this case, it can be observed that the increase in the number of samples used for the LHS method has no remarkable influence on the pdfs and on the corresponding normal distributions, which remain almost unchanged.   The comparison between the pdf histograms for the QoI in both cases of 1000 and 10,000 samples is reported in Figure 10. Even in this case, it can be observed that the increase in the number of samples used for the LHS method has no remarkable influence on the pdfs and on the corresponding normal distributions, which remain almost unchanged.   The results analysis suggests that in this case, the values of the QoI are spread over a wider range of values. Hence, a small variation in the input discharge pressure (with respect to r2 discharge pressure which is approximately 100 [kPa]) can change drastically the operative condition, and consequently, the nozzle performance.
The comparison between the pdf histograms for the QoI in both cases of 1000 and 10,000 samples is reported in Figure 10. Even in this case, it can be observed that the increase in the number of samples used for the LHS method has no remarkable influence on the pdfs and on the corresponding normal distributions, which remain almost unchanged.

Sensitivity Analysis on Gas Composition
In the previous part of the paper, an important operating parameter like the nozzle discharge pressure was considered as the input variable to quantitatively appreciate its influence on convergent-divergent performance using UQ approach.
The purpose of this section is to quantify the effect of a change in natural gas composition on the nozzle performance with fixed boundary conditions (total pressure, static pressure and total temperature). The natural gas case is of particular interest for industrial applications; the actual chemical composition of the gas delivered to a site is not constant, but can change significantly on a daily or seasonal base according to the pipeline used and the geographical extraction region. Table 7 summarizes some reference data of the most interesting cases for the Italian market: four different real natural gas compositions were selected and their properties were calculated. The variables obtained are the specific heat at constant pressure C p and the molecular weight MW of the mixture. The choice of these variables is related to the option of the fluid specification inside the CFD code that allows an easy automation of the simulation process through scripting. Four CFD simulations were performed, keeping constant the nozzle pressure ratio close to the r2 gas-dynamic condition (shock wave at the nozzle exit): Figure 11 shows Mach number contours for the different working conditions.

Sensitivity Analysis on Gas Composition
In the previous part of the paper, an important operating parameter like the nozzle discharge pressure was considered as the input variable to quantitatively appreciate its influence on convergent-divergent performance using UQ approach.
The purpose of this section is to quantify the effect of a change in natural gas composition on the nozzle performance with fixed boundary conditions (total pressure, static pressure and total temperature). The natural gas case is of particular interest for industrial applications; the actual chemical composition of the gas delivered to a site is not constant, but can change significantly on a daily or seasonal base according to the pipeline used and the geographical extraction region. Table 7 summarizes some reference data of the most interesting cases for the Italian market: four different real natural gas compositions were selected and their properties were calculated. The variables obtained are the specific heat at constant pressure Cp and the molecular weight MW of the mixture. The choice of these variables is related to the option of the fluid specification inside the CFD code that allows an easy automation of the simulation process through scripting. Four CFD simulations were performed, keeping constant the nozzle pressure ratio close to the r2 gas-dynamic condition (shock wave at the nozzle exit): Figure 11 shows Mach number contours for the different working conditions.
Comparing the results to those related to the discharge pressure (ps) sensitivity analysis, it is evident that operating fluid properties (Cp and MW) have a lower influence on the flow structure and on the Mach number values inside the nozzle.  Comparing the results to those related to the discharge pressure (p s ) sensitivity analysis, it is evident that operating fluid properties (C p and MW) have a lower influence on the flow structure and on the Mach number values inside the nozzle. This is also confirmed by the charts in Figure 12, where the Mach number and PR variation along the nozzle axis are displayed. A slight change of the shock wave position can be observed, caused by the different gas properties. This is also confirmed by the charts in Figure 12, where the Mach number and PR variation along the nozzle axis are displayed. A slight change of the shock wave position can be observed, caused by the different gas properties. On the other hand, due to gas composition change, a significant variation of the total pressure losses can be noticed (see Table 8): more than 10 [kPa] difference within the examined dataset.

Uncertainty Quantification Analysis
The UQ analysis for the gas composition was performed with the same approach of the previous test. The operating range for the variables Cp and MW (domain space) was identified from the preliminary sensitivity analysis: the specific heat at constant pressure varies from 2031 [J/(kgK)] to 2214 [J/(kgK)], whereas the MW is in the range 16.1÷19.5 [kg/kmol]. The considered QoI are again the nozzle exit Mach number (Mexit) and the total pressure drop (Δpt-Equation (3)).
The DoE was generated with 121 samples according LHS methodology. The response functions scattered distributions are reported in Figure 13 (MW* and Cp* are in normalized form, as in Equations (7) and (8)).  On the other hand, due to gas composition change, a significant variation of the total pressure losses can be noticed (see Table 8): more than 10 [kPa] difference within the examined dataset.

Uncertainty Quantification Analysis
The UQ analysis for the gas composition was performed with the same approach of the previous test. The operating range for the variables C p and MW (domain space) was identified from the preliminary sensitivity analysis: the specific heat at constant pressure varies from 2031 [J/(kgK)] to 2214 [J/(kgK)], whereas the MW is in the range 16.1 ÷ 19.5 [kg/kmol]. The considered QoI are again the nozzle exit Mach number (M exit ) and the total pressure drop (∆p t -Equation (3)).
The DoE was generated with 121 samples according LHS methodology. The response functions scattered distributions are reported in Figure 13 (MW* and Cp* are in normalized form, as in Equations (7) and (8)).
Algorithms 2019, 12, x FOR PEER REVIEW 12 of 18 This is also confirmed by the charts in Figure 12, where the Mach number and PR variation along the nozzle axis are displayed. A slight change of the shock wave position can be observed, caused by the different gas properties. On the other hand, due to gas composition change, a significant variation of the total pressure losses can be noticed (see Table 8): more than 10 [kPa] difference within the examined dataset.

Uncertainty Quantification Analysis
The UQ analysis for the gas composition was performed with the same approach of the previous test. The operating range for the variables Cp and MW (domain space) was identified from the preliminary sensitivity analysis: the specific heat at constant pressure varies from 2031 [J/(kgK)] to 2214 [J/(kgK)], whereas the MW is in the range 16.1÷19.5 [kg/kmol]. The considered QoI are again the nozzle exit Mach number (Mexit) and the total pressure drop (Δpt-Equation (3)).
The DoE was generated with 121 samples according LHS methodology. The response functions scattered distributions are reported in Figure 13 (MW* and Cp* are in normalized form, as in Equations (7) and (8)).  Cp * = Cp − Cp min Cp max − Cp min (8) A surrogate model based on the DoE (121 samples) was obtained using Kriging Response Surface Methodology. The resulting response surfaces for the QoI are shown in Figure 14.  The cross-validation results reported in Table 9 confirm the model reliability. Even in this second application, the validation with three random points of the database was performed.
The results of the challenge set and the comparison in terms of relative percentage error with respect to the CFD results are reported in Tables 10 and 11 respectively.  The cross-validation results reported in Table 9 confirm the model reliability. Even in this second application, the validation with three random points of the database was performed. The results of the challenge set and the comparison in terms of relative percentage error with respect to the CFD results are reported in Tables 10 and 11 respectively. The low values of the errors certify the reliability of the response surfaces. The UQ analysis was implemented according LHS method (1000 samples) in order to evaluate the component response to input variables statistical distributions. A Gaussian probability distribution function was considered for both variables C p and MW. The characteristics of the pdfs are reported in Table 12 and the corresponding trends are shown in Figure 15.  The low values of the errors certify the reliability of the response surfaces. The UQ analysis was implemented according LHS method (1000 samples) in order to evaluate the component response to input variables statistical distributions. A Gaussian probability distribution function was considered for both variables Cp and MW. The characteristics of the pdfs are reported in Table 12 and the corresponding trends are shown in Figure 15.  The mean values m1 and m2 of the pdfs were chosen as the arithmetic mean of the selected ranges. The standard deviations σ1 and σ2 are about 1% of the corresponding mean values in order to simulate the effects of a small perturbation of the involved variables.
The uncertainties propagation through the meta-model determines the output distributions for the QoI (see Figure 16); the statistical moments of the resulting pdfs are summarized in Table 13. The mean values m 1 and m 2 of the pdfs were chosen as the arithmetic mean of the selected ranges. The standard deviations σ 1 and σ 2 are about 1% of the corresponding mean values in order to simulate the effects of a small perturbation of the involved variables.
The uncertainties propagation through the meta-model determines the output distributions for the QoI (see Figure 16); the statistical moments of the resulting pdfs are summarized in Table 13.  The 'shape' of the response functions pdfs is very different from the corresponding reference Gaussian distribution: both of them are sharper (positive Kurtosis) and slightly shifted from the mean. The range of M exit variation, with the selected input, is very narrow. The Mach number is slightly influenced by the gas composition, confirming what was previously observed in the sensitivity analysis.
Regarding total pressure losses, it can be noticed that even with a small variation in Cp and MW, the range of the possible total pressure drop values (Δpt) is quite wide, and consequently, the probability to fall into a non-optimal operating condition is high.
Finally, to test the results with a higher number of samples, a second UQ analysis was performed with 10,000 LHS samples. The comparison between the statistical distributions of the QoI for the two cases are reported in Figure 17. The very small differences in the pdf trends confirm that 1000 samples are sufficient to accurately detect the system response to input data.  The 'shape' of the response functions pdfs is very different from the corresponding reference Gaussian distribution: both of them are sharper (positive Kurtosis) and slightly shifted from the mean. The range of M exit variation, with the selected input, is very narrow. The Mach number is slightly influenced by the gas composition, confirming what was previously observed in the sensitivity analysis.
Regarding total pressure losses, it can be noticed that even with a small variation in Cp and MW, the range of the possible total pressure drop values (∆p t ) is quite wide, and consequently, the probability to fall into a non-optimal operating condition is high.
Finally, to test the results with a higher number of samples, a second UQ analysis was performed with 10,000 LHS samples. The comparison between the statistical distributions of the QoI for the two cases are reported in Figure 17. The very small differences in the pdf trends confirm that 1000 samples are sufficient to accurately detect the system response to input data.  The 'shape' of the response functions pdfs is very different from the corresponding reference Gaussian distribution: both of them are sharper (positive Kurtosis) and slightly shifted from the mean. The range of M exit variation, with the selected input, is very narrow. The Mach number is slightly influenced by the gas composition, confirming what was previously observed in the sensitivity analysis.
Regarding total pressure losses, it can be noticed that even with a small variation in Cp and MW, the range of the possible total pressure drop values (Δpt) is quite wide, and consequently, the probability to fall into a non-optimal operating condition is high.
Finally, to test the results with a higher number of samples, a second UQ analysis was performed with 10,000 LHS samples. The comparison between the statistical distributions of the QoI for the two cases are reported in Figure 17. The very small differences in the pdf trends confirm that 1000 samples are sufficient to accurately detect the system response to input data.