Source Term Modelling of Vane-Type Vortex Generators under Adverse Pressure Gradient in OpenFOAM

An analysis of the generation of vortices and their effects by vane-type vortex generators (VGs) positioned on a three-dimensional flat plate with a backward-facing ramp and adverse gradient pressure is carried out. The effects of a conventional vortex generator and a sub-boundary layer vortex generator are implemented by using a source term in the corresponding Navier-Stokes equations of momentum and energy according to the so-called jBAY Source Term Model. The influence of the vortex generator onto the computational domain flow is modelled through this source term in the Computational Fluid Dynamics (CFD) simulations using the open-source code OpenFOAM. The Source Term Model seems to simulate relatively well the streamwise pressure coefficient distributions all along the flat plate floor as well as certain parameters studied for vortex characterization such as vortex path, decay and size for the two vane-type vortex generators of different heights studied. Consequently, the implementation of the Source Term Model represents an advantage over a fully Mesh-Resolved Vortex Generator Model for certain applications as a result of a meaningful decrease in the cell number of the computational domain which implies saving computational time and resources.


Introduction
The need for bigger capacities and installed power of wind generation systems has allowed the implementation of larger rotor wind turbines. Nowadays, wind power generators of 5-7 MWatt may present rotor blades of 60 m long or even larger.
This meaningful increase in the rotor size and weight of wind turbines has made that the techniques used in the past years to control them have to be improved. Johnson et al. [1] summarized a big part of the most relevant flow control techniques applied in wind turbines to work in an optimal and safe way under diverse atmospheric conditions. Several flow control devices have been developed in recent decades, Gad-el-Hak [2]. Most of them were firstly designed for aeronautical applications and consequently the initial starting point of research was aeronautics, Taylor [3]. According to Liu et al. [4], flow control devices are regularly also applied in turbo machinery systems. In recent years, researchers have been trying to optimize and install these devices in multi-megawatt wind turbines. Wood [5] and Johnson et al. [1] developed a four layer scheme to conceptually classify these flow control devices.
These control devices can be classified as active or passive ones depending on the control techniques applied as proposed by Aramendia et al. [6]. Passive flow control techniques enhance the turbine's efficiency and load reduction without external energy consumption. Active control techniques need a supplementary energy source to obtain the effect on the air flow and, unlike vortex generators (VGs) and other passive devices, active flow control requires elaborate algorithms to get the most favorable gain (Becker et al. [7]). Johnson et al. [1] carried out an investigation on around 15 control devices for wind turbines. Some of these flow control devices are still checked nowadays on full-scale wind turbines in working conditions. Vortex generators (VGs) are passive flow control devices that modify the boundary layer fluid motion exchanging momentum from the outer into inner region. Due to this energy transport, the velocity of the inner region increases at the same time as the boundary layer thickness decreases and in turn the flow separation is delayed, Rao et al. [8]. Moreover, Lin et al. [9] proved the effect of drag reduction and lift increase induced by sub-boundary layer vortex generators (VGs).
According to Doerffer et al. [10] and Steijl et al. [11], the streamwise vortices generated by vane-type VGs produce an endless momentum to counter the momentum decrease of the natural boundary layer and the growth of its thickness induced by viscous friction and adverse gradients of pressure. These vane-type VGs are able to reduce or remove flow separation when there is a limited adverse pressure gradient as commented by Velte et al. [12,13]. If the flow separation even appears for cases of large adverse gradients of pressure, the vortices action of mixing will limit the reversed flow region in the shear layer and maintain, to some extent, pressure recovery all along the detached flow.
Fernández-Gamiz et al. [14] and Urkiola et al. [15] analysed the flow effects of a vane-type VG on a flat plate under conditions of negligible streamwise pressure gradient and the primary vortices generated to examine how the physics of the wake past the VGs can be simulated by means of Computational Fluid Dynamics (CFD) computations. In those studies, the CFD simulations were compared with experimental observations. For instance, Gao et al. [16] and Baldacchino et al. [17] focused on the increase of the maximum value of the lift coefficient of a 30% thick DU97-W-300 airfoil designed by the Delft University of Technology by implementing passive VG devices. When increasing the angle of attack, the lift and drag coefficients raised up to values higher than the ones obtained under steady-state conditions.
The so-called micro vortex generators were introduced by Keuthe [18]. These flow control devices were also named as sub-boundary layer vortex generators (SBVGs) according to Holmes et al. [19], submerged vortex generators in Lin et al. [20], low-profile vortex generators by Martinez-Filgueira et al. [21] and micro vortex generators in Lin et al. [22]. Kenning et al. [23] summarized in his study the potential applications of all these vortex generators, including the control of shock-induced separation, smooth surface separation or leading edge separation.
Vortex generators are fixed onto the blades of wind turbines in order to prevent or remove the flow separation as well as to decrease the blade sensitivity. VGs are usually mounted on the suction side of the blade in a spanwise array for instance when a wind turbine does not perform as expected. Accordingly, the implementation of VGs onto the blade could eventually be a valid way to enhance the efficiency of a wind turbine rotor, Schubauer et al. [24] and Bragg et al. [25].
In order to design and optimize the position of the vortex generator on a wind turbine blade, CFD tools can be used as for instance Fernández-Gamiz et al. [26] or Troldborg et al. [27] do. Nonetheless, modelling a full rotor by using a fully-meshed VG computation becomes excessively expensive. As a matter of fact, the size of the vortex generators is generally comparable to the local boundary layer thickness and this implies that a meaningful number of cells around the device in the computational domain is needed in order to accurately model the flow. An alternative way of modelling VGs in the CFD computations is to model the device effect on the boundary layer by using body forces.
Bender et al. [28] established a source term model, the so-called BAY model, based on the Joukowski Lift Theorem and the Thin Airfoil Theory which used body forces. This model was designed for simulating vane-type Vortex Generators in a finite volume Navier-Stokes code and ignores the condition to fully define the geometry in the mesh. In order to calibrate the model, Bender [28] created a test case to compare the results obtained for a source term modelled VG with those obtained for a fully-meshed VG. The test case was based on 24 VGs circumferentially mounted on a pipe in a co-rotating configuration and the results obtained were satisfying. For instance, Dudek [29] successfully evaluated the BAY model for different types of flows in a rectangular duct with a single vane-type VG.
Recently, a better version of the BAY model was designed by Jirasek [30], namely the jBAY model. The improved version was based on the Lift Force Theory based on Bender et al. [28] and brought a more suitable technique for simulating the flow when using VG arrays. Jirasek [30] applied a reduced method for determining the model control points to implement in a effortless way the model and obtain more correct results. The new version of the model was checked by reproducing a single vane-type VG on a flat plate in an S-duct air intake and a high-lift wing configuration. The results showed very good agreement between experimental data and the CFD simulations. Thus, Florentie et al. [31] studied the potential of the jBAY model to reproduce the effects of rectangular VGs positioned on a flat plate and proposed a model optimization.
The main goal of the current study is firstly the implementation of the so-called jBAY Source Term Model into OpenFOAM [32] and secondly to investigate how well the open source CFD simulations are able to mimic the physics of the flow behind modelled rectangular conventional and sub-boundary layer vortex generators (VGs) mounted on a three-dimensional flat plate with a backward-facing ramp under adverse pressure gradient conditions. The satisfactory implementation of the Source Term Model in OpenFOAM would be eventually able to bring new model applications of this open source CFD code. The theoretical background of this study is described in Section 2; Section 3 presents the experimental data; Section 4 describes in detail the computational configuration; the results and their discussion are covered in Section 5 and the conclusions are presented in Section 6.

