0179_CALCULATION OF AERODYNAMIC FORCE

CALCULATION OF AERODYNAMIC FORCE ON HORIZONTAL AXIS WIND TURBINE ROTOREXERTED BY TOWER EFFECTYutaka HASEGAWA *1, Junsuke MURATA *2, Hiroshi IMAMURA *3, Shinya HOTTA *4and Koji KIKUYAMA *5*1Ecotopia Science Institute, Nagoya UniversityFuro-cho, Chikusa-ku, Nagoya 464-8603, Japan, phone: +81-52-789-2711, fax: +81-52-789-5123e-mail: hasegawa@mech.nagoya-u.ac.jp *1*2,4Department of Mechanical Science, Nagoya UniversityFuro-cho, Chikusa-ku, Nagoya 464-8603, Japan, phone: +81-52-789-2712, fax: +81-52-789-5123e-mail: h054101d@mbox.nagoya-u.ac.jp *2 e-mail: h044123m@mbox.nagoya-u.ac.jp *4*3College Master Hands Inc. e-mail: ima@ *3*5Nagoya Sangyo UniversityABSTRACT: The existence of tower in horizontal axis wind turbine system has a large effect on surrounding fluid flow. In case of upwind wind turbines, the unsteady aerodynamic loads occur when the blades of a wind turbine are passing in front of the tower. The fluctuation of load causes a variation of power output and determines a fatigue life of the wind turbine, directly concerning the cost problems. The present paper describes a numerical method for calculation of aerodynamic loads on the rotor blades, an inviscid aerodynamic model based on the asymptotic acceleration potential method is adopted, which is capable of calculating the unsteady aerodynamic loads. In the model, the rotor blades are represented by spanwise and chordwise pressure distributions. The tower is represented by pressure distribution, in order to adapt the numerical model to that for calculation of aerodynamic loads on the blades.Keywords : Wind turbine, Tower effect, Numerical methods, Acceleration potential method, Load, Wind shear1. INTRODUCTIONThe existence of tower in horizontal axis wind turbine system has a large effect on surrounding fluid flow. In case of upwind wind turbines, the unsteady aerodynamic loads occur when the blades of a wind turbine are passing in front of the tower. The fluctuation of load causes a variation of power output and determines a fatigue life of the wind turbine,directly concerning the cost problems. Clarification of the tower effects on the aerodynamic loads, the resulting performance and fatigue damage of the horizontal axis windturbine rotor is imperative in the design process for the rotor blades. The present paper describes a numerical method for calculating the effect of tower existence on the aerodynamicloads, and shows the calculated results.For the calculation of aerodynamic loads on the rotorblades, an inviscid aerodynamic model based on theasymptotic acceleration potential method is adopted [1][2],which is capable of calculating the unsteady aerodynamicloads caused by the dynamic inflow [3], the yawed inflowsituations [4] and the combined inflow situations [5]. Underthe assumption of incompressible, inviscid and irrotational flow around the wind turbine, it can be derived that the pressure perturbation satisfies the Laplace equation and acts as an acceleration potential function. In the model the rotor blades are represented by spanwise and chordwise pressure distributions, which are composed of analytical first-orderasymptotic solutions for the Laplace equation. The tower is also represented by pressure distribution, in order to adapt the numerical model to that for the calculationof aerodynamic loads on the blades. 2. CALCULATION METHOD2.1 Load calculationBasic equations An inviscid aerodynamic model based on the asymptotic acceleration potential method is adopted for the calculation of aerodynamic loads on the blades. Assuming incompressible and inviscid flow around a wind turbine rotor, the governing equations are the Euler equation and the continuity equation. The Euler equation can be linearized by assuming that the induced velocities are small enough against the undisturbed velocity. Taking divergence ofthe linearized Euler equation and with the assumption of irrotational flow, the Laplace equation of the pressureperturbation is derived, using the Cartesian coordinate systemx ,y ,z () in Fig. 1, as follows;.02222222=++=∇zpy p x p p ∂∂∂∂∂∂ (1)The pressure perturbation acts on air particles as an acceleration potential function. The rotor blades can be represented by spanwise and chordwise pressure distributions, which are composed of analytical first-order asymptotic solutions for the Laplace equation.Integration of the acceleration experienced by the particles of air, which travel from far upstream to the rotorFigure 1: Coordinate Systems for Rotorblade, determines the velocities including the induced velocities in the rotor plane. With these velocities the aerodynamic loads can be calculated for the moment at which the air particles arrive at the rotor plane.Approximate solution methodThe Laplace equation of the pressure has been solved by adopting approximate solution techniques listed below:(i) asymptotic expansion method of the first-order for the pressure field around the rotor(ii) spanwise lift distributions expressed by finite series of associate Legendre functions.At considerable distance from the rotor blade, the experienced accelerations are equivalent to those caused by a lifting line, on which the blade load is concentrated. Therefore in the far field, the blade can be modeled with a pressure dipole line which gives the three-dimensional pressure distribution around the rotor blade, i.e.,.p far Close to the blade, the experienced accelerations are dominated by the chordwise load distributions. Thus the accelerations can be modeled with two-dimensional pressure distributions around the blade section in the near field.p near These two asymptotic approximations are combined with the use of the third term (a common field term ), which cancels the far field term near the blade, as well as the near field term far away from the blade. The sum of the three terms of the pressure perturbation gives an expression that is valid in the whole field.p common p composite =p near +p far −p common (2) In the application of asymptotes, only the first order term is considered for each field and the solutions are adopted in the present code.Using a local coordinate system shownin Fig. 1, the boundary conditions for the problem arex b ,y b ,z b () (3) p =0for x b2+y b 2+z b 2→∞∂p∂y b =0on the rotor blade (4) wΩr=θp on the rotor blade (5) where θp denotes a spanwise distribution of the blade pitch angle.The first condition, Eq. (3), is obvious from the definition of the function as the pressure perturbation. The second one, Eq. (4), means that the flow on a flat plate airfoil located in the p x b −z b plane does not experience the acceleration normal to the airfoil. The third one, Eq. (5), where w is the axial velocity including the induced velocity and r the radialposition, is the kinematic condition on the rotor blade, i.e.,Kutta’s condition.The instantaneous velocities of the air particles can be calculated by integrating the accelerations of the particles, which arrive at each collocation point, along their paths. The collocation points, where the kinematic boundary condition must be satisfied, are located on half chord lines of each blade at different spanwise position. However, the pressure field and the paths of the particles are unknown a priori, thus an iterative calculation for the pressure and the paths is performed at each time step.Viscous effect considerationIn order to take into account the viscous effects of the flow around the blade on the load distributions, a viscous routine is included in the preceding iterative process. The viscous routine makes use of two dimensional airfoil characteristics as a basis, and determines the local lift and drag coefficients according to the inflow angle against the local blade section.Tower effect considerationThe tower is represented by pressure distributions, in order to adapt the numerical model to that for the calculation of aerodynamic loads on the blades. One of the basic solution for the Laplace equation (Eq. (1)) is a pressure quadrupole, which is expressed using the coordinate system (x t , y t , z t ) in Fig. 2, as follows,⎟⎟⎠⎞⎜⎜⎝⎛−=52234t t t r r z W p πµρ (6) where µ is the strength of the pressure quadrupole , r t is the distance from the pressure quadrupole. It corresponds to doublet in velocity potential method, which can represent the flow around a cylinder by superposing on the uniform flow. In this work, the tower is represented as a distribution of pressure quadrupole located on the center axis of the tower.In the region close to the tower, the linear assumption which derives Laplase equation (1) is not valid, because the induced velocities are not small enough against the undisturbed velocity. Therefore, a corrective method is adopted for the calculation of particle path, as follows;.tower near t v W v t v tower ⎟⎠⎞⎜⎝⎛∂∂=⎟⎠⎞⎜⎝⎛∂∂ (7) By using a numerical method similar to the load calculation on the rotor blades, the instantaneous velocities of the air particles can be calculated by integrating the accelerations of the particles, which arrive at the stagnation points, along their paths. The boundary condition for the problem is taken as follows;point.stagnation the at 0w = (8)2.2 Inflow conditiontTwo types of inflow model have been adopted for the present paper. The first is uniform inflow perpendicular to rotor plane. The other is wind shear model which is represented by the power low [6], as follows;().shub H H y W y W α⎟⎠⎞⎜⎝⎛+= (9)where W hub is the wind speed at the hub height H and s αis the wind shear exponent. In the present calculation s α is fixed as 0.2.Figure 2: Coordinate Systems for Tower2.3 Calculation conditionThe number of collocation points on each blade is taken 6 and the time step for the time marching process is [s], which corresponds to the period of 1/180 rotation of the rotor.31065.4−×The calculation has been performed for the geometry of the NREL test turbine, which has two blades of radius R =5.029 [m] with the blade sections consisting of NREL S809 airfoil series, as summarized in Table 1. Its rotational speed in operation is [][rad/s 5.760/2rpm 72=]×=πΩ. The rotor configuration is shown in table 1.The tower consists of two different diameter cylinders and short conical section. The conical section base is 3.4 [m] and top is 3.9 [m] above the tower base. The radius of the upper cylinder is 0.2032 [m] and the radius of the lower cylinder is 0.3048 [m]. In this calculation, the radius of the tower is taken uniform as 0.2032 [m] for simplicity. The distance between the rotor plane and the central axis of the tower is 1.401 [m].The present calculation deals only with the aerodynamic loads, and the effects from the gravity force, the centrifugal force, and the force due to blade vibration are not considered.3. CALCULATION RESULTS AND DISCUSSIONS 3.1 Tower effects with uniform windIn this section, the inflow is taken to be uniform wind in order to clarify the effects only from the tower existence on the aerodynamic loads on the blades.Figure 3 shows the fluctuation of induced velocity experienced at three different radial positions, that is, near the root (r /R =0.35), the middle (r /R =0.60), and the tip (r /R =0.85) of the blade. The existence of the tower causes the reduction of wind speed in front of the tower. The increase of induced velocity can be found in the region close to the tower. Since the tip of blade has large relative velocity to the tower, it passes in front of the tower in short period. The fluctuation amplitude of the induced velocity is almost the same at the three radial positions.Figure 4 shows the difference of lift coefficient from no-tower situation at the three radial positions. The increase in of the induced velocity causes the decrease of the attack angle, as shown in Fig. 5. It brings about the reduction of lift coefficient about the rotor blade section. Even for almost the same fluctuation amplitude of the induced velocity, the reduction of lift coefficient at root is larger than that at tip. It is because that the relative velocity due to blade rotor rotation (Ω r in fig. 5) is smaller at blade root.Figure 6 and 7 show the fluctuation of flapwise moment (see Fig. 8) and torque in comparison with the experimentalFigure 4: Local Lift Coefficient at Three Radial PositionsTable 1: Specification of wind turbine rotor blade901802703600.080.100.120.140.160.180.200.220.24W hub =5.0 [m/s]I n d u c e d V e l o c i t y w i n d /W h u b [-]Azimuth Angle [deg]r/R =0.35 r/R =0.60 r/R =0.85Figure 3: Induced V elocity at Three Radial Positions901802703600.350.400.450.500.550.60W hub =5.0 [m/s]L o c a l L i f t C o e f f i c i e n t [-]Azimuth Angle [deg]r/R =0.35 r/R =0.60 r/R=0.85W 090180270360140145150155160165W hub =5.0 [m/s]T o r q u e [N m ]Azimuth Angle [deg] Calculation ExperimentFigure 7: Fluctuation of Torque (Comparing with Experimental Results)901802703608809009209409609801000102010401060W hub =5.0 [m/s]F l a p w i s e M o m e n t [N m ]Azimuth Angle [deg]Calculation ExperimentFigure 6: Fluctuation of Flapwise Moment (Comparing with Experimental Results) Figure 5: Relative V elocity and Attack Angle for Bladeresults [7], respectively. Extreme reduction of the flapwise moment and the torque can be found in the region close to the tower. These periodic fluctuation result in a fatigue load on the blade. Both the flapwise moment and the torque fluctuation components show good agreement with measurements in the region close to the tower.3.2 Effect of distance between rotor plane and towerThe distance between the rotor plane and the tower is the key parameter for considering the tower effects. Figure 9 shows the fluctuation of the flapwise moment at the blade root for various distances between the rotor plane and the central axis of the tower, ∆z (see Fig. 8) for no-tower situation. The tower effect increases with the decrease in the distance, and the slight increase in the flapwise moment can be found when the blade is approaching to the tower for ∆z /R t =4.0.In order to evaluate the effect on the blade load from the distance, the ∆M flp is defined as difference between minimum flapwise moment and that of no-tower situation, as shown in Fig. 9. Figure 10 plots the value of ∆M flp against non-dimensional distance ∆z /R t for various wind speed with approximated curves. It is shown that ∆M flp decreases with the distance approximately as (∆z /R t )-α. The exponent α is about 1.5-1.8. Generally, the average value of flapwise moment has almost linear relation with the inflow wind speed. Considering the flow field around a two-dimensional cylinder in uniform flow, the reduction ratio of wind speed to the undisturbed velocity can be written as (∆z /R t )-2. However, calculation results show different relation. It is because of the effects of interaction between the rotor and the tower as well as the effects of unsteadiness.3.3 Effect of tower radiusThe tower radius is also the key parameter for the tower effect along with the distance between the rotor plane and the tower. In order to evaluate the effect on the flapwise moment fluctuation from tower radius, the non-dimensional tower radius R t /R t0 is defined as the ratio of tower radius to the actual one (R t0=0.2032 [m]). Figure 11 shows fluctuation of flapwise moment for various tower radii. The fluctuation amplitude and the azimuthal extent of tower effect increase with the tower radius.Figure 12 plots the value of ∆M flp against non-dimensional tower radius R t /R t0 for various wind speed with approximated curves. It shows that ∆M flp increases with the tower radius approximately as (R t /R t0)β. The exponent β is about 2.1-2.2.90180270360950100010501100∆z =1.401[m]W hub =5.0 [m/s]uniform inflow F l a p w i s e M o m e n t [N m ]Azimuth Angle [deg]R t /R t0=0.5 Rt /R t0=1.0 R t /R t0=1.5Figure 11: Fluctuation of Flapwise Momentfor V arious Tower Radii50100150200F l u c t u a t i o n A m p l it u d e d u e t o T o w e r E f f e c t ∆M f l p [N m ]Nondimensional Distance ∆z /R t [-]Figure 10: Relation between Tower Effectand DistanceFigure 9: Fluctuation of Flapwise Momentfor V arious DistanceFigure 8: Flapwise Moment and Non-dimensional DistanceF l a p w i s e M o m e n t [N m ]Azimuth Angle [deg]050100150200250F l u c t u a t i o n A m p l i t u d e d u e t o T o w e r E f f e c t ∆M f l p [N m ]Non-dimensional radius R t t0/R [-]Figure 12: Relation between Tower Effectand Tower Radius50100150200F l u c t u a t i o n A m p l i t u d e d u e t o T o w e r E f f e c t ∆M f l p [N m ]Nondimensional Distance ∆z /R t [-]90180270360750800850900950100010501100115012001250W hub =5.0 [m/s]with wind shearF l a p w i s e M o m e n t [N m ]Azimuth Angle [deg]∆z/R t =4.0 ∆z/R t =6.0 ∆z/R t =8.0 without tower3.4 Tower effect with wind shearAn actual wind field experienced by a wind turbine installed in the open air has vertically sheared profile, namely wind shear. In this section, the inflow is taken to be sheared profile expressed by Eq. (9). Figure 13 shows the fluctuation of the flapwise moment for various distance ∆z together with the no-tower situation. When wind shear is considered, extreme reduction of the flapwise moment also can be found in the region close to the tower, resulting in an additional fatigue load.To compare with the uniform inflow case, ∆M flp is defined in the same manner as the uniform inflow case. It means the additional fluctuation amplitude of the flapwise moment due to the existence of the tower. Figure 14 plots the value of ∆M flp against non-dimensional distance ∆z /R t for various wind speed with approximated curves. The values of ∆M flp are slightly smaller than these of the uniform inflow case, because the lower half of rotor plane, where the tower effects on the surrounding flow, is lower wind speed region in the case of the wind shear. It shows that ∆M flp decreases with the distance approximately as (∆z /R t )-α. The exponent α is about 1.5-1.7. The relation is almost the same as the uniform inflow case.3.5 Tower effect for tilted rotorMany wind turbines have tilt angle for their rotor axes. Tilt angle stands for the angle between the inflow direction and the axis of rotation of rotor in vertical plane, as shown in Fig. 15. Figure 16 shows the fluctuation of flapwise moment for various tilt angles. The inflow is taken to be sheared. Thetilting of the rotor reduces the fluctuation amplitude of the flapwise moment due to the tower effect. It is because that the distance between the rotor blade and the tower increases with the tilt angle in the lower half of the rotor plane. 4. CONCLUSIONS The effects of the tower on the aerodynamic behavior for the upwind turbine rotor have been examined by performing the numerical calculation. The major results obtained are summarized as follows.(1) The existence of tower causes the reduction of wind speed in front of the tower. It brings about decreases of attack angle and lift coefficient of the rotor blades.(2) The Flapwise moment on the rotor blade is extremely reduced in the region close to the tower.(3) The fluctuation amplitude of the flapwise moment decreases with the distance approximately as (∆z /R t )-α. The exponent α is about 1.4-1.8.(4) The fluctuation amplitude of the flapwise moment increases with the tower radius approximately as (R t /R t0)β. The exponent β is about 2.1-2.2.(5) The fluctuation amplitude of the flapwise moment due to the tower effect is decreases when the wind shear is considered.(6) The tilting of the rotor reduces the fluctuation amplitude of the flapwise moment due to the tower effect.Figure 14: Relation between Tower Effectand Distance for Sheared InflowFigure 13: Fluctuation of Flapwise Momentfor Sheared InflowFigure 15 Definition of Tilt AngleFigure 16 Fluctuation of Flapwise Moment90180270360750800850900950100010501100115012001250with wind shearW hub =5.0 [m/s]F l a p w i s e M o m e n t [N m ]Azimuth Angle [deg]tilt angle=0 [deg] tilt angle=3 [deg] tilt angle=6 [deg]for V arious Tilt AngleACKNOWLEDGEMENTThe National Wind Technology Center (NWTC) is gratefully acknowledged for provision of experimental data. REFERENCES1. van Holten, Th., The computation of aerodynamic loads on helicopter blades in forward flight using the method of the acceleration potential, Report VTH-189, Delft Univ. Tech., Delft, The Netherlands, 19752. van Bussel, G. J. W., The use of asymptotic acceleration potential methods for horizontal axis wind turbine rotor aerodynamics, Proc. EWEC’91, Amsterdam, pp. 18-23, 19913. Hasegawa, Y. and van Bussel, G. J. W., Application of the Asymptotic Acceleration Potential Method for the Calculation of Dynamic Inflow Effects on Horizontal Axis Wind Turbines, 4th Int. Symp. Transport Phenomena Dyn. Rotat. Mach., Honolulu, 186-193 (1992).4. Hasegawa, Y., et al., Numerical Analysis of Yawed Inflow Effects on a HAWT Rotor, 3rd ASME/JSME Joint Fluids Engineering Conference, San Francisco, FEDSM99-S295-9 (1999).5. Hasegawa, Y., et al., Aerodynamic Loads Calculation of a Horizontal Axis Wind Turbine Rotor in Combined Inflow Condition, EWEC04, London, Proceeding CD-ROM (2004).6. IEC 61400-1, 2nd ed., Wind Turbine Generator Systems - Part 1: Safety Requirements, International Electrotechnical Commission (1999).7. NREL’s National Wind Technology Center Unsteady Aerodynamics Experiment Database.。