Analytical modeling of asymmetric multi-segment rotor – bearing systems with Timoshenko beam model including gyroscopic moments

Computers and Structures 144 (2014) 119–126

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Analytical modeling of asymmetric multi-segment rotor – bearing systems with Timoshenko beam model including gyroscopic moments O. Özsßahin a,⇑, H.N. Özgüven a, E. Budak b a b

Department of Mechanical Engineering, Middle East Technical University, 06531 Ankara, Turkey _ 81474, Turkey Manufacturing Research Laboratory, Sabanci University, Istanbul

a r t i c l e

i n f o

Article history: Received 13 May 2014 Accepted 11 August 2014 Available online 16 September 2014 Keywords: Rotor dynamics High speed rotors Gyroscopic moments Timoshenko beam Multi-segment rotor Asymmetric bearing

a b s t r a c t In this study, analytical modeling and an analysis approach for asymmetric multi-segment rotor – bearing systems are presented. Timoshenko beam model which includes the effect of gyroscopic moments is employed for modeling rotor segments. Instead of applying FEM, sub-segment Frequency Response Functions (FRFs) are obtained analytically, and sub-segment FRFs obtained are coupled by using receptance coupling method. Bearing properties are included into system dynamics by employing structural modification techniques. The proposed analytical model is verified by using FEM approach. It is shown that using analytical model and receptance coupling, compared with FEM, reduces computational time drastically without losing accuracy. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Accurate determination of the dynamic response of rotating machinery plays a crucial role in the design state and failure prevention of rotating machinery. In the rotor dynamics research area, in order to determine the dynamic response of the system, Nelson [1] applied finite element method (FEM) and constructed the element matrices for rotating Timoshenko beam, and FEM became the main analysis tool in the rotor dynamics research area in the last two decades [1–8]. In addition to the FEM, Frew and Scheffer [9] applied numerical methods based on modified Euler equations, and Chena et al. [10] applied pseudo mode shape method for rotor – bearing foundation identification and integrated the identified foundation dynamics to the FEM. However; in the literature there exist limited numbers of studies that concentrate on the analytical modeling of rotor systems. For the continuous beam model of rotor systems, Lee et al. [11–13] proposed modal analysis solution based on non-self-adjoint system characteristics. Similarly, Wang and Kirkhope [14,15] proposed eigensolution and modal analysis for the undamped rotor systems and applied perturbation analysis for the damped cases. Parker and Sathe [16] proposed an exact solution for free and forced vibration of rotating disk spindle systems. In these analytical studies, rotors are considered uniform and modeled using Rayleigh beam model. However, for low slen⇑ Corresponding author. Tel.: +90 542 454 80 00. E-mail address: [email protected] (O. Özsßahin). http://dx.doi.org/10.1016/j.compstruc.2014.08.001 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

derness ratios, shear deformation becomes important at high frequencies and Rayleigh beam model which includes rotary inertia effects but neglects shear deformation does not provide accurate results. Therefore, for accurate modeling of the system, Timoshenko beam model should be used. Another pitfall of the rotor dynamics is the multi segment characteristics of the rotor geometry. For instance; in the classical eigensolution of an m-segment beam, EVP problem requires the solution of a characteristic equation expressed in terms of a 4m  4m matrix, where m is the number of segments of the rotor system which increases the size of the matrices. To overcome this complexity, Schmitz and Donaldson [17] and Ertürk et al. [18] proposed receptance coupling procedure to determine end point FRFs of multi segment beams. However; their study deals with non rotating shafts and neglects gyroscopic effects. For modeling of multi segment rotors Hong and Park [19] offers a method based on distributed transfer function synthesis (DTFS) which requires the use of Laplace transformations and construction of global matrix in the same manner as in FEM. In this present study, an analytical modeling approach for asymmetric multi-segment rotor systems is presented. Timoshenko beam model which includes the effect of gyroscopic moments is employed for the modeling of rotor segments. Instead of applying FEM, eigensolution of rotating Timoshenko beam is performed based on the method proposed by Aristizabal-Ochoa [20] which is derived for non-rotating Timoshenko beams, and thus the sub-segment Frequency Response Functions (FRFs) are

120

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126

obtained analytically using non-self-adjoint system characteristics. To couple the sub-segments dynamically and include the bearing dynamics to the system, FRF coupling is used along with structural modification technique as successfully used in the study of Ertürk et al. [18]. Being different from the approach proposed in [18], cross FRFs between two orthogonal planes are also included in the receptance matrices of rotor sub-segments. The whole rotor is modeled using symmetric beam elements with free–free end conditions and asymmetric bearing properties are added to the rotor employing structural modification techniques. Thus, response of an asymmetric rotor – bearing system is obtained using symmetric elements. The proposed analytical modeling approach is verified by using FEM. A spindle–holder –tool system is used in order to illustrate the effect of including gyroscopic moment in the analysis. The main advantage of using analytical modeling compared with FEM is the considerable saving in computational time, which is illustrated with a numerical example. 2. Mathematical modeling

Consider a beam element with generalized coordinates as shown in Fig. 1. Here uz is the axial displacement, ux and uy are the lateral displacements and /x , /y , and /z are the rotations with respect to x, y and z axes, respectively. Based on the generalized coordinates given in Fig. 1, the equations of motion for the rotating Timoshenko beam can be written as

  @/ @ / @ / @u qI 2x þ 2XIq y  EI 2x þ kAG /x þ y ¼ 0 @t @z @z @t   2 2 @ / @ / @/ @u qI 2y  2XIq x  EI 2y þ kAG /y  x ¼ 0 @t @z @z @t 2

ð1Þ ð2Þ

Here, q is the density, E is Young’s modulus, G is the shear modulus, A is the cross sectional area, I is the cross sectional area moment of inertia, k is the shear coefficient, and X is the rotational speed of the beam. The boundary conditions of the rotating Timoshenko beam in y direction for free–free end conditions can be expressed as follows:

ux ðz; tÞ ¼ U x ðzÞeixt ;

/x ðz; tÞ ¼ hx ðzÞeixt

ð8Þ

Note that due to the symmetricity, for the free vibration analysis of the beam element, the mode shapes of the beam in two orthogonal planes will be related with each other’s by the following relations [21]: f

b

U fx ðzÞ ¼ iU y ðzÞ U bx ðzÞ ¼ iU y ðzÞ

ð9Þ

where superscripts f and b stand for forward and backward modes respectively. If the expressions given by the Eqs. (8) and (9) are substituted into Eq. (7), differential equations representing motions in two orthogonal planes can be decoupled and the equation of motion can be written as ordinary differential equations for the backward motion as



  2 b d Uy EqI 2 x  2 q I X x 4 2 kG dz dz  2  2 qI 4 qI þ x  2 Xx3  qAx2 U by ¼ 0 kG kG 4

ð3Þ ð4Þ

ð7Þ

Note that, the same equation of motion can be written for the x direction as well. As seen from the Eqs. (3), (4) and (7), due to the gyroscopic effects, motion in two orthogonal planes are coupled. Thus classical solution methods cannot be applied for the solution of rotating Timoshenko beam. For the solution of rotating Timoshenko beam, harmonic vibration can be assumed and separation of variables can be applied such that

EI

2

    @/y @/y M y ð0; tÞ ¼ EI ¼ 0; My ðL; tÞ ¼ EI ¼0 @z z¼0 @z z¼L       @ux @ux Sy ð0; tÞ ¼ kAG  /y ¼ 0; Sy ðL; tÞ ¼ kAG  /y ¼0 @z @z z¼0 z¼L

  4 @ 4 uy @ 2 uy E @ uy q2 I @ 4 uy þ q A  q I 1 þ þ y 2 2 4 2 kG @z @t kG @t 4 @z @t !  2  3 @ @ux q @ ux  þ 2qIy X ¼0 @z2 @t kG @t 3

EIx

uy ðz; tÞ ¼ U y ðzÞeixt ;

2.1. Modeling of rotating Timoshenko beam

! @/y @ 2 ux @ 2 ux qA 2  kAG  þ 2 ¼0 @z @z @t ! @ 2 uy @/x @ 2 uy qA 2  kAG þ 2 ¼0 @z @z @t

where L is the length of the beam. A similar set of equations can be written in x direction. By eliminating the bending rotations /x and /y , equation of motion of the rotating Timoshenko beam in y direction will take the following form

d U by

qI þ

þ

ð10Þ

and, for the forward motion as



  2 f d Uy EqI x2 þ 2qIXx 2 kG dz dz  2  qI 4 q2 I þ x þ 2 Xx3  qAx2 U fy ¼ 0 kG kG 4

EI

d U fy 4

þ

qI þ

ð11Þ

2.2. Free vibration solution for rotating Timoshenko beam ð5Þ ð6Þ

The eigensolution for rotating Timoshenko beam can be obtained by using the solution procedure given by AristizabalOchoa [20] for non-rotating Timoshenko beam. For free–free end conditions, characteristic equation of rotating Timoshenko beam for backward motion can be written as

  D11  D 21

 D12  ¼ D11 D22  D12 D21 ¼ 0 D 

ð12Þ

22

where

D11 ¼ ða  kÞ cos ðaLÞ þ ðk  aÞ cosh ðbLÞ ka D12 ¼ ðk  aÞ sin ðaLÞ þ ðb  dÞ sinh ðbLÞ db ka sinh ðbLÞ  ka sin ðaLÞ D21 ¼ db bd D22 ¼ kaðcosh ðbLÞ  cos ðaLÞÞ Fig. 1. Rotating Timoshenko beam generalized coordinates.

a2  K k¼ a

b2 þ K d¼ b

K¼

qAx2 kAG

ð13Þ ð14Þ ð15Þ ð16Þ ð17Þ

121

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126

pffiffiffiffiffiffiffiffiffiffiffi

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi

a ¼ g þ e b ¼ g þ e pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 b b  4d g¼ e¼ h2 2 qI I þ EkG

q

b¼

ð18Þ ð19Þ

x2  2qIXx

i

h d¼

EI

q2 I kG

q2 I

x4  2 kG Xx3  qAx2

ð20Þ

EI

U by ðzÞ ¼ Ar ðC 1 sinðar zÞ þ C 2 cosðar zÞ þ C 3 sinhðbr zÞ þ C 4 coshðbr zÞÞ ð21Þ hbx ðzÞ ¼ Ar ½kr ðC 1 cosðar zÞ  C 2 sinðar zÞÞ þ dr ðC 3 coshðbr zÞ þ C 4 sinhðbr zÞÞ

ð22Þ

where

C3 ¼ C1

ar  kr dr  br

C 4 ¼ C 1

;

ar kr D11 dr br D12 ð23Þ

For forward motion, the free vibration solution can be obtained in a similar way. The only difference will be in the speed dependent terms expressed by Eq. (20). For forward motion b and d will be expressed as:

h

d¼

0 M



M M _ ¼ fwg G 0



0 fwg þ fQ ðtÞg K

ð29Þ

where

fwg ¼

 q_

fQ ðtÞg ¼

q



0

ð30Þ

fFðtÞg

Eq. (29) can also be written as:

_ ¼ ½Dfwg þ fQ ðtÞg ½Cfwg

ð31Þ

Note that in Eq. (31), ½D is self adjoint and ½C is non-self adjoint, therefore system will have left eigenvectors and their adjoint right eigenvectors. Thus, Eigen Value Problem (EVP) of rotating Timoshenko beam for left and right eigenvectors can be written as follows:

n o n o klr ½C Ulr ðzÞ ¼ ½D Ulr ðzÞ ; n o  Wl ðzÞ ¼ ½Df  Wl ðzÞg; kl ½C s s s

l ¼ b; f

ð32Þ

l ¼ b; f

ð33Þ

where

D11 C 2 ¼ C 1 ; D12

C1 ¼ C1;



qI þ EkGqI x2 þ 2qIXx

h



i

As can be seen from Eq. (20), rotating Timoshenko beam model has additional speed dependent terms, compared to the non-rotating case. Thus by using the characteristic equation given by Eq. (12), natural frequencies of each elastic mode and corresponding frequency numbers ar and br can be determined. Finally, the eigenfunction expressions for the dynamic transverse deflection and bending rotation can be obtained as follows:

b¼



i

EI q2 I

2

x4 þ 2 qkGI Xx3  qAx2 kG

i ð24Þ

EI

fWs ðzÞgT ¼

^ ly ¼ klr uly u r r

qA 0 6 0 qI 6

0

0

07 7 7; 05

½M ¼ 6 4 0

0

qA

0

0

0

2

@2 kAG @z 2

3

0

2

0

0

60 6 ½G ¼ 6 40

0

0 2qIX 7 7 7 0 0 5

@ kAG @z

0

0 2qIX 0 0

0

0

n

¼ K r /lxr

U lyr

U lxr

hlxr

o hlyr ;

l ¼ b; f

a ^l U yr

a ^l h xr

a ^l U xr

a ^l hyr

a

U lyr

a l hxr

a

U lyr

a l hxr

o

; l ¼ b; f

a l uyr

¼ K r ulyr

ð36Þ

Eigenvectors can be simplified using the relation given below:

^ l ¼ kl /l / r xr xr

ð37Þ

For the non self adjoint EVP, left and right eigenvectors are biorthonormal and biorthonormality can be defined using the inner product operator as

hC Ur ; Ws i ¼ drs

3

7 6 @2 7 6 kAG @ kAG  EI @z 0 0 2 7 6 @z ½K ¼ 6 7 2 @ @ 7 6 0 0 kAG kAG 5 4 @z @z2 @ @2 0 0 kAG @z kAG  EI @z2 8 9 8 9 uy > F y ðtÞdðz  zn Þ > > > > > > > > > > > > > u > 0 > > > > x> > > > : : ; ; /y 0

fag ¼ f a1

3

0

qI

^hl yr

hDUr ; Ws i ¼ kr drs

ð38Þ

In Eq. (38) h; i is the inner product operator and defined as follows:

where M; G, and K are defined as follows:

2

^l U xr

In Eq. (33) bar denotes complex conjugate, and in Eqs. (34) and (35) superscript a denotes adjoint and superscript b andf denotes the backward and forward modes respectively. For non-self adjoint systems, left eigenvectors and complex conjugate of their adjoint right eigenvectors are related with a constant K r which is a pure imaginary number [2,11]

2.3. Forced vibration response of rotating Timoshenko beam

ð25Þ

^hl xr

ð35Þ

a l /xr

_ þ ½Kfqg ¼ fFðtÞg ½Mfqg þ ½Gfqg

^l U yr

ð34Þ

Also note that the constant Ar can be obtained by the normalization of the eigenfunctions using biorthonormality condition which will be discussed in Section 2.3.

To obtain the forced response of a rotating Timoshenko beam to an external force applied on the beam at location z ¼ zn in y direction, Eqs. (1)–(4) can be expressed in the following form

n

T

fUlr ðzÞg ¼

ð26Þ

T

a2 gT ; fbg ¼ f b1

ha; bi ¼ ha1 ; b1 i þ ha2 ; b2 i ¼

Z

L

b1 a1 dz þ

0

0

ð27Þ

Note that, due to the gyroscopic effects, the system is non-self adjoint. Therefore, in order to obtain forced response of a rotating Timoshenko beam, equation of motion should be represented in state space as

Z

L

b2 a2 dz

ð40Þ

0

To determine the forced response of the system, backward and forward modes can be treated separately. Thus if we apply the biorthonormality, the constant Ar in the eigenfunction expressions of backward and forward modes can be determined. For backward modes the constant Ar can be determined as follows:

A2r ¼ ð28Þ

ð39Þ

b2 g

1 K r kr C r  K r Er

ð41Þ

where Cr ¼

Z L  0

Er ¼

 2    2    2    2  2qA U byr þ 2qI hbxr þ 2qA U bxr þ 2qI hbyr dz ð42Þ

Z L  2  2  dz i2qIX hbxr þ i2qIX hbyr

ð43Þ

0

Similarly, Ar in the eigenfunction expressions of forward modes can be determined.

122

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126

Response of the system to an externally applied excitation can be expressed by using expansion theorem as follows:

XX

Ulr ðzÞqlr ðtÞ

wðz; tÞ ¼

ð44Þ

r

l¼b;f

If the expression given by Eq. (44) is substituted into Eq. (31) and biorthonormality is applied, Eq. (31) reduces to:

q_ lr ðtÞ

l

¼ q ðtÞkr þ

Plr ðtÞ

ð45Þ

where,

D n oE Z Pır ðtÞ ¼ fQ ðtÞg; Wls ¼

L



0

 l fQðtÞgdz W s

ð46Þ

Thus, solution in terms of modal coordinates can be obtained as:

qlr ðtÞ ¼

htrans x

Plr ðtÞ ix 

ð47Þ

klr

uy ðz; tÞ ¼

1 U l ðzÞa U  l ðz ¼ zn ÞF y ðtÞ XX yr yr

ð48Þ

ix  klr

l¼f ;b 1

Since in the state space, the order of the EVP is doubled, eigenvalues and corresponding eigenvectors will be in the form of complex conjugates. Thus, the transverse response of the system to the applied force F y ðtÞ can be written as

uy ðz; tÞ ¼

1 XX l¼b;f r¼0

! U lyr ðzÞU lyr ðz ¼ zn ÞF y ðtÞ 2ixr  2 x2 þ xlr

ð49Þ

For the rotating Timoshenko beam, the receptance functions in y–z plane can be defined as follows:

 j uy F iy

;

Niy;jx ¼

ð/x Þj F iy

;

Lix;jy ¼

 j uy M ix

;

Pix;jx ¼

ð/x Þj M ix

¼ 0;

! U lyr ðzj ÞU lyr ðzi Þ Hiy;jy ðxÞ ¼ 2ix  2 x2 þ xlr ð1 þ icÞ l¼b;f r¼0 ! 1 XX hlxr ðzj ÞU lyr ðzi Þ Niy;jx ðxÞ ¼ 2ixlr  2 x2 þ xlr ð1 þ icÞ l¼b;f r¼0 ! 1 XX U lyr ðzj Þhlxr ðzi Þ l Lix;jy ðxÞ ¼ 2ixr  2 x2 þ xlr ð1 þ icÞ l¼b;f r¼0 ! 1 XX hlxr ðzj Þhlxr ðzi Þ 2ixlr Pix;jx ðxÞ ¼  2 x2 þ xlr ð1 þ icÞ l¼b;f r¼0 l r

¼A

x 2

L



ð55Þ

rot

ð56Þ

Atrans ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ; ixr C trans  Etrans

Arot ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ixr C rot  Erot

ð57Þ

where

C trans ¼ 4qAL; Etrans ¼ 0  Z L   2  2 dz; 2qA U rot þ 2qA U rot C rot ¼ y x rot

E

¼

Z

0

L

  2   2 dz i2qIX hrot þ i2 q I X hrot x y

ð58Þ

ð59Þ

In addition to the point and cross FRFs in each orthogonal plane, due to the gyroscopic effects, there exist cross coupling between the orthogonal planes. Therefore; for the prediction of the complete dynamics of the rotating Timoshenko beam, cross FRFs between the orthogonal planes should be considered. These FRFs can be determined using the same approach presented in this paper but unlike the FRFs given by Eqs. (51)–(54), the FRF expressions for cross coupling terms eigenfunctions of the two orthogonal planes will appear in equations. However; using the relations given by Eq. (9), cross FRFs between the two orthogonal planes can be expressed by the eigenfunctions of each plane as follows: 0 1 ! 1 X U byr ðzi ÞU byr ðzj Þ U fyr ðzi ÞU fyr ðzj Þ B C þ Hix;jy ðxÞ ¼ 2x @2x Að60Þ  2  2 x2 þ xbr ð1 þ icÞ 0 0 x2 þ xfr ð1 þ icÞ 0 1 ! 1 1 X X hbxr ðzj ÞU byr ðzi Þ hfxr ðzj ÞU fyr ðzi Þ B C þ Niy;jy ðxÞ ¼ 2x @2x Að61Þ  2  2 f x2 þ xbr ð1 þ icÞ 2 0 0 x þ xr ð1 þ icÞ 0 1 ! f f 1 1 X X U byr ðzj Þhbxr ðzi Þ U ðz Þh ðz Þ j i x y B C r r þ Lix;jx ðxÞ ¼ 2x @2x A ð62Þ  2  2 x2 þ xbr ð1 þ icÞ 0 0 x2 þ xfr ð1 þ icÞ 0 1 ! 1 1 X X hbxr ðzj Þhbxr ðzi Þ hfxr ðzj Þhfxr ðzi Þ B C þ 2x Pix;jy ðxÞ ¼ @2x A ð63Þ  2  2 x2 þ xbr ð1 þ icÞ 0 0 x2 þ xfr ð1 þ icÞ 1 X

ð50Þ 2.4. FRFs of multi segment rotors by receptance coupling

Note that structural damping can be included into the model by adding modal damping terms to the denominators in Eq. (49) [21], and then point and cross receptance functions in y–z plane can be obtained as 1 XX

hrot x

rot U rot y ¼ A

Here, Atrans and Arot are constants and by applying biorthonormality can be determined as

0

Using the solution in terms of modal coordinates, given by Eq. (47), the time response of the system to an applied excitation can be determined. For instance, response of the system to an applied force in y direction at the location z ¼ zn can be obtained as follows:

Hiy;jy ¼

U trans ¼ Atrans ; y

ð51Þ

ð52Þ

ð53Þ ð54Þ

Note that U yr and hxr in Eqs. (51)–(54) are the normalized backward and forward eigenfunctions obtained using biorthonormality condition. Since the beam has free–free end conditions, there exist two rigid body modes where beam can translate or rotate without any elastic deformation. Thus the summation term in Eqs. (51)– (54) start from zero, including the rigid body mode contributions. Translational and rotational rigid body modes can be expressed as

In this study the end point FRFs of a multi segment shaft are calculated by using receptance coupling method, rather than classical approaches in which matrix sizes are increased by assembling the individual matrices. Graphical representation of receptance coupling of two beam segments is shown in Fig. 2. In order to apply receptance coupling to a rotating Timoshenko beam, unlike in non-rotating case, cross FRFs between two orthogonal planes should be considered as well. Thus based on the FRF definitions given by Eq. (50), and considering the response of the rotating Timoshenko beam at A1 end to applied force and moment excitations at the same point, receptance matrix can be defined as follows:

8 9 2 HA1x;A1x uA1x > > > > > = 6N A1y 6 A1x;A1y ¼6 > uA1y > > 4 HA1x;A1y > > > : ; NA1x;A1x hA1x

LA1y;A1x

HA1y;A1x

PA1y;A1y

NA1y;A1y

LA1y;A1y

HA1y;A1y

PA1y;A1x

NA1y;A1x

9 38 LA1x;A1x > f A1x > > > > > < = PA1x;A1y 7 7 MA1y 7 LA1x;A1y 5> f A1y > > > > > : ; PA1x;A1x MA1x

ð64Þ

In Eq. (64), subscripts indicate the location and direction of the excitation or response. For instance, in HA1y;A1x the first subscript A1x indicates that the force is applied at A1 location in y-direction and the second subscript A1y indicates the response considered is at A1 location in x direction.

123

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126

Fig. 2. Receptance coupling of two beams. Fig. 3. Including bearing dynamics to the rotor segment by structural modification.

Using the response and excitation relations given in Eq. (64), receptance matrices for the subassemblies A and B can be written as

 ½A ¼

½A11  ½A12  ½A21  ½A22 



 ½B ¼

½B11  ½B12 

 ð65Þ

½B21  ½B22 

where,

2

3

HA1x;A1x

LA1y;A1x

HA1y;A1x

LA1x;A1x

6N 6 A1x;A1y ½A11  ¼ 6 4 HA1x;A1y

PA1y;A1y LA1y;A1y

NA1y;A1y HA1y;A1y

PA1x;A1y 7 7 7 LA1x;A1y 5

N A1x;A1x

PA1y;A1x

NA1y;A1x

PA1x;A1x

ð66Þ

Also, ½A12  and ½A22  matrices can be constructed similar to ½A12  expression given in Eq. (66). Using the compatibility and continuity conditions at the connection points of subassemblies A and B, receptance coupling procedure can be applied and FRFs of the coupled structure C can be obtained as follows: 1

½C 11  ¼ ½A11   ½A12 ½½A22  þ ½B11  ½A21 

ð67Þ

½C 12  ¼ ½A12 ½½A22  þ ½B11 1 ½B12 

ð68Þ

1

½C 11  ¼ ½B21 ½½A22  þ ½B11  ½A21 

ð69Þ 1

½C 22  ¼ ½B22   ½B21 ½½A22  þ ½B11  ½B12 

ð70Þ

2

0 60 6 6 60 6 60 6 ½D ¼ 6 60 6 60 6 6 40 0

2

½HC1C1  ½HC1C2 

½LC1C1 

½LC1C2 

3

6 ½H 7 6 C2C1  ½HC2C2  ½LC2C1  ½LC2C2  7 ½aC  ¼ 6 7 4 ½NC1C1  ½NC1C2  ½PC1C1  ½PC1C2  5 ½NC2C1  ½NC2C2  ½PC2C1  ½PC2C2 

ð71Þ

where

½HC1C1  ¼ and



HC1x;C1x

HC1y;C1x

HC1x;C1y

HC1y;C1y

 ð72Þ

0 0 0 0 0 0 0 0

3 0 0 0 7 0 0 0 7 7 7 0 0 0 7 7 0 0 0 7 7 7 0 0 0 7 7 0 0 0 7   7 5 0 0 kh y þ i x c h y ðkhx þ ixchx Þ 0 0

ð73Þ

Then by using the structural modification method suggested by Ozguven [22], the receptance matrix of the modified system represented as C 0 can be obtained as

½aC0  ¼ ½½I þ ½aC ½D1 ½aC 

ð74Þ

However, as the nonzero part of the modifying matrix is smaller than the total dof, the size of the matrix that needs to be inverted can be reduced by partitioning the matrices [14] so that

h







 12 T

h



11 a11 D11 C 0 ¼ ½I þ aC



 21



ii1 

h 21

a11 C

h

 i

ð75Þ i 11

¼ aC 0 ¼ aC ½I  D11 aC 0  22   22   21 h 11 i 21  aC0 ¼ aC  aC D aC0

aC 0

h

i

"h

i

ð76Þ ð77Þ

where D11 is the submatrix of the modification matrix which can be represented as:

2.5. Including bearing dynamics In this study, the bearing dynamics is included into the system dynamics by using the structural modification method suggested by Ozguven [22]. In this method, the receptance matrix of a modified system is obtained from the receptance matrix of the original system and the modification matrix. Bearing dynamics can be added to the left end of the rotor shaft as shown in Fig. 3, and the remaining shaft segments can be coupled to the modified system with the receptance coupling procedure given in Section 2.4. In Fig. 3 ky is the translational stiffness, cy is the translational damping, kh is the rotational stiffness and ch is the rotational damping of the bearing. The receptance matrix of the unmodified system (shown in Fig. 3 as part C), and the modification matrix representing dynamic stiffness matrix of the bearing can be written in partitioned form as follows:

0 0 0 0 0 0   0 0 kx þ ixcy   0 0 ky þ ixcy 0 0 0 0 0 0 0 0 0 0 0 0

½D ¼

D11

½0

#

½0 ½0 3 2 ðkx þ ixcx Þ 0 0 0   7 6 0 ky þ ixcy 0 0 7 6   ½D11  ¼ 6 7 5 4 0 0 khy þ ixchy 0 ðkhx þ ixchx Þ 0 0 0

ð78Þ

ð79Þ

Note that to be able to partition ½D as given in Eq. (78), it will be necessary to renumber the dofs. 3. Numerical case study In this section a numerical case study is presented to verify the analytical model and the solution approach, as well as to demonstrate the effect of including gyroscopic moments on rotor dynamics. 3.1. Verification of the model Accurate determination of dynamic response of a rotor–bearing systems plays a crucial role in many applications such as jet engines, turbines and machine tools. For example; in machining centers consisting of multi-segment spindle, holder, tool and bearings, to avoid chatter vibrations which result poor surface finish, reduced productivity and even tool breakage, stability diagram of the cutting operation is needed. In order to generate stability

124

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126

As seen from Figs. 5–7, the analytical model proposed predicts the dynamics of the system accurately. The maximum frequency difference between the resonance frequencies calculated with ANSYS and the analytical model developed occurs in the sixth mode, and the difference is only 0.9%. For the remaining modes in the frequency range of interest, differences are less than 0.1%. In addition to its high accuracy, the proposed analytical approach is computationally very efficient as well. The time required for ANSYS solution is 1500 s whereas using the proposed analytical model the end point FRF of the same assembly can be determined in 102 s. 3.2. The effect of gyroscopic moment

Fig. 4. Multi-segment rotor (spindle–holder–tool assembly) supported by asymmetric bearings.

diagrams, FRF of the machining center at the tool tip should be determined. In a machining center, cutting operations can be performed with different spindle, holder and tool combinations and for each combination, tool point FRF should be determined. It is not practical to perform experimental modal analysis or employ FEA for the machining center for all holder and tool combinations. However; with the proposed analytical model, dynamics of the rotating system for each combination can be predicted with increased computational efficiency. First, in order to verify the accuracy and computational efficiency of the proposed analytical modeling approach, a multi segment spindle–holder–tool assembly with asymmetric bearings (Fig. 4) is modeled and analyzed with ANSYS, as well as with the proposed approach. For the ANSYS model, beam element BEAM 188 which is based on Timoshenko beam theory is used for the rotor segments, and asymmetric front and rear bearings are modeled using COMBIN 214 element. Frequency increment in the FRF calculation is selected as 0.5 Hz in both ANSYS and in the analytical modeling approach. Geometry of the assembly and bearing properties are given in Tables 1 and 2, respectively. Material properties of the assembly are taken as follows; mass density q = 7860 kg/m3, Young’s modulus E = 200 GPa, Poisson’s ratio m = 0.3 and structural damping ratio c ¼ 0:06. In Fig. 4, superscripts f and r correspond to the front and rear bearings, respectively. In order to check the accuracy of the proposed analytical model, end point FRFs of the assembly are calculated for the x–z and y–z planes with both ANSYS and the proposed model for 40 000 rpm spindle speed. In addition to the point FRFs, the cross FRFs of the assembly between two orthogonal planes are also calculated. The results are given in Figs. 5–7.

In order to study the effects of gyroscopic moment on the end point FRF of the spindle–holder–tool assembly system, the end point FRF is calculated at 40 000 rpm and also at idle conditions. The results are shown in Fig. 8. As seen from Fig. 8, gyroscopic effects cause the separation of the idle modes into backward and forward modes as expected, and the amount of separation between backward and forward modes increases with spindle speed. Although the effect of gyroscopic moment on frequency response may not be much, its further effects may be more important in some applications. In order to investigate the effects of the gyroscopic moments on the stability of the cutting process in this application, the end point FRF of the spindle–holder–tool assembly given in Table 1 is calculated at 20 000 rpm spindle speed, as well as at idle state. Then stability limits of the cutting process are calculated using the FRFs in operating and idle conditions by using the analytical milling chatter stability model proposed by Budak and Altintas [23]. In this application the workpiece is an aluminum alloy, and the radial depth of cut is 3 mm. The cutting force coefficients are taken as K t ¼ 625 MPa and K r ¼ 100 MPa [24]. The effect of the gyroscopic moments on the stability diagram is shown in Fig. 9 which demonstrates that rotational effects may increase the stability limits of the cutting operation so that an unstable point on the diagram may actually be stable. Note that stability diagram given in Fig. 9 is for 20 000 rpm spindle speed. The tool point FRFs are also calculated for spindle speeds 30 000 rpm and 40 000 rpm and corresponding stability limits are determined. The change in the stability limits for different spindle speeds are tabulated in Table 3. As can be seen from Table 3, if gyroscopic moment is not considered and FRFs of the idle state were used, the stability limits predicted would be significantly underpredicted. 3.3. Computational efficiency of the method In this section the computational efficiency of the proposed approach is demonstrated for the case study considered. In

