titre de la communication

1ère Conférence Franco-Syrienne sur les énergies renouvelables

Damas, 24-28 Octobre 2010

WIND TURBINE FARM OPTIMIZATION USING ACTUATOR DISC AND ACTUATOR LINE MODELS

Munif Jourieh

Patrick Kuszla

Ecole Nationale Supérieure d’Arts et Métiers 151 Bd de l’Hôpital, 75013 Paris FRANCE [email protected]


Ivan Dobrev

Fawaz Massouh



Abstract The interaction between wind turbines in Wind Park is a crucial problem to be taken into account in preparation for the park construction. Due to important advances in numerical simulation, the study of this problem is achievable. However, it is very difficult to simulate many wind turbines with real geometry in a reasonable time. Now, for this reason it is necessary to use simplified hybrid models able to represent one or more wind turbines. Two hybrid models are developed to predict the aerodynamic performance of a horizontal axis wind turbine. A tridimensional simulation is thus coupled to a blade element method. This coupling uses two separate models: a model of actuator disc and a model of actuator line. These models are favorably compared with NREL S809 experiments and other simulations. These models allow us to access to the performance of several wind turbines with interacting wakes. Key words : Hybrid models, wind turbine, wake interaction, wind park.

1. Introduction Wind energy is one of the cleanest energy in the world. Since several decades, this kind of energy began to take its real place among the others. The need for this clean energy was lead to fast development of wind energy technology. To minimize the cost of energy production high power units are required and nowadays the largest wind turbines have 5 MW rated power and rotor diameter of 126 m. Also, to optimize power production, several machines must be installed on the same wind farm. Generally the wind farm development is a very complex problem and variety of factors influenced the wind turbines positioning. However in all cases it is needed to avoid the negative effect of aerodynamic interference between them. Therefore, it is preferable to investigate the wake development downstream of each wind turbine in the wind farm. Numerical simulation of the flow around even one wind turbine is computationally very expensive. Accurate results can only be obtained using very fine mesh around the blades and in their wake. As a consequence, modeling a wind turbine farm with more than 2 or 3 machines is technically impossible. Hence, to overcome this limitation, the rotors must be represented by numerical models involving less grid points. These kinds of models, without compromising the important features of the flow structure behind the rotor, are so called hybrid models. The hybrid models are frequently used in interactional aerodynamic of helicopters to obtain the interference fuselage-rotor. In the case of wind turbine aerodynamics the hybrid models are also applied [1]-[5]. The hybrid model comprises two modules. In the first module a CFD solver computes the velocity field of the domain of investigation. Here the presence of the rotor is modeled with body or surface forces. To calculate the rotor forces, the second module uses a conventional method based usually on the blade element method (BEM) or vortex methods. These methods use only the local blade element lift and drag coefficients and it is not needed to model flow around the blades. As a consequence, the needs for computational power decrease significantly. In this paper two simple hybrid models are studied. The first hybrid model replaces all the complicated geometry of the wind turbine rotor by a thin cylindrical volume; this is the so called actuator disc model. The CFSER 2010

Page 1



second hybrid model replaces each blade by an equivalent cylinder placed along the blade axis; this is called the actuator line model. These models are implemented in the commercial CFD code “Fluent” 6.3, where the wind turbine rotor is replaced with relevant body forces. These forces are obtained by means of the blade section airfoil performances and the rotor inflow. Then they are implicated in the computational domain by means of a user-defined function (UDF). 2. Hybrid models The first model considers replacing the rotor geometry by a thin cylinder with the same diameter; this model is a kind of the actuator disc model (AD). The second one is so called the actuator line model (AL). Here, circular cylinders along the blade axis replace the blades. The first model is simple to apply. The second is more complicated but is able to represent the tri-dimensional wake structure. However in case of wind farm modeling, unsteady solver with “moving mesh” capability must be used. In the two cases, we compute the forces exerted by the blades and then we distribute these forces on the corresponding geometry. 2.1. Actuator disc model Here the radius of the actuator disc is the same as the wind turbine radius R and its axial thickness is equal to e, this axial thickness is calculated depending on blade chord and blade pitch angle[13]. To replace the blades in the CFD solver, the actuator disc must exert the same aerodynamic forces as the real rotor. To obtain these results, we use the BEM [4]. According to the velocity triangle shown on Fig. 1, the elementary lift and drag forces created by one blade element with chord c and located at the radius r, are

1  ρ cCL ( ) Wref2 dr   2  1 dD  ρcCD ( ) Wref2 dr   2  dL 

where CL CD 

(1)

blade section lift coefficient, blade section drag coefficient, air density.

