the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Forcedmotion simulations of vortexinduced vibrations of wind turbine blades – a study of sensitivities
Felix HoutinMongrolle
Niels Nørmark Sørensen
Georg Raimund Pirrung
Pim Jacobs
Aqeel Ahmed
Bastien Duboc
Vortexinduced vibrations on wind turbine blades are a complex phenomenon not predictable by standard engineering models. For this reason, higherfidelity computational fluid dynamics (CFD) methods are needed. However, the term CFD covers a broad range of fidelities, and this study investigates which choices have to be made when wanting to capture the vortexinduced vibration (VIV) phenomenon to a satisfying degree. The method studied is the socalled forcedmotion (FM) approach, where the structural motion is imposed on the CFD blade surface through mode shape assumptions rather than fully coupled twoway fluid–structure interaction. In the study, two independent CFD solvers, EllipSys3D and Ansys CFX, are used and five different turbulence models of varying fidelities are tested. Varying flow scenarios are studied with low to high inclination angles, which determine the component of the flow in the spanwise direction. In all scenarios, the crosssectional component of the flow is close to perpendicular to the chord of the blade. It is found that the lowinclinationangle and highinclinationangle scenarios, despite having a difference equivalent to up to only a 30^{∘} azimuth, have quite different requirements of both grid resolution and turbulence models. For high inclination angles, where the flow has a large spanwise component from the tip towards the root, satisfying results are found from quite affordable grid sizes, and even with unsteady Reynoldsaveraged Navier–Stokes (URANS) k–ω turbulence, the result is quite consistent with models resolving more of the turbulent scales. For low inclination, which has a high degree of natural vortex shedding, the picture is the opposite. Here, even for scaleresolving turbulence models, a much finer grid resolution is needed. This allows us to capture the many incoherent vortices, which have a large impact on the coherent vortices, which in turn inject power into the blade or extract power.
It is found that a good consistency is seen using different variations of the higherfidelity hybrid RANS–large eddy simulation (LES) turbulence models, like improved delayed detached eddy simulation (IDDES), stressblended eddy simulation (SBES) and k–ω scaleadaptive simulation (SAS) models, which agree well for various flow conditions and imposed amplitudes.
This study shows that extensive care and consideration are needed when modeling 3D VIVs using CFD, as the flow phenomena, and thereby solver requirements, rapidly change for different scenarios.
 Article
(7550 KB)  Fulltext XML
 BibTeX
 EndNote