Vortex Generator Setup
The model consists of a single rectangular VG positioned on a flat plate with a backward-facing ramp. The three-dimensional computational domain modelled in OpenFOAM represents the extended test section of the wind tunnel experiment performed by Lin [33] and is shown in the illustration of Figure 1. Note that the single VG is placed upstream the backward-facing ramp. The domain has been designed following a previous simulation study by Konig et al. [34]. The geometry dimensions of the vane-type VG determined by a length L two times the VG height H are depicted in Figure 2a for two different heights H 1 = 0.8δ and H 2 = 0.2δ where δ is the local boundary layer thickness just upstream edge of the ramp. The two vortex generators, namely, conventional VG and sub-boundary layer VG of heights H 1 = 0.8δ and H 2 = 0.2δ respectively are positioned at distances of 5δ and 2.5δ from the upstream edge of the ramp. Figure 2b illustrates a classical configuration example of a pair of counter-rotating VGs. The rectangular vane has a constant thickness and no sharp edges. The parameter λ represents the spacing between VGs and lines A, B and C are measurement lines along the floor past the VG. The geometry and parameters of the two VG cases are summarized in the table of Section 4. Note that the computational domain has only one vane-type VG oriented to the main flow instead of two vanes in order to reduce the meshing and the computational resources needed. A boundary layer region develops over the flat plate induced by the viscous interaction between the wall and the air flow. The angle of incidence between the main flow direction and the VG vane is β = 15 • or β = 25 • depending on the VG case (see Figure 3) and the Shear Stress Transport (SST) turbulence model Menter [35] has been elected because of its capacity to solve whirling flows according to Liu et al. [4].   The Reynolds number for the CFD simulations is Re θ = 9100 in order to mimic the experiments explained in the next section and is based on the local boundary local (BL) momentum thickness θ and computed by the following expression: where U ∞ = 40.23 m·s −1 is the free stream velocity, ν = 1.9439 × 10 −5 m 2 ·s −1 the kinematic viscosity and the momentum thickness θ defined as: where u is the streamwise velocity component and y the vertical coordinate normal to the wall. Thus, the momentum thickness of the local boundary layer δ at the upstream edge of the ramp calculated by the previous equation is θ = 4.3967 mm.

