The Improved WNOFRFs Feature Extraction Method and Its Application to Quantitative Diagnosis for Cracked Rotor Systems

During its operation, a rotor system can be exposed to multiple faults, such as rub-impact, misalignment, cracks and unbalancing. When a crack fault occurs on the rotor shaft, the vibration response signals contain some nonlinear components that are considerably tougher to be extracted through some linear diagnosis methods. By combining the Nonlinear Output Frequency Response Functions weighted contribution rate (WNOFRFs) and Kullback–Leibler (KL) divergence, a novel fault diagnosis method of improved WNOFRFs is proposed. In this method, an index, improved optimal WNOFRFs (IOW), is defined to represent the nonlinearity of the faulty rotor system. This method has been tested through the finite element model of a cracked rotor system and then verified experimentally at the shaft crack detection test bench. The results from the simulation and experiment verified that the proposed method is applicable and effective for cracked rotor systems. The IOW indicator shows high sensitivity to crack faults and can comprehensively represent the nonlinear properties of the system. It can also quantitatively detect the crack fault, and the relationship between the values of IOW and the relative depth of the crack is approximately positively proportional. The proposed method can precisely and quantitatively diagnose crack faults in a rotor system.


Introduction
Rotor systems can be considered the most common and critical mechanical components in rotating machinery. They are widely used in various engineering fields, such as aerospace, transportation and electric power [1]. However, due to the more demanding requirements for the mechanical equipment, the structure of the rotor is more complicated, so it is prone to breakdowns in operation due to assembly errors or failure of its critical components [2,3]. Crack faults in rotor shafts are a common phenomenon that can lead to catastrophic failure and great economic loss if not identified in time [4,5]. Therefore, it is essential to detect and identify crack fault in rotor systems.
Vibration signals play an important role in reflecting the health status of the rotor system [6][7][8][9]. At present, many studies have focused on fault diagnosis methods for cracked rotor systems through vibration signals. R. Gradzki et al. [10] presented a rotor fault detection method through the auto-correlation and power spectral density functions of the rotor response signals and confirmed the sensitivity of this method, as well as its reliability. L. Xiang et al. [11] employed the orbit morphological characteristics to diagnose cracks in rotor systems. Its effectiveness at detecting crack information was verified through the dynamic model analysis and experimental study. Lu et al. [12] proposed a crack identification method for breathing crack identification of rotors based on the Kriging surrogate model by extracting the super-harmonic nonlinear characteristics. The results of the finite element method and experiments verify that this method is effective, accurate and robust for breathing crack identification in rotors. A. Hajnayeb et al. [13] created a feature space including the features of the largest Lyapunov exponent, approximate entropy and correlation dimension for fault diagnosis in rotors due to cracks. Saeed N. A. et al. [14] utilized the bifurcation diagram, Poincaré map, frequency spectrum and whirling orbit to investigate the Jeffcott rotor system. This strategy can not only inspect the existence of a crack but also predict the crack size based on the qualitative change of the system's whirling motions when the crack size increases. Yang et al. [15] proposed the hybrid MVA method and used the cracked bench rotor to validate its effectiveness. These methods can identify the crack fault of the rotor system, but all of them cannot quantitatively describe the crack to detect its depth.
The Frequency Response Function (FRF) is the basic theory in the analysis of a linear system. The FRF represents the relationship of the input and output for a linear system in the frequency domain. The theory of the FRF has been widely used in engineering practice [16]. However, most engineering systems are complicated and nonlinear, so they cannot be simply described through the linear model. In the late 1950s, the Generalized Frequency Response Functions (GFRFs) were developed in the nonlinear case based on the theory of the Volterra series. The GFRFs were defined as the multidimensional Fourier transformations of the Volterra kernel functions. They are used to represent the nonlinear characteristics in the frequency domain [17]. However, the calculation of high-order GFRFs is quite tedious because it will involve a large amount of algebra or symbolic manipulations [18]. Therefore, the GFRFs actually do not have a broad range of applications in engineering practice.
To tackle the above problem, Lang put forward the Nonlinear Output Frequency Response Functions (NOFRFs) on the basis of the Volterra series [19]. NOFRFs are a significant extension of the FRF theory for linear systems to nonlinear situations. It reveals the rationale of nonlinear phenomena well and provides a certain theoretical basis for analyzing nonlinear systems in the frequency domain [20]. One of the most attractive points of the NOFRFs is its one-dimensional nature, which has many advantages, so it is extensively used in engineering practice to solve some difficult problems, such as in the field of structural health monitoring and fault diagnosis [21]. Peng et al. [22] used the NOFRFs to detect cracks in beams, as well as the nonlinear component's position of periodic structures [23]. Based on the Nonlinear Auto Regressive with eXegenous Input (NARX) model and the NOFRFs, Bayma et al. [24] designed a structural damage detection method for detecting cracks in beam structures. By introducing the evidence theory to the NOFRFs, Cao et al. [25] put forward a fault diagnosis method for complicated nonlinear systems by combining the NOFRFs and evidence theory. This approach was employed to diagnose the fault of the transmission system of a numerical control machine tool. Zhao et al. [26] put forward a method for detecting early damage to nonlinear systems by combining the NOFRFs with the convolutional neural network and long short-term memory network (CNN-LSTM) model, which was verified through a cantilever beam with a breathing crack.
Some studies are devoted to NOFRF-based indexes to enhance the effect of the NOFRFs for defect detection and fault diagnosis. Peng et al. [27] put forward an index, Fe(n), toward the Nonlinear Auto-Regressive Moving Average with eXegenous Inputs (NARMAX) model in damage detection. The experiments of aluminum plates with holes and cracks indicated that the index Fe(n) of the inspected structure is considerably different from the damage-free structure. Huang [28] put forward the index Ne, which integrates all orders of the NOFRFs by introducing information entropy. The experiments of the plate specimens were conducted to reveal the higher sensitivity of the index Ne. Mao et al. [29] introduced Kullback-Leibler (KL) divergence to the NOFRFs to propose the index NOFRFs-KL (NKL) for fatigue damage identification. The effectiveness of NKL was verified by an experiment on a train wheelset. For improving the sensitivity of the Fe and Ne indexes, Liu et al. [30,31] put forward the method of the Nonlinear Output Frequency Response Functions weighted contribution rate (WNOFRFs). The index of the second-order optimal weighted contribution rate Rm was defined for feature extraction of the rub-impact of the rotor system. Another index, K, which is based on Rm, was also defined. Its superiority was verified through the mathematical model and experiments on a cracked rotor system. Li et al. [32] proposed the index J by introducing the Clenshaw-Curtis numerical integration method and the minimum mutual entropy principle to the WNOFRFs for diagnosing the rotor rub-impact.
Further studies on the feature extraction method of the WNOFRFs and the associated index K indicated that this method can improve the performance of crack detection. However, in this method, only the second-order WNOFRFs are focused on, rather than comprehensively considering the other orders of the WNOFRFs. The increments of the index K are a bit small when the degree of the crack fault increases, so it is not sensitive enough for crack fault detection. The values of K and the degree of the crack are in a linear relationship to some extent rather than in a positive proportion, as might be expected. In the present study, a fault diagnosis method known as the improved WNOFRFs is proposed, which is employed to diagnose rotor cracks. An associated index denoted as IOW is proposed to represent the nonlinear characteristics of the system. The effectiveness of the method was verified through a mathematical model of the cracked Jeffcott rotor system and an experimental rotor system with crack faults.

NOFRFs Theory
For a stable linear system, its output frequency response can be expressed as where U(jω) and Y(jω) are the input and output spectra, respectively, which are the Fourier transforms of the input function u(t) and the output function y(t) of the system. For linear systems, this expression shows that the possible output frequencies are the same as the frequencies in the input. Considering the relationship in the time domain of the input and the output of a nonlinear system which are stable at zero equilibrium, this can be represented by the Volterra series: where u(t) denotes the input function of the nonlinear system at discrete time t, y(t) is the corresponding output function, h n (τ 1 , · · · , τ n ) is the nth discrete time Volterra kernel, N denotes the maximum order of the system nonlinearity and τ is the time delay. The output frequency response of the nonlinear system under general excitation can be represented by where Y(jω) represents the out spectrum of the system and Y n (jω) represents the nth order output frequency response of the nonlinear system, which is defined as Equation (4) is a natural extension of the well-known linear relationship (expressed in Equation (1)) to non-linear cases, which is much more complicated than in linear system cases. It reveals how the nonlinear mechanisms operate on the input spectra to produce the system output frequency response. ω 1 +···+ω n =ω H n (jω 1 , · · ·, jω n ) under the condition of ω 1 + · · · + ω n = ω in the n-dimensional hyperplane. H n (jω 1 , · · · , jω n ) is the Fourier transform of h n (jω 1 , · · · , jω n ) and denotes the nth GFRFs of the nonlinear systems [33]. Its expression is Under the condition of U n (jω) = This is substituted into Equation (4), and the nth output frequency response of the nonlinear system can be rewritten as Therefore, the relationship between the output spectrum and the input spectrum can be rewritten as [34] The first four-order NOFRFs are generally sufficient to reflect the nonlinear characteristics of a nonlinear system, so only the first four-order NOFRFs are considered.
Therefore, Equation (8) can be further expressed as [30] As shown in Equations (9)-(12), if a nonlinear system was excited twice by two harmonic input signals that had different amplitudes, the first four-order NOFRFs could be calculated through the least square algorithm.
The Fourier transforms of these two harmonic inputs are A The output spectra are expressed as G 1 (jω F ) and G 3 (jω F ) can be calculated as follows: The remaining values of the NOFRFs can be calculated in a similar way.

The Principle of WNOFRFs
On the basis of the associated index Fe proposed by Peng [27], Liu [30,31] introduced the weighting coefficient n ρ and proposed the feature extraction method known as the NOFRFs weighted contribution rate (WNOFRFs). Its expression is where the weighting coefficient n ρ amplifies the high-order G n (jω) of the nonlinear system. Further studies found that the second-order weighted contribution rate Rn2(ρ) is a good index for indicating the degree of the fault, and Rn2(ρ) is defined as In Equation (16), the maximum value of the curve Rn2(ρ) was denoted as the secondorder optimal WNOFRFs Rm, and the Newton-Raphson iterative method was applied to solve the x-coordinate corresponding to Rm, which is named as the optimal fitness factor ρ max . Its expression is as follows: On the basis of Rm, an associated index K was defined for solving the problem that the sensitivity of Rm changes significantly at different stages of crack growth [31]. The index K can be calculated through This index is used to represent the system's nonlinearity, and the severity of the faults can be distinguished by Rm to some extent. However, for the index K, only the second-order WNOFRFs are focused on, rather than comprehensively considering the other orders of the WNOFRFs. The increments of K are a bit small when the degree of the crack fault increases, so it is not sensitive enough to crack faults, and the values of K and the degree of the crack are in a linear relationship to some extent rather than in a positive proportion, as might be expected.

Improved WNOFRFs Fault Diagnosis Method
To explore an index that has a high enough sensitivity and its values having a certain relationship with the degree of faults to precisely and quantitatively diagnose faults, KL divergence is combined with WNOFRFs.
The matching degree between the referenced probability density distribution f (x i ), i = 1, · · · , N and the target probability density distribution g(x i ), i = 1, · · · , N can be measured by calculating their KL divergence. If the difference between f (x i ) and g(x i ) is larger, their KL divergence D KL ( f g) will be larger [29]. The expression of the KL divergence is given by According to Equation (19), the discrepancy between two WNOFRFs, such as the WNOFRF of the healthy rotor system and the WNOFRF of the cracked rotor system, can be represented through calculating the KL divergence. Therefore, to explore a superior fault diagnosis method, a fault diagnosis method of the improved WNOFRFs is proposed based on the WNOFRFs and the KL divergence. Its expression is as follows: In Equation (20), Rn h (n)(ρ) represents the nth-order WNOFRFs of the healthy rotor system and Rn t (n)(ρ) represents the nth-order WNOFRFs of the target rotor system.
If the target rotor system does not have any crack faults, the curve Rn t (n)(ρ) approximately coincides with Rn h (n)(ρ), so for a specific ρ = ρ 0 , log Rn h (n)(ρ 0 ) Rn t (n)(ρ 0 ) is near zero. At the same time, the value of IW(ρ 0 ) is close to zero. If the target rotor system contains crack faults, Rn t (n)(ρ 0 ) is bigger than Rn h (n)(ρ) for a specific ρ = ρ 0 , so log Rn h (n)(ρ 0 ) Rn t (n)(ρ 0 ) < 0. When the crack fault gets worse, the value of Rn t (n)(ρ 0 ) increases, so for a specific ρ = ρ 0 , the value of log Rn t (n)(ρ 0 ) increases. Therefore, the value of IW(ρ 0 ) gradually increases as the crack fault gets worse.

