Sensitivity Study of Ice Accretion Simulation to Roughness Thermal Correction Model

Abstract

The effects of atmospheric icing can be anticipated by Computational Fluid Dynamics (CFD). Past studies show that the convective heat transfer influences the ice accretion and is itself a function of surface roughness. Uncertainty quantification (UQ) could help quantify the impact of surface roughness parameters on the reliability of ice accretion prediction. This paper aims to quantify ice accretion uncertainties and identify the key surface roughness correction parameters contributing the most to the uncertainties in a Reynolds-Averaged Navier-Stokes (RANS) formulation. Ice accretion simulations over a rough flat plate using two thermal correction models are used to construct a RANS database. Non-Intrusive Polynomial Chaos Expansion (NIPCE) metamodels are developed to predict the convective heat transfer and icing characteristics of the RANS database. The metamodels allow for the computation of the 95% confidence intervals of the output probability distribution (PDF) and of the sensitivity indexes of the roughness parameters according to their level of influence on the outputs. For one of the thermal correction models, the most influential parameter is the roughness height, whereas for the second model it is the surface correction coefficient. In addition, the uncertainty on the freestream temperature has a minor impact on the ice accretion sensitivity compared to the uncertainty on the roughness parameters.

1. Introduction

Icing represents a major threat to in-flight safety. Ice accretion leads to increased aircraft weight, aerodynamic efficiency loss and a potential increase of up to more than 60% in drag [1]. In the 1997–2006 period, 13% of weather-related accidents were caused by icing [2]. The aerodynamic efficiency loss is determined by the geometry of the ice shape, mainly the thickness and extension [3]. Developments in numerical simulations, especially RANS-based simulations (Reynolds Averaged Navier-Stokes), have enabled the development of icing codes, useful for the effective anticipation of the icing process [4]. Even if icing codes share similar mathematical models and numerical methods, the predicted ice shapes can vary greatly for the same atmospheric conditions [5]. The discrepancies between shapes can be explained by various reasons, both mathematical and numerical, among them the discrepancy between the convective heat transfer models. The fraction of the shape discrepancy that could be attributed to the thermal modeling remains to be quantified.

Most of the icing codes used to predict the ice accretion shape are based on the Messinger model [6]. The model uses a mass and energy balance involving the convective heat transfer to compute the ice growth on the airplane surface [7]. Past studies have shown that the convective heat transfer, modeled by a convective heat transfer coefficient, has a major impact on the final shape of the ice accretion [8,9]. The convective heat transfer coefficient does play a key role in the geometry of the ice shape, which needs to be quantified.

The increase in surface roughness due to early ice formation has an impact on the convective heat transfer in the near-wall region. Experimental studies were conducted by NASA to measure the convective heat transfer coefficient, or, equivalently, the Stanton number, on a roughened NACA0012 airfoil [10,11]. More recently, by Liu & Hu [12], Hosni et al. [13] and Dukhan et al. [14,15] also carried out related experiments, but on a roughened flat plate, highlighting the influence of the roughness characteristics on the Stanton number. An increase in the heat transfer coefficient was observed for a rough surface as compared to a smooth one, in both airfoil and flat plate cases. From the foregoing, it is therefore reasonable to expect that the initial surface roughness of the wing is influencing the complete icing process by acting on the convective heat transfer at the surface [8,16]. Experimental studies show the variability of roughness characteristics, which depend mainly on the airfoil geometry [17] and atmospheric conditions [14]. The incorporation of surface roughness in the RANS model is then subject to uncertainty due to the prevailing lack of knowledge on roughness characteristics in real in-flight situations. The RANS model’s roughness parameters are thus epistemically uncertain variables.

To take into account the surface roughness, Aupoix [18] worked on improvements of the thermal boundary layer calculations by considering a local increment of the Prandtl number. The thermal correction model from Aupoix, denoted as the HAX model, uses three parameters to incorporate surface roughness in simulations: the roughness height, the equivalent sand grain roughness and the surface correction coefficient. Another thermal correction model, called the 2PP model, was developed and is also used in the present study [19]. This model relies on two parameters: the roughness height and the equivalent sand grain roughness. The parameters involved in both the HAX and 2PP models are denoted as the roughness parameters for the rest of the paper. To further improve these two models, sensibility studies can play a key role in establishing the most sensitive parameters. Once identified, it will be possible to vary these parameters in experiments or high-fidelity simulations in order to further calibrate and improve the thermal models.

For uncertain input parameters, the uncertainty quantification (UQ) enables the evaluation of the probability distribution function (PDF) of the output of interest and a complete sensitivity analysis [20]. A relatively small sample of RANS simulations is needed, and a metamodel is generated based on these sample results. This metamodel allows us to quickly compute a large number of new heat transfer coefficients and ice accretion results. The Non-Intrusive Polynomial Chaos Expansion (NIPCE) is commonly used for metamodel generation [21]. The NIPCE is an efficient approach to estimate measurement uncertainties when complex systems are involved, thus explaining its widespread use in UQ [22,23], especially in aeronautics applications [24]. The NIPCE approach is also used in turbulence modeling for UQ, as carried out in [25] to assess the Spalart-Allmaras uncertainty due to the epistemically uncertain model coefficients.