Source Term Model
Computational Fluid Dynamics (CFD) is a commonly method used to reproduce the air flow downstream the VGs and to estimate the blade efficiency. Furthermore, the combination of this technique of simulation and the corresponding experiments in wind tunnels can be considered as a reliable tool for parametric studies on VG layout. However, these computational methods are time-consuming when designing high-quality meshes and solving the corresponding flow equations.
In the current work, the so-called jBAY Source Term Model presented by Jirasek [30] and based on the widely used BAY Model developed by Bender et al. [28] is implemented. The BAY model was designed for simulating vane-type vortex generators into finite volume Navier-Stokes codes and allows substituting the VG geometry by a subdomain of similar size at the original VG location where a specific body force distribution is then applied. The model is integrated into the CFD code as a so-called source term in the Navier-Stokes momentum and energy equations.
In the model, a force normal is applied to the direction of the local flow and at the same time parallel to the surface. This normal force simulates the side force generated by a vane-type vortex generator. Bender et al. [28] designed this Source Term Model based on the Joukowski Lift Theorem and the Thin Airfoil Theory in order to model the VG influence. Taking into account a rectangular vortex generator (VG), the lift forces on the VG can be computed by the following expression: where − → L is the lift force on the VG, ρ the local density, u the local velocity, b the unit vector defined as b = n× t with n and t the unit vectors normal and tangential to the VG respectively (see Figure 3), h VG the VG height and Γ the circulation defined as: where β is the angle of incidence between the flow direction and the VG and l VG the VG length. As expressed in Equation (3), the direction of the lift force − → L on the VG is obtained by means of the product of the local velocity u and the unit vector b. According to the Joukowski Lift Theorem for 2D airfoils: where S VG is the vortex generator area ( l VG × h VG ) and the lift force on a single cell − → L cell is computed by the following expression: where V cell is the volume corresponding to a single cell and V s the total volume of the grid cells where the model is applied. According to Bender et al. [28] a new term ( u· t u ) is introduced in the previous equation in order to obtain the final lift force on a single cell: When applying the BAY Model, a calibration process is required and a Mesh-Resolved VG Model can be used as reference for model calibration. Particularly for this model, the empirical constant C VG of the Equation (7) is a relaxation parameter. This parameter is selected so that the simulated estimations for the lift force distribution fit the results corresponding to the computations of the mesh-resolved VG. Bender et al. [28] proposed that this Source Term Model can operate in two modes, the so-called asymptotic and linear ones. If V s , defined as the total volume of the grid cells where the force term is applied, differs considerably from the volume of the vane type vortex generator, the model behaves in linear mode and is dependent on the constant C VG . However, the model of the present work has been locally applied, such that the volume V s is similar to the vortex generator volume. Therefore, the model operates in the asymptotic mode and is non sensitive to C VG . Commonly values for this constant are C VG ∼ = 10, see the study by Jirasek [30]. In the present work, instead of applying forces in all the subdomain cells as the BAY Model specifies, the forces are applied in cells placed just in the outline of the VG geometry, see Figure 4.  Taking into account that the Source Term Model was designed to be user-friendly, any vane size of height H and length L and angle of incidence β can be modelled. When using the OpenFOAM CFD code, only three parameters have to be specified to model the vortex generator (VG): the grid cells where the Source Term Model will be applied, the vortex generator area S VG and the angle of incidence β of the main flow direction with respect to the VG orientation. In Figure 5, two examples of subdomains of the selected mesh cells which represent the two vortex generators (VGs) of heights H 1 = 0.8δ and H 2 = 0.2δ in the Source Term Model are shown.  The reader should finally note that the Source Term Model can replace the Mesh-Resolved VG Model when modelling VGs in certain cases. In a fully Mesh-Resolved VG Model, a fine mesh in the vicinity of the vortex generator (VG) must be designed in order to capture properly the VG effect on the domain flow. Thus, Fernández-Gamiz et al. [26] designed a fine mesh of approximately eight million cells for simulating the primary vortex induced by a rectangular sub-boundary layer VG on a flat plate in a fully meshed-resolved model. When using the BAY Model, a coarser mesh can be used to capture the VG influence on the flow saving computational time and resources as will be indicated later. However, attention must be paid when defining the mesh because the BAY Model accuracy can eventually be dependent on mesh resolution according to some authors, e.g., Florentie et al. [31].

