Evaluation of tip loss corrections to AD/NS simulations of wind turbine aerodynamic performance

: The Actuator Disc / Navier-Stokes (AD / NS) method has played a signiﬁcant role in wind farm simulations. It is based on the assumption that the ﬂow is azimuthally uniform in the rotor plane, and thus, requires a tip loss correction to take into account the e ﬀ ect of a ﬁnite number of blades. All existing tip loss corrections were originally proposed for the Blade-Element Momentum Theory (BEMT), and their implementations have to be changed when transplanted into the AD / NS method. The special focus of the present study is to investigate the performance of tip loss corrections combined in the AD / NS method. The study is conducted by using an axisymmetric AD / NS solver to simulate the ﬂow past the experimental NREL Phase VI wind turbine and the virtual NREL 5MW wind turbine. Three di ﬀ erent implementations of the widely used Glauert tip loss function F are discussed and evaluated. In addition, a newly developed tip loss correction is applied and compared with the above implementations. For both the small and large rotors under investigation, the three di ﬀ erent implementations show a certain degree of di ﬀ erence to each other, although the relative di ﬀ erence in blade loads is generally no more than 4%. Their performance is roughly consistent with the standard Glauert correction employed in the BEMT, but they all tend to make the blade tip loads over-predicted. As an alternative method, the new tip loss correction shows superior performance in various ﬂow conditions. A further investigation into the ﬂow around and behind the rotors indicates that tip loss correction has a signiﬁcant inﬂuence on the velocity development in the wake.


Introduction
Wind energy is nowadays an important and increasing source of electric power.It has been the biggest contributor of renewable electricity except for hydropower, sharing about 5.5% of the global electricity production in 2018 [1].As a primary subject of wind turbine technology, aerodynamics [2] largely determines the efficiency of wind energy extraction of an individual wind turbine or a wind farm.Along with the extensive development of high quality wind resources onshore, there are several trends in wind power industry: A large number of newly installed wind turbines have to be installed in areas with lower wind speeds and complex terrain; the layout optimization of wind turbine array becomes very important for wind farm design; offshore wind power development is accelerating, and the rotor size is continuously increasing.These trends need to be supported by more advanced and refined aerodynamic tools.More accurate aerodynamic load prediction is required for the design of a new generation of wind turbines with high efficiency and relatively low weight.Furthermore, a wind farm with dozens of wind turbines needs to be studied as a whole in order to take into account the complex interference between wind turbines.
Blade Element Momentum Theory (BEMT) [3,4] and Computational Fluid Dynamics (CFD) [5,6] are two essential methods of wind turbine aerodynamic computation.BEMT is no doubt the key method for rotor design [7], while CFD gives the most refined data of aerodynamic loads and flow parameters.A full CFD simulation with resolved rotor geometry is usually employed for an individual wind turbine [8,9].However, full CFD is not suitable for a wind turbine array, due to the huge computational cost.The Actuator Disc/Navier-Stokes (AD/NS) method [10][11][12][13] was developed for this situation, in which the rotor geometry is not resolved, and thus, the number of mesh cells is greatly reduced.AD/NS is based on a combination of the blade-element theory and CFD.The flow field is still solved by CFD, while the rotor entity is replaced by a virtual actuator disc on which an external body force is acted.If the blade entities are represented by virtual actuator lines, it is called the Actuator Lines/Navier-Stokes (AL/NS) method [14 -17].As compared to AD/NS, the flow around the rotor solved by AL/NS is closer to the reality, since the rotor vortices are simulated.However, AL/NS requires more grid cells to describe the actuator lines, and thus, is seldom employed in wind farm simulations.
The AD/NS method has played a significant role in wind farm simulations involving wake interaction [18,19], complex terrain [20,21], atmospheric boundary layer [22,23], noise propagation [24,25], etc.Nevertheless, the AD model assumes that the number of blades is infinite, and thus, no tip loss is simulated, causing an over-prediction of the blade tip loads and power output.In order to take tip loss into account, a reliable engineering model has to be embedded into the numerical solver, which is usually called tip loss correction.However, all tip loss corrections were originally proposed for BEMT, and there is no tip loss correction specially developed for AD/NS.Most of the literature about AD/NS simulations either did not mention tip loss or declared that the Glauert tip loss correction [26] was employed.It is worth mentioning that AD/NS and BEMT solve the momentum of the flow past the rotor by using two completely different approaches.The former employs CFD, while the latter applies the momentum theory, though they share the actuator disc assumption and the blade-element theory.That leads to different implementations of tip loss correction for the two methods.In the BEMT, tip loss correction is realized by applying a correction factor F into the induced velocity through the rotor.In the AD/NS, the velocity is naturally found by the NS solver, and factor F can only be used to modify the external body force.The question is whether a tip loss correction has a good global performance when it is transplanted into AD/NS.Even for the BEMT itself, evaluations indicate that accurate tip loss correction is still not well achieved [27], and the development of new correction models is still going on [28][29][30][31][32][33].In contrast with the massive study in the BEMT, tip loss corrections applied to AD/NS lack a comprehensive evaluation.
In the present work, the 2-Dimensional (2D) axisymmetric AD model is employed, and steady-state simulations are performed.The code solves the incompressible axisymmetric NS equations for the experimental NREL Phase VI wind turbine [34] and the virtual NREL 5MW [35] wind turbine under various axial inflow conditions.The main purpose of the numerical study is to evaluate the performance of the tip loss corrections applied to AD/NS.Three different implementations of the Glauert tip loss correction are discussed.In addition, a new tip loss correction recently proposed by Zhong et al. [33] is introduced and compared.The normal and tangential forces of blade cross-sections are chosen as the key parameters for the present evaluation.Because of the long arm of force, computational errors of the forces acted on the tip region are most likely to cause non-ignorable errors of the blade bending moment and the rotor torque (power generation).A BEMT study on the NREL Phase VI wind turbine by Branlard [27] showed that the power generation would be overestimated by 15% if no tip loss correction was made and a deviation of about 5% exists when various existing tip loss corrections were applied.That highlights the significance of the present study on tip loss correction.
The innovation of the present study lies in the following items.(a) Different implementations of the Glauert tip loss factor F are gathered, discussed and compared, and finally, one of them is recommended according to its best performance.Such kind of work has never been reported in the existing literature.(b) The Glauert correction is found to have a similar performance when it is transplanted from BEMT to AD/NS.(c) The new tip loss correction of Zhong et al. [33] is for the first time employed in AD/NS simulations.It is found to be generally superior to the Glauert-type corrections.That provides an alternative choice for a more accurate prediction of the blade tip loads.(d) Tip loss correction is found to have a significant influence on the velocity field, which highlights the importance of an accurate tip loss correction for not only the blade loads, but also the wake development.
The paper is organized as follows: In Section 2, the axisymmetric AD/NS method is introduced, including the governing equations of the NS approach and the formulae of the AD model; In Section 3, the Glauert and the new tip loss corrections are introduced; In Section 4, the implementations of the tip loss corrections used in AD/NS are described; In Section 5, the numerical setup, as well as the involved simulation cases, are introduced; In Section 6, all simulation results are presented and discussed; Final conclusions are made in the last section.

