jBAY Modeling of Vane-Type Vortex Generators and Study on Airfoil Aerodynamic Performance

The increased demand for wind power is related to changes in the sizes of wind turbines and the development of flow control devices, such as vortex generators (VGs). In the present study, an analysis of the vortices generated by a vane-type VG is performed. To that end, the aerodynamic performance of a DU97W300 airfoil with and without VG is evaluated. The jBAY source term model was implemented for simulation of a triangular-shaped VG and the resolution of the fully meshed computational fluid dynamics (CFD) model. Reynolds-averaged Navier–Stokes (RANS) based simulations were used to calculate the effect of VGs in steady state, and the detached eddy simulation (DES) method was used for angles of attack (AoAs) around the stall situation. All jBAY based numerical simulations were carried out with a Reynolds number of Re = 2 × 106 to analyze the influence of VGs with AoAs between 0 and 20° and were validated versus experimental wind tunnel results. The results show that setting up a VG device on an airfoil benefits its aerodynamic performance and that the use of the jBAY model for simulation is accurate and efficient.


Introduction
The increasing need to obtain more energy from renewable sources and develop efficient machines is among the reasons for the growth in the wind power industry, becoming the generating technology with the highest rate of installations (see Aramendia et al. [1]). The use of wind power has increased by 9% for the last 15 years, as presented by Houghton et al. [2], which is a clear example of the importance of this technology in the overall power supply. Nowadays, in order to achieve greater power values, rotor blades of 60 m or longer are used, which can generate between 5 and 10 MW, according to Errasti et al. [3] and Baldacchino et al. [4]. This increase in power demand is related to an increase in power consumption, leading to bigger blades being continuously designed with the aim of generating more power, achieving blades larger than 80 m [4]. Changes in the size and complexity of wind turbines means that modifications in wind blade design must be made in order to maintain good efficiency and performance of the system (see Fernandez-Gamiz et al. [5] and Martínez-Filgueira et al. [6]).
Flow control needs to be considered when talking about large wind turbine blades. These big airfoils lead to flow separation, which becomes a bigger problem with larger rotors because of the larger available wind field (see Baldacchino et al. [4]). This premature and undesirable flow separation in the wings is generated by the adverse pressure gradient generated by the shape of the airfoil at some angles of attack. Thus, having thicker airfoils is undesirable due to the higher adverse pressure gradients and should be avoided as much as possible. Therefore, different flow control devices have been developed in order to delay the flow separation of the whole system with some angles of attack, as shown by Johnson et al. [7] and Timmer et al. [8].
Despite the fact that these flow control devices were initially designed for aeronautic purposes in the works of Taylor [9] and Bragg et al. [10], they have also been used in advanced technologies (see [11][12][13]). Classifications of these control devices were made by Johnson et al. [7] and Zaman et al. [14], dividing them into active or passive. The different control techniques are the reason for this division, with active controllers, which need auxiliary energy, and passive controllers, which improve the efficiency of the elements in which they are working without any independent energy consumption. Astolfi et al. [15] presented three cases to evaluate wind turbine power curve upgrades under different possible scenarios through operational data. A multivariate linear method for selecting the most appropriate input for modeling was proposed by Terzi et al. [16]. Recently, it was applied to a multi-megawatt wind turbine with different passive flow control items installed. These active and passive devices are used to improve the performance of the wings. Among the variety of flow control devices available, as presented by Johnson et al. [7], passive vortex generators (VGs) are the most frequently used, as indicated by Florentie et al. [17]. They have proved to be able to improve the power generated by a wind turbine nearly 24% compared to a turbine without them (see Øye [18] and Miller et al. [19]).
Vortex generators are passive flow control devices that consist of small items with several different shapes and are placed with a certain angle respect to the flow (see Johnson et al. [7]). Schubauer et al. [20] showed that these devices can be used to delay the flow separation caused by the adverse pressure gradients generated by the shape and angle of attack in each situation. With the aim to delay the flow separation as much as possible and completely comprehend the VG influence, theoretical and experimental research has been carried out (see Lin [21]). A larger lift coefficient was obtained and the stall occurred at a greater angle of attack, which improved the efficiency. Furthermore, according to Gao et al. [22] and Timmer et al. [8], the airfoil performance can be enhanced with delayed flow separation.
Betterton et al. [23] stated that two main configurations of vortex generators are used nowadays, counter-rotational and co-rotational. In the co-rotational configuration, all VGs have the same angle with respect to the flow, and in the counter-rotational layout, opposite angles are designed (see Figure 1). Taking into account the research done by Godard and Stanislas [24], Jirasek [25], and Bray [26], it was proven that the counter-rotational configuration is a more effective way to suppress flow separation for most applications, see Bur et al. [27].
devices have been developed in order to delay the flow separation of the whole system with some angles of attack, as shown by Johnson et al. [7] and Timmer et al. [8].
Despite the fact that these flow control devices were initially designed for aeronautic purposes in the works of Taylor [9] and Bragg et al. [10], they have also been used in advanced technologies (see [11][12][13]). Classifications of these control devices were made by Johnson et al. [7] and Zaman et al. [14], dividing them into active or passive. The different control techniques are the reason for this division, with active controllers, which need auxiliary energy, and passive controllers, which improve the efficiency of the elements in which they are working without any independent energy consumption. Astolfi et al. [15] presented three cases to evaluate wind turbine power curve upgrades under different possible scenarios through operational data. A multivariate linear method for selecting the most appropriate input for modeling was proposed by Terzi et al. [16]. Recently, it was applied to a multi-megawatt wind turbine with different passive flow control items installed. These active and passive devices are used to improve the performance of the wings. Among the variety of flow control devices available, as presented by Johnson et al. [7], passive vortex generators (VGs) are the most frequently used, as indicated by Florentie et al. [17]. They have proved to be able to improve the power generated by a wind turbine nearly 24% compared to a turbine without them (see Øye [18] and Miller et al. [19]).
Vortex generators are passive flow control devices that consist of small items with several different shapes and are placed with a certain angle respect to the flow (see Johnson et al. [7]). Schubauer et al. [20] showed that these devices can be used to delay the flow separation caused by the adverse pressure gradients generated by the shape and angle of attack in each situation. With the aim to delay the flow separation as much as possible and completely comprehend the VG influence, theoretical and experimental research has been carried out (see Lin [21]). A larger lift coefficient was obtained and the stall occurred at a greater angle of attack, which improved the efficiency. Furthermore, according to Gao et al. [22] and Timmer et al. [8], the airfoil performance can be enhanced with delayed flow separation.
Betterton et al. [23] stated that two main configurations of vortex generators are used nowadays, counter-rotational and co-rotational. In the co-rotational configuration, all VGs have the same angle with respect to the flow, and in the counter-rotational layout, opposite angles are designed (see Figure  1). Taking into account the research done by Godard and Stanislas [24], Jirasek [25], and Bray [26], it was proven that the counter-rotational configuration is a more effective way to suppress flow separation for most applications, see Bur et al. [27]. In Figure 1, is the free stream velocity of the oncoming flow, and h is the height and l the length of the VG in the chord direction. The distance between two devices in the co-rotating and counter-rotating VG configuration is represented by λ and L, respectively, defining the distance between two devices in the trailing edge in the counter-rotating configuration. In Figure 1, U ∞ is the free stream velocity of the oncoming flow, and h is the height and l the length of the VG in the chord direction. The distance between two devices in the co-rotating and counter-rotating VG configuration is represented by λ and L, respectively, defining the distance between two devices in the trailing edge in the counter-rotating configuration.
The shape of the vortex generator and the location on the airfoil have a great impact on the flow control. Regarding the shape, extensive studies have been carried out with a triangular-shaped VG in order to determine its boundary-layer separation control effectiveness (see Lin [21] and Gad-el-Hak et al. [28]). Moreover, Zhang et al. [29] studied the effects of semi-delta VG on the aerodynamic performance of airfoils of thick wind turbines. Hansen et al. [30] presented an engineering model to calculate the generated vortex strength of rectangular and triangular VGs. In order to obtain the best solution for each case, parametric studies are necessary (see Jirasek [25]). Usually, these studies are conducted by means of fully meshed numerical simulations, which are very computationally demanding due to the high number of cells required to mesh the VG (see Fernandez-Gamiz et al. [31]).
A new source term model was settled by Bender et al. [32], called the BAY model, based on the Joukowski lift theorem and the thin airfoil theory, which uses forces to obtain a solution. Vane-type vortex generators were used for simulation by the finite volume Navier-Stokes code, while ignoring the condition to totally define the geometry in the mesh. For calibration, an analysis was conducted and the results achieved by this source term model were compared with those obtained for a fully meshed VG, as done by Errasti et al. [3]. Eventually, the BAY model for different flows with vane-type vortex generators was successfully evaluated by Dudek [33], changing the landscape of simulation of this type of aerodynamic element, enabling much faster analysis. Recently, an updated version called the jBAY model was developed and published by Jirasek [25]. This new version was based on the previous BAY model of Bender et al. [32], for which the technique used for flow simulation was more suitable.
The main objective of the current study was to implement a triangular vane-type vortex generator on a DU97W300 based on the jBAY source term model. This airfoil is part of the reference 5 MW wind turbine developed by NREL presented in Jonkman et al. [34]. Figure 2a represents the airfoil outline and Figure 2b the position of the airfoil on the blade. The aerodynamic performance of the airfoil is defined in the current study by the lift-to-drag ratio, and it was evaluated with and without a VG. Additionally, vortex trajectory, vortex decay, wall shear stress parameter, pressure coefficient distribution along the airfoil surface, and vortex visualization downstream of the VG were studied. The main advantage of the jBAY to model a triangular VG compared with the fully meshed model is that it is easy to implement, since no geometry of the VG is required. Additionally, this method is very flexible in case the user wants to change the size or geometry. The results obtained from this study provide specific information about VG performance by using the jBAY source term model. The shape of the vortex generator and the location on the airfoil have a great impact on the flow control. Regarding the shape, extensive studies have been carried out with a triangular-shaped VG in order to determine its boundary-layer separation control effectiveness (see Lin [21] and Gad-el-Hak et al. [28]). Moreover, Zhang et al. [29] studied the effects of semi-delta VG on the aerodynamic performance of airfoils of thick wind turbines. Hansen et al. [30] presented an engineering model to calculate the generated vortex strength of rectangular and triangular VGs. In order to obtain the best solution for each case, parametric studies are necessary (see Jirasek [25]). Usually, these studies are conducted by means of fully meshed numerical simulations, which are very computationally demanding due to the high number of cells required to mesh the VG (see Fernandez-Gamiz et al. [31]).
A new source term model was settled by Bender et al. [32], called the BAY model, based on the Joukowski lift theorem and the thin airfoil theory, which uses forces to obtain a solution. Vane-type vortex generators were used for simulation by the finite volume Navier-Stokes code, while ignoring the condition to totally define the geometry in the mesh. For calibration, an analysis was conducted and the results achieved by this source term model were compared with those obtained for a fully meshed VG, as done by Errasti et al. [3]. Eventually, the BAY model for different flows with vanetype vortex generators was successfully evaluated by Dudek [33], changing the landscape of simulation of this type of aerodynamic element, enabling much faster analysis. Recently, an updated version called the jBAY model was developed and published by Jirasek [25]. This new version was based on the previous BAY model of Bender et al. [32], for which the technique used for flow simulation was more suitable.
The main objective of the current study was to implement a triangular vane-type vortex generator on a DU97W300 based on the jBAY source term model. This airfoil is part of the reference 5 MW wind turbine developed by NREL presented in Jonkman et al. [34]. Figure 2a represents the airfoil outline and Figure 2b the position of the airfoil on the blade. The aerodynamic performance of the airfoil is defined in the current study by the lift-to-drag ratio, and it was evaluated with and without a VG. Additionally, vortex trajectory, vortex decay, wall shear stress parameter, pressure coefficient distribution along the airfoil surface, and vortex visualization downstream of the VG were studied. The main advantage of the jBAY to model a triangular VG compared with the fully meshed model is that it is easy to implement, since no geometry of the VG is required. Additionally, this method is very flexible in case the user wants to change the size or geometry. The results obtained from this study provide specific information about VG performance by using the jBAY source term model. The paper is structured as follows: Section 1 provides background information about the application of flow control devices in wind turbines and, in particular, the implementation of vortex generators. Section 2 gives a detailed description of the numerical model defined by means of computational fluid dynamics (CFD), while the results obtained and the effects on the airfoil The paper is structured as follows: Section 1 provides background information about the application of flow control devices in wind turbines and, in particular, the implementation of vortex generators. Section 2 gives a detailed description of the numerical model defined by means of computational fluid dynamics (CFD), while the results obtained and the effects on the airfoil aerodynamic performance are discussed in Section 3. Finally, Section 4 presents the more noteworthy conclusions obtained in the present study.