Experimental Data
Experimental data for this study has been taken from the experiments carried out by Lin et al. [33]. The experiments were performed in the NASA Langley 20 by 28-inch Shear-Flow Control Tunnel and free stream velocity of U ∞ = 40.23 m·s −1 . Baseline flow separation was determined on a 25 • sloped backward facing curved ramp of radius R = 20.32 cm with a behaviour similar to a classical two-dimensional flow. The ramp spanned the wind tunnel test section width. The local boundary layer thickness δ was approximately 32.51 mm and the Reynolds number was Re θ = 9100 at the upstream edge of the sloped backward facing curved ramp which was placed approximately at 1.98 m from the wind tunnel test section entrance. In the experiments, different types of passive flow control devices were analysed and between them vortex generators (VGs). Furthermore, sets of vortex generators (VGs) were simulated upstream of the backward-facing ramp. In the present study, only the pressure distribution data corresponding to a counter-rotating conventional vane-type vortex generator ( Figure 2) of height H 1 = 0.8δ and a sub-boundary layer vane-type vortex generator of height H 2 = 0.2δ where δ represents the local boundary layer thickness at the upstream edge of the ramp will be analysed for comparison with the data obtained from the OpenFOAM simulations. The higher vortex generator of height H 1 was positioned at a distance of 5δ upstream of the ramp and the smaller vortex generator of height H 2 at a distance of 2.5δ.

