Progress in the validation of rotor aerodynamic codes using field data

Within the framework of the fourth phase of the International Energy Agency (IEA) Wind Task 29, a large comparison exercise between measurements and aeroelastic simulations has been carried out featuring three simulation cases in axial, sheared and yawed inflow conditions. Results were obtained from more than 19 simulation tools originating from 12 institutes, ranging in fidelity from blade element momentum (BEM) to computational fluid dynamics (CFDs) and compared to state-of-the-art field measurements from the 2 MW DanAero turbine. More than 15 different variable types ranging from lifting-line variables to blade surface pressures, loads and velocities have been compared for the different conditions, resulting in over 250 comparison plots. The result is a unique insight into the current status and accuracy of rotor aerodynamic modeling.

For axial flow conditions, a good agreement was found between the various code types, where a dedicated grid sensitivity study was necessary for the CFD simulations. However, compared to wind tunnel experiments on rotors featuring controlled conditions, it remains a challenge to achieve good agreement of absolute levels between simulations and measurements in the field. For sheared inflow conditions, uncertainties due to rotational and unsteady effects on airfoil data result in the CFD predictions standing out above the codes that need input of sectional airfoil data. However, it was demonstrated that using CFD-synthesized airfoil data is an effective means to bypass this shortcoming. For yawed flow conditions, it was observed that modeling of the skewed wake effect is still problematic for BEM codes where CFD and free vortex wake codes inherently model the underlying physics correctly. The next step is a comparison in turbulent inflow conditions, which is featured in IEA Wind Task 47.

Doing this analysis in cooperation under the auspices of the IEA Wind Technology Collaboration Program (TCP) has led to many mutual benefits for the participants. The large size of the consortium brought ample manpower for the analysis where the learning process by combining several complementary experiences and modeling techniques gave valuable insights that could not be found when the analysis is carried out individually.

1 Introduction

Wind turbine design codes are essential for the industry to assess lifetime and energy production before the investment is made to build a turbine prototype. The aerodynamic model is then one of the most challenging components of these codes, because every aerodynamic process, in its basis, is described by means of the so-called Navier–Stokes equations. These equations cannot be solved in an analytical way, where also a numerical solution of the Navier–Stokes equations is out of reach in design due to extreme computational demands. The difficulty of accurate aerodynamic modeling is perhaps most convincingly illustrated by the fact that solving the Navier–Stokes equations (as a matter of fact only proving that a smooth solution exists) is one of the seven Millennium Prize Problems as formulated by the Clay Mathematics Institute in 2000 (Fefferman, 2000). As such, every aerodynamic model inherently suffers from simplifications. For wind turbine aerodynamics, an additional difficulty arises from the fact that the computational effort for design calculations is more extreme than it is for most other applications (e.g., fixed wing aerospace); see Schepers (2012). This necessitates the use of engineering models that are not only very efficient but also very simplified aerodynamic models based on the blade element momentum (BEM) method. Obviously more advanced methods like computational fluid dynamic (CFD) codes are applied too, but their use is, due to the computational demands, restricted to specific studies and load cases. From a practical point of view, the simplifications in engineering methods inevitably go together with a large uncertainty band, which is even larger for modern MW-scale wind turbines. The larger uncertainties with increased rotor size are partly a result of unknown (high) Reynolds effects, where moreover the more flexible blades will lead to larger deflections and more pronounced non-linear aeroelastic behavior with unknown aerodynamic implications. Other uncertainties result from the thick(er) airfoils that are applied on large rotors, which are very difficult to model and measure accurately in an aerodynamic sense. Last but not least, the changed relation between the scales in atmospheric inflow to the scales of the turbine (blades) lead to larger uncertainties for increased rotor size.

