Assessment of a Spalart–Allmaras Model Coupled with Local Correlation Based Transition Approaches for Wind Turbine Airfoils

Abstract

This paper present recent advances in the development of local correlation based laminar–to–turbulent transition modeling relying on the Spalart–Allmaras equation. Such models are extremely important for the flow regimes involved in wind energy applications. Indeed, fully turbulent flow models are not completely reliable to predict the aerodynamic force coefficients. This is particularly significant for the wind turbine blade sections. In this paper, we focus our attention on two different transitional flow models for Reynolds–Averaged Navier–Stokes (RANS) equations. It is worth noting that this is a crucial aspect because standard RANS models assume a fully turbulent regime. Thus, our approaches couple the well–known $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ technique and $log\gamma$ equation with the Spalart–Allmaras turbulence model in order to overcome the common drawbacks of standard techniques. The effectiveness, efficiency, and robustness of the above-mentioned methods are tested and discussed by computing several flow fields developing around airfoils operating at Reynolds numbers typical of wind turbine blade sections.

1. Introduction

The aim of this work is to provide a contribution to the emerging research area related to the development of CFD techniques for the prediction of flows’ laminar–to–turbulent transition. This issue is particularly relevant in the study of the flow field developing around wind turbine airfoils. In this regard, the prediction of the flow’s laminar–to–turbulent transition is extremely important, especially when Laminar Separation Bubbles (LSBs) occur near the leading edge of the airfoil [1].

Reynolds–Averaged Navier–Stokes (RANS) approaches can be considered very appealing for the aerodynamic design of modern wind turbines due to their limited computational resource requirements. However, these methods are often not reliable since they are developed under the hypothesis of a fully turbulent regime. In order to overcome this drawback, some locally formulated transition models have been included in the RANS equations [2,3,4].

This paper focuses on Local Correlation based Transition Models (LCTM), which have demonstrated acceptable accuracy in a wide range of applications. LCTM methods integrate CFD, compatible transport equations, and empirical correlations for transition modeling purposes. In this context, it is essential to underline that transition models compatible with the RANS approach have become popular methods to deal with flows at moderate Reynolds numbers. Indeed, under these specific flow conditions, the transition from a laminar to a turbulent regime can have a significant impact, not only on the design of wind turbine blades, but also in several aerodynamic applications [5,6].

Standard LCTM approaches use two additional transport equations, one for the turbulence intermittency factor and the other for the transition onset correlation. The transition model, called $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$, provides correlations for the sake of initiating the transition process in the equations. More recently, Menter et al. [7] proposed a new version of the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ model, named the $\gamma$ model, which avoids the $\widetilde{{\mathrm{Re}}_{\theta ,t}}$ transport equation by employing an algebraic scheme.

It is important to put into evidence that the LCTM methods were initially coupled with the Shear–Stresses Transport (SST) k–$\omega$ turbulence model by developers. However, several papers concerning the the adoption of the Spalart–Allmaras (SA) turbulence model in the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ context have appeared in recent years [8,9,10,11,12,13,14]. The main reason behind this choice is due to the fact that the SA equation has produced very good results in the computation of external flows. Furthermore, the model has also shown a lower computational cost than the SST k–$\omega$ model. Despite this evident attention to the LCTM–SA approach, only one paper focused on the $\gamma$ model coupling with the SA equation was published by Liu et al. [13]. They proposed an implementation for compressible flows developed within the open–source SU2 library.

In this paper, we present recent advances in the development of laminar–to–turbulent transition modeling based on the SA equation. In particular, starting from our implementation, introduced in [11,12,15], we discuss the impact of a low Reynolds number correction for the turbulent variable. Furthermore, we also introduce and assess a new implementation of the $\gamma$–SA model derived from the contribution of Liu et al. [13]. Specifically, we discuss the impact of a positivity preserving implementation for the intermittency transport equation, confirming the key role of the numerical solution procedure adopted for LCTM techniques [16]. The reliability and effectiveness of the proposed approaches are also discussed in the paper. Several benchmark cases, spanning a wide range of Reynolds numbers, are faced.

This paper is organized as follows. Section 2 presents the governing equations considered here. Section 3 briefly describes the discretization technique. Section 4 is devoted to the presentation of the numerical results, and Section 5 contains the conclusions.

2. Governing Equations

The complete set of our flow governing equations can be written as:

