Accurate load prediction by BEM with airfoil data from 3D RANS simulations

In this paper, two methods for the extraction of airfoil coefficients from 3D CFD simulations of a wind turbine rotor are investigated, and these coefficients are used to improve the load prediction of a BEM code. The coefficients are extracted from a number of steady RANS simulations, using either averaging of velocities in annular sections, or an inverse BEM approach for determination of the induction factors in the rotor plane. It is shown that these 3D rotor polars are able to capture the rotational augmentation at the inner part of the blade as well as the load reduction by 3D effects close to the blade tip. They are used as input to a simple BEM code and the results of this BEM with 3D rotor polars are compared to the predictions of BEM with 2D airfoil coefficients plus common empirical corrections for stall delay and tip loss. While BEM with 2D airfoil coefficients produces a very different radial distribution of loads than the RANS simulation, the BEM with 3D rotor polars manages to reproduce the loads from RANS very accurately for a variety of load cases, as long as the blade pitch angle is not too different from the cases from which the polars were extracted.


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 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 3D spanwise effects. It is generally agreed that this is not a good assumption on a rotating blade. In order to add 3D influence to the BEM, various correction methods for 2D polars have been proposed (e.g. [1,2,3,4]). In literature, usually a lot of attention is paid to the 3D effects at the inner part of the blade (rotational augmentation, stall delay), whereas less attention is dedicated to the region close to the blade tip. Although some authors have found considerable 3D effects in the tip region (e.g. [4,5]), and a number of tip loss models have been proposed (e.g. [6], [7]), Prandtl's tip loss factor is still widely used in current BEM codes (cf. e.g. [8]).
For the rotor investigated within this work [9], a large difference in the loads was found between BEM and RANS simulations, in particular at the outboard part of the blade. In order to improve the predictions of the BEM, lift and drag polars were extracted from several steady 3D RANS simulations of the rotor. The basic concept has already been described e.g. in [10,11]. Two methods for the evaluation of the induction factors and the angle of attack were tested and compared. The 3D rotor polars obtained by this method are then applied in a simple BEM code and the results are compared to those of a BEM simulation with 2D airfoil data.
The objective of this work is to treat BEM as a kind of reduced order modeling, in the sense that it is tried to reproduce results from RANS by a drastically simpler and faster model. 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. 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.

Basic relations
We want to obtain airfoil coefficients from steady 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 the rotation frequency, but not on the wind velocity. One steady simulation per wind velocity (or rotational frequency, respectively) is required. Each simulation produces one angle of attack at each position of the blade and hence produces one point in the lift and drag polars.
The 3D rotor polars are evaluated from steady 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 obviously not 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. [12]). From a RANS simulation, the surface forces in the axial and tangential direction F ax and F t at each blade element are known. These forces could be normalized in order to obtain the normal and tangential 2 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 l and C d , defined by (where L and D are the lift and the drag force), a rotation by the inflow angle φ is required (cf. fig. 2): 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 (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 hints on how to do this can be found in literature. Shen et al. [13] obtained the induced velocities and local angles of attack from the pressure distribution or the bound circulation around the blade. Guntur and Sørensen [10] 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, we use two different methods 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 [10]; this had already been described before e.g. by Johansen and Sørensen [11]). The second one is called the inverse BEM method (method 1 from [10], already described in principle by Lindenburg [14]).

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 [11]. 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.  Johansen and Sørensen [11] neglected the tangential induction factor because they found it to have little effect on the result. In the present work, we use the tangential induction factor as well, as it takes considerable values in the inboard sections of the blade (cf. fig. 3b). In addition to what was proposed in [10] and [11], we tried to keep a maximum amount of consistency with 1D momentum theory and extended the approach such that the mass flow is conserved along the axial direction, which requires that the annular sections extend in a way that their borders follow the streamlines. In order to achieve this, a reference mass flow is calculated at the first slice  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.  4a), whereas the velocity in the wake is more uniform at lower tip speed ratios ( fig. 4b). Figure  3b shows the axial and tangential induction factors, calculated from the averaged velocities by eq. (2), as a function of the axial coordinate for the same case as shown in fig. 3a.
From the averaged velocities, the airfoil coefficients C l and C d and the angle of attack α are calculated by the following steps: (i) 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 (ii) Calculate the axial and tangential induction factors a and a by their definition eq. (2) (iii) 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) (iv) Calculate the local inflow angle φ by eq. (5) and the local angle of attack α by eq. (6) (v) Compute the local inflow velocity V rel by eq. (3) (vi) Compute the normal and tangential force coefficients C n and C t by eq. (1) (vii) Obtain the lift and drag coefficients C l and C d from eq. (4)

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 is very similar to the description by Guntur and Sørensen [10]. The inverse BEM can be described by the following algorithm: (i) Initialize a and a , e.g. a = a = 0 (ii) Compute the local inflow angle φ by eq. (5) (iii) Compute the local inflow velocity V rel by eq. (3) (iv) 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 (v) Calculate the induction factors a and a by the same equation as used in BEM. If we take the classical BEM algorithm as described by Hansen [12], we need to use where F is Prandtl's tip loss factor (cf. e.g. [12] or [6]) and σ = cB 2πr is the so-called solidity of the rotor with B the number of blades, c the chord length of the airfoil and r the radial position of the blade element. This equation is 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. [12] (vi) If a and a have changed more than a certain tolerance, go back to step (ii), else, continue (vii) 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 Glauert correction for high tip speed ratios is applied in step (v) of the extraction. This is essential for the inverse BEM as well as for the BEM itself. Glauert correction accounts for the fact that the 1D momentum theory, on which the BEM is basically built (cf. e.g. [12], [15]) breaks down on a heavily loaded rotor (i.e. operating at high tip speed ratio λ = ΩR V 0 , with R the rotor radius). When not using Glauert correction, the axial induction factor approaches unity at high tip speed ratios, which leads to unphysical 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 force 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. On the other hand, the inverse BEM method provides less physical insight than the annular sections method (for instance, plots as shown in fig. 3 can only be obtained by the annular sections method).

