Performance prediction of a horizontal axis wind turbine using BEM and CFD methods

In this study, a combination of CFD (Computational Fluid Dynamics) and BEM (Blade Element Momentum Method) methods are used to simulate the flow field around a wind turbine rotor with horizontal axis. The main objectives are to predict the aerodynamic performances such as forces and torque imposed on the rotor blades which are essential to its structure or design. This approach requires much less computing time and memory than three-dimensional simulation flow around the wind turbine rotor with simple CFD method. The flow is assumed unsteady, incompressible and fully turbulent. This work consists of two parts: -The Calculation of aerodynamic coefficients by the CFD method, using Ansys / Fluent software. -The Simulation of 3D flow field through the rotor of the wind turbine using the BEM method. The obtained results are in good agreement with the experimental measurements.


Introduction
The study of aerodynamic performance of horizontal axis wind turbines has attracted many researchers in recent years because of potential applications in the production of electricity based on wind kinetic energy that falls into the category of clean renewable energy to the environment.
This study is based on CFD simulation (Computation Fluid Dynamic) using the software ANSYS/Fluent [1] that solves the equations of fluid mechanics in 2D and 3D flows, incorporating highly advanced turbulence models.This software uses the finite element method, appropriate digital techniques and features a mesh adapted to very complex configurations.
Many research studies have been conducted to determine the fine structure of the flow through the blades of horizontal axis wind turbines.As examples, we cite some work done in recent years, such as that of Chandrala et al. ( [2], 2012), which deals with the aerodynamic analysis of a wind rotor of the NACA 4402 type by comparing the calculation results with CFD measurements for different angles of attack, the Working Yilei He et al. ( [3] 2014) which focuses on optimizing the S809 airfoil, using a multi-objective algorithm MOGA, and the work of Bekhti et al. ([4], 2010) worked on the simulation of the airflow around a NREL S809 airfoil to predict the dynamic stall phenomenon and understand the physical behavior of the flow around a rotor wind.This work is based on solving the Navier-Stokes equations averaged RANS using advanced turbulence model.
In the present work, our study was focused on the numerical simulation using Ansys/Fluent software to understand and analyze the fine structure of the flow around wind turbine S809, and also determine the characteristics for aerodynamic optimization and the improvement of the net power.In addition, CFD simulation is used to calculate the aerodynamic forces and to identify the necessary correlations that will be incorporated into a global model for the calculation of aerodynamic performance equivalent wind turbines.The choice of a three-bladed wind rotor of the S809 type was made because of the availability of geometric and functional data provide a benchmark study.

CFD method 2.1 mathematical modelling
In this study, we are interested in the numerical simulation of the flow around airfoil wind turbine S809 exposed to the wind speed at a different angle of incidence α (Figure 1).The flow is assumed axisymmetric two-dimensional, unsteady and fully turbulent.The flow is therefore governed by the equations of conservation of mass and momentum, which are written respectively: Where t represents the time (in s), ρ is the density of the fluid, V is the free stream velocity of a fluid particle, p is the static pressure, τ is the viscous stress tensor, f is the resultant of the inertial forces acting on a fluid particle.The equation of state for a perfect gas is also added to express the density function of the static pressure.
The SST k-ω turbulence model is used in turbulent zones due to good prediction in separated flow simulation.The model uses the standard k-ω model near the wall, but switches to a k-ε model away from the wall.
For the flow analysis of the wind turbine blade mesh is created in the Ansys/workbench software.Figure 2 shows the meshing drawing of the airfoil below.

Mesh topology and boundary conditions
The mesh quality and domain size of the 2D CFD model do affect the computation accuracy and convergence time.A good mesh should be large enough to avoid boundary effects.There are many mesh topologies for airfoil CFD modeling.The most popular mesh topology is the C mesh which is designed to have a C-type topology around the airfoil as shown in figure 2.The mesh used here is the Ctopology structured mesh contains 181200 quadricelements.A large number of grids around the airfoil surface are used to capture the pressure gradient accurately at the boundary layer.In the far-field, the mesh resolution can become progressively coarser since the flow gradients approach zero.The dimensions of the domain of the following calculations: A height of 12 chord length, a length of 18 including 6 upstream of the trailing edge 12 and downstream.Momentum (Second order upwind)

Results analysis
In first part, we calculated the coefficients of lift and drag of S809 airfoil wind turbine.The calculation results are presented in Figures 4 and 5.In this simulation, the incidence angle α of the flow varies from 0 to 20°.
In Figure 4, we see that the curve is a straight line until α=9°, then undergoes a shift beyond this value.This inflection point is expressed by a detachment of the flow from the upper surface at the trailing edge airfoil.and 5 show the variations of the lift and drag coefficients as a function of incidence angle for S809 airfoil, these results are obtained for optimal mesh.We see that lift coefficient increases with incidence up to 8°.In this angle range, the curve is similar to a straight, after the lift increasing with incidence slowly, and then it undergoes a sharp fall: the static stall.
Figure 5 show that the drag coefficient sudden increase in the stall angle.These results correlated well with the experimental data.