Governing Equations
The governing equations of the present AD/NS simulation are the incompressible axisymmetric NS equations.In a cylindrical coordinate system as shown in Figure 1, the axial direction is defined along the z-axis, the radial direction is represented by r, and the tangential direction is represented by θ, the continuity equation is written as the axial and radial momentum equations read and the tangential momentum equation is

Force on Actuator Disc
The conceptual idea of the AD/NS method is to solve the aerodynamic force of rotor blades by using the blade-element theory and then applying its counterforce as an external body force into the momentum equations.In the above equations, u z /u r /u θ is the axial/radial/tangential velocity, p is the static pressure, ρ is the air density, ν is the kinematic viscosity coefficient of air, f z , f r , f θ are the source term of the axial, radial, tangential external body forces, respectively.

Force on Actuator Disc
The conceptual idea of the AD/NS method is to solve the aerodynamic force of rotor blades by using the blade-element theory and then applying its counterforce as an external body force into the momentum equations.
At each time step of solving the NS equations, the velocity passing through the actuator disc is detected and used to determine the axial induced velocity W a and the tangential induced velocity W t , The axial interference factor a and tangential interference factor a are then determined by, where V 0 is the wind speed, and Ω is the rotating speed of the wind turbine rotor.
The above interference factors are azimuthally unchanged according to the axisymmetric condition, which implies an infinite number of blades and zero tip loss.In order to estimate the tip loss of the realistic rotor with a finite number of blades, the interference factors need to be corrected as ã = f corr (a ), (10) where ã and ã denote the corrected axial and tangential interference factors, f corr and f corr represent correction functions.The axial and tangential induced velocities after the correction are written as According to the velocity triangle shown in Figure 2, the flow angle φ, angle of attack α, and relative velocity V rel are then determined by where β is the local pitch angle of a blade cross-section.According to the relationship of force projection, the normal and tangential force coefficients are determined by cos sin The axial force z F and tangential force F θ per radial length of the rotor are then given by ( ) ( ) where c is the local chord length, and B is the number of blades.Ignoring the effect of the radial flow, the vector form of the counterforce acting on the air is ( ) Rather than distributing the force only on the disc, the above force is regularized by the following Gaussian distribution along the axial direction in order to avoid a numerical singularity.The Gaussian distribution has been widely used and proven to be proper for AD/NS and AL/NS simulations [10,[14][15][16][17].
where 0 z is the axial position of the disc.The parameter ε serves to adjust the concentration of the regularized force, which is in the present study set to ε = 0.02R where R is the rotor radius.The source terms of the external body force in the momentum Equations (2-4) need to be replaced with the above force per unit volume.In the axisymmetric coordinate system, a 2D grid With the determined angle of attack and relative velocity, the lift coefficient C l and drag coefficient C d of each cross-section of the blade can be obtained from tabulated airfoil data.According to the relationship of force projection, the normal and tangential force coefficients are determined by The axial force F z and tangential force F θ per radial length of the rotor are then given by where c is the local chord length, and B is the number of blades.Ignoring the effect of the radial flow, the vector form of the counterforce acting on the air is Rather than distributing the force only on the disc, the above force is regularized by the following Gaussian distribution along the axial direction in order to avoid a numerical singularity.The Gaussian distribution has been widely used and proven to be proper for AD/NS and AL/NS simulations [10,[14][15][16][17].
where z 0 is the axial position of the disc.The parameter ε serves to adjust the concentration of the regularized force, which is in the present study set to ε = 0.02R where R is the rotor radius.
The source terms of the external body force in the momentum Equations ( 2)-( 4) need to be replaced with the above force per unit volume.In the axisymmetric coordinate system, a 2D grid cell with an axial length of ∆z and a radial height of ∆r represents a volume of 2πr∆r∆z, and thus, the force per unit volume is Using Equations ( 20) and ( 21), it reads

Glauert Correction
The tip loss correction of Glauert [26] was developed from the study of Prandtl [36].A function F, which was later recognized as the first tip loss factor, was derived by Prandtl making Betz's optimal circulation [37] go to zero at the blade tip, where B is the number of blades, r is the local radial location, R is the rotor radius, and λ is the tip speed ratio.The original derivation of Prandtl was written very briefly and was elaborated more in [4,27,38] as for a more detailed derivation.Later on, Glauert made further contributions.First, he interpreted the physical meaning of factor F as the ratio between the azimuthally averaged induced velocity and the induced velocity at the blade position, leading to where a = 1 2π 2π 0 adθ and a = 1 2π 2π 0 a dθ are the azimuthally averaged axial and tangential interference factors, and a B and a B are the interference factors at the blade position.Second, the local inflow angle φ was introduced into the function in order to make it consistent with the local treatment of the BEMT, leading to the following new formula of factor F, The factor F was then applied into the momentum theory through a straightforward way of multiplying the axial velocity at the rotor plane with factor F, resulting in the following thrust and torque for an annular element, dT = 4πrρV 2 0 a(1 − a)Fdr, Using another two equations of dT and dM derived from the blade-element theory, and the velocity triangle at the rotor plane, the equations for the interference factors in the BEMT approach was derived to be: The above two equations lead to the final iterative formulae of the BEMT approach with the Glauert tip loss correction.Their difference from the baseline BEMT approach is only the appearance of factor F in the equations.Obviously, the application of the Glauert correction in BEMT is very simple, which is also a great advantage.There are several variations of the Glauert correction in which the way of applying factor F is different [39,40], but the original Glauert tip loss correction is the most commonly used form till today.

A Newly Developed Correction
A new tip loss correction was recently proposed for BEMT by Zhong et al. [33], based on a novel insight into tip loss.In contrast with the Prandtl/Glauert series corrections that begin with an actuator disc and estimate the effect of the finite number of blades, the new correction begins with a non-rotating blade and estimates the effect of rotation on tip loss.It has been validated in BEMT computations and showed superior performances in the cases involving flow separation or high axial interference factor.
The correction was realized by using two factors of F R and F S that treat the rotational effect and the 3D effect, respectively.
where c denotes a geometric mean chord length, and S t is the projected area of the blade between the present cross-section and the tip.The purpose of introducing c is to deal with tapered blades and those with sharp tips.The factor F R was applied to the BEMT approach by using: in which the employed lift and drag coefficients were corrected by factor F S , rather than the direct use of the airfoil data. where where m is the slope of the lift-curve of the airfoil before flow separation, α i is called the downwash angle, and α e is called the effective angle of attack.By comparing with the Glauert tip loss correction, the new tip loss correction appears more complicated in application.It is simplified in the present study: First, Equations ( 36) and (37) will not be employed naturally in the AD/NS method; Second, Equations ( 38)-( 41) are further simplified to The simplification is derived by using the fact that α i and C d are relatively small.More detailed applications are shown in the next section.32) and (33) where the Glauert tip loss factor F is applied are not employed in the AD/NS method because the flow is simulated by the NS solver instead of the momentum theory.As a result, an alternative way has to be employed correctly to apply the tip loss factor F. Nevertheless, among the literature studies, little literature mentions the detail of how the tip loss factor is applied to the AD/NS simulations.After an extensive literature review, we have found three representative documents in which Sørensen et al. [41], Mikkelsen [42], and Shen et al. [43] described their implementations.In order to facilitate the distinction, the implementations adopted by the three studies are denoted as Glauert-A, Glauert-B and Glauert-C in the present paper, respectively.

Glauert-A Correction
Sørensen et al. [41] performed a tip loss correction by applying factor F to modify the aerodynamic force that determines the external body force in the NS equations.They replaced the lift coefficient C l by C l /F (there was no need to modify the drag in their study as the drag was assumed not to produce the external body force).In the present Glauert-A correction, C n and C t , in Equations ( 18) and ( 19), are replaced by: That is equivalent to replacing C l by C l /F and C d by C d /F, according to Equations ( 16) and ( 17).Sørensen et al. [41] did not explain the reason why factor F could be directly used to modify the force.They might be inspired by Equations ( 32) and (33) in which the existence of F can be looked as corrections to C n and C t , although from a physical point of view it is a correction to the interference factors.It is important to note that the corrected force is only used for determining the external body force, so as the flowfield, while the force acting on the blade should be regarded as the original one.In addition, the interference factors in Equations ( 9) and (10) should no longer be corrected, leading to ã = a and ã = a .This implementation involves a division operation with denominator F, which causes unreasonable big values of Cn and Ct at the extreme tip where F → 0 .However, there is no exact criterion for defining what a big value is unreasonable because the Cn and Ct themselves are introduced as an engineering correction rather than a physical concept.Sørensen et al. [41] did not mention this problem, and no obvious numerical fluctuation is observed in his result.That is possible because the problem is limited to a very small area at the tip, and thus, the integral effect of the resulted unreasonable body force is also small.

Glauert-B Correction
Mikkelsen [42] compared Equations (32) and (33) to the corresponding baseline equations without factor F. The baseline equations are Appl.Sci.2019, 9, 4919 9 of 21 In order to distinguish from the parameters in the baseline equations, we here rewrite Equations ( 32) and (33) to The implementation of Mikkelsen [42] adopts the following equations which implies an assumption of φ = φ which in fact does not accurately hold, The corrected interference factors were, thus, determined by The tip loss correction was completed as the above functions for ã and ã were used to replace Equations ( 9) and (10).

Glauert-C Correction
Shen et al. [43] performed their correction by directly using Equation (26) which represents the physical meaning of factor F. Using the azimuthally uniform condition (a = a), the correction was completed by replacing Equations ( 9) and (10) Similar to the Glauert-A correction, this implementation involves a division operation with denominator F. That causes an unphysical big value of ã at the extreme tip where F → 0 .A limiter is set in the present study to force the result to be ã = 1 when it is greater than 1.The value of ã is usually much less than ã and is not limited here.

Application of New Correction
The application of the new correction consists of two steps: The first is using factor F R determined by Equation ( 34) and the second is using factor F S determined by Equation (35).
In the first step, an implementation similar to the Glauert-C correction is adopted, which is performed by replacing Equations ( 9) and (10) A limiter of ã = 1 is used when a is greater than 1.The angle of attack α can then be determined by calculations using Equations ( 11)-(14).The second step is to correct the lift and drag coefficients by using Equations ( 42) and (43), the result of which is used to replace the C l and C d in Equations ( 18) and (19).

Flow Solver and Mesh Configuration
The 2D axisymmetric NS equations are solved by using an in-house code EllipSys2D [44,45] developed at Technical University of Denmark (DTU).The code is a general incompressible flow solver with multi-block and multi-grid strategy.The equations are discretized with a second-order finite volume method.In the spatial discretization, a central difference scheme is applied to the diffusive terms and the QUICK (Quadratic Upstream Interpolation for Convective Kinematics) upwind scheme is applied to the convective terms.The SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations) scheme is used for the velocity-pressure decoupling.The turbulence flow is simulated using the method of Reynolds-averaged Navier-Stokes equations (RANS) in which the k-ω turbulence model of Menter [46] is employed with a modification for a better simulation of the turbulence quantities in the free-stream flow [47,48].
The coordinate for the AD/NS simulations is defined in the z-y plane.A computational mesh is generated, as shown in Figure 3 where the inflow, outflow and axisymmetric boundaries are indicated, the length of the computational domain is 30, which is nondimensionalized with the rotor radius.The blade is positioned at z = 0 and the grid cells are clustered around z = 0 and y = 1 to ensure a better resolution near the blade tip.The mesh is composed of four blocks where 64 × 64 grid points are used for each block.There are 64 grid points on the AD along the radial direction, and 20 points in the range of [-0.02, 0.02] in the axial direction where the Gaussian distribution of the body force plays a significant role.Since the axisymmetric boundary condition is applied, the 2D flow solution is regarded as an azimuthal slice of a full 3D field.

