Multicrack Localization in Rotors Based on Proper Orthogonal Decomposition Using Fractal Dimension and Gapped Smoothing Method

Multicrack localization in operating rotor systems is still a challenge today. Focusing on this challenge, a new approach based on proper orthogonal decomposition (POD) is proposed for multicrack localization in rotors. A two-disc rotor-bearing system with breathing cracks is established by the finite element method and simulated sensors are distributed along the rotor to obtain the steady-state transverse responses required by POD. Based on the discontinuities introduced in the proper orthogonal modes (POMs) at the locations of cracks, the characteristic POM (CPOM), which is sensitive to crack locations and robust to noise, is selected for cracks localization. Instead of using the CPOMdirectly, due to its difficulty to localize incipient cracks, damage indexes using fractal dimension (FD) and gapped smoothing method (GSM) are adopted, in order to extract the locations more efficiently. The method proposed in this work is validated to be effective for multicrack localization in rotors by numerical experiments on rotors in different crack configuration cases considering the effects of noise. In addition, the feasibility of using fewer sensors is also investigated.


Introduction
Rotors are one of the most important components of rotating machines, which are widely used in many engineering fields, such as turbines, generators, and aeroengines.Cracks in rotors are the most critical and fundamental damage which may lead to a sudden and catastrophic failure of equipment.So it is of vital significance to identify these cracks, in order to reduce maintenance cost and avoid failure of a rotating machine.In view of the importance, crack identification in rotors has been the focus of many investigations in recent decades and numerous papers have been published [1][2][3][4][5][6][7][8][9], but crack localization is still a challenge for operating rotor systems.
A brief review of the relevant studies of crack localization in rotors is given firstly.The localization methods can be classified as model-based and non-model-based methods.It should be noted that model-based methods defined here are the methods which require a mathematical representation of the system under study, for example, a partial differential equation of motion of a rotor as a beam or mass and stiffness matrices of a finite element model of a rotor.
When it comes to model-based methods in crack identification, approaches based on equivalent crack forces will be mentioned at the first beginning.Equivalent crack force methods consider the effects of cracks as equivalent forces applied in intact systems.They have been adopted to identify the location and depth of cracks in rotors by many researchers, such as Pennacchi et al. [10], Lees et al. [11], and Sekhar [12].There are also some other model-based methods.Dong et al. [13] presented a method based on intersection of the contour curves of the first three natural frequencies obtained from a rotor modelled by wavelet finite element method to determine the crack location and depth.
Seibold and Weinert [14] proposed a method in time domain based on a bank of Extended Kalman Filters to realize the crack location and depth identification for an operating rotor.Methods based on pattern recognition can be model-based when the samples for training are obtained from mathematical models rather than experiments.Their main ideas are extracting features sensitive to crack parameters, which then are trained by artificial intellectual methods or machine learning methods to obtain the relationship between the crack parameters and features, and then matching the measured features with the established relationship to determine the crack parameters.Zapico-Valle et al. [15] adopted the artificial neural network which trained by samples gathered from a finite element model to identify crack location and depth of rotors.Söffker et al. [16] compared a modern model-based technique based on a proportional-integral observer with a signal-based technique based on support vector machine using features extracted from wavelets to identify a crack in an operating rotor.Methods based on optimization are also often model-based, because a large number of iterations are required, and it is almost impossible if there is no model.Genetic algorithm was adopted by Saridakis et al. [17], Xiang et al. [18], and He et al. [19] to minimize the difference between real outputs and model outputs to determine the location and depth of a crack in a rotor-bearing system.Cavalini Jr. et al. [20] put forward a crack identification methodology using external diagnostic forces at certain frequencies to obtain the nonlinear combinational resonances which were used as the objective function of a differential evolution optimization code to determine the crack location and depth minimizing the difference between the measured and modelled rotor system.
In contrast with model-based methods, some crack localization methods do not need a mathematical model of the system under study and are thus called non-model-based methods which often just need inputs and outputs of the system, or even outputs only.Rubio et al. [21] used changes in resonant and antiresonant frequencies to detect crack locations in a two-cracked torsional shaft.Rahman et al. [22] utilized the changes in phase angle of frequency response function to identify the location and depth of a rotor with an open crack.Seo et al. [23] proposed a method for open crack localization by comparing the map of the modal constants of the reverse directional frequency response functions with the reference map of the uncracked model.These methods are based on changes of natural dynamic characteristics; though no mathematical model of the system is needed, they often require information from intact systems which sometimes is not convenient to obtain.There are also some non-model-based methods which do not need reference information from intact systems.ODS measured by sensors distributed along rotors was used for crack localization in rotors by Saravanan and Sekhar [24] and Zhang et al. [25].Babu and Sekhar [26] proposed a modified ODS method called amplitude deviation curve to identify cracks in rotorbearing systems.A residual ODS based method considering higher harmonic components of exciting frequency was developed to localize cracks in a rotor by Asnaashari and Sinha [27].Singh and Tiwari [28] proposed an algorithm for crack localization based on the fact that cracks cause slope discontinuities in the shaft deflection.Due to the lower sensitivity of ODS to incipient cracks, some after-treatment techniques were developed, such as wavelets [29], FD [30,31], and GSM [32].
Both model-based and non-model-based methods have their advantages and disadvantages.In this paper, a nonmodel-based method based on POD is proposed to realize multicrack localization for operating rotors.The POD is a multivariate statistical method which aims at obtaining a compact representation of vibration data [33,34].Galvanetto and Violaris [35] investigated the feasibility of POD to localize a crack in beams.Shane and Jha [36] applied POD to detect damage in composite beams.POD combined with radial basis functions was proposed by Benaissa et al. [37] to identify a crack in plate structures.However, to the authors' best knowledge, the POD has not been used for multicrack localization in operating rotor systems.
Due to the breathing phenomenon of cracks during rotation and the difficulty to generate breathing cracks in rotors, a model that reflects the essential behaviour of a crack is vitally important to get the response of cracked rotors more accurately.There are many methods to model a crack.A nonlinear 3D finite element method was adopted to model a breathing crack in a rotor in [38], which may be the most accurate model, but the computation workload is very heavy.In [39], a rigid finite element method was put forward to model a cracked rotor, which also had good accuracy to model a breathing crack.Papadopoulos' review paper [4] elaborated the approach for modelling cracks in rotors based on SERR and showed that the model put forward by Darpe [40] which was also based on SERR could characterise a breathing crack in rotors more accurately, and it had the advantages of allowing general excitations without assuming that the gravitational force was dominant, and the behaviour of the breathing crack was response-dependent instead of being rotation-dependent.So, considering the complexity in computation and accuracy in modelling in relation to other methods, Darpe's method is adopted to model a cracked rotor in this investigation.
In this work, the feasibility of multicrack localization based on POD using FD and GSM in operating rotor systems is validated by numerical investigation.The proposed method is a kind of non-model-based approach and it does not need the knowledge of the undamaged rotors.The rest of the paper is organized as follows.In Section 2, the model of a two-disc rotor-bearing system considering the static unbalance of discs with response-dependent breathing cracks is established by the finite element method.Section 3 presents the multicrack localization method based on POD using FD and GSM.In Section 4, numerical experiments are carried out for the multicracked rotor with different crack configuration cases.Finally, conclusions are drawn.

