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

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 e^N 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 $\mathit{\gamma }-{\stackrel{\mathrm{̃}}{Re}}_{\mathrm{\Theta }}$ equation models (Menter et al., 2006; Langtry et al., 2006; Sørensen, 2009; Menter et al., 2015; Langtry et al., 2015); other ones rely on stability analysis, such as the e^N method (Smith and Gamberoni, 1956; van Ingen, 1956). 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 $\mathit{\gamma }-{\stackrel{\mathrm{̃}}{Re}}_{\mathrm{\Theta }}$ 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 Dallmann, 1992), 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 (Tollmien, 1929; Schlichting, 1933), 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 (Bippes, 1999; 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 e^N 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 (Michelsen, 1992, 1994; Sørensen, 1994). Transition prediction within this solver is obtained through the semiempirical e^N method of Drela and Giles (Drela and Giles, 1987; Ö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 (Drela, 1989) 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 x₁ follows the wing contour along a circular arc with radius r₀, the second coordinate direction x₂ is perpendicular to the x₁ direction in the plane tangent to the wing surface, and the third coordinate direction x₃ is defined to be in the direction normal to the surface. Hence, x₁, x₂, and x₃ describe an orthogonal, curvilinear coordinate system. The error committed by assuming that the x₁ and x₂ 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 x₂ and spanwise directions oscillates between 1 and 4^∘.

Figure 1 Coordinate system on the wind-turbine blade. Ω is the rotation speed, u_∥ the mean-flow velocity vector projected in the x₁–x₂ plane, $\mathit{k}=(\mathit{\alpha },\mathit{\beta },\mathrm{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 Selig, 2000; Dumitrescu and Cardos, 2011; Drela, 2013; 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 Garrel, 2004). 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 (Warsi, 1999)

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\left(\mathit{\rho }{h}_{\mathrm{2}}{h}_{\mathrm{3}}{u}_{\mathrm{1}}\right)+{\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\left(\mathit{\rho }{h}_{\mathrm{1}}{h}_{\mathrm{3}}{u}_{\mathrm{2}}\right)+{\displaystyle \frac{\partial }{\partial x_{\mathrm{3}}}}\left(\mathit{\rho }{h}_{\mathrm{1}}{h}_{\mathrm{2}}{u}_{\mathrm{3}}\right)=\mathrm{0},\text{(2)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }& \left({\displaystyle \frac{u_{\mathrm{1}}}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{1}}}}+{\displaystyle \frac{u_{\mathrm{2}}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{u_{\mathrm{3}}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{3}}}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}{h}_{\mathrm{2}}}}\left({\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial x_{\mathrm{2}}}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}-{\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}{u}_{\mathrm{2}}^{\mathrm{2}}\right)\right)\\ & =-{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial p}{\partial x_{\mathrm{1}}}}+{\displaystyle \frac{\mathrm{1}}{Re}}{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial }{\partial x_{\mathrm{3}}}}\left({\displaystyle \frac{\mathit{\mu }}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{3}}}}\right)\\ & +\mathit{\rho }\left(\mathrm{2}{\mathrm{\Omega }}_{\mathrm{3}}{u}_{\mathrm{2}}+{\displaystyle \frac{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}{h}_{\mathrm{1}}}}{\displaystyle \frac{\partial r^{\mathrm{2}}}{\partial x_{\mathrm{1}}}}\right),\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }& \left({\displaystyle \frac{u_{\mathrm{1}}}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}+{\displaystyle \frac{u_{\mathrm{2}}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{u_{\mathrm{3}}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}{h}_{\mathrm{2}}}}\left({\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}-{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial x_{\mathrm{2}}}}{u}_{\mathrm{1}}^{\mathrm{2}}\right)\right)\\ & =-{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial p}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{Re}}{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial }{\partial x_{\mathrm{3}}}}\left({\displaystyle \frac{\mathit{\mu }}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}}\right)\\ & +\mathit{\rho }\left(-\mathrm{2}{\mathrm{\Omega }}_{\mathrm{3}}{u}_{\mathrm{1}}+{\displaystyle \frac{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}{h}_{\mathrm{2}}}}{\displaystyle \frac{\partial r^{\mathrm{2}}}{\partial x_{\mathrm{2}}}}\right),\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }& c_{\mathrm{p}}\left({\displaystyle \frac{u_{\mathrm{1}}}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial T}{\partial x_{\mathrm{1}}}}+{\displaystyle \frac{u_{\mathrm{2}}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial T}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{u_{\mathrm{3}}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial T}{\partial x_{\mathrm{3}}}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{RePr}}{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial }{\partial x_{\mathrm{3}}}}\left({\displaystyle \frac{\mathit{\kappa }}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial T}{\partial x_{\mathrm{3}}}}\right)\\ & +\left(\stackrel{\mathrm{‾}}{\mathit{\gamma }}-\mathrm{1}\right)M^{\mathrm{2}}\left\{{\displaystyle \frac{u_{\mathrm{1}}}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial p}{\partial x_{\mathrm{1}}}}+{\displaystyle \frac{u_{\mathrm{2}}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial p}{\partial x_{\mathrm{2}}}}\right.\\ & \left.+{\displaystyle \frac{\mathit{\mu }}{Re}}\left[{\left({\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{3}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}}\right)}^{\mathrm{2}}\right]\right\}.\end{array}\end{array}$$

