Revisiting the assumptions and implementation details of the BAY model for vortex generator flows

Today, Vortex Generators (VGs) are becoming an integral part of a Wind Turbine blade design. However, the challenges involved in the computation of the flow around VGs are yet to be dealt with in a satisfactory manner. A large number of VG models for Reynolds Averaged Navier Stokes (RANS) solvers has been proposed and, among them, the BendereAndersoneYagle (BAY) model (ASME Pap. FEDSM99-6919) is one of the most popular, due to its ease of use and relatively low requirements for user input. In the present paper a thorough investigation on the performance and application of the BAY model for aerodynamic VG flows is presented. A fully resolved RANS simulation is validated against experiments and then used as a benchmark for the BAY model simulations. A case relevant to wind turbines is examined, which deals with the flow past a wind turbine airfoil at Reynolds number 0.87e6. When the grid related errors are excluded, it is found that the generated vortices are weaker in the BAY model simulations than in the fully resolved computation. The latter finding is linked to an inherent deficiency of the model, which is first found in this study and which is explained in detail. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Vortex Generators (VGs) are one of the most popular means of passive flow control, due to their effectiveness, simplicity and minimal cost for production and installation. Their operating principle is rather simple, as they generate streamwise vortices, which bring high momentum fluid from the free stream closer to the solid surface, energizing the local boundary layer (BL) and thus delaying separation [1]. (see Figs. 16e21) In contrast to the operating principle, the flow in the VG region is highly complex. The interaction between the generated streamwise vortices and the underlying boundary layer takes place in much smaller scales than those of the overall flow. For example, on wind turbine blades, the VG height is in the order of the local BL height, O (10 À3 m), while the rotor radius which is largest flow length scale can be up to five orders of magnitude larger, O (10 2 m). In order to accurately resolve such complex flows dense computational grids and high-fidelity turbulence models of the "Large Eddy Simulation" family are required, which ends up having prohibitively high cost for engineering usage. In fact, even with lower fidelity turbulence modeling, such as Reynolds Averaged Navier Stokes (RANS) solvers, resolving the VG geometry remains prohibitively expensive. In response to the requirement for realistic RANS CFD modeling, approximate VG models have been proposed that do not resolve the actual VG geometry, but model its effect on the flow.
VG modelling in itself is considered a challenging task [2]. In general, VG models can be categorized into three groups, (a) those modelling the effect of the shed vortex in a statistical manner (e.g. Refs. [3,4]); (b) those modelling the profile of the vortex shed by the VG (e.g. Refs. [5e9] and (c) those modelling the effect of the VG shape on the flow [10e13]. The interested reader is directed to Refs. [14,15] for a detailed literature review.
The fact that research on VG modelling remains highly active until today, e.g. Refs. [16,17], indicates that no clear winner has emerged among the different VG models. It is arguable however, that the BAY model [13], which models the effect of the VG on the flow by means of body forces, has been the most popular one in both academia and industry [14,18e24]. Its popularity comes mainly from its effectiveness, ease of implementation and relatively limited requirement for user input.
Despite its widespread use, the knowledge around its details of application and actual dependence on specific parameters remains limited, as very recent studies show [19,24]. In addition to this knowledge gap, all studies up to now compare BAY model results with fully resolved simulations performed on different computational meshes. As a result, it has not been possible to distinguish between the model error and the grid related error or to highlight inherent model weaknesses.
The present paper is focused on the industrial need for an accurate way to model VG flows using RANS solvers and presents a thorough investigation of the most popular VG model (BAY) and its application details. The main two objectives are the following: ATo investigate possible inherent model deficiencies that had been overlooked until now. BTo highlight the importance of critical model implementation details.
To this end, a comprehensive comparison exercise is performed, where BAY model simulations are benchmarked against wind tunnel experiments and Fully Resolved (FR) RANS CFD simulations. To the best of our knowledge this is the first time an inherent deficiency of the model is highlighted and that different implementations of the model are discussed in detail.
In the next section, a short outline of the BAY model derivation is given for the sake of completeness, followed by the description of the methodology adopted in the present investigation. The paper continues with the presentation and discussion of results and closes with a concluding section, where possible paths for further model development are suggested.

source term derivation
The model was originally presented in 1999 [13] and is named after the authors of the original publication, Bender, Anderson and Yagle (BAY). The VG effect on the flow is introduced via the force exerted on the flow by the VG, L ! , which is equal and opposite to the lifting force acting on the VG. In this sense, the VG is replaced by a thin lifting surface, the presence of which leads to generation of streamwise vorticity [10,25]. According to thin airfoil theory, the lift coefficient for the VG will be 2pb, where b is the angle between the VG and the local velocity vector. Hence the force acting on the fluid is expressed by Eqn. (1): where r is the local density, u ! the local velocity, S VG the VG planform area and b l is the unit vector in the direction of the force. In the BAY model formulation, b l is expressed based on the unit vectors describing the orientation of the VG surface. As shown in Fig. 1, left, and Fig. 2 these are b t; which is tangent to the VG chord, b b, which is parallel to the VG span and b n, which is normal to the VG surface. The unit vector b l is then defined as normal to the local velocity unit vector, b u, and the spanwise vector b b: Using the small angle approximation, the angle b is: To approximate the loss of lift at higher angles of attack of the VG with respect to the free stream the force source term is also multiplied by a factor of b u$ b t. In a finite volume context, the total force acting on the fluid will be where L ! i is the source term added to the momentum and energy equations at the cells where the model is applied. L ! i is given by Eq. (5): where c VG is the BAY model constant and V i is the grid cell volume.
The momentum and energy equations at the cells where the model is applied become: where t is time, E is the total internal energy, F ! Mj ; F Ej represent the sum of the convective and diffusive fluxes of the momentum and energy equations and S j is the cell face area.
In its jBAY variation, which is an elegant simplification of the model, the source term is only applied to the cells engulfing the VG geometry (VG cells), as first described in Ref. [10]. An example of VG cell selection is given schematically in Fig. 1, right.
The constant c VG is the only model parameter and acts as a relaxation parameter, controlling the strength of the force term. According to published literature, the model is insensitive to the value of c VG for sufficiently large values of the constant, provided that the total VG cell volume P V i is close to the real VG volume. However, there is currently a disagreement between researchers regarding what is large enough for c VG values. Anything from c VG > 1 to c VG ¼ 15 has been used and despite the apparent disagreement no satisfactory explanation has been given to date [10,14,18,19,21,22,26]. For the rest of this paper, and unless otherwise stated, a value of c VG ¼ 10 is used, in order to be in-line with current research practice. The effect of the model constant is examined in detail in Section 4.5.

Methodology
The BAY model offers a simplified way to model the flow past vane VGs in RANS simulations and, hence, it can only be as good as a fully resolved VG RANS simulation, at best. Therefore, in order to quantify the accuracy of the BAY model (BAY model error), BAY and FR simulations with the same RANS solver on the same grid are compared. Using the same grid, excludes grid related errors while using the same RANS solver excludes errors due to turbulence modeling. Furthermore, in order to establish FR simulations as the basis of comparison (RANS error), they are first compared to experimental measurements and similar FR simulations from the literature. In the last stage of the present study the effect of grid on the BAY results (Grid error) is examined within a more realistic context. The identified types of error are detailed in Fig. 3.
It is noted at this point that the grid selection for the BAY model error identification is so dense that it is not a realistic scenario for an industrial investigation, e.g. a VG parametric study, as it cancels out the cost-reducing benefits of using the model in the first place. However, it is a valuable research exercise that highlights model deficiencies excluding grid related issues. To the best of the authors knowledge, all relevant BAY model investigations up to now always compared FR and BAY simulations performed on different grids and thus were unable to distinguish between the grid related and the model related errors.
The grid dependence of the model is then further examined using more realistic, extruded meshes of different spanwise density. The effect of grid architecture and density is discussed based on these results.
The experimental data used for the validation include pressure and Stereo Particle Image Velocimetry (PIV) measurements [27] and have already been used for similar validation purposes [14,16,21]. Numerical data from a previously published FR simulation of the same flow case [14] are also used to validate the present FR computation. In that case simulations were independently performed by researchers from the Danish Technical University using the incompressible in-house finite volume RANS solver EllipSys3D [28].
All present simulations are performed using the in-house finite volume RANS solver MaPFlow [29] and the k-u SST turbulence model [30]. MaPFlow solves the compressible RANS equations, is second order accurate in time and space and is equipped with preconditioning for low Mach flows. Riemann invariant boundary conditions were used for the far-field, symmetry conditions at the sides and no-slip condition on the airfoil (and the VG for the FR case). In all cases a steady state, fully turbulent flow was

considered.
Although the accuracy of RANS simulations in vortical or highly separated flows is debatable, using higher fidelity methods (e.g. LES) can be computationally prohibitive for the Reynolds number considered in this work. In the context of the present study, where the implementation details and inherent shortcomings of industrially relevant options are examined, only RANS simulations are considered.

The examined case
The main case considered is that of a wind turbine airfoil (NTUA t18, see Ref. [27] and Fig. 4

top left) equipped with VGs at chord
Re ¼ 0.87e6 and a ¼ 10 . The VGs are located at 30% chord and their height, h, is equal to the local BL height, d. The set-up is counterrotating with common flow up. The VG dimensions are given in Fig. 4 for completeness. The more demanding case of a ¼ 20 is also examined using all available grids while force coefficient polar data for 10 < a < 20 are also discussed.

Fully resolved simulation validation
The FR simulation is validated against the published wind tunnel data. Numerical mesh details are given in Table 1, while details of the grid are shown in Fig. 5. A single vane VG was considered with symmetry boundary conditions at the sidewalls of the computational domain. Fig. 6 and Fig. 7 show vorticity and velocity contour comparisons   between the present FR predictions, the experimental data (Stereo PIV [27]) and the EllipSys3D FR data (EllipSys3D, [14]), respectively. The plane considered is at 60% chord, or 30 VG heights downstream of the VG LE, the wing angle of attack is a ¼ 10 and the Re number The agreement is good, especially between the two FR solutions. This comparison confirms that the MaPFlow simulation can be used as a benchmark case for the present BAY model investigation. As explained in Ref. [14], the disagreement in the absolute vorticity values in Fig. 6 is due to the different fact the quantity is calculated on different grids that are used in the PIV and CFD data sets. (SPIV measurement grid vs. CFD grid).

BAY model vs. Fully resolved simulation on the same grid
Before comparing the BAY predictions to the FR ones, the FR flow evolution along the VG is further examined in order to highlight the vortex generation mechanism. It is acknowledged at this point that RANS turbulence models are not ideal for vortical flow prediction, since they have not been calibrated for such flows. On the other hand, RANS errors are not dominant in the vortex formation, which is an inviscid process [2]. Fig. 8 shows vorticity contours with Q criterion [31] isolines to highlight vortex regions. In total, four vortices are observed:  The main vortex is shed off the VG tip, the horseshoe vortex is the roll-up of the incoming boundary layer and the other two vortices are induced by the main vortex [32]. The only vortex that is upstream of the VG is the horseshoe vortex, while all other vortical structures are downstream of the VG lee side. Identifying the distinct vortices is important in order to pin point the differences in the vorticity generated by the BAY model, as discussed below. Fig. 9 compares the early stages of the vortex generation mechanism between the FR and the BAY simulations. The vorticity contours and Q criterion isolines reveal that in the BAY model simulation all vortices are weaker in the region close to the VG  (from x ¼ 0.335 to x ¼ 0.355). However, further downstream, at x ¼ 0.8, the situation reverses as indicated in Fig. 10. The vortex in the BAY modeling appers larger, stronger and closer to the airfoil surface compared to that in the FR simulation. Fig. 11 depicts vorticity isosurfaces in a three-dimensional view that better pictures the different vorticity evolution in the BAY and FR results. In the BAY modeling, the horseshoe and the two induced vortices are significantly weaker compared to those in the more accurate FR simulation. This results in different flow orientation as shown by the surface streamlines shown in Fig. 12. Contours in these plots refer to the spanwise velocity (W). The most pronounced differences are found on the left of the VG (circled by a red double line), in the VG wake (circled by a black dotted line) and at the junction of the VG lee side with the wing surface (indicated by a white arrow).
As discussed in detail in the next section, in the BAY model case, the flow is deflected because of the body forces, but is never perfectly aligned with the VG direction. The fact that the flow is significantly deflected leads to the creation of a vortex system that is qualitatively similar to that of the FR case. The quantitative agreement is not possible because the flow is not perfectly aligned with the VG direction.
The fact that in the BAY model simulation the main vortex is weaker closer to the VG (Fig. 9), but stronger towards the wing TE ( Fig. 10) is counterintuitive. Since grid, turbulence modeling or solver related errors cannot be considered as the cause, a plausible explanation is that in the FR case all four vortices are stronger and that strong interactions between them result in momentum loss. Additionally, the increased strength of the main vortex closer to the VG justifies its higher position, as it is pushed upwards by its "mirror" vortex from the symmetry condition on the computational domain sidewalls. Fig. 13 shows the integral force imposed by the VG on the flow as calculated in the two simulations. The FR force is the sum of pressure and friction forces on the VG. In agreement with [19], the pressure force contribution is dominant and the friction force is significantly smaller, equal to 3% of the pressure force. Despite being small, this friction contribution is sufficient to modify the integral force orientation which is not normal to the VG, like the pressure force, but at 88 .
In the BAY results, the force is the sum of all force terms from all VG cells. The disagreement is obvious in both direction and magnitude. The BAY model force is stronger by 53% and at a different angle (100 ). This difference leads to the different vortex strengths close to the VG surface as further discussed in the next section.
One apparent cause for this disagreement is that the BAY model does not consider the friction forces from the VG surface, which are small, but quantifiable. BAY model only attempts to model the pressure force, by adding a term normal to the local velocity vector.    A second cause for this discrepancy is that, due to an inherent model deficiency, this force never has the right direction or distribution. This is further discussed in Section 4.4.

BAY model on extruded grids
The BAY model was further investigated with regard to its grid dependence. Three different grids were considered, which differed only in the spanwise resolution. The same 2D airfoil grid with~50k cells was extruded 10, 20 and 70 times until the end of the computational domain to generate the Coarse (0.5e6 cells), Medium (1.0e6 cells) and Fine (3.5e6 cells) grids. It is noted that the latter grid has a very similar total number of cells with the FR grid (3.6e6 cells) examined ion the previous section. However, the meshing strategy is significantly different with the FR grid resolving the VG geometry, whereas the Fine grid is a simple extrusion of a 2D grid. Table 2 presents the integral force results with respect to the VG force, as computed from the FR simulation, for two airfoil incidence angles, 10 and 20 . Results show that the BAY integral force direction and magnitude depend heavily on grid density and meshing strategy. Fig. 14 graphically displays this dependence for the 10 case, while Fig. 15 shows surface streamlines on the airfoil surface and BAY model cells for the three extruded grid cases. For a ¼ 10 , the force magnitude increases for coarser grids, as the total VG volume, i.e. the volumes sum of the cells that the BAY term is added to, is significantly higher, see also Fig. 15.
Surface streamlines in the vicinity of the VG also significantly change over the different grids. In the Coarse and Medium grids, the flow seems to pass through the VG close to airfoil surface, while in the Fine and the FR grids, the flow patterns of the BAY and FR simulations are very similar (compare Figs. 12 and 15, top). This suggests that for the Coarse and Medium grids, the local velocity close to the airfoil boundary deviates significantly from the VG direction. Fig. 16 displays the contour of the local velocity angle (angle b in Fig. 2) with respect to the VG direction for the Medium grid case. Indeed, close to the wing surface the local flow deviates from the VG direction by approximately 20 . In fact, even at the centre of the VG surface, the flow misalignment ranges between 2 and 4 , while at the VG tip the angle difference reaches 10 .
Possible explanations for this discrepancy are given below. Close to the wall, the local velocity magnitude is very small, and so is the j u ! j 2 term in Eq. (6), resulting in a small lift force which is not strong enough to align the flow with the VG. At the tip cells, it is conceivable that the flow deviation is caused by the local acceleration above the VG and the relative increase of momentum in the region.
The differences in the region close to the airfoil surface and at the tip is what leads to the generation of a weaker vortex system, as discussed in the previous section. Because the force is not sufficiently large locally, part of the flow goes through the VG cells and the BL roll-up that leads to the formation of the horseshoe vortex is not as intense. Also, since the BAY force is always normal to the local velocity, see Eq. (5), the flow misalignment means that the BAY force is not normal to the VG surface, which in turn explains the differences in the integral force angle shown in Figs. 13 and 14. According to Bender at al [13]. "… the model source term starts to dominate the other terms in the discrete finite volume equations. When this happens, the flow tends to align itself with the orientation of the vortex generator, such that the local angle of attack approaches zero, i.e. b ¼ b u$b nz0, thus allowing the discrete momentum equations to maintain equilibrium." According to the presented results however, what happens is that angle b gets small, but non-zero values (Fig 16), such that b ¼ b u$b nz0 and the applied force on the flow is never equal to zero, jL i js0. In fact, there is a dynamic equilibrium where in one step the force turns the flow towards the VG direction (reducing b) and in the following step the force magnitude is reduced as a result. Because of the reduced force magnitude, the flow now turns away from the VG direction and b is now increased. In the following step, because of increased b, the force magnitude is also increased and the flow is again turned towards the VG direction and so on. The result is that the flow is never completely aligned to the VG direction at the VG cells. Also, this oscillating behavior is very small in magnitude and only concerns a limited number of cells, so it does not affect solution convergence, which is satisfactory (density residual at the order of eÀ11).
In fact, the FR simulation suggests that the integral VG force is not normal to the surface and that its orientation is also angle dependent. This is because in addition to the lift force, the VG also experiences drag. As explained in Section 2.1, the BAY model source term only considers the lift force according to thin airfoil theory and completely disregards the drag force.
Investigation of the effect of the model constant, C VG .
Researchers agree that, provided the sum of all VG cell volumes is close to the actual VG volume, the BAY model is insensitive to the value of c VG for sufficiently large values of the constant. Nonetheless, there is a disagreement on the definition of what sufficiently large is and values from c VG > 1 to c VG ¼ 15 have been proposed [10,14,18,19,21,22,26].
In fact, if the model did align the flow with the VG direction locally when the solution converged, then it would indeed be independent of the constant value as the force term would be zero. However, as explained in Section 4.4, the force term never becomes  zero since the flow is never aligned with the VG direction. As a result, the constant value will have an effect on the final force value.
Simulations for a range of c VG values (0:1 c VG 100) were performed for the Medium and the Fine grids and the results are presented in Figs. 17 and 18. The integral force for each case is scaled with the integral force for the maximum c VG value, c VG ¼ 100. It is evident that convergence is reached for high c VG values, with the force reaching 95% of the converged value for c VG ¼ 10. On the other hand, for smaller c VG values the force is significantly underestimated. For very high values (c VG ¼ 100) the force angle is almost normal to the VG (b ¼ 91 ), but it never becomes zero as the model would predict. It is noted that the local velocity angle distribution is similar to that shown in Fig. 16, even for the highest c VG value, with the greatest misalignment observed at the base and the tip of the VG. Fig. 18 is a graphical representation of the data presented in Fig. 17, showing that the force and its angle with respect to the VG increase as c VG increases.
It is important to note at this point that for high c VG values (c VG > 20Þ, numerical convergence is harder to reach. Increasing c VG results in a stronger source term and this hinders convergence, especially during the initial stages of the solution procedure, when the flow is completely misaligned with the VG. To overcome this, CFL values ten times smaller had to be used for c VG > 20, which made the computations even more expensive. A possible solution to this could be a gradual increase of c VG in the course of iterations to facilitate convergence.
The effect of the various grid options on the aerodynamic coefficients of the examined airfoil is discussed in this section. Fig. 19 shows the lift and drag coefficients variation with angle of attack for the FR simulation and the different BAY simulations. Only the region around C lmax is shown as this is of the greatest interest.
The BAY simulation on the FR grid has a constant offset with respect to the FR polar, having higher lift and lower drag values. It follows the same trend as the FR simulation and has the same a Clmax , at a ¼ 14 . Skin friction contours and surface streamlines for the 10 case (Fig. 20) show that the lift offset is because the flow is less separated in the BAY simulation than in the FR results. Conceivably, this happens because the BAY vortex is stronger and closer to the wing surface towards the wing TE, as explained in Section 4.2. As a matter of fact, however, the FR simulation results agree better with flow visualization data from Ref. [27]. The reduced drag values are also justified by the smaller disturbance that the BAY model causes to the flow (weaker vortices close to the VG), as discussed earlier.
With regards to the extruded grids, Figure 20 shows that the Coarse grid causes excessive separation even at 10 , while at the more challenging 20 case (Fig. 21) even the Medium grid leads to highly separated flow unlike the FR simulation. This also explains the force discrepancies in Table 2 for the 20 case.

Conclusions and future research
A thorough investigation on the application of BAY model for Vortex Generator flows has been performed. A Fully Resolved RANS simulation, validated against experiments and other FR simulations from the literature, has been used as the benchmark case.
It is found that even on the same grid as the FR simulation, the BAY model produces weaker vortices in the VG region. The main vortex, the horseshoe vortex and the induced vortices are all weaker, arguably because the flow is never fully aligned with the VG direction. This is an inherent deficiency of the model as explained in detail in Section 4.4.
Despite the fact that initially the vortices are weaker for the BAY model simulation, the main vortex is stronger towards the airfoil TE compared to the benchmark case, due to weaker vortex interaction. This leads to reduced separation and increased lift compared to the FR simulations, which are found to agree very well with the experiments.   With regards to the grid dependence it is found that the integral BAY force depends on both grid density and grid architecture and that as pressure recovery requirements increase, denser grids are required to achieve accuracy.
Finally, a study of the BAY model constant, c VG , showed that, in agreement with current understanding, for high enough constant values, model performance is independent of the actual constant value. It was also shown that current in current practise, values below that critical limit are used, possibly because this facilitates convergence.
The present study aims to act as a base for future modifications and improvements of the BAY model. Based on the results presented herein, future directions to be explored could include: A different selection of the velocity magnitude term (j u ! j 2 ) in Eq. (5) (e.g. the average VG velocity instead of the local cell velocity).
ο This could reduce the flow misalignment close to the wing surface. ο A distributed c VG along the VG height, with an increase close to the wing surface and at the tip. ο This could also increase the force magnitude closer to the wing and reduce the flow misalignment. ο The addition of a drag component to the lift BAY force.
A The FR simulations showed that the resulting force is not normal to the VG surface. The BAY force would potentially be closer to the FR results if a calibrated drag component was added.
ο A calibrated angle offset to the direction of the force. ο This modification could also possibly improve the alignment of the integral force.
A gradual increase to the c VG value along the course of the iterative solution procedure.   ο This would allow the integral force to reach higher values, and hence better align the flow with the VG, without hindering convergence.