Computational Configuration
Non-commercial open source code OpenFOAM [32] has been used for simulating the vortex and its effects in the present study. The open source CFD code is an object-oriented library package designed in C++ programming language to analyse systems of computational continuum mechanics.
The current computational domain ( Figure 6) partially based on a previous simulation study by Konig et al. [34] consists of a flat plate with a backward-facing ramp and adverse pressure gradient where selected cells corresponding to a single rectangular vortex generator according to the Source Term Model namely jBAY previously presented in Jirasek [30] and positioned at a point where the VG height is 80% and 20% the local boundary layer thickness δ at the upstream edge of the ramp for the cases H 1 and H 2 respectively. The incidence angle of the vane with respect to the oncoming flow is β = 15 • for the H 1 = 0.8δ vane-type VG case and β = 25 • for he H 2 = 0.2δ vane-type VG case as shown in Figure 3. The modelled VGs have an aspect ratio (AR) defined as the relationship between the height H and the length L of AR = 0.5 (Table 1). A representation of the mesh subdomain and selected cells representing the vortex generator is shown in Figure 5. The origin of the computational domain was placed at the ramp center.  The variables of particular interest which require boundary conditions are the fluid quantities, velocity and pressure. Non-slip boundary conditions were applied to the flat plate and ramp of the computational domain to allow for boundary layer growth (see Figure 1). The sides of the computational domain have been defined as symmetry condition. The domain consists of only one vane and the symmetry assumption used in the present computations can be justified by previous studies as Gutierrez-Amo et al. [36] and Sørensen et al. [37]. The atmospheric conditions were selected in such a way that a boundary layer thickness of approximately δ = 32.51 mm just upstream edge of the ramp according to Lin [33] was attained in the simulation of the test section of the wind tunnel. As the k-ω Shear Stress Transport (SST) turbulence model designed by Menter [35] was applied, a requirement when setting boundary conditions for the turbulent viscosity ν t , the turbulent kinetic energy k and the dissipation rate ω was needed. Nearby the flat plate and ramp surfaces, wall functions were not used in the computations. Slip-wall boundary conditions were also implemented at the beginning and end of the flat plate ( Figure 1) according to Konig et al. [34].
Two CFD solvers were applied in this three-dimensional simulation. The potentialFoam solver was firstly used to create initial fields to accelerate the convergence process and secondly the simpleFoam solver for the steady-state, incompressible and turbulent flow using the RANS (Reynolds Average Navier-Stokes) equations. This second solver is based on the k-ω SST turbulence model by Menter [35]. This turbulence model combines the Wilcox k-ω model for the close-to-wall regions and the k-model for the outer regions. The normal forces are applied in the selected cells of the computational domain according to the Source Term Model explained in Section 2 using the topoSet tool.
The computational domain analysed was discretized with a structured-type mesh and hexahedral faces of around 5 × 10 5 cells for the H 1 = 0.8δ vane-type VG case and 2 × 10 6 cells for the H 2 = 0.2δ vane-type VG case. Full second order linear-upwind schemes for the discretization were implemented for all the simulations. Two quality mesh parameters such as mesh orthogonality and mesh skewness for the two VG cases have been studied in order to analyse the mesh quality and the results are shown in Table 2. Mesh non-orthogonality is an indicator of how orthogonal the pairs of neighboring cells which share a face are and mesh skewness of how optimum the cell shape is in relation with the corner angles. The non-orthogonality parameter should be close to 0 degrees and the maximum skewness should not exceed 0.85 as suggested by Gutierrez-Amo et al. [36] for hexahedral cells to obtain an accurate solution. Table 2. Quality meshing parameters corresponding to the mesh of the two vane-type VG cases.

VG Case Non-Orthogonality (Average) Maximum Skewness (-)
According to Richardson and Gaunt [38] Extrapolation Method, a grid dependency study was performed for three different mesh resolutions: coarse, medium and fine. The method was applied to a vortex parameter, namely the normalized peak vorticity (ω x max δ)/U ∞ analysed in the next section at a plane normal to the streamwise direction located at x/δ = −3 past the vane for the H 1 = 0.8δ VG case and x/δ = −1 for the H 2 = 0.2δ VG case. Three parameters were calculated following the Richardson Extrapolation Method: the extrapolated solution RE, the order of accuracy p and the error ratio R. Table 3 exhibits the results obtained in the mesh independency study where a monotonic convergence was fulfilled. The higher resolution of the fine mesh was used in the computations for the two vane-type VG cases analysed. The iterative solution process was carried out in the simulations of the two VG cases until the residual errors dropped below 10 −4 for the pressure p and 10 −5 for the velocities U x , U y and U z and the turbulence quantities used, e.g., the turbulent kinetic energy k and the specific dissipation rate ω. The values obtained for the dimensionless distance of the first wall cell for the meshes corresponding to the two VG cases is summarized in Table 4. In both cases the value of y+ is lower than 1 (y+ < 1) as required by the turbulence model adopted.  Figure 7 shows the vortex visualization based on the velocity field distribution at four planes normal to the streamwise direction and located at normalized distances x/δ respect to the boundary layer thickness δ from the backward-facing ramp where the domain origin is placed. The scale in the y and z axis for the four snapshots in this figure is approximately 13 mm × 13 mm. Thus, the vortex can be clearly visualized for the H 1 = 0.8δ vane-type VG case as expected. Note that the larger the normalized distance x/δ is, the closer the vortex is to the VG trailing edge (TE). The four planes depicted are located between the selected cells representing the VG and the ramp.   Figure 8 shows again the vortex visualization based on the velocity field distribution for the H 2 = 0.2δ vane-type VG case at four planes normal to the main flow direction. Now, the scale for the four snapshots is approximately 12.2 mm × 12.2 mm. The formation of the vortex is slightly observed but not as clear as it was observed in the previous H 1 = 0.8δ vane-type VG case shown in Figure 7. This seems to indicate that the vortex generated by the sub-layer boundary VG has a smaller strength than that corresponding to the conventional and higher VG. For instance, the vortex shape at plane x/δ = +0.35 is not as easy to observe in the figure as the vortex shape is for the H 1 = 0.8δ vane-type VG case.

