Unsteady Blade Element-Momentum Method Including Returning Wake Effects

ABSTRACT: The wind energy research has grown substantially in the past few years, considerably fostered by the pursuit for a clean and sustainable energy source. Improvements on the design methods are increasingly needed. The purpose of this research is to investigate the use of the Loewy ́s lift deficiency function (LDF), also named Returning Wake Model, coupled with a non-stationary Blade Element-Momentum Method (BEM). The LDF simulates the influence of the wake behind the wind turbine on its capacity to generate power. It is expected that this model reduce the dependency of the several empirical parameters necessary in other wake models which are currently used. Aiming to validate the results obtained in this new approach they are compared with those provided by commercial computational software and they have proven to be very consistent. It is concluded that the method is feasible to be used as an efficient design and optimization tool of upwind horizontal axis wind turbine blades.


INTRODUCTION
Wind turbines operate in a hostile environment in which strong flow fluctuations, mainly due to the nature of the wind, can produce high and variable loads on its components.These loads, combined with the elastic behavior of the turbine structural components, compromise the overall energy generation efficiency and cannot be neglected during the design phase.The need for experimental and computational approaches to investigate the behavior of unsteady loads produced on a wind turbine blade has grown in proportion to the growing nominal power and size of the actual horizontal axis wind turbines.
The objective of this work is to evaluate the mathematical model of the Loewy´s lift deficiency function (LDF) to be included in a computational package enabling the optimal design of horizontal axis wind turbine blades.This model allows the inclusion of non-stationarities as well as an approximation of the effects of the downwind wake dynamics on turbine aerodynamic performance.This inclusion of non-stationarities and the effects of the downwind wake allow more realistic results for the aerodynamic loads than the simple Blade Element-Momentum (BEM) method, and consequently leads to better designs of structural components.
The overwhelming majority of computer packages and procedures actually use semi-empirical models to represent the influence of the wake in the overall aerodynamic performance of a wind turbine blade.The LDF (Loewy, 1957) is an analytical solution based on the classical Theodorsen theory (Theodorsen, 1935).The use of an analytical model reduces the dependence on experimental tests for the adjustment of empirical parameters needed to calibrate Silva, C.T. and Donadon, M.V.
these semi-empirical engineering models.Furthermore, an analytical resolution results in a more elegant mathematical solution, computationally feasible and sufficiently fast to provide blade design optimization.
The induced effects generated by the shed wake have been modeled using two general approaches: dynamic inflow methods and vortex wake methods.
The principles of the dynamic inflow approach are attributed to Carpenter and Freidovich (1953).The idea is to consider the unsteady aerodynamic lag of the inflow development over the rotor disk in response to changes in blade pitch inputs of changes in rotor thrust.They are written in the form of ordinary differential equations, with a time constant (or constants) representing the dynamic lag in the build-up of the inflow.One of its less satisfying aspects is that time constants must be obtained by experimental calibrations.
Vortex wake models are based on the assumption of an incompressible potential flow, with all vorticity being assumed concentrated within vortex filaments (which in the case of rotors require a coupling to the blade lift distribution).The induced velocity field can be determined through the application of the Biot-Savart low.Different approaches are encompassed ranging from prescribed to free vortex techniques.The prescribed wake models are strictly applicable when the operating conditions are nominally steady-state, i.e., in a steady wind.The free vortex methods have fewer potential limitations.They have been widely developed to be used in helicopters rotor analyses.Free vortex methods are based on discretized, finite-difference representation of the governing equations for the wake, and when solved, they track the evolution of discrete vortex elements through the flow.The number of discrete elements per vortex filaments can be very large, making the tracking process memory intensive and computationally demanding.Leishman (2002) presents more details about dynamic inflow models and vortex wake models.