Simulation Cases
Simulations are performed for two different wind turbines of the NREL Phase Ⅵ [34] and the NREL 5MW [35] that represent a small and a large rotor size, respectively.Additionally, the NREL Phase Ⅵ rotor has a very blunt blade tip shape as compared with the NREL 5MW rotor.
The operational conditions of the two rotors are listed in Table 1.Cases 1, 2 and 3 are set to be consistent with the NREL UAE experiment in axial inflow conditions [34].Various wind speeds of 7 m/s, 10 m/s and 13 m/s are considered, while the rotating speed remains unchanged.The flow is fully attached on the blade surface at a wind speed of 7 m/s, but is partly separated at 10 m/s and 13 m/s (Higher wind speed leads to heavier flow separation) [8].Measured pressure distributions (from which the force can be obtained by pressure integral) at five blade sections (r/R = 0.30, 0.47, 0.63, 0.80, 0.95) are available for these cases [50].Cases 4 and 5 are two typical points on the designed power curve of the NREL 5MW wind turbine [35].The wind speed and rotating speed of Case 5 are the rated parameters of this wind turbine.Case 4 represents a condition with a lower wind speed of 8 m/s at which the rotating speed is reduced to 9.22 rpm for tracking the optimum tip speed ratio.
These cases cover the following multiple situations: Wind turbines with remarkably different rotor sizes and tip shapes; flow with fully attached and separated conditions; operating conditions

