the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Aeroelastic analysis of wind turbines under turbulent inflow conditions
Giorgia Guma
Galih Bangga
Thorsten Lutz
Ewald Krämer
The aeroelastic response of a 2 MW NM80 turbine with a rotor diameter of 80 m and interaction phenomena are investigated by the use of a highfidelity model. A timeaccurate unsteady fluid–structure interaction (FSI) coupling is used between a computational fluid dynamics (CFD) code for the aerodynamic response and a multibody simulation (MBS) code for the structural response. Different CFD models of the same turbine with increasing complexity and technical details are coupled to the same MBS model in order to identify the impact of the different modeling approaches. The influence of the blade and tower flexibility and of the inflow turbulence is analyzed starting from a specific case of the DANAERO experiment, where a comparison with experimental data is given. A wider range of uniform inflow velocities are investigated by the use of a blade element momentum (BEM) aerodynamic model. Lastly a fatigue analysis is performed from load signals in order to identify the most damaging load cycles and the fatigue ratio between the different models, showing that a highly turbulent inflow has a larger impact than flexibility, when low inflow velocities are considered. The results without the injection of turbulence are also discussed and compared to the ones provided by the BEM code AeroDyn.
 Article
(4877 KB)  Companion paper
 BibTeX
 EndNote
The current design trend of wind turbines is leading to rotor diameters becoming larger and larger, but they have to be light in order to decrease the cost of wind power generation in terms of leveling energy costs (USD per kWh) and make wind power generation a competitive resource in comparison to other electric generation systems. A lot of research is being carried out to investigate materials and construction techniques in order to allow lighter designs with the consequence that the rotor blades are becoming more and more flexible, which leads to large deformations with associated nonstationary loads and oscillations, resulting in unexpected changes in performances or even flutter if the damping is negative. Additionally, large rotor wind turbines are in reality subjected to diverse inflow conditions, such as shear, turbulence and complex terrain, leading to higher load fluctuations. Moreover, the aeroelastic instabilities strongly affect the operational life of wind turbines (M. O. L. Hansen et al., 2006). Most of the available simulation tools for wind turbines aeroelasticity are based on engineering models like blade element momentum (BEM) for the aerodynamics and 1D multibody simulation (MBS) for the structural response, like for example in Riziotis et al. (2008) and Jeong et al. (2011). These models are cheap but rely on different correction models to take unsteadiness and 3D effects into account (Madsen et al., 2012). In recent years, highfidelity fluid–structure interaction (FSI) has been frequently used for wind turbine applications. Sayed et al. (2016) implemented a coupling of the CFD (computational fluid dynamics) solver FLOWer to the CSD (computational structure dynamics) solver Carat$++$, where only the blades have been coupled to either a 1D beam or a 2D shell structural model. Yu and Kwon (2014) used a loose CFD–CSD coupling with an incompressible CFD solver and nonlinear Euler–Bernoulli beam elements for the structure in order to investigate the aeroelastic response of the generic NREL 5 MW rotor. The communication in this case was only once per revolution. The same turbine was also used by Bazilevs et al. (2011) and Hsu and Bazilevs (2012) by means of FSI between a loworder arbitrary Lagrangian–Eulerian variational multiscale (ALEVMS) flow solver and a structural solved based on a nonuniform rational basis spline (NURBS). For the same turbine, Heinz et al. (2016) compared the coupling of the flow solver EllypSys3D with the aeroelastic solver HAWC2 to the BEM results of HAWC2 alone. While they considered uniform inflow, Li et al. (2017) additionally considered a turbulent inflow synthetically generated by the use of a Mann box (Mann, 1994). Dose et al. (2018) presented a method to couple the flow solver OpenFOAM to the FEMbased beam solver BeamFOAM. A CFD–MBS coupling between the unsteady Reynoldsaveraged Navier–Stokes (URANS) solver TURNS and the MBS solver MBDyn was used by Masarati and Sitaraman (2011) to investigate the NREL Phase VI rotor.
Wind turbines are especially susceptible to fatigue damage, due to the oscillating characteristic of the affecting loads. Fatigue analyses are normally performed by manufacturers for certification purposes, and therefore such analyses are mostly BEMbased. In the EU project AVATAR (Schepers, 2016) it was shown that BEMbased calculations against highfidelity calculations led to a 15 % error in the computation of fatigue. This error motivated the TKI WoZ VortexLoads project (Boorsma et al., 2019), where starting from turbulent inflow conditions BEMbased and CFDbased calculations have been compared with each other and to experimental results.
Within the scope of the present study, a highly accurate CFDbased aeroelastic model of a 2 MW wind turbine was created and applied to study unsteady load characteristics. The objective was to identify the impact of the modeling of the individual turbine components and the occurring interactions on the transient loads. To achieve this goal, numerical models of successively increasing complexity are introduced. Starting from a onethird model of the blade in uniform inflow, over a complete rotor up to a complete flexible turbine in turbulent inflow, the transient loads were analyzed and compared. The aim was to analyze the main drivers for the load fluctuations and the damage equivalent loading (DEL) using highly accurate models. The different CFD configurations have been analyzed in detail because their computational costs vary enormously. It is therefore of interest, especially for the industry, to know limitations and differences within the highfidelity approaches. For the uniform inflow case, a comparison with BEMbased calculations is given and two additional inflow conditions are computed, because of its cheapness, in order to determine the generalization level of the results. The ability of BEM in predicting reliable fatigue values changing the computational settings is discussed. In Sect. 2 of this paper, the highfidelity framework (as presented in Klein et al., 2018) is described for fluid–structure interaction coupled simulations on the NM80 2 MW wind turbine rotor, also known as the DANAERO rotor (DANAERO project, 2020). The inflow conditions and setup for the different cases are described. Furthermore, the BEM model of the turbine is described with its validation, based on the usage of 3D CFD polars in order to ensure consistency with the highfidelity model. In Sect. 3, the aeroelastic response of the reference turbine is shown and the difference between the modeling approaches is exposed. Lastly, DEL calculation is performed in postprocessing of the different simulations, using two different timevarying input variables.
2.1 DANAERO wind turbine
The DANAERO wind turbine rotor is used for this paper. This is the reference wind turbine in the IEA Task 29 IV, also known as Mexnext IV (IEA Task 29, 2020). In this project different institutions and universities around the world compare their own codes and approaches, using them for the calculations organised into different subtasks of the same project. The results are compared not only to each other but also to experimental results provided by the DANAERO experiment (Aagaard Madsen et al., 2010). The experiments were conducted between 2007 and 2010 in cooperation between the Technical University of Denmark and the industrial partners Vestas, Siemens LM and DONG Energy, and then they were postprocessed and calibrated in the followup project DANAEROII (Troldborg et al., 2013). In this way it is possible not only to understand limitations and problematics of the different approaches but also to improve them. The turbine has a rotor diameter of around 80 m, a tilt angle of 5^{∘} and an around 1.4 m prebend. The hub, nacelle and tower have been modeled within the present study as cylinders, based on the available diameter distribution provided in the structural model.
2.2 CFD model and inflow conditions
The simulations are performed with the CFD code FLOWer (Raddatz, 2009). First developed at the German Aerospace Center (DLR), FLOWer has been expanded for many years now at the Institute of Aerodynamics and Gas Dynamics (IAG) for helicopter and wind turbine applications. It is a URANS and detachededdysimulation (DES) finitevolume solver for structured meshes. The present simulations are run using the shearstresstransport (SST) k–ω model according to Menter (1994), using a fully turbulent boundary layer. Two different spatial discretization schemes are available, a secondorder central cellcentered Jameson–Schmidt–Turkel (JST) (Jameson et al., 1981) and a fifthorder weighted essentially nonoscillatory (WENO) (Kowarsch et al., 2013) scheme. The second one is applied in the present study on the background mesh in order to reduce the dissipation of the vortices. The timestepping scheme is an artificial fivestage Runge–Kutta scheme, and multigrid level 3 is applied to accelerate the convergence of the solution. The time integration scheme is an implicit procedure called dual time stepping where at the beginning of each time step t an estimation of the solution is guessed. The closer this is to the final value, the smaller the necessary number of inner iterations to reach convergence. Independent grids need to be created for each single component, combined and overlapped by the use of the chimera technique.
The CFD model of the blade is created from the provided CAD file, where a “watertight” outer surface is extracted. For the hub, nacelle and tower, surface databases are recreated (cylinderbased) from provided geometrical properties. Meshes are generated by the use of the commercial software Pointwise in combination with inhouse scripts. All components have been meshed ensuring ${y}^{+}\le \mathrm{1}$ in the boundary layer region. The blades are meshed in an Omesh topology with 257 points over the profile and 201 points in the radial direction, for a total of around 9 million cells for each blade. The background mesh consists of hanging grid nodes in which the component meshes are embedded with the chimera technique. Three different CFD models have been created for the turbine, with increasing fidelity:

onethird model (BMU) of the rotor (only one blade) suited for uniform inflow conditions,

full model of the turbine (FMU) including nacelle and tower suited for uniform inflow conditions,

full model of the turbine (FMT) including nacelle and tower suited for turbulent inflow conditions.
The differences between the three models consist in the backgrounds that were used. Model 1 has no ground, because it is just a 120^{∘} model of the turbine. Model 2 has no friction on the ground in order to avoid the generation of a wind profile. Finally, model 3 has friction on the ground in order to consequently propagate the sheared turbulent inflow and is much more expensive in comparison to model 2 (87 million cells against 58 million), because an additional refinement is added upwind where the turbulence is injected, and different boundary conditions need to be applied in order to ensure a correct propagation of the turbulence. The 120^{∘} model is much cheaper than the other two, because it uses the periodic characteristic of a threebladed wind turbine, but of course it considers neither tilt angle nor tower influence. The different boundary conditions and CFD models are depicted in Fig. 1. In the following the meaning of the different boundary conditions is clarified:

NAVIER–STOKES and EULER wall represent the ground with and without friction, respectively;

FARFIELD represents the uniform inflow boundary condition;

PERIODIC and PERIODIC ROT represent the symmetrical boundary condition for the full and 120^{∘} model, respectively;

GUST is the Dirichlet boundary condition, by which arbitrary unsteady inflow can be applied;

