Analysis of the DANAERO wind turbine field database to assess the importance of different state‐of‐the‐art blade element momentum (BEM) correction models

Aerodynamic loads of wind turbine blades are often predicted by manufacturers using the blade element momentum (BEM) theory, for which many corrections have been proposed in the literature. The physical impacts of such corrections on field measurements have seldom been assessed because of the relative unavailability of dedicated measurements. Based on the unique full‐scale database of the DANAERO project, available through the IEA (International Energy Agency) Task 29, this work incrementally applies on aerodynamic field measurement improvements of the BEM theory: atmospheric boundary layer vertical velocity gradient, neighboring wake, yaw misalignment, wind inflow location, tower shadow effect, cone angles modeling, blade aeroelastic deformation, and dynamic wake. This is performed using the iBEM method (inverse Blade Element Momentum), which back‐calculates the aerodynamic coefficients (lift—CL and drag—CD) using aerodynamic loads from field tests.

inflow. 9 Recently, 3D Computational Fluid Dynamics (CFD) polars rather than 2D wind tunnel measurements have been used showing better agreement with high fidelity tools results. [10][11][12] The physical impacts of such corrections on field measurements have seldom been assessed because of the relative unavailability of dedicated in situ measurements.
Nowadays, the design is more and more helped by more advanced and accurate models: free vortex wake (FVW) and computational fluid dynamics (CFD). The benefit of these tools, compared with BEM, is to provide a greater accuracy in the physics analyzed but at a higher computational cost. 13 It is also possible to couple to the aforementioned tools structural solvers and/or controllers in order to simulate the wind turbines in a so-called aero-servo-elastic tool. Consequently, this coupling enhances even more the model accuracy along with the computational power needed. The use of the correction models, higher fidelity tools, and new methods allows the wind turbine design to be more cost-effective, increases its capacity factor, and decreases the levelized cost of energy (LCOE). 14 However, to achieve grid parity the LCOE needs to decrease even further. To help predicting realistic and more accurate local aerodynamic blade loads and develop low cost computation models, complex full-scale field experiments are necessary. However, such field measurements are scarce because of the significant cost involved in such complex experiments. To the authors knowledge, only in the DANAERO project a full-scale field test campaign was performed, spanning over three months, between July and September 2009 15 including many different operating conditions (wind direction and intensity, misalignment, wake effects from neighboring turbines). Alongside with these field tests, wind tunnel tests of blade sections (aerofoils) were performed. Also, 2D and 3D CFD simulations of the identical aerofoils and rotor were performed. This important database was available through the IEA group Task 29 Phase IV led by G. Schepers since 2010 and used in many studies. 7,[15][16][17][18][19][20][21][22][23][24][25][26] The focus of IEA Task29 Phase IV Work Package 3 was the validation, improvement, and understanding of aerodynamic models as implemented in wind turbine design codes by analyzing detailed aerodynamic measurements from the DANAERO experiment. Different phenomena were targeted: the impact of the turbulent inflow, the yawed conditions, the wake inflow, 2D/3D aerofoil characteristics, aeroelastic effects, the transition, and acoustics. Despite the significant learning brought by the IEA Task29 Phase IV, the accomplished work on the DANAERO database cannot identify the major phenomena responsible for the observed lift and drag dispersion.
The present work objective was to perform further analysis of this database to help in the identifications of the BEM hypotheses with the most significant impacts on aerodynamic loads calculations by comparing the wind tunnel measured aerodynamic coefficients with the ones calculated from field measurements. As a first approach, this work is restricted to the mean data binned by angle of attack. The origin of the instantaneous dispersion being another important challenge that is out of the scope of the present paper. Compared with the work of Troldborg et al, 27 the main difference is on the transfer function used between the 2D database in ideal inflow condition (low turbulence intensity, 2D axi-symmetric conditions) and field measurements. In the work of Troldborg et al, 27 the transfer function is estimated by using an optimization tool and a 1-minute time series. The present study is adding incrementally to the steady BEM formulation, using a 30 minutes long time series, the effect of atmospheric boundary layer vertical velocity gradient, neighboring wake, yaw misalignment, and wind inflow location. Then, the unsteady BEM formulation is assessed by adding to the model: tower shadow effect, yaw and cone angles modeling, blade aeroelastic deformation, and dynamic wake. Doing so, it is possible to highlight and quantify the model corrections that have the most impact on the local aerodynamic loads. This is performed using the iBEM (Inverse Blade Element Momentum), which back-calculates the aerodynamic coefficients (lift-C L and drag-C D ) using a set of aerodynamic loads, operating conditions, atmospheric conditions, and blade geometry. This paper is separated into four distinct sections: first, the presentation of the field data used; second, a presentation of the iBEM method and its verification; then, the presentation of the results after analyzing three time series of 10 minutes each; and finally, a conclusion.

| DANAERO DATABASE
The results presented in this paper are based on the large field database generated by the Danish project DANAERO 15,16 and available through the membership to the IEA (International Energy Agency) Wind Task 29. The IEA Task 29 Phase IV is the continuation of previous tasks (I to III) which focused on wind turbine measurements both in the field and in wind tunnels. Only details useful for the present study will be given in the present paper, more details on the DANAERO database are available. 28 The section will be separated in two parts: a first part dedicated to a brief description of the field test site with its instrumentation and a second part presenting the motivations of the chosen test cases.

