Computational Characterization of a Rectangular Vortex Generator on a Flat Plate for Different Vane Heights and Angles

Vortex generators (VG) are passive flow control devices used for avoiding or delaying the separation of the boundary layer by bringing momentum from the higher layers of the fluid towards the surface. The Vortex generator usually has the same height as the local boundary layer thickness, and these Vortex generators can produce overload drag in some cases. The aim of the present study was to analyze the characteristics and path of the primary vortex produced by a single rectangular vortex generator on a flat plate for the incident angles of β = 10◦, 15◦, 18◦ and 20◦. A parametric study of the induced vortex was performed for six VG heights using Reynolds average Navier–Stokes equations at Reynodls number Re = 27,000 based on the local boundary layer thickness, using computational fluid dynamics techniques with OpenFOAM open-source code. In order to determine the vortex size, the so-called half-life radius was computed and compared with experimental data. The results showed a similar trend for all the studied vortex generator heights and incident angles with small variations for the vertical and the lateral paths. Additionally, 0.4H and 0.6H VG heights at incident angles of β = 18◦ and β = 20◦ showed the best performance in terms of vortex strength and generation of wall shear stress.


Introduction
In the last few decades, several flow control devices have been developed, and some of them were installed onto wings in order to improve aerodynamic efficiency. At present, this type of elements is being optimized by researchers to be introduced in wind turbines. A four-layer scheme was developed by Wood [1] for the classification of the different types of flow control devices. Wood classified the control devices as active or passive, geometric or fluidic, steady or unsteady, and attached or separated devices.
Flow control devices can be categorized as active or passive; this classification depends on their operating principle, see Aramendia et al. [2,3] and Barlas et al. [4]. Passive control techniques would represent an improvement in turbine efficiency and load reduction without external energy consumption. In contrast to passive devices as VGs, more complex algorithms are needed to get efficiency advantages by power output between 3% to 9% due to another type of passive flow control devices, such as Gurney flaps and VG rows implementation on the blades.
Some interesting advantages are shown by VGs, such as their small size which allows distributing a big number of them along the profile. They can also be installed at any time after the construction of the blade or replaced without any substantial cost.
The principal aim of this study was to characterize the generated primary vortex by means of computational techniques and, in this way, obtain the VG size and orientation fitting for their application in blades and airfoils. In order to develop a parametric study that covers a big range of cases, numerical simulations were carried out evaluating four different vane incident angles, β = 10 • , β = 15 • , β = 18 • , and β = 20 • and six different VG heights of 0.2H, 0.4H, 0.6H, 0.8H, H, and 1.2H. The CFD simulations in this study were carried out with OpenFOAM , see Weller et al. [29], open-source code; the vortex size was computed for all these cases. Section 2 explains the numerical setup; Section 3 shows the results and their discussion divided in four subsections, the primary vortex center path, vortex size, the wall shear stress, and the vortex strength; and finally, in Section 4, conclusions are presented.

