Aeroelastic analysis of wind turbines under turbulent inflow conditions

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 high-fidelity model. A time-accurate unsteady fluid–structure interaction (FSI) coupling is used between a computational fluid dynamics (CFD) code for the aerodynamic response and a multi-body 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.

1 Introduction

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 non-stationary 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 multi-body 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, high-fidelity 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 low-order arbitrary Lagrangian–Eulerian variational multiscale (ALE-VMS) flow solver and a structural solved based on a non-uniform 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 FEM-based beam solver BeamFOAM. A CFD–MBS coupling between the unsteady Reynolds-averaged 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 BEM-based. In the EU project AVATAR (Schepers, 2016) it was shown that BEM-based calculations against high-fidelity 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 BEM-based and CFD-based calculations have been compared with each other and to experimental results.

Within the scope of the present study, a highly accurate CFD-based 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 one-third 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 high-fidelity approaches. For the uniform inflow case, a comparison with BEM-based 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 high-fidelity 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 high-fidelity 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 post-processing of the different simulations, using two different time-varying input variables.

2 Methodology

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 post-processed and calibrated in the follow-up 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 detached-eddy-simulation (DES) finite-volume solver for structured meshes. The present simulations are run using the shear-stress-transport (SST) k–ω model according to Menter (1994), using a fully turbulent boundary layer. Two different spatial discretization schemes are available, a second-order central cell-centered Jameson–Schmidt–Turkel (JST) (Jameson et al., 1981) and a fifth-order weighted essentially non-oscillatory (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 time-stepping scheme is an artificial five-stage Runge–Kutta scheme, and multi-grid 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 (cylinder-based) from provided geometrical properties. Meshes are generated by the use of the commercial software Pointwise in combination with in-house scripts. All components have been meshed ensuring $y^{+}\le \mathrm{1}$ in the boundary layer region. The blades are meshed in an O-mesh 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:

1. one-third model (BMU) of the rotor (only one blade) suited for uniform inflow conditions,

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

3. 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 three-bladed 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.

Figure 1 Details of the meshes and boundary conditions for BMU (a) and FMU and FMT (b).

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⁻¹ 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 detached-eddy 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 multi-body 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 multi-body 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 (three-dimensional beam structures in Simpack, described by a node-based nonlinear finite-differences 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.

Table 1 Comparison of natural frequencies between the measured ones and those computed by Simpack; single blade above and full turbine below.

Figure 2 Normal (on the left) and tangential (on the right) sectional load in comparison for a single rigid blade 3D CFD vs. AeroDyn.

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 BEM-based simulations using the aeroelastic tool HAWC2 to the high-fidelity 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.

Table 2 Computed cases with uniform inflow in BEM. The first case is the one also computed with CFD.

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⁻¹, 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.

Figure 3 Explicit coupling strategy.

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 SuperMUC-NG supercomputer at the Leibniz-Rechenzentrum 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.

Table 3 Computed cases with inflow condition, CFD modeled structures and flexibility.

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 time-varying 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

$$\begin{array}{}\text{(1)}& \mathrm{DEL}=S_{\mathrm{r},\mathrm{eq}}={\left(\sum _{i=\mathrm{1}}^n{\displaystyle \frac{{\left(S_{\mathrm{r},i}\cdot \frac{S_{\mathrm{u}}-S_{\mathrm{m},\mathrm{eq}}}{S_{\mathrm{u}}-S_{\mathrm{m},i}}\right)}^m}{\mathrm{Neq}}}\right)}^{\frac{\mathrm{1}}{m}},\end{array}$$

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 Results

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⁻¹); 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.

Figure 4 Comparison of experimental normal to the chord loads (F_N) in (a)–(c) and tangential to the chord loads (F_T) in (d)–(f) for three different radial sections (r=13, r=19, r=37 m) over the blade. The blue line represents the full turbine with flexible blades. The red line represents a rigid rotor without tower but a turbulent inflow with the same TI as in the experiments. Thin grey and pink lines represent the data per revolution for the experiments and “CFD Turb”, respectively.

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 only-rotor 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.

Figure 5 Tip deformations calculated with CFD at 6.1 m s⁻¹ comparing BMU coupled (C) vs. RMU coupled (C). Out-of-plane deformation in (a), in-plane deformation in (b) and torsion in (c).

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.

Figure 6 Thrust and torque calculated with CFD at 6.1 m s⁻¹ comparing BMU vs. RMU, both rigid (R) and coupled (C).

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.

Figure 7 Sectional loads for BMU vs. RMU: comparison between rigid (R) and coupled (C) calculated with CFD at 6.1 m s⁻¹.

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⁻¹ 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⁻¹. 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.

Figure 8 Aeroelastic calculations using BEM as aerodynamic model. Tip deformations in flapwise direction BMU vs. RMU: 6.1 m s⁻¹ in (a), 9.0 m s⁻¹ in (b) and 13 m s⁻¹ in (c). Torque (M_x) generated by one-blade BMU vs. RMU: 6.1 m s⁻¹ in (d), 9.0 m s⁻¹ in (e) and 13 m s⁻¹ in (f).

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 out-of-plane deflection.

Figure 9 Tip deformations calculated with CFD at 6.1 m s⁻¹ comparing RMU vs. FMU.

Figure 10 Thrust and torque calculated with CFD at 6.1 m s⁻¹ comparing RMU vs. FMU, both rigid (R) and coupled (C).

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.

Figure 11 Sectional loads for FMU: comparison between rigid (R) and coupled (C) calculated with CFD at 6.1 m s⁻¹.

Figure 12 Pressure distributions for FMU rigid and coupled in comparison calculated with CFD at 6.1 m s⁻¹.

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⁻¹, 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 time-dependent 3D effects that BEM cannot capture.

Figure 13 Aeroelastic calculations using BEM as aerodynamic model. Tip deformations in flapwise direction RMU vs. FMU: 6.1 m s⁻¹ in (a), 9.0 m s⁻¹ in (b) and 13 m s⁻¹ in (c). Torque (M_x) generated by one-blade RMU vs. FMU: 6.1 m s⁻¹ in (d), 9.0 m s⁻¹ in (e) and 13 m s⁻¹ in (f).

3.1.3 FMU vs. FMT

Figure 14 shows iso-surfaces of the λ₂ criterion for both inflow cases. The interaction can be seen between the near-wake 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.

Figure 14 Visualization of the λ₂ criterion.

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.

Figure 15 Tip deformations comparing FMU vs. FMT calculated with CFD.

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.

Figure 16 Tower top deformation in flapwise direction calculated with CFD.

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.

Figure 17 Spectra of the deformations comparing FMU vs. FMT.

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 above-described 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.

Figure 18 Global loads in FMT: comparison between rigid (R) and coupled (C) calculated with CFD.

Figure 19 Sectional loads in FMT: comparison between rigid (R) and coupled (C) calculated with CFD.

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⁵, $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.

Figure 20 DEL calculation based on CFD for the different cases using in (a) M_x and in (b) M_y. Comparison between rigid (R) and coupled (C).

Figure 21 Comparison of number of cycle counts to load ranges using M_x from CFD as input.

Figure 22 DEL calculation using BEM: results for M_x in (a) and for M_y in (b). Comparison between rigid (R) and coupled (C). Cycle count in comparison to load ranges for FMU using as input M_x in (c).

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.

4 Conclusions

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 multi-body 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 pre-processing 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:

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

2. 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.

3. 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.

4. 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).

5. 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.

6. 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.

7. 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 single-blade 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.