Data-Driven Calibration of Rough Heat Transfer Prediction Using Bayesian Inversion and Genetic Algorithm

Abstract

The prediction of heat transfers in Reynolds-Averaged Navier–Stokes (RANS) simulations requires corrections for rough surfaces. The turbulence models are adapted to cope with surface roughness impacting the near-wall behaviour compared to a smooth surface. These adjustments in the models correctly predict the skin friction but create a tendency to overpredict the heat transfers compared to experiments. These overpredictions require the use of an additional thermal correction model to lower the heat transfers. Finding the correct numerical parameters to best fit the experimental results is non-trivial, since roughness patterns are often irregular. The objective of this paper is to develop a methodology to calibrate the roughness parameters for a thermal correction model for a rough curved channel test case. First, the design of the experiments allows the generation of metamodels for the prediction of the heat transfer coefficients. The polynomial chaos expansion approach is used to create the metamodels. The metamodels are then successively used with a Bayesian inversion and a genetic algorithm method to estimate the best set of roughness parameters to fit the available experimental results. Both calibrations are compared to assess their strengths and weaknesses. Starting with unknown roughness parameters, this methodology allows calibrating them and obtaining between 4.7% and 10% of average discrepancy between the calibrated RANS heat transfer prediction and the experimental results. The methodology is promising, showing the ability to finely select the roughness parameters to input in the numerical model to fit the experimental heat transfer, without an a priori knowledge of the actual roughness pattern.

1. Introduction

The numerical simulation of airflow above a rough wall differs from a classical smooth wall situation. The surface roughness creates a different near-wall behaviour, especially for the skin friction and the heat transfers. This led to the modification of the classical turbulence models to account for the roughness elements, like in [1] for the Spalart–Allmaras model. The main characteristics of the rough extension of the Spalart–Allmaras model reside in setting a non-zero turbulence viscosity at the wall. Usually, these adapted models are specifically designed to predict the correct skin friction coefficients based on experimental benchmark cases. This attention given to the skin friction coefficients alone has a drawback: the predicted heat fluxes are higher than the experimentally measured ones. An additional thermal correction model is thus required to adjust heat fluxes. To correct the heat flux behaviour, [2] suggested increasing the turbulent Prandtl number close to the wall. Further thermal correction model developments carried out by [3] and recently by [4] continued to aim at a turbulent Prandtl increase.

The thermal correction models take as input various parameters describing the roughness pattern of the surface. These parameters can be physically measurable, such as the roughness height, or numerical, such as the equivalent roughness. For instance, the thermal correction of [3] takes three input parameters and the two-parameter Prandtl correction (2PP) of [4], takes two. The latter 2PP thermal correction model is the one used in the present study. Previous works highlighted that the heat flux is sensitive to those roughness parameters and to the thermal correction model itself [4,5]. Usually, the equivalent roughness is computed from the roughness shape using empirical relations such as [6]. Those relations were built using manufactured regular roughness patterns, such as regularly spaced cones, hemispheres or pyramids. In many real-life situations, for example in aircraft icing [7], roughness patterns are irregular and present several uncertainties. Therefore, classical empirical relations can fail to correctly describe a roughness pattern in such situations.

To avoid the dependence on empirical relations in rough heat transfer CFD simulations, a data-driven approach can be used, provided that experimental data are available. This data-driven approach allows the recovery of the roughness parameters to input in the numerical model by processing the experimental data, without an a priori knowledge of the real roughness pattern. Machine learning and data-driven techniques have already been successfully applied in CFD, for example for the calibration of the turbulence model’s constants in [8], highlighting its potential.

