Articles | Volume 7, issue 6
Wind Energ. Sci., 7, 2201–2213, 2022
Wind Energ. Sci., 7, 2201–2213, 2022
Research article
07 Nov 2022
Research article | 07 Nov 2022

Multiple limit cycle amplitudes in high-fidelity predictions of standstill wind turbine blade vibrations

Multiple limit cycle amplitudes in high-fidelity predictions of standstill wind turbine blade vibrations
Christian Grinderslev, Niels Nørmark Sørensen, Georg Raimund Pirrung, and Sergio González Horcas Christian Grinderslev et al.
  • Department of Wind Energy, Technical University of Denmark, Risø Campus, 4000, Roskilde, Denmark

Correspondence: Christian Grinderslev (


In this study, vortex-induced vibrations (VIVs) on the IEA 10 MW blade are investigated using two methodologies in order to assess strengths and weaknesses of the two simulation types. Both fully coupled fluid–structure interaction (FSI) simulations and computational fluid dynamics (CFD) with forced motion of the blade are used and compared. It is found that for the studied cases with high inclination angles, the forced-motion simulations succeed in capturing the power injection by the aerodynamics, despite the motion being simplified. From the fully coupled simulations, a dependency on initial conditions of the vibrations was found, showing that cases which are stable if unperturbed might go into large VIVs if provoked initially by, for instance, inflow turbulence or turbine operations. Depending on the initial vibration amplitudes, multiple limit cycle levels can be triggered, for the same flow case, due to the non-linearity of the aerodynamics. By fitting a simple damping model for the specific blade and mode shape from FSI simulations, it is also demonstrated that the equilibrium limit cycle amplitudes between power injection and dissipation can be estimated using forced-motion simulations, even for the multiple stable vibration cases, with good agreement with fully coupled simulations. Finally, a time series generation from forced-motion simulations and the simple damping model is presented, concluding that CFD amplitude sweeps can estimate not only the final limit cycle oscillation amplitude, but also the vibration build-up time series.

1 Introduction

As wind turbines become larger and more flexible, the risk of serious instabilities between flow and structure increase. An example of these instabilities is the phenomenon known as vortex-induced vibrations (VIVs). Here, the shedding of separated flow leads to vibrations of the structure. This does not necessarily lead to any structural issues; however, in some cases the so-called lock-in phenomenon can occur. In that situation, the vortex shedding frequency locks in with the structural frequency in which the blade is vibrating. If this occurs over a sufficient part of the blade span, the corresponding vibrations can become violent and threaten the structural integrity of the blades. The phenomenon of VIVs has been extensively studied for simple geometries such as cylinders, but the number of studies relating to wind turbine blades is quite limited. Main publications of VIVs on wind turbine blades are the ones of Horcas et al. (2019, 2020) and Heinz et al. (2016b), who have studied VIVs for large wind turbine blades (AVATAR 10 MW, Ceyhan et al.2016; IEA 10 MW, Bortolotti et al.2019; and DTU 10 MW RWT, Bak et al.2013). All the aforementioned works were strictly numerical and relied on fluid–structure interaction (FSI) simulations, involving the communication of a fluid solver and a structural solver during runtime, through a loose coupling scheme. These previous studies have identified VIV risk zones depending on the angles between flow and blade as well as the velocity magnitude. It is found that the risk of VIVs is very dependent on the so-called pitch P and inclination I angles, as defined in Fig. 1. Assuming that the blade points upwards, the pitch angle P is here defined as the relative horizontal angle between flow and the root section chord. P=0 when the horizontal wind projection strikes directly the leading edge of the untwisted blade, P=90 when striking the pressure side perpendicularly, and P=180 when striking the trailing edge. I is the relative vertical angle between the inflow wind and the plane intersecting the root section (i.e., the xz plane in Fig. 1). I is positive, as shown in Fig. 1, when the spanwise component of the wind flow is from tip to root, and it is zero when the wind strikes the blade perpendicularly. In particular, inclination angles of  50–60 and pitch angles of  80–100 are in the risk zone of high-amplitude VIVs for the studied blade. On top of this, studies have been conducted on the impact of geometry alterations on VIVs such as flap deformations (Horcas et al.2019) and alterations of the blade tip (Horcas et al.2020). These have shown that VIVs are sensitive to the geometrical alterations; however, these seem to only move the risk zones rather than alleviate the phenomenon. In all cases, the first edgewise vibration mode was found to be excited for these kind of vibrations, at least for the studied conditions. This mode is likewise the focus on the present study.

Figure 1Definition of inclination and pitch angles. Reproduced from Horcas et al. (2022).

The main objective of the present study is the exploration of the forced-motion approach as a complementary analysis tool for the study of wind turbine blade VIVs. The forced-motion technique has been extensively applied to academic configurations such as cylinders (Williamson and Govardhan2004). In this approach the surface displacements are imposed, usually through a simple harmonic function based on the eigenvalues and eigenvectors of the structural problem, rather than being the result of an aeroelastic interaction computed during runtime (Bertagnolio et al.2002). For a given inflow, if the imposed deflection or frequency of motion do not correspond to a physical aeroelastic response, the computed flow and transferred power could substantially differ from the results of a FSI simulation. For instance, it is well known that the lock-in region becomes narrower when simulating low amplitudes of motion, both for cylinders (Koopmann1967; Blevins2001; Placzek et al.2009) and airfoils (Meskell and Pellegrino2019; Hu et al.2021). Additionally, the forced motion disregards the structural damping, which can have an important effect on the VIV phenomenon (Lee and Bernitsas2011; Blevins and Coughran2009). However, if assumptions on the loading and on the response can be made, the forced-motion strategy is able to estimate the amplitude of vibrations of a given inflow (Bearman2011), which is often the targeted unknown in a design process. This is usually done by performing a series of forced-motion simulations of varying amplitude and subsequently building a reduced-order model of the aeroelastic system (Lupi et al.2021; Zhang et al.2022). Additionally, the ability to choose the vibration pattern and amplitudes in the forced motion yields opportunities to investigate the aerodynamics for vibration cases outside those found by a coupled approach. Another benefit of the forced-motion approach is the fact that neither time-marching structural solver nor coupling framework is needed, simplifying the simulations. As the motion of the structure is prescribed, the needed simulation time decreases significantly as no time for the motion build-up is needed.

The current study is part of the PRESTIGE project (EUDP2020), looking into stability of wind turbines in various aspects. In this particular work, the forced-motion and FSI simulations of the IEA 10 MW (Bortolotti et al.2019) blade are performed and compared for several conditions. First, simulations are devoted to assess the performance and particularities of both approaches. This initial comparison is followed by two detailed studies, aiming at analyzing the impact of (i) the structural damping and (ii) the initial conditions.

2 Methodology

2.1 Flow solver

To calculate the fluid flow, the computational fluid dynamics (CFD) code EllipSys3D (Michelsen1994, 1992; Sørensen1995) is used. This code solves the incompressible Navier–Stokes equations through Reynolds-averaged Navier–Stokes (RANS) simulations, large eddy simulations (LESs) or hybrids of these. In this particular study, the IDDES model (Gritskevich et al.2012) is utilized to effectively resolve the shed vortices but without having to resolve the viscous sub-layer in LES. The code utilizes structured grids in curvilinear coordinates using the finite-volume method. This enables highly scalable computations through MPI with multi-block decomposition and grid sequencing. Convection is solved through various schemes such as central difference, second-order upwind, QUICK or combinations of said schemes. Here, a combination of the QUICK and a fourth-order central difference scheme is used dependent on the considered region being in a LES or RANS region of the IDDES model (Shur et al.2008). For the velocity–pressure coupling, the SIMPLE algorithm or variations hereof are used along with Rhie–Chow interpolations to avoid odd–even pressure decoupling. Deformations imposed on the CFD surface mesh are propagated into the volume grid trough a blend factor. This factor is based on the grid line normal distance to the surface and ensures that mesh cells close to the surface are moved as rigid cells, while further from the surface the cells are deformed based on the distance. The solver has been extensively used and validated through multiple studies of wind turbines; see, for instance, the Mexico project (Bechmann et al.2011; Sørensen et al.2016), the Phase VI NREL rotor simulations (Sørensen and Schreck2014; Sørensen et al.2002), and recent comparisons with wind tunnel tests of a curved blade tip (Barlas et al.2021) and with the full-scale DanAero measurements (Madsen et al.2018; Grinderslev et al.2020b; Schepers et al.2021).

2.1.1 Mesh

The structured grid used for the CFD solver is similar to the mesh used in Horcas et al. (2020), where a sensitivity study of mesh resolution and time step was conducted. It is a body-fitted grid that, for simplicity, only accounts for the geometry of a single blade (and not for the full rotor). The surface mesh consists of 160 cells in spanwise and 256 in chordwise direction (see Fig. 2). The volume mesh is extruded from the grid surface into a sphere surrounding the blade using the mesh tool HypGrid3D (Sørensen1998). The diameter of the sphere corresponds to  15 blade lengths. A total of 256 cells are used in the normal direction from the surface with an initial grid size of 1×10-6 to ensure a y+ of less than 1. The mesh is refined such that cells do not expand rapidly in the region close to the blade where vortices shed from the blade are of interest, and a good cell resolution is needed by the IDDES turbulence model. Further out, the cell size is rapidly expanded toward the boundaries of the sphere. These boundaries are depicted in Fig. 3, together with a detail of the near-wall volume mesh.

Figure 2Surface grid. Only every second grid line is shown.


Figure 3Volume grid near the surface at 50 % blade span and sphere domain (intersected). Only every second grid line is shown.


2.2 Coupling framework

To capture the interaction between fluid and structure response, two methodologies are investigated: forced-motion (FM) and fully coupled fluid–structure interaction (FSI) simulations. These methodologies are also commonly known as one-way and two-way coupling, respectively. For both simulation types, it is here for simplicity assumed that the blade is clamped at the root. This means that singe-blade modes can be captured, but full rotor/turbine modes are omitted.

2.2.1 Structural solver

To solve the structural response for two-way coupling HAWC2 (Larsen and Hansen2007; Madsen et al.2020; Pavese et al.2015) is used. HAWC2 is a multi-physics code widely used as an aeroelastic solver using multi-body dynamics (MBD) finite elements for the structure and blade element momentum (BEM) theory for aerodynamics. Here, only the structural module is utilized. The blade is represented through a number of finite-element Timoshenko beams, which are connected through constraint equations to enable considerations of non-linearities. Structural damping is included through Rayleigh damping using the proportionality parameters of the official model of the IEA 10 MW wind turbine HAWC2 model (Bortolotti et al.2019). The blade in the present setup is divided into 18 beam elements. Aerodynamic loads are imposed through so-called aerosections distributed along the blade. In this setup 101 aerosections, uniformly distributed along the blade, are used.

2.2.2 Fully coupled simulations

For the two-way-coupled simulations (FSI), the DTU coupling (Heinz et al.2016a; García Ramos et al.2020) is used, which through a Python framework communicates loading and structural response between the flow- and structural solvers. Deformations predicted in HAWC2 are imposed in the EllipSys3D CFD where the corresponding load is found and communicated back to HAWC2 to correct the structural response before moving to the next time step. This is done in a loose coupling, meaning that only one communication step is done per global time step, while each solver separately uses inner iterations to converge during each time step. Usually, the CFD solver will set the necessary time steps of the simulations, which in the present setup was set to 0.006 s per time step, which corresponds to 250 steps per motion cycle for the edgewise mode shape with 0.67 Hz frequency. The framework has been used for various studies looking at VIVs (Heinz et al.2016b; Horcas et al.2019, 2020) and for wind turbine operation in various flow conditions (Grinderslev et al.2021a, b), in which the framework is described in more detail. Additionally, the framework was validated with experimental tests in Grinderslev et al. (2020a).

2.2.3 Forced-motion simulations

Forced-motion (FM) simulations are conducted directly through the CFD solver, EllipSys3D. This is done by imposing the deformation through analytical expressions of the mode shape investigated, which is the first edgewise mode. Vibrations are imposed by multiplying the mode shape with a sine function based on the first edgewise mode frequency and the amplitude desired at the tip of the blade. Along with this, a mean deflection of the blade is imposed based on the deflection found from a FSI simulation. The first edgewise mode shape has been fitted through a seventh-order polynomial to the mode shape obtained from the aeroelastic stability tool HAWCSTAB2 (Hansen2004). Both edgewise and flapwise components of the first edgewise mode are considered. The torsional component, which amounts to less than half a degree tip torsion for 1 m edgewise deflection, is omitted. This is possible because the influence of small amounts of torsion on the aerodynamics is negligible at around a 90 angle of attack. The mode shape polynomials along with the corresponding output from HAWCSTAB2 are shown in Fig. 4, and the static mean flap deflection subtracted from FSI simulations and fitted to a seventh-order polynomial are shown in Fig. 5. From earlier work (Horcas et al.2022; Riva et al.2022), it is known that, for the low wind speeds experienced here, the structural mode shape obtained in HAWCSTAB2 is very similar to the deflection shape seen in FSI.

Figure 4First edgewise mode shape of blade from HAWCSTAB2 along with polynomial fit used in forced-motion simulations.


Figure 5Static deflection imposed in forced-motion simulations and mean deformation extracted from FSI simulation with flow case P100-I50.


2.3 Cases studied

As the computations presented here are computationally heavy, three flow scenarios are picked with basis in the study of Horcas et al. (2022). One case which shows high vibrations in the mapping study of Horcas et al. (2022) is the P=100, I=50 case. Two cases at the vicinity of this case, P=95 and P=90, both with I=50, which both show no vibration in the aforementioned mapping are also included. These three cases will be identified as P90-I50, P95-I50 and P100-I50 in the remainder of this paper. All cases are simulated with a uniform inflow of 18 m s−1, and the angles are obtained by varying the inflow direction. Thereby, the blade is always pointing upwards with gravity pulling from tip to root. In this study, positive aerodynamic power is defined as injection of power into the structure, meaning that a positive value will increase the vibration amplitudes. The average aerodynamic power PAERO is found by the time-averaged cell-wise integration of the power along each of the cells of the blade surface mesh.

(1) P AERO = 1 T n cycles t = t 0 t 0 + T n cycles i = 1 N F i u ˙ i d t ,

where T is the time for one motion period, i is the cell index, N is the total number of cells, and Fi and u˙i are the force and velocity at every cell.

3 Results

The three aforementioned flow scenarios have been tested for various forced-motion amplitudes to investigate the impact on the aerodynamic power injection. In Fig. 6, the resulting aerodynamic power depending on amplitude is depicted for all three cases. As seen, the power in situation P100-I50, where vibrations were found using FSI in previous work, follows a quite smooth bell shape with positive power injection from small amplitudes until reaching an amplitude of 1.6 m. The blade would reach a limit cycle oscillation at this amplitude if no structural damping is considered. For the two other cases, which were found to be stable in FSI simulations, the initial slope of the power curves is much flatter, giving rise to very little power injection in the first 0.5 m amplitude. Therefore, structural damping will give rise to the blade being stable or having only low-amplitude vibrations. Going further up in amplitude, the power injection increases rapidly, showing that a range of more than 1 m actually has significant positive power injection, which could be higher than what is dissipated by the structural damping. The maximum amplitude with positive power injection is close to the former case around 1.6 m.

Figure 6Aerodynamic power obtained at various forced amplitudes for flow cases with 50 inclination and pitch angles of 90, 95 and 100.


Case P100-I50 was further tested in the coupled FSI framework with the baseline structural damping and without any structural damping. This leads to limit cycle oscillations (LCOs) in the edgewise direction, as presented in Fig. 7. As can be seen in the figure, both cases develop LCOs after some time but stabilize at different amplitudes, as expected. The undamped case obtains an amplitude of approximately 1.63 m, which corresponds well with the intersection with 0 W in aerodynamic power injection found for the forced-motion simulation presented in Fig. 6. The structurally damped case leads to amplitudes of  1.18 m amplitude, meaning that the power dissipated by structural damping reduces the vibration amplitude by  27 % for the studied case. If vibration magnitude is the desired outcome of a simulation, the figure above shows a limiting factor for the FSI, as the development time of the LCO is 60 simulated seconds for the baseline case and even longer for the undamped case.

Figure 7Edgewise tip deflection from coupled FSI simulation of flow case P100-I50 with and without considering structural damping.


A clear benefit of the forced-motion methodology is that the limit cycle amplitude is set in advance, meaning that no vibration build-up needs to be simulated. Thereby, the flow shedding just needs to stabilize with the prescribed motion, which is found to only take a few motion periods. This is depicted in Fig. 8, showing the convergence of the total injected aerodynamic power. This benefit, however, diminishes if data for many amplitudes are needed in order to, for instance, estimate how the vibrations build up.

Figure 8Development of aerodynamic power injection (moving average) in FM simulation. Case P100-I50 with 1.0 m amplitude.


Taking a close look at the average aerodynamic power injection along the blade span for the coupled and forced-motion cases, Fig. 9 shows the results for the flow case P100-I50. The spanwise power per cycle Pcycle is here found as the injected work averaged over the time span of n cycles as presented in Horcas et al. (2022). The work ws over n cycles is found as the integration over the time range of the sectional loads Fs and the velocity of the considered section u˙s.

(2) P cycle = w s , t 1 t 1 - t 0 = 1 t 1 - t 0 t = t 0 t 1 F s u ˙ s d t ,

where t0 is chosen to be after the establishment of the LCO and the integration is performed over n periods of motion:

(3) t 1 = t 0 + n f e .

For the FSI simulation results, the baseline damping (BL) and a zero structural damping (ZD) case were chosen and compared to the forced-motion cases with amplitudes closest to that obtained in FSI, being 1.2 and 1.57 m, respectively. The general picture shows a very good agreement in both trend and magnitude between the spanwise power injections of the two methods, which indicates that the simplified motion imposed in FM simulations are sufficient. The small differences that occur are likely due to the omission of torsion in FM simulations along with the fact that the amplitudes are not exactly alike between the FSI and FM cases. The main contributor to the total power injection is in the edgewise direction where the motion is largest and the loads interlock with the motion velocity. In the considered cases, especially the range between 40 and 60 m span of the blade introduces energy to the system, pointing to that range being the critical area where the phase between the flow shedding and blade velocity lock together.

Figure 9Averaged power injected in flapwise and edgewise directions by aerodynamics along span over one motion cycle for FSI and FM simulations. (a, b) FSI: baseline damping; FM: 1.2 m amplitude. (c, d) FSI: zero damping; FM: 1.57 m amplitude.


3.1 Structural damping fitting

In general, the coupled FSI simulations conducted in this study use the Rayleigh damping parameters of the official model of the IEA 10 MW wind turbine HAWC2 model (Bortolotti et al.2019). However, it is also interesting to assess the effect of other damping coefficients on the VIVs observed and to fit a simple relation between the amplitudes obtained and the power dissipated by damping. A set of simulations of flow case P100-I50 was conducted through FSI, using various scalings of the structural damping coefficients. A scale factor of δDAMP of [0, 0.5, 1, 2, 5, 10, 15] was used on all Rayleigh damping coefficients for the blade structural model. All cases except for δDAMP= 15 lead to VIVs, but with varying amplitudes as expected since the equilibrium between aerodynamic power injection and dissipating structural damping will depend on the amplitude. For the high structural damping of 15 times the baseline, the vibrations are removed while small vibrations still occur when having 10 times the baseline damping. Scaling the damping to half or twice the baseline damping yields only marginal effects on the steady-state amplitude with an increase (half damp) or decrease (double damp) of  0.07 m (for this reason these are not shown in Fig. 10). This finding is a result of the highly negative gradient of aerodynamic power around the seen amplitude, when inspecting Fig. 6. Since the aerodynamic power increases rapidly when lowering the amplitude only slightly, and vice versa, the doubling of power dissipated structurally will reach an equilibrium with the aerodynamics at an amplitude only a few centimeters lower. The power injected/dissipated is, however, doubled. Note also that the development time of the LCOs increases as the damping is increased. This could be beneficial as critical flow conditions need to remain for longer periods in order to develop the VIVs, which at the same time become less critical in terms of amplitudes with the increased damping.

Figure 10Edgewise tip displacements (x) for various scaling of the structural damping Rayleigh parameters.


A simple model for the power dissipated by structural damping PSTRUC based on tip motion amplitude of the first edgewise mode has been obtained based on the aforementioned FSI simulations with varying structural damping. The power dissipated by the modal structural damping can be computed from the integral of force times velocity, similarly to the power transferred by the aerodynamic forces in Eq. (2). The force due to the modal damping is proportional to the velocity, which means that the dissipated power (force times velocity) is proportional to the velocity squared. The velocity is proportional to the amplitude; thus, the dissipated power due to structural damping is proportional to the amplitude squared. The proportionality constant FSTRUC is here determined from the FSI simulations, where the power dissipated by damping is found using the equilibrium of structurally dissipated power and the aerodynamic power injection when the amplitude of the vibration has converged to a steady state.

To test various amplitudes using the coupling, the aforementioned simulations with varying damping coefficients are used. This means that the power dissipated is also scaled accordingly with the same scale δDAMP. Thereby, the fit of a factor FSTRUC can be based on multiple cases. The power dissipated by structural damping PSTRUC can be described as

(4) P STRUC = F STRUC δ DAMP A 2 ,

where A is the tip displacement amplitude in edgewise direction at steady-state vibrations. With the various cases of δDAMP being [0, 0.5, 1, 2, 5, 10, 15], a factor of FSTRUC=540 W m−2 is found suitable for all cases as depicted in Fig. 11.

This simple model only estimates the structural damping dissipation when vibrating in the first edgewise mode as this is what it is fitted for here. Another mode would yield another value of FDAMP as the motion differs. This model can be used to assess the dissipated power from structural damping for various amplitudes, making it possible to estimate an equilibrium amplitude by use of forced-motion simulations with prescribed amplitudes.

Figure 11Simple damping model of dissipated power compared with aerodynamic power obtained through FSI simulations.


3.1.1 Amplitude estimations from forced-motion simulations

Having obtained a simple relationship between dissipated power and amplitude from the structural damping, the FM simulations can be used to not only assess the amplitude where the aerodynamics stop injecting power, but also assess the limit cycle amplitude when including the structural damping. This amplitude can be found at the intersection between the curve of the structural damping and the aerodynamic power. An example is given in Fig. 12 for flow case P100-I50 with the FM aero-power curve along with the curve of the damping model along with the result of the FSI simulation. As seen, the aerodynamic power and motion amplitude match well between the coupled FSI simulation and the cross section between the forced-motion curve and the damping model. For FSI, the equilibrium state is found at  1.18 m amplitude with an aerodynamic power of  775 W. The equilibrium for the forced-motion case and the damping model is found at  1.24 m with a corresponding aerodynamic power of  880 W. This difference in aerodynamic power injection might seem high, but as indicated before, this intersection happens at a very steep range of the power injection curve, which is also seen by the only 6 cm difference between the FSI resulting amplitude and the FM estimated amplitude. Looking at the difference between aerodynamic power and the power dissipated by structural damping, the largest effective input of power happens around 0.8 m amplitude. This indicates that the development of VIV amplitudes will happen the fastest here while the amplitude development over time will slow down as the LCO establishes.

Figure 12Forced-motion results for various amplitudes along with FSI result and damping model for baseline structural damping. Case P100-I50.


3.2 Multiple-amplitude cases and initial condition dependencies

Originally in the work of Horcas et al. (2022), cases P90-I50 and P95-I50 were found to not cause VIVs, as found in the neighboring case P100-I50. However, from the FM amplitude sweeps depicted in Fig. 6, it is seen that a region of positive power injection occurs, when forced amplitudes are higher than  0.1–0.4 m depending on flow case. This result is also found through FSI simulations, if initial motions are triggered through an external forcing/motion. Examples of this could be buffeting motion from turbulent gusts or motion from inertial forces when moving the blade through a yaw, pitch or azimuthal change. In this work, various ways of triggering the VIVs branches have been tested, these being turbulent buffeting inflow, a wind speed ramp and finally externally prescribed synthetic initial forcing of the mode through the FSI coupling. All these methods are found to successfully trigger the VIVs as the amplitude of vibration merely needs to reach a certain value before the power injection takes care of building up the VIVs. The latter synthetic trigger method allows most control and is therefore chosen in this study. Initially forcing the blade harmonically in the first edgewise frequency until an amplitude point is reached, or merely pulling the blade in the edgewise direction by a static loading and releasing it, triggers VIVs as depicted in Fig. 13. The figure shows the low-equilibrium (LE) LCO being triggered by a harmonic forcing in the first 50 s and the high-equilibrium (HE) LCO being triggered by pull and release. After the period of synthetic forcing, the coupling continues as normal with the CFD forcing, and the structural response calculations until a steady-state response is found. The forced triggering of VIVs also allows a speed-up of two-way-coupled simulations, if the final LCO amplitude is the desired result, since the CFD does not need to be active in the initial forcing.

Figure 13High-equilibrium (HE) and low-equilibrium (LE) LCOs triggered by synthetic initial conditions for flow case P90-I50. HE is triggered by a pull and release, while LE is triggered by harmonic forcing of the blade for the first 50 s.


Considering again the power curves of the three studied flow cases but this time subtracting the power dissipated by the structural damping, a stability figure can be obtained as presented in Fig. 14. Here, the effective power injection is found as the difference between aerodynamic power injection and the structural damping dissipation (Eq. 4). As can be seen, the FM simulations and the FSI simulations do not agree entirely at low amplitudes, as the FM estimates slightly positive effective power injection for case P95-I50 in the whole range of amplitudes from  0.1 up to  1.4 m. As also depicted in Fig. 15, the FSI finds three equilibrium LCOs for this case, depending on the initial conditions. The amplitudes of motion corresponding to these equilibrium points are  0.1 m (as also found in Horcas et al.2022), 0.33 and 1.4 m. The LCOs, despite having quite different amplitudes, all vibrate harmonically in the first edgewise mode with  0.67 Hz. Reasons why the FM and FSI not to agree perfectly here could be the simple estimate of damping made by a fit or the assumptions done in FM simulations, for instance, neglecting torsion or using the static deflection from P100-I50 as a baseline. The simplified damping model fit is likely the reason for the FSI simulations at high amplitudes for P90-I50 and P95-I50 not being at precisely 0 W power despite being at a stable LCO level. For the P90-I50 flow case, the conclusions are similar, but with better agreement between intersection points and FSI results for low amplitudes.

Figure 14Effective power PEFF=PAERO-PSTRUC for cases P100, P95 and P90, all with inclination I50. Magenta and green areas, respectively, indicate unstable and stable regions. PSTRUC is here defined as Eq. (4), which is why FSI results are not exactly at 0 W despite being at stable LCO.


Figure 15Edgewise tip displacement of flow case P95-I50, with three equilibrium states.


3.3 Time series generation of vibrations

By knowing the relation between aerodynamically injected power, structurally dissipated power and the amplitude of the mode, one can generate a time series of the development of vibration. This can be useful in terms of assessing whether a vibration is worrisome or not, as a long development time will need long periods of steady inflow conditions, while quickly developing vibrations can happen over a short period, making the risk here more critical.

Considering only small relative deflections in the first edgewise mode, the harmonic motion of the damped linear one-degree-of-freedom system in free vibration can be described as

(5) x ( t ) = A 0 e - ζ ω t sin ( ω ̃ t + ψ ) = A ( t ) sin ( ω ̃ t + ψ ) ,

with the initial amplitude A0, the damping ratio ζ, and the phase ψ, which is assumed to be zero in the following. Further, for the small damping of the first edgewise mode, the damped frequency is almost equal to the undamped natural frequency ω1e=ω̃ω. The time derivative of the amplitude can be found to be a function of the amplitude itself, the damping ratio and the natural frequency:

(6) d A d t = - ζ ω A 0 e - ζ ω t = - ζ ω A ( t ) .

For a damped free vibration the only power dissipation is caused by the structural damping, with the power proportional to the damping ratio ζ. For the vibrations considered in this article, the aerodynamic forces are another source of power injection, resulting in the effective power:

(7) P EFF ( A ) = P AERO ( A ) - P STRUC ( A ) = P AERO ( A ) - F STRUC A 2 ,

with the assumption again that the frequency and mode shape are approximately those of the first edgewise mode. From this assumption follows that the proportionality constant between the structural power and the structural damping ratio has the same magnitude as the proportionality constant between the effective power and an effective damping ratio ζeff:


With this effective damping ratio, the change in tip amplitude at a specific time step can be found as

(10) d A d t = - ζ EFF ( A ( t ) ) ω 1 e A ( t ) = P EFF ( A ) ζ STRUC P STRUC ( A ) ω 1 e A ( t ) .

With these assumptions, a simple backwards differencing scheme can be used to find the amplitude variation over time using interpolated values of aerodynamic power found in FM simulations: the amplitude is initialized with a certain value A0; the time derivative of the amplitude is determined from Eq. (10), with the effective power as shown in Fig. 14; and the amplitude in the following time step is found as A(t+Δt)=A(t)+(dA/dt)Δt. The amplitude time series A(t) multiplied with sin (ω1et) will generate a time signal of the vibration development. The fact that the aerodynamic power is obtained from interpolation between a finite number of FM simulations is naturally a limitation of the method, as the amplitude discretization will determine the accuracy of the method.

Figure 16Time series of tip displacement in edgewise direction for FSI simulation, and generated signal based on FSI amplitude at 20 s as initial condition. Flow case P100-I50.


Figure 17Generated time series for P90-I50 along with corresponding FSI simulations, for low LCO amplitude (a) and high LCO amplitude (b). In FSI, the high-amplitude LCO was obtained by pulling the blade to 1.75 m edgewise amplitude and releasing it to vibrate in the flow.


Figure 16 shows the generated signal along with the obtained FSI signal for flow case P100-I50. The time evolution of the vibration amplitude is in good agreement between the two signals, despite the discretization of the aerodynamic power in the given range being based on only seven FM simulations.

The two-amplitude case of P90-I50 can also be generated, since FM simulations here also predict two LCO amplitudes, as also depicted in Fig. 17. This is done by defining an initial amplitude below or above the stable region depicted in Fig. 14. Since the time generation method is based on averaged power per amplitude, from FM simulations with steady inflow, the generation will not be able to capture triggering of amplitude branches by, for instance, turbulence buffeting.

For P95-I50, the FM simulations predict a slightly higher aerodynamic power than damping for the low amplitudes, meaning that here the generation will predict only one LCO stabilized amplitude, whereas the FSI simulation as mentioned predicted three. This might not be a critical point, as the low-amplitude LCOs are of less importance when designing the wind turbine blade.

4 Concluding remarks

The main objective of the current study has been to conduct and compare simulations of vortex-induced vibrations (VIVs) by the use of two methodologies both based on CFD. Simulations fully coupled to a structural solver are compared with simulations with the blade motion imposed by an analytical expression of the first edgewise mode instead. Fully coupled simulations are beneficial as no assumptions are made of the structural response but limited by the fact that long time series are needed in order to reach a limit cycle oscillation (LCO). Forced-motion simulations, on the other hand, have a very limited initial transient as the LCO is prescribed, and only the flow needs to converge to a periodic state. However, these simulations need assumptions on the motion shape and frequency, and various amplitudes need to be simulated.

For the studied cases of high inclination angles and pitch angles close to perpendicular to the chord, a very good agreement is seen between coupled and forced-motion simulations. This is found by comparing power injection along the span of the blade, which is a parameter sensitive to both loads, motion and flow shedding frequencies, making it difficult to predict.

By fitting a simple model for the blade structural damping for the first edgewise mode shape, it is possible not only to make good assessments of the final LCO amplitude, but also estimate the development of the motion by using data from multiple forced-motion simulations with varying amplitude.

Finally, the study reveals a VIV dependency on the initial conditions of the blade motion. Cases which at steady inflow and without initial motion do not build up VIVs might build up large VIVs if perturbed by initial motion or unsteady loading. Examples of such perturbations could be turbulent inflow or actions from the turbine controller such as yawing or pitching. An example of this was demonstrated for the flow case P95-I50, where three different stable LCO amplitudes can be found if the initial amplitudes are varied in the FSI simulations.

5 Future studies

Vortex-induced vibrations for wind turbines are still an understudied field, with many possible discoveries to be made. Within single-blade VIVs, the studied parameter space can be expanded by looking at different modes, a broader range of inflow angles and velocities, inflow turbulence, and much more. However, expanding the parameter space rapidly increases computational efforts needed, which is why the forced-motion approach, examined in the current study, might become relevant in order to reduce the cost.

Looking further, investigations of full rotors and entire wind turbine setups are naturally relevant. This would allow different modes to be excited for the same wind speed range, again expanding the parameter space. Here, the blades will interact in terms of injection and dissipation of power, depending on their positions, and the flow might be altered by upstream blades for some positions.

Appendix A: Abbreviations
BEM Blade element momentum
BL Baseline damping
CFD Computational fluid dynamics
DTU Technical University of Denmark
EUDP Energy Technology Development and
Demonstration Program
FM Forced motion
FSI Fluid–structure interaction
IDDES Improved Delayed Detached Eddy Simulation
IEA International Energy Agency
LCO Limit cycle oscillation
LES Large eddy simulation
MBD Multi-body dynamics
MPI Message Passing Interface
NREL National Renewable Energy Laboratory
RANS Reynolds-averaged Navier–Stokes
RWT Reference wind turbine
VIV Vortex-induced vibration
ZD Zero damping
Data availability

The geometry and the structural model of the considered wind turbine are publicly available. However, the structural and fluid solvers are licensed. The data that support the findings of this study are available upon reasonable request to the corresponding author.

Author contributions

CG conducted FSI simulations, performed the analysis of FSI and FM simulations, and primarily wrote the paper. NNS conducted FM simulations and supported FSI analysis. GRP supported the analysis of FSI simulations and developed the time series generation method after an initial idea from SGH. SGH performed the initial FSI study and meshing and supported the FSI analysis. All authors contributed in editing the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work is part of the PRESTIGE project, and computational resources were provided by the DTU Risø cluster Sophia (DTU Computing Center2021).

The authors would like to thank Pim Jacobs, Aqeel Ahmed, Bastien Duboc and the rest of the team from Siemens Gamesa Renewable Energy for many fruitful discussions regarding the presented work in the PRESTIGE project.

Financial support

This research has been supported by Innovation Fund Denmark (grant no. 9090-00025B).

Review statement

This paper was edited by Roland Schmehl and reviewed by Vasilis A. Riziotis and one anonymous referee.


Bak, C., Zahle, F., Bitsche, R., Kim, T., Yde, A., Henriksen, L., Hansen, M., Blasques, J., Gaunaa, M., and Natarajan, A.: Description of the DTU 10 MW Reference Wind Turbine, technical report, DTU Wind Energy Report_I_0092, DTU Wind Energy, 2013. a

Barlas, T., Pirrung, G. R., Ramos-García, N., Horcas, S. G., Mikkelsen, R. F., Olsen, A. S., and Gaunaa, M.: Wind tunnel testing of a swept tip shape and comparison with multi-fidelity aerodynamic simulations, Wind Energ. Sci., 6, 1311–1324,, 2021. a

Bearman, P. W.: Circular cylinder wakes and vortex-induced vibrations, J. Fluid. Struct., 27, 648–658,, 2011. a

Bechmann, A., Sørensen, N., and Zahle, F.: CFD simulations of the MEXICO rotor, Wind Energy, 14, 677–689, 2011. a

Bertagnolio, F., Gaunaa, M., Hansen, M., Sørensen, N., and Rasmussen, F.: Computation of aerodynamic damping for wind turbine applications, in: CD-Rom proceedings, edited by: Tsahalis, D., Greek Association of Computational Mechanics, 4th GRACM Congress on Computational Mechanics, Patras, Greece, 27–29 June 2002, (last access: 1 November 2022), 2002. a

Blevins, R.: Flow-induced vibration, Krieger Publishing Company, Malabar, Florida, USA, 2nd edn., ISBN 1-57524-183-8, 2001. a

Blevins, R. D. and Coughran, C. S.: Experimental investigation of vortex-induced vibration in one and two dimensions with variable mass, damping, and reynolds number, J. Fluid Eng., 131, 1012021–1012027,, 2009. a

Bortolotti, P., Canet Tarrés, H., Dykes, K., Merz, K., Sethuraman, L., Verelst, D., and Zahle, F.: Systems Engineering in Wind Energy – WP2.1 Reference Wind Turbines, Tech. rep., National Renewable Energy Laboratory (NREL), (last access: 1 November 2022), 2019. a, b, c, d

Ceyhan, J. S. O., Boorsma, K., Gonzalez, A., Munduate, X., Pires, O., Sørensen, N., Ferreira, C., Sieros, G., Madsen, J., Voutsinas, S., Lutz, T., Barakos, G., Colonia, S., Heißelmann, H., Meng, F., and Croce, A.: Latest results from the EU project AVATAR: Aerodynamic modelling of 10 MW wind turbines, J. Phys.-Conf. Ser., 753, 022017,, 2016. a

DTU Computing Center: DTU Computing Center resources, Technical University of Denmark,, 2021. a

EUDP: Predicting wind turbine stability in general inflow PRESTIGE, (last access: 1 November 2022), 2020. a

García Ramos, N., Sessarego, M., and Horcas, S. G.: Aero-hydro-servo-elastic coupling of a multi-body finite-element solver and a multi-fidelity vortex method, Wind Energy, 24, 481–501,, 2020. a

Grinderslev, C., Belloni, F., Horcas, S. G., and Sørensen, N. N.: Investigations of aerodynamic drag forces during structural blade testing using high-fidelity fluid–structure interaction, Wind Energ. Sci., 5, 543–560,, 2020a. a

Grinderslev, C., Vijayakumar, G., Ananthan, S., Sørensen, N., Zahle, F., and Sprague, M.: Validation of blade-resolved computational fluid dynamics for a MW-scale turbine rotor in atmospheric flow, J. Phys.-Conf. Ser., 1618, 052049,, 2020b. a

Grinderslev, C., Horcas, S., and Sørensen, N.: Fluid–structure interaction simulations of a wind turbine rotor in complex flows, validated through field experiments, Wind Energy, 24, 1426–1442,, 2021a. a

Grinderslev, C., Sørensen, N. N., Horcas, S. G., Troldborg, N., and Zahle, F.: Wind turbines in atmospheric flow: fluid–structure interaction simulations with hybrid turbulence modeling, Wind Energ. Sci., 6, 627–643,, 2021b. a

Gritskevich, M. S., Garbaruk, A. V., Schütze, J., and Menter, F. R.: Development of DDES and IDDES formulations for the kω shear stress transport model, Flow Turbul. Combust., 88, 431–449,, 2012. a

Hansen, M. H.: Aeroelastic stability analysis of wind turbines using an eigenvalue approach, Wind Energy, 7, 133–143,, 2004. a

Heinz, J., Sørensen, N., and Zahle, F.: Fluid-structure interaction computations for geometrically resolved rotor simulations using CFD, Wind Energy, 19, 2205–2221, 2016a. a

Heinz, J., Sørensen, N., Zahle, F., and Skrzypinski, W.: Vortex-induced vibrations on a modern wind turbine blade, Wind Energy, 19, 2041–2051,, 2016b. a, b

Horcas, S., Madsen, M., Sørensen, N., and Zahle, F.: Suppressing Vortex Induced Vibrations of Wind Turbine Blades with Flaps, in: Recent Advances in CFD for Wind and Tidal Offshore Turbines, edited by: Ferrer, E. and Montlaur, A., Springer, 11–24,, 2019. a, b, c

Horcas, S., Barlas, T., Zahle, F., and Sørensen, N.: Vortex induced vibrations of wind turbine blades: Influence of the tip geometry, Phys. Fluids, 32, 065104,, 2020. a, b, c, d

Horcas, S. G., Sørensen, N. N., Zahle, F., Pirrung, G. R., and Barlas, T.: Vibrations of wind turbine blades in standstill: Mapping the influence of the inflow angles, Phys. Fluids, 34, 033306,, 2022. a, b, c, d, e, f, g

Hu, P., Sun, C., Zhu, X., and Du, Z.: Investigations on vortex-induced vibration of a wind turbine airfoil at a high angle of attack via modal analysis, J. Renew. Sustain. Ener., 13, 054105,, 2021. a

Koopmann, G. H.: The vortex wakes of vibrating cylinders at low Reynolds numbers, J. Fluid Mech., 28, 501–512, 1967. a

Larsen, T. and Hansen, A.: How 2 HAWC2, the user's manual, Tech. rep., Risø National Laboratory, (last access: 1 November 2022), ISBN 978-87-550-3583-6, 2007. a

Lee, J. H. and Bernitsas, M. M.: High-damping, high-Reynolds VIV tests for energy harnessing using the VIVACE converter, Ocean Eng., 38, 1697–1712,, 2011. a

Lupi, F., Höffer, R., and Niemann, H. J.: Aerodynamic damping in vortex resonance from aeroelastic wind tunnel tests on a stack, J. Wind Eng. Ind. Aerod., 208, 104438,, 2021. a

Madsen, H. A., Larsen, T. J., Pirrung, G. R., Li, A., and Zahle, F.: Implementation of the blade element momentum model on a polar grid and its aeroelastic load impact, Wind Energ. Sci., 5, 1–27,, 2020. 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,, 2018. a

Meskell, C. and Pellegrino, A.: Vortex shedding lock-in due to pitching oscillation of a wind turbine blade section at high angles of attack, Int. J. Aerospace Eng., 2019, 6919505,, 2019. a

Michelsen, J.: Basis3D – A Platform for Development of Multiblock PDE Solvers., Tech. rep., Risø National Laboratory, (last access: 1 November 2022), 1992. a

Michelsen, J.: Block Structured Multigrid Solution of 2D and 3D Elliptic PDE's., Tech. rep., Technical University of Denmark, 1994. a

Pavese, C., Wang, Q., Kim, T., Jonkman, J., and Sprague, M.: HAWC2 and BeamDyn: Comparison Between Beam Structural Models for Aero-Servo-Elastic Frameworks, in: Proceedings of the EWEA Annual Event and Exhibition 2015, European Wind Energy Association (EWEA), paper for poster presentation, EWEA Annual Conference and Exhibition 2015, Paris, France, 17–20 November 2015, NREL/CP-5000-65115, (last access: 1 November 2022), 2015. a

Placzek, A., Sigrist, J. F., and Hamdouni, A.: Numerical simulation of an oscillating cylinder in a cross-flow at low Reynolds number: Forced and free oscillations, Comput. Fluids, 38, 80–100,, 2009. a

Riva, R., Horcas, S. G., Sørensen, N. N., Grinderslev, C., and Pirrung, G. R.: Stability analysis of vortex-induced vibrations on wind turbines, J. Phys.-Conf. Ser., 2265, 042054,, 2022. a

Schepers, J. G., Boorsma, K., Madsen, H. A., Pirrung, G. R., Bangga, G., Guma, G., Lutz, T., Potentier, T., Braud, C., Guilmineau, E., Croce, A., Cacciola, S., Schaffarczyk, A. P., Lobo, B. A., Ivanell, S., Asmuth, H., Bertagnolio, F., Sørensen, N. N., Shen, W. Z., Grinderslev, C., Forsting, A. M., Blondel, F., Bozonnet, P., Boisard, R., Yassin, K., Hoening, L., Stoevesandt, B., Imiela, M., Greco, L., Testa, C., Magionesi, F., Vijayakumar, G., Ananthan, S., Sprague, M. A., Branlard, E., Jonkman, J., Carrion, M., Parkinson, S., and Cicirello, E.: IEA Wind TCP Task 29, Phase IV: Detailed Aerodynamics of Wind Turbines, Zenodo,, 2021. a

Shur, M. L., Spalart, P. R., Strelets, M. K., and Travin, A. K.: A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities, Int. J. Heat Fluid Fl., 29, 1638–1649,, 2008.  a

Sørensen, N. and Schreck, S.: Transitional DDES computations of the NREL Phase-VI rotor in axial flow conditions, J. Phys.-Conf. Ser., 555, 012096,, 2014. a

Sørensen, N., Zahle, F., Boorsma, K., and Schepers, G.: CFD computations of the second round of MEXICO rotor measurements, J. Phys.-Conf. Ser., 753, 022054,, 2016. a

Sørensen, N.: General purpose flow solver applied to flow over hills, PhD thesis, (last access: 1 November 2022) 1995. a

Sørensen, N.: HypGrid2D. A 2-d mesh generator, Tech. rep., Risø National Laboratory, (last access: 1 November 2022), ISBN 87-550-2368-1, 1998. a

Sørensen, N. N., Michelsen, J. A., and Schreck, S.: Navier–Stokes predictions of the NREL phase VI rotor in the NASA Ames 80 ft × 120 ft wind tunnel, Wind Energy, 5, 151–169,, 2002. a

Williamson, C. and Govardhan, R.: Vortex-Induced Vibrations, Annu. Rev. Fluid Mech., 36, 413–55,, 2004. a

Zhang, M., Song, Y., Wang, J., Abdelkefi, A., Yu, H., and Yu, H.: Vortex-induced vibration of a circular cylinder with nonlinear stiffness : prediction using forced vibration data, Nonlinear Dynam., 108, 1867–1884,, 2022. a

Short summary
As wind turbines increase in size, the risk of flow-induced instabilities increases. This study investigates the phenomenon of vortex-induced vibrations (VIVs) on a large 10 MW wind turbine blade using two high-fidelity methods. It is found that VIVs can occur with multiple equilibrium states for the same flow case, showing an dependence on the initial conditions. This means that a blade which is stable in a flow can become unstable if, e.g., a turbine operation provokes an initial vibration.