The sensitivity of an output parameter to the input parameter(s) can be evaluated by using a variance-based decomposition [26,27]. In the variance-based family, the Sobol indexes allow for an estimation of the main (first-order indexes) and interaction (second- and higher-order indexes) effects of an input parameter on the model output [28]. The Sobol indexes classify the input parameters according to their level of influence on a specific output parameter. UQ coupled to NIPCE and Sobol indexes has been used for simulations in aeronautics [29,30] and heat transfer [31]. Despite UQ’s wide use in aeronautics and heat transfer applications, no published paper has yet thoroughly characterized, from a statistical point of view, the sensitivity of ice accretion to surface roughness parameters. The closest related work is from Raj et al. [32], where the sensitivity of the ice accretions to the equivalent sand grain roughness is compared to the sensitivity to ice density, droplet distribution and number of ice layers. However, no comparison is made for the roughness parameters between them.

The present work is a first attempt to use UQ in aircraft icing to relate directly heat transfer model parameters to ice shape predictions. In addition, UQ allows the identification of the key parameters. The technical challenge in this case is double. First, an efficient way of generating a large sample set had to be used: making a sensitivity study requires a large database. The metamodel approach is well suited to this challenge, avoiding the computational effort that computational fluid dynamics (CFD) simulations would have required. Secondly, a reliable gauge of the importance of each roughness parameter was needed. The Sobol indexes give us the opportunity of directly reading on a bar chart the classification of the parameters depending on their level of influence on ice accretion. The objective is to quantify ice accretion uncertainty and identify the key roughness parameters contributing the most to the uncertainty of the thermal correction models for RANS equations. The importance of roughness parameters uncertainty compared to the freestream temperature uncertainty on the ice accretion shape is also investigated. To reduce the number of parameters, avoid curvature effects on ice growth and pressure gradients, the airfoil is replaced by a flat plate.

Metamodels generated from a sample of RANS simulations, alternatively using the HAX [18] and the 2PP [19] thermal correction models, are constructed to predict the convective heat transfer and ice accretions on a flat plate. For the ice accretion part, the flat plate has the droplet collection efficiency of a NACA0012, to be representative of an airfoil. The NIPCE metamodels are then used for UQ and to evaluate the Sobol sensitivity indexes associated with the roughness height, the equivalent sand grain roughness, the surface correction coefficient in the case of the HAX model and, finally, the freestream temperature. The novelty of the current work resides in the application of the NIPCE metamodeling method to rough icing simulations. The NIPCE has low sample requirements compared to the traditional Monte-Carlo methods used in the past [33]. In similar applications, Raj et al. [32] use the radial basis function (RBF) metamodeling. The performance of polynomial-type metamodels, like the NIPCE, compared to RBF methods, is better for small sample designs, giving higher regression coefficients [34]. This enhanced accuracy, which combined with the fast response due to the low sample sizes led to the choice of this tool to perform the metamodeling.

The paper starts with a presentation of the physical set-up where the flat plate geometry used is presented and contextualized, along with the HAX and 2PP RANS thermal correction models and the ice accretion model. Then, the methods retained for UQ and sensitivity analysis, including the NICPE metamodeling and the Sobol indexes, are described. Finally, the last section focuses on the results obtained from the simulations and the metamodeling process by highlighting the PDFs of the outputs and the Sobol indexes associated with the roughness parameters. The results allow for the determination of the most influential roughness parameter for both thermal correction models, and show the relative importance of roughness uncertainty compared to the freestream temperature uncertainty on the ice accretion shape.

2. Physical Set-Up

The goal of this section is, first, to define the test case geometrical configuration in the context of the RANS equations coupled to a Messinger model for ice accretion prediction. In this study, RANS equations are solved using a thermal correction model—either the HAX model or the 2PP model, which is presented in the second part of the section.

2.1. Convective Heat Transfers over a Rough Flat Plate

For an iced wing, the quasi-random and unknown characteristics of surface roughness introduce uncertainties in the RANS model. The initial uncertainty on roughness characteristics propagates throughout the thermal correction and ice accretion models, leading to uncertain ice accretion characteristics. Similarly to the previous work of Dukhan et al. [14], a rough flat plate geometry is used to model an airfoil, neglecting the pressure gradient effects on the flow. The flat plate represents the upper part of the airfoil, that was unrolled and retains the curvilinear impinging water droplet collection efficiency of the airfoil, as shown on Figure 1a. These assumptions allow for the computing of an ice accretion similar to what is seen on the airfoil, but neglect the curvature effects. Neglecting the curvature effects will produce a slightly underestimated ice thicknesses compared to a curved airfoil. In addition, the absence of positive pressure gradient will increase the boundary layer thickness. This flat plate simplification should not impact the study conclusions, since the sensitivity is ultimately monitored, and not the absolute ice thickness. The epistemic uncertainty is illustrated by two different roughness characteristics (see Figure 1b). Given their geometrical differences, these patterns lead to different convective heat transfer coefficients (see Figure 1c), which in turn cause differences in the characteristics of ice accretions (see Figure 1d).