In these equations, c_p, $\stackrel{\mathrm{‾}}{\mathit{\gamma }}$, κ, μ, 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 l₀, and Prandtl number, respectively. Moreover, ρ, p, and T denote density, pressure, and temperature, whereas u and Ω represent velocity and rotation, respectively. h_i denotes the Lamé coefficients, where $h_i^{\mathrm{2}}=g_{ii}$ and g_ij is the metric tensor. Note that since the coordinate system is orthogonal g_ij=0 for j≠i. The subscripts 1, 2, and 3 indicate components in the respective x₁, x₂, and x₃ 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

$$\begin{array}{}\text{(5)}& h_{\mathrm{1}}={\displaystyle \frac{x_{\mathrm{2}}+r_{\mathrm{0}}}{r_{\mathrm{0}}}},h_{\mathrm{2}}=\mathrm{1},h_{\mathrm{3}}=\mathrm{1}.\end{array}$$

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 l₀ is set to c₀, the chord of the airfoil at the radial position r₀, where the analysis is performed.

Table 1 Reference values. ∞ denotes freestream values.

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 x₂ direction, instead of exact expressions, one can obtain a quasi-three-dimensional model requiring discretization in the x₁ and x₃ 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 x₂ 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 x₁ direction can be assumed to depend on the x₂ coordinate linearly (Greenspan, 1968; Hernandez, 2011). This approximation is employed in the present work, together with the further assumption that the velocity in the x₂ direction, pressure, and temperature does not depend on x₂. Thus,

$$\begin{array}{}\text{(6)}& u_{\mathrm{1}}=u_{{\mathrm{1}}_{\mathrm{0}}}{\displaystyle \frac{x_{\mathrm{2}}+r_{\mathrm{0}}}{r_{\mathrm{0}}}},u_{\mathrm{2}}=u_{{\mathrm{2}}_{\mathrm{0}}},p=p_{\mathrm{0}},T=T_{\mathrm{0}}.\end{array}$$

The subscript 0 denotes evaluation at the radial location r₀. This choice can result in a momentum imbalance in the x₂ 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(x₃) 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

$$\begin{array}{}\text{(7)}& A=\mathit{\rho }{u}_{{\mathrm{1}}_{\mathrm{e}}}{\displaystyle \frac{\partial u_{{\mathrm{2}}_{\mathrm{e}}}}{\partial x_{\mathrm{1}}}}-{\displaystyle \frac{\mathit{\rho }{u}_{{\mathrm{1}}_{\mathrm{e}}}^{\mathrm{2}}}{r_{\mathrm{0}}}}-\mathit{\rho }\left(-\mathrm{2}{\mathrm{\Omega }}_{\mathrm{3}}{u}_{\mathrm{1}}+{\mathrm{\Omega }}^{\mathrm{2}}{r}_{\mathrm{0}}\right),\end{array}$$

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.,

$$\begin{array}{}\text{(8)}& f\left(x_{\mathrm{3}}\right)={\displaystyle \frac{x_{\mathrm{3}}}{x_{{\mathrm{3}}_{\mathrm{e}}}}}.\end{array}$$

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 x₁ and x₃ directions. The derivatives in the x₃ direction are evaluated using a second-order central finite-difference scheme, whereas the derivatives in the x₁ direction are evaluated using a second-order backward Euler finite-difference scheme.

The BLEs can be expressed as

$$\begin{array}{}\text{(9)}& {\mathbf{A}}_{\mathrm{1}}\mathbf{\Phi }+{\mathbf{A}}_{\mathrm{2}}{\displaystyle \frac{\partial \mathbf{\Phi }}{\partial x_{\mathrm{3}}}}+{\mathbf{A}}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathbf{\Phi }}{\partial x_{\mathrm{3}}^{\mathrm{2}}}}+{\mathbf{A}}_{\mathrm{4}}{\displaystyle \frac{\partial \mathbf{\Phi }}{\partial x_{\mathrm{1}}}}={\mathbf{A}}_{\mathrm{5}},\end{array}$$

