An integral boundary layer engineering model for vortex generators implemented in XFOIL

To assess and optimize vortex generators (VGs) for flow separation control, the effect of these devices should be modelled in a cost and time efficient way. Therefore, it is of interest to extend integral boundary layer models to analyse the effect of VGs on airfoil performance. In this work, theturbulentboundarylayerformulationismodifiedusingasourcetermapproach.Anadditional term is added to the shear-lag equation, to account for the increased dissipation due to streamwise vortex action in the boundary layer, forcing transition at the VG leading edge where appli-cable. The source term is calibrated and a semi-empirical relation is set up and implemented in XFOIL.ThemodifiedcodeiscapableofaddressingtheeffectoftheVGheight,length,inflowangle, and chordwise position on the airfoil's aerodynamic properties. The predicted polars for airfoils with VGs show a good agreement with reference data, and the code robustness is demonstrated by assessing different airfoil families at a wide range of Reynolds numbers.


INTRODUCTION
For the next generation of wind turbines, manufacturers aim to design multimegawatt rotors to improve the competitiveness of wind energy technology. To up-scale wind turbines, novel technologies are required and new design challenges will appear. One such aerodynamic challenge is the management of separated flow. Preventing or at least delaying separation over the blades can positively affect the annual energy production (AEP).
On top of that, the magnitude and severe variations of the aerodynamic loads associated with separating flows can be reduced, mitigating structural fatigue issues. In the wind energy industry, separation control is often realized by using passive vortex generators (VGs).
Vortex generators improve the resistance of a boundary layer against flow separation by re-energizing the flow close to the surface. The streamwise vortices shed from the free tips of the VGs enhance mixing between the high-energy flow in the outer part of the boundary layer with the low-energy regions near the walls (see Schubauer and Spangenberg. 1 ) The physics involved is nontrivial and poses a number of modelling challenges.
To achieve a cost-effective scale up of current turbines, it will become necessary to evolve towards a multidisciplinary design process where VGs are already incorporated early in the design phase. To assess and optimize the use of VGs, there will be a need to effectively model these devices in a cost and time efficient manner.

Background
Numerous VG modelling techniques have been explored in literature, most of which use computational fluid dynamics (CFD). The most direct and intuitive approach is to model the effect of VGs by including them as a local geometrical protrusion in the domain. This approach requires a fully 3-dimensional treatment and a dense grid due to the complexity of the flow and the small size of the VGs (see Nikolaou and Politis, 2 Jirasek, 3 and von Stillfried et al 4 ).
In 1999, Bender et al 5 developed a simplified technique in which the effect of VGs was modelled without including the geometry explicitly in the computational mesh. This approach is now commonly known as the BAY model. The method relies on the addition of a side force as a source term in the momentum equation. This essentially introduces and simulates swirl in the flow, giving rise to the formation of vortices produced by the VGs. Jirasek 3,6 suggested further improvements to the original model by refining the definition of the source term. Besides adding a source term based on a lateral lift force, other researchers constructed source terms based on a predetermined vortex circulation, or on the velocity induced by this circulation according to the Biot-Savart law. This method is referred to as a vortex-source model. Research based on this method is found in the spanwise averaged CFD method of Nikolaou and Politis 2 and the 3-dimensional approach by Zhang. 7 As Jirasek 3 noted, the key problem of these kind of models is that the initial circulation introduced by the VG needs to be known beforehand. Later, Tornblom and Johansson 8 presented a method in which a Reynolds stress approach is used in a statistical sense. They explicitly introduced the vortex-added velocity and associated stresses into the differential Reynolds stress transport equations to mimic the increased mixing due to the VGs. 8,9 Despite considerable advances with CFD-based methods, such techniques are generally more time and resource consuming, compared with the industry workhorses that are integral boundary layer (IBL) tools and blade element momentum methods. This makes them impractical for use in iterative design processes. Therefore, it is of interest to extend integral boundary layer models to analyse the effect of flow mixing devices on airfoil performance. Limited research is done in such a modelling approach incorporating VGs, and it is useful to consider the physics involved from an integral boundary layer perspective. Drela 10 argues that VGs promote increased dissipation by introducing streamwise vortices into the boundary layer, which consequently increases the sustainable adverse pressure gradient. Lengani et al 11 demonstrated this by analysing the dissipation mechanisms in a separating duct controlled by vane-type vortex generators. Dissipation in this sense refers to the drain of mean flow kinetic energy through the action of the shear stress with the mean strain rate, essentially, through enhanced mixing.
In essence, these observations are captured in the IBL method proposed by Michael Kerho and Brian Kramer. 12 In their approach, the VG-added mixing effect was incorporated into the XFOIL code through an enhancement of the turbulence production in the turbulent boundary layer formulation by modifying the stress transport formulation. The authors suggested modifying the boundary layer dissipation through the shear stress coefficient by introducing a stepwise increase in the production term, which decreases exponentially downstream. An obvious limitation of this approach is that the fully 3-dimensional flow field induced by an array of vortex generators is represented and modelled by a 2-dimensional integral boundary layer formulation. However, it is possible for 2D design codes to handle these complex situations and offer a global or spanwise-averaged effect of VGs for design direction.

