CC BY
6MODELING OF WIND TURBINES USING HYBRID MODELS
I.K. Dobrev*, F. Massouh**
ENSAM bvd. L'Hôpital, 151, Paris, 75013, France +33 1 44246256/+33144246266 *ivan.dobrev@paris.ensam.fr; **fawaz.massouh@paris.ensam.fr
Received: 24 Sept 2007; accepted: 25 Oct 2007
A hybrid model of kind "actuator surface" has been developed to represent the flow past a wind turbine rotor. The model uses a Navier-Stokes solver and permits calculation of rotor power and wake if the aerodynamic properties of the blade sections are known. The rotor geometry is simplified; the blades are replaced by their mean surfaces and a "pressure jump" boundary condition is applied to these surfaces. Thus, the proposed model is economic compared to full geometry simulation, because it is not needed to have fine grid around the blade. The hybrid model couples the Navier-Stokes solver with a blade element method (BEM). The solving process is iterative: at the beginning of each iteration a BEM determines the pressure discontinuities on the blade by means of rotor inflow and blade section performances. Then the CFD solver applies this pressure discontinuity in order to model the blade forces and calculate the flow past the rotor. The obtained velocity field is compared with results of previous iterations and if the required precision is attained, the calculation stops. The proposed hybrid model is tested in the case of a horizontal axis wind turbine (HAWT). The obtained results for rotor power and axial trust are satisfactory. Thus, this model can be employed for simulation of aerodynamic interaction between the wind turbines installed in the wind farm.
Keywords: wind energy equipment
Organization(s): Ecole Nationale Supérieure d'Arts et Métiers, Researcher. Education: Tech. Univ.-Sofia, Faculty of Energetic Machines (1978-1983). Experience: Tech. Univ.- Sofia, assistant (1983-2003). ENSAM, researcher (2004 to now). Main range of scientific interests: wind turbine, aerodynamics. Publications: 2 papers in international scientific journals.
Ivan Dobrev
Organization(s): Ecole Nationale Supérieure d'Arts et Métiers, Researcher, Assoc. Prof. Education: PhD - Paris - VI University (1984). Experience: ENSAM (1979 to now).
Main range of scientific interests: fluid mechanics, wind energy. Publications: 7 papers in international scientific journals.
Fawaz Massouh
Introduction
The proximity between the wind turbines installed in wind farm creates problems of aerodynamic interactions. Generally the wind farm development is complex problem and multiplicity of factors comes in play when the wind turbines are positioned. To optimize the energy production and the operation costs, engineers use
software tools developed especially for wind farm design. These software tools take into account wind turbine data, wind speed and direction, site topography, etc. However, in all cases it is needed to avoid the negative effect of aerodynamic interference between the wind turbines.
The simplified aerodynamic models used for wind farm design are not well adapted and cannot describe correctly
International Scientific Journal for Alternative Energy and Ecology № 5 (61) 2008
© Scientific Technical Centre «TATA», 2008
the behavior of wind turbine rotor. These models cannot obtain with sufficient precision the velocity field downstream the rotor and therefore they are not capable to evaluate the development of the wind turbine wake. To obtain numerical results with sufficient quality, a CFD simulation with an appropriate fine grid mesh is needed, but the solution is computationally very expensive. As result, modeling the interaction between more than two machines is impossible in practice. To reduce the computational cost, it is possible to use a simplified equivalent representation of the real rotor blades geometry that needs less grid points. This representation must be able to describe the behavior of the wind turbine rotor without modeling the exact blade geometry in the CFD computations. This kind of numerical modeling belongs to what is called hybrid modeling.
The hybrid models comprise two modules. In the first module a CFD solver computes the velocity field around the wind turbine rotor. Here, the presence of the rotor is modeled with source terms, pressure or velocity discontinuity. To prescribe these source terms or discontinuities, one second module uses a conventional method based usually on the blade element method (BEM). Here, the forces applied to the blades are calculated using the upstream velocity field and also the drag and lift coefficients of the blade sections. Thus, there is no need to model the flow around real blade geometry and the grid around the rotor may be coarsened. As results the need of computer power is reduced significantly.
In the field of wind turbine aerodynamics, [1] presents a comprehensive review on wake aerodynamics of wind turbines and several hybrid models are discussed. It is shown that many hybrid models use an actuator disk with the application of pressure or source terms. In these axisymmetric models, the source terms are distributed uniformly in the azimuthal direction and as a result the individual presence of blades is lost. To overcome this limitation and to represent more realistically the flow field downstream the rotor, a three-dimensional representation of the rotor blades is developed in [2]. In this model named as "actuator line", the geometry of real blades is replaced in CFD by source terms distributed radially along the blade axis. Here, the blade forces are determined by means of two dimensional airfoil data and the results of CFD computations are used to obtain the relative velocity and angle of attack. Compared to actuator disk, this model permits to represent individually each blade with its tip and root vortices and thus to improve rotor wake representation. The comparison of the actuator line with experimental data reveals the effectiveness of this proposed model in case of yaw and non-yaw compared to actuator disk model [3].
This paper is intended to develop the model of "actuator surface" proposed by authors in [4] and [5]. Compared to actuator line model the actuator surface model goes further in the blade representation. Here, each blade is
replaced by a surface of pressure discontinuity. The distribution of this discontinuity varies along the span but also along the chord. Thus the actuator surface model improve the blade representation and therefore the initial conditions of wake development compared to active line. Finally, to validate the proposed model, the results of hybrid calculation in the case of wind turbine will be presented and will be compared with experimental data.
Hybrid modeling
Hybrid models like actuator disc and actuators line, are presented in details in [3] and therefore no additional explanations are required. However, the active surface model has some differences and this paper is intended to present them. Actuator surface model like other hybrid models also combines a blade element method with a Navier-Stokes solver. In the CFD domain, the rotor geometry is simplified and the blades are replaced by thin surfaces. The specified boundary condition on these surfaces is "pressure discontinuity". Hence, the imposed surface forces replace the rigid blade wall and the number of nodes is significantly reduced, as there is no need to model the blade boundary layer. In the beginning of each iteration from the upstream velocity, the blade geometry and the airfoil data, the BEM module calculates the pressure jump distribution on the surface replacing the blade. Then the CFD module computes the flow velocity field, using as boundary condition the pressure distribution previously obtained from the BEM module. The solution is carried out iteratively, exchanging data between the BEM and CFD modules; it stops after convergence is reached. The calculation of pressure discontinuity is based on the blade element approach. At the blade radius r, the elementary forces acting in the normal and tangential directions on a blade element with span dr and chord c are:
and
dFn = 2 pW 2Cn (a)cdr
dFt = - pW 2Ct (a)cdr.
(1)
(2)
In the above formulas the force coefficients Cn and Ct are determined using the aerodynamic blade sections performances Cn = Cn(a) and Ct = Ct(a). The angle of attack a may be expressed as:
a = ф - ß,
(3)
where p is the blade section pitch angle and 9 is the angle between the plane of rotation and the reference relative velocity W. In the vortex line methods or BEM the flow angle is easy to evaluate because the axial induced velocity wia and tangential velocity w it are known explicitly:
Международный научный журнал «Альтернативная энергетика и экология» № 5 (61) 2008 © Научно-технический центр «TATA», 2008
9 = arctan
Vn
■ w,-
Qr-
(4)
where V0 is upstream velocity and Q is rotor angular velocity.
However, in actuator surface model the angle of attack a cannot be calculated explicitly because in the case of CFD modeling there is no means to separate the induced velocities in equation (4) from the rest of velocity field. Also the exact location, where flow angle must be calculated is not possible to define. Thus a different approach is needed to obtain the angle of attack. The flow around the wind turbine may be presented as the sum of a non-perturbed flow and another flow induced by the rotor blades. Then, the induced velocity field by the rotor also may be presented as a sum of two components:
- Local induced flow, created by the presence of the blade airfoils.
- Global induced flow, due to the presence of the rotor like an actuator that extracts kinetic energy from the wind and that slows down the velocity of the mass of air, which passes through the disk.
Upstream of the blades sections, at a distance of some chord lengths, the flow is slightly perturbed by the presence of local blade section. Therefore for the reference place, where the velocity and angle of attack must be determined, it is acceptable to use one line located upstream of the blades, where the flow is slightly perturbed by the presence of blade sections. However, this line must be sufficiently close to the rotor plane of rotation to have the same global induced field. Obviously, the proposed line is closer to the airfoil compared to the appropriate plane usually specified in "infinity". Therefore, for the same airfoil force coefficients the reference angles of attack are different and the airfoil performances must be corrected. However, this approach is very advantageous when airfoil performances are known from experiment in case of rotating blades [6] or from numerical simulations. In this case the normal and tangential force coefficients are obtained after normalizing the blade forces with velocities measured close to the blade. Thus the results may be used immediately without any correction if the same reference place is used.
In the BEM module the obtained normal force from equation (1) permits to calculated the force applied to blade element. Usually, in the case of HAWT, the blade section tangential forces are low compared to normal forces and therefore may be neglected. Then normal force is used to calculate the pressure discontinuity along the chord. In order to make the velocity induced by this discontinuity more adequate, it is preferable to use a chordwise distribution proposed in [4]. The pressure shape is close to the thin flat plate pressure distribution but at leading edge the pressure have no singularity. The pressure is linear in the intervals between leading edge and 1/4 chord and also between 1/4 chord and trailing
edge. The pressure value is 4p at the leading edge, p at the 1/4 chord and zero at the trailing edge. Hence, the normal force Fn is equal to pc and also the moment of pressure forces with respect to the point at 1/4 of the chord is equal to zero.
The hybrid model proposed, Fig. 1, is based on the CFD code Fluent 6.3 and the solution is obtained iteratively. In the CFD model, which computes flow field around the wind turbine, the blades are replaced by surfaces defined as "fan" boundary condition. This boundary condition corresponds to an imposed pressure difference between adjacent cells, located at the opposite sides of the boundary. Once at the beginning of the current iteration the CFD code executes a user-defined function UDF in C language. This function plays the role of the BEM solver and calculates the pressure distribution from blade geometry, rotor inflow and aerodynamic data of blade sections.
Fig. 1. Hybrid modeling
Fig. 2. Normal force coefficient along the blade radius
depending of local flow angle
The UDF function has access to all flow variables, which are needed to calculate the relative velocity vector along one reference line. Consequently, the normal force coefficient for all blade section along the span is interpolated from the blade geometry and airfoil data, Fig. 2. Then the resulting pressure is imposed on the
w,Ar
International Scientific Journal for Alternative Energy and Ecology № 5 (61) 2008
© Scientific Technical Centre «TATA», 2008
blade equivalent surfaces using the simplified pressure distribution. The solution is obtained after thousands of iterations depending of the number of nodes, demanded residuals, etc.
Numerical results
The studied wind turbine is the NREL Phase VI case. The turbine has two-bladed rotor of 10-meter diameter with blade sections of S809 airfoil [6]. The choice of NREL wind turbine was made because of its large experimental database obtained in NASA Ames wind tunnel. There are available inflow measurements at five blade radii for different upstream velocities between 5 and 25 m/s. For these five blade sections, there are also measurements of the pressure distribution both on suction and pressure side. These measurements permit to calculate for each upstream velocity the coefficients of normal force reduced with measured local velocity. For each upstream velocity corresponds different local flow angle (LFA) and using techniques of interpolation a two-dimensional function Cn = Cn (LFA, r/R) can be created, Fig. 2.
It must be noted that the reference velocity is measured only at 0.5 chords upstream of the blade section leading edge. Thus, the local perturbation caused by the blade is relatively important. However, in [4] it is shown that the velocities induced by a line with pressure discontinuity are very close to the velocities induced by an airfoil upstream of leading edge for distances greater than 0.5 chords. This similarity between the induced velocities permits to apply the obtained experimental relations between the normal force coefficient and local flow angle directly without any correction. During calculation in order to obtain the pressure distribution, for each node of pressure discontinuity surface the velocity, the LFA, the relative distance s/c from the leading edge and the relative radius r/R are determined, Fig. 3. Using the r/R and LFA, the force coefficients can be obtained from the experimental data, Fig. 2. Then using the distance s/c and the shape of chordwise pressure distribution it is possible to calculate the node value of pressure discontinuity that Fluent will apply as "fan" boundary condition.
The simulation model presents case "H" of the NREL test [6]. Here the rotor is in upwind position with blade tip pitch of 3° and average section Reynolds number about 1106. A cylinder that has a radius of 3.316R represents the flow field around the wind turbine rotor with radius R. It has also the length of 6R upstream and of 25R downstream the rotor plane. Here the value 3.316 is used because it leads to the same ratio as between the wind tunnel cross-section and the area of the wind turbine rotor. Hence, experiments and numerical simulations have the same coefficient of blockage. The actuator surface replacing the blade is represented by 4.000 nodes from the 500.000 ones used for the whole model. This surface is divided in the chordwise direction into 40 intervals, which are refined near the leading edge, where the pressure discontinuity gradient is strongest. In spanwise direction the blade is divided into 100 intervals equally spaced. The initial cells size in normal direction is 0.01 chords and the growth factor of 1.5 is used. To improve the wake calculation the mesh enclosed in a cylinder with length of 10R and diameter of 1.25R downstream the rotor is refined. The applied interpolation function, Fig. 2, for the normal force coefficients is based on the inflow measurements at five spanwise stations (r/R = 0.30, 0.47, 0.63, 0.80, 0.95) obtained for the case "H" [6]. To apply these data without any correction, the reference line passes through the same geometrical points where experimental coefficients are obtained.
Calculation is carried out iteratively. After numerous iterations, according to the number of the nodes and the required residual value, the convergence process is achieved. Usually, the computed rotor power reaches a constant value quickly, but additional iterations are needed to obtain the wake development. The comparison between wind turbine performances obtained experimentally and numerically by means of actuator surface is presented in Fig. 4 for the power and Fig. 5 for the axial trust. Some incertitude is involved because it is difficult to extrapolate the force coefficients near tip and root regions. To overcome this problem, more experimental or additional CFD results are needed. It must be noted that the numerical result for the power is close to measured values for upstream velocities up to 12 m/s. The disagreement in case of high velocities is due to the fact that the effect of "centrifugal pumping" is not modeled. In this case the flow in vicinity of the blade ends becomes highly three-dimensional and the hybrid model is non-adequate, similarly of all models that use the airfoil performances. However, it is useful to employ a correction for the circulation distribution near the blade tip and root region and also for the induced velocity angles. The comparison between the measured and the computed angles of local flow is shown in Fig. 6. The obtained results are good but there exist some discrepancy near the root, where the effect of "centrifugal pumping" is important.
Международный научный журнал «Альтернативная энергетика и экология» № 5 (61) 2008 *~j © Научно-технический центр «TATA», 2008
12
10
0) 6 o
4
-
)
" ; \
C ) o experiment - calculation
- , , ,
5 10 15
Velocity, m/s
Fig. 4. Wind turbine power
20
Fig. 5. Axial force
50
40
30
20
10
1 7m/s 11m/s -
A O O A 7 m/s 11 m /s 15 m/. exp exp exp
- -o— D A'
: -- - - o - - - - 0
: O - o- - .
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 r/R
Fig. 6. Local flow angle along the blade span Conclusion
In this paper the new hybrid model based on the actuator surfaces is proposed for simplified CFD calculation of the flow around wind turbine rotors. The objective is to validate the feasibility of this model to calculate the
aerodynamic performance of rotor. The rotor is simplified; the blades are replaced by thin surfaces constituted by the blade mean surfaces. Thus, the surface that replaces blade has the same pitch angle and chord as the original. On these surfaces, a pressure discontinuity is applied. This discontinuity is similar to thin flat plate chordwise pressure distribution, but without singularity in the leading edge point. The applied pressure discontinuity is calculated from aerodynamic properties of the blade section airfoils, the local flow angle and local velocity. Compared to active line model the present model improve the initial conditions of the wind turbine wake development because chordwise pressure is variable. Also the pressure discontinuity is comparable to that created by the real blade. This permits to create a velocity field similar to that around a blade section with the same normal force coefficient. In the case of wind turbine simulation the actuator surface model is capable to reproduce rotor mechanical power and forces, but for high wind velocity some disagreement with experimental results is revealed. This is due to fact that in case of high wind velocity the flow is detached. Therefore the flow is highly three-dimensional and the actuator surface model is non-adequate.
The suggested model here has two significant advantages compared to the other hybrid models. The first advantage is the possibility of using the 3-D airfoil data without applying any correction. The second advantage is that the velocity field downstream the blades is closer to reality compared to the models that use actuator disk or actuator line. Compared to the CFD methods, which use the complete three-dimensional rotor geometry, this model have the advantage of using a limited number of nodes. Hence, the size of the model is now suitable for studying wind farm design.
References
1. Vermeer L.J., Sorensen, J.N., Crespo A. Wind turbine wake aerodynamics // Progress in Aerospace Sciences. 2003. 39. P. 467-510.
2. Sorensen J., Shen W. Numerical modeling of wind turbine wakes // Journal of Fluids Engineering. 2002. 124. P. 393-399.
3. Mikkelsen R. Actuator disc methods applied to wind turbines. 2003 Ph.D. Thesis Technical University of Denmark.
4. Massouh F., Dobrev I., Rapin M. Numerical simulation of wind turbine performance using a hybrid model // AIAA paper 2006-0782, 44th AIAA, Reno, NV, Jan. 9-12.
5. Dobrev I., Massouh F., Rapin M. Actuator surface hybrid model The Science of making Torque from Wind // TWIND 2007, DTU-Lyngby.
6. Hand M., Simms D., Fingersh L.J. et al. Unsteady aerodynamics experiment phase VI: wind tunnel test configurations and available data campaigns. 2001. NREL/TP-500-29955.
CXTt - TATA —
LXJ
International Scientific Journal for Alternative Energy and Ecology № 5 (61) 2008
© Scientific Technical Centre «TATA», 2008
