## Abstract

An analytical solution of the dynamic response of offshore wind turbines under wave load with nonlinear Stokes’s wave theory and wave–structure and soil–foundation interactions is developed. Natural frequencies and the corresponding modes are obtained. The effect of the wave–structure interaction, the added mass, the foundation stiffness, and the nacelle translational and rotational inertia on the motion of the structure is investigated. The nonlinear loading provided by the drag term of Morison’s equation is successfully handled. A parametric study to examine the effect of the structural parameters on the dynamic response is conducted, and the results of the proposed analytical solution are compared to numerical ones. The proposed method has the following advantages: (a) it is accurate and straightforward because of its analytical nature, (b) it does not ignore the drag term in the wave loading by keeping its nonlinearity nature, (c) the structure of the wind turbine is modeled as a continuous system, (d) it takes into account the effect of the rotational and translational inertia of the nacelle on the dynamic response, and (e) it provides an interpretation of the effect of the sea level variation in changing the natural frequencies.

## 1 Introduction

Passing from carbon emission to carbon neutral energy sources motivates many countries to invest in renewable energies. Among them, offshore wind turbines (OWTs) have been chosen as the main source of renewable energy for many European countries [1]. Being installed in the offshore and nearshore mostly by fixed support structures poses challenges in optimum designing. Furthermore, the prediction of the long-term structural behavior of these structures required for structural integrity assessment is not available because the majority of the OWTs have been installed in the past 10 years, and the data acquisition is very costly.

OWTs are subjected to coupled dynamic phenomena due to the interaction of wind and wave loads and the rotor’s vibration. Therefore, dynamic stability and vibration control [2,3] is a topic that attracts the attention of researchers. Pioneering works on vibrational control of large real-life structures under dynamic loading have been performed by Adeli and Saleh [4–7]. Reviews of advances in vibration control algorithms for smart structures up to 2017 are presented by Refs. [8–10]. Since the coupling phenomena are complicated to be modeled by regular deterministic techniques, machine learning (ML) algorithms seem to be a promising alternative tool. Adeli and associates developed intelligent control algorithms employing neural networks and ML techniques as far back as 2008 [11,12]. More recent work on intelligent control of large real-life buildings and bridge structures was presented by Refs. [9,13–19] to discuss and advance the novel concept of integration of vibration control, health monitoring of structures [20,21], and energy harvesting [22] for smart cities of the future.

A key technology in the field of structural engineering in recent years is automated structural health monitoring (SHM). There are two fundamentally different approaches to SHM, one based on vibration [23] and the other based on imaging and computer vision [24]. Vibration-based SHM technology requires adroit integration of vibration theory, signal processing such as wavelets [25], and machine learning [26,27]. SHM technology has been used successfully, and there is a significant body of research on health monitoring of building structures and bridges [28], dams [29], railways [30], pavements [31], retaining structures along highways [32], and tunnels [33], but little work has been reported on health monitoring of offshore structures and wind turbines. A method of vibration-based SHM is based on computing and monitoring the structural properties such as natural frequencies and mode shapes [34]. Besides, the remaining fatigue life is also necessary to be evaluated and monitored in the case of offshore structures. One of the main challenges in fatigue life estimation, apart from selecting a novel accumulating damage model [35–37], is the availability of the stress history resulting from dynamic response of the structure in the hot spots.

The stress history required for fatigue life evaluation can be obtained from two main sources. One source would be simply by measuring the stress history of the real operating structure. Although the field data are extremely valuable because it reflects the real behavior of the structure, harvesting field data are expensive and in some cases, impossible. An alternative source would be generating data through the mathematical models by simulating the real operating situation. An immediate method in providing the response or stress history of an OWT structure, including all its complexities, is to utilize the numerical models by discretizing the structure via finite elements methods. Many researchers [38–48] are using commercial software packages or numerical platforms to simulate these structures due to their accessibility and ability to create high-fidelity models. However, numerical solutions are still time consuming, and their performance is affected by stability issues especially in very complicated time domain problems. The other alternative to generating response data is to develop the analytical solutions to achieve the response as a single function by which the response is evaluated for arbitrary loading parameters. Analytical solutions are straightforward, provide always-true and reliable results as opposed to the numerical ones but challenging to achieve because of complex mathematics.

Unlike a large number of existing numerical models, very few analytical models of the OWTs have been published. Scientists such as Graff [49] and Meirovitch [50] started early works on the classical methods of solving the wave equation in a cantilever flexural beam in the 60s and 70s. They introduced some methods for solving these equations. For instance, Graff [49] listed five methods: (1) finite Fourier transform, (2) expansion in the natural modes (in the spatial domain), (3) Laplace transform, (4) Laplace transform-natural mode expansion, and (5) solution by the natural modes (in both spatial and temporal domain). In the first two methods, it is assumed that the solution contains two separate parts in time and space. While for the last three ones, the solution begins directly from the equation of motion with an arbitrary function. Among these methods, expansion in the natural modes relies on the structure’s natural modes, which can be found from the homogeneous form of the equation of motion with respect to its boundary conditions. This method also provides the result in the form of a single function as opposed to numerical solutions where a set of numbers is obtained as the response. Expansion in the normal modes is also flexible in the external load situations. Complex loading can be expanded by using the Fourier series and analyzed. Selecting other methods requires more effort and sometimes impossible to find the solution. For instance, the method of Laplace transforms requires inverse Laplace transform, which is nearly impossible to be determined analytically for complicated boundary conditions. Pavlou [51] developed an analytical solution for the evaluation of the response of the OWT under the linear waves. In this work, the translational and rotational inertial effect of the nacelle, the hydrodynamic damping, and the soil–foundation interaction have been analytically investigated for gravity-based supported structures. The achieved analytic inversion of Laplace transforms was very challenging in this analysis. Apart from mentioned methods, very few analytical solutions have been proposed in the past years. In a rare case, Wang et al. [52] developed a mathematical model for dynamic analysis of an onshore wind turbine by using the thin-walled theory [45,53–55] to simulate the comprehensive behavior of the wind turbine.

Expansion of the response in the natural modes requires having an accurate and reliable estimation of the natural frequencies of the OWT structure. In the past few years, researchers have attempted to provide models in which the realistic situation of an OWT is included. Most of the works have been focused on the soil and foundation situation. In a study conducted by Arany et al. [56], an analytical model was developed by simulating foundation flexibility using three springs. The effect of boundary conditions on the natural frequencies has been parametrically studied by defining some nondimensional parameters. Having compared the analytical natural frequencies with the measured ones from the actual OWTs, a slight inaccuracy was reported, which is not improved by modeling the tower with the Timoshenko beam theory. In another study, Arany et al. [57] proposed a simplified methodology to have a quick hand calculation of the first natural frequency of an OWT. In their works, the effect of the fluid–structure interactions of the added mass has not been considered, which may be the reason for the inaccuracy they have reported. Amar Bouzid et al. [58] established a nonlinear finite element model [59] to obtain the head stiffness of the monopile support structure at the mudline. Their purpose was to improve the accuracy of the results obtained by Arany et al. [56,57]. Recently, Alkhoury et al. [60] established a full 3D model for the DTU 10 MW OWT including all the details of the nacelle, blades, rotor, as well as full 3D modeling of the soil inside and outside of the monopile and cross-sectional variation of the tower structure. More details about this work will be presented later in this article for the sake of verification.

The natural frequencies of the OWT were also measured in the real operating systems. Damgaard et al. [61] have reported the cross-wind modal properties of an OWT. They have reported that the first natural frequency is time dependent, which might be because of erosion of the soil around the monopile or soil scouring. Later, variation of the natural frequency in time attracts the attention of Prendergat et al. [62,63], resulting in two publications. In their first work [62], the scouring effect on the natural frequency was investigated without considering the effect of added mass in the system. In their second work [63], however, they considered other factors such as water added mass influencing the dynamic properties of the system. Another investigation on the measured data conducted by Dong et al. [64] also reported the time-dependent dynamic properties of an OWT. Moreover, a precise analysis was performed on the measured data by Norén-Cosgriff et al. [65]. The first natural frequency was separately plotted versus the wave height and wind speed based on 20 min measured data in the calm sea condition, while the wind velocity was under the cut in speed. The trend reveals that the first natural frequency is reduced as the wave height and wind speed increase. This phenomenon raises suspicion about the effect of added mass on the natural frequencies by sea level variation.

Attempts toward providing an accurate response require accurate and realistic inputs in the analysis alongside accurate and realistic models and solutions. In the case of bottom-fixed OWTs, which are mostly installed in shallow to intermediate water depth, measurements and studies reveal that nonlinear wave theory should be implemented to simulate the realistic sea states. A study conducted by Natarajan [66] showed that using the second-order wave theory significantly increases the extreme loading on the monopile support structure of OWTs. Wang [67] utilized a transformed linear method to simulate the second-order irregular wave to obtain the wave load by including the sea bottom effect. Wang et al. [68] conducted a case study to investigate the ultimate wave loads on a 10 MW OWT. All of them have reported a significant increase in wave load on the structure.

In the present work, a new analytical solution for the modal analysis of OWT structures is presented. Nonlinear waves and wave–structure and soil–foundation interactions are accounted for in the solution. The consideration of the second-nonlinear wave kinematics improves the ability of the method to cover a higher range of wavelength and height, providing more realistic loading on the structure. OWTs are subjected to different types of environmental loading such as wind, waves, ocean currents, earthquakes, and ship collisions [69]. Among them, dealing with the wave loads is still a challenging task due to its complexities and uncertainties [70]. Therefore, this study focuses on the wave load. Moreover, the analytical modeling of the translational and rotational inertia effect of the nacelle and the fluid–structure and soil–foundation interaction improves the reliability of the dynamic simulation. The novelty of the presented work beyond the published literature is based on the following advantages: The solution is accurate and straightforward because of its analytical nature, it does not ignore the drag term in the wave loading, the structure of the wind turbine is modeled as a continuous system by including its geometrical discontinuities, it takes into account the effect of the rotational and translational inertia of the nacelle, and it provides an interpretation of the effect of the sea level variation on the natural frequencies. The research presented in this article can be extended for health monitoring of offshore structures and wind turbines, which is intended as part of future research by the authors.

## 2 Formulation of the Problem

A typical horizontal-axis OWT consists of nacelle and blades systems mounted on the top of a tower fixed to the seabed by a monopile, as illustrated in Fig. 1(a). A transition part connects the tower and monopile at sea level. In this article, an OWT is modeled as a cantilever column supported by a set of springs at one end and free at the other end (see Fig. 1(b)). This cantilever beam, which is called the system for the rest of this article, consists of two parts separated from the platform level at the top of the transition part. This is because their dimensions and properties can be significantly different. The cross-sectional properties of each part are assumed constant. The definition of symbols used in this article is presented in Table 1. The monopile under the seabed is also modeled by a set of four springs representing the lateral, rotational, cross-coupled, and vertical stiffnesses with constants of *K*_{L}, *K*_{R}, *K*_{LR}, *and**K*_{z}, respectively [71].

Symbol | Structural properties |
---|---|

L | Nacelle level from the seabed (m) |

L_{Tow} | Tower length (m) |

D_{Tow} | Tower average diameter (m) |

t_{Tow} | Tower average thickness (m) |

m_{Tow} | Tower mass of unit length (kg/m) |

E_{Tow} | Tower Young's modulus (GPa) |

EI_{Tow} | Flexural rigidity of the tower, i.e., E_{Tow}I_{Tow} (GPa · m^{4}) |

M_{n} | Nacelle-Rotor assembly mass (kg) |

J_{p} | Nacelle-Rotor assembly rotational inertia (kg · m^{2}) |

L_{Plat} | Platform level from the seabed (m) |

D_{Mon} | Monopile average diameter (m) |

t_{Mon} | Monopile average thickness (m) |

A_{Mon} | Monopile cross-sectional area (m^{2}) |

m_{Mon} | Monopile mass of unit length (kg/m) |

E_{Mon} | Monopile Young's modulus (GPa) |

EI_{Mon} | Flexural rigidity of the monopile, i.e., E_{Mon}I_{Mon} (GPa · m^{4}) |

ρ_{s} | Material density (kg/m^{3}) |

Symbol | Support stiffness |

K_{L} | Lateral stiffness (GN/m) |

K_{LR} | Cross stiffness (GN) |

K_{R} | Rotational stiffness (GN · m) |

Symbol | Hydrodynamic loading properties |

d | Water depth (m) |

C_{D} | Drag coefficient |

C_{A} | Added mass coefficient |

C_{M} | Inertia coefficient |

ρ_{w} | Sea water density (kg/m^{3}) |

λ | Ocean wavelength (m) |

H | Ocean wave height (m) |

Symbol | Structural properties |
---|---|

L | Nacelle level from the seabed (m) |

L_{Tow} | Tower length (m) |

D_{Tow} | Tower average diameter (m) |

t_{Tow} | Tower average thickness (m) |

m_{Tow} | Tower mass of unit length (kg/m) |

E_{Tow} | Tower Young's modulus (GPa) |

EI_{Tow} | Flexural rigidity of the tower, i.e., E_{Tow}I_{Tow} (GPa · m^{4}) |

M_{n} | Nacelle-Rotor assembly mass (kg) |

J_{p} | Nacelle-Rotor assembly rotational inertia (kg · m^{2}) |

L_{Plat} | Platform level from the seabed (m) |

D_{Mon} | Monopile average diameter (m) |

t_{Mon} | Monopile average thickness (m) |

A_{Mon} | Monopile cross-sectional area (m^{2}) |

m_{Mon} | Monopile mass of unit length (kg/m) |

E_{Mon} | Monopile Young's modulus (GPa) |

EI_{Mon} | Flexural rigidity of the monopile, i.e., E_{Mon}I_{Mon} (GPa · m^{4}) |

ρ_{s} | Material density (kg/m^{3}) |

Symbol | Support stiffness |

K_{L} | Lateral stiffness (GN/m) |

K_{LR} | Cross stiffness (GN) |

K_{R} | Rotational stiffness (GN · m) |

Symbol | Hydrodynamic loading properties |

d | Water depth (m) |

C_{D} | Drag coefficient |

C_{A} | Added mass coefficient |

C_{M} | Inertia coefficient |

ρ_{w} | Sea water density (kg/m^{3}) |

λ | Ocean wavelength (m) |

H | Ocean wave height (m) |

In this study, the motion of the system is considered as a lateral deflection due to the wave load applied up to the mean sea level (MSL). As shown in Fig. 1(c), the deflection, *x*(*z*, *t*), is defined as a continuous function of time, *t*, and space, *z*, implying to represent the motion of the system with an infinite degree-of-freedom.

### 2.1 The Governing Equation.

*q*(

*z*,

*t*) is the external forces acting perpendicular to the monopile’s longitudinal axis.

*EI*(

*z*) and

*M*(

*z*) are the flexural rigidity and mass per unit length of the system, respectively. For the system defined in Fig. 1, they are defined as follows:

*H*

_{i}(

*z*),

*i*= 1, 2, 3 are the step functions defined as follows:

*Du*(

*z*,

*t*)/

*Dt*= ∂

*u*/∂

*t*+

*u*∂

*u*/∂

*x*+

*w*∂

*u*/∂

*z*, in the inertial term can be reduced to $u\u02d9(z,t)=\u2202u(z,t)/\u2202t$ by neglecting the advocative terms, which are reported to slightly increase the load when they are included [70]. The inertia term reveals that the $\u2212\rho wCAAMonx\xa8(z,t)$ provides an additional mass to the system affecting its oscillating properties. This added mass can represent itself in the equation of motion, Eq. (1). By substituting Eq. (7) into Eq. (1), the result will be

Equation (8) represents a partial differential equation governing the motion of the system subjected to the wave load. The added mass is included on the right side of the equation. Thus, the left side is not dependent on the motion of the tower. Parameter *a*_{MU} is the mass of the system underwater including added mass. Therefore, the wave–structure interaction is included in the equation of motion.

In the aforementioned three equations, *x*_{Tow}(*z*, *t*), *x*_{MA}(*z*, *t*), and *x*_{MU}(*z*, *t*) stand for the lateral motion of the system in the tower, monopile above water and monopile underwater, respectively. Equations (14) and (15) are homogeneous partial differential equations implying that the motion above water is a kind of free vibration, whereas the underwater motion, represented in Eq. (16), is a forced vibration due to the external load. Therefore, the above water motion is activated by the motion of the underwater part via a series of boundary conditions, which will be introduced later.

### 2.2. The Boundary and Initial Conditions.

*K*

_{z}, can be neglected because the vertical motion of the system is negligible. The three remaining springs can be collected in the matrix form to obtain the shear force,

*F*, and bending moment,

*M*, at the seabed by the following equation [71]:

The rotational motion of the nacelle produces the momentum proportional to the tower’s rotational acceleration $x\xa8Tow\u2032(L,t)$ at the top of the system, which should be in equilibrium with the total moment of the system producing boundary condition in the form of Eq. (20). Besides, the translational acceleration of the nacelle, $x\xa8Tow(L,t)$, creates an inertial force that should be equal to the internal shear force of the tower at the nacelle level. This equilibrium is represented by Eq. (21).

The last eight boundary conditions are based on the system’s continuity on deflection and slope of the motion at sea level as well as the internal shear force and bending moment continuity of the system at a point where two systems are linked together.

## 3 Solution for the Equation of Motion

The method that has been chosen to solve Eq. (8) is to expand the response in the natural modes of the system. The natural frequencies of the system and consequently the natural mode shapes will be obtained. Then, they will be utilized in the solution for the forced vibration.

### 3.1 Natural Modes.

*n*in quantity indicates that it belongs to the

*n*th natural mode. So,

*x*

_{Tn}(

*z*,

*t*),

*x*

_{MAn}(

*z*,

*t*), and

*x*

_{MUn}(

*z*,

*t*) are the natural mode shapes of the tower, monopile above water, and monopile underwater, respectively. The solution for the aforementioned equations is proposed in the form of following equations:

*β*

_{Tn},

*β*

_{MAn}, and

*β*

_{MUn}are the wavenumbers for the tower and monopile above water and underwater, respectively. Equations (40)–(42) can be separated to form the following equations:

*T*

_{Tn}=

*T*

_{MAn}(

*t*) =

*T*

_{MUn}(

*t*). By comparing Eqs. (46) and (47), one can conclude that $aTow2\beta Tn4=aMU2\beta MUn4=aMA2\beta MAn4$ because

*T*

_{Tn}=

*T*

_{MAn}(

*t*) =

*T*

_{MUn}(

*t*). Therefore,

Since *T* is a periodic function, it will oscillate with the cyclic frequency of *ω*_{n}, which is a natural frequency of the system.

*α*

_{i},

*i*= 1…6 are the solution variables introduced to implement the soil–structure interactions, nacelle-blades mechanical properties, and the systems mechanical properties at the sea level.

*T*

_{i},

*A*

_{i}, and

*U*

_{i}are the constant coefficients for each natural mode shape. Substituting the proposed solutions expressed by Eqs. (57)–(59) into the 12 boundary conditions represented by Eqs. (51)–(56) will yield a system of 12 linear equations. In the matrix form, they can be represented as follows:

**is the matrix containing trigonometrical and hyper trigonometrical functions and**

*P***is the constant coefficients vector. The concept of natural modes oscillation is that the oscillation should be independent of the constant coefficients in Eqs. (57)–(59). Therefore, the determinant of matrix**

*D***should be zero to yield a singular matrix. In matrix**

*P***, there are three variables**

*P**β*

_{Tn},

*β*

_{MAn},

*β*

_{MUn}, and the singularity condition of matrix

**. Two more conditions are needed to find them with one equation. As mentioned earlier, the cyclic frequency of the system is unique, so recalling the definitions of the wavenumber from Eq. (50) and rewriting them yields**

*P**γ*

_{UT}and

*γ*

_{AT}are defined to introduce the effect of added mass as well as the rigidity changes at sea level and platform level to the natural mode shapes, respectively. To find the coefficient matrix

**for each natural mode, Eq. (60) should be solved by substituting the wavenumbers obtained from the singularity of matrix**

*D***P**. Finally, the natural mode shapes of the system can be found by substituting variables found for each mode in Eqs. (57)–(59). By merging them, it can be represented in a single function as follows:

### 3.2 The Solution for an External Load.

*X*

_{n}(

*z*) is the natural mode of the system, which satisfy

*X*

_{m}(

*z*) and integrating over the length of the tower, it yields

### 3.3 Solution for the Wave Load.

*u*(

*z*,

*t*), which is defined based on the second-order wave theory [72] as follows:

Equation (77) is the wave load based on Morison’s formula rewritten from Eq. (13). *H*_{t}(*z*, *t*) is a step function representing the absolute value function in his formula. To find the response of an OWT under this load, it needs Eq. (70) to be solved after substituting Eq. (77) in it. The resultant will be an ordinary nonhomogeneous second-differential equation. The solution to this differential equation can be found by solving Eq. (70) for *F*_{1}(*z*, *t*) and *H*_{t}(*z*, *t*) *F*_{2}(*z*, *t*) separately and adding them by using the superposition principle. Substituting *F*_{1}(*z*, *t*) as *Q*(*z*, *t*) into the right side of Eq. (70) and solving the integration with respect to the *z* variable analytically will lead to the trigonometrical functions depending on the temporal variable left on the right side of Eq. (70). The solution of it is pretty straightforward. So, the solution for *F*_{1}(*z*, *t*) can be proposed in the form of Eq. (65), where *T*_{n}(*t*) is found by solving Eq. (70) for *F*_{1}(*z*, *t*).

*H*

_{t}(

*z*,

*t*)

*F*

_{2}(

*z*,

*t*), is substituted as

*Q*(

*z*,

*t*) in the right side of it will be challenging. This term consists of

*F*

_{2}(

*z*,

*t*) in Eq. (79) in which temporal and spatial functions are squared and

*H*

_{t}(

*z*,

*t*) in Eq. (80), which is a step function depending on the sign of

*u*(

*z*,

*t*). To obtain a solution for Eq. (70), it needs to work on these two parts to transfer them into the conventional form of functions with a combination of spatial functions and linear trigonometrical terms. To start, the squared term of Eq. (79) is expanded to obtain

*u*(

*z*,

*t*) depends on the values of cos(

*ωt*) and

*f*(

*z*). To evaluate the sign of

*u*(

*z*,

*t*), one needs to find when and where it becomes zero. From Eq. (76), it can be found that

*f*(

*z*) is continuously decreasing when

*z*is increasing because the numerator is always smaller than the denominator. Also, cosh(

*kz*) is always positive and sinh(

*kd*) is positive as long as

*kd*is positive. Therefore,

*f*(

*z*) is always a positive quantity and does not influence the sign of

*u*(

*z*,

*t*). For this reason, the only term that governs the sign of

*u*(

*z*,

*t*) is

*f*(

*z*)cos(

*ωt*) + cos(2

*ωt*). For

*f*(

*z*) > 1, there are two positive roots in the [0, 2

*π*/

*ω*] domain. They are given as follows:

The aforementioned equation represents the value of *H*_{t}(*z*, *t*) in the time domain of [0, 2*π*/*ω*]. For the time domain beyond it, the value of *H*_{t}(*z*, *t*) can be evaluated by considering the time variable relative to a one-period time frame since *F*_{2}(*z*, *t*) is a periodic function.

*T*

_{n}(

*t*) satisfies

*T*

_{1n}(

*t*),

*T*

_{2n}(

*t*), and

*T*

_{3n}(

*t*), where

The solution for the aforementioned equation can be represented as follows:

*H*

_{t}(

*z*,

*t*) poses challenges in solving the integration analytically. To solve the double integration in the aforementioned equation, it requires the removal

*H*

_{t}(

*z*,

*t*) from inside the integrations to have a conventional double integration. As defined in Eq. (88),

*H*

_{t}(

*z*,

*t*) governs the sign of

*v*

_{2n}. It can be either positive or negative, depending on which time frame it is evaluated.

*H*

_{t}(

*z*,

*t*) is positive in the time frame of [0,

*t*

_{1}(

*z*)], negative in the time frame of [

*t*

_{1}(

*z*),

*t*

_{2}(

*z*)], positive in the time frame of [

*t*

_{2}(

*z*),

*t*

_{2}(

*z*) + 2

*t*

_{1}(

*z*)], and so on. By introducing a new integer variable,

*m*, by which

*H*

_{t}(

*z*,

*t*) is positive when

*m*is an odd integer and an even number when it is negative. Therefore,

*H*

_{t}(

*z*,

*t*) can be redefined as follows:

*t*∈ [

*t*

_{e(m−1)}(

*z*),

*t*

_{em}(z)], one can conclude

*t*

_{e0}(z) = 0,

*t*

_{e1}(z) =

*t*

_{1}(

*z*),

*t*

_{e2}(z) =

*t*

_{2}(

*z*),

*t*

_{e3}(z) =

*t*

_{2}(

*z*) + 2

*t*

_{1}(

*z*), and so on. Substituting the aforementioned equation into Eq. (107) and reversing the integration order yields

## 4 Parametric Study and Numerical Example

In this section, the effect of the three solution variables, *α*_{3}, *α*_{5}, and *α*_{6}, introduced in Eqs. (52), (55), and (56) representing foundation rotational stiffness, nacelle rotational mass, and nacelle mass, respectively, as well as the water depth on the natural wavenumbers, and the effect of added mass on the response are investigated.

A numerical example is represented for a reference OWT. The geometry of the system has been chosen from the DTU 10 MW three-bladed OWT presented by Bak et al. [74]. The structural properties are summarized in Table 2. Note that the density is considered approximately 8% more than the regular steel density to take into account the mass of the components such as paint, bolts, flanges, and stiffeners [75]. The average tower diameter is the average of the tower diameter along its height, and the average thickness is calculated from the actual tower mass [57].

Symbol | Value | |
---|---|---|

Tower length (m) | L_{Tow} | 119 |

Tower average diameter (m) | D_{Tow} | 9.6 |

Tower average thickness (m) | t_{Tow} | 0.0295 |

Tower Young's modulus (GPa) | E_{Tow} | 210 |

Nacelle-rotor assembly mass (kg) | M_{n} | 676,723 |

Nacelle-rotor assembly rotational inertia (kg · m^{2}) | J_{p} | 1.7 × 10^{8} |

Platform level from mudline (m) | L_{Plat} | 45 |

Monopile average diameter (m) | D_{Mon} | 8.3 |

Monopile average thickness (m) | t_{Mon} | 0.09 |

Monopile Young's modulus (GPa) | E_{Mon} | 210 |

Material density (kg/m^{3}) | ρ_{s} | 8500 |

Symbol | Value | |
---|---|---|

Tower length (m) | L_{Tow} | 119 |

Tower average diameter (m) | D_{Tow} | 9.6 |

Tower average thickness (m) | t_{Tow} | 0.0295 |

Tower Young's modulus (GPa) | E_{Tow} | 210 |

Nacelle-rotor assembly mass (kg) | M_{n} | 676,723 |

Nacelle-rotor assembly rotational inertia (kg · m^{2}) | J_{p} | 1.7 × 10^{8} |

Platform level from mudline (m) | L_{Plat} | 45 |

Monopile average diameter (m) | D_{Mon} | 8.3 |

Monopile average thickness (m) | t_{Mon} | 0.09 |

Monopile Young's modulus (GPa) | E_{Mon} | 210 |

Material density (kg/m^{3}) | ρ_{s} | 8500 |

Table 3 presents the hydrodynamic loading parameters used in this study. The coefficients are chosen by the recommendations provided by DNV-RP-C205 [72].

Symbol | Value | |
---|---|---|

Water depth (m) | d | 35 |

Drag coefficient | C_{D} | 0.65 |

Added mass coefficient | C_{A} | 1 |

Inertia coefficient | C_{M} | 2 |

Sea water density (kg/m^{3}) | ρ_{w} | 1025 |

Symbol | Value | |
---|---|---|

Water depth (m) | d | 35 |

Drag coefficient | C_{D} | 0.65 |

Added mass coefficient | C_{A} | 1 |

Inertia coefficient | C_{M} | 2 |

Sea water density (kg/m^{3}) | ρ_{w} | 1025 |

The values for the coupled spring model are presented in Table 4 by the work presented by Alkhoury et al. [60]. They calculated these values for the loose sand from the finite element model created for their study in which the same DTU 10 MW OWT is modeled and studied.

Symbol | Value | |
---|---|---|

Lateral stiffness (GN/m) | K_{L} | 2.48 |

Cross stiffness (GN) | K_{LR} | −20.7 |

Rotational stiffness (GN · m) | K_{R} | 412 |

Symbol | Value | |
---|---|---|

Lateral stiffness (GN/m) | K_{L} | 2.48 |

Cross stiffness (GN) | K_{LR} | −20.7 |

Rotational stiffness (GN · m) | K_{R} | 412 |

For this case, the solution variables are calculated and represented in Table 5. By using the procedure described in Sec. 3.1, the natural frequencies of the tower, *f*_{n}, are calculated for the first six modes and represented in Table 4. The corresponding mode shapes are illustrated in Fig. 2, which is to satisfy Eq. (71).

d/l | a_{Tow} | a_{MA} | a_{MU} | α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | α_{6} | γ_{AT} | γ_{UT} |
---|---|---|---|---|---|---|---|---|---|---|---|

0.227 | 12,125 | 14,585 | 7501 | 5.84E-04 | −4.88E-03 | 0.097 | 0.188 | 1.59E-04 | 8.48E-07 | 0.911 | 1.271 |

d/l | a_{Tow} | a_{MA} | a_{MU} | α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | α_{6} | γ_{AT} | γ_{UT} |
---|---|---|---|---|---|---|---|---|---|---|---|

0.227 | 12,125 | 14,585 | 7501 | 5.84E-04 | −4.88E-03 | 0.097 | 0.188 | 1.59E-04 | 8.48E-07 | 0.911 | 1.271 |

In Table 6, the natural frequencies are calculated for both cases of considering the added mass and without the effect of added mass. As represented, including the effect of added mass in the system decreases the natural frequencies. This effect in the first mode is not as significant as in higher modes. The reason can be explained by using Fig. 2 where the displacement of the underwater section in first mode is remarkably less compared to the higher modes. Besides, the presence of the weighty nacelle-rotor assembly mass dominates the motion of the first mode. Therefore, the motion of the system in the higher modes is less than in the first mode due to the substantial inertial force at the top of the system.

Mode number | f (Hz)_{n} | Difference (%) | |
---|---|---|---|

No added mass | Added mass | ||

1 | 0.166561 | 0.166393 | −0.1 |

2 | 1.13463 | 1.0322 | −9.0 |

3 | 2.3888 | 1.98416 | −16.9 |

4 | 4.3686 | 3.8174 | −12.6 |

5 | 8.025 | 6.593 | −17.8 |

6 | 12.198 | 9.8905 | −18.9 |

Mode number | f (Hz)_{n} | Difference (%) | |
---|---|---|---|

No added mass | Added mass | ||

1 | 0.166561 | 0.166393 | −0.1 |

2 | 1.13463 | 1.0322 | −9.0 |

3 | 2.3888 | 1.98416 | −16.9 |

4 | 4.3686 | 3.8174 | −12.6 |

5 | 8.025 | 6.593 | −17.8 |

6 | 12.198 | 9.8905 | −18.9 |

### 4.1 Comparison With the Finite Element Model and the Degree of Accuracy.

To evaluate the accuracy of the results, they are compared with the study conducted by Alkhoury at el. [48,60]. They created a detailed 3D FE model within abaqus/standard to compute the natural frequencies of the DTU 10 MW OWT. They used shell elements to model the tower, including the diameter variation in length and solid elements for monopile. They also fully modeled the soil inside and outside the monopile to investigate the soil structure interaction. They also performed a parametric study on the first natural frequency by varying the water depth and monopile’s diameter and thickness. The first natural frequencies of the system are calculated and compared with the values they calculated for the loose sand that are represented in Table 7. Note that the values of the coupled springs used in this study are also calculated by them which are obtained from the FE model. The differences between the results obtained by the proposed model and the FE model reveal that the proposed model underestimated the first natural frequency for every water depth in the range between 13% and 16.8%. Alkhoury et al. [60] also compared the results of the full 3D model with the one in which the tower cross section is constant for a water depth of 25 m. They found that simplifying modeling by considering the tower’s cross section constant reduced the first natural frequency by 11% for the monopile with 8.3 m in diameter and 9 cm in thickness. The findings of this article also verify this underestimation with a 13.8% deviation. Therefore, this simplification underestimates the first natural frequency that requires using more complicated equations of motion to improve the accuracy of the natural frequency estimation.

Water depth (m) | Platform level (m) | Monopile outer diameter (m) | Monopile thickness (cm) | First natural frequency (Hz) | Deviation (%) | |
---|---|---|---|---|---|---|

Alkhoury et al. [60] | Proposed model | |||||

25 | 35 | 8.3 | 9 | 0.2009 | 0.1731 | −13.8 |

10 | 0.202 | 0.1743 | −13.7 | |||

12 | 0.2056 | 0.1771 | −13.9 | |||

9 | 10 | 0.2074 | 0.177 | −14.7 | ||

12 | 0.2112 | 0.1795 | −15.0 | |||

14 | 0.2158 | 0.1813 | −16.0 | |||

10 | 11 | 0.2177 | 0.1809 | −16.9 | ||

13 | 0.2186 | 0.1828 | −16.4 | |||

15 | 0.2207 | 0.1837 | −16.8 | |||

35 | 45 | 8.3 | 9 | 0.1909 | 0.1663 | −12.9 |

10 | 0.1928 | 0.1682 | −12.8 | |||

12 | 0.197 | 0.171 | −13.2 | |||

9 | 10 | 0.1992 | 0.1719 | −13.7 | ||

12 | 0.2038 | 0.1734 | −14.9 | |||

14 | 0.2092 | 0.176 | −15.9 | |||

10 | 11 | 0.2096 | 0.1768 | −15.6 | ||

13 | 0.2129 | 0.1784 | −16.2 | |||

15 | 0.2155 | 0.1795 | −16.7 | |||

45 | 55 | 8.3 | 9 | 0.1811 | 0.1592 | −12.1 |

10 | 0.1836 | 0.1614 | −12.1 | |||

12 | 0.189 | 0.1648 | −12.8 | |||

9 | 10 | 0.1909 | 0.1658 | −13.1 | ||

12 | 0.1962 | 0.1686 | −14.1 | |||

14 | 0.2023 | 0.1708 | −15.6 | |||

10 | 11 | 0.203 | 0.1717 | −15.4 | ||

13 | 0.207 | 0.1736 | −16.1 | |||

15 | 0.2101 | 0.175 | −16.7 |

Water depth (m) | Platform level (m) | Monopile outer diameter (m) | Monopile thickness (cm) | First natural frequency (Hz) | Deviation (%) | |
---|---|---|---|---|---|---|

Alkhoury et al. [60] | Proposed model | |||||

25 | 35 | 8.3 | 9 | 0.2009 | 0.1731 | −13.8 |

10 | 0.202 | 0.1743 | −13.7 | |||

12 | 0.2056 | 0.1771 | −13.9 | |||

9 | 10 | 0.2074 | 0.177 | −14.7 | ||

12 | 0.2112 | 0.1795 | −15.0 | |||

14 | 0.2158 | 0.1813 | −16.0 | |||

10 | 11 | 0.2177 | 0.1809 | −16.9 | ||

13 | 0.2186 | 0.1828 | −16.4 | |||

15 | 0.2207 | 0.1837 | −16.8 | |||

35 | 45 | 8.3 | 9 | 0.1909 | 0.1663 | −12.9 |

10 | 0.1928 | 0.1682 | −12.8 | |||

12 | 0.197 | 0.171 | −13.2 | |||

9 | 10 | 0.1992 | 0.1719 | −13.7 | ||

12 | 0.2038 | 0.1734 | −14.9 | |||

14 | 0.2092 | 0.176 | −15.9 | |||

10 | 11 | 0.2096 | 0.1768 | −15.6 | ||

13 | 0.2129 | 0.1784 | −16.2 | |||

15 | 0.2155 | 0.1795 | −16.7 | |||

45 | 55 | 8.3 | 9 | 0.1811 | 0.1592 | −12.1 |

10 | 0.1836 | 0.1614 | −12.1 | |||

12 | 0.189 | 0.1648 | −12.8 | |||

9 | 10 | 0.1909 | 0.1658 | −13.1 | ||

12 | 0.1962 | 0.1686 | −14.1 | |||

14 | 0.2023 | 0.1708 | −15.6 | |||

10 | 11 | 0.203 | 0.1717 | −15.4 | ||

13 | 0.207 | 0.1736 | −16.1 | |||

15 | 0.2101 | 0.175 | −16.7 |

### 4.2 Parametric Study on Natural Wavenumbers.

The effect of the water depth and solution parameters, *α*_{3}, *α*_{5}, and *α*_{6}, on the tower wavenumber, *β*_{Tn}, have been parametrically studied for the first five modes, and the results have been presented in Figs. 3–6. In these figures, the tower wavenumbers are normalized to the values of *β*_{Tn} when *α*_{3} = *α*_{5} = *α*_{6} = *d* = 0. In this study, the ratio of water depth and tower length, *d*/*L*, varies from 0 to 1, implying *d* = 0 and *d* = *L*, respectively. The variation of the solution variables as well as their corresponding variation of the parameters used in the parametric study is presented in Table 8. It should be noted that the effect of the support’s lateral stiffness, *K*_{L}, and cross-coupled stiffness, *K*_{LR}, which are used in variables *α*_{1} and *α*_{2}, respectively, are well investigated in the literature in Ref. [56]. So, their effect is not included in the parametric study. Besides, the section of the system is kept constant throughout the length for simplicity.

Solution variable | Equivalent to |
---|---|

α_{1} = ∞ | K_{L} = ∞ |

α_{2} = 0 | K_{LR} = 0 |

0.2 < α_{3} < ∞ | 271 (GN · m) < K_{R} < ∞ |

α_{4} = 1 | EI_{Mon} = EI_{Tow} = 1357 GN · m^{2} |

0 < α_{5} < 10^{−3} | 0 < J_{P} < 1.35 × 10^{9} kg · m^{2} |

0 < α_{6} < 10^{−6} | 0 < M_{n} < 1357 ton |

γ_{AT} = 1 | A_{Mon} = A_{Tow} |

γ_{UT} = 1.36 | See Eq. (62) |

Solution variable | Equivalent to |
---|---|

α_{1} = ∞ | K_{L} = ∞ |

α_{2} = 0 | K_{LR} = 0 |

0.2 < α_{3} < ∞ | 271 (GN · m) < K_{R} < ∞ |

α_{4} = 1 | EI_{Mon} = EI_{Tow} = 1357 GN · m^{2} |

0 < α_{5} < 10^{−3} | 0 < J_{P} < 1.35 × 10^{9} kg · m^{2} |

0 < α_{6} < 10^{−6} | 0 < M_{n} < 1357 ton |

γ_{AT} = 1 | A_{Mon} = A_{Tow} |

γ_{UT} = 1.36 | See Eq. (62) |

#### 4.2.1 The Effect of Water Depth.

As illustrated in Fig. 3, the wavenumber for the first five modes decreases by increasing the water depth regardless of the boundary conditions properties because of the presence of added mass to the system for the case of the properties introduced in Table 8. The value of the wavenumber in the first mode remains almost constant by increasing the water depth up to 0.4*L* and drops up to *L* for the case of *α*_{6} = 0, while the variation is almost constant for the case when *α*_{6} = 10^{−6}. This can be because of the domination of the heavy mass in the motion of the system in the first mode. The decrease of the wavenumber for the second mode starts at 0.2*L* and decreases approximately the same amount as the first mode at the sea level equal to the tower length for the cases when *α*_{6} = 0. The initiation of the drop of wavenumber for higher modes is almost half of the previous modes. It can be concluded that the effect of shallow water compared to the tower length and consequently the added mass in the lower modes is not significant as opposed to higher modes where it drops immediately by increasing the sea level. This phenomenon may be important by the fact that the free sea level varies in each wave period. Therefore, the water surface variation can significantly change higher natural wavenumbers in the shallow water proportional to the tower length for each period of wave load. However, for higher values of *d*/*L*, the lower natural modes are also influenced by sea level variation. This variation of the natural frequencies is also reported in the literature for the first natural modes based on the measured data in Refs. [61,65].

The pattern of the variation of wavenumbers by varying the water depth shown in Fig. 3 reveal a wavy-shape decrease in which the reduction rate changes in different mode numbers. It can be seen that the number of crests in this pattern is equal to the mode number. For instance, for the second mode, two crests at around *d*/*L* = 0.2 and *d*/*L* = 0.8 can be seen, while five distinguished crests are visible in the figure for mode 5. Therefore, it can be concluded that the variation of wavenumbers versus *d*/*L* is converging to a linear reduction rate.

#### 4.2.2 The Effect of Support Rotational Stiffness, α_{3}.

Figure 4 illustrates the variation of the wavenumber against the different values of *α*_{3} in first and fifth modes at *d*/*L* = 0.227, for instance. As expected, by increasing *α*_{3}, or decreasing the support’s rotational stiffness, the value of the wavenumber decreases linearly in the first mode and nonlinearly in the fifth mode. By increasing the rotational softness of the support, the wavenumbers for all natural modes decrease. This is because the higher rotational softness provides higher rotation in the support resulting in a reduction of the wavenumbers. In the first mode, the effect of *α*_{3} is more than in the fifth mode. Also, the impact of the *α*_{3} is less in higher values of the *α*_{5} and *α*_{6}.

#### 4.2.3 The Effect of Nacelle-Rotor Assembly Rotational Moment of Inertia, α_{5}.

The parametric study on *α*_{5} has been illustrated in Fig. 5 for the first and fifth modes for two values of *α*_{3} = ∞ and 0.2 and *α*_{3} = 0 and 10^{−6} for *d*/*L* = 0.6 when the values of *α*_{5} varying between 0 and 10^{–3}. As shown in Fig. 5(a), the wavenumber decreases by increasing *α*_{5} for the first mode. For the sixth mode, the natural wavenumber is almost insensitive to the variation of *α*_{5} despite a rapid decrease of natural wavenumber at small values of *α*_{5}. The decrease of wavenumber can be explained by the fact that the rotational moment of inertia at the top of the tower increases the mass momentum of the system at the top. Therefore, the moment of inertia due to the nacelle-rotor assembly causes the system to oscillate slower, yielding to the lower values of the natural wavenumbers.

#### 4.2.4 The Effect of the Nacelle Mass, α_{6}.

Figure 6 illustrates the variation of the wavenumber against *α*_{6} ranging from 0 to 10^{−6} for two values of *α*_{3} = ∞ and 0.2 and *α*_{5} = 0 and 10^{−3} in first and fifth modes when *d*/*L* = 0.6. A reduction is shown in Fig. 6(a) for all cases of *α*_{3} and *α*_{5} in the first mode. By increasing *α*_{6}, which is the increase of the top mass with respect to the tower’s flexural rigidity, the motion of the tower becomes slower yielding to the smaller values of the natural wavenumbers.

The reduction effect of *α*_{6}, which is proportional to the nacelle-rotor assembly mass, can be explained by the whipping effect of the tower. The inertia of the heavy mass at the top of the tower may cause a delayed motion relative to the mid-section of the tower in the same direction. When the mid-section of the tower reaches its maximum displacement, the top section is still moving imposing extra shear force to the mid-section pushing it to move further, consequently increasing the oscillation period, decreasing the frequency, and decreasing the wavenumber of the system, see Eq. (50). Physically speaking, the heavy mass at the top of the tower produces an inertia force in the opposite direction of the motion, which slows down the motion of the tower.

In higher modes, as illustrated in Fig. 6, the variation of the natural wavenumber of the system is a smooth reduction despite a rapid reduction in the small values of *α*_{6}. The reason for that can be explained by the fact that the translational acceleration at the top of the tower in higher modes is relatively smaller than that of the lower modes due to the whipping effect explained earlier. This makes higher modes less sensitive to the variation of the nacelle mass. Furthermore, this effect is shown in Fig. 3 where the curves are gathering together by increasing the mode number implying that the effect of the boundary conditions at both sides of the system is fading out of the natural wavenumbers.

### 4.3 The Response of the Reference Tower to the Wave Load.

The wave load introduced in Eq. (89) is applied to the system, and the response is evaluated by Eq. (97). The properties of Stokes’s wave kinematics applied to Morison’s formula are briefly represented. The proposed model's ability to deal with the difficulties posed by the drag term is explained. Besides, the effect of the added mass on the response is investigated for a certain wave load configuration. Finally, a comparison is performed between the responses obtained by the proposed solution and the numerical one.

#### 4.3.1 The Application of the Wave Load in the Proposed Solution.

The inherent properties of wave load based on Morison’s formula with Stokes's wave kinematics raise some difficulties in evaluating the response of the tower. Before presenting how the proposed solution deals with these properties, one needs to discuss the properties of the wave load with Morison’s equation introduced in Eq. (89). To this end, an ocean wave with a height of 5.1 m and a length of 132 m is chosen. This ocean wave represents a normal sea state of a water depth of 35 m in a wind speed of 26 m/s, which is reported as a nonlinear ocean wave state in Ref. [76]. The properties of this ocean wave are presented in Table 9. The wave horizontal particle velocity, *u*(*z*, *t*), introduced in Eq. (72), and the corresponding ocean wave load are drawn in Figs. 7(a) and 7(b), respectively. As expected, the wave horizontal particle velocity is a periodic but nonsymmetric function, which is the property of the second-order wave theory. The wave load is also not started from zero as shown in Fig. 7 because of the presence of the cosine function in the drag term of the wave load formula.

Height H (m) | Ocean wavelength λ (m) | Water depth d (m) | Frequency f (Hz) | Period | $HgT2$ | $dgT2$ | $Hd$ | Ursell number Ur |
---|---|---|---|---|---|---|---|---|

5.1 | 132 | 35 | 0.105 | 9.53 | 0.00572 | 0.0392 | 0.145 | 2.072 |

Height H (m) | Ocean wavelength λ (m) | Water depth d (m) | Frequency f (Hz) | Period | $HgT2$ | $dgT2$ | $Hd$ | Ursell number Ur |
---|---|---|---|---|---|---|---|---|

5.1 | 132 | 35 | 0.105 | 9.53 | 0.00572 | 0.0392 | 0.145 | 2.072 |

Moreover, Fig. 7(a) reveals that the wave horizontal particle velocity, *u*(*z*, *t*), does not become zero at the same time for all values of *z*. This is the second property of the second-order wave theory that causes difficulties in obtaining the response of a system loaded with it. Therefore, the situation in which *u*(*z*, *t*) = 0, depends on temporal and special variables. This is shown more precisely in Fig. 8 for the variation of *u*(*z*, *t*) in the z-direction for some instant of time around the first zero value. Therefore, two functions of *t*_{1}(*z*) and *t*_{2}(*z*) represented by Eqs. (86) and (87) are defined to evaluate the time when *u*(*z*, *t*) is zero. The importance of defining these two functions is to evaluate the absolute value function in the drag term of the wave load, *u*(*z*, *t*)|*u*(*z*, *t*)|, in Eq. (13). To compensate for the absolute value function, *H*_{t}(*z*, *t*) is defined as a function *t*_{1}(*z*) and *t*_{2}(*z*) in Eq. (80). This poses difficulties in evaluating Eqs. (100) and (101). By defining *H*_{t}(*z*, *t*) as Eq. (109) and carrying out the integration, it becomes possible to remove *H*_{t}(*z*, *t*). In addition, changing the integration order by which the temporal integration is taken first in the domain as a function of *t*_{1}(*z*) and *t*_{2}(*z*) results in trigonometrical functions. By introducing *t*_{1}(*z*) and *t*_{2}(*z*) and their combinations in periods, the special outer integration can be evaluated analytically since those are functions of reversed trigonometrical functions as shown in Eqs. (86) and (87). Therefore, the response can be obtained as a function without any need for numerical evaluation.

The presence of *H*_{t}(*z*, *t*) in the wave load is quite essential. To show this, an imaginary wave load based on Eq. (77) is defined as *F*_{1}(*z*, *t*) + *F*_{2}(*z*, *t*). By drawing it together with the wave load in Eq. (77), the effect of the *H*_{t}(*z*, *t*) reveals itself. Figure 9 illustrates this comparison at the sea level, *z* = 35. It can be seen that *H*_{t}(*z*, *t*) causes significant changes in the wave load when *u*(*z*, *t*) is negative between *t*_{1}(*z*) and *t*_{2}(*z*). Therefore, its presence cannot be neglected in evaluating the response. However, the severity of this effect may be different in other wave configurations.

#### 4.3.2 The Response of the Reference Tower to the Second Wave Load.

The response history of the system is illustrated in Fig. 10 at different heights in which the first five natural modes are participating. The selected wave height and length correspond to a wave frequency of 0.1 Hz, *T* = 9.52 s, which is way below the first natural frequency of the reference tower, 0.1663 Hz, represented in Table 9. As shown in Fig. 10, the maximum deflection occurs at the top of the towers reaching 0.13 m. In addition, the time history deflection curve in Fig. 10 reveals that the responses are a nonperiodic vibration even though the loading is periodic. This can be because of the indirect presence of the term *H*_{t}(*z*, *t*) in the response in Eqs. (110) and (111). The indirect presence of *H*_{t}(*z*, *t*) in the response shows itself by being displayed in different time frames when the direction of the wave horizontal particle velocity changes. This triggers the transient responses at the beginning of each stage causing the response to be nonperiodic. Moreover, Fig. 10 shows that the natural frequencies originating from the transient response are carried by the steady-state response. This also causes the response to be nonperiodic.

The deflection of the reference tower at the early stages of motion is illustrated in Fig. 11. The motion starts from the rest initial condition and follows by imposing the wave load up to sea level, in this case, 35 m, causing the lower sections of the tower to move, while the movement of the upper section is delayed because of the nacelle mass and tower softness. This can be referred to as the *whipping effect* when the motion of the upper sections is amplified by the motion of the lower sections of the tower. This effect can be seen clearly in the motion of the reference tower, which has low stiffness or heavy nacelle mass at the top. Moreover, by looking at the deflection of the reference tower in Fig. 11 for time instants from *t*/*T* = 0.08 up to 0.1, *T* being the period of the wave which is 9.6 s, it implies that the motion of the tower in low sections is slowing down, while the upper section is still moving toward the negative deflections representing an instance for the whipping effect. It should be mentioned that the wave load at zero time is a positive value, as shown in Fig. 7(b), decreasing in the early stages of the loading and reversing its direction as time passes.

#### 4.3.3 The Effect of Added Mass on the Response.

The inertia term of Morison’s formula adds an extra mass to the system up to the sea level. To represent its effect on the response of the reference tower, it is evaluated by considering *C*_{A} = 1 and *C*_{A} = 0 to simulate a system with and without the presence of the added mass, respectively, and the time history responses are illustrated in Fig. 12 at the nacelle level. As mentioned in Sec. 4.2.1 and illustrated in Fig. 3, the value of the system’s natural wavenumbers decreases by increasing the sea level, resulting from the decrease in natural frequencies based on Eq. (50). This leads the system with added mass oscillating with higher natural periods falling forward than the one without added mass, as shown in Fig. 12. Furthermore, as presented in Table 6, the added mass is significantly influenced in the second and higher natural frequencies in the case of the system of this study. The differences between the natural frequencies in higher modes are shown in Fig. 12.

The effect of the added mass up to the sea level can also be seen in the early stages of the loading in Fig. 13. The added mass increases the inertia force of the system. The higher the inertia force is, the slower the vibration results. Therefore, the system's motion with added mass included delays compared to the one without the added mass as shown in Fig. 13. It is also revealing that the added mass is appended to the system up to sea level, providing lower deflection at the lower sections of the system, while the upper sections have almost the same deflection. The overall interpretation from Fig. 13 reveals that the effect of the added mass up to the sea level is successfully simulated by the proposed solution.

#### 4.3.4 Comparison With Numerical Results.

The derived formulation of the response of the system is based on the expansion in the natural modes. The accuracy of the solution depends on the number of natural modes participating in the solution. For a continuous system, an infinite number of modes is expected. The higher the number of modes participating in the solution, the more accurate the response obtained. It is worthwhile to mention that the Bernoulli–Euler beam theory in the form of Eq. (1) is applicable for the lower natural modes only. In higher modes, the shearing deformation and rotatory inertia significantly affect the natural frequencies [50]. Therefore, reaching higher natural modes requires the equation of motion with those effects considered. However, as far as the wave load by using Morison’s formula matters, its validity requires that the ocean wavelength should be five times [72] higher than the diameter of the slender structure. This limitation causes the frequency of the ocean wave load to be lower than the system’s first natural frequency. Therefore, only the first few natural modes may be enough to evaluate the tower's response.

To confirm the accuracy of the proposed solution, a numerical evaluation is performed by using the standard commercial software mathematica^{®} [77]. The partial differential numerical solver, *NDsolve*, is chosen to solve Eq. (8) by introducing the wave velocity from Eq. (72). The properties of the system chosen for the numerical evaluation are presented in Table 10.

Symbol | Value | |
---|---|---|

Structural properties | ||

Tower length (m) | L_{Tow} | 70 |

Tower average diameter (m) | D_{Tow} | 6 |

Tower average thickness (m) | t_{Tow} | 0.045 |

Tower Young's modulus (GPa) | E_{Tow} | 210 |

Nacelle-rotor assembly mass (kg) | M_{n} | 0 |

Nacelle-rotor assembly rotational inertia (kg · m^{2}) | J_{p} | 0 |

Platform level from mudline (m) | L_{Plat} | 45 |

Monopile average diameter (m) | D_{Mon} | 6 |

Monopile average thickness (m) | t_{Mon} | 0.045 |

Monopile Young's modulus (GPa) | E_{Mon} | 210 |

Material density (kg/m^{3}) | ρ_{s} | 7820 |

Support stiffness | ||

Lateral stiffness (GN/m) | K_{L} | ∞ |

Cross stiffness (GN) | K_{LR} | 0 |

Rotational stiffness (GN · m) | K_{R} | 205.72 |

Hydrodynamic loading properties | ||

Water depth (m) | d | 30 |

Drag coefficient | C_{D} | 0.65 |

Added mass coefficient | C_{A} | 1 |

Inertia coefficient | C_{M} | 2 |

Sea water density (kg/m^{3}) | ρ_{w} | 1020 |

Symbol | Value | |
---|---|---|

Structural properties | ||

Tower length (m) | L_{Tow} | 70 |

Tower average diameter (m) | D_{Tow} | 6 |

Tower average thickness (m) | t_{Tow} | 0.045 |

Tower Young's modulus (GPa) | E_{Tow} | 210 |

Nacelle-rotor assembly mass (kg) | M_{n} | 0 |

Nacelle-rotor assembly rotational inertia (kg · m^{2}) | J_{p} | 0 |

Platform level from mudline (m) | L_{Plat} | 45 |

Monopile average diameter (m) | D_{Mon} | 6 |

Monopile average thickness (m) | t_{Mon} | 0.045 |

Monopile Young's modulus (GPa) | E_{Mon} | 210 |

Material density (kg/m^{3}) | ρ_{s} | 7820 |

Support stiffness | ||

Lateral stiffness (GN/m) | K_{L} | ∞ |

Cross stiffness (GN) | K_{LR} | 0 |

Rotational stiffness (GN · m) | K_{R} | 205.72 |

Hydrodynamic loading properties | ||

Water depth (m) | d | 30 |

Drag coefficient | C_{D} | 0.65 |

Added mass coefficient | C_{A} | 1 |

Inertia coefficient | C_{M} | 2 |

Sea water density (kg/m^{3}) | ρ_{w} | 1020 |

The wave height, *H*, and length, *λ*, are selected to be 3 and 100 m, respectively. By setting *MaxStepSize* equal to 1.4, and *AccuracyGoal* and *PrecisionGoal* to 6, *NDsolve* solves the partial differential equation (PDE) by using the Hermite method in orders of 7 and 3 in *z* and *t* variables, respectively, in the domain of *z* ∈ [0, 115] and *t* ∈ [0, 34]. Besides, the first six natural modes of the tower with the properties and loading the same as the numerical one are selected for participating in the proposed solution results. The comparison between the response obtained by the numerical solution and the proposed one verifies the perfect agreement between the two methods, as illustrated in Fig. 14.

## 5 Conclusions

An analytical solution for the modal analysis of offshore wind turbine structures has been developed. The solution includes the wave–structure interaction by appending an extra mass to the system underwater. In an effort to propose a more accurate solution based on the classical analytical methods, the flexibility of the foundation as well as the inertial forces induced by the nacelle-rotor assembly translational and rotational inertia are assigned to the boundaries of the system. Besides, the considerable cross-sectional changes at the platform level where the monopile is connected to the tower by a transition part are taken into account in the solution. Overall, a system of three partial differential equations consisting of 12 boundary conditions and two initial conditions has been solved using the expansion theorem.

The effect of water depth, foundation rotational flexibility, nacelle mass, and nacelle-blades rotational inertia on the system’s natural wavenumber were studied parametrically for the first five natural modes of the system. The results reveal the following:

The natural wavenumber decreases by increasing the water depth to the tower–length ratio,

*d*/*L*, for all natural modes producing a wavy pattern based on the modal number. The effect of the foundation rotational flexibility, nacelle-blades rotational inertia, and nacelle mass on the natural wavenumber decreases by increasing the modal number for all sea level values. More importantly, the variation of the natural wavenumbers by variation of the sea level implies that the natural wavenumbers or natural frequencies of the system vary during sea level variation during a period of wave load. Therefore, the system's natural frequencies can be considered a time-dependent quantity. It can be essential in the assessment of the ringing-type resonance of the system and fatigue loading.The system’s natural frequencies decrease by increasing the foundation rotational flexibility, mass, and the rotational inertial of the nacelle-rotor assembly. This pattern has been seen at all water depths.

The proposed model, based on the simplification of considering the constant cross section for the tower, underestimates the first natural frequency of the system between 13% and 16.8%. Reaching higher accuracy requires establishing more complicated equations of motions by accounting for the cross-sectional variation of the tower. However, the proposed method is straightforward and agile in calculating cost-efficient natural frequencies.

The solution for the undamped response of the tower under the wave load with the second-order Stokes’s wave theory based on Morison’s formula has been developed as an analytical function. Two major contributions are as follows:

The drag term of Morison’s formula, neglected by many researchers, is successfully included in finding the response of the system by defining

*H*_{t}(*z*,*t*) to remove the absolute value function in the drag term.The comparison made between the responses of the system with and without added mass showed that the presence of added mass up to the sea level changes the shape of the response. This also can be important in fatigue evaluation of the system by providing a more realistic estimation of the stress status in the structure.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

No data, models, or code were generated or used for this paper.