Dynamic analysis of wind turbines including nacelle–tower–foundation interaction for condition of incomplete structural parameters

Abstract

Considering the complexity structure of wind turbines, it is difficult to establish an accurate model for wind turbines. This article proposed an improved model based on the vibration signal analysis to investigate the onshore wind turbine dynamic behavior. The model is formulated using the Euler–Lagrangian approach in which the dynamic interaction effects between the nacelle and the tower and the effects between the tower and the foundation are considered. The nacelle is assumed as a mass with 2 degrees of freedom at the top of a flexible tower. And the tower is considered as a beam with elastic end support. The foundation stiffness can influence the dynamic response of wind turbines due to the soft soils surrounding the tower base. A rotational spring and a lateral spring derived from the foundation are used to reflect the interaction between the tower and the foundation. Finally, a comparison with existing National Renewable Energy Laboratory 5-MW baseline wind turbine model is performed. And a field experiment has been implemented to validate the model. The vibration signal of the nacelle and the tower is analyzed to determine the model parameters. Results show that the proposed model is accurate and practical for dynamic property analysis of wind turbine.

Introduction

Wind energy is one of the fastest growing renewable energies, according to the collected statistical data by the Global Wind Energy Council (GWEC), the global total installations wind power capacity was 369.6 GW at the end of 2014. ¹ Along with it, the reliability issues have been emerging. Vibration characteristics of wind turbines may essentially influence the stability, performance, and lifetime of wind turbine systems. There are many studies that have been devoted to improve the availability and survivability of wind turbines.^2–4 Recently, condition monitoring and fault detection systems have been proved an effective way to ensure that the wind turbine works more efficient, more robust, and less costly. ⁵ However, the dynamic characteristic of wind turbine is also an important issue in wind turbine design, operation, and maintenance but is rarely considered for a condition monitoring system or a fault diagnosis system.

Much attention has been paid to the dynamic analysis of wind turbine.^6–8 Dynamic analysis for wind turbines includes experimental modal analysis and theoretical modal analysis. The theoretical modal analysis can be applied more easily and widely. The main challenge in analyzing the dynamic properties of wind turbines with the theoretical model is estimating the parameters of the model. For experimental analysis, it considers the dynamic properties of a kinetic system by measuring and analyzing the input and response signal. ⁹ But there are limits to implement experimental test for large-scale wind turbines because of the inaccessibility and complicated construction of wind turbine. As vibration analysis is widely used in wind turbine condition monitoring and fault diagnosis system, the vibration signals of some components of the wind turbine can be obtained by the condition monitoring system and developed to improve the dynamic analysis of the wind turbine.

It is difficult to model the entire wind turbine in detail, but the structure of the wind turbine can be simplified as several key subsystems. Many subsystems may contribute to the vibration of the wind turbine, such as the flapping of blades, the rotation of the shaft and driven chain, and the oscillation of the tower. A comprehensive aeroelastic powerful simulator capable of predicting both the extreme and fatigue loads of two- and three-bladed horizontal-axis wind turbines (HAWTs), named fatigue, aerodynamics, structures, and turbulence (FAST), was developed by National Renewable Energy Laboratory (NREL).^10,11 The modules of FAST (AeroDyn, HydroDyn, etc.) rely on different physical domains of the coupled aero-hydro-servo-elastic solution. Therefore, FAST requires detailed parameters of the wind turbine, including aerodynamic shape, material properties, control strategies, and so on. However, these parameters are not always available, which will limit the use of FAST.

The blade characteristics in the low-order vibration conditions were investigated by Ye et al. ⁹ Murtagh et al.^12,13 investigated the dynamic interaction between the rotor blades and the tower and developed an analytical model for passive control of wind turbine vibrations. A hybrid multibody system method was proposed for wind turbine structure analysis, which could capture all the rigid and elastic mechanical components, couplings, and important motions based on nonlinear kinetics. ¹⁴ Kessentini et al. ¹⁵ used a mathematical model with flexible tower and blades to estimate the effects of the vibrations of the tower and blades. They found that the effects of the pitch angle and the blade orientation on the natural frequencies and the mode shapes of the wind turbine can be considered insignificant factors, while small vibrations of the tower may induce severe blade deflections. Staino and Basu ¹⁶ also developed a multi-modal mathematical model describing the dynamics of flexible rotor blades and their interaction with the turbine tower, taking the variable rotor speed into account. However, the above works ignored the foundation factor due to different soil types. Nevertheless, the soil–structure interaction (SSI) is also essential for the dynamic properties of wind turbine. The effect of foundation has been considered for the free vibration analysis of wind turbines in Adhikari and Bhattacharya ¹⁷ and Harte et al. ¹⁸ For the purpose of wind turbine condition monitoring, Liu ¹⁹ established a blade–cabin–tower coupling system to analyze the vibration and force in wind turbines. Recently, Liu performed a numerical vibration analysis by further considering the wind turbine system as a blade–cabin–tower coupling system. However, the model parameters are difficult to discern, and the essential property of the wind turbine such as the natural frequency is not included in the mentioned papers. For this reason, appropriate modeling of the relevant dynamics is necessary for reliable wind turbine vibration analysis applications.

In this article, a dynamic model of the vibrations for a three-bladed HAWT is formulated according to the Euler–Lagrangian ²⁰ approach. The model contains the effect of foundation and nacelle system and the effect of interaction of foundation and tower. Moreover, blade flapping influence is considered as a part of the external force effect, and the mass of the blades is lumped with the mass of the nacelle. At first, a mathematical model for a HAWT structure is characterized by a set of ordinary differential equations and four partial differential equations with the boundary conditions. Then, the vibration analysis and an experiment are performed to validate the proposed model, and the effects of a variety of model parameters are discussed. The result shows that the proposed model is able to represent the dynamic properties of the wind turbine. The approach may be used in preliminary quantitative design and building a condition monitoring system with the dynamic behaviors of the wind turbine.