jBAY Model
Computational fluid dynamics (CFD) is the most common method used to simulate airflow and estimate airfoil efficiency. The combination of this tool and experiments in wind tunnels can be considered as a reliable method to conduct any parametric flow study. However, the fully meshed numerical simulations of VGs are very computationally demanding to carry out due to the high number of cells necessary to capture the VG effects.
In the current work, the jBAY source term model presented by Jirasek [25] and based on the BAY model developed by Bender et al. [32] was implemented in the commercial CFD code STAR CCM+ ® [35]. The BAY model was modeled to simulate vane-type vortex generators in finite volume Navier-Stokes codes and allows substituting the VG geometry by a subdomain of similar geometry at the same location. Then, a body force distribution normal to the local flow is applied to simulate the effect of the vane (see Figure 3). The normal force represented by Equation (1) is introduced as → L in Figure 3b. This normal force simulates the force generated by the VG without having to design and mesh it (see Fernandez-Gamiz et al. [36]). The force for each cell in the VG domain is calculated with Equation (1): where → L cell is the force of the corresponding single cell; C VG is a relaxation factor, which, according to Jirasek [25], usually has a value near 10; ρ is the local density; t (normal and tangential vectors to the VG); S VG is the area of VG; V cell is the volume of a single cell; and V S is the total volume of grid cells where the model is applied (see Errasti et al. [3]).
Energies 2020, 13, x FOR PEER REVIEW 7 of 15 aerodynamic performance are discussed in Section 3. Finally, Section 4 presents the more noteworthy conclusions obtained in the present study.

