Fault Simulation and Online Diagnosis of Blade Damage of Large-Scale Wind Turbines

Abstract

Damaged wind turbine (WT) blades have an imbalanced load and abnormal vibration, which affects their safe and stable operation or even results in blade rupture. To solve this problem, this study proposes a new method to detect damage in WT blades using wavelet packet energy spectrum analysis and operational modal analysis. First, a wavelet packet transform is used to analyze the tip displacement of the blades to obtain the energy spectrum. The damage is detected preliminarily based on the energy change in different frequency bands. Subsequently, an operational modal analysis method is used to obtain the modal parameters of the blade sections and the damage is located based on the modal strain energy change ratio (MSECR). Finally, the professional WT simulation software GH (Garrad Hassan) Bladed is used to simulate the blade damage and the results are verified by developing an online fault diagnosis platform integrated with MATLAB. The results show that the proposed method is able to diagnose and locate the damage accurately and provide a basis for further research of online damage diagnosis for WT blades.

1. Introduction

Blades are the most important component of a wind turbine (WT) and their operating status is an important factor to ensure the normal and stable operation of WTs. Since WTs are mostly located in harsh areas and the wind conditions are complex and changeable, blade faults occur more commonly with increasing operating hours. Blade faults manifest as blade damage and the reasons include fatigue due to alternating loads, lightning strokes, freezing and external impacts. If the damage is not detected and repaired in time, the blade will break or the WT may collapse and other serious accidents may occur. Therefore, online blade damage diagnosis is of great importance in terms of research value and applications [1].

The supervisory control and data acquisition (SCADA) information reflects the operating status of the WT and is easy to obtain. However, there are no data that directly reflect the state of the blades in the SCADA data. Therefore, it is important to investigate the health monitoring and fault diagnosis of the blades using intelligent algorithms such as artificial neural networks, expert systems, fuzzy logic systems and support vector machines and the WT SCADA data or vibration monitoring data. In Reference [2], a WT blade breakage monitoring method using SCADA data was proposed. A deep automatic coder (DA) model was proposed to determine impending blade damage from the SCADA data using a discriminant index. By analyzing fuzzy fault features, Yang et al. [3] detected blade faults by interpreting the data collected by the WT SCADA system; the authors used the conventional 10-minute average data in SCADA, which not only affected the diagnostic accuracy but also disregarded much of the fault training data [4].

On the other hand, it is more common to install sensors on the WT blades to monitor blade damage by using operational signals. The data directly reflect the structural damage to the blades and data processing is easy and the results are accurate. Therefore, much research has focused on the choice of sensors and signal processing. In Reference [5], an acoustic emission (AE) technique was used to obtain blade damage information and the location of the damage. Non-destructive AE methods were applied during a series of blade certification tests on a set of small WT blades [6]. However, the method requires the installation of an acoustic emitter and several receivers, which is technically difficult; in addition, the AE signals may suffer from interference from the signals from different blades or mechanical noise. In Reference [7,8], optical fiber sensors were embedded into WT blades and the blade damage was detected by processing the measured strain signal. However, optical fiber sensors can only be installed during blade manufacturing and this method is not suitable for WTs already in operation. A ceramic piezoelectric sensor mounted on the blades was used to perform a tensile test and a wavelet packet transform was used to diagnose and locate the damage [9]. Polyvinylidene fluoride resin (PVDF) film- based strain sensors were installed during a test of a full-scale WT [10]; the experimental results showed that the PVDF film-based sensors were effective for detecting damage at the trailing edge of the WT blade. In order to increase the accuracy of the damage diagnosis, the authors in Reference [11] used a hybrid sensor network consisting of capacitive film strain sensors and optical fiber sensors and using a sensor information fusion method based on a neural network. The above-mentioned experiments have achieved very good results and thin film sensors are easy to install; however, a sensor network consisting of a large number of sensors increases the installation difficulty and reduces the reliability of the system. Furthermore, it is still unknown whether the sensors affect the aerodynamic characteristics of the WT blades.

