Force-partitioning analysis of vortex-induced vibrations of wind turbine tower sections

Vortex-induced vibration (VIV) of wind turbine towers during installation is an aero-structural problem of significant practical relevance. Vibrations may happen in the tower structure, especially when the rotor-nacelle assembly is not yet attached to the tower or if the rotor blades are not yet connected to the tower-nacelle assembly. The complexity of aeroelastic phenomena involved in VIV makes modelling and analysis challenging. Therefore, the aim of the current research is to investigate the fundamental mechanisms causing the onset and sustenance of vortex-induced vibrations. To gain more understanding of the nature of vibrations, a methodology is established that distinguishes between different components of the forces at play. This approach allows for identifying how various force components impact the oscillation of a rigid body. The method is executed using the OpenFOAM open-source software. Numerical simulations are conducted on a two-dimensional smooth cylinder at both subcritical and supercritical Reynolds numbers to establish a correlation between wind turbine tower vibrations and the force mechanism. The analysis involves performing unsteady Reynolds-averaged Navier–Stokes (URANS) simulations using the modified pimpleFoam solver with the k–ω shear stress transport (SST) turbulence model. Both fixed and free-vibrating cases are studied for smooth cylinders. For the high-Reynolds-number cases, a setup matching the tower top segment of the IEA 15 MW reference wind turbine was chosen. Studying the flow around a cylinder at a subcritical Reynolds number reveals that the primary force involved is the vortex-induced force. The combined force due to viscosity, added mass, and vorticity contributes most to the overall force. For a freely vibrating cylinder with a single degree of freedom in the crossflow direction, the analysis indicates that the force component associated with the cylinder's motion is crucial and significantly affects the total force. Moreover, analysing the energy transfer between the fluid and the structure, a positive energy contribution by the vortex-induced force is observed on or before the dominant Strouhal velocity. This confirms observations at low Reynolds numbers in the literature that the vortex shedding predominantly contributes to the initiation of oscillations during VIV. The kinematic force contributes to the energy transfer of the system, but the mean energy transfer per cycle is negligible.

1 Introduction

Wind turbines continue to grow in size for energy production at a lower cost of energy. One of the challenges that follows this trend is the scaling of wind turbines, as taller wind turbines have longer and more flexible blades and towers (Hoen et al., 2023; Veers et al., 2019). Many physical phenomena must be taken into consideration during the design process for such large turbines. Vortex-induced vibration (VIV) is one such aerodynamic phenomenon that must be understood in detail during the design of wind turbines. Both wind turbine blades and towers can suffer from VIV, for which many studies have been carried out in the past (Livanos, 2018; Derksen, 2019; Horcas et al., 2020; Riva et al., 2022). However, the complex aerodynamic interactions that give rise to VIV are not fully understood, especially for the complex flow conditions experienced by wind turbines. The towers are prone to oscillations or vibrations due to vortex shedding, particularly when the rotor-nacelle assembly or the rotor blades are not yet installed. The load cases in which the wind turbine towers are susceptible to VIV are when the towers are (i) standing on the quayside, (ii) transported in a vessel, and (iii) standing on the foundation (e.g. monopile or jacket for offshore turbines) either without the rotor-nacelle assembly or before the installation of the blades on the tower-nacelle assembly. Figure 1 shows a load case when a tower is standing on the foundation during installation at the site. Understanding VIV and its impact on wind turbine towers is crucial for the design and operation of these structures. Therefore, further studies are needed to analyse VIV and its mitigation for wind energy systems.

Figure 1 Example of a load case where a wind turbine tower is susceptible to VIV, with the tower standing on the foundation during installation.

Extensive research has been conducted in the past to investigate VIV of circular cylinders, particularly for subcritical and critical Reynolds numbers (Williamson and Roshko, 1988; Bishop and Hassan, 1964; Brika and Laneville, 1993; Williamson and Govardhan, 2004). These studies investigated the vortex patterns and associated forces at various Reynolds numbers and the vibration pattern with respect to amplitude and direction. However, VIV is a complex phenomenon that depends on many other structural and fluid properties, like mass, natural frequency, damping ratio, inflow velocity, turbulence intensity, viscosity, etc. Whilst most of these studies have been carried out for low-Reynolds-number regimes, significantly fewer studies have been done for super- and transcritical Reynolds numbers. This is mainly due to the difficulties of performing experiments at these high Reynolds numbers and the high computational expense of numerical simulations. Experiments carried out by Belloli et al. (2012, 2015) analyse the lift force behaviour and the wakes formed behind the cylinder. The cylinder is covered with a rough mesh to increase roughness and an effective Reynolds number to perform these experimental analyses. Upon further investigation, the motion experienced at a high mass ratio (m^*) and low structural damping shows that the vibration amplitudes are much higher than those according to the Griffin plot (Griffin, 1980).

In addition to the experimental studies, several numerical studies were carried out at transcritical Reynolds numbers. These studies were performed using large-eddy simulations (LESs), unsteady Reynolds-averaged Navier–Stokes (URANS) simulations, and a hybrid RANS–LES approach (Catalano et al., 2003; Squires et al., 2008; Travin et al., 2000). These studies analysed the force and pressure around the cylinder at Reynolds numbers ranging from 1 to 8 million. Some observations indicated significant disparities between the simulations conducted using LES (with wall functions) and those performed using the URANS methodology. With the developments in hybrid RANS–LES (see Fröhlich and von Terzi, 2008, for a review), the simulations have become simpler to execute as the resolution in the boundary layer on the cylinder could be significantly reduced over wall-resolved LESs. Detached eddy simulation (DES) and its variants (see Spalart, 2009, for a review) is one such method that showed great promise for cylinder flows. This is shown in the work by Squires et al. (2008), where DES and delayed DES (DDES) for flow over cylinders at transcritical Reynolds numbers are performed and compared with experimental data. Nevertheless, the more recent study by Viré et al. (2020) gives confidence in performing URANS simulations for transcritical Reynolds numbers as envisioned here. Moreover, the studies were also carried out for oscillating cylinders under free and forced vibrations. Their results showed sufficient agreement for the purpose of the study between their URANS and work in the literature, for both experiments and high-fidelity simulations. Most previous works analysed the numerical simulations for flow over cylinders at supercritical and/or transcritical Reynolds numbers for the total force and surface pressure over the cylinder. The vorticity and wake patterns are deduced qualitatively from their various mode shapes and shedding frequencies. But to better understand the fundamental mechanisms of initiation and sustenance of VIV, it is important to look into the details of the various force components and their relations with the cylinder motions.

By quantifying different force mechanisms acting on the cylinder, more insight into the behaviour of oscillations can be gained. The work of Quartapelle and Napolitano (1983) proved to be a fundamental base for the force-partitioning method (FPM). A generic formula was derived for the force and moment of a rigid body immersed in incompressible flow. It was achieved by solving the Navier–Stokes equations where the momentum equation was projected on the gradient of a harmonic function, which satisfies appropriate boundary conditions. Chang (1992) further showed the force contribution due to potential flow, vorticity within the flow, and surface vorticity on the finite body, subsequently providing a further understanding of different force components.

The influence of different forces on the oscillation of a rigid body provides more insight into the nature of the force mechanism. Morse and Williamson (2009) showed that the energy extracted in forced oscillations is closely related to amplitude response in free oscillations. This is found to be useful in complicated nonlinear responses and has been shown for oscillating cylinders by Williamson and Govardhan (2004). This was combined along with the force-partitioning approach to analyse the force mechanisms in the sustenance of VIV by Menon and Mittal (2021). The latter found that the vortex shedding behind the cylinder contributes only during the initiation of VIV, whilst the vorticity in the shear layer leads to the sustenance of VIV. Hence, the vortex-induced force was found to be the most crucial for oscillations during VIV compared to the vortices in the wake behind the cylinder. While the previous study focuses on the flow around a cylinder at Reynolds number ≤ 250, the importance of the force-partitioning method at transcritical Reynolds numbers is yet to be understood, as the flow around wind turbines is characterised by a high Reynolds number ($Re>\mathrm{3.5}\times {\mathrm{10}}^{\mathrm{6}}$).

Even though VIV of circular cylinders has been extensively studied in the past, the study of VIV of cylinders under the influence of high-Reynolds-number flows and a large mass ratio is still limited. Moreover, as investigated by Viré et al. (2020), there is an interchange between self-exciting and self-limiting oscillation when the cylinder is allowed to vibrate freely. It is interesting to see how various force mechanisms contribute to such behaviour during VIV. Hence, the aim of the current work is twofold: (a) to implement the force-partitioning method in the open-source computational fluid dynamics (CFD) software OpenFOAM (Bainbridge et al., 2012) and (b) to determine the nature of the force mechanism that leads to initiation and sustenance of VIV of wind turbine towers in the transcritical Reynolds number regime. Additionally, it is intended to derive a correlation between the oscillations of the wind turbine tower and different force mechanisms using the above-mentioned approaches.

2 Methodology

2.1 Force-partitioning method

The force-partitioning method aims to separate the forces experienced by a body in a fluid into various components in order to better understand the physical mechanisms between structural motions and the surrounding fluid flow. The method used in the current study to analyse the various forces is derived from the previous works of Quartapelle, Mittal, and Zhang (Quartapelle and Napolitano, 1983; Menon and Mittal, 2021; Zhang et al., 2015). This section explains the methodology used for the partitioning of forces in detail.

Figure 2 shows the schematic representation of the domain. The surface of the body is represented as B, where the body is placed in the fluid domain V_f, which is bounded by an outer boundary ∑. The unit normal of every surface point of the object is represented as n. The method starts by introducing a harmonic function ϕ⁽ⁱ⁾, where $i\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}$ represents the direction in which the force is partitioned. The harmonic function corresponds to the potential associated with the flow around the body for a shape and position for time t=τ as shown in Eq. (1):

$$\begin{array}{}\text{(1)}& {\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\varphi }}^{\left(i\right)}=\mathrm{0},\mathrm{with}\mathit{n}\cdot \mathbf{\nabla }{\mathit{\varphi }}^{\left(i\right)}=\left\{\begin{array}{ll}{n}_i,& \mathrm{on}B\left(\mathit{\tau }\right)\\ \mathrm{0},& \mathrm{on}\mathrm{\Sigma }\end{array}\right..\end{array}$$

The derivation starts with the momentum conservation equation written in Lamb–Gromeka form, as shown in Eq. (2):

$$\begin{array}{}\text{(2)}& \mathit{\rho }{\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\mathit{\rho }(\mathit{\omega }\times \mathit{u})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\mathrm{\nabla }(\mathit{u}\cdot \mathit{u})=-\mathrm{\nabla }p-{\displaystyle \frac{\mathrm{1}}{Re}}\mathit{\mu }\mathrm{\nabla }\times \mathit{\omega },\end{array}$$

where t represents the time, and p, ρ, and μ correspond to pressure, density, and dynamic viscosity, respectively. The partitioned forces can be obtained by projecting this momentum equation on the gradient of the harmonic function, ϕ⁽ⁱ⁾, and integrating it over the fluid domain V_f; i.e.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathit{\rho }\underset{V_{\mathrm{f}}}{\int }& \left({\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\mathit{\omega }\times \mathit{u}\right)\cdot \mathrm{\nabla }{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}V+\underset{V_{\mathrm{f}}}{\int }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\mathrm{\nabla }(\mathit{u}\cdot \mathit{u})+\mathrm{\nabla }p\right)\\ & \cdot \mathrm{\nabla }{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}V+\mathit{\mu }\underset{V_{\mathrm{f}}}{\int }(\mathrm{\nabla }\times \mathit{\omega })\cdot \mathrm{\nabla }{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}V=\mathrm{0}.\end{array}\end{array}$$

Upon further mathematical simplifications, Eq. (3) yields the final decomposition of the forces acting on the body F_B. This has been explained in detail in Menon and Mittal (2021). The final force is the summation of the kinematic force (F_K), the vortex-induced force (F_ω), the viscous force (F_σ), the force that results purely from potential flow (F_ϕ), and the force due to the effects of the outer boundary (F_ϵ), as shown in Eq. (4). The coefficient of each force component can be obtained by dividing the force component by the free-stream dynamic pressure, q_∞, and the area, A, namely $C_i=F_i/\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{2}}A\right)$.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& F_{\mathrm{B}}=F_{\mathrm{K}}+F_{\mathit{\omega }}+F_{\mathit{\sigma }}+F_{\mathit{\varphi }}+F_{\mathit{ϵ}}\\ & F_{\mathrm{B}}=\underset{\mathrm{kinematic}\mathrm{force}}{\underbrace{-\mathit{\rho }\underset{\mathrm{B}}{\int }\mathit{n}\cdot {\displaystyle \frac{\mathrm{d}{U}_{\mathrm{B}}}{\mathrm{d}t}}{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}S-\mathit{\rho }\underset{\mathrm{B}}{\int }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_{\mathrm{B}}^{\mathrm{2}}\mathit{n}\cdot \mathrm{\nabla }{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}S}}\\ & +\underset{\mathrm{vortex}\text{-}\mathrm{induced}\mathrm{force}}{\underbrace{\mathit{\rho }\underset{V_{\mathrm{f}}}{\int }\mathrm{\nabla }\cdot (\mathit{\omega }\times \mathit{u}){\mathit{\varphi }}^{\left(i\right)}\mathrm{d}V+\mathit{\rho }\underset{V_{\mathrm{f}}}{\int }{\mathrm{\nabla }}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{u}}_{\mathit{\nu }}^{\mathrm{2}}+{\mathit{u}}_{\mathit{\varphi }}\cdot {\mathit{u}}_{\mathit{\nu }}\right){\mathit{\varphi }}^{\left(i\right)}\mathrm{d}V}}\\ & +\underset{\mathrm{viscous}\mathrm{force}}{\underbrace{\mathit{\mu }\underset{\mathrm{B}}{\int }(\mathit{\omega }\times \mathit{n})\cdot \mathbf{\nabla }\left(\mathit{\varphi }-\widehat{e_i}\right)\mathrm{d}S}}+\underset{\mathrm{potential}\mathrm{flow}\mathrm{force}}{\underbrace{\mathit{\rho }\underset{V_{\mathrm{f}}}{\int }{\mathrm{\nabla }}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{u}}_{\mathit{\varphi }}^{\mathrm{2}}\right){\mathit{\varphi }}^{\left(i\right)}\mathrm{d}V}}\\ & -\underset{\mathrm{force}\mathrm{due}\mathrm{to}\mathrm{flow}\mathrm{and}\mathrm{vortices}\mathrm{at}\mathrm{outer}\mathrm{boundary}}{\underbrace{\mathit{\rho }\underset{\mathrm{\Sigma }}{\int }{\displaystyle \frac{\mathrm{d}\mathit{U}}{\mathrm{d}t}}\cdot \mathit{n}{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}S+\mathit{\mu }\underset{\mathrm{\Sigma }}{\int }(\mathit{\omega }\times \mathit{n})\cdot \mathrm{\nabla }{\mathit{\varphi }}^{\left(i\right)}\mathrm{d}S}}\end{array}\end{array}$$

The velocity of the body at every point is represented as U_B, also shown in Fig. 2. The velocity components u_ϕ and u_ν are obtained using the Helmholtz decomposition to expand the velocity vector into curl-free and divergence-free components as shown in Eq. (5) (Batchelor, 2000):

$$\begin{array}{}\text{(5)}& \mathit{u}={\mathit{u}}_{\mathit{\varphi }}+{\mathit{u}}_{\mathit{\nu }}=\mathrm{\nabla }\mathrm{\Phi }+\mathrm{\nabla }\times \mathit{A},\end{array}$$

where Φ and A are scalar and vector potentials, respectively.

Figure 2 Schematic representation of a fluid domain containing a cylindrical rigid body.

2.2 OpenFOAM implementation

The numerical simulations for the present study are carried out using the open-source computational fluid dynamics (CFD) software OpenFOAM-v2012 (Bainbridge et al., 2012). OpenFOAM consists of many solvers and libraries, making it useful for various applications like fluid dynamics, structural dynamics, electromagnetics, etc. Modifications made in the libraries of OpenFOAM for the implementation of the force-partitioning method (FPM) are explained in the current section.

OpenFOAM consists of an incompressible solver named pimpleFoam, which solves for pressure and velocity according to the PIMPLE algorithm. PIMPLE is a combination of the Pressure Implicit with Splitting of Operator (PISO) and Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithms. The incompressible solver is modified, named fppimpleFoam, to calculate the velocity potential ϕ and the harmonic function ϕ⁽ⁱ⁾. However, to solve ϕ⁽ⁱ⁾, a boundary condition called fppotentialGradient is made, which satisfies Eq. (1). The variable i, being equal to 1, 2, and 3, partitions the force in the x, y, or z direction, respectively. The boundary condition is then given as input for the CFD simulation to consider the direction in which the force is partitioned. The auxiliary function ϕ⁽ⁱ⁾ is introduced in the 0 folder, where the surface can be given the boundary type fppotentialGradient and the direction as (1 0 0) in order to partition the force in the x direction. Finally, a customised function object, fpforceCoeff, is derived from the force class. This function object calculates the different force coefficients for every time step according to Eq. (4), which can later be accessed for post-processing. The equations mentioned above work for an immersed boundary method. If a moving mesh is used to simulate the motion of the object, then a correction term is to be introduced as shown in Eq. (6):

$$\begin{array}{}\text{(6)}& F_{\mathrm{vorticity}\mathrm{correction}}=-\mathit{\rho }{\mathit{\varphi }}^{\left(i\right)}\underset{\mathrm{B}}{\int }\mathit{n}\cdot \left(\mathit{\omega }\times {\mathit{U}}_{\mathrm{B}}\right)\mathrm{d}S.\end{array}$$

This term corrects for the influence of mesh motion in the Lamb vector as seen in Eq. (2), thereby calculating the rotational characteristics purely arising from the fluid flow. The complete flowchart of the implementation is shown in Fig. 3. The source code of the implementation can be found in the repository (VimalKumar, 2023), where the fundamental description of running the OpenFOAM simulation is explained in detail.

Figure 3 Flowchart for the implementation of the force-partitioning method in OpenFOAM.

3 Numerical setup

3.1 Governing equations and numerical methods

The numerical simulations in the current study are performed using fppimpleFoam together with the function object fpforceCoeff. As it is a relatively simple geometry, a mesh with O-grid topology is made to carry out the simulations. The unsteady Reynolds-averaged Navier–Stokes (URANS) equations are solved using the following equations:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{\mathit{U}}& =\mathrm{0},{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{U}}}{\partial t}}+\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \left(\mathrm{\nabla }\stackrel{\mathrm{‾}}{\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}}^{\prime }{\mathit{U}}^{\prime }}\right).\end{array}\end{array}$$

The corresponding velocity and velocity fluctuations are represented by U and ${\mathit{U}}^{\prime }=\mathit{U}-\stackrel{\mathrm{‾}}{\mathit{U}}$, where the overbar represents a suitable Reynolds average. For laminar simulations, ${\mathit{U}}^{\prime }=\mathrm{0}$ and $\mathit{U}=\stackrel{\mathrm{‾}}{\mathit{U}}$, and the equations revert to the traditional Navier–Stokes equations. The time derivatives in the above equation are solved using an implicit second-order backward scheme, and the spatial derivatives are computed using a second-order Gauss linear scheme. For the transcritical Reynolds numbers, the two-equation k–ω shear stress transport (SST) turbulence model (Menter, 1994) is used to compute the turbulent kinetic energy, k; its specific dissipation rate, ω; and from these the eddy viscosity, ν_T. These are then used to close Eq. (7) with

$$\begin{array}{}\text{(8)}& \stackrel{\mathrm{‾}}{{\mathit{U}}^{\prime }{\mathit{U}}^{\prime }}=\mathrm{2}{\mathit{\nu }}_{\mathrm{T}}(\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{U}}-\mathrm{1}/\mathrm{3}\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{\mathit{U}})-\mathrm{2}/\mathrm{3}k\mathbf{I},\end{array}$$

where I represents the identity matrix, which is equivalent to Kronecker delta. A two-dimensional cylindrical section with a spring–mass–damper system, as a way to represent wind turbine towers, is considered to analyse the fluid–structure interaction in this study. The spring and the damper are attached to the centre of mass of the cylinder. The spring has a spring constant of k_s, and the viscous damper has the damping coefficient c. The spring–mass–damper in this study gives the cylinder a single degree of freedom in the crossflow direction. The schematic representation of the system and the mesh are shown in Fig. 4. As OpenFOAM is a code designed for three-dimensional space, the mesh is extruded in the z direction with one cell thickness.

Figure 4 Schematic representation of the computational domain: computational mesh (a), domain characteristics (b).

The governing equations for the motion of this rigid cylinder are represented as follows:

$$\begin{array}{}\text{(9)}& M\text{ÿ}+c\dot{y}+k_{\mathrm{s}}y=F_{\mathrm{B}},\end{array}$$

where ÿ, $\dot{y}$, and y represent the acceleration, velocity, and displacement in the y direction, respectively. The natural frequency of the system is determined by $\mathit{\nu }=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\sqrt{\frac{k_{\mathrm{s}}}{m}}$ for a damping ratio defined as $\mathit{\zeta }=\frac{c}{\mathrm{2}\sqrt{k_{\mathrm{s}}m}}$. The fluid–structure interaction (FSI) simulations are carried out across various reduced velocities U^* for a single mass ratio m^*. The reduced velocity is calculated as $U^{*}=\frac{U_{\mathrm{\infty }}}{\mathit{\nu }D}$ and mass ratio as $m^{*}=\frac{m}{\mathit{\rho }\mathit{\pi }{D}^{\mathrm{2}}L/\mathrm{4}}$.

The structural properties of wind turbine towers must be calculated accurately to perform FSI simulations, which are carried out using modal approximation. With a given variation in mass density and stiffness, the natural frequency of the tower can be calculated. This is calculated from Rayleigh's method using the principle of equating maximum potential energy and maximum kinetic energy in the vibration mode (Clough and Penzien, 1993). The wind turbine is considered to be a stepped tower with n segments to calculate the natural frequency, as shown in Fig. 5. The displacement of the tower at any height z and time t can be written as

$$\begin{array}{}\text{(10)}& v(z,t)=\mathit{\psi }\left(z\right)\cdot y_{\mathrm{0}}\sin \left(\mathit{\omega }t\right),\end{array}$$

where ψ(z) is the shape function, which represents the ratio of the displacement at the point z to the reference displacement y(t). The reference displacement y(t) is represented as a sinusoidal function with amplitude y₀. As mentioned above, the maximum potential energy is equated to the maximum kinetic energy to obtain the natural frequency. The potential energy of the system is given by

$$\begin{array}{}\text{(11)}& V={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}EI\left(z\right){v^{\prime \prime }}^{\mathrm{2}}\mathrm{d}z,\end{array}$$

where EI(z) is the flexural stiffness of the tower and ^′ represents the space derivative. The kinetic energy of the non-uniformly distributed mass is

$$\begin{array}{}\text{(12)}& T={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}m\left(z\right){\dot{v}}^{\mathrm{2}}\mathrm{d}z,\end{array}$$

where $\dot{}$ represents the derivative in time. Upon substituting the displacement from Eq. (10), the maximum potential and kinetic energy can be written as

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& V_{\max }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{y}_{\mathrm{0}}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{L}{\int }}EI\left(z\right)\left({\mathit{\psi }}^{\prime \prime }\right(z){)}^{\mathrm{2}}\mathrm{d}z,\\ & T_{\max }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{y}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{L}{\int }}m\left(z\right)\left(\mathit{\psi }\right(z){)}^{\mathrm{2}}\mathrm{d}z.\end{array}\end{array}$$

Figure 5 Illustration of a stepped tower with element j, length of an element l_j, height Z_j, mass of the element m_j, and moment of inertia I_j.

The natural frequency of the system can be obtained by equating V_max and T_max as

$$\begin{array}{}\text{(14)}& {\mathit{\omega }}^{\mathrm{2}}={\displaystyle \frac{\underset{\mathrm{0}}{\overset{L}{\int }}EI\left(z\right)\left({\mathit{\psi }}^{\prime \prime }\right(z){)}^{\mathrm{2}}\mathrm{d}z}{\underset{\mathrm{0}}{\overset{L}{\int }}m\left(z\right)\left(\mathit{\psi }\right(z){)}^{\mathrm{2}}\mathrm{d}z}}.\end{array}$$

The exact mode shape of a uniform beam of length L undergoing vibration is taken as

$$\begin{array}{}\text{(15)}& \mathit{\psi }\left(z\right)=\left(\mathrm{1}-\cos \left({\displaystyle \frac{\mathit{\pi }z}{\mathrm{2}L}}\right)\right).\end{array}$$

By considering the different mass density m(z) over different turbine height l(z) for n sections and substituting Eq. (15) in Eq. (14), we get

$$\begin{array}{}\text{(16)}& {\mathit{\omega }}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{4}}}{\mathrm{16}{L}^{\mathrm{4}}}}\cdot E\cdot {\displaystyle \frac{\sum _{j=\mathrm{1}}^n{I}_j{l}_j{\cos }^{\mathrm{2}}\left(\frac{\mathit{\pi }{z}_j}{\mathrm{2}L}\right)}{\sum _{j=\mathrm{1}}^n{m}_j{l}_j{\left(\mathrm{1}-\cos \left(\frac{\mathit{\pi }{z}_j}{\mathrm{2}L}\right)\right)}^{\mathrm{2}}}},\end{array}$$

where the moment of inertia of a section of the tower is calculated as

$$\begin{array}{}\text{(17)}& I_j={\displaystyle \frac{\mathit{\pi }\left(D_{\mathrm{out}}^{\mathrm{4}}-{\left(D_{\mathrm{out}}-D_{\mathrm{in}}\right)}^{\mathrm{4}}\right)}{\mathrm{64}}},\end{array}$$

where D_out is the outer diameter of the cylinder, and D_in is the inner diameter of the cylinder. Defining an equivalent moment of inertia and mass per unit length of the stepped beam,

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& I_{\mathrm{eq}}={\displaystyle \frac{\sum _{j=\mathrm{1}}^n{I}_j{l}_j{\cos }^{\mathrm{2}}\left(\frac{\mathit{\pi }{z}_j}{\mathrm{2}L}\right)}{L}},\\ & m_{\mathrm{eq}}={\displaystyle \frac{\sum _{j=\mathrm{1}}^n{m}_j{l}_j{\left(\mathrm{1}-\cos \left(\frac{\mathit{\pi }{z}_j}{\mathrm{2}L}\right)\right)}^{\mathrm{2}}}{L}}.\end{array}\end{array}$$

Considering the tower top mass to be 0, Eq. (16) then condenses to

$$\begin{array}{}\text{(19)}& {\mathit{\omega }}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{4}}}{\mathrm{16}}}{\displaystyle \frac{E{I}_{\mathrm{eq}}}{m_{\mathrm{eq}}{L}^{\mathrm{4}}}}.\end{array}$$

A strong coupling algorithm is utilised for the present study to give accurate results for modelling VIV. The structural parameters are computed for each iteration of the fluid solver, where the fluid solvers are sub-iterated for multiple iterations until the convergence criteria are met. The computed displacement of the cylinder moves the mesh according to a spherical linear interpolation scheme (SLERP).

3.2 Grid convergence study

The structured grid for the computational domain is generated using the blockMesh module of OpenFOAM and is shown in Fig. 4 (left). Richardson's extrapolation method (Roache, 1994) is used for calculating a grid convergence index (GCI) and quantifying the errors from the spatial discretisation of the mesh. This is a standardised method to analyse the numerical solution's dependence on spatial discretisation. Three different grids are compared by having doubled the number of elements in the domain and carrying out h refinement. The first cell height of the mesh, y_wall, is $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m, satisfying the near-wall dimensionless distance (y⁺) to be less than or equal to 1 for every mesh. The domain consists of an O grid of 5 D diameter with the same degree of refinement throughout the meshes with a different number of elements. The ratio of the number of elements in a dimension between two meshes (coarse to medium and medium to fine) is adjusted to 1.26 in order to achieve an approximate total number of points of around two between them. The details of the meshes are given in Table 1.

Table 1 Parameters of three different grids used in the grid convergence study.

The pressure coefficient around the cylinder section is considered for the GCI study for a Reynolds number of 3.6 million, as the paper mainly focuses on the VIV at transcritical Reynolds number regimes. Figure 6b shows the coefficient of pressure along the cylinder surface of the medium grid with the error bars obtained using the GCI. The normalised x coordinate x^* is the x coordinate normalised to the chord of the cylinder, as $x/c$. The maximum error from the GCI study is found to be 2.45 %, but the mean error is calculated to be 0.68 %. The figure also shows the variation in pressure along the cylinder surface for various meshes in comparison to the study from Viré et al. (2020). It can be seen that there is hardly any difference in the pressure along different meshes until the flow separates. However, the coarse mesh underpredicts the pressure coefficient after the separation point in comparison to other meshes. Furthermore, the medium and fine meshes predict the separation point to be 112.2° ($x/c=\mathrm{0.6889}$) and 111.8° ($x/c=\mathrm{0.6857}$), respectively, in comparison to 111° ($x/c=\mathrm{0.6792}$) from the literature. Hence, the meshes are discretised with acceptable error, and further simulations are carried out with the medium mesh, as this mesh balances the complexity of the mesh motion and computational costs.

Figure 6 (a) Pressure coefficient along the surface of the cylinder section for coarse, medium, and fine meshes in comparison to Viré et al. (2020). (b) Pressure coefficient along the surface of the cylinder section with error bars, using the GCI on the medium grid.

4 Results

Simulations are conducted for both laminar and turbulent Reynolds number regimes. Initially, simulations in the laminar regimes are performed to validate and verify the code implemented in OpenFOAM. Subsequently, simulations are carried out to replicate realistic wind conditions encountered by wind turbines. The focus of these simulations is primarily on free vibrations, mirroring real-world scenarios.

Figure 7 Variation in (a) maximum lift coefficient and (b) mean drag coefficient for a fixed cylinder at different subcritical Reynolds numbers compared against the literature: Park et al. (1998), Placzek et al. (2009), and Posdziech and Grundmann (2007).

4.1 Flow in the laminar regime

A cylinder with a diameter of 2 m is considered for the simulations with subcritical Reynolds numbers. The simulations are performed for Reynolds numbers ranging from 80 to 200. The cylinder is considered to be rigid and steady, and hence no structural properties are defined. The time step used for all the simulations maintains a Courant–Friedrichs–Lewy (CFL) number equal to 0.7. The force coefficients are then validated with results from the literature as shown in Fig. 7, where numerical simulations are performed using second-order finite-difference or volume methods (Park et al., 1998; Placzek et al., 2009) and a spectral element method (Posdziech and Grundmann, 2007). The lift and drag coefficients calculated from the sum of the various force coefficients are marked as the force-partitioning method, and those obtained directly from the OpenFOAM simulation are marked accordingly. The lift and drag coefficients calculated from the OpenFOAM post-processing function tool are explained in Eq. (20), where u and v are the velocity components in the x direction and y direction, respectively. The lift and drag forces are defined as F_l and F_d, respectively, and the surface area of the cylinder is defined as S. It can be seen that the sum of the force coefficients matches well with both the data from the literature and the present OpenFOAM simulation results.

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& C_{\mathrm{l}}={\displaystyle \frac{\mathrm{2}{F}_{\mathrm{l}}}{\mathit{\rho }{v}^{\mathrm{2}}S}},\\ & C_{\mathrm{d}}={\displaystyle \frac{\mathrm{2}{F}_{\mathrm{d}}}{\mathit{\rho }{u}^{\mathrm{2}}S}}\end{array}\end{array}$$

The flow around the cylinder subjected to motion in the crossflow direction is also validated and is shown in Fig. 8. For a mass ratio of 2 and a damping ratio of 0.007, the cylinder is subjected to a motion with various reduced velocities between 2 and 15 at a constant Reynolds number of 150. It can be seen that there is a jump in the amplitude at a reduced velocity of 4, where the frequency of vortex shedding matches the natural frequency of the cylinder. The maximum amplitude attenuates for higher velocity ratios and approaches an almost-constant amplitude from a velocity ratio of 10 onwards. This matches well with results in the literature (see e.g. Carmo et al., 2011).

Figure 8 Comparison of the variation in the non-dimensional maximum amplitude of a two-dimensional cylinder for various reduced velocities with results from Carmo et al. (2011); Re=150, $m^{*}=\mathrm{2}$, and c=0.007.

The implementation of the force-partitioning method is verified against the calculation of the total force coefficient calculated by OpenFOAM. The solvers in OpenFOAM are well-validated and benchmarked, which is one of the main motivations for verifying the code implementation with the OpenFOAM solver itself. Figure 9 shows the simulation results along with the partitioned forces. The forces are partitioned in the y direction, corresponding to the direction of the lift force on the cylinder. As shown in the figure, the lift force matches well with the sum of the force coefficients. The slight difference between the total lift coefficient and the sum of the force coefficient is due to the discretisation errors in computing the forces with different methods (Viré and Knaepen, 2009).

Figure 9 (a) Comparison of the lift force coefficient against the sum of the force coefficients from the force-partitioning method. (b) A zoomed-in snippet is shown, which breaks down the total force coefficient into partitioned force coefficients.

When the total force is separated based on different coefficients, it becomes evident that the predominant factor influencing the force is the vortex-induced force. The force due to the vortices formed behind the cylinder is in phase with the total force, and it drives the total contribution. The other major contributors are the forces due to the motion of the cylinder and the viscous force. Depending on the oscillation frequency, the force could be in phase with or out of phase from the motion of the cylinder, which is explained in more detail in Sect. 4.2.1. The results obtained above are in accordance with the results observed by Menon and Mittal (2021) for subcritical Reynolds numbers.

4.2 Flow in the turbulent regime

The laminar flow regime is far from a realistic scenario for wind turbine towers. Therefore, this section analyses the motion due to VIV of a section representative of the tower top of the IEA 15 MW reference wind turbine (Gaertner et al., 2020). The characteristics of the top section of the IEA 15 MW reference wind turbine are chosen as the analysis of a three-dimensional structure in a two-dimensional cylinder is not straightforward, and it simplifies the overall analysis. The tower top section is chosen as it experiences the maximum displacement in the scenario of VIV and is not due to being a critical section during a lock-in.

Following the procedure explained in Sect. 3.1, the natural frequency of the tower is calculated to be 0.48 Hz. The parameters used for the simulation are mentioned in Table 2. The fppimpleFoam solver is employed to simulate and analyse the flow behind the two-dimensional cylinder, which gives insights into the various forces acting on the cylinder for the transcritical Reynolds number regime. The cylinder has 1 degree of freedom in the y direction. For the mass–damper system assumed for the motion of the cylinder, the parameters are tuned to obtain the required structural properties compiled in Table 2. The height parameter shows the single cell thickness in the z direction as mentioned in Sect. 3.1.

Table 2 Simulation parameters used for the free-vibration simulations of the IEA 15 MW reference wind turbine tower section.

In order to carry out the simulations for the turbulent Reynolds number regime, the k–ω SST turbulence model was chosen, as it is reported in the literature to be a good choice for flows with flow separation. The model also calculates the turbulent kinetic energy k and the turbulent dissipation rate ω at every time step. The initial value of the turbulent kinetic energy is calculated as

$$\begin{array}{}\text{(21)}& k={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\left(\mathrm{TI}\cdot U_{\mathrm{\infty }}\right)}^{\mathrm{2}},\end{array}$$

where the turbulence intensity, TI, is taken as 0.03, and the inlet velocity is calculated depending on the reduced velocity used in the simulation. The turbulent dissipation rate is calculated as

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

where $\frac{{\mathit{\mu }}_t}{\mathit{\mu }}$ is taken as 10, and μ is the dynamic viscosity of the fluid. The mesh parameters are explained in Sect. 3.2. The turbulence parameters at the cylinder boundary are taken as $k={\mathrm{10}}^{-\mathrm{10}}$ m² s⁻², and ω is calculated as

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

where β₁=0.075 (Menter, 1992). The boundary conditions for k and ω are Dirichlet boundary conditions at the inlet with values equal to the initial values as calculated above. A Neumann boundary condition is imposed for both the variables at the outlet.

As explained in Sect. 3.2, the simulations performed at Reynolds number $Re=\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{6}}$ serve as a validation for the simulations at the turbulent regime. The force coefficient, pressure distribution, and separation point align with the results observed from the work of Viré et al. (2020).

4.2.1 Free vibration

As mentioned in Sect. 3.1, the simulations at turbulent Reynolds numbers are carried out with a single degree of freedom with translational motion in the y direction. The cylinder has structural properties as mentioned in Table 2, and the simulations are performed for 12 different reduced velocities. The Reynolds number varies from 3 million to 17.2 million.

Figure 10 shows the variation in the force coefficients (Fig. 10a) and the maximum amplitude (Fig. 10b) for various velocity ratios. The Strouhal velocity is represented as U_st, which is defined as the wind velocity for which the vortex shedding frequency equals the natural frequency of the cylinder. The lift coefficient experiences a sharp rise when the velocity ratio $U/U_{\mathrm{st}}$ reaches a value of unity. This also results in a peak in the normalised maximum amplitude of the cylinder, where it almost reaches 0.7. The normalised maximum amplitude becomes significantly less when $U/U_{\mathrm{st}}$ is 1.2 or higher.

Figure 10 Variation in (a) the root mean square of the lift coefficient and mean drag coefficient and (b) the normalised maximum amplitude for various velocity ratios and corresponding Reynolds numbers.

Figure 11a and b show the behaviour of various partitioned forces at a reduced velocity of 2.24. It can be seen from Fig. 11a that the main contributor to the lift force is the force due to vorticity. However, for the laminar (low-Reynolds-number) flow case, other significant contributors to the total force were present, i.e. the force due to added mass and the viscosity of the fluid. The zoomed-in figure from Fig. 11b shows that, for the turbulent (high-Reynolds-number) case without cylinder motion, the influence of the vortex-induced force on the total force is the primary or rather the only significant contributor to the total force.

Figure 11 (a, c, e) Comparison of the lift force coefficient with the partitioning coefficient at U_red=2.24 ($U/U_{\mathrm{st}}=\mathrm{0.71}$), U_red=2.88 ($U/U_{\mathrm{st}}=\mathrm{0.997}$), and U_red=3.20 ($U/U_{\mathrm{st}}=\mathrm{1.05}$), respectively. (b, d, f) Zoomed-in snippets of the time series of (a) ,(c), and (e), respectively.

As the free-stream velocity reaches the Strouhal velocity, the behaviour changes drastically, as seen from Fig. 11c and e. As the lift force develops a beating pattern, the vortex-induced force follows a similar behaviour. But the kinematic force is the dominating force overall. The magnitude of the lift force reduces mainly due to the out-of-phase behaviour between the kinematic force and the vortex-induced force. This is a similar behaviour for both values of the velocity ratio: $U/U_{\mathrm{st}}=\mathrm{0.997}$ and $U/U_{\mathrm{st}}=\mathrm{1.05}$. The kinematic force is equally dominating as the oscillation continues, but any reduction in the total lift force is due to the out-of-phase trend of the two major force contributors, as can be seen in Fig. 11d and f.

Figure 12 shows the vorticity formed behind the cylinder when the lift force develops a beating pattern. Figure 12a and b show the wake pattern at maximum positive displacements in a cycle, where Fig. 12a depicts the time when the vortex-induced force is in phase with the kinematic force. However, Fig. 12b illustrates the wake pattern when the two forces have a significant phase difference. It can be observed that vortices are shed from the bottom part of the cylinder when the forces are in phase. In contrast, vortices are shed from the top of the cylinder when the kinematic force and vortex-induced force are out of phase. The wake pattern development observed here is similar to that reported in previous studies by Viré et al. (2020).

Figure 12 Non-dimensional vorticity magnitude for $U/U_{\mathrm{st}}=\mathrm{1.05}$ at the maximum positive displacement observed when the lift force develops a beating pattern and when the kinematic force and vortex-induced force (a) are in phase and (b) have a phase difference between them.

Figure 11 illustrates the changes in partitioned forces as the velocity ratio approaches and surpasses 1. The force coefficient amplitudes exhibit heightened values as the ratio nears unity, primarily attributed to the vortex-induced force. However, reductions in amplitude are influenced by the kinematic force. This pattern reduces as the ratio exceeds 1.2.

As seen in Fig. 11, the simulations were stopped once the vortex-induced force became completely transient in nature. The system becomes self-exciting, and the fluid imparts energy to the structure if the structural damping coefficient is not enough to suppress the aerodynamic excitement. This can be analysed with a longer simulation time, which is out of the scope of the current study due to increased time and cost. Moreover, the simulations are carried out with 20 cores, where one of the expensive simulations required approximately 4950 CPU hours for a velocity ratio of $U/U_{\mathrm{st}}=\mathrm{0.997}$ and 2520 CPU hours for $U/U_{\mathrm{st}}=\mathrm{0.71}$, using the Intel^® Xeon^® CPU E5-2670 v2 processor.

Figure 13 shows the frequency ratio of forces observed at various normalised velocities. It can be observed that the kinematic force fluctuates at the natural frequency of the structure, even before the structure starts to oscillate at its natural frequency. This provides the reason for the erratic fluctuations in the lift force. This behaviour continues even after the natural frequency is reached until almost a frequency ratio of 1.2. This suggests that it may be important to analyse the phase difference of the forces from the displacement of the structure. Understanding the nature of the phase differences along with the energy mapping can give more insights into the nature of the oscillations experienced by wind turbine tower sections, as is shown in the next section.

Figure 13 Frequency ratio of the forces as a function of normalised velocity $U/U_{\mathrm{st}}$.

4.2.2 Phase difference and energy transfer

The phase difference gives a first look into how energy transfers between the fluid and the structure. From Fig. 14, it can be seen that there is a considerable phase difference between the lift force and the displacement of the cylinder as the velocity ratio reaches unity. By breaking down the contribution from the different forces, it can be seen that the vortex-induced force is the driving factor during this scenario. However, the phase difference between the kinematic force and the cylinder displacement is slightly negative when the velocity ratio is approximately unity. This could be interpreted as energy being added to the system primarily from the force due to the vortex shedding. This shows the influence of vortices in driving the oscillation, especially around the natural frequency of the system. Meanwhile, the viscous force has a constant negative phase difference with a sink around the velocity ratio of unity. This shows the damping nature of the force throughout, and it is significantly influenced when the system oscillates at its natural frequency.

Figure 14 (a) Phase difference between dominant forces and the displacement of the cylinder section for various velocity ratios. (b) Zoomed-in snippet between $\mathrm{0.5}<U/U_{\mathrm{st}}<\mathrm{2}$.

A deeper analysis of the energy transfer between the fluid and the system provides more insights into the nature of each partitioned force. According to the literature (Morse and Williamson, 2009; Menon and Mittal, 2021; Kumar et al., 2016), the non-dimensional energy extracted by the cylinder can be written as

$$\begin{array}{}\text{(24)}& E^{*}=\underset{t}{\overset{t+\mathrm{d}t}{\int }}{C}_i{y}^{*}\mathrm{d}t,\end{array}$$

where C_i is the partitioned force, y^* is the non-dimensional velocity of the cylinder, and dt is the time step. A positive energy transfer means that the system takes energy from the fluid. Furthermore, this provides better insight into the major contributor to the motion of the cylinder, along with the initiation and sustenance of the oscillations. The energy transfer is analysed for the velocity ratios of interest ($U/U_{\mathrm{st}}=\mathrm{0.997}$, 1, 1.05), which are close to the Strouhal velocity. Figures 15–17 show the energy transfer per time step and the mean transfer during the oscillation cycles for the lift force, the kinematic force, the vortex-induced force, and the viscous force.

Figure 15 Non-dimensionalised energy transfer from the lift force and partitioned forces between the fluid and the structure at a reduced velocity of 2.88 ($U/U_{\mathrm{st}}=\mathrm{0.997}$). The continuous dark line shows the moving average of the instantaneous normalised energy, while the shaded curve shows the instantaneous value.

Figure 16 Non-dimensionalised energy transfer from the lift force and partitioned forces between the fluid and the structure at a reduced velocity of 3.04 ($U/U_{\mathrm{st}}=\mathrm{1}$). The continuous dark line shows the moving average of the instantaneous normalised energy, while the shaded curve shows the instantaneous value.

Figure 17 Non-dimensionalised energy transfer from the lift force and partitioned forces between the fluid and the structure at a reduced velocity of 3.20 ($U/U_{\mathrm{st}}=\mathrm{1}/\mathrm{05}$). The continuous dark line shows the moving average of the instantaneous normalised energy, while the shaded curve shows the instantaneous value.

Figure 15 shows the energy transfer when the reduced velocity is 2.88 and the velocity ratio is 0.997. The major energy contributor to the lift force during the onset of VIV is the energy from the vortex-induced force. There is a pattern that is similar to the beating pattern observed in the time series of the lift force. However, the mean non-dimensionalised energy from the vortex-induced force is positive, $E^{*}=\mathrm{1.27}\times {\mathrm{10}}^{-\mathrm{4}}$. As the oscillations become significantly larger, the energy from the added mass contributes to the overall energy in the lift force. Even though there is an overall contribution of energy from the kinematic force, the mean energy transfer as a result of this is almost zero. As the fluid flow considered during the simulations is at a high Reynolds number, the energy from the viscous force is negligible as well.

At a velocity ratio of unity, it can be observed from Fig. 16 that the transient nature of the flow is attained much more quickly than for the other velocity ratios. The non-dimensionalised energy transfer happens from the structure to the fluid until a peak is observed, at t≈320 s, before the system becomes completely transient in nature. Eventually, there is a net mean positive energy transfer from the fluid to the structure. This change in the nature of the energy transfer may be due to the fact that the energy extracted from the shear layer increases as the cylinder oscillates at the natural frequency and that dissipation to the wakes is comparatively less, as explained in Menon and Mittal (2021).

For the velocity ratio $U/U_{\mathrm{st}}=\mathrm{1.05}$, the net energy from the vortex-induced force is largely negative, which implies that it starts to have a damping nature after the vibrations have been initiated. The rest of the oscillations are sustained by the energy from the added mass. As seen from previous figures, the mean energy contribution from the kinematic force is also negligible in this scenario. When the lift force is completely transient, it can be observed that the vortex-induced force increases as seen in previous conditions. Generally, the energy transition between the fluid and the structure shifts relatively swiftly to a fully transient state before reaching the Strouhal velocity.

5 Conclusions

In this work, vortex-induced vibrations of a circular cylinder at a transcritical Reynolds number are analysed using the force-partitioning method. This method dissects the total force into various force components, such as the force due to added mass, vorticity, the viscosity of the fluid, forces associated with the potential flow, and the force due to the boundary of the domain. The methodology is implemented in the CFD software OpenFOAM, where pimpleFoam is modified to fppimpleFoam and the force coefficients are calculated using fpforceCoeff. Unsteady RANS simulations were carried out to analyse the oscillation in the crossflow direction using the k–ω SST turbulence model. The structural parameters corresponded to the tower top section of the IEA 15 MW reference wind turbine tower.

Upon analysing the decomposed forces at Reynolds numbers greater than 3×10⁶, the major forces that contribute to the lift force are the vortex-induced force and the kinematic force. When the cylinder oscillation is sufficiently far from the natural frequency ($U/U_{\mathrm{st}}<\mathrm{0.8}$), the dominant force contributor is the vortex-induced force. But as the frequency ratio increases, the amplitude of the oscillation increases, which leads to increased contribution of the kinematic force. The patterns observed in the total lift force are mainly due to the phase difference between the kinematic force and the vortex-induced force. The oscillations are initiated by the vortex formed behind the cylinder; consequently, the major force contributor in the initial time periods is the vortex-induced force. The kinematic force oscillates at the natural frequency even before the system does, which explains the transient fluctuations in the lift force.

The energy transfer between the system and the fluid gives a better idea of the nature of the force contribution. The mean energy contribution by the vortex-induced force is positive on or before the Strouhal velocity. Therefore, the net energy transfer occurs from the fluid to the system. Even though there is a high energy contribution from the kinematic force, the mean energy transfer between the system and the fluid is negligible. The energy contribution from the viscous force is negligible as the Reynolds numbers are high. This is in contradistinction to the low-Reynolds-number validation case investigated before.

Similar to the studies by Menon and Mittal (2021), even at transcritical Reynolds numbers, the main contributor to the initiation and sustenance of the oscillations is the vortex-induced force. Even though force-partitioning analyses have been carried out in previous studies, investigating the phenomenon of VIV at such high Reynolds numbers with this method is the first of its kind. In the present research, the analysis of energy derived from vortices spans the entire domain, while the impact of vortices in both the near field and the far field has not been considered and was explored in the previous study; see Menon and Mittal (2021). As the shape and size of the cylinder matter predominantly in the oscillation, care must be taken in designing the wind turbine towers. The change in taper ratio and tower top diameter will also significantly influence VIV, making it important to analyse the critical section regarding a potential lock-in. The current work can be extended to analyse how the various force components vary for a three-dimensional tower. It is critical to understand how the vortices in the near field and far field affect the oscillations for a complete tower. Furthermore, the effect of the taper ratio could be investigated with the force-partitioning method, as a tapered tower potentially reduces the occurrence or strength of VIV. As mentioned in the previous section, surface roughness is imparted to mimic flows with high Reynolds numbers at a low-Reynolds-number range. Here, this methodology can shed light on the contribution of surface roughness to VIV at lower Reynolds numbers.