Interactive comment on “ Updating BEM models with 3 D rotor CFD data

The manuscript is focused on the extraction of three-dimensional airfoil data from blade-resolved RANS analyses of an in-house code. The authors have conducted a great deal of work with respect to the number of CFD RANS analyses, extraction of airfoil tables from RANS results, and comparison of sectional blade forces among results obtained by RANS, inverse BEM, and standard BEM methods. The primary conclusion of the work is that airfoil tables (that account for 3D effects) appear to be dependent on the blade pitch angle when extracted with an inverse BEM methodology that does not include any tip or inboard stall correction. In comparison to a TORQUE paper presented by the authors, the present manuscript also includes an introductory treatment of airfoil data extraction for unsteady cases.


Introduction
Despite the availability of high-fidelity CFD solvers, the Blade Element Momentum method (BEM) is still the state of the art in wind turbine design.In particular in the certification process, where hundreds or thousands of simulations are required, the use of CFD simulations would be too expensive with respect to time and computational resources.Instead, it is often tried to improve the BEM by empirical models in order to obtain more realistic results.BEM usually uses two-dimensional airfoil data, which were obtained for cross sections of the blade by wind tunnel experiments or numerical simulations.It assumes that the aerodynamic behaviour of each blade element can be completely described by these 2D airfoil data.It does not account for any three-dimensional spanwise effects.In order to add 3D influence to the BEM, various correction methods for 2D airfoil coefficients have been proposed: For a better description of the rotational augmentation of aerodynamic forces and delayed stall close the the blade root, models were developed e.g. by Du and Selig (1997), Chaviaropoulos and Hansen (2000) and Dowler and Schmitz (2014).For the decrease of aerodynamic forces close to the blade tip (tip loss), the traditional model by Prandtl/Glauert (described e.g. in Hansen (2008)) is still widely used, although improved and new tip loss model have been proposed (e.g.Shen et al. (2005) or Sørensen et al. (2016)).Bak et al. (2006) present a more general approach which corrects 2D airfoil data based on semi-empirical relations for the pressure distribution.Within this work, a set of airfoil coefficients is determined from a number of steady-state three-dimensional RANS simulations of a rotor.This requires the local angle of attack and the lift and drag coefficients to be extracted from RANS solutions.
Concepts to do so have already been described e.g. by Johansen and Sørensen (2004) or Guntur and Sørensen (2014).In the present work, two methods for the evaluation of the induction factors and the local angles of attack at the blade were tested and compared.The airfoil coefficients obtained by this procedure will be called 3D rotor polars.These polars are applied for steady-state aerodynamic simulations in a simple BEM code and the results are compared to those of BEM simulations with 2D airfoil data.
Figure 1 shows the traditional concept of BEM with two-dimensional airfoil coefficients as input data (possibly somehow corrected by empirical models; this approach will be called BEM-2D), and, in contrast, the BEM with 3D rotor polars extracted from RANS simulations as input, which will be called BEM-3D.
The objective of this work is to treat BEM as a kind of reduced order modelling, in the sense that it is tried to reproduce results from RANS simulations by a drastically simpler and faster model (the BEM-3D).This means that the issue of validation of the RANS solution is not addressed.Instead, the results of the BEM with 3D rotor polars are validated against a number of RANS simulations in the sense of a code-to-code validation.
A similar work has already been presented at "The Science of Making Torque from Wind" (TORQUE) 2016 in Munich (Schneider et al., 2016).In contrast to the original paper, the present work uses the more renowned reference wind turbine from the European project INNWIND.EU, developed at DTU, with a rotor diameter of about 178 meters and a rated power of 10 MW (cf.Bak et al. (2013)).In addition, section 4 contains some new material in which the difference between steady states is evaluated in order to assess the accuracy of the slope of the 3D rotor polars and their potential for application in unsteady BEM simulations.
2 Extraction of 3D rotor polars from RANS simulations and their use in BEM

