#### 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 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 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 index *RI* 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 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 nth-order output response of the nonlinear system in the frequency domain, and *H*_{n}(*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 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*) = *A*sin (*ω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}(*j*3*ω*), *G*_{4}(*j*2*ω*), and *G*_{4}(*j*4*ω*). 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 NOFRFs-based 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 nth-order 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 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. 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, 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: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 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

The high-order 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 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 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 *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. 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 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 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 force *F*_{n} and the tangential rub-impact force *F*_{r} are considered when the phenomenon of rub-impact appears in rotor-stator, 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 rub-impact force *F*_{rub} between the rotor and the stator can be decomposed into the normal rub-impact force *F*_{n} and the tangential rub-impact 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 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 aswhere the Heaviside function *H*(*δ*_{r}) is defined as

Therefore, the governing equation of the rotor rub-impact 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 rub-impact 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 high-order 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 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.

**(a)**

**(b)**

**(c)**

**(d)**

In addition, Figure 3(a) shows that the two indexes *Fe*1(*ρ* = 0) and *P*_{Y}1(*ρ* = 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 *P*_{Y} decrease when the rub-impact 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 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 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. Second-Order 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.

Rub-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 index *RI*. 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.

In order to verify that the second-order weighted contribution rate integral index *RI* 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*, *Fe*1, *Fe*2, *Fe*3, and *Fe*4 under different rub-impact clearance is shown in Table 5. From this, we can see that the indexes *RI*, *Ne*, and *Fe*1 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 *Fe*1 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.

#### 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 *Fe*1 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.

#### 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 index *RI* 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.

#### 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).