Forced-motion simulations of vortex-induced vibrations of wind turbine blades – a study of sensitivities

Vortex-induced vibrations on wind turbine blades are a complex phenomenon not predictable by standard engineering models. For this reason, higher-fidelity 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 vortex-induced vibration (VIV) phenomenon to a satisfying degree. The method studied is the so-called forced-motion (FM) approach, where the structural motion is imposed on the CFD blade surface through mode shape assumptions rather than fully coupled two-way 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 cross-sectional component of the flow is close to perpendicular to the chord of the blade. It is found that the low-inclination-angle and high-inclination-angle 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 Reynolds-averaged 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 scale-resolving 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 higher-fidelity hybrid RANS–large eddy simulation (LES) turbulence models, like improved delayed detached eddy simulation (IDDES), stress-blended eddy simulation (SBES) and k–ω scale-adaptive 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.

1 Introduction

Vortex-induced 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 lock-in 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, high-fidelity 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 so-called 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 }$.

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

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 forced-motion 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 forced-motion (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 well-established 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.

2 Methodology

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 forced-motion 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 Forced-motion method

In order to undertake high-fidelity modeling of VIVs without using a structural coupling framework, the forced-motion (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 low-speed 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⁻¹. 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.

Figure 2 First edgewise mode shape of the blade from HAWCStab2 along with the polynomial fit used in forced-motion simulations. The effect of the torsional component, which was less than 0.5^∘ at the tip for 1 m edgewise deflection, was found to be negligible in Grinderslev et al. (2022). The maximum edgewise deflection to the leading edge occurs together with the maximum flapwise tip deflection towards the pressure side.

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 finite-volume code based on structured grids that solves the incompressible Navier–Stokes equations with Reynolds-averaged 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. Second-order implicit backward iterative time stepping (dual time stepping) is used as the temporal discretization scheme.

In this study, simulations are based on unsteady Reynolds-averaged Navier–Stokes (URANS) k–ω shear stress transport (SST) (Menter, 1993) along with the higher-fidelity 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 fourth-order 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 in-house 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.

Figure 3 Surface grid and volume mesh hyperbolically grown from the surface for EllipSys3D. Only every second line is shown for clarity.

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 block-structured 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 non-orthogonal grids at the surface, the normal grid lines are rotated according to the present surface-normal 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 far-field 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 steady-state computation. The blending function for translatoric deformations is $B_{\mathrm{def}}=\mathrm{1}-\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}-\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 multi-physical fluid dynamics. Only the incompressible Navier–Stokes equations with URANS and hybrid URANS–LES formulation are used in this work. Ansys CFX is a finite-volume code based on structured, unstructured and hybrid grids. Ansys CFX uses a coupled solver (CFX-Solver, 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 scale-adaptive simulation (SAS) (Egorov et al., 2010) and stress-blended 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 “close-to” LES formalism is used.

For k–ω and SAS, the numerical schemes used are second-order in space and second-order backward Euler in time. For SBES, the spatial numerical scheme is changed to a bounded central-difference scheme (Leonard, 1991), switching between a second-order central-difference scheme and first-order 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 vertex-centered orthogonality (VCO, area-weighted 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.

Figure 4 (a) Near-blade volume grid in the CFX setup. The structured mesh is grown from the blade surface and switched to unstructured at ≈20 m from the surface. (b) Diffusion of the mesh displacement norm emanating from the blade surface at a given time step.

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 unstructured-grid 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 fourth-order accurate in EllipSys3D and second-order 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 higher-fidelity 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.

Figure 5 Discrepancies of the blade surface at the tip between the two setups: EllipSys3D (rounded tip) and CFX (flat tip).

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.

$$\begin{array}{}\text{(1)}& {\displaystyle }{P}_{A,\mathrm{SPAN}}\left(y\right)={\displaystyle \frac{\mathrm{1}}{Tn}}\underset{t=t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+Tn}{\int }}F\left(y\right)\dot{u}\left(y\right)\mathrm{d}t,\text{(2)}& {\displaystyle }{P}_{A,\mathrm{TOT}}=\int P_{A,\mathrm{SPAN}}\left(y\right)\mathrm{d}y,\end{array}$$

where t₀ 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 single-degree-of-freedom systems, the energy E_D dissipated by damping during one cycle of harmonic vibration of frequency ω is given by Paz (2012):

$$\begin{array}{}\text{(3)}& E_{\mathrm{D}}=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/\mathit{\omega }}{\int }}\left(c\dot{u}\right)\mathrm{d}u=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/\mathit{\omega }}{\int }}c{\dot{u}}^{\mathrm{2}}\mathrm{d}t.\end{array}$$

Given $\dot{u}\left(t\right)=A\mathit{\omega }\cos (\mathit{\omega }t-\mathit{\varphi })$, A being the amplitude, Eq. (3) becomes

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{E}_{\mathrm{D}}& =c{A}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/\mathit{\omega }}{\int }}{\cos }^{\mathrm{2}}(\mathit{\omega }t-\mathit{\varphi })\mathrm{d}t=\mathit{\pi }c\mathit{\omega }{A}^{\mathrm{2}}\\ & \approx \mathrm{2}\mathit{\pi }\mathit{\zeta }k{A}^{\mathrm{2}}\approx \mathrm{2}\mathit{\pi }\mathit{\zeta }{\mathit{\omega }}^{\mathrm{2}}m{A}^{\mathrm{2}},\end{array}\end{array}$$

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 }$.

