Unbalance response of rotors: A modal approach with some extensions to damped natural systems

Journal of Sound and Vibrarion

UNBALANCE

(1990) 140(l), 129- 153

RESPONSE

APPROACH

WITH

DAMPED

OF ROTORS:

SOME

EXTENSIONS

NATURAL

G. GENTA AND

A MODAL TO

SYSTEMS

F. DE BONA:

Dipartimento di Meccanica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10100 Torino, Italy (Received 10 March 1989, and injnalform

13 September 1989)

A modal procedure for the computation of the unbalance response of rotors modelled as systems with many degrees of freedom (possibly by using the finite element method) is suggested. If damping is neglected, the procedure is quite similar to that usually followed for a non-rotating structure, the only difference being the possibility of obtaining also modes with negative modal masses which cannot in general be neglected. Following a similar procedure damped systems can be also studied; if the distribution of damping is not “proportional”, an iterative approach that is cost-effective from the computational point of view is suggested both for natural and non-natural systems. The modal procedure is extended also to the case of non-isotropic machines with a non-symmetric stator or non-symmetric rotor. A computational procedure which can be used for the computation of the “backbone” and the stable branches of the response of non-linear systems is shown. An example of application to a rotor with many degrees of freedom is also reported.

1. INTRODUCTION

The evaluation of the unbalance response is one of the most important steps in the dynamic analysis of rotors. By resorting to models based on the finite element method (FEM), there is no difficulty in computing the displacements, stresses and strains caused by an arbitrary unbalance distribution; however, when the number of degrees of freedom is high, this approach can be quite expensive in terms of computation time and cost. A noticeable reduction of computation time could be achieved by adapting the modal techniques currently used in structural dynamics to the study of rotors. The possibility of employing a modal approach to rotating structures was first suggested before the techniques based on the FEM gained large popularity. Child [l] obtained the equations of motion of complex rotors modeled as a number of rigid bodies that were elastically connected, the deflections of which were defined in a rotating frame of reference. The eigenproblem related to the “natural” part of the system was solved. Modal coordinates were introduced and modes with natural frequencies higher than the maximum spin speed of the system were neglected. Obviously, this procedure contains approximations due to the fact that the influence of the gyroscopic terms is not taken into account in the formulation of the eigenproblem. Child [l] proceeded to obtain the transient response of a turbopump by means of numerical integration of the equations of motion; the model also takes into account viscous and hysteretic damping and non-linear bearings. An interesting method which allows one to evaluate the dynamic behaviour of a gyroscopic system by means of a modal approach was suggested by Meirovitch [2]; in t Now at Sincrotrone

Trieste,

Scientific

Division,

Trieste,

Italy

129 0022-460X/90/

130129 + 25 $03.00/O

@ 1990 Academic

Press Limited

130

G.

GENTA

AND

F. DE

BONA

this case complex co-ordinates are not used and therefore only the mass and the stiffness matrices are symmetric, while the gyroscopic matrix is skew symmetric. With the equations in phase space form, the set of n equations of second order is reduced to a set of 2n equations of the first order, the relevant matrices of which are respectively symmetric and skew symmetric. Obviously, in this case the conventional modal approach is not applicable, and a method to transform the partially skew symmetric complex eigenproblem into two real symmetric eigenproblems was suggested. In the latter case the orthogonality properties of the eigenvectors hold and modal uncoupling can be performed. In this way both the steady state and the transient [3] analysis can be performed. More recently, with the introduction of FEM codes, analysis of complex rotor models with many degrees of freedom has become normal computational practice and several techniques have been suggested in order to reduce the size of the problem. Some authors suggested using dynamic reduction techniques [4,5], while others [6,7] have preferred the “component mode synthesis” approach. It is well known that the latter is particularly suitable when the dynamic behaviour of structures which have different geometrical configurations must be studied (i.e., robot arms); the present authors believe that this method does not hold real advantages in the study of systems the configuration of which is fixed. The use of “component mode synthesis” in rotor dynamics was suggested by Geradin et al. [6]; by following the method and neglecting the “normal modes” the frequencies of which are higher than the maximum spin speed, it is possible to reduce the number of equations of motion. By resorting to phase space a non-symmetric coupled eigenproblem is then obtained. In reference [6] it is stated that in the case of undamped systems with no stiffness coupling between degrees of freedom related to the lateral inflection planes (the xz and yz planes in Figure 2 of section 2.2), if the actual gyroscopic matrix (which is skew symmetric) is replaced by its associated symmetric one, the form of the eigenproblem is very similar to that of a natural system. Following a similar procedure, Glasgow and Nelson [7] applied “component mode synthesis” to study the stability of rotor-bearing systems. The disadvantages linked to the non-symmetric formulation were overcome by using biorthogonality relations, which allowed uncoupling of the equations of motion in the case of skew symmetric systems; nevertheless it must be pointed out that in this case the transposed eigenproblem must also be solved. All of the above-mentioned methods were formulated in order to simplify the evaluation of the response to transient loads such as that of “blade loss” type or for the study of stability of damped systems, but they do not introduce actual advantages in the case of the steady state unbalance response. The modal approach suggested here allows one to achieve a very noticeable reduction of the number of computations needed to obtain the unbalance response, and consequently is particularly suitable to models including a large number of degrees of freedom, the response of which at many spin speeds has to be computed.

2. UNBALANCE 2.1,

RESPONSE

OF UNDAMPED

AXISYMMETRIC

ROTORS

EQUATIONS

The equation describing the dynamic behaviour of an axisymmetric undamped rotor in terms of complex co-ordinates is [8] (a list of symbols is given in the Appendix) [M]{cj}-iw[G]{c.j}+[K]{q}={F}w2ei”’.

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

APPROACH

131

If the steady state unbalance response is required, the solution is (2)

(4) = {qOIe’“‘.

By introducing the solution (2) into equation (I), the following algebraic equation is found:

[-~2([~l-[GI)+~~ll~~,~=~~~~‘.

(3)

The corresponding homogeneous equation can be used in order to obtain the critical speeds. Although equation (1) is the mathematical model of a non-natural system [9], equation (3) can be interpreted as the equation of motion of a natural system with a mass matrix equal to ([Ml - [G]) and consequently can be dealt with by resorting to numerical modal analysis. In the present case, however, the relevant mass matrix ([Ml -[G]), although symmetrical, can be non-positive definite, owing to the presence of the gyroscopic terms. As a consequence, as is well known, not all the eigenvalues w2 which can be obtained from the homogeneous equation (4) are positive and yield a real value of the critical speed. The critical speeds are fewer in number than the natural frequencies of the system at any speed. Owing to the symmetry of all the relevant matrices, the eigenvectors {QI} are still orthogonal with respect to both the stiffness matrix [K] and matrix ([Ml -[G]) and consequently can be used to perform a modal uncoupling. If the modal co-ordinates defined as (40) = are introduced,

[@1{17cJ

(4)

equation (3) can be written in modal form as (-u’[r]+[wf]){~}={F)wz,

(5)

where the eigenvalue matrix [wf] and the modal forces {F} are respectively [wf] = diag {o:, o:, . . . , wf,},

~=[{~,}T{F}I{~i}T([~l-[Gl){~;}l.

(67)

The only difference between equation (5) and the usual modal formulation is that in the present case some of the terms W: can be negative. The excitation is also proportional to w’, as in the case in which a non-rotating structure is excited by motion of the constraints (with the obvious difference that in rotor dynamics o means angular velocity, while when dealing with non-rotating structures it is the angular frequency). The ith equation of motion can be directly solved, yielding the value of the ith modal response as T&F, = l/(+0’-

l),

(8)

which coincides with the unbalance response of a “Jeffcott rotor”. The solution (8) is plotted in Figure 1, for both positive and negative values of the eigenvalue w f . Once the eigenproblem has been solved, the computation of the modal forces and then of the modal responses is straightforward and requires little computational work. The response of the system in terms of “physical” co-ordinates {qo} can be easily obtained through equation (4). As in most cases only a few of the responses qoi are required, it is not necessary to compute the whole vector {qo} for all values of the speed. Also, the power of the modal formulation resides in the well known possibility of obtaining a good approximation of the response by considering only a limited number of modes, thus reducing the size of the problem to be solved, both searching a limited number of eigenvalues and eigenvectors and recombining only few modal responses. In

132

Figure

G. GENTA

1. Modal

response

of the system

for modes

AND

F. DE BONA

with both positive

(-)

and negative

(- - -) eigenvalues.

this case, however, some of the modes with a negative modal mass must also be considered. In fact these modes cannot be neglected even in the case in which the corresponding modal forces are small. It must be stressed that in the procedure for normalizing the eigenvectors which was used above unit values for the modal masses are assumed. In this way the modes with negative eigenvalues are characterized by negative modal stiffnesses. The physical behaviour of the system suggests, however, considering such modes as modes with negative modal mass. Actually, the “negative” modes correspond to imaginary critical speeds, which can be interpreted as “virtual” critical speeds at which no resonance is possible owing to the stabilizing effect of the gyroscopic moments. After all, they are negative because the mass matrix of the system is not positive definite, and the stiffness matrix has nothing to do with it. At one of the actual critical speeds one of the modal responses becomes infinitely large, while at a “virtual” critical speed all responses are limited and that of the modes directly involved takes the value noi/ Fi = 4. Even if no resonance is caused by the “negative” modes, they can have a certain influence on the global behaviour of the system, and a criterion for the selection of the relevant modes can be that of neglecting all the modes characterized by a value [w:] well in excess of the square of the maximum spin speed to be considered. 2.2. EXAMPLE: ROTOR WITH FOUR DEGREES OF FREEDOM In order to show how the various modes interact in an actual case, consider the simple rotor with four degrees of freedom sketched in Figure 2. The equation of motion (3) can be written in the non-dimensional form (9) where the non-dimensional spin speed w* is W* = wJm/(k, + k,). The response of the system is then a function of the non-dimensional parameters (Y,p and 6, of the eccentricity E and of the angle x between the spin axis and the principal axis of inertia of the rotor. For k, = k2 and a = l/4, the non-dimensional parameters (Yand /3 are (Y= O-3125 and p = -0.25. Two cases will be dealt with: 6 = 2 and 6 = -2. In the first case, J, > Jp and consequently the rotor behaves as a natural system. The following eigenvalues and eigenvectors,

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

133

APPROACH

Figure 2. Simple rotor with two complex degrees of freedom. Sketch and geometrical definitions. (.I,-J,,)/ml’, (Y=(k,a2+k,b2)/1’(k,+kz); p =(k,h-k,a)/l(k,+kl).

normalized

in

w

iT=

such a way that the modal masses have unit values, are immediately found:

T’= 1.0355,

[@I= 0.6932 0.1971

UT’= 0.12071,

-0.1394 0.9804 I ’

Both eigenvalues are positive, as expected. The response of the system to a unit static unbalance is shown in Figure 3(a) (translational degree of freedom) and Figure 3(b) (rotational degree of freedom). It can be seen clearly that self-centering is due mainly to the second mode, which involves essentially

4

2

3

1.5

2

1

1 \’ 9

X O5 : *o O -0.5

0 -1

-1

-2

-1.5

-3 -4

I

-2

I

(b)

_

(dl

_

Figure 3. Unbalance response of the system of Figure 2, in the case J, > J,,. (a) Response to static unbalance, translational degree of freedom; (b) response to static unbalance, rotational degree of freedom; (c) response to couple unbalance, translational degree of freedom; (d) response to couple unbalance, rotational degree of freedom. -, Total: ----, mode I; ---, mode 2.

134

G. GENTA

AND

F. DE

BONA

the translational motion. As no couple unbalance was assumed, the rotational motion tends to vanish with increasing speed, owing to the combined effect of the two modes. The response to a couple unbalance (Figures 3(c), (d)) is similar, but the two modes exchange their role: here self-centering is dominated by the first mode, while translational motion tends to vanish owing to the interaction between the modes. For the other case with 6 = -2, i.e., with J,
I.0125 0.1121

[@I= 0.7159

1.

The response to a unit static unbalance is shown in Figure 4(a), (b). The response of the translational degree of freedom z. is completely due to the second mode; the other one, cpo, is strongly influenced by the first mode, which is “negative”. At very high speed, the response vanishes, as the two modes cancel each other, as was predictable from the absence of couple unbalance.

1

-1 -2 -3 -4

ly

2

4-

1.5

3-

1

2-

I

I (d)

-

0.5

$

O -0.5 -1 -1.5 -21

0

I 0.5

1 1

1 1-5

Figure

I 2

4. As Figure

-4

L 0

_ ,0~5

II 1

I 1.5

I z

lb*

3, but with J, CJ,,.

The response to a unit couple unbalance is shown in Figure 4(c), (d). Here the two modes cancel each other at high speed in the response related to the translational degree of freedom, while the rotational response is mainly influenced by the first mode. The tendency to self-centering starts at speeds which are in excess of the “virtual” critical speed, and is only marginally disturbed by the growth of the amplitude near the actual critical speed.

UNBALANCE

3. AN ITERATIVE

RESPONSE

APPROACH

OF

ROTORS:

A MODAL

FOR THE SOLUTION

135

APPROACH

OF DAMPED

PROBLEMS

3.1. ANALYSIS Consider a generic damped system governed by the linear equation of motion

[Ml{~}+[Cl{~}+[KI{x}={F(t)}.

(10)

It is well known that it is not possible to uncouple the equations of motion except in some cases. In order for it to be possible to diagonalize by modal transformation, the damping matrix [C] must verify the Caughey relation [ 10, 11, 121 [aw-‘[~I

= [KI[Mlr’[Cl.

(11)

As this condition is necessary and sufficient, other conditions, such as that of Rayleigh (i.e., the usual formulation of “proportional damping”, [C] = a[ M]+P[K]), shall be considered as a particular case of equation (11). In the present paper damping matrices verifying this condition will be referred to as “pseudo-proportional”. Nevertheless, actual damped systems are hardly ever of either Caughey or Rayleigh type, and therefore these conditions are of little help in practice. The usual way to solve equation (10) is that of working in phase space, following the traditional method of Duncan [13], but this means doubling the dimension of the eigenproblem, which involves, in the case of systems of many degrees of freedom, high computational costs. In a general case, if the damping matrix does not satisfy condition (1 l), the equation of motion is changed by the modal transformation into

~~l~;i~+~~l’~~l~~l~ri~+~~I~77~ = VW,

(12)

where [k], [K ] and {F} are the usual modal mass and stiffness matrices and force vector, while the non-diagonal modal damping matrix [@]‘[c][0]=[C]=[C*]+[C**]

(13)

can be considered as the sum of the diagonal matrix [C*] containing all terms of the diagonal of [C] and by a symmetric matrix [C**] containing all other terms. By applying the inverse modal transformation, the original damping matrix [C] can be considered as constituted by

[cl=rcl,+caI,~

(14)

where [Cl, = ([ @J-‘)‘[ C*][ @I-’ and [Cl,, = ([ @] -I)‘[ C**][ @I-’ are respectively the pseudo-proportional component of the original damping matrix and the component expressing its deviation from “pseudo-proportionality”. If the system is lightly damped, damping introduces only a slight coupling into the system and affects the response of each mode only at frequencies which are near to the relevant natural frequency. In this case the deviation from “pseudo-proportionality” can be neglected, and the system can be treated as proportionally damped. This is well known, but the approach suggested here allows one to obtain in a simple way the diagonal part of the modal damping matrix. It also shows that it is not actually necessary that the system is lightly damped in order to uncouple the modes: it is sufficient that the nonproportional part [Cl, of the damping is small. The fact that the solution obtained is a good approximation of the exact one around the natural frequencies of the undamped system can be easily shown. It is well known [14] that, if the excitation forces are harmonic, are all in phase (the vector {FO} is real) and follow a particular distribution, it is possible to find a set of eigenvectors {r}i which are orthogonal with respect to the matrices [Ml, [K] and [C]

136

C.

GENTA

AND

F. DE

BONA

and therefore uncouple the equations of motion. They are the so-called forced normal modes which depend on the excitation frequency and differ from the eigenvectors of the undamped system {@}i. When the excitation frequency A is equal to the ith eigenvalue A, of the undamped system, {Z}i coincides with {~}i. As a consequence, in this case the ith eigenvector of the undamped system is orthogonal also with respect to the damping matrix; i.e.,

wTc1w1,

= c,.

(15)

This means that for A = A, the ith modal damping obtained from equation (15) coincides with the exact one. At frequencies which are close to the ith natural frequency, the behaviour is strongly influenced just by the ith modal damping, and consequently a very close agreement between the approximated solution and the exact one is to be expected. There are, however, cases in which the very nature of the problem causes heavy coupling between the equations of motion in modal form. Also in this case it is possible to avoid having to resort to Duncan’s phase space approach. Equation (12) can be written in the form

r~l{;i}+r~*1{7j}+[K1177}= -K**l{7i~+m)I.

(16)

In order to compute the response to a harmonic excitation {F(t)}={F,}e’“‘,

(17)

the following algebraic equation can be used: [-A*[&Z]+iA[~*]+[K]]{n,,}= By normalizing the eigenvectors equation (18) simplifies to

-iAIC**]{nO}+{FO}.

(18)

in such a way that the modal masses have unit values,

[-A2[Z]+2iA[A~~~]+[Af]]{~0}~-2iA[Ai~~*]{~o}+{~b},

(19)

~51=~5*l+r5**l=tr~lTr~l~~l~~~lTr~lr~ll-’, m = r@l’~m@l K@lTr~lE@ll-‘.

(20) (21)

where

Equation (19) can be solved in an iterative way, following the Jacobi or Gauss-Seidel iterative technique [ 151. Note that the matrix on the left side is diagonal and consequently the solution at each iteration is very fast. If the Jacobi scheme is followed, the first iteration is obtained by setting the right side {no} = 0; the first approximation solution is then coincident with that of a system with “pseudo-proportional” damping. If the Gauss-Seidel technique is used (i.e., while solving the jth equation all values of ni with i
F:-i(A,/A)

xjr: 2&njk’-i(Aj/A) [(Af/A*)+2i&(A,/A)-II

Cy=i+l 2&~77~k-‘)

(22)

In the case of systems which are only slightly non-proportional, the convergence of the iterative procedure can be easily assessed: a sufficient condition for convergence is that the matrix of coefficients of equation (10) is diagonally dominant [ 151. For highly damped systems with strong non-proportional characteristics, particularly when excited in the vicinity of a resonance, such condition, which at any rate is not necessary, cannot be met.

UNBALANCE

RESPONSE

OF ROTORS:

A MODAL

137

APPROACH

In several tests run by the authors convergence was always achieved, usually with a very small number of iterations. The computation time was always far shorter than that required by the conventional approach. One of the most important features of the modal approach to the computation of the dynamic behaviour of a system is that of allowing one to obtain a very good approximation with only a limited number of modes. The procedure suggested here retains this helpful characteristic but, owing to the coupling between the modes introduced by the nonproportional part of the damping matrix, it can be difficult to choose which modes are to be neglected. In theory, all modes can be strongly involved in the response at any frequency. This does not mean that a very good approximation cannot be obtained with a limited number of modes, but only that it is not possible to give a general rule for the selection of the relevant modes. In highly damped systems, it can be expedient to compute the response at a small number of frequencies with the conventional technique (which is fast if only a few computations are required) and to use the results so obtained to evaluate the accuracy of the modal computation. 3.2. EXAMPLE: THE QUARTER-CAR MODEL An example of a simple system which is usually quite heavily damped and the damping of which is far from being proportional is the so-called “quarter-car” model with two degrees of freedom, used to study the dynamic behaviour of motor vehicle suspensions (see Figure 5(a)). The model, in which masses m, and m, simulate respectively the sprung and unsprung masses of the vehicle, is excited by the motion of the supporting point with a known law h(t), in order to simulate the motion on uneven ground. 2.5

/

1

I

2.5

I

0

2

4

6

8

Figure 5. Sketch of the “quarter-car” model “sprung mass”; (b) response of the “unsprung (- - -) after different numbers of iterations.

,

I

,

I (b)

(0)

10 and its frequency mass” computed

response. (a) Sketch and response of the directly (-), and iterative computations

The equation for the response of the system to a harmonic excitation with frequency A and amplitude h,, is [

-**2[; ;]+[_;l;;]+2i@[ _: -;]](:;]={;},

(23)

where u = m,/ m,,, p = k,/ k, l= c/2fi, and the non-dimensional frequency and amplitudes are defined as h * = A/m, x$, = x,,/ ho and xZO= x,,,/ h,, . The following parameter

Ci. CENTA

138

AND

F. DE

BONA

values are assumed: p = 0.1; p = 4; 5 = 0.43. The value of { chosen is close to that which optimizes riding comfort, at least under widely accepted assumptions: i.e., < = -/2fi. The amplitude of the motion of the two masses can be obtained by solving directly equation (23), to obtain the curves shown in Figure 5. Note that, owing to the fairly high value of the damping, the response shows only one “resonance peak”. The eigenvalues, the eigenvectors (normalized in such a way that the modal masses have unit value) of the undamped system and the modal forces are respectively AT”= 0.796761, A;* = 50.203239, 0.997941

[@I= 0.202820 1

F,=

0.064137 -3.155767

1’

0.811281 and F2 = - 12.623067. The modal damping yielding

matrix can be computed,

By separating the proportional and non-proportional parts of the modal damping matrix and by applying the inverse modal transformation, the following damping matrices are readily obtained:

[Cl, =

1,

[Cl,=[

0.2819 -0.6905

-0.6905 -0.0282

1’

It is easy to verify that not only does matrix [Cl, satisfy condition (1 l), but in this case it is actually proportional (which is obvious, as the system has only two degrees of freedom), being a linear combination of the mass and stiffness matrices with coefficients 0.4087 and O-1695 respectively. As the damping in this case is neither small nor proportional, a very rough approximation can be expected from neglecting the elements of the modal damping matrix lying outside the principal diagonal. The iterative procedure suggested was consequently applied. The results obtained by using a “Jacobi” approach are plotted in Figure 5, after having been transformed back from modal to “physical” co-ordinates. The convergence, which is already quite fast, can be further hastened by resorting to the Gauss-Seidel iterative scheme. The example shows that even in a case in which the damping is very far from being small or proportional, a small number of iterations are required. In the case of lightly damped systems, the iteration “of order zero” (i.e., the results obtained with the modal damping matrix assumed to be diagonal) yields already very good results.

4. COMPUTATION

OF THE

RESPONSE

OF

DAMPED

ROTORS

The equation of motion of a linear rotor with viscous damping associated to both the rotating and non-rotating elements is

[~1~ij~-~~~~1~~~+~~~~1+~~~1~~~~+~~~1-~~~~~1~~q~=~~~~~~~‘“’+~~~~. (24) For the computation of the unbalance response, the solution (2) can be directly introduced into equation (24), yielding

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

139

APPROACH

As in the case of undamped systems, equation (25) can be regarded as the equation of motion of a natural system with mass matrix ([Ml - [G]), excited by a harmonic forcing function with frequency w and amplitude proportional to w’. The computational procedure described in the preceding section for damped natural systems can be used. As the non-rotating damping matrix is not, in general, pseudoproportional, the iterative approach may be needed. The damped response of the modes with positive modal mass is the usual response of the damped “Jeffcott rotor”, with the characteristic feature of having a peak displaced towards the high frequencies, if compared with the usual harmonic response of nonrotating damped systems (due to the fact that the excitation is proportional to w’) (see Figure 6(a)).

a B

-60

x -120

-60

(d)

-180

-90 0

1

2

3 w2/lu(

0

I

L

1

2

Figure6. Modalresponseof dampedsystems.Amplitude(a) and phase(cl of amplitude

(b) and phase

(d) of a mode with negative

3

I a mode with positive

eigenvalue;

eigenvalue.

The response of the modes with negative modal mass is shown in Figure 6(b), for different values of the modal damping. In this case, if w = m, the absolute value of the modal response and the phase are respectively voi/Fi = 1/2m and cp= arctg ( -l&O. Th e amplitude at a “virtual” critical speed is then always smaller than O-5 and the phase is smaller than 45”, if modal damping is smaller than the critical damping. 5. MACHINES WITH NON-SYMMETRIC STATOR Consider a rotor running on supports, the characteristics of which are non-isotropic in the rotation plane. The unbalance response is more complicated than that of a rotor belonging to a machine in which all parts are axially symmetric, mainly because a backward-whirl motion with frequency w is superimposed on the forward-whirl characteristic of the unbalance response. As the possibility of modal uncoupling is strictly linked to the fact that all the relevant matrices are symmetric, a formulation based on the use of complex co-ordinates is required [ 81.

140 The

G. GENTA

AND

F. DE

BONA

equation allowing one to compute the amplitude of the unbalance response, {q}={q,}e’“‘+(q,}e~““‘,

(26)

is equation (58) in reference [8]. It can be written in the form

-~2[~l~qo~+i~[Sl~qo~+[~l~qo~=~~:’~~2,

(27)

where the relevant matrices and vectors, all of order 2n, are

(28) Here a double bar indicates the conjugate of a complex matrix or vector. All submatrices included in equation (27) are symmetric. Consequently the condition for which a modal uncoupling of the undamped problem can be performed is that the mass and stiffness matrices, both in their mean and deviatoric part, are real. For the mass matrix, they are defined as

~~l,=t~~M,l+~M,~l~+~t~~~,,l-~~x,l~, [ML =~([M,l-[M,,I)+i~([M,,.~l+[M,,I).

(29)

The condition for real mass matrices is that axes x and y coincide with the principal axes of inertia of all elements. Similarly, for the stiffness matrix,

[Kim =~~~~.~l+~~~l~+~S~~~~,l-~~~,,l~, [Kid =f([K,l-[K,,l)+il([K,,,l+[K,.l).

(30)

Here the conditions for real stiffness matrices are two: the axes x and y coincide with the principal axes of elasticity of all elements, and all stiffness matrices, written with reference to axes x and y, are symmetrical. This obviously excludes stiffness matrices deriving from the linearization of the behaviour of hydrodynamic bearings. Consider as an example the rotor of Figure 2, in which the supports have different stiffness in the x-z and y-z planes. Upon assuming k, = k+Ak and kV= k - Ak, and introducing the non-dimensional spin speed w* defined in equation (9), the relevant matrices and vectors to be introduced into equation (27) are

where 6’= (JP + J,)/ml’ and y = Ak/ k. The amplitude of the response to a unit static unbalance, computed with the same numerical values of the parameters as for Figure 4 and S’= 8, y = O-2, are plotted in Figure 7(a). In the figure q,( 1) and q2( 1) are respectively the amplitudes of the forwardand backward-whirl components of the orbit of the centre of gravity of the rotor. They are shown in Figure 7(b), for some selected values of the spin speed, in order to obtain the orbits of the centre of gravity. In Figure 7(a) q,(2) and q2(2) are respectively the amplitudes of the forward and backward components of (p&l). The eigenvalues are wf’ = 0.0300, w$* = -0.1275, wf2 = 0.7919 and wq*’- 1.1885. Only one of the eigenvalues is negative and consequently three critical speeds exist. However, from the response

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

141

APPROACH

6-

-6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

X/E

Figure 7. Unbalance response of the system of Figure 2 but with non-symmetric backward components of the response as functions of the spin speed; --, q,(l); - - -, q*(2); (b) orbits of the centre of gravity at some selected values of speed.

supports.

----,

la) Forward

ql(l ); -. -,

and q,(2);

plotted in Figure 7(a) it is clear that the first one influences the behaviour of the system only to a very small extent. If matrices [R] and [T] are not real or non-symmetrical and if the system is damped, it is possible to proceed in the same way as seen for damped axisymmetrical rotors. The following four steps allow one to solve the problem. (i) The real parts of the mass and stiffness matrices [R] and [T] are made symmetrical just by adding their transpose and dividing by 2. The related eigenproblem is solved. A set of eigenvectors which can be used to define the modal co-ordinates and play the same role as the modes of the undamped system for natural systems is so obtained. (ii) The modal mass, stiffness and damping matrices and the modal load vector are computed, starting from the original complex, unsymmetrical matrices. All modal matrices can consequently be complex and non-diagonal. (iii) The complex, diagonal part of the matrices is extracted and the non-diagonal part is carried to the right side of the equation. (iv) The solution of the equations is sought iteratively and the results are recombined in order to compute the response in “physical” co-ordinates. Generally speaking, the convergence of the iterative procedure cannot be guaranteed. It is obviously faster if the deviations from symmetry and the imaginary parts of all matrices are small, and if damping is either small or not too far from being proportional. It must be stressed that all the above-mentioned effects act in such a way as to couple the modal responses; it thus can become difficult to assess which modes must be retained and which can be dropped. 6. MACHINES

WITH

NON-SYMMETRIC

ROTOR

If the rotor itself is non-isotropic in the plane of rotation but the supports are axially symmetric the response to unbalance is simpler: the whirl motion is again circular, and a solution of the type of equation (2) can be assumed. In this case it is expedient to resort to a formulation based on the matrices descr.ibing the behaviour of the system in the reference frame which rotates at the spin speed.

142

G.

GENTA

AND

F. DE

BONA

Equation (67) in reference [8] yields, for the case of unbalance response, [-w2[R]+w[S]+[T]]{x,,}={F:‘}w’,

(31)

where the relevant matrices and vectors, all of order 2n, are

Again all submatrices are symmetric. Consequently, the condition for which modal uncoupling of the undamped problem can be performed is that the mass and stiffness matrices with xy and yx subscripts are equal. Generally speaking, this is the case for non-symmetric rotor elements and, consequently, the undamped problem reduces to a form which is equal to that seen for the axisymmetric machine, the only difference being that the size of all matrices is doubled. Also, in the case of damped systems there is not much difference from the case of axisymmetric systems, the only difference being that a real equation of motion is obtained. A case which is particularly simple is a rotor in which no inertial or elastic coupling exists between the x-z and y-z planes. The free behaviour of the system can be studied by solving the two separate eigenproblems [-w2([Ml,

- [GI) + [Kl,l

[-~‘([W~-[G1)+[~1,1

Re ({xJ)

= 0,

Im({xJ)=O.

(33)

With subscripts 1 and 2 respectively to denote their eigenvalues and eigenvectors, the following equations of motion in modal co-ordinates are obtained for the damped system: WI*

- w2[U

[ ~[@l:[cn1r@l2

-;y$;plq

pi;}

= d{

$}.

(34)

If damping is proportional (which is most unlikely, as matrix [C,] should satisfy condition (11) for both eigenproblems) or is small enough to neglect the elements outside the main diagonal in the modal damping matrices, each pair of equations (34) related to the two degrees of freedom noil and Toi is uncoupled from all the other equations. If this does not happen, an iterative approach corresponding to that seen in section 3 can be used. Note that if the characteristics in the x-z and y-z planes are equal, equation (31) can be written in the complex form (25).

7. AXISYMMETRIC 7.1.

EQUATIONS

OF

NON-LINEAR

ROTORS

MOTION

equation of motion of an axisymmetric non-linear rotor in which the non-linearities are involved only in the “elastic” part of the system and damping can be assumed to be linear, can be written [16] by adding a term {f(qi)} expressing the non-linear behaviour of the system to equation (24): The

=

{Fn}+{Fr}02 e’“‘.

(35)

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

APPROACH

143

It is well known that it is not possible to separate the study of the “free” behaviour of the system from that of the response to unbalance or the deformation under static loads. An approximate solution of equation (35) can be of the type

{9~={9J+{q0~eiw’+i {qkleiAL’. Such a solution is an approximation, as for each response it takes into account only the fundamental harmonic of a time history which is certainly non-harmonic, owing to the non-linear nature of the problem. However, at least in the case of a rotor with four degrees of freedom with “hardening” non-linearities, it has been shown that the self-excited vibrations start to interact heavily with the unbalance response only at speeds in excess of the threshold of instability of the linearized system [17]. As rotors are usually designed in order to work at speeds not exceeding such a threshold, and as most non-linear elements, such as rolling elements bearings, show a hardening behaviour, the unbalance response will be studied independently. Static forces will also be neglected, which is obviously an approximation that in some cases can limit the application of the present approach. A solution of the type of equation (2) is consequently assumed as {91= {qO]eiw’.

12a)

Owing to the axial symmetry of the system, equation (2) certainly expresses an exact solution of equation (35), even if this solution is not the only one possible and is not necessarily stable. If the considerations in reference [17] are intuitively generalized to systems with many degrees of freedom, such a solution can be considered to be the governing one, at least in the case of hardening systems, up to the threshold of instability of the linearized rotor. Equation (2a) allows one to transform the non-linear differential equation (35) into the non-linear algebraic equation (37) Equation (37) can be solved directly by using an iterative technique, which is, generally speaking, very difficult. As often the number of degrees of freedom directly involved in the non-linearities is small, condensation techniques can be used to simplify the solution [ 16,181. If the number of degrees of freedom which remain after condensation is greater than two, however, the iterative solution is a difficult task. Moreover, the condensation has to be performed for each value of the speed, and this leads to long and costly computations. The extension of the modal approach seen in section 2 for linear rotors can simplify the solution. 7.2.

SOLUTION

VIA

MODAL

CO-ORDINATES

The modal co-ordinates of the undamped linearized system defined by equation (4) can be used in order to express the solution of the non-linear equation (37). The modal transformation yields

(38) The equation for the undamped

system can be written in the form

144

G.

GENTA

AND

F. DE

If damping is not neglected, it is possible to extract modal damping matrix as in equation (13), obtaining

RONA

the diagonal

part

[c*]

from the

~-~‘[~l+~~l+~~[C~l~~~0~=-~~~~~*1~~~,,~-[~lT~f~q,i~~+{F,,~w’.

(40)

Equation (39) or equation (40) cannot be solved by using a Jacobi or Gauss-Seidel scheme. Even in the case of systems with a single degree of freedom, or unlikely systems in which the non-linearities are such that each equation (39) contains only one unknown Toi (this situation could be referred to as “proportional non-linearities”), it is possible to show that no convergence is obtained, at least not in a speed range starting at the critical speed of the relevant mode and extending to the speed at which multiple solutions are found [ 181. If the approximated solution {qj,“} is known, it is possible to write the following series qo,, for functions {f(qoi)} for limited variations of parameters

{f(qo;)}={f(46:))}+rSl({qo>-{qb”}), or, in terms

of modal

(41)

co-ordinates,

(42)

{.f(‘loi)}={f(46:‘)}+[SI[~1({770}-{77bl’}). The equation which allows one to compute from the values at the ith iteration is

the displacements

at the (i+

l)th

iteration

(-w2[A2]+[K]+iw[C~]+[0]T[S][@]){~~+”} =[~]TIS]{q~‘}-[~]T{f(q~‘)}-iw[~~*]{~~)}+{Fo}w’.

(43)

The computation of the Jacobian matrix [S] is extremely complicated if the complex co-ordinates are maintained. It is consequently expedient to separate the real part from the imaginary part of equation (43), which yields [K]-w2[M]

4Cl

(44) where subscripts x and y indicate, respectively, the real (x-z plane) and imaginary plane) parts of the complex co-ordinates and forces. The matrix [S’] is linked with the Jacobian of the functions f by the equation

(y-z

where {x} and {y} are the real and imaginary parts of {qo}. The solution of equation (45) can be made simpler if the elements of the four relevant submatrices in [S] which lie outside the main diagonal are neglected. This obviously does not affect the solution of equation (44), but only the number of iterations needed to reach convergence. Moreover, in some cases convergence cannot be obtained. Several tests performed by the authors did show good convergence in most cases; the simplification which can be obtained in this way is at any rate so effective that it is certainly worthwhile to investigate whether or not this procedure is applicable to the particular cases under study.

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

145

APPROACH

The solution of equation (44) reduces to

&+“=[ai(Iti

-w2tii+

Dir,.)+ b,(d:-

D,,.,)]/R

(46)

77j~+‘)=[bi(Ki-w2n),+D,~i)-a;(W~~+D,,?,)]/R

where

R=(K,-o2~ii,+Di,,)(~,-w’n;ii+D,,.,.)-(uC”+D,,,.)(-“C”+D

,,,,),

If no damping is present, a complete uncoupling between equations (44) is reached and the solution simplifies to n:,;+‘) = a,/( K, - w2ri?; + D,,,) n::,+“= b;/(K,-w’n;i,+ D ,,,.) ’

(47)

where

If the unbalance is contained in a single plane, i.e., one of the two vectors {FX}or {F,.} vanishes, then the deformed shape is contained in the same plane and one of the two vectors (77,) or {n,,} vanishes. The number of equations to be solved iteratively is halved. The computation of the response can be performed by using the following scheme: (i) computation of matrices [K], [Ml, [G] and [C,] and vector {F,,}; (ii) application of a static condensation procedure, in order to reduce the number of degrees of freedom, if needed (it must be remembered that this introduces approximations in the results); (iii) solution of the eigenproblem related to the linearized undamped system and computation of the modal matrices; (iv) iterative computation of the response curves, starting at a speed which is well below the first critical speed of the linearized system; in this way the starting solution for the first computation can be the undeformed shape and convergence is granted. Each iteration is very fast, as the equations are uncoupled and the solution is in closed form. At each iteration, however, a transformation from modal to physical co-ordinates is needed, as the non-linear part is function of the latter. After convergence, the solution can be computed in terms of physical co-ordinates and then expanded back in order to obtain the values of the co-ordinates which were condensed in step (ii). The solution used to start the computation at the following speed can be that obtained for the previous one. If the new value of the speed is not too far from the previous one, convergence is usually fast. Obviously, only one solution is found also in cases in which multiple solutions do exist. If damping is low enough to prevent the jump from taking place, the amplitude grows indefinitely after the first critical speed. The results obtained in the way suggested correspond to those experimentally achievable in a test performed with slowly increasing speed.

G.

146

GENTA

AND

F. DE

BONA

If self-centering takes place, it is sufficient to repeat the computation in the downward direction starting from the maximum speed, in order to find other solutions, when they exist. If damping is too low and no jump takes place, the computation can be started from the maximum speed, with the undeformed configuration used as the starting solution. There is little chance of obtaining the unstable branches of the response curve, as the iterative procedure usually does not converge on them. These branches have a certain practical interest, even if the actual system cannot be made to work on them, as unstable solutions are located on the borders of the attraction domains of the stable solutions. 7.3.

COMPUTATION

OF

THE

“BACKBONE”

OF

THE

RESPONSE

In the study of the dynamic behaviour of all non-linear systems the knowledge of the so-called “backbone” of the response is very important. It could be said that the backbone plays, for non-linear systems, the same role which is characteristic of the natural frequencies in the case of linear systems. In the case of a Jeffcott rotor (two real degrees of freedom), the possibility for the “jump” phenomenon to actually take place and consequently to achieve self-centering in the high supercritical range can be assessed simply by searching the intersections between the backbone and the “limit envelope” [19]. Unfortunately, in the case of multi-degree-of-freedom systems the limit envelope can be determined only in a few special cases. The backbone of the response can be obtained by solving the equations of motion of the undamped system after having cancelled the forcing term {F,}:

Using the modal approach for the solution of equation (48) yields (49) The trivia1 solution {no} = 0 exists for all values of the spin speed. If the system has hardening characteristics, then after starting from the first critical speed of the linearized rotor a non-trivial solution, the backbone, starts being present. If the characteristic of the non-linear element is not hardening in the whole field, the backbone can be found also at speeds which are lower than the first critical speed. If difficulties are experienced in order to shift from the trivial to the non-trivia1 solution, a very small unbalance can be introduced. If it is small enough, it does not affect the results, apart from making the procedure shift to the desired solution. The computation of the backbone requires only a small fraction of the time required for the computation of the response and gives useful information on the behaviour of the system. 7.4. A NON-LINEAR SPRING ELEMENT The expressions for the relevant matrices and vectors obtained for a spring element in which the force is expressed as a function of the complex displacement z by the law

f(z) = 41 +&l*)z,

(50)

are to be presented in this section. This law is particularly well suited for modelling some rolling element bearings, in particular preloaded angular contact ball bearings [16], but has a more general application. It can be regarded as the general characteristic of any

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

APPROACH

147

non-linear element with an “odd” behaviour, when all terms of the series but the first two are neglected. It constitutes a sort of second order approximation, the linearized equations being the first order one. Equation (50) has been widely used in non-linear dynamics, starting from the work of Duffing. The linear part of equation (50) can be introduced directly into the stiffness matrix. If the non-linear spring connects the ith degree of freedom to the jth, it follows that

J;(qo)

I

{ “US”)

=

1 I (4i-$1

kP(qi-qj)2

(51)

’

(9,-q,)

or, in real co-ordinates,

_M40) {

txt -xj) cxi_xiJ

1 I

I

=kcL[(xi-xj)2+(Y~-Y;)21

f;x(c?o>

(Yi -y,) r (Yj-Yt) I ’

f;?(%) =k~[(x,-xj)2+(Y,-Y,)21 1 f;;J%J I The Jacobian matrices are [ afx(qo)

1

9

(52)

,afffij =~f~~~~~~i~~)~ -:],

[ _ai

2

9

I

I

1

[

(40)1 = k/L[(Xi-Xj)2+3(yi-yi)2]

aY

[

_,1

-11 . 1

(53)

Obviously the relevant degrees of freedom can be rotational as well as translational. The dimensions of k and CLare defined accordingly. The matrices related to the whole structure are assembled from those computed for the single elements by using the same procedure as in the finite element technique. 7.5.

EXAMPLE:

ROTOR

WITH

TWO

COMPLEX

DEGREES

OF

FREEDOM

Consider, as an example, the rotor studied in section 2. The two supports have now a non-linear behaviour, which can be modelled by using equation (50). The linearized system is the same as already studied. In order to simplify the statement of the problem, it is expedient to choose as generalized co-ordinates the displacements of the points in which the non-linear elements are connected. With symbols z, and z2 designating these complex co-ordinates, the equation of motion reduces to

where the non-dimensional spin speed w* has been defined in equation (9). The results obtained for the undamped system are shown in Figure S(a). The amplitude of the orbit of the centre of gravity and the angle rp, are shown in non-dimensional form for a value of the non-dimensional non-linearity coefficient given by E*P = O-1. As expected, there is no jump in the computation performed with increasing values of the speed, while an abrupt increase in amplitude is encountered when the speed is reduced. Convergence is always very fast, with usually only two or three iterations needed for each value of the speed, with the exception of the speed at which the jump is experienced. In the latter case a larger number of iterations, more than 20, are required. For the study of the damped system a linear viscous damper is added to each support, to obtain consequently a proportional damping matrix. A value of relative modal damping of 0.2 was assumed for both modes.

148

G. GENTA

AND

F. DE

BONA

Figure 8. Unbalance response of the system of Figure 2, but with non-linear supports. (a) Radius of the orbit of the centre of gravity and angle &, for the undamped system; (b) as (a), but for the damped system. . .. linearized w,. . -, Iz&E; - - - -, /lq,$~; - - -, backbone;

The results are plotted in non-dimensional form in Figure 8(b). As the damping is high enough, self-centering takes place and the solutions for increasing and decreasing speed are completely superimposed (no jump occurs). 8. NON-SYMMETRIC

NON-LINEAR

ROTORS

In case of non-linear rotors, no exact solution can be easily found for the case of non-symmetric machines. The very fact that no circular whirling is possible renders a simple harmonic response impossible and causes the onset of subharmonics and superharmonies of the fundamental frequency, as in the general case of vibrating non-linear systems. Approximate solutions for the fundamental harmonic can be obtained by resorting to the usual methods, such as harmonic balance or the Ritz averaging technique (see, e.g., reference [9,20]). Once the relevant non-linear differential equations have been transformed into algebraic non-linear equations, the present modal technique can be applied, following the procedure seen for linear systems. The non-linear terms are added to the modal equations obtained in sections 5 and 6, and an iterative technique of the type seen in section 7 can be applied. 9. EXAMPLE

OF

AN

APPLICATION

9.1. MODEL In order to show an application of the computational procedure suggested here, a rotor with many degrees of freedom has been studied. A good test case is a rotor which has been used by several authors for FEM rotor dynamics computations [6,21,22] and is sometimes referred to in the literature as the “Nelson-McVaugh” rotor. The system was modelled here by using 17 beam elements, one mass element and two damped elastic supports, with a total number of 36 complex degrees of freedom. The DYNROT code, described in detail in reference [23], was used for the construction of the relevant matrices, the solution of the eigenproblem and the computation of the unbalance responses with the conventional approach. Two different choices were made for modelling the bearings. First a linearized computation was performed with the following values assumed for the stiffness and the damping: k = 4.38 x 10’ N/m; c = 2630 Ns/m. A non-linear computation was then attempted with

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

149

APPROACH

use of the model expressed by equation (50) with the following values of the parameters: k = 4.38 x 10’ N/m; p = 1.0 x lOI N/m’. The force-displacement curve so obtained is similar to the characteristic of preloaded angular contact ball bearings. A sketch of the model used is shown in Figure 9.

Figure 9. Sketch of rotor used for the application

9.2.

LINEAR

example.

E = 2,078 x IO” N/m’;

p = 7806 kg/m3;

I = 355

mm.

SOLUTION

Both the damped and undamped responses were computed by using the modal procedure here described. The results obtained for the orbits at nodes 5, 11 and 14 (the centre of the disc and bearing locations) when the mass element in node 5 has a unit static unbalance, are plotted in Figure 10. A total of six modes, one of which has a negative eigenvalue, were used. The results obtained here are very close to those obtained in reference [22], and also shown in Figure 10, even if the logarithmic scale used tends to magnify the errors in the zone in which the response is small. The reduction of the computation time is very noticeable as shown in Table 1, particularly if the computation of the critical speeds has to be performed. In this case the computation time for the damped system has been reduced to l/10 of that of the conventional method. This result, although depending on the hardware and software used, can be considered as typical. 9.3. NON-LINEAR SOLUTION The static unbalance at node 5 used for the computation of the non-linear response is rn~ = 3 x 10S5 kgm. The response is shown in Figure 11. As only two degrees of freedom are directly involved in non-linearities, it is possible to compute the unbalance response by using the procedure developed in references [ 16,181. The results so obtained are also shown in Figure 11. While the iterative modal solution allows one to compute only the stable branches of the response curve, the non-modal solution also yields the unstable branches, which are shown in the figure. As expected for undamped hardening systems, no “downward” jump is found.

150

G. GENTA

/ 2000

AND

F. DE

I

BONA

I

4000

6000

I

6000

10 000

w bad/s)

Figure 10. Response of the rotor of Figure 9, damped linear case. Node 5: -, 11: - - -, modal; 0, non-modal; node 14: -. .-, modal; 0 non-modal.

TABLE

modal;

?,?non-modal;

node

1

Computation time (in seconds) for the modal and conventional solution in the case of the rotor of Figure 9 (damped and undamped cases); hardware HP 9836; code written in HP BASIC

(3) (1) Conventional Undamped Damped Jacobi Damped GaussSeidel

(2) Eigenproblem

Modal solution

(2)+(3)

[(2)+(3)1/(l) (3)/(l)

735 3255

148 148

194 340

332 488

45% 15%

3255

148

312

460

14%

26% 10% 9.6%

Convergence was always quite fast except for values of the speed in excess of the linear critical speed (for which the backbone of the response already exists) but lower than that at which a “self-centred” solution can be found. In some cases in those fields violent oscillations of the results were encountered, and it was found necessary to introduce a “mathematical damping” in the convergence process in order to allow the solution to convergence smoothly. The authors ascribe this behaviour to the intrinsic difficulties of Newton-Raphson methods to avoid searching for a non-existent solution in a field in which the sign of the derivatives makes it possible to suspect the existence of a solution (see Figures l-10 in reference [15]). In at least one case, the authors did find that under the critical speed the relevant function y = f (x) of the equation f(x) = 0 yielding the solution is monotonic and has only one solution; when three solutions exist it has a maximum and a minimum crossing the axis y = 0, three times, while in the field at which difficulties are encountered it again has a wavy pattern, but with only one crossing of the x-axis. The difficulty in applying the Newton-Raphson method in this case is obvious. A possible cure is to introduce a large correction of the solution when strong oscillations between two constant values are

UNBALANCE

RESPONSE

OF

ROTORS:

A MODAL

APPROACH

151

300 250 200 150 100 50 n

105 607

60-

,r 40-

60

c

4P

t

I

!A

60

w (rad /s)

11. Response of the rotor of Figure 9, undamped non-linear case. (a) Orbit at node 5 (centre of disc); node 11 (front bearing); (c) orbit at node 14 (rear bearing). ---, Stable (modal); ?? ma, stable (non-modal); A A A, unstable (non-modal); - - -, backbone. Figure

(b) orbit at

encountered, in order to exit from the troublesome field. Further tests aimed at solving this problem are planned. All the 18 modes of the condensed model were retained. Some tests did show that the strong modal coupling introduced by non-linearities lowers the precision of the solution, if some modes are neglected. Also in this case, the saving of computation time over the non-modal solution is quite noticeable, even if a quantitative statement cannot be given, as no standard procedures which can be taken as a reference exist for the non-modal computation. 10. CONCLUSIONS The use of modal techniques in the analysis of the dynamic behaviour of structures is very useful as it allows one to uncouple the equations of motion and to reduce the size of the problem. The approach proposed here allows one to extend those advantages to the computation of the unbalance response of axisymmetric rotors even when the gyroscopic effects are

152

G. GENTA

AND

F. DE: BONA

not neglected and consequently the system is “non-natural”. Consequently, the unbalance response becomes an easily obtainable “by-product” of the computation of the critical speeds of the undamped system. In the case of axisymmetric undamped linear rotors, as well as in some non-axisymmetric rotors, the modes are completely uncoupled and a further simplification of the problem comes from the possibility of taking into account only a limited number of modes. In the case of damped systems with strongly non-proportional features, the response can be obtained through a fast-converging iterative technique. Obviously, the same iterative approach can be extended to any “natural”, highly damped system. If the rotor has a non-linear behaviour, the steady state circular whirling due to unbalance can be computed by using a modal procedure, based on the modes of the linearized undamped system. Also in this case, the modes are not uncoupled, owing to both the non-linear terms and damping, if present. An iterative procedure based on the Newton-Raphson approach has been developed for solving both the damped and undamped cases. A simplified form of Jacobian matrix can be used in order to obtain, at each iteration, uncoupled linear equations. It has been shown that this simplification does not influence the results of the computation but only the number of iterations required. Owing to the coupling of the modes, the number of modes required to achieve the required precision is generally higher than in the case of the linear system, and it is not possible to predict simply which modes can be neglected. Although this limits the advantages of the modal formulation, the reduction of computation time is nevertheless noticeable. The possibility of resorting to dynamic reduction procedures can, in this case, be very effective in reducing the size of the problem. Some simple examples, and an example of an application, demonstrated the good results and computation time savings which can be obtained by resorting to the modal co-ordinates in the various cases. Although it was not possible to demonstrate the convergence of the iterative procedure in all cases, the examples did show that the convergence is indeed very fast, in most cases just two or three iterations being required. In the case of non-linear rotors some spin speed fields in which convergence problems can occur have, however, been identified. Further work is required in order to understand this issue better.

REFERENCES 1. D. W. CHILD 1974 Journal of Engineering for Industry (May), 659-669. A rotor fixed modal simulation model for flexible rotating equipment. 2. L. MEIROVITCH 1974 American Institute of Aeronautics and Astronautics Journal (October), 1337-1342. A new method of solution of the eigenvalue problem for gyroscopic systems. 3. L. MEIROVITCH 1975 Journal of Applied Mechanics (June), 446-450. A modal analysis for the response of linear gyroscopic systems. 4. K. E. ROUCH and J. S. KAO 1980 Journal of Mechanical Design 102, 360-368. Dynamic reduction in rotor dynamics by the finite element method. 5. G. GENTA 1985 Meccanica 20, 235-248. Consistent matrices in rotor dynamics. 6. M. GERADIN and N. KILL 1984 Engineering Computing 1, 52-64. A new approach to finite element modelling of flexible rotors. 7. D. A. GLASGOW and H. D. NELSON 1980 Journal of Mechanical Design (April), 352-359. Stability analysis of rotor bearing systems using component mode synthesis. 8. G. GENTA 1988 Journal ofSound and Vibration 124, 27-53. Whirling of unsymmetrical rotors: a finite element approach based on complex coordinates. 9. L. MEIROVITCH 1970 Methods ofAnalytical Dynamics. New York: McGraw-Hill. 10. T. K. CAUGHEY 1960 Journal ofApplied Mechanics (June), 269-271. Classical normal modes in damped linear dynamic systems.

UNBALANCE

RESPONSE

OF ROTORS: A MODAL

APPROACH

153

11. T. K. CAUGHEY and M. E. J. O'KELLY 1965 Journal of Applied Mechanics (September), 583-588. Classical normal modes in damped linear dynamic systems. 12. I. FAWZY 1977 Journal of Applied Mechanics (March), 132-134. A theorem on the free vibration of damped systems. 13. R. A. FRAZER, W. J. DUNCAN and A. R. COLLAR 1965 Elementary Matrices. Cambridge University Press. 14. K. ZAVERI and M. PHIL 1985 Modal Analysis of Large Structures. Naerum, Denmark: Briiel & Kjaer. 15. G. F. CURTIS 1978 Applied Numerical Analysis. Reading, Massachusetts: Addison-Wesley. 16. G. GENTA and A. REPACI 1987 11th Biennial ASME Design Engineering Conference on Vibration and Noise, Boston, 441-448. Circular whirling and unbalance response of nonlinear rotors. 17. A. TONDL 1976 Some Problems of SelfExcited Vibrations qf Rotors. Becovice: National Research Institute for Machine Design. 18. G. GENTA and C. SCHIPS 1986 VIII Conuegno Nazionale AIMETA, Torino, 239-252. Studio static0 e dinamico di rotori su cuscinetti a caratteristica nonlineare. University Press. 19. G. SCHMIDT and A. TONDL 1986 Non-linear Vibrations. Cambridge 20. S. TIMOSHENKO, D. H. YOUNG and W. WEAVER 1974 Vibration Problems in Engineering, New York: John Wiley. 21. H. D. NELSON and J. M. MCVAUGH 1976 Journal ofEngineering.for Industry (May), 593-600. The dynamic of rotor-bearing systems using finite elements. 22. G. GENTA 1985 Current Advances in Mechanical Design and Production, Third Cairo University MDP Conference, Cairo, 115-122. A consistent matrix approach for finite element method dynamic modeling of rotors. 23. G. GENTA, A. GUGLIOTTA and G. BRUSSINO 1984 VII Congress0 Nazionale AIMETA. 239-251. Un modello per lo studio dinamico dei rotori. APPENDIX:

SYMBOLS

i

imaginary unit (i = J-i) mass vector of the (complex) co-ordinates ;nql vector of the (real) co-ordinates ;?I matrix [A] in modal form conjugate of complex matrix [A] [Al damping matrix excitation vector ::; gyroscopic matrix [Gl identity matrix 111 lm( ) imaginary part polar moment of inertia J/J transverse moment of inertia Jr stiffness matrix [Kl mass matrix [Ml Re ( ) real part Jacobian matrix rs1 eccentricity relative damping ; modal co-ordinates la} A frequency ith natural frequency A, angle between the spin axis and the principal X w spin speed ? ith eigenvalue ith eigenvector ?i} matrix of the eigenvectors [@I Subscripts d deviatoric nr mean n non-rotating r rotating

axis of inertia