On the computation of stationary deterministic behaviour of
non-linear dynamic systems with applications to rotor-bearing
structures
Citation for published version (APA):
Crooijmans, M. T. M. (1987). On the computation of stationary deterministic behaviour of non-linear dynamic systems with applications to rotor-bearing structures. Technische Universiteit Eindhoven.
https://doi.org/10.6100/IR274906
DOI:
10.6100/IR274906
Document status and date: Published: 01/01/1987 Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at: openaccess@tue.nl
providing details and we will investigate your claim.
DETERMINISTIC BEHAVIOUR
OF NON-LINEAR DYNAMIC SYSTEMS
WITH APPLICATIONS TO
DETERMINISTIC BEHAVIOUR
OF NON-LINEAR DYNAMIC SYSTEMS
WITH APPLICATIONS TO
ROTOR-BEARING STRUCTURES
PROEFSCHRIFT
ter verkrijging van de graad van doctor aan de
technische universiteit eindhoven, op gezag van de rector magnificus,
prof. dr. F. N. Hooge, voor een commissie aangewezen door het college van
dekanen in het openbaar te verdedigen op vrijdag 16 oktober 1987 te 16.00 uur
door
MARCUS THEODORUS MARIA CAOOkiMANS GEBOREN TE DEN BOSCH
Prof.ir. W.L.Egmeijer en
Prof.dr.ir. D.H. van Campen Co-promotor:
de Stichting voor de Technische Wetenschappen (STW).
The re~earch, reported in thi:; di~:;ertation, wa~ :;upported by the Netherlands Teelmoloen Foundaticm.
Chapter Cbapter 2 2. 1 2.2 2. 2. 1 2.2.2 2.2.3 2.2.4 2.2.5 2.2.6 2.3 2.3.1 2.3.2 2.3.3 2.1.4 2.3.5 2.4 2. 4. 1 2.1.2 llbst;rac:t Li:;t of ::;ymbols General introduetion
Numerical method for the determination of ::;table and un::;table ::;tatic, periodic and quasi-periadie ::;olution::;
Introduetion
The determination of :;tatic, periodic and quasi-periadie :;olutions
Some general definit.ions
De::;cription of solutions with a column of parameter::; J Di:;cretization of solutions and resulting equation~ Numerical ::;olution method
Some examples
Discu::;::;ion of the method Stability of the solutions
The eigenvalue problem
Discretization of the eigenvalue problem Solution method
Examples
Discussion of the method
Initia! equilibrium solulions in regular and singular points Regular starting condition
Singular points or bifurcation points
Chapter 3
3. 1
3.2
3.3
A rotor with rubbing Introduetion 3. 3.1 3.3.2 3.4 3. 4. 1 Rotor geomet:ry
Discu::;sion of ::;ome of the de::;ign variable::;
The rotor support, mass eccentricity and constant load
*
The rubbing parameters Fr' a, and f Forced response
ccJeffic:ient
3.4.3 Influence of the boundary for large eccentricities 3.5 Conclusions
Chapter 4 Bearing rnodels for plain journal Learings 4.1 Introduetion
4.2 Some aspect::; of a dimensionless model for plain journal
4.3 4. 3. 1 4.3.2 4.4
bearinqs
Comparison between numerical and existing analytica! solutions Finite difference methad
Analytica! salution methods An alternative bearing model
Chapter 5 5. 1 5.2 5.3
Dynamic behaviour of a rotor-Learing mc~el Introduetion 5. 3. 1 5.3.2 5.3.3 5.4 5. 4. 1 5.4.2 5.5
The rotor model
Dynamic' behaviour of the rotor-bearing model without excitation Quasi-static equilibrium position:;
Stability threshold tor quasi-static s!Jaft positions Periadie solutions
The behaviour of t!Je rotor~odel wit!J excitation Excitation with arbitrary frequencie~
Excitation with ::;ynchronous frequency Discussion and conclusions
Appendix A: Discretization schemes used for discretizing equilibrium solutirm:>
Appendix B: Discretization schemes used for solving the dimensionless Reynolds equation
Re terences Samenvattinq Nawoord Levensbericht
In linear dynamic models, the internal and external forces will depend linearly un ditiplacements, veloeities and accelerations. For som~ pdrticular constructions or loading conditions, such linear dependency does not lead to accurate predict;ion:; for dynamic behaviour. Furthermore, some effects, :.;uch as super- and sub-harmonie resonances and jumping from one motion to another, just cannot be predicted with linear mudels. In t:hi:; thesi:.;, we concentrnted un the stationary motiuns of nun-linear dynamic systems with applications to rotor-hearing structures.
The mathematica! m(•del of a non-linear dynamic :;y5tem cunsists of sets of nun-linear ordinary differentlal equations. Deterministic stationary solut.ions from these equations were time-independent motions, periodic motion::;, ur quasi-perioctic motiuns. A numerical rnethod was introduced for det.ermininy these !:>Ulutiuns and their stability. The met.lwd is based on a discrete time presentation. The combination of thi::; method with an are cunt.inuat.ion met.hod formed a flexible tuol fur inve:.;ti•Jating dynamic behaviuur as a tunetion uf various de::;ign variables.
A typical non-linear model resulted frum the st.udy of a rutur rubbing a.gainst its huusing. The dynamic hehaviour of the model is reported and discussed.
Fluid-filrn bearings are often fuund in rotor-hearing structures and the model:;, generally, are strictly nun-linear. ScJme analytica! and numerical rnodels of plain cylindrical bearings which are frequently used were discussed. A new semi-analytica! model was introduced tltat tm>k into account the effect of cavitation pressure without requiring much rnure computing time t:han t;he exi:.;ting models.
The variuus mudels uf plain cylindrical hearings were used to
investigate the dynamic behaviuur of a model, with two degrees of free~>rn, represent.ing a rotor-bearing st.ructure. The effects of variuus design variables on the linear stability threshold were analysed. The periodic motions that bifurcate at the linear 5tability threshold were determined to show the influence of nun-linear characteristics of the hearing models. Alsu, for a harmonically excit.ed rot.or, their influence on mut.icm and :;tability hehaviour was di::;cussed.
-1 (i) ACM DOf ODE A 0 B
ê
êd
c
D f "'5 Frad Ftan ~·tot Fxba tilde denotes a column matrix an underscore denotes a square matrix
a minus one superscript denotes the inverse uf the square matrix a dot denoteu the first derivative with reupeet to time
a double dot denotes the second derivat.ive with respect to time an index between brackets denotes the relationship with the i-th independent frequency of an m-dimensional equilibrium furrn
a:rc continuation metbod degree of treedom
ordinary differentlal equation
maximum amplitude of rotor response predicted by linear mudel:> bearing length (m)
damping matrix
damping matrix after discretization
clearance between rutating parts and the fixed hou:;ing (m) rutor diameter (m)
fixed cartesian triad rotated cartesian triad
dimensionless rotor eccentricity
column with m independent frequencies of an m-dimensiunal
~tationary deterministic motion
column with ns forced frequencieu column with m-n5 free frequencies column with generalized torces dimensionless constant force
dimensionless coefficient in the expreusion uf the restoring force
dimensionle:;s bearing force in
..
e -directionE
..
dimensionless hearing force 1n e -di reetion
Ijl
dimensionless bearing force amplitude beuring in
..
-directiun dimensionless force e xset of algebraic equa ti ons resul ting aft er discret;izat.i.on
(j "'
h interval distance of an equidistant grid (Ch.2) film th.i.chne:~s
H dimensionless film thickness
K stiffness matrix
~d stiffne~s matrix after discretization
~Jac the Jacubian matrix of a set uf nonlinear algebraic equations m dimeosion of deterministic stationary motion (Ch.2)
dimensionless rotor mass (Ch.J,Ch.Sl md number of discretizatiun points
~ mass matrix
~d mass matrix after discretization nb total number of discretization points np dimeosion uf auymented state space nr number uf design variables
n
5 number of independent periuds in a ~et of ODE's nv number of generalized co-ordinates
P dimensionless pressure
Pcav dimen~ionles~ cavitation preusure P
9 geometrical centre uf disc Pm centre of disc mass
~ generalized co-ordinates
~u stationary deterministic solution of the set ut ODE's with respect to time t
~ design variables
R rotor radius (m)
S Sommerfeld numher
t time (s)
y column with perturbatiuns of the generalized co-ordinates Ye complex eigenform
~ statiunury deterministic sulution of the set uf ODE's with respect tu parameters stured in column ~
V squeeze velocity amplitude (eq.(4.11)) Vrad dimensiunless squeeze velocity
x,y
* *
x ,y z ....z
8 V adi:;placements of t:he rotor centre (m)
dimensionless di~placements of the rotor centre column wi th unkn,)wns af ter discretizii t i1)n dimensionless axial co-ordinate
exponent in expres~ion of the restoring force length diameter ratio
parameter that indicate~ whether the rotor touche~ the twu~ing
dimensionless mass eccentricity
dimensionless damping of the rotor support viscosity (Ch.4)
dimensionless co-ordinnte of the rob•r centre (Ch.5) anyle defined by eq. (3.8) (Ch.3)
circumferential co-ordinate (Ch.4l
design variable in the Van der Pol equation design variable in the Duffing equation cumplex.eiyenvalue
dimensionless excitation frequency
dimensionless co-ordinates of the rotor centre modified Sommerfeld number
dimensionless time
column with m parameters that describe a m-dimensional deterministic stationary motion
radial frequency (Ch.2)
dimensionless rotor speed (Ch.3) dimensionless frequency at A0 rotational speed (=0.5Q)
constant rotational rotor speed (rad/s) dimensionless rotor speed
Chapter 1 General introduetion
The non-linear dynamic analysis of rotating machinery such as pumps,
generators and fans is of increasing importance. As early as in 1965, Tonül publisbed some experimental results in which he ob~erved buth jumping from one stationary motion to another and sub-harmonie resonance. He succes:1fully explained his results with analytica! solutions of simple non-linear models.
Due to increased running speeds and the need for further optimization of the structure, we more often have to deal with problems that can only be explained by analyzinq non-linear models. Also the increased numerical facilities offer the possibility to carry out analy~es that are more
detailed and ar:curate t:han the analytica! techniques. At this point, we will start by giving a few characteristic examples selected from the extensive rotordynamics literature that illustrate the necessity of non-linear dynamic analyses fur rotor-hearing structures.
In rotordynamics, fluid film bearings often have the strongesl:
non-linear characteristics of the rotor-benring structure. The prediction of a static position of a rotor mounted in fluid film bearings is unreliable if we do nut take into account the non-linear characterjstics. Furthermore, experiment:; with rutor-bearing structures [Kirk,Gunter (1974),
Choy,Gunter,Allaire (1978)] very often showed resunances at fractions or multiples of the excitational frequencies which are typical of non-linear
system~:~.
Fluid film hearings can cause the so-called half-omega whirl
instability of the structure as is clearly described by Crandall (1984). The influence of misalignment of the hearings [Hori,Uematsu (1980),
Sprinqer,Ecker,Gunter (1984)] and the magnitude of the excitatiunal force:; [Dannister,Makdissy (1980)] on the stability threshold can only be predicted by nun-linear rnudels. Finally, experiment:; with fluit.! film hearings
(Choy,Gunter,Allaire (1978), Holmes (1980)] showed that the occurrence of the well-knuwn half omega whirl instability did nut always irnply that the dynarnic behaviour is not acceptable. Whether this is the case or not depends
strongly onthefluid film bearing type (Gunter,Humphris,Springer (1983)]. Thè interaction between rotating and non-rotating parts is indicated by rubbing. Fractions and multiples of the excitational frequency were experimentally observed [Muszynska (1984)) and analytica! analyses
[Muszynska ( 1984), Childs ( 1982)] predicted these re::mlts. These phenomena have been discussed and classified by Beatty (1985).
Our last example is that for the calculation of a steady state motion of a structure with squeeze film dampers. Squeeze film dampers do not have stiffness properties. Only because of the non-linear damping characteristics will the shaft vibrate close toa centred position [Greenhill,Nelson (1981), Holmes (1983)].
In this thesis, we will stress the importance of non-linear dynamic analysis and show new calculation techniques and mudels that can be used for the dynamic analyses of rotating machinery. As the set-up of the numerical methods is quite genera!, they can also he applied to general type:; of non-linear dynamic systems.
In this introduction, we will review some general characteristics of linear and non-linear systems and discuss the typical non-linear phenomena that cause linear prediction methods to fail. Furthermore, some remarks are made on the availability'of non-linear models for rotor-hearing structures. Finally, the contents of this thesis are summarized.
In non-linear analysis, the torces acting on or between structural parts can be modelled mathematically by arhitrary non-linear functions of the
generalized co-ordinates g and their derivatives, 9 and g, with respect to time t. Furthermore, an arbitrary number of design variahles, collected in the column~· and timet can play an explicit role in the dynamic behaviour. In genera!, the mathematica! model of a rotor-hearing structure will lead to a set of non-linear second order ordinary differential equations (ODE':;):
( 1 . 1 )
Here, the equations have been divided into a part
!
2 that depends explicitly on time (and nut on any of the generalized co-urdinates), and a part&
1 with all other contributions. The right-hand side is taken into account when we investigate forced vibrations of the structure.In a linear analy~is, we assume that the forces acting on or between structural parts linearly depend on the generalized co-ordinates, their veloeities and accelerations. For real structures, this is nut true but such linear models are often sufficiently accurate in the range considered here. The modelling ~tage then results in a mathematica! model that consists of a set of linear second order ordinary differentlal equations (ODE's):
( 1.2)
This equation remains linear even if the matrices ~.
ê
and ~ are time-dependent; huwever the treatment of such an equation in rotordynamics i:;less common.
The structural parts are modelled in terms of their contributions to the influence coefficients in the so-called mass matriK ~. the damping matriK ê and the stiffness matrix !S of these equaticms. The linear model can be derived from the non-linear model but will only be valid in a certain range of generalized co-ordinate values.
We will now make a few remarks on non-linear analysis that show it is quite different from linear analysis.
1) In the absence of eKplicit time-dependent forces, it is possible to find equilibrium forms in which time has an influence.
These equilibrium forms can include phenomena such as limit-cycles or periodic solutions, quasi-periodic solutions and stranye attractors (Thompson,Stewart (1986)].
The equilibrium forms can be distinyuished by means of their dimensions. In Lichtenberg and Lieberman (1982) a common definition of the dimension is given. In this definition, a static uolution has dimension zero, a periadie solution has dimension one and a quasi-periodic solution has an integer dimension of two or more. Typical of strange attractors is that they have a fractional dimenHion in this definitiun.
2) More than one equilibrium farm may elÜst for the same set of design variables.
In a linear analysis, we usually have one static solution; but in non-linear analysis only rarely can it be proved that the calculated solution i:; the only possible solution.
3) non-linear models can be made valid for a larger range of the design variables .
. In linear rotordynamics, the coefficients in the mass-, damping- and stiffness matrices are generally determined in an assumed steady state equilibrium of the rotor-hearing system. For a rotor mounted in fluid film bearing:;, such an estimation must be quite accurate. Small chanqes of the steady state position can alter the coefficients quite a lot. Generally linear analyses are justified if the rotor is slender and tlle bt:!ar.ings have some distance apart.
In a non-linear analysis the non-linear equations (1. 1) are related to a set of globally generalized co-ordinates. A static t:!quilibrium position i~
found when we sulve (1.1) and neglect all parts that depend, directly or indirectly, on time" t. This means that accelerations, veloeities and excitational furces are made zero. For investigating the influence of misalignment of hearings on the ~tability of a rotor-bearing system, such a non-linear analysis is necessary [Hori,Uematsu (1980), Springer,Ecker,Gunter (1984)].
In order to illustrate the last two remarks we will consider the Van der Pol equation:
.. 2
q + ~ (q -1) q + q 0 ( 1. 3)
The solutions of this equation can shown in the phase :;pace a:; a tunetion of the system parameter ~- It can be easily seen that for all values of IJ the salution q=O is valid. Additionally, for ~=0, the pha:;e :;pace i5 filled with periadie solutions, all with the same shape and period but with tlif(erent
amplitudes (figure 1a). For all other values of IJ, we can only find one additional periodic salution (figure 1b).
6.o r 6.0
MU.O.O
-6.0 -6.0
figure 1a figure 1b
Fig.1.1: Periadie solutions of the Van der Pol equations
MUoO.l8<12
3.00
-MU<!.l827
Non-linear dynamic model:; of structural parts are more complex, both in formulati.on and in application, than linear ones. However, when it is not possible to obtain a linear model that is valid over the entire range of interest or when typically non-linear phenornena have to be predicted, we must use non-linear model:;. At present, the generation of non-linear models based on experimental data iu difficult and time conuurning. Due to this fact, non-linear rn•>dels generally originate from analytica! or numerical studies. Verification of these models with experirnental results can give poor quantative results. However, with the application of non-linear models, it is very often possible to predict the obeserved phenomena qualitatively.
We will discus:; some of the parts of rotating machinery for which non-linear modelu are available and often necessary. The rotating part of the structure is modelled very often with flexible beams and rigid discs. The influence of geometrical non-linearitles due to deformations of.a Timoshenko beam element were investigated by Heeren et al.(1986). They concluded that for a bearn made of steel with a linear elastic material behaviour only for large torques and axial force:;, this influence becarne significant. The influence of non-linearities in the material behaviour on the dynamical respons was investigated by Muszynska (1974). If we restriet ourselves to structures made of steel, we can conclude that the influence of non-linearities in the material behaviour is small.
If we shift our attention from the rotor towards the rotor connections to its casing, we will find several parts that behave Btrongly non-linearly. These elements are the hearings, squeeze film dampers and seals. For plane annular seal~ in pumps, if the eccentricity of the rotor does not exceed 60 to 70\ of the clearance, the restoring farces will depend linearly on displacements, veloeities and accelerations (Childs (1985)]. For larger eccentricities, this is no longer true. For squeeze film dampers, many numerical investigations were carried out by Holmes (1984); many
experimental investigations were performed by Tonneson (1976). The models that were available coincided with models for plain journal hearings for which it was assumed that rotational velocity equalled zero. In tltat case, only the squeezing effect was taken into account.
Bearing models are generally based on the Reynolds equation (with a tew exceptions [Wegman (1980)]). Solving the Reynolds equation with a finite difference or a finite element technique is very time-consuming. For plain journal hearings analytica! models have been developed like the Ockvirk salution and the Sommerfeld salution (Booker (1%5)] that are rapid but les:; accurate. For other kinds of hearings, we usually need to solve the Reynolds equation.
Ball hearings behave very differently from the liquid hearings mentioned above. The clearance between the balls and the inner- and outer ringscan cause trouble and it was investigated by Saito (1985). Rubbing of the rotor against the housing sametimes gives rise to fractional whirl of the rotor and has been investigated analytically by Childs (1982).
Techniques for solving the dynamic analysis problems of rotating machinery are no different from those that are suitable for other structures. For time dependent behaviour such as that resulting from a shock loading or from changes in the rotor speed, numerical time integration is most frequently used. This technique is suitable for finding stable statie, perioctic and quasi-periodic behaviour. However, to get an idea of the sensitivities of these solutions with respect to rotor speed and other design var iable:;, thi::; metbod is laborous. This becomes worse, when we think of the numerous
rotordynamic calculation~ in which periadie and qua::;i-periodic behaviour play a central role, such as for unbalanced behaviour or whirling behaviour.
In order to aave computer time and to create a better understanding with re:; peet. to t.he non-linear b.ehaviour of rotor systerns, an alternative methad has been developed. In chapter 2 the non-linear calculation technique for
the analy~ing of statie, periodic and quasi-periadie behaviuur is
introduced. We :;hall see that in rotordynamics, many of the non-linear phenomena can be explained from the knawledge about periadie solutions and their stability.
The salution metbod facilitates the analysis of ~table as well as unstable model solutions. Based on an extension of the floquet theury
[Muller,Schielen (1985)], a metbod for evaluating the stability of the three equilibrium forms mentioned above is presented, it makes use of the
eigenvalues and eigenforma corre~ponding to these equilibrium forms.
The metbod is capable of carrying out a sensitivity analysis. In rotordynamics, interest in the dynamic behaviour a:; a functiun uf one <>f the design variables is common. A problem that can arise from changing a design variable is that the system can loose its robustness [Newhouse (1978)]. In such a situation, the solutions of the system bifurcat.e and we show how to use the results af tlte stability analyses and bifurcation theory for finding other parts of the family of solutions.
In real rotor-hearing systems, close tolerances between the rotor and the fixed housing are required, but misalignment or excessive unbalanced forces may cause the rotor to rub against the housing. In chapter 3, a simple model for this rubbing phenomenon is introduced that consists of a Laval rotor and a weight.less houuing with a non-linear restoring force. The rotor is forced against the housing by a constant force combined with unbalanced forces.
Friction forceu were modelled with eaulomb's law and writing the equations of motion in a dimensionless form left seven independent dimensionle:;:; design variables.
The non-linear behaviour of this systern was investigated for a rotor that rubbed only slightly against the housing. For two different static deflections, periodic motions of the shaft centre have been shown; the
influence of the design variables on the observed non-linear phenomena is discussed and europared with results from literature.
When the unbalanced forces increased, it was shown that the dynamic behaviour was determined increasingly by the stiffness of the housing. Only
for a small band of frequencies did the rotor hit against the housing several times.
Models of hearings often are the ~trongest non-linear component of a rotor-hearing model. In chapter 1, we derive the minimum number of dimen~ionless variables that is needed fora hearing model. The exi~ting models are based on evaluations of the Reynolds equation and can be divided in analytica! models and numerical models ba~ed on the finite difference metbod or finite element method. We show that one of the shortcomings of the existing
analytica! bearing models is that cavitation pre~~ure is not one of the model's variables. We derive a new model that includes this extra variabie and show that the model's application is not much slower. A comparison with finite difference calculations has been made
In· chapter 5, a very simple rotor-hearing system i~ introduced and it consists of a rotor with two degrees of freedom that i~ mounted in a plain journal hearing. For this model, the influence of the design variables on the linear stability threshold was determined. A compari~on has been made between the various analytica! models. The behaviour of a bearing in the neighbourhood of the linear ~tability threshold ha~ been compared with existing literature. A few of the assumptions found in the literature can be explained with the new method.
In measurements, the non-linearity of the hearing shows itself often as a peak at half the rotational speed in the amplitude-frequency spectrum. We verified this behaviour numerically. We show that. the amplitude frequency response is not only a function of the non-linearity of the bearinq, but may depend on the amount of unbalance present in the rotor too.
Chapter 2 Numerical methad for the determination of stabie and un5table statie, periadie and quasi-periodic solutiuns
2.1 Introduetion
In genera!, a dynarnic model for a rotating structure can consist af a set uf non-linear ODE's tltat can be used to predict the dynamic behaviour as a tunetion of time t and as a tunetion af the initia! utate of the system at
t=t0. The influence of the initia! state generally decreases with time.
The resultiny :.;a-called 'long term' salutians resembie reality depending on the quality of the model, consequently, effective methods for calculating these solutians are important.
A part of the long term solutions are the statie, periadie and quasi-periadie solutions. These three solutions or equilibria cover the majority of the experimentally observed dynamic phenomena for rotor-hearing
structures. Mathematically, these solutions hoive an integer dimen:;ion in common and their properties recur.
In this chapter, the problem of determining these sulutions is divided into two parts. The first part concerns the calculation of stable and
unstable solutions and the second part is about the evaluation of their stability. The numerical technique that will be used in the first part needs a sufficiently accurate starting point for the solution. In the last part, we will show how we can use solutions of the linearized problem and general bifurcat.ion properties for generating these starting points.
2.2 The determination of statie, periodic and quasi-periodic solutions
In thi:.; parayraph, we will introduce a numerical tedmique that is able to
determine equilibrium furms with an integer dimeosion such as statie, periodic and qua:.;i-periodic sulutions. Tt is essential that we give a discrete approximation of the time-dependent equilibrium forms. The basic eperation is the transformation of thesetof non-linear ODE's (1.1) into a
set of non-linear algebraic equations. These can be solved with standard
numerical procedures. Introducing an are continuatien method, we arealso able to carry out sen:;itivity analyses with re:;pect to one of the design variables. As an example, we will determine periodic solutions of the
Duffing equation with damping and the homogeneaus Van der Pol equation. Finally, the metbod is discussed and compared with other methods.
2.2.1 Some general definitions
We will consider the non-linear equations of motion mentioneu in chapter 1. We must assume that we have nv equations with nv generalized co-onlinates and nr design variables. The equations of motion may have periodic
coefficients with a finite number of independent frequencies, ns. Those frequencies can be regardedas special members of the nr design variables. The state of this second order system can be described in the augmented state space. The augmented state space is the np-dimensional space of all qeneralized co-ordinates, their velocities, and the design variables (np= 2•nv+nr). A subspace of the augmented statespace is the 2*nv-dimen~i•Jnal
phase space for all qeneralized co-ordinates and their velocities.
Each solution of the equations of motion can be mapped in the
augmented state space. A special class of solutions are those that have
recurrency propertiès [Lichtenberg,Lieberman (1982)]. From the physical point of view, these solutions are the most interesting, because they represent kinds of steady state solutions. The state of the system returns to the vicinity of an earlie~ state after a finite time. If we assume that the map of the l:lolution in the augmented state space covers a subspace 9 that is closed and invariant [ Lichtenberg, Lieberman ( 1982)], we can de fine the dimeosion (lf the solution as:
11(9) lim (ln(M(t)))/(ln(1/~))
t-tO
(2. 1)
HereMis the minimum number of np-dimensional volumes with side lengtbs c that covers the subspace 9 entirely. From this definition it foll(lws that a static solution will have a dimeosion of zero, a periodic solution will have a dimeosion of unity and a quasi-periodic solution, where two period times are involved, will have a dimeosion Df two.
2.2.2 Description of solutions with a column of parameters l
A periadie salution is fully determined if the generalized co-ordinatea as a tunetion of time are known during a period. For qua:.;i-periodic motions, it is not possible to find a finite period in which the salution repeats .. itself. Therefore, rather than descrihing the salution as a tunetion of
time, the salution can bedescribed by a minimum number of parameters. We will distinquish between deterministic solutions with an integer dimeosion and chaotic solutions with a fractal dimension, and restriet ourselves to the deterministic m-dimensional solutions including both static and periadie solutions which play an important role in rotordynamics.
If we look at the solutions mapped in the augmented state :;pace, each point of this m-dimensional subspace can be uniquely identified by m independent parameters t(i) (i=1, ... m), stored in a column l· A trajectory in the augmented state :;pace of an m-dimensional equilibrium form will remain in the corresponding subspace. As time proceeds, the salution points of the trajectory will always be elements of the m-dimensional :.;ubspace and in the end, the complete subspace will be covered. If this is not true, the dimension of the :;olution will not be equal to m. If we describe the
trajectory in terros of the column of parameters
l•
we can choose the parameters in such a way that each of them increa:;es or decreaseacontinuously. If this is not true, one of the parameters t(i) is unnecessary and the dimension of the equilibrium will be les:; than m. Due to the
assumption tbat the subspace is covered and closed, we can conclude that the salution ~must be periadie in each of the parameters l(i). We will define each parameter r(i) in such a way that
~
is periodjc with period one for each of the parameters. The relationship between the real time t and each of the cyclic parameters t(i) can be yiven by a fixed frequency f(i) such that:(2.2)
In the m-dimensional column
!•
the frequencies of the corresponding column of parameters J. are :;tored.A one- or a two-dimensional equilibrium form can be drawn in a three-dimensional augmented state :;pace represented as a closed curve and a torus, as is drawn in tigure 2.1. Fora closed curve, we only need one cyclic
,.<
1'=o
.
s
L""o.2
Fiq.2.1 A one- and a two-dimen~ional equilibrium mapped in
augmented state space with one and two cyclic parameters.
parameter t(1) along the curve in order to identify each point. For a torus,
weneed two cyclic parameters, t(l) and r(2 l, in order to give a unique
identification.
We have already mentioned that the equations of motion (1.1) may
contain periodic coefficients on both the left-hand and right-hand uide:;
with n5 independent frequencies. Therefore, the solution can contain ns
independent frequencies as well. We will assume at least n
5 different cyclic
parameters T(i) which are chosen in such a way that the first n frequencies . . s
f~l)
of!
are identical to the ns independent frequencies of the periodic coefficients of the equations of motion.These n
5 frequencies occur both in the coefficients of the equations
of motion and in the solution, they will be called forced frequencies and
stored in the column!:;· Oue to their system-determined character, it is
better that these frequencies are design variables. All other independent
frequencies in the solution will be called free frequencies and collected in
the column !f· The m-ns frequencies are unknowns and will have to be solved.
2.2.3 Discretization of solutions and the resulting equations
The equations of motion must be written with re~pect to the column of
parameters r. By taking discrete values of the parameters t(i) (i=1, ... m),
the ~et of non-linear ODE's will be transformed into a set of algebraic equation:.;.
For m-dimensional solutions with m greater or equal to one, we can
expre~~ their accelerations and veloeities as derivatives with re~pect to the cyclic parameters. Therefore, we must consiuer a salution 9(tl=9u(t) of the equations of motion corre~ponding to the m-dimensional equilibrium. We
~(~). The veloeities 9u(t) and the accelerations 9u(t) can be written as:
9u(t) m oU(l) m 1: fi :r"'(i) := r f.u:(t) i=1 V i=1 1"'1 "' m m
o
2U(t) m m 9u ( t l=
r f. ( r f . I ( . "') "' ( . l J J :=
r f. ( r f . u'.' . ( 1)) i=1 1 j=1 ) or 1 OT ) i=1 1 j=1 ]"'1 ] "'We can nuw refurmulate the ~et of non-linear ODE's in terms of their derivatives with respect to the cyclic parameters. Generally, we can write:
(2.4)
where the column 1 contains the n cyclic parameters conesponding with the
• '-"S 5
frequencies f(1). The m-n cyclic parameters 1f that correspond to the free
s s "'
frequencies
!t
do not occur explicitly in the equatiuns of motion.Therefore, at any point of the solution, we are tree to assume that Jt=Ç and to u~e that point a~ a reference point for de~cribing the equilibrium. In order to obtain a unique solution, we put the derivative of one element of
g<J,l
with respect to each of the cyclic parameters contained in the column 1t equal to zero atJ.t=9 .
This leads to m-ns additional equations:[ilU(k)/ot(i)]r=O 0 (2.5)
"'"'
Rather than describe the solution U as a continuous tunetion of each of the
"'
parameters r(i) , we will describe the salution as a discrete number of points. Each parameter T(i) is discretized in mái) points t~il, .... , 1(i) in
md (i)
The solution ~IJ) can be described now by a column matrix ~ in which
the values of
!!
at all discretization points are stored. The number ofunknowns nb in
X
ar.e:n = n
* :
m(i)b V i=1 d (2.7)
If wedetermine the equations of motion (2.4), at each discretization point, we will formulate nb equations. The derivatives of this solutiun with
respect to the cyclic parameters,
Yl
andYij
can be expressed by means ofcentral differences. We will apply schemes with a cunsistency of 0(112) and
O(h4) that both use an equidistant discretization grid with width of h. In appendix A both of the discretization schemes are shown. With the examples given in this paragraph, we will compare the two schemes.
We can now formulate a set of nth (nth=nb+m-ns) algebraic equatiuns
for nth unknowns:
The first nb equations consist of the discretized equations of motion (2.4)
while the last m-ns equations will be obtained by discretizing (2.5). The
solution column ~ is defined as:
2.2.4 Numerical solution method
Solutions of the set of equations will be investigated by means of existing salution procedureli. A Newton-Raphson procedure was used fur finuinq a
single solution of the set of non-linear algebraic equations. An are
continuatien method was used for performinq a sensitivity analysi:; of the
salution for changes of one of the design variables. Both the iterative
Newton-Raphson procedure and tl1e are continuatien method are discus:;e<l in
this paragraph.
If values for all the design variables stored in the column
5
arealgebraic equati0ns with nth unknowns !· The initia! estimate, ~
0
, fur a sulution of (2.8) can be improved by the Newton-Raphson proces:; for multiple directions. In this iterative pracess, the (i+1)-th appruximation !itl isfound from the data generated at an approximatiun ~i:
(2. 10)
Calculation of the inverse of the Jacobian KJ - ac fur a solution z. was
~1
carried out at eac:h iteration step. The numerical effort needed to decompose the Jacobian was reduced; a part of the matrix had a diagonal structure. After rearranging the equations (interchanging rows and columns in ~Jac),
*
the .Jacobian was f:ransformed to K :*
K A ~1
~2
c
( 2. 11)In (2.11), the matrixBis a banded uquare matrix whose dimensions are considerably larger than those of the square matrix Ç. Matrix ~
1
and matrix•
~
2
are rectangular matrices that complete the square matrix ~ . The banded matrix can be formed because the set of equations at one discretization point only uses displacement data of a few discretization points backward or forward. Thus both the storage order of the ::~olution ~ and the equations ~1
ands;
2 must be changedinorder to obtain the form (2.11).The accuracy of the solutions can be checked trom the ratio E of the norm of the residual values of the algebraic equations and the norm of the Jacobian. The latter was calculated as the root of its square elements. In genera!, we assumed that this ratio should be smaller than 10-8 .
Although we can solve the problem for one set of design variables,
frequently, we are intere::;ted in the sen:;itivities of the equilibrium form with respect to the design variables. If we choose values for all the design variable::; but one, say ri' then we can reformulate the problem as:
( 2. 12)
The salution ~of (2.12) is now a tunetion of the design variabie ri. This kind of problem can be solved by an are continuation metbod (ACM) which will
be discussed only briefly here. For a more detailed description of the method, we refer to Van der Varst (1982) or Allgower (1980).
The ACH is a second order solution technique for problems whose Jacobians are continuous and known. The method assumes that we know a
solution, say ~
0
, at ri=riO' The method has two parts: 1) a prediebon partthat gives an approximation ~ for the next point ~
1
of the series ofsolutions at ri=ri 1 and 2) a correction part that improves ~ until an
acceptable accuracy is reached.
The way of calculating the predietien direction (öz,ör. ), where öz and ~ l ~ 6r1. are the increments of zand r., is described in Allgower (1980) and
~ l
Fried (1984). The direction is supplemented by the elliptical constraint
where Às is the ~tepsize:
( 2. 13)
In the second part, we have ta calculate a number of corrections ö~i and
iterate until an acceptable accuracy is obtained. fried (1984) contains a survey of some existing correction methods. We will uue his orthogonal
trajectory accession, because it appears to be the most robust metbod
without requiring excessive computer time when compared with the other
aethods. The iterative correction procedure has to be stopped when the
solution satisfie~ an accuracy criterion which is the same as that used in
the Newton-Raphson procedure.
Allgower (1980) described a variable stepsize criterium that indicates
when the stepsize ~s should be increased or decreased. lt is based on the
ratio of the first two corrections ö~
2
and ö~1 in the correction part. Ifthis ratio is smaller than ub' then ~he stepsize should be increased; if it
is greater than at, the stepsize should be decreased. It sometimes happens
that after increasing the stepsize, the convergence rate is so poor that the
stepsize has to be decreased aqain. In that. case, the smallest step:;ize has to be used. The common values for ab and at are respectively 0.2 and 0.8.
2.2.5 Some examples
We will illustrate some of the properties and potentials of the methodology
Befare actuálly treating the examples, we will discuss how bJ judge the quality of a solution. For a set of design variables, we assumed that ~O
was a salution of the equations (2.8) with m
0 discretization points. We determined the solution ~
1
with twice as many discretization points for the same set of design variables, because of the equidistant discretization grid, ~1
contained improved estimates for the m0 discretization values of~o· In the same way, further imprave tl1ese estimates via determining the salution
&
2 with tour times as many discretization points. For each co-ordinate, now we can define a convergence ratio re as:mo (4i) (2i) [ z - z1 i=1 2 ( 2. 14al r c mo (2i) (i) r. z - zo i=1 1
Theoretically, r is equal to 1/4 for the O(h 2J scheme and reis equal to 1/16 for the O(hij) scheme. The global error öa of one of the co-ordinates was calculated as:
·~a lim n-+ .. with a 1 (2i) ~1 - z "'0 (i)
Bath the convergence ratio and the global error can be used to evaluate the accuracy of a solution.
Example Duffing equation with damping
A simple example that shows many of the characteristic difficulties of non-linear dynamic problems is the Duffing equation with damping:
The design variables are the damping coefficient D, the cubic stiffness coefficient ~ (mu) and the excitatiunal frequency f. For most of the values of the design variables, this system has periodic solutions (m=1) with a forced frequency f (n
5=1).
In tbis example, we take 0=0.01 and ~=0.04. Duriny tbe are
continuation proces::; we varied the frequency f from 0.005 to 0.56. In figure 2.2, the amplitudes of the periodic solutions are drawn ayainst the design variabie f. The characteristic Duffing response can be observed, a5 well as, the superharmonie resonance top at f=0.055 Herz (1/3 f0J. Using the are continuatien method, we were able to follow the sharp resonance peak. easily.
Figure 2.2 is based upon about 140 calculations; hence, we used about 140 calculation steps. Toeach solution a unique number was assigned, called the stepnumber. We needed this extra number because we might have multiple solutions at one value rJf the design variabie. Some of t:he stepnumbers are shown in the figure. For description of the periodic solution, 20
discretization points were used. The results of both t:he O(h 2J and the O(h 4J difference scheme are shown but a difference cannot be seen from figure 2.2. In figure 2.3, 'three periodic solutions that were calculated with the O(h2J discretization scheme have been drawn in the phase ::;pace. The symbols in this tigure indicate the discretization points. We can see that the number of points seems adequate to describe the periodic solution.
If we focus on the superharmonie resonance peak in tigure 2.2, we will obtain fiyure 2.4. Here, we can observe a difference between the results obtained with the O(h2) and O(h4J discretization schemes. For three stepnumbers of fiqure 2.4, we have drawn the solutions obtained with the O(h2) discretization scheme and mapped these in the phase space as shown in fiqure 2.5. Here, we can see the influence of the triple harmonie in the solutions. It can be seen easily that the curves are point symmetrie with respect to the origin.
In figures 2.6 and 2.7, the ylobal errors and the convergence ratios of both the O(h2J and the O(h4J discretization schemes are drawn. The convergence ratios that were calculated, closely match the theoretica! one::;.
"'
0 lflO 15.0 0.\\fPI~G OI Clll ST~ /kl ;.> 10.0s
~
5.0 .1L'10RO.DS. u.o 0.00 1.00 ~00 5.00 6.00Fig.2.2: Amplitude vs. frequency of the Duffing equation with D=0.01 and ~=0.04. The results obtained with the O(h 2l and O(h 4) discretization schemes have both been drawn.
60. DA~!Pl'IG .OI ctli.Sn< .C» 40. 20.
,.
-ZO.O~
ZO.O,.
-20. EQPO!Yf'"'l 6 6 6-dil. FJ>POI:'"fo60 Q
"'
"'
F.QJ'OI!'o> -.?O
-60. 0 9 9
Fig.2.3: Periadie solutions of the Duffing equation mapped in phase :;pace at three solution steps.
1.600 32. .tH 1..100 CLll.S"P.i .()l 1.400 34 1.300
"'
Q5
~ uro 1.100 I.OCO ~-JORO 05. 6 6 6 4tb ORO.OS. 0.900 3.00 4.00 l.OO 6.00 7.00 8.00 91)J FREQUENCY HzFig.2.4: Amplitude v:;. frequency of the Ouffing equation ohtained with
both
an
0(1/) and an O(h4l di:;cretization :;cheme (D=0.01,~=0.04). 1.()0 DAMPI:'G .01 CliB.STN .()1 -2.00 2.00 EQ.POLYT·lO 6 6 6 EQPOf:<f,J2 ? Q Q EQroc-cr•J·• -1.00 0 8 8
Fig.2.5: Periadie ~olution~ of the Duffing equation fur three ~ulution
.J •JO 400 )j() 3.00 2.50 2.00 1~0 2nd ORD.DS. 6 6 6 1.00 niEOREllCAL 7 9 9· .W.ORD.DS. 0~0 0 0 0 nJEO~EllCAL 0.00 dOO lO.O 60.0 70.0 80.0 90.0 STEPNUMBER
Fig.2.6: The theoretical and calculated convergence ratios for the
re:;ults :;howed in fjgure 2.2.
1.00 0~0 :te: 0.60 0 ~ ..J <
3
'-' OAO 0.20 0.00 1\0.0 50.0 W.O 70.0 Sll'PNUMD~R 80.0 90.0 2nd ORO.OS. "' 6 6 4lhORO.OS. '<1'- --9--...q..Fig.2.7: <.aobal errors vs. stepnumbers for the results showed in figure 2.2.
Example 2 The Van der Pol equatiun
The next example wil! concern the Van der Pol equation. Some of the salution forms were given in figures 1.1 and 1.2, but we will concentrate on the periodic solution (m=1) withafree frequency (n =0). The equation has only
!)
one design variabie ~. In tigure 2.8, the amplitude of the periadie
solutions vs. the design variabie ~ is drawn. Theoretically, the amplitude ought to be near the value 2 and it was found for the solutions obtained with 20 discretization points that the numerical error became significantly large when ~>1.8. The period time vs. ~is shown in tigure 2.9. A large
difference can be observed between the solutions obtained from the 011121 and O(h4) discretization schemes. In order to see what happens to the solutions mapped in the phase space when ~ increased, tigure 2.10 was prepared. From this figure, it is clear that 20 discretization points are not sufficient tt>
describe accurately the periodic 5olution for large values of u. lf we increase the number of discretization points, it can be seen from tigure 2.11 that the approxirnation impraves.
The convergencè ratio and the qlobal errors corresponding tu the results shown in tigure 2.8 have been drawn in figures 2.12 and 2.13. Convergence ratios above ~~3 becorne worse which indicates paor results.
Obviously, the global error obtained with the O(h2l discretization scheme is larger than that obtained with the O(h4) discretization scheme.
2.2.6 Discussion of the metbod
Althouqh a general theory for integer m-dimensional equilibrium forms has
been given, the examples presented there concerned only one-dimensional equilibria. We have two good reasons for this. Firstly, the non-linear periadie solution is a first step in non-linear dynamics and its physical meaning is clear. Secondly, zero- and one-dimensional equilibrium forms seem to be of primary interest for rotordynamics.
Although we have not yet determined the stability, frum literature, it is known that same of the solutions that were found are unstablc. The
stability of the perudie solutions that were found, seemed to have no influence on the salution process.
2.010 2.!XXl 1.950 1.900 1.750 1.71.0 1.650 0.00 0.50 1.00 1.50 2.00 Or'-'i.VAIU1U V 2.50 3.00 ).50 2nd ORD.lJS. 6 .!> .!> 4th ORD.DS. 9----'J----
SJ-Yig.2.8: Amplitude of the periadie solutions of the Van Der Pul equation vs. design variable ~. ealeulated with 20
diseretization points and 2 different diseretizatiun sehemes.
1.7UO 1.650 1.600 l.450 !.4(tl 0.00 0.50 100 150 'V: --'i' ---9- -- --'i' - - - -9 - ---V 2.00 Df...S.VI\R.MU 2.50 ).00 3.50 ~d ORD.DS. 6 .!> 6 4lh ORD.DS.
9----<:J.---SJ-Fig.2.9: Frequeney of periadie solutions of the Van Der Pol equation vs. design variable ~. ealculated with 20 discretization points and 2 different discretization schemes
----
.. ',9.·~-~-:~···...
\
~''\-,V ~11:=0.7)197) t. 6 6 Mli•l.91974 9---~----<r MU•l.47291 G-···0···0F~q.2.10: Periodic solutions of the Van Der Pol equation mappP.d in pha!:le space calculated with 20 discreUzation points and for three·different values of~.
6.0 2D -6.0 ··Q ... '<\~-v
..
' --~:\
:Ç 2.00 3.47291 20PNT (', 6 ... 40PKT 9---~----<r 80PNT G-···0···0Fiq.2.11: Periadie solutions of the Van Der Pol equation at ~=3.47291 mapped in phase :;pace, calculated with 20, 40 and 60
5.00 4.00 ~ < " ~ 3.00 8 1.00 ((I() 0.00 0.00 0.50 1.00 150 2.00 DF.S.VAR.MU 2.50 3.00
'
41.h ORD.I.>S. < ) . - - -_ " . - - - "J-3.50Fig.2.12: Convergence ratio for reaulta showed in figure 2.8.
"
0 ~ "' ~ <..
g 0 350 3.00 250 2.00 150 1.00 050 0.00 0.00 0.50 1.00-
"
,1?' 150 2.00 DES.VAR.MU 2.50 ? ' 3.00 3.50Fig.2.13: Global error fur results showed ~n fig.2.8.
2nd ORD.DS.
4th ORD.DS < ) . - - - _ " . - - -
"J-In the Duffing equation, we were able to find superharmonie re:;onance. This phenomenon cannot bedescribed by a linear model. The :;harp peaks indicating low damping terms were described adequately by the present lay-out of the are continuatien method (ACH), even when we used a fixed stepsize.
In the example of the periodic solutions of the Van der Pol equation we can see that the stronger the relaxation vibration, the less able is the metbod to fulfill accuracy demands. Large changes in the motion·tend to occur in an ever smaller part of the periodic vibration. In order to
maintain accuracy, the number of discretization points must become extremely large. Local refinement of the discretization grjd may Le a solution, but it has not been applied yet.
As mentioned before, the ACM is a second order :;olution methad that demands a continuous Jacobian. This means that a sudden change in the stiffness can not be treated which is a big disadvantaye when we want to describe the rubbing phenomena when bouncing of the rotor occura thal hita the housing.
When we applied the method to the Van Der Pol equation, no implicit warning was given for"the bad reaults obtained. The convergence ratio and global error had to be checked separately. This means that, in general, we have to be caretul with obtained results if global errors are not checked. From the examplea, it becomes clear that the O(h 4J diacretization scheme is preferred. The good convergence rate diminishes the global error much faster when the number of discretization points increaaes. Furthermore,. the global errors of the results obtained with the O(h 4l diacretization scheme and 20 discretization points were smaller than those obtained with the O(h 2J scheme and 40 discretization points.
At least two other methods are available for calculatiny m-dimensional equilibria.
1: The shooting method is basedon numerical time integration. By
making small changes to the starting solution, a relationship can be
obtained between these small changes and the changes after a period of time. Using this relationship an iteration process can be developed that chooses the starting conditions so that it fulfils the boundary condition of the system. The numerical atability of the method was shown to be good by Ling
disadvantage of the metl1od is tl1at it is difficult to extend to :>olutions with a larger uimension than unity.
2: Tl u~ harmonie balance method is based on minimizing the error function. Solutions are described as Fourier series. The dimension of a St>lution may be larger than unity according to Lau, Cheung and Wu ( 1982, 1983). For each Fourier coefficient in the solution, the non-linearities were described in term:; of Fourier coe(ficients by mean:; of integrals. A uisadvantage is that these integrals must be calculated for each Fourier coefficient at each i teration step in the solutitm process. The metlwd ha:; some very important advantages too; each coefficient may be calculated separately and each adds to the accuracy of the solution. After each step, we can merely chouse which Fourier coefficient should be added to the solution until we are satisfied and stop the process.
2.3 Stability of the solutions
From knowledge of the statie, periodic and qua~i-peri<,dic ~olutions alone, we do not know whether these solutions are actually important. Therefore, we
mu~t con~ider the behaviour of the ~y~tem in the auymented state ~pace for small perturbations with respect to the equilibrium solution. If the solution converge~ to the equilibrium ~olution for t~- and for all deviations that are sufficiently small, we can call the solution locally stable. Otherwise, we will call the ~olution unstable. More precisely, we will use the :;tability definitions according to Lyapunov (Muller, Schiehlen
(1985)].
Stable solutions are called attractors too. Factually they are the most significant because they persist even after a long time. A more detailed definition of an attractor that include~ statie, periodic and quasi-periodic solutions was given by Lichtenberg,Lieberman (1982) and Kreuzer (1985).
The local ~tability of a ~olution is determined by the linearized properties of a systèm in the equilibrium state. The stability area, the area of the augmented state space from which trajectorie~ converge to a certain equilibrium, is even more important than the local stabilily itself.
Sometimes, it is possible to find the ~tability area ~uch a~ for the Van der Pol equation. In genera!, the stability area is difficult to determine.
The ~tability of the equations could be investigated by means of a numerical integration. Starting with a small deviation at time t=t0, the trajectory near an unstable equilibrium will diverge from it after a finite period of time. A trajectory near a ~table equilibrium will converge toward~ it. This is a general metbod and easy to apply. However, the choice of the initia! state of the system at t=t0 may influence any conclusions that are drawn. Therefore, for a ~ufficient. number of initia! st.ates, the
trajectorie:; have to be calculated. For low damping rates, the period of time over which the~e trajectorie~ have to be calculat.ed will be very large.
Furthermor~. it is difficult to define the speed that a trajectory converges to the equilibrium ~olution.
In this paragraph, we will ~h<~ that the local stability of periodic and quasi-perioctic ~>tions can be found in much the same way as for static equilibrium points using Lyapunov coefficients. We will use a discretized
form of an eiyenvalue pwblem that links to the discretized form oJf th~ ~alution~ in paragraph 2.2. Thi~ adju~ted eigenvalue problem can be ~alved with the well-known eigenvalue algorithms. An example of calculat;ing
~tability af periodic salutions af the Van der Pal equation i~ given. The
method is discu:;:;ed in the last part of thi:; parayraph.
2.3.1 The eigenvalue problem
The m-dimensional :>olution 3 can be considered as a tunetion of the m parameter~ af the column J. and tlae m carre~ponding frequenc:ies af the column
!·
We will formulate an eigenvalue problem for this ~alutian that links to the discretizatian methodolagy described before.The linearized properties of the system in this ~u can be found after linearizing the equations of motion with respect to diuplacements,
velocitie~ and ac:celerations relative to the equilibrium form. For a time
dependent equilibrium solution, the reslultant mass-, damping- and
~tiffnes:.; mat.rice:.; will depend on time a~ well. Int.roducing a column with urnall perturbations y of the generalized co-ordinates relating to the equilibrium, we can build a set of linear ODE's:
~(!<!.l~ + B(f,1)U +
-
""' ""' "" !5(!<!.'~g
with M .. lJ oF 1 1 .foqJ . B .. lJ oF1l .foq. J K .. lJ=
oF11 .foqJ .The ruw index i and column inûex j indicate an element of the mas~-.
· damping- ur sliffnes:; matrix. The values of the m independent frequencies
stored in column
!
require knowledge of the equilibrium.The stability of this linear set of ODE's and thus the local ~tability
of the curresponding equilibrium will be determined by means of eigentorros and eigenvalues. Both for zero and ane-dimensional equilibria, the theory is well-developed and the coefficients for the curresponding of linear ODE's are con~tant and periodic, resp. For the periadie caefficients in (2.16), we refer to the Fluquet-theory [Cesari (1963), Muller,Schiehlen (1985)].
In case of either constant or periadie coefficients, the funJamental
motionsof equation (2.16) are de:;cribed by a (periodic) complex. matrix
Q
and a con~tant complex matrix B:
u(t) = Re[ n(t) e~t c]
~ ~ ~ (2. 17)
In column c, n complex constants are stored witll arbitrary value~. For
"' V
~onstant coefficients in equation (2.16), the complex matrix Q have to be
constant and for periodic coefficients with frequency f, tl1e eomplex matrix
Q
has to be periadie with frequency f too.The multiplicity of an eigenvalue is given by the number uf the
eigenvalues with that same value. The nullity of an eit]envalue 1~ yiven by
the minimum number of independent eigenfarms that belong tu eigenvalues with
the ~ame value. If the nullity equals the multiplicity for all ei•Jenvalue:;,
the matrix ~ can be transformed into a diagonal matrix. We will assume that
this is true for the systems here.
In that case,
ê
can be written au B=
V Ay-
1 where ~ is a diagonalmatrix containing the eigenvalueu of ~ and
Y
is a matrix containing thecorresponding eiqenvectors of ~- We then qet:
( 2. 18)
We can divide this fundamental salution (2.18) into 2'nv fundamental ~1tion~
that look very much alike:
( 2. 19)
For the k-th fundamental motion, the either constant or periodic, complex
•
column u (t) will correspond to the k-th column of the matrix Q (=Q '{) and ~
"'e
is the k-th diagonal element of the matrix 6·
In the case of periadie coefficients, the complex constants ~ in
(2.19) are not defined uniquely. The :;ame fundamental motion yltl can be
defined with an infinite number of combinations of );!ètl and ~. The
combination that we choose will depend on the accuracy with which we can
determine a certain.cornbination. It may be that the eigenvalue ~ is
