Extended Glauert Tip Correction to Include Vortex Rollup Effects

Wind turbine loads predictions by blade-element momentum theory using the standard tip-loss correction have been shown to over-predict loading near the blade tip in comparison to experimental data. This over-prediction is theorized to be due to the assumption of light rotor loading, inherent in the standard tip-loss correction model of Glauert. A higher- order free-wake method, WindDVE, is used to compute the rollup process of the trailing vortex sheets downstream of wind turbine blades. Results obtained serve an exact correction function to the Glauert tip correction used in blade-element momentum methods. It is found that accounting for the effects of tip vortex rollup within the Glauert tip correction indeed results in improved prediction of blade tip loads computed by blade-element momentum methods.


Introduction
Though significant advances have been made in high-fidelity computational fluid dynamics (CFD), wind turbine design and analysis continues to rely heavily on efficient blade-element momentum theory (BEMT) methods; however, the wind energy community has been challenged with a persistent overprediction of blade tip loads computed by BEMT, thus adding uncertainty to the anticipated performance of full-scale wind turbine rotors. Reference [1] includes a complete discussion of the state-of-the-art in BEMT tip modeling.
The inherent difficulty of BEMT in predicting blade tip loads is fundamentally rooted in its striptheory assumption, with the attempt of including all relevant three-dimensional effects into a single Glauert tip correction as used in most of today's BEMT methods. One has to be aware, though, that the Glauert tip correction [2] is a modified version of Prandtl's earlier suggestion of a correction factor obtained by replacing the helicoidal vortex sheet downstream of rotating wings with an infinite series of planar and rigid vortex sheets. Hence the Glauert tip correction does not account for tip vortex rollup and wake distortion, and its applicability is strictly only justified for lightly-loaded rotors, a limitation that stays in contrast to modern highly-loaded wind turbines.
Several attempts have been made to finding either a new tip loss factor or extending the classical Glauert tip correction. For example, the work of Shen et al. [3] suggested a correction function applied to sectional normal and tangential force coefficients in the tip region. More recently, Sørensen suggested a refined tip correction based on the decambering effect experienced by tip airfoils as a result of flow induction from wake vortices [4]. A highly innovative approach has been followed by Branlard et al. [5,6] who conducted a systematic parameter sweep of vortex wake solutions, with the intent to populate a database of blade circulation distributions that can be interrogated during a BEMT solution process. The challenge with all approaches to improve BEMT tip load predictions is that results obtained have to be consistent with variations in blade pitch angle, tip-speed ratio, blade tip shape, properties of tip airfoils, etc., without adding notable computational overhead to BEMT.
The basic idea being pursued in the current work is to support the hypothesis that vortex rollup is the primary cause of the observed overprediction of blade tip loads by today's BEMT methods. The paper is organized as follows: A mathematical theory is presented in section 2 that allows to extract analytically vortex rollup effects from a higher-fidelity computational method by means of an additional function added to the Glauert tip correction [7,8]. Section 3 introduces the higher-fidelity computational method employed in this work, the higher-order free-wake method WindDVE [9], and demonstrates that a BEMT code, XTurb [10], can be easily modified to obtain the exact same tip loads as the higher-fidelity method. Section 4 presents variations in blade pitch angle, tip-speed ratio, blade tip shape, and tip airfoils that are valuable data towards generalizing the method for a simple extension to the Glauert tip correction.

Analytical Tip Loss Factor from BEMT Equations for Given Blade Flow Angle, Φ
The ultimate goal of this work is to develop a simple addition to the standard Glauert correction that accounts for tip-vortex rollup effects on highly-loaded wind turbine blades. The approach is to introduce a function g(r) into the Glauert correction as with g(r) < 1 in the outer blade region to account for tip vortex rollup, and g = 1 to reduce to the classical Glauert correction. In Equation (1), B is the blade number and Φ is the blade flow angle. Functional forms of g(r) are computed analytically from a higher-order free-wake method and revisited BEMT equations [7,8]. The fundamental question arises as to how g(r) can be obtained from a higher-fidelity solution, for example a free-wake method capable of computing the rollup process of trailing wake vortex sheets. The difficulty lies in the fact that a suitable g(r) is expected to result in essentially the same blade tip loading, though computed by an iterative BEMT solution process. The theoretical solution to this problem can be obtained by first noting that the blade flow angle, Φ, determines a unique BEMT solution state [11]. Given Φ from a higher-fidelity solution, the BEMT equations can be inverted to solve exactly for the corresponding tip loss factor, F(r), and hence g(r). The detailed analysis is given below, with separate consideration of the 'windmill state' (a < 0.4) and the 'turbulent wake state' (a > 0.4) for highly-loaded blade elements.

