Crack identification for rotating machines based on a nonlinear approach

Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Crack identification for rotating machines based on a nonlinear approach A.A. Cavalini Jr a,n, L. Sanches a, N. Bachschmid b, V. Steffen Jr a a LMEst – Laboratory of Mechanics and Structures, Federal University of Uberlândia, School of Mechanical Engineering, Av. João Naves de Ávila, 2121, Uberlândia, MG 38408-196, Brazil b Department of Mechanical Engineering, Politecnico di Milano, Via La Masa, 1, Milano 20156, Italy

a r t i c l e i n f o

abstract

Article history: Received 22 February 2016 Accepted 23 February 2016

In a previous contribution, a crack identification methodology based on a nonlinear approach was proposed. The technique uses external applied diagnostic forces at certain frequencies attaining combinational resonances, together with a pseudo-random optimization code, known as Differential Evolution, in order to characterize the signatures of the crack in the spectral responses of the flexible rotor. The conditions under which combinational resonances appear were determined by using the method of multiple scales. In real conditions, the breathing phenomenon arises from the stress and strain distribution on the cross-sectional area of the crack. This mechanism behavior follows the static and dynamic loads acting on the rotor. Therefore, the breathing crack can be simulated according to the Mayes' model, in which the crack transition from fully opened to fully closed is described by a cosine function. However, many contributions try to represent the crack behavior by machining a small notch on the shaft instead of the fatigue process. In this paper, the open and breathing crack models are compared regarding their dynamic behavior and the efficiency of the proposed identification technique. The additional flexibility introduced by the crack is calculated by using the linear fracture mechanics theory (LFM). The open crack model is based on LFM and the breathing crack model corresponds to the Mayes' model, which combines LFM with a given breathing mechanism. For illustration purposes, a rotor composed by a horizontal flexible shaft, two rigid discs, and two self-aligning ball bearings is used to compose a finite element model of the system. Then, numerical simulation is performed to determine the dynamic behavior of the rotor. Finally, the results of the inverse problem conveyed show that the methodology is a reliable tool that is able to estimate satisfactorily the location and depth of the crack. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Rotating machine Crack identification Combinational resonances Multiple scales method

1. Introduction Visual examination, ultrasonic tests, and dye penetrant inspection are some examples of nondestructive techniques widely used for crack detection in rotors. These methods have proved to be costly, since satisfactory results rely on detailed and periodic inspections. Significant research effort has been directed in recent years to online monitoring techniques, i.e., based on vibration

n

Corresponding author. E-mail address: [email protected] (A.A. Cavalini Jr).

http://dx.doi.org/10.1016/j.ymssp.2016.02.041 0888-3270/& 2016 Elsevier Ltd. All rights reserved.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

2

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

signals measured during rotor operation for the detection or identification of cracks location and severity (size and depth of cracks) [1–3]. The techniques based on vibration measurements are recognized as promising SHM tools (Structural Health Monitoring). According to Dimarogonas [4], over 500 papers were published on this subject between the 1980s and 1990s and this number is still very significant in the present days. There are several Structural Health Monitoring techniques, the so-called SHM techniques, proposed in the literature for crack detection in rotating machines. Among them, the ones based on vibration measurements are recognized as useful tools because they lead to satisfactory results even when the damage location is not accessible or even unknown [5]. About these techniques, two accepted rules are employed for detecting a crack. The first one is based on the monitoring of the synchronous vibration amplitude and phase. According to [6], changes in 1X amplitude and phase are the primary indicators of crack presence. The second rule relies on 2X vibrations, where [6] states that if a cracked rotor has a steady unidirectional radial load, then a strong 2X response may appear when the rotor is turning at half of any balance resonance speed. However, although widely used in industry, when applied in non-ideal conditions such techniques can detect cracks that eventually have already spread significantly by the cross section of the shaft, usually above 40% of its diameter. Therefore, currently, the researchers' attention is turning to more sophisticated methods capable of identifying incipient cracks (cracks that spread up to 25% of shaft diameter), which represent a type of damage that are hardly observable in vibration analysis [7–9]. Based on previous simulations, [10] and [11] have experimentally detected a crack by applying a specified harmonic force on the rotor by means of a single magnetic bearing. The rotor used in the work is composed of a shaft supported by two active magnetic bearings and one disc located at the shaft middle span. In the analyses, a small notch was cut by using a wire electric discharge machine close to the disc in order to obtain a behavior similar to the breathing crack mechanism. The presence of the damage led to spectral responses with additional peaks at frequencies that are combinations of the rotor speed, its critical speed, and the frequency of the diagnostic force. The method of multiple scales was used to determine the conditions required to create the combinational resonance. In [12], a preliminary experimental application of a crack identification methodology based on combinational resonances was proposed (i.e., technique based on [10,11,13]). The technique uses external applied diagnostic forces at combinational frequencies, together with a pseudo-random optimization code known as Differential Evolution, in order to characterize the signatures of the crack in the spectral responses of the flexible rotor. The rotor used is composed by a horizontal flexible shaft, two rigid discs, and two roller bearings. A small notch was also cut on the shaft in order to obtain breathing crack behavior. The methodology proposed by this work (i.e. preliminary experimental application) proved to be an interesting tool to detect, locate and estimate the depth of transverse cracks in shafts of rotors. Regarding the crack simulation, two models are currently used to represent the breathing behavior with weight dominance, namely the models proposed by Gasch [14] and Mayes and Davies [15]. In both, the mechanism for opening and closing the crack is described by simple mathematical functions. The Gasch's model considers the crack opening and closing abruptly, while the Mayes' model allows for a smooth transition between the fully opened and fully closed crack. It is worth mentioning that when static loads prevail over the dynamic ones (i.e. in a heavy horizontal shaft with relatively low unbalance) the crack opens and closes gradually every revolution of the system. Thus, the breathing assumes weight dominance. Although lead to reliable results, there are other models that can accurately represent the breathing phenomenon. Among them are the 3D model (developed at the research center of EDF) and the FLEX model [7]. However, the computational cost of the Mayes' model is significantly smaller than the FLEX. In [16], an analytical method devoted to the modeling the crack behavior was carried out (i.e., finite shaft element with transverse crack). In [17], the authors present a number of different simplified approaches to modeling cracks in rotating shafts. The influence of the crack models on the dynamic response of a Jeffcott rotor is demonstrated by using the method of harmonic balance. New breathing functions are proposed by [18] in order to formulate the time-varying stiffness matrix of the cracked finite element. The results were compared with typical formulas for representing the breathing mechanism. According to the authors, the proposed functions give more accurate results regarding the dynamic behavior of the cracked rotor. Additionally, [19] proposes a new mathematical model for the cracked rotating shaft based on the rigid finite element method. The crack is represented as a set of spring-damping elements connecting the crack faces. The proposed model is validated by numerical and experimental evaluations. In this context, the present contribution comprises the numerical application of a SHM technique previously proposed by [12]. The open and breathing crack models are compared regarding its dynamic behavior and the efficiency of the proposed identification technique. Based on the nonlinear behavior introduced by transverse cracks in rotating shafts (i.e., the combinational resonance), the crack presence, location, and severity is identified by using an evolutionary optimization method combined with the mathematical models encompassing the rotor and the open crack. Additionally, the dynamic response of the rotating shaft considering the breathing crack is compared with the results obtained from the open crack model (i.e., the small notch used to represent the crack behavior). Therefore, the Fourier expansion of the periodic stiffness matrix truncated at the third harmonic component is performed on the crack finite element. According to [20], cracks where firstly found in shafts of some steam turbines around the 1950s. Thus, research on the vibration behavior of cracked rotating shafts began to be reported in order to prevent major accidents and to develop diagnostic systems. Simultaneously, a few scientific works on nonlinear subharmonic resonances due to ball bearings were published. However, the nonlinear treatment dedicated to cracked rotors is more recent. Research conducted by [10], and other experts, shows an interesting similarity between experimental and numerical results. Thus, it is expected that the methodology proposed here (i.e. Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

3

preliminary experimental application) becomes an interesting tool to detect, locate and estimate the depth of transverse cracks in shafts of rotating machines. The performance of optimization techniques for solving inverse problems of high complexity has been proven in the literature, including the characterization of cracks [21–23].

