Verifizieren der Lösung eines Anfangswertproblems

Aufgabe

Skizze zur Aufgabe

Zwei Massen m1 und m2 können auf einer vertikalen bzw. einer horizontalen Führung reibungsfrei gleiten. Sie sind durch eine starre Stange (Masse mS und Massenträgheitsmoment JS bezüglich des Schwerpunktes S) gekoppelt. Das gesamte System wird aus der skizzierten Lage ohne Anfangsgeschwindigkeit freigelassen.

Gegeben:

Gegebene Größen

Das nur auf numerischem Wege zu gewinnende Bewegungsgesetz soll punktuell mit dem Energiesatz kontrolliert werden.

Bewegungs-Differenzialgleichung, Energiesatz


Koordinaten

Die Aufgabe wird hier für verschiedene Fragestellungen mit Matlab gelöst. Zur Analyse des Bewegungsvorgangs muss zunächst das Anfangswertproblem gelöst werden (die Koordinate φ verfolgt die Bewegung der Verbindungsstange, siehe nebenstehende Skizze, die Herleitung der Differenzialgleichung wird im Kapitel "Kinetik starrer Körper" ausführlich behandelt):

Nichtlineares Anfangswertproblem

Wegen der vielfältigen Fehlerquellen wird im Kapitel "Kinetik starrer Körper" empfohlen, die Ergebnisse punktuell mit dem Energiesatz zu überprüfen.

Bewegungszustände für Energiebilanz

Die nebenstehende Skizze zeigt zwei Lagen des Systems, für die die Energiebilanz besonders einfach zu formulieren ist. Als Zustand 1 wird der Start der Bewegung und als Zustand 2 die horizontale Lage der Stange betrachtet. Mit dem Null-Potenzial in Höhe der horizontalen Führung hat das System nur im Zustand 1 potenzielle Energie. Im Zustand 2 befindet sich Masse 2 im Umkehrpunkt ihrer Bewegung, hat also keine kinetische Energie. Da sich damit auch der Bolzen 2 der Stange in Ruhe befindet, ist er der Momentanpol für die Stange, deren kinetische Energie damit besonders einfach formuliert werden kann. Dies wird im Kapitel "Kinetik starrer Körper" ausführlich beschrieben und führt über die Energiebilanz zu einer sehr einfachen Formel für die Winkelgeschwindigkeit der Verbindungsstange beim Passieren der Lage φ = 0:

Winkelgeschwindigkeit beim Nulldurchgang

Die numerische Integration des Anfangswertproblems (hier mit erweiterter Aufgabenstellung ausführlich demonstriert) liefert das Bewegungsgesetz in Abhängigkeit von der Zeit als Wertetabelle. Die benötigte Nullstelle von φ(t) könnte man daraus durch Interpolation annähern, um danach aus der ebenfalls als Wertetabelle vorliegenden Funktion ω(t) auch durch Interpolation den speziellen ω-Wert (auch näherungsweise) zu gewinnen.

Für diese Aufgabe bietet Matlab allerdings mit der Event-Strategie eine effektive und natürlich viel bessere Lösungsmöglichkeit. Man kann den Algorithmen der numerischen Integration "Events" vorgeben, deren Eintreten bei der Integration überprüft wird, hier also zum Beispiel: Nulldurchgang der Funktion φ(t), und es werden der genaue Zeitpunkt des "Events" und die zugehörigen Funktionswerte φ und ω berechnet und zusätzlich abgeliefert. Dies wird nachfolgend gezeigt.

Lösung mit Matlab, Kontrolle mit dem Energiesatz

Das Anfangswertproblem 2. Ordnung muss für die Lösung mit Matlab in ein System 1. Ordnung umgeschrieben werden. Dafür wird eine neue Variable ω als Ableitung der Winkelkoordinate φ nach der Zeit eingeführt:

Anfangswertproblem 1. Ordnung