Pressure Distribution
The effect generated by the H 1 = 0.8δ vane-type VG placed at x = −5δ upstream the ramp on the streamwise pressure coefficient distribution c p along the measurement lines A, B and C (Figure 2b) on the flat plate floor is appreciated in Figure 9. The three curves correspond to the pressure distribution (jBAY model, red rectangles) obtained from the CFD simulation with the jBAY Source Term Model implemented, the distribution data (Baseline, blue circles) obtained from the CFD simulation for a baseline case with no vortex generator positioned and the experimental pressure distribution data (EXP, green triangles) according to Lin et al. [33]. The error between the modelled (jBAY model) and the experimental (EXP) data is also shown in black vertical lines. The reader should note that the center of the backward-facing ramp is placed at the normalized distance x/δ = 0 and the VG is positioned in a negative value of x/δ.  According to Lin et al. [33] and Konig et al. [34], the pressure distribution variations and the two peaks observed in the three plots at the beginning and end of the ramp correspond to accelerations and decelerations as expected. The influence of the VGs on the pressure distribution is meaningful when comparing with the pressure distribution for the baseline case. The results seem to be in agreement with the effects corresponding with those generated by vane-type VGs which considerably increase the pressure gradient between the beginning and the end of the ramp. As expected, the vortex generator is going to bring about variations in the flow structures. When observing the simulated pressure distribution along the measurement line A (Figure 9a), the initial suction peak at the beginning of the ramp is slightly overestimated and underestimated at the end of the ramp. The gradient of the pressure is quite well matched. Small deviations are appreciated between the experimental (EXP) and the simulated pressure distribution (jBAY model) before and behind the pressure peaks which could be related to the evidence that those regions are placed well within the three-dimensional flow structures induced by the VGs as indicated by Konig et al. [34]. When observing the pressure distribution along the measurement line B (Figure 9b) defined as the floor line starting from the VG trailing edge (TE) and parallel to the main flow (Figure 2b), the initial suction peak at the beginning of the ramp is underestimated although the gradient of the pressure as well as the pressure maximum at the end of the ramp are quite well matched. Similarly, Figure 9c shows that the simulated pressure distribution (jBAY model) along the measurement line C has a smaller magnitude between the two peaks at the beginning and end of the ramp than the amplitudes obtained along lines A and B. This distribution is shifted towards positive distances downstream the ramp. Figure 10 depicts the effect generated by the H 2 = 0.2δ vane-type VG placed at x = −2.5δ upstream the ramp on the streamwise pressure distribution c p along the measurement lines A, B and C on the flat plate floor. The error (black vertical lines) between the modelled (jBAY model) and the experimental (EXP) data is again indicated. The amplitude between the peaks is smaller than the amplitude obtained in the previous VG case of height H 1 . On the other hand, the simulated pressure distribution (jBAY model, red rectangles) over the ramp as well as at the end of the ramp and for the three plots (Figure 10a-c) is slightly overestimated when comparing with the experimental distribution (EXP, green triangles).
In summary, the comparison of the simulated pressure distribution (jBAY model) with the experiment pressure distribution (EXP) for the H 1 = 0.8δ and H 2 = 0.2δ vane-type VG cases depicted in the two previous figures shows a relatively good agreement although relatively small deviations are appreciated, mostly for the H 2 = 0.2δ vane-type VG case. Thus, the mean absolute percentage error between the modelled and experimental data in the H 1 VG case is 7.57%, 8.10% and 7.10% for the three pressure distributions respectively measured in A, B and C lines. In the H 2 VG case, the mean absolute percentage error is 9.83%, 9.58% and 9.76%. The reason of these little bit larger deviations could be related to the fact that the VG height H 2 is within the sub-buffer zone of the local boundary layer where the shear and viscous stresses are predominant.

