Microtab Design and Implementation on a 5 MW Wind Turbine

Microtabs (MT) consist of a small tab placed on the airfoil surface close to the trailing edge and perpendicular to the surface. A study to find the optimal position to improve airfoil aerodynamic performance is presented. Therefore, a parametric study of a MT mounted on the pressure surface of an airfoil has been carried out. The aim of the current study is to find the optimal MT size and location to increase airfoil aerodynamic performance and to investigate its influence on the power output of a 5 MW wind turbine. Firstly, a computational study of a MT mounted on the pressure surface of the airfoil DU91W(2)250 has been carried out and the best case has been found according to the largest lift-to-drag ratio. This airfoil has been selected because it is typically used on wind turbine, such as the 5 MW reference wind turbine of the National Renewable Energy Laboratory (NREL). Second, Blade Element Momentum (BEM) based computations have been performed to investigate the effect of the MT on the wind turbine power output with different wind speed realizations. The results show that, due to the implementation of MTs, a considerable increase in the turbine average power is achieved.


Introduction
In the field of energy production, wind energy is a key issue in order to reduce fossil fuel dependency. In addition, the solar-wind hybrid energy system has become very popular (Bouzelata et al. [1]). Wind energy is an essential resource among the other clean energy production methods. The search for an energy power policy that is local, sustainable, and environmentally friendly, and optimizes resources has become a requirement. Therefore, models that include factors such as emissions reduction, minimization of imported energy, and even social acceptance are proposed in many studies, such as Novosel et al. [2] and Kumu et al. [3]. Lately research, has been focused on wind turbine blade improvements to optimize rotor dynamic behavior (Jaume et al. [4] and Vaz et al. [5]).
The yearly 9% increase of installed wind energy in Europe in the last fifteen years shows the significance of research in the field of flow control for large wind turbines (Houghton et al. [6]). The considerable growth of wind turbine rotor size and weight in the last decade has made it impossible to control as they were controlled 30 years ago. Rotors of 120 meters or even more are now a reality.
Johnson et al. [7] compiled some of the most important load control techniques that could be used in wind turbines to assure a safe and optimal operation under a variety of atmospheric conditions. The improvements of present wind turbines should be pointed to minimize the fatigue from the wind turbine rotor and other structural components due to changes in wind direction, speed and turbulence, as well as start-stop operations of the wind turbine and, of course, to maximize the energy production.
In the last decades, many different flow control devices have been designed and developed (see Chow et al. [8]). Most of them were intended for aeronautical applications (Taylor [9]), but also frequently used in turbo machinery (see Liu et al. [10] and Xu et al. [11]). Currently, researchers are working to optimize and introduce this type of devices in horizontal axis wind turbines (HAWT). Wood [12] developed a four layer scheme which allows classifying the different concepts that are part of all flow control devices. In the study of Shires et al. [13], 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. In addition, the dynamic stall control was investigated by Xu et al. [14] on a S809 airfoil by the implementation of a co-flow jet.
Depending on their operating principle, they can be classified as actives or passives (see Aramendia et al. [15]). Passive control techniques would represent an improvement in the turbine's efficiency and in loads reduction without external energy consumption. Active control techniques need an additional energy source to get the desired effect on the flow and, unlike microtabs and other passive devices, active flow control needs intricate algorithms to get the maximum benefit (see Becker et al. [16] and Macquart et al. [17]). Johnson et al. [7] made an analysis and discussed fifteen different devices for wind turbine control. Some of them are still being tested on full-scale turbines.
The microtabs (MTs) consist of a small tabs situated close to the trailing edge (TE) of an airfoil, which projects perpendicular to the surface of the airfoil a few percent of the chord length c (usually 1-2%) corresponding to the boundary layer thickness. The potential of MTs was first investigated by van Dam et al. [18]. Baker et al. [19] carried out broad study dedicated to the S809 airfoil with MTs. These MTs jets the flow in the boundary layer away from the blade's surface, bringing a recirculation zone behind the tab, as can be observed in Figure 1. The MTs affect the airfoil aerodynamics shifting the point of flow separation and, therefore, providing changes in lift. Lift improvement is obtained by implementing the MT downwards (on the pressure side) and lift reduction is obtained by deploying the MT upwards (on the suction side).
Appl. Sci. 2017, 7, 536 2 of 19 now a reality. Johnson et al. [7] compiled some of the most important load control techniques that could be used in wind turbines to assure a safe and optimal operation under a variety of atmospheric conditions. The improvements of present wind turbines should be pointed to minimize the fatigue from the wind turbine rotor and other structural components due to changes in wind direction, speed and turbulence, as well as start-stop operations of the wind turbine and, of course, to maximize the energy production.
In the last decades, many different flow control devices have been designed and developed (see Chow et al. [8]). Most of them were intended for aeronautical applications (Taylor [9]), but also frequently used in turbo machinery (see Liu et al. [10] and Xu et al. [11]). Currently, researchers are working to optimize and introduce this type of devices in horizontal axis wind turbines (HAWT). Wood [12] developed a four layer scheme which allows classifying the different concepts that are part of all flow control devices. In the study of Shires et al. [13], 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. In addition, the dynamic stall control was investigated by Xu et al. [14] on a S809 airfoil by the implementation of a co-flow jet.
Depending on their operating principle, they can be classified as actives or passives (see Aramendia et al. [15]). Passive control techniques would represent an improvement in the turbine's efficiency and in loads reduction without external energy consumption. Active control techniques need an additional energy source to get the desired effect on the flow and, unlike microtabs and other passive devices, active flow control needs intricate algorithms to get the maximum benefit (see Becker et al. [16] and Macquart et al. [17]). Johnson et al. [7] made an analysis and discussed fifteen different devices for wind turbine control. Some of them are still being tested on full-scale turbines.
The microtabs (MTs) consist of a small tabs situated close to the trailing edge (TE) of an airfoil, which projects perpendicular to the surface of the airfoil a few percent of the chord length c (usually 1-2%) corresponding to the boundary layer thickness. The potential of MTs was first investigated by van Dam et al. [18]. Baker et al. [19] carried out broad study dedicated to the S809 airfoil with MTs. These MTs jets the flow in the boundary layer away from the blade's surface, bringing a recirculation zone behind the tab, as can be observed in Figure 1.The MTs affect the airfoil aerodynamics shifting the point of flow separation and, therefore, providing changes in lift. Lift improvement is obtained by implementing the MT downwards (on the pressure side) and lift reduction is obtained by deploying the MT upwards (on the suction side). The implementation of this device near the airfoil (TE) provokes changes in the flow, causing modifications in the circulation of the flow on the airfoil. The effective camber of the airfoil is modified, promoting changes in the lift and drag forces. Placing a MT at the pressure surface the lift is increased, however if it is placed on the suction, the opposite effect is reached. The main advantages of MTs are: small size, low cost of manufacturing, low power requirements for its activation, and simplicity of the device design. Multiple studies into this topic were made by van Dam et al. [20] and Yen et al. [21], including wind-tunnel experiments in order to determine their optimal height and location. The implementation of this device near the airfoil (TE) provokes changes in the flow, causing modifications in the circulation of the flow on the airfoil. The effective camber of the airfoil is modified, promoting changes in the lift and drag forces. Placing a MT at the pressure surface the lift is increased, however if it is placed on the suction, the opposite effect is reached. The main advantages of MTs are: small size, low cost of manufacturing, low power requirements for its activation, and simplicity of the device design. Multiple studies into this topic were made by van Dam et al. [20] and Yen et al. [21], including wind-tunnel experiments in order to determine their optimal height and location.
The aim of the current study is to find the optimal MT position to increase DU91W(2)250 airfoil aerodynamic performance and to investigate its influence on the average power output of a 5 MW Wind turbine. Firstly, a parametric study of a MT mounted on the pressure surface of the airfoil DU91W(2)250 has been carried out and best case has been selected. This airfoil has been selected because it is used on the 5 MW reference wind turbine of the NREL, as described in Jonkman et al. [22]. Second, Blade Element Momentum (BEM) based computations have been performed to investigate the effect of this MTs on the wind turbine power output. The results on the rotor thrust and blade root bending moment are also presented.

