Unsteady Double Wake Model for the Simulation of Stalled Airfoils ()
1. Introduction
Airfoil aerodynamics is one of the key subjects in the wind energy field. With the fast growing wind turbine dimensions new challenges have arisen for scientists and engineers in the sector; thick airfoils, high Reynolds numbers and sub-sonic flow conditions are some of the new obstacles that the wind energy community is facing. However, besides the new challenges, flow separation, which has been one of the main difficulties faced by aerodynamicist in the last decades, remains a hot-topic nowadays. An all-times challenge in aerodynamics has been to predict the airfoil characteristics at high angles of attack, in the stalled region, where the flow partially or completely leaves the airfoil forming a separated flow region. It is well-known that even highly demanding Navier-Stokes codes with sophisticated turbulence models fail in accurately predicting the correct behavior of airfoils in the deep stall region.
With the aim of capturing the overall physics of separated flows over airfoils using vortex theory, Maskew and Dvorak [1] developed a simplified model based in an inviscid flow solver, or panel method, which could accurately simulate steady flows around airfoils at high angles of attack. Following this idea, Marion et al. [2] developed a similar model focusing in the deep stall region. An unsteady version of the model, capable of taking into account the dynamic effects of the separated wake was first developed by Vezza and Galbraith [3] three decades ago. This approach has been recently revisited by two different research groups [4] [5] which included an integral boundary solver to compute the separation location. All three approaches modelled the separated region with a set of time-updated point vortices, while different approaches were used to represent the airfoil contour as well as the handling of the first released wake panels.
In the present work, a combination of the different approaches has been implemented, making possible to carry out the flow simulations with a small time step and a large number of vortex singularities. Moreover, the separation point is calculated externally with a viscous-inviscid strong coupling procedure, improving the accuracy of its prediction. In what follows the main characteristics of the model will be outlined. For more detailed information about the USDWM the reader is kindly referred to the work carried by Antoine [6] during his stay at DTU Wind Energy.
2. Numerical Model
2.1. Potential Flow Background
If the flow around a body with a free stream velocity is considered, at a point in the domain surrounding the solid body the velocity can be defined as the sum of the free stream velocity and the induced velocity u,
(1)
If the flow is considered to be incompressible, inviscid and irrotational, u can be expressed as the gradient of the potential, ∇Φ, with Φ satisfying the Laplace equation. The general solution to the Laplace equation can be obtained through a source and vorticity distribution around the body contour. Moreover, the shedding wake can be modelled with downstream converging point vortices. The velocity induced at any point in the domain can therefore be written as the superposition of the body and the wake singularities as follows,
(2)
2.2. Problem Formulation
In practice, as sketched in Figure 1, the airfoil contour is discretized in N panels, with each of the panels represented by a linear distribution of vorticity, γ1, …, γN+1. Where γ1 is the trailing vorticity of the lower trailing edge panel and γN+1 can be seen as the leading vorticity of the upper trailing edge panel. The vorticity released at the trailing edge and the separation points is first taken into account as two uniform panel distributions, γw SEP and γw TE, which later are convected into the separated wake region as point vortices. Finally a constant source distribution is applied over the airfoil, σ.
Hence a total number of N + 4 unknowns have to be determined. The equations necessary in order to solve the system are presented in what follows,
Equation (1) to N: Neumann non-penetration condition applied at the center of each one of the panels.
(3)
Equation N + 1: Kutta condition [7].
(4)
Equation N + 2: vortex strength at the separation panel.
(5)
Equation N + 3: Kelvin’s theorem implementation.
Figure 1. Sketch of the singularity distributions in the USDWM model.
(6)
Equation N + 4: zero vortex strength at N + 1
(7)
The system of equations is closed and the unknowns can be determined. It is important to note that the influence of the point vortices in the wake is included in the right hand side of the system, and their strength remains constant in time until they are merged by coalescence or enter the far wake model. Once the system of equations has been solved, the length and angle of the trailing edge and the separation panels is updated as follows,
(8)
(9)
(10)
(11)
Convergence is obtained once the residual value of the above parameters lTE, lSEP, θTE, θSEP as well the total airfoil circulation, ΓB, is lower than a threshold, in this case set to 10−4.
Finally, the surface pressure coefficient, Cp is calculated using the unsteady Bernoulli equation and accounting for the increase in total pressure over that at infinity in the separated region,
(12)
where Δh is the total pressure jump at the interface of the two wake sheets, which only adopts a non-zero value inside the wake region. The head jump can be expressed as the difference between the heads on both sides of the wake sheets,
(13)
3. Simulations
To first validate the implemented model, fixed angle of attack simulations have been carried out. In what follows we will analyze both, the average aerodynamic coefficients as well as the instantaneous values driven by the vortex shedding. In all the simulations presented here 80 panels have been used to represent the airfoil surface, the time step is set to 0.004 s, the free-stream velocity is 1 m/s and the chord length is fixed to 1m. Moreover, the maximum number of point vortices in the wake has been fixed to 6000 and a coalescence criterion has been used in order to reduce the number of point vortices in the far wake. The coalescence radius between particles is fixed to 0.04 and all vortices located more than four chords downstream from the trailing edge are candidates for the merging. The flow past the airfoil has been simulated a total time of 50 seconds, therefore a total of 25,000 vortices have been released in the wake.
The FFA-W3-211 airfoil has been chosen for the study. The 21% thick airfoil was tested in the low speed wind tunnel L2000 [8], in the Royal Institute of Technology, Stockholm. Measurements were performed at a Reynolds number of 1.8 × 106 and a turbulence intensity of 0.15%. Simulations with the USDWM have been carried out for flows past the FFA airfoil at four different angles of attack: 12, 15, 19 and 23 degrees. The separation location, which is a necessary input to the USDWM is computed using the VI solver Q3UIC [9]. Q3UIC solves the integral form of the boundary layer equations using a strong viscous inviscid coupling procedure. All boundary layer calculations in Q3UIC have been carried out with natural transition at a Reynolds number of 1.8 × 106 and a turbulence intensity of 0.15%.
The instantaneous flow visualization at the final simulated time of 50 seconds is shown in the left side of Figures 2-5. The lift coefficient variation with time is presented on the right side of the figures. The length scales of the vortex shedding are directly related with the angle of attack, α, and therefore with the size of the separated region. It is seen how the deviations of the instantaneous Cl from its average value are closely related to α as well.
Figure 6 shows the calculated lift and drag coefficients in comparison with Q3UIC and the experimental data. Generally a good agreement is obtained in terms of average values of Cl and Cd. Low angles of attack computations present a small standard deviation and vice versa, clearly relating the size of the shedding structures to the instantaneous variations in Cl. There is a tendency to slightly over predict the average Cl at high angles of attack.
The energy spectra of the vertical fluctuating velocity component calculated at a point in the airfoils wake is shown in Figure 7 and briefly analyzed for two different angles of attack. The point is situated at the airfoils
Figure 2. USDWM simulation of the flow past the FFA-W3-211 airfoil at an angle of attack of 12 deg and a Reynolds number of 1.8 × 106. (left) Instantaneous flow visualization at t = 50 s. (right) Instantaneous and average value of Cl.
Figure 3. USDWM simulation of the flow past the FFA-W3-211 airfoil at an angle of attack of 15 deg and a Reynolds number of 1.8 × 106. (left) Instantaneous flow visualization at t = 50 s. (right) Instantaneous and average value of Cl.
Figure 4. USDWM simulation of the flow past the FFA-W3-211 airfoil at an angle of attack of 19 deg and a Reynolds number of 1.8 × 106. (left) Instantaneous flow visualization at t = 50 s. (right) Instantaneous and average value of Cl.
Figure 5. USDWM simulation of the flow past the FFA-W3-211 airfoil at an angle of attack of 23 deg and a Reynolds number of 1.8 × 106. (left) Instantaneous flow visualization at t = 50 s. (right) Instantaneous and average value of Cl.
Figure 6. (left) Lift force coefficients and (right) drag force coefficients for flows past the FFA-W3-211 airfoil at a Reynolds number of 1.8 × 106.
Figure 7. Energy spectra of the vertical fluctuating velocity at a location 2 chords dowstream the leading edge of the FFA-W3-211 airfoil and at its c/4 height.
quarter-chord height and a downstream distance of two chords from the leading edge. The USDWM is capable of capturing the variation in the shedding frequency with the change in angle of attack. Simulations at high angles of attack, where the airfoil presents large separated flow regions, exhibit low shedding frequencies. In massive separated cases, the distance between the two shedding layers is larger, making the interaction slower. At lower angles of attack, with small separated bubbles, the interaction between both shedding layers is much faster, and therefore the shedding frequency is larger.
4. Conclusion
The implemented Unsteady Double Wake Model is capable of capturing accurately the average aerodynamic coefficients of a typical wind turbine airfoil under stall conditions. Moreover the model is able to predict the change in the shedding frequency as well as in the size of the shed vortex structures as a function of the angle of attack. Further work has to be done to validate the accuracy of the predicted unsteady aerodynamic quantities as well as the vortex shedding characteristics.
Acknowledgements
The authors would like to acknowledge the support from the Danish Council for Strategic Research for the project Off Wind China (Sagsnr. 0603-00506B) and the European Union’s Seventh Programme for research, te- chnological development and demonstration for the project “AVATAR: AdVanced Aerodynamic Tools for large Rotors” (FP7-ENERGY-2013-1/no. 608396).