 Open Access
 Total Downloads : 1366
 Authors : Mohammad A. Hossain, Ghizlane Zemmouri, Ziaul Huque, Raghava R. Kommalapati
 Paper ID : IJERTV2IS2161
 Volume & Issue : Volume 02, Issue 02 (February 2013)
 Published (First Online): 01032013
 ISSN (Online) : 22780181
 Publisher Name : IJERT
 License: This work is licensed under a Creative Commons Attribution 4.0 International License
Discrete Optimization of Wind Turbine Blade Airfoil
Mohammad A. Hossain1,2, Ghizlane Zemmouri 1,2, Ziaul Huque1,2, Raghava R. Kommalapati 2,3
1Department of Mechanical Engineering 2Center for Energy and Environmental Sustainability 3Department of Civil and Environmental Engineering
Prairie View A&M University TX 77446, USA
Abstract
This work is focused on optimization of wind turbine blade profile by two dimensional flow field analysis. A multiobjective response surface technique with computational fluid dynamics (CFD) was performed on two dimensional wind turbine blade airfoil. Based on the data of National Renewable Energy Laboratory (NREL) phase VI wind turbine rotor, six different airfoil (NACA xxxxxx, Fx and E Series) were used to calculate different aerodynamic loads (lift & drag) and their effects. Commercial software ANSYS Fluent was used to evaluate the lift coefficient (CL) and drag coefficient (CD) at different angle of attack () and three different Reynolds number (Re). Statistical code JMP was used to perform the response surface and finally a discrete optimization technique was developed to find an optimal airfoil that gives satisfactory performance in a wide range of design conditions.
Keywords Wind turbine, CFD, optimization, response surface

Introduction
According to the US Department of Energy the combustion of fossil fuels results in a net increase of
10.65 billion ton of atmospheric carbon dioxide every year [1]. The field of wind energy start to develop in 1970s after the oil crisis, with a large infusion of research money in the United States, Denmark and Germany to fine alternative resource of energy especially wind energy [2]. Blade is the most crucial part of wind turbine. Total performance, power output and efficiency depends on the design of the blade. Seeking a low cost, highly efficient blade design method has been an important problem required to be solved during wind turbine development. And the design of the blade completely depends on the selection of airfoil. The aerodynamic of general aviation airfoil has been fully studied in last few decades. The traditional wind turbine blade is using the aviation airfoil [3]. At present numerous research is focused on how to improve the
performance of the existing airfoil. Holten and Gyatt showed that using small flap in horizontal axis wind turbine could increase the output power [4,5]. Most wind turbine blades were adaptations of airfoils developed for aircraft and were not optimized for wind turbine uses. In recent years development of wind turbine blade airfoil has been ongoing. That may have modifications in order to improve performance for special application and wind conditions. To gain efficiency the blade should have both twisted and tapered. The taper, twist and airfoil characteristics should all be combined in order to give the best possible energy capture for the rotor speed and site conditions [6]. In this paper optimization of airfoil is focused. For this purpose six different airfoil is selected and different aerodynamic simulation is performed. Huque and Zemmouri has showed different optimum condition for six different airfoil [7] but they did not consider airfoils as a variable. In this work airfoils are considered as a discrete variable and with the help of MATLAB we successfully find out a single airfoil that gives the optimum aerodynamic performance in a wide range of design condition.

Airfoil Selection
Six different airfoil is selected for this work based on the previous work and literature. The six airfoils are NACA 63218, NACA 63421, NACA 64421, NACA 65421,FX63128,E387[7].
Figure.1. S809 Airfoil profile
For the comparison of the CFD data with the experimental results S809 airfoil is also considered [8]. These airfoil are created from the set of vertices generated from the University of Illinois at Urbana Champagne (UIUC) airfoil database [9]. These vertices are connected with a smooth curve creating the surface of the airfoil.

CFD Simulation

CFD Modeling
We considered three Reynolds number Re = 68,421, Re = 479,210, Re = 958,422 and a range of 00 to 120 angle of attack (). The CFD data of the 15 simulated cases for each airfoil were used to generate a response surface. The response surfaces were fit using standard leastsquare regression with quadratic polynomial using JMP. These response surfaces are obtained between design variable (Re, and Airfoil AF) and objective functions (CL and CD ) for each airfoil profile. All the design variable and the objective functions are normalized between 0 and 1 based on their maximum and minimum values in order to determine the response surface. Grid generation is done by ANSYS ICEM CFD algorithm. In this work approximately 86,000 unstructured quadrilateral elements Fig.2.were used to generate the mesh.
Figure 2. Mesh Domain
Figure 3. Mesh around airfoil
(a)
(b)
Figure 4. Mesh around (a) Leading edge, (b) Trailing edge
In order to have a stable and reliable solution, the mesh has minimum number of elements in the airfoil wall and grid points are clustered near the leading edge and trailing edge Fig. 3 in order to capture the flow separation and boundary layer of the airfoil wall.
Figure 5. Comparison between CFD data and NREL data of lift coefficient (CL) Vs angle of attack ()
In order to solve 2D Navierstokes equation, correct boundary condition plays very important role for appropriate results. In our model we considered noslip boundary condition in the wall and Inlet
velocity varies from 1ms1 to 7ms1 . Outlet pressure is considered as atmospheric pressure. Realizable k turbulence model along with second order upwind method is used in order to get more realistic result. To validate the CFD model we compared the experimental data of National Renewable energy laboratory (NREL) with our simulated data Fig.5.

