Hybrid modelling for analysis and identification of rotors

Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

www.elsevier.com/locate/cma

Hybrid modelling for analysis and identi®cation of rotors M. Aleyaasin, M. Ebrahimi * Department of Mechanical and Medical Engineering, University of Bradford, West-Yorkshire BD7 1DP, UK Received 17 October 1998; received in revised form 7 January 1999

Abstract In this paper a rotating shaft-disc system, modelled as a series of interconnected, distributed and lumped elements is considered. The overall transfer matrix of the system formed by series ¯exible distributed elements and rigid disc elements is determined, then applying supporting end boundary conditions the characteristic determinant of the system is computed. The damped natural frequency of the system is computed by a direct determinant search method. Thereafter a novel method is proposed for the determination of the frequency response function. The overall system is identi®ed as a linear multivariable dynamic system with second order transfer function elements. Future applications of this hybrid technique are discussed. Ó 2000 Elsevier Science S.A. All rights reserved.

Notations q A E I s C L L Yi L ,Y i R

Yi R ,Y i L hi L ,hi

R hi ,hi L Mi L ,M i R Mi R ,M i L Qi L ,Qi R Qi R ,Qi R

T …i† ‰Ul Š …i† ‰Ud Š X U,V W,K

*

density of the distributed shaft area of the distributed shaft cross section distributed shaft modulus second moment of area in bending Laplace transform variable auxiliary Laplace variable length of each distributed part displacement function in the left of the ith lumped and distributed element displacement function in the right of the ith lumped and distributed element slope function in the left of the ith lumped and distributed element slope function in the right of the ith lumped and distributed element bending moment function in the left of the ith lumped and distributed element bending moment function in the right of the ith lumped and distributed element shear force function in the left of the ith lumped and distributed element shear force function in the right of the ith lumped and distributed element overall transfer matrix transfer matrix of the ith lumped element transfer matrix of the ith distributed element rotational speed of the shaft decomposed transfer matrix inverse of the decomposed transfer matrices

Corresponding author. Tel.: +44-1274-3845-26; fax: +44-1274-3845-25. E-mail address: [email protected] (M. Ebrahimi).

0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 0 9 0 - 0

164

Gij f m,c,k x,y

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

frequency response function external force function equivalent mass, damping, sti€ness, respectively real and imaginary part of the frequency response function

1. Introduction The analysis of free ¯exural vibrations in rotating shaft-disc and bearing systems is an important aspect in the design of rotating machinery. Accurate methods for computing the damped natural frequencies and modal responses in gas turbine and other high-speed rotor systems, for example, provides typical application studies. Nelson's [1] contribution to the dynamic rotor bearing system analysis, using ®nite element methods, resulted in a high number of eigenvalues which obtained through domain discretization. This led to development of methods whereby erroneous eigenvalues and modes could be neglected. Many techniques for this process exist now such as dynamic reduction by Rouch [2], Modal Transformation by Kim and Lee [3] and Component Mode Synthesis by Glasgow and Nelson [4] etc. As an alternative for ®nite element method (FEM), the Transfer Matrix Method (TMM) which is analogous to Myklestad±Prohl method, has been used by engineers [5] to avoid such high number of eigenvalues obtained through FEM. The more exact results could be achieved by using TMM, through computation of the successive transfer matrices for each part a system. For checking the results of the FEM, these could be compared with TMM results [6,7]. For the determination of the modal responses, the transfer matrix method [8] is limited in application especially for the cases when vibration control of the rotor bearing system is considered [9]. Here frequency response data is needed for the estimation of the rotor-bearing system parameter [10]. Recent strategies for the active control of rotor vibration [11], has demonstrated that an accurate knowledge of the system parameters is needed for controller design studies. Therefore the identi®cation of rotor-bearing systems and validation via simulation is a necessary prelude to controller design. Some e€orts have been made enabling the accurate analysis of the ¯exural vibration of rotors. For example Yang and Pilkey [12] analyzed rotor-bearing systems using the governing partial di€erential equations. A truncated analysis method using frequency dependent mass, sti€ness and damping matrices was thereafter presented [13]. In this paper the approach employed is based on the theory of Distributed-Lumped parameter systems as shown in Whalley [14] which has also been applied to transmission line modeling. Further investigations by Barttlet et al. [15,16] applied this method to unsteady gas ¯ow behavior in pipes and tunnels. More recently the transient, torsional vibrations in power transmission systems is analyzed in which the e€ect of second order wave propagation, partial di€erential equations was evaluated by introducing Distributed, Lumped and Termination elements. Herein a rotor system is considered as a series of ¯exible Distributed Elements which are connected together by rigid discs forming the Lumped Elements. Laplace transform techniques are applied to the fourth order partial di€erential equations representing the ¯exural vibration of the system. Transfer matrices for each distributed element are derived. Upon multiplied by the transfer matrix for the lumped elements, in a manner analogous to Myklestad±Prohl method, the damped natural frequency of the system may be computed. A novel method is proposed for the determination of the frequency response of the system through an overall transfer matrix decomposition. After obtaining the frequency response data, the identi®cation procedure can be executed enabling the mass, sti€ness and damping coecients to be determined using a peak and eigenvalue matching algorithm. A numerical example based on a rotor consisting of three discs connected by four distributed elements on simply supported bearings is considered. The natural frequency of the system is computed and the results are compared with classical methods. After determination of the frequency response data the rotor model is identi®ed in the form of a multivariable dynamic system, having a 3 ´ 3 transfer function matrix

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