where Φ=(u₁, u₂, 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 A₁, A₂, A₃, A₄, and A₅ are found by collecting terms in Eqs. (1) to (4).

The solution is computed by space marching in the x₁ direction. Uniform boundary conditions are assumed at the inflow. The attachment-line equations (Cebeci, 1999) are solved at the first inflow node, since the BLEs are ill-conditioned when u₁ is equal to zero. Because of the boundary-layer singularity (Goldstein, 1948), 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 x₂ 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 $u_{{\mathrm{2}}_{\mathrm{e}}}$ is devised with inspiration from the conical-wing approximation (Cebeci, 1999; Sturdza, 2003). An approximation for $u_{{\mathrm{2}}_{\mathrm{e}}}$ is obtained by combining the Euler equation in the x₂ direction with an approximation for the variation of the pressure coefficient in this direction. The Euler equation in the x₂ direction can be written as (Warsi, 1999)

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathit{\rho }& \left[{\displaystyle \frac{u_{\mathrm{1}}}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}+{\displaystyle \frac{u_{\mathrm{2}}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{u_{\mathrm{3}}}{h_{\mathrm{3}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}}\right.\\ & +{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}{h}_{\mathrm{2}}}}\left({\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}-{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial x_{\mathrm{2}}}}{u}_{\mathrm{1}}^{\mathrm{2}}\right)\\ & \left.+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}{h}_{\mathrm{3}}}}\left({\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial x_{\mathrm{3}}}}{u}_{\mathrm{2}}{u}_{\mathrm{3}}-{\displaystyle \frac{\partial h_{\mathrm{3}}}{\partial x_{\mathrm{2}}}}{u}_{\mathrm{3}}^{\mathrm{2}}\right)\right]=-{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial p}{\partial x_{\mathrm{2}}}}+F_{{\mathrm{rot}}_{\mathrm{2}}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{F}_{{\mathrm{rot}}_{\mathrm{2}}}& =\mathit{\rho }\left[\mathrm{2}{u}_{\mathrm{3}}{\mathrm{\Omega }}_{\mathrm{1}}-\mathrm{2}{u}_{\mathrm{1}}{\mathrm{\Omega }}_{\mathrm{3}}-\left({\mathrm{\Omega }}_{\mathrm{2}}{x}_{\mathrm{3}}-{\mathrm{\Omega }}_{\mathrm{3}}{x}_{\mathrm{2}}\right){\mathrm{\Omega }}_{\mathrm{3}}\right.\\ & \left.+\left({\mathrm{\Omega }}_{\mathrm{1}}{x}_{\mathrm{2}}-{\mathrm{\Omega }}_{\mathrm{2}}{x}_{\mathrm{1}}\right){\mathrm{\Omega }}_{\mathrm{1}}\right].\end{array}\end{array}$$

We assume that $\frac{u_{\mathrm{2}}}{h_{\mathrm{2}}}\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}\approx \mathrm{0}$, based on the fact that the flow and the variations in the x₂ direction have a small magnitude. A second hypothesis is that $\frac{u_{\mathrm{3}}}{h_{\mathrm{3}}}\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}\approx \mathrm{0}$, built on the evidence that the flow and variations in the normal direction at the boundary-layer edge are small. Since u₃≈0 and Ω₁≈0, the term 2u₃Ω₁ is neglected in Eq. (11). However, the terms $\frac{u_{\mathrm{3}}}{h_{\mathrm{3}}}\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}$ and 2u₃Ω₁ may be relevant close to the stagnation point because $u_{\mathrm{3}}\approx \left|\right|\mathit{u}\left|\right|$ and ${\mathrm{\Omega }}_{\mathrm{1}}\approx \left|\right|\mathbf{\Omega }\left|\right|$. Therefore, Eq. (10) should be valid only after a slightly downstream distance from the stagnation point. Moving all terms except the one containing $\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}}$ to the right-hand side, dividing both sides of the equation by $\mathit{\rho }\frac{u_{\mathrm{1}}}{h_{\mathrm{1}}}$, and including the scale factors given by Eq. (5) yield

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}={\displaystyle \frac{h_{\mathrm{1}}}{\mathit{\rho }{u}_{\mathrm{1}}}}\left(-{\displaystyle \frac{\partial p}{\partial x_{\mathrm{2}}}}+F_{{\mathrm{rot}}_{\mathrm{2}}}+\mathit{\rho }{u}_{\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial x_{\mathrm{2}}}}\right).\end{array}$$

All terms on the right-hand side are known except for the x₂ 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.,

$$\begin{array}{}\text{(13)}& p=C_{\mathrm{p}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{\left(\mathrm{\Omega }{r}_{\mathrm{0}}\right)}^{\mathrm{2}}+p_{\mathrm{\infty }},\end{array}$$

and assuming that

$$\begin{array}{}\text{(14)}& C_{\mathrm{p}}=C_{{\mathrm{p}}_{\mathrm{0}}}{\displaystyle \frac{r^{\mathrm{2}}}{r_{\mathrm{0}}^{\mathrm{2}}}}{\displaystyle \frac{\mathit{\alpha }}{{\mathit{\alpha }}_{\mathrm{0}}}},\end{array}$$

where $C_{{\mathrm{p}}_{\mathrm{0}}}$ is the pressure coefficient at the radial position r₀ and $r=x_{\mathrm{2}}+r_{\mathrm{0}}$. Equation (14) models the variation in C_p due to the change in the reference velocity with r, as well as a first-order variation in C_p due to the change in the angle of attack α. The latter is defined as

$$\begin{array}{}\text{(15)}& \mathit{\alpha }={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{w_{\mathrm{\infty }}}{\mathrm{\Omega }{r}_{\mathrm{0}}}}\right)+\mathit{\theta }\left(x_{\mathrm{2}}\right),\end{array}$$