Numerical Setup
In order to obtain some of the main features of the MT, Computational Fluid Dynamic (CFD) techniques have been employed. Currently, non-commercial and proprietary CFD codes are used to reproduce relatively well any physical problem. In this work, the open source code OpenFOAM has been used for simulating the effects of a MT on a DU91W (2)250. This open source CFD code is an object-oriented library written in C++ to solve computational continuum mechanics problems. One of its advantages is that the user can modify the code to create new solvers and applications as well as freely share the code developed.
The SIMPLE algorithm was employed for the pressure-velocity coupling. The convective terms were discretized with a second order linear-upwind scheme. The discretization of the viscous terms was achieved by means of a second order central-differences linear scheme. The simulations were run fully turbulent. Steady state simulations were carried out and performed with a structured finite-volume flow solver using Reynolds averaged Navier-Stokes (RANS) equations. For these computations the k-ω SST shear stress turbulence model developed by Menter [23] was used due to its superior separated flow performance, as reported by Kral [24] and Gatski [25]. The model is a combination of two models: Wilcox's k-ω model for near wall regions and the k-ε model for the outer region and in free shear flows. The SST model departs from existing k-ω and k-ε models by way of a modified eddy viscosity definition that results in improved prediction of separated flows. Reynolds-averaged Navier-Stokes calculations with SST turbulence model for various tab configurations applied to an airfoil are presented in Mayda et al. [26]. Figure 2 illustrates the computational setup with the current setting consisting of a DU airfoil. An O-mesh type was designed for the computations with a computational domain radius of 32 times the airfoil chord length R = 32c, which is in the order of the computational size recommended by Sørensen et al. [27] for this type of simulations. The aim of the current study is to find the optimal MT position to increase DU91W(2)250 airfoil aerodynamic performance and to investigate its influence on the average power output of a 5 MW Wind turbine. Firstly, a parametric study of a MT mounted on the pressure surface of the airfoil DU91W(2)250 has been carried out and best case has been selected. This airfoil has been selected because it is used on the 5 MW reference wind turbine of the NREL, as described in Jonkman et al. [22]. Second, Blade Element Momentum (BEM) based computations have been performed to investigate the effect of this MTs on the wind turbine power output. The results on the rotor thrust and blade root bending moment are also presented.

