# Gearbox Fault Diagnosis of Wind Turbine by KA and DRT

Authors: Mohammad Heidari
Journal: Journal of Energy
 (2016)

---

## Abstract

The spectral kurtosis analysis (KA) is used to select the filter parameters (FPs) combined with the application of the demodulation resonance technique (DRT) for a gearbox fault diagnosis (FD) of wind turbine. Based on the proposed method, the FPs can be selected automatically according to the kurtosis maximization principle. By changing of the shaft speed under the variable loads conditions, the natural frequency (NF) of the gearbox will be shifted and will affect the accuracy of the detection of the faults. So, the effect of the external loads on the NF of the gearbox is examined based on the simulation of the gearbox. In addition, the fast kurtogram (FK) combined with the demodulated resonance technology is used to process the simulated faulty signal of a gearbox. The results show that the FD of the gearbox is modified by correcting the NF shifts due to the variation of the operating loads.


---

## Body

## 1. Introduction

For a fault signal with wide frequency band, the inherent vibration of the system can be induced. That is, the fault signal with wideband (WB) will resonate the structure and the sensor itself in its natural vibration (NV) mode; it must also include the fault source signal, the NV signal of the tested structure system, and the actual vibration signal of the sensor itself. Analyzing the signals needs to choose a high frequency (HF), NV as the work goal by using the band-pass filter (BPF) to separate NB and then through the envelope detection to separate the fault information, so the system faults can be diagnosed by spectrum analysis. In the process of selecting the BPF, the center frequency (CF) of the BPF should be equal to the corresponding NF. However, this artificial method has some limitations in the real applications. Zhang et al. [1] studied the characteristics of demodulated resonance technique and applied it on train and off train failure diagnosing on power car and passenger trains bogie. Barszcz and Sawalhi [2] presented a method for FD of rolling bearings based on combining EEMD adaptive denoising with adaptive demodulated resonance. In the analysis of a real model and nondamping model, the NF and the vibration mode will be affected by the variation of mass and stiffness of the system. For the FD of gear, the changes of loads and speeds will produce dynamic rotating prestress, and this will lead to the shift of the nature frequency in the gear system. In kurtogram research, simulation and experimental verification of FD for bearing were carried out by Wang [3]; they applied the FK algorithm to the DR successfully. McDonald et al. [4] solved the problem of the difficulty of choosing the parameters in traditional DR by the method of acoustic emission in which the SVD and FK algorithm were applied to bearing FD. Heidari et al. [5] applied the spectral kurtosis (SK) method based on LMD to gear fault diagnosis; during the process, LMD has been used to obtain the different time domain distribution of signal, followed by determining the maximum kurtosis of the different channels in the time domain, according to the kurtosis maximization principle to determine the FPs. However, the effect of operating conditions in the diagnostic process is rarely considered. The diagnosis accuracy will be improved by considering the influence of the operating conditions. In this study, the influence of varying load and rotational speed on the NF of wind turbine gearbox is studied by calculating the prestress mode of the gearbox. The automatic acquisition method of BPF in CF and bandwidth (BW) in the DRT is validated by using the SK analysis method. Signal spectrum kurtosis index can reflect the signal frequency and indicate the transient impact strength, so it can achieve the effect of self-adaptation; the envelope analysis is simplified; and the diagnosis result is more accurate.


## 2. The Influence of Prestress on the Inherent Mode of Gear

