Discrete Optimization of Wind Turbine Blade Airfoil

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 multi-objective 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 xx-xxxx, 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

1. 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.

2. Airfoil Selection

   Six different airfoil is selected for this work based on the previous work and literature. The six airfoils are NACA 63-218, NACA 63-421, NACA 64-421, NACA 65-421,FX63-128,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.

3. CFD Simulation

   1. 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 least-square 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 Navier-stokes equation, correct boundary condition plays very important role for appropriate results. In our model we considered no-slip boundary condition in the wall and Inlet

      velocity varies from 1ms-1 to 7ms-1 . 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.

   2. CFD Result

      Figure 6. shows a static pressure distribution of S809 airfoil at zero angle of attack and 4.7ms-1 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.7ms-1, =0

      Figure 7. Velocity around S809 airfoil at V=4.7ms-1 and = 0

      Fig. 8 shows the Cp distribution for NACA 63-421 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 63-421 airfoil at Re = 479,210

      Figure 9. Integrated Pressure Coefficient of NACA 63-421 airfoil at different Reynolds number

      Fig.9 represents the overall integrated pressure coefficient (Cp) as a function of angle of attack () of NACA 63-421 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 63-218 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 63-218, NACA 63-421, NACA 64-421, and NACA 65-421 airfoils

      have small integrated Cp (Cp around 1.3 or 1.4) at stall condition which are much smaller than FX 63-137 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 polynomial-based 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 63-421 airfoil as a function of angle of attack () at different Reynolds number (Re)

      Figure 10. shows the plot of lift coefficient of NACA 63-421 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 63-218 at all Re did not reach the stall conditions; that is, the stall condition will be reached at much higher than = 12.

      • Both FX 63-137 and E 387 indicate reaching stall condition at around = 12.

      • Both FX 63-137 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.

4. Optimization Approach

   1. 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.

   2. 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 least-square 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)

+(Re-0.523)*(Re-0.523)*0.0897+(Re-0.523)*(-0.5)*

(-0.217)+(-0.5)*(-0.5)*0.2158+(Re-0.523)*(AF-

0.5833)*(-1.152)+(-0.5)*(AF-0.5833)*(0.2269)+(AF-

0.5833)*(AF-0.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.10-12 shows the response of different sets of variables with CL.

Table 1. Quality parameters of response surface of six airfoil combined

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

configurations and available data campaigns, Tech. Rep. NREL/TP-500-29955, 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 64-421 airfoil. From that result we can conclude that NACA 64-421 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

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

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

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

4. Van. Holten , " Windmill with Diffuser Effect Induced

   by Small Tip Vanes", Proc. Int. System Wind Energy Syst,

   Cambridge, U.K.,1976

5. Gyatt G.W, Lissaman, P.B.S, " Development and

   Testing of Tip Devices for Horizontal Axis Wind Turbine",

   Aero Vironment Inc, DOE/NASA/0341-1, NASACR-174991.

   AV-FR-85/802

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

7. 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

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

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

10. 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.

11. R. H. Myers and D. C. Montgomery, Response Surface Methodology-Process and Product Optimization Using Designed Experiment, Wiley-Inter science, 1995.

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