Simulation02


Thema Simulation
Veranstaltung ADP Laufrobotik
Semester WS 2012/13
Namen Burbach, J.-N.; Erler, P.; Hoitz, F.; Stuhlenmiller, F.; Xiaoguang, Z.; Zimmermann, C.; Zwetsch, F.
Bearbeitungsdauer 60 min
Author/Verantwortlich Philipp Erler


Die in folgenden Text erstellten Modelle sind in folgender Datei zu finden: modelle.rar

Das Hopper-Modell

Abbildung 1 zeigt das Hopper-Modell mit einer Freischnitt-Skizze. Die Kinematik, die die Masse m2 über den Muskel hält wird hier nicht dargestellt, sondern die Masse m2 direkt an die Feder SEE drangeheftet. Neben dem Fuß, der durch die Masse m1, einer Feder mit der Steifigkeit k1 und eine Dämpfung mit b1 abgebildet wird, sind an den Verbindungsstellen der Elemente, bzw. für jede Koordinate, Massen (m3 und m4)eingefügt worden. Diese dienen als Platzhalter und können später verwendet werden, wenn die einzelnen Komponenten des realen Aufbaus ausgelegt sind und somit auch deren Massen. Als Zwischenschritt zu diesem Modell wurde ein Hopper-Modell ohne m3 und m4 entwickelt. Da sich durch das Modell mit den Massen nahezu das gleiche Systemverhalten einstellt, wenn für diese sehr geringe Werte verwendet werden, wird an dieser Stelle nicht weiter auf das „Zwischen-Modell“ eingegangen. Da in der Realität jede Feder eine steifigkeitsproportionale Dämpfung besitzt, wird an deren Stellen eine Dämpfung in Höhe von b = 0.02 k vorgesehen. Im Folgenden werden die aus den Kräftegleichgewichten und Zwangsbedingungen hergeleiteten Gleichungen vollzählig aufgestellt. Aus dem Kräftesatz an den 4 Massen ergeben sich folgende Gleichungen:

(20 - 23)


Diese werden nach den Koordinaten aufgelöst und in Simulink implementiert. Für jede Koordinatenberechnung wird ein Subsystem erstellt. Desweiteren ergeben sich durch die Dämpfer und Federn folgende zum Teil schon bekannte Koppelgleichungen (Berechnung dPDE(FCE) siehe Formel 5):

(24, 25)


<imgcaption image9 |Hopper-Modell mit Freischnitt Skizze> </imgcaption>

(26 - 32)


Wie bereits erwähnt, repräsentieren die Steifigkeit k1 und die Dämpfungskonstante b1 den Fuß und entsprechen wertemäßig denen des Marco-Hopper-Designs [Kalveram u. a. (2012)]. Um die Übersicht im Simulink-Modell zu bewahren wird für jede dieser Kräfte ein Subsystem erstellt. Außerdem wir für jede Koordinate eine absolute Position festgelegt. Diese hat keinen Einfluss auf die Simulation, wird jedoch später für die Abschätzung der Federlängen verwendet (siehe Tabelle A.2 im Anhang)

Validierung

Zuerst werden nochmals die statischen Auslenkungen an den drei Koordinaten berechnet und mit den Ausgangssignalen der Simulation verglichen.

(33 - 35)


Wie in Abbildung 2 zu sehen ist, stimmen die berechneten Werten mit denen der Simulation überein. Die geringen Abweichungen ab der dritten Nachkommastelle sind auf numerische Fehler und die nicht vollständig abgeklungene Schwingung zurückzuführen.

<imgcaption image10 |Statische Auslenkung des Hopper-Modells. Verwendete Parameter siehe Tabelle A.3> </imgcaption>

Neben den statischen Auslenkungen wird auch das Sprungverhalten betrachtet, um abzuschätzen ob das dynamische Verhalten plausibel ist. In Abbildung 3 ist der Verlauf der Kraft FAE und der zugehörige Verlauf der Koordinaten dargestellt. Man sieht deutlich, dass die Koordinate ym1, an der die Kraft

<imgcaption image11 |Dynamisches Verhalten Hopper-Modell. Verwendete Parameter siehe Tabelle A.3> </imgcaption>

FAE die größte Auslenkung erfährt, gefolgt von einer Auslenkung von ym2, auf welche wiederum eine Auslenkung von y2 folgt. y1 wird zu dieser Phase aufgrund der hohen Steifigkeit lediglich leicht abgesenkt. Nach einer Aktivierung von tact = 0.25s ist die Maximalkraft FAE,max erreicht. Die Dauer der Aktivierung tac t wurde willkürlich festgelegt. Durch die Verzögerung von tverzoeg = 50ms fällt die Kraft nicht sofort auf Null herab. Nach diesem Zeitpunkt erreicht ym1 die maximale Auslenkung und kurz darauf ym2 ebenfalls. y2 steigt weiter an aufgrund der in der Masse m2 gespeicherten kinetischen Energie. Diese kinetische Energie wird von der Feder PEE abgestützt und auf die Masse m1 übertragen. Daher folgt kurz nach dem Erreichen der maximalen Auslenkung von y2 ein starker Anstieg von y1. Die unterkritisch gedämpften Koordinaten ym1, ym2 und y2 führen daraufhin eine der Eigendynamik des Systems entsprechenden Schwingung aus. Die Frequenz dieser Schwingung wird im Folgenden näher betrachtet. Die Verläufe werden mit weiteren Parameter Kombinationen betrachtet, um Erkenntnisse über den Einfluss einzelner Parameter zu erhalten und auf Korrektheit zu überprüfen. Da es dabei zu keinen unplausiblen Verhaltensweisen kam, wird an dieser Stelle nicht weiter darauf eingegangen. Die Struktur des Hoppers ist an dieser Stelle vollständig und weißt keine Fehler auf. Daher wird mit diesem Modell eine vorläufige Parametrierung durchgeführt. Dazu wird eine Auslegungsmethode verwendet, die in Rücksprache mit Häufle festgelegt wurde: Die Federsteifigkeit des SEE wird so gewählt, dass die Eigenfrequenz der Masse m2 der Frequenz entspricht, mit der gehüpft wird. Als „Hüpf-Frequenz“ wird 2Hz festgelegt, diese liegt unter der maximalen Hüpf-Frequenz des Marco-Hoppers. Daraus ergibt sich für die Steifigkeit

(35)


Zur Vereinfachung wird kSEE = 230Nm festgelegt. Die Steifigkeit kPEE soll nach Aussage von Häufle möglichst gering gewählt werden und deutlich kleiner als kSEE. Daraus ergibt sich jedoch das Problem, dass die Differenz der Auslenkungen y1 und ym1 sehr groß wird. Gerade diese soll jedoch klein werden, um beim realen Muskel einen Dämpfer einsetzen zu können, der den entsprechenden maximal Hub besitzt. Bei einem kPEE = 100Nm ergibt sich eine Hub von fast 0.5m. Bei höheren Werten wird dieser Hub geringer, jedoch zeigt sich auch, dass die Sprunghöhe des Hoppers abnimmt. Aus diesem Grund wird vorerst ein Wert von kPEE = 150Nm festgelegt, bei dem der Hub von ca. 0.4m immer noch zu hoch ist. Stattdessen wird der Einfluss der anderen Parameter untersucht. Eine Möglichkeit zur Verringerung des Hubs ist die Erhöhung des maximalen Dämpferkraft dPDE,max. Hier zeigt sich jedoch, dass die Sprunghöhe stark reduziert wird. Die gleichen Effekte entstehen durch eine Verringerung der Aktivierungszeit tact und der Verringerung der Maximalkraft FAE,max. Das Problem könnte gelöst werden, indem von die Dämpfung nicht durch einen realen Dämpfer, sondern ledlich durch die Motorsteuerung emuliert wird. Damit wäre jedoch nicht der Zielkonflikt zwischen realistischen Parametern und der Sprunghöhe beseitigt. Zudem erfordert dies je nach Auslegung des Antriebs, dass eine Dämpfung nur stattfinden kann, wenn das CE konzentrisch kontrahiert wird, aufgrund der Nichtlinearen Federkennlinie des SE. Da Häufle keine Aussage darüber macht, wie genau sich das CE bei exzentrischer Kontraktion verhält, wird für das SE eine Feder mit linearer Kennlinie angenommen. Bei der Untersuchung, wie sich diese Änderung auf das Schwingungsverhalten auswirkt, zeigt sich anhand eines Leistungsdichtespektrums (Abbildung 4) der Koordinate y2, dass diese nicht wie anfänglich angenommen mit der Frequenz von ca. 2Hz schwingt, sondern mit einer Frequenz die aus der seriellen Verschaltung von PEE und SEE approximiert werden kann:

(36)


<imgcaption image12 |Leistungsdichtespektrum der Koordinate y2 des Hopper-Modells. Verwendete Parameter siehe Tabelle A.3, Federdämpfungen auf 0 gesetzt> </imgcaption>

Der Wert fy2,approx stimmt fast mit der niedrigsten, tatsächlich auftretenden Frequenz überein. Das Problem ist jedoch, dass sich durch Hinzufügen der Dämpfung an den Stellen SE, PEE und SEE diese Frequenzen verschieben und mittels eines Leistungsdichtespektrums nicht mehr ermittelt werden können. Diese sollte jedoch möglichst exakt ermittelt werden, um daraus die optimale Aktivierungsdauer und Frequenz für das Hüpfen ableiten zu können. Dazu ist eine genaue strukturdynamische Untersuchung nötig, wie im folgenden Kapitel gezeigt wird.



Zustandsraumdarstellung des linearisierten Hopper-Modells

Bisher ist es aufgrund der vorhandenen Nichtlinearitäten in SE und DPE nicht möglich gewesen, dass Hopper-Modell analytisch zu untersuchen. Da diese in den zuvor beschriebenen Entwicklungsschritten herausgenommen wurden, verhält sich das Hopper-Modell in den Phasen des Bodenkontaktes linear. Daher wird vom Hopper-Modell eine Zustandsraumdarstellung (ZRD) erzeugt und mittels dieser die Eigenfrequenzen des Systems bestimmt. Aus den Gleichungen (20), (21), (22) und (23) wird die ZRD (38) gewonnen.



(38)


Die Gravitation ist eine von außen auf das System wirkende Erregung und hat keinen Einfluss auf die Eigendynamik des Systems. Die Kräfte FAE und FPDE werden ebenfalls als Eingangsgrößen betrachtet. Im Falle des Dämpfers PDE ist dies zulässig, wenn dieser nur bei Aktivierung Energie dissipiert. Mittels der Matlab-Funktion eig() werden die Eigenwerte des Systems bestimmt und aus deren Imaginärteil der konjugiert komplexen Eigenwerten die Eigenfrequenzen feig,n ermittelt. Für die in Tabelle A.3 eingetragenen Parameter ergeben sich folgende Werte.

Tabelle A.3 eingetragenen Parameter ergeben sich folgende Werte.

(39)


Ein Blick auf Abbildung 4 zeigt, dass diese Werte mit denen des Leistungsdichtespektrums übereinstimmen. Im folgenden können mittels dieser Analyse auch bei größeren Dämpfungen die Eigenfrequenzen exakt bestimmt werden.



Parameterermittlung

Abweichend von der von Häufle vorgeschlagenen Auslegungsmethode der Parameter wird nun SEE und PEE soweit erhöht, bis die niedrigste Eigenfrequenz des Gesamtsystems bei ungefähr 1,8Hz liegt. Die Eigenfrequenz wird etwas geringer als die vorherigen 2Hz gewählt, um eine etwas größere Sprunghöhe zu erreichen und die gesetzten Anforderungen zu erfüllen (niedrigere Sprungfrequenz erlaubt längere Aktivierung). Nach Ausprobieren verschiedener Parameterkombinationen werden die Parameter kSEE = 800Nm und kPEE = 300Nm festgelegt. Mit diesen Werten ergibt sich eine Eigenfrequenz von feig,1 = 1.7962. An dieser Stellen sollen die Parameter des CE anhand von in der Literatur angegebenen Werten ermittelt werden. Dazu werden die Berechnunggleichungen aus [Haeufle u. a. (2012)] verwendet:

(40, 41)


FAE,max = 60N wird wie zuvor so gewählt, dass es ungefähr der vierfachen Massenkraft von m2 entspricht. Die Länge LCE,opt ist die Länge, bei der der Muskel die Maximalkraft FAE,max erzeugt. Da in unserem Fall das Kraft-Längenverhältnis nicht beachtet wird, wird ein realistischer Wert von LCE,opt = 0,2m festgelegt (z. B. M. semitendinosus [Guenther (1997)]). Zur Berechnung der Muskelkonstanten werden Werte von van Soest [van Soest (1992)] (zitiert in [Guenther (1997)]) verwendet:

(42, 43)


Mit diesen Werten ergeben sich die Dämpfungsparameter RPDE = 0,2908 und DPDE,max = 81,3462 ≈ 80N. Zur Steifigkeit des SE konnten keine Informationen gefunden werden, daher wird dieser Wert unverändert beibehalten (kSE = 2000N ). Somit sind alle Parameter des CE festgelegt. Abbildung 5 und 6 zeigen das Kraft-Geschwindigkeit-Verhältnis mit den ermittelten Parametern. Ein Versuch nachzuweisen, dass es sich hierbei um ein realistisches Muskelverhalten handelt wird aus zeitlichen Gründen nicht vorgenommen. Ob ein einfach Abgleich dieser Kurven mit den zuvor ermittelten aus Abschnitt 4.1 möglich ist, ist schwer zu beurteilen. Dies liegt daran, dass Häufle für seine Versuche die Parameter eines Rattenmuskels verwendet und daher geringe Kräfte und Geschwindigkeiten zustande kommen. Folglich müssten die Kurven skaliert werden und es ist nicht bekannt, ob ein Vergleich der Ergebnisse daraufhin noch als sinnvoll zu erachten ist. Mit der Bestimmung der exakten Eigenfrequenzen feig,n und der Festlegung der Parameter ist es möglich die optimale Aktivierungsdauer und Frequenz für Hüpfen bzw. periodische Anregung zu ermitteln. Dazu wird angenommen, das die optimale Aktivierungsdauer einem Viertel der Periodendauer beträgt, die sich aus der niedrigsten Eigenfrequenz des Systems ergibt:

(44)


<imgcaption image13 |Kraft-Geschwindigkeits-Verhältnis des CE aus QRM, verwendete Parameter siehe Tabelle A.4> </imgcaption>

<imgcaption image14 |Kraft-Geschwindigkeits-Verhältnis des CE aus QRF, verwendete Parameter siehe Tabelle A.4> </imgcaption>

Um die Annahme auf Korrektheit zu überprüfen, wird die Simulation mit dieser Aktivierungsdauer und mit um den berechneten Wert variierende Dauern ausgeführt und die Sprunghöhen verglichen. Als Periodendauer wird immer das vierfache der jeweiligen Aktivierungsdauer verwendet. Diese Versuche werden einmal ohne eine Verzögerung der Kraft FAE und einmal mit einer Verzögerungzeit von tverzoeg = 30ms ausgeführt. Als Sprunghöhe wird der Wert genommen, der sich nach längerem Hüpfen einstellt. Die Ergebnisse dieser beiden Versuchsreihen sind in den Tabellen 7 und 8 zu sehen. Wie angenommen werden die größten Sprunghöhen mit dem berechneten Wert takt,opt erreicht. Der Verlauf der Koordinaten bei periodischer Aktivierung mit takt = 1,7962/4 und tverzoeg = 30ms ist in Abbildung 9.

<imgcaption image15 |Sprunghöhe in Abhängigkeit der Aktivierungsdauer t act , t verzoeg = 0, restliche Parameter siehe A.4> </imgcaption>

<imgcaption image16 |Sprunghöhe in Abhängigkeit der Aktivierungsdauer t act , t verzoeg = 30, restliche Parameter siehe A.4> </imgcaption>

Energiebetrachtung

An dieser Stelle wird eine Betrachtung der Energieverteilung bei einem Sprung des Systems durchgeführt. Dies dient der Überprüfung, ob die Eigenschaft von Sehnen, Energie zu speichern realistisch abgebildet wird. Beim Menschen konnte experimentell nachgewiesen werden, dass z. B. in der Achillessehne ca. 25% der beim Springen genutzen Energie gespeichert werden können [Lichtwark und Wilson (2005)]. Für das vorliegende System soll dies untersucht werden. Dazu wird die Simulation so umgeändert, dass der gesamte Hopper aus einer willkürlich festgelegten Höhe von 0,2m fallen gelassen wird. Beim Bodenkontakt wird ein Teil dieser potentiellen Energie in den Federn gespeichert. Der Rest wird durch die steifigkeitsproportionale Dämpfung der Federn dissipiert oder liegt als kinetische Energie vor. Der in den Federn gespeicherte Teil sollte am größten sein, wenn die Masse m2 den niedrigsten Punkt erreicht hat. Ziel ist es zu diesem Zeitpunkt die Energieanteile zu bestimmen und abzuschätzen, ob das System zu stark oder schwach gedämpft ist, im Vergleich zu einem menschlichen Muskel. Zudem kann diese Betrachtung dazu genutzt werden, das System zu validieren. So muss z. B. durch eine Erhöhung einzelner Dämpfergrade die dissipierte Energie immer ansteigen.

<imgcaption image17 |Dynamisches Verhalten Hopper-Modell. Verwendete Parameter siehe Tabelle A.3> </imgcaption>

Die Simulation wird mit denen aus dem vorherigen Abschnitt bestimmten Parametern durchgeführt (siehe Tabelle A.4). Die Abtastrate wird auf eine feste Schrittweite von tabt = 0,00001 gesetzt, um bei der Integration der durch die Dämpfer dissipierten Energie nur ausreichend kleine Fehler zu erhalten. Während der Simulation werden alle Verschiebungen, Geschwindigkeiten und Kräfte festgehalten. Wie bereits erwähnt, sind die Zustände zu zwei Zeitpunkten bzw. die dazwischen dissipierte Energie von Interesse: Der Zeitpunkt t = 0 und der Zeitpunkt, an dem die Masse m2 den niedrigsten Punkt erreicht tmin(y2). Zu Beginn, d. h. t = 0, liegt nur in den Massen gespeicherte potentielle Energie vor. Die Gesamtenergie des Systems berechnet sich zu:

(45)


Bis zum Zeitpunkt tmin(y2) ist folgender Anteil dieser potentiellen Energie umgewandelt:

(46)


Die in den Federn gespeicherte Energie zu diesem Zeitpunkt berechnet sich zu:

(47)


Somit wird folgender Anteil der potentiellen Energie in den Federn gespeichert :

(48)


Um zu kontrollieren, ob die berechneten Werte korrekt sind, werden die kinetischen Energien der Massen zum Zeitpunkt tmin(y2) berechnet und die bis zu diesem Zeitpunkt dissipierte Energie. In der Summe müssen diese und die in den Federn gespeicherte Energie EFedern,tmin(y2) der zugeführten Energie Epot,zu entsprechen. Für die kinetische Energie ergibt sich folgender Betrag:

(49)


Für die dissipierte Energie wird zuerst aus den Dämpferkräften und den Geschwindigkeiten die Leistungs berechnet:

(50)


Diese Leistungen der Dämpfer werden aufaddiert und über der Zeit integriert:

(51)


Die Summe aus dissipierter, kinetischer und in den Federn gespeicherter Energie ergibt:

(52)


Die Werte stimmen bis auf die zweite Nachkommastelle überein. Die Ursache für die geringe Abweichung ist nicht bekannt, evtl. handelt es sich um numerische Ungenauigkeiten. Die Berechnung der Energien kann somit als korrekt angenommen werden. Da ein Anteil von 62,44% der zugeführten Energie in den Federn gespeichert wird, müsste die Dämpfung noch weiter erhöht werden, um das Verhalten des Systems an das eines realen Muskels anzupassen. Da sich dies wahrscheinlich negativ auf die Sprunghöhe auswirken würde, wird keine Änderung an den Parametern vorgenommen. Eine Alternative bietet die Verwendung einer veränderten Aktivierung (siehe z. B. [Geyer u. a. (2003)]).

Danksagung

An dieser Stelle möchten wir uns für die freundliche Unterstützung von Daniel Häufle bedanken.

Quellen

Siehe Literaturverzeichnis





Weiter zu der Konzeptbildung

adp_laufrobotik/adp_2012_ws_group1/simulation02.txt · Zuletzt geändert: 25.04.2013 18:27 von Fabian Hoitz
GNU Free Documentation License 1.3
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0