Research objective
Kerho and Kramer's 12 source term approach was proposed for a specific case with a single row of co-rotating VGs. Their focus was on natural-laminar-flow airfoils at relatively low operational Reynolds numbers with VGs installed far downstream for trailing edge separation control.
Therefore, the objective of this article is to model the effect of VGs in an IBL code (ie, XFOIL) by extending the source term implementation and overall code robustness for thicker airfoils over a wider operational envelope. This objective will be approached by the following: 1. Modifying the closure relations using a source term approach.
2. Setting up a relation between the source term, VG geometry, and boundary layer properties, using CFD and experimental airfoil reference data.
3. Verifying and validating the implementation and assessing the code robustness.
The article proceeds in Section 2 with a description of the methodology, modelling rationale, and boundary layer modifications in the context of vortex generators. Validation results are presented in Section 3 through comparisons with flat plate experiments and airfoil measurements. Further discussion is presented in Section 4 and concluding remarks in Section 5.

Integral boundary layer code
The integral boundary layer code XFOIL is a viscous-inviscid interaction method designed for predicting airfoil flows and performance. 13 In this framework, the flow is decomposed into 2 regions: the inviscid outer flow where viscosity can be neglected and the thin, viscous shear layer, that is, the inner flow, where the boundary layer plays an important role.
The outer flow is solved using a linear-vorticity streamfunction panel method. A 2-dimensional inviscid airfoil flowfield is composed of a free stream flow, a vortex sheet of strength on the airfoil surface, and a source sheet of strength on the airfoil surface and in the trailing wake. 13 The system is closed with the Kutta condition applied at the trailing edge.
The viscous boundary layer solution is obtained using the so-called numerical integral method. The partial differential equations for continuity, momentum and energy can be reworked and transformed into the familiar ordinary differential equations in terms of integral quantities: This method therefore determines the integral thickness and key shear quantities, namely, displacement thickness ⋆ , momentum thickness , friction coefficient C f , and the dissipation coefficient C D . 10 The integral momentum and kinetic energy equations are consequently combined with a chain of laminar and turbulent closure relations in order to make the problem determinate. 14 Finally, both solutions are coupled using a fully simultaneous coupling scheme described by Drela et al. 15 The entire nonlinear equation set is solved simultaneously as a fully coupled system using a global Newton-Raphson method.

Boundary layer formulation modification
The effect of VGs can be introduced in a number of ways. In this case, the boundary layer formulation is modified through the dissipation term. To do so, a source term approach is used, in which an additional production term is added to the turbulence closure relations. The choice of closure relation should be physically consistent, but also facilitate implementation, without compromising the code convergence behaviour. Thus, 2 sub-goals were to determine (1) in which equation to incorporate the source term and (2) to explore the behaviour of this source term over the streamwise domain.