The database needed to obtain the metamodels relies on RANS simulations. The geometry used in the present work is a 2-m-long horizontal flat plate. This geometry is commonly used in CFD validation [35], and the mesh size is selected according to a grid study presented in [19]. The 2D mesh of the domain consists of a structured domain of 137 by 97 cells with 114 points on the flat plate surface. The growth rate in the direction perpendicular to the plate is 1.3, and the y⁺ value is between 0.45 and 1.35. In a rough configuration, with this mesh, the difference with the experimental results of Hosni et al. [13] for the Stanton number evaluation was less than 7% for both the HAX and 2PP thermal correction models. In addition, the flat plate configuration was extensively studied both numerically [36] and experimentally [14] in the context of aircraft ice accretion studies. The complete 2.33 m by 1 m domain is sketched in Figure 2. The leading edge of the flat plate is located 0.33 m downstream of the domain inlet.

2.2. Air Mathematical Models

The compressible flow field coming through the inlet boundary of Figure 2 is at Mach 0.2. The Reynolds number based on the length of the plate is 9 million. The Boeing correction for roughness [37] and the thermal correction models are implemented in the Spalart-Allmaras turbulent model of the SU2 6.2.0 solver [19,35]. To obtain the convective heat transfer, the surface temperature is set to 273.15 K, while the ambient air is at 262.04 K. In accordance with the experimental set-up of Dukhan et al. [14], the first 5.2 cm of the plate are unheated, adiabatic and smooth. This starting length allows the boundary layer at the beginning of the rough zone to have a thickness of 1.83 mm.

Two RANS thermal correction models are alternatively used: the HAX model [18] and the 2PP model [19]. Mainly, these models quantify the shift ΔT⁺ of the nondimensional wall temperature in the thermal boundary layer induced by surface roughness. This shift was first observed for the velocity profile, further denoted as Δu⁺ [38].

The temperature shift ΔT⁺ is modelled by increasing the turbulent Prandtl number Pr_t by an additional ΔPr_t in the near-wall region. For the HAX model [18], the Prandtl number correction ΔPr_t is a function of the velocity shift Δu⁺, the normal coordinate to the wall y, the roughness height k, the equivalent sand grain roughness k_s and a non-dimensional corrected surface ratio S_corr [39]. This computation only applies in the region located fewer than five times the roughness height away from the wall (i.e., y/k < 5).

$$\Delta {Pr}_t=\left(A{\left(\Delta u^{+}\right)}^2+B\Delta u^{+}\right)\exp {\left(-y/k\right)}_{ }$$

(1)

with:

$$A=\left(0.0155-0.0035S_{corr}\right)\left(1-e^{-12\left(S_{corr}-1\right)}\right)$$

(2)

$$B= -0.08+0.25e^{-10\left(S_{corr}-1\right)}$$

(3)

$$\Delta u^{+}= \frac{1}{\kappa ·\log \left(1+\frac{R{e}_s}{\exp \left(1.3325\right)}\right)}$$

(4)

The 2PP model [19] depends only on two parameters, k and k_s:

$$\Delta P{r}_t= g\times 0.07083\times R{e_s}^{0.45}\times P{r}^{0.8}\times \exp \left(-\frac{y}{k}\right)$$

(5)

In Equations (4) and (5), the roughness Reynolds number Re_s is defined as u_τk_s/v. Finally, g = 1 if Re_s ≥ 70; $g=\frac{\ln \left(R{e}_s\right)-\ln \left(5\right)}{\ln \left(70\right)-\ln \left(5\right)}$ if 5 < Re_s < 70 and g = 0 if Re_s ≤ 5.

2.3. Ice Accretion Model

Once the heat transfer along the flat plate is known, the ice accretion simulations can be performed. To compute the 2D ice accretions, the Messinger approach is used [6,7]. The basis of this approach is to compute the ice growth rate in a control volume using a mass and an energy balance. The control volume and the contributions to the balances are illustrated in Figure 3.

The mass balance (see Equation (6)) and the energy balance (see Equation (7)) are given by the following set of two equations:

$${\dot{m}}_{in}+{\dot{m}}_{imp}-{\dot{m}}_{ice}-{\dot{m}}_{es}-{\dot{m}}_{out}=0$$

(6)

$${\dot{Q}}_{conv}+{\dot{Q}}_{es}+{\dot{Q}}_{imp}-{\dot{Q}}_{kin}-{\dot{Q}}_{ice}-{\dot{Q}}_{in}=0$$

(7)

The subscripts appearing in Equations (6) and (7) are defined in Figure 3. More details on the terms involved and the resolution procedure are given in the literature [8,40]. The convective heat transfer, whose coefficient is obtained by the models involved in the previous two subsections, has the general expression given by Equation (8):

$${\dot{Q}}_{conv}= h_c\left(T_w-T_{rec}\right)$$

(8)

with the recovery temperature given by:

$$T_{rec}= T\left(1+P{r}^{1/3}·\frac{\gamma -1}{2}M2\right)$$

(9)

The collection efficiency over a NACA0012 airfoil was applied on the flat plate surface and is shown in Figure 4, where s represents the curvilinear abscissa along the airfoil surface, and c represents the chord. The air freestream velocity is 64.9 m/s, the liquid water content is 1.3 g/m³, the droplet diameter is 20 μm and the angle of attack is 0°. Only the upper part of the airfoil, studied through the flat plate representation, is plotted. Note that the plot of Figure 4 is translated for the simulations to have s/c = 0 at the very beginning of the rough zone (at x = 5.2 cm), to have a nonzero boundary layer thickness representative of an airfoil stagnation point.

