Detailed analysis of the blade root flow of a horizontal axis wind turbine

The root flow of wind turbine blades is subjected to complex physical mechanisms that influence significantly the rotor aerodynamic performance. Spanwise flows, the Himmelskamp effect, and the formation of the root vortex are examples of interrelated aerodynamic phenomena that take place in the blade root region. In this study we address those phenomena by means of particle image velocimetry (PIV) measurements and Reynolds-averaged Navier–Stokes (RANS) simulations. The numerical results obtained in this study are in very good agreement with the experiments and unveil the details of the intricate root flow. The Himmelskamp effect is shown to delay the stall onset and to enhance the lift force coefficient Cl even at moderate angles of attack. This improvement in the aerodynamic performance occurs in spite of the negative influence of the mentioned effect on the suction peak of the involved blade sections. The results also show that the vortex emanating from the spanwise position of maximum chord length rotates in the opposite direction to the root vortex, which affects the wake evolution. Furthermore, the aerodynamic losses in the root region are demonstrated to take place much more gradually than at the tip.


Introduction
The aerodynamic design of wind turbine blades is subjected to important levels of uncertainty.As a matter of fact, not only transient operational states can pose a challenge to the wind turbine designer but also seemingly simple cases involving steady operation under axisymmetric, uniform inflow conditions (Leishman, 2002;Schepers, 2012).This is especially true for the tip and root regions of the blades, where the flow is three-dimensional and strongly influenced by the trailing vortices (Micallef, 2012).

Spanwise flows and Himmelskamp effect
At the root of the blade, the angle of attack (AoA) is usually considerably higher than at the tip.This increases the complexity of the flow since it often leads to flow separa-tion, which in this part of the blade generally gives rise to the Himmelskamp effect (Himmelskamp, 1947).This effect delays the stall onset and enhances the lift force as compared to non-rotating blades operating at the same AoA.The Himmelskamp effect, also known as stall delay or rotational augmentation, has been studied by many authors both experimentally (Schreck and Robinson, 2002;Sicot et al., 2008;Ronsten, 1992) and numerically (Guntur and Sørensen, 2014;Herráez et al., 2014;Schreck et al., 2007), although it still remains far from being well understood and characterized.It mainly affects the blade root region and is known to be closely related to the existence of spanwise flows in the boundary layer.Snel et al. (1993) were the first to propose a correction model to be applied to 2-D airfoil characteristics in order to account for this effect in blade element momentum (BEM) and other engineering tools that rely on 2-D airfoil data.More correction models have been developed since then (e.g.Chaviaropoulos and Hansen, 2000;Bak et al., 2006;Raj, 2000;Corrigan and Schillings, 1994).However, Breton et al. (2008) and Guntur et al. (2011) proved that their accuracy is still a critical issue.Currently, a major impediment in the development of accurate correction models is the incomplete understanding of the physical mechanisms.It is worth noting that, up to now, the study of the Himmelskamp effect has been mostly focused on post-stall conditions.Consequently, very little is known about its onset at moderate angles of attack.

The root vortex
One fundamental feature of the root (and tip) flow is the formation of trailing vorticity that rolls up into a discrete vortex.Several authors have attempted to capture experimentally the root vortex in the near wake of a wind turbine.However, as Vermeer et al. (2003) highlighted, this can be extremely difficult to achieve due to the fact that the near wake usually does not present a distinctive, well-defined root vortex (in opposition to the tip vortex).Many wind tunnel experiments with model wind turbines confirmed this.For instance, Massouh and Dobrev (2007) and Haans et al. (2008) also came to that conclusion after studying a wind turbine rotor wake with particle image velocimetry (PIV) and hot film wake measurements, respectively.Furthermore, Ebert and Wood (2001) and Sherry et al. (2013) observed by means of PIV (among other measurement techniques) that the root vortex diffuses very rapidly.The PIV measurements performed by Akay et al. (2012) on two different rotors demonstrated that the evolution and strength of the root vortex highly depends on the blade root geometry and the spanwise distribution of circulation.In a subsequent work also based on PIV measurements, Akay et al. (2014) suggested that the flow in the root region is driven by the bound vorticity.
The study of the root (and tip) vortices can also be addressed by means of numerical simulations.For this purpose, large eddy simulations (LESs) are commonly combined with actuator line models (Ivanell et al., 2007;Troldborg et al., 2007;Nilsson et al., 2015).This technique is very useful for analysing the evolution of the trailing vortices in the wake.
However, it implies a very strong simplification of the blade geometry, which makes it unsuitable for studying the origin of the root and tip vortices.This is well exemplified in van Kuik et al. (2014), where it is concluded that the fact that actuator line models disregard the chordwise bound circulation at the blade tip prevents them from correctly computing the tip vortex trajectory in the vicinity of the blade.The same article also shows that full blade Reynolds-averaged Navier-Stokes (RANS) simulations as well as panel code computations allow a much more realistic study of the tip vortex formation mechanism.Indeed, the use of a panel code allowed Micallef et al. (2012) to study the origin of the tip vortex on a wind tunnel model rotor, unveiling the complex distribution of bound vorticity at the blade tip.However, to the best of our knowledge, the formation of the root vortex has not been addressed until now.