Modeling the nacelle–tower–foundation system

A typical HAWT is shown in Figure 1, which includes several components such as blades, gearbox, generator, nacelle, tower, and foundation. We develop a mathematical model for an HAWT consisted of a flexible tower, a mechanical idealization flexible foundation, a rigid blade–rotor–gearbox–generator (BRGG) system with a translational nacelle, as shown in Figure 2. The multi-degree of freedom (DOF) superstructure model is formulated using the Euler–Lagrangian equations, ²¹ as follows

$$\frac{d}{\mathrm{dt}}\frac{\partial \mathcal{T}}{\partial {\stackrel{·}{q}}_i}-\frac{\partial \mathcal{T}}{\partial q_i}+\frac{\partial \mathcal{V}}{\partial q_i}=F_i$$

(1)

where $\mathcal{T}$ is the total kinetic energy of the system and $\mathcal{V}$ is the total potential energy of the system; $q_i$ is the generalized coordinate displacement; ${\stackrel{·}{q}}_i$ is the velocity at generalized coordinate $q_i$ ; and $F_i$ is the generalized force, the portion of the total forces not related to the potential energy of the system, such as gravity and spring forces appear in the potential energy expression and are not included here.

Figure 1.

The main structure of a horizontal-axis wind turbine.

Figure 2.

Wind turbine model.

The model is developed for the vibration analysis of the system on the top of tower. Here, we consider the effect of the vibration of blades as a part of the outside force $P(t)$ . As a result, a three-mass system is derived from the wind turbine model, where the three BRGGs were integrated as a rigid mass denoted as $M_0$ . The rigid mass has 2 DOFs translating along the horizon and vertical directions. The generalized stiffnesses are $k_h$ and $k_v$ , while $c_h$ and $c_v$ are the structural damping associated with the tower and BRGG in the horizontal and vertical directions, respectively. Furthermore, the model consisted of a flexible tower of uniform hollow cross section and a flexible foundation. The mass of the earth far outweighs the wind turbine; therefore, the motion of the foundation could be neglected.

The lumped mass of BRGG is modeled as a 2-DOF system which could be represented by $q_h$ in the horizontal direction and $q_v$ in the vertical direction, respectively. The tower of wind turbine is considered as an idealized Euler–Bernoulli beam. The variable bending stiffness of the tower is $\mathrm{EI}(x)$ , where x is the spatial coordinate described in Figure 2. The interaction between the tower and the foundation is modeled using two springs, the lateral spring with spring stiffness $k_l$ , and the rotational spring with spring stiffness $k_r$ . The tower has a top mass with rotary inertia J which is considered the idealized to the BRGG. The mass per unit length of the tower is $m(x)$ , and H is the height of the tower. The transverse deflection per unit height of the tower is $y(x,t)$ , and t is the time. The tower is modeled as a flexible beam with N modes of vibration in the horizontal direction, and the mode is described by the generalized variable $q_i(t)$ , $i=1,\dots ,N$ . And the mode of vibration is associated with the corresponding mode shape ${\mathrm{\Phi }}_i(x)$ . The total displacement along the tower is given by

$$y(x,t)=\sum _{i=1}^N{\mathrm{\Phi }}_i\left(\mathrm{x}\right)q_i\left(t\right)$$

(2)

The wind turbine vibration model is, therefore, described by $N+2$ generalized coordinates, which could be defined as

$$\mathbf{q}\left(t\right)=\left[\begin{array}{c}{q}_1\left(t\right)\hfill \\ \multicolumn{1}{c}{q_2\left(t\right)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{q_N\left(t\right)}\\ \multicolumn{1}{c}{q_h}\\ \multicolumn{1}{c}{q_v}\end{array}\right]\in {\mathbb{R}}^{N+2}$$

(3)

The total kinetic energy of the model can be written as

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{T}=\frac{1}{2}{M}_0({\stackrel{·}{q}}_h^2\left(t\right)+{\stackrel{·}{q}}_v^2\left(t\right))+\frac{1}{2}\underset{0}{\overset{H}{\int }}m\left(x\right){\stackrel{·}{y}}^2(x,t)\mathrm{dx}}\\ \multicolumn{1}{c}{=\frac{1}{2}{M}_0({\stackrel{·}{q}}_h^2\left(t\right)+{\stackrel{·}{q}}_v^2\left(t\right))+\frac{1}{2}\sum _{i=1}^N{M}_i{\stackrel{·}{q}}_i^2\left(t\right)}\end{array}$$

(4)

$\mathrm{where}(\stackrel{•}{•})$ denotes derivative with respect to time and $M_i$ is the equivalent mass corresponding to the generalized coordinates.

The total potential energy of the system is given by

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{V}=\frac{1}{2}{k}_h{(q_h\left(t\right)-y(H,t))}^2+\frac{1}{2}{k}_v{q}_v^2\left(t\right)}\\ +\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N\left[(k_{e,\mathrm{ij}}+k_{p,\mathrm{ij}}+k_{g,\mathrm{ij}})q_i\left(t\right)q_j\left(t\right)\right]\hfill \\ +\frac{1}{2}{k}_l{y}^2(0,t)+\frac{1}{2}{k}_{r}y{\prime }^2(0,t)\hfill \end{array}$$

(5)

where the terms $k_{e,\mathrm{ij}}$ , $k_{p,\mathrm{ij}}$ , and $k_{g,\mathrm{ij}}$ represent contributions to the elastic stiffness, geometric stiffness, and gravitational stiffness, respectively ²²

$$k_{e,\mathrm{ij}}=\underset{0}{\overset{H}{\int }}\mathrm{EI}\left(x\right)\left[{\mathrm{\Phi }}_i^{\prime \prime }\left(x\right){\mathrm{\Phi }}_j^{\prime \prime }\left(x\right)\right]\mathrm{dx}$$

