Consistent 3D CFD and BEM simulations of a research turbine considering rotational augmentation

The present studies evaluate the importance of the three-dimensional (3D) effects in a research turbine blade by comparing the Blade Element Momentum (BEM) simulation results to the fully resolved 3D Computational Fluid Dynamics (CFD) computations. The turbine has a diameter of around 50 m producing a rated power of 750 kW and is designed to be erected at Stöttener Berg in the southern part of Germany within the framework of the WINSENT (Wind Science and Engineering in complex Terrain) project. In the first part of the studies, 3D CFD simulations employing the URANS approach are performed for several pitch angles at a constant freestream inflow condition. The effective angles of attack are evaluated using two different approaches for 6 defined sections along the blade radius to extract the 3D polar characteristics. Then, the extracted 3D polar data are applied in the BEM simulations and compared to the results employing 2D polar with a 3D correction model in the second part of the work. The results show that the use of the 3D polar data improves the BEM prediction significantly while the 2D corrected polar data generate a strong deviation compared to the 3D CFD results.


Introduction
Driven by the ambition for a power production based on renewable energy sources, the installed wind capacity is increasing constantly every year. More place is needed with the requirement to be far enough from human settlements in order to avoid disturbance. Now complex further terrains have become of interest, but their complex inflow conditions require deep investigations to determine their practical effects. At the same time, wind turbine blades are becoming not only larger, but especially lighter and this second condition is leading to a higher importance of the aeroelastic effects. Rigid body CFD simulations are not sufficient anymore because they do not take into account twist deformations and tip displacement, that are directly connected to different effective angles of attack (AoA), forces distribution and power estimation.
From these issues the WINSENT project [1] is born. Within the cooperation of different institutes in the south of Germany (Universität Stuttgart, TU Munich, Hochschule Esslingen, Karlsruher Institut für Technologie, Universität Tübingen and Zentrum für Sonnenenergie-und Wasserstoff-Forschung Baden Württemberg), the effects of aeroelasticity on two research wind turbines in complex terrain [2] will be analyzed. New numerical models at different levels will be assessed and compared within each other and with field results, that is why, in order to have a sound and fair comparison, it is necessary to ensure consistency within the models. In order to make BEM computations, a simplified blade model of the geometry is necessary (starting from the twist and chord distributions and shape of the forming airfoils). Then, polars at some defined sections need to be provided. The geometry properties are defined and unambiguous, while different polars can have different origins. For example, they can be evaluated experimentally, with Xfoil, with CFD, etc. In CFD, it needs to be distinguished between 2D and 3D simulations. Then, in the case of 2D polars (independently from their origin), these can be corrected in order to take into account the so called 3D effects. In contrast, 3D CFD does not need for any kind of correction, because of the resolution of the entire blade model in the domain.
The root area of the blade is usually subjected to separation because of the large airfoil thickness associated with strong adverse pressure gradient, while at the blade tip recirculation and vortices shedding (tip loss effect) do usually occur. In the root area the 3D effect is called the Himmelskamp effect [3,4] or rotational augmentation. This increases the local lift coefficient compared to the 2D conditions due to the Coriolis and centrifugal forces. On the other hand, the lift coefficient reduces near the tip area due to the tip loss effect. BEM usually uses 2D airfoil data adding tip loss and rootloss corrections like the Prandtl one [5], or the new tip loss models described by Shen et al. [6], and 3D corrections (see Snel et al. [7]). Analysis on the BEM performances with CFD input data and influence of the AoA extraction methods have been performed by Schneider et al. [8], Johansen et al. [9] and Guntur et al. [10]. In this work, the strategy used in order to ensure consistent numerical models of the turbine from the geometry and aerodynamic point of view is presented. Starting from the CAD model, 2D and 3D CFD simulations using the URANS solver FLOWer [11] are performed and based on that, the BEM model is defined in order to perform steady computations. Polars are extracted from 3D and 2D results (the last ones are considered both with and without corrections), in order to show the effects on the aerodynamic results of Qblade [12], a BEM solver developed at the Hermann Föttinger Institute of TU Berlin.  [11]. This was firstly developed by the DLR (German Aerospace Center) and now, since many years it is extended for wind turbine applications by the IAG. It is a finite volume URANS solver and simulations have been run using the Shear-Stress-Transport (SST) k-omega according to Menter [13], assuming fully turbulent boundary layer. The spatial discretization scheme used is a central cell-centered Jameson-Schmidt-Turkel (JST) [14] finite volume formulation. This scheme gives a second order accuracy in space on smooth meshes. An artificial hybrid 5-stage Runge-Kutta time-stepping scheme and multi-grid level 3 is used to accelerate the convergence of the simulations.