3D rotor polars
By plotting the local force coefficients C l and C d versus the local angle of attack α, polar diagrams 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 5, which will be discussed below, shows an example. In principle, any combination of C n , C t , C l or C d with either α or φ is possible, but we think that C l and C d vs. α is the most meaningful description for a 3D rotor polar as it describes what happens locally at the airfoil. Using C n and C t instead of C l and C d and φ instead of α leads to trouble when the rotor blades are pitched.

The wind turbine
RANS simulations were performed for the rotor of the IWT-7.5-164 wind turbine [9], which was developed in the German project "Smart Blades". It is an IEC-IA turbine with a rotor diameter of 164 m and a rated power of 7.5 MW. The blade has quite a large pre-bent (4.5 m at the tip) and a blunt trailing edge at the inner part of the blade.

RANS setup
For the steady 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.
Only the rotor has been simulated (including hub, but without tower and nacelle) in a uniform, constant, axial wind field. All simulations were performed on a rigid grid with about 16 million cells, but the elastic deformations during operation have been accounted for approximately by static grid deformation using deflection data from a previous aeroelastic simulation.

Polar extraction from RANS results
In order to obtain simulations with different local angles of attack on the blade, the rotation frequency of the rotor was fixed at 10 RPM and the wind speed was varied from 7 to 13 m/s in steps of 1 m/s (rated conditions are 10 RPM and 11 m/s). The blade pitch was zero, unless stated otherwise. 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. Figure 5 shows the lift and drag coefficients as a function of the angle of attack for three blade elements, one in the inboard part of the blade (r/R = 0.23), one in the central part (r/R = 0.66) and one close to the tip (r/R = 0.95). Each symbol corresponds to one RANS simulation at one wind speed. For comparison, we also show 2D polars which were created by two-dimensional RANS simulations in TAU for the single airfoils of which the blade consists, and in addition, for r/R = 0.23 the 2D polar plus rotational augmentation as predicted by the models of Du and Selig [1] for lift and Eggers [2] for drag are shown. In the inboard sections, the increased lift and drag due to rotational augmentation as well as a delayed stall is clearly visible, although it should be noted that a steady-state RANS simulation is not necessarily able to accurately capture the highly unsteady three-dimensional flow at the root of the blade. In the central section, the difference between 2D and 3D polar is small. In the outboard section, the lift curve has a peak close to the rated conditions of the turbine, which indicates flow separation at higher wind speeds. Indeed, when looking closer into the RANS solution, one can find local flow separation with a small vortex system parallel to the trailing edge at the tip for wind speeds around and above rated velocity (11 m/s). Drag is considerably larger at the outer sections in the 3D than in the 2D polars.

The BEM code
A simple BEM algorithm as described by Hansen [12], implemented as a Python script, was used for all BEM simulations within this work. Three common empirical corrections were accounted for: Glauert correction in order to reduce the induction factor on a heavily loaded rotor (high tip speed ratio), Prandtl's tip loss factor in order to account for 3D effects in the vicinity of the tip (downwash by tip vortex), and the stall delay model by Du and Selig [1] for lift and of Eggers [2] for drag are implemented as empirical polar correction models. The out-of-plane blade deflection was not considered (neither was it considered during the evaluation of the 3D rotor polars from RANS results). The airfoil data were supplied as tables for C l and C d vs. angle of attack at the center of each blade element. 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 severely). The 2D polars were calculated for the seven different airfoils of which the blade consists, and were linearly interpolated to the blade element centers prior to application in the BEM code. Figures 6 and 7 show the distributions of the local axial and tangential forces over the radius for different load cases. BEM-2D denotes the result of the BEM with the airfoil coefficients from two-dimensional RANS simulations plus the empirical corrections as mentioned in the previous section. BEM-3D denotes the BEM with 3D rotor polars with no empirical corrections except for the Glauert correction, which is an essential ingredient of any kind of BEM at high tip speed ratios. The RANS curve was directly evaluated from the RANS solution using the same blade element definition as in the BEM. Figure 6 includes the results of BEM-3D with both polar extraction methods (annular sections and inverse BEM). The load case corresponds to the rated conditions of the rotor (tip speed ratio λ = 7.8, wind 11 m/s, rotation with 10 RPM). In the axial force, there is only a small difference between the two methods, and both are very close to the RANS result. The tangential force shows a significant offset from the RANS values for the annular sections method, although For instance, the curves are almost identical when the reference geometry (undeformed shape of the rotor) is used, whereas there is a significant offset if the loaded and elastically deformed rotor geometry is used (as shown in the figures). For some load cases, there is also a significant offset in the axial force. The reason 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 for all load cases.  6 shows a case which has been used for the extraction of the 3D rotor polars. For validation, it is more interesting to look at cases which have not been used to obtain the polars. Therefore, figure 7a shows a case with 9 m/s wind speed and a reduced rotation frequency of 8.17 RPM (λ = 7.8). In figure 7, BEM-3D uses the polars from inverse BEM. Obviously, the force from the RANS is still quite accurately reproduced by BEM-3D. Figure 7b shows a wind speed of 13 m/s and 10 RPM (λ = 6.6). In contrast to the other figures, the blade is pitched by an angle of 8.33 • . The deviation of BEM-3D with the rotor polars obtained at zero pitch (orange star symbols) to RANS is large. In addition, the BEM-3D result with polars extracted from simulations at a pitch angle of 10 • is shown (red cross symbols), which meets the RANS result quite accurately. A closer analysis of the influence of the pitch angle on the 3D rotor polars is provided in the next section.