CFD Result
Figure 6. shows a static pressure distribution of S809 airfoil at zero angle of attack and 4.7ms1 velocity. Fig.7 shows the velocity distribution of the same condition as the previous.
Figure 6. Static Pressure distribution of S809 airfoil at V= 4.7ms1, =0
Figure 7. Velocity around S809 airfoil at V=4.7ms1 and = 0
Fig. 8 shows the Cp distribution for NACA 63421 for the Reynolds number Re = 479,210. In the figure and for each angle of attack (uniform color), the bottom line represents the Cp distribution at the top surface of the airfoil, indicating lower pressure, and the top line represents the Cp distribution on the bottom surface of the airfoil indicating higher pressure. As the angle of attack increases from 0 to 12 for any Re, the area under the Cp curve increases indicating larger pressure difference between the bottom and the top surfaces. Similar trend is observed for different Re with the same
angle of attack. These are expected trends for any airfoil.
Figure 8. Cp distribution around NACA 63421 airfoil at Re = 479,210
Figure 9. Integrated Pressure Coefficient of NACA 63421 airfoil at different Reynolds number
Fig.9 represents the overall integrated pressure coefficient (Cp) as a function of angle of attack () of NACA 63421 airfoil at the three different Reynolds numbers. As expected, as we increase the angle of attack, the overall pressure coefficient increases for all six airfoils. However, within the same airfoil, Cp has little change as we move from a lower Reynolds number (Re = 68, 459) to a higher Reynolds number (Re = 95, 422).
The Cp of NACA 63218 airfoil increases continuously as we increase the angle of attack which indicates that it has not reached the stall condition yet, while the Cp plot of the other airfoils starts to flatten at around 11 to 12 of angle of attack which indicates that it is close to its stall condition. In addition, NACA 63218, NACA 63421, NACA 64421, and NACA 65421 airfoils
have small integrated Cp (Cp around 1.3 or 1.4) at stall condition which are much smaller than FX 63137 and
E387 airfoils (Cp around 1.6 or 1.8). Thus, we can conclude that the stall conditions could vary significantly between various airfoil profiles.
The response surface method fits an approximate function to a set of experimentally or numerically evaluated design data points [10]. There are various response surface approximation methods available in the literature. The polynomialbased approximations is being the most popular. In this technique, an appropriate ordered polynomial is fitted to a set of data points, such that the adjusted RMS error a is
adj
adj
minimized and quality parameter R2
as possible to one [11]. The a and R2
is made as close
are defined as
follows [7].
adj
Let N be the number of data points and let Np be the number of coefficients, and error ei at any point i is defined:
Figure 10. Lift coefficient of NACA 63421 airfoil as a function of angle of attack () at different Reynolds number (Re)
Figure 10. shows the plot of lift coefficient of NACA 63421 airfoil. CL is plotted as functions of angles of attack () and Reynolds number (Re). The general trends of all the plots are similar as expected; that is, CL increases with increasing and Re. Some of the observations from the plots are as follows.

The variations of CL between different Re are not significant.

The difference are more significant at higher
for NACA 65 421 and E 387 airfoil

CL for NACA 63218 at all Re did not reach the stall conditions; that is, the stall condition will be reached at much higher than = 12.

Both FX 63137 and E 387 indicate reaching stall condition at around = 12.

Both FX 63137 and E 387 show smaller variation with Re and reach higher values of CL = around 1.8 and 1.6 at = 12.
It is obvious from the previous observations that different airfoils behave differently with angle of attack and Reynolds numbers in the case of lift and drag coefficient (CL & CD). Hence there is an optimum combination of and Re for the maximum ratio of CL by CD for each airfoil. These optimal conditions are presented in the next sections.



Optimization Approach

Methodology
where fi a is the actual value of the function at the design
point and fi p is the predicted value. Hence,
Where,
adj
adj
adj
adj
The number of data N has to be greater than the number of coefficients Np so that the denominator of (2) is always positive and well posed. Since R2 needs to be as close as possible to 1 to represent a good fit, the terms in the numerator of (3) (a)2(NP 1) should be less than or equal to the denominator so that R2 will
always be positive.
In this study, the response surface method is applied with two objectives, namely, to generate response surface from the CFD simulation results and Reynolds number (Re), angle of attack () and airfoil (AF) are considered as design variable.

Response Surface