The icing model presented was successfully verified in previous studies [41], and the ice accretion dependency on the roughness parameters was preliminarily highlighted [8]. The ice thickness is linked to the ice accretion mass, through the ice density of 917 kg/m³, and the exposure time was fixed to 480 s for this work.

3. Uncertainty Quantification and Sensitivity Analysis Method

UQ is used to evaluate the ice accretion uncertainties caused by the heat transfer model parameters in the presence of surface roughness and to identify the most influential parameters. In this section, the uncertain parameter sampling and the NIPCE metamodeling are first depicted, followed by the Sobol indexes definition. In this work, the expression “output of interest”, denoted as Y_i, refers to either one of the following features: (1) the local value of the convective heat transfer coefficient at a specified point or (2) the ice accretion characteristics: maximum thickness or ice extension.

3.1. Uncertain Parameter Sampling and Metamodel

The UQ analysis is used to determine how the roughness’ epistemic uncertainty can impact the PDFs of the ice accretion characteristics. The uncertain inputs of the thermal correction models’ parameters are the roughness height k, its ratio with the equivalent sand grain roughness k_s/k and the corrected surface ratio S_corr for the HAX model [18] and k, k_s/k for the 2PP model [19]. The choice of the ratio k_s/k instead of k_s alone comes from the fact that k and k_s are linked, for example, with the Dirling relation [25].

The roughness height k was experimentally studied in the literature. For example, Dukhan et al. [14] give a range of variations between 0.41 mm and 4.32 mm on a flat plate measuring 0.458 m by 0.458 m, at a Reynolds number of 2,000,000. These values correspond to the roughness elements observed in different atmospheric conditions, and which led to different surface roughness types. The parameter k_s is taken from Fortin [17], who reviewed the validation test cases used for the LEWICE software [42] and concluded on a range of variations for k_s between 0.309 mm and 1.247 mm for airfoils with a chord varying between 53.3 cm and 91.4 cm. Both of these ranges for k and k_s are maintained in the present study and lead to a k_s/k ratio ranging between 0.288 and 0.753. The last input parameter, S_corr, is taken from Aupoix and has a value between 1 and 2.5 [18]. For all the parameters, a uniform PDF between the margins given above is assumed. Although other values are possible in practice,

Table 1. Summary of the distributions for the input parameters.

Parameter Xi	Minimum	Maximum	Distribution
k (mm)	0.41	4.32	Uniform
Ratio ks/k	0.288	0.753	Uniform
Scorr	1	2.5	Uniform

gives the distributions used in this study.

With the range of distribution of each parameter defined, the input parameters can be sampled by Latin Hypercube Sampling (LHS), which is one of the methods commonly used to obtain random, but homogeneously spread, sample points [43]. The HAX model has 190 sample simulations and the 2PP model 120, for their respective analyses, as shown in Figure 5. These sample size ranges are in agreement with the literature, which gives a sample size between 160 and 200 for a three-input model [25]. For instance, Solaï et al. used a sample of 115 simulations for their analysis with four input parameters [44].

The sampling presented gives triplets (HAX model) or doublets (2PP model) of roughness parameters that are used as inputs for the RANS simulations. The database of convective heat transfer coefficients and ice accretion characteristics allows for the construction of the metamodels.

The framework used for the UQ analysis, as well as for the subsequent sensitivity analysis, is the UQLab suite [20]. To be able to generate a large number of results (i.e., up to several thousand for the sensitivity analysis) without the need to run additional RANS simulations, three metamodels constructed from the database are used, as shown in

Table 2. Metamodels generated.

Metamodel	Input Parameters Xj		Output of Interest Yi
	HAX	2PP
M1	k, ks/k, Scorr, x	k, ks/k, x	Convective heat transfer coefficient value at position x
M2	k, ks/k, Scorr	k, ks/k	Maximum ice accretion thickness
M3	k, ks/k, Scorr	k, ks/k	Ice accretion extension

. A metamodel can be represented as a function M estimating the output of interest Y_i from the input vector X, whose components are X_j.

The method retained to create the metamodels is the NIPCE [21]. This existing metamodeling tool is frequently used in UQ, and the present work specifically applies it to the UQ for roughness study in aircraft icing. This method is suitable for models with two and three input parameters. The general expression of a NIPCE metamodel is given by Equation (10):

$$Y_i=M\left(X\right)= {\displaystyle \sum }_{\alpha }{y}_{\alpha }\times {\psi }_{\alpha }\left(X\right)$$

(10)

When generating a metamodel, the leave-one-out (LOO) error value allows an estimation of the metamodel quality. The LOO value, based on the mean square error between the observed and predicted values, quantifies the ability of the metamodel to predict an output value close to what the RANS computation would have calculated [21]. The computation of the LOO error relies on the cross-validation technique. For an N-sample design, the main aspect of the LOO calculation is to create N metamodels Mⁱ, which are constructed on the CFD database where the ith simulation was dropped, leaving N-1 sample points. The difference between the predicted value of the dropped sample by Mⁱ and its real CFD value allows the computation of the LOO error [45] as in Equation (11).

$$LOO= \frac{{\sum }_{i=1}^N\left(Y^{\left(i\right)}-M^i\left(X^{\left(i\right)}\right)\right)2}{{\sum }_{i=1}^N\left(Y^{\left(i\right)}-\overline{Y}\right)2}$$