BEM method
BEM theory is a compilation of both momentum theory and blade element theory.Momentum theory, which is for the determination of forces acting on the rotor to produce the motion of the fluid.This theory has no connection to the geometry of the blade, thus is not able to provide optimal blade parameters.Blade element theory determines the forces on the blade as a result of the motion of the fluid in terms of the blade geometry.By combining the two theories, BEM theory, also known as strip theory, relates rotor performance to rotor geometry.The mathematical code in this work for wind turbine design is based on BEM method.By applying axial and angular momentum conservation, the incremental axial force dT and angular torque dM on blade sectors can be obtained: Where ρ is the air density, V is the mean air flow velocity, r is the local radius, ω is the blade angular rotational speed, a and ' a are the axial and the tangential induction factors, respectively.
According to the blade element theory and by referring to Figure 7: Thus, the total torque (M) is the summation of dM for all the blades elements, and the wind turbine power is given by P=M.ω.Eventually, the wind turbine power prediction lies in solving the axial induction factor a and the tangential induction factor ' a .By combining Equations 3, 4, 5, and 6, one can solve the induction factors a and ' a as given below: The axial thrust coefficient C T is expressed as: By substituting the two induction factors into Equation 4, the wind turbine performance can be assessed.

Corrections for BEM method
The BEM method was developed by Glauert [5], it was not to give reliable results, and many authors have introduced corrections on the method.
These corrections relate to tip and Losses factor, to correct the detached from tip and hub effects of tourbillions to model the velocity distribution induced at the rotor.An approximate method of estimating the effect of tip losses has been given by Prandtl, given by equation ( 10 ) 1 ( 4 c a is approximately 0.2.

Analysis results
BEM method is programming in MATLAB, the results of the optimal distribution of the chord length and the twist angle are summarized in the graphs below.Figure 7 shows the optimum distribution of chord length airfoil S809 as a function of the local radius.In figure 7 and 8, the distribution calculated for the two parameters: the chord length and twist angle are not linear.In practice, these distributions are not respected because of manufacturing difficulties.The compromise would be to use linear distribution close to those calculated without sacrificing too much energy performance.In figure 9, we see that the distribution of the thrust force is an increasing linear function until r = 5.7, but the torque remains constant along the blade.

Conclusion
The S809 airfoil with different attack angles throughout the length is analyzed to find the maximum L/D ratio.This is done at the angle ranging from 0° to 20° for the velocity 14.6 m/s.The maximum C L /C D ratio is achieved at 6° of the attack angle.The airfoil optimal design, for a given site, must be based on the wind speed corresponding to the maximum average power of wind.
The BEM theory is very successful in horizontal axis wind turbine blade design.The modification factor and model were also combined into the BEM theory to predict the blade performance and there are good results of torque and thrust in each section between the improved BEM theory and CFD simulation.
The most difficult issues for the BEM theory are the mathematical representation of the correct lift and drag coefficient values, and the evaluation of axial and tangential induction factors.

Figure 3 .Table 1 :
Figure 3. Grid mesh around trailing edge For the boundary conditions, at the inlet, the velocity components are calculated for each case of angle of attack as follows.The x-component and y component of velocity is calculated respectively by D cos V x V

Figure 6 :
Figure 6: Blade section for a given radius r and C d are the lift and drag coefficients respectively, c is the local chord length.local solidity, here B is the number of blades and I is the angle of relative wind.
factor, R is the radius of blade.The application of this equation for the losses at the blade tips is to provide an approximate correction to the system of equations for predicting rotor performance and blade design, we find expressions of thrust force and torque: new formulas of axial and tangential induction factors and thrust coefficient becomes respectively:

Figure 7 .
Figure 7. Chord length distribution versus local radius The chord line decrease with increase of local radius r, and reach to maximum value at r=0.65 m, then decrease steeply.

Figure 8 .
Figure 8. Twist angle distributions and optimal attack angle versus a local radius Figure 8 represent the relationship between twist angle θ and local radius r, the curve shows the twist angle decrease toward blade's tip and there are large twists in regions of close to blade's root.In figure7and 8, the distribution calculated for the two parameters: the chord length and twist angle are not linear.In practice, these distributions are not respected because of manufacturing difficulties.The compromise would be to use linear distribution close to those calculated without sacrificing too much energy performance.

Figure 9 .
Figure 9. Distributions of thrust and torque versus a local radius