Windmill State, a ≤ 0.4
The BEMT relations for sectional thrust, dT, and sectional torque, dQ, are given below:

Momentum Theory
Blade-Element Theory (2) Note that V0 refers to the incoming wind speed, Vrel to the relative velocity magnitude at a section of blade, c to the local chord length, Cl and Cd to sectional lift and drag coefficients, B to the blade number, ρ to the fluid density, and Ω to the rotor angular speed. Instead of solving for a and a', the system of BEMT Equations in (2) can be written with Vrel/V0 and F being the unknowns for a given blade flow angle, Φ, as: = Ω / 0 is the local tip-speed ratio, 1 = ′ ( cos( ) + sin( )), and 2 = ′ ( cos( ) − s( )), with ′ = /(2 ) being the local blade solidity. After some algebra, equation set (3) can be solved analytically for / 0 to obtain

Turbulent Wake State, a > 0.4
In the turbulent wake state, the momentum-theory relation for the incremental thrust, dT, is altered by an empirical relation obtained from comparison with data taken by Glauert. A number of relations are in use, among which are those documented by Wilson [12], Burton [13], and Buhl [14]. To simplify the mathematics, a linear model is used, similar to [12,13], resulting in

Analytical Solutions for ( ) and ( )
For all axial induction factors, a, the second relation in either equation sets (3) or (6) can be solved for the total loss factor, F, where / 0 is either determined by Eqns. (4) or (8) depending on the local axial induction factor, a. The total loss factor, F, can be written as the product F = FT *. FR, with FR being a standard root loss factor and FT * being the adjusted tip loss factor in Eqn. (1). The correction function g(r) is then obtained by inverting Eqn. (1) as: For a given higher-fidelity solution capable of resolving vortex rollup, e.g. a free-wake solution, the mathematical theory provided above allows to solve analytically for a correction function ( ) that can be implemented into a BEMT solution algorithm through Eqn. (1), thus implicitly enabling BEMT to account for vortex rollup.

Numerical Methods
This section introduces the free-wake method, WindDVE, used as the higher-fidelity method to compute vortex rollup and a series of g functions according to the theory derived in section 2. A brief introduction is given to the BEMT code, XTurb, followed by a proof-of-concept study demonstrating that a simple extension can be added to the Glauert tip correction such that XTurb can recover the WindDVE solution in the blade tip region.

