On the Fidelity of RANS-Based Turbulence Models in Modeling the Laminar Separation Bubble and Ice-Induced Separation Bubble at Low Reynolds Numbers on Unmanned Aerial Vehicle Airfoil

Abstract

The operational regime of Unmanned Aerial Vehicles (UAVs) is distinguished by the dominance of laminar flow and the flow field is characterized by the appearance of Laminar Separation Bubbles (LSBs). Ice accretion on the leading side of the airfoil leads to the formation of an Ice-induced Separation Bubble (ISB). These separation bubbles have a considerable influence on the pressure, heat flux, and shear stress distribution on the surface of airfoils and can affect the prediction of aerodynamic coefficients. Therefore, it is necessary to capture these separation bubbles in the numerical simulations. Previous studies have shown that these bubbles can be modeled successfully using the Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS) but are computationally costly. Also, for numerical modeling of ice accretion, the flow field needs to be recomputed at specific intervals, thus making LES and DNS unsuitable for ice accretion simulations. Thus, it is necessary to come up with a Reynolds-Averaged Navier–Stokes (RANS) equation-based model that can predict the LSBs and ISBs as accurately as possible. Numerical studies were performed to assess the fidelity of various RANS turbulence models in predicting LSBs and ISBs. The findings are compared with the experimental and LES data available in the literature. The structure of these bubbles is only studied from a pressure coefficient perspective, so an attempt is made in these studies to explain it using the skin friction coefficient distribution. The results indicate the importance of the use of transition-based models when dealing with low-Reynolds-number applications that involve LSB. ISB can be predicted by conventional RANS models but are subjected to high levels of uncertainty. Possible recommendations were made with respect to turbulence models when dealing with flows involving LSBs and ISBs, especially for ice accretion simulations.

1. Introduction

Turbulent flow often dominates the operational regime of manned aircraft traveling at high Reynolds numbers ($Re$) of $O\left({10}^6\right)$ and higher. Under such circumstances, viscous force causes flow separation to be delayed. The flow field around a UAV flying at a low $Re$ of ${O(10}^5$) or lower is primarily of laminar nature and is sensitive to even minor pressure differences. Previous studies have demonstrated the formation of LSBs at these Reynolds numbers [1]. An extensive review of the flows involving LSBs was performed by Tani [2] in 1964. Ice accumulation on the leading side of an airfoil can lead to the creation of an ISB. Bragg et al. [3] primarily studied the production of such bubbles and the related deterioration of aerodynamic performance. Most of his studies are based on manned aircraft Ice accretion with Reynolds numbers in the order of ${10}^6$ and the observations were not verified for low Reynolds number ranges. Manaf [4,5] conducted a comprehensive review of the impact of ice accumulation on UAVs and their aerodynamics in 2022–2023. The appearance of LSBs and ISBs can affect the aerodynamic coefficients by influencing the pressure and shear stress distributions considerably [6]. Further, it can also affect the heat flux distribution of the airfoil. The numerical ice accretion solvers demand accurate modeling of the flow field around the airfoil because heat flux and shear stress are the two main parameters that determine the water transport and ice accretion phenomena on the surface of the airfoil [7]. The capability of DNS and LES in predicting LSBs and ISBs can be found in the literature [1,8]. Implementing DNS/LES models for ice accretion simulations is expensive and is not feasible for multi-shot simulations (an iterative process to predict ice accretion with recalculation of the flow field at specific intervals). Loth [9] estimated the computational cost of various numerical approaches on the basis of the Reynolds number and found values of $R{e}_L^{2.3}, R{e}_L^{1.8},\mathrm{a}\mathrm{n}\mathrm{d} R{e}_L^{0.2}$ for DNS, LES, and RANS simulations, respectively. Therefore, using DNS and LES for an ice accretion prediction numerical solver is not feasible, especially for multi-shot icing simulations, which demand a flow field recalculation at every shot. Thus, it is necessary to come up with an effective RANS-based model that can predict LSBs and ISBs as accurately as possible. This study attempts to compare nine distinct RANS-based turbulence models to assess their capability in predicting LSBs and ISBs on Rg-15 airfoil at low values of $Re$. The effectiveness of existing turbulence models based on RANS equations in predicting LSBs and ISBs is studied by comparing their predictions with data available in the literature. Oo et al. [6] performed experimental and LES studies on LSBs and ISBs in 2020 around an Rg-15 airfoil at $Re=O\left({10}^5\right)$, and the separation and reattachment locations were identified. Further, an attempt was also made to understand the dynamics and structure of these bubbles from the perspective of the skin friction coefficient, which was studied earlier but only from the perspective of pressure coefficients. A brief literature review of the current wisdom about LSBs and ISBs is presented below.

1.1. Laminar Separation Bubble (LSB)

During the transition of flow from laminar to turbulent flow, the surface of the laminar boundary layer that was previously attached begins to detach from the surface because of unfavorable pressure gradients. It then reattaches itself at a certain downstream location, thus resulting in the creation of LSB [2]. These separation bubbles have the potential to modify the flow field and the effective geometry of the airfoil, where they are more likely to occur on the suction surface. Even slight modifications in the geometry and operational parameters can significantly impair the performance of low $Re$ airfoils. Those separation bubbles are more likely to be observed on the suction side of the airfoil with the potential to modify both the effective curvature of the airfoil and the flow field’s behavior. Airfoils designed for low Reynolds numbers are highly sensitive, and even minor modifications to their geometry and operating conditions can result in a mortification of their performance. LSB significantly influences the distribution of pressure, shear stress, and heat flux, thereby influencing the moment values, lift, drag, and nature of stalling of airfoils [10,11]. Jones [11] documented the presence of LSBs in 1934 as part of his experimental inquiries into the stalling characteristics of aircraft wings. Subsequent researchers conducted numerous experimental investigations in an effort to comprehend the behavior and form of LSB; a comprehensive compilation of these research studies is available in reference [12]. The research conducted by Horton [13] in 1968 yielded valuable insights into the physical mechanisms underlying the expansion and rupture of the bubble. As shown in Figure 1, he also postulated a fundamental structure for the LSB. At point $S$, the laminar boundary layer becomes detached from the surface of the airfoil. At point $T^{\prime },$ the flow transition from laminar to turbulent occurs, which takes place at the maximal height of the separation bubble. Eventually, turbulent mingling eliminates the reverse flow in the vicinity of the wall, and at point $R$, the flow reattaches. These are referred to as ‘short laminar separation bubbles’, and their impact on the pressure distribution is confined to a local area. A ‘short’ bubble can, however, expand and ‘burst’ into a ‘long laminar separation bubble’ in response to an increase in the angle of incidence or a velocity decrease. Long bubbles may persist as an ‘unattached free shear layer’ or reattach at a considerable downstream location, both of which significantly impact the distribution of pressure on a global scale [2]. The airfoil’s stalling behavior is determined by the nature of its bursting.

1.2. Ice-Induced Separation Bubble (ISB)

The adverse effect of ice accretion on the aerodynamic performance of an airfoil was initially reported in the studies of Gray [15] on NACA airfoils in 1953. Ice accretion on the leading side of the airfoil leads to significant performance losses like a decrease in lift and an increase in drag and is widely studied in the literature [16,17,18,19]. These studies were limited to the measurement of lift, drag, and moment on the iced airfoils only, and they lack information about the flow over iced airfoils. Later, from 1953 onwards, Bragg et al. [20,21,22,23,24,25] conducted a sequence of studies to understand the nature of the flow field in the vicinity of an iced airfoil. Measurement of the surface pressure was performed to understand the regions of flow separation, focusing on the location of separation, the reattachment point, and the length of the bubble [22,25]. The presence of an ISB was visualized for the first time by Bragg et al. [20] in 1984 using the oil flow visualization technique. In 1985, he [26] measured the time-averaged vertical velocity profiles inside the glaze ISB using split hot film measurements. The stagnation streamlines and dividing streamlines are traced using the measured velocity profiles, thereby providing the profile of ISBs at different angles of attack. In 1986, he [21] extended the studies using densely distributed surface pressure probes to study ISBs in detail. In the bubble zone, a nearly constant but slightly dropping pressure is followed by a pressure recovery region. The region of constant pressure can be considered the separated flow zone and the flow reattachment is expected to occur in this pressure recovery region. In 1992, he [24] concluded that the behavior of ISBs is quite similar to that of the LSBs in that the boundary layer transitions into turbulent flow with reattachment. Gurbacki [27] observed in 2002 that the reattachment location of the ISB varies as a function of time due to the inherent unsteadiness in the flow, leading to a reattachment zone, as shown in Figure 2. The length of the reattachment zone increases with an increase in the angle of attack, attributed to higher levels of flow unsteadiness.

Numerical studies about the aerodynamics of iced airfoils were extensively reviewed by Stebbins [29] in 2019 and are not intended to be repeated in detail here. Most of the numerical studies were performed on wings with leading-edge horn Ice accretions [30,31,32,33,34,35,36,37,38,39,40] and spanwise ridge ice [41,42,43]. Baldwin–Lomax (BL) and the Spalart–Allmaras (SA) model for turbulence were the widely employed turbulence models in studies with horn ice. The separation extent is overpredicted by RANS models, especially at angles of attack equal to and above 6 degrees, indicating that the accuracy of prediction decreases with an increase in the length of the separation bubble [29]. The numerical simulation predicted very high local suction peaks in the $Cp$ distribution, whereas the experimental results are smooth. Further, there were discrepancies in the prediction of surface pressure recovery on the pressure side. For the spanwise ridge ice, all the RANS-based studies were performed using the SA model. The distribution of pressure on the pressure side is accurately predicted, whereas the predictions of the same on the suction side were poor. For both ice shapes, RANS models predict the aerodynamic coefficients with reasonable accuracy in the linear-lift region but fail at higher angles of attack due to massive flow separation. Most of these studies were performed at Reynolds numbers in the order of ${10}^6$ and Mach numbers above 0.1. However, the operating Reynolds numbers of UAVs range from the higher end of ${10}^4$ to the lower end of ${10}^6$, and the flow field is predominantly laminar at these Reynolds number ranges. Not many studies were performed at these Reynolds number ranges to study the behavior of ISBs. The existing studies were limited to SA and BL turbulence models, and the capability of other RANS-based models has not been studied. The behavior of the ISB at low-Reynolds-number flows and the need for transition-based models in these flow regimes also need to be evaluated.

1.3. RANS-Based Turbulence Model for Prediction of LSBs and ISBs

The laminar-to-turbulent transition of the boundary layer has a significant impact on the flow characteristics. The distribution of shear stress and heat flux on the wall can also be thought of as a function of the location of turbulence onset and the length of transition regions. The transition process can also affect boundary layer separation behavior and the aerodynamic behavior of airfoils [44]. A major limitation of the existing numerical models is their inability to incorporate all those transition effects into a single model. Another issue stems from the fact that the transition process involves both linear as well as nonlinear effects, but the averaging technique used in Reynolds Average Navier Stokes equations (RANS) eliminates the effects of non-linear disturbance growth. Further, most transition models are based on non-local formulations as they use information about the flow variables outside the boundary layer and the integral thickness of the boundary layer. The dependence of transition models on non-local flow variables precluded its compatibility with existing CFD models. Modern CFD algorithms are incapable of computing integral boundary layer parameters due to their formulation based on unstructured grid and domain decomposition methods. Even though stability-based methods like the $e^n$ method [45,46] can eliminate these limitations, the requirement of prior knowledge of the grid and geometry also limits their application.

High-Reynolds-number turbulence models can often be used to predict the transition using two unique approaches. The incorporation of additional low-Reynolds-number model approximations is one option for managing turbulence formation during transitions. In order to address the implications of diminished turbulence intensity within a viscous sublayer, the equations incorporate damping functions. Despite the development of numerous low$-Re$-turbulence models throughout the years [47,48,49,50], their capacity to predict the transition phenomenon appears to be coincidental and, hence, untrustworthy [44]. The second approach is to model the transition process as a mixture of laminar and turbulent flow in proportions that vary with the intermittency factor ($\gamma$). The value of $\gamma$ ranges from 0 at the start of the transition to 1 at the end, and the overall flow could be $\gamma$ times turbulent and ($1-\gamma$) times laminar [51]. This method necessitates empirical information regarding the position of transition initiation and the exact streamwise variation of $\gamma$. These correlations relate the transition Reynolds number based on momentum thickness ($R{e}_{\theta t}$) with local variables like the turbulence intensity $\left(T_u\right)$ and the pressure gradient $\left({\lambda }_{\theta }\right)$. Dhawan [52] in 1958 made such correlations as a function of the momentum thickness $\left(\theta \right)$. In 1980, Abu-Ghannam and Shaw [53] proposed the empirical correlations for the prediction of the start and end of the transition by conducting experiments at different values of turbulence intensity $\left(T_u\right)$ and pressure gradient $\left(\lambda \right)$. Mayle [54] proposed correlations for zero-pressure-gradient flows in 1991, which are functions of $T_u$ and $\theta$. In 1994, Gostelow et al. [55] proposed an empirical correlation for the streamwise distribution of intermittency, based on $\lambda$, $T_u$, and $\theta$. Steelant and Dick [56] proposed a transport equation for intermittency ($\gamma$) in 1996 and used it in conjunction with a conditioned Navier Stokes equation. The transition criteria are often associated with a two-equation turbulence model by an adjustment of the turbulent production term using an intermittency equation. The intermittency transport equation was incorporated into the $k-ϵ$ turbulence model by Cho and Chung [57] in 1992. In 2000, Suzen [58] coupled his intermittency transport equation with the SST $k-\omega$ turbulence model. Most of these transition models are based on non-local formulations as they use information about the flow variables outside the boundary layer and the integral thickness of the boundary layer. To address these issues, Menter [44] proposed basic formulations of transition models based on local variables in 2002. In 2004, Menter [59] added a transport equation for the $R{e}_{\theta t}$ to the SST $k-\omega$ turbulence model and called it the $\gamma -R{e}_{\theta t}$ transition model. He introduced a local property called the vorticity Reynolds number ($R{e}_V$). The value of the maximum vorticity Reynolds number in the profile of a boundary layer is directly linked to the momentum thickness Reynolds number. To induce the transition, the vorticity Reynolds number is used rather than the momentum thickness Reynolds number directly.

2. Numerical Methodology

In this study, the 2D-incompressible RANS equations were solved with a finite volume formulation called SIMPLE for the pressure–velocity coupling. For the spatial discretization of fluxes, a second-order accurate upwind scheme is used. The various turbulence models employed in this study are briefly discussed below.

2.1. Turbulence Model

The Spalart–Allmaras (SA) model [60] is essentially a low-Reynolds-number model that uses a one-equation transport model for solving kinematic eddy viscosity $\left(\tilde{v}\right)$. The transport equation for the SA model is given below in Equation (1).

$$\frac{\partial }{\partial t}\left(\rho \tilde{v}\right)+\frac{\partial }{\partial x_i}\left(\rho \tilde{v}{u}_i\right)=G_v+\frac{1}{{\sigma }_v}\left[\frac{\partial }{\partial x_j}\left\{\left(\mu +\rho \tilde{v}\right)\left(\frac{\partial \tilde{v}}{\partial x_j}\right)\right\}+C_{b2}\rho {\left(\frac{\partial \tilde{v}}{\partial x_j}\right)}^2\right]-Y_v+S_{\tilde{v}}$$

(1)

where $\mu$ is the turbulent viscosity and $G_v$ and $Y_v$ are the turbulent viscosity production and destruction terms.

The turbulent kinetic energy $\left(k\right)$ and the specific dissipation rate $\left(\omega \right)$ transport equations are solved in the standard $k-\omega$ turbulence model as given in Equations (2) and (3). Because this model is sensitive to freestream conditions, an SST $k-\omega$ turbulence model was proposed in 1994 by Menter [61]. This model combines the advantages of the $k-\epsilon$ model [62] in the far field with that of the $k-\omega$ model in the region near the wall [63]:

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left({\mathsf{\Gamma }}_k\frac{\partial k}{\partial x_j}\right)+G_k-Y_k+S_k+G_b$$

(2)

$$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_i}\left(\rho \omega u_i\right)=\frac{\partial }{\partial x_j}\left({\mathsf{\Gamma }}_{\omega }\frac{\partial \omega }{\partial x_j}\right)+G_{\omega }-Y_{\omega }+S_{\omega }+G_{\omega b}$$

(3)

where $G_k$ and $G_{\omega }$ represent the generation of $k$ and $\omega$, respectively, due to the mean velocity gradients. ${\mathsf{\Gamma }}_k$ and ${\mathsf{\Gamma }}_{\omega }$ represent the effective diffusivity of $k$ and $\omega$, respectively. $Y_k$ and $Y_{\omega }$ represent the dissipation of $k$ and $\omega$ due to turbulence. $S_k$ and $S_{\omega }$ are user-defined source terms. $G_b$ and $G_{\omega b}$ account for buoyancy terms.

The interaction of the transition model with the SST turbulence model occurs by modifying the terms $G_k$ and $Y_k$ in the $k$-equation as shown in Equations (4) and (5), respectively.

$$G_k^{*}={\gamma }_{eff}{\tilde{G}}_k$$

(4)

$$Y_k^{*}=\min \left(\max \left({\gamma }_{eff},0.1\right),1.0\right)$$

(5)

where the original production and destruction terms for the SST model are represented by ${\tilde{G}}_k$ and $Y_k$. ${\gamma }_{eff}$ physically represents the percentage of time the turbulent fluctuations are present in the boundary layer, and it varies between 0 and 1 by damping the production of turbulence in the laminar and transitional regions of the boundary layer. The damping limiter ensures that the dissipation rate does not drop below 10% of the turbulent value. The terms $G_k$ and $Y_K$ are modified as shown in Equations (6)–(8), such that:

$$G_k^{*}=0 \mathrm{and} Y_k^{*}=0.1 Y_k \mathrm{in} \mathrm{the} \mathrm{laminar} \mathrm{boundary} \mathrm{layer}$$

(6)

$$0<G_k^{*}<1 \mathrm{and} Y_k^{*}={\gamma }_{eff} Y_k (0.1<{ \gamma }_{eff}<1) \mathrm{in} \mathrm{the} \mathrm{transitional} \mathrm{boundary} \mathrm{layer}$$

(7)

$$G_k^{*}=1 \mathrm{and} Y_k^{*}=Y_k \mathrm{in} \mathrm{the} \mathrm{turbulent} \mathrm{boundary} \mathrm{layer}$$

(8)

The algebraic transition model formulates the intermittency as an algebraic equation instead of the transport equation for intermittency [63]. The algebraic equation is integrated into the $k-\omega$ based turbulence models by multiplying intermittency ($\gamma$) into the production term $p_k$. The empirical correlation is made in such a way that it fits the experimental transition locations.

$$\gamma =f(R{e}_{\theta C}, R{e}_v, Tu, {\lambda }_{\theta })$$

(9)

where $R{e}_{\theta ,c}$ is the critical Reynolds number where the intermittency first starts to increase in the boundary layer, $Tu$ is the turbulence intensity, and ${\lambda }_{\theta }$ is the non-dimensional pressure gradient parameter. $R{e}_v$ is the vorticity Reynolds number defined as a function of the wall distance ($d_w$) and strain rate magnitude (S).

$$R{e}_v=\frac{\rho d_w^{2}S}{\mu }$$

(10)

The intermittency ($\gamma$) transition model [64] solves one transport equation for intermittency (see Equation (11)), where the production term controls the length of transition and the dissipation term allows the boundary layer to re-laminarize by dissipating the intermittency fluctuations. The momentum thickness Reynolds number is algebraically computed using local variables [64].

$$\frac{\partial }{\partial t}\left(\rho \gamma \right)+\frac{\partial }{\partial x_j}\left(\rho u_j\gamma \right)=\frac{\partial }{\partial x_j}[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }}\right)\frac{\partial \gamma }{\partial x_j} ] +P_{\gamma }-E_{\gamma }$$

(11)

The intermittency production term $P_{\gamma 1}$ is given by Equation (12):

$$P_{\gamma 1}=F_{length} \rho S\gamma \left(1-\gamma \right)F_{onset}$$

(12)

where $F_{length}$ is the empirical correlation that controls the length of the transition region and $F_{onset}$ switches on the production of intermittency.

$$F_{onset}=f\left(R{e}_{\theta c}, R{e}_v,d_w, \rho , k,\frac{1}{\omega },\frac{1}{\mu },s, Tu, {\lambda }_{\theta }\right)$$

(13)

$$F_{length}=100$$

(14)

In order to account for the transition behavior, the transition $k\omega -\mathrm{S}\mathrm{S}\mathrm{T}$ turbulence model, also called the $\gamma -R{e}_{\theta } transition model$ [65], is coupled to the $\mathrm{S}\mathrm{S}\mathrm{T} k-\omega$ model with a transport equation for intermittency and the transition Reynolds number. The former is given in Equation (15).

$$\frac{\partial }{\partial t}\left(\rho \gamma \right)+\frac{\partial }{\partial x_j}\left(\rho u_j\gamma \right)=\frac{\partial }{\partial x_j}[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }}\right)\frac{\partial \gamma }{\partial x_j} ]+P_{\gamma 1}-E_{\gamma 1}+P_{\gamma 2}-E_{\gamma 2}$$

(15)

The intermittency production term $P_{\gamma 1}$ is given by:

$$P_{\gamma 1}=F_{length}C{a}_1\rho S{\left(\gamma F_{onset}\right)}^{C_{\gamma 3}}$$

(16)

$$E_{\gamma 1}=P_{\gamma 1}\gamma$$

(17)

Similarly, $F_{onset}$ switches on the production of intermittency and $F_{length}$ controls the length of the length of the transition region.

$$F_{onset}=f\left(R{e}_{\theta c}, R{e}_v,d_w, \rho , k,\frac{1}{\omega },\frac{1}{\mu },s, Tu, {\lambda }_{\theta }\right)$$

(18)

$$F_{length}=f\left({\overline{Re}}_{\theta t}\right)$$

(19)

$$R{e}_{\theta C}=f\left({\overline{Re}}_{\theta t}\right)$$

(20)

$$R{e}_{\theta t}=f(Tu,\lambda )$$

(21)

$R{e}_{\theta ,t}$ is the transition Reynolds number where the transition occurs, and a transport equation is solved for ${\overline{Re}}_{\theta t}$

$$\frac{\partial }{\partial t}\left(\rho {\overline{Re}}_{\theta t}\right)+\frac{\partial }{\partial x_j}\left(\rho U_j{\overline{Re}}_{\theta t}\right)=\frac{\partial }{\partial x_j}[{\sigma }_{\theta t}\left(\mu +{\mu }_t\right)\frac{\partial {\overline{Re}}_{\theta t}}{\partial x_j} ]+P_{\theta t}$$

(22)

Transport equations for turbulent kinetic energy ${(k}_T)$, laminar kinetic energy ($k_L$), and the inverse turbulent time scale ($\omega$) is solved in the $k-kl-\omega$ transition model [66], which is classified as a three-equation eddy-viscosity type.

$$\frac{D{k}_T}{Dt}=P_{k_T}+R+R_{NAT}-\omega k_T-D_T+\frac{\partial }{\partial x_j}\left[\left(\upsilon +\frac{{\alpha }_T}{{\alpha }_k}\right)\frac{\partial k_T}{\partial x_j}\right]$$

(23)

$$\frac{D{k}_L}{Dt}=P_{k_L}-R-R_{NAT}-D_L+\frac{\partial }{\partial x_j}\left[\upsilon \frac{\partial k_T}{\partial x_j}\right]$$

(24)

