Accurate loads and velocities on low solidity wind turbines using an improved Blade Element Momentum model

The accurate prediction of loadings and velocities on a wind turbine blades is essential for the design and optimization of wind turbines rotors. However, the classical BEM still suffer from an inaccurate prediction of induced velocities and loadings, even if the classical correction like stall delay effect and tip loss correction are used. For low solidity rotors, the loadings are generally over-predicted in the tip region, since the far wake expansion is not accurately accounted for in the one-dimensional 5 (1D) momentum theory. The 1D dimensional momentum theory supposes that the far wake axial induction is equal to twice the axial induction in the rotor plane, which results in an under-estimation of the axial induction factor in the tip region. Considering the complex nature of the flow around a rotating blade, the accurate estimation of 3D effects is still challenging, since most stall delay models still often tend to under-predict or over-predict the loadings near the root region. As for the solution method for the classical BEM equation, the induced velocities are computed accounting for the drag force. However, according to the 10 Kutta-Joukowski theorem, the induced velocities on a blade element are only created by lift force. Accounting for drag force when solving the BEM will result in an over-estimation of the axial induction factor, while the tangential induction factor is under-estimated. To improve the accuracy of the BEM method, in this paper, the 1D momentum theory is corrected using a new far wake expansion model to take into account the radial flow effect. The blade element theory is corrected for threedimensional effects through an improved stall delay model. An improved solution method for the BEM equations respecting 15 the Kutta-Joukowski theorem is proposed. The improved BEM model is used to estimate the aerodynamic loads and velocities on the National Renewable Energy Laboratory Phase VI rotor blades. The results of this study show that the proposed BEM model gives an accurate prediction of the loads and velocities compared to the classical BEM model.