| Field test site and instrumentation
Field tests of the DANAERO project took place in the Tjaereborg wind farm located in Western Denmark (see Figure 1). The turbine instrumented within the wind farm is a NM80 2MW (80 m rotor diameter), labeled later in the text "Reference turbine" (see Figure 2). One of the LM38.8 blades (38.8  radial stations starting from the blade root as summarized in Table 1. It allows in particular the normal and tangential forces computation at each radial position. The pressure outputs were acquired simultaneously with the reference turbine operating information, including the electrical power production, the blade azimuthal position, the rotor speed, and the blade pitch. A meteorological mast was also installed South West (main wind direction) of the reference turbine to record the rotor inflow. For that purpose, met mast sensors (see Table 2) acquired atmospheric data simultaneously with the reference turbine operating information and blade instruments. The blade tip height ranges between 17 m and 96m. The met mast anemometers were installed at several heights to capture the wind speed across the rotor. Finally, the wind farm altitude is at sea level and does not vary significantly throughout the wind farm.

| Test cases choice
In this paper, we will focus on the normal operating conditions, 29 that is, the turbine facing the incoming wind, with limited wake effects from neighboring turbines. The pitch and RPM settings are driven by the turbine controller according to the normal operating conditions. The chosen window of analysis is the measurement campaign done on 16 July 2009 between 12:30 hours and 13:00 hours as it reaches the normal operating conditions. Figure 3 shows the wind rose of the times series selected for the present study, measured by the met mast, with superimposition of the Tjaereborg wind farm turbines. Highlighted in blue, the reference turbine is located downstream of the met mast (in red) recording the rotor inflow. The data were binned in 10 minutes during the measurement campaign. The data in Table 3 show the mean environmental values (wind speed, wind direction, and turbulence intensity) for each of the three 10min time series files analyzed. The mean wind speed is rather low (6.5 m/s), and the turbulence intensity is around 9%. With a low wind speed, one can expect high induction levels which will lead the turbine to operate in the so-called "turbulent wake state" (a > 0.5). The range of wind direction over the 30 minutes analyzed is contained within 20°, its standard deviation is 5°.

| Atmospheric stationarity and stability tests
Prior to performing the statistical analysis on the atmospheric dataset, the stationary assumption, in the Foken sense, has been checked for all time series measured by the met mast. According to Foken, 30 the wind stationarity is assessed when, within a 30min period, the difference between the complete time series average of W ′ ′ and the average of the sum of 5 minutes-pieces for W ′ i ′ i is <30% (see Equation 1).
Tjaereborg wind farm, the instrumented wind turbine is shown along with the other wind turbines. The met mast is also shown southwest to the instrumented turbine. 15 The wind turbines are circled in blue and the met mast in green F I G U R E 2 Tjaereborg wind farm layout including the reference turbine, the blue square, studied in the present paper. The red triangle indicates the location of the meteorological mast where W′ is the wind speed fluctuation in the vertical direction and θ′ is the temperature variation across the time period considered. The atmospheric stability was computed using the Monin-Obukhov similarity theory, 31,32 and it is to be noted that none of the time series analyzed show neutral conditions. 33 We can therefore expect a high level of turbulence due to the thermal exchanges. This atmospheric condition is obviously far from the wind tunnel conditions for measuring 2D aerofoil polars used by the BEM method.

| IBEM METHOD
In order to derive the aerodynamic coefficients from the field measurements, an iBEM method was developed, validated, and verified. This is detailed in the present section.

| iBEM method description
The iBEM method is the application in reverse of the classical BEM theory. BEM is used to design wind turbines, calculate associated loads, and compute the turbine AEP (Annual Energy Production). iBEM allows conversely, to back-calculate from the aerodynamic loads the aerodynamic coefficients for a given inflow time series. In this article, the aerodynamic loads are obtained from two distinct sources: numerical results from CFD for verification purposes and from field measurements for analysis. By adapting the BEM procedure, 34 it is possible to calculate for a given turbine operating condition: the local velocity, the C L and C D as well as the induction and the angle of attack from integrated pressure (producing the normal F n and tangential forces F t ). The main steps are described below, and a graphical format with detailed equations is shown in Figure 4: 1. The axial induction factor a and the tangential induction factor a′ are first estimated (typically a = a′ = 0). 2. Then, the inflow angle is estimated from the instantaneous velocity inflow V w , the rotor rotational speed , and the local radius r. 3. The angle of attack, α, is computed using the Blade Element Theory (BET) with θ the local twist angle and β the blade pitch angle. 4. The relative velocity V rel is computed using the velocity triangle from the BET 34 (see Figure 5). 5. To account for the finite blade span, the Prandtl's tip correction factor is calculated. 6. The initial induction coefficients, a and a′, are updated accounting for highly loaded rotors (eg, using Glauert's correction or Spera equation). 7. The unsteady BEM equations can be applied: yaw models, dynamic wake model, blade acceleration due to its deflection and tower shadow effect. 8. A convergence criterion, , is defined and the iteration process restarts from step 2 until the convergence criterion is reached. It is defined here as a j − a j−1 < a and a � j − a � j − 1 < a � with a and a ′ equal to at least several order of magnitude lower than the expected induction value. Taking a = a � = 10 − 10 , the process takes on average less than 10 iterations to converge for each time step, where the subscript j denotes the iteration for a chosen time step. 9. After convergence, the time step (subscript i) is incremented and the process is repeated starting from step 1.
In the present iBEM approach, the aerodynamic coefficients are extracted from the forces measured thanks to the pressure taps along the instrumented blade (see Table 1). The integration of the pressure gives the normal force and the pressure part of the tangential force. The viscous part is not measured which, at low angle of attack, adds large uncertainty in the drag coefficient (cosine in front of the tangential force in Equation 16) and only small uncertainties in the lift force (sine in front of the tangential force in Equation 16).