jBAY Model
Computational fluid dynamics (CFD) is the most common method used to simulate airflow and estimate airfoil efficiency. The combination of this tool and experiments in wind tunnels can be considered as a reliable method to conduct any parametric flow study. However, the fully meshed numerical simulations of VGs are very computationally demanding to carry out due to the high number of cells necessary to capture the VG effects.
In the current work, the jBAY source term model presented by Jirasek [25] and based on the BAY model developed by Bender et al. [32] was implemented in the commercial CFD code STAR CCM+ ® [35]. The BAY model was modeled to simulate vane-type vortex generators in finite volume Navier-Stokes codes and allows substituting the VG geometry by a subdomain of similar geometry at the same location. Then, a body force distribution normal to the local flow is applied to simulate the effect of the vane (see Figure 3). The normal force represented by Equation (1) is introduced as in Figure  3b. This normal force simulates the force generated by the VG without having to design and mesh it (see Fernandez-Gamiz et al. [36]). The force for each cell in the VG domain is calculated with Equation (1): where L is the force of the corresponding single cell; C is a relaxation factor, which, according to Jirasek [25], usually has a value near 10; ρ is the local density; u is the local velocity; b is a unit factor defined as b = n x t (normal and tangential vectors to the VG); S is the area of VG; V is the volume of a single cell; and V is the total volume of grid cells where the model is applied (see Errasti et al. [3]).  For this work, the forces, instead of applying them in all the subdomain cells as specified by BAY model, have been applied in those cells that represent the VG, as shown in Figure 4. aerodynamic performance are discussed in Section 3. Finally, Section 4 presents the more noteworthy conclusions obtained in the present study.