Mathematical Model of the Cracked Rotor System
The mathematical model in this study was considered to employ the Jeffcott rotor system [35] with a crack fault, which is shown in Figure 1.  Table 1. The dynamic equation of the Jeffcott rotor system can be expressed as where M is the mass matrix, C is the damping matrix, G is the gyro matrix, K is the stiffness matrix and F is the excitation vector.
There is a crack in a certain shaft segment, the radius of which is R and the length of which is L. The cracked shaft element is illustrated in Figure 2. Figure 2a is a diagram of the cracked shaft element. Figure 2b is the cross-sectional view of the crack. In Figure 2a, the element is loaded with axial forces P 3 and P 9 , shear forces P 1 , P 2 , P 7 and P 8 , bending moments P 4 , P 5 , P 10 and P 11 and torques P 6 and P 12 . In Figure 2b, a represents the depth of the crack, and b is its width. The relative depth of the crack is represented as a/R. The stiffness matrix of the crack-free shaft element is where the transition matrix T = According to the Castigliano theorem, the general strain energy of the cracked shaft element is expressed as where U 0 is the strain energy of the crack-free shaft element and U c is the strain energy of the cracked shaft element. Based on the fracture mechanics, U c can be expressed as where J(A) is the strain energy function. According to the Paris theory, the flexibility of the cracked shaft element caused by P i can be expressed as where u 0 represents the flexibility of the crack-free shaft element and u c represents the additional flexibility due to cracking [36]. The flexibility coefficients of the cracked shaft element can be expressed as Under the assumption that the direction of the crack is on the same side as the mass eccentricity, only the mixed situation of a Mode I crack caused by the bending moment P 4 and a Mode II crack caused by the shear force P 2 is considered, as the direction of P 2 and P 4 is the same as the direction of the crack. Based on this, the strain energy of the crack-free shaft element is expressed as Thus, the flexibility matrix of the crack-free shaft element is given as The stress intensity factors K I4 and K II4 , corresponding to the bending moment P 4 and the shear force P 2 , are defined as where σ 4I = 2(P 2 α−P 4  and h = 2 √ R 2 − w 2 . Therefore, the additional strain energy caused by the crack with depth a and width b in Figure 2b is Thus, the local additional flexibility coefficients c 44 and c 22 caused by the bending moment P 4 and the shear force P 2 are expressed as The additional flexibility matrix caused by the crack is given by Thus, the flexibility matrix of the full-open crack is expressed as The stiffness matrix of the full-open crack is where the transition matrix The flexibility matrix of the semi-open crack is expressed as The stiffness matrix of the semi-open crack is The stiffness matrix of the cracked shaft element varies with the time rotation angles due to the opening and closing characteristics of the crack. This variation [37] is expressed as where The Newmark-β algorithm is employed to solve the dynamic equation of the Jeffcott rotor system at different crack faults. The rotor system is excited twice using different unbalanced forces. The depth of the cracks is set to 1 mm, 2 mm, 3 mm or 4 mm, so the relative depth (representing the ratio between the crack depth and the radius of the shaft) is 0, 0.2, 0.4, 0.6 or 0.8, and the dynamic responses at node 8 of the crack-free rotor system and cracked rotor systems with different degrees of crack faults at the rotation speed of 1600 rpm are shown in Figure 3. The time domain responses are shown in Figure 3a, and the frequency domain responses of the mathematical model are shown in Figure 3b. It can be observed from Figure 3a that the domain response of the crack-free rotor system was a standard sine signal. With the degree of the crack fault getting worse (the relative depth of the crack ranges from 0 to 0.8), the amplitude of the time domain responses gradually increased both at the peak and at the trough. This change was not distinct at the peak but visible at the trough. As crack faults became gradually more serious, the peak values of the 1× component in the frequency domain response increased from 0.0257 mm to 0.02721 mm (0.0257 mm, 0.02579 mm, 0.02596 mm, 0.02596 mm and 0.02721 mm). The change was not visible, and the peak values of the 2× component in the frequency domain response increased from 6.325 × 10 −6 mm to 1.196 × 10 −3 mm (6.325 × 10 −6 mm, 3.54 × 10 −5 mm, 1.226 × 10 −4 mm, 3.672 × 10 −4 mm and 1.196 × 10 −3 mm). The relative change was distinct. In Figure 3b, it is found that the 3×, 4× and 5× components also appeared in the condition of the most serious crack fault, which was due to changes in the local flexural stiffness of the cracked rotor [38]. However, these changes were insufficient to recognize the depth of the cracks through simple spectral analysis and some linear methods. Therefore, a novel nonlinear method needs to be explored to evaluate the cracked rotor system.

Improved WNOFRFs for the Cracked Rotor System
Dynamic responses from the mathematical model of the cracked Jeffcott rotor system under different severities of the crack faults and corresponding excitations were used for the fault diagnosis through the method of improved WNOFRFs. First, the WNOFRFs were calculated according to Equation (15). The curves of the WNOFRFs (Rn) under different crack faults are shown in Figure 4. It can be found from Figure 4 that the changing trends of the first four orders of WNOFRFs were relatively distinct. Rn1 monotonically increased from 0 to 1 with the increase in ρ. By contrast, Rn4 monotonically decreased from 1 to 0. When ρ was constant, the values of Rn1 only slightly decreased as the crack depth became larger, while Rn4 only had very small changes. Figure 4b indicates that Rn2 first rose from 0 to the maximum, which corresponds with the optimal fitness factor. Then, it declined to zero. The maximum values of Rn2 under different crack fault severities increased from zero, with the crack fault getting worse. As shown in Figure 4c, Rn3 followed the same trend as Rn2. However, its maximum values manifested a reverse trend as the relative depths of the crack faults deepened.
Next, IW1-IW4 were further obtained according to Equation (20). The curves of IW1-IW4 under different crack depths are presented in Figure 5. It can be seen from Figure 5 that the changing trend of IW1, IW2, IW3 and IW4 was similar to that of Rn2. The maximum values of each order of improved WNOFRFs under different conditions were also similar to that of Rn2.
To integrate all orders of WNOFRFs, IW1, IW2, IW3 and IW4 were summed to yield the IW, the curves of which are shown in Figure 6. It can be obviously observed from Figure 6 that the changing trend of IW was first from 0 to the maximum. Then, it declined to zero. The maximum values of IW under different relative depths of crack faults increased from zero, with the fault getting worse.

The Index of IOW for the Cracked Rotor System
In this subsection, a novel index representing the nonlinear characteristics of the system, the improved optimal WNOFRFs (IOW), is proposed. Its expression is where ρ opt represents the optimal fitness factor. In Figure 6, when ρ = ρ opt , the corresponding maximum value of every curve is denoted as IOW.
In this paper, the index IOW was used to detect the cracked rotor system. The IOW and K of the finite element model under different crack conditions were calculated according to Equations (41) and (18). The values of IOW and K are shown in Figure 7. From Figure 7, a conspicuous increasing trend could be found for the values of IOW and K (IOW was from 0 to 3.410, and K was from 0.160 to 0.604), with the relative depth of the cracks increasing. These two indexes both can recognize crack faults in the rotor system. When comparing these two indexes, the values of IOW were larger than K under the same degree of cracking, so the index IOW was more sensitive.
Then, we explored the relationship between the values of the indexes and cracks. The relative depths were set from 0 to 0.8 in increments of 0.01. The values of the indexes were calculated, and then linear fitting was performed on the values. The fitting results of indexes IOW and K are illustrated in Figures 8 and 9. Figure 8a is the fitted line of IOW, and Figure 8b is its corresponding residuals plot. Figure 9a is the fitted line of K, and Figure 9b is its corresponding residuals plot. The fitted line of index IOW was a/R = 0.229IOW + 0.0301, the RMSE was 0.01868, and the R-square was 0.9938. Meanwhile, the fitted line of index K was a/R = 1.772K − 0.3176, the RMSE was 0.02234, and the R-square was 0.9911. According to the rule that the value of the R-square is closer to one, the better the result is, it is concluded that the fitting result of IOW is better than that of K. The RMSE of IOW was smaller than that of K. Therefore, IOW could diagnose crack faults more accurately.

Process of the Improved WNOFRFs
According to the principle of improved WNOFRFs, the process of the improved WNOFRFs for diagnosing crack faults of the rotor system is summarized as follows and shown in Figure 10. Step 1: Collect the vibration response signals of the crack-free rotor system being subjected to different magnitudes of unbalanced excitations.
Step 2: Calculate the NOFRFs of the crack-free rotor system G h n (jw), n = 1, 2, 3, 4 through the excitation signals and response signals. Then, calculate its NOFRFs weighted contribution rate Rn h (n)(ρ), n = 1, 2, 3, 4 through the following: Step 3: Collect the vibration response signals of the cracked rotor system being subjected to different magnitudes of unbalanced excitations.
Step 6: Determine the index IOW through Equation (41). Then, conduct fault diagnosis of the inspected system.

Experiment Study
To evaluate the performance of the fault diagnosis method of improved WNOFRFs, experiments were carried out to diagnose the cracked rotor system. The experimental test bench of the cracked Jeffcott rotor system is illustrated in Figure 11. Different magnitudes of unbalanced forces were implemented by installing different numbers of bolts on the disk. Different degrees of crack faults were implemented by setting different crack depths at the same location of the shaft by a fatigue machine. The depths of the crack faults were 1 mm, 2 mm, 3 mm and 4 mm, as shown in Figure 12. Therefore, the relative depths were 0.2, 0.4, 0.6 and 0.8. Then, the dynamic responses of the experimental rotor system operating at 1600 rpm under various degrees of cracked faults were collected by displacement transducers.  The curves of IW1-IW4 under different crack faults are illustrated in Figure 13. The improved WNOFRFs (IW) were obtained through summing IW1-IW4. Figure 14 shows the curves of IW.  When comparing Figures 13 and 14 with Figures 5 and 6, it was concluded that the performance of the improved WNOFRFs in the experimental study was similar to that in the simulation analysis. It can be seen from Figure 13 that each of IW1-IW4 rose from 0 to the maximum, which corresponded to the optimal fitness factor. Then, they declined to zero. The peaks of IW1-IW4 were proportional to the crack depths, and the curves of IW had the same trend with the curves of IW1-IW4. Thus, the changing trend of IW in the experiment was also similar to that in the simulation. Figure 15 shows the values of K and IOW. It can be seen that the index K increased from 0.5594 to 0.9138, and the index IOW increased from 0 to 1.9610. Obviously, the index IOW was greater than the index K when the rotor system was in the same crack fault condition. Therefore, the result of the experiment was consistent with that of the simulation. It can be concluded that the index IOW for the crack fault was more sensitive than K. Therefore, the method of improved WNOFRFs could effectively diagnose the cracked rotor system.

Conclusions
A new fault diagnosis method known as improved WNOFRFs was developed by combining the WNOFRFs and KL divergence. A new index, denoted as IOW, was proposed to detect a cracked rotor system. The effectiveness of the improved WNOFRFs was proven through the simulation model and experiments of the Jeffcott rotor system. The following conclusions can be drawn: (1) When the relative depth increased from 0 to 0.8, K increased from 0.160 to 0.604, and IOW went from 0 to 3.410, so IOW had the higher sensitivity to crack faults; (2) IOW could comprehensively represent the nonlinear properties of the system by integrating all orders of the WNOFRFs; (3) IOW could quantitatively represent the crack fault, the minimum relative depth of which could be detected with increments of 0.01; (4) In terms of the simulation model used in this study, the relationship between the values of IOW and relative depths of the crack faults was a/R = 0.229IOW + 0.0301. Therefore, the values of the index IOW had an approximately positive proportional relationship with the relative depth of the crack.
Therefore, compared with the original WNOFRFs, which mainly focus on the second weighted contribution rate, the improved method was superior in the fault diagnosis of the cracked rotor system. The index IOW manifests many advantages for detecting crack in a rotor system. Future studies will be focused on the effectiveness of the improved WNOFRFs in identifying other kinds of faults of rotor systems, as well as the relationship between the values of IOW and the degree of the fault. Moreover, fault diagnosis of other nonlinear engineering systems will be considered.