Scope and outline
This article aims at gaining insight both experimentally and numerically into the root flow of a horizontal axis wind turbine operating under design conditions.The focus is put on two important and interrelated aspects of the root flow that, as above explained, are insufficiently understood so far: 1. spanwise flows and onset of the Himmelskamp effect at moderate angles of attack (design operating conditions); 2. the formation of the root vortex.

Experimental setup
The scope of the experimental campaign is to measure the three components of the flow over the root region of a wind turbine blade.This is achieved by means of stereoscopic PIV.
The measurements are carried out in the Open Jet Facility of the Faculty of Aerospace Engineering at the Delft University of Technology.This wind tunnel has an octagonal open jet with an equivalent diameter of 3 m.The studied wind turbine consists of a two-bladed rotor with a diameter of 2 m.The chord and twist distributions are shown in Fig. 1.Table 1 shows the airfoil type distribution along the span.
The measurement campaign includes both a spanwise and a chordwise configuration of the PIV windows.The spanwise measurements are carried out at different azimuth angles for capturing the evolution of the near wake.In the chordwise configuration the PIV windows are orthogonal to the blade axis around the blade chord.This configuration, which included 40 different radial positions, offers the best insight into the flow around the blade surface and is the one presented in this work.A diagram of the corresponding experimental setup with the global coordinate system is displayed in Fig. 2.
The measurements are phase-locked and phase-averaged with the azimuthal position of the rotor blade rotation.This allows for reconstruction of the flow over each blade section after measuring the pressure and suction sides separately.
The rotor operated at rated conditions with a freestream wind speed U ∞ = 6 m s −1 and a rotational speed  Further details about the experimental setup can be found in Akay et al. (2014).

Numerical method and computational mesh
The simulations presented in this work are based on the RANS method, and they are performed with the open-source code OpenFOAM (2015).The computational model solves the incompressible Navier-Stokes equations using a finitevolume approach for the spatial discretization.The convective terms are discretized with a second-order linear-upwind scheme.For the viscous terms a second-order, centraldifferences linear scheme is employed.The use of a noninertial reference frame and the addition of the Coriolis and centrifugal forces to the momentum equations allows to account for the rotation of the system.The SIMPLE algorithm is employed for enforcing the pressure-velocity coupling.The turbulence in the boundary layer is modelled by means of the k − ω shear-stress transport (SST) model proposed by Menter (1993).This model has been proved to be suitable for the simulation of wind turbine blades (Bechmann et al., 2011;Johansen and Sørensen, 2004;Le Pape and Lecanu, 2004;Sørensen et al., 2002).However, the implicit assumption of fully turbulent flow might be a source of uncertainty since the existence of laminar-to-turbulent transition can not be completely ruled out.
The grid is generated with the software Pointwise (2015).The hub and nacelle geometries are disregarded in order to keep the mesh as simple as possible.This approach, which is based on the assumption that the hub and nacelle do not influence substantially the blade root flow, is usually followed when structured meshes are used for simulating wind turbine blades (Johansen et al., 2002;Sørensen et al., 2002;Le Pape and Lecanu, 2004;Schreck et al., 2007;Bechmann et al., 2011).The mesh exploits the symmetry of the rotor by modelling only one half of it and using periodic boundary conditions.The computational domain is shown in Fig. 3 and it consists of two independent block-structured grids connected by means of a so-called arbitrary mesh interface.The outer grid is a semi-sphere with a radius of 22R, where R is the blade radius.The inner grid, which contains the blade, is a cylinder with a radius of 1.1R and height of 1.1R.The motivation for using two structured grids connected by an interface is to independently control the mesh resolution in the proximity of the blade and in the far field.The total number of cells is 9.8 × 10 6 .The blade surface mesh (see Fig. 4) contains 130 cells along the chord, while 210 cells are used in the spanwise direction.In order to properly resolve the boundary layer, the height of the first cell in the normal direction to the blade surface is set to 5 × 10 −6 m, which ensures that Y + is smaller than one along the whole blade.The semi-spherical outer boundary employs a boundary condition that changes its behaviour depending on the direction of the flow: in regions where the flow goes in, it works like a Dirichlet boundary condition assuming a predefined value of the velocity field; in regions where the flow goes out, it enforces a zero-gradient condition (Neumann condition).The   and 7, and they will be just neglected in the interpretation of the results.
Figure 5 shows the azimuthal velocity component in an inertial reference frame for the radial stations r = 0.26R, r = 0.35R (the position of maximum chord length) and r = 0.45R.The results are normalized with the free-stream wind speed U ∞ .The agreement between experimental and numerical results is fairly good for all the studied radial positions.At r = 0.35R and r = 0.45R, the azimuthal velocities are positive over the whole suction side except in the trailing edge region.However, at r = 0.26R the suction side presents negative velocities from the mid-chord until the trailing edge.This does not necessarily imply flow separation and recirculation, though, since the relative velocity might still remain positive along the whole suction side.Indeed, at that radial position the local circumferential velocity caused by the blade rotation is 1.82/U ∞ , which implies that the flow remains attached up to U θ = −1.82/U∞ .In Sect.3.3 the lack of separation is demonstrated by means of the wall shear stress.The axial flow component is shown in Fig. 6.The axial velocity over the second half of the suction side becomes smaller with increasing radial position.This is in fact just a geometrical effect that occurs as a consequence of the twist of the rotor blade.At r = 0.45R the orientation of the second half of the suction side surface is slightly upstream, so the axial velocity becomes negative if the flow is attached.The effect is smaller for lower radial positions because of the larger twist angle, which neutralizes the mentioned geometrical effect.The numerical results are consistent with the experiments, although at r = 0.45R the axial velocity over the suction side is slightly overpredicted.
Figure 7 displays the radial velocity component for the three considered spanwise positions.At r = 0.26R, this velocity is very strong on the suction side.However, at r = 0.35R and r = 0.45R it becomes much weaker.Hence, the presence of spanwise flows seems to be limited to the innermost region of the blade.The agreement between experiments and simulations is again very good for all three stations.The velocity field 10 mm off the suction side has been extracted from both the numerical results and the available PIV data (including 40 different radial positions between r = 0.17R and r = 0.65R) in order to study the flow in the proximity of the blade surface in more detail.Figure 8 shows that the azimuthal component is always positive outboard of the position of maximum chord length (r = 0.35R).Below that position, a significant region of the blade presents negative azimuthal velocities close to the trailing edge.At r ≈ 0.3R, this effect is stronger in the PIV measurements than in the numerical results.Apart from that, the numerical results are in very good agreement with the experimental results.
The axial flow velocity component is displayed in Fig. 9. Outboard of the radial position of maximum chord length, the axial velocity becomes negative from the mid-chord towards the trailing edge.The same effect has already been discussed in relation to Fig. 6.The fact that the numerical results somewhat underpredict this effect, which is caused by the relative position of the suction side to the rotor plane, might indicate a possible small deviation in the pitch angle or some uncer-  tainty in the PIV fields.For a detailed discussion of the experimental uncertainties, the interested reader is referred to Akay et al. (2014).The agreement between PIV and computational fluid dynamics (CFD) in the root region is very good, which is a clear indication that the wake and blade inductions are correctly predicted with the current CFD model.
Figure 10 presents the distribution of radial velocity along the blade suction side.The experimental results show a substantial spanwise flow in the leading edge region from r = 0.45R outwards.This is rather surprising, since the flow in that region is fully attached (as shown, for instance, in Fig. 8) and it is far away from the tip and root, where spanwise flows are usually expected.The numerical results show much smaller radial velocities in that region.At present, the authors do not have a solid explanation for this discrepancy since both experimental and numerical uncertainties might play a role in the mentioned discrepancy.In the root region both PIV and CFD show evidence of strong spanwise flows in the proximity of the trailing edge, although the simulation tends to underpredict the spanwise flow in the region 0.3R < r < 0.35R, as also happened with the azimuthal velocity (Fig. 8).As stated earlier, this might be caused by a slight deviation in the pitch angle.Other than this, the consistency between PIV and CFD is again very good, which gives confidence in the reliability of the numerical model in predicting the complex flows of the root region.

The source of the spanwise flows
Two different explanations have been proposed in the literature for explaining the origin of the spanwise flows: 1. Spanwise pressure gradients: the dynamic pressure over the blade surface is inversely proportional to the radial position.Hence, the air is assumed to travel from the root towards outer positions as a consequence of spanwise pressure gradients (Schreck and Robinson, 2002;Schreck et al., 2010).
2. Centrifugal force: the centrifugal force that acts on the bottom of the boundary layer (i.e. the region where the flow is not detached from the surface) pushes the flow towards the tip (Du and Selig, 1999;Lindenburg, 2003;Guntur and Sørensen, 2014).
The numerical results help to elucidate which is the actual source of the spanwise flows.In Fig. 11 the computed isobars of the blade suction side are compared with the limiting streamlines obtained from the wall shear stress.As can be seen, the surface pressure does not present significant spanwise gradients.It is worth remarking that the same observation was made in the analysis of the MEXICO wind tunnel experiment (Herráez et al., 2014).Therefore, we conclude that the centrifugal force is the main source of spanwise flows.
Figure 11 also shows how the Coriolis force progressively redirects the spanwise flow coming from the root towards the trailing edge, which causes the flow to follow a curved trajectory.

Onset of the Himmelskamp effect
Figure 12 compares the pressure coefficient C p distribution obtained from the blade at r = 0.26R with the C p distribution extracted from a 2-D simulation at the same Reynolds number (Re ≈ 1×10 5 ) and same angle of attack (AoA ≈ 13 • , computed using the method proposed by Shen et al., 2009).The 2-D simulation is a RANS computation performed with the k−ω SST turbulence model.Also, the 2-D mesh is equivalent to the 3-D mesh except for the third dimension.Experimental results of the same 2-D airfoil with Re = 1 × 10 6 are displayed as well.The 2-D experimental and numerical results present some disparity in the region of the suction peak, but apart from that they are very similar in spite of the difference of Reynolds number.However, the 3-D results exhibit some important differences.The slope of the adverse pressure gradient is substantially reduced, which leads to a delay of the separation point.The separation point can be approximately identified as the point where the adverse pressure gradient meets the region with zero pressure gradient (i.e. the region where the flow is separated).In the 2-D airfoils the separation point is located at x/c ≈ 0.4.In the 3-D case, the adverse pressure gradient presents a kink at x/c ≈ 0.5, but it stays negative for the whole chord length, which seems to indicate that the flow remains attached.However, Sicot et al. (2008) concluded that rotating blades can present separation even in regions of adverse pressure gradient.In order to verify if there is separation in the 3-D case, the skin friction coefficient C fx on the suction side is displayed in Fig. 13 for both the 2-D and 3-D simulations.In the 2-D case, C fx becomes positive at x/c = 0.39, indicating that the flow separates exactly at that point (in good agreement with the estimation from the C p distribution).In the 3-D case, C fx becomes zero at x/c = 0.52, but it recovers directly afterwards, remaining always negative.This confirms that the flow stays completely attached all along the chord.The point where C fx becomes zero is actually the place where the chordwise flow is deflected towards the spanwise direction.The same happens for all other radial positions at the root.Therefore, the transition between the chordwise and spanwise flows in edge.This resembles the behaviour of the 2-D case in the region with zero pressure gradient.Finally, it is worth highlighting that the 3-D case presents a smaller suction peak than the 2-D case.
The resulting lift and drag coefficients (C l and C d , respectively) for the 2-D and 3-D cases are presented in Table 1.C l is increased by approximately 9 % as a consequence of the Himmelskamp effect, whereas C d does not seem to be influenced at all.This is also in agreement with our observations from the MEXICO turbine, where the Himmelskamp effect had a very limited influence on the drag (Herráez et al., 2014).
Wind Energ.Sci., 1, 89-100, 2016 Estimating the validity of the above described results for larger wind turbines is, however, not so straightforward.On the one hand, the analysed turbine operates at the tip speed ratio λ = ωR/U ∞ = 7, which is also realistic for full-scale wind turbines working at nominal conditions.The NREL 5 MW wind turbine, for instance, also presents the same rated tip speed ratio (Jonkman et al., 2009).This is important because it contributes to maintain the same balance between the centrifugal and Coriolis forces, which is fundamental for the Himmelskamp effect.Lindenburg (2003), for instance, estimated that the change of aerodynamic lift and drag due to the Himmelskamp effect is proportional to the square of the local speed ratio ωr/U rel .Dowler and Schmitz (2015) also included a similar parameter, namely 2U rel /ωr (obtained directly from the ratio between the Coriolis and centrifugal forces) in their stall delay model.Interestingly, they also estimated the change of the lift force to be proportional to the square of the mentioned parameter.
On the other hand, the local blade solidity c/r, which has also been identified as a fundamental parameter for the Himmelskamp effect (see, e.g., Snel et al., 1993;Chaviaropoulos and Hansen, 2000;Lindenburg, 2003;Dowler and Schmitz, 2015), differs substantially between the TU-Delft and the NREL 5 MW turbines: at r = 0.26R (i.e. the radial position studied in Figs. 12 and 13), c/R ≈ 0.15 for the TU-Delft turbine and c/R ≈ 0.07 for the NREL 5 MW one.For some authors (e.g.Chaviaropoulos and Hansen, 2000;Lindenburg, 2003), the change in lift force is linearly proportional to the c/r parameter, whereas for other authors (e.g.Snel et al., 1993;Dowler and Schmitz, 2015) it is proportional to the square of the mentioned parameter.In any case, it can be inferred that the large discrepancy in the local blade solidity between both turbines would lead to a weaker Himmelskamp effect in the NREL 5 MW turbine.This conclusion is also valid for other wind turbines of the same size since they usually present a similar local blade solidity.

The origin of the root vortex
The bound vorticity γ can be computed as the difference in the velocity outside the boundary layer of the pressure and suction sides.γ can then be decomposed into a radial γ radial and a chordwise γ chordwise component.Figure 14 shows both components side by side.
γ radial is concentrated around the 1/4 chord position, as might be expected.The radial circulation radial can be computed from γ radial as its integral along the chord: where le and te are the leading and trailing edges, respectively.The use of the Kutta-Jukowski theorem then allows for the sectional lift to be computed: where ρ is the air density and U rel is the relative velocity in the plane perpendicular to the blade axis.Owing to the γ radial distribution from Fig. 14 and the strong link between γ radial and the lift, it can be concluded that the lift force is generated almost exclusively in the first half of the chord all along the blade.The decay of γ radial , and hence the decay of the lift force, is much more sudden at the tip than at the root.As a consequence, the root losses take place much more gradually than the tip losses.This should be taken into account by the correction models used, for example, in the BEM and actuator line methods.γ radial is transformed into γ chordwise at the tip and root before becoming trailed free vorticity, which gives rise to the tip and root vortices.This is evidenced in Fig. 14b, where it can be seen that the tip and root regions present substantial γ chordwise in the proximity of the trailing edge.γ chordwise is distributed over a larger spanwise range at the root than at the tip, which is in agreement with the gradual root losses earlier described.Van Kuik et al. (2014) obtained very similar results at the tip of a different rotor but the root was not studied.In the innermost region of the blade the sign of γ chordwise at the trailing edge is opposite to that of the tip (as one would expect from a horseshoe vortex model).However, in the region of maximum chord length, γ chordwise at the trailing edge presents the same sign as the tip vortex.The negative γ chordwise at the root implies an outward motion of the flow over the part of the suction side where the azimuthal velocity is slow (see Fig. 8).However, the positive γ chordwise in the region of maximum chord leads to an inward flow motion.Akay et al. (2014) studied the wake of the same wind turbine with PIV and indeed observed the presence of an outward flow for r < 0.25R and the existence of an inward flow in the radial range 0.25R < r < 0.35R.Furthermore, Medici and Alfredsson (2006) did similar observations in their experimental wake study of a different wind turbine.The present results not only confirm the mentioned experimental observations but also explain the origin of this aerodynamic behaviour.
Figure 15 shows the bound vorticity vectors over the blade.From this figure it is evident how the direction of the bound vorticity γ changes at the root.As can be seen, at the most inboard part of the blade (r < 0.35R), γ chordwise dominates the flow over the second half of the chord, indicating that vorticity is trailed along that region.The fact that γ chordwise is distributed over such a large area of the root might explain that the root vortex does not present a well-defined, distinctive structure, as Vermeer et al. (2003), Massouh and Dobrev (2007), and Haans et al. (2008) reported in their experimental wake studies of different wind turbines.Furthermore, the existence of two adjacent root regions with counter-rotating γ chordwise might also explain the fast diffusion of the root vortex reported by Ebert and Wood (2001) and Sherry et al. (2013).