PRESSURE OUTLET defines the outflow based on pressure.
All simulations are run based on the conditions defined in the subtask 3.1 of the IEA Task 29; see IEA Task 29 (2020). Those require a rated inflow velocity of 6.1 m s^{−1} in the uniform case. For FMT, synthetic turbulence is generated by the use of a Mann box (Mann, 1994) and injected in the flow field at a plane four diameters (4D) upstream of the tower bottom. This is added using a momentum source term as prescribed in Troldborg et al. (2014) and superimposed onto the steady uniform inflow. The turbulence on this plane is updated every time step using Taylor's frozen turbulence hypothesis (Troldborg et al., 2014). A turbulence intensity (TI) of 20 %, a length scale of 0.59^{*} hub height (according to the IEC standard normative 61400) and a stretching factor Γ=3.9 to approximate the Kaimal spectral model (as prescribed in Kim et al., 2018) are preset. A mesh refinement of the background is applied from the inflow plane in order to allow a better propagation of the turbulence. The effective TI at the rotor is usually lower than the one prescribed in the Mann box, because it decays for both physical and numerical reasons. From an empty box calculation with a TI of 6.8 %, a turbulence decay of around 14 % was calculated, and therefore it is assumed for this case that the effective TI amounts to 17.2 %. Sheared inflow is superimposed by the use of a power law with α=0.025. Due to the low reference velocity considered during the DANAERO experiment, a very high TI was chosen in order to be able to identify distinctively the effects of a turbulent atmospheric boundary layer. Delayed detachededdy simulation (DDES) is used instead of URANS for the CFD solution, changing the boundary conditions accordingly.
2.3 MBS solver
2.3.1 Structural model
The multibody dynamics (MBD) simulation code Simpack is used to simulate the structural dynamics of the turbine (as in Jassmann et al., 2014, and Luhmann et al., 2017). The structural properties of the entire turbine have been modeled starting from the provided HAWC2 aeroelastic data. A multibody system consists of rigid or flexible bodies interconnected by force and joint elements that impose kinematic and dynamic constraints. Each body, represented by one or more markers, may then have three translational and rotational displacements as a result of deformations and motion. The body motion is described by a set of differential–algebraic equations (DAEs), a combination of differential motion equations and algebraic constraints. The blades are modeled as nonlinear SIMBEAM body types (threedimensional beam structures in Simpack, described by a nodebased nonlinear finitedifferences approach). These have been discretized into 22 Timoshenko elements in the radial direction, also taking into consideration gravitational and centrifugal forces. Structural damping is applied using the Rayleigh damping model with α=0.025 and β=0.014. Due to its small expected deflections, the tower has been modeled as a linear SIMBEAM discretized into 25 Euler–Bernoulli elements, the hub has been modeled with 2 linear Euler–Bernoulli elements and the nacelle is modeled with only one rigid node; i.e., it can move but not deform. Loads provided from the CFD are damped for the first 200 time steps (equivalent to 200 azimuth degrees) in order to avoid strong and fast deformations that can lead to numerical instabilities in the calculation. In order to validate the structural model, the natural frequencies of the single blade and turbine are compared to the measured ones from M. H. Hansen et al. (2006) in Table 1.
2.4 BEM model
A simplified aerodynamic model based on Blade Element Momentum (BEM) theory has been generated with the NREL code AeroDyn (AeroDyn Theory Manual, 2005). This has the advantage of being already incorporated in Simpack as additional module, and it can be therefore easily coupled to the structural model. In this case, the blade needs to be modeled aerodynamically with as many nodes as structurally, i.e., 21 for each blade. Polars have been extracted from 3D CFD calculations in order to avoid the use of any tip or hub correction model and ensure as much consistency as possible to the CFD calculations, as it was already shown in Guma et al. (2018). The 3D polars have been provided in a range of angles of attack (AOAs) of between around −30 and +30^{∘} and have been extracted from the CFD solution using the RAV method (Rahimi et al., 2018) and then extrapolated up to −180 to +180^{∘} using the Viterna method. Axial and tangential induction corrections have been taken into account. Tower shadow effect has been taken into account depending on the computed case (single blade or full turbine). The comparison of the sectional loads per unit length in the normal (F_{x}) and tangential (F_{y}) direction between BEM and CFD is depicted in Fig. 2. In this case only one blade, with no tower shadow and rigid conditions, has been taken into consideration, averaging the results of the last three revolutions. The curves show a good agreement, and therefore the BEM model of the turbine is validated. A discussion of limitations and capabilities of BEM under turbulent inflow conditions is out of the scope of this paper. This aspect has already been addressed by Madsen et al. (2018), who compared BEMbased simulations using the aeroelastic tool HAWC2 to the highfidelity code EllipSys3D and to experiments. A good agreement was found between the three, although CFD predicted an unforeseen stall in the inboard regions. In the present work only uniform inflow cases have been calculated using BEM as the aerodynamic model of the turbine. The chosen setups are shown in Table 2.
2.5 FSI setup and computed cases
In order to allow the communication between FLOWer and Simpack, moving, undeformed and reference system markers need to be defined as prescribed in Klein et al. (2018). In the present study no controller is taken into account, which is why each simulation is conducted with a fixed rotational speed and pitch. These have been set according to the inflow velocity of 6.1 m s^{−1}, which is at the same time as the chosen uniform inflow velocity and the average velocity at which the Mann box is generated. Even if a high TI is set, the resulting velocity is always far away from cutoff. Therefore, the controller would mainly change the rounds per minute (RPM) and not the pitch angle. The change in RPM has an influence on the full system natural frequencies (that is expected to be small), on the blade–tower passage frequency and on the thrust. This would increase with the RPM and therefore the flapwise tip deformations. The coupling algorithm used is explicit; i.e., deformations and loads are exchanged only once per physical time step. In particular, the loads at the end of the flow calculation time step are used to calculate deformations that are applied to the subsequent step; see Fig. 3. The chosen time step in this case corresponds to 1 azimuthal degree. An already converged rigid simulation of the turbine that has already run for at least 10 revolutions is used as the restart for the coupled simulation in order to speed up the calculation and save computational time. The DANAERO rotor has high induction; therefore it takes many revolutions for the wake to fully develop and for the loads to stabilize. In order to save computational time, turbulence is injected and flexibility is activated only after a cheaper simulation (FMU) reached a low residuum, a difference lower than 1 % in the averaged loads and deformations between two revolutions and a wake development long enough to avoid effects on the loads too.
For the BMU case it was sufficient to run the coupled simulation for only 6 further revolutions to achieve convergence and periodicity of the results. For the FMU, RMU and FMT at least 10 revolutions have been run, although periodicity cannot be reached in the FMT case, because the simulation time is much shorter than the length of the Mann box used. The elapsed time for the coupled simulations (starting from a rigid converged solution) varies from a minimum of 15 h with 1632 processors for the BMU to a maximum of 48 h with 4320 cores for the FMT case. All simulations are run on the SuperMUCNG supercomputer at the LeibnizRechenzentrum in Munich.
All the CFD–MBD computed cases and differences can be seen in Table 3. For each mentioned case a rigid and a coupled version is available, although RMU R (rigid) and FMU R (rigid) represent the same case.
2.6 Damage equivalent loading (DEL)
The DEL is a constant load that leads, when applied for a defined number of cycles, to the same damage as that caused by a timevarying load over the same period. With this method, two or more signals can be compared in order to obtain insight into the fatigue loadings that blades are facing during normal operation. The approach is based on the S–N curves (stress vs number of cycles) of the material on a log–log scale so that the material behavior is defined by the slope of a line. Additionally, a rainflow algorithm is applied to recognize the relative fatigue cycles in a load signal by filtering peaks and valleys. This algorithm allows us to estimate the amount of load change depending on the amplitude of the cycle. In this way closed stress hysteresis cycles can be identified, defining not only their amplitude but also how often they appear. The consequent damage is, in fact, dependent on the combination of the last two factors. The formulation used in this paper is the one from Hendrinks and Bulder (1995) in which the different load signals are compared on a quantitative basis, using not only the range but also the mean of the load cycles. According to this method, the final expression of the DEL resulting from a provided signal is
where n is the total number of cycles detected by the rainflow counting; S_{r,i} is the amplitude of the ith cycle; S_{u} is the ultimate load; S_{m,i} is the mean value of the ith cycle; Neq is the number of cycles corresponding to DEL; S_{m,eq} is the equivalent mean value of the cycle with amplitude DEL; and, finally, m is the slope of the S–N curve, considering a symmetric Goodman diagram with straight life lines.
S_{r,i} and S_{m,i} are direct output of the rainflow counting, meaning that they are an individual and inevitable characteristic of the spectrum itself. Differently, Neq, S_{u}, S_{m,eq} and m need to be chosen in advance. S_{u} and m are material dependent, where a log–log S–N curve is considered in order to have a straight line and a constant m, while S_{u} can be calculated in first approximation as 5 times the maximum load in the provided spectrum. Neq and S_{m,eq} are user dependent. It is then clear that the absolute value computed by the DEL strongly depends on the choice of the constants, but as long as the same constants are considered, the DEL values are consistent within themselves and, therefore, comparable.
3.1 Aeroelastic effects
In this first section, the effects of aeroelasticity on the reference wind turbine are analyzed. The considered DANAERO experiment was performed at a low inflow velocity (6.1 m s^{−1}); that is why it is expected to have small deformations and therefore especially a low tower effect. The structural model used is always the same, imposing opportunely the flexibility of the components as prescribed in Table 3. This means that the calculation of gravitational and centrifugal forces, which is made directly in Simpack, is always taking the tilt angle into account, even in the BMU case.
As validation of the results, the sectional normal (F_{N}) and tangential (F_{T}) loads according to the chord length for three different radial positions in comparison to experiments are shown in Fig. 4.
Results of different field tests have been considered and averaged (black line). As described in Sect. 2.2, turbulence has been generated in a stochastic way, and therefore the experimental and simulation time series of each revolution are not directly comparable but need to be averaged. For the validation, two different test cases have been compared: an entire CFD model of the turbine with flexible blades with uniform inflow conditions (RMU C or RMU Flex, blue line) and an onlyrotor CFD model that is completely rigid but with an inflow turbulence comparable to the experiments (CFD Turb, red line). It can be seen that in the outside region, although a correct modeling of the inflow provides results closer to the experiment, the shape of the experimental curve is mostly well matched by the RMU C curve. In the hub region, the two modeling approaches do not show much difference from each other, although the flexible case gives slightly better results.
3.1.1 BMU vs. RMU
The first considerations are made comparing BMU and RMU; the two differ from each other by the presence of a rigid tower and a tilt angle in the CFD model. Deformations in the flapwise, edgewise and torsion direction of the tip of the blade can be seen in Fig. 5. It can be noticed that, due to the inertia of the blade, the tip deformation starts its downturn by 180^{∘} but shows this local minimum with a delay of around 20^{∘} by 2.35 % of the rotor radius.
A clear sinusoidal trend can be seen in both cases, which leads to an oscillation of the tip deflection from around 2.3 % to 2.5 % of the blade radius for the BMU case and from around 2.2 % to 2.5 % for the RMU case. The reason for this is the presence of the tilt angle (5^{∘}) that leads the gravitational and centrifugal forces to produce an oscillating deformation component in the flapwise direction. On the contrary, the aerodynamic contribution remains almost constant in time, with an oscillation smaller than 1 %. As previously mentioned, the CFD model in BMU has no tilt, but the structural model does, which is why the resulting centrifugal and gravitational forces are accordingly affected. This leads to the oscillation in the response of BMU. This oscillation turns out to be stronger than the blade–tower passage for RMU; therefore after the minimum due to the blade–tower interaction, there is a recovery that immediately collapses in order to follow the sinusoidal trend. The difference in the maximum deflection between BMU and RMU is 2.4 % and is due to a higher oscillation of the affecting loads in the rigid version of RMU, as can be seen in Fig. 6, where the global thrust (F_{x}) and torque (M_{x}) in the rigid and coupled case on the blade are plotted.
The tip deformations in the edgewise direction are only dependent on the gravitational forces and show therefore almost no difference between BMU and RMU. The same happens for the torsion, whose minimum value is slightly lower in RMU with a very low maximum value of 0.075^{∘}.
Regarding the global thrust and torque in the BMU case for rigid and coupled conditions, it can be seen that M_{x} in Fig. 6 has an oscillatory trend, directly related to the sinusoidal oscillation of the blade. The global thrust is slightly shifted to higher values in the case of coupling, where the mean value increases 1 %. This is due to the deformation of a prebent blade, resulting in an increase in the effective rotor surface. Even if the torque oscillates more in the flexible case than in the rigid state, the average difference is lower than 0.1 % and therefore negligible. The RMU case shows a larger oscillation due to the tower passage, and as in the BMU case, the structural coupling leads to a shift of both thrust and torque curves to higher values. In particular, directly before the tower passage, the flexible blade reaches higher values of thrust (on average 1 % to 2 % more) with a consequently higher thrust in front of the tower (on average 2 % to 3 % more). The same effect, although less evident, can be seen for the torque. Averaging over three revolutions, the maximum difference in the produced power is up 2.3 % and can be seen between BMU R (rigid) and RMU C (coupled). Lastly, the difference in the sectional loads, averaged over the last revolution, is analyzed in Fig. 7. These are the sectional forces normal and tangential to the rotor plane.
The normal forces in coupled and rigid conditions show almost no difference. In the tangential loads, the ones responsible for the power at the shaft, a small increase (around 1 %) can be observed at between 40 % and 60 % of the blade radius, due to a local slightly higher angle of attack (around 0.8 % more), connected with the positive value of torsion shown before and due to the increase in the effective rotor area.
While the CFD calculations have been made based on the operating conditions of the DANAERO experiment, further simulations have been conducted using BEM in order to determine the generalization level of the results. Tip deformations in the flapwise direction can be seen in Fig. 8a–c. where an oscillation from 2.3 % to 2.5 % of the blade radius can be observed as in CFD. In these BEM calculations the tilt angle needs to be either in both aerodynamic and structural models or in neither of them; therefore the only difference between BMU and RMU is the blade–tower passage effect. Differently from CFD, where the impact was almost negligible, large oscillations occur due to the blade–tower passage, which already for the case with an inflow velocity of 6.1 m s^{−1} decreases up to 10 % (in comparison to no tower shadow). Increasing the inflow velocity and the RPM, these oscillations become strong enough to preclude the deformations to reobtain the same shape as in BMU. An overestimated blade–tower passage effect can be observed in the produced torque too; see Fig. 8d–f. In particular, while with CFD a reduction in this effect was observed when the structures were flexible (by low inflow velocity), this does not appear using BEM, which shows only an increase in it for high velocities of around 11 % (see Fig. 8f). At the same time, while flexibility shows almost no effect on the average torque at low velocities, an up to 6 % difference can be observed at 13 m s^{−1}. Especially in this case it can be seen that the RMU C case converges back to the sinusoidal form of BMU C after a time equivalent of 150^{∘} at which this oscillation is damped out.
3.1.2 RMU vs. FMU
As mentioned in Sect. 2.5, the difference between RMU and FMU consists of the flexibility of the tower and nacelle. The flapwise, edgewise and torsion deformations comparing between RMU and FMU can be seen in Fig. 9. Due to the low inflow velocity, the tower deflection contributes only 0.1 % of the blade radius to the total blade outofplane deflection.
Considering the edgewise deflection, the average value increases from 0.43 % of the blade length for RMU to 0.65 % for FMU due to the additional contribution of the tower top deformation. For the same aforementioned reasons, the torsion deflection has on average the same value, but due to the tower's torsion contribution, it shows a higher amplitude of the oscillation that increases in the FMU case by up to 17 % more. The global thrust (F_{x}) and torque (M_{x}) can be seen for the RMU and FMU rigid and coupled conditions in Fig. 10, where almost no difference is shown between RMU and FMU coupled, due to the small deflections of the tower top.
As for the difference in FMU between rigid and coupled conditions, it can be seen that the decay due to the tower passage decreases by 6 % (difference in M_{x} between rigid and coupled at 180^{∘}). This has a direct effect on the maximum value reached directly after the recovery, which is also always higher than in the rigid case. It can also be observed that within one revolution the amplitude of the oscillation is higher in the coupled simulation. By averaging the results over the revolutions, it is found that the coupled case produces 3.5 % more power than the rigid case. In order to understand this behavior, the averaged sectional loads of the FMU rigid and coupled cases are compared; see Fig. 11. The area of interest is from 20 % of the blade radius, because near the hub the difference between the two curves is mostly due to the strong unsteadiness affecting the hub region, where separation is occurring. The loads in the normal direction F_{x} are not affected at all by the coupling. In contrast, the tangential loads F_{y}, the ones generating the torque M_{x} and therefore the power, show some difference in the range of between 40 % and 70 % of the blade radius (around 2 % more). This effect was also discussed by Sayed et al. (2016), who explained it with a slight increase in the angle of attack in this region that is confirmed in pressure distributions at 40 % and 50 % of the blade radius in Fig. 12. A maximum c_{p} difference of around 2.5 % in the pressure side can be noticed. Considering that differently from Sayed et al. (2016), no decrease in the AOA occurs in the outer region of the blade (for these inflow conditions), no compensation of this effect occurs and, together with the increase in the rotor disk area, the increment in produced power is explained.
As in Sect. 3.1.1, the simulations including the tower and its flexibility have been repeated using BEM and two more cases at higher inflow velocities have been added. As can be seen in Fig. 13a–c, almost no tower influence can be seen in the total blade deformation, because the predicted tower top deformation by AeroDyn is very low. Therefore, almost no difference can be noticed between FMU C and RMU C in the produced torque, but only the flexibility effect that increases with the inflow velocity is apparent, leading to up to 6 % less power produced in comparison to rigid. Again, no decrease in the blade–tower passage effect can be noticed by 6.1 m s^{−1}, rather only its increase at high velocity. Differently from CFD, the predicted torque using BEM in the flexible case is always lower than in the rigid case, and the curves show less oscillation than in CFD because of the lack of timedependent 3D effects that BEM cannot capture.
3.1.3 FMU vs. FMT
Figure 14 shows isosurfaces of the λ_{2} criterion for both inflow cases. The interaction can be seen between the nearwake vortices and the Kármán vortex street of the tower. The tower faces not only the turbulence of the flow but also the wake generated by the blades, resulting in a strongly turbulent flow and oscillations in the computed loads.
The comparison of the tip deformations in flapwise and edgewise directions and the torsion can be seen in Fig. 15. The FMU case already reaches a periodic steady state after two revolutions, oscillating flapwise with an average of 2.45 % of the blade length. The same convergence trend can be seen for the edgewise deformation and for the torsion, both of them almost negligible. All three are oscillating according to the rotational frequency.
The flap and torsion deformations are mostly affected by the presence of turbulence. Especially in the flap direction, five major peaks in 10 revolutions can be observed where the maximum deformation is around 3.1 % of the blade length, which is 47 % higher than the maximum in the uniform case. At the same time, the minimum flapwise displacement, which is not due to the tower passage, is 30 % lower than in the uniform case. For the torsion deformations, the turbulence is mostly affecting the minimum, which for FMU is −0.008^{∘}, while it is −0.09^{∘} for FMT. In the defined coordinate system, a negative torsion moves the trailing edge more downwind. The edgewise displacement, although in both cases oscillating around a mean value of 0.22 %, has higher values for the first eight minima of FMT.
This can be explained by the tower top deformations in the flapwise direction in Fig. 16. In FMT the tower displacement is always smaller than in the FMU, and the tower deflection has an additional tilting effect on the rotor and consequently on the gravitational forces. After the eighth revolution, the tower top shows larger peaks in FMT than in FMU, leading to the opposite effect of a smaller peak in the edgewise deformation.
The spectra of the deformations are depicted in Fig. 17, where the rotor frequency together with the higher harmonics is marked by a symbol. High amplitudes of the harmonics of the rotor frequency can be seen in the flapwise direction, where the first one is particularly strong. Additionally, it can be recognized that due to the inflow turbulence in FMT, the higher harmonics of the rotor frequency are obscured in the broadband of the spectrum. In the edgewise direction, which is mostly influenced by gravitation and not by aerodynamics, no strong increase can be seen for the rotor frequency, and the same happens for the torsion. On the other hand, the broadband has higher amplitudes in FMT than in FMU.
The effect of the tower can again be recognized in both FMU and FMT with a delay of around 20^{∘}, where a sudden drop in the tip deformations can be seen in Fig. 15. Nevertheless this drop is almost negligible in comparison to the total affecting oscillation.
The loads resulting from the abovedescribed deformations of the FMT case are shown in Fig. 18 (the FMU case has already been discussed in Sect. 3.1.2). Independently of the rigidity of the structure, the turbulence leads to a much higher amplitude in the oscillation of the loads in comparison to FMU as seen in Fig. 10. In fact, the torque M_{x} fluctuates between 140 and 10 kNm, while in FMU it ranges between 86 and 72 kNm. Due to this high oscillation, the blade–tower passage can hardly be recognized. Unlike in the FMU case, the addition of flexibility does not have marked consequences either in thrust or in torque. Some peaks are increased in the flexible case, e.g., in both thrust and torque at 250, 315, 700 and 1000^{∘}. Averaging the result in time, the torque is increased by 2.5 % (against 3.5 % in the uniform case) due to flexibility. As for the blade–tower passage, the fluctuation inducted by the turbulence is the predominant source of oscillation; the flexibility represents only a secondary cause. This is valid only for the present case, where the inflow velocity and therefore the consequent deformations are small. In Fig. 19, the sectional loads averaged over the same revolution for both rigid and coupled conditions are plotted. It can be seen that although the shape of F_{y} has changed between 30 % and 70 % of the blade length due to the strong oscillation brought by the turbulence, almost no difference is observed by the inclusion of flexibility in comparison to the uniform case as in Fig. 11.
3.2 DEL analysis
For the fatigue loading study of the different cases considered, the necessary constants described in Sect. 2.6 have been set to Neq=10^{5}, ${S}_{\mathrm{m},\mathrm{eq}}=\mathrm{0}$ and m=11, where the last one is material dependent. The first two, as described in Hendrinks and Bulder (1995), do not influence the results, because when making fatigue comparison, it is not the absolute value but the ratio between the output from two signals that is of interest. In order to consistently compare the cycle counts, the last three revolutions of each simulation case have been considered. The chosen input signals for the following analysis are the flapwise and edgewise blade root moment, M_{y} and M_{x}, respectively. The first signal represents an unwanted action of the wind on the blade, while the second one is responsible for the power production.
The results are shown in Fig. 20, and switching in BMU from rigid (R) to coupled (C) doubles the DEL, independently of the input variable used. It is observed in Fig. 21a that the flexibility mainly increases the number of small cycles of the signal (fluctuations) and adds a few cycles with higher amplitude. In the case of FMU, already in rigid, DEL is increased by 7 times in comparison to BMU, due to the tower passage and this effect is more pronounced using M_{y} as input. It is interesting to observe that in this case, the coupling has almost no effect on the total damage. This is because, as shown in Sect. 3.1, the flexibility has two opposing influences on the loads: on the one hand the increase in the oscillations and their mean value and on the other hand the decrease in the blade–tower passage effect. These two effects almost counteract each other leading in total to a comparable value of fatigue.
Switching the FMT case from rigid to flexible increases the DEL, because, as seen in Fig. 21c, the flexibility adds a few more small cycles but no big cycles that are completely dominated by the impact of turbulence. Independently from the chosen input, the addition of turbulence drastically increases the fatigue. Far fewer cycles are detected by the rainflow counting, but they all have an amplitude larger than the largest cycles in FMU and BMU.
Finally, the ability of BEM in predicting the fatigue loading for the BMU and FMU cases is discussed. As can be seen in Fig. 22a, BEM predicts slightly higher fatigue for BMU using M_{x} as input signal than in CFD, and that is because, as prescribed in Sect. 3.1.1, the BEM model also presents a tilt angle in the BMU case (differently from CFD), leading to a sinusoidal oscillation of the forces. That means that although the CFD calculations present many more smaller cycles due to unsteady 3D effects, the DEL is mostly affected by the big ones. The same impact, but more pronounced, can be seen in BMU using M_{y} as input signal. This shows that modeling the turbine as a single blade in CFD when a tilt is given can lead to a high underevaluation of the fatigue.
The FMU case is different (no tilt modeling problem occurs), where, for both rigid and coupled and for both chosen input signals, BEM predicts higher fatigue than CFD. The difference between the rigid and coupled case remains the same as that predicted by CFD (so almost none), but the single values are almost 2 times the ones from CFD. The reason for this can be explained by looking at the cycle count in Fig. 22c. Although BEM predicts a smaller number of short cycles than CFD, cycles with around 25 kNm appear, influencing mostly the fatigue calculation. Those cycles represent the blade–tower passage, whose effect is shown to be overestimated by AeroDyn in comparison to CFD and therefore leads to higher DEL values.
In the present work, different computational fluid dynamics (CFD) models ranging from a single blade to the complete turbine including the nacelle and tower of the DANAERO turbine rotor were generated and coupled to a multibody dynamics (MBD) structural model of the same turbine, by means of a loose (explicit) coupling. The aeroelastic response of the reference turbine was calculated by the use of models increasing their complexity and fidelity in order to recognize differences and deviations connected to modeling approaches in which computational and preprocessing costs strongly differ. The effects of turbulent inflow conditions were analyzed in comparison to uniform inflow, considering both a rigid and a completely elastic wind turbine model. Additionally, a blade element momentum (BEM) model of the turbine was consistently generated and assessed against the CFD results. In this way it was possible to consider additional uniform inflow cases to determine the generalization level of the results. The objective of this study was to identify the impact and interaction of the different components and modeling approaches on the transient loads and on the damage equivalent loading (DEL) of the blade only. This was evaluated taking into account the flapwise and edgewise blade root moment at the rotor center. The major results of this study can be summarized in the following:

A highfidelity fluid–structure interaction (FSI) model of the DANAERO wind turbine has been generated and validated in comparison to experimental results.

Modeling the turbine as a single blade instead of entirely leads to only around 1 % to 2 % difference in the average quantities (sectional loads, average torque and deformations). Differently, the resulting DEL increases from BMU (blade only in uniform inflow) to RMU (entire turbine with flexible blades in uniform inflow) by up to 12 times due to the additional large cycles induced by the tower passage and because of the consideration of the tilt angle that leads to a sinusoidal oscillation of the loads, as shown by the BEM calculations.

The introduction of flexibility in BMU increases the DEL because of more load oscillations, which in FMU (entire turbine with uniform inflow) are balanced by a reduction in the tower effect. That is why the DEL was shown to not be affected by flexibility in this case.

When the entire turbine is computed as flexible, a slight increase in the torque is found in comparison to the rigid case at the computed low inflow velocity, due to the increase in the rotor disk area and a slightly increase in the angle of attack (AOA).

BEM shows in general a good agreement with CFD in evaluating the average quantities, although an overestimated tower effect is predicted (with the standard tower model implemented in the AeroDyn version coupled to Simpack) with a direct impact on the DEL evaluation. Additionally, CFD shows a decrease in the tower effect with the introduction of flexibility, which BEM does not show.

