Investigation of damage detection in blade root joints of a 100 kW wind turbine using frequency tracking

1. Introduction

During the past decades, wind turbine technology has continuously improved and now it is the main alternative energy resource to fossil fuels and it is expected that 10 % of world electricity could be supplied by wind power by 2020 [1]. In a horizontal axis wind power system, turbine blades are used to convert wind energy to electricity. In general, turbine blades are required to be light and have small mass moment of inertia so that they can start to rotate at low wind speed. On the other hand, a wind power system has to survive severe environment for at least 20 years. During its lifetime, a wind blade may experience strong wind with speed as high as 60 m/s and in such condition, the blade should have enough strength to sustain the extreme wind load. In view of these requirements, the reliability of wind turbine blade has become an important topic of research. Although damage can occur to any component or part of the wind turbine, the most common type is rotor or blade damage and tower damage. Thus, special attention should be given to the structural health of blades because they are the elements of a wind power generation system most likely to be damaged and the cost of the blades can account for 15–20 % of the total turbine [2]. Wind turbine blades are vulnerable to cumulative fatigue damage owing to the periodic nature of the loading. Detecting damage to blades before failure is crucial. All Forces and torques applied on blades are transmitted to hub and main shaft through root joint, so blade root and consequently root joint bear maximum stress in blade. Root joints are usually made from aluminum, steel or cast iron and link blade root to hub bearing. Significant research has been carried out into the area of blade design and failure characteristics [3-5]. In many of investigation on blades damages, the root joint and blade root have been identified as the first point of initiating faults like cracks or delaminations.

Cardenas et al. [6] in their study on coupled aeroelastic damage progression for blades indicated that in their all cases, damage had been initiated in the blade root section. Mari’n et al. [7] in their study of fatigue damage detected in blades of a 300 kW wind turbine reported that Most of damages were located in the joining zone of the blade with the root where fatigue loads are more severe. Jui-Sheng et al. [8] in their analysis of wind turbine blade under critical wind loads recognized the box spar cover between the root section webs as the most likely failure location. Chen and Kam [9] in their analysis of small composite sandwich turbine blade subjected to extreme wind loads concluded that in extreme loading, failure location was at the blade root where the fibers of the composite skin ruptured.

The aim of this paper is analyzing changes in oscillatory responses of wind turbine due to damages occurred in blades root joints. All simulations carried out into the model of a 100 kW wind turbine under construction in the Sun-Air Research Institute of Mashhad (SARI) [10] according to the wind data of Binaloud site [11] and dynamic and aerodynamic characteristics of the turbine. The key properties of the turbine are presented in Table 1.

The stiffness of the root joints was reduced in real-time during these simulations. This is essentially modeling the case where the blade root joints lose stiffness due to delamination of metallic and composite layers. The responses obtained from the model simulations were then analyzed by short term Fourier transform to identify blade damage. Nagarajaiah and Varadarajan [12] have developed algorithms to track the dominant frequencies of a system using STFTs. The STFT algorithm splits up the signal into shorter time segments and an FFT is performed on each segment to identify the dominant frequencies present in the system during the time period considered. Combining the frequency spectra of each of these time segments, results in the time frequency distribution of the system over the entire time history. Reducing the stiffness of the root joint to simulate its damage will alter the frequency of vibration of the blade and nacelle and thus allow the damage to be detected.

2. Dynamic modeling

In dynamic modeling of wind turbine, flexibility of different parts and interaction between blades and tower are important. Ahlstrom [13] carried out research into the effects of increased flexibility in turbine blades and found that it can lead to a significant drop in the power output of the turbine. However, it is only over the last few years that research has started to focus on the dynamic behavior of the turbine blades and the interaction that occurs between the blades and the tower. Two main types of vibration occur in wind turbine blades, flap-wise and edgewise. Edgewise vibrations are vibrations occurring in the plane of the rotor axis. Flap-wise vibrations are vibrations occurring out of the plane of rotation of the blades. Flap-wise vibrations are considered in this paper. Ronold and Larsen [14] studied the failure of a wind turbine blade in flap-wise bending during normal operating conditions of the turbine, while Murtagh and Basu [15] studied the flap-wise motion of wind turbine blades and included their dynamic interaction with the tower. They found that the inclusion of the blade-tower interaction could lead to significant increases in the maximum blade tip displacement [16].