Cracked Rotor Modelling
A finite element model of the cracked rotor considering bending-torsion coupling introduced by static unbalance is established in this work.A generalized breathing crack model which can represent any crack angles and any types of excitations applied to the rotor is adopted.To model the cracked rotor, the key point is to simulate the crack appropriately and calculate the stiffness matrix of the crack element.After that, through assembling the cracked and uncracked elements, the finite element model of the rotor can be obtained.

Model of a Cracked Shaft Element.
Figure 1 shows a cracked shaft element of length  and radius . 1 - 12 are the loads acting on the 12 degrees of freedom of the two nodes in the element coordinate system --.The local coordinate system   -  -  is defined on the flat crack face to describe the crack cross-section. is the crack angle between the crack face and the shaft centre line (formed by the negative  axis turning to the negative -axis in the counter-clockwise direction);  L is the location of the crack centre in the element coordinate system.CCL is an imaginary line that separates the open and closed parts of the crack which will be used to simulate the breathing of crack.The hatched area corresponds to the open area of the crack.
The flexibility matrix (G 0 ) 6×6 of the uncracked element and the additional flexibility matrix (G c ) 6×6 of the cracked element can be derived based on SERR theory [4], and the detailed expressions can be obtained from [40].
Thus, from ( 2) and ( 4) So the stiffness matrix of the cracked element K ce and the stiffness matrix of the uncracked element K uce are

