**UFRJ**

# Aerofólio S809 - 1 - s2 0 - s0306261911002595 - numerico

(Parte **1** de 2)

Numerical modeling of an S809 airfoil under dynamic stall, erosion and high reduced frequencies

Kobra Gharali, David A. Johnson⇑ Wind Energy Group, Department of Mechanical and Mechatronics Engineering, University of Waterloo, Waterloo, Canada N2L 3G1 article i nfo

Article history: Received 31 October 2010 Received in revised form 1 April 2011 Accepted 12 April 2011 Available online 1 July 2011

Keywords: Wind turbine Dynamic stall S809 Airfoil Numerical simulation Reduced frequency Erosion abstra ct

An oscillating freestream over a stationary S809 airfoil is simulated numerically using ANSYS Fluent 12.1. For comparison the laminar-turbulent transition is simulated with the realizable k– and SST k–x models, and ﬁnally the SST k–x model is chosen. Where possible, the resultant aerodynamic coefﬁcients presenting dynamic stall phenomena are compared with aerodynamic coefﬁcients from existing experimental and semi-empirical data for a pitch oscillating S809 airfoil. The simulation was extended for several Reynolds numbers, 104,1 05 and 106, for a large range of the reduced frequency, from 0.026 to 18. The level of the dependency of the results on the reduced frequency is examined based on wake velocity proﬁles, vorticity contours and aerodynamic coefﬁcients. Leading edge erosion, which is a common problem for in-service wind turbine blades, is simulated based on the range of erosion found in turbine operating conditions. For the eroded airfoil, the lift decrease, causing a severe impact on wind turbine performance, was investigated in detail. 2011 Elsevier Ltd. All rights reserved.

1. Introduction

Sometimes, natural wind has a strong wind shear with altered direction and velocity magnitude that causes dynamic stall (DS) phenomena in horizontal axis wind turbines (HAWT) which may result in dynamic blade loading, variable perfomance and failing stall regulation parameters. Since DS phenomena cannot be prevented, studying the unsteady aerodynamics is crucial to understanding the effect and assists in modifying common wind turbine designs. Experimental methods, semi-empirical strategies such as Beddoes–Leishman (BL) models, and numerical models have been applied to predict the aerodynamic loads and ﬂow conditions of wind turbines during DS phenomena. Examining both the steady and unsteady behavior of airfoils, wind tunnel studies were one of the ﬁrst methods used to deﬁne the ﬂow phenomena. An airfoil that has received much attention in investigating DS phenomena is the S809 airfoil. Ramsay et al. [1] tested a two dimensional (2D) S809 airfoil under stationary and dynamic conditions. In the dynamic case, the airfoil underwent pitch oscillations at different mean angles of attack and oscillating amplitude for where q and l are the air density and viscosity, U1 is the freestream velocity and C is the airfoil chord length. The effects of air- foil roughness were also considered in their study. Somers [2],a t the Delft low speed, low-turbulence wind tunnel, measured aerodynamic coefﬁcients for different Reynolds numbers. The results of these reported measurements have been used as two references for developing models, conducting new experimental techniques, and advancing computational ﬂuid dynamic (CFD) codes. The Beddoes–Leishman (BL) model, a very popular semiempirical DS model, has been used for modeling dynamic stall phenomena for HAWT although the BL model was originally developed for helicopters. So far, this model has been developed and used for predicting the unsteady aerodynamic characteristics of S809 airfoils by researchers such as Sheng et al. [3,4], Gupta and Leishman [5], Gonzalez and Munduate [6] and Johansen [7]. According to Sheng et al. [3], no DS model has been considered universal yet.

As well as dynamic stall there are other features of airfoils at high reduced frequencies and low Reynolds numbers that have been investigated. These features include thrust and lift generation due to oscillation. Gursul et al. [8] experimentally investigated the effects of the reduced frequency on the lift in an unsteady freestream. Drag reduction on an oscillating airfoil or an airfoil in an oscillating freestream is known as the Katzmayr effect [9]. If the two rows of vortices in the wake of an airfoil, a Kármán vortex pattern, are reversed, a thrust force on an airfoil is produced [10].A t higher reduced frequencies, the wake vortex pattern, a reverse Kármán street, generates thrust on the airfoil. Dependency of thrust generation on different parameters such as reduced frequency was investigated in detail by Sarkar and Venkatraman [1]. Although that work was based on the ﬂow ﬁeld over a harmonically pitching airfoil, the role of these parameters should also be important for an oscillating freestream over a stationary airfoil.