Fig. 1 Blade element velocities and forces The relative velocity Wref in (1) can be obtained using upstream velocity V0, rotor angular velocity Ω, axial wia and rotational wit induced velocities

Wref  (V0  wia )2  ( r  wit )2 .

(2)

The angle of attack is determined as

CFSER 2010

Page 2



    .

Here β is the blade pitch angle and  is the flow angle in the plane of rotation. The flow angle can be obtained from   arctan

V0  wia .  r  wit

(3)

Our method is different here from the classic BEM, i.e. the terms (V0-wia) and wit is obtained from the flow field computed by CFD code. To get the sectional lift and drag coefficients C L and CD, experimental data concerning aerodynamic properties of the blade section airfoil are usually used. The coefficients of aerodynamic forces; thrust and tangential are expressed as follows

CQ  CL  cos   CD  sin   . CT  CL  sin   CD  cos   

(4)

To evaluate the intensity of the body forces within the actuator disc, we divide it in n annular elements of radial thickness dr. Then the elementary blade element forces are averaged azimuthally for each annulus. For an elementary annulus of radius r the volume is dv  2 re dr .

(5)

Taking into account the blade number N and using (1), (4) and (5), the annulus thrust and torque body force intensities are expressed as

cCQWref2 N  dQ   dv 4re  . cCT Wref2 N  dT fT     dv 4re  fQ  

(6)

Due to CFD solver constraints, the body forces must be distributed in the actuator disc by referring the volume cells in a Cartesian coordinate system. Projecting these forces on X, Y and Z-axis, the force intensity at arbitrary point of the volume becomes

 sin  F  4re   N cCT Wref2  fy   F  4re   N cCQWref2 fz   cos  F  4re  fx 

N cCQWref2

(7)

Here the angle  is azimuthal coordinate of actuator disc control point, where the body forces are to be computed. In (7) F is the Prandtl tip loss factor and is applied to account the finite number of blades. The approximate formula for this factor may be found in [6][7] such as

F

 N Rr  2 arccosexp      2 r sin  

(8)

It must be noted that there are some inconsistencies in the obtained results involved by this correction. An improved tip loss correction method allowing a better prediction of the physical behavior of the flow in vicinity of the blade tip was presented recently by Sorensen and all [8] but has not been implement in the present work. The computing procedure is organized as follows. Starting from inflow obtained by CFD, the BEM module calculates the resultant velocity, flow angle and angle of attack. The airfoil lift and drag coefficients, for a given angle of attack, are interpolated using experimental or computational data. To obtain body force intensities (7), the torque and axial coefficients (4) are calculated and the Prandtl correction (8) is applied. Then the CFD computes again the flow field, applying intensity distribution obtained from the BEM on the actuator disc volume. Hence the solution is carried out using an iterative method, exchanging data between the BEM and CFD modules; it stops after convergence is reached.

CFSER 2010

Page 3



Fig. 2 Actuator disc model The important difference between the classical blade element momentum theory (BEMT) and the hybrid methods is in the calculation of the induced velocities. In BEMT they are calculated independently for each elementary annulus using momentum theory [6], [7]. As a consequence, the interference between aerodynamic forces in spanwise direction does not exist. Contrarily, in hybrid models the CFD takes into account the distribution of force intensity in spanwise direction in same time. Another advantage of the proposed method is the finite thickness e of actuator disc. In this volume two important flow characteristics as turbulence kinetic energy and turbulence dissipation can be involved. Their intensity may be evaluated as a function of the flow condition downstream the blade element; separated or non-separated flow. Therefore it is possible to influence significantly the rotor wake development. In this study, the influence of the actuator disc thickness e is not investigated. In all case of simulation it is equal to 0.10.2 of the blade chord [13]. As the body forces are averaged azimuthally, strong flow gradients do not exist. 2.2. Actuator line model This model replaces the geometry of each blade by a cylindrical volume positioned along the blade axis. The cylinder has the same length as the blade and the same cross-section A as cross section of the blade. In this cylinder the body forces, which replace the blade are distributed in the spanwise direction. Here at the same radius, the body forces contained in the elementary volume of the cylinder are equal to the force exerted by the corresponding elementary blade element. From an aerodynamic modeling point of view, it is necessary that flow perturbations created by the body forces are similar of those created by the blade. As a consequence we expect that the wake created by the cylinder will be comparable the wake created by the blade. To find the optimal cylinder cross-section it is easier to study the two-dimensional flow past an airfoil. Here it is needed to find the radius of a circle with body forces distributed inside. This circle should generate a velocity field as close as possible to that created really around the studied aerofoil. During all simulations the resulting field of body forces must be equivalent to the airfoil aerodynamic force. In this study we carry out numerical simulations in the case of a flow around a typical wind turbine airfoil S809. It was found that for three different angles of attack 6°, 12° and 18° the minimum differences were obtained in the case of a circle of radius equal to 0.05c. It is clear that it is not possible to model the boundary layer wakes or flow separation. However it is found that a circle of large radius creates a velocity field without sufficient velocity gradient. On the opposite, in the case of a too small radius, the body forces create strong velocity gradients that increase viscous damping. Also to model a circle with low radius and strong body force intensity, more grid points are needed inside. In all cases it seems that the approach applied here is more acceptable from the point of view of existing phenomena than that approach used in [2]-[4] where 2-D Gaussian distribution of body forces is applied in a plane normal to the blade axis. The coefficients of aerodynamic forces for an elementary blade section at each radius are still given by (1) and (4). The aerodynamic elementary forces are now distributed as body forces at the radius r inside an elementary cylindrical volume dv with thickness equal to dr