THE CLASSICAL BLADE ELEMENT MOMENTUM METHOD
The BEM method presented by Glauert (1935) enables to calculate the steady loads and also the thrust and power using different settings on wind speed, rotational speed and pitch angle.The method couples the momentum theory with local events taking place at the actual blades.The blade is analyzed as a number of independent stream tubes.In each one, the induced velocity is calculated by performing the conservation of momentum, and the aerodynamic forces are found with the 2D aerodynamic theory and airfoil data.The stream tubes are discretized into N annular elements.The lateral boundary of the elements does not admit any flow across them.Some assumptions are made for the annular elements: no radial dependence, that is, one element cannot be affected by the others; the forces from the blade on the flow are constant in each annular element, corresponding to a rotor with a number of blades.
A correction known as Prandtl´s tip loss factor (Glauert, 1935) is introduced to correct this latter assumption in order to compute a rotor with a finite number or blades.
A relative velocity V rel seen by a blade section is a combination of axial velocity V 0 (1-a), in which a is the axial induction factor, and the tangential velocity (1-a') ϖr, where a' is the radial induction factor at the rotor plane (Fig. 1).The angle θ is the local pitch angle of the blade element, i.e., the local angle between the chord and the plane of rotation.It is a combination of the pitch angle, measured between the tip chord, the rotor plane and the twist of the blade, relative to the tip chord.ϕ is the flow angle, measured between the plane of rotation and the relative velocity.The local angle of attack α is obviously found.

Rotor plane
Unsteady Blade Element-Momentum Method Including Returning Wake Effects The algorithm for the BEM model can be summarized as the sequence of steps that follows.Since different control volumes are assumed to be independent, each blade element can be treated separately and the solution at one radius can be computed before solving another radius.The following algorithm is applied for each control volume.
Step (3): Compute the local angle of attack α.
Step (4): Read off the lift C l (α) and drag C d (α) coefficients from a table.
Step (5): Compute C n and C t , respectively normal and tangential aerodynamic force coefficients.
Step (7): If a and a' have changed more than a tolerable amount, repeat step (2), or else finish.
Step (8): Compute local loads on the blade element.
All equations necessary to perform the above algorithm can are described by Burton (2001) and Hansen (2007).

UNSTEADY BEM MODEL
In order to obtain good estimates of the annual energy production of a wind turbine, a steady BEM method is adequate to compute the steady power curve.But in reality the rotor of a wind turbine feels the inherent unsteadiness of the wind caused by atmospheric turbulence, wind shear and the presence of the tower.It is necessary to use an unsteady BEM method to compute realistically this variable behavior of the wind.
One simple model and additional coordinate systems can by placed at the wind turbine and its blades so it is possible to know the relative position of any blade element at any time.This simple model is depicted at Fig. 2.
An inertial system of coordinates is placed at tower base and named System 1. System 2 is non-rotating and fixed in the nacelle.System 3 is solidary to the shaft turbine and rotates with it, and system 4 is aligned with one of the blades.The tilt and cone angles are shown and the azimuthal position of the blade is set by de wing angle θ wing , not depicted.
The undisturbed wind velocity seen by the blade is found by a simple coordinate transformation clearly detailed by Hansen (2007).The essence of the BEM method is to determine the induced velocity and thus the local angle of attack.This is achieved by a summation of vectors, V rel = V 0 + V rot +W, all written at the element blade coordinate system (System 4, Fig. 3), in which the induction velocity is the term W.
With the induced velocity known, the flow angle and angle of attack are found.Bramwell et al. (1976) states that Glauert's relation between thrust and this induced velocity for a gyrocopter in Silva, C.T. and Donadon, M.V.
forward flight (similar to the result of the lifting line for an elliptically loaded circular wing) is: in which n is the unit vector in the direction of the thrust, which in System 3 has the coordinates n=[0,0,1] T .It is assumed that only the lift contributes to the induced velocity, and that the induced velocity acts in the opposite direction to the lift.The force from this blade at radial position is assumed to affect the air in the area dA = 2πrdr/B, so that all B blades cover the entire annulus of the rotor disc at radius r.
The following expression can be derived for one blade, according to Hansen (2007), For the tangential component a similar expression is postulated, in which F is Prandtl's tip loss factor.If the rotor is yawed (and/or tilted), there will be an azimuthal variation of the induced velocity, so that it is lower when the blade is pointing upstream in relation to when the same blade, half a revolution later, is pointing downstream.The physical explanation for it is that a blade pointing downstream is deeper into the wake than a blade pointing upstream.This means that an upstream blade sees a higher wind speed and thus produces higher loads than the downstream blade, which produces a beneficial yawing moment that will try to turn the rotor more into the wind, thus enhancing yaw stability.The yaw model describes the distribution of the induced velocity.If a yaw model is not included, the BEM method will not be able to predict the restoring yaw moment, according to Hansen (2007): in which the wake skew angle, X, is defined as the angle between the wind velocity in the wake and the rotational axis of the rotor.θ 0 is the angle in which the blade is deepest into the wake.The skew angle can be found as: The skew angle is assumed to be constant with the radius and can be computed at a radial position close to r/R = 0.7.
The induced velocity is now known at the new azimuthal position at time t+∆t, θ wing (t+ ∆t)= θ wing (t)+ ω∆t.The angle of attack can thus be evaluated from the equation and the lift and drag coefficients can be looked up from a table.The normal, p z , and tangential, p y , loads can be determined from: The algorithm can be resumed like the following: Step (1): Initialize all necessary data (geometry and run parameters); Step (2): Initialize the position and velocity of blades; Step (3): Discretize the blades into N elements; Step (4): Initialize the induced velocity; -for n=1 to max time step (t=n∆t) -for each blade -for each element 1 to N Rotor plane Unsteady Blade Element-Momentum Method Including Returning Wake Effects Step (5): Compute relative velocity to the blade element using old values for induced velocity; Step (6): Calculate flow angle and thus the angle of attack (Eq.1); Step (7): Determine static drag and lift coefficients from tables; Step (8): Compute lift and drag for each blade element (Eq.8); Step (9): Compute loads (Eq.7); Step (10): Compute new equilibrium values for induced velocities (Eqs.3 and 4); Step (11) Calculate the azimuthal variation from Eq. 5 and compute the induced velocity for each blade; Step (12): Compute momentum, thrust and power; Step (13): Increment time step and repeat from step (5).
The equations of the BEM method must be solved iteratively.The flow angle and thus the angle of attack depend on induced velocity.But the described algorithm is unsteady, therefore time is used as relaxation.After blades have moved in one time step an azimuthal angle of ∆θ wing =ω∆t (for small ∆t´s), values from the previous time step are used on the right hand side of equations for W when updating new values for induced velocity.This can be taken into account since induced velocity changes relatively slowly in time.This eliminates the need of calculating the induction factors and use of tolerance.

DETERMINISTIC WIND MODEL
The exponential model used to simulate wind shear.It controls the wind speed according to the altitude, and the shear parameter used is 0.2.Detailed information about this wind shear model was described by Hansen (2007).
The wind is also influenced by the presence of the tower.The simple model used to simulate the tower shadow assumes potential flow.All details about this simple model were also described by Hansen (2007).This model is a bad approximation for a downwind machine, in which each blade passes the tower wake once every revolution.However, for an upwind machine, the object of this study, the model provides good estimation.Also, the turbulent part of the real atmospheric wind should be added for a realistic time simulation for a wind turbine.For this initial investigation no atmospheric turbulence is added to the simulation.

LOEWY´S LIFT DEFICIENCY FUNCTION
The problem of calculating the aerodynamic loading on an oscillating profile was first approached by Glauert (1929), but it was only properly solved by Theodorsen (1935).Theodorsen's approach gives the solution for unsteady aerodynamic loading on a 2D oscillating airfoil in an inviscid and incompressible flow, and subject to the assumption of small disturbances.Theodorsen's problem is to obtain the solution for loading on the surface of the airfoil under the condition of forced harmonic oscillations.
For a simple harmonic motion of the airfoil the solution given by Theodorsen in a way that represents a transfer function relating the forcing input (angle of attack) and the aerodynamic response (pressure distribution, lift, and pitching moment).The approach is summarized by Bisplinghoff et al. (1955).See also Bramwell et al. (1976) and Johnson (1980) for a detailed exposition of the theory.
Theodorsen's theory is not suitable for studies involving rotors.In these types of problems sections of blades can find wake vorticity due to other rotor blades, as well as the returning wake from the blade in question.This fact was recognized by Loewy (1957) and Jones (1958).They built a two-dimensional model of a 2-D blade section with a returning shed wake, as shown in Fig. 4.
As in Theodorsen's model, the shed wake is modeled with 2-D flat surfaces of vortices, but now with a series of surfaces below the airfoil section with a vertical separation h, which depends on the speed induced by the rotor disc V and the number of blades N b .Loewy shows that, in this case, the lift on the blade can be expressed by replacing the function Theodorsen by the named Loewy's function.
( ) in which it is known for Loewy's function of Loewy's lift deficiency function.
For a rotor with N b blades the complex function W is written as: The wake spacing ratio h/b can be determined from the spacing of vortex sheets that are laid down below the rotor.If an average induced velocity v i =λΩR is assumed, then during a single rotor revolution the shed wake generated by a single blade will be at a distance h=(2π/Ω)v i below the rotor.For a multiple blade rotor the spacing is (2π)v i / ΩN b , i.e.
in which σ is the rotor solidity.
Representative results from Loewy's theory show that the main consequence of including shed vorticity below the blade is that it serves to amplify or attenuate the unsteady lift response, depending on the reduced frequency, wake spacing and wake phase.The most important effects are for lower reduced frequencies, with oscillations at the harmonics of the rotor rotational frequency (Leishman, 2000).
Several aspects related to the non-stationarity of the operating environment of a wind turbine need to be addressed during its project.Among them there are the main variations in wind speed (gust and wind shear), dynamic inflow, yaw and tower shadow, turbulence, wake dynamics and interactions blade/wake, and the dynamic stall.The adoption of a non-stationary BEM method using a dynamic inflow model is a solution that tends to enhance the results, proving to be a good option to introduce the study of non-stationarity in the design of a wind turbine blade.However the same problems related to the physics of the method remain.
Unsatisfactory aspects of the inflow theory are the so-called dynamic time constants employed in the methods.They are developed using the concept of apparent mass and inertia of the fluid surrounding the rotor (noncirculatory effect) as opposed to the delay of the dynamic evolution of wake vortices (circulatory effects).The concept of apparent mass applied to the rotor also assumes equivalence between the apparent force of the rotor disk accelerating in a stopped fluid and the force in a fluid accelerating through a permeable actuator disc, which certainly is not an accurate analogy.
Loewy proposes a solution to the problem of aerodynamics of rotors affected by non-stationarity generated by shed wake.It is based on Theodorsen's solution applying a suitable physical model for  (Leishman, 2000).
Unsteady Blade Element-Momentum Method Including Returning Wake Effects rotorcraft aerodynamics.The results obtained in the solution of problems related to helicopter hovers, and duly validated with experimental results, confirm the efficiency of the method.
The aerodynamics of a helicopter hover resembles in many aspects the aerodynamics of the blades of a wind generator.This similarity, coupled with the need for a mathematical model for the wind generator blade design that considers the conditions of non-stationarity of the phenomena involved in its aerodynamics, are the main motivators for this work.
The corrections using Loewy's model are applied for the Equation 2.8 in order to provide the lift deficiency described above.The corrected equations can be rewritten.No corrections are applied to the drag force.

VERIFICATION AGAINST AN INDUSTRY-STANDARD SOFTWARE
The commercial package used for result validation is the GH Bladed V4.1, distributed by GL Garrad Hassan (2012), an independent renewable energy consultancy.GH Bladed is an industry-standard integrated software package for the design and certification of onshore and offshore turbines.It provides user with a design tool that has been extensively validated against measured data from a wide range of turbines and enables the conduction of the full range of performance and loading calculations (Garrad Hassan & Partners, 2011).
Its manual postulates that GH Bladed uses the same methods employed at the present work.That is, combined blade element and momentum theory, wake rotation with radial induction, tip and hub loss models, which is suppressed for the present comparison, dynamic wake model and dynamic stall, also suppressed.
An educational version of de GH Bladed GH Bladed 4.1, with a limitation of 10 blade elements, is used.
The object of study is a 2MW wind turbine.The main used parameters and data are shown in Tables 1 and 2.
All necessary data, like angle of attack, lift, drag and moment coefficients are described by McGee and Beasley (1976).
Just for visualization, BEM NE has a graphical interface that shows a simplified illustration of the wind turbine and his main geometrical parameter.This is shown at Fig. 5.
The developed computational routine uses these same data to be compared with the results of both packages.It also uses de same model for wind shear, i.e., the exponential  vertical shear model with wind shear coefficient 0.2 and the potential flow model for the tower shadow.Details about the models for wind shear end tower shadow were described by Hansen (2007).
Figures 6 and 7 illustrate the wind condition for one blade element.Figure 6 is an illustration of the magnitude of the wind for a blade element.Figure 7 shows the wind speed distribution influenced by the presence of the tower (tower shadow).
Unsteady Blade Element-Momentum Method Including Returning Wake Effects

RESULTS AND DISCUSSION
The most important aerodynamic and performance results are shown and compared here.It is important to mention that the GH Bladed package has a pitch control module that changes the pitch angle of the blades in each time step according to the wind speed, rotor speed and output power.The developed computational routine, named here BEM NE, uses this same pitch angle control map in order to achieve the correct level of the shaft power, once it does not have a control module.This control map is depicted in Fig. 8b as also the shaft power (Fig. 8a) against different wind speeds.
The blades of the wind turbine used in this study are designed to produce 2MW power between wind speeds of 12 and 25 m/s, at the rotor speed 18 rpm.Results at 12 and 18 m/s wind speed are compared.
At first the results obtained for 12 m/s wind speed is shown.Figures 9 and 10 compare the geometrical results of flow angle and angle of attack, respectively.The results shown are obtained from the blade element located 26.407 m from the rotor center.This blade element is the most representative one.It is located in 70% length position of the blade.
During the 12 seconds of simulation, deviations on the results are minimal, as demonstrated in Figs. 9 and 10.Mean relative error at inflow angle is -1.47% and mean relative error at angle of attack is -2.71%.
Figures 11 and 12 show the comparison of the relative wind speed and the corresponding lift coefficient of the investigated blade element.For the BEM NE results for the lift coefficient, Loewy's lift deficiency function is applied to correct its value.
Again, the mean relative errors are small: respectively +0.28% e +7.68% for relative wind speed and lift coefficient.
The most important result, the output shaft power, is shown in Fig. 13.The relative mean error is now -0.66%.
Results obtained until this point show a good correlation between the GH Bladed and the BEM NE for the wind speed 12 m/s.Now the same comparison of results for the wind speed 18 m/s is shown in the middle of the operation wind speed range.Figures 14 until 18 show these comparison.For wind speed 18 m/s the correspondent pitch angle is 14.9 degrees, as shown in Fig. 8.
Relative mean errors are compatible with the ones obtained for wind speed 12 m/s.The most important parameter, the shaft power, has a relative error of -4.1%.
It is observed that the power curve obtained from BEM NE has more spikes then the one from GH Bladed.The discontinuities coincide with the passage of the blade through the tower shadow.Even though both packages use the same model for the tower shadow, these discrepancies can be justified by the absence of any kind of pitch control on the BEM NE.At the current stage of the present study, the model is simplified because its primary objective is to investigate the viability of LDF as a simulation for the influence of the wake shed behind the rotor.Therefore the pitch angle is maintained constant during all time steps simulation.On the contrary, Bladed has a variation of the pitch angle during the simulation, as shown in Fig. 19.

Figure 3 .
Figure 3. Relative and induced wind speeds in System 4.

Figure 8 .
Figure 8. Shaft Power (a) and Pitch angle (b) control map.

Table 1 .
General characteristics of rotor and turbine.