Breathing Crack Model.
To consider the breathing phenomenon, what matters the most is to describe the variation of crack section.In this paper, CCL method in [40] is adopted to model the breathing crack.This method assumes that the CCL is perpendicular to the crack edge and separates the open and closed parts of the crack which can be seen in Figure 1.And the position of CCL is determined by calculating the opening mode SIF  I by (7), which depends on the crack element nodal forces, so the crack is responsedependent nonlinear.A positive  I corresponds to the open crack state and a negative one to the closed state.And the CCL is located at the position where the sign of  I changes.Once the CCL is ascertained, the stiffness matrix of the crack element can be obtained.
here  I is the opening mode SIF contributed by   .

Equations of Motion of Cracked Rotor-Bearing System.
The rotor-bearing system considered in this work is shown in Figure 2. The rotor is discretized by two-node Timoshenko beam elements.The discs are considered rigid bodies which have three translational and three rotational inertias.And they are added to the mass matrix elements at the corresponding degrees of freedom.Gyroscopic effect of the two discs is also included.The ball bearings are simplified as stiffness and damping, one of which constrains one axial degree of freedom.The torsional degree of freedom in power input end of the rotor is also constrained.The rotating frequency of the rotor is Ω.By assembling the system matrix of cracked elements and uncracked elements, the finite element model can be established.Denote q  as displacement vector of node  having 6 degrees of freedom: The equations of motion in the stationary coordinate system can be written as follows: where M is the system mass matrix, D = M + K is system damping matrix considering the Rayleigh damping, D g is system gyroscopic matrix, K() is system stiffness matrix which will be updated as the crack breathes, F u is excitation due to static unbalance of discs, F g is excitation due to the gravitational force, and F ex is external excitation during operation.
As for the disc located at node , the gravitational excitation vector is The excitation due to static unbalance is where the elements in F u can be expressed as [41]  u = 0; As shown in Figure 2(b),  Ω () is the relative angular displacement between the rotating coordinate and the stationary coordinate, which will be Ω when the rotor is rotating at a constant speed Ω.   () is the torsional angle. is the unbalance orientation angle of the disc.From (12), one can see that the excitation introduced by unbalance is bending-torsion coupled.So, the equations of motion of the cracked rotor are response-dependent nonlinear and the excitation term is bending-torsion coupled.The Newmark method [42] is used to solve the equations numerically.The stiffness and damping matrices and the coupled excitation term are updated at each integration step, and the next time step will not start until the response in the current time step reaches convergence, which means the increment of displacement is less than the tolerance.

POD Based Multicrack Localization Method
3.1.Theory of POD.The mathematical formulation of POD was reviewed in [34] and will be briefly introduced in the following.
Let (, ) be a random field on Π, and it can be written as where () is the mean value part and (, ) is the time varying part.
The goal of POD is to obtain the most characteristic structure () of an ensemble of snapshots (a snapshot at time   is defined as   () = (,   )) of (, ).It is equivalent to find the basis function that maximizes the ensemble average of the inner products between   () and (): where Here ((), ()) = ∫ Π ()()d denotes the inner product in Π; |⋅| denotes the modulus; ⟨⋅⟩ is the averaging operator; ‖()‖ = ((), ()) To reach the maximum, the derivative of (()) should be zero, which is derived as [43] ∫ where ⟨  ()  ()⟩ is the averaged autocorrelation function.
The optimized solution is given by the orthogonal eigenfunctions   () of ( 16), called POMs.The corresponding eigenvalues   are POVs.
The mathematical formulation mentioned above is the continuous form of POD; however, in real practice, the data obtained are discretized in time and space, so discrete realization of POD by SVD is used in this work.
To start with POD, the system response matrix Y which is measured simultaneously by  sensors at different locations needs to be obtained: where  is the sample length.Y corresponds to the discretized form of field (, ) in (13).
As for the discretized data, the averaged autocorrelation function is replaced by covariance matrix which can be estimated by the sample covariance matrix C s ; then the POMs and POVs correspond to the eigenvectors and eigenvalues of C s , respectively.In particular, if the data have a zero mean, C s can be expressed as SVD of Y can be written as where U × is an orthogonal matrix containing the left singular vectors; S × is a pseudo-diagonal matrix with singular values at the diagonal entries; V × is an orthogonal matrix containing the right singular vectors.
According to (19), one can get Then, The idea of multicrack localization based on POD is that the characteristic structure of measured system response of a cracked system will be different from that of an uncracked system, and cracks will introduce local discontinuities in the POMs, while the POMs should be continuous for an uncracked system, where there is no other factor which introduces discontinuity, for example, a large lumped mass.