Introduction
Thanks to its simplicity and robustness, the blade element momentum (BEM) theory has been widely used in the wind industry 20 for the design and optimization of wind turbine rotors (Tangler, 2002). It was first introduced by Glauert (1935), as a combination of the 1D momentum theory (MT) and the blade element theory (BET). The 1D MT, developed by Glauert (1935), is describing the process of energy extraction by the wind turbine using the conservation laws of fluid mechanics. The BET, developed by Froude (1878), is describing the local loadings and velocities computation on a blade using airfoil theory. Some assumptions in the BEM theory are corrected to make the prediction more realistic. The assumption of an infinite number of However, in order to close the equation system, They used the semi-analytical AD equation of radial flow taken from the AD simulation by Madsen et al. (2010). Nevertheless, these BEM models are strongly coupled and require multiple iterative loops to be solved. As a result, they are computationally expensive compared to the classical BEM model, especially when used for the optimization of wind turbines rotors.
Considering the complexity of three-dimensional flow around a rotating blade, the non-dimensional 2D lift and drag, extracted from a wind tunnel for a given airfoil, fail to predict the aerodynamic loadings on a rotating blade at the inboard sections of the blades. In the last decades, extensive studies of this phenomenon were performed and several correction models were proposed (Spera (2009), Butterfield et al (1992), Lee and Wu (2013), Snel et al. (1993), and Breton et al. (2008)). These studies agree on the fact that for low wind speeds and angles of attack bellow the static stall angle, the flow remains attached to the blade and there is a minor difference between the lift and drag of a rotating blade and a non-rotating blade. however, when 70 the 3D stall occurs at hight angles of attack, the flow began to separate from the blade surface and rotate with the blade. As a result, flow is subjected to rotational forces namely centrifugal and Coriolis forces. The centrifugal force induces an outboard radial flow, thus, the radial flow induces Coriolis forces that act as a favorable pressure gradient in the chord-wise direction.
These forces lead to a delay in the separation process, resulting in an increase in the maximum lift coefficient at a higher angle of attack. This phenomenon is known as stall delay or rotational augmentation. However, based on the comparative study 75 performed by Breton et al. (2008) on six existing stall delay models, they concluded that these models fail to predict the 3D behavior compared to NREL's phase VI experiment mainly because of their lack of generality since this phenomenon is not yet fully understood. However, compared to other stall delay models, Bak et al. (2006) model was found to have superior accuracy thanks to the modeling approach that was based on estimating the pressure gradient instead of estimating the gradient of the force due to rotational effects. However, Bak et al. (2006) model is still over-estimating the loads. This original approach was 80 adopted by Wang et al. (2013) to propose a new stall delay model by solving simplified Navier-Stocks equations for the pressure. Nevertheless, these models that are based on estimating pressure gradient are strongly empirical and showing a relatively good agreement with experimental data. However, Bak et al. (2006) model is not yet extensively validated, particularly, the force distributions along the blades (Breton et al. , 2008).
Moreover, The classical BEM model computes the induced velocities by solving iteratively the BEM equations accounting 85 for the drag force. However, the induced velocities are produced by vortices (lift or circulation) according to vortex theory and particularly the Kutta-Joukowski theorem (De Vries , 1979). Thus, the induced velocities are only created by lift force. Lindenburg (2003) has shown that in the case of a non-swept rotor blade, the relative induced velocity on the airfoil is a result of the circulation around it. Thus, it is perpendicular to the local relative flow. Wilson and Lissaman (1974) pointed out that the drag should be excluded from BEM equations because the drag-based velocity deficit is only a characteristic of the wake 90 and it does not contribute to the induced velocity at the rotor disc. Recently, Shen et al. (2005) made a detailed analysis of BEM equations in the tip region. It was found that if we account for the drag force in the BEM equation, the relative velocity becomes zero at the tip, which is not physically true since the tip vortex is created at the blade tip and then converted into the wake. Given these arguments, it is still arguably if the drag should be excluded from the BEM equations when computing the induced velocities ( Burton et al. (2001) and Hansen (2015) ). In either case, the drag force should be accounted for once the induced velocities are computed; since the drag force has an impact on the performance of a wind turbine rotor.
This manuscript is organized as follows. In Section 2, the GMT is applied to low solidity wind turbine rotors, then a new far wake expansion model is proposed using dimensional analysis. In Section 3, the blade element model is described, then an improved stall delay model is proposed. In Section 4, An improved BEM model respecting the Kutta-Joukowski theorem is given, the improved BEM equations are solved using the guaranteed convergence algorithm. In Section 5, A detailed validation 100 of the proposed BEM model is performed based on the experimental results of the NREL Phase VI rotor (Hand et al , 2001).
At last, the conclusion and further improvements of the BEM model will be given.
2 Improved momentum theory 2.1 General momentum theory applied to low solidity wind turbines The general momentum theory (GMT) is based on a global description of the flow around the wind turbine rotor. The rotor is 105 modeled as an actuator disc with an infinite number of blades. The flow is considered steady, incompressible and axisymmetric.
The GMT take into account the radial, axial and azimuthal velocities. The axial and radial velocity are continuous, while the azimuthal velocity has a discontinuity at the actuator disc. A rotation of the flow is present in the wake, and it is absent upstream of the rotor disc. The pressure drop at the rotor disc is due to the change in azimuthal velocity at the rotor disc (Sørensen , 2016). Assuming that the air masse contained in the rotor disc is separate from the rest of the air in the rotor plane, as a result, 110 the air masse at the rotor disc can be contained in a circular boundary surface. It can be extended upstream and downstream of the rotor forming a long stream-tube of circular cross-sections. Since the air within the circular boundary of the stream-tube is assumed to be incompressible; the cross-sectional area of the stream-tube must expand to accommodate the drop in wind speed. Thus, the conservation laws of fluid mechanics can be applied globally along the stream tube or locally along an annular stream tube (see Fig. 1).