jBAY Model
Computational fluid dynamics (CFD) is the most common method used to simulate airflow and estimate airfoil efficiency. The combination of this tool and experiments in wind tunnels can be considered as a reliable method to conduct any parametric flow study. However, the fully meshed numerical simulations of VGs are very computationally demanding to carry out due to the high number of cells necessary to capture the VG effects.
In the current work, the jBAY source term model presented by Jirasek [25] and based on the BAY model developed by Bender et al. [32] was implemented in the commercial CFD code STAR CCM+ ® [35]. The BAY model was modeled to simulate vane-type vortex generators in finite volume Navier-Stokes codes and allows substituting the VG geometry by a subdomain of similar geometry at the same location. Then, a body force distribution normal to the local flow is applied to simulate the effect of the vane (see Figure 3). The normal force represented by Equation (1) is introduced as in Figure  3b. This normal force simulates the force generated by the VG without having to design and mesh it (see Fernandez-Gamiz et al. [36]). The force for each cell in the VG domain is calculated with Equation (1): where L is the force of the corresponding single cell; C is a relaxation factor, which, according to Jirasek [25], usually has a value near 10; ρ is the local density; u is the local velocity; b is a unit factor defined as b = n x t (normal and tangential vectors to the VG); S is the area of VG; V is the volume of a single cell; and V is the total volume of grid cells where the model is applied (see Errasti et al. [3]).

Computational Domain and VG Setup
The computational model is a DU97W300 with 0.65 m chord length (c) and 0.0175 m span length. A single VG was set up instead of a pair, with the assumption that the flow created behind the VG is symmetric. This strategy could contribute to reducing the computational time.
An O-shape computational domain was defined with a radius of 30 times the airfoil chord length. A symmetry boundary condition was implemented for the front and back sides and a nonslip wall boundary condition for the airfoil. A free stream condition was determined for the far-field domain. The domain was discretized by a structured mesh with a total of 6.64 million cells. Mesh refinement was defined in the region close to the leading and trailing edge, as well as at the VG position area, with a mesh growth rate of 1.1.
The configuration of the main VG variables in the airfoil is described in Figure 5. Table 1 shows the values of geometric parameters employed in this case. The VG was positioned at 30% of the chord length from the leading edge of the airfoil, x c = 0.3.

Computational Domain and VG Setup
The computational model is a DU97W300 with 0.65 m chord length (c) and 0.0175 m span length. A single VG was set up instead of a pair, with the assumption that the flow created behind the VG is symmetric. This strategy could contribute to reducing the computational time.
An O-shape computational domain was defined with a radius of 30 times the airfoil chord length. A symmetry boundary condition was implemented for the front and back sides and a nonslip wall boundary condition for the airfoil. A free stream condition was determined for the far-field domain. The domain was discretized by a structured mesh with a total of 6.64 million cells. Mesh refinement was defined in the region close to the leading and trailing edge, as well as at the VG position area, with a mesh growth rate of 1.1.
The configuration of the main VG variables in the airfoil is described in Figure 5. Table 1 shows the values of geometric parameters employed in this case. The VG was positioned at 30% of the chord length from the leading edge of the airfoil, x c ⁄ = 0.3.

Numerical Models
Two numerical models were set up for the numerical simulations in the current work. In an incompressible and steady flow, the k-ω shear stress transport (SST) model of Menter [37] was used with the Reynolds-averaged Navier-Stokes (RANS) method to get a solution for cases where the angle of attack (AoA) is lower than 12°. However, according to Zhong et al. [38], the simplification of turbulence in the RANS method increases the inaccuracy of CFD, especially for phenomena directly related to turbulence. For this reason, at higher AoAs where the clean airfoil would be in a stall situation, Detached eddy simulation (DES) was set up. For instance, Meana-Fernandez et al. [39] determined that the results obtained in a DU-06-W-200 airfoil by the use of DES predicted the aerodynamic forces with good accuracy and provided more realistic flow patterns. This model combines features of the standard Spalart-Allmaras RANS model in boundary layers and large eddy simulation in unsteady separated regions (see Guo et al. [40]). The solution was considered to be converged with a three-order-of-magnitude drop in the numerical residuals.
One of the most important variables to define in a DES simulation is the Courant number (C). It depends on cell size (∆x), defined time steps (∆t), and flow velocity (u) as shown in Equation (2). Briefly, this number indicates how a particle travels from one cell to an adjacent cell in a time step. If the Courant number is lower than or equal to 1, a flow particle will move from one cell to another in

