RODY II: rotordynamica binnen Matlab

RODY II

Citation for published version (APA):

Schie, van, C. (1991). RODY II: rotordynamica binnen Matlab. (DCT rapporten; Vol. 1991.052). Technische Universiteit Eindhoven.

Document status and date: Gepubliceerd: 01/01/1991

Document Version:

Uitgevers PDF, ook bekend als Version of Record

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.

RODY

I1

Rotordynamica binnen Matlab

TUE, Faculteit Werktuigbouwkunde

WFW Rapportnr. 91.052

Cees van § d i e

Instelling:

TU-Eindhoven

Faculteit:

Werktuigbouwkunde

Begeleider:

A. de Kraker

Plaats:

Eindhoven

Datum:

22

augustus

1991

Samenvatting

De analyse van het dynamisch gedrag van een rotor moet met een computer voor iedere

geïnteresseerde ingenieur (gebruiker) éénvoudig mogelijk

zijn. Een aantal functies moeten

daarom voor de gebruiker standaard aanwezig

zijn. Daarbij moeten niet ter zake doende

dingen, zoals opslag van gegevens, zoveel mogelijk worden vermeden.

Het bestaande programma

RODY

is uitgebreid, waarbij de boekhouding van

de

gegevens

door

het

p r q p ~ m a

wedt

afgekandeid. De gebrdkeï

kaì

via

een menusimeiuur

opdrachten en invoer opgeven, waarna de gewenste berekening wordt afgehandeld.

Het programma-pakket

RODY

11

is een gereedschap, waarmee een dynamische analyse

van een éénvoudige rotor kan worden uitgevoerd.

Het is in eerste instantie niet nodig Matlab goed te beheersen.

Als de berekeningen

complexer worden moet de gebruiker meer van deze programmeertaal kennen.

O.

Opdracht

Maak op een relatief éénvoudig wijze het mogelijk het dynamisch gedrag van een rotor te

bepalen. Schrijf een programma in Matlab, zodat een aantal standaard gegevens als:

-

Eigenfrequentie

-

Eigentrillingsvormen

-

Excitatie-responsie

van een rotor kunnen worden berekend.

Daarnaast een overzicht van de rotor en figuren zoals:

-

Campbell-diagram

-

Eigentri

11

i

ngsvorm

-

Bode- en Nyquist-figuur

belangrijk voor analyse van de rotor

Ga uit van het programma

RODY,

geschreven in

1989

door J.A.H.M. Jacobs voor Philips

Natuurkundig laboratorium.

1.

Symbolenlijst

Symbolen

systeemmatrix

