Comparative Aerodynamic Analysis of Wind Turbine Blade Profiles

Comparative Aerodynamic Analysis of Wind Turbine Blade Profiles

1. M. M. Talukder*, M. K. Islam

   Department of Mechanical Engineering, Chittagong University of Engineering & Technology,

   Chittagong -4349, Bangladesh.

   M. R. Rukan

   Department of Mechanical Engineering, Bangladesh University of Engineering and Technology,

   Dhaka-1000, Bangladesh.

   Abstract Nowadays renewable energy has become the prime source of clean form of energy. Among various renewable energy sources, wind energy has become one of the key solutions to the prevailing energy crisis. To economically produce the maximum power by using wind turbine the aerodynamic performance of the blade profile must be assessed. The key parameters in analyzing the wind-turbine blade performance are lift and drag coefficient. In order to gain the maximum power from the wind turbine the maximum sliding ratio is desired. In this probe a comparative aerodynamic performance analysis of NRELs S819 and S821 airfoils were performed based on finite volume approach using a computational fluid dynamics (CFD) method. Changing the angle of attack and wind speed different aerodynamic parameters such as the drag coefficient, lift coefficient and pressure distribution over the airfoils were determined computationally. The computational results were verified experimentally by testing the wooden models of the airfoils in a subsonic open circuit suction type wind tunnel. The comparison with the experimental data indicate that the CFD approach applied in this investigation can precisely predict the aerodynamic characteristics of the wind- turbine blades.

   Keywords Lift; Drag; Sliding ratio; CFD; Airfoil

   1. INTRODUCTION

      Among all types of renewable and sustainable energy technologies wind power has the biggest future prospects because of its omnipresent and environment friendly nature [1]. Bangladesh being a tropical country does have a lot of wind flow at different seasons of the year and due to its high population density it has also huge thirst for electrical power. So wind energy can be a viable solution to this crisis. As wind energy is a low density source of power and its applications are growing rapidly in the world energy market, efficiency of wind energy construction must be assured. So it is crucial to maximize the efficiency of converting wind energy into mechanical energy. To achieve this goal from the aerodynamic perspective , the flow around and downstream of a wind turbine needed to be analyzed comprehensively and computational fluid dynamics offers the best options to do it. Without understanding the complete flow physics of Horizontal-Axis- Wind-Turbine (HAWT) further efficiency improvement is quite impossible. To achieve this aim and augment the efficiency, CFD can outplay the conventional aerodynamic approaches [2]. Thus, the recent years have seen the rise of numerical investigation on all HAWT aerodynamic features, performed on many different levels using available CFD techniques. Aerodynamic analysis of different wind-turbine-

      blade profiles using finite-volume method have been carried out by H. A. Kandil et al. [3]. The simulations were developed using the wind speeds over different sites in Egypt to find the suitable blade profiles for Egypts wind conditions. Also, 2D CFD-RANS simulations have been carried out at high Reynolds number (Re > 106) for the wind turbine blade profiles NACA 0008 and NACA0012 by Hoogedoorn et al. [4]. So it is evident that research articles available on HAWTs aerodynamic performance analysis mostly focused on the range of the high wind speeds and there are very little researches that focuses on the low-speed range available in Bangladesh. With its huge population density Bangladesh is currently facing an immense energy crisis where the average wind speed in the windiest places is about 5 m/s [5]. So it has paramount importance to conduct a research to find a suitable blade profile that will convert low speed wind energy available in Bangladesh in to electrical power in an efficient way. The maximum power developed by the wind turbine is desired in order to increase the efficiency.

      In the present study by analyzing the flow around the turbine blades the power is determined. Two different wind turbine blade profiles are selected from the profiles developed by the National Renewable Energy Laboratory (NREL). For both the airfoils lift and drag coefficients are determined at different angles of attack. The goal of this study is to represent the aerodynamic loads of the wind-turbine blade profiles in order to find the suitable profile for the wind conditions in Bangladesh. In order to determine the aerodynamic loads on the selected blade profiles, ANSYS commercial software was used. Finally, the simulation results were benchmarked against wind tunnel measurements. The comparisons show that the CFD code used in this study can precisely predict the wind-turbine blades aerodynamic loads up to a satisfactory level.

   2. BACKGROUND

      There are two forces and one moment that act upon an airfoil and these are lift force, drag force and pitching moment. Lift is the force used to overcome the gravity and is defined to be perpendicular to the direction of the oncoming airflow [1]. It

      Fig. 1. Airfoil nomenclature [6].

      Fig. 2. Resultant aerodynamic force and its components to which it splits [6].

      is formed as a consequence of the unequal pressure on the upper and lower surfaces of an airfoil. The higher the lift force the higher the mass that can be lifted off the ground.

      The drag force is defined as a force parallel to the direction of oncoming airflow. The drag force is due to both viscous friction forces at the surface of the airfoil and unequal pressure on the airfoil surfaces facing toward and away from the oncoming flow. For an airfoil, Hansen [1] stated that the lift to drag ratio should be maximized. As a result, it can improve efficiency when wind turbine generates electricity. Lift and drag coefficients CL and CD are defined as follows.

   3. COMPUTATIONAL APPROACH

      1. Geometry and CAD Model Preparation

         The 2D and 3D CAD representation of S819 and S821 airfoils are shown in Fig. 3 and Fig. 4, respectively. The co- ordinates of the airfoils were taken from NRELs airfoil database [7]. The co-ordinates of the airfoil were imported to SolidWorks to create the 2D geometry of the airfoils which was then imported to ANSYS design modeler for further processing such as incorporating the fluid volume around the airfoil and meshing.

         Fig. 3. 2D and 3D cad representation of S819 airfoil.

         = (1)

         = (1)

         Lift coefficient, CL

         FL 12c

         Fig. 4. 2D and 3D cad representation of S821 airfoil.

         2 0

      2. Mesh Generation

         The discretization of the fluid volume was done according

         12

         12

         Drag coefficient, CD = (2)

         to the principle of finite element analysis method. The mesh

         2 0 was generated using the mesh tools of ANSYS design modeler.

         The meshing strategy selected was an extruded 2D mesh in order to generate a 2D mesh of one cell thickness. Fig. 5 and

         Therefore, the sliding ratio can be defined as,

         =

         (3)

         Fig. 6, shows the generated mesh of S819 and S821 airfoils respectively. For S819 airfoil, the number of nodes were 29593 and number of elements were 29300 and for S821 airfoils the number of nodes were 29290 and number of elements were 29000. Accuracy depends on the number of elements and

         where, is the air density, c is the chord length of the airfoil, V0

         Flow separation occurs when the boundary layer travels far enough against an adverse pressure gradient that the speed of the boundary layer relative to the object falls almost to zero. The fluid flow becomes detached from the surface of the object and takes the forms of eddies and vortices. In aerodynamics, flow separation can often result in increased drag. Separation of boundary layer depends on Reynolds number. The higher the Reynolds number, there will be greater tendency that the flow is turbulent. Reynolds number is defined by,

         0

         nodes. With the increase in the number of nodes and elements the finite element analysis can be performed more accurately but the computation time will increase immensely. So a balance is required between the computation capability of the hardware and the number of elements used in the simulation process.

         Fig. 5. Meshing of S819 airfoil.

         =

         (4)

         where, 0 is the free stream velocity of the fluid, is the density of the fluid and is the dynamic viscosity of the fluid.

         Fig. 6. Meshing of S821 airfoil.

      3. Governing Equations

         NavierStokes equations arise from applying Newton's second law to fluid motion, together with the assumption that the fluid stress is the sum of a diffusing viscous term (proportional to the gradient of velocity) and a pressure term- hence describing the viscous flow. In an inertial frame of reference, the general form of the equation is:

         ( + . ) = + . + (5)

         where, v is the flow velocity, is the fluid density, p is the pressure, T is the stress tensor, f represents body forces (per unit volume) acting on the fluid and is the Del operator.

         Numerical flow simulation is performed by solving Navier- Stokes equations, which are formulation of mass, momentum and energy conservation laws.

         The Spalart Allmaras turbulence model solved a modelled transport equation for kinematic eddy viscosity without calculating the length scale related to the shear layer thickness. The variable transported in the Spalart Allmaras model is which is assimilated, in the regions which are not affected by strong viscous effects such as the near wall region, to the turbulent kinematic viscosity. This equation has four versions, the simplest one is only applicable to free shear flows and the most complicated, which is written below, can treat turbulent flow past a body with laminar regions [8].

         (6)

         This transport equation brings together the turbulent viscosity production term, and the destruction term, . The physics behind the destruction of turbulence occurs in the near wall region, where viscous damping and wall blocking effects

         are dominants. No heat generation or transfer is considered. The other terms or factors are constants calibrated for each physical effect which needs to be modelled. This equation allows to determinate for the computation of the turbulent viscosity, from:

         = 1 (7)

         The cell centered finite volume method (FVM) is used to obtain conservative approximations of the governing equations on the locally refined rectangular mesh. The governing equations are integrated over a control volume which is a grid cell and then approximated with the cell centered values of the basic variables. The integral conservation laws may be represented in the form of the cell volume and surface integral equation:

         + . = (8)

         This is replaced by,

      4. Simulation

      At first the mesh was done for the fluid volume around the airfoil. For the prediction of wind turbine aerodynamics [9], the k-epsilon method had been chosen to capture the turbulence. The k-epsilon model was chosen for accurate boundary layer detection due to its ability to capture the influence of different factors that affect transition such as the free-stream turbulence and pressure gradients. For controlling and reducing the numerical solution errors, the upwind scheme method was selected. In near wall regions, predicting the velocity gradients produced by boundary layer phenomena need elements with high aspect ratios. The boundary conditions were applied to inlets, outlets, walls. It should be noted that the flow was always subsonic due to the low range of velocity used in the study. The boundary condition at the airfoil surface had been set to no-slip- solid wall boundary. Inlet turbulence was assumed with a medium intensity. Pressure boundary conditions were applied at the domain outlets, and the average static pressure method was used in order to allow the pressure to vary locally on the boundary.

   4. EXPERIMENTAL APPROACH

      1. Model Fabrication for Wind Tunnel Testing

         Wooden models of S819 and S821 airfoil were made for wind tunnel testing shown in Fig. 7, with the predefined chord length (200mm) and span length (120mm). Beech wood was used to construct the model.

         Fig. 7. Wooden models of S819 and S 821 airfoil.

         Inside of the models were kept hollow to incorporate the pressure tubes and twelve pressure tapings were introduced along the center line at the upper and lower surfaces of the airfoil maintaining equal distance with each other. Individual pressure tubes were glued with each hole and were connected with the multi tube manometer for the measurement of the airfoil surface pressure. In order to hold the model inside the wind tunnel test section and connect it with the load measuring device a cylindrical support of stainless steel was provided with the airfoil models.

      2. Wind Tunnel Specification

         The experiment was carried out in a 380mm×355mm open circuit wind tunnel as shown in Fig. 8, available at compressible fluid mechanics laboratory of Department of Mechanical Engineering, CUET. The wind speed is created by rotating a fan. The fan motor was powered by 220V-3-50Hz power supply through motor speed controller. Thus the wind speed in the tunnel could be varied by controlling the fan motor speed. Air was sucked through the test section by the fan located at the rear of the tunnel. A honeycomb was employed in front of the test section to reduce flow turbulence.

         ()

         + . =

         (9)

         Fig. 8. Open circuit suction type wind tunnel.

         TABLE I. WIND TUNNEL SPECIFICATIONS [10]

         Type	Open circuit wind tunnel (S HUNT)
         Test section	380 × 355 × 330mm
         Fan	Compatible capacity
         Motor Type	DC Motor 220v,10A,2800 rpm
         Speed controller	For variable speed
         Air velocity	Maximum 5.3 m/s (in test section)
         Number of Tube in Manometer	20
         Fluid inside the Manometer	Kerosene

      3. Test conditions and data collection procedures

      The flow of air was considered incompressible throughout the experiment. All the data of the experiment were taken at the room temperature of 25 but flow velocities were varied with the help of speed regulator to 3 m/s, 4 m/s and 5 m/s. The angle of attack was varied from -4° to 16° with an increment of 4°. The static pressure for different flow velocities at different AOA were measured from both the upper and lower surfaces of the airfoil models through pressure tapings by using a multi-tube kerosene manometer. Specific density of both air and kerosene corresponding to room temperature was assumed to be 1.145 kg/m3 and 790 kg/m3 respectively.

      Fig. 9. Experimental setup.

   5. RESULTS AND DISCUSSION Computations were performed for S819 airfoil to test the

      grid independency. So, for a S819 airfoil a finer mesh, made of 29300 quadrilateral cells, 29593 number of nodes was compared with a mesh of 8670 quadrilateral cells and 18981 numbers of nodes with the same input of AOA = 80 and M = 0.015. The lift co-efficient found for the finer mesh is 1.40798 and varies with that of the other mesh by 0.75% which is concurrent to the independency test.

      In Fig. 10 and Fig. 12, we find the coefficients of drag and lift for S819 airfoil for zero degree angle of attack. After 74th iterations convergence was obtained and the values of CL and CD as 0.15386 and 0.022071 were found respectively. On the other hand for S821 airfoil convergence was obtained after 80th iterations and the values of CL and CD were found as 0.21763 and 0.030105 respectively for zero degree angle of attack which is represented in Fig. 13 and Fig. 11, respectively. For a certain angle of attack at the initial stage, the coefficient value was changing with the increase of iteration and later the value was stabilized to a fixed value. From simulation, the coefficient of lift CL and coefficient of drag CD was calculated for different angle of attacks for both the airfoils at different flow velocities.

      Fig. 10. Coefficient of drag of S819 airfoil for zero degree angle of attack.

      Fig. 11. Coefficient of drag of S821 airfoil for zero degree angle of attack.

      Fig. 12. Coefficient of lift of S819 airfoil for zero degree angle of attack.

      Fig. 13. Coefficient of lift of S821 airfoil for zero degree angle of attack.

      Fig. 14. Pressure contour of S819 and S821 airfoils for zero degree angle of attack at 5 m/s flow velocity.

      In Fig. 14, the pressure distribution over the S819 and S821 airfoil is shown. It is evident from the figure that the maximum pressure is generated at the leading edge of both the airfoils. Another thing is clear from the pressure contour is that negative pressure is created at the top and bottom surface of the airfoils. The imbalance between these two negative pressures is the reason for the lift force which is the resultant of these two forces. With the increase of angle of attack up to a certain limit the negative pressure at the top surface increases and the negative pressure at the bottom surface decreases so the lift increases. But after certain angle of attack the opposite phenomenon occurs due to the separation of flow from the upper surface of the airfoil and formation of wake near the trailing edge which ultimately decreases the lift and increases the drag.

      Fig. 15. Velocity distribution over S819 and S821 airfoils for zero degree angle of attack at 5 m/s flow velocity.

      From Fig. 15, it is observed that the streamlines over the airfoils started to detach from the upper surface as the air passes over the trailing edge of the airfoil. This is the region where small vortices are created and thus it creates a small amount of negative pressure. These vortices and negative pressure consequently induces wakes near the separation point of the streamlines over the airfoil which increases the drag.

      Flow separation begins to occur at small angles of attack while attached flow over the wing is still dominant. As angle of attack increases, the separated regions on the top of the wing increase in size and hinder the wings ability to create lift. At the critical angle of attack, separated flow is so dominant that further increases in angle of attack produce less lift and vastly more drag.

      Fig. 16. Computational lift co-efficient of S819 and S821 airfoil with respect to angle of attack.

      Fig. 17. Experimental lift co-efficient of S819 and S821 airfoil with respect to angle of attack.

      From the computational lift coefficient data shown in Fig. 16, it is obvious that with the increase of angle of attack the lift is increasing for both the airfoils. The maximum lift for S819 is 1.79246 and for S821 is 1.60134 both obtained at 120 angle of attack which is termed as the critical angle of attack. The decrease in lift begins after the critical angle of attack which is termed as stall. Analogous pattern of lift coefficient is observed in experimental data collected from the wind tunnel testing of the airfoils. The maximum lift for S819 is 1.3546 and for S821 is 1.2051 obtained at 100 and 120 angle of attack respectively.

      Fig. 18. Computational drag co-efficient of S819 and S821 airfoil with respect to angle of attack.

      Fig. 19. Experimental drag co-efficient of S819 and S821 airfoil with respect to angle of attack.

      From Fig. 18 and Fig. 19, it is seen that the computational drag coefficient for both the airfoils increases with the increase of angle of attack. The drag coefficient of S819 airfoil is always higher than that of S821 airfoil for any particular angle of attack. Similar trend of drag coefficient is observed in the experimental data which confirms the validity of the investigation.

      From the computational sliding ratio versus angle of attack curves shown in Fig. 20 and Fig. 21, it is observed that the sliding ratio slightly increases as the wind speed increases at the same operating AOA for both the airfoils. The maximum sliding ratio of S819 airfoil is at angle of attack 2° while for the S821 airfoil the maximum sliding ratio occurs at 0° angle of attack. The maximum sliding ratio of S821 airfoil is greater than that of S819 airfoil.

      Fig. 20. Computational sliding ratio versus AOA curve of S819 airfoil for different wind speeds.

      Fig. 21. Computational sliding ratio versus AOA curve of S821 airfoil for different wind speeds.

      It is concluded that to achieve the maximum power from the wind turbine, the operating AOA should lay between 00 and 20 for both the profiles.

      Furthermore, for the S819 profile the best operating AOA is positive because it is almost symmetric profile and the maximum thickness occurs between 40% and 50% of the chord. Also, for the S821 profile the best operating AOA is negative because the maximum thickness occurs between 20% and 40% of the chord. The best profile at each wind speed can be selected according to the maximum sliding ratio from the optimum range of the AOA [3].

      Finally, the results in Fig. 20 and Fig. 21, show that at these wind speed, the S821 profile gives the maximum power as it has the maximum value of the sliding ratio. This profile is strongly recommended to be used for the wind turbine blade operating at the wind speeds under consideration.

      The simulation result was verified with experimental data collected from the wind tunnel testing of the airfoils. According to the Fig. 22 and Fig. 23, it is observed that the sliding ratio increases with the increase of the wind speed for both the airfoils for a certain angle of attack. At the beginning sliding ratios of both the airfoils were increasing with the increase of the angle of attack. In this trend the maximum sliding ratio of S819 airfoil was found at 20 angle of attack and the maximum sliding ratio of S821 airfoil was found at 00 angle of attack which is analogous to the simulation result. The maximum sliding ratio of S821 airfoil was also higher than that of S819 airfoil. After reaching the maximum value of the sliding ratio there was a gradual fall of the value of the sliding ratio with the increase of the angle of attack. The computational sliding ratio was higher than the experimental sliding ratio for both the airfoils which confirms the significance of the mesh quality i.e. selection of the shape, size and also the number of elements in the mesh.

      Fig. 22. Experimental sliding ratio versus AOA curve of S819 airfoil for different wind speeds.

      Fig. 23. Experimental sliding ratio versus AOA curve of S821 airfoil for different wind speeds.

      The lack of calibration facilities for the wind tunnel and different atmospheric conditions of air can also play a significant role in the data collection and evaluation process which ultimately could results in the discrepancy between the experimental and simulation data.

      The better airfoil always has a higher lift to drag ratio when it is compared with other airfoils. In this case, S821 airfoil demonstrated higher lift to drag ratio than S819 airfoil. Soit is explicit that S821 airfoil is best suited for aerodynamic applications than S819 airfoil with higher lift to drag ratio and less wake generation.

   6. CONCLUSIONS

It can be concluded from the investigation that sliding ratio plays a key role in determining the effectiveness and usability of a wind turbine blade profile. Beside this the AOA has a predominant effect on determining the optimum profile while the wind speed does not affect the optimum profile. The optimum operating AOA should lie between 20 and 00 to get the maximum sliding ratio and the maximum power extracted from the wind by using S819 and S821 airfoils as at these angle of attacks the airfoils have maximum sliding ratios. For the same operating wind speed the S821 has better sliding ratio than that of the S819. So for low speed application in a country like Bangladesh S821 is more suitable than S819. It should also be mentioned that as the AOA increases above the optimum range the sliding ratio decreases and the difference in the sliding ratios between the profiles becomes very small.

ACKNOWLEDGEMENT

The authors would like to express their sincere gratitude to the Department of Mechanical Engineering, CUET for providing the wind tunnel facilities and other technical supports to carry out the research activity.

REFERENCES NOMENCLATURE

1. Hansen, M.O.L., 2008. Chapter 2:2-D Aerodynamics, Aerodynamics of wind turbine Second Edition, Earthscan, UK and USA.

2. Johansen, J., 1997. Prediction of Laminar/Turbulent Transition in Airfoil Flows, Journal of Aircraft, Aerospace Research Central, Denmark, Volume 36, No. 4, pp. 731-734.