(6)

$$k_{p,\mathrm{ij}}=-\underset{0}{\overset{H}{\int }}{M}_{0}g\left[{\mathrm{\Phi }}_i^{\prime }\left(x\right){\mathrm{\Phi }}_j^{\prime }\left(x\right)\right]\mathrm{dx}$$

(7)

$$k_{g,\mathrm{ij}}=-\underset{0}{\overset{H}{\int }}g\left[\underset{x}{\overset{H}{\int }}m\left(\epsilon \right)d\epsilon \right]\left[{\mathrm{\Phi }}_i^{\prime }\left(x\right){\mathrm{\Phi }}_j^{\prime }\left(x\right)\right]\mathrm{dx}$$

(8)

where E is Young’s modulus of elasticity of the tower material, $I(x)$ is the second moment of area of the tower about the axis, and ${\mathrm{\Phi }}_i^{\prime }(x)$ and ${\mathrm{\Phi }}_i^{\prime \prime }(x)$ denote the first and second derivatives of the mode shape with respect to x.

The dissipation function of the system is given by

$$\mathcal{G}=\frac{1}{2}\sum _{i=1}^N{c}_i{{\stackrel{·}{q}}^2}_i\left(t\right)+\frac{1}{2}{c}_h{{\stackrel{·}{q}}^2}_h\left(t\right)+\frac{1}{2}{c}_v{{\stackrel{·}{q}}^2}_v\left(t\right)$$

(9)

where $c_i,i=1,2,\dots ,N$ , denotes the structural damping associated with the mode of the tower.

According to equation (1), the equations of motion for the wind turbine vibration model with N + 2 DOFs can be written as

$$\mathbf{M}\left(t\right)\stackrel{··}{\mathbf{q}}\left(t\right)+\mathbf{C}\left(t\right)\stackrel{·}{\mathbf{q}}\left(t\right)+\mathbf{K}\left(t\right)\mathbf{q}\left(t\right)=\mathbf{F}\left(t\right)$$

(10)

where $\mathbf{F}(t)$ is generalized outside force, including the aerodynamic load on the blades and the nacelle $P(t)$ .

The mass matrix $\mathbf{M}(t)$ of the wind turbine model is defined as

$$\mathbf{M}\left(t\right)=\left[\begin{array}{ccccc}{M}_1\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & M_N\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & M_0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & 0\hfill & M_0\hfill \end{array}\right]\in {\mathbb{R}}^{(N+2)\times (N+2)}$$

(11)

where $M_i$ is the equivalent mass, same as in equation (4), can be defined as

$$M_i=\underset{0}{\overset{H}{\int }}m\left(x\right){\mathrm{\Phi }}_i^2\left(x\right)\mathrm{dx},i=1,2,\dots ,N$$

(12)

The damping matrix $\mathbf{C}(t)$ could be written as

$$\mathbf{C}\left(t\right)=\left[\begin{array}{ccccc}{c}_1\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & c_N\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & c_h\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & 0\hfill & c_v\hfill \end{array}\right]\in {\mathbb{R}}^{(N+2)\times (N+2)}$$

(13)

The stiffness matrix $\mathbf{K}(t)$ is given by

$$\begin{array}{c}\mathbf{K}\left(t\right)=\left[\begin{array}{ccc}{k}_h\tilde{\mathbf{\Phi }}\left(H\right)+k_l\tilde{\mathbf{\Phi }}\left(0\right)+k_r\tilde{\mathbf{\Phi }}\prime \left(0\right)+\hat{\mathbf{K}}\hfill & -k_h\mathbf{\Phi }\left(H\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{-k_h{\mathbf{\Phi }}^T\left(H\right)}& k_h\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & k_v\hfill \end{array}\right]\hfill \\ \in {\mathbb{R}}^{(N+2)\times (N+2)}\hfill \end{array}$$

(14)

where

$$\hat{\mathbf{K}}=\left[\begin{array}{cccc}{k}_{e,11}+k_{p,11}+k_{g,11}\hfill & k_{p,12}+k_{g,12}\hfill & \cdots \hfill & k_{p,1N}+k_{g,1N}\hfill \\ \multicolumn{1}{c}{k_{p,21}+k_{g,21}}& k_{e,22}+k_{p,22}+k_{g,22}\hfill & \cdots \hfill & k_{p,2N}+k_{g,2N}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{k_{p,N1}+k_{g,N1}}& k_{p,N2}+k_{g,N2}\hfill & \cdots \hfill & k_{e,\mathrm{NN}}+k_{p,\mathrm{NN}}+k_{g,\mathrm{NN}}\hfill \end{array}\right]\in {\mathbb{R}}^{N\times N}$$

(15)

assuming orthogonal modes of vibration implies $k_{e,\mathrm{ij}}=0$ , $\forall i\ne j$ ; furthermore, $k_{p,\mathrm{ij}}=k_{p,\mathrm{ji}}\mathrm{and}{k}_{g,\mathrm{ij}}=k_{g,\mathrm{ji}}$ . Here, the Hadamard multiplication is introduced for the mode shape matrix

$$\mathbf{\Phi }\left(x\right)={[{\mathrm{\Phi }}_1\left(x\right),{\mathrm{\Phi }}_2\left(x\right),\dots ,{\mathrm{\Phi }}_N\left(x\right)]}^T$$

(16)

the Hadamard multiplication defined as

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\mathbf{\Phi }}\left(x\right)=\left[\begin{array}{cccc}{\mathrm{\Phi }}_1\left(x\right)\hfill & {\mathrm{\Phi }}_1\left(x\right)\hfill & \cdots \hfill & {\mathrm{\Phi }}_1\left(x\right)\hfill \\ \multicolumn{1}{c}{{\mathrm{\Phi }}_2\left(x\right)}& {\mathrm{\Phi }}_2\left(x\right)\hfill & \cdots \hfill & {\mathrm{\Phi }}_2\left(x\right)\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathrm{\Phi }}_N\left(x\right)}& {\mathrm{\Phi }}_N\left(x\right)\hfill & \cdots \hfill & {\mathrm{\Phi }}_N\left(x\right)\hfill \end{array}\right]}\\ •{\left[\begin{array}{cccc}{\mathrm{\Phi }}_1\left(x\right)\hfill & {\mathrm{\Phi }}_1\left(x\right)\hfill & \cdots \hfill & {\mathrm{\Phi }}_1\left(x\right)\hfill \\ \multicolumn{1}{c}{{\mathrm{\Phi }}_2\left(x\right)}& {\mathrm{\Phi }}_2\left(x\right)\hfill & \cdots \hfill & {\mathrm{\Phi }}_2\left(x\right)\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathrm{\Phi }}_N\left(x\right)}& {\mathrm{\Phi }}_N\left(x\right)\hfill & \cdots \hfill & {\mathrm{\Phi }}_N\left(x\right)\hfill \end{array}\right]}^T\in {\mathbb{R}}^{N\times N}\hfill \end{array}$$

(17)

where $\tilde{\mathbf{\Phi }}(0)$ and $\tilde{\mathbf{\Phi }}(H)$ are associated with the corresponding position at the tower bottom and top, respectively. $\tilde{\mathbf{\Phi }}\prime (0)$ denotes the derivative of the mode shape with respect to x for all elements in the matrix. $[{]}^T$ is the transpose operation of the matrix. According to equation (14), the DOF to top mass in the axial direction is decoupling to other DOFs in the proposed model. It can be considered as a simple vibration damping system, where the free vibration frequency ${\omega }_d$ has a relationship with the natural frequency ${\omega }_n$ of the system defined as

$${\omega }_d=\sqrt{1-{\zeta }^2},{\omega }_n=\sqrt{1-{\zeta }^2}\sqrt{\frac{k}{m}}$$

(18)

According to the formulation derived in this article, the tower–foundation interaction and the nacelle–tower interaction can be illustrated by numerical analysis.

Natural frequency of the wind turbine system can be obtained from the free vibration status by considering the outside force of the system is zero. Therefore, we obtain the differential equation of motion of the system

$$\mathbf{M}\left(t\right)\stackrel{··}{\mathbf{q}}\left(t\right)+\mathbf{C}\left(t\right)\stackrel{·}{\mathbf{q}}\left(t\right)+\mathbf{K}\left(t\right)\mathbf{q}\left(t\right)=0$$

(19)

The equation of motion of the beam can be given by ²³

$$\begin{array}{c}\frac{{\partial }^2}{\partial x^2}\left(\mathrm{EI}\left(x\right)\frac{{\partial }^{2}y(x,t)}{\partial x^2}\right)+\frac{\partial }{\partial x}\left(M\left(x\right)g\frac{\partial y(x,t)}{\partial x}\right)\hfill \\ -\frac{\partial }{\partial x}\left(m{r}^2\left(x\right)\frac{\partial \stackrel{··}{y}(x,t)}{\partial x}\right)+m\stackrel{··}{y}(x,t)=f(x,t)\hfill \end{array}$$

(20)

where $r(x)$ is the radius of gyration, $y(x,t)$ is the transverse deflection of the beam, t is the time, and $f(x,t)$ is the applied time-dependent load on the beam.

The model including only the fundamental mode of vibration captures the essential dynamics of the system. However, the exact closed-form solution is impossible to be obtained due to the complex nature of the resulting equations. For simplicity, here, we assume constant equivalent properties. Assuming that the properties are not changing with x, the tower can be simplified as an Euler–Bernoulli beam. As an Euler–Bernoulli beam, the shear transform and inertia can be without consideration

$$\mathrm{EI}\frac{{\partial }^{4}y(x,t)}{\partial x^4}+M_{0}g\frac{{\partial }^{2}y(x,t)}{\partial x^2}+m\stackrel{··}{y}(x,t)=0$$

(21)

Assuming harmonic solution, we have

$$y(x,t)=Y\left(x\right)e^{i\omega t}=C{e}^{(i\omega t+\beta x)}$$

(22)

Then, normalizing the variable x, equation (22) can be written as

$$y(\epsilon ,t)=C{e}^{(i\omega t+\beta \epsilon )},\epsilon =\frac{x}{H}\in [0,1]$$

(23)

Substituting this in the equation of motion

$$-m{\omega }^2+\mathrm{EI}{\beta }^4-M_{0}g{\beta }^2=0$$

(24)

Solving equation (24), we have

$${\beta }^2=\frac{-M_{0}g\pm \sqrt{{\left(M_{0}g\right)}^2+4\mathrm{EIm}{\omega }^2}}{2\mathrm{EI}}$$

(25)

Assuming $\gamma =M_{0}g/2\mathrm{EI}$ and $\delta =m/\mathrm{EI}$ , then equation (25) can be written as

$${\beta }^2=\sqrt{{\gamma }^2+\delta {\omega }^2}-\gamma \mathrm{or}{\beta }^2=-(\sqrt{{\gamma }^2+\delta {\omega }^2}+\gamma )$$

(26)

Therefore, the four roots of $\beta$ can be expressed as

$${\beta }_{1,2}=\pm {\lambda }_1\mathrm{or}{\beta }_{3,4}=\pm i{\lambda }_2$$

(27)

where

$${\lambda }_1={(\sqrt{{\gamma }^2+\delta {\omega }^2}-\gamma )}^{\frac{1}{2}}\mathrm{and}{\lambda }_2={(\sqrt{{\gamma }^2+\delta {\omega }^2}+\gamma )}^{\frac{1}{2}}$$

(28)

Such that a solution of equation (22) can be given as

$$y(x,t)=(C_1{e}^{{\lambda }_{1}x}+C_2{e}^{-{\lambda }_{1}x}+C_3{e}^{i{\lambda }_{2}x}+C_4{e}^{-i{\lambda }_{2}x})e^{i\omega t}$$

(29)

with

$$Y\left(x\right)=C_1{e}^{{\lambda }_{1}x}+C_2{e}^{-{\lambda }_{1}x}+C_3{e}^{i{\lambda }_{2}x}+C_4{e}^{-i{\lambda }_{2}x}$$

(30)

which can be rewritten as

$$Y\left(x\right)=C_1{e}^{{\beta }_{1}x}+C_2{e}^{{\beta }_{2}x}+C_3{e}^{{\beta }_{3}x}+C_4{e}^{{\beta }_{4}x}$$

(31)

where $C_i,i=1,2,3,4$ , is undetermined constant.

There are four boundary conditions associated with the tower, the bending moment $P^0=0$ and shear force $Q^0=0$ at $x=0$

$$\mathrm{EI}\frac{{\partial }^{2}y(x,t)}{\partial x^2}-k_r\frac{\partial y(x,t)}{\partial x}=0{|}_{x=0}$$

(32)

$$\mathrm{EI}\frac{{\partial }^{3}y(x,t)}{\partial x^3}-M_{0}g\frac{\partial y(x,t)}{\partial x}+k_{l}y(x,t)=0{|}_{x=0}$$

(33)

And the bending moment $M^H=0$ and shear force $Q^H=0$ at $x=H$

$$\mathrm{EI}\frac{{\partial }^{2}y(x,t)}{\partial x^2}-k_{h}y(x,t)L=0{|}_{x=H}$$

(34)

$$\mathrm{EI}\frac{{\partial }^{3}y(x,t)}{\partial x^3}+M_{0}g\frac{\partial y(x,t)}{\partial x}-k_{h}y(x,t)=0{|}_{x=H}$$

(35)

where L is the distance between the nacelle center and the tower top.

Substituting equation (29) in the four boundary conditions, equations (32)–(35) result in

