Methodology for the engineering calculation of flaps on Wind Turbines using BEM codes

Aeroelastic codes based on Blade Element Momentum theory are the standard used by many wind turbine designers. These codes usually include models and corrections for unsteady aerodynamics, tip and root effect, tower shadow and other effects. In general, this kind of codes does not include models to adequately simulate aerodynamic control devices. This paper presents a method to take into account the unsteady contributions due to the flap motion (based on indicial models) and the spanwise 5 effects (based on circulation theory), in order to simulate flaps on the blades. This method can be included in BEM codes in general and it could also be applied to another kind of control devices. The validation and verification show the accuracy of this method using experimental data for two-dimensional unsteady cases, and CFD for three-dimensional steady and unsteady cases.


Introduction
The classical implementation of the blade element momentum theory is based on steady flow (equilibrium state of wake) and annular ring independence.The proper simulation of the flaps requires an unsteady model for the angle of attack, the flap motion and to include the annular ring dependence due to the spanwise effect of the flaps.The aim of this work is to develop models for the BEM codes in order to correctly simulate the effects of aerodynamic control devices, and the study is focused on trailing edge flaps.The code used for the present study is FASTv8 (Ref.Jonkman et al. (2005) , Website), an aeroelastic code developed by NREL and based on the BEM theory.

Approach and methods
Most of aeroelastic codes consider the unsteady effect of the variation of the angle of attack, and some of them also consider the effect of the flap in the 2D calculation from a steady point of view.These codes can obtain the aerodynamic response with a direct interpolation of static polars through the flap angle.However, this calculation is not complete without the unsteady Significant research efforts related to the unsteady contributions of the flaps have been made by Bernhammer (Ref. Bernhammer (2015)), Baek (Ref. Baek et al. (2011)), Bergami (Ref. Bergami (2008)) and Gaunaa (Ref. Gaunaa (2006)).Similar modifications in order to simulate dynamic flaps have been adapted to the formulation of FASTv8 and explained in this paper.
The spanwise effects observed in wind turbines have been studied by Pirrung, and a near wake model has been developed in order to reduce the limitations of the BEM codes (Ref.Pirrung et al. (2016a), Pirrung et al. (2016b)), using a coupling between the models (not specifically for flaps).This work is based on the work of Beddoes (Ref. Beddoes (1987)), and the developments of Wang and Coton (Ref. Wang and Coton (1999)) and Madsen (Ref.Madsen and Rasmussen (2004)).There is an specific application for flaps developed by Andersen (Ref.Andersen et al. (2010)).In the present work, the study is focused on the spanwise effect of the flaps.This method is basically based on BEM (with the introduction of vorticity concepts) and leads to lower computational cost.

Unsteady method
The use of smart blade technologies, for example trailing edge flaps, is one of the options to control the loads on wind turbine blades.One of the main objectives is the reduction of fatigue loads and, consequently, this technology must adjust lift quickly and attenuate these loads.It is important to notice that the simulation of the unsteady contribution will be a critical issue for the validity of the flap modelling.Therefore, it is necessary to take into account the current state and also the speed of the changes.Currently, engineering codes usually include models for unsteady aerodynamics: Beddoes-Leishman (Ref.Leishman and Beddoes (July 1989)), ONERA (Ref.Tran et al. (1981)) or Øye (Ref. Øye (1991)).
The oscillations of a trailing edge flap are equivalent to a change in the angle of attack (based on thin airfoil theory, Ref. Leishman (1991), Bisplinghoff et al. (1955)).Therefore, it is possible to use a similar formulation in the calculation of the aerodynamic load coefficients, based on Beddoes-Leishman model (Ref.Leishman and Hariharan (Sept-Oct 1996)).
Considering the flap deflection (β), a first order delay is applied to calculate the effective value (β ef f n ).This delay takes into account the unsteady contribution with a formulation similar to the one used for the angle of attack for the n time step.X βn and Y βn are deficiency functions depending on the previous value of the flap deflection (β n−1 ).The effective value, β ef f n , is then used to get the force coefficients.
The parameters used in the unsteady formulation (i.e.C N α , slope of the 2D normal force coefficient curve and α 0 , zerolift angle of attack) are obtained through a direct interpolation of the input data of the different deflections, using the β ef f n previously calculated.For instance, the normal coefficient obtained with FAST will have a similar formulation to the calculation without flap, using different values of the parameters.
In the tangential coefficient, an additional term is added, in order to include the unsteady effect of the flap deflection.
The main novelty of this implementation with respect to previous works is the unsteady contribution of C M due to the variation of β.It depends not only on the variation of C L but also on the variation of the point of application of the resultant C L .In the case of flaps, the hysteresis observed in the C M cycles is smoother than the hysteresis of the C L cycles.As this variation is not known using the current methodology, the calculation of C M is interpolated directly with another β ef f n (β ef f CM n ), which is defined with a semi-empirical delay.
Currently, the unsteady contribution of the flaps has been included in Aerodyn 14 and in Aerodyn 15 , combined with an improved version of the Beddoes-Leishman model called DYSTOOL (Ref.González et al. (2014)), developed by CENER, which is available in the code FASTv8.The two-dimensional validation have been performed using Aerodyn 15, and the three-dimensional verification with Aerodyn 14.

Three-dimensional method
Pure BEM formulations do not take into account spanwise three-dimensional effects of the flow.However, this type of effects is included using corrections, for example root corrections.A significant difference is observed between the results for a finite flap in a three-dimensional simulation and the two-dimensional approach, therefore the modelling of a flap should include its effect along the blade.Some of the definitions used in other aerodynamic methods are very useful for this new methodology.
First, the definition of circulation and the Kutta-Joukowski theorem is considered.The circulation (Γ) for steady conditions is obtained with the lift coefficient calculated by the BEM code.
The circulation is a macroscopic measure of the rotation for a finite area of the fluid.The circulation could be presented as horseshoes elements (according to the Helmholtz Vortex Theorems) on every strip of the BEM (Figure 1), and the strength of the vorticity emitted from the extremes of each strip is equal to the circulation of that strip.Studying separately these extremes, it is clear that the total vorticity presented at these points is the difference of circulation between the adjacent strips (Figure 2).
The current methodology has been only applied to the extremes of the flap (Figure 3).At the flap extremes, there is a break of the circulation, with a step of the lift coefficient and the emitted vorticity cannot be neglected.
The strength of the vortex filaments is the difference of the circulations in the adjacent regions to the extreme of the flap, whose value is significant.The vortex filaments generate an induced velocity, according to the Biot-Savart Law.Then, the emitted vortex line generates a three-dimensional effect affecting the induced velocity and, consequently, the angle of attack and the force distribution in the closest regions.The vorticity obtained at the extremes of the flaps is simplified to a straight semi-infinite vortex filament, whose strength is obtained with the following formulation.Where γ is the strength of the vortex filament and ∆Γ is the difference of circulation between the adjacent sections and the studied extreme.The parameters c and V are the chord and the relative wind speed of the profile at that section; and ∆C L the step of the lift coefficient between the adjacent sections to the extreme of the flap.
The induced velocity is calculated with the Biot-Savart law.The induced velocity at a section is presented as v i , and r represents the distance between that section and the extreme of the flap.
In the core region (r → 0), the effect of the vorticity is attenuated in order to avoid unrealistic induced velocities when the distance approaches zero.The increment of angle of attack (∆α) due to the induced velocity is approximately: Where ω is the rotational speed and R is the radius of the studied section.
The final equation implemented in the aerodynamic code is the following simplification.
The filament vortex has been simplified to a straight vortex line with constant strength in order to avoid the additional computational cost of the wake induction calculation.However, the straight vortex can be modified in order to simulate more approximately the helical shape of the wake.
This methodology has also been applied to unsteady cases, neglecting the unsteady effect of the circulation.For an oscillating flap case (k ≈ 0.15), this simplification involves deviations close to 6% in terms of ∆Γ max and a delay of approximately 0.3 rad.