0306-2619/$ - see front matter 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.apenergy.2011.04.037

⇑ Corresponding author. Tel.: +1 519 8 4567. E-mail address: da3johns@uwaterloo.ca (D.A. Johnson).

Applied Energy 93 (2012) 45–52 Contents lists available at ScienceDirect

Applied Energy journal homepage: w.elsevi er.com/locat e/apenergy

Compared to restricted experimental methods based on test speed ranges, measuring equipment, and accuracy and size of the model, well-established CFD methods are progressing to investigate complex aerodynamic problems attributed to unsteady, transient and dynamic ﬂow. Due to developments in commercial computational software it is feasible to compute unsteady airfoil ﬂows.

This study is organized to give an overview of the numerical simulation using a commercial CFD package, ANSYS Fluent 12.1, as an accurate, time efﬁcient and economic way of simulating an oscillating freestream over a stationary S809 airfoil for different Re and a range of reduced frequencies. Results for dynamic stall are compared with prior experimental and semi-empirical pitch oscillating S809 airfoil studies. Finally, the behavior of the vorticity ﬁelds, wake velocity proﬁles and aerodynamic coefﬁcients are provided in detail to demonstrate the connection of the results to Re and reduced frequency. The study also seeks to understand the effect of blade erosion on the unsteady performance of the airfoil, where the erosion was simulated based on real cases from in-service blades.

2. Dynamic stall

Varying natural wind conditions result in altered directions and magnitude of velocities near turbine blades. These inhomogeneous velocity proﬁles cause unsteady operating conditions for wind turbine rotors such as unsteady angle of attack changes. A combination of stall phenomena, inherently unsteady, and various types of unsteady angle of attack (a) motion create dynamic stall (DS) phenomena. DS results in the development of a dynamic stall vortex and delay of stall, meaning a higher DS angle compared to a static stall angle. Because of these two signiﬁcant features of DS, the dynamic lift coefﬁcient exceeds the static lift coefﬁcient. DS inﬂuences the loading of turbine blades that use stall regulation to control maximum loads. Thus, analysis based on static stall conditions may not be valid.

3. Erosion

Based on the application of the S809 airfoil in wind turbines, it was designed to be roughness insensitive [2], but for unsteady conditions, the effects of the roughness are considerable compared with those of steady conditions and reduce the aerodynamic performance of the S809 [2]. Ramsay et al. [1] considered contamination effects by applying leading edge grit roughness and Sheng et al. [4] used two sand strips on the upper and lower surfaces of the airfoil to evaluate the effects of roughness on a S809 airfoil. In this study, the effects of the leading-edge local erosion on performance of the S809 airfoil for dynamic cases were considered. Blade erosion is a common occurrence with wind turbine airfoils installed in service. Erosion is caused by impacts from airborne dirt and debris and insects on primarily the leading edge surfaces of the airfoil. The impacts wear the surface and create rough and pock marked leading edge surfaces. A closeup of a turbine blade leading edge near the blade tip as removed from a large in-service wind turbine is shown in Fig. 1. The maximum extent of the measured erosion is about 23% of the maximum airfoil thickness covering the majority of the leading edge area, not shown in Fig. 1. For this blade the leading edge is the only region that shows any indication of erosion and will be the focus of the modiﬁcations in this numerical study. Here, three models of a leading edge eroded blade were developed; Erosion A, B and C Fig. 2. To examine a range of erosion characteristics, erosions were simulated with two thicknesses 12% areas commonly have sharp edges which have been considered in the modeling.

4. Computational method

In this study, the sinusoidal oscillating freestream over a stationary airfoil has a reduced frequency k,

where X =2 pf and f is the oscillation frequency. The oscillation system is a ¼ amean þ aampsinðXtÞð 2Þ where amean is the mean angle of the wind with relative velocity and aamp is the pitch oscillation amplitude. For this simulation, the airfoil is aligned with the horizontal axis.

Fig. 1. Photograph of the leading edge of a wind turbine blade removed from service. Leading edge dimensions shown are 4 m horizontally and 20 m vertically. Largest scale divisions are 5 m.

Fig. 2. S809 Airfoil with model erosion A, B and C.

46 K. Gharali, D.A. Johnson/Applied Energy 93 (2012) 45–52

4.1. Boundary conditions

Here, a code was added as a user-deﬁned function to the commercial CFD code, Fluent [12]. The code alters the direction of the far-ﬁeld ﬂow over the stationary airfoil at each time step based on the sinusoidal equation to simulate a proper a or wind direction for the boundary conditions as shown in Fig. 3.

