CFD investigation on the maximum lift coefficient degradation of rough airfoils

1.   Introduction

Ice protection systems reduce risks of accident but increase aircraft energy consumption. To save energy, deicing systems are the best choice, but in some moments, aircrafts must fly with small amounts of ice. The small amount of ice is modeled as roughness, which degrades the aerodynamic performance. Therefore, a better understanding of roughness on wing performance can help design energy-efficient deicing systems. Experiments proved that even in its very earlier stages (increased roughness), icing could cause a reduction of 25% in the maximum lift, and an increase of 90% in drag of an aspect ratio 6 wing ^[1]. In the same way, experiments show that standard leading-edge roughness creates similar performance degradation on an airfoil, with an increase of around 60% of the drag for a NACA 0015 ^[2].

With the evolution of computers, numerical simulation gained popularity due to its reduced cost and satisfactory results ^[3]. In addition, only numerical methods are able to explore, safely, all icing conditions required for certification ^[4]. In this way, different approaches have been proposed such as "highest fidelity", "discrete element", and "equivalent sand grain" ^[5]. As reported in ^[6], the first approach consists of computing the flow around the roughness elements. However, due to its high computational cost, this approach is presently done only at a research level. In the second approach, roughness effects are averaged at a higher level from the reference surface, which requires the modification of the mean flow equations to account for blockage effects caused by the roughness elements. However, Aupoix emphasizes that this modification precludes the use of this technique in commonly used RANS (Reynolds-Averaged Navier-Stokes) solvers. The last approach, proposed by schlichting ^[7], is the most popular and, according to Aupoix, it is the only affordable one for industrial applications. In this technique, the equivalent sand-grain roughness (${k_s}$) represents the height of a sand grain from Nikuradse's experiments, and reproduce the drag increase in the fully rough regime. With the purpose of reducing the uncertainty in the determination of the ${k_s}$, dirling proposed a correlation ^[8] in which topological features are translated into the ${k_s}$. This correlation became very popular, and it is frequently cited in the literature ^[9]. Along the years, different roughness implementations based on the ${k_s}$ have been proposed. For the well-established and validated Spalart-Allmaras (SA) turbulence model, two different extensions were independently developed by ONERA and Boeing ^[10]. Both models utilize the equivalent sand grain approach to reproduce the drag increase in the turbulent model.

Scaled-down airfoils are commonly used to evaluate ice formation in wind-tunnel tests, which poses a concern on whether the results produced would be useful to represent full scale flight conditions. Another concern is the computational cost of CFD simulations, which is increased at higher Reynolds numbers due to wider range of eddy sizes to be solved as well as finer grids requirements ^[11]. Therefore, lower Re simulations are desired for industrial applications. On the other hand, when roughness is considered, result obtained at reduced Reynolds number are not meaningful or useful for higher Reynolds number applications ^[12]. In this way, it is desired to find the minimum Reynolds number that are useful for higher Reynolds number simulations in which roughness is considered.

Moreover, the process of in-flight ice formation might result in different roughness heights. As seen in the literature, the maximum lift coefficient $C_{l(\max )}$ is usually reduced as the roughness height and density are increased ^[13,14]. Thus, an evaluation of roughness effects considering different equivalent sand grain roughness heights would help to build a correlation to estimate aircraft performance degradation translated by the deterioration of the maximum lift coefficient.

In this work, we propose a correlation to predict the degradation of the maximum lift coefficient caused by evenly distributed roughness effects on flows over two airfoils, a NACA 0012 and a model 5–6. In addition, a second correlation is proposed to find the minimum Reynolds number that are useful for higher Reynolds number applications when roughness is considered.

The paper is structured as follows. Section 2 introduce the mathematical and numerical models used herein. First, the widely used RANS model is presented followed by the turbulence model applied in our simulations, Spalart-Allmaras (SA). Next, the extension of SA model that considers roughness is presented. Finally, the numerical model utilized in this paper is described. In Section 3, the verification and validation of the numerical approach is presented. Two airfoils are used in this process, a NACA 0012 and a model 5–6. First, the behavior of the flow over a smooth NACA 0012 is investigated to confirm whether the implementation has no influence on the original model when roughness is not activated. Then, roughness is activated, and the flow over a rough NACA 0012 is investigated. Estimations of lift coefficients and velocity profiles inside the boundary layer are evaluated and compared to numerical and experimental results. In Section 4 (Results and discussion), investigations on the maximum lift coefficient reduction caused by different equivalent sand grain roughness heights are performed. Then, the influence of Reynolds numbers on the maximum lift coefficient reduction is evaluated, and correlations are proposed.

2.   Materials and method

2.1.   Reynolds-averaged navier-stokes