3D setup
In order to create a CFD model, the extraction of a smooth and "water proof" outer surface of the components from the provided CAD model is necessary. For these computations, a "one third model" is used in order to use the periodic characteristic of a 3bladed wind turbine and in this way, to save computational time. In fact, considering that only the loads on the blade are of interest, there is no need for a full model for the polar evaluation process. In this way, only one blade and one third of the hub and nacelle have been meshed. The blade mesh has been created using a GRIDGEN script developed at the IAG, namely Automesh [15].This is not only able to create the CH mesh for the blade, but also to extract the local chord and twist that are necessary to make comparable 2D and BEM calculations. The hub, nacelle and background meshes have been created using the mesh generator Pointwise. All the components have been meshed ensuring y + ≈1 in the boundary layer region. Two background meshes have been created in order to optimize the cells number, a finer one (≈ 6.5R long) at the blade region and a coarser one (≈ 13R long) in the outer region. The sum of all the mesh components and the backgrounds led to a total number of ≈ 22 millions cells, see Figure 1.
A mesh refinement study of the blade mesh has been performed in advance. Three different blade grid resolutions are evaluated, namely the "normal" mesh with ≈ 22 million cells, a "fine mesh" with ≈ 12 million cells, and a "coarse" one that is the "normal" one with level 2 of multi-grid. The grid resolutions for the other grid components remain unchanged. The blade is in fact the most sensitive component. Note that the meshes of the background, hub and nacelle have already fine resolutions around the rotor area. The studies are carried out at a constant wind speed of 11 m/s and a rotational speed of 26.5 rpm. Steady simulations are performed for these two cases and the resulting sectional loads are presented in Figure 2. It can be seen clearly that the "fine" and "normal" meshes show a similar prediction for the radial fraction greater than 0.2, but a slight deviation is shown in the root area. The "coarser" one, even if completely converged, shows a much higher deviation from the other two meshes. Considering the computational efficiency, the "normal" mesh is chosen for further simulations.  For the generation of the 3D polars, different CFD simulations have been run varying the blade pitch between -55 and +20 degrees, and after reaching the steady state, AoA, c l and c d have been extracted by the use of the RAV method [16] and the LineAve (or Circave) method [17]. The difference between the two methods consist mainly in the concept used to average the flow velocity around the blade. In the RAV method, two annular elements upstream and downstream equidistant from the blade are used, while in the LineAve, a symmetric closed line (in this case a circle centered in the quarter chord position) around the blade is used. In both cases the bound circulation influence is eliminated and the AoA is calculated. In [18] the use of different methods to extract the AoA including RAV and LineAve has been evaluated. It was observed that some differences between the methods can occur in the root area due to the flow separation. In fact here the bound circulation has a chordwise component that is not taken into account by methods employing only 2D sections along the blade. Anyhow really good agreement between different approaches was observed for r/R > 0.3.
The same operating conditions as the grid studies are employed but the URANS approach is used instead of RANS with a time step size of 1 • . Then, 6 sections at different radial positions have been chosen with the project partners as a first configuration to test the procedure. Due to the blockage effect, it is not possible to evaluate c l and c d for high angles of attack from CFD, that is why the polars have been extracted only in the range ±30 degrees. In order to ensure consistency between 2D and 3D simulations, the local blade operating conditions have been evaluated and used as flow setup for the 2D CFD simulations. The 3D correction scheme used in the 2D computations is the one by Chaviaropoulus [19], where the constants a, h and n are set equal to 2.2, 1.3 and 4 respectively.
In this way three polar sets are prepared: • 3D polars based on a "one third" CFD model; • 2D polars from 2D CFD with and without empirical 3D corrections; • Xfoil polars (without corrections);  [12]. This an open source software package with several usages, for example, an extension the Xfoil functionality allowing Xfoil polar calculation and extrapolation within the software package, blade design and optimization, rotor performances analysis, turbulent inflow generation according to Veers, etc. For these studies, steady BEM computations are employed to simulate the blade loads. Further detail is given in the subsequent sections.