4.2. Grid generation

Using a C grid layout, different cell sizes for a quadratic mesh have been used to make the aerodynamic forces independent of the grid size Fig. 3. The mesh contained 60,0 cells. Around the airfoil 390 nodes were distributed with high resolution on the leading and trailing edges. In the wake of the airfoil 170 nodes and 80 nodes are placed in a horizontal and in vertical lines, respectively. In addition, 10,0 nodes are located in the eroded area for those simulations. The thickness of the cells around the airfoil is adapted for Y+ 1, Y+ is the non-dimensional wall distance for a turbulent boundary layer deﬁned as Yþ ¼ qu yl and uf is the friction velocity. The eroded area, modeling erosion, was intro- duced as connected blocks to the main block with triangular meshes. The boundary of the computational domain was located approximately 20C from the airfoil, a distance that is recommended for the pressure far ﬁeld boundary condition.

4.3. Viscous model

For the range of Re in this study, modeling of the laminar-turbulent transition was necessary. To simulate laminar-turbulent transition, one technique utilized a ﬁxed transition point on the lower and upper sides of the airfoil resulting in two zones to simulate transition effects. These two laminar and turbulent zones were meshed separately with ﬁne and coarse grids for the turbulent and laminar zones, respectively. The forward ﬂow zone was set to the laminar zone manually. Then, the rest of the domain, the turbulent zone, was solved based on the turbulent model. One of the drawbacks of using a ﬁxed transition point model was that the transition locations have to be available for each a and Re. Since the angle of attack changes in dynamic conditions, the location of the transition point also changes. Therefore, this method has been used just for the steady case.

Advanced near-wall treatment models, written in the C language, can be integrated with the code as user-deﬁned wall functions. In Fluent, the k– model can have a user-deﬁned wall function for the near-wall treatment. It should be noted that for the ﬂow conditions in this study, the SST k–x and realizable k– models are two of the suggested models [12]. Therefore, in this study, for the SST k–x model, the Low-Re correction option was activated and for the realizable k– model, a user-deﬁned wall function has been used [12].

4.4. Simulation setup

The 2D domain geometry simulation used the segregated solver and implicit and absolute velocity formulation to set up the numerical simulation in Fluent. The Simple algorithm was chosen for coupling the momentum pressure equations. For spatial discretization, a second-order upwind differencing scheme was applied. Since the time step size is a crucial parameter for unsteady cases and is a function of amplitude, frequency and the far ﬁeld velocity, a number of time-step reﬁnements have been employed to ensure the temporal accuracy of the results and to get the details of the ﬂow [13]. The time step is proportional to the characteristic s = 0.001 provided the best results. The selected time steps were small enough to make the results independent of the time step size. As the time step decreases, signs of unsteadiness will appear especially for separated ﬂows. On the other hand, the selected time steps are large enough to damp out the unsteadiness.

5. Results and discussion

All the results are provided based on the instantaneous data from ANSYS Fluent 12.1. It should be noted that the lift coefﬁcient is deﬁned as CL ¼ LqU C where L is the lift force. The drag coefﬁcient is similarly deﬁned as CD ¼ DqU C where D is the drag force.

5.1. Numerical approach for turbulence modeling

5.1.1. Static case

Fig. 4 shows the lift and drag coefﬁcients for Re =1 06 in the static case for the realizable k– and SST k–x methods and published experimental as well as numerical data. Two ﬁxed transition points for the lower and upper sides of the airfoil were considered for each a. The transition locations were provided in previous studies by Drela and Giles [14]. The simulation results are compared with experimental data reported by Somers [2] and numerical data for fully turbulent ﬂow by Johansen [7]. The experimental data shows a change in slope of the coefﬁcient of lift at approximately 6 until approximately 9 where the lift coefﬁcient begins to decrease. At 13 the lift begins to increase again. Compared with the experimental data, the SST k–x has the least discrepancy as it is able to capture the variation in the lift coefﬁcient from 9 to 13 . The numerical results match the experimental curve well for small angles of attack, but for high angles of attack, the SST k–x method was preferred. The realizable k– method agrees with the fully turbulent simulation for the drag coefﬁcient, but for the lift coefﬁcient, at very high angles of attack, the results of this method follow the trend of the experimental results rather than that of the fully turbulent results. The numerical results from ﬁxed

Fig. 3. Sinusoidal freestream and C-type mesh around S809 airfoil.

(Parte **1** de 2)