| Alternative models for the induction
The iBEM method developed by the authors compares several existing formulations to model the induction (axial and tangential): the BEM method as documented by Hansen 34 (see dotted-line box in Figure 4) and the one from Madsen. 7 After comparing both methods, the authors decided to merge both implementations to create a third modeling. Thereafter, the results will be named as follow and detailed below:

| Madsen modeling
In Madsen et al, 7 the authors defined an alternative approach to the axial induction modeling for highly loaded rotors. Rather than using the conditional approach, 34 they have fitted a polynomial curve following the curve shape for low induction values and the highly loaded correction (see Figure 6). The blue circles represent the evolution of the induction following the BEM equation: and the tangential induction is derived with respect to the axial induction through:

| Hybrid modeling
The Hybrid model combines the axial induction model derived by Madsen et al 7 (see Section 3.2.1) (Equation 18) and the tangential induction described by Hansen (Equation 7). The tangential induction is therefore described independently from the axial induction. It is possible to merge those two modelings because both are build on independent empirical formulations.

| iBEM verification
To verify the iBEM solver, benchmark cases in the form of CFD outputs of a wind turbine immersed in a constant and uniform flow for several wind speeds. The CFD cases were performed in RANS (Reynolds Averaged Navier-Stokes) with a Shear Stress Tensor (SST) k − turbulence model. The IEA Tasks 29 members extracted aerodynamic results using the AAT (Average Azimuthal Technique). 36,37 The CFD data are obtained from G. Bangga (personal communication). The results are presented in Table 4. The blade is rigid, the inflow is uniform and constant (without shear or yaw), and the rotor operating conditions are constants (RPM, pitch). The numerical environment is therefore controlled and within the limits of the BEM theory. The local velocity (V rel ), lift force (F L ), and drag force (F D ) are calculated thanks to the other given inputs using, respectively, the Equation (4), and the inverted Equation (17) (where the forces are the unknown). The axial induction modeling differences from BEM solvers and CFD computations (visible in Figure 6 for instance) can then be evaluated similar to the work of Rahimi et al. 38

| Verification results
The iBEM verification outputs are first analyzed by comparing the outcome between the iBEM model output and the CFD computations at section r∕R = 49%, for V w = 6.1 m∕s . The forces (F L and F D ), the aerodynamic coefficients (C L , C D ), the relative velocity (V rel ), the axial induction (a), the tangential induction (a′), and the angle of attack (α) are analyzed and plotted. Figure 7 shows that the lift and drag force calculations are well predicted for all three induction modelings (ie, Hansen, Madsen, and Hybrid). The aerodynamic coefficients (C L , C D ), the relative velocity (V rel ) are within ± 0.1% of the reference, depending on the induction model. The Hansen model seems to slightly underpredict the relative velocity leading to a marginally higher lift coefficient. However, the axial induction (a), tangential induction (a′), and the angle of attack (α) results show large differences between the models. Madsen and Hybrid results show an axial induction error <1%, while Hansen model is closer to 10% error. The tangential induction results are similar (6% error across all models) despite the slightly better prediction from Hansen. Finally, the angle of attack prediction shows 10% error when using Hansen model and less than 1% when using either Madsen or Hybrid.
Along the blade span, results from the Madsen modeling and Hybrid modeling are overlapping each other (see Figure 8 left) and close to the reference case (within 1%). The worse results for the Hansen model highlighted in Figure 7 are confirmed for the entire blade span. The axial induction starts to deviate from the other models as early as r∕R = 35% . It is to be noted that all models fail to capture the axial induction behavior calculated by CFD between r∕R = 70% and r∕R = 90% (see Figure 8 left), similarly found in Rahimi et al. 38 The tangential induction modelings show differences with the CFD reference up to r∕R = 45%, after this radial position the agreement is better. The independence of the tangential modeling used in Hansen and Hybrid is captured and is beneficial toward the tip of the blade (see Figure 8 right). Madsen modeling for tangential induction is underpredicting compared with the CFD reference from r∕R = 82% until the tip of the blade (see Figure 8 right). The calculated aerodynamic coefficients, presented in Figure 9 show very good F I G U R E 4 Diagram of the inverse Blade Element Momentum (iBEM) and Unsteady inverse Blade Element Momentum (UiBEM) method algorithm (for definition of variables used in Tower Shadow, Yaw and Dynamic wake model, please refer to the general nomenclature at the beginning of the document). The process is repeated until convergence is reached for each time step of the acquired inflow velocity and for each blade section. The axial induction and tangential induction in the convergence loop are marked with subscript j. The black path shows the use of BEM equations, the blue path shows the use of unsteady BEM equations, and the red path shows outputs of the present solver, which comprise of the measured forces projected in the aerofoil coordinates, C L and C D , the angle of attack α and the converged axial and tangential induction F I G U R E 5 BEM velocity triangle showing the inflow angle , twist and pitch angle + and angle of attack . The normal force (F n ), tangential force (F t ), lift force (F L ), and drag force (F D ) are also represented showing the relationship between rotor and local coordinates agreement between the CFD and all BEM models, albeit the Hansen results are slightly underpredicted.
By combining the tangential induction model from Hansen and the more accurate axial induction from Madsen, the Hybrid model yields results marginally better than the other two models. The verification is repeated for all wind speeds listed in Table 4 (results not showed here), and it appears that the Hybrid model predicts results closer to the input for all wind speed. Following the presented results and the previous outcomes from Rahimi et al, 38 the iBEM method is considered validated and verified. For the rest of the study, the results presented will be using the Hybrid model.