115
The annular element of the actuator disc is taken at a radial position r and the actuator disc rotating with a rotational speed noted Ω. The annular element at the far wake is taken at a radial position r w and rotating with a rotational speed noted ω w . Due to the discontinuity in the rotational motion of the flow at the rotor plane, The tangential velocity at the rotor disk is taken as the average of the upstream and downstream plane V t = −rΩa . Since there is no rotation of the flow in the upstream plan, the tangential velocity is zero upstream of the rotor V t+ = 0 , the tangential velocity straight after the rotor V t− = −2rΩa , and 120 V tw = −r w ω w in the far-wake. In contrast, the axial and radial velocities are continuous, thus the radial and axial velocities straight after and before the rotor are identical. The radial velocity in the far upstream and downstream plans is zero since there is no expansion of the far upstream and downstream wake. In the rotor plane, the radial velocity is noted U rad . The axial (1) Where a, a , and b are: the axial induction at the rotor plan, tangential induction at the rotor plane, and axial induction at the far wake, respectively.
Applying the conservation laws under the assumption of the general momentum theory, the governing equations can be 130 found. The thrust (Eq. (3)) and torque (Eq. (4)) at a given radius of the rotor disc are given by (Sørensen , 2016). The GMT accounts for the pressure drop due to wake rotation, the wake expansion, and the pressure contribution to the lateral force component. Where, dT p,Side is the pressure contribution to the lateral force component, dA is the elementary disk area, and ρ is the air density. p w and p 0 are the pressure in the far wake and the atmospheric pressure respectively. λ r = rΩ/U 0 is the speed ratio.
The thrust coefficient in the GMT does not form a closed and solvable system because the unknowns number is too high compared to the number of the equations. Consequently, additional assumptions need to be introduced. The lateral force component, dT p,Side is difficult to determine and can be found using CFD, however, this term is generally considered to be 140 small (Sørensen , 2016). Thus, the annular element are independent of each other. Additionally, modern horizontal axis wind turbines (HAWT) are operating at high rotational speed in order to minimize the wake rotation losses and thus extract more power (Gupta and Leishman , 2005). As a result of high rotational speed, the rotor will have less solidity than a rotor operating at low rotational speed (Burton et al., 2011). Additionally, low solidity wind turbines ( typically 7% or less) will result in a low manufacturing cost of the blades (Tangler (2000) and Burton et al. (2011)). As a result, the pressure drop due to wake 145 rotation can also be neglected (p 0 = p w ). The thrust coefficient becomes similar to 1D MT except for the expansion effect that is parameterized by the far wake expansion ratio (χ = b/a) to be determined instead of taking b = 2a as the 1D MT. The hypothesis of an infinite number of blades is corrected using the Prandtl tip and hub loss model (Prandtl, 1921) . The thrust and torque are given by Eq. (5) and Eq. (6), respectively.
Prandtl tip and hub loss model F = F tip F hub is given by Eq. (7) and Eq. (8), respectively.
Where φ is the inflow angle, R is the tip radius, r h is the hub radius, and N is the number of blades.

155
Note that the radial flow is absent in the thrust coefficient equation because of the absence of radial flow in the far upstream and downstream to the actuator disc and it is only present in the near wake, reaching a maximum at the rotor disc when following a streamline (Van Kuik , 2017a). This was proved by Van Kuik (2017a) using vorticity-based arguments. As a result, When applying the conservation of momentum, the radial flow is canceled out, but its effect remains present in the far wake by the far wake expansion ratio.