Usually, the data-driven approach relies on metamodels. One of the most common metamodel families is the polynomial chaos expansion (PCE) [9]. PCE metamodels are commonly used in problems where uncertainty propagation is studied, like in [10] and [11]. Most of the time, the PCE approach allows building efficient metamodels from a CFD database, and these are used in lieu of the time-consuming CFD model [12] for complex analyses. The PCE metamodeling was successfully applied in heat transfer applications, for instance in [13] or [14], highlighting its suitability for the present study. Globally, the coupling between PCE metamodeling and CFD has been extensively developed and applied during the last two decade [15,16]. Once the metamodels are established, the sensitivity analysis can be performed, usually using the Sobol indices [17]. This combination of PCE metamodels with Sobol sensitivity indices was previously used in a related work [5]. The Sobol indices allows the classification of the uncertain inputs according to their influence on the particular output of interest of a study using a variance-based approach. These commonly used indices are regularly employed in CFD sensitivity analysis [18]. Following the definition of the metamodels, the calibration of uncertain roughness parameters can be performed. In the present work, two calibration methods will be assessed and compared: the Bayesian inversion and the genetic algorithm. The Bayesian inversion is a commonly used tool allowing model calibration and fine-tuning of numerical parameters [19]. The Bayesian inversion transforms the initial uncalibrated input’s distribution into a refined posterior distribution, whose peak corresponds to the calibrated value [20]. This approach, using probability distributions instead of fixed values, allows flexibility in the interpretation of the calibration results. This calibration methodology was used in the context of the constant calibration of turbulence models [21]. Fine tuning and optimization tasks suit the application of Bayesian inversion. The authors of [22] used Bayesian inversion for the shape optimization of a wall to obtain a prescribed pressure gradient distribution. This application presents similarities with the present work aiming at optimizing the surface roughness to obtain a given heat flux distribution. On the other side, classical calibration methods (i.e., non-Bayesian), such as the genetic algorithm, usually work in a non-probabilistic way, making their calibrated results strictly defined, without controlled flexibility. Classical calibration methods imply more intervention from the human modeller in post-processing, making the final results sensitive to the modeller’s experience, expertise and biases [23]. A genetic algorithm calibrates the problem using an approach similar to biological evolution [24]. The final calibrated parameters are estimated by performing a method using selections, mutations and crossings between “generations” of parameters samples, like in the natural selection process. This type of calibration with genetic algorithm was already used in CFD applications [25], and more specifically for optimization, like in [26] or [27]. This use of genetic algorithms in the CFD field highlights the suitability of the method for the current application. The coupling between the PCE metamodeling and the optimization by Bayesian inversion or genetic algorithm was already used in numerical heat transfer applications [28], showing the potential of the methodology. Additionally, in case of a low trustworthiness metamodel, coupling it with a Bayesian analysis allows one to consider the metamodel’s uncertainty in the whole study [29].

The objective of this paper is to set up a methodology to calibrate the roughness parameters to input into the 2PP thermal correction model to match the experimental results of a given test case. More specifically, a rough curved channel geometry, inspired by [2], is used to set up the design of experiment (DOE) for heat transfers. Next, the DOE allows the generation of PCE metamodels predicting key features of the heat transfers on the wall. The next step is the sensitivity study, using the Sobol indices to identify the most sensitive parameters. Finally, the PCE metamodels are used in a calibration procedure to establish the roughness input parameters with the aim of retrieving the experimental heat transfer distribution in the curved channel. The calibration is done successively with a Bayesian approach and a genetic algorithm to compare both methods on the results discrepancies and practical application.

First, the test case geometry and setup will be described. Second, the 2PP thermal correction model is detailed, highlighting the role of the roughness parameters in the simulation. Next, the DOE construction is depicted prior to the description of the PCE metamodeling, the sensitivity study and the calibration procedures. Finally, the results are shown and compared, highlighting between 4.7% and 5.4% of average discrepancy after the Bayesian calibration compared to the experimental results. On the other side, the genetic algorithm approach gives between 5.7% and 10% of discrepancies.

2. Test Case Geometry and Setup

A specific curved channel test case is used for the application of the depicted methodology. This section will give details about the geometry used in the present study, along with the RANS setup to perform the flow simulation. The geometry and computational domain is detailed, prior to the domain mesh and numerical RANS configuration.

2.1. Physical Geometry and Boundary Conditions

The geometry used to apply the calibration work discussed in this paper is a curved channel experimentally studied by [30]. This test case was numerically studied by [2] few years later. The channel has a straight section followed by a curved one and another straight section. The channel height is 0.07 m while the overall curvilinear length is about 1.25 m. More specifically, the zone of interest is the rough bottom floor of the channel where the heat transfer coefficient will be monitored. Figure 1 illustrates the computational domain.

The curved part of the channel’s floor, whose curvilinear length is 0.398 m, is an arc of radius 1.200 m and an angle of 19°. The freestream values are a velocity magnitude of 40 m/s, an inlet total pressure of 102,304 Pa and a total temperature of 288.95 K. The boundary conditions are also displayed on Figure 1. The floor is divided into three main zones: an initial smooth and unheated wall, the rough and heated study zone (including the curved portion), and a downstream smooth unheated zone. The top wall is entirely smooth and adiabatic. The heated zones are isothermal at 303.15 K and their roughness parameters are the objective of the calibration. The roughness parameters’ distribution will be detailed in Section 4.

2.2. Mesh and Numerical Setup

The geometry is discretized using quadrilateral elements in a structured mesh. The rough zone of interest is composed of 499 nodes in the stream-wise direction, while the channel height has 399 elements. Coarser meshes were tested, giving similar results on the benchmark test cases. The finest mesh was retained to cope with the various untested roughness patterns planned to be run during the sampling of the DOE. The entire computational grid has 274,512 quadrilaterals. The first cell height is about 3 µm, which allows a y⁺ below 1 for all the roughness ranges tested. The growth rate normal to the wall is 1.1. A close-up on the mesh near the floor is shown on Figure 2. The main outcomes of the mesh convergence study are detailed in the Appendix A.

