A Novel Technique for the Calculation of Post-Stall Aerodynamic Coefficients for S809 Airfoil

In the autumn of 2000, the National Renewable Energy Laboratory (NREL) organised a 'blind comparison', in which wind turbine research centres worldwide were asked to provide stationary aerodynamic loads for the wind turbine Phase VI and aerodynamic coefficients for the S809 airfoil. These institutes participated in this benchmark with their wind tunnel measurements and numerical calculation programs. The results furnished by the participants showed a large discrepancy and scatter, in particular in deep stall. Based on the differences between the provided aerodynamic properties, the need for improvements on the current design codes and for implementation of new numerical method was imposed. The purpose of the present work is to reproduce the aerodynamic coefficients for S809 airfoil in attached flow and in-stall development, by using a coupled potential method and the boundary layers theory. For the fully stalled region, a novel technique was introduced, based on the AERODAS model developed by NASA. This model has been improved to consider the effect of Reynolds number, and to prevent numerical instability, which might be produced by the transition from pre-stall to post-stall regions. The results of our calculations were approximately in the mid-range of those from the 'blind comparison' participants while they were not highly disparate to the measurements.


INTRODUCTION
In the spring of 2000, NREL carried out a series of measurements for the wind turbine Phase VI, in the world's biggest wind tunnel at NASA-Ames. This wind turbine has two blades with S809 airfoils, 10.058 m rotor diameter, with twisted and tapered blades, 3 deg tip pitch angle, developing a rated power of 19.8 kW, and was mainly operated at 72 rpm. Based on the high quality of these measurements, a 'blind comparison' was arranged by NREL, in which wind turbine research organisations worldwide (such as ECN in Netherlands, and Risø in Denmark) were called to provide the aerodynamic loads for some stationary conditions, without knowledge of the measured results [1]. The research institutes participated with different design codes based on the various aerodynamic models, and with some of their wind tunnel measurements. The results furnished by the participants in this benchmark showed huge differences, particularly in stall. These differences were found to depend partly on the choice of aerodynamic coefficients and on the stall delay models used for the consideration of rotational effects. So, many investigations were conducted, based on the measured data for the S809 airfoil in Delft [2], Ohio State [3], and Colorado State [4] Universities, to improve the calculation codes for stationary and rotational aerodynamic coefficients. From these wind tunnel measurements, it has been concluded that the tunnel-wall effect leads to a deviation of the slope of the lift and drag curves.
Somers calculated the aerodynamic coefficients for S809 airfoil by using the design code Eppler [5]. Bernadett and Van Dam used the MSES code [6] to calculate the aerodynamic coefficients for a Reynolds number of 10 6 . The results of MSES code demonstrated a good agreement with the measurements in the TU-Delft wind tunnel, except for the larger angles of attack when an under-estimation of the drag coefficients is obtained. Zhang et al. [7] and Wolfe and Ochs [8] calculated the aerodynamic coefficients of the S809 airfoil using the CFD method (which requires a large amount of computational power, extensive time and storage space) using different turbulence models respectively for a Reynolds number of 0.65·10 6 and 2·10 6 . The aerodynamic coefficients are precise through approximately 9° angle of attack. Above this angle, the calculations do not pick up the airfoil's stall behaviour and, therefore, over-predicts the lift and drag coefficients. Lanzafame and Messina [9] used a fifth-order logarithmic polynomial to fit experimental data for attached flow and stall development regimes, and mathematical functions from the flat plate theory to implement the post-stall aerodynamic coefficients for S809 airfoil in their simulation code. Spera [10] from Jacobs Technology, Inc., under the sponsorship of the NASA Glenn, applied an empirical approach in which the trends of measurements data are modelled by a set of mathematical equations. The developed model, called AERO-DAS, allows the calculation of the lift and drag coefficients in unstalled and stalled aerodynamic regimes. The greatest criticalities for the flat-plate theory and the AERODAS model are represented by the non-consideration of the Reynolds number effect, and the need for some input parameters in the equations describing these models, which should be derivable from a minimum of test data, due to the reduced time available for wind tunnel test and increasing test costs [11].
Lindenburg from ECN also formulated an empirical tool, named StC, for the deep-stall aerodynamic coefficients at different angles of attack [12], comparison of the deep-stall aerodynamic coefficients of StC with measured coefficients still revealed differences, but with a slight improvement in the normal-force and tangential-force coefficients.
From the follow-up analysis, it was concluded that determining the aerodynamic coefficients for an airfoil is still an important challenge, especially in deep stall.
In previous works [13,14], authors showed how to mathematically reproduce the aerodynamic coefficients for the S809 airfoil in the pre-stall region using an interacted viscous-inviscid approach. For this approach, the major assumption suggests that the flow can be separated in two sectors: within the boundary layer where the fluid is treated as viscous and outer this layer where the flow is treated as potential [15]. This approach is computationally efficient and yields reliable solutions for attached flow, but it cannot be exploited for approximate poststall calculations [16]. The inviscid and viscous flows are respectively computed employing the potential theory of wing sections and the boundary layers theory. The resulting velocity distribution from the potential flow calculation is used as input in the boundary layer equations to find a so-called blowing velocity distribution. To include the viscous effects, this blowing velocity is fed back to the potential flow formulation to determine the viscous velocity distribution around the airfoil. This operation is replicated until the convergence of the calculation.
In the present paper, particular attention was gave to the modelling of the aerodynamic coefficients for S809 airfoil in the fully stalled region using an improved AERODAS model. This model is modified to consider the Reynolds number effect, and to prevent the numerical instability, which might be caused by the transition from pre-stall to post-stall regions.