To determine the initial stress using the changing in load when the prestressed model analysis is carried out and the static structure analysis is required before the mode analysis of the prestress [7], the equation can be expressed as (1)KX=F.In (1), [K] is the stiffness matrix (SM) and [F] is an external load matrix to determine the stress SM for structural analysis σ0→[S]; the modal analysis equation with prestress can be expressed as follows [8]:(2)K+S-ω12Mφ1=0.σ0 and [S] are the stress displacement based on static analysis and the prestress effect matrix, respectively. In the gearbox, the input speeds will be varied in the gear transmission system under variable load conditions. For the wind turbine, the impeller is driven by the wind energy to obtain the mechanical energy to drag the gearbox and then drive the generator shaft to rotate at high speed [9]; additional mass and rotational speed changes that were caused by wind impact will generate rotational prestress and change the NF [10]. So, for a gearbox, the transmission model of the gearbox is set up; the effects of rotational speed and other factors on the NF and vibration modes are studied under the premise of considering the inertia and the prestress of the system; the NFs and vibration modes of gear at different rotational speeds are calculated, that is, 100 RPM, 500 RPM, 800 RPM, 1200 RPM, 1500 RPM, 5000 RPM, and 10000 RPM. Table 1 shows the first 10 NFs of the gear at different speeds. Table 1 shows that the NF of each order is also varied alongside the changing of the operating speeds.Table 1
NF of the rotating system at different rotating speeds.
Speed (RPM) 0 100 500 1500 2000 2500 5000 10000 NF of each orderf/Hz 1 316.1 356.8 368.8 386.6 399.8 415.5 425.8 489.5 2 324.5 369.7 389.7 388.5 407.5 439.6 458.8 509.5 3 332.1 372.1 421.9 424.5 435.6 479.8 501.9 578.4 4 355.2 383.2 463.1 463.6 486.5 512.6 555.7 599.7 5 376.4 409.5 475.3 478.7 498.8 556.8 603.5 605.5 6 426.5 434.8 527.2 519.6 588.3 637.1 678.6 751.4 7 468.7 601.5 668.6 638.8 691.1 732.8 783.7 868.6 8 529.4 612.3 706.8 660.4 703.3 789.1 808.2 932.6 9 580.3 633.3 743.3 693.7 756.3 827.1 896.5 1031.3 10 604.4 698.3 772.5 731.5 792.5 894.1 935.5 1072.3Table1 demonstrates that the dynamic rotating prestress will lead to the variation of NF. Therefore, the influence of the variable loads on the NF of the system should be considered for the accuracy of FD results obtained. For the gearbox system, because of the random characteristics of the system, the input shaft speed and load of the gearbox are dynamic and will affect the characteristics of the NFs, so the changes of actual frequency during the analysis need to be taken into account.


## 3. RD and FK

### 3.1. RDT

RDT is a method developed from the vibration detection and analysis technology, which is based on the principle of DR for the impact FD [11]. In general, for the fault signal with WB, the intensity will be greatly improved at the NF of the structure, and the HF can be separated from the NV through a filter. Through analyzing the envelope of the signal after band-pass filtering, it can get a pulse whose frequency is consistent with the impact of failure. After the envelope processing, the signal is filtered to remove the residual HF interference signal, and the fault signal component with low frequency is reserved. Compared with the calculation of the characteristic frequency of the components, the specific reason for the failures can be judged using power spectrum analysis. The basic principle of DR is [12] as follows: (1) separating the HF vibration, because the wide frequency band signal is greatly enhanced at the resonant frequency of the acceleration sensor, with the design of a BPF whose CF is equal to the NF to separate the HF NB; (2) getting a pulse train from the filtered resonance signal by envelope demodulation whose frequency is equal to the fault signal; (3) filtering the envelope detection signal by low-pass filter (LPF), removing HF noise signal, and retaining the low frequency signal; (4) analyzing the power spectrum and fault characteristics, so as to extract the fault features. For the BPF design, the center frequency and BW need to be determined properly to obtain an accurate analysis result. However, the artificial method of determining the parameters of BPF has a great challenge and limitation. The FK algorithm can be adopted to realize the optimized parameters for the BPF design and hence could greatly improve the diagnostic efficiency and accuracy.


### 3.2. FK