In dem folgenden Matlab-Script ist das Differenzialgleichungsssystem in der Function Dgl realisiert, die dem Verfahren zur numerischen Integration (ode45) beim Funktionsaufruf angezeigt wird. Bemerkungen zur Event-Strategie findet man nach dem Listing des Scripts.

function BolzenkraefteEv

global JSS          m1S         m2S       gS 
       JSS = 1/12 ; m1S = 0.5 ; m2S = 2 ; gS = 29.43 ;

maxstep = 0.01 ;                                                 % ... fuer MaxStep-Option
phi0 = [asin(1/3) ; 0] ; tspan = [0 2] ;                         % Anfangswerte, Zeitintervall
options = odeset ('MaxStep' , maxstep , 'Events' , @Ereignis) ;  % Event-Behandlung wird angekuendigt 

[t phiphip te phiphipe ereig] = ode45 (@Dgl , tspan , phi0 , options) ;
phi = phiphip(:,1) ; phip = phiphip(:,2) ;                       % Sortieren der Ergebnisse

subplot(1,2,1) ; plot (t , phi)  , grid on , title  ('\phi(t)')
subplot(1,2,2) ; plot (t , phip) , grid on , title  ('\omega(t)')
 
disp([' ']) ;
disp(['Rechnung mit MaxStep = ' , num2str(maxstep) , ':']) ;
for i=1:length(te)                                                         % Auswertung der "Events"
     disp(['omega(t = ' , num2str(te(i)) , ') = ' , num2str(phiphipe(i,2))]) ;
end

disp(['phi  (t_end = ' , num2str(t(end)) , ') = ' , num2str(phi(end))])  ; % Zusaetzliche Kontrolle:
disp(['omega(t_end = ' , num2str(t(end)) , ') = ' , num2str(phip(end))]) ; % Werte am Ende des Intervalls

% ===================================================================
function [value,isterminal,direction] = Ereignis (t,phiphip) 

phi  = phiphip(1) ;
value(1)      = phi ; % Ereignis 1: "Nulldurchgang"
isterminal(1) = 0 ;   % Kein Abbruch, aber registrieren ...
direction(1)  = 0 ;   % ... in jedem Fall (Hin- und Rueckweg) 

% ============================================================
% Function, die die "rechte Seite" des Dgl.-Systems definiert: 
function f = Dgl (t , phiphip)

global JSS          m1S         m2S       gS 

phi   = phiphip(1) ;
omega = phiphip(2) ;

phiP  = omega ;
zaehl = -((m2S-m1S)*omega^2*sin(2*phi)+(2*m1S+1)*gS*cos(phi)) ;
nenn  =  (0.5+2*m2S*(sin(phi))^2+2*m1S*(cos(phi))^2+2*JSS) ;

f = [phiP ; zaehl/nenn] ;

Die Passagen im vorstehenden Script, die sich auf die Matlab-Event-Strategie beziehen, sind farblich gekennzeichnet (hellblau die Ankündigung und Definition der Events, weiß die Auswertung der abgelieferten Informationen):

Rechnungen mit mehreren MaxStep-Options

Die Rechnung wurde mehrmals mit unterschiedlichen MaxStep-Options durchgeführt. Die nebenstehend zu sehende Ausgabe in das Command Window zeigt, dass im Integrationsbereich drei Ereignisse registriert werden. Neben den "Ereignis-Winkelgeschwindigkeiten" werden auch die Endwerte für φ und ω ausgegeben. Man erkennt, dass die MaxStep-Option 0,01 offensichtlich für ausreichende Genauigkeit sorgt. Sowohl die Winkelgeschwindigkeiten als auch die Endwerte ändern sich bei feinerer MaxStep-Option nicht mehr.

