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 [47]. Reviews of advances in vibration control algorithms for smart structures up to 2017 are presented by Refs. [810]. 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,1319] 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 [3537], 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 [3848] 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,5355] 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 KL, KR, KLR, andKz, respectively [71].

Fig. 1
(a) A typical OWT configuration, (b) schematic of the model, and (c) wave loading direction and system coordinates
Fig. 1
(a) A typical OWT configuration, (b) schematic of the model, and (c) wave loading direction and system coordinates
Close modal
Table 1

Symbols definition

SymbolStructural properties
LNacelle level from the seabed (m)
LTowTower length (m)
DTowTower average diameter (m)
tTowTower average thickness (m)
mTowTower mass of unit length (kg/m)
ETowTower Young's modulus (GPa)
EITowFlexural rigidity of the tower, i.e., ETowITow (GPa · m4)
MnNacelle-Rotor assembly mass (kg)
JpNacelle-Rotor assembly rotational inertia (kg · m2)
LPlatPlatform level from the seabed (m)
DMonMonopile average diameter (m)
tMonMonopile average thickness (m)
AMonMonopile cross-sectional area (m2)
mMonMonopile mass of unit length (kg/m)
EMonMonopile Young's modulus (GPa)
EIMonFlexural rigidity of the monopile, i.e., EMonIMon (GPa · m4)
ρsMaterial density (kg/m3)
SymbolSupport stiffness
KLLateral stiffness (GN/m)
KLRCross stiffness (GN)
KRRotational stiffness (GN · m)
SymbolHydrodynamic loading properties
dWater depth (m)
CDDrag coefficient
CAAdded mass coefficient
CMInertia coefficient
ρwSea water density (kg/m3)
λOcean wavelength (m)
HOcean wave height (m)
SymbolStructural properties
LNacelle level from the seabed (m)
LTowTower length (m)
DTowTower average diameter (m)
tTowTower average thickness (m)
mTowTower mass of unit length (kg/m)
ETowTower Young's modulus (GPa)
EITowFlexural rigidity of the tower, i.e., ETowITow (GPa · m4)
MnNacelle-Rotor assembly mass (kg)
JpNacelle-Rotor assembly rotational inertia (kg · m2)
LPlatPlatform level from the seabed (m)
DMonMonopile average diameter (m)
tMonMonopile average thickness (m)
AMonMonopile cross-sectional area (m2)
mMonMonopile mass of unit length (kg/m)
EMonMonopile Young's modulus (GPa)
EIMonFlexural rigidity of the monopile, i.e., EMonIMon (GPa · m4)
ρsMaterial density (kg/m3)
SymbolSupport stiffness
KLLateral stiffness (GN/m)
KLRCross stiffness (GN)
KRRotational stiffness (GN · m)
SymbolHydrodynamic loading properties
dWater depth (m)
CDDrag coefficient
CAAdded mass coefficient
CMInertia coefficient
ρwSea water density (kg/m3)
λOcean wavelength (m)
HOcean 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.

For the system introduced in Fig. 1(c), the Bernoulli–Euler beam theory can be applicable for small displacements [49]. Therefore, the equation of motion along the height of the system with the origin from the seabed can be written as follows:
EI(z)x(iv)(z,t)+M(z)x¨(z,t)=q(z,t)H3(z)
(1)
where 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:
M(z)=mTowH1(z)+mMonH2(z)+mMonH3(z)
(2)
EI(z)=EITowH1(z)+EIMonH2(z)+EIMonH3(z)
(3)
where Hi(z), i = 1, 2, 3 are the step functions defined as follows:
H1(z)={1,forLPlatzL0,ford<zLPlat0,for0<zd
(4)
H2(z)={0,forLPlatzL1,ford<zLPlat0,for0<zd
(5)
H3(z)={0,forLPlatzL0,ford<zLPlat1,for0<zd
(6)
For a slender structure, Morison’s equation from Refs. [72,73] can be adopted for dynamic modeling. Also, the relative velocity formulation is applicable for a moving slender structure subjected to the wave loads. Since the bottom-fixed support structure is of interest, it is expected that the underwater motion of the system to be way below its diameter. Therefore, the relative velocity in the drag term can be reduced to the wave horizontal particle velocity. For an ocean wavelength larger than five times, the monopile diameter, and the small displacement, the wave load on the monopile can be represented by the relative velocity formulation [72] as follows:
q(z,t)=ρwCAAMonx¨(z,t)+ρwCMAMonDu(z,t)Dt+12ρwCDDMonu(z,t)|u(z,t)|
(7)
In Eq. (7), the two first terms inside the curl bracket are inertial and the third one is the drag term. The total derivative of the wave horizontal particle velocity, Du(z, t)/Dt = ∂u/∂t + uu/∂x + wu/∂z, in the inertial term can be reduced to u˙(z,t)=u(z,t)/t by neglecting the advocative terms, which are reported to slightly increase the load when they are included [70]. The inertia term reveals that the ρwCAAMonx¨(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
x(iv)(z,t)+A(z)x¨(z,t)=1EIMonQ(z,t)H3(z)
(8)
where
A(z)=1aTow2H1(z)+1aMA2H2(z)+1aMU2H3(z)
(9)
aTow2=EITowmTow
(10)
aMA2=EIMonmMon
(11)
aMU2=EIMonmMon+ρwCAAMon
(12)
Q(z,t)=ρwCMAMonu˙(z,t)+12ρwCDDMonu(z,t)|u(z,t)|
(13)

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 aMU is the mass of the system underwater including added mass. Therefore, the wave–structure interaction is included in the equation of motion.

It can be seen in Eq. (8) that there are three separate systems, the tower, monopile above water, and monopile underwater, acting together to govern the motion of the system. Therefore, it can be separated in the form of three independent equations as follows:
xTow(iv)(z,t)+1aTow2x¨Tow(z,t)=0,forLPlatzL
(14)
xMA(iv)(z,t)+1aMA2x¨MA(z,t)=0,ford<zLPlat
(15)
xMU(iv)(z,t)+1aMU2x¨MU(z,t)=1EIMonQ(z,t),for0<zd
(16)

In the aforementioned three equations, xTow(z, t), xMA(z, t), and xMU(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.

To accommodate the motion of the system with the equation of motion in Eq. (8), two sets of conditions at the two ends of the system are needed. The motion of the system at the seabed is governed by a set of four springs as illustrated in Fig. 1(b). Vertical stiffness, Kz, 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]:
[F(t)M(t)]=[KLKLRKLRKR][x(0,t)x(0,t)]
(17)
By substituting the shear force and bending moment from the beam theory into the aforementioned equation and expanding, it yields the following two boundary conditions:
EIMonxMU(0,t)=KLxMU(0,t)+KLRxMU(0,t)
(18)
EIMonxMU(0,t)=KLRxMU(0,t)+KRxMU(0,t)
(19)
At the top of the tower, a heavy nacelle is mounted providing lump mass to the system. The effect of the translational and rotational inertia of the mass of the nacelle can be simulated as the two following conditions:
EITowxTow(L,t)=JPx¨Tow(L,t)
(20)
EITowxTow(L,t)=MNx¨Tow(L,t)
(21)

The rotational motion of the nacelle produces the momentum proportional to the tower’s rotational acceleration x¨Tow(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¨Tow(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).

As mentioned earlier, Eqs. (15) and (16) govern the motion of the monopile. Therefore, it is necessary that these equations are linked together at the sea level via some boundary conditions. They are as follows:
xMU(d,t)=xMA(d,t)
(22)
xMU(d,t)=xMA(d,t)
(23)
xMU(d,t)=xMA(d,t)
(24)
xMU(d,t)=xMA(d,t)
(25)
Also, Eqs. (14) and (15) are connected by following boundary conditions:
xMA(LPlat,t)=xTow(LPlat,t)
(26)
xMA(LPlat,t)=xTow(LPlat,t)
(27)
EIMonxMA(LPlat,t)=EITowxTow(LPlat,t)
(28)
EIMonxMA(LPlat,t)=EITowxTow(LPlat,t)
(29)

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.

It is assumed that the tower motion starts from the position when it is at the rest. Therefore, the initial deflection and the velocity of each part of the tower are zero. It yields the following initial conditions:
x(z,0)=x˙(z,0)=0
(30)

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.

The natural modes will be evaluated by using the homogenous form of Eq. (16) as well as Eqs. (14) and (15), which are expressed as follows:
xTn(iv)(z,t)+1aTow2x¨Tn(z,t)=0,forLPlatzL
(31)
xMAn(iv)(z,t)+1aMA2x¨MAn(z,t)=0,ford<zLPlat
(32)
xMUn(iv)(z,t)+1aMU2x¨MUn(z,t)=0,for0<zd
(33)
It should be noted that the subscript n in quantity indicates that it belongs to the nth natural mode. So, xTn(z, t), xMAn(z, t), and xMUn(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:
xTn(z,t)=XTn(z)TTn(t),forLPlatzL
(34)
xMAn(z,t)=XMAn(z)TMAn(t),ford<zLPlat
(35)
xMUn(z,t)=XMUn(z)TMUn(t),for0<zd
(36)
Substituting Eqs. (34)(36) into Eqs. (31)(33) will yield
XTn(iv)(z)TTn(t)+1aTow2XTn(z)T¨Tn(t)=0,forLPlatzL
(37)
XMAn(iv)(z)TMAn(t)+1aMA2XMAn(z)T¨MAn(t)=0,ford<zLPlat
(38)
XMUn(iv)(z)TMUn(t)+1aMU2XMUn(z)T¨MUn(t)=0,for0<zd
(39)
After some algebraic manipulations, they are transformed to
XTn(iv)(z)XTn(z)=βTn4=1aTow2T¨Tn(t)TTn(t),forLPlatzL
(40)
XMAn(iv)(z)XMAn(z)=βMAn4=1aMA2T¨MAn(t)TMAn(t),ford<zLPlat
(41)
XMUn(iv)(z)XMUn(z)=βMUn4=1aMU2T¨MUn(t)TMUn(t),for0<zd
(42)
where β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:
XTniv(z)βTn4XTn(z)=0,forLPlatzL
(43)
XMAniv(z)βMAn4XMAn(z)=0,ford<zLPlat
(44)
XMUniv(z)βMUn4XMUn(z)=0,for0<zd
(45)
and
T¨Tn(t)+aTow2βTn4TTn(t)=0,forLPlatzL
(46)
T¨MUn(t)+aMU2βMUn4TMUn(t)=0,ford<zLPlat
(47)
T¨MAn(t)+aMA2βMAn4TMAn(t)=0,for0<zd
(48)
In the aforementioned six equations, the temporal and spatial variables are separated. Therefore, they can be solved independently. Since every section of a continuous system should vibrate with the same natural frequency in each mode shape, so TTn = TMAn(t) = TMUn(t). By comparing Eqs. (46) and (47), one can conclude that aTow2βTn4=aMU2βMUn4=aMA2βMAn4 because TTn = TMAn(t) = TMUn(t). Therefore,
T¨n(t)+ωn2Tn(t)=0,for0zL
(49)
where
ωn2=aTow2βTn4=aMU2βMUn4=aMA2βMAn4
(50)

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

The solutions of the motion Eqs. (34)(36) can be imposed into the boundary conditions defined by Eqs. (18)(29) to yield the boundary conditions independent from the time variable. They are
XMUn(0)+α1XMUn(0)+α2XMUn(0)=0,α1=KLEIMon,α2=KLREIMon
(51)
XMUn(0)α2XMUn(0)α3XMUn(0)=0,α3=KREIMon
(52)
XMUn(d)XMAn(d)=0XMUn(d)XMAn(d)=0XMUn(d)XMAn(d)=0XMUn(d)XMAn(d)=0
(53)
XMAn(LPlat)XTn(LPlat)=0XMAn(LPlat)XTn(LPlat)=0XMAn(LPlat)α4XTn(LPlat)=0XMAn(LPlat)α4XTn(LPlat)=0,α4=EITowEIMon
(54)
XTn(L)α5ωn2XTn(L)=0,α5=JpEITow
(55)
XTn(L)+α6ωn2XTn(L)=0,α6=MnEITow
(56)
where α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.
The solution for Eqs. (43)(45) are in the form of
XTn(z)=T1cos(βTnz)+T2cosh(βTnz)+T3sin(βTnz)+T4sinh(βTnz)
(57)
XMAn(z)=A1cos(βMAnz)+A2cosh(βMAnz)+A3sin(βMAnz)+A4sinh(βMAnz)
(58)
XMUn(z)=U1cos(βMUnz)+U2cosh(βMUnz)+U3sin(βMUnz)+U4sinh(βMUnz)
(59)
where Ti, Ai, and Ui 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:
P×D=0,D={U1U2U3U4A1A2A3A4T1T2T3T4}T
(60)
where P is the matrix containing trigonometrical and hyper trigonometrical functions and D 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 P should be zero to yield a singular matrix. In matrix P, there are three variables βTn, βMAn, βMUn, and the singularity condition of matrix P. 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
βMAn=γATβTn,γAT=α4mMonmTow4
(61)
βMUn=γUTβTn,γUT=α4(mMon+ρwCAAMonmTow)4
(62)
where γ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 D for each natural mode, Eq. (60) should be solved by substituting the wavenumbers obtained from the singularity of matrix 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:
Xn(z)=XT(z)H1(z)+XMA(z)H2(z)+XMU(z)H3(z),0zL
(63)

3.2 The Solution for an External Load.

The solution of an equation of motion in the form of
x(iv)(z,t)+A(z)x¨(z,t)=1EIMonH3(z)Q(z,t)
(64)
can be obtained by using the expansion theorem to represent the motion of the tower in the form of
x(z,t)=n=1Xn(z)Tn(t)
(65)
where the Xn(z) is the natural mode of the system, which satisfy
Xn(iv)ωn2A(z)Xn(z)=0
(66)
The aforementioned equation is obtained by substituting Eq. (50) into Eqs. (43) and (44) and merging them by using Eq. (9). Substituting Eq. (65) into Eq. (64) yields
n=1Xn(iv)(z)Tn(t)+A(z)n=1Xn(z)T¨n(t)=1EIMonH3(z)Q(z,t)
(67)
Multiplying both sides by Xm(z) and integrating over the length of the tower, it yields
n=1Tn(t)0LXm(z)Xn(iv)(z)dz+n=1T¨n(t)0LA(z)Xm(z)Xn(z)dz=1EIMon0LH3(z)Q(z,t)Xm(z)dz
(68)
Substituting Xn(iv) from Eq. (66) into Eq. (68) results in
n=10LA(z)Xm(z)Xn(z)dz(T¨n(t)+ωn2Tn(t))=1EIMon0LH3(z)Q(z,t)Xm(z)dz
(69)
The natural modes are orthogonal and normalized. Therefore, the aforementioned equation is simplified to
T¨n(t)+ωn2Tn(t)=1EIMon0LH3(z)Q(z,t)Xm(z)dz
(70)
It should be noted that the natural modes are normalized such that
0LA(z)Xn2(z)dz=1
(71)

By finding Tn(t) from solving Eq. (70), it can be substituted into Eq. (65) to obtain the response of the tower.

3.3 Solution for the Wave Load.

The wave load acting on the tower was introduced in Sec. 2.1. Eq. (8) governs the motion of the tower with the external load in the form of Eq. (13). It is a function of horizontal particle wave velocity, u(z, t), which is defined based on the second-order wave theory [72] as follows:
u(z,t)=f1cosh(kz)cos(ωt)+f2cosh(2kz)cos(2ωt)
(72)
where
f1=ωH2sinh(kd)
(73)
f2=316ωkH2sinh4(kd)
(74)
Equation (72) can be rewritten as follows:
u(z,t)=f2cosh(2kz)(f(z)cos(ωt)+cos(2ωt))
(75)
where
f(z)=83sinh3(kd)kcosh(kz)cosh(2kz)
(76)
Substituting Eq. (75) into Eq. (13) yields
Q(z,t)=F1(z,t)+Ht(z,t)F2(z,t)
(77)
where
F1(z,t)=ρwCMAωf2cosh(2kz)(f(z)sin(ωt)+2sin(2ωt))
(78)
F2(z,t)=12ρwCDDf22cosh2(2kz)(f(z)cos(ωt)+cos(2ωt))2
(79)
Ht(z,t)={+1,u(z,t)>01,u(z,t)<0
(80)

Equation (77) is the wave load based on Morison’s formula rewritten from Eq. (13). Ht(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 F1(z, t) and Ht(z, t) F2(z, t) separately and adding them by using the superposition principle. Substituting F1(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 F1(z, t) can be proposed in the form of Eq. (65), where Tn(t) is found by solving Eq. (70) for F1(z, t).

However, the solution for Eq. (70) when the second term of Eq. (77), Ht(z, t) F2(z, t), is substituted as Q(z, t) in the right side of it will be challenging. This term consists of F2(z, t) in Eq. (79) in which temporal and spatial functions are squared and Ht(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
F2(z,t)=12ρwCDDMonf22cosh2(2kz)(f2(z)cos2(ωt)+2f(z)cos(ωt)cos(2ωt)+cos2(2ωt))
(81)
By using the following trigonometrical relationships
cos2(θ)=12(cos(2θ)+1)
(82)
cos(θ)cos(2θ)=12(cos(θ)+cos(3θ))
(83)
and substituting them into Eq. (81), it yields
F2(z,t)=14ρwCDDMonf22cosh2(2kz)×(f2(z)+1+2f(z)cos(ωt)+f2(z)cos(2ωt)+2f(z)cos(3ωt)+cos(4ωt))
(84)
which is a function in which the trigonometrical terms are linear.
The value of Ht(z, t) can be determined by the sign of u(z, t). Recalling u(z, t) from Eq. (75) and using the trigonometric relationship in Eq. (82) yields
u(z,t)=f2cosh(2kz)(2cos2(ωt)+f(z)cos(ωt)1)
(85)
The aforementioned equation reveals that the sign of 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:
t1(z)=1ωcos1(f(z)+f2(z)+84)
(86)
t2(z)=2πωt1(z)
(87)
Therefore,
Ht(z,t)={+1,0t<t1(z)ort2(z)<t2πω1,t1(z)<t<t2(z)
(88)

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

Therefore, Eq. (13) can be written in the following form:
Q(z,t)=i=12P1i(z)sin(iωt)+Ht(z,t)P2(z)+Ht(z,t)j=14P3j(z)cos(jωt)
(89)
where
P11(z)=ρwCMAMonωf2cosh(2kz)f(z)
(90)
P12(z)=2ρwCMAMonωf2cosh(2kz)
(91)
P2(z)=0.25ρwCDDMonf22cosh2(2kz)(f2(z)+1)
(92)
P31(z)=0.25ρwCDDMonf22cosh2(2kz)2f(z)
(93)
P32(z)=0.25ρwCDDMonf22cosh2(2kz)f2(z)
(94)
P33(z)=P31(z)
(95)
P34(z)=0.25ρwCDDMonf22cosh2(2kz)
(96)
The solution for Eq. (89) can be found by using the superposition principle since the properties of the system is linear. The solution for Eq. (89) can be expanded in the natural modes as follows:
x(z,t)=n=1Tn(t)Xn(z)
(97)
where Tn(t) satisfies
T¨n(t)+ωn2Tn(t)=i=12v1nisin(iωt)+v2n+j=14v3njcos(jωt)
(98)
v1ni=1EIMon0dP1i(z)Xn(z)dz
(99)
v2n=1EIMon0dHt(z,t)P2(z)Xn(z)dz
(100)
v3nj=1EIMon0dHt(z,t)P3j(z)Xn(z)dz
(101)
Solving Eq. (98) is possible by using the superposition principle in three parts, namely, T1n(t), T2n(t), and T3n(t), where
T¨1n(t)+ωn2T1n(t)=i=12v1nisin(iωt)
(102)
T¨2n(t)+ωn2T2n(t)=v2n
(103)
T¨3n(t)+ωn2T3n(t)=j=14v3njcos(jωt)
(104)
Therefore,
Tn(t)=T1n(t)+T2n(t)+T3n(t)
(105)
For Eq. (100), substituting into Eq. (103) results in
T¨2n(t)+ωn2T2n(t)=1EIMon0dHt(z,t)P2(z)Xn(z)dz
(106)

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

T2n(t)=1ωn1EIMon0t0dHt(z,τ)P2(z)Xn(z)sin(ωn(tτ))dzdτ
(107)
In the aforementioned equation, the function Ht(z, t) poses challenges in solving the integration analytically. To solve the double integration in the aforementioned equation, it requires the removal Ht(z, t) from inside the integrations to have a conventional double integration. As defined in Eq. (88), Ht(z, t) governs the sign of v2n. It can be either positive or negative, depending on which time frame it is evaluated. Ht(z, t) is positive in the time frame of [0, t1(z)], negative in the time frame of [t1(z), t2(z)], positive in the time frame of [t2(z), t2(z) + 2t1(z)], and so on. By introducing a new integer variable, m, by which Ht(z, t) is positive when m is an odd integer and an even number when it is negative. Therefore, Ht(z, t) can be redefined as follows:
Ht(z,t)={+1,m=1,3,5,1,m=2,4,6,
(108)
Since t ∈ [te(m−1)(z), tem(z)], one can conclude
Ht(z,t)=(1)m1
(109)
where te0(z) = 0, te1(z) = t1(z), te2(z) = t2(z), te3(z) = t2(z) + 2t1(z), and so on. Substituting the aforementioned equation into Eq. (107) and reversing the integration order yields
T2n(t)=1ωn1EIMon0dP2(z)Xn(z)×(0t1(z)sin(ωn(tτ))dτt1(z)t2(z)sin(ωn(tτ))dτ+t2(z)t2(z)+2t1(z)sin(ωn(tτ))dτ+(1)m1te(m1)(z)tsin(ωn(tτ))dτ)dz
(110)
Similarly, for Eq. (104), the solution is expressed as follows:
T3n(t)=j=141ωn1EIMon0dP3j(z)Xn(z)(0t1(z)cos(jωt)sin(ωn(tτ))dτt1(z)t2(z)cos(jωt)sin(ωn(tτ))dτ+t2(z)t2(z)+2t1(z)cos(jωt)sin(ωn(tτ))dτ+(1)m1te(m1)(z)tcos(jωt)sin(ωn(tτ))dτ)dz
(111)

Therefore, a method of solving the double integration in Eq. (107) is proposed by removing the function Ht(z, t), or better to say the staging of the integration domain, as represented in Eqs. (110) and (111).

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].

Table 2

DTU 10 MW OWT structural properties [74]

SymbolValue
Tower length (m)LTow119
Tower average diameter (m)DTow9.6
Tower average thickness (m)tTow0.0295
Tower Young's modulus (GPa)ETow210
Nacelle-rotor assembly mass (kg)Mn676,723
Nacelle-rotor assembly rotational inertia (kg · m2)Jp1.7 × 108
Platform level from mudline (m)LPlat45
Monopile average diameter (m)DMon8.3
Monopile average thickness (m)tMon0.09
Monopile Young's modulus (GPa)EMon210
Material density (kg/m3)ρs8500
SymbolValue
Tower length (m)LTow119
Tower average diameter (m)DTow9.6
Tower average thickness (m)tTow0.0295
Tower Young's modulus (GPa)ETow210
Nacelle-rotor assembly mass (kg)Mn676,723
Nacelle-rotor assembly rotational inertia (kg · m2)Jp1.7 × 108
Platform level from mudline (m)LPlat45
Monopile average diameter (m)DMon8.3
Monopile average thickness (m)tMon0.09
Monopile Young's modulus (GPa)EMon210
Material density (kg/m3)ρs8500

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

Table 3

Hydrodynamic loading properties

SymbolValue
Water depth (m)d35
Drag coefficientCD0.65
Added mass coefficientCA1
Inertia coefficientCM2
Sea water density (kg/m3)ρw1025
SymbolValue
Water depth (m)d35
Drag coefficientCD0.65
Added mass coefficientCA1
Inertia coefficientCM2
Sea water density (kg/m3)ρw1025

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.

Table 4

The values of the coupled springs [60]

SymbolValue
Lateral stiffness (GN/m)KL2.48
Cross stiffness (GN)KLR−20.7
Rotational stiffness (GN · m)KR412
SymbolValue
Lateral stiffness (GN/m)KL2.48
Cross stiffness (GN)KLR−20.7
Rotational stiffness (GN · m)KR412

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, fn, 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).

Fig. 2
The normalized natural mode shapes of the system with the added mass (MSL = mean sea level)
Fig. 2
The normalized natural mode shapes of the system with the added mass (MSL = mean sea level)
Close modal
Table 5

Solution variables

d/laTowaMAaMUα1α2α3α4α5α6γATγUT
0.22712,12514,58575015.84E-04−4.88E-030.0970.1881.59E-048.48E-070.9111.271
d/laTowaMAaMUα1α2α3α4α5α6γATγUT
0.22712,12514,58575015.84E-04−4.88E-030.0970.1881.59E-048.48E-070.9111.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.

Table 6

The natural mode frequencies of the system with and without the effect of added mass

Mode numberfn (Hz)Difference (%)
No added massAdded mass
10.1665610.166393−0.1
21.134631.0322−9.0
32.38881.98416−16.9
44.36863.8174−12.6
58.0256.593−17.8
612.1989.8905−18.9
Mode numberfn (Hz)Difference (%)
No added massAdded mass
10.1665610.166393−0.1
21.134631.0322−9.0
32.38881.98416−16.9
44.36863.8174−12.6
58.0256.593−17.8
612.1989.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.

Table 7

Comparison of the first natural frequency between the proposed method and full 3D FE-based model by Ref. [60]

Water depth (m)Platform level (m)Monopile outer diameter (m)Monopile thickness (cm)First natural frequency (Hz)Deviation (%)
Alkhoury et al. [60]Proposed model
25358.390.20090.1731−13.8
100.2020.1743−13.7
120.20560.1771−13.9
9100.20740.177−14.7
120.21120.1795−15.0
140.21580.1813−16.0
10110.21770.1809−16.9
130.21860.1828−16.4
150.22070.1837−16.8
35458.390.19090.1663−12.9
100.19280.1682−12.8
120.1970.171−13.2
9100.19920.1719−13.7
120.20380.1734−14.9
140.20920.176−15.9
10110.20960.1768−15.6
130.21290.1784−16.2
150.21550.1795−16.7
45558.390.18110.1592−12.1
100.18360.1614−12.1
120.1890.1648−12.8
9100.19090.1658−13.1
120.19620.1686−14.1
140.20230.1708−15.6
10110.2030.1717−15.4
130.2070.1736−16.1
150.21010.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
25358.390.20090.1731−13.8
100.2020.1743−13.7
120.20560.1771−13.9
9100.20740.177−14.7
120.21120.1795−15.0
140.21580.1813−16.0
10110.21770.1809−16.9
130.21860.1828−16.4
150.22070.1837−16.8
35458.390.19090.1663−12.9
100.19280.1682−12.8
120.1970.171−13.2
9100.19920.1719−13.7
120.20380.1734−14.9
140.20920.176−15.9
10110.20960.1768−15.6
130.21290.1784−16.2
150.21550.1795−16.7
45558.390.18110.1592−12.1
100.18360.1614−12.1
120.1890.1648−12.8
9100.19090.1658−13.1
120.19620.1686−14.1
140.20230.1708−15.6
10110.2030.1717−15.4
130.2070.1736−16.1
150.21010.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. 36. 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, KL, and cross-coupled stiffness, KLR, 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.

Fig. 3
Variation of system’s wavenumber, βTn, for different values of the water depth ratio, d/L, in the first five modes
Fig. 3
Variation of system’s wavenumber, βTn, for different values of the water depth ratio, d/L, in the first five modes
Close modal
Fig. 4
Variation of normalized wavenumber of the system for the above water section, βTn, for different values of α3 in the first and fifth modes
Fig. 4
Variation of normalized wavenumber of the system for the above water section, βTn, for different values of α3 in the first and fifth modes
Close modal
Fig. 5
Variation of normalized wavenumber of the system for the above water section, βTn, for different values of α5 in the first and fifth modes
Fig. 5
Variation of normalized wavenumber of the system for the above water section, βTn, for different values of α5 in the first and fifth modes
Close modal
Fig. 6
Variation of normalized wavenumber of the system for the above water section, βTn, for different values of α6 in the first and fifth modes
Fig. 6
Variation of normalized wavenumber of the system for the above water section, βTn, for different values of α6 in the first and fifth modes
Close modal
Table 8

The range of solution variables used for the parametric study and the equality to the parameters

Solution variableEquivalent to
α1 = ∞KL = ∞
α2 = 0KLR = 0
0.2 < α3 < ∞271 (GN · m) < KR < ∞
α4 = 1EIMon = EITow = 1357 GN · m2
0 < α5 < 10−30 < JP < 1.35 × 109 kg · m2
0 < α6 < 10−60 < Mn < 1357 ton
γAT = 1AMon = ATow
γUT = 1.36See Eq. (62)
Solution variableEquivalent to
α1 = ∞KL = ∞
α2 = 0KLR = 0
0.2 < α3 < ∞271 (GN · m) < KR < ∞
α4 = 1EIMon = EITow = 1357 GN · m2
0 < α5 < 10−30 < JP < 1.35 × 109 kg · m2
0 < α6 < 10−60 < Mn < 1357 ton
γAT = 1AMon = ATow
γUT = 1.36See 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.4L 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.2L 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.

Fig. 7
(a) The wave horizontal particle velocity and (b) wave load obtained from Morison’s equation, for H = 5.1 m and λ = 132 m
Fig. 7
(a) The wave horizontal particle velocity and (b) wave load obtained from Morison’s equation, for H = 5.1 m and λ = 132 m
Close modal
Table 9

Specification of the ocean wave applied to the reference tower

Height H (m)Ocean wavelength λ (m)Water depth d (m)Frequency f (Hz)PeriodHgT2dgT2HdUrsell number Ur
5.1132350.1059.530.005720.03920.1452.072
Height H (m)Ocean wavelength λ (m)Water depth d (m)Frequency f (Hz)PeriodHgT2dgT2HdUrsell number Ur
5.1132350.1059.530.005720.03920.1452.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 t1(z) and t2(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, Ht(z, t) is defined as a function t1(z) and t2(z) in Eq. (80). This poses difficulties in evaluating Eqs. (100) and (101). By defining Ht(z, t) as Eq. (109) and carrying out the integration, it becomes possible to remove Ht(z, t). In addition, changing the integration order by which the temporal integration is taken first in the domain as a function of t1(z) and t2(z) results in trigonometrical functions. By introducing t1(z) and t2(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.

Fig. 8
The variation of the wave horizontal particle velocity at around t1(z), for wave with H = 5.1 m and λ = 132 m
Fig. 8
The variation of the wave horizontal particle velocity at around t1(z), for wave with H = 5.1 m and λ = 132 m
Close modal

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

Fig. 9
The effect of the term Ht(z, t) in Morison’s formula with the second-order wave kinematics for H = 5.1 m and λ = 132 m at z = 35 m
Fig. 9
The effect of the term Ht(z, t) in Morison’s formula with the second-order wave kinematics for H = 5.1 m and λ = 132 m at z = 35 m
Close modal

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 Ht(z, t) in the response in Eqs. (110) and (111). The indirect presence of Ht(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.

Fig. 10
The response of the reference tower at different heights
Fig. 10
The response of the reference tower at different heights
Close modal

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.

Fig. 11
The response of the reference tower at the early stages of the motion
Fig. 11
The response of the reference tower at the early stages of the motion
Close modal

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 CA = 1 and CA = 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.

Fig. 12
Comparison between the response of the tower at the hub level with and without the added mass
Fig. 12
Comparison between the response of the tower at the hub level with and without the added mass
Close modal

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.

Fig. 13
Comparison between the response of the reference tower with and without the added mass at early stages of the motion
Fig. 13
Comparison between the response of the reference tower with and without the added mass at early stages of the motion
Close modal

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.

Table 10

The properties of an OWT for the numerical comparison

SymbolValue
Structural properties
Tower length (m)LTow70
Tower average diameter (m)DTow6
Tower average thickness (m)tTow0.045
Tower Young's modulus (GPa)ETow210
Nacelle-rotor assembly mass (kg)Mn0
Nacelle-rotor assembly rotational inertia (kg · m2)Jp0
Platform level from mudline (m)LPlat45
Monopile average diameter (m)DMon6
Monopile average thickness (m)tMon0.045
Monopile Young's modulus (GPa)EMon210
Material density (kg/m3)ρs7820
Support stiffness
Lateral stiffness (GN/m)KL
Cross stiffness (GN)KLR0
Rotational stiffness (GN · m)KR205.72
Hydrodynamic loading properties
Water depth (m)d30
Drag coefficientCD0.65
Added mass coefficientCA1
Inertia coefficientCM2
Sea water density (kg/m3)ρw1020
SymbolValue
Structural properties
Tower length (m)LTow70
Tower average diameter (m)DTow6
Tower average thickness (m)tTow0.045
Tower Young's modulus (GPa)ETow210
Nacelle-rotor assembly mass (kg)Mn0
Nacelle-rotor assembly rotational inertia (kg · m2)Jp0
Platform level from mudline (m)LPlat45
Monopile average diameter (m)DMon6
Monopile average thickness (m)tMon0.045
Monopile Young's modulus (GPa)EMon210
Material density (kg/m3)ρs7820
Support stiffness
Lateral stiffness (GN/m)KL
Cross stiffness (GN)KLR0
Rotational stiffness (GN · m)KR205.72
Hydrodynamic loading properties
Water depth (m)d30
Drag coefficientCD0.65
Added mass coefficientCA1
Inertia coefficientCM2
Sea water density (kg/m3)ρw1020

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.

Fig. 14
Comparison of the results of the proposed solution with the numerical one at z = 115 m
Fig. 14
Comparison of the results of the proposed solution with the numerical one at z = 115 m
Close modal

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:

  1. 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.

  2. 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.

  3. 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:

  1. The drag term of Morison’s formula, neglected by many researchers, is successfully included in finding the response of the system by defining Ht(z, t) to remove the absolute value function in the drag term.

  2. 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.

References

1.
Wind Europe
,
2019
, “
The European Offshore Wind Industry—Key Trends and Statistics 2019
,” Wind Europe, https://windeurope.org/intelligence-platform/product/wind-energy-in-europe-in-2019-trends-and-statistics/
2.
Adeli
,
H.
, and
Kim
,
H.
,
2009
,
Wavelet-Based Vibration Control of Smart Buildings and Bridges
,
CRC Press, Taylor & Francis
,
Boca Raton, FL
.
3.
Kim
,
H.
, and
Adeli
,
H.
,
2004
, “
Hybrid Feedback-Least Mean Square Algorithm for Structural Control
,”
J. Struct. Eng.
,
130
(
1
), pp.
120
127
.
4.
Adeli
,
H.
, and
Saleh
,
A.
,
1997
, “
Optimal Control of Adaptive/Smart Bridge Structures
,”
J. Struct. Eng.
,
123
(
2
), pp.
218
226
.
5.
Saleh
,
A.
, and
Adeli
,
H.
,
1994
, “
Parallel Algorithms for Integrated Structural/Control Optimization
,”
J. Aerosp. Eng.
,
7
(
3
), pp.
297
314
.
6.
Saleh
,
A.
, and
Adeli
,
H.
,
1996
, “
Parallel Eigenvalue Algorithms for Large-Scale Control-Optimization Problems
,”
J. Aerosp. Eng.
,
9
(
3
), pp.
70
79
.
7.
Saleh
,
A.
, and
Adeli
,
H.
,
1997
, “
Robust Parallel Algorithms for Solution of Riccati Equation
,”
J. Aerosp. Eng.
,
10
(
3
), pp.
126
133
.
8.
El-Khoury
,
O.
, and
Adeli
,
H.
,
2013
, “
Recent Advances on Vibration Control of Structures Under Dynamic Loading
,”
Arch. Comput. Meth. Eng.
,
20
(
4
), pp.
353
360
.
9.
Gutierrez Soto
,
M.
, and
Adeli
,
H.
,
2017
, “
Recent Advances in Control Algorithms for Smart Structures and Machines
,”
Expert Syst.
,
34
(
2
), p.
e12205
.
10.
Ghaedi
,
K.
,
Ibrahim
,
Z.
,
Adeli
,
H.
, and
Javanmardi
,
A.
,
2017
, “
Invited Review: Recent Developments in Vibration Control of Building and Bridge Structures
,”
J. Vibroeng.
,
19
(
5
), pp.
3564
3580
.
11.
Jiang
,
X.
, and
Adeli
,
H.
,
2008
, “
Neuro-Genetic Algorithm for Non-Linear Active Control of Structures
,”
Int. J. Numer. Methods Eng.
,
75
(
7
), pp.
770
786
.
12.
Jiang
,
X.
, and
Adeli
,
H.
,
2008
, “
Dynamic Fuzzy Wavelet Neuroemulator for Non-Linear Control of Irregular Building Structures
,”
Int. J. Numer. Methods Eng.
,
74
(
7
), pp.
1045
1066
.
13.
Li
,
Z.
, and
Adeli
,
H.
,
2016
, “
New Discrete-Time Robust H2/H∞ Algorithm for Vibration Control of Smart Structures Using Linear Matrix Inequalities
,”
Eng. Appl. Artif. Intell.
,
55
, pp.
47
57
.
14.
Gutierrez Soto
,
M.
, and
Adeli
,
H.
,
2017
, “
Many-Objective Control Optimization of High-Rise Building Structures Using Replicator Dynamics and Neural Dynamics Model
,”
Struct. Multidiscipl. Optim.
,
56
(
6
), pp.
1521
1537
.
15.
Soto
,
G.
, and
Adeli
,
H.
,
2017
, “
Multi-Agent Replicator Controller for Sustainable Vibration Control of Smart Structures Mariantonieta
,”
J. Vibroeng.
,
19
(
6
), pp.
4300
4322
.
16.
Gutierrez Soto
,
M.
, and
Adeli
,
H.
,
2018
, “
Vibration Control of Smart Base-Isolated Irregular Buildings Using Neural Dynamic Optimization Model and Replicator Dynamics
,”
Eng. Struct.
,
156
, pp.
322
336
.
17.
Li
,
Z.
, and
Adeli
,
H.
,
2018
, “
Control Methodologies for Vibration Control of Smart Civil and Mechanical Structures
,”
Expert Syst.
,
35
(
6
), p.
e12354
.
18.
Gutierrez Soto
,
M.
, and
Adeli
,
H.
,
2019
, “
Semi-Active Vibration Control of Smart Isolated Highway Bridge Structures Using Replicator Dynamics
,”
Eng. Struct.
,
186
, pp.
536
552
.
19.
Javadinasab Hormozabad
,
S.
, and
Soto
,
M. G.
,
2021
, “
Real-Time Damage Identification of Discrete Structures via Neural Networks Subjected to Dynamic Loading
,”
Proceedings Volume 11593, Health Monitoring of Structural and Biological Systems XV; 115932O
,
Virtual Online
.
20.
Azimi
,
M.
,
Eslamlou
,
A. D.
, and
Pekcan
,
G.
,
2020
, “
Data-Driven Structural Health Monitoring and Damage Detection Through Deep Learning: State-Ofthe- Art Review
,”
Sensors
,
20
(
10
), p.
2778
.
21.
Ngeljaratan
,
L.
,
Moustafa
,
M. A.
, and
Pekcan
,
G.
,
2021
, “
A Compressive Sensing Method for Processing and Improving Vision-Based Target-Tracking Signals for Structural Health Monitoring
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
36
(
9
), pp.
1203
1223
.
22.
Long
,
J.
, and
Büyüköztürk
,
O.
,
2020
, “
Collaborative Duty Cycling Strategies in Energy Harvesting Sensor Networks
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
35
(
6
), p.
534
548
.
23.
Sajedi
,
S.
, and
Liang
,
X.
,
2021
, “
Dual Bayesian Inference for Risk-Informed Vibration-Based Damage Diagnosis
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
36
(
9
), pp.
1168
1184
.
24.
Sajedi
,
S. O.
, and
Liang
,
X.
,
2021
, “
Uncertainty-Assisted Deep Vision Structural Health Monitoring
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
36
(
2
), pp.
126
142
.
25.
Huang
,
C. S.
,
Le
,
Q. T.
,
Su
,
W. C.
, and
Chen
,
C. H.
,
2020
, “
Wavelet-Based Approach of Time Series Model for Modal Identification of a Bridge With Incomplete Input
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
35
(
9
), pp.
947
964
.
26.
Gil-Gala
,
F. J.
,
Mencía
,
C.
,
Sierra
,
M. R.
, and
Varela
,
R.
,
2021
, “
Learning Ensembles of Priority Rules for Online Scheduling by Hybrid Evolutionary Algorithms
,”
Integr. Comput. -Aided Eng.
,
28
(
1
), pp.
65
80
.
27.
Sørensen
,
R. A.
,
Nielsen
,
M.
, and
Karstoft
,
H.
,
2020
, “
Routing in Congested Baggage Handling Systems Using Deep Reinforcement Learning
,”
Integr. Comput.-Aided Eng.
,
27
(
2
), pp.
139
152
.
28.
Ni
,
F. T.
,
Zhang
,
J.
, and
Noori
,
M. N.
,
2020
, “
Deep Learning for Data Anomaly Detection and Data Compression of a Long-Span Suspension Bridge
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
35
(
7
), pp.
685
700
.
29.
Ren
,
Q.
,
Li
,
M.
,
Li
,
H.
,
Song
,
L.
,
Si
,
W.
, and
Liu
,
H.
,
2021
, “
A Robust Prediction Model for Displacement of Concrete Dams Subjected to Irregular Water-Level Fluctuations
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
36
(
5
), pp.
577
601
.
30.
Ghofrani
,
F.
,
Pathak
,
A.
,
Mohammadi
,
R.
,
Aref
,
A.
, and
He
,
Q.
,
2020
, “
Predicting Rail Defect Frequency: An Integrated Approach Using Fatigue Modeling and Data Analytics
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
35
(
2
), pp.
101
115
.
31.
Tong
,
Z.
,
Yuan
,
D.
,
Gao
,
J.
, and
Wang
,
Z.
,
2020
, “
Pavement Defect Detection With Fully Convolutional Network and an Uncertainty Framework
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
35
(
8
), pp.
832
849
.
32.
Kalenjuk
,
S.
,
Lienhart
,
W.
, and
Rebhan
,
M. J.
,
2021
, “
Processing of Mobile Laser Scanning Data for Large-Scale Deformation Monitoring of Anchored Retaining Structures Along Highways
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
36
(
6
), pp.
678
694
.
33.
Zhu
,
M.
,
Zhu
,
H.
,
Guo
,
F.
,
Chen
,
X.
, and
Ju
,
J. W.
,
2021
, “
Tunnel Condition Assessment via Cloud Model-Based Random Forests and Self-Training Approach
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
36
(
2
), pp.
164
179
.
34.
Amezquita-Sanchez
,
J. P.
,
Park
,
H. S.
, and
Adeli
,
H.
,
2017
, “
A Novel Methodology for Modal Parameters Identification of Large Smart Structures Using MUSIC, Empirical Wavelet Transform, and Hilbert Transform
,”
Eng. Struct.
,
147
, pp.
148
159
.
35.
Pavlou
,
D.
,
2022
, “
A Deterministic Algorithm for Nonlinear, Fatigue-Based Structural Health Monitoring
,”
Comput.-Aided Civ. Infrastruct. Eng.
,
37
(
7
), pp.
809
831
.
36.
Bjørheim
,
F.
,
Siriwardane
,
S. C.
, and
Pavlou
,
D.
,
2022
, “
A Review of Fatigue Damage Detection and Measurement Techniques
,”
Int. J. Fatigue
,
154
, p.
106556
.
37.
Bjørheim
,
F.
,
Pavlou
,
D. G.
, and
Siriwardane
,
S. C.
,
2022
, “
Nonlinear Fatigue Life Prediction Model Based on the Theory of the S-N Fatigue Damage Envelope
,”
Fatigue Fract. Eng. Mater. Struct.
,
45
(
5
), pp.
1480
1493
.
38.
Luan
,
M.
,
Qu
,
P.
,
Jeng
,
D. S.
,
Guo
,
Y.
, and
Yang
,
Q.
,
2008
, “
Dynamic Response of a Porous Seabed–Pipeline Interaction Under Wave Loading: Soil–Pipeline Contact Effects and Inertial Effects
,”
Comput. Geotech.
,
35
(
2
), pp.
173
186
.
39.
Li
,
X.-J.
,
Gao
,
F.-P.
,
Yang
,
B.
, and
Zang
,
J.
,
2011
, “
Wave-Induced Pore Pressure Responses and Soil Liquefaction Around Pile Foundation
,”
Int. J. Offshore Polar Eng.
,
21
(
3
), pp.
233
239
.
40.
Chang
,
K. T.
, and
Jeng
,
D. S.
,
2014
, “
Numerical Study for Wave-Induced Seabed Response Around Offshore Wind Turbine Foundation in Donghai Offshore Wind Farm, Shanghai, China
,”
Ocean Eng.
,
85
, pp.
32
43
.
41.
Chen
,
L. F.
,
Zang
,
J.
,
Hillis
,
A. J.
,
Morgan
,
G. C. J.
, and
Plummer
,
A. R.
,
2014
, “
Numerical Investigation of Wave-Structure Interaction Using OpenFOAM
,”
Ocean Eng.
,
88
, pp.
91
109
.
42.
Sui
,
T.
,
Zhang
,
C.
,
Guo
,
Y.
,
Zheng
,
J.
,
Jeng
,
D.
,
Zhang
,
J.
, and
Zhang
,
W.
,
2015
, “
Three-Dimensional Numerical Model for Wave-Induced Seabed Response Around Mono-Pile
,”
hips Offshore Struct.
,
11
(
6
), pp.
667
678
.
43.
Zhang
,
C.
,
Zhang
,
Q.
,
Wu
,
Z.
,
Zhang
,
J.
,
Sui
,
T.
, and
Wen
,
Y.
,
2015
, “
Numerical Study on Effects of the Embedded Monopile Foundation on Local Wave-Induced Porous Seabed Response
,”
Math. Probl. Eng.
,
2015
, pp.
1
13
.
44.
Lin
,
Z.
,
Pokrajac
,
D.
,
Guo
,
Y.
,
Jeng
,
D.-s.
,
Tang
,
T.
,
Rey
,
N.
,
Zheng
,
J.
, and
Zhang
,
J.
,
2017
, “
Investigation of Nonlinear Wave-Induced Seabed Response Around Mono-Pile Foundation
,”
Coastal Eng.
,
121
, pp.
197
211
.
45.
Bazeos
,
N.
,
Hatzigeorgiou
,
G. D.
,
Hondros
,
I. D.
,
Karamaneas
,
H.
,
Karabalis
,
D. L.
, and
Beskos
,
D. E.
,
2002
, “
Static, Seismic and Stability Analyses of a Prototype Wind Turbine Steel Tower
,”
Eng. Struct.
,
24
(
8
), pp.
1015
1025
.
46.
Murtagh
,
P. J.
,
Basu
,
B.
, and
Broderick
,
B. M.
,
2005
, “
Along-Wind Response of a Wind Turbine Tower With Blade Coupling Subjected to Rotationally Sampled Wind Loading
,”
Eng. Struct.
,
27
(
8
), pp.
1209
1219
.
47.
Lavassas
,
I.
,
Nikolaidis
,
G.
,
Zervas
,
P.
,
Efthimiou
,
E.
,
Doudoumis
,
I. N.
, and
Baniotopoulos
,
C. C.
,
2003
, “
Analysis and Design of the Prototype of a Steel 1-MW Wind Turbine Tower
,”
Eng. Struct.
,
25
(
8
), pp.
1097
1106
.
48.
Alkhoury
,
P.
,
Soubra
,
A.-H.
,
Rey
,
V.
, and
Aït-Ahmed
,
M.
,
2022
, “
Dynamic Analysis of a Monopile-Supported Offshore Wind Turbine Considering the Soil-Foundation-Structure Interaction
,”
Soil Dyn. Earthquake Eng.
,
158
, p.
107281
.
49.
Graff
,
K. F.
,
2012
,
Wave Motion in Elastic Solids
,
Dover Publications
,
Mineola, New York.
50.
Meirovitch
,
L.
,
1967
,
Analytical Methods in Vibrations
,
Macmillan
,
New York
.
51.
Pavlou
,
D. G.
,
2021
, “
Soil–Structure–Wave Interaction of Gravity-Based Offshore Wind Turbines: An Analytical Model
,”
ASME J. Offshore Mech. Arct. Eng.
,
143
(
3
), p.
032101
.
52.
Wang
,
J.
,
Qin
,
D.
, and
Lim
,
T. C.
,
2010
, “
Dynamic Analysis of Horizontal Axis Wind Turbine by Thin-Walled Beam Theory
,”
J. Sound Vib.
,
329
(
17
), pp.
3565
3586
.
53.
Adeli
,
H.
, and
Karim
,
A.
,
1997
, “
Neural Network Model for Optimization of Cold-Formed Steel Beams
,”
J. Struct. Eng.
,
123
(
11
), pp.
1535
1543
.
54.
Karim
,
A.
, and
Adeli
,
H.
,
1999
, “
Global Optimum Design of Cold-Formed Steel Hat-Shape Beams
,”
Thin-Walled Struct.
,
35
(
4
), pp.
275
288
.
55.
Tashakori
,
A.
, and
Adeli
,
H.
,
2002
, “
Optimum Design of Cold-Formed Steel Space Structures Using Neural Dynamics Model
,”
J. Constr. Steel Res.
,
58
(
12
), pp.
1545
1566
.
56.
Arany
,
L.
,
Bhattacharya
,
S.
,
Adhikari
,
S.
,
Hogan
,
S. J.
, and
Macdonald
,
J. H. G.
,
2015
, “
An Analytical Model to Predict the Natural Frequency of Offshore Wind Turbines on Three-Spring Flexible Foundations Using Two Different Beam Models
,”
Soil Dyn. Earthquake Eng.
,
74
, pp.
40
45
.
57.
Arany
,
L.
,
Bhattacharya
,
S.
,
Macdonald
,
J. H. G.
, and
John Hogan
,
S.
,
2016
, “
Closed Form Solution of Eigen Frequency of Monopile Supported Offshore Wind Turbines in Deeper Waters Incorporating Stiffness of Substructure and SSI
,”
83
, pp.
18
32
.
58.
Amar Bouzid
,
D.
,
Bhattacharya
,
S.
, and
Otsmane
,
L.
,
2018
, “
Assessment of Natural Frequency of Installed Offshore Wind Turbines Using Nonlinear Finite Element Model Considering Soil-Monopile Interaction
,”
J. Rock Mech. Geotech. Eng.
,
10
(
2
), pp.
333
346
.
59.
Adeli
,
H.
,
Gere
,
J. M.
, and
Weaver
,
W.
,
1978
, “
Algorithms for Nonlinear Structural Dynamics
,”
ASCE J. Struct. Div.
,
104
(
2
), pp.
263
280
.
60.
Alkhoury
,
P.
,
Soubra
,
A. H.
,
Rey
,
V.
, and
Aït-Ahmed
,
M.
,
2021
, “
A Full Three-Dimensional Model for the Estimation of the Natural Frequencies of an Offshore Wind Turbine in Sand
,”
Wind Energy
,
24
(
7
), pp.
699
719
.
61.
Damgaard
,
M.
,
Ibsen
,
L. B.
,
Andersen
,
L. V.
, and
Andersen
,
J. K. F.
,
2013
, “
Cross-Wind Modal Properties of Offshore Wind Turbines Identified by Full Scale Testing
,”
J. Wind Eng. Ind. Aerodyn.
,
116
, pp.
94
108
.
62.
Prendergast
,
L. J.
,
Gavin
,
K.
, and
Doherty
,
P.
,
2015
, “
An Investigation Into the Effect of Scour on the Natural Frequency of an Offshore Wind Turbine
,”
Ocean Eng.
,
101
, pp.
1
11
.
63.
Prendergast
,
L. J.
,
Reale
,
C.
, and
Gavin
,
K.
,
2018
, “
Probabilistic Examination of the Change in Eigenfrequencies of an Offshore Wind Turbine Under Progressive Scour Incorporating Soil Spatial Variability
,”
Mar. Struct.
,
57
, pp.
87
104
.
64.
Dong
,
X.
,
Lian
,
J.
,
Wang
,
H.
,
Yu
,
T.
, and
Zhao
,
Y.
,
2018
, “
Structural Vibration Monitoring and Operational Modal Analysis of Offshore Wind Turbine Structure
,”
Ocean Eng.
,
150
, pp.
280
297
.
65.
Norén-Cosgriff
,
K.
, and
Kaynia
,
A. M.
,
2021
, “
Estimation of Natural Frequencies and Damping Using Dynamic Field Data From an Offshore Wind Turbine
,”
Mar. Struct.
,
76
, p.
102915
.
66.
Natarajan
,
A.
,
2014
, “
Influence of Second-Order Random Wave Kinematics on the Design Loads of Offshore Wind Turbine Support Structures
,”
Renewable Energy
,
68
, pp.
829
841
.
67.
Wang
,
Y.
,
2020
, “
Bottom Effects on the Tower Base Shear Forces and Bending Moments of a Shallow Water Offshore Wind Turbine
,”
Mar. Struct.
,
70
, p.
102705
.
68.
Wang
,
S.
,
Larsen
,
T. J.
, and
Bredmose
,
H.
,
2021
, “
Ultimate Load Analysis of a 10 MW Offshore Monopile Wind Turbine Incorporating Fully Nonlinear Irregular Wave Kinematics
,”
Mar. Struct.
,
76
, p.
102922
.
69.
Hirdaris
,
S. E.
,
Bai
,
W.
,
Dessi
,
D.
,
Ergin
,
A.
,
Gu
,
X.
,
Hermundstad
,
O. A.
,
Huijsmans
,
R.
,
Iijima
,
K.
,
Nielsen
,
U.D.
,
Parunov
,
J.
,
Fonseca
,
N.
,
Papanikolaou
,
A.
,
Argyriadis
,
K.
, and
Incecik
,
A.
,
2014
, “
Loads for Use in the Design of Ships and Offshore Structures
,”
Ocean Eng.
,
78
, pp.
131
174
.
70.
Bachynski
,
E.
,
Thys
,
M.
, and
Delhaye
,
V.
,
2019
, “
Dynamic Response of a Monopile Wind Turbine in Waves: Experimental Uncertainty Analysis for Validation of Numerical Tools
,”
Appl. Ocean Res.
,
89
, pp.
96
114
.
71.
Darvishi-Alamouti
,
S.
,
Bahaari
,
M. R.
, and
Moradi
,
M.
,
2017
, “
Natural Frequency of Offshore Wind Turbines on Rigid and Flexible Monopiles in Cohesionless Soils With Linear Stiffness Distribution
,”
Appl. Ocean Res.
,
68
, pp.
91
102
.
72.
Veritas
,
D. N.
,
2010
, “
Environmental Conditions and Environmental Loads,” DNV-RP-C205
.
73.
Morison
,
J. R.
,
Johnson
,
J. W.
, and
Schaaf
,
S. A.
,
1950
, “
The Force Exerted by Surface Waves on Piles
,”
J. Pet. Technol.
,
2
(
5
), pp.
149
154
.
74.
Bak
,
C.
,
Zahle
,
F.
,
Bitsche
,
R.
,
Kim
,
T.
,
Yde
,
A.
,
Henriksen
,
L. C.
,
Natarajan
,
A.
, and
Hansen
,
M.
,
2013
, "Description of theDTU 10 MW Reference Wind Turbine.", DTU Wind Energy Report-I-0092:1-138. https://dtu-10mw-rwt.vindenergi.dtu.dk.
75.
Zuo
,
H.
,
Bi
,
K.
, and
Hao
,
H.
,
2018
, “
Dynamic Analyses of Operating Offshore Wind Turbines Including Soil-Structure Interaction
,”
Eng. Struct.
,
157
, pp.
42
62
.
76.
Søren
,
P. H. S.
, and
Ibsen
,
L. B.
,
2013
, “
Assessment of Foundation Design for Offshore Monopiles Unprotected Against Scour
,”
Ocean Eng.
,
63
, pp.
17
25
.
77.
Wolfram Mathematica
,
2021
,
2021
, “
The World’s Definitive System for Modern Technical Computing
,” https://www.wolfram.com/mathematica