| FIELD MEASUREMENTS ANALYSIS
As explained in Section 2.3, the normal and tangential forces are measured in different radial positions along the blade. When applying the verified iBEM method to the selected time series, it is possible to calculate the lift and drag, thus evaluating the phenomena modeled by the BEM theory for each blade section and each time step. First, the comparison with the raw data is presented (Section 4.1), and then, the different filtering and correction models developed with their impact on the average aerodynamic coefficients will be shown. Only the lift coefficient will be presented as the pressure sensors cannot measure the complete drag force. The turbine is operating at high induction level (a ≥ 1∕3), and the calculations may contain errors due to possible inaccurate high induction modeling.
In their work, Ozçakmak et al 26,39 performed the analysis of field transition measurements for the tip section, focusing on the boundary layer transition laminar to turbulent. They showed that the turbulent inflow and surface roughness were the main contributors in the boundary layer transition. Transition on the suction side occurs, in the field, at approximately 5% of the aerofoil chord from the leading edge with little to no fluctuation of the mean location. The measurement is very well correlated to 3D CFD results where the transition is also predicted to occur at this location (ie, 5% from the leading edge at the tip section). The mid-span section, r∕R = 49%, showed that the transition on the suction side was at approximately 10% from the leading edge. It is close to the wind tunnel tripped position and rather far away from the predicted 2D natural transition which varies between 50% and 30% (XFOIL, Re c = 5 × 10 6 ). The wind tunnel clean conditions are therefore nonrealistic operating conditions. It was expected as experiments performed in the wind tunnel are done at very low levels of the turbulent intensity (0.1%) while the turbulent intensity in field measurements is around 9% F I G U R E 6 Different induction modelings assuming the Prandtl's tip correction factor F = 1 and the critical induction factor a crit = 0.2

Inputs Symbol Value Units
Wind speed V w 3, 6.
Tangential force coefficient Angle of attack (see Table 3). Also, as already known, the inboard area of the wind turbine blade is sensitive to rotational effects where the lift characteristics tends to be enhanced by the Coriolis forces. 2,37,[40][41][42] Therefore, for a fairer comparison with field data, different corrected values of the 2D wind tunnel tests are included to the C L plots including: • wind tunnel tests with tripped conditions at the blade surface, available from IEA task 29. • rotational effect, also called "3D corrections," of the 2D wind tunnel tests using the model derived by Chaviaropoulos et al. 4

| Raw comparison
The expected dispersion between 2D wind tunnel measurements, 2D CFD simulation (EllipSys 43 data) and field instantaneous C L using a steady iBEM method without correction or filter, is shown in Figure 10. Figure 11 shows the averaged C L iBEM output for the test cases described in Section 2.2. The results are compared with wind tunnel measurements for blade sections from r∕R = 34% , closest to the hub, to r∕R = 95%, closest to the blade tip. Using the calculated local velocity, the local Reynolds number for the different sections is approximately 4.5 × 10 6 . This value is close enough from the wind tunnel measurement of 5 × 10 6 to allow meaningful comparison since different wind tunnel test types will capture the possible effect seen in the field. All sections show a decreasing C L trend for an increasing α. Hence, even when selecting appropriate wind direction/speed and normal rotor operation, it is still not suitable to compute the aerodynamic coefficients using steady iBEM directly from field measurements.

| Filters and corrections
The next subsections will show the different steps the authors followed in order to improve the iBEM results accuracy and assess the importance of each phenomena in the BEM formulation. It is performed by the evaluation of the wind tunnel measurements and the back computed C L using iBEM for the radial section r∕R = 49%. Only one section is presented, because this mid-span section is less likely to be affected by the flow three dimensionality: root and tip effects, thus ensuring the most appropriate comparison. The study will evaluate various physical phenomena that may violate BEM hypothesis: the location of the wind inflow measurement (time lag), atmospheric boundary layer vertical gradient, wakes of neighboring turbines and yaw misalignment. Then, the application of the unsteady BEM equations will be evaluated and the tangential forces measurement approximation will be assessed.
The different steps can be classified in two main categories: filtering applied to the time series and correction of inputs for the iBEM method. The filters and correction methods developed are cumulative, that is, the results presented in each subsection will also account for any previously correction introduced. For instance, the neighboring wake filter presented in 4.2.2 will also include the Mean vertical gradient filter presented in 4.2.1. The unsteady BEM results will only use the wake filter and the vertical gradient correction as the other filters are addressed in the equations themselves. Each analyzed case is numbered from 0 to 7.