Simulation Cases
Simulations are performed for two different wind turbines of the NREL Phase VI [34] and the NREL 5MW [35] that represent a small and a large rotor size, respectively.Additionally, the NREL Phase VI rotor has a very blunt blade tip shape as compared with the NREL 5MW rotor.
The operational conditions of the two rotors are listed in Table 1.Cases 1, 2 and 3 are set to be consistent with the NREL UAE experiment in axial inflow conditions [34].Various wind speeds of 7 m/s, 10 m/s and 13 m/s are considered, while the rotating speed remains unchanged.The flow is fully attached on the blade surface at a wind speed of 7 m/s, but is partly separated at 10 m/s and 13 m/s (Higher wind speed leads to heavier flow separation) [8].Measured pressure distributions (from which the force can be obtained by pressure integral) at five blade sections (r/R = 0.30, 0.47, 0.63, 0.80, 0.95) are available for these cases [50].Cases 4 and 5 are two typical points on the designed power curve of the NREL 5MW wind turbine [35].The wind speed and rotating speed of Case 5 are the rated parameters of this wind turbine.Case 4 represents a condition with a lower wind speed of 8 m/s at which the rotating speed is reduced to 9.22 rpm for tracking the optimum tip speed ratio.These cases cover the following multiple situations: Wind turbines with remarkably different rotor sizes and tip shapes; flow with fully attached and separated conditions; operating conditions with various wind speeds and rotating speeds.It is, therefore, interesting to see the performance of the tip loss corrections in these cases.
Experimental data are preferred as the reference for the computational results of the NREL Phase VI rotor, as well as the full CFD data are also displayed in some results.For the cases of the NREL 5MW rotor, there is no experimental data available such that the full CFD data are the only reference.The related full CFD simulations were previously conducted by the authors, see Zhong et al. [8,33].