New far Wake expansion model
The far wake expansion radius can be found by applying the conservation of the axial flow between the rotor disc and the far wake. It is given by Eq. (9).
In case of the 1D MT, the far wake axial induction is taken as twice the axial velocity in the rotor plan (b=2a), the 1D far 165 wake expansion radius is given by Eq. (10).  (2013), Sørensen and Kock (1995), and Zahle and Sørensen (2007)).Thus for a given axial induction factor, the far wake expansion radius for low solidity rotors will be further over-estimated. Consequently, the 1D MT will be valid only for a very low axial induction factor which corresponding to a low tip speed ratio. This is generally because the far wake expansion ratio in 1D MT is independent of operating conditions and geometry of the wind turbine rotor. To take into account the operating conditions and rotor's geometry effects on the far wake expansion ratio. First, the fractional decrease in axial wind speed between the rotor plane and the far wake is defined by Eq. (11). It is a normalized and dimensionless quantity between 0 and 1. In case of no expansion of the wake ξ = 0, and the limiting case ξ = 1 corresponds to the maximum wake expansion ratio ( b=2a) since b/a < 2 according to GMT. Second, the classical dimensional analysis (Graebel , 2001) is used to identify all relevant dimensional grouping to the far wake expansion, then applying the Buckingham 180 pi theorem to estimate the fractional decrease in axial wind speed ξ. Hence, the far wake expansion ratio can be computed by All dimensional parameters that are relevant to the far wake expansion are identified. They are shown in Table 1. The actuator disc geometry parameters include the blade tip radius R, the local radial station r, and the blade number N . The flow conditions 185 are the free stream speed U 0 , the rotor rotational speed Ω, and the additional velocity induced by the wake ν (m/s). The classical Buckingham pi theorem states that the number of dimensionless groupings is equal to the number of dimensional parameters subtracted by the number of independent dimensions (length, time). Note that the blade number is already dimensionless. thus, there are a total of three dimensionless groupings: the normalized blade radius, the tip speed ratio, and additional velocity induced by the wake normalized by the free stream. The blade number is the fourth dimensionless number (Table 2).
190 Table 1. Parameters relevant to far wake expansion ratio.
Blade geometry Flow conditions r R N U0 Ω ν Table 2. Dimensionless groupings relevant to far wake expansion ratio.
Blade geometry Flow conditions The above dimensionless grouping can be further reduced to one dimensionless group that characterizes the wake, which is the apparent helical pitch (Eq.12) relating the rotor parameters (N, λ, r/R) to the far wake parameters (h,ν).
Where, h is the helical pitch of the wake, andh is the normalized apparent helical pitch of the wake.
(r, R, N, λ, C T ) by an expansion model (Okulov et al. (2015)). The thrust coefficient is excluded in this work since it is already depending on the wake expansion and also to avoid a strong coupling of BEM equations. Thus, the dimensionless expansion radius r w /r depends on the apparent helical pitchh, dimensionless velocity induced by the wakeν. However, according to Okulov et al. (2015), for small expansion of the wake, the additional velocity induced by the wake of the tip vortices equals half the averaged induced axial velocity with a correction of a small expansion ( = r w − r) as shown in Eq.13. As a result, the 200 velocity induced by the wake has a small variation (0 ≤ν ≤ 0.25(1 + )), when the rotor is operating in the momentum region (0 ≤ a ≤ 0.5), compared to the apparent helical pitch varying from zero to infinity. Thus, the additional velocity induced by the wake can be neglected and the apparent helical pitch is considered as the most important dimensional parameter.
Before applying the Buckingham pi theorem to estimate the fractional decrease in axial wind speed, the normalized apparent 205 helical pitch should be normalized by the relative distance √ r 2 + h 2 instead of the blade radius only. According to the pi theorem, the fractional decrease in axial wind speed is given by Eq. (14). As a result, the far wake expansion ratio is given by Eq. (15).
Where a 0 and b 0 are empirical parameters. In this paper, these parameters are taken as a 0 = 1 and b 0 = 1 for the NREL Phase VI rotor. Once the far wake expansion ratio is found, the new far wake expansion radius can be computed using Eq.16 as follows: 3 Blade element model accounting for stall delay effect 215

Blade element theory
The blade element theory assumes that the blades are made up of a number of blade elements arranged in the spanwise direction, so the airfoil theory can be used to compute the local loads and velocities acting on each blade element. The BET assumes that the blade elements are aerodynamically independent and do not have any interference between them. The loads can be obtained from the 2D lift, drag and moment coefficients of the airfoil at any radial position. The inflow is known at the blade 220 and a lifting-line assumption is used (Branlard , 2017). As a result of these assumptions, the global loads on a blade can be integrated over the total blade span incorporating the velocity terms, to obtain the thrust, the torque, and power developed by the blade. This is further multiplied by the number of blades to get the total rotor thrust, torque, and power. A blade element is an infinitesimal fraction dr of the blade radius R at a radial position r, so it can be considered as an airfoil. An airfoil is defined with by its chord c(r), its twist θ(r) and its shape. The relative wind applied to a blade element, 225 noted U rel (r), is decomposed into a normal component U n (r) and a tangential component U t (r) to the rotor plane. Three characteristic angles defined by the velocity axis and the chord axis: the flow angle φ(r) between the tangential and relative velocity and it is assumed to be known, the airfoil twist θ(r) about to the rotor plane, and the angle of attack α(r) between the relative velocity and the chord axis and it is related to aerodynamic loads. Since the lifting-line assumption is used, the velocity triangle can be defined in terms of the axial and tangential inductions factors a and a (Eq. 17 and Eq.18). The relationship 230 between these angles and the velocity components are given in Equations 19 and 20. The blade element, velocity triangle, and aerodynamic forces are shown in Figure 3.
The lift and drag forces per unit of length applied on the blade element of the surface are defined as follows: 240 Where (C l2D (α), C d2D (α)) are the 2D lift and drag coefficients respectively. The normal and tangential forces per unit of length are defined as the projection of the lift and drag on the longitudinal plane and rotor plane (Fig. 3). They can also be defined in the function of normal and tangential coefficients.