| Vertical gradient filter
The mean vertical gradient from the atmospheric boundary layer flow is not taken into account by the BEM theory. By analyzing the vertical gradient of the mean wind speed, it has been decided to remove from the dataset any time step whose inflow gradient, measured by the mast, is above a threshold. The limit has been imposed based on the average value of the steepest part of the vertical mean gradient of several time series analyzed. Its value is 0.0025 s −1 , which corresponds to 50 m height, close to the hub height, a region where the speed gradient is the most pronounced. A filter is therefore applied to exclude the data when the blade is passing below the hub. Analyzing only the upper part of the rotor reduces the influence of the mean vertical gradient and its associated turbulence and uncertainties due to the tower shadow effects (ie when blades are passing in front of the wind turbine tower). It should be pointed out that, obviously, this operation removes approximately 50% of the available data points. Figure 12 shows the effect of removing these data points in the time series. The raw data are shown in blue for reference, and this selection method is shown to have a negligible impact on the C L ( ) curve, since the remaining data are overlapping with the raw time series. Figure 13 shows that there is an overlap between the reference wind turbine and the wake from its neighbors during the 30 minutes analyzed. To ensure that the met mast is not reading any wake data and the reference wind turbine is not in the wake of a neighboring one, a dedicated filter was developed. Figure 14 shows the wind rose changes when the wake is accounted for. The wind farm layout is over imposed with the wind direction and intensity. The length of each colored line represents the probability of occurrence during the time series. Also, to ensure that the wake is well-developed, the Frandsen 44 wake development model is used and represented by black dotted line circles. The wind farm turbines are outside the black dotted circle, meaning the wake is well developed when it reaches the next turbine in the wind farm.

| Neighboring wake filter
The Jensen model is therefore a valid assumption. Using the Jensen wake model, 45 it is possible to calculate the wind speed deficit contribution of each wind turbine at any position in the wind farm. The circular red curve shows the wake F I G U R E 7 Comparison between the CFD reference and the different axial induction modelings at 6.1 m/s for r∕R = 49% F I G U R E 8 a and a ′ -Full blade comparison between the reference data and the different axial induction modelings at V w = 6.1 m/s. The square represents the reference data, the circle represents the Hansen modeling, the triangle represents the Madsen modeling, and the diamond represents the Hybrid modeling | 1487 contribution of all the wind turbines seen by the met mast, and the circular blue curve shows the same wake contribution seen by reference wind turbine, but graphically centered on the mast to allow for a visual comparison of the wake and instrumented wind turbine wake situation. Having both curves tangent would mean that there is no wind speed deficit due to the wake of neighboring wind turbines. For instance, it can be seen in Figure 14 left, in the sector between 240° and 250° that both curves are not tangent. The reference turbine is operating in the wake of another turbine. After application of the method, the inflow data are updated and only the wind direction where both circular curves are tangent is plotted in Figure 14 right. After filtering out the selected data samples, the inflow data are clean of any wake effect. Figure 15 shows the effect of removing the data points from the time series on C L when either the reference turbine or the met mast is in the wake of a neighboring turbine. Only a slightly positive impact is observed: a small lift increase and small standard deviation reduction.

| Yaw misalignment filter
An important limitation of the BEM theory is the assumption of rotor alignment with the inflow. In reality, large deviations can occur within 10 minutes (see Figure 16). The blue curve represents the turbine yaw position and the orange curve represents the wind direction measured by the met mast. It can be seen that the turbine orientation does not change throughout the 10 minutes time series. The horizontal black dotted lines represent the yaw error allowed, in this example set to be ±5°, around the turbine position. The choice of ±5° allows to deviate slightly from the BEM theory while keeping enough data points for the analysis. The allowed yaw error removes 50% of the entire dataset.
As seen in Figure 17, removing from the time series the data points where the yaw error is above an arbitrary threshold (here ±5°) seems to have no impact on the averaged aerodynamic coefficients. On the other hand, allowing an even larger yaw error will violate the BEM theory assumption of inflow alignment. A yaw correction model in the BEM model would certainly be more appropriated to keep all points for statistics purposes.

| Vertical gradient correction
Rather than filtering out data within the vertical gradient from the atmospheric boundary layer flow (see Section 4.2.1) and using a single height measurement as inflow, the present section corrects the wind inflow using known laws describing the vertical gradient velocity profile from atmospheric boundary layer flows (eg, Monin-Obukhov, power law). It will associate to the analyzed section the calculated inflow based on the azimuthal position. While the Monin-Obukhov method is more accurate than the power law method to rebuild the wind F I G U R E 9 C L and C D -Full blade comparison between the reference data and the different axial induction modelings at 6.1 m/s. The square represents the reference data, the circle represents the Hansen modeling, the triangle represents the Madsen modeling, and the diamond represents the Hybrid modeling profile, it requires 10 minutes averaged values to build a velocity profile. In our case, it is not fast enough since we are building the wind profile for every time step measured at a sampling frequency of 35 Hz. Thus, the power law approach is preferred, the difference between both methods is anyway well within the standard deviation of the measurements, for the analyzed time series (see Figure 18). The power law exponent, , is derived for each time step i based on the mast anemometer measurement, available at 90 m and 57 m above ground: The blade position is discretized per time step in a Cartesian coordinate system whose origin is centered at the bottom of the tower (x is parallel to the ground, and y is parallel to the tower): Finally, the wind inflow is reconstructed for each time step using the y i position of the blade along with the i exponent.
Thanks to Equation (22), the wind inflow seen by each section is now dependent of the blade azimuthal position. The inflow velocity is assumed homogeneous in the x direction, so that the power law applies across the rotor width. A similar solution was implemented in a BEM codes comparison. 46 Figure 19 shows the iBEM results when the vertical wind gradient correction is taken into account and consequently the combined effect of all previously introduced selection methods (ie, vertical gradient filter, neighboring wake filter, yaw misalignment filter). The postprocessed results (pink circles) tend to decrease slightly the C L . The raw data (blue square) are co-plotted for reference.