Glauert-A/B/C Corrections
The blade loads for Cases 1, 2 and 3 are gathered in Figure 4, which represent results obtained from the NREL Phase VI rotor.The NREL Phase VI blade has a linear change of chord distribution starting from r = 1.3 m with a constant slope of 0.1.The chord length is 0.358 m at the tip, which is about 48% of the largest chord length in the blade inboard part.Therefore, the loads do not converge to zero without a tip loss correction.The blade loads for Cases 1, 2 and 3 are gathered in Figure 4, which represent results obtained from the NREL Phase Ⅵ rotor.The NREL Phase Ⅵ blade has a linear change of chord distribution starting from r = 1.3 m with a constant slope of 0.1.The chord length is 0.358 m at the tip, which is about 48% of the largest chord length in the blade inboard part.Therefore, the loads do not converge to zero without a tip loss correction.It is noticed that there is a certain degree of difference between the results of the Glauert-A/B/C corrections in the tip region (e.g., at r/R > 0.8).In general, Glauert-C leads to lower loads which are closer to the experimental data, Glauert-A gives results close to, but slightly higher than, those of Glauert-C, while Glauert-B results in relatively higher loads near the tip.Taking the n F in Figure 4a as an example of quantitative comparison, the relative difference between the results of Glauert-B and Glauert-C is about 4% at r/R = 0.95.The difference in other cases is not larger than this value.
At a wind speed of 7 m/s, all the curves of n F generally agree with the five experimental data except that at r/R = 0.95 where an overestimation is observed.This overestimation becomes much more remarkable as the wind speed increases to 10 m/s and 13 m/s.At 10 m/s, there is no much difference around r/R = 0.9 between the results with no correction and with the Glauert-type corrections, indicating that almost no effective correction is made here.At 13 m/s, the curves of the Glauert-type corrections even exceed the uncorrected curve in a r/R range of about [0.66, 0.88].It is clear that the corrections perform much worse at the higher wind speeds as compared with 7 m/s.The overestimation at r/R = 0.95 is also observed in the results for t F , although it looks better to some extent, especially at 7 m/s.Considering the fact that the Glauert correction was derived based on the potential hypothesis, it is reasonable to believe that the unusually poor performance at the higher wind speed is caused by the flow separation.According to the results of the NREL UAE Phase Ⅵ experiments [50], the angle of attack at most cross-sections of the blade exceeds the linear range of their aerodynamic polar when the wind speed is increased to be higher than 10m/s.Taking the cross-section of r/R = 0.63 as an example, the angles of attack are 5.9°, 11.8°, 16.9° at a wind speed of 7 m/s, 10 m/s, 13 m/s, respectively.The exact operating points for three cases on the lift polar of the S809 airfoil (the airfoil of all cross-sections of the NREL Phase Ⅵ blade) are depicted in Figure 5   The normal and tangential forces along a blade of the NREL 5MW rotor are plotted in Figure 6.As a widely studied virtual rotor, full CFD solutions [33] are often used as the reference data.It is noticed that there is a certain degree of difference between the results of the Glauert-A/B/C corrections in the tip region (e.g., at r/R > 0.8).In general, Glauert-C leads to lower loads which are closer to the experimental data, Glauert-A gives results close to, but slightly higher than, those of Glauert-C, while Glauert-B results in relatively higher loads near the tip.Taking the F n in Figure 4a as an example of quantitative comparison, the relative difference between the results of Glauert-B and Glauert-C is about 4% at r/R = 0.95.The difference in other cases is not larger than this value.
At a wind speed of 7 m/s, all the curves of F n generally agree with the five experimental data except that at r/R = 0.95 where an overestimation is observed.This overestimation becomes much more remarkable as the wind speed increases to 10 m/s and 13 m/s.At 10 m/s, there is no much difference around r/R = 0.9 between the results with no correction and with the Glauert-type corrections, indicating that almost no effective correction is made here.At 13 m/s, the curves of the Glauert-type corrections even exceed the uncorrected curve in a r/R range of about [0.66, 0.88].It is clear that the corrections perform much worse at the higher wind speeds as compared with 7 m/s.The overestimation at r/R = 0.95 is also observed in the results for F t , although it looks better to some extent, especially at 7 m/s.
Considering the fact that the Glauert correction was derived based on the potential hypothesis, it is reasonable to believe that the unusually poor performance at the higher wind speed is caused by the flow separation.According to the results of the NREL UAE Phase VI experiments [50], the angle of attack at most cross-sections of the blade exceeds the linear range of their aerodynamic polar when the wind speed is increased to be higher than 10m/s.Taking the cross-section of r/R = 0.63 as an example, the angles of attack are 5.9 • , 11.8 • , 16.9 • at a wind speed of 7 m/s, 10 m/s, 13 m/s, respectively.The exact operating points for three cases on the lift polar of the S809 airfoil (the airfoil of all cross-sections of the NREL Phase VI blade) are depicted in Figure 5.It clearly shows that Case 2 and Case 3 are in the nonlinear lift region corresponding to flow separation, while Case 1 is in the linear lift region corresponding to attached flow, which implies that flow separation occurs at 10 m/s and 13 m/s.
(c) Case 3: V0 = 13 m/s, Ω = 72 rpm  The normal and tangential forces along a blade of the NREL 5MW rotor are plotted in Figure 6.As a widely studied virtual rotor, full CFD solutions [33] are often used as the reference data.(The Reynolds number of the airfoil is 1 × 10 6 which is close to that of the cross-section).
The normal and tangential forces along a blade of the NREL 5MW rotor are plotted in Figure 6.As a widely studied virtual rotor, full CFD solutions [33] are often used as the reference data.Although the typical rotating speed of the NREL 5MW rotor is a few times less than the NREL Phase VI rotor, the resulted tip speed ratio is higher, which is closer to the optimal design point.The blade tip of the NREL 5MW rotor is smoothly sharpened where the slope at the last 4 m is 0.46, and the chord length at the tip is only 4.3% of the largest chord length in the blade inboard part.It can be seen in Figure 6 that the loads at the tip naturally approach to zero even without a tip loss correction.Although the typical rotating speed of the NREL 5MW rotor is a few times less than the NREL Phase Ⅵ rotor, the resulted tip speed ratio is higher, which is closer to the optimal design point.The blade tip of the NREL 5MW rotor is smoothly sharpened where the slope at the last 4 m is 0.46, and the chord length at the tip is only 4.3% of the largest chord length in the blade inboard part.It can be seen in Figure 6 that the loads at the tip naturally approach to zero even without a tip loss correction.
As a large wind turbine with pitch control, there is no flow separation on the most area of the blade, including the tip.That means the flow pattern ideally keeps the same on the blade as the wind speed increases, which is significantly different from the situation of the NREL Phase Ⅵ rotor.
As seen in Figure 6, by increasing the wind speed from 8 m/s to 11.4 m/s, the forces are greatly increased, whereas, their distributions near the tip region remain similar.All the Glauert-A/B/C corrections over-predict the normal forces near the blade tip, whereas, better agreement with the reference data is found for the tangential forces.Such a result is to some extent similar to that of the NREL Phase Ⅵ rotor at 7 m/s.Slight discrepancies are also noticed between the three Glauert-type corrections, and the Glauert-C performs better among them.