In order to reduce the uncertainty band of aerodynamic models and to make them reliable enough for the design of cost-effective turbines, aerodynamic models need to be improved and validated with good measurements. Conventional wind turbine measurements of for example power and blade root bending moment lack sufficient detail for that purpose. More detailed sectional load information is necessary for a better validation and understanding. Historically, progress on this topic has taken advantage of international cooperation in research tasks under the auspices of the IEA Wind Technology Collaboration Program (TCP) (IEA, 2021), leading to many mutual benefits for the participants. IEA Wind Task 14 and 18 contributed to this objective where field measurements from all over the world (some including sectional pressure measurements) were studied by an international consortium of wind energy researchers (Schepers et al., 1997, 2002). However, the main conclusion was that constant, uniform and controlled inflow conditions are necessary to make progress in this field, which led to several wind tunnel experiments with rotating rotors. Amongst these was Phase VI of the Unsteady Aerodynamics Experiment (UAE), testing a two-bladed 10 m diameter wind turbine in the wind tunnel of NASA's Ames Research Center, featuring a test section of 80 ft (24 m) by 120 ft (37 m) (Hand et al., 2001). The measurements taken in this experiment have been the subject of investigation in IEA Wind Task 20 (Schreck, 2008), and a blind comparison to simulations has been carried out (Simms et al., 2001). A follow-up from these experiments was the European Union project MEXICO (Model rotor EXperiments In Controlled cOnditions) in which 10 institutes from six countries cooperated in doing experiments on an instrumented, three-bladed wind turbine of 4.5 m diameter placed in the open section of the large low-speed facility (LLF) of German-Dutch Wind Tunnels (DNW) in the Netherlands. These experiments, which also featured extensive flow field measurements using particle image velocimetry (PIV), featured campaigns in 2006 and 2014 and were the subject of analysis in IEA Wind Task 29 phases I to III (Schepers et al., 2012, 2014; Boorsma et al., 2018). A more detailed summary of dedicated wind tunnel experiments has recently been published online (Boorsma, 2021). Although these wind tunnels have been used successfully to validate rotor aerodynamic models, translating these results to “real life” flexible turbines in turbulent inflow conditions remains a challenge. Therefore the comparison rounds of IEA Task 29 Phase IV (Schepers et al., 2021) have focused on newly released field measurements on a 2 MW turbine from the DanAero experiment (Bak et al., 2010; Madsen et al., 2010b). The unique data from this experiment, including a description of the turbine, were made available to the participants of IEA Task 29 Phase IV so that it could form a basis for a thorough analysis. This paper presents the progress of this task, where many participants from different countries simulated the same experiment. The studies may serve as a benchmark for performing code-to-code comparison involving many participants. Suggestions will be given to improve the agreement between codes.

Section 2 presents the methodology of the comparison round, including a description of the measurements and the setup of the comparison. Sections 3 and 4 give the results of the comparison for the two cases under investigation together with a discussion, which is followed by conclusions.

2 Methodology

Firstly, a description is given of the DanAero experiment. Then the setup for the comparison rounds is given, including a summary of the simulation codes.

2.1 DanAero field measurement campaign

Detailed aerodynamic measurements on MW-scale wind turbines are scarce and open publications about them even more so. An exception lies in the DanAero experiment (Bak et al., 2010; Madsen et al., 2010b), which was carried out in atmospheric field conditions on an NM80 2.3 MW turbine in a Danish project by the Danish Technical University (DTU) and four industrial partners (LM Glasfiber, Siemens WindPower, Vestas and Dong Energy) in two periods from 2007 until 2010 and from 2010 until 2013. At the initiation of IEA Task 29 Phase IV in 2018, the DanAero partners agreed that the measurements, as well as the model data for aerodynamic and aeroelastic modeling of the NM80 turbine, could be shared with the partners participating in this task. The pitch-controlled turbine features three LM38.2 blades, resulting in an 80 m rotor diameter at a 57 m hub height. The level of detail from the instrumentation lies far above the level of conventional wind turbine measurements. In addition to conventional power and loads measurement using strain gauges, not only surface pressures at four sections along a blade (at 13, 19, 30 and 37 m from the rotor center) were measured but also inflow velocities using pitot tubes, and a row of surface flush-mounted microphones was installed at the outer part of the blade. Integrating the measured surface pressures around the sectional airfoil contour resulted in the chord-normal and tangential pressure forces at these four stations. Meteorological measurements were performed using a mast located 313 m (3.9 D – diameters) in a southwesterly direction (237^∘) upwind from the turbine. The meteorological mast included cup and sonic anemometers at seven heights up to 93 m. An overview is given in Fig. 1. More details about instrumentation can be found in the dedicated report (Bak et al., 2013).

