Table of Contents Author Guidelines Submit a Manuscript

A corrigendum for this article has been published. To view the corrigendum, please click here.

Journal of Renewable Energy
Volume 2013 (2013), Article ID 653103, 9 pages
Research Article

Propagation of Shock on NREL Phase VI Wind Turbine Airfoil under Compressible Flow

1Department of Mechanical Engineering, Prairie View A&M University, Prairie View, TX 77446, USA
2Center for Energy and Environmental Sustainability, Prairie View A&M University, Prairie View, TX 77446, USA
3Department of Civil and Environmental Engineering, Prairie View A&M University, Prairie View, TX 77446, USA

Received 19 March 2013; Revised 1 July 2013; Accepted 1 July 2013

Academic Editor: Cheng-Xian Lin

Copyright © 2013 Mohammad A. Hossain et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The work is focused on numeric analysis of compressible flow around National Renewable Energy Laboratory (NREL) phase VI wind turbine blade airfoil S809. Although wind turbine airfoils are low Reynolds number airfoils, a reasonable investigation of compressible flow under extreme condition might be helpful. A subsonic flow (mach no. ) has been considered for this analysis and the impacts of this flow under seven different angles of attack have been determined. The results show that shock takes place just after the mid span at the top surface and just before the mid span at the bottom surface at zero angle of attack. Slowly the shock waves translate their positions as angle of attack increases. A relative translation of the shock waves in upper and lower face of the airfoil are presented. Variation of Turbulent viscosity ratio and surface Y+ have also been determined. A k-ω SST turbulent model is considered and the commercial CFD code ANSYS FLUENT is used to find the pressure coefficient (Cp) as well as the lift (CL) and drag coefficients (CD). A graphical comparison of shock propagation has been shown with different angle of attack. Flow separation and stream function are also determined.

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] which has an enormous impact on environmental imbalance. As a result more focus on conversion of energy from alternate source has been given for the last few decades. In near future wind will be the most reliable green energy in the history of mankind. The field of wind energy started to develop in 1970s after the oil crisis, with a large infusion of research money in the United States, Denmark, and Germany to find alternative resource of energy especially wind energy [2]. To design the blade of a wind turbine proper assessment of aerodynamic characteristics of airfoil plays the most important role. The most effective way to design the blade is to have accurate experimental data set for the correct airfoil. But such data set are not always available and the designer must rely on calculated data such as simulated data generated by large-scale CFD code.

Recent applications of CFD to solve the Navier Stokes equations for wind-turbine airfoils are reflected in the works of Chang et al. [3]. They used their in-house code to solve the 2D flow field around the S805 and S809 airfoils with attached flow and the S809 airfoil with separated flow. Computations were made with the Baldwin-Lomax, Chein’s low-Reynolds-number k-ε [4], and Wilcox’s low-Reynolds-number k-ω turbulence models [5]. Unsteady compressible flow over airfoils was extensively studied by Chen et al. [6]. They studied various fundamental flow mechanisms dictating the intricate flow phenomena, including moving shock wave behaviors, turbulent boundary layer characteristics, kinematic of coherent structures, and dynamical process in flow evolution as shown in Figure 1.

Figure 1: Contour of turbulent kinetic energy studied by Li Chen et al. [6].

They also studied the moving shock wave characteristics and moving shock wave-boundary layer interaction. Shock generates sudden fluctuation of pressure and velocity. Detailed study of compressible flow has also been studied by Tijdeman and Seebass [7]. They characterized three types of moving shock waves. Lee [8] investigated self-sustained shock wave on airfoil experimentally and proposed a feedback model to estimate the frequency of oscillation of the shock motion. Wang et al. [9] studied unsteady boundary layer separation and vortex shedding in the trailing edge region for high-intensity turbulent flow. Marvin et al. [10] and Rumsey et al. [11] performed numerical simulation using the time-dependent two-dimensional Reynolds-average Navier Stokes (RANS) equation with turbulence model. Deck [12] studied a zonal detached eddy simulation (DES) method to predict buffet phenomenon on a super critical airfoil. Interaction between turbulent boundary layer and sustained moving shock was numerically studied by Smits and Dussauge [13]. Recent advancement on numeric simulation has been developed by Spalart [14]. His explicit study on DSE has become a powerful tool investigating high Reynolds number compressible flow. Moreover, three-dimensional Favre-averaged compressible Navier-Stokes are solved numerically by Lu et al. [15] where he used finite volume method with the combination of shock capture technique.