Comparison with BEMT Results
The Glauert tip loss correction was originally proposed for BEMT.As applied to AD/NS simulations, the correction factor F works in a different way.It is worth to mention that the specific difference in the final results would be caused by the transplantation from BEMT to AD/NS.In  As a large wind turbine with pitch control, there is no flow separation on the most area of the blade, including the tip.That means the flow pattern ideally keeps the same on the blade as the wind speed increases, which is significantly different from the situation of the NREL Phase VI rotor.As seen in Figure 6, by increasing the wind speed from 8 m/s to 11.4 m/s, the forces are greatly increased, whereas, their distributions near the tip region remain similar.All the Glauert-A/B/C corrections over-predict the normal forces near the blade tip, whereas, better agreement with the reference data is found for the tangential forces.Such a result is to some extent similar to that of the NREL Phase VI rotor at 7 m/s.Slight discrepancies are also noticed between the three Glauert-type corrections, and the Glauert-C performs better among them.

Comparison with BEMT Results
The Glauert tip loss correction was originally proposed for BEMT.As applied to AD/NS simulations, the correction factor F works in a different way.It is worth to mention that the specific difference in the final results would be caused by the transplantation from BEMT to AD/NS.In order to perform a comparison, BEMT computations under the same conditions of the present study are conducted.
The comparisons were shown by percentages, which represent the relative change of blade loads before and after using a tip loss correction, where F n,No and F n,Gl are the normal forces obtained from computations with no tip loss correction and with the Glauert correction (the standard Glauert correction in BEMT and the Glauert-type corrections in AD/NS), respectively, F t,No and F t,Gl are the corresponding tangential forces.The above parameters are calculated independently using BEMT and AD/NS.The results are shown in Figure 7 where the Glauert-BEMT denotes the results for the standard Glauert correction used in BEMT, and the Glauert-A/B/C denotes the results for the Glauert-type corrections used in AD/NS.A certain degree of difference between the Glauert-BEMT and the Glauert-A/B/C results is observed.The Glauert-C again shows the best performance, since it is generally most consistent with the Glauert-BEMT, which means the transplantation of the tip loss correction from BEMT to AD/NS does not notably influence its performance.The results for Cases 3 and 4 also agree with the above conclusion, which is not displayed here for simplicity.