165

description. For each element of the transfer function matrix, the mass, sti€ness and damping coecients due to the gyroscopic e€ects, generated by the lumped discs, are determined.

2. Analysis In Fig. 1, a rotating shaft representation consisting of series of distributed and lumped elements is shown. The number of lumped elements is n and there are n + 1 distributed elements. The Distributed-Lumped parameter model for the system is shown in Fig. 2. For the distributed part there are left and right terminations each of which have four parameters. These are vertical displacement, L R L R L R L R Y i ; Y i slopes hi ; hi , bending moments M i ; M i and shear forces Qi ; Qi for the left and right boundary terminations, respectively. For the lumped part there are also left and right terminations each speci®ed by the above four parameters. The input to the lumped element is the output from the distributed element. The procedures for the derivation of an Euler beam model are exercised for the distributed part while the equations for rigid body motion are considered for the lumped [17] model. The following relations exist for the input and output of each distributed, and lumped element the detailed analysis of which are mentioned in Ref. [17]. The governing distributed element equation is 2 R 3 2 L3 Yi Yi 7 6 R 7 6 6 h 7 h i 6 hL 7 6 i 7 6 i 7 …i† …1† 6 R 7 ˆ Ud
6 L 7 ; 6M 7 6M 7 4 i 5 4 i 5 R

Qi

L

Qi

h i …i† where Ud is the 4 ´ 4 matrix

…2†

Fig. 1. A Distributed-Lumped element rotating shaft.

166

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

Fig. 2. The Distributed-Lumped parameter model.

For the lumped part 2 R 3 2 L3 Yi Yi 6 hR 7 h i 6 hL 7 6 i 7 6 7 …i† 6 R 7 ˆ Ul
6 i L 7; 4 Mi 5 4 Mi 5 QR QLi i 2 1 h i 6 0 …i† 6 Ul ˆ 4 0 ms2 where Cˆs

0 1 sXJp 0

0 0 1 0

…3†

3 0 07 7; 05 1

…4†

p L0 C0 ;

…5†

qA ˆ L0 ;

…6†

1 ˆ C0 ; EI

…7†

In Eqs. (1)±(7) s is the Laplace transform variable, m mass of each lumped element, X the rotational speed of the shaft and q the density of the shaft material, A the cross sectional area, E the modulus of elasticity and I the moment of inertia of the cross section in bending, for each distributed element. The polar moment of inertia of the lumped disc is Jp and L is the length of each distributed element. L R L R L L R R According to Fig. 2, Y 1 ; h1 ; M 1 ; Q1 are the values for the left end and Y n‡1 ; hn‡1 ; M n‡1 ; Qn‡1 the values for …i† the right end. If the transfer matrix of the ith distributed element Di is ‰Ud Š and the transfer matrix of the …i† ith lumped element Li is ‰Ul Š, then the transfer matrix of the overall system T relating the above mentioned values for the left and right end could be computed as follows. 2 R 3 2 L3 Y1 Y n‡1 6 R 7 6 L 7 7 7 6h 6 6 n‡1 7 ˆ ‰T Š
6 h1 7; …8† 7 6 6 L7 4 M n‡1 5 4 M1 5 L Qn‡1 Q 1