(11)

In Equation (11), X⁽ⁱ⁾ and Y⁽ⁱ⁾ represent the input parameters of the ith CFD sample simulation and the corresponding CFD output, respectively. Finally, $\overline{Y}$ is the mean value of the output of interest on the sample dataset. The LOO error is a common tool frequently used in metamodeling application. For instance, Jung et al. extensively used the LOO error for anti-icing performance metamodeling for engine intake [46].

3.2. Sobol Indexes Definition

The purpose of a sensitivity analysis is to identify the uncertain parameters that contribute the most to the variations observed on Y_i [26]. Based on the results of the metamodels, the sensitivity of each Y_i to the X_i is quantified. The Sobol sensitivity indexes are based on the decomposition of the total variance of Y_i. The first-order Sobol index for X_i and the second-order Sobol index for X_i and X_j are defined by Equations (12) and (13), respectively:

$$S_i= \frac{V\left(Y\right)-E\left(V\left(Y|X_i\right)\right)}{V\left(Y\right)}$$

(12)

$$S_{i,j}= \frac{V\left(Y\right)-E\left(V\left(Y|X_i, X_j \right)\right)}{V\left(Y\right)}$$

(13)

Finally, for a generic study with three input parameters (i, j, k), the total Sobol index for the parameter i is:

$$S_{Ti}=1-\left(S_j+S_k+S_{j,k}\right)$$

(14)

According to Chan et al. [28], a parameter with a total Sobol index greater than 0.8 can be considered as “very important”, with the total Sobol index between 0.5 and 0.8 as “important”, between 0.3 and 0.5 as “unimportant” and “negligible” with an index below 0.3. Consequently, in the present icing application, a total Sobol index close to 1 characterizes an input parameter having a big effectiveness and influence on the ice accretion.

To study the impact of the temperature on the ice accretion sensitivity, the Sobol indexes are computed in two distinct sets of tests. First, the Sobol indexes are estimated on metamodels M₁, M₂ and M₃ as described in Table 2, with the roughness parameters as the only uncertain inputs. For the second set of tests, the temperature is added as an uncertain input parameter. Along with the roughness uncertain parameters of ice accretion metamodels M₂ and M₃, the freestream temperature T is added as an epistemically uncertain parameter varying uniformly between 261.5 K and 262.5 K. This last series of tests with an uncertain temperature is used to determine if ice accretion is more sensitive to the roughness uncertainty than to the temperature uncertainty.

4. Results and Discussion

With the RANS and the ice accretion models defined, the UQ and Sobol indexes can be applied to study the ice accretion uncertainties over a flat plate. For the outputs of interest, the goals are to evaluate the PDFs and the level of influence of the roughness parameters. For each metamodel and each thermal correction model, the following aspects retain attention: (i) the quality of the metamodel, quantified by the LOO error; (ii) the PDF of the outputs of interest, with the 95% confidence intervals; and (iii) the sensitivity of the outputs of interest to the input parameters, described by the Sobol indexes. For concision, one relevant figure is shown for each aspect listed, and the remaining data for all possible metamodel/thermal correction model combinations are summarized in tables.

The outputs of interest listed in Table 2 are illustrated on Figure 6. Figure 6a shows the outputs of interest on a plot of the convective heat transfer coefficient versus the x position on the flat plate surface. Figure 6b shows the outputs of interest on a graphical representation of an ice accretion. Note that the zone of interest is only the rough part of the flat plate, beginning at x = 5.2 cm.

All the sample simulations will give the same global trends, as observed on Figure 6, with variations induced by the roughness from one sample to another. The convective heat transfer coefficient (Figure 6a) will present a maximum value at the stagnation point before decreasing when traveling downstream along the flat plate. For the ice accretion shape (Figure 6b), the maximum thickness will be at the stagnation point, corresponding to the location of maximum collection efficiency (see Figure 4). The thickness will decrease along the plate with the decrease of the collection efficiency. In addition, the runback water will end up freezing, giving the observed step, located, in the example of Figure 6b, at x = 0.08 m. The step is due to the fact that no more liquid water runs back further downstream.

4.1. Precision of the Metamodels Generated

Prior to looking at the LOO error value, the precision of the metamodel may be visualized by plotting the NIPCE prediction for the set of input parameters used for the RANS simulations versus the actual RANS calculated value. Figure 7 shows this scatter plot for the metamodel M₁ of the HAX model. Y_RANS represents the RANS-calculated values of the local convective heat transfer coefficient at x = 11.36 cm (with x = 0 taken at the very beginning of the flat plate), while Y_NIPCE represents the corresponding values predicted by the metamodel. The points for the 190 sample simulations are plotted and show that the metamodel M₁ for the HAX model has a good precision, with an absolute mean error between the RANS-calculated values and the NIPCE-predicted values of 0.90 W/m²/K and a linear regression coefficient of 0.9999. This regression coefficient value is higher than the ones commonly accepted in the literature for sensitivity analysis [32]. The LOO error for this metamodel, M₁, at x = 11.36 cm is 4.25 × 10⁻⁴. The local position x = 11.36 cm was chosen among the 114 possible surface grid points to illustrate the behavior of the metamodel.

