GH Bladed Theory Manual

GH Bladed Theory Manual Document No Classification Issue no. Date

282/BR/009 Commercial in Confidence 11 July 2003

Author: E A Bossanyi Checked by: D C Quarton Approved by: D C Quarton

DISCLAIMER Acceptance of this document by the client is on the basis that Garrad Hassan and Partners Limited are not in any way to be held responsible for the application or use made of the findings of the results from the analysis and that such responsibility remains with the client.

Key To Document Classification Strictly Confidential

:

Recipients only

Private and Confidential

:

For disclosure to individuals directly concerned within the recipient’s organisation

Commercial in Confidence

:

Not to be disclosed outside the recipient’s organisation

GHP only

:

Not to be disclosed to non GHP staff

Client’s Discretion

:

Distribution at the discretion of the client subject to contractual agreement

Published

:

Available to the general public

© 2003 Garrad Hassan and Partners Limited

Garrad Hassan and Partners Ltd

Document: 282/BR/009

ISSUE:011

FINAL

CONTENTS

1.

Introduction 1.1 1.2 1.3 1.4 1.5

2.

1

Purpose Theoretical background Support Documentation Acknowledgements

1 2 3 3 3

AERODYNAMICS

4

2.1

Combined blade element and momentum theory 2.1.1 Actuator disk model 2.1.2 Wake rotation 2.1.3 Blade element theory 2.1.4 Tip and hub loss models 2.2 Wake models 2.2.1 Equilibrium wake 2.2.2 Frozen wake 2.2.3 Dynamic wake 2.3 Steady stall 2.4 Dynamic stall

3.

STRUCTURAL DYNAMICS

13

3.1

Modal analysis 3.1.1 Rotor modes 3.1.2 Tower modes 3.2 Equations of motion 3.2.1 Degrees of freedom 3.2.2 Formulation of equations of motion 3.2.3 Solution of the equations of motion 3.3 Calculation of structural loads

4. 4.1

4.2

4.3 4.4 4.5

5. 5.1 5.2

4 4 5 6 8 9 9 9 9 11 11 13 14 15 16 16 16 17 18

POWER TRAIN DYNAMICS

19

Drive train models 4.1.1 Locked speed model 4.1.2 Rigid shaft model 4.1.3 Flexible shaft model Generator models 4.2.1 Fixed speed induction generator 4.2.2 Fixed speed induction generator: electrical model 4.2.3 Variable speed generator 4.2.4 Variable slip generator Drive train mounting Energy losses The electrical network

19 19 19 19 20 20 21 22 23 24 24 25

CLOSED LOOP CONTROL

27

Introduction The fixed speed pitch regulated controller

-i-

27 27



5.2.1 Steady state parameters 5.2.2 Dynamic parameters 5.3 The variable speed stall regulated controller 5.3.1 Steady state parameters 5.3.2 Dynamic parameters 5.4 The variable speed pitch regulated controller 5.4.1 Steady state parameters 5.4.2 Dynamic parameters 5.5 Transducer models 5.6 Modelling the pitch actuator 5.7 The PI control algorithm 5.7.1 Gain scheduling 5.8 Control mode changes 5.9 Client-specific controllers 5.10 Signal noise and discretisation

6.

SUPERVISORY CONTROL

ISSUE:011

FINAL

28 28 28 28 30 31 31 32 33 33 36 37 38 38 39

40

6.1 6.2 6.3 6.4 6.5 6.6

Start-up Normal stops Emergency stops Brake dynamics Idling and parked simulations Yaw control 6.6.1 Active yaw 6.6.2 Yaw dynamics 6.7 Teeter restraint

40 41 41 42 42 42 42 43 44

MODELLING THE WIND

45

Wind shear 7.1.1 Exponential model 7.1.2 Logarithmic model Tower shadow 7.2.1 Potential flow model 7.2.2 Empirical model 7.2.3 Combined model Upwind turbine wake 7.3.1 Eddy viscosity model of the upwind turbine wake 7.3.2 Turbulence in the wake Time varying wind 7.4.1 Single point time history 7.4.2 3D turbulent wind 7.4.3 IEC transients Three dimensional turbulence model 7.5.1 The basic von Karman model 7.5.2 The improved von Karman model 7.5.3 The Kaimal model 7.5.4 Compatibility with IEC 1400-1 7.5.5 Using 3d turbulent wind fields in simulations

46 46 46 46 46 47 47 47 48 50 51 51 51 52 53 53 55 59 59 59

MODELLING WAVES AND CURRENTS

61

7. 7.1 7.2

7.3 7.4

7.5

8. 8.1 8.2

Tower and Foundation Model Wave Spectra

61 62 - ii -



ISSUE:011

8.2.1 JONSWAP / Pierson-Moskowitz Spectrum 8.2.2 User-defined Spectrum 8.3 Upper Frequency Limit 8.4 Wave Particle Kinematics 8.5 Wheeler Stretching 8.6 Simulation of Irregular Waves 8.7 Simulation of Regular Waves 8.8 Current Velocities 8.8.1 Near-Surface Current 8.8.2 Sub-Surface Current 8.8.3 Near-Shore Current 8.9 Total Velocities and Accelerations 8.10 Applied Forces 8.10.1 Relative Motion Form of Morison’s Equation 8.10.2 Longitudinal Pressure Forces on Cylindrical Elements

9.

POST-PROCESSING

62 62 63 63 64 64 66 67 68 68 68 69 69 69 69

71

9.1 9.2 9.3 9.4 9.5 9.6

Basic statistics Fourier harmonics, and periodic and stochastic components Extreme prediction Spectral analysis Probability, peak and level crossing analysis Rainflow cycle counting and fatigue analysis 9.6.1 Rainflow cycle counting 9.6.2 Fatigue analysis 9.7 Annual energy yield 9.8 Ultimate loads 9.9 Flicker

10.

FINAL

References

71 71 72 75 75 76 76 77 78 79 79

80

- iii -



ISSUE:011

FINAL

1. INTRODUCTION

1.1 Purpose GH Bladed is an integrated software package for wind turbine performance and loading calculations. It is intended for the following applications: • Preliminary wind turbine design • Detailed design and component specification • Certification of wind turbines With its sophisticated graphical user interface, it allows the user to carry out the following tasks in a straightforward way: • Specification of all wind turbine parameters, wind inputs and load cases. • Rapid calculation of steady-state performance characteristics, including: Aerodynamic information Performance coefficients Power curves Steady operating loads Steady parked loads • Dynamic simulations covering the following cases: Normal running Start-up Normal and emergency shut-downs Idling Parked Dynamic power curve • Post-processing of results to obtain: Basic statistics Periodic component analysis Probability density, peak value and level crossing analysis Spectral analysis Cross-spectrum, coherence and transfer function analysis Rainflow cycle counting and fatigue analysis Combinations of variables Annual energy yield Ultimate loads (identification of worst cases) Flicker severity • Presentation: results may be presented graphically and can be combined into a word processor compatible report.

1 of 82



ISSUE:011

FINAL

1.2 Theoretical background The Garrad Hassan approach to the calculation of wind turbine performance and loading has been developed over many years. The main aim of this development has been to produce reliable tools for use in the design and certification of wind turbines. The models and theoretical methods incorporated in GH Bladed have been extensively validated against monitored data from a wide range of turbines of many different sizes and configurations, including: • • • • • • • • • • • • • • • • • • • • • • • • • • • •

WEG MS-1, UK, 1991 Howden HWP300 and HWP330, USA, 1993 ECN 25m HAT, Netherlands, 1993 Newinco 500kW, Netherlands, 1993 Nordex 26m, Denmark, 1993 Nibe A, Denmark, 1993 Holec WPS30, Netherlands, 1993 Riva Calzoni M30, Italy, 1993 Nordtank 300kW, Denmark, 1994 WindMaster 750kW, Netherlands, 1994 Tjaereborg 2MW, Denmark, 1994 Zond Z-40, USA, 1994 Nordtank 500kW, UK, 1995 Vestas V27, Greece, 1995 Danwin 200kW, Sweden, 1995 Carter 300kW, UK, 1995 NedWind 50, 1MW, Netherlands, 1996 DESA, 300kW, Spain 1997 NTK 600, UK, 1998 West Medit, Italy, 1998 Nordex 1.3 MW, Germany, 1999 The Wind Turbine Company 350 kW, USA, 2000 Windtec 1.3 MW, Austria, 2000 WEG MS-4, 400 kW, UK, 2000 EHN 1.3 MW, Spain, 2001 Vestas 2MW, UK, 2001 Lagerwey 750 Netherlands, 2001 Vergnet 200, France 2001

This document describes the theoretical background to the various models and numerical methods incorporated in GH Bladed.

2 of 82



ISSUE:011

FINAL

1.3 Support GH Bladed is supplied with a one-year maintenance and support agreement, which can be renewed for further periods. This support includes a ‘hot-line’ help service by telephone, fax or e-mail: Telephone: Fax: E-mail

+44 (0)117 972 9900 +44 (0)117 972 9901 [email protected]

1.4 Documentation In addition to this Theory Manual, there is also a GH Bladed User Manual which explains how the code can be used.

1.5 Acknowledgements GH Bladed was developed with assistance from the Commission of the European Communities under the JOULE II programme, project no. JOU2-CT92-0198.

3 of 82



ISSUE:011

FINAL

2. AERODYNAMICS The modelling of rotor aerodynamics provided by Bladed is based on the well established treatment of combined blade element and momentum theory [2.1]. Two major extensions of this theory are provided as options in the code to deal with the unsteady nature of the aerodynamics. The first of these extensions allows a treatment of the dynamics of the wake and the second provides a representation of dynamic stall through the use of a stall hysteresis model. The theoretical background to the various aspects of the treatment of rotor aerodynamics provided by Bladed is given in the following sections.

2.1 Combined blade element and momentum theory At the core of the aerodynamic model provided by Bladed is combined blade element and momentum theory. The features of this treatment of rotor aerodynamics are described below. 2.1.1 Actuator disk model To aid the understanding of combined blade element and momentum theory it is useful initially to consider the rotor as an “actuator disk”. Although this model is very simple, it does provide valuable insight into the aerodynamics of the rotor. Wind turbines extract energy from the wind by producing a step change in static pressure across the rotor-swept surface. As the air approaches the rotor it slows down gradually, resulting in an increase in static pressure. The reduction in static pressure across the rotor disk results in the air behind it being at sub atmospheric pressure. As the air proceeds downstream the pressure climbs back to the atmospheric value resulting in a further slowing down of the wind. There is therefore a reduction in the kinetic energy in the wind, some of which is converted into useful energy by the turbine. In the actuator disk model of the process described above, the wind velocity at the rotor disk Ud is related to the upstream wind velocity Uo as follows: U d = ( 1 a )U o

The reduced wind velocity at the rotor disk is clearly determined by the magnitude of a, the axial flow induction factor or inflow factor. By applying Bernoulli’s equation and assuming the flow to be uniform and incompressible, it can be shown that the power P extracted by the rotor is given by : P = 2 AU o3a( 1 a )3

where

is the air density and A the area of the rotor disk.

The thrust T acting on the rotor disk can similarly be derived to give:

4 of 82



ISSUE:011

FINAL

T = 2 AU o2 a( 1 a )

The dimensionless power and thrust coefficients, CP and CT are respectively: CP = P / ( 1 2 AU o3 ) = 4a( 1 a )2

and: CT = T / ( 1 2 AU o2 ) = 4a( 1 a )

The maximum value of the power coefficient CP occurs when a is 1 /3 and is equal to 16/27 which is known as the Betz limit. The thrust coefficient CT has a maximum value of 1 when a is 1 /2. 2.1.2 Wake rotation The actuator disk concept used above allows an estimate of the energy extracted from the wind without considering that the power absorbed by the rotor is the product of torque Q and angular velocity of the rotor. The torque developed by the rotor must impart an equal and opposite rate of change of angular momentum to the wind and therefore induces a tangential velocity to the flow. The change in tangential velocity is expressed in terms of a tangential flow induction factor a’. Upstream of the rotor disk the tangential velocity is zero, at the disk the tangential velocity at radius r on the rotor is ra’ and far downstream the tangential velocity is 2 ra’. Because it is produced in reaction to the torque, the tangential velocity is opposed to the motion of the blades. The torque generated by the rotor is equal to the rate of change of angular momentum and can be derived as: Q=

R 4 (1 a )a ,U o

5 of 82



ISSUE:011

FINAL

2.1.3 Blade element theory Combined blade element and momentum theory is an extension of the actuator disk theory described above. The rotor blades are divided into a number of blade elements and the theory outlined above used not for the rotor disk as a whole but for a series of annuli swept out by each blade element and where each annulus is assumed to act in the same way as an independent actuator disk. At each radial position the rate of change of axial and angular momentum are equated with the thrust and torque produced by each blade element. The thrust dT developed by a blade element of length dr located at a radius r is given by: dT = 1 2 W 2 ( CL cos + CD sin )cdr

where W is the magnitude of the apparent wind speed vector at the blade element, is known as the inflow angle and defines the direction of the apparent wind speed vector relative to the plane of rotation of the blade, c is the chord of the blade element and CL and CD are the lift and drag coefficients respectively. The lift and drag coefficients are defined for an aerofoil by: CL = L / ( 1 2 V 2 S )

and CD = D / ( 1 2 V 2 S )

where L and D are the lift and drag forces, S is the planform area of the aerofoil and V is the wind velocity relative to the aerofoil. The torque dQ developed by a blade element of length dr located at a radius r is given by: dQ = 1 2 W 2 r( CL sin

CD cos )cdr

In order to solve for the axial and tangential flow induction factors appropriate to the radial position of a particular blade element, the thrust and torque developed by the element are equated to the rate of change of axial and angular momentum through the annulus swept out by the element. Using expressions for the axial and angular momentum similar to those derived for the actuator disk in Sections 2.1.1 and 2.1.2 above, the annular induction factors may be expressed as follows: a = g1 / ( 1 + g1 )

and a , = g2 / ( 1 g 2 )

where 6 of 82



g1 =

Bc ( CL cos + CD sin ) H 2 r 4 F sin 2

g2 =

Bc ( CL sin CD cos ) 2 r 4 F sin cos

ISSUE:011

FINAL

and

Here B is the number of blades and F is a factor to take account of tip and hub losses, refer Section 2.1.4. The parameter H is defined as follows: for a

0.3539, H = 10 .

for a > 0.3539, H =

4a (1 a ) (0.6 + 0.61a + 0.79a 2 )

In the situation where the axial induction factor a is greater than 0.5, the rotor is heavily loaded and operating in what is referred to as the “turbulent wake state”. Under these conditions the actuator disk theory presented in Section 2.1.1 is no longer valid and the expression derived for the thrust coefficient: CT = 4a( 1 a )

must be replaced by the empirical expression: CT = 0.6 + 0.61a + 0.79a 2

The implementation of blade element theory in Bladed is based on a transition to the empirical model for values of a greater than 0.3539 rather than 0.5. This strategy results in a smoother transition between the models of the two flow states. The equations presented above for a and a’ can only be solved iteratively. The procedure involves making an initial estimate of a and a’, calculating the parameters g1 and g2 as functions of a and a’, and then using the equations above to update the values of a and a’. This procedure continues until a and a’ have converged on a solution. In Bladed convergence is assumed to have occurred when: ak

ak

1

tol

a'k

a'k

1

tol

and

where tol is the value of aerodynamic tolerance specified by the user.

7 of 82



ISSUE:011

FINAL

2.1.4 Tip and hub loss models The wake of the wind turbine rotor is made up of helical sheets of vorticity trailed from each rotor blade. As a result the induced velocities at a fixed point on the rotor disk are not constant with time, but fluctuate between the passage of each blade. The greater the pitch of the helical sheets and the fewer the number of blades, the greater the amplitude of the variation of induced velocities. The overall effect is to reduce the net momentum change and so reduce the net power extracted. If the induction factor a is defined as being the value which applies at the instant a blade passes a given point on the disk, then the average induction factor at that point, over the course of one revolution will be aFt,, where Ft is a factor which is less than unity. The circulation at the blade tips is reduced to zero by the wake vorticity in the same manner as at the tips of an aircraft wing. At the tips, therefore the factor Ft becomes zero. Because of the analogy with the aircraft wing , where losses are caused by the vortices trailing from the tips, Ft is known as the tip loss factor. Prandtl [2.2] put forward a method to deal with this effect in propeller theory. Reasoning that, in the far wake, the helical vortex sheets could be replaced by solid disks, set at the same pitch as the normal spacing between successive turns of the sheets, moving downstream with the speed of the wake. The flow velocity outside of the wake is the free stream value and so is faster than that of the disks. At the edges of the disks the fast moving free stream flow weaves in and out between them and in doing so causes the mean axial velocity between the disks to be higher than that of the disks themselves, thus simulating the reduction in the change of momentum. The factor Ft can be expressed in closed solution form: Ft = 2 arccos[exp(

s )] d

where s is the distance of the radial station from the tip of the rotor blade and d is the distance between successive helical sheets. A similar loss takes place at the blade root where, as at the tip, the bound circulation must fall to zero and therefore a vortex must be trailed into the wake, A separate hub loss factor Fh is therefore calculated and the effective total loss factor at any station on the blade is then the product of the two: F = Ft Fh

The combined tip and hub loss factor is incorporated in the equations of blade element theory as indicated in Section 2.1.3 above.

8 of 82



ISSUE:011

FINAL