Table 1 Multi-segment rotor (spindle–holder–tool assembly) dimensions. Spindle dimensions Segment number Length (mm) Outer diameter (mm)

1 25 65

2 115 100

3 57 120

4 13 110

5 20 100

6 230 90

7 30 70

8 33 60

9 57.5 55

10 6.5 50

Holder dimensions Segment number Length (mm) Outer diameter (mm)

1 25 70

2 110 55

Tool dimensions Segment number Length (mm) Outer diameter (mm) Inner diameter (mm)

1 45 18 0

2 55 22 0

125

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126 Table 2 Front and rear bearing properties. Front bearing

Rear bearing

kx (N/m)

7

5:5  10

5:5  107

ky (N/m)

7:5  107 10 10

7:5  107 10 10

cx (N s/m) cy (N s/m)

Analytical model

-4

H 1x,1x ( ω) [m/N]

10

ANSYS

Fig. 8. End point FRF of the multi segment beam in the x–z plane calculated by the analytical model for the idle state and for 40 000 rpm spindle speed.

-6

10

60

Idle condition Rotating condition

-8

0

500

1000

1500

2000

2500

3000

3500

4000

4500

Frequency [Hz] Fig. 5. End point FRF of the multi segment beam in the x–z plane calculated using ANSYS and the analytical model.

Axial depth of cut [mm]

10

50 40 30 20 10

-2

10

1

Analytical model ANSYS

1.5

H1y,1y(ω ) [m/N]

-4

2.5

3

3.5 4

x 10

10

Fig. 9. Stability diagrams predicted using the tool tip FRFs of the idle state and 20 000 rpm spindle speed.

-6

10

-8

10

Table 3 Variations in the stability limits for different spindle speeds.

-10

10

0

500

1000

1500

2000

2500

3000

3500

4000

4500

Frequency [Hz]

Spindle Speed (rpm)

Fig. 6. End point FRF of the multi segment beam in the y–z plane calculated using ANSYS and analytical model.

Analytical model ANSYS

H1x,1y(ω) [m/N]

2

Spindle speed [rpm]

20 25 35 40

000 000 000 000

Maximum axial depth of cut using (mm) Using FRFs at idle state

Using FRFs for rotating spindle

55 12 10 12

63 19 16 20

Difference (%)

14.5 58 60 67

-6

10

Table 4 Tool dimensions.

-8

10

0

500

1000

1500

2000

2500

3000

3500

Segment number Length (mm) Outer diameter (mm) Inner diameter (mm)

1 40 10 0

2 57 12 0

Frequency [Hz] Fig. 7. Cross FRFs between the two orthogonal planes obtained with ANSYS and analytical model.

engineering applications where modifications exist, FEM should be reconstructed and the dynamics of the whole system should be analyzed for each case. However, in the proposed approach, receptances of the unmodified part of the rotor bearing system can be stored and the system dynamics can be determined by analyzing the modified part of the system only. For example, in machine tool applications, the spindle remains the same in the analyses but different tool holders and cutting tools are clamped to the spindle during the operation. In order to demonstrate the advantage of

the proposed approach, the same case study is used and the spindle–holder–tool assembly given in Table 1 is reanalyzed by using a different tool whose properties are given in Table 4. Tool point FRFs are calculated both with ANSYS and the proposed analytical model. Note that in the analytical model, only the receptances of the tool segments are calculated and these receptances are coupled with the previously determined spindle–holder subassembly receptances. Therefore; while the determination of the modified system dynamics with ANSYS takes 1500 s, proposed method requires only 21 s for the computation of the end point FRFs of the modified system.

126

O. Özsßahin et al. / Computers and Structures 144 (2014) 119–126

4. Conclusion

Acknowledgement

In this study, for the dynamic response analysis of asymmetric multi-segment rotor – bearing systems including gyroscopic effects an analytical modeling procedure is proposed. The objective of the study is to determine FRFs of rotor segments by using analytical methods, and then employ FRF coupling and structural modification techniques to calculate the FRFs of the assembly. Rotor system is modeled using Timoshenko beam model and eigenfunctions of the rotating Timoshenko beam are determined for backward and forward motions analytically. Then, forced response of the beam is obtained using non-self adjoint system characteristics, and thus receptance matrix of the rotating Timoshenko beam which includes the effect of gyroscopic moment is determined. In order to calculate the FRFs of a multi segment beam, receptance coupling is applied and receptance matrix of the multi segment beam for the free–free end condition is obtained. Finally, bearing dynamics are included into the system by using a structural modification technique. The method proposed is verified by comparing the results obtained from a FEM of a rotor system using ANSYS with those obtained employing the method suggested. The results show that the method developed can provide highly accurate predictions with a considerable reduction in computational time. In addition to the accuracy and computational speed of the method proposed, the further advantage of using receptance coupling and modification techniques is demonstrated with a machine tool application in which modifications only in certain parts of the system are required. In such applications computational time is reduced further drastically. Finally, proposed analytical method is applied on a machining center by using the proposed method in order to determine tool point FRF and then stability of the cutting process is obtained for the rotating and idle conditions. Effects of the spindle rotation on the stability of the cutting process are examined, and it is illustrated that although the effect of gyroscopic moment on FRF may not be much, its effect on stability can be considerable. The contribution of this study to the literature can be summarized in three categories. First, forced response of a rotating Timoshenko beam is determined analytically in terms of frequency response functions. In addition to being more accurate, it requires much less computational time and effort compared to FEM. Second, instead of multi segment solution procedure, receptance coupling is used in the proposed response model. Therefore, the size of the matrices involved is very small (4  4) and the solution procedure is much easier. Third, when the bearings in two orthogonal planes are not symmetrical the solution procedure in classical approaches becomes complicated. However, in this study, bearing dynamics are included into the analysis by employing structural modification to the free–free solution of a symmetric beam. Thus, dynamics of an asymmetric rotor – bearing system can be determined by reducing the computational cost drastically, without losing accuracy.