Figure 1 DanAero turbine and instrumentation (Bak et al., 2010).

2.2 Setup

Several cases are under investigation focusing on axial and sheared–yawed inflow, which are described in more detail below. The comparison rounds mainly focus on the data obtained from the pressure measurements, i.e., pressure distributions and the derived normal and tangential force at the four instrumented stations. For the CFD modelers, the blade geometry was not only made available by means of an IGES (Initial Graphics Exchange Specification) file but also a pre-processed multi-block surface and volume mesh were distributed amongst the participants. For the lifting-line codes, airfoil data were prescribed. Hereto data were obtained from dedicated wind tunnel testing of the scanned sectional airfoil geometries of the real blade, which were 3D-corrected afterwards based on field measurements of generator power and blade bending moments (Bak and Fuglsang, 2004). A wide variety of over 19 different codes have been used by the participants ranging from BEM to CFD models for the rotor aerodynamics. For CFD models, both Reynolds-averaged Navier–Stokes (RANS), as well as detached eddy simulation (DES) formulations, have been used. Additionally, medium fidelity tools such as panel codes, actuator line (AL) and lifting-line free vortex wake (LL-FVW) models are part of the comparison. An attempt to summarize the different codes has been given in Table 1. For a more detailed description of the different simulation tools used by the participants of the comparison round, the reader is referred to the code description appendix of the final report of IEA Task 29 Phase IV (Schepers et al., 2021).

Collier (2019)Schwamborn et al. (2006)Michelsen (1992)Sørensen and Shen (2002)Meyer Forsting et al. (2019)Madsen et al. (2010a)Madsen et al. (2020)Pirrung et al. (2016, 2017)Le Cunff et al. (2013)Perdrizet et al. (2013)Bozonnet et al. (2017)Blondel et al. (2018)Greco and Testa (2021)Moriarty and Hansen (2005)Shaler et al. (2020)Sprague et al. (2019)Mudry (1982)Cambier (2008)Boorsma et al. (2011)van Garrel (2003)Bauchau et al. (2001)Schwamborn et al. (2006)Kroll et al. (2000)

Table 1 High-level summary of participant codes and settings.

The distribution of chord-normal and tangential forces along the blades is supplied by all participants, together with rotor axial force and torque. For the CFD and panel codes, also pressure distributions are compared to the measurements. The lifting-line codes that use airfoil data also supplied the distribution of so-called “lifting-line variables” (angle of attack (AoA), effective wind speed and induced velocities), which are compared between the simulations only. To improve the comparison of the results, also non-dimensionalized values of normal and tangential force are compared. These are non-dimensionalized using undisturbed local dynamic pressure (determined from wind and local rotational speed) to allow for a solid comparison of airfoil coefficients between experiment and simulations along the span without the complications of uncertainty in rotor-induced velocities. Hence the definition for the non-dimensionalized normal force, which is equivalent for the tangential force, can be given as

$$\begin{array}{}\text{(1)}& \mathrm{Fnc}\mathrm{\_}\mathrm{qc}={\displaystyle \frac{\mathrm{Fn}}{\mathrm{0.5}\mathit{\rho }\left(U_{\mathrm{\infty }}^{\mathrm{2}}+(\mathit{\omega }r{)}^{\mathrm{2}}\right)c}},\end{array}$$

where Fnc_qc [–] is the non-dimensionalized chord-normal force, Fn [N m⁻¹] is the chord-normal force, ω [rad s⁻¹] is the rotor speed, ρ [kg m⁻³] is the air density, U_∞ [m s⁻¹] is the wind speed, r [m] is the local radius and c [m] is the local chord.

The supplied axial force and torque have been post-processed to thrust and power coefficients Cd_ax and Cp using

$$\begin{array}{}\text{(2)}& {\mathrm{Cd}}_{\mathrm{ax}}={\displaystyle \frac{\mathrm{Fax}}{\mathrm{0.5}\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{2}}\mathit{\pi }{R}^{\mathrm{2}}}},\mathrm{and}\mathrm{Cp}={\displaystyle \frac{\mathrm{Torque}\mathit{\omega }}{\mathrm{0.5}\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{3}}\mathit{\pi }{R}^{\mathrm{2}}}},\end{array}$$

