Abstract
This paper proposes a vibration-based damage identification method based on 6-dof Jeffcott rotor system, which is based on harmonic balance and Newton-Raphson methods. First, the equations of motion are derived by using energy method and Lagrange principle. The crack model is based on strain energy release rate (SERR) in fracture mechanics and modified to accommodate 6-dof Jeffcott rotor model. Then, Gear’s method is used to solve the vibration responses of nominal and damaged rotor systems. By processing vibration responses, the transfer function shifts between nominal and damaged systems are taken as the input of damage identification algorithm. Finally, damage severity can be correlated with the damage parameter estimated via developed damage identification model. Numerical examples are shown to demonstrate the effectiveness in identifying the breathing crack in the rotor system.
1. Introduction
Fatigue crack is one of the most common failure modes in the rotordynamic system. It may have an impact on the rotor’s stability and dynamical response at the initial stage and subsequently cause catastrophic accident if untended. Hence, early crack detection and identification are very important especially for high speed rotating machinery. Currently, damage identification methods can be classified into several categories [1]: vibration-based identification, model-based identification [2, 3], signal-based identification [4, 5], artificial intelligence based identification [6, 7], and statistics-based identification. Among these methods, vibration-based identification method is widely explored by many researchers due to its feasibility during operation. Structural damage will result in the change of structural parameters as well as modal parameters and the frequency response function (FRF). Hence, the modal parameters (frequencies, mode shapes, and FRF) can be seen as a sign of structural damage. First, due to easy measurement and high accuracy of low-order frequencies, inherent frequency has long been used to detect structural cracks. Cawley and Adams [8] used the changes of natural frequencies to detect the structural damage and Lele and Maiti [9] presented the forward (determination of frequencies of beams knowing the crack parameters) and inverse problem (determination of crack location knowing the natural frequencies) in a Timoshenko beam. Zhao and DeSmidt [10] proposed the Floquet based damage identification methods to identify the location and severity of open cracks in the rotor shaft. However, the measurement difficulty of high orders of frequency may restrict the application of this method. Second, mode is one of the common techniques to identify the crack. Karthikeyan et al. [5] identify the crack in a beam based on free and forced response measurements. Gounaris and Papadopoulos [11] used a method based on the basic observation that the eigenmodes of any cracked structure are different from those of the healthy one. Another characteristic used for damage identification is FRF. Ishida et al. [12] used external excitations of rotating cracked shafts, thus causing the excitation of nonlinear characteristics of the crack in order to identify it. Hwang’s and Kim’s idea was to minimize the measured and calculated FRF data [13] in order to identify the damage.
Besides aforementioned damage identification literatures, a number of studies have been reported regarding the dynamics of cracked Jeffcott rotor system and its indication (parameter change) for damage detection [14–17]. However, not much literature is found regarding how to identify the crack (inverse problem) by using parameter change of Jeffcott rotor. Here, it is important to emphasize the distinction between damage detection and identification. The former involves nondestructively determining if structural damage has occurred and the latter involves not only detecting but also estimating the severity and location of the damage within a structure. To advance the state of the art, this paper modifies a FRF-based damage identification method [18] to identify the severity of breathing crack in Jeffcott rotor system. First, the damaged rotor system is modeled by synthesizing 6-dof Jeffcott rotor and modified breathing crack model. The stiffness reduction caused by breathing crack is updated at every time step and vibration responses of rotor system are solved by Gear’s method [19]. Then, damage identification model is developed and FRF-based damage identification method is modified by using harmonic balance and Newton-Raphson methods. Finally, the damage identification method is demonstrated via several numerical examples to quantify the severity of crack in the rotor system.
2. Jeffcott Rotor Model
A horizontally placed Jeffcott rotor is employed to model the vibration of rotor system and thereafter damage identification of breathing crack on the rotor. Figure 1 shows the rotor system with the circular shaft (length ; radius ) and the disk in the midspan (thickness ; outer radius ) rotating about its longitudinal axis . is the fixed coordinate system and is the rotating system. The relation between the two coordinate systems is shown in
The symbols , , and denote the rotational deflection in , , and direction, respectively. The system is subjected to the gravity force, unbalance force, external torque, and axial force if applicable. Several assumptions are made to reduce the complexity of the system: (1) the shaft is massless and the rotor disk does not affect the stiffness of the massless shaft. (2) The bearings are considered to be infinitely stiff. (3) The aerodynamics or fluid-film cross coupling forces generated at the seals of the rotor due to the circumferential difference in the flow are not modeled in this analysis.
To derive the equations of motion, we start from writing out the displacement of central disk shown in Figure 2. Point and point are the geometry center and mass center of the disk, respectively, at the initial time. At time , they move to positions and , respectively. The eccentricity is defined by the distance between geometry center and mass center as . Hence, the displacement of disk is expressed in fixed frame as [20]:where , , and are the translational displacements in , , and direction, respectively; note that the tilt deflections are included here so that the gyroscopic effect can be added into the more complex model.
The kinetic energy of the system can be written as follows:where “·” means the first derivative with respect to time.
The potential energy of the system is expressed as follows:By substituting the kinetic and potential energy into Lagrange equation, the equations of motion in 6 degrees of freedom can be established [20]:withwhere , , , and are mass, damping, gyroscopic, and nominal stiffness matrix, respectively; is the stiffness reduction matrix induced by breathing crack; is the force vector (detail can be found in Appendix); is the shaft rotating speed, is the rotation angle induced by external torsion, and is the eccentricity phase. The rotating frame is in line with fixed frame at the initial time.
3. Breathing Crack Model
The model is based on the concept of energy release in fracture mechanics which was first proposed by Griffith in 1920. Figure 3(a) shows a crack under different types of loading on a shaft and Figure 3(b) shows the cross section of the shaft at the crack location.
(a)
(b)
The breathing crack model has been widely used in a number of literatures either in 4-dof Jeffcott rotor model [15, 16] or finite element rotor model [21]. The model is modified to accommodate 6-dof Jeffcott rotor model [20].
The flexibility matrix induced by the crack can be obtained bywhere is released strain energy due to crack [1]. Then, the total compliance matrix of the rotor system can be obtained by adding the nominal flexibility matrix to the flexibility matrix induced by crack:Note that the total compliance matrix is derived in the crack coordinate system . Therefore, the stiffness matrix obtained by inversing the compliance matrix is required to perform coordinate transformation to get the stiffness matrix in fixed coordinate system: with
Assuming that the stiffness matrix and vibration response at time step have been obtained, the force vector in the fixed frame at this time step can be calculated by
The SIF for crack opening mode iswithNote that the force should be transformed from fixed frame to rotating frame to calculate the stress intensity factors. The forces , , and are the components in rotating frame corresponding to , , and in fixed frame. The points in the crack area are swept to determine the Crack Closure Line (CCL) [15] by finding the points with zero . CCL is calculated for every time step and the crack stiffness matrix is updated to solve for the vibration response of next time step. Thus, the crack open area at time step can be used to integrate the compliance coefficients and update the stiffness matrix. The vibration response at time step can then be solved by Gear’s method.
4. Damage Identification
The damage identification of small crack is more interesting since the early damage detection is more important to prevent the catastrophic accident. In this section, the damage identification model is built and FRF-based method is modified to identify the severity of the breathing crack. The idea is to utilize the transfer function shifts induced by the crack to estimate the stiffness reduction, which can reflect the crack depth.
The transverse crack induces the stiffness reduction in diagonal terms and stiffness coupling terms as shown in
Figure 4 shows the dimensionless crack stiffness in lateral, torsional, and longitudinal directions in terms of crack depth. The subscript “” stands for the rotating coordinate system. The maximum values of diagonal terms are the nominal stiffness of system which is used to nondimensionalize the parameters. In lateral direction, the stiffness reduction in and direction of rotating frame is increasing with the crack depth and so does the coupling stiffness. The torsional and longitudinal stiffness have not much reduction for small crack (). In addition, the cross stiffness between longitudinal and lateral direction is trivial in comparison with longitudinal stiffness.
(a)
(b)
(c)
4.1. Damage Identification Model
For the small crack, we assume that the lateral stiffness , torsional stiffness , and longitudinal stiffness are constant because of their negligible variation. All the coupling terms are neglected due to their small magnitudes with respect to the diagonal terms. Therefore, the crack stiffness matrix in rotating frame can be treated as the nominal stiffness matrix except 1 stiffness element: lateral stiffness . The damaged stiffness matrix for the identification is expressed as
In damage identification model, a sinusoid function as shown in (17) is used to approximate the actual variation of the lateral stiffness in y direction :
The stiffness matrix in fixed coordinate system is
4.2. Damage Identification Algorithm
The equations of motion are recast into the following state-space form with state vector :with The nonlinear damage perturbation matrix is approximated as the summation of Fourier series and the homogeneous solution of the linear periodically time-varying system can be expressed aswhere is the number of harmonics for damage perturbation matrix and is the number of system harmonics. After performing harmonic balance, the following hyperdimensional problems whose solutions characterize the nominal and damaged system transfer function matrices are given as follows:with and are the hyperdimensional nominal system matrix and damage perturbation matrix with details in Appendix.
Using the Taylor expansion on (23), the damage identification is formulated in terms of the following sensitivity matrix :where is transfer function shifts between nominal and undamaged systems; is the vector of unknown damage parameters which can be estimated in the least-squares sense asUsing (26) as an initial damage estimate, subsequent estimates can be iteratively improved by Newton-Raphson method. The th estimate of the damage parameters can be determined fromwith where is Jacobian matrix.
4.3. Damage Identification Workflow
Figure 5 shows the damage identification workflow which contains 2 parts. The first part simulates the vibration responses of nominal and damaged rotor system via 6-dof Jeffcott rotor model and breathing crack model. The second part numerically tests the damage identification method with 6-dof Jeffcott model, which utilizes the transfer function shifts to estimate crack depth. The transfer function shifts in time domain can be approximated by harmonic Fourier expansion: is the number of harmonics; , , and are Harmonic Fourier Coefficients (HFCs) of transfer function shifts. Finally, crack depth can be estimated by looking up the damage estimates in the crack lookup curve as shown in Figure 6.
(a)
(b)
Figure 6 shows the maximum and minimum stiffness in direction of rotating frame in terms of the crack depth. Obviously, the deeper crack shows the larger stiffness variation. The crack lookup curve plots the actual ratios between the stiffness variation and the nominal stiffness in direction of rotating frame for different crack depths.
5. Results and Discussion
In this section, the damage identification workflow is applied to identify the shallow crack with depths 0.1R, 0.2R, and 0.3R. The breathing cracked rotor system with the parameters listed in Table 1 is used to demonstrate the damage identification results. The rotor is assumed to have 0.01 mm eccentricity which is about 10% of static deflection. The damping ratio can vary between 0.03 and 0.3 for traditional bearings, although a minimum of 0.1 is normally considered as needed for the safe operation of the machine [22]. There is no external torsional or axial excitation applied on the rotor while rotating.
Figures 7, 8, and 9 show the iteration history of damage estimates and its correlation to crack depth for 3 cases, respectively. In this paper, crack location has been assumed to be midpoint of shaft and only crack depth is needed to be estimated to characterize the damage. It can be seen that the damage parameter is estimated by DC terms, 1st and 2nd Harmonic Fourier Coefficients individually. By fitting the damage parameters into the lookup curve the estimated crack depth is close to the actual one for all three cases in numerical simulations. However, the noise effect needs to be taken into account during experimental validation because the high order of harmonics will be disrupted especially when the crack is very small. The method may need further optimization to improve signal-to-noise ratio.
(a)
(b)
(a)
(b)
(a)
(b)
If multicracks exist on the rotor system, Jeffcott rotor model is not able to describe the problem anymore and other rotor models (i.e., finite element rotor model) are needed [21]. It would be another challenge to develop corresponding FRF based damage identification method for such models to identify crack location and depth and phase if necessary.
6. Conclusion
The 6-dof Jeffcott rotor system is built based on the Lagrange principle. The breathing behavior of crack is captured by modified model based on existing breathing crack model. The damage identification model is developed and FRF-based damage identification method was developed to identify the depth of breathing crack by utilizing the Harmonic Fourier Coefficients of transfer function shifts obtained from vibration simulation of rotor model. The case study shows that the crack can be identified by proposed methodology which steps forward from the damage detection/diagnosis to the damage identification.
The paper has proposed a theoretical framework for the damage identification of breathing cracked rotor system. There still exist many challenges between this theoretical work and real industrial applications. The Jeffcott rotor assumes the massless shaft and lump mass at disk which may not be accurate enough for the real structure. The noise generated during the measurement may contaminate the transfer function signals and thus lead to estimation error. These issues deserve further investigation in the future work.
Appendix
Mass matrix is as follows:
Damping matrix is as follows:
Gyroscopic matrix is as follows:
Nominal stiffness matrix is as follows:
Force vector is as follows: where is the mass of disk, and represent the damping coefficient and stiffness in th DOF, is the coupling stiffness between th and th DOF induced by the crack, and and are the mass moment of inertia in principle and tangential axis, respectively.
Stress intensity factors are as follows: with
Compliance coefficients are as follows:
Hyperdimensional nominal system matrix is as follows:
Hyperdimensional damage perturbation matrix is as follows:
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.