2.2 Wake models 2.2.1 Equilibrium wake The use of blade element theory for time domain dynamic simulations of wind turbine behaviour has traditionally been based on the assumption that the wake reacts instantaneously to changes in blade loading. This treatment, known as an equilibrium wake model, involves a re-calculation of the axial and tangential induction factors at each element of each rotor blade, and at each time step of a dynamic simulation. Based on this treatment the induced velocities along each blade are computed as instantaneous solutions to the particular flow conditions and loading experienced by each element of each blade. Clearly in this interpretation of blade element theory the axial and tangential induced velocities at a particular blade element vary with time and are not constant within the annulus swept out by the element. The equilibrium wake treatment of blade element theory is the most computationally demanding of the three treatments described here. 2.2.2 Frozen wake In the frozen wake model, the axial and tangential induced velocities are computed using blade element theory for a uniform wind field at the mean hub height wind speed of the simulated wind conditions. The induced velocities, computed according to the mean, uniform flow conditions, are then assumed to be fixed, or “frozen” in time. The induced velocities vary from one element to the next along the blade but are constant within the annulus swept out by the element. As a consequence each blade experiences the same radial distribution of induced flow.. It is important to note that it is the axial and tangential induced velocities aUo and a’r not the induction factors a and a’ which are frozen in time.

and

2.2.3 Dynamic wake As described above, the equilibrium wake model assumes that the wake and therefore the induced velocity flow field react instantaneously to changes in blade loading. On the other hand, the frozen wake model assumes that induced flow field is completely independent of changes in incident wind conditions and blade loading. In reality neither of these treatments is strictly correct. Changes in blade loading change the vorticity that is trailed into the rotor wake and the full effect of these changes takes a finite time to change the induced flow field. The dynamics associated with this process is commonly referred to as “dynamic inflow”. The study of dynamic inflow was initiated nearly 40 years ago in the context of helicopter aerodynamics. In brief, the theory provides a means of describing the dynamic dependence of the induced flow field at the rotor upon the loading that it experiences. The dynamic inflow model used within Bladed is based on the work of Pitt and Peters [2.3] which has received substantial validation in the helicopter field, see for example Gaonkar et al [2.4]. The Pitt and Peters model was originally developed for an actuator disk with assumptions made concerning the distribution of inflow across the disc. In Bladed the model is applied at blade element or actuator annuli level since this avoids any assumptions about the distribution of inflow across the disc. 9 of 82



ISSUE:011

FINAL

For a blade element, bounded by radii R1 and R2 , and subject to uniform axial flow at a wind speed Uo, the elemental thrust, dT, can be expressed as: dT = 2U o am + U o m A a&

where m is the mass flow through the annulus, mA is the apparent mass acted upon by the annulus and a is the axial induction factor. The mass flow through the annular element is given by: m = U o (1 a )dA

where dA is the cross-sectional area of the annulus. For a disc of radius R the apparent mass upon which it acts is given approximately by potential theory, Tuckerman, [2.5]: mA = 8

3

R3

Therefore the thrust coefficient associated with the annulus can be derived to give: C T = 4a (1 a ) +

16 (R 32 3 U o (R 22

R 13 ) R 12 )

a&

This differential equation can therefore be used to replace the blade element and momentum theory equation for the calculation of axial inflow. The equation is integrated at each time step to give time dependent values of inflow for each blade element on each blade. The tangential inflow is obtained in the usual manner and so depends on the time dependent axial value. It is evident that the equation introduces a time lag into the calculation of inflow which is dependent on the radial station. It is probable that the values of time lag for each blade element calculated in this manner will under-estimate somewhat the effects of dynamic inflow, as each element is treated independently with no consideration of the three dimensional nature of the wake or the possibly dominant effect of the tip vortex. The treatment is, however, consistent with blade element theory and provides a simple, computationally inexpensive and reasonably reliable method of modelling the dynamics of the rotor wake and induced velocity flow field.

10 of 82



ISSUE:011

FINAL

2.3 Steady stall The representation and to some extent the general understanding of aerodynamic stall on a rotating wind turbine blade remain rather poor. This is a rather extraordinary situation in view of the importance of stall regulation to the industry. Stall delay on the inboard sections of rotor blades, due to the three dimensionality of the incident flow field, has been widely confirmed by measurements at both model and full scale. A number of semi-empirical models [2.6, 2.7] have been developed for correcting two dimensional aerofoil data to account for stall delay. Although such models are used for the design analysis of stall regulated rotors, their general validity for use with a wide range of aerofoil sections and rotor configurations remains, at present, rather poor. As a consequence Bladed does not incorporate models for the modification of aerofoil data to deal with stall delay, but the user is clearly able to apply whatever correction of the aerofoil data he believes is appropriate prior to its input to the code.

2.4 Dynamic stall Stall and its consequences are fundamentally important to the design and operation of most aerodynamic devices. Most conventional aeronautical applications avoid stall by operating well below the static stall angle of any aerofoils used. Helicopters and stall regulated wind turbines do however operate in regimes where at least part of their rotor blades are in stall. Indeed stall regulated wind turbines rely on the stalling behaviour of aerofoils to limit maximum power output from the rotor in high winds. A certain degree of unsteadiness always accompanies the turbulent flow over an aerofoil at high angles of attack. The stall of a lifting surface undergoing unsteady motion is more complex than static stall. On an oscillating aerofoil, where the incidence is increasing rapidly, the onset of the stall can be delayed to an incidence considerably in excess of the static stall angle. When dynamic stall does occur, however, it is usually more severe than static stall. The attendant aerodynamic forces and moments exhibit large hysteresis with respect to the instantaneous angle of attack, especially if the oscillation is about a mean angle close to the static stall angle. This represents an important contrast to the quasi-steady case, for which the flow field adjusts immediately, and uniquely, to each change in incidence. Many methods of predicting the dynamic stall of aerofoil sections have been developed, principally for use in the helicopter industry. The model adopted for inclusion of unsteady behaviour of aerofoils is that due to Beddoes [2.8]. The Beddoes model was developed for use in helicopter rotor performance calculations and has been formulated over a number of years with particular reference to dynamic wind tunnel testing of aerofoil sections used on helicopter rotors. It has been used successfully by Harris [2.9] and Galbraith et al [2.10] in the prediction of the behaviour of vertical axis wind turbines.

11 of 82



ISSUE:011

FINAL

The model used within Bladed is a development of the Beddoes model which has been validated against measurements from several stall regulated wind turbines. The model utilises the following elements of the method described in [2.8] to calculate the unsteady lift coefficient • The indicial response functions for modelling of attached flow • The time lagged Kirchoff formulation for the modelling of trailing edge separation and vortex lift The use of the model of leading edge separation has been found to be inappropriate for use on horizontal axis wind turbines where the aerofoil characteristics are dominated by progressive trailing edge stall. The time lag in the development of trailing edge separation is a user defined parameter within the model implemented in Bladed. This time lag encompasses the delay in the response of the pressure distribution and boundary layer to the time varying angle of attack. The magnitude of the time lag is directly related to the level of hysteresis in the lift coefficient. The drag and pitching moment coefficients are calculated using the quasi-steady input data along with the effective unsteady angle of attack determined during the calculation of the lift coefficient.

12 of 82



ISSUE:011

FINAL

3. STRUCTURAL DYNAMICS In the early days of the industry, wind turbine design was undertaken on the basis of quasistatic aerodynamic calculations with the effects of structural dynamics either ignored completely or included through the use of estimated dynamic magnification factors. From the late 1970’s research workers began to consider more reliable methods of dynamic analysis and two basic approaches were considered: finite element representations and modal analysis. The traditional use of standard, commercial finite element analysis codes for dealing with problems of structural dynamics is problematic in the case of wind turbines. This is because of the gross movement of one component of the structure, the rotor, with respect to another, the tower. Standard finite element packages are only used to consider structures in which motion occurs about a mean undisplaced position and for this reason the finite element models of wind turbines which have been developed have been specially constructed to deal with the problem. The form of wind turbine dynamic modelling most commonly used as the basis of design calculations is that involving a modal representation. This approach, borrowed from the helicopter industry, has the major advantage that it offers a reliable representation of the dynamics of a wind turbine with relatively few degrees of freedom. The number and type of modal degrees of freedom used to represent the dynamics of a particular wind turbine will clearly depend on the configuration and structural properties of the machine. At present, largely because of the very extensive computer processing requirements associated with the use of finite element models, the state of the art in the context of wind turbine dynamic modelling for design analysis is based squarely on the use of limited degree of freedom modal models. The representation of wind turbine structural dynamics within Bladed is based on a modal model.

3.1 Modal analysis Because of the rotation of the blades of a wind turbine relative to the tower support structure, the equations of motion which describe its dynamics contain terms with periodic coefficients. This periodicity means that the computation of the modal properties of an operating wind turbine as a complete structural entity is not possible using the standard eigen-analysis offered by commercial finite element codes. One solution to this problem is to make use of Floquet analysis to determine the modal properties of the periodic system. However, the mode shapes obtained by such calculations are complex and not directly useful for a forced response analysis. An alternative solution is based on the use of “component mode synthesis”. Here the modal properties of the rotating and non-rotating components of the wind turbine are computed independently. The component modes are then coupled by an appropriate formulation of the equations of motion of the wind turbine in the forced response analysis. This approach has been adopted for Bladed.

13 of 82



ISSUE:011

FINAL

3.1.1 Rotor modes The vibration of the tapered and twisted blades of a wind turbine rotor is a complex phenomenon. A classical method of representing the vibration is by means of the orthogonal, uncoupled “normal” modes of the structure. Each mode is defined in terms of the following parameters: • Modal frequency,

i

• Modal damping coefficient, • Mode shape,

i

i

(r )

where the subscript i indicates properties related to the ith mode. The modal frequencies and mode shapes of the rotor are calculated based on the following information: The mass distribution along the blade. The mass distribution is defined as the local mass density (kg/m) at each radial station in addition to the magnitude and location of any discrete, lumped masses. The bending stiffnesses along the blade. The bending stiffnesses are defined in local flapwise and edgewise directions at each radial station. The twist angle distribution along the blade. The mode shapes are computed in the rotor in-plane and out-of-plane directions and hence the flapwise and edgewise stiffnesses at each radial station are resolved through the local twist angle. The blade pitch and setting angles. The mode shapes are computed in the rotor in-plane and out-of-plane directions and hence the flapwise and edgewise stiffnesses at each radial station are resolved through the blade pitch and setting angles. The user of Bladed may select a series of different pitch angles for which the modal analysis is carried out. During subsequent dynamic simulations, the modal frequencies appropriate to the instantaneous blade pitch angle are therefore obtained by linear interpolation of the results of the modal analyses. The presence or otherwise of a hub teeter hinge for a two bladed rotor. For a two-bladed rotor the hub can be rigid or teetered. The presence of a teeter hinge will introduce asymmetric rotor modes involving out-of-plane rotation of the rotor about the teeter hinge. The presence or otherwise of a flap hinge for a one-bladed rotor. For a one-bladed rotor the hub can be rigid or have a flap hinge. The presence of a flap hinge will introduce rotor modes involving out-of-plane rotation of the rotor about the teeter hinge. The counter-weight mass and moment of inertia about the flap hinge for a one-bladed rotor. Whether the hub can rotate. 14 of 82



ISSUE:011

FINAL

Rotation of the hub will affect the frequencies and mode shapes of the in-plane rotor modes. With the shaft brake engaged and the rotor locked in position, the in-plane modes will include both symmetric and asymmetric cantilever-type modes. With the rotor free to rotate, the cantilever-type asymmetric modes will be replaced by asymmetric modes involving rotation about the rotor shaft. The rotational speed of the rotor. The frequencies and mode shapes of both in-plane and out-of-plane modes will be dependent on the rotational speed of the rotor. This dependence is explained by the additional bending stiffness developed because of centrifugal loads acting on the deflected rotor blades. The user of Bladed may select different rotational speeds for which the modal analysis is carried out. During subsequent dynamic simulations, the modal frequencies appropriate to the instantaneous rotational speed are therefore obtained by quadratic interpolation of the results of the modal analyses. The frequencies and mode shapes of the rotor modes are computed from the eigen-values and eigen-vectors of a finite element representation of the rotor structure. The finite element model of the rotor is based on the use of two-dimensional beam elements to describe the mass and stiffness properties of the rotor blades. The outputs from the modal analysis of the rotor are the modal frequencies and mode shapes defined in the rotor in-plane and out-of-plane directions. The modal damping coefficients are an input defined by the user and may be used to represent structural damping. 3.1.2 Tower modes The representation of the bending dynamics of the tower is based on the modal degrees of freedom in the fore-aft and side-side directions of motion. As for the rotor, the tower modes are defined in terms of their modal frequency, modal damping and mode shape. The modal frequencies and mode shapes of the tower are calculated based on the following information: The mass distribution along the tower. The mass distribution is defined as the local mass density (kg/m) at each tower station height in addition to the magnitude and location of any discrete, lumped masses. The bending stiffness along the tower. The tower is assumed to be axisymmetric with the bending stiffness therefore independent of bending direction. The mass, inertia and stiffness properties of the tower foundation. The influence of the foundation mass and stiffness properties on the tower bending modes may be taken into account. The model takes account of motion of the foundation mass and inertia against both translational and rotational stiffnesses. The mass and inertia of the nacelle and rotor For calculation of the tower modes, the nacelle and rotor are modelled as lumped mass and inertia located at the nacelle centre of gravity and rotor hub respectively. For one and twobladed rotors, the influence of the rotor inertia on the tower modal characteristics depends on the rotor azimuth and this may therefore be defined by the user. The variation of the tower modal frequencies with rotor azimuth is normally small and the assumption of a single rotor

15 of 82



ISSUE:011

FINAL

azimuthal position for the modal analysis is therefore a reasonable approximation. The user can, of course, determine the extent of the azimuthal variation in the tower modal frequencies by undertaking the modal analysis at a series of different rotor azimuths. The frequencies and mode shapes of the tower modes are computed from the eigen-values and eigen-vectors of a finite element representation of the tower structure. The finite element model of the tower is based on the use of two-dimensional beam elements to describe the mass and stiffness properties of the tower. The outputs from the modal analysis of the tower are the modal frequencies and mode shapes defined in the fore-aft and side-side directions. The modal damping coefficients are an input defined by the user and may be used to represent structural damping.

3.2 Equations of motion Because of the complexity of the coupling of the modal degrees of freedom of the rotating and non-rotating components, the algebraic manipulation involved in the derivation of the equations of motion for a wind turbine is a complicated problem. In the case of the dynamic model within Bladed, the derivation has been carried out using energy principles and Lagrange equations by means of a computer algebra package. 3.2.1 Degrees of freedom The degrees of freedom involved in the equations of motion for the structural dynamic model for Bladed are as follows: • • • • •

Rotor out of plane including teeter, maximum six modes Rotor in-plane, maximum six modes Nacelle yaw Tower fore-aft, maximum three modes Tower side-side, maximum three modes

In addition, a sophisticated representation of the power train dynamics is offered as described in Section 4 of this manual. 3.2.2 Formulation of equations of motion The equation of motion for a single modal degree of freedom, assuming no coupling with other degrees of freedom, is as follows:

q&&i + 2

i

q& +

i i

2 i

q = Fi / Mi

where: qi is the time dependent modal displacement,

Mi =

m(r )

2 i

(r )dr is the modal mass,

rotor

16 of 82



ISSUE:011

FINAL

and:

Fi =

f (r ) i (r )dr is the modal force. rotor

Here f(r) is the distributed force over the rotor or tower component. The modal degrees of freedom are, of course, coupled and the formulation of the equations of motion within Bladed is as follows: && + [ C]q& + [ K ]q = F [ M ]q

where [M], [C] and [K] are the modal mass, damping and stiffness matrices, q is the vector of modal displacements and F the vector of modal forces. The system matrices are full due to the coupling of the degrees of freedom and contain periodic coefficients because of the time dependent interaction of the dynamics of the rotor and tower. Because of their complexity, the equations of motion are not presented in this manual. The following key comments are, however, provided: • Although the equations of motion are based on a linear modal treatment of the structural dynamics, the model does contain non-linear terms associated primarily with gyroscopic coupling. • The rotor teeter degree of freedom is provided through the first out-of-plane mode and the equation of motion includes representation of mechanical damping, stiffness and pre-load restraints as specified by the user. • The equation of motion for the nacelle yaw degree of freedom is based on the inertia of the wind turbine about the yaw axis with mechanical restraints provided through yaw damping and stiffness as specified by the user. • The aeroelasticity of the wind turbine is taken into account in the equations of motion by consideration of the interaction of the total structural velocity vector with the wind velocity vector at each element along the rotor blades. The total structural velocity vector at each element on the rotor blades is composed of the appropriate summation of the velocities associated with each structural degree of freedom. In addition to the feedback of the structural velocities into the rotor blade aerodynamics, the structural displacement associated with the rotor teeter and nacelle yaw is also taken into account. 3.2.3 Solution of the equations of motion The equations of motion are solved by time-marching integration of the differential equations using a variable step size, fourth order Runge Kutta integrator.

17 of 82



ISSUE:011

FINAL

3.3 Calculation of structural loads The structural loads acting on the rotor, power train and tower are computed by the appropriate summation of the applied aerodynamic loads and the inertial loads. The inertial loads are calculated by integration of the mass properties and the total acceleration vector at each station. The total acceleration vector includes modal, centrifugal, Coriolis and gravitational components.

18 of 82



ISSUE:011

FINAL

4. POWER TRAIN DYNAMICS The power train dynamics define the rotational degrees of freedom associated with the drive train, including drive train mountings, and the dynamics of the electrical generator. The drive train consists of a low speed shaft, gearbox and high speed shaft. Direct drive generators can also be modelled.

