Paper received: 14.06.2007 Paper accepted: 07.07.2008 Numerical and Experimental Analysis of a Shaft Bow Influence on a Rotor to Stator Contact Dynamics Sanjin Braut - Roberto Zigulic1 - Mirko Butkovic2 1 University of Rijeka, Faculty of Engineering, Rijeka, Croatia 2 Polytechnic of Karlovac, Karlovac, Croatia The shaft bow problem presents a real situation especially in case of slender rotors. This paper investigates the shaft bow influence on the rotor-stator contact dynamics. For this purpose the rotor is described as a simple Jeffcott model and the stator as an elastically suspended rigid ring. To test the numerical model, except a usual run down analysis, an emergency shut down after the sudden rotor unbalance increase is also analyzed. Numerical integration is carried out by the fourth-order Runge-Kutta method. Two different normal force models for the rotor-stator interaction are analyzed. For analyzed parameters, both linear and nonlinear (Hunt and Crossley), normal force models gave similar rotor and stator responses. To confirm some of the results and to tune the numerical model, the experimental investigation on the test rig was conducted. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: rotor-stator contact, shaft bow, numerical simulations, rotor vibration, vibration measurements 0 INTRODUCTION Contact or rub between rotor and stator is one of the most intensive research subjects in rotor dynamics. Most papers which are considering the rotor-stator contact phenomena can be classified in rigid rotor disc-rigid stator contact [1] to [3], bladed disc-stator contacts [4] and rotor-stator contact in retainer bearings when the rotor is supported by active magnetic bearings [5] and [6]. Choy and Padovan [1] developed a general analytical rub model using the following assumptions: simple Jeffcott rotor model, linear stiffness and damping characteristics, rigid casing supported by springs acting in the radial direction, mass inertia of the casing small enough to be neglected, simple Coulomb friction and onset of rub caused by unbalance. Bartha [2] performed an extensive numerical and experimental research of the backward whirl of rotors considering the rigid and elastically suspended stator. Von Grol and Ewins compared measurements and simulations for a windmilling imbalance in aero-engines which is very similar to the classical rotor to stator contact. They have revealed the rotor response rich in different subharmonics. Influence of torsion on the rotor-stator contact was introduced by Edwards et al. [7]. Attention was paid to the effects of torsion on a steady state response of a system experiencing rotor to stator contact. The analyzed Jeffcott rotor model had only 3 degrees of freedom - d.o.f. because the stator was assumed to be rigid. Zigulic et al. [8] considered nonlinear dynamics of a rotor supported by two dry-friction bearings. The rotor was modeled with FEM-Finite element method and the obtained nonlinear system was integrated with Hilbert-Hughes-Taylor - HHT a method. Karpenko et al. [9] have shown how the preloading of the snubber ring could stabilize the dynamic response of the rotor. Theoretical predictions of the two-degrees-of-freedom Jeffcott rotor model with the preloaded snubber ring subjected to out-of-balance excitation was studied in [10], while an experimental verification with special attention to the analysis of the shaft rotational speed (excitation frequency) as well as different eccentricities and preloading was given in [11]. This paper considers the contact between the rotor, with a bowed shaft, and the stator at the position of rotor disc (seal), after the appearance of sudden unbalance due to blade loss. Linear and nonlinear models of the normal contact force are compared for different shaft bow - mass unbalance combinations. In this study, not only lateral but also torsional d.o.f. are taken into account allowing the additional extension of the *Corr. Author's Address: University of Rijeka, Faculty of Engineering, Vukovarska 58, HR-51000, Rijeka, Croatia, sanjin.braut@riteh.hr model with mechanical model of the induction motor. Results of numerical simulations and experiments are presented and compared using the rotor and stator responses in time domain, spectral maps and rotor orbits 1 MODELING OF THE ROTOR-STATOR SYSTEM A rather simple mechanical model for the rotor-stator interaction description is used. The rotor is described by Jeffcott model with a shaft bow while the stator is modeled as an elastically suspended rigid ring. Except traditional consideration of the rotor lateral d.o.f. (xr, yr), torsional d.o.f. (&r) is taken into account allowing an extension of the model with torsional d.o.f. of the induction motor (&m) as well as appropriate mechanical model. To keep consistency of the described consideration the stator had additional three d.o.f., two lateral (xs, ys) and one torsional (&s). Fig. 1 shows the considered rotor-stator system. Next sections briefly discuss the differential equations of motion, the contact model and mechanical model of the induction motor. Motor Fig.1. Extended Jeffcott rotor model 1.1. Differential equations of motion The equations of motion for the system described above can be derived using the Lagrange equations and have the following form: mrxr + crxr + kr xr = mreip2 cos&+mre&sin&+kr xr0 - FCx my + cryr + kryr = m^eqf sin&-mr eq>oosp+ky - Fey J&& + Crt (r ) + krt (r ) = Mfr1 -MC ms xs + Cs xs + ks xs = FCx ms ys + cy + ks ys = FCy JM + Cst& + kst& = MC Jm&m + Crt (m)+ krt (m) = Mf2 +Mm (1) The Equation (1) is taken from literature [12] and presents a further development of the models found in [2], [7] and [13], where mr, ms, c, cs, k, ks are mass parameters, lateral damping coefficients and stiffnesses of the rotor and stator respectively, crt, cst, krt and kst are torsional damping coefficients and stifnesses of the rotor shaft and stator, e is the rotor mass eccentricity, Jr, Js and Jm represent mass moments of inertia of the rotor, stator and induction motor respectively. Mfr1 and Mfr2 represent the torques of total losses (friction in bearings and losses in motor fan). These torques were identified from the measurements conducted on the test rig. To simplify the already complicated model, rotor and stator are assumed as isotropic regarding lateral stiffness and damping coefficient, what is still common practice in literature [14]. Vector rr0 from the Fig. 2 represents the static equilibrium position of the bowed shaft [13], where xr0 and yr0 are its X,Y components and are given by xr0 = r0x cos - r0y sin yr0 = r0x sin + r0x cos (2) Components r0x and r0y represent the static equilibrium position of the bowed shaft in the rotor fixed x,y,z coordinate system. yA Fig.2. Jeffcott rotor with a bowed shaft Contact forces FCx and FCy as well as the contact moment MC are defined, according to Fig. 3, by the following expression: FCx = FNcos7- FTsinY FCy = Fn sin Y + Ft cos y Mc = ft r (3) where the normal FN and the tangential FT contact forces have positive value when they are acting on the stator. Fig. 3 presents the geometry of the rotor - stator contact model with force definition. 1.2. Contact Models Traditionally, for the rotor-stator impact modeling, two different methods have been applied, namely the Newton's restitution coefficient model and the Contact forceindentation model. The first model is based on the simplifying assumption of perfectly rigid bodies. Actual physical objects are compliant and hence the impact duration is strictly greater than zero. This more realistic view of impact phenomena led many researchers to consider the continuous-dynamics models of collision where bodies deform during impact and the collision dynamics is treated as a continuous-time dynamic phenomenon. In its general form, the forceindentation relationship looks like [15]: Fn = Fk (¿) + F (¿, Cr ft = Mc fn Sgn(vsl) (16) where ¿uc is the coefficient of friction and vsl is the sliding velocity between rotor and stator surfaces in contact. The angle y from equation (3) represents the normal force angle or direction of the smallest gap when rotor and stator are not in contact i.e. tany = yr— ys (17) To determine S and vsl further vector calculation should be taken. Fig. 6 shows the velocity definition for the rotor and stator in contact. If the normal direction is known, n0 = M , r = ^ ry] = [ - — yr — ys] (18) as well as the velocities of rotor and stator centers Vr = [Vxr,Vyr] and Vs = [vXs,VyS], one can easily obtain the relative normal rotor-to-stator velocity, S = V — v„ (19) where vrn = vr ■ n0 and vsn = vs ■ n0 represent the rotor and stator velocities in normal direction. Fig.6. Velocity definition for rotor and stator in contact The sliding velocity Vsl can be further defined as, where Mmax and smn are the maximum motor torque or breakdown torque and slip for maximum motor torque, s is the actual slip of motors rotor with respect to the synchronous electromagnetic field and p = fsf /f is the ratio of actual frequency given from frequency inverter and nominal frequency given by electrical network (f =50 Hz). If we replace Mmax in Equation (21) with Mmaxb = -Mmax ( + smn )/(1-smn ) the KloSS expression for regenerative braking is obtained. The regenerative braking starts when actual speed of the motor's rotor is greater than the synchronous speed. The motor is then behaving as an electromagnetic brake. Fig. 7 shows the speed torque characteristics of an induction motor. Fig.7. Speed torque characteristics of an induction motor where vit = vrt 0 +pr R and vst = vst0 +ps (R + Cr) represent the rotor and stator velocities on its periphery in tangential direction. R is the radius of the rotor while pr, ps are angular velocities of the rotor and stator. 1.3. Mechanical Model of the Induction Motor The driving moment of the induction electric motor has been modeled according to the Kloss expression [18] and [19] Mm = Mm 2 + 2sm sp mn sp ■ 2sm (20) 2 NUMERICAL SIMULATIONS The numerical integration has been carried out by the fourth-order Runge-Kutta method. Time step of At = 2 10-5 s, has been applied in all simulations where the contact between rotor and stator was reasonably expected. Otherwise At = 1-10-4 s has been used. In some earlier simulations a greater time step was used but it turned out that it couldn't describe some specific rotor behaviors. Time step during one simulation has been fixed, so the subsequent frequency analysis could be done. The parameters used in all calculations are as follows: rotor mass mr = 4.25 kg, stator mass ms = 3.838 kg, mass moments of inertia Jr = 4.532-10-3 kgm2, Jm = 2.115-10-3 kgm2, Js = (21) 2.11-10-2 kgm2, rotor disc radius R = 0.06 m, stiffnesses kr = 131540 N/m, krt = 1080 Nm/rad, ks = 1.237-106 N/m, kst = 7917 Nm/rad, damping coefficients cr = 22.43 Ns/m, crt = 0.0374 Nms/rad, cs = 11.96 Ns/m, cst = 0.0258 Nms/rad, radial clearance Cr = 0.4 mm, motor breakdown torque Mmax = 10.25 Nm and motor speed of max. breakdown torque (at power network frequency f = 50 Hz) nm = 2400 rpm. The identified dependence between angular deceleration and angular velocity of the rotor, obtained from free run down tests (induction motor switched off), had the following form e = -4-10-6®3 +0.002® -0.5376®-10.513 rad/s2 (22) 2.1. Numerical Simulation of Motor Controlled Run-Down (RD) Analysis The numerical RD analysis has been performed to identify the basic dynamic characteristics of the rotor-stator system and to correlate them with the experimental results. Knowing, from the modal testing (rotor response on the impulse force excitation at standstill) that the first rotor natural bending frequency is equal to f = 28.0 Hz and according to the Jeffcott (or Laval) theory it is the only natural frequency, RD analysis has been focused on speed range n = 40 to 0 Hz. The second rotor natural bending frequency obtained by the modal testing has been f2 = 148.5 Hz hence it has been far away enough from the first natural frequency and its influence has been negligible. The basic goal of this part of numerical analysis was to discover the appropriate relation of the mass unbalance eccentricity vector e and the static equilibrium position of the bowed shaft rr0 according to the measured maximum lateral displacement amplitude of the rotor disc and the qualitative shape of the rotor lateral response in time domain, Fig. 8. Final parameters have been e = 0.046 mm, rr0 = 0.05 mm and phase lag of the e in relation to rr0 i.e. a0 = n rad (see Fig. 2). RD analysis has been controlled by a linear frequency ramp (speed law) of the synchronous electromagnetic field of induction motor, fsf = 40 to 0 Hz in the time period of 34 s, while its rotor behaved according to the presented mechanical model. At the beginning of the numerical simulation, 2 s have been taken additionally for disappearance of any transients caused by initial conditions. Fig.8. Simulated lateral rotor response in horizont direction for RD analysis and speed law Fig. 9 shows the spectral map of the lateral rotor response in the horizontal direction where as expected only the first order can be seen. On the same figure at the beginning of the simulation, the rotor critical speed is excited but only because of the numerical transient due to initial conditions. 200 150 I 100 Fig.9. Spectral map of the lateral rotor response in horizontal direction xr 2.2. Numerical Simulations of Sudden Rotor Unbalance Increase (SRUI) Analysis The basic assumption has been that the rotor is running at a nominal speed above the first rotor critical speed and that the rotor is well balanced. To achieve the initial condition which is similar to the steadily running speed and to avoid any transients, the sudden unbalance was introduced 1 second after the start of simulation. After sudden appearance of the additional unbalance madd it has been assumed that the permanent monitoring system of the turbo machine has recorded unallowable vibrations and as a precaution measure it starts an immediate shut down procedure. This was accomplished by switching off the induction motor. Then the whole rotor starts to decelerate due to given torques of losses, expressed by the identified rotor angular deceleration mentioned above. The rotor is then approaching to its critical speed and because of additional unbalance muadd it contacts the stator. When the additional unbalance is big enough, the rotor contacts the stator immediately after appearance of the additional unbalance i.e. before switching off the motor. The time period for response of the shut down procedure after sudden unbalance appearance has been set to 1 s. Two different contact normal force models have been compared i.e. linear and nonlinear. Linear model has contact stiffness kC = 7.7-107 N/m and contact damping cC = 3.2-103 Ns/m and nonlinear model has contact toughness kC = 1.8-108 Pa/m12 and parameter a = 3 s/m. The output inverter frequency before unbalance appearance was set to nsf = 40 Hz. It turns out that the frictional losses in test rig are quite big so the rotor decelerates, with small unbalance and without contact, in some 5 to 6 seconds from 40 Hz to standstill. In Figs. 10 to 12 the rotor response without contact to stator is presented, without any additional unbalance, madd = 0 (or presented via rotor mass eccentricity eadd = 0). If someone compares Fig. 9 and Fig. 12, except differences in the shape of first order (linear in Fig. 9 and nonlinear in Fig. 12), on Fig. 12 exists additional horizontal line at 28 Hz. The reason for this mainly lies in fast rotor speed change while passing though its critical speed. -0.5 — -0.5 0 , mm 0.5 Fig. 11. Simulated lateral rotor response in horizontal direction for free RD analysis with speed law eadd = 0 Fig. 12. Spectral map of rotor displacements in horizontal direction, free RD, eadd = 0 Fig. 10. Rotor center orbit during free RD (red dashed line represents a clearance circle) Fig.13. Rotor center orbit during RD, rotorstator contact appeared, eadd = 0.124 mm Fig. 14. Lateral rotor and stator responses in horizontal direction, contact appeared, eadd = 1.24-10-4 m, rr0 = 0.05 mm, a0 = n rad Fig. 15. Rotor speed laws for simulations without and with contact 200 150 I 100 Fig. 17. Spectral map of rotor relative torsional angles qm - qr, contact appeared 200 150 % 100 50 Fig. 18. Spectral map of stator lateral velocities in horizontal direction, contact appeared 200 150 I 100 Fig. 16. Spectral map of rotor displacement in horizontal direction, contact appeared 200 150 % 100 Fig.19. Spectral map of stator torsional angles, contact appeared Figures 13 to 19 show rotor and stator response after the appearance of the additional mass unbalance (presented in rotor mass eccentricity) eadd = 0.124 mm, so the total rotor mass eccentricity is equal to e = em + eadd =1.17 mm. The magnitude of the vector rr0 was rr0 = 0.05 mm and the phase lag was a0 = n rad. In this simulation the linear model of normal contact force has been applied and coefficient of friction in contact was /uc = 0.18, according to Bartha [2]. In Fig. 11 we had the supercritical self-balancing effect because the mass eccentricity e was smaller than the magnitude of the shaft bow vector ro. On the contrary in Fig. 14 the self balancing effect is placed subcritically at the t1 = 3.8 s of simulation because after sudden unbalance increase, the total mass unbalance eccentricity e became greater than rr0. Fig. 15 shows the difference between speed laws for simulations with and without the rotor-stator contact appearance. As expected a rotor deceleration for simulation where the rotorstator contact appeared, was more intensive especially in a time period with established rotorstator contact (see contact indication in Fig. 15). After introduction of SRUI, the rotor made several impacts to the stator and then a permanent contact was established. In Figs. 16 to 19 there is a narrow time period around t1 = 1 s of simulation where multiple harmonics can be seen. This is due to intermittent rotor-stator contact. It is interesting to see after the rotor-stator separation (ti = 2.75 s) that the rotor and stator continued to vibrate with their own flexural natural frequencies, f = 28 Hz (Fig. 16), f = 90 Hz (Fig. 18) and torsional frequencies ft = 137.9 Hz (Fig. 17), ft = 97.7 Hz (Fig. 19). In Fig. 20, the influence of the phase lag of the mass unbalance eccentricity vector e in relation to the static equilibrium position of the bowed shaft rr0 i.e. a0 has been analyzed. Presented contact forces are related to a0 = 180°, 135° and 90°. In the same time normal contact forces of the linear and nonlinear models have been compared, so the diagrams in the left column correspond to the linear and in the right column to the nonlinear model of the normal contact force. The main difference between the linear and nonlinear model of the contact normal force, presented in Fig. 20, can be seen in the region immediately after appearance of SRUI. Nonlinear model generally had a greater force response for the first impact, while in permanent contact both linear and nonlinear model had almost the same values. Details about the differences between analyzed normal force models and influence of the phase lag a0 are presented in Table 1. The normal force for the permanent rotorstator contact (with sliding) FCperm has showed expected increasing trend with decreasing phase lag a0, because a0 smaller than 180° means greater Table 1. Results of contact normal force analysis regarding phase lag of the mass eccentricity vector e in relation to shaft bow radius vector rr0 «0, fcmax / fcperm , N Linear Nonlinear 180° 327.7 / 96.8 611.9 / 97.1 135° 279.8 / 99.3 462.2 / 99.5 90° 203.7 / 104.8 277.1 / 105.0 unbalance excitation. The reason for contradictory trend of decreasing of the FCmax (max. value of the normal force for the first impact) with decreasing a0 (increasing unbalance excitation) lies in the fact that a0 of the em and eadd is the same. This caused an increasing tendency of the radius of initial rotor orbit and therefore the decrease of the maximum normal velocity at the first rotor-stator impact. 3 EXPERIMENTAL ANALYSIS The analysis has been performed on the test rig which can be seen on Fig. 21. The test rig has been specifically designed for the rotor-stator contact measurements. It consists of the following main components: robust foundation, mounting plate, roller bearings with their supports, rotor, stator, flexible coupling, induction motor with speed controller, measuring system based on Bruel and Kjaer non-contacting displacement sensors and National Instruments data acquisition card PCI NI 4472 with adequate software (LabVIEW, Matlab). Measurements showed that the first critical speed of the rotor was 28 Hz, while the first natural frequency of the stator was 90 Hz. Radial clearance between the rotor and stator was 0.4 mm. Before the experiments the rotor had been balanced. a) b) c) Linear model Non-linear model Fig. 20. Normal forces in contact between rotor and stator; left column - linear model, right column nonlinear model; a) a0 = 180°, b) a0 = 135°, c) a0 = 90° Also before every single experiment the contact surfaces of the rotor and stator were lubricated with WD40 spray [2] to decrease the coefficient of friction to the value of approximately ju = 0.18 and to minimize the possibility of destruction of sensors and other structural parts of the test rig. 3.1 Experimental, Motor Controlled, RunDown (RD) Analysis To verify the numerical results presented in chapter 2.1, an experimental RD analysis has been performed. Fig.21. Test rig for rotor-stator contact investigation The speed change has been controlled via frequency inverter Micromaster 440 from Siemens. The speed rate has been the same as in the numerical simulation only the starting speed has been 70 Hz. Figs. 22 and 23 show measured lateral responses of the rotor in horizontal x direction with speed law and spectral map of the rotor lateral response in the horizontal direction. This figures can be compared to Figs. 8 and 9 respectively, for a speed range n = 40 to 0 Hz. In Fig. 23 except for the first harmonic, other higher harmonics can be also seen due to various imperfections of the experimental model like radial and angular misalignment etc. Fig. 24 shows a spectral map of stator displacements in horizontal direction. Although there wasn't direct contact between the rotor and stator, the stator was excited as it can be seen on Fig. 24. This is due to the fact that rotor bearing supports and stator supports are rigidly connected on the mounting plate, so this situation allows an indirect stator excitation through the mounting plate. Measured flexural natural frequencies were, fr = 28 Hz (Fig. 23), fs = 90 Hz (Fig. 24) and stator torsional frequency fst = 102.5 Hz (Fig. 24). To excite effectively both stator natural frequencies (with second harmonic), the 70 Hz as a starting speed of this RD analysis was chosen. 0.5 80 60 N -L 40 c 20 0 it 0 10 20 30 40 50 60 30 t. s Fig.22. Measured lateral rotor xr response in horizontal direction and speed law . ■ . '■■■■ , v-'.V Awl Mi Fig.23. Measured spectral map of rotor displacement in horizontal direction I HI ■ .... ' . , v': ': ■ VSTy .-