Numerical Models
Two numerical models were set up for the numerical simulations in the current work. In an incompressible and steady flow, the k-ω shear stress transport (SST) model of Menter [37] was used with the Reynolds-averaged Navier-Stokes (RANS) method to get a solution for cases where the angle of attack (AoA) is lower than 12 • . However, according to Zhong et al. [38], the simplification of turbulence in the RANS method increases the inaccuracy of CFD, especially for phenomena directly related to turbulence. For this reason, at higher AoAs where the clean airfoil would be in a stall situation, Detached eddy simulation (DES) was set up. For instance, Meana-Fernandez et al. [39] determined that the results obtained in a DU-06-W-200 airfoil by the use of DES predicted the aerodynamic forces with good accuracy and provided more realistic flow patterns. This model combines features of the standard Spalart-Allmaras RANS model in boundary layers and large eddy simulation in unsteady separated regions (see Guo et al. [40]). The solution was considered to be converged with a three-order-of-magnitude drop in the numerical residuals.
One of the most important variables to define in a DES simulation is the Courant number (C). It depends on cell size (∆x), defined time steps (∆t), and flow velocity (u) as shown in Equation (2). Briefly, this number indicates how a particle travels from one cell to an adjacent cell in a time step. If the Courant number is lower than or equal to 1, a flow particle will move from one cell to another in one time step; if the Courant number is bigger than 1, particles will move two or more cells in one time step. A Courant value of 20 was used in all unsteady cases studied, with a time step of 10 −5 . C = u· ∆t ∆x (2) In this study, the Reynolds number is Re = 2 × 10 6 , based on airfoil chord length c = 0.65 m; ν = 1.51 × 10 −5 is the kinematic viscosity; and the free stream velocity is U ∞ = 46.52 m/s (see Equation (3)):

Results
In the current study, the effects of VG implementation on an aerodynamic airfoil was investigated, with the following aerodynamic properties: lift and drag coefficients, vortex trajectory, vortex decay and wall shear stress, pressure coefficient distribution along the airfoil surface, and vortex visualization downstream of the VG for different angles of attack.

Lift and Drag Coefficients
According to Sørensen et al. [41], Fernandez-Gamiz et al. [42], and Timmer et al. [8], placing the VG over an airfoil results in a higher lift coefficient and a slight increase of residual drag. In all cases, this device allows the airfoil to be in steady state for larger AoAs than without it, getting higher lift values above all at stall angles of the clean airfoil. The lift and drag coefficients of the airfoil can be calculated with Equation (4): where C L and C D are the lift and drag coefficients, respectively; f pressure f and f shear f are the force vectors on the f surface face; ρ re f , V re f , and a re f are the reference density, velocity, and airfoil area, respectively; and n D is the specified direction vector depending on each force. The pressure force vector on a surface is computed by means of Equation (5): where p f is the face static pressure; p re f is the reference pressure; and a f is the face area vector. The shear force on surface face f can be computed by Equation (6): where T f is the stress tensor. The analyzed angles of attack (AoAs) are in the range between 0 and 20 • . The lift coefficients increase linearly until stall angles. These stall angles are located around 12 • in the case of clean airfoil, according to Timmer et al. [8], and around 18-20 • in the case of airfoil with VG according to the study by Gao et al. [22]. At those angles, effects related to flow separation appear and it begins to be affected by a worse lift coefficient, as shown in Figure 6a. The higher lift value appears just before the stall point, and it is greater in both value and AoA when using VGs. With the aim of validating the numerical results, they were compared with those of Gao et al. [22] in lift and drag coefficients and with the experimental data in the clean airfoil case [8]. Some differences are noted in the clean case between the numerical and experimental data at the beginning of stalling angles. Usually, CFD overestimates the lift coefficient at high AoAs (see the study by Li et al. [43]). Moreover, Raffel et al. [44] found some differences in light-stall situation between CFD and experimental data. Figure 6b shows the exponential increase of drag coefficient that arises beyond the stall angle.

Vortex Trajectory and Decay
The study performed in this work on peak vorticity in the VG was divided into two parts, investigating the values obtained for AoA in each case. The first part was focused on analyzing the path following peak vorticity values, whether horizontally or vertically. The second part studied the decay of peak vorticity values due to their distance in each moment from the VG device.
As concluded by Zhen et al. [45], vortices can be assumed as circular and it is in their center where the maximum peak value appears. Thanks to this peak localization, it is possible to determine the vortex path of each AoA downstream of the VG. The path of peak values is an important parameter to reach a more complete characterization of the created vortex and optimize the VG device. Even though several simulations and studies were done around this characterization and the detection of these peak values proved to be critical, Fernandez-Gamiz et al. [42] addressed successfully the characterization of the vortex path.
Meanwhile, in the case of a flat plane with VG devices, the vortex vertical path does not change its trajectory and is usually horizontal (see Ibarra-Udaeta et al. [46]). In the case of an airfoil profile, it was concluded that the vortex tends to follow parallel to the wall. This behavior of the vortex was concluded in Errasti et al. [3], where a VG device was set up on a flat plate and before going down a

