Aufgabe

Ein Doppelpendel wird definiert durch die beiden Pendelmassen m₁ und m₂, die auf die jeweiligen Schwerpunkte bezogenen Massenträgheitsmomente J_S1 und J_S2, die Schwerpunktabstände von den Drehpunkten s₁ und s₂ und den Abstand l₁ der beiden Drehpunkte voneinander.

Die Bewegung soll durch die Funktionen φ₁(t) und φ₂(t) beschrieben werden, die für das Zeitintervall t = 0 ... 10 s zu berechnen sind.

Die Aufgabe wird in den Kapiteln "Prinzipien der Mechanik" und "Verifizieren von Computerrechnungen" behandelt.

Die gegebenen Werte gelten für zwei schlanke Stäbe gleicher Masse und gleicher Länge. Sie sollen aus der nebenstehend skizzierten Anfangslage ohne Anfangsgeschwindigkeiten freigelassen werden, so dass folgende Anfangsbedingungen gelten:

Bewegungs-Differenzialgleichungen, Vorbereitung der Simulink-Rechnung

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

Damit kann das Anfangswertproblem so formuliert werden, wie man es in Simulink modellieren kann:

Entsprechend der Empfehlung im Kapitel "Verifizieren von Computerrechnungen" wird eine Kontrollfunktion berechnet und graphisch dargestellt. Dort wird für die Gesamtenergie, die im System zu jedem Zeitpunkt gespeichert ist, folgender Ausdruck hergeleitet (hier in dimensionsloser Form aufgeschrieben, das Null-Potenzial liegt auf der Höhe des Aufhängepunktes): Dieser komplizierte Ausdruck muss für jeden Satz der Bewegungsgrößen φ₁, φ₂, ω₁, ω₂ eines beliebigen Zeitpunkts t den gleichen Wert ergeben. Es ist die Gesamtenergie, die bereits am Anfang der Bewegung (als potenzielle Energie) vorhanden ist: Diese Funktion wird mit den ermittelten Werten für φ₁(t), φ₂(t) und deren Ableitungen berechnet und graphisch dargestellt.

Simulink-Blockschaltbild

Das nachfolgend zu sehende Simulink-Blockschaltbild bildet das oben angegebene Anfangswertproblem ab:

• Links werden die beiden Winkelbeschleunigungen hineingeführt (wo die herkommen, wird nachfolgend erläutert). Sie werden jeweils durch zwei Integratoren geschickt: Der erste Integrator (phi1pp->phi1p bzw. phi2pp->phi2p) macht aus den Winkelbeschleunigungen Winkelgeschwindigkeiten, nach Durchlaufen des zweiten Integrators (phi1p->phi1 bzw. phi2p->phi2) sind die Funktionen φ₁(t) und φ₂(t) entstanden. Diese sind bereits das Ziel der Berechnungen und werden deshalb sofort den "Scope"-Blöcken mit den Bezeichnungen phi1(t) bzw. phi2(t) zugeführt, die die Funktionen graphisch darstellen.

• Die Integratoren benötigen für ihre Arbeit jeweils die zugehörige Anfangsbedingung. Nach Doppelklick auf einen Integrator öffnet sich ein Fenster, nebenstehend zu sehen für den Integrator phi1p->phi1. Im Feld "Initial condition" muss der Anfangswert für φ₁ eingetragen werden (π/2). Die 0 ist die Standardvoreinstellung, und weil alle anderen Anfangswerte bei dieser Aufgabe diesen Wert haben, ist also mit der einen Eintragung alles erledigt.
• Alle Ausgangswerte der Integratoren werden abgezweigt und in einen "Mux-Block" gesteckt, der daraus einen Vektor u zusammensetzt, der dann am Ausgang dieses Blocks verfügbar ist. Damit stehen (neben den Parametern, die die Aufgabenstellung vorgibt) die Funktionen bereit, mit denen alle Ausdrücke, die die Differenzialgleichungen definieren, gebildet werden können. Entsprechend der Reihenfolge am Eingang des Mux-Blocks sind die Werte dann im Vektor u positioniert:
  

• Zunächst wird der Weg verfolgt, auf dem der Ausdruck im Zähler auf der rechten Seite der ersten Differenzialgleichung erzeugt wird: In dem Fcn-Block ("Function-Block") ganz oben im Block-Schaltbild (mit der Bezeichnung "b1") wird b₁ erzeugt. Nach Doppelklick auf diesen Block öffnet sich ein Fenster, nebenstehend ist zu sehen, wie dort die Funktion eingetragen wird. Dabei müssen die u-Werte verwendet werden, hier steht z. B. u(2) für die Variable φ₂.
• Für den zweiten Term im Zähler wird (wieder startend am Ausgang des Mux-Blocks) zunächst in einem Fcn-Block mit dem Namen "a12" der Ausdruck a₁₂ erzeugt. Dieser wird in einen Product-Block mit dem Namen "a12*phi2pp" eingespeist. Der andere in diesen Block einfließende Faktor ist die Winkelbeschleunigung, die aus dem Ergebnis des (bisher noch nicht betrachteten) unteren Zweigs abgezweigt wird.
• Das Ergebnis dieses Product-Blocks wird mit dem Ergebnis des Fcn-Blocks, der den Ausdruck b₁ erzeugte, in einem Sum-Block zum Zähler des Beschleunigungs-Ausdrucks zusammengeführt. Dessen Ergebnis wird schließlich noch in einem Gain-Block mit dem reziproken Wert von a₁₁ multipliziert, und damit schließt sich der Kreis: Das Ergebnis (die Winkelbeschleunigung des oberen Pendels) wird den Integratoren zugeführt.
• Auf entsprechendem Wege entsteht der Ausdruck für die Winkelbeschleunigung des unteren Pendels, wobei einige bereits berechnete Ausdrücke verwendet werden können.
• In dem großen Fcn-Block "Kontrollfunktion Tges" im unteren Bereich des Blockschaltbilds entsteht die zur Kontrolle zu erzeugende Gesamtenergie des Systems. Das Ergebnis wird unmittelbar einem Scope-Block (Name: "Tges(t)") zugeführt.
• Bevor die Rechnung gestartet werden kann, müssen noch die Parameter der Aufgabenstellung als Zahlenwerte bekanntgemacht werden, indem sie z. B. (natürlich mit den im Blockschaltbild verwendeten Namen) über das Command Window eingegeben werden. Sie sind am unteren Rand des Blockschaltbilds aufgelistet und können mit "Cut and Paste" in das Command Window übertragen werden.

