Study on the diagnosis of rub-impact fault based on finite element method and envelope demodulation

2.1. Modal analysis

Since the rub-impact occurred between the blade and the casing can lead to the unsteady motion of the blade and even arose the high frequency vibration modes of casing, it is significant to conduct the modal analysis of casing. For achieving the finite element analysis about the casing, meshes are created utilizing the ANSYS software. Finite element meshing of the casing surface and the pedestal are done with three-dimensional, 8-node tetrahedral structural solid elements (SOLID 185).

Casing material is considered with the isotropic elasticity. The casing material properties are assigned as follows.

Young’s Modulus: 2.09×10⁵ MPa;

Poisson’s ratio: 0.269; Density:7890 kg/m³.

Boundary condition of meshed model are assigned as all nodes of the whole pedestals’ bottoms are fixed and the contact surfaces between the casing and the pedestals are bonded, and there are a total of 41533 elements, as shown in Fig. 1.

After assigning both bottom faces of the pedestal and casing material properties for the analysis, the calculation has been done through FE based model analysis solver of ANSYS.

The following natural frequencies and mode shapes obtained from the FEM analysis are shown in Table 1 and Fig. 2, respectively. Only the first six orders frequencies are extracted in this paper, since the vibration amplitude of the high frequency is relatively small.

Fig. 1Finite element mesh chart of casing and pedestal

Fig. 2Natural frequency and mode shapes for casing and pedestal

a) Mode 1

b) Mode 2

c) Mode 3

d) Mode 4

e) Mode 5

f) Mode 6

2.2. The establishment of the contact force model

In the section, we fully consider that elastically deformable radial blades of outer radius ‘$R$’, are installed on a rigid disk of hub radius ‘r’, as illustrated in Fig. 3. The casing interaction forces upon the blade tip is acquired by means of Coulomb friction model. Based on the experimental conditions, it is supposed that the rub-impact only occurs in the radial direction. The individual blades behave like a cantilever beam-column of span length ‘$L$’ and are subjected to a centripetal force field generated due to the rotor spin velocity ‘$\omega$’ with the blade tip rub load $F$ (+sign: tension and – sign: compression) acting along the beam axis. The tip rub force $F$ denotes a dynamic contact load, which is non-zero only under rub-impact status and will always have a (−) sign because of the contact load being compressive in nature. Let μ denotes the friction coefficient.

Fig. 3Rub forces F and μF on the rotor blade

Fig. 4Schematic view of the rub zone on the casing

Pulse loads [11-16] refer to forces that are applied over a short amount of time, less than 100 ms or so. Although that seems like a very short time window, it can easily represent 30 or more revolutions in some rotating machines, certainly long enough to include an entire rub event. The stationary global Cartesian unit vectors $(i,j,k)$ is used to describe the position of the blades. The parameter $∆t$ is defined as the time width in time, as shown in Fig. 4. The backward whirl induced by dry friction is ignored because the rotating speed is relative small and the casing is typical thin-walled structure. It is assumed that every blade has the same physical and geometric parameters, and rubbing occurs in the same position, the pulse force can be expressed as:

where $A$ denotes the amplitude of pulse force, n is the number of blade, $\omega$ represents the angular velocity of rotor.

The finite element method (FEM) is chosen to perform the time-accurate or transient analysis because it includes an explicit based solver that is suitable for short-duration loading events such as impacts.

2.3. Transient dynamics analysis of casing

Transient dynamics analysis is used to identify the structure dynamic characteristics in response of the load changing over time. The same boundary condition is adopted when the transient dynamic analysis is carried out. Periodic very short pulse force will be imposed on the casing when the rubbing fault occurs between the rotor and the casing. According to the experiment conditions, the amplitude of rubbing force acting on inner surface of the casing is set to 50 N and the rubbing cycle is defined to 0.0025 s, as shown in the left side of Fig. 5. The number of node which the rubbing force is imposed on is 11045, as shown in the right side of Fig. 5. The 11045th node’s $Y$-component of vibration acceleration response shown in Fig. 6 is obtained by using ANSYS to solve.

It can be seen from Fig. 6 that pulse rubbing force generated by rubbing can arose the casing’s high-frequency inherent vibration mode, which is taken as the carrier. The amplitude of carrier is subjected to periodic rubbing pulse modulation, so the casing vibration mode shows complex amplitude modulation phenomenon. Meanwhile, the contact and separation process often generated impacts between the rotor and casing. Therefore, the dynamic response of casing under rubbing have obvious impact characteristics. It can be observed from Fig. 6 that the amplitudes under rub-impact condition roughly present a downward trend, which is due to the attenuated vibration of corresponding node on casing in the short time after rub-impact. In order to extract the rub-impact feature, further investigate the amplitude-modulated and amplitude variation characteristics, it is essential to carry out the blade-casing rubbing experiments and collect the vibration acceleration signals from the casing.

Fig. 5a) Periodic rubbing force; b) force diagram of the casing

a)

b)

Fig. 6Y-component of vibration acceleration response

3.1. Experiment arrangement of the rotor system

In order to find out the laws of blade-casing rubbing fault, a simplified rotor-bearing-casing test bed system is carried out for the sake of simulating the rubbing fault. The schematic of the rotor-bearing-casing test bed is shown in Fig. 7. The picture of the rotor-bearing-casing test rig and the sketch of data acquisition system are illustrated in Fig. 8(a) and 8(b), respectively. The framework of the rotor-bearing-casing system are composed of a casing, one disk mounted on a shaft, and the shaft is supported by double row angular contact ball bearing. A motor drives the shaft to rotate by a flexible coupling.

Fig. 7Schematic of the rotor-bearing-casing test bed

The rotor-bearing-casing system is driven by the motor with the speed of 2100 rpm. The radial rubbing status is simulated by adopting a rub screw to hit the blade tip. The vibration acceleration signals of casing are measured by the acceleration sensors. An eddy current sensor is used to obtain the rotational speed by monitoring the key phase. The data acquisition system is utilized to capture the time series. The sampling frequency of the data acquisition device is set to be 4 kHz.

Fig. 8Experiment setup of rotor-bearing-casing system

a) Test rig for rubbing between the blade and casing

b) Data acquisition software and test rig monitoring interface

4.1. Extraction of signal envelopes

A large amount of variation component in signals could be embodied in the envelopes of signals, and the carrier wave of envelopes is usually taken by the high-frequency components of vibration signal. Hilbert technique is an effective tool to obtain the envelopes of signals in fault diagnosis [20, 21].

As for a real signal $x\left(t\right)$, $\mathcal{H}\left\{x\left(t\right)\right\}$ represents the Hilbert transform of $x\left(t\right)$, which is defined using the Cauchy principal value (denoted by $P$) as:

The analytic signal $g\left(t\right)$ of the real signal can be obtained which is defined as:

where $\widehat{x}\left(t\right)$ denotes Hilbert transform of the real signal $x\left(t\right)$, and the amplitude of analytic signal$$is expressed as:

This is the envelope of $x\left(t\right)$.

4.2. The application of envelope demodulation method to the rubbing fault diagnosis of the rotor and the casing

Based on the basic characteristics of the casing vibration signal caused by rubbing, the envelope demodulation method is utilized to separate the high frequency inherent vibration signal of the casing. Specific procedures are shown in Fig. 11.

The envelope analysis results of original casing acceleration signal are shown in Fig. 12, respectively. It can be found that the rubbing fault information can be extracted from the strong interference signal, such as background signal and noise signal, when the envelope analysis of the original signal under rub-impact status is carried out.

Fig. 11Specific procedures

Fig. 12Envelope analysis of casing signal under rubbing

a) Time domain after bandpass filter

b) Envelope of time domain

c) Time domain after low-pass filter

d) Spectral analysis after filtering

4.3. Rubbing features extraction based on EMD

As aforementioned, the responses of the system present nonlinear phenomenon and some faults almost have similar frequency characteristics, which to some extent increases the complexity of fault detect just depend on tradition features. Thus, the method of empirical mode decomposition (EMD) is applied to further process the nonstationary and nonlinear vibration acceleration signals and excavate more fault information that cannot be detected just relying on statistical parameters in time domain and frequency domain. The novel tool, EMD, first developed for data analysis by Huang et al [22], is a kind of adaptive and unconstrained decomposition of the complex signal into a finite number of almost orthogonal intrinsic mode functions (IMFs) following a recursive process and has been widely used for various applications [23-26]. An IMF is supposed to satisfy the following two conditions:

1) In the whole data set, the number of extrema and the number of zero crossings should be equal or different no more than one.

2) At any point, the mean value of the envelope defined by the local minima and local maxima is zero.

Sifting is implemented iteratively for extracting IMFs from original signal using following algorithm:

1. Let $m_1$ represent the mean of upper and lower envelope value of original signal $x\left(t\right)$, which are determined from a cubic-spline interpolation of local maxima and minima. The first component, $h_1\left(t\right)$ is calculated as shown in Eq. (5):

2. In next step, $h_1\left(t\right)$ is considered as the original signal, and $m_{11}$ is the mean of upper and lower envelope value of $h_1\left(t\right)$ and $h_{11}\left(t\right)$ is obtained:

3. Above procedure is repeated $n$ times, until $h_{1n}$ conforms to the conditions of an IMF. Then it is designated as first IMF, $c_1\left(t\right)=h_{1n}\left(t\right)$, it is then separated from rest of the data Eq. (7):

4. Now $r_1\left(t\right)$ is considered as main signal and steps 1-3 are repeated for calculating second IMF.

5. The number of IMFs that can be extracted depends on the original signal. The stopping criterion is that the $r_m\left(t\right)$ becomes monotonic function from which no more IMF can be extracted.

Therefore, based on the definition of EMD approach, any complicated signal could be theoretical expressed as:

where $c_k\left(t\right)$ denotes the $k$th IMF of the original signal $x\left(t\right)$ and $r_m\left(t\right)$ represents the ultimate residue value, which could be the signal mean trend or a constant. On account of the two conditions of each IMF as mentioned above, the IMF may be a single component of frequency-modulated (FM) or amplitude-modulated (AM) signal.

The acceleration signal shown in Fig. 9 is also used in the section. EMD is used to decompose the original signal and the result is presented in Fig. 13. The $x$ axis signifies the time and the $y$ represents amplitude in Fig. 14.

Fig. 13IMFs obtained by using EMD

Fig. 14The first six order IMFs of acceleration signal decomposed by EMD and their corresponding spectrums

It can be observed from Fig. 14 that the complicated acceleration vibration data is effectively decomposed into ten IMFs and residue, which reflect the information of different frequency bands. It should be noted from Fig. 14 that the first IMF1 contains the amplitude-modulated features, which is generated due to the periodic rub-impact occurred between the rotor and the casing. Furthermore, Fig. 14 also clearly shows that there exist the 1/2×, 1×, and some sup-× components etc. in the spectrum. These components are identified as the indicators of rub-impact in most cases, in other words, rubbing fault could be detected with EMD. The EMD approach successfully extracts all the mentioned signatures once more. This application, thus further demonstrates the validity of EMD in diagnosing rubbing fault.

In order to verify the effectiveness, some comparisons are performed. In the past, cepstral analysis has been applied in the rotating machinery condition monitoring since it can extract excitation information from the vibration signal [3, 7]. However, when the cepstral is employed to rotor vibration signals analysis, the impulse response function and the excitation should meet certain conditions, which the rotor vibration signals may not satisfy [27]. Wavelet owns the variable time-frequency window and could decompose a complicated signal into a number of components, which enable wavelet analysis to distinguish the rub-impact component from the rotor system vibration signal [6, 28]. However, wavelet analysis is necessarily an adjustable window Fourier transformation. Due to the limited length of wavelet base function, energy leakage will happen in wavelet transformation of signal. Furthermore, once wavelet base function and decomposition scales were confirmed, the result of wavelet transform would be the signal under a fixed scale, whose frequency components relate only to the sampling frequency of original signal and have no relation to the original signal itself. A demodulation method on the basis of the local mean decomposition (LMD) was proposed for rub-impact fault diagnosis [29], which can extract the transient fluctuations of the instantaneous frequency of the fundamental harmonic component. However, it should be noted that using LMD, the mode mixing problem cannot be solved.