Basic relations
The aim is to obtain airfoil coefficients from steady-state 3D RANS simulations of a rotor.These 3D rotor polars are meant to be used in a BEM code in order to substitute the commonly used polars obtained from non-rotating 2D wings in a wind tunnel or from 2D airfoil simulations.In order to obtain lift and drag coefficients over a range of angles of attack, it is possible to either fix the rotational frequency of the rotor and vary the wind velocity, or vice versa.Note, however, that these two possibilities are not strictly equivalent, as the centrifugal and Coriolis forces on the blades depend on rotation frequency, but not on wind velocity.One steady-state simulation per wind velocity (or rotational frequency, respectively) produces one angle of attack at each position of the blade and hence one point on the lift and drag polars for each blade element.
The 3D rotor polars are evaluated from steady-state simulations of an unyawed rotor in a uniform wind field.During the polar extraction, the same assumptions as in the BEM are used: The rotor is assumed to be located in one plane which is perpendicular to the inflow, the rotor blades are discretized into a number of blade elements in the radial direction, and in each blade element, the chord and twist are assumed to be constant.These assumptions are not strictly valid for a realistic 3D geometry, however, the same assumptions will be made later in the BEM in which the airfoil coefficients will be applied.
First of all, some basic relations from BEM theory are recalled (cf.e.g.Hansen (2008)).From a RANS simulation, the surface forces in the axial and tangential direction F ax and F t on each blade element are known.These forces could be normalized in order to obtain the normal and tangential force coefficients C n and C t : where V rel is the relative inflow velocity at the blade element, is mass density, c is the chord length and ∆r the length of the blade element in radial (spanwise) direction.
However, from a RANS simulation, V rel is not directly known.Important quantities in the classical BEM theory are the axial and tangential induction factors a and a : where V 0 is the wind velocity, Ω the angular frequency of the rotor, r the radial position of the blade element and V ax , V t are the axial and tangential velocities in the rotor plane (cf.fig.2).If the induction factors a and a were known, the inflow velocity V rel could be calculated by The coefficients C n and C t refer to the forces normal and tangential to the rotor plane.The lift and drag forces are by definition perpendicular and tangential to the inflow velocity V rel .In order to obtain the lift and drag coefficients C where the inflow angle φ, in turn, can be calculated from the induction factors by (cf.fig.2) When φ is known, the local angle of attack α can be calculated by where θ t is the twist angle and θ p the pitch angle (positive for nose-down; in fig.2, it is θ = θ t + θ p ).
In summary, if the induction factors a and a are available, then C l and C d can be obtained from eqs. (3), (1), ( 5) and (4), and the corresponding angle of attack α from eqs. ( 5) and ( 6).
The remaining task is now to determine the induction factors a and a from the RANS simulation.Some methods for that purpose can be found in literature.Shen et al. (2009) obtained the induced velocities and local angles of attack from the pressure distribution or the bound circulation around the blade.Guntur and Sørensen (2014) compared four different methods for the evaluation of the angle of attack from 3D CFD simulations, three of which provide the induction factors as well.
In the present work, two different methods are used to obtain the induction factors: The first one is based on the evaluation of the circumferentially averaged axial and tangential velocity in annular sections in a number of slices upstream and downstream of the rotor (method 2 from Guntur and Sørensen (2014); this had already been described before e.g. by Johansen and Sørensen (2004) and is sometimes referred to as azimuthal averaging technique (AAT)).The second method is called the inverse BEM method (method 1 from Guntur and Sørensen (2014), already described in principle by Lindenburg (2003)).

The annular sections method
In this method, the circumferentially averaged axial and tangential velocities are evaluated in annular slices through the computational domain upwind and downwind of the rotor, similar as described by Johansen and Sørensen (2004)  also referred to as azimuthal averaging technique).Then the value of the induction factor in the rotor plane is interpolated for each blade element from the velocities in the last slice upstream and the first slice downstream of the rotor.
From the averaged velocities, the airfoil coefficients C l and C d and the angle of attack α are calculated by the following steps: 1. From the solution of the RANS simulation, calculate the circumferentially averaged axial and tangential velocity component V ax and V t in each annular section 2. Calculate the axial and tangential induction factors a and a by their definition eq. ( 2) 3. Interpolate the induction factors to the rotor plane from the last slice upstream and the first slice downstream of the rotor (these slices need to be sufficiently far away from the rotor for not to cut through the walls or the boundary layers of the blades) 4. Use eqs.( 5) and ( 6) to compute the angle of attack α, and eqs.( 3), ( 1) and ( 4) to compute the lift and drag coefficients C l and C d

