Prediction of Airfoil Stall Based on a Modified $\mathit{k}-\overline{{\mathit{v}}^2}-\mathit{\omega }$ Turbulence Model

Abstract

The accuracy of an airfoil stall prediction heavily depends on the computation of the separated shear layer. Capturing the strong non-equilibrium turbulence in the shear layer is crucial for the accuracy of a stall prediction. In this paper, different Reynolds-averaged Navier–Stokes turbulence models are adopted and compared for airfoil stall prediction. The results show that the separated shear layer fixed $k-\overline{v^2}-\omega$ (abbreviated as SPF $k-\overline{v^2}-\omega$) turbulence model captures the non-equilibrium turbulence in the separated shear layer well and gives satisfactory predictions of both thin-airfoil stall and trailing-edge stall. At small Reynolds numbers ($Re~{10}^5$), the relative error between the predicted $C_{L,max}$ of NACA64A010 by the SPF $k-\overline{v^2}-\omega$ model and the experimental data is less than 3.5%. At high Reynolds numbers ($Re~{10}^6$), the $C_{L,max}$ of NACA64A010 and NACA64A006 predicted by the SPF $k-\overline{v^2}-\omega$ model also has an error of less than 5.5% relative to the experimental data. The stall of the NACA0012 airfoil, which features trailing-edge stall, is also computed by the SPF $k-\overline{v^2}-\omega$ model. The SPF $k-\overline{v^2}-\omega$ model is also applied to a NACA0012 airfoil, which features trailing-edge stall and an error of $C_L$ relative to the experiment at $C_L>1.0$ is smaller than 3.5%. The SPF $k-\overline{v^2}-\omega$ model shows higher accuracy than other turbulence models.

1. Introduction

The stall phenomenon often appears in airplane wings and empennages, compressor blades and helicopter rotors [1]. The stall of the main wing can cause the aircraft to lose lift. Flow separation associated with stall can compromise control surface performance or even make an aircraft uncontrollable. The aerodynamic performance of a deeply stalled airfoil (i.e., airfoil operating at angles of attack larger than the stall angle of attack $AO{A}_{stall}$) is also important in some applications, such as rotorcraft [2,3] and wind turbines. Consequently, the characteristics of the airfoil stall must be accurately predicted in aerodynamic design.

The stalling characteristics of airfoils have been widely studied experimentally. McCullough and Gault [4] conducted detailed wind tunnel tests and classified airfoil stalls into three types: trailing-edge stall, leading-edge stall and thin-airfoil stall. The locations and evolution of separation bubbles are different for each stall type. Trailing-edge stall is preceded by the turbulent boundary layer separation point gradually moving forward as the angle of attack increases. This movement eventually results in a mild loss of the lift coefficient beyond $AO{A}_{stall}$. Leading-edge stall results from the burst of a laminar separation bubble that has formed near the leading edge of the airfoil. When $AOA$ increases, the flow separates from the upper surface of the airfoil abruptly and causes a sharp loss in the lift. In a thin-airfoil stall, a laminar separation bubble forms at the leading edge, and the airfoil stall is preceded by the gradual growth of the separation bubble as the angle of attack increases. Bragg [5] experimentally studied the stalling characteristics of NACA64A010 and found that the airfoil exhibited thin-airfoil stall features. At $AOA=4°$, the slope of the $C_L-AOA$ curve decreased significantly. Bragg attributed the loss in slope to the appearance of a leading-edge separation bubble at $AOA=4°$. Tani [6] suggested that this is because the short separation bubble becomes a long bubble. At $AOA=10°$, the lift coefficient attained its maximum and started to decrease. At $AOA=15°$, the lift started to recover, and bluff body vortex shedding was observed by Bragg [5].

Accurate numerical simulation of airfoil stall is challenging. Complex phenomena, including an adverse pressure gradient, boundary layer separation and transition, are involved in the flow past a stalled airfoil [7]. Many unsteady simulations have numerically predicted airfoil stall performance. Rodriguez et al. [8] conducted a direct numerical simulation (DNS) study to learn the flow details around the NACA0012 airfoil at a large angle of attack. The results showed that the flow had three main characteristics: boundary layer separation at the leading edge, Kelvin–Helmholtz (K-H) instability in the separated shear layer, which eventually leads to transition, and large eddies, which induce strong momentum exchange after the transition point. Mary and Sagaut [9] conducted a large eddy simulation (LES) of the turbulent flow past an A-Airfoil. This study showed that LES could give relatively accurate results for the velocity profiles and the pressure distribution at large $AOAs$. At $AOA={13}^{\circ }$, the velocity profiles given by LES at $x/c\in \left[0.85,1.00\right],y/c<0.03$ reflected the flow separation successfully and differed from the experimental data by only 5%. Additionally, the relative error of the computed pressure coefficient and experimental data under the separation bubble ($x/c\in \left[0.85,1.00\right]$) was less than 5%. Im and Zha [10] used a delayed, detached eddy simulation (DDES) with the Spalart–Allmaras turbulence model to predict the force coefficient of the NACA0012 airfoil. The lift coefficient and the drag coefficient were accurately predicted at angles of attack up to $60°$. Kawai and Fujii [11] also used the Reynolds-averaged Navier–Stokes/large eddy simulation (RANS/LES) hybrid method to simulate the flow past the NACA64A006 airfoil at $AOAs$ up to $11°$. These unsteady simulations successfully predicted the growth of the laminar separation bubble. The results of the high-fidelity unsteady simulation methods were reported to be satisfactory. However, the large computational cost makes these methods unsuitable for the daily design process.