Higher-Order Free-Wake Method, WindDVE
The WindDVE code is a free-wake method specifically designed to eliminate vortex-filament singularities through the use of distributed-vorticity elements (DVEs) [9,[15][16], thus enabling the computation of vortex and wake rollup, see Fig. 1. Using the same airfoil tables as a respective BEMT analysis, the WindDVE code can then be used to compute functions g(r) through Eqn. (10) and the blade flow angle, Φ. The blade flow angle was computed in WindDVE through an equivalent angle-of-attack, which is based on the blade circulation. The circulation across a lifting surface vortex panel is found by assuming a parabolic spanwise distribution of circulation across each blade element, with the circulation equal to zero at the blade tip and by enforcing that the velocity normal to a blade element's control point is zero [16]. The lift of a blade element is then computed by summing the analytic spanwise integration of the freestream kinematic force with the numerical integration of the induced force as computed at the element trailing edge [9]. Finally, the equivalent angle of attack of the zero-lift chord line is computed by dividing the lift coefficient by the theoretical thin-airfoil theory lift curve slope of 2π per radian, as typically done in potential-flow based methods.

BEMT Analysis, XTurb
The XTurb code is a blade design and analysis code developed at The Pennsylvania State University [10]. The code can be run in two primary modes, i.e. a classical BEMT analysis or using a prescribed wake method [17,18]. Furthermore, the XTurb code employs automatic airfoil data correction for rotational augmentation according to Du and Selig or a novel concept of solution-based stall delay (SBSD) as proposed by Dowler and Schmitz [19]. In this work, XTurb is used exclusively in BEMT mode and using the SBSD model to account for rotational augmentation.

Proof-Of-Concept -Extended Glauert Tip Correction for BEMT Analysis
An initial proof-of-concept study was conducted to verify the theory presented in section 2 by demonstrating that a simple (approximate) g function can be implemented into a BEMT analysis code, XTurb, through Eqn. (1). A more detailed study for the NREL Phase VI rotor can be found in Reference [7,8]. Figure 2a shows a tip-speed ratio (TSR) sweep for the J-Sequence of the NREL Phase VI rotor as an example of WindDVE computed g functions. In general, one expects g < 1 in the tip region in the presence of vortex rollup effects. This can be seen indeed in Fig. 2a; however, it is interesting to note that depending on TSR, the g function may actually exceed a value of one and reduce to close-to zero at the blade root. This behavior is attributed to the interplay between the Glauert tip correction and root loss factor. The primary observation from Fig. 2a is that there is an approximately linear decrease in g in the tip region. In addition to that, the value of g at the blade tip, g Tip, is practically independent of the tip-speed ratio (TSR), while the spanwise location at which an extrapolated g function from the tip reaches a value of one, (r/R)g=1, shows some dependence on tipspeed ratio.
These first observations give rise to the possibility that a general, though simple, g function may be found as an extension to the Glauert tip correction by means of two parameters, i.e. g Tip and (r/R) g=1 . The case of a tip-speed ratio of 5.42 (NREL Phase VI Rotor, J -Sequence) is taken as an example of approximating the WindDVE computed g function with a linear function in the tip region (solid black line in Fig. 2a). The approximated g function was implemented into XTurb through Eqn. (1). Computed results for the sectional normal force coefficient, C n , are shown in Fig. 2b in comparison to WindDVE and NREL data. It can be seen that a standard analysis performed by XTurb using the classical Glauert tip correction and a solution-based stall delay model (SBSD, [19]) results in the wellknown overprediction of blade tip loads. On the other hand, results obtained by WindDVE show  excellent agreement with NREL data in the tip region, and so do results obtained by XTurb (Schmitz&Maniaci, SBSD) when using the approximated g function shown in Fig. 2a. It is interesting to note that a simple approximated g function results in a nearly perfect match between a BEMT and a higher-fidelity method capable of computing vortex rollup. Of further note is that inboard of (r/R) g=1 , both XTurb solutions are the same as g = 1. Comparison against WindDVE and NREL data shows the benefits of the SBSD model to account for rotational augmentation, which is not included in WindDVE that essentially represents a XTurb solution that would be obtained without using the SBSD model. Furthermore, it is worth mentioning that the addition of a g function to XTurb results in the exact same number of iterations at every blade element; hence, the g function does practically not add any additional computational expense.

Results and Discussion
The proof-of-concept study presented in the previous section calls for extended parametric studies on WindDVE computed g functions that include variations in tip shape, blade tip pitch angle, and tip airfoil. The following subsections present parameter variations for the NREL Phase VI Rotor and the NREL 5-MW Turbine and the resulting g Tip and (r/R) g=1 extracted from the WindDVE results. Additionally, selected XTurb analyses and comparison with available data, WindDVE, and previous CFD work are included [20,21].  Fig. 2a. The WindDVE results were computed with a minimum wake length of 7 rotor diameters, resulting in a wake length convergence error for rotor coefficient of power of less than 1.5% for TSR less than 8 and less than 4% for TSR greater than 8. The NREL Phase VI rotor was computed at a constant rotor speed of 72 rpm with tip pitch angles of 0, 2, 3, and 6 degrees. As for the NREL 5-MW turbine, a complete power curve was computed, i.e. nearly constant tip-speed ratio (TSR) and 0 degree blade tip pitch in Region II (constant pitch, variable speed), and constant rotor speed and varying tip pitch angle in Region III (variable pitch, constant speed). Note that both the rotor planforms and the airfoils differ in the tip region of both rotors: The NREL Phase VI rotor is equipped exclusively with the S809 airfoil (α 0 = -1.23deg) and 7.3% c/R at r/R = 0.975; the NREL 5-MW turbine has the NACA 64_A17 as the tip airfoil (α 0 = -3.84deg) and only 2.4% c/R at r/R = 0.975. For future work, it is anticipated that the effects of tip airfoil and tip shape (or solidity) be investigated separately.  As observed in section 3.3, it seems that, at least to a good approximation, the g function parameters g Tip and (r/R) g=1 are primarily either functions of the tip-speed ratio (TSR) or the blade tip pitch angle (PITCH), with the latter having an additional uncertainty associated with differences in the zero-lift angle, α 0 . It is very interesting to note in Fig. 3 that data align linearly, which gives rise that a universal g function may be found, though subject to lower/upper constraints. The shift in (r/R) g=1 between rotors in Fig. 3a at constant slope and the difference in slope in Fig. 3b are conjectured to be a function of the airfoil zero-lift angle, α 0 , and the tip shape (or a representative local solidity parameter). These are the subject of ongoing and future studies. The blue circles in Fig. 3 are drawn around data points where quantitative comparisons are performed in the following subsections of XTurb analyses against measured data or available computational fluid dynamics (CFD) results. The objective of these comparisons is to investigate whether or not a simple extension to the Glauert tip correction according to Eqn. (1) consistently results in appreciable improvements in the prediction of blade tip loads. Figure 4 shows spanwise distributions of XTurb computed sectional normal and tangential forces, F n and F t , compared to measured NREL data. The parameters for g Tip and (r/R) g=1 were taken from Fig. 3 and implemented into XTurb through Eqn. (1). It can be seen in Fig. 4 that results compare very well with measured NREL data and that an extension to the classical Glauert tip correction by means of a simple g function (Schmitz&Maniaci) leads to improved prediction of sectional normal forces, F n , in the blade tip region. As far as the tangential forces, F t , in Fig. 3b are concerned, however, the tip loads are not as well predicted. One has to keep in mind, though, that i) tangential forces are an order of magnitude smaller than the normal forces and ii) tangential forces are subject to integration error from pressure tap measurements that are typically distributed using a quadrature method optimized for normal (or lift) force integration and also do not directly include surface friction. On absolute scales, the error in tangential forces, F t , remains very small; however, could be the subject of future work.