The trend shown in Figure 7, with all the scatter points homogeneously packed close to the identity line, suggests a metamodel making accurate predictions. Checking the precision of all the metamodels leads to the LOO errors presented in

Table 3. Leave-one-out (LOO) errors for all the metamodels.

Metamodel	Output of Interest Yi	LOO Error
		HAX	2PP
M1	Heat transfer coefficient value at 11.36 cm	4.25 × 10−4	2.15 × 10−4
M2	Maximum ice accretion thickness	4.12 × 10−3	1.05 × 10−4
M3	Ice accretion extension	1.52 × 10−2	4.49 × 10−4

.

The values presented in Table 3 show that the convective heat transfer coefficient is well predicted by the metamodels: the LOO error for both HAX and 2PP is in the neighborhood of 10⁻⁴, which is an acceptable value that leads to a precision similar to the one of Figure 7. Nevertheless, the 2PP model has a more predictable behavior than its HAX counterpart. The values show a LOO error two to four times smaller than the one of the HAX model for a given metamodel. For the ice accretion characteristics, the HAX model shows a loss of predictability since the LOO errors rise to the 10⁻³ to 10⁻² range of values. On the other hand, the 2PP model still presents a well-predictable icing behavior with the LOO error once again in the 10⁻⁴ magnitude.

4.2. Outputs of Interest PDFs

The metamodel-generated databases allow the computation of the outputs of interest PDFs. Figure 8 shows all the ice accretion shapes obtained for (a) the HAX model and (b) the 2PP model. Additionally, for each of the models, one ice accretion chosen among the set is plotted alone for clear visualization. Each of the colored lines on the top charts of Figure 8 represents one specific accretion obtained with one particular sample run. For the single plotted accretions, the specific roughness parameters for the HAX model (a) are k = 1.14 mm, k_s/k = 0.599 and S_corr = 1.17, and for the 2PP model (b), k = 2.07 mm and k_s/k = 0.385. The general trend illustrated in Figure 6, with a step observed in the accretion, is visible for all the sample accretions. Figure 8 shows that varying the roughness parameters changes the position of the step, which characterizes the ice extension along the surface.

Using all 190 (HAX model) or 120 (2PP model) data sets allows for the computing of the NIPCE-predicted PDF of the outputs of interest. The NIPCE allows for an increase in sample size, but the PDF is similar to the PDF predicted, with a sample size of 190 or 120. To verify this for the 2PP model, a sample of 10,000 doublets (k, k_s/k) was generated with the Latin Hypercube Sampling, and the NIPCE metamodel applied to the sample. The PDFs of the 120 and 10,000 NIPCE predictions are shown together in Figure 9, where the bars’ heights are scaled by the total number of predictions for each case (i.e., 120 and 10,000, respectively).

Figure 9 shows a similar trend. A closer analysis allows for the determination of the PDF of the data in the figure. The fitted normal distributions, associated to the histograms of Figure 9, have their parameters μ and σ summarized in

Table 4. PDF of normal distribution parameters of Figure 9’s histograms.

120 Samples NIPCE Prediction	10,000 Samples NIPCE Prediction
μ = 0.0522 m [0.0512; 0.0531] σ = 0.0052 [0.0046; 0.0059]	μ = 0.0521 m [0.0520; 0.0522] σ = 0.0047 [0.0046; 0.0047]

, with the 95% confidence intervals in square brackets. The first column in Table 4 gives the parameters for the NIPCE prediction histogram of Figure 9 obtained with 120 samples, and the second column gives the PDF parameters as extracted from the prediction with 10,000 samples. Note that the NIPCE metamodel obtained, defined by a polynomial as in Equation (10), in both the sample sizes (120 and 10,000), is of degree 7.

The parameters in Table 4 show that the NIPCE prediction M₃ with 120 samples matches closely the prediction obtained with 10,000 samples. The PDF parameters of the NIPCE response only change by a magnitude of 10⁻⁴ for μ and 0.0005 for σ. This aspect shows that the sample size of 120 was already enough to estimate the PDF of the output with a good precision. Nevertheless, the confidence intervals are reduced with the increase of the sample size: the confidence interval on μ for 120 samples is 1.9 mm wide, while it is 0.3 mm wide for the 10,000 samples. By applying the relation from [25], with an oversampling ratio of 3 for two uncertain parameters and a NIPCE of degree 7, the minimum sample size required for the problem is 108. This confirms, in addition to the graphical visualization, that the sample size of 120 is enough. This characteristic of having a high enough precision with the RANS sample size is also verified for all the other metamodels. To ensure a small confidence interval, the results obtained for a 10,000 sample size are kept. The normal distributions obtained for each metamodel with the NIPCE prediction with 10,000 samples are summarized in

Table 5. Output PDF for each NIPCE prediction on the 10,000 sample size.

Metamodel	Output of Interest Yi	Output PDFs’ Parameters
		HAX	2PP
M1	Heat transfer coefficient at x = 11.36 cm	μ = 544.2 W/m2 K [542.1; 546.3] σ = 106.0 W/m2 K [104.8; 107.7]	μ = 430.8 W/m2 K [430.2;431.3] σ = 27.8 W/m2 K [27.4; 28.2]
M2	Maximum ice accretion thickness	μ = 0.0264 m [0.0263; 0.0265] σ = 0.0046 m [0.0046; 0.0047]	μ = 0.0218 m [0.0218; 0.0218] σ = 0.0017 m [0.0017; 0.0017]
M3	Ice accretion extension	μ = 0.0512 m [0.0511; 0.0513] σ = 0.0061 m [0.0061; 0.0062]	μ = 0.0521 m [0.0520; 0.0522] σ = 0.0047 m [0.0046; 0.0047]

