Question? Leave a message!




How wind Turbines produce Electricity

how wind turbines work and how wind turbines are installed and how wind turbines affect the environment and how wind turbine blades work
3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 1 1 General Introduction to Wind Turbines Before addressing more technical aspects of wind turbine technology, an attempt is made to give a short general introduction to wind energy. This involves a very brief historical part explaining the development of wind power, as well as a part dealing with economy and wind turbine design. It is by no means the intention to give a full historical review of wind turbines, merely to mention some major milestones in their development and to give examples of the historical exploitation of wind power. Short Historical Review The force of the wind can be very strong, as can be seen after the passage of a hurricane or a typhoon. Historically, people have harnessed this force peacefully, its most important usage probably being the propulsion of ships using sails before the invention of the steam engine and the internal combust- ion engine. Wind has also been used in windmills to grind grain or to pump water for irrigation or, as in The Netherlands, to prevent the ocean from flooding low-lying land. At the beginning of the twentieth century electricity came into use and windmills gradually became wind turbines as the rotor was connected to an electric generator. The first electrical grids consisted of low-voltage DC cables with high losses. Electricity therefore had to be generated close to the site of use. On farms, small wind turbines were ideal for this purpose and in Denmark Poul la Cour, who was among the first to connect a windmill to a generator, gave a course for ‘agricultural electricians’. An example of La Cour’s great foresight was that he installed in his school one of the first wind tunnels in the world in order to investigate rotor aerodynamics. Gradually, however, diesel engines and steam turbines took over the production of electricity and only during the two world wars, when the supply of fuel was scarce, did wind power flourish again.3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 2 2 Aerodynamics of Wind Turbines However, even after the Second World War, the development of more efficient wind turbines was still pursued in several countries such as Germany, the US, France, the UK and Denmark. In Denmark, this work was undertaken by Johannes Juul, who was an employee in the utility company SEAS and a former student of la Cour. In the mid 1950s Juul introduced what was later called the Danish concept by constructing the famous Gedser turbine, which had an upwind three-bladed, stall regulated rotor, connected to an AC asynchronous generator running with almost constant speed. With the oil crisis in 1973, wind turbines suddenly became interesting again for many countries that wanted to be less dependent on oil imports; many national research programmes were initiated to investigate the possibilities of utilizing wind energy. Large non-commercial prototypes were built to evaluate the economics of wind produced electricity and to measure the loads on big wind turbines. Since the oil crisis, commercial wind turbines have gradually become an important industry with an annual turnover in the 1990s of more than a billion US dollars per year. Since then this figure has increased by approximately 20 per cent a year. Why Use Wind Power? As already mentioned, a country or region where energy production is based on imported coal or oil will become more self-sufficient by using alternatives such as wind power. Electricity produced from the wind produces no CO 2 emissions and therefore does not contribute to the greenhouse effect. Wind energy is relatively labour intensive and thus creates many jobs. In remote areas or areas with a weak grid, wind energy can be used for charging batteries or can be combined with a diesel engine to save fuel whenever wind is available. Moreover, wind turbines can be used for the desalination of water in coastal areas with little fresh water, for instance the Middle East. At windy sites the price of electricity, measured in /kWh, is competitive with the production price from more conventional methods, for example coal fired power plants. To reduce the price further and to make wind energy more competitive with other production methods, wind turbine manufacturers are concen- trating on bringing down the price of the turbines themselves. Other factors, such as interest rates, the cost of land and, not least, the amount of wind available at a certain site, also influence the production price of the electrical energy generated. The production price is computed as the investment plus the discounted maintenance cost divided by the discounted production measured in kWh over a period of typically 20 years. When the character-3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 3 General Introduction to Wind Turbines3 istics of a given turbine – the power for a given wind speed, as well as the annual wind distribution – are known, the annual energy production can be estimated at a specific site. Some of the drawbacks of wind energy are also mentioned. Wind turbines create a certain amount of noise when they produce electricity. In modern wind turbines, manufacturers have managed to reduce almost all mechanical noise and are now working hard on reducing aerodynamic noise from the rotating blades. Noise is an important competition factor, especially in densely populated areas. Some people think that wind turbines are unsightly in the landscape, but as bigger and bigger machines gradually replace the older smaller machines, the actual number of wind turbines will be reduced while still increasing capacity. If many turbines are to be erected in a region, it is important to have public acceptance. This can be achieved by allowing those people living close to the turbines to own a part of the project and thus share the income. Furthermore, noise and visual impact will in the future be less important as more wind turbines will be sited offshore. One problem is that wind energy can only be produced when nature supplies sufficient wind. This is not a problem for most countries, which are connected to big grids and can therefore buy electricity from the grid in the absence of wind. It is, however, an advantage to know in advance what resources will be available in the near future so that conventional power plants can adapt their production. Reliable weather forecasts are desirable since it takes some time for a coal fired power plant to change its product- ion. Combining wind energy with hydropower would be perfect, since it takes almost no time to open or close a valve at the inlet to a water turbine and water can be stored in the reservoirs when the wind is sufficiently strong. The Wind Resource A wind turbine transforms the kinetic energy in the wind to mechanical energy in a shaft and finally into electrical energy in a generator. The maxi- mum available energy, P , is thus obtained if theoretically the wind speed max · 2 3 · could be reduced to zero: P = 1/2 mV = 1/2 AV where m is the mass flow, o o V is the wind speed,  the density of the air and A the area where the wind o speed has been reduced. The equation for the maximum available power is very important since it tells us that power increases with the cube of the wind speed and only linearly with density and area. The available wind speed at a given site is therefore often first measured over a period of time before a project is initiated.3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 4 4 Aerodynamics of Wind Turbines In practice one cannot reduce the wind speed to zero, so a power coefficient C is defined as the ratio between the actual power obtained and p the maximum available power as given by the above equation. A theoretical maximum for C exists, denoted by the Betz limit, C = 16/27 = 0.593. p Pmax Modern wind turbines operate close to this limit, with C up to 0.5, and are p therefore optimized. Statistics have been given on many different turbines sited in Denmark and as rule of thumb they produce approximately 2 1000kWh/m /year. However, the production is very site dependent and the rule of thumb can only be used as a crude estimation and only for a site in Denmark. Sailors discovered very early on that it is more efficient to use the lift force than simple drag as the main source of propulsion. Lift and drag are the components of the force perpendicular and parallel to the direction of the relative wind respectively. It is easy to show theoretically that it is much more efficient to use lift rather than drag when extracting power from the wind. All modern wind turbines therefore consist of a number of rotating blades looking like propeller blades. If the blades are connected to a vertical shaft, the turbine is called a vertical-axis machine, VAWT, and if the shaft is horizontal, the turbine is called a horizontal-axis wind turbine, HAWT. For commercial wind turbines the mainstream mostly consists of HAWTs; the following text therefore focuses on this type of machine. A HAWT as sketched in Figure 1.1 is described in terms of the rotor diameter, the number of blades, the tower height, the rated power and the control strategy. D H Figure 1.1 Horizontal-axis wind turbine (HAWT)3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 5 General Introduction to Wind Turbines5 The tower height is important since wind speed increases with height above the ground and the rotor diameter is important since this gives the area A in the formula for the available power. The ratio between the rotor diameter D and the hub height H is often approximately one. The rated power is the maximum power allowed for the installed generator and the control system must ensure that this power is not exceeded in high winds. The number of blades is usually two or three. Two-bladed wind turbines are cheaper since they have one blade fewer, but they rotate faster and appear more flickering to the eyes, whereas three-bladed wind turbines seem calmer and therefore less disturbing in a landscape. The aerodynamic efficiency is lower on a two- bladed than on a three-bladed wind turbine. A two-bladed wind turbine is often, but not always, a downwind machine; in other words the rotor is downwind of the tower. Furthermore, the connection to the shaft is flexible, the rotor being mounted on the shaft through a hinge. This is called a teeter mechanism and the effect is that no bending moments are transferred from the rotor to the mechanical shaft. Such a construction is more flexible than the stiff three-bladed rotor and some components can be built lighter and smaller, which thus reduces the price of the wind turbine. The stability of the more flexible rotor must, however, be ensured. Downwind turbines are noisier than upstream turbines, since the once-per-revolution tower passage of each blade is heard as a low frequency noise. The rotational speed of a wind turbine rotor is approximately 20 to 50 rpm and the rotational speed of most generator shafts is approximately 1000 to 3000 rpm. Therefore a gearbox must be placed between the low-speed rotor shaft and the high-speed generator shaft. The layout of a typical wind turbine can be seen in Figure 1.2, showing a Siemens wind turbine designed for offshore use. The main shaft has two bearings to facilitate a possible replace- ment of the gearbox. This layout is by no means the only option; for example, some turbines are equipped with multipole generators, which rotate so slowly that no gearbox is needed. Ideally a wind turbine rotor should always be perpendicular to the wind. On most wind turbines a wind vane is therefore mounted somewhere on the turbine to measure the direction of the wind. This signal is coupled with a yaw motor, which continuously turns the nacelle into the wind. The rotor is the wind turbine component that has undergone the greatest development in recent years. The aerofoils used on the first modern wind turbine blades were developed for aircraft and were not optimized for the much higher angles of attack frequently employed by a wind turbine blade. Even though old aerofoils, for instance NACA63-4XX, have been used in the light of experience gained from the first blades, blade manufacturers have3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 6 6 Aerodynamics of Wind Turbines With permission from Siemens Wind Power. Figure 1.2 Machine layout now started to use aerofoils specifically optimized for wind turbines. Different materials have been tried in the construction of the blades, which must be sufficiently strong and stiff, have a high fatigue endurance limit, and be as cheap as possible. Today most blades are built of glass fibre reinforced plastic, but other materials such as laminated wood are also used. It is hoped that the historical review, the arguments for supporting wind power and the short description of the technology set out in this chapter will motivate the reader to study the more technical sections concerned with aerodynamics, structures and loads as applied to wind turbine construction.3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 7 2 2-D Aerodynamics Wind turbine blades are long and slender structures where the spanwise velocity component is much lower than the streamwise component, and it is therefore assumed in many aerodynamic models that the flow at a given radial position is two dimensional and that 2-D aerofoil data can thus be applied. Two-dimensional flow is comprised of a plane and if this plane is described with a coordinate system as shown in Figure 2.1, the velocity component in the z-direction is zero. In order to realize a 2-D flow it is necessary to extrude an aerofoil into a wing of infinite span. On a real wing the chord and twist changes along the span and the wing starts at a hub and ends in a tip, but for long slender wings, like those on modern gliders and wind turbines, Prandtl has shown that local 2-D data for the forces can be used if the angle of attack is corrected accordingly with the trailing vortices behind the wing (see, for example, Prandtl and Tietjens, 1957). These effects will be dealt with later, but it is now clear that 2-D aerodynamics is of practical interest even though it is difficult to realize. Figure 2.1 shows the leading edge stagnation point present in the 2-D flow past an aerofoil. The reacting force F from the flow is decomposed into a direction perpendicular to the velocity at infinity V and ∝ to a direction parallel to V . The former component is known as the lift, L; ∝ the latter is called the drag, D (see Figure 2.2). y z x Figure 2.1 Schematic view of streamlines past an airfoil3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 8 8 Aerodynamics of Wind Turbines Figure 2.2 Definition of lift and drag If the aerofoil is designed for an aircraft it is obvious that the L/D ratio should be maximized. The lift is the force used to overcome gravity and the higher the lift the higher the mass that can be lifted off the ground. In order to maintain a constant speed the drag must be balanced by a propulsion force delivered from an engine, and the smaller the drag the smaller the required engine. Lift and drag coefficients C and C are defined as: l d L C = ––––––––– (2.1) l 2 1 / ρV c 2 ∝ and: D C = ––––––––– (2.2) d 2 1 / ρV c 2 ∝ where ρ is the density and c the length of the aerofoil, often denoted by the chord. Note that the unit for the lift and drag in equations (2.1) and (2.2) is force per length (in N/m). A chord line can be defined as the line from the trailing edge to the nose of the aerofoil (see Figure 2.2). To describe the forces completely it is also necessary to know the moment M about a point in the aerofoil. This point is often located on the chord line at c/4 from the3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 9 2–D Aerodynamics9 Figure 2.3 Explanation of the generation of lift leading edge. The moment is positive when it tends to turn the aerofoil in Figure 2.2 clockwise (nose up) and a moment coefficient is defined as: M C = –––––––––– (2.3) m 2 2 1 / ρV c 2 ∝ The physical explanation of the lift is that the shape of the aerofoil forces the streamlines to curve around the geometry, as indicated in Figure 2.3. From basic fluid mechanics it is known that a pressure gradient, ∂p/∂r = 2 V /r, is necessary to curve the streamlines; r is the curvature of the stream- line and V the speed. This pressure gradient acts like the centripetal force known from the circular motion of a particle. Since there is atmospheric pressure p far from the aerofoil there must thus be a lower than atmospheric o pressure on the upper side of the aerofoil and a higher than atmospheric pressure on the lower side of the aerofoil. This pressure difference gives a lifting force on the aerofoil. When the aerofoil is almost aligned with the3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 10 10 Aerodynamics of Wind Turbines flow, the boundary layer stays attached and the associated drag is mainly caused by friction with the air. The coefficients C , C and C are functions of α, Re and Ma. α is the angle l d m of attack defined as the angle between the chordline and V; Re is the ∝ Reynolds number based on the chord and V, Re = cV /ν, where ν is the ∝ ∝ kinematic viscosity; and Ma denotes the Mach number, in other words the ratio between V and the speed of sound. For a wind turbine and a slow ∝ moving aircraft the lift, drag and moment coefficients are only functions of α and Re. For a given airfoil the behaviours of C , C and C are measured l d m or computed and plotted in so-called polars. An example of a measured polar for the FX67-K-170 airfoil is shown in Figure 2.4. C increases linearly with , with an approximate slope of 2/rad, until a l certain value of α, where a maximum value of C is reached. Hereafter the l aerofoil is said to stall and C decreases in a very geometrically dependent l manner. For small angles of attack the drag coefficient C is almost constant, d but increases rapidly after stall. The Reynolds number dependency can also be seen in Figure 2.4. It is seen, especially on the drag, that as the Reynolds number reaches a certain value, the Reynolds number dependency becomes small. The Reynolds number dependency is related to the point on the aerofoil, where the boundary layer transition from laminar to turbulent flow occurs. The way an aerofoil stalls is very dependent on the geometry. Thin aerofoils with a sharp nose, in other words with high curvature around the leading edge, tend to stall more abruptly than thick aerofoils. Different stall behaviours are seen in Figure 2.5, where C (α) is compared for two different l aerofoils. The FX38-153 is seen to lose its lift more rapidly than the FX67- K-170. The explanation lies in the way the boundary layer separates from the upper side of the aerofoil. If the separation starts at the trailing edge of the aerofoil and increases slowly with increasing angle of attack, a soft stall is observed, but if the separation starts at the leading edge of the aerofoil, the entire boundary layer may separate almost simultaneously with a dramatic loss of lift. The behaviour of the viscous boundary layer is very complex and depends, among other things, on the curvature of the aerofoil, the Reynolds number, the surface roughness and, for high speeds, also on the Mach number. Some description of the viscous boundary is given in this text but for a more elaborate description see standard textbooks on viscous boundary layers such as White (1991) and Schlichting (1968). Figure 2.6 shows the computed streamlines for a NACA63-415 aerofoil at angles of attack of 5° and 15°. For α = 15° a trailing edge separation is observed. The forces on the aerofoil stem from the pressure distribution p(x)3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 11 2–D Aerodynamics11 Figure 2.4 Polar for the FX67-K-170 airfoil3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 12 12 Aerodynamics of Wind Turbines Figure 2.5 Different stall behaviour Figure 2.6 Computed streamlines for angles of attack of 5° and 15°3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 13 2–D Aerodynamics13 and the skin friction with the air  = µ (∂u/∂y) . (x,y) is the surface w y=0 coordinate system as shown in Figure 2.7 and µ is the dynamic viscosity. The skin friction is mainly contributing to the drag, whereas the force found from integrating the pressure has a lift and drag component. The drag component from the pressure distribution is known as the form drag and becomes very large when the aerofoil stalls. The stall phenomenon is closely related to separation of the boundary layer (see next paragraph); therefore rule number one in reducing drag is to avoid separation. In Abbot and von Doenhoff (1959) a lot of data can be found for the National Advisory Committee for Aeronautics (NACA) aerofoils, which have been extensively used on small aircraft, wind turbines and helicopters. Close to the aerofoil there exists a viscous boundary layer due to the no-slip condition on the velocity at the wall (see Figure 2.7). A boundary layer thickness is often defined as the normal distance δ(x) from the wall where u(x)/U(x) = 0.99. Further, the displacement thickness δ (x), the momentum thickness θ(x) and the shape factor H(x) are defined as: δ u δ (x) = (1– —)dy, (2.4) ∫ 0 U δ u u θ(x)= — (1– —)dy, and (2.5) ∫ 0 U U δ . H(x)= – (2.6) θ The coordinate system (x,y) is a local system, where x = 0 is at the leading edge stagnation point and y is the normal distance from the wall. A turbulent boundary layer separates for H between 2 and 3. The stagnation streamline (see Figure 2.1) divides the fluid that flows over the aerofoil from the fluid that flows under the aerofoil. At the stagnation point the velocity is zero and the boundary layer thickness is small. The fluid which flows over the aerofoil accelerates as it passes the leading edge and, since the leading edge is close to the stagnation point and the flow accelerates, the boundary layer is thin. It is known from viscous boundary layer theory (see, for example, White, 1991) that the pressure is approximately constant from the surface to the edge of the boundary layer, i.e. ∂p/∂y = 0. Outside the boundary layer the Bernoulli equation (see Appendix A) is valid and, since the flow accelerates, the pressure decreases, i.e. ∂p/∂x0. On the lower side the pressure gradient is much smaller since the curvature of the wall is small compared to the leading edge. At the trailing edge the pressure must be the same at the upper and lower side (the Kutta condition) and therefore the pressure must rise, ∂p/∂x0, from a minimum value somewhere on the upper side to a higher3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 14 14 Aerodynamics of Wind Turbines Figure 2.7 Viscous boundary layer at the wall of an airfoil Figure 2.8 Schematic view of the shape of the boundary layer for a favourable and an adverse pressure gradient value at the trailing edge. An adverse pressure gradient, ∂p/∂x0, may lead to separation. This can be seen directly from the Navier-Stokes equations (see Appendix A) which applied at the wall, where the velocity is zero, reduces to: 2 ∂ u 1 ∂p –— = – — (2.7) 2 ∂y ∂x The curvature of the u-velocity component at the wall is therefore given by the sign of the pressure gradient. Further, it is known that ∂u/∂y = 0 at y = δ.3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 15 2–D Aerodynamics15 From this can be deduced that the u velocity profile in an adverse pressure gradient, ∂p/∂x 0, is S-shaped and separation may occur, whereas the curv- ature of the u velocity profile for ∂p/∂x 0 is negative throughout the entire boundary layer and the boundary layer stays attached. A schematic picture showing the different shapes of the boundary layer is given in Figure 2.8. Since the form drag increases dramatically when the boundary layer separ- ates, it is of utmost importance to the performance of an aerofoil to control the pressure gradient. For small values of x the flow is laminar, but for a certain x the laminar trans boundary layer becomes unstable and a transition from laminar to turbulent flow occurs. At x the flow is fully turbulent. In Figure 2.9 transition from a T laminar to a turbulent boundary layer is sketched. The transitional process is very complex and not yet fully understood, but a description of the phenomena is found in White (1991), where some engineering tools to predict x are also given. One of the models which is sometimes used in trans aerofoil computations is called the one-step method of Michel. The method predicts transition when: 0.4 Re = 2.9 Re (2.8) θ x where Re = U(x)⋅θ(x)/ν and Re = U(x)⋅ x/ν. For a laminar aerofoil (see later), θ x however, the Michel method might be inadequate and more advanced 9 methods such as the e method (see White, 1991) should be applied. Figure 2.9 Schematic view of the transitional process Turbulent flow is characterized by being more stable in regions of adverse pressure gradients, ∂p/∂x0, and by a steeper velocity gradient at the wall, ∂u/∂y . The first property is good since it delays stall, but the second y = 0 property increases the skin friction and thus the drag. These two phenomena3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 16 16 Aerodynamics of Wind Turbines are exploited in the design of high performance aerofoils called laminar aerofoils. A laminar aerofoil is an aerofoil where a large fraction of the boundary layer is laminar and attached in the range it is designed for. To design such an aerofoil it is necessary to specify the maximum angle of attack, where the boundary layer to a large extent is supposed to be laminar. The aerofoil is then constructed so that the velocity at the edge of the boundary layer, U(x), is constant after the acceleration past the leading edge and downstream. It is known from boundary layer theory (see White, 1991, and Schlichting, 1968) that the pressure gradient is expressed by the velocity outside the boundary layer as: dp dU(x) –– = –U(x) ––––– (2.9) dx dx At this angle the pressure gradient is therefore zero and no separation will occur. For smaller angles of attack the flow U(x) will accelerate and dp/dx becomes negative, which again avoids separation and is stabilizing for the laminar boundary layer, thus delaying transition. At some point x at the upper side of the aerofoil it is, however, necessary to decelerate the flow in order to fulfil the Kutta condition; in other words the pressure has to be unique at the trailing edge. If this deceleration is started at a position where the boundary layer is laminar, the boundary layer is likely to separate. Just after the laminar/ turbulent transition the boundary layer is relatively thin and the momentum close to the wall is relatively large and is therefore capable of withstanding a high positive pressure gradient without separation. During the continuous deceleration towards the trailing edge the ability of the boundary layer to withstand the positive pressure gradient diminishes, and to avoid separation it is therefore necessary to decrease the deceleration towards the trailing edge. It is of utmost importance to ensure that the boundary layer is turbulent before decelerating U(x). To ensure this, a turbulent transition can be triggered by placing a tripwire or tape before the point of deceleration. A laminar aerofoil is thus characterized by a high value of the lift to drag ratio C /C below the design l d angle. But before choosing such an aerofoil it is important to consider the stall characteristic and the roughness sensitivity. On an aeroplane it is necessary to fly with a high C at landing since the speed is relatively small. If the pilot l exceeds C and the aerofoil stalls, it could be disastrous if C drops as l ,max l drastically with the angle of attack as on the FX38-153 in Figure 2.5. The aeroplane would then lose its lift and might slam into the ground. If the aerofoil is sensitive to roughness, good performance is lost if the wings are contaminated by dust, rain particles or insects, for example. On a wind turbine this could alter the performance with time if, for instance, the3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 17 2–D Aerodynamics17 turbine is sited in an area with many insects. If a wind turbine is situated near the coast, salt might build up on the blades if the wind comes from the sea, and if the aerofoils used are sensitive to roughness, the power output from the turbine will become dependent on the direction of the wind. Fuglsang and Bak (2003) describe some attempts to design aerofoils specifically for use on wind turbines, where insensitivity to roughness is one of the design targets. To compute the power output from a wind turbine it is necessary to have data of C (α,Re) and C (α,Re) for the aerofoils applied along the blades. l d These data can be measured or computed using advanced numerical tools, but since the flow becomes unsteady and three-dimensional after stall, it is difficult to obtain reliable data for high angles of attack. On a wind turbine very high angles of attack may exist locally, so it is often necessary to extrapolate the available data to high angles of attack. References Abbot, H. and von Doenhoff, A. E. (1959) Theory of Wing Sections, Dover Publications, New York Fuglsang, P. and Bak, C. (2003) ‘Status of the Risø wind turbine aerofoils’, presented at the European Wind Energy Conference, EWEA, Madrid, 16–19 June Prandtl, L. and Tietjens, O. G. (1957) Applied Hydro and Aeromechanics, Dover Publications, New York Schlichting, H. (1968) Boundary-Layer Theory, McGraw-Hill, New York White, F. M. (1991) Viscous Fluid Flow, McGraw-Hill, New York3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 18 3 3-D Aerodynamics This chapter describes qualitatively the flow past a 3-D wing and how the spanwise lift distribution changes the upstream flow and thus the local angle of attack. Basic vortex theory, as described in various textbooks (for example Milne-Thomsen, 1952), is used. Since this theory is not directly used in the Blade Element Momentum method derived later, it is only touched on very briefly here. This chapter may therefore be quite abstract for the reader with limited knowledge of vortex theory, but hopefully some of the basic results will be qualitatively understood. A wing is a beam of finite length with aerofoils as cross-sections and therefore a pressure difference between the lower and upper sides is created, giving rise to lift. At the tips are leakages, where air flows around the tips from the lower side to the upper side. The streamlines flowing over the wing will thus be deflected inwards and the streamlines flowing under the wing The wing is seen from the suction side. The streamline flowing over the suction side (full line) is deflected inwards and the streamline flowing under (dashed line) is deflected outwards. Figure 3.1 Streamlines flowing over and under a wing3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 19 3–D Aerodynamics19 will be deflected outwards. Therefore at the trailing edge there is a jump in the tangential velocity (see Figures 3.1 and 3.2). A jump in the tangential velocity is seen, due to the leakage at the tips. Figure 3.2 Velocity vectors seen from behind a wing Because of this jump there is a continuous sheet of streamwise vorticity in the wake behind a wing. This sheet is known as the trailing vortices. In classic literature on theoretical aerodynamics (see, for example, Milne- Thomsen, 1952), it is shown that a vortex filament of strength Γ can model the flow past an aerofoil for small angles of attack. This is because the flow for small angles of attack is mainly inviscid and governed by the linear Laplace equation. It can be shown analytically that for this case the lift is given by the Kutta-Joukowski equation: L = ρV Γ. (3.1) ∝ An aerofoil may be thus substituted by one vortex filament of strength Γ and the lift produced by a 3-D wing can be modelled for small angles of attack by a series of vortex filaments oriented in the spanwise direction of the wing, known as the bound vortices. According to the Helmholtz theorem (Milne- Thomsen, 1952), a vortex filament, however, cannot terminate in the interior of the fluid but must either terminate on the boundary or be closed. A complete wing may be modelled by a series of vortex filaments, Γ, i = i 1,2,3,4,..., which are oriented as shown in Figure 3.3. In a real flow the trailing vortices will curl up around the strong tip vortices and the vortex system will look more like that in Figure 3.4.3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 20 20 Aerodynamics of Wind Turbines Figure 3.3 Asimplified model of the vortex system on a wing The model based on discrete vortices, as shown in Figure 3.3, is called the lifting line theory (see Schlichting and Truckenbrodt, 1959 for a complete description). The vortices on the wing (bound vortices) model the lift, and the trailing vortices (free vortices) model the vortex sheet stemming from the three dimensionality of the wing. The free vortices induce by the Biot-Savart law a downwards velocity component at any spanwise position of the wing. For one vortex filament of strength Γ the induced velocity at a point p is (see Figure 3.5): Γ r ds   ∫ w =— ——— . (3.2) 3 4 r Figure 3.4 More realistic vortex system on a wing3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 21 3–D Aerodynamics21 Figure 3.5 Induced velocity from a vortex line of strength Γ The total induced velocity from all vortices at a section of the wing is known as the downwash, and the local angle of attack at this section is therefore reduced by α , since the relative velocity is the vector sum of the wind speed i V and the induced velocity w. α , α and α denote the geometric, the induced ∝ g i e and the effective angles of attack respectively. The effective angle of attack is thus: α = α – α . (3.3) e g i In Figure 3.6 the induced velocity w, the onset flow V and the effective ∝ velocity V are shown for a section on the wing together with the different e angles of attack α , α and α . It is assumed that equation (3.1) is also valid g i e for a section in a 3-D wing if the effective velocity is used. The local lift force R, which is perpendicular to the relative velocity, is shown in Figure 3.6. The global lift is by definition the force perpendicular to the onset flow V and the ∞ resulting force, R, must therefore be decomposed into components perpendicular to and parallel to the direction of V . The former component is ∞ thus the lift and the latter is a drag denoted by induced drag D . At the tips of i the wing the induced velocity obtains a value which exactly ensures zero lift. An important conclusion is thus: For a three-dimensional wing the lift is reduced compared to a two- dimensional wing at the same geometric angle of attack, and the local lift has a component in the direction of the onset flow, which is known as the induced3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 22 22 Aerodynamics of Wind Turbines Figure 3.6 The effective angle of attack for a section in a wing and the resulting force R, lift L and induced drag D i drag. Both effects are due to the downwash induced by the vortex system of a 3-D wing. In the lifting line theory it is assumed that the three-dimensionality is limited to the downwash, in other words that the spanwise flow is still small com- pared to the streamwise velocity and 2-D data can therefore be used locally if the geometric angle of attack is modified by the downwash. This assump- tion is reasonable for long slender wings such as those on a glider or a wind turbine. One method to determine the value of the vortices quantitatively and thus the induced velocities is Multhopp’s solution of Prandtl’s integral equation. This method is thoroughly described in, for example, Schlichting and Truckenbrodt (1959) and will not be shown here, but it is important to understand that the vortex system produced by a three-dimensional wing changes the local inflow conditions seen by the wing, in other words that although the flow is locally 2-D one cannot apply the geometric angle of attack when estimating the forces on the wing. This error was made in early propeller theory and the discrepancy between measured and computed performance was believed to be caused by wrong 2-D aerofoil data. On a3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 23 3–D Aerodynamics23 rotating blade Coriolis and centrifugal forces play an important role in the separated boundary layers which occur after stall. In a separated boundary layer the velocity and thus the momentum is relatively small compared to the centrifugal force, which therefore starts to pump fluid in the spanwise direction towards the tip. When the fluid moves radially towards the tip the Coriolis force points towards the trailing edge and acts as a favourable pressure gradient. The effect of the centrifugal and Coriolis force is to alter the 2-D aerofoil data after stall. Considerable engineering skill and exper- ience is required to construct such post-stall data – for example to compute the performance of a wind turbine at high wind speeds – in order to obtain an acceptable result (see also Snel et al, 1993 and Chaviaropoulos and Hansen, 2000). Figure 3.7 shows the computed limiting streamlines on a modern wind turbine blade at a moderately high wind speed (Hansen et al, 1997). Limiting streamlines are the flow pattern very close to the surface. Figure 3.7 shows that for this specific blade at a wind speed of 10m/s the flow is attached on the outer part of the blade and separated at the inner part, where the limiting streamlines have a spanwise component. Figure 3.7 Computed limiting streamlines on a stall regulated wind turbine blade at a moderately high wind speed Vortex System behind a Wind Turbine The rotor of a horizontal-axis wind turbine consists of a number of blades, which are shaped as wings. If a cut is made at a radial distance, r, from the rotational axis as shown in Figure 3.8, a cascade of aerofoils is observed as shown in Figure 3.9. The local angle of attack α is given by the pitch of the aerofoil θ; the axial velocity and rotational velocity at the rotor plane denoted respectively by V a and V (see Figure 3.9): rot α = φ – θ (3.4)3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 24 24 Aerodynamics of Wind Turbines Figure 3.8 Rotor of a three-bladed wind turbine with rotor radius R Figure 3.9 Radial cut in a wind turbine rotor showing airfoils at r/R3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 25 3–D Aerodynamics25 Source: Wilson and Lissaman (1974), reproduced with permission Figure 3.10 Schematic drawing of the vortex system behind a wind turbine where the flow angle φ is found as: V a tan φ = –— (3.5) V rot Since a horizontal-axis wind turbine consists of rotating blades, a vortex system similar to the linear translating wing must exist. The vortex sheet of the free vortices is oriented in a helical path behind the rotor. The strong tip vortices are located at the edge of the rotor wake and the root vortices lie mainly in a linear path along the axis of the rotor, as shown in Figure 3.10. The vortex system induces on a wind turbine an axial velocity component opposite to the direction of the wind and a tangential velocity component opposite to the rotation of the rotor blades. The induced velocity in the axial direction is specified through the axial induction factor a as aV, where V is o o the undisturbed wind speed. The induced tangential velocity in the rotor wake is specified through the tangential induction factor a’ as 2a’ω r. Since the flow does not rotate upstream of the rotor, the tangential induced velocity in the rotor plane is thus approximately a’ω r. ω denotes the angular velocity of the rotor and r is the radial distance from the rotational axis. If a and a’ are known, a 2-D equivalent angle of attack could be found from equations (3.4) and (3.5), where: V = (1– a)V , (3.6) a o3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 26 26 Aerodynamics of Wind Turbines and: V = (1+ a’)ωr. (3.7) rot Furthermore, if the lift and drag coefficients C (α) and C (α) are also known l d for the aerofoils applied along the blades, it is easy to compute the force distribution. Global loads such as the power output and the root bending moments of the blades are found by integrating this distribution along the span. It is the purpose of the Blade Element Momentum method, which will later be derived in detail, to compute the induction factors a and a’ and thus the loads on a wind turbine. It is also possible to use a vortex method and construct the vortex system as shown in Figure 3.10 and use the Biot-Savart equation (3.2) to calculate the induced velocities. Such methods are not derived in this book but can found in, for example, Katz and Plotkin (2001) or Leishman (2006). References Chaviaropoulos, P. K. and Hansen, M. O. L. (2000) ‘Investigating three-dimensional and rotational effects on wind turbine blades by means of a quasi-3D Navier-Stokes solver’, Journal of Fluids Engineering, vol 122, pp330–336 Hansen, M. O. L, Sorensen, J. N., Michelsen, J. A. and Sorensen, N. N. (1997) ‘A global Navier-Stokes rotor prediction model’, AIAA 97-0970 paper, 35th Aerospace Sciences Meeting and Exhibition, Reno, Nevada, 6–9 January Katz, J. and Plotkin, A. (2001) Low-Speed Aerodynamics, Cambridge University Press, Cambridge Leishmann, J. G. (2006) Principles of Helicopter Aerodynamics, Cambridge University Press, Cambridge Milne-Thomson, L. M. (1952) Theoretical Aerodynamics, Macmillan, London Schlichting, H. and Truckenbrodt, E. (1959) Aerodynamik des Flugzeuges, Springer-Verlag, Berlin Snel, H., Houwink, B., Bosschers, J., Piers, W. J., van Bussel, G. J. W. and Bruining, A. (1993) ‘Sectional prediction of 3-D effects for stalled flow on rotating blades and comparison with measurements’, in Proceedings of European Community Wind Energy Conference 1993, H. S. Stephens & Associates, Travemunde, pp395–399 Wilson, R. E. and Lissaman, P. B. S. (1974) Applied Aerodynamics of Wind Power Machines, Technical Report NSF-RA-N-74-113, Oregon State University3212 J&J Aerodynamic Turbines 15/11/07 1:42 PM Page 27 4 1-D Momentum Theory for an Ideal Wind Turbine Before deriving the Blade Element Momentum method it is useful to examine a simple one-dimensional (1-D) model for an ideal rotor. A wind turbine extracts mechanical energy from the kinetic energy of the wind. The rotor in this simple 1-D model is a permeable disc. The disc is considered ideal; in other words it is frictionless and there is no rotational velocity component in the wake. The latter can be obtained by applying two contra-rotating rotors or a stator. The rotor disc acts as a drag device slowing the wind speed from V far upstream of the rotor to u at the rotor plane and to o u in the wake. Therefore the streamlines must diverge as shown in Figure 1 4.1. The drag is obtained by a pressure drop over the rotor. Close upstream of the rotor there is a small pressure rise from the atmospheric level p to p o before a discontinuous pressure drop ∆ p over the rotor. Downstream of the rotor the pressure recovers continuously to the atmospheric level. The Mach number is small and the air density is thus constant and the axial velocity must decrease continuously from V to u . The behaviour of the pressure and o 1 axial velocity is shown graphically in Figure 4.1. Using the assumptions of an ideal rotor it is possible to derive simple relationships between the velocities V, u and u, the thrust T, and the absorbed o 1 shaft power P. The thrust is the force in the streamwise direction resulting from the pressure drop over the rotor, and is used to reduce the wind speed from V to u : o 1 T = ∆ pA, (4.1) 2 where A= πR is the area of the rotor. The flow is stationary, incompressible and frictionless and no external force acts on the fluid up- or downstream of the rotor. Therefore the Bernoulli equation (see Appendix A) is valid from far upstream to just in front of the rotor and from just behind the rotor to far downstream in the wake:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 28 28 Aerodynamics of Wind Turbines Figure 4.1 Illustration of the streamlines past the rotor and the axial velocity and pressure up- and downstream of the rotor 1 1 2 2 p + —V = p + — u , (4.2) o o 2 2 and: 1 1 2 2 p – ∆ p + — u = p + — u . (4.3) o 1 2 2 Combining equation (4.2) and (4.3) yields: 1 2 2 ∆ p = — (V – u ). (4.4) o 1 2 The axial momentum equation in integral form (see Appendix A) is applied on the circular control volume with sectional area A drawn with a dashed cv line in Figure 4.2 yielding: ∂ — u(x, y, z)dxdydz + ∫∫ u(x, y, z)V·dA = F + F . (4.5) ∫∫∫ cv cs ext pres ∂t3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 29 1-D Momentum Theory for an Ideal Wind Turbine29 dA is a vector pointing outwards in the normal direction of an infinitesimal part of the control surface with a length equal to the area of this element. F pres is the axial component of the pressure forces acting on the control volume. The first term in equation (4.5) is zero since the flow is assumed to be stationary and the last term is zero since the pressure has the same atmospheric value on the end planes and acts on an equal area. Further, on the lateral boundary of the control volume shown in Figure 4.2, the force from the pressure has no axial component. Figure 4.2 Circular control volume around a wind turbine Using the simplified assumptions of an ideal rotor, equation (4.5) then yields: 2 2 · 2 u A + V (A – A ) + m V – V A = –T. (4.6) 1 1 o cv 1 side o o cv · m can be found from the conservation of mass: side · A u + (A – A )V + m =  A V . (4.7) 1 1 1 o side cv o cv yielding: · m =  A (V – u ). (4.8) side 1 o 1 The conservation of mass also gives a relationship between A and A as: 1 · m = uA = u A . (4.9) 1 13212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 30 30 Aerodynamics of Wind Turbines Combining equations (4.8), (4.9) and (4.6) yields: · T = uA(V – u ) = m (V – u ). (4.10) o 1 o 1 If the thrust is replaced by the pressure drop over the rotor as in equation (4.1) and the pressure drop from equation (4.4) is used, an interesting observation is made: 1 u = – (V + u ). (4.11) o 1 2 It is seen that the velocity in the rotor plane is the mean of the wind speed V o and the final value in the wake u . 1 An alternative control volume to the one in Figure 4.2 is shown in Figure 4.3. Figure 4.3 Alternative control volume around a wind turbine The force from the pressure distribution along the lateral walls F of the press, lateral control volume is unknown and thus so is the net pressure contribution F . pres On this alternative control volume there is no mass flow through the lateral boundary, since this is aligned with the streamlines. The axial momentum equation (4.5) therefore becomes: T = uA(V – u ) + F . (4.12) o 1 pres Since the physical problem is the same, whether the control volume in Figure 4.2 or that in Figure 4.3 is applied, it can be seen by comparing equations3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 31 1-D Momentum Theory for an Ideal Wind Turbine31 (4.10) and (4.12) that the net pressure force on the control volume following the streamlines is zero. The flow is assumed to be frictionless and there is therefore no change in the internal energy from the inlet to the outlet and the shaft power P can be found using the integral energy equation on the control volume shown in Figure 4.3: p p · 1 2 o 1 2 o P= m (– V + — – – u – —). (4.13) o 1   2 2 · and since m = ρuA the equation for P simply becomes: 1 2 2 P= – ρuA(V – u ). (4.14) o 1 2 The axial induction factor a is defined as: u = (1 – a)V. (4.15) o Combining equation (4.15) with (4.11) gives: u = (1 – 2a)V, (4.16) 1 o which then can be introduced in equation (4.14) for the power P and into equation (4.10) for the thrust T, yielding: 3 2 P = 2ρV a(1 – a) A (4.17) o and: 2 T = 2ρV a(1 – a)A. (4.18) o The available power in a cross-section equal to the swept area A by the rotor is: 1 3 P = – ρAV (4.19) avail o 2 The power P is often non-dimensionalized with respect to P as a power avail coefficient C : p P C = —–—– (4.20) p 3 1 ρV A – o 23212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 32 32 Aerodynamics of Wind Turbines Similarly a thrust coefficient C is defined as: T T C = —–—– (4.21) T 2 1 ρV A – o 2 Using equations (4.17) and (4.18) the power and thrust coefficients for the ideal 1-D wind turbine may be written as: 2 C = 4a(1 – a) (4.22) p and: C = 4a(1 – a). (4.23) T Differentiating C with respect to a yields: p dC P — — = 4(1 – a)(1 – 3a). (4.24) da It is easily seen that Cp, max = 16/27 for a = 1/3. Equations (4.22) and (4.23) are shown graphically in Figure 4.4. This theoretical maximum for an ideal wind turbine is known as the Betz limit. Figure 4.4 The power and thrust coefficients C and C as a function of the axial induction p T factor a for an ideal horizontal-axis wind turbine3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 33 1-D Momentum Theory for an Ideal Wind Turbine33 Experiments have shown that the assumptions of an ideal wind turbine leading to equation (4.23) are only valid for an axial induction factor, a, of less than approximately 0.4. This is seen in Figure 4.5, which shows measurements of C as a function of a for different rotor states. If the T momentum theory were valid for higher values of a, the velocity in the wake would become negative as can readily be seen by equation (4.16). Source: Eggleston and Stoddard (1987), reproduced with permission. Figure 4.5 The measured thrust coefficient C as a function of the axial induction factor T a and the corresponding rotor states As C increases the expansion of the wake increases and thus also the T velocity jump from V to u in the wake, see Figure 4.6. o 1 The ratio between the areas A and A in Figure 4.6 can be found directly o 1 from the continuity equation as: A o — — = 1 – 2a. (4.25) A 1 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 34 34 Aerodynamics of Wind Turbines Figure 4.6 The expansion of the wake and the velocity jump in the wake for the 1-D model of an ideal wind turbine For a wind turbine, a high thrust coefficient C , and thus a high axial T induction factor a, is present at low wind speeds. The reason that the simple momentum theory is not valid for values of a greater than approximately 0.4 is that the free shear layer at the edge of the wake becomes unstable when the velocity jump V – u becomes too high and eddies are formed which transport o 1 momentum from the outer flow into the wake. This situation is called the turbulent-wake state, see Figures 4.5 and 4.7.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 35 1-D Momentum Theory for an Ideal Wind Turbine35 Figure 4.7 Schematic view of the turbulent-wake state induced by the unstable shear flow at the edge of the wake Effects of Rotation For the ideal rotor there is no rotation in the wake; in other words a’ is zero. Since a modern wind turbine consists of a single rotor without a stator, the wake will possess some rotation as can be seen directly from Euler’s turbine equation (see Appendix A) applied to an infinitesimal control volume of thickness dr, as shown in Figure 3.8: · 2 dP = mωrC = 2r uωC dr, (4.26) θ θ where C is the azimuthal component of the absolute velocity C = (C ,Cθ,C ) r a θ after the rotor and u the axial velocity through the rotor. Since the forces felt by the wind turbine blades are also felt by the incoming air, but with opposite sign, the air at a wind turbine will rotate in the opposite direction from that of the blades. This can also be illustrated using Figure 4.8, where the relative velocity upstream of the blade V is rel,1 given by the axial velocity u and the rotational velocity V . For moderate rot angles of attack the relative velocity V downstream of the rotor rel,2 approximately follows the trailing edge. The axial component, C , of the a absolute velocity equals u due to conservation of mass, and the rotational3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 36 36 Aerodynamics of Wind Turbines speed is unaltered. The velocity triangle downstream of the blade is now fixed and, as Figure 4.8 shows, the absolute velocity downstream of the blade, C, has a tangential component C in the opposite direction of the blade. θ Figure 4.8 The velocity triangle for a section of the rotor From equation (4.26) it is seen that for a given power P and wind speed the azimuthal velocity component in the wake Cθ decreases with increasing rotational speed ω of the rotor. From an efficiency point of view it is therefore desirable for the wind turbine to have a high rotational speed to minimize the loss of kinetic energy contained in the rotating wake. If we recall that the axial velocity through the rotor is given by the axial induction factor a as in equation (4.15) and that the rotational speed in the wake is given by a’ as: C = 2a’ωr. (4.27) θ Equation (4.26) may then be written as: 2 3 dP = 4ρω Va’(1 – a)r dr. (4.28) o The total power is found by integrating dP from 0 to R as: 2 R 3 P = 4ρω V ∫ a’(1 – a)r dr. (4.29) o 03212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 37 1-D Momentum Theory for an Ideal Wind Turbine37 or in non-dimensional form as: 8 λ 3 C = –– ∫ a’(1 – a)x dx, (4.30) p 0 2 λ where λ = ωR/V is the tip speed ratio and x = ωr/V is the local rotational o o speed at the radius r non-dimensionalized with respect to the wind speed V. o It is clear from equations (4.29) and (4.30) that in order to optimize the power it is necessary to maximize the expression: f (a, a’) = a’(1 – a). (4.31) If the local angles of attack are below stall, a and a’ are not independent since the reacting force according to potential flow theory is perpendicular to the local velocity seen by the blade as indicated by equation (3.1). The total induced velocity, w, must be in the same direction as the force and thus also perpendicular to the local velocity. The following relationship therefore exists between a and a’: 2 x a’(1 + a’) = a(1 – a). (4.32) Equation (4.32) is directly derived from Figure 4.9 since: a’ωr tan φ = ––––– (4.33) aV o and: (1 – a)V o tan φ = –––– –––––– (4.34) (1 + a’)ωr x = ωr/V denotes the ratio between the local rotational speed and the wind o speed. For local angles of attack below stall a and a’ are related through equation (4.32) and the optimization problem is thus to maximize equation (4.31) and still satisfy equation (4.32). Since a’ is a function of a, the expression (4.31) is maximum when df/da = 0 yielding: df da’ — = (1 – a) — – a’ = 0, (4.35) da da3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 38 38 Aerodynamics of Wind Turbines Note that for small angles of attack the total induced velocity w is perpendicular to the relative velocity. Figure 4.9 Velocity triangle showing the induced velocities for a section of the blade which can be simplified to: da’ (1 – a) — = a’ (4.36) da Equation (4.32) differentiated with respect to a yields: da’ 2 (1 + 2a’) — x = 1 – 2a. (4.37) da If equations (4.36) and (4.37) are combined with equation (4.32), the opti- mum relationship between a and a’ becomes: 1–3a a’ = ––––– . (4.38) 4a – 1 A table between a, a’ and x can now be computed. a’ is given by equation (4.38) for a specified a and then x is found using (4.32). It can be seen that as the rotational speed ω and thus also x = ωr/V is o increased the optimum value for a tends to 1/3, which is consistent with the simple momentum theory for an ideal rotor. Using the values from the table,3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 39 1-D Momentum Theory for an Ideal Wind Turbine39 the optimum power coefficient C is found by integrating equation (4.30). p This is done in Glauert (1935) for different tip speed ratios λ = ωR/V. o Glauert compares this computed optimum power coefficient with the Betz limit of 16/27, which is derived for zero rotation in the wake a’= 0 (see Table 4.2). In Figure 4.10, Table 4.2 is plotted and it can be seen that the loss due to rotation is small for tip speed ratios greater than approximately 6. Table 4.1 The numerical relationships between a, a’ and x aa’ x 0.26 5.5 0.073 0.27 2.375 0.157 0.28 1.333 0.255 0.29 0.812 0.374 0.30 0.500 0.529 0.31 0.292 0.753 0.320.143 1.15 0.33 0.031 2.63 0.333 0.00301 8.58 Table 4.2 Glauert’s comparison of the computed optimum power coefficient including wake rotation with the Betz limit λ = ωR/V 27C /16 o p 0.5 0.486 1.0 0.703 1.5 0.811 2.0 0.865 2.5 0.899 5.0 0.963 7.5 0.983 10.0 0.9873212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 40 40 Aerodynamics of Wind Turbines The efficiency is defined as the ratio between C , including wake rotation, and the Betz p limit C = 16/27. p, Betz Figure 4.10 The efficiency of an optimum turbine with rotation References Eggleston, D. M. and Stoddard, F. S. (1987) Wind Turbine Engineering Design, Van Nostrand Reinhold Company, New York Glauert, H. (1935) ‘Airplane propellers’, in W. F. Durand (ed) Aerodynamic Theory, vol 4, Division L, Julius Springer, Berlin, pp169–3603212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 41 5 Shrouded Rotors It is possible to exceed the Betz limit by placing the wind turbine in a diffuser. If the cross-section of the diffuser is shaped like an aerofoil, a lift force will be generated by the flow through the diffuser as seen in Figure 5.1. Figure 5.1 Ideal flow through a wind turbine in a diffuser As shown in de Vries (1979), the effect of this lift is to create a ring vortex, which by the Biot-Savart law will induce a velocity to increase the mass flow through the rotor. The axial velocity in the rotor plane is denoted by V, and 2 ε is the augmentation defined as the ratio between V and the wind speed V, 2 o i.e. ε = V /V. A 1-D analysis of a rotor in a diffuser gives the following 2 o expression for the power coefficient: P T·V 2 C = –––––––– = ––––––––––– = C ε. (5.1) p,d T 1 3 V – ρV A o o 2 2 1 –– ρV V A – o 2 2 V 23212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 42 42 Aerodynamics of Wind Turbines For an ideal bare turbine equations (4.22) and (4.23) are valid yielding: C = C (1 – a). (5.2) p,b T Combining equations (5.1) and (5.2) yields: C p,d ε ––– – = ––––– – (5.3) C p,b (1 – a) Further, the following equations are valid for the mass flow through a bare · · turbine m and the mass flow through a turbine in a diffuser m : b d · m (1 – a)V A o b –– ––– = = 1 – a (5.4) VA VA o o · m V A 2 d –––– – – = ––––– – = ε. (5.5) V A VA o o Combining equations (5.3), (5.4) and (5.5) yields: · C p,d m d ––– – = ––– (5.6) · m C p,b b Equation (5.6) states that the relative increase in the power coefficient for a shrouded turbine is proportional to the ratio between the mass flow through the turbine in the diffuser and the same turbine without the diffuser. Equation (5.6) is verified by computational fluid dynamics (CFD) results as seen in · · Figure 5.2, where for a given geometry the computed mass flow ratio m /m d b is plotted against the computed ratio C /C . The CFD analysis is done on a p,d p,b simple geometry, without boundary layer bleed slots, as suggested by Gilbert and Foreman (1983). The diffuser was modelled using 266,240 grid points with 96 points around the diffuser aerofoil section, and a turbulence model was chosen which is sensitive to adverse pressure gradients (see Hansen et al, 2000). The rotor was modelled by specifying a constant volume force at the position of the rotor. To check this approach some initial computations were made without the diffuser and in Figure 5.3, which shows the relationship between the thrust and power coefficients, it is seen that this approach gave good results compared to the following theoretical expression, which can be derived from equations (4.22) and (4.23): —— — 1 C = —C (1 + 1 – C ) (5.7) p,b T T 23212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 43 Shrouded Rotors43 Figure 5.2 Computed mass flow ratio plotted against the computed power coefficient ratio for a bare and a shrouded wind turbine, respectively In Figure 5.3 it is also seen that computations with a wind turbine in a diffuser gave higher values for the power coefficient than the Betz limit for a bare turbine. The theoretical relationship equation (5.7) for a bare rotor is also compared in this figure. Figure 5.3 Computed power coefficient for a rotor in a diffuser as a function of the thrust coefficient C T3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 44 44 Aerodynamics of Wind Turbines The results are dependent on the actual diffuser geometry, in other words on the amount of lift which can be generated by the diffuser. An adverse pressure gradient is present for the flow in the diffuser, and the boundary layer will separate if the ratio between the exit area and the area in the diffuser becomes too high. To increase the lift, giving a higher mass flow through the turbine and thus a higher power output, any trick to help prevent the boundary layer from separating is allowed, for example vortex generators or boundary bleed slots. The computations in Figure 5.3 and the wind tunnel measurements of Gilbert and Foreman (1983) show that the Betz limit can be exceeded if a device increasing the mass flow through the rotor is applied, but this still has to be demonstrated on a full size machine. Furthermore, the increased energy output has to be compared to the extra cost of building a diffuser and the supporting structure. References Gilbert, B. L. and Foreman, K. M. (1983) ‘Experiments with a diffuser-augmented model wind turbine’, Journal of Energy Resources Technology, vol 105, pp46–53 Hansen, M. O. L, Sorensen, N. N. and Flay, R. G. J. (2000) ‘Effect of placing a diffuser around a wind turbine’, Wind Energy, vol 3, pp207–213 de Vries, O. (1979) Fluid Dynamic Aspects of Wind Energy Conversion, AGARDograph No 243, Advisory Group for Aeronautical Research and Development3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 45 6 The Classical Blade Element Momentum Method All definitions and necessary theory to understand the Blade Element Momentum (BEM) method have now been introduced. In this chapter the classical BEM model from Glauert (1935) will be presented. With this model it is possible to calculate the steady loads and thus also the thrust and power for different settings of wind speed, rotational speed and pitch angle. To calculate time series of the loads for time-varying input some engineering models must be added, as will be shown in a later chapter. In the 1-D momentum theory the actual geometry of the rotor – the number of blades, the twist and chord distribution, and the aerofoils used – is not considered. The Blade Element Momentum method couples the momentum theory with the local events taking place at the actual blades. The stream tube introduced in the 1-D momentum theory is discretized into N annular elements of height dr, as shown in Figure 6.1. The lateral boundary of these elements consists of streamlines; in other words there is no flow across the elements. Figure 6.1 Control volume shaped as an annular element to be used in the BEM model3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 46 46 Aerodynamics of Wind Turbines In the BEM model the following is assumed for the annular elements: 1 No radial dependency – in other words what happens at one element cannot be felt by the others. 2The force from the blades on the flow is constant in each annular element; this corresponds to a rotor with an infinite number of blades. A correction known as Prandtl’s tip loss factor is later introduced to correct for the latter assumption in order to compute a rotor with a finite number of blades. In the previous section concerning the 1-D momentum theory it was proven that the pressure distribution along the curved streamlines enclosing the wake does not give an axial force component. Therefore it is assumed that this is also the case for the annular control volume shown in Figure 6.1. The thrust from the disc on this control volume can thus be found from the integral momentum equation since the cross-section area of the control volume at the rotor plane is 2πrdr: · dT = (V – u )dm = 2ru(V – u )dr. (6.1) o 1 o 1 The torque dM on the annular element is found using the integral moment of momentum equation on the control volume (see Appendix A) and setting the rotational velocity to zero upstream of the rotor and to C in the wake: θ 2 · dM = rC dm = 2r uC dr. (6.2) θ θ This could also have been derived directly from Euler’s turbine equation (4.26), since: dP = ωdM (6.3) From the ideal rotor it was found that the axial velocity in the wake u1 could be expressed by the axial induction factor a and the wind speed V as u = (1 o 1 – 2a)V, and if this is introduced into equations (6.1) and (6.2) together with o the definitions for a and a’ in equations (4.15) and (4.27) the thrust and torque can be computed as: 2 dT = 4rV a(1 – a)dr (6.4) o and:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 47 The Classical Blade Element Momentum Method47 3 dM = 4r Vω(1 – a)a’dr. (6.5) o The left hand sides of equations (6.4) and (6.5) are found from the local flow around the blade. It is recalled that the relative velocity V seen by a section rel of the blade is a combination of the axial velocity (1 – a)V and the tangential o velocity (1 + a’)ωr at the rotorplane (see Figure 6.2). Figure 6.2 Velocities at the rotor plane θ is the local pitch of the blade, in other words the local angle between the chord and the plane of rotation. The local pitch is a combination of the pitch angle, θ , and the twist of the blade, β, as θ = θ + β, where the pitch angle p p is the angle between the tip chord and the rotorplane and the twist is measured relative to the tip chord. φ is the angle between the plane of rotation and the relative velocity, V , and it is seen in Figure 6.2 that the local angle rel of attack is given by:  = φ – θ. (6.6) Further, it is seen that: (1 – a)V o tan φ = –––––––– – . (6.7) (1 + a’)ωr It is recalled from the section concerning 2-D aerodynamics that the lift, by definition, is perpendicular to the velocity seen by the aerofoil and the drag is parallel to the same velocity. In the case of a rotor this velocity is V due rel to arguments given in the section about the vortex system of a wind turbine.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 48 48 Aerodynamics of Wind Turbines Further, if the lift and drag coefficients C and C are known, the lift L and l d drag D force per length can be found from equations (2.1) and (2.2): 1 2 L = – ρV cC (6.8) rel l 2 and: 1 2 D = – ρV cC . (6.9) rel d 2 Since we are interested only in the force normal to and tangential to the rotorplane, the lift and drag are projected into these directions (see Figure 6.3): p = Lcos φ + Dsin φ (6.10) N and: p = Lsin φ – Dcos φ (6.11) T 1 2 The equations (6.10) and (6.11) are normalized with respect to –ρV c yielding: 2 rel R is the vector sum of the lift L and the drag D; p and p are the normal and tangential components of N T R respectively. Figure 6.3 The local loads on a blade3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 49 The Classical Blade Element Momentum Method49 C = C cos φ + C sin φ (6.12) n l d and: C = C sin φ – C cos φ (6.13) t l d where: p N C = –––––––– (6.14) n 1 2 – ρV c 2 rel and: p T C = ––––––– (6.15) t 1 2 – ρV c 2 rel From Figure 6.2 it is readily seen from the geometry that: V sin φ = V (1 – a) (6.16) rel o and: V cos φ = ωr(1 + a’) (6.17) rel Further, a solidity σ is defined as the fraction of the annular area in the control volume which is covered by blades: c(r)B σ(r) = –––– – (6.18) 2r B denotes the number of blades, c(r) is the local chord and r is the radial position of the control volume. Since p and p are forces per length, the normal force and the torque on N T the control volume of thickness dr are: dT = Bp dr (6.19) N and: dM = rBp dr. (6.20) T3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 50 50 Aerodynamics of Wind Turbines Using equation (6.14) for p and equation (6.16) for V , equation (6.19) N rel becomes: 2 2 V (1 – a) o 1 dT = – ρB–––––––– –– cC dr. (6.21) n 2 2 sin φ Similarly, if equation (6.15) is used for p and equations (6.16) and (6.17) are T used for V , equation (6.20) becomes: rel V(1 – a)ωr (1 + a’) 1 o _______________ dM = – ρB cC rdr. (6.22) t 2 sin φ cos φ If the two equations (6.21) and (6.4) for dT are equalized and the definition of the solidity equation (6.18) is applied, an expression for the axial induction factor a is obtained: 1 a = ––––––––––––––– . (6.23) 2 4 sin φ ––––––––– 1 C n If equations (6.22) and (6.5) are equalized, an equation for a’ is derived: 1 a’ = ––––––––––––––– . (6.24) 4 sin φ cos φ –––––––––– 1 C t Now all necessary equations for the BEM model have been derived and the algorithm can be summarized as the 8 steps below. Since the different control volumes are assumed to be independent, each strip can be treated separately and the solution at one radius can be computed before solving for another radius; in other words for each control volume the following algorithm is applied. Step (1) Initialize a and a’, typically a = a’ = 0. Step (2) Compute the flow angle φ using equation (6.7). Step (3) Compute the local angle of attack using equation (6.6). Step (4) Read off C (α) and C (α) from table. l d Step (5) Compute C and C from equations (6.12) and (6.13). n t Step (6) Calculate a and a’ from equations (6.23) and (6.24). Step (7) If a and a’ has changed more than a certain tolerance, go to step (2) or else finish. Step (8) Compute the local loads on the segment of the blades.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 51 The Classical Blade Element Momentum Method51 This is in principle the BEM method, but in order to get good results it is necessary to apply two corrections to the algorithm. The first is called Prandtl’s tip loss factor, which corrects the assumption of an infinite number of blades. The second correction is called the Glauert correction and is an empirical relation between the thrust coefficient C and the axial induction T factor a for a greater than approximately 0.4, where the relation derived from the one-dimensional momentum theory is no longer valid. Each of these corrections will be treated in separate sections. After applying the BEM algorithm to all control volumes, the tangential and normal load distribution is known and global parameters such as the mechanical power, thrust and root bending moments can be computed. One has to be careful, however, when integrating the tangential loads to give the shaft torque. The tangential force per length p is known for each segment at T,i radius r and a linear variation between r and r is assumed (see Figure 6.4). i i i+1 The load p between r and r is thus: T i i+1 pT = A r + B (6.25) i i where: p – p T,i + l T,i A = (6.26) –––––––– i r – r i + l i and: p r – p r T,i i + l T,i + l i B = (6.27) –––––––––––– i r – r i + l i The torque dM for an infinitesimal part of the blade of length dr is: 2 dM = rp dr = (A r + B r)dr (6.28) T i i and the contribution M to the total shaft torque from the linear tangential i,i+1 load variation between r and r is thus: i i+1 r 1 1 i+1 1 1 3 2 3 3 2 2 M = — A r + —B r = —A (r – r )+ —B (r – r ). (6.29) i,i 1 i i i i+1 i i i+1 i r i 3 2 3 2 The total shaft torque is the sum of all the contributions M along one blade i,i+1 multiplied by the number of blades: N–1 M = B M . (6.30) tot i,i+1 13212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 52 52 Aerodynamics of Wind Turbines Figure 6.4 Alinear variation of the load is assumed between two different radial positions r and r i i+1 Prandtl’s Tip Loss Factor As already mentioned, Prandtl’s tip loss factor corrects the assumption of an infinite number of blades. For a rotor with a finite number of blades the vortex system in the wake is different from that of a rotor with an infinite number of blades. Prandtl derived a correction factor F to equations (6.4) and (6.5): 2 dT = 4rV a(1 – a)Fdr (6.31) o and: 3 dM = 4r Vω(1 – a)a’Fdr. (6.32) o F is computed as: 2 –1 –f F = — cos (e ), (6.33) 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 53 The Classical Blade Element Momentum Method53 where: B R – r f = — – (6.34) 2 rsinφ B is the number of blades, R is the total radius of the rotor, r is the local radius and φ is the flow angle. Using equations (6.31) and (6.32) instead of equations (6.4) and (6.5) in deriving the equations for a and a’ yields: 1 a = ––––––––– (6.35) 2 4F sin φ ––––––– +1 C n and: 1 a’ = –––––––––––– (6.36) 4F sin φ cos φ ––––––––––– –1 C t Equations (6.35) and (6.36) should be used instead of equations (6.23) and (6.24) in step 6 of the BEM algorithm and an extra step computing Prandtl’s tip loss factor F should be put in after step 2. Deriving Prandtl’s tip loss factor is very complicated and is not shown here, but a complete description can be found in Glauert (1935). Glauert Correction for High Values of a When the axial induction factor becomes larger than approximately 0.4, the simple momentum theory breaks down (see Figure 4.5, where the different states of the rotor are also shown). Different empirical relations between the thrust coefficient C and a can be made to fit with measurements, for example: T 1 4a(1 –a)Fa ≤ — 3 (6.37) C = 1 T  1 — 4a(1 – (5–3a)a)Fa — 4 3 or: 4a(1 –a)Fa ≤ a c (6.38) C = T 2  4(a +(1–2a )a)Fa a c c c The last expression is found in Spera (1994) and a is approximately 0.2. F c is Prandtl’s tip loss factor and corrects the assumption of an infinite number3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 54 54 Aerodynamics of Wind Turbines of blades. In Figure 6.5 the two expressions for C (a) are plotted for F = 1 T and compared to the simple momentum theory. Figure 6.5 Different expressions for the thrust coefficient C versus the axial induction T factor a From the local aerodynamics the thrust dT on an annular element is given by equation (6.21). For an annular control volume, C is by definition: T dT __________ . (6.39) C = T 1 2 –ρV 2rdr 2 o If equation (6.21) is used for dT, C becomes: T 2 (1 – a) C n C = ––––––––– – – (6.40) T 2 sin φ This expression for C can now be equated with the empirical expression T (6.38).3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 55 The Classical Blade Element Momentum Method55 If a a : c 2 (1 – a) σC n 4a(1 – a)F = –––––––––––– (6.41) 2 sin φ and this gives: 1 a = ––––––––– (6.42) 2 4F sin φ ––––––– +1 C n which is the normal equation (6.35). If a a : c 2 (1 – a) σC n 2 4(a + (1 – 2a )a)F = ––––––––– (6.43) c c 2 sin φ and this gives: 1 2 2 a = –2 + K(1 – 2a ) – (K(1 – 2a ) +2) + 4(Ka – 1) (6.44) c c c 2 where: 2 4Fsin φ K = –––––––– (6.45) σ C n In order to compute the induced velocities correctly for small wind speeds, equations (6.44) and (6.42) must replace equation (6.35) from the simple momentum theory. Annual Energy Production The BEM method has now been derived and it is possible to compute a power curve, in other words the shaft power as a function of the wind speed V. In order to compute the annual energy production it is necessary to o combine this production curve with a probability density function for the wind. From this function the probability, f (V V V ), that the wind speed i o i+1 lies between V and V can be computed. Multiplying this with the total i i+1 number of hours per year gives the number of hours per year that the wind speed lies in the interval V V V . Multiplying this by the power (in kW) i o i+1 produced by the wind turbine when the wind speed is between V and V i i+13212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 56 56 Aerodynamics of Wind Turbines gives the contribution of the total production (in kWh) for this interval. The wind speed is discretized into N discrete values (V, i = 1,N ), typically with i 1m/s difference (see Figure 6.6). Figure 6.6 Probability f(V V V ) that the wind speed lies between V and V and a i o i+1 i i+1 power curve in order to compute the annual energy production for a specific turbine on a specific site It must be noted that the production must be corrected for losses in the generator and gearbox, which have a combined efficiency of approximately 0.9. Typically the probability density function of the wind is given by either a Rayleigh or a Weibull distribution. The Rayleigh distribution is given by the mean velocity only as:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 57 The Classical Blade Element Momentum Method57 2  V  V o o h (V ) = — — exp – — — (6.46) R o – – 2 2 V 4 V In the more general Weibull distribution, some corrections for the local siting (for example for landscape, vegetation, and nearby houses and other obsta- cles) can be modelled through a scaling factor A and a form factor k: k–l k kV V o o h (V ) = — — exp – — (6.47) w o AA A The parameters k and A must be determined from local meteorological data, nearby obstacles and landscape. Help in doing this can be obtained from the European Wind Atlas (Troen and Petersen,1989). From the Weibull distri- bution, the probability f (V V V ) that the wind speed lies between V and i o i+1 i V is calculated as: i+1 k k V V i i+1 o f (V V V ) = exp – — – exp – — (6.48) i i +1 AA The total annual energy production can thus be evaluated as: N–1 1 AEP = – (P(V ) + P(V)) · f (VVV )·8760. (6.49) i +1 i i o i +1 2 i=1 Example Having derived the BEM method and shown how annual energy production can be calculated, it is time for a simple example to illustrate the accuracy that can be obtained for a real turbine. The following example is of a Nordtank NTK 500/41 wind turbine. The turbine is stall regulated (fixed pitch) and the main parameters are listed below: Rotational speed: 27.1 rpm; 3 Air density: 1.225kg/m ; Rotor radius: 20.5m; Number of blades: 3; Hub height: 35.0m; and Cut-in wind speed: 4m/s; cut-out wind speed: 25m/s.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 58 58 Aerodynamics of Wind Turbines Blade description r m twist degrees chord m 4.5 20.0 1.63 5.5 16.3 1.597 6.5 13.0 1.540 7.5 10.05 1.481 8.5 7.45 1.420 9.5 5.85 1.356 10.5 4.85 1.294 11.5 4.00 1.229 12.5 3.15 1.163 13.5 2.60 1.095 14.5 2.02 1.026 15.5 1.36 0.955 16.5 0.77 0.881 17.5 0.33 0.806 18.5 0.14 0.705 19.5 0.05 0.545 20.3 0.02 0.265 Since the power depends directly on the air density ρ, according to the 3 standards the computations must be performed for ρ = 1.225kg/m . The difficult part is to find reliable aerofoil data C () and C () for the different l d aerofoils applied along the span. The data available in the literature are for thin aerofoils not much thicker than 20 per cent of the chord and for angles of attack only slightly above C . For structural reasons it is desirable to use l,max very thick aerofoils of approximately 40 per cent of the chord at the root of the blades in order to absorb the high bending moments. Further, the boundary layers on the rotating blades are influenced by centrifugal and Coriolis forces, which alter the post-stall lift and drag coefficients from those measured in a wind tunnel. It is therefore clear that it requires significant engineering skill and experience to construct good aerofoil data for thick aerofoils at high angles of attack including 3-D effects. In Snel et al (1993) and Chaviaropoulos and Hansen (2000) some guidelines are given to correct3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 59 The Classical Blade Element Momentum Method59 for the rotational effects. When the power curve from an actual wind turbine is known from measurements, it is common to calibrate the aerofoil data afterwards in order to achieve better agreement between measurements and computations. If a new blade is constructed which is not too different from blades for which the aerofoil data have been calibrated, it is possible to predict the power curve very well. But if a new blade is to be designed with completely new aerofoils, one has to be very careful in using the computed results. The actual geometry and aerofoil data on the blade of the Nordtank NTK 500/41 is not shown here, but an educated guess has been applied to extrapolate the data into high angles of attack. The power curve has been measured (Paulsen, 1995) and the comparison between the computed and measured power curve is shown in Figure 6.7 to give an idea of the accuracy of the BEM model. Figure 6.7 Comparison between computed and measured power curve, i.e. mechanical shaft power as a function of the wind speed It is seen in Figure 6.7 that except for at very high wind speeds the BEM method captures the measurements very well. The power curve is often shown in non-dimensional form as in Figure 6.8 or 6.9. Figure 6.8 shows that this particular wind turbine has a maximum efficiency of approximately C = p 0.5 for a tip speed ratio λ between 9 and 10. The advantage of plotting the -1 power coefficient as a function of the inverse tip speed ratio is that λ increases linearly with the wind speed V. o3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 60 60 Aerodynamics of Wind Turbines Figure 6.8 Power coefficient C as a function of the tip speed ratio  = ωR/V p o –1 Figure 6.9 Power coefficient C as a function of the inverse tip speed ratio  =V /ωR p o3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 61 The Classical Blade Element Momentum Method61 If the turbine were erected at a site where the Weibull parameters are k = 1.9 and A = 6.8m/s, the annual energy output from the mechanical power would 6 be 1.09·10 kWh, which corresponds to the consumption of approximately 250 households. The actual number it could supply would be less, however, since the losses from the generator and gearbox have not been taken into account. It should also be remembered that the Weibull parameters vary from site to site and have to be evaluated for each individual siting. The example shown here is for a stall regulated wind turbine, but the BEM method can also be used to predict the necessary pitch setting of a pitch regulated wind turbine. When the pitch is referred to for an entire wind turbine, it means the angle between the chord line of the tip aerofoil and the rotor plane. A pitch regulated wind turbine may operate at a fixed pitch until a certain nominal power is generated. For higher wind speeds the blades are pitched normally with the leading edge into the wind in order to keep this nominal power. Therefore the power curve of a pitch regulated wind turbine is absolutely flat after the nominal power has been reached. More runs are required with the BEM method to predict the power curve and the pitch setting for different wind speeds. One procedure for computing a pitch regulated wind turbine is sketched in Figure 6.10. Figure 6.10 Sketch of a power curve for a pitch controlled wind turbine3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 62 62 Aerodynamics of Wind Turbines At point A the nominal power is reached and it is necessary to change the pitch. BEM calculations are made for a pitch of θ and θ . It is seen that for 2 3 a wind speed of V(B) the pitch much be changed to θ in order to obtain the o 2 nominal power. The BEM method as derived in this chapter is steady, so it is not possible to compute the transient from point A to point B and from point B to point C. This requires an extended BEM method as shown in a later chapter. References Chaviaropoulos, P. K. and Hansen, M. O. L. (2000) ‘Investigating three-dimensional and rotational effects on wind turbine blades by means of a quasi-3D Navier-Stokes solver’, Journal of Fluids Engineering, vol 122, pp330–336 Glauert, H. (1935) ‘Airplane propellers’, in W. F. Durand (ed) Aerodynamic Theory, vol 4, Division L, Julius Springer, Berlin, pp169–360 Paulsen, U. S. (1995) ‘Konceptundersøgelse Nordtank NTK 500/41 Strukturelle Laster’ (in Danish), Risø-I-936(DA), Risø National Laboratory Snel, H., Houwink, B., Bosschers, J., Piers, W. J., van Bussel, G. J. W. and Bruining, A. (1993) ‘Sectional prediction of 3-D effects for stalled flow on rotating blades and comparison with measurements’, in Proceedings of European Community Wind Energy Conference 1993, pp395–399 Spera, D. A. (1994) Wind Turbine Technology, ASME Press, New York Troen, I. and Petersen, E. L. (1989) European Wind Atlas, Risø National Laboratory3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 63 7 Control/Regulation and Safety Systems The control or regulation system ensures that the turbine operates within the design range, in other words that it: • keeps the rotational speed within a certain range; • yaws the turbine; • keeps the power output within a certain range; or • starts and stops the turbine. Further, the control system can ensure a smooth power output P(t) and/or may optimize the power output at lower wind speeds. To limit the power at high wind speeds, the following three strategies may be used, the first two being by far the most common: 1 stall regulation; 2pitch regulation; and 3 yaw control. Stall Regulation Stall regulation is mechanically the most simple, since the blades are fixed to the hub and cannot be pitched. A stall regulated wind turbine is normally operated at an almost constant rotational speed and thus the angle of attack increases as the wind speed increases (see Figure 6.2). Eventually, as the local angles of attack are increased, the blades stall, causing the lift coefficient to decrease and the drag coefficient to increase, yielding a lower tangential load according to equation (6.13). The power decrease depends on the pitch angle (angle between rotor plane and tip chord), the twist and chord distributions, and the aerofoils used for the blades. If a site test shows that the3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 64 64 Aerodynamics of Wind Turbines power is not limited sufficiently, it is necessary to unbolt the blades and change their fixed pitch setting. On a stall regulated wind turbine an asynchronous generator is often used whereby the rotational speed is almost constant and determined by the torque characteristic of the generator, in other words the shaft torque into the generator, M , as a function of the rotational G speed of this shaft, n. A typical torque characteristic is illustrated in Figure 7.1, where it is seen that the asynchronous generator can act both as a motor and as a generator. The motor mode can be used to start the turbine. The sign here is deemed positive when the generator is producing electricity. The rotational speed of the generator will be between n and n and the torque o nom will equal the torque produced by the rotor blades M at the generator shaft. R The rotational speed of the generator for zero shaft torque, n , for an asyn- o chronous generator is: n = 60f /p, (7.1) o grid where f is the frequency of the grid (in Europe 50Hz) and p denotes the grid number of pole pairs. n is thus 1500 rpm for 4 poles and 1000 rpm for 6 o poles. The rotational speed of the generator is higher than the rotational speed of the rotor and therefore there is a gearbox between the generator and the rotor. The relationship between the rotational speed of the rotor, ω, and the rotational speed of the generator, n, is given through the transmission factor r as ω = n/r. The relative difference between the actual rotational speed n and n is called the slip SL = (n – n )/n and for a normal stall regulated wind o o o turbine the slip is about 1–3 per cent. This means that the rotational speed of the rotor is almost constant and the possibility of using the rotor as a flywheel to store energy – from a gust, for example – is small. Changes in the rotor torque M from, for example, turbulence in the wind are thus almost imme- R diately transferred to the generator torque M and thus to the produced G electrical power P = M 2n/60. Consider a wind turbine operating at point EL G A on Figure 7.1 and the wind speed increasing. In this case the torque, M , R from the rotor blades will also increase and the rotor accelerate according to equation (7.2): dω I ––– = M – M , (7.2) R G dt until M again equals M at point B. I denotes the moment of inertia of the G R rotor about the rotational axis and here M and M are the rotor and generator R G torque at the rotor shaft.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 65 Control/Regulation and Safety Systems65 Figure 7.1 Atypical torque characteristic for an asynchronous generator If M exceeds the maximum point on the torque characteristic, M , or R Break down the generator is disconnected from the grid, the term M – M on the R G right–hand side of equation (7.2) is always positive and the rotor will start to accelerate; in this case the rotational speed can become so high that a risk of breakdown exists. The safety system must detect this and ensure that the rotor is stopped. On a stall regulated wind turbine it is common that the outer parts of the blades are activated by centrifugal force to turn 90° and thus act as an aerodynamic brake limiting M , (see Figure 7.2). R An example of a time history of a start-up at a high wind speed for a stall regulated wind turbine is shown in Figure 7.3. This figure shows actual measurements performed on the Elkraft 1MW demonstration wind turbine, sited near Copenhagen. This turbine can run as a stall regulated machine as well as a pitch regulated machine (see later). The first curve in Figure 7.3 shows the wind speed at hub height, the second curve shows the rotational speed of the generator n and the last curve shows the corresponding power as a function of time. When starting the turbine, at t = 420s, the generator is switched off; in other words M is zero and the rotor starts to accelerate until G3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 66 66 Aerodynamics of Wind Turbines n = n at t = 445s. Then the generator is connected instantaneously, giving o rise to an overshoot in the power, as seen directly in the power time history. Figure 7.2 Turnable tip used as an aerodynamic brake and activated by the centrifugal force Pitch Regulation (Constant Rotational Speed) On a pitch regulated machine it is possible to actively pitch the entire blade and thus to change simultaneously the angles of attack along its entire length. One way of controlling the pitch is sketched in Figure 7.4, where a piston placed within the main shaft changes the pitch. The position of the piston is determined by applied hydraulic pressure. If the hydraulic pressure is lost a3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 67 Control/Regulation and Safety Systems67 Figure 7.3 Starting a stall regulated wind turbine at high wind speed: Elkraft 1MW demonstration wind turbine at Avedøre, Denmark spring will retract the piston, twisting the leading edges of the blades up against the wind. It is worth mentioning here that this is not the only way to pitch the blades: each blade could be fitted with a small electrical motor, meaning that they could be pitched independently. A pitched blade can act as an aerodynamic brake and it is no longer necessary to include tip brakes as on a stall regulated machine. By pitching the entire blade it is possible to control the angles of attack and thus the power output. Normally, the power is reduced by decreasing the angles of attack by pitching the leading edge of the blades up against the wind, in other words by increasing θ in the expression for the angle of attack = φ –θ, where φ is the flow angle. Alternatively, one could reduce the power output by increasing the angle of attack, thus forcing the blades to stall. This is called active stall. Due to the turbulent nature of the wind, the instantaneous power output of a pitch regulated machine will often exceed the rated power, the timescales of these fluctuations being smaller than the time it takes to pitch the blades (see Figure 7.5). Figure 7.5 shows a start-up of the same3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 68 68 Aerodynamics of Wind Turbines Figure 7.4 Sketch of mechanism to change the pitch of the blades through a piston placed inside the main shaft turbine as in Figure 7.3 but now running as a pitch controlled machine. It is seen that the start-up occurs much more smoothly since the blades are gradually pitched from 50° to about 15° (second curve in Figure 7.5). Comparing the time history of the power output after the start-up for the turbine running as a pitch regulated wind turbine (Figure 7.5) with the time history of the power output for the same turbine running as a stall regulated wind turbine (Figure 7.3), it is seen that the peaks are smaller when the turbine is running as a stall regulated turbine. Since turbulent fluctuations occur much faster than the time it takes to pitch the blades, the power from a pitch regulated wind turbine follows for some time the stationary power curve for a fixed pitch turbine (see Figure 6.10). For high wind speeds the3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 69 Control/Regulation and Safety Systems69 stationary power curve for a pitch regulated machine at a fixed pitch has a much higher slope dP/dV than the corresponding stationary power curve for o a stall regulated wind turbine and thus a larger variation ∆ P = dP/dV ·∆ V for o o the wind speed interval ∆ V . This is the reason why the power fluctuations o are lower for a stall regulated wind turbine than for a pitch regulated machine at high wind speeds. The control diagram when the pitch is used to control the power is shown in Figure 7.6. The controller of a classical pitch regulated wind turbine is not, however, responding to the wind speed but reacts directly to the power as: dθ KI(P– P ) p ref ––– = –––––––––, (7.3) dt θ (t) p ––––– 1+ KK where KI is an integration constant and KK is the gain reduction that reduces the pitch rate at high values of the pitch angle itself. High values of the pitch angle correspond to high wind speeds, where the loads are also very big and thus very sensitive to the pitch angle. Figure 7.5 Starting a pitch regulated wind turbine at high wind speed: Elkraft 1MW demonstration wind turbine at Avedøre, Denmark3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 70 70 Aerodynamics of Wind Turbines Figure 7.6 Control diagram for a pitch regulated wind turbine controlling the power To overcome the problem of the large peaks in power and loads on a pitch ® regulated wind turbine at high wind speeds, a system called OptiSlip is used ® by VESTAS. OptiSlip utilizes the fact that the torque characteristic; in other words the slip, for an asynchronous generator can be altered by changing an inner resistance R in the generator. A torque characteristic with a constant torque, after a certain value, M = M = const, can thus be obtained as G nom indicated by the thick curve in Figure 7.7, which shows the effective characteristic obtained with different resistances, R , R , R ,… . In high winds 1 2 3 ® the power output from a wind turbine using OptiSlip will be almost constant, with very small fluctuations around the nominal value. When the rotor torque M exceeds M the rotor will start to accelerate in R nom accordance with equation (7.2). The control system detects this and starts pitching the blades. The advantage of this is that the time needed to accelerate the rotor is longer than that needed to pitch the blades, so there is enough time available to physically move the blades, and at the same time a much smoother power output is obtained since the torque on the generator shaft is almost constant. Further, the loads are also reduced, which increases the fatigue lifetime. The pitch system is thus controlling the rotational speed and not the power.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 71 Control/Regulation and Safety Systems71 Figure 7.7 Torque characteristic for an asynchronous generator with variable slip There follows an example illustrating stall regulation, pitch regulation and active stall, based on the previously introduced stall controlled NTK 500/41 wind turbine. Using the same blades, the same rotational speed and the same Weibull parameters, the potential increase from a pitch mechanism in the annually captured energy is investigated. It is decided to use the same generator; in other words the nominal mechanical power remains unaltered at approximately 580kW. For each wind speed an optimum pitch angle is sought by varying the pitch in a BEM calculation (see Figure 7.8). It is seen that the mechanical power for a constant wind speed has an optimum value for the global pitch, but that for this turbine the variation is small. For higher wind speeds the optimum value of the power exceeds the nominal power, as shown in Figure 7.9. In this figure it is also seen that there exist two choices of pitch yielding the nominal power. The smaller value, θ p = 0.16°, corresponds to active stall, since the local angles of attack are higher than the limit for unseparated flow. The higher value, θ = 20.2°, corresponds p to classical pitch regulation, whereby the local angles of attack and thus the loads on the blades are reduced.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 72 72 Aerodynamics of Wind Turbines Figure 7.8 Optimum pitch angle for different wind speeds Figure 7.9 Variation of the mechanical power with the pitch for a wind speed of 20m/s on an NTK 500/41 wind turbine3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 73 Control/Regulation and Safety Systems73 The computed mechanical power curves for the NTK 500/41 wind turbine running as a stall and pitch controlled wind turbine are plotted in Figure 7.10. For the variable pitch machine, the optimum values for the pitch have been used for the lower wind speeds. These optimum pitch values are shown in Figure 7.11 as a function of the wind speed. For wind speeds slightly higher than 14m/s, the blades must be pitched to ensure a power below the nominal value. The lower branch in Figure 7.11 shows the pitch setting on a machine controlled by active stall and the upper branch shows the pitch setting on a pitch controlled wind turbine. Figure 7.10 Computed power curve for the NTK 500/41 wind turbine running as a stall controlled or pitch regulated machine Assuming the same Weibull parameters as in the example concerning the NTK 500/41 machine, i.e. k = 1.9 and A = 6.8m/s, the annual energy 6 production from the turbine running as a pitch regulated machine is 1.11·10 kWh. In this example an increase in the annual energy production of approximately 2 per cent has thus been achieved by changing from a stall to a pitch regulated machine. (It should be noted that the annual energy production in both cases has been computed using mechanical power; in other words the losses in the generator and in the gearbox have not been taken into account.) The main contribution to this increase comes from the3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 74 74 Aerodynamics of Wind Turbines Figure 7.11 Optimum values for the pitch and the necessary pitch to avoid exceeding the nominal power shape of the power curve just before reaching the nominal power, where the power curve of the pitch regulated wind turbine is steeper than that of the stall regulated machine. Due to the turbulent characteristics of the wind, however, the power curve on a real pitch regulated machine is not so steep close to P as the theoretical curve indicates in Figure 7.10, so the increase nom in the annual production is probably slightly lower than estimated from the theoretical power curves. Yaw Control Instead of limiting the power output using pitch or stall regulation, it is possible to control the yaw of the turbine. On normal pitch or stall regulated machines, it is common to have a yaw drive, which is constantly trying to rotate the nacelle to minimize the yaw misalignment in order to get as much air through the rotor disc as possible. On a yaw controlled machine, the rotor is turned out of the wind in high winds to limit the airflow through the rotor and thus the power extraction. Yaw control was used on the old Western mills. For larger machines, only the 1.5MW Italian prototype called GAMMA 60 is, to the author’s knowledge, yaw controlled. Therefore yaw control will not be treated further in this text.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 75 Control/Regulation and Safety Systems75 Variable Speed The power coefficient is a function of the tip speed ratio and the pitch angle, C (θ ,λ), and by applying variable speed on a pitch regulated rotor it is p p possible to run the turbine at the optimum point C occurring at θ and p, max p,opt λ . From Figure 6.8, which shows the C -λ curve for the NTK 500/41 opt p machine, it is seen that the turbine for this pitch angle runs most efficiently at λ approximately equal to 9.8. From a C -λ curve, a plot showing the power p P as a function of the rotational speed ω for different wind speeds can be derived as: 1 3 P = —ρV AC (λ) (7.4) o p 2 and: 1 λV o ω = ––– . (7.5) R This has been done for the NTK 500/41 machine and the result is shown in Figure 7.12. The turbine is equipped with an asynchronous generator forcing the blades to rotate at 27.1 rpm, indicated by the vertical line. It is seen that the turbine is running most efficiently at a wind speed V of approximately o 7m/s. Some stall and pitch regulated wind turbines using asynchronous generators, in other words running at a fixed rotational speed, therefore have two generators, one which is efficient at lower wind speeds and one which is efficient at higher wind speeds. If another type of generator had been used, one which is able to run at different rotational speeds, the turbine could be operated at the optimum rotational speed for each wind speed, as indicated in Figure 7.12 by the operational line that intersects all the top points in the curves for the different wind speeds. All points on this line correspond to the highest C that can be obtained for the applied pitch angle. It is noted that at p a fixed value of the tip speed ratio λ, the angular velocity ω and thus the tip speed will increase proportionally with the wind speed according to equation (7.5). Due to noise emission, the tip speed is limited to approximately ωR = 70m/s and therefore the optimum tip speed ratio can only be obtained for lower wind speeds. The alternating current produced by such a generator will have a frequency different from the frequency of the grid (50Hz in Europe). Therefore the alternating current (AC) is first transformed into direct current (DC) and then back to alternating current with the correct frequency using a so-called ACDC/DCAC device.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 76 76 Aerodynamics of Wind Turbines Figure 7.12 Constant speed versus variable speed Combining the fact that power is the torque multiplied by the angular velocity, P= Mω, with equations (7.4) and (7.5) yields: 1 2 3 3 M = – ρω R AC (λ,θ )/λ . (7.6) p p 2 The generator characteristic for the highest obtainable power coefficient C (θ , λ ) is thus: p,max p,opt opt 1 2 3 3 2 M = – ρω R AC (θ ,λ )/λ = const · ω . (7.7) opt p,max p,opt opt opt 2 This characteristic should be used until the maximum allowable tip speed is reached. Thereafter the torque should be constant, such as the characteristic shown in Figure 7.7. If the rotor torque exceeds this maximum generator torque it will accelerate according to equation (7.2), but since the inertia of the rotor is large it will take some time for the rotor to build up speed. The electrical output will be almost constant since the generator torque is constant, and the pitch system should control the rotational speed and not the power. Since the change in rotational speed is slow due to the high inertia of the rotor, there is plenty of time for the pitch mechanism to act. This will overcome the problem of the large spikes/fluctuations from a pitch regulated3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 77 Control/Regulation and Safety Systems77 wind turbine operating at constant rotational speed as seen in the power output shown in Figure 7.5. The control diagram for a pitch regulated variable speed machine is shown in Figure 7.13. Figure 7.13 Schematic drawing of a control diagram for a pitch regulated variable speed machine One way of obtaining variable speed with an asynchronous generator is to apply a so-called double feed induction generator (DFIG). Here a control device is connected to the generator that effectively controls the net frequency seen by the rotor in the generator. With a DFIG it is thus possible to create a variable speed that can reduce fatigue loads and improve power quality and that for lower wind speeds can produce slightly more power by running at the optimum C . However, the extra cost of the control system and p the necessary converters also has to be considered.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 78 8 Optimization Having derived all the necessary equations to compute a given wind turbine, one should be in a position to use these to compute an optimum design. First, an optimum design must be defined, but since the purpose of a wind turbine is to produce electricity, and this should be at a competitive cost, the object function is a design that can last for a typical design lifetime of 20 years and at a given site minimize energy production cost (in /kWh). To do this it is necessary to estimate a cost function for every component of the wind turbine and a maintenance cost. From a purely technical point of view, the optimum design could be a wind turbine which for a given rotor diameter captures as many kWh/year as possible. If the turbine is sited where the wind speed V is constant in time, it is obvious to optimize o the power coefficient at this wind speed. Since in practice wind speed is not constant, however, an optimum design could have lower C as p,max sketched in the second design in Figure 8.1. As shown in the previous chapter, the annual energy production is a combination of the wind distribution and the power curve. Thus the optimum design also depends on the actual siting. A BEM code as described earlier can be coupled with an optimization algorithm with appropriate constraints to optimize the geometry of the blades, for example. Of course, it is imperative afterwards to verify that the calculated optimum design will also survive the entire design period, taking both extreme and fatigue loads into account. It is possible to compute analytically the geometry of design 1 in Figure 8.1 with the already derived equations as described in, for example, Glauert (1935). Such a design is more interesting for an aeroplane propeller, which is mainly operating at cruise speed, but it could perhaps also be interesting for a wind turbine with variable rotational speed. Ideally such a machine can be kept by a control mechanism at the optimum tip speed ratio, λ , and pitch angle, θ , as described in Chapter 7. opt p,opt First a good aerofoil is chosen; this must be relatively roughness insensitive and possess an acceptable stall characteristic. Noise considerations might3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 79 Optimization79 Figure 8.1 Two different designs: design 1 has a high C but C drops off quickly at p,max p different tip speed ratios; design 2 has a lower C but performs better over a range of tip p,max speed ratios also influence the choice of aerofoil. A possible choice is the NACA63-415 aerofoil, which has proven to possess good stall characteristics on stall regulated wind turbines. Since it is designed to operate at one point only, it is ensured that the effective angle of attack has an optimum value along the entire span. The optimum value is where the ratio between the lift and the drag is highest. 6 In Figure 8.2 it is seen that the NACA63-415 aerofoil at Re = 3⋅10 has a maximum value of C /C of approximately 120 at an angle of attack  of 4°. l d It is also seen that this maximum value drops to approximately 67 when standard roughness is added to the aerofoil at the same optimum angle of attack of 4°. In the following example the values with roughness are used: C (4°) = 0.8 and C (4°) = 0.012. Further, the number of blades B is chosen l,opt d,opt as 3 and the design point as λ = ωR/V = 6. Since an angle of attack of 4° is o chosen in the design point, the flow is attached to the blades and equations (4.32) and (4.38) are valid and can be combined to give an optimum3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 80 80 Aerodynamics of Wind Turbines Source: Abbot and von Doenhoff (1959), reproduced with permission. Figure 8.2 Airfoil data for the NACA63-415 airfoil3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 81 Optimization81 Source: Abbot and von Doenhoff (1959), reproduced with permission. Figure 8.2 continued3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 82 82 Aerodynamics of Wind Turbines relationship between x and a: 3 2 2 2 16a – 24a + a(9 – 3x ) – 1 + x = 0 (8.1) The optimum value of a’ is found using (4.38) and the optimum local pitch angle can then be computed as: θ = φ –  , (8.2) opt opt since it is recalled that the flow angle is computed as: (1 – a)V (1 – a) o tan φ = –––––––– = ––––––– . (8.3) (1 a’)ωr (1 a’)x The optimum chord distribution is found from equation (6.23) using the optimum values for a and a’: Figure 8.3 Optimum pitch distribution (neglecting Prandtl’s tip loss factor) for λ = 6,  = 4, C = 0.8, C = 0.012 and B = 3 l,opt d,opt opt3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 83 Optimization83 2 c(x)8ax sin φ ––– = ––––––––– , (8.4) R (1–a)BC λ n where: C = C cos(φ) + C sin(φ) (8.5) n l,opt d,opt For λ = 6, α = 4, C = 0.8, C = 0.012 and number of blades B = 3, the opt l,opt d,opt optimum chord and pitch distribution can now be computed; the solution is shown graphically in Figures 8.3 and 8.4. Note, however, that in these simple expressions Prandtl’s tip loss factor has not been taken into account. Taking Prandtl’s tip loss factor into account will change the optimum pitch and chord distribution close to the tip. Figure 8.4 Optimum chord distribution (neglecting Prandtl’s tip loss factor) for λ = 6,  = 4, C = 0.8, C = 0.012 and B = 3 opt l,opt d,opt 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 84 84 Aerodynamics of Wind Turbines References Abbot, H. and von Doenhoff, A. E. (1959) Theory of Wing Sections, Dover Publications, New York Glauert, H. (1935) ‘Airplane propellers’, in W. F. Durand (ed) Aerodynamic Theory, vol 4, Division L, Julius Springer, Berlin, pp169–3603212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 85 9 Unsteady BEM Model To estimate the annual energy production from a wind turbine at a given site with known wind distribution it is sufficient to apply a steady BEM method as described in Chapter 6 to compute the steady power curve. But due to the unsteadiness of the wind seen by the rotor caused by atmospheric turbulence, wind shear and the presence of the tower it is necessary to use an unsteady BEM method to compute realistically the aeroelastic behaviour of the wind turbine. To do this a complete structural model of the wind turbine is also required; this must be coupled with the unsteady BEM method since, among other things, the velocity of the vibrating blades and the tower change the apparent wind seen by the blades and thus also the aerodynamic loads. Since the wind changes in time and space it is important at any time to know the position relative to a fixed coordinate system of any section along a blade. The fixed or inertial coordinate system can be placed at the bottom of the tower. Depending on the complexity of the structural model a number of additional coordinate systems can be placed in the wind turbine. The following example illustrates a very simple model, with the wind turbine described by four coordinate systems as shown in Figure 9.1. First, an inertial system (coordinate system 1) is placed at the base of the tower. System 2 is non-rotating and placed in the nacelle, system 3 is fixed to the rotating shaft and system 4 is aligned with one of the blades. Note that due to the orientation of coordinate system 2, the tilt angle θ must be tilt negative if the shaft is to be nose up as sketched in Figure 9.1. A vector in one coordinate system X = (x , y , z ) can be expressed in A A A A another coordinate system X = (x , y , z ) through a transformation matrix a : B B B B AB X = a X . (9.1) B AB A The columns in the transformation matrix a express the unit vectors of AB system A in system B. Further, the transformation from system B to system T A can be found as a = a . The rules will now be applied to the coordinate BA AB systems shown in Figure 9.1. First the transformation matrix a will be 12 constructed. To begin with, system 1 and 2 are identical with the exception3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 86 86 Aerodynamics of Wind Turbines Figure 9.1 Wind turbine described by four coordinate systems of the position of the origins. System 2 is first rotated about the x-axis with the angle θ . This gives a transformation matrix: yaw 10 0 a = 0 cosθ sinθ . (9.2) 1 yaw yaw   0 –sinθ cosθ yaw yaw Hereafter system 2 is rotated θ about the y-axis, yielding a transformation tilt matrix:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 87 Unsteady BEM Model87 cosθ 0 –sinθ tilt tilt a = 0 1 0 . (9.3) 2   sinθ 0 cosθ tilt tilt Since system 2 is not rotated about the z-axis, this transformation matrix then becomes: 100 a = 0 1 0 . (9.4) 3   001 The total transformation matrix, X = a X , between system 1 and system 2 2 12 1 is found as a = a a a . 12 3 2 1 · · Since the shaft in this simple model is assumed to be stiff, the only transform- ation between system 2 and system 3 is a rotation about the z-axis: cosθ sinθ 0 wing wing a = –sinθ cosθ 0 . (9.5) 23 wing wing   00 1 θ is the rotation of blade 1 as shown in Figure 9.2. wing System 4 is only rotated θ about the y-axis and the transformation matrix cone is thus: cosθ 0 –sinθ cone cone a = 0 1 0 . (9.6) 34   sinθ 0 cosθ cone cone Note that in order for the blades to cone as shown in Figure 9.1 the cone angle must be negative. A point along blade 1 is described in coordinate system 4 as r = (x,0,0), where x is the radial distance from the rotational axis 4 to the point on the blade. To transform this vector from system 4 to the inertial system 1, the following transformations are thus required:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 88 88 Aerodynamics of Wind Turbines Figure 9.2 Rotor seen from downstream First, the vector is transformed to system 3, which again can be transformed to system 2 and finally to system 1 as: r = a r 3 43 4 · r = a r 2 32 3 · r = a r 1 21 2 · or directly as: T T T r = a a a r = a a a r . 1 21 32 43 4 12 23 34 4 · · · · · · To find the coordinates of the point on the blade in system 1 the vector addition shown in Figure 9.3 is applied: xp r = yp = r + r + r . (9.7) t s b   zp3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 89 Unsteady BEM Model89 Figure 9.3 Apoint on the wing described by vectors All the vectors in equation (9.7), r , r and r , must be given in the same t s b coordinate system. For a position on a blade the natural coordinate system is 1 because the wind velocity is given in the fixed system. The undisturbed wind velocity seen by the blade is found by transforming this velocity V to system 4: 1 V x V = V = a · a · a V = a V. (9.8) y o 34 23 12 1 14 1   V z To find the relative velocity seen by the blade, V , the rotational velocity, V , rel rot plus the induced velocity, W, must be added as vectors to V in system 4 as o shown in Figure 9.4.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 90 90 Aerodynamics of Wind Turbines V = V + V + W ⇒ rel o rot V V –ωxcosθ W rel,y y cone y = + + (9.9)         V V 0 W rel,z z z Figure 9.4 Velocity triangle seen locally on a blade The angle of attack, , can be computed if the induced velocity, W, is known:  = φ – ( + θ ), (9.10) p where: V rel,z tan φ = ––––– – – . (9.11) – V rel,y Knowing , the lift and drag coefficients can be looked up in a table. The essence of the BEM method is to determine the induced velocity W and thus the local angles of attack. From a global consideration the rotor acts as a disc with a discontinuous pressure drop across it. The thrust generated by this pressure drop induces a velocity normal to the rotor plane, W, that deflects the wake as shown in n Figure 9.5. From simple momentum theory it is known that the induced velocity in the far wake is two times the induced velocity in the rotor plane. Bramwell (1976) states that Glauert’s relation between the thrust and this induced velocity for a gyrocopter in forward flight (similar to the lifting line result for an elliptically loaded circular wing) is:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 91 Unsteady BEM Model91 Figure 9.5 The wake behind a rotor disc seen from above T W = n·W = ––––––– , (9.12) n 2ρAV’ where V’ = V + n(n·W). o n is the unit vector in the direction of the thrust, which in system 3 has the coordinates n = (0,0,–1). Further, equation (9.12) reduces to the classical BEM theory for zero yaw misalignment. A gyrocopter in forward flight corresponds to 90° yaw misalignment and it is therefore postulated that equation (9.12) is valid for any yaw angle, which in fact is far from obvious. Figure 9.6 shows the local effect close to a blade section. 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 this radial position is assumed to affect the air in the area dA = 2rdr/B, so that all B blades cover the entire annulus of the rotor disc at radius r (see Figure 9.7). From equation (9.12) and Figure 9.6 the following can be derived for one blade: –L cos φ·dr –BLcos φ W = W = –––––––––––––––––––––– =––––––––––––––––– . (9.13) n z 2rdr 2ρ –––––– F V + f n(n·W) 4ρrFV +f n(n·W) o g o g B3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 92 92 Aerodynamics of Wind Turbines Figure 9.6 The local effect on an airfoil Figure 9.7 Annular strip3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 93 Unsteady BEM Model93 For the tangential component a similar expression is postulated: –BLsin φ W = W = ––––––––––––––––– (9.14) t y 4ρrFV + f n(n·W) o g F is Prandtl’s tip loss factor. It is noted that the equations for the induced velocity are identical with the classical BEM method in the case of zero yaw misalignment and can be reduced to the well-known expression: C = 4aF(1 – f · a), (9.15) T g where by definition for an annual element of infinitesimal thickness, dr, dT C = –––––––– . (9.16) T 1 2 –ρV dA 2 o f , usually referred to as the Glauert correction, is an empirical relationship g between the thrust coefficient C and the axial induction factor a in the T turbulent wake state and if equation (6.38) is used assumes the form: 1 for a ≤ a c f = (9.17) a a g c c ––(2 – ––) for a a .  c aa a is around 0.2. Further, it is noted that the equations must be solved c iteratively since the flow angle and thus also the angle of attack depend on the induced velocity itself. But the method described here is unsteady and thus time is used as relaxation; in other words after the blades have moved in one time step an azimuthal angle of ∆θ = ω∆ t (assuming a small ∆ t), values wing from the previous time step are used on the right-hand side of equations (9.13) and (9.14) for W when updating new values for the induced velocity. This can be done since the induced velocity changes relatively slowly in time due to the dynamic wake model. Dynamic Wake Model To take into account the time delay before equations (9.13) and (9.14) are in equilibrium with the aerodynamic loads, a dynamic inflow model must be applied. In two EU-sponsored projects (Snel and Schepers, 1995; Schepers and Snel, 1995) different engineering models were tested against3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 94 94 Aerodynamics of Wind Turbines measurements. One of these models, proposed by S. Øye, is a filter for the induced velocities consisting of two first order differential equations: dW dW qs int W + τ –––– – = W + k · τ ––––– – (9.18) int 1 qs 1 dt dt dW W + τ –––– – = W . (9.19) 2 int dt W qs is the quasi-static value found by equations (9.13) and (9.14), W an int intermediate value and W the final filtered value to be used as the induced velocity. The two time constants are calibrated using a simple vortex method as: 1.1 R τ = –––––––– · –– (9.20) 1 (1–1.3a) V o and: r 2 τ = (0.39 – 0.26 –– )·τ · (9.21) 2 1   R R is the rotor radius, k = 0.6 and a is the axial induction factor defined for zero yaw as a = W /V or more generally estimated as: n o V–V’ o a = ––––––– (9.22) V o Using equation (9.20), however, a is not allowed to exceed 0.5. Equations (9.18) and (9.19) can be solved using different numerical techniques. The one suggested here is to assume that the right-hand sides are constant, which allows them to be solved analytically, yielding the following algorithm: i 1 calculate W using equations (9.13) and (9.14); qs 2estimate right hand side of equation (9.18) using backward difference, i i–l W –W qs qs i H = W + k·τ qs l ∆ t i i–1 –∆ t 3 solve equation (9.18) analytically, W = H + (W – H)exp(–––); and int int τ l i i i–1 i –∆ t 4 solve equation (9.19) analytically, W =W int + (W – W int)exp(–––). τ 2 Applying a dynamic filter for the induced velocity is necessary in order to capture the time behaviour of the loads and power when the thrust is changed, by pitching the blades, for example. Figure 9.8 shows for a 2MW3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 95 Unsteady BEM Model95 machine the computed and measured response in the rotor shaft torque for a sudden change of the pitch angle. At t = 2s the pitch is increased from 0 to 3.7°, decreasing the local angles of attack. First the torque drops from 260kNm to 150kNm and not until approximately 10s later do the induced velocities and thus the power settle at a new equilibrium. At t = 32s the pitch is changed back to 0° and a similar overshoot in the torque is observed. The decay of the spikes seen in Figure 9.8 can only be computed with a dynamic inflow model; such a model is therefore of utmost importance for a pitch regulated wind turbine. Figure 9.8 Comparison between measured and computed time series of the rotor shaft torque for the Tjaereborg machine during a step input of the pitch for a wind speed of 8.7m/s Dynamic Stall The wind seen locally on a point on the blade changes constantly due to wind shear, yaw/tilt misalignment, tower passage and atmospheric turbulence. This has a direct impact on the angle of attack, which changes dynamically during the revolution. The effect of changing the blade’s angle of attack will not appear instantaneously on the loads but will take place with a time delay proportional to the chord divided by the relative velocity seen at the blade3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 96 96 Aerodynamics of Wind Turbines section. The response of the aerodynamic load depends on whether the boundary layer is attached or partly separated. In the case of attached flow the time delay can be estimated using the Theodorsen theory for unsteady lift and aerodynamic moment (Theodorsen, 1935). For trailing edge stall, in other words when separation starts at the trailing edge and gradually increases upstream at increasing angles of attack, so-called dynamic stall can be modelled through a separation function, f , as described in Øye (1991) (see s later). The Beddoes–Leishman model (Leishman and Beddoes, 1989) further takes into account attached flow, leading edge separation and compressibility effects, and also corrects the drag and moment coefficients. For wind tur- bines, trailing edge separation is assumed to be the most important phenomenon in terms of dynamic aerofoil data, but effects in the linear region may also be important (see Hansen et al, 2004). It is shown in Øye (1991) that if a dynamic stall model is not used, one might compute flapwise vibrations, especially for stall regulated wind turbines, which are non- existent on the real machine. For stability reasons it is thus highly recom- mended to at least include a dynamic stall model for the lift. For trailing edge stall the degree of stall is described through f as: s C = f C () + (1 – f )C (), (9.23) l s l,inv s l,fs where C denotes the lift coefficient for inviscid flow without any separat- l,inv ion and C is the lift coefficient for fully separated flow, for example on a l, fs flat plate with a sharp leading edge. C is normally an extrapolation of the l,inv st static aerofoil data in the linear region; one way of estimating C and f is l,fs s st shown in Hansen et al (2004). f is the value of f that reproduces the static s s aerofoil data when applied in equation (9.23). The assumption is that f will s always try to get back to the static value as: st df f – f s s s ––– = –––––– – , (9.24) τ dt which can be integrated analytically to give: st st ∆ t f (t + ∆ t) = f +( f (t) – f )exp(– –– ). (9.25) s s s s τ τ is a time constant approximately equal to A·c/V , where c denotes the local rel chord and V is the relative velocity seen by the blade section. A is a constant rel that typically takes a value of about 4. Applying a dynamic stall model, the aerofoil data is always chasing the static value at a given angle of attack that is also changing in time. If, for example, the angle of attack is suddenly3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 97 Unsteady BEM Model97 increased from below to above stall, the unsteady aerofoil data contains for a short time some of the inviscid/unstalled value, C , and an overshoot l,inv relative to the static data is seen. It can thus be seen as a model of the time needed for the viscous boundary layer to develop from one state to another. Figure 9.9 shows the result using the dynamic stall model for  = 15 + 2sin(12.57·t) with a time constant τ = 0.08s and the initial condition st f (0) = f . It is seen that the mean slope of the lift curve, dC /d, becomes s s l positive for the dynamic aerofoil data in stall, which is beneficial for stability. Figure 9.9 Example of the result using a dynamic stall model Yaw/Tilt Model If the rotor is yawed (and/or tilted), as shown in Figure 9.10, there will be an azimuthal variation of the induced velocity, so that the induced velocity is smaller when the blade is pointing upstream than when the same blade half a revolution later is pointing downstream. The physical explanation of this 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 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 is3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 98 98 Aerodynamics of Wind Turbines Figure 9.10 Yawed rotor disc in wind field not included, the BEM method will not be able to predict the restoring yaw moment. The following yaw model proposed by Glauert is also found in Snel and Schepers (1995): χ r W = W (1 + – tan(–)cos(θ –θ )), (9.26) o wing o R 2 where the wake skew angle, χ, is defined as the angle between the wind velocity in the wake and the rotational axis of the rotor (see Figure 9.10). W o is the mean induced velocity found by equations (9.13) and (9.14) and followed by (9.18) and (9.19). θ is the angle where the blade is deepest into o the wake. The skew angle can be found as: n·V’ cosχ = ––––– (9.27) nV’ ’ where n is the normal vector in the direction of the rotational axis (see Figure 9.10). 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, θ (t+∆ t) = θ (t)+ω·∆ t. The angle of attack can thus be evaluated wing wing from equation (9.10) and the lift, drag and moment coefficients can be looked up from a table. The normal, p , and tangential, p , loads can be determined z y from: p = Lcos φ + Dsin φ (9.28) z3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 99 Unsteady BEM Model99 and: p = Lsin φ – Dcos φ (9.29) y where: 1 2 L = –ρV cC (9.30) rel l 2 and: 1 2 D = –ρV cC (9.31) rel d 2 To summarize the unsteady BEM model: • read geometry and run parameters; • initialize the position and velocity of blades; • discretize the blades into N elements; • initialize the induced velocity – for n = 1 to max time step (t=n·∆ t); – for each blade; – for each element 1 to N; • compute the relative velocity to the blade element from equation (9.9) using old values for the induced velocity; • calculate the flow angle and thus the angle of attack from equations (9.10) and (9.11); • determine static C and C from a table; l d • determine dynamic aerofoil data using a dynamic stall model; • calculate lift using equation (9.30); • compute the loads p and p using equations (9.28) and (9.29); z y • compute new equilibrium values for the induced velocities W and W z y using equations (9.13) and (9.14); • find the unsteady induced velocities, W and W, using a dynamic wake z y model; and • in the case of yaw, calculate the azimuthal variation from equation (9.26) and compute the induced velocity for each blade.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 100 100 Aerodynamics of Wind Turbines In the case of an aeroelastic computation, where the structure is not considered stiff, the loads p and p are used to determine a local blade z y velocity V . This blade velocity must be taken into account when computing b the relative velocity and equation (9.9) should be extended to: V V –ωxcosθ W V rel,y y cone y b,y = + + – (9.32)           V V 0 W V rel,z z z b,z Deterministic Model for Wind The time averaged atmospheric boundary layer shown in Figure 9.11 can be modelled as: v x V (x) = V (H) –– , (9.33) o o  H where H is the hub height, x the distance from the surface and v a parameter giving the amount of shear. ν is in the range between 0.1 and 0.25. Figure 9.11 Deterministic wind velocity shear3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 101 Unsteady BEM Model 101 The wind is also influenced by the presence of the tower. A simple model for the influence of the tower is to assume potential flow (see Figure 9.12). Figure 9.12 The effect of the tower Coordinate system 1 is oriented so that the undisturbed wind velocity, V, is o aligned with the z-axis. In a polar coordinate system as shown in Figure 9.12, the radial and tangential velocity components around the tower, assuming potential flow, can be computed as: a 2 V = V (1–(—) )cosθ, (9.34) r o r and: a 2 V = –V (1+(—) )sinθ, (9.35) θ o r where a is the radius of the tower. Transforming the velocity to the cartesian coordinate system 1 yields: V = V cos θ – V sin θ z r θ (9.36) V = –V sin θ – V cos θ y r θ The following relation between r, θ and y, z can be used: z cosθ = – r y sinθ = – – . (9.37) r 2 2 r =  z + y3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 102 102 Aerodynamics of Wind Turbines It is noted that assuming potential flow is a bad approximation for a downwind machine, where each blade passes the tower wake once every revolution. Further, the turbulent part of the real atmospheric wind should be added for a realistic time simulation of a wind turbine. References Bramwell, A. R. S. (1976) Helicopter Dynamics, Edward Arnold Ltd, London Hansen, M. H., Gaunaa, M. and Madsen, H. A. (2004) ‘A Beddoes–Leishman type dynamic stall model in state-space and indicial formulations’, Risoe-R-1354(EN), Roskilde, Denmark Leishman, J. G. and Beddoes, T. S. (1989) ‘A semi-empirical model for dynamic stall’, Journal of the American Helicopter Society, vol 34, no 3, pp3–17 Øye, S. (1991) ‘Dynamic stall, simulated as a time lag of separation’, in K. F. McAnulty (ed) Proceedings of the 4th IEASymposium on the Aerodynamics of Wind Turbines , ETSU-N-118, Harwell Laboratory, Harwell, UK Schepers, J. G. and Snel, H. (1995) Dynamic Inflow: Yawed Conditions and Partial Span Pitch Control, ECN-C- -95-056, Petten, The Netherlands Snel, H. and Schepers, J. G. (1995) Joint Investigation of Dynamic Inflow Effects and Implementation of an Engineering Method, ECN-C- -94-107, Petten, The Netherlands Theodorsen, T. (1935) ‘General theory of aerodynamic instability and the mechanism of flutter’, NACA report no 496, National Advisory Committee for Aeronautics, pp413–4333212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 103 10 Introduction to Loads and Structures Having described in detail how to calculate the aerodynamic loads on a wind turbine, the following text concerns the structural issues that need to be addressed to ensure that the construction will not break down during its design lifetime of typically 20 years. Normally a breakdown is caused by an inadequate control system, extreme wind conditions, fatigue cracks or a defective safety system. A very dangerous breakdown may occur if power to the generator is lost. In this case there is no braking torque on the rotor, which, in the absence of a safety system such as mechanical or aerodynamic emergency brakes, is free to accelerate. Because the aerodynamic forces increase with the square of the rotor speed, the blades will bend more and more in the downwind direction and might end up hitting the tower or flying off due to centrifugal forces. It has been estimated (Sørensen, 1983) that torn- off blades from an overspeeding wind turbine could land up to about 300 m from the tower. Fortunately, violent failures are extremely rare and no humans have, to the author’s knowledge, ever been reported to have suffered injuries from debris flying off a wind turbine. Safety standards such as IEC 61400 (2004) exist to ensure that wind turbines operate safely. The standards define load cases, such as extreme gusts, which a wind turbine must be able to survive. Lightning is also known to have caused disintegration of blades. Fatigue is a very important issue in a wind turbine construction, which is 9 built to run for a minimum of 20 years and thus performs in the order of 10 revolutions. To estimate the loads on a wind turbine throughout its lifetime, the loads and hence the stresses in the material must either be computed using an aeroelastic code in a realistic wind field including turbulence or be measured directly on an already built turbine. Once the dynamic stresses are known, it is possible to calculate the fatigue damage using standard methods such as the Palmgren-Miner rule. 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 104 104 Aerodynamics of Wind Turbines Description of Main Loads In order that the following text can be better understood, a short description of the main loads on a horizontal-axis wind turbine is given. To extract energy from the wind it is necessary to slow down the wind speed using a force pointing in the upwind direction. This force is called the thrust and is caused by a jump in the pressure over the rotor, induced by the flow past the individual rotor blades. The total load has not only a component normal to the flow but also a tangential component in the rotational direction of the blades. The tangential load component delivers the shaft torque that turns the rotor. To characterize these loads it is common to state the flapwise and the edgewise bending moments at a position close to the root of the blades together with the tilt rotor moment and yaw rotor moment in the shaft between the rotor and the first bearing. The flapwise bending moment M flap (see Figure 10.1) stems from the normal forces (thrust), which tend to deflect the blades out of the rotor plane in the downwind direction: R M (r ) = ∫ rp (r)dr, (10.1) flap pos N r pos where p is the normal force per length, r the local radius, R the total radius N of the rotor and dr an incremental part of the blade. The edgewise bending moment is the bending moment in the rotor plane from the tangential force distribution. The edgewise bending moment is sometimes referred to as the lead-lag moment. Figure 10.1 Some main loads on a horizontal-axis wind turbine3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 105 Introduction to Loads and Structures 105 Figure 10.2 Example of computed loads using FLEX4 for a mean wind speed of 11m/s3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 106 106 Aerodynamics of Wind Turbines The tilt rotor moment in the main shaft, as shown in Figure 10.1, tries to tilt the nacelle over the tower. The yaw rotor moment tries to turn the nacelle on the tower. Sometimes the two bending moments at the root of the tower are also stated. Figure 10.2 shows an example of computed loads using the aero- elastic code FLEX4 at a mean wind speed of 11m/s. The wind speed at hub height and the power output are seen on the first graphs. The last four graphs show the corresponding time histories of the flapwise and edgewise bending moments and the tilt rotor and yaw rotor moments. The effect of gravity is clearly seen on the edgewise bending moment as a dominant sinusoidal variation, upon which is superimposed some small high frequency signal stemming from atmospheric turbulence. The flapwise bending is mostly influenced by the aerodynamic loads that vary with the turbulent wind field, and this signal is therefore more stochastic. References IEC 61400-1 (2004) ‘Wind turbines. Part 1: Design requirements’, CD, edition 3, second revision, IEC TC88-MT1 Sørensen, J. N. (1983) ‘Beregning af banekurver og kastelængder for afrevne vindmøllevinger’ (in Danish), AFM83-06, Department of Fluid Mechanics (AFM), Technical University of Denmark (DTU) 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 107 11 Beam Theory for the Wind Turbine Blade This section describes how a blade, whose outer contour is designed from aerodynamical considerations, is built to be sufficiently strong and stiff. In the past materials like wood, steel, aluminium, glass-fibre-reinforced plastics (GRPs) and carbon fibre reinforced plastics (CFRPs) have been used. The choice depends on many parameters such as strength, weight, stiffness, price and, very important for wind turbines, fatigue properties. The majority of wind turbine blades today are built using GRPs, and therefore a short description of a manufacturing process using this material is given. A negative mould for the upper part (suction side) and lower part (pressure side) of the blade is made. A thin film of so-called gelcoat is first laid in the moulds. The gelcoat gives a smooth white finish to the blades and therefore it is not necessary to paint the blades afterwards. Then a number of glass fibre mats are laid in. On each mat a layer of epoxy or polyester is rolled on to bind the mats into a hard matrix of fibres. The number of mats gives the thickness of the shell; typically a thin shell is made around the leading and trailing edges and a thick shell is made in the middle of the aerofoil. A section of such a blade is shown in Figure 11.1. To make the blade stronger and stiffer, so-called webs are glued on between the two shells before they are glued together. To make the trailing edge stiffer, foam panels can also be glued on before assembling the upper and lower parts. Because such a construction consists of different layers, it is often called a sandwich construction; a sketch is given in Figure 11.2 which can readily be compared to the real section in Figure 11.1. It is seen that the thick layer of mats and epoxy in the middle of the skin and the webs form a box-like structure. For structural analysis, the box-like structure, which is the most important structural part of the blade, acts like a main beam on which a thin skin is glued defining the geometry of the blade. Fixing a thin skin on a main beam (so the skin is not carrying loads, but merely gives an outer aerodynamic shape) is an alternative way that is sometimes used to construct a blade.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 108 108 Aerodynamics of Wind Turbines Figure 11.1 Section of an actual blade A blade can thus be modelled as a beam, and when the stiffnesses EI and GI v at different spanwise stations are computed, simple beam theory can be applied to compute the stresses and deflections of that blade. E is the modulus of elasticity, G is the modulus of elasticity for shear and I denotes different moments of inertia. In the next section a more elaborate explanation is given for the moments of inertia and it is shown how the stiffnesses can be computed for a wind turbine blade such as the one shown in Figure 11.1. The simple beam theory described here is found in almost any basic book on mechanics of materials (for example Gere and Timoshenko, 1972). Further, it is outlined how to compute the important structural parameters shown in Figure 11.3. Values of these parameters are necessary to compute the deflection of a blade for a given load or as input to a dynamic simulation using an aeroelastic code. Figure 11.2 Schematic drawing of a section of a blade3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 109 Beam Theory for the Wind Turbine Blade 109 Figure 11.3 Section of a blade showing the main structural parameters EI – bending stiffness about first principal axis; 1 EI – bending stiffness about second principal axis; 2 GI – torsional stiffness; v X – the distance of the point of elasticity from the reference point; E X – the distance of the centre of mass from the reference point; m X – the distance of the shear centre from the reference point; s  – the twist of the aerofoil section measured relative to the tip chord line; v – angle between chord line and first principal axis; +v – angle between tip chord line and first principal axis; The point of elasticity is defined as the point where a normal force (out of the plane) will not give rise to a bending of the beam. The shear centre is the point where an in-plane force will not rotate the aerofoil. If the beam is bent about one of the principal axes, the beam will only bend about this axis. As will be seen later, it is convenient to use the principal axes when calculating the blade deflection. Before continuing, some necessary definitions must be introduced. The following quantities are defined in terms of the reference coordinate system (X ,Y ) in Figure 11.4: R R • Longitudinal stiffness: EA = ∫ EdA. A X • Moment of stiffness about the axis X : ES = ∫ EY dA. R R A R Y • Moment of stiffness about the axis Y : ES = ∫ EX dA. R R A R3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 110 110 Aerodynamics of Wind Turbines Figure 11.4 Section of a blade 2 • Moment of stiffness inertia about the axis X : EIX = ∫ EY dA. R R A R 2 • Moment of stiffness inertia about the axis Y : EIY = ∫ EX dA. R R A R • Moment of centrifugal stiffness: EDXY = ∫ EX Y dA. R A R R From these definitions, the point of elasticity P = (X ,Y ) can be computed in E E E the reference system (X ,Y ) as: R R ESY R X = ––––– (11.1) E EA and: ESX R Y = ––––– (11.2) E EA For E and ρ constant the point (X ,Y) equals the centre of mass for the E E section, where ρ denotes the density of the material used. Now the moments of stiffness inertia and the moment of centrifugal stiffness are moved to the coordinate system (X’,Y’), which is parallel to the reference system (X ,Y ) R R and has its origin in the point of elasticity, using the formulaes: 2 2 EI ’ = ∫ E(Y’) dA = EIX –Y EA (11.3) X A R E3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 111 Beam Theory for the Wind Turbine Blade 111 2 2 EI = ∫ E(X’) dA = EIY –X EA (11.4) Y’ A E R ED = ∫ EX’Y’dA = EDXY –X Y EA. (11.5) X’Y’ A R E E It is now possible to compute the angle  between X’ and the first principal axis and the bending stiffness about the principal axes. The second principal axis is perpendicular to the first principal axis: 2ED X’Y’ 1 –1 –  = tan ––––––––––– (11.6) 2   EI – EI Y’ X’ EI = EI – ED tan  (11.7) 1 X’ X’Y’ EI = EI ED tan  (11.8) 2 Y’ X’Y’ The stress σ (x,y) in the cross-section from the bending moments about the two principal axes M and M and the normal force N is found from: X Y σ (x,y) =E(x,y)ε(x,y), (11.9) where the strain ε is computed from: M M N 1 2 ε(x,y) = –––– y – –––– x + –––– . (11.10) EIEIEA 1 2 σ, ε and N are positive for tension and negative for compression. The bending moments M and M and the normal force N are computed from the 1 2 loading of the blade, as is shown later. The main structural data are now determined. Since a wind turbine blade is very stiff in torsion, the torsional deflection is normally not considered. A complete description of how to compute the shear centre and the torsional rigidity is, however, given in Øye (1978). An example of results from Øye (1988), for the 30m blade used at the 2MW Tjæreborg wind turbine, is listed in Table 11.1. Table 11.1 shows that the position of the first principal axis, described by the angle +v, varies with the radius r. It is also seen that the position of the first principal axis is almost identical with the chord line since the angle v is small for most of the blade. 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 112 112 Aerodynamics of Wind Turbines Table 11.1 Main structural parameters for the Tjaereborg blade rEAEI EI GI mX X X v +v 1 2 v E m s 2 2 2 m GN MNm MNm MNm kg/m mm mm mm ° ° 1.8 36.0 12,000 12,000 7500 1700 0 0 0 0 0 3.0 6.14 1630 1725 362 330 2 2 0 5.4 14.4 4.5 5.821080 19408 389 54 159 11 0.94 9.44 6.0 5.10 623 1490 207 347 59 165 13 1.30 9.30 9.0 4.06 255 905 92.8 283 63 170 18 1.09 8.09 12.0 3.33 129 557 47.7 235 58 158 15 0.86 6.86 15.0 2.76 64.8 349 24.7 196 51 137 15 0.86 5.86 18.0 2.33 32.4 221 12.9 166 45 121 16 0.91 4.91 21.0 1.83 15.2 131 6.23 172 41 110 17 0.83 3.83 24.0 1.21 6.04 65.7 2.57 90.3 40 102 16 0.63 2.63 27.0 0.63 1.82 28.1 0.84 52.6 47 108 14 0.16 1.16 30.0 0.21 0.32 9.5 0.18 24.2 82 136 10 –0.52 –0.52 Deflections and Bending Moments A wind turbine blade such as that in Figure 11.5 can be treated as a technical beam as sketched in Figure 11.6. Note that the coordinate system used for the blade is different from the one shown in Figure 11.4. Figure 11.5 Wind turbine blade 323212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 113 Beam Theory for the Wind Turbine Blade 113 Figure 11.6 Technical beam If the external loadings, p and p , are known along the blade, the shear forces, y z T and T , and bending moments, M and M , can be found as: z y y z dT z ––– = –p (x) + m(x)u ¨ (x) (11.11) z z dx dT y ––– = –p (x) + m(x)u ¨ (x) (11.12) y y dx dM y ––– ––– = T (11.13) z dx dM z ––– ––– = –T. (11.14) y dx Equations (11.11–11.14) can be derived using Newton’s second law on an infinitesimal piece of the beam as shown in Figure 11.7. u ¨ is the acceleration and if the blade is in equilibrium, the last term (inertia term) on the right hand side of equations (11.11) and (11.12) is zero. Figure 11.7 Infinitesimal piece of the beam3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 114 114 Aerodynamics of Wind Turbines The bending moments can now be transformed to the principal axes as (here the y-axis is aligned with the tip chord): M = M cos( + v) – M sin( + v) (11.15) 1 y z M = M sin( + v) + M cos( + v), (11.16) 2 y z where +v is the angle between the y-axis and the first principal axis as shown in Figure 11.8. All other things being equal, it can be assumed that the first principal axis lies along the chord line, which is only true for a symmetric aerofoil. Note that  is negative for a normally twisted blade, but (+v) is assumed positive in equations (11.15), (11.16), (11.19) and (11.20). Figure 11.8 Orientation of principal axes The curvatures about the principal axes are from simple beam theory: M 1 κ = –— (11.17) 1 EI 1 M 2 κ = –—. (11.18) 2 EI 2 These curvatures are then transformed back to the y-axis and z-axis as: κ = –κ sin( + v) + κ cos( + v) (11.19) z 1 2 κ = κ cos( + v) + κ sin( + v) (11.20) y 1 23212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 115 Beam Theory for the Wind Turbine Blade 115 The angular deformations and thus the deflections are now calculated from: dθ y — = κ (11.21) y dx dθ z — = κ (11.22) z dx du z — –– – = – θ (11.23) y dx du y — –– – = θ . (11.24) z dx If the loads are given in discrete points as shown in Figure 11.9 and the loads are assumed to vary linearly between two points i and i 1, then equations (11.11) to (11.14) and equations (11.21) to (11.24) can be integrated to give the following numerical algorithm to calculate bending moments and deflections. Figure 11.9 Discretized cantilever beam Numerical Algorithm for Determining the Bending Moments and Deflection N N Boundary conditions: T = 0 and T = 0 y z For i = N to 2: 1 i–1 i i–1 i i i–1 – T = T + ( p + p )(x – x ) (11.25) y y y y 2 1 i–1 i i–1 i i i–1 – T = T + ( p + p )(x – x ) (11.26) z z z z 2 N N Boundary conditions: M = 0 and M = 0 y z3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 116 116 Aerodynamics of Wind Turbines For i = N to 2: i–1 i i i i–1 1 i–1 1 i i i–1 2 M = M – T (x – x ) – (–p + – p )(x – x ) (11.27) y y z 6 3 z z i–1 i i i i–1 1 i–1 1 i i i–1 2 M = M T (x – x ) (–p + –p )(x – x ) . (11.28) z z y 6 y 3 y For all points calculate M and M using equations (11.15) and (11.16). Then 1 2  and  are computed using equations (11.17–11.20). y z 1 1 Boundary conditions: θ = 0 and θ = 0 y z For i = 1,N–1: i 1 i 1 i+1 i i+1 i θ = θ + – (κ + κ )(x – x ) (11.29) y y 2 y y i+1 i 1 i+1 i i+1 i θ = θ + – (κ + κ )(x – x ). (11.30) z z 2 z z 1 1 Boundary conditions: u = 0 and u = 0 y z For i = 1,N–1: i+1 i i i+1 i 1 i+1 1 i i+1 i 2 u = u θ (x – x ) (–κ + –κ )(x – x ) (11.31) y y z 6 z 3 z i+1 i i i+1 i 1 i+1 1 i i+1 i 2 u = u – θ (x – x ) – (–κ + –κ )(x – x ) . (11.32) z z y 6 y 3 y The boundary conditions are for a cantilever beam and the inertia terms in equations (11.11) and (11.12) have been neglected for simplicity, but must be added for an unsteady computation. A Method to Estimate the First Flapwise, First Edgewise and Second Flapwise Eigenmodes The equations (11.25) to (11.32) can also be used to estimate the first few eigenmodes. An eigenmode is a free vibration without the presence of external loads; in other words equations (11.11) and (11.12) become: dT z —– = m(x)ü (x) (11.33) z dx dT y —– = m(x)ü (x). (11.34) y dx3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 117 Beam Theory for the Wind Turbine Blade 117 Since for an eigenmode the deflection is of the type u = Asin(ωt), the acceleration is proportional to the deflection as: 2 ü = –ω u, (11.35) where ω is the associated eigenfrequency. Using (11.35), equations (11.33) and (11.34) become: dT z 2 —– = –m(x)ω u (x) (11.36) z dx dT y 2 —– = –m(x)ω u (x). (11.37) y dx Comparing equations (11.11–11.12) with (11.36–11.37), it is seen that an eigenmode can be found using the static beam equations applying the external loads as: 2 p = m(x)ω u (x) (11.38) z z 2 p = m(x)ω u (x). (11.39) y y Since the deflections in equations (11.38) and (11.39) are not known, the equations must be solved iteratively, which will converge to the mode with the lowest eigenfrequency, known as the first flapwise mode. First, an initial deflection can be found using, for example, a constant loading for both the z and y directions. With this deflection the eigenfrequency is estimated at the tip as: N p 2 z ω = –––– (11.40) N N u m z and a new loading is computed in all discrete points as: i u z i 2 i p z = ω m ––––––––––––– – – . (11.41) –––––––––– N 2 N 2 (u ) + (u ) z y i u y i 2 i p y = ω m ––––––––––––– – – . (11.42) –––––––––– N 2 N 2 (u ) + (u ) z y Note that the loading is normalized with the tip deflection to ensure that the tip deflection becomes 1 in the next iteration. A new deflection is now3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 118 118 Aerodynamics of Wind Turbines computed using the loading from equations (11.41) and (11.42). The proced- ure is repeated a few times until the eigenfrequency ω becomes constant. 1f 1f Now the deflection shape, u and u , of the first flapwise eigenmode as z y sketched in Figure 11.10 is known. To find the first edgewise deflection the procedure can be used again. However, some modification is needed to make sure the method does not converge towards the first flapwise mode. Every time a new deflection, u (x) y and u (x), is computed, it is necessary to subtract the part that contains the z first flapwise mode : le lf u = u – const ·u (11.43) z z 1 z le lf u = u – const ·u . (11.44) y y 1 y To determine the constant, the following orthogonality constraint must be satisfied: R R lf le lf le ∫ u mu dx+∫ u mu dx = 0 (11.45) z z y y 0 0 Combining equations (11.43) and (11.44) with equation (11.45) yields an expression for the constant: R R lf lf ∫ u mu dx+∫ u mu dx (11.46) z z y y 0 0 const = –––––––––––––––––– . 1 R R lf lf lf lf ∫ u mu dx+∫ u mu dx z z y y 0 0 Thus to force the method to converge to the first edgewise eigenmode it is necessary every time a new deflection has been computed to remove the part from the first flapwise eigenmode using equations (11.43), (11.44) and (11.46). After a few iterations the deflection shape of the first edgewise le le eigenmode, u and u , as sketched in Figure 11.11, is known. z y Finally, the second flapwise eigenmode can be found similarly. Now it is also necessary to substract not only the part that contains the first flapwise eigenmode but also the part from the first edgewise eigenmode using the following orthogonality constraint: R R le 2 f le 2 f ∫ u mu dx+∫ u mu dx = 0. (11.47) z z y y 0 0 Every time a new deflection has been computed, equations (11.43), (11.44) and (11.46) are used to remove the part from the first flapwise eigenmode. Next equations (11.48–11.50) are applied to remove the part from the first3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 119 Beam Theory for the Wind Turbine Blade 119 edgewise eigenmode fulfilling constraint (11.47): 2f 1e u = u – const ·u (11.48) z z 2 z 2f 1e u = u – const ·u (11.49) y y 2 y R R le le ∫ u mu dx + ∫ u mu dx z z y y 0 0 const = ––––––––––––––––––– (11.50) 2 R R le le le le ∫ u mu dx + ∫ u mu dx z z y y 0 0 After a few iterations the deflection shape of the second flapwise eigenmode, 2f 2f u and u , as sketched in Figure 11.12, is known. z y Figure 11.10 First flapwise eigenmode (1f) Figure 11.11 First edgewise eigenmode (1e) The eigenmode estimation can be checked by comparing with the analytical solutions for a cantilever beam having constant properties along the beam (see, for example, Craig, 1981):3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 120 120 Aerodynamics of Wind Turbines Figure 11.12 Second flapwise eigenmode (2f) 1 / 3.516 EI 2 ω = ––––– ––– (11.51) 1 2  L m 1 / 22.03 EI 2 ω = ––––– ––– (11.52) 2 2  L m ω and ω are the first and second eigenfrequencies respectively and L is the 1 2 length of the beam. Applying this algorithm on a beam with constant stiffness 2 EI = 1Nm , a mass distribution m = 1kg/m, a length L = 1m and using 11 points yields ω = 3.513rad/s, which is close to the analytical solution of 1 3.516 rad/s. Using 11 points, the second eigenfrequency becomes 22.273rad/s and 22.044rad/s for 51 points, which should be compared to the analytical value of 22.03rad/s. The solution for 51 points is shown in Figure 11.13. Figure 11.13 Second eigenmode using 50 elements3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 121 Beam Theory for the Wind Turbine Blade 121 In an accurate computation the effect of the centrifugal force on a deflected blade must be included. In Figure 11.14, the total force F from the x centrifugal acceleration at r = x is shown. Figure 11.14 Total centrifugal force F at a spanwise position x x F can be computed by integrating the incremental contribution dF from x to x x the tip of the blade R. R R 2 F = ∫ dF = ∫ x’ω mdx’, (11.53) x x x x3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 122 122 Aerodynamics of Wind Turbines where ω is the angular velocity of the rotor, m the mass per length, dx’ an incremental part of the blade and the x’ radial distance to the incremental part. The projection of F normal to the rotor blade is (see Figure 11.14): x F = F sinθ  Fθ. (11.54) c x x This force corresponds to a loading pc (force per length): dF dθ dF c x p =––– = F ––– + θ ––– (11.55) c x dx dx dx 2 F decreases for increasing x, thus dF /dx = –xω m and equation (11.55) x x becomes: dθ 2 p = F ––– – mxω θ (11.56) c x dx p is the aerodynamic load computed by BEM method and p +p is the aerodynamic plus aero aero c the centrifugal load. Figure 11.15 The computed load about the first principal axis for the 2MW Tjæreborg machine at V = 10m/s with and without the centrifugal loading o 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 123 Beam Theory for the Wind Turbine Blade 123 Figure 11.16 The computed bending moment about the first principal axis for the 2MW Tjæreborg machine at V = 10m/s with and without the centrifugal loading o The extra loading from equation (11.56) should be added to the load calculated from the aerodynamics and the net result is a redistribution of the loading, which results in a reduction in the flapwise bending moment. Since pc is a function of the angular deformation θ and the curvature dθ/dx, which again is a function of the loads, it is necessary to iterate a few times in order to have the deflection including the effect of the centrifugal acceleration. In other words, first the deflections and the angular deformations are computed by using only the aerodynamic loads. Then the extra loading is computed using equations (11.53) and (11.56). A new loading is found by adding this extra loading to the original aerodynamic loading and a new set of deflect- ions and angular deformations are computed. Hereafter the centrifugal loading is calculated again using equation (11.56) and the procedure is repeated a few times until a stationary solution is found. Figure 11.15 shows the effect of including the centrifugal loading on the 2MW Tjæreborg machine at a wind speed of 10m/s. It is seen that the centrifugal loading reduces the loads at the tip but increases the loads at the root. The corres- ponding reduction in the bending moments is seen in Figure 11.16.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 124 124 Aerodynamics of Wind Turbines References Craig, R. R. Jr (1981) Structural Dynamics, John Wiley & Sons, Hoboken, NJ Gere, J. M. and Timoshenko, S. P. (1972) Mechanics of Materials, van Nostrand Reinhold Company, New York Øye, S. (1978) ‘Beregning af tværsnitskonstanter for bjælker med inhomogent tværsnit’ (in Danish), AFM notat VK-45-781201, Department of Fluid Mechanics, Technical University of Denmark (DTU) Øye, S. (1988) ‘Projekt K 30m Glasfibervinge, Teknisk beskrivelse’ (in Danish), AFM 88- 12, Department of Fluid Mechanics, Technical University of Denmark (DTU)3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 125 12 Dynamic Structural Model of a Wind Turbine The main purpose of a structural model of a wind turbine is to be able to determine the temporal variation of the material loads in the various com- ponents in order to estimate the fatigue damage. Further, a dynamic system is used when analysing the stability of the wind turbine design, including perhaps the control system. To calculate the deflections and velocities of the various components in the wind turbine in the time domain, a structural model including the inertia terms is needed. Then the dynamic structural response of the entire construction can be calculated subject to the time dependent load found using an aerodynamic model, such as the BEM method. For offshore wind turbines, wave loads and perhaps ice loads on the bottom of the tower must also be estimated. One way of setting up the structural model, based on the principle of virtual work, is presented here in detail. However, more formal ‘finite element’ methods have also been used in different aeroelastic codes. The velocity of the vibrating wind turbine construction must be subtracted when calculating the relative velocity seen locally by the blade as shown in equation (9.32). The loads therefore depend on the deflections and velocities of the structure, which again depend on the loads. The structural and aerodynamic models are therefore highly coupled and must be solved together in what is known as an aeroelastic problem. Principle of Virtual Work and Use of Modal Shape Functions The principle of virtual work is a method to set up the correct mass matrix, M, stiffness matrix, K, and damping matrix, C, for a discretized mechanical system as: · M¨ x + Cx + Kx = F , (12.1) g3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 126 126 Aerodynamics of Wind Turbines where F denotes the generalized force vector associated with the external g loads, p. Equation (12.1) is of course nothing but Newton’s second law, assuming linear stiffness and damping, and the method of virtual work is nothing but a method that helps in setting up the correct mass, stiffness and damping matrices for a multi-body system, which is especially well suited for a chain system. Knowing the loads and appropriate conditions for the velocities and the deformations, equation (12.1) can be solved for the accelerations, from which the velocities and deformations can be estimated for the next time step. The number of elements in x is called the number of degrees of freedom, DOF, and the higher this number the more compu- tational time is needed in each time step to solve the matrix system. Use of modal shape functions is a tool to reduce the number of degrees of freedom and thus reduce the size of the matrices to make the computations faster per time step. A deflection shape in this method is described as a linear combination of a few but physically realistic basis functions, which are often the deflection shapes corresponding to the eigenmodes with the lowest eigenfrequencies. For a wind turbine such an approach is suited to describe the deflection of the rotor blades and the assumption is that the combination of the power spectral density of the loads and the damping of the system do not excite eigenmodes associated with higher frequencies. In the com- mercially available and widely used aeroelastic simulation tool FLEX (see, for example, Øye, 1996), only the first 3 or 4 (2 flapwise and 1 or 2 edge- wise) eigenmodes are used for the blades and results are in good agreement with measurements, indicating the validity of the assumption. First one has to decide on the DOFs necessary to describe a realistic deformation of a wind turbine. For instance, in FLEX4 17–20 DOFs are used for a three- bladed wind turbine: 3–4 DOFs per blade as described above; the deformation of the shaft is described using 4 DOFs (1 for torsion, 2 hinges just before the first bearing with associated angular stiffness to describe bending and 1 for pure rotation); 1 DOF to describe the tilt stiffness of the nacelle; and finally 3 DOFs for the tower (1 for torsion, 1 mode in the normal direction of the rotor and 1 in the lateral direction). The values in the vector x describing the deformation of the construction, x , are known as the general coordinates. To each generalized coordinate is i associated a deflection shape, u , that describes the deformation of the con- i struction when only x is different from zero and typically has a unit value. i The element i in the generalized force corresponding to a small displacement in DOF number i, dx , is calculated such that the work done by the general- i ized force equals the work done on the construction by the distributed external loads on the associated deflection shape:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 127 Dynamic Stuctural Model of a Wind Turbine 127 F dx = ∫p·u dS, (12.2) g,i i i s where S denotes the entire system. Please note that the generalized force can be a moment and that the displacement can be angular. All loads must be included, in other words also gravity and inertial loads such as Coriolis, centrifugal and gyroscopic loads. The non-linear centrifugal stiffening can be modelled as equivalent loads calculated from the local centrifugal force and the actual deflection shape as shown in the previous chapter. The elements in the mass matrix, m , can be evaluated as the generalized force from the i,j inertia loads from an unit acceleration of DOF j for a unit displacement of DOF i. The elements in the stiffness matrix, k , correspond to the generalized i,j force from an external force field which keeps the system in equilibrium for a unit displacement in DOF j and which then is displaced x = 1. The elements i in the damping matrix can be found similarly. For a chain system the method of virtual work as described here normally gives a full mass matrix and diagonal matrices for the stiffness and damping. For one blade rigidly clamped at the root (cantilever beam), it is relatively easy to estimate the 1f 1e lowest eigenmodes (first flapwise u , first edgewise u and second flapwise 2f u ), using, for example, the iterative procedure from the previous chapter. This description of the principle of virtual work might seem very abstract, but will hopefully become clearer after using it on two examples, one of two discrete masses connected by springs and dampers and another for an isolated blade (see later). If the structural system comprises a system of continuous mass distri- butions, such as a system of beams, equation (12.1) is the result of discretizing the system, since in reality such a system has an infinite number of DOFs. The elements in the mass, stiffness and damping matrices depend on the system and, in the case of a continuous system, also of the discre- tization. If the right hand side of equation (12.1) is 0 the system is said to perform its natural motion. · Provided that the deflections, x, and velocities, x, are known, equation (12.1) can alternatively be written as: · · M¨ x = F – Cx – Kx = f(x, x, t), (12.3) g where the function f in general is non-linear. Non-linarity can come, for example, from non-linear loads p or from aerodynamic damping (see later). A non-linear system can be treated as a linearized eigenvalue approach or as a full non-linear time domain approach. Only the latter method is treated in this text.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 128 128 Aerodynamics of Wind Turbines n Knowing the right hand side of equation (12.3) at time t = n∆ t, the n acceleration at time t is found solving the linear system of equations: n –1 n n n · x ¨ = M f(x , x , t ). (12.4) n n n n · Knowing the accelerations, x ¨ , the velocities, x , and positions, x , at time t , n+1 · a Runge-Kutta-Nyström scheme can be used to estimate the velocities, x , n+1 n+1 n+1 n+1 n+1 n+1 · and positions, x , at t . New loads, p (x , x , t ), can be calculated using, for example, an unsteady BEM method and thus equation (12.4) can be updated and a new time step can be performed. This can be continued until a sufficient time period has been simulated. · Box 12.1 The Runge-Kutta-Nyström integration scheme of x ¨ = g(t, x, x) ∆ t n A = –– x ¨ 2 ∆ t n 1 · – b = –– (x + A) 2 2 1 ∆ t – n+ n n · 2 B = –– g(t , x + b, x + A) 2 1 ∆ t – n+ n n · 2 C = –– g(t , x + b, x + B) 2 n · d = ∆ t (x + C) ∆ t n+1 n n · D = –– g(t , x + d, x + 2C) 2 and the final update: n+1 n t = t + ∆ t n+1 n n 1 · – x = x + ∆ t(x + (A + B + C)) 3 n+1 n 1 · · – x = x + (A+2B +2C+ D) 3 n+1 n+1 n+1 n+1 · x ¨ = g(t , x , x )3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 129 Dynamic Stuctural Model of a Wind Turbine 129 SDOF The simplest dynamic system is called SDOF (single degree of freedom) and is comprised of only one concentrated mass. Imagine, for example, a mass which is hung up in a spring as shown in Figure 12.1. Figure 12.1 SDOF system with no damping Statically the spring will be stretched until the spring force equals the weight Mg. The dynamic equation is: Mü + ku = 0, (12.5) where k is the spring constant and u the displacement from the equilibrium position. The well known analytical solution to this system is: · u o u = u cos(ωNt) + –– sin(ω t). (12.6) o N ω N ––– – · u is the pertubation at t =0, u is the velocity at t = 0 and ω = k/M is the o o N eigenfrequency of the system. The system performs a simple harmonic vibration and the time for one cycle is T = 2/ω . This simple SDOF system N can be used, for example, to check the implementation, accuracy and stability of the Runge-Kutta-Nyström method.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 130 130 Aerodynamics of Wind Turbines Aerodynamic Damping The aerodynamic forces from the flow past a structure can cause damping that may be negative and thus amplify vibrations. Below is a description of the aerodynamic damping using a simple SDOF system, but now including aerodynamic forces, as sketched in Figure 12.2. The system comprises an aerofoil in a wind tunnel, which is mounted with a geometrical angle of attack  on a spring that allows the aerofoil to move up and down. Imagine g · now that the blade is moving downwards with a velocity x. The aerofoil then · feels an extra velocity, equal to x but in the opposite direction, giving a component from below that, when added to the velocity in the wind tunnel V, gives the relative velocity (see Figure 12.2). Provided that the aerofoil o data, C (), are known, the flow angle φ, the angle of attack and the force in l the x-direction can be estimated as: · x tan φ = –– (12.7) V o  =  + φ (12.8) g 1 2 F = –ρV AC ()cos φ (12.9) x 2 rel l When the aerofoil moves downwards, the angle of attack increases, and when it moves upwards, the angle of attack is decreased according to equations (12.7) and (12.8). This changes the lift coefficient: ∂C l C = C + ––– – ∆ , (12.10) l l,o ∂ where C is the previous value of the lift coefficient and C the new lift l,o l coefficient at the higher angle of attack. If the aerofoil is moving downwards and the slope ∂C /∂ is positive, the lift coefficient and thus the aerodynamic l force is increased and, since this increased force is in the opposite direction to the motion, the vibration is damped. The same argument is valid when the blade is moving upwards. If the slope is negative, however, as in stalled conditions, the aerodynamic damping is negative. More formally, the work done by the aerodynamic forces on the section during one cycle can be calculated as: W = – F·dx. (12.11)   ∫ F is the aerodynamic force on the aerofoil and x the displacement. If the work is positive the section is positively damped, and if the work is negative the3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 131 Dynamic Stuctural Model of a Wind Turbine 131 Figure 12.2 SDOF with lift section is negatively damped. To have the correct slopes for the lift and drag coefficients is therefore essential for correctly predicting the stability. For these types of vibrations, dynamic stall models as described in Chapter 9 are thus very important. In some cases the dynamic stall model effectively increases the slope ∂C /∂, and thus the stability. In other words compu- l tations may over-predict the oscillations if static aerofoil data are used. Examples of Using the Principle of Virtual Work Two examples are given here of how to use the principle of virtual work to create the mass, stiffness and damping matrices in equation (12.1). The first example is a 2-DOF system comprising two masses connected by springs and dampers as shown in Figure 12.3. The second example is a wind turbine blade using modal shapes in order to reduce the number of DOFs. 2-DOF system The mass and stiffness matrices together with the generalized force vector for the 2-DOF system shown in Figure 12.3 will be set up using the principle of virtual work and generalized coordinates. 3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 132 132 Aerodynamics of Wind Turbines Figure 12.3 Example of a 2-DOF system The methodology is as follows: first, the generalized coordinates, (x ,x ), are 1 2 defined as the relative displacements between the two masses: x = u 1 1 (12.12) x = u – u 2 2 1 where u and u are the displacements of mass 1 and 2, respectively. The 1 2 generalized force vector, F , is found by the principle of virtual work as the g work done by the external forces, F and F , for a displacement of one of the 1 2 generalized coordinates keeping the other(s) zero. The first component is found for x = 1 and x = 0 as F = F +F since in this case both masses moves 1 2 g,1 1 2 a unit length. The second component is found putting x = 0 and x = 1 as F 1 2 g,2 = F since in this case only mass 2 moves a unit length. The mass matrix is 2 found by specifying a unit acceleration of one of the generalized coordinates keeping the other(s) zero and replacing the external forces by inertia forces, in other words mass times acceleration. For ¨ x = 1 and ¨ x = 0 both masses 1 1 2 and 2 accelerate with a unit acceleration, in other words the corresponding inertial forces becomes F = m and F = m2 and the generalized force based 1 1 2 on these become F = m +m and F = m . This gives the first column in the g,1 1 2 g,2 2 mass matrix as shown below: m m 1 m m + m 11 12 11 1 2 = = . (12.13) m m 0 m m 21 22 21 23212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 133 Dynamic Stuctural Model of a Wind Turbine 133 The second column is found by specifying ¨ x = 0 and ¨ x = 1; in other words 1 2 only mass 2 accelerates and the inertia forces become F = 0 and F = m . The 1 2 2 generalized force based on these inertia forces becomes F = m and F = m , g,1 2 g,2 2 yielding: m m 0 m m 11 12 12 2 = = . (12.14) m m 1 m m 21 22 22 2 The mass matrix therefore becomes: m + m m 1 2 2 M = (12.15) m m 2 2 The first column in the stiffness matrix can be found using the necessary generalized force to obtain a unit static displacement of the first generalized coordinate, i.e. x = 1 and x = 0. The necessary external forces for this unit 1 2 displacement are F = k and F = 0, where k is the spring constant for spring 1 1 2 1 1. The corresponding generalized force is F = F +F = k and F = F = 0 and g,1 1 2 1 g,2 2 the first column in the stiffness matrix becomes: k k 1 k k 11 12 11 1 = = . (12.16) k k 0 k 0 21 22 21 The second column is found as the necessary generalized force to obtain a unit static displacement of the second generalized coordinate, i.e. x = 0 and 1 x = 1. The necessary external forces for this unit displacement are F = –k 2 1 2 and F = k , where k is the spring constant for spring 2. The corresponding 2 2 2 generalized force is F = F +F = 0 and F = F = k , yielding: g,1 1 2 g,2 2 2 k k 0 k 0 11 12 12 = = . (12.17) k k 1 k k 21 22 22 2 The stiffness matrix therefore becomes: k 0 1 K = . (12.18) 0 k 2 The damping matrix can be found similarly. The first column is the generalized force necessary to obtain a unit velocity of the first generalized · · coordinate, i.e. x = 1 and x = 0. The result is similar to the stiffness matrix 1 2 and the damping matrix becomes:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 134 134 Aerodynamics of Wind Turbines d 0 1 C = . (12.19) 0 d 2 where d and d are the coefficients of viscous damping for the two dampers. 1 2 In a chain formulation like the above example the stiffness and damping matrices become diagonal matrices, whereas the mass matrix becomes full. The full system is thus determined as: · m + m m x ¨ d 0 x k 0 x F + F 1 2 2 1 1 1 1 1 1 2 ++ = (12.20) · m m x ¨ 0 d x 0 k x F 2 2 2 2 2 2 2 2 To verify the described method Newton’s second law is applied on each mass: · · · m ü = –k u –d u +k (u –u ) + d (u –u )+F (12.21) 1 1 1 1 1 1 2 2 1 2 2 1 1 · · m ü = –k (u –u ) –d (u –u )+F. (12.22) 2 2 2 2 1 2 2 1 2 Replacing u by x and u by x +x and replacing equation (12.21) with 1 1 2 1 2 equation (12.21) + equation (12.22) makes equations (12.21) and (12.22) identical to equation system (12.20). It is straightforward to add more masses. Dynamic system for blade The methodology of generalized coordinates is now applied to a wind turbine blade as shown in Figure 11.5. Assume that we know the normalized eigen- modes, with a maximum tip deflection of 1m. The first eigenmodes can then be calculated as shown in the previous chapter and sketched in Figures 11.10–11.12. It is assumed that a deformation of a blade can be described as a linear combination of these three modes as: lf le 2f u (x) = x ·u (x) + x ·u (x) + x ·u (x) (12.23) z 1 z 2 z 3 z and: lf le 2f u (x) = x ·u (x) + x ·u (x) + x ·u (x) (12.24) y 1 y 2 y 3 y The deflection shape can thus be described by three parameters only, which are denoted by the generalized coordinates x , x and x . Since the modes are 1 2 33212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 135 Dynamic Stuctural Model of a Wind Turbine 135 constant, the velocities and accelerations along the blade are: lf le 2f · · · · u (x) = x ·u (x) + x ·u (x) + x ·u (x) (12.25) z 1 z 2 z 3 z lf le 2f · · · · u (x) = x ·u (x) + x ·u (x) + x ·u (x) (12.26) y 1 y 2 y 3 y and: lf le 2f ü (x) = x ¨ ·u (x) + x ¨ ·u (x) + x ¨ ·u (x) (12.27) z 1 z 2 z 3 z lf le 2f ü (x) = x ¨ ·u (x) + x ¨ ·u (x) + x ¨ ·u (x) (12.28) y 1 y 2 y 3 y As indicated by equation (12.2), the generalized force for each mode is the work done on this mode by the external loads, p (x) and p (x), without z y contribution from the other modes, i.e.: 1f l f F = p (x)u (x)dx + ∫p (x)u (x)dx (12.29) ∫ g,1 z z y y 1e l e F = p (x)u (x)dx + ∫p (x)u (x)dx (12.30) ∫ g,2 z z y y and: 2f 2f F = p (x)u (x)dx + ∫p (x)u (x)dx (12.31) ∫ g,3 z z y y The first column of the mass matrix is found by evaluating the generalized force from external forces corresponding to the inertia forces for a unit acceleration of the first degree of freedom and the others set to 0, i.e. (x ¨ , x ¨ , 1 2 x ¨ ) = (1, 0, 0). Using equations (12.27) and (12.28), the inertia loads become 3 lf lf ( p , p ) = (mü , mü ) = (mu , mu ) and thus: y z y z y z lf lf lf lf u (x)m(x)u (x)dx + u (x)m(x)u (x)dx ∫ ∫ z z y y m GM 11 1 lf le lf le m = u (x)m(x)u (x)dx + u (x)m(x)u (x)dx = 0 . (12.32) ∫ ∫ 21 z z y y m lf 2f lf 2f 0 31 u (x)m(x)u (x)dx + u (x)m(x)u (x)dx ∫ ∫ z z y y The first element is sometimes denoted by the first generalized mass GM ; the 1 two other integrals are 0 due to the orthogonality constraints between eigen- modes. The second column of the mass matrix is found by evaluating the generalized force from external forces corresponding to the inertia forces for a unit acceleration of the second degree of freedom and the others set to 0,3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 136 136 Aerodynamics of Wind Turbines i.e. (x ¨ , x ¨ , x ¨ ) = (0, 1, 0). Using equations (12.27) and (12.28) the inertia loads 1 2 3 le le become ( p , p ) = (mü , mü ) = (mu , mu ) and thus: y z y z y z le lf le lf u (x)m(x)u (x)dx + u (x)m(x)u (x)dx ∫ ∫ z z y y m 0 12 le le le le m = u (x)m(x)u (x)dx + u (x)m(x)u (x)dx = GM . (12.33) ∫ ∫ 22 z z y y 2 m le 2f le 2f 0 32 u (x)m(x)u (x)dx + u (x)m(x)u (x)dx ∫ ∫ z z y y The second element is sometimes denoted the second generalized mass; the two other integrals are 0 due to the orthogonality constraints between eigenmodes. The third column of the mass matrix is found by evaluating the general- ized force from external forces corresponding to the inertia forces for a unit acceleration of the third degree of freedom and the others set to 0, i.e. (x ¨ , x ¨ , 1 2 x ¨ ) = (0, 0, 1). Using equations (12.27) and (12.28), the inertia loads become 3 2f 2f (p , p)=(mü , mü)=(mu , mu ) and thus: y z y z y z 2f lf 2f lf u (x)m(x)u (x)dx + u (x)m(x)u (x)dx ∫ ∫ z z y y m 0 13 2f le 2f le m = u (x)m(x)u (x)dx + u (x)m(x)u (x)dx = . (12.34) ∫ ∫ 23 z z y y 0 m 2f 2f 2f 2f 33 u (x)m(x)u (x)dx + u (x)m(x)u (x)dx GM ∫ ∫ z z y y 3 The third element is sometimes denoted the third generalized mass; the two other integrals are 0 due to the orthogonality constraints between eigenmodes. The first column in the stiffness matrix can be found using the generalized force necessary to obtain a unit static displacement of the first generalized coordinate, i.e. (x , x , x ) = (1, 0, 0). In this case the deflection, according 1 2 3 1f 1f to equations (12.23) and (12.24), is identical to u , u . The loads yielding y z this deflection, according to equations (11.38) and (11.39), are ( p , p)= y z 2 1f 2 1f (mω u , mω u ,) , where ω is the eigenfrequency associated with the first 1 y 1 z 1 flapwise eigenmode; applying this in the equations for the generalized force (12.29–12.31) yields: 2 lf lf 2 lf lf ω u mu dx + ω u mu dx 2 ∫ ∫ 1 z z 1 y y k ω GM 11 1 1 2 lf le 2 lf le k = ω u mu dx + ω u mu dx = 0 . (12.35) ∫ ∫ 21 1 z z 1 y y k 2 lf 2f 2 lf 2f 0 31 ω u mu dx + ω u mu dx ∫ ∫ 1 z z 1 y y The two last integrals are 0 due to the orthogonality constraint of the eigen- modes. Similarly it can be shown that:3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 137 Dynamic Stuctural Model of a Wind Turbine 137 k 0 12 2 k = ω GM (12.36) 22 2 2 k 0 32 and: 13 k 0 k = 0 . (12.37) 23 2 k ω GM 33 3 3 ω2 and ω3 are the eigenfrequencies associated with the first edgewise and second flapwise eigenmodes respectively. The generalized masses GM , GM 1 2 and GM are defined in equations (12.32–12.34). 3 Neglecting the structural damping, the equation for one blade becomes: 2 GM 00 x ¨ ω GM00 x F 1 1 1 1 1 g,1 2 0 GM 0 x ¨ +0 ω GM 0 x = F . (12.38) 2 2 2 2 2 g,2 2 00 GM x ¨ 00 ω GM x F 3 3 3 3 3 g,3 Structural damping terms can be modelled as: δ 1 ω GM–– 00 1 1 π δ 2 C = 0 ω GM–– 0 . (12.39) 2 2 π δ 3 00 ω GM–– 3 3 π where δ is logarithmic decrement associated with mode i. The system for the i beam in equation (12.38) comprises three uncoupled differential equations. This is a result of the orthogonality constraints on the eigenmodes used; it is not a general result of the principle of virtual work. For instance, when the method is used on a whole wind turbine construction the equation of motion becomes fully coupled. FEM Models Alternatively, a more formal finite element method (FEM) can be used to set up the dynamic structural model in the form of equation (12.1). In Ahlström (2005) and Schepers (2002) a long list of various aeroelastic codes are compiled, and many of the recently developed codes have used the FEM3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 138 138 Aerodynamics of Wind Turbines approach. However, the number of DOFs in a FEM discretization will be considerably larger than, for example, that from combining the principle of virtual work with mode shapes as shown in the text above for a single blade. Therefore the computational time required to calculate a time history of a certain length is much larger when a FEM approach is taken. In this book there will be no attempt to use the FEM approach for modelling a wind turbine, but more details and further reference can be found in Hansen et al (2006), for example. References Ahlström, A. (2005) ‘Aeroelastic simulation of wind turbine dynamics’, Doctoral thesis in structural mechanics, Royal Institute of Technology, KTH, Stockholm, Sweden Hansen, M. O. L., Sørensen, J. N., Voutsinas, S., Sørensen, N. and Madsen, H. Aa. (2006) ‘State of the art in wind turbine aerodynamics and aeroelasticity’, Progress in Aerospace Sciences, vol 42, no 4, pp285–330 Øye, S. (1996) ‘FLEX4 simulation of wind turbine dynamics’ in Proceedings of 28th IEA Meeting of Experts Concerning State of the Art of Aeroelastic Codes for Wind Turbine Calculations (available through International Energy Agency) Schepers, J. G. (2002) ‘Verification of European wind turbine design codes, VEWTDC; final report’, technical report ECN-C-01-055, Energy Research Centre of The Netherlands3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 139 13 Sources of Loads on a Wind Turbine The three most important sources of the loading of a wind turbine are 1 gravitational loading; 2inertial loading; and 3 aerodynamic loading. Gravitational Loading: The Earth’s gravitational field causes a sinusoidal gravitational loading on each blade, as indicated in Figure 13.1. Figure 13.1 The loading caused by the Earth’s gravitational field When the blade is in position 1 in Figure 13.1 (down-rotating) the blade root at the trailing edge side is exposed to tensile stress and the leading edge side of the blade root is exposed to compressive stress. In position 2 (up-rotating) the trailing edge side of the blade root is exposed to compressive stress and the leading edge side of the blade root is exposed to tensile stress. Thus gravity is responsible for a sinusoidal loading of the blades with a frequency3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 140 140 Aerodynamics of Wind Turbines corresponding to the rotation of the rotor often denoted by 1P. This loading is easily recognized in Figure 10.2 in the time series of the edgewise bending moment. Note that a wind turbine is designed to operate for 20 years, which means that a machine operating at 25 rpm will be exposed to 20 365 24 8 60 25 = 2.6 10 stress cycles from gravity. Since a wind turbine blade might weigh several tons and be more than 30m long, the stresses from the gravity loading are very important in the fatigue analysis. Inertial Loading Inertial loading occurs when, for example, the turbine is accelerated or decelerated. An example is the braking of the rotor, where a braking torque T is applied at the rotor shaft. A small section of the blade will feel a force dF in the direction of the rotation as indicated in Figure 13.2. Figure 13.2 Loading caused by braking the rotor3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 141 Sources of Loads on a Wind Turbine 141 The size of dF is found from: · dF = ωrmdr, (13.1) where m is the mass per length of the blade, r the radius from the rotational · axis to the section and dr the size of the small section; ω = dω/dt can be found from: dω I ––– = T, (13.2) dt where I is the moment of inertia of the rotor. The mü terms in equations (11.11) and (11.12) are also inertia loads stemming from local accelerations. Another inertial loading stems from the centrifugal force on the blades. In order to reduce the flapwise bending moment, the rotor can be coned backwards with a cone angle of θ as shown in Figure 13.3. cone Figure 13.3 Effect of coning the rotor3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 142 142 Aerodynamics of Wind Turbines The centrifugal force acting on the incremental part of the blade at a radius r 2 from the rotational axis as shown in Figure 13.3 is F = ω rmdr, where M is c the mass of the incremental part and ω the angular velocity of the rotor. Due to the coning the centrifugal force has a component in the spanwise direction of the blade, Fcosθ , and a component normal to the blade, Fsinθ , as c cone c cone shown in Figure 13.3. The normal component gives a flapwise bending moment in the opposite direction to the bending moment caused by the thrust and thus reduces the total flapwise bending moment. Aerodynamic Loading The aerodynamic loading is caused by the flow past the structure, in other words the blades and the tower. The wind field seen by the rotor varies in space and time due to atmospheric turbulence as sketched in Figure 13.4. Figure 13.4 Sketch of turbulent inflow seen by wind turbine rotor As also seen in Figure 13.4, the wind field is characterized by shear; in other words the mean wind speed increases with the height above the ground. For neutral stability this shear may be estimated as: V (x) 1n(x/z ) 10min o –––––– = –––––– . (13.3) V (h) 1n(h/z ) 10min o3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 143 Sources of Loads on a Wind Turbine 143 V (x) is the time averaged value for a period of 10 minutes at a height x 10min above the ground. V (h) is the time averaged value at a fixed height h and 10min z is the so-called roughness length. Alternatively the wind shear can be given o by an exponent as in equation (9.33). The roughness length depends on the –4 surface characteristics and varies from 10 m over water to approximately 1m in cities. Values of z can be found in Troen and Petersen (1989) and are o summarized in Table 13.1. Table 13.1 Roughness length table z m Terrain surface characteristics o 1.0 City 0.8 Forest 0.2Many trees and bushes 0.1 Farmland with closed appearance 0.05 Farmland with open appearance 0.03 Farmland with very few buildings, trees, etc. –3 5 10 Bare soil –3 1 10 Snow surface –4 3 10 Sand surface (smooth) –4 1 10 Water areas Source: Troen and Petersen (1989). Wind shear gives a sinusoidal variation of the wind speed seen by a blade with a frequency corresponding to the rotation of the rotor 1P. The turbulent fluctuations superimposed on the mean wind speed also produce a time variation in the wind speed and thus in the angle of attack. In order to simulate the behaviour of a wind turbine using an aeroelastic code, it is therefore necessary first to generate a realistic wind field as shown in Chapter 14. The tower also contributes to variation in the inflow; for an upwind rotor this can be calculated using equations (9.34–9.36). A wind turbine might operate in yaw if, for example, the direction of the wind is not measured correctly or in the event of a malfunctioning yaw system. In this case, the reduced wind speed u at the rotor plane has a component normal to the rotor, ucosθ , and tangential to the rotor, usinθ . yaw yaw If the blade at the top position in Figure 13.5, corresponding to θ = 0°, wing moves in the same direction as the wind, the relative rotational speed ωr(1+a’) is reduced by usinθ ; at the bottom position, θ = 180°, it will yaw wing move against the wind speed and the relative rotational speed is increased by usinθ . Further, the axial induced velocity, and thus the axial velocity u, is yaw not constant in an annular element of the rotor disc (see equation (9.26)).3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 144 144 Aerodynamics of Wind Turbines Figure 13.5 Rotor plane showing the azimuthal position of a blade and a yawed rotor plane seen from the top For pure yaw the relative rotational speed varies sinusoidally as ωr(1+a’) – usinθ cosθ . Therefore both the relative velocity seen by the blade and the yaw wing angle of attack vary with the frequency of the rotor, ω. The loads therefore also vary in yaw; this will therefore contribute to the fatigue loads and thus influence the expected lifetime of the rotor. It should now be clear that the rotor blades experience a variation of the angle of attack due to turbulence, wind shear, tower passage and yaw/tilt. The forces and moments on the blades and thus on the entire structure will therefore also vary in time. Since a wind turbine is designed to last at least 20 years, it is very important to quantify the loads in order to be able to perform a reliable fatigue analysis. To make a wind turbine last for the design period, one must also take situations like extreme wind speeds into account. In this case the blades are parked or idling and the Danish standard (DS 412, 1992), for example, describes how to compute the extreme loads as: p(r) = q C c(r), (13.4) 2s f 1 2 where C is a force coefficient, c(r) the chord and q = –ρV the dynamic f 2s 2 2s pressure from an extreme wind speed time averaged over 2 seconds. For a Danish homogeneous terrain, V can be computed using: 2s h V = Vk (1n(–) + 3) = V + 3, (13.5) 2s b t 10min z o where V = 27m/s is a basis wind speed, h the height above the ground (mini- b mum 4m), z the roughness length and k a terrain factor.  is the standard o t3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 145 Sources of Loads on a Wind Turbine 145 deviation in a 10 minute time series. The terrain factor k is related to the t roughness length z as: o z m k o t 0.004 0.16 0.010 0.17 0.050 0.19 0.300 0.22 The extreme wind speed time, averaged over 10 minutes, can be estimated from: h V = Vk 1n( —). (13.6) 10min b t z o The following simplified example shows how to quantify extreme loads on wind turbine blades. In the Danish standard (DS 472, 1992) for loads and safety in wind turbine construction, it is stated that for the blades the extreme wind speed V must be calculated at h = h +2/3R using C = 1.5. For a 2s hub f wind turbine with a hub height of 40m and a rotor radius of 20m, this corresponds to h = 40+2/3·20 = 53.3m. If the surrounding landscape has no nearby obstacles, such as houses, and very low vegetation, the roughness length is approximately 0.01m and the terrain factor k = 0.17. In this case V t 2s = 27·0.17(ln(53.3/0.01)+3) = 53.2m/s. Assuming that the chord, c, is a 3 constant 1.3m and the density is 1.28kg/m , the load according to equation 1 2 (13.4) is p = – .1.28.53.2 .1.3.1.5 = 3532N/m. The root bending moment at 2 r = 3m for the constant load then becomes: R 1 2 2 1 2 2 M = ∫ rp(r)dr = – p(R –r ) = – 3532(20 –3 ) = 690533Nm. (13.7) r 2 2 For a simplified structural cross section such as the one shown in Figure 13.6, the moment of inertia about the flapwise axis is: 1 3 3 2 3 3 I = — a((2b ) – (2b ) ) = –a(b – b ). (13.8) l 2 l 2 12 3 9 2 Further it is assumed that E = 14Gpa = 14 10 N/m and that the thickness of the aerofoil is t/c = 35%. The thickness is thus t =(t/c)·c = 0.35·1.3m = 0.455m and b = t/2 = 0.228m. 1 Given the extreme wind speed, the thickness of the aerofoil and E, it is possible to estimate the necessary shell thickness t = b –b so that the blade s 1 2 does not break.3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 146 146 Aerodynamics of Wind Turbines Figure 13.6 Simplified structural model Since the force in the tangential direction is small, equation (11.10) reduces to: M ε = –– y (13.9) EI and it is clear that the maximum strain occurs for y = b . From testing it has 1 been found that the material used in this example fails for ε = ε = 0.02. The fail necessary moment of inertia can now be estimated from equation (13.9) as: 690533 M –4 4 I = ––– b = –––––––––– 0.228 = 5.62·10 m . (13.10) 1 Eε 9 14·10 ·0.02 fail Using this necessary moment of inertia in equation (13.8), b can be found 2 for a = 0.5m: –4 3·5·62·10 3 3I 1 3 1 b = (b – — ) /3 = (0.228 – ––––––––––) /3 = 0.217m. (13.11) 2 1 2a 2·0.5 The necessary shell thickness is thus t = b Nb = 0.228-0.217 = 0.011m = s 1 2 1.1cm. In the standard (DS 472, 1992) it is further stated, however, that the above sketched method can only be applied if the construction is assumed not to auto-vibrate. References DS 472 (1992) Dansk Ingeniørforenings og ingeniørsammenslutningens Norm for ‘Last og sikkerhed for vindmøllekonstruktioner’ (in Danish), Danish standard Troen, I. and Petersen, E. L. (1989) European Wind Atlas, Risø National Laboratory3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 147 14 Wind Simulation To calculate realistic time series for the loads on a wind turbine construction as shown in Figure 10.2, for example, it is important to have as input a wind field with correct spatial and temporal variation, as sketched in Figure 13.4. This variation should include turbulence, wind shear and the effect on the wind from the tower. This chapter describes some of the statistical properties of atmospheric turbulence and how to model a 3-D wind field. Wind Simulation at One Point in Space A simple anemometer measures the wind speed at one point with a sample frequency of f =1/∆ t, where ∆ t is the time between two measurements. The s output is a list of numbers, ui i = 1,..,N, where the corresponding time is t = 1·∆ t,2·∆ t,...,N·∆ t The total time is T = ∆ t·N and the sample frequency is f = s N/T. An example of such a time history is shown in Figure 14.1. Figure 14.1 Time history of discrete sampled wind speed at one point3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 148 148 Aerodynamics of Wind Turbines Figure 14.1 shows that the highest frequency that can be resolved is f = f /2 h s = (N/2)/T, since three discrete points are needed as a minimum to describe a vibration. The lowest frequency that can be resolved is f = 1/T. If the signal low is assumed to be periodical, the time history can be decomposed using a discrete Fourier transformation as: N/2 2n u(t) = a + a cos(ω t) + b sin(ω t), ω = –––– . (14.1) o n n n n n T n=1 The coefficients are found as: N 1 a = — u (mean wind speed) (14.2) o i N n=1 N 2n 2 N a = — u cos( –––– i), n = 1, ..., — –1 (14.3) n i N N 2 i=1 N 2n 2 N b = – u sin( –––– i), n = 1, ..., — –1 (14.4) n i N N 2 i=1 N 1 a = – u cos(i) (14.5) N i — N 2 i=1 b = 0 (14.6) N — 2 Putting equation (14.1) into the definition of the variance, , the following relationship can be proved: N N/2–1 _ 2 1 2 1 2 2 2  = – (u –u) = – (a + b )+a . (14.7) i n n N — 2 N i=1 n=1 2 For the power spectral density function, PSD( f ), the variance is: ∞ 2  = ∫ PSD( f )df, (14.8) 0 where f is the frequency. Discretizing the integral (14.8), and given that the frequencies lie between f =1/T and f =(N/2)/T, yields (see Figure 14.2): 1ow N/2 N/2 ω 2 n n   PSD( f )∆ f, f = –– = – . (14.9) n n 2π T n=13212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 149 Wind Simulation 149 Comparing equations (14.9) and (14.7), it is seen that the power spectral density function PSD( f ) can be evaluated as: 1 2 2 T 2 2 PSD( f ) = —– (a + b ) = – (a + b ). (14.10) n n n n n 2∆ f 2 Figure 14.2 Power spectral density function From a measured time series, u ,...,u , the Fourier coefficients and thus the 1 N power spectral density function can be determined from equations (14.3), (14.4) and (14.10). This method is called discrete Fourier transform (DFT). Constructing a time series from a known power spectral density (PSD) function is known as the inverse DFT. For the atmospheric boundary layer, different analytical expressions to approximate measured PSD functions exist, for example a Kaimal spectrum as given in DS 472 (1992): 2 I V l 10min PSD( f ) = ––––––––––– (14.11) f·l 5 / (1+1.5 –––– ) 3 V 10mins I = σ/V is the turbulence intensity, f is the frequency (in Hz), V is the 10 10min 10min minute averaged wind speed, and l is a length scale, l = 20h for h 30m and3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 150 150 Aerodynamics of Wind Turbines l = 600m for h 30m, where h is the height above ground level. The inverse DFT should provide a Fourier transform, as shown below, that satisfies the prescribed PSD function: N/2 – u(t)= u + a cos(ω t)+b sin(ω t). (14.12) n n n n n=1 Equation (14.12) can be rewritten as: N/2 –––––– – 2 2 u(t) = u +  a +b cos(ω t–ϕ ), (14.13) n n n n n=1 – ––––– 2 2 where ϕ is the phase angle at frequency ω . Replacing the term a +b by n n n n 2PSD( f )/T in equation (14.13) using, for example, equation (14.11) for the n spectrum, yields a Fourier transform that exactly satisfies the prescribed PSD function: –––––––– N/2 2PSD(ω ) n – u(t)= u + –––––––– – cos(ω t–ϕ ). (14.14) n n  T n=1 t = i·∆ t for i = 1,...,N The phase angle ϕ is not reflected in the PSD function and can be modelled n using a random number generator yielding a value between 0 and 2π. Using equation (14.14) it is very easy to compute a discrete time series having exactly the prescribed PSD function. A note should be made on the PSD function such as, for example, equation (14.11). These functions assume that all frequencies between 0 and ∞ are present, but as described earlier only frequencies between 1/T and (N/2)/T are present for the discrete time series. It is therefore practical to scale the PSD function as: f =(N/2)/T PSD( f )df = 1, (14.15) ∫ f =1/T corresponding to a standard deviation of 1, in other words a turbulent inten- sity of I = 1/V . A standard deviation of S and thus a turbulent intensity of 10min I = S/V can be found simply as: 10min – – u – u = S(u – u) (14.16) i i or: – u = Su +(1– S)u (14.17) i i3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 151 Wind Simulation 151 More information on different spectras such as the Kaimal and von Karman can be found in Burton et al (2001) and Rohatgi and Nelson (1994). 3-D Wind Simulation To simulate a time history of the wind speed at two or more points in space, it must also be considered that the time histories are not independent. This dependency is of course related to the physical distance between two points, but also to the frequency. The high frequency content of the time series is a result of small vortices, which have small spatial influence. Similarly the low frequency part is related to large-scale vortices covering a bigger volume of the flow. A necessary coherence function must therefore take into account both the distance, L, between points j and k and the frequency; one possible choice is given in DS 472 (1992): coh (L, f ) = exp(–12(fL/V )). (14.18) jk 10min Veers (1988) presents a method that will generate a 3-D wind field with a prescribed PSD function and coherence function. In the following the method will not be proved but only given as an algorithm. First, a matrix, S , is created: jk –––––– Sjk = coh S ·S . (14.19) jk jj kk where S and S are the PSD functions of points j and k respectively. The off- jj kk diagonal terms are the magnitudes of the cross-spectras. If the number of points in space is NP, S is an NP NP matrix. jk Second, a lower triangular H matrix is computed through following recur- sive formulae: 1/2 H = S 11 11 H = S /H 21 21 11 2 1/2 H = (S – H ) 22 22 21 H = S /H 31 31 11 . (14.20) . . k–1 H = (S – H H )/H jk jk jl kl kk l=l k–1 1/2 2 H = (S – H ) kk kk kl l=l3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 152 152 Aerodynamics of Wind Turbines For each point indexed by k and for each discrete frequency, f = m/T, a m random number, ϕ , between 0 and 2π is found to represent the phase as in km equation (14.14). m varies between 1 and N/2, where N is the number of discrete points in the time histories (t = i⋅∆ t,i = 1,…,N). A vector with imaginary components of length equal to the number of points in space, V = V( f ), is now calculated as: j m j Re(V ( f )) = H cos(ϕ ) j m jk km k=l (14.21) j Im(V ( f ))= H sin(ϕ ) j m jk km k=l Re(V ( f )) and Im(V ( f )) is transformed to an amplitude, Amp ( f ), and a j m j m j m phase, Φ ( f ), as: j m 2 2 Amp ( f )= Re(V ( f )) + Im(V ( f )) (14.22) j j m m j m Im(V ( f )) j __ _ _______ __ _ m tan Φ ( f )= j m Re(V ( f )) j m Finally, the time histories at the points j =1,NP can be computed as: N/2 – U(t)= U + 2Amp ( f )cos(2πf ·t – Φ ( f )) (14.23) j j m j m m m =1 t = i·∆ t for i = l,...,N Figure 14.3 plots the time series of two points spaced 1m apart. The mean wind speed is 10m/s and the turbulence intensity is 0.1. It is clearly seen that the two curves are well correlated for the lower frequencies. Figure 14.4 plots the coherence function computed from the two time series shown in Figure 14.3 together with the input given by equation (14.18). Using equation (14.23), together with appropriate PSD and coherence functions, the time histories at all points are computed for each velocity component U = (u,v,w) independently. Therefore there is no guarantee of obtaining correct cross-correlations. In Mann (1998) a method ensuring this is developed on the basis of the linearized Navier-Stokes equations. In the future, wind fields are expected to be generated numerically from large eddy simulations (LES) or direct numerical simulations (DNS) of the Navier- Stokes equations for the flow on a landscape similar to the actual siting of a specific wind turbine. For an aeroelastic calculation of a wind turbine construction it is natural to construct the time histories in a series of points arranged as indicated in3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 153 Wind Simulation 153 The mean wind speed is 10m/s and the turbulence intensity is 0.1. Figure 14.3 Computed time series of wind speed in two points separated by 1m Figure 14.4 Comparison of actual coherence from the two time series shown in Figure 14.3 and specified by equation (14.18)3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 154 154 Aerodynamics of Wind Turbines Figure 14.5. The velocities on the blade sweeping through the grid must in general be found by spatial interpolation. It should be mentioned that the time history of the wind seen by a point on the blade is different from the time history of a point fixed in space. A time history for a point on the rotating blade is called rotational sampling; Veers (1988) shows how this can be directly calculated for a blade with a constant rotational speed. Figure 14.5 Sketch of point distribution for time histories of wind speed for aeroelastic calculation of wind turbine construction To illustrate rotational sampling a number of time series are created for a constant radius of r = 5m and azimuthally spaced by 7.5°. The velocity seen by a blade at the same radius and rotating with ω = 13.1 rad/s is found by interpolation in the time series. In Figure 14.6 the calculated PSD from the rotational sampling is plotted together with the PSD from a stationary point3212 J&J Aerodynamic Turbines 15/11/07 1:43 PM Page 155 Wind Simulation 155 Figure 14.6 Comparison between specified PSD and the PSD seen by rotating blade in space. The PSD for the rotating blade shows clear spikes at multiples of the rotational frequency 1P, 2P,… (the frequency in Figure 14.6 is non- dimensionalized with the rotational frequency of the rotor, f = ω/2π). These 0 spikes, which will make a major contribution to the fatigue damage of the wind turbine construction, are due to the spatial coherence, since for a frozen time there exist areas of the rotor disc where the velocity is relatively high or low (see, for example Figure 13.4), and a blade will pass these areas once per revolution. Finally, it should be noted that the simulated time series not only depend on the prescribed PSD and coherence functions, but also on the random numbers for the phases. Therefore it is recommended to calculate at least three different time series for each case and perform for each series a run with the aeroelastic code to get an idea of the uncertainty from the different time series. References Burton, T., Sharpe, D., Jenkins, N. and Bossanyi, E. (2001) Wind Energy Handbook, John Wiley & Sons, Chichester, UK DS 472 (1992) Dansk Ingeniørforenings og ingeniørsammenslutningens Norm for ‘Last og sikkerhed for vindmøllekonstruktioner’ (in Danish), Danish standard
Website URL
Comment