Abstract
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.
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.
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.
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.
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.
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.
(a)
(b)
(c)
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.
(a)
(b)
(c)
(d)
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 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 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 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.
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.
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.
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.
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.
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.