the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Multiple limit cycle amplitudes in highfidelity predictions of standstill wind turbine blade vibrations
Niels Nørmark Sørensen
Georg Raimund Pirrung
Sergio González Horcas
In this study, vortexinduced 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 forcedmotion 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 nonlinearity 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 forcedmotion simulations, even for the multiple stable vibration cases, with good agreement with fully coupled simulations. Finally, a time series generation from forcedmotion 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 buildup time series.
Please read the corrigendum first before continuing.

Notice on corrigendum
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Article
(7578 KB)

The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(7578 KB)  Corrigendum
 BibTeX
 EndNote
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 vortexinduced 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 socalled lockin 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 socalled 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=\mathrm{0}{}^{\circ}$ when the horizontal wind projection strikes directly the leading edge of the untwisted blade, $P=\mathrm{90}{}^{\circ}$ when striking the pressure side perpendicularly, and $P=\mathrm{180}{}^{\circ}$ 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 x–z 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 highamplitude 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.
The main objective of the present study is the exploration of the forcedmotion approach as a complementary analysis tool for the study of wind turbine blade VIVs. The forcedmotion technique has been extensively applied to academic configurations such as cylinders (Williamson and Govardhan, 2004). 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 lockin region becomes narrower when simulating low amplitudes of motion, both for cylinders (Koopmann, 1967; Blevins, 2001; Placzek et al., 2009) and airfoils (Meskell and Pellegrino, 2019; Hu et al., 2021). Additionally, the forced motion disregards the structural damping, which can have an important effect on the VIV phenomenon (Lee and Bernitsas, 2011; Blevins and Coughran, 2009). However, if assumptions on the loading and on the response can be made, the forcedmotion strategy is able to estimate the amplitude of vibrations of a given inflow (Bearman, 2011), which is often the targeted unknown in a design process. This is usually done by performing a series of forcedmotion simulations of varying amplitude and subsequently building a reducedorder 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 forcedmotion approach is the fact that neither timemarching 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 buildup is needed.
The current study is part of the PRESTIGE project (EUDP, 2020), looking into stability of wind turbines in various aspects. In this particular work, the forcedmotion 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.1 Flow solver
To calculate the fluid flow, the computational fluid dynamics (CFD) code EllipSys3D (Michelsen, 1994, 1992; Sørensen, 1995) is used. This code solves the incompressible Navier–Stokes equations through Reynoldsaveraged 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 sublayer in LES. The code utilizes structured grids in curvilinear coordinates using the finitevolume method. This enables highly scalable computations through MPI with multiblock decomposition and grid sequencing. Convection is solved through various schemes such as central difference, secondorder upwind, QUICK or combinations of said schemes. Here, a combination of the QUICK and a fourthorder 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 Schreck, 2014; Sørensen et al., 2002), and recent comparisons with wind tunnel tests of a curved blade tip (Barlas et al., 2021) and with the fullscale 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 bodyfitted 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ørensen, 1998). 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$\times {\mathrm{10}}^{\mathrm{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 nearwall volume mesh.
2.2 Coupling framework
To capture the interaction between fluid and structure response, two methodologies are investigated: forcedmotion (FM) and fully coupled fluid–structure interaction (FSI) simulations. These methodologies are also commonly known as oneway and twoway coupling, respectively. For both simulation types, it is here for simplicity assumed that the blade is clamped at the root. This means that singeblade modes can be captured, but full rotor/turbine modes are omitted.
2.2.1 Structural solver
To solve the structural response for twoway coupling HAWC2 (Larsen and Hansen, 2007; Madsen et al., 2020; Pavese et al., 2015) is used. HAWC2 is a multiphysics code widely used as an aeroelastic solver using multibody 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 finiteelement Timoshenko beams, which are connected through constraint equations to enable considerations of nonlinearities. 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 socalled 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 twowaycoupled 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 Forcedmotion simulations
Forcedmotion (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 seventhorder polynomial to the mode shape obtained from the aeroelastic stability tool HAWCSTAB2 (Hansen, 2004). 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 seventhorder 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.
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=\mathrm{100}{}^{\circ}$, $I=\mathrm{50}{}^{\circ}$ case. Two cases at the vicinity of this case, P=95 and $P=\mathrm{90}{}^{\circ}$, both with $I=\mathrm{50}{}^{\circ}$, which both show no vibration in the aforementioned mapping are also included. These three cases will be identified as P90I50, P95I50 and P100I50 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 P_{AERO} is found by the timeaveraged cellwise integration of the power along each of the cells of the blade surface mesh.
where T is the time for one motion period, i is the cell index, N is the total number of cells, and F_{i} and ${\dot{u}}_{i}$ are the force and velocity at every cell.
The three aforementioned flow scenarios have been tested for various forcedmotion 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 P100I50, 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 lowamplitude 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.
Case P100I50 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 forcedmotion 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.
A clear benefit of the forcedmotion methodology is that the limit cycle amplitude is set in advance, meaning that no vibration buildup 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.
Taking a close look at the average aerodynamic power injection along the blade span for the coupled and forcedmotion cases, Fig. 9 shows the results for the flow case P100I50. The spanwise power per cycle P_{cycle} is here found as the injected work averaged over the time span of n cycles as presented in Horcas et al. (2022). The work w_{s} over n cycles is found as the integration over the time range of the sectional loads F_{s} and the velocity of the considered section ${\dot{u}}_{\mathrm{s}}$.
where t_{0} is chosen to be after the establishment of the LCO and the integration is performed over n periods of motion:
For the FSI simulation results, the baseline damping (BL) and a zero structural damping (ZD) case were chosen and compared to the forcedmotion 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.
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 P100I50 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 steadystate 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.
A simple model for the power dissipated by structural damping P_{STRUC} 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 F_{STRUC} 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 F_{STRUC} can be based on multiple cases. The power dissipated by structural damping P_{STRUC} can be described as
where A is the tip displacement amplitude in edgewise direction at steadystate vibrations. With the various cases of δ_{DAMP} being [0, 0.5, 1, 2, 5, 10, 15], a factor of F_{STRUC}=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 F_{DAMP} 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 forcedmotion simulations with prescribed amplitudes.
3.1.1 Amplitude estimations from forcedmotion 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 P100I50 with the FM aeropower 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 forcedmotion 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 forcedmotion 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.
3.2 Multipleamplitude cases and initial condition dependencies
Originally in the work of Horcas et al. (2022), cases P90I50 and P95I50 were found to not cause VIVs, as found in the neighboring case P100I50. 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 lowequilibrium (LE) LCO being triggered by a harmonic forcing in the first 50 s and the highequilibrium (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 steadystate response is found. The forced triggering of VIVs also allows a speedup of twowaycoupled simulations, if the final LCO amplitude is the desired result, since the CFD does not need to be active in the initial forcing.
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 P95I50 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 P100I50 as a baseline. The simplified damping model fit is likely the reason for the FSI simulations at high amplitudes for P90I50 and P95I50 not being at precisely 0 W power despite being at a stable LCO level. For the P90I50 flow case, the conclusions are similar, but with better agreement between intersection points and FSI results for low amplitudes.
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 onedegreeoffreedom system in free vibration can be described as
with the initial amplitude A_{0}, 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} =$\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u0303}}{\mathit{\omega}}\approx \mathit{\omega}$. The time derivative of the amplitude can be found to be a function of the amplitude itself, the damping ratio and the natural frequency:
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:
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
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 A_{0}; 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+\mathrm{\Delta}t)=A\left(t\right)+(\mathrm{d}A/\mathrm{d}t)\mathrm{\Delta}t$. The amplitude time series A(t) multiplied with sin (ω_{1e}t) 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 16 shows the generated signal along with the obtained FSI signal for flow case P100I50. 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 twoamplitude case of P90I50 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 P95I50, 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 lowamplitude LCOs are of less importance when designing the wind turbine blade.
The main objective of the current study has been to conduct and compare simulations of vortexinduced 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). Forcedmotion 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 forcedmotion 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 forcedmotion 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 P95I50, where three different stable LCO amplitudes can be found if the initial amplitudes are varied in the FSI simulations.
Vortexinduced vibrations for wind turbines are still an understudied field, with many possible discoveries to be made. Within singleblade 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 forcedmotion 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.
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  Multibody dynamics 
MPI  Message Passing Interface 
NREL  National Renewable Energy Laboratory 
RANS  Reynoldsaveraged Navier–Stokes 
RWT  Reference wind turbine 
VIV  Vortexinduced vibration 
ZD  Zero damping 
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.
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.
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 Center, 2021).
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.
This research has been supported by Innovation Fund Denmark (grant no. 909000025B).
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., RamosGarcí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 multifidelity aerodynamic simulations, Wind Energ. Sci., 6, 1311–1324, https://doi.org/10.5194/wes613112021, 2021. a
Bearman, P. W.: Circular cylinder wakes and vortexinduced vibrations, J. Fluid. Struct., 27, 648–658, https://doi.org/10.1016/j.jfluidstructs.2011.03.021, 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: CDRom proceedings, edited by: Tsahalis, D., Greek Association of Computational Mechanics, 4th GRACM Congress on Computational Mechanics, Patras, Greece, 27–29 June 2002, https://backend.orbit.dtu.dk/ws/portalfiles/portal/148414970/paper.pdf (last access: 1 November 2022), 2002. a
Blevins, R.: Flowinduced vibration, Krieger Publishing Company, Malabar, Florida, USA, 2nd edn., ISBN 1575241838, 2001. a
Blevins, R. D. and Coughran, C. S.: Experimental investigation of vortexinduced vibration in one and two dimensions with variable mass, damping, and reynolds number, J. Fluid Eng., 131, 1012021–1012027, https://doi.org/10.1115/1.3222904, 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), https://www.osti.gov/biblio/1529216ieawindtcptasksystemsengineeringwindenergywp2referencewindturbines (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, https://doi.org/10.1088/17426596/753/2/022017, 2016. a
DTU Computing Center: DTU Computing Center resources, Technical University of Denmark, https://doi.org/10.48714/DTU.HPC.0001, 2021. a
EUDP: Predicting wind turbine stability in general inflow PRESTIGE, https://energiforskning.dk/en/node/15974 (last access: 1 November 2022), 2020. a
García Ramos, N., Sessarego, M., and Horcas, S. G.: Aerohydroservoelastic coupling of a multibody finiteelement solver and a multifidelity vortex method, Wind Energy, 24, 481–501, https://doi.org/10.1002/we.2584, 2020. a
Grinderslev, C., Belloni, F., Horcas, S. G., and Sørensen, N. N.: Investigations of aerodynamic drag forces during structural blade testing using highfidelity fluid–structure interaction, Wind Energ. Sci., 5, 543–560, https://doi.org/10.5194/wes55432020, 2020a. a
Grinderslev, C., Vijayakumar, G., Ananthan, S., Sørensen, N., Zahle, F., and Sprague, M.: Validation of bladeresolved computational fluid dynamics for a MWscale turbine rotor in atmospheric flow, J. Phys.Conf. Ser., 1618, 052049, https://doi.org/10.1088/17426596/1618/5/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, https://doi.org/10.1002/we.2639, 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, https://doi.org/10.5194/wes66272021, 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, https://doi.org/10.1007/s1049401193784, 2012. a
Hansen, M. H.: Aeroelastic stability analysis of wind turbines using an eigenvalue approach, Wind Energy, 7, 133–143, https://doi.org/10.1002/we.116, 2004. a
Heinz, J., Sørensen, N., and Zahle, F.: Fluidstructure 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.: Vortexinduced vibrations on a modern wind turbine blade, Wind Energy, 19, 2041–2051, https://doi.org/10.1002/we.1967, 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, https://doi.org/10.1007/9783030118877_2, 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, https://doi.org/10.1063/5.0004005, 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, https://doi.org/10.1063/5.0088036, 2022. a, b, c, d, e, f, g
Hu, P., Sun, C., Zhu, X., and Du, Z.: Investigations on vortexinduced vibration of a wind turbine airfoil at a high angle of attack via modal analysis, J. Renew. Sustain. Ener., 13, 054105, https://doi.org/10.1063/5.0040509, 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, https://orbit.dtu.dk/en/publications/how2hawc2theusersmanual (last access: 1 November 2022), ISBN 9788755035836, 2007. a
Lee, J. H. and Bernitsas, M. M.: Highdamping, highReynolds VIV tests for energy harnessing using the VIVACE converter, Ocean Eng., 38, 1697–1712, https://doi.org/10.1016/j.oceaneng.2011.06.007, 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, https://doi.org/10.1016/j.jweia.2020.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, https://doi.org/10.5194/wes512020, 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, https://doi.org/10.1088/17426596/1037/2/022011, 2018. a
Meskell, C. and Pellegrino, A.: Vortex shedding lockin due to pitching oscillation of a wind turbine blade section at high angles of attack, Int. J. Aerospace Eng., 2019, 6919505, https://doi.org/10.1155/2019/6919505, 2019. a
Michelsen, J.: Basis3D – A Platform for Development of Multiblock PDE Solvers., Tech. rep., Risø National Laboratory, https://backend.orbit.dtu.dk/ws/portalfiles/portal/272917945/Michelsen_J_Basis3D.pdf (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 AeroServoElastic 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/CP500065115, https://www.nrel.gov/docs/fy16osti/65115.pdf (last access: 1 November 2022), 2015. a
Placzek, A., Sigrist, J. F., and Hamdouni, A.: Numerical simulation of an oscillating cylinder in a crossflow at low Reynolds number: Forced and free oscillations, Comput. Fluids, 38, 80–100, https://doi.org/10.1016/j.compfluid.2008.01.007, 2009. a
Riva, R., Horcas, S. G., Sørensen, N. N., Grinderslev, C., and Pirrung, G. R.: Stability analysis of vortexinduced vibrations on wind turbines, J. Phys.Conf. Ser., 2265, 042054, https://doi.org/10.1088/17426596/2265/4/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, https://doi.org/10.5281/zenodo.4813068, 2021. a
Shur, M. L., Spalart, P. R., Strelets, M. K., and Travin, A. K.: A hybrid RANSLES approach with delayedDES and wallmodelled LES capabilities, Int. J. Heat Fluid Fl., 29, 1638–1649, https://doi.org/10.1016/j.ijheatfluidflow.2008.07.001, 2008. a
Sørensen, N. and Schreck, S.: Transitional DDES computations of the NREL PhaseVI rotor in axial flow conditions, J. Phys.Conf. Ser., 555, 012096, https://doi.org/10.1088/17426596/555/1/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, https://doi.org/10.1088/17426596/753/2/022054, 2016. a
Sørensen, N.: General purpose flow solver applied to flow over hills, PhD thesis, https://backend.orbit.dtu.dk/ws/portalfiles/portal/12280331/Ris_R_827.pdf (last access: 1 November 2022) 1995. a
Sørensen, N.: HypGrid2D. A 2d mesh generator, Tech. rep., Risø National Laboratory, https://backend.orbit.dtu.dk/ws/portalfiles/portal/7750949/RIS_R_1035.pdf (last access: 1 November 2022), ISBN 8755023681, 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, https://doi.org/10.1002/we.64, 2002. a
Williamson, C. and Govardhan, R.: VortexInduced Vibrations, Annu. Rev. Fluid Mech., 36, 413–55, https://doi.org/10.1146/annurev.fluid.36.050802.122128, 2004. a
Zhang, M., Song, Y., Wang, J., Abdelkefi, A., Yu, H., and Yu, H.: Vortexinduced vibration of a circular cylinder with nonlinear stiffness : prediction using forced vibration data, Nonlinear Dynam., 108, 1867–1884, https://doi.org/10.1007/s11071022073327, 2022. a
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(7578 KB)  Fulltext XML