4.1 Drive train models 4.1.1 Locked speed model The simplest drive train model which is available is the locked speed model, which allows no degrees of freedom for the power train. The rotor is therefore assumed to rotate at an absolutely constant speed, and the aerodynamic torque is assumed to be exactly balanced by the generator reaction torque at every instant. Clearly this model is unsuitable for start-up and shut-down simulations, but it is useful for quick, preliminary calculations of loads and performance before the drive train and generator have been fully characterised. 4.1.2 Rigid shaft model The rigid shaft model is obtained by selecting the dynamic drive train model with no shaft torsional flexibility. It allows a single rotational degree of freedom for the rotor and generator. It can be used for all calculations and is recommended if the torsional stiffness of the drive train is high. The acceleration of the generator and rotor are calculated from the torque imbalance divided by the combined inertia of the rotor and generator, making allowance for the gearbox ratio. Direct drive generators are modelled simply by setting the gearbox ratio to 1. The torque imbalance is essentially the difference between the aerodynamic torque and the generator reaction torque and any applied brake torque, taking the gearbox ratio into account. However, this is corrected to account for the inertial effect of blade deflection due to any edgewise blade vibration modes. To use the rigid shaft model, a model of the generator must also be provided, so that the generator reaction torque is defined. During a parked simulation, or once the brake has brought the rotor to rest during a stopping simulation, the actual brake torque balances the aerodynamic torque exactly (making allowance for the gearbox ratio if the brake is on the high speed shaft) and there is no further rotation. However, if the aerodynamic torque increases to overcome the maximum or applied brake torque, the brake starts to slip and rotation recommences. The rigid drive train model may be used in combination with flexible drive train mountings. In this case the equations of motion are more complex - see Section 4.3. 4.1.3 Flexible shaft model The flexible shaft model is obtained by selecting the dynamic drive train model with torsional flexibility in one or both shafts. It allows separate degrees of freedom for the rotation of the turbine rotor and the generator rotor. The torsional flexibility of the low speed and high speed shafts may be specified independently. As with the rigid shaft model, a model of the generator must be provided so that the generator reaction torque is specified.

19 of 82



ISSUE:011

FINAL

The turbine rotor is accelerated by the torque imbalance between the aerodynamic torque (adjusted for the effect of edgewise modes as explained in Section 4.1.2) and the low speed shaft torque. The generator rotor is accelerated by the imbalance between high speed shaft torque and generator reaction torque. The shaft torques are calculated from the shaft twist, together with any applied brake torque contributions depending on the location of the brake, which may be specified as being at either end of either the low or high speed shaft. During a parked simulation, or once the brake disk has come to rest during a stopping simulation, the equations of motion change depending on the brake location. If the brake is immediately adjacent to the rotor or generator then there is no further rotation of that component, but the other component continues to move and oscillates against the torsional flexibility of the shafts. If the brake is adjacent to the gearbox and both shafts are flexible, then both rotor and generator will oscillate. However, if the torque at the brake disk increases to overcome the maximum or applied brake torque, then the brake starts to slip again. The flexible drive train model may be used in combination with flexible drive train mountings. In this case the equations of motion are more complex - see Section 4.3. It should be pointed out that while the flexible shaft model provides greater accuracy in the prediction of loads, there is potential for one of the drive drain vibrational modes to be of relatively high frequency, depending on the generator inertia and shaft stiffnesses. The presence of this high frequency mode could result in slower simulations.

4.2 Generator models The generator characteristics must be provided if either the rigid or flexible shaft drive train model is specified. Three generator models are available: • A directly-connected induction generator model (for constant speed turbines), • A variable speed generator model (for variable speed turbines), and • A variable slip generator model (providing limited range variable speed above rated) 4.2.1 Fixed speed induction generator This model represents an induction generator directly connected to the grid. Its characteristics are defined by the slip slope h and the short-circuit transient time constant . The air-gap or generator reaction torque Q is then defined by the following differential equation:

Q& = 1 [h( where

0

) Q]

is the actual generator speed and

0

is the generator synchronous or no-load speed.

The slip slope is calculated as

Pr

h= r

(

r

0

)

20 of 82



ISSUE:011

where r is the generator speed at rated power output Pr , given by r = S is the rated slip in %, and is the full load efficiency of the generator.

0

FINAL

(1 + S/100) where

4.2.2 Fixed speed induction generator: electrical model A more complete model of the directly-connected induction generator is also available in Bladed. This model requires the equivalent circuit parameters of the generator to be supplied (at the operating temperature, rather than the ‘cold’ values), along with the number of pole pairs, the voltage and the network frequency. It is also possible to model power factor correction capacitors and auxiliary loads such as turbine ancillary equipment. The equivalent circuit configuration is shown in Figure 4.1. Rr/s

Rs xs

xr

xm

Ra C Xa

Rs = Stator resistance xs = Stator reactance Rr = Rotor resistance xr = Rotor reactance xm = Mutual reactance C = Power factor correction Ra = Auxiliary load resistance Xa = Auxiliary load reactance s = slip

Figure 4.1: Equivalent circuit model of induction generator The equivalent circuit parameters should be given for a star-connected generator. If the generator is delta-connected, the resistances and reactances should be divided by 3 to convert to the equivalent star-connected configuration. The voltage should be given as rms line volts. To convert peak voltage to rms, divide by 2. To convert phase volts to line volts, multiply by 3. Since this model necessarily includes electrical losses in the generator and ancillary equipment, it is not possible to specify any additional electrical losses, although mechanical losses may be specified - see Section 4.4. Four different models of the electrical dynamics of the system illustrated in Figure 4.1 are provided: • • • •

Steady state 1st order 2nd order 4th order

The steady state model simply calculates the steady-state currents and voltages in Figure 4.1 at each instant. The 1st order model introduces a first order lag into the relationship between the slip (s) and the effective rotor resistance (Rr/s), using the short-circuit transient time constant given by [4.1]:

=

X s X r x 2m XsR r s

where Xs = xs+xm, Xr = xr+xm, and

s

is the grid frequency in rad/s.

21 of 82



ISSUE:011

The 2nd order model represents the generator as a voltage source reactance X’ = Xs - xm2/Xr, ignoring stator flux transients:

FINAL

behind a transient

is (rs + jX’) = vs where is and vs are the stator current and terminal voltage respectively. The dynamics of the rotor flux linkage r may be written as

1 & r = rr i r + js ( + s ) 1 s

r

where s is the fractional slip speed (positive for generating) and ir is the rotor current. This can be re-written in terms of the induced voltage using xm r = j Xr to give

T0 & =

rs + jX s rs + jX

js

s

T0

+j

Xs X vs rs + jX

where T0 =

Xr . s rr

The 4th order model is a full d-q (direct and quadrature) axis representation of the generator which uses Park’s transformation [4.2] to model the 3-phase windings of the generator as an equivalent set of two windings in quadrature [4.3]. Using complex notation to represent the direct and quadrature components of currents and voltages as the real and imaginary parts of a single complex quantity, we can obtain

xsx r s

x 2m d i s dt i r

=

x r rs + jx 2m (1 + s) x m rs

jx m x s (1 + s)

x m rr + jx m xr (1 + s) x s rr

jx s x r (1 + s)

is ir

+

xr v xm s

where all the currents and voltages are now complex. Where speed of simulation is more important than accuracy, one of the lower order models should be used. The 4th order model should be used for the greatest accuracy, although in many circumstances the lower order models give very similar results. The lower order models do not give an accurate representation of start-up transients, however. 4.2.3 Variable speed generator This model should be used for a variable speed turbine incorporating a frequency converter to decouple the generator speed from the grid frequency. The variable speed drive, consisting of both the generator and frequency converter, is modelled as a whole. A modern variable speed drive is capable of accepting a torque demand and responding to this within a very short time to give the desired torque at the generator air-gap, irrespective of the generator speed (as long as it is within specified limits). A first order lag model is provided for this response:

22 of 82


Qg =


ISSUE:011

FINAL

Qd (1 + e s)

where Qd is the demanded torque, Qg is the air-gap torque, and e is the time constant of the first order lag. Note that the use of a small time constant may result in slower simulations. If the time constant is very small, specifying a zero time constant will speed up the simulations, without much effect on accuracy. A variable speed turbine requires a controller to generate an appropriate torque demand, such that the turbine speed is regulated appropriately. Details of the control models which are available with Bladed can be found in Section 5. The minimum and maximum generator torque must be specified. Motoring may occur if a negative minimum torque is specified. The phase angle between current and voltage, and hence the power factor, is specified, on the assumption that, in effect, both active and reactive power flows into the network are being controlled with the same time constant as the torque, and that the frequency converter controller is programmed to maintain constant power factor. An option for drive train damping feedback is provided. This represents additional functionality which may be available in the frequency converter controller which adds a term derived from measured generator speed onto the incoming torque demand. This term is defined as a transfer function acting on the measured speed. The transfer function is supplied as a ratio of polynomials in the Laplace operator, s. Thus the equation for the air-gap torque Qg becomes

Qg =

Qd Num(s) + (1 + e s) Den(s)

g

where Num(s) and Den(s) are polynomials. The transfer function would normally be some kind of tuned bandpass filter designed to provide some damping for drive train torsional vibrations, which in the case of variable speed operation may otherwise be very lightly damped, sometimes causing severe gearbox loads. 4.2.4 Variable slip generator A variable slip generator is essentially an induction generator with a variable resistance in series with the rotor circuit [4.3, 4.4]. Below rated power, it acts just like a fixed speed induction generator, so the same parameters are required as described in Section 4.2.1. Above rated, the variable slip generator uses a fast-switching controller to regulate the rotor current, and hence the air-gap torque, so the generator actually behaves just like a variable speed system, albeit with a limited speed range. The same parameters as for a variable speed system must therefore also be supplied (see Section 4.2.3), with the exception of the phase angle since power factor control is not available in this case. Alternatively, a full electrical model of the variable slip generator is available. The generator is modelled as in Section 4.2.2, and the rotor current controller is modelled as a continuoustime PI controller which adjusts the rotor resistance between the defined limits (with

23 of 82



ISSUE:011

FINAL

integrator desaturation on the limits), in response to the difference between the actual and demanded rotor current. The steady-state relationship between torque and rotor current is computed at the start of the simulation, so that the torque demand can be converted to a rotor current demand. The scheme is shown in Figure 4.2. Torque demand

Current demand

1 |I|

PI with limits

Rotor resistanc e

Measured current |I| Figure 4.2: Variable slip generator – rotor current controller

4.3 Drive train mounting If desired, torsional flexibility may be specified either in the gearbox mounting or between the pallet or bedplate and the tower top. This option is only allowed if either the stiff or flexible drive train model is specified, and it adds an additional rotational degree of freedom. In either case, the torsional stiffness and damping of the mounting is specified, with the axis of rotation assumed to coincide with the rotor shaft. The moment of inertia of the moving components about the low speed shaft axis must also be specified. In the case of a flexible gearbox mounting, this is the moment of inertia of the gearbox casing. In the case of a flexible pallet mounting, it is the moment of inertia of the gearbox casing, the generator stator, the moving pallet and any other components rigidly fixed to it. If either form of mounting is specified, the direction of rotation of the generator shaft will affect some of the internal drive train loads. If the low speed and high speed shafts rotate in opposite directions, specify a negative gearbox ratio in the drive train model. The effect of any offset between the low speed shaft and high speed shaft axes is ignored. Any shaft brake is assumed to be rigidly mounted on the pallet. Thus any motion once the brake disk has stopped turning depends on the type of drive train mounting as well as on the position of the brake on the low or high speed shaft. For example if there is a soft pallet mounting, then there will still be some oscillation of the rotor after the brake disk has stopped even if both shafts are stiff. As in the case of the flexible shaft drive train model, it should be pointed out that while modelling the effect of flexible mountings provides greater accuracy in the prediction of loads, there is potential for one or two of the resulting drive train vibrational modes to be of relatively high frequency, depending on the various moments of inertia and shaft and mounting stiffnesses. The presence of high frequency modes could result in slower simulations.

4.4 Energy losses Power train energy losses are modelled as a combination of mechanical losses and electrical losses in the generator (including the frequency converter in the case of variable speed turbines).

24 of 82



ISSUE:011

FINAL

Mechanical losses in the gearbox and/or shaft bearings are modelled as either a loss torque or a power loss, which may be constant, or interpolated linearly from a look-up table. This may be a look-up table against rotor speed, gearbox torque or shaft power, or a two-dimensional look-up table against rotor speed and either shaft torque or power. Mechanical losses modelled in terms of power are inappropriate if calculations are to be carried out at low or zero rotational speeds, e.g. for starts, stops, idling and parked calculations. In these cases, the losses are better expressed in terms of torque. The electrical losses may specified by one of two methods: Linear model: This requires a no-load loss LN and an efficiency , where the electrical power output Pe is related to the generator shaft input power Ps by: Pe =

(Ps - LN)

Look-up table: The power loss L(Ps) is specified as a function of generator shaft input power Ps by means of a look-up table. The electrical power output Pe is given by: Pe = Ps - L(Ps) Linear interpolation is used between points on the look-up table. Note that if a full electrical model of the generator is used, additional electrical losses in this form cannot be specified since the generator model implicitly includes all electrical losses.

4.5 The electrical network Provided either the detailed electrical model of the induction generator or the variable speed generator model is used, so that electrical currents and voltages are calculated, and reactive power as well as active power, then the characteristics of the network to which the turbine is connected may also be supplied. As well as allowing the voltage variations, and hence the flicker, at various points on the network to be calculated, the presence of the network may also, in the case of the directly connected induction generator, influence the dynamic response of the generator itself particularly on a weak network. The network is modelled as a connection, with defined impedance, to the point of common coupling (PCC in Figure 4.2) and a further connection, also with defined impedance, to an infinite busbar. Further turbines may be connected at the point of common coupling. These additional turbines are each assumed to be identical to the turbine being modelled, including the impedance of the connection to the point of common coupling. However they are modelled as static rather dynamic, with current and phase angle constant during the simulation. The initial conditions are calculated with the assumption that all turbines are in an identical state, and the ‘other’ turbines then remain in the same state throughout. Thus the steady state voltage rise due to all the turbines at the point of common coupling will be taken into account in calculating the performance of the turbine whose performance is being simulated .

25 of 82



ISSUE:011

FINAL

Other turbines (if required) Wind turbine

R1 + jX1 Windfarm interconnection impedance

PCC

Figure 4.2: The network model

26 of 82

R2 + jX2 Network connection impedance

Infinite busbar



ISSUE:011

FINAL

5. CLOSED LOOP CONTROL

5.1 Introduction Closed loop control may be used during normal running of the turbine to control the blade pitch angle and, for variable speed turbines, the rotor speed. Four different controller types are provided: 1. Fixed speed stall regulated. The generator is directly connected to a constant frequency grid, and there is no active aerodynamic control during normal power production. 2. Fixed speed pitch regulated. The generator is directly connected to a constant frequency grid, and pitch control is used to regulate power in high winds. 3. Variable speed stall regulated. A frequency converter decouples the generator from the grid, allowing the rotor speed to be varied by controlling the generator reaction torque. In high winds, this speed control capability is used to slow the rotor down until aerodynamic stall limits the power to the desired level. 4. Variable speed pitch regulated. A frequency converter decouples the generator from the grid, allowing the rotor speed to be varied by controlling the generator reaction torque. In high winds, the torque is held at the rated level and pitch control is used to regulate the rotor speed and hence also the power. For a constant speed stall regulated turbine no parameters need be defined as there is no control action. In the other cases the control action will determine the steady state operating point of the turbine as well as its dynamic response. For steady state calculations it is only necessary to specify those parameters which define the operating curve of the turbine. For dynamic calculations, further parameters are used to define the dynamics of the closed loop control. The parameters required are defined further in the following sections. Note that all closed loop control data are defined relative to the high speed shaft.

5.2 The fixed speed pitch regulated controller This controller is applicable to a turbine with a directly-connected generator which uses blade pitch control to regulate power in high winds. It is applicable to full or partial span pitch control, as well as to other forms of aerodynamic control such as flaps or ailerons. In the latter case, the pitch angle can be taken to refer to the deployment angle of the flap or aileron. From the optimum position, the blades may pitch in either direction to reduce the aerodynamic torque. If feathering pitch action is selected, the pitchable part of the blade moves to reduce its angle of attack as the wind speed (and hence the power) increases. If stalling pitch action is selected, it moves in the opposite direction to stall the blade as the wind speed increases. In the feathering case, the minimum pitch angle defines the pitch setting below rated, while in the stalling case the maximum pitch angle is used below rated, and the pitch decreases towards the minimum value (usually a negative pitch angle) above rated.

27 of 82


Wind


Electric power

Turbine

Blade pitch

ISSUE:011

FINAL

Power Measured transducer power Controlle

Pitch actuator

Pitch demand

Power set-point

Figure 5.1: The fixed speed pitch regulated control loop

Figure 5.1 shows schematically the elements of the fixed speed pitch regulated control loop which are modelled. 5.2.1 Steady state parameters In order to define the steady-state operating curve, it is necessary to define the power setpoint and the minimum and maximum pitch angle settings, as well as the direction of pitching as described above. The correct pitch angle can then be calculated in order to achieve the setpoint power at any given steady wind speed. 5.2.2 Dynamic parameters To calculate the dynamic behaviour of the control loop, it is necessary to specify the dynamic response of the power transducer and the pitch actuator, as well as the actual algorithm used by the controller to calculate a pitch demand in response to the measured power signal. Section 5.5 describes the available transducer and actuator models, while Section 5.6 describes the PI algorithm which is used by the controller.

5.3 The variable speed stall regulated controller This controller model is appropriate to variable speed turbines which employ a frequency converter to decouple the generator speed from the fixed frequency of the grid, and which do not use pitch control to limit the power above rated wind speed. Instead, the generator reaction torque is controlled so as to slow the rotor down into stall in high wind speeds. The control loop is shown schematically in Figure 5.2. 5.3.1 Steady state parameters The steady-state operating curve can be described with reference to a torque-speed graph as in Figure 5.3. The allowable speed range in the steady state is from S1 to S2. In low winds it is possible to maximise energy capture by following a constant tip speed ratio load line which corresponds to operation at the maximum power coefficient. This load line is a quadratic curve on the torque-speed plane, shown by the line BG in Figure 5.3. Alternatively a look-up table may be specified. If there is a minimum allowed operating speed S1, then it is no longer possible to follow this curve in very low winds, and the turbine is then operated at nominally constant speed along the line AB shown in the figure. Similarly in high wind speeds, once the maximum operating speed S4 is reached, then once again it is necessary to 28 of 82


Wind


ISSUE:011

Generator speed

Speed transducer

Electrical power

Power Measured Controlle transducer power

FINAL

Measured speed

Turbine

Generator torque demand

Desired power, torque, speed

Figure 5.2: The variable speed stall regulated control loop depart from the optimum load line by operating at nominally constant speed along the line GH. Once maximum power is reached at point H, it is necessary to slow the rotor speed down into stall, along the constant power line HI. If high rotational speeds are allowed, it is of course possible for the line GH to collapse so that the constant power line and the constant tip speed

Figure 5.3: Variable speed stall regulated operating curve

29 of 82



ISSUE:011

FINAL

ratio line meet at point J. Clearly the parameters needed to specify the steady state operating curve are: • The minimum speed, S1 • The maximum speed in constant tip speed ratio mode, S4 • The maximum steady-state operating speed. This is usually S4, but could conceivably be higher in the case of a turbine whose characteristics are such that as the wind speed increases, the above rated operating point moves from H to I, then drops back to H, and then carries on (towards J) in very high winds. This situation is somewhat unlikely however, because if rotational speeds beyond S4 are permitted in very high winds, there is little reason not to increase S4 and allow the same high rotor speeds in lower winds.) • The above rated power set-point, corresponding to the line HI. This is defined in terms of shaft power. Electrical power will of course be lower if electrical losses are modelled. • The parameter K which defines the constant tip speed ratio line BG. This is given by: K =

R5 Cp( ) / 2

3

G3

where = air density R = rotor radius = desired tip speed ratio Cp( ) = Power coefficient at tip speed ratio G = gearbox ratio Then when the generator torque demand is set to K 2 where is the measured generator speed, this ensures that in the steady state the turbine will maintain tip speed ratio and the corresponding power coefficient Cp( ). Note that power train losses may vary with rotational speed, in which case the optimum rotor speed is not necessarily that which results in the maximum aerodynamic power coefficient. As an alternative to the parameter K , a look-up table may be specified giving generator torque as a function of speed. 5.3.2 Dynamic parameters To calculate the dynamic behaviour of the control loop, it is necessary to specify the dynamic response of both power and speed transducers, as well as the actual algorithm used by the controller to calculate a generator torque demand in response to the measured power and speed signals. Section 5.5 describes the available transducer and actuator models. Two closed loop control loops are used for the generator torque control, as shown in Figure 5.4. An inner control loop calculates a generator torque demand as a function of generator speed error, while an outer loop calculates a generator speed demand as a function of power error. Both control loops use PI controllers, as described in Section 5.6. Below rated, the speed set-point switches between S1 and S4. In low winds it is at S1, and the torque demand output is limited to a maximum value given by the optimal tip speed ratio curve BG. This causes the operating point to track the trajectory ABG. In higher winds, the set-point changes to S4, and the torque demand output is limited to a minimum value given by the optimal tip speed ratio curve, causing the operating point to track the trajectory BGH. Once the torque reaches QR, the outer control loop causes the speed set-point to reduce along HI, and the inner loop tracks this varying speed demand. 30 of 82



PI controller

Power set-point

ISSUE:011

FINAL

PI controller

Speed demand

Measured power Measured speed Generator torque demand Figure 5.4: Stall regulated variable speed control loops

5.4 The variable speed pitch regulated controller This controller model is appropriate to variable speed turbines, which employ a frequency converter to decouple the generator speed from the fixed frequency of the grid, and which use pitch control to limit the power above rated wind speed. The control loop is shown schematically in Figure 5.5. 5.4.1 Steady state parameters The steady-state operating curve can be described with reference to the torque-speed graph shown in Figure 5.6. Below rated, i.e. from point A to point H, the operating curve is exactly as in the stall regulated variable speed case described in Section 5.3.1, Figure 5.3. Above rated however, the blade pitch is adjusted to maintain the chosen operating point, designated

Wind

Generator speed

Speed transducer

Measured speed

Turbine Controlle Blade pitch

Pitch actuator

Pitch demand

Generator torque demand

Desired torque and speed

Figure 5.5: The variable speed pitch regulated control loop L. Effectively, changing the pitch alters the lines of constant wind speed, forcing them to pass through the desired operating point.

31 of 82



ISSUE:011

FINAL

Figure 5.6: Variable speed pitch regulated operating curve

Once rated torque is reached at point H, the torque demand is kept constant for all higher wind speeds, and pitch control regulates the rotor speed. A small (optional) margin is allowed between points H (where the torque reaches maximum) and L (where pitch control begins) to prevent excessive mode switching between below and above rated control modes. However, this margin may not be required, in which case points H and L coincide. As with the stall regulated controller, the line GH may collapse to a point if desired. Clearly the parameters needed to specify the steady state operating curve are: • The minimum speed, S1 • The maximum speed in constant tip speed ratio mode, S4 • The speed set-point above rated (S5). This may be the same as S4. • The maximum steady-state operating speed. This is normally the same as S5. • The above rated torque set-point, QR. • The parameter K which defines the constant tip speed ratio line BG, or a look-up table. This is as defined in Section 5.3.1. 5.4.2 Dynamic parameters To calculate the dynamic behaviour of the control loop, it is necessary to specify the dynamic response of the speed transducer and the pitch actuator, as well as the actual algorithm used by the controller to calculate the pitch and generator torque demands in response to the measured speed signal. Section 5.5 describes the available transducer and actuator models.

32 of 82



ISSUE:011

FINAL

Figure 5.7 shows the control loops used to generate pitch and torque demands. The torque demand loop is active below rated, and the pitch demand loop above rated. Section 5.6 describes the PI algorithm which is used by both loops. Below rated, the speed set-point switches between S1 and S4. In low winds it is at S1, and the torque demand output is limited to a maximum value given by the optimal tip speed ratio curve BG. This causes the operating point to track the trajectory ABG. In higher winds, the set-point changes to S4, and the torque demand output is limited to a minimum value given by the optimal tip speed ratio curve, causing the operating point to track the trajectory BGH, and a maximum value of QR. When point H is reached the torque remains constant, with the pitch control loop becoming active when the speed exceeds S5. Above rated Speed set-point Below rated

Measured speed Blade

PI controller

PI controller

pitch

Generator torque demand Figure 5.7: Pitch regulated variable speed control loops

5.5 Transducer models First order lag models are provided in Bladed to represent the dynamics of the power transducer and the generator speed transducer. The first order lag model is represented by y& =

1 (x T

y)

where x is the input and y is the output. The input is the actual power or speed and the output is the measured power or speed, as input to the controller.

5.6 Modelling the pitch actuator The pitch actuator may be modelled as either a pitch position or pitch rate actuator, and either active or passive dynamics may be specified. The simplest model is a passive actuator, with the relationship between the input and the output represented by a transfer function. For the pitch position actuator, the input is the pitch demand generated by the controller and the output is the actual pitch angle of the blades. For the pitch rate actuator, the input is the pitch rate demand generated by the controller and the output is the actual pitch rate at which the blades move. The transfer

33 of 82



ISSUE:011

FINAL

function may be a first order lag, a second order response, or a general transfer function, up to 8th order. The first order lag model is represented by y& =

1 (x T

y)

where x is the input and y is the output. The second-order model is represented by && y + 2! y& =

2

(x

y)

where is the bandwidth and ! the damping factor. The general transfer function model is represented by numerator and denominator polynomials in the Laplace operator. For detailed calculations, especially to understand the loads on the pitch actuator itself and the duty which will be required of it, it is possible to enter a more detailed model. This can take into account any internal closed loop dynamics in the actuator, and also the pitch motion resulting from the actuator torque acting on the pitching inertia, with or against the aerodynamic pitch moment and the pitch bearing friction. The bearing friction itself depends critically on the loading at the pitch bearing. Figure 5.8 shows the various options for controlling the pitch angle, starting from either a pitch position demand or a pitch rate demand. The pitch position demand may optionally be processed through a ramp control, shown in Figure 5.9, which smooths the step changes in demand generated by a discrete controller by applying rate and/or acceleration limits. Then the pitch position demand can act either through passive dynamics to generate a pitch position, or through a PID controller on pitch error to generate a pitch rate demand. Rate limits are applied to the output, with instantaneous integrator desaturation to prevent wind-up in the PID case. Thus the pitch rate demand may come either from here or directly from the controller. This rate demand can act either through passive dynamics to generate a pitch rate, or through a PID controller on pitch rate error to generate an actuator torque demand. In the latter case, the pitch actuator passive dynamics then generate an actual actuator torque, which acts against bearing friction and any aerodynamic pitching moment to accelerate the pitching inertia of the blades and the actuator itself. An optional first order filter on each PID input allows step changes in demand from the controller to be smoothed, and instantaneous integrator desaturation prevents wind-up when the torque limits are reached. Both PID controllers include a filter on the differential term to prevent excessive high frequency gain. Also there is a choice of derivative action, such that the derivative gain may be applied either to the feedback (i.e. the measured position or rate), the error signal, or the demand. The latter case represents a feed-forward term in the controller. If passive pitch rate dynamics are selected, the response will be subject to acceleration limits calculated from the aerodynamic pitching moment, bearing friction and the actuator toque limits acting on the pitching inertia. If the total pitching inertia is zero, no limits will be applied. The pitch bearing sliding friction torque is modelled as the sum of four terms: a constant, a term proportional to the bending moment at the bearing, and a terms proportional to the axial and radial forces on the bearing. Sometimes the actuator cannot overcome the applied

34 of 82



ISSUE:011

FINAL

torques and the pitch motion will stick. Before it can move again, the break-out or ‘stiction’ torque must be overcome. This is modelled as an additional contribution to the friction torque while the pitch is not moving. This additional contribution is specified as a constant torque, plus a term proportional to the sliding friction torque.

Pitch position demand from controller

Measured pitch position

Pitch rate demand from controller

Measured pitch rate

Bearing loads

Ramp control

+

Pitching moment PID controller

Actuator torque limits Pitch rate demand + Pitching inertia Acceleration limits

Passive dynamics

PID controller

Actuator torque demand

Passive dynamics

Passive dynamics Actuator torque Pitching inertia

Pitch rate

Actual pitch position

Figure 5.8: Pitch actuator options

35 of 82

Bearing friction



ISSUE:011

FINAL

1.2

1

Demand

0.8

0.6

0.4

Raw demand Rate limit

0.2

Acceleration limit Rate & acceleration limits

0 -0.2

0

0.2

0.4

0.6

0.8

1

1.2

Timesteps

The ramp is re-started each timestep. If the ramp is not completed by the end of the timestep and an acceleration limit is specified, the slope at the start of the next timestep will be nonzero. Figure 5.9: Ramp control for pitch actuator position demand

5.7 The PI control algorithm All the closed loop control algorithms described above use PI controllers to calculate the output y (pitch, torque or speed demand) from the input x (power or speed error). The basic PI algorithm can be expressed as

y& = K p x& + Ki x where Kp and Ki represent the proportional and integral gains. The ratio Kp/Ki is also known as the integral time constant. Calculation of appropriate values for the gains is a specialist task, which should take into account the dynamics of the wind turbine together with the aerodynamic characteristics and principal forcing frequencies, and should aim to achieve stable control at all operating points and a suitable trade-off between accuracy of tracking the set-point and the degree of actuator activity. Straightforward implementation of the above equation leads to the problem of ‘integrator wind-up’ if the output y is subject to limits, as is the case here. This means that the raw output calculated as above continues to change as a result of the integral (Ki) term even though the actual output is being constrained to a limit. When the direction of movement of y changes, it will then take a long time before it comes back to the limit so that the final (constrained) output starts to change. This is avoided in the continuous-time implementation of the PI controller by an additional term -"y/Td in the above equation, where "y is the

36 of 82



ISSUE:011

FINAL

amount by which the raw output y has gone beyond the limit, and Td is the desaturation time constant which must be supplied by the user. In practice the control algorithm is usually implemented in a digital controller working on a discrete timestep. In the Bladed model, the continuous implementation of the controller is an approximate representation, although the discrete timestep is usually fast enough for the approximation to be a very good one. Since the integrator desaturation in a discrete controller can be implemented by fully adjusting the raw integrator output at every timestep, a suitable approximation for the continuous case is to use a desaturation time constant approximately equal to the discrete controller timestep. Alternatively, perfect or instantaneous desaturation can be specified by setting the desaturation time constant to zero. 5.7.1 Gain scheduling Since the characteristics of the turbine, especially the aerodynamic characteristics, are not constant but will vary according to the operating point, and hence the wind speed, it may be necessary to adjust the controller gains as a function of the operating point in order to ensure that suitable control loop characteristics are achieved at all wind speeds. This is known as gain scheduling, and the gain scheduling model provided in Bladed allows both the proportional and integral gains of any control loop to be scaled by a factor 1/F, where F is a function of some variable V which is accessible to the controller and which is representative of the operating point in some way. The choices available are: • F = constant • F = F(V) as defined by a look-up table • F = F(V) as defined by a polynomial, but with minimum and maximum limits applied to F The choice of variable V depends on the particular control loop. The following choices are provided: Fixed speed pitch regulated controller: • Electrical power, pitch angle, wind speed. Variable speed below-rated torque controller: • Electrical power, generator speed, wind speed, and pitch angle (in the pitch regulated case). Variable speed stall regulated above-rated controller: • Electrical power, generator speed, wind speed. Variable speed pitch regulated above-rated controller: • Electrical power, generator speed, wind speed, pitch angle. The variables shown in bold are normally recommended. Gain scheduling is unlikely to be required for the variable speed below rated controllers. For the variable speed stall regulated above-rated controller, no general rule can be given. Gain scheduling on wind speed is not usually a practical proposition because of the difficulty of measuring a representative wind speed, and this option is only provided for research purposes. The wind speed used is the

37 of 82



ISSUE:011

FINAL

hub wind speed, which may differ from any wind speed measured by an anemometer mounted on the nacelle, especially in the case of an upwind turbine. Gain scheduling on pitch angle is recommended for the pitch regulation controllers, to compensate for the large changes in the sensitivity of aerodynamic torque to pitch angle over the operating range. The steady loads calculation may be used to calculate the partial derivative of aerodynamic torque with respect to pitch angle, and F may be set proportional to this. In many cases, simply setting F proportional to pitch angle is a good approximation, but a lower limit for F must be set to prevent excessive gains at small pitch angles.

5.8 Control mode changes The variable speed controllers, both stall regulated and pitch regulated, require the following mode changes: • Change of speed set-point from S1 to S4 (refer to Figures 5.3 and 5.6). This occurs when the measured speed crosses the threshold value (S1+S4)/2. This mode change is completely benign as the control action along the optimum tip speed ratio line BG is the same either side of the mode change point, so no hysteresis is required. • Change from below rated to above rated control. For the stall regulated case, the change from below rated to above rated is also benign. Making the switch in the middle of the section GH of Figure 5.3 causes no immediate change in control action. However, in the case of G and H coinciding, or being very close together, it may be necessary to modify the mode change strategy, depending on the turbine characteristics. For the pitch regulated case, the change to above rated control occurs when the torque demand is at maximum (QR) and the speed exceeds S5 (refer to Figure 5.6). The change to below rated occurs when the pitch demand is at fine pitch (minimum pitch for the feathering case, maximum pitch for pitch-assisted stall) and the speed falls below S4. While this strategy is usually suitable, it may be desirable to modify it depending on the turbine characteristics. The mode changes occur on a discrete timestep set to a default value of 0.1 seconds.

5.9 Client-specific controllers The control algorithms described above have been developed to be suitable for a wide range of cases. However, it is recognised that there is great variation in the design of controllers for wind turbines. In a number of specific cases, Garrad Hassan have enhanced these basic controller designs in various ways to suit particular turbine designs, further improving the control performance. In many cases, particularly for variable speed turbines, both the closed loop performance and the mode changing behaviour can be improved significantly with a small additional degree of sophistication.

38 of 82



ISSUE:011

FINAL

It is also recognised that simulations must be able to be adapted to use any particular controller design, both to allow algorithms different from the standard ones described above, and also to allow the modelling of discrete controllers, for example so that the effect of controller timestep can be investigated. For these reasons, Bladed offers the possibility of incorporating user-defined controllers in the dynamic simulations. Through a defined interface which makes use of a shared file, a user’s control program, written in any language, can be used to control the simulation. The user-defined controller may do any of the following: • Blade pitch angle or blade pitch rate control during any phase of operation including power production, stops, starts, idling etc. • Generator torque control for variable speed turbines • Control the generator contactor, allowing the generator to be switched on or off for simulating stops and starts • Control the shaft brake, to simulate transitions between parked, idling, starting, stopping, and power production states. • Control of nacelle yaw to simulate closed loop yaw control algorithms and/or yawing strategies for start-up, shutdown etc. The User Manual describes how to write a user-defined control program.

5.10 Signal noise and discretisation When a discrete external controller is used, Bladed offers the possibility of adding random noise to the measured signals sent to the controller, and also to discretise the signals to a specified resolution. The random noise may be Gaussian, in which case the standard deviation of the noise must be specified, or it may be from a rectangular distribution, in which case the half-width of the distribution should be given. The noise is added to the signal before it is discretised.

39 of 82



ISSUE:011

FINAL

6. SUPERVISORY CONTROL This section of the manual covers the modelling of the following aspects of turbine operation: • • • • • •

Start-up Shut-down (normal and emergency stops) Non-operational situations (rotor parked or idling) Operation of the shaft brake Teeter restraints Yaw control

The standard implementation of these features in the simulation model is described. As in the case of Closed Loop Control, alternative supervisory control logic can be incorporated in a user-defined controller - see Section 5.9.

6.1 Start-up Simulation of a wind turbine start-up begins with the rotor at a specified speed (usually but not necessarily zero) and the generator off-line. The brake is assumed to be released at the start of the simulation (i.e. at time zero). If blade pitch or aileron control is available, the initial pitch or aileron angle is specified, along with a constant rate of change which continues until either a specified angle is reached or the closed loop controller takes over. When a specified rotational speed is reached, the generator comes on line, and the closed loop controller begins to operate. The simulation continues until the specified simulation end time. In the case of a variable speed turbine, there may be a transition period after cut-in of the closed loop controller before the turbine is fully in the normal running state. There are two different cases: Variable speed pitch regulation: in the case when the pitch angle has not yet reached the normal operating value (‘fine pitch’) at the moment when the closed loop controller cuts in, then the pitch change rate for start-up continues to apply until either fine pitch is reached, or until the conditions of Section 5.8 for starting the closed loop pitch controller are satisfied. Variable speed stall regulation: when the closed loop controller cuts in, the above-rated control mode is assumed to apply initially. In practice this assumption does not affect the start-up since in low winds the operating point would be constrained by the quadratic optimum-Cp characteristic in any case.

40 of 82



ISSUE:011

FINAL

6.2 Normal stops A normal stop is initiated at a specified time after the start of the simulation. Normal operation in power production mode is assumed prior to this point, with full structural and control dynamics in effect if desired. The structural dynamics continue in effect during the entire simulation. The standard logic for a normal stop is to start pitching the blades (or moving the ailerons) at a specified rate from the moment that the stop is initiated, continuing until a final pitch angle is reached. The generator is taken off-line when the electrical power reaches zero in the case of a fixed-speed turbine, or when the minimum generator speed is reached in the case of a variable speed turbine. Once the rotational speed drops below a specified value, the shaft brake is applied to bring the rotor to rest. The simulation continues until the rotor comes to rest, or for a certain time longer if so desired in order that the transient loads can be simulated as the brake disk stops. However, the simulation end time overrides this, so it must be set long enough for the stop event to be completed. If there is no pitch control, the brake trip speed may be set high so that the shaft brake is applied immediately at the initiation of the stop. Section 6.4 describes the dynamic characteristics of the shaft brake itself.

6.3 Emergency stops An emergency stop is initiated at a specified time after the start of the simulation. Normal operation in power production mode is assumed prior to this point, with full structural and control dynamics in effect if so desired. The structural dynamics continue in effect during the entire simulation. Several options are available for simulating emergency stops. In all cases it is assumed that the generator load is lost at the initiation of the emergency stop, whether because of grid failure or some electrical or mechanical failure of the turbine. Pitch (or aileron) action is initiated either immediately or when the rotational speed exceeds a specified value. A fixed pitch rate then applies until a final pitch angle is reached. Provision is made for the pitch of one or more of the blades to ‘stick’ at a specified angle to simulate failure of a pitch bearing or actuator. The shaft brake can also be applied either at the initiation of the stop or when a specified overspeed is reached. Section 6.4 describes the dynamic characteristics of the shaft brake itself. There is also a rotational speed below which the shaft brake is applied for parking, in the event that it has not already been applied because of load loss or overspeed. The simulation continues until the rotor comes to rest, or for a certain time longer if so desired in order that the transient loads can be simulated as the brake disk stops. However,

41 of 82



ISSUE:011

FINAL

the simulation end time overrides this, so it must be set long enough for the stop event to be completed.

6.4 Brake dynamics When the shaft brake is applied, either during a normal or an emergency stop, the full braking torque is not available instantly. Instead, the torque builds up to the full value over a short period of time. This torque build-up may be modelled as either a linear torque ramp, or by specifying a look-up table giving achieved braking torque as a function of time.

6.5 Idling and parked simulations For simulations in the idling and parked states, a fixed pitch angle is specified, the generator is off line, and there is no pitch control action. In the case of a parked rotor the shaft brake is applied, and the rotor azimuth must be specified. The azimuth is measured from zero with blade 1 at top dead centre. All specified structural dynamics will be in effect during these simulations. This also allows for the possibility of the shaft brake slipping during a parked simulation if the shaft torque exceeds the specified brake torque.

6.6 Yaw control 6.6.1 Active yaw Active yaw movement may be specified in one of two ways: 1. One fixed-rate yaw manoeuvre may be specified, starting at a given point in any simulation. This represents a change in the nominal nacelle position through a given angle at a specified angular speed. 2. A user-defined controller (Section 5.9) may be used to specify either the yaw rate or the yaw actuator torque at any time. If active yaw is used to control the yaw rate, the effect of this is to change the ‘demanded nacelle angle’ in a specified way. The actual nacelle angle depends on the yaw dynamics see next section.

42 of 82



ISSUE:011

FINAL

6.6.2 Yaw dynamics Three options are available to define the yaw dynamics: 1. Rigid yaw: the actual nacelle angle exactly follows the ‘demanded nacelle angle’ 0. 2. Flexible yaw: a certain amount of flexibility is present, usually in the yaw actuation system, such that the actual nacelle angle may not follow the ‘demanded nacelle angle’ 0 exactly. The extreme case is free yaw, when the demanded nacelle angle does not have any effect. 3. Controlled yaw torque: this is available only with an external controller to define the yaw actuator torque demand

Demanded yaw rate Aerodynamic and inertial yaw torque

Tower Controlled torque

Yaw spring

Friction

Yaw control type None Rigid Flexible Controlled torque

Demanded yaw rate No Yes Yes No

Damper

Yaw spring and damper No No Yes No

Friction No No Yes Yes

Controlled torque No No No Yes

In the case of flexible or free yaw, the yaw damping Dy may be specified. This specifies a torque Qd which opposes the yaw motion, given by

Qd = Dy ( & 0

&)

In the case of flexible yaw, a yaw spring may be specified either as a linear spring or as a hydraulic accumulator system such as is often used to provide flexibility in hydraulic yaw drives. The hydraulic system is assumed to be double-acting, with one accumulator (or set of accumulators) on either side of the yaw motor. The torque opposing the motion is provided by compression of the gas in the accumulators. If the nominal gas volume is V0 and the instantaneous gas volumes either side of the yaw motor are v1 and v2 then the opposing torque Qk is given by

Q k = KP0

V0 v1

#

V0 v2

#

43 of 82



ISSUE:011

FINAL

where v1 = V0 + F( - 0 ) and v2 = V0 - F( - 0 ) and P0 is the equilibrium pressure in the hydraulic system. The constant K defines the relationship between the torque developed at the yaw bearing and the pressure difference across the yaw motor, while F the relationship between the volume of oil flowing through the yaw motor and the resulting angular movement at the yaw bearing. # is the gas law constant: PV# = RT. Putting # = 1 specifies isothermal conditions in the accumulators.

6.7 Teeter restraint Although not strictly a supervisory control function, the teeter restraint model available in Bladed for teetered rotors is described here. The model allows a linear variation of restoring torque with teeter angle, but also allows a free teeter range and an initial pre-load. Figure 6.1 defines the relevant parameters. Linear damping is also allowed, giving an additional torque contribution proportional to teeter rate. Restoring torque

Pre-load

Spring constant

free teeter angle

Figure 6.1: Teeter restraint model

44 of 82

Teeter angle



ISSUE:011

FINAL

7. MODELLING THE WIND The wind field incident on the turbine may be specified in a number of ways. For some simple calculations, a uniform, constant wind speed is assumed, such that the same incident wind speed is seen by every point on the rotor. For more detailed calculations however, it is important to be able to define both the spatial and temporal variations in wind speed and direction. The steady-state spatial characteristics of the wind field may include any combination of the following elements: • Wind shear: the variation of wind speed with height. • Tower shadow: distortion of the wind flow by the wind turbine tower. • Upwind turbine wake: full or partial immersion of the turbine rotor in the wake of another turbine operating further upwind. The wind direction must also be specified, both relative to the direction in which the nacelle is pointing (to define the yaw error), and relative to the horizontal plane (to define the upflow angle). The latter effect may be important for turbines operating in hilly terrain. For simulations, it is also important to be able to define how the wind speed and direction vary with time. The following alternative models are provided: • Constant wind: no variation with time. • Single point history: a time history of wind speed and direction, which is fully coherent over the whole rotor, is specified as a look-up table against time. Linear interpolation is used between the time points. • 3D turbulent wind: this option uses a 3-dimensional turbulent wind field with defined spectral and spatial coherence characteristics representative of real atmospheric turbulence. This option will give the most realistic predictions of loads and performance in normal conditions. • IEC transients: this option uses wind speed and direction transients as defined by the IEC 1400-1 standard [7.1, 7.7]. It is intended for evaluating specific load cases, for example during extreme gusts.

45 of 82



ISSUE:011

FINAL

7.1 Wind shear Wind shear is the variation of steady state mean wind speed with height. Two alternative models are provided, to relate the wind speed V(h) at height h above the ground to the wind speed V(h0) at some reference height h0.. 7.1.1 Exponential model This model is defined in terms of a wind shear exponent $:

h V (h) = V (h0 ) h0

$

Specifying the exponent as zero results in no wind speed variation with height. 7.1.2 Logarithmic model This model is defined in terms of the ground roughness length z0:

V (h) = V (h0 )

log(h / z0 ) log(h0 / z0 )

7.2 Tower shadow Tower shadow defines the distortion of the steady-state mean wind field due to the presence of the wind turbine tower. Three different models are available: a potential flow model for upwind rotors, an empirical tower wake model for downwind rotors, and a combined model which is useful if the rotor yaws in and out of the downwind shadow area. 7.2.1 Potential flow model This model is appropriate for rotors operating upwind of the tower. The longitudinal wind velocity component upwind of the tower (V0) is modified using the assumption of incompressible laminar flow around a cylinder of diameter D = F.DT where DT is the tower diameter at the height where the tower shadow is being calculated, and F is a tower diameter correction factor supplied by the user. For a point at a distance z in front of the tower centreline and x to the side of the wind vector passing through the centreline, the wind speed V is given by:

V ( x , z ) = AV0 where

A = 1+

D 2

2

( x2 z2 ) ( x2 + z2 )2

provided the point is at an azimuth within +60° from bottom dead centre relative to the hub centre. For azimuth within +60° of top dead centre it is assumed that V(x,z) = V0 , and to 46 of 82



ISSUE:011

FINAL

ensure a smooth transition between these two zones, for all other azimuths, the factor A is modified to A( 0.5 cos( )) + (0.5 + cos( )) where is the blade azimuthal position 7.2.2 Empirical model For rotors operating downwind of the tower, an empirical model is provided, based on the work of Powles [7.2] which uses a cosine bell-shaped tower wake. For a point at a distance z behind the tower centreline and x to the side of the wind vector passing through the centreline, the wind speed V is given by:

V ( x , z ) = AV0 where A = 1 " cos 2

x WDT

for azimuth angles within +60° of bottom dead centre. For other azimuth angles, the same correction is applied as for the potential flow model, Section 7.2.2. Here " is the maximum velocity deficit at the centre of the wake as a fraction of the local wind speed, and W is the width of the tower shadow as a proportion of the local tower diameter DT. These quantities are defined for a given downwind distance, also expressed as a proportion of DT . For other distances, W increases, and " decreases, with the square root of the distance. 7.2.3 Combined model The combined model simply uses the potential flow model at the front and sides of the tower, and whichever of the other models gives the larger deficit at any point downwind. To ensure a smooth transition, the product of the A factors of the two models is used in any small areas where the potential flow model gives accelerated flow and the empirical model gives a velocity deficit.

7.3 Upwind turbine wake If the turbine rotor being modelled is assumed to be wholly or partially immersed in the wake of another turbine operating further upwind, a model is provided to define the modification to the steady-state mean wind profile caused by that wake. A Gaussian profile is used to describe the wake of the upstream turbine. The local velocity at a distance r from the wake centreline (which may be offset from the hub position) is given by: r2

V = V0 1 "e

2W2

47 of 82



ISSUE:011

FINAL

where V0 is the undisturbed wind speed, " is the fractional centre line velocity deficit, and W is the width of the wake (the distance from the wake centre line at which the deficit is reduced to exp(-0.5) times the centre line value). Two options are provided for defining the velocity deficit " and the wake width W . They can be defined directly, or they can be calculated by Bladed by specifying the characteristics of the upwind turbine. In the latter case, an eddy viscosity model of the wake is used, developed by Ainslie [7.8,7.9] and described in the next section. 7.3.1 Eddy viscosity model of the upwind turbine wake The eddy viscosity wake model is a calculation of the velocity deficit field using a finitedifference solution of the thin shear layer equation of the Navier Stokes equations in axissymmetric co-ordinates. The eddy viscosity model automatically observes the conservation of mass and momentum in the wake. An eddy viscosity, averaged across each downstream wake section, is used to relate the shear stress term in the thin shear equation to gradients of velocity deficit. The mean field can be obtained by a linear superposition of the wake deficit field and the incident wind flow. An illustration of the wake profile used in the eddy viscosity model is shown in Figure 7.1.

Figure 7.1: Wake profile used in the eddy viscosity model The Navier Stokes equations with Reynolds stresses and the viscous terms dropped gives [7.10]:

U

&U &U 1 &( ruv) +V = &x &r r &r

The turbulent viscosity concept is used to describe the shear stresses with an eddy viscosity defined by [7.11]:

48 of 82



ISSUE:011

FINAL

(x) = L m (x).U m (x) and

uv =

&U &r

Lm and Um are suitable length and velocity scales of the turbulence as a function of the downstream distance x but independent of r. The length scale is taken as proportional to the wake width Bw and the velocity scale is proportional to the difference UI – Uc across the shear layer. Thus the shear stress uv is expressed in terms of the eddy viscosity. differential equation to be solved becomes:

U

The governing

&U &U &( r&U / &r ) +V = &x &r r &r

Because of the effect of ambient turbulence, the eddy viscosity in the wake can not be wholly described by the shear contribution alone. Hence an ambient turbulence term is included, and the overall eddy viscosity is given by [7.12]:

= FK 1 B w ( U i

Uc ) +

amb

where the filter function F is a factor applied for near wake conditions. This filter can be introduced to allow for the build up of turbulence on wake mixing. The dimensionless constant K1 is a constant value over the whole flow field and a value of 0.015 is used. The ambient eddy viscosity term is calculated by the following equation proposed by Ainslie [7.12]: 2

amb

= F. K k . I amb / 100

Kk is the von Karman constant with a value of 0.4. Due to comparisons between the model and measurements reported by Taylor in [7.13] the filter function F is fixed at unity. The centre line velocity deficit Dmi can be calculated at the start of the wake model (two diameters downstream) using the following empirical equation proposed by Ainslie [7.12]: D mi = 1

Uc = Ct Ui

0.05

[(16C t

0.5)I amb /1000]

Assuming a Gaussian wind speed profile and momentum conservation an expression for the relationship between the deficit Dm and the width parameter Bw is obtained as

49 of 82


Bw =


ISSUE:011

FINAL

3.56C t 8D m (1 0.5D m )

Using the above equations, the average eddy viscosity at a distance 2D downstream of the turbine can be calculated. The equations can then be solved for the centre-line deficit and width parameter further downstream. Assuming to the Gaussian profile, the velocity deficit a distance r from the wake centreline is given by: r 3.56 Bw

D m ,r = exp

2

Therefore the wake width W used by Bladed is given by: W = Bw

0 .5 3.56

7.3.2 Turbulence in the wake If the eddy viscosity wake model is used, it is also possible to calculate the additional turbulence caused by the wake. The added turbulence is calculated using an empirical characterisation developed by Quarton and Ainslie [7.14]. This characterisation enables the added turbulence in the wake to be defined as a function of ambient turbulence Iamb, the turbine thrust coefficient Ct, the distance x downstream from the rotor plane and the length of the near wake, xn. The characterisation was subsequently amended slightly by Hassan [7.15] to improve the prediction, resulting in the following expression:

I add = 5.7C t 0.7 I amb 0.68 ( x / x n )

0.96

in which all turbulence intensities are expressed as percentages. Using the value of added turbulence and the incident ambient turbulence the turbulence intensity Itot at any turbine position in the wake can be calculated as

2

I tot = I amb + I add

2

The near wake length xn is calculated according to Vermeulen et al [7.16,7.17]: in terms of the rotor radius R and the thrust coefficient Ct as

xn =

n r0 dr dx

50 of 82



ISSUE:011

FINAL

where

r0 = R

1

m=

n=

m +1 2

1 Ct

(1

(

0.214 + 0.144m 1 0.214 + 0.144m

)

0.134 + 0.124m

)

0.134 + 0.124m

and dr/dx is the wake growth rate:

dr = dx dr dx dr dx and

dr dx

dr dx

2

+ $

dr dx

2

+ m

dr dx

2

= 2.5I 0 + 0.005 is the growth rate contribution due to ambient turbulence, $

= m

(1

m ) 1.49 + m is the contribution due to shear-generated turbulence, (1 + m ) 9.76

= 0.012 B

the number of blades and

is the contribution due to mechanical turbulence, where B is is the tip speed ratio.

7.4 Time varying wind Various forms of temporal variation of wind speed and direction may be superimposed on the spatial variations described in Sections 7.1 to 7.3 above. 7.4.1 Single point time history A look-up table can be used to supply the wind speed and direction as a function of time, at a defined reference height. Linear interpolation between time points is used. For any particular point in space, the wind speed is then multiplied by the appropriate correction factors for wind shear, tower shadow and upwind turbine wake as defined above. 7.4.2 3D turbulent wind A 3-dimensional turbulent wind field is generated, with statistical properties representative of real atmospheric turbulence. Section 7.5 describes how the turbulence is generated. It consists of dimensionless wind speed deviations, defined as + = (V-Vo)/IV0 where V0 is the mean wind speed and I the turbulence intensity, at a number of grid points on a rectangular array large enough to encompass the rotor swept area in the vertical and lateral (cross-wind)

51 of 82



ISSUE:011

FINAL

directions, and long enough in the longitudinal (along-wind) direction to allow a simulation of the desired length as the whole wind field moves past the rotor at the mean wind speed. At any point in time, the position in the longitudinal direction is calculated. The position in the lateral and vertical directions is calculated depending on the radial (r) and azimuthal position ( ) of any particular point on the rotor at that time, and 3-dimensional linear interpolation is then used to calculate the appropriate wind speed deviation +. The actual wind speed is then given by V(r, ,t) = V0Fs0 (Fs + I.+(r, ,t)) .FT .FW where Fs0 is the wind shear factor from the reference height (for mean speed V0 ) to the hub height, Fs is the wind shear factor from the hub height to the point (r, ), FT is the tower shadow factor for the point (r, ), and FW is the upwind turbine wake factor for the point (r, ). 7.4.3 IEC transients The transient variations of wind speed, shear and wind direction defined in the international standard for the safety of wind turbine systems, IEC 1400-1 [7.1, 7.7], may be simulated with Bladed. Transient changes in each of the following quantities may be independently simulated, each with its own parameter values: • • • •

Wind speed Wind direction Horizontal shear (linear variation of wind speed from one side of the rotor to the other) Vertical shear (linear variation of wind speed from bottom to top of the rotor)

Each may be either a half-wave transient or a full-wave transient. The transients are sinusoidal, with a more complex shape defined in edition 2 of the standard [7.7]. The parameters needed to define each transient are the starting value Y0, the start time t0, the duration T, and the amplitude A. These parameters are illustrated in Figure 7.2. 12.5

Half wave

Y0 + 12 A

11.5

11

Full wave 10.5

Y0

10

IEC edition 2

9.5

9 -0.2

t00

0.2

0.4

0.6

Time

0.8

t0 1+ T

1.2

Figure 7.2: Definition of IEC sinusoidal transients

52 of 82


The actual wind speed at radius r, azimuth


ISSUE:011

FINAL

and time t is then given by:

V(r, ,t) = (V0Fs0 Fs + Vtrans) .FT .FW where V0 is the starting wind speed at the reference height, Vtrans is the combined effect of the wind speed and horizontal and vertical shear transients, and other parameters as defined in Section 7.4.2.

7.5 Three dimensional turbulence model The wind simulation method adopted in Bladed is based on that described by Veers [7.3]. The rotor plane is covered by a rectangular grid of points, and a separate time history of wind speed is generated for each of these points in such a way that each time history has the correct single-point wind turbulence spectral characteristics, and each pair of time histories has the correct cross-spectral or coherence characteristics. Calculations using such a turbulent wind field will take into account the crucially important 'eddy slicing' transfer of rotor load from low frequencies to those associated with the rotational speed and its harmonics. This 'eddy slicing', associated with the rotating blades slicing through the turbulent structure of the wind, is a significant source of fatigue loading. The wind speed time histories may, in principle, be generated from any user-specified autospectral density and spatial cross-correlation characteristics. A choice of two different models of atmospheric turbulence has been provided. These are the von Karman and the Kaimal models. Both models are generally accepted as good representations of real atmospheric turbulence, although they use slightly different forms for the autospectral and cross-spectral density functions. The von Karman model can be used either to generate just the longitudinal component of turbulence, or to generate all three components if required. Two versions of the von Karman model are available: the basic model, given in [7.4] and described in Section 7.5.1, and the improved model, described in Section 7.5.2, which is based on more up-to-date information [7.5, 7.6]. The Kaimal model in Bladed gives only the longitudinal component of turbulence. It should be remembered, of course, that all these models tend to be based largely on observations for flat land sites. 7.5.1 The basic von Karman model The autospectral density for the longitudinal component of turbulence, according to the von Karman model, is given in [7.4] as

nSuu (n)

, 2u

=

4n~u (1 + 70.8n~u2 )5/ 6

where Suu is the auto-spectrum of wind speed variation, n is the frequency of variation, , u is

the standard deviation of wind speed variation and ~ n u is a non-dimensional frequency parameter given by:

53 of 82



ISSUE:011

FINAL

n xLu n~u = U Here xLu is the length scale of longitudinal turbulence and U is the mean wind speed. If the three-component model is selected, the corresponding spectra for the lateral (v) and vertical (w) components are:

nSii (n)

, i2

4n~i (1 + 755.2n~i2 ) = (1 + 282.3n~i2 ) 11/ 6

where

n x Li n~i = U and i is either v or w. Associated with the von Karman spectral equations is an analytical expression for the cross-correlation of wind speed fluctuations at locations separated in both space and time, derived assuming Taylor's frozen turbulence hypothesis. Accordingly for the longitudinal component at points separated by a distance "r perpendicular to the wind direction, the coherence Cu ("r,n), defined as the magnitude of the cross-spectrum divided by the autospectrum, is:

Cu ( "r , n) = 0.994( A5/ 6 (- u )

1 2

- u 5/ 3 A1/ 6 (- u ))

Here Aj(x) = xj Kj(x) where K is a fractional order modified Bessel function, and