Damage Indexes from FD and GSM.
When there is no crack, the POMs will be continuous, but discontinuities will be introduced at the locations of cracks.In order to amplify the effect of the discontinuities in localization, damage indexes based on FD and GSM are used.

FD Based on POMs.
The FD of a curve defined by  points (O 1 , . . ., O  ) is estimated by [44] FD = log 10 ( − 1) log 10 ( − 1) + log 10 (/) , Here, dist(⋅, ⋅) denotes the distance between two points.So, the FD of a specific curve is definite and it is a measurement of the complexity of a curve.Generally speaking, places where discontinuities occur will show high complexity, which is the main idea to use FD as damage index of cracks.In order to detect discontinuities in a POM, a sliding window is used to truncate the curve.The FD in every window is calculated to represent the complexity of the local segment falling into the window and a proper window chosen can amplify the local discontinuities of the whole curve.
Let  be the width of the sliding window and  be the sliding step; then the FD in the th window can be expressed as FD () = log 10 ( − 1) log 10 ( − 1) + log 10 ( () / ()) , () = During the process of crack localization,  responds to the location of midpoint in the window.

GSM Based on POMs.
The GSM is a kind of polynomial curve fitting method.It is used to extract the discontinuities induced by cracks in the POMs in this paper.Its main idea is to fit the cracked POM using gapped polynomial to obtain the approximate uncracked POM and then to calculate the difference function between the actual POM and the fitted POM.Large differences indicate presence of cracks.Generally speaking, the order of gapped polynomial is chosen to be three, so the gapped polynomial function (GPF 3  ) at the gapped point O  (  ,   ) can be written as [32] GPF where  0 ,  However, for the crack localization in rotors, the gapped linear interpolation is found to be more efficient.In this case, the gapped polynomial function (GPF 1  ) can be expressed as where  0 and  1 are determined by O −1 , O +1 .Then, two damage indexes are put forward as the squared difference between the gapped polynomial function and the corresponding value of the actual POM:

Numerical Investigation
In order to investigate the multicrack localization methods, numerical experiments are carried out for the rotor-bearing system shown in Figure 2 and its detailed parameter values are given in Table 1, where  and  are calculated by assuming The rotor is discretized into 28 equivalent two-node twelve-degree-of-freedom Timoshenko beam elements, and cracks with different configurations are embedded using the cracked shaft elements.All the cracks considered are transverse ones and the cracks are assumed not to propagate during the short period of excitation while measurement is made.Newmark method is adopted to obtain the responses in time domain.The Newmark constants are 0.25 and 0.5, respectively, the sampling frequency is 5000 Hz or the integration step is 2 × 10 −4 s, and the accuracy of convergence for each step is set to 10 −11 .