There are quite a number of design codes and simulation software packages available which model the structural behaviour of a wind turbine system. Lee et al. [17] gives an outline of these and the various methods used to solve wind turbine dynamics. The most common approach used by researchers in this area, however, is to formulate the equations of motion by the Lagrangian approach. This allows the derivation of coupling of the numerous flexible bodies (blades, tower, rotor shaft, etc.) and rigid bodies (i.e. yaw mechanism and control actuators) in two or even three dimensional space, directly minimizing the total energy functions of the dynamical system [18, 19]. A Lagrangian formulated model considering only flap-wise vibration is developed for the tower and rotating blades in this paper and subjected to a turbulent wind loading acting solely in the out-of-plane direction of the rotor. The proposed flap-wise system is illustrated in Fig. 1. The model includes two coordinate frames of reference, a local co-rotating system for each blade ($x$, $y$, $z$) and a global system for the combined elements which includes the tower and nacelle ($X$, $Y$, $Z$). At the root of each blade exists the origin of the local blade system. Wind loading is considered in the global $Z$-direction as significantly less loading occurs in the other directions [20]. Blades are rotating composite beams with non-uniform cross section and varying density through their length. The key properties of the blades are listed in Table 2. Root joints are massless revolute joints in flap-wise direction with torsional stiffness of $k_{root}$ which connect blades to the hub. Hub is considered rigid and is fixed to the rotor shaft. Nacelle contains all parts for converting mechanical energy to electrical energy and is considered as a solid mass located on the top of the tower. Tower is a cantilever beam with non-uniform cross section. Details of the key tower properties are outlined in Table 3.

2.1. Kinetic and potential energy

The kinetic and potential energies of the model were first derived (Eqs. (1), (2)). Total kinetic energy includes kinetic energy of three blades, nacelle and tower. Total potential energy includes potential energy due to deformation of three blades and tower in addition to deformation of root joints:

1

$T=\frac{1}{2}{{\sum }_{i=1}^3{\int }_0^L{m}_b\left(x\right){V_{bi}}^{2}dx+\frac{1}{2}{M}_{nac}{\dot{q}}_{nac}^2+\frac{1}{2}{\int }_0^h{m}_t\left(y\right)\left(\frac{\partial \omega }{\partial t}\right)}^{2}dy,$

2

$V=\frac{1}{2}{E}_b{\sum }_{i=1}^3{\int }_0^L{I}_b\left(x\right){\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}^{2}dx+{\sum }_{i=1}^3{\int }_0^{L}p\left(x\right){\left(\frac{\partial u}{\partial x}\right)}^{2}dx+\frac{1}{2}{\sum }_{i=1}^3{K}_{root}{q}_{rooti}^2+\frac{1}{2}{E}_t{\int }_0^h{I}_t\left(y\right){\left(\frac{{\partial }^2\omega }{\partial y^2}\right)}^{2}dy,$

where $m_b$ is the mass of blade, $L$ is the length of the blade, $M_{nac}$ is the mass of nacelle, $E_b$ is the Young’s Modulus for the blade, $I_b$ is the second moment of area of blade, $K_{root}$ is the stiffness of the root joint, $q_{root}$ is the generalized displacement of root joint, $u_i$is the displacement of the blade $i$, ${\dot{q}}_{nac}$ is the velocity of the nacelle, $m_t\left(x\right)$ is the mass of tower, $\omega (y,t)$ is the displacement of tower, $E_t$ is the Young’s modulus for the tower, $I_t\left(y\right)$is the second moment of area of the tower, $h$is the tower height. Also $p\left(x\right)$ is centrifugal force on blades (Eq. (3)) and $V_{bi}$ is absolute velocity of blade $i$ including the nacelle motion that causes blade tip displacement, this is a function of both the position along the blade $x$, and time (Eq. (4)):

3

$p\left(x\right)={\int }_0^L{\mathrm{\Omega }}^2{m}_b\left(x\right)\xi d\xi ,$

4

$V_{bi}^2={\left({\left(\frac{\partial \omega }{\partial t}\right)}_{(h,t)}+\left(\frac{\partial u_i}{\partial t}\right)\right)}^2+{\mathrm{\Omega }}^2{x}^2,$

where $\mathrm{\Omega }$ is the rotary speed of the rotor.

Fig. 1Wind turbine model

2.2. Assumed mode method

The motion of every substructure may be approximated by a weighted superposition of admissible functions [21]. In the assumed modes method, the solution of the vibration problem of the continuous system is assumed in the form of a series composed of a linear combination of admissible functions ${\phi }_i$, witch are functions of the special coordinates, multiplied by time-dependent generalized coordinates $q_i$(t). So:

5

$\omega (y,t)={\sum }_{i=1}^2{\phi }_t^i\left(y\right)q_t^i\left(y\right),$

6

$u(x,t)={\sum }_{i=1}^2{\phi }_b^i\left(y\right)q_b^i\left(y\right),$

where ${\phi }_t^i\left(y\right)\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}{\phi }_b^i\left(x\right)$ are admissible functions of tower and blade respectively, $q_t^i\left(t\right)$ and $q_b^i\left(t\right)$ are generalized coordinates of tower and blade respectively. The expressions of kinetic energy (Eq. (1)), potential energy (Eq. (2)) and generalized conservative forces are expressed in terms of the assumed modes solution (Eqs. (8)-(10)) and the Lagrange’s equation (Eq. (7)) are used to derive the equation of motion of the equivalent discrete system of the continuous system:

7

$\frac{d}{d}\left(\frac{\partial T}{\partial {\dot{q}}_i}\right)-\frac{\partial T}{\partial q_i}+\frac{\partial V}{\partial q_i}=Q_{i,nc},$

where $T$ is the kinetic energy of the system, $V$ is the potential energy of the system, $q_i$ is the displacement of the generalized degree of freedom (DOF) $i$ and $Q_i$ is the generalized loading for degree of freedom $i$. Also:

8

$Q_b^i\left(t\right)={\int }_0^L{f}_i\left(t\right){\phi }_b^i\left(x\right)dx,$

9

$Q_r^i\left(t\right)={\int }_0^{L}x{f}_i\left(t\right)dx,$

10

$Q_t^i\left(t\right)={\int }_0^{L}g\left(y,t\right){\phi }_t^i\left(y\right)dy.$

In witch $Q_b^i\left(t\right)$ is the generalized loading for blades, $Q_t^i\left(t\right)$ is the generalized loading for tower, $Q_r^i\left(t\right)$ is the generalized loading for root joint. $g\left(y\right)$ is the aerodynamic force on tower and $f_i\left(x,t\right)$ is the aerodynamic force on blade $i$.

The blades and tower were assumed to be vibrating in their two first modes. Mode shapes were estimated by extracting discrete points on mode shapes from frequency analysis in Abaqus and fitting polynomial functions on Abaqus data. Fig. 2, 3 show two first mode shapes of tower and blades and the path set on their nodes used for extracting mode shapes. Blade mode shapes are affected by centrifugal forces due to rotation of the rotor. Thus before extracting blade mode shapes a rotational body force proportional to the rotor speed was applied on the blade and then frequency analysis carried out. Figs. 4, 5 show fitted polynomial functions on Abaqus discrete points. Eq. (11) presents mode shape functions of the tower and blades. Eq. (12) and Eq. (13) are corresponding to the blade and the tower cross section areas and second moment of areas respectively. Constant values of following equations are listed in Appendix:

11

${\phi }_t,{\phi }_b={\sum }_{i=1}^5{a}_i{x}^{5-i},$

12

$I_b,A_b={\sum }_{i=1}^3{a}_i{e}^{-{\left(\frac{x-b_i}{c_i}\right)}^2},$

13

$I_t,A_t={\sum }_{i=1}^4{a}_i{x}^{4-i},$

The damage detection procedure in this paper only needs dominant frequency of the nacelle vibrations which are caused by the first or second natural frequencies of the tower or the blades, thus tower and each blade are modeled as one element beam which is accurate enough for dynamic modeling of the wind turbine in low frequency vibrations [22].

The model consists of 11 DOF (two in each of the blades, two for tower and one for each root joints) and by using assumed mode method can be expressed in the standard form as in following equation:

14

$\left[M\right]\ddot{\overrightarrow{q}}+\left[K\right]\overrightarrow{q}=\overrightarrow{Q},$

where $\left[M\right]$ and $\left[K\right]$ are the mass and stiffness matrices of the system respectively. $\ddot{\overrightarrow{q}}$ and $\overrightarrow{q}$ are the acceleration and displacement vectors and $\overrightarrow{Q}$ is the loading vector.

Fig. 2Two first mode shapes of blades: a) first mode, b) second mode

Fig. 3Two first mode shapes of tower: a) first mode, b) second mode

Fig. 4Blade mode shape functions

Fig. 5Tower mode shape functions

3. Wind loading

As the wind passes through a wind turbine it generates both a lift and drag force. The lift force is utilized as the main driving force of the blades and acts in the plane of rotation of the blades whereas the drag force acts in the plane parallel to the primary wind flow and therefore in the flap-wise direction of the blades. For the purposes of the model presented here, only the drag force is considered.

There are different ways for modeling aerodynamic loads on blades and tower. In this article according to the characteristics of the turbine and the installation site and according to Risø’s Guidelines for Design of Wind turbines [23] the following equation is proposed to apply the aerodynamic forces on the blades:

where $A_b$ is the cross section area of blade, $\rho$ is the air density, $V_{nac}$ is the wind speed at nacelle height, $H_{nac}$ is the nacelle height, $C_d^b$ is the blade drag coefficient, $\alpha$ is wind shear ratio, $\mathrm{\Omega }$ is the rotational speed of rotor.

Also loading on the tower was assumed to be a function of tower height as following equation [23]:

where $D_b$ is the tower bottom diameter, $D_{top}$ is the tower top diameter and $C_d^t$ is the tower drag coefficient.

Table 4 outlines numerical values for some parameters of Eq. (15), (16). Other parameters were presented in Tables 2, 3. Fig. 6 shows wind speed at nacelle height for 140 seconds taken from Binaloud site:

Fig. 6Wind speed at nacelle height of the 100 kW wind turbine (courtesy with Binaloud site)

4. Results

To obtain the simulation results, equations modeled in Matlab2010/Simulink and solved with the Runge-Kutta method. Fig. 7 shows the response of Blade 1 to turbulent load and Fig. 8 shows the frequency content. Fig. 9 and Fig. 10 show the response of the nacelle to turbulent load and the nacelle response frequency content respectively.

The model was then changed. A loss in stiffness in one root joint was simulated. This is modeling the case where the root joint loses stiffness due to delamination of blade root composite layers and root joint metallic structure. A simulation was set up for 140 seconds. Two losses in root joint stiffness of 1 % and 2 % were applied after 40 and 86 seconds of the simulation respectively. Fig. 11 shows the time history response for Blade 1. Fig. 12 shows the dominant frequency of the blade with respect to time. The STFT was applied at 30 seconds. It is clear that before the root joint stiffness was reduced the dominant frequency in the blade was about 2.5 Hz. However, after the blade’s stiffness was reduced twice the dominant frequency of the blade also reduced for 8 % and 14 %. Fig. 13 plots the nacelle displacement time history. Fig. 14 shows the dominant frequency of the system with respect to time. After that root joint stiffness was reduced two times, the dominant frequency in the nacelle reduced for 6 % and 12 %. Thus reduction of 6 % or 12 % in nacelle dominant frequency or 8 % or 14 % in blade dominant frequency led to a reduction about 10 % or 20 % in blade root stiffness respectively. So dominant frequency tracking was successful in detecting losses in root joints stiffness and hence damage for the turbine. As a result by tracking dominant frequencies of the blades and the nacelle in certain duration in turbine life time it is possible to discover root joints damages by paying attentions to reduction in dominant frequency of the nacelle or the blades.

Fig. 7Blade 1 flap-wise time history response

Fig. 8Blade 1 flap-wise frequency spectrum

Fig. 9Nacelle time history response

Fig. 10Nacelle frequency spectrum

Fig. 11Blade 1 displacement, two losses in stiffness

Fig. 12Blade 1 dominant frequency tracking

Fig. 13Nacelle displacement, two losses in stiffness

Fig. 14Nacelle dominant frequency tracking

5. Conclusion

In this article an 11 DOF model for a 100 kW wind turbine under construction in Sun-Air Research Institute (SARI) was proposed. Using assumed modes method and Lagrange equation dynamic equations of the turbine were derived. All mode shapes also derived by frequency analysis on CAD model of tower and blades of the turbine in Abaqus. In addition aerodynamic loading on blades and tower according to the characteristics of the turbine and the installation site was considered. By reducing the stiffness value of the root joints, damage was simulated. This in turn had an effect on the natural frequencies of the blades and nacelle. By tracking the dominant frequencies of blades and nacelle over time, simulated damage for the turbine was successfully detected.