Vortex Path
The path or trajectory of the vortex in the vertical (y) and lateral (z) directions can be determined by computing the location of the vortex center generated by the VG all along the downstream axis x. The vortex center can be analytically defined as the point in a cross-stream plane which has the maximum value of the vorticity, the so-called peak vorticity ω max . In Figure 11, a comparison between the vortex paths corresponding to the H 1 = 0.8δ and H 2 = 0.2δ vane-type VG cases is indicated. Streamwise, lateral and vertical coordinates are normalized by the local boundary layer thickness δ to show the effects of VG scaling downstream of the vane.
The vertical path corresponding to the conventional H 1 = 0.8δ vane-type VG case represented in Figure 11a tends to be parallel to the flat plate and the flow direction far from the backward-facing ramp placed at x/δ ∼ =0, but when the vortex approaches the ramp, it starts moving downward probably due to the adverse gradient of pressure created by the ramp. As observed in the same figure for the sub-boundary layer H 2 = 0.2δ vane-type VG case, the vertical path has a fast downward deviation probably due to the adverse gradient of pressure created by the ramp and its proximity to the ramp.  The vortex paths for the H 1 = 0.8δ and H 2 = 0.2δ vane-type VG cases in the lateral direction are represented in Figure 11b. Both lateral trajectories seem to show a relatively small increasing trend in the same direction in which the vanes are pointed. As expected, the deviation of the lateral trajectory of the primary vortex described in the spanwise direction (z) by the highest VG is larger than the lateral deviation observed in the lowest VG. The reason for this different trajectory would lay on the fact that the lowest VG is within sub-buffer region of the boundary layer and consequently its effect on the oncoming flow is less significant. The starting point for both cases again depends on the location of the vane on the flat plate. This lateral increasing trend for the two trajectories is almost the same for both cases although for the H 2 = 0.2δ vane-type VG case the movement occurs in a lower streamwise distance. Thus, the mean positive slopes of the lateral trajectory for the H 1 and H 2 cases are around 0.005 and 0.006 respectively.

Vortex Decay
In order to illustrate how the vortex decays, the non-dimensional streamwise distribution of the normalized peak vorticity (ω x max δ)/U ∞ can be studied. In Figure 12 the non-dimensional streamwise distribution of the normalized peak vorticity is plotted function of the normalized downstream distance x/δ for the two vane-type cases. The figure depicts the vortex decay in the streamwise direction and ratifies that the peak vorticity ω x max quite fast weakens past the VG for the two cases analysed. As observed, the maximum of the streamwise peak vorticity is obtained downstream the vane at the position x/δ = −2.8 and x/δ = −1.8 for the H 1 = 0.8δ and H 2 = 0.2δ VG cases respectively. After the maximum, the peak vorticity seems to decay in a decreasing exponential way to x/δ which is in concordance with Fernández-Gamiz et al. [39]. As expected, the peak vorticity magnitude decreases when the VG height is reduced. Figure 12. Non-dimensional streamwise peak vorticity (ω x max δ)/U ∞ function of the normalized distance x/δ for the H 1 = 0.8δ (blue rectangles) and H 2 = 0.2δ (red circles) vane-type VG cases. VGs, flat plate and ramp also shown (not at scale).

