Variable cross-sectional effect on bi-directional blades–tower–soil–structure dynamic interaction on offshore wind turbine subject to wind–wave loads

Abstract

Background

This study introduces a numerical model designed to simulate interactions occurring between a wind turbine's tower and the surrounding soil, as well as between the nacelle, blades, and the surrounding environment. This simulation accounts for both fore–aft and side-to-side movements. To describe these interactions, the model leverages the Euler–Lagrange equations. It calculates wave loads utilizing the Morison equation, with wave data generated based on the JONSWAP spectrum. Furthermore, aerodynamic loads are determined using the blade element moment theory, and the wind spectrum is generated using the Von Karman turbulence model. The tower is represented as a variable cross-sectional beam, employing a two-noded Euler beam element with two degrees of freedom: transverse displacement and rotation, and utilizing Hermite polynomial shape functions.

Results

In a comparative analysis against experimental data, this modified model demonstrates significant enhancements in accurately reproducing the dynamic behavior of wind turbines with variable cross-sectional towers, outperforming models that approximate the tower with a constant cross section. Our findings reveal that the modified model achieves a remarkable improvement of 15% in replicating the tower's dynamic response when compared to the constant cross-sectional models. As a case study, a 5 MW monopile wind turbine with a flexible foundation, specifically the one provided by the National Renewable Energy Laboratory (NREL), is employed to simulate its dynamic response.

Conclusions

This research presents a robust numerical model for simulating wind turbine behavior in various environmental conditions. The incorporation of variable cross-sectional tower representation significantly improves the model's accuracy, making it a valuable tool for assessing wind turbine dynamics. The study's findings highlight the importance of considering tower flexibility in wind turbine simulations to enhance their real-world applicability.

1 Background

The present global political and economic landscape has led to a shift in the scientific community's focus toward renewable energy sources. This shift is primarily driven by escalating costs and the diminishing reserves of fossil fuels. Many nations have devised medium- and long-term strategic plans that prioritize the development and integration of green energy sources in the upcoming decades. Furthermore, recent international climate conferences have underscored the imperative of increased investments in clean energy technology to mitigate greenhouse gas emissions. Wind energy stands out as a remarkably efficient resource within the realm of renewable energy, chiefly due to its cost-effectiveness and wide availability. Among the various wind energy technologies, offshore wind turbines (OWT) have garnered significant attention due to their exceptional energy density. However, despite the evident advantages of OWT when compared to onshore turbines, there remains a crucial need for further research concerning their maintainability, as well as their impact on operational conditions and longevity. This research is vital to optimizing the performance of offshore wind turbines and ensuring the enduring sustainability of this technology. The study of the dynamic behavior of offshore wind turbines (OWT) is a multifaceted and complex subject that encompasses various research areas. These areas include the modeling of the entire wind turbine, the tower (utilizing continuous-multi step or uniform-variable cross-sectional methods), the interaction between the structure and the soil (involving springs, dampers, combinations thereof, or piles), the nacelle–blades system (considering multi-blades, nacelle with gearbox, or simplified mass representations, depending on the study's objectives), and the primary direction of analysis (uni-directional, bi-directional, or torsional). These research domains are essential for gaining a comprehensive understanding of OWT's performance and for optimizing its long-term sustainability. The process of modeling the dynamic behavior of OWT commences with the determination of the fundamental natural frequencies, which provides critical insights into the turbine's behavior. Using these fundamental characteristics as a foundation, simulations are conducted to predict the loads acting on the wind turbine. These loads encompass seismic, hydrodynamic, and aerodynamic forces. The response of various components of the wind turbine is obtained by solving the governing equations used in the modeling process. Based on the acquired responses, improved designs can be proposed through optimization techniques. These optimizations may involve adjustments to blade characteristics, tower dimensions, or the implementation of specific control methods to manage the turbine's response effectively. A substantial body of research is readily accessible concerning the dynamics of wind turbines. This extensive body of work encompasses various facets, such as the dynamics of individual components, their interactions, and the identification of pivotal factors that hold significance. Moreover, a wide array of mathematical and experimental techniques have been proposed to enhance the precision of modeling, all while maintaining computational and experimental efficiency. The overarching goal of these studies has been to furnish a more comprehensive comprehension of wind turbine dynamic behavior and to facilitate the optimization of its performance, ensuring its long-term sustainability.

Thresher [1] conducted a comprehensive review of various methods used to dynamically simulate horizontal-axis turbines in comparison to experimental data. Halfpenny [2] developed an innovative frequency domain model suitable for analyzing both offshore and onshore wind turbines with full flexibility, utilizing the finite element technique. These studies have provided valuable insights into the dynamic behavior of wind turbines and the simulation methods employed, which can contribute to the development of enhanced wind turbine system designs. Ahlstrom [3] delved into the aeroelastic dynamic responses of horizontal-axis turbines by developing a finite element model. Kessentini et al. [4] employed the differential quadrature method (DQM) to create a mathematical model for horizontal-axis turbines with elastic blades and towers, accounting for structural damping and nacelle pitch angles. As this current research specifically focuses on the dynamic interactions between the tower, blades, and nacelle in offshore wind turbines with elastic foundations, we will provide a brief overview of the literature that addresses various aspects related to this topic in the following sections.

Molenaar [5] conducted an extensive review of the theoretical foundations and design possibilities related to wind turbine dynamics, including the dynamics of offshore turbines influenced by waves. Oh et al. [6] conducted a review of various types of wind turbine foundations and the different modeling techniques used to analyze the interaction between the structure and the soil. These studies offer valuable insights and information on various aspects of wind turbine dynamics, including offshore wind turbine foundations and wave effects, which can contribute to the development of improved wind turbine system designs.

Xu et al. [7] utilized finite element analysis to investigate the dynamic and static behaviors of a tower supporting a 600 KW wind turbine, while considering the influence of flexible subgrade soils. Bhattacharya and co-authors [8] summarized findings from a series of 1:100 scale tests of a V120 Vestas turbine supported by two types of foundations: monopiles and tetrapod suction caissons. These outcomes provided valuable insights into long-term performance and issues related to structure–soil interaction. Adhikari and Bhattacharaya [9] presented a closed-form approximate model for estimating the fundamental frequency of wind turbine towers with elastic foundations. The proposed analytical design is based on the theory of Euler–Bernoulli beam columns with flexible end supports. Bhattacharaya and Adhikari [10] introduced innovative experimental procedures for obtaining the parameters required for the dynamic modeling of offshore wind turbines. Their results indicated significant variations in the damping factors and natural frequencies of wind turbine towers depending on the foundation–soil interaction model. Arany et al. [11] developed an analytical model for offshore wind turbines with elastic foundations to obtain a rapid and relatively accurate estimation of the fundamental frequency. This model is suitable for initial design or for confirming results obtained using finite element analysis. Alamo et al. [12] investigated the dynamic impact of foundation constraints on offshore wind turbines, highlighting the influence of soil–structure interaction on the accuracy of the assessed fundamental frequency and the corresponding soil–structure system damping. Kumar and Nasar [13] studied the dynamic behavior of offshore wind turbines subjected to rotor-induced and environmental vibration forces using an Euler–Bernoulli beam-column model with flexible end supports for the tower. Subsequently, Wang et al. [14] conducted an in-depth investigation of wind turbine dynamics by incorporating an enhanced model that considered vibration signal analysis and employed the Euler–Lagrangian approach. This approach allowed for the consideration of dynamic interactions between the nacelle–tower and the tower–foundation system. Bouzid et al. [15] developed a computational methodology that employed nonlinear finite element analysis to simulate the dynamic characteristics of monopile foundations for offshore wind turbines, with a particular focus on incorporating foundation flexibility as a crucial aspect of the model.

Negm and Maalawi [16] utilized an optimization model that conceptualizes the wind turbine tower as a series of fixed, uniform, round tubular segments, incorporating a top mass to simulate the nacelle–blade system. The forces, including aerodynamic, gravitational, and inertial forces, are distributed along the tower within an undeformed inertial coordinate system. Murtagh et al. [17] proposed a model for a wind turbine under prolonged wind loading, representing the turbine as a tower with a top mass encompassing the blade/nacelle system. The applied load is divided into two components: a drag force affecting the tower and a force affecting the top mass, derived from the rotation of the blades treated as cantilever beams subjected to gravity and centrifugal forces. However, this model does not consider soil–structure or blade–tower interactions. Mohammad AlHamaydeh and Saif Hussain [18] employed three-dimensional finite element analysis using SAP2000 software to model a pile–tower–foundation system constructed of reinforced concrete, incorporating a top mass. The model accounts for wind loads under both steady and unsteady conditions, and it is represented as a single-degree-of-freedom (SDOF) system. John Arrigan et al. [19] investigated the variation in blade natural frequency due to changes in blade rotation and stiffness to reduce flapwise vibration through the use of semi-active tuned mass dampers (STMD). The model considers three uniformly cantilevered beams with mass at the root, representing the nacelle and accounting for blade–tower interaction. Arturo Soriano et al. [20] conducted a comprehensive review of nonlinear and linear techniques for dynamic modeling of wind turbines, including actuators, and explored various control methodologies, including nonlinear ones. Liu [21] analyzed the tower–cabin–blade coupling system for wind turbines, establishing the coordinate system and motion equation as a single-degree-of-freedom (SDOF) system. Genov et al. [22] developed a model for a wind turbine experiencing random aerodynamic load, with the tower represented as a viscously damped fixed-bottom beam with a top mass for the nacelle/blade system. The model includes four equations accounting for various dynamic aspects. Fitzgerald and Basu [23] proposed a multi-degree-of-freedom (MDOF) model using the Euler–Lagrangian method, considering interactions between out-of-plane and in-plane blade vibrations and exploring the significance of structure–soil interaction in wind turbine control. Zuo et al. [24] introduced a 3D finite element model for a 5 MW turbine from the National Renewable Energy Laboratory (NREL) subjected to wind, sea wave, and seismic loads, focusing on the tower while excluding the blades from the scope. Staino and Basu [25] proposed a novel control technique merging active tendon and passive pitch controls, using a fixed tower model with bi-directional blade and tower vibrations and considering three types of loads: wind, gravity, and rotor acceleration. Mohammadi et al. [26] employed an adaptive control algorithm to cancel vibration signals after identifying their frequencies in a fixed tower structure, utilizing multiple software codes for different aspects of wind energy system simulation. Hussan et al. [27] summarized methods to reduce jacket vibration in offshore wind turbines exposed to wind, waves, and seismic loads, incorporating multiple tuned mass dampers (MTMDs). Sun [28] explored the control of a wind turbine exposed to wave, seismic, and wind loads using a semi-active tuned mass damper (STMD), with a focus on uni-directional tower vibration and bi-directional blade vibration. Hemmati and Oterkus [29] proposed a model for an OWT considering time-varying STMD and damage conditions, including nonlinear soil–pile interaction. Brodersen et al. [30] applied an active tuned mass damper (ATMD) to a fixed wind turbine tower. Hsu [31] utilized the differential quadrature method to represent the dynamic behavior of rotating wind turbine blades. Park et al. [32] investigated the modal characteristics of rotating blades and examined frequency variations with changes in blade rotation speed. Li et al. [33] developed a nonlinear model for wind turbine blade vibrations, considering longitudinal, in-plane bending, out-of-plane bending, and torsion vibrations, accounting for damping and moment and force distributions on the blade section. Kono et al. [34] studied the upwind influence on wind turbine towers using computational fluid dynamics (CFD) simulations, with a specific focus on tower behavior. El Mouhsine et al. [35] optimized wind turbine blade shapes based on aerodynamic and structural analysis. Júnior et al. [36] evaluated the performance of nonlinear large deformation blade models, including three-dimensional shell finite elements and exact beam models, for both static and dynamic analysis.

Expanding upon the research conducted by El Absawy et al. [37, 38], a multi-degree-of-freedom (MDOF) model has been developed for dynamic analysis of offshore wind turbine towers characterized by varying cross-sectional geometries and elastic foundations. This model takes into consideration the interplay between the tower, nacelle, and blades in both lateral (side-to-side) and longitudinal (fore-and-aft) directions. The proposed model represents the interaction between the structure and the soil by incorporating longitudinal, coupled, and rotational spring–damper elements in both lateral and longitudinal directions. Blade dynamics are modeled by accounting for centrifugal stiffness, bending stiffness, and gravitational effects. Wave loading is assessed using the Morison equation, with wave data generated based on the JONSWAP spectrum. Aerodynamic loads are calculated using the blade element moment theory, while the wind spectrum is generated through the Von Karman turbulence model. The nonlinear governing equations are derived using an Euler–Lagrangian energy-based approach and solved using the Newmark procedure. The impact of the varying cross-sectional geometry of the tower is compared with that of a tower with a constant cross-sectional shape for further analysis.

2 Methods

2.1 Physical model

The research focuses on a 3D wind turbine system comprising a tower, nacelle, and blades. The system has varying cross sections and an elastic foundation. It faces random wave and wind turbulence loads. Wave load is misaligned with the wind direction, represented by β. The tower is an Euler beam, and the nacelle contains rotating blades. The coordinate system is based at the tower's center line and the soil level intersection (Fig. 1).

Fig. 1

Wind turbine model

The three blades coordinate q₁–q₃ denotes the edgewise coordinates in Fig. 2a and q₄–q₆ denotes the flapwise coordinates in Fig. 2b. The rotating blades speed is Ω in which the blades angle ψ_i of the ith blade is given by:

$$\psi_{i} \left( t \right) = \Omega t + \frac{2\pi }{3} \left( {i - 1} \right),\quad i = 1,2,3$$

(1)

Fig. 2

Blade vibration in a edgewise direction and b flapwise direction

Nacelle's relative motion represented by q₇ (fore–aft) and q₈ (side–side). Soil–structure interaction (SSI) is modeled by springs (klx, kly, krx, kry, kxy) and dashpot dampers (clx, cly, crx, cry). OWT is MDOF with 12 DOFs. Blades have six DOFs (flapwise, edgewise). SSI is indicated by four DOFs: q₉–q₁₁ (translation), q₁₀ and q₁₂ (rotation).

2.2 Mathematical model (Euler–Lagrangian equation)

The model is expressed utilizing the Euler–Lagrangian energy-based method:

$$\frac{{\text{d}}}{{{\text{d}}t}} \frac{{\partial T\left( {t,q\left( t \right),\dot{q}\left( t \right)} \right)}}{{\partial \dot{q}_{i} \left( t \right)}} - \frac{{\partial T\left( {t,q\left( t \right),\dot{q}\left( t \right)} \right)}}{{\partial q_{i} \left( t \right)}} + \frac{{\partial U\left( {t,q\left( t \right)} \right)}}{{\partial q_{i} \left( t \right)}} = F_{i} \left( t \right)$$

(2)

where T accounts for the system's total kinetic energy; U is the system's total potential energy; t is the time; $\widetilde{q}$ is the generalized system degree of freedom (coordinates); $\dot{\widetilde{q}}$ is the first derivative w.r.t time of the generalized system degree of freedom (coordinates); and $F$ is the generalized force vector.

2.2.1 Kinetic energy

The kinetic energy of the OWT system is given by:

$$T = T_{{{\text{tow}}}} + T_{{{\text{nac}}}} + T_{{\text{f}}} + T_{{{\text{blades}}}}$$

(3)

where ${ }T_{{{\text{tow}}}}$ represents the kinetic energy of the tower; ${ }T_{{{\text{nac}}}}$ denotes the kinetic energy of the nacelle; ${ }T_{{\text{f}}}$ denotes the kinetic energy of the foundation; and **${ }T_{{{\text{tow}}}}$ is the blades' kinetic energy.

The tower total absolute velocity is:

$$v_{{{\text{tow}}}} = \sqrt {{v_{{{\text{tow}}}}^{{{\text{fa}}}}}^{2} + {v_{{{\text{tow}}}}^{{{\text{ss}}}}}^{2} } = \sqrt {\left( {\dot{q}_{7} \varphi_{{{\text{fa}}}} + \dot{q}_{9} + z \dot{q}_{10} } \right)^{2} + \left( {\dot{q}_{8} \varphi_{{{\text{ss}}}} + \dot{q}_{11} + z \dot{q}_{12} } \right)^{2} }$$

(4)

where z denotes the tower coordinate, $\varphi_{{{\text{fa}}}}$ is the tower normalized fundamental mode shape (first mode shape) at the fore–aft direction, and $\varphi_{{{\text{ss}}}}$ is the tower normalized fundamental mode shape (first mode shape) at the side–side direction.

Through taking an infinitesimal element dz and integrating along the tower height h, therefore $T_{tow}$ is given by:

$$T_{{{\text{tow}}}} = \frac{1}{2} \mathop \smallint \limits_{0}^{h} M_{{\text{tow,d}}} v_{{{\text{tow}}}}^{2} {\text{d}}z$$

(5)

where $M_{{\text{tow,d}}}$ represents the mass per length density of the tower.

The kinetic energy of the nacelle $T_{{{\text{nac}}}}$ is given by:

$$T_{{{\text{nac}}}} = \frac{1}{2} M_{{{\text{nac}}}} v_{{{\text{nac}}}}^{2}$$

(6)

where $M_{{{\text{nac}}}}$ denotes the total mass of nacelle and hub.

The foundation's kinetic energy is $T_{{\text{f}}}$ and given by:

$$T_{{\text{f}}} = \frac{1}{2} M_{{\text{f}}} \left( {\dot{q}_{9}^{2} \left( t \right) + \dot{q}_{11}^{2} \left( t \right)} \right) + \frac{1}{2} I_{{\text{f}}} \left( {\dot{q}_{10}^{2} \left( t \right) + \dot{q}_{12}^{2} \left( t \right)} \right)$$

(7)

where $M_{{\text{f}}}$ and $I_{{\text{f}}}$ are the mass of the foundation and its inertia moment, respectively.

The kinetic energy of the three blades $T_{{{\text{blades}}}}$ is given by:

$$T_{{{\text{blades}}}} = \frac{1}{2} \mathop \sum \limits_{i = 1}^{3} \mathop \smallint \limits_{0}^{R} m_{{{\text{blade}},d}} v_{{{\text{blade}}.i}}^{2} \left( {r,t} \right) {\text{d}}r$$

(8)

where $m_{{{\text{blade}},d}}$ denotes the mass per length density of the blades and R is the blade length. The total absolute velocity of the blade is:

$$v_{{{\text{blade}},i}} \left( {r,t} \right) = \sqrt {\dot{x}_{{{\text{blade}},i}}^{2} + \dot{y}_{{{\text{blade}},i}}^{2} + \dot{z}_{{{\text{blade}},i}}^{2} }$$

(9)

where

$$\begin{aligned} & v_{blade,i} = f\left( {\dot{x}_{blade,i} ,\dot{y}_{blade,i} ,\dot{z}_{blade,i} } \right) \\ & \dot{x}_{blade,i} = v_{nac}^{ss} + \Omega r\cos \left( {\psi_{i} } \right) + \dot{q}_{i} \varphi_{e} \cos \left( {\psi_{i} } \right) - \Omega q_{i} \varphi_{e} \sin \left( {\psi_{i} } \right) \\ & \dot{y}_{blade,i} = v_{nac}^{fa} + \dot{q}_{i + 3} \varphi_{f} \\ & \dot{z}_{blade,i} = - \Omega r\sin \left( {\psi_{i} } \right) - \dot{q}_{i} \varphi_{e} \sin \left( {\psi_{i} } \right) - \Omega q_{i} \varphi_{e} \cos \left( {\psi_{i} } \right) \\ \end{aligned}$$

(10)

where $\varphi_{{\text{e}}} \left( r \right)$ and $\varphi_{{\text{f}}} \left( r \right)$ represent the normalized fundamental mode shape (first mode shape) in both edgewise direction and flapwise direction.

2.2.2 Potential energy

The OWT potential energy is given by:

$$U = U_{{{\text{tow}}}} + U_{{\text{f}}} + U_{{{\text{blades}}}}$$

(11)

where ${ }U_{{{\text{tow}}}}$ represents the potential energy of the tower; ${ }U_{{\text{f}}}$ represents the potential energy of the foundation; and ${ }U_{{{\text{blades}}}}$ is the potential energy of the blades.

The potential energy of the tower $U_{{{\text{tow}}}}$ is determined by:

$$\begin{aligned} & U_{{{\text{tow}}}} = \mathop \smallint \limits_{0}^{h} \frac{1}{2} {\text{EI}}\left( z \right)\left[ {\frac{{\partial^{2} w_{fa} \left( {z,t} \right)}}{{\partial z^{2} }}} \right]^{2} {\text{d}}z + \mathop \smallint \limits_{0}^{h} \frac{1}{2} {\text{EI}}\left( z \right)\left[ {\frac{{\partial^{2} w_{ss} \left( {z,t} \right)}}{{\partial z^{2} }}} \right]^{2} {\text{d}}z \\ & U_{{{\text{tow}}}} = \frac{1}{2} k_{7} q_{7}^{2} + \frac{1}{2} k_{8} q_{8}^{2} \\ \end{aligned}$$

(12)

${ }U_{{\text{f}}}$ is the potential energy of the foundation given by:

$$U_{{\text{f}}} = \frac{1}{2} k_{x} q_{9}^{2} + \frac{1}{2} k_{xr} q_{10}^{2} + \frac{1}{2} k_{y} q_{11}^{2} + \frac{1}{2} k_{yr} q_{12}^{2}$$

(13)

$U_{{{\text{blades}}}}$ is the potential energy of the blades owing to bending, gravity, and centrifugal stiffness in both edgewise and flapwise direction.

$$\begin{aligned} U_{{{\text{blade}},i}} = & \mathop \smallint \limits_{0}^{R} \frac{1}{2} {\text{EI}}_{{\text{e}}} \left( r \right)\left[ {\frac{{\partial^{2} r_{ie} \left( {r,t} \right)}}{{\partial r^{2} }}} \right]^{2} {\text{d}}r + \mathop \smallint \limits_{0}^{R} \frac{1}{2} {\text{EI}}_{{\text{f}}} \left( r \right)\left[ {\frac{{\partial^{2} r_{if} \left( {r,t} \right)}}{{\partial r^{2} }}} \right]^{2} {\text{d}}r \\ & + \; \mathop \smallint \limits_{0}^{R} \frac{1}{2} \left[ { - W\left( r \right) } \right] \left( {\frac{{\partial r_{ie} \left( {r,t} \right)}}{\partial r}} \right)^{2} {\text{d}}r + \mathop \smallint \limits_{0}^{R} \frac{1}{2} \left[ { - W\left( r \right)} \right] \left( {\frac{{\partial r_{if} \left( {r,t} \right)}}{\partial r}} \right)^{2} {\text{ d}}r \\ & + \, \mathop \smallint \limits_{0}^{R} \frac{1}{2} \left[ {F_{r} \left( r \right)} \right] \left( {\frac{{\partial r_{ie} \left( {r,t} \right)}}{\partial r}} \right)^{2} {\text{d}}r + \mathop \smallint \limits_{0}^{R} \frac{1}{2} \left[ {F_{r} \left( r \right) } \right] \left( {\frac{{\partial r_{if} \left( {r,t} \right)}}{\partial r}} \right)^{2} {\text{d}}r \\ \end{aligned}$$

(14)

where $W\left( r \right)$ and $F_{r} \left( r \right)$ are the axial weight component and the centrifugal force along the blade, respectively.

By simplifying the blades potential energy, it can be written as:

$$U_{{{\text{blades}}}} = \frac{1}{2} \mathop \sum \limits_{i = 1}^{3} \left[ {k_{ie} q_{i}^{2} + k_{if} q_{i + 3}^{2} } \right]$$

(15)

where $k_{ie}$ and $k_{if}$ are the stiffness edgewise and flapwise force the ith blade, respectively.

2.3 Loading

The OWT is subjected to several forces as shown in the following.

2.3.1 Wave turbulence loading

The water wave loading on the tower (cylindrical structure) is governed by Morison’s equation [39]. To calculate the force acting on an infinitesimal element dz,

$${\text{d}}F_{{{\text{wave}}}} = \frac{{\pi D^{2} }}{4}C_{M} \rho \dot{u} {\text{d}}z + \frac{1}{2}C_{D} \rho D u \left| u \right| {\text{d}}z$$

(16)

C_D and C_M are the drag and inertia coefficients (C_D = 1, C_M= 1.2); ρ is the water density; D is the mean diameter of dz; and u and u̇ are the wave velocity and acceleration of the fluid particles, respectively. Wave is dominant, uni-directional, and governed by JONSWAP spectrum [40].

$$S\left( f \right) = 0.3125 H_{{\text{S}}}^{2} T_{{\text{p}}} \left( {\frac{f}{{f_{{\text{p}}} }}} \right)^{ - 5} \exp^{{ - \frac{5}{4} \left( {\frac{f}{{f_{{\text{p}}} }}} \right)^{ - 4} }} \left( {1 - 0.287\ln \gamma } \right) \gamma^{{\exp^{{\left( { - \frac{{\left( {\omega - \omega_{{\text{p}}} } \right)^{2} }}{{2 \sigma^{2} \omega_{{\text{p}}}^{2} }}} \right)}} }}$$

(17)

where $H_{{\text{S}}}$ denotes the significant wave height, $f_{{\text{p}}}$ denotes the wave peak frequency $\omega = 2\pi f$, and $T_{{\text{p}}}$ is the wave period $T_{{\text{p}}} = \frac{1}{{f_{{\text{p}}} }}$, in which $\sigma = 0.07$ for $f \le f_{{\text{p}}}$ and $\sigma = 0.09$ for $f > f_{{\text{p}}}$, where $\sigma$ is the input parameter.

$\gamma$ is the peak parameter of the JONSWAP spectrum as follows:

$$\gamma = \left\{ {\begin{array}{*{20}c} 5 & {\frac{{T_{{\text{p}}} }}{{\sqrt {H_{{\text{S}}} } }} \le 3.6} \\ {\exp^{{\left( {5.57 - \frac{{1.15 T_{{\text{p}}} }}{{\sqrt {H_{{\text{S}}} } }}} \right)}} } & {3.6 < \frac{{T_{{\text{p}}} }}{{\sqrt {H_{{\text{S}}} } }} \le 5} \\ 1 & {\frac{{ T_{{\text{p}}} }}{{\sqrt {H_{{\text{S}}} } }} < 5} \\ \end{array} } \right.$$

(18)

JONSWAP spectrum assumes multiple waves with different frequencies and amplitudes. To calculate fluid particle velocity u, sum up contributions for all frequencies.

$$A_{i} = \sqrt {2 S\left( {\omega_{i} } \right)\Delta \omega }$$

(19)

and

$$u\left( {z,t} \right) = \mathop \sum \limits_{i = 1}^{N} \omega_{i} A_{i} \frac{{\cosh \left[ {k\left( {z + d_{{\text{w}}} } \right)} \right]}}{{T_{{\text{w}}} \sinh \left( {kd_{{\text{w}}} } \right)}} \sin \left( { \omega_{i} t - k_{i} x + \emptyset_{i} } \right)$$

(20)

where $\omega$ is the wave frequency; k is the wave number; z is the coordinate as mentioned before from MSL; ${\text{d}}_{w}$ is the water depth; t is the time; x is equal to zero; and $\emptyset_{i}$ is a random phase angle generating and taking values from 0 to 2π, in which N is the number of wave component of spectrum decomposition. k is calculated from the dispersion equation:

$$k \tanh \left( {kz} \right) = \sqrt {2 S\left( {\omega_{i} } \right) \Delta \omega }$$

(21)

Differentiate the velocity equation to get the acceleration $\dot{u}$:

$$\dot{u}\left( {z,t} \right) = \mathop \sum \limits_{i = 1}^{N} \omega_{i}^{2} A_{i} \frac{{\cosh \left[ {k\left( {z + d_{w} } \right)} \right]}}{{T_{w} \sinh \left( {kd_{w} } \right)}} \cos \left( { \omega_{i} t - k_{i} x + \emptyset_{i} } \right)$$

(22)

Using the virtual work principle, the wave force (dF_wave) acting on the virtual displacement (δu_tower) of the immersed tower part is calculated as δW_wave.

$$\begin{aligned} \delta W_{{{\text{wave}}}} = & \smallint dF_{{{\text{wave}}}} \delta u_{{{\text{tower}}}} \\ = & \smallint dF_{{{\text{wave}}}} \left[ {\cos \left( \beta \right)\left( {\varphi_{fa} \delta q_{7} + \delta q_{9} + z\delta q_{10} } \right) + \sin \left( \beta \right) \left( {\varphi_{ss} \delta q_{8} + \delta q_{11} + z\delta q_{12} } \right)} \right] \\ \end{aligned}$$

(23)

where $\beta$ is the angle between wind and wave.

Therefore, the wave force vector is as follows:

$$F_{{{\text{wave}}}} = \left[ {\begin{array}{*{20}c} 0 & 0 & 0 & 0 & 0 & 0 & {F_{7w} } & {F_{8w} } & {F_{9w} } & {F_{10w} } & {F_{11w} } & {F_{12w} } \\ \end{array} } \right]$$

(24)

where

$$\begin{aligned} & F_{7w} = \mathop \sum \limits_{i = 1}^{N} \cos \left( \beta \right)A\left( {z_{i} } \right)\varphi_{fa} \left( {z_{i} } \right) \& F_{8w} = \mathop \sum \limits_{i = 1}^{N} \sin \left( \beta \right)A\left( {z_{i} } \right)\varphi_{fa} \left( {z_{i} } \right) \\ & F_{9w} = \mathop \sum \limits_{i = 1}^{N} \cos \left( \beta \right)A\left( {z_{i} } \right) \& F_{10w} = \mathop \sum \limits_{i = 1}^{N} \cos \left( \beta \right)A\left( {z_{i} } \right)z_{i} \\ & F_{11w} = \mathop \sum \limits_{i = 1}^{N} \sin \left( \beta \right)A\left( {z_{i} } \right) \& F_{12w} = \mathop \sum \limits_{i = 1}^{N} \sin \left( \beta \right)A\left( {z_{i} } \right)z_{i} \\ & A\left( {z_{i} } \right) = \left[ {\frac{{\pi D^{2} \left( {z_{i} } \right)}}{4}C_{M} \rho \dot{u}\left( {z_{i} ,t} \right) \Delta z + \frac{1}{2}C_{D} \rho D\left( {z_{i} } \right) u\left( {z_{i} ,t} \right) \left| {u\left( {z_{i} ,t} \right)} \right| \Delta z} \right] \\ \end{aligned}$$

(25)

where $\Delta z$ is the increment length.

The wave load at the mudline in fore–aft (FA) and side–side (SS) directions is shown in Fig. 3.

Fig. 3

Hydrodynamic load and bending loads at mudline

2.3.2 Aerodynamic turbulence loading

BEM derives aerodynamic force using momentum theory and local blade events. Assume constant force on each element and radial independency. Blades are divided into elements. Reference [41] applies BEM to one rotating blade element and velocities in Fig. 4.

$$\begin{aligned} & \tan \varphi = \frac{{\left( {1 - a} \right)V_{o} }}{{\left( {1 + \dot{a}} \right)\Omega r}} \\ & \alpha = \varphi {-} \theta \\ \end{aligned}$$

(26)

Fig. 4

Annular element of blade element with rotational and axial velocities

In which V_o denotes the total wind speed; r represents the blade element radius; Ω denotes the rotating blades speed; a denotes the axial velocity induction factor; $\dot{a}$ represents the rotational velocity induction factor; V_rel is the relative velocity; θ denotes the blade local pitch angle; α is the angle of attack; and φ is the angle between the relative velocity V_rel and the plane of rotation (flow angle), where:

$$V_{{{\text{rel}}}} = \sqrt {\left( {\left( {1 - a} \right)V_{o} } \right)^{2} + \left( {\left( {1 + \dot{a}} \right)\Omega r} \right)^{2} }$$

(27)

in which the drag force D and lift force L are as follows:

$$L = \frac{1}{2} \rho V_{{{\text{rel}}}}^{2} c C_{l} \& D = \frac{1}{2} \rho V_{{{\text{rel}}}}^{2} c C_{d}$$

(28)

in which $C_{{\text{d}}}$ represents the drag coefficient and $C_{{\text{l}}}$ represents the lift coefficient. The plane of the blade element is concerning with tangential $F_{{\text{T}}}$ and normal $F_{{\text{N}}}$ forces, respectively.

$$\begin{aligned} & F_{{\text{T}}} = L\sin \varphi - D \cos \varphi \\ & F_{{\text{N}}} = L \cos \varphi + D \sin \varphi \\ \end{aligned}$$

(29)

According to the above equations, the axial a, rotational $\dot{a}$ velocity induction factor, and total wind speed V_o are required. Firstly, the induction factors are obtained from the algorithm stated in eight steps in Ref. [50] and through the definitions.

$$\begin{aligned} & a = \frac{1}{{\frac{{4 \sin \varphi^{2} }}{{\sigma C_{n} }} + 1}} \\ & \dot{a} = \frac{1}{{\frac{4 \sin \varphi \cos \varphi }{{\sigma C_{t} }} - 1}} \\ \end{aligned}$$

(30)

where $\sigma$ is the solidity and $C_{{\text{n}}}$ is the normal coefficient and $C_{{\text{t}}}$ is the tangential coefficient, which are expressed as:

$$\begin{aligned} & C_{{\text{n}}} = C_{{\text{l}}} \cos \varphi + C_{{\text{d}}} \sin \varphi \\ & C_{{\text{t}}} = C_{{\text{l}}} \sin \varphi - C_{{\text{d}}} \cos \varphi \\ \end{aligned}$$

(31)

Equations from (26) to (31) are obtained from Ref. [41].

Secondly, the total wind speed V_o is;

$$V_{o} \left( {z,t} \right) = V_{{{\text{mean}}}} \left( z \right) + V_{{{\text{turb}}}} \left( t \right)$$

(32)

where $V_{{{\text{mean}}}} \left( z \right)$ is the mean velocity which is constant and depending consequently on the height z according to the equation of Log Wind Profile:

$$V_{{{\text{mean}}}} \left( z \right) = u_{{{\text{hub}}}} \frac{{\log \frac{z}{{z_{o} }}}}{{\log \frac{{h_{{{\text{hub}}}} }}{{z_{o} }}}}$$

(33)

where $u_{{{\text{hub}}}}$ represents the wind velocity at reference length $h_{{{\text{hub}}}}$ which is the hub length or top tower length and $z_{o}$ is the roughness length and its value is 0.03.

The turbulent wind velocity $V_{{{\text{turb}}}} \left( t \right)$ term is obtained through von Karman spectrum model Ref. [42]. The spectral model is:

$$S_{u} \left( f \right) = \frac{{4 \sigma_{u}^{2} \frac{L}{{u_{{{\text{hub}}}} }}}}{{\left( {1 + 71\left( {\frac{f L}{{u_{{{\text{hub}}}} }}} \right)^{2} } \right)^{\frac{5}{6}} }}$$

(34)

and

$$S_{v,w} \left( f \right) = \frac{{2 \sigma_{v,w}^{2} \frac{L}{{u_{{{\text{hub}}}} }}}}{{\left( {1 + 71\left( {\frac{f L}{{u_{{{\text{hub}}}} }}} \right)^{2} } \right)^{\frac{11}{6}} }} \left( {1 + 189\left( {\frac{f L}{{u_{{{\text{hub}}}} }}} \right)^{2} } \right)$$

(35)

in which L, f, and σ are integral scale parameter; cyclic frequency, and standard deviation, respectively. The integral scale parameter is defined by the turbulence scale parameter $\Lambda_{\upsilon }$ by:

$$L = 3.5 \Lambda_{\upsilon }$$

(36)

where

$$\sigma_{u} = \sigma_{v} = \sigma_{w}$$

(37)

The relative velocity of the wind on a fixed position of one rotor blade, with a radial position of 48 m away from the hub center, is compared to HAWC2 code simulation [43] with a mean wind speed of 11 m/s as demonstrated in Fig. 5.

Fig. 5

The blade 1 relative velocity at a radial position of 48 m away from the hub center compared to HAWC2 code simulation

To determine the aerodynamic force vector, virtual work principle is applied in which virtual work $\delta W_{{{\text{Aero}}}}$ of wind force acted by wind force $dF_{{{\text{Aero}}}}$ on the virtual displacement $\delta u$ as follows:

$$\delta W_{{{\text{Aero}}}} = \mathop \sum \limits_{i = 1}^{3} \left[ {\mathop \smallint \limits_{0}^{R} F_{Ni} \left( {r,t} \right) \left\{ {\varphi_{fa} \delta q_{i + 3} + \delta r_{{{\text{nac}}}}^{{{\text{fa}}}} } \right\}{\text{d}}r + \mathop \smallint \limits_{0}^{R} F_{Ti} \left( {r,t} \right) \left\{ {\varphi_{ss} \delta q_{i} {\text{ cos}}\left( {\psi_{i} } \right) + \delta r_{{{\text{nac}}}}^{{{\text{ss}}}} } \right\}{\text{d}}r } \right]$$

(38)

Therefore, the aero force vector is as follows:

$$F_{{{\text{Aero}},i}} = \frac{{\partial \left( {\delta W_{{{\text{Aero}}}} } \right)}}{{\partial \left( {\delta q_{i} } \right)}}$$

(39)

$$F_{{{\text{Aero}}}} = \left[ {\begin{array}{*{20}c} {F_{1A} } & {F_{2A} } & {F_{3A} } & {F_{4A} } & {F_{5A} } & {F_{6A} } & {F_{7A} } & {F_{8A} } & {F_{9A} } & {F_{10A} } & {F_{11A} } & {F_{12A} } \\ \end{array} } \right]$$

(40)

where

$$\begin{aligned} & F_{iA} = \mathop \smallint \limits_{0}^{R} F_{Ti} \left( {r,t} \right) \varphi_{ss,i} {\text{d}}r \& F_{{\left( {i + 3} \right)A}} = \mathop \smallint \limits_{0}^{R} F_{Ni} \left( {r,t} \right) \varphi_{fa,i} {\text{d}}r,\quad i = 1,2,3 \\ & F_{7A} = F_{9A} = \mathop \sum \limits_{i = 1}^{3} \left[ {\mathop \smallint \limits_{0}^{R} F_{Ni} \left( {r,t} \right) {\text{d}}r} \right]\& F_{8A} = F_{11w} = \left[ {\mathop \sum \limits_{i = 1}^{3} \mathop \smallint \limits_{0}^{R} F_{Ti} \left( {r,t} \right) \sin \left( {\psi_{i} } \right) {\text{d}}r} \right] \\ & F_{10w} = h\mathop \sum \limits_{i = 1}^{3} \left[ {\mathop \smallint \limits_{0}^{R} F_{Ni} \left( {r,t} \right) {\text{d}}r} \right] \& F_{12w} = h\mathop \sum \limits_{i = 1}^{3} \left[ {\mathop \smallint \limits_{0}^{R} F_{Ti} \left( {r,t} \right) \sin \left( {\psi_{i} } \right) {\text{d}}r} \right] \\ \end{aligned}$$

(41)

2.3.3 Solution of the nonlinear equation of motion

Recalling equation (2) and substituting the kinetic and potential energy and differentiate, the following nonlinear equation is obtained:

$$\left[ M \right] \ddot{q} + \left[ C \right] \dot{q} + \left[ K \right] q = F_{{{\text{wave}}}} + F_{{{\text{Aero}}}}$$

(42)

Equation (42) is nonlinear with time-varying terms, involving stiffness matrices [K], mass matrix [M], damping matrix [C], and degrees of freedom vector q and its derivatives $\dot{q}$ and $\ddot{q}$. It also includes force vectors F_wave, F_Seismic, and F_Aero for wave, seismic, and aerodynamic forces, respectively. The Newmark integration scheme [44] is used to solve it, and the results are presented in the next section.

3 Results

3.1 Comparative study

The numerical model is verified with the NREL 5 MW wind turbine [43, 45] phase 2, with soil–structure interaction SSI properties according to Ref. [46], in which the longitudinal stiffness in both x and y directions is $K_{lx} = K_{ly} = 2.574E6{\text{ KN}}/{\text{m}}$; the rotational stiffness in both x and y directions is $K_{rx} = K_{ry} = 2.629E8{\text{ KN}}/{\text{rad}}$; stiffness in coupled direction is $K_{xy} = - 2.253E7{\text{ KNm}}/{\text{m}}$; the longitudinal damping ratio in both x and y directions is $\zeta_{lx} = \zeta_{ly} = 0.6{\text{ \% }}$; the rotational damping ratio in both x and y directions is $\zeta_{rx} = \zeta_{ry} = 0.6{\text{ \% }}$; the blade damping ratio in FA and SS is $\zeta_{bFF} = \zeta_{bSS} = 0.48$; and the tower damping ratio in FA and SS is $\zeta_{TFF} = \zeta_{TSS} = 1{\text{\% }}$.

The comparison of the natural frequencies of the presented nonlinear model with FAST in which a great agreement is obtained is given in Table 1.

Table 1 Natural frequencies of FAST and nonlinear model

Comparison First fundamental frequencies—tower and blade in static conditions (FA and SS). Tower: identical frequencies (0.1% error). Blade: around 1% higher than FAST program, consistent with other models [].

The hydrodynamic response is validated with all simulation models in the technical report [43] and matches with all of them in the following.

Table 2 compares the translation displacement at the mudline in meters and the rotational displacement in radians resulting from the hydrodynamic loads on the monopile. The results from two models, namely NREL FAST CS and CENER FAST CS, are compared to the results obtained from the modified model. The comparison indicates a strong alignment between the results from the modified model and those from the other two models.

Table 2 Monopile displacement of two different models and the presented one

Table 3 displays the results of the blade's steady-state response under aerodynamic loads, and these results are compared to those obtained from the FAST model. The comparison demonstrates a strong agreement between the findings of the two models.

Table 3 The blade flap wise steady-state response

3.2 Hydrodynamic (wave) and aerodynamic (wind) loads effect

For $\beta = 0^\circ$, wind speed of 11.4 m/s and the rotating blades with 12.1 rpm, the wind turbine responses for the acceleration and displacement are shown below, and the response is presented for blade in the FA and SS directions in Fig. 6.

Fig. 6

The blade acceleration and displacement response

Figure 6 focuses on the blade's acceleration and displacement response, in which the acceleration in the FA direction is greater than the acceleration in the SS direction.

Figure 7 illustrates the acceleration and displacement responses of the tower. Notably, in the fore–aft (FA) direction, the response dominates (details to be discussed later). In this FA direction, the acceleration oscillates around zero, while the displacement predominantly fluctuates in the positive direction. This indicates that the displacement continuously oscillates about a positive value due to the continuous airflow in this direction. Furthermore, in the FA direction, the acceleration is significantly larger than in the side–side (SS) direction. Additionally, the displacement and acceleration exhibit values of approximately the same magnitude in this direction.

Fig. 7

The Nacelle acceleration and displacement response

In Table 4, the mean displacement for blade and nacelle in the FF and SS directions for different $\beta$ is given, where $q_{1}$ is the displacement of the first blade in the SS direction; $q_{4}$ displacement of the first blade in the FA direction; $q_{7}$ displacement of the nacelle in the FA direction; $q_{8}$ displacement of the nacelle in the SS direction; $r_{{{\text{nac}}}}^{{{\text{fa}}}}$ the tower absolute displacement in FA direction; and$r_{{{\text{nac}}}}^{{{\text{ss}}}}$ the absolute displacement of the tower in side–side direction.

Table 4 The mean displacement for blade and nacelle in the FF and SS directions for different $\beta$

The various values of β represent the direction of wave forces. It is worth noting that the impact of the wave loads primarily influences the side–side (SS) direction. This explains why the values for the SS direction in the tower are generally of smaller magnitude compared to the fore–aft (FA) direction.

The calculated responses indicate that:

The maximum displacement of the tower in the FA direction is reduced by 42% due to the misalignment of forces.

The maximum displacement of the tower in the SS direction is significantly affected by changes in direction, with an increase of approximately 140% between no misalignment and a complete change in direction by 90° for wave loads.

For the blades, the maximum displacement in both the SS and FA directions shows slight changes between no misalignment and a complete change in direction by 90° for both wave loads.

The mean displacement remains unaffected by misalignment, while the maximum values are significantly altered due to misalignment for the tower only.

3.3 Effect of variable cross section

The following table shows the effect of variable cross of wind turbine tower for $\beta = 90^\circ$.

Table 5 shows the significant effect: variable cross-sectional tower vs. constant cross section. Nacelle FA displacement: 15–16% error at β = 0°, ~ 18% increased relative displacement. Nacelle SS displacement: 8–9% error. Blades: < 1% error in SS, almost zero error in FA. Constant cross-sectional overestimates displacement.

Table 5 Worthiness of variable cross-sectional consideration

4 Discussion

In this study, we have undertaken the development of a sophisticated multi-degree-of-freedom (MDOF) model tailored for offshore wind turbine towers. This model is distinct in its consideration of towers with non-uniform cross sections and flexible foundations. The core components of our model encompass the intricate interactions between the tower, nacelle, and blades, as well as the dynamic coupling between the structure and the seabed, a relationship we have represented using a combination of springs and dampers. We have also paid close attention to the complex dynamics of the turbine blades, taking into account centrifugal forces, bending stiffness, and the influence of gravity on their behavior. To simulate the impact of ocean waves, we have employed the renowned Morison equation, complemented by the Joint North Sea Wave Project (JONSWAP) spectrum for side–side (SS) direction loading.

Our results, as presented before, shed light on the mean displacement values for both the blade and nacelle in the fore–aft (FA) and SS directions across a range of β angles.

Upon meticulous examination of our findings, it becomes evident that the maximum displacement of the tower in the SS direction is profoundly influenced by variations in β. We observe a striking approximate 140% increase in displacement when transitioning from a scenario of no misalignment to a scenario involving a complete 90° change in direction with regard to wave loads.

In contrast, the blade's maximum displacement values in both the SS and FA directions exhibit only marginal variations when comparing the no misalignment scenario to a complete 90° change in direction for wave loads. The overall conclusion here is that blade displacement remains relatively stable in the face of misalignment.

It is essential to underline that while the mean displacement of the system remains comparatively consistent across various misalignment scenarios, the most significant shifts are concentrated within the tower's response. These fluctuations in tower displacement serve as a vital point of consideration for structural engineers and designers.

In summary, our research has offered a deep insight into the intricate dynamics of offshore wind turbine towers. It elucidates the profound impact that different orientations and misalignments can have on structural behavior. These findings hold considerable significance for the offshore wind energy sector, impacting tower design and structural integrity. This study underscores the need for further research aimed at devising strategies to mitigate the impact of misalignment on tower displacement, with a broader objective of enhancing the overall performance and safety of offshore wind turbines.

5 Conclusion

A modified multi-degree-of-freedom (MDOF) model has been developed to analyze offshore wind turbine towers with varying cross sections and elastic foundations. This model accounts for the interaction between the tower, nacelle, and blades in both the fore–aft (FA) and side–side (SS) directions. The interaction between the structure and the soil is represented using longitudinal and rotational springs and dampers in both FA and SS directions. Blade dynamics are incorporated in the model by considering centrifugal stiffness, bending stiffness, and gravity. Wave loads are calculated using the Morison equation, with wave data generated from the JONSWAP spectrum. Aerodynamic loads are determined using the blade element moment theory, and the wind spectrum is generated using the Von Karman spectrum. The nonlinear governing equations are derived using an Euler–Lagrangian energy-based method and solved using the Newmark method. The primary objective of this study is to quantitatively assess the impact of varying cross-sectional shapes of wind turbine towers on the simulated response values of the real NREL 5 MW offshore wind turbine. The results of this study indicate that simulating the tower with varying cross sections, as opposed to assuming a constant cross section, can have significant implications. The error increases by 20% when a constant cross-sectional assumption is applied. The study systematically examines the response of the offshore wind turbine to various load types and different misalignments of wave and seismic loads with aerodynamic loads, presenting a detailed analysis of the results.

6 Recommendations

• The tower can be modeled as viscoelastic material instead of flexible foundations.

• New numerical techniques can be used as mentioned in Ref. [47,48,49,50,51,52] in its stability as in Ref. [53].