(n,

,

n J

ingangsmatrix (n,

,

n,,)

uitgangsmatrix (ny

,

n J

ciempingsmatrix

demping

[Ns/m]

of [Nms/rad]

vector

met externe belasting [NI of

[Nm]

frequentie [Hz]

overdrachtsmatrix van

Q(s)

naar

x(s)

stij fheidsmatrix

stijfheid [N/m] of [Nm/rad]

massamatrix

vector- of matrixlengte, aantal vrij heidsgraden

ingangsvector in frequentiedomein

ingangsvector met

nu

elementen

eigenvector, eigentrillingsvorm

[m]

of [rad]

toestandsvector

in

frequentledomein

toestandsvector met n, elementen [m],

[rad],

[mls] of [rad/s]

uitgangsvector in frequentiedomein

uitgangsvector met

ny

elementen

[m]

of [rad]

vector van vrijheidsgraden [m] of [rad]

eigenwaarde [rad/s]

exponentiële dempingsfactor

'gedempte' eigenhoekfrequentie [Hz]

dimensieloze dempingsfactor

[-I

hoeksnelheid [rad/s]

Y ~ o ~ ~ ~ ~ ~ .(Ey S

,

rEB> l , ~ ~ ~ ~ ~ ~ 2 ~ ~ ~ Ä

Indexen

U

Y

mode, eigenwaarde

element (lokale systeem)

globale systeem

imaginair deel, rijnummer

kolomnummer

eigentrillingsvorm (mode)

reëel deel

systeem

ingang

uitgang

2.

Inhoudsopgave

Samenvatting

O.

Opdracht

1. Symbolenlijst

3. Inleiding

4. Algemene theorie

4.1 Lokale beschrijving

4.2 Globale beschrijving

4.3 Assemblage

4.4 Toestandsbeschrijving

4.5 Eigenwaarde-probleem

4.6 Systeembeschrijving in frequentiedomein

5,

Theorie

binnen Matlab

5-1 Eigenwaarden

5.2 ûverdrachtsmatrix

5.3 Whirl

6. Gebruik van

RODY

I1

6.1 Mogelijkheden

6.2 Algorithme voor het gebruik

6.3 Invoeren van gegevens

6.4 Functies van Matlab

7. Voorbeeld

7.1 Rotormodel

7.2 Modes

7.3 Excitatie-responsie

8.

Nabeschouwing

8.1 Gemaakte aannamen

8.2 Overwonnen moeilijkheden

8.3 Wat is nog niet zo handig

8.4 Nog te overwinnen moeilijkheden

8.5

Ondervonden problemen met Matlab

9. Programma-pakket

RODY

I1

10.

Conclusie

1

2

3

5

6

6

7

7

8

9

10

12

12

13

13

15

15

15

16

16

18

18

19

21

23

23

23

23

23

24

25

26

11.

Literatuur

27

Bijlage

I:

Basiselementen

1.1

As-element met

8 vrijheidsgraden:

’SHAFTS’

1.2

As-element met

4

vrijheidsgraden:

’SHAFT4’

1.3

Balk-element met

4

vrijheidsgraden:

’BEAM4’

1.4 Schijf-element met 4 vrijheidsgraden:

’DISK4’

1.5

Schijf-element met

2

vrijheidsgraden:

’DISK2’

I.%

Massa-element met

1

vdjheidsgraad:

91wASS9

11.1

Veer met

4

vrijheidsgraden:

’SPRING4’

11.2

Veer met

2

vrijheidsgraden:

’SPRING2’

11.3

Veer

met

2

vrijheidsgraden:

’SPRING2M’

11.4

Veer met

1

vrijheidsgraad:

’SPRINGl’

11.5

Demper met 4 vrijheidsgraden:

’DAMP4’

11.6

Demper met

2

vrijheidsgraden:

’DAMP2’

11.7

Demper met

2

vrijheidsgraden:

’DAMP2M’

11.8

Demper met

1

vrijheidsgraad:

’DAMBI’

Bijlage

11:

Verbindingselementen

Bijlage

111:

Rotorvoorbeeld

Bijlage IV: Programma’s

IV.

1

Eigenwaarde-probleem oplossen:

EIG-V0RM.m

IV.2 Globale variabelen:

GLOB VAR.m

IV.3 ûverdrachtsfuncties bepale;

OVERDRm

IV.4 Hoofdprogramma:

R0DY.m

IV.5 Dummy programma:

SHELL.m

IV.6

Modes tekenen:

TEKEN.m

V.l Gegevens invoeren:

IN-CHECK.m

V.2 Sorteren en controleren:

SAMEN.m

V.3 Elementen in tabel tonen:

D1SPLAY.m

V.4

Vrijheidsgraden controleren:

D0FCHECK.m

V.5 Assembleren systeemmatrix:

ABMDK.m

V.6

Eigentrillingsvormen tekenen:

TEKENF1G.m

Y.7

E!!@

erd delen:

VERnXEL2.m

V.8

3D-figuur transformeren:

TRANS3D.m

V.9 Whirl bepalen:

WH1RL.m

V.10

Plot voorbereiden:

PL0TPREP.m

V.ll Excitatie invoeren:

EXCF0RCE.m

V.12

ûverdrachtsfuncties tekenen:

RODYB0DE.m

Bijlage VI: Formaat matrices

VI.

f

ûverdrachtsmatrix

VI.2,3 Frequenties bij excitatie

VI.4

Rotorfrequenties

bij

eigenwaarden

VIS Ingangsvectoren

VI.6

Uitgangsvector

VI.7

Eigenvector

VI.8

Eigenfrequentie

Bijlage V: Functies

I.

1

I.

1

1.5

1.8

I.

10

1.12

I.

13

11.1

11.1

11.3

11.4

11.5

11.6

11.7

11.9

11.10

111.1

IV.

1

IV.

1

1v.3

1v.3

1v.6

1v.8

IV.

11

v.

1

v.l

v.2

v.4

V.8

v.9

v.ll

v.

12

V.

13

V.

13

V. 14

V.

14

V.16

VI.

1

v1.1

VI,

1

v1.2

v1.2

v1.2

v1.2

v1.3

3. Inleiding

Eén van de moeilijkheden in de rotordynamica is de gyroscopie van de rotorschijven. Dit

uit

zich

in een dempingsmatrix die afhankelijk is van de rotorhoeksnelheid, zodat de

systeemeigenschappen wijzigen met het rotortoerental.

Het programma RODY

is geschreven in PC-Matlab en gaat uit van een eindige

elementenmodel van de rotor. Echter, voor het

gebruik

van

RODY

in de

oude

v ~ r m

i

gmde

kennis

vaiì

de mzthem&che

pïûgïammeeïtaal

PC-lvlaiiab

ñûûuZakeiijk. Om het

programma éénvoudiger toegankelijk te maken is het uitgebreid met een aantal standaard-

functies en van menu’s voorzien. De toegevoegde functies van het RODY-pakket

zijn

globaal:

-

Rotor verdelen in groepen

-

Berekenen eigenfrequenties en -trillingsvormen

-

Tekenen eigentrillingsvormen

-

Berekenen overdrachtsfunctie

-

Tekenen van Bode- en Nyquist-figuur

-

Berekenen gedempte, ongedempte eigenfrequentie en dimensieloze demping

-

Bepalen van whirl

-

Eénvoudige excitatie (bijvoorbeeld onbalans)

Dit heeft tot gevolg dat de rotorberekeningen een grote boekhouding en veel rekentijd op

de computer vergt. Daarom is voor rotoren met meer dan 20 vrijheidsgraden een AT-

386sx

computer of krachtiger sterk aan te raden.

Hoofdstuk

4

behandelt de algemene theorie, waarna deze

in hoofdstuk

5

wordt toegepast

op het gebied van rotordynamica en de verwerking binnen Matlab. Paragraaf

5.3

behandeld het begrip whirl. In hoofdstuk

6

wordt de werkwijze met het programma-pakket

RODY

I1

doorgenomen, aangevuld met belangrijke Matlab-commando’s. Dit wordt

in

hoofdstuk 7 toegepast a.h.v. een rotorvoorbeeld. De nabeschouwing van hoofdstuk

8

heeft

betrekking op de (0n)mogelijkheden van RODY

11.

In hoofdstuk

9

zijn

de opzet en

beknopte uitleg van de Matlab-programma’s opgenomen.

4.

Algemene theorie

In

dit hoofdstuk wordt het werkelijke systeem omgeschreven naar

een

toestands-

beschrijving. Daarbij wordt het systeem opgedeeld in deelsystemen (elementen), waarna

het totale systeem via de eindige elementenmethode wordt samengesteld.

Er wordt ongeveer dezelfde schrijfwijze als in het programma-pakket

RODY

gebruikt.

Figuur

4.1:

Globale en lokale assenstelsel

Indien een rechtsdraaiend assenstelsel, figuur (4.1), wordt gedefiniëerd en daarin de

Y-as

als draaiingsas van de rotor wordt gekozen kan vector

9

de translaties

u resp. w en de

rotaties 0 resp.

cp

in

de positieve X-as resp.

Z-as richting bevatten. De dimensies voor

translatie resp. rotatie zijn in termen van

[ml

resp. [rad].

De rotaties 0 resp.

cp

kunnen worden

u&drukt

in partiële afgeleiden van de translaties

naar de coördinaat-assen:

Bij de elementenmethode wordt rotor opgedeeld

in

lineaire elementen. Daarbij wordt

gebruik gemaakt van de eindige elemententheorie. Het systeem wordt rond de

basiselementen opgebouwd en samengehouden met verbindingselementen. De elementen

worden verdeeld

in twee groepen, namelijk

-

Basiselementen

(As,

balk, schijf en massa)

-

Verbindingselementen (Demper en veer)

Voor deze elementen kunnen bewegingsvergelijkingen worden opgesteld. De element-

matrices en -vectoren zijn

in

bijlage

I

voor basiselementen en

in

bijlage

I1

voor

verbindingselementen beschreven. De bewegingsvergelijkingen (zie literatuur 3) van een

element

zijn

in

het algemeen:

Me

-

4,

+

De

-

q e +

Ke

- -

4e

=

f,

(4.2)

waarin:

9,

f

a

Me

De

K,

lokale vector met

ne vrijheidsgraden

element massamatrix

element dempingsmatrix

element stijfheidsmatrix

lokale

vector

met

exteme

belasting

4.2

Globale beschrijving

De bewegingsvergelijkingen beschrijven het systeem

in een set 2e-orde differentiaal

vergelijkingen. De globale vergelijking is:

M g S

+ D g p

+ K g g = f

(4.3)

wzarln:

9

vector met n

vrijheidspaden

vector met externe belasting

-

f

M,

globale massamatrix

D,

globale dempingsmatrix

$

globale stijfheidsmatrix

4.3 Assemblage

De lagrange beschrijvingen (4.2) van alle elementen worden samengevoegd tot een globale

beschrijving

(4.3)

door een zogenaamde assemblage. De basiselementen worden met de

eindige elementen methode aan elkaar of m.b.v. verbindingselementen gekoppeld.

De afzonderlijke elementmatrices worden geassembleerd door de volgende formules:

Mg = M e i(.) j(b) n-1 a=1 b=l K g =

c

Ke (u) j(b) n = l u=l b=l

(4.4)

(44.6)

waarin:

N

aantal basiselementen

Ne

i(a)

transformatie voor rijen

j(b) transformatie voor kolommen

aantal vrijheidsgraden van het betreffende element

Kolommen i en

j

transformeren de lokale

in globale vrijheidsgraden en zijn afhankelijk

van de keuze van de globale vrijheidsgraden en het basiselement.

De Lagrange bewegingsbeschrijving (4.3)

wordt getransformeerd naar de

toestandsbeschrijving, dit

is een algemene beschrijvingsmethode voor lineaire systemen

waarbij geen afbreuk aan de algemeenheid wordt gedaan (literatuur 4). Deze beschrijving

bestaat uit een toestandsvergelijking (4.7) en een uitgangsvergelijking

(4.8).

X = A x + B u

-

-

-

(4.7)

vraarin:

5

E

A

B

C

D

U

-

tcestandsvector met

E,

elementen

ingangsvector met

nu

elementen

uitgangsvector met

ny

elementen

systeemmatrix (n,

,

n$

ingangsmatrix (n,

,

n,,)

uitgangsmatrix (ny

,

n$

doorverbindingsmatrix (ny

,

n,,)

De bewegingsvergelijkingen

(4.3)

kunnen worden getransformeerd nadat de

toestandsvector bekend is. Een geschikte keuze voor deze toestandsvector:

z

=

k]

(4.9)

Schrijf

(4.3)

m.b.v. de toestandsvector om naar een set le-orde differentiaal vergelijkingen,

waaruit volgt:

Mi'

(M,

4

+

D,

4

+

K,

g)

=

Mi1

f

(4.10)

(4.11)

Indien voor de ingangsvector

combinatie met toestandsvergelijking

(4.7):

de externe belasting

f

wordt genomen volgt uit (4.11)

in

Na vergelijking van (4.7) met (4.12) volgt voor matrices

A

en

B:

(4.12)

(4.13)

Binnen RODY kan de rotor met meerdere externe belastingen worden doorgerekend, deze

belastingen worden na elkaar opgelegd. Daardoor is het aantal ingangen nu =

1.

Omdat de eerste

n

elementen van toestandsvector de vector

9

bevatten wordt de uitgang

y

van de uitgangsvergelijking geschreven in de vorm:

o

o o

o

1

..

o

o

1..

Y

=

I:]

=

[u:]

=

[

1

o

o

o o

..

o

o

...

;]i

+

[o

:::

;]u

(4.14)

In deze uitgangsvergelijking zijn als voorbeeld vrijheidsgraden qs en ql beschouwd. De rij

j

van matrix

C4

bevat alleen nullen, behalve één 1 op de positie i, waarbij de

ie-vrijheidsgraad als de je-uitgang wordt bekeken.

4.5

Eigenwaarde-probleem

Als

van de toestandsvergelijking (4.7) het homogene deel wordt beschouwd, dus =

0,

kunnen de eigenwaarden en eigenvectoren van het systeem worden bepaald. Deze

gegevens vertellen hoe het systeem zich bij vrije responsie gedraagt. De eigenwaarde resp.

eigenvector geven een eigenhoeksnel heid en dempingsfactor [rad/s] resp. een

eigentrillingsvorm

[m]

en [rad] weer. Dit wordt ook mode genoemd. De gedwongen

responsie worât bepaaid door een combinatie van eigenfiequenties en -trillingsvormen.

De toestandsvergelijking van het ongedwongen systeem is:

X

-

= A X

-

(4.15)

De oplossing van (4.15) is van de vorm:

x(t)

-

=

E

e h r

waarin:

1

eigenvector

h

eigenwaarde

(4.16)

Substitutie van (4.16)

in

(4.17) levert:

-

v

h

ent = A

-

v

ent

(hl

-

A)x ent

=

-

O

(4.17)

Voor niet-triviale oplossingen 3

# Q

moet de matrix A uit vergelijking (4.17) singulier

zijn, zodat moet gelden: det(M

-

A)

=

O.

Hieruit volgt de z.g.n. karakteristieke

eigenvector

3,

waarin k =

1 ,...,

n,.

vergelijkimg9

waarin voor .A.

de

h^^gffP

m x h t

E,

is.

Bij

elke

eigeE%raa:Ue

h,

behoc;r:

een

In

tegenstelling tot de eigenwaarden liggen de eigenvectoren

op

een complexe constante

a

na

vast, omdat

a*>

_evenals

3

de eigenvector kan zijn. Omdat dit niet éénduiding vastligt

wordt de modulus van de eerste n vrijheidsgraden van de eigenvector op

1

genormeerd.

4.6

Systeembeschrijving in frequentiedomein

Wanneer een gedwongen responsie wordt opgelegd kan bij de berekening worden

uitgegaan van de toestandsbeschrijving (4.7) en (4.8). Dit is een beschrijving in het

tijddomein.

Als

de toestandsbeschrijving wordt getransformeerd naar het frequentiedomein

kunnen

zaken

als amp~itildeveïhsudimg

en

fcisevers-Jekuiving meer informatie omtrent het

systeem geven,

De Laplace-getransformeerde van de impulsresponsie van de ingang op een uitgang wordt

de overdrachtsfunctie van het systeem genoemd.

In de toestandsbeschrijving ligt dit

verband vast. Met als beginvoorwaarde voor

g&(tJ =

g

met:

to

= O

De toestandsvergelijking (4.7) wordt Laplace-getransformeerd:

s ~ ( s )

-

2

=

A

Z(S)

+

B

E(s)

X(S)

-

= (SI

-

A)-'

x

-

+

(SI

-

A)-l

B

-

U(S)

waarin:

x(s)

toestandsvector

in

frequentiedomein

-

x

beginvoorwaarde

in

tijddomein

-

U(s) ingangsvector in frequentiedomein

De uitgangsvergelijking op overeenkomstige manier wordt getransformeerd tot:

-

Y($

=

c

Z(S)

+

D

-

U(S)

waarin: I(s) uitgangsvector in frequentiedomein

Substitutie van (4.18)

in

(4.19) levert:

-

Y(s)

=

C

(SI -

A)-'

-

x

+

(C

(SI -

A)-'

B

+

O)

-

U(S)

=

C

(SI -

A)-'

₊

_H(s)

_E(s)

met:

H(s)

=

C

(SI

-

A)-l

+ D

(4.18)

(4.19)

(4.20)

waarin: H(s) overdrachtsmatrix van

Q(s)

naar

I(s) met

(n,

,

n,,)

H(S)~

is de overdrachtsfunctie van de i"-uitgang door de je-ingang. Omdat de matrix

(SI -

A)-l

uit (4.20) dezelfde vorm heeft als in het eigenwaarde-probleem

(4.17)

kan deze

matrix ook worden omgeschreven naar:

(4.21)

De dimensies van de responsie voor translatie- resp. rotatie-vrijheidsgraden

zijn

[m/N]

of

[radmm] voor de amplitude resp. [graden] voor de faseverschuiving.

5.

Theorie binnen Matlab

In de rotordynamica blijven de matrices in de bewegingsvergelijkingen (4.3) niet constant.

Door de gyroscopie van de schijven is de dempingsmatrix

D, = D&w) afhankelijk van de

rotorhoeksnelheid,

in een ander geval kunnen dempings- en stijfheidsmatrices ook door

toerentalafhankelijke lagercoëfficiënten wijzigen. De bewegingsvergelijkingen worden:

Mg

ij

+

Dg(0)

4

+

Kg

4

=

f

(5.1)

Dit heeft tot gevolg dat de eigenfrequenties en daarmee de eigentrillingsvormen en

overdrachtsmatrix ook van de rotorhoeksnelheid afhankelijk zijn.

Binnen RODY wordt bij elke rotorhoeksnelheid de systeemmatrix

A

opnieuw

geassembleerd. Het levert geen beduidende tijdwinst op alleen de dempingsmatrix

D, te

assembleren, omdat de rekentijd van bijvoorbeeld eigenfrequenties veel meer vergt.

5.1 Eigenwaarden

In RODY worden alle eigenwaarden van de systeemmatrix uitgerekend. Hoewel alleen één

partitie van de systeemmatrix, zie (4.13) rechtsonder, afhankelijk is van de

rotorhoeksnelheid veranderen alle eigenwaarden.

Indien we te maken hebben met niet al te sterk gedempte systemen blijken de

eigenwaarden in toegevoegd complexe paren voor te komen.

' k

= &

+ j v k

xk

=

&

- j

v k

voor

k

=

1,

...,

n

waarin:

vk

pk

'gedempte' eigenhoeksnelheid voor de

k"

trillingsvorm (mode)

exponentiële dempingsfactor voor de

k"

mode

Na sortering wordt een selectie uitgevoerd welke eigenwaarden worden opgeslagen. Bij

een dynamische analyse zijn vaak alleen de laagste eigenfrequenties van belang, omdat

bijvoorbeeld een onbalans-excitatie

bij

het opstarten van de rotor het lage frequentiegebied

doorloopt. Het programma gaat ervan uit dat de eigenwaarden niet in toegevoegd

complexe paren voorkomen en worden alle eigenwaarden opgeslagen, dit gebeurt zonder

controle met variabele complexjaar, bij de waarde

1

resp. 2 worden toegevoegde

complexe paren wel resp. niet opgeslagen.

Uit de eigenhoeksnelheid

h k

voor de

k"

mode kunnen naast de 'gedempte' eigenhoek-

snelheid

Yk

nog worden bepaald (literatuur 2):

O O k = Iykl

(5.3)

waarin:

mok

'ongedempte' eigenhoeksnelheid voor k" mode

g k

dimensieloze dempingsfactor voor

k"

mode

5.2

Overdrachtsmatrix

Omdat de overdrachtsmatrix afhankelijk

is van de rotorhoeksnelheid is een linearisatie

rond de gemiddelde rotorhoeksnelheid uitgevoerd. Het frequentiegebied wordt opgedeeld

in

intervallen waarbij voor elk interval de systeemmatrix wordt geassembleerd. Voor het

interval worden de systeemeigenschappen nagenoeg constant verondersteld (linearisatie).

De eigenfrequenties en -trillingsvormen kunnen ook worden gebruikt voor het berekenen

van de overdrachtsmatrix, hiervan is geen gebruik gemaakt omdat:

-

De responsie vaak later en bij een andere rotorfrequentie wordt berekend.

-

Alle eigenfrequenties en -trillingsvormen moeten worden opgeslagen.

-

Responsie-berekening vergt minder rekentijd dan voor berekenen van eigentrillings-

vormen nodig is.

De eigentrillingsvormen geven de knooppuntsverplaatsing van de rotor weer en

beschrijven de whirl. De baan van het knooppunt, waarbij

=

O

wordt gehouden, wordt

getekend door oplossing (4.16) en toestandsvector (4.9) alsvolgt om te schrijven:

waarin:

Yi

reële deel vector

imaginaire deel vector

De whirl geeft de draairichting van het knooppunt aan t.o.v. de draairichting van de rotor.

De draairichting kan gelijkgericht resp. tegengesteld

zijn, dit wordt 'forward whirling'

(FW)

resp. 'backward whirling'

(BW)

genoemd. Uit de toestandsvector

en de oplossing

wordt een algorithme ontwikkeld voor de bepaling van de whirl uit de translaties. De

draairichting van het knooppunt kan op elk willekeurig tijdstip worden bepaald door het

uitwendig produkt van verplaatsing en snelheid te bepalen (literatuur

1).

@*4

,

y)

=

@*A

4

9

y)

waarin:

9

vector van knooppunt

e

eenheidsvector in +Y-richting

-Y

Figuur

5.1:

Forward whirling

(FIV)

(5.6)

Figuur

5.2:

backward whiriing

(SW)

Zoals uit figuren

(5.1)

en (5.2) blijkt kan de draairichting worden opgemaakt

uit de

translaties

qi,

q2

en hun afgeleiden

qi,

q2.

De whirl wordt bepaald, door substitutie van qi

en

q2

in

(5.6)

en het uitwendig produkt op t=O te bepalen. Indien deze vector in de

+Z-richting wijst is de draairichting gelijkgericht, forward whirling.

@*a

,

y)l,,

=

reëeZ(ql) reëel@

qJ -

reëel@ ql) reëeZ(qJ

_(5.7)

Een getal kleiner resp. groter dan

nul

geeft forward resp. backward whirling aan. Als het

inwendig produkt precies

nul

is de baan van het knooppunt een rechte.

6. Gebruik

van

RODY I1

In

dit hoofdstuk wordt puntsgewijs aangehaald wat de mogelijkheden van het programma

zijn en een handleiding gegeven

van het programma-pakket RODY

11.

In dit hoofdstuk

wordt een algorithme gegeven hoe het kan worden gebruikt, waarna

in hoofdstuk

7

een

voorbeeld wordt behandeld. Met de beperkingen uit hoofdstuk

8

moet rekening worden

gehouden.

6.

i

Mogeiijkheden

Bij de elementenverdeling van een rotor kan van onderstaande mogelijkheden worden

gebruik gemaakt:

-

Groepen, elementen en vrijheidsgraden kunnen in willekeurige volgorde worden

opgegeven, hoewel wel rekening met de lokale vrijheidsgraden binnen de elementen

moet worden gehouden, zodat de lokale vrijheidsgraden elkaar niet gedeeltelijk

overlappen

-

Omdat voor de basiselementen de lengte wordt opgegeven

ligt

de Y-coördinaat van de

vrijheidsgraad vast

-

Door deze opzet kan een (lokale) elementenverfijning en dergelijke éénvoudig worden

doorgevoerd, zonder dat alle vrij heidsgraden en Y-coördinaten opnieuw moeten worden

bepaald

-

Indien het maximale vrijheidsgraadnummer groter is dan het aantal vrijheidsgraden

wordt met een grotere systeemmatrix gerekend, zodat de berekening onnodig veel tijd

vergt. Dit wordt voorkomen door alle vrijheidsgraden 'te bezetten'

6.2

Algorithme voor het gebruik

Het algorithme, zoals ik het voorstel, ziet er alsvolgt uit:

3)

4)

7)

9)

