On Modeling of Vibration and Crack Growth in a Rotor for Prognostics

Prognostics and health management (PHM) include comprehensive engineering approaches that evaluate the realtime health condition of an asset and predict its future states under the actual operating conditions. There are two main concepts in PHM: diagnostics and prognostics. Detection, isolation and identification of faults are done by diagnostics, while prognostics is concerned with predicting future fault progression of a system by assessing the extent of deviation from its expected normal operating conditions. Mechanical fatigue phenomenon that causes crack initiation and propagation is considered to be a common reason for failure in mechanical parts. The current paper studies mutual effects between the fatigue crack propagation and vibrations of an unbalanced rotor. At the first step, the coupled equations of rotor motion and crack growth are obtained. The Paris– Erdogan law is used for crack growth modeling and the Jeffcott model is used for the rotor. The coupled equations are solved numerically using the Runge–Kutta method. The mutual effects between cycles of loading, excitation frequency and crack growth rate are demonstrated for high cycle and low cycle loadings using numerical simulations. The outputs of the model show the importance of considering coupling between the dynamic response and crack growth for fatigue degradation modeling. The proposed model is capable of calculating the instantaneous crack depth for a given rotary system within a range of excitation frequency and loading cycles. In addition, it can model rotor degradation evaluation that is of significance in future state estimation, prognostics, and lifetime prediction of rotating systems.


INTRODUCTION
Planning an effective maintenance strategy is crucial for safety of mechanical systems. The maintenance strategies have historically changed from post-failure repair to preventive maintenance to Condition Based Maintenance (CBM). Unlike the first two strategies, CBM is a costeffective maintenance approach that makes maintenance decisions only when needed. Recently, CBM has been suggested for many advanced mechanical systems with high reliability requirements (Kim, An, & Choi, 2016).
Prognostics and Health Management (PHM) are the primary concepts to aid CBM. PHM includes engineering approaches for real-time health assessment of a mechanical asset under real operational conditions and predicts the system states evolution based on online data captured form the system. Diagnostics, as the first step in PHM, is concerned with detection, isolation and identification of faults, while prognostics is concerned with predicting the future fault progression and degradation evolution when the captured data shows deviation from ideal operating conditions (Goebel, Daigle, Saxena, Sankararaman, Roychoudhury, & Celaya, 2017).
Mechanical fatigue phenomenon that causes crack initiation and propagation is considered to be one of the most common reasons for failure in engineering systems. Rotating machineries like compressors, turbines, pumps and expanders are under high risk of such fatigue failures due to long runs with rotating motion. Crack initiation and propagation can strongly affect the nonlinear dynamics of a rotating system due to major structural changes that it brings about (Schijve, 2001). Hence, modeling the dynamics of a rotating system under the effects of crack propagation is arguably a crucial step in its future state estimation, prognostics and life prediction. Notable research in this area are addressed in the following. 2 A literature review by Wauer (1990) on dynamics of cracked rotating systems can be considered as an early work in this area. Dimarogonas (1996) reviewed dynamics of cracked structures with focus on rotors. Nelson and Nataraj (1986) analyzed the dynamics of a rotor-bearing system with a transversely cracked rotor. The effects of crack presence on dynamics of the system and consequent stiffness variation and parametric excitation were studied. Papadopoulos (2008) proposed the strain energy release method for obtaining the stiffness coefficients in cracked rotating systems. Kumar and Rostagi (2009) presented a review on approaches applied to modeling dynamics of cracked rotating systems. The major part of related research focused on modeling the cracked rotating systems under negligible crack growth assumption (Darpe, Gupta, & Chawla, 2004;Ebrahimi, Heydari, & Behzad, 2017;Palacios-Pineda, Gómez-Mancilla, Martínez-Romero, & Elías-Zúñiga, 2017;Papadopoulos & Dimarogonas, 1992;Sinou & Lees, 2005). Under this assumption, a wide range of research tried to study the effects of crack on dynamics of a rotor using numerical analyses such as bifurcation, phase portraits and Poincaré maps (Patel & Darpe, 2008;Weiyang Qin, Chen, & Ren, 2004;Qin, Meng, & Zhang, 2003;Saeed & Eissa, 2018;Yiming, Yufang, & Shijian, 2003). Although these studies made notable improvement in modeling of cracked rotors, the negligible crack growth assumption undermines their effects on PHM advances. These approaches include analyses like modeling fault propagation, degradation evaluation and future state estimation when the system goes under a long period of operation. In such conditions, assuming that crack growth is negligible cannot be valid. Hence, both short term and long term behaviors of the system should be taken into account in modeling. This emphasizes the importance of modeling long time degradation and short term vibration in a coupled form for analyzing cracked rotating systems.
In the case of coupling nonlinear vibration and fatigue induced degradation several studies are worthy of mention. Initial research concerned with modeling simple vibratory systems under effects of different excitation types and fatigue induced stiffness degradation (Sobczyk, Perros, & Papadimitriou, 2010;Sobczyk & Trebicki, 2000;Trebicki & Sobczyk, 2008). Oppenheimer and Loparo (2002) presented a physics-based approach for diagnostics and prognostics using integrated observers and life models in a rotating system. Using observers for shaft cracking and imbalance, the number of machine fault strengths and corresponding remaining machine life were determined. Niu and Yang (2018) considered a coupled rotor vibration and crack growth model and studied the effects of the crack growth on stability of the system using numerical simulations.
Based on the literature review, the majority of research studies focused on modeling the cracked rotor under negligible crack growth assumption, so modeling the coupled vibration and crack propagation seems still far from mature.
This paper tries to model precisely the mutual effects between crack growth and nonlinear vibrations of a rotor under excitation of an unbalance force for different scenarios of loading cycle. Also, the degradation evolution is modeled by considering these mutual effects. The interdependency between cycles of loading, excitation frequency and crack growth rate are demonstrated for high cycle and low cycle loadings using numerical simulations. It is shown that the proposed model can calculate the instantaneous crack depth for a given rotating system within a range of excitation frequency and loading cycles.
The remaining parts of this paper are organized as follows. In Section 2, a mathematical model of the cracked rotor is derived. In Section 3, results of numerical simulations under different loading scenarios are discussed. The conclusion is presented in Section 4. Fig. 1a shows a flexible Jeffcott rotor including a massless elastic shaft and a rigid disk with transverse surface crack located at the mid span of the rotor. In order to model the dynamics of the system, rotating and fixed coordinates are considered. and are axes of the fixed coordinate and and are the rotating coordinate axes positioned at the crack cross-section according to Fig. 1b. It is assumed that the shaft is under excitation of an unbalanced force with eccentricity of ε at an angle β respect to axis. The rotational speed is notated by ω and ( ) is the instantaneous rotation angle. For obtaining the equations of motion the following assumptions are made (Abbasi, Khadem, Bab, & Friswell, 2016;Darpe, Chawla, & Gupta, 2002).

MATHEMATICAL MODELING
1. Amplitude of torsional and axial vibrations are negligible, so only lateral vibration and corresponding transverse crack growth are taken into account. 2. Rotor damping is of linear viscous type. 3. The only source of excitation is the eccentricity of the unbalance force. 4. The initial crack front is a straight line and propagates in the same form. 5. The applied stress range is well below the yield stress, so linear elastic fracture mechanics is valid in crack modeling. 6. The disk thickness is negligible, so the displacement at the crack location is the same as the disk displacement. 7. Only plane strain relations are considered for the crack front propagation under the excitation force. 8. It is assumed that the magnitude of the gravity force dominates the unbalance force.

Rotor Dynamic Equations
The equations of motion of the cracked rotor using Newton's second law of motion is derived as (Patel & Darpe, 2008): where and show the stiffness and damping coefficients of the rotor. and are displacements along and directions, respectively. and represent direct stiffness coefficients along and directions in the fixed coordinate, respectively and and are the cross-stiffness coefficients in this coordinate. The stiffness coefficients in the fixed coordinate can be related to those in the rotating coordinate using a proper transformation matrix as (Patel & Darpe, 2008): where is the transformation matrix; y and represent the direct stiffness coefficients along and directions in the rotating coordinate, respectively and yz and are the cross-stiffness coefficients in the rotating coordinate. The transformation matrix can be written as: Also, displacement in the rotating coordinate and fixed coordinate can be related to each other using the same transformation matrix as: Given the assumption on the dominancy of the gravity force over the unbalance force, it is assumed that the open part of the crack changes continuously with shaft rotation. Thus, a harmonic breathing function (F(θ)) for modeling the crack opening and closing versus shaft rotation is considered as: The relation between stiffness coefficients of the uncracked rotor with those of the cracked one using the breathing function can be written as: where 0 is stiffness of uncracked rotor and ̂y and ̂z are stiffness coefficients in the rotating coordinate when the crack is fully open along and directions, respectively. For calculating these stiffness coefficients, a cross section of the cracked rotor according to Fig. 2 is considered. In this figure is the crack depth; is the crack depth; is the rotor diameter and ′ = √ 2 − (2 ) 2 is the depth of the crack where is equal to . Using the strain energy method the flexibility of the cracked cross section and the corresponding stiffness, coefficients can be obtained as (Chasalevris & Papadopoulos, 2006;Patel & Darpe, 2008 For obtaining the stiffness coefficients in the rotating coordinates when the crack is fully open (̂ and ̂z ) the following relations can be used (Chasalevris & Papadopoulos, 2006;Patel & Darpe, 2008 where ̂y , ̂, ̂y z and ̂ are the flexibility coefficients when the band of integration is set to the depth of the crack when it is fully open (Chasalevris & Papadopoulos, 2006;Papadopoulos & Dimarogonas, 1992;Patel & Darpe, 2008).

Crack Propagation Equation
Forced vibration of the rotor under effects of the unbalance force and consequent fatigue stress can result in slow propagation of the rotor transverse crack. Experimental tests on fatigue process show that the range of the stress intensity factor range ( ) of the crack front is the most effective parameter in crack propagation. When the stress intensity factor is higher than a threshold value, the crack growth can be modeled using Paris-Edgard equation as (Shih & Chen, 1997): where is number of loading cycles; and are empirical constants that are dependent on material properties and environmental effects. For obtaining the stress intensity factor range, the bending force applied on the crack cross section is calculated as (Chasalevris & Papadopoulos, 2006;Patel & Darpe, 2008): where and are bending forces along and directions in the rotating coordinate, respectively. The corresponding stress relations can be obtained as (Chasalevris & Papadopoulos, 2006;Patel & Darpe, 2008): where and z are stress terms along and directions in the rotating coordinate, respectively. is the crack cross section moment of inertia. So, the total stress intensity of the crack front ( I ) will be According to Eq. 10 it can be seen that forces applied on the crack cross section and dynamic responses are interdependent. Also, the mutual effects between the dynamic response and crack growth can be seen in Eqs. 9 and 13, in which the crack depth affects the cross-section stiffness while the dynamic response changes the stress intensity range. Using the transform matrices in Eq. 3 and the relations found for stiffness of the crack cross section in Eq. 9, the equations of motion can be rewritten as (Chasalevris & Papadopoulos, 2006;Patel & Darpe, 2008): The equations of motion can be written in dimensionless form using the following parameters (Patel & Darpe, 2008):

Solving Coupled Equations of Rotor Vibration and Crack Propagation
The rate of crack propagation and corresponding degradation process is low compared to the vibration time history. Hence, in order to couple the rotor vibrations and crack propagation, additional assumptions on numerical simulation time scales should be made. It is assumed here that for low cycles of loading, the crack growth is negligible, so its depth is constant. High cycle loadings can be separated to low cycle loading steps. The length of low cycle loading is assumed to be 1000. Hence, the crack depth changes after each 1000 cycles of loading and it is assumed to be constant meanwhile. Based on this assumption, the second order differential equations governing the motion equations (Eq. 16) are transformed to the first order ones in the state space; they are then solved numerically step by step using the Runge-Kutta method. Each step lasts for 1000 cycles and the maximum stress intensity factor range for each step is obtained. Then, the crack depth at the end of each step is calculated using Eq.
10. This depth is set as the initial crack depth of the next step of 1000 cycles and this process continues until loading is over. It should be noted that at each step, in order to eliminate the transient behavior of the system, the numerical results of the first 100 cycles of loading are eliminated. In the numerical simulation, the following parameters are used (Patel & Darpe, 2008

RESULTS AND DISCUSSION
In this section, mutual effects between the crack growth and dynamic responses are analyzed under different loading scenarios.

Low Cycle Loading with Negligible Crack Growth
In this section, it is assumed that the crack growth is negligible, so the initial crack depth remains constant (1mm), and the rotor undergoes 1000 cycles of loading. Figs. 3 and 4 show the frequency response of the cracked rotor along the and directions, respectively. The horizontal axis is ratio of the excitation frequency to the natural frequency of the uncracked rotor increased by a constant step. At each frequency ratio step, the motion equations (Eq.16) that represent a classical two degree of freedom forced vibration equations are solved numerically using Runge-Kutta method. Then, the amplitude of response within the range of 1000 cycle is considered as the frequency response at this step. This procedure is repeated for a range of frequency ratio with lower and upper bands of 0.001 and 2, respectively. The outputs of this procedure form Figs. 3 and 4. The diagrams show that the motions of the rotating system in the and directions are periodic and similar but not identical. The positive direction of the vertical displacement has been assumed to be upward. Given the direction of the gravity force and effects of the static deflection, the vertical frequency response is negative in the low frequency range. Also, it can be seen that the peak amplitude of response of the cracked rotor happens in resonance condition in which the excitation frequency gets close to the natural frequency of uncracked rotor and the corresponding ratio is close to 1. It can be inferred that, for the crack depth of 1mm, the variation in cracked rotor stiffness is not significant, so the frequency response and resonance conditions of the cracked and uncracked rotor are very similar.
6 Figure 3. Vertical frequency response of the rotor with constant crack depth of 1mm Figure 4. Horizontal frequency response of the rotor with constant crack depth of 1mm The effects of the crack on dynamics of the system are intensified for higher values of the initial crack depth. Figs. 5 and 6 show the effects of initial crack depth on the response amplitude of the cracked rotor along the and directions, respectively. The results are obtained for several values of the initial crack ranging from 1mm to 14mm. The system goes through low cycle range of loading (1000 cycles) at the resonance frequency of the uncracked rotor condition. It can be seen that for higher values of the initial crack depth, the amplitude of the rotor response will be lower in both and directions. It shows that increasing the initial crack depth escalates the cracked rotor stiffness degradation and increases the equivalent flexibility of the rotor, so the resonance frequency reduces. For example, considering the case 1 and 2 with the crack depth of 1mm and 8mm, respectively. The resonance frequency of the case 1 is the post-resonance frequency in case 2. Hence, a drop in amplitude of response between two cases can be seen in Figs. 5 and 6.

High Cycle Loading and Effect of Crack Growth
The crack propagation is highly slower than exposure rate of nonlinear behaviors of the rotor. Upon crack initiation, each cycle of loading increases the crack depth a small amount. In the low cycle loading the summation of this changes can be neglected, while it can cause fatigue failure over the high cycle loadings. Majority of industrial rotating machines like gas turbines and compressors undergo continuous and longtime loadings. Given the wide application range of these machines from power generation to aircrafts, failure can bring about many major consequences. In this section, the interdependency between the crack growth and nonlinear behaviors of the rotor over high cycle loadings is investigated. As mentioned earlier, for high cycle loading (here, over 1000 cycles) the loading history can be separated to low cycle loading steps with the length of 1000 cycles. In the low cycle loading steps the crack depth is assumed to be constant. Upon calculating the maximum stress intensity factor range in each step, the crack depth is updated and is considered as the initial depth of the consecutive step. Fig. 7 shows crack growth trajectories versus cycles of loading for different values of the frequency ratio, which is defined as the ratio of excitation frequency to the natural frequency of the uncracked rotor. The initial crack depth is 1mm and the simulation is over when the crack depth reaches 6mm. Within the frequency ratio range of 1 to 2.1, it can be seen that as the system comes close to the resonance condition, where the ratio is equal to 1, the crack growth rate increases. It is clear that the most rapid crack growth happens when the frequency ratio is equal to 1 that represents the resonance condition. Figure 7. Crack growth trajectories for different frequency ratio values The mutual effects between the crack growth and dynamics of rotor can be illustrated by comparing appropriate indices. With this aim, Fig. 8 demonstrates the relationship between normalized vibration peak amplitude and normalized maximum stress intensity factor range. In this figure, horizontal axis shows the ratio of excitation frequency to natural frequency of the uncracked rotor and results are obtained in a range of loading cycle in which the crack depth grows from 1mm to 6mm. These two normalized parameters can be considered as indices showing the magnitude of the dynamic response and the crack growth rate, respectively. It can be seen that the variation of these two indices are in accordance with each other. The maximum normalized stress intensity factor range happens where the peak amplitude of the response happens. This confirms maximum crack growth happens when the system is in the resonance condition. This coupled model shows how the dynamic response can affect the crack growth and fatigue degradation process and vice versa. This model can precisely follow the degradation evolution of the rotor according to Figs. 7 and 8 by considering the interdependence between crack growth and dynamic responses. Calculating the crack depth at a specific excitation frequency and cycle of loading is a primary step in prediction of crack growth that can be done by this model. In addition, another notable output is stiffness degradation evolution modeling due to fatigue phenomena (see Fig. 7). The evolution data can be used in future state estimation of the system that is so important from failure prediction and reliability aspects. Hence, the outputs of this model are of significance to prognostics and lifetime prediction of rotating systems.

CONCLUSION
This paper has been concerned with modeling the coupling between the fatigue degradation due to crack propagation and dynamic response of an unbalanced cracked rotor. Given the significant difference between the time rate of fatigue degradation and dynamic response, loading history has been separated to high cycle and low cycle loading steps. It has been assumed that high cycle loadings are consecutive steps of low cycle loadings in which the crack depth is constant. The mutual effects between the dynamic response and crack propagation rate in high cycle loadings have been demonstrated using numerical simulations. The proposed model is capable of calculating the instantaneous crack depth for a given rotary system within a range of excitation frequency and loading cycles. This model can provide detailed insight on degradation evaluation and future state of the system. The outputs of this model are of significance to prognostics and lifetime prediction of rotating systems. Our ongoing research will focus on improvements of this model in the following aspects.
The length of the low cycle loading step can be variable using an adaptive approach. Significant variations in amplitude of the response can undermine the validity of negligible crack growth assumption. Hence, an adaptive approach can set the length of low cycle loading step based on dynamic response variation rendering the corresponding stress intensity factor range more realistic.
The crack growth trajectories of high cycle loadings have been obtained under the assumption that all of them start with the same initial crack depth (1mm). Not making this assumption and enabling the model to follow the crack growth for any arbitrary initial crack depth would be a notable improvement in the generalization of the model. Upon this step, the model will be able to provide a dataset containing initial crack depth, excitation frequency, cycle of loading and instantaneous crack depth. This dataset can be used in combination with data driven techniques in a hybrid prognostic approach.
The capability of this model in tracking the degradation evaluation can be used along with a diagnostics approach. Further, a health management strategy can be designed in which at the first step a diagnostics approach can detect the fault in the form of a crack in the rotor; then, the proposed model of this paper can predict the fault propagation (crack growth in this study) and degradation evaluation in the prognostic step.