Vortex Trajectory and Decay
The study performed in this work on peak vorticity in the VG was divided into two parts, investigating the values obtained for AoA in each case. The first part was focused on analyzing the path following peak vorticity values, whether horizontally or vertically. The second part studied the decay of peak vorticity values due to their distance in each moment from the VG device.
As concluded by Zhen et al. [45], vortices can be assumed as circular and it is in their center where the maximum peak value appears. Thanks to this peak localization, it is possible to determine the vortex path of each AoA downstream of the VG. The path of peak values is an important parameter to reach a more complete characterization of the created vortex and optimize the VG device. Even though several simulations and studies were done around this characterization and the detection of these peak values proved to be critical, Fernandez-Gamiz et al. [42] addressed successfully the characterization of the vortex path.
Meanwhile, in the case of a flat plane with VG devices, the vortex vertical path does not change its trajectory and is usually horizontal (see Ibarra-Udaeta et al. [46]). In the case of an airfoil profile, it was concluded that the vortex tends to follow parallel to the wall. This behavior of the vortex was concluded in Errasti et al. [3], where a VG device was set up on a flat plate and before going down a linear ramp. In this study, the vortex tended to continue the bottom wall path going down the ramp. The simulations performed for this purpose in the current work show that the generated vortex axes had similar vertical behavior in the first 20 mm after VG at all AoAs. From that distance onwards, the vortex vertical path slightly evolves in different ways, decaying faster at higher AoAs due to the length of the vortex, which is shorter at higher AoAs. Near stall state, the generated vortex in the boundary layer tends to separate from the surface of the airfoil, spreading in turbulent wakes that introduce major fluctuations in forces and increase the deviations of drag and lift coefficients. Thanks to VG devices, the boundary layer gets attached to the wall for larger degrees, but gets worse attitude at higher AoAs until a stall point is detected at 20 • . The vertical trajectories can be seen in Figure 7, where the black line represents the airfoil surface.
Energies 2020, 13, x FOR PEER REVIEW 11 of 15 linear ramp. In this study, the vortex tended to continue the bottom wall path going down the ramp. The simulations performed for this purpose in the current work show that the generated vortex axes had similar vertical behavior in the first 20 mm after VG at all AoAs. From that distance onwards, the vortex vertical path slightly evolves in different ways, decaying faster at higher AoAs due to the length of the vortex, which is shorter at higher AoAs. Near stall state, the generated vortex in the boundary layer tends to separate from the surface of the airfoil, spreading in turbulent wakes that introduce major fluctuations in forces and increase the deviations of drag and lift coefficients. Thanks to VG devices, the boundary layer gets attached to the wall for larger degrees, but gets worse attitude at higher AoAs until a stall point is detected at 20°. The vertical trajectories can be seen in Figure 7, where the black line represents the airfoil surface. As concluded in Ibarra-Udaeta et al. [46], the vortex peak path tends to deviate laterally in the direction marked by the device. In Figure 8 the lateral displacement can be noticed, which has a similar deviation for all AoAs in steady angles and stall angles, but is less displaced at high degrees. The trend of its deviation is higher at first distances in all cases, with a decay of deviation at 80 mm from the device. At 160 mm onwards, it is possible to see how the vortices at the highest AoA tend to go back to the device zero coordinate. Vortex decay is a well-studied parameter that consists of the strength force loss of a vortex. In a study by Martínez-Filgueira et al. [6], the peak decay was analyzed over a flat plate. In that study, the decay development was exponential; the further it was from the device, the smaller the peak. In As concluded in Ibarra-Udaeta et al. [46], the vortex peak path tends to deviate laterally in the direction marked by the device. In Figure 8 the lateral displacement can be noticed, which has a similar deviation for all AoAs in steady angles and stall angles, but is less displaced at high degrees. The trend of its deviation is higher at first distances in all cases, with a decay of deviation at 80 mm from the device. At 160 mm onwards, it is possible to see how the vortices at the highest AoA tend to go back to the device zero coordinate.
Energies 2020, 13, x FOR PEER REVIEW 11 of 15 linear ramp. In this study, the vortex tended to continue the bottom wall path going down the ramp. The simulations performed for this purpose in the current work show that the generated vortex axes had similar vertical behavior in the first 20 mm after VG at all AoAs. From that distance onwards, the vortex vertical path slightly evolves in different ways, decaying faster at higher AoAs due to the length of the vortex, which is shorter at higher AoAs. Near stall state, the generated vortex in the boundary layer tends to separate from the surface of the airfoil, spreading in turbulent wakes that introduce major fluctuations in forces and increase the deviations of drag and lift coefficients. Thanks to VG devices, the boundary layer gets attached to the wall for larger degrees, but gets worse attitude at higher AoAs until a stall point is detected at 20°. The vertical trajectories can be seen in Figure 7, where the black line represents the airfoil surface. As concluded in Ibarra-Udaeta et al. [46], the vortex peak path tends to deviate laterally in the direction marked by the device. In Figure 8 the lateral displacement can be noticed, which has a similar deviation for all AoAs in steady angles and stall angles, but is less displaced at high degrees. The trend of its deviation is higher at first distances in all cases, with a decay of deviation at 80 mm from the device. At 160 mm onwards, it is possible to see how the vortices at the highest AoA tend to go back to the device zero coordinate. Vortex decay is a well-studied parameter that consists of the strength force loss of a vortex. In a study by Martínez-Filgueira et al. [6], the peak decay was analyzed over a flat plate. In that study, the decay development was exponential; the further it was from the device, the smaller the peak. In Vortex decay is a well-studied parameter that consists of the strength force loss of a vortex. In a study by Martínez-Filgueira et al. [6], the peak decay was analyzed over a flat plate. In that study, the decay development was exponential; the further it was from the device, the smaller the peak. In this study, values of the peak vorticity for every AoA in relation to the range toward VG are shown in Figure 9. It can be observed that the value varies for each angle, with different behavior between steady and stall angles in the near wake behind the VG. Higher peak vorticity values correspond to steady degrees, decreasing at higher AoAs. After the first 40 mm, the value of vorticity starts converging, with a similar value until the last 175 mm from the device, where it is less than 5000 s -1 , and always lower for higher AoA.
this study, values of the peak vorticity for every AoA in relation to the range toward VG are shown in Figure 9. It can be observed that the value varies for each angle, with different behavior between steady and stall angles in the near wake behind the VG. Higher peak vorticity values correspond to steady degrees, decreasing at higher AoAs. After the first 40 mm, the value of vorticity starts converging, with a similar value until the last 175 mm from the device, where it is less than 5000 s -1 , and always lower for higher AoA.  Figure 10 shows the wall shear stress distribution along the first distances from the VG. The measurements were taken in a chordwise parallel line passing from the device trailing edge to 175 mm downstream. The results obtained agree with those obtained by Godard and Stanislas [24]. As in that case, the implementation of VGs increased the wall shear stress values downward in comparison with the clean airfoil. In the cases with a VG set up over the airfoil, all plots show an increasing trend up to 54 mm downstream from the leading edge of the VG device. From this point onward, the values decay gradually up to a value between 42% and 70% at the last 173 mm from the VG (at higher AoAs, the difference between maximum and minimum is smaller). A maximum wall shear stress value of 14.86 s is observed with an AoA of 15.25°.  Figure 10 shows the wall shear stress distribution along the first distances from the VG. The measurements were taken in a chordwise parallel line passing from the device trailing edge to 175 mm downstream. The results obtained agree with those obtained by Godard and Stanislas [24]. As in that case, the implementation of VGs increased the wall shear stress values downward in comparison with the clean airfoil.

