INVESTIGATION OF THE YAWED OPERATION OF WIND TTURBINES BY MEANS OF A VORTEX PARTICLE METHOD

S.G.Voutsinas, M.A.Belessis and K.G.Rados

National Technical University of Athens,

Department of Mechanical Engineering, Fluids Section,

P.O. Box 64070, 15710 Zografou, Athens, Greece


CONTENT

1. SUMMARY
2. DYNAMICS OF WIND TURBINES IN YAW AND MODELLING
3. THE NUMERICAL METHOD
3.1. GENUVP
• FLOW CHART OF GENUVP
• THE CALCULATIONS OF LOADS
3.2. THE AEROLASTIC COUPLING
4. RESULTS
5. CONCLUSIONS
ACKNOWLEDGEMENTS
REFERENCES
FIGURES

1. SUMMARY

A fully three-dimensional non-linear aeroelastic numerical investigation of the response of horizontal axis wind turbines during yawed operation was carried out. The numerical tool used, consists of a time-marching method based on the coupling of an unsteady free-wake vortex particle model and a 3D beam-type structural model. The investigation led to a complete data base of numerical results concerning the Tjareborg wind turbine for which extensive full scale measurements of very good quality exist. Among the points that were given particular attention are: the effect of the root vortex, the coupling with the shear of the inflow and the tower effect on the dynamics of the blades. Herein the most significan results are presented and discussed.


2. DYNAMICS OF WIND TURBINES IN YAW AND MODELLING

The operation of wind turbines in yaw constitutes a special case of dynamic excitation that consists of a combination of different mechanisms. As a whole this combination leads to a non-linear response that is not fully understood basically because of its close connection to three-dimensional effects on one hand and dynamic stall on the other. Due to its frequent appearance especially in cases of non-flat terrain applications, yawed operation is of particular importance from the design point of view. It certainly affects the performance (energy production) as well as the availability of the wind turbines.

As in all physical problems, the most appropriate methodology for a better understanding, is a combination of measurements and theoretical calculations. Because the corresponding experiments are difficult to realize, the number of available complete data bases is small. Even smaller is the number of campaigns that someone could rely on during his theoretical/numerical investigations. Fortunately, in Denmark an extensive campaign of high quality full scale measurements taken from the 2MW wind turbine at Tjaereborg. This campaign covered a wide range of cases with dynamic inflow effects. Some of them were made available to us. Based on these data, we carried out our numerical study (Ref 1).

In most cases the existing theoretical work on the subject, concerned mainly the evaluation and improvement of engineering tools of rather simplified structure such as the blade element model (Ref 2). The reason is clear: In practice only aeroelastic calculations leading to load spectra of wide range make sence. This type of calculations are extremely time consuming so any complication in the modeling would render the work more a research than a design proccedure. More elaborate theoretical models however do exist (Ref 3-7). They were developed rather recently as tools forapplied research. In their basic features, they are all variants of the panel or boundary element method, differing the one from the other in what concerns the modeling of the wakes of the blades. Following closely the developments made in the past for helicopters in forward flight, the aerodynamic models for wind turbines are classified into prescribed and free-wake models. No doubt that free-wake models are more suitable to give a detailed understanding of the different excitation mechanisms related to yaw. This is due to the complexity of the case of yawed operation. In particular, as the blade rotates fluctuations of the local angle will result load vibrations that can be significant especially if the incidence enters the stalling envelope of the corresponding airfoil. These fluctuations of the incidence are due to: (a) the lateral component of the inflow velocity which has an alternating effect on the magnitude of the relative with respect to the blade velocity, (b) the shear of the inflow that acts on the axial velocity component, (c) the tower that decreases again the axial velocity relative to the blade when it passes in front of it and (d) the structural flexibility that contributes to the velocity diagramm the rate of deformation as an additional velocity term. These fluctuations of the velocity result a time variation of the circulation around the airfoil considered. In its turn this variation will lead to extra vortex emission in the wake that will feed back to the velocity diagramm different induced velocity. Moreover if the amount of the shed vorticity varies, also the velocity field in the wake will vary which means that geometry of the wake becomes important. Several important cases could be mentioned: (a) Case of uniform inflow, no tower effects, no flexibility effects, (b) Case of sheared inflow, no tower effects, no flexibility effects, (c) Case of uniform inflow, with tower effects, no flexibility effects, (d) Case of sheared inflow, tower effects included, no flexibility effects, and (e) Case of sheared inflow, tower effects and flexibility effects included.