Numerical Setup
In order to obtain some of the main features of the MT, Computational Fluid Dynamic (CFD) techniques have been employed. Currently, non-commercial and proprietary CFD codes are used to reproduce relatively well any physical problem. In this work, the open source code OpenFOAM has been used for simulating the effects of a MT on a DU91W (2)250. This open source CFD code is an object-oriented library written in C++ to solve computational continuum mechanics problems. One of its advantages is that the user can modify the code to create new solvers and applications as well as freely share the code developed.
The SIMPLE algorithm was employed for the pressure-velocity coupling. The convective terms were discretized with a second order linear-upwind scheme. The discretization of the viscous terms was achieved by means of a second order central-differences linear scheme. The simulations were run fully turbulent. Steady state simulations were carried out and performed with a structured finite-volume flow solver using Reynolds averaged Navier-Stokes (RANS) equations. For these computations the k-ω SST shear stress turbulence model developed by Menter [23] was used due to its superior separated flow performance, as reported by Kral [24] and Gatski [25]. The model is a combination of two models: Wilcox's k-ω model for near wall regions and the k-ε model for the outer region and in free shear flows. The SST model departs from existing k-ω and k-ε models by way of a modified eddy viscosity definition that results in improved prediction of separated flows. Reynolds-averaged Navier-Stokes calculations with SST turbulence model for various tab configurations applied to an airfoil are presented in Mayda et al. [26]. Figure 2 illustrates the computational setup with the current setting consisting of a DU airfoil. An O-mesh type was designed for the computations with a computational domain radius of 32 times the airfoil chord length R = 32c, which is in the order of the computational size recommended by Sørensen et al. [27] for this type of simulations.  The Reynolds number based on the airfoil chord length of c = 1 m is Re = 7 × 10 6 . The computational setup of the simulations consists of a structured mesh with the first cell height ∆z/c of 1.45 × 10 −6 normalized by the airfoil chord length. The stretching in the chord-wise and normal directions is accomplished by double sided tanh stretching functions based on Vinokur [28] and Thompson et al. [29]. The mesh domain was designed to have a dimensionless distance less than 1 (y+ < 1) on the airfoil wall. An optimized mesh plays a major role in the CFD simulations as the tool that will help the user to discretize the domain. It is important to identify the mesh regions where the results have to be quite accurate as well as to establish a balance between the accuracy of the simulations and the computational cost. Figure 3 shows the cell distribution around the MT and in the near wake of the trailing edge. The wake and the regions where high gradients were expected were accordingly refined. There are certain regions close to the MT in which the velocity gradient changes drastically, which is the reason why those areas are so important. directions is accomplished by double sided tanh stretching functions based on Vinokur [28] and Thompson et al. [29]. The mesh domain was designed to have a dimensionless distance less than 1 (y+ < 1) on the airfoil wall. An optimized mesh plays a major role in the CFD simulations as the tool that will help the user to discretize the domain. It is important to identify the mesh regions where the results have to be quite accurate as well as to establish a balance between the accuracy of the simulations and the computational cost. Figure 3 shows the cell distribution around the MT and in the near wake of the trailing edge. The wake and the regions where high gradients were expected were accordingly refined. There are certain regions close to the MT in which the velocity gradient changes drastically, which is the reason why those areas are so important.  In OpenFoam, both the velocity and the pressure conditions have to be defined at all boundaries, since the velocity-pressure coupling is based on a collocated grid approach. Non-slip boundary condition was set for the airfoil and MT walls. Computational simulations of the DU91(2)250 airfoil without any MT have been carried out and validated against the data obtained by Xfoil from DOWEC project of Kooijman et al. [30] and Lindenburg [31]. The Lift-to-drag ratio was calculated for ten angles of attack (AoA) from α = 0 to α = 9 degrees. Figure 4 shows the results of the CFD computations against the Xfoil results for all angles of attacks of the airfoil. A mesh independency study was carried out to verify enough grid resolution with three mesh sizes using a refinement ratio of 2. The coarse mesh contains 72,500 cells. For the medium and fine mesh the number of cells is 145,000 and 290,000, respectively. Drag and lift results obtained for the finer mesh are compared with the results of a regular and a coarse mesh. Less than 4% mesh dependency has been found for both drag and lift. The simulations were converged until a satisfactory residual convergence was achieved on the velocities, pressure and turbulence quantities.  In OpenFoam, both the velocity and the pressure conditions have to be defined at all boundaries, since the velocity-pressure coupling is based on a collocated grid approach. Non-slip boundary condition was set for the airfoil and MT walls. Computational simulations of the DU91(2)250 airfoil without any MT have been carried out and validated against the data obtained by Xfoil from DOWEC project of Kooijman et al. [30] and Lindenburg [31]. The Lift-to-drag ratio was calculated for ten angles of attack (AoA) from α = 0 to α = 9 degrees. Figure 4 shows the results of the CFD computations against the Xfoil results for all angles of attacks of the airfoil. directions is accomplished by double sided tanh stretching functions based on Vinokur [28] and Thompson et al. [29]. The mesh domain was designed to have a dimensionless distance less than 1 (y+ < 1) on the airfoil wall. An optimized mesh plays a major role in the CFD simulations as the tool that will help the user to discretize the domain. It is important to identify the mesh regions where the results have to be quite accurate as well as to establish a balance between the accuracy of the simulations and the computational cost. Figure 3 shows the cell distribution around the MT and in the near wake of the trailing edge. The wake and the regions where high gradients were expected were accordingly refined. There are certain regions close to the MT in which the velocity gradient changes drastically, which is the reason why those areas are so important.  In OpenFoam, both the velocity and the pressure conditions have to be defined at all boundaries, since the velocity-pressure coupling is based on a collocated grid approach. Non-slip boundary condition was set for the airfoil and MT walls. Computational simulations of the DU91(2)250 airfoil without any MT have been carried out and validated against the data obtained by Xfoil from DOWEC project of Kooijman et al. [30] and Lindenburg [31]. The Lift-to-drag ratio was calculated for ten angles of attack (AoA) from α = 0 to α = 9 degrees. Figure 4 shows the results of the CFD computations against the Xfoil results for all angles of attacks of the airfoil. A mesh independency study was carried out to verify enough grid resolution with three mesh sizes using a refinement ratio of 2. The coarse mesh contains 72,500 cells. For the medium and fine mesh the number of cells is 145,000 and 290,000, respectively. Drag and lift results obtained for the finer mesh are compared with the results of a regular and a coarse mesh. Less than 4% mesh dependency has been found for both drag and lift. The simulations were converged until a satisfactory residual convergence was achieved on the velocities, pressure and turbulence quantities.  A mesh independency study was carried out to verify enough grid resolution with three mesh sizes using a refinement ratio of 2. The coarse mesh contains 72,500 cells. For the medium and fine mesh the number of cells is 145,000 and 290,000, respectively. Drag and lift results obtained for the finer mesh are compared with the results of a regular and a coarse mesh. Less than 4% mesh dependency has been found for both drag and lift. The simulations were converged until a satisfactory residual convergence was achieved on the velocities, pressure and turbulence quantities. The CFD results follow reasonable good the trend of DOWEC results. Lift and drag coefficients were calculated according to the Equations (1) and (2), respectively: The air density was defined by ρ = 1.204 kg/m 3 and the free stream velocity, far ahead of the airfoil, corresponds to U ∞ = 10.66 m/s. L and D represent the lift and drag forces per unit of area, since the simulations are in 2D.