The CFD data of 15 cases were used to generate a response surface for each of the two objective functions for each airfoil shape. The response surfaces were fit using standard leastsquare regression with quadratic polynomial using JMP [12]. The following response surfaces for each of the objective function were obtained as a function of the three design variables (Re, , AF) of six airfoils combined:
Lift Coefficient response
CL= 0.0393+(0.0455*Re)+(1.2097 * )+ (0.5987*AF)
+(Re*Re*0.0200)+(Re**0.2175)+(alp*alp*0.3863)+ (Re*AF* 0.0803)+(*AF * 0.0290)+ (AF*AF*0.5279)
(5)
Drag Coefficient response
CD= 0.355+(0.260*Re) + (0.3616* )+ (0.1324*AF)
+(Re0.523)*(Re0.523)*0.0897+(Re0.523)*(0.5)*
(0.217)+(0.5)*(0.5)*0.2158+(Re0.523)*(AF
0.5833)*(1.152)+(0.5)*(AF0.5833)*(0.2269)+(AF
0.5833)*(AF0.5833)*(0.6689))
(6)
The quality of the response surface of this airfoil is shown in Table 2. The response surface for the entire objective had very high adjusted coefficient of both CL and CD which indicate good capabilities for this airfoil. Fig.1012 shows the response of different sets of variables with CL.
Table 1. Quality parameters of response surface of six airfoil combined
Observation 
CL 
CD 
R2 
0.95911 
0.96423 
adj 

Root Mean Square Error 
0.05803 
0.03248 
Mean of Response 
0.61775 
0.62735 
Figure 10. Response of angle of attack (alp) and airfoil (AF) in CL
Figure 11. Response of angle of attack (alp) and Reynolds number (Re) in CL
Figure 12. Response of Reynolds number (Re) and airfoil (AF) in CL
In the previous work [7], optimization is done by considering two design variable (Re & ). The result of that optimization is shown in Fig.13. After that we considered airfoil a discrete variable and optimized the objective function assigning the normalized value Table 2 of the Airfoil.
Figure 13. Optimum CL Vs angle of attack for all six airfoil
Table 2. Assigned value for the Airfoil variable
Airfoil Name 
Discrete value 
Normalized value 
NACA 63218 
1 
0.166 
E387 
2 
0. 333 
FX 63 128 
3 
0. 500 
NACA 63421 
4 
0. 666 
NACA 64421 
5 
0. 833 
NACA 65421 
6 
1.000 
Airfoil Name 
Discrete value 
Normalized value 
NACA 63218 
1 
0.166 
E387 
2 
0. 333 
FX 63 128 
3 
0. 500 
NACA 63421 
4 
0. 666 
NACA 64421 
5 
0. 833 
NACA 65421 
6 
1.000 
configurations and available data campaigns, Tech. Rep. NREL/TP50029955, NREL, 2001.
By using fmincon optimization tool of MATLAB we tried to maximize the objective function CL considering CD as an inequality constrain. After several iterative process the optimum value for airfoil we got is 0.825 which is very close to 0.833 or NACA 64421 airfoil. From that result we can conclude that NACA 64421 gives the optimum aerodynamic performance among all these six airfoil.
Acknowledgment
This work is supported by the National Science Foundation (NSF) through the Center for Energy and Environmental Sustainability (CEES), a CREST Center, award no. 1036593.
Reference

US Government, Department of Energy, US Department of Energy on Green House Gases, 2009 http://www.eia.doe.gov

Promod jain, "Wind Energy Engineering", McGraw Hill 2011, Chapter 1, pp 1,

Huang Jxioxg, "Research of New Wind Turbine Airfoils and Their Aerodynamic Performance", shantou, shantou university, 2001

Van. Holten , " Windmill with Diffuser Effect Induced
by Small Tip Vanes", Proc. Int. System Wind Energy Syst,
Cambridge, U.K.,1976

Gyatt G.W, Lissaman, P.B.S, " Development and
Testing of Tip Devices for Horizontal Axis Wind Turbine",
Aero Vironment Inc, DOE/NASA/03411, NASACR174991.
AVFR85/802

Anders Ahlstrom, " Aeroelastic Simulation of Wind Turbine Dynamics", PhD Thesis, Royal Institute of Technology, Department of Mechanics, Sweden, April 2005

Ziaul Huque, Ghizlane Zemmouri, Donald Harby,Raghava Kommalapati, "Optimization of Wind Turbine Airfoil Using Nondominated Sorting Genetic Algorithm and Pareto Optimal Front", International Journal of chemical Engineering,
Volume 2012, Article ID 193021

M. Hand, D. Simms, L. Fingersh et al., Unsteady aerodynamics experiment phase vi: wind tunnel test

University of Illinois at UrbanaChampaign Airfoil Coordinates Database, http://www.ae.illinois.edu/m selig/ads/coord database.html.

T. Goel, R. Vaidyanathan, R. T. Haftka,W. Shyy, N. V. Queipo, and K. Tucker, Response surface approximation of pareto optimal front in multi objective optimization, in Proceedings of the 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, pp. 2230 2245, Albany, NY, USA, September 2004.

R. H. Myers and D. C. Montgomery, Response Surface MethodologyProcess and Product Optimization Using Designed Experiment, WileyInter science, 1995.

JMP, The statistical discovery software, 19892002, Version 10, SAS Institute Inc., Cary, NC, USA.