Research on vibration response of a multi-faulted rotor system using LMD-based time-frequency representation

Unbalance, fatigue crack and rotor-stator rub are the three common and important faults in a rotor-bearing system. They are originally interconnected each other, and their vibration behaviors do often show strong nonlinear and transient characteristic, especially when more than one of them coexist in the system. This article is aimed to study the vibration response of the rotor system in presence of multiple rotor faults such as unbalance, crack, and rotor-stator rub, using local mean decomposition-based time-frequency representation. Equations of motion of the multi-faulted Jeffcott rotor, including unbalance, crack, and rub, are presented. By solving the motion equations, steady-state vibration response is obtained in presence of multiple rotor faults. As a comparison, Hilbert-Huang transformation, based on empirical mode decomposition, is also applied to analyze the multi-faults data. By the study some diagnostic recommendations are derived.


Introduction
In essence, fault diagnosis can be looked on as a typical task of pattern recognition, in which feature extraction is inevitably involved and plays a crucial role [1]. Therefore, it is important for fault diagnosis of rotor system to extract representative features from time-domain vibration response belonging to different fault patterns (including normal condition, different faults and their different combinations). In the present day, rotors become more flexible and operate under tight clearances and harsh environment. Thus, fault diagnosis of rotor system becomes more complicated and difficult, especially when there exist more faults than one in the system. The actual rotor system can have a varying degree of nonlinearity and could often be a multi-fault system; the assumption of single fault rotor system may be an idealistic situation. For example, a rotor system can always have some amount of unbalance, misalignment in the driveline, temporary bow, etc. However, these can be within the permissible limits. Furthermore, these initial drawbacks or faults may excite one or more faults in the system under such circumstances. For example, rotor might develop contact with the stationary parts under tighter clearances or rotor could develop fatigue crack under severe thermal and mechanical stresses [2]. These two important rotor faults, i.e., fatigue crack and rotorstator rub, could lead to catastrophic failure, if undetected. Many researches have been made on fault diagnosis of crack and rotor-stator rub. Zhou et al. [3] studied the nonlinear dynamical behavior of a horizontal Jeffcott rotor system with a transverse crack. They used a new acceleration scheme to reduce shocks to the rotor and diagnose the existence of cracks. Sabnavis et al. [4] made extensive reviews of the literature on vibration characteristics of cracked rotor shaft. Muszynska [5,6] made excellent researches and reviews on rotor-stator contacting phenomenon and rub-related diagnostics. She found contacting interaction between rotor and stator is manifested by one or more of the physical phenomena, i.e., friction, impacting, and stiffening, depending upon system parameters and operating conditions, different types of rubbing conditions prevail. Feng and Zhang [7] studied vibration response of a rotor rubbing caused by an initial perturbation. They discussed influence of various parameters on the vibration phenomena of the rotor system. Chu and Lu [8] observed very rich form of periodic and chaotic vibrations in their experimental study. They found that besides harmonic components 2× and 3×, the 1/2 fractional harmonic components such as 1/2×, 3/2×, etc., and the 1/3 fractional harmonic components such as 1/ 3×, 2/3×, etc., are also present in some of the cases. Patel and Darpe [2] pointed out that higher levels of vibration due to crack in a rotor may develop rotor-stator rubbing in tight clearance situation. Therefore, this study is aimed to examine one multi-faults case in which three faults, i.e., unbalance, transverse fatigue crack, and rotor-stator rub, are simulated together, and the effect of presence of multiple faults is studied. In recent past, only a few research studies have addressed this issue. Based on the model-based diagnostic approaches, Bachschmid and Pennacchi [9] and Platz et al. [10] studied multiple faults identification of rotor system. Darpe [11] applied the signal-based approach for analysis on the vibration response of the rotor system with crack and asymmetry faults together. Recently, Patel and Darpe [2,12] have numerically and experimentally investigated the rotor vibration characteristics of the unbalance, crack, and rotor-stator rub faults, using both classical Jeffcott rotor model and Timoshenko beam element, considering both lateral and torsional degrees of freedom.
In this study, steady-state vibration response of a multifaults (including unbalance, crack, and rotor-stator rub) Jeffcott rotor supported on simple rigid bearings is investigated, considering the fact that rotors become more flexible in the present day. Our primary objective is to investigate nonlinear and non-stationary characteristics of the rotor vibration in presence of multi-faults. This study utilizes the crack model and the rotor-stator rub model presented in [12]. The model of the rotor-stator system thus allows us to investigate the influence of unbalance, crack, and rotor-stator rub on the vibration response. Since the rub fault is highly nonlinear and transient, timefrequency representation (TFR) based on local mean decomposition (LMD) is extensively used for the study. As comparison, Hilbert-Huang transformation (HHT), based on empirical mode decomposition (EMD), is also used for analysis of the multi-faults data.
In Section 2, an outlined presentation is made on solution to unbalance vibration response of a Jeffcott rotor, in presence of crack, rotor-stator rub, or their combination. LMD-based techniques for TFR along with HHT, based on EMD, are interpreted in brief in Section 3. In Section 4, several numeric simulations are made to analyze the steady-state vibration response of the faulted rotor, and simulation results are discussed. Conclusions are given in Section 5.