Microtab Lay-Out
The MT position in the airfoil is sketched in Figures 5 and 6. Dimension x represents the position from the LE and y represents the height of the MT, both in percentage of c. Twelve cases have been established depending on the distance measured respective to the LE in %c (see Table 1). These cases are: 93%c, 94%c, 95%c and 96%c. The MT height relative to the chord length measured in percentage is 1%c, 1.5%c and 2%c. These series of cases has been designed according to the previous studies of Standish et al. [32], Mayda et al. [26] and Yen et al. [21], where the maximum translation was estimated in the order of the boundary layer thickness at the device position: 1-2% of the chord length and the optimal location for a lower surface tab in terms of lift and drag was found to be around 95% of c. The MTs are placed on the pressure surface and have been studied for ten different angles of attack, from 0 • to 9 • . The combination of all these positions for the MTs gives 120 different cases to study (Ayerdi-Zaton et al. [33]). The airfoil DU91W(2)250 without any flow control device was also simulated to study the influence of the previously described MTs. The CFD results follow reasonable good the trend of DOWEC results. Lift and drag coefficients were calculated according to the Equations (1) and (2), respectively: The air density was defined by = 1.204 kg/m 3 and the free stream velocity, far ahead of the airfoil, corresponds to U∞ = 10.66 m/s. L and D represent the lift and drag forces per unit of area, since the simulations are in 2D.

Microtab Lay-Out
The MT position in the airfoil is sketched in Figures 5 and 6. Dimension x represents the position from the LE and y represents the height of the MT, both in percentage of c. Twelve cases have been established depending on the distance measured respective to the LE in %c (see Table 1). These cases are: 93%c, 94%c, 95%c and 96%c. The MT height relative to the chord length measured in percentage is 1%c, 1.5%c and 2%c. These series of cases has been designed according to the previous studies of Standish et al. [32], Mayda et al. [26] and Yen et al. [21], where the maximum translation was estimated in the order of the boundary layer thickness at the device position: 1-2% of the chord length and the optimal location for a lower surface tab in terms of lift and drag was found to be around 95% of c. The MTs are placed on the pressure surface and have been studied for ten different angles of attack, from 0° to 9°. The combination of all these positions for the MTs gives 120 different cases to study (Ayerdi-Zaton et al. [33]). The airfoil DU91W(2)250 without any flow control device was also simulated to study the influence of the previously described MTs.    The CFD results follow reasonable good the trend of DOWEC results. Lift and drag coefficients were calculated according to the Equations (1) and (2), respectively: The air density was defined by = 1.204 kg/m 3 and the free stream velocity, far ahead of the airfoil, corresponds to U∞ = 10.66 m/s. L and D represent the lift and drag forces per unit of area, since the simulations are in 2D.

Microtab Lay-Out
The MT position in the airfoil is sketched in Figures 5 and 6. Dimension x represents the position from the LE and y represents the height of the MT, both in percentage of c. Twelve cases have been established depending on the distance measured respective to the LE in %c (see Table 1). These cases are: 93%c, 94%c, 95%c and 96%c. The MT height relative to the chord length measured in percentage is 1%c, 1.5%c and 2%c. These series of cases has been designed according to the previous studies of Standish et al. [32], Mayda et al. [26] and Yen et al. [21], where the maximum translation was estimated in the order of the boundary layer thickness at the device position: 1-2% of the chord length and the optimal location for a lower surface tab in terms of lift and drag was found to be around 95% of c. The MTs are placed on the pressure surface and have been studied for ten different angles of attack, from 0° to 9°. The combination of all these positions for the MTs gives 120 different cases to study (Ayerdi-Zaton et al. [33]). The airfoil DU91W(2)250 without any flow control device was also simulated to study the influence of the previously described MTs.

Computational Results
A parametric study has been carried out in order to find the optimal position of the MT on the airfoil DU91W(2)250. Table 1 illustrates the 13 cases studied in the current work. The first case studied is the one with no MT implemented and the other 12 cases with different sizes (y) and position of the MT from the leading edge of the airfoil (x). Each case has been studied for ten different degrees of angle of attack α, in the range from 0 • to 9 • . Figure 7 illustrates the lift-to-drag ratio C L /C D evolution for every angle α and all MT cases. In the left column the evolution of the C L /C D vs. the location of the MT from the airfoil leading edge x is represented. Right column plots illustrate the C L /C D evolution against the three different heights of the MT y = 1, 1.5 and 2%c. Both x and y parameters are represented in terms of percent of the airfoil chord length c. The Lift-to-drag ratio of the airfoil without MT for each AoA is represented by the black continuous line in each plot.