forces.
The above parameters are calculated independently using BEMT and AD/NS.The results are shown in Figure 7 where the Glauert-BEMT denotes the results for the standard Glauert correction used in BEMT, and the Glauert-A/B/C denotes the results for the Glauert-type corrections used in AD/NS.A certain degree of difference between the Glauert-BEMT and the Glauert-A/B/C results is observed.The Glauert-C again shows the best performance, since it is generally most consistent with the Glauert-BEMT, which means the transplantation of the tip loss correction from BEMT to AD/NS does not notably influence its performance.The results for Cases 3 and 4 also agree with the above conclusion, which is not displayed here for simplicity.

Results of New Tip Loss Correction
The simplified form (see Section 3.2) of the new correction developed by Zhong et al. [33] is used in the present simulations.The results of blade loads for Cases 1, 2 and 5 are displayed in Figure 7. Distributions of δ n and δ t a rotor blade.

Results of New Tip Loss Correction
The simplified form (see Section 3.2) of the new correction developed by Zhong et al. [33] is used in the present simulations.The results of blade loads for Cases 1, 2 and 5 are displayed in Figure 8a-c, respectively.Case 1 represents a typical condition of the NREL Phase VI wind turbine with the fully attached flow on its blades, while Case 2 represents another condition of flow separation (stall).Case 5 represents the rated operating condition of the NREL 5MW wind turbine.The results for Cases 3 and 4 are not displayed here for simplicity.(The performance of the correction in Case 3 is similar to that in Case 2, while the performance in Case 4 is similar to that in Case 5.) The Glauert-C correction, which performs best in Glauert-A/B/C, is compared with the new correction.Full CFD results [33] are employed as the reference data for the NREL 5MW wind turbine.As compensation for the sparse points of the experimental data, full CFD data [8,33] are also employed for the NREL Phase VI wind turbine.

