Abstract

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 CPOM directly, 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.

1. 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 [19], 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 rotor-bearing 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 non-model-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.

2. 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.

2.1. Model of a Cracked Shaft Element

Figure 1 shows a cracked shaft element of length and radius . 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); 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 of the uncracked element and the additional flexibility matrix of the cracked element can be derived based on SERR theory [4], and the detailed expressions can be obtained from [40].

According to the assumption of linear elasticity, the total flexibility matrix of the crack element is the sum of and :With the flexibility matrix, the stiffness matrix can be derived considering static equilibrium of crack element (as shown in Figure 1) and definition of the stiffness.

The nodal force vector of an element can be expressed in terms via transformation matrix asHere,According to Hooke’s law, the left term of (4) is the vector of displacement along .Thus, from (2) and (4)So the stiffness matrix of the cracked element and the stiffness matrix of the uncracked element are

2.2. 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 by (7), which depends on the crack element nodal forces, so the crack is response-dependent nonlinear. A positive 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 changes. Once the CCL is ascertained, the stiffness matrix of the crack element can be obtained.here is the opening mode SIF contributed by .

2.3. 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 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 is the system mass matrix, is system damping matrix considering the Rayleigh damping, is system gyroscopic matrix, is system stiffness matrix which will be updated as the crack breathes, is excitation due to static unbalance of discs, is excitation due to the gravitational force, and is external excitation during operation.

As for the disc located at node , the gravitational excitation vector isThe excitation due to static unbalance is where the elements in can be expressed as [41]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.

3. 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 aswhere 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 :Here denotes the inner product in ; denotes the modulus; is the averaging operator; denotes the norm of a function.

By introducing Lagrange multiplier, the optimization problem can be expressed asTo 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 which is measured simultaneously by sensors at different locations needs to be obtained:where is the sample length. 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 ; then the POMs and POVs correspond to the eigenvectors and eigenvalues of , respectively. In particular, if the data have a zero mean, can be expressed asSVD of can be written aswhere is an orthogonal matrix containing the left singular vectors; is a pseudo-diagonal matrix with singular values at the diagonal entries; is an orthogonal matrix containing the right singular vectors.

According to (19), one can getThen,So, one can see that the eigenvectors or POMs of are the left singular vectors of , and the eigenvalues or POVs of are the squares of singular values of divided by .

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.

3.2. 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.

3.2.1. FD Based on POMs

The FD of a curve defined by points () is estimated by [44]Here, 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 asDuring the process of crack localization, responds to the location of midpoint in the window.

3.2.2. 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 at the gapped point can be written as [32]where , , , and are determined by , , , and .

However, for the crack localization in rotors, the gapped linear interpolation is found to be more efficient. In this case, the gapped polynomial function can be expressed aswhere and are determined by , .

Then, two damage indexes are put forward as the squared difference between the gapped polynomial function and the corresponding value of the actual POM:

4. 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 modal damping ratios of the first two modes being 0.005 and 0.01. And the first critical speed is calculated in the no-crack condition.

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.

4.1. 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/s2 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.

4.2. 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 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 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 can be expressed aswhere is the length of y. is the constant noise level within (0, 1). is the mean value of . 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 double-cracked rotor without and with noise.

4.2.1. 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 after-treatment 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) and 8(b), one can see that multicrack localization using GSM by 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 915).

From Figures 815, 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.

4.2.2. 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 cracks (this number is unknown). For GSM by linear gapped interpolation, to cover the worst situation there should be at least sensors, shown as Figure 17(a); for FD method with window width of 3, at least 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.

5. 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 single-crack 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 good-quality vibration data, the proposed method would not work well.

Abbreviations

CCL:Crack closure line
CPOM:Characteristic proper orthogonal mode
FD:Fractal dimension
GSM:Gapped smoothing method
ODS:Operational deflection shape
POD:Proper orthogonal decomposition
POM:Proper orthogonal mode
POV:Proper orthogonal value
SERR:Strain energy release rate
SIF:Stress intensity factor
SVD:Singular value decomposition.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.

Acknowledgments

This study is partly supported by the National Natural Science Foundation of China (51405399) and the Fundamental Research Funds for the Central Universities (DUT16RC(3)027) and carried out by the first author during his visit to the University of Liverpool, sponsored by the China Scholarship Council.