Vibration response of cracked rotor with unbalance and/or rotor-stator rub
A Jeffcott rotor is considered on rigid-bearing supports having a central disk of mass m on the shaft of length L and diameter D. Transverse surface crack of depth a is assumed at the mid-span of the rotor. Figure 1a shows the coordinates (Y and Z are the fixed coordinates and ξ and h are the rotating coordinates) in the crack cross section. The unbalance eccentricity ε is assumed to be at an angle b with the weak crack direction. θ(t) is the instantaneous angle of rotation and ω is the rotor speed. Rotor-stator interaction forces are considered at disk location, as shown in Figure 2. The equations of motion for the cracked Jeffcott rotor can be written in fixed coordinate Y-Z system ( Figure 1a) as [12] mŸ + cẎ + k y Y + k yz Z = F y + mεω 2 cos(θ + β) − mg, The forces F y and F z are the nonlinear rub forces. k y and k z are direct and k yz and k zy are cross-stiffness coefficients of the cracked rotor, defined in the fixed coordinate system. These coefficients are calculated from the stiffness coefficients k ξ , k h , k ξh , and k hξ defined in the rotational frame in a cross section containing the transverse fatigue crack. In breathing crack model considered for this study, amount of open part of the crack continuously changes with shaft rotation, thereby accounting partial open/close state of the crack. Due to this partial opening and closing of the crack, cross-stiffness terms k ξh and k hξ appear in the equations of motion (Equation 1, as k ξh is related to k yz ). However, k ξh comes out to be equal to the k hξ [12]: where T is transformation matrix. The equations of motion (Equation 1) can be rewritten as [12] mŸ + cẎ 2.1. Calculation of stiffness coefficients k ξ , k hξ , and k ξh (k hξ ) The stiffness coefficients k ξ , k h , and k ξh (k hξ ) is found out from the direct and cross-flexibility coefficients g ξ ,, g h , and g ξh (g hξ ) of the cracked shaft, using linear elastic fracture mechanics theory: Castingliano's theorem, as follows [11,12]: The flexibility coefficients can be found using Castingliano's theorem using the strain energy of the cracked shaft. The strain energy of the rotor in presence of crack can be found by using the strain energies of the uncracked shaft and the additional strain energy due to the presence of crack, i.e., [13] where U 0 is the strain energy of the uncracked shaft and U c is the strain energy due to crack.  Let u i and P i are displacement and force, respectively, along the ith co-ordinate. Thus, [13] which can be written as [13] Thus, total flexibility of the cracked shaft made up of two parts: one is the flexibility of uncracked shaft and the second one is the additional flexibility introduced by the crack. The flexibility introduced by crack changes with the amount of the open part of crack. As the rotor rotates, the crack breathes. The amount of open part of the crack constantly changes, thereby changing the flexibility of the cracked rotor. The flexibilities are computed using [12] where the term L 3 /48EI represents the flexibility of the uncracked shaft, and functions F and F' are given by [12] The crack breathing condition has significant influence on the variation of stiffness coefficients, which is found using the sign of the overall stress intensity factor (SIF) [2]. Integration limits in expressions of Equation 8 are from 0 to a for depth of crack and the limits for width are not specified (Figure 1b). The limits of integration on width depend on the open part of the crack, which can be obtained from the sign of SIF along the crack edge. The negative sign of SIF indicates compressive stress and the closed crack condition and positive sign indicates tensile state of stress and open crack condition. Accordingly, depending on the instantaneous state of forces acting on the crack element, the SIF is influenced, which influences the stiffness coefficients. The detailed derivation of the flexibility coefficients of the cracked rotor can be found in [12].