$$\begin{array}{c}\nabla ·\mathbf{u}=0,\hfill \\ \frac{\partial \mathbf{u}}{\partial t}+\nabla ·\left(\mathbf{u}\otimes \mathbf{u}\right)+\nabla p-\nabla ·\left(\left(\nu +{\nu }_t\right)\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\right)=0,\hfill \\ \frac{\partial \tilde{\nu }}{\partial t}+\nabla ·\left(\mathbf{u}\tilde{\nu }\right)-{\mathrm{P}}_{\tilde{\nu }}+{\mathrm{D}}_{\tilde{\nu }}-\frac{c_{b2}}{\sigma }\nabla \tilde{\nu }·\nabla \tilde{\nu }-\frac{1}{\sigma }\nabla ·\left(\left(\nu +\tilde{\nu }\right)\nabla \tilde{\nu }\right)=0,\hfill \end{array}$$

(1)

where $\mathbf{u}$ is the velocity vector, $p=P/\rho$ is the pressure divided by the density, d is the distance form the nearest wall, and $\nu$ is the kinematic viscosity. The turbulent viscosity, ${\nu }_t$, is computed according to the $\tilde{\nu }$ variable as:

$${\nu }_t=f_{v1}\tilde{\nu },$$

(2)

the following closure functions are introduced to close the $\tilde{\nu }$ equation (Equation (1)):

$$\begin{array}{cccc}\hfill f_{v1}& =\frac{{\chi }^3}{\left({\chi }^3+c_{v1}^3\right)},\hfill & \hfill f_{v2}& =1-\frac{\chi }{\left(1+\chi f_{v1}\right)},\hfill \\ \hfill g& =r+c_{w2}\left(r^6-r\right),\hfill & \hfill f_w& =g{\left[\frac{1+c_{w3}^6}{g^6+c_{w3}^6}\right]}^{\frac{1}{6}},\hfill \\ \hfill \tilde{S}& =\left[\mathsf{\Omega }+\min \left(0,S-\mathsf{\Omega }\right)\right]+\frac{\tilde{\nu }}{k^2{d}^2}{f}_{v2},\hfill & \hfill r& =\left\{\begin{array}{cc}{r}_{max}\hspace{1em}\hfill & \frac{\tilde{\nu }}{\tilde{S}{k}^2{d}^2}<0\hfill \\ min\left(\frac{\tilde{\nu }}{\tilde{S}{k}^2{d}^2},r_{max}\right)\hspace{1em}\hfill & \frac{\tilde{\nu }}{\tilde{S}{k}^2{d}^2}\ge 0\hfill \end{array}\right.,\hfill \end{array}$$

(3)

where $\chi =\tilde{\nu }/\nu$ is the dimensionless turbulent variable, $\mathsf{\Omega }=\sqrt{2\mathbf{W}:\mathbf{W}}$ is the vorticity tensor module, $S=\sqrt{2\mathbf{D}:\mathbf{D}}$ is the strain rate tensor module, and $\tilde{S}$ is a function of both the vorticity magnitude $\mathsf{\Omega }$ and $\tilde{\nu }$. Finally, to complete the SA model, the following (standard) closure constants are adopted:

$$\begin{array}{cccccc}\hfill c_{b1}& =0.1355,\hfill & \hfill \hspace{1em}{c}_{b2}& =0.622,\hfill & \hfill \hspace{1em}{c}_{v1}& =7.1,\hfill \end{array}$$

(4)

$$\begin{array}{cccc}\hfill \sigma & =2/3,\hfill & \hfill \hspace{1em}{c}_{w1}& =\frac{c_{b1}}{k^2}+\frac{(1+c_{b2})}{\sigma },\hfill \end{array}$$

(5)

$$\begin{array}{cccccc}\hfill c_{w2}& =0.3,\hfill & \hfill \hspace{1em}{c}_{w3}& =2,\hfill & \hfill \hspace{1em}k& =0.41.\hfill \end{array}$$

(6)

2.1. $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA Transition Model

In this approach, we use two transport equations to model the transition:

$$\begin{array}{c}\frac{\partial \gamma }{\partial t}+\nabla ·\left(\mathbf{u}\gamma \right)={\mathrm{P}}_{\gamma }-{\mathrm{D}}_{\gamma }+\nabla ·\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_f}\right)\nabla \gamma \right],\hfill \\ \frac{\partial \widetilde{{\mathrm{Re}}_{\theta ,t}}}{\partial t}+\nabla ·\left(\mathbf{u}\widetilde{{\mathrm{Re}}_{\theta ,t}}\right)={\mathrm{P}}_{\theta ,t}+\nabla ·\left[{\sigma }_{\theta ,t}\left(\nu +{\nu }_t\right)\nabla \widetilde{{\mathrm{Re}}_{\theta ,t}}\right].\hfill \end{array}$$

(7)