In this paper the aerodynamic characteristics of wind turbine airfoil (S809) under compressible flow condition have been studied because, to the best of the author’s knowledge, very little work has been done in this field of wind turbine blade airfoil due to lack of available experimental data. In recent years development of wind turbine blade airfoil has been ongoing and has many modifications in order to improve performance for special application and wind conditions. To gain efficiency the blade should have both twist, and taper. 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 [16]. Moreover, the blade length of a wind turbine blade is increasing day by day. In recent years NREL 5MW offshore wind turbine has blade length of 90 m. In near future the length of the blade will increase more and more. For a 100 m blade rotating at 25 rpm the tip speed will approximately reach a Mach no. of 0.78. Therefore a compressible flow analysis might be helpful for the designer in order to avoid unexpected incidents.

2. k-ω SST Turbulent Model

The SST k-ω turbulence model (Menter 1993) is a two equation eddy-viscosity model. SST k-ω model can also be used as a low-Re turbulent model without considering any extra damping function [17]. This model can produce a large turbulence levels with regions of large normal strain like stagnation region and regions with strong acceleration [18]. The original k-ω model can be defined as

The Shear Stress Transport (SST) formulation combines the two equations. The shear stress boundary layer and kinematic eddy viscosity can be defined as

3. Airfoil Selection

National Renewable Energy Laboratory (NREL) has developed different airfoils specially for horizontal axis wind turbine [19]. Some of the airfoils are S801, S805, S809, S8012, and so forth. Among them we considered S809 as this airfoil was used in NREL phase VI wind turbine experiments. The airfoil for simulation is created from the set of vertices obtained from the University of Illinois at Urbana Champagne (UIUC) airfoil database [20]. These vertices are connected with a smooth curve creating the surface of the airfoil as shown in Figure 2.

Figure 2: S809 airfoil profile.

4. CFD Simulation

4.1. CFD Modeling

We considered a subsonic flow where represents the mach number and a range of 0° to 10° angle of attack (). Grid generation is done by ANSYS ICEM CFD algorithm. In this work approximately 0.2 million unstructured triangular elements were used to generate the mesh. Computational domain consists of a smooth parabola for better resolution of results as shown in Figure 3. The maximum domain length is 125 m and the width is 90 m.

Figure 3: Mesh domain.

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 (Figures 4 and 5) in order to capture the flow separation and boundary layer of the airfoil wall. The simulation is modeled in such a way that it would be mesh independent. The quality of the mesh largely depends on the orthogonal quality parameter. This orthogonal quality varies from 0 to 1. Values close to 0 correspond to low-quality mesh [21]. In this CFD simulation the orthogonal quality is 0.95.

Figure 4: Mesh around airfoil.
Figure 5: Mesh around trailing edge.

A pressure-based solver is set and ideal gas approximation is considered for all the CFD simulations. In order to solve 2D Navier-stokes equation, correct boundary condition plays a very important role for appropriate results. k-ω SST turbulent model with no slip boundary condition at the wall has been considered. Outlet pressure is considered as atmospheric pressure. Coupled second-order upwind method is used as a solving method. The turbulent viscosity ratio is considered 10% and operating temperature is assumed 300 K. The operating condition is zero gage pressure or an absolute pressure of 101325 Pa. Sutherland’s viscosity law which is the relation between the dynamic viscosity (μ) and the absolute temperature (T) is considered. Sutherland’s law is based on kinetic theory of ideal gases and an idealized intermolecular-force potential [22] which is being used for many advanced CFD simulations. ANSYS FLUENT is used with Semi-Implicit Method for Pressure Linked Equation (SIMPLE) solution method. It is a steady state iterative technique. Here velocity field is obtained by solving the momentum equation, and pressure distribution is calculated by the initial condition. SIMPLE solver algorithm resolves the pressure-velocity coupling [23]. To validate the model, a simple incompressible flow is considered at a velocity of 4.45 m/s at an angle of attack of 5.12°, and the pressure distribution around the airfoil is determined. The pressure distribution has been compared with the experimental data as shown in Figure 6. The experimental data are collected from Delft University low-turbulence wind tunnel data. [24]. Simulation results show good agreement with the experiment. Although the analysis is for compressible flow, there is no such experimental result for compressible flow over wind turbine. But ample amount of experiments have been conducted for incompressible flow over wind turbine blade. Because of the availability of the experimental data, the CFD model is validated for incompressible flow, and then the velocity is increased successively in order to resolve the compressible flow.

Figure 6: CFD model validation at velocity 4.45 m/s and .
4.2. CFD Result

The main objective of this study is to find the flow behavior and the shock propagation around the airfoil in compressible flow condition. In order to do that the static pressure, the mach number, the turbulent viscosity, and the temperature variation around the airfoil have been determined. The coefficient of pressure distribution around the airfoil and the lift and drag () coefficients at different angles of attack have also been determined.

According to NREL Phase VI wind turbine blade, A 5 m blade rotating at 70 rpm would have a flow approximately mach 0.11 at the tip for the corresponding upstream wind velocity of  m/s. Figure 7(a) shows a static pressure contour of S809 airfoil in monochromic form at incompressible flow condition with velocity  m/s (mach 0.11) and . On the other hand Figure 7(b) shows the pressure contour for the same airfoil and same angle of attack but in compressible flow (mach 0.8). Figure 7(c) also shows the same pressure contour as in Figure 7(b) but in color to illustrate the shock generation at the top and the bottom surfaces of the airfoil.

Figure 7: (a) Pressure contour at mach 0.11 and , (b) pressure contour at mach 0.8 and , and (c) static pressure distribution of S809 airfoil at .

Figure 8 shows the pressure coefficient of S809 airfoil under incompressible flow condition which is presented in order to illustrate the variation of pressure coefficient with the compressible flow. The pressure contour shows that there is a shock on both top and bottom walls of the airfoil. As angle of attack increases shock shifts its positions as shown in Figures 9(a)9(c) and at above 8 degree the shock has a remarkable change at the lower surface (Figure 9(d)). It is observed that at compressible flow condition pressure suddenly changes both in upper faces and lower face of the airfoil and its position changes with the change of angle of attack.

Figure 8: Coefficient of pressure along airfoil at and .
Figure 9: (a) Coefficient of pressure along airfoil at , (b) coefficient of pressure along airfoil at , (c) coefficient of pressure along airfoil at , and (d) coefficient of pressure along airfoil at .

Figure 10 shows the translation of the shock at the top surface, and Figure 11 shows the translation of the shock at the bottom surface. As the distribution illustrates, the shock changes its position as angle of attack increases. For top surface the shock moves towards the leading edge and for bottom surface the shock moves towards the trailing edge. A combined plot of translation of shock at both the top and bottom surfaces is presented in Figure 12 which indicates the relative change in position of the shock as angle of attack increases.

Figure 10: Translation of shock on the top surface of airfoil.
Figure 11: Shock propagation at the bottom surface of airfoil.
Figure 12: Shock propagation comparison.

Figure 13 shows the contour of mach number variation at . As it is expected the mach number changes drastically in position of the shock and at the downstream it decreases. This happens due to the change of stagnation density and stagnation pressure. At the point of shock the sudden release of pressure energy converts it to kinetic energy and increases the mach number. Figure 14 shows the variation of mach number at the top surface of the airfoil which has the same explanation as the contour does. A relative comparison between the variation of mach number at the top and the bottom surfaces of the airfoil is shown in Figure 15.