with w_∞ and θ representing the incoming-flow velocity and the geometric twist angle, respectively. Note that Eq. (14) is singular for α₀=0 and may not be very accurate for small values of α₀. Therefore, some other approximations may be more suitable for these cases. With inspiration from the conical-wing approximation (Cebeci, 1999; Sturdza, 2003), $C_{{\mathrm{p}}_{\mathrm{0}}}$ 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 2 Conical parameters. O and A are the center of rotation and the cone apex, respectively. Lines of constant β₁ are the conical lines.

With this assumption, the derivative of $C_{{\mathrm{p}}_{\mathrm{0}}}$ in the x₂ direction can be related to its derivative in the x₁ direction by

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial C_{{\mathrm{p}}_{\mathrm{0}}}}{\partial x_{\mathrm{2}}}}=-\tan \left({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{0}}\right){\displaystyle \frac{\partial C_{{\mathrm{p}}_{\mathrm{0}}}}{\partial x_{\mathrm{1}}}}.\end{array}$$

The angles β₁ and β₀ are defined as

$$\begin{array}{}\text{(17)}& {\mathit{\beta }}_{\mathrm{1}}={\sin }^{-\mathrm{1}}\left({\displaystyle \frac{x_{{\mathrm{1}}_{\mathrm{c}}}-x_{\mathrm{1}}}{r_{\mathrm{1}}}}\right),{\mathit{\beta }}_{\mathrm{0}}={\sin }^{-\mathrm{1}}\left({\displaystyle \frac{x_{{\mathrm{1}}_{\mathrm{c}}}-x_{\mathrm{1}}}{r_{\mathrm{0}}}}\right),\end{array}$$

where $x_{{\mathrm{1}}_{\mathrm{c}}}$ denotes the x₁ coordinate of point C, where the line connecting the center of rotation O and the cone apex A intersects the arc with radius r₀. These assumptions lead to an expression for the pressure derivative, given by

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\displaystyle \frac{\partial C_{\mathrm{p}}}{\partial x_{\mathrm{2}}}}& =-\tan \left({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{0}}\right){\displaystyle \frac{\partial C_{{\mathrm{p}}_{\mathrm{0}}}}{\partial x_{\mathrm{1}}}}{\displaystyle \frac{r^{\mathrm{2}}}{r_{\mathrm{0}}^{\mathrm{2}}}}{\displaystyle \frac{\mathit{\alpha }}{{\mathit{\alpha }}_{\mathrm{0}}}}\\ & +C_{{\mathrm{p}}_{\mathrm{0}}}\left({\displaystyle \frac{\mathit{\alpha }}{{\mathit{\alpha }}_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{2}r}{r_{\mathrm{0}}^{\mathrm{2}}}}+{\displaystyle \frac{r^{\mathrm{2}}}{r_{\mathrm{0}}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathit{\alpha }}{\partial x_{\mathrm{2}}}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(19)}& u_{{\mathrm{2}}_{\mathrm{e}}}=\left(\mathrm{2}\mathrm{\Omega }{r}_{\mathrm{0}}-u_{{\mathrm{1}}_{\mathrm{e}}}\right)\tan \left({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{0}}\right),\end{array}$$

where 2Ωr₀ is a reference velocity. However, Eq. (19) is not very accurate if $u_{{\mathrm{1}}_{\mathrm{e}}}$ is small, as is the case near the attachment line. Thus, it is advisable to start the integration at a position $x_{{\mathrm{1}}_{\mathrm{0}}}$ downstream of the attachment line, where $u_{{\mathrm{1}}_{\mathrm{e}}}$ has a value that is comparable to the freestream velocity. An approximate initial value for $u_{{\mathrm{2}}_{\mathrm{e}}}$ at $x_{{\mathrm{1}}_{\mathrm{0}}}$ can be found from

$$\begin{array}{}\text{(20)}& u_{{\mathrm{2}}_{\mathrm{e}}}\left(x_{{\mathrm{1}}_{\mathrm{0}}}\right)=\left[\mathrm{2}\mathrm{\Omega }{r}_{\mathrm{0}}-u_{{\mathrm{1}}_{\mathrm{e}}}\left(x_{{\mathrm{1}}_{\mathrm{0}}}\right)\right]\left(x_{{\mathrm{1}}_{\mathrm{c}}}-x_{{\mathrm{1}}_{\mathrm{0}}}\right){\displaystyle \frac{r_{\mathrm{0}}+r_{\mathrm{1}}}{r_{\mathrm{0}}{r}_{\mathrm{1}}}}.\end{array}$$

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 x₂ direction. Moreover, rotation terms are added to the momentum equations.

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{u}\right)=\mathrm{0},\text{(22)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }& \left[{\displaystyle \frac{\partial \mathit{u}}{\partial t}}+(\mathit{u}\cdot \mathrm{\nabla })\mathit{u}\right]=-\mathrm{\nabla }p\\ & +{\displaystyle \frac{\mathrm{1}}{Re}}\mathrm{\nabla }\left[\mathit{\lambda }\right(\mathrm{\nabla }\cdot \mathit{u}\left)\right]+{\displaystyle \frac{\mathrm{1}}{Re}}\mathrm{\nabla }\cdot \left[\mathit{\mu }\right(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T\left)\right]+{\mathit{F}}_{\mathrm{rot}},\end{array}\text{(23)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }& c_{\mathrm{p}}\left[{\displaystyle \frac{\partial T}{\partial t}}+(\mathit{u}\cdot \mathrm{\nabla })T\right]={\displaystyle \frac{\mathrm{1}}{RePr}}\mathrm{\nabla }\cdot \left(\mathit{\kappa }\mathrm{\nabla }T\right)\\ & +(\stackrel{\mathrm{‾}}{\mathit{\gamma }}-\mathrm{1})M^{\mathrm{2}}\left[{\displaystyle \frac{\partial p}{\partial t}}+(\mathit{u}\cdot \mathrm{\nabla })p+{\displaystyle \frac{\mathrm{1}}{Re}}\mathrm{\Phi }\right],\end{array}\text{(24)}& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\gamma }}{M}^{\mathrm{2}}p=\mathit{\rho }T,\text{(25)}& {\displaystyle }{\mathit{F}}_{\mathrm{rot}}=-\mathit{\rho }[\mathrm{2}\mathbf{\Omega }\times \mathit{u}+\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{x}\left)\right],\text{(26)}& {\displaystyle }\mathrm{\Phi }=\mathit{\lambda }(\mathrm{\nabla }\cdot \mathit{u}{)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\mu }{\left(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T\right)}^{\mathrm{2}},\end{array}$$

where $\mathit{\lambda }=-\frac{\mathrm{2}}{\mathrm{3}}\mathit{\mu }$ 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

$$\begin{array}{}\text{(27)}& \mathit{q}\left(x_{\mathrm{1}},x_{\mathrm{3}},t\right)=\stackrel{\mathrm{‾}}{\mathit{q}}\left(x_{\mathrm{1}},x_{\mathrm{3}}\right)+\mathit{ϵ}\stackrel{\mathrm{̃}}{\mathit{q}}\left(x_{\mathrm{1}},x_{\mathrm{3}},t\right),\end{array}$$

where t denotes time, q=(u₁, u₂, u₃, 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

$$\begin{array}{}\text{(28)}& \stackrel{\mathrm{̃}}{\mathit{q}}\left(x_{\mathrm{1}},x_{\mathrm{3}},t\right)=\widehat{\mathit{q}}\left(x_{\mathrm{1}},x_{\mathrm{3}}\right)e^{i\mathrm{\Theta }},\end{array}$$

where $\widehat{\mathit{q}}(x_{\mathrm{1}}$, x₃) denotes the slowly varying part of the perturbation, i denotes the imaginary unit, and Θ is

$$\begin{array}{}\text{(29)}& \mathrm{\Theta }=\underset{x_{\mathrm{0}}}{\overset{x_{\mathrm{1}}}{\int }}\mathit{\alpha }\left(x^{\prime }\right)d{x}^{\prime }+\mathit{\beta }{x}_{\mathrm{2}}-\mathit{\omega }t,\end{array}$$

where α and β are the wavenumber in the x₁ and x₂ directions, respectively, whereas ω denotes the temporal angular frequency of the disturbance. x₀ is the chordwise coordinate of the initial point of analysis. Including these relations in Eqs. (21) to (24), assuming that the variation in the x₁ direction is weak compared to the variation in the x₃ one (there is a scale of $\mathrm{1}/Re$ between them), neglecting terms of order ϵ², and collecting the terms, we obtain a system of the form

$$\begin{array}{}\text{(30)}& {\mathbf{B}}_{\mathrm{1}}\widehat{\mathit{q}}+{\mathbf{B}}_{\mathrm{2}}{\displaystyle \frac{\partial \widehat{\mathit{q}}}{\partial x_{\mathrm{3}}}}+{\mathbf{B}}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\widehat{\mathit{q}}}{\partial x_{\mathrm{3}}^{\mathrm{2}}}}+{\mathbf{B}}_{\mathrm{4}}{\displaystyle \frac{\partial \widehat{\mathit{q}}}{\partial x_{\mathrm{1}}}}=\mathbf{0}.\end{array}$$

In addition, the following normalization condition is used

$$\begin{array}{}\text{(31)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\widehat{\mathit{q}}}^{*}{\displaystyle \frac{\partial \widehat{\mathit{q}}}{\partial x_{\mathrm{3}}}}d{x}_{\mathrm{3}}=\mathrm{0},\end{array}$$

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

$$\begin{array}{}\text{(32)}& \left\{\begin{array}{ll}{\widehat{u}}_{\mathrm{1}}={\widehat{u}}_{\mathrm{2}}={\widehat{u}}_{\mathrm{3}}=\widehat{T}=\mathrm{0},& \mathrm{for}{x}_{\mathrm{3}}=\mathrm{0},\\ {\widehat{u}}_{\mathrm{1}},{\widehat{u}}_{\mathrm{2}},{\widehat{u}}_{\mathrm{3}},\widehat{T}\to \mathrm{0},& \mathrm{for}{x}_{\mathrm{3}}\to \mathrm{\infty }\end{array}\right..\end{array}$$

Notice that the far-field condition ${\widehat{u}}_{\mathrm{3}}\to \mathrm{0}$ can be replaced by $\widehat{\mathit{\rho }}\to \mathrm{0}$. The derivatives in the x₃ direction are computed with a fourth-order compact finite-difference scheme, whereas the derivatives in the x₁ direction are computed with a second-order compact finite-difference scheme. Given initial values of α and β, the growth of the disturbances along x₁ is evaluated by marching Eq. (30) in the x₁ direction.

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

$$\begin{array}{}\text{(33)}& N=\ln \left(A/A_{\mathrm{0}}\right)=\underset{x_{\mathrm{I}}}{\overset{x}{\int }}\mathit{\sigma }\left(x^{\prime }\right)d{x}^{\prime },\end{array}$$

where A is the amplitude of the perturbations (A₀=A(x₀)), x_I 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)

$$\begin{array}{}\text{(34)}& \begin{array}{rl}& \mathit{\sigma }={\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}}}\left[-Im\left(\mathit{\alpha }\right)+Re\left({\displaystyle \frac{\mathrm{1}}{E}}{\displaystyle \frac{\partial E}{\partial x_{\mathrm{1}}}}\right)\right],\\ & E=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\stackrel{\mathrm{‾}}{\mathit{\rho }}\left({\widehat{u}}_{\mathrm{1}}^{\mathrm{2}}+{\widehat{u}}_{\mathrm{2}}^{\mathrm{2}}+{\widehat{u}}_{\mathrm{3}}^{\mathrm{2}}\right)d{x}_{\mathrm{3}}.\end{array}\end{array}$$

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 Casalis, 2000).

Figure 3 Wind-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 (Menter, 1993), and the transition prediction is performed using an e^N method (Drela and Giles, 1987) combined with a model for the turbulence intermittency factor γ (Özçakmak et al., 2020). The intermittency function is defined as

$$\begin{array}{}\text{(35)}& \mathit{\gamma }=\mathrm{1}-\exp \left\{-{\left(x-x_{\mathrm{tr}}\right)}^{\mathrm{2}}{\left({\displaystyle \frac{U_{\mathrm{e},\mathrm{tr}}}{\mathit{\nu }}}\right)}^{\mathrm{2}}\widehat{n}\mathit{\sigma }\right\},\mathrm{for}x\ge x_{\mathrm{tr}},\end{array}$$

where x is the chordwise position (measured from the stagnation line), x_tr is the chordwise position of the transition onset, ν is the kinematic viscosity, σ is the spot propagation rate, $\widehat{n}$ is the nondimensional spot formation rate, and U_e,tr is the edge velocity at the chordwise position of the transition onset (Mayle, 1999). For laminar flow, i.e., x<x_tr, γ=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 $\stackrel{\mathrm{‾}}{R}=\mathrm{41}$ and 15 for geometries 1 and 2, respectively, where $\stackrel{\mathrm{‾}}{R}=(u_{\mathrm{\infty }}{R}_{l_{\mathrm{e}}}\sin \mathit{\varphi }\tan \mathit{\varphi }/(\mathrm{2}\mathit{\nu }){)}^{\mathrm{1}/\mathrm{2}}$, u_∞ is the incoming infinite velocity, $R_{l_{\mathrm{e}}}$ is the curvature radius of the leading edge, and ϕ is the sweep angle. This is well below the threshold of 250 for contamination (Poll, 1978).

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⁻³, dynamic viscosity of $\mathrm{1.784}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻¹ s⁻¹, ratio of specific heats of 1.4, and gas constant of 287 J kg⁻¹ K⁻¹. The meshes used for the RANS computations of geometries 1 and 2 have 15.5×10⁶ nodes, of which 118×10³ 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 2 Physical parameters of the wind turbines.

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

Table 3 Abbreviations of the employed methods.

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 $r_{\mathrm{0}}/R=\mathrm{0.26}$ and Geometry 2 at $r_{\mathrm{0}}/R=\mathrm{0.89}$, XFOIL results also indicate small separation bubbles at x₁≈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 C_p 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 $r_{\mathrm{0}}/R=\mathrm{0.26}$.

Figure 4 Comparison 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 $u_{{\mathrm{2}}_{\mathrm{e}}}$ obtained with RANS simulations and the edge velocity model (EVM). The analyses are performed at three radial locations r₀ in the inner ($r_{\mathrm{0}}/R=\mathrm{0.26}$ and 0.40), middle ($r_{\mathrm{0}}/R=\mathrm{0.58}$), and outer ($r_{\mathrm{0}}/R=\mathrm{0.89}$) parts of the blade, where R is the radius of the rotor. The inner section for Geometry 2 ($r_{\mathrm{0}}/R=\mathrm{0.40}$) is chosen after the location of the maximum chord at $r_{\mathrm{0}}/R=\mathrm{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 $u_{{\mathrm{2}}_{\mathrm{e}}}$ 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 $C_{{\mathrm{p}}_{\mathrm{0}}}$ 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 5 Spanwise edge velocity.

The results for Geometry 2 are presented in Fig 5b, d, and f. At the inner radial location, $r_{\mathrm{0}}/R=\mathrm{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 $r_{\mathrm{0}}/R=\mathrm{0.36}$ was observed. Moreover, a separation bubble is also present from the root to almost $r_{\mathrm{0}}/R=\mathrm{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 $r_{\mathrm{0}}/R$ and x₁. This is particularly true at $r_{\mathrm{0}}/R=\mathrm{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 C_p 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 $u_{{\mathrm{2}}_{\mathrm{e}}}$ 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 x₃ 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 x₁=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 x₁=0.30 for the inner radial position and at approximately x₁=0.40 for the middle and outer radial locations. Therefore, an earlier transition may be expected for Geometry 2 at $r_{\mathrm{0}}/R=\mathrm{0.40}$.

Figure 6 Boundary-layer profiles for Geometry 1.

Figure 7 Boundary-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 Selig, 2000). 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 $r_{\mathrm{0}}/R=\mathrm{0.58}$ of geometries 1 and 2 (Figs. 6d and 7c, d) and the RANS results at $r_{\mathrm{0}}/R=\mathrm{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 $u_{{\mathrm{2}}_{\mathrm{e}}}$ 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 $u_{{\mathrm{2}}_{\mathrm{e}}}$. Nonetheless, the shapes of the BLX profiles agree with that of the other methods, indicating that the mismatch is only due to the $u_{{\mathrm{2}}_{\mathrm{e}}}$ 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 $r_{\mathrm{0}}/R=\mathrm{0.58}$ and x₁=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⁻¹ 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 −x₂ 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ρu₁Ω₃ in Eq. (11), pointing in the −x₂ 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 $\mathit{\rho }{\mathrm{\Omega }}_{\mathrm{3}}^{\mathrm{2}}{x}_{\mathrm{2}}$ in Eq. (11) and the inertial term with $\mathit{\rho }{u}_{\mathrm{1}}^{\mathrm{2}}$ in Eq. (12) overcome the Coriolis force and accelerate the flow in the +x₂ 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 −x₂ 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 +x₂ direction.

Figure 8 presents the BLX spanwise velocity profiles for Geometry 1. Compared to an almost translatoric situation (0.032 rad s⁻¹), rotation tends to accelerate the flow in the x₂ direction, driven by the centrifugal and inertial forces. Considering $r_{\mathrm{0}}/R=\mathrm{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⁻¹ 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 8 Spanwise velocity profiles for Geometry 1 for several rotation speeds.

Figure 9 Spanwise 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 −x₂ at the inner locations and in +x₂ directions at the outer sections. An increase in the rotation speed tends to accelerate the flow in the −x₂ direction at $r_{\mathrm{0}}/R=\mathrm{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 N_crit. This state corresponds to the appearance of the first turbulent spots. Although not representative of all atmospheric conditions, it is assumed N_crit=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 e^N method of Drela and Giles (1987) indicates that N_crit 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 Selig, 2000). PSER and RANS transition locations agree from $r_{\mathrm{0}}/R=\mathrm{0.68}$ to the tip of the blade. For $r_{\mathrm{0}}/R<\mathrm{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 $r_{\mathrm{0}}/R=\mathrm{0.26}$, where the former indicates transition slightly downstream. Concerning the PSEX results, earlier transition locations are obtained for $r_{\mathrm{0}}/R\ge \mathrm{0.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 C_p 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 10 Transition locations.

Figure 11 N-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 $\mathrm{\Psi }=-\mathrm{58}$, −24, and −6^∘ at $r_{\mathrm{0}}/R=\mathrm{0.26}$, 0.58, and 0.89, respectively. The lower $\mathrm{\Psi }=-\mathrm{58}$^∘ at $r_{\mathrm{0}}/R=\mathrm{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 (x₁=0.23, compared to x₁=0.34 and 0.37 at $r_{\mathrm{0}}/R=\mathrm{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 $\left|\mathrm{\Psi }\right|$ (Ψ=17, 5^∘, and 4^∘ for $r_{\mathrm{0}}/R=\mathrm{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 $\mathrm{\Psi }=-\mathrm{12}$, −16, and −12^∘ for $r_{\mathrm{0}}/R=\mathrm{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 $\left|\mathrm{\Psi }\right|$ 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 $r_{\mathrm{0}}/R=\mathrm{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 12 N-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 u₁ 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 $r_{\mathrm{0}}/R=\mathrm{0.26}$, these modes have a single peak, located at $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{0.2}$, associated with their high $\left|\mathrm{\Psi }\right|$ 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 $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{0.1}$, and having a second lobe for $x_{\mathrm{3}}/\mathit{\delta }>\mathrm{0.7}$. At $r_{\mathrm{0}}/R=\mathrm{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 $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{0.2}$, and a second lobe for $x_{\mathrm{3}}/\mathit{\delta }>\mathrm{0.7}$. The PSER and PSEX modes finally converge to a 2D mode at $r_{\mathrm{0}}/R=\mathrm{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 $r_{\mathrm{0}}/R=\mathrm{0.58}$. The appearance of near-wall peaks in the PSER and PSEX modes at $r_{\mathrm{0}}/R=\mathrm{0.58}$ and 0.89 as well as the close agreement between these modes and the PSER 2D ones at $r_{\mathrm{0}}/R=\mathrm{0.89}$ can be related to the amplification of 2D TS waves due to an adverse pressure gradient.

Figure 13 PSE 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 $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{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 $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{0.2}$, where the PSER and PSEX modes have more pronounced peaks. The presence of a peak at $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{0.1}$ for all radial locations is related to a strong adverse pressure gradient in Geometry 2. The second peak, at $x_{\mathrm{3}}/\mathit{\delta }=\mathrm{0.2}$, seems to be associated with the obliqueness of the mode, having a larger amplitude for larger values of $\left|\mathrm{\Psi }\right|$. 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 $r_{\mathrm{0}}/R=\mathrm{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 14 Transition locations for several rotation speeds.

Figure 15 N-factor contours from PSEX at $r_{\mathrm{0}}/R=\mathrm{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⁻¹ leads to transition 37 % earlier. The case corresponding to 5 % of the RANS rotation speed (not shown) did not present any mode reaching N_crit, 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 $r_{\mathrm{0}}/R=\mathrm{0.89}$, but transition occurs earlier regardless. There is a delay in transition for increasing radius up to $r_{\mathrm{0}}/R=\mathrm{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 N_crit. 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⁻¹. 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 $r_{\mathrm{0}}/R=\mathrm{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 16 PSEX results for the mode leading to transition for several rotation speeds.

The PSEX contours of the N factor at $r_{\mathrm{0}}/R=\mathrm{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 x₁. This region lies mostly in the −Ψ half plane, meaning that the critical waves propagate towards the root of the blade. These modes present $\mathrm{\Psi }=-\mathrm{25}$, $-{\mathrm{24}}^{\circ }$, and −25^∘ for Ω=0.32, 0.64, and 0.96 rad s⁻¹. For Geometry 2, in Fig. 15d–f, we also observe the displacement of the critical region to lower x₁ with the increase in Ω. Moreover, the flat critical region extending from $\mathrm{\Psi }=-\mathrm{60}$ to 40^∘ obtained with Ω=1.35 rad s⁻¹ 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 $\mathrm{\Psi }=-\mathrm{16}$, −15, and −13^∘ for Ω=0.45, 0.9, and 1.35 rad s⁻¹. The analysis of the full geometry indicates that the increase in Ω reduces the critical $\left|\mathrm{\Psi }\right|$ in the region $\mathrm{0}\le r_{\mathrm{0}}/R\le r$, 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 $r_{\mathrm{0}}/R=\mathrm{0.58}$ and 0.89, the modes for Ω=0.64 and 0.96 rad s⁻¹ are in close agreement. However, the mode for Ω=0.32 rad s⁻¹ 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-$\left|\mathrm{\Psi }\right|$ modes occurring at locations of inflectional spanwise velocity and tending to have a single peak like those at $r_{\mathrm{0}}/R=\mathrm{0.26}$. Figure 16d–f shows the results for Geometry 2. At $r_{\mathrm{0}}/R=\mathrm{0.40}$, the increase in Ω reduces $\left|\mathrm{\Psi }\right|$ and makes double-peak modes such as those for Ω=0.45 and 0.9 rad s⁻¹ become a 2D, single-peak mode like the one for Ω=1.35 rad s⁻¹. At $r_{\mathrm{0}}/R=\mathrm{0.58}$, all modes collapse and present double peaks. At the outer radial location, the mode for Ω=0.45 rad s⁻¹ is nearly 2D, and the rise in Ω increases its obliqueness (i.e., increases $\left|\mathrm{\Psi }\right|$). 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 e^N 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.