The flow is modeled with the compressible Reynolds-Averaged Navier-Stokes equations ^[15]. According to ^[16], when the density is not constant (compressible flows), it is recommended to apply Reynolds averaging for density and pressure, whereas Favre averaging is used for the other variables, which results in:

where $\bar \rho$ represents the mean density, ${\tilde \upsilon _i}$ is the mean velocity, $\bar p$ is the mean pressure, h is the enthalpy, $k$ is the thermal conductivity coefficient, and K is the turbulent kinetic energy. Here, the bar denotes the Reynolds-averaged values, the tilde represents the Favre-averaged values, and the double prime is used for fluctuating parts. The viscous stress tensor, ${\tilde \tau _{ij}}$ is defined as:

where $\lambda$ represents the second viscosity coefficient, ${\mu _{dyn}}$ is the dynamic viscosity, and ${\delta _{ij}}$ is the Kronecker symbol. The dynamic viscosity can be calculated based on the Sutherland law and reference values, ${\mu _{ref}} = 1.716 \times {10^{ - 5}}kg/ms$, ${T_{ref}} = 273.15K$ and $Su = 110.4K$. Hence:

where the value of $T$ is obtained from the energy equation and the enthalpy definition, $H = \rho {C}_{p}T$.

The transfer of momentum caused by turbulent fluctuations is based on the Boussinesq eddy-viscosity hypothesis, in which a linear relationship between the turbulent shear stress and the mean strain rate is assumed ^[16]:

where $\tau _{ij}^F$, $\mu {}_{tur}$, ${\tilde S_{ij}}$ and $\tilde K$ represent the Favre-averaged Reynolds stress tensor, eddy viscosity, the Favre-averaged strain rate and the Favre-averaged turbulent kinetic energy, respectively.

The development of RANS equations yields six extra unknowns called Reynolds stresses, for which, the solution requires additional turbulent models. In this work, the one-equation Spalart-Allmaras (SA) turbulence model is used. The SA model, as implemented into the open-source code SU2 ^[17], is presented by ^[18] as follows:

where $\hat \nu$ is the eddy-viscosity variable ^[16].

For the turbulent viscosity:

where the production term is $\hat S = \left| {\vec \omega } \right| + \frac{{\hat \nu }}{{{\kappa ^2}{d^2}}}{f_{\upsilon 2}}$, $\vec \omega = \nabla \times \vec \upsilon$ is the fluid vorticity, $d$ is the distance to the nearest wall, and ${f_{\upsilon 2}} = 1 - \frac{\chi }{{1 + \chi {f_{\upsilon 1}}}}$; ${f_\omega } = g{\left[{\frac{{1 + c_{\omega 3}^6}}{{{g^6} + c_{\omega 3}^6}}} \right]^{{\raise0.7ex\hbox{$1$} \!/ \!\lower0.7ex\hbox{$6$}}}}$, where $g = r + c_{\omega 2}^{}\left({{r^6} - r} \right)$ and $r = \frac{{\hat \nu }}{{\hat S{\kappa ^2}{d^2}}}$. The closure constants are:

2.2.   Extension of the spalart-allmaras model to consider roughness

The technique used to evaluate roughness effects in this paper is based on the equivalent sand-grain roughness, proposed by schlichting ^[7], which is an idealized roughness height that represent the height of a sand grain from Nikuradse's experiments ^[19]. The commonly used correlation proposed by ^[8] to translate topological features into ${k_s}$ is summarized below:

where ${k_s}$ represents the equivalent sand-grain roughness height, $k$ is the mean roughness height, and $\lambda$ is the roughness shape parameter, which is defined as:

where ${r_0}/k$ represents the spacing parameter, ${r_0}$ is the mean distance between roughness elements; ${A_p}$ is the projected area in the direction of the freestream velocity vector; and ${A_s}$ is the windward surface area of the element seen by the flow.

The roughness extension is designed to mimic roughness by increasing the eddy viscosity ${\nu _t}$ at the wall which results in increased skin friction and heat flux. An increased wall distance is also needed to model the roughness induced velocity shift, which is obtained by an imposed offset ${d_{new}} = d + {d_0}$, where ${d_0}$ is a length dependent of ${k_s}$. This offset in linked to the equivalent sand grain roughness based on the velocity profiles presented by Nikuradse ^[19]. Assuming that ${\nu _t} = \hat \nu$ ^[10], spalart reduced the momentum equation to:

from which, the solution reads:

which yields:

The extension of the SA model to consider roughness effects ^[10] has been implemented into an open-source code called Stanford University Unstructured, SU2 ^[17]. For uniform ${k_s}$, the new distance (${d_{new}}$) replaces every $d$ in the original model. The $\chi$ is changed to achieve good predictions for smaller roughness^[10]:

with ${c_{R1}} = 0.5$. As $\hat S$ is not affected, ${f_{\upsilon 2}}$ becomes:

Then, the Dirichlet wall boundary condition (${\hat \nu _{wall}} = 0$) is replaced by:

where $\vec n$ represents the outside normal to the wall vector. This is a Robin boundary condition where the ${\hat \nu _{wall}}$ is obtained at previous iteration by solving Eq 7 with the imposed flux$\frac{{\partial \hat \nu }}{{\partial \vec n}}$.

2.3.   Numerical method

The 3D flows are modeled by the compressible, turbulent Navier-Stokes equations, which are expressed in the conservative form ^[20]. The vector of conservative variables is represented as:

where ρ is the air density, E is the total energy per unit mass, and $\vec \upsilon = {\upsilon _1}, {\upsilon _2}, {\upsilon _3} \in {\mathbb{R}^3}$ is the flow velocity. Then, the Navier-Stokes equations are presented as a general convection diffusion equation:

in which $\vec U$ represents the vector of state variables, ${\vec F^c}(\vec U)$ the convective fluxes, ${\vec F^\upsilon }(\vec U)$ the viscous fluxes and $\vec Q(\vec U)$ a generic source term. Thus, the convective and viscous fluxes are respectively:

where $P$ represents the static pressure, which under the assumption of a perfect gas can be determined from $P = (\gamma - 1)\rho [E- 0.5(\vec \upsilon.\vec \upsilon)]$, $H$ is the fluid enthalpy, ${C_P} = \gamma R/(\gamma - 1)$ is the specific heat of air at constant pressure, $T = P/\rho R$ is the temperature, ${\delta _{ij}}$ is the Kronecker delta function, the viscous stresses are written as:

where ${\Pr _d}$ and ${\Pr _t}$ represent dynamic and turbulent Prandtl numbers. The dynamic viscosity, ${\mu _{dyn}}$, is assumed to satisfy Sutherland's law ^[21], and the turbulent viscosity, ${\mu _{tur}}$, is obtained from the one equation SA turbulent model ^[18]. Moreover, in the framework of the general convection diffusion equation, the roughness SA equation can be rewritten as:

The boundary conditions at farfield is:

and the rough wall boundary conditions are defined as presented in Eq 17.

The system of equations, RANS and turbulence models, are solved iteratively using an implicit time marching approach without multi-grid. For the convective flux interpolation, the Roe scheme is used with a second order MUSCL approach and the Venkatakrishnan limiter. The linear system of equations is solved with a GMRES solver and the Lower-Upper Symmetric Gauss-Seidel (LU-SGS) preconditioner ^[16].

3.   Verification and validation of roughness extension

In this section, the roughness model implementation performed into SU2 code is verified and validated using RANS simulations to estimate lift coefficients and nondimensional velocity profiles of flows over smooth and rough NACA 0012 airfoils. To facilitate the identification of the model active during the simulations throughout this document, the following nomenclature is used: SU2-RANS_smooth is used when RANS simulations are performed with roughness deactivated, while SU2-RANS_rough is used when roughness is activated. SU2-RANS is used for the general case in which both cases are treated. The smooth and rough simulations use exactly the same grids, farfield boundary conditions, and numerical methods. For the rough airfoils, the Robin boundary condition is used at the wall instead of the Dirichlet boundary condition. Also, as the roughness is uniformly distributed on the airfoil, the wall distance is increased proportionally to the sand grain roughness value $d + 0.03{k_s}$.

3.1.   Numerical method smooth and rough NACA 0012 airfoils

First, the behavior of the flow over a smooth NACA 0012 is investigated to confirm whether the implementation has no influence on original model when roughness is not activated. Then, roughness is activated, and the flow over a rough NACA 0012 is investigated. SU2-RANS estimations of lift coefficients and velocity profiles inside the boundary layer are evaluated and compared to numerical ^[22,23,24] as well as experimental ^[25,26,27] results.

For the smooth NACA 0012 airfoil, SU2-RANS_smooth is used to estimate the lift coefficient, and results are compared to estimations of two different codes, CFL3D and OVERFLOW, obtained from ^[22]. The validation is performed using experimental results from ^[25,27]. The parameters are kept as the ones utilized by reference ^[22], being $Ma{\rm{ }} = {\rm{ }}0.15$, $\operatorname{Re} {\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$ and $T{\rm{ }} = {\rm{ }}300K$. The grid used for the smooth airfoil was obtained from NASA's turbulence modeling resource webpage ^[28], and it has 897 × 257 grid points. This grid was selected after a grid studies involving five grids, 113 × 33,225 × 65,449 × 129,897 × 257, and 1793 × 513. The minimum space from the wall to the grid is $y = 8.14 \times {10^{ - 7}}$m, and the average ${y^ + }$ around the airfoil is 0.2. The farfield boundary condition is 500c (c: chord of the airfoil) away from the airfoil to minimize issues associated with the boundaries. This higher than usual values for the farfield boundary condition ensure minimal influence on the solution and enable comparison with reference results ^[28]. However, a value of 100c at AoA 15˚ would have increased the drag coefficient values by 1% and the lift coefficient by less than 0.1%.

For the rough NACA 0012 airfoil, two separated evaluations are performed, similar to what was done by ^[24]. First, SU2-RANS_rough estimations of lift coefficient are compared to numerical results produced by the WMB code used by ^[24], and to experimental data from ^[29]. Mendez uses RANS simulations with two different turbulence models to take roughness into account ^[24], the extension of the k--SST ^[30] and the k- models ^[31]. For the experimental data, Abbott et al. use a 24-inch-chord airfoil, in which 0.011-inch carborundum grains are applied on the first 0.08c of the leading edge of the airfoil ^[29]. The parameters used in our simulations are $\operatorname{Re} {\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$, $Ma{\rm{ }} = {\rm{ }}0.15$ and ${k_s}/c{\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, as the ones used by ^[24].

For the evaluation of the velocity profiles within the boundary layer, the experimental results obtained from ^[26], and the numerical estimations produced by ^[23] are used for comparisons. In the experiments, the measurements were performed using a 0.5334-meter-chord NACA 0012 airfoil with a span of 0.8573 m. The authors simulate roughness using strips tapes of 12.7 × 101.6 mm (0.5 × 4 in). The tapes were 0.35 mm high, and the roughness element's center-to-center distance was kept 1.3 mm. The roughness was spread between $x/c = 0.00612$ ($x/s = 0.0150$, where $s$ is the distance on the upper camber) and $x/c = 0.0258$ ($x/s = 0.388$), which results in a length of 12.7 mm (1/2 in). For the numerical estimations of ^[23], a modified version of the Langtry-Menter $\gamma - \tilde R{e_{\theta t}}$ transition model to account for roughness was used. The model was implemented into OVERFLOW-2.2k, and a ${k_s}/c{\rm{ }} = {\rm{ }}6.56{\rm{ }} \times {\rm{ }}{10^{ - 4}}$ was chosen for their simulations.

Our numerical simulations are performed at a Reynolds number of 1.25 × 10⁶, as used by references ^[24] and ^[23]. The computational domain is the same as the one used for the smooth NACA 0012, and the heat flux is set to zero at the wall. The farfield boundxary is 500c away from the airfoil. Roughness is considered along the whole airfoil. The grid was obtained from NASA's turbulence modeling resource webpage ^[28], and it has 897 × 257 grid points, in which the closest point to the wall is $y = 8.14 \times {10^{ - 7}}$m, and the maximum $y +$ is 0.105. Dirling's correlations ^[8] is used to calculate the equivalent sand-grain roughness height, and the parameter $({A_p}/{A_s})$ was assumed to be 0.5, which resulted in a ${k_s}/c{\rm{ }} = {\rm{ }}6.94{\rm{ }} \times {\rm{ }}{10^{ - 4}}$.

Figures 1 and 2 present the velocity streamlines of flows over smooth and rough NACA0012 airfoils at an AoA of 12 degrees. In both figures, the streamlines are colored by the streamwise velocity, and the range of the legends are based on the maximum and minimum values of the smooth case. When comparing both figures, the first noticeable difference caused by roughness effects is the thickening of the boundary layer, which causes trailing-edge separation as well as flow recirculation, as seen in Figure 2. These effects caused by roughness were also pointed out by references ^[12,32]. In addition, Bragg et al. highlight that these effects are manifested through modified skin friction and pressure distribution into performance degradation ^[32].

Figures 3 and 4 present the variation of the lift coefficient as the AoA of the airfoil is increased. In addition, numerical and experimental results are added to the figures for comparison. As shown in Figure 3, SU2-RANS_smooth, CFL3D and OVERFLOW almost perfectly match, and the maximum discrepancies between the results are inferior to 1%. When compared to experimental results, the maximum error is 2%, obtained with OVERFLOW at an AoA of 15 degrees, and the minimum error presented is 1.18%, achieved by CFL3D. SU2-RANS_smooth presents a good estimation of the maximum lift coefficient, ${C_{l(\max)}} = {\rm{ }}1.66$, at 17 degrees, against ${C_{l(\max)}} = {\rm{ }}1.64$ in the experimental data, which results in an error of 1.2%. For the rough airfoil, Figure 4 shows that SU2-RANS_rough estimations are globally closer to experimental data ^[29] than Mendez's estimations using the WMB code ^[24]. When estimating the ${C_{l(\max)}}$, SU2-RANS_rough and WMB reach the same error of 7%, as compared to experiments.

Figures 5a to 5d depict the profiles of the nondimensional velocities ($\upsilon /{\upsilon _\infty }$) over the rough NACA 0012 (AoA = zero degrees). The cuts are taken in the normal direction on the upper surface of the airfoil. Therefore, the perpendicular-to-the-wall velocity profiles are measured from point 0 (at the wall) up to $7{\rm{ }} \times {\rm{ }}{10^{ - 3}}$m in the normal direction ($y/c$). The black dots represent the experimental reference ^[26], the red lines represent the numerical reference ^[23], and the blue lines are our results. In the first region of the airfoil (Figure 5a), SU2-RANS_rough estimations achieved better agreement with experimental results than the numerical reference. At $x/c{\rm{ }} = {\rm{ }}0.05$, SU2-RANS_rough and the numerical reference agree on estimations for values obtained at a distances superior to $y/c = 1.3 \times {10^{ - 3}}$ from the airfoil. In the region below this distance, SU2-RANS_rough performed better with a maximum error of 10.6% for velocity estimations, against an error of 15.7% of the numerical reference, when considering the same distance $y/c$ of $9.2{\rm{ }} \times {\rm{ }}{10^{ - 4}}$. At $x/c{\rm{ }} = {\rm{ }}0.075$, SU2-RANS_rough and the numerical reference achieved reduced errors globally, as shown in Figure 5b. SU2-RANS_rough results of nondimensional velocities are slightly underestimated with a maximum error of 5.5% at $y/c{\rm{ }} = {\rm{ }}1.4{\rm{ }} \times {\rm{ }}{10^{ - 3}}$, whereas the numerical reference's results are overestimated with a maximum error of 6.6% at the same $y/c$. For the next cuts, the numerical reference estimations are slightly better than SU2-RANS_rough's, as seen in Figures 5c and 5d. Nevertheless, the maximum difference between SU2 and the numerical reference along these cuts is inferior to 11%. The reason for this underestimation is twofold. First, the choice of roughness in the entire airfoil resulted in a somewhat higher thickening of the boundary layer. Second, Langel used a modified version of the Langtry-Menter $\gamma - \tilde R{e_{\theta t}}$ transition model, based on the k-$\mathrm{\omega }$ tubulence model, to account for roughness, which presented a better accuracy in this case ^[23].

Considering the uncertainties involved in the transitional roughness regime, and the fact that the model implemented is not a laminar-turbulent transition model, SU2-RANS_rough predictions are satisfactory when compared to Kerho's experimental data ^[26] and with numerical results from ^[23]. Therefore, these results show that SU2 roughness implementation behaves as expected, and that the code is applicable for the estimation of lift coefficient and velocity profiles over airfoils with rough surfaces.

4.   Maximum lift coefficient predictions

The NACA 0012 and the model 5–6 used for this study are illustrated in Figure 6. The NACA0012 airfoil is a symmetric airfoil often used for model validation in aeronautics. The model 5–6 airfoil is a modern airfoil typical of an outer section of a transport aircraft. The thickness to cord ratio is 12% for the NACA0012 and 10.6% for the model 5–6.

The NACA 0012 airfoil has a chord of 1 m. We use $Ma = 0.15$ and $T = 300K$. The grid used was obtained from NASA's turbulence modeling resource webpage ^[28], and it has 230,529 grid points with the farfield boundary condition at 500c away, as shown in Figure 7. The model 5–6 airfoil has a chord of 1 m. The farfield boundary surface is a circle of diameter 44.39c (Figure 8) defined in a grid with 52,462 vertices and 37,007 cells. The aerodynamic parameters are set as $Ma{\rm{ }} = {\rm{ }}0.2$ and $T{\rm{ }} = {\rm{ }}288.15{\rm{ }}K$. The computational domain and the grid are similar to the ones used in ^[33]. The grid gives almost grid independent results and good predictions of the lift and drag coefficient for several CFD codes, including our modified version of SU2 ^[34].

4.1.   Reynolds number sensitivity analysis

The increase in Reynolds number can result in improvement or degradation of performance. According to ^[12], airfoils with a thickness ratio of < 8% had little improvement in the maximum lift coefficient (${C_l}_{(\max)}$). On the other hand, for geometries with thickness ratios of 9% or higher, a very definite reduction in ${C_l}_{(\max)}$ (from 15 to nearly 40%) was noted. As both airfoils studied herein present thickness ratios higher than 9% (NACA 0012 $t/c{\rm{ }} = {\rm{ }}12\%$, model 5–6 $t/c{\rm{ }} = {\rm{ }}10.6\%$), only performance degradation is expected.

Figure 9 presents estimations of the lift coefficients in function of the AoA for different Reynolds numbers using the NACA 0012 airfoil. The AoA is varied from 11 to 15 degrees, and the Reynolds number is varied from $1{\rm{ }} \times {\rm{ }}{10^6}$ to $18{\rm{ }} \times {\rm{ }}{10^6}$ with increments of $3{\rm{ }} \times {\rm{ }}{10^6}$. The Re reductions are calculated based on the first Re chosen, $\operatorname{Re} {\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^6}$, herein called baseline Re. Hence, $({C_{l(}}_{\max)\_\operatorname{Re} } - {C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}})/{C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}}$, where ${C_{l(}}_{\max)\_\operatorname{Re} }$ is the maximum lift coefficient estimated using a ${k_s}/c = 6 \times {10^{ - 4}}$ at a selected Reynolds number, and ${C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}}$ is maximum lift coefficient with the baseline Re, $\operatorname{Re} {\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^6}$. As a matter of space, $({C_{l(}}_{\max)\_\operatorname{Re} } - {C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}})/{C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}}$is represented as $\Delta {C_{l(}}_{\max)\operatorname{Re} }/{C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}}$ in the graphs. As shown, ${C_l}_{(\max)}$ suffers a decrease of 0.082 units (reduction of 6.7%) when the Reynolds number is increased from $\operatorname{Re} {\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^6}$ to $\operatorname{Re} {\rm{ }} = {\rm{ }}3{\rm{ }} \times {\rm{ }}{10^6}$. In the next increments, the reduction achieves 9.5% at $\operatorname{Re} {\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$, then, 10% at $\operatorname{Re} {\rm{ }} = {\rm{ }}9{\rm{ }} \times {\rm{ }}{10^6}$.

Similar to the previous study case, estimations of the lift coefficient using the model 5–6 airfoil in function of the AoA for different Reynolds numbers are presented in Figure 10. The AoA of the model 5–6 is varied from 10 to 17 degrees, and the Reynolds number from $1{\rm{ }} \times {\rm{ }}{10^6}$ to $18{\rm{ }} \times {\rm{ }}{10^6}$ with increments of $3{\rm{ }} \times {\rm{ }}{10^6}$. As shown in Figure 10, a reduction in ${C_l}_{(\max)}$ of 5% is reached when the Reynolds number is changed from $1{\rm{ }} \times {\rm{ }}{10^6}$ to $3{\rm{ }} \times {\rm{ }}{10^6}$. Then, at $6{\rm{ }} \times {\rm{ }}{10^6}$, a reduction of 6.2% is obtained as compared to the baseline $\operatorname{Re} {\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^6}$. The difference between the $\operatorname{Re} {\rm{ }} = {\rm{ }}3{\rm{ }} \times {\rm{ }}{10^6}$ and $\operatorname{Re} {\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$ is 1.2%, and from$\operatorname{Re} {\rm{ }} = 6{\rm{ }} \times {\rm{ }}{10^6}$ to $\operatorname{Re} {\rm{ }} = {\rm{ }}18{\rm{ }} \times {\rm{ }}{10^6}$, no significant reduction is noticed.

Although practice in aircraft icing favors the ${\mathrm{k}}_{\mathrm{s}}/\mathrm{c}$ parameter, the boundary layer theory relates the friction coefficient to the Reynolds number and the roughness Reynolds number, ${\mathrm{R}\mathrm{e}}_{\mathrm{k}} = {\mathrm{k}}_{\mathrm{s}}{\mathrm{u}}_{\mathrm{\tau }}/\mathrm{\nu }$ with ${\mathrm{u}}_{\mathrm{\tau }}$ the shear velocity. The friction coefficient increases with ${\mathrm{R}\mathrm{e}}_{\mathrm{k}}.$ For the airfoils, even if the ${\mathrm{k}}_{\mathrm{s}}/\mathrm{c}$ is kept constant, ${\mathrm{R}\mathrm{e}}_{\mathrm{k}}$ growth almost linearly with the Reynolds number, as shown in Figure 11. This creates an increase in the friction coefficient, removes more energy from the flow, and triggers an earlier boundary layer separation on the airfoil upper surface.

4.2.   Equivalent sand-grain roughness heights sensitivity analysis

For the evaluation of the decrease in the maximum lift coefficient caused by roughness, a range of equivalent sand-grain roughness heights ${k_s}$, was used. The same two airfoils from the previous analysis were used. Therefore, results of the NACA 0012 airfoil are presented in Figure 12, whereas results of the 5–6 model airfoil are shown in Figure 13.

In Figures 12 and 13, the increase of AoA of the airfoil is represented over the x-axis, and the variation of the lift coefficient is represented over the y-axis. Each curve is a representation of the lift coefficients using a specific equivalent sand-grain roughness height. The reductions are calculated based on the $\left(C_{l(\max )_{-} s m o o t h}-C_{l(\max )}\right) / C_{l(\max )_{-} s m o o t h}$, which is represented as $\Delta C_{l(\max ) k_{s}} / C_{l(\max )_{-} s m o o t h}$ in the graphs.

The first equivalent sand-grain roughness heights chosen for the NACA 0012 airfoil, ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 7}}$, ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 6}}$ and ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 5}}$, cause no significant reduction on $C_{l(\max )}$, as shown in Figure 12. The reason for this lack of roughness effects can be better understood when the roughness Reynolds number ${\operatorname{Re} _k}$, is investigated. Figure 14 illustrates the reduction in lift coefficient as ${k_s}/c$ is increased. In addition, the flow regime based on the ${\operatorname{Re} _k}$ is highlighted in red, green and blue, for aerodynamically smooth, transitional roughness, and fully rough surface, respectively. As seen in Figure 14, the first roughness heights chosen for the NACA 0012 airfoil fall in the aerodynamically smooth regime. According to ^[35], only in the transitional regime that some effects of roughness become noticeable. Among the equivalent sand grain roughness heights chosen for the simulations over the NACA 0012, roughness effects become significant from ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, which marks the beginning of the transitional ${\operatorname{Re} _k}$ regime. Thus, from ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, $C_{l(\max )}$ decreases rapidly, as illustrated in Figure 12. However, as ${k_s}/c{\rm{ }} = {\rm{ }}9{\rm{ }} \times {\rm{ }}{10^{ - 4}}$ is achieved, the changes in the reduction of $C_{l(\max )}$ begins to decrease, as illustrated by the gentler slope of the curve shown in Figure 14. Therefore, even for the increased ${k_s}$ the changes in reduction were mild.