The flow simulation is performed using the compressible RANS solver SU2 6.2 [31]. The Reynolds number based on the total floor length is $3.2\times {10}^6$ and the Mach number is 0.118. To accelerate the convergence, the CFL number is set to 10. Finally, the convective fluxes are discretized using a Roe scheme with MUSCL reconstruction [32]. The solver used includes in-house implementations of the rough modification of the Spalart–Allmaras turbulence model [1] and the addition of the 2PP [4] thermal correction model. The next section will detail the mathematical model making up the thermal correction model.

3. 2PP Thermal Correction Model

This section will give an overview of the equations of the 2PP thermal correction model. The 2PP thermal correction model was implemented into SU2, and aims at increasing the Prandtl number in the near-wall region to reduce the heat flux when the rough version of the Spalart–Allmaras turbulence model is used.

The 2PP thermal correction model achieves the goal of increasing the turbulent Prandtl number by computing an increment ΔPr_t. The model takes two roughness input parameters: the roughness height k (m) and the equivalent roughness k_s (m). Equation (1) allows the computation of ΔPr_t. In Equation (1), Pr is the laminar Prandtl number and d is the distance to the wall.

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

(1)

The roughness Reynolds number and the parameter g in Equation (1) are detailed in Equations (2) and (3).

$$R{e}_S=\frac{u_{\tau }{k}_s}{\nu }$$

(2)

$$\left\{\begin{array}{c}g=1 \mathrm{if}R{e}_s\ge 70\\ g=\frac{\ln \left(R{e}_s\right)-\ln \left(5\right)}{\ln \left(70\right)-\ln \left(5\right)} \mathrm{if}5<R{e}_s<70\\ g=0 \mathrm{if}R{e}_s\le 5\end{array}\right\}$$

(3)

In Equation (2), $u_{\tau }$ is the friction velocity, and ν is the air kinematic viscosity. Finally, Equations (1)–(3) highlighted that the roughness parameters have a direct impact on the Prandtl number correction, and thus on the predicted heat flux.

4. PCE Metamodeling

This section will depict the metamodeling process retained, using polynomial chaos expansion (PCE) metamodels. This step allows the creation of metamodels to predict the heat flux behaviour above any roughness pattern without the need for a complete CFD run. The first step of the metamodeling task is the creation of a numerical DOE. The PCE metamodels are then generated, and finally their accuracy is checked to ensure they are reliable enough for the study.

4.1. Design of Experiment (DOE)

To prepare the metamodeling and calibration steps, a numerical DOE of heat transfers for various roughness patterns is needed. First, the distribution of the input parameters k and k_s/k are defined to set the ranges for the sampling. Note that the ratio k_s/k is used instead of k_s alone, since it will allow one to directly evaluate the relation between the roughness height and the equivalent roughness. The present work is included in the broader scope of in-flight aircraft icing. Therefore, the typical ranges of variation of k and k_s/k are obtained from the icing literature [33,34]. These ranges are wide enough to ensure that they are suitable for the current study, which is not specifically a simulation in icing conditions. The compilation of the distribution of all the input parameters is given in

Table 1. Distribution of the input parameters.

Parameter	Minimum	Maximum	Distribution
k (mm)	0.41	4.32	Uniform
Ratio ks/k	0.2	6.5	Uniform

. The distributions are chosen as uniform since there is no a priori knowledge of the experimental roughness pattern.

Following the distributions of Table 1, a sampling is done using the Latin hypercube sampling [35]. The sample size is defined according to the literature [36] for a response surface and gives, with an oversampling factor, 120 samples. Following the sampling, 120 CFD simulations of the curved channel are run. This allows the construction of a numerical heat transfer database, needed for the next metamodeling step.

4.2. Metamodels Generation

Once the heat transfer DOE is set up and ran, the metamodeling tool uses its output database to estimate a mathematical relation between the roughness parameter inputs and the heat flux output. For this study, polynomial chaos expansion (PCE) metamodels were chosen [9] and generated using the UQLab tool [9]. The choice of this type of metamodel is motivated by its wide use in uncertainty quantification for CFD and aerodynamic applications [37]. The general form of a PCE metamodel is given by Equation (4).

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

(4)

In Equation (4), Y_i is the output of interest, X = (X₁, X₂) is the input parameter vector, M_i is the corresponding PCE metamodel defined by its coefficients $y_{\alpha }$ and the multivariate polynomials ${\psi }_{\alpha }$ of the decomposition. α = (α₁, α₂) is the multi-index with two components (since there are two input parameters; see Table 1). The multivariate polynomials ${\psi }_{\alpha }$ are obtained as the tensor product of the two (in the present application) univariate basis polynomials φ (Equation (5)).