Comparing uniform and turbulent inflow, the spectra of the blade tip deformations show that the turbulence increases the amplitude of the broadband while obscuring the higher harmonics of the rotor frequency.

Independently of the rigidity of the turbine, turbulence leads to a much higher amplitude in the load oscillations, in which the tower passage becomes only a neglectable effect. This has a direct impact on the DEL of the blade that increases by up to 11 times in comparison to FMU. Flexibility is indeed additionally increasing the fatigue but much less so than in comparison to what turbulence does, showing that this is the main factor influencing the DEL calculation.
In general it can be concluded that, in the computed cases, turbulence is shown to be the most important factor influencing the DEL of the single blade, more than flexibility, which played in comparison only a marginal role for this specific case where the rotor radius is only 40 m long. Note that when the rotor size increases, the effect of flexibility may play a greater role. Also, the modeling of the turbine as a single blade strongly underestimates the DEL even if CFD is used. On the other hand, a singleblade model (that is much cheaper than a full CFD model of the turbine) is realized to give valid results when just the averaged deformations and loads in uniform inflow are of interest and the predicted tower top deformations are low (as for the low inflow velocity studied in this paper). AeroDyn overestimates the blade–tower effect in comparison to CFD, leading to higher fatigue values, but excluding this overestimated tower effect, BEM can be employed to give useful conclusions regarding the effect of flexibility on fatigue for the uniform inflow conditions under which it has been used.
The raw data of the simulation results can be provided by contacting the corresponding author.
GG generated a part of the CFD model, generated the MBD model, ran the coupled simulations, and performed the postprocessing and analysis. GB generated a part of the CFD model. TL and EK supported the research, defined and supervised the work, and revised the manuscript.
The authors declare that they have no conflict of interest.
The authors gratefully acknowledge the DANAERO Consortium for providing the geometry and structural data. They additionally acknowledge Simpack for providing the user licenses and the funders of the project WINSENT (code number 0324129); the Federal Ministry for Economic Affairs and Energy (BMWi); and the Ministry of the Environment, Climate Protection and the Energy Sector BadenWürttemberg under the funding number L75 16012, under which project improvements on the simulation chain were performed. Computer resources were provided by the Gauss Centre for Supercomputing and Leibniz Supercomputing Centre under grant pr94va. Additionally particular thanks are given to the DLR and SWE, University of Stuttgart, for the productive discussions that helped in improving the structural model of these simulations. Finally the authors would like to acknowledge Aimable Uwumukiza for his effort in correcting the language and presentation of this work in its first submission.
This article is part of the special issue “Wind Energy Science Conference 2019”. It is a result of the Wind Energy Science Conference 2019, Cork, Ireland, 17–20 June 2019.
This openaccess publication was funded
by the University of Stuttgart.
This paper was edited by Katherine Dykes and reviewed by David Verelst and two anonymous referees.
Aagaard Madsen, H., Bak, C., Schmidt Paulsen, U., Gaunaa, M., Fuglsang, P., Romblad, J., Olesen, N. A., Enevoldsen, P., Laursen, J., and Jensen, L.: The DANAERO MW experiments, Final report, Denmark, 2010. a
AeroDyn Theory Manual: National Renewable Energy Laboratories, available at: https://www.nrel.gov/docs/fy05osti/36881.pdf (last access: 17 September 2020), 2005. a
Bazilevs, Y., Hsu, M. C., Kiendl, J., Wüchner, R., and Bletzinger, K. U.; 3d simulation of wind turbine rotors at full scale. prt ii: Fluid–structure interaction modeling with composite blades, Int. J. Numer. Meth. Fluids, 65, 236–253, 2011. a
Boorsma, K., Wenz, F., Aman, M., Lindenburg, C., and Kloosterman, M.: TKI WoZ VortexLoads, Final report TNO 2019 R11388, TNO, Petten, available at: http://resolver.tudelft.nl/uuid:f4d6218af19c4246a3b0d0ac64af7618 (last access: 11 January 2021), 2019. a
DANAERO project: Experimental Rotor and Airfoil Aerodynamics on MW Wind Turbines, available at: http://www.danaero.dtu.dk/, last access: 28 January 2020. a
Dose, B., Rahimi, H., Herráez, I., Stoevesandt, B., and Peinke, J.: Fluidstructure coupled computations of the nrel 5 mw wind turbine by means of cfd, Renew. Energy, 129, 591–605, 2018. a
Guma, G., Bangga, G., Jost, E., Lutz, T., and Kraemer, E.: Consistent 3D CFD and BEM simulations of a research turbine considering rotational augmentation, J. Phys.: Confer. Ser., 1037, 022024, https://doi.org/10.1088/17426596/1037/2/022024, 2018. a
Hansen, M. H., Thomsen, K., Fuglsang, P., and Knudsen, T.: Two methods for estimating aeroelastic damping of operational wind turbine modes from experiments, Wind Energy, 9, 179–191, 2006. a
Hansen, M. O. L., Sørensen, J. N., Voutsinas, S., Sørensen, N., and Madsen, H. A.; State of the art in wind turbine aerodynamics and aeroelasticity, Prog. Aerosp. Sci., 42, 285–330, 2006. a
Heinz, J. C., Sørensen, N. N., and Zahle, F.: Fluid–structure interaction computations for geometrically resolved rotor simulations using cfd, Wind Energy, 19, 2205–2221, 2016. a
Hendriks, H. B. and Bulder, B. H.: Fatigue Equivalent Load Cycle Method, ECN Publicaties, ECN, available at: https://www.osti.gov/etdeweb/biblio/191445 (last access: 11 January 2021), 1995. a, b
Hsu, M. C. and Bazilevs, Y.: Fluid–structure interaction modeling of wind turbines: simulating the full machine, Comput. Mech., 50, 821–833, 2012. a
IEA Task 29: Task 29 Work Plan and Objectives, available at: https://community.ieawind.org/task29/29workplan, last access: 28 January 2020. a, b
Jameson, A., Schmidt, W., and Turkel, E.: Numerical solutions of the Euler equations by finite volume methods using Runge–Kutta timestepping schemes, AIAA Paper 811259, available at: https://www.researchgate.net/publication/247576565_Solutions_of_the_Euler_Equations_by_Finite_Volume_Methods_Using_RungeKutta_TimeStepping_Schemes. Last access: 15/01/2020 (last access: 11 January 2021), 1981. a
Jassmann, U., Berroth, J., Matzke, D., Schelenz, R., Reiter, M., Jacobs, G., and Abel, D.: Model predictive control of a wind turbine modelled in Simpack; J. Phys.: Conf. Ser., 524, 012047, https://doi.org/10.1088/17426596/524/1/012047, 2014. a
Jeong, M. S., Yoo, S. J., and Lee, I.: Aeroelastic analysis for large wind turbine rotor blades, in: 52nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Structures, Structural Dynamics, and Materials and Colocated Conferences, AIAA, 7 April 2011 Denver, Colorado, 9–14, 2011. a
Kim, Y., Weihing, P., Schulz, C., and Lutz, T.: Do turbulence models deteriorate solutions using a nonoscillatory scheme?, J. Wind Eng. Indust. Aerodynam., 156, 41–49, 2016. a
Klein, L., Gude, J., Wenz, F., Lutz, T., and Krämer, E.: Advanced computational fluid dynamics (CFD)–multibody simulation (MBS) coupling to assess lowfrequency emissions from wind turbines, Wind Energ. Sci., 3, 713–728, https://doi.org/10.5194/wes37132018, 2018. a, b
Kowarsch, U., Keßler, M., and Krämer, E.: High order CFDsimulation of the rotorfuselage interaction, in: 39th European Rotorcraft Forum, 3–6 September 2013, Moscow, Russia, 2013. a
Li, Y., Castro, A. M., Martin, J. E., Sinokrot, T., Prescott, W., and Carrica, P. M.: Coupled computational fluid dynamics/multibody dynamics method for wind turbine aeroservoelastic simulation including drivetrain dynamics, Renew. Energy, 101, 1037–1051, 2017. a
Luhmann, B., Seyedin, H., and Cheng, P.W.: Aerostructural dynamics of a flexible hub connection for load reduction on twobladed wind turbines, Wind Energy, 20, 521–535, 2017. a
Madsen, H. A., Riziotis, V., Zahle, F., Hansen, M., Snel, H., Grasso, F., Larsen, T., Politis, E., and Rasmussen, F.: Blade element momentum modeling of inflow with shear in comparison with advanced model results, Wind Energy, 15, 63–81, 2012. a
Madsen, H. A., Sørensen, N. N., Bak, C., Troldborg, N., and Pirrung, G.: Measured aerodynamic forces on a full scale 2 MW turbine in comparison with EllipSys3D and HAWC2 simulations, J. Phys.: Conf. Ser., 1037, 022011, https://doi.org/10.1088/17426596/1037/2/022011, 2018. a
Mann, J.: The spatial structure of neutral atmospheric surfacelayer turbulence, J. Fluid Mech., 273, 141–168, 1994. a, b
Masarati, P. and Sitaraman, J.: Coupled cfd/multibody analysis of nrel unsteady aerodynamic experiment phase vi rotor, in: AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper, Orlando, FL, 2011. a
Menter, F. R.: Two equation eddyviscosity turbulence models for engineering applications, AIAA J., 32, 1598–1605, 1994. a
Raddatz, J.: The blockstructured RANS solver FLOWer, DLR, Institute of Aerodynamics and Flow Technology, in: MEGAFLOW – Numerical Flow Simulation for Aircraft Design, Springer, Berlin, Heidelberg, 27–44, 2009. a
Rahimi, H., Schepers, J. G., Shen, W. Z., García, N. R., Schneider, M. S., Micallef, D., Simao Ferreira, C. J., Jost, E., Klein, L., and Herráez, I.: Evaluation of different methods for determining the angle of attack on wind turbine blades with CFD results under axial inflow conditions, Renew. Energy, 125, 866–876, 2018. a
Riziotis, V. A., Voutsinas, S. G., Politis, E. S., Chaviaropoulos, E. S., Hansen, A. M., Madsen Aagaard, H., and Rasmussen, F.: Identification of structural nonlinearities due to large deflections on a 5 MW wind turbine blade, in: EWEC, Vindenergi, Scientific proceedings, European Wind Energy Conference and Exhibition, 31 March–3 April 2008, Brussels, 9–14, 2008. a
Sayed, M., Lutz, T., Krämer, E., Shayegan, S., Ghantasala, A., Wüchner, R., and Bletzinger, K.U.: High fidelity CFDCSD aeroelastic analysis of slender bladed horizontalaxis wind turbine, J. Phys.: Conf. Ser., 753, 042009, https://doi.org/10.1088/17426596/753/4/042009, 2016. a, b, c
Schepers, J.: Latest results from the EU project AVATAR: Aerodynamic modelling of 10 MW wind turbines, J. Phys.: Conf. Ser., 753, 022017, https://doi.org/10.1088/17426596/753/2/022017, 2016. a
Troldborg, N., Bak, C., Madsen, A. H., and Skrzypinski, W. R.: DANAERO MW: Final Report, DTU Wind Energy, 2013. a
Troldborg, N., Sørensen, J. N., Mikkelsen, R., and Sørensen, N. N.: A simple atmospheric boundary layer model applied to large eddy simulations of wind turbine wakes, Wind Energy, 17, 657–669, https://doi.org/10.1002/we.1608, 2014. a, b
Yu, D. O. and Kwon, O. J.: Timeaccurate aeroelastic simulations of a wind turbine in yaw and shear using a coupled CFDCSD method, J. Phys.: Conf. Ser., 524, 012046, https://doi.org/10.1088/17426596/524/1/012046, 2014. a