CFSER 2010

Page 4


dv  Adr .


(9)

Then we obtain the volume force intensity 2 dQ cCQWref    dv 2A  . 2  dT cCT Wref  fT     dv 2A 

fQ  

(10)

Fig. 3. Actuator line model Projecting these forces on X, Y and Z-axis, the force components for blades becomes

 sin   2A   cCT Wref2 fy   . 2A   cCQWref2 fz   cos   2A  fx 

cCQWref2

(11)

In the case of the calculation of a single wind turbine without ground influence it is more practical to use the rotor rotational symmetry and to model only one blade period. Then if the blade axis is coincident with one of coordinate system axis it is possible to use (10). In all other cases the rotational symmetry does not exist, the simulations are unsteady and (11) must be applied to obtain the body forces intensity inside the cylinders. 3. Reference values and cases To validate the results obtained by the hybrid models developed here, numerical simulations are carried out in the case of the NREL Phase II wind turbine [9]. The tested machine is a three-bladed horizontal axe wind turbine. The rotor has a diameter of 10.1 meters and untwisted blades with a constant chord of 0.458 meters. The rotational speed of the rotor is 72 rpm. The blade pitch angle is fixed at 12°. The wind turbine tower has a 0.458 m diameter and a height of 16.8 meters above the ground. To compute the blade aerodynamic forces we use the S809 airfoil data obtained in Delft University of Technology low-turbulence wind tunnel [10]. Usually the BEM are based on 2-D airfoil data and give realistic results in case of angle of attack below the stall. When this condition is not respected the flow past the blades becomes tri-dimensional. Here the detached flow, which is rotated with blade, creates radial flow in spanwise direction. As a result, the Coriolis force affects the flow in vicinity of the blade suction surface. The influence of these forces leads to significantly higher lift coefficients than those obtained during wind tunnel tests, in particular at the blade root. Actually there are several studies concerning the methods of correcting airfoil 2-D data. Also, for S809 airfoil there are experiments, which allow extracting airfoil data on the base of measured blade pressure distribution.

CFSER 2010

Page 5



For example, the comparison of the lift and the drag coefficients obtained for different r/R and the wind tunnel data can be found in [11]. Because our study concerns only the applicability of hybrid modeling so we will use the classical pre and post stall values of Viterna, given in [11]. Our models are implemented in Fluent 6.3 using User-Defined-Functions. All numerical results are obtained in case of incompressible fluid. We used a k- SST turbulence model. 4. Numerical results and comparisons Experimental data from NREL Combined Experiment Phase II is plotted against the obtained numerical results, Fig. 4. The actuators disc model is computed with and without Prandtl correction. The power predicted by this model with Prandtl correction is in good agreement with the NREL for all wind speed conditions, from well-attached flow to the heavy stalled cases. We also compare the numerical results of the actuator line model and the NREL experimental data. As can be seen on the Fig.4, there is a good agreement between them, for low and high wind speeds. It is interesting to compare the results of actuator disc and actuator line models. For low wind speed the same results are obtained; but the difference becomes more important after the velocity of 12 m/s. Here the results obtained by the actuator disc model are more precise. The possible reason of this is the difficulty to determine correctly the reference plane where the vector of velocity is obtained in the case of the actuator line model. Finally we compare our two hybrid models with the results obtained in [12] for the same wind turbine using the Navier-Stokes code of Shu and Shankar. From the presented results on Fig. 4, it can be seen that this model, especially developed for wind turbine performance analysis, gives good results, comparable to those presented here.

power (kW)

20

15

10

Actuator disk Actuator line NREL Experiment N-S Solver, Xu and Sankar