$${\psi }_{\alpha }\left(X\right)={\displaystyle \prod }_{i=1}^2{\phi }_{{\alpha }_i}(X_i)$$

(5)

The input parameters having a uniform distribution, the univariate basis polynomials are Legendre polynomials of indices α_i [9].

Three metamodels are generated, as listed in

Table 2. The metamodels created.

Metamodel	Output(s) of Interest
M1	hc at the starting point of the rough zone, W/m2K
M2	Mean relative error with experimental hc (%)
M3	hc values at N equally spaced locations along the rough zone (multi-output), W/m2K

where h_c is the heat transfer coefficient in W/m²K.

The second metamodel, M₂, predicts the mean relative error compared to the experimental results of [2]. The mean relative error ε is the average of the relative errors taken on every mesh point (Equation (6)).

$$\epsilon =\frac{1}{N_{points}}{\displaystyle \sum }_{i=1}^{N_{points}}\left|\frac{h_{c,CFDi}-h_{c, EXPi}}{h_{c, EXPi}}\right|$$

(6)

In Equation (6), N_points is the number of mesh points in the study zone, and h_c,CFDi and h_c,EXPi are the CFD predicted h_c at point i and the experimental h_c at point i, respectively.

M₃ is a multi-output metamodel evaluating in one estimation N values of h_c along the wall. Setting N to a high value allows one to predict the complete h_c distribution on the entire rough wall. The metamodel M₁ is separated from the metamodel M₃: this separation allows one to compare the cases where the starting h_c is calibrated alone (M₁) and where several locations are simultaneously calibrated (M₃).

Once a metamodel is created, its evaluation on the inputs of the DOE allows the comparison of its prediction with the actual CFD prediction obtained when setting up the DOE. Doing a linear regression between the PCE prediction and the CFD output allows the computation of the R² coefficient (Equation (7)).

$$R^2=1-\frac{{\sum \left(Y_{CFD}-Y_{PCE}\right)}^2}{{\sum (Y_{CFD}-{\overline{Y}}_{CFD})}^2}$$

(7)

In Equation (7), Y_CFD and Y_PCE are the CFD and PCE predictions, respectively. ${\overline{Y}}_{CFD}$ is the mean value of the output of interest (CFD). An R² coefficient close to 1 ensures a PCE metamodel with a good accuracy, since it predicts outputs close to what the full CFD simulation gives. Once the metamodels are established, the next step is to use them for the sensitivity study and the calibration purpose.

5. Sensitivity Study

The sensitivity study allows the identification of the most sensitive parameter(s) in the model. Using the previously described PCE metamodels, the sensitivity analysis computes the Sobol indices, described in [17]. These indices allow the classification of the input parameters from the most to the least sensitive. The Sobol index is a ratio of variances, taking values between 0 and 1. For input parameters i and j, the first and second order Sobol indices are defined by Equations (8) and (9), respectively.

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

(8)

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

(9)

In Equations (8) and (9), Y is the output of interest, V is the variance operator, E the mean value, and the notation Y|X_i denotes the output of interest when the ith input parameter is fixed. The total Sobol index, which is monitored in the present study, for the ith input parameter for a generic three parameters study is given by Equation (10).

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

(10)

The total Sobol index gives the contribution of the variability of the ith input and its interactions with the other inputs (here the jth and kth) on the response sensitivity. A total Sobol index close to zero means that the ith parameter does not contribute much to the variability of the studied output. According to the classification criteria suggested by [38], an input parameter with a total Sobol index:

• above 80% is very important;
• between 50% and 80% is important;
• between 30% and 50% is unimportant;
• below 30% is negligible.

6. Model Calibration

Observing experimental results without an a priori knowledge of the roughness pattern is not helpful to precisely extract the roughness parameters. This task is even more non-trivial since the roughness parameters to input in the model vary depending on the thermal correction model chosen. The calibration is intended to estimate those roughness parameters by working on the PCE metamodels previously created. The Bayesian inversion calibration is first described prior to the genetic algorithm approach.

6.1. Bayesian Inversion Calibration

Based on the Bayes theorem (Equation (11)), the principle of the Bayesian inversion is to get the posterior distributions π(θ|X_i) of input X_i, based on the assumed prior distributions (see Table 1) π(θ) and on the information provided by the experimental data. Here, θ denotes the distribution parameters and the | symbol denotes the conditional dependence. The posterior distributions are the distributions of the input parameters knowing the information brought by the experimental data.

$$\pi \left(\theta |X_i\right)=\frac{\pi \left(X_i|\theta \right)·\pi \left(\theta \right)}{\pi \left(X_i\right)}$$

(11)

