Laufkatze mit pendelnder Last - Berechnung mit der Matlab-Event-Strategie
Matlab (einfach)  |   Matlab-Event-Strategie  |   Simulink  |   Runge-Kutta  |   Animation

Aufgabe


Laufkatze mit pendelnder Last

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

Generalisierte Koordinaten

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:

Zeitabhängige Antriebskraft, zeitabhängige Federzahl

Unter Verwendung der nebenstehend skizzierten Koordinaten gelten folgende Bewegungs-Differenzialgleichungen (die ausführlich kommentierte Herleitung findet man im Kapitel "Prinzipien der Mechanik"):

Bewegungs-Differenzialgleichungen

Die Matlab-Solver für die numerische Integration erwarten ein Differenzialgleichungssystem 1. Ordnung. Deshalb werden die beiden neuen Variablen v und ω eingeführt:

Definition von zusätzlichen Variablen

Damit kann das Anfangswertproblem so formuliert werden, wie man es Matlab anbieten kann:

Komplettes Anfangswertproblem

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:

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" xa. 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

Null-Potenzial

An zwei Stellen (gekennzeichnet durch graue Schriftfarbe) sind Kontrollmöglichkeiten im Matlab-Script vorgesehen:

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 von 4 Rechnungen

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) = −mg 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.

Graphikausgabe des Matlab-Scripts

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:

 
Bewegungsgesetz der Laufkatze Animation der Bewegung

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.