3. Mohamed, A. S., Kandil, H. A., Shaltot, A., 2012. Aerodynamic analysis of different wind-turbine-blade profiles using finite- volume method, Volume 64, pp. 541-550

4. Hoogedoorn, E., Jacobs, G.B., Bey, A., 2010. Aero-elastic behavior of a flexible blade for wind turbine application: a 2D computational study, Elsevier Ltd. J Energy; 35(2):77885.

5. Khan, M.G. K., Rahman,T. and Alam, M.M., 2004. Wind energy in Bangladesh: prospects and utilization initiatives, Third ICECE 2004, 28-30,December 2004, Dhaka, Bangladesh

6. Anderson, J. D., 2011. Fundamentals of Aerodynamics, Fifth Edition, McGraw-Hill Companies, Inc.

7. https://wind.nrel.gov/airfoils/Documents/S819,S820,S821_Design Date of access : 23/05/2015

8. Azim, R., Hasana, M. M., Alib, M., 2014. Numerical investigation on the delay of boundary layer separation by suction for NACA 4412, 6th BSME International Conference on Thermal Engineering (ICTE 2014).

9. Hossain, S., 2014. A comparative flow analysis of NACA 6409 and NACA 4412 airfoil, IJRET, Volume 3, Issue 10 pp. 342-350.

