Articles | Volume 6, issue 3
Wind Energ. Sci., 6, 715–736, 2021
Wind Energ. Sci., 6, 715–736, 2021

Research article 25 May 2021

Research article | 25 May 2021

A simplified model for transition prediction applicable to wind-turbine rotors

A simplified model for transition prediction applicable to wind-turbine rotors
Thales Fava1, Mikaela Lokatt1, Niels Sørensen2, Frederik Zahle2, Ardeshir Hanifi1, and Dan Henningson1 Thales Fava et al.
  • 1Department of Mechanics, Linné Flow Centre, SeRC, KTH Royal Institute of Technology, Stockholm, Sweden
  • 2Department of Wind Energy, Technical University of Denmark, Risø Campus, Roskilde, Denmark

Correspondence: Thales Fava (


This work aims to develop a simple framework for transition prediction over wind-turbine blades, including effects of the blade rotation and spanwise velocity without requiring fully three-dimensional simulations. The framework is based on a set of boundary-layer equations (BLEs) and parabolized stability equations (PSEs), including rotation effects. An important element of the developed BL method is the modeling of the spanwise velocity at the boundary-layer edge. The two analyzed wind-turbine geometries correspond to a constant airfoil and the DTU 10-MW Reference Wind Turbine blades. The BL model allows an accurate prediction of the chordwise velocity profiles. Further, for regions not too close to the stagnation point and root of the blade, profiles of the spanwise velocity agree with those from Reynolds-averaged Navier–Stokes (RANS) simulations. The model also allows predicting inflectional velocity profiles for lower radial positions, which may allow crossflow transition. Transition prediction is performed at several radial positions through an “envelope-of-envelopes” methodology. The results are compared with the eN method of Drela and Giles, implemented in the EllipSys3D RANS code. The RANS transition locations closely agree with those from the PSE analysis of a 2D mean flow without rotation. These results also agree with those from the developed model for cases with low 3D and rotation effects, such as at higher radial positions and geometries with strong adverse pressure gradients where 2D Tollmien–Schlichting (TS) waves are dominant. However, the RANS and PSE 2D models predict a later transition in the regions where 3D and rotation effects are non-negligible. The developed method, which accounts for these effects, predicted earlier transition onsets in this region (e.g., 19 % earlier than RANS at 26 % of the radius for the constant-airfoil geometry) and shows that transition may occur via highly oblique modes. These modes differ from 2D TS waves and appear in locations with inflectional spanwise velocity. However, except close to the root of the blade, crossflow transition is unlikely since the crossflow velocity is too low. At higher radial positions, where 3D and rotation effects are weaker and the adverse pressure gradient is more significant, modes with small wave angles (close to 2D) are found to be dominant. Finally, it is observed that an increase in the rotation speed modifies the spanwise velocity and increases the Coriolis and centrifugal forces, shifting the transition location closer to the leading edge. This work highlights the importance of considering the blade rotation and the three-dimensional flow generated by that in transition prediction, especially in the inner part of the blade.

1 Introduction

In wind-turbine design, accurate determination of aerodynamic loads is of importance as they are related to properties, such as performance and structural loads. Since aerodynamic loads can be influenced by the boundary-layer character, an accurate determination of the transition location can be significant to obtain a successful wind-turbine design. This has long been recognized by aerodynamicists, and significant efforts have been devoted to the development of transition models.

There are several transition models available (for a review see, e.g., Saric et al.2003; Langtry et al.2006; Pasquale et al.2009; Colonia et al.2017). Some of these are based on the transport equations, such as the γ (Colonia et al.2017) and γ-RẽΘ equation models (Menter et al.2006; Langtry et al.2006; Sørensen2009; Menter et al.2015; Langtry et al.2015); other ones rely on stability analysis, such as the eN method (Smith and Gamberoni1956; van Ingen1956). These models are compatible with modern RANS solvers. In particular, the models of natural and bypass transition coupled with RANS solvers have shown good agreement with experiments on wind turbines (Özçakmak et al.2020). The γ-RẽΘ has also been used for prediction of transition dominated by crossflow instability (Guerrero et al.2018). More accessible measurement techniques such as ground-based thermographic imaging (Reichstein et al.2019) have offered further data for the development, calibration, and comparison of transition models. The methods mentioned above can provide transition predictions at a relatively low computational cost, being common in engineering applications. While their accuracy has been validated for a number of two- and three-dimensional flows, further knowledge about their performance for rotating wind-turbine blades would be beneficial.

There are also more advanced transition-prediction methods, such as those based on direct numerical simulations (DNSs) and parabolized stability equations (PSEs) (Bertolotti et al.1992; Simen and Dallmann1992), which can provide accurate transition prediction in three-dimensional flows. DNSs aim at exactly resolving the flow field, and they can thus provide detailed information about velocity fluctuations within the boundary layer, based on which results about transition and turbulence characteristics can be derived. At this moment, only a few studies of the transition process on wind-turbine blades using high-resolution simulations are available (Jing et al.2020). The DNS approach for transition prediction provides accurate results, but it implies a high computational cost. With the current available computational power, simulations at Reynolds numbers corresponding to those of real wind turbines are not possible. The PSE analysis has a much lower computational cost compared to DNSs (Özçakmak et al.2020), but it provides more accurate transition predictions than the RANS approach with an algebraic-integral or transport model. However, there are limitations in the linear PSE approach, which are the inability to predict (i) transition in strongly non-parallel flows with rapid variation in the streamwise direction; (ii) transition in strongly three-dimensional flows; and (iii) transition caused by global instability, as in the case of strong separation bubbles.

In two-dimensional flow fields, the waves causing instability are typically of the Tollmien–Schlichting (TS) type (Tollmien1929; Schlichting1933), whereas in three-dimensional flow fields, waves of crossflow type are also common (Saric et al.2003). The former is more prone in wings with small sweep angles and very weak or adverse chordwise pressure gradients, while the latter generally takes place for large sweep angles and favorable chordwise pressure gradients. Borodulin et al. (2019) showed a good agreement between linear stability results and experiments for TS waves developing over a swept wing. There were similarities between the TS waves found experimentally and those for the Blasius boundary layer, such as the shape of the eigenfunctions and phase speed. However, the waves observed over the swept wing could propagate at a broader range of angles relative to the inviscid streamline, being more unstable at propagation angles between 25 and 70. Unlike the TS instability, crossflow instability has an inviscid origin, caused by the inflection of the crossflow velocity profile (Saric et al.2003). Unstable crossflow modes can be triggered by noise or even microscopic surface roughness (Bippes1999; Gaponenko et al.2002). The crossflow instability can manifest as stationary vortices in environments with low turbulence intensity and as traveling modes in cases with high turbulence intensity/low surface roughness. These waves can propagate at a narrower range of angles compared to TS waves and are more unstable for directions nearly perpendicular to the inviscid flow direction.

The present work aims to develop a simple model for transition prediction applicable for wind-turbine blades and to understand the effects of blade rotation on the boundary-layer flow and its stability. Firstly, a model to compute the boundary-layer profiles over the wind-turbine blades is developed. This model is based on the quasi-three-dimensional boundary-layer equations (BLEs) and accounts for effects of the blade rotation and the three-dimensional outer flow. A technique to obtain an approximation for the spanwise velocity is also provided, such that the only required inputs are the chordwise distribution of pressure or streamwise velocity and the blade geometry. Secondly, the eN method is employed to predict the transition locations. The N factors are obtained using an existing PSE code (Hanifi et al.1994; Hein et al.1994) to which rotation effects are added. The developed framework is applied to two different full-scale wind-turbine geometries and the results are compared with the mean-flow and transition data from EllipSys3D RANS simulations (Michelsen1992, 1994; Sørensen1994). Transition prediction within this solver is obtained through the semiempirical eN method of Drela and Giles (Drela and Giles1987; Özçakmak et al.2020). This transition model does not account for effects of the blade rotation or the three-dimensional flow. The PSE results may also indicate accuracy of the RANS prediction. Finally, effects of the rotation speed and spanwise velocity on the transition location are analyzed, and the suitability of XFOIL (Drela1989) data as the input to the developed model is assessed.

2 Boundary-layer model

This section describes the boundary-layer (BL) model developed in this work.

2.1 Coordinate system

The coordinate system of the BL model is illustrated in Fig. 1. The blade rotates around a vertical axis at a constant angular velocity Ω, and the coordinate system is fixed to the blade. Therefore, centrifugal and Coriolis forces need to be included in the fluid-dynamic equations (Kundu et al.2016). The first coordinate direction x1 follows the wing contour along a circular arc with radius r0, the second coordinate direction x2 is perpendicular to the x1 direction in the plane tangent to the wing surface, and the third coordinate direction x3 is defined to be in the direction normal to the surface. Hence, x1, x2, and x3 describe an orthogonal, curvilinear coordinate system. The error committed by assuming that the x1 and x2 directions are respectively the chordwise and spanwise directions is low. That is because the chord-to-radius ratio and the sweep angle are small in the analyzed wind-turbine blades. For instance, the angle between the x2 and spanwise directions oscillates between 1 and 4.

Figure 1Coordinate system on the wind-turbine blade. Ω is the rotation speed, u the mean-flow velocity vector projected in the x1x2 plane, k=(α,β,0) the wave vector, and Ψ the perturbation propagation angle relative to the outer streamline.


2.1.1 Boundary-layer equations

There are several integral formulations of the boundary-layer equations (BLEs) (Du and Selig2000; Dumitrescu and Cardos2011; Drela2013; Garcia et al.2014). However, a differential formulation is expected to be more accurate than its integral counterpart because the latter requires closure relations which are found through empirical relations (van Garrel2004). For this reason, a differential formulation is selected in the present case. When expressed in the coordinate system described in Sect. 2.1, the differential form of the BLEs can be written as (Warsi1999)


In these equations, cp, γ, κ, μ, M, Re, and Pr denote specific heat capacity at constant pressure, ratio of specific heats, thermal conductivity, dynamic viscosity, Mach number, Reynolds number based on a reference length l0, and Prandtl number, respectively. Moreover, ρ, p, and T denote density, pressure, and temperature, whereas u and Ω represent velocity and rotation, respectively. hi denotes the Lamé coefficients, where hi2=gii and gij is the metric tensor. Note that since the coordinate system is orthogonal gij=0 for ji. The subscripts 1, 2, and 3 indicate components in the respective x1, x2, and x3 directions. r is the radial position.

In the BL model, the chordwise curvature of the wing model is neglected, while the radial curvature is considered. Thus, the metric vector becomes

(5) h 1 = x 2 + r 0 r 0 , h 2 = 1 , h 3 = 1 .

Since the code is intended for analysis of laminar flows, turbulent fluctuations and statistics need not be considered. In order to obtain a well-conditioned system whose solution is compatible with the subsequent PSE analysis, the terms in the system of Eqs. (1) to (4) are normalized by the reference quantities given in Table 1. The value of l0 is set to c0, the chord of the airfoil at the radial position r0, where the analysis is performed.

Table 1Reference values.  denotes freestream values.

Download Print Version | Download XLSX

2.1.2 Approximations of the spanwise derivatives

As they stand, the BL equations are dependent on all three coordinate directions so that their numerical solution requires a full volume discretization. A three-dimensional discretization can result in a solution procedure that is costly in terms of computational capacity and CPU time. By employing approximate models for the derivative terms in the x2 direction, instead of exact expressions, one can obtain a quasi-three-dimensional model requiring discretization in the x1 and x3 directions only. The reduced dimension of the discretization typically results in significant savings in computational cost and meshing effort. Furthermore, a judicious selection of the model for the x2 derivative can provide accurate mean flows. These beneficial properties lead a quasi-three-dimensional model to be employed in the present work.

Similarity solutions for rotating flows suggest that the velocity in the x1 direction can be assumed to depend on the x2 coordinate linearly (Greenspan1968; Hernandez2011). This approximation is employed in the present work, together with the further assumption that the velocity in the x2 direction, pressure, and temperature does not depend on x2. Thus,

(6) u 1 = u 1 0 x 2 + r 0 r 0 , u 2 = u 2 0 , p = p 0 , T = T 0 .

The subscript 0 denotes evaluation at the radial location r0. This choice can result in a momentum imbalance in the x2 direction at the boundary-layer edge, as pointed by Sturdza (2003) for swept-wing flows. Sturdza argued that the imbalance could be compensated for by defining an additional source term A that accounts for the momentum difference. The extra source term is then multiplied by a blending function f(x3) and added to the right-hand side of the spanwise momentum equation (Eq. 3). A is found by considering momentum balance at the boundary-layer edge. With the current approximation of spanwise derivatives and curvature terms, A becomes

(7) A = ρ u 1 e u 2 e x 1 - ρ u 1 e 2 r 0 - ρ - 2 Ω 3 u 1 + Ω 2 r 0 ,

where the subscript “e” denotes evaluation at the boundary-layer edge. The blending function is selected to linearly depend on the wall-normal distance inside the boundary layer, i.e.,

(8) f x 3 = x 3 x 3 e .

2.1.3 Discretization of BLEs

The spanwise approximations described in Sect. 2.1.2 make the system of the BLEs (Eqs. 1 to 4) include only derivatives in the x1 and x3 directions. The derivatives in the x3 direction are evaluated using a second-order central finite-difference scheme, whereas the derivatives in the x1 direction are evaluated using a second-order backward Euler finite-difference scheme.

The BLEs can be expressed as

(9) A 1 Φ + A 2 Φ x 3 + A 3 2 Φ x 3 2 + A 4 Φ x 1 = A 5 ,

where Φ=(u1, u2, T)T denotes the vector of primary variables. The density is calculated from the temperature and pressure using the equation of state and the BL approximation of pressure being constant inside the boundary layer. The components of the matrices A1, A2, A3, A4, and A5 are found by collecting terms in Eqs. (1) to (4).

The solution is computed by space marching in the x1 direction. Uniform boundary conditions are assumed at the inflow. The attachment-line equations (Cebeci1999) are solved at the first inflow node, since the BLEs are ill-conditioned when u1 is equal to zero. Because of the boundary-layer singularity (Goldstein1948), the system of equations can become strongly ill-conditioned if flow separation is encountered. However, the present code is intended to be used for transition prediction, and separation within a laminar-flow region typically causes transition. Therefore, the separation point can be taken as a reasonable approximation of the transition location, and the issue is circumvented.

2.2 Edge velocity model

The velocity in the x2 direction at the boundary-layer edge is required as input to the quasi-three-dimensional BL model. In order to avoid the necessity of a costly simulation to obtain it, a model for u2e is devised with inspiration from the conical-wing approximation (Cebeci1999; Sturdza2003). An approximation for u2e is obtained by combining the Euler equation in the x2 direction with an approximation for the variation of the pressure coefficient in this direction. The Euler equation in the x2 direction can be written as (Warsi1999)

(10) ρ u 1 h 1 u 2 x 1 + u 2 h 2 u 2 x 2 + u 3 h 3 u 2 x 3 + 1 h 1 h 2 h 2 x 1 u 1 u 2 - h 1 x 2 u 1 2 + 1 h 2 h 3 h 2 x 3 u 2 u 3 - h 3 x 2 u 3 2 = - 1 h 2 p x 2 + F rot 2 ,


(11) F rot 2 = ρ 2 u 3 Ω 1 - 2 u 1 Ω 3 - Ω 2 x 3 - Ω 3 x 2 Ω 3 + Ω 1 x 2 - Ω 2 x 1 Ω 1 .

We assume that u2h2u2x20, based on the fact that the flow and the variations in the x2 direction have a small magnitude. A second hypothesis is that u3h3u2x30, built on the evidence that the flow and variations in the normal direction at the boundary-layer edge are small. Since u3≈0 and Ω1≈0, the term 2u3Ω1 is neglected in Eq. (11). However, the terms u3h3u2x3 and 2u3Ω1 may be relevant close to the stagnation point because u3||u|| and Ω1||Ω||. Therefore, Eq. (10) should be valid only after a slightly downstream distance from the stagnation point. Moving all terms except the one containing u2x1 to the right-hand side, dividing both sides of the equation by ρu1h1, and including the scale factors given by Eq. (5) yield

(12) u 2 x 1 = h 1 ρ u 1 - p x 2 + F rot 2 + ρ u 1 2 h 1 x 2 .

All terms on the right-hand side are known except for the x2 pressure gradient. An approximation for this term can be found by rewriting the definition of the pressure coefficient with the reference speed equal to the rotational one, i.e.,

(13) p = C p 1 2 ρ Ω r 0 2 + p ,

and assuming that

(14) C p = C p 0 r 2 r 0 2 α α 0 ,

where Cp0 is the pressure coefficient at the radial position r0 and r=x2+r0. Equation (14) models the variation in Cp due to the change in the reference velocity with r, as well as a first-order variation in Cp due to the change in the angle of attack α. The latter is defined as

(15) α = tan - 1 w Ω r 0 + θ x 2 ,

with w and θ representing the incoming-flow velocity and the geometric twist angle, respectively. Note that Eq. (14) is singular for α0=0 and may not be very accurate for small values of α0. Therefore, some other approximations may be more suitable for these cases. With inspiration from the conical-wing approximation (Cebeci1999; Sturdza2003), Cp0 is assumed to be constant along conical lines. These lines as well as other parameters related to the conical-wing approximation are illustrated in Fig. 2.

Figure 2Conical parameters. O and A are the center of rotation and the cone apex, respectively. Lines of constant β1 are the conical lines.


With this assumption, the derivative of Cp0 in the x2 direction can be related to its derivative in the x1 direction by

(16) C p 0 x 2 = - tan β 1 + β 0 C p 0 x 1 .

The angles β1 and β0 are defined as

(17) β 1 = sin - 1 x 1 c - x 1 r 1 , β 0 = sin - 1 x 1 c - x 1 r 0 ,

where x1c denotes the x1 coordinate of point C, where the line connecting the center of rotation O and the cone apex A intersects the arc with radius r0. These assumptions lead to an expression for the pressure derivative, given by

(18) C p x 2 = - tan β 1 + β 0 C p 0 x 1 r 2 r 0 2 α α 0 + C p 0 α α 0 2 r r 0 2 + r 2 r 0 2 1 α 0 α x 2 .

Inserting Eqs. (13), (14), and (18) in Eq. (12) provides an expression that can be integrated along x1 to obtain the distribution of u2e in this direction. However, it is necessary to obtain an approximation for u2e at the initial point of integration. In order to do that, we use as inspiration the swept-wing approximation (Cebeci1999) and assume that u2e can be approximated by the velocity over a conical line (see Fig. 2). This approximation yields

(19) u 2 e = 2 Ω r 0 - u 1 e tan β 1 + β 0 ,

where r0 is a reference velocity. However, Eq. (19) is not very accurate if u1e is small, as is the case near the attachment line. Thus, it is advisable to start the integration at a position x10 downstream of the attachment line, where u1e has a value that is comparable to the freestream velocity. An approximate initial value for u2e at x10 can be found from

(20) u 2 e x 1 0 = 2 Ω r 0 - u 1 e x 1 0 x 1 c - x 1 0 r 0 + r 1 r 0 r 1 .
3 PSEs

The coordinate system employed in the PSE analysis is the one in Fig. 1. The PSEs are derived from the continuity, momentum, energy, and state equations (Hanifi et al.1994; Kundu et al.2016), as shown in Eqs. (21) to (24). Because of the complexity of performing a full three-dimensional analysis, periodicity is assumed in the x2 direction. Moreover, rotation terms are added to the momentum equations.


where λ=-23μ denotes the second viscosity coefficient under the Stokes hypothesis. The quantities in these equations have been normalized with the reference values given in Table 1.

The flow can be decomposed as

(27) q x 1 , x 3 , t = q x 1 , x 3 + ϵ q ̃ x 1 , x 3 , t ,

where t denotes time, q=(u1, u2, u3, T, ρ)T. Here, pressure is eliminated using the equation of state. The bar denotes the mean-flow variables from the BL model or RANS, tilde denotes the perturbation quantities, and ϵ≪1 (Hanifi et al.1994; Hein et al.1994). The perturbation part has the form

(28) q ̃ x 1 , x 3 , t = q ^ x 1 , x 3 e i Θ ,

where q^(x1, x3) denotes the slowly varying part of the perturbation, i denotes the imaginary unit, and Θ is

(29) Θ = x 0 x 1 α ( x ) d x + β x 2 - ω t ,

where α and β are the wavenumber in the x1 and x2 directions, respectively, whereas ω denotes the temporal angular frequency of the disturbance. x0 is the chordwise coordinate of the initial point of analysis. Including these relations in Eqs. (21) to (24), assuming that the variation in the x1 direction is weak compared to the variation in the x3 one (there is a scale of 1/Re between them), neglecting terms of order ϵ2, and collecting the terms, we obtain a system of the form

(30) B 1 q ^ + B 2 q ^ x 3 + B 3 2 q ^ x 3 2 + B 4 q ^ x 1 = 0 .

In addition, the following normalization condition is used

(31) 0 q ^ * q ^ x 3 d x 3 = 0 ,

where the superscript * denotes the complex conjugate (Hanifi et al.1994). The following boundary conditions are employed

(32) u ^ 1 = u ^ 2 = u ^ 3 = T ^ = 0 , for x 3 = 0 , u ^ 1 , u ^ 2 , u ^ 3 , T ^ 0 , for x 3 .

Notice that the far-field condition u^30 can be replaced by ρ^0. The derivatives in the x3 direction are computed with a fourth-order compact finite-difference scheme, whereas the derivatives in the x1 direction are computed with a second-order compact finite-difference scheme. Given initial values of α and β, the growth of the disturbances along x1 is evaluated by marching Eq. (30) in the x1 direction.

In the eN method, transition location is predicted based on the amplification of disturbances presented by the N factors computed as

(33) N = ln A / A 0 = x I x σ ( x ) d x ,

where A is the amplitude of the perturbations (A0=A(x0)), xI the location where the perturbation first starts to grow, and σ the growth rate of the perturbation kinetic energy E defined as (Hanifi et al.1994)

(34) σ = 1 h 1 - I m ( α ) + R e 1 E E x 1 , E = 0 ρ u ^ 1 2 + u ^ 2 2 + u ^ 3 2 d x 3 .

Here, consistent with the PSE framework, we use an “envelope-of-envelopes” approach, meaning that transition is predicted based on the envelope of the amplification curves computed for fixed values of ω and β (see, e.g., Arnal and Casalis2000).

Figure 3Wind-turbine blades with radial sections of analysis. The surface is colored with a normalized measure of the axial position of the mesh point. The radial coordinate r is given in meters. R is the radius of the wind-turbine rotor.


4 Results

The results of the proposed approach are compared to those from the EllipSys3D RANS code. This solver is based on the incompressible Navier–Stokes equations and employs a block-structured, finite-volume discretization, including a second-order upwind scheme for the discretization of convective terms and a central difference scheme for the discretization of the viscous ones. Turbulence is modeled using the shear stress transport (SST) kω turbulence model (Menter1993), and the transition prediction is performed using an eN method (Drela and Giles1987) combined with a model for the turbulence intermittency factor γ (Özçakmak et al.2020). The intermittency function is defined as

(35) γ = 1 - exp - x - x tr 2 U e , tr ν 2 n ^ σ , for x x tr ,

where x is the chordwise position (measured from the stagnation line), xtr is the chordwise position of the transition onset, ν is the kinematic viscosity, σ is the spot propagation rate, n^ is the nondimensional spot formation rate, and Ue,tr is the edge velocity at the chordwise position of the transition onset (Mayle1999). For laminar flow, i.e., x<xtr, γ=0, and for fully turbulent flow, γ=1.

4.1 Test cases

Two different full-scale wind-turbine rotors are investigated. Both have three blades, and their geometries are illustrated in Fig. 3. The shaded colors show a normalized measure of the axial position of each mesh point on the blade surface. The first geometry (Geometry 1) has a tapered and twisted blade with a symmetric NACA 63-018 airfoil profile along its entire span. It was mainly designed to allow the investigation of the accuracy of the conical-wing-based edge velocity model when applied to a geometry respecting its geometrical assumptions. The second geometry (Geometry 2) corresponds to the blade of the DTU 10-MW Reference Wind Turbine (Bak et al.2012). It has a tapered and twisted blade with spanwise-varying cross-sectional properties. This enables the evaluation of our quasi-three-dimensional model when applied to a general wind-turbine blade geometry. It is assumed that the flows over the three blades are similar so that it is sufficient to analyze one blade. We focus on the suction side of the blade since transition often occurs earlier there. Attachment-line transition is not expected to occur as the attachment-line Reynolds number R=41 and 15 for geometries 1 and 2, respectively, where R=(uRlesinϕtanϕ/(2ν))1/2, u is the incoming infinite velocity, Rle is the curvature radius of the leading edge, and ϕ is the sweep angle. This is well below the threshold of 250 for contamination (Poll1978).

The main parameters of the two cases are given in Table 2. Both were computed using a temperature of 287.5 K, density of 1.225 kg m−3, dynamic viscosity of 1.784×10-5 kg m−1 s−1, ratio of specific heats of 1.4, and gas constant of 287 J kg−1 K−1. The meshes used for the RANS computations of geometries 1 and 2 have 15.5×106 nodes, of which 118×103 are surface ones. The boundary layer is discretized with approximately 50 nodes in the wall-normal direction. The corresponding meshes for the BL and PSE models have 200 and 500 points in this direction, respectively. This level of discretization provided spatially converged results for test cases. However, a lower number of grid points could be used for increased performance when computing the envelope of N factors with the PSEs.

Table 2Physical parameters of the wind turbines.

Download Print Version | Download XLSX

For the benefit of the reader, the abbreviations of the methods used in the following sections are summarized in Table 3.

Table 3Abbreviations of the employed methods.

Download Print Version | Download XLSX

4.2 Pressure distributions

The pressure distributions from RANS and XFOIL are shown in Fig. 4. Close agreement is obtained for the middle and outer radial locations of Geometry 1. For Geometry 2 and the inner radial location of Geometry 1, XFOIL results indicate a less severe pressure drop along the airfoil, although RANS and XFOIL pressure gradients are close to each other for the initial chordwise extent of the airfoils. For Geometry 1 at r0/R=0.26 and Geometry 2 at r0/R=0.89, XFOIL results also indicate small separation bubbles at x1≈0.45, which are not present in RANS distributions. A possible source of those differences is the mismatch between the angles of attack (AoAs) of XFOIL and RANS. The XFOIL computations are for an AoA calculated based on the inflow velocity and that generated by the blade rotation, which may differ from the actual AoA in the RANS simulation. Moreover, XFOIL Cp distributions were obtained for a two-dimensional section of the wing, without considering its spanwise variation and the three-dimensionality of the flow present in the RANS results. Those effects are particularly important for Geometry 1 at r0/R=0.26.

Figure 4Comparison between XFOIL and RANS pressure distributions for the suction side of the airfoils of geometries 1 and 2 at three radial positions.


4.3 Spanwise edge velocity

Here, we compare the chordwise distributions of spanwise velocity at the edge of the boundary layer u2e obtained with RANS simulations and the edge velocity model (EVM). The analyses are performed at three radial locations r0 in the inner (r0/R=0.26 and 0.40), middle (r0/R=0.58), and outer (r0/R=0.89) parts of the blade, where R is the radius of the rotor. The inner section for Geometry 2 (r0/R=0.40) is chosen after the location of the maximum chord at r0/R=0.30.

Figure 5a, c, and e present the results for Geometry 1. The spanwise velocity is of the order of 1 % of the freestream velocity, except close to the stagnation point, where it can reach higher values. EVMR and RANS results agree for the middle and outer radial locations after 10 % of the chord. The differences between EVMX and RANS results are also small for these locations. The small overestimation of u2e of the EVMX method compared to RANS/EVMR is related to the smaller flow acceleration predicted by XFOIL compared to its RANS counterpart (see Eq. 12). The differences between the EVM and RANS results are larger at the inner radial position and close to the stagnation point. The reason is that the approximation for the spanwise pressure gradient given by Eq. (16) is more accurate at large radii and chordwise positions. This approximation relies on the assumption of Cp0 being constant over conical lines, which may not be respected at the mentioned locations due to the strong variation of the geometry in the radial direction and the flow three-dimensionality.

Figure 5Spanwise edge velocity.


The results for Geometry 2 are presented in Fig 5b, d, and f. At the inner radial location, r0/R=0.40, EVMR and EVMX results indicate a higher spanwise velocity than RANS, similarly to Geometry 1. In previous analysis of Geometry 2 (Zahle et al.2014), a region of three-dimensional flow radially pumped from the root to r0/R=0.36 was observed. Moreover, a separation bubble is also present from the root to almost r0/R=0.40 (Horcas et al.2017). These factors increase the flow three-dimensionality at the inner radial part of the blade, making it more difficult for the quasi-three-dimensional BL model to capture the flow features correctly. However, the agreement between EVM and RANS results improves with r0/R and x1. This is particularly true at r0/R=0.58 and 0.89 after 15 % of the chord. The differences between EVM and RANS velocity distributions were expected to be higher for Geometry 2 because the spanwise variation of the airfoil spurs changes in the Cp along conical lines. The higher spanwise velocity of Geometry 1, especially at the inner radial location, associated with inflectional spanwise velocity profiles, as shown in the next section, indicates a higher potential for crossflow instability. These results suggest that the edge velocity model can provide a reliable approximation for u2e for radial positions not too close to the root of the blade and stagnation point. The results are expected to be more accurate for geometries respecting the assumptions of the model and generating a less three-dimensional flow, such as Geometry 1.

4.4 Velocity profiles

We present the chordwise and spanwise velocity profiles obtained with RANS simulations and the boundary-layer model as a function of the normal coordinate x3 nondimensionalized by the BL thickness δ. Two chordwise positions are analyzed for each radial location. Figure 6 presents the results for Geometry 1. The BLR, BLX, and BLR 2D profiles of chordwise velocity are in close agreement with the RANS results for all locations. They resemble the Falkner–Skan type of profiles for an accelerating flow and seem to be little affected by three-dimensionality since they agree with the BLR 2D solution. Further downstream, around x1=0.40, the flow starts to decelerate (see Fig. 4), which may allow the appearance of a viscous instability of the Tollmien–Schlichting (TS) type. These conclusions also apply to Geometry 2, whose results are shown in Fig. 7. The qualitative behavior of the chordwise velocity profiles is similar. However, the flow starts to decelerate earlier at around x1=0.30 for the inner radial position and at approximately x1=0.40 for the middle and outer radial locations. Therefore, an earlier transition may be expected for Geometry 2 at r0/R=0.40.

Figure 6Boundary-layer profiles for Geometry 1.


Figure 7Boundary-layer profiles for Geometry 2.


The spanwise velocity at the inner radial position of Geometry 1 is directed towards the root of the blade as portrayed in Fig. 6a and  b. This reverse flow supports the hypothesis of considerable three-dimensionality at radial locations closer to the root of the blade (Du and Selig2000). Although the BLR and BLX profiles of spanwise velocity are close to each other, they indicate a positive velocity (flow towards the tip of the blade), whereas the spanwise velocity profile from RANS is only positive in the near-wall region. The RANS, BLR, and BLX spanwise velocity profiles present inflection points. Therefore, they are susceptible to an inviscid instability of the crossflow type. Other cases with inflection of the spanwise velocity profile are the RANS and BLR results at r0/R=0.58 of geometries 1 and 2 (Figs. 6d and 7c, d) and the RANS results at r0/R=0.40 of Geometry 2 (Fig. 7b).

The BLR and RANS spanwise velocity profiles are in close agreement at the middle and outer radial positions of Geometry 1 as presented in Fig. 6c–f. The higher values obtained with the BLX approach in those cases are caused by the larger u2e predicted with the edge velocity model (EVMX). The same occurs at the outer radial location of Geometry 2, as shown in Fig. 7e and f, in which BLR and RANS spanwise velocity profiles agree, but the result from BLX overestimates u2e. Nonetheless, the shapes of the BLX profiles agree with that of the other methods, indicating that the mismatch is only due to the u2e values.

The BLR and BLX results for the spanwise velocity at the inner and middle radial parts of Geometry 2 (Fig. 7a–d) in general do not follow the trend of the RANS results. An exception is the BLR spanwise velocity profile at r0/R=0.58 and x1=0.25. As shown in Fig. 7a and b, the RANS profile presents an inversion of direction between 10 % and 20 % of the chord. This also occurs in a smaller extent at the inner radial position of Geometry 1 (Fig. 6a and b), where, at the near-wall region, the spanwise velocity profile presents an inversion of direction. The fact that the inversion of the spanwise velocity profile only occurs at the inner radial position of geometries 1 and 2 may confirm the three-dimensional character of the flow at smaller radii.

The effects of rotation on the spanwise velocity are investigated using the approach of Du and Selig (2000), in which the rotation speed is varied while the angle of attack is kept constant. This allows for segregating the effects of the variation of the spanwise velocity as well as Coriolis and centrifugal forces from those caused by the variation of the angle of attack. The selected rotation speeds are 5 %, 50 %, 100 %, and 150 % of that used in RANS (0.64 and 0.9 rad s−1 for geometries 1 and 2, respectively). The 5 % and 50 % cases account for the accelerating phase of the wind turbine, whereas the 150 % case is not in the normal operating range of most turbines but offers insight into how overspeed may impact transition.

Analysis of the data for Geometry 1 shows that the inviscid flow is accelerated in the x2 direction near the stagnation point due to a negative spanwise pressure gradient and the Coriolis force to a lesser extent. The dominant term of the latter is −2ρu1Ω3 in Eq. (11), pointing in the x2 direction. After roughly 10 % of the chord, where the flow reaches its maximum streamwise velocity, the spanwise pressure gradient decreases substantially. Hence, the centrifugal force with leading term ρΩ32x2 in Eq. (11) and the inertial term with ρu12 in Eq. (12) overcome the Coriolis force and accelerate the flow in the +x2 direction. For small radii, the Coriolis force tends to increase faster with the rotation speed than the centrifugal and inertial ones, impelling the flow in the x2 direction. The centrifugal and inertial forces tend to grow faster with Ω at the middle and outer parts of the blade, forcing the flow in the +x2 direction.

Figure 8 presents the BLX spanwise velocity profiles for Geometry 1. Compared to an almost translatoric situation (0.032 rad s−1), rotation tends to accelerate the flow in the x2 direction, driven by the centrifugal and inertial forces. Considering r0/R=0.58 and 0.89, the spanwise velocity increases with Ω since the centrifugal and inertial forces grow faster at larger radii. At the inner radial position, the spanwise velocity decreases when Ω increases from 0.32 to 0.96 rad s−1 because the Coriolis force grows faster than its counterparts. These velocity profiles present inflection points, indicating the potential of crossflow instability. Inflectional profiles can also be observed at the inner radial position of Geometry 2.

Figure 8Spanwise velocity profiles for Geometry 1 for several rotation speeds.


Figure 9Spanwise velocity profiles for Geometry 2 for several rotation speeds.


The boundary-layer profiles for Geometry 2 are presented in Fig. 9. The airfoils of Geometry 2 sustain negative chordwise and spanwise pressure gradients over a larger chordwise extent compared to Geometry 1. Therefore, it is not possible to decouple a region where the pressure gradient is dominant from another in which rotation effects are preponderant. This fact makes the effects of rotation less clear than in the previous geometry. However, one can still observe the trend described in the theoretical analysis. At the downstream chordwise stations, the flow accelerates with Ω in the x2 at the inner locations and in +x2 directions at the outer sections. An increase in the rotation speed tends to accelerate the flow in the x2 direction at r0/R=0.58. This fact indicates that the pressure gradient and Coriolis forces are more important than the inertial and centrifugal ones at this location. This trend remains for the downstream chordwise station since the negative spanwise pressure gradient does not vanish.

4.5 Transition prediction

The quasi-three-dimensional PSE model is applied to analyze the disturbance growth inside the boundary layer. The onset of transition is assumed to occur when the amplification factor N based on the integral disturbance energy (Hanifi et al.1994) reaches Ncrit. This state corresponds to the appearance of the first turbulent spots. Although not representative of all atmospheric conditions, it is assumed Ncrit=9 in the current work to have a larger region of laminar flow in the RANS results, allowing a more detailed comparison between the developed model and RANS. In the EllipSys3D code, when the eN method of Drela and Giles (1987) indicates that Ncrit was reached, the onset of transition is detected and the intermittency factor γ starts to grow from zero in the laminar region to one in the fully turbulent flow (Özçakmak et al.2020). As the transition location is not directly stored in RANS data, we choose to select a small value for this parameter (γ=0.01 is selected) to indicate the transition location.

The transition locations for Geometry 1 as a function of the radial position are presented in Fig. 10a. Transition is delayed as the radial position increases, which agrees with previous works that observed stabilizing effects of rotation for increasing radii (Du and Selig2000). PSER and RANS transition locations agree from r0/R=0.68 to the tip of the blade. For r0/R<0.68, PSER results indicate an earlier transition compared to RANS. This is due to the effects of the spanwise velocity and rotation, which are not considered in the EllipSys3D transition model. As shown in Sect. 4.4, the spanwise velocity reaches higher values at lower radii. Moreover, the presence of a laminar separation bubble at the inner part of the blade increases the rotation effects because the Coriolis force passes to act in the same direction of the centrifugal one. Therefore, differences between transition locations from RANS and the developed model were expected to be larger at lower radii. Another conclusion is that considering three-dimensional and rotation effects leads to the prediction of earlier transition locations. The PSER 2D transition locations, which do not consider 3D and rotational effects, are in close agreement with the RANS results, except at r0/R=0.26, where the former indicates transition slightly downstream. Concerning the PSEX results, earlier transition locations are obtained for r0/R0.58 compared to RANS and PSER. This is likely due to the higher spanwise velocity found at these locations with the PSEX method. PSEX and PSER transition locations are close to each other for lower radial positions, probably because the differences between their predicted spanwise velocity profiles are smaller.

Figure 10b presents the transition locations for Geometry 2. PSER and PSER 2D results are in close agreement. This indicates that three-dimensional effects and rotation are likely not very important for this blade. As discussed in Sect. 4.4, the pressure gradient seems to be more important than rotation effects in Geometry 2. PSER and PSER 2D present slightly downstream transition locations when compared to RANS. The PSEX transition locations are downstream of the PSER ones, possibly due to the weaker adverse pressure gradient in the Cp distributions from XFOIL. The transition delay due to increasing radius is less significant in Geometry 2, probably because of the lower influence of rotation effects.

Figure 10Transition locations.


Figure 11N-factor contours from PSER for three radial positions. The white line indicates the critical region, and the red dashed line shows the transition location.


The PSER contours of the N factor as a function of the chordwise position and propagation angle Ψ are shown in Fig. 11. Ψ is the angle between the inviscid streamline and the perturbation wave vector (see Fig. 1). The dashed red line indicates the transition location. Considering Geometry 1 in Fig. 11a–c, the region of critical N factor is displaced in the −Ψ direction, and it is less symmetrical at the inner radial location. The critical modes have Ψ=-58, −24, and −6 at r0/R=0.26, 0.58, and 0.89, respectively. The lower Ψ=-58 at r0/R=0.26 is possibly related to the stronger and inflectional spanwise velocity occurring at this location, which makes transition more susceptible to oblique and crossflow modes. Transition occurs significantly earlier at this position (x1=0.23, compared to x1=0.34 and 0.37 at r0/R=0.58 and 0.89, respectively). The PSER 2D contours of the N factor, shown in Fig. 12a–c, are more symmetrical around Ψ=0, with the critical modes having lower |Ψ| (Ψ=17, 5, and 4 for r0/R=0.26, 0.58, and 0.89, respectively). This shows that the oblique critical modes obtained in the PSER results are caused by three-dimensionality and rotation.

Figure 11d–f shows that the PSER critical regions are more elongated in the Ψ direction for Geometry 2, indicating transition susceptibility to a broader range of waves. The critical modes have Ψ=-12, −16, and −12 for r0/R=0.40, 0.58, and 0.89. These waves are less oblique than those for Geometry 1, particularly at the inner radial location. Notice that the BL profiles of spanwise velocity at this location (Fig. 7b) do not present an inflection point, making transition via lower |Ψ| modes more likely. Regarding the PSER 2D results, in Fig. 12d–f, the regions of critical N factors are more centered around Ψ=0, with the critical modes for r0/R=0.40, 0.58, and 0.89 presenting Ψ=0. This means that disregarding 3D and rotation effects in the mean flow leads to 2D critical modes for Geometry 2.

Figure 12N-factor contours from PSER 2D for three radial positions. The white line indicates the critical region, and the red dashed line shows the transition location.


Figure 13a–c presents the profiles of the perturbation of u1 velocity of the modes leading to transition in Geometry 1. The PSER and PSEX modes are in close agreement for the three radial positions, indicating that they predict the same transition mechanism. At r0/R=0.26, these modes have a single peak, located at x3/δ=0.2, associated with their high |Ψ| and the inflectional spanwise velocity (Fig. 6b). This indicates that transition may be triggered by oblique TS or crossflow modes. The PSER 2D critical mode differs from the previous ones by presenting a near-wall peak, at x3/δ=0.1, and having a second lobe for x3/δ>0.7. At r0/R=0.58, the PSER and PSEX modes approach the PSER 2D one by developing a near-wall peak, although less important than the one at x3/δ=0.2, and a second lobe for x3/δ>0.7. The PSER and PSEX modes finally converge to a 2D mode at r0/R=0.89, where they are in close agreement with the PSER 2D one. The latter is similar to a 2D TS wave, as also observed for r0/R=0.58. The appearance of near-wall peaks in the PSER and PSEX modes at r0/R=0.58 and 0.89 as well as the close agreement between these modes and the PSER 2D ones at r0/R=0.89 can be related to the amplification of 2D TS waves due to an adverse pressure gradient.

Figure 13PSE results for the mode leading to transition.


The results for Geometry 2 are presented in Fig. 13d–f. As occurs for Geometry 1, the PSER and PSEX modes agree for the three radial positions. They indicate double-peak modes, with maxima at x3/δ=0.1 and 0.2. The former has a larger or similar magnitude compared to the latter. These modes are close to the PSER 2D ones except around x3/δ=0.2, where the PSER and PSEX modes have more pronounced peaks. The presence of a peak at x3/δ=0.1 for all radial locations is related to a strong adverse pressure gradient in Geometry 2. The second peak, at x3/δ=0.2, seems to be associated with the obliqueness of the mode, having a larger amplitude for larger values of |Ψ|. A 2D TS mechanism seems to be more important in Geometry 2 because the critical modes are closer to the PSER 2D ones, and the adverse pressure gradient is stronger. However, a mechanism related to oblique TS waves, engendered by 3D and rotation effects, appears to be more important for transition in Geometry 1. This is due to its larger sweep angle and region of favorable pressure gradient. Although the crossflow velocity profiles are inflectional, the magnitude of this velocity component is very low, of the order of 0.1 % of the freestream velocity, except for the inner radial location of Geometry 1, where it reaches 3.5 %. Thus excluding Geometry 1 at r0/R=0.26, a crossflow transition mechanism is unlikely. Nevertheless, the effect of the spanwise velocity on transition cannot be neglected as it allows transition through oblique modes.

Figure 14Transition locations for several rotation speeds.


Figure 15N-factor contours from PSEX at r0/R=0.58 for several rotation speeds. The white line indicates the critical region, and the red dashed line shows the transition location.


In the following, we analyze the effects of rotation on the transition locations. Figure 14a presents the PSEX transition locations as a function of the radial position and rotation speed for Geometry 1. The displayed trend indicates that an increase in the rotation speed shifts the transition location closer to the nose. In particular, the rise in Ω from 0.32 to 0.96 rad s−1 leads to transition 37 % earlier. The case corresponding to 5 % of the RANS rotation speed (not shown) did not present any mode reaching Ncrit, further indicating the destabilizing effect of rotation. These effects occur through the Coriolis and centrifugal forces acting on the disturbances as well as through the variation of the spanwise velocity. The former seems to be preponderant since there is no significant variation in the spanwise velocity with Ω at r0/R=0.89, but transition occurs earlier regardless. There is a delay in transition for increasing radius up to r0/R=0.47, where the Coriolis force is prevalent. Further increases in radius do not significantly change the transition locations, indicating a balance between the rotation effects. The presence of a laminar separation bubble for radial positions closer to the root can make the Coriolis force act in the same direction as the centrifugal one. For higher radial positions and in the absence of separation, these two forces tend to balance each other.

Figure 14b portrays the results for Geometry 2. The increase in Ω plays a destabilizing role. This observation is supported by the fact that the case with 5 % of the RANS rotation speed (not shown) presented no mode reaching Ncrit. However, the variation of Ω does not play a role as important as for Geometry 1. For instance, transition occurs 8 % earlier on average for an increase in Ω from 0.45 to 1.35 rad s−1. The transition location moves less with the rotation speed for Geometry 2 because this blade maintains a non-negligible pressure gradient over a larger chordwise extent, overtaking rotation effects. The fact that the spanwise velocity in Geometry 2 varies more with Ω than in Geometry 1 with a smaller effect on transition corroborates this claim. Transition is delayed when increasing the radius up to r0/R=0.58, a range along which the Coriolis force is dominant. Only slight variations in transition locations occur after this radial position, pointing to a balance in the rotation effects.

Figure 16PSEX results for the mode leading to transition for several rotation speeds.


The PSEX contours of the N factor at r0/R=0.58 for geometries 1 and 2 are shown in Fig. 15. In the case of Geometry 1, as shown in Fig. 15a–c, the increase in Ω forces the critical region towards lower x1. This region lies mostly in the −Ψ half plane, meaning that the critical waves propagate towards the root of the blade. These modes present Ψ=-25, -24, and −25 for Ω=0.32, 0.64, and 0.96 rad s−1. For Geometry 2, in Fig. 15d–f, we also observe the displacement of the critical region to lower x1 with the increase in Ω. Moreover, the flat critical region extending from Ψ=-60 to 40 obtained with Ω=1.35 rad s−1 shows that the higher rotation velocity allows transition through a broader range of disturbances. The critical regions are mostly located in the −Ψ half plane, indicating stronger transition susceptibility to waves traveling to the inner blade part. The critical modes present Ψ=-16, −15, and −13 for Ω=0.45, 0.9, and 1.35 rad s−1. The analysis of the full geometry indicates that the increase in Ω reduces the critical |Ψ| in the region 0r0/Rr, where r=0.58 and 0.5 for geometries 1 and 2. For larger r, the opposite occurs, i.e., rising Ω leads to increasingly oblique critical modes.

Figure 16a–c shows the PSEX profiles of the critical modes for Geometry 1. All modes collapse at the inner radial location, indicating that Ω does not alter the transition mechanism. The inflectional spanwise velocity profiles at this location (Fig. 8a and d) seem to render the transition mechanism, through oblique modes, quite robust to changes in Ω. At r0/R=0.58 and 0.89, the modes for Ω=0.64 and 0.96 rad s−1 are in close agreement. However, the mode for Ω=0.32 rad s−1 differs from the previous ones by the presence of a near-wall peak. As already discussed, the mode shapes are closely related to their propagation angles, with higher-|Ψ| modes occurring at locations of inflectional spanwise velocity and tending to have a single peak like those at r0/R=0.26. Figure 16d–f shows the results for Geometry 2. At r0/R=0.40, the increase in Ω reduces |Ψ| and makes double-peak modes such as those for Ω=0.45 and 0.9 rad s−1 become a 2D, single-peak mode like the one for Ω=1.35 rad s−1. At r0/R=0.58, all modes collapse and present double peaks. At the outer radial location, the mode for Ω=0.45 rad s−1 is nearly 2D, and the rise in Ω increases its obliqueness (i.e., increases |Ψ|). The modes for higher Ω are in close agreement at this location. In Geometry 2, the adverse pressure gradient is more important, and transition is more susceptible to modes closer to 2D TS waves with near-wall peaks. The increase in the rotation tends to prompt these 2D modes at low radial locations, while it makes the critical modes more oblique at higher radii.

5 Conclusions

A framework for transition prediction applicable to flows over wind-turbine blades is developed. The method, which comprises a boundary-layer model and the PSE, accounts for effects of the quasi-three-dimensional flow and the blade rotation. It aims to provide more reliable transition predictions without requiring three-dimensional simulations. Using the developed method, we have analyzed the role of flow three-dimensionality and rotation on the transition onset over two geometries.

The developed method provides accurate chordwise velocity profiles and, for locations not too close to the root of the blade and stagnation point, spanwise velocity. The flow is highly three-dimensional close to the root of the blade, reducing the accuracy of a quasi-three-dimensional approach. The spanwise velocity obtained with the model better agrees with RANS for geometries respecting the conical-wing approximation. Some of the spanwise velocity profiles contain inflection points, which may allow crossflow instability, not considered in two-dimensional transition models. Rotation was shown to accelerate the flow towards the tip of the blade in the developed flow region, while the opposite occurs near the stagnation point.

Transition locations from the eN method implemented in the EllipSys3D RANS code closely agree with those from the PSE analysis of a 2D mean flow without rotation. RANS transition locations are close to those from the model developed in this work in places where 3D and rotation effects are low. This occurs for Geometry 2 and higher radial positions in Geometry 1. However, results of the RANS transition model and the 2D approach deviate from those from the new approach for locations from the root to approximately 58 % of the radius of Geometry 1, where 3D and rotation effects are important. At these locations, the combined influence of three-dimensionality and rotation leads to earlier transition onsets. These effects make the transition occur through oblique modes, which have single peaks and are not predicted with the 2D approach. The oblique modes appear in locations where the spanwise velocity profile is inflectional, raising the possibility of being related to crossflow instability. However, except for the inner radial location of Geometry 1, the magnitude of the crossflow velocity seems to be too low to trigger crossflow transition. The single-peak modes may be very oblique TS waves. For larger radial positions, the flow tends to be more two-dimensional, and the adverse pressure gradient is more important. Thus the critical modes become less oblique and develop features of 2D TS waves, such as a second peak near the wall. Finally, it is also shown that the increase in the rotation speed, through the modification of the spanwise velocity and the increase in the Coriolis and centrifugal forces, seems to shift the transition location closer to the leading edge.

In order to better understand the transition process over the rotating blades and validate the prediction of the presented approach, in-depth investigation through DNS simulations and detailed experimental works are desired.

Code and data availability

Part of the codes and data employed/developed is available upon direct request to the corresponding author.

Author contributions

TF implemented the models, performed the analysis, and wrote the final version of the manuscript. ML developed the model, obtained part of the results, and wrote the first version of the document. NS and FZ performed the RANS simulations using the EllipSys3D code. AH and DH developed the NOLOT PSE code (among other researchers), supported the analysis, provided useful discussions, and contributed with critical feedback. All authors reviewed the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


This research has been supported and funded by the Strategic Research Area initiative STandUP for Energy.

Financial support

This research has been supported by the STandUP for Energy.

Review statement

This paper was edited by Sandrine Aubrun and reviewed by three anonymous referees.


Arnal, D. and Casalis, G.: Laminar-turbulent transition prediction in three-dimensional flows, Prog. Aerosp. Sci., 36, 173–191,, 2000. a

Bak, C., Bitsche, R. D., Yde, A., and Kim, T.: Light Rotor: The 10-MW reference wind turbine, in: Proceedings of EWEA 2012 – European Wind Energy Conference & Exhibition, 16–19 April 2012, Copenhagen, Denmark, 2012. a

Bertolotti, F. P., Herbert, T., and Spalart, P. R.: Linear and nonlinear stability of the Blasius boundary layer, J. Fluid Mech., 242, 441–474, 1992. a

Bippes, H.: Basic experiments on transition in three-dimensional boundary layers dominated by crossflow instability, Prog. Aerosp. Sci., 35, 363–412,, 1999. a

Borodulin, V. I., Ivanov, A. V., Kachanov, Y. S., Mischenko, D. A., Örlü, R., Hanifi, A., and Hein, S.: Experimental and theoretical study of swept-wing boundary-layer instabilities. Three-dimensional Tollmien–Schlichting instability, Phys. Fluids, 31, 1–18,, 2019. a

Cebeci, T.: An engineering approach to the calculation of aerodynamic flows, 1st Edn., Springer, Long Beach, USA, 1999. a, b, c, d

Colonia, S., Leble, V., Steijl, R., and Barakos, G.: Assessment and Calibration of the γ-Equation Transition Model at Low Mach, AIAA J., 55, 1126–1139,, 2017. a, b

Drela, M.: XFOIL: An analysis and design system for low Reynolds number airfoils, in: Lecture Notes in Engineering, vol. 54, Springer-Verlag, Cambridge, USA, 1989. a

Drela, M.: Three-Dimensional Integral Boundary Layer Formulation for General Configurations, in: Proceedings of the 21st AIAA Computational Fluid Dynamics Conference, 24–27 June 2013, San Diego, California, 2013. a

Drela, M. and Giles, M. B.: Viscous-Inviscid Analysis of Transonic and Low Reynolds Number Airfoils, AIAA J., 25, 1347–1355,, 1987. a, b, c

Du, Z. and Selig, M. S.: The effect of rotation on the boundary layer of a wind turbine blade, Renew. Energ., 20, 167–181,, 2000. a, b, c, d

Dumitrescu, H. and Cardos, V.: Inboard Stall Delay Due to Rotation, in: Fundamental and Advanced Topics in Wind Power, chap. 3, edited by: Carriveau, R., IntechOpen, Rijeka, Croatia,, 2011. a

Gaponenko, V. R., Ivanov, A. V., Kachanov, Y. S., and Crouch, J. D.: Swept-wing boundary-layer receptivity to surface non-uniformities, J. Fluid. Mech., 461, 93–126,, 2002. a

Garcia, N. R., Sørensen, J. N., and Shen, W  Z.: A quasi-3D viscous-inviscid interaction code: Q3UIC, J. Phys.: Conf. Ser., 555, 012041,, 2014. a

Goldstein, S.: On laminar boundary layer flow near a position of separation, Q. J. Mech. Appl. Math., 1, 43–69,, 1948. a

Greenspan, H. P.: The theory of rotating fluids, 1st Edn., Cambridge University Press, Cambridge, UK, 1968. a

Guerrero, A. M., Gutierrez, L. M. G., and Remola, A. O.: On the influence of transition modeling and crossflow effects on open water propeller simulations, Ocean. Eng., 156, 101–119,, 2018. a

Hanifi, A., Henningson, D. S., Hein, S., Bertolotti, F. P., and Simen, M.: Linear nonlocal Instability Analysis – the linear NOLOT code, Technical report FFA-TN 1994-54, The Aeronautical Research Institute of Sweden, Bromma, Sweden, 1994. a, b, c, d, e, f

Hein, S., Bertolotti, F. P., Simen, M., Hanifi, A., and Henningson, D. S.: Linear nonlocal instability analysis – the linear NOLOT code, Tech. Rep. DLR-IB 223-94 A43, DLR, Göttingen, Germany, 1994. a, b

Hernandez, G. G. M.: Laminar-Turbulent transition on Wind Turbines, PhD thesis, Department of Mechanical Engineering, Technical University of Denmark, Lyngby, Denmark, 2011. a

Horcas, S. G., Debrabandere, F., Tartinville, B., Hirsch, C., and Coussement, G.: Rotor-tower interactions of DTU 10 MW reference wind turbine with a non-linear harmonic method, Wind Energy, 20, 619–636,, 2017. a

Jing, Z., Ducoin, A., and Braud, C.: Large eddy simulation of transitional boundary layer on horizontal axis wind turbine blade, PhD thesis, LHEEA, Ecole Centrale Nantes, Nantes, 2020. a

Kundu, P. K., Cohen, I. M., and Dowling, D. R.: Fluid Mechanics, 6th Edn., Academic Press, Oxford, UK, 2016. a, b

Langtry, R. B., Gola, J., and Menter, F. R.: Predicting 2D Airfoil and 3D Wind Turbine Rotor Performance using a Transition Model for General CFD Codes, in: Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, 9–12 January 2006, Reno, Nevada, 2006. a, b

Langtry, R. B., Sengupta, K., Yeh, D. T., and Dorgan, A. J.: Extending the γ-ReΘt Local Correlation based Transition Model for Crossflow Effects, in: Proceedings of the 45th AIAA Fluid Dynamics Conference, 22–26 June 2015, Dallas, Texas, 2015. a

Mayle, R. E.: A Theory for Predicting the Turbulent-Spot Production Rate, J. Turbomach., 121, 588–593,, 1999. a

Menter, F.: Zonal Two Equation kω Turbulence Models For Aerodynamic Flows, in: Proceedings of the 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, 6–9 July 1993, Orlando, Florida, 1993. a

Menter, F. R., Langtry, R., and Völker, S.: Transition Modelling for General Purpose CFD Codes, Flow Turbul. Combust., 95, 277–303,, 2006. a

Menter, F. R., Smirnov, P. E., Liu, T., and Avancha, R.: A One-Equation Local Correlation-Based Transition Model, Flow Turbul. Combust., 77, 583–619,, 2015. a

Michelsen, J.: Basis3D – a Platform for Development of Multiblock PDE Solvers, Technical report AFM 92-05, Technical University of Denmark, Lyngby, Denmark, 1992. a

Michelsen, J.: Block structured Multigrid solution of 2D and 3D elliptic PDE's, Technical Report AFM 94-06, Technical University of Denmark, Lyngby, Denmark, 1994. a

Özçakmak, O. S., Madsen, H. A., Sørensen, N. N., and Sørensen, J. N.: Laminar-turbulent transition characteristics of a 3-D wind turbine rotor blade based on experiments and computations, Wind Energ. Sci., 5, 1487–1505,, 2020. a, b, c, d, e

Pasquale, D. D., Rona, A., and Garrett, S. J.: A selective review of CFD transition models, in: Proceedings of the 39th AIAA Fluid Dynamics Conference, 22–25 June 2009, San Antonio, Texas, 2009. a

Poll, D. I. A.: Some aspects of the flow near a swept attachment line with particular reference to boundary layer transition, CoA Report 7805, Cranfield University, Cranfield, 1978. a

Reichstein, T., Schaffarczyk, A. P., Dollinger, C., Balaresque, N., Schülein, E., Jauch, C., and Fischer, A.: Investigation of Laminar–Turbulent Transition on a Rotating Wind-Turbine Blade of Multimegawatt Class with Thermography and Microphone Array, Energies, 12, 2102–2122,, 2019. a

Saric, W. S., Reed, H. L., and White, E. B.: Stability and transition of three-dimensional boundary layers, Annu. Rev. Fluid Mech., 35, 413–440,, 2003. a, b, c

Schlichting, H.: Zur Entstehung der Turbulenz bei der Plattenströmung, Nachr. Ges. Wiss., Göttingen, 181–208, 1933. a

Simen, M. and Dallmann, U.: On the instability of hypersonic flow past a pointed cone – comparison of theoretical and experimental results at Mach 8, in: Deutscher Luft- und Raumfahrtkongress/DGLR-Jahrestagung, 29 September–2 October 1992, Bremen, 1992. a

Smith, A. M. O. and Gamberoni, N.: Transition, pressure gradient and stability theory, Tech. Rep. ES 26388, Douglas Aircraft Co., El Segundo, USA, 1956. a

Sørensen, N. N.: General Purpose Flow Solver Applied to Flow over Hills, Technical Report Risø-R- 827-(EN), Risø National Laboratory, Roskilde, Denmark, 1994. a

Sørensen, N. N.: CFD modeling of Laminar-Turbulent Transition for Airfoils and Rotors using the γ-RẽΘ Model, Wind Energy, 12, 715–733,, 2009. a

Sturdza, P.: An aerodynamic design method for supersonic natural laminar flow aircraft, PhD thesis, Department of Aeronautics and Astronautics, Stanford University, Stanford, 2003. a, b, c

Tollmien, W.: Über die Entstehung der Turbulenz, English translation NACA TM 609, 1931, Nachr. Ges. Wiss., Göttingen, 21–44, 1929. a

van Garrel, A.: Integral Boundary Layer Methods for Wind Turbine Aerodynamics – A Literature Survey, Technical report ECN-C–04-004, Energy Research Centre of the Netherlands, Petten, the Netherlands, 2004.  a

van Ingen, J. L.: A suggested semiempirical method for the calculation of the boundary layer transition region, Technical report VTH-74, Faculty of Aerospace Eng., Delft University of Technology, Delft, 1956. a

Warsi, Z. U. A.: Fluid Dynamics, 2nd Edn.,CRC Press, Boca Raton, Florida, 1999. a, b

Zahle, F., Bak, C., Guntur, S., Sørensen, N. N., and Troldborg, N.: Comprehensive Aerodynamic Analysis of a 10 MW Wind Turbine Rotor Using 3D CFD, in: Proceedings of the 32nd ASME Wind Energy Symposium, 13–17 January 2014, National Harbor, Maryland, 2014. a

Short summary
This work develops a simplified framework to predict transition to turbulence on wind-turbine blades. The model is based on the boundary-layer and parabolized stability equations, including rotation and three-dimensionality effects. We show that these effects may promote transition through highly oblique Tollmien–Schlichting (TS) or crossflow modes at low radii, and they should be considered for a correct transition prediction. At high radii, transition tends to occur through 2D TS modes.