In Equation (11), π(X_i|θ) is called the likelihood function and π(X_i) is seen as a normalization constant called the marginal likelihood. The new posterior distributions are generally not uniform anymore and allow identifying the input parameter values that will produce an output that best fits the experimental data. Additionally, a discrepancy between the PCE output and the CFD prediction is given to the solver, along with the experimental observations themselves. For instance, a discrepancy between 0 W/m²K and 15 W/m²K means that up to 15 W/m²K of difference between the PCE and CFD predictions is expected. Numerically speaking, the Bayesian inversion is performed with the UQLab tool [19]. The computation of the posterior distribution is made using a Markov chain Monte Carlo (MCMC) algorithm. The samplers used in the study are the affine invariant ensemble algorithm (AIES) or the Metropolis–Hastings (MH) algorithm. For the purposes of the work, the MCMC solver is tuned to perform 70,000 iterations and generates 15 chains. The Bayesian solver setup for each metamodel is summarized in

Table 3. Bayesian solver settings.

Metamodel Calibrated	Objective of the Calibration	Sampler	Discrepancy	Experimental Observation Supplied
M1	Recovering the same starting value of hc	AIES	Uniform [0; 15] W/m2K	255.1 W/m2K
M2	Having a mean relative error with experimental hc of 0%	AIES	Uniform [0; 5]%	0%
M3	Recovering the same hc values at the N equally spaced locations along the rough zone	MH	Uniform [0; 15] W/m2K	Experimental hc at the N locations (W/m2K)

. For the metamodel M₂ (predicting the mean relative error with experimental data), the experimental observation is 0% of mean relative error. This means that the objective of the calibration is to obtain a mean relative error close to 0%.

After having obtained the posterior distributions, the mean value or the maximum a posteriori (MAP) of these distributions can be estimated. These two statistics can be used to estimate which single values of input parameters lead to a model response close to the experimental data. After evaluating the PCE metamodels and the complete CFD simulation with the calibrated parameters, it is possible to assess the success—or not—of the calibration procedure. One can estimate if the model response is close to the experimental data. Additionally, the distributions of the model predictions using the posterior distributions should present less uncertainty (reduced variance) when using the prior distributions.

6.2. Calibration Using a Genetic Algorithm

The Matlab framework has a built-in genetic algorithm that is used in the study [39,40]. Its purpose is to find the set of input parameters that minimize an objective function. The default settings provided in the Matlab framework provide a straightforward application for the present case. The objective functions supplied to the genetic algorithm are summed up in

Table 4. Objective functions for calibration with the genetic algorithm.

Metamodel Calibrated	Objective of the Calibration	Objective Function Used
M1	Recovering the same starting value of hc	|M1 − 255.1|
M2	Having a mean relative error with experimental hc of 0%	|M2|
M3	Recovering the same hc values at the N equally spaced locations along the rough zone	|M3[i] − hc[i]| i = 1…N

.

Note that since the metamodel M₃ outputs multiple (N) values, the genetic algorithm is executed N times. Each one of the executions gives a set of calibrated parameters and the final retained parameters are obtained by doing an average or a weighted average over the N solutions. Most of the time, the simple averaging produces unsatisfactory results. The weighted average, with appropriate weight(s) on the most relevant output(s), allows one to improve the results, as shown in the next section.

7. Results

This section will display the results obtained all through the process described in the previous sections. First, the CFD results obtained initially before calibration are shown to illustrate the baseline results. Following this, the characteristics and precision of the metamodels are shown prior to the sensitivity indices. Finally, the calibration results for both Bayesian inversion and genetic algorithm are displayed, prior to the comparison of both methods.

7.1. CFD Results before Calibration

Prior to the metamodeling and calibration process, the curved channel test case, as described in Section 2.1 and Section 2.2, is run alone to verify the CFD settings and mesh. This initial simulation, called here the baseline simulation, is the same as the one run for the mesh study with the fine mesh. The roughness height is k = 0.5 mm and the equivalent roughness is k_s = 1.55 mm. This baseline simulation also provides useful information about the convergence of the CFD configuration used. The iterative convergence residual curves are displayed on Figure 3.

Figure 3 shows that after the prescribed number of iterations (125,000), the residuals are stable and decreased by at least three orders of magnitude. The configuration used has a satisfactory convergence behaviour which confirms the choice of the numerical setup. The heat transfer results obtained with these settings are plotted on Figure 4.

From Figure 4, it is possible to notice a relatively poor agreement with the experimental data. Finding the best roughness parameters is non-trivial, and it legitimates the use of the data-driven calibration approach used to match the literature.

7.2. Visualization of the DOE outputs

After a sampling of the roughness parameters with the Latin hypercube method, 120 CFD simulations were run. The heat transfer databases obtained from the DOE are plotted in Figure 5, where each one of the black curves is the heat transfer coefficient of one simulation of the DOE. On Figure 5, the experimental data, in red, from the literature are shown.