Source term implementation
The source term is implemented into the rate equation (or shear-lag equation), an ordinary differential equation for the shear stress level inside the boundary layer. The equation models turbulence history effects that dominate turbulence intensity (see Drela and Giles 16 ) Additionally, the rate equation depends on an empirical constant K c , controlling the reactivity of the boundary layer, and on the equilibrium shear stress coefficient C EQ . The latter represents the shear stress level that would exist if the local boundary layer would be in equilibrium, in this sense, meaning that the boundary layer profile for a turbulent flow exhibits a behaviour analogous to the similar flows of laminar boundary layers. 17 In slowly changing flows, the flow is almost in equilibrium and thus C closely follows C EQ . However, this does not hold for rapidly changing flows, and thus, the rate equation plays an important role in modelling the lag effect. A case in point are vortical flows, in which the Reynolds stresses are known to lag the mean strain field. 18 A simple expression for C EQ can be obtained from the well-known G-locus of equilibrium boundary layers proposed by Clauser. 16 The source term implementation is presented in Equation 3 where S VG represents the source term itself: Because C EQ is seen as the fundamental parameter of the rate equation, the VG source term is added explicitly to the equilibrium shear stress coefficient C EQ , effectively enhancing the production term of the turbulent kinetic energy budget. In this sense, the non-equilibrium nature of VG flow increases the tendency of the boundary layer to depart from its equilibrium state.
This implementation was also favourable considering the convergence speed and numerical stability of XFOIL. Qualitatively speaking, these modifications promote an increased momentum thickness but reduce the displacement thickness. The shape factor decreases and the skin friction coefficient will increase. 10,11 It is important to note that these trends may be skewed by the physical presence of the devices, which according to Schubauer and Spangenberg, 1 cause local jumps in * and that are amplified with the developing boundary layer. This will become evident in the following section.

Source term shape function
The source term shape function is selected so as to mimic the streamwise vortex strength decay from mixing devices (see, eg, Lögdberg et al 19 and Baldacchino et al 20 ). Such studies have shown that vortex strength decay is roughly exponential and arises because of the presence of wall shear, boundary layer turbulence and mutual interference between adjacent vortices. On the other hand, vortices shed by VGs do not appear abruptly, and generally require around one vane chord length to develop before reaching their full rolled-up intensity. Therefore, the selected source term shape is given by Equation 4 where the source term develops gradually before decaying exponentially.
As shown in Figure 1, the shape function differs from the expression used by Kerho and Kramer, 12 who proposed instead a simple step function.
The new source term expression consists of 3 variables: the source term strength 0 , decay rate , and the location of the VG with respect to the leading edge. Assuming that the source term is related to the strength of the shed circulation, it is expected that the source term strength 0 and decay rate will be a function of the VG configuration and the local boundary layer properties.

Laminar-turbulent transition
XFOIL uses the e N method for predicting natural transition. 21 This method assumes that transition occurs when the most unstable Tollmien-Schlichting wave in the boundary layer has grown by a given factor e N , where N = 9 for natural transition. 13 Vortex generators that are typically sized on the order of the boundary layer thickness will likely lead to bypass-type transition within a short region. Research has shown that small VGs can actually delay transition, by attenuating critical perturbations in the boundary layer. 22 However, this is out of scope of the applications envisaged in this work. A simplified transition definition is used that assumes that VGs promote flow transition at their leading edge, independent of their configuration: where x tr represents the transition location with respect to the airfoil leading edge. The validity of this assumption will depend somewhat on the VG configuration and the local pressure gradient. However, in many practical applications, VGs are mounted on a strip that itself can enforce transition, as illustrated in Figure 2.