Vibration signal analysis is widely used in the field of mechanical fault diagnosis because it is a mature technology and vibration sensors are easy to install [12]. Therefore, vibration or acceleration sensors have been used commonly for damage diagnosis of blades. The authors in Reference [13,14,15] used the finite element analysis software ANSYS to simulate blade damage. The changes in the modal vibration shapes before and after blade damage were determined using vibration signal analysis. Because it is difficult to simulate WT operation processes in the ANSYS software, online blade damage diagnosis is not possible. However, the simulation results represent a good theoretical reference for damage diagnosis of blades using vibration sensors. In Reference [16], the characteristics of the vibration data before and after blade damage were determined under different environmental conditions and a principal component analysis was used for detecting the damage. Small-scale WT blades were used for damage experiments in Reference [17,18]. The change in the modal parameters was determined by analyzing the vibration signals of the blades. In Reference [19], a neural network was used to diagnose blade damage; the vibration signal data were obtained from excitation experiments. These types of methods are more common than simulation methods using ANSYS but they are still in the experimental stage. In many experiments, the vibration sensors were not actually installed on the blades. The use of a large number of sensors not only increases the installation difficulty but also affects the system reliability. In Reference [20], an exciter and several vibration sensors were installed on the blade of a Vestas V27 WT to monitor the operating conditions. This method was effective for monitoring the blade damage but requires an exciter for creating blade vibrations. In addition, harsh environmental conditions affect the reliability of the sensors.

In China, the most effective blade damage diagnosis method used in actual operations is still human observation. Therefore, a reliable and accurate method for online damage diagnosis is urgently needed. In this study, the GH Bladed software is used to simulate blade damage, which is then initially detected by analyzing the wavelet packet energy spectrum of the blade tip displacement. Unlike the SCADA system, Bladed is a reliable wind turbine simulation software that simulates blade faults accurately and obtains fault information without requiring a large amount of fault data. It has been shown that the wavelet packet method is more effective than the Fourier transform method. Moreover, the initial damage determination of the blades is also simpler and easier to implement using the wavelet packet decomposition method and the analysis does not complicate the diagnosis. When the blade is damaged, the modal parameters of the blade sections are determined using an operational modal analysis and subsequently, the modal strain energy (MSE) change ratio (MSECR) is calculated. The MSECR is used as an index to locate the damage. Computing the MSECR only requires the vibration signal of the blade, which can be measured by installing a small number of patch vibration sensors on the tip or embedding vibration sensors in the blade. Unlike acoustic emission technology and other methods of fault diagnosis that use sensors, an external signal excitation source does not have to be installed, which reduces the complexity and enhances the reliability of the system. Finally, a dynamic link library (DLL) is used to interface the Bladed and MATLAB software (2016a, MathWorks, Natick, MA, USA) to create an online diagnostic tool for WT blade damage. The simulation results show that this method can diagnose blade damage accurately. Moreover, the method can be applied to actual WTs by using a small number of patch-type vibration sensors on the blade tip.

2. Analysis and Simulation of Blade Damage Using GH Bladed

WT blades are typically constructed using fiber-reinforced polymeric composites and sandwich structures. Moreover, their geometries (e.g., the aerofoil chord length) gradually change along the pitch axis direction. It is, therefore, a challenging task to develop an accurate analytical model of the structural damage of WT blades. For simplification, in the present study, a WT blade is regarded as a multi-degree-of-freedom system consisting of 9 sections (U₁–U₉), as shown in Figure 1. The mass and stiffness can be adjusted in each section, as shown in

Table 1. Blade geometry and mass/stiffness information.

Section	1	2	3	4	5	6	7	8	9	10
Distance along pitch axis (m)	0	1.15	3.44	5.74	9.19	16.07	26.41	35.59	38.23	38.75
Chord (m)	2.07	2.07	2.76	3.44	3.44	2.76	1.84	1.15	0.69	0.03
Aerodynamic twist (deg)	0	0	9	13	11	7.8	3.3	0.3	2.75	4
Thickness (%)	100	100	64	40	30	22	15	13	13	13
Mass/unit length(kg/m)	1084.77	369.81	277.36	234.21	209.56	172.58	103.55	55.47	40.68	24.65
Edgewise stiffness (N·m2)	7.47 × 109	2.61 × 109	2.09 × 109	1.43 × 109	1.29 × 109	5.65 × 108	1.22 × 108	2.43 × 107	4,518,000	8167.51
Flapwise stiffness (N·m2)	7.47 × 109	2.43 × 109	1.41 × 109	8.34 × 108	5.56 × 108	2.09 × 108	2.95 × 107	2,259,000	113,824	3127.98

.

Due to the large size and weight of large-scale WT blades, it is difficult to conduct damage experiments and the cost is high. GH Bladed provides a simulation platform that creates an approximation of an actual WT; the software contains a wealth of blade information, including length, thickness, stiffness, quality, airfoil and other data. It can also simulate the blade icing by setting the ice position and ice density. Therefore, it is convenient to use GH Bladed for fault analysis and diagnosis of blade icing. In some cases, blade icing and damage appear to be similar but the results are different. Blade icing mainly changes the quality of the blades. However, blade damage does not change the quality of blades but causes changes in the characteristic parameters such as blade stiffness. In addition, blade icing affects the operational parameters of the WT, such as the power and motor speed [21], whereas blade damage usually does not. Blade icing faults mostly occur in the cold season and cold regions and WT maintenance personnel detect blade icing using weather and WT SCADA data. Many studies have focused on blade icing diagnosis and deicing [22,23,24]. However, the connection between blade damage and weather is insignificant and there is no effective method for its detection.

In this study, we focus on changing the section stiffness and damping to simulate blade damage and we obtain the vibration information from GH Bladed. As the blade stiffness decreases, the vibration signal gradually increases along the major axis and the response is most pronounced at the tip of the blade. Therefore, the signal obtained from the tip is used for the initial damage detection of the blade damage. Figure 2 shows the vibration displacement of section 3, section 6 and the blade tip with severe damage in section 2; the simulation results demonstrate that the vibration signal increases gradually in the direction of the main axis. Figure 3 shows the blade tip displacement of a normal blade and a blade with severe damage.

3. Signal Analysis of the Blade Tip Displacement

3.1. Fast Fourier Transform (FFT) Analysis of Blade Tip Displacement

Fourier analysis is commonly used in traditional signal analysis; it uses a fixed window function and does not reflect the non-stationary, time-domain and frequency-domain characteristics of the signals. Figure 4 shows the results of an FFT analysis of a blade tip displacement signal with different damage degrees in section 2. In the blade damage simulation, the blade parameters are assessed before and after the blade damage and the damage is divided into four classes, namely, slight damage, moderate damage, severe damage and extreme damage. For example, slight damage is defined as 80% of the original stiffness and moderate, severe and extreme damage are defined as 70%, 60% and 50% of the original stiffness, respectively. When there is little damage, the FFT of the blade tip displacement has no apparent influence on the low-frequency part and it is difficult to determine the blade damage using this method. The information in the high-frequency part is not accurate because of the limitations of the FFT analysis method. The wavelet transform is a method that uses a fixed area but a variable window size. Wavelet packet decomposition can adaptively select the frequency bands that match the signal spectrum based on certain signal characteristics [25].

3.2. Wavelet Packet Energy Spectrum Extraction

A damaged blade exhibits abnormal vibration, which can be detected in the energy spectrum. The band energy reflects the operating status of the blade. By determining the change in the energy in different frequency bands, blade damage can be initially diagnosed. The energy-based wavelet packet decomposition results are called the wavelet packet energy spectrum of the blades.

We use a 3-layer wavelet packet decomposition as an example; the wavelet packet energy spectrum extraction method consists of the following steps:

(1) Obtaining the decomposition coefficients. The blade vibration signals are decomposed using a 3-layer wavelet packet decomposition; 8 decomposition coefficients of layer 3 from low frequency to high frequency are obtained $\left(X_3^0,X_3^1,X_3^2,\cdots ,X_3^7\right)$.

(2) Reconstruction of the wavelet coefficients. We extract each sub-band range signal $S_3^j\left(j=0,1,\cdots ,7\right)$; the total signal is expressed as:

$$S=S_3^0+S_3^1+S_3^2+\cdots +S_3^7$$

(1)

(3) Calculating the signal energy of each sub-band. The reconstructed signal of each layer 3 node is expressed as: $S_3^j\left(j=0,1,\cdots ,7\right)$; the corresponding band energy is $E_3^j\left(j=0,1,\cdots ,7\right)$, which is defined as:

$$E_3^j=\int {\left|S_3^j\left(t\right)\right|}^{2}dt={\displaystyle \sum }_{k=1}^N{\left|x_{jk}\right|}^2$$

(2)

where $x_{jk}\left(j=0,1,\cdots ,7;k=0,1,\cdots ,n\right)$ represents the amplitude of the signal’s discrete points.

The blade tip displacement signals were decomposed using the wavelet packet method to obtain the characteristic information of the tip displacement energy band. The details for the energy bands are shown in

Table 2. Tip displacement energy spectrum data.

Damage Degree	Band 1	Band 2	Band 3	Band 4	Band 5	Band 6	Band 7	Band 8
Extreme damage	94.8525	2.1198	0.5638	0.1149	0.0585	0.0259	0.0139	0.0043
Severe damage	94.7138	2.1173	0.5634	0.1148	0.0585	0.0259	0.0139	0.0043
Moderate damage	94.6368	2.1126	0.5631	0.1146	0.0585	0.0258	0.0139	0.0044
Slight damage	94.5633	2.1203	0.5573	0.1141	0.0584	0.0258	0.0138	0.0043
Normal	94.2946	15.7002	1.8677	7.7063	0.4163	0.3263	0.9087	3.8432

. The results show that the values of Band 1 increase with increasing degree of damage, although the differences are not large. The same is observed for the values of Bands 2 to 8 but the difference between the damaged blade and the normal blade is large enough to determine if the blade is damaged. Unlike the FFT analysis, this method has more significant eigenvalue changes and it is easier to determine the blade damage.

This method only requires vibration sensors at the blade tips to measure the displacement signal, which is easier to achieve than sensor installation and signal analysis for other kinds of sensors. If patch-type wireless sensors are used, there is less impact on the aerodynamic characteristics of the blade, making this method very suitable for blade health monitoring in currently used commercial WTs. Depending on the blade monitoring requirements, a number of sensors can be installed in different sections during blade manufacturing to detect the presence of damage and also the damage location. We report on the blade damage location using multi-sensor vibration signals in the next section of this paper.

4. Blade Damage Location Based on Operational Modal Testing

In an experimental modal analysis, the modal parameters of a structure are obtained by the parameter identification of the system input and output signals collected under experimental conditions. Nevertheless, during the actual operation of a WT, no man-made excitation can be applied; therefore, an ambient excitation method is used, that is, the external wind condition is used as an excitation source. Subsequently, the operational mode test theory is applied to analyze the blade vibration signals. After a pretreatment using a random decrement technique, the output signal time series was analyzed using an autoregressive moving average (ARMA) model to solve for the modal parameters of each section. Finally, the blade damage location was determined by calculating the MSECR.

4.1. Operational Modal Test Method Applicable to WT Blades

The random decrement technique refers to removing or reducing random components from one or more stationary random response samples of a linear vibration system to obtain a free response signal under a certain initial excitation [26]. The following describes the basic principle of obtaining free vibration response signal data from a structure response signal using the random reduction method. For a linear system structure, the forced vibration response of a measuring point under any excitation can be expressed as:

$$y\left(t\right)=y\left(0\right)D\left(t\right)+\dot{y}\left(0\right)V\left(t\right)+{{\displaystyle \int }}_0^t\mathrm{h}\left(t-\tau \right)f\left(\tau \right)d\tau$$

(3)

where D(t) is the free vibration response at an initial displacement of 1 and an initial velocity of 0. V(t) is the free vibration response at an initial displacement of 1 and initial velocity of 0. y(0) is the initial displacement and $\dot{y}\left(0\right)$ is the initial velocity of the system vibration. h(t) is the system unit impulse response function. f (t) is the external excitation. We selected a suitable constant A to intercept a measured random vibration structure response signal y(t); a series of different intersection moments t_i (i = 1, 2, ..., N) can be obtained and the response y(t − t_i) from the moment t_i can be seen as a linear superposition of three parts: the free vibration response caused by the initial displacement at t_i, the free vibration response caused by the initial velocity at t_i and the forced vibration response caused by the random excitation f(t). Therefore:

$$y\left(t-t_i\right)=y\left(t_i\right)D\left(t-t_i\right)+\dot{y}\left(t_i\right)V\left(t-t_i\right)+{{\displaystyle \int }}_{t_i}^{t}h\left(t-\tau \right)f\left(\tau \right)d\tau$$

(4)

Since the excitation f(t) is stationary and the starting point does not affect its random characteristics, a series of starting points t_i of y(t − t_i) can be moved to coordinate the origin to obtain the subsample function x_i(t) (i = 1, 2, ..., N). That is:

$$x_i\left(t\right)=AD\left(t\right)+\dot{y}\left(t_i\right)V\left(t\right)+{{\displaystyle \int }}_0^{t}h\left(t-\tau \right)f\left(\tau \right)d\tau$$

(5)

The statistical average of x_i(t) is:

$$\begin{array}{l} x\left(t\right)=\frac{1}{N}{\displaystyle {\displaystyle \sum }_{i=1}^N}{X}_i\left(t\right)\approx E\left[AD\left(t\right)+\dot{y}\left(t_i\right)V\left(t\right)+{{\displaystyle \int }}_0^{t}h\left(t-\tau \right)f\left(\tau \right)d\tau \right]\\ \approx AD\left(t\right)+E\left[\dot{y}\left(t_i\right)\right]V\left(t\right)+{{\displaystyle \int }}_0^{t}h\left(t-\tau \right)E\left[f\left(\tau \right)\right]d\tau \end{array}$$

(6)

If the excitation f(t) is a stationary pure random vibration with a mean value of 0 and the system vibration response y(t) and $\dot{y}\left(t_i\right)$ is also a stationary random vibration with a mean of 0, then:

$$\begin{gathered} E\left[f\left(t\right)\right]=0 \\ E\left[\dot{y}\left(t_i\right)\right]=0 \end{gathered}$$

(7)

According to the above:

$$x\left(t\right)\approx AD\left(t\right)$$

(8)

where x(t) is called the free vibration response obtained by the random decrement method. Generally, turbulent wind is a natural random excitation in WTs and the mean vibration velocity of a blade tip should be zero, that is, $E\left[\dot{y}\left(t_i\right)\right]$ is zero. However, the mean value response of the blade tip displacement is not zero; therefore, the random decrement method cannot be directly applied to the original tip displacement signal [27]. Under the excitation of an average wind speed, the tip displacement is:

$${{\displaystyle \int }}_0^{t}h\left(t-\tau \right)E\left[f\left(\tau \right)\right]d\tau ={{\displaystyle \int }}_0^{t}h\left(t-\tau \right)g\left(\overline{v}\right)d\tau =B$$

(9)

where $\overline{v}$ is the average wind speed from time 0 to t, B is the mean value of the blade tip displacement,$g\left(\overline{v}\right)$ is the blade excitation function for a wind speed of $\overline{v}$. In order to ensure that this method is applicable to the tip displacement signal, the signal should be pre-processed by subtracting the mean value B so that the mean value is zero. Therefore:

$$x^{\prime }\left(t\right)=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\left[X_i\left(t\right)-B\right]\approx E\left[AD\left(t\right)+\dot{y}\left(t_i\right)V\left(t\right)\right]=AD\left(t\right)$$

(10)

According to Equation (10), a free vibration response with an initial displacement of A and an initial velocity of 0 is obtained. The response is determined by using an ARMA model time series analysis, which is a method of using parametric models to process ordered random vibration response data for modal parameter identification [28]. The relationship between the linear system excitation and the response for N degrees of freedom can be described by a higher-order differential equation, which becomes a differential equation represented by several time series at different times, the ARMA temporal model equation, in a discrete time domain:

$${\displaystyle \sum }_{k=0}^{2N}{a}_k{x}_{t-k}={\displaystyle \sum }_{k=0}^{2N}{b}_k{f}_{t-k}$$

(11)

Equation (11) describes the relationship between the response data sequence x_t and the historical value x_t−k, where 2N is the order of the autoregressive model and the sliding mean model, a_k, b_k denote the autoregressive coefficient and the sliding mean coefficient to be identified respectively and f_t denotes white noise excitation. When k = 0, let a₀ = b₀ = 1. The ARMA equation {x_t} has a unique smooth solution:

$$x_t={\displaystyle \sum }_{i=0}^{\infty }{h}_t{f}_{t-i}$$

(12)

where h_t is the impulse response function and f_t is white noise. Therefore:

$$E\left[f_{t-i}{f}_{t+\tau -k}\right]=\{\begin{array}{c}{\sigma }^2\left(k=\tau +i\right)\\ 0\left(others\right)\end{array}$$

(13)

where δ² is the white noise variance. By substituting the result of Equation (13) into Equation (12), the following is obtained:

$$R_{\tau }={\sigma }^2{\displaystyle \sum }_{i=0}^{\infty }{h}_i{h}_{i+\tau }$$

(14)

Since the linear system impulse response function h_t is the system output response when excited by a pulse signal δ_t, the expression defined by the ARMA process is:

$${\displaystyle \sum }_{k=0}^{2N}{a}_k{h}_{t-k}={\displaystyle \sum }_{k=0}^{2N}{b}_k{\sigma }_{t-k}=b_t$$

(15)

After calculating the autoregressive coefficient a_k and the sliding mean coefficient b_k, the system modal parameters can be calculated by using the expression of the ARMA model transfer function:

$$H\left(z\right)=\frac{{{\displaystyle \sum }}_{k=0}^{2N}{b}_k{z}^{-k}}{{{\displaystyle \sum }}_{k=0}^{2N}{a}_k{z}^{-k}}$$

(16)

The root of the denominator polynomial equation is solved using a high-order algebraic equation solving method and the obtained root is the pole of the transfer function. Their relationship with the system modal frequency ω_k and the damping ratio ξ_k is:

$$\{\begin{array}{c}{z}_k=e^{s_k\Delta t}=e^{\left(-{\xi }_k{\omega }_k+j{\omega }_k\sqrt{1-{\xi }_k^2}\right)\Delta t}\\ z_k^{*}=e^{s_k^{*}\Delta t}=e^{\left(-{\xi }_k{\omega }_k-j{\omega }_k\sqrt{1-{\xi }_k^2}\right)\Delta t}\end{array}$$

(17)

The modal frequency ω_k and the damping ratio ξ_k can be obtained from Equation (18), that is:

$$\{\begin{array}{c}{R}_k=\ln z_k=s_k\Delta t\\ {\omega }_k=\frac{\left|R_k\right|}{\Delta t}\\ {\xi }_k=\sqrt{\frac{1}{1+{\left(\frac{Im\left(R_k\right)}{Re\left(R_k\right)}\right)}^2}}\end{array}$$

(18)

Suppose that the k-order residue of H_pq (s), which is the transfer function of the p-point response excited at point q, is A_kpq, then the residue can be calculated as follows:

$$A_{kpq}=\underset{z\to z_k}{lim}{H}_{pq}\left(z\right)\left(z-z_k\right)=\frac{{{\displaystyle \sum }}_{k=0}^{2N}{b}_k{z}^{-k}}{{{\displaystyle \sum }}_{k=0}^{2N}{a}_k{z}^{-k}}\left(z-z_k\right)|z=z_k$$

(19)

The modal vector can be obtained by processing the residue obtained from a series of measured response points. For a structure with n response points, the first task is determining the measurement point with the largest absolute value from the residue of n corresponding modes of the same order. Assuming that the point is the measurement point m, the normalized complex modal vector corresponding to the k-order mode can be obtained by the following formula:

$$\left\{{\varphi }_k\right\}={\left[A_{K1q}{A}_{k2q}\cdots A_{kmq}\right]}^T/A_{kmq}$$

(20)

In this way, the modal parameters such as the modal frequency, modal damping and modal vibration mode can be obtained.

4.2. Blade Damage Location Analysis Based on the MSECR

Because a blade can be simplified as a hollow cantilever beam structure, a model of the blade sections was developed based on the structural characteristics and material properties (Figure 5). Three neighboring sections are denoted as U_n−1, U_n and U_n+1 with the masses of m_n−1, m_n and m_n+1, respectively. The sections U_n−1 and U_n are connected via stiffness k_n−1,n and damping c_n−1,n and the sections U_n and U_n+1 are connected via k_n,n+1 and c_n,n+1. Consequently, when an external load is applied to the blade, the dynamic response can be expressed by the following equation:

$$M\ddot{x}+C\dot{x}+Kx=F\left(t\right)$$

(21)

where x represents the vector of the displacement responses along the blade; M, C and K denote the equivalent structural mass, the equivalent structural damping and the equivalent structural stiffness matrices, respectively; F is the matrix of the external forces. It can be inferred that when a local defect occurs in section n, the values of c_n−1,n, c_n,n+1, k_n−1,n and k_n,n+1 change correspondingly, whereas the damping and stiffness in the other sections may not change.

This configuration was the inspiration for developing a damage location diagnosis method for WT blades using the modal analysis of the sections. Blade icing has a significant influence on the mass characteristic parameters of the blade but has little effect on stiffness and damping. Since the values of the damping and stiffness are dependent only on the blade’s structural integrity, the proposed damage location diagnosis method responds only to changes caused by structural damage. Therefore, false alarms due to ice on the blade surfaces can be avoided.

The damage to the blade structure has nothing to do with its mass, that is to say, $\left[\Delta M\right]=0$. Therefore, structural damage is equivalent to a change in stiffness. Equation (22) shows the relationship between the structural stiffness and the modal vibration modes before and after the damage has occurred (superscript d indicates the damaged section):

$$\begin{array}{c}\left[K^d\right]=\left[K\right]+\left[\Delta K\right]=\left[K\right]+{{\displaystyle \sum }}_{j=1}^L{\partial }_j\left[K_j\right]\\ \left\{{\varphi }_i^d\right\}=\left\{{\varphi }_i\right\}+\left\{\Delta {\varphi }_i\right\}=\left\{{\varphi }_i\right\}+{\displaystyle {\displaystyle \sum }_{j=1}^n}{c}_{ij}\left\{{\varphi }_j\right\}\end{array}$$

(22)

In Equation (22), $-1<{\partial }_j<0$ and $0<c_{ij}<1$. In addition, the MSE is related to the mode shape. The i-order MSE of the j-th section before and after structural damage is defined as follows:

$$\begin{array}{c}MS{E}_{ij}={\left\{{\varphi }_i\right\}}^T\left[K_j\right]\left\{{\varphi }_i\right\}\\ MS{E}_{ij}^d={\left\{{\varphi }_i^d\right\}}^T\left[K_j\right]\left\{{\varphi }_i^d\right\}\end{array}$$

(23)

Usually, only a low-order modal term is calculated and the high-order modal terms are ignored; the MSE of the section before and after structural damage is defined as follows:

$$MSE{C}_{ij}=MS{E}_{ij}^d-MS{E}_{ij}=2{\left\{{\varphi }_i\right\}}^T\left[K_j\right]\left\{\Delta {\varphi }_i\right\}$$

(24)

This analysis indicates that the damage location can be diagnosed by the index vectors obtained from the MSE. However, this is not sufficient to determine the damage degree. The MSECR is defined as the indicator of the damage degree:

$$MSEC{R}_{ij}=\frac{|MS{E}_{ij}^d-MS{E}_{ij}|}{MS{E}_{ij}}$$

(25)

where MSECR_ij is the MSECR for the j-th section with respect to the i-order mode.

In order to reduce the influence of the experimental modal random noise, multiple different order modalities can be used to diagnose the location of structural damage.

$$MSEC{R}_j=\frac{1}{m}{\displaystyle \sum }_{i=1}^{m\sum }\frac{MSEC{R}_{ij}}{MSEC{R}_{imax}}$$

(26)

where m is the number of used modalities and MSECR_imax is the maximum MSECR [25].

The theoretical analysis shows that when the change in the stiffness matrix and the modal parameters before and after structural damage are used as input diagnostic information, the location of the blade damage can be obtained using the MSECR.

5. Simulation Verification

The GH Bladed software possesses high-precision WT modeling and simulation functions but the data post-processing capability is limited. In contrast, the MATLAB software has strong data analysis and post-processing capability and a rich algorithm toolbox. However, the nonlinear model of the WT created in MATLAB has low accuracy and the simulation function is also very limited. In order to verify the proposed series of blade damage diagnosis methods and perform online diagnostic simulation, the two software applications have to be combined into a fault simulation platform. Because there is no direct data communication interface between the two applications, the data interaction between GH Bladed and MATLAB was implemented by using the Bladed external DLL interface file. The schematic diagram of the integration is shown in Figure 6. The fault data were exported to MATLAB through the DLL file and the signal analysis and damage diagnosis were performed. Because the Bladed external DLL data interface allows for setting a data transmission time interval, the data transmission mode is very similar to the data acquisition process between the SCADA system and the actual sensors. Moreover, Bladed has a sensor simulation function, in which the sensor characteristics such as the signal noise, delay and fault can be specified. This greatly increases the simulation accuracy and credibility of the online monitoring and diagnosis simulation performed in this study.

The experimental WT data used in GH Bladed are shown in

Table 3. Parameters of the experimental WT.

Number	Parameter Name	Value
1	Length of blade	38.75 m
2	Rated wind speed	12 m/s
3	Rated power	2 WM
4	Above-rated generator speed set-point	1500 rpm
5	Transmission ratio	83.33
6	Minimum generator speed	850 rpm
7	Pitch angle range	−2–90°
8	Height of tower	60 m
9	Rated generator torque	13,403 Nm
10	Maximum generator torque	14,400 Nm
11	Air density	1.225 kg/m3
12	Cut-in wind speed	4 m/s
13	Cut-out wind speed	25 m/s
14	Number of blades	3
15	Rotor diameter	80 m

. Figure 6 shows the flowchart of the blade damage diagnosis. The specific steps were as follows: First, the wavelet packet decomposition of the blade tip displacement signal was conducted to obtain the characteristic information of the energy bands. The blade damage was preliminarily determined based on these characteristics. For the damaged blade, the modal parameters of each section were calculated using the operational modal method and the MSECR was obtained. The MSECR represented the index to determine the location of the damage.

When the wavelet packet energy of the blade tip displacement signals is abnormal, the displacement signals of each section are analyzed using the operational modal analysis method to obtain the modal parameters of every section and the MSECR. Figure 7 shows the flowchart of the blade damage diagnosis.

Table 4. Modal parameters of each unit.

Blade Unit	First-Order Frequency	First-Order Mode	Second-Order Frequency	Second-Order Mode
1	0.3201	−0.005207	3.8343	0.000059
2	0.2986	−0.002834	3.7893	−0.000078
3	0.2794	0.017348	3.9138	−0.000144
4	0.2677	0.073219	4.0817	0.000130
5	0.2472	0.304202	4.3765	−0.000204
6	0.2036	1.419367	4.3478	−0.001643
7	0.2349	1.096382	3.3818	−0.030174
8	0.2445	1.043385	3.0982	−0.011586
9	0.2464	1.026799	3.1026	−0.013261

shows the modal parameters of section 2 with severe damage. The MSECR was calculated using the data in Table 4 and the damage location was determined using MSECR as the discriminant index. Figure 8 shows the MSECR histograms of section 2 before and after the severe damage occurred at a mean wind speed of 12 m/s.

In order to further illustrate the effectiveness of this method, a large number of simulation experiments were conducted using blade damage in different locations and different degrees of damage. Figure 9 shows the MSECR histograms for the different locations and degrees of damage for different wind speeds. The blade damage location can be clearly detected in Figure 9. Figure 10 shows the MSECR histograms of when two sections were damaged.

Figure 9 shows the MSECR of each blade section when a single section is damaged. The simulation results show that the MSECR of the damaged section has a prominent peak, which is significantly higher than those of the undamaged sections. The MSECR of the damage section in Figure 9 is about five times as large as those of the normal sections. In addition, the MSECRs of the adjacent sections were also slightly higher than (about 50% higher than the normal sections) and there was no significant change in the MSECR of the distant sections. Therefore, the damaged section can be determined based on the MSECR. Figure 10 shows the MSECR of each blade section when two sections are damaged. It is observed that the MSECR of the damaged section is higher than that of the normal section. It is noteworthy that the MSECR provided a good indication of the degree of damage if the damage degree of the two damaged sections was large (as shown in Figure 10a). However, when the damage degree was significantly different for the two sections, the sections with the smaller damage sections may be difficult to detect (as shown in Figure 10b). As a result, it is possible to miss some minor damage when using this method for blade damage location but this does not affect the diagnosis and location of the major damage when the blade is damaged in multiple sections. Generally, blade damage consists mostly of single-section damage. And the more serious parts should be mainly considered in multiple-section damage. Environmental conditions affect the structural characteristics of the blade and the calculated value of the wavelet energy spectrum will differ under different weather conditions; however, the location of the damage is not affected. Therefore, this method can be used to locate blade damage effectively.

6. Conclusions

The proposed method for blade damage diagnosis and damage location based on the tip displacement signal, wavelet packet decomposition and operational modal analysis has the following characteristics:

(1) The wind power simulation software Bladed was used to simulate the blade damage fault of a wind turbine by determining the change in the structural characteristic parameters before and after blade damage. An integrated MATLAB and Blade simulation platform was developed for real-time online damage diagnosis.

(2) An initial on-line blade damage assessment was conducted using Fourier analysis and wavelet packet energy spectrum analysis of the tip displacement. The simulation results showed that the wavelet packet energy spectrum analysis was not only easy to implement but also provided significantly better results than the traditional Fourier analysis.

(3) A method for identifying the working modal parameters of the blades was proposed by combining the random decrement method and ARMA model parameter identification. The damage was accurately located by calculating the MSECR of each blade section without requiring additional excitation signals and the method proved effective for different degrees of damage occurring simultaneously.

Although the results of this study are based on software simulation results, the Blade software provides high accuracy and is suitable for preliminary assessment prior to empirical research. In a future study, we will build a small-blade experimental wind turbine and develop a wireless patch sensor for empirical research. The results of this study provide methodological guidance for system implementation.