Selected Results -NREL Phase VI Rotor (Comparison against measured data)
A second example is given in Fig. 5 for a different tip-speed ratio (TSR) and blade tip pitch angle (PITCH). Again, values for g Tip and (r/R) g=1 are taken from Fig. 3. Shown in Fig. 5 are normal and tangential force coefficients, C n and C t , rather than corresponding forces. Note that the normal force coefficient, C n , is defined as the sectional force perpendicular to the local chord line, F n , and   normalized by the local dynamic pressure; the tangential force coefficient, C t , is defined as the sectional force parallel to the local chord line and facing to the leading edge, F t , normalized by the local dynamic pressure. This was simply done to demonstrate that both forces and force coefficients scale accordingly, which is important to show due to the respective change in local induction factors and associated scaling with (V rel /V 0 ) 2 . Results shown in Fig. 5 are consistent with the observations in Fig. 4; hence, this is supporting data that a simple, though general, extension to the Glauert tip correction may be found from the theory and analyses presented in this work.   Figure 6 shows comparisons of XTurb results against CFD results using the OVERFLOW code [20,21]. Note again that XTurb was run using the classical Glauert tip correction and using values for g Tip and (r/R) g=1 from Fig. 3. Shown are one representative case each of Region II and Region III operation. The NREL 5-MW turbine is a highly-loaded rotor in Region II, as can be seen in Fig. 6a. It is interesting to note that all observations made thus far are also true for a utility-scale turbine at representative wind speeds of the power curve. It is clear from Fig. 6 that the classical Glauert tip correction overpredicts tip loads as observed and known for many other rotors. On the other hand, XTurb results accounting implicitly for vortex rollup by means of a simple g function lead to a notable improvement in the prediction of blade tip loads. Of further note is that differences in rotor power amount to 1.8% in Region II at V 0 = 8m/s and 3% in Region III at V 0 = 15m/s. These are small but appreciable changes that, depending on the wind distribution, can be expected to be about 2.5% in turbine AEP. This is potentially important as the classical Glauert tip correction overpredicts turbine AEP, thus leading to higher-than-expected gains.

Conclusions and Future Work
This paper presented a novel idea of finding a simple extension to the classical Glauert tip correction that leads to improved BEMT predictions of blade tip loads on horizontal-axis wind turbines. At the core of the idea is the hypothesis that tip vortex rollup is a primary contributor to the seemingly universal overprediction of blade tip loads via BEMT methods that use the classical Glauert tip correction. It is oftentimes forgotten that the Glauert tip correction was designed for lightly-loaded rotors, an assumption that is practically invalid for most of the operating regime of modern horizontalaxis wind turbines.
A mathematical theory was presented that enables the extraction of an exact tip loss factor and corrective g function from a higher-fidelity method, provided that this method is capable of computing/extracting the spanwise distribution of the blade flow angle, Φ. In theory, this can be done by a CFD code; however, uncertainties in airfoil data and Φ are unavoidable. Instead, free-wake methods are a natural choice for this task due to their capability of computing tip vortex rollup and the fact that the exact same airfoil table data can be provided to the free-wake method as for the corresponding BEMT analyses. Some parameter variations were performed using the higher-order free-wake method, WindDVE, with the intent to find general relations for the parameters g Tip and (r/R) g=1 towards a simple extension to the Glauert tip correction. Though the present work is not complete, it provides good evidence that such a simple extension in the form of a g function can be found with future work and parameter sweeps that include different tip airfoils and tip shapes. Quantitative comparisons against measured data and available CFD analyses of tested g functions implemented into the BEMT code XTurb are very promising and, probably most importantly, consistent with one another in the improved predictive capability of blade tip loads. This provides good support and motivation for future research on BEMT tip modeling