Wall Shear Stress
Energies 2020, 13, x FOR PEER REVIEW 12 of 15 this study, values of the peak vorticity for every AoA in relation to the range toward VG are shown in Figure 9. It can be observed that the value varies for each angle, with different behavior between steady and stall angles in the near wake behind the VG. Higher peak vorticity values correspond to steady degrees, decreasing at higher AoAs. After the first 40 mm, the value of vorticity starts converging, with a similar value until the last 175 mm from the device, where it is less than 5000 s -1 , and always lower for higher AoA.  Figure 10 shows the wall shear stress distribution along the first distances from the VG. The measurements were taken in a chordwise parallel line passing from the device trailing edge to 175 mm downstream. The results obtained agree with those obtained by Godard and Stanislas [24]. As in that case, the implementation of VGs increased the wall shear stress values downward in comparison with the clean airfoil. In the cases with a VG set up over the airfoil, all plots show an increasing trend up to 54 mm downstream from the leading edge of the VG device. From this point onward, the values decay gradually up to a value between 42% and 70% at the last 173 mm from the VG (at higher AoAs, the difference between maximum and minimum is smaller). A maximum wall shear stress value of 14.86 s is observed with an AoA of 15.25°. In the cases with a VG set up over the airfoil, all plots show an increasing trend up to 54 mm downstream from the leading edge of the VG device. From this point onward, the values decay gradually up to a value between 42% and 70% at the last 173 mm from the VG (at higher AoAs, the difference between maximum and minimum is smaller). A maximum wall shear stress value of 14.86 s −1 is observed with an AoA of 15.25 • .

Wall Shear Stress
In the case of the clean airfoil, two simulation types were evaluated: steady simulations at AoAs lower than 12.25 • and DES based simulations at higher degrees. On the one hand, steady clean airfoil simulations present a maximum value closer to the device position than flow-controlled simulations, at 37 mm. From this point onward, the wall shear stress value starts decaying gradually. All steady simulations performed with the clean airfoil show a similar decaying gradient, with higher values at lower AoAs. On the other hand, the stall simulations present a maximum value at the first sampling point 20 mm from the VG leading edge and a subsequent decay until values reach next to zero. From this point onward, the wall shear stress values present a small increase that tends to be more stable. At higher AoAs, the distance of the null wall shear stress value from the VG device is closer. The difference between values at the pressure side and suction side represents the capacity to generate lift. The VG implementation barely affects the Cp distribution at low AoAs between the clean and flow-controlled airfoil. However, the increment in the pressure coefficient due to VG implementation is very noticeable at high angles of attack. It is from 12 • onward where the device demonstrates high performance capacity. It is noteworthy that there exists a cut value in the suction side case due to the presence of the VG.