5

0

5

10

15

20

25

wind speed (m/s)

Fig. 4. Comparison of power as predicted by computation methods and experimental results.

5. interaction between 6 machines in cascade Here we present the axial velocity field for six wind turbines replaced by the actuator disc hybrid model DA. The distance between two consecutive machines is equal to 5D; upstream wind velocity is equal to 10 m/s. The turbines type is NREL phase II. On Fig 5, after the first wind turbine, we observe that each machine is located in the wake of another one, and the axial velocity deficit behind each machine gradually increases with the distance from the field entrance. The deficit between the first and second one is greater than between the second and third and so on. The power results from these machines are presented in Fig 6; we observe that the loss between two consecutive machines decreases downstream. In this figure, we plot the power of the i th machine divided by the power of the first one. The power difference between the first and second machine is 10%, the difference between the third and fourth machines is 5%, the difference between the fourth and fifth is 1.5% and the difference between the fifth and sixth machine is less than 1%.

CFSER 2010

Page 6



+ Fig. 5. Axial velocity field for six machines in cascade.

Fig. 6. Power evolution for six wind turbines located at the same rotation axis,

CFSER 2010

Page 7



6. Conclusions and perspectives

We have implemented two rotor models in an industrial CFD code. These models are based on the blade element theory. They produce promising results, which are validated by a comparison with the experimental data obtained in the case of NREL S809 phase II wind turbine. Despite the good agreement between our calculations and the experiments, there are some differences. These differences are probably due to the model simplifications introduced to avoid a complete and expensive representation of the real geometry. It must be noted that the aerodynamic forces distributed on the control volumes defining the equivalent rotors are evaluated using 2D characteristics. Better results should be obtained using the characteristics of rotating airfoils. Finally it must be noted that calculation speed and memory usage allows us to compute the aerodynamic performance of a wind farm with several machines. Further work will be now aimed at the validation of the downstream wake of the rotor. Some computations should be also made in more complicated situations taking into account the ground geometry and some more complex interactions between the wind turbines. 7. References

[1] L.J. Vermeer, J.N. Sørensen and A. Crespo, Wind turbine wake aerodynamics, Progress in Aerospace Sciences, Vol. 37, Aug.-Oct. 2003, pp 467-510 [2] J.N. Sørensen and W.Z. Shen, Numerical Modeling of Wind Turbine Wakes, Journal of Fluids Engineering, Vol. 124, Jun. 2002, pp.393-399 [3] R. Mikkelsen, Actuator Disc Methods Applied to Wind Turbines, Ph.D. Thesis, Technical University of Denmark, [4] S. Ivanell, Numerical Computations of Wind Turbine Wakes, Technical Reports from KTH Mechanics, Royal Institute of Technology, Stockholm, Sweden, April 2005 [5] F. Massouh, I. Dobrev and M. Rapin, Numerical Simulation of Wind Turbine Performance Using a Hybrid Model, 44th AIAA Aerospace Science Meeting and Exhibit, Reno, Nevada, Jan. 2006 [6] Glauert H. Airplane Propellers, Vol. IV, Div. L, in Aerodynamic Theory, edited by Durand W.F., Dover Pub. Inc, NY, 1963, p. 434 [7] T. Burton, D. Sharpe, N. Jenkins, E. Bossanyi, Wind Energy Handbook, 2001, Wiley & Sons, p. 617 [8] W.Z. Shen, R. Mikkelsen, J.N. Sorensen and C.Bak, Tip loss Corrections for Wind Turbine Computations,, Wind Energy, Vol. 8, No. 4, Oct.-Dec. 2005, pp. 457-475 [9] E.P.N. Duque, W. Johnson, C.P. van Dam, R. Cortes and K. Yee, Numerical Predictions of Wind Turbine Power and Aerodynamic Load for the NREL Phase II Combined Experiment Rotor, AIAA-2000-0038 [10] D.M., Somers, Design and Experimental Results for the S809 Airfoil. NREL/SR-440-6918. Golden, Colorado: NREL, 1997 [11] J. Tangler, J.D. David Kocurek, Wind Turbine Post-Stall Airfoil Performance Characteristics Guidelines for Blade-Element Momentum Methods, NREL/CP-500-36900, Oct., 2004 [12] G. Xu, L.N. Sankar, Computational Study of Horizontal Axis Wind Turbines, Journal of Solar Energy Engineering [13] M.JOURIEH, “Développement d’un modèle representative d’une éolienne afin d’étudier l’implantation de plusieurs machines sur un parc éolien’’, Décembre 2007

CFSER 2010

Page 8