where in Eq. (8) h i …n‡1†
‰T Š ˆ Ud

nÿ1 h Y iˆ0

…nÿi† Ud

! ih i …nÿi† : Ul :

…9†

Finally the Eq. (8) could be expanded according to the elements of the T matrix yielding. R

L

L

L

L

…10†

R

L

L

L

L

…11†

Y n‡1 ˆ T11 Y 1 ‡ T12 h1 ‡ T13 M 1 ‡ T14 Q1 ; hn‡1 ˆ T21 Y 1 ‡ T22 h1 ‡ T23 M 1 ‡ T24 Q1 ;

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176 R

L

L

167

L

L

M n‡1 ˆ T31 Y 1 ‡ T32 h1 ‡ T33 M 1 ‡ T34 Q1 ; R

L

L

…12†

L

L

Qn‡1 ˆ T41 Y 1 ‡ T42 h1 ‡ T43 M 1 ‡ T44 Q1 :

…13† L

L

For the con®guration in Fig. 2 the values of M 1 ˆ 0; Y 1 ˆ 0 are for the left end distributed part and the R R values of M n‡1 ˆ 0; Y n‡1 ˆ 0 are for the right end distributed part. Therefore Eqs. (10) and (12) are simply. L

L

…14†

L

L

…15†

T12 h1 ‡ T14 Q1 ˆ 0; T32 h1 ‡ T34 Q1 ˆ 0:

In order to obtain non-trivial solution for Eqs. (14) and (15) the following determinant must be zero. T12 T14 ˆ 0: …16† R ˆ det T32 T34 The damped natural frequency of the system could be computed by determining the roots of the irrational characteristic Eq. (16). This procedure will be shown via a worked example in the next section. 3. Determination of the frequency response The overall transfer matrix of the rotor T in Eq. (9) is decomposed into two parts, the pth lumped station Lp , so that the ®rst part of the transfer matrix U is de®ned as follows: Up ˆ

pÿ1 h Y iˆ0

…pÿi†

Ul

ih i …pÿi† : Ud :

Then we can write the 2 R3 2 Yp Up 6 hR 7 6 U 11 6 p 7 6 p21 6 R7ˆ4 Up31 4 Mp 5 R Up41 Qp

…17†

following relation for the matrix U 2 L3 3 Y Up12 Up13 Up14 6 L1 7 6 7 Up22 Up23 Up24 7 7
6 h1 7 : L7 Up32 Up33 Up34 5 6 4 M1 5 L Up42 Up43 Up44 Q1

…18†

The second part of the transfer matrix is de®ned for the rest of the system as ! h i nÿpÿ1 Y h …nÿi† i h …nÿi† i …n‡1† Ul
: Ud : Vp ˆ U d

…19†

iˆ0

By de®ning the inverse of the matrix Vp further relations could be written Wp ˆ Vpÿ1 ; 2

L

Y p‡1

3

…20† 2

Wp11 6 L 7 6 h 7 6W 6 p‡1 7 6 p21 6 L 7ˆ4 Wp31 6 M 7 4 p 5 Wp41 L Q0 p

Wp12 Wp22 Wp32 Wp42

Wp13 Wp23 Wp33 Wp43

2 R 3 3 Y n‡1 Wp14 6 R 7 7 6 Wp24 7 6 hn‡1 7 7
R 7: Wp34 5 6 4 M n‡1 5 R Wp44 Qn‡1

…21†

L

L

For the systems of Eqs. (18) and (19) for the simply supported end conditions the values M 1 ˆ 0; Y 1 ˆ 0, R R M n‡1 ˆ 0; Y n‡1 ˆ 0 so we can write the following eight equations. L

L

YpR ˆ Up12 h1 ‡ Up14 Q1 ;

…22†

168

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176 L

L

hR p ˆ Up22 h1 ‡ Up24 Q1 ; L

…23†

L

MpR ˆ Up32 h1 ‡ Up34 Q1 ;

…24†

L

…25†

L

QR p ˆ Up42 h1 ‡ Up44 Q1 ; L

R

R

…26†

hp‡1 ˆ Wp22 hn‡1 ‡ Wp24 Qn‡1 ;

L

R

R

…27†

M p‡1 ˆ Wp32 hn‡1 ‡ Wp34 Qn‡1 ;

L

R

R

…28†

L

R

R

…29†

Y p‡1 ˆ Wp12 hn‡1 ‡ Wp14 Qn‡1 ;

Q0 p‡1 ˆ Wp42 hn‡1 ‡ Wp44 Qn‡1 :

According to Fig. 1 by leaving the pth lumped element at the left we will enter the p + 1th distributed element at the right, so we can write. L

L

hR p ˆ hp‡1 ;

YpR ˆ Y p‡1 ;

L

MpR ˆ M p‡1 : L

instead of Upon substituting, the external force fp into transfer matrix we can write Qp‡1 ˆ QR p , but L considering the external force fp at the pth Lumped element transfer matrix, we introduce Q0 p‡1 6ˆ QR p and the external force fp could be expressed as follows, this enables, driving an easy algorithm for determination of the frequency response. L

0 fp ˆ QR p ÿ Q p‡1 :

…30†

The Eqs. (22)±(30) could be written in the matrix form as follows 3 0 607 6 7 607 6 7 6 7 07 iT   6 6 7 L 7 6 0 ˆ K
Q p‡1 p 6 0 7: 607 6 7 6 7 607 6 7 405 fp

…31†

3ÿ1 0 0 7 7 0 7 7 7 0 7 7 0 7 7 : 0 7 7 7 0 7 7 0 5 ÿ1

…32†

2

h

L

Y p‡1

L

hp‡1

L

M p‡1

QR p

L

h1

L

Q1

R

hn‡1

R

Qn‡1

The 9 ´ 9 matrix Kp is as the following form: 2

1 60 6 60 6 6 60 6 Kp ˆ 6 61 60 6 6 60 6 40 0

0 1 0 0 0 1 0 0 0

0 0 1 0 0 0 1 0 0

0 0 0 1 0 0 0 1 1

ÿ Up12 ÿ Up22 ÿ Up32 ÿ Up42 0 0 0 0 0

ÿ Up14 ÿ Up24 ÿ Up34 ÿ Up44 0 0 0 0 0

0 0 0 0 ÿ Wp12 ÿ Wp22 ÿ Wp32 ÿ Wp42 0

0 0 0 0 ÿ Wp14 ÿ Wp24 ÿ Wp34 ÿ Wp44 0

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

169

The input force fp is applied in the pth lumped element and the displacement Yp results there, so the desired transfer function could be obtained Gpp …s† ˆ

Yp …C† ˆ Kp19 …C†: fp …C†

…33†

For computation of the irrational transfer function given by Eq. (33), ®rst we use relations (2) and (5), then following Eqs. (17), (19) and (20) and ®nally inversion according to Eq. (32) is undertaken. By putting s ˆ ix in Eq. (33) the desired frequency response function can be computed. In Section 4 this procedure is executed numerically. The displacement in the qth Lumped element due to the external force applied at pth element could be computed by de®ning the matrix. nÿq h i h i Y …nÿi† …nÿi† q Ud
Ul : …34† ‰HŠp ˆ iˆnÿpÿ1

The H matrix in Eq. (34) has the 2 R3 2 Yq Hpq Hpq Hpq 12 13 6 R 7 6 11 6 h 7 6 Hpq Hpq Hpq 6 q 7 6 21 22 23 6 R 7 ˆ 6 pq pq 6 Mq 7 4 H31 Hpq H 32 33 5 4 Hpq 41

QR q

Hpq 42

Hpq 43

following form 2 L 3 Y p‡1 pq 3 H14 6 L 7 7 7 6 7 6h Hpq 24 7 7
6 p‡1 7: pq 7 6 L 7 H34 5 6 M 7 4 p‡1 5 Hpq L 44 Q0 p‡1

…35†

The YqR that represents the frequency response could be computed by expanding Eq. (35), and further substitution from Eq. (31), enabling the relevant in¯uence transfer function to be expressed as. pq pq pq Gpq …s† ˆ Kp19 …C†
Hpq 11 …C† ‡ Kp29 …C†
H12 …C† ‡ Kp39 …C†
H13 …C† ‡ Kp49 …C†
H14 …C†:

…36†

4. Frequency domain identi®cation In Eqs. (33) and (36), by putting s ˆ ix, the frequency response data could be obtained from each transfer function, so that a vector of frequency response data as a function of x could be written as: Gpp …i
x† ˆ xpp ‡ i ypp ;

…37†

Gpq …i
x† ˆ xpq ‡ i ypq :

…38†

Although there are some techniques for the identi®cation of models from frequency response data [18], state space realisations are usually employed. If each set of the frequency response data is considered separately, then we can estimate any rational transfer function [19,20]. The existence of unstable poles in the resulting high order transfer function may occur. Since the results obtained through simulation do not include noise, the authors have investigated the time response from the frequency response data using the inverse Fourier transform [21]. The results show the desired response similar to a second order system model. This leads to the application of a second order identi®cation algorithm [21]. Therefore be employing the identi®cation technique, the actual transfer functions from the frequency response data in Eqs. (37) and (38) can be determined, which have the following form: G0pp …s† ˆ G0pq …s† ˆ

1 ; ‡ cpp s ‡ kpp

…39†

1 : mpq s2 ‡ cpq s ‡ kpq

…40†

mpp

s2

170

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

For the mass, damping, and sti€ness coecients to have the same roots as the determinant in Eqs. (16), (39) and (40) can be written as: G0pp …s† ˆ

1 ; mpp …s2 ‡ 2as ‡ a2 ‡ b2 †

…41†

G0pq …s† ˆ

1 ; mpq …s2 ‡ 2as ‡ a2 ‡ b2 †

…42†

where s ˆ ÿa  ib are the zeros of the characteristic determinant expressed in Eq. (16). …max† …max† and ypq are the maximum absolute value of the frequency response function of the actual If ypp system, in Eqs. (37) and (38), then the second order systems models described in Eqs. (39) and (40) have the same maximum absolute value, subject to the following equations [23]: …max† ˆ ypp

2mpp q ; cpp 4kpp mpp ÿ c2pp

…43†

…max† ˆ ypq

2mpq q : cpq 4kpq mpq ÿ c2pq

…44†

By comparing Eqs. (39) and (40) with Eqs. (41) and (42) the following equations can be obtained. cpp ˆ 2ampp ;

…45†

kpp ˆ mpp …a2 ‡ b2 †;

…46†

cpq ˆ 2ampq ;

…47†

kpq ˆ mpq …a2 ‡ b2 †:

…48†

By substituting Eqs. (45)±(48) into Eqs. (43) and (44) the following relations can be obtained. mpp ˆ

mpq ˆ

1 …max†

2abypp 1

…max†

2abypq

;

…49†

:

…50†

The parameters of the second order models have the same maximum absolute value as the frequency response function and also the poles are equal to the roots of the determinant in Eq. (16). As we will see in the next section this method of identi®cation leads to the optimal matching of the actual and estimated frequency response data. 5. Numerical example A rotating shaft on simply supported bearings has a 1.6 m span and rotational speed of 5000 rpm and carries three equi-spaced discs, each with mass of 100 kg The polar moment of each disc is Jp ˆ 5:1011 kgm2 and the modulus of the shaft is E ˆ 200 GN=m2 with a diameter of 8 cm . As shown in Fig. 1, the three lumped discs and four distributed ¯exible elements each with length of L ˆ 40 cm, form the complete system. The following procedure is used for computing the natural frequency of the system. 2 =4 1. Compute I ˆ pD4 =64 and then C0 ˆ 1=EI and also L0 ˆ qA ˆ pqDp  for each distributed part. 2. Choose a value for x and put s ˆ ix ‡ r in Eq. (5) then compute C for each distributed part.

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

3. 4. 5. 6.

171

Using the values of L; C; C0 compute the matrix ‰Ud Š for each distributed part from relation (2) Using the values of X; ix; Jp ; Jd ; m compute the matrix ‰Ul Š for each lumped part from relation (4) Find the elements of matrix T from Eq. (9) and then compute the value of R from the relation (16). Continue the steps 2±5 by choosing further values x and thereafter compute R. When R is minimized the natural frequency has been determined. For the above problem, according to Fig. 3, the minimum value of R in Eq. (16) is obtained by increasing the value of x to x ˆ 165 rad/s where the ®rst natural frequency of the system occurs. Here the value of R is minimized and the gyroscopic forces prevent R becoming zero [22]. The minimum value of the determinant is R ˆ 0.3898 for x ˆ 165 rad/sec where there is a peak in the frequency response curve. Due to gyroscopic effects the actual root of the determinant has a negative real

Fig. 3. Values of R versus angular frequencies.

Fig. 4. Flexural damped natural frequency.

172

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

part. This is equal to the damped natural frequency of the system. In Fig. 4 the 3-D graph, shows the value of the determinant R versus real and imaginary part of the root, r and x. The minimum value of the determinant is R ˆ 0.0075 for s ˆ ÿ18 + 179i which determines the damped natural frequency of the system so that x ˆ 179 ad/s. There is no peak in the frequency response curve at this frequency. This is due to the high inertia of the attached discs, which introduce gyroscopic moments, providing a damping ratio of approximately f ˆ 0:1. Generally for determination of the zeros in an irrational function, like the characteristic determinant R in Eq. (16), various optimization techniques could be implemented [24]. In this article in order to observe that the minimum of the determinant is global or local, we have drawn the 3-D graph like Fig. 4. By using this method, determination of the frequency in which the global minimum for the value of R in Eq. (16), occurs, is not dif®cult when a 2-D array for the absolute value of the R is available. For example in Fig. 4 the minimum value of R, and also the location of it, is founded by direct search from an 400 ´ 400 array, some modern softwareÕs have such in built commands [25]. The computational eciency is higher with this method than with Finite Element methods where many additional eigenvalues are introduced [1]. The natural frequency of the above system can be computed through Dunkerely`s method which excludes the e€ect of the gyroscopic moments. The natural frequency is computed from this to be xn ˆ 152.5 rad/s. when, we compare with 165 rad/s in Fig. 3 and 179 rad/s in Fig. 4, evidently the di€erence exists because DunkerelyÕs method has higher approximation than our accurate method, and also our approach considers the gyroscopic e€ects. Furthermore , even if we increase the range of the frequency, we can see only one minimum in the determinant as shown in Fig. 4, so the ability of this new method to capture the fundamental modes of vibration could be concluded. This also explains, why the higher natural frequency values have not been reported in experimental frequency response tests on aircraft gas turbines [6,7]. For frequency response determination, after the computation of ‰Ud …ix†Š for each distributed element and ‰Ul …ix†Š for each lumped element, the computation procedure is as follows: Step 1: compute Up …ix† from Eq. (17), Vp …ix† from Eq. (19) and Wp …ix† from Eq. (21). Step 2: compute Kp …ix† from Eq. (32) and then Gpp …ix† from Eq. (33). Step 3: when response in the lumped part q due to the excitation in lumped part p is desired, compute q ‰ H Šp …ix† from Eq. (34) and then Gpq …ix† from Eq. (36). For the system mentioned the frequency response curves can be drawn following the above steps for the range of 0±300 rad/s as in Figs. 5 and 6. The results obtained through Fig. 5 show that when x ˆ178 rad/s there is a peak in the frequency response curve, which con®rms the computed damped natural frequency of the system. The functions G11 …ix†; G12 …ix†; G13 …ix† are shown in Fig. 5. The ®rst sux in G is for excitation, and the second one is for

Fig. 5. Frequency response plots for G11, G12 and G13.

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

173

Fig. 6. Frequency response plot for G22 and G23.

the response. The three curves indicate that the system appears to be second order, as Fig. 4 shows, with only one pair of eigenvalues. As mentioned, by using inverse Fourier transform [21], the results show the desired response of a second order mechanical system. This leads to the application of the second order identi®cation algorithm to the frequency domain data in Figs. 5 and 6. Due to symmetry only the frequency response curves G22 …ix†; G23 …ix†, are drawn in Fig. 6, and again at x ˆ178 rad/s, the peak in the FRF curves have second order shape also. In these curves the vertical axis shows the displacement in microns due to a 1N sinusoidal force. Frequency domain identi®cation for the above example is based on the transfer functions G0pp …s†; G0pq …s† as a second order system estimation with the eigenvalues and peak frequency response values mentioned in Eqs. (39) and (40). The system is identi®ed as a three input, three output multivariable system. The identi®cation is performed by using the system eigenvalues as shown in Fig. 4 for which the parameters are a ˆ 18 and b ˆ 179, and thereafter by ®nding the maximum absolute value of the frequency response curves as indicated in Figs. 5 and 6. Finally mpp ; cpp ; kpp and mpq ; cpq ; kpq can be computed from Eqs. (45)±(50). With f1 ; f2 ; f3 as the excitation forces at discs 1, 2, 3 the displacement response, from the transfer function matrix is given by 2 3 2 0 3 3 2 G11 …s† G012 …s† G013 …s† F1 …s† Y1 …s† 4 Y2 …s† 5 ˆ 4 G021 …s† G022 …s† G023 …s† 5
4 F2 …s† 5; …51† G031 …s† G032 …s† G033 …s† Y3 …s† F3 …s† where G0pq …s† ˆ

mpq

s2

1 ‡ cpq s ‡ kpq

p; q ˆ 1; 2; 3:

Each transfer function of the matrix in Eq. (51) is as follows. G011 …s† ˆ G033 …s† ˆ G022 …s† ˆ

163s2

1 ; 309:9s2 ‡ 10493s ‡ 9837900

1 ; ‡ 5868s ‡ 5275700

G012 …s† ˆ G021 …s† ˆ

218:4s2

1 ; ‡ 7865:6s ‡ 7071400

…52†

174

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

G012 …s† ˆ G021 …s† ˆ G013 …s† ˆ G031 …s† ˆ G023 …s† ˆ G032 …s† ˆ

218:4s2

1 ; ‡ 7865:6s ‡ 7071400

240:4s2

1 ; ‡ 8654:5s ‡ 7780600

1 : 233:9s2 ‡ 8423:7s ‡ 7573100

The estimated, second order transfer functions, not only have the same eigenvalues as the system model, but also have the closest frequency response curve. In Fig. 7 the actual frequency response G11 …ix† and the estimated second order frequency response G011 …ix† are shown. These traces are close to each other and the actual response curve overlaps the estimated function plot. The third curve or G0011 …ix† is also drawn using a

Fig. 7. Actual, estimated and time delay frequency response plot for G11.

Fig. 8. Actual and estimated frequency response plot of for G12.

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

175

Fig. 9. Bode diagram for actual G11 and estimated 2nd order G11.

Fig. 10. Phase diagram for actual G11 and estimated 2nd order G11.

second order transfer function estimation with an optimal time delay. Regardless of the computing technique, this transfer function is mismatched with the actual frequency response function and G011 …ix† is very close to G11 (ix) only. In Fig. 8 G012 …ix†and G12 …ix† are shown and the solid line indicates the actual transfer function, con®rming the ef®ciency of this algorithm. In Figs. 9 and 10 the Bode diagram of the actual system and the estimated second order system are shown, for a range of operational frequencies of (0±400) rad/s. The response of the both actual and estimated systems are the same.

6. Conclusion The frequency response of the system model indicated in Figs. 7 and 8, has been computed. Thereafter, a second order modelling approach based on its frequency response function has been closely matched with

176

M. Aleyaasin, M. Ebrahimi / Comput. Methods Appl. Mech. Engrg. 182 (2000) 163±176

the actual model frequency response. This is important in the case of rotor-disc dynamics and ¯uid ®lm bearing systems in which the amplitude of vibration could be controlled using magnetic bearings. The transfer function matrix given by Eqs. (51) and (52) now enables a control strategy to be formulated and speci®c design strategies considered. References [1] H.D. Nelson, T. McVaugh, The dynamics of rotor bearing systems using ®nite elements, Trans. ASME J. Engrg. Industry 98 (1976) 593±600. [2] K.E. Rouch, J.S. Kao, Dynamic reduction in rotor dynamics by the ®nite element method, Trans. ASME J. Mech. Design 102 (1980) 360±367. [3] Y.D. Kim, C.W. Lee, Finite element analysis of rotor bearing systems using a modal transformation matrix, J. Sound Vibration 111 (3) (1986) 441±456. [4] D.A. Glasgow, H.D. Nelson, Stability analysis of rotor-bearing systems using component mode synthesis, Trans. ASME J. Mech. Design 102 (1980) 353±359. [5] J.W. Lund, Stability and damped critical speeds of a ¯exible rotor in ¯uid-®lm bearings, Trans. ASME J. Engrg. Industry May (1974), 509±517. [6] J.M. Vance, A.C. Royal, High speed rotor dynamics-an assessment of current technology for small turboshaft engines, J. Aircraft, 12 (4) (1975). [7] D.G., Hibner, Dynamic response of viscous-damped multishaft jet Engines, J. Aircraft 12 (4) (1975). [8] J.W. Lund, Modal response of a ¯exible rotor in ¯uid-®lm bearings, Trans. ASME J. Engrg. for Industry, May (1974), 525±533. [9] C.R. Burrows, M.N. Sahinkaya, Vibration control of multi-mode rotor-bearing systems, Proc. Roy. Soc. London A 386 (1983) 77±94. [10] C.R. Burrows, M.N. Sahinkaya, Parameter estimation of multi-mode rotor-bearing systems, Proc. Roy. Soc. London A 379 (1982) 367±387. [11] J.Y. Huang, Magnetic bearing control using fuzzy logic, IEEE Trans. Industry Appl. 31 (6) (1995) 1492±1497. [12] B.S. Yang, W.D. Pilkey, Accurate free vibration analysis for rotating shafts, ASME DE Rotating Machinery and Vehicle Dynamics 35 (1991) 133±138. [13] B.S. Yang, W.D. Pilkey, N.J. Fergusson, Frequency dependent element matrices for rotor dynamics analysis, J. Sound Vibration 159 (2) (1992) 339±351. [14] R. Whalley, The response of Distributed-Lumped parameter systems, I Mech E. Proc. 202 (C6) (1988) 421±429. [15] H. Bartlett, R. Whalley, Gas ¯ow in pipes and tunnels, Proc. Inst. Mech. Engrg. Part I 209 (1995) 41±52. [16] H. Bartlett, R. Whalley, S.S.I Rizvi, Hybrid modelling of marine power transmission systems, in: Proceedings of MBDMST Conference, MEP Publication, UK, 1997, pp. 97±108. [17] Hybrid analysis of the ¯exural vibration of a rotating shaft internal report, Department of Mechanical Engineering, University of Bradford, January, 1998. [18] C.W. Chen, J.N. Juang, G. Lee, Frequency domain state-space system identi®cation, Trans. ASME J. Vibration Acoustics 116 (1994) 523±528. [19] Multivariable Frequency Domain Toolbox For Use with Matlab, Cambridge Control, 1998. [20] T. McKelvey, Identi®cation of State Space Models from Time and Frequency data, PhD thesis 380, Linkoping University, 1995. [21] The Response of Distributed-Lumped Systems Using Inverse Fourier Transforms internal report, Department of Mechanical Engineering University of Bradford, April, 1998. [22] C.W. Lee, A.G. Katz, R.A. Scott, Modal analysis of a distributed parameter rotating shaft, J. Sound Vibration 122 (1988) 119±130. [23] C.F. Beards, Engineering Vibration Analysis with Application to Control Systems, Edward Arnold, London, 1995. [24] Masanao, A, Introduction to optimization techniques: Fundamentals and applications of nonlinear programming, Macmillan Series in Applied Computer Sciences, 1971. [25] Matlab, The Mathworks, Natick, MA, 1997.