The inverse BEM method
In this approach, the induction factors are obtained from the RANS surface forces iteratively using the same equations as in the BEM.The approach in the present work is very similar to the description by Guntur and Sørensen (2014).The inverse BEM can be described by the following algorithm: 1. Initialize a and a , e.g. a = a = 0 2. Compute the local inflow angle φ by eq. ( 5) 3. Compute the local inflow velocity V rel by eq. ( 3) 4. Calculate the normal and tangential force coefficient by eq. ( 1), where F ax and F t are the sum of all point forces in the RANS solution in the axial or tangential direction, respectively 5. Calculate the induction factors a and a by the same equation as used in BEM.If the classical BEM algorithm as described by Hansen is used (Hansen, 2008), then it is where σ = cB 2πr is the solidity of the rotor with B the number of blades, c the chord length of the airfoil, r the radial position of the blade element, and F is the tip loss factor according to Prandtl and Glauert (cf.e.g.Hansen (2008) applied only as long as a is below a certain threshold.For higher values of a, the value is artificially reduced by Glauert correction, which is an empirical relation, for which different formulations are possible, cf.e.g.Hansen (2008) 6.If a and a have changed more than a certain tolerance, go back to step (2), else, continue 7. Use eqs.( 5) and ( 6) to compute the angle of attack α, and eqs.( 3), ( 1) and ( 4) to compute the lift and drag coefficients Glauert correction for high tip speed ratios is applied in step 5 of the extraction.This is essential for the inverse BEM as well as for the BEM itself.Glauert correction (also referred to as turbulent wake state correction) accounts for the fact that the 1D momentum theory, on which the BEM method is built (cf.e.g.Hansen (2008), Leishman (2006)), breaks down on a rotor operating at high tip speed ratio λ = ΩR V0 (with R the rotor radius).Physically, the reason for this is that the rotor enters turbulent wake state (cf.e.g.Leishman (2006)).When not using Glauert correction, the axial induction factor approaches unity at high tip speed ratios, thus leading to non-physical results.Equation ( 7) is exactly the same as in the classical BEM, but it is used the other way round: The axial and tangential forces on the blade element are known from RANS, the induction factors are calculated from these forces (in order to do that, the force coefficients C n and C t need to be updated in every iteration using the current value of V rel ), and in the end, the local angle of attack α and the lift and drag coefficients C l and C d are obtained.
Using the same relation for the induction factors as in the BEM algorithm (eq.( 7)) instead of their definition (eq.( 2)) assures a maximum amount of consistency with the BEM algorithm, in which the 3D rotor polars are supposed to be used later.Another advantage of the inverse BEM method compared to the annular sections method is that only the surface forces from the RANS simulation are required.There is no need to do a post-processing of the volume solution of the RANS, as there is in the annular sections method for the calculation of the averaged axial and tangential velocities.

The wind turbine
RANS simulations were performed for the 10 MW rotor from the INNWIND.EU project, developed at DTU (Bak et al. (2013)).
Its diameter is about 178 m.The configuration without prebend, cone, tilt and structural deflections, as it is provided in the description of the turbine (Bak et al. (2013)), has been used.

RANS setup
For the steady-state RANS simulations, the fully compressible, unstructured, multigrid-accelerated, viscous Navier-Stokes solver TAU, developed at the German Aerospace Center (DLR), was applied.Turbulence was simulated with the Spalart-Allmaras one-equation turbulence model.The computational grid consists of all three blades, the spinner and the nacelle, but not the tower.It is an unstructured mesh of tetrahedra with prism layers built on the basis of a mostly structured surface mesh,  which contains 415 cells in the spanwise direction and 100 cells along the chord on the upper and the lower side.The mesh contains about 38 million cells in total.It was created with zero pitch.For different pitch angles, TAU's grid deformation mechanism has been used to pitch the blades.All simulations were conducted on a rigid mesh, and a uniform, constant, axial wind field was assumed.In order to obtain simulations with different local angles of attack on the blade, the rotation frequency of the rotor was fixed at 9.6 RPM (rated speed) and the wind speed was varied from 6 to 16 ms −1 in steps of 1 ms −1 at zero pitch (rated wind is 11.4 ms −1 ).In addition, pitch angles of 10 • and 20 • with wind velocities from 9 to 23 or 14 to 26 ms −1 , respectively, in steps of 2 ms −1 were simulated.

Polar extraction from RANS results
The evaluation was done as described in section 2, using the mass-consistent annular sections method as well as the inverse BEM method.The blade was discretized into 30 radial sections (blade elements) of equal length.The simulations for 0 • , 10 • and 20 • pitch were evaluated separately, which means that three sets of airfoil coefficients were obtained.
In the inverse BEM method, the tip loss factor F was set to unity in order to include the tip loss directly into the 3D rotor polars (this means that when these 3D rotor are applied in a BEM code, the tip loss model should be switched off as well).The critical induction value a c in the Glauert correction model (cf.step 5 in section 2.3) was set to 0.2.
For the annular sections method, the volume solution of the RANS simulation was evaluated in ParaView by Python macros in order to evaluate the induction in the annular sections upstream and downstream of the rotor.For the evaluation of the 3D rotor polars, 30 annular sections within the rotor diameter and only a few slices in the axial direction were used.The values of the axial and tangential induction factor in the rotor plane were linearly interpolated from the last upstream and the first downstream values, which were located 4 m upstream and downstream of the rotor.
The annular sections method calculates a circumferentially averaged value of the axial and tangential velocity components are obtained in a cross section through the computational domain with the azimuthal coordinate eliminated by averaging.
Figure 3 shows the averaged axial velocity in a longitudinal cut plane through the rotor (which is located at x = 0 from r = -89.16 to 89.16 m) for some of the RANS results which were used for polar extraction.The bold black lines are the bounds of the stream tube as assumed by 1D momentum theory, i.e. they are parallel to the streamlines and there is the same mass flow between those lines at every axial position (they were calculated by a modified version of the annular sections ParaView macro, which adapts the radial bounds of the annular sections in such a way that the mass flow in each blade element remains approximately constant along the axial direction).It is visible that the stream tube expands while passing through the rotor plane.The expansion is larger at high tip speed ratios (i.e.low V 0 for constant Ω, cf.fig.3(a)), and smaller at low tip speed ratios (fig.3(c)).Furthermore, these figures show that the axial velocity in the wake of the rotor is not uniform in every cross section as assumed by 1D momentum theory.It seems that the root vortex severely disturbs the velocity profile in the wake.This effect is very strong at high tip speed ratios (fig.3(a)), whereas the velocity in the wake is more uniform at lower tip speed ratios (fig.3(c)).This visualization is an additional benefit of the annular sections method, because it provides insight in the flow situation around the rotor and enables a comparison of the RANS result to the assumptions of 1D momentum theory.

3D rotor polars
By plotting the local force coefficients C l and C d versus the local angle of attack α, lift and drag polars can be obtained which have the same form as those from wind tunnel tests or from simulations for a 2D airfoil, but additionally contain the influence of the 3D rotational effects in the RANS solution.Figure 4 shows the polars obtained by inverse BEM for the rotor at zero pitch, with the rotational frequency kept fixed at 9.6 RPM (rated speed) and the wind speed varied from 6 to 16 ms −1 in steps of 1 ms −1 (each marker on the lines in figure 4  In the inboard sections, the increased lift and drag due to rotational augmentation is clearly visible.It is noted that at the three innermost sections (up to r/R = 0.1), the lift and drag polars had to be smoothed because the steady-state RANS solution showed some convergence problems which led to wiggly curves.In the central section, the difference between 2D and 3D polar is small.In the outboard section, lift is decreased by 3D effects, while drag is considerably larger than in the 2D polar.The 2D polars at r/R = 0.60 and 0.97 are very similar because the airfoil is very similar at these radial positions of the blade, however, the 3D rotor polars differ significantly due to 3D influence at the tip.
As already noted, it is also possible to obtain 3D rotor polars by keeping the wind velocity fixed and varying the rotational frequency.The polars obtained by doing so are very similar to those in figure 5, however, it is noted that the drag polars differ at lower angles of attack.In particular, the negative slope in the drag polars, which is observed at the inboard and central section in fig.4, does not appear when using simulations at constant wind speed and variable rotational frequency covering the same range of the tip speed ratio.This indicates that the actual value of the rotational frequency (and not only the tip speed ratio) does influence the drag polar.This might be explained by the centrifugal and Coriolis forces on the flow, which depend on rotational frequency, but not on wind velocity.
Prior to application in the BEM code, the 3D airfoil polars were extrapolated to the whole angle of attack range from -180 • to 180 • using the well-known extrapolation formulas by Viterna and Corrigan (1982).

The influence of pitch angle
The polar extraction has been repeated for RANS simulations with pitch angles of 10 • and 20 • (cf.section 3.2).The comparison of the 3D rotor polars obtained at different pitch angles shows that there is a dependency on the blade pitch angle.The polars at 0 • and 10 • pitch are compared in figure 5 (for the sake of clarity, the polars for 20 • are not included in the figure, but it is noted that they show the same qualitative difference to the 10 • polars than the 10 • to the 0 • polars).
At the inner section in fig. 5 (r/R = 0.23), it is obvious that different lift and drag values are obtained at the same angle of attack in the 10 • than in the 0 • case.At the outer section (r/R = 0.97), a different range of the angle of attack is obtained in the two configurations, and in the overlapping part, the values and the slope of the two polars do not match.In the middle of the blade (r/R = 0.60), the polars for the different pitch angles match rather well.This difference between the polars is caused by a totally different flow state at the rotor (compare fig. 6  fig.6 does not compare cases with equal local angles of attack (which is not possible anyway because the radial distribution of the angle of attack depends on the pitch angle as well), the different radial distribution of induction might explain the pitch dependency of the 3D rotor polars.The consequence is that the relation between the angle of attack α, the pitch angle θ and the inflow angle φ is not as easy as assumed in the BEM algorithm (cf.eq. ( 6)).For more accurate results over a range of pitch angles, the polar extraction should be performed for cases with different pitch angles, and interpolation of the polars can be applied for pitch angles in between (it will be demonstrated in section 3.8 that steps of 10 • for the pitch angle are sufficient).
It is noted that in the common BEM with 2D airfoil coefficients, the lift and drag polars are always assumed to be independent of the pitch angle.However, the findings of the present work indicate that this is wrong for 2D polars as well as for 3D rotor polars.

2D Polars
In the BEM simulations, not only the 3D rotor polars described in the previous section were applied, but also three different kinds of 2D airfoil polars for comparison.For complete consistency with the 3D rotor simulations, 2D airfoil polars were used which were obtained from two-dimensional RANS simulations in TAU using the same turbulence model (Spalart-Allmaras, fully turbulent) and -as far as possible -the same numerical setup as in the 3D RANS simulations.These 2D TAU polars were calculated for the five different airfoils of which the blade consists, including the Gurney flaps at the FFA-W3-600, -480 and -360 airfoils at the inner part of the blade (cf.Bak et al. (2013)).Just like the 3D rotor polars, the 2D TAU polars were extrapolated to the whole angle of attack range from -180 • to 180 • using the formulas by Viterna and Corrigan (1982).
Furthermore, the two sets of two-dimensional airfoil polars which DTU provides for the reference turbine were used.The first set (called 2D DTU in the figures below) are the purely two-dimensional airfoil coefficients for the five airfoils which define the blade, which were obtained by DTU with the code EllipSys2D (cf.Bak et al. (2013)).The second one is a 3D-

The BEM code
A simple BEM algorithm as described by Hansen (2008), implemented as a Python script, was used for all BEM simulations within this work.The tip loss factor by Prandtl/Glauert was used for simulations with 2D polars in order to account for the load reduction in the vicinity of the tip, but not for simulations with 3D rotor polars, because the tip loss is included directly in the 3D rotor polars.Glauert correction (or turbulent wake state correction; implemented as proposed by Hansen (2008) chapter 6) is used in order to reduce the induction factor on a heavily loaded rotor (i.e. at high tip speed ratio).This correction is necessary because of the limitations of 1D momentum theory, which fails if the rotor enters turbulent wake state (cf.Leishman (2006)).
Therefore, Glauert correction must be applied in any kind of BEM, no matter if 2D airfoil polars or 3D rotor polars are used.
The 2D airfoil polars were linearly interpolated to the radial positions of the blade element centers by the BEM code.The 3D rotor polars had been evaluated in the same blade sectioning before, so no interpolation was necessary (although interpolation of the 3D rotor polars is not a problem, as long as the element definitions in the polar extraction and in the BEM do not differ too much; if they differ significantly, then information from the 3D rotor polars may be lost due to the interpolation).can be considered a kind of verification).The other load cases (P1 to P3) have not been used for extraction of the polars and can be considered to be prediction.For comparison, additional RANS simulation were performed for these cases.Comparing the BEM predictions to the RANS results for cases P1 to P3 should therefore be considered a code-to-code validation.

Comparison between RANS, BEM-3D and BEM-2D
Figures 7 to 13 show the distributions of the local axial and tangential forces over the radius for different load cases.The curves labelled by BEM-2D DTU, BEM-2D corr DTU and BEM-2D TAU denote the results of the BEM with the 2D airfoil coefficients as described in section 3.6.BEM-3D invBEM and BEM-3D anSecs denote the BEM results obtained with 3D rotor polars extracted by the inverse BEM method and the annular sections method, respectively.The RANS curve was directly evaluated from the forces in the RANS solution using the same blade element definition as in the BEM.
Figures 7 shows load case E2, which is close to the rated conditions of the rotor (wind 11 ms −1 , rotation with 9.6 RPM, tip speed ratio λ = 8.15).In the axial force, there is a significant offset between the BEM-2D results and the RANS result.
The tangential force is predicted too low by BEM-2D DTU and BEM-2D TAU in the inboard and central parts of the blade, whereas the values are too large at the blade tip.The BEM-3D invBEM result meets the RANS result almost exactly, whereas the BEM-3D anSecs result shows a significant deviation from RANS.The same offset can be seen in the load cases E1 and E3 (at different wind speeds; figures 8 and 9).It is noted that the kinks in the tangential force in figure 9 do probably not represent actual physics.This feature from RANS, although probably not very realistic, is also reproduced by BEM-3D.
The reason for the deviation of the BEM-3D anSecs curves from RANS is most likely the inconsistency in the calculation of the induction factors between the annular sections method, which uses eq. ( 2), and the BEM, which uses eq. ( 7).With the inverse BEM method (which uses the original BEM relation eq. ( 7)), the BEM-3D meets the RANS solution almost exactly.This is true for all load cases which have been used for 3D polar extraction.It is noted that the difference between the annular sections method and the inverse BEM method is also visible in the lift and drag polars (the annular sections result has been omitted from figure 5).However, only the comparison of the BEM-3D forces with the original RANS forces, from which the 3D rotor polars were obtained, allows to judge which method produces the correct lift and drag polars (in this context, the "correct" polars are those who are able to reproduce the RANS forces from which they were obtained).Figure 10 shows the axial induction factor as a function of the radial coordinate for load case E2 (there is no RANS curve, as there is no unique way to determine the induction factor from RANS; it can only be obtained by e.g. the annular sections method or the inverse BEM method, but these results differ).It is clearly visible that the induction factor behaves differently in BEM-3D than in BEM-2D.BEM-3D does not use any tip loss model and the load reduction close to the tip is directly produced by appropriate values for the lift and drag coefficients close to the blade tip (lift is decreased, drag is increased  reduction at the tip by increasing the induction factor, which is physically wrong.Shen et al. (2005) have shown analytically that the induction factor goes to unity at the blade tip when using the Prandtl/Glauert tip loss factor (for non-zero drag), which would mean that there was zero axial velocity at the tip, which is obviously wrong and can be considered an inconsistency in the tip loss model by Prandtl/Glauert.This inconsistency is present in most tip loss models which are based on the Prandtl/Glauert model according to Shen et al. (2005).When using 3D rotor polars in BEM, the tip loss model is not required, and the BEM  For validation, it is important to look at load cases which actually represent prediction, i.e. which have not been used to obtain the polars.Therefore, figure 11 shows the tangential forces for the above-rated cases P2 (rated rotation speed 9.6 RPM, 14 ms −1 wind, pitch angle 9.292 • , tip speed ratio 6.40).The curve coloured in cyan and labelled BEM-3D invBEM 0 • was produced with the 3D rotor polars obtained by inverse BEM from RANS simulations at zero pitch.Obviously, the black RANS curve is not met at all, and the result is even worse than for BEM-2D at most of the locations.In contrast, when the 3D rotor polars obtained at 0 and 10 • pitch are both used as input to the BEM and the BEM code interpolates these polars linearly to the current pitch angle of 9.292 • (red colour and label BEM-3D invBEM in fig.11), then the RANS result is met almost exactly.
One reason for the deviation of the BEM-3D with 0 • pitch polars is that the simulation produces angles of attack which are not contained in the 3D rotor polars, so that the extrapolated range of the polars is used (the different ranges of the angle of attack can be seen in figure 5; this mostly causes the difference in the outboard half of the blade).Moreover, as can be seen from figure 5, in some parts of the blade there are different values for C l and C d for the same angle of attack at 10 • than at 0 • pitch, i.e. different 3D rotor polars are obtained for different pitch angles (which mostly causes the difference in the inboard half of the blade).Figure 12 shows the above-rated case P3 with a pitch angle of 17.618 • (wind 20 ms −1 , 9.6 RPM, λ = 4.48) with the pitch-interpolated 3D rotor polars represented by the red line (in this case, interpolated between the 3D rotor polars for 10 • and 20 • ).This case also shows very good agreement between BEM-3D and RANS.This demonstrates the necessity to account for the dependency of the 3D rotor polars on the pitch angle.The pitch angle should be considered as an additional parameter for the 3D rotor polars.Linear interpolation between 3D rotor polar sets with steps of 10 • pitch in between seems to be sufficient.
Figure 13 shows the below-rated case P1 with wind speed 9 ms −1 and 7.229 RPM (tip speed ratio λ = 7.5).It is obvious that the axial force matches quite well between RANS and BEM-3D, whereas the tangential force does not match very well.The results for the axial force match very well, as the axial force is hardly influenced by drag.This difference in drag is probably caused by the fact that the 3D rotor polar has been obtained obtained at a higher rotation speed than the RANS simulation for this load case, so the forces on the flow due to rotation (centrifugal and Coriolis forces), and hence the 3D influence, are different.This is confirmed by the observation that the difference between BEM-3D and RANS result at a given load case becomes smaller if the rotation speed in the RANS simulations from which the polars are extracted is closer to the investigated load case.This indicates that, strictly spoken, not only the pitch angle should be treated as an additional parameter for 3D rotor polars, but also the rotation speed.However, as this would drastically increase the effort required to obtain a set of 3D rotor polars, and as the BEM-3D results presented here are still a significant improvement compared to the BEM-2D results, this dependency was not further investigated.

Quasi-steady comparison
It was shown so far that the principle of 3D rotor polars is feasible in the context of steady-state BEM simulations.Of course, it is desirable to use 3D rotor polars also for unsteady simulations.Finding suitable descriptions of e.g.dynamic inflow models and dynamic stall models to describe the unsteady aerodynamic behaviour is a separate task, which is not directly connected  Figure 14 shows the result for a load case with wind 11 ms −1 , 8.836 RPM, pitch 0 • , which is very similar to load case E2 from the previous section, for a pitch change of 1 • (i.e. the difference between the local forces at 1 • and 0 • pitch is shown).
Figure 15 shows 17 ms −1 , 9.6 RPM, pitch 13.896 • with a pitch change of 1 • .The polar datasets shown in these figures are the 3D rotor polars from inverse BEM, the 3D-corrected 2D airfoil polars from DTU, and the 2D airfoil polars obtained with TAU.The figures show only the difference in the local forces between the two steady states (deformed minus undeformed).
The RANS evaluation and the BEM simulations were both performed with 20 blade elements of equal length.
The difference in the local axial forces is quite similar for the all sets of airfoil coefficients.Just like for the steady-state comparison, the BEM-3D result is closer to the RANS result, in particular close to the root and the tip of the blade, although the relative deviation tends to be larger than in the steady-state figures.
In the local tangential forces, the deviation between the curves is more pronounced.The 2D polars show implausible kinks close to the blade root.These kinks are caused by very different slopes in the interpolated 2D airfoil polars at these radial positions.The BEM-3D result, in contrast, behaves quite smoothly.In the tip region for the 17 ms −1 case (fig.15), the RANS result for F t is met quite accurately by BEM-3D and by BEM-2D, except for the outermost section, where the deviation is larger, but still acceptable.For the 11 ms −1 case (fig.14), both BEM-2D and BEM-3D show a strong deviation from the RANS result, however, the BEM-3D result is still slightly closer to RANS.As this large deviation from RANS does appear in the tangential force only, and as the tangential force is significantly influenced by drag whereas the axial force is usually not,  it is deduced that the deviation of the BEM-3D from the RANS result is mostly caused by the slope of the drag polars at this operating point.
In the 11 ms −1 case, there is also an unrealistic kink in the RANS curve, which is most likely caused by a convergence problem of the steady-state RANS simulation at this position of the blade (this kind of problem is rather common at the most inboard sections because the flow is separated at the upper side of the blade; it is questionable whether a steady-state RANS simulation can lead to accurate results under these circumstances).
On the whole, the 3D rotor polars seem to capture not only the value, but also the slope of the polars, and seem to be better suited for the description of the difference between steady states than the 2D airfoil coefficients used for comparison.
Therefore, it is deduced that 3D rotor polars can be used as the basis for unsteady modelling.The remaining task and the main 18 5 Summary and conclusions 3D rotor polars were extracted from a number of steady-state RANS simulations using two different methods for the calculation of the induction factors.These 3D rotor polars were applied in simulations with a simple BEM code without using any empirical corrections except for Glauert correction.Four main findings could be deduced: -The BEM with 3D rotor polars (BEM-3D) can produce results which are very similar to RANS results for a variety of load cases (also for many load cases which differ significantly from those used for polar extraction).
-The inverse BEM method is more reliable than the annular sections method, as the 3D rotor polars obtained by inverse BEM method are able to reproduce all RANS results used for polar extraction almost exactly.The annular sections method deviates from RANS results because of an inconsistency to the BEM algorithm in the calculation of the induction factors.The inverse BEM method, in contrast, can be made completely consistent to the BEM algorithm by using the same algorithm with the same assumptions on geometry and the same submodels (e.g.tip loss, turbulent wake correction).Furthermore, it is faster and easier to use.Therefore, it is recommended to use the inverse BEM method if airfoil coefficients are extracted from RANS simulations for the purpose of application in BEM simulations.
-The flow around the rotor blades and the structure of the rotor wake change significantly if the blades are pitched.This leads to different 3D rotor polars for different blade pitch angles.Therefore, 3D rotor polars should be extracted for multiple pitch angles.Only few different pitch angles are necessary.Steps of 10 • between the pitch angles with linear interpolation of the lift and drag polars for pitch angles in between were shown to be feasible.
-The comparison of the difference between two steady states (quasi-steady comparison), which is a way to assess the slope of the polars, showed that 3D rotor polars can be used as a basis for CFD-based unsteady modelling, although for one load case relatively large errors were found for the tangential force due to wrong slope in the drag polars.

Data availability
The data files associated with the simulations and the figures from this paper are archived at DLR. Requests for access to the data may be sent to the corresponding author.Data concerning the reference wind turbine are available from DTU via http://dtu-10mw-rwt.vindenergi.dtu.dk.

Figure 1 .
Figure1.The traditional BEM approach with 2D airfoil data (BEM-2D) and the BEM with 3D rotor polars from RANS simulations (BEM-3D) l and C d , defined by C l = L 0.5 V 2 rel c∆r and C d = D 0.5 V 2 rel c∆r (where L and D are the lift and the drag force), a rotation by the inflow

(
and induction factors) in each annular section.If the annular sections method is applied not only within the rotor diameter, Wind Energ.Sci.Discuss., doi:10.5194/wes-2016-51,2016 Manuscript under review for journal Wind Energ.Sci.Published: 20 December 2016 c Author(s) 2016.CC-BY 3.0 License.but in a certain surrounding of the rotor, and if a sufficiently fine resolution is used, then a visualization of the flow around the rotor in a circumferentially averaged sense is obtained, i.e., the values of the velocity components (or induction factors) corresponds to one RANS simulation).Polars for three blade elements are shown (lines with symbols): one in the inboard part of the blade (r/R = 0.23), one in the central part (r/R = 0.60) and one close to the tip (r/R = 0.97).Each symbol corresponds to one RANS simulation at one wind speed.For comparison, fig. 4 also contains the 2D polars which were created by two-dimensional RANS simulations in TAU for the single airfoils of which the blade consists (bold lines), and in addition, the 2D polars plus correction by the method of Bak et al. (2006) (as provided by DTU with the reference turbine) are shown (thin lines).

Figure 4 .
Figure 4. Airfoil coefficients at three radial positions, (a) lift and (b) drag coefficient vs. angle of attack, 2D airfoil polars from TAU (bold lines) vs. 3D-corrected 2D polars from DTU (thin lines) vs. 3D rotor polars obtained from RANS by the inverse BEM method (lines with symbols; each symbol represents one RANS simulation) ): If the blades are pitched (fig.6(b) and (c)), the axial induction in the wake is drastically lower than in the zero pitch case (fig.6(a)).Although

Figure 5 .
Figure 5.Comparison of the 3D rotor polars extracted at zero pitch (points only) and at 10 • collective pitch (lines with points)

Figure 12 .
Figure 12.Load case P3: Wind 20 ms −1 , 9.6 RPM, pitch 17.618 • , axial (a) and tangential (b) forces on the blade elements vs. radial coordinate (sum over all blades); case not used for polar extraction to the concept of 3D rotor polars.A tabulated set of airfoil coefficients cannot describe transient behaviour without further 16 Wind Energ.Sci.Discuss., doi:10.5194/wes-2016-51,2016 Manuscript under review for journal Wind Energ.Sci.Published: 20 December 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 14 .
Figure 14.Quasi-steady comparison: Difference between pitched and unpitched case in axial (a) and tangential (b) force, load case 11 ms −1 , 8.836 RPM, pitch 0 • , pitch deflection 1 • Wind Energ.Sci.Discuss., doi:10.5194/wes-2016-51,2016 Manuscript under review for journal Wind Energ.Sci.Published: 20 December 2016 c Author(s) 2016.CC-BY 3.0 License.effort for the description of unsteady behaviour is to characterize the unsteady transient from one steady state to another (this is usually referred to as dynamic inflow modelling in the context of BEM).The authors are currently working on CFD-based unsteady aerodynamic models for use in BEM.
).If possible, the tip loss correction should be switched of in the target BEM code and F should be set to unity; if this is not possible, then at least the same tip loss model should be used in the inverse BEM as in the BEM.This equation is c Author(s) 2016.CC-BY 3.0 License.

Table 1 .
Table 1 lists the load cases for which figures are shown within this section.This includes some of the load cases used for extraction of the 3D rotor polars (E1 to E3), where a very good agreement should be expected, as this simply means that the 3D rotor polars, when applied in BEM, reproduce the force distributions from which they were obtained by polar extraction (this Wind Energ.Sci.Discuss., doi:10.5194/wes-2016-51,2016 Manuscript under review for journal Wind Energ.Sci.Published: 20 December 2016 c Author(s) 2016.CC-BY 3.0 License.The load cases shown in the figures; upper part: Load cases used for polar extraction (only those shown in the figures); lower part: