Investigations on bending-torsional vibrations of rotor during rotor-stator rub using Lagrange multiplier method

Journal of Sound and Vibration 401 (2017) 94–113

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Investigations on bending-torsional vibrations of rotor during rotor-stator rub using Lagrange multiplier method Md Asjad Mokhtar, Ashish Kamalakar Darpe n, Kshitij Gupta Department of Mechanical Engineering, Indian Institute of Technology Delhi, New Delhi 110016, India

a r t i c l e in f o

abstract

Article history: Received 10 September 2016 Received in revised form 22 March 2017 Accepted 26 March 2017

The ever-increasing need of highly efficient rotating machinery causes reduction in the clearance between rotating and non-rotating parts and increase in the chances of interaction between these parts. The rotor-stator contact, known as rub, has always been recognized as one of the potential causes of rotor system malfunctions and a source of secondary failures. It is one of few causes that influence both lateral and torsional vibrations. In this paper, the rotor stator interaction phenomenon is investigated in the finite element framework using Lagrange multiplier based contact mechanics approach. The stator is modelled as a beam that can respond to axial penetration and lateral friction force during the contact with the rotor. It ensures dynamic stator contact boundary and more realistic contact conditions in contrast to most of the earlier approaches. The rotor bending-torsional mode coupling during contact is considered and the vibration response in bending and torsion are analysed. The effect of parameters such as clearance, friction coefficient and stator stiffness are studied at various operating speeds and it has been found that certain parameter values generate peculiar rub related features. Presence of sub-harmonics in the lateral vibration frequency spectra are prominently observed when the rotor operates near the integer multiple of its lateral critical speed. The spectrum cascade of torsional vibration shows the presence of bending critical speed along with the larger amplitudes of frequencies close to torsional natural frequency of the rotor. When 1 m × n X frequency component of rotational frequency comes closer to the torsional natural frequency, stronger torsional vibration amplitude is noticed in the spectrum cascade. The combined information from the stator vibration and rotor lateral-torsional vibration spectral features is proposed for robust rub identification. & 2017 Published by Elsevier Ltd.

Keywords: Rotor-stator rub Torsional vibration Lagrange multiplier method Contact mechanics

1. Introduction Among various research issues, fault diagnosis for a rotating machine is an important topic. Once the rotor bearing system is properly designed, its continuous performance without untimely failure can be ensured by way of monitoring its health during service. Among various faults, the rotor stator rub has always been of great significance, particularly in high speed flexible rotors. Reducing the clearance for efficiency increases the undesired rotor stator contact possibilities. The rub is termed as the contact between the rotor and the stator that leads to significantly changed rotordynamic behaviour. It is one of very few causes that influence both lateral and torsional rotordynamics. The detection of torsional vibrations is also n

Corresponding author. E-mail addresses: [email protected] (M.A. Mokhtar), [email protected] (A. Kamalakar Darpe), [email protected] (K. Gupta).

http://dx.doi.org/10.1016/j.jsv.2017.03.026 0022-460X/& 2017 Published by Elsevier Ltd.

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Nomenclature aT Ai cg Cr, Cs Dr Es f g gN , gN

gṪ , ġ T Gr GL GN, GT GNT Kr, Ks Mr, Ms nT Nr, Ns pNi

unit vector along tangential direction contact area for the ith nodal pair ¼ g/max(X(ω)), clearance ratio effective damping matrices: rotor, stator damping matrix for rotor elastic modulus of stator global force vector of body forces and surface traction gap or clearance between contacting bodies normal gap function and the array of all the normal gap functions time derivative of tangential gap function and the array of all the gṪ gyroscopic matrix for rotor local contact matrix for nodal contact pair global contact matrix for normal and tangential direction ¼ GN þ mGT, global normal-tangential contact matrix stiffness matrices: rotor, stator mass matrices: rotor, stator unit vector along normal direction degree of freedom: rotor, stator normal contact pressure for ith contact pair

95

PA, PB position vector for point A, B q number of active nodal contact pair r ¼ ω/ωL0, speed ratio rd radius of the disc Rr (t), Rs (t) external force vectors scalar form of tangential surface traction tTi for i th contacting pair displacement update without contact u * uc corrective displacement vector ur, u̇ r, ür displacement, velocity and acceleration vectors for rotor us, u̇ s, üs displacement, velocity and acceleration vectors for stator X° vector defining the initial gap between contacting nodal pairs X(ω) response as function of speed δ□ variation of the parameter □ normal and tangential contact force λN , λT normal and tangential contact force λN , λT vectors m friction coefficient Ω speed of rotation ωL0, ωT 0, ωS0 1st natural frequency: rotor's lateral, torsional and stator's lateral direction Π Potential energy Energy due to contact contribution ΠCLM

important in practical rotors where unnoticed large torsional vibrations can lead to fatigue induced cracks and failure. The rotor stator rub has been under investigations for over three decades. Childs [1,2] in his early works has shown that the rubbing causes parametric excitation of half speed whirl at a rotor's natural frequency. He showed increase in radial stiffness and tangential force due to rubbing. In subsequent studies, third order sub-synchronous motion is also found due to partial rubbing. During rotor operation under rub, Beatty [3] emphasized the presence of second and third harmonics of the synchronous frequency. He also concluded that the frictional cross-coupling term may lead to instability during full contact rub. Choy and Padovan [4] investigated the event of rub as a nonlinear phenomenon and studied the transient motion of the rotor due to rub. The process is divided into four stages as non-contact stage, rub initiation, rub interaction and separation. They studied effects of casing stiffness, friction co-efficient, unbalance load and system damping. Muszynska [5] has provided a detailed review of the literature on the theoretical and experimental works. Muszynska and Goldman [6] in their analytical, numerical, and experimental studies, have concluded that the vibrational behaviour can be characterized 1

by regular periodic vibrations of synchronous (1X), and sub synchronous ( 2 X ,

1 X ...) 3

orders. Chaotic vibration patterns were

found accompanied by higher harmonics. Further, it was observed that chaotic vibration zones decrease with increase in damping. In a study of relatively new problem of wind-milling imbalance in aero engines, Groll and Ewins [7] found very rich frequency spectrum consisting of super-harmonics upto order of 9 and sub-harmonics of the order of 32. Nonlinear vibrations related to rub have been studied by many researchers. Chu and Zhang [8] have used non-linear model with piecewise linear stiffness and have shown that rub impact between rotor and stator exhibits periodic, quasiperiodic and chaotic vibrations. In a similar study, Sun et al. [9] have shown that motion of the rotor system alternates among the periodic, chaotic and quasi-periodic vibrations, as the rotating speed increases. It was concluded that the proper increase of damping can make the chaotic motion return to periodic one, which is contrary to the observations made in [8]. Chu and Lu [10] have reported an experimental work to investigate nonlinear vibrations in a rub impact rotor system and observed a very rich form of periodic and chaotic vibrations with presence of fractional harmonic components such as 1 3 5 X , 2 X , 2 X , etc. Most of the work on the rub is based on conventional friction model. Only a few attempts are made to 2 address the rub phenomenon in an alternate way such as impact energy model as reported by Cong et al. [11]. The stiffening effect due to rub is observed and analysed by Chu and Lu [12] and Patel and Darpe [13]. Dynamic stiffness identification for determining rub location was studied by Chu and Lu [14]. Patel and Darpe [13] have shown the possibility of detecting rub at its initiation stage. They also investigated in detail the directional nature of rub fault along with the shift in resonance speed. Advanced signal processing tools have been used for time-frequency-energy mapping to estimate the exact time of occurrence of rub and associated excitation frequencies. Similarly, most of the rotordynamic models for rub investigations have considered lateral vibrations only. Presence of tangential force due to contact and friction leads to torsional vibrations in the system. Edwards et al. [15] have investigated

96

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

the torsional vibration response due to rubbing and showed that the system is more stable with torsional effect considered. Later, Sun et al. [16] have found that the periodic and chaotic responses are comparable with and without the bendingtorsion coupling. Al-Bedoor [17,18] had carried out a transient analysis with torsional degree of freedom and found irregular rubbing orbits and split in response at resonance. He observed super-harmonics of order 3, with possible instability, in case of torsional vibrations. Mohiuddin and Khulief [19] proposed an FE formulation using Lagrangian approach for coupled bending-torsional model which provides an efficient tool for dynamic analysis of a complex rotor system. Patel and Darpe [20] proposed a model for the coupled torsional and lateral vibrations of unbalanced rotors that accounts for the rotor-stator rub along with the crack. A nonlinear dynamic study is done by Khanlo et al. [21] showing disappearance of quasi periodic or chaotic behaviour beyond certain speed ratio when bending-torsional coupling is included in the system. Lately, contact mechanics based approach dealing with two or more bodies in contact is considered for rotor stator interaction. This approach has an advantage of modelling the stator as a part of the rotor bearing system and thus includes its dynamic effect on the rotor response. Roques et al. [22] have used this approach for constraint enforcement and did a transient analysis for the rotor stator interaction caused by accidental blade-off imbalance in a turbo-generator. Sensitivity analysis done in this work proves the significance of friction coefficient and the stator flexibility. Behzad et al. [23] investigated a Jeffcott rotor system under partial rub condition using contact mechanics formulation. The formulation techniques used in [22,23] are based on Lagrange multiplier method. Ma et al. [24] used an FE model of two disc rotor bearing system and considered augmented Lagrangian method for contact formulation. They presented a qualitative comparison of the response with respect to speed, gap and contact stiffness. Subsequently, Ma et al. [25,26] carried out similar studies with two different loading conditions and varied stator configuration. A comparative study of various constraint enforcement techniques on rotor stator rub simulation is done by Mokhtar et al. [27]. A literature review highlighting the phenomenology involved during the blade-casing contact and the rotor-stator contact is presented by Jacquet-Richardet et al. [28]. Various modelling approaches, the major experimental facilities and the results are summarised. They have delved on the complexity of the problem and identification of major influencing parameters on the dynamics. The recent works [23–27] with contact mechanics based approaches have considered only lateral motions of the rotor system. The torsional vibration is not given its due when it comes to investigations on rotor-stator rub phenomenon and related diagnostic feature extraction. It seems worthwhile to investigate the torsional vibrations of such contacting rotor with flexible stator including the stator dynamics that, to the best of author's knowledge, was neglected in the past. The paper also establishes the correlation between lateral and torsional vibration response in terms of frequency contents of the vibratory signal. The present work uses an FE model of a single disc rotor system and the Lagrange multiplier method as constraint enforcement technique. A node to node contact strategy is used as it accounts for the interaction in case of local deformation at the contact location. The dynamic motion of the stator is considered by modelling it as a flexible beam that can respond to both axial and lateral forces during contact. Most of the previous work considered only axially deformable stators. In addition, the conventional approach of modelling the stator as stationary spring element has a limitation that the dynamics of stator element is not considered. In fact, the contact between the rotor and stator depends not only on the rotor excursion but also on the instantaneous stator vibratory response. Thus, the dynamic stator response affects the rotor stator interaction and must be included in the mathematical model of rotor-stator contact system. Further, an approach of using the stator vibratory response for rub diagnosis is explored in the present work, instead of only relying on the rotor lateral and torsional vibration signature. The next section gives details of the formulation of governing equations along with the contact model and gap function definition, while the results are discussed in Section 3. Section 4 highlights conclusions of the work.

2. Formulation of governing equations The mathematical formulation of the rotor bearing system along with the stator is briefed in this section. The incorporation of constraint imposed by stator on the rotor while formulating the governing equations is also discussed. 2.1. Rotor stator model The problem taken up for the rub simulation consists of simple rotor bearing system with a disc centrally located and the bearings at the two ends of the shaft as shown in Fig. 1. The rotor shaft is discretized into 14 finite beam elements with each node consisting of three translational and three rotational degrees of freedom. The finite element includes effect of rotary inertia, gyroscopic moment and transverse shear effect for the shaft. A proportional damping is considered in the system

Stator y

Disc B1

Rotor shaft

B2

y x

Fig. 1. Rotor bearing stator system.

z

Gap

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

97

and the model assumes material linearity. The two bearings that support the rotating shaft are modelled using orthogonal spring elements and are included in the global stiffness matrix at its appropriate location. The equations of motion for the rotor system are written as

⎡⎣ M r ⎤⎦{ u¨ r } + ⎡⎣ C r ⎤⎦{ u̇ r } + ⎡⎣ K r⎤⎦{ u r } = { R r(t )} where

⎫ ⎬ ⎭

⎪

⎪

C r = Dr + ωGr

(1)

where Mr, Dr, Gr and Kr are mass, damping, gyroscopic and stiffness matrices for rotor. Rr and ur are generalized force and displacement vectors. The matrix Kr includes the support stiffness of the bearing and Cr is the effective damping matrix for rotor. Similarly, the stator is modelled as discretised beam with its inertial and elastic properties and placed vertically with a specified clearance with the disc. The equations of motion for the stator are written as

⎡⎣ Ms⎤⎦{ u¨ s } + ⎡⎣ Cs⎤⎦{ u̇ s } + ⎡⎣ K s⎤⎦{ us } = { R s(t )},

(2)

where subscript ‘s’ stands for stator. 2.2. Contact model For modelling the interaction of rotor and stator parts, Lagrange multiplier method is used for constraint enforcement. The method ensures perfect enforcement of constraint at the rotor stator contact. During the contact, the Lagrange multiplier is introduced and the governing equation of motion is added with an additional term that ensures to satisfy the constraints in the system. Modelling of the contact dynamics requires combining the rotor and stator equations of motion. The combined rotor stator system will act as a coupled system during the interaction and uncoupled otherwise. An appropriate gap function is defined that gives the instantaneous gap between the two contacting bodies. The gap function defined in normal direction leads to the contact condition that is the constraint for the system to ensure nonpenetration. Equations of motion are derived using variational mechanics, principle of virtual work, and Lagrange multiplier method [29]. The governing equations are then solved through a numerical integration technique. 2.2.1. Gap function definition For the two bodies coming into contact, the normal and tangential gap function between them needs to be defined. Consider points A and B as one nodal contacting pair (Fig. 2). PA and PB are position vectors of nodal points on the two contacting bodies. The corresponding normal gap function can be written as [23]

⎯⎯⎯→ gN = nT ( PB−PA) = nT AB

( )

(3)

where gN is the normal gap function, nT is the unit vector along AB. The matrix form for the gap function is

gN = ⎡⎣ −nx −ny nx

⎡ xA⎤ ⎡ xA⎤ ⎢ ⎥ ⎢ ⎥ A ⎢ y A⎥ ⎢y ⎥ ny ⎤⎦⎢ ⎥ = ⎡⎣ GL⎤⎦ ⎢ ⎥ 1×4 B B ⎢x ⎥ ⎢x ⎥ ⎢⎣ y B ⎥⎦ ⎢⎣ y B ⎥⎦

(4) A

B

A

B

GL is the local contact matrix for a pair of contacting nodes. x , x , y , y are the components of position vectors PA and PB.

y B a

PA

n

A

g

PB x Fig. 2. Gap function definition.

98

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Depending upon the other possible locations of interaction in the model, one can generate corresponding local contact matrices and then it can be assembled into a global contact matrix which has the information of all the nodal pairs in contact. With Nr and Ns as number of the degrees of freedom for rotor and stator respectively and q as the no. of active nodal contact pairs, the gap function can be generalized as [23]

(

gN(u) = GN X 0 + u

)

(5)

where T T T⎤ ⎡ T ⎡⎣ GN⎤⎦ = ⎣ G1L G2L ... G q ⎦ q × 2( Nr + Ns) L

X0 is the vector defining the initial gap between contacting nodal pairs at t¼0, u the displacement vector and the sum, X þu is spatial coordinate vector. Similarly, for relative tangential velocity, the time derivative of tangential gap function will be 0

gṪ = aT ( u̇ B −u̇ A),

hence

ġ T = G T u̇

(6)

where a is the unit vector in tangential direction, gṪ is the time derivative of tangential gap function ( gT ) and ġ T is the array of tangential gap functions. GT is global tangential contact matrix and u̇ corresponds to velocity vector. In classical contact mechanics, it is assumed that the reaction force between the contacting bodies is always compressive in nature making it either negative or zero. Similarly, the gap between contacting bodies are either positive or zero. It can be expressed mathematically as Hertz-Signorini-Moreau conditions [29]

⎧ gi ≥ 0 ⎪ N ⎪ i ⎨ p ≤ 0 for i = 1, 2.....q ⎪ N i i ⎪ ⎩ pN gN = 0

(7)

where pN is the normal contact pressure. The average rotor angular velocity remains constant as the energy dissipated due to friction is assumed to be compensated by the driving torque. Thus the rotor boundary steadily slides over the stator boundary. The Coulomb's law is used to consider the frictional effects into the constrained equations of motion. Thus, tangential surface traction (tT) is directly proportional to the normal reaction and the scalar form can be written as

tTi = μpNi

(8)

If the contact area for ith nodal contact pair is Ai

⎧ ⎪ λNi = pNi Ai ⎨ for i = 1, 2.....q i i i ⎪ ⎩ λT = tT A

(9)

where λN and λT are the normal and tangential contact force respectively. Eqs. (7) and (8) can be rewritten in terms of contact forces as

⎧ gi ≥ 0 N ⎪ ⎪ ⎨ λNi ≤ 0 and λTi = μλNi , ⎪ i i ⎪ ⎩ λN g N = 0

for i = 1, 2....q (10)

2.2.2. Constrained equation of motion The Lagrange multiplier is used to add constraints to a weak form of the equilibrium equation during interaction. The energy contributions related to interaction between rotor and casing is [29]

ΠCLM =

∫Γ ( pN gN + tT gT )dA

(11)

c

Variation of this contact contribution will be

δΠCLM =

∫Γ ( pN δgN + tT δgT )dA + ∫Γ ( δpN gN + δtT gT )dA c

c

(12)

In the Eq. (12), the first integral is associated with the virtual work of the Lagrange multipliers along the variation of the gap functions in normal and tangential directions. The second integral describes the enforcement of the constraints [29]. δgN is the variation of the normal gap function. The terms tT.δgT and δtT.gT are associated with tangential stick. In the case of

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

99

sliding, a tangential stress vector tT is determined from the constitutive law for frictional slip. Therefore, the normal contact pressure, pN acts as the sole Lagrange multiplier. The term related to variation of tangential surface traction may disappear; as it is not required to be zero in case of slip motion. Using node to node contact strategy in FE discretization, Eq. (12) is simplified for the q activated nodal contact pairs and can be written in matrix form as [23]

δΠCLM =

∫Γ ( pN δgN + tT ⋅δgT )dA + ∫Γ δpN gN dA = δgNT λN + δgTT λT + ( δλN )T gN c

where λN =

[λN1 ,

c

λN2......,

λNq ]T ,

λT =

[λT1,

λT2......,

λTq]T

and gN =

[gN1 ,

gN2 ......,

(13)

gNq ]T

In terms of virtual displacements (δu), T

δΠCLM = δ uT GTNλN + δ uT G TT λT + ( δ λN ) gN = δ uT GT λ + δ λTN gN

(14)

where λT = [λTN λTT ] and GT = [GTN G TT ] By applying FE discretization procedure to the linear elastic bodies of the rotor and stator, the total potential energy due to static deformation is given by

Π (total) = Π (system) + ΠC (contact )

(15)

Taking variation of Eq. (15) with respect to displacements gives the weak form of the equilibrium equation.

δΠ LM = 0 i. e . T

T

δΠ + δΠCLM = 0 , T

T

δ u Ku−δ u f + δ u G λ +

δ λTN gN

leads to

=0

(16)

where u is the vector of nodal displacements, K is the stiffness matrix of rotor-stator combined system and f is the vector of body forces and surface tractions. These variables contain information of both the rotor and stator in which deformations of rotor and stator are not coupled. Eq. (16) must hold for any arbitrary value of δu and the unilateral contact constraints (Eq. (10)) can be expressed as complementarity condition such that the contact force vector is used in expressing the governing equations of motion as [30]

⎧ T ⎪ Ku + G λ = f ⎨ ⎪ ⎩ gN ≥ 0 ; λN ≤ 0 ;

λNi gNi = 0

for i = 1, .....q

(17)

The governing equation of motion is second order in time with constraints holding on the displacement only. Due to unilateral constraint, the velocity can be discontinuous, and there is in general no strong solution to this problem [31]. Hence for the solution to be unique, the complementarity conditions should be supplemented with an impact law governing the dynamics during impact [30]. Thus, an additional condition, specifying the velocity after impact is required that is defined as u̇ þ ¼ -e  u̇ - where e is the coefficient of restitution (e ϵ [0 1]) [31]. In a continuous framework problem, where the contact boundary is reduced to a point, uniqueness of a solution do exist without the impact law [[31] and Refs. [11, 28] therein]. The present case looks for a weak solution where displacement u is continuous in time and velocity u̇ , acceleration ü and contact force λ is evaluated (containing impulses). The problem formulated here is to evaluate displacement u and contact force λ. With a node-to-node contact strategy and the matrix K being positive definite, the Eq. (17) can have a unique solution. Now, considering inertia forces and damping forces to study the contact dynamics between rotor and stator, Eq. (17) can be written as

⎧ ⎪ Mu¨ + Cu̇ + Ku + GT λ = f ⎨ 0 ⎪ ⎩ GN u + X −rd ≥ 0 ; λN ≤ 0 ;

(

)

λNi gNi = 0

for i = 1, .....q

(18)

where ü and u̇ are vectors of nodal acceleration and velocity respectively and rd is radius of the disc. This is the constrained equation of motion for the rotor and stator to be solved by explicit numerical integration scheme. The matrix M, C and K are mass, damping and stiffness matrices of the combined rotor-stator system. Using the expression λTi = μλNi , Eq. (18) can be expressed as

⎧ Mu¨ + Cu̇ + Ku + GT λ = f ⎪ NT N ⎨ 0 ⎪ G u X r 0 ; + − ≥ λN ≤ 0 ; N d ⎩

(

)

λNi gNi = 0

for i = 1, .....q

where GNT is the normal-tangential-contact matrix, GNT = GN + μG T

(19)

100

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

2.3. Numerical simulations A good contact treatment depends strongly on the integration method being used for solving the equations of motion. The implicit solutions are good for static and slow transient problems; however it behaves badly when inertial forces are relatively large [32]. Implicit methods with automatic control of time step size may also be used. Conditionally stable explicit methods appear to be more relevant to non-differentiable functions such as contact [22]. The solution procedure adopted here is based on forward increment Lagrange multiplier method with central difference integration scheme. It has been established that the Lagrange-explicit results are accurate enough and could further improve with mesh refinements [32]. The discretization in the present work is based on the wave speed in rotor's lateral and stator's longitudinal direction and the frequency of interest. Further, the convergence study over a no. of element size and time step is carried out and the time step of 0.5 ms is used for the present study. 2.3.1. Solution procedure The central difference scheme is a well-known solution procedure for numerical integration. The steps involved during the solution are as follows. The incremental solution of equations proceeds by first calculating the usual explicit displacement update u n + 1 for all degrees of freedom, including surface contactor nodes, assuming no contact i.e. λN ¼0 [23,32] *

˜ −1f˜ n un+ 1 = M * where ˜ = 1M+ 1 C , M 2h h2

⎛ ⎛ 1 n 2 ⎞ 1 ⎞ n− 1 f˜ = f n − ⎜ K − 2 M⎟un − ⎜ 2 M − C⎟u ⎝ ⎝h 2h ⎠ h ⎠

and

0

(20)

n+1

Following it, the spatial co-ordinates (X + u ) are processed by a search algorithm to identify all surface contactor * nodes that have penetrated the target surfaces. Based on this search information a constraint matrix is assembled and the Lagrange multiplier vector λN can be determined and included into the computations to ensure non-penetration. For actual displacement that ensures non penetration, the explicit displacement has to be corrected using

un + 1 = u n + 1 + u c *

(21) 0

where uc represents the required correction. Adding X and multiplying gives

(

)

(

GnN+ 1

to Eq. (21), then subtracting ‘rd’ from both sides

)

GnN+ 1 X 0 + un + 1 − rd = GnN+ 1 u n + 1 + X 0 − rd + GnN+ 1u c *

(

)

n+ 1

= gN u *

+

GnN+ 1u c

(22)

For the non-penetrability, the Eq. (22) must be equal to zero. Referring to Eq. (19) and Eq. (20), the corrective displacement vector can be evaluated by T ˜ n + 1 = f˜ n= − Gn + 1 λn Mu c c NT N

(

)

(23)

Substituting Eq. (23) into Eq. (22) and the condition of instantaneous gap function to be zero at time tn þ 1, the Lagrange multiplier and the corrective displacement vector can be calculated by following equations [23] −1 T ⎛ ˜ −1 Gn + 1 ⎞⎟ g u n + 1 λnN = ⎜ GnN+ 1M NT ⎝ ⎠ N *

(

)

T

˜ −1 Gn + 1 λn uc = M NT N

(

)

(

)

⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭

(24)

The explicit displacement calculated in the first place can now be updated with corrective displacement (Eq. (21)) and the other parameters are evaluated at time tn þ 1 before moving to next time step.

3. Results and discussion In this section, three aspects of the present work are discussed. Including the lateral motion of stator in addition to the axial one is one aspect. The other related aspect of this improvement of stator model is to explore the torsional vibration feature of the rotor when the problem formulation is based on contact theory. Finally, effect of parameters such as rotational speed, rotor-stator clearance, friction coefficient and the stator flexibility are investigated on the rotor response. In practical cases of rotor stator rub, the stator get excited as a general perturbation from the contacting spinning rotor and may vibrate in all its modes in particular axial and bending. In the past [23–26], most stator models are considered as an axial bar, disregarding any lateral degree of freedom. Such models thus assume only compressive movement of the stator under the

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

101

Fig. 3. Orbit plots for speed ratio, (a) r ¼0.75 (b) r¼ 1.5 (c) r¼ 3.

contact forces. However, the contact force between the rotor and stator does provide significant excitation tangentially due to contact friction. Inclusion of the lateral degree of freedom of stator makes it more realistic and an appropriately modelled stator is expected to correctly and comprehensively represent the dynamics of rotor-stator motion during the interaction. The torsional vibration waveform is discussed along with lateral vibration of the rotor during the interaction. 3.1. Stator model with lateral degree of freedom The contact formulation discussed in the previous sections allows inclusion of the FE model of the stator as a part of the global system. For the validation purpose, the configuration of rotor bearing system used by Behzad et al. [23] is considered and investigated for the vibration response with and without the lateral DOF of the stator. The orbit plots for three speed ratios are shown in Fig. 3. The dashed line corresponds to the one simulated without considering the lateral (bending) degree of freedom for the stator as used in their work [23]. The horizontal lines crossing the orbits represent the initial position of the stator node. While there is an overall agreement in the orbit shape, a variation in the orientation of orbit is noticed between the two models. The difference depends on the speed of rotation; the orbit plot at r ¼1.5 shows significant deviation. The orbits have tilted from vertical for the model with the inclusion of lateral degree of freedom of stator. For the validation of present rotor-stator contact formulation with Lagrange multiplier approach with torsional degree of freedom, unbalance vibration response is computed for the rotor model of Patel et al. [20]. The contact model used in [20] is a conventional friction contact. Torsional oscillations and vertical vibration waveforms are compared for assessing the response using the two different contact formulation techniques (Fig. 4). In terms of amplitudes of vibration response, the vertical responses are in close agreement with each other whereas the torsional vibration amplitudes show some difference. In the conventional model the vertical response waveform (Fig. 4a) has slight variation of the amplitude in successive cycles. The torsional vibration amplitudes are relatively less with the present contact formulation as compared with the conventional model as observed in Fig. 4d. Some difference in consecutive rotor-stator interaction (observed in lateral vibration response shown encircled regions in Fig. 4a) is noticed in the conventional model. It may be mentioned that in the conventional stator model, the stator boundary is assumed to be fixed, while in the present model involving contact formulation, it is dynamic. With the current model and for the given rotor system configuration and parameters such cycle to cycle variation of amplitudes at the interactions is not observed (as noticed in the encircled interaction regions in Fig. 4b). In addition, the response (torsional and vertical) appears to be periodic with rotational frequency, whereas in the case of conventional model, it is periodic with twice the period of a rotation for the given operating parameters and rotor bearing stator system configuration. 3.2. Simulation of rotor-stator interaction Having compared the rotor lateral and torsional vibration response, using the rotor bearing system configurations of respective past work, the investigations are now focussed on the rotor bearing stator system described in Section 2.1. The material and geometric properties of the rotor stator system along with FE information are given in Table 1. The external force vector ‘Rr’ (Eq. (1)) mainly constitutes an unbalance excitation applied at the disc location. The resulting rotor deflection once exceeds the specified clearance, induces rub that is studied over a wide range of rotational speed. The lateral rotor vibrations have been investigated extensively in the past, but the torsional vibrations are relatively less explored. Using Lagrange multiplier based contact formulation technique, the torsional vibration features in addition to the rotor lateral vibration for a rotor-stator system is now investigated. The rotor and stator vibration responses obtained are analysed using vibration waveforms of rotor and stator, orbit plots of shaft centre, time frequency plot, contact reaction force and corresponding frequency spectrum. Simulations are carried out over the speed range of r ¼0.75 to 6 and the orbit plots are observed. As the operating speed changes, the shape and size of the orbit is found to change as shown in Fig. 5. It has been shown in the past that the critical speed with large amplitude occurs well beyond the theoretical resonance speed without stator contact. This is due to the

102

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Fig. 4. Vertical and torsional response (a, c) conventional friction model [20] (b, d) present contact model. Table 1 Geometrical and material properties of the rotor stator system. Shaft

Disc

Stator

Unbalance

Length

80 cm

Diameter Elastic modulus Shear modulus Density Diameter/thickness Density Disc mass Elastic modulus Length Area of c/s Mass Radius

2 cm 210 GPa 79.6 GPa 7800kg-m  3 12 cm/3 cm 7800kg-m  3 2.64 kg 3 GPa (Es) 4 cm 25 mm2 0.008 kg 4.5 cm

1e6 N-m-1

Bearing linear translational stiffness Friction coefficient, m *Clearance ratio, cg Speed ratio, r¼ ω/ωL0

0.3 0.70 & 0.85 0.75, 1.5, 2.0 and 3.0

Natural frequency

Rotor

No. of elements

Stator Rotor Stator

31.56 Hz (ωL0, Lateral) 128 Hz (ωT0, Torsional) 803 Hz (ωS0, Lateral) 14 6

Disc Bearings

8 1, 15 (B1 and B2)

Nodal location

n Clearance ratio is defined as cg ¼ g/max(X(ω)), where g is the radial gap between rotor and stator at stationary condition. X(ω) represents response of the rotor (for given rotor configuration) as a function of speed ‘ω’ in absence of rotor-stator contact. The value of cg thus represents the fraction of allowable rotor motion without contacting the stator. For the rotor stator interaction to happen, the value of cg must be less than 1 which means that the initial gap ‘g’ specified should be less than the possible rotor deflection without any constraint.

fact that the rotor gets supported by the compliant stator at the disc location and hence has higher effective bending natural frequency. The resonance at operating speed close to this modified higher natural frequency is referred to as pseudo critical speed. At speed ratio r ¼1.7 the orbit shape is observed to attain maximum size and have a single loop which corresponds to 1 the vibration synchronous to its pseudo critical speed. In addition, at higher spin speed of r ¼3.5, the 2 X component due to 1 rotor stator rub, generates resonance in the rotor stator system. Similarly, at speed r ¼5.3, the 3 X component generates resonance which is evident from the change of orbital response of the rotor across this speed. Comparison of these orbital response with the experimental results reported by Muszynska (Fig. 5.6.8 in [5]) and theoretical results reported by Ehrich (Fig. 10(b) in [33]), shows that the present results are in good agreement. Fig. 6 shows the vibration response for speed ratio of r ¼0.75 and the rotor-stator clearance of cg ¼0.85. Recall that the cg represents the fraction of allowable rotor motion without contacting the stator. The vertical response as shown in Fig. 6a repeats itself after every two cycles. A close look at the vertical vibration response reveals significant deviation in the period 1

3

between successive peaks. This leads to the presence of 2 X and 2 X combined with fundamental period of 2 T (where T is time for a cycle of rotation). The orbit plot (Fig. 6h) shows two different lobes for successive cycles and repeats itself

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

103

Fig. 5. Orbit plots for speed ratio near (a) r ¼ 1.7 (b) r ¼ 3.5 (c) r ¼5.3.

afterwards. This kind of orbit is generally found in low speed rotor operation with rub. A similar feature is observed in horizontal response as shown in Fig. 6c where the repetition of two successive cycles is present. The corresponding fre1

3

3

quency spectra shown in Figs. 6b and 6d depicts 2 X , 2 X and 2X in vertical and 2 X component in horizontal response in addition to synchronous frequency component. Frequency domain parameters are evaluated with large number of cycles for better resolution. It has been ensured that the initial transient part of the simulation response is ignored and the steady state response is used to find the spectral content. The torsional vibration time domain waveform and its frequency spectrum are shown in Figs. 6e and 6f. The frictional contact between rotor and stator induces the torsional excitation to the rotor, which can be observed in the time domain torsional vibration response. The torsional vibration signal decays over time until another contact of the rotor with stator generates next excitation through a transient short duration contact. It is important to note that there is cycle to cycle variation in the level of contact reaction force generated between stator and rotor (as noticed later in Fig. 7a). This variation in turn generates cycle to cycle variation in the torsional excitation and corresponding torsional vibration response. For example, the engagement of rotor with stator at t¼ 7.855 s results in higher contact reaction force (Fig. 7a) than that generated during next engagement cycle at t¼7.901 s. This results in proportional frictional torque on the rotor and corresponding different level of torsional vibration amplitude immediately after the end of respective contacts (peak amplitude of 7.31e-5 rad after t¼7.855 s as against peak amplitude of 6.101e-5 rad after t¼ 7.901 s as noticed in Fig. 6e). This cyclic 1

modulation in amplitude reflects in the presence of 2 X component in the frequency spectrum of torsional vibration (Fig. 6f).

Fig. 6. Time response and frequency spectrum plots for r ¼0.75, Es ¼ 3e9, m ¼ 0.3, cg ¼ 0.85 (g ¼ 100 mm); (a, b) Vertical, (c, d) Horizontal, (e, f) Torsional, (g) Torsional acceleration (h) Orbit.

104

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Fig. 7. Time response plots for r ¼0.75, Es ¼ 3e9, m ¼ 0.3, cg ¼0.85; (a) Contact reaction force, (b) Stator lateral deflection and (c) Closeup view of (b) for t ¼7.985–7.995.

The dominant frequency in the torsional vibration spectrum is 130.2 Hz which is 11th harmonic of half the rotational 11

frequency ( 2 X ) and closest to the torsional natural frequency of the rotor (128 Hz). As the rotor contacts the stator, it is expected that there would be sudden change in the torsional response at the moment when contact is established and when it is lost. The torsional acceleration response (Fig. 6g) shows substantial change in the slope of the waveform indicating the instant of contact that is not clearly visible in torsional displacement plot (Fig. 6e). A closeup view of time domain plot itself shows the start and end of engagement. The frequency spectrum of torsional acceleration is found to have similar frequency content as observed for torsional displacement (not shown in Fig. 6). The duration of rotor stator contact for this case is 3 ms, which can also be observed by the duration of non-zero contact reaction force shown in Fig. 7a. The duration of the contact varies with different set of system parameters and configurations, as will be shown through various simulations later in Section 3.3. The stator deforms in axial and lateral directions due to contact between the rotor and stator. The axial deflection pattern of the stator is proportional to the reaction force that is shown in Fig. 7a. The plot shows the increase and decrease of the contact force during the interaction and remains zero otherwise. The amplitude of contact reaction force gives an idea about the extent of interaction or the intensity of impact. Fig. 7b shows the lateral vibration response of the stator beam where the amplitude of lateral deflection is initially high due to impact. The response amplitude quickly decays as free vibration response, after disengaging from the rotor until the next impact occurs. The frequency of vibration is observed to be the natural frequency of the stator in lateral mode, which is 803 Hz (ωs0). A closeup view of one interaction span is shown for the stator lateral response in Fig. 7c. The instant of start and end of engagement is also pointed out. The frequency spectrum of stator response shows harmonics of rotational frequency as the impulsive contact is repeated every cycle. The torsional vibrations are repeated transient perturbations due to impulsive excitations from rotor stator contact. Once the contact is lost, the ensuing torsional vibration response is a free vibration response at torsional natural frequency. Due to contact conditions and engagement of rotor with stator, there is a localised deviation in harmonic motion of the torsional vibration that results in generation of spurious frequencies in the torsional vibration spectra and the torsional natural frequency does not appear in the spectrum. Due to transients during the contact, the signal becomes non-stationary and in such cases, the time-frequency domain analysis is a very useful way in revealing temporal variation in frequency composition. The time-frequency domain analysis using Hilbert Huang transform (HHT) is carried out [34]. The instantaneous frequency is obtained by differentiating the analytical signal from Hilbert transform using first order difference method (forward differencing). Using the start and end instances of contact duration of the rotor motion, the torsional response data has been windowed using Hanning window for the non-engaged duration. The engaged duration in the response has not been included for windowing. The windowing helps in isolating the effect of localised deviations in time domain data during the engagement period. The processed time domain torsional vibration response is shown in Fig. 8a, which has been used to find the HHT. The Hilbert spectrum corresponding to the first IMF is plotted as it is expected to capture the higher frequencies present in the response. The darker spots in the map represents signal with high energy level. The energy distribution corresponding to the first IMF shown in Fig. 8b indicates the consistent presence of torsional natural frequency (ωT0 ¼128 Hz) during its non-engaged condition. Hilbert transform is also applied to the stator lateral response and the corresponding HHT spectrum for first IMF is shown in Fig. 8c. The ‘1–2’ and ‘4–5’ portion of plot is the contact duration which has higher energy distribution across wider range of frequencies. The spans ‘2–4’ and ‘5–7’ represent the nonengaged condition of rotor and stator. The frequency-time-energy distribution shows distinct presence of stator lateral natural frequency (803 Hz) for the non-engaged duration of the stator in lateral direction (as shown in spans ‘2–3’ and ‘5–6’ in Fig. 8c), once it is excited at the end of engagement of the rotor and stator (at ‘2’ and ‘5’). The torsional vibration response is found for three different cases of rotational speed such that one of the harmonics (6th in these cases) is close to the torsional natural frequency, ωT0 ¼128 Hz. (a) ω ¼20 Hz, where 6X¼120 Hz (6Xo ωT0) (b) ω ¼21.34 Hz, where 6X ¼128 Hz (¼ ωT0) and (c) ω ¼22 Hz, where 6X¼132 Hz (4 ωT0). For the three rotational speed cases considered, the peak amplitude in frequency spectra of torsional vibration response appears at 120 Hz, 128 Hz and 132 Hz as

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

x 10

(b)

−5

5

Contact duration

Instantaneous frequency (Hz)

Windowed torsional response

(a)

105

0

300 Non engaged duration

200 100

−5 7.8

7.85

(c)

7.9 Time (s)

Instantaneous frequency (Hz)

128Hz(ωT0)

128Hz(ωT0)

0 0

7.95

1,000

Contact duration

0.05

0.1 Time (s)

0.15

0.2

803 Hz (ωS0)

803 500 1 2

0 0

3

4 5

0.05

6

7

0.1

0.15

0.2

0.25

Time (s) Fig. 8. (a) Windowed signal of torsional vibration, (b) HHT of windowed signal and (c) HHT for stator lateral response.

shown in Fig. 9. It is noted that the frequency component of the largest amplitude is closest to the torsional natural frequency ωT0. The peak to peak amplitude of torsional oscillations in three cases are found to be 181.7 mil, 513.7 mil and 423.9 mil respectively for (a), (b) and (c). It is clear that when one of the harmonics of rotational speed matches with the torsional natural frequency (i.e. when the repetition rate of torsional excitation or its harmonic is closer to the natural frequency in torsion) the torsional vibration amplitude is maximum. A detailed spectrum cascade of torsional vibration is discussed later in Section 3.4. 3.3. Influence of system parameters on the rotor response In this section, influence of system parameters on the dynamics of rotor vibration is discussed. Apart from the rotor speed, variation in parameters such as friction coefficient m between the contacting bodies, clearance ratio cg, and the stator flexibility (indicated by Es) are studied for their influence on the dynamic response. 3.3.1. Friction coefficient The influence of friction coefficient on the overall response of the rotor system is analysed. Figs. 10a and b show the torsional vibration response of the rotor and lateral vibration response of the stator respectively with friction coefficient value varied as 0.1, 0.3 and 0.5. It is clear that the amplitude of vibration changes significantly with change in friction coefficient. The amplitude of torsional response and stator lateral response increases with increase in friction coefficient due to greater tangential force. The response character in terms of its frequency content, however, is unchanged. In case of higher speed ratio such as r ¼3, the torsional response amplitude increases with coefficient of friction value as shown in Fig. 11a. The orbit shape (Fig. 11b) has not changed significantly, although slight change in the orientation of the orbit with higher m is observed. Increase in tangential force due to higher friction coefficient allows the locus of rotor centre to have larger amplitude in horizontal direction. A considerable increase in the amplitudes of torsional and horizontal vibration is observed with increase in friction coefficient from 0.1 to 0.5.

120Hz

0.05

0 0

100 200 300 Frequency (Hz)

0.2

(c) 128Hz

0.1

0 0

100 200 300 Frequency (Hz)

Amplitude (mm)

(b)

0.1

Amplitude (mm)

Amplitude (mm)

(a)

0.15

132Hz

0.1 0.05 0 0

100 200 300 Frequency (Hz)

Fig. 9. FFT of torsional response for rotational speed (a) ω¼20 Hz, (b) ω¼ 21.34 Hz and (c) ω¼ 22 Hz.

106

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Fig. 10. Response waveform for r ¼0.75, Es ¼ 3e9, cg ¼ 0.85 (a) Torsional and (b) Stator lateral. 1

Frequency spectrum obtained from horizontal response is shown in Fig. 12a that depicts the strong presence of 2 X component. As explained earlier in Section 3.2 (Fig. 5), the pseudo critical speed in case of rub occurs at r¼1.7 for the present rotor 1

configuration. The rotational speed now is at r¼3 which is close to the twice the pseudo lateral critical speed. The 2 X frequency component generates resonance and exhibits a double looped orbit (Fig. 11b). The Hilbert spectrum is obtained for the torsional vibration response and the plot is shown in Fig. 12b. Although the energy content is low, the consistent appearance of torsional natural frequency 128 Hz (ωT0) in the spectrum indicates the torsional excitation of the rotor due to contact. 3.3.2. Gap or clearance between the rotor and the stator To investigate the influence of gap between rotor and stator on the rotor response, the gap is reduced to 84 mm (cg ¼ 0.7) from 100 mm (cg ¼0.85) used earlier (Fig. 6). The rotor vibration response is shown in Fig. 13. Unlike the case of larger clearance (100 mm), the vertical vibration response here shows periodicity of rotational cycle time (Fig. 13a). In case of larger gap (100 mm), the extent of interaction, in terms of the stator axial deflection or the reaction force, observed during contact is less as compared to the case of smaller gap (84 mm). The stator axial deflection and the reaction force increases roughly by 35% with 15% reduction in the clearance. The smaller clearance causes larger contact reaction force and correspondingly larger deflection in stator response. It may also be recalled that the periodicity of 2T (T ¼ period of 1 cycle rotation) in the vertical response resulted in the appearance of sub-harmonic component

( X ) in its frequency spectrum (Fig. 6b). In case of 1 2

smaller gap, the reaction forces due to contact in successive cycles are found to be same. This results in single loop orbit and harmonic vertical vibration response. In the torsional vibration response shown in Fig. 13b, it can be observed that the response amplitude has decreased by 10% compared to the case of larger clearance (Fig. 6e). With reduction in clearance, the torsional vibration waveform clearly shows the instants of engagement and disengagement. With reduced clearance this sudden change in waveform appears more prominent and the start and end of rotor stator contact interaction in time domain plot can be identified from the torsional acceleration plot (Fig. 13f). It is also noticed that the period of repetition of contact reaction force is more consistent in the case of reduced clearance (Fig. 13d) compared to the case of larger clearance (Fig. 7a). 3 The frequency spectrum shown in Fig. 14a for vertical response shows 2X and 3X frequency components. The 2 X and 1 X frequency components observed for higher clearance case (Figs. 6b and 6d) disappears as the response is periodic 2 with rotational frequency. In torsional vibration response, the frequency spectrum (Fig. 14b) shows harmonics of rotational speed and the largest amplitude is noticed for the 5th harmonic at 118.3 Hz, which is closest harmonic component to the torsional natural frequency (ω T0 ¼128 Hz). The time-frequency domain analysis using Hilbert Huang

Fig. 11. Response waveform for r ¼3, Es ¼ 3e9, cg ¼ 0.7 (76 mm) (a) Torsional and (b) Orbit.

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

2

(b)

−4

x 10

Amplitude (m)

47.25Hz (1X/2)

1

94.5Hz (1X)

3X/2

0 0

50

100

150

200

Instantaneous frequency (Hz)

(a)

107

350 300 250 200 150 100

0

128Hz

128Hz (ωT0)

50 0

0.02

0.04

Frequency (Hz)

128Hz

0.06

0.08

Time (s)

Fig. 12. (a) FFT for horizontal vibration response and (b) HHT plot for torsional vibration response; (r¼ 3, Es ¼3e9, m ¼ 0.3, cg ¼ 0.7).

7.92

7.94

7.96

7.98

Stator lateral response (mm)

0.5 0 7.96

7.98

8

Time (s)

0

−0.1 −0.3

7.9

(f)

0.1

−0.2

−0.1

0

0.1

Horizontal, z (mm)

0.2

0.3

start end

−0.05

(d)

1

7.94

0

8

Time (s)

7.92

0.05

Contact reaction force (N)

−0.1 7.9

7.9

Vertical, y (mm)

Torsion, θ (mil)

0

(c)

(e)

(b)

0.1

Torsional accl. (rad−s−2)

Vertical, y (mm)

(a)

7.92

7.94

7.96

7.98

8

7.98

8

7.98

8

Time (s) 0

−5

−10 7.9

7.92

7.94

7.96

Time (s) 50

end start

0

−50 7.9

7.92

7.94

7.96

Time (s)

Fig. 13. Time domain plots for r ¼0.75, Es ¼ 3e9, m¼ 0.3 and cg ¼ 0.7 (84 mm); (a) Vertical (b) Torsional (c) Stator lateral (d) Contact reaction force (e) Orbit and (f) Torsional acceleration.

transform is done for the torsional vibration data and the instantaneous frequency-time plots are shown in Fig. 14c. The procedure followed to obtain the HHT of torsional response at cg ¼0.85 is repeated for this case (cg ¼0.7). The HHT for first IMF indicates the consistent presence of torsional natural frequency (ω T0 ¼128 Hz) during its non-engaged condition. It is noted that the torsional vibration energy content is high and uniform during the non-engaged conditions of rotor and stator. The Hilbert spectrum for stator lateral response (Fig. 14d) indicates the presence of stator natural frequency (803 Hz) for the span ‘2–3’ during which more than 80% of the amplitude decays. The span ‘1–2’ is contact duration where the energy is spread over a wider band of frequency. In the span ‘3–4’, the stator comes closer to stationary condition and remains so until the next contact happens. Fig. 15 shows a comparative time domain plots for three different values of clearance: 110 mm, 90 mm and 70 mm with corresponding cg values equal to 0.85, 0.70 and 0.54 respectively. The vertical and torsional response amplitude decreases with increase in gap (Figs. 15a and c). The extent of interaction increases by reducing the gap therefore the reaction force and the stator response (lateral and axial direction) increases in amplitude. The horizontal response amplitude does not change significantly whereas it can be observed (Fig. 15b) that the mean has shifted in the direction of tangential force exerted on the rotor due to contact. The orbit plots (Fig. 15f) also indicates the variation of vertical and horizontal response with respect to change in gap. The sudden discontinuity in the waveform of torsional acceleration shown in Fig. 15d is noticed at the instant of impact as the extent of interaction has increased because the stator considered is stiffer (Es ¼ 30 GPa). The reaction force amplitude (Fig. 15e) decreases with increase in gap.

108

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Amplitude (m)

8

x 10

−5

1X

6 4 2X

2

3X

0 0

50 100 Frequency (Hz)

3

x 10

118Hz (5X) 142Hz (6X)

2 1 1X

0 0

150

(c)

100 200 Frequency (Hz)

300

(d) 1 2

300

Instantaneous frequency (Hz)

Instantaneous frequency (Hz)

−5

(b) Amplitude (m)

(a)

Contact duration

200

Non engaged duration

100 0 0

128Hz

0.05

0.1 Time (s)

128Hz

0.15

0.2

3

4

1,000

803Hz

803 500

0 0

0.05 Time (s)

0.1

Fig. 14. Frequency spectra for (a) Vertical response, (b) Torsional response and Hilbert spectra for (c) Torsional response, (d) Stator lateral response; (r¼ 0.75, Es ¼3e9, m ¼ 0.3, cg ¼ 0.7).

Fig. 15. Time domain response plots for r¼ 2, Es ¼ 30e9, m¼ 0.3 (a) Vertical, (b) Horizontal, (c) Torsional, (d) Torsional acceleration, (e) Reaction force and (f) Orbit. (Legends as in ‘e’ and ‘f’).

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

−3

4 3

126.3Hz (2X)

2 1X

1 0 0

−3

(b)

x 10

50

100 150 200 Frequency (Hz)

4 Amplitude (rad)

Amplitude (rad)

(a)

250

109

x 10

3 126.3Hz (2X)

2 1

1X

0 0

300

50

100 150 200 Frequency (Hz)

250

300

Instantaneous frequency (Hz)

(c) 128Hz (ωT0)

300 200 100 0 0

0.01

0.02

0.03

0.04

0.05 Time (s)

0.06

0.07

0.08

0.09

Fig. 16. Frequency domain plots for torsional response at r ¼2 (a, b) FFT for cg ¼ 0.55 and 0.85 and (c) HHT for cg ¼ 0.55.

The frequency spectra of torsional vibration response for two extreme cases of gap 70 mm and 110 mm is shown in Figs. 16a and b. The dominating frequency is 126.3 Hz which is harmonic of rotational speed (2X) and very close to the torsional natural frequency, ωT0. The Hilbert spectrum shown in Fig. 16c is for gap 70 mm where the spectrum for first IMF contains the torsional natural frequency, ωT0 (128 Hz) during non-contact conditions. 3.3.3. Stator stiffness In order to investigate the effect of stator stiffness, the elastic modulus of the stator material is varied. The values considered for the simulations are 3 GPa, 15 GPa and 30 GPa. The effect of increase in stator stiffness is depicted for two speed ratio's r ¼0.75 and r ¼1.5 in Fig. 17 and Fig. 18 respectively. The time domain response for the two extreme cases (Es ¼ 3 GPa and Es ¼30 GPa) for speed r ¼0.75 is shown in Fig. 17. The orbit plot (Fig. 17a) indicates that the amplitude variations in horizontal and vertical vibration response are not significant; however the distortion in waveform and the amount of penetration reduces significantly with increase in stator stiffness. The contact duration as observed from reaction force plot (Fig. 17d) reduces from 3ms to 1.2ms. Similarly, change in orbital response with stator stiffness is observed for higher speed ratio, r ¼1.5 (Fig. 18a). The orbit plot shows elongated large size orbit for flexible rotor and as the stiffness increases, the orbit size shrinks with reduction in vertical and horizontal vibration response. With increase in stator stiffness, the amplitude of stator axial deflection in both the cases (Fig. 17c and Fig. 18c) drops significantly. However the drop for speed r¼1.5 is larger (Fig. 18c) as the operating speed is near the pseudo critical speed of rotor. The contact reaction force is proportional to the stator axial deflection and its stiffness. For speed r¼0.75, the stator stiffness dominates and increases the contact reaction force (Fig. 17d), however for r¼1.5 the drop in axial deflection of stator influences the contact reaction force (Fig. 18d) that is found to change marginally. Corresponding to the generated contact reaction force, the tangential friction force changes and it governs the torsional vibration amplitudes in respective cases (Fig. 17b and Fig. 18b). The instants of engagement and disengagement distinctly appear in the waveform of torsional acceleration data. The stator lateral response shown in Fig. 19 for both the speed cases (r ¼0.75 and r ¼1.5) indicates smaller response amplitude and fast decaying curve for stiffer stator. The amplitude decays rapidly because of higher damping (proportional to mass and stiffness) present in the system. The foregoing parametric study shows following influence of various parameters on response features: Effect of friction coefficient: Amplitude of torsional vibration response of rotor (Fig. 10a) and lateral vibration of stator (Fig. 10b) increases significantly with increase in friction coefficient due to higher friction torque. Effect of clearance: With reduced clearance between rotor and stator, the contact reaction force increases (Fig. 13d) and hence more prominent indication of contact instances appear in the torsional vibration response data (Fig. 13b). In terms of 1

frequency content, the larger clearance case (cg ¼ 0.85) generates sub-harmonic response of order 2 ( 2 X component in frequency plot (Fig. 6b)) which is not observed in the frequency spectra with lower clearance cg ¼0.7 (Fig. 14a). At a higher operating speed (r ¼2), increase in clearance decreases the amplitudes for vertical and torsional vibration responses and the contact reaction force (Fig. 15).

110

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

Fig. 17. Time response plots for r ¼0.75, m ¼ 0.3, cg ¼0.7 (g ¼ 84 mm) (a) Orbit, (b) Torsional, (c) Stator lateral vibration and (d) Contact reaction force.

Fig. 18. Time response plots for r ¼1.5, m ¼ 0.3, cg ¼ 0.7 (g ¼ 150 mm) (a) Orbit, (b) Torsional, (c) Stator axial response and (d) Contact reaction force.

Effect of stator stiffness: Increasing the stator stiffness results in significant enhancement of contact reaction force (Fig. 17d) thus a considerable increase in the amplitude of torsional response is observed (Fig. 17b). The dynamic response analysis has shown the reasons and conditions for existence of sub-harmonic frequency components in the rotor vibration response. These frequencies have been proposed as general rub indicators in the past literature; however, these are not necessarily always observed in the response, unless the operating speed is favorable for their presence.

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

111

Fig. 19. Stator lateral response for (a) r ¼0.75 and (b) r ¼1.5.

3.4. Torsional vibration spectral features For the overall assessment of torsional vibration features, the spectrum cascade shown in Fig. 20 is obtained for the steady state torsional vibration response over the speed range from r ¼ 0.9 to r ¼ 6.1. The other simulation parameters include m ¼ 0.3, cg ¼ 0.7 and Es ¼ 30 GPa. One of the prominent features of the spectrum is the presence of large amplitude resonance near the torsional natural frequency (ωT0 ¼ 128 Hz). Depending on the operating speed, larger amplitudes of 1

1 X) 3

torsional vibrations are noticed when one of the sub-harmonic ( 2 X ,

or super-harmonic (2X, 3X) of the rotational

frequency component gets closer to ωT0. For example, when rotating at r ¼ 1.9, the 2X frequency component is closer to ωT0 and generates large amplitude (marker #d), while at r ¼ 1.3, the 3X frequency component generates resonant amplitude (marker #e). Similar behaviour is observed at speed r ¼ 5.7 where 14th harmonic of 14 (5.7  20

1 X 20

frequency component

 31.56 ¼ 126 HzE ωT0) attains largest amplitude (marker #f). At r ¼ 4.1, 1X frequency component exhibits re1

sonant amplitude in the spectrum (marker #g). In essence, when m × n X frequency component of rotational frequency comes closer to the torsional natural frequency ωT0, stronger amplitude is noticed in the spectrum cascade. The other important observation in the torsional vibration spectrum is the appearance of larger vibration amplitudes at pseudo bending natural frequency (53.5 Hz), a strong indication of coupling of lateral and torsional vibration modes. As explained in the Section 3.2 (ref. Fig. 5), the pseudo bending vibration resonance phenomenon is observed at r ¼1.7. Thus, if the operating speed is close to the pseudo bending critical speed, prominent torsional vibration amplitude is observed in the spectrum cascade. Even when the 1/nth harmonic of rotational frequency comes closer to the pseudo bending critical speed, the resonance condition in the lateral vibration prevails and it couples with the torsional mode to generate larger

1 X n

frequency component in the torsional vibration. These sub-harmonic frequency components are indicated by marker #a for n ¼ 1, marker #b for n ¼ 2 and marker #c for n ¼ 3 in Fig. 20. The peaking up of amplitudes across these markers indicates the passage through sub-harmonic resonances. There are other spectral features in the cascade such as larger amplitudes of 1

3rd and 4th harmonics of 3 X for speed near r ¼ 5.3 (marker #h and #i). These are in general due to large overall vibration f h

X/2

X/3

c

1X

i (4X/3)

k

3.5

d 3X

a

Pseudo critical speed (53.5 Hz)

25

50

2.9 2.3 1.7

e

0

4.1

2X

2X/3

75

1.1

Torsional natural frequency, 128Hz (ωT0)

100

125 150 Frequency (Hz)

175

200

225

250

Fig. 20. Spectrum cascade for torsional vibration response over speed range from r¼ 0.9 to r ¼6.1 in step of Δr ¼0.2.

Rotational speed ratio, r

j b

5.9 5.3 3X/2 4.7

g

112

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

1

amplitudes caused by sub-harmonic resonance with pseudo bending natural frequency. Similarly, higher harmonics of 2 X frequency are seen prominently in the spectrum near operating speed of r ¼ 3.5 (marker #j and #k). In this case, the sub1

harmonic resonance (marker #b) generates larger amplitudes for multiples of 2 X frequency components. One may also note the resonance build up across the ωT0, when one follows the 3 X , 1X, 2X, 3X excitation lines in the spectrum cascade. This is 2

similar to that observed across the passage through the pseudo-bending critical speed. In practice, when the sensors are mounted on both stator and rotor supports, the vibration signals from both the sensors in appropriate mode (stator structural vibration, rotor lateral vibration and rotor torsional vibration) can be used in conjunction to infer about the possible rub phenomenon in a more comprehensive manner based on these vibration features.

4. Conclusions The present work investigates the rotor stator rub phenomenon in an FE framework using Lagrange multiplier based contact mechanics approach. The stator is modelled more realistically and the complex rotor-stator rub phenomenon is analysed considering the stator dynamics including its axial and lateral degrees of freedom. The contact strategy based on instantaneous and independent motion of rotor and stator ensures that the system dynamics is adequately and realistically captured. The lateral–torsional coupling of rotor during contact is studied and vibration response is analysed. In addition to the rotor lateral vibration response and stator dynamic response, the present study focusses on the investigations of torsional vibration features of the rotor which is the first attempt in light of the contact mechanics based approach. The effect of parameters like rotor-stator clearance, stator stiffness and friction coefficient are also studied. The orbits around pseudo bending critical speed show period-1 motion and exhibits predominantly second and third harmonics in the frequency spectra. However sub harmonics may appear in the frequency domain at subcritical speeds with 1 lower contact intensity such as the case of larger gap (cg ¼ 0.85). Further at supercritical speeds, strong sub-harmonics of 2 X 1 and 3 X can be observed near the second and third multiple of pseudo bending critical speeds of the rotor respectively. A strong influence of clearance on the dynamic response is noticed; lower level of clearance suppresses the presence of subharmonic vibration frequencies. A strong coupling between lateral and torsional natural frequency is clearly observed in the spectrum cascade of the torsional vibration due to friction contact. The distinct presence of harmonics of rotational frequency (in absence of any explicit torsional excitation) in frequency spectra of torsional vibration due to the strong coupling indirectly indicates the presence of rub phenomenon. Further, Hilbert Huang spectrum indicates the intermittent appearance of torsional natural frequency, particularly during non-contact intervals, when the rotor is free to oscillate at its torsional natural frequency. This intermittent rub is effectively captured using HHS. Interesting resonance phenomenon related to both torsional and pseudo bending natural frequency are noticed due to the presence of submultiples and integer multiple of rotational speed components. The stator is excited periodically due to rotor interaction through a series of impulses and hence generates vibratory response comprising mainly the structural resonance frequencies. It is proposed that the measured data from stator structural vibration and rotor lateral-torsional vibration can be correlated and used for rub identification.

Acknowledgment This research is supported by Aeronautics Research and Development Board, Defence Research and Development Organization, Government of India (Project ref. no. DARO/08/1041661/M/I). The authors gratefully acknowledge the financial support.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

D.W. Childs, Fractional-frequency rotor motion due to nonsymmetric clearance effects, J. Eng. Power, Trans. ASME 104 (1982) 533–541. D.W. Childs, Rub-induced parametric excitation in rotors, J. Mech. Des., Trans. ASME, 101, , 1979, 640–644. R.F. Beatty, Differentiating rotor response due to radial rubbing, J. Vib. Acoust. Stress Reliab. Des. 107 (1985) 151–160. F.K. Choy, J. Padovan, Non-linear transient analysis of rotor-casing rub events, J. Sound Vib. 113 (1987) 529–545, http://dx.doi.org/10.1016/S0022-460X(87)80135-9. A. Muszynska, Rotordynamics, CRC Taylor and Francis Group, NewYork, 2005. P. Goldman, A. Muszynska, Rotor-to-stator, Rub-related, Thermal/Mechanical Effects in Rotating Machinery, Chaos, Solitons Fractals 5 (1995) 1579–1601. G.V. Groll, D.J. Ewins, A mechanism of low subharmonic response in rotor/stator contact-measurements and simulations, J. Vib. Acoust. 124 (2002) 350–358, http://dx.doi.org/10.1115/1.1467648. F. Chu, Z. Zhang, Periodic, quasi-periodic and chaotic vibrations of a rub impact rotor system supported on oil film bearings, Int. J. Eng. Sci. 35 (1997) 963–973. Z. Sun, J. Xu, T. Zhou, Analysis on complicated characteristics of a high-speed rotor system with rub-impact, Mech. Mach. Theory 37 (2002) 659–672. F. Chu, W. Lu, Experimental observation of nonlinear vibrations in a rub-impact rotor system, J. Sound Vib. 283 (2005) 621–643.

M.A. Mokhtar et al. / Journal of Sound and Vibration 401 (2017) 94–113

113

[11] F. Cong, J. Chen, G. Dong, K. Huang, Experimental validation of impact energy model for the rub-impact assessment in a rotor system, Mech. Syst. Signal Process. 25 (2011) 2549–2558. [12] F. Chu, W. Lu, Stiffening effect of the rotor during the rotor-to-stator rub in a rotating machine, J. Sound Vib. 308 (2007) 758–766. [13] T.H. Patel, A.K. Darpe, Study of coast-up vibration response for rub detection, Mech. Mach. Theory 44 (2009) 1570–1579. [14] F. Chu, W. Lu, Determination of the rubbing location in a multi-disk rotor system by means of dynamic stiffness identification, J. Sound Vib. 248 (2001) 235–246. [15] S. Edwards, A.W. Lees, M.I. Friswell, The Influence of Torsion on Rotor/Stator Contact in Rotating Machinery, J. Sound Vib. 225 (1999) 767–778. [16] Z. Sun, J. Xu, T. Zhou, N. Tan, Study on the influence of bending-torsional coupling in an impacting-rub rotor system, Appl. Math. Mech. 24 (2003) 1316–1323. [17] B.O. Al-Bedoor, Modeling the coupled torsional and lateral vibrations of unbalanced rotors, Comput. Methods Appl. Mech. Eng. 190 (2001) 5999–6008. [18] B.O. AL-Bedoor, Transient torsional and lateral vibrations of unbalanced rotors with rotor-to-stator rubbing, J. Sound Vib. 229 (2000) 627–645. [19] M.A. Mohiuddin, Y.A. Khulief, Coupled bending torsional vibration of rotors using finite element, J. Sound Vib. 223 (1999) 297–316. [20] T.H. Patel, A.K. Darpe, Coupled bending-torsional vibration analysis of rotor with rub and crack, J. Sound Vib. 326 (2009) 740–752. [21] H.M. Khanlo, M. Ghayour, S. Ziaei-Rad, The effects of lateral-torsional coupling on the nonlinear dynamic behavior of a rotating continuous flexible shaft-disk system with rub-impact, Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 1524–1538. [22] S. Roques, M. Legrand, P. Cartraud, C. Stoisser, C. Pierre, Modeling of a rotor speed transient response with radial rubbing, J. Sound Vib. 329 (2010) 527–546. [23] M. Behzad, M. Alvandi, D. Mba, J. Jamali, A finite element-based algorithm for rubbing induced vibration prediction in rotors, J. Sound Vib. 332 (2013) 5523–5542. [24] H. Ma, C. Shi, Q. Han, B. Wen, Fixed-point rubbing fault characteristic analysis of a rotor system based on contact theory, Mech. Syst. Signal Process. 38 (2013) 137–153. [25] H. Ma, Z. Wu, X. Tai, B. Wen, Dynamic characteristics analysis of a rotor system with two types of limiters, Int. J. Mech. Sci. 88 (2014) 192–201. [26] H. Ma, Q. Zhao, X. Zhao, Q. Han, B. Wen, Dynamic characteristics analysis of a rotor-stator system under different rubbing forms, Appl. Math. Model 39 (2015) 2392–2408. [27] M.A. Mokhtar, K. Gupta, A.K. Darpe, A comparative study of various constraint enforcement techniques for rotor stator rub, in: Proceedings of the IMechE’s 12th International conference on Vibrations in Rotating Machinery- VIRM 11, September 13–15, 2016. [28] G. Jacquet-Richardet, M. Torkhani, P. Cartraud, F. Thouverez, T. Nouri Baranger, M. Herran, C. Gibert, S. Baguet, P. Almeida, L. Peletan, Rotor to stator contacts in turbomachines. Review and application, Mech. Syst. Signal Process. 40 (2013) 401–420. [29] P. Wriggers, Computational Contact Mechanics, Springer-Verlag, Berlin Heidelberg, 2006. [30] M. Legrand, S. Junca, S. Heng, Nonsmooth modal analysis of a N-degree-of-freedom system undergoing a purely elastic impact law, hal-01185980, 2015, pp. 1–31. [31] D. Doyen, A. Ern, S. Piperno, Time-Integration schemes for the finite element dynamic signorini problem, SIAM J. Sci. Comput. 33 (2011) 223–249. [32] N.J. Carpenter, R.L. Taylor, M.G. Katona, Lagrange constraints for transient finite element surface contact, Int. J. Numer. Methods Eng. 32 (1991) 103–128. [33] F.F. Ehrich, High order subharmonic response of high speed rotors in bearing clearance, J. Vib. Acoust. Stress Reliab. Des. 110 (1988) 9–16. [34] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. Lond. A 454 (1998) 903–995.