Numerical Setup
CFD simulations were carried out in the present study to investigate and understand the development of the primary vortex generated by VGs. The vortex generated by a single vortex generator (VG) was simulated by the OpenFOAM open-source code. OpenFOAM is a library written in C++ that allows the manipulation and development of numerical solvers for problems of continuum mechanics. One of the main advantages of this noncommercial CFD code is that it can be improved and tailored to reproduce different physical phenomena depending on user interests.
A semi-implicit method for pressure-linked equations (SIMPLE) algorithm wasn applied for turbulent, steady-state, and incompressible flow solving the Reynolds average Navier-Stokes (RANS) equations. With the purpose of speeding up the convergence process, a potentialFoam solver was used to initialize the fields. This solver is appropriate to set initial conditions for the next step when a more advanced solver is applied. After that, a simpleFoam solver was used for solving RANS equations, and the turbulence model used in the present study was the k-ω shear-stress transport (SST) turbulence model presented by Menter [30]. Two transport equations were solved by this turbulence model, the one for the turbulent kinematic energy k and the specific dissipation ω. The used computational technique was justified on the steady Reynolds average Navier-Stokes (RANS) simulations made by Allan et al. [31] for a single VG in a flat plate for two incident angles computed with OVERFLOW, a code developed by the National Aeronautics and Space Administration (NASA) . These simulations were compared to experimentally obtained data, resulting on a good agreement for CFD results except for short distances downstream of the VG. In that study, the streamwise peak vorticity and vortex trajectory were better predicted by the shear-stress transport (SST) compared to the Spalart-Allmaras (SA) turbulence model. Second order discretization schemes were employed regarding the numerical schemes, following similar studies developed by Urkiola et al. [17].
The computational setup consists of a flat plate with a single VG whose leading edge (LE) is placed at 10.25 m from the inlet boundary, see Figure 1. This passive device consists of a rectangular vane with a length twice its height (H), which varies to make possible the parametric study. The normalization of the domain is done with respect to the conventional VG height (H); this height is equal to the boundary layer thickness δ = 0.25 m at the leading edge (LE) of the VG. Different boundary conditions have been set for the domain limits. According to Gutierrez-Amo et al. [23], the ceiling and the two lateral sides were set as slip conditions. Additionally, nonslip conditions were set for the vane and flat plate surfaces, while the upstream part of the domain was specified as inflow conditions with prescribed velocity according to the undisturbed velocity. The downstream part was set as an outlet assuming a fully developed flow.
The strategy carried out to save computational cost was to divide the mesh in blocks because the cell density in the mesh is different in order to capture more accurate results in those places where the vortex is generated, around and downstream of the VG. Each block has a particular stretching factor in the three spatial coordinates, as shown in Figure 2.  Two types of structured meshes in the computational domain were designed according to the different VG heights used for the parametric study. A first mesh with around 12 million hexaedral cells and 64 blocks for the 0.2H and 0.4H VG height cases, and a second mesh sized 11 million hexaedral cells and 44 blocks for the other three VG height cases were used with a normalized height of the closest cell to the wall of ∆y δ = 1.5 × 10 −6 . Both meshes have a block-based structure and are divided in regions which allow a different grid resolution. The mesh generation was performed by an OpenFoam block-mesh utility that reads a dictionary where points, faces, cells, and boundary are defined. A similar meshing strategy as that in Urkiola et al. [17] was followed in order to achieve a high quality grid.
The computation domain was normalized with respect to the conventional VG height H, which it is equal to the local boundary layer thickness at the leading edge of the VG. Thus, its width is 64 times, its height 10 times, and its length 75 times the local BL thickness (Figure 1).
In the present study, four incident angles were analyzed, β = 10 • , 15 • , 18 • , and 20 • Six VG heights were considered, 0.2H, 0.4H, 0.6H, 0.8H, H and 1.2H; for every case, the VG length is 2H. The free stream velocity of the fluid is U ∞ = 20 m · s −1 , the kinematic viscosity has a value of ν = 1.85185 × 10 −04 m 2 s −1 , and the Reynolds number Re = 27,000. The Reynolds number definition can be seen in Equation (1), where U ∞ represents the free stream velocity, δ the boundary layer thickness where the VG leading edge is located, and ν the kinematic viscosity.
The domain was automatically decomposed in eight subdomains as shown in Figure 3 in order to solve in parallel with all the workstation cores available, making the solving process faster and more efficient. Simulations were performed on a personal server-clustered parallel workstation with Intel R Xeon R E5-2630 v3 CPU 2.40GHzx16 cores and 32 GB RAM on 64 bits Linux. The simulations were completed when the residuals of the velocity, pressure, and turbulence fields reached previously defined threshold values. A qualitative convergence was considered for residuals decay to 10 −5 during the simulation. Taking into account the six VG heights and the four incident angles considered in this parametric study, twenty-four simulations were carried out. The average computational cost of each simulation took approximately 194 h. All variables needed to estimate the vortex features were extracted from eleven spanwise planes placed downstream of the VG ( Figure 1) in each simulation. The planes were normal to the main flow direction and located from 1 to 25 times the local boundary layer thickness δ downstream of the VG trailing edge (TE).