10. Pal, S., Kabir, S.M.H., Talukder, M. M. M., 2015. Aerodynamic analysis of a concept car model, Proceedings of the third ICMERE, 26 29, November 2015, Chittagong, Bangladesh.

11. Prince, S.A., Khodagolian, V., Singh, C., Kokkalis, T., 2009. Aerodynamic Stall Suppression on Airfoil Sections Using Passive Air-Jet Vortex Generators, AIAA Journal, India, Volume 47, No. 9, pp. 2232-2242.

12. Shih, T.H., Liou, W.W., Shabbir, A., Zhu, J., 1995. A New k- eddy Viscosity Model for High Reynolds Number Turbulent Flows, National Aeronautics and Space Administration, USA, pp. 1-32.

13. Bensiger, S.S.B., Prasanth, N., 2012. Analysis of Bi-Convex Airfoil Using CFD Software at Supersonic and Hypersonic Speed, Elixir International Journal, Elixirpublishers, India, Volume 53, pp. 11695-11698

14. Sutherland, W., 1893. The Viscosity of Gases and Molecular Force, Encyclopaedia of Australian Science, Australia, Volume 36, pp. 507-531.

15. Kevadiya, M., Hemish, A.,Vadiya, 2013. 2D Analysis of NACA 4412 Airfoil, International Journal of Innovative Research in Science Engineering and Technology, India, Volume 2, No. 5 , pp. 168-169.

16. Mashud, M., Ferdous, M., Omee, S. H., 2012. Effect of Spoiler Position on Aerodynamic Characteristics of an Airfoil, International Journal of Mechanical and Mechatronics Engineering, Bangladesh, Volume 12, No. 6.

17. Sayed, M.A., Kandil, H.A., Morgan, EII. , 2011. Computational fluid dynamics study of wind turbine blade profiles at low reynolds number for various angles of attack. , The 4th international meeting on advances in thermo fluids, Melaka, Malaysia.

18. Cao, H., May 2011. Aerodynamics Analysis of Small Horizontal Axis Wind Turbine Blades by Using 2D and 3D CFD Modeling

    M. Sc. Thesis, University of the Central Lancashire, Preston, England.

19. Gulzar, O., Gulzar, S., Bhatele, S., Soni, N., 2014. Impact of Variation in Angle of Attack on NACA 7420 Airfoil in Transonic Compressible Flow Using Spalart-Allamaras Turbulence Model, IJRMET, India, Volume 4, No. 2, pp. 35-39.

20. Launder, B.E., Spalding, D.B., 1974. The Numerical Computation of Turbulent Flows, Computer Methods in Applied Mechanics and Engineering, Science Direct, UK, Volume 3, No. 2, pp. 269-289.