245
Since the wind turbine rotor has identical blades, the loading on the blade elements of the blades at the same radius will be also identical. The local thrust, torque and power coefficients applied to an elementary annulus at a radial position r of the rotor disc are defined as follows: where σ = N c/(2πr) is the local solidity.

3D stall delay modeling
The stall delay phenomenon is still challenging current 1D models as it was shown by Breton et al. (2008) by performing a 255 comparative study of several stall delay model. However, Bak et al. (2006) model was found to have superior accuracy thanks to the modeling approach that was based on estimating the pressure gradient due to rotation instead of estimating the loading gradient. Although, Bak model is still over-estimating the loads. In addition, it is not yet fylly validated for the spanwise distribution of the loadings. In this work, The stall delay model originally developed by Bak et al. (2006) will be improved based on recent advances in this subject.

Bak model
Bak et al. (2006) was the first to develop a stall delay model that is based on estimating the pressure gradient instead of the force gradient due to the blade rotation. The pressure gradient is defined as the product of two functions: An amplification and a shape. The amplification was estimated as the ratio of the sum of centrifugal and Coriolis force to the pressure force, these forces were estimated through an order of magnitude analysis of NS equations. The shape was estimated through a simple 265 empirical equation based on fitting the 3D experimental data of NREL phase VI at 30% of the blade radius. The pressure gradient is given by Eq. 28.
Where α 0 is the AoA at which the flow begins to separate and α 1 is the AoA at which the flow is completely separated. For the S809 airfoil used in the NREL phase VI rotor. α 1 and α 2 can be approximately taken as 6.2 o and 21 o . x/c is the normalized 270 chord-wise position.
The normal and tangential forces gradient is found by integrating the pressure gradient over the full chord.

Improved Bak model
The stall delay model proposed by Bak assume that the chordwise velocity is zero and the spanwise velocity is dominant when the flow is separated. This hypothesis was also adopted by Corten (2001) based on an experimental study on the flow 280 separation on the wind turbine blade. As a result, this model is only valid inside the separated area where the radial flow is constant. Assuming that the pressure gradient induced by the spanwise flow in the separation area mainly affects the normal force of the local airfoil, since the pressure gradient depends only on the chord-wise position. This is similar to the analysis of (Lindenburg , 2003). The normal force gradient can be computed by integrating the pressure gradient in the separated area is given by Eq. 33.
Where f = 1 − x/c is the trailing edge separation factor. It can be estimated using the Kirchhoff-Helmholtz model (Eq. 34) that is commonly used in the dynamic stall computation (Leishman and Beddoes , 1986).
Note that the improved bak model is in good agreement with Wang et al. (2013) findings. Since, Wang et al. (2013) solved 290 the inviscid stall delay model developed by Corten (2001). They found that the shape of the normal force gradient to be a third-order plynomial as a function of the trailing edge separation factor, which is consistent with Eq. 33.

Performance tip loss correction
The blade element theory does not take into account the finite span of the blade; it treats all sections of the blade the same way.
In reality, this is not true because there is always aerodynamic loss at the tip of the blade (Shen et al. , 2005). The performance 295 tip loss is used to correct the 3D aerodynamic force coefficients; since the forces should converge to zero at the tip to equalize the pressure between the upper and lower surface of the blade. Shen's tip loss model (Shen et al. , 2005) is used in this paper.
It is given by Eq. (35).
The correction function h 0 is depending only on the number of blades and the tip speed ratio. The constants c 1 and c 2 are determined from experimental data. In case of the NREL phase VI rotor: c 1 = 0.125 and c 2 = 21.
The 2D aerodynamic data used in this work are taken from Lindenburg (2003) for the non-rotating blade. The performance tip loss model will be used in both Bak and improved Bak models. The 3D correction of the lift and drag are given by Eq. (37) and Eq. (38).