Preliminary description of the results
The unsteady two-dimensional validation was performed using experimental data of a flap oscillation (Ref.Leishman and Hariharan (Sept-Oct 1996)).The first case is the oscillation of a 25% chord flap in a NACA64A006 profile, at Reynolds Number ≈ 2 M .The amplitude of the flap movement is ±2.5 o , with a reduced frequency k = 0.098 and the angle of attack is 0 o .In order to obtain the two-dimensional result using the current aeroelastic code, the simulation for this case considers        The comparison of the lift coefficient with respect to the angle of attack is included in figures 12 and 13 for both cases.The shape of the first comparison is quite similar to the comparison with respect to the flap deflection.There are again two different regions and the analysis is similar to the previous comparison of the case.The position of the experimental turning point is very close to the turning point obtained in the simulation, which demonstrates that the differences observed in the previous comparison of the case (with respect to the flap deflection) is due to the delay observed between the flap motion and the pitch oscillation.The second case presents a reasonable good agreement in the clockwise loop.The most significant difference in this case is observed for the maximum lift coefficient (previously mentioned in the comparison with respect to the flap oscillation).
After the validation of the unsteady method implemented, a verification of the three-dimensional method for the flap integration is presented.The accuracy of this new methodology is shown in figures 14 and 15 for a rotating case using the INNWIND.EU rotor without tilt, precone and tower shadow (ω = 8.97 rpm, pitch = 4.97 o ) and a deformable static 10%c flap, centred at 75%R, width 14%R and β = 10 o (where R is the distance between the tip and the axis of the wind turbine).
These figures show the comparison with respect to CFD results (Ref.Gomez-Iradi (2009)).The results obtained using FAST without the spanwise effects of the flaps have also been included.The increment of the C N due to the effect of the flap is very similar for both codes along the blade, while the C T values present higher differences, related to differences of the static polars.
However, the C T trend of the CFD results is clearly captured in the sections close to the extremes of the flap.The second case for the three-dimensional validation has been performed with a comparison between different aerodynamic codes, studied within the AVATAR Project (Ref.González et al. (2016)).This case simulates the AVATAR rotor without tilt,    Finally, an unsteady three-dimensional case was performed with the modified aeroelastic code (Figures 18 and 19) within the European Project AVATAR (Ref.Ferreira et al. (2015)).This case simulates the AVATAR rotor without tilt, precone and tower shadow (ω = 9.6 rpm, pitch = 2.06 o ) with an oscillating 10%c flap (f β = 0.96 Hz, centred at 75%R, width 10%R and It is compared to the results obtained with a vortex method and with two different CFD codes.Figure 18 shows the axial force distribution along the blade for the maximum and the minimum flap deployments and presents similar results of the force along the blade compared to the rest of the codes.In general, the force obtained with FAST is slightly lower than for the rest of the codes.This deviation is related to a difference of the result for the case without flap that could be due to small differences in the static polars.Figure 19 presents the comparison at the mid flap section, and it can be observed a similar slope with respect to the vortex result.The maximum and minimum values present the same deviation previously described.

Conclusions
This paper demonstrates the possibility to properly simulate flaps in wind turbine blades with a BEM code.The method developed has been presented, including the theories related to this development and the validation (in this case, the aeroelastic code used is FASTv8, Ref. Jonkman et al. (2005), Website).The validation and verification show the accuracy of this method using experimental data for two-dimensional unsteady cases, and for three-dimensional steady and unsteady cases.The results of the validation are very encouraging and present similar trends to CFD results and experimental data.
It is also important to highlight that one of the main advantages presented by the BEM codes with respect to CFD or vortex methods is the computational cost required for the simulations, which is maintained in the same order with the current modifications.
The implementation of the spanwise 3D effect detailed in this work can also be applied to the analysis of blades with abrupt changes of the geometry (i.e.transition between different profiles, important step of the twist or chord and others).
The limitations of this methodology are going to be analyse in future validations.The suitability of this method (including the unsteady method) for another type of aerodynamic control device are going to be studied in future works (i.e.vortex generators, leading edge flaps).

Figure 1 .
Figure 1.Two adjacent strips of the blade, shown as horseshoes defined by their circulations
straight, parked blades, rigid components and effects (tip and root, tower shadow) are disabled.The agreement between the lift and moment coefficients with respect to the experimental data can be observed in figures 4 and 5.Both cycles present similar behaviour to the experimental data, in terms of maximum and minimum values; and the width of the loops.Another validation case has been studied, using experimental data carried out by CENER (in the DTU -Lyngby Campus, Red Wind Tunnel).It is an oscillating 15% chord flap on a NACA64418 profile, with a reduced frequency k = 0.100 at Reynolds Number 0.5 M .The amplitude of the flap oscillation is ±5 o , and the angle of attack of this case is 0 o .The simulation of this case presents the same considerations explained in the previous validation case (straight, parked blades, rigid components and effects are disabled).The agreement of the lift and moment coefficients is shown in the figures 6 and 7, where the Cl cycle is similar to the experiments, while the Cm cycle shows similar loop with a slight deviation towards lower values.Some of the cases of the previously mentioned Wind Tunnel Campaing (carried out by CENER in the DTU -Lyngby Campus, Red Wind Tunnel) include combined motions of the flap and the oscillation of the pitch of the airfoil.The studied airfoil is Wind Energ.Sci.Discuss., doi:10.5194/wes-2016-50,2016 Manuscript under review for journal Wind Energ.Sci.Published: 22 December 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 14 .
Figure 14.∆CN for the INNWIND.EU rotor and β=10 o with respect to the blade without flap.

Figure 15 .
Figure 15.∆CT for the INNWIND.EU rotor and β=10 o with respect to the blade without flap.

Figure 16 .
Figure 16.∆CN for the AVATAR rotor and β=10 o with respect to the blade without flap.

Figure 17 .
Figure 17.∆CT for the AVATAR rotor and β=10 o with respect to the blade without flap.

Figure 18 .
Figure18.CN for the AVATAR rotor, unsteady case.Figure19.CN at 75%R for the AVATAR rotor, unsteady case.

Figure 19 .
Figure18.CN for the AVATAR rotor, unsteady case.Figure19.CN at 75%R for the AVATAR rotor, unsteady case.