The source terms in the $\gamma$ equation are defined as:

$$\begin{array}{c}{\mathrm{P}}_{\gamma }=c_{a1}S{\left[\gamma {\mathrm{F}}_{\mathrm{onset}}\right]}^{0.5}\left(1-c_{e1}\gamma \right){\mathrm{F}}_{\mathrm{length}},\hfill \\ {\mathrm{D}}_{\gamma }=c_{a2}\mathsf{\Omega }\gamma {\mathrm{F}}_{\mathrm{turb}}\left(c_{e2}\gamma -1\right),\hfill \end{array}$$

(8)

and in ${\mathrm{P}}_{\gamma }$, the term ${\mathrm{F}}_{\mathrm{onset}}$ is computed as:

$${\mathrm{F}}_{\mathrm{onset}}=\max \left({\mathrm{F}}_{\mathrm{onset},2}-{\mathrm{F}}_{\mathrm{onset},3},0\right)$$

(9)

with:

$$\begin{array}{c}{\mathrm{F}}_{\mathrm{onset},2}=\min \left(\max \left({\mathrm{F}}_{\mathrm{onset},1},{\mathrm{F}}_{\mathrm{onset},1}^4\right),4\right),\hfill \\ {\mathrm{F}}_{\mathrm{onset},3}=\max \left(2-{\left(\frac{{\mathrm{R}}_{\mathrm{T}}}{2.5}\right)}^3,0\right),\hfill \\ {\mathrm{F}}_{\mathrm{onset},1}=\frac{{\mathrm{Re}}_{\nu }}{2.193{\mathrm{Re}}_{\theta ,c}}.\hfill \end{array}$$

(10)

In Equation (10), the terms ${\mathrm{Re}}_{\nu }$ and ${\mathrm{R}}_{\mathrm{T}}$ are obtained as follows:

$${\mathrm{Re}}_{\nu }=\frac{S{d}^2}{\nu },\hspace{1em}{\mathrm{R}}_{\mathrm{T}}=\frac{{\nu }_t}{\nu }.$$

(11)

The aspects concerning the terms ${\mathrm{F}}_{\mathrm{length}}$ and ${\mathrm{Re}}_{\theta ,c}$ are described in Section 2.1.1. As regards ${\mathrm{D}}_{\gamma }$, the coefficient ${\mathrm{F}}_{\mathrm{turb}}$ is defined as:

$${\mathrm{F}}_{\mathrm{turb}}=exp{\left(-\frac{{\mathrm{R}}_{\mathrm{T}}}{4}\right)}^4.$$

(12)

For the source terms in the transport equation for $\widetilde{{\mathrm{Re}}_{\theta ,t}}$, ${\mathrm{P}}_{\theta ,t}$, the following equation is adopted:

$${\mathrm{P}}_{\theta ,t}=\frac{c_{\theta ,t}}{\mathrm{T}}\left({\mathrm{Re}}_{\theta ,t}-\widetilde{{\mathrm{Re}}_{\theta ,t}}\right)\left(1-{\mathrm{F}}_{\theta ,t}\right).$$

(13)

In Equation (13), the last term ${\mathrm{F}}_{\theta ,t}$ is defined as:

$${\mathrm{F}}_{\theta ,t}=\min \left(\max \left(exp{\left(-\frac{{\left|\mathbf{u}\right|}^2}{375\mathsf{\Omega }\nu \widetilde{{\mathrm{Re}}_{\theta ,t}}}\right)}^4,1-{\left(\frac{\gamma -1/c_{e2}}{1-1/c_{e2}}\right)}^2\right),1.0\right).$$

(14)

The term $\mathrm{T}$, which appears in ${\mathrm{P}}_{\theta ,t}$, is defined as follows: $500\nu /{\left|\mathbf{u}\right|}^2$. Finally, the handling of ${\mathrm{Re}}_{\theta ,t}$ in Equation (13) is further discussed, together with the ${\mathrm{F}}_{\mathrm{length}}$ coefficient, in the next subsection.

For the transition model, the following closure constants are adopted:

$$\begin{array}{cccccc}\hfill c_{a1}& =2.0,\hfill & \hfill c_{a2}& =0.06,\hfill & \hfill c_{e1}& =1.0,\hfill \end{array}$$

(15)

$$\begin{array}{cccccc}\hfill c_{e2}& =50,\hfill & \hfill c_{\theta ,t}& =0.03,\hfill & \hfill {\sigma }_f& =1.0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\sigma }_{\theta ,t}& =2.0.\hfill \end{array}$$

(17)

2.1.1. Empirical Correlations in the Model