Computational Results
A parametric study has been carried out in order to find the optimal position of the MT on the airfoil DU91W(2)250. Table 1 illustrates the 13 cases studied in the current work. The first case studied is the one with no MT implemented and the other 12 cases with different sizes (y) and position of the MT from the leading edge of the airfoil (x). Each case has been studied for ten different degrees of angle of attack α, in the range from 0° to 9°. Figure 7 illustrates the lift-to-drag ratio CL/CD evolution for every angle α and all MT cases. In the left column the evolution of the CL/CD vs. the location of the MT from the airfoil leading edge x is represented. Right column plots illustrate the CL/CD evolution against the three different heights of the MT y = 1, 1.5 and 2%c. Both  For the range of angles of attack, from 0° to 9°, the best values of the CL/CD ratio are reached by the cases MT9510, MT9515 and MT9520. However, as can be seen in the left column plots of Figure 7, the highest values are reached around the x = 95% of the chord length. On the right column plots, it is clear again that the best values of the Lift-to-drag ratio are reached by x = 95% of c, and the highest ones by the MT height of y = 2% of c. Therefore, in the range of angles of attack studied in the present work, the best case in terms of Lift-to-drag ratio is the one defined by: x = 95% and y = 2% of c. These results are in concordance previous studies of Yen et al. [21] where it was found that the best place to situate the lower surface tab with respect to lift and drag was around 95%c. According to the classification of Table 1, that case corresponds to DU912250MT9520. The CL/CD ratio achieved in that case, keeps up to the ratio of the clean airfoil represented in Figure 7 with a continuous black line. Figure 8 illustrates the Lift-to-drag ratio values of the case DU912250MT9520 in comparison with the ratios obtained for the clean airfoil.
At low AoAs, the increase in the CL/CD ratio due to the microtab MT9520 implementation is clearly visible. However, at the AoA of 9° the CL/CD ratio stays below the CL/CD ratio of the clean airfoil DU91W(2)250. Note that this behavior is repeated by the other cases with different MTs configurations, as shown in Figure 7e. The reason could be found in the fact that at 9° of AoA the airfoil is working near the stall conditions and the CD increases considerably. Since the computations were performed in 2D, three-dimensional effects were neglected. In the studies of Mayda et al. [26] and Zahle et al. [34], a detailed study on the three-dimensional shedding and spanwise flow can be found. Figure 9 represents the streamwise velocity distribution around the MT of the case with the best Lift-to-drag ratio: DU91W(2)250 MT9520. The presence of the tab changes the trailing edge flow development, the so-called Kutta condition, and consequently the effective camber of the airfoils is modified, providing in this case lift enhance. The MT jets the flow in the BL away from the airfoil surface, producing a circulation region behind the tab. For the range of angles of attack, from 0 • to 9 • , the best values of the C L /C D ratio are reached by the cases MT9510, MT9515 and MT9520. However, as can be seen in the left column plots of Figure 7, the highest values are reached around the x = 95% of the chord length. On the right column plots, it is clear again that the best values of the Lift-to-drag ratio are reached by x = 95% of c, and the highest ones by the MT height of y = 2% of c. Therefore, in the range of angles of attack studied in the present work, the best case in terms of Lift-to-drag ratio is the one defined by: x = 95% and y = 2% of c. These results are in concordance previous studies of Yen et al. [21] where it was found that the best place to situate the lower surface tab with respect to lift and drag was around 95%c. According to the classification of Table 1, that case corresponds to DU912250MT9520. The C L /C D ratio achieved in that case, keeps up to the ratio of the clean airfoil represented in Figure 7 with a continuous black line. Figure 8 illustrates the Lift-to-drag ratio values of the case DU912250MT9520 in comparison with the ratios obtained for the clean airfoil.
At low AoAs, the increase in the C L /C D ratio due to the microtab MT9520 implementation is clearly visible. However, at the AoA of 9 • the C L /C D ratio stays below the C L /C D ratio of the clean airfoil DU91W(2)250. Note that this behavior is repeated by the other cases with different MTs configurations, as shown in Figure 7e. The reason could be found in the fact that at 9 • of AoA the airfoil is working near the stall conditions and the C D increases considerably. Since the computations were performed in 2D, three-dimensional effects were neglected. In the studies of Mayda et al. [26] and Zahle et al. [34], a detailed study on the three-dimensional shedding and spanwise flow can be found. Figure 9 represents the streamwise velocity distribution around the MT of the case with the best Lift-to-drag ratio: DU91W(2)250 MT9520. The presence of the tab changes the trailing edge flow development, the so-called Kutta condition, and consequently the effective camber of the airfoils is A comparison of the pressure distribution on the surface of the clean airfoil DU91W(2)250 and the airfoil with the best aerodynamic performance in terms of lit-to-drag ratio DU91W(2)250MT9520 is presented in Figure 10. The presence of the MT considerably increases the aft loading of the airfoil and a positive gap between the clean and the microtabed airfoil is clearly visible at all angles of attack. The airfoil shape is sketched by a continuous black line. A comparison of the pressure distribution on the surface of the clean airfoil DU91W(2)250 and the airfoil with the best aerodynamic performance in terms of lit-to-drag ratio DU91W(2)250MT9520 is presented in Figure 10. The presence of the MT considerably increases the aft loading of the airfoil and a positive gap between the clean and the microtabed airfoil is clearly visible at all angles of attack. The airfoil shape is sketched by a continuous black line. A comparison of the pressure distribution on the surface of the clean airfoil DU91W(2)250 and the airfoil with the best aerodynamic performance in terms of lit-to-drag ratio DU91W(2)250MT9520 is Appl. Sci. 2017, 7, 536 9 of 18 presented in Figure 10. The presence of the MT considerably increases the aft loading of the airfoil and a positive gap between the clean and the microtabed airfoil is clearly visible at all angles of attack. The airfoil shape is sketched by a continuous black line.