Mesh Quality Study
For all meshes used in this study, two different mesh quality parameters were studied, such as the mesh orthogonality and the skewness. The obtained values for these two parameters are shown in Table 1 for the six VG meshes. The non-orthogonality (nO) measures the angle between a line that connects the center of two adjacent cells and the normality of their common face; any non-orthogonality value satisfying the condition of nO < 70 is a safe value and for all the VG cases, the used meshes nO values lie in this safe zone. The skewness represents the distance between the point where the common face of two cells is crossed by the line that connects the center of those two cells and the center of this common face. The smaller the value for skewness, the better, and any value bigger than 20 may suggest lack of accuracy. Both meshing parameters, non-orthogonality and skewness, are in the range explained ahead for all meshes.

Results
According to Martinez-Filgueira et al. [32], the vortex center path, vortex size, and peak vorticity of the vortex are the most interesting parameters to show the characteristics of the progress of the streamwise vortexes generated downstream in rectangular VGs. Yao et al. [33] proposed the point with the highest peak streamwise vorticity as the center of the vortex core. The peak vorticity is the point which has the largest vorticity of the vortex in a spanwise plane past the trailing edge of the VG and normal to the streamwise direction. Thus, the vortex trajectory can be obtained looking at the evolution of the vortex center in the different spanwise planes downstream of the VG. The primary vortex lateral and vertical path produced by VGs have a key role in the momentum mixing performance.

Vortex Path
Vertical y and lateral z coordinates of the vortex center were calculated in spanwise planes to the streamwise direction placed at different distances x downstream of the VG. With these calculations, the vortex center deviation can be obtained for each simulated case corresponding to the six different heights.  As shown in Figure 4, a lateral evolution of the vortex center can be observed; this movement evolves along the distance x from the VG trailing edge at the considered different spanwise planes and for all the simulations. This lateral evolution shows a deviation of the vortex center in the direction indicated by the VG trailing edge. The nondimensional vortex center lateral paths observed in the figure follow a similar trend for the incident angles β = 15 • , 18 • , and 20 • , having higher lateral deviation for the highest VG height, 1.2H.
When looking at the 0.2H and 0.4H VG height cases, the deviations observed are smaller than the deviations corresponding to higher VGs. When analyzing the β = 10 • incident angle, the deviation observed is also smaller than the deviation corresponding to bigger incident angles.
Regarding the vertical path shown in Figure 5, the vortex centers calculated just past the VG trailing edge initially have a vertical position which is similar to the VG height for all the VG cases. As observed in the four plots corresponding to the four incident angles in Figure 5, there is not a meaningful vertical deviation. For the 0.8H, H, and 1.2H height cases, a relatively small downward deviation can be appreciated. On the other hand, for the two smallest height cases, the vertical path is pointing upwards. The reason for this upward deviation can be found in the fact that the lower VG vanes are within the inner part of the local boundary layer and, consequently, the viscous shear could prevail and have a determining influence when interacting with the wall. The cases with similar VG height to the local boundary layer thickness follow the conclusions of Fernandez-Gamiz et al. [24], who concluded that there is probably no vertical deviation of the vortex center position lengthwise in the streamwise direction.
In the four figures of Appendix A, Figures A1-A4, the primary vortex evolution can be observed. These four figures depict the axial velocity distribution U x at four different planes normal to the streamwise direction and located at distances proportional to the local boundary layer thickness δ for the six VG height and four incident angles.
The formation and evolution of the primary vortex is clear for all the simulations and even more evident in the 0.6H VG height case. Furthermore, the vortex size increases as the planes are farther away from the VG trailing edge. This size increase indicates a lower gradient of velocities from the vortex center to the outer part of the vortex. It can be qualitatively observed in the four figures of Appendix A that vortexes of increasing size are obtained for bigger incident angles β and higher VGs.
The lateral deviation of the vortex path can be also observed in the figures as the distance from the vortex generators' TE increases. This displacement in the z-axis occurs towards the side where the VG is rotated and is depicted for all the VG cases and is more evident for the highest incident angles β = 18 • and 20 • .
Summarizing the qualitative analysis made by means of the four figures shown in the Appendix A, the primary vortex can be clearly identified; it is simulated and captured relatively well with the exception of the smallest VG height case. where the vortex is hardly appreciated, perhaps because of its position near to the inner part of the local boundary layer.