where Cd_ax [–] is the axial force coefficient, Cp [–] is the power coefficient, Fax [N] is the rotor axial force, Torque [N m⁻¹] is the rotor torque and R [m] is the rotor radius.

For the axial flow case, the aerodynamic flatwise moment is deduced using a script that linearly integrates the simulated force distribution along the blade span. For the yawed and sheared case, the flatwise moment was directly supplied by the participants. The displayed experimental values have been obtained from the post-processed strain gauge measurements by removing gravity and inertial contributions.

3 Case IV.1: axial uniform inflow

The first case, IV.1, is summarized in Table 2, based on a measurement data point in the summer of 2009 with relatively steady and uniform inflow and constant operational conditions. It can be observed that the turbine is highly loaded with an average value of 0.38 for the axial induction factor, which is generally considered to be the turbulent wake state. Unfortunately no measurement data were available with the turbine operating at a lower induction factor. For this case, the tilt angle (5^∘) and tower shadow effects are neglected by the modelers, but blade pre-bend is included. The case is subdivided into Case IV.1.1 where the turbine (i.e., tower, blades etc.) should be modeled rigid and Case IV.1.2 where flexibility is included. It is noted that although a comparison is made with the measurements (for which the blades are obviously flexible), the CFD simulations were mostly performed for a rigid blade.

Table 2 DanAero comparison cases (axial flow).

^* Estimate.

Figure 2 Lifting-line variables for Case IV.1.2.

3.1 Lifting-line codes

Selected lifting-line variable comparison plots for Case IV.1.2 are given in Fig. 2. These variables are calculated for all codes that need the input of sectional airfoil data.

The effective velocity U_eff (a composite of wind-, motion- and rotor-induced velocities) is in good agreement between the codes, indicating that the inputted operational conditions are consistent between the codes. Small differences can be observed in the inboard region, where induction starts to play a role over the elsewhere dominant rotational velocity. The angle of attack (AOA) shows larger variations, caused by the differences in axial- and tangential-induced velocities U_i and V_i, as shown in Fig. 2c and d. A closer look at the axial-induced velocity U_i shows that, especially apparent for the participants that delivered both BEM and free vortex wake results, the free vortex wake codes (usually depicted with a dashed line) feature a roughly 10 % lower induced velocity. A similar observation was made in the final report of Mexnext-III (Boorsma et al., 2018) for the New Mexico case in the turbulent wake state, which featured a high axial induction factor similar to the case under investigation here.

The resulting loads in terms of non-dimensionalized sectional normal and tangential forces Fn and Ft are given in Fig. 3. Consistent with the induced velocities and angles of attack, the agreement in Fn is fair between the computational results, with a spread of around 5 %. A closer look at the levels in Fig. 3a again shows a difference between vortex wake codes and BEM codes, where the latter generally feature a lower loading than vortex code results. This is also reflected in the integral performance from Fig. 3c and d (blade root flatwise moment, axial force and power). Acknowledging the earlier observed variations in axial-induced velocities, this discrepancy seems to be in contradiction with momentum theory, where a higher loading is accompanied with more induction. However this could also be attributed to uncertainties in the engineering extensions used for the turbulent wake state in which the turbine is operating for this case. Grouping the results by code type, as will be shown in Sect. 3.3, allows us to better observe the loading differences.

Figure 3 Loads, performance and deformations for Case IV.1.2. Note that for the bar plots the order of the bars agrees with the order of the legend entries.

The comparison with measured values (obtained from the integration of the measured pressure distributions) shows the normal force to be consistently overpredicted by approximately 10 % for all radial positions. It could then be expected that this overprediction is reflected in the measured flatwise moments from Fig. 3c as well. However, the observed differences for this moment are around ±2 %, where it should be noted that the measured flatwise moments are obtained from the strain gauge rather than the pressure measurements. Detailed investigations of possible causes for the low experimental normal force were carried out during IEA Task 29 Phase IV, unfortunately without success. Contrary to the normal force, the tangential force (Fig. 3b) reveals an under- instead of overprediction. Here it must be noted that the measured values only contain the contribution of pressure, whereas the simulations are based on airfoil coefficients, which also contain the friction forces. For the tangential direction, this contribution is significant (as will be shown in Sect. 3.2), and taking this into account will result in a better agreement between measurements and simulations. In addition to that, it is known that obtaining tangential force from the integration of a finite number of pressures is very sensitive to the distribution of the taps, especially around leading and trailing edges. In summary, tangential forces obtained from pressure distributions should be interpreted with care.