Comparison of Velocity Field
The blade loads are affected by the local flow around the blade, and the forces consequently modify the flow around the blade and further in the wake.The flowfield simulated by the AD/NS It is seen in the result of Case 1, as shown in Figure 8a, that the full CFD result achieves the best agreement with the experimental data and the new correction is superior to the Glauert-C correction.As compared to the Glauert-C correction, the overestimation of F n in the tip region is properly corrected towards the reference data by the new correction.A similar tendency is seen from the tangential force distribution.In Case 2, as shown in Figure 8b, the new correction maintains its accuracy, while the Glauert-C correction produces larger positive error, due to the occurrence of flow separation.The performance of the two corrections is similar to those in BEMT computations [33].
For the NREL 5MW rotor, as shown in Figure 8c, the new correction again makes a better correction to the blade loads in the tip region.The corrected curve of F n for the new correction is closer to the reference data compared to that for the Glauert-C correction.The difference between the two corrections in tangential force (F t ) is small, and both the corrections can be considered to be doing well.

Comparison of Velocity Field
The blade loads are affected by the local flow around the blade, and the forces consequently modify the flow around the blade and further in the wake.The flowfield simulated by the AD/NS demonstrates how the tip loss corrections affect the velocity field.In order to clearly show the influence of using different tip loss corrections, instead of directly showing the velocity field, the following velocity difference is defined, where u z,new and u z,Gl are the axial velocity obtained from the simulations using the new tip loss correction and the Glauert-C tip loss correction, respectively.Figure 9 displays the contours of δu z in the flow field of the studied two rotors, which highlights the influence of the tip loss corrections in terms of changing the velocity field not only around the blade tip, but also in the wake.The influence of δu z is more concentrated near the blade tip in the rotor plane, whereas, the inner part is hardly affected.In the wake, the influence range of δu z tends to cover a wider radial space downstream.The magnitude is even increased in the wake of the NREL 5MW rotor.

)of 21 Figure 1 .
Figure 1.Definition of coordinates and the stream tube through an actuator disc.

Figure 1 .
Figure 1.Definition of coordinates and the stream tube through an actuator disc.

C
of each cross-section of the blade can be obtained from tabulated airfoil data.

Figure 2 .
Figure 2. Illustration of the velocity and the aerodynamic force for a blade cross-section.

Figure 2 .
Figure 2. Illustration of the velocity and the aerodynamic force for a blade cross-section.

Figure 3 .
Figure 3. Mesh configuration for the Actuator Disc/Navier-Stokes (AD/NS) simulations.The above computational setup has been proven to perform well in AD/NS simulations, as shown in the previous work of Cao et al. [49] where both blade loads and wake flows were validated against experiments.

Figure 4 .
Figure 4. Normal force n F and tangential force t F along the NREL Phase Ⅵ blade.

Figure 5 .
Figure 5. Operating points of the cross-section of r/R = 0.63 on the lift polar of the S809 airfoil.(The Reynolds number of the airfoil is 1 × 10 6 which is close to that of the cross-section).

Figure 4 .
Figure 4. Normal force F n and tangential force F t along the NREL Phase VI blade.

Figure 4 .
Figure 4. Normal force n F and tangential force t F along the NREL Phase Ⅵ blade.

Figure 5 .
Figure 5.Operating points of the cross-section of r/R = 0.63 on the lift polar of the S809 airfoil.(The Reynolds number of the airfoil is 1 × 10 6 which is close to that of the cross-section).

Figure 5 .
Figure 5. Operating points of the cross-section of r/R = 0.63 on the lift polar of the S809 airfoil.(TheReynolds number of the airfoil is 1 × 10 6 which is close to that of the cross-section).

Figure 6 .
Figure 6.Normal force Fn and tangential force Ft along the NREL 5MW blade.

Figure 6 .
Figure 6.Normal force F n and tangential force F t along the NREL 5MW blade.

Figure 7 .
Figure 7. Distributions of n δ and t δ along a rotor blade.

Figure
Figure 8a, b and c, respectively.Case 1 represents a typical condition of the NREL Phase Ⅵ wind turbine with the fully attached flow on its blades, while Case 2 represents another condition of flow separation (stall).Case 5 represents the rated operating condition of the NREL 5MW wind turbine.

Figure 8 .
Figure 8. Normal force Fn and tangential force Ft along the NREL 5MW blade.

Figure 8 .
Figure 8. Normal force F n and tangential force F t along the NREL 5MW blade.

Table 1 .
Simulation cases for two different wind turbines.