Source term calibration
Of the 3 source term shape function variables, the source strength 0 and the decay rate are unknown and expected to depend on the airfoil/VG configuration. The required value of both parameters is determined here based on a calibration process using reference data obtained from high-fidelity numerical simulations and measurements. Data sets are adopted from studies of Timmer and van Rooij, 23 Baldacchino et al, 20,24 Manolesos and Voutsinas, 25 and the public database provided by the AVATAR project. [26][27][28] These cover a variety of VG configurations, airfoil families, and inflow conditions. A summary of the data sets is given in the Appendix. For different VG configurations in various flow conditions, the source term is determined in 3 steps: 1. Analyse the effect of the source term parameters on lift and drag.
2. Define the target aerodynamic properties.
3. Determine the required source term parameters.
This procedure is summarized in Figure 3. Throughout the remainder of this article, cases without and with vortex generators are referred to as the clean case and VG case.

Part 1-source term effect
The effect of different source terms with different strength and decay parameters is presented in Figure 4 for 2 different angles of attack. Increasing the source term strength and decay at low angles of attack (ie, below stall) results in a decreasing lift coefficient and increasing drag. However, at low angles, the sensitivity is small because the boundary layer is less receptive to the added dissipation. At high angles of attack, the opposite is true and the increased mixing gradually takes effect with increasing angle of attack.
One may ascertain that multiple combinations of source term strength and decay predict the same lift and drag coefficient. It has been established that all these combinations of strength 0 and decay constants share the same value of the source term integral, as shown in Figure 5. This source term integral I ST , is the area enclosed by the source term function: The integral is numerically determined based on the airfoil panelling using the trapezoidal rule. In keeping with the IBL modelling ethos, the definition of I ST reduces the number of unknowns (originally the strength 0 and decay ) to a single variable.

Part 2-target lift polar
To account for the discrepancies of the clean polar predicted by the airfoil code, corrected (target) polars expected to be found by the airfoil code are introduced. With reference to Steps 1 and 2 in Figure 3, the lift slope discrepancy between clean XFOIL predictions and reference data is first calculated. The slopes were evaluated in the interval = [0 • , 4 • ]. Finally, in Step 3, the slope correction factor is applied to the reference lift polar for the VG case. This procedure avoids the source term accounting for inherent inaccuracies of the integral boundary layer code. The corrected lift polar is referred to as the target lift curve and is thus thought to be a more representative reference with which to calibrate the empirical function.

Part 3-source term selection
The determination of the optimal source term is based primarily on the lift coefficient. The drag coefficient is not found suitable because integral boundary layer codes significantly underestimate drag in clean case and because the source term approach does not account for the additional profile drag of the VG itself. Boundary layer measurements with vortex generators providing reliable integral properties are sparsely available and thus also unsuitable for the present purpose.

Source term semi-empirical relation
To  (4) the VG height, h VG ∕l VG , scaled as the vane aspect ratio. These variables, schematically presented in Figure 6, have been previously identified as being directly correlated to the strength of the shed vortex (see Wendt and Reichert, 29 Angele and The reliability of such a calibration depends of course on the quality of the predictions in clean conditions. To account for this, the calibration procedure uses a weighting function taking into account the error in predicting the lift slope, maximum lift coefficient, and stall angle compared with the reference data. Additionally, focus is placed on higher angles of attack since at lower angles, the effect of the source term on lift and drag is smaller than the expected accuracy of the code. This process is also shown in Figure 3.
The VG parameters h ⋆ VG and l ⋆ VG represent the non-dimensional vortex generator height and length with respect to the airfoil chord. To retrieve the flow velocity at the VG tip (U VG ), the Swafford boundary layer velocity profile as used in XFOIL is reconstructed using the momentum thickness Reynolds number Re , the kinematic shape parameter H, the edge velocity U e , and the boundary layer thickness . 13,32 Figure 7 provides the error distribution between the calibrated source term, S VG | cal and the empirical fit, S VG | emp . The residuals are approximately normally distributed.
The source term expression (Equation (8)) is implemented in XFOIL in such a way as to minimize the required user interaction. Therefore, the only inputs required are the VG geometry: the VG height, length, inflow angle, and position with respect to the airfoil leading edge. The inflow velocity U VG is calculated internally using an iterative process that is consequently used to update the value of the source term according to the empirical relation.

RESULTS
The VG modelling capabilities are assessed in this section. Results are presented for a canonical case of a flat plate with and without vortex generators and subsequently for 3 airfoil sets controlled using vortex generators. The former allows a more straightforward assessment of integral boundary layer properties, whereas the latter provides basis for global comparison.

Boundary layer properties
The polynomial that was merged with the middle section using a square root blending function. The trailing 20% chord was modelled as the aft 70% of a NACA0010 symmetric airfoil (further details of this method are found in Sanders. 33 ) A uniform pressure distribution was obtained at = 0 • .
The predicted evolution of the boundary layer properties is shown in Figure 8 for the clean and VG case. The relatively short range of the measurements downstream limits the extent of the comparisons. However, the measured and predicted trends compare reasonably well. Notable discrepancies exist at the location of the VG itself: the measurements demonstrate somewhat of a step increase in * and . The vortex action gradually enhances the boundary layer state, seen as an eventual decrease in the shape factor ( Figure 8C). However, a passive device is itself a source of drag that manifests as a mass and momentum deficit in the developing boundary layer, over and above the effect of the developing vortex. Schubauer and Spangenberg 1 originally demonstrated this aft of passive devices in a boundary layer wind tunnel. They observed that the evolution of the relative momentum deficit remained rather constant. This implies that the device drag could potentially be modelled as a stepwise perturbation that would subsequently be convected and modulated with the developing boundary layer. Since the present model implementation does not account for this penalty drag, this additional deficit cannot be captured.

Global performance assessment
This part of the validation aims to (1) validate the source term expression and (2) demonstrate the code robustness. This will be achieved by validating the source term sensitivity with VG location, height, and inflow angle by comparison with data sets outside the reference database.

Source term expression
In Figure 9A  At low angles of attack, the closer the VGs are to the leading edge, the higher the drag. Within the current modelling framework, this occurs primarily because of the earlier forced transition to turbulence, minimizing the lower-drag laminar flow region over the airfoil. These observations are evident in the experimental data, which are reasonably well predicted by XFOIL. Figure 10 shows a similar comparison using the same reference data set, comparing the effect of the VG height. The Delta-shaped VGs are located at 40% chord from the airfoil leading edge and have a height of 5 and 10 mm, all other parameters dictated through geometric similarity. Increasing the VG height at this location improves their effectiveness in separation delay but also introduces more drag at pre-stall angles of attack.
The effect of the inflow angle is presented in Figure 11. The reference data originates from a numerical parametric study performed within the AVATAR project. 27 The synthetic data were generated using fully turbulent, geometry-resolving CFD computations on a 33% thick FFAW3333 airfoil with a chord length of 5.84 m. The use of synthetic data permitted a comparison for a high Reynolds number (in this case 14 × 10 6 ) for which experimental data with VGs is not readily available. The airfoil was equipped with 30-mm high and 90-mm long counter-rotating Delta-shaped VGs, mounted at 40% chord. Evidently, increasing the vane inflow angle from 15 • to 25 • increases the maximum lift. The drag at low angles of attack increases for higher vane angles and decreases at high angles of attack as separation is delayed. Note that while the drag increases at a similar rate after maximum lift is attained, the numerical computations predict an earlier merging point than XFOIL. This is the angle at which the VGs, in all cases considered, are rendered completely ineffective because of large scale separation, resulting in a collapse of all polars.

Code robustness
Additional comparisons are made with the measurements reported in Fouatih et al. 34 Figure 12A to 12C presents the sensitivity of the relative change in maximum lift (ΔC l,max ) to the VG height h VG ∕c, VG location x VG ∕c, and the inflow angle VG . Except for the VG location, the trends are generally captured well, along with the order of magnitude.

Source term expression
The regression coefficients of the expression (Equation (8)) reveal that the source term integral increases with increasing VG length, VG height, inflow angle, and flow velocity. The capability of the airfoil code to model VGs significantly depends on the quality of the source term expression.
This in turn is limited by the quality of the calibration database, and how well represented the key variables are in those data sets. These limitations can be summarized as follows: • Error uncontrolled airfoil: For the uncontrolled performance (without VGs), XFOIL mostly overestimates the lift slope as well as the maximum lift coefficient and underestimates the drag. Despite efforts to anticipate for this using correct target polars, the calibrated source term will to some extent also account for the discrepancies in uncontrolled airfoil predictions.
• Reference data: The performance of the code is limited by the availability, consistency, and quality of the reference data. Not all VG parameters (eg, pair spacing and VG shape) were sufficiently represented in the database. Improving this representation is seen as a potential continuous improvement of the presented method.
• Boundary layer properties: The velocity at the VG is calculated using the edge velocity, integral boundary layer properties, and expressions for the boundary layer thickness and velocity profile. This includes additional uncertainties and affects the quality of the source term expression.

Validation
Comparing all available reference cases with the XFOIL predictions, we can establish goodness measures, as shown in Figure 13A to13C. These show the error distribution for the maximum lift coefficient, the stall angle, and the drag coefficient at zero angle of attack. The error is defined as the difference between the reference and predicted XFOIL values, normalized with respect to the reference value. The following observations can be made: • Maximum lift coefficient: The error for the maximum lift coefficient is normally distributed with a mean value near zero. This result was expected since the composition of the source term expression was based on the maximum lift coefficient. With a 90% confidence interval, the error of the predicted maximum lift coefficient with respect to the reference value will be within ±12%.
(A) (B) (C) FIGURE 13 Error distribution between the reference data and the XFOIL calculations [Colour figure can be viewed at wileyonlinelibrary.com] • Stall angle: While the modal value is roughly zero centred, the mean error is approximately −10%, meaning that XFOIL mostly overpredicts the stall angle. However, this is a tendency of the present code for predictions even without VGs.
• Drag coefficient: The drag coefficient is consistently underpredicted. The drag at zero angle of attack found by XFOIL is more than 25% lower than the reference value. The net increment in drag by adding VGs is lower than given by the reference data since the parasitic VG drag has not been accounted for in the present model.
The methodology presented demonstrates a potentially useful engineering approach for modelling vortex generators in integral boundary layer codes. Reasonable predictions have been demonstrated. Some of the limitations and assumptions discussed prior can be better managed with the following suggestions: 1. Reference data: An expanded reference database to improve the representation of important variables within the calibration process.

CONCLUSION
This article demonstrates a method for modelling the effect of VGs using the integral boundary layer code XFOIL. The methodology builds on the source term approach introduced by Kerho and Kramer 12 . This was realized by adding an additional term to the equilibrium shear stress coefficient of the shear-lag equation, accounting for the increased dissipation due to streamwise vortex action in the boundary layer. A gradual step input followed by exponential decay was introduced at the VG location, mimicking the inception and evolution of an embedded streamwise vortex. This perturbation was described using 3 variables: the source term strength, the decay rate, and the location of the VGs.
Synthesizing these parameters, the source term integral was defined and correlated with the VG geometry and inflow conditions. This was set up using a least-square regression to calibrate the expression parameters with an extensive reference database. The resulting expression gives the required value of the source term integral as a function of the VG height, length, inflow angle, and flow velocity at the tip of the VG.
The modified XFOIL code can predict the aerodynamic behaviour of an airfoil equipped with VGs. The new code is capable of addressing the effect of the VG height, length, inflow angle, and chordwise position but is not expected to capture subtle differences in VG geometry. The error on the maximum lift coefficient is normally distributed with a mean value near zero. The stall angle is mostly overpredicted, but this is in line with the clean airfoil predictions. The drag at zero angle of attack found by XFOIL is more than 25% lower than the reference value, which can be explained because of the lack in modelling the parasitic drag of the VGs itself.
Overall, the predicted polars for airfoils with VGs showed a good agreement with the reference data and the code is demonstrated to be robust and able to model different airfoil families at a wide range of Reynolds numbers. The source term approach is proven to be promising and can be elaborated further considering some of the suggested improvements.