Figure 3e and f then displays the predicted deformations between the aeroelastic codes. Here it is shown that besides some outliers, most codes agree within 1 % for the flapwise deformation (d_x), whereas there is a bit more scatter in the edgewise deformation (d_y). However, in an absolute sense, we are comparing differences of less than a centimeter. As some of the participants did not provide results for the flexible case (IV.1.2), several colors are visible on the zero line. Although not depicted here, the differences due to the flexibility between Case IV.1.1 and IV.1.2 are mostly noticeable in the outboard region of the blades. Because of the flapwise pre-bend, which is “flattened” for the flexible case, the wind faces the blades more head-on, and the radius is slightly increased, resulting in slightly larger loads (and larger effective velocity, induced velocities and angle of attack) toward the tip.

3.2 CFD and panel codes

Starting with a comparison of the sectional geometries used by the codes, Fig. 4a and b shows the shapes of the most in- and outboard section. Theoretically these shapes should be identical, but some offsets appear and in some cases even local twist differences are apparent. The question remains whether these deviations are due to different ways of post-processing (e.g., a different orientation or location of the section) or if the geometries are not identical. As the blade geometry was prescribed by an IGES file, the latter should not be the case. It is acknowledged that retrieving a sectional blade slice in a 3D pre-bended rotor geometry can result in differences due to small inconsistencies like definition of radial coordinate, angular orientation of the sectional plane, and direction of the chord line. What adds to that is the fact that different software packages are being used together with the corresponding “human factor”. Hence it can be considered evident that differences appear, which is representative for what happens in real-life CFD applications.

Figure 4 Sectional geometry and pressure distributions (Case 1.1).

The pressure distribution comparison plots from Fig. 4 show the models tend to agree well for all radial positions, although the inviscid panel code is unable to predict separated flow in the most inboard station in Fig. 4c. But at this station, the CFD codes also struggle to match the measured pressure distribution on the suction side just before the trailing edge. For the other three radial positions, the suction levels are 5 % to 10 % higher in comparison to the experiment. Acknowledging the fact that most CFD codes model the boundary layer fully turbulent (whereas the experiment features free transition), this is contrary to expectation. At r=30.20 m (Fig. 4e), the linear shape of the measured pressure distribution on the pressure side just before the chordwise position of maximum thickness differs from the rounded shape, as predicted by all simulations. A separate sectional study was performed by DTU to compare the performance between design and measured geometries, of which the first was used for this comparison round. This however revealed very small differences in comparison to the differences shown here. Besides these small deviations, one can conclude that generally speaking the pressure distributions are in good agreement between simulations and between simulations and experiment.

Table 3 Coarse and fine mesh properties.

To further dive into the observed discrepancies between the CFD codes, a sensitivity study was carried out by DLR, DTU, UAS Kiel and the IAG University of Stuttgart. Hereto their CFD results were post-processed using a common tool at IAG to improve the consistency between the CFD simulations. Additional simulations were run using a finer mesh resolution to assess the impact of the grid resolution. The properties of the coarse and fine mesh are summarized in Table 3. Here it is noted that DTU used a fine mesh with full rotor topology, while the others used a periodic boundary condition.

Although for the periodic fine mesh the chordwise resolution is the same as the coarse mesh, and the first cell normal size is even larger, the total number of cells for a full rotor equivalent of this mesh would roughly differ by a factor of 6 ($\mathrm{26}\times \mathrm{3}/\mathrm{14}$). Hence, the majority of the extra cells for this fine mesh are placed in the wake and the background mesh.

In addition to the grid sensitivity study, the results of the coarse and fine mesh were post-processed to observe the contribution of friction to the integrated forces. The following results can then be visualized:

• Cf. It includes coarse mesh, friction and pressure forces.

• Cp. It only includes coarse mesh and pressure forces.

• Ff. It includes fine mesh, friction and pressure forces.

• Fp. It only includes fine mesh and pressure forces.