Vortex Size
The vorticity on a perpendicular plane to the streamwise direction is a function of three variables, the distance from the vortex center, the peak vorticity of that plane, and the radius size.
The calculation of the vortex radius is a difficult task for experimental and CFD data. Martinez-Filgueira et al. [32] and Fernandez-Gamiz et al. [25] showed that the vorticity profile generated by a vortex generator in a plane normal to the streamwise direction and in a line passing horizontally through the vortex center has a Gaussian distribution-type shape and that the vortex center defined by the peak vorticity is located in the center of this Gaussian distribution.
The vorticity profile fits the Gaussian distribution better in the peak values than in the tails. A way to overcome this fitting problem in the tails comes by a innovative definition of the vortex size by means of the concept of the half-life radius (R 05 ). The R 05 is defined as the radial distance from the vortex center to the point where the half-value of the peak vorticity is located, as proposed by Bray [34].
In Figure 6, the vorticity distribution in a spanwise plane normal to the free stream velocity at x distance from the VG trailing edge can be observed. In order to obtain the vorticity profile associated to the vortex center, the vorticity values of a horizontal line (shown in white color) passing through the vortex center are extracted. For each VG case, the vorticity profile that crosses the vortex center of 12 spanwise planes normal to the streamwise direction was extracted. These planes are placed from 3δ to 25δ distances from the trailing edge of the VG.
Equations (2) and (3) represent the Gaussian distribution of the vorticity in a plane normal to the streamwise velocity, where r represents the radial coordinate, ω is the vorticity expressed in s −1 , R 0.5 the half-life radius in m and the dimensionless constant k with value of 1. According to Martinez-Filgueira et al. [32], the R 0.5 can be determined by calculating the half-value of the peak voticity ω peak .
(3) Figure 6. Vorticity distribution in a spanwise plane normal to the free stream velocity at x distance from the VG trailing edge. Primary and secondary vortexes with a vorticity field are shown in a normal plane to the free stream velocity. The white horizontal is parallel to the z-axis and passes through the vortex center. Figure 7 shows the CFD results for the R 0.5 considering the four different incident angles, β = 10 • , 15 • , 18 • , and 20 • , and some of the experimental (EXP) measurements obtained by Bray [34]. Represented values are normalized with respect to the local boundary layer thickness (δ) in the VG leading edge.
The vortex size considering the criteria of the R 0.5 shows an increasing trend along the planes placed normal to the streamwise direction. This increasing trend obtained in the CFD simulations is as well observed in the measurements of Bray [34] for the incident angles. As it gets farther from the vortex generator TE, the half-life radius gets higher values. Although the slope representing the evolution of the vortex size increases as the incident angle does, it is more dependent on the distance from the trailing edge of the VG.
Considering the VG height, bigger values of half-life radius are obtained for higher VG vanes. Even if the highest vane follows a similar trend, with each incident angle, it takes bigger values in small distances from the TE for the 10 • incident angle. The 10 • incident angle also shows a bigger difference between the half-life radius values of the highest VG case than other incident angles.
A qualitative study of the vortex size can also be done in Figures A1-A4 of the Appendix A. The velocity fields show for every VG case the vortex formation for several spanwise planes downstream of the vane and for a different incident angle in each figure. The interaction with the surface is clearly bigger for 0.2H and 0.4H VG heights, but also in these both cases, an evolution is identified where the vortex size seems bigger as the plane is further from the VG. The vortex formation gets more obvious for β = 18 • and β = 20 • . Vortexes are also generated, as expected, at a bigger height as the VGs are higher. Therefore, quantitative and qualitative results show the same conclusions about the development of the vortex size. Thus, the main vortex characteristics of this complex 3-D flow were captured and described in spite of the viscous interactions between the fluid and secondary structures and the wall. The evaluation of the CFD techniques performance of the characterization of this type of flows is important, since these flows are the most broadly used in order to solve this kind of problem.