By observing the experimental results, it is possible to notice oscillations resulting from the experimental uncertainty. Additionally, all the numerical simulations exhibit a sudden bump at x = 0.46 m, where the curved part of the channel starts. This local increase in the heat transfer coefficient is due to the flow acceleration in the curve portion. Studies in [30] quantified this increase in the heat transfer coefficient value at the start of the curve zone to be about 3%. Figure 5 shows that the target experimental results are included in the envelope defined by the sampling CFD simulations. This observation shows that the initial range of roughness parameters chosen contains the a priori unknown experimental values.

7.3. Characteristics of the Metamodels and Accuracy

The previous DOE is used as a basis for the PCE metamodel generation. Three metamodels are generated: one to predict the initial value of h_c along the rough zone, one to predict the mean relative error with experimental results and a third one to predict N values of h_c along the wall (see Table 2). For each metamodel, the resulting PCE is described with the polynomial degree p_PCE and the number of terms in the expression (Equation (4)). Additionally, the R² coefficient is calculated to assess the metamodel accuracy (Equation (5)). The values obtained are gathered in

Table 5. R² coefficient for each metamodel.

Metamodel	Output of Interest	PCE Degree pPCE	Number of Terms in Equation (4)	R2 Coefficient
M1	hc at the starting point of the rough zone, W/m2K	10	66	0.99994
M2	Mean relative error with experimental hc (%)	10	66	0.99962
M3	hc values at 6 equally spaced locations along the rough zone (multi-output), W/m2K	Y1:10	66	0.99994
		Y2:9	55	0.99993
		Y3:12	91	0.99999
		Y4:10	66	0.99992
		Y5:12	91	0.99999
		Y6:12	91	0.99997

. For the present study and for the rest of the paper, the metamodel M₃ is used with N = 6, meaning it computes six values of h_c regularly spaced along the wall. Note that the number of terms in Equation (4), for a two input parameters study, is equal to $\frac{(p_{PCE}+2)!}{p_{PCE}!2!}$.

Table 5 shows that the PCEs obtained have a degree between 9 and 12, with the corresponding polynomial expressions having between 55 and 91 terms. All regression coefficients are above 99.9%, meaning an excellent agreement between the CFD results and the PCE-predicted results on the same sample. For comparison, [41] performed the same type of uncertainty quantification analysis with R² coefficients as low as 94.6%. For graphical visualization, the regression for M₂ (with R² of 99.96%) is plotted on Figure 6 where Y_PCE and Y_CFD are the mean errors in percentage with the literature, as predicted by PCE and CFD, respectively. For synthesis and concision, only the regression for M₂ is displayed, since it is the worst among the three (see Table 5).

Figure 6 shows that the results are close to the identity line, visually confirming the good value of the R² regression coefficient. This R² assessment shows that the metamodels generated are accurate and reliable enough to be used in the present application.

7.4. Sensitivity Study

Using the PCE metamodels, the sensitivity study allows the Sobol indices to be computed (see Equation (8)).

Table 6. Total Sobol indices.

Metamodel	Total Sobol Indices
M1	k: 0.1445 ks/k: 0.8868
M2	k: 0.3061 ks/k: 0.9772
M3	k: 0.1167 ks/k: 0.9061

gathers the total Sobol index values for each metamodel.

Table 6 shows that the roughness height k is the least sensitive parameter in each case, with a total index between 11% and 30%. The ratio k_s/k is the most dominant parameter of influence, with 88% to 97% sensitivity. According to the classification made by [38], k is a negligible parameter in the model sensitivity while k_s/k is a very important parameter. Therefore, the relation between k_s and k is more critical for the model output than their absolute values in millimetres.

7.5. Bayesian Inversion Calibration

The metamodels built are used as forward models for the MCMC algorithms, evaluated at each iteration of the chain. The target of the calibrations is to retrieve the features of the experimental results by estimating the best roughness parameters. These targets were listed earlier in Table 3. The target values (last column of Table 3) are inputted into the Bayesian module and the computation is carried out on the corresponding metamodel. For metamodels M₁ and M₂, the calibrated roughness parameters retained are the mean values of the posterior distribution. For metamodel M₃, the maximum a posteriori (MAP) is retained. The calibrated roughness parameters retained are listed in

Table 7. Calibrated roughness parameters (Bayesian inversion).

Calibrated Metamodel	Values Retained	Final Calibrated Roughness Parameters
M1	Mean	k = 2.2 mm ks = 6.4 mm
M2	Mean	k = 1.6 mm ks = 4.2 mm
M3	MAP	k = 1.8 mm ks = 5.0 mm

, along with the value (mean or MAP) chosen.