MATHEMATICAL MODEL
The experimental lift coefficient results for the typical angle of attack range (from -20° to 90°) and for distinct Reynolds numbers, show that there are two different regions " Figure.  For the pre-stall region, the interacted viscous-inviscid approach is used for the calculation of the lift and drag coefficients for the S809 airfoil. The objective is to calculate the air flow around the S809 airfoil using the potential flow formulation, which is based on the classical assumption of irrotational flow, coupled with the boundary layers method. For the potential flow formulation, the potential theory of wing sections is implemented in previous works by the authors [13,14] and deals with the determination of the velocity distribution around airfoil, where V is the undisturbed velocity at infinity, α the angle of attack, ε T the zero-lift angle, ψ represents a real function of the angle θ (θ varies from zero to 2π), and Ψ 0 is a constant quantity. The equation (1) yields U e = ∞ near the leading edge. To determine the velocity near the airfoil nose, we consider the leading edge as analogous locally to a circle; consequently, the tangential velocity on the surface can be written as [17], To consider the effect of the angle-of-attack variation, we have implemented the coefficient cosα, we get then, To determine the velocity near the tail, we can write θ = π + Δθ and ε = ε T + ε'Δθ and, where Δθ is a small angle.
From the equation (1), we get the velocity potential at the tail, The Prandtl-Glauert rule [18] was used for compressibility correction.

