# Research on Fault Feature Extraction Method Based on NOFRFs and Its Application in Rotor Faults

Authors: Yang Liu; Yulai Zhao; Jintao Li; Fangquan Xi; Shuanghe Yu; Ye Zhang
Journal: Shock and Vibration
 (2019)

---

## Abstract

Rub-impact between the rotating and static parts is a more common fault. The occurrence of faults is often accompanied by the generation of nonlinear phenomena. However, it is difficult to find out because the nonlinear characteristics are not obvious at the beginning of the fault. As a new frequency domain-based method, nonlinear output frequency response functions (NOFRFs) use the vibration response to extract the nonlinear characteristics of the system. This method has a better recognition rate for fault detection. Also, it has been applied in structural damages detection, but the high-order NOFRFs have the characteristics that the signals are weak and the features are difficult to extract. On this basis, the concept of the weighted contribution rate of the NOFRFs is proposed in this paper. The variable weighted coefficients with orders are used to amplify the influence of high-order NOFRFs on the nonlinearity of the system so as to extract its fault characteristics. The new indexRI is proposed based on Clenshaw–Curtis quadrature formula to eliminate the effect of artificially selected weighted coefficients on sensitivity. Especially in the early stage of the fault, the new index varies greatly with the deepening of the fault. Both simulation and experimental results verify the validity and practicability of the new index. The new index has certain guiding significance in the detection of mechanical system faults.


---

## Body

## 1. Introduction

In the design process of modern mechanical products, in order to meet the requirements of small size, the distance between the components is getting smaller and smaller. For mechanical products containing rotating parts, this increases the risk of rub-impact faults between the rotor and stator.Scholars have done a lot of research on the diagnosis of rub-impact faults in rotating machinery [1–8]. The physical model established by the finite element method or the lumped mass method combined with the time domain frequency domain and the axial trajectory is the earliest developed method of detecting rub-impact faults. This method mainly relies on the accuracy of the model [2, 3]. Since then, a large amount of research work has been done to make the model more accurate [4]. In the further development, the detection method of the rub-impact faults based on the Poincaré diagram, the bifurcation diagram, and the Lyapunov exponent is gradually applied, which has great influence on the development of the rub-impact faults detection [5, 6]. In recent years, with the continuous improvement of signal processing technology, some special signal processing methods have also been applied to the diagnosis of rub-impact faults [7–10]. Chandra and Sekhar [9] compared the effects of three kinds of signal processing methods, short time Fourier transform (STFT), continuous wavelet transform (CWT), and Hilbert–Huang transform (HHT), on the rub-impact faults diagnosis. Zhou et al. [10] studied the intrawave frequency modulation phenomenon and extracted the fault characteristics of the rub-impact from the vibration signal by using the M-VNCMD method. Most of the above studies only qualitatively describe the fault characteristics of the rub-impact. The feature extraction of rub-impact faults, especially the weak fault features at the beginning of the rub-impact, is not well developed.Researchers found that the frequency domain output of the system will have nonlinear phenomena [11, 12], such as harmonics, intermodulation, and so on when the mechanical system is damaged. However, the knowledge of the linear method cannot explain the generation of nonlinear phenomena [13]. Based on this, scholars began to focus on nonlinear methods to find out the mechanism of the system faults. They have developed a number of nonlinear methods for the diagnosis of system faults, such as the nonlinear decomposition and regulation method [14] and the second-order output spectrum approach [15, 16]. Among many nonlinear methods, the nonlinear method based on Volterra series has been accepted by many scholars. The Volterra series can explain the nonlinear components of the frequency domain output of system. Bai and Zhang [17] proposed a blind identification method based on third-order cumulant and an inverse recursive method and used this method to identify the simplified secondary Volterra system of hydraulic turbine shafting. Tang et al. [18] used the basic theory of Volterra series to model the rotating machinery and diagnose the rotating machinery. The results show that the mechanical system fault is identified by the change of the generalized frequency response function under different operating conditions of the mechanical system. Xia et al. [19] proposed a key kernels-PSO (KK-PSO) method that identifies the Volterra kernel and used the kernel function of the Volterra series to establish a rotor-bearing fault diagnosis system. The results indicate that the KK-PSO method has better recognition ability.However, the Volterra series is quite complex and can only be applied to some simple nonlinear problems. The generalized frequency response functions (GFRFs) is a concept based on the Volterra series. It can be regarded as a case where the frequency response function (FRF) of a linear system is generalized to a nonlinear system. Scholars [20–22] used GFRFs to study the nonlinear components of the frequency domain output of the system. Since GFRF is a multidimensional function, it is difficult to calculate it accurately in practice. Therefore, Lang and Billings [23] proposed the concept of nonlinear output frequency response functions (NOFRFs), which are a one-dimensional function that allows the calculation of the frequency domain output of a nonlinear system in a similar way to calculation of the frequency domain output of a linear system. Peng et al. [11] used NOFRFs to detect cracked beams and found an input-independent relationship between NOFRFs and GFRFs when the system is subjected to harmonic excitation. The experimental results show that the input energy of the system will be converted into superharmonic components when the excitation frequency ω is close to the resonance frequency of the associated NOFRFs. This also implies that NOFRFs can be used to analyze the behavior of nonlinear systems. Xia et al. [24] proposed a new method for online fault recognition of hydroelectric generators and studied the fault mechanism under different conditions. The results prove that this method is simple and efficient. Huang et al. [25] established the NARMAX model identification with PSO-adaptive lasso algorithm. A damage detection method combining the NARMAX model with a rectangular pulse is proposed. The results show that the method can effectively detect the accumulated fatigue damage of used parts.Under the development of nonlinear output frequency response functions, researchers have begun to focus on relevant indexes based on NOFRFs [26, 27]. Peng et al. [26] proposed a method based on NOFRFs related index Fe. The feasibility of the method was verified by the detection of damaged aluminum plates under the NARMAX model. Huang et al. [25] proposed a new index Ne. The NOFRFs of the nonlinear system are calculated by the NARMAX model, and then the index Ne is calculated. Comparing with the index of nondamaged test pieces, the results reveal that the newly proposed method can be used to detect fatigue accumulation damage of used parts. Mao et al. [27] proposed a nonlinear scale index NL and a divergence index NDI based on the basis of NOFRFs. The experimental results show that these indexes can well identify the morphology and extent of microcracks before they appear.In [25–27], the structural damage of test piece was tested at static state (nonworking state). The system excitation used in these references is a method of pulse excitation. However, this method is no longer applicable for online monitoring of rotating machinery faults. Therefore, the authors did the following research studies. It is pointed out in [26] that the NOFRFs of a nonlinear system can be calculated under the condition that the two amplitudes are different but the excitation frequencies are the same. Due to manufacturing defects, the disc is more or less imbalanced. This will cause the disc to generate centrifugal force during rotation. The existence of centrifugal force just makes the rotor system meet the conditions for calculating NOFRFs in [26] so that the NOFRFs of the system can be obtained. During the research, it was found that the indexes proposed in [25, 26] did not change significantly in the rotor rub-impact fault system. In order to find out a suitable index for online monitoring of rotor rub-impact faults, a new index based on Clenshaw–Curtis integration method [28] is used in this paper, which is the second-order weighted contribution rate integral index RI. The simulation and experimental results show that the index is more suitable for detecting the rotor rub-impact faults than these indexes in [25, 26] and has certain efficiency.The contribution of this paper is to combine the weighted contribution rate of NOFRFs with Clenshaw–Curtis integration method for the first time. A new system nonlinear detection indexRI is proposed, which is used to detect the rub-impact faults between rotor and stator. The organization of the rest paper is summarized as follows. In Section 2, we introduce the method of calculating NOFRFs under harmonic excitation. In Section 3, we combine the second-order optimal weighted contribution rate of NOFRFs with the Clenshaw–Curtis integration method and propose the index RI. After that, we use the simulation to verify the effectiveness of the index. In Section 4, the effectiveness of this index in the detection of rotor rub-impact fault is verified by experiments. Finally, the conclusion is given in Section 5.


## 2. Theory of NOFRFs

### 2.1. NOFRFs under Harmonic Input

Transfer functions are often used to describe the characteristics of linear systems. It can be obtained from the inputs and outputs of a linear system. Based on the similar view, GFRFs were seen as the extension of the transfer functions of the linear system to the nonlinear system situation. At this point, the output of the nonlinear system can be expressed as(1)Yjω=∑n=1NYnjω,where(2)Ynjω=1/n2πn−1×∫ω1+⋯+ωn=ωHnjω1,…,jωn∏i=1nUjωidσnω,where N is the highest order of the system nonlinearity, U(jω) and Y(jω) are the frequency domain expressions of the system time domain input u(t) and the output y(t) after Fourier transform, Yn(jω) represents the nth-order output response of the nonlinear system in the frequency domain, and Hn(jω) is the nth-order GFRFs. It can be seen from its expression that GFRF is a multidimensional function, and its dimensions are the same as orders. As the number of dimensions increase, the amount of calculations of GFRFs become very large, which limit the application of GFRFs.Based on this, the concept of NOFRFs is proposed.On the condition that(3)∫ω1+⋯+ωn=ω∏i=1nUjωidσnω≠0,its specific expression is defined as(4)Gnjω=∫ω1+⋯+ωn=ωHnjω1,…,jωn∏i=1nUjωidσnω∫ω1+⋯+ωn=ω∏i=1nUjωidσnω.Its specific solving method will be mentioned in [12]. On the basis of the system NOFRFs acquired, the frequency domain output of the system can be rewritten as(5)Yjω=∑n=1NYnjω=∑n=1NGnjωUnjω.This is similar to the form of the linear system transfer function, which decomposes the nonlinear system into a summation form of multiple linear subsystems, and the nonlinear characteristics of the system will be contained in the transfer function of each linear subsystem. This greatly reduces the amount of calculations and avoids dimensionality disasters. It is advantageous for the application of the method in the field of engineering as well as further theoretical research studies.When the system is excited by harmonics of the formu(t) = Asin (ωt), the system frequency domain output can be simplified [29] as(6)Yjkω=∑n=2N−k+1/2Gk+2n−1jkωUk+2n−1jkω,k=0,1,2,…,N.System frequency domain input is defined as(7)Unj−n+2kω=12nn!k!n−k!An.System NOFRFs can be simplified to(8)Gnj−n+2kω=Hnjω,…,jω︷k,−jω,…,−jω︷n−k.According to equations (6)–(8), under the condition that the first four orders of NOFRFs are sufficient to represent the nonlinear characteristics of the system, the relationship between the frequency domain output response and the input is(9)YjωF=G1jωFU1jωF+G3jωFU3jωF,Yj2ωF=G2j2ωFU2j2ωF+G4j2ωFU4j2ωF,Yj3ωF=G3j3ωFU3j3ωF,Yj4ωF=G4j4ωFU4j4ωF.According to the method mentioned in [12], when two harmonics with different intensities but the same frequency are used to excite the system, the first four orders of NOFRFs can be obtained, namely, G1(jω), G2(jω), G3(jω), G3(j3ω), G4(j2ω), and G4(j4ω). They can be used directly for fault detection, and they can be further extracted from fault features.


### 2.2. Feature Extraction Using Traditional Methods

#### 2.2.1. IndexFe

NOFRFs can be used directly to detect damages to mechanical system or structure. In addition, the nonlinear features in the NOFRFs can be further extracted to detect faults in the mechanical system.Fe [15] is the earliest proposed NOFRFs-based index. Its specific expression is as follows:(10)Fen=∫−∞+∞Gnjωdω∑i=1N∫−∞+∞Gijωdω,1≤n≤N.Equation (10) shows that the index Fe integrates all the NOFRFs values of the nonlinear system. It has the following characteristic:(11)∑n=1NFen=1.Fe(n) can reflect the ratio of the nth-order NOFRFs in the system. For example, Fe(1) reflects the proportion of G1(jω) in the system. If Fe(1) ≈ 1, it means that G1(jω) is dominant, and high-order system frequency response NOFRFs can be ignored; that is, the system can be approximated as a linear system. In addition, if Fe(5) ≈ 1, it indicates that the fifth-order NOFRFs dominate the system behaviors. Therefore, the index Fe(n) can be used to describe the nonlinear behaviors of the system and to detect damages to the system.


#### 2.2.2. IndexNe

The information entropyNe [16] fuses all orders of NOFRFs, and its specific expression is as follows:(12)NePY=−∑n=1NPYnlogPYn,where(13)PYn=∫−∞+∞Gnjω/ndω∑i=1N∫−∞+∞Gijω/idω,1≤n≤N.SinceNe integrates all the orders of NOFRFs, it becomes easier to apply Ne, which is more suitable for fault detection in engineering than Fe. Yet, Ne may be integrated into useless information. Assume that Fe(3) cannot be used to detect system faults, but Ne also fuses Fe(3), which will reduce the sensitivity of the index to system faults detection. Therefore, the authors attempted to identify a more sensitive index to detect system faults and damages.


### 2.3. Weighted Contribution Rate of NOFRFs (Rn)

During the operation of the mechanical system, the system is prone to fault due to abominable working conditions and alternating loads. At this moment, the system will express nonlinearity and the nonlinearity will become apparent as the fault increase. When the fault reaches a certain level, it will lead to an accident. Therefore, in engineering practice, online monitoring of mechanical systems is of great significance. NOFRFs and related indexes can be used to detect system faults. In the actual solution process, high-order NOFRFs become smaller and smaller as the order increases, which is difficult to extract. High-order NOFRFs have a significant effect on the detection of system faults. Therefore, it is necessary to extract the nonlinear features of high-order NOFRFs. Based on this, a new method is proposed in this paper, which is to weight the high-order NOFRFs and increase their contribution rate to the system. The specific weighted process is as follows:(14)Tnjω=Gnjωnρ,where Tn(jω) represents Gn(jω) after weighted, nρ represents a weighted coefficient, and ρ is an indefinite constant. It can be seen that the weighted method adopted in this paper is the variable weighted method. When ρ < 0, the higher the order of NOFRFs and the greater the weighting, as follows:(15)1=T1jωG1jω<T2jωG2jω<⋯<TnjωGnjω.It will increase the contribution of high-order NOFRFs to the system, thus achieving the purpose of feature extraction. These are the basic theory of the weighted contribution rate of the NOFRFs, which is expressed as(16)Rnn=∫−∞+∞Tnjωdω∑i=1N∫−∞+∞Tijωdω=∫−∞+∞Gnjω/nρdω∑i=1N∫−∞+∞Gijω/iρdω,1≤n≤N,ρ∈−∞,0.The high-orderGnjωn>1 nonlinear characteristics of the system will be amplified by the weighted contribution rate of NOFRFs using feature weighted method. The weighted contribution rate index Rn introduces the order n of the NOFRFs, which causes the weighted coefficients to increase as the order of the NOFRFs increases. The contribution rate of the high-order Gnjωn>1 will be magnified under the condition ρ<0. Only in this way, the weighted coefficients at the denominator position is less than 1, namely, nρ<1. Therefore, the weighted contribution rate Rn solves the problem that the high-order NOFRFs and their contribution are too small.It can be seen from equation (16) that when the indefinite constant ρ = 0, the weighted contribution index Rn(1) is equal to Fe(1), so the feature of Rn(1) is the same as Fe(1); when the indefinite constant ρ = 0.5, Rn is equal to PY, and PY will be replaced the index Rn. From the above analysis, it can be discovered that the weighted contribution rate index Rn is a generalized form of the two indexes of Fe and PY. There is significant effect on the detection index when different indefinite constants were chosen. Therefore, how to choose an optimal indefinite constant or even how to remove its effect is very important for the accuracy of the detection index. If the final index is not affected by the indefinite constant selection, this will greatly improve the accuracy of the detection index and can be more conducive to the application of the NOFRFs in practical engineering.


## 2.1. NOFRFs under Harmonic Input

Transfer functions are often used to describe the characteristics of linear systems. It can be obtained from the inputs and outputs of a linear system. Based on the similar view, GFRFs were seen as the extension of the transfer functions of the linear system to the nonlinear system situation. At this point, the output of the nonlinear system can be expressed as(1)Yjω=∑n=1NYnjω,where(2)Ynjω=1/n2πn−1×∫ω1+⋯+ωn=ωHnjω1,…,jωn∏i=1nUjωidσnω,where N is the highest order of the system nonlinearity, U(jω) and Y(jω) are the frequency domain expressions of the system time domain input u(t) and the output y(t) after Fourier transform, Yn(jω) represents the nth-order output response of the nonlinear system in the frequency domain, and Hn(jω) is the nth-order GFRFs. It can be seen from its expression that GFRF is a multidimensional function, and its dimensions are the same as orders. As the number of dimensions increase, the amount of calculations of GFRFs become very large, which limit the application of GFRFs.Based on this, the concept of NOFRFs is proposed.On the condition that(3)∫ω1+⋯+ωn=ω∏i=1nUjωidσnω≠0,its specific expression is defined as(4)Gnjω=∫ω1+⋯+ωn=ωHnjω1,…,jωn∏i=1nUjωidσnω∫ω1+⋯+ωn=ω∏i=1nUjωidσnω.Its specific solving method will be mentioned in [12]. On the basis of the system NOFRFs acquired, the frequency domain output of the system can be rewritten as(5)Yjω=∑n=1NYnjω=∑n=1NGnjωUnjω.This is similar to the form of the linear system transfer function, which decomposes the nonlinear system into a summation form of multiple linear subsystems, and the nonlinear characteristics of the system will be contained in the transfer function of each linear subsystem. This greatly reduces the amount of calculations and avoids dimensionality disasters. It is advantageous for the application of the method in the field of engineering as well as further theoretical research studies.When the system is excited by harmonics of the formu(t) = Asin (ωt), the system frequency domain output can be simplified [29] as(6)Yjkω=∑n=2N−k+1/2Gk+2n−1jkωUk+2n−1jkω,k=0,1,2,…,N.System frequency domain input is defined as(7)Unj−n+2kω=12nn!k!n−k!An.System NOFRFs can be simplified to(8)Gnj−n+2kω=Hnjω,…,jω︷k,−jω,…,−jω︷n−k.According to equations (6)–(8), under the condition that the first four orders of NOFRFs are sufficient to represent the nonlinear characteristics of the system, the relationship between the frequency domain output response and the input is(9)YjωF=G1jωFU1jωF+G3jωFU3jωF,Yj2ωF=G2j2ωFU2j2ωF+G4j2ωFU4j2ωF,Yj3ωF=G3j3ωFU3j3ωF,Yj4ωF=G4j4ωFU4j4ωF.According to the method mentioned in [12], when two harmonics with different intensities but the same frequency are used to excite the system, the first four orders of NOFRFs can be obtained, namely, G1(jω), G2(jω), G3(jω), G3(j3ω), G4(j2ω), and G4(j4ω). They can be used directly for fault detection, and they can be further extracted from fault features.


## 2.2. Feature Extraction Using Traditional Methods

### 2.2.1. IndexFe

NOFRFs can be used directly to detect damages to mechanical system or structure. In addition, the nonlinear features in the NOFRFs can be further extracted to detect faults in the mechanical system.Fe [15] is the earliest proposed NOFRFs-based index. Its specific expression is as follows:(10)Fen=∫−∞+∞Gnjωdω∑i=1N∫−∞+∞Gijωdω,1≤n≤N.Equation (10) shows that the index Fe integrates all the NOFRFs values of the nonlinear system. It has the following characteristic:(11)∑n=1NFen=1.Fe(n) can reflect the ratio of the nth-order NOFRFs in the system. For example, Fe(1) reflects the proportion of G1(jω) in the system. If Fe(1) ≈ 1, it means that G1(jω) is dominant, and high-order system frequency response NOFRFs can be ignored; that is, the system can be approximated as a linear system. In addition, if Fe(5) ≈ 1, it indicates that the fifth-order NOFRFs dominate the system behaviors. Therefore, the index Fe(n) can be used to describe the nonlinear behaviors of the system and to detect damages to the system.


### 2.2.2. IndexNe

The information entropyNe [16] fuses all orders of NOFRFs, and its specific expression is as follows:(12)NePY=−∑n=1NPYnlogPYn,where(13)PYn=∫−∞+∞Gnjω/ndω∑i=1N∫−∞+∞Gijω/idω,1≤n≤N.SinceNe integrates all the orders of NOFRFs, it becomes easier to apply Ne, which is more suitable for fault detection in engineering than Fe. Yet, Ne may be integrated into useless information. Assume that Fe(3) cannot be used to detect system faults, but Ne also fuses Fe(3), which will reduce the sensitivity of the index to system faults detection. Therefore, the authors attempted to identify a more sensitive index to detect system faults and damages.


## 2.2.1. IndexFe

NOFRFs can be used directly to detect damages to mechanical system or structure. In addition, the nonlinear features in the NOFRFs can be further extracted to detect faults in the mechanical system.Fe [15] is the earliest proposed NOFRFs-based index. Its specific expression is as follows:(10)Fen=∫−∞+∞Gnjωdω∑i=1N∫−∞+∞Gijωdω,1≤n≤N.Equation (10) shows that the index Fe integrates all the NOFRFs values of the nonlinear system. It has the following characteristic:(11)∑n=1NFen=1.Fe(n) can reflect the ratio of the nth-order NOFRFs in the system. For example, Fe(1) reflects the proportion of G1(jω) in the system. If Fe(1) ≈ 1, it means that G1(jω) is dominant, and high-order system frequency response NOFRFs can be ignored; that is, the system can be approximated as a linear system. In addition, if Fe(5) ≈ 1, it indicates that the fifth-order NOFRFs dominate the system behaviors. Therefore, the index Fe(n) can be used to describe the nonlinear behaviors of the system and to detect damages to the system.


## 2.2.2. IndexNe

The information entropyNe [16] fuses all orders of NOFRFs, and its specific expression is as follows:(12)NePY=−∑n=1NPYnlogPYn,where(13)PYn=∫−∞+∞Gnjω/ndω∑i=1N∫−∞+∞Gijω/idω,1≤n≤N.SinceNe integrates all the orders of NOFRFs, it becomes easier to apply Ne, which is more suitable for fault detection in engineering than Fe. Yet, Ne may be integrated into useless information. Assume that Fe(3) cannot be used to detect system faults, but Ne also fuses Fe(3), which will reduce the sensitivity of the index to system faults detection. Therefore, the authors attempted to identify a more sensitive index to detect system faults and damages.


## 2.3. Weighted Contribution Rate of NOFRFs (Rn)

During the operation of the mechanical system, the system is prone to fault due to abominable working conditions and alternating loads. At this moment, the system will express nonlinearity and the nonlinearity will become apparent as the fault increase. When the fault reaches a certain level, it will lead to an accident. Therefore, in engineering practice, online monitoring of mechanical systems is of great significance. NOFRFs and related indexes can be used to detect system faults. In the actual solution process, high-order NOFRFs become smaller and smaller as the order increases, which is difficult to extract. High-order NOFRFs have a significant effect on the detection of system faults. Therefore, it is necessary to extract the nonlinear features of high-order NOFRFs. Based on this, a new method is proposed in this paper, which is to weight the high-order NOFRFs and increase their contribution rate to the system. The specific weighted process is as follows:(14)Tnjω=Gnjωnρ,where Tn(jω) represents Gn(jω) after weighted, nρ represents a weighted coefficient, and ρ is an indefinite constant. It can be seen that the weighted method adopted in this paper is the variable weighted method. When ρ < 0, the higher the order of NOFRFs and the greater the weighting, as follows:(15)1=T1jωG1jω<T2jωG2jω<⋯<TnjωGnjω.It will increase the contribution of high-order NOFRFs to the system, thus achieving the purpose of feature extraction. These are the basic theory of the weighted contribution rate of the NOFRFs, which is expressed as(16)Rnn=∫−∞+∞Tnjωdω∑i=1N∫−∞+∞Tijωdω=∫−∞+∞Gnjω/nρdω∑i=1N∫−∞+∞Gijω/iρdω,1≤n≤N,ρ∈−∞,0.The high-orderGnjωn>1 nonlinear characteristics of the system will be amplified by the weighted contribution rate of NOFRFs using feature weighted method. The weighted contribution rate index Rn introduces the order n of the NOFRFs, which causes the weighted coefficients to increase as the order of the NOFRFs increases. The contribution rate of the high-order Gnjωn>1 will be magnified under the condition ρ<0. Only in this way, the weighted coefficients at the denominator position is less than 1, namely, nρ<1. Therefore, the weighted contribution rate Rn solves the problem that the high-order NOFRFs and their contribution are too small.It can be seen from equation (16) that when the indefinite constant ρ = 0, the weighted contribution index Rn(1) is equal to Fe(1), so the feature of Rn(1) is the same as Fe(1); when the indefinite constant ρ = 0.5, Rn is equal to PY, and PY will be replaced the index Rn. From the above analysis, it can be discovered that the weighted contribution rate index Rn is a generalized form of the two indexes of Fe and PY. There is significant effect on the detection index when different indefinite constants were chosen. Therefore, how to choose an optimal indefinite constant or even how to remove its effect is very important for the accuracy of the detection index. If the final index is not affected by the indefinite constant selection, this will greatly improve the accuracy of the detection index and can be more conducive to the application of the NOFRFs in practical engineering.


## 3. Second-Order Weighted Contribution Rate Integral Index

### 3.1. Estimation of NOFRFs of Rotor-Stator Rub-Impact System

In order to eliminate the influence of the indefinite constants selection on the detection index, we first look at a practical example to realize how to affect the weighted contribution rate of the NOFRFs by the selection of the indefinite constants. The Jeffcott rotor fixed-point rub-impact system will be used as an example in this paper. The schematic diagram of the Jeffcott rotor fixed-point rub-impact system is shown in Figure1. The rotor system is divided into 11 shaft segments, and the specific parameters of each shaft segment are shown in Table 1. Each shaft segment has two nodes, so the rotor system has a total of 12 nodes.Figure 1
Jeffcott rotor fixed-point rub-impact system schematic.Table 1
Specific unit parameters of segments.
Segment no.
1
2
3
4
5
6
7
8
9
10
11
Length (mm)
50
45
45
45
45
27
45
45
45
45
50
Diameter (mm)
10
10
10
10
10
70
10
10
10
10
10Each shaft segment of the rotor system is simulated with Timoshenko beam, and the generalized coordinates of the beam element isu = (xA, yA, θA, θA, xB, yB, θB, θB)T. A and B represent the endpoints of the beam element, respectively, where xA, yA, θA, θA are the displacement and rotation angle of the x and y directions of the A node, respectively, and xB, yB, θB, θB are the displacements and rotation angles of the x and y directions of the B node, respectively.The fixed-point rub-impact of rotor-stator is mainly simulated in this paper; in other words, rub-impact will appear between a certain point on the rotating shaft and a stator. The friction torque generated by the rub-impact is ignored, and only the normal rub-impact forceFn and the tangential rub-impact force Fr are considered when the phenomenon of rub-impact appears in rotor-stator, as shown in Figure 2. In Figure 2, O1 represents the centroid of the rotor when it is static, O1′ represents the centroid of the rotor when it is dynamic, O2 denotes the center of mass of the rotor when it is static, δ0 represents the initial clearance between rotor and stator, and e represents the eccentricity. The friction coefficient and the contact stiffness between rotor and stator are denoted by kc and fc, respectively.Figure 2
Mechanics model of fixed-point rub-impact.In this paper, the Kuhn–Tucker law of Coulomb friction is used to detect the contact between rotor and stator, and the contact condition is solved by the augmented Lagrangian method of friction [30]. These ensure that the disc and the elastic rod meet the contact conditions and do not penetrate each other. In order to meet the penetration depth is within the allowable range, a suitable force will be applied to the normal contact in this method. The rub-impact force Frub between the rotor and the stator can be decomposed into the normal rub-impact force Fn and the tangential rub-impact force Fr, as follows:(17)Fn=0,r<δ0,r−δ0kc+λri+1,r≥δ0,Fr=fcFn,where r is the distance between O1 and O2 and λri+1 is the Lagrangian multiplier of the i + 1th step in the iterative process, which is expressed as(18)λri+1=λri+kcδr,δr≥ς,λri,δr<ς,where δr = r − c, ς is prescriptive permeability-tolerance. When ς exceeds the maximum allowable permeation amount after a given equilibrium iteration, the Lagrangian multiplier is used to enlarge the contact stiffness of the contact units. This process does not end until the magnitude of penetration is less than the maximum allowable value.Then, the rub-impact force can be expressed as(19)FrubxFruby=−Hδrδrδrkc+λri+1r1−fcfc1xy,where the Heaviside function H(δr) is defined as(20)Hδr=1,δr≥0,0,δr<0.Therefore, the governing equation of the rotor rub-impact system is expressed as(21)Mq¨+C+Gq˙+Kq=Fe+Frub,where, M, C, G, and K are the mass matrix, damping matrix, gyroscopic matrix, and stiffness matrix, respectively; Fe is the imbalance force vector; and Frub is the rub-impact force vector of the rotor system [31].The rest of simulation parameters are shown in Table2. The Newmark-β and Newton–Raphson iterative methods are used to solve equation (21), and the vibration response of each node at the excitation frequency of 20 Hz can be obtained. Then, two vibration responses with different excitation strengths but the same excitation frequency can be obtained only by changing the imbalance and solving again. According to equation (9) and [12], the first four orders of NOFRFs of the system can be obtained by the two inputs and outputs of the system. The specific values are shown in Table 3.Table 2
Simulation parameters of the rotor system.
Parameters
Values
Material
45#steel
Elastic modulus of shaft Ef (Pa)
2.1 × 1011
Density (kg/m3)
7850
Poisson ratioμv
0.3
Normal rub-impact stiffness ratio
2/4/6 × 104
Friction coefficient
0.1/0.2/0.3
Bearing stiffness (N/m)
2 × 106
Imbalance (kg·m)
167.4/236.7 × 10−6
Rotational speed (rpm)
1200Table 3
The NOFRFs of rub-impact rotor system.
NOFRFs
Normal
Clearance = 90μm
Clearance = 30μm
G
1
j
ω
1.5078 × 10−4
1.5069 × 10−4
1.5051 × 10−4
G
2
j
ω
2.1286 × 10−8
4.8293 × 10−7
1.2208 × 10−6
G
3
j
ω
2.5299 × 10−9
1.1921 × 10−8
3.9750 × 10−9
G
3
j
3
ω
8.9336 × 10−10
5.7439 × 10−8
1.0286 × 10−7
G
4
j
2
ω
2.7172 × 10−9
2.2947 × 10−9
3.5165 × 10−8
G
4
j
4
ω
3.8508 × 10−10
1.7238 × 10−8
3.0804 × 10−8


### 3.2. Feature Extraction Based on Weighted Contribution Rate of NOFRFs

According to the analysis in Section2.3, the contribution rate of the high-order NOFRFs Gnjωn>1 will be increased by the weighted contribution index Rn of the NOFRFs to the system. However, the choice of indefinite constants can have a large effect on the results. The authors attempt to eliminate these effects. For this result, the authors research the effect of indefinite constants on the weighted contribution rate. First, the obtained NOFRFs values are taken into equation (16) to calculate Rn, the range of indefinite constant ρ is [−30, 5], and a curve of Rn respect to ρ is then made under three different rub-impact clearances (normal, clearance A = 90 μm, and clearance B = 30 μm).When the indefinite constantρ changes within the interval [−30, 5], the trend of the weighted contribution index Rn is as shown in Figure 3. It can be seen from the trend in Figure 3 that Rn is almost equal and tends to 0 under different rub-impact clearances when the indefinite constant range is (−∞, −30] or [5, +∞), which means that the weighted contribution rate Rn will fail to diagnose unceasingly the severity of rub-impact of the system when ρ belongs to at this range.Figure 3
Rn change with indefinite constant (A = 90 μm, B = 30 μm). (a) Rn(1). (b) Rn(2). (c) Rn(3). (d) Rn(4).
(a)
(b)
(c)
(d)In addition, Figure3(a) shows that the two indexes Fe1(ρ = 0) and PY1(ρ = 0.5) increase with the decrease of the rub-impact clearance, but the weighted contribution rate increases from 0 to 1 in the process of the indefinite constant increasing from negative infinity to positive infinity; from Figures 3(a)–3(d), it can be seen that Fe(n) and PY decrease when the rub-impact clearance decreases. The value of Fe(n) and PY(n) n=2,3,4 is also increasing. As seen in Figure 3, the sensitivity of the index Fe is higher than that of PY, but from this, we can also know that in the indefinite constant ρ = 0 or ρ = 0.5, Fe and PY are not the best choice. As shown in Figure 3(b), the peak-weighted contribution rate sensitivity is the highest in the case of different rub-impact clearances, but the indefinite constant corresponding to the peaks changes as the rub-impact clearances change. It can be seen from Figure 3(b) that the trend of the weighted contribution rate Rn(2) is increased first and then decreased, and the curve and the horizontal axis constitute a closed graph when the indefinite constant change. A new index is proposed by using the integral of the indefinite constant ρ of the weighted contribution rate Rn(2). Geometrically, this new index corresponds to the area of the weighted contribution rate Rn(2) and the horizontal axis. This index can better detect the change of rub-impact. Based on this idea, the index will be theoretically derived.The specific expression ofRn(2) can be expressed from equation (16), as follows:(22)Rn2=∫−∞+∞T2jωdω∑i=1N∫−∞+∞Tijωdω=∫−∞+∞G2jω/2ρdω∑i=1N∫−∞+∞Gijω/iρdω,ρ∈−∞,0.Defineψi as ψi=∫−∞+∞Gijωdω, and equation (22) can be written as(23)Rn2=ψ2/2ρ∑i=1Nψi/iρ,ρ∈−∞,0.Equation (23) implies that the index Rn(2) is a function of the indefinite constant ρ, and R2(ρ) is regarded as(24)R2ρ=ψ2/2ρ∑i=1Nψi/iρ,ρ∈−∞,0.According to the analysis in Section2.1, only the first four order NOFRFs of the system are considered here, so equation (24) can be expressed as(25)R2ρ=12ρψ1/ψ2+2/3ρψ3/ψ2+1/2ρψ4/ψ2+1,ρ∈−∞,0.In the practical engineering case,ψi/ψ2, i = 1, 3, 4 in the function R2(ρ) are constants. It can be seen from equation (25) that the function R2(ρ) gradually approaches zero when ρ tends to infinity. According to the analysis mentioned above, the function R2(ρ) is bounded and continuous. So, the function R2(ρ) can be integrated.


### 3.3. Clenshaw–Curtis Integration

In engineering practice, researchers often use numerical integration methods when they solve the integrations of unknown analytic functions. There are many methods for numerical integration. A method with higher precision will be adopted in this paper, namely, Clenshaw–Curtis quadrature formula [28]. Using this method to numerically integrate the index Rn(2), the basic theory is that for the integrand function f(x), its definite integral expression on the interval [−1, 1] is(26)If=∫−11fxdx.The Clenshaw–Curtis quadrature formula is defined as(27)Hnf=∑k=0nθ•kfζk,where ∑• indicates that the first and last items in the formula are halved. ζk=coskπ/n, θk=4/n∑j=0n1/1−j2•cosjkπ/n, k = 0, 1, 2, ..., n, j = 1, 2, 4, ..., n, ζk, θk are the quadrature node and the quadrature coefficient, respectively. Regarding the accuracy of Hn(f), the following theorem is given.Theorem 1.
Iff satisfies the condition of f∈Cr−1,1,r>1, the error is expressed as(28)enf=Hnf−If=on1−r.On a given finite integration interval [a, b], there is an integral of f(x),(29)If=∫abfxdx,where a, b are noninfinite constants and f(x) has a certain smooth degree. The following processing is often done before using Clenshaw–Curtis integration for numerical integration:(30)x=a+b2+b−a2t,that is, the integral interval is converted to [−1, 1] by variable substitution.However, in practical engineering problems, the researchers found that the results obtained did not achieve the desired by linear transformation in most cases. The reason for this is that the trend of the function at the endpoints is different. Considering a class of function,f(x) has a significant change in the vicinity of a in the interval [a, b], and the change in the vicinity of b is slow. Such as f(x) = e−x, when a = 1 and b = 50, the ideal numerical integral is to take more integral nodes near a, and at the same time, in order to reduce the amount of calculation, there are fewer integral nodes near b. Considering the Clenshaw–Curtis numerical integration in the interval [−1, 1], the integral nodes are symmetric about the origin, and there are many nodes at the ±1 and fewer nodes near the origin. With Clenshaw–Curtis numerical integration, linear transformations cannot meet our requirements. Therefore, in this case, a new method of variable substitution will be considered, namely, nonlinear variable substitution.Extending the linear transformation to nonlinearity, equation (30) is rewritten as(31)x=git=b−a2it+1i+a,−1≤t≤1,where i ≥ 1 and i ∈ Z. It can be seen from equation (31) that git can convert the interval [a, b] to [−1, 1], and then perform Clenshaw–Curtis numerical integration on f(x). When i = 1, git is a linear variable substitution; when i = 2, 3, 4, …, git is a nonlinear variable replacement. Under nonlinear transformation, rewrite equation (29) to(32)If=∫−11fgitgi′tdt=b−ai2i∫−11fb−a2it+1i+at+1i−1dt.As can be seen from equation (32), hit=fb−a/2it+1i+at+1i−1 and f(x) have the same smooth order. Then, for equation (32), the Clenshaw–Curtis quadrature formula is used to obtain(33)Inif=b−ai2i∑k=0nf•b−a2iζk+1i+aθkζk+1i−1.In the case of nonlinear variable substitution, the quadrature nodeξi=b−a/2iζk+1i+a,k=0,1,2,...,n is denser near the end point a and sparse at the end point b. This is what we expected.


### 3.4. Second-Order Weighted Contribution Rate Integral Index (RI)

According to the analysis in Section3.2, the function R2(ρ) is integral. The integral expression is(34)IR=∫abR2ρdρ.From Figure3(a) in Section 3.2, it can be concluded that the value of Rn(2) is close to 0 when the indefinite constant is within (−∞, −30] or [5, +∞). So, the integral interval of equation (34) is taken [−30, 5], that is, a = −30, b = 5. Substituting equation (25) into equation (33), equation (34) is rewritten as(35)IniR=35i2i∑k=0nR2•352iζk+1i−30θkζk+1i−1.DefineRI=IniR.In order to verify the reliability of the Clenshaw–Curtis quadrature formula, the authors use the int infinite integral function in MATLAB to solveR2(ρ), and the integration interval is negative infinity to positive infinity. The results are compared with the Clenshaw–Curtis integration method used in this paper, as shown in Table 4. The results show that the selection of the integral interval in the Clenshaw–Curtis method is reasonable, and the results obtained by the Clenshaw–Curtis quadrature method are very accurate.Table 4
Comparison of Clenshaw–Curtis integral and int integral function results.
Clearance (μm)
Clenshaw–Curtis integral
Int infinite integral
I
Clenshaw
−
Curtis
−
I
int
Normal
0.0670145365
0.0670145578
2.1289 × 10−8
100
0.3580265528
0.3580266286
7.5884 × 10−8
80
0.5712722586
0.5712723441
8.5445 × 10−8
60
0.6289452527
0.6289452534
7.0059 × 10−8
40
0.6855516020
0.6855515858
1.6181 × 10−8
20
0.7381875132
0.7381875624
4.9269 × 10−8Rub-impact stiffness and friction coefficient are the two main parameters affecting rub-impact. The change of these two parameters will affect the change of detection indexRI. In order to study the influence of these two parameters on the results, the authors studied the variation of the second-order weighted contribution rate integral index RI with the rub-impact clearances when these two parameters change, as shown in Figures 4 and 5. They illustrated the trend of the index RI of the system under different rub-impact stiffness and friction coefficient. The results show that the index RI will increase with the decrease of the rub-impact clearance when the rub-impact stiffness and friction coefficient of the system increase. The variation of the index RI will gradually increase, which means that the nonlinearity of the system will be increased.Figure 4
The trend ofRI under different rub-impact clearances and rub-impact stiffnesses.Figure 5
The trend ofRI under different rub-impact clearances and friction coefficients.In order to verify that the second-order weighted contribution rate integral indexRI is more suitable for detecting the rub-impact faults of the rotor, according to equation (33), the index RI under different rub-impact clearances is calculated. The amount of change of RI, Ne, Fe1, Fe2, Fe3, and Fe4 under different rub-impact clearance is shown in Table 5. From this, we can see that the indexes RI, Ne, and Fe1 have the largest change with the reduction of the rub-impact clearance. In order to more intuitively show the advantages of the index RI, we compare the three indexes in the form of a histogram, as shown in Figure 6. The results show that the indices RI and Ne increase with the increasing rub-impact clearance, but the index Fe1 decreases as the rub-impact clearance decreases. However, the change of the second-order weighted contribution rate integral index RI is more obvious. Especially in the initial stage of the rub-impact fault, this index has strong superiority. These indicate that the second-order weighted contribution rate integral index RI proposed in this paper is more suitable for detecting the rub-impact fault of the rotator-stator.Table 5
The amount of change ofRI, Ne, Fe1, Fe2, Fe3, and Fe4 under different rub-impact clearances (μm).
Indexes
100
80
60
40
20
ΔRI
0.2910
0.2132
0.0577
0.0566
0.0526
ΔNe
0.0142
7.388 × 10−3
7.413 × 10−3
7.041 × 10−3
6.754 × 10−3
ΔFe1
2.771 × 10−3
1.700 × 10−3
1.810 × 10−3
1.807 × 10−3
1.803 × 10−3
ΔFe2
2.242 × 10−3
1.619 × 10−3
1.616 × 10−3
1.613 × 10−3
1.610 × 10−3
ΔFe3
3.949 × 10−4
8.183 × 10−5
8.164 × 10−5
8.150 × 10−5
8.135 × 10−5
ΔFe4
1.349 × 10−4
6.657 × 10−7
1.119 × 10−4
1.117 × 10−4
1.116 × 10−4Figure 6
Comparison of simulation results of indexesRI, Fe1, and Ne under different clearances.


## 3.1. Estimation of NOFRFs of Rotor-Stator Rub-Impact System

In order to eliminate the influence of the indefinite constants selection on the detection index, we first look at a practical example to realize how to affect the weighted contribution rate of the NOFRFs by the selection of the indefinite constants. The Jeffcott rotor fixed-point rub-impact system will be used as an example in this paper. The schematic diagram of the Jeffcott rotor fixed-point rub-impact system is shown in Figure1. The rotor system is divided into 11 shaft segments, and the specific parameters of each shaft segment are shown in Table 1. Each shaft segment has two nodes, so the rotor system has a total of 12 nodes.Figure 1
Jeffcott rotor fixed-point rub-impact system schematic.Table 1
Specific unit parameters of segments.
Segment no.
1
2
3
4
5
6
7
8
9
10
11
Length (mm)
50
45
45
45
45
27
45
45
45
45
50
Diameter (mm)
10
10
10
10
10
70
10
10
10
10
10Each shaft segment of the rotor system is simulated with Timoshenko beam, and the generalized coordinates of the beam element isu = (xA, yA, θA, θA, xB, yB, θB, θB)T. A and B represent the endpoints of the beam element, respectively, where xA, yA, θA, θA are the displacement and rotation angle of the x and y directions of the A node, respectively, and xB, yB, θB, θB are the displacements and rotation angles of the x and y directions of the B node, respectively.The fixed-point rub-impact of rotor-stator is mainly simulated in this paper; in other words, rub-impact will appear between a certain point on the rotating shaft and a stator. The friction torque generated by the rub-impact is ignored, and only the normal rub-impact forceFn and the tangential rub-impact force Fr are considered when the phenomenon of rub-impact appears in rotor-stator, as shown in Figure 2. In Figure 2, O1 represents the centroid of the rotor when it is static, O1′ represents the centroid of the rotor when it is dynamic, O2 denotes the center of mass of the rotor when it is static, δ0 represents the initial clearance between rotor and stator, and e represents the eccentricity. The friction coefficient and the contact stiffness between rotor and stator are denoted by kc and fc, respectively.Figure 2
Mechanics model of fixed-point rub-impact.In this paper, the Kuhn–Tucker law of Coulomb friction is used to detect the contact between rotor and stator, and the contact condition is solved by the augmented Lagrangian method of friction [30]. These ensure that the disc and the elastic rod meet the contact conditions and do not penetrate each other. In order to meet the penetration depth is within the allowable range, a suitable force will be applied to the normal contact in this method. The rub-impact force Frub between the rotor and the stator can be decomposed into the normal rub-impact force Fn and the tangential rub-impact force Fr, as follows:(17)Fn=0,r<δ0,r−δ0kc+λri+1,r≥δ0,Fr=fcFn,where r is the distance between O1 and O2 and λri+1 is the Lagrangian multiplier of the i + 1th step in the iterative process, which is expressed as(18)λri+1=λri+kcδr,δr≥ς,λri,δr<ς,where δr = r − c, ς is prescriptive permeability-tolerance. When ς exceeds the maximum allowable permeation amount after a given equilibrium iteration, the Lagrangian multiplier is used to enlarge the contact stiffness of the contact units. This process does not end until the magnitude of penetration is less than the maximum allowable value.Then, the rub-impact force can be expressed as(19)FrubxFruby=−Hδrδrδrkc+λri+1r1−fcfc1xy,where the Heaviside function H(δr) is defined as(20)Hδr=1,δr≥0,0,δr<0.Therefore, the governing equation of the rotor rub-impact system is expressed as(21)Mq¨+C+Gq˙+Kq=Fe+Frub,where, M, C, G, and K are the mass matrix, damping matrix, gyroscopic matrix, and stiffness matrix, respectively; Fe is the imbalance force vector; and Frub is the rub-impact force vector of the rotor system [31].The rest of simulation parameters are shown in Table2. The Newmark-β and Newton–Raphson iterative methods are used to solve equation (21), and the vibration response of each node at the excitation frequency of 20 Hz can be obtained. Then, two vibration responses with different excitation strengths but the same excitation frequency can be obtained only by changing the imbalance and solving again. According to equation (9) and [12], the first four orders of NOFRFs of the system can be obtained by the two inputs and outputs of the system. The specific values are shown in Table 3.Table 2
Simulation parameters of the rotor system.
Parameters
Values
Material
45#steel
Elastic modulus of shaft Ef (Pa)
2.1 × 1011
Density (kg/m3)
7850
Poisson ratioμv
0.3
Normal rub-impact stiffness ratio
2/4/6 × 104
Friction coefficient
0.1/0.2/0.3
Bearing stiffness (N/m)
2 × 106
Imbalance (kg·m)
167.4/236.7 × 10−6
Rotational speed (rpm)
1200Table 3
The NOFRFs of rub-impact rotor system.
NOFRFs
Normal
Clearance = 90μm
Clearance = 30μm
G
1
j
ω
1.5078 × 10−4
1.5069 × 10−4
1.5051 × 10−4
G
2
j
ω
2.1286 × 10−8
4.8293 × 10−7
1.2208 × 10−6
G
3
j
ω
2.5299 × 10−9
1.1921 × 10−8
3.9750 × 10−9
G
3
j
3
ω
8.9336 × 10−10
5.7439 × 10−8
1.0286 × 10−7
G
4
j
2
ω
2.7172 × 10−9
2.2947 × 10−9
3.5165 × 10−8
G
4
j
4
ω
3.8508 × 10−10
1.7238 × 10−8
3.0804 × 10−8


## 3.2. Feature Extraction Based on Weighted Contribution Rate of NOFRFs

According to the analysis in Section2.3, the contribution rate of the high-order NOFRFs Gnjωn>1 will be increased by the weighted contribution index Rn of the NOFRFs to the system. However, the choice of indefinite constants can have a large effect on the results. The authors attempt to eliminate these effects. For this result, the authors research the effect of indefinite constants on the weighted contribution rate. First, the obtained NOFRFs values are taken into equation (16) to calculate Rn, the range of indefinite constant ρ is [−30, 5], and a curve of Rn respect to ρ is then made under three different rub-impact clearances (normal, clearance A = 90 μm, and clearance B = 30 μm).When the indefinite constantρ changes within the interval [−30, 5], the trend of the weighted contribution index Rn is as shown in Figure 3. It can be seen from the trend in Figure 3 that Rn is almost equal and tends to 0 under different rub-impact clearances when the indefinite constant range is (−∞, −30] or [5, +∞), which means that the weighted contribution rate Rn will fail to diagnose unceasingly the severity of rub-impact of the system when ρ belongs to at this range.Figure 3
Rn change with indefinite constant (A = 90 μm, B = 30 μm). (a) Rn(1). (b) Rn(2). (c) Rn(3). (d) Rn(4).
(a)
(b)
(c)
(d)In addition, Figure3(a) shows that the two indexes Fe1(ρ = 0) and PY1(ρ = 0.5) increase with the decrease of the rub-impact clearance, but the weighted contribution rate increases from 0 to 1 in the process of the indefinite constant increasing from negative infinity to positive infinity; from Figures 3(a)–3(d), it can be seen that Fe(n) and PY decrease when the rub-impact clearance decreases. The value of Fe(n) and PY(n) n=2,3,4 is also increasing. As seen in Figure 3, the sensitivity of the index Fe is higher than that of PY, but from this, we can also know that in the indefinite constant ρ = 0 or ρ = 0.5, Fe and PY are not the best choice. As shown in Figure 3(b), the peak-weighted contribution rate sensitivity is the highest in the case of different rub-impact clearances, but the indefinite constant corresponding to the peaks changes as the rub-impact clearances change. It can be seen from Figure 3(b) that the trend of the weighted contribution rate Rn(2) is increased first and then decreased, and the curve and the horizontal axis constitute a closed graph when the indefinite constant change. A new index is proposed by using the integral of the indefinite constant ρ of the weighted contribution rate Rn(2). Geometrically, this new index corresponds to the area of the weighted contribution rate Rn(2) and the horizontal axis. This index can better detect the change of rub-impact. Based on this idea, the index will be theoretically derived.The specific expression ofRn(2) can be expressed from equation (16), as follows:(22)Rn2=∫−∞+∞T2jωdω∑i=1N∫−∞+∞Tijωdω=∫−∞+∞G2jω/2ρdω∑i=1N∫−∞+∞Gijω/iρdω,ρ∈−∞,0.Defineψi as ψi=∫−∞+∞Gijωdω, and equation (22) can be written as(23)Rn2=ψ2/2ρ∑i=1Nψi/iρ,ρ∈−∞,0.Equation (23) implies that the index Rn(2) is a function of the indefinite constant ρ, and R2(ρ) is regarded as(24)R2ρ=ψ2/2ρ∑i=1Nψi/iρ,ρ∈−∞,0.According to the analysis in Section2.1, only the first four order NOFRFs of the system are considered here, so equation (24) can be expressed as(25)R2ρ=12ρψ1/ψ2+2/3ρψ3/ψ2+1/2ρψ4/ψ2+1,ρ∈−∞,0.In the practical engineering case,ψi/ψ2, i = 1, 3, 4 in the function R2(ρ) are constants. It can be seen from equation (25) that the function R2(ρ) gradually approaches zero when ρ tends to infinity. According to the analysis mentioned above, the function R2(ρ) is bounded and continuous. So, the function R2(ρ) can be integrated.


## 3.3. Clenshaw–Curtis Integration

In engineering practice, researchers often use numerical integration methods when they solve the integrations of unknown analytic functions. There are many methods for numerical integration. A method with higher precision will be adopted in this paper, namely, Clenshaw–Curtis quadrature formula [28]. Using this method to numerically integrate the index Rn(2), the basic theory is that for the integrand function f(x), its definite integral expression on the interval [−1, 1] is(26)If=∫−11fxdx.The Clenshaw–Curtis quadrature formula is defined as(27)Hnf=∑k=0nθ•kfζk,where ∑• indicates that the first and last items in the formula are halved. ζk=coskπ/n, θk=4/n∑j=0n1/1−j2•cosjkπ/n, k = 0, 1, 2, ..., n, j = 1, 2, 4, ..., n, ζk, θk are the quadrature node and the quadrature coefficient, respectively. Regarding the accuracy of Hn(f), the following theorem is given.Theorem 1.
Iff satisfies the condition of f∈Cr−1,1,r>1, the error is expressed as(28)enf=Hnf−If=on1−r.On a given finite integration interval [a, b], there is an integral of f(x),(29)If=∫abfxdx,where a, b are noninfinite constants and f(x) has a certain smooth degree. The following processing is often done before using Clenshaw–Curtis integration for numerical integration:(30)x=a+b2+b−a2t,that is, the integral interval is converted to [−1, 1] by variable substitution.However, in practical engineering problems, the researchers found that the results obtained did not achieve the desired by linear transformation in most cases. The reason for this is that the trend of the function at the endpoints is different. Considering a class of function,f(x) has a significant change in the vicinity of a in the interval [a, b], and the change in the vicinity of b is slow. Such as f(x) = e−x, when a = 1 and b = 50, the ideal numerical integral is to take more integral nodes near a, and at the same time, in order to reduce the amount of calculation, there are fewer integral nodes near b. Considering the Clenshaw–Curtis numerical integration in the interval [−1, 1], the integral nodes are symmetric about the origin, and there are many nodes at the ±1 and fewer nodes near the origin. With Clenshaw–Curtis numerical integration, linear transformations cannot meet our requirements. Therefore, in this case, a new method of variable substitution will be considered, namely, nonlinear variable substitution.Extending the linear transformation to nonlinearity, equation (30) is rewritten as(31)x=git=b−a2it+1i+a,−1≤t≤1,where i ≥ 1 and i ∈ Z. It can be seen from equation (31) that git can convert the interval [a, b] to [−1, 1], and then perform Clenshaw–Curtis numerical integration on f(x). When i = 1, git is a linear variable substitution; when i = 2, 3, 4, …, git is a nonlinear variable replacement. Under nonlinear transformation, rewrite equation (29) to(32)If=∫−11fgitgi′tdt=b−ai2i∫−11fb−a2it+1i+at+1i−1dt.As can be seen from equation (32), hit=fb−a/2it+1i+at+1i−1 and f(x) have the same smooth order. Then, for equation (32), the Clenshaw–Curtis quadrature formula is used to obtain(33)Inif=b−ai2i∑k=0nf•b−a2iζk+1i+aθkζk+1i−1.In the case of nonlinear variable substitution, the quadrature nodeξi=b−a/2iζk+1i+a,k=0,1,2,...,n is denser near the end point a and sparse at the end point b. This is what we expected.


## 3.4. Second-Order Weighted Contribution Rate Integral Index (RI)

According to the analysis in Section3.2, the function R2(ρ) is integral. The integral expression is(34)IR=∫abR2ρdρ.From Figure3(a) in Section 3.2, it can be concluded that the value of Rn(2) is close to 0 when the indefinite constant is within (−∞, −30] or [5, +∞). So, the integral interval of equation (34) is taken [−30, 5], that is, a = −30, b = 5. Substituting equation (25) into equation (33), equation (34) is rewritten as(35)IniR=35i2i∑k=0nR2•352iζk+1i−30θkζk+1i−1.DefineRI=IniR.In order to verify the reliability of the Clenshaw–Curtis quadrature formula, the authors use the int infinite integral function in MATLAB to solveR2(ρ), and the integration interval is negative infinity to positive infinity. The results are compared with the Clenshaw–Curtis integration method used in this paper, as shown in Table 4. The results show that the selection of the integral interval in the Clenshaw–Curtis method is reasonable, and the results obtained by the Clenshaw–Curtis quadrature method are very accurate.Table 4
Comparison of Clenshaw–Curtis integral and int integral function results.
Clearance (μm)
Clenshaw–Curtis integral
Int infinite integral
I
Clenshaw
−
Curtis
−
I
int
Normal
0.0670145365
0.0670145578
2.1289 × 10−8
100
0.3580265528
0.3580266286
7.5884 × 10−8
80
0.5712722586
0.5712723441
8.5445 × 10−8
60
0.6289452527
0.6289452534
7.0059 × 10−8
40
0.6855516020
0.6855515858
1.6181 × 10−8
20
0.7381875132
0.7381875624
4.9269 × 10−8Rub-impact stiffness and friction coefficient are the two main parameters affecting rub-impact. The change of these two parameters will affect the change of detection indexRI. In order to study the influence of these two parameters on the results, the authors studied the variation of the second-order weighted contribution rate integral index RI with the rub-impact clearances when these two parameters change, as shown in Figures 4 and 5. They illustrated the trend of the index RI of the system under different rub-impact stiffness and friction coefficient. The results show that the index RI will increase with the decrease of the rub-impact clearance when the rub-impact stiffness and friction coefficient of the system increase. The variation of the index RI will gradually increase, which means that the nonlinearity of the system will be increased.Figure 4
The trend ofRI under different rub-impact clearances and rub-impact stiffnesses.Figure 5
The trend ofRI under different rub-impact clearances and friction coefficients.In order to verify that the second-order weighted contribution rate integral indexRI is more suitable for detecting the rub-impact faults of the rotor, according to equation (33), the index RI under different rub-impact clearances is calculated. The amount of change of RI, Ne, Fe1, Fe2, Fe3, and Fe4 under different rub-impact clearance is shown in Table 5. From this, we can see that the indexes RI, Ne, and Fe1 have the largest change with the reduction of the rub-impact clearance. In order to more intuitively show the advantages of the index RI, we compare the three indexes in the form of a histogram, as shown in Figure 6. The results show that the indices RI and Ne increase with the increasing rub-impact clearance, but the index Fe1 decreases as the rub-impact clearance decreases. However, the change of the second-order weighted contribution rate integral index RI is more obvious. Especially in the initial stage of the rub-impact fault, this index has strong superiority. These indicate that the second-order weighted contribution rate integral index RI proposed in this paper is more suitable for detecting the rub-impact fault of the rotator-stator.Table 5
The amount of change ofRI, Ne, Fe1, Fe2, Fe3, and Fe4 under different rub-impact clearances (μm).
Indexes
100
80
60
40
20
ΔRI
0.2910
0.2132
0.0577
0.0566
0.0526
ΔNe
0.0142
7.388 × 10−3
7.413 × 10−3
7.041 × 10−3
6.754 × 10−3
ΔFe1
2.771 × 10−3
1.700 × 10−3
1.810 × 10−3
1.807 × 10−3
1.803 × 10−3
ΔFe2
2.242 × 10−3
1.619 × 10−3
1.616 × 10−3
1.613 × 10−3
1.610 × 10−3
ΔFe3
3.949 × 10−4
8.183 × 10−5
8.164 × 10−5
8.150 × 10−5
8.135 × 10−5
ΔFe4
1.349 × 10−4
6.657 × 10−7
1.119 × 10−4
1.117 × 10−4
1.116 × 10−4Figure 6
Comparison of simulation results of indexesRI, Fe1, and Ne under different clearances.


## 4. Experimental Verification

The experimental equipment is shown in Figure7. In order to obtain the first four orders of NOFRFs of the system, as described in Section 2.1, the system requires two excitations of the same frequency but different intensities. In this paper, the two different strength excitations were obtained by changing the imbalance of the rotor system. In the rotor system, if the imbalance of the system is different, the centrifugal force of the system will also change, which will ultimately affect the vibration amplitude of the system. These two excitations of the same frequency can be obtained at the same speed. Therefore, when the data were collected in the experiment, it was necessary to change the imbalance of the rotor system at the same speed to meet the requirements. A different number of balanced bolts were used to change the imbalance in the rotor system. The rotor speed in the experiment was the same as the simulation, namely, 1200 rpm.Figure 7
Experimental device: 1, eddy current sensor; 2, shaft; 3, rub screw; 4, locknut; 5, motor; 6, dc power; 7, coupling; 8, bearing seat; 9, magnetic base; 10, disc.It can be seen from Figure8 that the index Fe1 changes slowly and the index Ne’s magnitude and the relative variation are small when the rub-impact degree is increased in the rotor rub-impact system. Meanwhile, the variation in the index RI is very obvious. Therefore, RI is a better index to detect system rub-impact faults. Figure 8 also shows that the experimental data and the simulation results basically coincide, which indicates that the simulation results are credible.Figure 8
Comparison of experimental results ofRI, Fe1, and Ne under different severities.


## 5. Conclusion

A new method of detection of the rub-impact fault of the rotor system is introduced in this paper. It introduces the concept of the order of the nonlinear output frequency response functions and uses the NOFRFs to amplify the signal characteristics. Therefore, the impact of higher-order NOFRFs on the system is amplified. On the basis, a new index called the second-order weighted contribution rate integral indexRI is proposed in this paper. After comparative analysis, it is known that the sensitivity of this index to the rub-impact faults of the rotor system is much higher than that of the other two indexes (Fe and Ne). In other words, RI is a better index to detect the rub-impact faults of the rotor system. On the other hand, simulation and experimental results verify that the efficiency of the index RI is in detecting system faults. The second-order weighted contribution rate integral index RI is of great significance in the study of fault diagnosis. The performance of this index in other mechanical faults and structural damage detection [32] will be discussed in future studies.



---

*Source: 3524948-2019-07-02.xml*