Response Characteristics of a Cracked
Rotor.In order to identify the crack locations in the rotor, response characteristics will be studied first.Figure 3 gives the vertical steady-state responses of the rotor in time and frequency domains without any crack, with one crack and with two cracks, respectively, measured from a single sensor located in the 14th element of the rotor.And  represents the frequency corresponding to the rotating speed.The vertical responses of rotor run-up with angular acceleration 10 rad/s 2 in time domain without any crack, with one crack, and with two cracks are shown in Figure 4.The response characteristics are consistent with the results in [45], so one can believe that the model established and its solution are correct.
From Figures 3 and 4 one can see that the presence of superharmonic components (or subharmonic resonances) generated by nonlinearity introduced by breathing of cracks is a clear indicator of a crack, but there is no qualitative difference between responses of a rotor with a single crack and a rotor with double cracks, whether in steady or unsteady state.
Because the crack number cannot be known in advance, crack detection results could be misleading, which shows the difficulties of multicrack localization by measuring the responses just from a single sensor, since no space information is produced.And it can also be concluded that those methods suitable for a single-crack rotor are not always suitable for a multicrack rotor.In view of the difficulties of multicrack localization in rotors using methods without space information, POD is introduced for the operating rotor in the same situation as the rotor which was used to get the responses in Figures 3 and 4, and the first and second normalized POMs are shown in Figure 5.
As it can be seen from Figure 5 that the one-crack and two-crack cases can be identified by POM1 and POM2, while the two cases are difficult or impossible to be distinguished without space information, as previously shown in Figures 3  and 4.And from Figure 5, one can see that the cracks will introduce discontinuities in POMs.Therefore, multicrack localization can be realized by detecting the discontinuity locations in POMs.2, where the relative phase angle is defined as the angle between the positive normal lines of the two-crack tips which is shown as  phi in Figure 2(b), FD and GSM with the CPOM are used, respectively, to localize the cracks.In order to simulate measurement errors, white Gaussian noise is added to the original response y, so the noise-polluted response y N can be expressed as

Localization of a Double-Cracked Rotor Using FD and GSM with the CPOM. Focusing on the cases of the rotor with double cracks of varying depths and relative phase angles at different locations as shown in Table
where  is the length of y.NL is the constant noise level within (0, 1). is the mean value of y. r is an N-length vector of normally distributed random numbers with zero mean and variance equal to 1. Figure 6 is a typical response of a doublecracked rotor without and with noise.