2. Crack breathing model As mentioned, two models are currently used to represent the breathing behavior with weight dominance, namely the models proposed by Gasch [17] and Mayes and Davies [18]. However, these models are not able to correlate the increased flexibility of a shaft element according to the depth of the crack. Thus, in this work the additional flexibility is calculated by using the linear fracture mechanics theory. This formulation is interesting since it gives the stiffness matrix of the element containing the crack explicitly in terms of its depth. Different formulations able to correlate the flexibility of the shaft element with the crack depth are found in the literature. In [24] a simple mathematical expression is presented. Another approach consists in the cohesive zone model [25], which describes the material failure without considering the geometric parameters. 2.1. Linear fracture mechanics Fig. 1a presents a rotor segment, represented by a shaft element of length L containing a transverse crack of depth α located at the middle position between the nodes 1 and 2. The six-degree-of-freedom element is loaded with axial forces, P1 and P7, shear forces, P2, P3, P8, and P9, torsion moments, P4, and P10, and bending moments, P5, P6, P11, and P12. Details about the cross-section of the element of diameter D (D ¼2R; R is the shaft radius) at the location of the crack are given in Fig. 1b. The cross sectional area without crack is represented by the dashed area (Ac represents the field not hatched). Using the Castigliano's second theorem, the cracked shaft displacement qi in the direction of the load Pi can be determined as shown in Eq. (1), as described by [2]. qi ¼ qUi þqCi ¼

∂U U ∂U C þ ∂P i ∂P i

ð1Þ

where UU is the elastic strain energy of the shaft element under pristine condition and UC is the additional strain energy due to the crack presence; qUi and qCi are displacements associated to the strain energies, respectively. The concepts of fracture mechanics show that the additional strain energy UC is given by the integration of the strain energy density function over the cracked area Ac, as shows Eq. (2), 2 !2 !2 !2 3 Z 6 6 6 X X X 0 4 K li þ K lIi þ ð1 þ νÞ K lIIi 5 dAc ð2Þ UC ¼ E A

c

i¼1

i¼1

i¼1

where E ¼(1  ν)/E and m¼1 þ ν. E is the Young's modulus and ν is the Poisson's ratio. KIi, KIIi and KIIIi are the stress intensity factors (SIF) corresponding to the opening, sliding, and shear modes of the crack displacement, respectively, according to the scheme shown in [26]. The SIF are determined according to the stress distribution on the cross-sectional area of the crack and a shape function, which depend on the size and shape of the body under consideration (i.e. the shaft) [27]. It is important to point out that in 0

Fig. 1. (a) Cracked shaft element. (b) Details about the crack element cross-section.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

practice, the mode I is predominant [7]. The modes II and III demonstrate little influence on the flexibility growth. Eq. (3) shows the SIF for a shaft portion having a crack with depth α.   pffiffiffiffiffiffiffi α ð3Þ K mi ¼ σ i πα F n hX where Kmi is the SIF (m ¼I, II or III), σi is the stress distribution that acts on the crack due to the applied loads Pi (i¼1, 2, …, 6), and Fn are the shape functions (hX is shown in Fig. 1b). The additional flexibility cij introduced in the shaft element due to the crack presence is obtained with the definition of compliance, as given by Eq. (4): cij ¼

∂2 U C ∂P i ∂P j

ð4Þ

where the resulting integrals were calculated by using the procedure described in [28]. Therefore, the additional flexibility matrix due to the crack is: 2 3 0 0 0 0 6 7 0 0 0 7 6 7 cCR ¼ 6 ð5Þ 6 c55 c56 7 4 5 SYM: c66 where the terms located outside the diagonal indicate the coupling between the orthogonal directions [4]. According to [29], acceptable values for the flexibility coefficients are obtained for α/Dr0.8 (due to the accuracy of the shape functions). In this work α/D r0.5 was considered. The axial forces P1 and P7, and the torsion moments, P4 and P10, acting along the Y direction (Fig. 1), are disregarded in the finite element considered for the cracked shaft (Timoshenko's beam element with 4 degrees-of-freedom per node, as shown in Fig. 2). Also, the effect of crack modes II and III are small and disregarded in the analysis. Consequently, the additional flexibility of the shaft due to the crack is represented by using only the coefficients c55, c56, and c66, which are shown in Fig. 3. 2.2. Mayes' breathing model The Mayes' model is based on a cosine function that allows a smooth transition between the fully open and fully closed crack [25]. The stiffness along the η and ξ directions (Fig. 4; Ω is the rotating speed of the rotor and t is the time) according to the Mayes' model (kηMayes and kξMayes, respectively, in the rotating coordinates) are defined by Eq. (6).       kξMayes θ ¼ 12 ko þ kξ þ 12 ko  kξ C 1   1  1  kηMayes θ ¼ ko þkη þ ko  kη C 1 2 2

ð6Þ

where C1 ¼cos θ (Ci ¼cos iθ; i¼1, 2, 3, …), ko is the stiffness of the shaft without crack; kη and kξ are the stiffness of the shaft with crack along the directions η and ξ, respectively. If C1 ¼1, i.e., for θ ¼0°, kξMayes (0°)¼ko and kηMayes (0°) ¼ko. Thus, the crack remains fully closed. However, if C1 ¼  1, i.e., for θ ¼180°, kξMayes (180°) ¼kξ and kηMayes (180°) ¼kη. In this case, the crack remains fully opened [27]. The stiffness of the shaft with crack in rotating coordinates (kRMayes in a matrix representation) according to the Mayes' model is given by: " # kM ξ þ kD ξ C 1 0 kRMayes ¼ ð7Þ 0 kM η þ kD η C 1 where kMξ ¼ (ko þkξ)/2, kMη ¼ (ko þkη)/2, kDξ ¼(ko–kξ)/2, and kDη ¼(ko–kη)/2.

Fig. 2. Timoshenko's shaft element with 4 degrees-of-freedom per node.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 3. Additional flexibility c55 (

), c56 (

), and c66 (

5

) versus crack depth.

Fig. 4. Fixed and rotating coordinates for the shaft with a crack.

Fig. 5. Stiffness behavior in fixed coordinates (

α/D ¼0.20;

α/D ¼ 0.35;

On the other hand, the stiffness of the shaft with crack in determined from a mathematic transformation, as shown in Eq. " # " # " kFMayes ð11Þ C1 S1 T C1 S1 kRMayes ¼ kFMayes ¼ kFMayes ð12Þ  S1 C 1 S1 C 1

α/D¼ 0.50); (a) kFMayes(11); (b) kFMayes(22); (c) kFMayes(12).

fixed coordinates (kFMayes in a matrix representation) is (8). # kFMayes ð12Þ ð8Þ kFMayes ð22Þ

where Si ¼sin iθ (i¼1, 2, 3, …). Thus, the terms of kFMayes are given by:         kFMayes ð11Þ ¼ 12 kM ξ kM η þ 14 3kDξ  kDη C 1 þ 14 kξ  kη C 2  18 kξ kη C 3  1    1 1 1 k þkM η þ kDξ 3kDη C 1  kξ  kη C 2  kξ kη C 3 2 Mξ 4 4 8    1 1 1 k kD η S1 þ kM ξ  kM η S2 þ kD ξ  kD η S3 kFMayes ð12Þ ¼ 4 Dξ 2 4 kFMayes ð22Þ ¼

ð9Þ

According to [29], in fixed coordinates the Mayes' model generates a constant value added to the components 1X, 2X, and 3X, the angular rotation of the rotor; considering the direct stiffness kFMayes(11) and kFMayes(22). The term kFMayes(12) does not have the constant value. Fig. 5 shows the behavior of the stiffness terms described in Eq. (9) considering cracks of different Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

depths located at the node #18 of the finite element model of the rotor (FE model) that will be adopted in the analysis of this work (Fig. 7b). Note that the term kFMayes(12) is well below the direct ones. The stiffness matrix of the shaft element with the transverse crack KCE (a 8  8 matrix in agreement with the element shown in Fig. 2) is determined firstly by adding the flexibility matrix of the shaft without crack co to the additional flexibility matrix due to the crack cCR. A similar procedure is adopted by [30]. cCE ¼ co þ cCR

ð10Þ

where cCE is the flexibility matrix of the element with crack. 1 The stiffness kξ and kη used in Eq. (9) are obtained from the inverse of cCE (kCE ¼cCE ; kξ ¼kCE(11) and kη ¼kCE(22)). Similarly, 1 ko is determined from the inverse of the flexibility matrix co (ko ¼ co ; ko ¼ko(11) ¼ko(22)). Thus, Eq. (9) is resolved and the stiffness kFMayes(11) and kFMayes(22) are substituted into Eq. (11) in order to obtain KCE (IY is the inertia moment with respect to the Y axis and ϑY is a constant that takes into account the shear effect of the Timoshenko's element). The stiffness matrix KCE is clearly a combination of the matrices shown in Eq. (11). 2 3 1 0 2 3" # L 6L 1 L 1 0 1 7 6 7 kFMayes ð11Þ 2 12EI Y 4 5 6 7 KCEXY ¼ L3 1 þ ϑ 6 ð4 þ ϑY Þ 2 ð 0 1 0 1 0 7 YÞ L 41 5 2L 12 0 1 2 3 1 0 2 3" # L 6 L 1  L 1 0 1 7 7 kFMayes ð22Þ  2 12 EI Y 6 5 74 KCEYZ ¼ 3  ð11Þ 6 ð 4 þ ϑY Þ 2 6 0 1 0 1 0 7 L L 1 þ ϑY 4 1 5  2L 12 0 1 The stiffness kFMayes(12), Eq. (9), is disregarded in Eq. (11) since it is much smaller than the stiffness along the X and Z directions (as previously mentioned).

3. Nonlinear technique Fig. 6 shows a flowchart to illustrate the proposed crack identification methodology. The method begins by defining a set of diagnostic forces with different combinational frequencies Ωdiag (i.e. the frequency of the external force applied to the structure) that will induce the rotor to combination resonances. The frequency of each force is determined by using the method of multiple scales. According to [11], the conditions required for a combination resonance occurs for:

ωΩ ¼ 7 Ωdiag þ n Ω

ð12Þ

where ωΩ is a critical speed (forward whirl) of the rotor, Ω is the rotation speed, and n ¼ 71, 72, 73, …. The forward whirl speeds were adopted in this work instead the natural frequencies (a more appropriated choice). The determined diagnostic forces Ωdiag are applied to the analyzed system and the corresponding spectral responses are stored (FFTexp; Fast Fourier Transform). The same forces are applied to a reliable finite element model (Rotor model) of the structure so that the optimizer is responsible for adding the model of a crack with different depths and locations, randomly generated. It is worth mentioning that the unbalance condition of the rotating machine has to be equally applied to the rotor model. The obtained spectral responses (FFTmod) are compared with the experimental ones by means of a given objective function DEOF, as shows

Fig. 6. Crack identification flowchart.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

7

Eq. (13). If the procedure converges to a minimum value of the objective function, the crack is identified. If this is not the case, the optimization procedure will propose a new crack configuration (depth and location). The optimization process continues iteratively until convergence is reached. DEOF ¼

n X ‖FFTexp; i  FFTmod ; i ‖ i¼1

ð13Þ

‖FFTexp; i ‖

where n is the number of FFTs used in the minimization procedure. Thus, the proposed SHM technique is based on crack signatures that appear on the spectral responses of the rotor due to the nonlinear effect caused by the crack presence allied to the application of the diagnostic forces whose frequency is properly determined by the method of multiple scales.

4. Rotor test rig Fig. 7a shows the rotor test rig used in the numerical simulations shown in this work, which is represented by a model with 33 finite elements (FE model; Fig. 7b). It is composed of a flexible steel shaft with 860 mm length and 17 mm of diameter (E¼ 205 GPa, ρ ¼7850 kg/m3, υ ¼0.29), two rigid discs D1 (node #13; 2637 kg; according to the FE model) and D2 (node #23; 2649 kg), both of steel and with 150 mm of diameter and 20 mm of thickness (ρ ¼7850 kg/m3), and two roller bearings (B1 and B2, located at the nodes #4 and #31, respectively). In this case, the vibration responses of the shaft are collect on the orthogonal directions of the nodes #8 (S8X and S8Z) and #28 (S28X and S28Z). Eq. 14 presents the equation of motion that governs the dynamic behavior of the cracked flexible rotor supported by roller bearings [31].   ð14Þ M q€ þ Dþ Ω Dg q_ þ KðΩtÞq ¼ W þ Fu þ Fe where M is the mass matrix, D is the damping matrix, Dg is the gyroscopic matrix, and K(Ωt) is the stiffness matrix with variable values due to the crack (i.e., Ωt stands to the angular position of the shaft). All these matrices are related to the rotating parts of the system, such as couplings, discs, and shaft. W stands for the weight of the rotating parts, Fu represents the rotating unbalance forces, and Fe represents the diagnostic force applied in the rotor (force fixed in space), and q is the generalized displacement vector (see Fig. 2).

kX/B1 = 8.551 x 105 kZ/B1 = 1.198 x 106 dX/B1 = 7.452 dZ/B1 = 33.679

kX/B2 = 5.202 x 107 kZ/B2 = 7.023 x 108 dX/B2 = 25.587 dZ/B2 = 91.033

γ = 2.730 β = 4.85 x 10-6 kROT = 770.442

Fig. 7. Rotor used in the numerical simulations of the SHM technique; (a) test rig; (b) FE model. *k: stiffness [N/m]; d: damping [N s/m].

Fig. 8. Updated model (

) and experimental (

) FRFs of the rotor; (impact along the X direction of D1 and sensor S8X).

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

A model updating procedure was used in order to obtain a representative FE model for the numerical simulations presented on this work (Fig. 7b). In this sense, a heuristic optimization technique (Differential Evolution [32]) was used to determine the unknown parameters of the model, namely the stiffness and damping coefficients of the bearings, the proportional damping added to D (coefficients γ and β; Dp ¼ γMþ βK), and the angular stiffness kROT due to the coupling between the electric motor and the shaft (added around the orthogonal directions X and Z of the node #1). The entire frequency domain process was performed 10 times, considering 100 individuals in the initial population of the optimizer. The objective function adopted is similar to the one presented in Eq. (13). However, in this case only the regions close to the peaks associated with the natural frequencies were taken into account. The experimental FRF (frequency response functions) were measured on the test rig for the rotor at rest by applying impact forces along the X and Z directions of both discs, separately. The response signals were obtained by the two proximity probes positioned along the same directions of the impact forces, resulting 8 FRF. The measurements were performed by the analyzer Agilents (model 35670A) in a range of 0– 200 Hz and steps of 0.25 Hz. Fig. 7 summarizes the parameters determined in the end of the minimization process associated with the smaller fitness value. Fig. 8a compares one simulated and experimental FRF considering the parameters shown in Fig. 7; validating the updating procedure performed. Fig. 8b compares the experimental orbit measured by the measure plan S8 (full acquisition time of 4 s in steps of 0.002 s) with the one determined by the updated FE model. The operational spin speed of the rotor Ω was fixed to 1200 rev/min and an unbalance of 487.5 g mm/0° applied to the disc D1 was considered (Fig. 9).

5. Numerical application The crack identification methodology was tested for the rotor under three different structural conditions. The first one comprises the shaft without crack (pristine condition). The second test was performed for the shaft with a crack located at the element #10 (between the nodes #10 and #11) with 20% depth. The last test was performed for the shaft with a crack located at the element #20 (between the nodes #20 and #21) with 20% depth. Simulations verified that the results were not sensitive to the length of the crack element of 45 mm. Additional numerical evaluations can be seem in [33,34]. For all the three structural conditions treated by this contribution, an unbalance of 637.5 g mm/  90° applied to the disc D1 was considered. Additionally, in all analyzes the operational spin speed of the system Ω was fixed to 1200 rev/min. The frequencies of the diagnostic forces were chosen in order to be distributed along the frequency bandwidth of the four first natural frequencies of the rotor. Thus, considering forward whirl critical speeds, for ωΩ ¼ ωfw1 ¼28.57 Hz, the diagnostic frequencies are the following: Ωdiag ¼48.57, 8.57, 68.57, 11.43, 51.43, and 71.43 Hz. They were obtained, respectively, for n ¼ 1, 1,  2, 2, 4, and 5. Using now ωΩ ¼ ωfw2 ¼97.66 Hz, the diagnostic frequencies are the following: Ωdiag ¼77.66, 57.66, 37.66, 17.66, 22.34, and 42.34 Hz. They were obtained, respectively, for n ¼1, 2, 3, 4, 6, and 7. This procedure gives a total of 12 diagnostic frequencies, resulting 24 FFTs responses, as only the two sensors along the horizontal direction were used (Ωdiag close to the critical speeds, Ωdiag 4110 Hz, and Ωdiag o5 Hz were discarded). All the diagnostic forces have the same amplitude. Tests were performed for different amplitudes and it was possible to observe that the amplitude has to be changed according to the level of noise. Thus, for the results that will be presented in the following, a satisfactory amplitude value for the diagnostic forces was found to be 25 N (forces applied along the Z direction in the node #4 of the FE model; see Fig. 7b). Fig. 10a shows the FFTs obtained using Ωdiag ¼ 48.57 Hz (sensor S8Z) for the shaft in its pristine condition, then for the case with a 15% depth breathing crack located at the element #20, and finally for a 30% depth breathing crack located at the same element. Note that the amplitudes of the peaks related with each rotor condition are different (i.e., the peaks associated with the combination frequencies: Ωdiag–2Ω, –Ωdiag þ3Ω,  Ωdiag þ4Ω, and so on). Fig. 10b shows the FFTs obtained using the same diagnostic force (sensor S8Z), but for the shaft with a 30% depth breathing crack at the element #5 and a 30% depth breathing crack at the element #20. Observe that the peaks associated with the combination resonances also present different amplitudes. For others Ωdiag, similar behavior is noted. As mentioned, the proposed methodology uses a set of

Fig. 9. Updated model (

) and experimental (

) orbits of the rotor determined at the measure plan S8.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Ω

– Ωdiag + 4Ω

4Ω

5Ω

Ωdiag 2Ω

3Ω

Ωdiag + 2Ω

3Ω

Ωdiag – 2Ω –Ωdiag + 3Ω

2Ω

ωfw1 Ωdiag + 3Ω

Ωdiag

Ωdiag + 2Ω

ωfw1

Ωdiag + Ω

–Ωdiag + 3Ω

Ωdiag – 2Ω

Ω

4Ω

Ωdiag + 3Ω

b

Ωdiag + Ω

a

9

5Ω

– Ωdiag + 4Ω

Fig. 10. FFTs obtained using Ωdiag ¼ 48.57 Hz for the shaft in two structural conditions ( □ pristine condition); (a) 15% depth breathing crack at element #20; 30% depth breathing at element #20; (b) 30% depth breathing crack at element #5; 30% depth breathing at element #20.

diagnostic forces with different frequencies to identify cracks. It is possible due to the behavior evidenced by Fig. 10. The application of several diagnostic forces with different frequencies makes the spectral signature of the crack specific for its characteristic (i.e. location and depth), which allows the correct identification. As mentioned, when static loads prevail over the dynamic ones the crack opens and closes gradually every revolution of the system. Thus, the breathing assumes weight dominance [7]. The crack behavior becomes function of the shaft angular position Ωt only, lead to a periodical stiffness variation. Therefore, the stiffness matrix of the shaft K(Ωt) (see Eq. (14)) can be written according to Eq. (15). KðΩtÞ ¼ Km þ

3 X

ΔKj ei jΩt

ð15Þ

j¼1

where Km is considered the mean stiffness matrix of the shaft with the crack and ΔKj (j ¼1, 2, and 3) are the stiffness variation related to the components 1X, 2X, and 3X of the rotor speed. The Fourier expansion of the periodic stiffness is truncated at the third harmonic component. Fig. 11 shows the behavior of the stiffness terms described in Eq. (15) considering a 20% depth crack located at the element #20 of the rotor FE model. Both open and breathing crack models are presented regarding the direct stiffnesses kFMayes(11) and kFMayes(22) (refer to Fig. 5). Note that the open crack behavior is represented only by the mean stiffness matrix and the component ΔK2. The breathing crack is affected by all the expansion terms considered in Eq. (15). Fig. 12 shows the FFTs obtained using Ωdiag ¼48.57 Hz (Fig. 12a – sensor S8X; Fig. 12b – sensor S8Z) for the shaft with a 20% depth crack located at the element #20, considering both the breathing and open models. Note that the amplitudes of the peaks related with each crack model are different. The peaks associated with the frequencies Ωdiag–2Ω, –Ωdiag þ3Ω, ωfw1, – Ωdiag þ4Ω, and Ωdiag þ3Ω, are not so evident for the case of open crack when the sensor S8X is considered. Similar behavior can be observed for the harmonic components 2Ω and 4Ω, which is not the case for the 3Ω and 5Ω vibration peaks. Regarding the sensor S8Z, the peaks associated with the frequencies –Ωdiag þ3Ω, ωfw1, Ωdiag þ Ω, and Ωdiag þ3Ω, are not evident for the case of open crack. For others Ωdiag, similar behavior is noted. However, the peaks associated with the combination resonances are still present on the spectral response considering the shaft with an open crack (e.g., the combination resonance Ωdiag þ2Ω in Fig. 12a and b), which allows the correct crack identification. It is worth mentioning that the presence of the peaks associated with the combination frequencies is interesting, since the amplitudes of the harmonic components Ω, 2Ω, 3Ω, and so on, are frequently affected by unbalance, misalignment, and other machinery malfunctions [34]. Therefore, the FE model of the rotor test rig shown in Fig. 7 considering the adjusted parameters is used in the numerical applications of the proposed SHM technique. The open crack model is adopted since many contributions try to represent the crack behavior by machining a small notch on the shaft instead by using the fatigue process. Here, only the proximity probes placed only along the vertical direction of the shaft is considered (the sensors S8Z and S28Z). Random noise is added to the dynamic responses of the rotating machine in order to simulate the experimental environment found in a real plant. Eq. (16) indicates how noise was added to the responses. Remember that the diagnostic forces were applied along the Z direction in the node #4 of the FE model (25 N). n o1=2 qnoise ¼ q þ P noise R noise E½q  EðqÞ2

ð16Þ

where qnoise is the response with noise, Pnoise is the parameter that defines the amount of noise to be added (Pnoise ¼1% in the present case), and Rnoise is the random noise. E [  ] indicates the expected of the value [  ]. Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

Fig. 11. Stiffness behavior in fixed coordinates ( K(Ωt); K m; ΔK1; ΔK2; ΔK3); (a) breathing crack model – kFMayes(11); (b) open crack model – kFMayes(11); (c) breathing crack model – kFMayes(22); (d) open crack model – kFMayes(22).

4Ω

5Ω

ωfw1 2Ω

Ωdiag 3Ω

Ωdiag + 2Ω

3Ω

–Ωdiag + 3Ω

2Ω

Ωdiag – 2Ω

ωfw1

Ωdiag + 2Ω

Ω

Ωdiag

Ωdiag – 2Ω –Ωdiag + 3Ω

Ω

4Ω

5Ω

Ωdiag + 3Ω

b

Ωdiag + Ω

a

– Ωdiag + 4Ω

Fig. 12. FFTs obtained using Ωdiag ¼ 48.57 Hz for the shaft with a 20% depth crack at element #20;( (b) sensor S8Z.

breathing crack;

open crack); (a) sensor S8X;

Table 1 shows the results obtained in the crack identification processes, driven by the optimizer (i.e. Differential Evolution). Twenty individuals were used in the initial population to identify the structural condition of the rotor. Five optimization processes were carried out for each analysis in order to avoid local minima solutions. It is important to point out Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

11

that Table 1 presents only the best optimization process (i.e., associated with the smallest fitness value). The element with crack and the crack depth were used as design variables, being the candidates shown in Table 2 (reduced possibilities to provide a faster minimization process). Note that the rotor conditions were satisfactorily identified. In [33] different rotor conditions were successfully identified considering the breathing crack (i.e., the Mayes' model). The results obtained lead to the numerical validation of the proposed crack identification technique. A preliminary experimental application of the proposed methodology is presented in [12]. Fig. 13 compares FFTs obtained from the rotating machine with those found in the end of the identification process associated with the smallest fitness value for the pristine rotor (see Table 1) using the sensor S8Z (i.e., node #8 along the Z Table 1 Results of the crack identification process. Structure configuration

Smallest fitness value

Results Element

Depth [%]

Pristine condition Element #10 depth 20%

0.0146 0.0145

#31 #10

0 20

Element #20 depth 20%

0.0148

#20

20

Table 2 Possible elements and depths used. Elements with crack

Crack depth [α/D]

Elements #5 to #11, #14 to #21, and #24 to #33.

0, 2.5%, 5%, 7.5%, 10%, 12.5%, 15%, …, and 50%.

Fig. 13. Crack identification process; pristine rotor ( (d) Ωdiag ¼77.66 Hz.

rotor model;

identified). (a) Ωdiag ¼48.57 Hz; (b) Ωdiag ¼ 11.43 Hz; (c) Ωdiag ¼ 17.66 Hz;

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

a

b Ωdiag Ω

2Ω

3Ω 4Ω

Ωdiag + 2Ω

Ωdiag

2Ω 5Ω

Ωdiag + 2Ω

Ω

3Ω 4Ω

5Ω

–Ωdiag + 2Ω

c

d Ω

2Ω

3Ω

4Ω

Ωdiag + 4Ω

Ω Ωdiag + 2Ω

Ωdiag

Ωdiag

2Ω 5Ω

Fig. 14. Crack identification process; element #10 with 20% depth ( (c) Ωdiag ¼ 17.66 Hz; (d) Ωdiag ¼ 77.66 Hz.

rotor model;

3Ω

5Ω

identified).(a) Ωdiag ¼ 48.57 Hz; (b) Ωdiag ¼11.43 Hz;

direction). The goal here was simply to check if the proposed methodology does not lead to false positive alarm. In these FFTs, one can observe only the peaks associated to the excitation frequencies and diagnostic forces, indicating the healthy state of the beam (i.e., no extra peaks are found in the spectra). Note that the identification process reproduced satisfactorily the FFTs obtained from the FE model of the rotating machine and is able to indicate correctly the pristine condition. The remaining FFTs are similar to the previous ones. Fig. 14 shows the results of the identification process associated with the smallest fitness value when the rotor was affected by a crack located in the element #10 with 20% depth, while Fig. 15 presents the results of the identification process associated with the smallest fitness value regarding a crack located in the element #20 with 20% depth. Now, it is possible to observe the diagnostic peaks. Note that the peaks are found for different amplitudes, depending on the frequency of the diagnostic force. The diagnostic peaks show up due to the fact that a deeper crack changes more intensely the dynamic behavior of the beam under the application of the diagnostic forces, i.e., the nonlinear effect becomes more evident in the system responses. The remaining FFTs are similar to the previous ones.

6. Conclusion The presented numerical results show that the developed inverse problem methodology is a promising tool for estimating satisfactorily the existence, location, and depth of a transverse crack found in shafts of rotating machines. However, a representative mathematical model of the mechanical system and the crack are required (i.e., the open model is also based on the linear fracture mechanics). Additionally, it was shown that the diagnostic peaks observed in the FFTs depend on the breathing behavior of the crack. The proposed methodology shows that the diagnostic forces generate diagnostic peaks in the region of the lower vibration modes, where the frequency bandwidth is associated with high vibration energy. Thus, the identification of the crack becomes easier. It is necessary to apply several diagnostic forces to the system to identify the crack. Thus, a time-consuming procedure is expected for a future experimental testing that will be performed. However, in certain cases the diagnostic forces can be applied to the structure under operating conditions. The amplitude of the diagnostic forces can be adjusted to keep the system on a safe vibration level. Additionally, as the amplitude of the diagnostic Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

a

13

b

2Ω

Ωdiag – 2Ω

Ωdiag Ω

Ωdiag + 2Ω

Ωdiag 3Ω 4Ω

2Ω Ωdiag + 2Ω

Ω

5Ω

3Ω 4Ω

5Ω

–Ωdiag + 2Ω

c

d Ω

Ω

4Ω

5Ω

Fig. 15. Crack identification process; element #20 with 20% depth ( (c) Ωdiag ¼ 17.66 Hz; (d) Ωdiag ¼ 77.66 Hz.

Ωdiag Ωdiag – 2Ω

3Ω

Ωdiag + 4Ω

2Ω Ωdiag + 2Ω

Ωdiag

rotor model;

2Ω 3Ω

5Ω

identified). (a) Ωdiag ¼48.57 Hz; (b) Ωdiag ¼11.43 Hz;

forces can be adjusted so that difficulties that arise from a faster crack growth due to the effect of adding additional forces to the system can be minimized. The results presented demonstrate the efficiency of the methodology conveyed. The positive results conveyed encourage the authors on doing further research efforts dedicated to experimental tests for validation purposes.

Acknowledgments The authors are thankful to the Brazilian Research Agencies FAPEMIG and CNPq (INCT-EIE / 5740012008-5) and also to CAPES for the financial support provided for this research effort.

References [1] S.W. Doebling, C.R. Farrar, M.B. Prime, A summary review of vibration-based damage identification methods, Shock Vib. Dig. 30 (2) (1998) 91–105. [2] A.K. Darpe, K. Gupta, A. Chawla, Coupled bending, longitudinal and torsional vibrations of a cracked rotor, J. Sound Vib. 269 (1) (2004) 33–60. [3] J.T. Sawicki, M.I. Friswell, A.H. Pesh, A. Wroblewski, Condition monitoring of rotor using active magnetic actuator, in: Proceedings of ASME Turbo Expo 2008: Power for Land, Sea and Air, Berlin, Germany, June 9–13, 2008. [4] A.D. Dimarogonas, Vibration of cracked structures: a state of the art review, Eng. Fract. Mech. 55 (5) (1996) 831–857. [5] E.P. Carden, P. Fanning, Vibration based condition monitoring: a review, Struct. Health Monit. 3 (5) (2004) 355–377. [6] D.E. Bently, C.T. Hatch, Fundamentals of Rotating Machinery Diagnostics, Bently Pressurized Bearing Company, Minden, NV, USA, 2002. [7] N. Bachschmid, P. Pennacchi, E. Tanzi, Cracked Rotors: A Survey on Static and Dynamic Behaviour Including Modelling and Diagnosis, Springer-Verlag, Berlin, Germany, 2010. [8] A.S. Sekhar, Crack identification in a rotor system: a model-based approach, J. Sound Vib. 270 (1) (2003) 887–902. [9] Z. Kulesza, J.T. Sawicki, Auxiliary state variables for rotor crack detection, J. Vib. Control 17 (6) (2010) 857–872. [10] J.T. Sawicki, D.L. Storozhev, J.D. Lekki, Exploration of NDE properties of AMB supported rotors for structural damage detection, J. Eng. Gas Turbines Power 133 (1) (2011) 1–9. [11] G. Mani, D.D. Quinn, M. Kasarda, Active health monitoring in a rotating cracked shaft using active magnetic bearings as force actuators, J. Sound Vib. 294 (1) (2006) 454–465.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i

14

A.A. Cavalini Jr et al. / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

[12] A.A. Cavalini Jr, L. Sanches, V. Steffen Jr, Experimental analysis of a model based crack identification approach for rotating machines, in: Proceedings of the ASME 2014 International Design Engineering Technical Conferences, Buffalo, NY, USA, 2014. [13] J.E.T. Penny, M.A. Friswell, The dynamics of cracked rotors, in: Proceedings of the IMAC XXV Orlando Florida, USA, 2007. [14] R. Gasch, Dynamic behaviour of a simple rotor with a cross sectional crack, in: Proceedings of the IMechE International Conference on Vibrations in Rotating Machinery, Cambridge, UK, 1976. [15] I.W. Mayes, W.G.R. Davies, Analysis of the response of multi-rotor-bearing system containing a transverse crack in a rotor, J. Vib. Acoust. Stress Reliab. Des. 106 (1) (1984) 139–145. [16] H.D. Nelson, C. Nataraj, The dynamics of a rotor system with a cracked shaft, J. Vib. Acoust. 108 (2) (1986) 189–196. [17] J.E.T. Penny, M.I. Friswell, Simplified modeling of rotor cracks, in: Proceedings of the ISMA 27, Leuven, Belgium, 2002. [18] M.A. Al-Shudeifat, E.A. Butcher, New breathing functions for the transverse breathing crack of the cracked rotor system: approach for critical and subcritical harmonic analysis, J. Sound Vib. 330 (1) (2011) 526–544. [19] Z. Kulesza, J.T. Sawicki, Rigid finite element model of a cracked rotor, J. Vib. Control 331 (18) (2012) 4145–4169. [20] Y. Ishida, T. Yamamoto, Linear and Nonlinear Rotordynamics, Wiley-VCH, Weinheim, Germany, 2012. [21] Y. He, D. Guo, F. Chu, Using genetic algorithms and finite element methods to detect shaft crack for rotor-bearing system, Math. Comput. Simul. 57 (1) (2001) 95–108. [22] A.S. Sekhar, Detection and monitoring of cracks in a coast-down rotor supported on fluid film bearings, Tribol. Int. 37 (3) (2004) 279–287. [23] P. Pennacchi, N. Bachschmid, A. Vania, A model-based identification method of transverse cracks in rotating shafts suitable for industrial machines, Mech. Syst. Signal Process. 1 (1) (2006) 2112–2147. [24] I.W. Mayes, W.G.R. Davies, The vibration behavior of a rotating shaft system containing a transverse crack, in: Proceedings of the IMechE International Conference on Vibrations in Rotating Machinery, Cambridge, UK, 1976. [25] R.T. Liong, Application of the cohesive zone model to the analysis of a rotor with a transverse crack, in: Proceedings of the 8th IFToMM International Conference on Rotor Dynamics, Seoul, Korea, 2010. [26] T.L. Anderson, Fracture Mechanics: Fundamentals and Applications, Taylor & Francis, USA, 2005. [27] C.R. Burbano, V. Steffen Jr, Diagnosis of cracked shafts by monitoring the transient motion response, in: Proceedings of the DINAME XII International Symposium on Dynamic Problems of Mechanics, Ilhabela, Brazil, 2007. [28] C.A. Papadopoulos, Some comments on the calculation of the local flexibility of cracked shafts, J. Sound Vib. 278 (1) (2004) 1205–1211. [29] M.I. Friswell, J.E.T. Penny, Crack modeling for structure health monitoring, Struct. Health Monit. 1 (2) (2002) 139–148. [30] A.S. Sekhar, Multiple cracks effects and identification, Mech. Syst. Signal Process. 22 (1) (2008) 845–878. [31] M. Lalanne, G. Ferraris, Rotordynamics Prediction in Engineering, John Wiley & Sons, INC., Baffins Lane, Chichester, West Sussex PO19 1UD, England, 1998. [32] R. Storn, K. Price, Differential evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces, Int. Comput. Sci. Inst. 12 (1) (1995) 1–16. [33] A.A. Cavalini Jr, V. Steffen Jr, J. Mahfoud, Crack identification approach for rotating machines based on combinational resonances, in: Proceedings of the 22nd COBEM International Congress of Mechanical Engineering, Ribeirão Preto, Brazil, 2013. [34] F.F. Ehrich, Handbook of Rotordynamics, Mcgraw-Hill, New York, 1991.

Please cite this article as: A.A. Cavalini Jr, et al., Crack identification for rotating machines based on a nonlinear approach, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.041i