Setup
The Qblade model of the blade has been created, starting from the chord and the twist distribution extracted from the blade CAD model. Linear interpolation is applied between the sections. It was then necessary to load the polars files for each desired section. The blade bottom is modeled as a circular foil with c d =0.5, so that the root area is a linear interpolation up to the first provided section at r/R=21.7%. Qblade, like in general BEM tools, needs the entire 360 degrees polar range, that is why extrapolation methods have been applied. Qblade allows the use of the the Montgomerie [20] and Viterna method [21], which have both as general assumption that the airfoil behaves as a flat plate at high angles of attack. The main difference between the two is the treatment of the spanwise effects. Although the extrapolation can consistently differ according to the parameter settings, considering that the range of the interested AoAs lays in the provided one within the reference blade sections, there is no relevant influence in the use of the two extrapolation methods. It is important to underline that the simulation in Qblade have been performed at a very low pitch, leading to angles of attack mostly in the provided range of AoA. Only in the root region, where the highest angles of attack can be found, can the extrapolation method have a strong influence. Viterna showed smoother results, that is why, it has been chosen and applied for all the sections.
The simulation parameters have been set, where the Prandtl tip correction [5] has been applied, in which the wind turbine is modeled through vortex sheets. Additionally, the foil interpolation option is used, that is not a correction of the BEM algorithm, but applies a linear interpolation between the polars at the provided sections. The same set up has been used for all the polar sets, in order to maintain consistency in the results.
Qblade allows to use in combination with the Prandtl Tip Loss also the Prandtl Root Loss, accounting for the vortex shedding at the blades root that results in a reduction of the AoA computing by BEM in this area. The affected area is below the first provided section, and so a correction here is not of interest for this study. Another possibility is the combination of the New Tip Loss and New Root Loss correction factors, that are the correction methods proposed  Shen [6] that change the c N and c T considering that the forces should reduce to zero moving to the tip because of pressure equalization. The relaxation factor was set equal to 0.2 and different number of inner iterations and tolerances have been tested to ensure the convergence of the results.

2D and 3D polars
As it can be expected, due to a strong radial flow component that occurs in the hub region, 3D CFD leads to a higher estimation of c l , and this effect gets smaller moving to the blade tip, see Figure 4. It was described by Bangga et al. [4,22] that the centrifugal force causing the underlying radial flow reduces the boundary layer thickness that improves the flow stability. The Coriolis force acts in the chordwise direction as an inertial force creating an additional acceleration term that delays the occurrence of separation [4,22]. Both effects were observed to be the dominant sources for the lift augmentation in the blade root area. As it can be seen, the correction of the 2D polars leads to an increase of c l , approaching the 3D CFD polars. The linear region shows differences smaller and smaller between 2D and 3D moving to the tip, but, as will be shown in the following section, the effective angles of attack over the blade radius lay mostly after the linear region. This means that also far away from the root region, a big difference in the use of different polars can be observed. The two methods RAV and LineAve have also been compared in Figure 3 where c l and c d at the section r/R=21.7% and 61.6% are compared. They show mostly a good agreement so that the difference is not enough to create visible differences in the BEM calculations.