Rotor-stator rub forces
When rub occurs as shown in Figure 2, the radial impact force F N and tangential rub force F T can thus be expressed as [12] where e r = (Y 2 + Z 2 ) 1/2 is the radial displacement of the rotor, δ the clearance between rotor and stator, μ the coefficient of friction, and k s the stator stiffness. The forces F y and F z can be written as (11) where ψ f is the function that decides the direction of frictional forces where R is the disc radius and υ t is the tangential velocity at disk location.

TFR based on LMD
To analyze vibration signals, a large number of methods have been proposed so far. The well-known Fast Fourier Transformation (FFT) is worthy of the most widely used tool for spectral analysis. Many FFT-based methods, for example, power spectral decomposition (PSD), cepstrum, and so on, have been applied to signal analysis in faulted rotor system. Unfortunately, these conventional techniques are often not suitable for nonlinear and nonstationary signal analysis because they always hypothesize that analyzed signal systems are linear and stationary. However, rotor faults like crack and rotor-stator rub show highly nonlinear and/or short-time transient characteristic. Such characteristics often make these conventional methods for signal analysis insufficient and sometimes theoretically inappropriate and inaccurate to reveal the nonlinear and non-stationary nature of the faults in a rotor system. Although researchers have developed some advanced time-frequency techniques, for example, short-time Fourier transformation, Wigner-Ville Distribution, and wavelet transformation, to represent the vibration signal in the form of time-frequencyamplitude (or energy) maps, all of them suffer from certain weaknesses, mainly the leakage of energy in the neighboring modes due to unequal time-frequency resolution, or the linear hypothesis on data segmentation [2].

HHT
The HHT technique, consisting of EMD and the wellknown Hilbert transform (HT), was proposed by Norden et al. [14]. Using the EMD, one can decompose any complicated vibration signal into finite number of intrinsic mode functions (IMFs). Such decomposition is completely adaptive and the extracted IMFs are almost orthogonal. Then, the HT is applied to every IMF component to produce time-frequency-energy distribution of the vibration data together, i.e., Hilbert-Huang spectrum (HHS). Theoretically, HHT technique is adaptive and able to provide equal resolution at all frequencies and time instants, which makes use of HHT more meaningful for transient vibration signals such as rotor coast up vibration response. As a powerful tool for analysis of nonlinear and non-stationary signal, the HHT technique has successfully been used in signal analysis of faulted rotor system [15][16][17].
The original signal can be represented as Here, the residue r n is left out as it is either a monotonic function or a constant. For one IMF c i (t), we can have its HT as where P is Cauchy principal value. With this definition, an analytic signal is given as in which One can write the instantaneous frequency (IF) as After performing the HT to each IMF component, the original signal can be expressed as the real part (RP) in the following form: Equation (18) gives both amplitude and frequency of each component as functions of time. Thus, time-frequency-energy map of the original signal y(t) can be obtained, which is known as HHS.

