Experimental wear parameters identification in hydrodynamic bearings via model based methodology

Author’s Accepted Manuscript Experimental Wear Parameters Identification in Hydrodynamic Bearings via Model Based Methodology Ricardo U. Mendes, Tiago H. Machado, Katia L. Cavalca www.elsevier.com/locate/wear

PII: DOI: Reference:

S0043-1648(16)30735-9 http://dx.doi.org/10.1016/j.wear.2016.12.002 WEA101843

To appear in: Wear Received date: 14 June 2016 Revised date: 18 October 2016 Accepted date: 5 December 2016 Cite this article as: Ricardo U. Mendes, Tiago H. Machado and Katia L. Cavalca, Experimental Wear Parameters Identification in Hydrodynamic Bearings via Model Based Methodology, Wear, http://dx.doi.org/10.1016/j.wear.2016.12.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

Experimental Wear Parameters Identification in Hydrodynamic Bearings via Model Based Methodology Ricardo U. MENDES1, Tiago H. MACHADO1*, Katia L. CAVALCA1 1

Laboratory of Rotating Machinery, Faculty of Mechanic Engineering, UNICAMP; 200, Rua Mendeleyev, Postal Box 6122, Campinas, CEP: 13083-970, SP, Brazil. * Corresponding author: [email protected] Tel.: +55 19 3521 3178 Abstract When a rotor operates in very low rotating speeds, the fluid-film in hydrodynamic bearings is not completely formed, which allows direct metal to metal contact between the bearing and the journal. Numerous start/stop cycles also may lead to bearing wear. In this context, the present work aims to evaluate a model-based method for identification of hydrodynamic bearing wear parameters: more specifically, wear depth and angular position. The proposed method is based on minimizing an objective function, which compares the directional Frequency Response Function (dFRF) of the physical system with the response of the model developed here. It contemplates the numerical model of the rotor and the worn bearings, in this case, represented by linearized coefficients of damping and stiffness. Through experimental results, a deeper evaluation of the identification method is presented, such as the sensitivity to the number of points of the dFRF used in the objective function, different meshes of starting points, and the region of the dFRF used – around the resonance peak or around the oil-whirl peak. The analyses are conducted for different levels of bearing wear. The precision obtained in the results indicates that the presented method is a promising tool in the identification of bearing wear. Keywords: bearings wear, parameters identification, directional frequency response. 1. Introduction Rotating machines are used in a wide range of industrial processes. To avoid sudden failures during operation that can seriously affect the productivity, early fault detection is essential. It is crucial for effective condition monitoring and to complement the predictive maintenance policies, that are known to reduce the cost considerably. Vibration analysis is widely used as an effective tool for fault monitoring and diagnosis of many ill conditions in rotating machinery, such as misalignment, cracks, shaft bow and bearing faults. In this context, the vibration analysis or modal analysis has always gained attention in the industry and research community. Rotating machines present very particular characteristics, such as non-symmetric stiffness and damping matrix due to gyroscopic effects and anisotropy of journal bearings and supporting structures. The implication of such peculiarities is that classical modal analysis methods, in most cases, are not suitable for the modal parameters identification of rotating systems. In 1975, Gasch and Pfützner [1] showed that when a rotating system is represented by the use of complex coordinates, the rotor precession movement can be

2 described in terms of two rotating vectors, one in the same direction of the angular velocity of the shaft, called forward, and other in the opposite direction, called backward. Nordmann [2] introduced some improvements in the classical modal analysis, and therefore, modal analysis for rotating machines could be treated in an appropriate way. The main improvement was the investigation of the modal parameters, such as eigenvalues and natural modes of a simple shaft, considering some effects peculiar to rotating systems (non symmetry of matrices and speed dependence of modal parameters, for instance). A combined experimental and analytical method was used to identify these parameters and a sensitivity analysis was applied to study the influence of the parameters to the dynamics of rotors. In 1991, Lee [3] presented the complex modal analysis for rotating machines used to distinguish between forward and backward modes shapes, and compared it with the classical modal analysis. With the use of complex notation, not only the forward and backward modes were easily distinguished, but were also separated in frequency domain, facilitating the identification of the modal parameters, and enabling the analysis of the degree of anisotropy of the system. Later, Kessler [4] used complex notation to represent directional vectors (forward and backward) and combined them to represent the dynamic of rotating machineries and the planar forces. Therefore, an association of a direction (forward or backward) to a frequency response function (directional FRF) was possible. In the same field, Mesquita et al. [5] presented a comparison between the use of the traditional Frequency Response Function and the directional Frequency Response Function in the rotordynamics analysis using simulations of rotors supported by isotropic and anisotropic bearings. Latter, Cavalca et al. [6] studied the foundation effect on the rotor behavior inside the concept of mixed coordinates. The degrees of freedom related to the most significant foundation modes were transformed to principal coordinates, allowing the direct representation of these modes. The global rotorbearings foundation system, hence, was represented by an equation of motion with mixed coordinates: physical coordinates for the rotor and principal coordinates for the foundation. With the use of directional notation, Cavalca and Okabe [7] continued the work and evaluated the influence of the supporting structure on the directional response of a rotor-bearing-foundation system. Still in the context of rotating machines, journal bearings play an important role in the dynamics of rotor-bearings systems, being responsible for supporting the machine and transmit forces between the rotor and the foundation. For this reason, in order to perform a dynamic analysis of the complete rotating system, it is necessary to know the dynamic characteristics of bearings, understanding the main phenomena related to this component, allowing predictions of malfunctions or even a proper predictive maintenance. Consequently, evaluation and monitoring of bearings conditions are critical issues in rotor dynamics. Under this perspective, the wear of the bearing material affects the bearing clearance, which results in changes in the dynamic characteristics of the bearing and, consequently, of the rotating system. Katzenmeier [8] was one of the first authors to analyze the effects of wear in hydrodynamic bearings, analyzing quantitatively and qualitatively several examples of damage due to wear in bearings after long periods of operation of rotating machines. Duckworth and Forrester [9] and Mokhtar et al. [10] conducted the first experimental studies related to the wear in hydrodynamic bearings. They investigated the evolution of wear in hydrodynamic bearings over several start/stop cycles, evaluating changes in

3 diametrical clearance, surface finish and roundness of the bearing. The authors concluded that the wear caused during the rotor starting, and the stopping process does not significantly contribute to wear. From experiments conducted on a steam turbine, Dufrane et al. [11] proposed two geometrical models for the bearing worn region, one based on the concept that the shaft makes a 'print' in the bearing and other based on the assumption of an abrasive wear. Hashimoto et al. [12] investigated theoretically and experimentally the effects of geometric changes due to wear, validating the second model proposed by Dufrane et al. [11]. Latter in the nineties, Scharrer et al. [13] presented a research on the dynamic coefficients of journal bearings, showing that small wear depths have small influence on bearing performance. In the same decade, Suzuki and Tanaka [14] analyzed the stability of a worn bearing. They concluded that the wear decreases the stability of the bearing when it is submitted to a light load. The stability of worn bearings was also studied by Kumar and Mishra [15] and [16] in non-laminar lubrication regimes. Regarding the most recent works, Fillon and Bouyer [17] studied the performance of worn plain journal bearings, taking into account the local thermal effects. Later, Bouyer et al. [18] presented and discussed some experimental data about the influence of wear on the behavior of a journal lobed bearing subjected to numerous starts and stops. Finally, Muzakkir et al. [19] experimentally investigated the effectiveness of axial and circumferential grooves in minimizing the wear of journal bearing operating in mixed lubrication regime. Considerable research has also been carried out for the development of various techniques for bearing fault detection and diagnosis. As described by Machado et al. [20], these techniques can be mainly classified into two categories: time domain and frequency domain techniques. For the time domain techniques, Wu et al. [21] proposed an experimental system to monitor the friction conditions in hydrodynamic bearings called "On-Line Visual Ferrograph". With the use of this system, it was concluded that the predominant wear mechanisms were the 'micro-plowing' and 'micro-cutting' induced by the initial roughness of the surfaces in the startup of the rotor. Gertzos et al. [22] developed a practical methodology to identify the wear depth in hydrodynamic bearings, using the model of Dufrane et al. [11]. The method was based on measurements of the basic bearing characteristics to obtain real time (online) wear identification. Later, Okac et al. [23] developed a new scheme based on wavelet packet decomposition and hidden Markov modeling (HMM) for tracking the severity of bearing faults, such as wear. Finally, Chasalevris et al. [24] presented an investigation on the emergence of additional harmonic components in the transient response of a continuous rotor mounted on worn hydrodynamic bearings. The harmonics were more sensitive to wear especially on the frequency of 1/2X and during the crossover through resonance. In the frequency domain techniques context, Papadopoulos et al. [25] presented a theoretical identification method for the wear in hydrodynamic bearings by means of measurements of the rotor response at a given point (usually the rotor midpoint). Least square technique was used to the identification at the predefined point. Machado and Cavalca [26] presented a numerical approach to identify the bearing wear parameters using the unbalance response of the system. The analysis took into account the frequency response due to unbalance of the rotor-bearing system in directional

4 coordinates. A searching method was used to minimize an objective function, considering the rising of the backward component due to the increasing of bearing anisotropy degree. Finally, Machado et al. [20] numerically applied the concept of wear identification, using the directional Frequency Response Function (dFRF). The identification procedure was tested using a simulated experimental response obtained by a noise formulation applied in the numerical response. The results showed that the dFRF matrix is sensitive to the wear parameters and more suitable for bearing wear identification, based on the cross-coupled terms of the dFRF matrix, which presented high sensitivity to bearing wear. The method proposed by Machado et al. [20] overcomes the problem of operating the rotor near the critical speed, from the method presented in the work of Machado and Cavalca [26]. However, in order to measure the dFRF, it is necessary to have an external excitation source, such as magnetic bearing or magnetic actuator (as suggested by the authors). It is important to highlight that many authors (Baumann, [27]; Moore et al. [28]; Sorokes et al. [29]; Bidaut et al. [30]; Pettinato et al. [31]) have already shown that this device can be used as an external excitation source for stability measurements in commissioning tests; replacing the well-known pedestal mounted electromechanical exciters. Usually, the magnetic bearing is temporarily mounted on the free end of the compressor’s shaft. In this scenario, the paper gives continuity to the analyses of the influence of wear in hydrodynamic bearings and the identification of the wear parameters using the directional frequency response of a rotating system. In fact, it is an experimental validation of the method proposed by Machado et al. [20]. Initially, the numerical model of the journal bearing supported rotor is described, followed by a brief introduction about directional analysis of rotors and a wear parameters identification method. The test rig used to validate the method is presented and, as an initial part of the identification procedure, its numerical model is adjusted. Next, a set of 5 different configurations (i.e. 5 different bearings – a healthy bearing and four bearings with different levels of wear) were used for the sensitivity analyses. The sensitivity analyses begin by comparing the wear parameter identification performance when using the Frequency Response Function in physical coordinates versus in directional coordinates; and the use of the four terms of the frequency response versus the use of only the cross terms. In following analysis, the use of points of the dFRF around the natural frequency peak and around the oil-whirl peak were compared, and the sensitivity of the method regarding the number of points of the dFRF used in the objective function was also analyzed. Finally, the sensitivity of the identification to the mesh of starting points used by the search algorithm is investigated. The results are presented by a set of tables; and the reconstructed frequency response for a selected case is graphically presented along with the experimental measurement. The main contribution of the work lies on the description of the wear parameters identification method, and of the details regarding its experimental application. The results have shown that the method is a promising tool for the identification of bearing wear, presenting a valuable contribution to the fault identification analysis in rotating systems.

5 2. Methodology This section presents how the complete system is modeled, including the rotor and hydrodynamic bearings, followed by a brief description of the theory behind the conversion of the FRF from physical to directional coordinates. The optimization problem to be solved is then addressed, and the full identification procedure is described. Finally, the test rig used to validate the identification method is presented. 2.1. System model The model of the rotor was obtained using the well-known finite element method (FEM) as described in the work of Nelson [32], which considers the Timoshenko beam theory. The final equation of motion of the system is given by Eq. 1, where [M] is the mass matrix containing translational and rotational inertias, [K] is the stiffness matrix containing both the shaft stiffness and the added stiffness of the bearings, [G] is gyroscopic matrix and ȳ is the rotational speed. The damping matrix [C] contains the proportional structural damping matrix [Cs] regarding the shaft, Eq. 2, and the damping contribution of the bearings. The system model has two translational and two rotational degrees of freedom per node in the generalized coordinate vector {q}, and the external forces due to the unbalance and the magnetic actuator are inserted in the force vector {F}. ሾ‫ܯ‬ሿሼ‫ݍ‬ሷሽ ൅ ሺሾ‫ܥ‬ሿ ൅ ȳሾ‫ܩ‬ሿሻሼ‫ݍ‬ሶሽ ൅ ሾ‫ܭ‬ሿሼ‫ݍ‬ሽ ൌ ሼ‫ܨ‬ሽ

(1)

ሾ‫ܥ‬௦ ሿ ൌ ߙሾ‫ܯ‬௦ ሿ ൅ ߚሾ‫ܭ‬௦ ሿ

(2)

ሾ‫ܪ‬ሿሼ‫ݍ‬ሽ ൌ ሼ‫ܨ‬ሽ  ื  ሼ‫ݍ‬ሽ ൌ ሾ‫ܪ‬ሿିଵ ሼ‫ܨ‬ሽ

(3)

Assuming that the rotating system described by Eq. 1 is excited by a harmonic force {F} of frequency ω and presents a complex response {q}, the FRF of the system ([H]-1 matrix) can be obtained as presented in Eqs. 3 and 4.

ሾ‫ܪ‬ሿିଵ ൌ ൣെ߱ଶ ሾ‫ܯ‬ሿ ൅ ݆߱ሺሾ‫ܥ‬ሿ ൅ ȳሾ‫ܩ‬ሿሻ ൅ ሾ‫ܭ‬ሿ൧

ିଵ

(4)