Comparison BEM and CFD
In order to evaluate the effects of the different polars sets on the BEM based calculations, the tangential loads (F T ) and the normal loads (F N ) in function of the normalized blade radius, and the Power in function of the blade pitch are computed using the different sets of polars. These curves have been then compared with the same ones extracted by 3D CFD computations, after normalization with the maximum value along r/R and pitch angle, respectively. It can be seen in Figure 5 that the BEM prediction gets closer to the CFD results as the polar data obtained from the 3D CFD simulations are applied. Despite that, the predictions in the root area are delicate. It is shown that even though the BEM results using the 3D CFD polar have a fairly good agreement with the 3D CFD data, the tangential force is not accurately simulated especially for r/R < 0.4. This shows that the underlying rotational augmentation is complex and its influence cannot be easily captured by BEM simulations even though consistent polar data are used. The 2D polar data show an even stronger deviation compared to the 3D CFD results, indicating that their usage is limited only to the outer blade area. The applied 3D correction model improves the prediction up to r/R = 0.35 of F T but results in higher values for F N further inboard in comparison to 3D CFD results.
For r/R within 0.4 and 0.8 there is no visible improvement in the F T estimation using the One of the first explanation for the difference in F T can lay on the AoA evaluation from BEM. In fact, looking at the Figure 6, it can be seen that BEM computes lower AoAs in this region with a difference of ≈ 6%. This means that BEM is referring in this area to a smaller c l , with a difference that is increasing with the inclination of the polar in the linear region.
For r/R > 0.8, BEM using 3D CFD polars results in higher values of F T and the reason for this has been associated with the low number of sections used for the polar extraction procedure. In fact as shown in table 1, no sectional polar is provided for almost 30% of the blade length, leading to the use of interpolated polars in between.
Finally, a further problem that can lead to such a difference in the F T estimation is the induction factor a. As known, one of the limitation of BEM is that the simple momentum theory breaks down when the induction factor becomes bigger than around 0.5, because the wind velocity behind the blade becomes negative. As can be seen in Figure 6b the induction factor calculated with both CFD and BEM is already higher than 0.4. In order to overcome this problem, empirical models should be applied to correct the thrust coefficient [23] (and this is not possible in Qblade) or more accurate models like the Vortex model should be used.
Observing the power curve in reference to the pitch angles, it can be seen that there is a good agreement between BEM and 3D CFD; in particular it can be seen, as expected, that the power curve calculated from the 2D corrected polars is in between the one with 2D polars without corrections and the one with 3D polars. This shows clearly the importance of the 3D effects and in particular the impact of the use of different models to take them into account. The BEM simulation using Xfoil uncorrected polars, from the other side, results in higher values for the Power. Considering that the Xfoil polars have been used without 3D corrections, and these lead to a higher c l , the resulting power and forces estimation would be even higher and completely not fitting with the 3D CFD results.

Conclusions
In this paper the influence of the 3D effects on BEM based simulations has been analyzed in comparison to high fidelity 3D CFD. The open source software Qblade has been used, creating models differing only in the input polars. CFD and BEM are two completely different methods, but in order to be compared and ensure the consistency of the comparison, a defined strategy needs to be developed. The origin of the polars is a crucial input that can make strong difference in the results. It has been shown that, even if consistent polars are used for BEM simulations and performances are in general improved, the rotational augmentation effects in the root area lead still to differences between BEM and CFD. At the same time, this is also the area in which the AoA is really sensitive to the chosen extrapolation method, see section 3.1.
The second important aspect is the number of the provided polars. BEM based models fully rely on the interpolation methods within the provided polars, but, when possible, a good strategy would be to increase the number of the input polars as much as possible. In the same way it is important to provide polars in ranges of AoA as big as possible, in order to avoid dependency of the results to BEM extrapolation methods.
The extraction of the effective AoA with BEM showed in some blade ranges differences with the one extracted in 3D CFD, as well as for the induction factor. This one has been shown to be a critical point for BEM calculations. For wind turbines operating in turbulent wake state, and therefore with high a, BEM calculations are not reliable. The differences in AoA and the high values of a are, in the authors' opinion, the main reason of discrepancies between 3D CFD and BEM calculations.
Further investigations can lead to the comparison of these models with field results, in order to validate also the 3D CFD results. Additionally, the polars that have been used, have been extracted from CFD, and it would be interesting to compare also polars from experiments, in order to give a complete overview of all the possibilities and results.