The basic idea of SK is to calculate the kurtosis value for each frequency line so as to find the impact of the frequency band. The bigger the absolute value of kurtosis, the more serious the fault.The calculation equation is as follows [13]:(3)SKf=Ht,f4Ht,f22-2f≠0,(4)Ht,f=∫-∞+∞xτr∗τ-te-j2πfτdt.In (4), H(t,f) is the complex envelope of vibration signal in time-frequency domain, r(t) is the time window, SK index is zero when x(t) is a stationary signal, and x(t) is a nonstationary signal when there is noise signal b(t) interference. Equation (5) shows the calculation of SK: (5)Kx+bf=Kxf1+ρf2.In (5), ρ(f) is the reciprocal of signal-to-noise ratio (SNR) in the equation. Equation (5) shows that larger ρ(f) will lead to smaller SK index that cannot reflect the features of shock. So, the kurtogram is introduced for the analysis. The basic idea is that because Kx+bf is determined by f and the frequency resolution, this is the same as the principle of how to choose the CF and BW of the BPF in the DR. Therefore, the kurtogram algorithm is applied to the selection of the BPF parameters of the resonant demodulation. For the kurtogram algorithm, constructing a filter model whose BW is 1/4,(6)h0n=hnejnπ/4f∈0,14,(7)h1n=hnej3nπ/4f∈14,12.In (6) and (7), h0(n), h1(n), and h(n) are LPF, high-pass filter (HPF), and LPF whose cutoff frequency is 1/8, respectively. By filtering the signal x(n) by the LPF and HPF, cki(n) is the short time Fourier transform coefficient of the first i filter in the k layer that is used as input again for filtering. cki(n) is the short time Fourier transform of x(n) where CF is fi=(i+2)-12-(k+1) and frequency resolution is (Δf)k=2-(k+1) and i=(0,1,…,2-k+1).


### 3.3. Applied FK to DR

Usingcki(n) instead of H(t,f) in (1), (8)SKf=ckin4ckin2-2i=0,1,…,2k-1.The SK index becomes maximum through the FK calculation to obtain the CF and frequency resolution. SK index is used as the CF and frequency resolution of the BPF in DR, and then the fault characteristic frequency can be obtained by filtering and analysis [14].


## 3.1. RDT

RDT is a method developed from the vibration detection and analysis technology, which is based on the principle of DR for the impact FD [11]. In general, for the fault signal with WB, the intensity will be greatly improved at the NF of the structure, and the HF can be separated from the NV through a filter. Through analyzing the envelope of the signal after band-pass filtering, it can get a pulse whose frequency is consistent with the impact of failure. After the envelope processing, the signal is filtered to remove the residual HF interference signal, and the fault signal component with low frequency is reserved. Compared with the calculation of the characteristic frequency of the components, the specific reason for the failures can be judged using power spectrum analysis. The basic principle of DR is [12] as follows: (1) separating the HF vibration, because the wide frequency band signal is greatly enhanced at the resonant frequency of the acceleration sensor, with the design of a BPF whose CF is equal to the NF to separate the HF NB; (2) getting a pulse train from the filtered resonance signal by envelope demodulation whose frequency is equal to the fault signal; (3) filtering the envelope detection signal by low-pass filter (LPF), removing HF noise signal, and retaining the low frequency signal; (4) analyzing the power spectrum and fault characteristics, so as to extract the fault features. For the BPF design, the center frequency and BW need to be determined properly to obtain an accurate analysis result. However, the artificial method of determining the parameters of BPF has a great challenge and limitation. The FK algorithm can be adopted to realize the optimized parameters for the BPF design and hence could greatly improve the diagnostic efficiency and accuracy.


## 3.2. FK

The basic idea of SK is to calculate the kurtosis value for each frequency line so as to find the impact of the frequency band. The bigger the absolute value of kurtosis, the more serious the fault.The calculation equation is as follows [13]:(3)SKf=Ht,f4Ht,f22-2f≠0,(4)Ht,f=∫-∞+∞xτr∗τ-te-j2πfτdt.In (4), H(t,f) is the complex envelope of vibration signal in time-frequency domain, r(t) is the time window, SK index is zero when x(t) is a stationary signal, and x(t) is a nonstationary signal when there is noise signal b(t) interference. Equation (5) shows the calculation of SK: (5)Kx+bf=Kxf1+ρf2.In (5), ρ(f) is the reciprocal of signal-to-noise ratio (SNR) in the equation. Equation (5) shows that larger ρ(f) will lead to smaller SK index that cannot reflect the features of shock. So, the kurtogram is introduced for the analysis. The basic idea is that because Kx+bf is determined by f and the frequency resolution, this is the same as the principle of how to choose the CF and BW of the BPF in the DR. Therefore, the kurtogram algorithm is applied to the selection of the BPF parameters of the resonant demodulation. For the kurtogram algorithm, constructing a filter model whose BW is 1/4,(6)h0n=hnejnπ/4f∈0,14,(7)h1n=hnej3nπ/4f∈14,12.In (6) and (7), h0(n), h1(n), and h(n) are LPF, high-pass filter (HPF), and LPF whose cutoff frequency is 1/8, respectively. By filtering the signal x(n) by the LPF and HPF, cki(n) is the short time Fourier transform coefficient of the first i filter in the k layer that is used as input again for filtering. cki(n) is the short time Fourier transform of x(n) where CF is fi=(i+2)-12-(k+1) and frequency resolution is (Δf)k=2-(k+1) and i=(0,1,…,2-k+1).