.

The data in Table 5 show that each thermal correction model predicts slightly different output PDFs for the same metamodel. The μ values in Table 5 show that, statistically, changing from the HAX to the 2PP thermal correction model leads to a 21% variation in the heat transfer coefficient. This change in the heat transfer coefficient then leads to a 17% variation in the ice maximum thickness and only 1.8% of change in the ice extension. The standard deviation σ allows us to estimate the relative uncertainty and how the observed values spread around the mean μ. For a normal distribution, 99.7% of the data are in the range [μ − 3σ; μ + 3σ]. The σ values in Table 5 suggest a larger uncertainty on the outputs for the HAX model compared to the 2PP model. For the ice extension, the relative uncertainty in the 99.7% range [μ − 3σ; μ + 3σ] is about 36% for the HAX model and 27% for the 2PP model. The same calculation for the maximum ice thickness leads to a relative uncertainty of 52% for the HAX model and 23% for the 2PP model. As expected, these uncertainties on the accretion find their origin in the heat transfer uncertainties, which are 58% and 19% for the HAX and 2PP models, respectively. As an example, Figure 10 shows the difference in the distributions predicted by metamodel M₂ for the maximum ice thickness. The difference between the distributions predicted for the HAX model (orange) and the 2PP model (blue) is noticeable. The larger standard deviation for the HAX model is clearly visible, with a larger spread of the values around the mean compared to the 2PP model.

Figure 10 shows that the two thermal correction models predict different maximum ice thickness ranges: the 2PP model gives ice accretions with thicknesses of up to about 23 mm, while the HAX model predicts thicknesses of up to 33 mm.

4.3. Sensitivity Analysis: Calculation of the Sobol Sensitivity Indexes

4.3.1. Sensitivity to the Roughness Parameters

With the metamodels generated, the Sobol sensitivity indexes are evaluated for each input parameter. The 2PP model, with two input parameters, allows the computation of the Sobol indexes up to order 1, while the HAX model, with three input parameters, allows Sobol indexes of order 2. Figure 11 shows the Sobol indexes for the metamodel M₂, evaluating the maximum thickness of the ice accretion and the HAX thermal correction model. The values of the total Sobol indexes obtained for each metamodel are shown in

Table 6. Total Sobol sensitivity indexes for each metamodel.

Metamodel	Output of Interest Y	Total Sobol Sensitivity Indexes
		HAX	2PP
M1	Heat transfer coefficient at x = 11.36 cm	k:  0.239 ks/k: 0.051 Scorr: 0.787	k:  0.820 ks/k: 0.271
M1	Heat transfer coefficient at x = 14.29 cm	k:  0.239 ks/k: 0.050 Scorr: 0.787	k:  0.802 ks/k: 0.291
M1	Heat transfer coefficient at x = 27.48 cm	k:  0.248 ks/k: 0.052 Scorr: 0.779	k:  0.756 ks/k: 0.341
M2	Maximum ice accretion thickness	k:  0.326 ks/k: 0.072 Scorr: 0.682	k:  0.894 ks/k: 0.202
M3	Ice extension	k:  0.166 ks/k: 0.088 Scorr: 0.822	k:  0.912 ks/k: 0.288

. In Table 6, three positions have been chosen for the metamodel M₁, to highlight the differences of sensitivity when traveling on the flat plate.

The results in Figure 11 and Table 6 show, for the HAX model, that the parameter S_corr has the largest impact on every output of interest: S_corr contributes up to 68% to the variability of the maximum ice thickness and 82% for the ice extension. For the 2PP model (see Table 6), the roughness height k has the greatest influence on the outputs of interest as compared to the ratio k_s/k, going up to about 91%. The results observed for the 2PP model, where k has a greater influence than the ratio k_s/k, are also observed for the HAX model. According to the classification given by Chan et al. [28], S_corr is an important or very important parameter for every metamodel, while k and the ratio k_s/k are negligible or unimportant for the HAX thermal correction model. For the 2PP thermal correction model, k is very important compared to the ratio k_s/k, which is negligible. These results show that in the case of the HAX model, S_corr needs to be carefully quantified prior to launching a simulation, since the outputs of interest are the most sensitive to its value. The same conclusion can be drawn for the parameter k when using the 2PP model. Looking at the values in Table 6, it can be observed that the roughness height k sees its level of influence decreasing as it gets away from the leading edge with the 2PP thermal correction. This observation is done by comparing the Sobol indexes from metamodel M₁ at x = 11.36 cm and further away at x = 14.29 cm and x = 27.48 cm, where the Sobol index for k decreases. Meanwhile, the Sobol indexes for k_s/k increase while traveling away from the leading edge. For the HAX model, the Sobol indexes for all three uncertain parameters remain quite steady when traveling downstream.

4.3.2. Sensitivity of the Ice Accretion to the Roughness Parameters and the Freestream Temperature