Boundary layers formulation
The viscous effect and the contribution of the skinfriction component in the calculation of the total drag coefficient are predominant. So, the boundary layers modelling is required. For this purpose, authors used the momentum integral equation [14], is the skin-friction coefficient, τ w is the wall shear stress and V s expressed the velocity over the wall as there is the possibility of flow due to suction passing through a porous wall, in the habitual case of an impermeable wall V s = 0. This momentum integral equation must be solved to yield the variations of the boundary-layer parameters, such as the boundary layer thickness δ, the displacement thickness δ * , the momentum thickness θ * and most importantly the skin-friction coefficient C f along the surface of the airfoil.
To calculate the characteristics of the boundary layer in the laminar region, up to the onset of transition, it is more stable and easier to use the Thwaites' method [19]. Thwaites has suggested the use of an empirical correlation to provide closure for the momentum integral equation in laminar boundary layer, With b = 0.45 and a = -6. The equation (6) can be solved analytically, we get, The shape factor H can be evaluated by the Cebeci and Bradshow semi-empirical formula [20] To investigate the turbulent boundary layer, it is crucial to predict the location of the laminar-to-turbulent transition, and to determine the starting values of θ * and H. For this purpose, Cebeci and Smith relation based on Michel's criterion [21], is used. According to this criterion, the transition takes place when, The Head's method [22], which is a typical integral method and especially fast, is used for solving the turbulent boundary layer.
This method uses an entrainment velocity parameter E, which is given by, The shape factor H 1 and the momentum thickness θ * , in turbulent boundary layer, can be determined by simultaneous solution of the following equations, Separation in the boundary layer occurs when the tangential flow velocity changes sign and recirculation occurs, a point of inflexion develops in the velocity profile and the shear stress at the wall vanishes, Reverse flow occurs and the boundary layer separates from the airfoil surface at some point.
The boundary layer equations are only valid up to the streamwise location of flow separation. The turbulent boundary layer separation is predicted when H > 2.4 or H 1 > 3.3 [23], or when the value of the skin-friction coefficient drops below C f = 0. After the separation point, the model assumes that the pressure remains constant on the profile. This pressure is assumed to be equal to the pressure at infinity in the undisturbed flow.
The momentum thickness and the shape factor will not be developed in the separation regime, and are supposedly diminished, contrary to the boundary layer thickness. The drag coefficient can be calculated by the Young Square formula, where the boundary layer characteristics are considered at the trailing edge [

Viscous/Inviscid coupling
Separate formulations for the external potential flow and the viscous boundary layer are mixed interactively to afford a composite solution to the entire flow field ( Figure.2). For this purpose, the strong interaction method [14] which employs the blowing velocity concept as an interacting agent is adopted. This blowing velocity distribution is used as a boundary condition in the inviscid method in order to include the effects of boundary layer in the potential flow solutions. This interactive procedure is reiterated until convergence is achieved. Figure 3 summarises, under schematic form, the structure of the numerical flow model.

Post-stall behavior
AERODAS is a mathematical model developed to meet the requirements for equations to calculate the aerodynamic coefficients of a blade airfoil, operating in unstalled and stalled aerodynamic regions. As mentioned above, the AERODAS model [10] is adopted to predict the airfoil lift and drag coefficients in the fully stalled region.
The lift coefficient in the post-stall region is given by the following equations, The drag coefficient in the post-stall region is given by, where, RCL = 1.632 -CL max , N = 1 + CL max /RCL, ACL: angle of attack at maximum pre-stall lift, CL max : maximum post-stall lift coefficient, A: angle of attack at which pre-stall lift coefficient worth 0, ACD: angle of attack at maximum pre-stall drag, CD 1max : maximum pre-stall drag coefficient, CD 2max : maximum post-stall drag coefficient.
The AERODAS model input parameters for the S809 airfoil, based on TU-Delft wind tunnel measurements are listed in table 1.
The greatest criticality of the AERODAS model is represented by the non consideration of Reynolds number. A single curve is fitted for all wind speeds.
To overcome this shortfall, the authors have introduced a trick that consists the consideration of the first point of the post-stalled AERODAS curve, which starts from the last point of the pre-stalled interacted viscousinviscid curve, which corresponds to CL atrans and CD atrans for the lift and drag coefficients respectively ( Figure. 4). Therefore, all the points forming the poststalled AERODAS curve are shifting in the same way as the first point of this curve and then a specific curve is obtained for each number of Reynolds. For this, the following equations are implemented, Also, this trick has helped to avoid the numerical instability which might be caused by the transition from pre-stall to post-stall regions.

RESULTS AND DISCUSSIONS
In this section, some illustrative results are presented to study the capability of the developed computer code (named AEROC) to evaluate the two dimensional S809 aerodynamic coefficients and the characteristics of the boundary layer, such as, displacement thickness, momentum thickness, the location of the separation point from the airfoil surface, and the inception point of the transition phenomenon.
In Figure 5, the location of the onset of transition phenomenon for the upper surface of the S809 airfoil at Re = 2·10 6 is sketched, and the beginning and end of this phenomenon are demonstrated. The comparison with experiments shows good agreement.
In Figure 6, the location of the inception of the separation phenomenon is presented. The boundary layer equations used in AEROC can determine the position of the separation phenomenon. In this figure, it is clear that the separation surface at the upper surface of the airfoil increases with the increasing angle of attack. In other words, the separation phenomenon relocates toward the leading edge. In Figures 7 and 8, the displacement thicknesses of S809 airfoil at the angle of attack of 0° and for two different Reynolds numbers are illustrated. Based on these plots, we can conclude that the increase in the Reynolds number leads to the reduction of the boundary layer thickness. Figure 9 illustrates the momentum thickness over the S809 airfoil at angle of attack of 5°. The position of the transition point on the upper surface of the airfoil has been captured by a red ellipse. With each alteration in the Reynolds number, based on the momentum thickness, the location of the transition point changes.
Indeed, the increase in momentum thickness leads to the displacement of the transition point forward to the leading edge of the airfoil.     On the other hand, at the lower surface of the airfoil, the location of the transition point depends only on edge velocity distribution and is not affected by the momentum thickness.
In Figures 10 and 15, plots of lift and drag coefficients for S809 airfoil are respectively sketched for different Reynolds numbers. The accord between AEROC code and experimental results can clearly be seen. Figures 10 and 11 show a comparison between experimental and calculated lift coefficient for the respective Reynolds number of 0.65·10 6 and 1·10 6 . A good accord over the pre-stall region is obtained for the two Reynolds numbers.
In addition, the lift coefficients provided by AEROC concede well with the experimental results at the beginning of the post-stall region, until 35° angle of attack for 0.65·10 6 and 26° for Re = 1·10 6 .
It must be noted that the post-stall modelling using the modified AERODAS model is highly dependent on the latest pre-stall 2D aerodynamic coefficients. Tables 2 and 3 compare the lift coefficients for the same Reynolds numbers. In the pre-stall region, the lift coefficients predicted by AEROC (resp. AERODAS) are between 1% and 18% (resp.4% and 39%) and between 1% and 25% (resp.6% and 27%) far the experiment results for the respective Reynolds numbers of Re = 0.65·10 6 and Re = 1·10 6 .
For the starting part of the post-stall region, the relative errors given by AEROC (resp. AERODAS) are within 10% (resp.35%) and 18% (resp.18%) for respectively Re = 0.65·10 6 and Re = 1·10 6 . The CFD model of Zhang et al. gives relative errors for Re = 0.65·10 6 within 52% and 47% for respectively the pre-stall and post-stall regions. Figures 12 and 13 show the predicted drag coefficients curve of AEROC code, the predicted drag coefficients curve of AERODAS model, and the experimental drag coefficients curve for the respective Re = 0.65·10 6 and 1·10 6 . Notice how the predicted drag coefficients curve, given by AEROC code, is closer to the experimental results, which is better than AERODAS model for large ranges of the angles of attack belonging to both pre-stall and post-stall regions. Therefore, the improved AERODAS model has given rise to the maximisation of the correlation between predicted and experimental data in the post-stall regime. Tables 4 and 5 contain a comparison between experimental and numerical drag coefficients for Re = 0.65·10 6 and Re = 1·10 6 .
Analysis of the effect of the Reynolds number on these coefficients of the airfoil is presented in Figures 14  and 15. As evidenced in these graphs, an increase in Reynolds number involves an increase in lift coefficient while it's would bring about a decrease in the drag coefficient.

CONCLUSION
This paper focuses on the characterisation of the S809 airfoil with accepted accuracy and low computational costs. A numerical code, based on a strong viscousinviscid interaction approach, has been developed for the pre-stall region. The flow is analysed using the potential theory of wing sections, coupled with the boundary layers theory to take into account the viscous effect. In the fully stalled region, the AERODAS model is adopted to calculate the lift and drag coefficients. This model is improved to consider the Reynolds number effect, and to avoid the numerical instability, which might be caused by the transition from pre-stall to post-stall regions.
The numerical results have been compared with wind tunnel measurements found in the literature and in the work reported by the Netherlands institute ECN under number 74189.
Generally, a good agreement is obtained against experiments for the predicted aerodynamic lift and drag coefficients for large ranges of attack angles and different Reynolds numbers. The approach exposed in this paper represents an interesting contribution that responds to a global aerodynamic challenge, and that can be used to provide a good compromise between accuracy and computational costs, compared to the CFD methods and the empirical models generally based on wind tunnel tests, such as the Flat plate and AERODAS models.