Wall Shear Stress
The wall shear stress τ ω is a significant parameter for the understanding of the flow separation. According to Godard et al. [35], higher values of wall shear stress are directly related with the delay or avoidance of the separation of the flow .
The wall shear stress of all VG height cases and four incident angles, β = 10 • , 15 • , 18 • , and 20 • , was calculated. For higher values of the shear stress, the boundary layer gets better attached to the surface. In Figure 8, the wall shear stress τ ω is shown for four incident angles depending on the distance from the TE of the VG.   The wall shear stress shows the highest value for the 0.6H VG height case for all the incident angles. For the incident angle β = 10 • , the 0.4H VG height case starts with a big value and decreases as the cutting planes are farther from the VG trailing edge. The wall shear stress varies with the vane heights and the different incident angles. This variation is even more obvious for incident angles of β = 18 • and β = 20 • , where higher values of wall shear stress are observed for the 0.4H and mainly the 0.6H VG height cases, and as long as the downstream distance from the TE is bigger, the wall shear stress curve decreases. However, for the 0.8H, H, and 1.2H VG height cases, the wall shear stress is almost constant for the four incident angles and all along the streamwise direction, with a small variation for the incident angle β = 10 • .

Vortex Strength
As commented on at the beginning of Section 3, the vortex peak vorticity is one of the main characteristics of a vortex, so the maximum value of the streamwise vorticity was studied for cases and incident angles at a plane placed to 5δ distance from the TE of the VG, see Figure 9a. Additionally, the strength of the vortex was evaluated in the present study by the circulation (Γ 5 ) defined by Equation (4), see Yao et al. [33]. It was calculated by finding the peak vorticity over the area surrounding the vortex core in cross-flow plane positions normal to the wall and at the distance of 5δ downstream of the VG trailing edge. Figure 9b presents the values of the circulation against the device heights.
The circulation measured at 5δ roughly varies for cases where device heights are larger than 0.4H. Furthermore, the highest values are obtained for β = 18 • and β = 20 • incident angle cases, following a similar trend in these two incident angles. In agreement with the values obtained for wall shear stress in Figure 8c

Conclusions
In the present study, the primary vortex generated by a single vane-type vortex generator (VG) has been characterized. Numerical simulations were carried out using open-source OpenFOAM code for four incident angles, β = 10 • , 15 • , 18 • , and 20 • and six VG heights, 0.2H, 0.4H, 0.6H, 0.8H, H, and 1.2H, where H is the value of the BL thickness δ at the VG position. The vortexes generated by these VGs were simulated considering turbulent, steady, state and incompressible flow with a Reynolds number of Re = 27,000 based on the local BL thickness δ at the VG position.
The parameter used to measure the vortex size was the half-life radius R 0.5 , and the results showed that this value is significantly dependent on the VG height and vane incident angle. In general, the higher the VG height, the larger the size of the primary vortex. For the lowest incident angle β = 10 • studied, the largest size of the vortex is achieved by the highest VG case of 1.2H. However, as the incident angle increases, the differences in the size of the vortex in comparison with smaller vanes are not relevant.
As expected, the vortex trajectory is affected by the VG height, as well as by the incident angle β, and this path variation is more evident for the larger angles. In each simulated case, the vertical trajectory of the vortex is almost constant, with the exception of the 0.2H VG height case. This behavior may be due to the low VG height, which is within viscous sublayer where viscous effects are dominant and, therefore, an interaction with the wall appears.
Regarding the wall shear stress, two conclusions can be drawn. Firstly, a variation of the wall shear stress with the incident angle was confirmed. Larger wall shear stress values were obtained for higher incident angles β = 18 • and 20 • , and these values decrease along the streamwise direction. Secondly, the VG sizes 0.4H and 0.6H at the incident angle of β = 20 • reached the largest wall shear stress along the streamwise direction. Additionally, the circulation shows the highest values are for 0.4H and 0.6H VG heights and β = 18 • and 20 • at 5δ distance from the TE, results that agree with the obtained wall shear stress conclusions. Therefore, those vane heights and incident angles seem to be the best option for flow separation control.
CFD techniques are the most commonly methods used to improve the optimization of VGs on wind turbine blades, so the results and information presented by this study can be used as a model and guidance for effective implementation of VGs. Due to the high variety of parameters that have influence in the problem, such as VG height and orientation, a precise understanding of the flow gains importance for successful application of VGs. A further step would be the simulation of different VG sizes and orientation in typically used airfoils for wind turbine blades at different Reynolds numbers. Funding: This research was partially funded by the Government of the Basque Country through the ELKARTEK program and by the Fundation VITAL Fundazioa.

Acknowledgments:
The authors are grateful to the University of the Basque Country UPV/EHU through the EHU 12/26 research program.

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

Abbreviations
The following abbreviations are used in this manuscript: