Effect of scour on the fatigue life of offshore wind turbines and its prevention through passive structural control

Offshore wind turbine (OWT) support structures are exposed to the risk of fatigue damage and scour, and this risk can be effectively mitigated by installing structural control devices such as tuned mass dampers (TMDs). However, time-varying scour altering OWTs' dynamic characteristics has an impact on the TMD design and fatigue life, which has rarely been studied before. In this paper, a simplified modal model is used to investigate the influence of scour and a TMD on the fatigue life evaluation of a 5 MW OWT's support structure, and a traditional method and a newly developed optimization technique are both presented to obtain TMD parameters. This optimization technique aims at finding optimal parameters of the TMD which maximize the fatigue life of a hotspot at the mudline, and the effect of time-varying scour can be considered. This study assumes that the TMD operates in the fore–aft (FA) direction, while the vibration in the side–side (SS) direction is uncontrolled. Results show that scour can decrease the fatigue life by about 24.1 % and that the TMD can effectively suppress vibration and increase the fatigue life. When the scour depth reaches 1.3 times the pile diameter, the TMD with a mass ratio of 1 % can increase the fatigue life of an OWT's support structure by about 64.6 %. Further, it is found that the fatigue life can be extended by 25 % with the TMD optimized by the proposed optimization technique rather than using a traditional design method which does not take the change in dynamic characteristics into account.

1 Introduction

With the continuous development of large-size fixed-bottom OWTs, local scour and scour protection of pile foundation have become a common issue (Wang et al., 2020; Wang et al., 2019; Zhang et al., 2022). Scour has a significant impact on the dynamic characteristics, vibration magnitudes, and thus fatigue life of OWTs under wind and wave loads. On the one hand, the action of currents and waves causes local scour pits around pile foundations, which reduces the burial depth of pile foundations. This phenomenon usually causes a reduction in the natural frequencies of OWTs and changes in other dynamic characteristics. This can potentially lead to resonance, large-amplitude stress cycles, and fatigue damage when one of the natural frequencies is close to the rotational frequency of the blades (Sørensen and Ibsen, 2013). On the other hand, current scour protection measures cannot completely avoid scour and have their own shortcomings. For example, armouring protection has the disadvantage that the projectile cannot be accurately cast in complex sea conditions and is easy to be washed away (Wang et al., 2023; Zhang et al., 2023). Flow-altering protection has the disadvantages of high cost and changing the dynamic characteristics of the foundation (Tang et al., 2023). As offshore structures, wind turbines are vulnerable to corrosion from seawater, which makes the fatigue problem worse (Amirafshari et al., 2021). Thus, the scour-induced changes in dynamic characteristics and the risk in resonance inevitably induce a further increase in fatigue damage and deserve in-depth research (Mayall et al., 2018).

Many researchers have studied the effect of scour on fatigue damage accumulation in OWTs. For instance, van der Tempel  (2006) investigated the frequency and fatigue of piles under different scour depths and concluded that scour has little effect on the natural frequencies but a major effect on fatigue damage. Zhang et al. (2021) found that scour depth has a significant influence on monopile impedance. Rezaei et al. (2018) showed that scour leads to an increase in the maximum bending moment of the monopile and a shortening of the fatigue life. Installing structural control devices is an effective way to mitigate the fatigue damage to OWTs. It was demonstrated that TMDs have a positive effect on reducing vibration amplitudes of wind turbine systems (Lackner and Rotea, 2011a; Dinh and Basu, 2015; Lu et al., 2023; Aydin et al., 2023). K. Dai et al. (2021) conducted a shaker experiment using a scaled wind turbine model and showed that the installed TMD can suppress the vibration of the structure more effectively considering soil–structure interaction (SSI).

In the previously mentioned studies, researchers have individually investigated the effect of scour on the structural vibration and fatigue and the structural control of TMDs for OWTs. However, in practice, the effect of scour combining structural control via TMDs could have a significant impact on OWTs' fatigue life. Moreover, whether considering scour could influence the design of TMDs and TMDs with different parameters can also have an impact on fatigue damage accumulation.

The purpose of this study is to explore the effect of scour on the fatigue life of wind turbine structures and the control effect of a TMD on the fatigue life of wind turbine structures under scour conditions. The authors present a case study with a 5 MW single-pile wind turbine to carry out related research. In this study, the software Abaqus is used to establish a detailed SSI model with different scour depths. A finite-element (FE) model which considers wind loads and TMD is established in MATLAB, and the scour effect is considered by establishing a relationship with the Abaqus model by means of the equivalent stiffness matrix. The finite-element model is also simplified to a modal model in order to achieve fast prediction of fatigue life. The TMD operates in the fore–aft (FA) direction and does not work in the side–side (SS) direction. This study investigates the effect of different scour depths on the performance of the TMD and the fatigue life of a 5 MW OWT's support structure, including a tower and a monopile foundation, and the optimization of the TMD's parameters considering time-varying scour depths to maximize fatigue life is also presented. This study provides some knowledge of the effects of the time-varying scour and the TMD on the fatigue life of wind turbines and a new TMD design method targeting the enhancement of fatigue resistance. The rest of the paper is organized as follows: Section 2 introduces the numerical models used in the research. Section 3 introduces the traditional TMD design method and the newly developed parameter optimization method. Section 4 describes the load cases for the fatigue analysis, the analysis results of this study, and the TMD parameter optimization results. Section 5 concludes the study.

2 Model description

2.1 Finite-element model and implementation of tuned mass damper

An FE model of a monopile-supported OWT installed with a TMD is established in MATLAB. This model contains a flexible tower, a rotor–nacelle assembly (RNA), and an external TMD, which consider the foundation flexibility. As shown in Fig. 1, the model is based on the widely used National Renewable Energy Laboratory (NREL) 5 MW reference OWT, and its detailed properties are shown in Table 1. Three-dimensional beam elements are used to create the FE model, and the theoretical basis is the standard Euler–Bernoulli beam theory. The wind turbine tower is divided into 18 beam elements, and the monopile between the mudline and the mean sea level (MSL) is divided into 4 beam elements. A convergence test performed by comparing the first natural frequencies shows that 22 beam elements are sufficient. Each element node has 6 degrees of freedom (DOFs) corresponding to translational and rotational motions in various directions. The mass matrix and stiffness matrix in the equation of motion of the OWT structure can be obtained using the material properties. The damping matrix is applied by means of Rayleigh damping, and the combined damping ratio of soil damping and structural damping is assumed to be 1 % (Chen and Duffour, 2018). The Rayleigh mass and stiffness coefficients, α₁ and α₂, are defined by ${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\frac{{\mathit{\zeta }}_C}{\frac{\mathrm{1}}{\mathrm{2}\mathit{\omega }}+\frac{\mathit{\omega }}{\mathrm{2}}}$. The variable ω is the natural frequency of the first FA mode, and ζ_C is the combined damping ratio. The RNA is represented by a lumped mass at the tower top.

Table 1 Basic properties of the NREL 5 MW reference OWT (Jonkman et al., 2009; Rezaei, 2017).

Figure 1 Schematic of NREL 5 MW wind turbine and scour effect.

The TMD is mounted on the top of the tower, and the effect of the TMD is considered by incorporating its mass, damping, and stiffness terms at relevant positions in the local mass, damping, and stiffness matrices of the beam element representing the top of the tower. The equation of motion of the OWT main structure is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathbf{M}}_{\mathrm{s}}{\ddot{\mathit{U}}}_{\mathrm{s}}& +{\mathbf{C}}_{\mathrm{s}}{\dot{\mathit{U}}}_{\mathrm{s}}+{\mathbf{K}}_{\mathrm{s}}{\mathit{U}}_{\mathrm{s}}+{\mathbf{C}}_{\mathrm{T}}({\dot{\mathit{U}}}_{\mathrm{s}}-{\dot{\mathit{U}}}_{\mathrm{T}})+{\mathbf{K}}_{\mathrm{T}}({\mathit{U}}_{\mathrm{s}}-{\mathit{U}}_{\mathrm{T}})\\ & ={\mathit{F}}_{\mathrm{wind}}+{\mathit{F}}_{\mathrm{wave}},\end{array}\end{array}$$

where M_s, C_s, and K_s are the mass, damping, and stiffness matrices of the main structure. C_T and K_T are matrices with the same dimensions containing c_T and k_T; ${\mathbf{C}}_{\mathrm{T}}=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& c_{\mathrm{T}}\end{array}\right];{\mathbf{K}}_{\mathrm{T}}=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& k_{\mathrm{T}}\end{array}\right]$. U_s is the displacement vector of the main structure; ${\mathit{U}}_{\mathrm{s}}=\left[\begin{array}{c}{\mathit{u}}_{\mathrm{s}-\mathrm{1}}\\ \mathrm{⋮}\\ {\mathit{u}}_{\mathrm{s}-top}\end{array}\right]$. U_T is the displacement vector containing u_T; ${\mathit{U}}_{\mathrm{T}}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{⋮}\\ u_{\mathrm{T}}\end{array}\right]$. F_wind and F_wave are the aerodynamic and wave load vectors. The equation of motion for the TMD can be represented by

$$\begin{array}{}\text{(2)}& m_{\mathrm{T}}{\ddot{u}}_T+c_{\mathrm{T}}({\dot{u}}_{\mathrm{T}}-{\dot{u}}_{\text{s-top}})+k_{\mathrm{T}}(u_{\mathrm{T}}-u_{\text{s-top}})=\mathrm{0},\end{array}$$

where m_T,c_T, and k_T are the mass, damping, and stiffness values of the TMD and u_T and u_s-top are the displacement of the TMD and the displacement of the top node. The modelling of the SSI is realized by an equivalent stiffness matrix, which will be introduced in detail subsequently in Sect. 2.3.

Wind loads were calculated using modified unsteady blade element momentum (BEM) theory (Branlard, 2017; Jonkman and Buhl, 2006) with Prandtl–Glauert corrections. Ignoring the iterative loop (Chen et al., 2021b) in the steady-state BEM code, the instantaneous aerodynamic forces were calculated for each time step within the time integration. The turbulent wind field was generated using the Kaimal spectrum according to the wind field parameters of International Electrotechnical Commission (2019), assuming moderate turbulence intensity. It should be noted that the aerodynamic loads from the rotor applied at the tower top were calculated using an aerodynamic force linearization technique previously developed by the authors (Chen et al., 2021a; Chen et al., 2020). This technique divides the aerodynamic loads into two parts. The first part is the quasi-steady aerodynamic force calculated by the BEM theory, which does not consider the influence of tower top motion. The second part considers the effect of aerodynamic damping by introducing an additional aerodynamic damping matrix. The adoption of this technique aims to facilitate the development of a simplified modal model for rapid fatigue calculation, which will be introduced in detail in Sect. 2.4. To represent the influence of the controller in the OWT, a standard relationship (Jonkman et al., 2009) between the mean wind speed, rotor rotation speed, and blade pitch angles, which represent the OWT's normal operational conditions, is adopted throughout the wind loading calculations.

Wave loads were calculated using the Morison equation, which includes viscous drag and inertial forces:

$$\begin{array}{}\text{(3)}& {\mathit{F}}_{\mathrm{wave}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{w}}{D}_{\mathrm{pile}}{C}_{\mathrm{d}}\left|{\dot{\mathit{u}}}_{\mathrm{w}}\right|{\dot{\mathit{u}}}_{\mathrm{w}}+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{\mathit{\rho }}_{\mathrm{w}}{D}_{\mathrm{pile}}^{\mathrm{2}}{C}_{\mathrm{m}}{\ddot{\mathit{u}}}_{\mathrm{w}},\end{array}$$

where ${\dot{\mathit{u}}}_{\mathrm{w}}$ and ${\ddot{\mathit{u}}}_{\mathrm{w}}$ are the velocity and acceleration vectors of water particles, C_d is the drag coefficient, D_pile is the diameter of the monopile between the mean sea level and the mudline, C_m is the inertia coefficient, and ρ_w is the density of water. C_d and C_m were chosen as 1 and 2 respectively as the recommended values in Shirzadeh et al. (2013). The wave profiles were obtained through the superposition of wave components, combining linear wave theory and the Joint North Sea Wave Project (JONSWAP) spectra (Hasselmann et al., 1973). The application of wind and wave loads is shown in Fig. 2.

Figure 2 Schematic of wind turbine load application.

2.2 Scour modelling in Abaqus

Using solid elements to model pile–soil interaction (S. Dai et al., 2021; Fard et al., 2022; Ma and Chen, 2021; Zdravković et al., 2015) is usually considered more accurate than the p–y-curve method (Liang et al., 2018; Song and Achmus, 2023) and the equivalent embedding method (Shahmohammadi and Shabakhty, 2020; Bergua et al., 2022). The solid-element method can also reduce the influence of an empirical formula on the results. Therefore, the solid-element method is used to establish the wind turbine scour model. The wind turbine scour model established in Abaqus contains soil, pile foundation, and tower, and the RNA is replaced by a concentrated mass located at the top of the tower. The diameter of the soil body is selected as 20 times the pile diameter, the soil under the pile foundation is selected as 2.5 times the pile diameter, and the total height of the soil body is 60 m. The soil body is made of homogeneous dense sandy soil, and the pile and tower are made of steel. The material parameters of the soil body, pile, and tower are shown in Table 2.

Table 2 Soil, pile, and tower material parameters.

The Mohr–Coulomb model is used for the soil, and the pile, tower, and nacelle are assumed to be elastic since they are much stiffer than the soil and do not deform plastically under normal operational conditions. The pile and tower are connected with a binding relationship. The normal contact between the pile and soil adopts the hard contact, and the tangential contact adopts the friction penalty function. The relative sliding friction factor at the interface, μ, is equal to tan (0.75 φ), where φ is the internal friction angle. The pile–soil contact is in the form of frictional contact, where mutual contact pairs are established between the pile and the soil, including the contacts between the pile bottom surface and the soil, the outside surface of the pile and the soil, and the inside surface of the pile and the soil core. The frictional contact between the pile bottom surface and soil is omitted due to the small area of the contact surface. These frictional contacts all adopt the face-to-face contact, and the contact discretization method adopts the face-to-face discretization method, considering the large stiffness of the main surface and the small stiffness of the slave surface. The perimeter of the soil body is translationally constrained, and the bottom surface of the soil adopts a fixed constraint. The eight-node linear brick element (C3D8R) is used to model the pile and soil, and the mesh division is realized by arranging seeds as shown in Fig. 3. The whole model is set up by adopting the modelling method of “element birth and death”, which makes the initial soil stress balance operation possible and sets up contacts and other related steps by killing and activating relevant elements.

Figure 3 Pile–soil interaction modelled in Abaqus.

The scour conditions can be represented by a deep conical pit around the pile under the long-term action of the waves and currents. According to the specification by Det Norske Veritas (DNV) (Det Norske Veritas AS, 2014b), the radius of the pit surface formed by scour, R, can be related to the depth of the scour pit by

$$\begin{array}{}\text{(4)}& R={\displaystyle \frac{D}{\mathrm{2}}}+{\displaystyle \frac{S}{\tan \mathit{\phi }}},\end{array}$$

where D is the diameter of the pile, S is the scour depth, and φ is the angle of the internal friction of the soil.

2.3 Equivalent stiffness matrix method

It is necessary to consider the effect of scour in the FE model in MATLAB. An equivalent stiffness matrix method is adopted in the FE model to consider the flexibility induced by the SSI. The 6 DOFs of nodes at the mudline are assumed to be constrained by a series of coupled springs, and the stiffness values of the coupled springs form a 6×6 stiffness matrix. For a specific stiffness term used in the FE model – for instance, the one relevant to the lateral displacement in the FA direction – the value of the stiffness term can be found from the relationship between the reaction force at the mudline and the pile top displacement (Jung et al., 2015). The equivalent stiffness schematic of the pile–soil interaction in the FA direction for the OWT is shown in Fig. 4.

Figure 4 Equivalent stiffness schematic of the pile–soil interaction in the FA direction.

According to the principle of virtual displacement and with the DOFs in other directions constrained, a unit displacement or rotation is first applied in one direction, and then the reaction force in that direction can be determined. The equivalent stiffness in that direction can be subsequently calculated by the relationship between the displacement and reaction force. Using the same approach, the stiffness terms corresponding to the remaining 5 DOFs are calculated. The stiffness terms in all 6 DOFs together form all the diagonal terms of the soil stiffness matrix. With the diagonal terms known, the off-diagonal stiffness terms can be found by applying a unit displacement in one direction and looking at the reaction force in the other concerned direction, with the other 4 DOFs constrained. Using the same principle, the off-diagonal terms can also be found from the relationship between the displacements and reaction forces, which ultimately results in a 6×6 stiffness matrix (Bergua et al., 2021; Pedersen and Askheim, 2021):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{F}}_{\mathrm{soil}}& =\left\{\begin{array}{c}\begin{array}{c}{F}_x\left(t\right)\\ F_y\left(t\right)\\ F_z\left(t\right)\end{array}\\ M_x\left(t\right)\\ M_y\left(t\right)\\ M_z\left(t\right)\end{array}\right\}\\ & =\left[\begin{array}{c}\begin{array}{lll}{k}_{xx}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{yy}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{zz}\end{array}\\ \begin{array}{ccc}\mathrm{0}& k_{\mathit{\theta }xy}& \mathrm{0}\\ k_{\mathit{\theta }yx}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\begin{array}{c}\begin{array}{lll}\mathrm{0}& k_{x\mathit{\theta }y}& \mathrm{0}\\ k_{y\mathit{\theta }x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ \begin{array}{lll}{k}_{\mathit{\theta }x\mathit{\theta }x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathit{\theta }y\mathit{\theta }y}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathit{\theta }z\mathit{\theta }z}\end{array}\end{array}\right]\\ & \left\{\begin{array}{c}\begin{array}{l}{u}_x\left(t\right)\\ u_y\left(t\right)\\ u_z\left(t\right)\end{array}\\ {\mathit{\theta }}_x\left(t\right)\\ {\mathit{\theta }}_y\left(t\right)\\ {\mathit{\theta }}_z\left(t\right)\end{array}\right\}={\mathbf{K}}_{\mathrm{soil}}{\mathit{u}}_{\mathrm{soil}},\end{array}\end{array}$$

where K_soil is the equivalent soil stiffness matrix, u_soil is the displacement vector, and F_soil is the reaction force vector. The equivalent soil stiffness matrix ignores the nonlinearity in the force–displacement relationship. This approach is suitable for fatigue analysis, as in normal operation conditions, the deformation of the soil around the monopile is relatively small and the nonlinearity in soil stiffness is very weak. The 6×6 soil stiffness matrix obtained from Abaqus is imported into the FE model in MATLAB. This modelling method combines the accuracy enhancement provided by the scour model in Abaqus using solid elements with the fast calculation speed and convenience of applying wind and wave loads using the FE model in MATLAB.

2.4 Rapid fatigue evaluation method

The established FE model in MATLAB can generate dynamic responses of the OWT, considering the wind and wave loads and the scour effect. However, a comprehensive fatigue life prediction in time domain needs to consider a large number of environmental states and load cases, so simulation efficiency is very important. Moreover, the TMD design optimization requires much more dynamic response time series. The FE model is not fast enough in this case. Therefore, a simplified modal model is developed from the FE model in MATLAB following the method developed in Chen et al. (2021b). The total aerodynamic forces from the rotor applied on the tower top node are linearized as the sum of a term corresponding to the forces for an assumed rigid tower and a term proportional to the tower top linear and angular velocities. The hydrodynamic forces are linearized by ignoring the relatively small monopile vibrations. The details of the force linearization can be found in a previous study by the authors (Chen et al., 2021a). Since the dynamic responses of the OWT are mainly dominated by the first two bending vibration modes, the FE model is reduced to a 4-DOF simplified modal model by considering only the first two bending modes in the FA and SS directions respectively. The development of the simplified 4-DOF modal model is briefly introduced as follows. Denoting the mass matrix and stiffness matrix of the OWT as M and K, including the TMD and the lumped soil stiffness matrix, the undamped vibration mode matrix, Ψ, can be obtained directly through eigen analysis. According to the relationship u=Ψα and multiplying the transpose of the undamped vibration matrix Ψ^T with the equation of motion, the following equation is obtained:

$$\begin{array}{}\text{(6)}& {\mathbf{\Psi }}^{\mathbf{T}}\mathbf{M}\mathbf{\Psi }\ddot{\mathit{\alpha }}+{\mathbf{\Psi }}^{\mathrm{T}}\mathbf{C}\mathbf{\Psi }\dot{\mathit{\alpha }}+{\mathbf{\Psi }}^{\mathrm{T}}\mathbf{K}\mathbf{\Psi }\mathit{\alpha }={\mathbf{\Psi }}^{\mathrm{T}}\mathit{F}.\end{array}$$

Then, the above equation can be rewritten as

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{‾}}{\mathbf{M}}\ddot{\mathit{\alpha }}+\stackrel{\mathrm{‾}}{\mathbf{C}}\dot{\mathit{\alpha }}+\stackrel{\mathrm{‾}}{\mathbf{K}}\mathit{\alpha }=\stackrel{\mathrm{‾}}{\mathit{F}},\end{array}$$

where α is the general coordinate vector, $\stackrel{\mathrm{‾}}{\mathbf{M}}$ is the modal mass matrix, $\stackrel{\mathrm{‾}}{\mathbf{C}}$ is the modal damping matrix, $\stackrel{\mathrm{‾}}{\mathbf{K}}$ is the modal stiffness matrix, and $\stackrel{\mathrm{‾}}{\mathit{F}}$ the modal load vector. Truncating Eq. (6) by only considering the first two bending modes, the FE model is reduced to a 4-DOF modal model, which can be used for a rapid fatigue analysis. The dynamic responses of the OWTs can be obtained by modal superposition after solving the general coordinate vector by time integration. In the 4-DOF simplified modal model, the cross section stress at any height can be calculated from the calculated node displacements. According to the dynamic stress extraction method provided by Pelayo et al. (2015), the cross section stress, σ_Z(t), at any moment at a given location can be obtained by

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_Z\left(t\right)=-E\left({\mathit{N}}^{e^{\prime \prime }}\right(z\left){\mathit{u}}_x^e\right(t)x+{\mathit{N}}^{e^{\prime \prime }}(z\left){\mathit{u}}_y^e\right(t\left)y\right),\end{array}$$

where u^e is the nodal displacement vector at the cross section, E is the material elastic modulus, and N^e is the elemental shape function vector of the FE model; subscripts x and y are the positions within the section at height z of the tower. After the cyclic counting of the stress time series using the rainfall counting method, the fatigue damage at the hotspot can be evaluated by utilizing the Palmgren–Miner rule based on the S–N fatigue calculation method. According to the DNV (Det Norske Veritas AS, 2014a) and considering the thickness effect, the S–N curve for steel under water can be obtained by the following equation:

$$\begin{array}{}\text{(9)}& \mathrm{logN}=\log \stackrel{\mathrm{‾}}{a}-m\cdot \log \left[\mathrm{\Delta }\mathit{\sigma }{\left({\displaystyle \frac{t}{t_{\mathrm{ref}}}}\right)}^k\right],\end{array}$$

where N is the number of cycles to failure and Δσ is the stress range. The variable Δσ is calculated from the nominal stress Δσ_nominal by the equation $\mathrm{\Delta }\mathit{\sigma }=\mathrm{SCF}\cdot \mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{nominal}}$, and SCF is the stress concentration factor. The variable m is the negative inverse slope of the S–N curve, $\log \stackrel{\mathrm{‾}}{a}$ is the intercept between the log N axis and the S–N curve, t_ref is the reference thickness for welded joints, and t is the thickness at which cracks may grow. Finally, t=t_ref is used for a thickness lower than t_ref. When t is larger than t_ref, t is the actual thickness of the pile. The variable k is the thickness exponent of the fatigue strength. For pile joints, the reference thickness, t_ref, is 25 mm. According to the DNV code, a bilinear S–N curve is usually used for offshore structures subjected mainly to typical wind and wave loads, using the class-E structural detail S–N curve shown in Table 3.

Table 3 Class-E structural detail S–N curves.

For stresses with variable amplitudes, the fatigue damage index is calculated using the Palmgren–Miner summation rule:

$$\begin{array}{}\text{(10)}& D_k=\sum _{i=\mathrm{1}}^{N_c}{\displaystyle \frac{n_i}{N_i}},\end{array}$$

where N_c is the total number of bins, n_i is the number of cycles in ith stress bin, N_i is the number of cycles to failure for the ith stress range, and D_k is the total fatigue damage index. The rainflow function in MATLAB is used for rainflow counting. When a stress time history is provided, this function can automatically calculate the ith stress range and the corresponding cycle number, n_i; N_c is the total number of stress ranges. Fatigue failure occurs at the hotspot when the fatigue damage index reaches 1.

3 Damper design and optimization method

Installing damping devices can efficiently reduce the vibration amplitudes of OWTs so that their service life can be greatly prolonged. Using TMDs as passive control devices is the most widely used method of controlling the vibration of OWT support structures. Usually, most TMDs are designed based on the dynamic characteristics of the OWTs determined in the preliminary design stage, without considering the changes in dynamic properties that may be caused by scour and soil degradation. In the real environment, scour can alter the dynamic characteristics of OWTs, potentially reducing the effectiveness of installed dampers or rendering them completely ineffective. Therefore, it is of great significance to consider the changes in dynamic properties caused by scour on the TMD design. The following two subsections first introduce the traditional TMD design method considering constant dynamic characteristics in the initial state. Then, an optimal parameter search method for the design of TMDs is presented, taking into account the effects of scour and the fatigue life evaluation.

Figure 5 Schematic diagram of TMD arrangement in the tower tube.

3.1 TMD design in initial state

As the dominant vibration mode of the OWT structure in operation is the first bending mode, the largest vibration amplitude occurs at the top of the tower, making the installation of the TMD at the tower top most effective. Therefore, the TMD is installed inside the steel tube at the tower top mainly to control the vibration in the FA direction, as shown in Fig. 5. Also, the TMD can be aligned with the FA direction by rotating the damper. Accordingly, the initial design of the TMD is mainly created based on the dynamic properties of the first-order mode. The initial design is used based on the assumption that the monopile foundation is not scoured.

Numerous studies have shown that a TMD can effectively suppress the vibration of a main structure when the mass ratio of the TMD to the main structure is 1 %–2 % (Lackner and Rotea, 2011b; Zhang et al., 2019). After determining the mass ratio, according to the classic TMD optimization theory proposed by Den Hartog (1957), the optimal frequency ratio of the TMD to the OWT structure is

$$\begin{array}{}\text{(11)}& {\mathit{\alpha }}_{\mathrm{opt}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathit{\mu }}}.\end{array}$$

The optimal damping ratio for the TMD can be calculated by

$$\begin{array}{}\text{(12)}& {\mathit{\xi }}_{\mathrm{opt}}=\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\mu }}{\mathrm{8}\left(\mathrm{1}+\mathit{\mu }\right)}}},\end{array}$$

where μ is the mass ratio of the TMD to the OWT structure, α_opt is the optimal frequency ratio of the TMD to the OWT structure, and ξ_opt is the optimal damping ratio of the TMD.

Considering that excessive mass would lead to increased construction costs, difficulties, and changes in the inherent characteristics of the original structure, the mass ratio of the TMD system to the main structure is initially set to 1 %. Moreover, previous studies have found that a tuned liquid column damper (TLCD) with a mass ratio of 1 % and a TMD with a mass ratio of 2 % can effectively suppress vibration (Colwell and Basu, 2009; Lackner and Rotea, 2011b; Zhang et al., 2019). According to Den Hartog's optimization theory for the initial TMD design, it can be determined that the optimal frequency ratio of the TMD to the main structure is 0.99 and that the optimal damping ratio of the TMD is 0.061. When the OWT support structure is not scoured, the first-order modal mass of the structure is 440 350 kg and the first-order modal frequency is 0.265 Hz. Therefore, according to the initial design parameters, the mass, stiffness coefficient, and damping coefficient of the TMD system are 4403.5 kg, 11 952 N m⁻¹, and 885 N s m⁻¹ respectively.

3.2 Fatigue-based damper optimization technique

After scouring occurs around the monopile foundation, the burial depth of the monopile and the natural frequencies of the OWT gradually change. The vibration mitigation effect of the TMD designed based on the dynamic parameters in the initial state can be reduced, which may lead to an increase in fatigue damage to the OWT support structure. Therefore, when designing the TMD, considering the influence of time-varying scour can enhance the performance of the TMD and result in a longer fatigue life of the support structure.

Figure 6 Flowchart of the TMD fatigue-life-based optimization technique.

Here, a fatigue-life-based optimization technique (FOT) to find optimal parameters of the TMD is developed in MATLAB, as shown in Fig. 6. In this technique, the frequency ratio, mass ratio, and damping ratio of the TMD are set as the optimal parameters to be searched, with the optimization objective being the fatigue life. When considering the time-varying scour process, the time-varying scour depth curve is first divided into a number of scour depths in increments of 0.1 D. For each scour depth, the fatigue damage is calculated and the total fatigue damage in a particular duration is then summarized. When the scour pit becomes deeper, the fatigue damage accumulates and finally reaches 1, which denotes the end of fatigue life. The simplified 4-DOF modal model incorporating scour modelling is used to generate the stress time series. The optimization problem is formulated to determine the optimal parameters of the TMD that maximize the fatigue life of the OWT support structure. The GlobalSearch function in MATLAB is used to solve the optimization problem. In the TMD optimization process, the mass ratio of the TMD is initially set to 1 %, and only the parameters of the frequency ratio and the damping ratio are optimized. Subsequently, in order to understand the optimization effect of the TMD when the value of the TMD mass ratio is not fixed, a mass ratio optimization interval is provided, making the mass ratio a variable within the optimization range.

4 Results

4.1 Environmental states and load cases

In this study, fatigue analyses are performed under 22 environmental states provided by van der Tempel (2006), taking into account both operational and parked conditions. These 22 environmental states are shown in Table 4. In operating conditions, the wind turbine withstands the aerodynamic load of the rotating rotor and the wind load of the tower. The wind load on the rotor is calculated using the BEM theory. In parked conditions, the wind turbine mainly withstands the aerodynamic load on the tower, and the aerodynamic damping is very small. The aerodynamic loading on the blades is calculated by directly looking at the aerodynamic loading coefficient table based on the local attack angles. The wind and wave loads are assumed to always act in the same direction to simplify the analysis. When the mean wind speeds are above the cut-in wind speed and below the cut-out wind speed, a 95 % wind turbine availability is assumed, following the setting in Velarde et al. (2020), meaning that the OWT does not produce power 5 % of the time when the mean wind speeds are in the operating range. For a specific combination of mean wind speed, wave period, and wave height, six different random seed numbers are used to produce various wind fields and wave profiles, minimizing the impact of randomness. To obtain the stress time histories at the mudline, a 700 s simulation is conducted for each random seed, and the response in the first 100 s is subtracted to eliminate the effect of initial transient vibration (Capaldo and Mella, 2023; Stieng and Muskulus, 2020).

Table 4 Environmental states, adopted from van der Tempel (2006).

In Table 4, V_w is the wind speed, T_z is the zero-crossing wave period, H_s is the wave height, and P_state is the probability of environmental state. To investigate the effect of scour and installation of the TMD on the fatigue damage accumulation, six load cases (LCs) are selected, as shown in Table 5. LC 1 is used as the reference case, and other cases are distinguished by different scour and TMD settings. For LC 4 to LC 6, the initial design of the TMD with the mass ratio of 1 % is used.

Table 5 Load case definition.

When considering the time-varying scour depth, for a particular time, t, the time-varying scour depth S can be predicted by the equation provided by Nakagawa and Suzuki (1976):

$$\begin{array}{}\text{(13)}& S={\left({\displaystyle \frac{t}{t_{\mathrm{1}}}}\right)}^{\mathrm{0.22}}D,\end{array}$$

where D is the diameter of the monopile and t₁ is the reference time and can be calculated by

$$\begin{array}{}\text{(14)}& t_{\mathrm{1}}=\mathrm{29.2}\cdot {\displaystyle \frac{D}{\sqrt{\mathrm{2}}\cdot u}}\cdot {\left({\displaystyle \frac{\sqrt{\mathrm{\Delta }\cdot g\cdot d_{\mathrm{50}}}}{\sqrt{\mathrm{2}}\cdot u-u_c}}\right)}^{\mathrm{3}}\cdot {\left({\displaystyle \frac{D}{d_{\mathrm{50}}}}\right)}^{\mathrm{1.9}}.\end{array}$$

The variable u is the tidal velocity and taken as 0.5 m s⁻¹, u_c is the critical shear velocity and taken as 0.37 m s⁻¹, g is the acceleration of gravity and taken as 9.8 m s⁻², and d₅₀ is the grain size of sea sand and taken as 0.2 mm. The parameter $\mathrm{\Delta }=\frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}}-\mathrm{1}$, where ρ_s is the density of sand and taken as 2.65 g cm⁻³ and ρ_w is the density of water and taken as 1 g cm⁻³. Rudolph et al. (2004) provided information on the sea state and measured the scour depth in the North Sea where the monopile N7 is situated. The measured scour depth was fitted well for the first 5 years based on the time-varying scour depth prediction equation shown in Eq. (13). As shown in Fig. 7, the data from the North Sea site can represent a typical ocean environment with time-varying scour and are used for the correlation analysis in this study.

Figure 7 Time-varying scour depth curve for pile N7 in the North Sea.

When conducting an analysis with the time-varying scour, an increment of scour depth equal to 0.1 D is utilized. At a specific scour depth, the fatigue damage is calculated, and the total fatigue damage during a longer period with varying scour depths can then be obtained by damage accumulation. According to the specification by DNV, the maximum depth of a local scour pit formed around a pile foundation is 1.3 times the diameter of the pile. Therefore, it is assumed that the local scour pit has a maximum scour depth of 1.3 D, at which the scour process achieves equilibrium.

4.2 Scour influence on natural frequencies

The scour of the soil around the monopile has an important effect on the natural frequencies of the OWT. For different scour depths, the first natural frequencies obtained by the models in Abaqus and MATLAB are compared in Fig. 8. The figure shows that the increase in the scour depth leads to a decrease in the first natural frequency of the OWT. The first natural frequency is 0.265 Hz when no scour occurs, and the natural frequency is reduced to 0.248 Hz when the depth of the scour pit reaches the maximum depth. The first natural frequency is reduced by 6.42 % due to the maximum scour depth. It shows that the natural frequency nearly monotonically decreases with the increase in the scour depth. The installation of a TMD also influences the natural frequency of the OWT main structure. The TMD with a mass ratio of 1 % makes the first natural frequency of the OWT main structure reduce to 0.251 Hz when no scour occurs, meaning that the natural frequency is reduced by 5.28 %.

Figure 8 Relationship between the wind turbine natural frequency and scour depth.

Figure 9 Displacement of the tower top and TMD in the 22nd environmental state.

In the TMD design process, the feasible displacement should be considered. The smaller the mass ratio of TMD, the larger the feasible displacement required. The 22nd environmental state corresponds to the greatest vibration responses of the wind turbine tower top due to large wind speed variations and lower aerodynamic damping, and the stroke of the TMD can be the largest. As shown in Fig. 9, the relative displacement between the TMD and the tower top is much smaller than the inner diameter of the wind turbine tower top in the 22nd environmental state. It shows that the stroke of the TMD is sufficient when the mass ratio of the TMD is 1 %.

Figure 10 Dynamic response of the wind turbine under wind–wave coupled loads for four operating conditions.

4.3 Dynamic response analysis

When the OWT in the operating state is under the ninth environmental state, which corresponds to the rated wind speed of 12 m s⁻¹, a comparison of the tower top displacements is made for LC 1, LC 3, LC 4, and LC 6, as shown in Fig. 10. These displacements are obtained from the FE model in MATLAB described in Sect. 2.1. By comparing the displacements from 300 to 420 s for LC 1 and LC 4, it can be found that the vibration amplitude of the tower top slightly decreases when the TMD is installed. Moreover, by comparing the displacement responses for LC 1 and LC 3, it can be found that the average of the displacement at the tower top increases when the scour depth reaches 1.3 D. This is because scour makes the OWT support structure become more flexible. It is known that the aerodynamic damping is large when the OWT is operating under the rated wind speed, so it is normal that the vibration mitigation effect of the TMD is less significant in this case. The effect of the TMD is more prominent in parked conditions with less aerodynamic damping. As shown in Fig. 11, the vibration mitigation effect of the TMD is more significant in the parked condition with a wind speed of 3 m s⁻¹.

Figure 11 The displacement response of the wind turbine tower in the parked condition with a wind speed of 3 m s⁻¹.

4.4 Fatigue calculation results

In Sect. 2.4, it is mentioned that in the process of fatigue life analysis, the 4-DOF simplified modal model is used to greatly decrease the calculation time. The accuracy test of the 4-DOF modal model in generating dynamic responses is first presented in this subsection. Under the turbulent wind field with a turbulence intensity of 11.9 % and an average wind speed of 12 m s⁻¹, the FE model and the 4-DOF simplified modal model are used to calculate the stress responses at the mudline for 10 min.

Figure 12 Comparison of stresses at the mudline from the FE model and the 4-DOF model in the time domain (a) and frequency domain (b).

As shown in Fig. 12, the stress responses from these two models are very close, confirming the accuracy of the 4-DOF modal model. The fatigue damage caused by the FE model in 10 min is $\mathrm{2.108}\times {\mathrm{10}}^{-\mathrm{7}}$, and the fatigue damage caused by the 4-DOF model in 10 min is $\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{7}}$, with an error of 0.05 %. Moreover, the calculation time of the 4-DOF simplified modal model is only about $\mathrm{1}/\mathrm{55}$th of that of the FE model, which shows that the 4-DOF simplified modal model is adequate to replace the FE model when conducting the fatigue life prediction.

The 4-DOF simplified modal model is used to conduct the fatigue life prediction for the OWT support structure for LC 1 to LC 6. A 10 min simulation is conducted for each of the six different random seed numbers to obtain stress time histories at the mudline. The location of the hotspot selected to assess fatigue damage is where the maximum stress occurs, specifically in the support structure cross section at the mudline. Although the location in the monopile where the moment reaches its maximum value can be below the mudline, the location at the mudline is chosen for simplicity. Further, since the SSI is modelled in the FE model by an equivalent soil stiffness matrix, it is not straightforward to determine the internal forces at the cross sections below the mudline. Given the stress time series at the selected hotspot, the corresponding fatigue damage is calculated. The fatigue damage for the combination of the mean wind speed, wave period, and wave height over 10 min is calculated by averaging the fatigue damage across six random seeds. For all 22 environment states, the 10 min fatigue damage is calculated. The fatigue life is predicted according to the Palmgren–Miner sum rule by combining these calculated fatigue damage values and the probabilities of the environmental states.

Figure 13 Fatigue life of the wind turbine with different scour depths.

For various scour depths, the fatigue life of the OWT considering both operating and parked conditions is predicted with or without TMD installation, and the results are shown in Fig. 13. It is shown that an increase in scour depth leads to a decrease in fatigue life and that an increasing fatigue life reduction rate can be observed when the scour depth increases. When no scour occurs and the TMD is not installed on the OWT, the OWT's fatigue life is 59.3 years, and the fatigue life drops to 45.0 years when considering the maximum scour depth of 1.3 D. There exist some uncertainties in the fatigue life prediction process due to the generation of a random wind field and wave profile. It should be noted that the predicted fatigue life is much longer than the normally adopted OWTs' design fatigue life of 25–30 years. This can be explained by the following reasons. First, the maximum moment of the OWT support structure is not at the cross section at the mudline where the selected hotspot is located. Second, the complex wind and wave directionality during the OWT's lifetime is simplified, which can influence the fatigue calculation results. Third, many other operational conditions such as start-up and shutdown phases are not considered in this study, which can also impact fatigue damage accumulation. Moreover, the installation of the TMD greatly extends the OWT support structure's fatigue life by 51.8 % when the pile foundation is not scoured.

Figure 14 Fatigue life of the wind turbine in six operating conditions.

Table 6 Optimization of TMD parameters.

The fatigue life prediction results of the OWT are obtained for all six LCs, as shown in Fig. 14. The fatigue life of the reference case LC 1, which is 59.3 years, is regarded as the reference fatigue life. It shows that the fatigue life decreases by 14.3 years, or about 24.1 %, when the scour depth is set to the maximum value of 1.3 D without applying the TMD compared to the reference fatigue life. When considering the time-varying scour, the fatigue life decreases by about 22.1 % from the reference value. When comparing the results for LC 1 and LC 4, it shows that the installation of the TMD results in a significant increase in the fatigue life of the OWT, with an increase in fatigue life of about 30.7 years, which is about 51.8 %. In LCs with the TMD installed, the fatigue life in LC 6 decreases by about 17.7 % when the scour depth reaches 1.3 D compared to the result for LC 4. But the fatigue life in LC 6 is still 1.25 times the reference fatigue life, which indicates that the imposition of TMD can effectively increase the fatigue life of the OWT by reducing vibration amplitudes.

4.5 Fatigue calculation with optimized TMD

To compare the optimization effect and speed up the optimal parameter search process, the mass ratio of TMD is first kept at 1 %. Before the optimization, the parameter ranges of the frequency ratio and the damping ratio need to be defined. The optimal frequency of the TMD is usually close to the frequency of the main structure, so the range of the frequency ratio is chosen to be from 0.8 to 1.1 for optimization. As the optimal damping ratio can vary in a relatively larger range, the range of the damping ratio for optimization is chosen to be from 1 % to 30 %. The optimization of the TMD is also conducted, with the mass ratio not fixed. A range of the mass ratio from 0.001 to 0.1 is used to optimize the TMD so that the influence of the mass ratio can be evaluated.

The optimal parameters obtained by FOT and the predicted fatigue life are listed in Table 6. The fatigue life for LC 5 and the parameters of the initially designed TMD are also shown in Table 6 for comparison. It shows that when the mass ratio is fixed at 1 %, the optimal frequency ratio is 0.94, the optimal damping ratio is 5 %, and the final fatigue life is 93.2 years. Compared to the fatigue life with the initially optimized TMD using the traditional method without considering scour, the fatigue life increases by 18.6 years, which is approximately 25 %. It indicates that the parameter search in the optimization process is correct and that it is better to use the TMD parameter search method to design the TMD after obtaining the time-varying scour curve. When the mass ratio range is taken from 0.1 % to 10 %, the optimal mass ratio of the TMD is 9.7 %, the frequency ratio is 0.92, the damping ratio is 15 %, and the final fatigue life is 133.2 years. In this case, the fatigue life of the OWT is significantly longer primarily because of the large mass ratio. However, in practice, it may be uneconomical to implement a TMD with such a large mass ratio.

5 Conclusions

This study establishes a rapid numerical model that can consider the effect of scour and the installation of a TMD, and the TMD operates only in the FA direction. The model is simplified by concentrated mass instead of RNA and ignores the nonlinearity of the equivalent stiffness matrix. The established model is used to investigate the influence of scour and the installed passive structural control device on the OWT's natural frequencies and fatigue life in 22 environmental states. An optimization technique was also developed to find the optimal parameters of the TMD, considering time-varying scour. Moreover, it shows that the vibration amplitude of the OWT can be effectively reduced by the TMD. On the one hand, the results show that the TMD reduces the vibration amplitude of the tower top. On the other hand, when the scour depth reaches 1.3 D, the wind turbine support structure becomes more flexible, with the displacement of the tower top increasing without TMD.

In addition, the fatigue calculation results show that the installation of the TMD significantly extends the fatigue life of the OWT, but scour can cause a reduced performance of the TMD. It is found that when scour and the scour-induced natural frequency reduction are considered during the OWT's lifetime, the performance of the initially designed TMD is not as good as the TMD optimized by the developed FOT in terms of fatigue life enhancement. Given a mass ratio of 1 %, the fatigue life can be extended by 25 % with the TMD optimized by FOT. This is because the FOT can consider the effect of time-varying scour. This study focuses solely on analysing scour, but it is important to note that other factors, like soil degradation, can also impact the dynamic characteristics of OWTs. These factors may influence on structural control devices and the evaluation of fatigue life. Additionally, during OWTs' lifetime, the properties of installed TMDs can also change, making the evaluation of TMDs' performance and OWTs' fatigue life more complicated. These factors are worth investigating in the future.