TFR based on LMD
The LMD method can be used to analyze a wide variety of natural signals such as vibration, sound, electrocardiograms, functional magnetic resonance imaging data, and earthquake data. Smith [18] thought the instantaneous amplitude (IA) and IF calculated using LMD are more stable and precise than those obtained by EMD because LMD uses smoothed local means and local magnitudes which facilitate a more natural decomposition than that using the cubic spline approach of EMD. In this way, he also illustrated how LMD gives a more localized time-frequency estimate for EEG. Wang et al. [19] declared that more accurate IF and more meaningful interpretation of the signals can be acquired by LMD than by EMD, by comparable numeric simulation. In recent past, the LMD-based technique for TFR has been used for signal analysis in fault detection of rotating machinery [19][20][21][22], and is evaluated as a more effective tool than the EMD-based method.
Essentially, the LMD scheme involves progressively separating a frequency-modulated signal from an amplitude-modulated envelope signal. This separation is achieved by smoothing the original signal, subtracting the smoothed signal from the original signal, and then amplitude demodulating the result using an envelope estimate. For an original signal y(t), the first step of the decomposition involves calculating the mean of the maximum and minimum points of each half-wave oscillation of the signal. So, the ith mean value m i of each two successive extrema n i and n i+1 is given by [18] A corresponding envelope estimate can also be derived. The local magnitude of each half-wave oscillation is given by Both the local means and the local magnitudes are then smoothed using moving averaging, weighted by the time-lapse between the successive extrema of the original signal, to form a smoothly varying continuous local mean and envelope functions m(t) and a(t), respectively. If the resulting signal is not a purely frequency modulated signal (i.e., a signal with a flat envelope), the procedure is repeated until a frequency-modulated signal is obtained.
Let the initial envelope estimate be denoted by a 11 (t), and the initial mean by m 11 (t). m 11 (t) is subtracted from the original data y(t), Then, h 11 (t) is amplitude demodulated by dividing it by a 11 (t), The envelope a 12 (t) of s 11 (t) can then be calculated. If a 12 (t) ≠ 1, the procedure needs to be repeated for s 11 (t). A smoothed local mean m 12 (t) is calculated for s 11 (t), subtracted from s 11 (t), and the resulting function amplitude demodulated using a 12 (t). This iteration process continues n times until a purely frequency-modulated signal s 1n (t) is obtained. So, where s 11 (t) = h 11 (t) a 11 (t), The corresponding envelope is given by a 1 (t) = a 11 (t)a 12 (t) · · · a 1n (t) = where the objective is that The IF can then be calculated from the frequencymodulated signal s 1n (t). The frequency-modulated signal can be written as where 1 (t) is the instantaneous phase ϕ 1 (t) = arccos s 1n (t) .
It should be noted that when calculating the phase, it is important that -1 ≤ s 1n (t) ≤ 1. From a practical point of view, the iteration process can be stopped when a 1n (t) ≠ 1. Then, any extrema of s 1n (t) which do not exactly equal ± 1 can be set equal to ± 1. Having derived the instantaneous phase, which should be set to range between ± π, the IF1 ω 1 (t) is defined from the unwrapped phase as Multiplying s 1n (t) by the envelope function a 1 (t) gives a product function PF 1 (t),

PF 1 (t) = a 1 (t)s 1n (t).
(30) This product function PF 1 (t) is then subtracted from the original data y(t) resulting in a new function u 1 (t), which represents a smoothed version of the original data since the highest frequency oscillations have been removed from it. u 1 (t) now becomes the new data and the whole process is repeated k times until u k (t) is a constant or contains no more oscillations The scheme is complete in the sense that the original signal can be reconstructed according to We can also calculate a demodulated signal spectrum analogous to the Fourier spectrum, similar to marginal spectrum of HHT, as follows where T is the length of the data and DS (ω, t) is TFR of the demodulated signal. Finally, the resulting IF and envelope values can be plotted together in the form of a demodulated signal TFR. Intuitively, we name this kind of representation as LMD-based TFR by direct solution to instantaneous frequencies (DLMD-TFR).
Different from direct computing instantaneous frequencies ω p (t), p = 1, ..., k, from the purely frequency modulated s pq (q is the iteration number of the pth times decomposition) like Equations (27) to (29), we can also have the following HT, for each PF component PF i (t), just as Equation (14): The IF can be written as After performing the HT to each PF component, the original signal can be expressed as the RP in the following form: Thus, TFR of the original signal y(t) can be, alternatively, obtained. The above representation, in this article, is named as LMD-based TFR by IDLMD-TFR.

Numeric results and discussion
Values of the system parameter in our numerical simulations are selected below, according to [12]: Mass of the disk, m = 6 kg; Shaft diameter, D = 2; R = 0.025 m; Crack depth a/D = 0.15; Shaft length, L = 0.7 m; External damping, c = 182.56 Ns/m; Stator stiffness, k s = 140 × 10 6 N/m; Coefficient of friction, μ = 0.2; Unbalance eccentricity, ε = 1.0597 × 10 -5 m. Bending natural frequency of the uncracked rotor, ω cr = 48.42 Hz; Rotating frequency of rotor shaft, ω = 0.5; ω cr = 24.21 Hz. In this section, unbalance vibration responses of the cracked rotor in presence of rotor-stator rub is discussed. Timefrequency characteristic of the steady-state rotor vibration response is identified, when unbalance, rub, and crack coexist, using both EMD and LMD-based TFR.

Vibration response of the rotor system with multifaults
For the multi-faults case (i.e., unbalance, crack, and rotor-stator rub coexist), the time-domain response in vertical direction and its FFT spectral are shown in Figure 3.
From the vibration waveform and FFT plot (Figure 3), it may be noted that the harmonics of 1× order (24.21 Hz) exists in the response. The FFT spectral plot also shows that the higher frequencies, for example 2×, 3×, 5×, etc., get excited as rub takes place. Weak subharmonic resonance at 1/2 of the bending critical speed can be observed in crack-rubbing rotor. It should be noted that the 2×, 3×, and 5× harmonics are mainly due to crack and harmonics at 4×, 6×, etc., arise due to rotorstator rub, previously illustrated in [12]. These harmonics are clearly seen in FFT plot (Figure 3b).

Time-frequency analysis of multi-faulted vibration response using EMD-and LMD-based methods
In this section, time-frequency characteristic of vibration response of the multi-faulted rotor is identified using both EMD-and LMD-based methods for TFR. The multi-faulted data are first analyzed using the EMDbased method for TFR. Figure 4 shows waveform plots of the IMFs extracted by EMD, the IFs by HT, marginal spectral, and time-frequency distribution (TFD), respectively.
Five IMF components are extracted from the multifaulted data by EMD (Figure 4a). However, it is evident that the decomposition by EMD is not stable enough, and some irregular and imprecise fluctuations exist in the extracted IMF ( Figure 4a) and IF components (Figure 4b). The main reason is that EMD extracts IAs using the cubic spline approach, rather than using more smooth local means and local magnitudes in LMD.
From the IF components, marginal spectral and HHS distribution (Figure 4b-d), the rotor-stator contact phenomenon can be found, although there exist some trivial interferences components (e.g., the embraced part by red dash line in Figure 4b), introduced by the EMD procedure. The higher frequencies are excited by contact rubbing. Accumulated energy distribution of these higher frequencies can be observed from the marginal spectral in Figure 4c. These higher frequencies disappear once the rotor loses the contact with the stator. With shaft rotation, approximately periodical occurrences of the rub can roughly be observed from Figure 4b, d. However, it can clearly be noted that there exist considerable mismatches on both time and period of rub occurrences between the IF (Figure 4b) and the HHS (Figure 4d), due to the inherent drawbacks in EMD. Thus, it is apparently impossible to determine the precisely equal ranges of rub occurrence T 1 , T 2 , T 3 , and T 4 , which may make some later procedures, such as feature extraction from the time-frequency data, difficult.
We then use two LMD-based methods for TFR, i.e., DLMD-TFR and IDLMD-TFR, to analyze the same multi-faulted data. The analysis results are shown in Figure 5 together.
Three PF components are extracted from the multifaulted data by LMD (Figure 5a). They can be further decomposed into IA and frequency components ( Figure  5b, c). Apparently, more stable decomposition is obtained by LMD, compared with that by EMD ( Figure  4a). Two marginal spectral distributions by LMD ( Figure  5e, f) also show accumulated effect of energy of different frequencies, similar to that by EMD (Figure 4c). It can be seen that the PF components show regular and precise variations of amplitude and frequency embedding in the multi-faulted data, excited by the coexisting faults.
Especially, symptoms of different faults, i.e., unbalance, crack, and rub, are clearly separated on both the IF and the TFD, computed by the DLMD-TFR method ( Figure  5c, g). According to the results by Patel and Darpe [12], three fault classes are identified and isolated with respect to different characteristic frequency bins. Unbalance is mainly localized at 1× order harmonic, crack at the 2×, 3×, and 5× frequency bins (between two dot lines in Figure 5c, g), along with the 1/2× subharmonic, and rotor-stator at 4×, 6×, etc., frequency bins (between two dash lines in Figure 5c, g). Furthermore, equal occurrence ranges of rub, i.e., T 1 , T 2 , T 3 , and T 4 , can precisely be divided into 90-260°, 450-620°, 810-980°, and 1170-1340°of shaft rotation angular positions in both the IF and the TFD, computed by the IDLMD-TFR method (Figure 5d, h). The higher frequencies disappear when there is no rub, for example, at 260-450°, 620-810°, and 980-1170°of shaft rotation, which can be isolated as occurrence ranges of unbalance and crack. The equally periodical occurrences of higher harmonics, observed in the IFs and TFDs (Figure 5c, d, g, h), may be a typical and unique rub feature and are proposed for its identification.
It should be also noted that higher level of vibrations, due to presence of crack, leads to more much interactions between rotor and stator, compared with the case of uncracked rotor with rub in [12]. In turn, contact rubbing between rotor and stator also increases level of vibration of cracked rotor without rub. Interactions among different faults, such as unbalance, fatigue crack, and rotor-stator rub, is a complex cross-coupling issue to be extensively studied in future.

Conclusions
In this article, one multiple faults case, i.e., unbalance, crack, and rotor-stator rub, in a rotor system is considered. Equations of motion of the multi-faulted Jeffcott rotor are first derived. Numerical simulation is then carried out to solve steady-state vibration response, and to comparably analyze vibration signals using the both EMD-and LMD-based methods for TFR. The study brings out the nonlinear and non-harmonic nature of the spectral components for the rotor faults. The nonlinear dynamic response due to rotor-stator interaction that appears in vibration waveform is also revealed by two LMD-based TFR methods, i.e., DLMD-TFR and IDLMD-TFR, with complete information in the form of time, frequency, and associated energy. In this way, the non-stationary nature of the rotor-stator contact is accurately captured in both IFs and TFD. Distinguishing features between rub and crack are established. This study shows that the LMD-based TFR is a more effective tool, than the EMD-based TFR, to reveal the nonlinear and nonstationary nature of the excitation frequencies produced by some rotor faults such as crack and rotor-stator rub. It is also illustrated how LMD gives a more localized time-frequency estimate for multiple rotor faults. From the investigations, relatively new diagnostic recommendations are proposed. It has been shown that rub in uncracked rotor excites higher frequency components with almost equally periodical mode and this is suggested as fault feature typical of rotor rub. When the rotor-stator rub occurs, periodical occurrences of rub can clearly be observed, and their ranges can precisely be divided equally, with respect to shaft rotation angular positions. All of these can be used for isolation of multifaults such as unbalance, crack, and rub in a rotor system.