Improved BEM model
In order to compute the BEM solution (a, a , φ), the improved BEM model needs to be completed by respecting the KJ theorem and using a wake state model for axial induction factor a > a c . Then, the improved BEM equations will be solved using the 310 guaranteed convergence algorithm to avoid any convergence issues.

Induced velocities computational method
The conceptual model of Kutta-Joukowski for the wind turbine rotor is identical to the BET, except for the exclusion of drag force. Since, the induced velocities are produced only by lift force according to vortex theory and particularly the Kutta-Joukowski theorem (Okulov et al. , 2015). As a result, The orthogonality condition must be respected between the total lift 315 force and the relative velocity (see Fig. 4). By applying The Kutta-Joukowski theorem to a blade element of length dr, the total force is given by Eq. (39).
Where, Γ is the circulation vector defined by Eq. (40).
The momentum theory is valid up to a critical value of the axial induction factor (a < a c ). For higher induction the wind turbine operating in the turbulent wake state. Thus, an empirical relationship between the thrust coefficient and the axial induction factor is needed. In this paper, Buhl and Marshall (2005) approach is used because it does not suffer from numerical discontinuity between the theoretical and empirical thrust coefficients and also it shows a good agreement with experimental data. The proposed correction model is developed following the approach of Buhl and Marshall (2005). The thrust coefficient 330 for the axial induction factor larger than a critical value ( a > a c = 0.4) is given by Eq. (43).
It is clear that this equation is very similar to the original Buhl and Marshall (2005) correction. Equation 43 can be reduced to that of Buhl if relation the far wake expansion ratio χ = 2.

Algorithm solution 335
Historically, the solution method for solving the BEM equations is based on taking the axial and tangential induction factors as the unknowns, then solve the fixed-point problem (a, a ) = f (a, a ) using iterative method. This solution method is simple and fast but suffer from some convergence problems and does not always converge to the right solution (Maniaci , 2011) . Recently Ning (Ning , 2014) proposed a solution method that guaranteed the convergence which is based on reducing the two governing equations of the BEM model to one residual function Eq. (44), then solve it by a root-finding algorithm. The inflow angle is 340 the unknown and not the induction factors used as the unknown in the classical iterative algorithms.
The momentum/empirical region corresponds to local inflow angles in the range φ ∈ C(]0, π[). In this case, there is twoequation for the axial induction factor, a criterion is needed to distinct the momentum region (a < a c ) and the empirical region (a > a c ). The criterion is given by (η 0 = 2/3). The axial and tangential induction factors can be computed by equating the 345 thrust and torque in the MT and Kutta-Joukowski theory. The axial and tangential induction factors are given by Eq. (45) and Eq.(46), respectively.

355
Once the induced velocities are converged, the performances of the wind turbine rotor can be computed using BET, since the drag force is accounted for when computing the performances of the wind turbine rotor. This is illustrated by the following algorithm.
classical BEM. However, for briefness and relevance, the validation of the proposed improvements to the BEM method will be done by considering six different variations of the BEM method (Table 3)