$$\{\begin{array}{c}\sum _{j=1}^4(\mathrm{EI}{\beta }_j^2-k_r{\beta }_j)C_j=0\hfill \\ \multicolumn{1}{c}{\sum _{j=1}^4(\mathrm{EI}{\beta }_j^3-M_{0}g{\beta }_j+k_l)C_j=0}\\ \multicolumn{1}{c}{\sum _{j=1}^4(\mathrm{EI}{\beta }_j^2-k_{h}L)e^{{\beta }_j}{C}_j=0}\\ \multicolumn{1}{c}{\sum _{j=1}^4(\mathrm{EI}{\beta }_j^3+M_{0}g{\beta }_j-k_h)e^{{\beta }_j}{C}_j=0}\end{array}$$

(36)

Arranging equation (36) with a matrix format

$$\mathbf{WC}=0$$

(37)

where

$$\begin{array}{c}\mathbf{W}=[\begin{array}{cc}\mathrm{EI}{\beta }_1^2-k_r{\beta }_1\hfill & \mathrm{EI}{\beta }_2^2-k_r{\beta }_2\hfill \\ \multicolumn{1}{c}{\mathrm{EI}{\beta }_1^3-M_{0}g{\beta }_1+k_l}& \mathrm{EI}{\beta }_2^3-M_{0}g{\beta }_2+k_l\hfill \\ \multicolumn{1}{c}{(\mathrm{EI}{\beta }_1^2-k_{h}L)e^{{\beta }_1}}& (\mathrm{EI}{\beta }_2^3-k_{h}L)e^{{\beta }_2}\hfill \\ \multicolumn{1}{c}{(\mathrm{EI}{\beta }_1^3+M_{0}g{\beta }_1-k_h)e^{{\beta }_1}}& (\mathrm{EI}{\beta }_2^3+M_{0}g{\beta }_2-k_h)e^{{\beta }_2}\hfill \end{array}\begin{array}{cc}\mathrm{EI}{\beta }_3^2-k_r{\beta }_3\hfill & \mathrm{EI}{\beta }_4^2-k_r{\beta }_4\hfill \\ \multicolumn{1}{c}{\mathrm{EI}{\beta }_3^3-M_{0}g{\beta }_3+k_l}& \mathrm{EI}{\beta }_4^3-M_{0}g{\beta }_4+k_l\hfill \\ \multicolumn{1}{c}{(\mathrm{EI}{\beta }_3^2-k_{h}L)e^{{\beta }_3}}& (\mathrm{EI}{\beta }_4^2-k_{h}L)e^{{\beta }_4}\hfill \\ \multicolumn{1}{c}{(\mathrm{EI}{\beta }_3^3+M_{0}g{\beta }_3-k_h)e^{{\beta }_3}}& (\mathrm{EI}{\beta }_4^3+M_{0}g{\beta }_4-k_h)e^{{\beta }_4}\hfill \end{array}]\hfill \end{array}$$

(38)

$$\mathbf{C}={[C_1,C_2,C_3,C_4]}^T$$

(39)

Considering the constant vector $\mathbf{C}$ cannot be zero, hence, the natural frequency can be obtained by

$$|\mathbf{W}|=0$$

(40)

This transcendental equation can be solved numerically. For the purpose of illustration, a numerical analysis and an experimental test have been compared for demonstrating the performance of the proposed model for the system.

Numerical simulation and analysis

In this article, the parameters are obtained from a 1.5-MW wind turbine YDF-1500-77. The wind turbine is a horizontal-axis upwind type, mainly consists of a three-bladed rotor which has an electric pitch-regulation, gearbox, yaw system, doubly fed induction generator (DFIG), and a partially rated converter. The wind turbine parameters used in the numerical example and field test are provided in Table 1.

Table 1.

Wind turbine parameters.

Name	Value	Name	Value
Rated power (P0)	1500 kW	Tower height (H)	73 m
Rated wind speed (vs)	11.5 m/s	Tower mean external diameter (D)	3.2 m
Cut-in wind speed (vin)	3 m/s	Tower shell thickness (th)	0.018 m
Cut-out wind speed (vout)	25 m/s	Tower mass density (ρ)	7850 kg/m3
Blade length (lb)	36.5 m	Tower’s modulus of elasticity (E)	2.1 × 1011 N/m2
Blade mass (mb)	5600 kg	Rated generator speed (vg)	1800 r/min
Blade number (Nb)	3	Generator mass (mg)	7500 kg
Hub height (Hh)	75 m	Rotor diameter (Dr)	77 m
Hub mass (mh)	33,000 kg	Structural-damping ratio (ζ)	1%
Nacelle mass (mn)	58,000 kg	Grid frequency (F0)	50 Hz

The second moment of area of tower can be obtained as

$$I=\frac{1}{64}\pi (D^4-{(D-2t_h)}^4)$$

(41)

The mass density per unit length of the tower can be obtained as

$$m=\frac{1}{2}\rho \pi ({\left(\frac{D}{2}\right)}^2-{\left(\frac{D-2t_h}{2}\right)}^2)$$

(42)

The top mass of wind turbine can be obtained as

$$M_0=m_h+m_n+3\times m_b$$

(43)

The distance between the nacelle center and the tower top can be obtained as

$$L=H_h-H$$

(44)

In general, at the fixed support, the values of the foundation stiffness parameters $k_r$ and $k_l$ will be infinity. Although the value is not given, a reasonable range can be determined. In this example, we will use $k_r$ and $k_l$ as variable parameters, where $k_r\in [1,100]\mathrm{GN}\mathrm{m}/\mathrm{rad}$ and $k_l\in [1,10]\mathrm{GN}/\mathrm{m}$ . On the same consideration, the stiffness between the top mass and the tower $k_h\in [1,100]\mathrm{GN}/\mathrm{m}$ and $k_v$ are supposed to be [0.01, 10] MN/m. Substitute these derived constants in equations (25)–(28) and (40) to obtain the first natural frequency ${\omega }_0$ . Then, the nature frequency in Hertz can be obtained as

$$f_0=\frac{{\omega }_0}{2\pi }$$

(45)

As a result, how these stiffness parameters affect the overall behavior of the system will be analyzed, and the determination values of them will be derived from the test experiment in the next section.

The proposed model has been implemented in MATLAB with the numerical solution. In order to illustrate the effects of the variable parameters in the proposed model, different values of foundation lateral stiffness $k_l\in [1,10]\mathrm{GN}/\mathrm{m}$ and rotational stiffness $k_r\in [1,100]\mathrm{GN}\mathrm{m}/\mathrm{rad}$ , and four fixed values of the top mass horizontal stiffness $k_h\in \{1.16,3.28,35.35,100\}\mathrm{GN}/\mathrm{m}$ are analyzed and shown in Figure 3. A reference natural frequency f₀ = 0.69 Hz is plotted in the subplot as a constant plane.

Figure 3.

The first natural frequency of the wind turbine with respect to the foundation lateral stiffness $k_l$ and rotational stiffness $k_r$ ; four fixed values of the top mass horizontal stiffness $k_h$ are considered, and the static response natural frequency f₀ = 0.69 Hz as a plane is also shown in the subplots: (a) k_h = 1.16 GN/m, (b) k_h = 3.28 GN/m, (c) k_h = 35.35 GN/m, and (d) k_h = 100.00 GN/m.

In Figure 3, the first natural frequency of the wind turbine is enlarged with the increasing values of the stiffness parameters. The increasing stiffness will enhance the stability of the wind turbine and get higher values of the natural frequency as the fact. The parameter $k_l$ has high sensitivity for the natural frequency. There is a fast change range for the natural frequency when the parameter $k_r$ has small values. However, the natural frequency will tend to stabilize, when the stiffness parameter $k_r$ is more than50 GN m/rad. The value of the natural frequency is not sensitive to the parameter $k_h$ . The natural frequency as a whole will be heightened when the value of $k_h$ is increasing. In Figure 4, according to equation (18), the variation of the first natural frequency of nacelle in the vertical direction with respect to the vertical stiffness $k_v$ is shown. It can be seen that higher values of $k_v$ correspond to higher values of the natural frequency.

Figure 4.

The first natural frequency of nacelle in the axial direction under various vertical stiffnesses k_v which range from 10⁴ N/m to 10⁷ N/m (with reference frequency f₀ = 5.4 Hz when k_v = 1.58 MN/m).

To validate the proposed model, we use the FAST^10,11 model as a reference. The parameters used in the FAST numerical simulation must match those specified parameters in Table 1. We build the wind turbine model (YDF-1500-77) in this article by modifying the baseline wind turbine WindPACT 1.5 MW in FAST. The natural frequencies are calculated by performing an eigenanalysis on the first-order state matrix created from the linearization analysis. The natural frequency of tower fore-aft is 0.6642 Hz, and the tower side to side is 0.6557 Hz. The results show that the agreement among the proposed model, FAST, and the experiment are quite good.

In addition, specifications of the NREL offshore 5-MW baseline wind turbine ²⁴ have been considered for validation of the proposed model. The details of the 5-MW wind turbine are provided in Table 2. The results showed that the natural frequency consistent with the FAST when the stiffness parameters of the proposed model had a reasonable range, ¹⁷ as shown in Figure 5. In particular, while the rotational stiffness k_r = 0.50 GN m/rad, horizontal stiffness k_h = 0.60 GN/m, and the lateral stiffness $k_l\in [0.5,1]\mathrm{GN}/\mathrm{m}$ , the natural frequency carried out by the proposed model was close to 0.32 Hz (the same values as FAST ²⁴ ).

Table 2.

Properties of NREL 5-MW baseline HAWT. ²⁴

Name	Value	Name	Value
Rated power (P0)	5000 kW	Tower height (H)	87.6 m
Rated wind speed (vs)	11.4 m/s	Tower mean external diameter (D)	4.875 m
Cut-in wind speed (vin)	3 m/s	Tower shell thickness (th)	0.023 m
Cut-out wind speed (vout)	25 m/s	Tower mass density (ρ)	8500 kg/m3
Blade length (lb)	61.5 m	Tower’s modulus of elasticity (E)	2.1 × 1011 N/m2
Blade mass (mb)	17,740 kg	Rated generator speed (vg)	1173.7 r/min
Blade number (Nb)	3	Rotor diameter (Dr)	126 m
Hub height (Hh)	90 m	Structural-damping ratio (ζ)	1%
Hub mass (mh)	56,780 kg
Nacelle mass (mn)	240,000 kg

Figure 5.

The first natural frequency of NREL 5-MW wind turbine under various stiffnesses: (a) k_l = 0.50 GN/m and (b) k_r = 0.50 GN m/rad.

Field test and analysis

In order to determine the natural frequency of the wind turbine used in the numerical analysis, an onsite test and the data analysis were implemented. The vibration signal of the wind turbine was acquired by four vibration transducers of the condition monitoring system. The vibration transducers were mounted on the top of the tower and the nacelle at the axial direction and the radial direction, respectively, as shown in Figure 6. The response frequency range of the vibration transducers is 0.04–1000 Hz, and the sensitivity is 960 mV/g. The sample frequency is 8 KHz, and sample length is 58 s for static status and 26 s for shock response, respectively.

Figure 6.

The experimental sketch map (four acceleration transducers mounting on the top of tower and the nacelle, two in the axial direction and the other two in the radial direction, respectively. The shock applying point on the same horizontal plane with the radial direction vibration transducer).

The detail of the wind turbine is illustrated in Table 1. At first, the static response of the wind turbine was tested when the wind speed is lower than 1 m/s, and wind turbine was stopped with rotor brake. The vibration signals for the static situation are shown in Figure 7. The subplots show the different vibration responses with respect to the nacelle and the tower. All the vibration signals’ intensity is very small in the stationary state as expected. The tower axial direction has uniform and stable vibration signal.

Figure 7.

Vibration signals of wind turbine under static response: (a) tower radial direction, (b) tower axial direction, (c) nacelle radial direction, and (d) nacelle axial direction.

Then, an impact was performed at the shock applying point, and the vibration response of the tower was recorded by the acceleration transducers as shown in Figure 8. The impact phenomenon in the vibration signal is significantly recorded by all the subplots. The tower radial direction receives the strongest shock response, and the signal’s amplitude exceeds the sensor maximum range and is recorded at the beginning of the shock activation. The vibration signal attenuates to zero rapidly in the axial direction of the tower, and the amplitude is about half of the radial direction. The impact action on the nacelle is weaker and decays faster than the tower.

Figure 8.

Vibration signals of the wind turbine under shock response: (a) tower radial direction, (b) tower axial direction, (c) nacelle radial direction, and (d) nacelle axial direction.

After that the fast Fourier transform (FFT) was utilized to analyze the frequency spectrum of the vibration signals. The result will be analyzed with the numerical model in the next section.

In the static response experiment, the vibration signals of the wind turbine in the radial direction and the axial direction are acquired under the static status. The power spectrums of the vibration signals are given in Figure 9. Note that the wind turbine first natural frequency is typically in the low frequency range, a low-pass filter with cutoff frequency of 100 Hz has been used to process the vibration signals which mean only the frequency lower than 100 Hz are analyzed. The frequency of the five largest peaks in the power spectrum is extracted. These peaks can be essentially considered corresponding to the natural frequency of the wind turbine.

Figure 9.

Power spectrums of the static response vibration signals: (a) tower radial direction, (b) tower axial direction, (c) nacelle radial direction, and (d) nacelle axial direction.

In the radial direction of the tower, most power focuses on f₁ = 0.69 Hz in the power spectrum. In the axial direction, the highest frequency is f₄ = 5.48 Hz which take more power energy in the power spectrum. Hence, the natural frequency of the tower is 0.69 Hz in the radial direction and 5.48 Hz in the axial direction based on the static response experiment. Similarly, the natural frequency of the nacelle is 5.38 and 5.43 Hz in the radial direction and axial direction respectively, as shown in Figure 9(c) and (d).

For the purpose of reducing the influence of high frequency noise, the low-pass filtering is applied to the raw signals in the shock response experiment. The power spectrums of the shock response vibration signals are illustrated in Figure 10. There is only one distinguished peak frequency in the radial direction of the tower, that is f₁ = 0.57 Hz. Hence, the natural frequency of the radial direction is 0.57 Hz from the shock response signal. The frequency of the highest peaks in the axial direction of the tower is f₂ = 5.46 Hz. The natural frequency of the nacelle is 0.68 Hz in both directions based on the shock response.

Figure 10.

Power spectrums of shock response vibration signals: (a) tower radial direction, (b) tower axial direction, (c) nacelle radial direction, and (d) nacelle axial direction.

According to the mutual influence among the wind turbine components, the natural frequency of one direction may be represented in both direction signals. Peaks of the natural frequency based on the static response and the shock response are listed in Table 3. The first-order peak is lower than 1 Hz, which can be found in all vibration signals. This indicates the first natural frequency for the wind turbine in the radial direction. This low frequency is essentially affected by the stiffness of the tower, as the shear and bending moment action on one side supported beam. The vibration in the axial direction of the wind turbine is mainly caused by the interaction between the nacelle and the tower. The highest natural frequency is focused on about 5.40 Hz in the axial direction both on the tower and on the nacelle.

Table 3.

Peaks of natural frequency based on the static response and the shock response.

Response type	Position of vibration signal		Peaks of natural frequency (Hz)
			1	2	3	4	5	6	7	8	9
Static response	Radial	Tower	0.69	−	3.41	−	5.69	7.36	−	16.67	−
		Nacelle	0.69	−	3.47	−	5.38	6.88	8.88	−	−
	Axial	Tower	0.69	1.84	3.98	−	5.48	−	10.95	−	−
		Nacelle	0.67	2.29	3.69	−	5.43	6.88	−	−	−
Shock response	Radial	Tower	0.57	−	−	−	−	−	−	−	−
		Nacelle	0.68	2.01	3.33	−	5.70	−	10.61	−	−
	Axial	Tower	0.65	−	−	−	5.46	−	9.65	15.26	18.27
		Nacelle	0.68	2.02	3.44	4.36	5.41	−	−	−	−

The maximum energy value in the power spectrum is shown in bold.

Discussion

Compared to the static response and shock response experiments, the wind turbine study in this article has a radial direction natural frequency close to 0.57–0.69 Hz and an axial direction natural frequency close to 5.38–5.70 Hz. Without loss of generality, in this article we choose f₀ = 0.69 Hz in the radial direction and f₀ = 5.40 Hz in the axial direction as the first natural frequency of the wind turbine, respectively. As shown in Figure 4, the vertical stiffness of nacelle can be determined when the natural frequency of the nacelle in the axial direction is given. According to equation (18), the top mass stiffness of the vertical direction $k_v$ will be determined. An acceptable value of $k_v$ is 1.58 MN/m when the natural frequency is 5.40 Hz.

In order to examine the effects of the stiffness in the proposed model, different combinations of the parameters have been investigated. For example, a point of intersection where k_r = 35.35 GN m/rad, k_l = 1.45 GN/m, and k_h = 1.16 GN/m and the natural frequency f₀ = 0.69 Hz is chosen to analyze the different effects of the three stiffness parameters. The variation of $k_h$ with five different values of $k_r$ and fixed $k_l$ is shown in Figure 11(a), while fixed $k_r$ with five different values of $k_l$ is shown in Figure 11(b). The natural frequency is close to the experimental values f₀ = 0.69 Hz when k_l = 1.45 GN/m and k_r = 1.16 GN m/rad. Moreover, in this condition, increased value of $k_h$ has little influence on the natural frequency. This is because as the top mass horizontal direction stiffness is more than a particular value. It can be considered as a fixed mass on the top of tower. However, when higher value of k_r = 35.35 GN m/rad is fixed, the appropriate values of $k_l$ is smaller than k_l = 1.08 GN/m. This indicates that the foundation has a specified total stiffness of the wind turbine. In the same way, two similar diagrams with variations of $k_l$ and $k_r$ are shown in Figures 12 and 13, respectively.

Figure 11.

The first natural frequency of the wind turbine with respect to the foundation lateral stiffness $k_l$ , rotational stiffness $k_r$ , and the top mass horizontal stiffness $k_h$ . The experimental natural frequency f₀ = 0.69 Hz as a line is illustrated in the diagram: (a) fixed k_l = 1.45 GN/m and (b) fixed k_r = 35.35 GN m/rad.

Figure 12.

The first natural frequency of the wind turbine with respect to the foundation lateral stiffness $k_l$ , rotational stiffness $k_r$ , and the top mass horizontal stiffness $k_h$ . The experimental natural frequency f₀ = 0.69 Hz as a line is illustrated in the diagram: (a) fixed k_h = 1.16 GN/m and (b) fixed k_r = 35.35 GN m/rad.

Figure 13.

The first natural frequency of the wind turbine with respect to the foundation lateral stiffness $k_l$ , rotational stiffness $k_r$ , and the top mass horizontal stiffness $k_h$ . The experimental natural frequency f₀ = 0.69 Hz as a line is illustrated in the diagram: (a) fixed k_h = 1.16 GN/m and (b) fixed k_l = 1.45 GN/m.

Based on the numerical simulation solutions and the experimental results, the stiffness parameters can be determined when one of the three is given. According to this, the stiffness parameters of the proposing model will be determined. However, during the wind turbine design phase, the proposed model can be utilized to improve the performance of the wind turbine by choosing suitable materials and structures with different stiffnesses to avoid the resonant frequency range. In the condition monitoring system, the natural frequency can be quickly identified among the abnormal frequency of the wind turbine system.

Conclusion

This article provides an improved dynamic analysis model in order to investigate the effects among the tower–nacelle–foundation system in wind turbine in case of incomplete structural parameters. A mathematical model of the wind turbine has been formulated using Euler–Lagrangian approach. A closed-form expression of the characteristic equation is given. Numerical simulations have been carried out using realistic data from a 1.5-MW wind turbine and a 5-MW wind turbine. The natural frequencies of the system have been carried out by the present methods and FAST. The results indicate that the presented model is able to describe the dynamic behavior of the wind turbine under incomplete structural parameters. The vibration signal of wind turbine acquired by the condition monitoring system is used to validate the proposed method. Results show that the calculated and measured natural frequencies are reasonably close.

The main aim of this study was to show how the proposed simplified model could be successfully used for analyzing the dynamic characteristic of the tower–nacelle–foundation system in the wind turbine. The modal expressions derived in this article can be used for checking the model parameters verifying and improving analytical models contrasting experimental to analytic data. Moreover, it can be utilized for estimating the dynamic property of the wind turbine in the design phase.

Further investigations are needed for considering the behavior of the blade system and environmental condition (such as the wind speed, wind direction, and temperature) for onshore and offshore wind turbines.