## 3.3. Applied FK to DR

Usingcki(n) instead of H(t,f) in (1), (8)SKf=ckin4ckin2-2i=0,1,…,2k-1.The SK index becomes maximum through the FK calculation to obtain the CF and frequency resolution. SK index is used as the CF and frequency resolution of the BPF in DR, and then the fault characteristic frequency can be obtained by filtering and analysis [14].


## 4. Filtered CF Calibration

The CF of the filter which is automatically captured is close to the NF of the structure; modifying the filter CF according to the actual NF of the structure can make the results more accurate. Firstly, the quality and SM of the system are calculated when calculating the NF of the system, and SM is also divided into the time-varying SM and average SM. So, the average SM is used to find the NF of the system and the time-varying SM is used as the parameters of the system. Based on the above parameters, the NF can be calculated as in [15]. In order to measure the NF of the system more accurately, it is necessary to carry out the vibration test of the structure. Based on the measured amplitude-frequency curve, we determine the NF and damping ratio of the system. Figure 1 shows the test rig.Figure 1
The diagram of test rig.The exciter is applied to the excitation of different frequencies near the theoretical value of the NF, and the frequency is the NF of the system when the amplitude reaches the maximum value by referring to the theoretical value of NF in Figure2 [16]. The excitation device gradually increases the excitation frequency on the measured part, and, with the frequency increasing, the amplitude of the measured part is gradually increased between 300 Hz and 550 Hz, and the amplitude of the sample is gradually reduced when the excitation frequency is greater than 550 Hz. Therefore, the amplitude of the test piece reaches its maximum value at 550 Hz, which inferred the notion that the NF of the measured part is about 550 Hz.Figure 2
Vibration curve [6].The NF of the system can be measured, and we extract the HF of NV with maximum value of kurtosis as the CF for the BPF and then extract the fault characteristic frequency to judge the feature of the faults. However, the NF of the system is measured under the condition of no stress state, and the data of resonance demodulation is measured under a certain working condition with stress; it is necessary to exclude the frequency increment when choosing the FPs by resonance demodulation to get more accurate results. With the increasing of the order, the compensation value of NF is increased, and the compensation value of each step also increased with the increasing of the prestress. The greater the prestress, the greater the difference that needs to be compensated.


## 5. Simulation Study

### 5.1. Establish Fault Model

The model should be built to meet the following relationship with the typical case of 1.5 MW wind turbine gearbox: the modulus of the ring, the planet gear, and the sun gear must be the same in order to ensure that the gear train can have a correct assembly relationship [6]. In the simulation study, cut off one tooth of the planet gear to simulate the whole tooth broken fault.


### 5.2. Signal Simulation

The model is imported into ADAMS and the material properties are 45 steel, the density is7.801E-006 kg/mm3, the elastic modulus is 2.07E+005 N/mm2, and Poisson’s ratio is 0.29. Contact parameters are set as follows: IMPACT-function-based contact is used to determine the contact force and ADAMS/Solver using the IMPACT function in ADAMS library function is used to calculate the contact force; the principle of parameters is set as follows: (1) material stiffness is specified; the higher the stiffness value, the more difficult the integral solution; (2) exponent force is used to calculate the value of the contribution of the material stiffness to the instantaneous stress. Taking 1.5 or more than this value, the range becomes greater than or equal to 1. A value of 2 or even 3 is desirable for rubber, whereas for metals values are often taken from 1.3 to 1.5; (3) define the damping properties of the contact material whose range is greater than or equal to 0, usually taking from 0.1% to 1% of the value of the stiffness; (4) the penetration values are defined for full damping. The damping coefficient is zero when penetration is zero; ADAMS/Solver uses the three STEP functions to solve the damping coefficient between the two points. Its range is greater than or equal to 0. Simulation is carried out under no-load conditions, where the input speed is set to 60 RPM, the simulation time is 1 s, the step size is 0.001, the average value of the output shaft speed is 2125 RPM, theoretical output average speed is 2130 RPM, and the error rate is 0.23%. This proved the accuracy of the simulation model. The impact signal in the time domain can be obtained by the optimization process as shown in Figure 3.Figure 3
Impulse signal.The fault simulation signal is processed through adding random noise with SNR of −9 db. This presents the notion that the background noise energy is much greater than the failure impact signal, and the fault signal is submerged by the noise signal, and the actual impact of the fault line is very difficult to observe. It is more similar to the actual measured signals in the real applications. The waveform of vibration signal with noise in time domain is shown in Figure4.Figure 4
Vibration signal with noise.