Maak m.b.v. tekeningen en meetgegevens een beschrijving van de rotor

in

elementen, door knooppunten en vrij heidsgraden te kiezen.

De bepaalde gegevens kunnen worden ingevoerd

in een invoer-file, daarbij kan

bijvoorbeeld

SHELL.m

als uitgangspunt worden genomen. De naam van deze file

De berekening start, wat gebeurt met

RODY.

Er verschijnt een hoofdmenu.

Als

de rotor voor het eerst wordt doorgerekend moet de invoer-file worden

doorlopen (optie 2), als er reeds gegevens bekend zijn kunnen deze uit een uitvoer-

file worden gelezen (optie 3).

De elementen worden gesorteerd en geassembleerd (optie

10).

Het gebruik van vrijheidsgraden wordt gecontroleerd (optie

ll),

tevens wordt een

overzicht van elementen met elementgegevens gepresenteerd

Nu is het verstandig, de berekende startgegevens op te slaan

in een uitvoer-file

(optie

4).

Dit kan ook na elke tijdrovende berekening worden gedaan.

De eigenfi-equenties en e.v.t. -trillingsvormen worden berekend (optie 20) nadat een

aantal gegevens

zijn

ingevoerd, daarom verschijnt een ander menu.

Nadat de eigenfrequenties zijn berekend kunnen deze worden getekend (optie

21)

in

een Campbell-diagram. Indien eigentrillingsvormen bekend zijn kan een animatie

van 'SHAFTS' elementen en kunnen ook de eigentrillingsvormen worden getekend.

wordt

door

het programma

ook elders

gebruih.

10)

Een excitatie-responsie wordt berekend (optie

30)

nadat een aantal gegevens, zoals

excitatie, zijn ingevoerd. De invoer gebeurt in een ander menu.

11)

Na de berekening kan de responsie worden afgebeeld

in een Bode-

of

Nyquist-

diagram (optie

31).

12)

Matlab-commando’s (optie

O)

kunnen ook worden ingevoerd, zie paragraaf 6.3.

13)

De figuren door het programma gemaakt kunnen onder dezelfde naam als van de

invoer-file worden opgeslagen (optie

12).

14)

Het pïoqramma-pakket

RODY

wordt veïlaten met optie 99.

De punten

(1)

t/m

(3) zijn noodzakelijk voor de dynamische analyse. Punt

(99)

sluit het

programma af, alle variabelen blijven dan in het geheugen. Gebleken is dat sommige

versies van Matlab fouten vertonen (zie hoofdstuk

9,

nabeschouwing), dit is extra reden de

berekende gegevens vaker op te slaan (met optie

4).

6.3

Invoeren van gegevens