This project is funded by the Scientific and Technological Research Council of Turkey (TUBITAK) under Project numbers 110M52 and 108M340 which is gratefully acknowledged. References [1] Nelson D. A finite rotating shaft element using Timoshenko beam theory. ASME J Mech Des 1980;102:793–803. [2] Lee CW. Vibration analysis of rotors. Dordrecht: Kluwer Academic Publishers; 1993. [3] Kang Y, Chang YP, Tsai JW, Mu LH, Chang YF. An investigation in stiffness effects on dynamics of rotor–bearing–foundation systems. J Sound Vib 2000;231:343–74. [4] Jang GH, Lee SH. Free vibration analysis of a spinning flexible disk spindle system supported by ball bearing and flexible shaft using the finite element method and substructure synthesis. J Sound Vib 2002;251(1):59–78. [5] Chatelet E, D’Ambrosio F, Richardet GJ. Toward global modelling approaches for dynamic analyses of rotating assemblies of turbomachines. J Sound Vib 2005;282:163–78. [6] Cavalca KL, Cavalcante PF, Okabe EP. An investigation on the influence of the supporting structure on the dynamics of the rotor system. Mech Syst Signal Process 2005;19:157–74. [7] Combescure D, Lazarus A. Refined finite element modelling for the vibration analysis of large rotating machines: application to the gas turbine modular helium reactor power conversion unit. J Sound Vib 2008;318:1262–80. [8] Friswell MI, Penny JET, Garvey SD, Lees AW. Dynamics of rotating machines. Cambridge aerospace series; 2010. [9] Frew DA, Scheffer C. Numerical modelling of a high-speed rigid rotor in a single-aerostatic bearing using modified Euler equations of motion. Mech Syst Signal Process 2008;22:133–54. [10] Chena YS, Chenga YD, Yang T, Koai KL. Accurate identification of the frequency response functions for the rotor–bearing–foundation system using the modified pseudo mode shape method. J Sound Vib 2010;329:644–58. [11] Lee CW, Katz R, Ulsoy AG, Scott RA. Modal analysis of a distributed parameter rotating shaft. J Sound Vib 1988;122:119–30. [12] Lee CW, Jei YG. Modal analysis of continuous rotor–bearing systems. J Sound Vib 1988;126:345–61. [13] Jei YG, Lee CW. Modal analysis of continuous asymmetric rotor–bearing systems. J Sound Vib 1992;152:245–62. [14] Wang W, Kirkhope J. New eigensolutions and modal analysis for gyroscopic/ rotor systems. Part 1: Undamped systems. J Sound Vib 1994;175:159–70. [15] Wang W, Kirkhope J. New eigensolutions and modal analysis for gyroscopic/ rotor systems. Part 2: Perturbation analysis for damped systems. J Sound Vib 1994;175:171–83. [16] Parker RG, Sathe PJ. Exact Solutions for the free and forced vibration of a rotating disk–spindle system. J Sound Vib 1999;223:444–65. [17] Schmitz T, Donaldson R. Predicting high-speed machining dynamics by substructure analysis. Ann CIRP 2000;49(1):303–8. [18] Ertürk A, Özgüven HN, Budak E. Analytical modeling of spindle–tool dynamics on machine tools using Timoshenko beam model and receptance coupling for the prediction of tool point FRF. Int J Mach Tools Manuf 2006;46:1901–12. [19] Hong SW, Park JH. Dynamic analysis of multi-stepped, distributed parameter rotor–bearing systems. J Sound Vib 1999;227:769–85. [20] Aristizabal-Ochoa JD. Timoshenko beam-column with generalized end conditions and nonclassical modes of vibration of shear beams. J Eng Mech 2004;130:1151–9. [21] Ewins DJ. Modal testing: theory, practice and application. Baldock (Hertfordshire, England): Research Studies Press; 2000. [22] Ozguven HN. Structural modifications using frequency response functions. Mech Syst Signal Process 1990;4:53–63. [23] Budak E, Altintas Y. Analytical prediction of chatter stability in milling – Part I: General formulation; Part II: Application to common milling systems. Trans ASME J Dyn Syst Measure Control 1998;120:22–36. [24] Altintas Y. Manufacturing automation: metal cutting mechanics, machine tool vibrations and CNC design. Cambridge University Press; 2000.