Vortexinduced vibrations (VIVs) on wind turbine blades are a phenomenon gaining relevance as wind turbines become larger and more flexible. When the turbine is not in operation, for instance due to grid loss, maintenance or storm conditions or during erection, the blades can experience wind from various directions, which can result in large angles of attack that are close to perpendicular to the chord. In this range of wind directions, deep stall with a high degree of vortex shedding can occur, meaning that a risk of lockin between structural modes and shedding frequencies increases.
As VIVs are directly depending on vortex shedding frequency and phase between the corresponding loads and motion velocity, engineering models struggle to compute the phenomenon. It becomes especially problematic as the blade shape, by twist and chord length, changes over the blade span, making a simple Strouhal relationship analysis difficult. For this reason, highfidelity methods such as computational fluid dynamics (CFD) are needed. Examples of this are the works of Horcas et al. (2022a, b, 2020), who studied wind turbine blade VIVs through fluid–structure interaction (FSI) simulations coupling CFD with a structural solver. It was shown that branches of VIVs can be found for various flow angles, defined by socalled pitch and inclination angles. This flow direction definition is also used in the present study and is depicted in Fig. 1. As shown, the pitch angle, P, is the angle between flow and the chord of the root airfoil section. The pitch angle is defined similarly to a standard geometrical angle of attack, reaching $P=\mathrm{90}{}^{\circ}$ when the wind vector strikes the pressure side perpendicularly to the chord at the root. Inclination, 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 when the spanwise component of the wind flow is from the tip to the root and zero when the wind strikes the blade perpendicularly to the span. It is important to notice that various combinations of wind turbine settings, i.e., blade pitch, yaw and azimuth, can result in identical inclination and pitch angles, meaning these parameters are more general than a single turbine setting. In the following, flow cases will be named based on their pitch and inclinations angles; i.e., P100I30 means a pitch angle of $P=\mathrm{100}{}^{\circ}$ and an inclination angle of $I=\mathrm{30}{}^{\circ}$.
The positions of the VIV branches depend on the blade shape, structural properties and flow velocity. As shown in Horcas et al. (2022a, 2020), changing the shape of the blade by tip and/or flap altercations moves the regions of VIVs; however, it does not seem plausible, within realistic alternations, to remove the VIV risk entirely, as the branches seem rather to shift towards other flow angles and/or velocities.
In a recent study, Grinderslev et al. (2022) used the same setup as Horcas et al. (2022a, 2020) and showed that it is feasible to omit the coupling with the structural solver and replace it with an analytical imposed motion of the structural mode, at least when considering single wind turbine blades which are rigidly clamped at the root. The approach is to simulate the forced motion for various defined amplitudes, revealing a picture of the power injection by the aerodynamic loads. By comparison to the dissipated power from structural damping, an equilibrium vibration state can then be estimated. The benefit of such an approach is that no coupling framework is needed, and for a specific simulation the approach is likewise faster, as no time for the buildup of vibration is needed. This latter benefit, however, disappears when multiple simulations with various amplitudes are needed if the vibration development is of interest.
The forcedmotion method is not a new concept and has been used extensively, especially on cylinder VIVs (Placzek et al., 2009; Viré et al., 2020) but also for airfoils in 2D (Skrzypiński et al., 2014; Hu et al., 2021).
In the present study, the approach of forcedmotion (FM) CFD analysis of VIVs is studied further in various aspects. The influence of modeling schemes, turbulence models, grids and more is studied using two wellestablished CFD solvers: one used and developed by the Technical University of Denmark (DTU), EllipSys3D (Sørensen, 1995; Michelsen, 1992), and one commercial code used by Siemens Gamesa Renewable Energy (SGRE), Ansys CFX (v. 2021, R1; CFX, 2021). The aim is to provide knowledge about good practices when simulating VIVs for wind turbine blades. The present work uses the IEA 10 MW blade (Bortolotti et al., 2019), also studied in the aforementioned FSI studies by Horcas et al. (2022b) and Grinderslev et al. (2022).
The study shows that the chosen modeling approach has large effects on the computed power input for cases with low to medium inclination angles, where uncorrelated natural shedding occurs. For cases with high inclination angles, the sensitivity is found to be much lower, as the defined blade motion controls the flow pattern more in this region.
As two different codes are used, various combinations of grid methods, convective schemes, turbulence models and much more can be studied. The setups of the two solvers are based on the experiences, common practices and computational resources of the users (DTU uses EllipSys3D; SGRE uses Ansys CFX), albeit with a desire of being able to capture the same physics. Common in all simulations is the use of forcedmotion CFD simulations as described below in Sect. 2.1. In the following subsections, the fluid solver codes will be described along with the chosen models. Finally, the analysis methods will be described.
2.1 Forcedmotion method
In order to undertake highfidelity modeling of VIVs without using a structural coupling framework, the forcedmotion (FM) method is used. Here, it is assumed that the structural response of the blade seen during VIVs can be simplified to being purely the structural modes. This assumption works well when the mode shapes and natural frequencies are not altered significantly by the surrounding flow, which is the case for the current study with lowspeed airflow. In these simulations, the first edgewise mode has been chosen, as this is the mode being provoked by the investigated flow scenarios when using fully coupled FSI simulations (Horcas et al., 2022b; Grinderslev et al., 2022). In these studies the assumption of having modes that are close to purely structural has been validated for a wind speed of 18 m s^{−1}. For high wind speeds, the assumption might not hold as the aeroelastic mode shape moves away from the structural one.
The aeroelastic model of the IEA 10 MW wind turbine is publicly available from Bortolotti et al. (2019), and in the present study the aeroelastic stability tool HAWCStab2 (Hansen, 2004) was used on said model to extract the modes. The first edgewise mode has a frequency of 0.67 Hz, and the edge and flapwise motion components of the mode shape are depicted in Fig. 2, along with the corresponding polynomials fitted for use with the CFD solvers. The phase between the flapwise and edgewise amplitude is such that the maximum tip deformation towards the pressure side of the blade occurs at the same time as the maximum deformation towards the leading edge. The amplitude used in the present study is 1 m in the edgewise direction, except for amplitude sweeps presented in Sect. 3.3.
For the specific study, some assumptions are made to enable the FM approach. Firstly, as mentioned, it is assumed that the first structural edgewise mode shape is the only motion present. This means that no contribution from static loads or buffeting loads is included in the motion. This assumption aligns well with what was found in previous studies using fully coupled FSI (Horcas et al., 2022b; Riva et al., 2022). Another assumption made is the disregard of the torsional part of the mode shape. This assumption is made for practical reasons in terms of imposing motion in the CFD solvers. The effect has again been tested with FSI simulations which include torsion, and it was found that the consequence of not considering torsion is marginal. In the present case, the torsional component is less than half a degree at the tip for a 1 m amplitude. This is, however, not a general conclusion and is something that should be assessed for the specific blade and flow scenario considered. One reason that torsion has little effect in the present study is the angles of attack (AoAs) studied, which are close to perpendicular to the chord. Here, the aerodynamics are less sensitive to small changes in the AoA than for instance at stall onset.
2.2 EllipSys3D setup
The EllipSys3D CFD solver (Sørensen, 1995; Michelsen, 1994, 1992) is a finitevolume code based on structured grids that solves the incompressible Navier–Stokes equations with Reynoldsaveraged Navier–Stokes (RANS), large eddy simulation (LES) or hybrid turbulence equations. The solution algorithm is based on the SIMPLE algorithm in combination with Rhie–Chow interpolation to avoid pressure decoupling. Secondorder implicit backward iterative time stepping (dual time stepping) is used as the temporal discretization scheme.
In this study, simulations are based on unsteady Reynoldsaveraged Navier–Stokes (URANS) k–ω shear stress transport (SST) (Menter, 1993) along with the higherfidelity k–ωbased improved delayed detached eddy simulation (IDDES) (Gritskevich et al., 2012) for better resolution of turbulent structures created in the near wake of the blade.
For the URANS simulations the quadratic upstream interpolation for convective kinematics (QUICK) convective scheme is used, while for IDDES a combination of QUICK (in the RANS region) and fourthorder central difference (in the LES region) is used as described by Strelets (2001).
2.2.1 EllipSys3D grids
Various grids have been tested in the present study. All surface grids are based on the DTU inhouse Parametric Geometry Library (PGL) tool (Zahle, 2022), and volume grids are created through hyperbolic extrusion from the surface grid to a spherical domain with a radius of ≈700 m (≈ seven blade lengths), using the mesh tool HypGrid3D (Sørensen, 1998). Multiple grid refinements have been tested to study the sensitivity of the VIVs to the resolved vortices. The baseline mesh used, if not otherwise stated, has 512 cells spanwise along the blade, 256 cells chordwise and 256 cells normal to the surface. The total number of cells for the baseline mesh is 35.6 million cells. This mesh is finer than that used in previous publications (Horcas et al., 2020, 2022b; Grinderslev et al., 2022) and in this study was found to be necessary for certain flow scenarios; see Sect. 3.1. The first cell size normal to the surface is set to $\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ m, which ensures a y+ value of much less than 1. This cell height is a common choice in EllipSys3D for operational cases and is kept here, despite being unnecessarily low for the wind speed used, as the additional computational cost is low.
The grid deformation procedure in the EllipSys3D simulations is based on an explicit algebraic algorithm, transferring the deformation of the surface grid into the volume grid by a blending approach that exploits the blockstructured nature of the computational grid.
The deformed grid is computed by enforcing the Cartesian translation and deformation of the surface grid points along the grid lines normal to the surface. To avoid generating highly nonorthogonal grids at the surface, the normal grid lines are rotated according to the present surfacenormal direction. Using blending functions in the direction normal to the surface, it is assured that the grid translation and rotation are only enforced in the proximity of the surface of the geometry. This ensures that the original grid quality is conserved at the surface while preserving the original grid far away from the surface. In between the surface proximity region and the farfield region, a blending region is present where the grid quality risks deterioration in the case of large deformations if the blending is not adequately tuned. Typically, the surface deflections are enforced far away from the surface, while the rotations are limited to a region close to the surface. The blending is based on hyperbolic tangent functions, using the normalized curve length along the grid lines normal to the surface. The procedure can easily be tuned for specific cases by calibrating the blending function constants for a severe static deformation using a steadystate computation. The blending function for translatoric deformations is ${B}_{\mathrm{def}}=\mathrm{1}\mathrm{tanh}\left(a\right(\frac{s}{S}{)}^{b})$, with s being the curve length from the surfaces of the specific grid line and S being its total length. Factors a and b are the tuning parameters, which in this specific study were set as 1500 and 3.0 respectively. The orthogonality blending is done similarly: ${B}_{\mathrm{orth}}=\mathrm{1}\mathrm{tanh}\left(c\right(\frac{s}{S}{)}^{d})$, with tuning factors c and d here set to 15 000 and 1.05 respectively. To assure that the grid is not degenerating, a simple check for negative cell volumes is performed after each grid deformation. No checks for orthogonality are performed.
2.3 Ansys CFX setup
The Ansys CFX library (CFX, 2021) gathers a set of solvers to resolve different multiphysical fluid dynamics. Only the incompressible Navier–Stokes equations with URANS and hybrid URANS–LES formulation are used in this work. Ansys CFX is a finitevolume code based on structured, unstructured and hybrid grids. Ansys CFX uses a coupled solver (CFXSolver, 2021) combined with Rhie–Chow interpolation, which solves the hydrodynamic equations as a single system differing from the SIMPLE algorithm. The solver uses a multigrid accelerated incomplete lower–upper (Raw, 1996) factorization technique for solving the discrete system of linearized equations. It relies on an iterative process to approach the exact solution.
In this study, the simulations are based on URANS with a k–ω SST (Menter, 1993), SST scaleadaptive simulation (SAS) (Egorov et al., 2010) and stressblended eddy simulation (SBES) (Menter, 2018). SAS and SBES aim for a better resolution of the turbulent structures by either decreasing the added turbulence modeling or relying on LES turbulence models. They both rely on shielding functions to delimit volumes where a “closeto” LES formalism is used.
For k–ω and SAS, the numerical schemes used are secondorder in space and secondorder backward Euler in time. For SBES, the spatial numerical scheme is changed to a bounded centraldifference scheme (Leonard, 1991), switching between a secondorder centraldifference scheme and firstorder upwind scheme, based on the local convection boundedness criterion (Jasak et al., 1999).
2.3.1 Ansys CFX grids
For the CFX grid, a combination of structured and unstructured grids is used to keep more control near the blade while exploiting the unstructured expansion of the grid further away. The different meshes are generated using Pointwise v22.2, allowing the control of the structured mesh. The first cell size normal to the surface is set to $\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ m with a growth rate of 1.07, which ensures a y+ of 1 or less. Several grid resolutions are investigated in this work where only the structured part is refined; i.e., the background unstructured cell size is kept constant. The baseline structured mesh used, if not otherwise stated, has 500 cells spanwise, 544 cells chordwise and 152 cells normal to the surface. This adds up to ≈50 million cells total: 48 million in the structured region and 2 million in the unstructured region. The mesh quality is evaluated based on the cell length ratios in the chordwise and spanwise direction, keeping them below 1.5. In the structured region, near the blade, the vertexcentered orthogonality (VCO, areaweighted average of the orthogonality angles associated with each bounding face of the dual mesh control volume around the vertex; a 90^{∘} VCO represents perfect orthogonality) is kept higher than 20^{∘}. For the unstructured region, a Delaunay triangulation algorithm is used and a smooth transition from structured to unstructured is achieved by using a growth rate of 1.07 for the tetrahedron. An overview of the mesh is given in Fig. 4a. The domain is square with 500 m side lengths (≈ five blade lengths), and the blade is placed in the center. Side boundaries are velocity inlets and pressure outlets.
To take into account the motion of the selected mode shape, the mesh is deformed periodically at the mode frequency. The mesh deformation is computed only during the initialization step as the displacement is imposed. This deformation is computed by diffusing the displacement registered on the blade boundary to the neighboring mesh cells. To prevent any cell from folding over, a mesh stiffness is defined. This stiffness is set to increase near the blade boundaries at a cubic rate and after a distance to the blade boundary of 1 m. The obtained mesh displacement at a given time step is depicted in Fig. 4b. The VCO and negative cell volumes are monitored to ensure that the grid remains suited for resolving the flow.
2.4 Setup differences
The main difference between the two used CFD codes is their discretization methods, with EllipSys3D being a structured solver, while CFX uses unstructured grids. Both of these have pros and cons; the unstructuredgrid approach is more flexible in terms of grid manufacturing but often results in a slower performance. In this work, the grid close to the surface was chosen as structured for the CFX solver as well to avoid too rapid a dissipation of the shed vortices, which is found when using an unstructured approach. Further from the surface, the unstructured grid rapidly expands, limiting the cell count, i.e., ensuring faster computations. For the structured grid in EllipSys3D, an expansion of cells also happens when moving far from the surface. The number of cells used also varies between the two setups based on grid sensitivity studies; see more in Sect. 3.1. This is partly a consequence of the convective schemes having different orders of accuracy, being fourthorder accurate in EllipSys3D and secondorder accurate in CFX.
Domain shapes differ between the two, being spherical for EllipSys3D and square for the CFX setup. This should have no impact, as boundary conditions are far from the considered blade.
The turbulence models implemented in the two solvers differ but should have similar capabilities of capturing the vortex shedding with varying degrees of accuracy from URANS k–ω SST (Menter, 1993) in both solvers to the higherfidelity hybrid RANS–LES models like IDDES (Gritskevich et al., 2012; Menter et al., 2003; Shur et al., 2008) in EllipSys3D and SAS (Egorov et al., 2010) and SBES (Menter, 2018) in CFX.
Finally, the blade surface shows discrepancies at the tip. As the meshing methodology differs between the two setups, the tip cap surface used in CFX is flat, while it is rounded in EllipSys3D. This introduces an 8 cm difference corresponding to less than 0.1 % of the total blade length, which showed a low impact on the aerodynamic spanwise power distribution that is introduced in the following section.
2.5 Analysis
2.5.1 Aerodynamic power
In this paper, aerodynamic power will be defined as positive when injecting power into the structural system and negative when damping the structural response. When calculating aerodynamic power, the mean power over n full motion cycles is considered; see Eqs. (1) and (2). Through one motion period, there might be both a positive and a negative contribution of power, and so the total power over the full cycle needs to be considered to assess whether the structural response is excited or damped. The power, P_{A,TOT}, is total power for the full blade, meaning that power is found spanwise, P_{A,SPAN}, and integrated over the length of the blade.
where t_{0} is the start time of the integration, T is the time for one motion period and n is the number of full cycles considered. In this work, a minimum of five full cycles, after convergence is obtained, is used to analyze the aerodynamic power. F(y) and $\dot{u}\left(y\right)$ are the spanwise force distribution and structural velocity along the blade span y.
When considering the risk of VIVs, it is important to realize that the total power P_{A,TOT} is the driving factor. If this is negative, the vibration is damped, despite there being positive values of P_{A,SPAN} at some spanwise positions.
2.5.2 Assessment of VIV amplitudes
As with the aerodynamic power injected into the structure for a given amplitude, the corresponding power dissipated by structural damping can be found and compared to the aerodynamic power to assess whether the operating point is stable or not.
Structural damping is estimated using modal analysis. For singledegreeoffreedom systems, the energy E_{D} dissipated by damping during one cycle of harmonic vibration of frequency ω is given by Paz (2012):
Given $\dot{u}\left(t\right)=A\mathit{\omega}\mathrm{cos}(\mathit{\omega}t\mathit{\varphi})$, A being the amplitude, Eq. (3) becomes
where $\mathit{\zeta}=c/{c}_{\mathrm{cr}}\ll \mathrm{1}$ is the damping ratio, $\mathit{\omega}=\sqrt{k/m}\sqrt{\mathrm{1}{\mathit{\zeta}}^{\mathrm{2}}}\approx \sqrt{k/m}$ is the damped natural frequency and ${c}_{\mathrm{cr}}=\mathrm{2}\sqrt{km}$ is the critical damping. Parameters m, k and c are mass, stiffness and structural damping respectively.
To obtain the power dissipation, the energy needs to be divided by the period $T=\mathrm{1}/f=\mathrm{2}\mathit{\pi}/\mathit{\omega}$.
where F_{STRUC} is a structural damping dissipation factor as also described in Grinderslev et al. (2022).
An equivalent singledegreeoffreedom system for a given mode can be constructed using modal analysis; thus, the stiffness and damping are replaced by modal stiffness and modal damping for a given mode to compute the energy dissipated when the structure moves by a unit amplitude with a certain mode shape.
With eigenmatrix V containing in each column the eigenvector V_{i} for a given mode i,
where the first edgewise eigenvector has been scaled to correspond to 1 m blade tip deflection in the edgewise direction. Thus the effective mass, damping and stiffness values to be used in Eq. (5) are
It is also evident from Eq. (5) that structural damping is a quadratic function of the amplitude of the displacement.
Using the method above and HAWCStab2 (Hansen, 2004), the structural damping dissipation factor F_{STRUC} of the blade undergoing the first edgewise mode was found to be 495.8 W m^{−2}. This deviates from the value of 540 W m^{−2} that was found in previous FSI studies (Grinderslev et al., 2022). The reason for this deviation is likely a combination of the calculation method and deviations of the aeroelastic deflection shapes observed in the FSI simulations from the purely structural mode shapes investigated here. Finally, numerical damping in the loosely coupled FSI framework could play a role; however, this effect was found to be low by Heinz et al. (2016) when developing the framework. Considering the uncertainties in the structural damping of wind turbine blades, the agreement within 10 % is found to be acceptable.
3.1 Grid and time dependency
3.1.1 Grid sensitivity
Various grid configurations have been tested in the present study using varying turbulence models. In the higherfidelity turbulence models (SAS, SBES and IDDES), the resolved length scale in the LES region depends on the grid cell size itself, meaning that large changes in the grid can lead to large changes in the resulting flow. Two separate grid studies were conducted: first, the EllipSys3D solver using the IDDES turbulence models with different inclination angles and largely varying grid resolutions. The grid setups tested are defined in Table 1 and go from the coarsest case, EA, of ≈12 million cells to the finest case, EE, of 2281 million cells. Obviously, a grid with more than 2×10^{9} cells is more an academic case than a practical case due to the corresponding immense computational cost for both running and postprocessing. Luckily, it is also found to be unnecessarily fine in the current study, as depicted in Fig. 6, which shows the total power for various flow inclinations as result of the grid. For all the cases, a 1 m amplitude has been imposed. As seen, the resulting power becomes close to stable from grid case C and finer, and the sensitivity seems to be highest for the lowerinclination cases with 30^{∘} and 40^{∘} inclination. For the higherinclination cases $I=\mathrm{50}{}^{\circ}$ and $I=\mathrm{60}{}^{\circ}$, the sensitivity is in general low, which was also found in previous work by Horcas et al. (2022b).
Note that for the sake of visibility, the figure has the number of cells spanwise on the x axis, despite cases D and E having varying numbers of cells in the normal and chordwise directions as well, for the sake of grid quality.
Physically, it makes sense that lowerinclination cases are more sensitive to the grid and turbulence model than higherinclination cases, as the amount of chaotic natural shedding in lowinclination cases is quite high. For higher inclinations, there seems to be much more shedding that is correlated with the motion of the blade, meaning that larger, more ordered vortices are resolved. For the flow case of $I=\mathrm{30}{}^{\circ}$, the importance of considering the grid is evident, as the CFD predicts positive power injection for the coarsegrid setups A and B but negative power for the finer cases. In the specific case, this is conservative, as one would “overdesign” the turbine if considering it in the design. However, it is not a given that the opposite cases could not exist, where the positive injection of power would only be captured for finergrid setups as is almost the case for $I=\mathrm{40}{}^{\circ}$.
This finding indicates that the previously found VIV risk mapping from Horcas et al. (2022b) shows false positives in the lowinclination region around 30^{∘}, as the grid in that study was coarser than what is here found to be necessary. However, the main risk region found in the mapping at higher angles of inclination is valid with these findings.
Secondly, a grid study using various turbulence models was conducted using the CFX setup and the turbulence models URANS, SAS and SBES for the inclination angle of $I=\mathrm{30}{}^{\circ}$, which as mentioned was found to be quite grid sensitive for the EllipSys3D setup. The grid setups for the CFX cases are described in Table 2. The results are shown in Fig. 7 along with the corresponding IDDES cases from EllipSys3D.
As seen, the behavior of the higherfidelity models, SAS, SBES and IDDES, is very similar, and a large dependency on the grid setups is found. The URANS case, however, does not see this dependency but appears to overshoot the aerodynamic power injection for all grids considered.
3.1.2 Time step sensitivity
The sensitivity to time step size for the EllipSys3D simulations was studied for the flow case P100I50 with a 1 m amplitude. The baseline time step was set to $\mathrm{6}\times {\mathrm{10}}^{\mathrm{3}}$ s, which was found to be suitable in earlier work by Horcas et al. (2022b) and results in 250 time steps per motion cycle for the first edgewise mode. The time step was varied to half and double the baseline, and the resulting power injection distribution along the span is given in Fig. 8. It has in general been found that a deep convergence is needed to capture the power injection well, as it is directly dependent on both force amplitude and the phase between motion and force. In the EllipSys3D setup, each time step has between 5 and 20 inner iterations, which dynamically change based on the convergence of the flow residuals. This means that when the time step is low enough to capture the physics, even lower time steps would typically reduce the number of subiterations needed.
For the CFX setup, a similar sensitivity study is performed. The baseline time step is less strict, reaching $\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{2}}$ s, resulting in 100 time steps per motion cycle for the same first edgewise mode. It is then reduced by 2 and 4, where the resulting power injection along the span is shown in Fig. 8. In a similar way to EllipSys3D, the time step inner iterations fluctuate according to the flow residuals' convergence. A limit of five inner iterations was used for CFX. As seen in Fig. 8, the two solvers obtain very similar results when using time steps half the size of the baseline choices. When considering power as the objective, the total convergence of results is hard to obtain, as it is extremely sensitive to the phase between the found forces and the motion velocity as well as the force amplitude. In this case the difference between the accumulated power over a full cycle is less than 7 % between the two smallest time steps for the EllipSys3D setup and even less for the CFX setup. This is deemed acceptable for the purpose of this investigation, as the uncertainty in the corresponding dissipated power from structural damping is likely much higher. As seen in the middle panels of Fig. 8, the main difference in accumulated power stems from the range between 40 and 70 m along the span, where the majority of aerodynamic power is inputted. Closer to the tip, the power aligns for all cases.
3.2 Turbulence model dependency
3.2.1 High inclination
The highinclination case has flow coming with an inclination angle I of 50^{∘} and pitch angle P of 100^{∘}. In this scenario, the shedding is quite correlated for all turbulence models, as depicted in Fig. 10, but with the highest degree of correlation for URANS turbulence modeling, as this cannot resolve the smallstructure vortices. The spanwise power distribution is quite similar, no matter the turbulence model, and the accumulated power injected into the blade is also close between all methods, yielding a minimum of 1280 W from EllipSys3D URANS and a maximum of 1450 W from CFX SAS; see Fig. 9.
3.2.2 Low inclination
The lowinclination case has flow coming with an inclination angle I of 30^{∘} and pitch angle P of 100^{∘}. Here, a large degree of natural, more chaotic, shedding occurs, which is seen to be quite different between the various turbulence models; see Figs. 11 and 12. For IDDES, SAS and SBES, a high number of smallscale vortices are created, without largescale spatial or temporal correlation.
These smallscale vortices result in negative accumulated power along the span. The spanwise power distribution is therefore also much less in agreement between the lower and higherfidelity turbulence models than in the highinclination case. The lowerfidelity URANS turbulence models predict high spanwise correlation, resulting in a high power injection between 1560 and 2000 W. The higherfidelity models, however, predict the situation to be positively damped with an accumulated power of ≈ −300 W; see Fig. 13.
It is important to note that for IDDES, SBES and SAS, this flow case resulted in a positive power injection for lower grid resolution – presented in Figs. 6 and 7, which depict the evolution of the total aerodynamic power injection according to the grid resolution and turbulence models for this low inclination. In Fig. 7, URANS results seem to have converged numerically yet present results in direct opposition to those of SBES, SAS and IDDES once grid convergence is achieved. This is likely due to the fact that these higherfidelity models are blending URANS and LES models and, with lower grid resolution, the URANS region is larger, meaning the spanwise vortex shedding becomes more correlated for coarser grid resolutions. This shows that URANS, irrespective of the grid resolution, is not suited for simulating the VIV of blades with lowinclination flow.
Figures 11 and 12 show the vortices developed in the near wake captured by the various turbulence models. In Fig. 12 the outer approximately 40 % of the blade is shown and it is clear how the URANS turbulence model creates much more coherency between the vortices than what is seen for IDDES and SAS, which both create more incoherent natural shedding.
3.3 Resulting vibrations
Using both CFD setups, sweeps of amplitudes of up to 2 m were conducted using the various turbulence models. URANS and IDDES simulations were conducted using the EllipSys3D setup, while SBES and SAS simulations were performed on the CFX setup. By these sweeps an approximation of the vibration level can be given, using the structural damping of the considered blade. As presented in Sect. 2.5.2, the power dissipated by structural damping is proportional to the square of the vibration amplitude. For the considered blade a relation between dissipated power and amplitude was found: ${P}_{\mathrm{STRUC}}=\mathrm{495.8}\cdot {{A}_{\mathrm{tip}}}^{\mathrm{2}}$. By this relation, an effective power ${P}_{\mathrm{EFF}}={P}_{A,\mathrm{TOT}}{P}_{\mathrm{STRUC}}$ can determine whether the blade is in a stable or unstable situation, as depicted in green and red regions respectively of Fig. 14.
As seen for the highinclination case with $I=\mathrm{50}{}^{\circ}$, a similar trend of effective power is found between the various models, varying mostly at higher amplitudes. The equilibrium points of said simulations lie between ≈1.25 m for SAS and ≈1.6 m for IDDES. For the lowinclination case, however, it is again clear that the URANS model results in farfromrealistic scenarios due to the artificial vortex coherency created. For this flow scenario, all higherfidelity simulations (IDDES, SBES and SAS) lie well within the stable region, meaning that no vibrations should occur. URANS simulations, however, show a high injection of power all the way to 1.25 m amplitude. One could state that this is acceptable as the result is conservative, but note that this was not the case for the highinclination case. A scenario in which higherfidelity turbulence will predict vibrations where URANS turbulence does not is plausible. However, this scenario, if present, is still to be revealed to the authors.
A comprehensive study has been conducted, investigating the impacts of various simulation choices for vortexinduced vibrations of wind turbine blades. Common for all studies was the forcedmotion CFD approach, where the structural first edgewise mode was imposed as a motion in the CFD simulations, which in earlier work has been shown to be feasible. Two independent CFD methodologies were used: the DTU inhouse EllipSys3D solver and the commercial Ansys CFX solver used at Siemens Gamesa. Various grid strategies and turbulence models were tested and compared, showing a high degree of sensitivity for lowinclination flow, meaning spanwise flow closer to perpendicular to the span rather than along the span. It is found that for these inclinations, care is needed regarding the selection of the turbulence model and grid. The observed differences are due to the artificial coherency in the vortex structures created by unsteady RANS models, leading to a high input of aerodynamic power. For coarse grids, the URANS region of the higherfidelity DESlike turbulence models becomes too big, leading to similar results to those for pure URANS simulations. For finer grids, the higherfidelity models resolve the more chaotic smallerscale vortices, which breaks the coherence and power injection. For higherinclination cases, the sensitivities to the grid and turbulence models are much lower, as the degree of chaotic natural shedding is low compared to the coherent structures, which can be resolved fairly well even using URANS turbulence.
This leads to the main conclusion of the present study: a lot of care needs to be taken when simulating vortexinduced vibrations of wind turbine blades. Various conditions will need separate sensitivity investigations in order to ensure the accuracy of the results. This is important since it is otherwise a risk that computations that are much too heavy are conducted for cases that do not need it. Even worse, the computations that were found to be well resolved for one case may fail to predict the VIVs in other cases.
The topic of VIVs is becoming increasingly relevant with the increasing sizes of wind turbines, and much more research is needed. As a continuation of the current study, an expansion on the parametric space is needed to make final conclusions on turbulence models and grid requirements. This study shows that the necessity of highfidelity turbulence and fine grids is highly dependent on flow scenario. As of now, no general rule of thumb about how and when to use various models can be justified. This would need a larger mapping of flow cases and rotor designs. However, the tendency seems to be that the need for higher accuracy increases with the degree of natural shedding along the blade, meaning that low inclinations are more difficult to compute correctly.
Even though forcedmotion simulations have the possibility of being more efficient than fully coupled FSI simulations, the simulation time needed for broad mappings is still high, especially if various amplitudes or flow velocities are needed. By use of reducedorder modeling, the number of simulations needed could possibly be reduced significantly.
The current study only considers clamped single blades undergoing the first edgewise blade mode vibration. As wind turbines are coupled systems, the coupled rotor modes should likewise be studied. This highly increases the complexity as the modes of a wind turbine are many, and the motion of the individual blades will then depend on the azimuth position.

