Zur
Startseite
Laufkatze mit pendelnder Last - Berechnung mit der Matlab-Event-Strategie
Aufgabe
Eine Laufkatze (Masse mK)
trägt eine Last (Masse einschließlich
Anhängevorrichtung: mL,
Massenträgheitsmoment bezüglich des
Schwerpunktes S:
JL). In der skizzierten Ruhelage beginnt für
eine kurze Zeit Δt die konstante
Antriebskraft F0 zu
wirken, die danach wieder abgeschaltet wird. Nach dem Zurücklegen
der Strecke a stößt die Laufkatze auf einen
elastischen Puffer
(Federzahl c).
Die Bewegung von Laufkatze und Last soll, beginnend aus der
Ruhelage, für die ersten 10 Sekunden analysiert werden.
Gegeben:
mK = 100 kg ; JL = 400 kgm2 ; lS = 4 m ; F0 = 2000 N ;
mL = 500 kg ; c = 200000 N/m ; Δt = 1 s ; a = 5 m .
Die Aufgabe wird in den Kapiteln "Prinzipien der Mechanik"
und "Verifizieren von Computerrechnungen" behandelt.
Bewegungs-Differenzialgleichungen, Vorbereitung der Matlab-Rechnung
Die Besonderheit dieser Aufgabe besteht in dem Eintreten
von "Ereignissen"
(Abschalten der Antriebskraft, Zu- und Abschalten einer Feder).
Man erfasst sie, indem
an die Stelle der Kraft F0
die zeitabhängige Kraft Ft
tritt und die Federkonstante
c durch ct
ersetzt wird:
Unter Verwendung der nebenstehend skizzierten Koordinaten
gelten folgende Bewegungs-Differenzialgleichungen (die ausführlich
kommentierte Herleitung findet man im Kapitel
"Prinzipien der Mechanik"):
Die Matlab-Solver für die numerische Integration erwarten ein
Differenzialgleichungssystem 1. Ordnung. Deshalb werden die beiden
neuen Variablen v und
ω eingeführt:
Damit kann das Anfangswertproblem so formuliert werden, wie
man es Matlab anbieten kann:
Matlab-Script
Es wird der Matlab-Standard-Solver ode45 verwendet, dem die Informationen
übergeben werden müssen, die das Anfangswertproblem definieren. Diese
werden u. a. in drei Functions zusammengestellt:
- Die "rechte Seite" des Differenzialgleichungssystems (in der oben
angegebenen Zusammenstellung ist das der Vektor
r) wird in einer Function
RechteSeite (Name ist frei wählbar)
definiert. Der Name dieser Function wird als erster Parameter
des ode45-Aufrufs bekanntgegeben.
- Weil das Bewegungs-Differenzialgleichungssystem in den Beschleunigungsgliedern
gekoppelt ist, entsteht nach dem Umschreiben auf ein
System 1. Ordnung auf der linken Seite des Differenzialgleichungssystems
eine "Massenmatrix" M. Matlab kann
damit umgehen, man muss nur über die "options", die dem
Solver übergeben werden, ankündigen, dass
eine gesonderte Function diese Matrix beschreibt (hier wird
für diese Function - willkürlich - der
Name Massenmatrix gewählt, in den
"options" wird dieser Name nach dem Keyword "Mass" bekanntgemacht).
- Die Werte für die Antriebskraft Ft
und die Federzahl ct
müssen beim Eintreten der entsprechenden Ereignisse ("Abschalten
des Antriebs nach Δt" bzw. "Anschlag
an die Feder" und "Lösen von der Feder") geändert werden. Dafür wird hier
die "Matlab-Event-Strategie" genutzt (eine wesentlich einfachere, aber
"nicht ganz saubere" Lösung wird
hier gezeigt). Die
Ereignisse müssen definiert werden, wofür hier die Function
Ereignis (Name ist auch frei
wählbar) geschrieben wurde. Der Name wird in den
"options" nach dem Keyword "Events" bekanntgemacht.
Zusätzlich zu diesen Informationen werden der ode45-Function
noch die Anfangsbedingungen und der Zeitbereich mitgeteilt, über den integriert
werden soll. Der über die "options" transportierte
"MaxStep-Parameter" dient der Steuerung der Schrittweite für die
numerische Integration. Die Strategie der Berechnung und der
Verifizierung der Ergebnisse wird nach dem Listing des
Matlab-Scripts erläutert:
function LaufkatzeEvKontr ()
% Parameter (global definiert, damit Wertaenderung nur an einer Stelle erforderlich):
global mK mL JL ct lS Dt Ft a g
mK = 100 ; mL = 500 ; JL = 400 ; c = 200000 ; lS = 4 ; Dt = 1 ; F0 = 2000 ; a = 5 ; g = 9.81 ;
x0 = [0 ; 0 ; 0 ; 0] ; ct = 0 ; Ft = F0 ; % Anfangswerte
tend = 10 ;
tspan = [0 tend] ; % Zeitintervall
% Integration des Anfangswertproblems:
options = odeset ('Mass' , @Massenmatrix , 'MaxStep' , 0.001 , 'Events' , @Ereignis) ;
tout = [] ; xvout = [] ; % Auf tout und xvout werden die Ergebnisse gesammelt
cont = 1 ;
while cont == 1
[t xv te xve ereig] = ode45 (@RechteSeite , tspan , x0 , options) ;
tout = [tout ; t] ; % Die gerade abgelieferten Ergebnisse werden an ...
xvout = [xvout ; xv] ; % ... die bereits gespeicherten Ergebnisse angehaengt
lent = length (t) ;
if (t(lent) >= tend)
cont = 0 ; % Ende erreicht, Abbruch der while-Schleife
else % ... war ein "Ereignis" der Grund fuer den Abbruch
telen = length (te) ;
disp(['Ereignis ' , num2str(ereig(telen)) , ' nach t = ' , ...
num2str(te(telen)) , ' s bei x = ' , num2str(xve(telen)) , ' m']) ;
if ereig(telen) == 1 % Ereignis "Eine Sekunde nach dem Start"
Ft = 0 ; % Kraft abschalten
elseif ereig(telen) == 2 ct = c ; % ... Feder zuschalten
elseif ereig(telen) == 3 ct = 0 ; % ... Feder abschalten
end
tspan = [te(telen) tend] ; % ... weiter mit "Restintervall ...
x0 = [xve(telen,1) xve(telen,2) xve(telen,3) xve(telen,4)] ;
end % ... und neuen Anfangsbedingungen
end
% Grafische Ausgabe von x und phi, zunaechst "Sortieren der Ergebnisse":
x = xvout(:,1) ; % xvout ist die Matrix der Ergebnisse,
phi = xvout(:,2) ; % x, phi ...
v = xvout(:,3) ; % v und omega ...
omega = xvout(:,4) ; % ... sind Vektoren
Tges = 0.5*mK*v.^2 + 0.5*mL*(v.^2+lS^2*omega.^2+2*lS*v.*omega.*cos(phi)) ...
+0.5*JL*omega.^2-mL*g*lS*cos(phi) ; % Kontrollfunktion "Gesamtenergie"
subplot(3,1,1) ; plot (tout , x) , grid on , title ('Bewegung der Laufkatze: x(t)') , axis([0 5 -1 6])
subplot(3,1,2) ; plot (tout , phi) , grid on , title ('Pendeln der Last: \phi(t)')
subplot(3,1,3) ; plot (tout , Tges) , grid on , title ('Gesamtenergie: T_{ges}(t)')
disp(['Ergebnisse nach ' , num2str(length(tout)-1) , ' Zeitschritten:']) ; % "End"-Kontrolle
disp(['x (t_end = ' , num2str(tout(end)) , ') = ' , num2str(x (end))]) ;
disp(['phi(t_end = ' , num2str(tout(end)) , ') = ' , num2str(phi(end))]) ;
disp(['T_ges(t_end = ' , num2str(tout(end)) , ') = ' , num2str(Tges(end))]) ;
disp(['(zum Vergleich: T_ges(t = ' , num2str(0) , ') = ' , num2str(Tges(1)) , ')']) ;
% ===================================================================
% Function, die die Ereignisse definiert:
function [value,isterminal,direction] = Ereignis (t,xv)
global mK mL JL ct lS Dt Ft a g
x = xv(1) ; % Aktuelle Position der Laufkatze
value(1) = t-Dt ; % Ereignis 1: "Eine Sekunde nach dem Start"
isterminal(1) = 1 ; % Abbruch, ...
direction(1) = 0 ; % ... in jedem Fall
value(2) = x-a ; % Ereignis 2: "Punkt x=a erreicht"
isterminal(2) = 1 ; % Abbruch, ...
direction(2) = 1 ; % ... wenn Bewegung nach rechts
value(3) = x-a ; % Ereignis 3: "Punkt x=a erreicht"
isterminal(3) = 1 ; % Abbruch, ...
direction(3) = -1 ; % ... wenn Bewegung nach links
% ============================================================
% Function, die die "Massenmatrix" des Dgl.-Systems definiert:
function M = Massenmatrix (t , xv)
global mK mL JL ct lS Dt Ft a g
phi = xv(2) ;
M = [1 0 0 0 ;
0 1 0 0 ;
0 0 mK+mL mL*lS*cos(phi) ;
0 0 mL*lS*cos(phi) mL*lS^2+JL ] ;
% ===================================================================
% Function, die die "rechte Seite" des Dgl.-Systems definiert:
function r = RechteSeite (t , xv)
global mK mL JL ct lS Dt Ft a g
x = xv(1) ; phi = xv(2) ; v = xv(3) ; omega = xv(4) ;
b1 = mL*lS*omega^2*sin(phi)-ct*(x-a)+Ft ;
b2 = -mL*g*lS*sin(phi) ;
r = [v ; omega ; b1 ; b2] ;
Ereignisse (Events), Ankündigung und Definition
Die Option, die den Lösungsalgorithmus veranlasst,
auf Events zu achten, ist im Matlab-Script hellblau gekennzeichnet:
options = odeset ('Mass' , @Massenmatrix , 'MaxStep' , 0.005 , 'Events' , @Ereignis) ;
Damit wird auf eine Funktion Ereignis verwiesen,
die die Ereignisse definieren muss. In diesem Fall werden die drei
Ereignisse "Zeitpunkt Δt erreicht",
"Punkt x = a
erreicht, weiter in Richtung
x > a",
"Punkt x = a
erreicht, weiter in Richtung
x < a" definiert:
% Function, die die Ereignisse definiert:
function [value,isterminal,direction] = Ereignis (t,xv)
global mK mL JL ct lS Dt Ft a g
x = xv(1) ; % Aktuelle Position der Laufkatze
value(1) = t-Dt ; % Ereignis 1: "Eine Sekunde nach dem Start"
isterminal(1) = 1 ; % Abbruch, ...
direction(1) = 0 ; % ... in jedem Fall
value(2) = x-a ; % Ereignis 2: "Punkt x=a erreicht"
isterminal(2) = 1 ; % Abbruch, ...
direction(2) = 1 ; % ... wenn Bewegung nach rechts
value(3) = x-a ; % Ereignis 3: "Punkt x=a erreicht"
isterminal(3) = 1 ; % Abbruch, ...
direction(3) = -1 ; % ... wenn Bewegung nach links
Die Function Ereignis bekommt als Input die gleichen Werte
(Zeit t und
Vektor der Funktionswerte mit
x(t),
φ(t),
v(t),
ω(t))
wie die Functions für die Definition der Differenzialgleichungen
(hier: Massenmatrix und
RechteSeite).
Sie muss für jedes Ereignis 3 Werte abliefern.
Eine "Nullfunktion" (hier auf Vektor value) definiert
durch ihre Nullstelle das Ereignis, außerdem:
Was soll beim Eintreten dieses Ereignisses geschehen
(hier auf isterminal: 1 bedeutet Abbruch, 0 würde "Weiterrechnen,
Ereignis nur registrieren"
bedeuten), und soll dies in jedem Fall (die 0 auf direction)
geschehen bzw. nur beim Nulldurchgang von negativen in positive Werte (1)
oder nur beim Nulldurchgang von positiven in negative Werte (-1)?
Die Ereignisse werden mit ihrer Position in den drei
Vektoren identifiziert. Alle Ereignisse sollen hier zum Abbruch
führen. Die Ereignisse 2 und 3 unterscheiden sich durch
die Richtung des Nulldurchgangs der "Ereignisfunktion"
x−a.
Es wird also jeweils nur eines dieser beiden Ereignisse
registriert.
Reaktion auf das Eintreten eines Ereignisses, Strategie der Berechnung
Wenn die "Events-Option" gesetzt ist, liefert die ode45-Function
(neben dem Vektor der Zeitschritte und der Matrix der
berechneten Funktionswerte) drei zusätzliche Ergebnisse ab. Alle zur Auswertung
der Ereignisse gehörenden Passagen sind im Matlab-Script
weiß gekennzeichnet:
[t xv te xve ereig] = ode45 (@RechteSeite , tspan , x0 , options) ;
Auf te werden die Zeiten abgeliefert, zu denen die
Ereignisse eingetreten sind, auf xve findet man die
Zustandsgrößen der Bewegung zu diesen Zeitpunkten
(hier:
x(t),
φ(t),
v(t),
ω(t)),
und der Vektor ereig enthält die Ereignisnummern
(hier 1, 2 oder 3). Nach jedem Abbruch werden diese
Informationen ausgewertet, um die Parameter
(hier: Ft
und ct)
zu aktualisieren, außerdem wird das Ereignis im Command Window registriert:
telen = length (te) ;
disp(['Ereignis ' , num2str(ereig(telen)) , ' nach t = ' , ...
num2str(te(telen)) , ' s bei x = ' , num2str(xve(telen)) , ' m']) ;
if ereig(telen) == 1 % Ereignis "Eine Sekunde nach dem Start"
Ft = 0 ; % Kraft abschalten
elseif ereig(telen) == 2 ct = c ; % ... Feder zuschalten
elseif ereig(telen) == 3 ct = 0 ; % ... Feder abschalten
end
Weil die Rechnung mehrfach (bei jedem eigetretenen Ereignis)
abgebrochen wird, ist die ode45-Funktion in eine while-Schleife
eingebettet. Die Ergebnisse, die jeder ode45-Aufruf auf dem
Vektor t und der Matrix xv abliefert, werden deshalb
auf tout und xvout kumuliert (grüne Zeilen):
tout = [] ; xvout = [] ; % Auf tout und xvout werden die Ergebnisse gesammelt
cont = 1 ;
while cont == 1
[t xv te xve ereig] = ode45 (@RechteSeite , tspan , x0 , options) ;
tout = [tout ; t] ; % Die gerade abgelieferten Ergebnisse werden an ...
xvout = [xvout ; xv] ; % ... die bereits gespeicherten Ergebnisse angehaengt
lent = length (t) ;
if (t(lent) >= tend)
cont = 0 ; % Ende erreicht, Abbruch der while-Schleife
else % ... war ein "Ereignis" der Grund fuer den Abbruch
telen = ...
tspan = [te(telen) tend] ; % ... weiter mit "Restintervall ...
x0 = [xve(telen,1) xve(telen,2) xve(telen,3) xve(telen,4)] ;
end % ... und neuen Anfangsbedingungen
end
Die rot gekennzeichneten Zeilen realisieren die Strategie:
Bei jedem Abbruch wird gefragt, ob das Ende des Integrationsintervalls
erreicht ist. Die Rechnung wird dann entweder beendet oder fortgesetzt,
wobei im letztgenannten Fall das Integrationsintervall auf den
"Rest der noch nicht abgearbeiteten Zeit" gesetzt wird, und für
die Anfangsbedingungen der Folgerechnung werden die Zustandsgrößen
zum Zeitpunkt des Ereignisses eingesetzt.
Strategie zum Verifizieren der Ergebnisse
An zwei Stellen (gekennzeichnet durch graue Schriftfarbe)
sind Kontrollmöglichkeiten im Matlab-Script vorgesehen:
- Es bietet sich an, die Energiekurve als Kontrollfunktion parallel
zur numerischen Integration der Bewegungs-Differenzialgleichungen zu
ermitteln.
Die potenzielle Energie soll nur für das bewegte System berücksichtigt
werden (während des Kontakts mit der Feder ist auch in dieser potenzielle
Energie gespeichert).
Im Kapitel "Verifizieren von Computerrechnungen"
wird gezeigt, dass mit dem Null-Potenzial entsprechend
nebenstehender Skizze und den Anteilen der kinetischen
Energie aus der Bewegung von Laufkatze und pendelnder Last
folgende Beziehung für die Gesamtenergie des bewegten Systems gilt:
Diese Funktion wird als
Tges(t)
mit den ermittelten Werten für
x(t),
φ(t),
v(t) und
ω(t)
berechnet und graphisch dargestellt. Außerdem wird
der Zahlenwert für
Tges(tend)
in das Command Window ausgegeben.
- Die typische "Endkontrolle" zur Einschätzung, ob die
Schrittweitensteuerung ausreichend feine Zeitschritte erzeugte,
wird mit der Ausgabe der Werte
x(tend) und
φ(tend)
in das Command Window durchgeführt. Der "MaxStep"-Parameter
wird von Rechnung zu Rechnung verkleinert, bis schließlich
die Kontrollwerte ausreichend genau übereinstimmen.
Tges = 0.5*mK*v.^2 + 0.5*mL*(v.^2+lS^2*omega.^2+2*lS*v.*omega.*cos(phi)) ...
+0.5*JL*omega.^2-mL*g*lS*cos(phi) ; % Kontrollfunktion "Gesamtenergie"
subplot(3,1,1) ; plot (tout , x) , grid on , title ('Bewegung der Laufkatze: x(t)') , axis([0 5 -1 6])
subplot(3,1,2) ; plot (tout , phi) , grid on , title ('Pendeln der Last: \phi(t)')
subplot(3,1,3) ; plot (tout , Tges) , grid on , title ('Gesamtenergie: T_{ges}(t)')
disp(['Ergebnisse nach ' , num2str(length(tout)-1) , ' Zeitschritten:']) ; % "End"-Kontrolle
disp(['x (t_end = ' , num2str(tout(end)) , ') = ' , num2str(x (end))]) ;
disp(['phi(t_end = ' , num2str(tout(end)) , ') = ' , num2str(phi(end))]) ;
disp(['T_ges(t_end = ' , num2str(tout(end)) , ') = ' , num2str(Tges(end))]) ;
disp(['(zum Vergleich: T_ges(t = ' , num2str(0) , ') = ' , num2str(Tges(1)) , ')']) ;
Ergebnisse der Berechnung
Das Matlab-Script gibt die Ergebnisse als graphische Darstellungen
und zusätzlich ausgewählte Zahlenwerte in das Command Window aus.
Bei Problemen dieser Art ist auf eine geeignete Schrittweitensteuerung
zu achten, obwohl der Matlab-ode45-Solver mit automatischer
Schrittweitensteuerung arbeitet.
Die Schrittweiten werden über die "MaxStep-Option" beeinflusst.
Man sollte sie eigentlich immer verwenden und normalerweise
den Wert nicht größer als 0.1 setzen. Hier wurden vier Rechnungen
mit den "MaxStep-Options" 0.1, 0.01, 0.001 und 0.0001 durchgeführt.
Das ist sicherlich übertrieben, zeigt aber sehr schön, wie sich
die Ergebnisse stabilisieren. Die völlige Übereinstimmung
der Ergebnisse (in allen ausgegebenen Dezimalstellen)
der beiden letzten Rechnungen mit mehr als
40000 bzw. 400000 Integrationsschritten zeigt zusätzlich,
dass sich die Rundungsfehler bei der gewaltigen Anzahl
von Rechenoperationen noch nicht zu einem erkennbaren Fehler
kumulieren.
Als Indikatoren können bei jedem beliebigen Anfangswertproblem
die berechneten Funktionswerte am Ende des Integrationsintervalls
dienen, hier also
x(tend) und
φ(tend).
Bei dieser Rechnung können zusätzlich die Werte
zu den Ereigniszeitpunkten und der Wert für
Tges(tend)
verglichen werden.
Alle Indizien sprechen für eine "gesunde" Rechnung, und man
erkennt auch, dass mit der Matlab-Event-Strategie die
Zeitpunkte und die Funktionswerte der Ereignisse exakt
ermittelt werden, so dass auch die Schaltpunkte für die
Änderung der Parameter genau den Vorgaben der Aufgabenstellung
entsprechen.
Im Kapitel "Verifizieren von Computerrechnungen" wird
zusätzlich noch eine Betrachtung der Zahlenwerte für die
Gesamtenergie des bewegten Systems Tges
empfohlen: Weil das System aus der Ruhelage startet, hat es
zum Startzeitpunkt nur potenzielle Energie. Mit der oben
skizzierten Lage des Null-Potenzials
berechnet man:
Tges(t=0) = −mL g lS = −19620 Nm .
Dass dieser Wert exakt ausgegeben wird, ist selbstverständlich
(bis zu diesem Zeitpunkt ist ja noch nichts geschehen).
Dem System wird nur durch die Antriebskraft Energie zugeführt, und
weil diese Kraft konstant ist, kann die Arbeit, die von ihr
geleistet wird, aus "Kraft · Weg" berechnet werden.
Der bis zum Ereignis 1 (Abschalten der Antriebskraft) zurückgelegte
Weg wird mit
x(t=Δt) = 3.6642 m ausgewiesen,
so das sich aus "Anfangsenergie + Arbeit infolge Antriebskraft"
die Gesamtenergie
Tges(t=tend) = Tges(t=0) + F0 · x(t=Δt)
ergibt, die am Ende des Integrationsintervalls noch vorhanden
sein muss. Genau diesen Wert
Tges(t=tend) = −19620 Nm + 2000 N · 3.6642 m = −12291.6 Nm
findet man im Command Window.
Das Matlab-Script liefert als graphische Darstellungen die beiden
Bewegungsgesetze
x(t) (Laufkatze) und
φ(t) (Pendeln der Last)
und die Kontrollfunktion
Tges(t).
In der nebenstehenden Graphik sieht man, dass die Laufkatze
offensichtlich zweimal an die Feder stößt, bevor sie ihren
Rückweg antritt und mehrmals (auch ohne äußere Krafteinwirkung)
ihre Bewegungsrichtung ändert. Dieses etwas überraschende Verhalten
wird hier ausführlich diskutiert.
Die Kurve der Kontrollfunktion
Tges(t)
liefert ein weiteres Indiz für die Richtigkeit der Rechnung:
Startend mit der (negativen) potenziellen Energie,
wächst sie infolge der Antriebskraft auf ihren Maximalwert,
den sie zum Zeitpunkt
t = Δt
erreicht.
Die beiden "Zacken" in der Kurve kennzeichnen die Energieabgabe
an die Feder nach dem Anschlag. Die Energie wird komplett an
das bewegte System zurückgegeben, so dass der Rest des
Kurvenverlaufs wieder eine horizontale Gerade ist.
Die nachfolgende Abbildung links zeigt einen Zoom in
das Bewegungsgesetz der Laufkatze, der den besonders
interessanten Teil der Bewegung verdeutlicht.
Was im Diagramm noch etwas überraschend erscheint, wird mit der
unten rechts zu sehenden Animation der Bewegung verständlich:
Man sieht, wie die schwingende Last auf die Laufkatze einwirkt
und Änderungen der Bewegungsrichtung auch ohne äußere
Krafteinwirkung erzwingt:
|
|
Die Animation zeigt die Bewegung von Laufkatze und Last
etwa für den Zeitbereich, der im links zu sehenden Diagramm
als Bewegungsgesetz der Laufkatze dargestellt ist.
|
Download
Das oben zu sehende Matlab-Script steht als
LaufkatzeEvKontr.m
zum Download zur Verfügung.