Like other $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ approaches available in the literature, the present model contains three empirical correlations needed to compute ${\mathrm{Re}}_{\theta ,t}$, ${\mathrm{Re}}_{\theta ,c}$, and ${\mathrm{F}}_{\mathrm{length}}$. In this paper, the correlation developed by Menter et al. [2] for ${\mathrm{Re}}_{\theta ,t}$ is adopted:

$${\mathrm{Re}}_{\theta ,t}=\left\{\begin{array}{c}\left(1173.51-589.428·Tu+0.2196/T{u}^2\right)F\left({\lambda }_{\theta }\right)\hspace{1em}Tu\le 1.3\hfill \\ 331.5{\left(Tu-0.5668\right)}^{-0.671}F\left({\lambda }_{\theta }\right)\hspace{1em}Tu>1.3\hfill \end{array}\right.,$$

(18)

$$F\left({\lambda }_{\theta }\right)=\left\{\begin{array}{c}1+\left[12.986{\lambda }_{\theta }+123.66{\lambda }_{\theta }^2+405.689{\lambda }_{\theta }^3\right]exp\left(-{\left(\frac{Tu}{1.5}\right)}^{1.5}\right)\hspace{1em}{\lambda }_{\theta }\le 0\hfill \\ 1+0.275\left[1-exp\left(-35{\lambda }_{\theta }\right)\right]exp\left(-\frac{Tu}{0.5}\right)\hspace{1em}{\lambda }_{\theta }>0\hfill \end{array}\right..$$

(19)

It is important to note that the correlations in Equations (18) and (19) contain the turbulence intensity $Tu$. In the framework of the k–$\omega$ model, $Tu$ can be computed using the solution to the k equation. Here, we adopt the approach introduced in [8]. Specifically, $Tu=T{u}_{\infty }$ is set for all the points of the flow field.

Moreover, ${\mathrm{Re}}_{\theta ,t}$ is computed by iterating on the value of ${\theta }_t$, since ${\mathrm{Re}}_{\theta ,t}$ is a function of ${\theta }_t$ itself because of the presence of ${\lambda }_{\theta }$. Differently, for ${\mathrm{Re}}_{\theta ,c}$ and ${\mathrm{F}}_{\mathrm{length}}$, we use the correlations introduced by Malan et al. [17]:

$${\mathrm{Re}}_{\theta ,c}=\min \left(0.615\widetilde{{\mathrm{Re}}_{\theta ,t}}+61.5,\widetilde{{\mathrm{Re}}_{\theta ,t}}\right),$$

(20)

$${\mathrm{F}}_{\mathrm{length}}=\min \left(exp\left(7.168-0.01173\widetilde{{\mathrm{Re}}_{\theta ,t}}\right)+0.5,300\right).$$

(21)

2.1.2. $\tilde{\nu }$ Equation Coupling with the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ Model

The production and destruction terms that appear in the $\tilde{\nu }$ transport equation are suitably defined as follows:

$$\begin{array}{c}{\mathrm{P}}_{\tilde{\nu }}={\gamma }_{\mathrm{eff}}{c}_{b1}\tilde{S}\tilde{\nu },\hfill \\ {\mathrm{D}}_{\tilde{\nu }}=c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2.\hfill \end{array}$$

(22)

The term ${\gamma }_{\mathrm{eff}}$ in Equation (22) is devoted to modeling the separation–induced transition, and it is defined as follows:

$${\gamma }_{\mathrm{eff}}=\max \left(\gamma ,{\gamma }_{\mathrm{sep}}\right)$$

(23)

with:

$${\gamma }_{\mathrm{sep}}=\min \left(2.0·\max \left[0,\left(\frac{{\mathrm{Re}}_{\nu }}{3.235{\mathrm{Re}}_{\theta ,c}}\right)-1\right]{\mathrm{F}}_{reattach},2.0\right){\mathrm{F}}_{\theta ,t}$$

(24)

and:

$${\mathrm{F}}_{reattach}=exp{\left(-\frac{{\mathrm{R}}_{\mathrm{T}}}{20}\right)}^4.$$

(25)

2.2. $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA20 Transition Model

This version of the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ transition model exactly replicates the framework described above. The only difference, compared to the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA model, is the $c_{w2}$ coefficient, which controls the strength of the near–wall destruction term in the $\tilde{\nu }$ equation [18]. In this case, $c_{w2}$ is replaced by the function:

$$c_{w2LRe}=c_{w4}+\frac{c_{w5}}{{\left(\frac{\chi }{40}+1\right)}^2}$$

(26)

with $c_{w4}=0.21$ and $c_{w4}=1.5$.

This approach was presented by Spalart and Garbaruk [19], to reduce the destruction term in the SA equation in the near–wall region to increase the skin friction. This technique is tested here, for the first time, to predict transitional flows.

2.3. log $\gamma$–SA Transition Model

As already introduced in Section 1, we also consider a second strategy to include transitional effects in the Spalart–Allmaras turbulence model. Specifically, in our approach, only one equation for $\tilde{\gamma }=log\gamma$, Equation (27), is solved to account for the laminar–to–turbulent transition. This method, formerly introduced by Ilinca and Pelletier [20], is very attractive because it allows guaranteeing the positivity of the intermittency. Moreover, the logarithmic distribution of a variable ensures a much smoother behavior with respect to the use of the primitive variable itself. This is a focal point of our approach, because the baseline model formulation exhibits a blow–up of the computations.

The transport equation for $\tilde{\gamma }$ is obtained by changing $\gamma =exp\left(\tilde{\gamma }\right)$ in the intermittency equation reported in Liu et al. [13]:

$$\frac{\partial \tilde{\gamma }}{\partial t}+\nabla ·\left(\mathbf{u}\tilde{\gamma }\right)={\mathrm{P}}_{\tilde{\gamma }}-{\mathrm{D}}_{\tilde{\gamma }}+\nabla ·\left[\left(\nu +{\nu }_t\right)\nabla \tilde{\gamma }\right]+\left(\nu +{\nu }_t\right){\left|\nabla \tilde{\gamma }\right|}^2.$$

(27)

Production and destruction terms for $\tilde{\gamma }$ equation have the following expressions:

$${\mathrm{P}}_{\tilde{\gamma }}={\mathrm{F}}_{\mathrm{length}}S\left(1-e^{\tilde{\gamma }}\right){\mathrm{F}}_{\mathrm{onset}}\hspace{1em}{\mathrm{D}}_{\tilde{\gamma }}=c_{a2}\mathsf{\Omega }{\mathrm{F}}_{\mathrm{turb}}\left(c_{e2}{e}^{\tilde{\gamma }}-1\right){\mathrm{F}}_{\mathrm{onset}}.$$

(28)

The ${\mathrm{F}}_{\mathrm{onset}}$ term in Equation (28) is defined as in the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ framework. Differently, ${\mathrm{F}}_{\mathrm{onset},2}$ and ${\mathrm{F}}_{\mathrm{onset},3}$ are re–calculated as follows:

$$\begin{array}{c}{\mathrm{F}}_{\mathrm{onset},2}=\min \left({\mathrm{F}}_{\mathrm{onset},1},2.0\right),\hfill \\ {\mathrm{F}}_{\mathrm{onset},3}=\max \left(1-{\left(\frac{{\mathrm{R}}_{\mathrm{T}}}{3.5}\right)}^3,0\right).\hfill \end{array}$$

(29)

As regards ${\mathrm{F}}_{\mathrm{turb}}$, the following relation is adopted for this LCTM model:

$${\mathrm{F}}_{\mathrm{turb}}=exp{\left(-\frac{{\mathrm{R}}_{\mathrm{T}}}{2}\right)}^4.$$

(30)

We want to also remark that the constants $c_{a2}$ and $c_{e2}$ assume the same values reported in Equation (15).

The parameter ${\mathrm{F}}_{\mathrm{length}}$ is fixed at $0.5$ as in Lit et al. [13] in order to achieve a smooth growth of the intermittency. On the other hand, ${\mathrm{Re}}_{\theta ,c}$ is defined here as:

$${\mathrm{Re}}_{\theta ,c}=f_{local}\left(T{u}_l\right)f_{far}\left(T{u}_{\infty }\right),$$

(31)

where:

$$\begin{array}{c}{f}_{local}\left(T{u}_l\right)=803.73{\left(T{u}_l+0.6067\right)}^{-1.027},\hfill \\ f_{far}\left(T{u}_{\infty }\right)=-3.162T{u}_{\infty }^2-0.4565T{u}_{\infty }+1.7.\hfill \end{array}$$

(32)

Looking at Equations (31) and (32), it is very easy to note that both the local turbulence intensity and the far-field turbulence information are involved. This is an interesting feature of the $log\gamma$–SA correlation based transition model. Indeed, a local turbulence intensity controlling term, $f_{local}$, is adopted to reflect the effect of local turbulence intensity variation on ${\mathrm{Re}}_{\theta ,c}$. The formulation of Cakmakcioglu [21] is used for the local turbulence intensity, $T{u}_l$.

The integration of the $log\gamma$ model and the Spalart–Allmaras turbulence model is handled using a scheme similar to the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$ one. The production term of the modified turbulent viscosity, ${\mathrm{P}}_{\tilde{\nu }}$, is treated exactly as in Equation (22). However, for the coefficient ${\gamma }_{sep}$, we use a correlation similar to the one suggested by Liu et al. [13]:

$${\gamma }_{sep}=min\left(8.0max\left[0.0,\left(\frac{{\mathrm{Re}}_{\nu }}{2{\mathrm{Re}}_{\theta ,c}}\right)-1\right],{\gamma }_{lim}\right)$$

(33)

In the above presented equation, ${\gamma }_{lim}$ is a user–defined parameter imposed equal to $2.5$. As regards the destruction term appearing in the turbulence equation, ${\mathrm{D}}_{\tilde{\nu }}$, we use the same approach proposed by Liu et al. [13]:

$$D_{\tilde{\nu }}=min\left(max\left(\gamma ,0.5\right),1.0\right)\left[c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2\right].$$

(34)

2.4. Boundary Conditions

The boundary conditions for $\tilde{\nu }$ are standard: ${\tilde{\nu }}_{\infty }=3\nu$ at the free stream and $\tilde{\nu }=0$ at the wall, whereas the boundary condition for $\gamma$ at the wall is a zero normal gradient. At the inlet, the intermittency is set to one. The boundary condition for $\widetilde{{\mathrm{Re}}_{\theta ,t}}$ at the wall is zero flux; differently, a fixed value condition is adopted at the inlet for $\widetilde{{\mathrm{Re}}_{\theta ,t}}$. This value is calculated from the specific empirical correlation based on the inlet turbulence intensity.

All the grids used in the following are able to guarantee a viscous sub–layer scaled first cell height, $y^{+}$, of approximately one [3]. The value of $y^{+}$ is estimated as $y^{+}=\frac{u_{\tau }}{\nu }{y}_c$, where $u_{\tau }=\sqrt{{\tau }_w/\rho }$ is the friction velocity, ${\tau }_w$ is the viscous stress component measured at the wall, and $y_c$ is the height of the cells next to the wall.

3. Numerical Solution

The governing equations described in Section 2 are spatially discretized by adopting an unstructured colocated FVM. In particular, our solution strategy relies on the well–known OpenFOAM library [22], which is an open–source package for continuum mechanics released under the GNU Public License (GPL).

It is important to remark that numerical solutions are obtained through simpleFoam, which is the steady solver for incompressible flows available in OpenFOAM official releases. simpleFoam uses the well–established SIMPLE algorithm [23] for pressure–velocity decoupling; the Rhie–Chow correction is adopted to remove oscillations in the solutions [24].

For all the computations presented in this paper, the diffusive terms and pressure gradients were approximated with second–order accurate central schemes. The convective terms for the momentum, turbulence, and transition equations were handled with a second–order accurate linear–upwind scheme. As regards the linear solvers, a PBiCG with the DILU preconditioner was used to solve the discretized momentum, $\tilde{\nu }$, $\gamma$, and $\widetilde{{\mathrm{Re}}_{\theta ,t}}$ equations. On the other hand, a PCG with a diagonal incomplete Cholesky preconditioner was adopted for the pressure. Lastly, a local accuracy of ${10}^{-7}$ was established for the pressure, whereas other linear systems were considered as converged when the residuals reached the machine precision.

4. Results

In this section, we discuss the prediction capabilities of the different transition models of interest in this paper. Three different airfoils were identified as proper benchmarks: SD7003 at $\mathrm{Re}=6\times {10}^4$, E387 at $\mathrm{Re}=2\times {10}^5$, and S809 at $\mathrm{Re}=2\times {10}^6$. The first two airfoils were not designed for wind energy applications, so it is worth justifying this choice. The authors participated in the GARTEUR (Group for Aeronautical Research and Technology in Europe) Action Group: Improving the Modeling of Laminar Separation Bubbles (IMoLa). This paper represents a part of our effort in this research field. The preliminary agreement of the project members is on the basis of two fundamental points: (i) the reliability of literature data available for the specific flow problem; (ii) the flow problem difficulty. Thus, the airfoils mentioned above are suitable choices to test the analyzed transitional models, and the reliability of the reference results used for benchmarking is ensured. Moreover, it was our intention to test our approaches against a wide range of $\mathrm{Re}$ numbers, i.e., a very important issue in wind turbine blade aerodynamics.

4.1. The $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA and $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA20 Models’ Results

4.1.1. SD7003 Airfoil

The Selig–Donovan (SD) 7003 airfoil was investigated at $\mathrm{Re}=6\times {10}^4$. A C–topology grid with 768 (96 in the wake) cells in the stream–wise direction and 176 cells in the normal–to–the–wall direction were used. This grid was generated at “Centro Italiano Ricerche Aerospaziali” (CIRA)research center and extensively tested in the computation of this kind of flow field [25]. It was made available to the authors as participants of the previously mentioned GARTEUR project, AG59, which started in 2019.

Here, we compare the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA and $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA20 transition models with the LES and k–$\omega$ results of Catalano and Tognaccini [25,26]. At $\alpha =4^{°}$ the pressure coefficient, $c_p=2\left(p-p_{\infty }\right)/\rho u_{\infty }^2$, the obtained distributions put in evidence a good agreement between our RANS results and the literature ones; see Figure 1a. A negligible impact of the Spalart and Garbaruk correction can be noted. The skin friction coefficient behavior, $c_f=2{\tau }_w/\rho u_{\infty }^2$, at the same angle of attack is shown in Figure 1b. It is quite evident that the RANS separation and reattachment points are in good agreement with the LES ones. The $C_f$ peak is sufficiently reproduced in the downstream (transitional part) of the bubble. Evident discrepancies remain past the turbulent reattachment. In this case, the Spalart and Garbaruk correction slightly improves the overall model behavior.

In Figure 2 pressure and skin friction coefficient distributions at $\alpha =6^{°}$ are reported. At this angle of attack, it is possible to note a situation similar to the previous case. Results for $\alpha =8^{°}$ are depicted in Figure 3; for the skin friction coefficient, we obtained a similar behavior in comparison with the other angles of attack presented above. On the other hand, the $c_p$ coefficient prediction gets clear benefits from the $c_{w2LRe}$ expression. Indeed, it is very easy to note that the pressure wiggles, which appear near the transition point, are completely suppressed; see Figure 3a. Nevertheless, $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA–20 agrees rather well with the k–$\omega$ solutions published by Catalano and Tognaccini [25], rather than the literature LES data. At $\alpha ={10}^{°}$, the impact of the Spalart and Garbaruk correction is even more marked, see Figure 4. The $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA computation produces the pressure distribution of a stalled airfoil, whereas the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA20 model results are in satisfactory agreement with the k–$\omega$ data.

The force coefficients are compared, in Figure 5, to three sets of experimental data and to the numerical results obtained by Catalano and Tognaccini [25]. The lift coefficient, $C_L=2L^{\prime }/\rho u_{\infty }^{2}c$, and the drag coefficient, $C_D=2D^{\prime }/\rho u_{\infty }^{2}c$, are considered for comparison. The measurements depicted in Figure 5 were obtained from Selig et al. [27] at the University of Princeton in 1989, from Selig et al. [28] at the University of Illinois in 1996, and from Ol et al. [29] at the Horizontal Fee–Surface Wind Tunnel (HFWT) of the Air Force Research Laboratory in 2005. At low–medium angles of attack, the computed force coefficients are well in accordance with the experimental data. Airfoil stall is differently predicted by the models presented here. In particular, $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA shows the stall at $\alpha ={10}^{°}$. On the contrary, $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA20 seems to approach k–$\omega$ of Catalano and Tognaccini [25], but unsteadiness appears in the solution for $\alpha >{10}^{°}$.

4.1.2. S809 Airfoil

The airfoil S809 is a laminar–flow airfoil with a $21\%$ thickness and was specially designed for horizontal axis wind turbines [30]. Several wind tunnel tests were conducted for this airfoil at Colorado State University, Ohio State University, and Delft University of Technology. Moreover, in the literature, several papers treating CFD simulations of S809 [31,32,33,34] have been published. This airfoil has a crucial feature related to the very thin LSB appearing on its surface. Thus, the computation of the flow field developing over this airfoil can be considered a reliable and challenging benchmark for the last generation of laminar–turbulent transition models.

In the following, a chord line length based Reynolds number equal to $2\times {10}^6$ was considered. The solutions presented here were computed using an O–shaped structured grid having about $7.5\times {10}^5$ cells. The inflow/outflow boundaries were placed at about 18 c from the airfoil. The first cell height was arranged in order to obtain $O\left(y^{+}\right)\simeq 1$.

The pressure distribution for some angles of attack are reported in Figure 6, Figure 7 and Figure 8. The obtained results underline the good agreement with the experimental data published by Somers [30]. The force coefficient curves are given in Figure 9a. These data put in evidence the effectiveness of the transition model. Indeed, the force coefficient curve prediction capability of our approach is certainly better than the fully turbulent SA technique of Xu et al. [35]. Furthermore, transition locations are satisfactory well captured by the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA model, as reported in Figure 9b.

In our experience, we noted that the impact of the $c_{w2LRe}$ correlation was negligible for this high $\mathrm{Re}$ benchmark problem; hence, for this reason, we do not present the results obtained from the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA20 model.

4.2. log $\gamma$–SA Model Results

In this subsection, we discuss the $log\gamma$–SA transition model. It is worth noting that in this paper, we developed and tested an implementation for incompressible flows within the OpenFOAM library. In such a context, we had to introduce the $log\gamma$ formulation, following Ilinca and Pelletier [20], in order to overcome the blow–up of the computations. Moreover, we also verified that, for $\alpha \ge 6^{°}$, the Spalart and Garbaruk correction improves the model behavior for medium $\mathrm{Re}$ numbers. On the other hand, for low $\mathrm{Re}$ configurations, this approach does not produce relevant effects; see Figure 10a.

Liu et al. [13] did not report the implementation issues discussed above. However, it is important to note that they handled compressible flows through the SU2 code.

In Figure 10, we represent the $c_p$ distribution computed on SD7003’s surface at $\mathrm{Re}=6\times {10}^4$. The same computational grid and numerical settings described in Section 4.1.1 were employed. The results obtained from the $log\gamma$–SA model were not satisfactory for $\alpha =4^{°}$ (see Figure 10a) and even nonphysical for $\alpha =8^{°}$ (see Figure 10b). A similar behavior was obtained for other angles of attack, but not reported for compactness. Therefore, we cannot consider the $log\gamma$–SA approach sufficiently reliable for low $\mathrm{Re}$ applications.

The Eppler 387 airfoil at $\mathrm{Re}=2\times {10}^5$ was studied for comparison with Liu et al. [13]. The solutions reported in the following were computed using a structured C–shaped computational grid having about $2.22\times {10}^5$ cells. It was built by the authors and already successfully tested for the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA model in [11,12]. The $c_p$ data presented here (Figure 11 and Figure 12) are globally very consistent with the experimental findings of McGhee et al. [36]. A particular mention has to be made for $\alpha =8^{°}$. McGhee et al. [36] observed an LSB near to the leading edge at $\alpha =8.5^{°}$, while a natural transition in the boundary layer was evidenced at $\alpha =8^{°}$. This behavior may contribute to explaining the discrepancy between the numerical and experimental data, already discussed in D’Alessandro et al. [11]. The force coefficients, shown Figure 13, confirm the very good consistency between $log\gamma$–SA and experimental data prior to the stall region. However, it is important to highlight that the stall was not correctly predicted by $log\gamma$–SA, different from what happens for the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA technique. Furthermore, the $\gamma$–SA results published in Liu et al. [13] reveal a dissimilar behavior if compared with the $log\gamma$–SA implemented here.

Lastly, the S809 airfoil outputs are not shown here since the $log\gamma$–SA technique proved to be unable to correctly predict the transitional effect on this high $\mathrm{Re}$ problem.

5. Conclusions

In this work, we present transitional RANS models for the SA equation. Flow regimes relevant for wind energy applications are considered. In particular, starting from a previous implementation of the authors, we carefully assess the impact of a low–$\mathrm{Re}$ number correction for the $c_{w2}$ term. Furthermore, a new implementation for a one–equation ($\gamma$) transition model for incompressible flows is presented.

The obtained results put in evidence that the $\gamma$–$\widetilde{{\mathrm{Re}}_{\theta ,t}}$–SA model is rather mature for predicting laminar separation bubbles. In fact, good results are highlighted from low–$\mathrm{Re}$ numbers up to high–$\mathrm{Re}$ numbers. In this framework, the approach introduced by Spalart and Garbaruk [19] is a key ingredient for improving the performance of the flow model at low–$\mathrm{Re}$ numbers. On the other hand, $log\gamma$–SA shows fairly encouraging results. For this specific model, we proposed to switch to the $log\gamma$ formulation, since unstable computations were evidenced in the standard configuration. It is worth noting that for the SD7003 airfoil, we obtain a satisfactory feedback. However, the $c_p$ distribution is still far from the reference data. Only for the E387 airfoil, we obtain results prior to the stall region. On the contrary, for S809 at $\mathrm{Re}=2\times {10}^6$, this approach is not able to detect the laminar–to–turbulent transition. Finally, the one–equation LCTM approach is very attractive because it allows predicting complex transitional flows with only a further transport equations. However, in our experience, it can be considered to be only as a premature stage for the SA equation. For this reason, future work will be devoted to its further development.