$$\begin{array}{cc}\frac{D_{\omega }}{D_t}=C_{\omega 1}\frac{\omega }{k_T}{P}_{k_T}& +\left(\frac{C_{\omega R}}{f_W}-1\right)\frac{\omega }{k_T}\left(R+R_{NAT}\right)-C_{\omega 2}{\omega }^2+C_{\omega 3}{f}_{\omega }{\alpha }_T{f}_W^2\frac{\sqrt{\left(k_T\right)}}{d^3}\hfill \\ & +\frac{\partial }{\partial x_j}\left[\left(\upsilon +\frac{{\alpha }_T}{{\alpha }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right]\hfill \end{array}$$

(25)

The effective diffusivity term in Equations (2) and (3) depends on the coefficient ${\alpha }^{*}$. A low $Re$ correction damps the turbulent viscosity by modifying the coefficient ${\alpha }^{*}$.

$${\alpha }^{*}={\alpha }_{\infty }^{*}\left(\frac{{\alpha }_0^{*}+\frac{R{e}_t}{R_k}}{1+\frac{R{e}_t}{R_k}}\right)$$

(26)

where $R{e}_t=\frac{\rho k}{\mu \omega }$ and the remaining terms are constants.

The author intended to provide a very concise explanation of the turbulence models considered in this study, and many terms are not defined in detail. Reference [63] provides a full discussion of the numerous source terms, empirical relations, and constants employed in these turbulence models. Since all these turbulence models are based on many empirical equations and constants, researchers have shown that the finetuning of such constants can improve the predictions. However, it will remain case-specific and cannot be considered universal for industrial applications. Therefore, the numerical simulations in this research have been performed using the standard empirical constant recommended by Ansys [63] to find the best generic turbulence model for low-Reynolds-number studies.

2.2. Computational Model

Computational simulations were conducted on the low-Reynolds-number airfoil Rg-15, which is well-suited for unmanned aerial vehicle (UAV) applications. The reference airfoil had a chord length of $C=210 \mathrm{m}\mathrm{m}$, resulting in $1.07\times {10}^5$ Reynolds number at $7.25 \mathrm{m}/\mathrm{s}$ velocity. The ice shape was adopted from the studies of Williams [67], and the experimental data about LSBs and ISBs are derived from the research of Oo et al. [6]. Numerical simulations were run using the Ansys FLUENT program v18.2, a widely used commercial tool. Multiblock structured 2-D C-grid was made using the software called ICEM CFD, following the guidelines outlined in reference [63] for transition models. The flow domain stretches from the leading edge in all directions for ten times the chord length. To guarantee that the Y-plus value stays below 1, the first cell with a height of $1\times {10}^{-6}$ is used, and the computation mesh is massed in the streamwise directions. Grid independence studies are conducted to verify that the separation bubble is unaffected by changes in the grid resolution, as shown in Figure A1 of Appendix A. The computation utilizes a numerical grid consisting of 231,032 cells (371 grid points on the airfoil surface), as depicted in Figure 3. The simulations were conducted for a significant number of iterations to guarantee the convergence of force coefficients, precision in LSB location, and the fall of residuals to a value significantly below $1\times {10}^{-6}$. The research was conducted for $0, 3, \mathrm{a}\mathrm{n}\mathrm{d} 6^{\mathrm{o}}$ angles of attack, and its findings match the literature data [6]. Turbulence intensity can be considered a crucial aspect influencing the LSB. A wind tunnel with a turbulence level of 1.2% is used for experimental testing and, therefore, the same value is also employed for numerical research. A comparison of the predictions of LSBs using steady-state and time-averaged simulations is shown in Figure A2 of Appendix A. Detailed comparisons can also be found in the earlier research works of the author [68]. At lower angles of attack, the differences in the predictions of steady-state and time-averaged simulations are negligible. Therefore, considering the computational cost, steady-state simulations are performed here. The experiments were performed in the literature on a 2D infinite wing, and therefore, the simulations are also performed in 2D for one-to-one comparisons. The above-mentioned simplifications are adopted since the goal of this study is to assess the capacity of different models to predict LSBs and ISBs. However, from a physics perspective, these separations are three-dimensional and unsteady.

3. Results and Discussion

3.1. Laminar Separation Bubbles (LSBs)

Figure 4 depicts the traits of LSBs as the angle of attack on the RG-15 airfoil surface increases. At an angle of attack of zero degrees, the boundary layer detaches from the surface close to the airfoil trailing edge and does not reconnect to the surface of the airfoil, thus resulting in an open-type LSB. An upstream movement of the separation onset is observed with an increase in the angle of attack to 3 degrees. Also, the flow reattaches to the airfoil surface, resulting in a closed-type LSB. The separation onset travels upstream with a decrease in length and thickness when the angle of attack changes to six degrees. In experiments, the presence of an LSB is detected utilizing the $C_p$ distribution curve [24,69] as shown in Figure 5. Due to the same pressure within the separated flow regions, an LSB is denoted by a flattening of the $C_p$ distribution curve. The initiation of separation is denoted by the beginning of $C_p$ flattening, and the termination of flattening is designated as the transition point. The instant at which the $C_p$ value returns to its inviscid state is considered the juncture of reattachment.

The value of coefficient of skin friction ($C_f$) is negative inside the separation bubble due to the reverse flow, and it can be used to distinguish the regions of flow separation [70]. However, the $C_f$ curve is not examined in detailed like the $C_p$ distribution curve. Thus, an attempt is made here to study the geometry of LSBs from the $C_f$ curve. Figure 6 shows the LSB, $C_p$ distribution, and $C_f$ distribution over Rg-15 at a 6-degree angle of attack predicted by the $\gamma -R{e}_{\theta }$ transition model. The $\gamma -R{e}_{\theta }$ transition model is chosen because of its closer agreement with experimental results as discussed in Section 3.2. The point where the $C_f$ curve initially crosses zero is considered the point of separation onset and, therefore, the location denoted by ‘$S$’ in Figure 6 could be considered the location of separation onset. After the separation onset, the thickness of the bubble gradually increases and the magnitude of negative $C_f$ also increases up to some point before it remains flat until point ‘$T$’ (X/C = 0.378). $C_f$ is a function of the gradient of velocity in the direction normal to the wall, and an increase or decrease in the velocity gradient is responsible for the sudden variations in $C_f$, as shown in Figure 6. Compared to a laminar boundary layer, velocity increases more rapidly in the direction normal to the wall within a turbulent boundary layer. Consequently, the profile of the turbulent boundary layer is exponential, whereas the profile of the laminar boundary layer is parabolic instead. In laminar flow, the free-flowing fluid’s kinetic energy is transferred to the slower-moving fluid at the surface solely by viscosity, i.e., frictional shear stresses. There is no entrainment of energy from the mean flow into the boundary layer observed in the laminar region of the separation bubble [71]. Therefore, the fluid inside the laminar part of the separation bubble is recirculating at a reduced velocity. In a turbulent flow, the kinetic energy of the mean free stream is entrained into the boundary layer and, therefore, the fluid gains more energy and its velocity increases. The velocity profiles along the airfoil surface at specific chordwise locations are shown in Figure 7. The velocity magnitudes are normalized by the local velocity of the main flow, and the chordwise positions are normalized by chord length. For the positions upstream of the maximum thickness point in LSBs (T, X/C = 0.378), the reverse flow velocity profile is parabolic in nature, and after point ‘T’, it becomes exponential. Thus, during the onset of LSB, there is a small velocity gradient due to the development of reverse flow, and because the velocity gradient is smooth in a parabolic profile, the $C_f$ remains relatively flat. At point ‘$T$’, the LSB has attained the maximum thickness, and it can be considered the start of the transition process. The laminar-to-turbulent transition can be considered as a process that occurs over a finite distance and, therefore, it cannot be defined at a specific point [70]. The “Transition zone” could be considered a better expression than the transition location or transition point when speaking about flow transitions. The sudden fall in the $C_f$ value (sudden increase in the negative $C_f$ value) is an indication of a full transition and is denoted by the letter ‘$T^{\prime }$’ in Figure 6. From Figure 7, it can be observed that the profile of the reverse flow starts to become exponential at point ‘T’ (X/C = 0.378), and the gradient becomes sharper at point ‘$T^{\prime }$’ (X/C = 0.393). Therefore, the region ‘${T-T}^{\prime }$’ can be considered the transition zone. The entrainment of energy from the mean flow develops a greater velocity gradient in the turbulent region of the bubble, and the $C_f$ value increases rapidly up to point ‘$R^{\prime }$’. After this point, the velocity gradient decreases until the flow reattaches at Point ‘$R$’.

The slope of the ‘$R-R^{\prime }$’ curve depends on the rate of entrainment of energy from the mean flow. The entrained fluid rolls up into vortices and is then propagated downstream by shedding vortices, so no apparent reverse flow layer is observed [71]. This entrainment, combined with the favorable pressure gradient and the shear layer’s dominance over the region of “dead fluid”, allows the boundary layer to quickly reattach and form a bubble. The fluid flow appears to be very smooth and orderly after reattachment, with vortex structures stretching as they are connected downstream [72]. It could be noted from the $C_p$ distribution curve that the region inside the bubble is not completely flat and is increasing at a slower rate. Also, it can be observed that, at the beginning of the LSB, the $C_p$ distribution is not affected much, and the effect starts to be revealed after the LSB has gained some thickness. Therefore, using the $C_p$ distribution curve to measure the separation and reattachment point has its limitations. In the $C_f$ graphs presented in this manuscript on Section 3.2 and Section 3.4, the black dashed horizontal line represents the $C_f=0$ line and the black dash-dotted and dotted vertical lines represent the experimental and LES predictions of the separation onset, respectively. In a similar way, the purple lines represent the reattachment locations. The LSB is predicted only on the suction side of the Rg-15 airfoil in both the experiments and numerical simulations.

3.2. Fidelity of RANS Models in Predicting Laminar Separation Bubbles (LSBs)

The efficiency of various turbulence models based on RANS equations in modeling the LSB is studied by conducting numerical simulations using nine different RANS-based turbulence models. The characteristics of the LSB over the RG-15 airfoil were studied by Oo et al. [6] in 2020, and large eddy simulations were performed to predict the location of LSBs and validated with experimental studies in a wind tunnel. The Cp distribution over the RG-15 airfoil at various angles of attack (0, 3, and 6 degrees) is reported in his work along with the separation onset and reattachment location of LSBs. The fidelity of each turbulence model in predicting the LSB is analyzed using three different parameters: the $C_p$ distribution, separation onset, and reattachment point from the $C_f$ distribution and aerodynamic coefficients. The separation onset and reattachment are calculated in the literature using $C_p$ distribution curves; therefore, it has limitations in predicting the exact locations as discussed in section A. The aerodynamic coefficients ($C_L$ and $C_D$) values are based on the shear stress and pressure distribution on the airfoil surface. $C_L$ values mainly depend on the difference in mean pressure on the pressure side and suction side of the airfoil and, therefore, it is relative. Therefore, using $C_L$ and $C_D$ values in assessing the prediction of a turbulence model is trivial. As a result, in this study, weightage is given to the $C_p$ distribution more than other factors as a mandate to assess the performance of the turbulence model.

The $C_p$ profile over Rg-15 for nine different RANS-based turbulence models at 0-degree angles of attack is shown in Figure 8. Only the transition turbulence models are capable of capturing the $C_p$ distribution accurately. Also, it can be observed that a negative value of $C_f$ (indicative of separation) is predicted only by the transition turbulence models. Therefore, the conventional RANS turbulence models failed to predict the LSBs, while transition-based models succeeded. Upon comparing the $C_p$ and $C_f$ distributions predicted by the numerical models with experiments, it can be observed that all the models predict the early onset of separation when compared to experiments. It should also be noted that even the LES simulation results (obtained from the literature) also predicted the early onset of the separation bubble. The $\gamma -R{e}_{\theta }$ transition model, the $\gamma$ transition model, and the $K-Kl-\omega$ model predict the $C_p$ distribution accurately on the pressure side of the airfoil. However, these models overpredicted the $C_p$ predictions on the suction side of the airfoil. Also, it should be noted that the $C_p$ plateau predicted by numerical simulation is very distinct, but the $C_p$ predictions in the separation region from experiments have a slight positive slope. This may be due to the fact that the separation bubble predicted by these transition models is stronger (large velocity gradients) than the actual LSBs observed in experiments. Introducing the low-Reynolds-number correction to the $\gamma$ transition model corrected this problem, but then the pressure side $C_p$ was overpredicted. In contrast with other transition models, the algebraic transition model predicted an LSB of low strength and, therefore, the $C_p$ distribution on the suction side has a better match with experiments. However, this model overpredicted the $C_p$ on the pressure side. Introducing the low-Reynolds-number correction to this model improved the overall $C_p$ predictions, but still the $C_p$ on the pressure side is overpredicted. This is because these models predicted a trailing edge separation bubble on the pressure side of the airfoil, which is not observed in the experimental and LES studies. Overall, the $C_p$ predictions of the $\gamma -R{e}_{\theta }$ transition model and the $\gamma$ transition model are closer to experimental prediction provided the $C_p$ on the suction side is overpredicted.

Table 1. Separation and reattachment locations of LSBs on the surface of Rg-15 airfoil.

AOA	0 Deg	% Diff from Exp	% Diff from LES	3 Deg		% Diff from Exp		% Diff from LES		6 Deg		% Diff from Exp		% Diff from LES
Turbulence Model	Seperation	Seperation	Seperation	Seperation	Reattachment	Seperation	Reattachment	Seperation	Reattachment	Seperation	Reattachment	Seperation	Reattachment	Seperation	Reattachment
Experiments	0.78	0.00	11.11	0.54	0.88	0.00	0.00	−1.10	0.00	0.24	0.49	0.00	0.00	22.45	11.36
LES	0.70	−10.00	0.00	0.55	0.88	1.11	0.00	0.00	0.00	0.20	0.44	−18.33	−10.20	0.00	0.00
SA	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
K-W	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
low re KW	-	-	-	-	-	-	-	-	-	0.20	0.37	−16.67	−23.67	2.04	−15.00
Trans KW	0.73	−6.41	3.99	0.53	0.87	−1.85	−0.68	−2.93	−0.68	0.20	0.46	−16.67	−6.73	2.04	3.86
gamma	0.69	−11.54	−1.71	0.51	1.00	−5.56	13.64	−6.59	13.64	0.19	0.47	−21.25	−4.69	−3.57	6.14
lowre gamma	0.65	−16.67	−7.41	0.51	1.00	−5.56	13.64	−6.59	13.64	0.18	0.47	−25.00	−4.49	−8.16	6.36
k-kl-w	0.68	−12.82	−3.13	0.47	1.00	−12.96	13.64	−13.92	13.64	0.16	0.50	−33.33	2.86	−18.37	14.55
algebraic	0.75	−3.85	6.84	0.55	0.78	1.85	−11.02	0.73	−11.02	0.20	0.40	−16.67	−19.18	2.04	−10.00
low re algebraic	0.77	−1.28	9.69	0.48	1.00	−11.11	13.64	−12.09	13.64	0.19	0.47	−20.83	−4.69	−3.06	6.14

shows the comparison of separation and reattachment positions of each transition model with experimental and LES observations. The separation onset is modeled more closely using the experimental predictions of the algebraic transition models with and without low $Re$ correction. However, as mentioned before, these models predicted non-physical separation bubbles on the pressure side. The $\gamma -R{e}_{\theta }$ transition model has the next closest prediction when compared with experiments with an error of 6.41%.

Figure 9 shows the $C_p$ and $C_f$ distribution over Rg-15 at a 3-degree angle of attack. In contrast with the observations made at a 0-degree angle of attack, the overall $C_p$ predictions of both the conventional and transitional turbulence models show a close match with the experimental predictions. All the models slightly overpredicted the $C_p$ on the pressure side and it shows deviations at the region of LSB. Similar to the 0-degree case, the $C_p$ predictions of conventional RANS models on the suction side are in line with the experimental measurements forward of the separation zone, but they deviate inside the recirculation zone due to its inability to predict LSBs. It can be observed from the $C_f$ distribution that all $\gamma -R{e}_{\theta }$ only transition models and algebraic transition models predict a separation bubble that reattaches to the surface of the airfoil (the same is also observed in experiments). Similar to the zero-degree case, the $\gamma -R{e}_{\theta }$ transition model predicted a strong separation bubble and, therefore, $C_p$ is slightly overpredicted locally in the LSB zone. The origin of the separation is accurately predicted by both models, but the reattachment point is predicted accurately by the $\gamma -R{e}_{\theta }$ transition model only.

When the angle of attack is increased to 6 degrees, the size of the LSB is reduced and its effect on overall $C_p$ distribution is more localized as shown in Figure 10. Except inside the LSB, the $C_p$ predictions of the numerical simulations and experiments are quite comparable. The conventional models failed to predict the LSB, and all of the transitional models predicted a stronger separation bubble than the one in the experiments. Therefore, $C_p$ is overpredicted by transition models in the LSB zone but is slightly underpredicted by conventional RANS models. From the $C_f$ distribution shown in Figure 10 and the information given in Table 1, it can be observed that all the transition models predict a closed LSB that reattaches to the airfoil surface. Most of the transition models predict the separation much earlier, when compared to experiments. The conventional $k\omega -\mathrm{S}\mathrm{S}\mathrm{T}$ turbulence model with low $Re$ correction and the algebraic transition model have better predictions of the separation onset location, but these models predict much earlier reattachment. Other transition models predict a reattachment location that is in between the predictions of experiments and LES, except the $k-kl-\omega$ model, which has a reattachment location aft of the experimental reattachment point. The legend for Figure 8, Figure 9, Figure 10 and Figure 11 is given in Figure 12.

Figure 8. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil at 0-degree angle of attack (See Figure 12 for legend).

Figure 8. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil at 0-degree angle of attack (See Figure 12 for legend).

Figure 9. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil at 3-degree angle of attack (See Figure 12 for legend).

Figure 9. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil at 3-degree angle of attack (See Figure 12 for legend).

Figure 10. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil at 6-degree angle of attack (See Figure 12 for legend).

Figure 10. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil at 6-degree angle of attack (See Figure 12 for legend).

Figure 11. $C_L$ (left) and $C_D$ (right) distribution over Rg-15 airfoil (See Figure 12 for legend).

Figure 11. $C_L$ (left) and $C_D$ (right) distribution over Rg-15 airfoil (See Figure 12 for legend).

Figure 12. Legend for Figure 8, Figure 9, Figure 10 and Figure 11.

Figure 12. Legend for Figure 8, Figure 9, Figure 10 and Figure 11.

Upon comparing the $C_p$ distribution at the three angles of attack considered in this study, it can be observed that the occurrence of laminar separation bubble has altered the $C_p$ distribution on the entire airfoil surface (both pressure and suction sides). For the 3-degree case, the influence of LSBs on the $C_p$ distribution outside the LSB is considerably reduced and is negligible for the 6-degree case. Thus, the open separation bubble formed in the 0-degree case has a global effect on the $C_p$ distribution and can be considered a “long LSB”. The bubbles formed at 3-degree and 6-degree angles of attack have only a local influence on the $C_p$ distribution (mainly in the separated flow zone) and can be considered “short LSB”. The thickness of the separation bubble at the onset of separation is very small and increases quite gradually, thus increasing the chances of experimental uncertainties in predicting the separation onset of an LSB. Previous researchers also reported that the LSB is quite unsteady and the reattachment point fluctuates, leading to a reattachment zone [27]. Therefore, considering the prediction of the reattachment location and $C_p$ distributions, the performance of the $\gamma -R{e}_{\theta }$ transition model classified it as a better model compared to the various transition-based models for the study of LSBs at lower angles of attack.

Figure 11 shows the $C_L$ and $C_D$ predictions of the nine different turbulence models considered in this study for all three angles of attack. The experimental $C_L$ and $C_D$ values are obtained from the studies of Selig [73] for $Re={10}^5$. It can be observed that LSB has a considerable influence in the prediction of aerodynamic coefficients. For zero-degree angle of attack, the algebraic transition models with and without low $Re$ correction have the closest match with experimental predictions followed by the $\gamma$ and $\gamma -R{e}_{\theta }$ transition models. For 3-degree angle of attack $\gamma -R{e}_{\theta }$ transition models have the best prediction and for 6- degrees the algebraic transition model, $\gamma$ and $\gamma -R{e}_{\theta }$ transition models predicts closely with experiments. The disparity in the predictions decreases with an increase in the angle of attack, and this may be due to the local and global effect of LSB on $C_p$ as discussed in previous paragraphs. Further, the difference in the $C_L$ and $C_D$ predictions among the transition models indicate the importance of predicting the LSB accurately. Considering the predictions at all three angles of attack, the predictions of the $\gamma -R{e}_{\theta }$ transition model seems more reliable.

3.3. Ice-Induced Separation Bubbles (ISBs)

The nature of ISBs with an increase in the angle of attack on the surface of an RG-15 airfoil is shown in Figure 13. The point of onset of separation seems unaffected by the angle of attack as observed by previous researchers [6,24]. When the angle of attack increased from 3 degrees to 6 degrees, it can be observed from Figure 13 that a small, separated flow region exists ahead of the main separation bubble. Such separation bubbles are not reported in the experimental and LES studies, where the separation point remains the same for all angles of attack [24]. However, careful observation of the $C_f$ distribution and velocity contour presented by Oo et al. [6] reveals the existence of small separated flow regions aft of the main separation bubble. The separation and reattachment points are obtained from the $C_p$ plots as shown in Figure 14. Like LSBs, $C_p$ flattening and negative $C_f$ values are observed inside the separation bubble. Interestingly, it should be noted that, unlike LSBs, ISBs are observed on both the suction side and pressure side of the iced airfoil. However, the experimental data are only available for the ISBs formed on the suction side. For ISBs, the flow separation originates from the tip of the leading-edge ice horn as shown in Figure 2, so the reattachment location is the major focus when studying ISBs, in addition to the $C_p$ distribution. Further, the aerodynamic lift coefficients of the iced Rg-15 airfoil predicted by the numerical solvers are compared with the LES simulation data provided by Oo et al. [6]. It is important to note that the conventional RANS turbulence models are also quite capable of predicting ISBs, unlike LSBs. Figure 15 shows the $C_p$ distribution over the ISB region predicted by the $K-\omega \mathrm{S}\mathrm{S}\mathrm{T}$ model (representative of conventional RANS) and the $\gamma -R{e}_{\theta }$ transition model (representative of transition models). It can be observed that, unlike the LSB case, the experiments predicted a very distinct plateau in the region of ISBs. The conventional models are not able to predict a distinct plateau like the transition models. Also, the length of the plateau region increases with the angle of attack.

3.4. Fidelity of RANS Models in Predicting Ice Induced Separation Bubbles (ISBs)

In this section, the capability of various RANS turbulence models in predicting an ISB is studied. Oo et al. [6] in 2020 performed experimental and LES studies on an iced Rg-15 airfoil to study the behavior of ISBs. The $C_p$ distribution over an RG-15 airfoil with mixed ice accretion and the separation onset and reattachment location of ISB is reported in his work. The $C_p$ and $C_f$ distributions over the Rg-15 airfoil with mixed ice accretion at a 0-degree angle of attack for nine different RANS-based turbulence models are shown in Figure 16. The leading edge of the clean Rg-15 airfoil starts at $\frac{X}{C}=0$ and the region forward of this point represents the geometry of the ice. The onset of separation starts on the surface of the ice at $\frac{X}{C}=-0.0183$ for all the turbulence models considered in this study, but the predictions of the reattachment location vary. Similar to the LSB case, the $C_p$ values are slightly overpredicted in the numerical studies, compared to experiments. Also, the values of the suction peak predicted by numerical simulations are quite high compared to experimental predictions. From Figure 15, it can be observed that the $C_p$ distribution in the ISB zone is well predicted by the transition models in comparison to the conventional models. However, from the $C_f$ distribution, it can be observed that the $\gamma$-based transition models and the $K-Kl-\omega$ model predict a trailing edge separation bubble like LSB, which is not observed in experiments and LES simulation cases. Therefore, the $C_p$ values outside the ISB zone show deviations from experiments. From

Table 2. Separation and reattachment locations of ISBs on the surface of Rg-15 airfoil.

AOA	0 Deg				3 Deg				6 Deg
Turbulence Model	Seperation	Reattachment	% Diff from Exp	% Diff from LES	Seperation	Reattachment	% Diff from Exp	% Diff from LES	Seperation	Reattachment	% Diff from Exp	% Diff from LES
Experiments	0.00	0.10	0.00	3.09	0.00	0.10	0.00	−15.97	0.00	0.14	0.00	6.92
LES	0.00	0.10	−3.00	0.00	0.00	0.12	19.00	0.00	0.00	0.13	−6.47	0.00
SA	−0.02	0.05	−46.00	−44.33	−0.02	0.06	−41.00	−50.42	−0.05	0.07	−51.44	−48.08
K-W	−0.02	0.06	−37.50	−35.57	−0.02	0.07	−30.50	−41.60	−0.05	0.07	−51.44	−48.08
low re KW	−0.02	0.08	−20.00	−17.53	−0.02	0.08	−25.00	−36.97	−0.05	0.09	−32.37	−27.69
Trans KW	−0.02	0.08	−22.50	−20.10	−0.02	0.08	−25.00	−36.97	−0.05	0.08	−45.32	−41.54
gamma	−0.02	0.07	−27.00	−24.74	−0.02	0.08	−18.00	−31.09	−0.05	0.09	−35.25	−30.77
lowre gamma	−0.02	0.07	−30.00	−27.84	−0.02	0.08	−23.00	−35.29	−0.05	0.09	−34.53	−30.00
k-kl-w	−0.02	0.08	−23.00	−20.62	−0.02	0.11	12.00	−5.88	−0.05	0.08	−43.17	−39.23
algebraic	−0.02	0.07	−35.00	−32.99	−0.02	0.07	−30.50	−41.60	−0.05	0.08	−43.88	−40.00
low re algebraic	−0.02	0.06	−44.50	−42.78	−0.02	0.06	−35.60	−45.88	−0.05	0.08	−45.68	−41.92

, it can be observed that the reattachment location predicted by the $K-\omega \mathrm{S}\mathrm{S}\mathrm{T}$ model with low $Re$ correction gives the closest match to experiments, followed by the $\gamma -R{e}_{\theta }$ transition model.

For the 3-degree angle-of-attack case, except the $K-Kl-\omega$ model, all other turbulence models predicted similar $C_p$ distributions outside the ISB zone as shown in Figure 17. This is because the $K-Kl-\omega$ turbulence model predicts a long trailing edge separation bubble as evident from the $C_f$ distribution shown in Figure 17. On the pressure side, the $C_p$ values are slightly overpredicted by numerical simulations. Inside the ISB, the $C_p$ distribution does not match the experimental predictions well, but downstream of the bubble, the predictions perform quite well. The reattachment location predicted by the $\gamma$ transition model is closer to the experimental prediction with a difference of 18% forward. The $C_p$ predictions are improved at a 6-degree angle of attack, except at regions inside the suction-side ISB as shown in Figure 18. An almost identical $C_p$ distribution is predicted by all the turbulence models considered in this study, except in regions inside the suction-side ISB. Likewise, the 0-degree case with the $K-\omega \mathrm{S}\mathrm{S}\mathrm{T}$ model with low $Re$ correction has the closest prediction of the reattachment point compared to experiments. However, it should be noted that the deviation from experiments is considerably higher by 32%, as given in Table 2.

Figure 19 shows the numerical $C_L$ and $C_D$ values predicted by different turbulence models for the angles of attack considered in this study. The $C_L$ values are compared with the LES data available in the literature by Oo et al. [6]. The $C_D$ data from LES simulation were not available in the literature and are therefore not compared with tests. It should be noted that the SA turbulence model predicts $C_L$ much better than any other model considered in this study for all angles of attack. Such a conclusion about the effectiveness of the SA model was also made by Costes et al. [42] during his studies on the NACA 23,012 airfoil with spanwise ridge ice at a 2-degree angle of attack. Similar observations were also made by Fajt [74]. The geometry of ISBs may affect the drag coefficient more evidently than lift due to the increase in shear stress, and the lack of experimental/LES drag data is a hindrance to assessing the capability of these models in detail. The legend for Figure 16, Figure 17 and Figure 18 is given in Figure 20.

Figure 16. $C_p$ (left) and $C_f$ (right) distributions over Rg-15 airfoil with mixed ice at 0-degree AOA (See Figure 20 for legend).

Figure 16. $C_p$ (left) and $C_f$ (right) distributions over Rg-15 airfoil with mixed ice at 0-degree AOA (See Figure 20 for legend).

Figure 17. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil with mixed ice at 3-degree AOA (See Figure 20 for legend).

Figure 17. $C_p$ (left) and $C_f$ (right) distribution over Rg-15 airfoil with mixed ice at 3-degree AOA (See Figure 20 for legend).

Figure 18. $C_p$ (left) and $C_f$ (right) distributions over Rg-15 airfoil with mixed ice at 6-degree AOA (See Figure 20 for legend).

Figure 18. $C_p$ (left) and $C_f$ (right) distributions over Rg-15 airfoil with mixed ice at 6-degree AOA (See Figure 20 for legend).

Figure 19. $C_L$ (left) and $C_D$ (right) distributions over Rg-15 airfoil with mixed ice accretion (See Figure 20 for legend).

Figure 19. $C_L$ (left) and $C_D$ (right) distributions over Rg-15 airfoil with mixed ice accretion (See Figure 20 for legend).

Figure 20. Legend for Figure 16, Figure 17, Figure 18 and Figure 19.

Figure 20. Legend for Figure 16, Figure 17, Figure 18 and Figure 19.

The RANS simulations underpredicted the reattachment location for all the angles of attack considered, and the growth of the separation bubble is limited to less than 15% as we increase the angle of attack from 0 to 6 degrees. However, the experimental studies predict 40% growth in the bubble length. Analogous observations were made by Bragg [24] during his experimental studies on the NACA 0012 airfoil as explained in Section 1.2. Thus, the RANS-based turbulence models predicted the reattachment much earlier than the actual scenario in the case of ISBs. This may be because the turbulence level predicted by the numerical solvers is higher than the experiments. These observations demand an improvement in the existing models for low-Reynolds-number ISB problems. Unlike LSBs, conventional turbulence is also capable of predicting ISBs, therefore the shear stress distribution does not vary considerably among the models except in cases where the use of transition models can lead to the prediction of nonphysical flow separation near the leading edge of the airfoil. Therefore, the low computational cost and simplicity associated with conventional RANS models make then a viable option to study ISBs with compromised accuracy.

4. Conclusions

The ability of nine distinct RANS-based turbulence models to predict the laminar separation bubble over the RG-15 airfoil at a Reynolds number of $1.07\times {10}^5$ was evaluated through comparison. Turbulence models with low Reynolds numbers that do not incorporate transition modeling are incapable of forecasting laminar separation bubbles; instead, they predict a fully attached flow at all angles of attack taken into account. Transition-based turbulence models are capable of predicting LSBs, but the predictions (location and size) differ among the models. Among the different transition models considered, the $\gamma -R{e}_{\theta }$ transition model has a better match with the experimental results for all of the angles of attack considered. The LSB migrates in the opposite direction of the flow as the angle of attack increases, while the separation bubble becomes smaller in both length and thickness. Conventional RANS-based models are also capable of predicting the ISB. Thus, unlike LSBs, the study of ISBs does not demand the use of any transition-based turbulence modes, while in some cases, the use of such models can lead to the prediction of nonphysical flow separation near the leading edge of the airfoil. The SA turbulence model can be considered the most favored turbulence model to study the aerodynamic coefficient of Iced airfoils but is not quite reliable enough to estimate the pressure distribution and bubble geometry. Also, the predictions of separation and reattachment locations by RANS models are subjected to uncertainty. Thus, the use of a single turbulence model that can be reliably used for the entire ice accretion simulation is difficult to select from the observations made in this study. Transition-based models can be used for the initial shots of the ice accretion simulation, but they need to be changed to conventional RANS models once an ISB forms. This helps save computational time and prevents the use of non-physical shear stress and heat flux values in the ice accretion calculations.