Kompendium Der Maschinendynamik
Kompendium Der Maschinendynamik
Kompendium Der Maschinendynamik
Dieter Kraft
Fachbereich Maschinenbau
Fachhochschule München
ii
Dieter
c Kraft. 5. Auflage, WS 1999/2000.
Inhaltsverzeichnis
Vorwort viii
1 Einführung 1
1.1 Modellbildung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Newton-Euler . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Harmonische Bewegung . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Begriffe . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 Beschreibung . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.1 Zeitbereich . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.2 Frequenzbereich . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Physikalische Ersatzsysteme . . . . . . . . . . . . . . . . . . . . . 8
1.4.1 Federkonstanten . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.1.1 Addition . . . . . . . . . . . . . . . . . . . . . . 8
1.4.1.2 Balken und Stäbe . . . . . . . . . . . . . . . . . 9
2 Ein-Freiheitsgrad-Systeme 11
2.1 Ungedämpfte Eigenbewegung . . . . . . . . . . . . . . . . . . . . 11
2.1.1 Lösung der Bewegungsgleichung . . . . . . . . . . . . . . 11
2.1.2 Bewegungsformen . . . . . . . . . . . . . . . . . . . . . . 13
2.1.2.1 Anfangsbedingungen . . . . . . . . . . . . . . . . 13
2.1.2.2 Winkeladdition . . . . . . . . . . . . . . . . . . . 14
2.1.2.3 Phasenportrait . . . . . . . . . . . . . . . . . . . 15
2.1.3 Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.3.1 Manometer-Schwingung . . . . . . . . . . . . . . 15
2.1.3.2 Strukturschwingung . . . . . . . . . . . . . . . . 16
2.1.3.3 Torsionsschwingung . . . . . . . . . . . . . . . . 16
2.2 Gedämpfte Eigenbewegung . . . . . . . . . . . . . . . . . . . . . 17
2.2.1 Lösung der Bewegungsgleichung . . . . . . . . . . . . . . 17
2.2.2 Bewegungsformen . . . . . . . . . . . . . . . . . . . . . . 18
2.2.2.1 Ungedämpfte Schwingung . . . . . . . . . . . . . 18
2.2.2.2 Kritische Dämpfung . . . . . . . . . . . . . . . . 19
2.2.2.3 Unterdämpftes System . . . . . . . . . . . . . . 19
2.2.2.4 Überdämpftes System . . . . . . . . . . . . . . . 22
iv INHALTSVERZEICHNIS
2.2.3 Charakteristika . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.3.1 Stabilität . . . . . . . . . . . . . . . . . . . . . . 23
2.2.3.2 Logarithmisches Dekrement . . . . . . . . . . . . 24
2.2.3.3 Energieverluste . . . . . . . . . . . . . . . . . . . 25
2.2.4 Alternative Dämpfungsmodelle . . . . . . . . . . . . . . . 25
2.2.4.1 Trockene Reibung . . . . . . . . . . . . . . . . . 26
2.3 Harmonische Anregung . . . . . . . . . . . . . . . . . . . . . . . 28
2.3.1 Harmonische Anregung des ungedämpften Systems . . . . 29
2.3.1.1 Resonanz . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1.2 Schwebung . . . . . . . . . . . . . . . . . . . . . 32
2.3.2 Harmonische Anregung des gedämpften Sytems . . . . . . 34
2.3.2.1 Lösung der Bewegungsgleichung . . . . . . . . . 34
2.3.2.2 Diskussion der Lösung . . . . . . . . . . . . . . . 35
2.3.2.2.1 Vergrößerungsfunktion I . . . . . . . . . 35
2.3.2.2.2 Phasenwinkel I . . . . . . . . . . . . . . 37
2.3.2.2.3 Kraftübertragung I . . . . . . . . . . . 38
2.3.2.2.4 Bandbreite . . . . . . . . . . . . . . . . 38
2.3.3 Komplexe Erregerfunktion . . . . . . . . . . . . . . . . . . 40
2.3.3.1 Lösung der Bewegungsgleichung . . . . . . . . . 41
2.3.3.2 Frequenzgang . . . . . . . . . . . . . . . . . . . . 41
2.3.3.3 Frequenzgang und Übertragungsfunktion . . . . 42
2.3.3.4 Geometrische Darstellung . . . . . . . . . . . . . 43
2.3.4 Alternative harmonische Anregungen . . . . . . . . . . . . 43
2.3.4.1 Harmonische Anregung der Basis . . . . . . . . . 44
2.3.4.1.1 Vergrößerungsfunktion II . . . . . . . . 46
2.3.4.1.2 Relativbewegung . . . . . . . . . . . . . 48
2.3.4.1.3 Kraftübertragung II . . . . . . . . . . . 49
2.3.4.2 Rotatorische Unwucht . . . . . . . . . . . . . . . 50
2.3.5 Alternative Dämpfungsmodelle II . . . . . . . . . . . . . . 53
2.3.5.1 Coulomb-Dämpfung . . . . . . . . . . . . . . . . 53
2.3.5.2 Hysterese-Dämpfung . . . . . . . . . . . . . . . . 53
2.3.5.3 Luft-Dämpfung . . . . . . . . . . . . . . . . . . 56
2.3.5.4 Kombinierte Dämpfungen . . . . . . . . . . . . . 57
2.3.6 Strömungs-induzierte Schwingungen . . . . . . . . . . . . 57
2.3.7 Beispiele und Übungen . . . . . . . . . . . . . . . . . . . . 58
2.4 Nicht-harmonische Anregungsfunktionen . . . . . . . . . . . . . . 66
2.4.1 Sprung-Funktion . . . . . . . . . . . . . . . . . . . . . . . 67
2.4.1.1 Sprungfunktion im Zeitbereich . . . . . . . . . . 67
2.4.1.2 Sprungfunktion im Frequenzbereich . . . . . . . 68
2.4.1.3 Anwendung . . . . . . . . . . . . . . . . . . . . . 68
2.4.2 Impuls-Funktion . . . . . . . . . . . . . . . . . . . . . . . 71
2.4.2.1 Impulsfunktion im Zeitbereich . . . . . . . . . . 72
2.4.2.2 Impulsfunktion im Frequenzbereich . . . . . . . 73
2.4.2.3 Anwendung . . . . . . . . . . . . . . . . . . . . . 74
2.4.3 Allgemeine Anregung . . . . . . . . . . . . . . . . . . . . 77
2.4.3.1 Impulsantwort . . . . . . . . . . . . . . . . . . . 77
INHALTSVERZEICHNIS v
3 Zwei-Freiheitsgrad-Systeme 81
3.1 Eigenbewegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.1.1 Ungedämpfte Eigenbewegung . . . . . . . . . . . . . . . . 81
3.1.1.1 Modellbildung . . . . . . . . . . . . . . . . . . . 81
3.1.1.2 Lösung der Matrix-Differentialgleichung . . . . . 82
3.1.1.2.1 Lösungsansatz . . . . . . . . . . . . . . 82
3.1.1.2.2 Eigenfrequenzen . . . . . . . . . . . . . 83
3.1.1.2.3 Eigenformen . . . . . . . . . . . . . . . 84
3.1.1.2.4 Allgemeine Lösung . . . . . . . . . . . . 85
3.1.1.3 Modal-Analyse . . . . . . . . . . . . . . . . . . . 87
3.1.1.3.1 Eigenwerte und Eigenvektoren . . . . . 87
3.1.1.3.2 Transformation auf Diagonalform . . . 90
3.1.1.3.3 Entkopplung der Systemgleichungen . . 91
3.1.1.4 Beispiele . . . . . . . . . . . . . . . . . . . . . . 95
3.1.1.4.1 Werkzeugmaschine . . . . . . . . . . . . 96
3.1.1.4.2 Doppelpendel . . . . . . . . . . . . . . . 99
3.1.2 Gedämpfte Eigenbewegung . . . . . . . . . . . . . . . . . 102
3.1.2.1 Modell-Gleichungen . . . . . . . . . . . . . . . . 102
3.1.2.2 Modal-Analyse . . . . . . . . . . . . . . . . . . . 103
3.2 Erzwungene Bewegung . . . . . . . . . . . . . . . . . . . . . . . . 104
3.2.1 Bewegungsgleichungen . . . . . . . . . . . . . . . . . . . . 104
3.2.2 Modal-Analyse . . . . . . . . . . . . . . . . . . . . . . . . 104
3.2.3 Drei Integrale . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.2.3.1 Harmonische Anregung des gedämpften Systems 105
3.2.3.2 Sprungförmige Anregung des gedämpften Systems105
3.2.3.3 Harmonische Anregung des ungedämpften Sys-
tems . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.2.3.4 Lösungen der Integrale . . . . . . . . . . . . . . 105
3.2.4 Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.2.4.1 KFZ-Aufhängung . . . . . . . . . . . . . . . . . 106
3.2.4.2 KFZ-Nickbewegung . . . . . . . . . . . . . . . . 109
4 Viel-Freiheitsgrad-Systeme 111
Anhang 154
A Integral-Transformationen 155
A.1 Fourier-Transformation . . . . . . . . . . . . . . . . . . . . . . . 155
A.1.1 Fourier-Reihen . . . . . . . . . . . . . . . . . . . . . . . . 155
A.1.2 Diskrete Fourier-Transformation . . . . . . . . . . . . . . 156
A.1.3 Fenster-Funktionen . . . . . . . . . . . . . . . . . . . . . . 156
A.1.4 Eigenschaften und Berechnung . . . . . . . . . . . . . . . 156
A.1.5 Beispiel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
A.2 Laplace-Transformation . . . . . . . . . . . . . . . . . . . . . . . 158
A.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
A.2.2 Eigenschaften . . . . . . . . . . . . . . . . . . . . . . . . . 161
A.3 Partialbruchzerlegung . . . . . . . . . . . . . . . . . . . . . . . . 165
INHALTSVERZEICHNIS vii
B Matrixalgebra 171
B.1 Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
B.2 Algebraische Operationen . . . . . . . . . . . . . . . . . . . . . . 172
B.3 Iterative Operationen . . . . . . . . . . . . . . . . . . . . . . . . 173
C Linearisierung 179
C.1 Taylor-Entwicklung . . . . . . . . . . . . . . . . . . . . . . . . . . 179
C.2 Standard-Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
C.3 Beispiele und Übungen . . . . . . . . . . . . . . . . . . . . . . . . 181
Bibliographie 209
viii INHALTSVERZEICHNIS
Vorwort
2. Das Schreiben dieses Manuskripts hat mir sehr viel Freude gemacht,2
3. Das Manuskript enthält wahrscheinlich viele Fehler. Helfen Sie mir durch Meine Sprechstunde:
Rückkopplung: Nennen Sie Fehler, machen Sie Verbesserungsvorschläge, Di 13–14 Uhr,
fragen Sie bei Unklarheiten, etc.3 Raum W-1.07,
Tel. 1265-1177.
4. Der jeweils letzte Stand des Manuskripts ist als Postscript-File unter http://www.fm.fh-
muenchen.de:8080/˜dkraft verfügbar.
5. Dank sei gesagt den Vielen, die nach den Grundsätzen der Free Software
Foundation das Gnu-Projekt vorantreiben, das Betriebssystem Linux pro-
movieren, das geniale Textsatzsystem TEX geschrieben haben, das einfache
CAD-System Xfig managen, u.v.a.m. Stellvertretend für viele hier nur
einige Namen: Richard Stallman, Linus Torvalds, Donald Knuth, Brian
Smith, Tomas Rokicki, Karl Berry, Thomas Esser.
6. Ein dringender Hinweis:Vollziehen Sie die Beispiele nach; fragen Sie, wenn !
ich Zwischenschritte zu eilig ausgelassen habe, rechnen Sie die Aufgaben
mit Papier und Bleistift durch,4 reines Durchlesen ist überflüssige Zeitver-
schwendung.
Einführung
Die Anzahl der unabhängigen Koordinaten, die zur Beschreibung eines Systems
genügt, wird als der Freiheitsgrad des Systems bezeichnet. So ist ein inertial auf- !
gehängtes, ebenes Pendel ein System mit einem Freiheitsgrad, da der momen-
tane Pendelwinkel relativ zu einer Referenzlage die Pendelbewegung eindeutig
beschreibt. Ein analog geartetes Doppelpendel ist ein System mit zwei Freiheits-
graden, eine schwingende Saite ein solches mit unendlich vielen Freiheitsgraden.
Ein räumlich ausgedehnter starrer Körper hat sechs Freiheitsgrade, drei der
Translation und drei der Rotation.
Diese Einführung beschränkt sich auf lineare oder abschnittsweise lineare Bewe-
gung1 , bei denen das Superpositionsprinzip gilt, und für die die abgeschlossene
Theorie linearer Differentialgleichungen zur Verfügung steht.
1.1 Modellbildung
1.1.1 Newton-Euler
(unabhängig von Gottfried Wilhem Leibniz (1646–1716)) die Differential- und Integralrech-
nung. Bahnbrechende Beiträge zur Entwicklung der klassischen Mechanik.
3 Leonhard Euler (1707–1783), schweizer Mathematiker, der in St. Petersburg und Berlin
wirkte. Beiträge auf fast allen Gebieten der Mathematik und Physik: Differential- und Dif-
ferenzengleichungen, komplexe Analysis, spezielle Funktionen, Reihen, Variationsrechnung,
Mechanik, Fluidmechanik.
4 In Gl. (1.1.1) bedeutet m die Masse des Systems, und in Gl. (1.1.2) ein von außen auf das
Modell(Differential)gleichung
Sie wird noch durch Anfangs- und/oder Randwerte, z.B. x(0) = x0 und ẋ(0) = ẋ0
festgelegt.
Die Differentialgleichung ist homogen, aber nichtlinear. Eine Entwicklung der
Sinus-Funktion sin θ = θ − θ3 /3 + und Abbruch der Entwicklung nach dem
ersten Term führt auf folgende lineare Differentialgleichung, s.a. Anhang C,
g
θ̈(t) + θ(t) = 0.
l
1.1.2 Lagrange
Die Modellbildung nach Newton und Euler hat den Nachteil, daß Sie bei Syste-
men mit mehreren Freiheitsgraden die Teilsysteme freischneiden und ihre Reak-
tionskräfte bestimmen müssen, die als innere Kräfte und Momente dann in den
Bewegungsgleichungen garnicht explizit auftreten. Dieser Umstand haftet der
Methode nach Lagrange5 nicht an. Es wird die kinetische Energie T = T (ẋ, x)
und die potentielle Energie U = U(x) und die Lagrange-Funktion L = L(ẋ, x)als
deren Differenz gebildet:
und dann erhält man die Bewegungsgleichungen aus den Gleichungen 2. Art von
Lagrange:
d ∂L ∂L
− = f(t), i = 1, . . . , n, (1.1.3)
dt ∂ẋi (t) ∂xi (t)
mit n der Anzahl der beteiligten Teilsysteme.
Beispiel 2 Für das in Abb. 1.1 dargestellte mechanische System mit den beiden
Freiheitsgraden x1 , der Position der Masse m, und x2 , der Winkelauslenkung
des Schwungrades, seien die Bewegungsgleichungen nach der Methode von La-
grange herzuleiten.
Die kinetische Energie des Systems ist
1 1
T= mẋ21 + Jẋ22 ,
2 2
die potentielle Energie
1 2 1
U= kx + k(rx2 − x1 )2 ,
2 1 2
5 Joseph Louis Lagrange (1736–1813), lebte und wirkte in Turin, Berlin und Paris. Bedeu-
x1 x2
11
00
00
11
00
11 k k
00
11 m
00
11
00
11 r
00
11 11111111111
00000000000 J
111111
000000
00
11
00
11 00000000000
11111111111 00
11
000000
111111
m(t)
d ∂L ∂L ∂F
− + = f(t), i = 1, . . . , n. (1.1.4)
dt ∂ẋi (t) ∂xi (t) ∂ẋi (t)
x(t)
A sin ω t
A A
ωt
ωt
2π
In Abb. 1.3 ist eine abklingende (gedämpfte) harmonische Bewegung einer Masse
m, die über eine elastische Feder und einen viskosen Dämpfer inertial fest auf-
gehängt ist und die eine Lichtquelle trägt, auf ein lichtempfindliches Papier
projiziert, das an der Anordung mit konstanter Geschwindigkeit vorbeiläuft[48].
11111
00000
k c
x(t)
A
m
t
1.2.1 Begriffe
Die harmonische Bewegung wird nach Abb. 1.3 durch die Auslenkung x(t)
der Masse beschrieben, wobei A die Amplitude und ω die Kreisfrequenz, ge-
messen in rad/s, ist. Die Dauer einer Schwingungsperiode ist T , und der Zusam-
6 Einführung
1.2.2 Beschreibung
Wir sind nicht nur an der Auslenkung x sondern auch an der Geschwindigkeit ẋ
der Masse interessiert, und wir beobachten jeweils eine Phasenverschiebung von
π/2 zwischen diesen Größen. Zwischen Auslenkung und Beschleunigung besteht
folgender Zusammenhang:6
ẍ = −ω2 x, (1.2.5)
die Beschleunigung ist der Auslenkung proportional und zeigt in Richtung auf
den Ursprung der Bewegung, die Gleichgewichtslage.
Mit Hilfe der Euler-Identität
z(t) = Aejωt
= A cos ωt + jA sin ωt (1.2.7)
= x + jy
ausgedrückt werden.
Aufgabe 3 Zeigen Sie, daß die Beziehung (1.2.7) die Differentialgleichung (1.2.5)
erfüllt.
2
Die zu z konjugiert komplexe Größe ist
z£ = Ae−jωt ,
1.3 Definitionen
1.3.1 Zeitbereich
1.3.2 Frequenzbereich
Dezibel: Das Dezibel8 ist ein logarithmisches Maß des Verhältnisses zweier Leistungen:
P1
db = 10 log10
P2
oder mit einem Ergebnis aus Beispiel 6, das besagt, daß die Leistung proportio-
nal der Amplitude ist,
x21 x1
db = 10 log10 2
= 20 log10 . (1.3.4)
x2 x2
Oktave: Die Oktave ist eine Bezeichnung für die Bandbreite einer Schwingung. Wenn
die obere Grenzfrequenz doppelt so groß ist wie die untere, dann deckt das
Frequenzspektrum eine Oktave ab.
Beispiel: Der Frequenzbereich erstrecke sich von 250 Hz bis 500 Hz, dann ist die
Bandbreite 500 − 250 = 250 Hz oder eine Oktave.
f1 f1
= 2n oder n = log2 (1.3.5)
f2 f2
gilt.
1.4.1 Federkonstanten
1.4.1.1 Addition
Zwei oder mehrere (n) parallel geschaltete Federn mit den Steifigkeiten ki haben
als äquivalente Federkonstante k die Summe der ki:
n
k= ki , (1.4.1)
i=1
8 Das
Dezibel ist nach G. A. Bell, dem Gründer der Bell Company, später AT& T, benannt.
9 Dieser
Abschnitt wird im Laufe des Semesters ergänzt. Einige Formeln sind schon im Text
verwendet, ohne hier zu erscheinen, s.z.B. Gl. (2.1.12).
1.4 Physikalische Ersatzsysteme 9
1 1
n
= , (1.4.2)
k ki
i=1
oder n
ki
k = i=1
n . (1.4.3)
i=1 ki
Ein-Freiheitsgrad-Systeme
111111
000000
k
entspannte k∆
∆ k( ∆ +x) stat.
Feder
m x Gleichgew.
m
m
. ..
W x x
W
ms2 + k = 0, (2.1.3)
k
die sogenannte charakteristische Gleichung, aus der
k
s1,2 = − = j (2.1.4)
m m
folgt. Der Ausdruck6 k
ωn = (2.1.5)
m
! wird als die natürliche (ungedämpfte) Frequenz des Systems bezeichnet. Diese
wird experimentell wie folgt ermittelt:
Aus der statischen Auslenkung ∆ der Feder, s.a. Abb. 2.1, und der Masse m
wird zunächst die Federkonstante k bestimmt:
mg
k=
∆
→ Experiment und diese dann in Gl. (2.1.5) eingesetzt:
g
ωn = .
∆
Damit erhält man auch die Periode Tn = 1/fn der Eigenschwingung:
∆
Tn = 2π .
g
5 Die Lösung für C = 0 wird triviale Lösung genannt und interessiert uns nicht.
6 Oft findet man in der Literatur auch ω0 statt ωn .
2.1 Ungedämpfte Eigenbewegung 13
Die beiden Lösungen (2.1.4) der charakteristischen Gleichung (2.1.3) heißen cha-
rakteristische oder Eigenwerte des Systems. Mit ihnen ist die allgemeine Lösung
der Dgl. (2.1.1)
x(t) = C1 ejωn t + C2 e−jωn t , (2.1.6)
mit aus den Anfangswerten zu bestimmenden konjugiert komplexen Konstanten
C1 und C2 , oder unter Verwendung der Euler-Identität e¦jα = cos α j sin α
x(t) = A1 cos ωn t + A2 sin ωn t, (2.1.7)
wobei A1 und A2 wiederum aus den Anfangswerten wie folgt zu bestimmen
sind:
x(0) = A1 = x0 , ẋ(0) = ωn A2 = ẋ0
und damit
ẋ0
A1 = x0 , A2 = .
ωn
In die Lösung (2.1.7) eingesetzt wird
ẋ0
x(t) = x0 cos ωn t + sin ωn t (2.1.8)
ωn
als allgemeine Lösung erhalten.
Zwischen den Konstanten Ai und Ci besteht folgende Beziehung:
A1 = C1 + C2 , und A2 = j(C1 − C2 ).
Aufgabe 5 Leiten Sie diese Beziehung her; und schreiben Sie sich Ci = f(Ai )
in den Margin.
2
2.1.2 Bewegungsformen
Die allgemeine Lösung (2.1.8) stellt eine harmonische Schwingung dar, deswegen
ist das durch Gl. (2.1.1) beschriebene System ein harmonischer Oszillator.
2.1.2.1 Anfangsbedingungen
Gl. (2.1.8) beschreibt die Bewegung der Masse m um ihre statische Gleichge-
wichtslage. Eine (Eigen-)Bewegung kommt nur zustande, wenn mindestens eine
der beiden Anfangsbedingungen von Null verschieden sind.
Für ẋ0 = 0 ergibt sich
x(t) = x0 cos ωn t
und für x0 = 0
ẋ0
x(t) = sin ωn t.
ωn
Der Einfluß der Anfangswerte auf die Lösung x(t) ist in Abb. 2.3 für ωn = 1
dargestellt. Sie erkennen deutlich die unterschiedlichen Amplituden der Schwin-
gungen.
14 Ein-Freiheitsgrad-Systeme
1.5
0.5
←x =0,xp =1/2
0 0
←x0=1,xp0=0
x(t)
−0.5
←x =−1/2,xp =1
0 0
−1
−1.5
0 1 2 3 4 5 6 7
t[s]
2.1.2.2 Winkeladdition
2.1.2.3 Phasenportrait
Ein- bzw. zweimalige Differentiation der Lösung (2.1.9) ergeben die Geschwin-
digkeit
x(t)
cos (ωn t − φ) =
A
und mit
ẋ(t)
sin (ωn t − φ) = −
Aωn
das Phasendiagramm einer Ellipse
x2 ẋ2
+ = 1.
A2 (Aωn )2
2.1.3 Beispiele
2.1.3.1 Manometer-Schwingung
mẍ = fx ,
d.h.
LAρẍ = −2Axρg
oder vereinfacht
2g
ẍ + x=0
L
16 Ein-Freiheitsgrad-Systeme
2.1.3.2 Strukturschwingung
2.1.3.3 Torsionsschwingung
Ein Rad hängt an einer eingespannten Welle mit Durchmesser d und Länge l.
Das Rad wird aus seiner statischen Gleichgewichtslage ausgelenkt und vollführt
in t1 Sekunden n Schwingungen. Bestimmen Sie das polare Trägheitsmoment J
des Rades.
Lösung: Die Bewegungsgleichung ist
Jθ̈ + kt θ = 0.
Die Eigenkreisfrequenz ωn berechnet sich aus der Anzahl der Schwingungen in der
gegebenen Zeit:
n
ωn = 2π .
t1
Die Torsionssteifigkeit der Welle ist
Gπd4
kt = (2.1.12)
32l
mit G dem Schermodul des Wellenmaterials.
kt Gd4 t21
J= 2
= .
ωn 128πln2
2.2 Gedämpfte Eigenbewegung 17
lautet.
Zur Lösung der Bewegungsgleichung (2.2.1) wird wieder der Ansatz (2.1.2)
x(t) = Cest herangezogen. Die charakteristische Gleichung ergibt sich hier nach
(ms2 + cs + k)Cest = 0 zu
ms2 + cs + k = 0, (2.2.2)
c
c 2 k
s1,2 =− − . (2.2.3)
2m 2m m
und damit ist die allgemeine Lösung der Bewegungsgleichung (2.2.1) !
− 2m
c
+
( 2m
c
)2 − m
k
t
− 2m
c
−
( 2m
c
)2 − m
k
t
x(t) = C1 e + C2 e , (2.2.4)
geschrieben werden, wobei nun die drei Parameter m, c, k durch nur zwei, ωn , ζ,
ersetzt sind.
Das Dämpfungsverhältnis wird in die Beziehung (2.2.3) eingesetzt:
s1,2 = −ζ ζ 2 − 1 ωn , (2.2.7)
und damit wird als Lösung (2.2.4)
−ζ+ ζ2 −1 ωn t −ζ− ζ2 −1 ωn t
x(t) = C1 e + C2 e (2.2.8)
erhalten. Setzt man das Dämpfungsverhältnis ζ und die Eigenkreisfrequenz ωn
in die Systemgleichung (2.2.1) ein, dann erhält man die häufig verwendete Dars-
tellung
ẍ(t) + 2ζωn ẋ(t) + ω2n x(t) = 0 (2.2.9)
und entsprechende Anfangsbedingungen.
2.2.2 Bewegungsformen
Der Fall ζ = 0 führt auf ungedämpfte Schwingungen des Abschnitts (2.1). Hier
soll die Lösung (2.1.8) durch Anwendung der Laplace-Transformation[36], s.a.
Anhang A.2, gewonnen werden. Die Laplace-Transformation der homogenen Dif-
ferentialgleichung mẍ(t) + kx(t) = 0 oder, nach Division durch m,
ẍ(t) + ω2n x(t) = 0
Merke: x(t) c s X(s) lautet
s2 X(s) − sx0 − ẋ0 + ω2n X(s) = 0,
und nach X aufgelöst:
s 1
X(s) = x0 + 2 ẋ0 .
s2 + ω2n s + ω2n
Die inverse Laplace-Transformation dieses Ausdrucks ergibt mit Zeile 23 und 24
der Tab. A.1
ẋ0
x(t) = x0 cos ωn t + sin ωn t,
ωn
wie nach Ergebnis (2.1.8) zu erwarten war.
2.2 Gedämpfte Eigenbewegung 19
Im Falle von ζ = 1, der kritischen Dämpfung, fallen die beiden Wurzeln der
charakteristischen Gleichung (2.2.2) zusammen:
ck
s1,2 = − = −ωn ,
2m
dementsprechend lautet die Systemgleichung
ẍ(t) + 2ωn ẋ(t) + ω2n x(t) = 0.
Im Fall des unterdämpften Systems 0 < ζ < 1 wird die Diskriminante der Lösung
der charakteristischen Gleichung (2.2.7) negativ, deshalb kann diese Lösung auch
als
s1,2 = −ζ j 1 − ζ2 ωn (2.2.11)
geschrieben werden. Damit wird die Lösung (2.2.8) der Systemgleichung
−ζ+j −ζ−j
1−ζ2 ωn t 1−ζ2 ωn t
x(t) = C1 e + C2 e
= e−ζωn t C1 ej 1−ζ2 ωn t
+ C2 e−j 1−ζ2 ωn t
= e−ζωn t (C1 + C2 ) cos 1 − ζ 2 ωn t
+ j(C1 − C2 ) sin 1 − ζ 2 ωn t
= e−ζωn t A1 cos 1 − ζ2 ωn t + A2 sin 1 − ζ 2 ωn t
A1 = x0
und
ẋ0 + ζωn x0 ẋ0 + ζωn x0
A2 = = ,
1 − ζ 2 ωn ωd
wobei, mit
ωd = 1 − ζ 2 ωn (2.2.14)
der Frequenz der gedämpften Schwingung, die allgemeine Lösung
ẋ0 + ζωn x0
−ζωn t
x(t) = e x0 cos ωd t + sin ωd t (2.2.15)
ωd
ẋ 2
0 + ζωn x0
A= A21 + A22 = x20 +
ωd
←x =1,xp =0
0 0
0.8
0.6
←x0=0,xp0=1
0.4
←x =−1/2,xp =0
0 0
x(t)
0.2
−0.2
−0.4
−0.6
0 2 4 6 8 10 12 14 16 18 20
t[s]
x0=-1/2;
x0p=0;
A1=x0;
A2=(x0p+z*wn*x0)/wd;
x3=exp(-z*wn.*t).*(A1.*cos(wd.*t)+A2.*sin(wd.*t));
plot(t,x1,’k-’,t,x2,’k--’,t,x3,’k-.’)
grid
xlabel(’t[s]’)
ylabel(’x(t)’)
gtext(’\leftarrowx_0=1,xp_0=0’)
gtext(’\leftarrowx_0=0,xp_0=1’)
gtext(’\leftarrowx_0=-1/2,xp_0=0’)
Abb. 2.4 zeigt den Einfluß der Dämpfung, ζ {1/5, 1/3, 1/2}, auf das System-
verhalten.
Aufgabe 8 Untersuchen Sie anhand der Abb. 2.4 den Einfluß der Dämpfung auf
die Schwingungsdauer T . Begründen Sie Ihre Erfahrung durch die Theorie.
2
22 Ein-Freiheitsgrad-Systeme
0.8
0.6
0.4
←ζ=1/5
x(t)
0.2
←ζ=1/2
0
−0.2
←ζ=1/3
−0.4
−0.6
0 2 4 6 8 10 12 14 16 18 20
t[s]
reel und unterschiedlich, und es gilt s2 < s1 < 0. Die Lösung der Bewegungs-
gleichung ist
−ζ+
ζ2 −1 ωn t
−ζ−
ζ2 −1 ωn t
x(t) = C1 e + C2 e , (2.2.16)
A := 1 1
−ζ˜ + ζ˜2 − 1 w ˜ −ζ˜ − ζ˜2 − 1 w ˜
> b:=vector([x0,xp0]);
b := [ x0 xp0 ]
> C:=linsolve(A,b);
1 xp0 + w ˜ ζ˜ x0 + w ˜ ζ˜2 − 1 x0
C :=
2 w˜ ζ˜2 − 1
1 w ˜ ζ˜ x0 − w ˜ ζ˜2 − 1 x0 + xp0
−
2 ζ˜2 − 1 w ˜
2
Aufgabe 9 Überprüfen Sie, ob Maple richtig ‘gerechnet’ hat!
2
2.2.3 Charakteristika
2.2.3.1 Stabilität
jω
ζ=0
00
11
00ω n
11
s1
0<ζ<1 1
0 1−ζ 2 ω n
ωn
(si ) < 0, i,
wobei (s) der Realteil der komplexen Größe s ist. D.h. die Nullstellen der cha-
rakteristischen Gleichung8 müssen negativen Realteil haben. Dieser Sachverhalt
ist in Abbildung 2.5 dargestellt für zwei Polpaare s1,2 , von denen eines auf dem
Halbkreis mit Radius ωn liegt und eines auf der negativen reellen Achse.
Das Verhältnis
x(t1 ) ζ !
δ = ln = ζωn Td = 2π (2.2.17)
x(t2 ) 1 − ζ2
wird als logarithmisches Dekrement bezeichnet. Für kleine ζ 1 gilt approxi-
mativ
δ 2πζ.
Wenn durch Gleichung (2.2.17) das logarithmische Dekrement gegeben ist, dann
kann daraus das Dämpfungsverhältnis berechnet werden:
δ
ζ= ,
4π2 + δ2
wobei für kleine δ wieder die approximative Linearisierung
δ
ζ
2π
gilt. Mit diesen Beziehungen kann aus dem Graphen x(t) der Lösung das Dämp- → Experiment
fungsverhältnis ζ experimentell durch Auslesen zweier konsekutiver Werte x(t1 )
und x(t2 ) gewonnen werden.
2.2.3.3 Energieverluste
Die zeitliche Änderung der Energie W ist das Produkt aus Kraft und Geschwin-
digkeit
d0
tW(t) = fv = −cẋ2 (t), (2.2.18)
dt
und die dissipierte Energie ∆W während einer Periode Td ist
Td
d0
∆W = tW(t)dt. (2.2.19)
0 dt
Für eine gedämpfte harmonische Bewegung sei x(t) = X sin ωd t wird das Inte-
gral (2.2.19) zu
Td
∆W = cX2 ω2d cos2 ωd tdt
0
2π
= cX2 ωd cos2 ωd td(ωd t)
0
= πcωd X2 .
Der Energieverlust ist proportional zum Quadrat der Bewegungsamplitude X;
weiterhin ist er abhängig von sämtlichen Systemparametern (ωn , ζ).
Die bisherige Modellbildung der Dämfungskraft geht aus von der Vorstellung
der Geschwindigkeitsproportionalität fc (t) = cẋ(t). Der Vorteil dieses Ansatzes
liegt in der Tatsache, daß die Bewegungsgleichung linear ist, der Nachteil ist
die geringere Realitätsnähe. Herrscht trockene Reibung vor mit einem Reibung-
skoeffizienten µ, dann ist die Dämpfungskraft auch richtungsabhängig:
⎧
⎨ µN ẋ < 0,
fc = 0 ẋ = 0,
⎩
−µN ẋ > 0,
dieser Dämpfungsansatz heißt Coulomb-Dämpfung, und mit ihm wird die Be-
wegungsgleichung (1.1.1) nichtlinear. Unter Einführung der Signum-Funktion
⎧
⎨ 1 f ür x > 0,
sign(x) = 0 f ür x = 0,
⎩
−1 f ür x < 0,
Diese Gleichung ist inhomogen, die Nichtlinearität ist einfach: die Gleichung ist
abschnittsweise linear:
mẍ + kx = 0 für ẋ = 0
und
mẍ + kx = µmg für ẋ < 0.
Damit wird auch die Lösung in Anlehnung an Gleichung (2.1.7) einfach
µmg
x(t) = A1 cos ωn t + A2 sin ωn t − für ẋ > 0, (2.2.21)
k
und
µmg
x(t) = A3 cos ωn t + A4 sin ωn t + für ẋ < 0, (2.2.22)
k
was durch Einsetzen in Gl. (2.2.20) nachprüfbar ist. Die Konstanten Ai , i =
1, . . . , 4 werden wieder aus Anfangsbedingungen gewonnen. Diese seien zum
Zeitpunkt t = 0
x(0) = x0 > 0,
und
ẋ(0) = 0,
d.h. die Bewegung startet von ‘rechts’ nach ‘links’ und damit muß Gleichung (2.2.22)
verwendet werden, aus der
µmg
A3 = x0 −
k
2.2 Gedämpfte Eigenbewegung 27
und
A4 = 0
erhalten werden, und x(t) in der ersten Halbperiode wird
µmg µmg
x(t) = x0 − cos ωn t + . (2.2.23)
k k
Bei t = T/2 = π/ωn findet Richtungsumkehr der Masse m statt; jetzt gilt
Gl. (2.2.21) mit Anfangsbedingungen aus Gl. (2.2.23) und t = π/ωn
π µmg
x1 = x( ) = −x0 + 2
ωn k
und
π
) = 0,
ẋ(
ωn
da bei Richtungsumkehr die Geschwindigkeit der Masse gerade Null ist.
Damit werden in der zweiten Halbperiode für π/ωn t 2π/ωn
µmg
A1 = x0 − 3
k
und
A2 = 0
und µmg µmg
x(t) = x0 − 3 cos ωn t − . (2.2.24)
k k
Am Ende der zweiten Halbperiode ist
µmg
x2 = x(2π/ωn ) = x0 − 4
k
und ẋ(2π/ωn ) wiederum Null. Dies sind die Anfangsbedingungen für den dritten
Halbzyklus, für den wieder Gl. (2.2.22) zu wählen ist.
Und so fort, bis die Bewegung zum Stillstand kommt, was dann geschieht, wenn
xn µmg/k ist, da dann die Rückstellkraft der Feder nicht mehr ausreicht, um
die Reibungskraft µmg zu überwinden, d.h.
2µmg µmg
x0 − n
k k
oder nach
µmg
x0 −
n= k
2µmg
k
Halbperioden.
In jeder Halbperiode reduziert sich die Amplitude X der Schwingung um 2µmg/k,
d.h.
4µmg
Xi+1 = Xi − ,
k
d.h. sie reduziert sich linear mit der Zeit, im Gegensatz zur viskosen Dämpfung,
bei der die Amplitude exponentiell abfällt, s.a. Gl. (2.2.15).
28 Ein-Freiheitsgrad-Systeme
1
x(t)
−1
−2
−3
−4
−5
0 5 10 15 20 25 30 35 40
t[s]
f(t) = F0 ej(ωt+φ)
oder
f(t) = F0 cos(ωt + φ) (2.3.1)
formuliert wird, mit F0 der Amplitude, ω der Frequenz und φ dem Phasenwin-
kel. Wenn zum Zeitpunkt t = 0 auch f(t) = 0 ist, dann ist auch φ = 0.
Wir werden nachweisen, daß harmonische Anregung immer eine harmonische
Bewegung x(t) bewirkt, die sich jedoch i.a. in Amplitude und Phase von der
Anregung unterscheidet. Wenn die anregende Frequenz ω mit der Eigenfrequenz
ωn des angeregten Systems zusammenfällt, dann überhöht sich die Amplitude
X sehr stark, und es kann zu schweren Schäden am System führen. Man spricht
von Resonanz.10
Die Bewegungsgleichung des fremd-erregten Systems lautet:
Deren allgemeine Lösung wird superponiert aus der Lösung xh (t) des homogenen
Systems
mẍh (t) + cẋh (t) + kxh (t) = 0
und der partikulären Lösung xp (t) des angeregten Systems (2.3.2):
Darin ist xh für den transienten Teil der Lösung und xp für deren stationären
Anteil verantwortlich.
Für das ungedämpfte System wird die Bewegungsgleichung (2.3.2) mit der An-
regung (2.3.1)
oder
ẍ(t)
+ x(t) = f0 cos ωt,
ω2n
wobei
F0
f0 = =∆ (2.3.4)
k
10 Resonanz von resonare (lat.) widerhallen.
30 Ein-Freiheitsgrad-Systeme
F0
A0 =
k − mω2
F0
x(t) = A1 cos ωn t + A2 sin ωn t + cos ωt. (2.3.5)
k − mω2
F0
A1 = x0 −
k − mω2
und
ẋ0
A2 =
ωn
und damit wird Gleichung (2.3.5) schließlich zu
F0
ẋ0 F0
x(t) = x0 − cos ωn t + sin ωn t + cos ωt. (2.3.6)
k − mω2 ωn k − mω2
Es ist leicht nachzuvollziehen, daß mit dem Verhältnis von Erreger- und Eigen-
kreisfrequenz
η = ω/ωn
V(η) =
A0
f0
= 1ω 2
=
1
1 − η2
(2.3.7)
1−
ωn
gilt; in Abb. 2.7 ist V über η abgebildet. Das Amplitudenverhältnis V wird auch
als Vergrößerungsfunktion bezeichnet.
Die Abbildung zeigt zwei Bereiche, die durch einen Pol bei ω = ωn getrennt
sind.
2.3 Harmonische Anregung 31
ω=ω
n
V(η)
−1
−2
−3
−4
−5
0 0.5 1 1.5 2 2.5 3
η
2.3.1.1 Resonanz
es gilt η < 1, und x(t) berechnet sich nach (2.3.8) und V ist im linken Ast der
Abb. 2.7 gezeichnet. Der Systemausgang x(t) ist in Phase mit dem Systemein-
gang f(t).
Jetzt ist η > 1, und die partikuläre Lösung wird xp (t) = −A0 cos ωt, hier mit ωn < ω < ∞
f0
A0 =
η2 − 1
Das Amplitudenverhältnis ist im rechten Ast abgebildet. Für sehr große ω ten-
diert A0 gegen Null. A ist negativ, d.h. für eine positive ‘statische’ Amplitude
f0 ist die ‘dynamische’ Amplitude A0 negativ; man sagt, Eingang und Ausgang
Æ
sind (180 ) phasenverschoben.
In diesem Fall kann der Ansatz xp (t) = A0 cos ωt nicht zum Ziel führen, da er ωn = ω
32 Ein-Freiheitsgrad-Systeme
auch die homogene Gleichung löst. Die modifizierte Ansatzfunktion lautet nun
[28]
xp (t) = A0 t sin ωt,
und damit die Lösungsfunktion x(t) = xh (t) + xp (t)
f0
A0 = .
2ω
Damit wird die Lösung für ω = ωn
f0
x(t) = A1 cos ωt + A2 sin ωt + t sin ωt, (2.3.9)
2ω
wobei die Amplituden A1 und A2 wieder aus den beiden Anfangsbedingungen
x0 und ẋ0 gewonnen werden, die in Gleichung (2.3.9) und deren Ableitung für
t = 0 eingesetzt werden. Damit erhält man schließlich
ẋ0 f0
x(t) = x0 cos ωt + sin ωt + t sin ωt (2.3.10)
ω 2ω
als Lösung. In ihr kommt ωn nur scheinbar nicht vor, da es gleich der anregenden
Kreisfrequenz ω ist. Man erkennt leicht, warum Resonanz im ungedämpften Fall
katastrophale Auswirkungen hat: Der dritte Summand in der Lösung (2.3.10)
wächst für t → ∞ über alle Schranken, s.a. Abb. 2.8
2.3.1.2 Schwebung
Wir kehren zurück zur Gleichung (2.3.6) und setzen dort die Anfangsbedingun-
gen zur Vereinfachung — aber ohne Beschränkung der Allgemeinheit — zu Null.
Dann erhalten wir nach kurzer Umrechnung
F0
x(t) = (cos ωt − cos ωn t),
m(ω2n − ω2 )
was unter Verwendung eines Theorems der Winkelfunktionen als Produkt dars-
ω ω
tellbar ist:
2F0 n +ω n −ω
x(t) = sin t sin t . (2.3.11)
m(ω2n − ω2 ) 2 2
Nun sei die erregende Frequenz ω nur wenig kleiner als die Eigenfrequenz ωn :
ωn − ω = 2
ωn + ω 2ω
2.3 Harmonische Anregung 33
50
40
30
20
10
x(ω t)
−10
−20
−30
−40
−50
0 5 10 15 20 25 30 35 40 45 50
ωt
ω2n − ω2 4ω.
F0
x(t) = sin t sin ωt.
2mω
F0
A(t) = sin t.
2mω
Die Periode der Amplitude ist TA 2π/. Diese Schwingungsform ist Ihnen aus
dem Physikunterricht als Schwebung bekannt. Eine typische Schwebung ist in
Abb. 2.9 dargestellt.
34 Ein-Freiheitsgrad-Systeme
0.8
0.6
0.4
0.2
x(t)
−0.2
−0.4
−0.6
−0.8
−1
0 5 10 15 20 25 30 35
t[s]
Die partikuläre Lösung dieser Gleichung wird wegen des Dämpfungsterms cẋ(t)
nicht mehr in Phase mit der Anregung schwingen, sondern wir müssen einen
Phasenwinkel φ0 in den Ansatz einführen:
xp (t) = A0 cos(ωt − φ0 ).
Die Konstanten A0 und φ0 werden wieder durch Einsetzen von xp und deren
Ableitungen in Gleichung (2.3.12) bestimmt:
Dies ist Gleichungssystem mit den beiden Unbekannten A0 und φ0 mit der
Lösung
F0
A0 = (2.3.13)
(k − mω2 )2 + (cω)2
und
cω
φ0 = arctan . (2.3.14)
k − mω2
Übung 1 Überprüfen Sie die Lösung mit Bleistift und Papier.
Die Gesamtlösung der Bewegungsgleichung ergibt sich wieder aus der Überla-
gerung von homogener und partikulärer Lösung x = xh + xp , mit xh z.B. aus
Abschnitt 2.2.2.3 Gleichung (2.2.12)
A0 1
V(ζ, η) = = (2.3.16)
f0 (1 − η2 )2 + (2ζη)2
und
2ζη
φ0 = arctan . (2.3.17)
1 − η2
Übung 2 Vollziehen Sie bitte die Ausdrücke (2.3.16) und (2.3.17) mit Bleistift
und Papier nach.
4.5
3.5
←ζ=0.1
3
V(η)
2.5
←ζ=0.2
2
←ζ=0.3
1.5
←ζ=0.5
1
←ζ=0.707
0.5
0
0 0.5 1 1.5 2 2.5 3
η
• Die Steigung des Graphen von V im Ursprung ist für alle ζ Null,
• Für ζ Ô1
2
nimmt V mit wachsenden η monoton ab,
• Für ζ < Ô1
2
erreicht V ein Maximum bei
η= 1 − 2ζ2 . (2.3.18)
Übung 3 Leiten Sie die Beziehung (2.3.18) ab. Wie groß ist der Wert dieses
Maximums?
2
2.3 Harmonische Anregung 37
180
160
120
100
φ0(η)
80
60
40
←ζ=0.2
20
←ζ=0.1
0
0 0.5 1 1.5 2 2.5 3
η
Æ
• Für große Werte von η nähert sich φ0 → 180 , d.h. Anregung und Antwort
sind außer Phase.
38 Ein-Freiheitsgrad-Systeme
Sie wird im stationären Fall aus der partikulären Lösung xp (t) = A0 cos(ωt −
φ0 ) und deren Ableitung ẋp (t) = −ωA0 sin(ωt − φ0 ) zu
berechnet. Die zweite Zeile dieser Gleichung zeigt, daß die beiden Vektorsum-
manden orthogonal zueinander sind (π/2!), und der Betrag der Summe ist
FB = A0 (k)2 + (cω)2
woraus das Verhältnis der Amplitude der anregenden Kraft F0 zur übertragenen
Kraft FB gebildet werden kann:
FB 1 + (2ζη)2
= , (2.3.19)
F0 (1 − η2 )2 + (2ζη)2
1
Vmax (2.3.20)
2ζ
approximiert werden kann. Dieser Wert wird auch als Qualitätsfaktor
Q be-
zeichnet. Die Frequenzwerte ω1 und ω2 , die den Werten Q2 = Q/ 2 auf dem
V-Graphen sind Frequenzen halber Leistung, da die durch Dämpfung dissipierte
Leistung ∆W = πcωd X2 proportional zum Quadrat der Amplitude X ist, s.a.
Abschnitt 2.2.3.3.
Der Sachverhalt ist in Abb. 2.12 skizziert.
2.3 Harmonische Anregung 39
Q2
∆ω
η1 1 η2 η
Abbildung 2.12: Bandbreite
die
η21 = 1 − 2ζ2 − 2ζ 1 + ζ2
und
η22 = 1 − 2ζ2 + 2ζ 1 + ζ2
sind und für ζ 1 als
und außerdem
ω2 + ω1 2ωn ,
ωn 1
Q (2.3.22)
ω2 − ω1 2ζ
d.h.
f0 = 0.02
2ζ = 0.04ζ.
Andererseits ist
A0
=
1
f0 (1 − η2 )2 + (2ζη)2 η=0.75
0.01 1
=
0.04ζ 0.1914 + 2.25ζ2
Die Euler-Identität
e¦jωt = cos ωt j sin ωt
2.3 Harmonische Anregung 41
erlaubt es, den Realteil F0 cos ωt und den Imaginärteil F0 sin ωt einer harmoni-
schen Anregung zusammenzufassen zur Exponentialfunktion f(t) = F0 ejωt und
damit die Bewegungsgleichung (2.3.12) in folgende Form zu bringen
Der Realteil der Lösung x(t) dieser Gleichung entspricht dabei der physikalischen
Lösung.
Wie in Abschnitt 2.3.2.1 wird zur Lösung der — jetzt komplexen — Bewegung-
sgleichung die partikuläre Lösung als
angesetzt. Dieser Ansatz und dessen erste und zweite Ableitungen werden in die
Bewegungsgleichung (2.3.24) eingesetzt. Unter der Bedingung, daß diese für alle
Zeitpunkte t gelten muß, d.h. daß ejωt nicht identisch verschwinden darf, wird
für die Ausgangsamplitude A0 folgender Ausdruck erhalten:
F0
A0 =
k − mω2 + jcω
k − mω2 cω
= F0 − j F0
(k − mω ) + (cω)
2 2 2 (k − mω2 )2 + (cω)2
F0
= e−jφ , (2.3.26)
(k − mω2 )2 + (cω)2
wobei die letzte Formelzeile aus der Äquivalenz der Darstellung komplexer Va-
riablen in kartesischen und in Polarkoordinaten x + jy = rejφ mit r = x2 + y2
und φ = arctan(y/x) folgt und für
cω
φ = arctan
k − mω2
resultiert. A0 in den Ansatz (2.3.25) eingesetzt ergibt für die partikuläre Lösung
F0
xp (t) = ej(ωt−φ) (2.3.27)
(k − mω2 )2 + (cω)2
2.3.3.2 Frequenzgang
H(jω) = |H(jω)|e−jφ
2ζη
φ = arctan
1 − η2
beträgt. Damit kann die partikuläre Lösung (2.3.27) auch folgendermaßen ges-
chrieben werden:
xp (t) = f0 |H(jω)|ej(ωt−φ) . (2.3.29)
Spezialfälle dieser Gleichung sind
f0
xp (t) = f0 |H(jω)|ej(ωt−φ) = cos(ωt − φ) (2.3.30)
(1 − η2 )2 + (2ζη)2
lautet
ms2 X(s) + csX(s) + kX(s) = F(s),
und heißt als Verhältnis von Ausgang X zu Eingang F geschrieben Übertragung-
sfunktion
X(s) 1
G(s) = = .
F(s) ms2 + cs + k
2.3 Harmonische Anregung 43
X(jω) 1
k = .
F(jω) 1 − η + 2jζη
2
kG(jω) = H(jω).
Die Analyse der Lösung (2.3.29) läßt eine anschauliche geometrische Deutung in
der komplexen Ebene zu. Xp wird zweimal differenziert: Es wird damit zunächst
die Geschwindigkeit ẋp und dann die Beschleunigung ẍp gewonnen:
und
ẍp (t) = j2 ω2 f0 |H(jω)|ej(ωt−φ) = −ω2 xp (t).
In vielen Fällen wirkt nicht eine unmittelbare äußere Kraft als Anregung auf das
dynamische System, sondern solche Kräfte ergeben sich mittelbar, z.B. durch
harmonische Erregung der Systembasis (Fußpunkterregung) oder eine rotato-
rische Unwucht im System.
44 Ein-Freiheitsgrad-Systeme
Im
x(t)
f(t)
kx(t)
φ
.
cx(t)
ωt ..
mx+kx
Re
..
mx(t)
Eine harmonische Erregung der Systembasis ist in Abb. 2.14 dargestellt, in der
ein stark vereinfachtes Fahrzeug über eine sinusförmige Fahrbahn fährt.
Die freigeschnittene Fahrzeugmasse mit der an ihr angreifenden Kräften ist dort
ebenfalls angegeben, und damit ergibt die Modellbildung folgende Bilanz:
woraus
folgt. Die Anregungsfunktion sei y(t) = Y sin ωt, und in Gl. (2.3.32) eingesetzt
ergibt sich
Dies ist ein dynamisches System mit zwei Anregungen, also müssen auch zwei
partikuläre Lösungen angesetzt werden, die allerdings — wegen der Linearität
2.3 Harmonische Anregung 45
..
mx
x(t)
m
k c
y(t)
. .
k(x-y) c(x-y)
11
00
00
11
00
11
00
11
1 + (2ζη)2
= Y cos(ωt − φ1 − φ2 ), (2.3.33)
(1 − η2 )2 + (2ζη)2
mit
cω 2ζη
φ1 = arctan = arctan
k − mω2 1 − η2
und
k 1
φ2 = arctan = arctan ,
cω 2ζη
wobei wieder das Frequenzverhältnis η = ω/ωn und das Dämpfungsverhältnis
ζ = c/ck sind.
Übung 6 Berechnen Sie die beiden zu cωY cos ωt und kY sin ωt gehörenden par-
tikulären Lösungen x1p (t) und x2p (t) und damit die Lösung xp (t) = x1p (t)+x2p (t).
2
Übung 7 Fassen Sie mit Hilfe einer guten Formelsammlung[6]11 die beiden Pha-
11 Der Maple-Befehl combine(expression,arctan,’symbolic’) führt auch zum Erfolg.
46 Ein-Freiheitsgrad-Systeme
1 + (4ζ2 − 1)η2
φ0 = − arctan . (2.3.34)
2ζη3
180
160
140
←ζ=0.1
120
←ζ=0.2
100
φ0(η)
←ζ=0.3
80
←ζ=0.5
60
←ζ=0.707
40
20
0
0 0.5 1 1.5 2 2.5 3
η
wird ins Verhältnis zu der der Anregung Y gesetzt und ergibt dann den Auslen-
kungsfaktor12 V2 , eine weitere Vergrößerungsfunktion
A0 1 + (2ζη)2
V2 (η, ζ) = = , (2.3.35)
Y (1 − η2 )2 + (2ζη)2
die beschreibt, wie die Auslenkung der Basis auf das System in Abhängigkeit von
der Anregungsfrequenz übertragen wird. Sie ist in Abb. 2.16 für verschiedene
Dämpfungsverhältnisse dargestellt.
4.5
3.5
←ζ=0.1
3
V2(η)
2.5
2
←ζ=0.2
1.5 ←ζ=0.5
←ζ=0.3
1
←ζ=0.707
0.5
0
0 0.5 1 1.5 2 2.5 3
η
größerungsfunktionen ist uneinheitlich. Hier wird der Vorschlag von Wittenburg[49] übernom-
men.
48 Ein-Freiheitsgrad-Systeme
1 + 8 ζ2
V2max
1− 1
2
1 + 8 ζ2 − 1
+ 1 + 8 ζ2 − 1
4 ζ2
mω2
z(t) = Y sin(ωt − φ1 ) = Z sin(ωt − φ1 ),
(k − mω2 )2 + (cω)2
mit
cω 2ζη
φ1 = arctan = arctan .
k − mω 2 1 − η2
Die Vergrößerungsfunktion V3
Z mω2 η2
V3 = = =
Y (k − mω2 )2 + (cω)2 (1 − η2 )2 + (2ζη)2
ist in Abb. 2.17 für unterschiedliche Dämpfungsverhältnisse aufgezeichnet. Auch
hier zeigt die gestrichelte Linie den Graphen des geometrischen Ortes des Maxi-
mums von V3 . Der Graph des Phasenwinkels φ1 ist identisch mit dem in Abb. 2.11
Die Relativbewegung ist durch folgende Merkmale gekennzeichnet:
• Für kleine Frequenzen ω 1 ist Z 0, unabhängig von ζ,
• In der Resonanzregion, ω ωn , steigt die Vergrößerungsfunktion sehr
stark mit abnehmendem Dämpfungsverhältnis ζ; um diesen Effekt zu verk-
leinern, müssen zusätzliche Dämpfungsmaßnahmen ergriffen werden,
• Für große Frequenzen ω → ∞ tendiert die Vergrößerungsfunktion V3 → 1,
unabhängig von ζ,
• Für 0 < ζ < 2/2 erreicht die Vergrößerungsfunktion V3 ein Maximum
bei
1
η= ,
1 − 2ζ2
2.3 Harmonische Anregung 49
4.5
←ζ=0.1
4
3.5
3
V3(η)
2.5 ←ζ=0.2
2
←ζ=0.3
1.5
1 ←ζ=0.5
←ζ=0.707
0.5
0
0 0.5 1 1.5 2 2.5 3
η
• Für ζ > 2/2 steigt die Vergrößerungsfunktion monoton an.
Übung 9 Zeigen Sie, daß der Wert des Maximums der Vergrößerungsfunktion
1
V3max =
2ζ 1 − ζ2
ist.
2.3.4.1.3 Kraftübertragung II Aus der Abb. 2.14 ist ersichtlich, daß die bes-
chleunigte Masse m eine Kraft f
auf die Basis ausübt. Sie kann — im stationären Zustand — aus der partikulären
Lösung (2.3.33) durch zweimalige Differentiation gewonnen werden
1 + (2ζη)2
f(t) = mω2 Y cos(ωt − φ1 − φ2 )
(1 − η2 )2 + (2ζη)2
1 + (2ζη)2
= kYη 2
cos(ωt − φ1 − φ2 )
(1 − η2 )2 + (2ζη)2
= FT cos(ωt − φ1 − φ2 ),
4.5
3.5
←ζ=0.1
←ζ=0.707
3
←ζ=0.2
←ζ=0.5
V4(η)
2.5
2
←ζ=0.3
1.5
0.5
0
0 0.5 1 1.5 2 2.5 3
η
emr η2
xp (t) = A0 cos(ωt − φ0 ) = cos(ωt − φ0 ), (2.3.38)
M (1 − η2 )2 + (2ζη)2
mit dem Phasenwinkel φ0
2ζη
φ0 = arctan . (2.3.39)
1 − η2
52 Ein-Freiheitsgrad-Systeme
ωt
1111
0 01
0000
1 111
000
10
e
000
111 000
111 ξ ξ
0
1000
111
000
111
0 mr/2
1 m
000
111
000
111
000
111 1010 x(t)
0
1 10 η η
0
1
0
1 1010 µ µ
0
1 10
k c
.
kx cx
11111111111111
00000000000000
Abbildung 2.19: Rotatorische Unwucht
Die Vergrößerungsfunktion
MA0 η2
= ,
mr e (1 − η2 )2 + (2ζη)2
mX 1
= (2.3.40)
me max 2ζ 1 − ζ2
bei
1
ηmax = > 1, (2.3.41)
1 − 2ζ2
4. Für ζ > 2/2 wächst die Vergrößerungsfunktion monoton.
2.3 Harmonische Anregung 53
2.3.5.1 Coulomb-Dämpfung
wobei die Reibungskraft µN addiert wird, wenn ẋ > 0 gilt und sie subtrahiert
wird, wenn ẋ < 0 ist.
2.3.5.2 Hysterese-Dämpfung
F = kx + cẋ
Die stationäre Lösung sei von der Form x(t) = X sin ωt, dann wird
Nach Quadrieren und Umordnen ergibt sich die Gleichung einer rotierten Ellipse:
deren Graph in Abb. 2.20 dargestellt ist, und deren Fläche den oben berechneten
Wert ∆E einnimmt. Wenn keine viskose Dämpfung wirkt (c = 0), dann kollabiert
die Ellipse auf die eingezeichnete Gerade mit der Steigung k, der Kennlinie der
Feder.
54 Ein-Freiheitsgrad-Systeme
F(t)
cωX
X x(t)
∆E = πkβX2 = πhX2 ,
3.5
3
←ζ=0.3
2.5
V5(η)
←ζ=0.5
1.5
←ζ=0.707
1
←ζ=1
0.5 ←ζ=2
0
0 0.5 1 1.5 2 2.5 3
η
wobei
c^
ζ^ =
2 mk
und
ω ω
η= =
ωn k/m
sind.
In den Abb. 2.21 und 2.22 sind Amplitudengang und Phasengang für verschie-
dene Werte von ζ dargestellt.
Sie erkennen folgendes:
• Die Vergrößerungsfunktion V5 = X/f0 erreicht ihr Maximum
1
Vmax =
b
bei ω = ωn unabhängig von ζ,
• Der Phasenwinkel bei ω = 0 beträgt
φ = arctan β,
180
160
←ζ=0.707
140
←ζ=1
←ζ=2
120
100
Φ(η)
80
60
40
←ζ=0.5
←ζ=0.3
20
0
0 0.5 1 1.5 2 2.5 3
η
2.3.5.3 Luft-Dämpfung
Wenn ein System in Luft oder einem anderen viskosen Medium schwingt, dann
ist die Dämpfungskraft proportional dem Quadrat der Geschwindigkeit ẋ und
die Systemgleichung wird folgendermaßen formuliert:
Die Dämpfungskraft
1
Fd = αsignẋẋ2 = cw ρAsignẋẋ2 ,
2
die der Geschwindigkeitsrichtung entgegenwirkt, hängt vom Widerstandsbeiwert
cw , von der Dichte ρ des Mediums und der Querschnittsfläche A des Systems
ab. Die Berechnung des Energieverlustes geschieht nach dem Beispiel von [25]:
8
∆E = αω2 X3 ,
3
2.3 Harmonische Anregung 57
8 αωX
c^ = . (2.3.44)
3 π
Aufgabe 15 Versuchen Sie aus der Literatur, z.B. [1],[5] eine möglichst vollstän-
dige und übersichtliche Zusammenstellung der unterschiedlichen Dämpfungsfor-
men zu kompilieren.
auftritt, wobei ρ die Dichte des strömenden Mediums ist, v dessen Geschwin-
digkeit, η dessen dynamische Zähigkeit und d der Bauteildurchmesser ist. Dabei
ist die Frequenz f der abgehenden Wirbel durch die Strouhalzahl
fd
Sr = = 0.21 (2.3.47)
v
58 Ein-Freiheitsgrad-Systeme
1
f(t) = cA ρv2 A sin ωt, (2.3.48)
2
mit cA = 1, dem Auftriebskoeffizienten von Zylindern und A der Querschnitts-
fläche.
Beispiel 7 Für das in Abb. 2.23 dargestellte mechanische System einer an einer
Stange der Länge l und Masse m2 aufgehängten Punktmasse m1 , die mit der
harmonischen Kraft f(t) = F0 sin ωt angeregt wird, soll die Bewegungsgleichung
erstellt und die Amplitude Θ der stationären Lösung berechnet werden.
11111
00000 1
0
0
1
0
1
1
0 0
1
0 F(t)
1
0
1 0
1
k2
0
1
0
1 0
1
0
1
0
1
0
1
0
1 11
00
00
11
m2 0
1 m1
0
1
0
1
θ 0
1
0
1
0
1
0
1 0
1
0
1 1
0
0
1 k1 0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1 0
1
0
1 0
1
0
1
a 1111
0000
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
1
0 0
1
0
1 0
1
0
1 b
0
1 0
1
0
1 0
1
0
1
0
1
l 0
1
0
1
0
1
Abbildung 2.23: Dynamisches System
θp = Θ sin ωt
F0 l
Θ= .
k1 a2 + k2 b2 − (J0 + m1 l2 )ω2
500
Θ= = −8.572
10−4 rad.
5000(0.252 + 0.52 ) − (3.333 + 50)104.722
Besorgen Sie sich14 das Computer-Algebra-System MuPAD[19] und lassen Sie sich MuPAD
den Lösunggraphen zeichnen. Das geht ganz einfach:
Mit read( beispiel.mu ) lesen Sie am MuPAD-Prompt folgendes Beispiel-File
ein:
a:=0.25:
b:=0.5:
l:=1.0:
M:=50.0:
m:=10.0:
k:=5000.0:
F:=500.0:
w0:=200.0:
w:=2*PI*w0/60:
J:=m*l^2/3+M*l^2:
K:=k*(a^2+b^2):
f:=ode({J*diff(Q(t),t,t)+K*Q(t)=F*l*sin(w*t),Q(0)=0,D(Q)(0)=0},Q(t));
x:=solve(f);
graph:=[Mode=Curve,[t,op(x,1)],t=[0,2],Grid=[50],Smoothness=[20]];
plot2d(Axes=Origin,Ticks=5,Labeling=TRUE,Scaling=UnConstrained,graph);
Zunächst werden die Daten eingelesen, dann mit der Prozedur ode die Differen-
tialgleichung formuliert, in solve gelöst, und in plot2d wird der in Abb. 2.24
folgende Graph gezeichnet. Dabei variiert ω0 von 200 bis 800.
0.15
0.1
←ω0=200
0.05
←ω =400
0
←ω0=600
Θ(t)
0
←ω0=800
−0.05
−0.1
−0.15
0 0.5 1 1.5 2 2.5 3 3.5
t
Anregung: 1) Variieren Sie einige der Parameter des Beispiels, z.B. m1 und
wiederholen Sie dann die MuPAD-Sitzung.
2) Beantworten Sie dann wieder die o.a. Fragen.
2
Beispiel 8 Ein Beispiel für rotatorische Unwucht ist die Reifen/Rad-Kombination
an Kraftfahrzeugen. Ein vereinfachtes Modell ist in Abb. 2.25 dargestellt.
Die Masse M vereinigt alle festen Teile (Bremsen, Streben, etc.), die Masse m
alle beweglichen Teile (Felge, Reifen). Die Feder-Dämpfer-Kombination (c, k2 )
repräsentiert die Radaufhängung, die Feder k1 stellt die Reifensteifigkeit dar.
Der Schwerpunkt der Masse m habe eine (Unwucht-)Exzentrizität e von der
Drehachse, die bei bestimmten Drehzahlen des Rades (d.h. Fahrgeschwindigkei-
ten) unangenehme Resonanzerscheinungen verursacht.
Gesucht ist die maximale Unwucht-Amplitude und die Geschwindigkeit, bei der
sie auftritt. Folgende Zahlenwerte seien gegeben: M = 20 kg, m = 15 kg, k1 =
1.8
105 N/m, k2 = 2
104 N/m, c = 265 Ns/m (der Dämpfer sei beschädigt!),
e = 2 mm. Der Bodenabstand der Radnabe beträgt h = 0.28 m.
Zunächst ist die Eigenfrequenz des Gesamtsystemsωn = k/m= 200, 000/35 =
75.59 rad/s und das Dämpfungsverhältnis ζ = c/2 km = 265/2 200, 000
35 =
0.05 und damit erhält man die Amplitude X mit Gl. (2.3.41) und Gl. (2.3.38)
2.3 Harmonische Anregung 61
111111111111111111
000000000000000000
000000000000000000
111111111111111111
000000000000000000
111111111111111111
k2
c
M m ω
k1
11111111111111111
00000000000000000
00000000000000000
11111111111111111
00000000000000000
11111111111111111
e
bei ηmax = 1/ 1 − 2ζ2 = 1.0025 zu
em η2
X= = 8.57 mm.
M+m (1 − η2 )2 + (2ζη)2
Die Amplitude tritt bei ωmax = rωn = 75.78 rad/s auf, was folgender Fahrges-
chwindigkeit entspricht:
v = hωmax = 0.28
75.78 = 21.2 m/s.
Berechnen Sie den Koeffizienten β der Hysterese-Dämpfung und den der äqui-
valenten viskosen Dämpfung, wenn das System bei 10 Hz betrieben wird.
Hysterese-Dämpfung:
∆E = πkβX2
2.5 = π5
104 β0.0084
⇒β = 0.294
62 Ein-Freiheitsgrad-Systeme
Viskose Dämpfung:
kβ
c^ =
ω
5
104 0.294
=
20π
⇒ c^ = 198kg/s
2
Beispiel 10 Berechnung der Resonanz-Amplitude für quadratische Dämpfung:
Aus Beziehung (2.3.16) wird für ω = ωn , d.h. η = 1, und mit Ausdruck (2.3.44)
f0 f0 f0
X= = =
2ζω 2 c^ω 8/(3π)αω2 X
2
Beispiel 11 Stabilitätsbedingungen:
1. Die Stabilität eines Systems wird nur durch seine Eigenbewegung charak-
terisiert.
2. Ein System ist dann und nur dann stabil, wenn die Nullstellen des cha-
rakteristischen Polynoms negativen Realteil haben.
3. Das charakteristische Polynom q(s) ist der Nenner der Übertragungsfunk-
tion des Systems.
4. Im Fall des gedämpften mechanischen Systems mit einem Freiheitsgrad ist
c 1
c 2 k
s1,2 =− −4 .
2m 2 m2 m
Diese haben negativen Realteil wenn alle Polynomkoeffizienten (m, c, k) positiv
(oder negativ — was im physikalischen Kontext sinnlos ist) sind.
Aufgabe:Versuchen Sie diese Aussage zu beweisen.
2
Beispiel 12 Ein Stahlkamin mit der Höhe h = 20 m, dem Außendurchmesser
da = 0.80 m und dem Innendurchmesser di = 0.75 m werde von Wind angebla-
sen.
2.3 Harmonische Anregung 63
damit
(2.07
1011 )(0.00457)
ω1 = 1.875 2
= 12.415rad/s = 1.976Hz.
(7.9
103 )(0.06087)(20)4
Die Frequenz der Wirbelbildung ist nach Gleichung (2.3.47) fd = 0.21v und
damit
fd 1.976
0.8
v= = = 7.53m/s.
0.21 0.21
2
Übung 10 Stellen Sie die Bewegungsgleichung für das in Abb. 2.26 gezeichnete
mechanische System auf und bestimmen Sie die stationäre Lösung. Die Masse
m ist homogen verteilt. Das System werde harmonisch angeregt mit M(t) =
M0 sin ωt.
1111
0000 11111
00000
k k
1111111
0000000
M(t)
0000000
1111111
1
0
0000000
1111111
1111111111
0000000000
0000000
1111111
0
1
θ m
1
0
0
1 0
1 0
1
11111
00000
l/4 3l/4
Übung 11 Die Konfiguration der Basisanregung sei wie in Abb. 2.27 dargestellt,
d.h. die anregende Bewegung wird über einen Dämpfer auf die Masse übertragen.
Die obere Abstützung sei inertial fest angenommen.
111111111111111111111
000000000000000000000
000000000000000000000
111111111111111111111
000000000000000000000
111111111111111111111
111111111111111111111
000000000000000000000
x(t)
y(t)=Ysin ω t
11111
00000
111
000
000
111
000
111
Welche Kraft wird im stationären Zustand auf die Masse übertragen? Skizzieren
Sie auch den Amplitudengang.
Übung 12 Ein Kraftfahrzeug der Masse m kg fahre über eine unebene Straße,
die sinusförmig mit Wellenlänge l = 6 m und Amplitude Y = 0.1m angenähert
sei. Die KFZ-Aufhängung bestehe aus Feder (k = 4
105 N/m) und Dämpfer
(c = 20
103 Ns/m).
000000000
111111111
000000000
111111111
Lager
1010
1010
Welle
1010
1010
l=2m
d
1010
Rotor 1010 s=5mm
1010
000000000000000000000
111111111111111111111 10 01
A
000000000000000000000
111111111111111111111 1010 A
000000000000000000000
111111111111111111111
000000
111111
000000000000000000000
111111111111111111111 0000000
1111111 1010
000000
111111
000000000000000000000
111111111111111111111
000000
111111 B 0000000
1111111
0000000
1111111
B
10
000000
111111
1111111
0000000
000000
111111 1111111
0000000 0000000
1111111
0000000
1111111
Stator
Sie dann zwei Fahrzeuge mit den Massen m = {1000, 1500} kg bei Geschwindig-
keiten von v = {20, 80, 100, 150} km/h.
2
Übung 13 Eine Francis-Turbine, deren vereinfachter Schnitt in Abb. 2.28 dar-
gestellt ist, mit Rotor-Masse m = 250 kg hat eine nominelle Spaltbreite zwischen
Rotor und Stator von s = 5 mm. Es wirkt eine Unwucht von emr = 5 mmkg
auf den Rotor. Der Drehzahlbereich variiert von 600 n 6000 U/min.
Die Turbinen-Welle hat eine Länge von l = 2 m. Berechnen Sie den Wellen-
Durchmesser d der Turbine, damit diese im gesamten Drehzahlbereich berüh-
rungsfrei läuft.
Hinweis: Nehmen Sie die Welle wie einen im Lager eingespannten Balken an,
dessen freies Ende die Rotormasse belastet. In diesem Fall ist die äquivalente
Federsteifigkeit k gegeben durch
EI E
πd
4
k=3 =3 3 . (2.3.52)
l3 l 64
2
Übung 14 Selbsterregte Schwingungen
Das dynamische Eigen-Verhalten einer Flugzeugtragfläche sei sehr einfach durch
mẍ + c0 ẋ + kx = 0 (2.3.53)
ż = Az + bu (2.3.54)
1111
0000 11111
00000 11111
00000
k k
c
1111111
0000000
M(t)
0000000
1111111
1
0
0000000
1111111
0
1
θ
0000000 10
1111111 0
1
m
0
1
0
1 0
1 0
1
01
11111
00000 10101111111
0000000
l/4
l/4 3l/41 0
2.4.1 Sprung-Funktion
σ (t) σ (t)
1 1
0 t 0 a t
0 für t < 0
σ(t) = (2.4.1)
1 für t 1
15 Siewird auch Heaviside-Funktion genannt, nach dem englischen Physiker Oliver Heaviside,
1850–1925.
68 Ein-Freiheitsgrad-Systeme
und allgemein:
0 f ür t < a
σ(t − a) =
1 f ür t a
=
1
− e−st
s
∞
=
e−as
s
.
a
2.4.1.3 Anwendung
Es soll nun der Einfluß einer sprungförmigen Anregung f(t) auf ein gedämpftes
mechanisches System mit Hilfe der Methode der Laplace-Transformation unter-
sucht werden:
wird Laplace-transformiert
woraus
mẋ0 + (ms + c)x0 F(s)
X(s) = +
ms2 + cs + k ms2 + cs + k
= Xh (s) + Xp (s)
1 s + 2ζωn
= ẋ0 + 2 x0 (2.4.3)
s2 + 2ζωn s + ω2n s + 2ζωn s + ω2n
F(s)
+
m(s2 + 2ζωn s + ω2n )
folgt, wobei F(s) die Laplace-Transformierte der Anregungsfunktion f(t) ist. Die
Gesamtlösung X(s) setzt sich auch im Bild- oder Frequenzbereich mit der Va-
riablen s wieder additiv aus der homogenen Lösung Xh (s) und der partikulären
Lösung Xp (s) zusammen.
2.4 Nicht-harmonische Anregungsfunktionen 69
Wenn die Anfangsbedingungen Null sind, wird der Quotient aus Ausgangsfunk-
tion Xp (s) und Eingangsfunktion F(s) als Übertragungsfunktion G(s) bezeich-
net:
Xp (s) 1
G(s) = .
F(s) ms2 + cs + k
Zunächst wird also X(s) im Frequenzbereich gebildet, und daraufhin wird die
gesuchte Lösung x(t) im Zeitbereich durch inverse Laplace-Transformation ge-
wonnen:16
x(t) = −1 {X(s)}.
Dieser Schritt ist nun für X(s) in Gleichung (2.4.3) durchzuführen: Die erste
Zeile wird in Tab. A.1 ausgelesen, die zweite Zeile wird mit Hilfe von Gl. (A.2.6)
gelöst, nachdem
F(s)
m(s2 + 2ζωn s + ω2n )
in die Faktoren F1 (s) = F(s) und
1
F2 (s) =
m(s2 + 2ζωn s + ω2n )
faktorisiert wird, deren inverse Laplace-Transformation
t
−1 1
[F1 (s)F2 (s)] = f(τ)e−ζωn (t−τ) sin( 1 − ζ2 ωn (t − τ))dτ
m 1 − ζ 2 ωn 0
ist. Damit wird
x0
x(t) = e−ζωn t sin( 1 − ζ2 ωn t + φ)
1 − ζ2
ẋ0
+ e−ζωn t sin 1 − ζ2 ωn t (2.4.4)
1 − ζ 2 ωn
t
1
+ f(τ)e−ζωn (t−τ) sin 1 − ζ2 ωn (t − τ)dτ
m 1 − ζ 2 ωn 0
mit φ = arccos ζ. Der dritte Summand dieser Beziehung ist ein Duhamel- oder
Faltungsintegral.
Jetzt soll das Integral in der Lösung (2.4.4)
noch für die Sprungfunktion 0110 F(t)
f(t) =
F0 0 t < t0 1010
0 t0 t 1010
ausgewertet werden. Das geschieht wieder
im Bildbereich, also Gl. (2.4.3); die Laplace-
Transformierte [f(t)] ist
F0 e−t0 s
F(s) =
F0
s
−
s
. 11111111111111
00000000000000
00000000000000
11111111111111
11111111111111
00000000000000
16
11111111111111
00000000000000
Zumeist geschieht dies durch Nachschlagen in geeigneten Tabellen, z.B. Tab. A.1, u.U.
nach vorheriger Vereinfachung durch Partialbruchzerlegung, s.a. Anhang 31.
m
11111111111111
00000000000000
00000000000000 1010
11111111111111
11111111111111
00000000000000
00000000000000 1010x(t)
11111111111111
k/2 k/2
c
111111111111111
000000000000000
70 Ein-Freiheitsgrad-Systeme
F0
e −ζωn t
xp (t) = 1− sin( 1 − ζ2 ωn t + φ)
mω2n
F0
e 1−ζ 2
−ζωn (t−t0 )
− 1− sin( 1 − ζ2 ω n (t − t0 ) + φ)
mω2n
F0 1 − ζ2
= −e−ζωn t sin( 1 − ζ2 ωn t + φ)
m 1 − ζ2 ω2n
+ e−ζωn (t−t0 ) sin( 1 − ζ2 ωn (t − t0 ) + φ) ,
mit ⎧
⎨ 0 t<1
f(t) = 1 1t<2
⎩
0 2t
Es sind ωn = 2 und ζ = Ô
3
2 2
> 1 und mit
−s −2s
^F = e − e
s
erhält man
X(s) = G(s) e−s − e−2s
mit
1
G(s) =
s(s + 1)(s + 2)
und nach Partialbruchzerlegung, siehe Anhang A.3
1 1 1
G(s) = − + .
2s s + 1 2(s + 2)
1 1
g(t) = −1 [G(s)] = − e−t + e−2t
2 2
2.4 Nicht-harmonische Anregungsfunktionen 71
1
−1 −s
= G(s) − −1 e−2s G(s)
−(t−1)1
= −e + e−2(t−1) σ(t − 1)
21 2
1
− + e−(t−2) − e−2(t−2) σ(t − 2)
2 2
bedeutet.
Die Antwort des Systems, wieder mit MuPADmit dem nachfolgenden File gewon-
nen, ist in Abb. 2.31 dargestellt. Dort ist auch der Graph der Impulsantwort
(Pulsbreite=0) des Systems mit eingetragen.
f:=proc(t) begin
if t<1 then 0;
elif t<2 then
1/2-exp(-(t-1))+1/2*exp(-2*(t-1));
else
-exp(-(t-1))+1/2*exp(-2*(t-1))+exp(-(t-2))-1/2*exp(-2*(t-2));
end_if;
end_proc;
plot2d(Axes=Origin,Scaling=UnConstrained,Ticks=5,
Labeling=TRUE,Labels=["",""],[Mode=Curve,[t,hold(f(t))],
t=[0,10],Grid=[100],Style=[Lines]]);
Aufgabe 16 Wiederholen Sie Beispiel 13 für ein unterdämpftes System mit ζ < 1
mit z.B. q(s) = s2 + 2s
+ 2 und für ein kritisch gedämpftes System mit ζ = 1
mit z.B. q(s) = s2 + 2 2s + 2.
2.4.2 Impuls-Funktion
Eine Impulsfunktion liegt vor, wenn ein hoher Kraft- oder Momenten-Impuls
über eine sehr kurze Zeit erfolgt, z.B ein Hammerschlag auf ein massives Werkstück
oder der Stoß zweier Billiard-Kugeln. In elektrotechnischen Systemen geschieht
die schnelle Aufladung eines Kondensators durch eine Stromquelle impulsförmig.
72 Ein-Freiheitsgrad-Systeme
0.25
← Pulsbreite=0
0.2
← Pulsbreite=1
0.15
x(t)
0.1
0.05
0
0 1 2 3 4 5 6 7 8 9 10
t
2. ∞
F(t − τ)dt = ^F. (2.4.6)
−∞
2.4 Nicht-harmonische Anregungsfunktionen 73
F(t)
F
2ε
τ−ε τ+ε t
Ist ^
F = 1, dann spricht man vom Dirac-Impuls.17
1 1 − e−ks
[fk (t)] = e−t0 s − e−(t0 +k)s = e−t0 s . (2.4.7)
ks ks
Der Grenzwert
lim fk (t) = δ(t − t0 )
k→ 0
∞ f ür t = t0 ,
δ(t − t0 ) = (2.4.8)
0 sonst,
17 Paul Dirac (1902–1984), Englischer Physiker; erhielt 1933 den Nobel-Preis für seine Ar-
beiten zur Quanten-Mechanik.
74 Ein-Freiheitsgrad-Systeme
und ∞
δ(t − t0 )dt = 1. (2.4.9)
−∞
2.4.2.3 Anwendung
und abschnittsweise:
0 t < t0 ,
x(t) =
e−(t−t0 ) − e−2(t−t0 ) t0 < t.
1
0
0
1 J0
0
1
0
1
0
1 k2
11
00
00
11
r
00
11
0
1 00
11
0
1
00
11
O1
0
2r
111
000
1
0
0
1
0
1
0
1F(t)
0
1
0
1
0
1
0
1
0
1
0
1
m
1
0
0
1
0
1
0
1
0
1
0
1
0
1
x(t)
k1
1111111
0000000
0000000
1111111
Abbildung 2.33: Flaschenzug
∂
> deq:=m*diff(x(t),t,t)+k*x(t)=F*exp(-t);
2
deq := m˜ x( t ) + k ˜ x( t ) = F ˜ e( −t )
∂t2
> ini:=D(x)(0)=0,x(0)=0;
ini := D( x )( 0 ) = 0, x( 0 ) = 0
0.06
ungedämpft
0.04
← gedämpft
0.02
x(t)
−0.02
−0.04
−0.06
0 2 4 6 8 10 12 14 16 18 20
t
dsolve(deq,ini,x(t));
>
m˜ k ˜ t m˜ k ˜ t
F ˜ cos F ˜ m˜ sin
F ˜ e( −t ) m˜ m˜
x( t ) = − +
m˜ + k ˜ m˜ + k ˜ m˜ k ˜ ( m˜ + k ˜ )
Für die o.a. numerischen Werte ergibt sich dann der Graph in Abb. 2.34. In-
terpretieren Sie den Graphen!
2
Aufgabe 17 Lösen Sie das voranstehende Beispiel mit Bleistift und Papier. Wenn
Sie die Laplace-Transformation verwenden, dann werden Sie Ihr Ergebnis mit
dem von MuPAD vergleichen wollen. Verwenden Sie dazu die Funktionen trans-
form::laplace bzw. transform::ilaplace.
2
Aufgabe 18 Fügen Sie nun in Abb. 2.33 parallel zur Feder k2 einen Dämpfer
mit dem Dämpfungsbeiwert c = 100 Ns/m ein. Leiten sie die Bewegungsglei-
chung her und lösen Sie sie. Je nach Dämpfungsbeiwert könnte der Graph wie
in Abb. 2.34 eingetragen aussehen.
2
2.4 Nicht-harmonische Anregungsfunktionen 77
2.4.3.1 Impulsantwort
Die Impulsantwort h(t) ist die Antwort (Lösung der Systemgleichung) x(t) eines
linearen dynamischen Systems, z.B. mẍ + cẋ + kx = f(t), auf den Dirac-Impuls
f(t) = δ.
Im Frequenzbereich wird das System durch die Übertragungsfunktion G(s) =
Ê
X(s)/F(s) beschrieben. Da für f(t) = δ die Beziehung (2.4.11) F(s) = ∆ = 1 gilt,
stimmt formal die Impulsantwort H(s) = X(s) = G(s)F(s) im Frequenzbereich
mit ihrer Übertragungsfunktion überein: H(s) = G(s). Merke aber: Während
H(s) eine Lösung(sfunktion) ist, ist G(s) der Quotient zweier Systemgrößen.
Aus physikalischen Gründen darf für den Grenzübergang ∆ → 0 der zweite und
dritte Summand der linken Seite keinen Beitrag leisten, ansonsten müßte gelten
x(∆) → ∞, also:
m[ẋ(∆) − ẋ(0)] = 1,
und mit ẋ(0) = 0:
1
ẋ(∆) = .
m
Aus dieser Überlegung folgt, daß eine impulsförmige Systemanregung äquivalent
ist zur Eigenbewegung des Systems mit x(0) = 0 und ẋ(0) = 1/m.
Jetzt sei noch einmal die gedämpfte Eigenbewegung aus Abschnitt 2.2.2.3, Gl. (2.2.12), Eigenbewegung
betrachtet:
Für x(0) = 0 und ẋ(0) = 1/m erhält man nach kurzer Rechnung:
e−ζωn t
h(t) = x(t) = sin(ωd t),
mωd
welches wegen der gewählten Anfangsbedingungen auch die Impulsantwort h(t)
ist.
Nun sei die allgemeine Anregung f(t) betrachtet. Wir integrieren von t = 0 bis Anregung f(t)
78 Ein-Freiheitsgrad-Systeme
x(t) = f(0)∆h(t).
und weiter:
n−1
Ê
x(t) = f(i∆)h(t − i∆)∆.
i=0
Das Integral in Gl. (2.4.12) ist ein Faltungsintegral.20 Bisher waren beide An-
fangsbedingungen zu Null gesetzt, d.h. xp(t) = x(t).
Anwendung Mit der Impulsantwort für, z.B. das unterdämpfte System mẍ + cẋ + kx = δ,
1
h(t − τ) = e−ζωn (t−τ) sin ωd (t − τ) (2.4.13)
mωd
Beispiele mit Lösungen dieses Integrals für diverse f(t) finden Sie in Absch-
nitt 3.2.3.
Beispiel 15 Es soll das gedämpfte System mẍ + cẋ + kx = σ(t), mit σ(t) dem
Einheitssprung (2.4.1) als Anregung und mit Anfangsbedingungen zu Null ge-
setzt, berechnet werden.
Das Integral in Gl. (2.4.14) ist für t 0 auszuwerten, d.h. in diesem Fall ist
f = 1, und Maple löst uns das Integral wie folgt:21
>
> int(exp(zeta*omega*tau)*
> sin(sqrt(1-zeta^2)*omega*(t-tau)),tau=0..t);
t
1 − ζ2 e(ζωn t)
ζωn τ
e sin ωd (t − τ)dτ =
0 ωn
1 − ζ2 cos( 1 − ζ2 ωn t) + ζ sin( 1 − ζ2 ωn t)
−
ωn
20 Englisch:convolution integral.
21 Die weitere Umwandlung in die Form sin(ωd t + φ) muß man noch per Hand erledigen;
Anfrage bei Waterloo Software läuft.
2.4 Nicht-harmonische Anregungsfunktionen 79
1
1
x(t) = 1− e−ζωn t sin(ωd t + φ)
k 1 − ζ2
mit
1 − ζ2
φ = arctan .
ζ
2
Die Gesamtlösung x(t) = xh + xp des Systems Gesamtlösung
ist mit Gl. (2.2.12) für den unterdämpften Fall, 0 < ζ < 1,
ẋ 2
+ ζωn x0
e−ζωn t sin (ωd t + φ)
0
x(t) = x20 +
ωd
t (2.4.15)
1
+ e−ζωn t f(τ)eζωn τ sin ωd (t − τ)dτ,
mωd 0
mit
x0 ωd
φ = arctan .
ẋ0 + ζωn x0
In Abb. 2.35 ist eine allgemeine Anregungsfunktion abgebildet und dazu ihre
Approximation durch eine Treppenfunktion. Mit dieser Idee sind lineare dyna-
mische Systeme, die durch allgemeine Erregerfunktionen angetrieben werden,
abschnittsweise durch Laplace-Transformation und Anregung durch Sprung-
funktionen (jeweils die Differenz zur vorangegangenen Treppe, die natürlich auch
negativ sein kann) und nachfolgende inverse Laplace-Transformation lösbar.
Die Approximation der Anregungsfunktion durch eine (abschnittsweise kons-
tante) Treppenfunktion ist eine (grobe) Approximattion nullter Ordnung und
liefert nur dann brauchbare Ergebnisse, wenn die Zeitachse genügend fein dis-
kretisiert wird, was einen hohen rechnerischen Aufwand zur Folge hat.
Approximationen höherer Ordnung (Polygonzug, Spline-Funktionen, etc.) lassen
eine gröbere Diskretisierung der Zeitachse zu [12],[39].
F(t)
F(a) 00
11
11
00
00
11
00
11
00
11
00
11
00
11
00
11
00
11
00
11
00
11
00
11
11
00 t
0 a a+h
fangswertprobleme, für die eine Vielzahl von Verfahren, Methoden und Softwa-
repaketen zur Verfügung steht [23].22 Eine kurze Einführung ist in Anhang D
zusammengefaßt.
22 Im Internet ist die netlib mit der URL ‘http://www.netlib.org/liblist.html’ eine uner-
schöpfliche Quelle freier Software.
Kapitel 3
Zwei-Freiheitsgrad-Systeme
3.1 Eigenbewegung
3.1.1.1 Modellbildung
1
0
0
1
x1 x2
0
1
0
1
k1 k2
0
1
m1 m2
0
1
0
1
0
1
k1 x1 k 2 (x 2 -x 1 )
m1 m2
und umgeformt:
ist ein Spaltenvektor, der (um Raum zu sparen) im Text zum Zeilenvektor x =
¼
x1 , x2 transponiert wird. Die Anfangswertvektoren sind
x (0)
1
ẋ (0)
1
x0 = und ẋ0 =
x2 (0) ẋ2 (0)
Übung 16 Bestimmen Sie Massen- und Steifigkeitsmatrix für die in Abb. 3.2
dargestellte Beispielkonfiguration.
2
1
0
0
1
x1 x2
11
00
00
11
0
1 k1 k2 k3
00
11
0
1 00
11
0
1 00
11
m1 m2
0
1 00
11
0
1 00
11
0
1 00
11
Abbildung 3.2: Zwei-Freiheitsgrad-System
oder, da ejωt = 0, t,
(−ω2 M + K)u = 0. (3.1.5)
Da wir nicht an der trivialen Lösung u = 0 interessiert sind, darf die Inverse der
Matrix −ω2 M + K nicht existieren. Dies ist der Fall, wenn ihre Determinante
Null wird:
det(−ω2 M + K) = 0. (3.1.6)
Um dies einzusehen, wird Gl. (3.1.5) mit (−ω2 M + K)−1 von links vormultipli-
ziert:
(−ω2 M + K)−1 (−ω2 M + K)uejωt = (−ω2 M + K)−1 0
oder
Iu = 0,
was wieder nur die triviale Lösung u = 0 bedeutet. Die Inverse einer Matrix
A existiert nicht, wenn ihre Determinante |A| Null wird, was unmittelbar aus
Gl. (B.3.2) folgt. Dies sei gezeigt an folgendem:
Beispiel 16 Gegeben sei die Matrix A
a b
A=
c d
−1 1
d
−b
A =
|A| −c a
3.1.1.2.2 Eigenfrequenzen Aus der Gl. (3.1.6) läßt sich ω gewinnen. Dies sei
gezeigt am Beispiel des Abschn. 3.1.1.1:
Beispiel 17 Die Beziehung det(−ω2 M + K) = 0 ergibt mit den aktuellen Matri-
zen
−ω2 m1 + k1 + k2
−k2
−k2
−ω2 m2 + k2
=0
84 Zwei-Freiheitsgrad-Systeme
oder ausgewertet:
ω4 − 6ω2 + 8 = 0
oder
(ω2 − 2)(ω2 − 4) = 0
und daraus ω1 = 2 und ω2 = 2.
2
Die Gl. (3.1.7) wird als charakteristische Gleichung bezeichnet, die Frequenz
ω ist die natürliche Frequenz des Systems (3.1.2). Im nächsten Abschnitt wird
gezeigt daß zwischen dieser und der Eigenfrequenz λ einer transformierten Stei-
figkeitsmatrix folgender Zusammenhang besteht:
λ = ω2 . (3.1.8)
(−ω2i M + K)ui = 0.
In Übung 17 wurde gezeigt, daß die Matrix (−ω2i M + K) einen Rang-Abfall von
Zwei auf Eins hat, das heißt, daß die beiden Zeilen der Matrix linear abhängig
sind. Das bedeutet, daß von den beiden Komponenten des Vektors ui nur das
Verhältnis bestimmt werden kann, nicht ihre absolute Größe. Wir zeigen das
mit den Werten des Beispiels 17 am
Beispiel 18 Die Lösung der (homogenen) Matrixgleichung
(−ω2i M + K)ui = 0
9u11 − 3u12 = 0,
−3u11 + u12 = 0.
d.h. wir kennen nur die Richtung des Vektors u1 , nicht seine Größe. Für ω22 = 4
ergibt sich
−9
4 + 24 + 3 −3
u21 0
= ,
−3 −1
4 + 3 u22 0
oder
u21 1
=− .
u22 3
So können die beiden Vektoren u1 und u1 folgendermaßen geschrieben werden:
1
1
u1 = 3 = = ...
1 3
und − −1 1
1
u2 = 3 = = = ...
1 3 −3
1 1
x1 (t) = A1 sin( 2t + φ1 ) − A2 sin(2t + φ2 ), (3.1.11a)
3
3
x2 (t) = A1 sin( 2t + φ1 ) + A2 sin(2t + φ2 ). (3.1.11b)
Gl. (3.1.12a) und Gl. (3.1.14a) bilden ein Gleichungssystem mit den vier unbe-
kannten Integrationskonstanten A1 , A2 , φ1 und φ2 . Gelöst wird es mit einem
Paket für symbolische Algebra, z.B. Maple. Hier folgt der Ausdruck der Maple-
Sitzung:
> eqn1:=1=1/3*A1*sin(phi.1)-1/3*A2*sin(phi.2);
1 1
eqn1 := 1 = A1 sin(φ1) − A2 sin(φ2)
3 3
> eqn2:=0=A1*sin(phi.1)+A2*sin(phi.2);
> eqn3:=0=sqrt(2)/3*A1*cos(phi.1)-2/3*A2*cos(phi.2);
1 2
eqn3 := 0 = 2 A1 cos(φ1) − A2 cos(φ2)
3 3
> eqn4:=0=sqrt(2)*A1*cos(phi.1)+2*A2*cos(phi.2);
eqn4 := 0 = 2 A1 cos(φ1) + 2 A2 cos(φ2)
> solve({eqn1,eqn2,eqn3,eqn4},{A1,A2,phi.1,phi.2});
1 1 −3 −3
{φ2 = π, φ1 = − π, A2 = , A1 = },
2 2 2 2
1 1 3 −3
{φ1 = − π, φ2 = − π, A2 = , A1 = },
2 2 2 2
1 −3 1 3
{φ2 = π, A2 = , φ1 = π, A1 = },
2 2 2 2
1 3 1 3
{φ2 = − π, A2 = , φ1 = π, A1 = }
2 2 2 2
Von den vier Lösungsmengen wählen wir die dritte mit jeweils positiven Pha-
senwinkeln und erhalten damit als analytische Lösung der Differentialgleichung
3.1 Eigenbewegung 87
1 π 1 π 1 1
x1 (t) = sin( 2t + )+ sin(2t + )= cos 2t + cos 2t (3.1.15a)
2 2 2 2 2 2
3 π 3 π 3 3
x2 (t) = sin( 2t + )− sin(2t + )= cos 2t − cos 2t (3.1.15b)
2 2 2 2 2 2
2
← x2
← x1
x(t)
−1
−2
−3
0 2 4 6 8 10 12 14 16 18 20
t
3.1.1.3 Modal-Analyse
3.1.1.3.1 Eigenwerte und Eigenvektoren Aus der Linearen Algebra [21, 43]
wissen Sie, daß die Matrix −ω2 M + K oder — bekannter mit λ = ω2 als —
−λM + K dem allgemeinen Eigenwertproblem zugrunde liegt. Bevor wir uns
dessen Lösung zuwenden, sei zunächst das spezielle Eigenwertproblem
behandelt. Die Genese von Gl. (3.1.16) ist leicht nachzuvollziehen. Wie im Falle
der Systeme mit einem Freiheitsgrad, s.a. Gl. (2.2.9), wird die Systemgleichung
88 Zwei-Freiheitsgrad-Systeme
(−ω2 I + A)vejωt = 0
oder
Av = ω2 v = λv (3.1.21)
Wenn die Matrizen M und K, wie in unserem Fall, symmetrisch sind, dann ist
auch die Matrix A symmetrisch, und man spricht vom symmetrischen Eigen-
wertproblem. Es zeigt folgende Eigenschaften:
vi ¼ vj = 0, i = j, (3.1.22)
V = (v1 , v2 , . . . , vn ). (3.1.24)
Nun werden aus der Bedingung der Singularität der Matrix in der voranstehen-
den Gleichung ihre Eigenwerte bestimmt:
3 − λ
−1 −1
3−λ
= 0 =⇒ λ2 − 6λ + 8 = 0 =⇒ λ1,2 = {2, 4}.
1
1
v1 = .
2 1
3 − λ
Und nun
2 −1
v2 = 0 =⇒ −v21 − v22 = 0 =⇒ v21 = −v22 .
−1 3 − λ2
Und aus
v2
= 2v21 = 1 folgt v21 = 1/ 2 und damit wird
1
1
v2 = .
2 −1
1
Analog zur Transformation (3.1.17) x = M− 2 q gilt auch die Transformation
1 1
ui = M− 2 vi oder vi = M 2 ui . (3.1.25)
orthonormal ist.
Übung 18 Zeigen Sie die eben hergeleitete Aussage mit den Ergebnissen des
Beispiels 20.
2
V ¼ V = I. (3.1.27)
Übung 19 Zeigen Sie die Eigenschaft (3.1.27) anhand der Eigenvektoren aus
Beispiel 20.
2
Übung 20 Zeigen Sie, daß aus der Eigenschaft (3.1.27) folgende einfache Aus-
sage über die Inversion einer orthonormalen Matrix getroffen werden kann:
V −1 = V ¼ . (3.1.28)
2
Eine Transformation einer Matrix X mit einer nicht-singulären Matrix T
Y = T −1 XT (3.1.29)
wird als Ähnlichkeitstransformation bezeichnet. Die Matrix Y wird als zur Ma-
trix X ähnlich bezeichnet.
Folgende Eigenschaft bildet die Grundlage der Modalanalyse: Eine symmetrische
Matrix A ist durch die Ähnlichkeitstransformation mit der Modal-Matrix V
Λ = V −1 AV = V ¼ AV = diag(λ1 , λ2 , . . . , λn ) (3.1.30)
V p̈(t) + AVp(t) = 0
und diese Gleichung mit der Inversen, d.h. der Transponierten, von V von links
(vor-)multipliziert:
Damit ist die Entkopplung der Systemgleichungen vollendet, und jede dieser
kann für sich nach den Methoden des Kap. 2 gelöst werden. Die Behandlung der
Anfangsbedingungen wird in der Zusammenfassung des Algorithmus im nächs-
ten Paragraphen erläutert.
Beispiel 21 In Komponenten lautet Gl. (3.1.32c) für den Zwei-Freiheitsgrad-Fall
p̈ λ
1 1 0
p 0
1
+ =
p̈2 0 λ2 p2 0
oder
p̈1 + λ1 p1 = 0,
p̈2 + λ2 p2 = 0.
Mẍ(t) + Kx(t) = 0,
x(0) = x0 ,
ẋ(0) = ẋ0 .
1 1
1. Berechnung von M 2 und M− 2 ,2
1 1
2. Berechnung von A = M− 2 KM− 2 ,
3. Berechnung der Eigenwerte λi = ω2i und Eigenvektoren vi
der Matrix A,
4. Ordnen der Eigenvektoren in die Modal-Matrix
V = (v1 , v2 , . . . , vn ),
5. Berechnung von P = M− 2 V und P−1 = V ¼ M 2 ,
1 1
> S:=‘linalg/matpower‘(M,1/2);
S := result
> print(S);
3.1 Eigenbewegung 93
M =
1 0
0 1
K =
3 -2
-2 3
M1 =
1 0
0 1
M2 =
1
( 2a + 2d + 2 a2 − 2ad + d2 + 4b2 a2 − 2ad + d2 + 4b2
4
− a 2a + 2d − 2 a2 − 2ad + d2 + 4b2
+ a 2a + 2d + 2 a2 − 2ad + d2 + 4b2
+ d 2a + 2d − 2 a2 − 2ad + d2 + 4b2
− d 2a + 2d + 2 a2 − 2ad + d2 + 4b2
Æ
+ 2a + 2d − 2
a2 − 2ad + d2 + 4b2
1
a2 − 2ad + d2 + 4b2 )
a2 − 2ad + d2 + 4b2 , b(− 2a + 2d − 2 a2 − 2ad + d2 + 4b2
2
Æ
+ 2a + 2d + 2 a2 − 2ad + d2 + 4b2 ) a2 − 2ad + d2 + 4b2
1
b(− 2a + 2d − 2 a − 2ad + d + 4b
2 2 2
Æ a − 2ad + d + 4b , −1(
2
+ 2a + 2d + 2 a − 2ad + d + 4b )
2 2 2 2 2 2
4
− 2a + 2d + 2 a − 2ad + d + 4b a − 2ad + d + 4b
2 2 2 2 2 2
+ d 2a + 2d − 2 a − 2ad + d + 4b
2 2 2
− d 2a + 2d + 2 a − 2ad + d + 4b
2 2 2
− a 2a + 2d − 2 a − 2ad + d + 4b
2 2 2
+ a 2a + 2d + 2 a − 2ad + d + 4b
2 2 2
Æ
− 2a + 2d − 2 a − 2ad + d + 4b a − 2ad + d + 4b )
2 2 2 2 2 2
a2 − 2ad + d2 + 4b2
94 Zwei-Freiheitsgrad-Systeme
1 0
0 1
A =
3 -2
-2 3
V =
-0.7071 -0.7071
0.7071 -0.7071
D =
5 0
0 1
P =
-0.7071 -0.7071
0.7071 -0.7071
PI =
-0.7071 0.7071
-0.7071 -0.7071
x0 =
0
1
xp0 =
0
0
p0 =
0.7071
-0.7071
pp0 =
0
0
w1 =
2.2361
w2 =
1
und hier der Eingabe-File für die Sitzung:
M=[1,0;0,1]
K=[3,-2;-2,3]
M1=sqrt(M)
M2=inv(M1)
A=M2*K*M2
[V,D]=eig(A)
P=M2*V
PI=V’*M1
x0=[0;1]
xp0=[0;0]
r0=PI*x0
rp0=PI*xp0
t=linspace(0,20,2000);
w1=sqrt(D(1,1))
w2=sqrt(D(2,2))
r1=r0(1)*cos(w1.*t)+rp0(1)/w1*sin(w1.*t);
r2=r0(2)*cos(w2.*t)+rp0(2)/w2*sin(w2.*t);
3.1 Eigenbewegung 95
r1=sqrt(r0(1)^2+rp0(1)^2/w1^2)*sin(w1.*t+atan2(r0(1)*w1,rp0(1)));
r2=sqrt(r0(2)^2+rp0(2)^2/w2^2)*sin(w2.*t+atan2(r0(2)*w2,rp0(2)));
x1=P(1,1)*r1+P(1,2)*r2;
x2=P(2,1)*r1+P(2,2)*r2;
plot(t,x1,’-’,t,x2,’--’)
grid
xlabel(’t’)
ylabel(’x(t)’)
gtext(’\leftarrow x_1’)
gtext(’\leftarrow x_2’)
0.8
← x1
0.6
0.4
0.2
← x2
x(t)
−0.2
−0.4
−0.6
−0.8
−1
0 2 4 6 8 10 12 14 16 18 20
t
3.1.1.4 Beispiele
96 Zwei-Freiheitsgrad-Systeme
m J0
k1 k2
11111111111111111111111111111
00000000000000000000000000000
l1 l2
J0 θ̈ = f1 l1 − f2 l2 = k1 l1 (x − l1 θ) − k2 l2 (x + l2 θ).
Die Lösung ist wieder in der folgenden Matlab-Sitzung gegeben, wobei die Werte
aus Matlab heraus mit dem Befehl diary gespeichert werden:
1 m =
3.1 Eigenbewegung 97
l1 l2
θ (t)
f1
f2
2 1000
3 J =
4 300
5 k1 =
6 3000000
7 k2 =
8 2000000
9 l1 =
10 0.50000000000000
11 l2 =
12 0.80000000000000
13 M =
14 1000 0
15 0 300
16 K =
17 1.0e+06 *
18 5.00000000000000 0.10000000000000
19 0.10000000000000 2.03000000000000
20 M1 =
21 31.62277660168379 0
22 0 17.32050807568877
23 M2 =
24 0.03162277660168 0
25 0 0.05773502691896
26 A =
27 1.0e+03 *
28 5.00000000000000 0.18257418583506
29 0.18257418583506 6.76666666666666
30 V =
31 0.99481178318674 0.10173257115012
32 -0.10173257115012 0.99481178318674
98 Zwei-Freiheitsgrad-Systeme
33 D =
34 1.0e+03 *
35 4.98132939148435 0
36 0 6.78533727518232
37 U =
38 0.98301325847268 0.05592416811137
39 -0.18353455714641 0.99843501912796
40 E =
41 1.0e+03 *
42 4.98132939148435 0
43 0 6.78533727518232
44 Modeshape_1 =
45 -5.35601182554696
46 Modeshape_2 =
47 0.05601182554696
48 U =
49 0.03145871078044 0.00321706637060
50 -0.00587353273389 0.05743548508159
51 U_norm =
52 0.54488085412332 0.05572122405192
53 -0.10173257115012 0.99481178318674
%Werkzeugmaschine
format long
m=1000
J=300
k1=3000000
k2=2000000
l1=1/2
l2=8/10
M=[m,0;0,J]
K=[k1+k2,-k1*l1+k2*l2;-k1*l1+k2*l2,k1*l1^2+k2*l2^2]
M1=sqrt(M)
M2=inv(M1)
A=M2*K*M2
[V,D]=eig(A)
[U,E]=eig(K,M)
Modeshape_1=U(1,1)/U(2,1)
Modeshape_2=U(1,2)/U(2,2)
U=M2*V
U_norm=U/norm(U)
Der Unterschied der Zeilen 38,39 und 52,53 ist kein Rechenfehler, sondern beruht
auf der Definition der Matrix-Norm norm(U).
Aufgabe 19 Wiederholen Sie dieses Beispiel, wobei Sie nun die Koordinaten
nicht in den Schwerpunkt der Maschine legen, sondern um den Abstand e hori-
zontal nach links versetzt.
3.1 Eigenbewegung 99
3.1.1.4.2 Doppelpendel
11111111
00000000
l1
θ
1
S1
m1
Bestimmen Sie für das
W1 l2
θ nebenstehenden Doppel-
2 S2
pendel die Bewegungsglei-
x1 chung in den Koordi-
m2
naten x1 und x2 für
x2 W2 kleine Ausschläge. Berech-
nen Sie die natürlichen
Frequenzen, die Amplitu-
denverhältnisse und die
Lage der Schwingungskno-
ten für m1 = m2 = m
und l1 = l2 = l. Die Seil-
kräfte sind mit S1 und S2
bezeichnet.
Das Momentengleichgewicht um den Aufhängepunkt lautet:
Mit den Beziehungen θ1 = x1 /l1 und θ2 = (x2 − x1 )/l2 und der Approximation
S2 W2 erhalten wir
m l l1 + l2
+ xx = 00 . (3.1.33)
l1
1 1 0 ẍ W 1 1 + W2 −W2 1
l2 l2
0 m l ẍ
2 2 2 −W2 +W2 2
Dieses Problem wollen wir symbolisch, mit den Parametern m und l lösen.
Damit wir Körper der irrationalen Zahlen bleiben, wählen wir g = 10. Hier die
Maple-Sitzung:
> with(linalg):
M:=matrix([[m*l,0],[0,m*l]]);
ml
>
0
M :=
0 ml
> g:=10;
g := 10
> K:=matrix([[3*m*g,-m*g],[-m*g,m*g]]);
30 m −10 m
K :=
−10 m 10 m
> M1:=matrix([[sqrt(m*l),0],[0,sqrt(m*l)]]);
ml 0
M1 :=
0 ml
!
> M2:=inverse(M1);
""
ml
M2 :=
0
ml
ml
#
0
ml
> A:=evalm(M2&*K&*M2);
!
""
30 10
−
A := #
l l
10 10
−
l l
> V:=eigenvectors(A);
V :=
2+ 2 2− 2
[10 , 1, { −1 − 2, 1 }], [10 , 1, { 2 − 1, 1 }]
l l
> v1:=V[1];
2+ 2
v1 := [10 , 1, { −1 − 2, 1 }]
l
> v2:=v1[3];
3.1 Eigenbewegung 101
v2 := { −1 − 2, 1 }
v3:=v2[1];
>
v3 := −1 − 2, 1
> w1:=V[2];
2− 2
w1 := [10 , 1, { 2 − 1, 1 }]
l
w2:=w1[3];
>
w2 := { 2 − 1, 1 }
w3:=w2[1];
>
w3 := 2 − 1, 1
−1 −
> W:=matrix([v3,w3]);
2 1
W :=
2−1 1
!
> U:=evalm(M2&*W);
m l (−1 − 2)
""
ml
U := #
ml ml
m l ( 2 − 1) ml
ml ml
> r1:=U[1,1]/U[1,2];
r1 := −1 − 2
> r2:=U[2,1]/U[2,2];
r2 := 2−1
> rf1:=evalf(r1);
rf1 := −2.414213562
> rf2:=evalf(r2);
rf2 := .414213562
> node1:=0;
node1 := 0
> node2:=zero/1=(1-zero)/(-r1);
102 Zwei-Freiheitsgrad-Systeme
1 − zero
node2 := zero = −
−1 − 2
> solve(node2,zero);
1
2+ 2
> evalf(");
.2928932188
3.1.2.1 Modell-Gleichungen
abgekürzt
Mẍ + Cẋ + Kx = 0. (3.1.35)
1
0
0
1 x1 x2
0
1
k1 k2
0
1
0
1
0
1
m1 m2
0
1
0
1
0
1 c1 c2
Die Modalanalyse des Abschnitts 3.1.1.3 kann unmittelbar zur Lösung der Gl. (3.1.35)
verwendet werden, wenn proportionale Dämpfung
C = µM + κK (3.1.36)
3 Wir betrachten vorläufig nur die geschwindigkeits-proportionale Coulomb-Dämpfung.
3.1 Eigenbewegung 103
wird.
3.1.2.2 Modal-Analyse
zu
p̈ + (µI + κΛ)ṗ + Λp = 0. (3.1.39)
Komponentenweise geschrieben, sind dies entkoppelte Differentialgleichungen
der Form
p̈i (t) + 2ζi ωni ṗi (t) + ω2ni pi (t) = 0 (3.1.40)
mit ω2ni = λi und 2ζi ωni = µ + κω2ni oder
1
µ
ζi = + κωni . (3.1.41)
2 ωni
Gl. (3.1.40) ist ganz analog zu Gl. (2.2.9) gebaut, und demnach ist die Lösung
für den unterdämpften Fall 0 < ζi < 1
mit ṗ (0) + ζ ω 2
i i ni pi (0)
Ai = p2i (0) +
ωdi
und
pi (0)ωdi
φi = arctan
ṗi (0) + ζi ωni pi (0)
und mit
ωdi = 1 − ζ2i ωni .
Die zur Ermittlung von A und φ notwendigen Anfangswerte pi (0) und ṗi ste-
hen mit p(0) = V ¼ M1/2 x0 und ṗ(0) = V ¼ M1/2 ẋ0 zur Verfügung, und die
Rücktransformation der Lösung (3.1.42) geschieht wieder durch
Nachdem im ersten Abschnitt dieses Kapitels über Systeme mit zwei Freiheits-
graden die Eigenbewegung ohne und mit Dämpfung relativ ausführlich be-
handelt wurde, soll für die erzwungene Bewegung nur der Fall der Coulomb-
Dämfung beschrieben werden.
3.2.1 Bewegungsgleichungen
3.2.2 Modal-Analyse
mit
pi (0)ωdi
φi = arctan
ṗi (0) + ζi ωni pi (0)
und mii den Diagonalelementen der Massenmatrix M und pi (0) und ṗi (0) wie
in Paragraph 3.1.2.2. Auch die Rücktransformation in die Original-Koordinaten
x geschieht wie dort.
Für die Lösung (3.2.5) sind drei Integrale von Bedeutung, nämlich:
Integral:
t
sin(ωτ)eζi ωni τ sin 1 − ζ2i ωni (t − τ)dτ, (3.2.6)
0
Anregung: 1
Integral:
t
eζi ωni τ sin 1 − ζ2i ωni (t − τ)dτ, (3.2.7)
0
Integral:
t
sin(ωτ) sin ωni (t − τ)dτ, (3.2.8)
0
Die Lösungen der Integrale in der o.a. Reihenfolge präsentiert für Sie Maple:
> interface(labeling=false);
> int(sin(omega*tau)*exp(zeta[i]*omega[i]*tau)*
> sin(sqrt(1-zeta[i]^2)*omega[i]*(t-tau)),tau=0..t);
106 Zwei-Freiheitsgrad-Systeme
−e(t ζi ωi )
$
1 − ζi 2 ωi (2 ζi ωi cos(t ω) ω + sin(t ω) ω2 − sin(t ω) ωi 2 )
((ω + 2 ω 1 − ζ ω + ω ) (ω − 2 ω 1 − ζ ω + ω
2 2 2 2 2 2
)) + ω
i
i i i i i
(2 ζ ω cos( 1 − ζ ω t) 1 − ζ + 2 sin( 1 − ζ ω t) ζ
2 2 2 2 2
ωi 2
i i i
$ i i i i i
+ sin( 1 − ζ ω t) ω − sin( 1 − ζ ω t) ω ) (
2 2 2 2
i i
i i i
(ω2 + 2 ω 1 − ζi 2 ωi + ωi 2 ) (ω2 − 2 ω 1 − ζi 2 ωi + ωi 2 ))
> int(exp(zeta*omega[i]*tau)*sin(sqrt(1-zeta^2)
> *omega[i]*(t-tau)),tau=0..t);
1 − ζ2 e(t ζ ωi ) 1 − ζ2 cos( 1 − ζ2 ωi t) + ζ sin( 1 − ζ2 ωi t)
−
ωi ωi
> int(sin(omega*tau)*sin(omega[i]*(t-tau)),tau=
> 0..t);
sin(t ω) ωi sin(ωi t) ω
− +
(ω + ωi ) (ω − ωi ) (ω + ωi ) (ω − ωi )
3.2.4 Beispiele
3.2.4.1 KFZ-Aufhängung
In Abb. 3.8 ist eine einfache KFZ-Aufhängung und ihr mechanisches Ersatz-
modell dargestellt. Es sollen zuerst die Bewegungsgleichungen und daran an-
schließend die natürlichen Frequenzen des Systems bestimmt werden. Die Fahr-
zeugmasse sei munnamed.fig1 = 2000 kg, die Federsteifigkeit k1 = 103 N/m,
die Reifenmasse m2 = 50 kg und die Reifensteifigkeit k2 = 104 N/m.
Die natürlichen Frequenzen ergeben sich aus der Lösung der Polynomgleichung
det(−ω2n M + K) = 0
oder
−2000ω + 1000
−1000 2
n −1000
−50ω2n + 11000
= 105 ω4n − 2.205
107 ω2n + 107 = 0
und damit ω2n1,2 = λ1,2 = {0.4545, 220.0455} oder ωn1,2 = {0.6741, 14.8339}.
x1
m m1
1
k1 k1
x2
m2
m2
k2
1111111111
0000000000 11111111111
0
1
00000000000
k2
0000000000
1111111111 0
1
00000000000
11111111111
Abbildung 3.8: Prinzipskizze und Ersatzmodell einer KFZ-Aufhängung
Nun sei parallel zu den Federn je ein Dämpfungsglied mit den Dämpfungsve-
rhältnissen ζ1,2 = {0.01, 0.2} geschaltet. Außerdem sei das Rad von einer Kraft
f2 (t) = 10 sin 3t N angeregt. Es soll die Systemantwort x1 (t) der Fahrzeugmasse
bestimmt werden. Die Anfangsbedingungen seien Null.
M1 =
44.72135954999580 0
0 7.07106781186548
>> M2=inv(M1)
M2 =
0.02236067977500 0
0 0.14142135623731
>> A=M2*K*M2
1.0e+02 *
0.00500000000000 -0.03162277660168
-0.03162277660168 2.20000000000000
>> [V,D]=eig(A)
V =
0.99989628223203 -0.01440224907975
0.01440224907975 0.99989628223203
D =
1.0e+02 *
0.00454451365276 0
0 2.20045548634724
Die Koordinatentransformation nach p = V ¼ q führt auf die (zweimal) transfor-
mierte Systemgleichung
p̈(t) + diag(2ζωni )ṗ(t) + Λp(t) = V ¼ M−1/2 f(t),
deren rechte Seite
f(t) = (0.02036785597726 sin 3t, 1.41406688329898 sin 3t) ¼
ist. Die entkoppelten Gleichungen lauten somit
p̈1 + 0.0134826ṗ1 + 0.4544514p1 = 0.0203679 sin 3t
p̈2 + 5.9335729ṗ2 + 220.0455486p2 = 1.4140669 sin 3t,
> interface(labeling=false):
> dsolve({diff(p(t),t,t)+a*diff(p(t),t)+b*p(t)
> =c*sin(3*t),p(0)=0,D(p)(0)=0},p(t)):
> simplify(");
1 Ô
p(t) = (−6 e(−1/2 (a− a −4 b) t) b + 2 sin(3 t) b a2 − 4 b
2
2 Ô
+ 6 e(−1/2 (a+ a −4 b) t) b − 6 cos(3 t) a a2 − 4 b
2
Ô Ô
+ 54 e(−1/2 (a− a2 −4 b) t)
+ 3 e(−1/2 (a− a2 −4 b) t)
a2 − 4 b a
Ô Ô
(−1/2 (a+ a2 −4 b) t) (−1/2 (a+ a2 −4 b) t)
− 3e a + 3e
2
a2 − 4 b a
Ô Ô
Æ
(−1/2 (a− a2 −4 b) t) (−1/2 (a+ a2 −4 b) t)
+ 3e a − 54 e
2
Da nur die erste Koordinate x1 (t) gesucht ist, ist dies ein einfaches Skalarpro-
dukt, dessen Ausführung sei dem Leser überlassen.
3.2.4.2 KFZ-Nickbewegung
Viel-Freiheitsgrad-Systeme
Systeme mit mehr als zwei Freiheitsgraden sind mit Papier und Bleistift“mühsam
”
zu analysieren. Die Verallgemeinerung von zwei auf mehr als zwei Freiheitsgrade
unter Verwendung des Matrix-Kalküls und eines Computers (Matlab) ist dage-
gen trivial.
Die entsprechende Ausarbeitung erfolgt erst in den nächsten Semesterferien.
Diesmal ist die Materie kein Prüfungsstoff!
112 Viel-Freiheitsgrad-Systeme
Kapitel 5
Kontinuierliche Systeme
Nun sollen Systeme beschrieben und behandelt werden, die nicht mehr durch
diskret verteilte Massen sondern durch kontinuierlich verteilte Parameter cha-
rakterisiert sind. Dazu gehören Saiten, Stäbe, Balken, Platten, Schalen und
andere mechanische Systeme. Diese werden nun als partielle Differentialglei-
chungen modelliert: Die Variablen sind nun von Ort und Zeit abhängig, z.B.
bezeichnet w(x, t) die Durchbiegung eines Balkens am Ort x zum Zeitpunkt t.
Systeme zweiter Ordnung sind schwingende Saiten und längs- und torsionssch-
wingende Stäbe. Zunächst seien schwingende Saiten (Streichinstrumente, Brü-
ckenkabel, Überlandleitungen) behandelt.
In Abbildung 5.1 ist eine beidseitig eingespannte Saite mit der Länge l gezeigt.
Sie steht unter (Längs-)Spannung τ, hat eine homogene Dichteverteilung ρ und
wird quer zur Längenabmessung mit Last/Länge f(x, t) belegt. Die seitliche
Auslenkung der Saite ist w(x, t).
114 Kontinuierliche Systeme
11
00 f(x,t) 111
000
00
11
00
11
000
111
000
111
00
11
00
11
000
111
000
111
00
11
000
11
000
111
l111
000
00
11 x
000
111
f(x,t) w(x,t) τ2
θ2
x2 =x1 + ∆ x
θ1 x
τ1 1
Die Abbildung zeigt auch einen infinitesimalen Ausschnitt ∆x, anhand dessen
die Bewegungsgleichung modelliert wird. Das zweite Axiom von Newton ergibt
quer zur Saite:
∂2 w(x, t)
ρ∆x = f(x, t)∆x − τ1 sin θ1 + τ2 sin θ2 (5.1.1)
∂t2
und längs zur Saite
0 = −τ1 cos θ1 + τ2 cos θ2 . (5.1.2)
Mit der Approximation für kleine Winkel θ und der Steigung ∂w/∂x der Saite
∼ tan θi = ∂w(x, t)
sin θi =
∂x
xi
Bei Vernachlässigung der Terme höherer Ordnung O(∆x2 ) wird Gl. (5.1.3) mit
c = τ/ρ, der Ausbreitungsgeschwindigkeit der Welle, und f(x, t) = 0 die freie
Schwingung der Saite
∂2 w(x, t) ∂2 w(x, t)
2
= c2 . (5.1.4)
∂t ∂x2
5.1 Systeme zweiter Ordnung 115
Gleichung (5.1.4) ist die Gleichung der schwingenden Saite; sie wird in der ma-
thematischen Physik als Wellengleichung bezeichnet. Es ist eine lineare par-
tielle Differentialgleichung zweiter Ordnung in Raum und Zeit. Zu ihrer Lösung
werden zwei (räumliche) Randbedingungen und zwei (zeitliche) Anfangsbedin-
gungen benötigt: Man spricht von einem Anfangs-Randwert-Problem. Bei der
eingespannten Saite liegen beide Ränder fest:
Ganz analog zu den Systemen mit endlich vielen Freiheitsgraden sind Anfang-
swerte für den Ort und die Geschwindigkeit anzugeben:
Die Saiten- oder Wellengleichung ist eine Gleichung vom hyperbolischen Typ.
Bezüglich der Klassifizierung von partiellen Differentialgleichungen geben An-
hang E und Logan[30] Auskunft.
∂2 w(x, t) 2
2 ∂ w(x, t)
= c (5.1.7)
∂t2 ∂x2
und separieren die Variable w(x, t) in einen ortsabhängigen Anteil X(x) und
einen zeitabhängigen Anteil T (t):1
Dies bedeutet, daß X(x) die Form der Saite bestimmt und T (t) ihre Bewegung.
Damit wird Gl. (5.1.7)
T ¼¼ + λc2 T = 0. (5.1.11)
Wir lösen zunächst das Ortsproblem. Da w(0, t) = X(0)T (t) = 0 und w(l, t) =
X(l)T (t) = 0 sind, schließen wir auf X(0) = X(l) = 0, da die alternative An-
nahme T (t) = 0 zu w(x, t) = 0 also zur unbewegten Saite führt. Damit wird
Problem (5.1.10) zu einem Randwertproblem
Hierin ist λ ganz analog zu diskreten Systemen Eigenwert des Problems und
nichttriviale Lösungen X die zugehörigen Eigenfunktionen.
Wir untersuchen unterschiedliche Lösungstypen: a) Für λ = 0 ist X¼¼ = 0, was
nach zweifacher Integration zu X(x) = ax + b mit Konstanten a und b führt;
diese sind wegen X(0) = b = 0 und X(l) = al = 0 beide Null, weswegen
Problem (5.1.12) keine nichttriviale Lösung besitzt. b) Nun sei λ < 0, z.B.
λ = −k2 , was zu X¼¼ − k2 X = 0 führt mit der Lösung X(x) = aekx + be−kx .
Die Randwerte ergeben zunächst X(0) = a + b = 0 oder a = −b mit X(x) =
2a sinh kx und dann speziell X(l) = 2a sinh kl = 0, was nur für a = 0 möglich
ist: Auch dieser Fall führt nur auf triviale X(x) 0. Es bleibt c) λ > 0, z.B.
λ = k2 und X¼¼ + k2 X = 0 mit der Lösung
Der linke“ Randwert fordert X(0) = a = 0 und damit X(x) = b sin kx. Der
”
rechte“ Randwert bedingt X(l) = b sin kl = 0, was einerseits für b = 0 wiede-
”
rum nur die triviale Lösung X(x) 0 ergibt, andererseits auch für sin kl = 0
erfüllt ist, was in
kl = nπ, n = 1, 2, 3, . . . ,
oder
λn = k2 = n2 π2 /l2 , n = 1, 2, 3, . . . , (5.1.13)
resultiert mit entsprechenden Eigenfunktionen
x
Xn (x) = sin nπ . (5.1.14)
l
Die erhaltenen Eigenwerte werden nun in die Zeitgleichung eingesetzt:
n2 π2 c2
T ¼¼ + T = 0,
l2
5.1 Systeme zweiter Ordnung 117
Orts- und Zeitlösung werden zur Gesamtlösung zusammengeführt:
ct ct
x
wn (x, t) = Xn (x)Tn (t) = an cos nπ + bn sin nπ . sin nπ
l ll
(5.1.15)
Diese Lösung erfüllt die Wellengleichung für alle positiven n und die Randbe-
dingungen w(0, t) = w(l, t) = 0. Die noch unbestimmten Koeffizienten an und
bn in Gl. (5.1.15) werden durch Auswertung der Anfangsbedingungen w0 (x)
und v0 (x) gewonnen. Hier wird nur das Ergebnis mitgeteilt, eine ausführliche
Begründung findet der interessierte Leser in O’Neil[37].
∞
Auch die unendliche Summe
ct
ct
x
w(x, t) = an cos nπ + bn sin nπ sin nπ (5.1.16)
l l l
n=1
erfüllt. Ebenso wird mit der partiellen Ableitung der Lösung nach der Zeit
∞
c ct ct x
wt (x, t) = nπ −an sin nπ + bn cos nπ sin nπ
l l l l
n=1
der Anfangswert
∞
c x
wt (x, 0) = nπ bn sin nπ = v0 (x)
l l
n=1
Randbedingungen:
w(0, t) = 0, w(1, t) = 0, t 0,
Anfangsbedingungen:
⎧
⎪
⎪ falls0 x 1/4,
⎨x
w(x, 0) = w0 (x) = 1/4 falls1/4 x 3/4, ẇ(x, 0) = v0 (x) = 0, 0 x 1.
⎪
⎪
⎩1 − x falls3/4 x 1,
0.2
0.1
–0.1
–0.2
0
0.2 0
0.4 0.2
0.4
t 0.6
0.6 x
0.8 0.8
1
5.1 Systeme zweiter Ordnung 119
X (0)
die beiden Randbedingungen in Vektordarstellung formuliert:
¼ −ka0 kb1
0
= = .
X¼ (l) −ka sin kl kb cos kl 0
Da aus der ersten Gleichung b = 0 folgt, kann die Vektordarstellung umgeformt
werden in eine Matrixdarstellung:
0 k
a 0
= .
−k sin kl k cos kl b 0
Dieses lineare Gleichungssystem mit den unbekannten Koeffizienten a und b hat
nur dann eine nicht verschwindende Lösung, wenn die Determinante der Matrix
Null wird, d.h.
k2 sin kl = 0. (5.1.20)
Gleichung (5.1.20) ist die charakteristische Gleichung des Systems.
Aufgabe 21 Geben Sie die Lösung der Wellengleichung mit den Randbedingun-
gen dieses Abschnitts an.
2
120 Kontinuierliche Systeme
Nach der schwingenden Saite betrachten wir als nächstes Kontinuum die Läng-
sschwingungen von Stäben. Auch sie führen auf das Modell einer Differential-
gleichung zweiter Ordnung im Ort und in der Zeit.
Eine Prinzipskizze der Anordnung zeigt Abb. 5.3, die auch ein infinitesimales
Stabelement mit der Längenausdehnung dx enthält,2 auf das das zweite Axiom
von Newton angewendet, folgende Gleichgewichtsbeziehung ergibt:
∂2 w(x, t)
ρA(x)dx = F + dF − F.
∂t2
Hierin ist w(x, t) die Auslenkung des Stabes am Ort x zum Zeitpunkt t, ρ die
Dichte des Stabmaterials, A(x) die (variable) Querschnittsfläche, und F bzw. F+
dF sind die (Zug-)Kräfte, die in x−Richtung auf die beiden Querschnittsflächen
wirken.
11
00
00
11
00
11
00
11
00
11
00
11
00
11 0
00
11
l
11
00 x
w(x,t)
F F+dF
∂w(x, t)
F = EA(x) .
∂x
In die Gleichgewichtsbeziehung eingesetzt und mit dF = (∂F/∂x)dx erhält man
∂2 w(x, t) ∂
∂w(x, t)
ρA(x) 2
= EA(x) .
∂t ∂x ∂x
2 Das durchgezogene Element kennzeichnet den Gleichgewichtszustand, das gestrichelte den
deformierten.
5.1 Systeme zweiter Ordnung 121
∂2 w(x, t) E ∂2 w(x, t)
2
= . (5.1.21)
∂t ρ ∂x2
Die Lösung der Modellgleichung (5.1.21) wird wieder durch Separation der Va-
riablen w(x, t) = X(x)T (t) gewonnen. Sie hängt von den Randwerten
und Anfangswerten
ab. Wir besprechen zwei Beispiele und fassen Resultate für andere gebräuchliche
Randwerte in einer Tabelle zusammen.
wird differenziert
X¼ (x) = −ak sin kx + bk cos kx.
Aus w(0, t) = X(0)T (t) = 0 folgt X(0) = 0 und damit a = 0. Die charakteris-
tische Gleichung folgt aus
1 0
a 0
=
−k sin kl k cos kl b 0
und 1
−k sin kl 0
k cos kl
=0
zu cos kl = 0 und kl = 2(n + 1)π/2. Hieraus folgen die Eigenwerte
2(n + 1)π 2
λn = k2n = , n = 1, 2, 3, . . . .
2
122 Kontinuierliche Systeme
mit
2
l 2(n + 1)π x
cn = w0 (x) sin dx
l 0 2 l
und
4
l 2(n + 1)π x
dn = v0 (x) sin dx
(2n + 1)πc 0 2 l
2
Beispiel 25 Stab mit Endmasse
In diesem Beispiel trage der eingespannte Stab an seinem freien Ende eine Zu-
satzmasse m1 .
111
000
000
111
0000
111
000
111
l
000
111
000
111 ρ AE m1
000
111
000
111
000
111
x
Der Randwert am linken“ Rand ist wieder w(0, t) = 0, und am rechten“ Rand
” ”
ergibt er sich aus der Gleichgewichtsbedingung zwischen Zugkraft am Balkenende
und Trägheitskraft der Endmasse.
∂w(l, t) ∂2 w(l, t)
AE = −m1 .
∂x ∂t2
Die partiellen Ableitungen der Auslenkung w(x, t) erhält man aus dem Lösung-
sansatz
zu
∂wn (x, t)
= kn (−an sin kn x + bn cos kn x)(cn cos kn ct + dn sin kn ct)
∂x
5.1 Systeme zweiter Ordnung 123
bzw.
∂2 wn (x, t)
= −k2n c2 (an cos kn x + bn sin kn x)(cn cos kn ct + dn sin kn ct).
∂t2
Setzt man diese in die Randwertbedingung ein und berücksichtigt, daß an = 0
ist (was aus dem linken Randwert folgt), dann ergibt sich
AE cos kn l = m1 kn c2 sin kn l
oder
AEl m
kn l tan kn l = = .
c2 m1 m1
Dies ist eine transzendente Gleichung zur Bestimmung der Eigenfrequenzen kn .
Sie muß numerisch gelöst werden. Die rechte Seite ist das Verhältnis der Stab-
masse m zur Endmasse m1 .
11
00
00
11 00
11
11
00
00
11 00
11
0 l
00
11 00
11
00
11
00
11
ρ AE m1
00
11
00
11
00
11 x k 00
11
00
11 00
11
Abbildung 5.5: Schwingender Stab mit Endmasse und Federabstützung
EAkn
tan kn l = − (5.1.22)
m1 k2n c2 + k
Aufgabe 23 Beschreiben Sie den Unterschied zwischen Zeile 1 und 3 der Tabelle.
Erklären Sie die physikalische Bedeutung.
2
124 Kontinuierliche Systeme
nπc x
sin kn l = 0 kn = cos nπ
2n + 1 x
l l
11
00
00
11 (2n + 1)πc
00
11 cos kn l = 0 kn =
00
11 x sin π
00
11 2l 2 l
00
11
00
11
00
11
00
11 11
00
00
11 00
11
00
11 00
11 nπc x
00
11 x 00
11 sin kn l = 0 kn = sin nπ
00
11 00
11 l l
00
11 00
11
00
11 00
11
00
11 00
11
00
11 111
000
00
11 000
111
00
11 000
111 kn cot kn =
kl
kn =
nπc x
00
11 x 000
111 sin nπ
00
11 000
111 EA l l
00
11 000
111
00
11 000
111
00
11 000
111
Tabelle 5.1: Alternative Randwerte
Das mathematische Modell der Torsionsschwingung von Stäben wird unter Be-
zug auf Abb. 5.6 gewonnen. Zwischen dem Torsionsmoment T (x, t) und dem
Torsionswinkel θ(x, t) besteht folgender Zusammenhang[20]:
∂
T (x, t) = GJ(x) θ(x, t),
∂x
mit G dem Schermodul des Stabmaterials und J(x) dem polaren Trägheitsmo-
ment der Querschnittsfläche des Stabes.
Betrachtet werde ein infinitesimal kleiner Ausschnitt des Stabes der Länge dx,
auf den ein äußeres Moment f(x, t) pro Längeneinheit wirke; Newtons zweites
Gesetz ergibt
∂2
J0 (x)dx θ(x, t) = T (x, t) + dT (x, t) − T (x, t) + f(x, t)dx
∂t2
mit J0 (x) dem polaren Trägheitsmoment pro Längeneinheit dx. Mit dT (x, t) =
∂x T (x, t)dx erhalten wir
∂
∂2 ∂
∂
J0 (x)dx 2 θ(x, t) = GJ(x) θ(x, t) dx + f(x, t)dx,
∂t ∂x ∂x
und für homogenes Material mit konstantem Querschnitt, also J0 = ρJ, mit ρ
der Materialdichte:
∂2 ∂2
ρ 2
θ(x, t) = G 2 θ(x, t) + f(x, t). (5.1.23)
∂t ∂x
5.2 Systeme vierter Ordnung 125
dx
θ
l
x
f(x,t) T
θ
T+dT
θ+d θ
x
dx
x+dx
∂2 2 ∂
2
θ(x, t) = c θ(x, t). (5.1.24)
∂t2 ∂x2
Diese Beziehung ist strukturidentisch zur Saitengleichung (5.1.4) und zur Glei-
chung für die Longitudinalschwingung von Stäben (5.1.21).
Zur Modellbildung sei folgende vereinfachende Annahme getroffen, die mit guter
Näherung für schlanke Balken (l/b, l/h > 10) zutrifft[2, Kap. 11].
Annahme: Die Scher-Verformung sei klein gegenüber der Biegung w(x, t), d.h.
die Seiten des Elements dx verformen sich nicht und bleiben senkrecht zur Bie-
gelinie.
126 Kontinuierliche Systeme
f(x,t)
z
w(x,t)
x
dx
x
l
f(x,t)
M(x,t)+dM(x,t)
M(x,t)
w(x,t) V(x,t)+dV(x,t)
V(x,t)
dx
x
∂2
ρA(x)dx w(x, t) = V(x, t) − (V(x, t) + dV(x, t)) + f(x, t)dx,
∂t2
oder mit
∂V(x, t)
dV(x, t) = dx
∂x
∂2 ∂
ρA(x) w(x, t) + V(x, t) = f(x, t). (5.2.1)
∂t2 ∂x
Ebenso gilt für die Momente um O:3
1
−M(x, t) + (M(x, t) + dM(x, t)) − (V(x, t) + dV(x, t))dx + dxf(x, t)dx = 0
2
oder mit
∂M(x, t)
dM(x, t) = dx
∂x
und unter Vernachlässigung von Termen höherer Ordnung
∂M(x, t)
V(x, t) = .
∂x
Einsetzen dieser Beziehung in Gl. (5.2.1) ergibt
∂2 ∂2
ρA(x) 2
w(x, t) + 2 M(x, t) = f(x, t).
∂t ∂x
Die Festigkeitslehre lehrt die Proportionalität von Biegemoment M(x, t) und
Durchbiegung w(x, t) mit dem Elastizitätsmodul E als Proportionalitätsfaktor:
∂2
M(x, t) = EI(x) w(x, t).
∂x2
Beide voranstehenden Beziehungen ergeben
∂2 ∂2
∂2
ρA(x) 2 w(x, t) + 2 EI(x) 2 w(x, t) = f(x, t).
∂t ∂x ∂x
∂2 ∂4
2
w(x, t) + c2 4 w(x, t) = 0, (5.2.2)
∂t ∂x
wobei
EI
c= (5.2.3)
ρA
ist.
3 Wegen Annahme I ist wirkt kein Massenträgheitsmoment.
128 Kontinuierliche Systeme
w(x, 0) = w0 (x),
und
ẇ(x, 0) = v0 (x),
die nicht beide identisch verschwinden dürfen.
Die Lösung der Modellgleichung (5.2.2) wird wieder durch Separation der Va-
riablen w(x, t) = X(x)T (t) gewonnen. Einsetzen dieses Ansatzes in Gl. (5.2.2)
führt auf
T ¼¼ (t)X(x) + c2 X¼¼¼¼ (x)T (t) = 0,
welche Beziehung nach Umformung auf
nur für eine Konstante ω2 erfüllt sein kann, wie in Abschnitt 5.1 ausführlich
dargelegt wurde.
5.2 Systeme vierter Ordnung 129
und
T ¼¼ (t) + ω2 T (t) = 0, (5.2.6)
wobei abkürzend
ω2 ρAω2
β4 = 2
= (5.2.7)
c EI
steht.
X(x) = Ceλx ,
(λ4 − β4 )Ceλx = 0
oder
λ4 = β4
mit den Lösungen
λ = β, jβ
führt. Die Lösung ist damit
woraus sich mit Gl. (5.2.8) folgendes homogene Gleichungssystem für die ci
ergibt:
! ! !
"" "" ""
0 1 0 1 c1 0
−1
# # #
0 0 1 c2 0
=
sin βl cos βl sinh βl cosh βl c3 0
− sin βl − cos βl sinh βl cosh βl c4 0
und da c2 + c4 = 0 und −c2 + c4 = 0 sind, vereinfacht es sich zu
sin βl sinh βl
c 0
1
= . (5.2.9)
− sin βl sinh βl c3 0
Diese Gleichung führt nur dann zu nichttrivialen Lösungen, wenn ihre Deter-
minante gleich Null ist:
βn l = nπ, n = 1, 2, . . .
und damit
n
βn = π, n = 1, 2, . . . . (5.2.10)
l
n EI
Damit werden die Eigenfrequenzen mit der Beziehung (5.2.7)
2
ωn = π , n = 1, 2, . . . . (5.2.11)
l ρA
2
Aufgabe 24 Wiederholen Sie Beispiel 26 für den einseitig und den beidseitig
eingespannten Balken.
2
5.2.1.2.2 Zeitkomponente Nach Kap. 2.1 ist die Lösung der Zeitkomponente,
T ¼¼ (t) + ω2 T (t) = 0,
hier durch
Tn (t) = an cos ωn t + bn sin ωn t (5.2.13)
formuliert, wobei die Koeffizienten an und bn aus den Anfangswerten bestimmt
werden, und die ωn während der Berechnung der Ortskomponente X(x) ange-
fallen sind.
5.2 Systeme vierter Ordnung 131
Beispiel 27 Wir führen Beispiel 26 fort und geben als w0 (x) = sin 2π xl und
ẇ(x) = v0 (x) = 0 an.
Die letzte Bedingung führt wegen
d
Tn (0) = −an
0 + bn
1 = 0
dt
auf bn = 0 und damit wird die Gesamtlösung
∞
w(x, t) = X(x)T (t) = [(a1 c1 )n cos ωn t + (a1 c1 )n sin ωn t] sin βn x,
n=1
wobei die (a1 c1 )n = αn neue Paramater sind, die aus der verbleibenden An-
fangsbedingung w(x) wie folgt berechnet werden.
Es ist
∞
x
w0 (x) = sin 2π = w(x, 0) = αn sin βn x,
l
n=1
fallen O−T und O−N nicht mehr zusammen: Es ergibt sich ein Scherwinkel γ =
φ−α. Die Scherung hat nur Verformung, keine Rotation zur Folge. Die Rotation
ist Konsequenz der Verformung in Zusammenhang mit der Durchbiegung des
Balkenelements.
Für das Biegemoment M und die Scherkraft V übernehmen wir folgende Zu-
sammenhänge aus der Festigkeitslehre[3]:
∂
M(x, t) = EI(x) φ(x, t) (5.2.14)
∂x
und
∂
V(x, t) = kA(x)Gγ = kA(x)G(φ(x, t) − w(x, t)) (5.2.15)
∂x
mit G dem Schermodul und k einem Formfaktor, der die Querschnittsform
berücksichtigt. Dieser Faktor wird heißt Timoshenko-Scher-Koeffizient und ist
k = 56 für Rechteckquerschnitte und k = 10
9
für Kreisquerschnitte[10].
y,w
M+dM
N
V+dV
P φ
γ
T
Q
M O α
O´
R
w(x,t)
V
dx
Die Newton- und Eulergleichungen der Translation und Rotation des Balkene-
lement können nun wie folgt formuliert werden:
1. Translation in y−Richtung
∂2
ρA(x)dx w(x, t) = V(x, t) − (V(x, t) + dV(x, t)) + f(x, t)dx,
∂t2
5.2 Systeme vierter Ordnung 133
∂2 dx
ρI(x)dx φ(x, t) = −M(x, t)+(M(x, t)+dM(x, t))+(V(x, t)+dV(x, t))dx+f(x, t)dx .
∂t2 2
Einsetzen der Ausdrücke für Biegemoment und Scherkraft in diese beiden Glei-
chungen und Vernachlässigung von Termen höherer Ordnung in dx führt auf
∂2 ∂
∂
ρA(x) w(x, t) = − kA(x)G φ(x, t) − w(x, t) + f(x, t)
∂t2 ∂x ∂x
und
∂2 ∂
∂
∂
ρI(x) 2 φ(x, t) = EI(x) − kA(x)G φ(x, t) − w(x, t) .
∂t ∂x ∂x ∂x
∂
Für konstanten Querschnitt kann die erste dieser Gleichungen nach ∂x φ(x, t)
aufgelöst und in die zweite eingesetzt werden:
∂4
E
∂4 ρ2 I ∂4 ∂2
EI w(x, t) − ρ 1 + w(x, t) + w(x, t) + ρA 2 w(x, t)
∂x4 kG 2
∂x ∂t 2 kG ∂t 4 ∂t
EI ∂2 ρI ∂2
= f(x, t) + 2
f(x, t) − f(x, t),
kAG ∂x kAG ∂t2
und für den homogenen Fall f(x, t) = 0:
∂4
E
∂ w(x, t)
4
ρ2 I ∂4 ∂2
EI 4 w(x, t) − ρ 1 + + w(x, t) + ρA w(x, t) = 0.
∂x kG ∂x2 ∂t2 kG ∂t4 ∂t2
(5.2.16)
134 Kontinuierliche Systeme
Kapitel 6
Technische Anwendungen
6.1 Schwingungsdämpfung
6.2 Schwingungstilgung
Zunächst sei keine Dämpfung im System vorausgesetzt. Die Anregung der Primär-
masse sei f(t) = F0 sin ωt, dann gelten folgende Systemgleichungen:
m 1 0
ẍ k
1 1 + k2 −k2
x F
1 0 sin ωt
+ = , (6.2.1)
0 m2 ẍ2 −k2 k2 x2 0
Maschine
F sin( ω t)
0
m1
x1
k2
x2
k_1 k
_1
m2
2 2
Schwingungstilger
111111111111111111111
000000000000000000000
Abbildung 6.1: Schwingungstilgung
k2 − m2 ω2
X1 = F0 = 0, (6.2.4)
(k1 + k2 − m1 ω2 )(k2 − m2 ω2 ) − k22
woraus die Abstimmungsbedingung für den Isolator
k2
= ω2 (6.2.5)
m2
folgt. Diese Bedingung soll insbesondere die Vibrationen nahe der Systemreso-
nanz
k1
ω211 = ω
m1
absorbieren,1 was zusammen mit Gl. (6.2.5) zu
k1 k2
=
m1 m2
führt. Mit diesen Beziehungen und den Definitionen
ω ω
η1 = und η2 =
ω11 ω22
1 Beachten Sie, daß ω11 und ω22 die Eigenfrequenzen der beiden getrennten Systeme
(Maschine und Absorber) sind.
6.2 Schwingungstilgung 137
und der statischen Auslenkung ∆ = F0 /k1 wird die Lösung (6.2.3) als
X1
∆
= k2
1 − η21
k2
(6.2.6a)
1+ − η21 1 − η22 −
k1 k1
X2
∆
= k2
1 − k 2
(6.2.6b)
1+ − η21 1− η22
k1 k1
Der Absolutwert von (6.2.6a) ist für µ = m1 /m2 = 1/4 in Abb. 6.2 dargestellt.
1.8
µ=1/4
ω1=ω2
1.6
← ohne Isolator
1.4
1.2
|X1/∆|
0.8
0.6
0.4
0.2
a b
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
η
F2 = k2 X2 = −F0 .
Die beiden Eigenwerte der isolierten Maschine ergeben sich als Nullstelle der
Determinante det(A) mit folgenden Definitionen:
m2 ω22
µ= und ν =
m1 ω11
ω11
ω1,2 = 1 + ν2 (1 + µ) ν4 (1 + µ)2 − 2(1 − µ)ν2 + 1. (6.2.7)
2
In Abb. 6.3 sind die Eigenwerte als Funktion der Isolatormasse für unterschied-
liche Parameter ν dargestellt.
ν=1/2
2.5
2
ω2/ω11
ν=1
1.5
ω1/ω11
ν=2
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
µ
2. Die Vibrationen der Maschine werden nur für eine bestimmte Frequenz
getilgt; für andere (Betriebs-)Frequenzen können die Amplituden durchaus
störend hoch werden.
6.2 Schwingungstilgung 139
Nach Einführung einer viskosen Dämpfung entsprechend Abb. 6.4 wird Gl. (6.2.1)
Maschine
F sin( ω t)
0
m1
x1
k2 c2
x2
k_1 k
_1
m2
2 2
1
0
0
1 Schwingungstilger
1111111111111111111111
0000000000000000000000
0000000000000000000000
1111111111111111111111
Abbildung 6.4: Gedämpfte Schwingungstilgung
m ẍ c ẋ k x F
folgendermaßen ergänzt:
1 0 1 1+ c2 −c2 1 1 + k2 −k2 1 0
+ + = ejωt ,
0 m2 ẍ2 −c2 c2 ẋ2 −k2 k2 x2 0
(6.2.8)
oder allgemein:
Mẍ + Cẋ + Kx = Fejωt , (6.2.9)
in die nun wegen der durch die Dämpfung eingeführten Phasenverschiebung
zwischen Eingang und Ausgang der Lösungsansatz
xi = Xi ejωt , i = 1, 2, (6.2.10)
eingesetzt wird. Dann ergibt sich folgendes Gleichungssystem für die Amplituden
k X F
Xi :
X k F
und daraus
1 1 2 − m2 ω2 + jc2 ω k2 + jc2 ω 0
= ,
X2 det(A) k2 + jc2 ω k1 + k2 − m1 ω2 + j(c1 + c2 )ω 0
(6.2.12)
mit A = K − ω2 M + jωC und det(A)
k2 − m2 ω2 + jc2 ω
X1 = F0 , (6.2.14a)
|A|
k2 + jc2 ω
X2 = F0 . (6.2.14b)
|A|
Der Absolutwert der auf die statische Auslenkung ∆1 = F0 /k1 normierten Am-
plituden Xi wird für den wichtigen Fall c1 = 0 mit
m2 ω22 c2 ω
µ= , ν= , ζ= , η1 =
m1 ω11 2m2 ω1 ω11
berechnet als
X
∆ 1 = (2ζη ) (−1 + (1 +(2ζη
2
) + (η − ν )
)
1
2 2
1
µ η ) + (µν η + (1 − η )(η
2 2 2 2
2 2
2 2
− ν2 ))2
,
1 1
1 1 1 1
X (6.2.15a)
∆
2 = (2ζη ) (−1 + (1 + µ η(2ζη
2
) +ν
) 2 2
1
2
) + (µν η + (1 − η )(η
4
− ν2 ))2
2 2
. 2 2
1 1 1 1 1 1
(6.2.15b)
Die Maschinenauslenkung |X1 /∆1 | ist in Abb. 6.5 für µ = 1/4, ν = 1, η1 = 1 und
verschiedene ζ dargestellt. Es ist zu erkennen, was schon Gl. (6.2.14a) zeigt, daß
für kein ω die Maschinenamplitude X1 mehr Null ist. Dafür ist aber gegenüber
dem ungedämpften Fall in Abb. 6.2 der Bereich [a, b], in welchem X1 /∆1 1
gilt, aufgeweitet, d.h. das System ist unempfindlicher gegenüber Änderungen
der nominellen Drehzahl.
6.2.2.1 Optimierung
12
10
←ζ=1/10
8
|X1/∆|
←ζ=1/2
←ζ=1/4
2
←ζ=1/6
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
η
genügen.
Das zu vermessende System wird nicht mit einer Probe f(t) angeregt, sondern
es wird eine bestimmte, statistisch signifikante Anzahl von Experimenten durch-
geführt, und deren Resultate x(t) zu einem Ensemble zusammengefaßt. Wichtige
Zufallsfunktionen sind solche, deren statistische Eigenschaften nicht von der Zeit
t abhängen; man nennt sie stationäre Zufallsfunktionen. Z.B. ist ihr Mittelwert
x̄
1 T
x̄ = lim x(t)dt (6.3.1)
T→∞ T 0
ist ein Maß für die Schwankungen des Signals x(t). Die Quadratwurzel aus der
Varianz ist2
xrms = x̄2 .
Die Autokorrelationsfunktion Rxx (τ) ist definiert als
1 T
Rxx (τ) = lim x(t)x(t + τ)dt. (6.3.3)
T→∞ T 0
Sie ist ein Maß dafür, wie schnell sich das Signal x(t) mit der Zeit ändert. Die
Fourier-Transformierte der Autokorrelationsfunktion ist die spektrale Leistungs-
dichte Sxx (ω)
1 ∞
Sxx (ω) = Rxx (τ)e−jωτ dxτ (6.3.4)
2π −∞
oder rücktransformiert:3
∞
Rxx (τ) = Sxx (ω)ejωτ dxω. (6.3.5)
−∞
und
Sfx (ω) = S£xf (ω), (6.3.9)
mit S£xf (ω) der konjugiert-komplexen Funktion Sfx (ω).
Kreuzkorrelationen beinhalten Informationen über die Abhängigkeit unterschied-
licher Signale. Sie geben z.B. an, ob eine Schwingung an einem bestimmten Ort
einer Struktur eine Schwingung an einem anderen Ort beinflußt. Ist die Kreuz-
korrelation für einen bestimmten Wert von τ hoch, dann ist der wechselseitige
Einfluß groß, und die gegenseitige Übertragungszeit beträgt τ.
Unter Zuhilfenahme längerer Umrechnungen von Fourier-Integralen können für
gedämpfte mechanische Systeme folgende Beziehungen gezeigt werden[32, Seite
26–29]:
Sfx (ω) = H(ω)Sff (ω) (6.3.10)
und
Sxx (ω) = H(ω)Sxf (ω) (6.3.11)
2 Sie heißt im Englischen root-mean-square value.
3 Gl. (6.3.4) und Gl. (6.3.5) sind die Beziehuungen von Wiener und Khintchine.
6.3 Experimentelle Modalanalyse 143
und
Sxx (ω) = |H(ω)|2 Sff (ω). (6.3.12)
Damit wird die spektrale Leistungsdichte des Systemausgangs Sxx mit der des
Eingangs Sff über das Quadrat des Betrags des Frequenzgangs4 H(ω) verbun-
den
Da sowohl Gl. (6.3.10) als auch Gl. (6.3.11) eine Basis zur Berechnung des Fre-
quenzganges H(ω) bilden, kann sowohl die eine als auch die andere zur Über-
prüfung von H herangezogen werden. Am einfachsten geschieht dies durch Aus-
wertung der Kohärenzfunktion γ2
H1 (ω) Sfx (ω) Sxf (ω) Sfx (ω)S£fx (ω)
=
=
H2 (ω) Sff (ω) Sxx (ω) Sff (ω)Sxx (ω)
|Sfx (ω)|2
= =: γ2 (ω). (6.3.13)
Sxx (ω)Sff (ω)
Es gilt
0 < γ2 1,
und die Messungen sind kohärent für γ2 = 1, d.h. sie sind nicht durch Meßrau-
schen kontaminiert. Bei γ2 = 0 haben Sie reines Rauschen gemessen.
1
0
0
1
0
1
0
1
Beschleunigungsaufnehmer 6
0
1
5
0
1
4
0
1
elastische Struktur
2
0
1
1
0
1
0
0 0.5 1 1.5 2 2.5 3
0
1 Shaker Anregung
Analysator
Die i−te Eigenfrequenz ωi wird dort angenommen, wo das i−te Maximum des
Amplitudengangs auf der Frequenzabszisse gefunden wird, s.a. Abschnitt 2.3.2:
ηmax = 1 − 2ζ2 (6.3.14)
und
1
Vmax = . (6.3.15)
2ζ 1 − ζ2
Das Dämpfungsverhältnis ζi wird aus der Bandbreite ∆ωi = 2ζi ωi zu
∆ωi
ζi = (6.3.16)
2ωi
berechnet. Dabei wird die Bandbreite als Frequenz-Differenz der Schnittpunkte
der Geraden Ô12 Vmax mit dem Graphen des Amplitudengangs gewonnen, s.a.
Abb. 2.12 auf Seite 39.
6.3.2.2 Eigenvektoren
1. Transformation:
2. Diagonalmatrizen:
und
Λ = diag(ω2i ) = V T M−1/2 KM−1/2 V,
3. Einsetzen von
C = M1/2 V∆V T M1/2
und
K = M1/2 VΛV T M1/2
in Gl. (6.3.19) mit Q = M1/2 V
4. Rezeptanz-Matrix:
= Q−T diag(ω2i − ω2 + 2jζi ωi ω)−1 Q−1
1
= Q−T diag Q−1
ωi − ω + 2jζi ωi ω
2 2
n
ui uTi
= , (6.3.20)
i=1
ω2i − ω2 + 2jζi ωi ω
Uk
Ê
αkl (ω) = Hkl (jω) = . (6.3.21)
Fl
Für ω = ωi erhält man aus Gl. (6.3.20) für den Betrag des Elementes {kl} der
Frequenzgangsmatrix die Approximation6
Ê
und damit
|ui uTi |kl = 2ζi ω2i |Hkl (ωi )|. (6.3.23)
Das Vorzeichen von |ui uTi |kl wird demPhasendiagramm H(ωi ) entnommen.
5 Die
6 Wegen
restlichen Eingänge fi , i = 1, . .. , n, i = l werden Null gehalten.
der Resonanzüberhöhung bei ω = ωi sind die restlichen Summanden aus (6.3.20)
klein gegenüber (6.3.22).
146 Technische Anwendungen
6.3.2.2.1 Betrag der Eigenvektorelemente Wie werden nun aus der gemesse-
nen Rezeptanz-Matrix die Eigenvektoren ui , i = 1, . . . , n, gewonnen? Darüber
geben die Matrizen ui uTi Auskunft: Dies sind Dyaden der Form
u 2
u u
i,1 ui,1 ui,2 ui,1 ui,n
ui uTi = ...
i,2 i,1 u2i,2
..
.
..
.
ui,2 ui,n
..
.
(6.3.24)
|H
ui,n ui,1 ui,n ui,2 u2i,n
11 (ωi )| |H12 (ωi )| |H1n (ωi )|
|H
..
.
..
.
|Hn1 (ωi )| |Hn2 (ωi )| |Hnn (ωi )|
die bekanntlich den Rang Eins haben[46] und deren Zeilen oder Spalten alle von
einer linear abhängig sind. Es muß also bereits eine Zeile, z.B. die erste, alle
notwendigen Informationen liefern. Aus Gl. (6.3.24) erhält man zunächst:
ui,1 = 2ζi ω2i |H11 (ωi )|, i = 1, . . . , n, (6.3.25)
6.3.2.3 Beispiel
mit 1 0 2 −1
M= , K= , C = µM + κK
0 1 −1 2
mit µ = 1/10 und κ = 1/40.
Das folgende ist als Maple-Worksheet ausgeführt:
> restart;
> with(linalg):
Warning, new definition for norm
6.3 Experimentelle Modalanalyse 147
0
M :=
0 1
K:=matrix([[2,-1],[-1,2]]);
>
2 −1
K :=
−1 2
> mu:=1/10;
1
µ :=
10
> kappa:=1/40;
1
κ :=
40
> C:=mu*M+kappa*K;
1 1
C := M+ K
10 40
> A:=K-omega^2*M+I*omega*C;
1 1
A := K − ω2 M + I ω ( M+ K)
10 40
> H:=inverse(A);
40 − 20 ω + 3 I ω 2
40 + I ω
!
−16
H :=
−8 ""
%1
40 + I ω
%1
40 − 20 ω2 + 3 I ω #
−8 −16
%1 %1
%1 := −960 + 1287 ω − 176 I ω − 320 ω + 96 I ω3
2 4
> plot(abs(H[1,1]),omega=0..3):
> plot(argument(H[1,1]),omega=0..3):
> plot(abs(H[2,1]),omega=0..3):
> plot(argument(H[2,1]),omega=0..3):
> eiv:=eigenvectors(K);
> v1:=eiv[1,3,1];
148 Technische Anwendungen
4.5
Q
4
3.5
3
−1/2
Q2=2 Q
2.5
|H11(ω)|
1.5
0.5
ω ω
11 12
0
0 0.5 1 1.5 2 2.5 3
ω
v1 := [1, 1]
> v1n:=v1/norm(v1,2);
1
v1n := v1 2
2
> v2:=eiv[2,3,1];
v2 := [−1, 1]
> v2n:=v2/norm(v2,2);
1
v2n := v2 2
2
> V:=augment(v1n,v2n);
6.3 Experimentelle Modalanalyse 149
ω1 ω2
0
−0.5
−1
−1.5
∠ H11(ω)
−2
−2.5
−3
−3.5
0 0.5 1 1.5 2 2.5 3
ω
!
""
1 1
2 −
V :=
2
2
1
2
1 #
2 2
2 2
VT:=transpose(V);
!
>
""
1 1
VT :=
2 2
2
1
2
1 #
− 2 2
2 2
Delta:=multiply(VT,C,V);
!
>
""
1
∆ :=
0
8
7 #
0
40
> Lambda:=multiply(VT,K,V);
150 Technische Anwendungen
1 0
Λ :=
0 3
> zeta[1]:=Delta[1,1]/sqrt(Lambda[1,1])/2;
1
ζ1 :=
16
> zeta[2]:=Delta[2,2]/sqrt(Lambda[2,2])/2;
7
ζ2 := 3
240
Betrag und Argument von H11 und H21 sind in den Abbildungen 6.7, 6.8, 6.9
und 6.10 dargestellt und werden in Matlab nachbearbeitet.
6.3.2.3.1 Eigenfrequenz und Dämpfung Aus den Abb. 6.7 und 6.8 werden
folgende Werte entnommen und berechnet:
3.5
2.5
|H21(ω)|
1.5
0.5
0
0 0.5 1 1.5 2 2.5 3
ω
wobei die Vorzeichen nach Gl. (6.3.27) gefunden werden. Das Ergebnis ist in gu-
ter Übereinstimmung mit den exakten theoretischen Werten des Maple-Worksheet.
1 global F
x=0:0.001:3;
h11=-16*(40-20*x.^2+3*i*x)./(-960+1287*x.^2-176*i*x-320*x.^4+96*i*x.^3);
h21=-8*(40+i*x)./(-960+1287*x.^2-176*i*x-320*x.^4+96*i*x.^3);
5 %Eigenwert 1, Daempfung 1
x1=fmin(’h’,sqrt(1)-1/6,sqrt(1)+1/6,[1,1e-8,1e-8])
F=-h(x1)/sqrt(2)
x11=fzero(’h11’,x1-1/10)
x21=fzero(’h11’,x1+1/10)
10 zeta1=(x21-x11)/2/x1
y=0:0.01:x21;
z=F*ones(length(y),1);
y0=0:0.01:x1;
z0=-h(x1)*ones(length(y0),1);
15 y1=0:0.01:-h(x11);
z1=x11*ones(length(y1),1);
152 Technische Anwendungen
1
∠ H21(ω)
−1
−2
−3
−4
0 0.5 1 1.5 2 2.5 3
ω
y2=0:0.01:-h(x21);
z2=x21*ones(length(y2),1);
%Eigenwert 2, Daempfung 2
20 x2=fmin(’h’,sqrt(3)-1/6,sqrt(3)+1/6,[1,1e-8,1e-8])
F=-h(x2)/sqrt(2)
x12=fzero(’h11’,x2-1/10)
x22=fzero(’h11’,x2+1/10)
zeta2=(x22-x12)/2/x2
25 plot(x,abs(h11));grid;
xlabel(’\omega’);
ylabel(’|H_{11}(\omega)|’);
hold on
plot(y,z,’--’);plot(y0,z0,’--’);
30 plot(z1,y1,’--’);plot(z2,y2,’--’);
%gtext(’Q’);gtext(’Q_2=2^{-1/2}Q’);
%gtext(’\omega_{11}’);gtext(’\omega_{12}’);
hold off
pause
35 y1=0:-0.01:h12(x1);
y2=0:-0.01:h12(x2);
z1=x1*ones(length(y1),1);
6.3 Experimentelle Modalanalyse 153
z2=x2*ones(length(y2),1);
plot(x,angle(h11));grid;
40 xlabel(’\omega’);
ylabel(’\angle H_{11}(\omega)’);
hold on
plot(z1,y1,’--’);
plot(z2,y2,’--’);
45 %gtext(’\omega_1’);gtext(’\omega_2’);
hold off
pause
plot(x,abs(h21));grid;
xlabel(’\omega’);
50 ylabel(’|H_{21}(\omega)|’);
pause
plot(x,angle(h21));grid;
xlabel(’\omega’);
ylabel(’\angle H_{21}(\omega)’);
55 pause
%Probe der Daempfungen
xg1=fmin(’g’,sqrt(1)-1/6,sqrt(1)+1/6,[1,1e-8,1e-8])
xg2=fmin(’g’,sqrt(3)-1/6,sqrt(3)+1/6,[1,1e-8,1e-8])
F=-g(xg1)/sqrt(2)
60 x11=fzero(’h22’,xg1-1/10)
x21=fzero(’h22’,xg1+1/10)
zeta1=(x21-x11)/2/xg1
F=-g(xg2)/sqrt(2)
x12=fzero(’h22’,xg2-1/10)
65 x22=fzero(’h22’,xg2+1/10)
zeta2=(x22-x12)/2/xg2
%Betraege der Eigenvektoren
a1a1=-2*zeta1*h(x1)*x1^2
a1a2=-2*zeta1*g(xg1)*xg1^2
70 u11=sqrt(a1a1)
u12=a1a2/u11
a11=-2*zeta2*h(x2)*x2^2
a12=-2*zeta2*g(xg2)*xg2^2
u21=sqrt(a11)
75 u22=a12/u21
function y=h(x)
y=-abs(-16*(40-20*x.^2+3*i*x)./(-960+1287*x.^2-176*i*x-320*x.^4+96*i*x.^3));
80 function y=h11(x)
global F
y=abs(-16*(40-20*x.^2+3*i*x)./(-960+1287*x.^2-176*i*x-320*x.^4+96*i*x.^3))-F;
154 Technische Anwendungen
function y=h12(x)
85 y=angle(-16*(40-20*x.^2+3*i*x)./(-960+1287*x.^2-176*i*x-320*x.^4+96*i*x.^3));
Wenn Sie die Student-Edition von Matlab besitzen, dann ist eine symbolisch-
numerische Lösung durch Zugriff auf den Maple-Kern optimal.
2
Anhang A
Integral-Transformationen
A.1 Fourier-Transformation
A.1.1 Fourier-Reihen
Diese einfachen Beziehungen erhält man durch Multiplikation von Gl. (A.1.1)
mit cos nωt bzw. sin nωt und Integration über eine Periode T .
2
1 Jean-Baptiste Joseph Fourier (1768–1830), französischer Mathematiker und Physiker, be-
deutende Beiträge zur Thermodynamik und zu den Differentialgleichungen der Physik.
156 Integral-Transformationen
Wenn die Funktion x(t) mit der Periode T zu N diskreten Zeitpunkten tk als xk
gemessen wurde, dann läßt sie sich als diskrete Fourier-Transformation darstel-
len:
a0
N/2 2πitk 2πitk
x(tk ) = xk = + ai cos + bi sin , k = 1, . . . , N. (A.1.3)
2 T T
i=1
1
N
a0 = xk , (A.1.4a)
N
k=1
1
N
2πik
ai = xk cos , (A.1.4b)
N N
k=1
1
N
2πik
bi = xk sin . (A.1.4c)
N N
k=1
A.1.3 Fenster-Funktionen
Wenn die Abtastung des Signals nicht bei ganzzahligen Vielfachen der Periode
T erfolgt, kommt es zu Verfälschungen der Fourier-Approximation, da dann
eine Funktion mit einer veränderten Periode zugrunde liegt. Um diesen Effekt
abzuschwächen, wird das Signal mit einer Gewichtungsfunktion bewertet, die die
Ränder des Abtastintervalls gering gewichtet. Man spricht von einer Bewertung
mit einer Fenster-Funktion.2 Eine häufig verwendet Fenster-Funktion ist die
Hanning-Funktion:
1 k
w(k) = 1 − cos 2π , k = 1, . . . , n. (A.1.5)
2 n+1
Alternative Darstellung: Fourier-Reihen können als Sinus- oder Cosinus-Terme alleine beschrieben wer-
den, z.B.
∞
x(t) = d0 + dn cos(nωt − φn ),
n=1
mit
d0
= a0 /2, (A.1.6a)
= a2n + b2n ,
dn
b n
(A.1.6b)
φn = arctan . (A.1.6c)
an
2 Englisch: Windowing.
A.1 Fourier-Transformation 157
Gerade Funktionen genügen der Bedingung x(−t) = x(t); für diese fallen die Gerade und ungerade
Sinus-Terme der Fourier-Entwicklung weg, und es bleibt Funktionen:
∞
1
x(t) = a0 + an cos nωt.
2
n=1
Für ungerade Funktionen gilt x(−t) = −x(t), und in diesem Fall gilt
∞
x(t) = (bn ) sin nωt.
n=1
Die harmonischen Funktionen an cos(nωt und bn sin(nωt sind die harmoni- Frequenz-Spektrum:
schen Anteile der Ordnung n der periodischen Funktion x(t). Die Koeffizienten
der Gln. (A.1.2) und (A.1.6a) können als vertikale Linien über den zu ihnen
gehörenden Frequenzen nω aufgetragen werden und bilden dann das (diskrete
)Frequenzspektrum der Funktion x(t). Auf diese Weise kann jede Funktion im
Zeitbereich als x(t) oder im Frequenzbereich als z.B an (nω), n = 1, . . . , ∞,
dargestellt werden. Beide Darstellungsweisen haben je von Fall zu Fall ihre
Vorteile[29].
Zur numerischen Berechnung der Koeffizienten in Gl. (A.1.2) wird die schnelle Numerische Berechnung:
Fourier-Transformation3 von Cooley und Tuckey[29] verwendet. Dazu wird das
Gleichungssystem (A.1.4) in Matrix-Vektor-Darstellung
Aa = x
umformuliert, mit x dem Vektor mit den Komponenten xk der abgetasteten dis-
kreten Funktionswerte, a dem Vektor der zu berechnenden Fourier-Koeffizienten
a0 , an , bn und A der Koeffizientenmatrix mit den Eingängen cos 2πitk /T und
sin 2πitk /T . Der Algorithmus von Cooley und Tuckey beinhaltet eine effiziente
Inversion der Koeffizientenmatrix zur Berechnung der Fourier-Koeffizienten
a = A−1 x.
A.1.5 Beispiel
die Dreieckswelle folgendermaßen approximieren:
8 2π 1 6π 1 10π
f(t) = cos t + cos t+ cos t + ... . (A.1.7)
π2 T 9 T 25 T
Übung 24 Überprüfen Sie die Approximationsgüte der Gl. (A.1.7) indem Sie
sukzessive die Summanden addieren. Benützen Sie Matlab oder Gnuplot oder
ein ähnliches Ihnen verfügbares Programm.
A.2 Laplace-Transformation
A.2.1 Definition
die eine im Zeitbereich gegebene Funktion f(t) in die Funktion F(s) im Bildbe-
reich transformiert, der durch die komplexe Variable s = σ + jω charakterisiert
ist.
Lineare, zeitinvariante Differentialgleichungen lassen sich durch Integral-Trans-
formationen in rationale Funktionen ‘algebraisieren’ und lösen, wenn als Kern
der Transformation
K(s; t) = e−st
gewählt und die Integrationsgrenzen t0 = 0 und t1 = ∞ gewählt werden:
Regel 2
∞
F(s) = [f(t)] = 0
f(t)e−st dt.
∞
−st
1
−st
T
[f(t)] = [1] = e dt = lim − e
0 T→∞
1s −sT
0
1
1
= lim − e + e0 =
T→∞ s s s
No F(s) f(t), t 0
∞
1 F(s) 0
e−st f(t)dt
d
2 sF(s) − f(0) f(t)
dt
d2
3 s2 F(s) − sf(0) − f ¼ (0) f(t)
dt2
dn
4 sn F(s) − sn−1 f(0) − fn−1 (0) f(t)
dtn
F(s) t
5 0
f(τ)dτ
s
t t
8 F(s) G(s) 0
f(τ)g(t − τ)dτ = 0
g(τ)f(t − τ)dτ
160 Integral-Transformationen
9 1 δ(t)
e−as
10 σ(t − a)
s
1
11 σ(t)
s
n!
12 tn
sn+1
1
13 e−at
s+a
1 1
14 tn−1 e−at
(s + a)n (n − 1)!
a
15 1 − e−at
s(s + a)
s
16 (1 − at)e−at
(s + a)2
1 1
17 (1 − e−at − ate−at )
s(s + a)2 a2
1 e−at − e−bt
18
(s + a)(s + b) b−a
s+c (c − a)e−at − (c − b)e−bt
19
(s + a)(s + b) b−a
1 1
b a
20 1− e−at + e−bt
s(s + a)(s + b) b−a b−a
ab
s+c 1 b(c − a) −at a(c − b) −bt
21 c− e + e
s(s + a)(s + b) ab b−a b−a
1 e−at e−bt e−ct
22 + +
(s + a)(s + b)(s + c) (b − a)(c − a) (c − b)(a − b) (a − c)(b − c)
s+d (d − a)e−at (d − b)e−bt (d − c)e−ct
23 + +
(s + a)(s + b)(s + c) (b − a)(c − a) (c − b)(a − b) (a − c)(b − c)
a
24 sin at
s2 + a2
s
25 cos at
s2 + a2
A.2 Laplace-Transformation 161
a
26 e−bt sin at
(s + b)2 + a2
s+b
27 e−bt cos at
(s + b)2 + a2
ω2n ωn
28 e−ζωn t sin ωn 1 − ζ2 t
s2 + 2ζωn s + ω2n 1 − ζ2
s 1
29 − e−ζωn t sin(ωn 1 − ζ2 t − θ), θ = arccos ζ
s2 + 2ζωn s + ω2n 1 − ζ2
s + 2ζωn 1
30 e−ζωn t sin(ωn 1 − ζ2 t + θ), θ = arccos ζ
s2 + 2ζωn s + ω2n 1 − ζ2
ω2n 1
31 1− e−ζωn t sin(ωn 1 − ζ2 t + θ), θ = arccos ζ
s(s2 + 2ζωn s + ω2n ) 1 − ζ2
s + 2ζωn
32 e−ζωn t sin(ωn 1 − ζ2 t + θ), θ = arccos ζ
s(s2 + 2ζωn s + ω2n )
ω 1
33 arctan sin ωt
s t
0 falls t < 0,
σ(t) =
1 falls t 0,
σ(t) − σ(t − T )
δ(t) = lim
T→0 T
definiert werden kann. Zusammen mit der Einheitsrampenfunktion ρ(t)
ρ(t) = tσ(t)
A.2.2 Eigenschaften
Die Laplace-Transformation ist eine lineare Transformation mit folgenden Ei- Linearität:
genschaften einer linearen Transformation:
[kf(t)] = k [f(t)]
162 Integral-Transformationen
e−at
[f(t − t0 )] = e F(s),
f(t) = F(s + a),
(A.2.5)
(A.2.6)
t t
−1
[F1 (s)F2 (s)] = f1 (t − τ)f2 (τ)dτ = f1 (τ)f2 (t − τ)dτ. (A.2.7)
0 0
Bemerkung: Die erste Eigenschaft (A.2.1) des Satzes bildet die Grundlage zur Algebrai-
sierung linearer Differentialgleichungen mit konstanten Koeffizienten. Die Ei-
genschaften (A.2.3) und (A.2.4) bilden die sogenannten Grenzwertsätze. Die
Eigenschaft (A.2.5) zeigt, daß auch Systeme mit Totzeit t0 mit der Laplace-
Transformation behandelt werden können. Die letzte Eigenschaft (A.2.7) wird
als Faltungseigenschaft bezeichnet.
Beweis: Hier sollen als Beispiele für die Anwendung der Regel 2 einige Aussagen des
Satzes 1 bewiesen werden.
a) Differentiation, Gl. (A.2.1): Als Werkzeug aus der Integralrechnung wird die
partielle Integration verwendet[6].
∞
¼
{f (t)} = e−st f ¼ (t)dt
0
∞
= e−st f(t)
∞
+s
e−st f(t)dt
0 0
= −f(0) + s{f(t)}
= −f(0) + sF(s).
A.2 Laplace-Transformation 163
{f ¼¼ (t)} = s2 {f(t)} − sf(0) − f ¼ (0) wird bewiesen durch Ersetzen von f ¼ (t)
durch f ¼¼ (t) in der o.a. Gleichungskette, der allgemeine Beweis von Gl. (A.2.1)
geschieht durch Induktion. (Versuchen Sie es als Wiederholung der Grundvor-
lesung Mathematik!)
b) Verschiebung der t-Achse, Gl. (A.2.5): Zunächst soll gezeigt werden, daß mit
der Heaviside-Funktion σ(τ) folgendes gilt:
0 f ür t < t0
f̃(t) =
f(t − t0 ) sonst
Dann gilt {f̃} = e−t0 s F(s), wobei F(s) = {f(t)} ist, und damit ist die Bezie-
hung (A.2.8) hergeleitet.
Nun ist
∞
e−t0 s F(s) = e−t0 s e−sτ f(τ)dxτ
∞ 0
−s(t0 +τ)
= e f(τ)dxτ;
τ=0
c) Verschiebung der s-Achse, Gl. (A.2.6): Es gilt nach Definition der Laplace-
Transformation ∞
F(s) e−st f(t)dt
0
d) Faltungsintegral, Gl. (A.2.7): Die Definition von F2 (s) und Gl. (A.2.5) führen
auf ∞ ∞
e−sτ F2 (s) = e−st f2 (t − τ)σ(t − τ)dt = e−st f2 (t − τ)dt.
0 τ
Dies und die Definition von F1 (s) ergeben:
∞ ∞ ∞
F1 (s)F2 (s) = e−sτ f1 (τ)F2 (s)dt = f1 (τ) e−st f2 (t − τ)dtdxτ
0 0 τ
∞ t
−st
= e f1 (τ)f2 (t − τ)dxτdt,
0 0
A.3 Partialbruchzerlegung
Rationale Funktionen P(s)/Q(s), mit P(s) und Q(s) als Polynomen beliebiger
Ordnung, werden nicht immer in Korrespondenzentabellen zu finden sein. Um
von ihnen die inverse Laplace-Transformation bilden zu können, erweist sich die
Partialbruchzerlegung als wichtiges Hilfsmittel. Es sei
P(s) bm s m + + b1 s + b0
F(s) = = ,
Q(s) a n sn + + a 1 s + a 0
wobei ohne Beschränkung der Allgemeinheit m < n gesetzt werden darf. Nach
dem Fundamentalsatz der Algebra kann Q(s) in seine Nullstellen si faktorisiert
werden:
n
Q(s) = (s − si ) = (s − s1 )(s − s2 )
(s − sn ).
i=1
k1 k2 kn
F(s) = + + +
s − s1 s − s2 s − sn
zerlegen, wobei das Vorgehen durch den Typ der Nullstellen bestimmt wird.
Die Koeffizienten kj für einfache reelle Nullstellen sj werden wie folgt berechnet !
kj = (s − sj )F(s)|s=sj , j = 1, . . . , n.
166 Integral-Transformationen
gewonnen.
ki (s − α) + ki+1 β
F(s) = +
(s − α)2 + β2
Diese komplexe Gleichung kann durch Vergleich der Realteile und der Imaginär-
teile nach den beiden Unbekannten ki und ki+1 aufgelöst werden. Die inverse
Laplace-Transformation des betrachteten ersten Terms von F(s) ist unter Beach-
tung der No. 16 und 17 in Tabelle A.1 besonders einfach zu berechnen.
Y(s) 2s + 3
T (s) = = 2 .
U(s) s + 2s + 1
Es soll die Antwort y(t) des Systems im Zeitbereich auf einen Einheitssprung
u(t) = σ(t) als Systemeingang berechnet und deren Graph skizziert werden.
Y(s) 2s + 3
T (s) = = 2 .
U(s) s + 2s + 1
5 Mehrfache komplexe Polpaare kommen in der Praxis äußerst selten vor, deshalb wird auf
deren Behandlung verzichtet.
A.3 Partialbruchzerlegung 167
[(s − s ) Y(s)]
(l−j)
1 d
kij = i
l
, j = 1, . . . , l,
(l − j)! ds(l−j) s=si
d 2s + 3
[(s − s ) Y(s)]
d(2−1)
1
k21 = 2
2
= = −3,
(2 − 1)! ds(2−1)
ds s
2s + 3
s=s2 s=−1
k22 =
1 d(2−2)
(2 − 2)! ds(2−2)
[(s − s2 )2 Y(s)] =
s = −1.
s=s2 s=−1
Damit ergibt die Partialbruch-Zerlegung
2s + 3 3 3 1
Y(s) = = − − ,
s(s + 1)2 s s + 1 (s + 1)2
und unter Verwendung der Korrespondenzen-Tabelle wird die Sprungantwort im
Zeitbereich
y(t) = 3 − 3e−t − te−t .
2
6 In der Theorie komplexer Funktionen werden die kij als Residuen bezeichnet.
168 Integral-Transformationen
2
Übung 27 Wiederholen Sie Aufgabe 26 für folgende Funktion
5
Y(s) = .
s(s + 1)3 (s + 2)
dn y dn−1 y dy
an n
+ an−1 n−1 + + a1 + a0 y
dt dt dt
dm u dm−1 u du
= bm m
+ b m−1 m−1
+ + b1 + b0 u, m < n,
dt dt dt
deren sämtliche Anfangsbedingungen,
dn y(0)
= 0, ... y(0) = 0,
dtn
verschwinden mögen, was wegen der Linearität (Superponierbarkeit) der be-
trachteten dynamischen Systeme keine Einschränkung bedeutet.
! Anwendung der ersten Eigenschaft des Satzes 1, [di y/dti ] = si Y(s), i =
0, . . . , n, ergibt als Laplace-Transformation der Impulsantwort einer Differen-
tialgleichung deren Übertragungsfunktion
oder
(ms2 + ds + c)Y(s) = m[y(0)s + ẏ(0)] + dy(0) + U(s)
oder, mit allen Anfangsbedingungen zu Null gesetzt,
Y(s) 1
= .
U(s) ms + ds + c
2
2
Beispiel 33 Ein pneumatisches Drossel-Speicher-System hat die Übertragung-
sfunktion
Xa (s) 1
= .
Xe (s) Ts + 1
Erweiterung mit dem Nenner der rechten Seite ergibt
mit xe (t) dem zeitlichen Verlauf des Eingangsdruckes vor der Drossel und xa (t)
dem zeitlichen Verlauf des Ausgangsdruckes im Kessel. Die Zeitkonstante T =
RC ist das Produkt aus Drosselwiderstand R und Kesselkapazität C.
2
Übung 28 Betrachten Sie die Differentialgleichung
d2 x(t) dx(t)
2
+2 + 3x(t) = 4σ(t),
dt dt
mit σ(t) der Einheitssprungfunktion.
Zunächst seien die Anfangswerte x(0) = 0 und x ¼ (0) = 0. Lösen Sie die Diffe-
rentialgleichung durch Transformation in den Bildbereich, Partialbruchzerlegung
und Rücktransformation in den Zeitbereich.
Nehmen Sie nun x(0) = −1 und x ¼ (0) = 3 an und wiederholen Sie die Aufgabe.
Interpretieren Sie die unterschiedlichen Ergebnisse.
2
170 Integral-Transformationen
Anhang B
Matrixalgebra
B.1 Definitionen
Eine m
n Matrix ist eine rechteckige Anordnungen von (ganzen, reellen, kom- Matrizen:
plexen, etc.) Zahlen oder Funktionen:
a11 a12 a1,n−1 a1,n
a21 a22 a2,n−1 a2,n
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
A= .. .. .. .. .. .. .. ..
. . . . . . . .
a
.. .. .. .. .. .. .. ..
. . . . . . . .
m−1,1 am−1,2 am−1,n−1 am−1,n
am,1 am,2 am,n−1 am,n
mit m Zeilen und n Spalten. Die einzelnen Zahlen sind die Komponenten, Ele-
mente oder Eingänge der Matrix, das Element in der i-ten Zeile und j-ten Spalte
wird mit aij bezeichnet. Gilt m = n, dann heißt die Matrix quadratisch. Das
Tupel (m, n) wird Dimension der Matrix genannt.
Matrizen mit einer Spalte oder einer Zeile heißen Vektoren, genauer Spalten- Vektoren:
bzw. Zeilenvektoren. Im allgemeinen werden Vektoren als Spaltenvektoren an-
gegeben, Zeilenvektoren sind transponierte Spaltenvektoren: a ist ein Spalten-
vektor, aT ist derselbe Vektor zum Zeilenvektor transponiert. Vektoren sind in !
aller Regel als Spaltenvektoren notiert; alle Ausnahmen davon sind ausdrück-
lich erwähnt. Die ganze Zahl n der Elemente des Vektors heißt Dimension des
Vektors.
Matrizen werden hier und im folgenden in Großbuchstaben und fett geschrieben, Konventionen:
Vektoren in Kleinbuchstaben und ebenfalls fett. Komponenten oder Elemente
von Vektoren werden normal gedruckt und mit einem Index versehen, z.B. xi ,
Elemente von Matrizen werden ebenfalls normal gedruckt und mit zwei Indizes
172 Matrixalgebra
gekennzeichnet, z.B. aij , wobei der erste Index der Zeilenindex, der zweite der
Spaltenindex ist.
Spezielle Matrizen: Quadratische Matrizen, deren nicht-diagonale Elemente Null sind, heißen Dia-
gonalmatrizen, Diagonalmatrizen, deren Diagonalelemente alle Eins sind, heißen
Einheitsmatrizen. Die Einheitsmatrix wird I geschrieben. Will man ihre Dimen-
sion n noch angeben, dann indiziert man die Matrix In .
Eine quadratische Matrix, die durch Spiegelung der oberen rechten Hälfte an
der Diagonalen entsteht, heißt symmetrische Matrix.
Beispiel 34 Folgende Matrix A erfüllt die Definition und ist damit symmetrisch:
2
A= 1
1 3
4 0
3 0 7
Transponierung: Die Transponierte einer Matrix wird durch Vertauschen von Zeilen und Spal-
ten der Matrix erzeugt. Sie wird durch ein hochgestelltes T nach der Matrix
gekennzeichnet. Ein transponierter Spaltenvektor ist ein Zeilenvektor.
Beispiel 35
3 4 −2
8
A=
−6 17 , =
3 8 9 0
.
9 2 1 AT 4 −6
−2 17
2 4
1 7
0 4 7
2
Addition: Matrizen oder Vektoren gleicher Dimension werden elementweise addiert oder
subtrahiert.
Multiplikation Matrizen oder Vektoren werden elementweise mit einer skalaren Größe multipli-
mit einem Skalar: ziert oder dividiert.
Skalarprodukt: Das Skalarprodukt eines Zeilenvektors mit einem Spaltenvektor gleicher Dimen-
sion ist die Summe der elementweisen Produkte der einzelnen Vektorelemente,
d.h. explizit
2
Das Produkt einer m
n Matrix A mit einer n
p Matrix B ist eine m
p Matrix Matrixprodukt:
C, deren Elemente cij die Skalarprodukte der Zeilenvektoren ai , i = 1, . . . , m
der Matrix A mit den Spaltenvektoren bj , j = 1, . . . , p der Matrix B sind:
n
cij = ai bj = aik bkj , i = 1, . . . , m, j = 1, . . . , p.
k=1
Die Anzahl der Spalten von A muß gleich der Anzahl der Zeilen von B sein,
sonst ist die Matrixmultiplikation nicht definiert.
Beispiel 37 für ein Matrix
Matrix-Produkt:
1 7
4 = 1854
8 9
2 3 24 30
5 6
4 5 6 69 84
1 2 3
2
Es gilt folgendes:
(AB)T = BT AT .
Aber es gilt im allgemeinen nicht: !
AB = BA,
Die Determinante |A| einer Matrix wird auch als det(A) bezeichnet. Die zum Determinante:
Element aij definierte Unterdeterminante Aij ist die Determinante der Matrix,
die aus der Matrix A entsteht, wenn deren i-te Zeile und j-te Spalte gestrichen
werden. Diese wird auch Minore genannt.
Beispiel 38 Die Unterdeterminante A13 der Matrix
2
A= 1
1 1
1 −1 (B.3.1)
2 1 3
ist 1
= = −1.
1
A13
2 1
2
174 Matrixalgebra
2
Expansion der Die Ordnung einer Determinante kann um eine Ordnung reduziert werden, in-
Determinante: dem man sie um irgendeine Zeile oder Spalte in die Kofaktoren expandiert.
2
Beispiel 40 Folgende Determinante D werde um die zweite Spalte expandiert:
4 1 5
D = 2 2 1
0
4
3
2 2
= 1(−1)1+2
1
3
2 + 2(−1)
2+2 2 5
3
+ 0(−1) 3+2 4 5
1
= −10 − 8 = −18.
2
Adjungierte Matrix: Die zur Matrix A adjungierte Matrix adj(A) ist die Transponierte der Matrix
der Kofaktoren:
T
c11 c12 c13
adj(A) = c21 c22 c23 .
c31 c32 c33
Beispiel 41 Die Matrix A sei wieder wie im Beispiel (B.3.1), dann ist deren
Adjungierte
4 −5 −1
T
adj(A) = −2 4
0 .
−2 3 1
1 Sie sind sich bewußt, daß solch ein Überprüfen kein Beweis der Eigenschaften ist!
B.3 Iterative Operationen 175
Die Inverse A−1 einer quadratischen Matrix ist, falls sie existiert, die Lösung Inverse:
der Matrixgleichung
AX = I.
adj(A)
A−1 = . (B.3.2)
|A|
Beispiel 42 Die Matrix A sei wieder wie im Beispiel (B.3.1), dann ist deren
Inverse
2 −1 −1
A−1
= − 5
2 2 3
2
.
− 12 0 1
2
Welche Bedingung muß an die Koeffizienten aij gestellt werden, damit A inver-
tierbar ist?
AT A = AAT = I, (B.3.4)
und damit gilt auch folgende einfache Relation zur Ermittlung der Inversen einer
orthogonalen Matrix:
A−1 = AT . (B.3.5)
Die Eigenwerte λ einer Matrix sind die Lösungen der Polynomgleichung Eigenwerte:
|λI − A| = 0. (B.3.6)
176 Matrixalgebra
0 λ 0 − 1 1 −1
= 0 0 λ 3
λ − 2 −1 −1
2 1
−1 λ − 1 1
= −2 −1 λ − 3
= λ3 − 6λ2 + 9λ − 2 = 0,
2
Es gelten die Beziehungen
n
n
|A| = λi und spur(A) = λi ,
i=1 i=1
wobei spur(A), die Spur von A, die Summe der Diagonalelemente der Matrix A
ist.
Übung 30 Verifizieren Sie diese Beziehungen für das gerade behandelte Beispiel.
2
Übung 31 Bilden Sie die Eigenwerte der Matrix aus Übung 29.
2
Eigenvektoren: Ein Eigenvektor x zum Eigenwert λ einer Matrix A ist durch folgende homogene
Gleichung definiert:
(λI − A)x = 0. (B.3.7)
Eigenwerte und Eigenvektoren sind unerläßliche Hilfsmittel bei der Behand-
lung von Schwingungsproblemen mit mehreren Freiheitsgraden, s.a. Kap. 3 und
Kap. 4.
Rang einer Matrix: Die m-dimensionalen Vektoren x1 , x2 , . . . , xn heißen linear abhängig, wenn es
Zahlen c1 , c2 , . . . , cn gibt, die nicht alle Null sind und die Gleichung
c1 x1 + c2 x2 + + cn xn = 0
1
Beispiel 45 Die 3
2-Matrix
A= 2
2
4
3 7
hat den Rang min(3, 2) = 2, die 2
3-Matrix
1 2 3
B=
2 4 6
dagegen nur den Rang 1 = min(3, 2).
2
Quadratische Matrizen A, mit dim(A) = n, haben genau dann vollen Rang n,
wenn ihre Determinante nicht verschwindet
|A| = 0.
1
Übung 33 Hat die Matrix
A= 4
2 3
5 6
7 8 9
vollen Rang?
2
Eine quadratische Matrix, deren Determinante Null ist, heißt singuläre Matrix. Singuläre Matrix:
Wichtige Literatur zur Matrix-Algebra bilden die Monographien von Golub und
Van Loan[21] und Stewart[44]. In der Maschinendynamik sind die Matrizen i.A.
symmetrisch; dadurch spezialisiert sich das Eigenproblem auf das symmetrische.
Das Standardwerk dazu hat Parlett[38] verfaßt.
178 Matrixalgebra
Anhang C
Linearisierung
C.1 Taylor-Entwicklung
Für die Untersuchung mit unseren Methoden muß diese Gleichung linearisiert
werden.
Dazu wird eine nominale Lösung x0 (t) von (C.1.1) betrachtet, die durch den
Eingang u0 (t) und eine gewisse Anfangsbedingung x0 (t) erzeugt worden ist.
Um diese nominale Lösung werde die i-te Komponente der Differentialgleichung
(C.1.1) in eine Taylor-Reihe bis zum ersten Glied entwickelt:
ẋi (t) = fi (x0 , u0 ) +
n
∂fi (x, u)
∂x
(xj − x0j )
j=1
j x0 ,u0
∂f (x, u)
+
m
∂u
i (uj − u0j ),
j=1
j x0 ,u0
i = 1, . . . , n.
Nun seien
∆xi = xi − x0i
und
∆ui = ui − u0i ,
1 EinSystem wird nicht-autonom genannt, wenn in den rechten Seiten der Differentialglei-
chungen die Zeit t explizit auftritt: ẋ(t) = f(x(t),u(t), t).
180 Linearisierung
und daher
∆ẋi = ẋi − ẋ0i .
Es gilt
ẋ0i = fi (x0 , u0 ),
und damit kann Gleichung (C.1.2) folgendermaßen umgeformt werden:
∆ẋi (t) =
n
∂fi (x, u)
∂xj
∆xj (t) +
m
∂fi (x, u)
∂uj
∆uj (t),
j=1 x0 ,u0 j=1 x0 ,u0
C.2 Standard-Form
ẋ = Ax + Bu. (C.2.1)
∂f1 ∂f1
∂f1
∂u1 ∂u2 ∂um
∂f2 ∂f2 ∂f2
∂u1 ∂u2 ∂um
B = .. .. .. .. ..
(C.2b)
. . . . .
.. .. .. .. ..
. . . . .
∂fn ∂fn
∂fn
∂u1 ∂u2 ∂um
C.3 Beispiele und Übungen 181
ẍ + ẋ + ẋx = u,
lautet.
Die Jacobi-Matrizen werden als
∂f1 ∂f1
A=
∂x1 ∂x2 0 1
=
∂f2 ∂f2 −x2 −x1 − 1 x0 ,u0
∂x1 ∂x2
x0 ,u0
∂f1 0
B= ∂u =
∂f2 1
∂u x0 ,u0
berechnet, und mit x0 = (x01 x02 )T als Nominallösung erhält man folgende li-
neare Zustandsgleichungen in den Variablen ∆x und ∆u
∆ẋ1 = ∆x2 ,
∆ẋ2 = −x02 ∆x1 − (x01 + 1)∆x2 + ∆u, (C.3.1)
∆ẋ1 = ∆x2 ,
∆ẋ2 = −∆x1 − 2∆x2 + ∆u.
2
Beachten Sie, daß die linearen Terme der nichtlinearen Darstellung unverändert
in die lineare Darstellung übergehen. Dies ist plausibel, da die Linearisierung
der linearen Darstellung an dieser nichts ändern soll.
Übung 34 Untersuchen Sie den Fehler in Beispiel 46, wenn die Linearisierung
für
0 x1 2, 0 x2 2,
gültig sein soll.
182 Linearisierung
Hinweis: Vergleichen Sie für extreme x die lineare mit der nichtlinearen Zus-
tandsdarstellung.
2
Beispiel 47 Betrachten sie folgendes lineare System
u = (1 − eK|x1 | )sign(x1 ),
+1 if z > 0
sign(z) = ,
−1 if z < 0
Für x01 = 0 ergibt sich a = K und für sehr großes x01 strebt a → 0. Im ersten
Fall ist ∆ẋ2 = K∆x1 , im zweiten bleibt ∆x2 praktisch konstant.
2
Übung 35 Gegeben sei das nichtlineare System
1
ẋ1 = − ,
x22
ẋ1 = ux1 .
a) Bilden sie die linearisierte Zustandsdarstellung zunächst allgemein für belie-
biges x01 , x02 und u0 .
b) Nun nehmen Sie u0 (t) = 0 und x01 (t0 ) = x02 (t0 ) = 1 an. Zeigen Sie, daß für
diesen Eingang und mit diesen Anfangsbedingungen die allgemeine Lösung der
Differentialgleichungen
x01 (t) = 1 − t,
x02 (t) = 1,
C.3 Beispiele und Übungen 183
ist.
c) Bilden Sie für diese x0 (t) und u0 (t) die Jacobi-Matrizen. Hier tritt die Zeit
explizit auf !
2
184 Linearisierung
Anhang D
Lösung gewöhnlicher
Differentialgleichungen
In wenigen Fällen gering komplexer Systeme niedriger Ordnung werden die Zus-
tandsgleichungen
ẋ(t) = Ax(t) + Bu(t), x(t0 ) = x0 , (D.1a)
y(t) = Cx(t), (D.1b)
analytisch und exakt gelöst werden. In aller Regel werden die Lösungen mit dem
Computer numerisch und approximativ gefunden.
Wir werden auch den ersten Weg nicht vernachlässigen, da er Einsicht in die
Lösungsstruktur vermittelt. Für den zweiten Weg werden wir auch Computer-
Prozeduren angeben.
Es sei eine Zeile i der Gleichung (D.1a) und eine Zeile j der Gleichung (D.1b)
herausgegriffen:
ẋi (t) = ai1 x1 (t) + + ain xn (t) + bi1 u1 (t) + + bim um (t),
xi (t0 ) = xi0 , (D.1a)
yj (t) = cj1 x1 (t) + + cjn xn (t), (D.1b)
und
Yj (s) = cj1 X1 (s) + + cjn Xn (s)
oder umgeformt
Mit der Beziehung (B.3.2) kann die Inverse (sI − A)−1 in Gleichung (D.2a) als
adj(sI − A)
(sI − A)−1 = (D.1.3)
|sI − A|
Von der Darstellung (4.2) der Lösung der Zustandsgleichung gelangt man, wenn
die Anfangsbedingungen x0 = 0 gesetzt werden, sofort zur Übertragungsfunk-
tion G(s):
Y(s)
G(s) = = C(sI − A)−1 B. (D.1.4)
U(s)
Beispiel 48 Das Matrix-Tripel (A, B, C) der Strecke des balancierten Stabes ist
0 1
0
A= , b= , cT = (1 0).
g/l 0 −1
−1
1 s −1 0 −1
c (sI − A)
T
b = 0
−g/ls s 1 −1
1 s − g/l s − g/l 0
0 g/l −1
2 2
=
s
s2 − g/l s2 − g/l
−1
= .
s2 − g/l
2
Übung 36 Wiederholen Sie den Rechengang aus Beispiel 48 für folgende Daten:
0 1
0
A= , b= , cT = (1 0).
−2 −3 1
Überprüfen Sie Ihr Ergebnis durch Vergleich mit den Matrix-Elementen der Nor-
malform.
2
Zunächst sei die Eingangsfunktion u(t) Null, d.h. wir suchen die Lösung xh (t)
des homogenenen Systems
Die Lösung sei genügend glatt, so daß sie als Reihenansatz geschrieben werden
kann
xh (t) = α0 + α1 (t − t0 ) + α2 (t − t0 )2 + . (D.1.6)
Mit t = t0 erhält man sofort x0 = α0 . Differentiation von (D.1.6) und Substi-
tution in (D.1.5) ergibt
und daraus α1 = Ax0 für t = t0 . Dieser Vorgang wird fortgesetzt und führt zu
1 1
xh (t) = I + A(t − t0 ) + A2 (t − t0 )2 + A3 (t − t0 )3 + x0 .
2 6
188 Lösung gewöhnlicher Differentialgleichungen
Der Ausdruck in der eckigen Klammer ist die Reihendarstellung der Matrix-
Exponentialfunktion
∞
(t − t0 )k
eA(t−t0 ) = Ak , (D.1.7)
k!
k=0
die unter gewissen Voraussetzungen sehr schnell konvergiert [34], so daß nur eine
geringe Anzahl kmax der Summanden in (D.1.7) berechnet werden muß. Damit
errechnet sich die homogene Lösung zu
kmax
(t − t0 )k
xh (t) = Ak x0 . (D.1.8)
k!
k=0
Hinweis: Lösen Sie die homogene Gleichung zuerst von t0 nach t1 und dann
von t1 nach t2 .
2
Die partikuläre Lösung xp (t) der Zustandsgleichungen wird durch Variation der
Parameter erhalten. Folgender Lösungsansatz sei angenommen
oder aufgelöst
ξ̇(t) = e−A(t−t0 ) Bu(t),
und integriert
t
ξ(t) = e−A(τ−t0 ) Bu(τ)dτ.
t0
und unter Verwendung von (D.1.9) erhalten wir eine partikuläre Lösung der
Zustandsgleichung
t
xp (t) = eA(t−τ) Bu(τ)dτ. (D.1.11)
t0
D.1 Analytische Lösung 189
D.1.2.3 Gesamtlösung
Die homogene Lösung (D.1.8) und die Partikularlösung (D.1.11) ergeben die
Gesamtlösung
Wie schon die Lösung (D.2a) im Frequenzbereich zeigt auch die Lösung (D.12a)
im Zeitbereich zwei Summanden: der erste beschreibt den Einfluß des Anfang-
swertes x(t0 ) auf die Lösung (Anteil der Lösung der homogenen Differentialglei-
chung), der zweite den Einfluß einer Eingangsfunktion u(t).
Beachten Sie, daß eA(t−t0 ) in (D.12a) eine Matrix ist mit elementweise Expo-
nentialfunktionen. Diese Matrix wird als Transitionsmatrix
des dynamischen Systems (D.1a) bezeichnet, und mit ihr wird (4.12) zu
t
x(t) = Φ(t, t0 )x0 + Φ(t, τ)Bu(τ)dτ, (D.14a)
t0
t
y(t) = CΦ(t, t0 )x0 + C Φ(t, τ)Bu(τ)dτ. (D.14b)
t0
Ein Vergleich der beiden Lösungen (D.2a) und (D.15a) zeigt folgenden Zusam-
menhang zwischen der Laplace-Transformierten der Transitionsmatrix Φ(t) und
der Matrix (sI − A)−1 :
−1
{Φ(t)} = (sI − A) . (D.1.16)
190 Lösung gewöhnlicher Differentialgleichungen
Übung 38 Rechnen Sie dieses Ergebnis für den skalaren Fall durch.
Die Methode der unbestimmten Koeffizienten ist ein einfacher Spezialfall der
Methode der Variation der Parameter[28] und ist u.a. auf die Ansätze der Erre-
gerfunktionen f(t) in Tab. D.1 anwendbar, bei denen die Ableitungen von f(t)
der Funktion f(t) ähnlich sind.
f(t) xp (t)
aeωt Aeωt
atn A0 + A1 t + + An−1 tn−1 + An tn
a cos ωt A0 cos ωt + A1 sin ωt
a sin ωt A0 cos ωt + A1 sin ωt
1. Ist f(t) eine der Funktionen in Spalte 1 von Tab. D.1, dann wird der
korrespondierende Ansatz xp (t) aus Spalte 2 gewählt. Die unbestimmten
Koeffizienten des Ansatzes werden als Lösung eines linearen Gleichungs-
systems erhalten, das sich durch Einsetzen von xp und seiner Ableitungen
in Gl. (D.1.17) ergibt.
2. Löst f(t) den homogenen Teil von Gl. (D.1.17), dann wird der entspre-
chende Ansatz xp (t) mit t multipliziert, wenn die Lösung einer einfachen
Wurzel der charakteristischen Gl. von (D.1.17) entspricht oder mit t2 im
Falle einer doppelten Nullstelle usf..
3. Ist f(t) Summe von Komponenten aus Spalte 1, dann ist xp (t) Summe der
korrespondierenden Komponenten aus Spalte 2.
xp = A2 t2 + A1 t + A0 ,
D.1 Analytische Lösung 191
ẋp = 2A2 t + A1 ,
ẍp = 2A2 .
Der Koeffizientenvergleich der Polynome auf der rechten und linken Seite ergibt
folgendes lineare Gleichungssystem:
A2 = 2,
4A2 + A1 = −1,
4A2 + 2A1 + A0 = 0,
2
Übung 39 Bestimmen Sie für dieses Beispiel die Gesamtlösung für x0 = 1 und
ẋ0 = −1.
2
Übung 40 Warum funktioniert xp = A2 t2 + A1 t nicht?
2
Beispiel 50 Regel 2. Lösen Sie
ẍ + 5ẋ + 6x = e−2t .
Der natürliche“ Ansatz xp = Aet führt nicht zum Ziel. (Probieren Sie ihn und
”
beurteilen Sie den Mißerfolg!)
Wir wählen
xp = Ate−2t
und
xp = te−2t
ergibt.
2
Übung 41 Bestimmen Sie für dieses Beispiel die Gesamtlösung für x0 = 0 und
ẋ0 = 0. Vorsicht: Es ist nicht etwa xh 0!
2
Ein Maple-Worksheet vollzieht Beispiel und Übung nach:
> restart;
> ode:=diff(x(t),t,t)+5*diff(x(t),t)+6*x(t)=exp(-2*t);
2
(−2 t)
ode := ( ∂t
∂
2 x(t)) + 5 ( ∂t x(t)) + 6 x(t) = e
∂
> dsolve({ode,x(0)=0,D(x)(0)=0},x(t));
x(t) = −e(−2 t) + t e(−2 t) + e(−3 t)
> f:=exp(-2*t);
f := e(−2 t)
> a:=A*t*f;
a := A t e(−2 t)
> b:=diff(a,t);
b := A e(−2 t) − 2 A t e(−2 t)
> c:=diff(a,t,t);
c := −4 A e(−2 t) + 4 A t e(−2 t)
> leq:=c+5*b+6*a=f;
leq := A e(−2 t) = e(−2 t)
> A:=solve(leq,A);
A := 1
Beispiel 51 Regel 3. Lösen Sie
ẍ + 5ẋ + 6x = e2t + t2 .
xp (t) = A3 e2t + A2 t2 + A1 t + A0
1 2t 1 2 5 19
= e + t − t+ ,
20 6 18 108
und sie wird in folgender Maple-Sitzung erarbeitet.
> restart;
> ode:=diff(x(t),t,t)+5*diff(x(t),t)+6*x(t)=exp(2*t)+t^2;
2
(2 t)
ode := ( ∂t
∂
2 x(t)) + 5 ( ∂t x(t)) + 6 x(t) = e
∂
+ t2
D.1 Analytische Lösung 193
> dsolve({ode,x(0)=0,D(x)(0)=0},x(t));
1 37 (−3 t) 1 (−2 t)
x(t) = (27 e(5 t) + 90 e(3 t) t2 − 150 t e(3 t) + 95 e(3 t) ) e(−3 t) + e − e
540 135 2
> dsolve({ode},x(t));
1
x(t) = (27 e(5 t) + 90 e(3 t) t2 − 150 t e(3 t) + 95 e(3 t) ) e(−3 t) + C1 e(−3 t) + C2 e(−2 t)
540
> expand(%);
1 t 2 1 2 5 19 C1 C2
x(t) = (e ) + t − t+ + +
20 6 18 108 (et )3 (et )2
> f:=exp(2*t)+t^2;
f := e(2 t) + t2
Koeffizientenvergleich
> eq[1]:=op(1,op(1,lhs(leq)))*op(2,op(1,lhs(leq)))=1;
eq 1 := 20 A3 = 1
> eq[2]:=coeff(lhsc,t,2)=1;
eq 2 := 6 A2 = 1
> eq[3]:=coeff(lhsc,t,1)=0;
eq 3 := 6 A1 + 10 A2 = 0
> eq[4]:=coeff(lhsc,t,0)=0;
eq 4 := 5 A1 + 6 A0 + 2 A2 = 0
19 −5 1 1
{A0 = , A1 = , A2 = , A3 = }
108 18 6 20
2
Aufgabe 31 Lösen Sie
D.2.1 Diskretisierung
D.2.1.1 Lösungsalgorithmus
konstant (Halteglied nullter Ordnung). Dadurch kann u() aus dem Integranden
genommen werden. Änderung der Variablen η = kT + T − τ erleichtert die
Integration:
T
x(kT + T ) = eAT x(kT ) + eAη dηBu(kT ). (D.2.2)
0
Mit
AT A2 T 2
Ψ(T ) = I + + +
2! 3!
oder numerisch stabiler
AT AT AT AT
Ψ(T ) = I + I+ I+
2 3 N−1 N
definiert wird. Dann berechnet sich
und
∞
∞
Ak T k+1 Ak T k
Γ (T ) = B= T B = Ψ(T )T B. (D.2.6)
(k + 1)! (k + 1)!
k=0 k=0
10 %
11 % x[n+1] = Phi * x[n] + Gamma * u[n]
12 %
13 % assuming a zero-order hold on the inputs
14 % and sample time T.
15
16 error(nargchk(3,3,nargin));
17 error(abcdchk(a,b));
18
19 [m,n] = size(a);
20 [m,nb] = size(b);
21 s = expm([[a b]*t; zeros(nb,n+nb)]);
22 Phi = s(1:n,1:n);
23 Gamma = s(1:n,n+1:n+nb);
24
25
26 function e = expm(a)
27 % Matrix exponential via Taylor series.
28 % Scale A by power of 2 so that its norm is < 1/2 .
29 s = round(log(norm(a,1))/log(2)+1.5);
30 if s < 0,
31 s = 0;
32 end
33 a = a/2^s;
34
35 % Taylor series for exp(A)
36 [m,n] = size(a);
37 k = 1;
38 e = 0*a;
39 f = eye(m,n);
40 while norm(e+f-e,1) > 0,
41 e = e + f;
42 f = a*f/k;
43 k = k+1;
44 end
45
46 % Undo scaling by repeated squaring
47 for k = 1:s,
48 e = e*e;
49 end
D.2.2 Taylor-Entwicklung
D.2.2.1 Lösungsalgorithmus
die Funktion f genügend oft nach x, u und t differenzierbar ist. Wenn ∂f/∂x =: fx
und ∂f/∂u =: fu stetig sind, dann besitzt die Differentialgleichung (D.2.7) eine
eindeutige Lösung x(t). Diese entwickeln wir in eine Taylor-Reihe um den Punkt
t = t0 :
(t − t0 )2
x(t) = x0 + (t − t0 )ẋ(t0 ) + ẍ(t0 ) + , (D.2.8)
2!
deren Ableitungen wir nicht explizit kennen, da die Lösung nicht bekannt ist.
Wir beschaffen die Ableitungen durch totale Differentiation von (D.2.7)
Übung 43 Bilden Sie die dritte totale Ableitung der Differentialgleichung (D.2.7).
2
Diese Übung zeigt, daß die höheren totalen Ableitungen sehr komplexe Aus-
drücke werden. Daraus werden alternativ folgende Konsequenzen gezogen:
1. die Taylor-Entwicklung wird nach wenigen Termen abgebrochen; dann ist
(D.2.8) nur in einer kleinen Umgebung von t0 gültig,
2. die totalen Ableitungen in der Taylor-Entwicklung werden durch zusätzliche
Auswertung der Funktion f an geeigneten Zwischenwerten im Integration-
sintervall gewonnen; dann muß ein höherer numerischer Aufwand einge-
setzt werden, (D.2.8) ist aber in einer größeren Umgebung von t0 gültig.
Es sei nun angenommen, daß die Entwicklung (D.2.8) eine brauchbare Näherung
der Lösung für einen Schritt h = t − t0 darstellt. Dann können x und alle seine
totalen Ableitungen an der Stelle t = t0 + h ausgewertet werden, und mit ihnen
kann ein neuer Schritt t + h gemacht werden. Mit der Definition des Operators
T
h ¼ hk−1 (k−1)
Tk (x, u, t) = f(x, u, t) + f (x, u, t) + + f (x, u, t), (D.2.10)
2! k!
wird dieser Vorgang als Algorithmus 3 formalisiert
198 Lösung gewöhnlicher Differentialgleichungen
2. setze
ti = t0 + ih, i = 0, . . . , N,
h(k+1) (k)
R = f (x(τ), u(τ), τ),
(k + 1)!
h(k+1) (k+1)
= x (τ), ti < τ < ti + h,
(k + 1)!
zeigt, daß der Fehler der Approximation von der Ordnung k + 1, (hk+1 ), ist.
Algorithmus 4 Euler-Verfahren:
Schritt 3:
xi+1 = xi + hf(xi , u(ti ), ti ), i = 0, . . . , N − 1.
Das Restglied,
R = 1/2h2 ẍ(τ),
Algorithmus 5 Heun-Verfahren:
Schritt 1 und 2 wie in Algorithmus 3,
Schritt 3:
k1 = hf(xi , u(ti ), ti ),
xi+1 = xi + 1/2[k1 + hf(xi + k1 , u(ti+1 ), ti+1 )], i = 0, . . . , N − 1.
2
Übung 44 Bilden Sie das Restglied und daraus die Ordnung des Heun-Verfahrens.
2
k1 = hf(xi , u(ti ), ti ),
k2 = hf(xi + 1/2k1 , u(ti + 1/2h), ti + 1/2h),
k3 = hf(xi + 1/2k2 , u(ti + 1/2h), ti + 1/2h),
k4 = hf(xi + k3 , u(ti + h), ti + h),
xi+1 = xi + 1/6(k1 + 2k2 + 2k3 + k4 ), i = 0, . . . , N − 1.
2
Der Fehler ist hier von (h5 ).
D.2.2.2 Computer-Prozedur
Als Beispiel für die Anwendung der im Abschnitt D.2 besprochenen numerischen
Lösungsverfahren wollen wir den mit einem PD-Regler und einem Zustandsre-
gler geregelten Stab simulieren. Aus Stabilitätsgründen wird nicht die ungere-
gelte Strecke berechnet, sondern die geregelte mit der Übertragungsfunktion
Y(s) −(KD s + KP )
G(s) = = 2 . (D.3.1)
U(s) s − KD s − (g/l + KP )
a0 = −(g/l + KP ), a1 = −KD ,
und
b0 = −KP , b1 = −KD .
D.3 Beispiel: Balancierter Stab 201
Damit läßt sich die Regler-Normalform des geregelten balancierten Stabes for-
ẋ (t) x (t) 0
mulieren:
1 0 1 1
= + u(t), (D.2a)
ẋ2 (t) g/l + KP KD
x (t)
1
x (t) 2 1
- b1
u - m ẋ2 - 1 x2 1 x1 ? y
+ r- r- b0 - +m
I
@
6 s s
@
−a1
−a0
KD = −2 und KP = −11.
Übung 45 Bitte bestätigen Sie durch Vergleich der Koeffizienten des charakte-
ristischen Polynoms der Übertragungsfunktion (D.3.1) die numerischen Werte
der Regler-Parameter.
2
Die numerischen Lösungsverfahren gelten für beliebiges (beschränktes) Eingang-
ssignal u(t). Für einen Vergleich beschränken wir uns hier auf den Einheitss-
prung, d.h.
0 für t < 0
u(t) = ,
1 für t 0
dessen Laplace-Transformierte U(s) = 1/s ist.
Die exakte Lösung des Problems (D.3.1) ist dann → Kompendium: Abschnitt
y(t) =
−KP
1 − e−ζωn t sin(βωn t + θ)
2.2.2
−(KP + g/l)
−KD
+ e−ζωn t sin βωn t, (D.3.3)
−(KP + g/l)
mit β = 1 − ζ2 und θ = arccos ζ.
202 Lösung gewöhnlicher Differentialgleichungen
Übung 46 Bitte rechnen Sie als Wiederholung des Stoffs der Vorlesung Rege-
lungstechnik diese Lösung nach. Überlegen Sie sich, daß deren zweiter Summand
der differenzierende Anteil des Zählers der Übertragungsfunktion ist.
D.3.1 Diskretisierung
Die Lösung des Problems (4.2) mit dem Diskretisierungsansatz (D.4a) ist in
Abbildung D.2 in der linken Hälfte dargestellt. Dabei ist das Zeitintervall t
[0, 5] in k = 128, k = 8 und k = 1 Teilintervalle aufgeteilt. Die Graphen der
entsprechenden Lösungen sind durch ausgezogene, gestrichelte, bzw. punktierte
Linienzüge gekennzeichnet.
Für enge Diskretisierung (k = 128) fällt der Graph mit der exakten Lösung
(4.3) zusammen. Die Lösungen der beiden anderen Diskretisierungen stimmen
nur für t = iT, i = 0, . . . , k, mit T dem Diskretisierungsintervall, mit der exakten
Lösung überein.
→ Rechenzeitvergleich Soll die Lösung nur für gewisse diskrete Zeitpunkte erhalten werden, dann ist
in der Vorlesungsübung für lineare Differentialgleichungen der Diskretisierungsansatz sehr ökonomisch
bezüglich der aufzuwendenden Rechenleistung. Die Lösungen für Zwischenzeit-
punkte können durch verschiedene Interpolationsansätze gewonnen werden: li-
neare Interpolation oder Spline-Interpolation.
! Zur Rekapitulation: Dieser Diskretisierungsansatz eignet sich nur für lineare
bzw. linearisierte Differentialgleichungen, deren Parameter nicht von der unabhän-
gigen Variablen t abhängen.
D.3.2 Reihen-Entwicklung
Die Lösung des Problems (4.2) mit Hilfe der Taylor-Entwicklung des Algorith-
mus 3 zur Ordnung k = 1 ist in Abbildung D.2 in der rechten Hälfte dargestellt.
Dabei sind die Schrittweiten h = 0.05, h = 0.2 und h = 1 Sekunden gewählt.
Die Graphen der entsprechenden Lösungen sind wieder durch ausgezogene, ges-
trichelte, bzw. punktierte Linienzüge gekennzeichnet.
Charakterisierung partieller
Differentialgleichungen
Consider the first order linear partial differential equation in two variables x and
y:
a(x, y)ux + b(x, y)uy + c(x, y)u = f(x, y). (E.0.1)
where v(ξ, η) = u(x(ξ, η), y(ξ, η)), and equation (E.0.2) is a linear first order
ordinary differential equation in ξ with η as a parameter.
We want the co-ordinate changes to be one-to-one, which requires a non-zero
ϕ
determinant of the Jacobian
J =
x ϕy
= ϕx ψy − ϕy ψx = 0. (E.0.3)
ψ x ψy
a(vξ ξx + vη ηx ) + b(vξ ξy + vη ηy ) + cv = f
or
(aξx + bξy )vξ + (aηx + bηy )vη + cv = f. (E.0.4)
206 Charakterisierung partieller Differentialgleichungen
ηx b
=−
ηy a
for ηy = 0.
dy ηx b
=− = . (E.0.5)
dx ηy a
η(x, y) = k (E.0.6)
defines a family of curves in the (x, y) plane which are called the characteristics
of problem (E.0.1).
ξ = x.
or in terms of (ξ, η)
C(ξ, η) F(ξ, η)
vξ + v= (E.0.7)
A(ξ, η) A(ξ, η)
Solution of ODE
d(ξ,η)dξ
Now, for the solution of equation (E.0.7), we multiply it by e to get
e d(ξ,η)dξ
vξ + d(ξ, η)e d(ξ,η)dξ
v = g(ξ, η)e d(ξ,η)dξ
.
Charakterisierung partieller Differentialgleichungen 207
As the left hand side is (e d(ξ,η)dξ
v)ξ the equation can be integrated with
respect to ξ with the result
e d(ξ,η)dξ
v= g(ξ, η)e d(ξ,η)dξ
dξ + h(η)
or
v(ξ, η) = e− d(ξ,η)dξ
g(ξ, η)e d(ξ,η)dξ
dξ + h(η)e− d(ξ,η)dξ
, (E.0.8)
where h(η) is any differentiable function of one variable. This expression has to
be transformed back into the original variables (x, y)
u(x, y) = eα(x,y) [β(x, y) + h(ψ(x, y))], (E.0.9)
which is the desired solution.
Example 1 With a = 0, b, and c as real numbers find the solution of
aux + buy + cu = 0.
The characteristic equation is
dy b
= ,
dx a
which can be integrated as
bx − ay = k.
Now take
ξ=x and η = bx − ay,
to get
avξ + cv = 0
or
c
v = 0.
vξ +
a
By multiplying this equation by e (c/a)dξ = e(c/a)ξ we get
c
e(c/a)ξ vξ + ve(c/a)ξ = 0
a
or
∂ (c/a)ξ
(e v) = 0.
∂ξ
Integration yields
e(c/a)ξ v = h(η)
or
v(ξ, η) = e−(c/a)ξ h(η)
and back transformation
u(x, y) = e−(c/a)x h(bx − ay).
[14] D. J. Ewins. Modal Testing: Theory and Practice. Research Studies Press,
New York, 1984.
210 LITERATURVERZEICHNIS
[27] D. Kraft. Tomp — fortran modules for optimal control calculations. ACM
Trans. Math. Softw., 20(3):262–281, 1994.
[31] N. M. M. Maia and D. J. Ewins. A new approach for the modal identifi-
cation of lightly damped structures. Mech. Systems and Signal Processing,
3(2):173–193, 1989.
[33] C. Moler and et al. Using Matlab Version 5. The MathWorks, Natick, MA,
1996.
[34] C. B. Moler and C. F. van Loan. Nineteen dubious ways to compute the
exponential of a matrix. SIAM Review, 20:801–836, 1978.
[35] M. B. Monagan, K. O. Geddes, G. Labahn, and S. Vorkoetter. Maple V
Programming Guide. Springer, New York, 1996.
[36] K. Ogata. Modern Control Engineering. Prentice Hall, Englewood Cliffs,
1990.
[37] P. V. O’Neil. Beginning Partial Differential Equations. John Wiley, New
York, 1999.
[38] B. N. Parlett. The Symmetric Eigenvalue Problem. Prentice-Hall, Engle-
wood Cliffs, 1980.
[39] M. J. D. Powell. Approximation Theory and Methods. Cambridge Univer-
sity Press, Cambridge, 1981.
[40] S. S. Rao. Mechanical Vibrations. Addison-Wesley, Reading, 1995.
[41] D. Redfern. The Maple Handbook. Springer, New York, 1996.
[42] L. F. Shampine and M. W. Reichelt. The matlab ode suite. SIAM J. Sci.
Comput., 18(1):1–22, 1997.
[43] B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, Y. Ikebe, V. C.
Klema, and C. B. Moler. Matrix Eigensystem Routines - EISPACK Guide.
Springer, New York, 1976.
[44] G. W. Stewart. Introduction to Matrix Computations. Academic Press,
New York, 1973.
[45] J. Stoer and R. Bulirsch. Introduction to Numerical Analysis. Springer,
New York, 1980.
[46] G. Strang. Linear Algebra and its Applications. Academic Press, New York,
1976.
[47] TheMathWorks. Matlab User’s Guide. TheMathWorks, Natick, 1992.
[48] W. T. Thomson. Theory of Vibration with Applications. Prentice-Hall,
Englewood Cliffs, 1981.
[49] J. Wittenburg. Schwingungslehre. Springer, Berlin, 1996.
[50] S. Wolfram. Mathematica – A System for Doing Mathematics by Computer.
Addison-Wesley, Redwood City, CA, 1991.
212 LITERATURVERZEICHNIS