Two-dimensional numerical simulations of vortex-induced vibrations for a cylinder in conditions representative of wind turbine towers

Vortex-induced vibrations (VIVs) of wind turbine towers can be critical during the installation phase, when the rotor–nacelle assembly is not yet mounted on the tower. The present work uses numerical simulations to study VIVs of a two-dimensional cylinder in the transverse direction under flow conditions that are representative of wind turbine towers both from a fluid dynamics and structural dynamics perspective. First, the numerical tools and fluid–structure interaction algorithm are validated by considering a cylinder vibrating freely in a laminar flow. In that case, both the motion amplitude and frequency are shown to agree well with previous results from the literature. Second, VIVs are modelled in the turbulent supercritical regime using unsteady Reynolds-averaged Navier–Stokes equations. In this context, the turbulence model is first validated against flow past a stationary cylinder with a high Reynolds number. Then, the results from forced vibrations are validated against experimental results for a range of reduced frequencies and velocities. It is shown that the behaviour of the aerodynamic damping changes with the frequency ratio and can therefore lead to either self-limiting or self-exciting VIVs when the cylinder is left to freely vibrate. Finally, results are shown for a freely vibrating cylinder under realistic flow and structural conditions. While a clear lock-in map is identified and shows good agreement with published numerical and experimental data, the work also highlights the unsteady nature of the aerodynamic forces and motion under certain operating conditions.

1 Introduction

Vortex-induced vibrations (VIVs) are structural vibrations that can occur due to the shedding of flow vortices when a fluid flow passes around a structure. Due to this fluid–structure interaction (FSI) phenomenon, a synchronization (also called lock-in) of the vortex shedding and the structural motion can occur for certain flow conditions and/or structural properties, leading to premature fatigue failures. Vortex-induced vibrations occur in many engineering applications, such as suspension bridges, marine risers, and industrial chimneys. In the context of wind energy, such vibrations have been observed numerically both in the tower (Livanos, 2018; Derksen, 2019) and the wind turbine rotor (Heinz et al., 2016; Horcas et al., 2019). The increase in rotor height and diameter makes VIVs an increasingly problematic issue when designing the new generation of wind turbines. This work is motivated by the specific problem of VIVs of wind turbine towers during their installation. In that phase, when the tower does not yet support the rotor–nacelle assembly, it acts as a beam clamped at one of its ends and is subjected to wind flow. Because of the vortex shedding developing in the tower wake, the tower may start to oscillate. This is particularly critical at three different stages of the installation process, namely when the tower (i) stands on the quayside, (ii) is transported on a vessel offshore, and (iii) is installed on the offshore foundation in the absence of the rotor–nacelle assembly.

The design of wind turbine towers depends on a number of parameters (Dykes et al., 2018). The tower's natural frequency, which is important for the dynamic response of the structure to vortex-induced excitations, depends on the tower diameter, thickness, material, and geometry (tapered). Various combinations of these parameters are therefore critical for the occurrence of VIVs. For the current towers used by Siemens Gamesa Renewable Energy, it is observed that tower diameters in the range of $\mathrm{4.5}\mathrm{m}\le D\le \mathrm{6.5}\mathrm{m}$ and with a first bending frequency of $\mathrm{0.5}\mathrm{Hz}\le {\mathit{\omega }}_{bm\mathrm{1}}\le \mathrm{1.1}\mathrm{Hz}$ are most critical for the occurrence of VIVs. However, as mentioned above, this is thickness-, material-, and shape-dependent. Furthermore, the trend over the years has been to reduce the tower-clamped first-mode frequency. This brings challenges as the critical resonance velocities are at a lower wind speed, which results in higher VIV incidences.

In practice, wind turbine towers are usually tapered cylinders with sections of discrete diameters. This shape, together with wind shear, might influence the VIV behaviour of the structure compared to that of circular cylinders. This has already been investigated to some extent in the literature (Balasubramanian et al., 2001; Hover et al., 1998; Bourguet et al., 2013; Bourguet and Triantafyllou, 2016). In this paper, however, this effect is neglected, and circular cylinders are considered instead, with a focus on large Reynolds numbers. Despite this limitation, the present study is relevant for the wind energy industry because of the interest of wind turbine developers in reducing the tapering at the top sections. The phenomenological aspects of VIVs for circular cylinders have been extensively analysed in the literature for several decades. However, studies involving large Reynolds numbers and large mass ratios (defined as the ratio between the structure mass and the mass of the displaced fluid) are still quite rare. Experiments in both air (Feng, 1968; Brika and Laneville, 1993) and water (Govardghan and Williamson, 2000) have shown that cylinders with large mass ratios exhibit much smaller amplitudes of vibrations than those with low mass ratios. Belloli et al. (2012), however, showed that this is not necessarily true in the presence of large Reynolds numbers of the order of Re=50 000. In this context, the maximum non-dimensional oscillation amplitude can become larger than the cylinder diameter. This was also described in other studies with low mass ratios (Govardghan and Williamson, 2006). For this reason, the original Griffin plot (Griffin, 1980), which describes the maximum oscillation amplitude as a function of the mass damping parameter, needs to be modified for large Reynolds numbers.

It is worth noting that all these VIV studies were limited to the subcritical regime, in which the boundary layer remains fully laminar, and the drag coefficient is nearly constant with the Reynolds number. By contrast, wind turbine towers experience flows both in the transitional regime ($\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{5}}\times \le \mathit{Re}\le \times \mathrm{3.5}\times {\mathrm{10}}^{\mathrm{6}}$), where laminar separation bubbles can exist and reduce the drag coefficient, and in the supercritical regime ($\mathit{Re}\ge \mathrm{3.6}\times {\mathrm{10}}^{\mathrm{6}}$), in which the boundary layer is fully turbulent. Experiments in these regimes can be challenging due to the need for large wind tunnels and wind speeds. In the context of stationary cylinders, a series of experiments have been performed with these very large Reynolds numbers (Roshko, 1961; Achenbach, 1968; Jones et al., 1969; van Nunen, 1974; Schewe, 1983). These studies generally agree well with the value of the separation angle and the location of the minimum pressure coefficient. However, large discrepancies are found between the various experimental data for the values of the aerodynamic quantities, such as pressure coefficient, lift and drag coefficients, and Strouhal number. Jones et al. (1969) also looked at a circular cylinder under forced vibration. Their study showed that a lift amplification exists and is due to the cylinder motion, which is in agreement with other works. Numerical studies have also been performed for flow past 2D and 3D stationary cylinders. These include unsteady Reynolds-averaged Navier–Stokes (URANS) models (Travin et al., 2000; Catalano et al., 2003; Ong et al., 2009; Rosetti et al., 2012), large-eddy simulations (LESs; Breuer, 2000; Catalano et al., 2003; Singh and Mittal, 2005), and detached-eddy simulations (DESs; Travin et al., 2000; Squires et al., 2008). Scatter has been observed between the various URANS simulations and was attributed to either different wall function implementations or the fact that some simulations were done in 2D instead of 3D. In most cases, the numerical results also deviate from the experiments, especially with URANS simulations. The latter methodology indeed presents shortcomings, such as an isotropic eddy viscosity, homogeneous Reynolds stresses, and the modelling of the full range of turbulent eddies (Rosetti et al., 2012). Yet, it was found that URANS simulations can still lead to satisfactory engineering results in the supercritical flow regime (Ong et al., 2009). The present study will go a step further by presenting URANS results of a freely vibrating cylinder in the supercritical regime ($\mathit{Re}\ge \mathrm{3.6}\times {\mathrm{10}}^{\mathrm{6}}$). To the best of the authors' knowledge, this is the first numerical work to do so.

The paper is organized as follows. Section 2 presents the fluid dynamics and structural dynamics models used in this study. It also explains the algorithm used to couple these models and perform two-way coupled fluid–structure interaction (FSI) simulations of a freely vibrating cylinder. Section 3 shows the results for two different regimes. First, laminar flow is considered in order to validate both the models and the FSI algorithm. The results of a freely vibrating cylinder are presented and compared to the literature. Second, the turbulent supercritical regime is considered for three cases: stationary cylinder, cylinder undergoing forced vibrations, and cylinder undergoing free vibrations. Finally, conclusions are drawn in Sect. 4.

2 Computational approach

2.1 Fluid dynamics model

The computational fluid dynamics (CFD) model used in this study is the open-source code Open-FOAM-v1812. The model uses a finite-volume discretization of the incompressible Navier–Stokes equations for a Newtonian fluid. In the laminar regime, these equations are solved directly, without using a turbulence model, through direct numerical simulation (DNS). In the turbulent regime, the unsteady Reynolds-averaged Navier–Stokes (URANS) equations are solved using a k–ω shear stress transport (SST) turbulence model (Menter, 1994), namely

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{u}}}{\partial t}}+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{\mathit{u}\mathit{u}}\right)=-{\displaystyle \frac{\mathrm{\nabla }\stackrel{\mathrm{‾}}{p}}{\mathit{\rho }}}+\mathrm{\nabla }\cdot \left(\mathit{\nu }\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{u}}\right)-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\mathbf{\prime }}{\mathit{u}}^{\mathbf{\prime }}}\right),\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{\mathit{u}}=\mathrm{0},\end{array}$$

where the fluid velocity is decomposed into an averaged and a fluctuating part, i.e. $\mathit{u}(\mathit{x},t)=\stackrel{\mathrm{‾}}{\mathit{u}}(\mathit{x},t)+{\mathit{u}}^{\mathbf{\prime }}(\mathit{x},t)$ with $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\mathbf{\prime }}}=\mathrm{0}$. Furthermore, p is the pressure field, ρ is the fluid density, and ν is the fluid kinematic viscosity. Additionally, the turbulence model assumes that

$$\begin{array}{}\text{(3)}& -{\stackrel{\mathrm{‾}}{u_i^{\prime }{u}^{\prime }}}_j=\mathrm{2}{\mathit{\nu }}_t{S}_{ij}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\delta }}_{ij}k,\end{array}$$

with δ_ij denoting the Kronecker symbol, ν_t the eddy viscosity, S_ij the strain rate tensor defined as

$$\begin{array}{}\text{(4)}& S_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_j}{\partial x_i}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\delta }}_{ij}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_k}{\partial x_k}},\end{array}$$

and k the turbulent kinetic energy defined as

$$\begin{array}{}\text{(5)}& k={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_i^{\prime }}.\end{array}$$

The specific implementation of the model is that presented by Menter et al. (2003), with a revised turbulence-specific dissipation rate production term from Menter and Esch (2001). The turbulent-specific dissipation rate ω is computed by solving the following equation:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\omega }}{\partial t}}& +{\displaystyle \frac{\partial u_j\mathit{\omega }}{\partial x_j}}=\mathit{\gamma }{\displaystyle \frac{{\stackrel{\mathrm{̃}}{P}}_k}{{\mathit{\nu }}_t}}-\mathit{\beta }{\mathit{\omega }}^{\mathrm{2}}+{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mathit{\nu }+{\mathit{\sigma }}_{\mathit{\omega }}{\mathit{\nu }}_t\right){\displaystyle \frac{\partial \mathit{\omega }}{\partial x_j}}\right]\\ & +\mathrm{2}(\mathrm{1}-F_{\mathrm{1}}){\mathit{\sigma }}_{{\mathit{\omega }}_{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{\displaystyle \frac{\partial k}{\partial x_j}}{\displaystyle \frac{\partial \mathit{\omega }}{\partial x_j}},\end{array}\end{array}$$

where the blending function F₁ is given by

$$\begin{array}{}\text{(7)}& F_{\mathrm{1}}=\tanh \left[{\left[min\left[max\left({\displaystyle \frac{\sqrt{k}}{{\mathit{\beta }}^{*}\mathit{\omega }\mathit{\gamma }}},{\displaystyle \frac{\mathrm{500}\mathit{\nu }}{{\mathit{\gamma }}^{\mathrm{2}}\mathit{\omega }}}\right),{\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{{\mathit{\omega }}_{\mathrm{2}}}k}{{\mathrm{CD}}_{k\mathit{\omega }}{\mathit{\gamma }}^{\mathrm{2}}}}\right]\right]}^{\mathrm{4}}\right],\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\mathrm{CD}}_{k\mathit{\omega }}=max\left(\mathrm{2}{\mathit{\sigma }}_{{\mathit{\omega }}_{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{\displaystyle \frac{\partial k}{\partial x_j}}{\displaystyle \frac{\partial w}{\partial x_j}},{\mathrm{10}}^{-\mathrm{10}}\right).\end{array}$$

The transport equation for the turbulent kinetic energy k is given by

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial k}{\partial t}}+{\displaystyle \frac{\partial u_{j}k}{\partial x_j}}={\stackrel{\mathrm{̃}}{P}}_k-{\mathit{\beta }}^{*}k\mathit{\omega }+{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mathit{\nu }+{\mathit{\sigma }}_k{\mathit{\nu }}_t\right){\displaystyle \frac{\partial k}{\partial x_j}}\right],\end{array}$$

where the limited production term ${\stackrel{\mathrm{̃}}{P}}_k$ is given by

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{̃}}{P}}_k=min(P_k,\mathrm{10}{\mathit{\beta }}^{*}k\mathit{\omega }),\end{array}$$

and

$$\begin{array}{}\text{(11)}& P_k={\mathit{\nu }}_t{\displaystyle \frac{\partial u_i}{\partial x_j}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right).\end{array}$$

Furthermore, the values of the turbulence coefficients are given in Table 1.

Table 1 Turbulence coefficients for the present k–ω SST model.

If needed, the constants are blended (typically close to the boundary layer) by the following interpolation:

$$\begin{array}{}\text{(12)}& \mathit{\varphi }=F_{\mathrm{1}}{\mathit{\varphi }}_{\mathrm{1}}+(\mathrm{1}-F_{\mathrm{1}}){\mathit{\varphi }}_{\mathrm{2}},\end{array}$$

where ϕ represents any of the coefficients in Table 1. Once the two turbulence transport equations are solved, the eddy viscosity field is obtained as

$$\begin{array}{}\text{(13)}& {\mathit{\nu }}_t={\displaystyle \frac{a_{\mathrm{1}}k}{max(a_{\mathrm{1}}\mathit{\omega },\mathrm{\Omega }{F}_{\mathrm{2}})}},\end{array}$$

in which Ω is the magnitude of the strain rate tensor, a₁ is a coefficient defined in Table 1, and F₂ is the following blending function:

$$\begin{array}{}\text{(14)}& F_{\mathrm{2}}=\tanh \left[{\left[max\left({\displaystyle \frac{\mathrm{2}\sqrt{k}}{{\mathit{\beta }}^{*}\mathit{\omega }\mathit{\gamma }}},{\displaystyle \frac{\mathrm{500}\mathit{\nu }}{{\mathit{\gamma }}^{\mathrm{2}}\mathit{\omega }}}\right)\right]}^{\mathrm{2}}\right].\end{array}$$

The equations are discretized using second-order methods in both space and time. Only the equations for k and ω are solved using a first-order upwind scheme. The PISO-SIMPLE (PIMPLE) algorithm is used for the pressure–velocity coupling, where PISO is the pressure-implicit splitting of operators, and SIMPLE is the semi-implicit method for pressure-linked equations. This model represents the state of the art in URANS modelling. However, it does have the limitation that the model is two-dimensional and therefore cannot take vertical wind shear into account. Furthermore, it has some limitations for the boundary layer transition regime.

2.2 Structural dynamics model

The two-dimensional cylinder is modelled as a rigid body of mass m attached to a linear spring with spring constant k_s and a viscous damper with damping coefficient c. The motion of the rigid body is constrained so that it cannot rotate and can only translate in the y direction, transverse to the inflow. The response of the body is therefore governed by the equations of motion of a single degree of freedom oscillator, i.e.

$$\begin{array}{}\text{(15)}& m\ddot{y}+c\dot{y}+k_{\mathrm{s}}y=F_y,\end{array}$$

where $\ddot{y}$, $\dot{y}$, and y are the y components of the linear acceleration, velocity, and displacement, respectively, of the centre of gravity of the structure, and F_y is the resultant traction force exerted on the structure by the fluid. Equation (15) can be written as

$$\begin{array}{}\text{(16)}& \ddot{y}+\mathrm{2}\mathit{\zeta }{\mathit{\omega }}_{\mathrm{n}}\dot{y}+{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}y={\displaystyle \frac{F_y}{m}},\end{array}$$

where ${\mathit{\omega }}_{\mathrm{n}}=\sqrt{k_{\mathrm{s}}/m}$ is the natural frequency of the structure, and ζ is the damping coefficient defined as the ratio between the actual damping and the critical value, i.e.

$$\begin{array}{}\text{(17)}& \mathit{\zeta }={\displaystyle \frac{c}{\mathrm{2}\sqrt{k_{\mathrm{s}}m}}}.\end{array}$$

For a three-dimensional tower, the values for mass, stiffness, and damping vary along the tower. Here, a modal analysis was performed to lump these values into a single value representing the two-dimensional property of the cylinder. To do so, the three-dimensional tower is divided into 40 segments with 6 degrees of freedom at each node, leading to stiffness and mass matrices of dimensions 246×246. The natural frequency of the first mode is obtained by solving Eq. (15) in the frequency domain. The modal mass and stiffness are then obtained by premultiplying and postmultiplying the mass and stiffness matrices by the eigenvector of the first bending mode. Equation (16) is further non-dimensionalized as

$$\begin{array}{}\text{(18)}& {\ddot{y}}^{*}+\mathrm{2}\mathit{\zeta }{\dot{y}}^{*}+y^{*}=\mathrm{2}n{U}^{*}{C}_{F_{y^{*}}},\end{array}$$

denoting $y^{*}=y/D$ (D being the cylinder diameter), $t^{*}=t/{\mathit{\omega }}_{\mathrm{n}}$, and $m^{*}=\mathrm{2}m/\left(\mathit{\rho }D\right)$. Additionally, $C_{F_{y^{*}}}$ stands for the aerodynamic coefficient expressed in terms of y^*, $n=\mathit{\rho }{D}^{\mathrm{2}}/\left(\mathrm{4}m\right)$, and the reduced velocity is $U^{*}=U/\left(D{\mathit{\omega }}_{\mathrm{n}}\right)$.

2.3 Fluid–structure interaction coupling algorithm

In this work, fluid–structure interactions are modelled using a partitioned approach, whereby the fluid and structural dynamics are solved alternatively at every time step. Here, both weakly and strongly coupled schemes are considered. This is illustrated by Fig. 1. At the start of a new time step, Eq. (15) is solved for the structural dynamics. Once the structure displacement is known, the geometry is moved accordingly, and the mesh is diffused using a spherical linear interpolation scheme (SLERP) algorithm. With this approach, it is possible to specify a mesh region where the cells preserve their shape during motion. In this work, for the freely vibrating cases, the mesh at a distance of up to 25D around the cylinder is kept rigid throughout the simulations in order to keep its initial high-quality characteristics in the O-grid block. By contrast, the mesh deformation is applied to all cells located farther away from the cylinder. Once the mesh has moved, the fluid dynamics solver solves the Navier–Stokes equations, accounting for the mesh motion using an arbitrary Lagrangian–Eulerian approach. In the weakly coupled algorithm, the fluid dynamics equations are solved iteratively without solving again for the structural dynamics. The structural and mesh motion are then only solved once per time step. By contrast, in the strongly coupled approach, the structural dynamics equations are solved for each sub-iteration of the fluid solver. Thus, the structural dynamics and mesh movement are solved multiple times per time step.

Figure 1 Sketch of the weakly (explicit) and strongly (implicit) coupled algorithms between structure and fluid solvers. Figure taken from Derksen (2019).

In this work, a weakly coupled algorithm was found to give accurate results for modelling VIVs in the laminar regime. In that case, the structural dynamics equations were also solved using an explicit solver, the so-called symplectic integrator (Dullweber et al., 1997), which is based on the leapfrog method. The explicit character of this solver leads to a constant structural displacement for each iteration within a time step, while the acceleration and velocity may change at each iteration. The structural displacement of the current time step is only dependent on the acceleration computed in the previous time step. In the turbulent regime, a strongly coupled FSI scheme was adopted, and three to four sub-iterations were necessary between fluid and structural solvers to achieve accurate results. In that case, an implicit Newmark structural solver was further used with the so-called average constant acceleration method (Newmark, 1959).

3 Results

The results are shown for two distinct flow regimes. First, a freely vibrating cylinder is considered in a laminar flow in order to validate both the set-up and the fluid–structure interaction algorithm. Second, the results are presented in a turbulent regime relevant for wind turbines. In that case, stationary cylinders, forced vibrations, and free vibrations are considered. The data presented in this section are available at https://doi.org/10.4121/uuid:33f213fe-df8e-4b71-bec5-feada4e12a6d (Derksen et al., 2020).

3.1 Laminar flow

This section presents the simulation set-up and results for a freely vibrating cylinder in a laminar flow.

3.1.1 Simulation set-up

The computational domain used in the present work is sketched in Fig. 2. The cylinder is centred in a squared domain of length 50D, D being the cylinder diameter. A uniform streamwise velocity field is imposed as an initial condition and also as a Dirichlet boundary condition at the inlet of the domain. A zero-pressure Dirichlet condition is set at the outlet. The lateral far-field boundaries are set as slip walls. In this subsection, the free-vibration results are shown, whereby the dynamics of the cylinder are determined by solving Eq. (15). The domain is meshed with a structured hexahedral grid with a minimum cell height at the cylinder corresponding to $y^{+}=y{u}_{\mathit{\tau }}/\mathit{\nu }=\mathrm{1}$, where $u_{\mathit{\tau }}=\sqrt{{\mathit{\tau }}_w/\mathit{\rho }}$ is the friction velocity at the wall. As shown by Fig. 3, the mesh has an O-grid topology, which is commonly used for flow around cylinders. For the present laminar flow computations, the cell height at the cylinder boundary of Δy is 0.05. A mesh convergence analysis was performed and led to the conclusions that a mesh with 30 096 cells was sufficient to obtain accurate converged results. The time step size was further varied to keep the Courant–Friedrichs–Lewy (CFL) number equal to 0.7 for all the computations.

Figure 2 Sketch of the computational domain.

Figure 3 Computational mesh: full domain (a), zoom around the cylinder (b).

3.1.2 Free vibration

This section shows the results of a cylinder vibrating freely in the direction transverse to the flow. Here, the flow is laminar, and the results are compared with those available in the literature (Carmo et al., 2011), where a loosely coupled numerical FSI framework was also used. The simulations are conducted with a Reynolds number of Re=150, a mass ratio $m^{*}=\mathrm{4}m/\left(\mathit{\rho }\mathit{\pi }{D}^{\mathrm{2}}L\right)=\mathrm{2}$ (L being the length of the equivalent three-dimensional cylinder), and a damping factor of ζ=0.007. The reduced velocity $U^{*}=U/\left({\mathit{\omega }}_{\mathrm{n}}D\right)$ varies in the range of $\mathrm{2.5}<U^{*}<\mathrm{16}$. The non-dimensional maximum amplitude of the cylinder is shown in Fig. 4 for both the present simulations and the results from the literature. The lock-in region is clearly identified for $\mathrm{4}<U^{*}<\mathrm{7}$, and good agreement between the present results and the reference solution is obtained for the whole range of reduced velocities. Figure 5 shows the cylinder frequency divided by its natural frequency. Again, good agreement is obtained between the present results and the literature. As expected, the motion frequency equals the natural frequency in the lock-in region (i.e. $\mathit{\omega }/{\mathit{\omega }}_n=\mathrm{1}$), whilst the frequency ratio follows the Strouhal relation outside the lock-in band. These results demonstrate that the present weakly coupled FSI model is capable of predicting the dynamics of a light cylinder undergoing VIV in laminar flow conditions.

Figure 4 Non-dimensional maximum motion amplitude as a function of the reduced velocity for Re=150, $m^{*}=\mathrm{2}$, and ζ=0.007.

Figure 5 Ratio of the motion frequency to the natural frequency as a function of the reduced velocity for Re=150, $m^{*}=\mathrm{2}$, and ζ=0.007.

3.2 Turbulent flow

In this section, the cylinder is placed in a turbulent flow as typically encountered by wind turbine towers. First, the simulation set-up is presented. Then, the simulation results are presented for three conditions: stationary cylinder, forced vibrations, and free vibrations.

3.2.1 Simulation set-up

The computational domain and boundary conditions are identical to those used in the previous subsection. However, since a turbulence model is used here, additional boundary and initial conditions are needed for the turbulent quantities. In particular, the initial value of the turbulent kinetic energy is set to

$$\begin{array}{}\text{(19)}& k_{t=\mathrm{0}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\left(\mathrm{TI}{u}_{\mathrm{\infty }}\right)}^{\mathrm{2}},\end{array}$$

with the turbulence intensity being set at TI =0.03, and the inlet velocity is u_∞=1, whilst the initial value of the turbulence-specific dissipation rate is equal to

$$\begin{array}{}\text{(20)}& {\mathit{\omega }}_{t=\mathrm{0}}={\displaystyle \frac{\mathit{\rho }k}{\mathit{\mu }}}{\left({\displaystyle \frac{{\mathit{\mu }}_t}{\mathit{\mu }}}\right)}^{-\mathrm{1}},\end{array}$$

with ${\mathit{\mu }}_t/\mathit{\mu }=\mathrm{10}$, μ being the fluid dynamic viscosity. The initial value of the eddy viscosity is ν_t=0. Regarding the boundary conditions, Dirichlet conditions are imposed on both k and ω at the inlet, with values equal to the initial ones, respectively, and Neumann conditions are imposed at the outlet. At the cylinder, the mesh is fine enough to resolve the flow up to the viscous sub-layer ($y^{+}=\mathrm{1}$) using prism layers around the cylinder, with a mesh growth ratio equal to 1.2 and a minimum cell height at the cylinder equal to $\mathrm{\Delta }y=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{6}}$. Converged results were obtained for a mesh with 98 730 grid points. As for the laminar flow simulations, the maximum CFL number was kept constant at CFL =0.7. Since the viscous sub-layer is mostly laminar, the turbulent quantities are such that turbulence is suppressed at the wall. Thus, at the cylinder boundary, $k={\mathrm{10}}^{-\mathrm{10}}{m}^{\mathrm{2}}/s^{\mathrm{2}}$, ν_t=0, and

$$\begin{array}{}\text{(21)}& \mathit{\omega }={\displaystyle \frac{\mathrm{6}\mathit{\mu }}{{\mathit{\beta }}_{\mathrm{1}}{y}_{\mathrm{wall}}^{\mathrm{2}}}},\end{array}$$

with β₁=0.075 and y_wall being the mesh cell height at the wall (Menter, 1992).

3.2.2 Flow past a stationary cylinder

Before considering flow past moving cylinders, the accuracy of the k–ω SST turbulence model in the supercritical regime is first assessed on flow past a stationary cylinder with a Reynolds number $\mathit{Re}=\mathrm{3.6}\cdot {\mathrm{10}}^{\mathrm{6}}$. Table 2 compares the results of the present simulations with numerical and experimental data from the literature in terms of the Strouhal number St, mean drag coefficient C_D,mean, base pressure coefficient c_p,b at $\mathit{\theta }=\mathrm{180}{}^{\circ }$ (θ being the circumferential angle at a point on the cylinder surface starting from the stagnation point), and flow separation angle θ_sep. The values marked with a double asterisk are estimated indirectly based on the available pressure distribution where the constant pressure plateau is reached. Overall, it is observed that the present results agree well with the literature, especially when compared with other numerical results obtained with either URANS (Ong et al., 2009) or DES (Travin et al., 2000) turbulence models. The pressure coefficient at the cylinder surface is further plotted in Fig. 6. Again, it is clear that the present URANS results match well with the results from other numerical studies. Compared to experimental works (Achenbach, 1968; Jones et al., 1969), the present pressure coefficient agrees well in the front part of the cylinder for $\mathrm{0}{}^{\circ }\le \mathit{\theta }\le \mathrm{63}{}^{\circ }$. For larger values of θ, the numerical results show both an earlier pressure recovery (at $\mathit{\theta }\approx \mathrm{80}{}^{\circ }$) and a less pronounced negative peak for c_p. It was shown by Achenbach (1968) that the transition point of the boundary layer is at $\mathit{\theta }=\mathrm{65}{}^{\circ }$ for that Reynolds number. It seems that the pressure distribution starts to deviate around this point, which could indicate that the transition has an impact on the pressure distribution. The boundary layer in the experiment of Jones et al. (1969) continues to decelerate further aft (at $\mathit{\theta }=\mathrm{90}{}^{\circ }$), which explains why the separation also occurs later in their measurements. This might indicate that the numerical models do not fully predict the detached flow accurately. This is in line with previous observations made in the literature, although with different Reynolds numbers (Ferrand et al., 2006; Catalano et al., 2003). The friction coefficient at the cylinder surface is also plotted in Fig. 7. The present numerical results agree well with other numerical studies that use either LES (Catalano et al., 2003) or DES (Travin et al., 2000). The peak of the maximum friction coefficient is slightly lower than that obtained from DES. However, the overall agreement is good, and the separation location is also accurately captured and compares well with the experiment of Achenbach (1968). The numerical values of the friction coefficient c_f for $\mathrm{0}{}^{\circ }\le \mathit{\theta }\le \mathrm{90}{}^{\circ }$ are generally much higher than the experimental values. This can be explained by the fact that the numerical studies assume fully turbulent transition and therefore lead to significantly higher friction coefficient values. In particular, the overestimation is caused by the fact that the laminar flow at the beginning part of the cylinder wall is not taken into account in the numerical models. This is in line with the observations of Travin et al. (2000) and Catalano et al. (2003). Only the URANS simulations of Ong et al. (2009) show much smaller friction coefficient values. However, this study used a k–ϵ turbulence model, which is known to perform poorly with flow separation and strong pressure gradients. This could explain why these URANS results underestimate the skin friction coefficient compared to the other numerical studies. Despite the discrepancies in c_p and c_f between the numerical and experimental results, the present turbulence model is deemed to be adequate for the purpose of this work as it provides results that are in line with other trustworthy numerical studies in the supercritical regime.

(Ong et al., 2009)(Travin et al., 2000)(Achenbach, 1968)(Jones et al., 1969)

Table 2 Comparison of data for flow past a stationary cylinder in the supercritical regime. The values marked with ^* have been estimated indirectly based on the available pressure distribution where the constant pressure plateau was reached.

Figure 6 Time-averaged pressure distribution for flow past a stationary cylinder in the supercritical turbulent regime.

Figure 7 Time-averaged friction coefficient for flow past a stationary cylinder in the supercritical turbulent regime.

3.2.3 Flow past a cylinder under forced vibration

In this section, the transverse motion of the cylinder is prescribed for different frequencies ω_f. Seven different frequencies are considered, each corresponding to a different reduced velocity and Reynolds number, as indicated in Table 3. In the latter, St_f denotes the non-dimensional forcing frequency of the cylinder, while St is the Strouhal number associated with flow past a stationary cylinder with the same Reynolds number. The latter is obtained numerically with the same mesh as for the corresponding moving cases. The simulations are further run for seven different values of the prescribed maximum motion amplitude. The root mean square of the lift coefficient (non-dimensionalized by its value for a stationary cylinder) is shown in Fig. 8 for different ratios St_f∕St and prescribed motion amplitudes $y_{\max }^{*}$. As expected, the lift magnification due to the cylinder motion increases with both the non-dimensional frequencies and $y_{\max }^{*}$. The increase with frequency is rather linear except close to ${\mathit{St}}_{\mathrm{f}}/\mathit{St}=\mathrm{1}$. The peak in lift magnification around a unitary reduced frequency is believed to be a consequence of the wake lock-in on the cylinder motion. For ${\mathit{St}}_{\mathrm{f}}/\mathit{St}>\mathrm{1}$ and $y_{\max }^{*}\ge \mathrm{1}$, our simulation results indicate that the lift magnification becomes very large, which is believed to be due to the effective added mass of the fluid caused by the vorticity dynamics (Williamson and Govardhan, 2008).

Figure 8 Root mean square of the lift coefficient (non-dimensionalized by its value for a stationary cylinder) for forced vibrations in the supercritical turbulent regime with different frequency ratios St_f∕St and prescribed motion amplitudes $y_{\max }^{*}$.

Another observation is that the lock-in band width increases with increasing motion amplitudes, with a move towards lower reduced frequencies as the oscillation amplitude increases. Both these observations, namely increased lift magnification at high frequency ratios and wider lock-in band with larger motions, were observed in the experimental results obtained from NASA under similar conditions (Jones et al., 1969), as shown in Fig. 9. Although the experimental results show some scatter, a peak of lift magnification is obtained around ${\mathit{St}}_{\mathrm{f}}/\mathit{St}=\mathrm{1}$, decreasing for slightly larger Strouhal number ratios before increasing further as ${\mathit{St}}_{\mathrm{f}}/\mathit{St}>\mathrm{1}$ (see for example the results at $y_{\max }^{*}=\mathrm{0.0278}$ in Fig. 9).

Table 3 Range of frequency ratios and velocities considered for the forced vibration of a cylinder in the supercritical regime.

Figure 9 Root mean square of the lift coefficient (non-dimensionalized by its value for a stationary cylinder) for forced vibrations compared with NASA experimental data (Jones et al., 1969) at three values of the prescribed motion amplitude.

Figure 10 Aerodynamic damping for forced vibrations in the supercritical turbulent regime with different frequency ratios St_f∕St and prescribed motion amplitudes $y_{\max }^{*}$.

Although the cylinder motion is prescribed, it is interesting to look in more detail at the feedback mechanism between the cylinder and the fluid motions. In order to do this, the aerodynamic damping is analysed. It is defined as the part of the harmonic lift force that is in phase with the cylinder velocity. There are different ways to compute this damping. Here, an approach based on the time-averaged energy transfer between fluid and vibrating cylinder is considered (Bourguet et al., 2011; Gopalakrishnan, 1993) as it is believed to be accurate when signals exhibit multiple frequencies. This leads to the following expression for the aerodynamic damping:

$$\begin{array}{}\text{(22)}& C_{L,V}=-\sqrt{\mathrm{2}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{C_L^{\prime }{\dot{y}}^{\prime }}}{\sqrt{\stackrel{\mathrm{‾}}{{\dot{y}}^{\mathrm{2}}}}}},\end{array}$$

where $C_L=\mathrm{2}{F}_y/\left(\mathit{\rho }{U}^{\mathrm{2}}D\right)$ denotes the lift coefficient, a prime denotes the fluctuating part of a quantity, and an overline denotes its time average. The aerodynamic damping is plotted for the different simulation cases in Fig. 10. For ${\mathit{St}}_{\mathrm{f}}/\mathit{St}\ge \mathrm{1}$, it is observed that C_L,V increases with the motion amplitude. This indicates that the wake dynamics tend to oppose the cylinder motion as the latter increases, hence stabilizing the system. By contrast, for ${\mathit{St}}_{\mathrm{f}}/\mathit{St}\le \mathrm{0.75}$, the aerodynamic damping is negative, with an absolute value increasing as the cylinder displacement increases. In that case, the wake dynamics tend to further amplify the system motion, hence leading to an unstable system behaviour. In the lock-in band, for $\mathrm{0.8}\le {\mathit{St}}_{\mathrm{f}}/\mathit{St}\le \mathrm{1}$, the nature of the interaction between fluid and structural motion depends on the amplitude of the prescribed motion. However, a stabilizing fluid–structure interaction behaviour is observed for large motion amplitudes. These observations agree with those obtained from Jones et al. (1969). In particular, the overall behaviour of the experimental aerodynamic damping in the supercritical turbulent regime for varying frequency ratios was found to be similar to the present CFD results.

Figure 11 Non-dimensional maximum motion amplitude as a function of the reduced velocity for $m^{*}=\mathrm{29.6}$, ζ=0.003, and $\mathrm{1.88}\le U^{*}\le \mathrm{10.74}$. Insets show the time evolution of the lift coefficient at certain values of reduced velocities.

3.2.4 Flow past a freely vibrating cylinder

In this section, the cylinder is left to freely vibrate under the vortex-induced vibrations. The results are shown for a baseline value of mass ratio $m^{*}=\mathrm{29.6}$ and damping coefficient ζ=0.003, which are representative of a realistic wind turbine tower. The effect of these values on the results is also briefly discussed in this section. The reduced velocity is varied from $U^{*}=\mathrm{1.88}$ to $U^{*}=\mathrm{10.74}$, corresponding to Reynolds numbers between $\mathit{Re}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{6}}$ and $\mathit{Re}=\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{7}}$, respectively. The associated non-dimensional maximum motion amplitude under these conditions is shown in Fig. 11, in which the lock-in band can be clearly visualized. Additionally, the insets show the time evolution of the lift coefficient on the cylinder for certain values of reduced velocity. The associated values of reduced frequency, ${\mathit{\omega }}^{*}={\mathit{\omega }}_{\mathrm{n}}/{\mathit{\omega }}_{(\mathrm{shed},\mathrm{stat})}$, are also shown for each inset. These insets illustrate that for some values of U^*, a non-harmonic lift response (and consequently also motion response) is obtained. Outside the lock-in region, at low values of U^*, the vortex-shedding frequency is dominated by the stationary Strouhal relation and is very different from the natural frequency of the cylinder. This is also apparent from Fig. 12. In that case, the mean aerodynamic forces are similar to those of a stationary cylinder, and little cylinder motion is observed. Furthermore, the vorticity contours in the cylinder wake follow a 2S pattern similar to flow past a stationary cylinder. Furthermore, the phase angle ϕ between the lift force and the transverse displacement is constant and equals 0^∘. When U^* increases to values close to the lock-in band, the time histories of the aerodynamic forces and cylinder motion change. First, the cylinder displacement increases substantially because the vortex-shedding frequency gets closer to the natural frequency of the cylinder. Second, the cylinder motion shows two main frequency responses instead of one. The dominant frequency in the time evolutions of the lift and motion amplitude still equals the Strouhal frequency. However, performing a power spectral density on these time signals shows that a second (much smaller) frequency peak coincides with the natural frequency. This explains the observed increase in cylinder motion compared to cases at smaller reduced velocities. The phase angle between lift force and transverse displacement is also different than for smaller values of U^*. Instead of being constant, it periodically oscillates in time around a mean value of $\mathit{\varphi }=\mathrm{5}{}^{\circ }$. When the phase angle is larger than the mean value, the aerodynamic damping decreases, and the fluid amplifies the cylinder excitation. By contrast, when the phase angle is smaller than the mean value, the fluid damps the cylinder motion. When U^* is increased further inside the lock-in band, the time histories of the lift and displacement also show two frequencies. However, the dominant peak corresponds to the natural frequency of the cylinder, hence confirming that the wake behaviour becomes driven by the cylinder oscillations rather than the Strouhal relation. Furthermore, the lift force (and also motion displacement) alternates between growth and decay, which also corresponds to a change in sign of the aerodynamic damping. In particular, when the lift force increases, the fluid dynamics excite the cylinder. By contrast, when the lift force decreases, the fluid dynamics damp the cylinder motion. In doing so, the vortex-induced vibrations continuously alternate between self-exciting and self-limiting behaviours. This alternation can also be related to a drop in phase angle between lift force and cylinder displacement. Figure 13 shows the time evolution of the phase angle for $U^{*}=\mathrm{3.1}$. It is apparent that the phase angle changes sign as time evolves. When ϕ>0, the fluid excites the cylinder motion, and the magnitude of both the lift force and the transverse displacement increases. By contrast, when ϕ<0, the lift force lags behind the displacement, and the fluid dampens the cylinder motion. The magnitude of both the lift force and the cylinder displacement also decreases. The phase angle drop can also be related to the wake pattern development as shown in Fig. 14. The left figure corresponds to a positive phase angle at $t/T_s=\mathrm{143.5}$, while the right figure is taken at $t/T_s=\mathrm{146.5}$, when the phase angle is negative. When the phase angle is positive, the vortices are shed from the bottom side of the cylinder when the cylinder reaches its maximum positive displacement. The opposite is true when the phase angle switches sign. In that case, the vortices are shed from the same side as where the cylinder is oscillating to. It seems that when the cylinder motion becomes too large, the wake reorganizes itself in such a way that it starts damping the cylinder motion. The self-exciting behaviour becomes a self-limiting behaviour. This alternation of behaviours is influenced by both the mass ratio and the damping factor. If the mass ratio is reduced, the non-dimensional displacement amplitude increases. This is expected as a lower mass ratio leads to a lower structural inertia when compared to the fluid inertia. This makes the cylinder more susceptible to oscillations and possibly lock-in. By contrast, for very large mass ratios, the cylinder motion decreases and eventually tends to a stationary cylinder. There is thus a critical value of the mass ratio at which the system stops undergoing enhanced fluid–structure interactions. Non-harmonic behaviours of the lift and motion displacement are also less likely to occur at large mass ratios and large damping factors because the larger structural inertia prevents the cylinder from changing to a different dynamic state. When U^* is increased further, in a certain range, the lift force, drag force, and transverse displacement do not show a clear converged pattern (see inset in Fig. 11 at $U^{*}=\mathrm{4.3}$). It is not excluded that, if the simulation were run for even longer times, a self-limiting behaviour could appear. However, the observed behaviour is believed to have physical causes rather than being a numerical artefact. This is because the operation points for which this happens correspond to reduced frequencies between 0.5 and 0.7, for which the results of the forced-vibration cases led to unstable behaviours. By contrast, at smaller values of U^*, the reduced frequencies are large ($\mathrm{0.9}\le {\mathit{\omega }}^{*}\le \mathrm{1.2}$), and a stable forced-vibration behaviour was obtained. It is also worth noting that the large lock-in band observed in the present study is in line with previous observations made in the literature under similar conditions. For example, Guilmineau and Queutey (2004) and Assi (2009) showed large values of non-dimensional displacement for reduced frequencies as low as ${\mathit{\omega }}^{*}=\mathrm{0.7}$, whilst the wind tunnel experiment of Feng (1968) also showed a lock-in region extending to reduced frequencies of ${\mathit{\omega }}^{*}=\mathrm{0.6}$.

Figure 12 Ratio of the motion frequency to the natural frequency as a function of the reduced velocity for $m^{*}=\mathrm{29.6}$, ζ=0.003, and $\mathrm{1.88}\le U^{*}\le \mathrm{10.74}$.

Figure 13 Time evolution of the phase angle between lift force and transverse displacement at $U^{*}=\mathrm{3.1}$. Figure taken from Derksen (2019).

Figure 14 Instantaneous non-dimensional vorticity magnitude for $U^{*}=\mathrm{3.1}$ at the maximum positive displacement: $t/T_s=\mathrm{143.5}$ (left) and $t/T_s=\mathrm{146.5}$ (right). The vorticity is filtered to the range $\mathrm{0}\le \mathit{\xi }D/U_{\mathrm{\infty }}\le \mathrm{5}$. Figures taken from Derksen (2019).

4 Conclusions

In this work, numerical simulations of a two-dimensional cylinder undergoing vortex-induced vibrations were performed. First, the fluid–structure interaction methodology was validated in the laminar regime, for which direct numerical simulations (DNSs) were able to reproduce results from the literature. Second, supercritical turbulent regimes experienced by wind turbine towers were simulated using the k–ω SST unsteady Reynolds-averaged Navier–Stokes (URANS) model. The ability of our model to simulate flow past a stationary cylinder in the turbulent supercritical regime was demonstrated by comparing the present results with other URANS, large-eddy simulations (LESs), detached-eddy simulations (DESs), and experimental results from the literature. Additionally, the present tools were capable of simulating the dynamic behaviour of the cylinder (including lock-in) for a range of reduced velocities and Reynolds numbers. In particular, good agreement was found between the present results and those obtained from wind tunnel experiments under forced vibrations. When the cylinder was left free to oscillate under the effect of vortex shedding, the results highlighted a complex interplay between structural and fluid dynamics for values of reduced velocity close to or inside the lock-in band. In particular, there was a continuous alternation between self-exciting and self-limiting vortex-induced vibrations. The physical feasibility of this interplay was supported by both the role of the aerodynamic damping, which was shown to continuously change sign at these operating points, and the results from forced-vibration oscillations under similar conditions.

Traditionally, the designs of wind turbine towers are soft-stiff, with the tower's natural frequency being larger than the rotor's rotational frequency but smaller than the blade's passing frequency. Today, advanced wind turbine controllers can work through resonance conditions, allowing the tower's natural frequencies to overlap with the rotor's rotational frequency and the blade's passing frequency. This makes soft-soft designs possible, which is promising to enable a significant reduction in tower mass when the tower height is larger than 100 m (Dykes et al., 2018). During the installation phase of wind turbines, however, the rotor–nacelle assembly is absent, and wind turbine controllers cannot be used to damp vibrations. This is the case during the whole installation cycle of jackets and monopile foundations as well as their vertical storage and transport. During that time, towers might stand without rotors for several weeks. They act as a clamped beam and are subjected to supercritical Reynolds numbers. Over the past years, the tower-clamped first-mode frequency has been reduced, and this trend is continuing. This means that both the combined first eigenmodes and the second mode of vibration can interact with the wind climate and become critical for the dynamics of the structure and VIVs. The present work shows that, in the supercritical flow regime, a cylinder can experience episodes of sustained VIV, where the Strouhal relation is not valid. This means that the best mitigation strategy is to have a robust damping mechanism that is effective for the range of critical first bending frequencies. This work presents some limitations that should be addressed in future studies. For example, future work should look at three-dimensional tower sections as well as the influence of non-uniform wind flows on the present results. Although this has already been partly investigated in the literature, there is a lack of analyses for the supercritical regime. Furthermore, the use of a tapered cylinder instead of a circular one is expected to lower the risk of VIVs and move the lock-in region to lower reduced velocities. Nevertheless, circular cylinders are very relevant to new wind turbine tower designs as some manufacturers are currently looking into having minimal taper at the top sections. Finally, during the installation phase, wind turbine towers are usually transported in groups. Therefore, wake–tower interactions and their effect on the development of VIVs should be further investigated.