The contribution of the cylindrical hydrodynamic bearings is inserted in the finite element model of the rotor by means of linear equivalent coefficients of damping and stiffness. The coefficients are obtained through the linearization, via first order Taylor expansion, of the nonlinear hydrodynamic supporting forces in the lubricant film, as shown in Eq. 5 ([Machado and Cavalca [33]). In the forces equation (Eq. 5), W is the bearing load, Kij are the stiffness coefficients, Bij are the damping coefficients and ȟ‫ݕ‬ത,ȟ‫ݖ‬ҧ, ȟ‫ݕ‬തሶ and ȟ‫ݖ‬ҧሶ are small perturbations in the journal bearing displacements and velocities around its equilibrium position. The coefficients are partial derivatives and the resulting differential expressions are approximated by finite differences using centered formulation as presented in Eq. 6; where ଓҧ, ଔҧand ଓҧሶ, ଔҧሶare, respectively, small perturbations in the journal bearing displacements and velocities around its equilibrium position at each rotating velocity. ‫ܨ‬௬ ൌ ‫ܭ‬௬௬ ȟ‫ݕ‬ത ൅ ‫ܭ‬௬௭ ȟ‫ݖ‬ҧ ൅ ‫ܤ‬௬௬ ȟ‫ݕ‬തሶ൅ ‫ܤ‬௬௭ ȟ‫ݖ‬ҧሶ ‫ܨ‬௭ ൌ ܹ ൅ ‫ܭ‬௭௬ ȟ‫ݕ‬ത ൅ ‫ܭ‬௭௭ ȟ‫ݖ‬ҧ ൅ ‫ܤ‬௭௬ ȟ‫ݕ‬തሶ൅ ‫ܤ‬௭௭ ȟ‫ݖ‬ҧሶ

(5)

6

‫ܭ‬௜௝ ൎ

ி೔శ౴ഢҧିி೔శ౴ഢҧ ଶ୼ఫҧ

and ‫ܤ‬௜௝ ൎ

ி೔శ౴ഢҧିி೔శ౴ഢҧ ଶ୼ఫҧሶ

(6)

The forces are obtained integrating the pressure field distribution inside the bearings, which is the result of the Reynolds equation solution. The worn bearing has a geometric discontinuity in the radial clearance and consequently, in the oil film. The solution of the Reynolds equation when the oil-film thickness is discontinuous was first presented by Arghir et al. [34] in one dimension (circumferential), and later expanded to two dimensions (radial and circumferential) by Machado and Cavalca [35]. The solution procedure takes into account a uniform mesh discretization of the fluid-film by finite volumes. The wear region of the bearing assumes abrasive wear and is based on the geometry proposed by Dufrane et al. [11], and experimentally validated by Hashimoto et al. [12]. This geometry considers that the wear has a uniform thickness in the axial direction and is symmetrically located at the bearing bottom. Even though wear due to startup and stops is generally located at the bottom of the bearing, wear due to mixed regime operation does not present this symmetry. In fact, the wear is developed through these two sources: start/stop cycles and mixed regime operation. Therefore, practical evidence has shown that the wear is commonly located in the fourth quadrant of the bearing. In this sense, the wear model proposed by Machado and Cavalca [33] included an angular displacement (γ) of the center of the wear, as shown in Fig. 1, and the authors experimentally validated the model comparing to literature. Thereby, the model from Machado and Cavalca [33] was used in this work.

Fig. 1. Worn bearing geometry. In the schematic of the worn bearing (Fig. 1) it can be seen that the wear introduces an additional oil layer with depth ߜ௛ ሺߠሻ in the region delimited by ߠ௦ and ߠ௙ , assumed symmetrical with relation to the point of highest depth (d0), resulting in a fluidfilm thickness ݄ሺߠሻ given by Eq. 7, where ݄଴ ሺߠሻ is the film thickness due to the shaft eccentricity.

7

݄ሺߠሻ ൌ ቊ

݄଴ ሺߠሻǡͲ ൑ ߠ ൑ ߠ௦ ǡ ߠ௙ ൑ ߠ ൑ ʹߨ ݄଴ ሺߠሻ ൅ ߜ௛ ሺߠሻǡߠ௦ ൏ ߠ ൏ ߠ௙ 

(7)

According to Fig. 1, the terms of Eq. 7 can be written in the inertial reference system (y,z) as presented in Eq. 8, where ey and ez are the components of the shaft eccentricity in y and z coordinates, respectively, d0 is the maximum wear depth and Cr is the bearing radial clearance. ݄଴ ሺߠሻ ൌ ‫ܥ‬௥ െ ݁௭ …‘•ሺߠ െ ߨȀʹሻ ൅ ݁௬ •‹ሺߠ െ ߨȀʹሻ ߜ௛ ሺߠሻ ൌ ݀଴ െ ‫ܥ‬௥ ሺͳ ൅ …‘•ሺߠ െ ߨȀʹሻሻ

(8)

The wear depth ߜ௛ ሺߠሻ is considered to be zero in the edges, ߠ ൌ ߠ௦ and ߠ ൌ ߠ௙ , which angular position is given by Eq. 9 ߠ௦ ൌ ߨȀʹ ൅ …‘• ିଵሺ݀଴ Ȁ‫ܥ‬௥ െ ͳሻ ൅ ߛ ߠ௙ ൌ ߨȀʹ െ …‘• ିଵ ሺ݀଴ Ȁ‫ܥ‬௥ െ ͳሻ ൅ ߛ

(9)

Finally, Eqs. 7-9 are inserted into the finite volumes model for the evaluation of pressure distribution, hydrodynamics forces and bearings linear dynamic coefficients. The wear parameters to be identified are the maximum wear depth (d0) and the wear angular position (γ), which allow the wear geometry to be obtained by means of Eqs. 8 and 9. 2.2. Directional response analysis The rotor orbit is obtained by combination of lateral displacements in y and z directions through time, thus, representing the precession movement or whirl (ψ). The precession movement can be forward, if it is in the same direction of the shaft spin, Fig. 2a, or backward, if it is in the opposite direction, Fig. 2b. The same principle can be extended for frequency domain analysis by converting the Frequency Response Function (FRF), in physical coordinates y and z, to the directional Frequency Response Function (dFRF), in directional coordinates f (forward) and b (backward). The transformation was presented by Lee [3] and previously applied by Cavalca and Okabe [7]. In Time Domain, a real signal can be represented by a Fourier series, as presented in Eq. 10 (Proakis and Manolakis [36]), and re-written in terms of complex exponentials as in Eq. 11. ݂ሺ‫ݐ‬ሻ ൌ σஶ ௡ୀଵሾܽ௡ ܿ‫ݏ݋‬ሺ߱௡ ‫ݐ‬ሻ ൅ ܾ௡ ‫݊݅ݏ‬ሺ߱௡ ‫ݐ‬ሻሿ ௔೙ ି௝௕೙

݂ሺ‫ݐ‬ሻ ൌ σஶ ௡ୀଵ ቂቀ

ଶ

௔೙ ା௝௕೙

ቁ ݁ ௝ఠ೙௧ ൅ ቀ

ଶ

ቁ ݁ ି௝ఠ೙௧ ቃ

(10) (11)

8

Fig. 2. Precession motion: (a) forward; (b) backward. Therefore, the steady state response of a rotational system submitted to a harmonic excitation is described by its Fourier series expansion as shown in Eq. 12. ଵ ௝ఠ೙ ௧ ‫ݕ‬ሺ‫ݐ‬ሻ ൌ σஶ ൅ തതത ܻ௡ ݁ ି௝ఠ೙௧ ൧ ௡ୀଵ ଶ ൣܻ௡ ݁

ଵ ௝ఠ೙ ௧ ‫ݖ‬ሺ‫ݐ‬ሻ ൌ σஶ ൅ ܼҧ௡ ݁ ି௝ఠ೙௧ ൧ ௡ୀଵ ଶ ൣܼ௡ ݁

where

௬

௬

ܻ௡ ൌ ܽ௡ െ ݆ܾ௡ ܼ௡ ൌ ܽ௡௭ െ ݆ܾ௡௭

(12)

The rotor response in the plane yz can be re-written using the complex notation ‫ݎ‬ሺ‫ݐ‬ሻ ൌ ‫ݕ‬ሺ‫ݐ‬ሻ ൅ ݆‫ݖ‬ሺ‫ݐ‬ሻ, resulting in Eq. 13. ଵ ௝ఠ೙ ௧ ൅ ሺܻത௡ ൅ ݆ܼ௡ҧ ሻ݁ ି௝ఠ೙௧ ൧ ‫ݎ‬ሺ‫ݐ‬ሻ ൌ σஶ ௡ୀଵ ଶ ൣሺܻ௡ ൅ ݆ܼ௡ ሻ݁

(13)

Alternatively, each harmonic of the rotor whirl can be described in directional coordinates (Fig. 3), i.e. by the sum of two rotating vectors - forward (f) and backward (b) – as in Eq. 14. ௝ఠ೙ ௧ ‫ݎ‬ሺ‫ݐ‬ሻ ൌ σஶ ൅ ܾ௡ ݁ ି௝ఠ೙௧ ൧ ௡ୀଵ ൣ݂௡ ݁

(14)

Comparing Eqs. 13 and 14, it is possible to obtain a relation between physical coordinates and direction coordinates as presented in Eq. 15. ͳȀʹ ݂ ܻ ൜ ൠ ൌ ሾ‫ܣ‬ሿିଵ ቄ ቅ ; ሾ‫ܣ‬ሿିଵ ൌ ൤ ത ͳȀʹ ܼ ܾ

݆Ȁʹ ൨ െ݆Ȁʹ

(15)

9

Fig. 3. Precession motion represented in directional coordinates. By means of the transformation matrix (Eq. 15), the rotor displacement and the external force vectors in physical (ሼ‫ݍ‬ሽ and ሼ‫ܨ‬ሽ) and in directional (ሼ‫݌‬ሽ and ൛ܲ෠ൟ) coordinates are related, as in Eq. 16. The same relations allow the conversion of the FRF (Eq. 4) into the dFRF as presented in Eqs. 17 and 18. Finally, considering a single input and output, the dFRF matrix has four terms as shown in Eq. 19. ሼ‫ݍ‬ሽ ൌ ሾ‫ܣ‬ሿሼ‫݌‬ሽ and ሼ‫ܨ‬ሽ ൌ ሾ‫ܣ‬ሿ൛ܲ෠ൟ

ሾ‫ܣ‬ሿሼ‫݌‬ሽ ൌ ሾ‫ܪ‬ሿିଵ ሾ‫ܣ‬ሿ൛ܲ෠ൟ  ื  ሼ‫݌‬ሽ ൌ ሾ‫ܣ‬ሿିଵ ሾ‫ܪ‬ሿିଵ ሾ‫ܣ‬ሿ൛ܲ෠ൟ ሾ݀‫ܪ‬ሿିଵ ൌ ሾ‫ܣ‬ሿିଵ ሾ‫ܪ‬ሿିଵ ሾ‫ܣ‬ሿ ݂݂ ሾ݀‫ܪ‬ሿିଵ ൌ ቈ ܾത݂

݂ܾത ቉ ܾതܾത

(16) (17) (18) (19)

2.3. Procedure for wear parameters identification The wear parameters identification method consists of solving an optimization problem. The objective function to be minimized compares the directional Frequency Response Function (dFRF) experimentally obtained with the response of the rotating system model. The objective function used, Fmin, is the classical least mean square function, as given by Eq. 20. In the objective function, ‫ܨ‬௘௫௣೔ and ‫ܨ‬௠௢ௗ೔ are, respectively, the i-th complex components of the experimental response and model response (dependent of {Xp}). The number of points of the dFRFs (‫ܨ‬௘௫௣೔ and ‫ܨ‬௠௢ௗ೔ ) which were used in the objective function (Fmin) is represented by n. {Xp} is the vector containing the parameters to be minimized (maximum wear depth, d0, and wear angular position, γ), whose upper and lower limits are given, respectively, by {Ll} and {Lu}. One of the parameters analyzed in the results of this work is the influence of the number of points

10 (n) in the dFRF to be used in the objective function. It is important to highlight that, since the dFRF is a set of complex numbers, the result of the least mean square function will also be a complex number. Thus, the objective function is obtained by taken the modulus of the resultant complex number. By taking the modulus after the least mean square calculation, the dFRF phase is taken into account, and the objective function still gives a real number (‫ܨ‬௠௜௡ ‫ א‬Թ). ்

‫ܨ‬௠௜௡ ൌ σ௡௜ ቚ൫‫ܨ‬௘௫௣೔ െ ‫ܨ‬௠௢ௗ೔ ൯ ൫‫ܨ‬௘௫௣೔ െ ‫ܨ‬௠௢ௗ೔ ൯ቚ ሼ‫ܮ‬௟ ሽ ൑ ൛ܺ௣ ൟ ൑ ሼ‫ܮ‬௨ ሽ

(20)

However, as stated by Arruda [37], the convergence of searching algorithms tends to be difficult when working with objective functions presented in the form of Eq. 20 due to the abrupt variations of magnitude, especially around resonance regions. The author found out that the use of a logarithmic scale improves the convergence of the searching algorithms. Therefore, the objective function proposed by Arruda [37] was used as presented in Eq. 21, resulting in the restricted optimization problem stated in Eq. 22. ி೘೚೏೔

ଶ

௟௢௚

‫ܨ‬௠௜௡ ൌ σ௡௜ൣŽ‘‰൫ห‫ܨ‬௠௢ௗ೔ ห൯ െ Ž‘‰൫ห‫ܨ‬௘௫௣೔ ห൯൧ ൌ σ௡௜ ൤Ž‘‰ ൬ฬ ி ௟௢௚

ி೘೚೏೔

݉݅݊݅݉݅‫ܨ݁ݖ‬௠௜௡ ൌ σ௡௜ ൤Ž‘‰ ൬ฬ ி

೐ೣ೛೔

ଶ

ฬ൰൨ ǡ

‫݋ݐ݀݁ݐ݆ܾܿ݁ݑݏ‬ሼ‫ܮ‬௟ ሽ ൑ ൛ܺ௣ ൟ ൑ ሼ‫ܮ‬௨ ሽ

೐ೣ೛೔

ฬ൰൨

ଶ

(21)

(22)

In order to improve the finding of global minimum of the obtained optimization problem, a MATLAB routine was implemented using the MultiStart method with the Interior Point Algorithm, both from the Optimization Toolbox. The Multi-Start method solves the optimization problem using the Interior Point Algorithm for a given mesh of starting points and selects the one with the lower objective function. More details regarding the Interior Point Algorithm are found in the works of Byrd et al. [38 and 39] and Waltz et al. [40]. At every step of the optimization process, the algorithm calculates the journal bearing coefficients and the dFRF response to evaluate the objective function. The influence of the mesh of starting points was also investigated in the results presented in this work. For example, if the mesh is taken as [10,40] μm for wear depth and [-15,15]° for wear angular position, its grid is composed by the points: {10 μm, 15°}, {10 μm, 15°}, {40 μm, -15°}, {40 μm, 15°}; as presented in Fig. 4. The same notation was used in section 3.2.

11

Fig. 4. Representation of the mesh of starting points given by [10,40] μm for wear depth and [-15,15]° for wear angular position. The procedure for the wear parameters identification is summarized in the following six steps: 1 – Initially, the test rig is assembled with the healthy bearing; 2 – Next, the FRF in physical coordinates is obtained using the excitation applied by a magnetic actuator and the displacements measured in the bearing; 3 – The finite element model of the rotor is adjusted in order to better represent the real system. This adjustment involves the number of shaft and disc elements and adjustment of structural damping coefficients (α and β in Eq. 2); 4 – After the model calibration tuning, the worn bearings are assembled in the test rig, whose wear parameters will be identified; 5 – The FRF of the new set-up is measured and converted into directional coordinates (dFRF); 6 – Finally, the identification algorithm is applied using the adjusted model and the experimental dFRF. 2.4. Experimental test rig The test rig used to validate the wear parameters identification method is shown in Fig. 5. The SAE 1020 steel shaft, of 600 mm length and 12 mm diameter, is supported by two cylindrical hydrodynamic bearings placed at 400 mm from each other, being the first bearing at 120 mm from the flexible coupling with the motor. The bearings have a 90 μm radial clearance, 18 mm length and 31 mm of inner diameter. The lubricant used is the oil AWS 32 from Castrol (ISO VG 32) at 28 °C. The journal of the magnetic actuator, with 40 mm diameter and 80 mm length, is placed at 200 mm from the first bearing. The actuator constantly applies a down force of 80 N in order to make the bearings operate in a higher eccentricity condition, at high rotation speed (Machado and Cavalca [33]). The control of the actuator force is based on magnetic flux control using

12 hall-sensors placed at the coils poles; more details regarding the magnetic actuator are to be found in the work of Mendes et al. [41].

Fig. 5. Test rig at the Laboratory of Rotating Machinery (University of Campinas). Fig. 6 shows the test rig instrumentation and control system. The WEG motor is driven by a frequency inverter WEG CFW-08, which is controlled via software using a serial communication module WEG XC8. The displacements of the shaft inside the bearings are measured by four Bently Nevada sensors. The displacement signals pass through signal conditioners, a DC filter to remove the offset and a low-pass analogic filter before being acquired by the National Instruments acquisition board (NI PCIe6363); the acquisitions were performed at a sample rate of 1024 samples/s. The magnetic actuator force is controlled via the same board, where the force reference signal goes to the control module, which sends a signal proportional to the required magnetic field to the PWM power amplifiers (4-Q-DC-Servoamplifier ADS50/5 from MAXON). The power amplifiers provide the necessary current to the magnetic actuator coils, and the magnetic fields are fed back via Hall sensors.

13

Fig. 6. Scheme of the acquisition and data processing. 3. Results and Discussion Although the identification procedure has six steps, the results are divided here into two parts: sub-section 3.1 contemplates the first three steps of the identification, related to the rotor model tuning, and sub-section 3.2 addresses the results regarding the wear parameters identification (last three steps). 3.1. Model Adjustment In order to allow the correct identification of the bearing wear parameters, the numerical model needs to be as close as possible to the real system. Therefore, the first three steps of the identification procedure, presented in section 2.3, concern the rotor model calibration. Initially, the test rig was assembled with the healthy bearing and operated at 50 Hz for the measurement of the four terms of the FRF. Fig. 7a represents the direct term of the FRF related to the y direction response of the rotor, which represents how an excitation in the y direction affects the y component of the rotor response. Fig. 7d represents the other direct term of the FRF and its interpretation is analogous to Fig. 7a, but by means of the z component. Fig. 7b contain the cross term of the FRF which represents the effect of an excitation in the z direction in the y component of the response. Fig. 7c is analogous to Fig. 7b, but representing the response of the z component to a y direction excitation. The experimental FRF is represented by the black lines in Fig. 7. It can be seen the peaks of the first bending modes of the rotor around 82 Hz, a narrow peak at 50 Hz (the rotational speed) due to unbalance and the oil-whirl peak at 25 Hz (half the

14 rotational speed). The peak at 50 Hz did not appear in the simulated results (gray curves in Fig. 7) once the simulated FRF was obtained from the solution of Eq. 4 (harmonic response of the system). A time simulation was not needed and, consequently, there was no unbalance force applied. The effect of the rotational speed in the model FRF is accounted for in the gyroscopic matrix and in the calculus of the bearings equivalent coefficients of stiffness and damping. Since the rotor is almost symmetric, the gyroscopic effect is not representative, and the peak at 82 Hz, correspondent to the first bending mode of the shaft, is equally prominent in the entire operational range of the rotor (up to 100) Hz. As the shaft is more flexible than the bearings and the load only affects the bearings equivalent coefficients, the bending mode frequency is not affected by load variations. However, the damping of the mode (bending mode peak shape) is affected by the load variation, once the journal bearings are responsible for most of the damping of the system. The physical reasoning for the damping variation is presented next with the oil-whirl peak description. The oil-whirl peak, near 25 Hz, corresponds to the sub-synchronous phenomenon, occurring at near half the rotational speed. Thus, its frequency will change as the rotational speed changes; and the peak also gets narrower as the rotational speed increases because its damping decreases. Also, the lower the load, the more prominent the oil-whirl becomes. Both phenomena occurs because when the rotor operates in a low eccentricity condition, either due to the self-centering effect (at high rotational speed) or a low load condition, the thicker fluid-film is less stiff and causes a disequilibrium of the stabilizing forces inside the bearing. A complete discussion regarding the stability of lightly loaded and heavily loaded rotors is to be found in the work of Mendes and Cavalca [44]. Consequently, depending on the operation condition of the rotor (load and rotational speed) during the wear parameters identification, the points used in the objective function should be changed. In section 3.2, a complete discussion is presented regarding the choice of the dFRF points to be used in the objective function. The same behavior is observed in the presence of different levels of bearing wear, once the wear is related to a change in the fluid film thickness. This dependence of the bending mode damping and the oil-whirl damping with the operating condition of the bearing is the reason why both peaks are suitable for the bearing wear parameters identification, as it will be shown in section 3.2. In all the FRFs measured in this work, the magnetic actuator was used as an external excitation source and the displacements were measured at the second bearing. White noise excitations of 30 N amplitude and 5 s duration were applied. The signals were measured at a sample rate of 1024 samples/s, thus resulting in a 0.2 Hz resolution in the frequency response. The total duration of each FRF matrix measurement was about 14 minutes, considering 100 averages and Hanning window.

15

(a)

(b)

(c) (d) Fig. 7. Comparison between the Frequency Response Functions (FRF) obtained experimentally and from the adjusted model for the healthy bearing: (a) yy; (b) yz; (c) zy; (d) zz. Next, the finite element model was compared and tuned with the measured FRF. The model discretization is presented in Fig. 8. It consists of 15 shaft elements (16 nodes) with dimensions given in Table 1, and 3 discs elements, as shown in Table 2. A hypothetic disc located in node 1 represents the inertia added to the shaft due to the coupling connecting the shaft and the motor. The linear dynamic coefficients of damping and stiffness of the first bearing (black triangles) are added in the

16 correspondent global matrixes in the position of node 4, and the coefficients of the second bearing (gray triangles) are added in the position of node 14. As previously mentioned, the coefficients were evaluated according to Machado and Cavalca [35]. The fluid viscosity at 28 °C is 0.1044 Pa.s and the bearing loads are 35.38 N and 58.82 N, respectively, for the first and second bearings. It is important to recall that, besides the rotor weight, the magnetic actuator applies a constant down force of 80 N in the shaft.

Fig. 8. Finite element model of the shaft. The steel shaft is considered to have an elastic modulus of 2.1x1011 N/m2, 7800 kg/m density, poisson coefficient 0.3 and a 0.796x1011 N/m2 shear modulus. The aluminum of the coupling disc has 0.69x1011 N/m2 elastic modulus, density of 2700 kg/m3, poisson coefficient 0.33 and a 0.26x1011 N/m2 shear modulus. The coefficients which relate the proportional damping matrix of the shaft [Cs] with the shaft mass [Ms] and stiffness [Ks] matrixes are, respectively, ߙ ൌ Ͳ and ߚ ൌ ͳǤ͸‫ିͲͳݔ‬ସ, and were obtained in previous tests and simulations (Santana et al. [42]). The coefficient ߙ is usually considered nonzero only for composite materials, once the masses of the fibers and the matrix have a significant influence on the behavior of such materials. Therefore the damping matrix is only a weighting of the stiffness matrix. The coefficient ߚ was then adjusted by comparing the damping (peak shape) of the first forward and backward bending modes through the model response and the experimental dFRF (Fig. 7). 3

Table 1. Dimensions of the shaft elements from the finite element model of the test rig. Elements Diameter (mm) Length (mm) 1,2 12 60 3,4,13,14 31 10 5,6,7,8 12 55 9,10 30 40 11,12 12 50 15 12 30 Table 2. Dimensions of the disc elements from the finite element model of the test rig. Inner Outer Nodes Length (mm) Material Diameter (mm) Diameter (mm) 1 12 25 20 Aluminum 9,10 30 40 40 Steel 1020 10,11 30 40 40 Steel 1020

17 The journal of the magnetic actuator is represented between nodes 9 and 11, and the magnetic force is considered to be applied in the central node of the journal (node 10). Due to the actuator journal length, the stiffness properties of the shaft are changed in the region between nodes 9 and 11. In order to compensate this effect in the FEM model, Lalanne and Ferraris [43] suggest the journal (which is a long disc) to be divided in more than one element. Each element of the journal is modeled part as a shaft element (in white), and part as disc element (in gray). Through such division, the local stiffness of the shaft is increased due to the use of a shaft element, and the remaining inertia of the journal is compensated using the disc element. The adjustment performed in the model was made by comparing the frequency of its first bending mode (around 82 Hz in Fig. 7) with the experimental result (FRF). As previously mentioned, the final dimensions for the shaft and disc elements (after the adjustment) are found in Tables 1 and 2. The comparison between the experimental FRF (black lines) and the FRF of the adjusted model (gray lines) is presented in Fig. 7 – recalling that the model FRF is calculated between nodes 10 (external excitation from the magnetic actuator) and 14 (displacements of the second hydrodynamic bearing). The results indicated that the model is suitable for wear parameters identification, as seen from the good agreement with the experimental FRFs. The only exception is an unknown (thus, not modeled) factor in yz term (Fig. 7c) of the experimental curve in the range of 85 to 90 Hz. 3.2. Experimental wear parameters identification Once the model is tuned, four bearings with different levels of wear were assembled in the test rig in the position of the second bearing (see Fig. 5 and Fig. 8), and the FRFs were experimentally obtained. Therefore, at the experimental tests, five sets of FRFs were obtained, one being for the healthy bearing and four for the bearings with different wear parameters, as described in Table 3. In order to verify the algorithm, the healthy bearing was also evaluated in all the sensitivity analyses performed. It should be remarked that the wear in the bearings was artificially made.

Worn Bearing Wear 0 Wear 1 Wear 2 Wear 3 Wear 4

Table 3. Parameters of the worn bearings tested. Maximum Wear depth (μm) Wear Angular Position (°) 0 0 22.27 11 22.27 -11 40.19 8 48.26 17

The evaluation of the wear parameter identification method was based on the sensitivity analysis of the method to the following parameters: - The use of the FRF versus dFRF, using the four terms of the matrixes or only the two crossed terms in the objective function; - The use of dFRF around the natural frequency peak versus the oil-whirl peak in the objective function; - The mesh discretization of starting points used by the MultiStart algorithm;

18

- The number of points of the dFRF used in the objective function. As previously mentioned, according to practice, the wear depth (d0) identification can be useful up to 50% of the bearing radial clearance; after which the bearing is not useful anymore. The wear angular position (γ) ranges between ±20°, once the main cause of bearing wear are the start/stop cycles, when the fluid-film is not completely formed and the rotor is located around the bottom of the bearing. In this sense, the search for the variables to be optimized (wear parameters) in the identification method were limited inside the intervals: Ͳߤ݉ ൑ ݀଴ ൑ ͷͲߤ݉ and െʹͲι ൑ ߛ ൑ ʹͲι.

The first parameter to be analyzed in the identification process was the use of the terms of the frequency response in physical (FRF) and directional (dFRF) coordinates – considering only the cross terms or the full matrix (cross and direct terms). In these cases, 7 points in the proximities of the first natural frequencies were used in the objective function (78.1 to 87.7 Hz, at every 1.6 Hz). The mesh of starting points used in the MultiStart algorithm was [0, 15, 30, 45] μm and [-15, -7.5, 0, 7.5, 15]°. The results are summarized in Tables 4 (wear depth) and 5 (wear angular position), where the identified parameters and relative error (in parenthesis) can be found. The error is not presented in the case of bearing without wear, once the error relative to a null value is not defined; thus, only in these cases, absolute values should be compared. Table 4. Identification results of wear depth and relative errors (in parenthesis) regarding the real values – FRF versus dFRF. Wear 0 Wear 1 Wear 2 Wear 3 Wear 4 0.22 21.05 21.45 42.21 45.78 FRF (-%) (5.48%) (3.68%) (-5.03%) (5.14%) (full) 0.12 0.06 0.01 28.19 28.47 FRF (-%) (99.73%) (99.96%) (29.86%) (41.01%) (cross) 0.03 21.92 22.33 41.98 46.58 dFRF (-%) (1.57%) (-0.27%) (-4.45%) (3.48%) (full) 0.05 21.62 20.76 42.15 46.79 dFRF (-%) (2.92%) (6.78%) (-4.88%) (3.05%) (cross) Table 5. Identification results of wear angular position and relative errors (in parenthesis) regarding the real values – FRF versus dFRF. Wear 0 Wear 1 Wear 2 Wear 3 Wear 4 -6.51 7.89 -7.88 6.12 19.85 FRF (-%) (28.27%) (-28.36%) (23.5%) (16.76%) (full) -6.95 -1.16 -2.28 19.74 19.31 FRF (-%) (110.55%) (-79.27%) (-146.75%) (-13.59%) (cross) 0.07 13.99 -7.85 5.45 15.85 dFRF (-%) (-27.18%) (-28.64%) (31.88%) (6.76%) (full) 1.23 8.11 -11.32 6.12 15.41 dFRF (-%) (26.27%) (2.91%) (23.5%) (9.35%) (cross) The results presented in Tables 4 and 5 show a similar behavior than in the simulations previously performed by Mendes et al. [45], the full dFRF presented smaller errors than the full FRF. When using only cross terms, the dFRF is much superior to the FRF, which have failed in the identification of wear in all cases. By

19 comparing the cases where the dFRF was used, it can be seen that the use of the crossed terms presented levels of error close to the use of all four terms. This observation experimentally validates the simulated results presented by Machado et al. [20], which have shown that the cross terms of the dFRF are more sensitive to bearing wear than the direct terms. This phenomenon is intrinsically related to the anisotropy level inserted in the system by the journal bearings. The high sensitivity of the dFRF cross terms is due to the variation of the anisotropy level which occur when the wear parameters changes. Such anisotropy effects are a result of the fluid-film thickness variation, as discussed at the beginning of section 3.1. Thereby, all the following identifications were performed using the full dFRF matrix, once the direct terms are less sensitive then the cross terms but still contribute to a more accurate identification. The next set of tests analyzed the influence of different number of points of the dFRF used in the objective function (n – see section 2.3, Eq. 22, for a complete explanation) for two cases: points in the surrounding of the natural frequency peak and in the surroundings of the oil-whirl peak. The range of dFRF points used in each test is presented in Table 6, and the identification results are presented in Tables 7 (depth) and 8 (angular position). In order to determine the number of points necessary in the objective function, the half-power points can be used. Thus, at least three points are necessary: the central frequency (the peak value) and the two half-power points. It is recalled that the halfpower points are those for which the energy dissipated per cycle of vibration is half the energy from the central frequency (peak value). Note that, when using logarithmic amplitude scales, the half-power points amplitudes are exactly 3dB lower than the peak amplitude. From a modal analysis point of view, the mode frequency can be obtained from the central point, and the modal damping can be estimated from the half-power points. Therefore, from the information contained in such points, the optimization algorithm should be able to identify the wear parameters which better adjust the model dFRF with the measured dFRF. It is important to take into account the difference in frequency of the first forward and first backward bending modes; see the peaks around 80 Hz in the experimental curves (black) in Fig. 7a and Fig. 7d. In this sense, for the first bending frequencies, four points were used in the range 78.1 Hz to 87.7 Hz, in order to cover both peaks. Since the oil-whirl is only a forward phenomenon, represented by a single peak around 25 Hz, only three points were used in the range 15.9 Hz to 34.3 Hz. These results were than compared when adding extra points between the points used in the previous cases; thus resulting in 7 points for the natural frequency peak and 5 points for the oil-whirl peak.

20 Table 6. Frequency ranges used for the results presented on Tables 7 and 8. Points of the dFRF used 78.1 ~ 87.7 Hz. at every 3.2 Hz Natural Freq. – 4 points 78.1 ~87.7 Hz. at every 1.6 Hz Natural Freq. – 7 points 15.9 ~ 34.3 Hz. at every 9.2 Hz Oil-Whirl – 3 points 15.9 ~ 34.3 Hz. at every 4.6 Hz Oil-Whirl – 5 points Table 7. Identification results of wear depth and relative errors (in parenthesis) regarding the real values – different number of points of the dFRF in the surroundings of the natural frequency and oil-whirl peaks. Wear 0 Wear 1 Wear 2 Wear 3 Wear 4 0.02 21.67 22.50 37.82 46.58 Natural Freq. (-%) (2.69%) (-1.03%) (13.4%) (3.48%) 4 points 0.03 21.92 22.33 41.98 46.58 Natural Freq. (-%) (1.57%) (-0.27%) (-4.45%) (3.48%) 7 points 3.01 25.47 23.24 41.03 39.89 Oil-Whirl (-%) (-14.37%) (-4.36%) (-2.09%) (17.34%) 3 points 4.35 24.58 23.41 37.65 37.81 Oil-Whirl (-%) (-10.37%) (-5.12%) (6.32%) (21.65%) 5 points Table 8. Identification results of wear angular position and relative errors (in parenthesis) regarding the real values – different number of points of the dFRF in the surroundings of the natural frequency and oil-whirl peaks. Wear 0 Wear 1 Wear 2 Wear 3 Wear 4 0.01 8.52 -7.50 13.40 15.01 Natural Freq. (-%) (22.55%) (-31.82%) (-67.5%) (11.71%) 4 points 0.07 13.99 -7.85 5.45 15.85 Natural Freq. (-%) (-27.18%) (-28.64%) (31.88%) (6.76%) 7 points -8.14 9.86 -11.14 4.98 15.89 Oil-Whirl (-%) (10.36%) (1.27%) (37.75%) (6.53%) 3 points -12.91 10.62 -12.59 4.65 12.93 Oil-Whirl (-%) (3.45%) (14.45%) (41.88%) (23.94%) 5 points Tables 7 and 8 indicate that the increase in the number of points (from 4 to 7 points) reflected in a small reduction in the identification errors, when using the region of the first natural frequencies of the shaft. However, when points around the oil-whirl frequency were used, the increment in the number of points, from 3 to 5 points, presented a small increment in the identification errors. Thus, the obtained results show that the algorithm have different sensitivities to the number of points depending on the region of the dFRF used in the objective function (the use of 7 points presented the best results for the region of the first natural frequencies, and the use 3 points presented slightly better results when using points around the oil-whirl peak). This behavior can be explained due to different shapes (smoother or sharper - Fig. 9 and Fig. 10) in different peaks (in this case, the oil-whirl and the natural frequency). Therefore, the number of points and the spacing required should be adjusted depending on the peak used in the identification process. This choice must allow the identification algorithm to correctly interpret and be sensitive to the peaks changes (system modal parameters changes) regarding the wear parameters variation. As previously mentioned, the central frequency and the half-power points can be used. In the cases studied in Tables 7 and 8, the increase in the number of points only increased the computational time, but without

21 a significant change in the results, once that the half-power points presented accurate results. Finally, different meshes of starting points are evaluated in the MultiStart algorithm. Two meshes were compared: the first is composed of 4 points for depth ([0, 15, 30, 45] μm) and 5 points for angular position ([-15, -7.5, 0, 7.5, 15]°), called 4/5; and the second, by 7 points for each parameter ([0, 7.5, 15, 22.5, 30, 37.5, 45] μm and [15, -10, -5, 0, 5, 10, 15]°); called 7/7. The identification results for wear depth are presented in Table 9, and the identification results for wear angular position in Table 10. Table 9. Identification results of wear depth and relative errors (in parenthesis) regarding the real values – different meshes of starting points. Wear 0 Wear 1 Wear 2 Wear 3 Wear 4 0.02 21.67 22.50 37.82 46.58 Mesh with (-%) (2.69%) (-1.03%) (5.9%) (3.48%) 4/5 points 0.02 21.67 22.49 37.82 46.57 Mesh with (-%) (2.69%) (-0.99%) (5.9%) (3.5%) 7/7 points Table 10. Identification results of wear angular position and relative errors (in parenthesis) regarding the real values – different meshes of starting points. Wear 0 Wear 1 Wear 2 Wear 3 Wear 4 0.01 8.52 -7.50 13.40 15.01 Mesh with (-%) (22.55%) (-31.82%) (-67.5%) (11.71%) 4/5 points -4.45 8.51 -7.51 13.45 14.99 Mesh with (-%) (22.64%) (-31.73%) (-68.13%) (11.82%) 7/7 points In Mendes et al. [45], similar simulated tests have shown that more points in the mesh of starting points improved the identification results; however, it is important to remark that the authors used the same simulated model to generate what would be an experimental dFRF and to perform the wear parameters identification. In the case of the experimental results of this work, presented in Tables 9 and 10, the similarity between the results for different meshes suggests that the quality of the parameters identification asymptotically tends to a limit. Such limit will always exist even though the mesh size overcome the problem of local minima of the objective function, given that the model will never perfectly represent the real system. Regarding all the identification results presented in this section, it is noticeable that the identification of wear angular position presented higher difficulty than the identification of wear depth. This behavior occurred because, as shown in simulations in the work of Machado et al. [20], the dFRF is much more sensitive to changes in wear depth than in wear angular position. Due to the manufacturing process (the wear in the bearings was artificially made), the accuracy obtained in the wear angular position was much lower than the accuracy in the wear depth. Therefore, higher uncertainties and, consequently, higher errors in the identification of the angular position were obtained. In order to illustrate the quality in which the model can represent the system response, Fig. 9 and Fig. 10 present the reconstructions of the directional frequency response (dFRF) for the case of Wear 1 (full and cross) showed in Tables 4 and 5. Both experimental and simulated dFRFs, obtained with the identified wear parameters are

22 shown and compared. This case was chosen to show that even with a low wear depth the model can satisfactorily represent the real system

(a)

(b)

(c) (d) Fig. 9. Comparison between the directional Frequency Response Functions (dFRF) obtained experimentally and from the model with the adjusted parameters using the full matrix for the bearing with Wear 1: (a) ff; (b) fb; (c) bf; (d) bb.

23

(a)

(b)

(c) (d) Fig. 10. Comparison between the directional Frequency Response Functions (dFRF) obtained experimentally and from the model with the adjusted parameters using the cross matrix for the bearing with Wear 1: (a) ff; (b) fb; (c) bf; (d) bb. Comparing Fig. 9 and Fig. 10, the adjustment using the full and cross dFRFs matrices showed similar results, as expected, in agreement with the levels of error obtained. It is remarkable that the identification method is capable of identifying the wear parameters despite the differences between the rebuild response using the model and the experimental curve. These results reinforce the work of Machado et al. [20], which have shown that the dFRF (manly the cross terms) are more sensitive to the wear

24 parameters variation than the FRF, and consequently, this sensitivity is the key feature for the identification. The graphics of Fig. 9 and Fig. 10 are also in accordance with the results in Tables 7 and 8, showing that the natural frequency peak is better adjusted than the oil-whirl peak; and consequently, the natural frequency peak is more suitable for the identification. 3.3. Method Limitations and Future Improvements One limitation of the presented method comes from the approximation of the wear geometry as a function of the maximum depth d0 (see Eqs. 8b and 9). Even the artificially made wear in the bearings used in this work did not perfectly match such geometry. The combination of the geometry approximation with the low sensitivity of the dFRF to the wear angular position resulted in the high errors in this parameter identification. Another limitation, related to the wear geometry used in this work, is that the wear is considered constant in the axial direction; and so are the artificially made wear in the bearings. However, such assumption is not true for real worn bearings and, therefore, the robustness of this method must be tested in worn bearings whose wear is not constant in the axial direction. In this sense, the proposed method could be changed in two points: 1. The wear model could be extended to a three dimensional model (by including the axial direction), thus the wear would not be constant in the axial direction; 2. The number of parameters to be identified could be increased, so that the geometry would not be only a function of the maximum depth, but would have more parameters regarding its shape (e.g. a half ellipse with major semi-axis to be identified – the minor semi-axis would be the maximum depth). The inclusion of more wear parameters in the identification is expected to lead to a more complex optimization problem and, consequently, to more computational time, without a practical improvement in the result. This statement is based on this work results observation, where the parameter that most affects the bearing behavior is wear depth; hence, this is the parameter with higher potential to compromise the bearing performance. As previously mentioned, the dFRF is more sensitive to the wear depth than to the wear angular position, and the depth identification presented high precision (below 5% for the best case – Full dFRF) besides of the limitations presented. The first point could be feasible if the geometry in the axial direction was described by another equation depending on the maximum depth (similar to the equation for the circumferential direction). Therefore, without increasing the number of parameters, the identification of real worn bearings could be improved. Despite the method limitations, the results obtained have shown the proposed method robustness in the estimation of the wear parameters. Therefore, the method is a promising tool for the use in real rotors, which present a higher challenge in the wear identification due to a series of complications not present in the laboratory test rig; such as many bearings, high level of noise and influence of foundation and fluid dynamics.

25

4. Conclusions This work presented an experimental validation of a model-based method for identification of the parameters of cylindrical hydrodynamic bearing wear. The tests performed were also intended for sensitive analyses regarding several parameters of the method. The first identification results have shown that the dFRF is more suitable for the wear parameters identification than the FRF; and the use of the full dFRF matrix presented similar performance to the use of only the cross terms of the dFRF matrix. These results are in accordance with the findings of Machado et al. [20] that the cross terms of the dFRF are more sensitive to changes in the bearing wear parameters than the direct terms of the dFRF. The next set of tests have shown that both the oil-whirl and the natural frequency peaks are sensitive to changes in the bearing wear parameters and can be used in the identification method proposed. From the obtained results, it was concluded that, depending on the peak chosen and, obviously, on the characteristics of the system under analysis, different number of points with different frequency spacing are required in order to allow the method to correctly identify the wear parameters. The half-power points and the distribution of the peaks chosen for the identification (first forward and first backward bending modes peaks, or the oil-whirl peak) can give an important hint about how to tune these parameters. It is recalled that the points mentioned refer to the points of the dFRF used in the objective function of the optimization problem. The last analysis regards the size of the mesh of starting points. The MultiStart method was used along with the optimization method chosen (Interior Point Algorithm) to avoid local minima of the objective function. The experimental tests indicated that, since the model will never perfectly represent the system under test, the model uncertainties limit the increase in performance obtained from increase of mesh size of starting points. For real bearing users, the model adjustment can be performed during the rotor commissioning; thus, providing a representative model that can be used as a reference for the system with healthy bearings. Before the wear identification with the proposed algorithm, a previous identification step is necessary in order to identify that the changes in the rotor behavior are due to bearing wear. The next step of this research is to try to isolate bearing wear from other failures, such as misalignment. From the results, the presented method is a very promising tool for bearing wear parameters identification. As previously mentioned, an important advantage of the method is that it can obtain the necessary information for the wear parameters identification directly form usual vibration measurements at the bearings; in contrast with most of the wear identification methods based on non-trivial measurements of bearing characteristics, such as bearing locus, pressure distribution or Sommerfeld number.

26 5. Acknowledgment The authors would like to thank CAPES, CNPq and grant #2015/20363-6 from the São Paulo Research Foundation (FAPESP) for the financial support to this research. 6. Conflict of Interests The authors declare that there is no conflict of interests regarding the publication of this paper.

References [1] GASCH, R.; PFÜTZNER, H. Rotordynamik. Springer-Verlag, BerlinHeidelberg, 188p., 1975. [2] NORDMANN, R. Modal Parameter identification and sensitivity analisys in rotating machinery. University of Kaiserlautern, Departament of Mechanical Engineering, German Federal Republic, 1982. [3] LEE, C. W. A complex modal testing theory for rotating machinery. Mechanical Systems and Signal Processing, v. 5, p. 119-137, 1991. [4] KESSLER, C. L. Complex Modal Analysis of Rotating Machinery. Thesis, University of Cincinnati, Cincinnati, OH, United States, 1999. [5] MESQUITA, A. L.; DIAS, M. J.; MIRANDA, U. A. A Comparison between the Traditional Frequency Response Function (FRF) and the directional Frequency Response Function (dFRF) in Rotordynamic Analysis. Mecánica Computacional, v. 21, p. 2227-2246, MECOM – First South-American Congress on Computational Mechanics, Santa Fé - Paraná, Argentina, 2002. [6] CAVALCA, K. L.; CAVALCANTE, P. F.; OKABE, E. P. An investigation on the influence of the supporting structure on the dynamics of the rotor system. Mechanical Systems And Signal Processing, Cambridge, UK, v. 19, n. 1, p. 157-174, 2006. [7] CAVALCA, K. L.; OKABE, E. P. On the analysis of rotor-bearing-foundation systems. In: IUTAM Bookseries - Symposium on Emerging trends in rotor dynamics, 1ed. Dordrecht: Springer, v. 25, p. 89-101, 2010. [8] KATZENMEIER, G. The influence of materials and surface quality on wear behaviour and load capacity of journal bearings. In: Proceedings of the Institution of Mechanichal Engineers’ Tribology Convention, p. 80–4, 1972. [9] DUCKWORTH, W. E.; FORRESTER, P. B. Wear of lubricated journal bearings. Proceedings of the Institution of Mechanical Engineers Conference on Lubrication and Wear, London, p. 714-719, 1957.

27 [10] MOKHTAR, M.; HOWARTH, R.; DAVIES, P. Wear characteristics of plain hydrodynamic journal bearings during repeated starting and stopping. ASLE Transactions, v. 20, p. 191-194, 1977. [11] DUFRANE, K .F.; KANNEL, J. W.; MCCLOSKEY, T. H. Wear of steam turbine journal bearings at low operating speeds. Journal of Lubrication Technology, v. 105, p. 313–317, 1983. [12] HASHIMOTO, H.; WADA, S.; NOJIMA K. Performance characteristics of worn journal bearings in both laminar and turbulent regimes. Part I: Steady-state characteristics. ASLE Transactions, v. 29, p. 565-571, 1986. [13] SCHARRER, J. K.; HECHT, R. F.; HIBBS, R. I. The effects of wear on the rotordynamic coefficients of a hydrostatic journal bearing. ASME Journal of Tribology, p. 113:210-3, 1991. [14] SUZUKI, K.; TANAKA, M. Stability characteristics of worn journal bearing. In: Proceedings of the Asia–Pacific Vibration Conference, Kuala Lumpur, p. 296–301, 1995. [15] KUMAR, A.; MISHRA, S. Stability of a rigid rotor in turbulent hydrodynamic worn journal bearings. Wear, v. 193, p. 25-30, 1996. [16] KUMAR, A., MISHRA, S. Steady-state analysis of non-circular worn journal bearings in non-laminar lubrication regimes. Tribology International, p.29:493–8, 1996. [17] FILLON, M.; BOUYER, J. M. Thermohydrodynamic analysis of a worn plain journal bearing. Tribology International, p. 129-136, 2004. [18] BOUYER, J. M.; FILLON, M.; PIERRE-DANOS, I. Influence of wear on the behavior of a two-lobe hydrodynamic journal bearing subjected to numerous startups and stops. Journal of Tribology, p. 205-208, 2007. [19] MUZAKKIR, S. M.; HIRANI, H.; THAKRE, G. D. Experimental investigations on effectiveness of axial and circumferential grooves in minimizing wear of journal bearing operating in mixed lubrication regime. International Journal of Current Engineering and Technology, v. 5, p. 486-489, 2015. [20] MACHADO, T. H.; MENDES, R. U.; CAVALCA, K. L. Directional Frequency Response Applied to Wear Parameter Identification in Hydrodynamic Bearings. Mechanics Research Communications, v. 74, p. 60-71, 2016. [21] WU, T.; MAO, J.; DONG, G.; XU, H.; XIE, Y. Journal bearing wear monitoring via on-line visual ferrograph. Advanced Materials Research, v. 44, p. 189194, 2008. [22] GERTZOS, K. P.; NIKOLAKOPOULOS, P. G.; CHASALEVRIS, A. C.; PAPADOPOULOS, C. A. Wear identification in rotor-bearing systems by measurements of dynamic bearing characteristics. Computers and Structures, v. 89, p. 55-66, 2011.

28

[23] OKAC, H.; LOPARO, K. A.; DISCENZO, F. M. Online tracking of bearing wear using wavelet packet decomposition and probabilistic modeling: A method for bearing prognostics. Journal of Sound and Vibration, v. 302, p. 951-961, 2007. [24] CHASALEVRIS, A. C.; NIKOLAKOPOULOS, P. G.; PAPADOPOULOS, C. A. Dynamic effect of bearing wear on rotor-bearing system response. ASME Journal of Vibration and Acoustics, v. 135, Issue 1, 2013. [25] PAPADOPOULOS, C. A.; NIKOLAKOPOULOS, P. G.; GOUNARIS, G. D. Identification of clearances and stability analysis for a rotor-journal bearing system. Mechanism and Machine Theory, v. 43, p. 411-426, 2008. [26] MACHADO, T. H.; CAVALCA, K. L. Geometric discontinuities identification in hydrodynamic bearings. In: Proceedings of 9th IFToMM International Conference on Rotor Dynamics. Politecnico di Milano, v. 1, p. 1-10, Milan, Italy, 2014. [27] BAUMANN, U. Rotordynamic Stability Tests on High-Pressure Radial Compressors. Proceedings of the Twenty-Eighth Turbomachinery Symposium, Turbomachinery Laboratory, Texas A&M University, College Station, Texas, p. 115122, 1999. [28] MOORE, J. J.; WALKER, S. T.; KUZDZAL, M. J. Rotordynamic Stability Measurement During Full-Load Full-Pressure Testing of a 6000 PSI Reinjection Centrifugal Compressor. Proceedings of the Thirty-First Turbomachinery Symposium, Turbomachinery Laboratory, Texas A&M University, College Station, Texas, p. 29-38, 2002. [29] SOROKES, J. M.; SOULAS, T. A.; KOCH, J. M.; GILARRANZ, J. L. FullScale Aerodynamic and Rotordynamic Testing for Large Centrifugal Compressors. Proceedings of the Thirty-Eighth Turbomachinery Symposium, Turbomachinery Laboratory, Texas A&M University, College Station, Texas, p. 71-79, 2009. [30] BIDAUT, Y.; BAUMANN, U.; Al-HARTHY, S. M. H. Rotordynamic Stability of a 9500 PSI Reinjection Centrifugal Compressor Equipped with a Hole Pattern Seal— Measurement Versus Prediction Taking into Account the Operational Boundary Conditions. Proceedings of the Thirty-Eighth Turbomachinery Symposium, Turbomachinery Laboratory, Texas A&M University, College Station, Texas, p. 251259, 2009. [31] PETTINATO, B. C.; CLOUD, C. H.; CAMPOS, R. S. Shop acceptance testing of compressor rotordynamic stability and theoretical correlation. Proceedings of the thirty-ninth turbomachinery symposium, p. 31-42, 2010. [32] NELSON, H. D. A Finite Rotating Shaft Element Using Timoshenko Beam Theory. Journal of Mechanical Design, v.102, p.793-803, October 1980. [33] MACHADO, T. H.; CAVALCA, K. L. Modeling of Hydrodynamic Bearing Wear in Rotor-bearing Systems. Mechanics Research Communications, v. 69, p. 15-23, 2015.

29

[34] ARGHIR, M.; ALSAYED, A.; NICOLAS, D. The finite volume solution of the Reynolds equation of lubrication with film discontinuities. International Journal of Mechanical Sciences, v. 44, p. 2119-2132, 2002. [35] MACHADO, T. H.; CAVALCA, K. L. Dynamic Analysis of Cylindrical Hydrodynamic Bearings with Geometric Discontinuities. In: International Conference on Vibration Problems, ICoVP-2011, September 5-8, Prague, Czech Republic, 2011. [36] PROAKIS, J. G.; MANOLAKIS, D. G. Digital Signal Processing: Principles, Algorithms and Applications. Prentice Hall, 3rd edition, 1016p., 1995. [37] ARRUDA, J. R. F. Objective Functions for the Curve-Fit of Frequency Response Functions. AIAA Journal, v. 30, n. 3, p. 855-857, 1992. [38] BYRD, R. H.; HRIBAR, M. E.; NOCEDAL, J. An Interior Point Algorithm for Large-Scale Nonlinear Programming. SIAM Journal on Optimization, v. 9, n. 4, p. 877–900, 1999. [39] BYRD, R. H.; GILBERT, J. C.; NOCEDAL, J. A Trust Region Method Based on Interior Point Techniques for Nonlinear Programming. Mathematical Programming, v. 89, n. 1, p. 149-185, 2000. [40] WALTZ, R. A.; MORALES, J. L.; NOCEDAL, J.; ORBAN, D. An interior algorithm for nonlinear optimization that combines line search and trust region steps. Mathematical Programming, v. 107, n. 3, p. 391–408, 2006. [41] MENDES, R. U.; CAVALCA, K. L.; FERREIRA, L. O. S. Analysis of a complete model of rotating machinery excited by magnetic actuator system, Proc IMechE Part C: J Mechanical Engineering Science, pp. 1–17, 2012. [42] SANTANA, P. M. F.; CAVALCA, K. L.; OKABE, E. P.; MACHADO, T. H., Complex Response of a Rotor-Bearing-Foundation System. Proceedings of IFToMM – 8th International Conference on Rotor Dynamics, Seoul, Korea, pp. 231238, 2010. [43] LALANNE, M.; FERRARIS, G. Rotordynamics Prediction in Engineering. John Wiley & Sons Ltd., Chichester, England, 1990. [44] MENDES, R. U.; CAVALCA, K. L. On the Instability Threshold of Journal Bearing Supported Rotors. International Journal of Rotating Machinery, v. 2014, p. 117, 2014. [45] MENDES, R. U.; MACHADO, T. H.; CAVALCA, K. L. Evaluation of a Model Based Identification Method for Hydrodynamic Bearing Wear. Proceedings of the IMechE – Vibrations in Rotating Machinery, VIRM11, Manchester, England, 2016. Highlights

30 Experimental analysis of the wear in hydrodynamic bearings. Sensitivity of the directional frequency response function (dFRF) due to wear. Wear parameters identification in hydrodynamic bearings using the FRF and dFRF. Analysis of the wear identification method to the variation of several parameters. Advantage of the use of vibration measurements made directly on the rotor.

Test Rig

white noise applied by magnetic actuator

displacement measurements in the bearing

Graphical Abstract (for review)

f

b

FRF in Directional Coordinates (dFRF)

Rotor FEM model Worn Bearing FVM model

Optimization algorithm using system model and measured dFRF d0 Wear Parameters Identified: - Maximum Depth (d0) - Angular Position (γ)

γ