The values in Table 7 are different from the ones tested in the baseline simulation, what explains the poor agreement with the literature prior to the calibration. Before assessing the quality of the calibrations, it is interesting to look at the main power of the Bayesian analysis: transforming the prior probability distributions of the inputs (uniform, see Table 1) into the calibrated posterior distributions, which are closer to Gaussian distributions. Figure 7 and Figure 8 illustrate as an example the prior and posterior distributions for the M₃ calibration for both k and k_s/k roughness parameters, respectively.

In Figure 7 and Figure 8, it is possible to see the difference from the uncalibrated distributions (prior) to the posterior calibrated distributions. The calibrated values are clearly visible since the distributions are less spread and present peaks. Figure 7 and Figure 8 also allow visualizing graphically the meaning of the mean and MAP values: the mean is the mid value between the distribution’s lower and upper limits, while the MAP value corresponds to the highest peak. From the posterior distributions, it is possible to recover the calibrated values of Table 7 for metamodel M₃: k = 1.8 mm and k_s/k = 2.8 (i.e., k_s = 5.0 mm). The observations highlight that the ratio k_s/k (Figure 8b) is finely calibrated with a smaller uncertainty compared to the roughness height k (Figure 7b). The posterior distribution for the ratio k_s/k exhibits a narrower peak, typical of a small variance.

For the metamodel M₁ and M₂, the calibrated roughness parameters from Table 7 are inputted into the CFD solver and the simulation is run to verify the new heat transfer obtained after calibration. Figure 9 shows the heat transfer using the calibrated roughness parameters for M₁ and M₂.

Figure 9 shows that the calibrated results have a better agreement with the experimental results compared to the baseline simulation of Figure 4. The calibration for M₁ allows one to obtain an average relative error with the experimental data of 4.7%. The relative error is computed as the average for the entire rough zone among all the grid points. The experimental data being available at about only 25 locations, the experimental values are interpolated at the grid points to allow the relative error calculation. The calibration for M₂ presents similar errors with the experimental data with 4.8% of relative mean error. Globally, both calibrations present satisfactory results, showing less than 5% of error on average compared to the experimental results. The results for the metamodel M₁ are higher than the one for the metamodel M₂, which is due to the higher roughness parameters in the case of M₁. This observation is in accordance with the usual experimental observation, where higher roughness elements lead to an enhanced heat transfer [2].

Figure 10 illustrates the calibrated results using the results for the metamodel M₃. Since M₃ outputs h_c values at different locations, it is possible to use the calibrated inputs to plot both PCE metamodel and CFD predictions.

Figure 10 is informative on several aspects: first, the PCE prediction using the calibrated roughness parameters presents a satisfactory agreement with the experimental data. Second, and most importantly, the CFD prediction is close to the PCE prediction. It allows assessing whether the PCE metamodel is reliable enough to predict a solution, and makes it a less costly and time-consuming tool which can replace the CFD on that particular application. Finally, the calibrated CFD results present a good agreement with the experimental data. The mean relative error with the experimental curve is 5.4%.

7.6. Genetic Algorithm Calibration

The same three metamodels are now calibrated with the genetic algorithm. The objectives of the calibrations remain the same, and the objective functions were listed earlier in Table 4. The calibrated results obtained are displayed below in

Table 8. Calibrated roughness parameters (genetic algorithm).

Calibrated Metamodel	Values Retained	Final Calibrated Roughness Parameters
M1		k = 1.9 mm ks = 3.1 mm
M2		k = 4.3 mm ks = 8.2 mm
M3	Average among the N values	k = 3.0 mm ks = 8.3 mm
M3	Weighted average among the N values	k = 2.9 mm ks = 7.1 mm

. For the metamodel M₃, the genetic algorithm outputs N calibrated parameters corresponding to each of the N locations. The final retained parameters are computed by average or weighted average (where the initial point at the start of the rough zone has a weight of three).

Figure 11 shows the CFD results using the calibrated parameters for the metamodels M₁ and M₂.

Figure 11 highlights that the agreement with the experimental data is better at the beginning of the rough zone and tends to present higher discrepancy at the end of the rough zone. The mean relative error for the metamodel M₁ is 8.2% and is 5.7% for metamodel M₂.

Figure 12 displays the PCE and CFD predictions using the calibrated values of metamodel M₃. The calibration shows better agreement in the middle and at the end of the wall. This feature has a drawback: the agreement at the beginning is not as good as for the metamodels M₁ and M₂. Taking the average roughness parameters gives a large overestimation of h_c at the start of the rough zone. Performing a weighted average, giving a weight of three to the first point, tends to lower the starting h_c value and globally lowers the values everywhere.

7.7. Comparison of Both Calibration Methods

After calibrating the metamodels using both Bayesian inversion and genetic algorithm approaches, the first comparison is the capacity for predicting the experimental data. This comparison is made by computing the mean relative error compared with experimental data.