Based on the above remarks, it is easy to see that during yawed operation strong aerodynamic and aeroelastic interactions take place that justify the choice of elaborate models and in particular free-wake models. In this sub-class of free-wake models, we can distinguish the potential approach that considers the wake as a surface of discontinuity (Ref 6), the vortex-filament models that represent the wake by a set of thin vortex filaments (Ref 4,5) and the vortex particle model that represents the wake as an ensemble of blobs that carry vorticity (Ref 7). The model that was used in the present work is a later version of the GENUVP code also used and outlined in Ref 7. The difference consists of the aeroelastic coupling with a beam type structural model and a three-dimensional wind simulator. It corresponds to a computational environment that supports the research carried out on Aerodynamics at NTUA.

In brief GENUVP is defined along the following guidelines:

According to Helmholtz's decomposition theorem (Ref 8) the velocity field is made up of an irrotational part representing the disturbance due to the presence of the solid boundaries and a rotational part representing the disturbance due to the wakes. In order to determine the irrotational part a Neumann boundary value problem for the Laplace equation is solved. On the other hand the rotational part of the velocity field is determined directly as convolution of the vorticity contained in the wakes. As the fluid is assumed inviscid, there is no link between the irrotational and the rotational part of the velocity field. In order to bring these parts into contact, a link can be defined that models the vorticity emission process observed in real flows. Mathematically this link is based on the Kutta condition (more exactly on an appropriate formulation of Kelvin's theorem). The fulfilment of the Kutta condition, permits to define quantitatively the conversion of the bound-vorticity into freevorticity. This mechanism constitutes the inviscid analog of the vorticity production process already existing in all the conventional aerodynamic models.

In the discrete formulation, singularity distributions are defined on all solid boundaries (sources on the nacelle and tower, dipoles on the blades) in order to generate the potential part of the flow, whereas the wakes of the blades are approximated by vortex particles that are continuously generated along the edges of the blades. Next these vortex particles are convected downstream by the flow. This is performed by solving the vorticity transport equations by means of a Lagrangian procedure.

In the sequel, the most important findings of this campaign are presented together with a concise presentation of the model. As a general conclusion one could say that the comparisons between predictions and measurements on one hand showed the robustness of the model and on the other hand its ability to reproduce in a consistent way the reality.


3. THE NUMERICAL METHOD


3.1 GENUVP: A Free-Wake Aerodynamic Model Based on Vortex Particle Approximations

GENUVP consists of a GENeralized Unsteady Vortex Particle computational tool that can perform unsteady calculations around a configuration consisted of NB threedimensional bodies Bk, k=1...N_B that move independently the one from the other. Each component Bk can be either nonlifting (as the nacelle and the tower) or lifting (as the blades), according to its operational characteristics. For the purposes of the present work the configuration of a wind turbine is formed by the tower and the blades. Also in the same figure the convention regarding the basic operational characteristics such as the rotation speed of the turbine, the position of zero azimuthal angle and the inflow velocity.

Let D denote the flow field and S its boundary. According to Helmholtz's decomposition theorem, the velocity field

u_vec can be given the form (Ref 9):

(1)

where nabla phi represents the irrotational part of the flow,

u_inf is the inflow velocity and

u_omega is a divfree field induced by the free vorticity.

For an incompressible fluid, nabla phi is included basically to account for the non-entry boundary conditions. Because ö(•;t) is harmonic, it can be represented by source and/or dipole distributions defined on S (Ref 10). Dipole distributions are used to generate lift on solid boundaries and surface vorticity on vortex sheets whereas sources are used in order to give the displacement generated by the thickness of the bodies. In order to reduce the computational cost, the blades are considered to be lifting surfaces, i.e. of zero thickness. Therefore they are represented solely by dipole distributions. Source distributions were used only for the non-lifting tower.


Figure 1: The configuration of a Horizontal Axis Wind Turbine


In the case of a wind turbine, S will comprise the solid boundaries Sk, k=1...NB of the different components as well as the wakes SWk, k=1...NB-1 shed by the NB-1 blades (The tower is set last in the list of components. Moreover index w is used to designate quantities defined on the wakes). Thus,

(2)

For an inviscid fluid, a wake can be represented either by means of a dipole distribution, and thus included in the irrotational part of the flow field, or as a vorticity carrying (vortex) sheet in which case it will appear in u_omega (Ref 11, 8). This dual interpretation of a wake permits us to have the following hybrid formulation:

The numerical model is defined as a time marching algorithm. In fact the computations follow the gradual generation of the wakes, as part of the responce of the wind turbine to all the different external excitations. During a time step of the scheme, first the wakes that were generated up to that time, i.e the old or far parts SWk*, k=1...NB-1 , are convected downstream leaving a gap in the near vicinity of the sheding body. Next this gap is filled by the new or near part of the wakes SWkÄ (Fig 2).

In GENUVP the near-wakes are represented as strips on which dipole distributions are defined. On the contrary the far-wakes are represented as spatial distributions of vorticity. This means that as time passes, the model treats the wakes within different physical contexts. In particular, at the end of each time step as the near wake is convected it becomes part of the far-wake. In other words the dipoles o f the near-wakes are transformed into vorticity. That is why the whole scheme is termed hybrid. The passage from the "potential" interpretation to the "vortical" can be reasoned qualitatively, as being part of the deformation process that the wakes will undergo. A more rigorous justification is given bellow.


Figure 2: The hybrid modelling of a wake


The potential interpretation considers a wake Ó(t) to be a moving surface on which the potential of the flow and the tangential to Ó(t) velocity component are discontinuous (for incompressible fluids the normal to Ó(t) velocity remains continuous).

This means that Ó(t) can be represented by a dipole distribution ì(•;t). Let

x_vec=, ...

denote the Lagrangian representation of Ó(t). Then:

my=..,u_matrix=..

denotes jump on Ó(t) (Fig 3). The evolution of Ó(t) is determined by its

equation of motion as defined by the mean velocity u_vec_m

on Ó(t) and the zero pressure jump condition obtained by applying

Bernoulli's equation to the two faces of Ó(t):

(3a)

(3b)


Figure 3: Definition of surface vorticity on a vortex sheet


Condition (3b) can be recast in Lagrangian form to give the material conservation of ì(•;t), equivalent to Kelvin's theorem:

(4)

Note that the superficial material time derivative is defined with respect to the mean velocity u_vec_m.

(3b) also yields an evolution equation for the intrinsic surface vorticity g_vec=... is the unit normal to Ó and N its metric. Because this equation is similar to the well known Helmholtz evolution equation for the vorticity, g_vec is identified as part of the spatially distributed vorticity omega_vec (Ref 12). In case Ó is open, in addition to g_vec , there is also a line term defined along the boundary of Ó. Summarizing, for a vortex sheet Ó(t) , the generalized vorticity is defined as follows:

(5) where äÓ(.) and ä (.) denote the surface and line Dirac functions defined on the interior and the boundary of Ó(t) respectively and

tetra_vec denotes the unit tangential to vector (Fig 3). Note that if µ(.;t) is constant then there is no surface term and the line vorticity is reduced to a closed vortex filament along . This situation appears in the case of piecewize constant dipole distributions which is in fact the choice made in GENUVP for thin lifting surfaces.

According to the above analysis, ö(.;t) will be represented by means of: a source distribution ó(.;t) defined on the tower, bound dipole distributions µk(.;t) defined on the blades and free dipole distributions µWk(.;t) defined on the near-wakes of the blades SWkÄ. As explained in the sequel ó(.;t), µk(.;t) and µWk(.;t) constitute the unknown degrees of freedom at every time step. It is noted that also unknown is the exact geometry of the near-wakes SWkÄ.

In order to determine the unknown fields of the problem there are two types of conditions available:

(a) the kinematic conditions and more specifically the nonentry conditions on all the solid surfaces, and the conditions of material motion of the near-wakes SWkÄ(as in (3a))

(b) the dynamic conditions, i.e. the requirement of zero pressure jump throughout the vortex sheets SWk (as in (3b) or (4)).

Let u_vec_Bkdenote the rigid body velocity distribution defined on Sk and v_vec the unit normal to Sk. Then the non-entry condition for x_vec and k=1...`NB takes the form:

(6)

In the discrete formulation, ó(.;t), µk(.;t) and µWk(.;t) are defined with reference to surface grids of boundary elements. Let Ske, e=1...Ek and SWke, e=1...EWk , denote the elements of Sk and SWkÄ respectively (Ek and EWk are the number of elements that form Sk and SWkÄ ). The boundaries of these elements will be denoted by dSke and dSWke (Fig 4). Assuming that ó(.;t), µk(.;t) and µWk(.;t) are piecewize constant, the aerodynamic discrete degrees of freedom will be their intensities at the boundary elements óe(.;t), µek(.;t) and µWke(.;t) respectively.


Figure 4: The notations of the boundary element grids


The application of (6) to the centervectors xcpe(t) of the elements covering the solid boundaries Sk , gives the discrete equations wherefrom óe(.;t) and µek(.;t), e=1...Ek are determined. As regards µWke(.;t) , they are determined by applying (4) along the so called emission lines, i.e the parts of the border of the blades wherefrom the wakes originate. These lines usually include the trailing edge, the tip and the root of the blades. The leading edge is also included in the case of stall. This option was left apart throughout the work presented herein. By interpreting (4) as a Kutta condition, it results that at all boundary elements SWke of the near-wakes, µWke(.;t) must be set equal to the corresponding bound dipole intensity, i.e. the dipole intensity of the element on Sk that is adjacent to SWke. As regards the geometry of the near-wakes SWkÄ , this is defined by the Lagrangian kinematics. According to (3a), SWkÄ will be convected by the local velocity.

The last point to be considered, concerns the rotational part of the flow u_omega and in particular the process of transforming the dipole distribution of the near-wake into vorticity. Let, ùW(.;t) denote the vorticity contained in the far-wakes of all blades. Then

(7)

where denotes the support of the free vorticity

(8)

according to (5). In (8), the three summations run respectively the blades, the elements on each near-wake strip and the time steps already passed. Also note that only line-vorticity terms are included. This is due to the piecewize approximation of the dipole distributions µWk(.;t). Although it is not stricktly said, (8) suggests that the wakes are interpreted as a collection of panels that formed all the near-wake strips throughout the time interval [0,NT •Ät] (where NT is the current number of timesteps). There is enough evidence both theoretical and experimental, showing that the evolution of a vortex sheet is unstable in time. This is due to the singular character of the velocity field that it requires the calculation of Cauchy type integrals. Consequentely some kind of smoothening is required. A method to accomplish stable and consistent calculations is to introduce particle approximations for the vorticity (Ref 13,14,15). In this case,

(9)

where (Omega and Z) denote the intensities and positions of the vortex particles, J(t) the index set for the vortex particles and å(t) the cutoff function that ensure the neccessary smoothening. By selecting (Ref 16):

(10)

(Omega and Z) takes the form:

(11)

where

(12)

Thus instead of calculating the geometry of the vortex sheets and the dipole distributions they carry, we follow the evolution of the vortex particles defined by the following dynamic equations:

(13)

Equations (13) concern the evolution of the far-wakes. As the near parts still retain their character as vortex sheets their determination is different. Let U_em denote the mean velocity at a point X_em along the vorticity emission line of a lifting body. The geometry of the near part of the corresponding wake SWÄ is determined kinematically through the following relation:

(14)

where X_delta-X_em denotes the width of SWkÄ in vectorial form (Fig 5). Finally the intensity of the dipole distribution of SWkÄ is determined by means of the dynamic condition (4). According to (4) and because the model uses piecewize constant dipole distributions,ìWke(t) must be set equal to the bound dipole intensity of the element adjacent to the eWk near-wake element.

Due to the time dependent character of the problem, the wakes as well as the vortex particles they include in their far parts, will be constructed gradually. This means that vortex particles will be created as the near parts of the wakes evolve. In order to make this approach compatible with the dynamics of vortex sheets,

(Omega and Z) are defined as follows:

(15)

In the above relations, the integration covers for every vortex particle the surface of an element of the near part of the wake considered.


Figure 5: The mechanism of vortex particles generation


THE FLOW CHART OF THE GENUVP MODEL

For every time step ( t = n•Ät )

A. POTENTIAL CALCULATIONS

u=...

0. Initialize SWÄ and ÖWÄ (through the values of µW along the emission lines of the lifting bodies )

Iterative schemes for the near wake:

1. Calculate Ö (H-fulfilment of non-entry boundary conditions)

2. Calculate the emission velocities U_em along the emission lines

3.1 Correct SWkÄ

3.2 Correct matrix_PHI (H-fulfilment of the Kutta condition)

4. Check for convergence: delta_matrix_PHI

5. FIRST OUTPUT: Force and velocity calculations

B. VORTICITY CALCULATIONS

1. Create new vortex particles

2. SECOND OUTPUT: Wake structure and velocity profiles in the wake

3. Move and deform all the Vortex Particles

3.1 Calculate the velocities and deformations induced at the Vortex Particle locations

3.2 Check for Vortex-Solid surface interaction and correct accordingly

3.3 Produce the new far-wake


The calculations of loads. The calculation of the unsteady loads on the rotor blades was based on Bernoulli's equation. For a thin lifting surface this equation takes the form:

(16)

where -matrix_phi denotes the surface dipole distribution, U_m denotes the mean velocity, marix_phi=.. denotes the tangential velocity discontinuity and matrix_p the pressure jump. In the discrete problem equation (29) is applied to the control points. Due to the piecewize approximations used for the dipole distributions, marix_u is calculated through a zeroorder difference scheme [H-18]. Let matrix_p=..., denote the calculated values of the pressure jump. Integrating over the corresponding elements, the calculation of the loads is obtained as follows:

for every element,

F=..

for every body,

F=..

where [Ske] denotes the area of Ske. In addition to these loads and because of the thin-wing approximation that was made for the blades, the suction force is also added.

Due to the essentially inviscid character of the modelling adopted can only predict the "inviscid" part of the loading. The rest, i.e. the loads due to viscous effects must be superimposed. In this connection an approximate aposteriori scheme, similar to the one used by the classical strip theory was applied. This scheme is based on the calculation of the radial distribution of the effective angle of attack. Next the drag coefficient CD is estimated from given airfoil data. Finally the drag forces are integrated over radial strips and thus an account of viscous effects is obtained.

3.2 The Aeroelastic Coupling

Measurements as well as prior numerical investigations, clearly showed that flexibility effects are important in the case of yawed operation. In order to take them into account GENUVP was coupled with the structural dynamics of the rotor. To this end , the beam approximation was used to give the flexibility characteristics of the blades. The coupling between the structural and thr aeroelasticpart of the problem was realized in the time domain by performing consecutively aerodynamic and structural time steps. Two coupling interfaces were defined: the first communicates to the structure, the loadsb are calculated by the aerodynamic model, and the second communicates the the blade deformations together with their rates, to aerodynamic. Because of the different character that the aerodynamical and the structural part have , the stability restrictions for the two parts are not the same. In fact the structural part is much stiffer demanding a time step smaller than that used for the aerodynamic calculations. Consequently in order to cover the time interval corresponding to one aerodynamic time step, sevaral structural steps must be performed. In all cases reported herein, the calculations were realized by repeating the sequence: one aerodynamic step followed by several structural steps. Such a marching scheme signify the choice of an explicit time integration for the aeroelastic problem. More elaborate ad particular implicit schemes can be easily defined. They will require an iterative exchange of information between the aerodynamic and the structural part that should converge before the calculations proceed in time. There is no doubt that such a scheme would increase the computational cost a fact we took seriouslt into account because the number of runs that had to be made ( Futher details an be found in Ref 17).


Figure 6: Typical velocity diagrams for 0o and 180o azimuthal angles


4. RESULTS

The main body of the tests concern the Tjareborg wind turbine for which extensive full scale measurements exist of very good quality. The machine is upwind oriented and pitch regulated. The diameter of the 3-bladed rotor is 60 meters, rotating at 22rpm. In order to have a clear idea of the different effects, the following simulations of yawed operation were performed:

(a) uniform inflow, no tower effects, no flexibility effects,

(b) sheared inflow, no tower effects, no flexibility effects,

(c) uniform inflow, with tower effects, no flexibility effects,

(d) sheared inflow, tower effects included, no flexibility effects, and

(e) sheared inflow, tower effects and flexibility effects included.

In the sequel a representative part of the results obtained, is given.

As an introduction to the presentation of the results we thought useful to start with a brief discussion on the basics of yhe yawed operation of wind turbines. Let us assume that the rotation is clockwise and the yaw angle positive. This corresponds exactly to what Fig 1 shows. The first and most critical point is the role of the lateral component Uy of the inflow velocity. Fig 6 shows two typical velocity diagramms for azimuthal angles 0o and 180o, in the caseof uniform inflow. These two azimuth positions constitute the positions where the effect of lateral velocity due to yaw reaches its peaks. The differece is clear. At 0o azimuthal angle Uy is added to the circumferential compnent of the effective velocity UÖ . On the contrary at 180o azimuthal angle Uy is substracted. In fact this additive/subtractive role is con served respectively throughout the lower/upper half of the rotation of the blade. As a result, UÖ willl exhibit a periodic variation that results an increase of the axial force in most of the lower part and a decrease in the upper part of the rotation. Fig 7 and 8 clearly show this effect. As expected, there is a phase shift so that the two peaks do not appear at azimuthal angles 0o and 180o. The phase shift is about 45o and is due to the fact that the wake of the turbine is inclined with respect to its axis. In fact there are two interacting mechanisms of excitation with 90o lag the one from the other. The first is due to the lateral velocity at azimuthal angles 0o and 180o and the second is due to the wake with peaks at azimuthal angles 90o and 270o. The fine tuning of the 45o phase shift is determined by the details of the configuration (e.g. the pitch, the twist distribution e.t.c.). In Fig 8 the dependance on the value of the yaw angle is shown. As expected for larger yaw angles the amplitude of the fluctuation increases. For larger yaw angles, a second local fluctuation is seen around the azimuthal angle 270o as a result of a more intense interaction with the wake. In particular, most of the local fluctuation is due to the root vortex as shown in Fig 9. The effect of the root vortex is moderate and has to do with the fact that the design point of the turbine corresponds to in flow velocity equal to 9 m/s. At higer velocity the intensity of the root vortex is expected to equally increase together with its influence on the axial forces. Because of the special significance of the root vortex, its placement is given in Fig 10b.Compared to the configuration of the wake in case of zero yaw, the deviation with respect to the direction of the inflow velocity is remarked. In what concerns the calculations, this deviation is a direct consequence of the velocity deficit that changes the pitch angles along which the shedding of vorticity takes place. In our calculations, this is ensured by the freedom we give to the wake to deform in shape. It is noted that flow visualizations carried out by FFA revealed similar wake shapes.

Having theabove remarks in mind, we proceed with the discussion on the influence of the other parameters of the problem. Asummary of results is given in Fig 11, where for differenr operational conditions, the variation of the flapping moment at root due to aerodynamic forces is given for yaw angles 0o, +32o and -32o. In Fig 11a, the ideal case with no-shear , no-tower and no-flexibility effects is shown giving a simple 180o phase shift in the azimuthal angle. In Fig 11b, the blade sare treated as defrmable beams. As seen in the figure, the result of the flexibility is the decrease of the amplitudes. In Fig 11c, the effect of shear is examined. As expected for zero yaw, there is a peak at 180o azimuthal angle because the inflow in the upper part of the rotation is larger in comparision to the inflow in the lower part. When this effect is coupled with yaw then the peak is shifted , of course differently for positive and negative yaw angles respectively. Again we have the superposition of two excitation mechanisms with pahse lag of 180o. As a result not only a phase shift of the peak takes place but also a modification of the amplitude. More specifically, in the case of negative yaw, the sheared inflow acts as an amplifier for the peak whereas in the case of positive yaw it decreases it. In Fig 11d the flexibility effects are added. Again there is no significant remark to make. In Fig 11e the effect of the tower is examined. For a moreclear vision of the tower effect in Fig 12 the azimuthal variation of the included velocity in the axial direction is given for the whole span of one of the blades. As expected the tower will induce extra deficit around the zero azimuthal angle. Compared to the calculations without tower effects, it can be remarked that this effect is significat for a range about 90o. Due to the rotation of the blades ven at zero yaw, the tower effect is non-symmetric. This asymmetry becomes more pronounced when yaw is imposed. For positive yaw the minimum is shifted at an azimuthal angle of about 15o, whereas for negative yaw, thereis no real shift of the minimum but a distributed decrease ofor about 30o before the zero azimuthal angle. When flexibility effects are added to this case (Fig 11f), a distortion of the shape of the curves is noticed especially after the blade has passed the zero azimuthal angle. Fig 11g and 11hcomplete the results by adding shear to the inflow. There is no special comment to make except that the curves represent the result of the superposition of more than two excitation mechanisms on the rotor of the wind turbine.

In Fig 13 the predictions are compared against full scale measurement for four different yaw angles. The quantity considered is the flapping moment, which is important as regarda the assesssment of the structural robustness of wind turbines. In all cases the model reproduces accurately the form of the corresponding experimental curve. Also in terms of absolute values, there is a level difference which could be attributed to the difference between the processing followed for obtaining the series compared in the figure. More specifically, the inflow velocity used in the calculations is the result of the averaging over 1 min out of a 10 min time series. No spatial variations were input at the inflow plane of the turbine as happens in reality. Also approximate is the estimation of the shear. Nevertheless thre are some discrepancies that need further investigation. In three cases of positive yaw , experiments have a bump just after 270o azimuthal angle. The code onlt predicts the bump in the case of the extreme yaw of 54o. This could be due to the pronounced influence of the root vortex, which according to Fig 8 , should be small for much smaller angles. This could mean that the bumpsfor the two smaller yaw angles are of different origin. Indded the appearane of the bump for yaw angle 3o justifies such a possibility.

As regards the effect of the tower, it can be noticed that in the case of positive yaw, the low peak just after the 0o azimuthal angle, is predicted satisfactorily, even though there is a phase shift. On the contrary, at negative yaw the code fails to predict the peak as if the tower and yaw effects cancel each other.


5. CONCLUSIONS

During the investigation carried out on the yawed operation of wind turbines:

a) an advancedprediction methodology was established,

b) a detailed data-base was produced regarding the behaviour of a specific wind turbine.

According to the evauation perfformed the following conclusions can be made:

a) In general the predictions can be considered satisfactory both in quantative and qualitative sense: the model gives values comparable to measurements and the behaviour of the wind turbine as revealed by the model is similar to the experimental one.

b) However difference between measurements and calculations does exist giving hints for possible improvements. In this respect the following points are of interest:

  • The time integration scheme used for the aeroelastic coupling could be made implicit.

  • The blades of the rotor should be turned to thick wings so as to get a better flow field and also to prepare the neccesary input for viscous corrections.

  • in the case of thick blades , the coupling of the outer-inviscid region with in inner viscous layer could improve the flow predictions on the blades and consequently the prediction of the loads. This is true especially when stall occurs.


Acknowledgements

The present was was partially financed by the DG XII of the Commission of the European Union under the contract JOU2-CT92-0186 and JOU2-CT92-0113. The authors would like to thank Professor A. Panaras and Professor S. Tsaggaris for their assistance in diffusing the present work. Special thanks are addressed to Captain K. Gianniopoulos and Lieutenant A. Petropoulos of the Greek Air Force for the help they provided to K. Rados.


References

1. Oye, S.,"Tjaereborg wind turbine, First Dynamic InflowMeasurements", AFM Notak Vk-189, April 1991

2. Snel, H. and J.G., Schepers, "Investigation and Modelling of Dynamic Inflow Effects", in "European Wind Energy Conference Proccedings", March 1993, Paper H3, pp 371-375

3. Bussel, G.J.W, "The Use of the Asymptotic Acceleration Potential Method for Horizontal Axis Wind Turbine Rotor Aerodynamics", in "Wind Energy: Technology and Implementation", Proccedings of Amsterdam EWEC '91, Part I, October 1991, pp18-23

4. Gould, J. and Fiddes, S.P., "Computational Methods for the Performance Prediction of HAWTs", in "Wind Energy: Technology and Implementation", Proccedings of Amsterdam EWEC '91, Part I, October 1991, pp29-33

5. Pesmajoglou, S. and Graham, J.M.R., "Prediction of Yaw Loads on a Horizntal Axis Wind Turbine", in "European Wind Energy Conference Proccedings", March 1993, Paper H17, pp 420-423

6. Bereib. S. and Wagner, S., "The Free Wake / Hybride Wake Code ROVLM - A Tool for Aerodynamic Analysis of Wind Turbines", in "European Wind Energy Conference Proccedings", March 1993, Paper H18, pp 424-427

7. Voutsinas, S.G., Belessis, M.A. and Huberson, S., "Dynamic Inflow Effects and Vortex Particle Methods", in "European Wind Energy Conference Proccedings", March 1993, Paper H19, pp 428-431

8. Richardson, S.M. and Conrish, A.R.H., "Solution of threedimensional incompressible flow problems", J. Fluid Mech., 82, 2, 1977, pp. 309319

9. Batchelor, G.K., "An introduction to fluid dynamics", Cambridge University Press, 1967

10. Hunt, B., "The Panel Method for Subsonic AerodynamicFlows: A Survey of Mathematical Formulations and Numerical Models and an Outline of the New British Aerospace Scheme", VKI Lecture Notes, 1980

11. Hess, J.L., "Calculation of Potential Flow about Arbitary 3D Lifting Bodies", McDonnel Douglas Report, MDC J567901, 1973

12. Voutsinas, S.G., "Numerical Simulation of Three Dimensional Vortical Flows", NTUA report, 1994

13. Leonard, A., "Computing ThreeDimensional Incompressible Flows with Vortex Filaments", Ann. Rev. Fluid Mech., 17, 1985, pp. 523559

14. Hoeijmakers, "Numerical Simulation of Vortical Flows, Introduction to Vortex Dynamics", VKI Lecture Notes, 1986

15. Smith, J.H.B., "Vortex flows in Aerodynamics", Ann. Rev. Fluid Mech., 18, 1986, pp. 221242 (See also: Smith, J.H.B., "Modelling 3D Vortex Flows in Aerodynamics, Introduction to Vortex Dynamics", VKI Lecture).

16. Beale, J.T. and Majda, A., "Higher Order Accurate Vortex Methods with Explicit Velocity Kernels", J. Computational Physics, 58, 1985, pp.188208

17.Belessis, M.A., Kouroussis, D.A., Voutsinas, S.G., "Development , Testing and Evaluation of an Advanced Aeroelastic Numerical Tool for Horizontal Axis Wind Turbines", NTUA report, 1994


Figure 7: Radial distribution of the axial loading on the blade at several azimuthal positions (phi=0: blade oriented to the ground)


Figure 8: The effect of yaw angle in flapping moment variation


Figure 9: The influence of root vorticity emission


(a)


(b)

Figure 10: The influence of the yaw angle in the geometry of the wake


(a)


(b)


(c)


(d)


(e)


(f)


(g)


(h)

Figure 11: Effects of sheared inflow, tower existence and elastodynamic deformation on the loadings in yawed operation


(a)


(b)

Figure 12: Radial distribution of the induced velocities on the blade at several azimuthal positions with or without tower influence (phi=0: blade oriented tothe ground)


(a)


(b)


(c)


(d)

Figure 13: Loadings calculations for real cases of yawed operation