Five Megawatt Wind Turbine Power Output Improvements by Passive Flow Control Devices

The effects of two types of flow control devices, vortex generators (VGs) and Gurney flaps (GFs), on the power output performance of a multi-megawatt horizontal axis wind turbine is presented. To that end, an improved blade element momentum (BEM)-based solver has been developed and BEM-based computations have been carried out on the National Renewable Energy Laboratory (NREL) 5 MW baseline wind turbine. The results obtained from the clean wind turbine are compared with the ones obtained from the wind turbine equipped with the flow control devices. A significant increase in the average wind turbine power output has been found for all of the flow control device configurations and for the wind speed realizations studied in the present work. Furthermore, a best configuration case is proposed which has the largest increase of the average power output. In that case, increments on the average power output of 10.4% and 3.5% have been found at two different wind speed realizations. The thrust force and bending moment in the root of the blade have also been determined and compared with the values of the clean wind turbine. A residual increase in the bending moment of less than 1% has been found.


Introduction
The rise of installed wind power in Europe for the last fifteen years, along with the increasing significance of offshore wind energy, shows the relevance of research in the field of flow control for large wind turbines.The considerable increase of wind turbine rotor size and weight in recent years has made it impossible to control as they were controlled 20 years ago.Rotors of 120 meters, or even more, are now a reality.Johnson et al. [1] compiled some of the most important flow control techniques that could be used in wind turbines to assure a safe and most favourable operation under different atmospheric conditions.To maximize the lifetime energy captured by the wind turbine is a key factor in the wind energy field.
Many different flow control devices have been developed in the last decades.Most of them were created for aeronautical issues and this was its first research field and application [2].They are also frequently used in turbo machinery [3].Nowadays researchers are working to optimize and introduce these types of devices in multi-megawatt wind turbines.Wood [4] developed a four layer scheme which allows classifying the different concepts that are part of all flow control devices.In the study of Shires and Kourkoulis [5] a tangential air jet was used on a vertical axis wind turbine (VAWT) blade to control the separation of the flow and, therefore, to increase the aerodynamic performance.Additionally, the dynamic stall control was investigated by Xu et al. [6] on a S809 airfoil (NREL's S-Series airfoils, Golden, CO, USA) by the implementation of a co-flow jet.
Depending on their operating principle control techniques can be classified as active or passive [7].Passive control techniques represent an improvement in the turbine's efficiency and in load reduction without external energy consumption.Active control techniques need an additional energy source to obtain the desired effect on the flow and, unlike vortex generators (VGs) and other passive devices, active flow control needs complex algorithms to obtain the maximum benefit (see Becker et al. [8]).Johnson et al. [1] conducted an analysis and discussed 15 different devices for wind turbine control.Some of them are still being tested on full-scale turbines.
A VG is defined as a passive flow control system which modifies the boundary layer (BL) fluid motion bringing momentum from the outer flow part into the inner flow part of the wall-bounded flows.Its main goal is to delay the separation of the flow and increase the maximum lift coefficient C L,max .VGs are intended in order to re-energize the BL by transferring momentum between the free stream velocity and the near wall region.VGs have been investigated for more than fifty years for a wide range of applications in aerodynamics and airplane wings [9].They are small vanes, usually triangular or rectangular (Figure 1), inclined at an angle to the oncoming flow and placed as close as possible of the leading edge (LE).They are generally assembled spanwise on the suction side of the blade and present the advantage that they can be added as a post-production fix to blades that do not perform as expected.Its height is usually similar to the boundary layer thickness at the VG position.
Energies 2017, 10, 742 2 of 15 (VAWT) blade to control the separation of the flow and, therefore, to increase the aerodynamic performance.Additionally, the dynamic stall control was investigated by Xu et al. [6] on a S809 airfoil (NREL's S-Series airfoils, Golden, CO, USA) by the implementation of a co-flow jet.Depending on their operating principle control techniques can be classified as active or passive [7].Passive control techniques represent an improvement in the turbine's efficiency and in load reduction without external energy consumption.Active control techniques need an additional energy source to obtain the desired effect on the flow and, unlike vortex generators (VGs) and other passive devices, active flow control needs complex algorithms to obtain the maximum benefit (see Becker et al. [8]).Johnson et al. [1] conducted an analysis and discussed 15 different devices for wind turbine control.Some of them are still being tested on full-scale turbines.
A VG is defined as a passive flow control system which modifies the boundary layer (BL) fluid motion bringing momentum from the outer flow part into the inner flow part of the wall-bounded flows.Its main goal is to delay the separation of the flow and increase the maximum lift coefficient CL,max.VGs are intended in order to re-energize the BL by transferring momentum between the free stream velocity and the near wall region.VGs have been investigated for more than fifty years for a wide range of applications in aerodynamics and airplane wings [9].They are small vanes, usually triangular or rectangular (Figure 1), inclined at an angle to the oncoming flow and placed as close as possible of the leading edge (LE).They are generally assembled spanwise on the suction side of the blade and present the advantage that they can be added as a post-production fix to blades that do not perform as expected.Its height is usually similar to the boundary layer thickness at the VG position.Fernandez-Gamiz et al. [10] and Urkiola et al. [11] studied the behaviour of a rectangular VG on a flat plate and the streamwise vortices produced to investigate how the physics of the wake behind VGs in a negligible streamwise pressure gradient flow can be reproduced in CFD simulations.In the work carried by Gao et al. [12] on a 30% thick DU97-W-300 airfoil (Delft University family airfoils, Delft, The Netherlands), the maximum lift coefficient was increased from 1.5 to approximately 2 due to the implementation of passive VGs.When the angle of attack increases, both the lift and drag coefficients rise to values higher than the ones reached in the steady state conditions.A vortex structure grows and the airfoil goes into a stall situation, then the loss of lift is larger than in steady operation [13].The non-linear and unsteady aerodynamic behaviour of a multi-megawatt horizontal axis wind turbine (HAWT) is a problem from structural and electric generation points of view.The use of VGs can led to a reduction of periodic loads, thus improving power output and cyclic fatigue life, as described in Gebhardt et al. [14].Øye [15] compared the measured power curves with VGs and without them on a 1 MW wind turbine.Although quite rough methods were used for the VG design optimization, the experiment showed that, for a stall-regulated wind turbine, power increased nearly 24% by using VGs through field tests.Furthermore, Sullivan [16] conducted an experiment on a 2.5 MW wind turbine to test the effects of adding VGs on the power conversion performance.An increase of 11% in the annual energy production was found.
VGrs are usually mounted in a spanwise array on the suction side of the blade and have the advantage that they can be added as a post-production fix to blades that do not perform as expected.On the other hand, their main disadvantage is that drag increases (CD) due to the implantation of this device, an undesirable feature for this kind of application.Great care is also needed to be taken in Fernandez-Gamiz et al. [10] and Urkiola et al. [11] studied the behaviour of a rectangular VG on a flat plate and the streamwise vortices produced to investigate how the physics of the wake behind VGs in a negligible streamwise pressure gradient flow can be reproduced in CFD simulations.In the work carried by Gao et al. [12] on a 30% thick DU97-W-300 airfoil (Delft University family airfoils, Delft, The Netherlands), the maximum lift coefficient was increased from 1.5 to approximately 2 due to the implementation of passive VGs.When the angle of attack increases, both the lift and drag coefficients rise to values higher than the ones reached in the steady state conditions.A vortex structure grows and the airfoil goes into a stall situation, then the loss of lift is larger than in steady operation [13].The non-linear and unsteady aerodynamic behaviour of a multi-megawatt horizontal axis wind turbine (HAWT) is a problem from structural and electric generation points of view.The use of VGs can led to a reduction of periodic loads, thus improving power output and cyclic fatigue life, as described in Gebhardt et al. [14].Øye [15] compared the measured power curves with VGs and without them on a 1 MW wind turbine.Although quite rough methods were used for the VG design optimization, the experiment showed that, for a stall-regulated wind turbine, power increased nearly 24% by using VGs through field tests.Furthermore, Sullivan [16] conducted an experiment on a 2.5 MW wind turbine to test the effects of adding VGs on the power conversion performance.An increase of 11% in the annual energy production was found.
VGrs are usually mounted in a spanwise array on the suction side of the blade and have the advantage that they can be added as a post-production fix to blades that do not perform as expected.On the other hand, their main disadvantage is that drag increases (C D ) due to the implantation of this device, an undesirable feature for this kind of application.Great care is also needed to be taken in their integration with the blade so as not to deteriorate the performance of the wind turbine or the aeroelastic conditions.An in-depth review of the use of VGs for of boundary-layer flow-separation control is presented in [17].
The Gurney flap (GF) is a simple vane with a height between 1% and 2% of the airfoil chord length c, located perpendicular to the lower or upper side of the airfoil near the trailing edge.According to Liebeck [18], GFs increase the total lift of the airfoil while reducing the drag once appropriately sized.Storms [19] found a lift increase of around 13% for a GF size of 0.5% c with minimal or no drag penalties for low and moderate lift coefficient values.The studies of van Dam et al. [20,21], based on the research of Meyer et al. [22], studied the effects of serrated and split GFs to avoid the vortex shedding.
In order to study the impact of VGs and GFs on wind turbines, to optimize their position and distribution, CFD tools can be used such as in the work by Troldborg et al. [23] where computational simulations of a wind turbine rotor with and without a row of VGs on the blades were carried out.In those simulations, a 10 MW wind turbine rotor and airfoils were fully resolved by a grid, while the VGs were modeled using a modified version of the BAY model developed by Bender et al. [24].According to Troldborg et al. [23], it is not computationally feasible to directly simulate a full wind turbine blade equipped with many VGs due to the large range of scales present in a rotor equipped with VGs.Therefore, modelling a fully-meshed rotor with VGs and GFs becomes prohibitively expensive because of its small size compared with the rotor dimensions.

Blade Element Momentum Method
The blade element momentum (BEM) method is the most frequently-used tool to calculate the wind turbine rotor's aerodynamic performance since it is computationally economical and, consequently, very fast [25].In addition, the BEM method gives reasonable results provided that usually high-quality airfoil data are accessible for the lift, drag, and moment coefficients, and also at different Reynolds numbers.The method assumes that all sections of the rotor are independent of each other and, consequently, can be treated individually.The momentum loss is produced by the axial loads of the flow passing the blades, causing a pressure drop over the blade section.Since the induced velocity generated by the action of the loads is known, the local angle of attack at a given radial sector on a blade can be calculated.As described in Hansen [26], two corrections to the BEM method are required to obtain good-quality results.The Prandtl's tip loss factor correction is the first one, which solves the hypothesis of an infinite number of blades.The second correction is applicable when the axial induction factor achieves values greater than about 0.2-0.4,where the simple momentum theory fails.This is known as Glauert [27] correction and consists of an empirical relation between the thrust coefficient C t and the axial induction factor a.
The BEM-based algorithm was developed and programmed by the authors of the current study based on the numerical iterative approach of Hansen [26].All of the necessary equations were derived and computed based on the steps proposed by the classical BEM method.The usual basic steps for BEM calculations were as follows: Energies 2017, 10, 742 4 of 15 (7) State a tolerance for a and a and if it has changed more than that tolerance, go to (b), or else finish.(8) Compute the local loads.
The power coefficient C p calculation has been implemented as follows: For each r i station in the blade the working point is evaluated in order to calculate the torque and the power captured by this station.The input conditions are given by the wind velocity and the rotor speed.Once the working point defined by a i and a i is well known for each r i station, the torque contribution of this station can be calculated with the following equation: When the algorithm must calculate the working point (a i , a i ), usually an iterative equation system is applied.The stop condition is usually given by a maximum number of steps or enough low values of change in the working point (a i , a i ).The presented algorithm is a Gauss Seidel-like equation solver and consists of the following equations: The variables used are described in Table 1.In the present implementation, the solver has followed the flow charge described in Figure 2 for each iteration.Following the proposed scheme of Figure 2, for each value at iteration k, the previous equation system from Equations ( 2) to (10) gives a new value for iteration k + 1.It is simple to propose a grid of values of from 0 to 2π, and calculate the function ( + 1) defined by Equation (11).When this function crosses the line ( + 1) = ( ), it is possible to say that the solver has obtained the value of Φ that solves the equation system from Equations ( 2) to ( 10): In this way, for a given value of Φ, it is possible to obtain the values of and .Once these two parameters have been estimated with Equations ( 2) to (10), the torque contribution of this station is calculated.These calculations can be repeated for all of the stations in the blade and, finally, the total torque can be calculated by the numerical integration of these torque contributions for a given rotor speed, wind speed, and given pitch values.Once the total torque has been estimated for a given state, (wind speed, rotor speed, and pitch angle), the power captured by the wind turbine can be directly estimated by multiplying the rotor speed by the total torque.This algorithm allows the calculation of the working point of each station without any numerical oscillations or instabilities because all possible values of ( ) are calculated.Hence, there is no possibility of numerical convergence problems.

NREL 5 MW Reference Wind Turbine
The NREL 5 MW reference wind turbine is being widely used in research studies of the wind energy field since it represents a baseline of the modern and future offshore HAWT [28].Many investigations have been carried out based on this wind turbine concept, including studies about rotor aerodynamics, controls, offshore dynamics, and design code development.This conception of a 5 MW wind turbine is based on the data from the DOWEC study [29,30]; with a concept from the UpWind project [31].The airfoils and chord schedule used in the present work are the same as [28], also adopted from the DOWEC project.More exhaustive information about the DU Delft University family of airfoils used can be found in Timmer et al. [32].The reported NREL 5 MW airfoil distribution is shown in Table 2.
Table 2. Airfoil distribution along the blade of the NREL 5 MW wind turbine described in [28].Following the proposed scheme of Figure 2, for each Φ i value at iteration k, the previous equation system from Equations ( 2) to (10) gives a new Φ i value for iteration k + 1.It is simple to propose a grid of values of Φ from 0 to 2π, and calculate the function Φ(k + 1) defined by Equation (11).When this function crosses the line Φ(k + 1) = Φ(k), it is possible to say that the solver has obtained the value of Φ that solves the equation system from Equations ( 2) to (10):

Station
In this way, for a given value of Φ, it is possible to obtain the values of a i and a i .Once these two parameters have been estimated with Equations ( 2) to (10), the torque contribution of this station is calculated.These calculations can be repeated for all of the stations in the blade and, finally, the total torque can be calculated by the numerical integration of these torque contributions for a given rotor speed, wind speed, and given pitch values.Once the total torque has been estimated for a given state, (wind speed, rotor speed, and pitch angle), the power captured by the wind turbine can be directly estimated by multiplying the rotor speed by the total torque.This algorithm allows the calculation of the working point of each station without any numerical oscillations or instabilities because all possible values of Φ(k) are calculated.Hence, there is no possibility of numerical convergence problems.

NREL 5 MW Reference Wind Turbine
The NREL 5 MW reference wind turbine is being widely used in research studies of the wind energy field since it represents a baseline of the modern and future offshore HAWT [28].Many investigations have been carried out based on this wind turbine concept, including studies about rotor aerodynamics, controls, offshore dynamics, and design code development.This conception of a 5 MW wind turbine is based on the data from the DOWEC study [29,30]; with a concept from the UpWind project [31].The airfoils and chord schedule used in the present work are the same as [28], also adopted from the DOWEC project.More exhaustive information about the DU Delft University family of airfoils used can be found in Timmer et al. [32].The reported NREL 5 MW airfoil distribution is shown in Table 2.
Table 2. Airfoil distribution along the blade of the NREL 5 MW wind turbine described in [28].

Flow Control Devices and Experimental Data Description
All of the experimental data have been taken from the AVATAR (Advance Aerodynamic Tools for Large Rotors) European project.The results from wind tunnel tests at the Low-Speed Tunnel of TU Delft at a Reynolds number of Re = 2 × 10 6 are used in the current work.For a detailed description of the experimental techniques and analysis, refer to Timmer et al. [32].The data used in the present study refers to the following airfoils: DU97W300, DU91W(2)250, and DU93W210.The configuration of the flow control devices is described in Table 3 with the inclusion of a triangular, vane-type VG array on the upper airfoil surfaces data combining with GFs.The suffixes VG20 and VG30 represent that VGs are mounted on the airfoil suction side at a distance from the LE of 20% and 30% of c, respectively.The GF2 suffix corresponds to a GF with height of 2% of c implemented on the airfoil.The notation for each case is denoted in the third column of the Table 3.

Wind Speed Model
The wind speed realizations utilized in the BEM computations, as shown in Figure 3, have been calculated with the TurbSim tool (v1.06.00,NREL, Golden, CO, USA) [33].The wind speed series have been generated with the following parameters, and the turbulence model is the normal turbulence model (NTM) following the IEC 61400 norm [34]: The velocity spectra of the IECKAI model are assumed to be invariant across the grid.In practice, a small amount of variation in the u-component standard deviation occurs due to the spatial coherence model.Figure 3a,b represent the wind speed series used in the present study according to the NTM with a 5 m/s of average velocity and with a 10 m/s average velocity, respectively.These two different wind speed realizations were chosen to investigate the effects of the VGs and GFs, since it is a good way to evaluate the wind turbine power output at low and medium wind speeds.The wind speed dataset is sorted using the "method of bins".This method divides the wind speed range into 0.5 m/s contiguous bins centered on integer multiples of 0.5 m/s. Figure 4 shows the corresponding histograms of the percentage of occurrences per bin for both cases at 5 and 10 m/s averaged wind speed.

Methodology
The primary tools used in the current work to investigate the effects of the passive flow control devices on the NREL 5 MW baseline wind turbine are engineering models.The present procedure is described in the following steps: (1) First of all, BEM-based computations were carried out in order to characterize the dynamic behavior of the NREL 5 MW wind turbine, including BEM calculation improvements described in Section 2. TurbSim uses an adapted version of Veers [35] to generate the time series based on the spectral representation.The IECKAI (IEC Kaimal) model is defined in IEC 61400-1 2nd ed.[34] and 3rd ed.[36] and assumes neutral atmospheric stability.The spectra for the wind field are given by Equation ( 12): where f and L K are the cyclic frequency of turbulent wind field and the integral scale parameter, respectively, corresponding to a value of f = 50.000Hz and L K = 340.2m.The mean wind velocity at the hub height is indicated by u.
The velocity spectra of the IECKAI model are assumed to be invariant across the grid.In practice, a small amount of variation in the u-component standard deviation occurs due to the spatial coherence model.Figure 3a,b represent the wind speed series used in the present study according to the NTM with a 5 m/s of average velocity and with a 10 m/s average velocity, respectively.These two different wind speed realizations were chosen to investigate the effects of the VGs and GFs, since it is a good way to evaluate the wind turbine power output at low and medium wind speeds.
The wind speed dataset is sorted using the "method of bins".This method divides the wind speed range into 0.5 m/s contiguous bins centered on integer multiples of 0.5 m/s. Figure 4 shows the corresponding histograms of the percentage of occurrences per bin for both cases at 5 and 10 m/s averaged wind speed.
Energies 2017, 10, 742 7 of 15 where f and LK are the cyclic frequency of turbulent wind field and the integral scale parameter, respectively, corresponding to a value of f = 50.000Hz and LK = 340.2m.The mean wind velocity at the hub height is indicated by .
The velocity spectra of the IECKAI model are assumed to be invariant across the grid.In practice, a small amount of variation in the u-component standard deviation occurs due to the spatial coherence model.Figure 3a,b represent the wind speed series used in the present study according to the NTM with a 5 m/s of average velocity and with a 10 m/s average velocity, respectively.These two different wind speed realizations were chosen to investigate the effects of the VGs and GFs, since it is a good way to evaluate the wind turbine power output at low and medium wind speeds.The wind speed dataset is sorted using the "method of bins".This method divides the wind speed range into 0.5 m/s contiguous bins centered on integer multiples of 0.5 m/s. Figure 4 shows the corresponding histograms of the percentage of occurrences per bin for both cases at 5 and 10 m/s averaged wind speed.

Methodology
The primary tools used in the current work to investigate the effects of the passive flow control devices on the NREL 5 MW baseline wind turbine are engineering models.The present procedure is described in the following steps: (1) First of all, BEM-based computations were carried out in order to characterize the dynamic behavior of the NREL 5 MW wind turbine, including BEM calculation improvements described in Section 2.

Methodology
The primary tools used in the current work to investigate the effects of the passive flow control devices on the NREL 5 MW baseline wind turbine are engineering models.The present procedure is described in the following steps: (1) First of all, BEM-based computations were carried out in order to characterize the dynamic behavior of the NREL 5 MW wind turbine, including BEM calculation improvements described in Section 2. (2) Following the specifications of the utility scale multi megawatt wind turbine NREL offshore 5 MW baseline described in [28], all of the wind turbine rotor properties were introduced as input characteristics.The polar curves of airfoils with no flow control device were taken from [32] and [29].The data of the airfoil with the flow control devices are described in Section 3.1.(3) The surfaces of the power coefficient C p were calculated for all cases of the present study according to the matrix distribution described in Table 4. (4) Once the C p surfaces have been generated, BEM-based computations are run for all of the cases and the power curve vs. wind speed is calculated to compare the power curve of the clean turbine with the curve of the cases with the flow control devices implemented.(5) Afterwards, the wind speed realizations explained in Section 4 are introduced to calculate the average wind turbine power output for all cases.(6) The results of the average wind turbine power output for all cases and at two different wind speed realizations are compared with the mean power output of the clean wind turbine, the one without any flow control devices implemented.
Table 4 illustrates a matrix with the distribution of the cases.The clean wind turbine was taken as the baseline case, without any devices implemented.The cases are different depending on the blade span position where the passive devices were implemented, the airfoil type, and the device location from the LE of the airfoil.The suffix st means the blade station where the passive devices were introduced.According to the airfoil distribution described in Table 2, the stations 7-11 were selected with their corresponding airfoil types.For simplicity, a test ID number will be used in the legends of the plots and tables.As an example, the case ID1 represents the NREL 5 MW wind turbine with VGs implemented in the upper surface of the airfoil DU97W300 at 20% of c from the LE and corresponding to blade station 7.

Results and Discussion
In order to investigate the influence of the flow control devices on the power of the NREL 5 MW reference wind turbine, BEM-based computations have been carried out following the steps explained in the previous section.The BEM computations have been derived and the power has been computed.Table 5 shows the results of the average wind turbine power calculations for the clean wind turbine and for the cases with flow control devices.Equation (13) shows how this average power is calculated: V 0,j : is the instantaneous wind speed according to the realizations shown in Figure 3.
Nbins: number of bins per data.P(V 0,j ): power at the wind speed V 0,j .
N(V 0,j ): number of data at the wind speed V 0,j .
where the P(V 0,j ) has been determined by Equation ( 14): Firstly, the average power was calculated for the clean wind turbine at the two wind speed realizations presented in Figure 3, without any flow control device mounted on the blade.The results are presented in Table 5.At NTM5 wind speed realization the average power output was 4.85570 × 10 5 W and at NTM10 3.417917 × 10 6 W. The same procedure was carried out for the wind turbine cases with the flow control devices following the test cases described in Table 4.In all wind turbine cases an increment in the wind turbine power output is achieved, but it is more notable at low wind speed.The greatest mean power value defining the best combination case ID25) at the NTM5 wind speed realization was achieved by case ID6, with a value of 5.34667 × 10 5 W, which supposes an increase of 10.111% in comparison with the value obtained by the clean wind turbine.At NTM10 wind speed realization the largest average power value (before defining the best combination case ID25) is again reached by case ID6 with an increase with respect to the clean case of 3.388%.The other cases with the flow control devices mounted at different stations experience a similar increase.The symbol ∆ in Table 5 represents the increments in the average power output of all cases, from ID1 to ID25, in comparison with the clean wind turbine at NTM5 and NTM10 wind speed realizations.The case indicated in the last line of Table 5 has been determined as a combination of the best cases from ID1 to ID24 for each airfoil type and blade station.The largest value in the average power for the cases with the flow control devices implemented into the airfoil DU97 is achieved by the case ID2; DU97W300VG30st7, consequently, was selected to be included into the best case ID25.Following the same criterion of the maximum mean power output, ID6 and ID24 cases have been chosen for the best combination case.Therefore, the wind turbine of ID25 was defined by the cases ID2, ID6, and ID24, and the average power was calculated and compared with the average power values of the clean wind turbine.At the wind speed realization of NTM5, the ID25 case achieves an increase of 10.401% and, at NMT10, the increase is 3.541% in comparison with the clean wind turbine.Both increments are the largest ones in comparison with the other cases from ID1-ID24.Those results are in concordance with the previous study by Sullivan [16], where an increase of approximately 11% in the power production was achieved.

Best Combination Case Comparison vs. Clean Case
Figure 5a illustrates the power curves vs. the wind speed for the case with no passive devices implemented into the blade in comparison with case ID25.The curve of the wind turbine with the best combination case ID25 follows the trend of the curve of the clean wind turbine.However, at the wind speeds before the rated power is achieved, the power output increases slightly in the ID25 case, as shown in the enlargement view embedded in Figure 5a.The power coefficient C p has also been computed for both the clean wind turbine and the ID25 case and is represented in Figure 5b against the wind speed.The improvement of the C p for the ID25 case is observable in comparison with the clean case.In Figure 5a,b, the calculations of the wind turbine power output and the power coefficients provided by [28] have also been added.Those parameters have been determined by FAST which is an NREL computer-aided engineering (CAE) tool for horizontal axis wind turbines widely used in wind turbine aero-elasticity calculations.The BEM-based computations of the current study follow the trends of the calculations made by [28] for both power output and C p variables.
The wind turbine power output production has also been investigated along with the duration of the wind speed realizations described in Section 4. The upper part of Figure 6a,b illustrates the comparison between the clean case and the previously-denoted best case ID25 of the power distribution with time for both wind profiles NTM5 and NTM10.In the bottom of Figure 6a,b an enlargement view of the power profiles has been inserted to show the increments of the ID25 case power with respect to the clean case.At the NTM5 wind realization, the power produced by the ID25 case follows the profile formed by the clean case with time but with a positive gap.However, at the NTM10 wind profile the power of the ID25 case is constant and equal to the wind turbine rated power of 5 MW in contrast to the power profile achieved by the clean case, which is irregular and lower than the rated power.
The gain in the power output production due to the passive flow control devices implementated in the ID25 case is clearly visible for both wind speed realizations.After applying the BEM algorithm to all control volumes, the tangential and normal load distributions are known and global parameters, such as thrust and bending moment at the root of the blade, can be computed.The mean values of both the thrust and bending moment have been calculated for each wind speed realization.The thrust T has been calculated by Equation ( 15) taking into account the thrust distribution along the blade and the bending moment M is determined by Equation ( 16): The Prandtl's tip loss correction factor F has been determined by Equation ( 17), for both thrust and bending moment calculations: The variables used for thrust and bending moment estimations, and the corresponding dimensions, are shown in Table 1.All of the calculations are based on the 5 MW reference wind turbine described in [28], where R is the rotor radius.Figures 7 and 8 illustrate the comparison between the clean case and the ID25 case for the thrust and bending moment results at the wind speed realizations of NTM5 and NTM10, respectively.The thrust force is represented by the green   After applying the BEM algorithm to all control volumes, the tangential and normal load distributions are known and global parameters, such as thrust and bending moment at the root of the blade, can be computed.The mean values of both the thrust and bending moment have been calculated for each wind speed realization.The thrust T has been calculated by Equation ( 15) taking into account the thrust distribution along the blade and the bending moment M is determined by Equation ( 16): The Prandtl's tip loss correction factor F has been determined by Equation ( 17), for both thrust and bending moment calculations: The variables used for thrust and bending moment estimations, and the corresponding dimensions, are shown in Table 1.All of the calculations are based on the 5 MW reference wind turbine described in [28], where R is the rotor radius.Figures 7 and 8 illustrate the comparison between the clean case and the ID25 case for the thrust and bending moment results at the wind speed realizations of NTM5 and NTM10, respectively.The thrust force is represented by the green   After applying the BEM algorithm to all control volumes, the tangential and normal load distributions are known and global parameters, such as thrust and bending moment at the root of the blade, can be computed.The mean values of both the thrust and bending moment have been calculated for each wind speed realization.The thrust T has been calculated by Equation ( 15) taking into account the thrust distribution along the blade and the bending moment M is determined by Equation ( 16): The Prandtl's tip loss correction factor F has been determined by Equation ( 17), for both thrust and bending moment calculations: Energies 2017, 10, 742 12 of 15 The variables used for thrust and bending moment estimations, and the corresponding dimensions, are shown in Table 1.All of the calculations are based on the 5 MW reference wind turbine described in [28], where R is the rotor radius.Figures 7 and 8 illustrate the comparison between the clean case and the ID25 case for the thrust and bending moment results at the wind speed realizations of NTM5 and NTM10, respectively.The thrust force is represented by the green color line and the bending moment by the magenta line.An enlargement view has been included in the bottom of each plot in order to show a detailed view of the variation between the clean case and the ID25 case.At both wind speed realizations the increase in the thrust force is notable due to the flow control devices implemented in the ID25 case.As expected, the bending moment in the root of the blade is also increased, but at a lower scale.
Energies 2017, 10, 742 12 of 15 color line and the bending moment by the magenta line.An enlargement view has been included in the bottom of each plot in order to show a detailed view of the variation between the clean case and the ID25 case.At both wind speed realizations the increase in the thrust force is notable due to the flow control devices implemented in the ID25 case.As expected, the bending moment in the root of the blade is also increased, but at a lower scale.Table 6 represents the mean values of thrust and blade root bending moment calculations for the clean case and for the ID25 case.The thrust was calculated by integrating the thrust distribution along the blade by Equation (15) for each wind speed realization.Afterwards, the average thrust force was determined according to the wind speed realization duration in Figure 3.The increments in the thrust force of case ID25 with respect to the clean case were estimated to be 1.639% and 1.473% at the average wind speeds of NTM5 and NTM10, respectively.As expected, those increments of the average thrust in the ID25 case led to an augmentation of the bending moment in the root of the blade.The augmentations in the average bending moment for both wind speed realizations are represented in the last column of Table 6 and were estimated to be less than 1% in comparison against the clean case.color line and the bending moment by the magenta line.An enlargement view has been included in the bottom of each plot in order to show a detailed view of the variation between the clean case and the ID25 case.At both wind speed realizations the increase in the thrust force is notable due to the flow control devices implemented in the ID25 case.As expected, the bending moment in the root of the blade is also increased, but at a lower scale.Table 6 represents the mean values of thrust and blade root bending moment calculations for the clean case and for the ID25 case.The thrust was calculated by integrating the thrust distribution along the blade by Equation (15) for each wind speed realization.Afterwards, the average thrust force was determined according to the wind speed realization duration in Figure 3.The increments in the thrust force of case ID25 with respect to the clean case were estimated to be 1.639% and 1.473% at the average wind speeds of NTM5 and NTM10, respectively.As expected, those increments of the average thrust in the ID25 case led to an augmentation of the bending moment in the root of the blade.The augmentations in the average bending moment for both wind speed realizations are represented in the last column of Table 6 and were estimated to be less than 1% in comparison against the clean case.Table 6 represents the mean values of thrust and blade root bending moment calculations for the clean case and for the ID25 case.The thrust was calculated by integrating the thrust distribution along the blade by Equation (15) for each wind speed realization.Afterwards, the average thrust force was determined according to the wind speed realization duration in Figure 3.The increments in the thrust force of case ID25 with respect to the clean case were estimated to be 1.639% and 1.473% at the average wind speeds of NTM5 and NTM10, respectively.As expected, those increments of the average thrust in the ID25 case led to an augmentation of the bending moment in the root of the blade.The augmentations in the average bending moment for both wind speed realizations are represented in the last column of Table 6 and were estimated to be less than 1% in comparison against the clean case.

Conclusions
BEM-based computations have been carried out to investigate the effects of GFs and VGs on the power performance of a 5 MW wind turbine.A new BEM-based algorithm has been developed and all of the necessary equations have been computed based on the classical BEM method.This improved method has been designed to avoid numerical oscillations or instabilities.
An overall increase on the average wind turbine power output has been found in the current study due to the implementation of GFs and VGs at different blade stations.At lower average wind speeds of NTM5 the increment in the power output wanders around 10% for all the cases illustrated in Table 5.That increase decays to 3% at NTM10 of average wind, which is still significant.The best results, in terms of average power, are reached by the case denoted by ID25, which is a combination of the best cases found from ID1 to ID24 cases.The rise in the average wind turbine power production is 10.4%, in comparison with the clean wind turbine at a wind speed of NTM5 and 3.5% at NTM10.
The wind turbine rotor thrust has also been computed at two wind speed realizations for the ID25 case in comparison against the clean wind turbine.An increase of 1.6% at NTM5 and an increase of 1.5% at NTM10 in the ID25 case average thrust have been achieved in comparison with the thrust force of the clean wind turbine.The enhancement is more prominent at low average wind speeds than at the speeds close to the rated power.
As expected, the increase in the wind turbine power output due to the passive flow control device implementation leads to an augmentation of the average bending moment in the root of the blade of the ID25 case in comparison against the clean wind turbine.At the wind speed realization of NTM5 the increase was estimated as 0.9% and as 0.8% for NTM10.However, that enhancement in the bending moment is acceptable, taking into account the rise in the average power output production of 10.4% and 3.5% for both wind speed realizations, NMT5 and NTM10, respectively.
The results of the current study show that careful analysis of the flow control device type and location on the airfoil, combined with a selection of an appropriate spanwise location on the blade, can yield an effective device flow control system to increase the wind turbine power output with a low penalty in the root bending moment.

Figure 1 .
Figure 1.A prismatic airfoil equipped with vortex generators (VGs).(a) Row of triangular VGs; (b) A detailed view of the VG pairs.

Figure 1 .
Figure 1.A prismatic airfoil equipped with vortex generators (VGs).(a) Row of triangular VGs; (b) A detailed view of the VG pairs.

( 1 )
Initialization by guessing the values of a and a , and the axial and tangential induction factors, respectively.(2) Calculate the flow angle Φ i .(3) Calculate the local angle of attack α. (4) Read off C L (α) and C D (α).(5) Compute the normal C n and tangential C t load coefficients.(6) Re-calculate a and a .

Figure 2 .
Figure 2. Schematic flowchart of the improved BEM algorithm with the recursive function of Equation (11) for the flow angle variable instead of a and a'.

Figure 2 .
Figure 2. Schematic flowchart of the improved BEM algorithm with the recursive function of Equation (11) for the flow angle variable Φ instead of a and a .

Figure 3 .
Figure 3. Wind speed realization calculated with Turbsim.(a) NTM at 5 m/s of average wind speed; (b) NTM at 10 m/s of average wind speed.

Figure 3 .
Figure 3. Wind speed realization calculated with Turbsim.(a) NTM at 5 m/s of average wind speed; (b) NTM at 10 m/s of average wind speed.

Figure 3 .
Figure 3. Wind speed realization calculated with Turbsim.(a) NTM at 5 m/s of average wind speed; (b) NTM at 10 m/s of average wind speed.

Figure 5 .Figure 6 .
Figure 5.Comparison of the (a) power curves and (b) Cp calculations of the clean turbine with the best case ID25 and the calculations made by Jonkman [28].

Figure 5 .Figure 5 .Figure 6 .
Figure 5.Comparison of the (a) power curves and (b) C p calculations of the clean turbine with the best case ID25 and the calculations made by Jonkman [28].

Figure 6 .
Figure 6.Wind turbine power distribution with time.(a) The wind speed realization of NTM5; (b) The wind speed realization of NTM10.

Figure 7 .Figure 8 .
Figure 7. Wind turbine performance at NTM5 wind speed realization.(a) Thrust force, and (b) bending moment at the root of the blade

Figure 7 .
Figure 7. Wind turbine performance at NTM5 wind speed realization.(a) Thrust force, and (b) bending moment at the root of the blade.

Figure 7 .Figure 8 .
Figure 7. Wind turbine performance at NTM5 wind speed realization.(a) Thrust force, and (b) bending moment at the root of the blade

Figure 8 .
Figure 8. Wind turbine performance at NTM10 wind speed realization.(a) Thrust force; (b) Bending moment at the root of the blade.

Table 1 .
Description of current variables.

Table 3 .
Types of flow control devices and locations on the airfoils.

Table 4 .
Cases studied according to flow control device span distribution.St.: station.

Table 5 .
Average wind turbine power output for the clean wind turbine in comparison with the cases with flow control devices.Calculations have been made at two different wind speed realizations, NTM5 and NTM10.∆ indicates the average power variation of each ID case with respect to the clean case.

Table 6 .
Thrust force and bending moment average values for the clean and the ID25 cases at two different wind speed realizations.∆T and ∆M indicate the variation of thrust and bending moment, respectively, with respect to the clean case, in %.