### 5.3. Simulation and Verification

The signal is calculated in Figure4 by FK, and the calculated map is shown in Figure 5.Figure 5
The calculation of the kurtosis.It can be seen from Figure5 that the CF is 375.1 Hz where the spectrum kurtosis index is maximum, and the corresponding frequency resolution is 25.6 Hz. The analog signal is filtered using the 375.1 Hz and 25.6 Hz as the CF and BW of the resonant demodulation; the envelope of the filtered signal in the time domain is shown in Figure 6; and its spectrum is shown in Figure 7. Based on the spectrum of the resonance demodulation results, the fault frequency of the simulated signal can be obtained which is 27.6 Hz when the operating speed is 1500 RPM as shown in Figure 7.Figure 6
Time-frequency envelope.Figure 7
The spectrum of resonance demodulation.Figure8 shows the final demodulation spectrum at different rotation speeds of 2000 RPM, 2500 RPM, 5000 RPM, and 10000 RPM.Figure 8
FFT analysis of RD at 2000 RPM, 2500 RPM, 5000 RPM, and 10000 RPM.When the rotation speeds of the planetary gear are 2500 RPM, 5000 RPM, 2000 RPM, and 10000 RPM, respectively, the corresponding fault frequencies will be 36.5 Hz, 41.8 Hz, 86.9 Hz, and 172.3 Hz as shown in Figure8. Calculate the fault frequency again after correction of the inherent frequency of the FPs in the process of SK calculation. The fault characteristic frequency at each speed is obtained in Figure 9. It can be seen from Figure 9 that when the rotation speeds are 2000 RPM, 1500 RPM, 2500 RPM, 5000 RPM, and 10000 RPM, the corrected corresponding fault frequencies will be 25.4 Hz, 34.8 Hz, 41.9 Hz, 84.6 Hz, and 168.2 Hz. And then the theoretical value of the fault characteristic frequency at different speeds is obtained according to the fault characteristic frequency calculation equation, compared to the unmodified fault characteristic frequency, the modified fault characteristic frequency, and the theoretical value. From the results, it can be seen that the modified fault characteristic frequency curve is closer to the theoretical value. The error of the correction value and the theoretical value is 1.86%. From the above, to a certain extent, we can know that the diagnosis is more accurate, for the error between the modified frequency value and the theoretical value is smaller.Figure 9
FFT analysis of RD after frequency correction at 1500 RPM, 2000 RPM, 2500 RPM, 5000 RPM, and 10000 RPM, respectively.


## 5.1. Establish Fault Model

The model should be built to meet the following relationship with the typical case of 1.5 MW wind turbine gearbox: the modulus of the ring, the planet gear, and the sun gear must be the same in order to ensure that the gear train can have a correct assembly relationship [6]. In the simulation study, cut off one tooth of the planet gear to simulate the whole tooth broken fault.


## 5.2. Signal Simulation

The model is imported into ADAMS and the material properties are 45 steel, the density is7.801E-006 kg/mm3, the elastic modulus is 2.07E+005 N/mm2, and Poisson’s ratio is 0.29. Contact parameters are set as follows: IMPACT-function-based contact is used to determine the contact force and ADAMS/Solver using the IMPACT function in ADAMS library function is used to calculate the contact force; the principle of parameters is set as follows: (1) material stiffness is specified; the higher the stiffness value, the more difficult the integral solution; (2) exponent force is used to calculate the value of the contribution of the material stiffness to the instantaneous stress. Taking 1.5 or more than this value, the range becomes greater than or equal to 1. A value of 2 or even 3 is desirable for rubber, whereas for metals values are often taken from 1.3 to 1.5; (3) define the damping properties of the contact material whose range is greater than or equal to 0, usually taking from 0.1% to 1% of the value of the stiffness; (4) the penetration values are defined for full damping. The damping coefficient is zero when penetration is zero; ADAMS/Solver uses the three STEP functions to solve the damping coefficient between the two points. Its range is greater than or equal to 0. Simulation is carried out under no-load conditions, where the input speed is set to 60 RPM, the simulation time is 1 s, the step size is 0.001, the average value of the output shaft speed is 2125 RPM, theoretical output average speed is 2130 RPM, and the error rate is 0.23%. This proved the accuracy of the simulation model. The impact signal in the time domain can be obtained by the optimization process as shown in Figure 3.Figure 3
Impulse signal.The fault simulation signal is processed through adding random noise with SNR of −9 db. This presents the notion that the background noise energy is much greater than the failure impact signal, and the fault signal is submerged by the noise signal, and the actual impact of the fault line is very difficult to observe. It is more similar to the actual measured signals in the real applications. The waveform of vibration signal with noise in time domain is shown in Figure4.Figure 4
Vibration signal with noise.


## 5.3. Simulation and Verification

The signal is calculated in Figure4 by FK, and the calculated map is shown in Figure 5.Figure 5
The calculation of the kurtosis.It can be seen from Figure5 that the CF is 375.1 Hz where the spectrum kurtosis index is maximum, and the corresponding frequency resolution is 25.6 Hz. The analog signal is filtered using the 375.1 Hz and 25.6 Hz as the CF and BW of the resonant demodulation; the envelope of the filtered signal in the time domain is shown in Figure 6; and its spectrum is shown in Figure 7. Based on the spectrum of the resonance demodulation results, the fault frequency of the simulated signal can be obtained which is 27.6 Hz when the operating speed is 1500 RPM as shown in Figure 7.Figure 6
Time-frequency envelope.Figure 7
The spectrum of resonance demodulation.Figure8 shows the final demodulation spectrum at different rotation speeds of 2000 RPM, 2500 RPM, 5000 RPM, and 10000 RPM.Figure 8
FFT analysis of RD at 2000 RPM, 2500 RPM, 5000 RPM, and 10000 RPM.When the rotation speeds of the planetary gear are 2500 RPM, 5000 RPM, 2000 RPM, and 10000 RPM, respectively, the corresponding fault frequencies will be 36.5 Hz, 41.8 Hz, 86.9 Hz, and 172.3 Hz as shown in Figure8. Calculate the fault frequency again after correction of the inherent frequency of the FPs in the process of SK calculation. The fault characteristic frequency at each speed is obtained in Figure 9. It can be seen from Figure 9 that when the rotation speeds are 2000 RPM, 1500 RPM, 2500 RPM, 5000 RPM, and 10000 RPM, the corrected corresponding fault frequencies will be 25.4 Hz, 34.8 Hz, 41.9 Hz, 84.6 Hz, and 168.2 Hz. And then the theoretical value of the fault characteristic frequency at different speeds is obtained according to the fault characteristic frequency calculation equation, compared to the unmodified fault characteristic frequency, the modified fault characteristic frequency, and the theoretical value. From the results, it can be seen that the modified fault characteristic frequency curve is closer to the theoretical value. The error of the correction value and the theoretical value is 1.86%. From the above, to a certain extent, we can know that the diagnosis is more accurate, for the error between the modified frequency value and the theoretical value is smaller.Figure 9
FFT analysis of RD after frequency correction at 1500 RPM, 2000 RPM, 2500 RPM, 5000 RPM, and 10000 RPM, respectively.


## 6. Summary

In this paper, the effect of the rotating prestress on the NF of a wind turbine gearbox was studied, and the FK algorithm is applied to the resonance demodulation in the analysis of the impact of rotating machinery. It is shown that this method can obtain the central frequency and BW of the BPF automatically. Through the analysis, it is found that the prestress has a certain effect on the CF, which eliminated the effect of prestress on center frequency after it is automatically obtained; the analysis shows that the demodulation results are more accurate for the gearbox condition monitoring and FD.



---

*Source: 9451631-2016-10-19.xml*