It is also worth discussing at which AoA the $C_{l(\max )}$ occurs over the NACA 0012 airfoil. The AoA of the $C_{l(\max )}$ is the same as that of the smooth case for ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 7}}$, ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 6}}$ and ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 5}}$. Then, from ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 4}}$ even mild increments suffice to drop the AoA of the ${C_{l(}}_{\max)}$. Thus, at ${k_s}/c{\rm{ }} = {\rm{ }}1{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, the AoA of the ${C_{l(}}_{\max)}$ is 16 degrees; at ${k_s}/c{\rm{ }} = 3{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, the AoA drops to 14 degrees; and at ${k_s}/c{\rm{ }} = {\rm{ }}7{\rm{ }} \times {\rm{ }}{10^{ - 4}}$ it becomes 12 degrees. When further increased, the maximum reduction on the AoA achieved is 11 degrees.

For the investigation of the 5–6 model airfoil, a similar behavior is found, as shown in Figure 13. Hence, the maximum lift coefficient suffers no significant influence for small ${k_s}$, and, for sufficiently large ${k_s}$, mild changes in the reduction of $C_{l(\max )}$ is noticed. On the other hand, for intermediate heights of ${k_s}$ (from ${k_s}/c{\rm{ }} = {\rm{ }}2{\rm{ }} \times {\rm{ }}{10^{ - 4}}$ to ${k_s}/c{\rm{ }} = {\rm{ }}16{\rm{ }} \times {\rm{ }}{10^{ - 4}}$), the effect of roughness is elevated, and a significant loss in $C_{l(\max )}$ is seen. This behavior is also illustrated Figure 14, in which the steeper slope of the curve represents the highest influence of roughness.

The AoA of the maximum lift coefficient for ${k_s}/c{\rm{ }} = {\rm{ }}2{\rm{ }} \times {\rm{ }}{10^{ - 6}}$ and ${k_s}/c{\rm{ }} = {\rm{ }}2{\rm{ }} \times {\rm{ }}{10^{ - 5}}$ is the same as that of the smooth airfoil, 17 degrees. Then, for ${k_s}/c{\rm{ }} = {\rm{ }}2{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, ${k_s}/c{\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^{ - 4}}$ and ${k_s}/c{\rm{ }} = {\rm{ }}8{\rm{ }} \times {\rm{ }}{10^{ - 4}}$, a unitary reduction is achieved in each change, 16, 15 and 14 degrees, respectively. From this point on, mild variations in the AoA are seen, and the smallest value reached is 12 degrees for ${k_s}/c{\rm{ }} = {\rm{ }}40{\rm{ }} \times {\rm{ }}{10^{ - 4}}$.

4.3.   Correlation and discussion

As it can be seen in Figure 15 for both the NACA0012 and the model 5–6, a plateau in the reduction is achieved for Re above $6{\rm{ }} \times {\rm{ }}{10^6}$. This behavior results in no evolution on the degradation as Reynolds is increased, which can also be seen in the experimental results of ^[35].

The reduction of the maximum lift coefficient reached by both airfoils can be obtained by Eq 26, which is represented by the read lines shown in Figure 15.

where ${a_1}$, ${b_1}$ and ${c_1}$ are parameters that vary from one airfoil to the other. Table 1 gathers the values of these parameters for the NACA 0012 and the model 5–6 airfoils.

Our results show that although the reduction of the maximum lift coefficient caused by different equivalent sand-grain roughness heights and Reynolds number present similar behavior, they fall into different curve formats. In this way, the reduction of the maximum lift coefficient caused by different equivalent sand-grain roughness heights can be represented by Eq 27. Figure 16 presents the reduction of the maximum lift coefficient reached by both airfoils as ${k_s}$ is increased.

Table 2 gathers the parameters ${a_2}$ and ${b_2}$ for each airfoil.Table 2.

Results show that, for the two airfoils considered, moderate reduction of maximum lift coefficient is noticed at Reynolds numbers greater than $6 \times {10^6}$. This value is in accordance with the observations of ^[12], in which it is shown that Re values much below $5{\rm{ }} \sim {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$ are most likely not meaningful or useful for higher Re applications. Therefore, higher-Reynolds-number simulations could be performed at reduced Reynolds numbers ($\operatorname{Re} {\rm{ }} \geqslant {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$) when estimating lift coefficient curves.

As for the equivalent sand-grain roughness heights investigated, results suggest that variation of sufficiently small heights has no significant influence on the maximum lift coefficient degradation. Moreover, when roughness is continuously increased, a saturation point seems to be approached, in which the variation of the maximum lift coefficient degradation is reduced. On the other hand, in-between values of ${k_s}$ produce the maximum variation of the ${C_l}_{(\max)}$ during the increments.

5.   Conclusions

Using CFD and the Spalart-Allmaras roughness extension developed by Boeing, this work investigated the effects of equivalent sand grain roughness heights on the reduction of the maximum lift coefficient of two airfoils, a NACA 0012 and a model 5–6. Our results demonstrated that when roughness is continuously increased, a saturation point seems to be approached, in which the variation of the maximum lift coefficient degradation is reduced. As a Reynolds number investigation shows, lift degradation reach a plateau for values greater than $\operatorname{Re} {\rm{ }} = {\rm{ }}6{\rm{ }} \times {\rm{ }}{10^6}$. In this way, the maximum lift coefficient reduction caused by the equivalent sand-grain roughness height can be successfully represented by $\Delta {C_{l(\max ){k_s}}}/{C_{l(\max )\_smooth}} = \left( {{a_2}{k_s}} \right)/\left( {{k_s} + {b_2}} \right)$, whereas the reduction caused by Re is represented by $\Delta {C_{l(}}_{\max)\operatorname{Re} }/{C_l}_{(\max)\_\operatorname{Re} = 1 \times {{10}^6}}{\rm{ }} = {\rm{ }}{a_1}{\operatorname{Re} ^{{b_1}}} + {c_1}$.

We believe that those correlations can be useful since scaled-down airfoils are commonly used to evaluate ice formation in wind-tunnel tests. Another concern is the computational cost of CFD simulations, which is increased at higher Reynolds numbers. Therefore, lower Re simulations are desired for industrial applications.

Future research should consider hybrid models such as DES and DDES with the roughness extension implemented to improve the prediction of performance degradation due to icing. These models allow simulations of study cases with massive flow separation, for which RANS is not recommended. This will enable the study of post-stall aerodynamic performance.

Acknowledgment (All sources of funding of the study must be disclosed)

Computations were made on supercomputers managed by Calcul Québec and Compute Canada. The operation of these supercomputers is funded by the Canada Foundation for Innovation (CFI), Ministère de l'Économie, des Sciences et de l'Innovation du Québec (MESI) and le Fonds de recherche du Québec - Nature et technologies (FRQ-NT).

Some numerical simulations presented in this paper were carried out using the PLaFRIM experimental testbed. PLaFRIM cluster being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS, see https://www.plafrim.fr/.

Gitsuzo B.S. Tagawa would like to thank the support of CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico—Brazil).

Conflict of interest

The authors declare no conflict of interest.