To highlight the importance of roughness parameters on the ice accretion shape compared to the freestream temperature, the freestream temperature T is now added as an epistemically uncertain variable, varying randomly between 261.5 K and 262.5 K, and the Sobol indexes are computed. For these tests, the parameters are summarized in

Table 7. Metamodels used for the temperature impact study.

Metamodel	Input Parameters Xj		Output of Interest Yi
	HAX	2PP
M2T	k, ks/k, Scorr, T	k, ks/k, T	Maximum ice accretion thickness
M3T	k, ks/k, Scorr, T	k, ks/k, T	Ice accretion extension

.

The new metamodels, with an increased number of uncertain input parameters, have a subscript “T” to avoid confusion. For the sake of concision, the LOO errors and outputs’ PDFs of metamodels M_2T and M_3T are not displayed, but the precision is of the same order of magnitude as that of metamodels M₂ and M₃. The Sobol indexes of metamodel M_2T for the HAX thermal correction model are displayed in Figure 12.

Figure 12 shows that the temperature has a minor impact on the sensitivity of the ice accretion’s maximum thickness. When comparing with Figure 11, it is also possible to note that the level of influence of the roughness parameters remains almost the same despite the addition of the temperature as uncertain. This negligible influence of the temperature is also observed for the other combinations of metamodel/thermal correction model, as summarized in

Table 8. Total Sobol sensitivity indexes for the temperature impact study.

Metamodel	Output of Interest Y	Total Sobol Sensitivity Indexes
		HAX	2PP
M2T	Maximum ice accretion thickness	k:  0.302 ks/k: 0.067 Scorr: 0.667 T: 0.039	k:  0.750 ks/k: 0.170 T: 0.162
M3T	Ice extension	k:  0.190 ks/k: 0.050 Scorr: 0.784 T: 0.022	k:  0.681 ks/k: 0.255 T: 0.176

.

The results in Table 8 show that for every metamodel and thermal correction model, the temperature is the least influent parameter. For the HAX thermal correction model, T is completely negligible, while for the 2PP model it has a greater impact, closer to k_s/k, but still far behind the level of influence of the roughness height k and still below 0.3. These results highlight the importance of the roughness model in the simulation, whose parameters have to be carefully selected, since their variation leads to bigger changes to the ice accretion shape than the temperature variation within 1 K. This relative importance of roughness on the accretion sensitivity compared to other factors related to ambient conditions has also been noticed in the literature. In [32], Raj et al. showed that the roughness plays a bigger role in the accretion shape than the ice density, the droplet distribution and the evaporation.

The results show that the methodology retained is well suited to the ice accretion, since the NIPCE satisfactorily predicts the behavior of the output of interest. The NIPCE metamodeling gave fast responses and facilitated the analysis for the aircraft icing study. The results obtained using the NIPCE tool allow for a direct observation of the key features of the accretion’s characteristics when roughness parameters are varying. Following that, the sensitivity analysis allowed us to determine that the outputs of interest are very sensitive to S_corr when using the HAX model, and to the roughness height k when using the 2PP model.

The ranges of the uncertain roughness parameters come from the literature, and uniform PDFs have been chosen. Changing the ranges and/or the PDFs might impact the results. Other types of metamodels can also be investigated. A similar study can also be conducted around an airfoil to highlight the effects of curvature and pressure gradient effects on the results. Furthermore, k, unlike S_corr and k_s/k, is linked to the skin friction coefficient, which is seldom available and known for an iced surface.

5. Conclusions

This paper has introduced an efficient methodology to quantify the sensitivity of the ice accretion characteristics to the thermal correction parameters on a rough flat plate. The objective of the paper to quantify the ice shape uncertainty and identify the key roughness parameters leading to this uncertainty was reached. The main parts of the methodology were the metamodeling process, using the Non-Intrusive Polynomial Chaos Expansion (NIPCE), and the sensitivity study computing the Sobol sensitivity indexes. The NIPCE metamodeling tool gave entire satisfaction regarding accuracy and quickly computed results with up to several thousand inputs. These features highlighted NIPCE’s suitability for the study, despite its novel usage for aircraft icing with roughness analysis. The study allowed a comparison of two thermal correction models: the HAX model, with three input parameters, and the 2PP model, with two input parameters. The predicted ice accretion’s maximum thickness is very sensitive to the choice of the thermal correction model and its respective parameters. Switching from the HAX to the 2PP thermal correction induces a 17% change in the maximum ice thickness. This aspect is mainly due to the impact of the thermal correction on the convective heat transfer, which is changed up to about 20%. In addition, the standard deviation study highlighted a larger relative uncertainty on the accretion for the HAX model, which can reach 52% on the ice maximum thickness compared to the uncertainty of 23% of the 2PP model. It was also concluded from the Sobol sensitivity indexes that the most influential parameter for the HAX model is S_corr, which is responsible of 68% and 82% of variations on the ice accretion’s maximum thickness and extension, respectively. For the 2PP model, the roughness height k is the most sensitive parameter, responsible for about 90% of the changes in ice thickness and extension. To improve the rough turbulence models, it seems necessary to reinforce the experimental databases by precisely measuring the heat transfer variations induced by typical ice roughness. The present work shows that the precise knowledge of the heat transfers associated to a specific roughness pattern is key to ensuring reliable numerical icing predictions. Future investigations will be conducted to apply the same methodology to an airfoil.