3D BEM solution
The spanwise distribution of the converged angles of attack for different wind speed is shown in Fig. 5. The BEM-KJ-E,0 and BEM-KJ-E,1 models are used and compared with the estimated AoA from experimental data given by Sant et al. (2006).  In order to further highlight the effect of the far wake expansion model on the induced velocities. The relative velocity 390 and inflow angles, at r/R = 0.8 and U 0 = 7, predicted using the BEM-KJ-E,0 and BEM-KJ-E,1 models along with the ones estimated from experimental data are given in Figure 7. The axial induction factor predicted using the BEM-KJ-E,0 model is under-predicted due to the use of the 1D far wake expansion model. However, the new far wake expansion model improves the prediction of the axial induction factor, and it is in good agreement with the estimate axial induction from experimental data as shown in Fig.7. The tangential induction factor is unchanged for both models since the torque coefficients in both models 395 are identical. The inflow angle predicted using the BEM-KJ-E,1 model is also improved compared to the BEM-KJ-E,0 model, since the axial induction is under-estimated and the inflow angle is over-estimated. Thus, the new far wake expansion effect is canceled out when computing the relative velocity as shown in Fig.6. This behavior remains valid for other wind speed cases.
In order to further highlight the effect of excluding drag force on the induced velocities. The relative velocity and inflow angles, at r/R = 0.3 and U 0 = 20, predicted using the BEM-KJ-E,0 and BEM-KJ-E,1 models along with the ones from ex- perimental data are given in Figure 8. Using the velocity triangle, it is clear that the axial induced velocity is over-estimated in the BEM-KJ-E,0 model, while the tangential induction factor is under-estimated because the drag force is considered when solving the BEM equations, which contradict the KJ theorem. However, When the KJ condition is respected, the axial and tangential induction factors are improved and correlate well with the ones estimated from experimental data. Additionally, the inflow angle is slightly improved as shown in Fig. 5, and the relative velocity is significantly improved as shown in Fig. 6. This 405 behavior remains valid for other wind speed cases. that for low solidity wind turbines, the wake expansion radius converges after small distance downstream the rotor plan to an almost cylindrical wake. In the case of the NREL phase VI, the wake radius converges approximatively after one apparent helical pitch (Fig. 9). This was also confirmed by Madsen et al. (2010) showing that the axial velocity converges after just one diameter downstream of the rotor plane. Others like Sørensen and Kock (1995) and Zahle and Sørensen (2007) have also shown similar behavior.
415 Figure 9. Wake expansion radius at U0 = 7 5.2 3D stall delay model Figure 10 shows the 3D sectional loading at five radial locations (r/R = 30%, 47%, 63%, %80, and 95%). The BEM-S,0 and BEM-S,1 models are used to validated the improved stall delay model. It can be seen that both stall delay models are identical and predict relatively well the 3D behavior of the rotating blade, except in the separation area (α 0 < α < α 1 ). In the separated area, the improved Bak model predicts accurately the lift and drag coefficients, compared to Bak model that over-predict the 420 lift and drag due to integrating the pressure gradient over the whole chord. In fact, both stall delay models are only valid in the separated area where the radial flow is dominant compared to the chordwise flow. As a result, the pressure gradient is strongly dependent on the AoA, especially in the separation area, which was recently confirmed by Mauro et al. (2018) using CFD simulation. At r/R = 95%), the predicted loadings are in good agreement with experimental data, thanks to Shen tip loss model that used to correct the lift and frag coefficients in both stall delay models.

425
In both stall delay models, the amplitude of the pressure gradient depends on external forces at a given spanwise position, and the shape of the pressure gradient depends on the AoA and it is strongly empirical. Since it is taken from the 2D airfoil, then corrected using an empirical quadratic function. The amplitude is generally in good agreement with experimental data (Fig.12 and Fig.13), however, some over-estimation and under-estimation of sectional loadings are due to the shape of the pressure gradient that is strongly empirical and needs to be further studied. In fact, the shape of the pressure gradient can vary from 430 triangular shape to trapezoidal shape depending on the separation type: leading-edge separation or trailing edge separation. Schreck and Robinson (2002) has analyzed the pressure data on the NREL Unsteady Aerodynamics Experiment and found 22 https://doi.org/10.5194/wes-2020-43 Preprint. Discussion started: 12 February 2020 c Author(s) 2020. CC BY 4.0 License.
out that the 3D pressure coefficients have a trapezoidal as well as triangular shapes. The pressure shape variation was attributed to various mechanisms mediated by Coriolis and centrifugal forces and also to the presence of an energetic vortex structure. is due to the use of 1D MT expansion model, that results in an over-estimation of the AoA, resulting in an over-estimation of normal force coefficient. At low tip speed (U 0 ≥ 10m/s), the stall delay of Bak is responsible for the over-estimation of the spanwise distribution of both normal and tangential force coefficients at the hub region and at U 0 = 7, 10m/s, which is corresponding to the separation area. For the fully separated flow, both stall delay models are identical and they are in good