$$\begin{array}{}\text{(5)}& P=c/\mathrm{2}{\mathit{\omega }}^{\mathrm{2}}{A}^{\mathrm{2}}\approx \mathit{\omega }\mathit{\zeta }k{A}^{\mathrm{2}}\approx {\mathit{\omega }}^{\mathrm{3}}\mathit{\zeta }m{A}^{\mathrm{2}}=F_{\mathrm{STRUC}}{A}^{\mathrm{2}},\end{array}$$

where F_STRUC is a structural damping dissipation factor as also described in Grinderslev et al. (2022).

An equivalent single-degree-of-freedom 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,

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathbf{M}}_{\mathrm{modal}}={\mathbf{V}}^T\mathbf{MV},\text{(7)}& {\displaystyle }& {\displaystyle }{\mathbf{C}}_{\mathrm{modal}}={\mathbf{V}}^T\mathbf{CV},\text{(8)}& {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{modal}}={\mathbf{V}}^T\mathbf{KV},\end{array}$$

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

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{eff}}={\mathbf{M}}_{\mathrm{modal}}(i,i),\text{(10)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{eff}}={\mathbf{C}}_{\mathrm{modal}}(i,i),\text{(11)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{eff}}={\mathbf{K}}_{\mathrm{modal}}(i,i).\end{array}$$

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⁻². This deviates from the value of 540 W m⁻² 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 Results

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 higher-fidelity 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, E-A, of ≈12 million cells to the finest case, E-E, of 2281 million cells. Obviously, a grid with more than 2×10⁹ cells is more an academic case than a practical case due to the corresponding immense computational cost for both running and post-processing. 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 lower-inclination cases with 30^∘ and 40^∘ inclination. For the higher-inclination 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 lower-inclination cases are more sensitive to the grid and turbulence model than higher-inclination cases, as the amount of chaotic natural shedding in low-inclination 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 coarse-grid setups A and B but negative power for the finer cases. In the specific case, this is conservative, as one would “over-design” 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 finer-grid 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 low-inclination 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.

Table 1 Grid refinement cases for the EllipSys3D setup.

Table 2 Grid refinement cases for the CFX setup. The total number of cells is given for the structured part of the mesh.

Figure 6 Total aerodynamic power per cycle, P_A,TOT, for varying grid refinements and flow inclination angles. $P=\mathrm{100}{}^{\circ }$ for all cases.

Figure 7 Total aerodynamic power per cycle, P_A,TOT, for varying grid refinements and turbulence models. Flow case P100I30.

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 higher-fidelity 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.

Figure 8 Accumulated power in time along with spanwise distribution of average aerodynamic power injection for various time steps in the CFD solvers for flow case P100I50. For EllipSys3D the baseline mesh and IDDES were used. For CFX the baseline mesh and SBES were used.

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 sub-iterations needed.

Figure 9 Spanwise distribution of aerodynamic power for flow case P100I50. (a) Average power per cycle P_A,SPAN; (b) P_A,SPAN accumulated over span.

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.

Figure 10 Vorticity fields resulting from the turbulence models used along with isosurfaces of q criterion = 5.0. Flow case P100I50.

Figure 11 Vorticity fields resulting from the turbulence models used along with isosurfaces of q criterion = 5.0. Flow case P100I30.

3.2 Turbulence model dependency

3.2.1 High inclination

The high-inclination 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 small-structure 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 low-inclination 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 small-scale vortices are created, without large-scale spatial or temporal correlation.

Figure 12 Zoomed-in view of the outer part of the blade for URANS (DTU setup), IDDES and SAS turbulence; q criterion = 20.0 and isosurfaces are colored by vorticity between −1 and 1 rad s⁻¹. Flow case P100I30.

These small-scale vortices result in negative accumulated power along the span. The spanwise power distribution is therefore also much less in agreement between the lower- and higher-fidelity turbulence models than in the high-inclination case. The lower-fidelity URANS turbulence models predict high spanwise correlation, resulting in a high power injection between 1560 and 2000 W. The higher-fidelity 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 higher-fidelity 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 low-inclination flow.

Figure 13 Spanwise distribution of aerodynamic power for flow case P100I30. (a) Average power per cycle P_A,SPAN; (b) P_A,SPAN accumulated over span.

Figure 14 Effective power, $P_{\mathrm{EFF}}=P_{A,\mathrm{TOT}}-P_{\mathrm{STRUC}}$, for various amplitudes using SAS, SBES, IDDES and URANS turbulence modeling. Flow cases P100I30 (a) and P100I50 (b).

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 high-inclination 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 low-inclination case, however, it is again clear that the URANS model results in far-from-realistic scenarios due to the artificial vortex coherency created. For this flow scenario, all higher-fidelity 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 high-inclination case. A scenario in which higher-fidelity turbulence will predict vibrations where URANS turbulence does not is plausible. However, this scenario, if present, is still to be revealed to the authors.

4 Conclusions

A comprehensive study has been conducted, investigating the impacts of various simulation choices for vortex-induced vibrations of wind turbine blades. Common for all studies was the forced-motion 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 in-house 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 low-inclination 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 higher-fidelity DES-like turbulence models becomes too big, leading to similar results to those for pure URANS simulations. For finer grids, the higher-fidelity models resolve the more chaotic smaller-scale vortices, which breaks the coherence and power injection. For higher-inclination 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 vortex-induced 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.

5 Future studies

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 high-fidelity 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 forced-motion 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 reduced-order 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.

Appendix A: Acronyms and abbreviations

• CFD – computational fluid dynamics

• VIV – vortex-induced vibration

• FM – forced motion

• FSI – fluid–structure interaction

• LES – large eddy simulation

• IDDES – improved delayed detached eddy simulation

• URANS – unsteady Reynolds-averaged Navier–Stokes

• SAS – scale-adaptive simulation

• SBES – stress-blended 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 – vertex-centered orthogonality

• SST – shear stress transport