Conclusions
The use of PIV measurements and RANS simulations enabled the analysis of the flow in the root region of a wind turbine blade operating under design conditions with axisymmetric inflow.The following conclusions are drawn: -The RANS method is capable of capturing accurately the main features of the root flow of wind turbine blades operating under design conditions.
-The spanwise flows in the root region are caused by the centrifugal force and not, as some authors have suggested, by radial pressure gradients.
-Even at relatively moderate angles of attack (AoA ≈ 13 • ), the interaction of the centrifugal and Coriolis forces can give rise to the Himmelskamp effect.
-The influence of the Himmelskamp effect on the sectional C p distribution is twofold: on the one hand the suction peak is reduced, while on the other hand the separation point is delayed (indeed, in our case, the separation is completely avoided).As a consequence of both counteracting effects, the influence of the Himmelskamp effect on the loads is weaker than on the C p distribution.
-The reduction of the aerodynamic performance is more gradual at the root than at the tip.Tip/root loss correction models (as used, for example, in BEM simulations) should account for this effect.
-The trailing vorticity in the spanwise position of maximum chord length presents the opposite sign to that at the blade root.This contributes to the diffusion of the root vortex.
We recommend considering these points for a better characterization of the root flow of wind turbine blades.This can help to reduce the uncertainty in the blade design process, which would in turn contribute to making the turbines more cost-effective.

Figure 1 .
Figure 1.Chord and twist distribution along the blade span.

Figure 2 .
Figure 2. Experimental setup with the chordwise PIV measurement configuration.The coordinate system used in this work is also displayed.The azimuthal direction θ is opposite to the direction of rotation ω.

Figure 3 .
Figure 3. Schematic representation of the computational domain.The inner cylinder represents the arbitrary mesh interface.

Figure 4 .
Figure 4. Detail of the surface mesh in the blade root region.

Figure 5 .
Figure 5. Experimental and numerical results of the azimuthal velocity component at different blade spanwise positions.

Figure 6 .
Figure 6.Experimental and numerical results of the axial velocity component at different blade spanwise positions.

Figure 7 .
Figure 7. Experimental and numerical results of the radial velocity component at different blade spanwise positions.

Figure 8 .
Figure 8. Experimental and numerical results of the azimuthal velocity component 10 mm off the blade suction side.

Figure 9 .
Figure 9. Experimental and numerical results of the axial velocity component 10 mm off the blade suction side.

Figure 10 .
Figure 10.Experimental and numerical results of radial velocity component 10 mm off the blade suction side.

Figure 11 .
Figure 11.Isobars and limiting streamlines over the suction side of the blade root region (obtained from the numerical results).

Figure 14 .
Figure 14.Radial and chordwise components of the bound vorticity.

Figure 15 .
Figure 15.Bound vorticity vectors over the blade.