Wind Speed Model
The wind speed realizations used in the current study, as shown in Figures 11, have been calculated with the TurbSim tool (Kelley et al. [35]). The turbulence model is the Normal Turbulence Model (NTM) following the IEC 61400 norm and the wind speed series have been generated with the following parameters: TurbSim uses an adapted version of Veers [36] to generate time series based on spectral representation. The IECKAI (IEC Kaimal) model is defined in IEC 61400-1 2nd ed. [37] and 3rd ed. [38], and assumes neutral atmospheric stability. The spectra for the three wind components, K = u, v, w, are calculated by Equation (3): where f is the cyclic frequency and LK is an integral scale parameter defined in IEC 61400-1 standard 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 11 represents the wind speed series used in the present study according to the Normal Turbulence model with 5 m/s, 7.5 m/s and with 10 m/s average velocity, respectively. According to TurbSim user specifications, the first input value is a random seed that must be an integer between −2,147,483,648 and 2,147,483,647 (inclusive). In the current study, three different seeds have been chosen for each wind realizations. Figure 11a 0deg no MT 2deg no MT 4deg no MT 6deg no MT 9deg no MT 0deg MT9520 2deg MT9520 4deg MT9520 6deg MT9520 9deg MT9520

Wind Speed Model
The wind speed realizations used in the current study, as shown in Figure 11, have been calculated with the TurbSim tool (Kelley et al. [35]). The turbulence model is the Normal Turbulence Model (NTM) following the IEC 61400 norm and the wind speed series have been generated with the following parameters: TurbSim uses an adapted version of Veers [36] to generate time series based on spectral representation. The IECKAI (IEC Kaimal) model is defined in IEC 61400-1 2nd ed. [37] and 3rd ed. [38], and assumes neutral atmospheric stability. The spectra for the three wind components, K = u, v, w, are calculated by Equation (3): where f is the cyclic frequency and L K is an integral scale parameter defined in IEC 61400-1 standard. 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 11 represents the wind speed series used in the present study according to the Normal Turbulence model with 5 m/s, 7.5 m/s and with 10 m/s average velocity, respectively. According to TurbSim user specifications, the first input value is a random seed that must be an integer between −2,147,483,648 and 2,147,483,647 (inclusive). In the current study, three different seeds have been chosen for each wind realizations. Figure 11a-c illustrates the wind patterns generated for each wind speed. The values of Seeds 1-3 have been chosen to obtain different wind patterns and are the same for each average wind speed. These different wind speed realizations were chosen to investigate the effects of the MTs, since it is a good way to evaluate the wind turbine power output at low and medium wind velocities. Four different cases were considered in the current study. The clean wind turbine was taken as the baseline case, without any device implemented and named DU91W(2)250. According to the matrix presented in Table 2, the cases are different depending on the blade station where the passive devices were implemented. The suffix st means the blade station where the MTs were introduced. According to the airfoil distribution described in [22], stations 8 and 9 were chosen for the present study.

Methodology
The primary tools used in the current work to investigate the effects of the passive MTs on the NREL 5 MW Baseline Wind Turbine are engineering models. The NREL 5 MW reference wind turbine is widely used in research studies in the wind energy field since it represents a baseline of the modern and future offshore HAWT. 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 concept of a 5 MW wind turbine is based on the data from the These different wind speed realizations were chosen to investigate the effects of the MTs, since it is a good way to evaluate the wind turbine power output at low and medium wind velocities. Four different cases were considered in the current study. The clean wind turbine was taken as the baseline case, without any device implemented and named DU91W(2)250. According to the matrix presented in Table 2, the cases are different depending on the blade station where the passive devices were implemented. The suffix st means the blade station where the MTs were introduced. According to the airfoil distribution described in [22], stations 8 and 9 were chosen for the present study.

Methodology
The primary tools used in the current work to investigate the effects of the passive MTs on the NREL 5 MW Baseline Wind Turbine are engineering models. The NREL 5 MW reference wind turbine is widely used in research studies in the wind energy field since it represents a baseline of the modern and future offshore HAWT. 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 concept of a 5 MW wind turbine is based on the data from the DOWEC study [30,31], with a concept from the UpWind project [39]. The airfoils and chord schedule used in the present work are presented in Table 3 and are the same from NREL [22], also adopted from the DOWEC project. The blade airfoil locations, labeled as r (m) in Table 3, are directed along the blade-pitch axis from the rotor center to the blade cross sections. The DU25 airfoil corresponds to the DU91W(2)250. More detailed information on the DU family of airfoils used in the current work can be found in the study made by Timmer [40]. The reported NREL 5 MW airfoil distribution is shown in Figure 12. DOWEC study [30,31], with a concept from the UpWind project [39]. The airfoils and chord schedule used in the present work are presented in Table 3 and are the same from NREL [22], also adopted from the DOWEC project. The blade airfoil locations, labeled as r (m) in Table 3, are directed along the blade-pitch axis from the rotor center to the blade cross sections. The DU25 airfoil corresponds to the DU91W(2)250. More detailed information on the DU family of airfoils used in the current work can be found in the study made by Timmer [40]. The reported NREL 5 MW airfoil distribution is shown in Figure 12. 61.6333 NACA64XX Figure 12. Sketch of the airfoil distribution along the blade (not to scale) of the NREL 5 MW wind turbine according to [22].
Once the lift and drag coefficients are identified for the airfoils along the blades, it is feasible to compute the force distribution. Global loads such as the power output and the root bending moment of the blade can be found by integrating this distribution along the blade span. It is the principle of the BEM method, which will be derived to compute the induction factors a and a' and thus the loads on a wind turbine. The present procedure is described in the following steps: Figure 12. Sketch of the airfoil distribution along the blade (not to scale) of the NREL 5 MW wind turbine according to [22].
Once the lift and drag coefficients are identified for the airfoils along the blades, it is feasible to compute the force distribution. Global loads such as the power output and the root bending moment of the blade can be found by integrating this distribution along the blade span. It is the principle of the BEM method, which will be derived to compute the induction factors a and a' and thus the loads on a wind turbine. The present procedure is described in the following steps:

1.
First of all, BEM based computations were carried out in order to characterize the dynamical behavior of the NREL 5 MW wind turbine, including Prandtl's tip loss factor and Glauert's correction. The BEM solver was developed and programmed by the authors of the current study based on the numerical iterative approach of Hansen [41]. All the necessary equations were derived and computed based on the steps proposed by the classical blade element momentum method. The usual basic steps for BEM calculations was followed; this is a short schedule: a. Initialization by guessing values of a and a' values, axial and tangential induction factors, respectively. b.
Calculate the flow angle Φ. c.
Calculate the local angle of attack α. d.
Read off C L (α) and C D (α). e.
Compute the normal C n and tangential C t load coefficients. f.
Re-calculate a and a'. g.
State a tolerance for a and a' and if it has changed more than that tolerance, go to (b) or else finish. h.
Compute the local loads.

2.
Following the specifications of the utility scale multi megawatt wind turbine NREL 5 MW baseline described in [22], all the wind turbine rotor properties were introduced as input characteristics. The polar curves of the airfoil with the MT were taken from best case found in Section 2, which is the DU91W(2)250MT9520 with the MT position from the leading edge at 95% of c and the MT height of 2% of c.

3.
The surfaces of power coefficient Cp were calculated for all cases of the present study according to the matrix distribution described in Table 2.

4.
Once the Cp surfaces have been generated, BEM based computations are run for the four cases and the power curve vs. wind speed is calculated to compare the power curve of the clean turbine with the curve of the turbine with the MTs implemented.

5.
Afterwards, the wind speed realizations explained in Section 3 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 device implemented.