The resulting comparison plots are shown in Fig. 5. For the DTU result in Fig. 5a and b, it becomes clear that the mesh refinement hardly has an influence on the results, and the coarse mesh suffices. For the IAG University of Stuttgart results in Fig. 5c and d, this conclusion does not hold, and level differences up to 10 % can be observed. The same holds for the DLR results (which are not shown here for the sake of briefness). While the coarse mesh is already providing grid-independent results for the incompressible DTU EllipSys3D solver, the compressible solvers require further refinement. Although the underlying pressure distributions are not shown here, the compressible codes reveal a different pressure level on the suction side between the coarse and fine mesh results, which causes the differences in the shown sectional loads as obtained from the integrated pressure distributions. A similar finding was previously reported in the EU AVATAR project (Hansen et al., 2016). Figure 5e and f demonstrates that refining the mesh improves the consistency between the codes significantly. Apart from the blade root region inboard of a 10 m span, the results are close to identical.

Figure 5 Effect of CFD sensitivity studies on predicted loads for Case IV.1.1.

As expected, the effect of friction is hardly present in the normal force, which can be deduced from Fig. 5a and c. For the tangential force this is different, especially for the three most outboard stations that feature moderate angles of attack with attached flow. The friction reduces the tangential force by approximately 15 %, consistent between the two code results displayed here.

Figure 6 Loads comparison by model types for Case IV.1.1.

3.3 Model types comparison

In addition to displaying the loading results from the various codes, the supplied data also give the opportunity to calculate an average result for each code type. Here the following code types are distinguished:

• BEM, blade element momentum methods using the prescribed airfoil data set;

• FVW, lifting-line free vortex wake methods, also using the prescribed airfoil data; and

• CFD, computational fluid dynamics codes, which model the rotor blade geometry and the 3D space around it. In most cases, the blade boundary layer is modeled as fully turbulent to promote consistency between the results.

In addition to the code types listed above, there are also a panel code and actuator line codes that joined the comparison round, but they were excluded from the averaging as the number of codes for these types is deemed too few (<3) to obtain valuable statistics.