| Wind inflow location correction
As seen in Figure 2, the met mast is far from the turbine, it is installed approximately 300 m away. It leads to a discrepancy between the data measured by the mast at t 0 and the data measured by the turbine at t 0 . The wind measured by the met mast reaches the turbine after a convective time which certainly depends on the atmosphere state. However, a simple model 47 was found to significantly improve results. Based on a time shift, Δt, calculated using the distance between the mast and the turbine and the angle between the mean wind direction and the direction formed by the segment "mastturbine." The following time shift is then applied to the met mast data as well as the turbine data: where l is the distance between the mast and the turbine, is the horizontal angle difference between the mean wind speed direction and the orientation "mast-turbine" and V w mean is the mean wind speed at hub height over the time series. Applying this time lag correction to the previous case (all filters and vertical gradient correction used) does not show a significant impact on the C L obtained until now (see Figure 20). However, it seems to reduce the standard deviation. Also, in order to assess the pertinence of the convective time correction, the measured electrical power is analyzed. At low wind speed, the expected power curve evolution is to increase with an increasing wind speed. However, the raw data outcome (blue curve) in Figure 21   shows the opposite: the measured power is decreasing with wind speed. When including the convective time correction, the slope is reversed and thus closer to a more expected evolution. This clearly demonstrates the importance of the time lag correction.

| Unsteady Inverse Blade Element
Momentum (UiBEM) The application of steady BEM method did not yield satisfactory results, and the C L trend remains opposite to the expected behavior regardless of the filter or correction applied. In the linear region, the lift coefficient decreases while the angle of attack increases, all the while reducing the number of samples available (see Figure 26A). The filters and corrections developed, despite ensuring the validity of the BEM equations, do not account for the turbulent conditions and stochastic nature of the environment. For these reasons, the uses of the unsteady formulation seem evident. As shown in Figure 4, the UiBEM relies initially on steady iBEM computations. The essence of the Unsteady Inverse Blade Element Momentum (UiBEM) does not differ from the steady BEM, and the aim was to determine the local velocity accounting for unsteadiness. The instantaneous local velocity is now considered as a two coordinates vector rather than a scalar as seen in Equation (14) and Equation (24). The formulation now includes: aeroelastic deformation ( ��� ⃗ V b ), tower shadow effect ( �����⃗ V w ts ), yaw and cone angles ( ����� ⃗ V rot ), and dynamic wake ( ��� ⃗ W). 34 The tilt angle has been ignored as deemed insignificant. 13 This new formulation for the velocity yields an updated velocity triangle (see Figure 22).
The DANAERO blade was equipped with tri-dimensional accelerometers (see Table 1) positioned along the radius. By applying the Equation (25), it is possible to use a first approximation of the relative blade velocity based on its local acceleration.
where V b,y and V b,z are the blade local velocity; a r,y and a r,z are the blade local accelerations; and f is the sampling frequency.
Tower shadow effect is the alteration in uniform inflow due to the presence of the tower, using potential flow theory Das et al 35 modeled the tower impact as function of the blade azimuthal position and the blade radius analyzed (see Equation 9). The dynamic wake modeling has been developed by Glauert, to account for the time delay before reaching wake equilibrium introduced by the wake deflection behind the rotor in the local velocity definition. 49 Finally, when applying the wind inflow location correction (time lag correction), the vertical gradient correction and the unsteady BEM equations we obtain the C L ( ) evolution as seen in Figure 23. The impact is significant, and the overall trend and C L values follow now the wind tunnel measurements. It is interesting to note that the computed C L is falling between the tripped and clean 3D corrected wind tunnel data.