Table 9. Errors with experimental data for each calibration.

Calibrated Metamodel	hc Mean Relative Error with Experimental Data
	Bayesian Inversion	Genetic Algorithm *
M1	4.7%	8.2%
M2	4.8%	5.7%
M3	5.4%	Avg: 10% W-Avg: 7.0%

* For M₃, Avg: Average among the N outputs; W-Avg: Weighted average.

gathered together the mean relative errors for both approaches.

The Bayesian calibration allows an agreement between the numerical results and the experimental data with an error about 5%. For the genetic algorithm, the errors are higher: they are between 5.7% and 10%. These values show an agreement not as good as the one obtained with the Bayesian inversion, despite displaying interesting behaviour at the beginning of the rough zone for the metamodels M₁ and M₂ (see Figure 11). By purely looking at the values in Table 9, the Bayesian inversion provides more satisfying and consistent results compared to the genetic algorithm. Note that the uncertainty in the experimental results was not considered. Given the oscillations observed on the experimental curve, one can estimate that the uncertainty in the experimental data is between 6% and 10%.

Nevertheless, by comparing Figure 9, Figure 10, Figure 11 and Figure 12, both methods have their own strengths and weaknesses. The calibrations for the metamodels M₁ and M₂ present a better agreement at the beginning of the rough zone for the genetic algorithm calibration (Figure 11) compared to the Bayesian inversion (Figure 9). When calibrating the starting value of h_c with the metamodel M₁, the genetic algorithm performs better since the starting value of h_c after calibration agrees better with the experimental data compared to the Bayesian inversion. On the other hand, the Bayesian inversion presents a better agreement between the calibrated results and the experimental data in the middle and at the end of the rough zone. For the metamodel M₃, the trend is the opposite: the Bayesian inversion (Figure 10) performs better at the beginning of the rough zone while the genetic algorithm (Figure 12) exhibits a better agreement at the end of the rough zone.

On usage, the Bayesian inversion is more automated and gives a more flexible output. The posterior distribution gives an idea of the uncertainty of the calibrated parameters. In the meantime, the output of the genetic algorithm is more rigid, since apart from the calibrated value, no distribution is given. In the case of a multi-output metamodel (M₃), the Bayesian inversion calibrates the model by working simultaneously on all outputs, thus estimating correlated calibrated parameters. The genetic algorithm works on each output independently, and requires a human intervention to establish the correlation between the outputs, like in the case of the weighted-average metamodel M₃. Modeller intervention can have a large impact on the results. Table 9 shows that changing the weight of the first output of the metamodel M₃ from one to three in the averaging decreases the error from 10% to 7%. An optimization of the weights in the future could improve the results even further.

The novelty of this work resides in the joint application of PCE metamodeling with a calibration method in the field of rough heat transfer prediction. Usually, unknown roughness patterns lead to large uncertainty in the aerothermal behaviour observed. With the methodology presented here, it is possible to reduce the lack of knowledge about a given experimental roughness pattern. As an extension, the present methodology can be applied to aircraft icing, where aerothermal and ice accretion behaviour depend strongly on the initial and uncertain (and unknown) roughness pattern [5,42]. Finely calibrating the roughness parameters in an icing simulation will reduce one source of uncertainty that affects the final prediction of the ice shape.

8. Conclusions

The paper presents a methodology to perform a calibration of roughness parameters aiming at approaching experimental heat transfers. Two methods of calibration were compared: the Bayesian inversion and the genetic algorithm. Furthermore, this methodology was applied to a rough curved channel test case to illustrate its capacity. Starting with an unknown experimental roughness pattern, the procedure allowed the recovery of roughness parameters, providing between 4.7% and 5.4% of discrepancies with the experimental heat transfer when using Bayesian inversion calibration. In the case of the genetic algorithm, the agreement was above 5.7%, reaching 10%. The Bayesian inversion handled better than the genetic algorithm in the present application since its calibrated results fit the experimental data better. Furthermore, the Bayesian inversion behaves better in the case of multi-output metamodels, since it does not require manual tuning to account for each output contribution. This building of this methodology was achieved by combining polynomial chaos expansion (PCE) metamodeling and calibration techniques. The methodology showed high accuracy with respect to the PCE metamodels, confirming their suitability in the current application. The sensitivity analysis using the Sobol sensitivity indices highlighted that for this test case, the relation between the roughness height and the equivalent roughness plays a bigger role than the roughness height alone. The data-driven approach showed its suitability in CFD applications. In the case of unknown roughness patterns, it allows one to select the numerical parameters to input into a CFD simulation to retrieve the experimental data. Future extensions of the work will allow one to calibrate the roughness parameters not only to fit the heat transfer value but related quantities such as ice shape geometry in in-flight icing simulations.