CFD – computational fluid dynamics

VIV – vortexinduced vibration

FM – forced motion

FSI – fluid–structure interaction

LES – large eddy simulation

IDDES – improved delayed detached eddy simulation

URANS – unsteady Reynoldsaveraged Navier–Stokes

SAS – scaleadaptive simulation

SBES – stressblended eddy simulation

QUICK – quadratic upstream interpolation for convective kinematics

AoA – angle of attack

DTU – Technical University of Denmark

SGRE – Siemens Gamesa Renewable Energy

IEA – International Energy Agency

PGL – Parametric Geometry Library

VCO – vertexcentered orthogonality

SST – shear stress transport
The geometry and the structural model of the considered wind turbine are publicly available at https://github.com/IEAWindTask37/IEA10.0198RWT (IEA Wind Task 37, 2023). However, fluid solvers are licensed. The data that support the findings of this study are available upon reasonable request to the corresponding author.
CG conducted EllipSys3D simulations, did analysis of the results and did the main writing of the paper. FHM did the majority of the CFX calculations, analysis and methodology implementation in CFX, as well as participating in the redaction and reviewing of the paper. NNS conducted EllipSys3D simulations and took part in the analysis and the reviewing of the paper. GRP conducted HAWCStab2 analysis and damping estimations, supported on analysis of EllipSys3D simulations, and reviewed the paper. Valuable support was provided by PJ, AA and BD regarding the analysis of CFX results, methodology implementation in CFX, and reviewing the paper. All authors contributed to 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 made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Computational resources were provided by the DTU Risø cluster Sophia (DTU Computing Center, 2021).
This research is part of the PRESTIGE project supported by Innovation Fund Denmark (grant no. 909000025B).
This paper was edited by Johan Meyers and reviewed by two anonymous referees.
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: 24 October 2023), 2019. a, b
CFX, A.: Modeling guide, Release 21R2, https://dl.cfdexperts.net/cfd_resources/Ansys_Documentation/CFX/Ansys_CFXSolver_Modeling_Guide.pdf (last access: 24 October 2023), 2021. a, b
CFXSolver, A.: Theory guide, Release 21R2, https://dl.cfdexperts.net/cfd_resources/Ansys_Documentation/CFX/Ansys_CFXSolver_Theory_Guide.pdf (last access: 24 October 2023), 2021. a
DTU Computing Center: DTU Computing Center resources, https://doi.org/10.48714/DTU.HPC.0001, 2021. a
Egorov, Y., Menter, F., Lechner, R., and Cokljat, D.: The scaleadaptive simulation method for unsteady turbulent flow predictions. part 2: Application to complex flows, Flow, Turbulence and Combustion, 85, 139–165, https://doi.org/10.1007/s1049401092654, 2010. a, b
Grinderslev, C., Nørmark Sørensen, N., Raimund Pirrung, G., and González Horcas, S.: Multiple limit cycle amplitudes in highfidelity predictions of standstill wind turbine blade vibrations, Wind Energ. Sci., 7, 2201–2213, https://doi.org/10.5194/wes722012022, 2022. a, b, c, d, e, f, g
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, Turbulence and Combustion, 88, 431–449, https://doi.org/10.1007/s1049401193784, 2012. a, b
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, b
Heinz, J., Sørensen, N., and Zahle, F.: Fluidstructure interaction computations for geometrically resolved rotor simulations using CFD, Wind Energy, 19, 2205–2221, 2016. a
Horcas, S. G., 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., Madsen, M., Sørensen, N.N. Zahle, F., and Barlas, T.: Influence of the installation of a trailing edge flap on the vortex induced vibrations of a wind turbine blade, J. Wind. Eng. Ind. Aerod., 229, 105118, https://doi.org/10.1016/j.jweia.2022.105118, 2022a. a, b, c
Horcas, S. G., Sørensen, 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, 054105, https://doi.org/10.1063/5.0088036, 2022b. a, b, c, d, e, f, g, h, i
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, 033306, https://doi.org/10.1063/5.0040509, 2021. a
IEA Wind Task 37: IEA10.0198RWT, GitHub [code], https://github.com/IEAWindTask37/IEA10.0198RWT (last access: 26 October 2023), 2023. a
Jasak, H., Weller, H., and Gosman, A.: High resolution NVD differencing scheme for arbitrarily unstructured meshes, Int. J. Numer. Meth. Fl., 31, 431–449, 1999. a
Leonard, B.: The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection, Comput. Method. Appl. M., 88, 17–74, 1991. a
Menter, F.: Zonal Two Equation Kappa–Omega Turbulence Models for Aerodynamic Flows, in: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, 6–9 July 1993, Orlando, FL, USA, https://doi.org/10.2514/6.19932906, 1993. a, b, c
Menter, F.: StressBlended Eddy Simulation (SBES) – A New Paradigm in Hybrid RANSLES Modeling, in: Progress in Hybrid RANSLES Modelling, edited by: Hoarau, Y., Peng, S.H., Schwamborn, D., and Revell, A., 27–37, Springer International Publishing, Cham, https://doi.org/10.1007/9783319700311_3, 2018. a, b
Menter, F., Kuntz, M., and Langtry, R.: Ten Years of Industrial Experience with the SST Turbulence Model, in: Proceedings of the 4th International Symposium on Turbulence, Heat and Mass Transfer, Begell House Inc., West Redding, 625–632, 2003. 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: 24 October 2023), 1992. a, b
Michelsen, J.: Block Structured Multigrid Solution of 2D and 3D Elliptic PDE's., Tech. rep., Technical University of Denmark, 1994. a
Paz, M.: Structural dynamics: theory and computation, Springer Science & Business Media, https://doi.org/10.1007/9783319947433, 2012. 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, Computers and Fluids, 38, 80–100, https://doi.org/10.1016/j.compfluid.2008.01.007, 2009. a
Raw, M.: Robustness of coupled algebraic multigrid for the NavierStokes equations, in: 34th Aerospace sciences meeting and exhibit, American Institute of Aeronautics and Astronautics, p. 297, https://doi.org/10.2514/6.1996297, 1996. 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
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
Skrzypiński, W., Gaunaa, M., Sørensen, N., Zahle, F., and Heinz, J.: Vortexinduced vibrations of a DU96W180 airfoil at 90^{∘} angle of attack, Wind Energy, 17, 1495–1514, https://doi.org/10.1002/we.1647, 2014. a
Strelets, M.: Detached eddy simulation of massively separated flows, 39th Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, https://doi.org/10.2514/6.2001879, 2001. a
Sørensen, N.: General purpose flow solver applied to flow over hills, PhD thesis, Risø National Laboratory, https://backend.orbit.dtu.dk/ws/portalfiles/portal/12280331/Ris_R_827.pdf (last access: 24 October 2023), 1995. a, b
Sørensen, N.: HypGrid2D. A 2d mesh generator, Tech. rep., Risø National Laboratory, ISBN 8755023681, https://backend.orbit.dtu.dk/ws/portalfiles/portal/7750949/RIS_R_1035.pdf (last access: 24 October 2023), 1998. a
Viré, A., Derksen, A., Folkersma, M., and Sarwar, K.: Twodimensional numerical simulations of vortexinduced vibrations for a cylinder in conditions representative of wind turbine towers, Wind Energ. Sci., 5, 793–806, https://doi.org/10.5194/wes57932020, 2020. a
Zahle, F.: PGL – Parametric Geometry Library, https://gitlab.windenergy.dtu.dk/frza/PGL (last access: 24 October 2023), 2022. a