| Tangential force measurement correction
The drag measurement for thin, high lift aerofoils is already a complex task in wind tunnels. It is near impossible for any profile in an open field. In this field measurement, the viscous drag, a part of the drag force F D , is left out because of the use of pressure sensors. The tangential force measurement, F t , is also impacting the lift force F L through the projection of the rotor forces onto the local aerofoil reference frame (see Figure 5). The relationship is made through the inflow angle , and looking at Equation (26), it is then apparent that a tangential force measurement error will influence both the lift and drag. Albeit, the C L is less impacted than the C D at low values of inflow angle ( ) thanks to the C t sin in the C L transfer formulae: (26) C L = C t sin( ) + C n cos( ) C D = −C t cos( ) + C n sin( ) F I G U R E 1 2 Binned averaged C L comparison between raw time series and gradient selection for the radial location r∕R = 49%. The square represents the average coefficients for the raw data and the circle represents the average coefficients when the mean vertical gradient filter is used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series F I G U R E 1 3 Wake calculated for the extreme wind directions measured by the met mast (in red) in the selected period. The reference wind turbine (in blue) is partially subject to a wake from WT2 | 1491 POTENTIER ET al.
In order to estimate the impact on the results due to tangential force uncertainty, a 2D CFD unsteady simulation using an in-house solver ISIS-CFD 50 was run. ISIS-CFD, has been developed by Centrale Nantes and CNRS and available as a part of the FINE™/Marine computing suite, is an incompressible Unsteady Reynolds-Averaged Navier-Stokes (URANS) method. The solver is based on the finite volume method to build the spatial discretization of the transport equations. The unstructured discretization scheme is facebased, which means that cells with an arbitrary number of arbitrarily shaped faces are accepted. A second order backward difference scheme is used to discretize time. The solver can simulate both steady and unsteady flows. The velocity field is obtained from the momentum conservation equations and the F I G U R E 1 4 Wake interaction in the wind farm for the chosen time series F I G U R E 1 5 Binned averaged C L comparison between raw time series and wake effect selection for the radial location r∕R = 49%. The square represents the average coefficients for the raw data and the triangle represents the average coefficients when the neighboring wake filter is used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series pressure field is extracted from the mass equation constraint, or continuity equation, transformed into a pressure equation. In the case of turbulent flows, transport equations for the variables in the turbulence model are added to the discretization. The turbulence model used is SST k − , 51 and the flow characteristics are representing the air at sea level at a temperature of 15°C, ie: = 1.81 × 10 − 5 (dynamic viscosity [kg m −1 s −1 ]) and = 1.225 (air density[kg m −3 ]). Regarding the simulation boundary conditions, the aerofoil related surfaces were described as "no slip wall," which means that viscous forces on the aerofoil are taken into account. The free stream velocity condition was imposed on the Inlet, Outlet and Top/Bottom, and the Outlet was using the condition Prescribed pressure. Finally, y + = 0.15 was imposed on the aerofoil surfaces.
The CFD simulations allowed to define the ratio between viscous drag (missing) and pressure drag (measured). The ratio changed between low and high angle of attack, and indeed, at low angle of attack the viscous effect is dominant (65% viscous and 35% pressure) while at high angle of attack the pressure part is dominant (29% viscous and 71% pressure). For conservatism, the results at = 2 • were used, that is, the tangential forces measured were increased by 65% to simulate the viscous drag effect. As we can see in Figure 24, the impact of the missing viscous drag measurement on the results is negligible. The added viscous drag changes neither the overall trend nor the absolute values, thus comforting the current findings.

| Discussion
Upon looking at the results presented from  it is apparent that the steady BEM equations are not sufficient to model accurately the expected behavior of a wind turbine operating in the field. Figure 25 shows the comparison of all different cases with the corrections and filters developed in Sections 4.2.1-4.2.5 and the F I G U R E 1 6 Nacelle alignment discrepancy F I G U R E 1 7 Binned averaged C L comparison between raw time series and yaw error selection for the radial location r∕R = 49%. The square represents the average coefficient for the raw data and the diamond represents the average coefficient when the yaw misalignment filter is used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series | 1493 POTENTIER ET al.
reference wind tunnel data. Nevertheless, to ensure remaining within the BEM theory limits, the application of the neighboring wake filter (Section 4.2.2) is necessary. To account for more realistic inflow conditions, the mean vertical gradient velocity correction (Section 4.2.4) and the wind inflow location correction (Section 4.2.5) should be used.
The use of the UiBEM method highlights the importance of the unsteady part in a field data analysis. First, in a physical analysis point of view, the addition of dynamic wake changes the C L trend outcome thanks to the more accurate description of the inflow velocity. The lift coefficient evolution follows now the trend measured in a wind tunnel rather than the flat or even decreasing C L ( ) relationship seen when F I G U R E 1 8 Wind profiles based on measured values F I G U R E 1 9 Binned averaged C L comparison between the raw time series and the vertical gradient velocity correction for the radial location r∕R = 49%. The square represents the average coefficient for the raw data and the circle represent the average coefficient when the vertical gradient velocity correction is used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series using the steady BEM equations. We can also see that some rotational effects are present since the mean results hovers between the 2D and 3D corrected lift curves. The error bars show the dispersion in the results and most probably show the effect of transition location displacement, or the intermittent apparition of the flow separation, both phenomena missing in the state-of-the-art BEM models. However, available measurements are limited to explain further this dispersion, as instantaneous angles of attack and wind inflow conditions must be known locally in front of the studied aerofoil for that purpose. Second, in a statistical point of view, including the unsteady equations improves the confidence in the results. Thanks to the tower shadow model and the yaw model a larger set of samples can be kept compared with the steady BEM equations. Figure 26A shows the number of samples converged using iBEM solver with the progressive application of filters and different corrections (from case 0 to case 5) at two radial locations. The addition of the transparent and opaque shaded area represents the total available number of samples to be analyzed after the cumulative application of filters and correction models. The transparent shaded area is the nonconverged part of the total number of samples. The converged number of samples decreases with the application of correction and filters (down to 10 000 samples for r∕R = 49% ). It also highlights that the number of input samples and converged samples after the application of the unsteady BEM equations increase drastically. For the tip section, the number of converged samples, when using the UiBEM method, is even greater than the raw analysis (see Figure 26B). For both sections, when using UiBEM, a maximum of 72% of the available dataset is analyzed because of wake effect between neighboring turbines; out of which 85% have achieved convergence. Figure 27 shows the results after applying the UiBEM method to all the sections. The inboard sections show a good agreement between the tripped wind tunnel conditions with the application of the 3D correction model and the calculated UiBEM C L . The agreement is poorer when it comes to the F I G U R E 2 0 Average C L comparison between the raw time series and the azimuthal wind profile correction for the radial location r∕R = 49% . The square represents the average coefficient for the raw data and the hollow diamonds represent the average coefficient when the convective time correction is used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series outer sections, and it may be attributed to the weakness of the BEM theory to capture both the tip vortex and the aeroelastic effects at those locations and needs dedicated investigations.
Previous data were presented using an angle of attack binning method. The blade azimuthal position impact was therefore compensated. In the following, the lift hysteresis is now investigated. The rapid variation of angle of attack during the blade rotation does not enable the flow around the aerofoil to reach a steady state organization. In particular, dynamic stall can locally occur and may cause deviation from the mean C L ( ) curve. Figure 28 shows the C L ( ) evolution when the data are averaged by azimuthal position (ΔΦ = 5 • ). The dynamic stall behavior is modeled using the adapted Beddoes-Leishman model from DTU. 52 Some parameters need to be chosen for this model: the reduced frequency k has been defined by the blade rotational speed (since the pitch is constant throughout the time series), the mean angle of attack m is extracted from the UiBEM simulation of the considered time series, while the amplitude of the angle of attack oscillations a is the difference between the highest quadrant averaged angle of attack and lowest quadrant averaged angle of attack.
The quadrant average results follow the overall shape described by the instantaneous cloud of points. In Figure 28, and the numbers written on the plot show some azimuthal blade position (0° being the instrumented blade pointing upwards). We can clearly see the effect of the ABL vertical gradient, which decreases the angle of attack and the C L when the blade is close to the ground (between 90° to 270°). Then, the tower shadow effect is visible between 90° and 180°, because when the blade goes down, the angle of attack does not linearly change with the blade azimuth but loops back. Despite using a crude assumption for the reduced frequency, the dynamic stall model captures the essence of both the instantaneous data and the quadrant average data. The remaining dispersion of points is certainly to be attributed to the missing knowledge of the instantaneous wind inflow at the blade location. It may be different than the wind inflow at the mast location used in the present study and needs dedicated investigations.