Ergebnisse der Berechnung

Nachdem (z. B. wie oben beschrieben) die Problemparameter als Zahlenwerte eingegeben wurden, sollte man sich über Simulation|Configuration parameters ... die Standardeinstellungen ansehen, mindestens aber den Zeitbereich für die Berechnung einstellen:

Wenn dieses Fenster wie abgebildet aussieht, kann die Rechnung aus dem Fenster mit der Darstellung des Blockschaltbilds über das Menüangebot Simulation|Start gestartet werden. Man sieht zunächst nichts, die Ergebnisse landen in den "Scopes". Durch Doppelklick auf Scope "phi1(t)" wird die graphische Darstellung des Bewegungsgesetzes des oberen Pendels, durch Doppelklick auf Scope "phi2(t)" wird die graphische Darstellung des Bewegungsgesetzes des unteren Pendels sichtbar (es ist empfehlenswert, jeweils nach dem Öffnen des Graphik-Fensters auf das "Fernglas"-Symbol zu klicken, wodurch automatisch ein geeigneter Maßstab für die Darstellung gewählt wird):

Die Kurve der Kontrollfunktion T_ges(t) (nebenstehendes Bild, erscheint nach Doppelklick auf Scope "Tges(t)") liefert ein wichtiges Indiz für die Richtigkeit der Rechnung: Es müsste sich im Idealfall eine horizontale Gerade ergeben. Man sieht, dass die Kontrollfunktion davon einige Abweichungen zeigt. Dies sollte Anlass dafür sein, eine etwas genauere Rechnung zu erzwingen. Dazu kann man im Fenster "Configuration parameters" die Einstellungen ändern.

Die "Zacken" in der Kurve der Kontrollfunktion verschwinden jedoch nicht völlig, was vor allem daran liegt, dass nach Anklicken des "Fernglases" die optimale Einstellung des Bereichs der vertikalen Koordinatenachse gnadenlos auch die kleinsten Abweichungen visualisiert.

Für diese Aufgabe (und den vorgegebenen Zeitbereich) ändern sich die Bewegungsgesetze auch bei genauerer Rechnung aber nicht mehr nennenswert, so dass die oben zu sehenden Bewegungsgesetze als ausreichend genau angesehen werden dürfen. Aber die "nie ganz saubere Kontrollfunktion" sollte immer daran erinnern, wie riskant solche numerischen Berechnungen sind und dass man in jedem Fall mehrere Rechnungen, Verfahren und Softwareprodukte bemühen sollte, ehe man den Ergebnissen vertraut. Diese Simulink-Rechnung sollte eigentlich auch nur die Berechnung mit zwei unterschiedlichen Integrationsverfahren bestätigen.

Bemerkung zur Simulink-Berechnung

Die Benutzeroberfläche von Simulink, die im Vergleich zu vielen anderen so ganz anders aussieht, suggeriert, dass hier mit einem Blockschaltbild eine Art Regelkreis definiert wird. Tatsächlich ist das wirklich nur die Oberfläche, intern wird (natürlich!) numerisch gerechnet, und eigentlich passiert nichts anderes als nach dem Schreiben eines Matlab-Scripts, wie es zum Beispiel auf den Seiten "Doppelpendel - Berechnung mit Matlab" oder "Doppelpendel - Matlab (mit Kontrollen)" demonstriert wird.

Man kann sich das anzeigen lassen (und gegebenenfalls beeinflussen), was intern geschieht. In dem Fenster, das über das Menüangebot Simulation|Configuration Parameters ... erreicht wird (Ausschnitt ist hier zu sehen) sieht man, dass die Voreinstellung für den Solver genau der ode45-Algorithmus ist, der auch als Standard-Solver für Matlab-Scripts (und damit für die Berechnungen, auf die oben verwiesen wird) verwendet wird. Man kann sowohl diese Einstellung als auch die "MaxStep"-Einstellungen, auf die bei den Matlab-Berechnungen ausführlich eingegangen wird, verändern.

Fazit: Simulink ist nur eine andere Oberfläche für die Matlab-Funktionen.

Download

Das oben zu sehende Simulink-Blockschaltbild steht als doppelpendelsim.mdl zum Download zur Verfügung.