nLu ( "r , n) "r - u = 0.747 1 + 70.8 Lu ( "r , n) U

2

The local length scale Lu("r,n) is defined by: L u ("r, n ) = 2MIN(1.0,0.04n

2/3

)

( y L u "y) 2 +( z L u "z) 2 "y 2 + "z 2

where "y and "z are the lateral and vertical components of the separation "r, and yLu and zLu are the lateral and vertical length scales for the longitudinal component of turbulence. For the lateral and vertical components, the corresponding equations are:

Ci ( "r , n) =

[

0.597 4.781# i2 A5/ 6 (- i ) 2 2.869# i 1

54 of 82

A11/ 6 (- i )

]



ISSUE:011

FINAL

where

nLi ( "r , n) "r - i = 0.747 1 + 70.8 Li ( "r , n) U and

#i =

2

- i Li ( "r , n) "

for i = v or w. In this case the local length scales are given by: L v ("r, n ) = 2MIN(1.0,0.05n

2/3

( y L v "y / 2) 2 +( z L v "z) 2

)

"y 2 + "z 2

and L w ("r, n ) = 2MIN(1.0,0.2n

1/ 2

)

( y L w "y) 2 +( z L w "z / 2) 2 "y 2 + "z 2

The three turbulence components are assumed to be independent of one another. This is a reasonable assumption, although in practice Reynolds stresses may result in a small correlation between the longitudinal and vertical components near to the ground. 7.5.2 The improved von Karman model The improved von Karman model [7.5] attempts to rectify some deficiencies of the basic model at heights below about 150m. The autospectral density for the longitudinal component of turbulence is given by:

nS uu (n)

, u2

= .1

(

2.987n~u / a 2 1 + ( 2 n~u / a )

)

5/ 6

+ .2

1294 . n~u / a

(1 + (

2 n~u / a )

)

5/ 6

F1

where Suu is the auto-spectrum of wind speed variation, n is the frequency of variation, ,u is the standard deviation of wind speed variation and n~u is a non-dimensional frequency parameter given by:

n xLu n~u = U Here xLu is the length scale of longitudinal turbulence and U is the mean wind speed. If the three-component model is selected, the corresponding spectra for the lateral (v) and vertical (w) components are:

55 of 82


nSii (n)

, i2

= .1


2.987(1 + (8 / 3)(4 n~i / a ) 2 )(n~i / a )

(

2 1 + ( 4 n~i / a )

)

11/ 6

+ .2

ISSUE:011

(

1294 . n~i / a 2 1 + (2 n~i / a )

FINAL

)

5/ 6

F2

where

n x Li n~i = U and i is either v or w. The five additional parameters a, . 1, . 2, F1 and F2 are defined as follows:

[ ] F = 1 + 2.88 exp[ 0.218( n~ / a ) ] F1 = 1 + 0.455 exp 0.76(n~u / a )

0.8

0. 9

2

i

. 2 = 1 .1 .1 = 2.357a 0.761

a = 0535 . + 2.76(0138 . A) 0.68 where

A = 0115 . [1 + 0.315(1 z / h) 6 ]2 / 3 Here z is the height above ground, and h is the boundary layer height obtained from:

h = u * / (6 f )

f = 2 sin(

)

(the Coriolis parameter: the earth, and

u = ( 0.4U 34.5 f . z ) / ln( z / z0 ) z0 = surface roughness length

is the angular speed of rotation of

is the latitude)

*

The turbulence intensities of the three components of turbulence are also defined for the same choice of z, z0, U and , as follows:

- = 1 6 f . z / u* p = -16 ,u =

7.5-( 0538 . + 0.09 ln( z / z0 )) p u * 1 + 0156 . ln u * / f . z0

(

Iu = , u / U

I v = I u 1 0.22 cos4

)

(the longitudinal turbulence intensity)

z 2h

(the lateral turbulence intensity)

56 of 82



I w = I u 1 0.45 cos4

z 2h

ISSUE:011

FINAL

(the vertical turbulence intensity)

although these may be changed by the user for any particular simulation. The nine turbulence length scales are also defined, as follows:

x

Lu =

(

)

3

A1.5 , u / u * z 2.5Kz1.5 (1 z / h) 2 (1 + 5.75z / h)

( = 05 . L (1

( 0.68 exp(

y

Lu = 0.5x Lu 1 0.46 exp 35( z / h) 1.7

z

Lu

x

u

)) ))

35( z / h)1.7

Lv = 05 . x Lu (, v / , u ) 3 x Lw = 0.5x Lu ( , w / , u ) 3 y Lv =2 y Lu (, v / , u ) 3 z Lv = zLu ( , v / , u ) 3 y Lw = yLu (, w / , u ) 3 z Lw =2 z Lu (, w / , u ) 3 x

where

[

K z = 019 . ( 019 . K0 ) exp B( z / h)

N

]

K0 = 0.39 / R 0.11 B = 24 R 0.155 N = 124 . R 0.008 u* R= f . z0 Associated with the von Karman spectral equations is an analytical expression for the cross-correlation of wind speed fluctuations at locations separated in both space and time, derived assuming Taylor's frozen turbulence hypothesis [7.6]. Accordingly for the longitudinal component at points separated by a distance "r perpendicular to the wind direction, the coherence Cu ("r,n), defined as the magnitude of the cross-spectrum divided by the auto-spectrum, is:

Cu ( "r , n) = 0.994( A5/6 (- u )

1 2

- u 5/3 A1/ 6 (- u ))

Here Aj(x) = xj Kj(x) where K is a fractional order modified Bessel function, and

-i =

0.747 "r 2 Li

2

2 n"r + c U

2

for i = u

57 of 82



ISSUE:011

FINAL

The local length scale Lu("r,n) is defined by: ( y L u "y) 2 + ( z L u "z) 2

L u ("r, n ) =

"y 2 + "z 2

while

c = max(10 . , with

16 . ( "r / 2 Lu ) 0.13

- 0b

b = 0.35( "r / 2 Lu )

)

0.2

and

-0 =

0.747 "r 2 Lu

2

2 n"r + U

2

"y and "z are the lateral and vertical components of the separation "r, and yLu and zLu are the lateral and vertical length scales for the longitudinal component of turbulence. For the lateral and vertical components, the corresponding equations are:

Ci ( "r , n) =

[

0.597 4.781# i2 A5/ 6 (- i ) 2 2.869# i 1

A11/ 6 (- i )

]

for i = v,w

where -i is defined as above for i = v, w, and

#i =

- i 2 Li ( "r , n) "

In this case the local length scales are given by: L v ("r, n ) =

( y L v "y / 2) 2 +( z L v "z) 2 "y 2 + "z 2

and L w ( "r, n ) =

( y L w "y) 2 +( z L w "z / 2) 2 "y 2 + "z 2

The three turbulence components are assumed to be independent of one another. This is a reasonable assumption, although in practice Reynolds stresses may result in a small correlation between the longitudinal and vertical components near to the ground.

58 of 82



ISSUE:011

FINAL

7.5.3 The Kaimal model The autospectral density for the longitudinal component of turbulence, according to the Kaimal model, is:

nSuu (n)

, 2u

=

4n~u (1 + 6.0n~u )5/ 3

where Suu is the auto-spectrum of wind speed variation, n is the frequency of variation, , u is the standard deviation of wind speed variation and n~u is a non-dimensional frequency parameter given by:

n L1 n~u = U Here L1 = 2.329 xLu where xLu is the length scale of longitudinal turbulence, and U is the mean wind speed as before. A simpler coherence model is used in conjunction with the Kaimal model. With the same notation as in Section 7.5.1, the coherence is given by

C ( "r , n) = exp

8.8"r

n U

2

012 . + L( "r , n)

2

7.5.4 Compatibility with IEC 1400-1 The turbulence model defined in the IEC standard 1400-1 [7.1] assumes isotropic turbulence. In this case, we have xLu = 2 yLu = 2 zLu and n-2/3 high frequency modification to the local length scale L("r,n) is not applicable. The above relationships are in fact equivalent to the IEC 1400-1 definition for frequencies below 0.008 Hz provided xLu = 2 yLu = 2 zLu. 7.5.5 Using 3d turbulent wind fields in simulations The following points should be noted when using these turbulent wind fields for wind turbine simulations: • The length of the wind field, Lwind, must be sufficient for the simulation to be carried out. For a simulation of T seconds at a mean wind speed of U m/s, Lwind must be at least UT + D metres where D is the turbine diameter (the extra diameter is needed in case the turbine is yawed with respect to the mean wind direction). • The width and height of the wind field must evidently be sufficient to envelope the whole rotor, i.e. at least equal to the rotor diameter. • A grid of about 7x7 points to cover the rotor plane is generally sufficient. Clearly the number of points required to achieve suitable resolution of the spatial turbulent variations will depend on the ratio of the turbulence length scales used to the rotor diameter. • If a simulation uses only a part of a turbulent time history, the mean wind speed and turbulence intensity for that part of the time history may not be the same as for the whole time history, and therefore may not match the mean wind speed and turbulence intensity 59 of 82



ISSUE:011

FINAL

which was specified for the simulation since this assumes that the whole time history will be used. Note also that a turbulence history of length Lwind with an along-wind step size of "L would be calculated at Lwind/"L points in the along-wind direction. This must be a power of two for efficient calculation, since Fast Fourier Transform techniques can then be used. If it is not a power of two, then the spacing "L will automatically be decreased to make Lwind/"L a power of two. • Different time histories with the same turbulence characteristics can be generated by changing the random number seed. • A sinusoidal half- or full-wave wind direction transient as described in Section 7.4.3 may be superimposed on the turbulent wind field. This is intended for use with turbulent wind fields when only the longitudinal component has been generated, to ensure that some yaw error occurs during the simulation. Using all three components of turbulence should give a more realistic variation of yaw error.

60 of 82



ISSUE:011

FINAL

8. MODELLING WAVES AND CURRENTS For wind turbines sited offshore, the fatigue loads and extreme loads experienced by the tower are strongly dependent on the action of waves and currents on the tower base. For fatigue load calculations in particular it is important to couple the wind and wave load calculations so that both aerodynamic and hydrodynamic damping act together to moderate tower movement. For fatigue load calculations, Bladed creates a series of irregular waves based on linear Airy theory. The amplitude and frequency content of these waves are specified by the user in terms of a power spectral density function. This may be either: • the standard JONSWAP / Pierson-Moskowitz function, or • a user-defined function. For extreme load calculations, a regular wave train may be defined. The kinematics of this wave are calculated using stream function theory.

8.1 Tower and Foundation Model Offshore wind turbines are most likely to be installed in relatively sheltered inshore conditions, where the sea depth is in the range 5m to 25m. Bladed assumes that the tower is fixed to the sea bed as a simple monopile as shown in Figure 8.1 below. The tower may be defined over the full depth (Figure 8.1a) or above a rigid base (Figure 8.1b). In both cases, the turbine structure is regarded as being transparent to the waves, implying that both tower and base are slender in comparison to the wavelength.

a) Simple Monopile

b) Monopile with narrow base

Figure 8.1: Assumed base structures As for onshore cases, the tower is assumed to have a circular cross-section and may be tapered. Foundation translational and rotational stiffnesses may also be specified.

61 of 82



ISSUE:011

FINAL

8.2 Wave Spectra To create an irregular wave train for fatigue load calculations the user must specify a suitable wave spectral formula S! ( f ) . This function will depend on the location of the turbine being modelled and the prevailing meteorological and oceanographic conditions. Bladed allows the wave spectrum to be specified in one of two ways: as a JONSWAP / Pierson-Moskowitz spectrum or as a user-defined look-up table. 8.2.1 JONSWAP / Pierson-Moskowitz Spectrum There are several different versions of the JONSWAP formula. The version used is based on an expression by Goda [8.1]. S! ( f ) = $

2 2 H s Tp

5

f fp

f 125 . fp

exp

4

#.

where f is the wave frequency (in Hz), H s is the significant wave height, Tp is the peak spectral period, f p = 1 Tp , # is the JONSWAP peakedness parameter,

$2 =

0.0624

. = exp 0.5

and

0.185 1.9 + #

0.230 + 0.0336#

, = 0.07 for f

f fp

2

1

,

fp

, = 0.09 for f > f p The Pierson-Moskowitz spectral density function may be regarded as a special case of the JONSWAP spectrum with # = 1.0 : S! ( f ) =

0.3123H s2 Tp

f fp

5

exp

f 1.25 fp

4

If the JONSWAP / Pierson-Moskowitz option is selected, the user is required to enter values for H s , Tp and # . 8.2.2 User-defined Spectrum A user-defined spectrum may be entered in the form of a look-up table. Up to 100 pairs of S! ( f ) and f may be entered. The values of S! ( f ) at the lowest and highest frequencies entered should be zero. At frequencies between the specified values of f , values of S! ( f ) are linearly interpolated. 62 of 82



ISSUE:011

FINAL

8.3 Upper Frequency Limit Waves which have wavelengths much smaller than the diameter of the tower do not contribute to the net force because regions of positive and negative velocity are experienced by the tower at the same time. Applied forces are calculated from the wave particle kinematics at the tower centreline only and so the calculations specifically exclude these high frequency components. The frequency cut-off is based on experimental work by Hogben and Standing [8.2] which show that the applied force on a cylinder falls off rapidly when the wave number exceeds 1 / radius. Therefore: S! ( f ) = 0 for k >

1 radius

The radius is taken as the minimum tower radius between the sea bed and a height of 3 standard deviations of the wave elevation above the mean water level. At any instant, the wave elevation has a probability of 99.85% of being within this range.

8.4 Wave Particle Kinematics For both the fatigue and extreme wave load calculations, wave particle kinematics are based on linear Airy theory. The following equations describe the wave particle velocity vector u w = u wx , u wy , u wz , the corresponding acceleration vector u& w = u& wx , u& wy , u& wz , the hydrodynamic component of the pressure p and the water surface elevation ! for a regular wave of height H and period T at the point ( x , y , z ) :

(

(

)

u wx =

H cos µ w cosh k (d + z ) cos($ 2 sinh( kd )

t)

u wy =

H sin µ w cosh k (d + z ) cos($ 2 sinh( kd )

t)

u wz =

H sinh k (d + z ) sin($ 2 sinh( kd )

u& wx = u& wy = u& wz =

[

]

[

[

]

]

2

t)

H cos µ w cosh k (d + z ) sin($ 2 sinh( kd )

[

2

]

H sin µ w cosh k (d + z ) sin($ 2 sinh( kd )

[

2

]

H sinh k (d + z ) cos($ 2 sinh( kd )

[

]

p=

gH cosh k (d + z ) cos($ 2 cosh( kd )

!=

H cos($ 2

[

]

)

t) t)

t) t)

t)

where = 2 f is the angular wave frequency, f is the wave frequency, t is time, d is the water depth (assumed to be constant), is the water density, g is the acceleration due to gravity and 63 of 82


$ = kx cos µ w


ISSUE:011

FINAL

ky sin µ w

where µ w is the direction from which waves arrive at the tower. The wave number k is found as the solution to the dispersion relation: 2

= gk tanh kd

The co-ordinate system used for the wave and current calculations is a right-handed Cartesian system in which the xy plane is horizontal with the x-axis pointing to the North, the y-axis pointing to the West and the z-axis pointing vertically upwards. The origin of the co-ordinate system lies where the tower centre line intersects the mean water level. Angles are defined relative to the x-axis (North) and increase positively toward the East. For the calculation of regular extreme waves, the above equations are used directly to calculate the wave particle kinematics at each submerged tower station. For fatigue load calculations, however, it is necessary to calculate an irregular (i.e. random, non-repeating) series of waves. This is achieved using the filtered white noise ‘shift register’ procedure described in section 8.6 below.

8.5 Wheeler Stretching A limitation of Airy theory is that it only defines wave particle kinematics up to the mean water level (z = 0). The theory can be extended above the mean water surface, up to the level of the wave crest, by using the Airy formulae with positive values of z. However this approach causes calculation difficulties and is known to over-estimate particle velocities and accelerations in the crest region and to underestimate velocities and accelerations in the troughs. To avoid these difficulties, Bladed uses Wheeler stretching [8.3] to take account of the forces acting between mean water level and the instantaneous free surface. Experimental results by Gudmestad [8.4] indicate that Wheeler stretching provides satisfactory estimates of particle kinematics in the free surface zone in deep water. Wheeler stretching assumes that particle motions calculated using Airy theory at the mean water level should actually be applied at the instantaneous free surface. Airy particle motions calculated at locations between the sea-bed and mean free surface are shifted vertically to new locations in proportion to their height above the sea bed. Airy wave particle kinematics calculated at a vertical location z are therefore applied to a new location z defined by: z =

d + ! (t ) d

+ ! (t )

where ! (t ) is the surface elevation above the location in question.

8.6 Simulation of Irregular Waves During an offshore simulation in which waves are specified, the following records are synthesised:

64 of 82



ISSUE:011

FINAL

• Wave elevation at the tower centre-line, • Wave particle velocities, accelerations and dynamic pressures at various points on the structure, • Wave forces on the submerged tower. For irregular waves, these records are created by the digital filtering of pseudo-random white noise. A single white noise record is used, together with a different filter for each time history to be generated. Because each filter introduces the correct amplitude variation and phase shift, the resulting output time histories display the correct amplitude and phase relationships to each other. Unlike the generation of turbulent wind records, which are generated and written to a file before running the simulation, wave data are generated as the simulation proceeds. The relationship between the parameter of interest (i.e. the wave particle velocity at the first tower station, the particle acceleration at the sea bed etc.) and the water surface elevation is defined in terms of a complex function of the wave frequency known as a Response Amplitude Operator (RAO). It is represented as a complex number of the form: RAOr = Rr e i0 r

The filters used to process the pseudo-random white noise are Finite Impulse Response (FIR) filters and are defined in terms of their frequency transforms. The transformed filter for response r is given by: zm,r = Rr ( f m )

z

m ,r

S ( f m )"f 4N

[

]

exp i0 r ( f m )

= z m ,r

where f m = m"f "f = f max N

and m is in the range 0

m

N.

The filter weights are then obtained as the transform of the expression: wn ,r =

N

1

m = N +1

z m,r exp

imn N

Having generated the filter functions for each parameter at each required location, time histories are generated using a shift-register technique. Firstly an N-element array of normally-distributed random numbers is created. The random numbers are generated by converting the output of a simple random number generator to a normally distributed deviate with zero mean and unit variance using the Box-Muller method. For each filter function in turn, the N filter weights are multiplied by the values of the equivalent elements in the random number array and the N products are then summed to give the value of the property at one particular instant in time. To calculate the value of the property at the next time step, the elements of the random number array are ‘shifted’ one place higher in the array, a new

65 of 82



ISSUE:011

FINAL

random number is introduced at element 1 and the multiplication and summation process is repeated.

8.7 Simulation of Regular Waves If “Extreme Deterministic Waves” are selected from the Waves panel, the wave kinematics are calculated using stream function theory. This method is more accurate than linear wave theory in cases where the wave height is a significant proportion of the mean water depth. The method may even be used to model waves with amplitudes close to the breaking wave limit. In cases when currents are specified in addition to regular waves (see section 8.8), the wave calculation takes proper account of the influence of the current profile on the wave kinematics. The non-linear regular wave calculations within Bladed are based on original coding by Chaplin [8.5]. Regardless of whether current components are specified, Bladed first solves the wave equation using stream function theory for the case of no currents. Stream function theory was first applied to wave modelling by Dean [8.6 & 8.7] who developed the following form of stream function: X1 z+ T

( x, z ) =

where

1

N n=2

X n sinh(nk ( z + d )) cos(nkx)

X 1 = wavelength X n +1 =

n

and N is the order of the stream function solution. as defined above satisfies the requirements that (i) the shape of the The stream function free surface is compatible with the motion of the water just below it (the Kinematic Free Surface Boundary Condition), (ii) the flow is periodic, and (iii) the flow is compatible with the presence of a horizontal sea bed at the specified depth. The values of X n are determined by a least-squares method to satisfy the additional requirements that (i) the pressure on the free surface is uniform (the Dynamic Free Surface Boundary Condition), and (ii) the required wave height is obtained. As implemented within Bladed, the order of the solution, N, is automatically chosen based on the input values of wave height, period and mean water depth. Once the stream function solution has been obtained, the horizontal and vertical velocities (in the absence of a current) are calculated using the relations:

u=

& &y

and

v=

& &x

and the dynamic pressure is calculated using Bernoulli’s equation. In cases where a current profile is specified the flow is in general rotational and the wave solution must be modified. The method used follows the approach developed by Dalrymple [8.8 & 8.9] and is based on coding by Chaplin [8.5]. It is assumed that the relationship

66 of 82



ISSUE:011

FINAL

and stream function is the same for the combined flow as for the between the vorticity undisturbed current when viewed from a reference frame moving at such a speed that total flow rate is the same as that in the x,y frame. This requirement can be stated mathematically as:

&2 &2 + = &x 2 &y 2

= f( )

In Bladed, the stream function is computed at discrete points in the x,y plane using a finite-difference calculation scheme. The most difficult feature of this approach is that the location of the free surface is not known in advance. A regular grid of points in the x,y plane would therefore have awkward intersections with the free surface profile, which must itself be calculated as part of the computation. To overcome this difficulty, Dubreil-Jacotin’s method is used to transform the problem from the x,y plane to the x, plane, with y as the field variable. The position of the free surface is now defined along the upper boundary of a plane. Treating x and as the independent parameters, the rectangular grid in the x, velocity components are now given by:

u=

1 &y

&y and

&

v=

&y

&x

&

The accuracy of the solution relies on a sufficiently fine mesh in the x, plane to resolve the structure of the flow and to allow the evaluation of derivatives on the boundaries of the computational domain, particularly at the free surface. For this purpose a regular grid in the x, plane is rather inefficient and therefore a stretched grid is employed which is finer near the free surface than the sea bed. After solving the finite difference relations on this plane, the flow velocities are calculated using the equations above and dynamic pressures are calculated using Bernoulli’s equation. Reference [8.5] should be consulted for further details of this method.

8.8 Current Velocities Bladed allows current velocities to be calculated based on three current profiles, either separately or in combination: • a near-surface (wind/wave generated) current: u cw • a sub-surface (tidal and thermo-saline) current: u cs • a near-shore (wind induced surf) current: u cn These three velocity vectors have the form:

67 of 82



ISSUE:011

FINAL

u cw = ucw ( z ) (cos µ cw , sin µ cw ,0) u cs = ucs ( z ) (cos µ cs , sin µ cs ,0)

u cn = ucn ( z ) (cos µ cn , sin µ cn ,0)

where µ cw , µ cs and µ cn are the directions from which the three current components arrive at the tower. Components of the calculated current velocities are then combined linearly: u c = u cw + u cs + u cn

8.8.1 Near-Surface Current The near-surface current velocity profile is of the form: u cw ( z ) = 2 ( z ) us ( z10 )

where us ( z10 ) is an input parameter, representing the mean wind speed at a height 10m above the mean water surface. 2 ( z ) is given by the formulae: z if 15m 15 2 ( z ) = 0.0 if z < 15m

2 ( z ) = 0.01 1

z

0m

8.8.2 Sub-Surface Current The sub-surface current velocity profile is of the form: ucs ( z ) =

z+d d

17

u s 0 ( z = 0)

for 0 4 z 4 d , where d is the water depth and u s0 ( z = 0) is an input parameter equal to the velocity at the sea surface. 8.8.3 Near-Shore Current The near-shore current velocity has a uniform profile, independent of depth. The design velocity at the location of the breaking wave is defined as: ucn = 2 s gH B

where g is the acceleration due to gravity, s is the beach slope and H B is the breaking wave height given by: HB =

b a 1 + d B gTB2

68 of 82


where:

[


]

a = 44 1 exp( 19 s)

[

b = 1.6 1 + exp( 19 s)

ISSUE:011

FINAL

]

d B is the water depth at the location of the breaking wave and TB is the period of this wave. For very small beach slopes H B may be estimated using the formula H B = 0.8d B .

8.9 Total Velocities and Accelerations The wave particle velocity and acceleration vectors at a particular location ( x , y , z ) in the wave field at time t, obtained from the white noise filtering procedure, are denoted by u w and u& w . The total current velocity vector at the same location is u c and the velocity and acceleration of the tower structure itself are u s and u& s . The total velocity u t and acceleration u& t of the fluid relative to the structure at this location and time are therefore: ut = uw + uc u& t = u& w u& s

us

8.10 Applied Forces Having evaluated the total particle kinematics relative to the tower, the resulting forces are calculated as the sum of two components: • Drag and inertia forces calculated using the relative motion form of Morison’s equation, • Longitudinal pressure forces. These forces are then used to calculate the tower modal forces as described in Section 3.2.2. 8.10.1 Relative Motion Form of Morison’s Equation To calculate the forces on the tower, the monopile is approximated by 10 cylindrical subelements of equal height. The forces on each sub-element, acting normal to the cylinder axis, are calculated using the ‘relative motion’ form of Morison’s equation: F = ( Cm

1)

D2 u& t + 4

D2 1 Lu& w + Cd DLut ut 4 2

where F is the normal force on a segment of cylinder of length L and diameter D, water density, Cm is the inertia coefficient and Cd is the drag coefficient.

is the

8.10.2 Longitudinal Pressure Forces on Cylindrical Elements Morison’s equation gives the force on the cylindrical element normal to the element’s axis. In situations where the monopile is strongly tapered, pressures acting longitudinally on the changing cross-sectional area may cause a significant vertical force to act on the pile.

69 of 82



ISSUE:011

FINAL

The added mass acting in the longitudinal direction is very small and so the longitudinal forces are estimated using the hydrodynamic pressure in the ambient wave field acting over the change in cross-sectional area of the tower between the top and bottom faces of each subelement. For a tower with diameter Da at the top of a sub-element and diameter Db at the bottom, the longitudinal force acting on this portion of the tower is: F=

(D

2 a

4

Db2

)p

No pressure force is included where the end of the tubular member passes through the free surface or terminates at the sea-bed.

70 of 82



ISSUE:011

FINAL

9. POST-PROCESSING Bladed includes an integrated post-processing facility which allows the results of wind turbine calculations to be processed further in various ways. The theory behind these postprocessing calculations is described in this section.

9.1 Basic statistics The following basic statistical properties of a signal are calculated: Minimum Maximum Mean

MIN(x) MAX(x) x

Standard deviation

,=

Skewness

(x

x )3 / , 3

Kurtosis

(x

x )4 / , 4

(x

x )2

9.2 Fourier harmonics, and periodic and stochastic components Wind turbine loads consist of both periodic and random or stochastic components. The periodic components of loads result from effects which vary as a function of rotor azimuth, such as gravitational loads, tower shadow, yaw misalignment, wind shear etc. The stochastic components result from the random nature of wind turbulence. In understanding the loads on a wind turbine it is often useful to separate out the periodic and stochastic parts of a load time history, and future to analyse the periodic part in terms of the harmonics of the fundamental rotational frequency. The periodic part of a signal is obtained by binning the signal against rotor azimuth. The number of azimuth bins may be specified by the user, otherwise it is calculated from the first two azimuth values in the time history. These are used to define the azimuth bin width, which is then adjusted to an exact sub-multiple of a revolution. The number of azimuth bins must be compatible with the sampling interval of the time history. If too many bins are used, it is possible for some of them to be empty, in which case the calculation will not proceed. Having obtained the periodic component of the signal, the Fourier harmonics are obtained by means of a discrete Fourier transform, after first increasing the number of bins by two to four times using linear interpolation. The stochastic component of the signal is obtained for each time point by subtracting the periodic component calculated from the azimuth at that time point. Linear interpolation is used between azimuth bins.

71 of 82



ISSUE:011

FINAL

9.3 Extreme prediction The prediction of the extreme loads which are likely to be encountered by a wind turbine during its lifetime is clearly a crucially important part of the design process. It is common practice to base the prediction of these extreme loads on deterministic load cases, in which the wind turbulence is represented in terms of discrete gusts with amplitudes and rise times as specified by design standards and certification rules. Discrete gusts can be modelled with Bladed as described in Section 7.4.3. An alternative approach, which avoids the problem of the rather arbitrary nature of these discrete gusts, is based on probabilistic techniques, with the stochastic nature of the loads due to wind turbulence represented by means of a probability distribution. Although this approach has been used for many years for the evaluation of extreme loads on buildings and similar structures, its application to wind turbine loads is relatively rare. The analysis involved in applying it to an operational wind turbine is rather more complicated since the probability distribution of the combined stochastic and deterministic load components must be considered. Any particular wind turbine loading can be expressed as y(t) = z(t) + x(t) where z and x represent the periodic and stochastic parts of the load respectively (see Section 9.2). It is generally a good approximation to assume that the stochastic part of the load is Gaussian, so its probability distribution is:

p( x ) =

1

,x 2

e

x 2 / 2, 2x

where ,x is the standard deviation of x. For such a signal, Rice [9.1] has derived the probability distribution of signal peaks as:

p$ ( x ) =

1 #2

,x 2

e

- 2 / 2 (1 # 2 )

+

-# e 2, x

-2 / 2

1 + erf

2

#2

2

where

- = x /, x # = 50 /5m

50 =

M2 M0

(the zero up-crossing or apparent frequency)

5m =

M4 M2

(the frequency of peaks)

72 of 82



ISSUE:011

FINAL

6

Mi =

f i H ( f )df

(the ith spectral moment)

0

f = frequency (Hz) H ( f ) = power spectral density (see next section for calculation details), and erf () = error function. Knowing the probability distribution of peaks for such a process, the probability distribution of extremes can then be deduced. For the extreme of the signal in a given period to be x, one peak must have this value and all other peaks in the period must have a lesser value. The probability distribution can be written N p$$ (-) = Np$ (-)(1 Q(-))

1

where 6

Q(-) = p$ (-) d- , and -

N = Number of peaks in the period. Davenport [9.2] combined this with Rice’s equation to give the following analytical expression for the probability distribution of extremes:

p$$ (-) = - e where 2

= 5 0 T e - / 2 and T = time period. The mean of this distribution is

-ext = . +

$ .

where

. = 2 ln(5 0 T ) and $ = 0.5772 (Euler’s constant). As the term 5 0 T increases, the distribution of extremes has a larger mean and becomes very narrow. For an operational wind turbine whose loads are a combination of stochastic and periodic components, Madsen et al [9.3] proposed an approach based on Davenport’s model of the stochastic signal, with the assumption that the extremes in the total signal occurred at minima and maxima of the periodic component. This allows the periodic time history to be idealised as a square waveform as follows:

73 of 82



ISSUE:011

FINAL

z t1

zmax

t2

zmean

Time t3

zmin T0

The resulting expressions for the mean and standard deviation of the extreme distributions are: For extreme maxima:

ye max = z max + , x . 1 +

, e max = , x

$ .1

6. 1

where

. 1 = 2 ln( 15 0T ) 1 =

t1 = T0 ( zmax

, 2z zmean )( z max

zmin )

while for extreme minima:

ye min = zmin , e min = , x

, x .3 +

$ .3

6. 3

where

. 3 = 2 ln( 35 0T ) 3

t = 3 = T0 ( zmean

, 2z zmin )( z max

z min )

Here , z is the standard deviation of the periodic component z. The time period T should be taken as the total time for which the condition being modelled will be experienced during the lifetime.

74 of 82



ISSUE:011

FINAL

9.4 Spectral analysis Bladed allows the calculation of auto-spectral density, cross-spectral magnitude and phase, transfer functions and coherence functions. All calculations involving spectral analysis use a Fast Fourier Transform technique with ensemble averaging. To perform the spectral analysis, the signal is divided into a number of segments of equal length, each of which contains a number of points which must be a power of 2. The segments need not be distinct, but may overlap. Each segment is then shaped by multiplying by a ‘window’ function which tapers the segment towards zero at each end. This improves the spectrum particularly at high frequencies. A choice of windowing functions is available. Optionally, each segment may have a linear trend removed before windowing, which can improve the spectral estimation at low frequencies. The final spectrum is obtained by averaging together the resulting spectra from each segment, and scaled to readjust the variance to account for the effect of the window function. The information required is therefore as follows: Number of points: the number of datapoints per segment. This must be a power of 2: if it is not, it is adjusted by the program. The maximum allowed is 4096. The larger the number of points, the better will be the frequency resolution, which may be important especially at low frequencies. However, choosing fewer points may result in a smoother spectrum because there will be more segments to average together. If in doubt, 512 is a good starting point. Percentage overlap: the overlap between the segments. This must be less than 100%. 50% is often satisfactory, although 0% may be more appropriate if a rectangular window is used. Window: a choice of five windowing functions is provided: (a) rectangular (equivalent to not using a window) (b) triangular: 1 2 f 1 (c) Hanning: (d) Hamming:

(1 cos(2 f )) / 2 0.54 0.46 cos(2 f )

(e) Welch:

1 (2 f

1) 2

where f is the fractional position along the segment (0 at the start, 1 at the end). One of the last three windows (which are all quite similar) is recommended. Trend removal is usually desirable.

9.5 Probability, peak and level crossing analysis These calculations work by binning values. The range and size of the bins to be used are calculated by the program, unless they have been supplied by the user. The probability density analysis simply bins the signal values. From the probability density function it also calculates the cumulative probability distribution. Also a Gaussian distribution is calculated for comparison, which has the same mean and standard deviation as

75 of 82



ISSUE:011

FINAL

the signal. There is an option to remove the mean of the signal: this merely moves the mean of the calculated distribution to zero. The peak analysis bins only those signal values which are turning points of the signal. Peaks and troughs are binned separately, so that the probability distribution of each can be output. For the level crossing analysis, the number of up-crossings and down-crossings are counted at each of the bin mid-points. The number of crossings per unit time in each direction is output for each bin mid-point.

9.6 Rainflow cycle counting and fatigue analysis Bladed offers the possibility of rainflow cycle counting of a stress time history and of subsequent fatigue analysis based on the cycle count data. A suitable stress time history can be generated from one or more load time histories by use of the channel combination and factoring facility provided by the code. 9.6.1 Rainflow cycle counting Rainflow cycle counting is the most generally accepted method used as the basis of fatigue analysis of structures. The key advantage of the rainflow cycle counting method is that it is able to take proper account of stress or strain reversals in the context of a stress-strain hysteresis loop. The cycle counting procedure involves the following steps: • The stress history is searched to determine the successive peaks and troughs by identification of turning points. • The successive peaks and troughs are re-ordered so that the sequence begins with the highest peak value of the stress history. • The sequence of peaks and troughs is now scanned to determine the rainflow cycles. A rainflow cycle is only recorded when the range exceeds a user specified minimum range. The purpose of this user-specified minimum range is to filter out very small cycles where this is desired. • The mean and range of each rainflow cycle is recorded. • The count of rainflow cycles is binned according to the cycle mean and range values. The distribution of bins is defined by the user who is required to specify minimum and maximum values of stress and the number of bins to be used. The output from the rainflow cycle counting analysis consists of the two-dimensional distribution of the number of cycles binned on the means and ranges of the cycles. This calculation can also be extended to generate damage equivalent loads. The user specifies one or more inverse S-N slopes m (see next section) and a frequency f (typically 1P for fixed speed machines), and an equivalent load is calculated as the amplitude of a

76 of 82



ISSUE:011

FINAL

sinusoidal load of constant frequency f which would produce the same fatigue damage as the original signal. The equivalent load is therefore given by:

1 ni Si m

1 m

i

Tf where ni is the number of cycles in stress range Si and T is the duration of the original time history. 9.6.2 Fatigue analysis As is described above, a complex stress history can be represented in terms of constituent cycles by use of the rainflow cycle counting technique. The distribution of rainflow cycles is defined in terms of the number of cycles binned against stress range and mean value. The basis of the fatigue analysis provided in Bladed is that fatigue failure is predicted to occur according to the Palmgren-Miner [9.4] linear cumulative damage law. Failure will occur when the “accumulated fatigue damage number” is equal to 1.0 as follows:: ni

1N i

= 1.0

i

where ni is the number of rainflow cycles of the ith stress range and Ni is the corresponding number of cycles to failure. The summation is defined as the accumulated damage. For rainflow cycles of stress range Si, the number of cycles to failure Ni is given by the S-N curve for the material. The user of Bladed must supply the S-N curve in one of two ways. The first possibility is that the S-N curve is provided as a log-log relationship of the form: log S =

1 log k m

1 log N m

so that: N = kS

m

The user must specify the value of m, the inverse slope of the log S against log N relationship. The user must also specify the intercept of the log-log relationship, c. The parameter k above is related to the intercept c by: k = cm

The second option is for the user to specify the S-N curve as an arbitrary function through the use of a look-up table. For a material where the mean stress has an influence on the fatigue damage accumulated, Bladed offers the option of converting each cycle range to the equivalent range assuming a zero mean stress value. (A cycle with a zero mean value has a R-ratio of -1, where R is the 77 of 82



ISSUE:011

FINAL

ratio of minimum to maximum stress.) This conversion is performed by means of a Goodman diagram and the user is required to provide the ultimate tensile strength (UTS) of the material. Following the conversion, the fatigue analysis proceeds using the Palmgren-Miner law and user specified S-N curve as described above. The output from the fatigue analysis consists of the accumulated damage due to the stress history as well as the two-dimensional distribution of the proportion of the accumulated damage binned on the means and ranges of the rainflow cycles.

9.7 Annual energy yield The annual energy yield is calculated by integrating the power curve for the turbine together with a Weibull distribution of hourly mean wind speeds. The power curve is defined at a number of discrete wind speeds, and a linear variation between these points is assumed. The Weibull distribution is defined by:

F (V ) = 1 e

V cV

k

where F is the cumulative distribution of wind speed V. Thus the probability density f(V) is given by

Vk 1 f (V ) = k e (cV ) k

V cV

k

Here k is the Weibull shape factor, and c is the scale factor. For a true Weibull distribution, these two parameters are related by the gamma function:

c = 1/ 7 1+

1 k

Unless the user supplies a value for c, its value is calculated as above. Note that if a different value is supplied, the resulting distribution will have a mean value which is different from V . The annual energy yield is calculated as cutout

E =Y

P (V ) f (V )dV cutin

where

P (V ) = power curve, i.e. electrical power as a function of wind speed, Y = the length of a year, taken as 365 days. The result is further multiplied by the availability of the turbine, which is assumed for this purpose to be uncorrelated with wind speed.

78 of 82



ISSUE:011

FINAL

Frequently a steady state power curve is used, combined with a Weibull distribution of hourly mean wind speeds. For a more accurate calculation, it is desirable to use a dynamically calculated power curve given as the average power from a series of simulations based on a model of the turbulent wind field. It is common practice to use 10-minute simulations to capture the effects of turbine dynamics and wind turbulence. Strictly speaking, the appropriate Weibull distribution to use in this case would be one representing the distribution of 10-minute mean wind speeds in a year. This will typically have a slightly smaller shape factor than that for hourly means.

9.8 Ultimate loads The ultimate loads calculation, which is often required for certification calculations, is simple in concept: the results of a load case simulation are analysed to find the times at which each of a number of specified loads reaches its maximum and minimum values. The simultaneous values of all the loads at each of those instants is reported. A further calculation named ‘ultimate load cases’ further analyses the results of a number of ultimate loads calculations for different groups of load cases, and generates a histogram showing the load cases in which the maximum and minimum values of each load occurred within each group.

9.9 Flicker The Flicker calculation generates short-term flicker severity values (Pst), either from a voltage time history, or from time histories of active and reactive power. Such time histories are available from simulations with the full electrical model of the fixed speed induction generator, and also with the variable speed generator model. The flicker severity is a measure of the annoyance created by voltage variations through perception of the resulting flicker of incandescent lights. The calculation of flicker from a voltage time history is defined in [9.5]. An algorithm conforming to this standard is incorporated into the Bladed post-processor. It can also calculate flicker from a time history of active and reactive power. In this case a voltage time history is calculated first, and this can be calculated for any given network impedance to which the wind turbine might be connected. In fact the flicker for several different network impedances can be calculated in a single calculation. The network impedances are entered as a set of short circuit power levels and network angles, the network angle being arctan(X/R), where X and R are the network reactance and resistance respectively. The voltage is calculated as the solution of the following equation: U4 + U2(2{QX - PR} - U02) + (QX - PR)2 + (PX + QR)2 = 0 where U0 is the voltage at the infinite busbar, and P and Q are the active and reactive power respectively.

79 of 82



ISSUE:011

FINAL

10. REFERENCES 2.1

Glauert H, “An aerodynamic theory of the airscrew”, Reports and memoranda, AE. 43, No 786, January 1922

2.2

Prandtl L and Tietjens O G, “Applied Hydro and Aeromechanics”, Dover Publications, 1957

2.3

Pitt D M and Peters D A, “Theoretical prediction of dynamic inflow derivatives”, Vertica, Vol. 5, No 1, 1981

2.4

Gaonkar G H, Sastry VV,, Reddy T S R, Nagabhushanam J and Peters D A, “The use of actuator disc dynamic inflow for helicopter flap lag stability”, 8th European Rotorcraft Forum, France, Sept. 1982

2.5

Tuckerman L B, “Inertia factors of ellipsoids for use in airship design”, NACA Report No 210, 1925

2.6

Rasmussen F R, Petersen S M, Larsen G, Kretz A and Andersen P D, “Investigations of aerodynamics, structural dynamics and fatigue on Danwin 180 kW”, Risø M-2727, June 1988

2.7

Snel H, Houwink R, Bosschers J, Piers W J and van Bussel G J W, “Sectional prediction of 3D effects for stalled flow on rotating blades and comparison with measurements”, EWEC ‘93, Travemunde, March 1993

2.8

Leishman J G and Beddoes T S, “A semi-empirical model for dynamic stall”, Journal of the American Helicopter Society, July 1989

2.9

Harris A, “The role of unsteady aerodynamics in vertical axis wind turbines”, Recent developments in the aerodynamics of wind turbines, BWEA workshop, University of Nottingham, February 1990

2.10 Galbraith R A McD, Niven A J and Coton F N, “Aspects of unsteady aerodynamics of wind turbines”, Recent developments in the aerodynamics of wind turbines, BWEA workshop, University of Nottingham, February 1990 4.1

Ahmed-Zaid S and Taleb M, Structural modelling of small and large induction generators using integral manifolds, IEEE trans. Energy Conversion 6, 3, September 1991.

4.2

Park R H, Two-reaction theory of synchronous machines, Trans AIEE 48, 1929.

4.3

Bossanyi E A, Investigation of torque control using a variable slip induction generator, ETSU WN 6018, ETSU, 1991.

4.4

Pedersen T K, Semi-variable speed operation - a compromise? Wind Energy Conversion 1995, 17th BWEA Conference (Warwick), Mechanical Engineering Publications Ltd.

80 of 82



ISSUE:011

FINAL

7.1

IEC 1400-1, Wind turbine generator systems - Part 1: Safety requirements, First edition, 1994-12.

7.2

Powles S R J, “The effects of tower shadow on the dynamics of a HAWT”, Wind Engineering, 7, 1, 1983.

7.3

Veers P S, “Three dimensional wind simulation”, SAND88 - 0152, Sandia National Laboratories, March 1988.

7.4

Engineering Sciences Data Unit, “Characteristics of atmospheric turbulence near the ground. Part II: Single point data for strong winds”, ESDU 74031, 1974.

7.5

Engineering Sciences Data Unit, “Characteristics of atmospheric turbulence near the ground. Part II: Single point data for strong winds (neutral atmosphere)”, ESDU 85020, 1985 (amended 1993).

7.6

Engineering Sciences Data Unit, “Characteristics of atmospheric turbulence near the ground. Part III: Variations in space and time for strong winds (neutral atmosphere)”, ESDU 86010, 1986 (amended 1991).

7.7

IEC 1400-1, Wind turbine generator systems - Part 1: Safety requirements, Second edition, 1997.

7.8

Ainslie J F, “Development of an eddy viscosity model for wind turbine wakes”, Proceedings of 7th BWEA Wind Energy Conference, Oxford 1985.

7.9

Ainslie J F, “Development of an Eddy Viscosity model of a Wind Turbine Wake”, CERL Memorandum TPRD/L/AP/0081/M83, 1983.

7.10 H Tennekes and J Lumley, “A first course in turbulence”, MIT Press, 1980. 7.11 L Prandtl, “Bemerkungen zur Theorie der freien Turbulenz”, ZAMM, 22(5), 1942. 7.12 Ainslie J F, “Calculating the flowfield in the wake of wind turbines”, Journal of Wind Engineering and Industrial Aerodynamics, Vol 27, 1988. 7.13 Taylor G J, “Wake Measurements on the Nibe Wind Turbines in Denmark”, National Power, ETSU WN 5020, 1990. 7.14 Quarton D C and Ainslie J F, “Turbulence in Wind Turbine Wakes”, Wind Engineering, Vol. 14 No. 1, 1990. 7.15 U Hassan, “A Wind Tunnel Investigation of the Wake Structure within Small Wind Turbine Farms”, Department of Energy, E/5A/CON/5113/1890, 1992. 7.16 Vermeulen P and Builtjes P, “Mathematical Modelling of Wake Interaction in Wind Turbine Arrays, Part 1”, report TNO 81-01473, 1981. 7.17 Vermeulen P and Vijge J, “Mathematical Modelling of Wake Interaction in Wind Turbine Arrays, Part2”, report TNO 81-02834, 1981.

81 of 82



ISSUE:011

FINAL

8.1

Goda Y, “A Review on Statistical Interpolation of Wave Data”, Report of the Port and Harbour Research Institute, Vol. 18, No. 1, March 1979.

8.2

Hogben N and Standing R, “Experience in Computing Wave Loads on Large Bodies”, OTC 2189, Offshore Technology Conference, Houston, 1975.

8.3

Wheeler J D, “Method for Calculating Forces Produced by Irregular Waves”, J. Petr. Techn., pp.359-367, March 1970.

8.4

Gudmestad O T, “Measured and Predicted Deep Water Wave Kinematics in Regular and Irregular Seas”, Marine Structures, Vol. 6, pp.1-73, 1993.

8.5

Chaplin J R, “The Computation of Non-Linear Waves on a Current of Arbitrary NonUniform Profile,” Den Report OTH 90 326, HMSO, 1990.

8.6

Dean R G, “Stream Function Representation of Nonlinear Ocean Waves,” Journal of Geophysical Research, Vol.70, No. 18, Sept. 1965.

8.7

Dean R G, “Stream Function Wave Theory: Validity and Application,” Proceedings of the Santa Barbara Specialty Conference, Ch. 12, Oct. 1965.

8.8

Dalrymple R A, “A Finite Amplitude Wave on a Linear Shear Current,” Journal of Geophysical Research, Vol. 79, pp. 4498-4504, 1974.

8.9

Dalrymple R A, “A Numerical Model for Periodic Finite Amplitude Waves on a Rotational Fluid,” Journal of Computational Physics, Vol. 24, pp. 29-42, 1977.

9.1

Rice S O, “Mathematical analysis of random noise”, Selected papers on noise and stochastic processes, ed. N Wax, 1959.

9.2

Davenport A G, “Note on the distribution of the largest value of a random function with application to gust loading”, Proc. Inst. Civil Eng. 28, pp187-196, 1964.

9.3

Madsen P H, Frandsen S, Holley W E, Hansen J C, “Dynamics and fatigue damage of wind turbine rotors during steady operation”, Risø report R-512, 1984.

9.4

Miner M A, “Cumulative damage in fatigue”, Transactions of the American Society of mechanical Engineers, Vol. 67, A159-A164, 1945.

9.5

Flickermeter functional and design specification, BSEN60868, 1993, and evaluation of flicker severity, BSEN60868-0, 1993, equivalent to IEC 868-0, 1991.

82 of 82

GH Bladed Theory Manual

Recommend Documents