Vortex Size
The vortex size is a key parameter when modelling vortices. A suitable way to study the vortex size is by means of the half-life radius R 0.5 defined as the radial distance from the vortex center to the point where the vorticity is half the peak vorticity ω max captured in a cross-stream plane according to Bray [40]. At that point, the measurement errors when calculating the vortex size are negligible and the accuracy degree is high. Figure 13 shows the vortex size evolution expressed in terms of the normalized half-life radius R 0.5 for the H 1 = 0.8δ and H 2 = 0.2δ vane-type VG cases at different locations x/δ past the VG. This figure shows a increasing trend in the vortex size when moving downstream for both VG cases. This relationship observed between both magnitudes could be formulated by a mathematical expression and this could eventually allow establishing a prediction model between the vortex size and the streamwise distance. As expected, the vortex size generated by the conventional VG of height H 1 = 0.8δ is larger than the vortex generated by the sub-boundary layer VG of height H 2 = 0.2δ. Figure 13. Vortex size evolution expressed in terms of the normalized half-life radius R 0.5 at difference normalized distances x/δ for the H 1 = 0.8δ (blue rectangles) and H 2 = 0.2δ (red circles) vane-type VG cases. Figure 14 shows the wall shear stress distribution along the measurement line B on the flat plate floor for the two vane-type VG cases analysed in the present study. The wall shear stress distribution obtained from the implementation of the jBAY Source Term Model (red rectangles) is compared with that corresponding to the CFD simulation for the baseline case (blue circles). According to Godard and Stanislas [41], the implementation of VGs leads to an increase of the values of the wall shear stress downstream of the vanes in comparison with the cases where passive devices are not implemented. The results show a larger increment between the jBAY and baseline distributions for the H 1 = 0.8δ vane-type VG case than for the H 2 = 0.2δ vane-type VG case. This higher increase could be related to the fact that in the case of the conventional VG of height H 1 = 0.8δ the streamwise vortex induced is stronger and consequently its effects on the wall shear stress bigger than in the case of the sub-boundary layer VG of height H 2 = 0.2δ whose vortex induced is weaker and so are its effects, for instance, on the wall shear stress. However, these results obtained for the jBAY modelled wall shear stress seem to be again in accordance with those generated by rectangular VGs located on a flat plate which increase the wall shear stress downstream of the VG. According to Godard and Stanislas [41], this increase tends to delay the boundary layer detachment. Two wall shear stress peaks are appreciated at the beginning and end of the ramp which indicate the influence of the ramp on the shear stress. As observed, these two peaks are larger for the H 1 = 0.8δ vane-type VG case than for the H 2 = 0.2δ vane-type VG case.

Conclusions
The generation of vortices and their effects by a conventional vortex generator and a sub-boundary layer vortex generator positioned on a three-dimensional flat plate with a backward-facing ramp and adverse gradient pressure has been carried out by means of CFD simulations using the open-source code OpenFOAM. The influence of these two vane-type vortex generators (VGs) on the computational domain flow is implemented by using a source term in the corresponding Navier-Stokes equations according to the so-called jBAY Source Term Model. Steady-state, incompressible and turbulent flow is assumed and Reynolds Average Navier-Stokes (RANS) turbulence modelling is applied in the simulations.
The Source Term Model seems to simulate relatively well the streamwise pressure distributions along the floor of the flat plate for the two vane-type vortex generators studied of heights H 1 = 80% and H 2 = 20% the local boundary layer thickness δ at the upstream edge of the ramp. The results obtained for the jBAY modelled pressure distributions are in concordance with the experiments where the influence of these VG devices is measured. However, the streamwise pressure distribution has been slightly underestimated in certain regions far from the ramp and overestimated in other regions by the jBAY Model.
The jBAY modelled pressure distributions for the conventional vortex generator seem to be more accurate than those obtained for the sub-boundary layer vortex generator. The reason of this fact could be related to the height of the sub-boundary layer vortex generator which is within the sub-buffer zone of the local boundary layer where the shear and viscous stresses are predominant. Consequently, the results could potentially indicate that the jBAY Source Term Model presents more difficulties when modelling sub-boundary layer vortex generators.
Other parameters for vortex characterization analysed in the present study such as vortex path, vortex decay and vortex size obtained by means of the jBAY Source Term modelling, are in relatively good agreement with the results expected.
The implementation of the Source Term Model can represent an advantage over a fully Mesh-Resolved Vortex Generator Model for certain application cases due to a meaningful decrease in the cell number of the computational domain with the corresponding saving of computational time and resources. In the present work where the Source Term Model has been applied in the simulations of a conventional and a sub-boundary layer vortex generator of respective heights H 1 and H 2 , the total cell number of the mesh could be eventually small when compared with the conventionally used fully Mesh-Resolved Vortex Generator Model. In addition, the jBAY Source Term Model could be quite helpful to reduce the meshing time.
Further research should be done in order to model and characterize vortex generators (VGs) of different dimensions and geometries positioned in surfaces of interest such as high-lift airfoils or ducts by implementing the Source Term Model or an optimized version of the Source Term Model into OpenFOAM or other CFD codes.