Die theoretischen Vorgaben werden bestätigt: Die Winkelgeschwindigkeiten beim Passieren der Lage φ = 0 entsprechen genau dem nach dem Energiesatz zu erwartenden Wert (die Vorzeichen entsprechen der mit der Definition der Koordinate φ festgelegten positiven Drehrichtung). Die nachfolgend zu sehenden kinematischen Diagramme zeigen, dass die Ereignisse tatsächlich den Nullstellen der Funktion φ(t) zuzuordnen sind, und die ausgegebenen Winkelgeschwindigkeiten entsprechen ebenfalls den Werten, die das Diagramm anzeigt:

Kinematische Diagramme

Plausibilität der Kinematik, Animation


Die oben zu sehenden Diagramme φ(t) und ω(t) sind schwierig zu interpretieren. Immerhin zeigt sich bereits der zu erwartende periodische Verlauf, die Erfüllung der Anfangsbedingungen und die erwartete Zuordnung spezieller ω-Werte zu den Nulldurchgängen der Funktion φ(t). Eine Vorstellung vom zu erwartenden Bewegungsverlauf hat man allerdings eher für die beiden Gleitsteine, deren Bewegungsgesetze deshalb nachfolgend noch berechnet werden.

Außerdem ist aus der Kinematik bekannt, dass sich der Mittelpunkt der Verbindungsstange eines Doppelschiebers auf einer Kreisbahn bewegt. Auch dies wird mit der nachfolgenden Testrechnung überprüft.

Die für die Darstellung benötigten Koordinaten sind über die hier angegebenen Zwangsbedingungen aus der Funktion φ(t) leicht zu berechnen. Sie lauten in dimensionsloser Form (bezogen auf die Länge l der Verbindungsstange):

Bewegungskoordinaten der Gleitsteine und des Stangenmittelpunkts

Das Matlab-Script enthält nur noch die wesentlichen Passagen des oben gezeigten Scripts, mit denen die graphische Darstellung der gewünschten Funktionen gelingt:

function BolzenkraefteEv2

global JSS          m1S         m2S       gS 
       JSS = 1/12 ; m1S = 0.5 ; m2S = 2 ; gS = 29.43 ;

phi0 = [asin(1/3) ; 0] ; tspan = [0 2] ;                         % Anfangswerte, Zeitintervall
options = odeset ('MaxStep' , 0.01) ;

[t phiphip] = ode45 (@Dgl , tspan , phi0 , options) ;
phi = phiphip(:,1) ; phip = phiphip(:,2) ;                       % Sortieren der Ergebnisse
xS  = 0.5*cos(phi) ; yS  = 0.5*sin(phi)  ;                       % Mittelpunkt der Verbindungsstange
x2S = cos(phi)     ; y1S = sin(phi)      ;                       % Gleitsteine

subplot(1,3,1) ; plot (t  , x2S) , grid on , title  ('x_2^{\ast}(t)')
subplot(1,3,2) ; plot (t  , y1S) , grid on , title  ('y_1^{\ast}(t)')
subplot(1,3,3) ; plot (xS , yS ) , grid on , title  ('y(x)') , axis equal

% ============================================================
% Function, die die "rechte Seite" des Dgl.-Systems definiert: 
function f = Dgl (t , phiphip)

global JSS          m1S         m2S       gS 

phi   = phiphip(1) ;
omega = phiphip(2) ;

phiP  = omega ;
zaehl = -((m2S-m1S)*omega^2*sin(2*phi)+(2*m1S+1)*gS*cos(phi)) ;
nenn  =  (0.5+2*m2S*(sin(phi))^2+2*m1S*(cos(phi))^2+2*JSS) ;

f = [phiP ; zaehl/nenn] ;

Die Graphik-Ausgabe des Scripts zeigt die drei Diagramme:

Bewegungsgesetze der Gleitsteine, Bahnkurve des Stangemittelpunkts
Animation der Bewegung

Interpretation der Kontrollrechnung:

Download


Die oben zu sehenden Matlab-Scripts stehen als BolzenkraefteEv.m und BolzenkraefteEv2.m zum Download zur Verfügung.