Dynamics of Vertical Axis Wind Turbines ( Darrieus Type )

A computing package that combines finite element methods for evaluating the resonance frequencies and modes of turbine subcomponents (blade, tower and shaft) together with the aerodynamic calculations for forces and moments taking into 
consideration the dynamic stall as well as the dynamic response is developed. This method was applied to a realistic VAWT; namely; the PIONEER I built in the Netherlands by Fokker company. A reasonable agreement between the calculated and field results was predicted.


INTRODUCTION
o meet economic levels of power generation, the wind turbine designer is being forced towards compliant designs.As a consequence of this fact, more emphasis is being placed on both aeroelastic and dy- namic analyses in the design phase of producing a new product.Since a wind turbine blade is highly dynami- cally loaded, it requires intensive investigations before being put into production to ensure both safety and reliability.This paper is devoted towards the investiga- tion of the aerodynamics and stability of a Darrieus vertical axis wind turbine (VAWT).VAWT's have sev- eral advantages over the horizontal axis wind turbines; (HAWT's) like the indifference to wind direction, the proximity of the gearbox and generator to the ground and the non-reversal gravitational stresses in their troposkien shaped blades.However, the main disadvantage of VAWT's is their disabilities for self starting.
At the Vrije Universiteit Brussel (V.U.B) a substantial effort has been done to investigate both the aerodynamics and dynamics/aeroelasticity of both HAWT's and VAWT' s.Concerning the specific work with VAWT' s, a computing package is developed which combines the versatile finite element methods (FEM) for evaluating the natural frequencies and modes of the turbine's subcom- ponents; Bathe [1982], modal coupling technique; Ben- field, et al. [1971], to define the behavior of the whole structure, the aerodynamic calculations for loads on blades, eigenvalue problem solver for natural frequencies and modes evaluation as well as a stability investigation via a forced vibration analysis.This package is analogous to that developed by Garrad, et al. [1984].In this regard it can substitute using the well-known code NASTRAN (MSA/NASTRAN User's manuals 1981 ]), and also the European code ARLIS, Kirchgi[ner [1984].Extensive analytical and experimental research works regarding the aerodynamics, aeroelasticity and dynamics of Darrieus VAWT's are performed both in Sandia laboratories; Lobitz [1986] and Popelka [1982], in USA and in the Netherlands Energy Research Foundation (ECN) in Neth- erlands; Machielse [1984]  and Machielse, et al. [1986].
Most investigations elsewhere compare their results with the findings published by Sandia and/or ECN.In this paper, our model is applied to the 15 m VAWT manufac- tured by ECN and named Pioneer I, Ottens [1981].

MATHEMATICAL MODEL
A mathematical model capable for investigating both of the aeroelastic and dynamic responses of a Darrieus VAWT is presented here.The equations of motion for a small displacement (q) from a steady state for the discretized model may be derived either from the energy equation of continuum mechanics, or the principle of virtual displacements leading to a Lagrangian form of governing equation.Both techniques lead to the same equation.The following assumptions are made: 1.The blade has a troposkien shape; refer to Blackwell,  et al. [1974]. 2. Blade cross section is symmetric about its major principal axis and can have three distinct points; namely, the elastic, aerodynamic and tension centers.The elastic and mass centers are coincident. 3. Blade material is isotropic and linearly elastic. 4. Blade cross sections which are normal to the elastic axis before deformation remain plain after deforma- tion (Bernoulli-Euler hypothesis).Also torsional warping is negligible. 5. Blade can bend in two mutually perpendicular direc- tions normal to the elastic axis. 6. Deflections, strains and rotations are small. 7. Quasi-steady blade-element strip theory is applicable.

Frames of Reference
The first frame of reference is an inertial one So(Xo, Yo, Zo) fixed to the undeflected rotor.Next a series of transformations including rotations and translations are considered.S I(X1, Y1, Z1) has the same origin but rotated an angle 0 w.r.t, the inertial frame So.Next $2(X2, Y2, Z2) is shifted a distance g from the frame S 1.While $3(X3, Y3, Z3) and $4(X4, Y4, Z4) are inclined angles + and w.r.t, the Y2 and X3 respectively.Finally the frame $5(X5, Y5, Z5) locates any point on the blade itself whose tangent is inclined angle -,/ to the Z4 direction.The different frames of reference are related to each other by the equation.

Finite Element Method
Treatment of such flexible continuous body can be performed by considering its deformation as the sum of few number of its normal modes multiplied by general- ized coordinates.This method is known as the MODAL METHOD.This method is advantageous over the lumped mass deformation as in the former only few number of the generalized coordinates or in other words unknowns are assigned.Usually, vibration modes of non-rotating turbines are adopted.To calculate these vibration modes, finite element methods are applied.Vibration modes of the components of the whole stBc- ture; here the blades, tower and the shaft together with the transmission and generator, are calculated.This is described as the "component mode method".
The subcomponents of the wind turbine are discretized into finite number of beam elements.Each beam element is characterized by two nodes and 12 degrees of freedom (6 rotational and 6 translational).The coordinates of each node together with the cross-sectional area and moments of ineaia about three mutually pewendicular axes and the material density distribution are fed to a stBctural code.Consequently, the mass, flexural and torsional rigidities of each element are calculated.A powerful code developed by the stBctural depament of the Vfije Universiteit BBssel; Belgium, known as DAPST (which is analogous to NASTRAN) is used to calculate the natural frequencies and modal shapes of the three sub- components (refer to its user's guide [1985]).
Since blades are axisymmetfic, then the same mode shape can be used for all blades.The displacement anor rotation of a point on the kth blade and tower are denoted by Sik and ; where varies from to n and varies from to m.The generalized coordinates of the shaft are its angular displacement (0), and the torsional deflection of the mechanical transmission and rotating part of the generator (A02).

Equations of Motion in Lagrangian Form
Subsequent to these steps the kinetic and potential energies are derived as described in [12].Since the governing equations in the Lagrangian form have the form: where T and V are respectively the kinetic and potential energies of the whole turbine, while hi are the generalized coordinates (A0,A02,qj,Sik).Forces are divided into conservative and non-conservative forces.Here both of the aerodynamic and gravitational forces represent the conservative forces.These conservative forces are derived from the potential energy.The non-conservative forces (except for the viscous damping) are expressed by Qi.The viscous damping are expressed by 0F / 0hi where F is the well known Rayleigh's Dissipation Function (Strickland 1975]).

Aerodynamic Force
The aerodynamic forces are defined both theoretically (refer to Strickland [1975], Sharpe [1977], Templin [1974], Simhan [1984] and Mandal [1986]) and experi- mentally (refer to Worstell [1979], Muraca, et al. [1976] and South, et al. 1972]).The double-multiple stream tube theory [23, 24] is employed in this model.In such a model the Darrieus VAWT is represented by a pair of actuator disks in tandem at each level of the rotor as illustrated in Figure (2).The elements of normal and tangential forces (SNk and 8Tk) and the resulting mo- ment on a strip (SMk) through the principle of virtual work provide the generalized forces QA0, Qqi and Qski.Dynamic stall will be discussed later on.

Gravitational Force
The forces due to gravity are the second conservative force.The generalized forces due to gravity associated with the tower and blades are, GA0, Gqi and GSik respectively Modal Coupling To apply the modal coupling technique; Benfield, et al.  [1971], the generalized coordinates are assembled into

Equations of Motion
After a lengthy mathematical deductions, equation (2)   will have the form: where [k] are the structural mass, damping and stiffness matrices respectively.[g], [s] are the coriolis (gyroscopic) and geometric stiffness matrices.
[Gr], [c]  are the gravitational and centrifugal matrices.
Moreover, Q represents all the forces, which is now a function of both space and time (due to the aerodynamic forces on blades).Calculations are now performed in the time domain using the classical 4th order Runge-Kutta integration scheme, with automatic time step integration.

Aerodynamics with Dynamic Stall Calculations
As was pointed out earlier, the aerodynamic for- ces are calculated with the double multiple streamtube theory.At higher wind speed this theory fails to predict the forces on the turbine correctly due to dynamic stall effects.Therefore, the present theory has been extended, to account for these effects.The circular path of the blades of a vertical axis wind turbine cause periodic oscillations of the angle of attack at a frequency of once per revolution and at higher multiples of this frequency.This rapid variation introduces unsteady flow phenomenon, that are not accounted for by the classical theories.They are known as dynamic stall.Stall is a phenomenon that occurs at higher angles of attack (of the order of 10 degrees) when the flow separates from the body.As a result a pressure drop and thus a decrease in lift will occur.When this occurs at a relatively slow angle of attack change, the phenomenon is called "static stall".When the angle of attack changes more rapidly unsteady flow phenomenon introduces hysteresis effects as shown in figure 3. When the angle of attack increases rapidly above the static stall angle, the separation of the flow is retarded and stall will occur at higher angles of inci- dence.The lift forces will be higher than in the static case.In order to take into account dynamic stall effects, the Boeing-Vertol dynamic-stall model, modified by Strickland [1975], is used.This model assumes that the lift-curve slope and the zero-lift angle remain unchanged and that the only change due to dynamic effects is the angle of attack at which stall occurs.Thus a modified angle of attack is utilized for use in entering two- dimensional force coefficient data.
In this method dynamic stall corrections are consid- ered only for angles of attack above 5 degrees and below three times the static stall angle.Lift and drag coeffi- cients are calculated first for the static case and then adjusted for dynamic stall effects.A modified angle of attack om is calculated as: where: Ae K. [e.C/(2 U)] 5 + K K, K1 are empirical constants, depending on the Mach number and the thickness ratio of the profile.For the unstalled fully stalled "... NACA 0012 profile the actual values are K 74.5, K1 4.47.G is a weighing factor depending on which part of the dynamic stall hysteresis loop, the cycling value of am lies.(See Figure 3).[1 + (. Once the modified angle of attack is calculated, the lift force coefficient CLs is calculated for stalled conditions.Finally the modified lift force coefficient is calculated from the static value CLo and the stall value CLs as" Ce C0. (1 P) + Cs.P where: P is defined in figure 4. (6)

RESULTS
In this section the dynamic response of blades and tower of the Pioneer I to the external forces is presented.These external forces are the centrifugal, gravitational and aerodynamic forces.

CASE STUDY
The mathematical model described above is now applied to one of the most documented wind turbines, namely PIONEER I, built in the Netherlands by Fokker com- pany.Figure 5  discretized model.PIONEER I is a two-bladed VAWT with rotor height of 15 m and diameter of 14.92 m, has a troposkien blade shape with a blade chord of 0.75 m.
The main bearings on the top of the tower and in the rotor center are situated at the inside of the rotor tube.Both of the rotor tube and the turbine tower are made of two steel tubes of different length, diameter and wall thickness.
The blades have a NACA 0012 cross section and made of glass fiber reinforced polyester (GFRP) skins and a polyurethane foam core.Operational speed varies from 25 to 50 rpm at a wind speed of l0 m/s.Power output at 18 m/s is 65 kW.

Forces and Power
Figure 6 shows the power coefficient versus wind speed for 30 and 50 rpm.(without dynamic stall).At 30 rpm a maximal power is extracted at a wind speed of 7 m/s.This power is approximately 10 kW.At 50 rpm the power increases to 53 kW at a wind speed of 12 m/s.In figure 7 the power at a rotational speed of 40 rpm is shown, as calculated without and with dynamic stall.It may be seen that for higher wind speeds the aerodynamic power is considerably higher when dynamic stall is taken into account.
The normal and tangential forces on the blades as a function of azimuthal position at a rotational speed of 50 m/s and a wind speed of 10 m/s are presented in figures 8 and 9 respectively.At the upstream side these forces are considerably larger than at the downstream side.

Dynamic Response
In figures 10 to 12 the generalized forces and dynamic blade responses due to the centrifugal forces are shown for a two mode representation of the blades.As can be seen from figure 10 the centrifugal force is constant and thus leads to a constant generalized force.In principle, the centrifugal force is balanced by the gravitational forces, due to the special geometry of a troposkien shape.This leads to a general force, that should be equal to zero.However due to numerical rounding errors, a small force remains that will result in a very low dynamic response, as can be seen in figures 11 and 12.This response is characterized by a damped oscillation level at transient, and a constant level after some settling time.This type of response is typical for a damped oscillation to a constant force.
In figures 13 to 15 the dynamic response of the blade to the gravitational forces are shown.As the centrifugal forces and gravity introduce a constant force into the blades of the turbine and into the tower, it will introduce a damped oscillation with constant level in the flapping direction.In lagging direction a zero displacement level should be found.However, due to the modal approxima- tion and rounding errors, a small amplitude of the order of 10-6 remains.
In figures 16 to 18 the displacements of the midpoint of the blade and the generalized forces, due to aerodynamic loading are presented.As can be seen the aerody- namic loading is more or less periodic, with a maximal amplitude when the blade is in the upstream region.For a damping ratio of 0.1, the flapping displacement of the blade is symmetric, with a nearly zero amplitude at mid span.The lagging displacement follows the patt,ern of the aerodynamic forces with a period of 2 per revolution.
In figures 19 to 21 the generalized forces and displace- ments of the blade to all external loading are presented.As can be seen, the displacements of the midpoint of the blade in flapping motion are dominated by the gravita- tional forces, while the displacements in lagging motion are mainly due to the aerodynamic loading.CONCLUSIONS 1.The designed package which is described here predicted its validity and accuracy in evaluating the aerodynamic, aeroelastic and dynamic behaviors of Darrieus VAWT' s.
2. Dynamic stall has to be incorporated in the aerodynamic calculations since it drastically changes the extracted power compared to models ignoring it. 3. Centrifugal forces on blades are counteracted by gravitational ones.The transient behavior due to them shows a damped response.

REVDLUT
FIGU 17 Flapping displacement of the centre of the blade due to aerodynamic loading.FIGURE 18 Lagging displacement of the centre of the blades due to aerodynamic loading.
4. Aerodynamic loading is periodic with maximum amplitude when blades are in the upstream region. 5. Displacements in the lagging and flapping directions are dominated by aerodynamic and gravitational forces respectively.

ACKNOWLEDGEMENT
The first author acknowledges the Financial support given to him by Vrije Universiteit Brussel while he was visiting professor to the fluid mechanics department in 1986.
The support of Mr. Riad M. E1-Sharawi; Enppi Co., during the preparation of this paper is greatly acknowl- edged.AZIMUTH FIGURE 19 Generalized forces on the blades due to all loading.

Figure
FIGUDifferent frames of reference.

FIGURE 2
FIGURE 2 Darrieus rotor geometry and the double-multiple stream tube.

FIGURE 3
FIGURE 3 Lift force and weighing factor for hysteresis loop.
FIGURE4 Weighing factor E