Figure 13: Mach number contour at .
Figure 14: Mach number variation on the top surface at .
Figure 15: Mach number variation comparison on top and bottom surfaces at .

Figure 16 illustrates the turbulent viscosity contour at . As expected highly unsteady turbulent behaviors have been observed. Variation of nondimensional wall distance (Y+) that determines the turbulent behavior and the quality of mesh to resolve viscosity effect [25] has been shown in Figures 17 and 18. Figure 17 shows the variation of Y+ on the top surface of the airfoil which indicates that the Y+ value is little high to resolve the near-wall region.

Figure 16: Turbulent viscosity at .
Figure 17: Top surface Y+ variation at .
Figure 18: Average Y+ on top surface at .

Figure 18 illustrates the average Y+ between the airfoil chord length with interval which indicates that the mesh size near the first half of the airfoil surface needs to be refined to resolve better result as it is highly turbulent region. Some other fluid properties have also been studied during this investigation. Fluid density and temperature are some of them. Figure 19 shows the density contour at which clearly indicates the variation of density during compressible flow. The change of density is much intense at the stagnation point, and the stagnation density reaches approximately , while in the shock region, the density becomes much smaller from the actual air density. The pressure change due to shock causes this according to the law of compressibility.

Figure 19: Density contour.

On the other hand Figure 20 illustrates the total temperature contour at and Figure 21 shows the variation of temperature on the top surface of the airfoil. The operating temperature has been considered 300 K, but at the stagnation point, the temperature increases, and it reduces until the shock happens. At the position of the shock the temperature suddenly jumps due to sudden transform of energy. This energy comes from the pressure energy that is suddenly released due to shock.

Figure 20: Total temperature.
Figure 21: Temperature distribution at the top surface.

The change of and with respect to has also been studied. It is found that increases as increases (Figure 22). First it increases at a steady rate, but after 8° the lift coefficient increases rapidly. The same pattern has also been observed for drag coefficient as shown in Figure 23.

Figure 22: Integrated lift coefficient with respect to angle of attack ().
Figure 23: Integrated drag coefficient with respect to angle of attack ().

Flow separation has also been observed during the simulation. Observation indicates that after the shock the flow separation starts (Figure 24), and as increases flow separation also occurs more rapidly.

Figure 24: Separation of flow just after the shock at angle of attack .

5. Conclusion