| CONCLUSION
The field database of wind turbine blade aerodynamic data from the DANAERO project have been analyzed further in the present paper. The aim was to identify blade element momentum theory (BEM) hypotheses of interest, to match field test measurements with simulations in idealized conditions (from 2D blade sections in wind tunnels or CFD solvers). This was performed using the back computation of the aerodynamic forces from the field measurements through the application of the inverse BEM (iBEM) and the Unsteady iBEM (UiBEM) method. The solver was first developed and verified against a CFD computation dataset. Then, several filters and correction models were implemented, whose effects on the BEM hypotheses were analyzed. The main outcome of F I G U R E 2 3 Binned averaged C L comparison for the radial location r∕R = 49%. The square represents the average coefficient for the raw data and the crossed circles represent the average coefficients when the unsteady BEM equations are used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series F I G U R E 2 2 Updated velocity triangle 34 showing the new local velocity formulation. The blade deflection velocity is not pictured the paper is that applying only the steady BEM formulation does not proves sufficient enough to capture the dynamic behavior of the turbine in the field despite the filters developed. In order to calculate the mean C L from field measurements, the use of UiBEM is mandatory since it models the dynamic nature of the inflow and turbine response.
At last, applying a dynamic stall model to the 3D corrected wind tunnel data, with an input oscillation angle associated to the blade rotation in the atmospheric boundary layer vertical gradient, proved to capture the trend of the UiBEM calculated C L azimuth averaged. The instantaneous description of the lift, mostly associated to the different F I G U R E 2 4 Binned averaged C L comparison for the radial location r∕R = 49%. The crossed circles represent the average coefficients when the unsteady BEM equations are used, and the crossed squares represent the average coefficients when the tangential force measurement are increased by 65%. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series F I G U R E 2 5 Binned averaged C L comparison for the radial location r∕R = 49%. All the previously filters, corrections, and BEM methods are overimposed. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series | 1497 POTENTIER ET al.
aerodynamic boundary layer states (eg, laminar to turbulent transition, dynamic flow separation/reattachment) under fluctuations of the atmospheric boundary layer inflow, needs dedicated local measurements with in particular the knowledge of the wind inflow and the angle of attack at each instant and in front of the blade section of interest.

F I G U R E 2 6
Effect on the data points converged of the filters, the corrections, and unsteady BEM equations for different radial locations. The transparent shaded area represents the available data points after application of the filters, and the opaque area represents the converged number of samples through the iterative iBEM procedure F I G U R E 2 7 Binned averaged C L comparison for all radial location. The square represents the average coefficient for the raw data and the crossed circles represent the average coefficients when the unsteady BEM equations are used. The black crosses show the clean 2D wind tunnel data, and the gray crosses show the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series. The red diamonds show the 3D correction model applied to the clean 2D wind tunnel data, and the green diamonds show the 3D correction model applied to the tripped 2D wind tunnel data. The error bars represent the C L standard deviation in the analyzed time series

ACKNOWLEDGMENTS
We would like to thank the IEA Task 29 members for their time and advice during meetings and emails exchanges. A special thanks to: Gerard Schepers, Helge Madsen, Koen Boormsa, and Galih Bangga for their precious help.

NOMENCLATURE B number of blades [-]
normal projection of C L and C D in the rotor rotational plane [-] C t tangential projection of C L and C D in the rotor rotational plane [ F I G U R E 2 8 Quadrant averaged C L comparison for the radial location r∕R = 49%. The green crossed circles represent the instantaneous coefficients when the unsteady BEM equations are used. The orange crossed circles represent the quadrant coefficients when the unsteady BEM equations are used. The error bars represent the C L standard deviation in the analyzed time series. The black crosses show the tripped 2D wind tunnel data. The red crosses show the 3D correction model applied to the tripped wind tunnel data. The solid green line represents the dynamic data using the tripped 3D corrected wind tunnel data as input