Research Article  Open Access
Research on Fault Feature Extraction Method Based on NOFRFs and Its Application in Rotor Faults
Abstract
Rubimpact 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 domainbased 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 highorder 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 highorder NOFRFs on the nonlinearity of the system so as to extract its fault characteristics. The new index RI 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.
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 rubimpact faults between the rotor and stator.
Scholars have done a lot of research on the diagnosis of rubimpact 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 rubimpact 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 rubimpact 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 rubimpact 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 rubimpact 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 rubimpact faults diagnosis. Zhou et al. [10] studied the intrawave frequency modulation phenomenon and extracted the fault characteristics of the rubimpact from the vibration signal by using the MVNCMD method. Most of the above studies only qualitatively describe the fault characteristics of the rubimpact. The feature extraction of rubimpact faults, especially the weak fault features at the beginning of the rubimpact, 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 secondorder 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 thirdorder 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 kernelsPSO (KKPSO) method that identifies the Volterra kernel and used the kernel function of the Volterra series to establish a rotorbearing fault diagnosis system. The results indicate that the KKPSO 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 onedimensional 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 inputindependent 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 PSOadaptive 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 rubimpact fault system. In order to find out a suitable index for online monitoring of rotor rubimpact faults, a new index based on Clenshaw–Curtis integration method [28] is used in this paper, which is the secondorder weighted contribution rate integral index RI. The simulation and experimental results show that the index is more suitable for detecting the rotor rubimpact 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 index RI is proposed, which is used to detect the rubimpact 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 secondorder 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 rubimpact 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 aswherewhere 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, Y_{n}(jω) represents the nthorder output response of the nonlinear system in the frequency domain, and H_{n}(jω) is the nthorder 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 thatits specific expression is defined as
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
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 form u(t) = Asin (ωt), the system frequency domain output can be simplified [29] as
System frequency domain input is defined as
System NOFRFs can be simplified to
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
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, G_{1}(jω), G_{2}(jω), G_{3}(jω), G_{3}(j3ω), G_{4}(j2ω), and G_{4}(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. Index Fe
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 NOFRFsbased index. Its specific expression is as follows:
Equation (10) shows that the index Fe integrates all the NOFRFs values of the nonlinear system. It has the following characteristic:
Fe(n) can reflect the ratio of the nthorder NOFRFs in the system. For example, Fe(1) reflects the proportion of G_{1}(jω) in the system. If Fe(1) ≈ 1, it means that G_{1}(jω) is dominant, and highorder 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 fifthorder 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. Index Ne
The information entropy Ne [16] fuses all orders of NOFRFs, and its specific expression is as follows:where
Since Ne 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, highorder NOFRFs become smaller and smaller as the order increases, which is difficult to extract. Highorder NOFRFs have a significant effect on the detection of system faults. Therefore, it is necessary to extract the nonlinear features of highorder NOFRFs. Based on this, a new method is proposed in this paper, which is to weight the highorder NOFRFs and increase their contribution rate to the system. The specific weighted process is as follows:where T_{n}(jω) represents G_{n}(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:
It will increase the contribution of highorder 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
The highorder 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 highorder will be magnified under the condition . Only in this way, the weighted coefficients at the denominator position is less than 1, namely, . Therefore, the weighted contribution rate Rn solves the problem that the highorder 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 P_{Y}, and P_{Y} 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 P_{Y}. 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. SecondOrder Weighted Contribution Rate Integral Index
3.1. Estimation of NOFRFs of RotorStator RubImpact 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 fixedpoint rubimpact system will be used as an example in this paper. The schematic diagram of the Jeffcott rotor fixedpoint rubimpact system is shown in Figure 1. 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.

Each shaft segment of the rotor system is simulated with Timoshenko beam, and the generalized coordinates of the beam element is u = (x_{A}, y_{A}, θ_{A}, θ_{A}, x_{B}, y_{B}, θ_{B}, θ_{B})^{T}. A and B represent the endpoints of the beam element, respectively, where x_{A}, y_{A}, θ_{A}, θ_{A} are the displacement and rotation angle of the x and y directions of the A node, respectively, and x_{B}, y_{B}, θ_{B}, θ_{B} are the displacements and rotation angles of the x and y directions of the B node, respectively.
The fixedpoint rubimpact of rotorstator is mainly simulated in this paper; in other words, rubimpact will appear between a certain point on the rotating shaft and a stator. The friction torque generated by the rubimpact is ignored, and only the normal rubimpact force F_{n} and the tangential rubimpact force F_{r} are considered when the phenomenon of rubimpact appears in rotorstator, as shown in Figure 2. In Figure 2, O_{1} represents the centroid of the rotor when it is static, represents the centroid of the rotor when it is dynamic, O_{2} 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 k_{c} and f_{c}, respectively.
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 rubimpact force F_{rub} between the rotor and the stator can be decomposed into the normal rubimpact force F_{n} and the tangential rubimpact force F_{r}, as follows:where r is the distance between O_{1} and O_{2} and is the Lagrangian multiplier of the i + 1th step in the iterative process, which is expressed aswhere δ_{r} = r − c, ς is prescriptive permeabilitytolerance. 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 rubimpact force can be expressed aswhere the Heaviside function H(δ_{r}) is defined as
Therefore, the governing equation of the rotor rubimpact system is expressed aswhere, M, C, G, and K are the mass matrix, damping matrix, gyroscopic matrix, and stiffness matrix, respectively; F_{e} is the imbalance force vector; and F_{rub} is the rubimpact force vector of the rotor system [31].
The rest of simulation parameters are shown in Table 2. 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.


3.2. Feature Extraction Based on Weighted Contribution Rate of NOFRFs
According to the analysis in Section 2.3, the contribution rate of the highorder NOFRFs 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 rubimpact 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 rubimpact 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 rubimpact of the system when ρ belongs to at this range.
(a)
(b)
(c)
(d)
In addition, Figure 3(a) shows that the two indexes Fe1(ρ = 0) and P_{Y}1(ρ = 0.5) increase with the decrease of the rubimpact 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 P_{Y} decrease when the rubimpact clearance decreases. The value of Fe(n) and P_{Y}(n) is also increasing. As seen in Figure 3, the sensitivity of the index Fe is higher than that of P_{Y}, but from this, we can also know that in the indefinite constant ρ = 0 or ρ = 0.5, Fe and P_{Y} are not the best choice. As shown in Figure 3(b), the peakweighted contribution rate sensitivity is the highest in the case of different rubimpact clearances, but the indefinite constant corresponding to the peaks changes as the rubimpact 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 rubimpact. Based on this idea, the index will be theoretically derived.
The specific expression of Rn(2) can be expressed from equation (16), as follows:
Define as , and equation (22) can be written as
Equation (23) implies that the index Rn(2) is a function of the indefinite constant ρ, and R_{2}(ρ) is regarded as
According to the analysis in Section 2.1, only the first four order NOFRFs of the system are considered here, so equation (24) can be expressed as
In the practical engineering case, , i = 1, 3, 4 in the function R_{2}(ρ) are constants. It can be seen from equation (25) that the function R_{2}(ρ) gradually approaches zero when ρ tends to infinity. According to the analysis mentioned above, the function R_{2}(ρ) is bounded and continuous. So, the function R_{2}(ρ) 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
The Clenshaw–Curtis quadrature formula is defined aswhere ∑^{•} indicates that the first and last items in the formula are halved. , , k = 0, 1, 2, ..., n, j = 1, 2, 4, ..., n, , θ_{k} are the quadrature node and the quadrature coefficient, respectively. Regarding the accuracy of H_{n}(f), the following theorem is given.
Theorem 1. If f satisfies the condition of , the error is expressed as
On a given finite integration interval [a, b], there is an integral of f(x),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: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 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 aswhere i ≥ 1 and i ∈ Z. It can be seen from equation (31) that can convert the interval [a, b] to [−1, 1], and then perform Clenshaw–Curtis numerical integration on f(x). When i = 1, is a linear variable substitution; when i = 2, 3, 4, …, is a nonlinear variable replacement. Under nonlinear transformation, rewrite equation (29) to
As can be seen from equation (32), and f(x) have the same smooth order. Then, for equation (32), the Clenshaw–Curtis quadrature formula is used to obtain
In the case of nonlinear variable substitution, the quadrature node is denser near the end point a and sparse at the end point b. This is what we expected.
3.4. SecondOrder Weighted Contribution Rate Integral Index (RI)
According to the analysis in Section 3.2, the function R_{2}(ρ) is integral. The integral expression is
From Figure 3(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
Define .
In order to verify the reliability of the Clenshaw–Curtis quadrature formula, the authors use the int infinite integral function in MATLAB to solve R_{2}(ρ), 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.

Rubimpact stiffness and friction coefficient are the two main parameters affecting rubimpact. The change of these two parameters will affect the change of detection index RI. In order to study the influence of these two parameters on the results, the authors studied the variation of the secondorder weighted contribution rate integral index RI with the rubimpact 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 rubimpact stiffness and friction coefficient. The results show that the index RI will increase with the decrease of the rubimpact clearance when the rubimpact 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.
In order to verify that the secondorder weighted contribution rate integral index RI is more suitable for detecting the rubimpact faults of the rotor, according to equation (33), the index RI under different rubimpact clearances is calculated. The amount of change of RI, Ne, Fe1, Fe2, Fe3, and Fe4 under different rubimpact 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 rubimpact 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 rubimpact clearance, but the index Fe1 decreases as the rubimpact clearance decreases. However, the change of the secondorder weighted contribution rate integral index RI is more obvious. Especially in the initial stage of the rubimpact fault, this index has strong superiority. These indicate that the secondorder weighted contribution rate integral index RI proposed in this paper is more suitable for detecting the rubimpact fault of the rotatorstator.

4. Experimental Verification
The experimental equipment is shown in Figure 7. 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.
It can be seen from Figure 8 that the index Fe1 changes slowly and the index Ne’s magnitude and the relative variation are small when the rubimpact degree is increased in the rotor rubimpact system. Meanwhile, the variation in the index RI is very obvious. Therefore, RI is a better index to detect system rubimpact faults. Figure 8 also shows that the experimental data and the simulation results basically coincide, which indicates that the simulation results are credible.
5. Conclusion
A new method of detection of the rubimpact 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 higherorder NOFRFs on the system is amplified. On the basis, a new index called the secondorder weighted contribution rate integral index RI is proposed in this paper. After comparative analysis, it is known that the sensitivity of this index to the rubimpact 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 rubimpact 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 secondorder 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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant nos. 51875093 and U1708257), Basic Research Business Fee of the Central University of Education (grant no. N180304017), and China Postdoctoral Science Foundation (grant nos. 2014M551105 and 2015T80269).
References
 G. JacquetRichardet, M. Torkhani, P. Cartraud et al., “Rotor to stator contacts in turbomachines. Review and application,” Mechanical Systems and Signal Processing, vol. 40, no. 2, pp. 401–420, 2013. View at: Publisher Site  Google Scholar
 H. Ma, C. Shi, Q. Han, and B. Wen, “Fixedpoint rubbing fault characteristic analysis of a rotor system based on contact theory,” Mechanical Systems and Signal Processing, vol. 38, no. 1, pp. 137–153, 2013. View at: Publisher Site  Google Scholar
 J. Páez Chávez and M. Wiercigroch, “Bifurcation analysis of periodic orbits of a nonsmooth Jeffcott rotor model,” Communications in Nonlinear Science and Numerical Simulation, vol. 18, no. 9, pp. 2571–2580, 2013. View at: Publisher Site  Google Scholar
 N. Wang, D. Jiang, and K. Behdinan, “Vibration response analysis of rubbing faults on a dualrotor bearing system,” Archive of Applied Mechanics, vol. 87, no. 11, pp. 1891–1907, 2017. View at: Publisher Site  Google Scholar
 L. Hu, Y. Liu, L. Zhao, and C. Zhou, “Nonlinear dynamic response of a rubimpact rod fastening rotor considering nonlinear contact characteristic,” Archive of Applied Mechanics, vol. 86, no. 11, pp. 1869–1886, 2016. View at: Publisher Site  Google Scholar
 M. Behzad and M. Alvandi, “Unbalanceinduced rub between rotor and compliantsegmented stator,” Journal of Sound and Vibration, vol. 429, pp. 96–129, 2018. View at: Publisher Site  Google Scholar
 H. Cao, L. Niu, S. Xi, and X. Chen, “Mechanical model development of rolling bearingrotor systems: a review,” Mechanical Systems and Signal Processing, vol. 102, pp. 37–58, 2018. View at: Publisher Site  Google Scholar
 Y. F. Yang, Q. Y. Wu, Y. L. Wang, W. Y. Qin, and K. Lu, “Dynamic characteristics of cracked uncertain hollowshaft,” Mechanical Systems and Signal Processing, vol. 124, pp. 36–48, 2019. View at: Google Scholar
 N. H. Chandra and A. S. Sekhar, “Fault detection in rotor bearing systems using time frequency techniques,” Mechanical Systems and Signal Processing, vol. 7273, pp. 105–133, 2016. View at: Publisher Site  Google Scholar
 P. Zhou, M. Du, S. Chen, Q. He, Z. Peng, and W. Zhang, “Study on intrawave frequency modulation phenomenon in detection of rubimpact fault,” Mechanical Systems and Signal Processing, vol. 122, pp. 342–363, 2019. View at: Publisher Site  Google Scholar
 Z. K. Peng, Z. Q. Lang, S. A. Billings, and Y. Lu, “Analysis of bilinear oscillators under harmonic loading using nonlinear output frequency response functions,” International Journal of Mechanical Sciences, vol. 49, no. 11, pp. 1213–1225, 2007. View at: Publisher Site  Google Scholar
 H. Xu, N. Wang, D. Jiang, T. Han, and D. Li, “Dynamic characteristics and experimental research of dualrotor system with rubimpact fault,” Shock and Vibration, vol. 2016, Article ID 6239281, 11 pages, 2016. View at: Publisher Site  Google Scholar
 R. K. Pearson, Discrete Time Dynamic models, Oxford University Press, Oxford, UK, 1994.
 X. Jing and Q. Li, “A nonlinear decomposition and regulation method for nonlinearity characterization,” Nonlinear Dynamics, vol. 83, no. 3, pp. 1355–1377, 2016. View at: Publisher Site  Google Scholar
 Q. Li and X. Jing, “A secondorder output spectrum approach for fault detection of bolt loosening in a satellitelike structure with a sensor chain,” Nonlinear Dynamics, vol. 89, no. 1, pp. 587–606, 2017. View at: Publisher Site  Google Scholar
 Q. Li and X. Jing, “Fault diagnosis of bolt loosening in structures with a novel secondorder output spectrum–based method,” Structural Health Monitoring, Article ID 1475921719836379, 2019. View at: Publisher Site  Google Scholar
 B. Bai and L. Zhang, “HOC based blind identification of hydroturbine shaft Volterra system,” Shock and Vibration, vol. 2017, Article ID 6732704, 11 pages, 2017. View at: Google Scholar
 H. Tang, Y. H. Liao, J. Y. Cao, and H. Xie, “Fault diagnosis approach based on Volterra models,” Mechanical Systems and Signal Processing, vol. 24, no. 4, pp. 1099–1113, 2010. View at: Publisher Site  Google Scholar
 X. Xia, J. Zhou, J. Xiao, and H. Xiao, “A novel identification method of Volterra series in rotorbearing system for fault diagnosis,” Mechanical Systems and Signal Processing, vol. 6667, pp. 557–567, 2016. View at: Publisher Site  Google Scholar
 A. K. Swain, S. A. Billings, P. K. Stansby, and M. Baker, “Accurate prediction of nonlinear wave forces: part I (fixed cylinder),” Mechanical Systems and Signal Processing, vol. 12, no. 3, pp. 449–485, 1998. View at: Publisher Site  Google Scholar
 X. J. Jing, Z. Q. Lang, and S. A. Billings, “Parametric characteristic analysis for generalized frequency response functions of nonlinear systems,” Circuits, Systems & Signal Processing, vol. 28, no. 5, pp. 699–733, 2009. View at: Publisher Site  Google Scholar
 O. M. Boaghe, S. A. Billings, L. M. Li, P. J. Fleming, and J. Liu, “Time and frequency domain identification and analysis of a gas turbine engine,” Control Engineering Practice, vol. 10, no. 12, pp. 1347–1356, 2002. View at: Publisher Site  Google Scholar
 Z. Q. Lang and S. A. Billings, “Energy transfer properties of nonlinear systems in the frequency domain,” International Journal of Control, vol. 78, no. 5, pp. 354–362, 2005. View at: Publisher Site  Google Scholar
 X. Xia, J. Zhou, C. Li, and W. Zhu, “A novel method for fault diagnosis of hydro generator based on NOFRFs,” International Journal of Electrical Power & Energy Systems, vol. 71, pp. 60–67, 2015. View at: Publisher Site  Google Scholar
 H. Huang, H. Mao, H. Mao et al., “Study of cumulative fatigue damage detection for used parts with nonlinear output frequency response functions based on narmax modelling,” Journal of Sound and Vibration, vol. 411, pp. 75–87, 2017. View at: Publisher Site  Google Scholar
 Z. K. Peng, Z. Q. Lang, C. Wolters, S. A. Billings, and K. Worden, “Feasibility study of structural damage detection using NARMAX modelling and nonlinear output frequency response function based analysis,” Mechanical Systems and Signal Processing, vol. 25, no. 3, pp. 1045–1061, 2011. View at: Publisher Site  Google Scholar
 H. Mao, W. Tang, Y. Huang et al., “The construction and comparison of damage detection index based on the nonlinear output frequency response function and experimental analysis,” Journal of Sound and Vibration, vol. 427, pp. 82–94, 2018. View at: Publisher Site  Google Scholar
 D. B. Hunter and H. V. Smith, “A quadrature formula of ClenshawCurtis type for the gegenbauer weightfunction,” Journal of Computational and Applied Mathematics, vol. 177, no. 2, pp. 389–400, 2005. View at: Publisher Site  Google Scholar
 Z. K. Peng, Z. Q. Lang, and S. A. Billings, “Resonances and resonant frequencies for a class of nonlinear systems,” Journal of Sound and Vibration, vol. 300, no. 3–5, pp. 993–1014, 2007. View at: Publisher Site  Google Scholar
 V. Barzdaitis, M. Bogdevicius, and R. Didziokas, “Diagnostics procedure for identification of rubs in rotor bearings,” Journal of Vibroengineering, vol. 12, no. 4, pp. 552–565, 2010. View at: Google Scholar
 M. Jiang, Y. Kuang, J. Wu, and X. Li, “Rubimpact detection in rotor systems with pedestal looseness using a nonlinearity evaluation,” Shock and Vibration, vol. 2018, Article ID 7928164, 15 pages, 2018. View at: Publisher Site  Google Scholar
 P. F. Gao, Z. N. Lei, X. X. Wang, and M. Zhan, “Deformation in fatigue crack tip plastic zone and its role in crack propagation of titanium alloy with trimodal microstructure,” Materials Science and Engineering: A, vol. 739, pp. 198–202, 2019. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Yang Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.