Compressible flow analysis around NREL Phase VI S809 airfoil is the primary objective of this work. CFD simulations using SST k-ω turbulence model are performed. Flow mach no. 0.8 is used for the simulations. The results show the pressure distribution and effect of shock around the airfoil and the translation of shock as a function of angle of attack. The shocks are found to move towards the leading edge on the top surface and towards the trailing edge on the bottom surface. Different unsteady flow phenomena and change in fluid property have also been observed and showed. These variations are needed to be taken into account during design of large wind turbine blades where compressible flows are expected under normal wind speed conditions.


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.


  1. US Government, Department of Energy, US Department of Energy on Green House Gases, 2009,
  2. P. Jain, Wind Energy Engineering, chapter 1, McGraw Hill, New York, NY, USA, 2011.
  3. Y. L. Chang, S. L. Yang, and O. Arici, “Flow field computation of the NREL S809 airfoil using various turbulence models,” in Energy Week-96, Book VIII, Volume I: Wind Energy, pp. 172–178, ASME, 1996. View at Google Scholar
  4. K. Y. Chien, “Predictions of channel and boundary-layer flows with a low-Reynolds-number turbulence model,” AIAA Journal, vol. 20, no. 1, pp. 33–38, 1982. View at Google Scholar · View at Scopus
  5. D. C. Wilcox, Turbulence Modeling for CFD, DCW Industries, La Cañada, Calif, USA, 1994.
  6. L. Chen, C. Xu, and X. Lu, “Numerical investigation of the compressible flow past an aerofoil,” Journal of Fluid Mechanics, vol. 643, pp. 97–126, 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. H. Tijdeman and R. Seebass, “Transonic flow past oscillating airfoils,” Annual Review of Fluid Mechanics, vol. 12, pp. 181–222, 1980. View at Google Scholar
  8. B. H. K. Lee, “Self-sustained shock oscillations on airfoils at transonic speeds,” Progress in Aerospace Sciences, vol. 37, no. 2, pp. 147–196, 2001. View at Publisher · View at Google Scholar · View at Scopus
  9. M. Wang, J. B. Freund, and S. K. Lele, “Computational prediction of flow-generated sound,” Annual Review of Fluid Mechanics, vol. 38, pp. 483–512, 2006. View at Publisher · View at Google Scholar · View at Scopus
  10. J. G. Marvin, L. L. Levy, and H. L. Seegmiller, “Turbulence modelling for unsteady transonic flows,” AIAA Journal, vol. 18, no. 5, pp. 489–496, 1980. View at Google Scholar · View at Scopus
  11. C. L. Rumsey, M. D. Sanetrik, R. T. Biedron, N. D. Melson, and E. B. Parlette, “Efficiency and accuracy of time-accurate turbulent Navier-Stokes computations,” Computers and Fluids, vol. 25, no. 2, pp. 217–236, 1996. View at Publisher · View at Google Scholar · View at Scopus
  12. S. Deck, “Numerical simulation of transonic buffet over a supercritical airfoil,” AIAA Journal, vol. 43, no. 7, pp. 1556–1566, 2005. View at Google Scholar · View at Scopus
  13. A. J. Smits and J. P. Dussauge, Turbulent Shear Layers in Supersonic Flow, American Institute of Physics, College Park, Md, USA, 1996.
  14. P. R. Spalart, “Detached-eddy simulation,” Annual Review of Fluid Mechanics, vol. 41, pp. 181–202, 2009. View at Publisher · View at Google Scholar · View at Scopus
  15. X. Lu, S. Wang, H. Sung, S. Hsieh, and V. Yang, “Large-eddy simulations of turbulent swirling flows injected into a dump chamber,” Journal of Fluid Mechanics, vol. 527, pp. 171–195, 2005. View at Publisher · View at Google Scholar · View at Scopus
  16. A. Ahlstrom, Aeroelastic simulation of wind turbine dynamics [Ph.D. thesis], Royal Institute of Technology, Department of Mechanics, Stockholm, Sweden, 2005.
  17. F. R. Menter, “Zonal two equation k-ω turbulence models for aerodynamic flows,” AIAA Paper 93-2906, 1993. View at Google Scholar
  18. F. R. Menter, “Two-equation eddy-viscosity turbulence models for engineering applications,” AIAA Journal, vol. 32, no. 8, pp. 1598–1605, 1994. View at Google Scholar · View at Scopus
  19. J. L. Tangler and D. M. Somers, NREL Airfoil Families for HAWTs, Doc, AWEA, New York, NY, USA, 1995.
  20. UIUC Airfoil Coordinates Database,
  21. ANSYS ICEM CFD 12. 0 documentation, chapter 4: Meshing in FLUENT, 2011.
  22. W. Sutherland, “The viscosity of gases and molecular force,” Philosophical Magazine Series 5, vol. 36, no. 223, pp. 507–531, 1893. View at Publisher · View at Google Scholar
  23. S. V. Patankar and D. B. Spalding, “A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows,” International Journal of Heat and Mass Transfer, vol. 15, no. 10, pp. 1787–1806, 1972. View at Google Scholar · View at Scopus
  24. D. M. Somers, Design and Experimental Results for the S809 Airfoil, Airfoils, State College, Pa, USA, 1989.
  25. S. M. Salim and S. C. Cheah, “Wall Y+ strategy for dealing with wall bounded turbulent flow,” in Proceedings of the International MultiConference of Engineers and Computer Science (IMECS '09), vol. 2, March 2009.