In veel gevallen kan data-invoer worden overgeslagen als de waarde bekend is of een

invoer-loop moet worden afgesloten door op RETURN te drukken, dit is natuurlijk niet

mogelijk als invoer gewenst is.

In Matlab kunnen variabelen op een aantal manieren worden ingevoerd:

-

Eén waarde: De waarde kan normaal worden ingevoerd, waarna een RETURN

-

Een rij

(1):

Meerdere waarden worden per rij ingevoerd, dit gaat alsvolgt. [waarde,

waarde,

...

waarde,], waarna een RETURN

-

Een rij (2): Ais de waarden met een vast patroon op- of aflopen kan de invoer alsvolgt

-

worden opgegeven: [waarde,

:

increment

:

waarde,],

waarna een

RETURN

6.4

Functies van Matlab

De handleiding van Matlab (literatuur

5 )

behandeld alle commando’s van deze

mathematische programmeertaal. Hieronder staan de belangrijkste genoemd.

Vanuit Matlab moeten alle DOS-commando’s vooraf worden gegaan door een

uitroepteken, zodat bijvoorbeeld inhoudsopgave van de directory wordt opgevraagd met

!DIR

/w

Alle tekst die op het scherm verschijnt, zonder invoer, kan worden opgeslagen in een

bestand onder de naam ’filenaam’ door het commando diary filenaam Dit wordt

ongedaan gemaakt door diary off

Met optie

12

van het hoofdmenu worden alle plaatjes toegevoegd aan de z.g.n. META-

file genaamd ’filenaam.MET’. Dit bestand wordt na het afwerken van het programma

doorgenomen met het commando

GPP

filenaam.MET /dvga

/p

of worden

afgedrukt op een HP laserprinter met het commando

GPP

filenaam.MET /djet

/fprn of

via postscipt met GPP filenaam.MET /dps /fprn

pagina

16

Opmerking:

de diary-file en de META-file worden aangevuld, waardoor de bestanden

een grote omvang kunnen aannemen. Het is daarom verstandig ze

regelmatig te hernoemen

of

te verwijderen.

pagina

17

7.

Voorbeeld

Een kleine rotor voor hoge toerentallen moet aan een dynamische analyse worden

onderworpen. Het rotorfrequentiegebied van O

tlm

2000

Hz

wordt volgens de opzet

uit het

vorige hoofdstuk nader bekeken. De lagering is éénvoudig gemodelleerd. In bijlage

I11

is

de invoer-file

V0ORBEEL.m

opgenomen.

De rotor van figuur (7.1) moet

in eindige elementen worden opgedeeld. Voor een

dynamische analyse moet elk knooppunt

4

vrijheidsgraden bezitten, zodat ook de whirling

wordt bepaald.

i

taatslagers

1

Figuur

7.1:

Dwarsdoorsnede van rotor

Om het programma-pakket

RODY

I1

te kunnen gebruiken worden de afmetingen en

eigenschappen van rotor

en

lagering omgezet naar een eindige elementenmodel. De

topologie staat in figuur (7.2) afgebeeld.

Iii

r

r5

r9

41 3 41 7 10 41 4 41 8

i""

41 45 42 46 Z

T

ì

412 416 j 4 2 0

txy

Figuur

7.2

Vrijheidsgraden van rotor

De gegevens staan niet

in een tabel, omdat deze

in

bijlage

I11

staan of met RODY (optie

11)

zijn op te vragen. Het aantal vrijheidsgraden is zo klein mogelijk gekozen, hier 20.

De lagering is gemodelleerd met

in

zowel

X-

als Z-richting gelijke stijfheden zonder

demping. Doordat

in

het systeem geen demping zit komen de eigenwaarden

in

complexe

paren voor

(zie paragraaf

5.1),

dit verklaart de parameter complexgaar = 2 in de

invoer- file.

Allereerst worden de laagste 8 eigenfrequenties

en

-trillingsvormen (modes) berekend,

voor 8 discrete rotorfrequenties. Na de berekening op een AT-386sx/16 met PC-Matlab

worden modes

2,

4,

6 en 8

in

een Campbell-diagram (figuur

7.3)

weergegeven.

Campbell-diagram van: Voorbeeld

1300 / t

M

IOoo

t

8

' i

FW ...

j

... ~

r

I

... I B W

6

I

/ I I I I ' i

-

/ / I I J I 3 1 i I 200 400 600 8OOx 1000 1200

₌

1400 1600 1800 2ûuO

Figuur

7 3 :

Campbell-diagram van

4

laagste eigenfrequenties

De eigentrillingsvormen (modes)

6

en

8

worden nader toegelicht. Mode

6

vertoond een

backward whirl

(BW),

de eigenfrequentie loopt van 1192 Hz bij een stilstaande rotor tot

1171

Hz bij een rotorfrequentie van 2000 Hz. Mode 8 beschrijft daarentegen een fonvard

whirl

(FW).

De eigenfrequenties en dus ook de trillingsvormen (figuur 7.4) van mode

6

en

8 vallen samen voor de stilstaande rotor, waarbij geen whirling optreedt.

Eigentrillingsvormen van modes

6 en 8 staan

in

figuren (7.5)

en

(7.6).

De onderbroken

lijn

geeft de posities van de knooppunten op tijdstip t=O weer, de hartlijn wordt met een

streep-stiplijn aangegeven en is tevens Y-as.

Eigenfrequentie

= 1193

Hz,

L

I

I

Figuur

7.4:

Modes

6

en

8

bij rotorfrequentie

O

Hz

Eigenfrequentie-=

1180

Hz.

+

I 1

De

whirl kan uit de figuur worden afgeleid, omdat de baan van de knooppunt vanaf

t=O

door een (op het scherm) doorgetrokken

lijn

wordt weergegeven, waarbij het laatste

segment van de baan niet wordt getekend. Uit de lengte van de Y-as kan worden afgelezen

dat de amplitude van de

translatie-vrijheidsgraden

17

en

18

groot

is. De schaalfactor in het

tekenmenu stelt de vergroting van de mode t.o.v. de genormeerde modulus

1

in.

pagina

2Q

t

Eigenfrequentie=

1210

Hz.

I

Groep=

1,

+

Figuur 7.6:

Mode 8 bij rotorfrequentie

1000

Hz, forward whirl

(FW)

7.3

Excitatie-responsie

De rechterkant van de rotor wordt bijvoorbeeld gebruikt om gereedschappen in te spannen.

Als

dit gereedschap niet correct wordt ingespannen voelt de rotor een kracht t.g.v.

onbalans. De vrijheidsgraden

17

en

18 worden daarom geëxciteerd met een onbalanskracht

van

1

N die, natuurlijk, een forward whirl

(FW)

beschrijft. Doordat de willekeurige

beginhoek van de kracht t.o.v. vrijheidsgraad

17

gelijkgesteld is aan O" beschrijft de kracht

op q17 een cosinus en de kracht op qls een

sinus. De excitatie wordt op soortgelijke wijze

als in

(5.5)

voor de baan beschreven is omgezet naar complexe getallen, waaruit volgt:

1 r

UI7 = +I

*\uls

=

-j

u,, =

cos(0

t)

u18

=

sin(o t)

u =

_-

ur

cos(o

t)

+

j

_-

ui

sin(o

t)

-

waarin: u,

reële deel ingangsvector, excitatie

ui

imaginaire deel ingangsvector, excitatie

Bij de responsie worden de verplaatsingen

17

en 18 bekeken, omdat het gereedschap over

nader te bepalen frequentlegebied

aan

bepaalde nauwkeurigheidseisen moet voldoen.

Daartoe moet het frequentiegebied worden opgedeeld in

10

intervallen.

Voor elk interval

wordt de systeemmatrix opnieuw geassembleerd.

Binnen elk interval wordt aangenomen dat de systeemeigenschappen weinig veranderen,

waarbij voor een aantal punten

(10+1)

de overdracht van de excitatie op de geselecteerde

vrijheidsgraden wordt bepaald. Indien de linearisatie een onnauwkeurig resultaat geeft kan

het aantal intervallen alsnog worden vergroot. Het Bode-diagram van vrijheidsgraad

17

is

in figuur

(7.7)

opgenomen.

100- O -100

De onbalans-excitatie is als een doorgetrokken lijn

in

het Campbell-diagram (figuur

7.3)

opgenomen. De onbalans snijdt

in het frequentiegebied

O

t/m

2000

Hz

vier eigen-

frequenties (punten

A

t/m

D), deze frequenties zijn:

A)

780Hz

C)

1175

Hz

D) 1215Hz

B)

865 Hz

,

,;I,,

, !

. f

-

Door de linearisatie verschijnen de resonansiepieken

in

het Bode-diagram (figuur

7.7)

bij

867

en

1220

Hz,

deze frequenties komen overeen met de punten

P

en

Q

van het

Campbell-diagram. Punt

P

resp.

Q

stelt de eigenfrequentie van de FW-mode, berekend

bij

interval V resp.

VI1 voor. Eigenfrequenties uit het Bode-diagram (punten P en

Q)

wijken

iets af van de werkelijke eigenfrequenties uit het Campbell-diagram (punten

B

en D).

Bode

figuur van:

Voorbeeld

3

-

- _- - - o 200 400 600

800

'

i 1 0 0 0 , i 2 0 ~ 140Om16O0= 1800x 2000

I

=I O

Hoekfrequentie rotor

(Hz)

ì d w l j a l :

1 n

=

I I I n on o a on W 1-

.$

s

1

5

2

CIZ

Hoekfrequëntie rotor (HZ)

Ingangnr:

'1,

Uitgang:

17 I t

1

7

200

I

I

Figuur

7.7:

Onbalans-excitatie, gemeten aan vrijheidsgraad

17

8.

Nabeschouwing

In

dit hoofdstuk worden opmerkingen geplaatst over Matlab en het programma-pakket

RODY

11.

A l s

het programma fouten bevat of wordt uitgebreid kan men met de inhoud

van de onderstaande paragrafen hopelijk zijn voordeel doen.

8.1

Gemaakte aannamen

-

Lineair, tijdsonafhankelijk systeem

-

Massaloze, oneindige verbindingselementen

-

Excitatie

is alleen zinvol indien deze betrekking heeft op één groep. Bij een tandwiel-

overbrenging kan bijvoorbeeld tegelijk een moment aan de ingaande en uitgaande

as

worden opgelegd

-

De eigenfrequenties mogen elkaar niet snijden, omdat anders het Campbell-diagram

fout wordt getekend, dit heeft te maken met de sortering

8.2

Overwonnen moeilijkheden

Elementen kunnen

in groepen worden verdeeld, zodat deze groepen onderling m.b.v.

veren en/of dempers verbonden kunnen worden

De verhouding tussen hoeksnelheid van excitatie en hoofdrotor wordt opgegeven met

een Ie-orde polynoom: coexcibtie

=

al*~h,fhtor

+

a,

De elementen kunnen

in

een willekeurige volgorde staan, zodat binnen een groep de

basiselementen van links naar rechts moeten worden gerangschikt

Na het rangschikken kan binnen de groep de Y-coördinaat van elke vrijheidsgraad

worden berekend

Voor een éénvoudige excitatie is het alleen nodig voor de betreffende vrijheidsgraden

de amplitude, beginhoek en whirl op te geven. Waaruit de complexe excitatie wordt

berekend

De whirl wordt bepaald uit de translatie-vrijheidsgraden qi en q2 respectievelijk q5 en q6

van een as-element met

8 vrijheidsgraden

'SHAFT8'

8.3

Wat is nog

niet

zo

handig

-

Snelle uitbreiding van het aantal eigenhoekfrequenties en -trillingsvormen is niet

voorgeprogrammeerd, dit zou handig

zijn om bijvoorbeeld een gebied van

rotorfrequenties achteraf beter te analyseren. Dit kan door de bekende gegevens op te

slaan en de extra gegevens na berekening toe te voegen en daarna e.v.t. te sorteren naar

rotor frequentie

-

Als

in het Campbell-diagram de eigenfrequenties elkaar snijden worden, door de

sortering van eigenwaarden, de modes met elkaar verward

8.4

Nog

te overwinnen moeilijkheden

-

De efficiëntie van de eigenwaarde-routine

in Matlab laat voor dit probleem te wensen

over. Omdat alleen de laagste eigenwaarden van belang

zijn en de eigenwaarden van

pagina

23

het systeem bij de huidige hoekfrequentie van de rotor

in de omgeving liggen van de

eigenwaarden bij lagere hoekfrequentie kan een efficiëntere oplosmethode worden

toegepast

Misschien kan een routine worden geschreven voor de analyse van het opstarten

Voor tandwieloverbrengingen is een groepenverdeling gemaakt, er moeten alleen nog

bewegingsvergelijkingen worden geïmplementeerd, de tandwielen staan immers met

elkaar

in interactie (zie literatuur

1)

Momenten

op

tandwielen moeten ook als excitatie kunnen worden opgegeven

Bij groepen met verschillende hoekfrequenties moet excitatie eenvoudig mogelijk

worden gemaakt door de excitatie per groep af te handelen

soft_eeroutine voor

de

eigenfrequenfies vin

het

Cim@!!!-diiigriim

8.5

Ondervonden problemen met Matlab

-

386-Matlab loopt

(in

1991)

wel eens vast bij inlezen van variabelen uit een MAT-file.

Matlab wordt in zo’n geval afgebroken waarna de DOS-prompt verschijnt

-

Op één of andere manier kunnen de excitatie-responsies soms niet worden getekend na

berekening van eigenwaarden en responsie, het programma RODY wordt afgebroken

-

Computers met kleine geheugens kunnen snel de fout

OUT

of MEMORY geven, vooral

wanneer het aantal vrijheidsgraden groot

is

9.

Programma-pakket

RODY

I1

In

de bijlagen zijn alle (functie-)programma’s opgenomen. Hieronder staat een korte

omschrijving van de inhoud

in volgorde van aanroep, daarbij worden sommige functies

meerdere malen vanuit programmadelen aangeroepen.

1) SHELL

2)

RODY

3)

GLOB-VAR

4) IN-CHECK

5)

SAMEN

6)

DISPLAY

7) DOFCHECK

8) EIG VORM

9)

ABMDK

10) TEKEN

11) TEKENFIG

12)

VERDEEL2

13)

TRANS3D

14) WHIRL

15) PLOTPREP

16)

OVERDR

17) EXCFORCE

18) RODYBODE

Dummy hoofdprogramma wat gecopiëerd kan worden voor het

Hoofd-MENU

waaruit verschiliende functies van

RODY

worden

opgestart

(Her)definieer de globale variabelen

Functie voor het invoeren en controleren van de gegevens

Functie stelt de gebruikte elementen samen in een rij

en

bepaalt de

Y-coördinaat voor elke vrijheidsgraad. Tevens wordt gecontroleerd

of het systeem kan worden geassembleerd

Functie voor het controleren van vrijheidsgraden op dubbel gebruik.

Laat ook een overzicht zien van gebruikte elementen

Functie controleert of de elementen aan elkaar kunnen passen i.v.m.

translatie-, rotatie-vrij heidsgraden en hoeksnelheid

Bepaalt de eigenwaarden en e.v.t. de eigentrillingsvormen

Functie assembleert de massa-, dempings- en stijfieidsmatrix

waaruit de A en

€3 matrix van wordt bepaald

Teken de eigenwaarden of vectoren

in

een plaatje

Functie voor het tekenen van eigentrillingsvormen en bepalen

welke vrij heidsgraden

bij

een bepaalde groep betrokken zijn

Functie voor het berekenen van de sinus- en cosinuswaarden om de

cirkel op te delen

in

stukken, zodat de whirl kan worden getekend

Functie voor het transformeren van 3D-figuur

in

2D-figuur

m.b,v

een transformatiematrix. De tekening wordt isometrisch getekend

Functie voor het bepalen van de whirl a.h.v.

ql

en q2

Functie voor het berekenen van de schaal en het transformeren en

schalen van de figuur naar binnen een rechthoek

Bepaalt de excitatie-overdracht m.b.v. Nyquist

Functie voor het invoeren en berekenen van de excitatiekrachten

Functie voor het tekenen van BODE- of NYQUIST-diagrammen

aanmaken va::

een

nieKW hGGfdprGgr2rnrna

Bij uitbreiden van RODY moet vooral op de volgende punten worden gelet:

Voor nieuwe basiselementen moeten invoer SHELL, assemblage

SAMEN, controle

DISPLAY en de tekenroutine TEKEN worden uitgebreid.

Als

het verbindingselementen

zoals

lagennodellen betreft moeten invoer

SMELL, controle DI§PLAY

wsrdefi

uitgebreid

De nieuwe elementen worden meegerekend indien ze als globale variabelen

GLOB

VAR bekend zijn en ze

bij

het genereren van de systeemmatrix ABMDK

worden- meegenomen. Als bestaande lokale variabelen globaal worden gedefiniëerd

kunnen problemen met de aanroep van subroutines ontstaan

In bijlage VI is het formaat van de gebruikte arrays binnen RODY opgenomen.

10.

Conclusie

Met het programma-pakket

RODY

I1

is voor de ingenieur getracht een dynamische

analyse van rotoren hanteerbaar te maken.

In

het verslag

is het begrip whirl, volgens mij,

genoeg aandacht gegeven. Zoals vermeld moet men bij ingewikkelde excitaties nog

rekenen op trucs met matrices.

SarneEVatteEd -werkt

het

pgrarnrna

gûed

FGû:

ééEV=Udige

rûtG:en, als

CGrnpdterttyd

riet

11.

Literatuur

1) Lalanne, M., Rotordynamics Prediction in Engineering. John Wiley and Sons,

Chichester, 1990

2) Kraker, de,

A.,

Numerieke-Experimentele analyse van Dynamische Systemen.

Technische Universiteit Eindhoven, faculteit Werktuigbouwkunde, Eindhoven, 1990

3)

Campen, van, D.H., Het Dynamische Gedrag van Constructies. Technische

Universiteit Eindhoven, faculteit Werktuigbouwkunde, Eindhoven, 1989

4)

Kok,

J.J.,

Werktuigbouwkundige Regeltechniek

11.

Technische Universiteit

Eindhoven, faculteit Werktuigbouwkunde, Eindhoven, 1990

5 )

Moler,

C.,

PC-Matlab, for MS-DOS Personal Computers. The Mathworks Inc.,

Sherborn, 1987

6) Nelson, H.D., The Dynamics of Rotor-Bearing Systems Using Finite Elements,

Journal

of

Engineering for Industry, 1976

pagina 27

Bijlage

I:

Basiselementen

1.1

As-element met

8

vrijheidsgraden:

'SHAFTS'

i

q3

+

--

I

q5

I

Figuur I.1:As-element met

8

vrijheidsgraden

f,

-

=

v;

fi

m 2 f3

f4

m3 m4IT

Me = Met

+

M,

waarin:

PAL

M&

=

-

420

156

O

o

22L

54

o

O

O

156

-22L

O

O

54

13L

o

- 2 z

4L2

o

o

-13L

-3L2

22L

o

O

4L2

13L

O

O

54

o

O

13L 156

O

O

O

54 -13L O

O

156

22L

O

13L

-3L2 O

o

22L

4L2

-13L

O

o

-3L2

-22L

o

O

-

13L

O

O

-3L

-22L

O

O

4L2

pagina 1.1

36

O

O

3L

-36

O

O

3L

O

36

-3L

O

O

-36

-3L

O

o

-3L

4L2

o

o

3L -L2

o

I

P I

M,

=

-

30L

3L

o

o

4L2

-3L

o o

-L2

-36

O

O

-3L

36

O

O

-3L

-36

3é

0

o

36

3L

0

o

-3L

-L2

o o

3L

4L2

o

I:

o o

-L2

-3L

o o

4L2

EI

K,

=

-

L3

-Pp

0

De

=

30L

6L

O

O

4L2

-6L

O

O

2Lz

-12

O

O

-6L

12

O

O

-6L

'

O

-36

3L

O

O

36

3L

O

36

O

O

3L

-36

O

O

3L

-3L

o o

-4L2

3L

o o

L2

o

-3L

4L2

o

o

3L

-L2

o

O

36

-3L

O

O

-36

-3L

O

-36

O

O

-3L

36

O

O

-3L

-3L

o o

L2

3L

o o

-4E2

o

-3L

-L2

o

o

3L 4L2

o

O

-12

6L

O

O

12

6L

O

O

-6L

2L2

O

O

6L 4L2

O

6L

O

O

2L2

-6L

O

O

4L2

-

12

O

O

6L

-12

O

O

6L

O

12

-6L

O

O

-12

-6L

O

O

-6L

4L2

O

O

6L

2L2

O

I

Het

programma SHAFT8.m is alsvolgt:

function [M2,D2,K2]=shaft8(M,D,K,Sproper,Sto~l,omega,ndim,is)

% SHAFT8 Functie voor het vullen van massa-, dempings-, en stijfheidsmatrix

% voor een as-element met 8 _{vrijheidsgraden. Deze matrices worden} % toegevoegd aan de invoer matrices H, D en K.

% Aanroep: [M2,D2,K2]=shaft8(M,D,K,Sproper,Stoper,Stopol,omega,nd~,~s)

% %

% %

Invoer: M,D,K: resp. globale massa-,dempings- en stijfheidsmatrix, Sproperr eigenschappen as-element:

alle met dimensie (ndim,ndim)

% Sproper(is,l)r lengte (m)

% Sproper(is,Z)r buitendiameter (m)

% Sproper(is,3)r binnendiameter (m)

% Sproper(is,4)r elasticiteitsmodulus (Pa)

% Sproper(is,S)r dichtheid (kg/mA3)

% Stopolr vrijheidsgraden as-element:

% Stopol(is,l)r laagste vrijheidsgraad linkerzijde

% Stopol(is,Z)r laagste vrijheidsgraad rechterzijde

% Stopol(is,3)r elementgroep

% omegat werkelijke hoeksnelheid as (rad/s)

% ndimr aantal vrijheidsgraden totale systeem

% is: elementnummer

%

% uzcvoerr Eizr matrix bi pius bijcirage as-element, aimensie (nciim,ndim)

% D2r matrix D plus bijdrage as-element, dimensie (ndim,ndim)

% K2r matrix K plus bijdrage as-element, dimensie (ndim,ndim)

L=Sproper(is,l);

% Johan Jacobs Philips Natlab jan. 1989 versie 1.0

% Gewijzigd door Cees van Schie

Du=Sproper(is,Z); Di=Sproper(is,d); E=Sproper(is,4); Rho=Sproper ( is, 5 ) ; Q(l)=Stopol(is,l); Q(2)=Q(1)+1; Q(4)=Q(1)+3; Q(5)=Stopol(is,2); Q ( 6) =Q(5 +I; Q(7)=Q(5)+2; Q(8)=Q(5)+3; elementgroep=Stopol(is,3); cl=156; c2=22*L; c3-54; c4=13*L; c5=4*L*L; c6=3 *L*L; mu=Rho*pi*(Du*D~-Di*Di)/4; c7=mu*L/420;

met-zeros (ndim, ndim) ; MBer=MBet;

DBe-MBet ;

KBe=MBet ;

%

% vullen van massa-matrix

% MBet(Q,Q)=c7*[ cl O O c2 c3 O O -c4 i r z I _ - _ 17-4-91 4(3)=~(1)+2;

o

cl -c2

o

o c3 c4

o

O -c2 ~5 O O -c4 -c6 O ~2 O O ~5 ~4 O O -c6 6 3 O 9 CO cl O O -e2

o

c3 -c4

o

o cl c2 o O ~4 -c6 O O c2 c5 O -c4 O O -c6 -c2 O O ~5

1 ;

cl-36; c2=3 *L; c3=4*L*L; c4=L*L; it=mu*(Du*Du+Di*Di)/16; c5=it/(30*L) ; MBer(Q,Q)=c5*[ cl O O o cl -c2

o

-c2 c3 c2 o

o

-cl o

o

o

-cl c2

o

-c2 -c4 c2

o

o

MZ=M+MBet+MBer; %

% vullen van dempingsmatrix

% ip=2*it; c2 -cl

o

o

c2

o

o

-cl -c2

o

o

o c2 -c4

o

c3 -c2 o 0 -co -c2 cl

o

o

-c2

o

o

cl c2 o

o

o

c2 c3 o -e4 -e2 o

o

c3 I ;