Pressure Coefficient Distribution
Energies 2020, 13, x FOR PEER REVIEW 13 of 15 In the case of the clean airfoil, two simulation types were evaluated: steady simulations at AoAs lower than 12.25° and DES based simulations at higher degrees. On the one hand, steady clean airfoil simulations present a maximum value closer to the device position than flow-controlled simulations, at 37 mm. From this point onward, the wall shear stress value starts decaying gradually. All steady simulations performed with the clean airfoil show a similar decaying gradient, with higher values at lower AoAs. On the other hand, the stall simulations present a maximum value at the first sampling point 20 mm from the VG leading edge and a subsequent decay until values reach next to zero. From this point onward, the wall shear stress values present a small increase that tends to be more stable. At higher AoAs, the distance of the null wall shear stress value from the VG device is closer. Figure 11 shows the pressure distributions on the clean DU97W300 airfoil and the airfoil with VG at 0°, 4°, 6°, 10.27°, 15.25°, and 18.18° AoA. The airfoil profile and VG are illustrated with black lines. The difference between values at the pressure side and suction side represents the capacity to generate lift. The VG implementation barely affects the Cp distribution at low AoAs between the clean and flow-controlled airfoil. However, the increment in the pressure coefficient due to VG implementation is very noticeable at high angles of attack. It is from 12° onward where the device demonstrates high performance capacity. It is noteworthy that there exists a cut value in the suction side case due to the presence of the VG.   Figure 12 illustrates the axial velocity fields for a qualitative comparison between six angles, AoA = 0 • , 4 • , 6 • , 10.37 • , 15.25 • , and 18.18 • . These snapshots were taken at the plane position five device heights downstream of the VG. According to Urkiola et al. [47], at that position the vortex is totally developed, and this shows how the vortex shape varies with the angle of attack. As a qualitative observation, it denotes that the size of the vortex increases with the AoA. The effect of the VG device is lower at small AoAs in comparison with higher angles. This is in agreement with the results presented in Figure 11 about the pressure distribution along the airfoil.

Vortex Vizualization
Energies 2020, 13, x FOR PEER REVIEW 14 of 15 Figure 12 illustrates the axial velocity fields for a qualitative comparison between six angles, AoA = 0°, 4°, 6°, 10.37°, 15.25°, and 18.18°. These snapshots were taken at the plane position five device heights downstream of the VG. According to Urkiola et al. [47], at that position the vortex is totally developed, and this shows how the vortex shape varies with the angle of attack. As a qualitative observation, it denotes that the size of the vortex increases with the AoA. The effect of the VG device is lower at small AoAs in comparison with higher angles. This is in agreement with the results presented in Figure 11 about the pressure distribution along the airfoil.

Conclusions
The aerodynamic efficiency and vortex generated downstream of a delta VG device were investigated in this study. The geometry of the VG device is a semi-delta shape with a β = 18° incidence inclination, 5 mm height, and 17 mm length, and it was modelled by the so-called jBAY source term model. The distance between two devices is 10 mm measured from the leading edge and 20 mm from the trailing edge. The device was placed at ⁄ = 0.30 in a DU97W300 airfoil of 0.65 m of chord. Computational simulations for this purpose were performed with a Reynolds number Re = 2 × 10 and free stream velocity U = 46.52 m/s using the DES method. The main conclusions are as follows: First, setting up a VG device is an efficient way to improve the lift coefficient of the DU97W300 airfoil in spite of the residual drag increment due to shear forces on the suction side. The jBAY model could be a key tool for parametric studies to find the optimal sizes and configurations of vane-type VGs. In addition, the vertical vortex trajectory is conditioned by the wall. The vortex will tend to adhere to the airfoil suction side for more angles of attack (AoAs), delaying the stall situation

Conclusions
The aerodynamic efficiency and vortex generated downstream of a delta VG device were investigated in this study. The geometry of the VG device is a semi-delta shape with a β = 18 • incidence inclination, 5 mm height, and 17 mm length, and it was modelled by the so-called jBAY source term model. The distance between two devices is 10 mm measured from the leading edge and 20 mm from the trailing edge. The device was placed at x c = 0.30 in a DU97W300 airfoil of 0.65 m of chord. Computational simulations for this purpose were performed with a Reynolds number Re = 2 × 10 6 and free stream velocity U ∞ = 46.52 m/s using the DES method. The main conclusions are as follows: First, setting up a VG device is an efficient way to improve the lift coefficient of the DU97W300 airfoil in spite of the residual drag increment due to shear forces on the suction side. The jBAY model could be a key tool for parametric studies to find the optimal sizes and configurations of vane-type VGs.
In addition, the vertical vortex trajectory is conditioned by the wall. The vortex will tend to adhere to the airfoil suction side for more angles of attack (AoAs), delaying the stall situation and improving the lift coefficient at those degrees. It allows improvement of the aerodynamic performance of the airfoil. Furthermore, from the results obtained it can be concluded that the vortex peak decay has a similar downing gradient at any AoA, but with higher values at lower AoAs. Another remarkable fact is the difference between the first values of peak vorticity, which have a clear difference between AoA = 12.45 • and AoA = 15.25 • . It is in this degree range where the clean airfoil starts entering the stall situation. Furthermore, on the one hand, it is clear that the VG device increases the wall shear stress, which leads to a slight increment of the overall drag force at those AoAs where the clean airfoil would be in a stall regime, due to the circulation generated by the VG. On the other hand, implementing the VG has a noticeable effect on the increment of the pressure coefficient of the DU97W300 airfoil from an AoA of 12 • onward. At low AoAs, no clear difference exists between the Cp distribution on the clean and flow-controlled airfoils. Finally, it can be observed that the size of the vortex generated by implementing the VG increases with increasing angles of attack. Acknowledgments: The authors are thankful for the technical and human support provided by SGIker of UPV/EHU.

Conflicts of Interest:
The authors declare no conflict of interest.