Localization Results and Robustness of the Method.
Without losing generality, case 2 is chosen to determine the CPOM which is the most robust to noise and most sensitive to cracks.All the cases are measured by 29 sensors in the corresponding nodes, except when investigating the effects of sensor numbers.In order to investigate the effect of noise on POMs, a higher noise level of 5% is considered and the POMs from the first order to the forth order are compared with the corresponding unnoised ones in Figure 7. From Figure 7 one can see that the cracks will affect all the first four POMs, but the first two POMs are less sensitive to noise.In addition, the discontinuity locations in higher order POMs are dominated by one of the cracks; for example, POM2 and POM3 are dominated by crack 2, while POM4 is mainly influenced by crack 1. Therefore in view of the robustness to noise and sensitivity to cracks, POM1 is selected as the CPOM to identify multicrack locations for various cases of cracked rotors in Table 2.However, it is still not easy to identify crack locations from the CPOM directly, so aftertreatment methods which can amplify discontinuities in the CPOM are required.
In order to amplify discontinuities in the CPOM, further treatment is performed by GSM and FD.When GSM is used, the cubic and linear gapped interpolations are compared.And the width of sliding window  for FD in ( 23) is set to 3 which is determined by trial and error.From Figures 8(a cubic gapped interpolation can identify the locations roughly, but the resolution is lower and it is more sensitive to noise compared with GSM by linear gapped interpolation as shown Figure 8(a).In addition, multicrack localization result using FD is also quite good in Figure 8(c).So, in the following, GSM by linear gapped interpolation and FD will be used for multicrack localization (see Figures 9-15).From Figures 8-15, one can see that all the double-crack cases are identified correctly and the method based on CPOM using GSM with linear gapped interpolation and FD is robust to noise.In Figure 8, though the two cracks are located correctly, there are two more discontinuities apart from the crack locations which correspond to the locations of the two discs, but these discontinuities are relatively weak.And, fortunately, as the crack depth increases, the discontinuities induced by discs almost disappear.And from Figure 12, one can see that the method is still reliable even when a crack is located in the same element as the disc in case 5.So it can be concluded that crack locations can be identified regardless of the disc locations.Besides, cracks at different locations with different depths can be localized, and the deeper the crack, the larger the corresponding magnitude of the damage indexes, which can be seen in Figures 9, 11, and 13.And one can also see that even if a crack is near a bearing, it can also be localized correctly, as shown in Figure 13.From Figures 10, 14, and 15, one can see that, under the same crack depths and locations, the relative phase angle will change the values of damage indexes.Because the relative phase angle between two cracks will definitely influence the response of the rotor, thus the CPOM will be different.However, the localization results are still quite good, which means that the proposed method is suitable for cracks in rotors with any crack phase angles.

Effects of Sensor Numbers.
In order to investigate the feasibility to reduce sensor numbers, fewer sensors are used to measure the responses of the cracked rotor in case 3. Fifteen sensors are used and the location results using GSM and FD are shown in Figure 16.
As can be seen from Figure 16, the locations of the two cracks are identified correctly and also insensitive to noise, but with lower resolution.As a matter of fact, the number of sensors determines the spatial resolution and thus it will influence the accuracy of crack localization.So, the more sensors are used, the more accurate localization is, in theory.As for the minimal number of sensors, it can be assumed that there are  c cracks (this number is unknown).For GSM by linear gapped interpolation, to cover the worst situation there should be at least 3 c + 1 sensors, shown as Figure 17(a); for FD method with window width of 3, at least 3 c + 3 sensors are required, shown as Figure 17(b).
In practice, when a crack is localized using  sensors and if it is suspected that the accuracy is poor, all these  sensors can be placed around the damage location and the local responses are measured again.This will lead to a more accurate localization.

Conclusions
Numerical investigation is carried out for multicrack localization in rotors based on proper orthogonal decomposition (POD) using fractal dimension (FD) and gapped smoothing method (GSM).A two-disc rotor-bearing system with response-dependent breathing cracks at different locations of varying depths considering the static unbalance of the two discs is established by the finite element method.Through comparing response characteristics of the rotor with a single crack and two cracks, it is observed that it is very difficult or impossible to distinguish a multicrack case from a singlecrack case just based on the response from one sensor.So, proper orthogonal modes (POMs) are extracted by POD from the responses "measured" from sensors distributed along the rotor.Discontinuities are found to have been introduced by cracks at the corresponding locations in the POMs.Considering the sensitivity to cracks and noise, the characteristic POM (CPOM) is selected.Instead of utilizing the CPOM directly, after-treatment techniques of FD and GSM are used to amplify the discontinuities in the CPOM to realize the multicrack localization more effectively.All the localization results for the rotor with cracks at different locations of varying depths based on CPOM using FD and GSM are quite good.And the crack localization method is robust to noise and fewer sensors are still feasible to successfully locate the cracks.In addition, regardless of input excitations, only responses are needed by the proposed method.What is more, no prior knowledge about the model is demanded which is of great significance for rotors with complex structures and complicated boundaries that are difficult to model.Therefore, the method will be useful in real applications.However, vibration-based damage identification relies heavily on measurement technology.For some machines working in hostile environments, such as steam turbines, noncontact heat-and humidity-resistant sensors should be used.Without goodquality vibration data, the proposed method would not work well.

Figure 2 :
Figure 2: (a) Schematic diagram of cracked two-disc rotor-bearing system.(b) Definition of rotating and stationary coordinates.
So, one can see that the eigenvectors or POMs of C s are the left singular vectors of Y, and the eigenvalues or POVs of C s are the squares of singular values of Y divided by .

Figure 3 :
Figure3: Steady-state responses comparison of rotors without any crack, with one crack and with two cracks (the one crack located in the 12th element with depth of 0.2 and the two cracks located in the 11th and 17th elements both with depth of 0.2).

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Run-up responses comparison of rotors without any crack, with one crack and with two cracks.(a) No crack.(b) One crack located in the 12th element with depth of 0.2.(c) Two cracks located in the 11th and 17th elements both with depth of 0.2.

Figure 7 :
Figure 7: The first four POMs of double-cracked rotor in case 2.

Figure 8 :
Figure 8: Localization results of double-cracked rotor in case 1.(a) Localization using GSM by linear gapped interpolation.(b) Localization using GSM by cubic gapped interpolation.(c) Localization using FD.

Figure 9 :Figure 10 :
Figure 9: Localization results of double-cracked rotor in case 2. (a) Localization using GSM by linear gapped interpolation.(b) Localization using FD.

Figure 15 :Figure 16 :
Figure 15: Localization results of double-cracked rotor in case 8. (a) Localization using GSM by linear gapped interpolation.(b) Localization using FD.

Table 1 :
Parameters of the cracked rotor.