Results from BEM Computations
In order to investigate the influence of MTs on the power of the NREL 5 MW reference wind turbine, BEM computations have been carried out following the steps explained in the previous section. The BEM based computations has now been derived and the power has been computed versus the wind speed. Figure 13 illustrates the power curves along the wind speed for the clean case with no passive device implemented into the blade in comparison with the cases with MTs (see Table 2). The power curves of the wind turbine with MTs follow the trends 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 cases with the MT implemented, as shown in the enlargement view embedded in Figure 13. Additionally, the average wind turbine power output has been calculated for the clean wind turbine and for the cases with MTs at the wind speed realizations described in Section 3. Equation (4) shows how this average power is calculated:  vj: is the wind speed according to the realizations shown in Figure 11.  Nbins: number of bins per data.  P(vj): power at the wind speed vj.  N(vj): number of data at the wind speed vj. where the P(vj) has been determined from the data obtained in Figure 13. Table 4 shows the results of the power output calculations according to Equation (4). Firstly, the average power has been calculated for the clean wind turbine for the three wind speed realizations illustrated in Figure 11, without any flow control device mounted on the blade. Afterwards, the average wind turbine power was calculated for all cases of MTs distribution described in Table 2 and compared with the clean turbines power values. The symbols denoted by Δ represent the increment of average power in comparison with the clean turbine.  Additionally, the average wind turbine power output has been calculated for the clean wind turbine and for the cases with MTs at the wind speed realizations described in Section 3. Equation (4) shows how this average power is calculated: v j : is the wind speed according to the realizations shown in Figure 11. where the P(v j ) has been determined from the data obtained in Figure 13. Table 4 shows the results of the power output calculations according to Equation (4). Firstly, the average power has been calculated for the clean wind turbine for the three wind speed realizations illustrated in Figure 11, without any flow control device mounted on the blade. Afterwards, the average wind turbine power was calculated for all cases of MTs distribution described in Table 2 and compared with the clean turbines power values. The symbols denoted by ∆ represent the increment of average power in comparison with the clean turbine. At the wind speed realization of NTM5, the greatest increase in the average power was achieved by the case st8, with a value of 5.3218 × 10 5 W, which supposes an increase of 9.599% in comparison with the value obtained by the clean wind turbine. Moreover, the other cases with MTs mounted on st9 and st8st9 present a similar increase. At the NTM7.5 wind speed realization, the largest increment average power output was reached by the case with the MTs implemented in the station 8 of the blade, with an increase of the power in compared with the clean wind turbine of 4.425%. At the NTM10 wind speed realization the largest average power value is reached by the case with the MTs mounted on st8 with an increase with respect to the clean case of 3.282%. The largest increases for the three wind speed realization used in the present work have been performed by the case with the MTs in the blade station 8, corresponding to the case DU91W(2)250MT9520st8. Figure 14 illustrates in a bar plot the increments in the power output for every case in comparison with the clean case. The effect of mounting MTs on the blade stations studied in the present work is more significant at low wind speeds than at wind speed close to the nominal power. At those low wind speeds, the MTs can help to increase notably the wind turbine power output performance. Note that even thought the case DU91W(2)250MT9520st8st9 has more MTs mounted along the blade, its power performance in terms of average power output is quite similar in comparison with the other cases with MTs on stations 8 and 9. At the wind speed realization of NTM5, the greatest increase in the average power was achieved by the case st8, with a value of 5.3218 × 10 5 W, which supposes an increase of 9.599% in comparison with the value obtained by the clean wind turbine. Moreover, the other cases with MTs mounted on st9 and st8st9 present a similar increase. At the NTM7.5 wind speed realization, the largest increment average power output was reached by the case with the MTs implemented in the station 8 of the blade, with an increase of the power in compared with the clean wind turbine of 4.425%. At the NTM10 wind speed realization the largest average power value is reached by the case with the MTs mounted on st8 with an increase with respect to the clean case of 3.282%. The largest increases for the three wind speed realization used in the present work have been performed by the case with the MTs in the blade station 8, corresponding to the case DU91W(2)250MT9520st8. Figure  14 illustrates in a bar plot the increments in the power output for every case in comparison with the clean case. The effect of mounting MTs on the blade stations studied in the present work is more significant at low wind speeds than at wind speed close to the nominal power. At those low wind speeds, the MTs can help to increase notably the wind turbine power output performance. Note that even thought the case DU91W(2)250MT9520st8st9 has more MTs mounted along the blade, its power performance in terms of average power output is quite similar in comparison with the other cases with MTs on stations 8 and 9. After applying the BEM algorithm to all control volumes, the tangential and normal load distribution is known and global parameters such as thrust and bending moment at the root of the blade can be computed. In the current work, the mean values of both thrust and bending moment have been calculated for each wind speed realization. The thrust has been calculated by Equation (5) taking into account the trust distribution along the blade and the bending moment has been determined by Equation (6). After applying the BEM algorithm to all control volumes, the tangential and normal load distribution is known and global parameters such as thrust and bending moment at the root of the blade can be computed. In the current work, the mean values of both thrust and bending moment have been calculated for each wind speed realization. The thrust has been calculated by Equation (5) taking into account the trust distribution along the blade and the bending moment has been determined by Equation (6).
The Prandtl's tip loss correction factor F has been determined by the Equation (7), for both thrush and bending moment calculations: The variables used for thrust and bending moment estimations and the corresponding dimensions are shown in Table 5. All calculations are based on the 5 MW reference wind turbine described in [22].  Table 6 represents the mean values of thrust calculations for all MT cases according to the matrix described in Table 2. The thrust was calculated by integrating the thrust distribution along the blade by the Equation (5) for each wind speed realization. Afterwards, the average thrust was determined according to the wind realization duration of Figure 11. The increments in the average thrust of the microtabed blade cases have been illustrated in Table 6 in comparison with the clean case. The average thrust has been experienced an increase in every case with MTs implemented and once again, the case with the MTs mounted on the blade station 8 presents larger increments than the other cases, which is in concordance with the results of power output presented in Figure 14. The bending moment in the root of the blade has been determined by Equation (6) taking into account the bending moment distribution along the blade. The bending moment in the root of the blade was determined for each wind speed realization and computed along blade according to Equation (6). Afterwards, the mean bending moment was determined according to the wind realization duration of Figure 11. The increments in the average bending moment of the microtabed blade cases have been illustrated in Table 7 in comparison with the clean case. No extraordinary increment in the mean bending moment in the root of the blade has been found. All the increments in the bending moment due to the MTs implementation are acceptable taking into account the increase in the average power output production presented in Figure 14.

Conclusions
A parametric study for design and analysis of a MT on an airfoil has been carried out. To that end, 2D computational fluid dynamic simulations have been performed at Reynolds number of Re = 7 × 10 6 . The MT design attributes resulting from the simulations have allowed the sizing and positioning of the passive device based on aerodynamic performance. Comparisons of the CFD simulations and the DOWEC results have been made and verified the effectiveness of the MTs as flow control devices to increase the aerodynamics performance. The case DU91W(2)250MT9520 with the MT positioned at 95% of c and with the height of 2% is the one with the best aerodynamic performance in terms of lift-to-drag ratio. Afterwards, BEM based computations have been carried out to investigate the effects of that designed MT on the power performance of a 5 MW wind turbine. An increase on the average wind turbine power output has been found in the current study due to the implementation of MTs at different blade stations. That increase is more notable for the wind speed realizations with lower average wind speed NTM5. However, that increase is still significant with the wind speed realizations with average speeds of 7.5 m/s NTM7.5 and 10 m/s NTM10. In those cases, the increase in the power output is lower but still important. The best results in terms of average power are reached by the case denoted by DU91W(2)250MT9520st8 with the MTs implemented into the blade station 8. The largest increase in thrust has also been achieved by the case with the MTs mounted on the blade station 8. As expected, the increase in the wind turbine power output due to the MTs implementation leads to an augmentation in the bending moment in the root of the blade. However, this increase in the bending moment is acceptable taking into account the raise in the average power output production achieved by all cases. Moreover, no significant variation in the power increase has been found for the other MT locations st9 and st8st9. Because of the cheaper assembly of MT in only one blade station, the case with the MTs on the station 8 DU91W(2)250MT9520st8 is recommended.
The results of the current study shows that carefully analysis of MT height and location in the pressure surface from the airfoil leading edge, 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.