To obtain the loading averages, first the normal and tangential force are determined at the same spanwise four positions as the instrumented sections using linear interpolation from the supplied radial distributions. A simple average $\stackrel{\mathrm{‾}}{x}$ and standard error x_err between the supplied results of a code type are determined using

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{x}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left(x_i\right)\mathrm{and}{x}_{\mathrm{err}}=\sqrt{\sum _{i=\mathrm{1}}^n{\displaystyle \frac{{\left(x_i-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}}{n(n-\mathrm{1})}}},\end{array}$$

where x_i is the sample value, $\stackrel{\mathrm{‾}}{x}$ is the average, x_err is the standard error and n is the number of samples.

To give an indication of the variability between the results, a partly transparent band is plotted around the average of each code type illustrating the standard error x_err between the supplied results of a code type. The same procedure is applied to the flatwise moments, axial force and torque (average levels only). This process is only performed for the rigid results, because the largest number of CFD results is provided for this configuration.

Figure 7 Pressure distribution comparison at 0^∘ azimuth angle, Case IV.2.1. Animations of variation with an azimuth angle are available from the Supplement.

The results are illustrated in Fig. 6. One can conclude that the level of the normal force generally agrees within 5 % between the different code types. Acknowledging the fact that we are operating in the turbulent wake state, and BEM codes need an empirical relation to replace momentum theory, this difference between the code types is considered small. The same holds for the variability between the BEM codes, as different empirical relations exist that describe the turbulent wake state. Besides this observation, the tip region shows more of a falloff for CFD, and the root region features a sudden drop in level for the CFD results. Here it is noted that most BEM codes feature a Prandtl tip correction to correct for the non-uniformity between the blades (which is intrinsic to the FVW codes). It can be hypothesized that the larger load falloff in the tip region for CFD is related to the decambering effect on the tip sections of a finite blade, which is not properly accounted for by lifting-line codes. The higher loading for free vortex wake models as discussed in Sect. 3.1 is also apparent from the plots. For the tangential force, the differences between the code types are larger, but in an absolute sense these are rather small. The reason for the good trend overlay between CFD and FVW as opposed to BEM is not fully understood and could well be attributed to coincidence. The bands showing the variability between results from each code type are similar between the different types and amount to ±1 % for the normal force and ±5 % for the tangential force. In the previous phases of IEA Task 29, the human factor would result in a larger spread between CFD results. The dedicated efforts to minimize this spread as described in Sect. 3.2 seem to have paid off together with the fact that multiple iterations were performed to eliminate errors.

Comparing against the measurements, the normal force is overpredicted by approximately 10 % for almost all sections, while the flatwise moment agrees within ±2 % as was observed in the dedicated section discussing lifting-line code and CFD results. For the discussion on the tangential force discrepancy, the reader is referred to the lifting-line-code discussion in Sect. 3.1.

Table 4 DanAero comparison cases (shear and yawed flow).

^* Estimate of azimuth-averaged value.

Figure 8 Pressure distribution comparison at a 0^∘ azimuth angle, Case IV.2.2. Animations of variation with an azimuth angle are available from the Supplement.

Figure 9 Normal force (Fn) variation by model types, Case IV.2.1.

4 Case IV.2: sheared and yawed inflow

Two cases are defined featuring a significant vertical wind shear and yaw misalignment (also including vertical shear) as summarized in Table 4, carefully selected from the available measured time series. Again it can be observed that the turbine is highly loaded and operating in the turbulent wake state for both cases. For these cases, flexibility, tilt angle and tower shadow effects are included. Although the flexibility was accounted for in the BEM and FVW simulations, it was anticipated that for many CFD contributions this was not manageable (a rigid and tilted rotor without the tower effect will suffice for CFD). However, the main priority here is comparison to the measurements, and, fortunately, an investigation using a stiff and flexible aeroelastic model of the turbine in simulations indicated that the flexibility had a very small impact on the aerodynamic blade force variation. As we are studying the load variation due to shear and yaw, the results are presented as a function of the rotor azimuth angle. To better compare the load trends as a function of the azimuth angle, the mean over the rotor revolution has been subtracted from these results now showing a “delta”. For the measurements, because results over multiple revolutions are available, the standard deviation between these measured values gives an indication of the repeatability. The standard deviation is indicated in the graphs by a grey band around the mean value.

Figure 10 Normal force (Fn) variation by model types, Case IV.2.2.

4.1 Pressure distributions

Pressure distributions at a 0^∘ azimuth are given for all four radial stations in Figs. 7 and 8 for Case IV.2.1 and IV.2.2 respectively. Also note that animations featuring four different azimuth angles (0, 90, 180 and 270^∘) are available from the Supplement. Generally speaking the trend with azimuth is well captured, but the separated flow conditions for Case IV.2.1 (especially inboard) are a challenge for the panel code and some of the CFD codes. For both cases, the measured pressure distributions at 92 % R (Figs. 7d and 8d) feature a dip at the suction side around 20 % chord, which seem rather fierce to have be introduced by transition. It is noted that for these dynamic cases involving shear, most of the CFD modelers employed a different meshing strategy in comparison to the uni-axial case. More details can be found in the detailed code descriptions from the final report (Schepers et al., 2021).

4.2 Model types

Similar to the axial flow analysis described in Sect. 3.3, results between code types have been averaged to give a better overview of the differences between them. To obtain the loading averages and standard error, the same data reduction procedure is adopted as for axial flow. Because we are interested in the load variation as a function of azimuth, all supplied code results are linearly mapped onto an azimuth angle distribution with a 5^∘ step, prior to calculating the average and standard error for each code type.

The results in terms of normal force variation are illustrated in Figs. 9 and 10. Although the average level of the result is removed, since we are focusing on the azimuthal variation, these were mostly in line with the differences observed in Case IV.1 in axial inflow. The trends for the sheared case in Fig. 9 clearly illustrate the benefit of CFD simulations over the lifting-line codes. The variability within each code type (as illustrated by the colored band around the mean results) is of the same order for all code types for this case, although there seems to be a significant variability in the amplitude of the dip of around 180^∘ azimuth as predicted by CFD codes.

Figure 11 Normal force (Fn) variation and airfoil data using CFD-synthesized airfoil data, Case IV.2.1. USTUTT reference CFD data are plotted as well.

The yawed case shows different trends, as it is dominated by induction rather than airfoil aerodynamics. This is illustrated by the normal force in Fig. 10, with a maximum near 270^∘ azimuth for the outer part of the blade and near 90^∘ azimuth for the inner part. This is explained by induction effects from the skewed wake; see Schepers (2012): at the inner part of the blade the azimuthal load variation at yawed conditions is mainly determined by the induction from the skewed root vortex, whereas the load variation at the outer part of the blade is mainly determined by the induction from the skewed tip vortex, which gives a 180^∘ different phase shift to the induction (and resulting load) variation. This skewed wake effect is inherently modeled by CFD and vortex methods by which they almost completely fall within the experimental uncertainty band, especially for the outboard sections. BEM codes rely on uncertain engineering methods to model these skewed wake effects by which the variability between BEM codes is more than 2 times larger than vortex and CFD codes.

4.3 CFD-synthesized airfoil data

The poor agreement of lifting-line codes with measurements for the shear case together with the encouraging CFD results instigated further investigations. The IAG University of Stuttgart provided an airfoil data set synthesized from 3D rotational CFD computations employing the azimuthal averaging method as described in Hansen et al. (1997) and Bangga (2018). Two participants then re-simulated Case IV.2.1 with their BEM codes. The resulting normal force trends are given in Fig. 11.

It is clear that for Case IV.2.1 the results featuring the CFD-synthesized polars (dashed lines) outperform the results with the original polars. This indicates the earlier found discrepancies to be airfoil data related. It is noted that the airfoils are operating in the stall region, and apparently the stall delay due to the 3D rotational and unsteady effects is not correctly taken into account by the original airfoil data (despite the fact that a 3D correction was applied (Bak and Fuglsang, 2004)). Additionally it was found that application of a classical 3D-correction model like Snel Snel et al. (1993) also can not predict the stall delay correctly. Comparing the underlying polars at r=48 % R in Fig. 11f clearly shows the stall delay as predicted by CFD. This explains, in the case of a downward-pointing blade when using the original airfoil data, the loading and induction augmentation as the lift increases for a decreasing angle of attack (when coming back from stall around 12^∘). It is noted that, although this reversed trend is most apparent for the 48 % R station, the CFD-synthesized airfoil data clearly improve the agreement with measurements also for the outboard stations. This observation hints at the need for a stall delay model that works along the whole blade span instead of only the inboard part. For Case IV.2.2, which is dominated by induction aerodynamics and featuring angles of attack further away from stall, it can be shown that the CFD-synthesized polars hardly impact the trend variation.

5 Conclusions

A large comparison exercise has been performed featuring three simulation cases in axial, sheared and yawed inflow conditions. Results were obtained from more than 19 simulation tools originating from 12 institutes ranging in fidelity from BEM to CFD and compared to state-of-the-art field measurements from the DanAero turbine. More than 15 different variable types ranging from lifting-line variables to pressures, loads and velocities have been compared for the different conditions, resulting in over 250 comparison plots. The result is a unique insight into the current status and accuracy of rotor aerodynamic modeling.

For axial flow conditions, a good agreement was found between the various code types, where a dedicated grid sensitivity study was necessary for the CFD simulations. However, compared to wind tunnel experiments that feature controlled conditions (like New Mexico), it remains a challenge to achieve good agreement of absolute levels between simulations and measurements in the field. Considerable efforts were spent on investigating possible causes in the measurements for the deviations to the simulations without success. For sheared inflow conditions, uncertainties due to rotational effects on airfoil data result in the CFD results to stand out above the codes that need input of sectional airfoil data. However, it was demonstrated that using CFD-synthesized airfoil data is an effective means to bypass this shortcoming. For yawed flow conditions, it was observed that modeling of the skewed wake effect is still problematic for BEM codes, whereas CFD and FVW codes inherently model the underlying physics correctly. The next step is a comparison in turbulent inflow conditions, which is featured in IEA Wind Task 47.

Doing this analysis in cooperation with and under the auspices of IEA TCP Wind has led to many mutual benefits for the participants. The large size of the consortium brought ample manpower for the analysis where the learning process by combining several complementary experiences and modeling techniques gave valuable insights that could not be found when the analysis was carried out individually.