The influence of rotor pitch angle on the 3D rotor polars
One reason for the deviation of the BEM-3D with 0 • pitch polars in fig. 7b is that the simulation produces angles of attack which are not contained in the 3D rotor polars, so that the BEM code has to use extrapolation (the formulas by Viterna and Corrigan [16] have been used for this purpose, although linear interpolation yields almost the same result in this case because the angle of attack is not far outside of the range of the RANS polars). Moreover, it is noted that at a pitch of 10 • , we have different values for C l and C d for the same angle of attack α than in a zero pitch configuration in some parts of the blade, i.e. we obtain different 3D rotor polars for different pitch angles. Figure  This difference between the polars is caused by a different three-dimensional flow state at the rotor (compare fig. 3a to fig. 4c): The axial velocity in the wake is significantly higher and the expansion of the stream tube is hardly visible in the pitched case ( fig. 4c). Obviously, as we have the same tip speed ratio, but different pitch, the local angles of attack on the blade are different in figure 3a than in fig. 4c. However, the same qualitative difference is found if cases which lead to similar local angles of attack are compared. Therefore, we deduce that the 3D rotor polars are dependent on the blade pitch angle. For more accurate results over a range of pitch angles, the polar extraction can be performed for cases with different pitch angles (e.g. 0 • , 10 • , 20 • ), and interpolation of the polars can be applied for pitch angles in between.

Remarks on the difference between BEM-2D and RANS
It can be seen from the 3D polar data ( fig. 5), and also from a closer look into the solution of the RANS, that the simulation shows a flow separation at the tip around and above rated conditions. This could be a specific feature of this wind turbine, and it might depend on the numerical setup or the turbulence model whether this separation occurs or not. This feature contributes to the large difference between RANS and BEM-2D in the tip region. However, it is stressed that the difference between RANS and BEM-2D is large even in cases where the flow is fully attached at the tip (e.g. fig. 7a). The same analysis as discussed above for the IWT-7.5-164 turbine was done for the well-established NREL-5MW turbine with a diameter of 126 m (cf. [17]) for zero pitch. The results were qualitatively the same, although the maximum difference between RANS and BEM results was smaller for the NREL turbine. The largest difference was again found close to the tip. BEM-3D was again able to considerably improve the BEM prediction when compared  to RANS. By choosing the IWT-7.5-164 turbine for presentation in this paper, we have shown that even for cases in which there is a large difference between RANS and BEM results, we are able to bring the BEM prediction very close to the RANS results by using 3D rotor polars. The feasibility of using 3D rotor polars in the BEM for steady load prediction was shown.

Summary and Conclusions
3D rotor polars were extracted from a number of steady 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. Three main findings could be deduced: • The BEM-3D was able to reproduce the RANS results very accurately, even for load cases which were not used for polar extraction, as long as the pitch angle was not too different from the setting in which the polars were extracted. In particular, the local loads in the outer part of the blade are quite accurately reproduced by the BEM-3D, whereas BEM-2D shows very large deviations from RANS in the tip region. • Two methods for obtaining the induction factors during polar extraction were investigated.
The inverse BEM method was more reliable as it was able to reproduce all RANS results used for polar extraction almost exactly. The annular sections method sometimes produced an offset from the RANS values, which seems to be caused by an inconsistency to the BEM in the calculation of the induction factors. The inverse BEM method, in contrast, can be made completely consistent to the BEM algorithm, and furthermore, it is faster and easier to use. Therefore, it is recommended to prefer the inverse BEM method for such an application. Despite its worse performance, the annular sections method is valuable because it computes circumferentially averaged quantities which can be used to visually assess the flow situation around the rotor and provide physical insight. • 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.
Current research is focused on concepts for the improvement of unsteady BEM simulations by using high-fidelity unsteady RANS data.