Currently, the steady Reynolds-averaged Navier–Stokes (RANS) method is widely used in practical engineering. The computational cost of the RANS method is relatively low, and it can provide an acceptable force coefficient at a small angle of attack. Although the flow past a stalled airfoil is essentially unsteady, using the steady RANS method to obtain a relatively accurate lift coefficient or other averaged features of the flow field is still quite desirable for daily aerodynamic design. The turbulence model should be used in the RANS method to approximate the unclosed Reynolds stress in the averaged Navier–Stokes equations. Two of the most commonly used turbulence models for aeronautics applications are the SST turbulence model and the SA turbulence model. The SST turbulence model gives accurate results in the near-wall layers. The zonal formulation based on blending functions makes the model insensitive to the free stream $\omega$ value. The SST turbulence model gives a relatively accurate separation point, but it underpredicts the level of turbulence stresses in the separated shear layer. The position of reattachment also disagrees with the experimental data [12]. The SA turbulence model contains a single PDE. Consequently, it is easy to implement the SA turbulence model in the existing code. The SA turbulence model is also proven to be numerically well behaved. The single PDE of the SA turbulence model can be applied to both incompressible and compressible flow, which makes it convenient for aeronautic applications of a wide speed range. However, the model uses a prescribed turbulence scale. This simplification compromises the model’s ability to predict the separated flow [13]. In Section 3 of this paper, the SST and SA model are used to predict the airfoil stall. The results show that their ability in stall prediction is limited.

Indeed, most RANS turbulence models have poor performance in predicting airfoil stalling characteristics. Rizzetta and Visbal [14] applied the Baldwin–Lomax algebraic model and two-equation $k-ϵ$ model to compute the stall performance of the Sikorsky SSC-A09 airfoil and found that both models gave poor results. Genç [15] used the $k-\omega$ shear–stress transport (SST) model to simulate the stalling characteristics of the NACA64A006 airfoil and found that the model seriously underpredicted the lift coefficient when $AOA$ was large.

RANS methods with transition prediction may more accurately evaluate the airfoil’s stalling characteristics. Catalano and Tognaccini [16] developed an empirical method based on the analysis of the laminar separation bubble to approximate the transition point on an airfoil. The RANS simulation using a fixed transition point obtained by the empirical method gave satisfactory stall prediction of the SD7003 airfoil. Wokoeck et al. [17] introduced a transition prediction based on linear stability theory to RANS simulation and successfully predicted the laminar separation bubble on an airfoil surface. However, the maximum lift and the stalling angle of attack were not accurately computed. Genç [15] applied the $k_T-k_L-\omega$ transitional model [18,19,20] to predict the stalling characteristics of the NACA64A006 airfoil. The results showed that the size of the laminar separation bubble at small angles of attack agreed well with the experimental results. Morgado et al. [21] used the same $k_T-k_L-\omega$ transitional model to simulate the stall of E387 and S1223 airfoils. The results showed that the model could give a promising prediction of $C_L$ at large $AOAs$, provided that the $AOAs$ were less than $AO{A}_{stall}$. In Section 3, an alternative version of the $k_T-k_L-\omega$ transitional model, called the $k-\overline{v^2}-\omega$ model [20], is used to simulate the airfoil stall. The $k-\overline{v^2}-\omega$ model and the $k_T-k_L-\omega$ model have similar characteristics.

Nearly all turbulence models mentioned above are calibrated by wall-bounded flow without massive flow separation. The lack of calibration by the separated flow leads to poor capability in predicting massive separation in the flow past an airfoil operating at $AOA\ge AO{A}_{stall}$. Rumsey [22] proposed that most RANS turbulence models underpredict the eddy viscosity in the shear layer region of a separation bubble and thus lead to insufficient mixing and delayed reattachment after flow separation. On the other hand, the non-equilibrium turbulence (which is characterized by a large ratio between the rates of turbulence kinetic energy production and the rates of dissipation) is prominent in the region of flow separation, where the adverse pressure gradient is strong [23]. In the research conducted by Fang [24,25], it was also found that the non-equilibrium turbulence appears in the separated shear layer. Based on these observations, Rumsey [26] developed a modified $k-\omega$ model, which recognizes the shear layer region by its high $P_k/\epsilon$, enhances the destruction of $\omega$ in this region and produces a large eddy viscosity in the separated shear layer. This modification substantially improves the turbulence model’s ability to predict massive separation in some cases, such as the flow past a bump [26]. Li et al. [27,28,29] used a similar non-equilibrium modification function to improve the SST, $k-\overline{v^2}-\omega$ and Wilcox 2006 $k-\omega$ models. The fixed term in these modified models, based on the non-equilibrium features of turbulence, can augment the eddy viscosity in the region where the non-equilibrium turbulence and the shear are strong. The improved $k-\overline{v^2}-\omega$ model is called the separated shear layer fixed $k-\overline{v^2}-\omega$ model. The separated shear layer fixed $k-\overline{v^2}-\omega$ model successfully predicted the non-equilibrium turbulence in the separated shear layer produced by the ice accretion at the leading edge of the airfoil. The size of the separation bubble behind the ice corn also agreed well with the experiment. Zhang et al. [30] applied the separated shear layer fixed $k-\overline{v^2}-\omega$ model to study the multi-element airfoils and a complex three-dimensional JAXA standard model with high-lift devices. The results showed that the separated shear layer fixed $k-\overline{v^2}-\omega$ model performed well in predicting the separation bubble in the flow field. The lift coefficients of the multi-element airfoils and the three-dimensional JAXA standard model also agreed with the experiment results.

To summarize, the airfoil stall is an important phenomenon in wing design and turbine blade design. The stall performance is critical to the overall performance of an aircraft. An efficient method to predict airfoil stall is highly desirable. However, current high-fidelity methods (DNS and LES) are inefficient. RANS is significantly cheaper, but the mediocre performance of traditional turbulence models in separation flow prediction limits RANS’s application in stall prediction. The facts listed above motivated the authors to develop a turbulence model that is capable of stall prediction to serve for daily engineering design.

The authors believe that the SPF $k-\overline{v^2}-\omega$ model has great potential in stall prediction, in which the separated shear layer plays an important role. In this paper, the separated shear layer fixed $k-\overline{v^2}-\omega$ model’s (SPF $k-\overline{v^2}-\omega$ model) ability to simulate airfoil stall is verified. The model is applied to several airfoil stall cases, including the thin-airfoil stall of the NACA64A010 airfoil at both low and high Reynolds numbers, thin-airfoil stall of the NACA64A006 airfoil and trailing-edge stall of the NACA0012 airfoil. The results of the test cases show that the SPF $k-\overline{v^2}-\omega$ model performs better than traditional turbulence models in predicting the stall characteristics of different airfoil stall types. In thin-airfoil stall and trailing-edge stall, the separated shear layer where the non-equilibrium turbulence is strong develops and plays an important role in determining the lift coefficient of the airfoil. Since the modification term in the SPF $k-\overline{v^2}-\omega$ model is aimed at the separated shear layer, the test cases contain both airfoils exhibiting thin-airfoil stall (NACA64A006 and NACA64A010) and airfoils exhibiting trailing-edge stall (NACA0012). These specific airfoils (NACA64A006, NACA64A010 and NACA0012) are chosen because their experimental data are relatively abundant in the current literature.

The authors believe that the SPF $k-\overline{v^2}-\omega$ model can be used to evaluate the stall performance of a newly designed airfoil and can also be used to conduct numerical optimization of airfoils, for which the stall performance is important.

2. Governing Equations and Numerical Method

2.1. The Reynolds-Averaged Navier–Stokes Equations

Direct numerical simulation of fluid flow using Navier–Stokes equations is numerically expensive or even prohibitive due to the strong fluctuation and the multiple scales of turbulence. Averaging can be applied to Navier–Stokes equations to reduce the cost of simulation. The Reynolds-averaged equations can be written as follows [31]. Note that in all equations listed in this paper, the summation convention is used:

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\rho }\overline{u_j}\right)=0$$

(1)

$$\frac{\partial \left(\overline{\rho }\overline{u_i}\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\rho }\overline{u_i}\overline{u_j}\right)=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(2)

$$\frac{\partial \left(\overline{\rho }\overline{E}\right)}{\partial t}+\frac{\partial [\overline{\rho }\overline{u_j}(\overline{E}+\frac{\overline{p}}{\overline{\rho }})]}{\partial x_j}=\frac{\partial }{\partial x_j}\left({\overline{u}}_i{\tau }_{ij}+\frac{\mu }{\left(\gamma -1\right)Pr}\frac{\partial }{\partial x_j}(\frac{\gamma \overline{p}}{\overline{\rho }})+\overline{{\dot{q}}_j}\right)$$

(3)

where $Pr,\gamma$ are constants and

$$\overline{E}=\frac{\overline{p}}{\left(\gamma -1\right)\overline{\rho }}+\frac{1}{2}\left(\overline{u_i}\overline{u_i}\right)$$

(4)

$$\overline{{\dot{q}}_j}=\frac{\overline{\rho }{\nu }_T}{\left(\gamma -1\right)P{r}_T}\frac{\partial }{\partial x_j}\left(\frac{\gamma \overline{p}}{\overline{\rho }}\right),P{r}_T=const.$$

(5)

All variables with a bar denote averaged variables. The vector $\overline{{\dot{q}}_j}$ describes the heat diffusion caused by turbulent motion. The stress tensor ${\tau }_{ij}$ in Equations (2) and (3) can be written as

$${\tau }_{ij}=\mu \left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)-\frac{2}{3}\mu \frac{\partial \overline{u_k}}{\partial x_k}{\delta }_{ij}-\overline{\rho }\overline{u_i^{\prime }{u_j}^{\prime }},\mu =const.$$

(6)

where $u_i^{\prime }$ denotes the fluctuation velocity, which is the difference between the instantaneous velocity and the averaged velocity:

$$u_i^{\prime }=u_i-\overline{u_i}$$

(7)

After substituting Equations (4) and (5) into Equations (2) and (3), one can find that there are 11 dependent variables. (The vector $\overline{u_j}$ has 3 components, and the symmetric tensor $\overline{u_i^{\prime }{u_j}^{\prime }}$ called Reynolds stress tensor has 6 components.):

$$\overline{\rho },\overline{p},\overline{u_i},\overline{u_i^{\prime }{u_j}^{\prime }}$$

whereas there are only five partial differential equations in Equations (1)–(3). Consequently, the averaged Navier–Stokes equations are unclosed. To reduce the number of dependent variables, people use the Boussinesq approximation to express $\overline{u_i^{\prime }{u_j}^{\prime }}$:

$$-\overline{u_i^{\prime }{u_j}^{\prime }}={\nu }_T\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}-\frac{2}{3}\frac{\partial \overline{u_k}}{\partial x_k}{\delta }_{ij}\right)-\frac{1}{3}\overline{u_k^{\prime }{u_k}^{\prime }}{\delta }_{ij}$$

(8)

Equation (8) is an analog of the first two terms on the right-hand side of Equation (6). Physically speaking, the approximation in Equation (8) implies that the mechanism of the momentum transfer caused by the velocity fluctuation is the same as the momentum transfer caused by the molecular transport (which is depicted by the first two terms on the right-hand-side of Equation (6)).

The number of dependent variables is reduced to seven after substituting Equations (4), (5) and (8) into Equations (1)–(3):

$$\overline{\rho },\overline{p},\overline{u_i},{\nu }_T,\overline{u_k^{\prime }{u_k}^{\prime }}$$

where ${\nu }_T$ is called the eddy viscosity and $\frac{1}{2}\overline{u_k^{\prime }{u_k}^{\prime }}$ is called the turbulent kinetic energy. Then, the turbulence model is introduced to compute the turbulent kinetic energy and eddy viscosity. The turbulence model greatly affects the accuracy of the flow simulation and should be chosen or constructed based on the physical properties of the flow.

2.2. SPF $k-\overline{v^2}-\omega$ Turbulence Model

The SPF $k-\overline{v^2}-\omega$ turbulence model proposed by Li et al. [29] is used to predict the stall performance of airfoils in this paper. This is because the separated shear layer appears and develops in the thin-airfoil stall and trailing-edge stall, and the SPF $k-\overline{v^2}-\omega$ model is designed to accurately compute the separated shear layer. Furthermore, the successful application of the SPF $k-\overline{v^2}-\omega$ model to the complex separated flow around an iced airfoil [29] and multi-element airfoil [30] makes the choice of the SPF $k-\overline{v^2}-\omega$ model in this work more reasonable. The transport equations of the SPF $k-\overline{v^2}-\omega$ model include the equation of the fluctuation kinetic energy $k$ (including both turbulent kinetic energy and pre-transitional velocity fluctuation), the equation of the kinetic energy of full turbulent velocity fluctuation $\overline{v^2}$ and the equation of the specific dissipation rate $\omega$. These three equations correspond to Equations (9), (10) and (11), respectively. The eddy viscosity is determined by $k, \overline{v^2},\omega ,\overline{u_j}$ and $\overline{\rho }$ through the algebraic relation in Equation (12), and the turbulence kinetic energy $\frac{1}{2}\overline{u_l{u}_l}$ equals $\overline{v^2}$:

$$\frac{\partial k}{\partial t}+\overline{u_j}\frac{\partial k}{\partial x_j}=\frac{1}{\overline{\rho }}{P}_k-\min \left(\omega k,\omega \overline{v^2}\right)-\frac{1}{\rho }{D}_k+\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{\rho {\alpha }_T}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]$$

(9)

$$\frac{\partial \overline{v^2}}{\partial t}+\overline{u_j}\frac{\partial \overline{v^2}}{\partial x_j}=\frac{1}{\overline{\rho }}{P}_{\overline{v^2}}+R_{BP}+R_{NAT}-\omega \overline{v^2}-\frac{1}{\rho }{D}_{\overline{v^2}}+\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{\rho {\alpha }_T}{{\sigma }_k}\right)\frac{\partial \overline{v^2}}{\partial x_j}\right]$$

(10)

$$\frac{\partial \omega }{\partial t}+\overline{u_j}\frac{\partial \omega }{\partial x_j}=\frac{1}{\overline{\rho }}{P}_{\omega }+\left(\frac{C_{\omega R}}{f_W}-1\right)\frac{\omega }{\overline{v^2}}\left(R_{BP}+R_{NAT}\right)-f_{NE}{C}_{\omega 2}{\omega }^2{f}_W^2+2{\beta }^{*}\left(1-F_1^{*}\right){\sigma }_{\omega 2}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}+\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{\rho {\alpha }_T}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right]$$

(11)

$${\nu }_T=g\left(k, \overline{v^2},\omega ,\overline{u_j},\overline{\rho }\right)$$

(12)

$$\frac{1}{2}\overline{u_l{u}_l}=\overline{v^2}$$

(13)

The dependent variables of Equations (9)–(13) are

$$\overline{\rho }, \overline{u_j},k, \overline{v^2},\omega ,{\nu }_T,\overline{u_l{u}_l}$$

The definitions of other variables can be found in the Appendix A of this paper. Substitute Equations (4) (5) and (8) into Equations (1)–(3) and combine the resulting equation with Equations (9)–(13), and we obtain a system of equations that has eight dependent variables:

$$\overline{\rho },\overline{p},\overline{u_j},{\nu }_T, \overline{u_l{u}_l}, k, \overline{v^2},\omega$$

There are also eight equations in the system. Consequently, the system is closed and can be solved numerically.

The major difference between the SPF $k-\overline{v^2}-\omega$ turbulence model and the original $k-\overline{v^2}-\omega$ turbulence model is the modification factor $f_{NE}$ added in the SPF $k-\overline{v^2}-\omega$ model. The modification factor $f_{NE}$ is in the destruction term of $\omega$ and has the following expression:

$$f_{NE}=\{\begin{array}{c}1.0, R{e}_{\mathsf{\Omega }}<\frac{1.0}{Th{r}_{SSL}{\mathsf{\Gamma }}_{SSL}}\\ 300R{e}_{\mathsf{\Omega }}{\mathsf{\Gamma }}_{SSL},\frac{1.0}{Th{r}_{SSL}{\mathsf{\Gamma }}_{SSL}}\le R{e}_{\mathsf{\Omega }}<\frac{Mag}{Th{r}_{SSL}{\mathsf{\Gamma }}_{SSL}}\\ Mag, R{e}_{\mathsf{\Omega }}\ge \frac{Mag}{Th{r}_{SSL}{\mathsf{\Gamma }}_{SSL}}\end{array}$$

(14)

$${\mathsf{\Gamma }}_{SSL}=\frac{1}{1+e^{-10{\left(\frac{P_{\overline{v^2}}}{\epsilon }-C_{SSL}\right)}^2}}$$

(15)

$$R{e}_{\mathsf{\Omega }}=\frac{d^2\mathsf{\Omega }}{\nu }$$

(16)

$$\frac{P_{\overline{v^2}}}{\epsilon }=\frac{{\mu }_{T,s}{S}^2}{\rho \overline{v^2}\omega }$$

(17)

where $d$ is the distance between the point considered and the wall. The $f_{NE}$ is larger than one in regions where $\frac{P_{\overline{v^2}}}{\epsilon }-C_{SSL}$ and $R{e}_{\mathsf{\Omega }}>\frac{1}{Th{r}_{SSL}}$, while $f_{NE}$ remains one in other regions. The local value of $R{e}_{\mathsf{\Omega }}$ and $P_{\overline{v^2}}/\epsilon$ decides the magnitude of the modification factor $f_{NE}$. The diagram of the switch function ${\mathsf{\Gamma }}_{SSL}$ in Figure 1 illustrates that ${\mathsf{\Gamma }}_{SSL}$ approaches one when $\frac{P_{\overline{v^2}}}{\epsilon }>C_{SSL}$ and approaches zero when $\frac{P_{\overline{v^2}}}{\epsilon }<C_{SSL}$.

In the fluid region where $\frac{P_{\overline{v^2}}}{\epsilon }<C_{SSL}$, ${\mathsf{\Gamma }}_{SSL}$ approaches zero, and $\frac{1}{Th{r}_{SSL}{\mathsf{\Gamma }}_{SSL}}$ approaches $\infty$. According to Equation (14), for all $R{e}_{\mathsf{\Omega }}$ values that are physically meaningful, $f_{NE}=1.0$. In the region where $\frac{P_{\overline{v^2}}}{\epsilon }>C_{SSL}$, ${\mathsf{\Gamma }}_{SSL}$ approaches one, and $\frac{1}{Th{r}_{SSL}{\mathsf{\Gamma }}_{SSL}}$ approaches $\frac{1}{Th{r}_{SSL}}$. Equation (14) then indicates that $f_{NE}>1.0$ if $R{e}_{\mathsf{\Omega }}>\frac{1}{Th{r}_{SSL}}$.

A more physical description is as follows. The modification term $f_{NE}$ is increased in the region away from the wall where the shear presents ($R{e}_{\mathsf{\Omega }}>1/Th{r}_{SSL}$) and the non-equilibrium turbulence is strong ($\frac{P_{\overline{v^2}}}{\epsilon }>C_{SSL}$). The separated shear layer extending from the leading edge of a stalled airfoil has the characteristic mentioned above and $f_{NE}$ is increased there. If $f_{NE}$ is increased, the destruction term $f_{NE}{C}_{\omega 2}{\omega }^2{f}_W^2$ in Equation (16) will increase. Therefore, the local value of $\omega$ tends to decrease. The eddy viscosity ${\nu }_T$ is approximately proportional to $\frac{k}{\omega }$. Therefore, due to the decrease in $\omega$, the local value of ${\nu }_T$ tends to increase. Consequently, the modification term $f_{NE}$ can increase ${\nu }_T$ in the separated shear layer. This characteristic makes the SPF $k-\overline{v^2}-\omega$ model able to avoid the underpredicted ${\nu }_T$ in the separated shear layer suffered by many traditional turbulence models [22].

According to the calibration of Li et al. [27] and the physical properties of the separated shear layer [22], the parameters in the modification term are reported in

Table 1. The value of the parameters of the modification term $f_{NE}$.

Parameter	Value
$Th{r}_{SSL}$	300
$Mag$	1.8
$C_{SSL}$	2.1

.

2.3. Numerical Method

All the computations in the following sections are carried out using a finite volume solver (CFL3D version 6.7) [31]. Roe’s flux difference splitting method is used to calculate the inviscid flux at the cell interface:

$${\widehat{F}}_{i+\frac{1}{2}}=\frac{1}{2}\left[\widehat{F}\left(q_L\right)+\widehat{F}\left(q_R\right)-\left|\tilde{A}\right|\left(q_R-q_L\right)\right]$$

(18)

Here, $q_L$ and $q_R$ are reconstructed by a third-order upwind-biased scheme:

$$q_L=q_i+\frac{1}{6}\left(q_i-q_{i-1}\right)+\frac{1}{3}\left(q_{i+1}-q_i\right)$$

(19)

$$q_R=q_{i+1}-\frac{1}{6}\left(q_{i+2}-q_{i+1}\right)-\frac{1}{3}\left(q_{i+1}-q_i\right)$$

(20)

The scalar tridiagonal inversion method is used for matrix inversion [31]. Local time stepping and multigrid are used to accelerate the convergence. The implicit approximate factorization method is used to solve the three transport equations of the SPF $k-\overline{v^2}-\omega$ model.

The boundary conditions are as follows: the no-slip condition is imposed on the surface of the airfoil, and the free-stream condition is imposed on the outer boundary of the physical domain. For the three transport equations of the SPF $k-\overline{v^2}-\omega$ model, $k$, $\overline{v^2}$ and the normal gradient of $\omega$ are set to zero on the no-slip wall boundaries. At the outer boundary of the physical domain, $k$ is given, and $\omega$ is computed to match the specified freestream eddy viscosity.

2.4. Grid Convergence Study

Three C-type grids of the NACA64A010 airfoil with different cell numbers are used for the grid convergence study of the SPF $k-\overline{v^2}-\omega$ model. The schematic diagram of the computational domain and the boundary conditions are illustrated in Figure 2a. The parameters $r, d_1, d_2$ describe the scales of the computational domain and are equal to 40 times the chord length of the airfoil. The fine grid has 681 and 189 grid points in the circumferential direction and the wall-normal direction, respectively. The medium grid and coarse grid have $481\times 133$ and $341\times 97$ grid points, respectively. Figure 2b shows the fine grid used in this paper. The lift coefficients of the NACA64A010 airfoil at $Ma=0.12$, $Re=3\times {10}^5$ and $AOA=11°$ obtained on these grids are shown in Figure 2c. The lift coefficients are plotted as a function of $\sqrt{1/N}\times {10}^3$, where $N$ is the number of grid cells. In 2D cases, $\sqrt{1/N}$ is proportional to the averaged grid spacing. As shown in Figure 2c, the lift coefficient approached the experimental results [5] as the grid became finer. Consequently, the model exhibited good grid convergent behavior.

3. Numerical Results from Stalled Airfoils

Three test cases are presented in this section. In the first two cases, the SPF $k-\overline{v^2}-\omega$ model is used to predict the stall characteristics of the NACA64A010 and NACA64A006 airfoils. Both airfoils display thin-airfoil stall. In the third test case, the turbulence model is applied to calculate the aerodynamic performance of the NACA0012 airfoil at large $AOAs$, which undergoes trailing-edge stall at high Reynolds numbers.

3.1. Stall Prediction of the NACA64A010 Airfoil

3.1.1. Low Reynolds Number Case

Different turbulence models, including SA [32], SST [33], $k-\overline{v^2}-\omega$ [20] and SPF $k-\overline{v^2}-\omega$ (termed the SPF model later in this article), were used to predict the stall behavior of the NACA64A010 airfoil, which underwent thin-airfoil stall. The computation was carried out on the fine grid shown in Figure 2a. The freestream Mach number was 0.12, and the Reynolds number (based on the chord length) was $3\times {10}^5$. Figure 3 shows the $C_L-AOA$ curves obtained by different turbulence models. The SA and SST models underpredicted $C_L$ at relatively large $AOA$s. $C_{L,max}$ and $AO{A}_{stall}$ were also greatly underpredicted by the SST model. The $k-\overline{v^2}-\omega$ model performed better than the SA and SST models, but the predicted $C_L$ value was still much lower than the experimental value. Unlike the experimental data, the lift curve given by $k-\overline{v^2}-\omega$ did not rise at $AOA>15°$. The SPF model predicted $C_{L,max}$ and $AO{A}_{stall}$ more accurately than the other models, and the $C_L-AOA$ curve rose as $AOA$ increased at $AOA>15°$. Moreover, the relative error between the predicted $C_{L,max}$ and experiment was smaller than 3.5%.

Figure 4 shows the separation bubbles obtained by the SST, $k-\overline{v^2}-\omega$ and SPF at $AOA=17°$. The heights of the separation bubbles ($h/c$) given by the SST, $k-\overline{v^2}-\omega$ and SPF models were 0.39, 0.33 and 0.29, respectively. The SPF model gave the smallest separation bubble. Since the flow with higher momentum encouraged flow reattachment and made the separation bubble shrink, the smaller separation bubble obtained by SPF indicated a flow with greater momentum near the upper surface of the airfoil. The velocity profiles obtained at different stations on the airfoil in Figure 5 confirm that the flow between the airfoil surface and the separated shear layer obtained by the SPF model has a greater velocity compared to the other models. The greater momentum obtained by the SPF model induced a lower pressure coefficient $C_p$ on the upper surface of the airfoil, as shown in Figure 6. The longer interval with lower $C_p$ on the airfoil surface directly led to a larger $C_L$ obtained by the SPF model at $AOA=17°$, as shown in Figure 3.

The reason why the SPF model gave a flow with higher momentum between the airfoil surface and the separated shear layer is explained here. In the context of turbulence models based on the Boussinesq hypothesis, a larger ${\nu }_T$ can intensify the momentum exchange. Figure 7 shows the contour of ${\nu }_T$ given by the SST, $k-\overline{v^2}-\omega$ and SPF models. The SPF model obtained a larger region where ${\nu }_T$ was greater than 2000 compared with the SST and $k-\overline{v^2}-\omega$ models. Consequently, the momentum exchange was stronger near the separated shear layer in the flow field predicted by the SPF model. Eventually, the stronger momentum exchange will increase the momentum of the flow near the wall. Accordingly, the SPF model gave a shear layer separated from the leading edge, with greater momentum between the separated shear layer and the airfoil surface compared with the SST and $k-\overline{v^2}-\omega$ models.

The larger region with a high ${\nu }_T$ value obtained by the SPF model was due to the modification introduced to it. Figure 8 shows the contour of $P_{\overline{v^2}}/\epsilon$ at $AOA=17°$ obtained by the SST, $k-\overline{v^2}-\omega$ and SPF models. Typically, $P_{\overline{v^2}}/\epsilon$ is as high as 3–4 in some part of the separated shear layer [26]. Figure 8a shows that $P_{\overline{v^2}}/\epsilon$ was smaller than two along the separated shear layer, indicating a failure in the prediction of non-equilibrium turbulence. In Figure 8b, $k-\overline{v^2}-\omega$ predicted a small region in which $P_{\overline{v^2}}/\epsilon$ was greater than two in the separated shear layer. However, this region only extended to $x/c=0.05$. On the other hand, the SPF model predicted that $P_{\overline{v^2}}/\epsilon$ was greater than two along the separated shear layer, as shown in Figure 8c. There was also a region in the separated shear layer that extended to $x/c=0.25$, where $P_{\overline{v^2}}/\epsilon$ was as high as 3–5, which is in agreement with the observation of Rumsey [23]. Consequently, among the three turbulence models considered here, only the SPF model recognized the strong non-equilibrium turbulence in the separated shear layer. The recognition of non-equilibrium turbulence activated the modification in the SPF model through a rapid increase in ${\mathsf{\Gamma }}_{SSL}$ and then $f_{NE}$. As shown in Equation (11), the magnified $f_{NE}$ augmented the destruction term of $\omega$. Since ${\nu }_T\propto \frac{1}{\omega }$, the decrease in $\omega$ caused by $f_{NE}$ led to an increase in ${\nu }_T$. Finally, a higher ${\nu }_T$ induced a lower $C_p$ on the airfoil surface and a larger $C_L$, as mentioned earlier. Relying on the mechanism discussed above, the SPF model gave a more accurate prediction of $C_L$ at large $AOA$s.

3.1.2. High Reynolds Number Case

In this section, the NACA64A010 airfoil at the freestream Mach number $0.167$ and the Reynolds number $4.1\times {10}^6$ is studied. The computation was carried out on the fine grid shown in Figure 2a. This case was used to test the SPF model’s applicability to thin-airfoil stall problems at high Reynolds numbers. All numerical results were compared with the experiment conducted by Peterson [34]. Figure 9 shows the $C_L-AOA$ curves obtained by different turbulence models, including the SA, SST, $k-\overline{v^2}-\omega$ and SPF models. SA, SST and $k-\overline{v^2}-\omega$ substantially overpredicted $AO{A}_{stall}$ and $C_{L,max}$. The $AO{A}_{stall}$ and $C_{L,max}$ given by SPF were $9°$ and 0.98, respectively, and were close to $AO{A}_{stall}=9.5°$ and $C_{L,max}=1.02$, which were obtained by the experiment. The relative error between the computed $C_{L,max}$ using the SPF model and the experimental result was smaller than 4.0%. On the other hand, the decrease in the lift coefficient after the stall was slightly overpredicted by the SPF model, which led to a more abrupt stall than was observed in the experiment. However, the result given by the SPF model obviously showed the best agreement with the experimental data among the four turbulence models.

Figure 10 shows the pressure distributions on the airfoil at $AOA=4°$ and $AOA=6°$. The results obtained by SPF, SST and SA agreed well with the experimental results. The $k-\overline{v^2}-\omega$ model gave the pressure distribution with a zigzag pattern, which showed poor convergence in the steady RANS computation. Figure 11 shows the velocity profiles at $x/c=80\%$ obtained by the SPF model. The scatter diagram represents the experimental data, and the curves represent the results obtained by the SPF model. The computational results agreed well with the experiment at $AOA=4°$, $AOA=6°$ and $AOA=9°$. However, the velocity profile obtained by the SPF model was slightly different from the experiment at $AOA=9.5°$. This difference corresponded to the sharp stall in Figure 9 calculated by the SPF model. In conclusion, the SPF model performed better than the SST, SA and $k-\overline{v^2}-\omega$ models at predicting the stall characteristics of the NACA64A010 airfoil at high Reynolds numbers.

3.2. Stall Prediction of the NACA64A006 Airfoil

The SPF model was used to predict the stalling behavior of the NACA64A006 airfoil, which exhibited thin-airfoil stalling characteristics. The experimental results obtained by McCullough et al. [35] were used to test the accuracy of the numerical results. The freestream Mach number was $0.167$, and the Reynolds number (based on the chord length) was $5.8\times {10}^6$. A C-type grid like Figure 2a was adopted. This case was chosen to test whether the fixed SPF model had general applicability in thin-airfoil stall problems.

Figure 12 shows the $C_L-AOA$ curves obtained by different turbulence models, including SA, SST, $k-\overline{v^2}-\omega$ and SPF. $AO{A}_{stall}$ and $C_{L,max}$ were both seriously underpredicted by SA and SST. The computation using these two models quickly diverged after $AO{A}_{stall}$. The $k-\overline{v^2}-\omega$ model performed relatively better. However, the $C_L-AOA$ curve jumped after $AO{A}_{stall}$, and the computation diverged at $AOA=17°$. The SPF model still obtained the result closest to the experiment among all four turbulence models. For the SPF model, the relative error between the computed $C_{L,max}$ and the experimental result was smaller than 5.5%. The lift coefficient predicted by the SPF model rose as $AOA$ increased at large $AOAs$, which is a typical characteristic of a thin-airfoil stall.

Figure 13 shows the $C_p$ distributions on the airfoil surface obtained by the SPF, $k-\overline{v^2}-\omega$, SST and experimental data. All numerical results agreed well with the experiment at $AOA=4°$ and $AOA=5°$. At $AOA=10°$, the $k-\overline{v^2}-\omega$ model failed to predict a flat $C_p$ distribution on the upper surface, while the SST model gave a wavy $C_p$ distribution that was not observed in the experiment. The SPF model successfully captured the flat $C_p$ distribution, and the result agreed well with the experiment. Figure 14 shows the velocity profiles in the boundary layer at $AOA=4°$ given by the SPF, $k-\overline{v^2}-\omega$ and SST models and the experiment. The figure illustrates that at all stations selected, the prediction of the velocity profile given by the SPF model agreed well with the experiment, while $k-\overline{v^2}-\omega$ overestimated the velocity and SST underestimated the velocity at most locations. The accurate computation of the boundary layer velocity profile constituted a basis for the SPF model to correctly predict the flow separation associated with the airfoil stall. Considering the relatively precise prediction of $C_L$ and the details of the flow field, the fix in the SPF model had some general applicability in the prediction of the thin-airfoil stall.

3.3. Stall Prediction of the NACA0012 Airfoil

The stall characteristics of the NACA0012 airfoil are reported in this section. The freestream Mach number was $0.15$, and the Reynolds number (based on the chord length) was $8.86\times {10}^6$. All the numerical results were compared with the experiment carried out by Ladson [36]. The computation was carried out on a C-type grid that had $457$ grid points in the circumferential direction and $97$ grid points in the normal direction.

According to Gregory and Reilly [37], the NACA0012 airfoil only exhibits trailing-edge stalling characteristics at the Reynolds number $Re=8.86\times {10}^6$. The boundary layer separation point moves from the trailing edge to the leading edge, and the separation bubble grows in size as $AOA$ increases before the airfoil stall. The separation point is not fixed while the non-equilibrium turbulence is not that strong in the separated shear layer. Consequently, the flow pattern in this case differs from the previous patterns, except that the separated shear layer still presents itself.

Figure 15a,b shows the $C_L-AOA$ curve and drag polar curves obtained by the SST, SA and SPF models. The result of the $k-\overline{v^2}-\omega$ model is not presented here because the computation had difficulty converging at this case. The SST and SA models seriously underestimated $AO{A}_{stall}$ and $C_{L,max}$, while the SPF model underestimated them by only $1.5°$ and 0.06, respectively. The relative error between the computed $C_L$ and the experimental data was smaller than 3.5% when the experimental $C_L>1.0$. Because of underestimated $AO{A}_{stall}$, both the SST and SA models gave drag polar curves that dramatically deviated from the experiment. On the other hand, the drag polar curve obtained by SPF started to deviate from the experiment only at a larger $C_L$. Therefore, the SPF model predicted the force coefficient near the stall more accurately than the other two models.

Figure 16a,b shows the streamlines past the airfoil at $AOA=17°$ and $AOA=18°$ obtained by the SPF model. At both angles of attack, there was a separation bubble near the trailing edge. As the angle of attack increased from $17°$ to $18°$, the separation point moved forward, and the separation bubble grew, which is exactly the characteristic of trailing-edge stall. This indicates that the SPF model correctly captured the features of boundary layer separation in the trailing-edge stall. Consequently, the SPF model’s applicability was not limited to thin-airfoil stall problems. The model for trailing-edge stall prediction was also satisfactory.

4. Conclusions

Stall is a typical flow phenomenon for an airfoil operating at large angles of attack. To accurately predict airfoil stall, an accurate RANS model is expected to reflect the non-equilibrium characteristic of turbulence in the separated shear layer. In this paper, the SA, SST, $k-\overline{v^2}-\omega$ and SPF $k-\overline{v^2}-\omega$ models were used to compute the stall performance of the NACA64A010, NA64A006 and NACA0012 airfoils. The SA and SST models were mainly calibrated for equilibrium turbulence, and these models failed to predict the airfoil stall performance. The SPF $k-\overline{v^2}-\omega$ model supplemented by the separated shear layer fix terms captured the non-equilibrium characteristic in the shear layers successfully and predicted the stall performance accurately. The work of this paper can be summarized as follows:

(1) The SPF model was proven to be effective for the stall prediction of NACA64A006 and NACA64A010 airfoils at different Reynolds numbers, which all feature thin-airfoil stall. The relative errors between the computed $C_{L,max}$ and the experimental data at both cases were less than 5.5%.

(2) The SPF model was also applied to predict the stall of the NACA0012 airfoil, which exhibited a trailing-edge stalling feature. The relative error between the predicted $C_L$ and the experimental data was smaller than 3.5% when the lift $C_L>1.0$. The SPF model also showed good effectiveness on the trailing-edge stall problem.

The authors believe that the SPF $k-\overline{v^2}-\omega$ model can also be applied to other airfoils featuring thin-airfoil stall and trailing-edge stall if the subsonic flow is considered. Hence, the SPF $k-\overline{v^2}-\omega$ model has the potential to be applied to the verification of the stall performance and the numerical optimization of an airfoil’s stall performance.

We need to further explore the applicability of the model in the 3D cases and the 3D corrections that may be needed. This is because the stall performance of a 3D wing is also critical in aircraft design. In addition, we need to consider turbulence models that take fluid compressibility into account and are capable of stall prediction, given that many current aircraft are faster and may continue to get faster.