Analyse nichtglatter dynamischer Systeme mit mengenorientierten Methoden am Beispiel eines Ultraschall-Stoßbohrsystems zur Erlangung des akademischen Grades eines DOKTORS DER INGENIEURWISSENSCHAFTEN(Dr.-Ing.) der Fakult¨at fu¨r Maschinenbau der Universit¨at Paderborn genehmigte DISSERTATION von Dipl.-Ing. Nicolai Neumann aus Bad Vilbel Referent: Korreferent: Tag der Einreichung: Tag der mu¨ndlichen Pru¨fung: Prof. Dr.-Ing. J¨org Wallaschek Prof. Dr. Michael Dellnitz 11. Oktober 2007 2. Juni 2008

.

Vorwort Die Arbeiten fu¨r die vorliegende Dissertationsschrift fu¨hrte ich w¨ahrend meiner dreij¨ahrigen Stipendienzeit ab 1. September 2003 als Stipendiat des PaSCo-GKs(DFG Graduiertenkolleg Wissenschaftliches Rechnen des Paderborn Institute for Scientific Computation: Anwendungsorientierte Modellierung und Algorithmenentwicklung) und wissenschaftlicher Mitarbeiter des Lehrstuhls fu¨r Mechatronik und Dynamik am Heinz Nixdorf Institut (HNI) der Universit¨at Paderborn und einem weiteren, vierten Jahr durch. Zu zweierlei Gruppierungen durfte ich mich dazuz¨ahlen-- einerseits zu meiner der Fakult¨at fu¨r Maschinenbau angeh¨orenden Fachgruppe, andererseits zum Graduiertenkolleg mit seinem Schwerpunkt in Mathematik, Ingenieurwissenschaften und Informatik. So genoss ich Betreuung und Austausch, Beratung und Diskussion in beiden Disziplinen, der Mechanik und der Mathematik, die beide wissenschaftliche Handwerkszeuge meiner Arbeit sind. Dankend m¨ochte ich die DFG(Deutsche Forschungsgemeinschaft) nennen, die meine Forschung u¨ber das PaSCo mit einem Stipendium und der Finanzierung einer wissenschaftlichen Hilfskraft fu¨r das letzte halbe Jahr angenehm gemacht hat sowie zusammen mit der Fachgruppe Mechatronik und Dynamik des HNIs die Teilnahme und Mitwirkung an Fachkonferenzen erm¨oglicht hat. Umfangreiche Laborausstattung, Werkstatt und Rechnercluster wurden mir vom HNI zur Verfu¨gung gestellt. Fachliche Konsultation erhielt ich auf der Seite des Maschinenbaus durch meinen Doktorvater und Fachgebietsleiter Prof. Dr.-Ing. J¨org Wallaschek und Prof. Dr.-Ing. Thomas Sattel. Besonders m¨ochte ich hier Thomas Sattel erw¨ahnen, der mein Interesse auf die Forschungsarbeiten im Heinz Nixdorf Institut lenkte. Er hat mich wissenschaftlich mit großem Engagement bis zuletzt durch meine Dissertationszeit begleitet und oft wichtige, motivierende Impulse gegeben. Fachliche Unterstu¨tzung fu¨r die faszinierenden mathematischen Belange meiner Arbeit erhielt ich durch den Sprecher des PaSCo und meinen Zweitgutachter Prof. Dr. Michael Dellnitz und maßgeblich auch durch Prof. Dr. Oliver Junge, Mitbegru¨nder und Entwickler der mengenorientierten numerischen Methoden und ihrer Algorithmen. Meinen Dank m¨ochte ich auch an all jene richten, die mir w¨ahrend meiner Zeit in Paderborn begegnet sind und Diskussionspartner, Ideengeber, iii

Helfer oder Ratgeber waren, im Besonderen an Jens Twiefel, mit dem ich einige mechanische Modelle und L¨osungsans¨atze er¨orterte, meine sehr angenehmen Bu¨rokollegen Bj¨orn Richter und Dr. Jia Du, Christian Potthast, Stefan Sertl, der mir das Softwarepaket GAIO(Global Analysis of Invariant Objects) mehrfach kompilierte und mit mir geduldig oft Algorithmen und Ergebnisse besprach, Dr. Kathrin Padberg, Mirko Hessel von Molo, Dr. Marcus Post, Dr. Martin Ziegler, Dr. Sina Ober-Bl¨obaum, Dr. Tobias Hemsel, Martin Liekenbr¨ocker, Wilfried Br¨ockelmann, Dr. Walter Littmann, Carmelo Quirante Kneba, Marina Kassu¨hlke, Kerstin Hille, Begona Nasarre, Axel Keller fu¨r die freundliche Hilfe mit der Bedienung des Hochleistungsrechners Arminius des PC 2 s, Axel Betanski, Dr. Markus Hohenhaus vom Rechnerbetrieb des HNI, an meinen Freund und Kommilitonen Dr. Tobias Niemz, der ein weiterbringendes Dissertationstreffen initiierte. Abschließend gebu¨hrt meiner Freundin und Verlobten Rebekka Oeters ein sehr großer Dank fu¨r ihre qualitative, wertvolle, unermu¨dliche und ausdauernde Unterstu¨tzung, motivierende Hilfe und ihr Mitdenken und Korrekturlesen. Die wichtigsten Personen aber, denen ich unermesslich danke, ohne die diese Arbeit nie entstanden w¨are, sind meine Eltern Elisabeth NeumannBeuerle und Karl-Christoph Neumann, die mir durch ihre dreißigj¨ahrige, u¨berallhin begleitende Unterstu¨tzung zu dieser Wu¨rde der Promotion verhalfen. Baunatal, im Juli 2008 Nicolai Neumann iv

Inhaltsverzeichnis 1 Einleitung 1.1 Begriffsdefinition nichtlineares und nichtglattes System.... 1.2 Dynamische Systeme nichtglatter Charakteristik....... 1.3 Analyse nichtglatter dynamischer Systeme.......... 1.4 Zielsetzung........................... 1.5 Vorgehensweise......................... 1 1 3 4 5 6 2 Mengenorientierte numerische Methoden 2.1 Invariante Mengen und globale Attraktoren.......... 2.2 Approximation relativer globaler Attraktoren......... 2.3 Aufenthaltswahrscheinlichkeiten................ 7 8 10 13 3 Ultraschall-Stoßbohrer 17 3.1 Entwicklung des Ultraschall-Stoßbohrers........... 17 3.2 Aufbau des Ultraschall-Stoßbohrers.............. 18 4 Modellierung und Analyse nichtglatter Dynamik-Stand der Technik 21 4.1 Klassische Methoden zur Analyse nichtlinearer dynamischer Systeme ............................. 22 4.1.1 Zeitreihenanalyse.................... 22 4.1.2 Phasenportrait..................... 22 4.1.3 Periodische Punkte und Fixpunkte.......... 23 4.1.4 Ljapunov-Stabilit¨at und Ljapunov-Exponent..... 23 4.1.5 Poincar´e-Abbildung................... 24 4.1.6 Verzweigungs- oder Bifurkationsdiagramm...... 24 4.1.7 Einzugsgebiet...................... 24 4.2 Zellabbildungsmethoden.................... 25 4.3 Modellierung von Schwingstoßsystemen............ 27 5 Modellierung 5.1 Herleitung eines Modells 2. Ordnung............. 5.1.1 Zeitdiskrete Modellierung............... 5.1.2 Herleitung zeitdiskreter Bewegungsgleichungen.... 5.1.3 Notwendigkeit transzendenter Bewegungsgleichungen 5.1.4 Mehrfachst¨oße...................... 31 34 34 36 38 39 v

INHALTSVERZEICHNIS 5.2 Herleitung eines verfeinerten zeitdiskreten Modells...... 40 5.2.1 Aufstellen des verfeinerten Modells 4. Ordnung... 41 5.2.2 Berechnung des Stoßzeitpunktes t k+1 ......... 43 5.2.3 Berechnen der verbleibenden Zustandsgr¨oßen..... 46 5.2.4 Zusammenfassung des Stoßkontaktmodells 4. Ordnung 48 5.3 Algorithmus zum L¨osen transzendenter Stoßgleichungen... 49 6 Verhalten des Stoßbohrers 6.1 Experimentelle Untersuchungen am Schwingstoßpru¨fstand. 6.1.1 Aufbau des Schwingstoßpru¨fstands.......... 6.1.2 Messtechnik....................... 6.1.3 Stoßzahlen aus Labormessungen............ 6.2 Simulation des Modells 2. Ordnung.............. 6.2.1 Analytische Fixpunktbestimmung........... 6.2.2 Periodische und chaotische L¨osungen......... 6.2.3 Bifurkationsanalyse................... 6.3 Experimentelle Untersuchungen am Stoßbohrer....... 6.3.1 Pru¨fstand fu¨r Messungen am Stoßbohrer....... 6.3.2 Messergebnisse am Stoßbohrer............. 6.4 Simulation des Modells 4. Ordnung.............. 6.5 Vergleich von Modellen und Messungen............ 53 53 53 56 61 63 63 65 68 73 73 75 75 78 7 Analyse mit mengenorientierten Methoden 83 7.1 Analyse des Stoß-Modells 2. Ordnung............. 84 7.1.1 Der zylindrische Zustandsraum............ 84 7.1.2 Bestimmung eines relativen globalen Attraktors... 86 7.1.3 Berechnung der Aufenthaltswahrscheinlichkeiten... 87 7.1.4 Analyseergebnisse aus Aufenthaltswahrscheinlichkeiten 89 7.1.5 Analyseergebnisse aus Einzugsgebieten........ 94 7.2 Analyse des Stoß-Modells 4. Ordnung............. 97 7.2.1 Approximation des relativen globalen Attraktors... 97 7.2.2 Analyseergebnisse und Diskussion........... 99 8 Zusammenfassung, Diskussion und Ausblick 105 8.1 Zusammenfassung........................ 105 8.2 Diskussion der Ergebnisse................... 106 8.3 Ausblick............................. 108 vi

Kapitel 1 Einleitung Die vorliegende Dissertationsschrift betrachtet die Anwendung und den EinesaintezndBerei"trmaegn, gdeinesoeriMenettiherotdeennniunmdeernisicnhgeenniMeuertwhiosdseenns"c.hDafatldicuhrcehn leistet sie Bereichen st¨arker zum Einsatz zu bringen und ihren Bekanntheitsgrad zu erh¨ohen. Exemplarisch wird im Rahmen dieser Arbeit ein System mit nichtglatten Eigenschaften aus der realen Ingenieurwelt herangezogen. Hierbei handelt es sich um ein Ultraschall-Stoßbohrsystem, welches durch seine nichtlineare Dynamik, die durch die auftretenden Kontaktwechselvorg¨ange insbesondere nichtglatten Charakter besitzt, ein Beispiel aus einer herausragenden Klasse von Systemen darstellt. In Kapitel 1.1 erkl¨are ich zun¨achst die Begriffe linear/ nichtlinear und die Unterscheidung in glatt/ nichtglatt-- fu¨r Leser, denen diese Termini im Kontext dynamischer Systeme nicht gel¨aufig sind. Danach gebe ich einen U¨ berblick u¨ber die spezielle Gruppe nichtglatter dynamischer Systeme und gehe auf deren Analysem¨oglichkeiten ein. Schließlich pr¨asentiere ich in den weiteren Unterkapiteln die Zielsetzung und Vorgehensweise und erl¨autere sukzessive den Aufbau der folgenden Arbeit. 1.1 Begriffsdefinition nichtlineares und nichtglattes System Linear/ nichtlinear In einem linearen Modell, das ein physikalisches Ph¨anomen oder einen physikalischen Zusammenhang beschreibt, h¨angt eine Systemgr¨oße-z. B. die mechanische Spannung in einem Bauteil-- in linearer Weise, also proportional von einer zweiten Systemgr¨oße-- etwa der Dehnung -- ab(siehe Bild 1.1). Diese Abh¨angigkeit l¨asst sich in Form linearer Modellgleichungen mathematisch beschreiben. In einem kartesischen Koordinatensystem werden lineare Abh¨angigkeiten durch gerade Linien dargestellt. Nichtlinearit¨aten werden nicht durch lineare, sondern durch gekru¨mmte, also z. B. in ihrer Steigung zunehmende Kurven beschrieben. Zum Basiswissen der Ingenieurwissenschaften geh¨oren heute g¨angige 1

KAPITEL 1. EINLEITUNG Techniken, um das Verhalten solcher linearen Modelle untersuchen und simulieren zu k¨onnen. Vergleiche mit Labormessungen, die das Verhalten eines realen, physikalischen Systems wiedergeben, weisen jedoch in vielen F¨allen die Existenz von Ph¨anomenen auf, die mit linearen, die Realit¨at vereinfachenden Modellen, nicht abgebildet werden k¨onnen. Tragen diese Ph¨anomene wesentlich zur Funktion eines Mechanismus bei oder muss man ihr Auftreten gezielt beeinflussen oder gar verhindern k¨onnen, so wird es unerl¨asslich, die auf sogenannten Nichtlinearit¨aten beruhenden Eigenschaften in die Modellgleichungen mit aufzunehmen. Nichtlinear bedeutet allgemein, dass es keinen linearen Zusammenhang zwischen verschiedenen physikalischen Systemgr¨oßen und ihren Ableitungen gibt: eine Systemgr¨oße ist oder Integrale oder nicht einer "aenindfearcehn" das Vielfache Systemgr¨oße, einer ihrer Ableitungen sondern sie, eine ihrer Ableitungen oder Integrale kann in Differentialgleichungen als Argument einer nichtlinearen Funktion(z. B. Polynom zweiter oder h¨oherer Ordnung, Exponential-, Logarithmus- oder trigonometrische Funktion) auftreten oder z. B. in einem Produkt von Zustandsgr¨oßen mit deren Ableitungen oder Integralen enthalten sein. Als eine von vielen m¨oglichen Beispielen hierzu sei die Arbeit von Neumann[Neumann 2002] genannt, in der lange Zeit nicht beru¨cksichtigte, nichtlineare Effekte in den L¨angsschwingungen axial polarisierter, piezoelektrischer Keramiken experimentell nachgewiesen und in Modellgleichungen mit nichtlinearen Termen umgesetzt wurden. Bild 1.1: Beispiel fu¨r einen linearen Zusammenhang von Systemgr¨oßen: Spannungs-Dehnungs-Diagramm(elastischer Bereich) Glatt/ nichtglatt Auf einen Spezialfall nichtlinearer Systeme soll innerhalb dieser Arbeit besonderes Gewicht gelegt werden. Es sind sogenannte nichtglatte dynamische Systeme. In glatten Systemen l¨asst sich das Zusammenwirken von Systemgr¨oßen durch glatte Kurven veranschaulichen. Glatt bedeutet, die Kurven sind kontinuierlich, d. h. zusammenh¨angend, und stetig differenzierbar, d. h. ohne Knick. Die Charakteristika nichtglatter Systeme dagegen k¨onnen Knicke oder gar Spru¨nge enthalten. Mathematisch definierend ist eine Funktion innerhalb eines Intervalls dann glatt, wenn alle ihre Ableitungen innerhalb diesen Intervalls kontinuierlich sind, d. h. keinen Sprung besitzen. Glatte Funktionen geh¨oren somit 2

Position Geschwindigkeit 1.2. DYNAMISCHE SYSTEME NICHTGLATTER CHARAKTERISTIK der Differenzierbarkeitsklasse C an. Ein offensichtliches Beispiel fu¨r eine Vielzahl nichtglatter Systeme sind s¨amtliche Mechanismen, in denen Stoßprozesse stattfinden. W¨ahrend eines Stoßes ¨andert ein K¨orper abrupt seinen Bewegungszustand(siehe Bild 1.2). Seine Geschwindigkeit macht einen Sprung. Der Zeitverlauf der Position des K¨orpers erh¨alt dabei einen Knick. Zeit Zeit Bild 1.2: Sprunghafte A¨ nderung der Geschwindigkeit eines K¨orpers 1.2 Dynamische Systeme nichtglatter Charakteristik Nahezu alle Systeme, die heutzutage in Forschung und Wissenschaft-- besonders in den Ingenieurwissenschaften-- eine Rolle spielen, sind durch nichtlineare Eigenschaften gepr¨agt. In einigen F¨allen ist es legitim, diese Nichtlinearit¨aten zu vernachl¨assigen. H¨aufig ist es jedoch sinnvoll, sie in die Modellierung der Systeme mit aufzunehmen, wenn sie deren qualitatives Verhalten entscheidend mitbestimmen. Eine besondere Untergruppe der Systeme mit nichtlinearem Verhalten bilden solche Systeme, in denen nichtglatte Eigenschaften auftreten. Ursachen dafu¨r k¨onnen die oben erw¨ahnten Kontaktwechsel sein. Diese k¨onnen erwu¨nscht(z. B. in Stoßbohrern) oder unerwu¨nscht(z. B. Zahneingriffst¨oße in Getrieben) sein. Oder man denke an Tankschiffe, die an Hochseebojen vor O¨ lbohrplattformen vertaut werden. Periodische Anregung durch vorbeiziehende Wellenzu¨ge und die wechselnde Elastizit¨at der Taue k¨onnen zu unerwarteten und unerwu¨nschten Bewegungen fu¨hren. Die Steifigkeit der Taue h¨angt hierbei in sprunghafter Weise davon ab, ob sie unter Spannung stehen oder schlaff durchh¨angen. Der U¨ bergang zwischen diesen beiden benachbarten Zust¨anden ist durch einen Knick in der Elastizit¨atskurve der Vertauung (Ru¨ckstellkraft u¨ber Auslenkung) gekennzeichnet. 3

KAPITEL 1. EINLEITUNG Eine weitere Gruppe von Systemen mit sehr offensichtlichen nichtglatten Eigenschaften sind all jene, in denen Stoßvorg¨ange wesentlich zur Funktionsweise beitragen, wie etwa s¨amtliche h¨ammernden Schlagbohrprozesse. Hier findet w¨ahrend des Stoßes ein U¨ bergang einer freien(evtl. linearen) Bewegung hin zu einem Kontakt statt. Je nach Art der Modellbildung l¨asst sich dieses Verhalten wie beim Hochseetankschiff z. B. durch eine Steifigkeit nachbilden, die in Abh¨angigkeit einer Systemgr¨oße-- hier der Auslenkung -- ihren Wert sprunghaft ¨andert. H¨aufig findet man in der Literatur zur Beschreibung solcher St¨oße den Einsatz der Stoßzahlhypothese nach Newton, nach der die Geschwindigkeit des Stoßk¨orpers oder beider Stoßpartner im Stoßaugenblick schlagartig ihren Wert ¨andert. 1.3 Analyse nichtglatter dynamischer Systeme Im Prozess der Entwicklung, Gestaltung, Auslegung und Optimierung von Systemen aus der Mechanik, dem Maschinenbau und der Elektrotechnik ist stets ein eingehendes Verst¨andnis des Systemverhaltens unerl¨asslich. Hierfu¨r sind Schritte zum Modellieren des realen Verhaltens und schließlich zum Analysieren wesentliche und immer wiederkehrende Bestandteile. Zur Simulation nichtlinearen Verhaltens muss von Fall zu Fall abgewogen werden, welche der bew¨ahrten Methoden zur L¨osung nichtlinearer Differentialgleichungen zur Anwendung kommen k¨onnen oder ob es zul¨assig ist, Nichtlinearit¨aten durch Linearisieren zu vernachl¨assigen. Oft wird die letzte Vorgehensweise verwendet, da fu¨r die Analyse linearer Systeme robuste und bew¨ahrte Techniken vorliegen. Hierbei nimmt der Ingenieur jedoch in Kauf, wesentliche Eigenschaften seines Systems, die allein auf die Nichtlinearit¨aten zuru¨ckzufu¨hren sind, zu verlieren. Noch schwieriger wird es, wenn ein Extremfall von nichtlinearen Kennkurven das Systemverhalten beeinflusst, n¨amlich wenn nichtglatte Effekte mit im Spiel sind: einige Methoden, die fu¨r die Untersuchung nichtlinearen Verhaltens anwendbar sind, sind nur auf glatten Teilbereichen einsetzbar. Zwar besteht die M¨oglichkeit, nichtglatte Kennlinien z. B. durch die Arkustangens-Funktion 1 zu gl¨atten, doch solches Vorgehen kann zum Verlieren oder Verf¨alschen wesentlicher Effekte fu¨hren. Die Reihe der zur Verfu¨gung stehenden Methoden wurde in den letzten Jahren um die sogenannten mengenorientierten numerischen Methoden erweitert. Sie ¨ahneln den Zellabbildungsmethoden, die Hsu seit den 80er Jahren in Berkeley entwickelt[Hsu 1980]. Hsu behandelt die Zustandsgr¨oßen nicht als Kontinua, sondern als diskrete Gr¨oßen, die er in Zellen einteilt. Einen Unterraum des n dimensionalen Zustandsraumes diskretisiert er in 1 Die Funktion f ( x )= 2 lim a arctan ax konvergiert punktweise gegen die Vorzeichenfunktion signum( x ), die an der Stelle x = 0 von- 1 nach 1 springt. Fu¨r große a gibt 2 arctan ax n¨aherungsweise die Sprungfunktion in glatter Weise wieder und besitzt alle Ableitungen, geh¨ort also der Differenzierbarkeitsklasse C an. 4

1.4. ZIELSETZUNG viele, kleine n dimensionale Zellen. Die Zellabbildungsmethoden werden in Unterkapitel 4.2 vorgestellt. Die mengenorientierten Methoden liefern eine vielversprechende Weiterentwicklung des Konzepts von Hsu. Da sie durch einen effektiven Unterteilungsalgorithmus nur denjenigen Zustandsbereich fein diskretisieren, in dem die Dynamik des Systems stattfindet, verh¨alt sich ihr Rechenaufwand nicht wie die Anzahl der Zustandsraumdimensionen, sondern wird durch die (unter Umst¨anden fraktale) Dimensionsgr¨oße der Mannigfaltigkeiten 2 oder Attraktoren des Systems bestimmt. Dies kann zu einer erheblichen Rechenersparnis fu¨hren. Dennoch entfalten die mengenorientierten Methoden ihr Potential besonders mit dem Aufkommen immer leistungsf¨ahigerer Rechenmaschinen. In die Grundlagen dieser Methoden wird Kapitel 2 einfu¨hren. 1.4 Zielsetzung Exemplarisch fu¨r den oben vorgestellten Systemtypus mit nichtglattem Verhalten soll im Rahmen dieser Arbeit ein spezieller Ultraschall-Stoßbohrer untersucht werden. Das allein auf Stoßvorg¨angen beruhende Bohrkonzept ist auf eine große Klasse von Systemen mit Stoß- oder nichtglatten Vorg¨angen u¨bertragbar. Der Bohrer, der erst vor wenigen Jahren entwickelt wurde, wird in Kapitel 3 im Detail vorgestellt. Obgleich seine Funktion unter Laborversuchen nachgewiesen werden konnte, gibt es bislang nur wenige Arbeiten, die sich auf theoretischer Ebene mit seinem Funktionsprinzip besch¨aftigen[Badescu et al. 2005],[Neumann et al. 2007], [Potthast et al. 2007a]. Ein tiefgehendes Verst¨andnis der mechanischen Prinzipien, die die Funktion des Stoßbohrsystems ausmachen, ist bislang nicht bekannt. Mit mengenorientierten numerischen Methoden, die in Kapitel 2 beschrieben werden, soll ein Beitrag in Richtung eines tieferen Verst¨andnisses des Funktionsprinzips erreicht werden. Ein Bestandteil dieser Arbeit ist, neben der Einfu¨hrung und Vorstellung der relativ jungen Methoden, ihre Anwendung auf das vorgeschlagene, nichtglatte Stoßsystem zu demonstrieren und ihren Beitrag zum globalen Systemverst¨andnis eines Bohrprozesses zu betrachten. In Erg¨anzung der Analyse mit mengenorientierten Methoden werden dazu Vorgehensweisen angewandt, die in der Untersuchung nichtlinearer Modelle g¨angig sind. Hierzu wird eine Bifurkationsanalyse der Dynamik in Abh¨angigkeit von Systemparametern durchgefu¨hrt. Ziel ist es, Fixpunkte zu identifizieren sowie deren Periodizit¨at und Stabilit¨atsverhalten zu berechnen. Weiterhin ist zu pru¨fen, ob es Bereiche mit koexistierenden L¨osungen gibt, fu¨r die ggf. Einzugsgebiete zu berechnen sind. Fu¨r die angestrebten Untersuchungen ist die Dynamik einer freien Stoßmasse als wesentliche Komponente aus dem Inneren des Bohrsystems zu analysieren. In realen Laborversuchen sind die theoretisch gewonnenen Ergebnisse schließlich mit experimentellen Resultaten zu vergleichen. 2 die Mannigfaltigkeit: Der Begriff der Mannigfaltigkeit verallgemeinert den Begriff der Fl¨ache im Raum, welche man sich als verbogenes Stu¨ck eines deformierbaren Materials vorstellen kann, wie z. B. eine Kugel- oder Torusoberfl¨ache.[Brockhaus 1989] 5

KAPITEL 1. EINLEITUNG 1.5 Vorgehensweise Die weiteren Kapitel sind wie folgt aufgebaut: Kapitel 4 gibt zun¨achst einen U¨ berblick u¨ber klassische Analysetechniken fu¨r nichtlineare Systeme, die langj¨ahrige Forschungsthemen in Mathematik und Mechanik sind und ebenfalls fu¨r nichtglatte Systeme Bedeutung haben. Es enth¨alt daraufhin Erl¨auterungen zu Arbeiten, die der modernen Wissenschaft der Systemanalyse entstammen und solche, die speziell fu¨r nichtglatte Systeme entwickelt wurden. Schließlich beschreibt Kapitel 4 Arbeiten, die wichtige Ans¨atze zur modellhaften Formulierung von Schwingstoßsystemen geben. Bevor eine numerische, auf mathematischen Algorithmen basierende Methode auf ein reales System angewandt werden kann, ist es n¨otig, die Realit¨at bis auf ihre wesentlichen Bestandteile zu abstrahieren und mittels mathematischer Gleichungen zu formulieren. Dieser Schritt der Modellierung geschieht fu¨r die Dynamik im hier untersuchten Bohrprozess in Kapitel 5 . Es werden zwei Modelle unterschiedlicher Verfeinerungsstufen und Ordnungen vorgeschlagen. Modellsimulationen sind Thema in Kapitel 6 . Hier wird das dynamische Verhalten in der Theorie anhand von periodischen und chaotischen Zeitreihen und in Bifurkationsdiagrammen besprochen. Der Theorie werden danach experimentelle Studien an vergleichbaren Laborsystemen gegenu¨bergestellt, um die Gu¨te der Modellbildung bewerten und Annahmen validieren zu k¨onnen. In Kapitel 7 schließlich wird die eigentliche Analyse der Modelle vollzogen: die mengenorientierten numerischen Methoden demonstrieren ihr Potential im Einsatz auf nichtglatte Systeme. Dazu werden relative globale Attraktoren berechnet, die die statistische Wahrscheinlichkeitsverteilung fu¨r das Auftreten s¨amtlicher m¨oglicher Zust¨ande der Attraktoren wiedergeben. Dies vermittelt ein globales Bild der den Systemen zugrundeliegenden Dynamik. Enden wird die vorliegende Arbeit mit einer abschließenden Zusammenfassung und Bewertung des Beitrags der eingesetzten Methoden auf nichtglatte Systeme und einem Ausblick in Kapitel 8 . 6

Kapitel 2 Mengenorientierte numerische Methoden Ein wichtiger Bestandteil in Entwicklungsprozessen des Maschinenbaus ist gepr¨agt vom Verst¨andnis des dynamischen Verhaltens mechanischer Systeme und der M¨oglichkeit, diese mathematisch simulieren zu k¨onnen. Fru¨her gab es kaum andere M¨oglichkeiten, als durch Linearisieren vereinfachte Abbildungen der Wirklichkeit mathematisch zu untersuchen. Mitunter ist es jedoch entscheidend, Nichtlinearit¨aten zu betrachten. m¨ogEliicnhekeiRteenihefu¨retanbiclihetrltienrear"eklaSsyssistcehmeer" Analyse- und Betrachtungssind im Wesentlichen Zeitreihenanalysen, Fixpunktsuche, Stabilit¨atsbewertung, Phasenportraits, Poincar´e-Abbildungen und Bifurkationsanalysen. Sie werden in Unterkapitel 4.1 zusammenfassend dargestellt. Als Ergebnis der Analyse eines nichtlinearen Systems liefern sie wichtige und wesentliche Informationen u¨ber Gleichgewichtspunkte und periodische L¨osungen. Wenn Systeme irregul¨ares bzw. nichtperiodisches(also chaotisches) Verhalten aufweisen, ist es sinnvoll, statistische Aussagen treffen zu k¨onnen. Fragen wie Mit welcher Wahrscheinlichkeit wird eine theoretisch m¨ogliche L¨osungsform in der Realit¨at tats¨achlich angenommen? oder-- vorausgesetzt, die dynamische Entwicklung des Systems findet wissentlich innerhalb eines Attraktors statt-Mit welcher Wahrscheinlichkeit wird ein Systemzustand innerhalb des Attraktors angenommen? k¨onnen nicht allein durch die oben erw¨ahnten Methoden bearbeitet werden. Es kommt vor, dass zur alleinigen Erlangung von erstem Systemverst¨andnis oder gar als grundlegende Analysetechnik Langzeitsimulationen anhand eines Modells praktiziert werden. Ausgehend von einigen Anfangszust¨anden berechnet man einige wenige Trajektorien. In Systemen irregul¨aren oder chaotischen Verhaltens fu¨hren Langzeitintegrationen aufgrund begrenzter, 7

KAPITEL 2. MENGENORIENTIERTE NUMERISCHE METHODEN numerischer Rechengenauigkeiten jedoch m¨oglicherweise auf fehlerhafte Resultate. Eine vielversprechende Alternative und Erg¨anzung bieten die sogenannten mengenorientierten numerischen Methoden. Wie bei der Zellabbildungsmethode[Hsu 1987] basieren mengenorientierte Methoden auf Kurzzeitsimulationen sehr vieler Trajektorien. Hierin besteht ein wesentlicher Unterschied zu anderen Methoden, bei denen Langzeitintegrationen auf fehleranf¨allige Trajektorien fu¨hren. Hinsichtlich der Zellabbildungsmethoden unterscheiden sich mengenorientierte Methoden unter anderem durch algorithmische Erweiterungen, die die Recheneffizienz steigern: der Rechenaufwand h¨angt hier nicht wie sonst von der Dimension des zugrundeliegenden Zustandsraumes ab, sondern von der(unter Umst¨anden fraktalen, d. h. nicht-ganzzahligen) Dimension der zu bestimmenden Mannigfaltigkeit. Im Folgenden wird eine Einfu¨hrung in die Techniken und Algorithmen gegeben, die den Kern der mengenorientierten numerischen Methoden bilden. Entwickelt und implementiert wurden die den Methoden zugrundeliegenden Algorithmen unter anderem am Lehrstuhl fu¨r Angewandte Mathematik der Universit¨at Paderborn, siehe hierzu[Dellnitz & Junge 2002]. Ausfu¨hrliche Herleitungen und Beweise finden sich außerdem in[Dellnitz, Froyland, Junge 2001] und [Dellnitz & Hohmann 1997]. Zur Untersuchung dynamischer Systeme aus der Ingenieurspraxis wurden die Methoden zum Beispiel fu¨r Schienenfahrzeuge in[Goldschmidt 2003] verwendet. Die folgenden Abschnitte stellen die Konzepte der mengenorientierten Methoden vor. 2.1 Invariante Mengen und globale Attraktoren Im Zuge der modellbasierten Entwicklung eines nichtlinearen dynamischen Systems ist es wichtig, Kenntnisse u¨ber dessen globales Langzeitverhalten zu erlangen. Dies bezieht das Auffinden von Gleichgewichtszust¨anden, periodischen oder chaotischen L¨osungen sowie deren Einzugsgebiete mit ein. Im Falle eines schwach ged¨ampften Systems, wie es bei vielen Stoß-Schwing-Anwendungen zu finden ist, kann die Bestimmung von periodischen Punkten schwierig sein. Aus topologischer Sicht lassen sich sowohl Gleichgewichtspunkte als auch periodische und chaotische L¨osungen den invarianten h¨angen nur von Mengen zuordnen. Systemparametern ab "uInnvdarviearn¨atn"debrendesiucthetn, icdhite Mengen aufgrund der dynamischen Entwicklung des Systems mit der Zeit. Attraktor Der Begriff des Attraktors wird hier vorgestellt: ein Attraktor beschreibt einen Bereich oder eine Untermenge eines Zustandsraumes, zu dem sich 8

2.1. INVARIANTE MENGEN UND GLOBALE ATTRAKTOREN die Systemdynamik im Langzeitverhalten hinentwickelt. Bezu¨glich der zeitlichen Evolution eines Systems ist ein Attraktor eine invariante Menge: wenn der Systemzustand den Attraktor einmal erreicht hat, wird er ihn unter dem Einfluss der zugrundeliegenden Dynamik nie wieder verlassen. Daher stellt das Bestimmen von Attraktoren ein wesentliches Mittel dar, Aussagen u¨ber das Langzeitverhalten eines Systems zu treffen. Dynamisches System Mathematisch l¨asst sich ein autonomes 1 , zeitkontinuierliches, dynamisches System durch ein System von gew¨ohnlichen Differentialgleichungen beschreiben: d x ( t ) d t = g ( x ( t )) , g : R n R n , x : R R n . (2.1) Darin wird g als Vektorfeld bezeichnet, x ist der Zustandsvektor, der den Systemzustand zur Zeit t angibt und R n ist der n -dimensionale Zustandsraum stand (mauitchVo"rPahnascshenreriatuemn"degreZnaeintnet)n,twinicnkeerlhta. lIbstdGesls.e(n2.s1i)chandaelrytSiyscshtenmiczhutl¨osbar, so l¨asst sich die Zustandsentwicklung durch numerische Integration von(2.1) u¨ber die Zeit berechnen. Dann wird die zugrundeliegende Dynamik mittels eines zeitdiskreten, dynamischen Systems T der Form x k+1 = T ( x k ) , k = 0 , 1 , 2 ,..., T : R n R n , (2.2) angen¨ahert, wenn eine explizite Integrationsmethode gew¨ahlt wurde. Der Anfangszustand ist x 0 R n , und x k = x ( t k ) beschreibt den Zustand zum Zeitschritt t k . Invariante Menge, Einzugsgebiet und Attraktor Fu¨r eine Abbildung T kann der Begriff der invarianten Menge, des Einzugsgebietes und des Attraktors entsprechend[Dellnitz & Junge 2002] wie folgt definiert werden: Definition invariante Menge und Attraktor: · Eine Untermenge A R n wird invariant genannt, wenn T ( A )= T 1 ( A )= A . · Wenn zu A ein Einzugsgebiet D A existiert mit A D A und D A R n , dann wird A Attraktor genannt. · tDraaskt E or i s nz is u t g d sg ie eb M ie e t ng(eendgelr. je"nbiagseinn of attraction") eines Anfangsbedingungen, Atvon 1 autonom: die Anregung eines Systems(meistens auf der rechten Seite der Differentialgleichung, die die Dynamik eines Systems beschreibt, notiert) h¨angt nicht explizit von der Zeit t ab. 9

KAPITEL 2. MENGENORIENTIERTE NUMERISCHE METHODEN denen aus sich das Langzeitsystemverhalten in den Attraktor hinein bewegt. Formal bildet die Vereinigungsmenge aller Ru¨ckw¨artsabbildungen der fundamentalen Umgebung U das Einzugsgebiet des Attraktors: D A = k N T k ( U ) . · Falls D A = R n , dann wird A globaler Attraktor genannt. Zur Bestimmung eines Attraktors kann R n wegen seiner Unbegrenztheit numerisch nie vollst¨andig beru¨cksichtigt werden. Anstelle dessen wird der Zustandsraum auf eine interessierende Untermenge Q R n begrenzt und der [Dellnitz, relative globale Froyland, Junge Attraktor 2001] und (["DReGllnAit"z) A Q eingefu¨hrt & Junge 2002]): (siehe auch Definition relativer globaler Attraktor: A Q ( T )-- der globale Attraktor eines dynamischen Systems T relativ zu einer interessierenden, gegebenen Menge Q -- ist definiert als Schnittmenge aller Vorw¨arts-Abbildungen des Urbildes Q , bzw. A Q ( T ):= T k ( Q ) . k=0 (2.3) b"ImSietta,Fratobulgofexdn"adsgeenanlsawnAinrudts.gSdainieegssiemniteeenirngeesrsediecihreetAnecdlkgeiogrUibtnhetmgerreemnnezantnegsgee n w Qadn"idmAt enwnfeasrinodgnesna-"l,edsoieGdieenrden folgenden Abschnitten vorgestellt werden. 2.2 Approximation relativer globaler Attraktoren Dieser Abschnitt zeigt, wie der relative globale Attraktor zu einem gegebenen, zeitdiskreten, dynamischen System T numerisch ermittelt wird. Dazu wird die Idee des mengenorientierten Vorgehens erl¨autert, welche die Basis der mengenorientierten numerischen Methoden bildet. Ziel ist es, den Attraktor mit einer Menge von kleinen, rechteckigen sog.( n -dimensionalen) Boxen n¨aherungsweise zu u¨berdecken. Der erste Schritt des Vorgehens auf dem Weg zum relativen globalen Attraktor ist die Wahl eines rechteckigen Teilgebietes Q aus dem zugrundeliegenden Zustandsraum R n . Das Teilgebiet Q werden wir im weiteren Verlauf mit Startbox bezeichnen. Nur innerhalb dieser Startbox wird die Analyse stattfinden. Beginnt man eine Analyse ganz ohne Vorwissen u¨ber das System, so ist es ratsam, zun¨achst eine recht große Startbox zu w¨ahlen. Im Folgenden wird ein Box-Unterteilungsalgorithmus vorgestellt [Dellnitz, Froyland, Junge 2001]. Der Algorithmus besteht aus zwei Schritten, welche einige Male iteriert 2 werden. Begonnen wird mit der rechteckigen Startbox Q (0) = Q . Die zwei Schritte sind: 2 iterieren: wiederholen, wiederholt anwenden. 10

2.2. APPROXIMATION RELATIVER GLOBALER ATTRAKTOREN 1. Unterteilung: Gegeben sei der i -te iterative Unterteilungsschritt, nach dem N ( i) Boxen B k ( i) , k = 1 , 2 ,..., N ( i) vorliegen, und es gelte N ( i ) B k ( i) = Q ( i)  Q (0) R n . k=1 (2.4) Halbiere im n¨achsten Schritt i + 1 jede Box B k ( i) . Dies ergibt 2 N ( i) kleinere, verfeinerte Boxen fu¨r die wieder gilt: 2 N ( i ) B k ( i+1) = Q ( i) . k=1 (2.5) 2. Auswahl: Bestimme das Urbild T 1 ( B k ( i+1) ), k = 1 ,..., 2 N ( i) jeder der verfeinerten Boxen durch Anwendung der invertierten Abbildung 3 . L¨osche daraufhin all jene Boxen B k ( i+1) , deren Urbild die gegenw¨artige Menge Q ( i) an Boxen nicht schneidet. Die Menge Q ( i+1) fu¨r den n¨achsten Schritt i + 1 ergibt sich aus den u¨brig bleibenden Boxen: 2 N ( i ) Q ( i+1) = T 1 k=1 B k ( i+1)  Q ( i) . (2.6) Bei der Unterteilung der n -dimensionalen Boxen B k ( i) wird die Halbierungsebene entlang der n Dimensionen des Zustandsraumes zyklisch permutiert. Beispielsweise im Falle eines 2-dimensionalen Systems mit den Zustandsgr¨oßen x und y wird die Halbierung der rechteckigen Boxen zun¨achst in x -Richtung, dann in y -Richtung und im n¨achsten Iterationsschritt wieder in x -Richtung und so weiter vorgenommen. Die Anzahl der Iterationsschritte i sollte(muss aber nicht) ein Vielfaches der Dimensionszahl n betragen, so dass die Unterteilung in jede Koordinatenrichtung gleich oft geschieht. Im Gegensatz zu anderen Zellteilungsverfahren h¨angt der Rechenaufwand des hier vorgestellten Verfahrens nicht alleine von der Anzahl der Dimensionen ab, sondern kann in Abh¨angigkeit von der zugrundeliegenden Dynamik geringer sein. Ein Attraktor mit schlanker, fein ver¨astelter Struktur kann daher mit geringer Rechenzeit angen¨ahert werden, da nach jedem Unterteilungsschritt nur wenige Boxen behalten werden, die weiter unterteilt werden mu¨ssen. Bild 2.1 zeigt eine solche Anordnung von Boxen, die durch den Unterteilungsalgorithmus entstand. Diese Menge von Boxen u¨berdeckt den gesuchten relativen globalen Attraktor eines dynamischen Systems, welches sp¨ater in Kapitel 3 vorgestellt wird. Fu¨r die Feinheit der Rechteck-Boxen im gegebenen Beispiel waren 14 Unterteilungsschritte n¨otig. In diesem Beispiel wird 3 Die Existenz der inversen Abbildung T 1 stellt in der Praxis keine Voraussetzung dar. Die Auswahl der zu l¨oschenden Boxen wird vielmehr derart vorgenommen, dass zun¨achst die Bilder jeder der verfeinerten Boxen berechnet werden. Daraufhin werden all solche Boxen gel¨oscht, die sich mit keinem der Bilder schneiden. 11

KAPITEL 2. MENGENORIENTIERTE NUMERISCHE METHODEN der Rechenvorteil des mengenorientierten Ansatzes gegenu¨ber der Zellabbildungsmethode deutlich. In der Zellabbildungsmethode unterteilen die Boxen das gesamte Teilgebiet Q R n : fu¨r die gew¨ahlte Verfeinerung mu¨ssten 2 14 = 16384 Boxen untersucht werden, um den Attraktor zu bestimmen. Der mengenorientierte Unterteilungsalgorithmus bestimmt in diesem Falle nur 2925 Boxen-- n¨amlich nur solche, die zum Attraktor geh¨oren. Die Anzahl der Boxabbildungen liegt jedoch etwas h¨oher, da die Anzahl 2925 durch 14 sukzessive Verdoppelungen und Reduktionen der Boxanzahl erreicht wird. Bild 2.1: Ein relativer globaler Attraktor A Q , der nach 14 Unterteilungsschritten durch U¨ berdeckung mit 2925 Boxen angen¨ahert wurde. Die Bestimmung der Bilder der Boxen, wie es in Schritt 2 des Algorithmus vorgeschrieben ist, wird durch Abbilden einer großen Zahl von sog. Testpunkten vorgenommen. Die Testpunkte k¨onnen dazu in unterschiedlichen Weisen u¨ber jede zu traktierende Box verteilt werden; das Programm GAIO 4 sieht hierfu¨r vier Optionen vor. Die Testpunkte k¨onnen die Boxen ausfu¨llen, in dem sie auf ein regelm¨aßiges Gitter verteilt werden, wobei die R¨ander miteingeschlossen werden( Grid ) oder nur das Innere mit Testpunkten bedeckt wird( InnerGrid ). Alternativ kann es vorteilhaft sein, die Punkte zufallsverteilt zu platzieren( MonteCarlo ), oder es genu¨gt in gewissen F¨allen, nur die R¨ander der Boxen mit Testpunkten zu versehen( Edges ). Die letztgenannte M¨oglichkeit, die Rechenzeit einspart, ist sinnvoll, wenn die Abbildung in eine Dimension sehr stark kontrahierend wirkt, so dass rechteckige Boxen zu sehr schlanken Bildern transformiert werden. 4 G lobal A nalysis of I nvariant O bjects steht fu¨r ein Softwarepaket, welches s¨amtliche Algoritmen der mengenorientierten Methoden vereint und verfu¨gbar macht. Bild 2.1 wurde so mit dem oben eingefu¨hrten Unterteilungsalgorithmus erzeugt, der in der GAIO Komponente RGA implementiert ist. Details und Anwendungsbeispiele zu der Software sind ihrem Handbuch[Dellnitz, Froyland, Junge 2001] zu entnehmen. 12

2.3. AUFENTHALTSWAHRSCHEINLICHKEITEN 2.3 Aufenthaltswahrscheinlichkeiten Mit dem Box-Unterteilungsalgorithmus des vorigen Unterkapitels k¨onnen relative globale Attraktoren angen¨ahert werden. Kenntnis u¨ber die Lage und Form des Attraktors innerhalb des Zustandsraumes lassen jedoch keine Ru¨ckschlu¨sse auf die Bewegung des Systems innerhalb des Attraktors zu. Oft eru¨brigt sich dieser weitere Wissensbedarf: wenn der Attraktor n¨amlich einen periodischen Grenzzyklus(engl. limit cycle, kann bei zeitkontinuierlichen Systemen auftreten), einen Fixpunkt oder einen 1-periodischen bzw. k k periodische Punkte darstellt, so bleiben u¨ber die Dynamik des Systems, sobald sich sein Zustand auf dem Attraktor befindet, keine weiteren Fragen offen. Die geometrische Dimension solcher Attraktoren ist ganzzahlig(Fixpunkt: dim= 0; geschlossene Zykluskurve: dim= 1). Wenn dagegen der Attraktor wie in Bild 2.1 keine vollfl¨achige Form besitzt, so ist seine Dimension fraktal(1 < dim < 2), also eine gebrochene Zahl. In solchen F¨allen spricht man von einem seltsamen Attraktor. Deterministisches, also vorherbestimmtes Chaos beherrscht das zeitliche Verhalten der Systemzust¨ande. Hier ist es sinnvoll, die irregul¨aren Bewegungen des Systems tiefgehender zu studieren. Hilfreich ist die Berechnung von sog. invarianten Maßen, wie etwa die Aufenthaltswahrscheinlichkeit eines ist. Sie ist ein in der zeitlichen Entwicklung unver¨anderliches Maß, das statistische Aussagen u¨ber die Wahrscheinlichkeit trifft, mit der die im Attraktor m¨oglichen Zust¨ande vom System erreicht werden. Die folgenden Abschnitte erl¨autern die Berechnung von Aufenthaltswahrscheinlichkeitswerten fu¨r alle Boxen, die einen Attraktor approximieren, wie er z. B. in Bild 2.1 gezeigt ist. U¨ bergangswahrscheinlichkeitsmatrix Die Berechnung der Aufenthaltswahrscheinlichkeiten basiert auf der Matrix der U¨ bergangswahrscheinlichkeiten(engl.: transition probability matrix). Sie wird mit P bezeichnet und kann fu¨r einen gegebenen relativen globalen Attraktor A Q berechnet werden. Sie enth¨alt die Elemente P =[ p kl ] mit k, l = 1 ,..., N . N stellt die Anzahl der Boxen dar, die den relativen globalen Attraktor A Q bilden. Jede Komponente p kl der Matrix P steht fu¨r die U¨ bergangswahrscheinlichkeit, mit der das System nach einem Schritt einen innerhalb Box B l liegenden Zustand hin zu einem in Box B k liegenden Zustand ¨andert. Die U¨ bergangswahrscheinlichkeit p kl wird angen¨ahert, indem man viele Testpunkte von innerhalb Box B l abbildet und den Anteil an Punkten z¨ahlt, die in Box B k , k = 1 ,..., N landen. Die Illustration in Bild 2.2 verdeutlicht das Vorgehen der Berechnung aller Eintr¨age p kl der U¨ bergangswahrscheinlichkeitsmatrix P . Die siebte Spalte p k, 7 erh¨alt man, indem man(hier aus Demonstrationszwecken nur) vier Testpunkte x i  B 7 , die in Box B 7 liegen, abbildet und danach fu¨r jede Box die relative Zahl an Bildpunkten z¨ahlt. 13

KAPITEL 2. MENGENORIENTIERTE NUMERISCHE METHODEN 12 456 7 89 10 11 12 13 14 15 16 3 Bild 2.2: Fiktiver Attraktor zur Demonstration der Berechnung von U¨ bergangswahrscheinlichkeiten 0 P = ··· ··· ··· ... 0 1 / 4 1 / 4 0 ... 2 / 4 0 ··· ··· ··· p 6 ,l p 7 ,l p 14 ,l ... p k, 7 Vektor der Aufenthaltswahrscheinlichkeiten Im vorigen Abschnitt wurden die U¨ bergangswahrscheinlichkeiten des Systemzustands von jeder Box in jede Box definiert und berechnet. In diesem Abschnitt wird nun der Begriff der Aufenthaltswahrscheinlichkeit y k zur Box B k definiert: Definition Aufenthaltswahrscheinlichkeit: Gegeben sei eine Menge von N Boxen B k , k = 1 ,..., N , wobei N die Anzahl der Boxen darstellt, die den relativen globalen Attraktor A Q beschreiben. Man nehme an, dass sich der Systemzustand x auf dem Attraktor A Q befinde. Die Aufenthaltswahrscheinlichkeit y k der Box B k ist definiert als die Wahrscheinlichkeit, mit der der Systemzustand x innerhalb der Box B k liegt. Es folgt noch die Definition Aufenthaltswahrscheinlichkeitsvektor: Der Aufenthaltswahrscheinlichkeitsvektor ist definiert als eine Spaltenmatrix aus den Aufenthaltswahrscheinlichkeiten y k aller Boxen B k : y :=[ y 1 , y 2 ,..., y N ] T (2.7) Unter Verwendung der beiden Definitionen kann schließlich die Berechnung der Aufenthaltswahrscheinlichkeiten hergeleitet werden. Nimm dazu an, der Aufenthaltswahrscheinlichkeitsvektor y =[ y 1 , y 2 ,..., y N ] T w¨are bereits bekannt. y l bezeichnet die Wahrscheinlichkeit, mit der ein Systemzustand x innerhalb der Box B l liegt. Die Wahrscheinlichkeit, dass dieser Zustand x vorlag und im n¨achsten Schritt in Box B k hinein abgebildet wird, ist durch das Produkt y ¯ kl = p kl y l der U¨ bergangswahrscheinlichkeit p kl und der Aufenthaltswahrscheinlichkeit y l gegeben. Die Summe u¨ber alle Boxen B l ergibt die gesamte Aufenthaltswahrscheinlichkeit B k : N y k = y ¯ kl = p kl y l . l=1 l (2.8) 14

2.3. AUFENTHALTSWAHRSCHEINLICHKEITEN In Matrix-Notation kann Gl.(2.8) als y = P y (2.9) geschrieben werden. Verglichen mit dem Standard-Eigenwertproblem y = P y stellt Gl.(2.9) das spezielle Eigenwertproblem dar, n¨amlich eines mit Eigenwert = 1. Der Aufenthaltswahrscheinlichkeitsvektor y ist der Rechts-Eigenvektor der U¨ bergangsmatrix P zum Eigenwert 1. Da U¨ bergangsmatrizen stochastischer Natur sind, haben sie stets mindestens einen Eigenwert 1[Norris 1997]. 15

.

Kapitel 3 Ultraschall-Stoßbohrer In Kapitel 2 wurden die mengenorientierten numerischen Methoden vorgestellt. Ihre effiziente Vorgehensweise bei der Approximation niedrigdimensionaler Attraktoren bzw. bei der U¨ berdeckung von Mannigfaltigkeiten wurde eindrucksvoll z. B. an der H´enon-Abbildung, an Chuas Circuit oder am Lorenz-System bewiesen[Dellnitz & Junge 2002]. Anwendungen der Methoden auf weitere bzw. neuere Systeme aus dem Ingenieursalltag sind bislang nur wenige erfolgt. Insbesondere wurden sie nicht an Systemen mit ausgepr¨agter nichtglatter Charakteristik erprobt. Die vorliegende Arbeit m¨ochte diese Lu¨cke schließen und konzentriert sich daher auf den Einsatz der mengenorientierten Methoden an einer fu¨r die Klasse der nichtglatten Systeme beispielhaften Anwendung. Ausgew¨ahlt wurde das neuartige Konzept eines Stoßbohrers. Es besticht durch seinen einfachen, aus wenigen Komponenten bestehenden Aufbau. Gleichzeitig vereinigt es daher die M¨oglichkeit einer analytisch u¨berschaubaren Modellierung mit der Analyse einer sehr komplexen Dynamik. Im Folgenden wird dieses Bohrsystem vorgestellt. 3.1 Entwicklung des Ultraschall-Stoßbohrers Vor einigen Jahren befasste sich die kalifornische Raumfahrtbeh¨orde NASA mit der Planung einer Mars-Mission. Die Entnahme von Gesteinsproben sollte fu¨r deren Untersuchung auf der Marsoberfl¨ache m¨oglich sein. Der Bedarf entstand, einen Bohrer zu entwickeln, der vielen besonderen Anforderungen genu¨gte. Herk¨ommliche Gestein-Schlagbohrer schieden aus, weil sie einige der wichtigsten Anforderungen nicht erfu¨llen konnten: so waren ein niedriges Gewicht und geringe elektrische Leistungsaufnahme wichtige Bedingungen. Wartungsfreiheit war eine ebenso wichtige Voraussetzung. Da der Bohrer auf einem kleinen, mobilen Roboter-Fahrzeug, dem Mars-Rover, zu befestigen war, durfte er weder ein großes Halte-Drehmoment, noch eine hohe Anpresskraft auf die Bohrstelle ben¨otigen. Entwickelt wurde daraufhin ein neuartiges Stoßbohrsystem, welches den Namen USDC bekam, was fu¨r U ltrasonic/ S onic D riller/ C orer steht; sie17

KAPITEL 3. ULTRASCHALL-STOSSBOHRER he dazu[Bar-Cohen et al. 2001b]. Der Bohrer nutzt Schallwellen sowohl im Ultraschallbereich(ab 20kHz), als auch im h¨orbaren Schallbereich(100 Hz­ 1000 Hz) und ist in der Lage, Vollquerschnittsbohrungen wie auch Kernlochbohrungen zur Probenentnahme vorzunehmen. Gute Bohrergebnisse in Laborversuchen und die erfolgreiche Anwendung des Bohrkonzepts zur planetaren Gesteinsprobenentnahme(siehe [Bar-Cohen et al. 2001b]) machen dieses Bohrprinzip fu¨r terrestrische Anwendungen interessant. Vorstellbar sind leichte, ger¨auscharme Handbohrer mit geringer Leistungsaufnahme fu¨r spr¨odhartes Material wie Beton und Fels. Anwendung k¨onnten sie im Klettersport, im Baugewerbe, in der Geologie, in der Medizintechnik oder als Verschu¨ttetenbergungsger¨at fu¨r Katastropheneins¨atze finden. Allerdings besteht fu¨r denkbare Anwendungen noch Optimierungspotenzial: Probebohrungen zeigten, dass die Bohrzeiten noch relativ hoch liegen, siehe z. B.[Bar-Cohen et al. 2001a]. 3.2 Aufbau des Ultraschall-Stoßbohrers Bild 3.1 zeigt einen am Heinz Nixdorf Institut gefertigten Labor-Prototypen des Ultraschall-Stoßbohrers, mit dem experimentelle Messungen gemacht wurden; in Bild 3.2 ist dessen schematischer Aufbau dargestellt. Deutlich wird dort ein herausragendes Merkmal, n¨amlich sein einfacher Aufbau. Außer der elektrischen Ansteuerung besteht das Bohrsystem aus nur drei Teilen: · piezoelektrischer Oszillator · freie Stoßmasse · Bohrstab Bild 3.1: Prototyp des Stoßbohrers fu¨r Laborexperimente 18

3.2. AUFBAU DES ULTRASCHALL-STOSSBOHRERS Oszillator Der piezoelektrische Oszillator, fu¨r den auch die Namen Transducer oder Aktor benutzt werden, ist ein metallischer Zylinder, in den ein Stapel von (hier vier) piezoelektrischen Keramiken eingespannt ist. Im Betrieb wird dieser Piezostapel durch Anlegen eines harmonisch oszillierenden Spannungssignals an seine Elektroden zu Schwingungen angeregt. Als Anregefrequenz wird die unterste L¨angseigenfrequenz des Oszillators gew¨ahlt. Dadurch wird die erste L¨angseigenmode in Resonanz zu Schwingungen angeregt. Die erste L¨angseigenform eines Stabes beschreibt eine stehende sog. Lambda-halbeWelle. Das heißt, an den Stabenden bilden sich Schwingungsb¨auche, dazwischen ein Schwingungsknoten. Der fu¨r Messungen verwendete Oszillator hat eine L¨ange von l = 125 mm, was bei einer Schallwellengeschwindigkeit in Stahl f = c 2 l = 20 kHz von c = 5 entspricht. km s einer ersten L¨angseigenfrequenz von etwa Hülse(Handhalterung) ~ U Piezo-Stapelaktor Oszillator Flugmasse Bohrer Gestein Bild 3.2: Schematischer Aufbau des Ultraschall-Stoßbohrsystems USDC Freie Stoßmasse Der zweite Bestandteil des Bohrers ist eine wenige Gramm schwere, freifliegende Stoßmasse. Sie kann z. B. Unterlegscheiben- oder Kugelform besitzen. Die Stoßmasse fliegt in einem kleinen Zwischenraum zwischen dem schwingenden Ende des Ultraschall-Oszillators und der dem Bohrgut(Gestein) abgewandten Seite des Bohr-Stabes hin und her. Dadurch u¨bertr¨agt die Stoßmasse Stoßenergie vom Oszillator hin zum Bohrer. Ihre besondere 19

KAPITEL 3. ULTRASCHALL-STOSSBOHRER Aufgabe besteht darin, eine Transformation der Ultraschall-Schwingungen von kleiner Amplitude, wie sie am Ende des Oszillators auftreten, hin zu K¨orperschallimpulsen im Bohrstab von großer Intensit¨at vorzunehmen. Als weitere Funktion erh¨alt sich die Stoßmasse den fu¨r ihre Flugbewegung notwendigen Zwischenraum. Die vorliegende Arbeit setzt einen Fokus auf die Analyse und Beschreibung der Bewegung dieser Stoßmasse, von der der Bohrprozess im Wesentlichen abh¨angt. Bohrstab Der dritte Teil des Bohrsystems ist der Bohrstab. Er leitet die Stoßwellen, die wiederholt durch Abprallen der freien Stoßmasse an seiner systeminneren Seite erzeugt werden, hin zur Bohrer-Gestein-Kontaktstelle, wo sie im zu bearbeitenden Gestein dessen Bruchspannung u¨berschreiten und so den eigentlichen Bohrfortschritt erzielen. Obgleich die oben vorgestellte Bohrtechnologie vielversprechend ist, konnte ihr Verhalten auf dynamischer Ebene bislang kaum entschlu¨sselt und durchdringend verstanden werden. Beispielsweise ist es noch unklar, inwiefern die Gr¨oße, Gestalt und Materialzusammensetzung der Stoßmasse, die Anregefrequenz und die Bohrstabgeometrie den Materialabtrag beeinflussen. Ein hinreichendes Systemverst¨andnis ist allerdings unabdingbar, m¨ochte man das System-Design fu¨r weitere Anwendungsfelder nutzbar machen und Eigenschaften wie die Bohrgeschwindigkeit optimieren oder nur auf gr¨oßere Bohrquerschnitte und h¨ohere Leistungsaufnahmen skalieren. 20

Kapitel 4 Modellierung und Analyse nichtglatter Dynamik-Stand der Technik Im Fokus der vorliegenden Arbeit steht die Analyse von nichtglatten Systemen der Dynamik. Diese Systeme bilden eine besondere Untergruppe der nichtlinearen Systeme. Seit einigen Jahrzehnten kommt Nichtlinearit¨aten in der Wissenschaft immer gr¨oßere Aufmerksamkeit zu. Dementsprechend ist die Literatur, in der man sich mit nichtlinearen Ph¨anomenen auseinandersetzt, in dieser Zeit auf eine beachtliche Zahl angewachsen. In den Anf¨angen der Arbeit mit nichtlinearen Systemen entwickelten Mathematiker und Ingenieure mathematische Modelle von technischen Systemen, die von sehr grundlegender Natur waren. Zu den bekanntesten Modellen z¨ahlen Differentialgleichungen von Duffing 1 , van der Pol 2 und Lorenz 3 . Vor dem Zeitalter der maschinellen, numerischen Rechenunterstu¨tzung waren die Untersuchungsmethoden rein analytisch bzw. experimentell. Die zur theoretischen Analyse damals entwickelten wichtigsten Methoden sind 1 Der deutsche Ingenieur und Experimentalist Georg Duffing (1861­1944) formulierte 1918 die Duffing-Gleichung, indem er dem Potential eines linearen Oszillators eine kubische Ru¨ckstellkraft hinzufu¨gte. Das dynamische System entspricht einem einfachen, ged¨ampften Feder-Masse-Schwinger mit nichtlinearer Federkennlinie [Abraham & Shaw 2000]. 2 Der niederl¨andische Physiker und Philips-Ingenieur Balthasar van der Pol (1889­ 1959) untersuchte 1920 die van der Pol-Gleichung, die einen Feder-Masse-Schwinger mit nichtlinearer D¨ampfung modelliert[Van der Pol 1920]. 3 Der emeritierte Meteorologieprofessor des MIT Edward Norton Lorenz (geb. 1917 in den USA) stellte ca. 1960 sein Wettermodell auf, aus dessen starker Abh¨angigkeit von Anfangsbedingungen der sogenannte "Schmetterlingseffekt" hervorging[Lorenz 1963]. 21

KAPITEL 4. MODELLIERUNG UND ANALYSE NICHTGLATTER DYNAMIK-- STAND DER TECHNIK unter anderem den Arbeiten von Lagrange 4 , Poincar´e 5 , Ljapunov 6 und Birkhoff 7 zu verdanken, die heute immer noch Aktualit¨at genießen. Einige Methoden, die auf diese Pioniere zuru¨ckgehen, werden im folgenden Unterkapitel beschrieben. Im Unterkapitel 4.2 wird auf ju¨ngere Methoden eingegangen, die erst in den letzten Dekaden mit dem Aufkommen maschineller Rechenunterstu¨tzung entwickelt wurden. In dieser Arbeit werden die genannten Methoden auf zwei Schwingstoßsysteme angewendet. Auf Modellierungsans¨atze zu ¨ahnlichen Systemen wird in Unterkapitel 4.3 verwiesen. 4.1 Klassische Methoden zur Analyse nichtlinearer dynamischer Systeme Dieses Unterkapitel stellt die g¨angigsten Methoden zusammenfassend dar, die fu¨r Untersuchungen an nichtlinearen dynamischen Systemen in der Literatur regelm¨aßig Anwendung finden. Die nichtglatten dynamischen Systeme geh¨oren auch in diese Gruppe. Daher bilden die nachfolgenden Techniken auch fu¨r alle Schwingstoßsysteme eine wesentliche Grundlage. 4.1.1 Zeitreihenanalyse Die naheliegendste M¨oglichkeit, Einblick in die Bewegungsabl¨aufe eines Systems zu erhalten, ist die Zeitreihenanalyse. Dabei wird die zeitliche Entwicklung von Zustandsvariablen des Systems in einem Diagramm u¨ber die Zeit aufgetragen. Ein Beispiel dafu¨r zeigt Bild 6.11 auf Seite 66. Die zugrundeliegenden Daten k¨onnen experimentell oder durch Simulation eines mathematischen Modells gewonnen werden. An dieser Darstellungsweise lassen sich Eigenschaften wie Periodizit¨at, asymptotisches oder stabiles Verhalten, Stabilit¨at oder Chaos in der Regel nicht ablesen. Ihr Vorteil liegt in der einfachen Zug¨anglichkeit der Bewegungsabl¨aufe auch fu¨r Nichtwissenschaftler. 4.1.2 Phasenportrait Wgeenradnennt)dbieenZuetizttaulsndLaduiefpZaursatmanetdesrwiemrteZduesstaSnydsstreamusmin(achurcohn"oPlohgaisscehnerrauFmol"4 Der italienische Schu¨ler Leonhard Eulers(1707­1783) Joseph-Louis Lagrange (1736­ 1813) wurde als Begru¨nder der analytischen Mechanik bekannt[Wikipedia 2007]. Fu¨r die Stabilit¨atsuntersuchungen nichtlinearer diskreter und kontinuierlicher Systeme lieferte er einen Stabilit¨atsbegriff zur Klassifikation von L¨osungsverhalten [Thompson & Stewart 1986], der aber heute eher ungebr¨auchlich ist. 5 Nach dem franz¨osischen Mathematiker, Physiker und Philosophen Jules Henri Poincar´e (1854­1912) ist die Poincar´e-Abbildung benannt, die den Zustand eines Systems stroboskopartig zu diskreten, periodischen Zeitpunkten beleuchtet und dadurch besondere Strukturen innerhalb chaotischer Bewegungen enthu¨llt. 6 Aleksandr Mikhailovich Liapounov (1857­1918) wurde beru¨hmt fu¨r seine grundlegenden Ideen zum Begriff der Stabilit¨at von L¨osungen. 7 Der US-amerikanische Mathematiker George David Birkhoff (1884­1944) wurde 1913 durch den Beweis Poincar´es letzten Theorems beru¨hmt und arbeitete im Bereich der statistischen Mechanik. 22

4.1. KLASSISCHE METHODEN ZUR ANALYSE NICHTLINEARER DYNAMISCHER SYSTEME ge eingetragen, so entsteht eine Trajektorie(auch Orbit genannt). Fu¨r eine Darstellung n dimensionaler Trajektorien, wenn n> 3, muss die Trajektorie in den Raum hinein projiziert werden, der von zwei oder drei Zustandsvariablen aufgespannt wird. Diese Betrachtungsweise erm¨oglicht das Erkennen von sich wiederholenden Zust¨anden und periodischen Bewegungen. Bei zeitdiskreten Systemen, die durch Abbilden eines Zustands auf einen folgenden Zustand beschrieben werden, besteht eine Trajektorie aus einzelnen Zustandspunkten, die nicht miteinander verbunden sind. Beispiele fu¨r diskrete Orbits zeigt Bild 6.14 auf Seite 70. Periodische Bewegungen lassen sich erkennen, wenn Zustandspunkte wiederholt auf sich selbst oder einen ihrer Vorg¨anger abgebildet werden. Periodizit¨at in zeitkontinuierlichen Systemen ist durch sich schließende Orbits charakterisiert. Aus Phasenportraits lassen sich weiterhin auch Aussagen u¨ber das Langzeitverhalten eines Systems treffen. Z. B. k¨onnen die Zustandsbereiche, innerhalb derer sich ein durch D¨ampfung begrenztes System aufh¨alt, abgesch¨atzt werden. Die Daten, aus denen sich Phasenportraits erstellen lassen, sind die der Zeitreihe, k¨onnen also experimentell aufgezeichnet oder simulatorisch berechnet worden sein. 4.1.3 Periodische Punkte und Fixpunkte Fixpunkte im Zustandsraum von zeitdiskreten Systemen stellen stabile Gleichgewichtslagen dar. Fu¨r zeitdiskret dargestellte Systeme repr¨asentieren Fixpunkte sich periodisch wiederholendes Verhalten. Besteht ein Orbit aus m Fixpunkten, so odische Fixpunkte wlaisrsdendiseicBhenwuergusenlgte"n m a napleyrtiiosdchiscbhe"stgimenmanenn.t.InMPehhrapseerni-portraits sind sie aber einfach erkennbar. 4.1.4 Ljapunov-Stabilit¨at und Ljapunov-Exponent Eine weitere Charakterisierung des Systemverhaltens kann durch die Beschreibung der Stabilit¨at von L¨osungen vorgenommen werden. Es gibt mehrere Stabilit¨atsbegriffe. Eine L¨osung u ( t ) eines autonomen oder nichtautonomen Differentialgleichungssystems wird Ljapunov-stabil genannt, wenn fu¨r jede beliebige, kleine Zahl  > 0 eine Zahl = ( ) > 0 existiert, so dass jede andere L¨osung v ( t ), fu¨r die zur Zeit t = t 0 die Abstandsbeziehung u ( t 0 ) - v ( t 0 ) < gilt, fu¨r alle Zeiten t> t 0 die Beziehung u ( t ) - v ( t ) < erfu¨llt[Nayfeh 1995]. Anders ausgedru¨ckt wird eine L¨osung dann Ljapunov-stabil genannt, wenn sie sich ab einem bestimmten Zeitpunkt nie mehr aus den gedachten Schl¨auchen begrenzter Durchmesser, die um alle benachbarten L¨osungen herum verlaufen, herausbewegt. Ljapunov-Exponenten hingegen geben ein Maß fu¨r die Stabilit¨at einer L¨osung an. Dabei vergleicht man zwei Trajektorien, die durch im Zustandsraum sehr dicht benachbarte Zust¨ande laufen. In der weiteren zeitlichen Entwicklung verh¨alt sich der Abstand zwischen den Zust¨anden beider Trajektorien im Zustandsraum meistens exponentiell: der Ljapunov-Exponent gibt die St¨arke dieser exponentiellen Separation oder Ann¨aherung an. Entscheidend ist meistens nur das Vorzeichen des Ljapunov-Exponenten. Es gibt 23

KAPITEL 4. MODELLIERUNG UND ANALYSE NICHTGLATTER DYNAMIK-- STAND DER TECHNIK an, ob es sich um eine stabile oder um eine instabile Bewegung z. B. in Bezug auf einen Fixpunkt handelt. 4.1.5 Poincar´e-Abbildung Die Poincar´e-Abbildung ist eine besondere Darstellungsform der Dynamik, die die Zust¨ande von zeitkontinuierlichen Systemen nur zu diskreten Zeitpunkten stroboskopartig beleuchtet. Dieses Vorgehen erzeugt aus einem kontinuierlichen System eine diskrete Abbildungsvorschrift. Geschlossene (mehr-)periodische Trajektorien werden zu(mehreren) Fixpunkten und lassen sich einfach identifizieren. Besondere Relevanz hat die Poincar´eAbbildung bei chaotischer Bewegung: irregul¨ar verlaufende Trajektorien bilden im Phasenraum ein Liniengewirr. Blitzt man den mit der Zeit durch den Phasenraum wandernden Zustandspunkt dagegen in periodischen Abst¨anden gedanklich mit einem Stroboskop und friert die so entstehenden, einzelnen Punkte ein, so enthu¨llt eine deterministisch-chaotische Bewegung m(siietuhnetaeur cehinUenetresrtkaaupniltieclh2e.3S)t.ruAkltsuPr,erdiioedmenadnau"esreletsiganmeetnsiAchttbraeki tnoirc"htnaeuntnotnomen Systemen(also solchen mit expliziter Zeitabh¨angigkeit, etwa durch eine ¨außere Kraft- oder Weganregung) fu¨r die Abbildung die Anregeperiodendauer. Bei autonomen Systemen ist die Wahl der Periodendauer weniger klar. Es eignen sich meistens systemeigene Schwingungsdauern. 4.1.6 Verzweigungs- oder Bifurkationsdiagramm Mit einem Bifurkations- oder Verzweigungsdiagramm 8 l¨asst sich gut die Abh¨angigkeit der L¨osungstypen von einem Systemparameter untersuchen. Dazu wird der Wert einer Zustandsgr¨oße an bestimmten Zeitpunkten in Abh¨angigkeit des Systemparameters dargestellt. Bei zeitkontinuierlich definierten Systemen wertet man wie bei der Poincar´e-Abbildung den Systemzustand an periodischen, diskreten Zeitpunkten aus. So lassen sich im Bifurkationsdiagramm periodische L¨osungen mit ihrer Periodizit¨atszahl und chaotische L¨osungen gut voneinander in Abh¨angigkeit des Systemparameters unterscheiden. Bild 4.1 zeigt die fu¨r nichtlineare und Stoßsysteme typische Kaskade von Periodenverdoppelungen in Bifurkationsdiagrammen fu¨r beide Zustandsgr¨oßen V und des dynamischen Systems 2. Ordnung, das in Unterkapitel 5.1 hergeleitet wird. Mit gr¨oßer werdendem Systemparameter 1 leiten die Periodenverdoppelungen ins Chaos u¨ber, das an den rechten R¨andern in Bild 4.1 in Form gestreuter Punkte erkennbar ist. 4.1.7 Einzugsgebiet In besonderen F¨allen, die vom Zusammenwirken der Systemparameter abh¨angen, k¨onnen unterschiedliche, station¨are L¨osungen(das sind die Be8 Bifurkation : Gabelung; bifurkare(lateinisch), bifurcer(franz¨osisch), bifurcate(englisch): sich(in zwei A¨ ste) gabeln. 24

4.2. ZELLABBILDUNGSMETHODEN V [m/s] [ rad] 3.9496 3.9495 3.9494 3.9493 3.9492 3.9491 3.949 1.501 1.5 1.499 1.498 1.497 1.496 1.495 1.494 1.493 1.492 1.491 0.612100 0.612105 1 0.612110 0.612100 0.612105 1 0.612110 Bild 4.1: Periodenverdoppelungen im Bifurkationsdiagramm des zeitdiskreten Systems 2. Ordnung(Unterkapitel 5.1) fu¨r den Systemparameter 1 . Die Verzweigungsdiagramme werden beispielhaft fu¨r beide Zustandsgr¨oßen V und gezeigt. wegungen, die sich nach einiger Zeit einstellen und danach ihre Eigenschaften fu¨r alle Zeiten beibehalten) gleichzeitig in einem System existieren. Die Menge an Zustandspunkten, die von allen Trajektorien eingenommen werden, nachdem Anfangsbewegungen gungen") abgeklungen sind, werden (A"ttrtaranksiteonrteen" goednearnn"tin(sstiaethieonU¨anreteBrkeawpei-tel 2.1). Eine L¨osungstrajektorie wird nach einiger Zeit immer in einen der Attraktoren laufen. Von ne vom Anfangszustand welchem Attraktor abh¨angig, bei dem siihere"aBngereezcohgneun"ngwbiredg,ainsnt .aDlleaisEinzugsgebiet(engl. basin of attraction oder domain of attraction) eines ZFuixsptaunndkstreasumma, rvkoienrtdednieejnenaiugesnTAranjefakntogrsibeendiinngeuinngeennbaeustfimeimnetren"KAatrttrea"ktiomr hinein laufen. Bild 7.7 auf Seite 96 zeigt die beiden Einzugsgebiete zu zwei Attraktoren(weiß und schwarz). Gemeinsam fu¨llen sie den Zustandsraum komplett aus. 4.2 Zellabbildungsmethoden Seit etwa dem Ende des 20. Jahrhunderts stehen der Wissenschaft H¨ochstleistungsrechenzentren zur Verfu¨gung, deren gewaltige Rechenleistung Jahr fu¨r Jahr weiter nach oben schnellt. So entstehen auf dem Gebiet der Analyse nichtlinearer Systeme und der Chaosforschung M¨oglichkeiten zur Entwicklung immer neuer, numerischer Verfahren: Hsu aus Berkeley dis25

KAPITEL 4. MODELLIERUNG UND ANALYSE NICHTGLATTER DYNAMIK-- STAND DER TECHNIK kretisiert den Zustandsraum des Systems und entwickelt die Zellabbildungsmethoden, die in Fachkreisen unter dem Namen bekannt geworden sind. Publiziert wurden sie in "[HCseull-1t9o8-C0]e,ll[HMsuap1p9i8n1g]", [Hsu 1987] und[Hsu 1992]. Ein System kann damit global, also in seinem gesamten(aber technisch bedingt fu¨r die Berechnungen dennoch begrenzten) Zustandsraum betrachtet werden. Eingesetzt und erweitert wurden die Zellabbildungsmethoden unter anderen von Tongue in[Tongue 1987] und [White & Tongue 1994]. In diesem Unterkapitel soll die Funktionsweise der Zellabbildungsmethoden erl¨autert werden. Denn dieser Technik in der Grundidee ¨ahnlich, aber um entscheidende Verfeinerungen erg¨anzt, sind die recht jungen orientierten numerischen Methoden", die seit etwa einer Dekade e"nmtwenicgkeenltwerden. Detaillierte Informationen dazu sind Kapitel 2 zu entnehmen. Funktionsweise Nach der Idee der Zellabbildungsmethode betrachtet man den Phasenoder Zustandsraum nicht wie herk¨ommlich als Kontinuum, sondern als Ansammlung sehr vieler Zustandszellen. Jede Zelle soll einen Zustandsbereich repr¨asentieren. Hsu unterscheidet dabei zwei Arten, n¨amlich die Einfache Z Z ("u e G l s l t a ean b n b ed i r l as d bl u ize n er g edi(c"hCSeielmlinpMelreaCpjepedlielnngM"aZ).peplIlienngnd"ue)rruaneuidnffdaeiceihn G een eeni Z enr e za l i l lg a ies b i b e Z r il t ed e lule Z n e ga ll b a wg b ei b rb i d l i d ld u de n et g r. Diese Methode dient prim¨ar zur Darstellung periodischer Systemantworten. Falls die Abbildungsvorschrift in einer Dimension expandierend ist, so erstreckt sich der Bildbereich einer Zelle u¨ber mehrere Bildzellen. Weil die einfache Zellabbildung, die nur eine Bildzelle beru¨cksichtigt, so zu Verf¨alschungen der tats¨achlichen Systemantwort fu¨hren kann, entstand die Erweiterung der generalisierten Zellabbildung. Danach kann jede Zelle eine endliche Zahl von Bildzellen aufweisen. Die Bildzellen werden bestimmt, indem eine endliche Zahl von Anfangsbedingungen, die innerhalb der Ursprungszelle liegen, einmal integriert(d. h. abgebildet) werden. Bildzellen sind l¨asst jene, sich deieindeurU¨chbedrigeasenngsSwcahhrirtstch"geientlriocffhekne"itswvuerrtdeeinlu.nAguskodniesstermuieVreonrg. ehDeine U¨ bergangswahrscheinlichkeiten fu¨r alle N Zellen k¨onnen in einer N × N großen Matrix abgelegt werden, die danach die gesamte Information u¨ber die Dynamik des Systems enth¨alt. Ein Bild des globalen Verhaltens nichtlinearer Systeme entsteht. Seltsame Attraktoren, periodische Orbits und statistisch-deterministische Eigenschaften lassen sich lokalisieren. Aufbauend auf dem Konzept der Zellabbildung haben sich mittlerweile in verschiedene Richtungen Erweiterungen oder A¨ nderungen aufgetan. Zu nennen ist die Interpolierende Zellabbildungsmethode und die Zellabbildungsmethode fu¨r diskontinuierliche Systeme. Auf beide Erweiterungen wird in den n¨achsten Abschnitten eingegangen. 26

4.3. MODELLIERUNG VON SCHWINGSTOSSSYSTEMEN Interpolierende Zellabbildungsmethode Die Zellabbildungs-Methoden betrachten ein Kontinuum-- n¨amlich den Zustandsraum-- notwendigerweise als Diskontinuum. Um wegen der begrenzten Feinheit der Diskretisierung hieraus m¨oglicherweise resultierende Verf¨alschungen auszuschließen, schl¨agt Tongue die Interpolierende Zellabbild B u il n d g pu(n"IkntteerdpeorlaintteedgrCieerltleMn aApnpfianngg"s)bevdoirng[Tuonnggeunem1i9t8i7h].reHnie e r x i a n kt w en erdWenertdeine im Kontinuum gespeichert. Durch eine bi-lineare Interpolationsstrategie gelingt es, die Trajektorie zu jeder Anfangsbedingung zu ermitteln, solange sie im Feld der Zellen bleibt. Somit gewinnen die L¨osungen wieder ihren urspru¨nglich kontinuierlichen Charakter zuru¨ck. Zellabbildungsmethode fu¨r diskontinuierliche Systeme In[Kraker, van der Spek und van Campen 1999] werden einige Erweiterungen der Zellabbildungsmethoden fu¨r diskontinuierliche Systeme vorgeschlagen. Zus¨atzlich beschreiben die Autoren eine Parametervariationsmethode fu¨r Zellabbildungen und eine Erweiterung fu¨r Systeme mit vielen Freiheitsgraden. Die Erweiterung der Parametervariations-Zellabbildung baut auf der einfachen Zellabbildung auf. Mit ihr ist es m¨oglich, die Evolution-- also die Verschiebung und Verformung-- der Einzugsgebietsgrenzen zu bestimmen, wenn sich Systemparameter ver¨andern. Ein spezieller Algorithmus berechnet hierfu¨r nur die A¨nderung der Grenzen, anstelle sie fu¨r jeden Systemparameterwert neu zu berechnen. Dies fu¨hrt zu Einsparungen an Rechenoperationen. Die Zellabbildungsmethode fu¨r Systeme mit vielen Freiheitsgraden wurde durch den exponentiellen Anstieg der Zellenanzahl mit dem Zunehmen der Systemfreiheitsgrade motiviert. Dies erfordert ein hohes Speichervolumen. Zur Zeit der Entwicklung dieser Erweiterung war es kaum m¨oglich, Systeme mit mehr als zwei Freiheitsgraden mit Zellabbildungsmethoden zu untersuchen. Die Darstellung bzw. Auswertung von Ergebnissen aus Zustandsr¨aumen mit mehr als zwei Dimensionen auf Papier oder einem Monitor erfordert stets die Projektion auf einen zu w¨ahlenden 2-dimensionalen Unterraum des Zustandsraumes(also eine Ebene). Eine große Datenmenge bleibt ungenutzt. Die Idee der Erweiterung ist, die Wahl der Projektionsebene vorab zu treffen und nur fu¨r solche Anfangsbedingungen das Langzeitverhalten zu berechnen, die in der betrachteten Ebene liegen. Aufgabe ist es folglich, die Schnittmenge zwischen und dem Attraktor zu finden. 4.3 Modellierung von Schwingstoßsystemen In unz¨ahligen Bereichen-- nicht nur in den Ingenieurswissenschaften-wechseln Systeme ihre Bewegungsgleichung, sobald eine Zustandsgr¨oße-etwa eine Verschiebungs- oder Geschwindigkeitsgr¨oße-- einen Schwellen27

KAPITEL 4. MODELLIERUNG UND ANALYSE NICHTGLATTER DYNAMIK-- STAND DER TECHNIK wert u¨berschreitet. In vielen Anwendungsgebieten des Maschinenbaus treten Unstetigkeiten oder St¨oße auf. In der Elektrotechnik wechseln beispielsweise Dioden ihr Verhalten bei einer Durchlassspannung. Dies fu¨hrt dazu, dass sich die ¨außeren Einflu¨sse, die das Systemverhalten bestimmen, ¨andern. Wegen ihrer Relevanz und dem h¨aufigen Vorkommen entsprechender Anwendungen gewinnen nichtglatte dynamische Systeme in der Vergangenheit und heute große Aufmerksamkeit. Dieser Trend macht sich auch in der Literatur bemerkbar. mische Systeme" in [Leine & Nijmeijer drei Kategorien: 2004] unterteilen "nichtglatte dyna1. Systeme, die eine Ru¨ckstellkraft in ihrer Bewegungsdifferentialgleichung enthalten, die eine vom(z. B.) Verschiebungszustand abh¨angig nichtglatte, aber kontinuierliche Kennlinie besitzt. Beispiel: FederMasse-Schwinger, der mit einer zweiten Feder erst ab einer gewissen Verschiebung in Kontakt steht. Die zweite Federkraft ist eine kontinuierliche aber, nichtglatte Funktion der Verschiebung. 2. Systeme, die einen Kraftterm enthalten, der nichtkontinuierlich vom (z. B.) Geschwindigkeitszustand abh¨angt. Beispiel: Bei trockener Reibung oder Visko-Elastizit¨at h¨angt die Wirkrichtung der Reibungskraft vom Vorzeichen der Geschwindigkeit ab. Beim U¨ berschreiten der NullGeschwindigkeit springt die D¨ampferkraft. 3. Systeme, bei denen eine Zustandsgr¨oße nichtkontinuierliche Spru¨nge ausu¨ben kann. Beispiel: Bei einem Stoß zwischen zwei K¨orpern ¨andert der Geschwindigkeitszustand sprunghaft seinen Wert. Bei den in dieser Arbeit untersuchten Stoßsystemen handelt es sich um Systeme dieser dritten Art. Eine interessante Auswahl von Arbeiten aus dem Gebiet der nichtglatten Dynamik sind[Babitsky 1998],[Brogliato 1996], [Leine et al. 2000],[Leine& Nijmeijer 2004],[Neumann et al. 2007] und [Pavlovskaia & Wiercigroch 2003]. Babitsky baut in seinem Buch die Theorien u¨ber Schwingstoßsysteme aus und nutzt Methoden, die auf dem Konzept der sog. harmonischen, ¨aquivalenten Linearisierung basieren. Dies ist bei periodischem Systemverhalten m¨oglich. Schwingstoßsysteme definiert er als solche, bei denen St¨oße systematisch mit der Frequenz des Schwingers auftreten. Im Rahmen dieser Arbeit wird der Begriff des Schwingstoßsystems auf Systeme aufgeweitet, in denen eine schwingende Komponente in die Stoßdynamik involviert oder fu¨r diese verantwortlich ist. Leine entwickelt Pfadverfolgungstechniken zur Bifurkationsanalyse von periodischen L¨osungen in nichtglatten mechanischen Systemen weiter, in denen trockene Reibung eine Unstetigkeit in der Kraft erzeugt(siehe 2. in der Unterteilung oben). nes"E)i,nbeeAi drterSetionßebhoahrrmeronfu¨isrchErkdrbafotharnugnegreegnte(eMngals.se"gerrosutnndacmh oU¨libnegrwminacdheineines leeren Zwischenraumes Kontakt mit einem Feder-D¨ampfer-Element hat, untersucht Wiercigroch. Auch Schiffsvert¨auungen an Hochseebojen fu¨r 28

4.3. MODELLIERUNG VON SCHWINGSTOSSSYSTEMEN O¨ l- und Gastanker fallen in diese Kategorie. Die turmf¨ormigen Bojen sind mit ihrem Fuß am Meeresgrund verankert und oszillieren, harmonisch angeregt von vorbeilaufenden Meereswellen, durch die Auftriebskraft wie ein kopfstehendes Pendel. Ein daran vert¨autes Tankschiff verh¨alt sich im Vergleich zur Pendelschwingbewegung wie ein festes Objekt. Mit den Pendelbewegungen in die eine Richtung spannt sich das Tau und h¨angt in die andere Richtung schlaff durch. Die Ru¨ckstellkraft, die der Bojenturm durch das Tau erf¨ahrt, ist nicht kontinuierlich und somit nichtglatt. Diese nichtglatte Steifigkeitskennlinie kann in diesem Fall zu gef¨ahrlichen, unerwarteten, subharmonischen Resonanzen fu¨hren[Thompson & Stewart 1986]. Ein großes Anwendungsgebiet fu¨r aus St¨oßen resultierende Dynamik ist Getrieberattern und-klacken, das aufgrund der Zahneingriffst¨oße, von Teilungsfehlern, im lastfreien Betrieb oder bei Lastwechseln kaum vermeidbar ist. Nichtglatte Modellierung ist in allen weiteren Stoßbohruntersuchungen[Wiercigroch, Wojewoda, Krivtsov 2000],[Neumann et al. 2007], [Potthast et al. 2007a] oder auch bei Druckh¨ammern in Druckmaschinen notwendig. Eine weitere Anwendung von Ultraschallwellen in Bohrprozessen ist im DFG-Projekt Tieflochbohren zu finden, in dem h¨ohere Bohrraten und exaktere Bohrergebnisse durch ultraschallu¨berlagertes Bohren erm¨oglicht werden[Potthast et al. 2007b]. "Bouncing Ball" Systeme Im System des Stoßbohrers, der in[Bar-Cohen et al. 2001b], [Badescu et al. 2005] und Kapitel 3 vorgestellt wird, bewegt sich eine freie Stoßmasse zwischen einer oszillierenden Transducerspitze und der Ru¨ckseite eines Bohrstabes. Die Stoßmasse kann aufgrund ihrer geringen Gr¨oße und Masse als starres Element angesehen werden. Somit wird sie in Modellen als Massepunkt beschrieben. Ein vergleichbares System, das Mathematiker und Physiker schon lange besch¨aftigt, ist unter dem Namen d "b u o n u g n , c w in o g be b i al h l" ierbfeu¨kranmnitt. In der Partikel-Physik rein elastischen St¨oßen fand das System Anwen(Stoßzahl 1.0) gerechnet wird. Mittlerweile stellt das Bouncing Ball System ein Standard-Beispiel der Stoßdynamik und mathematischen Analyse nichtglatter Systeme dar. Sehr detailliert wurde es bisher untersucht, unter anderen von[Holmes 1982], [Guckenheimer & Holmes 1983],[Tufillaro, Mello, Choi, Albano 1986], [Mello, Tufillaro 1985]. Mitunter wird das Modell des Bouncing Ball Systveemresinsfoagcharen"ddeisAsinpnaatihvme eSgtaentrdoaffredn-,AdbibeilFdluunggh"¨ohgeenl¨aagnenut.mDGabr¨oeßi ewnoirrddnouftngdeine u¨ber der Schwingamplitude der Anregung[Holmes 1982]. Nur mit dieser Annahme ist es m¨oglich, die Abbildung explizit zu definieren. [Tufillaro, Mello, Choi, Albano 1986] stellen jedoch eine implizite, exakte Abbildung auf und weisen Periode-1-Fixpunkte samt deren Stabilit¨at durch Linearisieren und Vorzeichenauswerten der Eigenwerte analytisch nach. So der laAssnernegseifcrhequfu¨ernzd"ieGPreanrazmkuertveernebeznweis"cAhennregGeaebmieptleitnuduentu¨ebresrchQiedulaidchraert Periodizit¨at bestimmen, die erstaunlicherweise alle Geraden sind. 29

KAPITEL 4. MODELLIERUNG UND ANALYSE NICHTGLATTER DYNAMIK-- STAND DER TECHNIK Die A¨ hnlichkeiten des Bouncing Ball Systems mit dem UltraschallStoßbohrsystem aus Kapitel 3, das in den folgenden Kapiteln in unterschiedlichen Modellen beschrieben und untersucht wird, beschr¨anken sich auf die Modellierungsweise mit Newtons Stoßzahl und die Art der Phasenraumdarstellung.[Brogliato 1996] benennt die Abbildungsvorschriften, die als Differenzengleichung die Bewegungsgleichungen fu¨r Stoßsysteme ohne dspieeicPhoeirnncdaer´eE-Alebmbeinldtuendgarfsu¨terlledne,ra"r i t m ig p e ac S t y m st a e p m"eunddefihnaietrtd.amNietuglbeieczhu¨zgeliitcihg der nichtlinearen Systemtheorie ist, dass in der Poincar´e-Abbildung der "piemripoadcistchmearpZ"eidtipeunmk¨otgeligcehwer¨awheltiswe eurndreeng.eSlmom¨aßitigiesnt Stoßzeitpunkte anstelle die Poincar´e-Abbildung bereits durch die Differenzenbewegungsgleichung gegeben. In einem wesentlichen Punkt unterscheidet sich die Dynamik der SlettozßtemrearsseistimunUbletsrcahscrh¨analkl-tStbozßgbl.ohsreeinr esrtaFrklugvho¨onhde.erSdomesit"blo¨aussntcienignebahlolsh"e: Abfluggeschwindigkeit eine große Flugh¨ohe und somit eine lange Flugzeit folgen. Im Stoßbohrer gilt dieser Zusammenhang reziprok: die Flugh¨ohe ist beschr¨ankt; somit bedingt eine große Abfluggeschwindigkeit der Masse eine kurze Flugdauer. Vielfach finden die analytischen, numerischen oder einxpeeirniemmenutenltleernenUFnrteeqrsuuecnhzubnegreenichdsetsat"tb, oiunndceimngAbnarlelgs"epeinrioddeer Literatur und Stoßperiode gleiche Gr¨oßenordnung haben k¨onnen. Den Besonderheiten, die bei Frequenzen im Ultraschallbereich auftreten, ist noch nicht Rechnung getragen worden. So bedingen die im Vergleich zur Anregeoszillation m¨oglicherweise flachen Flugbahnen besondere Beru¨cksichtigung beim numerischen L¨osen der impliziten Bewegungsgleichungen. Auf einen speziell auf dieses Problem zugeschnittenen Algorithmus wird in Unterkapitel 5.3 eingegangen. Weiterhin ist es notwendig, sogenannte Mehrfachst¨oße (siehe dazu Fußnote auf Seite 37) gesondert zu beru¨cksichtigen: unter bestimmten Voraussetzungen kann es zwischen Oszillator und Stoßmasse zu wiederholtem Aufeinanderstoßen kommen, bevor ein vollst¨andiges L¨osen der Stoßmasse vom Oszillator m¨oglich ist. In den folgenden Kapiteln dieser Arbeit wird das oben beschriebene, dynamische System mit begrenzter Flugh¨ohe in zwei Varianten analytisch untersucht. Im ersten Modell(Unterkapitel 5.1) wird die begrenzte Flugh¨ohe auf einem festen Wert gehalten; im zweiten Modell(Unterkapitel 5.2) ist die Flugh¨ohe variabel und wird als Systemzustand aufgenommen. 30

Kapitel 5 Modellierung des Ultraschall-Stoßbohrers Notwendiger und wesentlicher Schritt hin zur analytischen oder numerischen Untersuchung eines realen Systems ist zun¨achst die mathematische Abbildung des Systems in einem Modell. Diesem ersten Schritt muss Beachtung geschenkt werden, da von der Gu¨te eines Modells nachweislich die Ergebnisse einer Simulation, die Optimierung von Systemeigenschaften oder die Systemauslegung abh¨angen. Ein Modell beschreibt stets ein Abbild der Natur. H¨aufig idealisiert man darin Eigenschaften, die fu¨r das Systemverhalten keine oder nur untergeordnete Rollen spielen und abstrahiert das Verhalten auf die wesentlichen, zu untersuchenden Merkmale. Modelle realer Systeme k¨onnen in bestimmten F¨allen als konservativ ­ also energieerhaltend ­ angen¨ahert werden, und unter Umst¨anden kann eine lineare Formulierung brauchbare L¨osungen liefern. Wichtige Ph¨anomene der meisten realen, mechanischen Systeme sind jedoch auf Nichtlinearit¨aten und Dissipation mechanischer Systemenergie zuru¨ckzufu¨hren, was oft eine · dissipative und · nichtlineare Modellierung erfordert. Dissipativ-- also energieverzehrend-- sind reale Systeme, weil die mechanische Systemenergie entweder durch mechanische Reibung, D¨ampfung oder durch inelastische St¨oße dissipiert 1 wird. Lineare Systeme sind zwar in bestimmten F¨allen in der Lage, einen ersten Eindruck des Systemverhaltens zu geben. Fu¨r eine Vorhersage s¨amtlicher tats¨achlicher Ph¨anomene eines realen Modells sind nichtlineare Terme in den Modell-Differentialgleichungen jedoch unerl¨asslich. Fu¨r die Analyse hat das Konsequenzen: w¨are ein dissipatives System linear, so bes¨aße es genau eine stabile L¨osung, zu der jegliches Langzeitverhalten hin tendiert. Nichtlineare Systeme beinhalten besondere Ph¨anomene wie die Koexistenz mehrerer mec 1 h D a i n ss is i c p h a e ti r on Enbeerdgeieu,tedtiewd¨oardtluicrhch"ZinerWstr¨aerumuengu"mugnedwabnedzeelitchwniertd damit die Zerstreuung und somit als fu¨r das mechanische System verloren gilt. 31

KAPITEL 5. MODELLIERUNG Attraktoren, nichtperiodisches(also chaotisches) Verhalten, selbst wenn das System periodisch angeregt wird, und unterschiedliche Arten von Bifurkationen, also pl¨otzliches A¨ ndern von L¨osungseigenschaften wie Vielfachheit, Form, Art, Ausdehnung und Stabilit¨at, wenn ein Systemparameter seinen Wert infinitesimal 2 variiert. Dergleichen Ph¨anomene nichtlinearer Systeme verlangen besondere Analysemethoden, die fu¨r lineare Systeme nicht n¨otig sind. Systemzustand und Parameter Das Verhalten eines Systems h¨angt von System- und Prozessparametern ab. Systemparameter bemessen Eigenschaften, die dem System fest anhaften. Prozessparameter sind Gr¨oßen, die sich w¨ahrend des Betriebes ergeben oder eingestellt und ver¨andert werden k¨onnen. Von den folgenden Parametern h¨angt das Verhalten des UltraschallStoßbohrers ab: Systemparameter sind · Materialeigenschaften wie Dichte oder Elastizit¨at des Transducers(Oszillators), des Bohrers und der freien Stoßmasse, · Oberfl¨achenbeschaffenheiten(Rauhigkeit, Kru¨mmung, H¨arte) der Kontaktfl¨achen, die in St¨oße involviert sind, · Geometriegr¨oßen(L¨angen und Querschnittsfl¨achen der K¨orper) und · mechanische Vorspannung der Piezokeramiken. Prozessparameter sind · elektrische Anregespannung, · Anregefrequenz, · Anpressdruck des Bohrers auf das Bohrgut(Gestein), · Flugabstand der Stoßmasse und · Umgebungstemperatur. Ein mathematisches Modell will das Verhalten eines Systems beschreiben. Der Begriff des Systemverhaltens bezeichnet die Angabe des gesamten Systemzustands zu allen Zeiten. Unter dem Systemzustand versteht man die zahlenwertige Angabe s¨amtlicher Zustandsgr¨oßen, die den Systemzustand eindeutig definieren. Zeitabh¨angige, ver¨anderliche Zustandsgr¨oßen k¨onnen z. B. · Schwerpunktsgeschwindigkeiten aller starren K¨orper in drei Dimensionen, · Schwerpunktspositionen aller starren K¨orper in drei Dimensionen, 2 infinitesimal: unendlich klein werdend 32

· Winkelgeschwindigkeiten aller starren K¨orper in drei Dimensionen, · Winkellagen aller starren K¨orper in drei Dimensionen oder · Verformungen der elastischen K¨orper in drei Dimensionen sein. Ein Modell mit dem Anspruch der Vollst¨andigkeit wu¨rde versuchen, in Abh¨angigkeit aller System- und Prozessparameter die zeitlichen Ver¨anderungen all dieser Zustandsgr¨oßen miteinander in Beziehung zu setzen. Der Aufwand, ein Systemverhalten derart vollst¨andig zu beschreiben, w¨are enorm. Es erg¨abe sich eine Komplexit¨at hohen Ausmaßes, die numerisches L¨osen kaum noch m¨oglich machte. Schließlich w¨are die Flut an Ergebnisdaten nicht mehr auswert- oder begreifbar. Daher ist es das Ziel des Wissenschaftlers zu erkennen, welche Systemgr¨oßen fu¨r eine Analyse des Verhaltens berechnet werden mu¨ssen und welche Parameter die Entwicklung der Systemgr¨oßen maßgeblich beeinflussen. Dabei ist es ein bew¨ahrtes Vorgehen, sich von sehr einfach gehaltenen Modellen, die reales Verhalten noch zu ungenau wiedergeben werden, hin zu aufw¨andigeren Modellen zu verbessern, um die Komplexit¨at so gering wie m¨oglich und nur so hoch wie n¨otig zu halten. Modellierungsvorgehen fu¨r den Ultraschall-Stoßbohrer In Kapitel 3 wurde das reale System des Ultraschall-Stoßbohrers vorgestellt. Das vorliegende Kapitel wird sich der modellhaften Beschreibung von dessen Eigenschaften in mathematischer Sprache widmen. Es werden zwei Modelle hergeleitet werden: zun¨achst ein einfaches, grundlegendes Modell in Unterkapitel 5.1. Es dient in erster Linie einer Einfu¨hrung in die Thematik der nichtglatten Systemanalyse. Es ist bewusst einfach gehalten und zeigt dem Leser so das Potential der mengenorientierten Numerik, mit der die Modelle in Kapitel 7 untersucht und studiert werden. Die Zahl der Zustandsgr¨oßen bleiben hier auf zwei beschr¨ankt; dieses Modell fokussiert ausschließlich auf die Bewegung der freien Stoßmasse. Jede Modellierung der Wirklichkeit abstrahiert das tats¨achliche Verhalten auf diejenigen Vorg¨ange, die fu¨r den Zweck der angestrebten Untersuchungen relevant erscheinen. Gleichzeitig stellt man Zusammenh¨ange idealisiert dar. So reduziert man z. B. einen K¨orper auf eine Punktmasse, wenn Wellenausbreitungen innerhalb des K¨orpers und Verformungen vernachl¨assigbar sind und man annehmen darf, dass er nur translatorische (und keine rotatorischen) Bewegungen ausfu¨hrt. In dieser Vernachl¨assigung stecken zwei Annahmen: erstens nimmt man an, dass sich das gesamte System, in dem dieser K¨orper eine Komponente ist, genau so verh¨alt, wie wenn es sich um einen starren, also unverformbaren K¨orper handeln wu¨rde. Zweitens geht man davon aus, dass Rotationen des K¨orpers mit endlicher Ausdehnung entweder nicht stattfinden oder so klein sind, dass sie das untersuchte Verhalten des Gesamtsystems nicht beeinflussen. Die Bewegungsgleichungen k¨onnen dann einfach gehalten werden, was eine 33

KAPITEL 5. MODELLIERUNG wichtige Voraussetzung fu¨r eine Analyse darstellt. Prinzipiell ist es wichtig, der Modellierung zugrundeliegende Annahmen zu validieren. Fu¨r die Validierung des einfach gehaltenen Modells aus Kapitel 5.1 wurde ein reales System im Labor untersucht. Der Pru¨fstand und die Messverfahren werden in Kapitel 6 beschrieben. Der Vergleich des Modells mit der Wirklichkeit wird im Unterkapitel 6.5 durchgefu¨hrt. Das zweite Modell, das anschließend in Unterkapitel 5.2 hergeleitet wird, stellt eine Verfeinerung des ersten Modells in einigen Aspekten dar: die im realen Bohrsystem freie Beweglichkeit des Oszillators wird in dieses Modell mit zwei zus¨atzlichen Zustandsgr¨oßen aufgenommen. Außerdem werden der Prozessparameter Anpresskraft und der Systemparameter Oszillatormasse zus¨atzlich in die Modellgleichungen einfließen. Unterkapitel 5.3 geht im Detail auf die Implementierung eines L¨osungsalgorithmus ein, der die fu¨r die Modellgleichung erforderliche Bestimmung der Schnittpunkte einer gegebenen Gerade mit einer Kosinusfunktion und einer gegebenen Parabel mit einer Kosinusfunktion vornimmt. Im Kontext dieser Anwendung ist es nicht von Belang, alle Schnittpunkte zu finden, sondern es ist wichtig, den numerisch kleinsten der Menge an Schnittpunkten zu identifizieren, da dieser den Zeitpunkt des Stoßkontaktes zwischen zwei Stoßk¨orpern des Stoßbohrsystems beschreibt. 5.1 Herleitung eines Modells 2. Ordnung Die Bewegung der frei-fliegenden Masse innerhalb des Stoßbohrers ist von zentraler Bedeutung fu¨r die Funktionsweise des Bohrsystems. In diesem Unterkapitel wird ein einfaches Modell in Form von mathematischen Bewegungsgleichungen aufgestellt. Es modelliert ein einfaches Schwingstoßsystem, das in den folgenden Kapiteln simuliert und schließlich unter Anwendung der mengenorientierten Methoden untersucht wird. 5.1.1 Zeitdiskrete Modellierung Dieser Abschnitt leitet ein System zeitdiskreter Bewegungsgleichungen 2. Ordnung ­ also mit zwei Zustandsvariablen-- her. Bild 5.1 zeigt eine Skizze des Modells. Die Bohrrichtung ist in der schematischen Skizze bewusst horizontal angenommen, um gedanklich den Einfluss eines Gravitationsfeldes auf die Bewegung der Masse zu vernachl¨assigen. Wegen der kurzen Flugzeiten der Masse zwischen zwei Richtungswechseln ist diese Annahme selbst fu¨r das Bohren vertikaler L¨ocher n¨aherungsweise gerechtfertigt bzw. gilt im Falle horizontalen Bohrens. Der Bohrstab ist hier mit seiner systeminneren Kontaktseite als starre, unbewegliche Wand(Prallplatte) dargestellt. Die Starrk¨orperbewegungen des gesamten Oszillators werden in diesem Modell vernachl¨assigt, um zun¨achst mit wenigen Freiheitsgraden die Komplexit¨at des Modells zu begrenzen und sich ganz auf die Dynamik der Masse zu konzentrieren. Demzufolge ist die Anregung durch die Oszillatorspitze allein als eine harmonische Fußpunkterregung dargestellt. 34

5.1. HERLEITUNG EINES MODELLS 2. ORDNUNG Prallplatte u Oszillatorspitze W k = - u 0  sin t k x - U k V k 1 h 2 Bild 5.1: Einfaches Modell 2. Ordnung des Stoßbohrers Von 5 Systemparametern h¨angt das Modellverhalten ab: Prozessparameter sind · u 0 : Anregeamplitude und · : Anregekreisfrequenz, Modellparameter sind · 1 , 2 : Stoßzahlen nach Newton und · h : Abstand zwischen Nulllage der harmonischen Fußpunkterregung und Prallplatte. W k , U k , V k bezeichnen die Geschwindigkeiten der Oszillatorspitze, der Stoßmasse vor und der Stoßmasse nach dem Stoß k ; dabei zeigt die positive Bewegungsrichtung nach rechts. u symbolisiert die Position der Oszillatorspitze. In vielen F¨allen werden Bewegungsgleichungen in zeitkontinuierlicher Form, also als Differentialgleichungssystem aufgestellt. Dies kann sinnvoll sein, wenn eine stetige, kontinuierliche und glatte Bewegung beschrieben werden soll. Bei nichtglatten Systemen wie dem Vorliegenden wechselt das System bei jedem Stoßkontakt die Bewegungsgleichungen. Die Dynamik l¨asst sich dann nicht mehr in Form eines geschlossenen Differentialgleichungssystems beschreiben. Es eignet sich hier eine zeitdiskrete Formulierungsweise. Der Systemzustand wird nicht mehr kontinuierlich zu jeder Zeit t definiert, sondern nur zu diskreten Zeitpunkten. Fu¨r den Stoßbohrer ist es sinnvoll, diese Zeitpunkte in die Augenblicke zu legen, in denen die Stoßmasse auf den Oszillator st¨oßt. In diesen Momenten wird der Systemzustand durch die Phasenlage des Stoßes in Bezug auf die harmonische Anregung und die Abfluggeschwindigkeit der Masse nach dem Stoß eindeutig definiert. Aus diesen zwei Zustandsgr¨oßen l¨asst sich problemlos der gesamte Bewegungsverlauf bis zum n¨achsten Masse-Oszillator-Stoß ableiten. Ein System von zwei Differenzengleichungen wird aufgestellt. Es definiert eine Abbildungsvorschrift, die vom Stoßzustand j auf den Stoßzustand j + 1 abbildet. A¨ hnliche Vorgehensweisen zur Modellierung von Stoßsystemen 35

KAPITEL 5. MODELLIERUNG sind aus[Holmes 1982],[Guckenheimer& Holmes 1983],[Brogliato 1996], [Mello, Tufillaro 1985] oder[Neumann& Sattel 2007] bekannt. Folgende vereinfachende Annahmen wurden fu¨r diese erste Modellbildung getroffen: · die harmonische Oszillatorschwingung l¨asst sich nicht durch die St¨oße der Stoßmasse beeinflussen, d. h. die Verschiebung u ( t ) der Oszillatorspitze wird als harmonische Fußpunkterregung fest vorgegeben. 3 · die Starrk¨orperbewegungen des Bohrstabes werden vernachl¨assigt, die Kontaktfl¨ache des Bohrstabes mit der Kugel wird als eine Prallplatte abgebildet. · die Starrk¨orperbewegungen des gesamten Oszillators werden in diesem Modell zun¨achst ignoriert. · das systeminnere Bohrstabende-- hier als Prallplatte modelliert-und die Mittelposition u = 0 der Oszillatorspitze befinden sich in konstantem Abstand h zueinander. · der Stoßvorgang zwischen dem Oszillator und der Stoßmasse wird als unendlich kurz angesehen und durch Newtons Stoßhypothese mit Stoßzahl 1 beschrieben. · der Stoßvorgang zwischen der Stoßmasse und dem systeminneren Bohrstabende(Prallplatte) wird ebenfalls durch die Stoßhypothese mit Stoßzahl 2 beschrieben. · Gravitation beeinflusst die Bewegungen nicht 4 . 5.1.2 Herleitung zeitdiskreter Bewegungsgleichungen Zwischen zwei St¨oßen fu¨hrt die modellierte Stoßmasse eine unged¨ampfte Bewegung konstanter Geschwindigkeit aus. Somit ist es angebracht, ihr Langzeitverhalten wie oben beschrieben als zeitdiskrete Entwicklung aufeinanderfolgender Kontakte mit dem Oszillator zu beschreiben. Anstelle von zeitkontinuierlichen Zustandsgr¨oßen x ( t ) , x  ( t ) werden zur eindeutigen Beschreibung des Systemzustands die Zustandsvariablen · t k , Zeitpunkt des k -ten Kontaktes zwischen Stoßmasse und Oszillator und · V k , Geschwindigkeit der Stoßmasse nach dem k -ten Stoß am Oszillator verwendet. Die Abbildung( t k , V k ) ( t k+1 , V k+1 ) wird als Differenzengleichungssystem definiert, das im Folgenden hergeleitet wird. 3 Zur Validierung dieser Annahme siehe Unterkapitel 6.1.2, Seite 56. 4 Weil die Gravitationskraft nicht auf die Stoßmasse wirkt, fließt die Gr¨oße der Stoßmasse in dieses Modell nicht ein. 36

5.1. HERLEITUNG EINES MODELLS 2. ORDNUNG Die Fußpunktverschiebung u ( t ) des harmonisch schwingenden Oszillators sei vorgegeben als u ( t )= u 0 cos t. (5.1) Fu¨r die Herleitung der Modellgleichungen fu¨hren wir drei Geschwindigkeiten als Hilfsgr¨oßen ein, die auch in Skizze 5.1 enthalten sind: U k := x  ( t k ) , V k := x  ( t k + ) , W k := u  ( t k ) , (5.2) U k und V k bezeichnen die Geschwindigkeit der Stoßmasse vor bzw. nach dem Stoß k am Oszillator, und W k ist die Oszillatorgeschwindigkeit zum Stoßzeitpunkt. Die St¨oße der Stoßmasse mit dem Oszillator bzw. der Prallplatte sind teilelastisch. Die Energiedissipation beim Stoß wird mittels Stoßzahlen nach Newton beru¨cksichtigt, die wie folgt definiert werden: 1 : Stoßzahl zwischen Oszillator und Stoßmasse: 1 := V k U k W k W k V k+1 =(1+ 1 ) W k+1 - 1 U k+1 (5.3) 2 : Stoßzahl zwischen Stoßmasse und Bohrstab: 2 := U k+1 V k U k+1 = - 2 V k (5.4) Mit dem unver¨anderlichen Abstand h zwischen der Mittelposition der Oszillatorspitze und der Prallplatte kann die Zeitspanne zwischen zwei aufeinanderfolgenden St¨oßen am Oszillator berechnet werden, wenn man von der M¨oglichkeit eines Mehrfachstoßes 5 absieht: t k+1 - t k = t Oszillator Prallplatte + t Oszillator Prallplatte = h- u ( t k ) V k + h- u ( t k+1 ) - U k+1 . (5.5) Einsetzen von U k+1 aus Gl.(5.4) und u ( t ) aus Gl.(5.1) und Umstellen fu¨hrt auf 0= t k+1 - t k 1 V k h 1+ 1 2 - u 0 cos t k u 0 2 cos t k+1 . (5.6) Die transzendente 6 Gleichung(5.6) kann explizit nicht nach der gesuchten Zustandsgr¨oße t k+1 aufgel¨ost werden. Die zweite zeitdiskrete Zustandsgleichung gibt die Folge der Abfluggeschwindigkeiten V k nach jedem MasseOszillator-Stoß an. Sie folgt aus Gl.(5.3) durch Einsetzen von(5.2) und (5.4): V k+1 = 1 2 V k (1+ 1 ) u 0  sin t k+1 . (5.7) 5 Ein sogenannter Mehrfachstoß meint in diesem Zusammenhang einen Mehrfachstoß am Oszillator und bezeichnet das mehrfache Aufeinandertreffen von freier Stoßmasse und Oszillator, ohne dass die Stoßmasse dazwischen Kontakt mit der Prallplatte hat. Dies geschieht immer dann, wenn V k < 0 oder wenn V k einen betragsm¨aßig kleinen Wert erreicht. Dann n¨amlich wird die Stoßmasse den Aktionsradius( u 0 ) der Oszillatorspitze nicht verlassen, bevor sie von der harmonischen Schwingung wiederholt eingeholt wird. 6 transzendent: siehe Fußnote auf Seite 45. 37

KAPITEL 5. MODELLIERUNG Somit ist das System von zwei Differenzengleichungen mit den Zustandsgr¨oßen x k =[ V k , t k ] T definiert durch die implizite Gleichung g ( x k+1 , x k )= 0 (5.8) mit g ( x k+1 , x k ):= t k+1 - t k 1 V k h 1+ 1 2 - u 0 cos t k u 0 2 cos t k+1 . V k+1 - 1 2 V k +(1+ 1 ) u 0  sin t k+1 (5.9) 5.1.3 Notwendigkeit transzendenter Bewegungsgleichungen Wu¨nschenswert ist es, die die Bewegung beschreibende Differenzengleichung (5.9) in einfacherer, also expliziter Form zu formulieren, um sie analytisch l¨osen und eine Zeitreihe berechnen zu k¨onnen. Guckenheimer und Holmes[Guckenheimer & Holmes 1983] bzw. Holmes[Holmes 1982] treffen im Bouncing-Ball Problem(siehe Unterkapitel 4.3) die Annahme, die Flugh¨ohe der Masse sei der Oszillatoramplitude um Gr¨oßenordnungen u¨berlegen, weshalb die sich ¨andernde Oszillatorspitzenposition vernachl¨assigt werden kann. Im realen Stoßbohrer stellt sich dieser Abstand jedoch variabel ein und nimmt mitunter sehr kleine Werte an. Fu¨r das vorliegende, einfach gehaltene Modell wollen wir uns jedoch auf feste Abst¨ande h in der Gr¨oßenordnung von Millimetern beschr¨anken, wie sie auch sp¨ater in einem Laborexperiment eingestellt werden sollen(Unterkapitel 6.1). Realistische Amplitudenwerte der Oszillatorspitze liegen im Bereich von ca. 10 µ m. Wir k¨onnten vereinfachend die Annahme treffen, dass die Masse nahezu auf einer H¨ohe von 0 m mit dem Oszillator kollidiert. Der relative Fehler, den wir in Kauf nehmen mu¨ssten, l¨age lediglich bei s s < u 0 h = 10 µ m 1mm = 1%(5.10) und erscheint hinnehmbar. Wird also in Gleichung(5.5) u ( t k )  u ( t k+1 ) 0(5.11) angenommen, folgt anstelle von Gleichung(5.6) fu¨r die Flugdauer t k+1 - t k = h V k 1+ 1 2 mit ( t k+1 - t k ) t k+1 - t k 1% . (5.12) Wie setzt sich dieser geringe Fehler nun aber in der unver¨andert geltenden Gl.(5.7) V k+1 = 1 2 V k (1+ 1 ) u 0  sin t k+1 (5.13) fort? Die Flugdauer t k+1 - t k mit einem relativen Fehler von 1% geht hier in die Sinus-Funktion ein. Wegen ihrer Periodizit¨at muss der absolute Fehler 38

5.1. HERLEITUNG EINES MODELLS 2. ORDNUNG des Arguments herangezogen und auf eine Periodendauer bezogen werden. In Laborversuchen wird sich sp¨ater zeigen, dass fu¨r h = 5 mm Flugdauern von bis zu 50 ms erreicht werden k¨onnen(siehe Bild 6.6). Bei einer Frequenz von 20 kHz entspricht ein 1%-Fehler der Flugdauer bereits dem 10-fachen einer Periodendauer. Somit k¨onnte bereits jedes Argument des Sinus innerhalb einer Periode fehlerhaft erreicht werden. Der relative Fehler der Geschwindigkeits-Zustandsgr¨oße kann wie folgt nach oben abgesch¨atzt werden: V k+1 V k+1 (1+ 1 ) u 0 (max[sin t k+1 ] 20 min[sin t k+1 ] 20 ) 1 2 V k (1+ 1 ) u 0  sin t k+1 = 2(1+ 1 ) u 0 1 2 V k (1+ 1 ) u 0  sin t k+1 2(1+ 1 ) u 0 1 2 V k (1+ 1 ) u 0 = 2 V k 2 u 0  1+1/ 1 1 (5.14) Fu¨r 1 [0 , 1] und 2 [0 , 1] ist der Term 2 1+ 1 / 1 [0 , 0 . 5] . Fu¨r = 2  · 20 kHz 1 . 3 · 10 5 rad/s, u 0 [1 , 10] µ m und V k [0 . 1 , 10] m/s ist der Term V k u 0 [0 . 08 , 1 . 3] . Somit ist der gesuchte relative Fehler V k+1 V k+1 [ 200 , 8]% . (5.15) Mit bis zu 200% Abweichung vom tats¨achlichen Wert fu¨r die Abfluggeschwindigkeit V k folgt die Notwendigkeit, fu¨r die genaue Abflugposition die momentane Auslenkung u ( t k ) der Oszillatorspitze zum Abflugzeitpunkt mit zu beru¨cksichtigen, was auf die transzendente, implizite Differenzengleichung(5.6) fu¨hrt. Diese Notwendigkeit k¨onnte bezu¨glich des realen Experiments allenfalls dadurch relativiert werden, dass eine pr¨azise Festlegung des Abstandes h kaum m¨oglich ist; weiterhin wird sp¨ater eine starke Empfindlichkeit des Verhaltens bezu¨glich der Stoßzahlen festgestellt(Kapitel 6.2.3), die im Modell als konstant angenommen werden. 5.1.4 Mehrfachst¨oße Die Abbildungsfunktion(5.9) beschreibt die Bewegung der Stoßmasse unter der Annahme, dass sie nach einem Stoß am Oszillator direkt auf die 39

KAPITEL 5. MODELLIERUNG Prallplatte zufliegt. Tats¨achlich kann es jedoch vorkommen, dass sie vorher mehrfach auf den Oszillator st¨oßt. Diesen Vorgang bezeichnet man als Mehrfachstoß. Sein Auftreten ist selten und wird von kleinen oder negativen Abfluggeschwindigkeiten V k verursacht. In solchen F¨allen verl¨asst die Masse den Aktionsradius des schwingenden Oszillators nicht, bevor dessen Spitze wiederholt auf die Masse trifft. Bild 5.2 zeigt m¨ogliche F¨alle eines doppelten und eines dreifachen Mehrfachstoßes. Nach jedem Stoß der Masse am Oszilx ( t ) u ( t ) t Bild 5.2: Zeitverlauf der Stoßmasse x ( t ) sowie des Oszillators u ( t ) w¨ahrend Mehrfachst¨oßen lator lautet ihre kontinuierliche Bewegungsgleichung bis zum n¨achsten Stoß x ( t )= V k ( t- t k )+ u 0 cos t k . (5.16) Ein Mehrfachstoß tritt genau dann auf, wenn ein t k+1 existiert mit x ( t k+1 )= u ( t k+1 ) , (5.17) was auf die Bedingung u 0 V k (cos t k+1 cos t k )= t k+1 - t k , t k < t k+1 (5.18) fu¨hrt. Numerisch zu bestimmen, ob ein t k+1 existiert, welches Gl.(5.18) erfu¨llt, stellt eine Herausforderung dar. Existiert solch ein t k+1 nicht, so findet kein Mehrfachstoß statt und es ist die Abbildungsgleichung(5.9) zu l¨osen, um t k+1 zu bestimmen. 5.2 Herleitung eines verfeinerten zeitdiskreten Modells Das in Unterkapitel 5.1 untersuchte einfache Modell 2. Ordnung gab einen konstanten Abstand zwischen Oszillatorschwerpunkt und Prallplatte(Bohrer) vor. Eine Eigenschaft des USDC-Prototypen ist jedoch ein sich selbst¨andig einstellender Abstand. Das folgende Unterkapitel leitet ein verfeinertes Modell her, das diesem Umstand gerecht wird. Es beru¨cksichtigt die Starrk¨orperbewegung des Aktors. Da in diesem Modell auch die Stoßmasse durch den Einfluss der Schwerkraft beschleunigt wird, geht hier auch 40

5.2. HERLEITUNG EINES VERFEINERTEN ZEITDISKRETEN MODELLS ihre Masse als Parameter in die Gleichungen mit ein. Zur Erinnerung-im ersten, einfacheren Modell spielte die Masse der Stoßmasse keine Rolle, da hier kein Einfluss der Schwerkraft vorausgesetzt wurde(wie es etwa bei horizontalem Bohren der Fall ist). Somit wurde die Bewegung allein durch kinematische 7 Gesetze bestimmt. 5.2.1 Aufstellen des verfeinerten Modells 4. Ordnung F A M x ( t ) u ¯( t )= u 0 cos t u ( t ) m g y ( t ) Oszillator Bild 5.3: Verfeinertes Modell 4. Ordnung Bild 5.3 zeigt den Oszillator, dessen unteres Ende infolge der angeregten ersten L¨angseigenmode(in der Skizze rechts vom Oszillator dargestellt) harmonisch schwingt, und die Stoßmasse(dargestellt als Kugel). Die Positionen x ( t ) des Oszillators sowie y ( t ) der Stoßmasse werden als H¨ohe u¨ber dem Boden gemessen. So kann die Auslenkung des Oszillator-St¨oßels u ( t ) als eine an den Oszillator gekoppelte Bewegung u ¯( t ):= u 0 cos t formuliert werden: u ( t )= x ( t ) - u ¯( t ) . (5.19) Um ein Eindringen der K¨orper ineinander auszuschließen, wird voraus7 Die Kinematik ist die Lehre von der Bewegung von K¨orpern. Nur geometrische U¨ berlegungen fließen in die kinematischen Gesetze ein. Im Ggs. dazu ist die Kinetik die Lehre von der Wirkung von Kr¨aften. 41

KAPITEL 5. MODELLIERUNG gesetzt: 0  y ( t )  t und y ( t )  u ( t )  t. (5.20) Auf den Oszillator wirke die Anpresskraft F A . Diese kann positiv, negativ oder Null sein. Auf beide K¨orper wirke außerdem ihre Gewichtskraft. Von 7 Systemparametern h¨angt das Modellverhalten ab: Prozessparameter sind · u 0 : Schwingamplitude der Oszillatorspitze, · : Anregekreisfrequenz, · F A /M : Oszillatorschwerpunktsbeschleunigung aus Anpresskraft und · g : Erdbeschleunigungskomponente in L¨angsrichtung; Modellparameter sind · 1 , 2 : Stoßzahlen nach Newton und · µ = M m : Massenverh¨altnis Oszillatormasse zu Stoßmasse. Vier Zustandsgr¨oßen Im Betrieb des Bohrers kommt es zu Stoßkontakten zwischen den drei Kontaktpartnern Oszillatorspitze, freie bezeichnet) und dem inneren Ende Stoßmasse(im des Bohrstabes. FDoilegseensdEenndaelsde"sKBuogherl"stabes soll im Starrheit und Folgenden Fixiertheit hminitzu"Pwreaisllepnl,adttiee"wbirezineidchienseetr werden, um auf seine Modellierung voraussetzen. Zwischen zwei aufeinanderfolgenden St¨oßen fu¨hren die beiden beweglichen K¨orper Oszillator und Kugel freie, beschleunigte Bewegungen durch, die dem zweiten Newtonschen Axiom(Kraft pro Masse= Beschleunigung) gehorchen und durch Anfangsgeschwindigkeit sowie Anfangsauslenkung eindeutig festgelegt sind. Somit genu¨gt zur eindeutigen Beschreibung der Bewegungen der beiden K¨orper allein die Angabe des Bewegungszustands beider K¨orper nach jedem Stoß sowie die Angabe des Zeitpunktes, an dem der jeweilige Stoß stattfand. Bis zum Auftreten des n¨achsten Stoßes k¨onnen die Bewegungen damit als gegeben angesehen werden. Folglich wird der dynamische Zustand des Systems nur zu diskreten Zeitpunkten beschrieben, n¨amlich genau zu den Zeitpunkten, an denen ein Stoß zwischen Kugel und Oszillator auftritt. Durch die folgenden vier Zustandsgr¨oßen wird der Systemzustand demnach vollst¨andig beschrieben: t k : x ( t k ): x  ( t + k ): y ( t + k ): Stoßzeitpunkt Oszillator-Position im Stoßzeitpunkt Oszillator-Geschwindigkeit direkt nach dem Stoß Kugel-Geschwindigkeit direkt nach dem Stoß 42

5.2. HERLEITUNG EINES VERFEINERTEN ZEITDISKRETEN MODELLS Anmerkung: Da Kugel und Oszillator im Stoßzeitpunkt in Kontakt sind, ist wegen x ( t k )= y ( t k ) die Kugel-Position y ( t k ) im Stoßzeitpunkt ebenso gegeben. Abku¨rzend werden die Symbole x k := x ( t k ) x + k := x  ( t + k ) y k + := y ( t + k ) eingefu¨hrt; damit lautet der Zustandsvektor x k := t k , x k , x + k , y k + T . (5.21) Wie im Fall des einfachen, dynamischen, 2-dimensionalen Modells aus Unterkapitel 5.1 wird die Dynamik also durch ein zeitdiskretes dynamisches System definiert. Das heißt, die Bewegungsgleichungen werden durch ein System von vier Differenzengleichungen gebildet, welche im Weiteren aufgestellt werden sollen. Die Differenzengleichungen bilden den Zustand des Systems im Zeitpunkt t k auf den Zustand im Zeipunkt t k+1 ab. Somit handelt es sich hier um die Abbildung: t k t k+1 x k x + k ---A b b il d u n g v o n -Stoß k nach Stoß k + 1 x k+1 x + k+1 y k + y k ++1 Die Abbildungsvorschrift wird in den folgenden Abschnitten hergeleitet. 5.2.2 Berechnung des Stoßzeitpunktes t k +1 Die Bewegung der Kugel y ( t ) nach Stoß k ist durch folgendes Anfangswertproblem(AWP) gegeben: y ¨= - g y ( t k )= y k + y ( t k )= u ( t k ) Bewegungsdifferentialgleichung Anfangsgeschw.(Abfluggeschw. v. Oszill.)(5.22) Anfangspos.(Kontakt mit Oszillator) Fu¨r die Anfangsposition gilt u ( t k )= x ( t k ) - u ¯( t k )= x k - u 0 cos t k . (5.23) L¨osen des AWPs(5.22) fu¨hrt auf die Bewegungsgleichung fu¨r die Kugel nach Stoß k : y ( t )= g 2 ( t- t k ) 2 + y k + ( t- t k )+ x k - u 0 cos t k , t k  t< t k+1 . (5.24) 43

KAPITEL 5. MODELLIERUNG Die Bewegung des Oszillators x ( t ) nach Stoß k ist gegeben durch das AWP x ¨= - g- F A /M x  ( t k )= x + k x ( t k )= x k Bewegungsdifferentialgleichung Anfangsgeschw.(Abfluggeschw. v. Kugel)(5.25) Anfangsposition L¨osen des AWPs(5.25) fu¨hrt auf die Bewegungsgleichung fu¨r den Oszillator nach Stoß k : x ( t )= g + F A /M 2 ( t- t k ) 2 + x + k ( t- t k )+ x k , t k  t< t k+1 . (5.26) Nun gilt es zu bestimmen, wann der n¨achstfolgende Stoß auftritt: gesucht ist t k+1 . Dabei muss zwischen zwei m¨oglichen F¨allen unterschieden werden: 1. M¨oglichkeit: Stoß Kugel-Prallplatte zum Zeitpunkt t W k mit t k < t W k < t k+1 Stoßbedingung: y ( t W k )= 0 2. M¨oglichkeit: Mehrfachstoß Kugel-Oszillator zum Zeitpunkt t k+1 Stoßbedingung: y ( t k+1 )= u ( t k+1 ) Unter der Annahme, nach dem Abflug der Kugel vom Oszillator geschieht der n¨achste Stoß der Kugel an der Prallplatte(1. M¨oglichkeit), wird t W k , also der Zeitpunkt des Stoßes zwischen Kugel und Prallplatte, aus folgender Kontaktbedingung berechnet: y ( t W k )= 0 , gesucht: t W k . (5.27) Die Bedingung fu¨hrt auf ein quadratisches Polynom in( t W k - t k ), das sich zu t W k 1 , 2 = t k + 1 g y k + ( y k + ) 2 + 2 g ( x k - u 0 cos t k ) = u( t k )> 0 (5.28) l¨ost. Die Minus-L¨osung interessiert nicht, da sie zu einem t W k mit t W k < t k fu¨hren wu¨rde, was physikalisch nicht m¨oglich ist. Unter der Annahme, es komme zum Mehrfachstoß(2. M¨oglichkeit), also einem wiederholten Kontakt zwischen Kugel und Oszillator ohne zwischenzeitlichen Kontakt mit der Prallplatte, wird der Kontaktzeitpunkt t M k+F1S gesucht, an dem dieser Mehrfachstoß(MFS) auftritt. Die Bedingung fu¨r einen Oszillator-Kugel Kontakt ist y ( t M k+F1S )= u ( t M k+F1S ) = x ( t M k+F1S ) - u ¯( t M k+F1S ) . 44

5.2. HERLEITUNG EINES VERFEINERTEN ZEITDISKRETEN MODELLS Hier werden die Bewegungs-Funktionen y ( t )(5.24), x ( t )(5.26) und u ¯( t ) eingesetzt:  g 2 ( t M k+F1S - t k ) 2 + y k + ( t M k+F1S - t k )+ x k - u 0 cos t k = g + F A /M 2 ( t M k+F1S - t k ) 2 + x + k ( t M k+F1S - t k )+ x k - u 0 cos t M k+F1S F A 2 M t M k + F 1 S - t k 2 + y k + - x + k t M k + F 1 S - t k = - u 0 (cos t M k + F 1 S cos t k )(5.29) Gesucht ist nun das kleinste, reelle t M k + F 1 S mit t k < t M k+F1S , welches die letzte Gleichung(5.29) l¨ost, also: t M k+F1S = min t M k+F1S | t k < t M k+F1S  y ( t M k+F1S )= u ( t M k+F1S ) . (5.30) Diese Aufgabe stellt kein triviales Problem dar, da Gleichung(5.29) transzendent 8 ist. Allgemein l¨asst sich der Problemtypus folgendermaßen e fo in rm em ulie K r o en si : nu"sF."ind D e ie s¨a A m n t z li a c h h l e Schnittpunkte zwischen einer der Schnittpunkte variiert in Parabel und Abh¨angigkeit von den Parametern zwischen Null und unendlich. Die Behandlung dieses Problems wird im folgenden Unterkapitel 5.3 beschrieben. Das aus Gleichung(5.29) gefundene t M k+F1S wird mit dem unter der Prallplattenstoß-Annahme(1. M¨oglichkeit) gefundenen t W k (Gleichung (5.28)) verglichen: Falls t M k+F1S < t W k , tritt ein Mehrfachstoß am Oszillator-- und kein Prallplattenstoß-- auf und t k+1 ist gefunden: t k+1 = t M k+F1S . Andernfalls( t W k < t M k+F1S  t M k+F1S = ) tritt ein Prallplattenstoß der Kugel zum Zeitpunkt t W k auf(1. M¨oglichkeit). In diesem Fall muss noch der nachfolgende Stoßzeitpunkt t k+1 der Kugel mit dem Oszillator bestimmt werden. Zun¨achst wird angenommen, dass der dem Prallplattenstoß folgende Stoß der Kugel mit dem Oszillator erfolgt. (Es gibt seltene Situationen, in denen es zu einem Mehrfachstoß an der Prallplatte kommen kann. Dies kann insbesondere dann passieren, wenn die (Anpress)kraft F A negativ ist. Siehe hierzu die Fußnote 13 auf Seite 77.) 8 transzendente Gleichung:(lateinisch trancendere: u¨berschreiten) Eine Gleichung, die keine algebraische Gleichung ist bzw. eine Gleichung, die unterschiedliche transzendente Funktionen(Logarithmus, Exponentialfunktion, trigonometrische Funktion) enth¨alt. Eine transzendente Gleichung l¨asst sich nur in impliziter Form darstellen. L¨osungen k¨onnen numerisch oder grafisch, jedoch im Allgemeinen nicht analytisch gefunden werden. 45

KAPITEL 5. MODELLIERUNG Die absorbierte Energie beim Stoß der Kugel mit der Prallplatte wird nach der Newtonschen Stoßhypothese u¨ber die Stoßzahl 2 modelliert, welche wie folgt definiert ist: 2 := | v Kugel nach Stoß mit Prallplatte | | v Kugel vor Stoß mit Prallplatte | = y W + k - y W k (5.31) PDPrraaarllillnpplliaasttttt y ee W . u k Wn=deg yy e W + (n t k W y d k W )i k ed < Giee0Gsc(ehFswclhuinwgdiiningdkniegeiktgeaditteivrdeeKrRuKigcuehlgteunlnavcgoh)r dem Stoß dem Stoß und | y W + k | mit der mit der <| y W k | (aufgrund Energiedissipation in jedem realen Stoß) gilt 2 (0 , 1). Aus Gleichung(5.31) folgt mit Gleichung(5.24) y W + k = 2 g ( t W k - t k ) - y k + . (5.32) Die Bewegung y 2 ( t ) der Kugel nach dem Stoß mit der Prallplatte wird durch das Anfangswertproblem y ¨ 2 ( t )= - g y 2 ( t W k )= y W + k y 2 ( t W k )= 0 (5.33) beschrieben, das von y 2 ( t )= g 2 ( t- t W k ) 2 + y W + k ( t- t W k )(5.34) gel¨ost wird. Fu¨r den Kontaktzeitpunkt t k+1 zwischen Oszillator und Kugel muss gelten: y 2 ( t k+1 )= u ( t k+1 )= x ( t k+1 ) - u 0 cos t k+1  y 2 ( t k+1 ) > 0 . (5.35) Werden fu¨r y 2 ( t ) und u ( t ) die jeweiligen Bewegungsfunktionen eingesetzt, folgt wieder eine transzendente Gleichung in t k+1 , das heißt, es mu¨ssen wieder alle(reellen) Schnittpunkte zwischen einer Kosinusfunktion und einer Parabel gefunden werden. Daraus folgt t k+1 = min { t k+1 | t W k < t k+1  y 2 ( t k+1 )= u ( t k+1 )  y 2 ( t k+1 ) > 0 }. (5.36) Also ist der kleinste, reelle Schnittpunkt t k+1 mit t W k < t k+1 der gesuchte Stoßzeitpunkt. Wenn er nicht existiert(weil y 2 = u keine L¨osung besitzt oder weil y 2 ( t k+1 ) < 0), so tritt der seltene Fall eines Mehrfachstoßes an der Prallplatte auf. Beachte hierzu die Fußnote auf Seite 77. 5.2.3 Berechnen der verbleibenden Zustandsgr¨oßen Bislang wurde fu¨r die beiden m¨oglichen F¨alle des Oszillator-PrallplatteOszillator Stoßes(1. M¨oglichkeit) und des Oszillator-Oszillator Mehrfachstoßes(2. M¨oglichkeit) die Stoßzeit t k+1 , an der die Kugel wieder 46

5.2. HERLEITUNG EINES VERFEINERTEN ZEITDISKRETEN MODELLS vom Oszillator abprallt, berechnet. Fu¨r die vollst¨andige Angabe des Systemzustands fehlt noch die Berechnung der drei weiteren Zustandsgr¨oßen Abfluggeschwindigkeit der Kugel vom Oszillator nach dem Stoß, Oszillatorschwerpunktsgeschwindigkeit nach dem Stoß sowie Oszillatorbzw. Kugel-Position beim Stoß. Zur kinematischen Beschreibung des Stoßes am Oszillator wird die Stoßzahl 1 eingefu¨hrt: 1 := v K+ugel v K ugel v O+szillatorspitze v O szillatorspitze = y k ++1 - u + k+1 y (2) ( t k+1 ) - u  ( t k+1 ) (5.37) mit u + k+1 = x + k+1 - u ¯( t k+1 ). Die Geschwindigkeit der Kugel vor dem Stoß v K ugel wird beim Mehrfachstoß am Oszillator durch y ( t k+1 ) und beim Prallplattenstoß durch y 2 ( t k+1 ) erhalten. Die beiden F¨alle wurden in der letzten Gleichung durch die eingeklammerte (2) im Index der Geschwindigkeit gekennzeichnet. Eine zweite Gleichung fu¨r die zwei gesuchten Geschwindigkeiten folgt aus der Impulserhaltung: m · v K ugel + M · v O szillatorschwerpunkt = m · v K+ugel + M · v O+szillatorschwerpunkt , bzw. mit den hier verwendeten Symbolen: m · y (2) ( t k+1 )+ M · x  ( t k+1 )= m · y k ++1 + M · x + k+1 . (5.38) Darin sind m = Masse der Kugel, M = Masse des Oszillators. Es sei erw¨ahnt, das die Kugelmasse m erstmals in der Impulsbilanz(5.38) ins Modell miteinfließt. Die Oszillatormasse M hatte bei endlicher Anpresskraft F A bereits Einfluss auf die Bewegungsgleichung des Oszillators. Es wird das Massenverh¨altnis µ eingefu¨hrt: µ := M m . (5.39) Die beiden in den gesuchten Geschwindigkeiten linearen Gleichungen (5.37) und(5.38) werden nun nach den Geschwindigkeiten aufgel¨ost: x + k+1 = 1 1+ µ (1+ 1 )( y (2) ( t k+1 )+ u ¯( t k+1 ))+( µ- 1 ) x  ( t k+1 ) (5.40) y k ++1 = 1 1+ µ µ (1+ 1 )( x  ( t k+1 ) - u ¯( t k+1 )) ( µ 1 1) y (2) ( t k+1 ) . (5.41) Schließlich kann die Position des Oszillators im Stoßmoment durch Einsetzen von t k+1 in die Bewegungsgleichung fu¨r x ( t )(5.26) angegeben werden: x k+1 = x ( t k+1 )= g + F A /M 2 ( t k+1 - t k ) 2 + x + k ( t k+1 - t k )+ x k . (5.42) 47

KAPITEL 5. MODELLIERUNG 5.2.4 Zusammenfassung des Stoßkontaktmodells 4. Ordnung Das Modell besitzt sieben Parameter, von denen sein Verhalten abh¨angig ist. Die Systemparameter des Stoßbohrermodells 4. Ordnung sind 1 , 2 und µ = M/m welche fu¨r die Stoßzahl zwischen Stoßmasse und Oszillator, die Stoßzahl zwischen Stoßmasse und Prallplatte und fu¨r das Verh¨altnis Oszillatormasse zu Stoßmasse stehen. Die Prozessparameter des Stoßbohrermodells sind u 0 ,, F A /M und g , die fu¨r die Schwingamplitude der Oszillatorspitze, die Anregekreisfrequenz, die Beschleunigung aus axialer Anpresskraft auf den Oszillator und den Anteil der Erdbeschleunigung in L¨angsrichtung stehen. Die vier diskreten Zustandsgr¨oßen des Stoßbohrermodells 4. Ordnung sind die Stoßkontaktzeit t k der Kugel am Oszillator, die Position des Oszillatorschwerpunktes x k im Stoßzeitpunkt, die Geschwindigkeit des Oszillatorschwerpunktes x + k unmittelbar nach dem Stoß und die Geschwindigkeit der Kugel y k + unmittelbar nach dem Stoß. Ausgehend von dem Zustand x k = t k , x k , x + k , y k + T , der einen Kontakt zwischen Stoßmasse und Oszillator beschreibt, wird die Abbildung des Zustands k auf den Zustand k + 1 als implizite Funktion f ( x k , x k+1 )= 0 (5.43) angegeben, die sich zusammensetzt aus den Gleichungen fu¨r t k+1 (5.30) bzw.(5.36), x k+1 (5.42), x + k+1 (5.40) und y k ++1 (5.41). Die Berechnung des Folgezustands x k+1 gliedert sich in zwei Schritte auf. Im ersten Schritt wird der Kontaktzeitpunkt t k+1 berechnet. Im zweiten Schritt werden die drei u¨brigen Zustandsgr¨oßen x k+1 , x + k+1 , y k ++1 berechnet. Um im ersten Schritt die Kontaktzeit t k+1 zu berechnen, wird zwischen zwei m¨oglichen F¨allen unterschieden, n¨amlich dem Oszillator-PrallplatteOszillator Stoß(1. M¨oglichkeit) und dem Oszillator-Oszillator Mehrfachstoß (2. M¨oglichkeit). Um festzustellen, welcher Fall auftritt, werden zwei Zeitpunkte berechnet und miteinander verglichen. Einerseits ist dies die theoretische Kontaktzeit t W k der Kugel mit der Prallplatte(oder " W and") t W k = t k + 1 g y k + + ( y k + ) 2 + 2 g ( x k - u 0 cos t k ) , (5.44) siehe(5.28). Andererseits muss u¨berpru¨ft werden, ob ein Mehrfachstoß zwischen Kugel und Oszillator auftrat, bevor der Kugel-Prallplatte Stoß t W k stattfinden wu¨rde. Hierfu¨r muss der Zeitpunkt t M k+F1S berechnet werden aus t M k+F1S = min t M k+F1S | t M k+F1S > t k  y ( t M k+F1S )= u ( t M k+F1S ) , (5.45) siehe(5.30). Die beiden Zeitpunkte t W k und t M k+F1S werden miteinander verglichen. Der Fall, in dem ein Kugel-Prallplatte Stoß direkt auf einen KugelOszillator Stoß folgt(1. M¨oglichkeit), tritt auf, falls t W k < t M k+F1S  t M k+F1S =  . (5.46) 48

5.3. ALGORITHMUS ZUM LO¨ SEN TRANSZENDENTER STOSSGLEICHUNGEN Fu¨r diesen Fall erh¨alt man den Zeitpunkt fu¨r den n¨achsten Kugel-Oszillator Kontakt t k+1 aus t k+1 = min { t k+1 | t k+1 > t W k  y 2 ( t k+1 )= u ( t k+1 )  y 2 ( t k+1 ) > 0 }, (5.47) wobei u ( t k+1 )= x ( t k+1 ) - u 0 cos t k+1 mit x ( t k+1 )= x k+1 aus Gleichung (5.42). Im anderen Fall, wenn t M k+F1S < t W k , wenn also ein Kugel-Oszillator Mehrfachstoß auftritt, d. h. zwei aufeinanderfolgende Kugel-Oszillator Kontakte ohne einen zwischenzeitlichen Kugel-Prallplatte Stoß erfolgen, erh¨alt man den gesuchten Zeitpunkt t k+1 = t M k+F1S (5.48) aus Gleichung(5.45). Nachdem der diskrete Zeitpunkt t k+1 bekannt ist, folgt der zweite der beiden Schritte, in dem die restlichen drei der vier diskreten Zustandsgr¨oßen in expliziter Weise aus den Gleichungen(5.40)­(5.42) mit µ = M/m resultieren: x + k+1 = 1 1+ µ (1+ 1 )( y (2) ( t k+1 )+ u ¯( t k+1 ))+( µ- 1 ) x  ( t k+1 ) , y k ++1 = 1 1+ µ µ (1+ 1 )( x  ( t k+1 ) - u ¯( t k+1 )) ( µ 1 1) y (2) ( t k+1 ) , x k+1 = g + F A M 2 ( t k+1 - t k ) 2 + x + k ( t k+1 - t k )+ x k . (5.49) Als Ergebnis der gesamten Prozedur kann schließlich der Zustandsvektor x k+1 =[ t k+1 , x k+1 , x + k+1 , y k ++1 ] T aufgeschrieben werden. 5.3 Algorithmus zum L¨osen transzendenter Stoßgleichungen Eine besondere Schwierigkeit stellt das Finden der Schnittpunkte zwischen einer linearen Funktion(einer Geraden) und einer harmonischen Funktion (z. B. der Kosinusfunktion) dar. Im Fall der horizontalen(gravitationsfreien) Bewegung im Schwingstoßsystem 2. Ordnung finden keine Beschleunigungen statt. Daher bewegt sich die Stoßmasse auf einer geraden Flugbahn. Ihre gerade Flugbahn beginnt-- falls ein sog. Mehrfachstoß auftreten wird (kein Erreichen der Prallplatte vor dem folgenden Oszillator-Stoß)-auf der Position der Oszillatorspitze im Stoßzeitpunkt(Kosinusfunktion), andernfalls an der Prallplatte mit jeweils gegebener Abfluggeschwindigkeit der Kugel vom Oszillator bzw. von der Prallplatte. Wird das Modell 2. Ordnung zum Modell 4. Ordnung erweitert, indem m¨ogliche Starrk¨orperbewegungen des Oszillators zugelassen werden und Anpress- und Gravitationskr¨afte wirken k¨onnen, so haben die Flugbahnen parabelf¨ormigen Verlauf. Im Zuge der L¨osungsfindung der Bewegungsgleichungen dieses Modells mu¨ssen die Schnittpunkte zwischen einer 49

KAPITEL 5. MODELLIERUNG Kosinusfunktion und einer Parabel-- also einem Polynom 2. Ordnung -- gefunden werden, was ein ungleich schwierigeres Problem als das mit linearer Flugbahn darstellt. Problembeschreibung Bei den genannten Aufgaben existieren(Null-Geschwindigkeiten ausgeschlossen) jeweils endlich viele L¨osungen. In allen F¨allen ist es entscheidend, genau eine bestimmte L¨osung zu finden, n¨amlich diejenige, die zeitlich dem Startpunkt der Flugbahn am n¨achsten liegt, also am fru¨hesten auftritt. Sei der Startpunkt mit x s bezeichnet, so ist also die kleinste L¨osung x c mit x s < x c gesucht. Im Folgenden sollen die m¨oglichen Verfahren beschrieben werden, durch die die gesuchte L¨osung numerisch angen¨ahert werden kann. Schnittpunkte zwischen Gerade und Kosinusfunktion Problem: Gesucht ist die kleinste L¨osung x c mit x c > x s und bx + c = cos( x ) , b, c, x s R , b = 0 . Schnittpunkte zwischen Parabel und Kosinusfunktion Problem: Gesucht ist die kleinste L¨osung x c mit x c > x s und ax 2 + bx + c = cos( x ) , a, b, c, x s R , a = 0 . L¨osungsalgorithmus Zun¨achst kann auf einfache Weise bestimmt werden, ob die Parabel y = ax 2 + bx + c den Wertebereich[ 1 , 1] der Kosinusfunktion schneiden kann: Fu¨r a< 0 ist die Parabel nach unten ge¨offnet. Die Parabel kann die Kosinusfunktion h¨ochstens dann schneiden, wenn fu¨r ihren Extremwert y e = - b 2 4 a + c folgendes gilt: y e > 1. Es k¨onnen daraufhin die Stellen des Ein- und Austritts der Parabel aus dem Wertebereich[ 1 , 1] der Kosinusfunktion berechnet werden. Dadurch ergibt sich fu¨r 1 < y e < 1 ein Intervall[ x 1 , x 2 ] mit x 1 = 1 2 a - b x 2 = 1 2 a - b + b 2 4 a ( c + 1) b 2 4 a ( c + 1) (5.50) (5.51) bzw. fu¨r y e > 1 ergeben sich zwei Intervalle[ x 1 , x 2 ] und[ x 3 , x 4 ] mit x 1 = 1 2 a - b + x 2 = 1 2 a - b + x 3 = 1 2 a - b x 4 = 1 2 a - b b 2 4 a ( c + 1) b 2 4 a ( c 1) b 2 4 a ( c 1) b 2 4 a ( c + 1) (5.52) (5.53) (5.54) (5.55) 50

5.3. ALGORITHMUS ZUM LO¨ SEN TRANSZENDENTER STOSSGLEICHUNGEN in denen x c h¨ochstens liegen kann. Ferner l¨asst sich zeigen, dass, wenn eine Funktion in[ x 1 , x 2 ] den Wertebereich[ 1 , 1] durchquert, ihre kleinste Schnittstelle mit der Kosinusfunktion auch in[ x 1 , x 1 + 2 ] liegen muss. Folglich liegt die kleinste Schnittstelle in [ x 1 , x 2 ] [ x 1 , x 1 + 2 ] bzw. falls y e > 1 in[ x 1 , x 2 ] [ x 1 , x 1 + 2 ] oder in [ x 3 , x 4 ] [ x 3 , x 3 + 2 ]. Also kann die gesuchte Schnittstelle auf ein Intervall [ x li , x re ] mit Breite < 2 eingegrenzt werden, dessen Mitte an der Stelle x mi = x li + x re 2 liegt. In einem zweiten Schritt werden nun die Schnittstellen numerisch ermittelt. Hierzu wird die Problemstellung zun¨achst in ein neues Koordinatensystem mit z = x- x mi transformiert. Durch diese Koordinatentransformation l¨asst sich das Problem wie folgt formulieren: Schnittpunkte zwischen Parabel und Kosinusfunktion Problem: Gesucht sind alle L¨osungen z c der Gleichung a ( z c + x mi ) 2 + b ( z c + x mi )+ c = cos( z c + x mi )(5.56) innerhalb des symmetrisch zu z = 0 liegenden Intervalls x re 2 x li , x re 2 x li mit x re - x li 2  . Mit der trigonometrischen Umformung cos( + )= cos  · cos  sin  · sin l¨asst sich Gl.(5.56) schreiben als a ( z c + x mi ) 2 + b ( z c + x mi )+ c cos x mi · cos z c + sin x mi · sin z c = 0 . (5.57) In Gl.(5.57) ersetzen wir die beiden harmonischen Funktionen cos z c und sin z c nun durch ihre Potenzreihen. Die gesuchten L¨osungen z c liegen wie oben erl¨autert innerhalb des Radius r mit r< um z = 0. In dieser Umgebung ist der Fehler zwischen dem Taylorpolynom 8. Ordnung der harmonischen Funktionen und ihrem tats¨achlichen Wert hinreichend klein (" < S1¨a0m t 8 l"ic)h. e Nullstellen des so gebildeten Polynoms lassen sich mittels Eigenwertberechnung der Companion-Matrix der Polynomialkoeffizienten numerisch berechnen. Die kleinste reelle Nullstelle z c > z s mit z s = x s - x mi wird bestimmt und zur x Koordinate x c zuru¨cktransformiert. Zur Erh¨ohung der Genauigkeit k¨onnen ausgehend von der gefundenen L¨osung x c noch einige Schritte mit dem Newton-Verfahren durchgefu¨hrt werden. 51

.

Kapitel 6 Verhalten des Stoßbohrers In diesem Kapitel wird das Bewegungsverhalten des Stoßbohrers beschrieben. Hierzu wurden im Labor Experimente mit einem Schwingstoßpru¨fstand, der dem Modell 2. Ordnung aus Unterkapitel 5.1 entspricht, und einem Prototypen des Ultraschall-Stoßbohrers, der in Kapitel 3 vorgestellt wurde, durchgefu¨hrt. Mit verschiedenen Messverfahren werden Starrk¨orperbewegungen der Stoßmasse und des Oszillators ermittelt. Das dynamische Verhalten beider in Kapitel 5 hergeleiteten Modelle wird simuliert und dem Verhalten der realen Laborsysteme gegenu¨bergestellt. Die Unterkapitel 6.1 und 6.2 beinhalten Experimente und Simulationen zum System 2. Ordnung. Die Unterkapitel 6.3 und 6.4 beinhalten Experimente und Simulationen zum System 4. Ordnung. Unterkapitel 6.5 vergleicht die theoretischen Simulationsergebnisse der Modelle mit den praktisch gewonnenen Daten der realen Laborversuche. 6.1 Experimentelle Untersuchungen am Schwingstoßpru¨fstand Dieses Unterkapitel beschreibt die Experimente, die am System 2. Ordnung durchgefu¨hrt wurden. Nachdem der Aufbau des Pru¨fstands im Labor erl¨autert wird, werden die Messverfahren, mit denen ein Zugang zur Bewegung der Kugel m¨oglich ist, vorgestellt. Der letzte Abschnitt gibt Messdaten fu¨r die Stoßzahlen wieder. 6.1.1 Aufbau des Schwingstoßpru¨fstands Zur experimentellen Analyse der Stoßdynamik wurde ein Laborversuch aufgebaut, der sich an dem in Abschnitt 5.1 aufgestellten Modell 2. Ordnung orientiert. Bild 5.1 auf Seite 35 zeigt dieses Modell schematisch. Den Oszillator bildet im Laborexperiment ein Ultraschalltransducer, wie er links im Photo auf Seite 18 zu sehen ist. Der Stapelaktor aus vier Piezokeramikscheiben wird mit einer sinusf¨ormigen Wechselspannung gespeist, die u¨ber einen Verst¨arker von einem Sinusgenerator erzeugt wird und in Frequenz und Amplitude eingestellt werden kann. Die Anregefrequenz wird auf 53

KAPITEL 6. VERHALTEN DES STOSSBOHRERS Waveform Generator Netzteil -15 V 0 V +15 V 1.2 k SpeicherOszilloskop Verst¨arker 1 k 1 k Piezo-Oszillator HochgeschwindigkeitsKamera FastCam VLiabsreormDeoteprpler Laser Vibrometer Controller Bild 6.1: Skizze des Schwingstoßpru¨fstands die erste L¨angseigenfrequenz eingestellt, die den Oszillator in seiner ersten L¨angseigenmode anregt. Fu¨r den verwendeten Oszillator betr¨agt diese Frequenz f = 20 . 19 kHz. Die erste L¨angseigenmode ist auf Seite 41 in Bild 5.3 rechts des Oszillators gezeigt: an den Enden des Oszillatorstabes befinden sich Schwingungsb¨auche, dazwischen ein Schwingungsknoten. Die Schwingungsb¨auche lassen die Enden in L¨angsrichtung sinusf¨ormig schwingen. Im Schwingungsknoten wird der Oszillatorstab gehalten und auf den Labortisch geschraubt. Im Abstand von 11 mm in L¨angsrichtung von der Oszillatorspitze wird als eine Gegenprallelement Hammerbahn 1 auf (d"ePnraLllapblaotrttei"sc;h in den Skizzen geklemmt. Die gru¨n dargestellt) Hammerbahn und die platte Schwingfl¨ache des Oszillators stehen parallel zueinander. Im Zwischenraum wird als frei fliegende Masse eine Stahlkugel mit 6 mm Durchmesser pendelnd aufgeh¨angt(in den Skizzen rot dargestellt). Sie kann sich horizonal in Schwingungsrichtung in einem Freiraum von h = 5 mm frei bewegen. Zum Start des Experiments werden Sinusgenerator und Verst¨arker eingeschaltet. Die frei pendelnde Kugel muss leicht angestoßen werden, gerade so stark, dass sie einmal mit dem Schwinger in Kontakt kommt. 1 Die Hammerbahn bezeichnet die flache Seite eines Hammerkopfes. 54

6.1. EXPERIMENTELLE UNTERSUCHUNGEN AM SCHWINGSTOSSPRU¨ FSTAND Nach einigen wenigen transienten St¨oßen stellt sich eine Bewegung der Kugel ein, bei der sie in ihrem Freiraum hin- und herfliegend abwechselnd St¨oße am schwingenden Ende des Oszillators und an der Hammerbahn ausfu¨hrt. Bild 6.2: Stabfu¨hrung gelochter Scheibe Bild 6.3: Fu¨hrung u¨ber r¨ohrchenf¨ormigen "Balken" Bild 6.4: Fu¨hrung durch Blattfeder Um die Kugelbewegung auf eine 1-dimensionale und nahezu geradlinige Flugbahn zu begrenzen, ist eine 5-wertige Lagerung der Kugel notwendig (eine Kugelfu¨hrung, die die Kugelbewegung in 5 ihrer 6 Freiheitsgrade verhindert): diese Fu¨hrung soll die Bewegung der Kugel in L¨angsrichtung des Oszillators nicht beeinflussen, jedoch eine translatorische Bewegung entlang der zwei weiteren Raumdimensionen und eine rotatorische Bewegung um alle drei Raumrichtungen unterbinden. Die Bilder 6.2­6.4 zeigen einige realisierte Fu¨hrungsm¨oglichkeiten. Als nachteilig im Sinne einer 1-dimensionalen Bewegung haben sich die L¨angsstabfu¨hrung(Bild 6.2) und die Blattfederfu¨hrung (Bild 6.4) erwiesen. Um die Stoßmasse auf einem Stab in L¨angsrichtung gleiten zu lassen, wurde eine kreisrunde, scheibenf¨ormige Metallscheibe im Zentrum gelocht und auf einem runden Stab gefu¨hrt. Die Reibung zwischen Stab und Loch der Scheibe l¨asst diese taumeln und sich verkanten. Die Bewegung der Scheibe entspricht anstelle eines geradlinigen Fluges einer komplizierten, 3-dimensional translatorischen und rotatorischen Bewegung um ihre Raumachsen. Somit ist diese Fu¨hrungsvariante fu¨r eine Bewegungsstudie von begrenzter Komplexit¨at ungeeignet. Eine klarere Bewegung fu¨hrt die in einer Blattfeder gehaltene Kugel aus, die weit entfernt von deren statischer Gleichgewichtslage betrieben wird. Die in Bewegungsrichtung sehr biegeweich gew¨ahlte Blattfeder ist jedoch auch um ihre eigene L¨angsachse torsionsweich, sodass die Kugel eine schlackernde Bewegung ausfu¨hren kann. Eine passable Fu¨hrung der Kugel konnte mit einem r¨ohrchenf¨ormigen Balken 2 geringer Eigenmasse erzielt werden(Bild 6.3), der drehbar an einer 2 Der Balken bezeichnet in der Kontinuumsmechanik ein elastisches Element, das Querund L¨angskr¨afte sowie Biege- und Torsionsmomente u¨bertragen kann. 55

KAPITEL 6. VERHALTEN DES STOSSBOHRERS Welle gelagert ist. Nur geringfu¨gig werden die Kugelbewegungen von den Biegeschwingungen des du¨nnen Metallr¨ohrchens beeinflusst. 6.1.2 Messtechnik Um die Bewegung der Kugel zu erfassen, ist eine beru¨hrungsfreie Messtechnik notwendig, die die Bewegung der Kugel nicht beeinflusst. M¨ogliche Messmethoden sind in Tabelle 6.1 aufgefu¨hrt. Messmethode Messziel Laser-Doppler-Vibrometrie Auslenkung der Oszillatorspitze Strommessung Piezoaktor Stoßzeitpunkte am Oszillator Hochgeschwindigkeitsfilmen Kontinuierliche Position der Stoßmasse Elektrische Kontaktdetektion Stoßzeitpunkte aller St¨oße Tabelle 6.1: Verschiedene Messmethoden zur Aufnahme der Kugelbewegung Laser-Doppler-Vibrometrie Der Einsatz eines Laser-Doppler-Vibrometers 3 zur Aufnahme der Kugelbewegung gestaltet sich schwierig, weil das Messobjekt im makroskopischen 4 Bereich seine Position ¨andert. So ist es nicht m¨oglich, den Laserstrahl auf die Kugel zu fokussieren. Auch kann die Kugel nicht l¨angs ihrer Bewegungsrichtung per Laserstrahl erfasst werden, da Oszillator und Gegenprallplatte dwiuerKdeugaebl earufeiinhgreerseFtzlut,gaucmhsed"ieumulstcrhalsicehßaenll"fr.eDquaesnLtaeseSrc-hDwoipnpgluenr-gVdibersoOmseztielrlators in Frequenz und Amplitude zu erfassen. Hierfu¨r wurde es unter einem Winkel von 45 auf die schwingende Oberfl¨ache gerichtet. So wurde der Strahl nicht durch die Flugbahn der Kugel unterbrochen. Bild 6.5 zeigt das Geschwindigkeitssignal der Oszillatorspitze und best¨atigt eine Annahme, die fu¨r die Modellierungen getroffen wurde: in den Modellen wurde der Oszillator als eine kinematische Fußpunktanregung formuliert, dessen harmonische Bewegungen im Modell durch die St¨oße am Oszillator unbeeinflusst bleiben. Das Messsignal zeigt, dass die Oszillatorschwingung durch die Stoßimpulse beeinflusst wird. holt" sich die Schwingbewegung jedoch sehr schnell Nach einem Stoß und befindet sich "veorrdem n¨achsten Stoß wieder in ihrer station¨aren, freien Schwingungsmode. 3 Die Laser-Doppler-Vibrometrie nutzt den Doppler-Effekt, nach dem die von einer bewegten Oberfl¨ache reflektierte Welle eine von der urspru¨nglich ausgesandten Welle unterschiedliche Frequenz besitzt. Dadurch l¨asst sich die Geschwindigkeit einer bewegten Oberfl¨ache messen. 4 Die Bezeichnung makroskopisch besagt(im Ggs. zu mikroskopisch), dass etwas mit bloßem Auge, ohne optisch-vergr¨oßernde Hilfsmittel erkennbar ist. 56

6.1. EXPERIMENTELLE UNTERSUCHUNGEN AM SCHWINGSTOSSPRU¨ FSTAND 8 KontaktSignal [V] 0 -8 Laser-Vibrometer Signal [m/s] 2 1 0 -1 -2 0 0.05 0.1 0.15 0.2 Zeit t[s] Bild 6.5: Kontaktsignal und Laservibrometersignal 0.25 Messung der Stromaufnahme Die Messung des Stromsignals ist die mit Abstand einfachste und somit sicherste M¨oglichkeit, auf indirekte Weise Auskunft u¨ber die ungef¨ahren Zeitpunkte der Kugel-Oszillator-St¨oße zu erhalten. Lediglich eine Strommesszange wird verwendet, um den Stromfluss vom Verst¨arker zum Piezotransducer zu messen. Das Stromsignal l¨asst sich mit einem Speicheroszilloskop aufzeichnen. Die indirekte Ermittlung der Stoßzeitpunkte ist m¨oglich, da die Stromaufnahme kurzzeitig ansteigt oder abf¨allt, wenn die Stoßwelle, die die Stoßmasse in den Transducerstab einbringt, die freie Resonanzschwingung des Piezotransducers u¨berlagert. Die Laufzeit der Stoßwelle durch den 60 mm langen Transducerstab(ca. 12 µ s) entspricht etwa einer Viertelschwingungsperiode und ist fu¨r eine grobe Zeitpunktsbestimmung daher weitgehend vernachl¨assigbar. Den exakten Zeitpunkt des Stoßes nur anhand des Stromsignales festzustellen ist jedoch schwierig, da die Abweichung des Stromsignals vom bisherigen harmonischen Verlauf bewertet werden muss. Die Messmethode der Stromaufnahme trifft keine Aussage u¨ber die Bewegung der Kugel zwischen zwei St¨oßen am Oszillator. Bild 6.6 zeigt das Messergebnis des Stromsignals zusammen mit dem u¨berlagerten Kontaktmesssignal. Einsatz einer Hochgeschwindigkeitskamera Aussagen u¨ber die Kugelbewegung zwischen ihren St¨oßen liefert allein die Hochgeschwindigkeitskamera, da sie die Kugelposition in direkter Weise optisch zu jedem Zeitpunkt erfasst. Ein weiterer, großer Vorteil dieser Messmethode ist, dass s¨amtliche makroskopischen Bewegungen zeitlich fein aufgel¨ost 57

Kontaktsignal in [8V] Stromverlauf in [6.4A] KAPITEL 6. VERHALTEN DES STOSSBOHRERS 1 0 -1 0.145 0.15 0.155 0.16 0.165 0.17 0.175 0.18 Zeit t[s] Bild 6.6: Kontaktsignal und Stromsignal: Deutlich zu sehen sind die Stromeinbru¨che und-spitzen, die durch die St¨oße am Oszillator hervorgerufen werden; vor dem jeweils folgenden Stoß befindet sich die Stromaufnahme wieder im harmonisch-station¨aren Zustand. wiedergegeben werden k¨onnen. Problemlos lassen sich die Schwingungen der umgebenden Struktur qualitativ analysieren und ihr Einfluss auf die Dynamik der Kugel beurteilen. So l¨asst sich schnell abw¨agen, welche Art der Kugelfu¨hrung inwieweit nachteilig erscheint. Auch zur quantitativen Analyse der Kugelbewegung k¨onnen fein aufl¨osende Bilder der Hochgeschwindigkeitskamera herangezogen werden: ein Formerkennungsalgorithmus in der Software Photron Motion Tools wurde genutzt, um die Position der Kugel auf jedem aufgenommenen Einzelbild zu quantifizieren. So entsteht eine genaue Zeitreihe der Kugelbewegung. Da pro Zeitschritt ein Bild erzeugt und zur sp¨ateren Auswertung digital gespeichert wird, ist der Speicherbedarf fu¨r diese Methode hoch und erm¨oglicht vergleichsweise kleine Aufnahmedauern. Ein Einzelbild des Hochgeschwindigkeitsfilms der Kugelbewegung zeigt Bild 6.7. Die kleinen roten + - Zeichen markieren den Verlauf des oberen Kugelrandes auf den ersten 63 Einzelbildern der Aufnahme, was bei einer Aufnahmerate von 5000 fps(frames per second) einer Aufnahmedauer von 12.4 ms entspricht. Die Wellenform der Spur ist ein Hinweis auf eine nichtgeradlinige oder leicht bogenf¨ormige Bewegung, wie man es von einem masselosen und unendlich steifen Fu¨hrungsbalken erwartet h¨atte. Statt dessen u¨bt die Biegeschwingung des Metallr¨ohrchens, das die Kugel u¨ber eine 80 mm entfernte Drehachse fu¨hrt, einen deutlichen Einfluss auf die Kugelbewegung aus. Die Positionsbestimmung der Stoßmasse in allen aufeinanderfolgenden Bildern der Hochgeschwindigkeitsaufnahme fu¨hrt zu einem kontinuierlichen Bewegungsablauf, der in Bild 6.8 u¨ber der Zeit aufgetragen ist. Die Auswertung von Hochgeschwindigkeitsdaten ist fu¨r einen derartigen Pru¨fstand die einzige M¨oglichkeit, beru¨hrungsfrei makroskopische Positionen aufzunehmen. Daru¨berhinaus ist es eine direkte Methode zur Positionsbestimmung. Die anderen hier genannten Messverfahren lassen nur indirekt auf die Bewegung der Stoßmasse schließen. Hierfu¨r war bisher die Annahme Voraussetzung, dass die Stoßmasse zwischen den St¨oßen eine unbeschleunigte Bewegung ausfu¨hrt. Die Analyse der Kameradaten in Bild 6.8 validiert nun diese Annahme. 58

Position [mm] Flugmasse 6.1. EXPERIMENTELLE UNTERSUCHUNGEN AM SCHWINGSTOSSPRU¨ FSTAND Bild 6.7: Einzelbild der Hochgeschwindigkeitskamera mit Spur(+++) der Positionsidentifizierung der Kugel w¨ahrend 12.4 ms; am unteren Bildrand ist die Oszillatorspitze zu sehen, am oberen Rand die Hammerbahn. ------------- Wand(Hammer-Bahn)------------3 2 1 0 ------------- Oszillator------------3.2 3.3 3.4 3.5 3.6 3.7 3.8 Zeit t[s] Bild 6.8: Position u¨ber Zeit der Stoßmasse aus der Positionserkennung der Hochgeschwindigkeitsaufnahme Elektrische Kontaktmessung Bei der elektrischen Kontaktmessung wird die Leitf¨ahigkeit der metallischen Stoßpartner ausgenutzt, um durch Schließen und O¨ ffnen elektrischer Kreise den Kontaktzustand zu u¨bermitteln. Hierfu¨r mu¨ssen der Oszillator, die Kugel sowie der Hammer elektrisch voneinander getrennt auf den Labortisch montiert werden. Dies l¨asst sich erreichen, indem man die Bauteile z. B. zwischen Pertinaxplatten einspannt oder sie mit Magnetfu¨ßen ausstattet, die u¨ber Kunststofffolien vom metallenen Labor-Spanntisch isoliert werden. Zu beachten ist, dass der Transducerstab des Oszillators u¨ber die Wechselspannungserregung vom Verst¨arker immer geerdet ist. Die u¨brigen zwei Teile mu¨ssen daher vom ebenfalls geerdeten Labortisch elektrisch isoliert werden. Die Skizze 6.1 zeigt die elektrische Verschaltung der Bauteile. Dem Hammerkopf als Gegenprallfl¨ache wird aus einem Gleichspannungsnetzteil ein Potential von+15 V, dem Oszillator ein Potential von-15 V auferlegt. Eine feine, auf die Kugel gel¨otete Einzellitze 5 und der 0 V-Ausgang des Netzteils werden an einen Messkanal eines Speicheroszilloskops angeschlossen. Besteht w¨ahrend eines Stoßes elektrisch leitender Kontakt zwischen Kugel und Oszillator, so wird ein Stromkreis geschlossen, der aus einer Reihenschaltung aus 5 Eine Litze ist ein sehr du¨nnes Metalldr¨ahtchen. Flexible, elektrische Leitungen bestehen aus einem Verbund vieler Metalllitze. 59

KAPITEL 6. VERHALTEN DES STOSSBOHRERS der 15 V­Spannungsquelle, einem 1 k­ und einem 1.2 k­Widerstand besteht. Der am gr¨oßeren Widerstand entstehende Spannungsabfall von 8.3 V wird von einem potentialgetrennten Eingang des Oszilloskops aufgezeichnet. Die 1 k­Widerst¨ande dienen lediglich dem Vermeiden von Kurzschlu¨ssen w¨ahrend Einstellarbeiten am Versuchsaufbau. Der 1.2 k­Widerstand beschleunigt das Entladen von oszilloskopseitigen Kapazit¨aten nach Trennung der Kontakte. W¨ahrend einer Messung wird der Spannungsverlauf dieses Messkanals im Oszilloskop mit einer hohen Abtastrate 6 von 5 MS/s aufgezeichnet, aus dem sich mit hoher zeitlicher Genauigkeit auf die Stoßzeitpunkte und die Kontaktdauern schließen l¨asst. Ein derart entstandener Kontaktspannungsverlauf ist in roter Farbe in Bild 6.5 und Bild 6.6 zu sehen; Bild 6.9 zeigt den Kontaktspannungsverlauf eines einzelnen Stoßes zusammen mit dem Auslenkungssignal der Oszillatorspitze, das aus Integrieren des Laservibrometersignals gewonnen wurde. Rot gezeichnet ist in Bild 6.9 die Zeitspanne, w¨ahrend der metallischer Kontakt zwischen beiden Stoßpartnern besteht. Beginn und Ende des Kontaktes lassen sich aus der Kontaktspannung ableiten. Deutlich sieht man den exponentiellen Anstieg und das Abklingen der Spannung, bedingt durch Induktivit¨aten und Kapazit¨aten in Kabeln und Oszilloskop. Als Kontaktbeginn und ­ende werden jeweils die Startzeitpunkte der Spannungsanstiege bzw. abf¨alle ausgewertet. 0 -2 -4 -6 -8 -8 0 11 20 30 40 50 -6 -3 0 3 6 8 0 11 20 30 40 50 Zeit t[ µ s] Bild 6.9: Messdaten w¨ahrend eines Stoßes(rot) am Oszillator: Kontaktspannungsverlauf(oben) und Auslenkung der Oszillatorspitze(unten) 6 Die Abtastfrequenz oder engl. Samplingrate wird in der Einheit samples per second (S/s) angegeben. 60 Auslenkung Oszillator [ µ m] Kontaktsignal [V]

6.1. EXPERIMENTELLE UNTERSUCHUNGEN AM SCHWINGSTOSSPRU¨ FSTAND Paralleler Einsatz verschiedener Messmethoden Die unterschiedlichen hier vorgestellten Messmethoden schließen sich gegenseitig nicht aus, sondern sind synchron nebeneinander einsetzbar. Allein die spannungsliefernden Messverfahren machen sich gegenseitig Konkurrenz bei der hochfrequenten zeitlichen Abtastung der Spannungssignale in den Oszilloskopeing¨angen. Bei einer hohen Abtastrate k¨onnen nicht alle vier Messsignale aufgezeichnet werden. Die Bilddaten der Kamera werden in einem kamerainternen 40 GB großen RAM(Speicher) erfasst. Zur Synchronisation aller Aufzeichnungen ist es jedoch n¨otig, auch einen Triggerkanal mit aufzuzeichnen, der so einen weiteren Datenfluss erzeugt. Die zeitsynchrone Auswertung der Messdatenverl¨aufe ergab, dass sich die unterschiedlichen Methoden gegenseitig verifizieren k¨onnen. Die Hochgeschwindigkeitsfilmaufnahmen haben gezeigt, dass die Annahme nichtbeschleunigter Kugelbewegungen gerechtfertigt ist. 6.1.3 Stoßzahlen aus Labormessungen Der Parametersatz des hier untersuchten Stoßschwingers 2. Ordnung umfasst die zwei Prozessparameter u 0 und sowie die drei Systemparameter h , 1 und 2 . Die Stoßzahlen 1 und 2 beschreiben die D¨ampfungseffekte w¨ahrend der St¨oße. Es ist nicht m¨oglich, den Wert der Stoßzahlen allein aus Kru¨mmungsradien der Oberfl¨achen, Oberfl¨achenbeschaffenheit, Materialzusammensetzung, H¨arte und Temperatur der Stoßpartner zu berechnen, vielmehr ist er in der Praxis auch von Stoß zu Stoß leicht schwankend. Aus einer Serie von Messreihen, die fu¨r die Flugweiten h = 3 mm und h = 5 mm unter nach M¨oglichkeit ansonsten konstanten Bedingungen durchgefu¨hrt wurden, wurden Abweichung und Mittelwert der Stoßzahl 2 fu¨r den Stoß zwischen der Hammerbahn frei fliegenden ("Prallplatte") Kexupgeerlim("eSnttoeßllmbaessstei"m)mutn. d der eingespannten Da fu¨r die Stoßzahl der Zusammenhang(5.4) (Seite 37) gilt, waren die Geschwindigkeiten vor und nach jedem Stoß zu bestimmen. Die Kugelgeschwindigkeiten sind nicht direkt messbar. Sie wurden indirekt aus den Laufzeiten zwischen dem Oszillatorstoß und dem Prallplattenstoß bestimmt. Bild 6.10 gibt im oberen Diagramm die Verteilung der Stoßzahlen aus 520 Messungen in einem Histogramm wieder. Das Maximum der H¨aufigkeitsverteilung liegt bei 2, max = 0 . 84, der Mittelwert bei 2, Mittelw. = 0 . 81. Physikalisch m¨oglich und sinnvoll sind nur Stoßzahlen im Intervall[0 , 1]. 4% der Messdaten liegen außerhalb dieses Intervalls im grau unterlegten Bereich. Eine Erkl¨arung fu¨r diese Messwerte liegt in der Natur des Messverfahrens. Liegen Stoßzahlen oberhalb von Eins, so l¨asst sich mit Gewissheit ein Laufzeitenverh¨altnis gr¨oßer Eins ableiten. Das heisst, die Kugel brauchte fu¨r den Hinweg zum Hammer weniger Zeit als fu¨r den Weg zuru¨ck zum Oszillator. Das Laufzeitenverh¨altnis mit dem Stoßgeschwindigkeitenverh¨altnis gleichzusetzen setzt aber eine konstante Geschwindigkeit der Kugel auf ihrem Weg vom Oszillator zum Hammer und zuru¨ck voraus. Eine Zeitlupenanaly61

KAPITEL 6. VERHALTEN DES STOSSBOHRERS se mittels Hochgeschwindigkeitskameraaufnahme enthu¨llte aber, dass diese Voraussetzung nicht ausnahmslos Gu¨ltigkeit beansprucht. Die Biegeschwingungen der r¨ohrchenf¨ormigen Aufh¨angung der Kugel u¨berlagern sich zur Flugbewegung der Kugel. Dieser Effekt ist in Bild 6.8 kaum ersichtlich, ist in der 2-Dimensionalit¨at der roten Positionsspur in Bild 6.7 jedoch deutlich erkennbar. Bild 6.10 zeigt im unteren Diagramm die gleichen Messdaten wie im Histogramm daru¨ber. Hier ist zus¨atzlich die Information enthalten, bei welchen Geschwindigkeiten und bei welchen Flugabst¨anden die Messdaten erzeugt wurden. Untersuchungen in der Literatur ist zu entnehmen, dass Häufigkeit 80 70 60 50 40 30 20 10 0 0 0.2 0.4 0.6 0.8 1 gemessene Stoßzahl 2 1.2 +: h= 3mm +: h= 5mm 1 1.2 1.4 Ankunftsgeschwindigkeit an Hammerbahn [m/s] 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 gemessene Stoßzahl 2 Bild 6.10: Experimentelle Messung der Stoßzahl 2 an der Hammerbahn 62

6.2. SIMULATION DES MODELLS 2. ORDNUNG Stoßzahlen auch von der Stoßgeschwindigkeit abh¨angig sind. Ein Zusammenhang zwischen Stoßzahl und Geschwindigkeit, u¨ber den unter Miteinbeziehung von Wellenausbreitungsvorg¨angen im stabf¨ormigen Stoßpartner in[Seifried, Schiehlen, Eberhard 2004] berichtet wurde, kann in den erhaltenen Messwerten nicht beobachtet werden. Die Modelle, die in Kapitel 5 aufgestellt wurden, beschreiben die kinematischen Zusammenh¨ange in den Stoßzeitpunkten durch die Stoßzahlenhypothese nach Newton. Hier wird eine zeitlich konstante Stoßzahl angenommen, die sich fu¨r einen Parametersatz nicht ¨andert. Dieses in der Literatur g¨angige Vorgehen(siehe Unterkapitel 4.3) vernachl¨assigt die in der Praxis beobachtete Streuung der Stoßzahlenwerte. Abschnitt 6.2.3 wird die qualitative Auswirkung einer Stoßzahlenstreuung auf den L¨osungstypus betrachten. 6.2 Simulation des Modells 2. Ordnung Dieses Unterkapitel betrachtet das Verhalten des Modells 2. Ordnung. Wir werden Fixpunkte analytisch bestimmen und Zeitverl¨aufe fu¨r sowohl periodische als auch chaotische L¨osungen sehen. Um die Parameterabh¨angigkeit von periodischen L¨osungen zu verstehen, werden wir eine Bifurkationsanalyse fu¨r die Systemparameter Flugabstand h sowie die Stoßzahl 1 durchfu¨hren. 6.2.1 Analytische Fixpunktbestimmung Die Analyseergebnisse in Kapitel 7.1 werden zeigen, dass fu¨r einen effizienten Betrieb des Stoßbohrsystems ein periodisches Verhalten der Stoßmasse angestrebt wird. Ein Fixpunkt der Abbildung, die das Modell 2. Ordnung definiert, stellt eine periodische Bewegung der Stoßmasse dar, denn bei jedem Stoß wird das System wiederholt den gleichen Zustand annehmen. Dies motiviert die Bestimmung von Fixpunkten. Aufgrund der impliziten Natur der Abbildungsgleichung(5.9) ist es nicht m¨oglich, Fixpunkte von mehrperiodischen Bewegungen analytisch zu bestimmen. Dies geschieht auf numerische Weise in Abschnitt 6.2.3 mittels einer Bifurkationsanalyse. Im Folgenden wird gezeigt, dass Fixpunkte [ t , V ] T zu 1-periodischen L¨osungen analytisch gefunden werden k¨onnen. Eine Anmerkung sei vorab zur Angabe der ersten Komponente eines Zustandswertes gegeben: bisher war die erste Zustandskomponente die Stoßzeit t k . Im Folgenden wird bei der Angabe eines Zustands oft anstelle der absoluten Stoßzeit t k deren relative Phasenlage k in Bezug zur harmonischen Anregung durch den Oszillator angegeben. Aus t k wird somit k durch die Transformation k =( t k mod T ) · . (6.1) Ein Zustand[ , V ] T heißt genau dann Fixpunkt der Abbildung T , wenn T [ , V ] T =[ , V ] T . 63

KAPITEL 6. VERHALTEN DES STOSSBOHRERS Fu¨r eine 1-periodische Bewegung folgt hieraus die Periodizit¨atsbedingung k+1 = k  V k +1 = V k  k N . Wegen t k der Transformation(6.1) folgt fu¨r die Zustandsgr¨oße "Stoßzeitpunkt" t k+1 = t k + n T, n N , weil der Zeitverzug eines nachfolgenden Stoßes dem Vielfachen einer Anregeperiodendauer T = 2 gleichen muss. Damit gilt fu¨r die Auslenkung u ( t )= u 0 cos t der Oszillatorspitze u ( t k )= u ( t k+1 ) . Unter Ausnutzung diesen Zusammenhangs folgt aus der Differenzengleichung(5.6) des Stoßzeitpunktes die Flugdauer t = t k+1 - t k zwischen zwei aufeinanderfolgenden St¨oßen t = 1+ 1 2 h- u 0 cos t V = nT. (6.2) Mit der Differenzengleichung(5.7) ergibt sich fu¨r die periodische Abfluggeschwindigkeit die Bedingung V = 1 2 V (1+ 1 ) u 0  sin t . (6.3) Umstellen von Gleichung(6.2) und Gleichung(6.3) nach V und Gleichsetzen fu¨hrt auf h u 0 = cos t - b sin t (6.4) oder 1= a sin t arctan b 1 (6.5) mit a = u 0 h 1+ b 2 b = 2 n (1 1+ 1 1 2 )(1+ 1 2 ) , Aus Gleichung(6.5) folgt 1, 2 (0 , 1) . t = arcsin 1 + arctan 1 ab und aus Gleichung(6.3) kann V berechnet werden: (6.6) V = 1+ 1 1 - 1 2 u 0  sin t . (6.7) Die Fixpunkte[ t , V ] T sind gefunden. Aufgrund der nicht-expliziten Formulierbarkeit der Abbildungsgleichungen ist es nicht m¨oglich, die Stabilit¨at der Fixpunkte analytisch zu bestimmen. In den folgenden Abschnitten wird sich mittels numerischer Zeitreihenintegration zeigen, dass fu¨r bestimmte Parameterbereiche einer der Fixpunkte stabil ist. 64

6.2. SIMULATION DES MODELLS 2. ORDNUNG 6.2.2 Periodische und chaotische L¨osungen In Abschnitt 5.1 wurde ein Modell 2. Ordnung aufgestellt, welches die Dynamik des zentralen Elements des Ultraschall-Stoßbohrers simulieren soll. Genauer gesagt interessieren die Bewegungen der Kugel, die durch St¨oße des schwingenden Oszillators angeregt und durch Kontakt mit einer Prallplatte in ihrem Aktionsradius begrenzt wird. Um einen ersten Eindruck von der Bewegung der Kugel zu erhalten, wollen wir ihre Dynamik in diesem Abschnitt simulieren. Das Modell 2. Ordnung wurde mit dem Differenzengleichungssystem (5.9) auf Seite 38 definiert. Der Zustandsvektor[ t k+1 , V k+1 ] T wurde dort iterativ aus einer Abbildung t k+1 V k+1 = T t k V k (6.8) erzeugt, wobei die Funktion T : R 2 R 2 nicht explizit definierbar ist. Um diese zeitdiskreten Bewegungsgleichungen u¨ber berechnet man eine Folge aufeinanderfolgender Zduesrt¨aZnedite.zDu ie"sinetFegorlgieerevno"n, Zustandspunkten ergibt einen sog. Orbit bzw. eine sog. Trajektorie. Man beginnt die Iterationen bei einem Anfangszustand[ t 0 , V 0 ] T bzw.[ 0 , V 0 ] T und iteriert diesen Zustand u¨ber die Abbildung(6.8). Wir w¨ahlen einen Parametersatz 7 Oszillator­ Ampl. Kreisfrequenz u 0 10 µ m 2  · 20 . 19 kHz Stoßzahl an Osz. Bohrer 1 2 0.6121 0.7879 Flugabstand h 5 mm (6.9) und iterieren eine Zeitreihe ausgehend von einer zuf¨allig gew¨ahlten Anfangsbedingung, etwa 0 V 0 = 0 . 265518 rad 17 . 349077 m/s . In Bild 6.11, obere H¨alfte, ist der so entstandene Bewegungsverlauf gezeichnet. Die iterierten Zustandswerte geben die Stoßzeiten t k bzw. die Abfluggeschwindigkeiten V k der Stoßmasse am Oszillator wieder, der in Bild 6.11 am unteren Bildrand in blauer Farbe dargestellt ist. Die kontinuierliche Flugbahn wurde unter Annahme beschleunigungsfreier Bewegung, also konstanter Geschwindigkeit aus der Geschwindigkeit V k fortgesetzt bis zum Erreichen der Prallplatte(gru¨ne Farbe, oberer Bildrand) im Abstand von h = 5 mm. Mit der Stoßgleichung(5.4) wird die Abfluggeschwindigkeit U k+1 7 Flugabstand, Anregefrequenz und Anregeamplitude wurden entsprechend den in Laborversuchen eingestellten Werten gew¨ahlt. Die Stoßzahlen liegen im Bereich der Werte, die im vorigen Kapitel experimentell ermittelt wurden und somit fu¨r den Pru¨fstand realistisch sind. Die hohe Genauigkeit, mit der die Stoßzahlen und die Anfangsbedingungen hier angegeben wurden, ist fu¨r die Praxis wenig sinnvoll, fu¨r die rein mathematische Untersuchung und die Erzeugung von Periodizit¨at hier jedoch zwingend notwendig und gilt auch fu¨r die u¨brigen Parameter. 65

KAPITEL 6. VERHALTEN DES STOSSBOHRERS Position Position h 0 0 0.05 0.1 0.15 h 0 0 0.05 0.1 0.15 Zeit t[s] Bild 6.11: Zeitreihe aus dem Modell 2. Ordnung; oben: chaotische Bewegung; unten: periodische Bewegung -mit chaotischer Bewegung -- zum Vergleich(Das rasche Auseinanderstreben der zwei Zeitverl¨aufe mit nahezu identischem Anfang, d. h. das Entstehen der chaotischen L¨osung ist erkennbar.) nach dem Stoß an der Prallplatte berechnet und ebenfalls als Fluggerade in den Zeitverlauf eingezeichnet. Aus dem Zeitverlauf l¨asst sich eine erste Vorstellung der simulierten Bewegung der Kugel fu¨r einen kleinen Zeitraum und fu¨r diesen einen Parametersatz und die speziellen Anfangsbedingungen gewinnen. Eine Regelm¨aßigkeit der Bewegung l¨asst sich nicht erkennen. Da sich diese Behauptung durch bloßes Betrachten eines Zeitverlaufes noch nicht mit Sicherheit aufstellen l¨asst, werden wir im n¨achsten Abschnitt in einer Bifurkationsanalyse eine andere Sichtweise auf die Art der Bewegung bekommen. Vorher wollen wir mit dem gleichen Parametersatz eine weitere Zeitreihe bilden, die von anderen Anfangsbedindungen ausgeht. Zur Wahl der Anfangswerte u¨bernehmen wir die obige Startbedingung und erh¨ohen lediglich die Anfangsgeschwindigkeit in der letzten bezifferten Stelle: 0 V 0 = 0 . 265518 rad 17 . 349078 m/s . Das Resultat der Bewegung, ausgehend von den vorgegebenen Bedingungen, ist in Bild 6.11, untere H¨alfte, dargestellt. Augenscheinlich stellt sich nach zwei Oszillatorst¨oßen eine ein-periodische Bewegung ein 8 . Tats¨achlich kann man zeigen, dass die Fixpunktbedingung T V = V von V = 1 . 49542278921781 rad 3 . 94968938856343 m/s erfu¨llt wird. Die Periodendauer der Stoßbewegung ist hier 2.87 ms bzw. entspricht einer Frequenz von 1/2.87 ms= 348 Oszillatorst¨oßen pro Sekunde. 8 Die Startwerte zu der periodischen L¨osung mit einer-- verglichen mit der Fixpunktgeschwindigkeit-- hohen Geschwindigkeit von 17 m/s wurden fu¨r dieses Beispiel ermittelt, indem ein L¨osungsalgorithmus fu¨r die inverse Abbildungsvorschrift T 1 aufgestellt und programmiert wurde. Ausgehend von dem zu diesem Zeitpunkt bereits bekannten Fixpunkt als Startbedingung wurden auf der instabilen Mannigfaltigkeit von T einige Schritte ru¨ckw¨arts iteriert. 66

6.2. SIMULATION DES MODELLS 2. ORDNUNG Plausibel erscheint dieses Ergebnis im Vergleich mit Untersuchungen des Ultraschall-Stoßbohrers am Jet Propulsion Laboratory in Pasadena: die Stoßfrequenz liegt im Bereich, der dort festgestellt wurde, siehe [Bar-Cohen et al. 2001a]. Bei exakt gleich gehaltenen Parametern fu¨hren zwei sich minimal unterscheidende Anfangsbedingungen nach kurzer Zeit auf zwei voneinander verschiedene L¨osungen. Diese fu¨r deterministisch chaotische Systeme typische Eigenschaft weist auf mindestens zwei koexistierende L¨osungen hin, die nebeneinander existieren[Strogatz 1994]. Jede der beiden m¨oglichen L¨osungsformen besitzt ihr eigenes Einzugsgebiet von Anfangsbedingungen. Wie sich die Einzugsgebiete u¨ber den Zustandsraum erstrecken, werden wir in Kapitel 7.1 untersuchen. Die rasche Separation der beiden dicht beieinander beginnenden Trajektorien aus Bild 6.11 verdeutlicht Bild 6.12, in dem die periodische L¨osung (durchgezogen --) und die chaotische(gestrichelt---) in das selbe Achsensystem u¨bereinander gelegt wurden. Der erste Vergr¨oßerungsausschnitt a hebt den dritten Stoß, der bei einer Phasenlage des Oszillators von etwa k =1.355 rad stattfindet, hervor. Die periodische L¨osung erreicht den Oszillator wenig sp¨ater, was der Oszillatormaximalgeschwindigkeit (bei Phase 1 . 5 ) geringfu¨gig n¨aher liegt. Daher ist der Impuls auf die Kugel, die sich auf der chaotischen Flugbahn bewegt, geringer. So erreicht sie den Oszillator im n¨achsten Zyklus bereits deutlich sp¨ater als die periodische Flugbahn(Ausschnitt b) und hinkt der Periodizit¨at einen weiteren Stoß sp¨ater bereits eine dreiviertel Phase hinterher(Ausschnitt c). Die Oszillatorgeschwindigkeit ist in diesem Stoßmoment aber nahezu Null, genau genommen sogar negativ, wodurch die Abfluggeschwindigkeit sogar leicht gebremst wird. Dies verursacht die lange sechste Flugdauer und begru¨ndet die Irregularit¨at der weiteren Bewegung, siehe Bild 6.11 oben. x 10 -3 4 Position [m] 2 Position [m] 0 0 x 10 -6 1 0 -1 0.005 x 10 -5 0.01 0.015x 10 -5 1 0 a -1 2 1 0 b -1 0.02 Zeit t[s] c Bild 6.12: Transiente Phase der periodischen Zeitreihe -mit u¨berlagerter chaotischer L¨osung -- ; darunter drei Detailvergr¨oßerungen. 67

KAPITEL 6. VERHALTEN DES STOSSBOHRERS 6.2.3 Bifurkationsanalyse Bevor wir die Dynamik des Modells 2. Ordnung detailliert unter Anwendung der mengenorientierten Methoden fu¨r einen ausgew¨ahlten Parametersatz im folgenden Kapitel 7 analysieren und interpretieren werden, soll dieser Abschnitt einen U¨ berblick u¨ber die Abh¨angigkeit der Dynamik von den Parametern des Modells geben. Wie wir sp¨ater genauer sehen werden, ist die am h¨aufigsten anzutreffende Form der Bewegung eine irregul¨are, d. h. chaotische. Unter gewissen Parameterkonstellationen koexistieren auch chaotische und periodische Bewegungen von unterschiedlicher Periodizit¨atszahl. Im vorigen Abschnitt wurde ein solcher Parametersatz(siehe (6.9)) gew¨ahlt und die sensible Abh¨angigkeit des L¨osungstypus von den Anfangsbedingungen anhand eines Beispiels gezeigt. Wie verh¨alt sich das System aber, wenn sich Prozess- oder Systemparameter ¨andern? In welchen Parameterbereichen existieren periodische L¨osungen neben chaotischen? Gibt es Parameter, fu¨r die periodisches Verhalten nicht m¨oglich ist? Wie verschieben sich die Fixpunkte periodischer L¨osungen, wenn sich ein Parameter ver¨andert? Diesen Fragen soll in diesem Abschnitt auf numerischem Wege nachgegangen werden. 3.95 3.9498 1.5 1.499 1.498 3.9496 1.497 1.496 V [m/s] [ rad] 3.9494 1.495 3.9492 1.494 1.493 3.949 1.492 1.491 3.9488 4.9995 5 h[m] 5.0005 x 10 -3 1.49 4.9995 5 h[m] 5.0005 x 10 -3 Bild 6.13: Bifurkationsdiagramm mit Bifurkationsparameter h : erkennbar ist die Kaskade von Pitchfork-Bifurkationen fu¨r minimal abnehmendes h . Der Ausschnitt zeigt einen Bereich von nur ca. 1 µ m Breite! Die Verzweigungs­ oder Bifurkationsanalyse entschlu¨sselt die Existenz von periodischer versus chaotischer Dynamik in Abh¨angigkeit eines Prozess- oder Systemparameters. Man erstellt ein sogenanntes Bifurkationsdiagramm, in welches man eine(beliebig) gew¨ahlte Zustandsvariable u¨ber 68

6.2. SIMULATION DES MODELLS 2. ORDNUNG den interessierenden Parameter auftr¨agt, dessen Abh¨angigkeit untersucht werden soll(siehe auch Abschnitt 4.1.6). Die Wahl der Zustandsvariablen ist weitestgehend frei. Die Abszisse 9 des Diagramms wird in feine ParameterSchritte unterteilt. Fu¨r jeden Parameterwert wird daraufhin die Trakjektorie von einem Anfangszustand aus integriert. Nachdem transiente Einschwingvorg¨ange abgeklungen sind, wertet man den Systemzustand x wiederholt an einzelnen Zeitpunkten aus, die im Abstand einer Periodendauer T aufeinander folgen, wobei t 0 das gesch¨atzte Ende der Einschwingzeit kennzeichnet: x ( t 0 ) , x ( t 0 + T ) , x ( t 0 + 2 T ) , x ( t 0 + 3 T ) ,.... Etwa 100 Zustandspunkte ermittelt man so. Die Werte der gew¨ahlten Zustandsvariablen aus jedem der Zustandspunkte plottet man als einzelne Punkte in das Diagramm u¨ber die Stelle des Systemparameters. Falls sich das System bei diesem Systemparameter 1-periodisch verh¨alt, so wird die Zustandsvariable nach jeder Periode T wieder den gleichen Wert annehmen. Im Bifurkationsdiagramm entsteht nur ein Punkt. Liegt eine 2periodische L¨osung vor, so werden zwei diskrete Punkte entstehen, da das System erst wieder nach stets 2 Periodendauern T zum ersten Zustandspunkt zuru¨ckkehrt. Findet dagegen chaotisches Verhalten statt, so liegen die Punkte verstreut 10 . Nachdem fu¨r den festgehaltenen Parameter die etwa 100 Zustandspunkte geplottet wurden, ru¨ckt man ihn ein kleines Inkrement oder Dekrement 11 weiter und wiederholt den Vorgang. Bei zeitkontinuierlichen, nichtautonomen 12 Systemen bietet sich bei periodischer Anregung als Zeitmaß T die Wahl der Anregeperiode an, zu dem die integrierte Trajektorie wiederholt ausgewertet wird. Im Falle des zeitdiskreten Modells liegt bereits eine Abbildungsvorschrift vor, die von einem Zustandspunkt in einem diskreten Schritt zum n¨achsten Punkt abbildet. Diese Schrittweite wird als ein Periodenzeitschritt( T ) angesehen. Genauso wie beim zeitkontinuierlichen System lassen sich so 1-periodische von 2-periodischen von mehr-periodischen von chaotischen ("unendlich-periodischen") Attraktoren unterscheiden. Bedingt durch die komplizierte Dynamik, die nichtglatte Stoßsysteme aufweisen, koexistieren oft mehrere Attraktoren nebeneinander. Die Analyse des vorliegenden Systems zeigte, dass fu¨r alle Parameterkombinationen stets ein chaotischer Attraktor existiert. Dies erzeugt eine Schwierigkeit bei 9 Abszisse: (lateinisch: die abgeschnittene Linie) horizontale Achse eines Koordinatensystems 10 In einer Poincar´e-Abbildung, in der sehr viele solcher Zustandspunkte im Zeitabstand T stroboskop¨ahnlich in ein Zustandsdiagramm eingezeichnet werden, erkennt man eine seltsame Struktur, die den Determinismus, also die Vorherbestimmheit, des nur oberfl¨achlich zuf¨allig traktor". wirkenden Chaos nachweist. Diese Struktur nennt man "seltsamen At11 Das Inkrement: (lateinisch Dekrement: Verminderung. "Zuwachs") Betrag, um den eine Gr¨oße zunimmt. Das 12 Ein nichtautonomes Differentialgleichungssystem verh¨alt sich von der Zeit abh¨angig, meist durch einen zeitabh¨angigen Erregerterm auf der rechten Seite. Beispiel: rechte Seite = u 0 cos t. 69

KAPITEL 6. VERHALTEN DES STOSSBOHRERS V [m/s] V [m/s] 3.9496 3.9494 3.9492 6 3.949 3 1.49 2 7 54 1 1.495 [ rad] 1.5 3.9496 3.9494 5 3.9492 10 6 2 7 11 9 4 1 3.949 8 3 1.49 1.495 [ rad] 1.5 Bild 6.14: Trajektorien fu¨r mehrperiodische L¨osungen; links: 2-periodisch ( h = 4 . 9996 mm), rechts: 4-periodisch( h = 4 . 99949 mm) der Bestimmung der Bifurkationen von periodischen L¨osungen: es ist notwendig, dass die Anfangsbedingungen zum Nachweis periodischer L¨osungen stets innerhalb des Einzugsgebietes der periodischen L¨osungen liegen. Mit weiterru¨ckendem Bifurkationsparameter ver¨andert sich auch das Einzugsgebiet geringfu¨gig. Fu¨r das vorliegende System wurde daher eine Strategie implementiert, die als Anfangsbedingung den Fixpunkt aus der Rechnung mit dem vorigen Parameter verwendet. Bei n -periodischen L¨osungen( n 2) wird derjenige der n Fixpunkte mit dem betragsm¨aßig gr¨oßten Zahlenwert im Geschwindigkeitszustand gew¨ahlt, da dieser maximal weit entfernt vom chaotischen Attraktor liegt, der typischerweise in Bereichen kleinerer Geschwindigkeiten residiert. Durch dieses Vorgehen kann bei genu¨gend kleiner Parameterschrittweite nahezu sichergestellt werden, dass die Anfangsbedingungen mit dem Einzugsgebiet der periodischen L¨osungen mitwandern. Bifurkationsparameter h Im vorigen Abschnitt 6.2.2 wurde eine 1-periodische L¨osung bei einem Flugabstand von h = 5 mm mit dem Fixpunkt V = 1 . 49542278921781 rad 3 . 94968938856343 m/s gefunden. Im Labor l¨asst sich der Flugabstand zwar einstellen, jedoch ist dies nur mit begrenzter Pr¨azision m¨oglich. In der folgenden Bifurkationsanalyse soll daher untersucht werden, welchen Einfluss h auf die Lage und die Art des Fixpunktes hat. Bild 6.13 zeigt zwei Bifurkationsdiagramme mit dem Bifurkationsparameter h . Das linke Diagramm wurde fu¨r die Zustandsvariable V k , das rechte fu¨r die Zustandsvariable k erzeugt. Fu¨r den Parameterwert h = 5 mm finden wir den zuvor gefundenen Fixpunkt[ k , V k ] T in den Diagrammen wieder. Ausgehend von h = 5 mm wurde der Flugabstand in feinen Schritten zun¨achst erh¨oht, dann verringert. W¨achst h an, so steigt die Fixpunktgeschwindigkeit, und die Phasenlage des Fixpunktes n¨ahert sich von unten dem Oszillatorgeschwindigkeitsmaximum bei k = 3 / 2 rad. Die periodische L¨osung existiert mit Periode eins weiterhin, lediglich der Fixpunkt 70

6.2. SIMULATION DES MODELLS 2. ORDNUNG V [m/s] 3.9500 3.9498 3.9496 3.9494 [ rad] 1.5 1.498 1.496 1.494 1.492 4.9995 3.95 5 h[m] 5.0005 x 10 -3 V [m/s] 3.9492 3.9495 3.9490 1.492 1.494 1.496 1.498 1.5 [ rad] 3.949 4.9995 5 h[m] 5.0005 x 10 -3 Bild 6.15: Ortskurve des Bifurkationsdiagramms fu¨r den Laufparameter h verschiebt sich mit steigendem h etwas im Zustandsraum. Sobald der Parameter jedoch den Flugabstand von 5.00066 mm u¨bersteigt, verschwindet die stabile Fixpunktl¨osung pl¨otzlich und weicht der alleine dominierenden chaotischen L¨osung. Fu¨r h> 5 . 00066 mm breitet sich die Punktewolke, die das Chaos beschreibt, bei Geschwindigkeitszust¨anden, die unterhalb des gezeigten Ausschnitts im Bild links liegen und sind im Diagramm daher nicht sichtbar. Verringert man den Parameter ausgehend von der 1-periodischen L¨osung bei h = 5 mm schrittweise, so passiert U¨ berraschendes bei h 4 . 9998 mm: der Ast periodischer L¨osungen verzweigt sich. Der Ast vormals stabiler 1-periodischer L¨osungen wechselt seine Stabilit¨at hin zu einem weiterlaufenden Ast von instabilen Fixpunkten. Da die instabilen L¨osungen durch das beschriebene Vorgehen nicht erhalten werden, wird der Ast 1-periodischer L¨osungen im Diagramm nicht fortgesetzt. Wie aus dem Nichts verzweigen sich aus dem Ast 1-periodischer L¨osungen zwei Zweige. Sie enthu¨llen die Existenz 2-periodischer L¨osungen. Dieses Ph¨anomen nennt man Periodenverdoppelung bzw. Heugabel-Verzweigung(engl. pitchfork bifurkation). Startet das System von einem Anfangszustand aus dem Einzugsbereich der 2-periodischen L¨osung, so konvergiert die Reihe in wenigen Schritten gegen die zwei Fixpunkte. Diesen Vorgang zeigt die Trajektorie in Bild 6.14, linke H¨alfte. Im Zustandsdiagramm ist die Reihenentwicklung als PunktReihe V k von der au¨ubseret w k aei1n5geSzcehicrhitnteet.itDerieieZrtiffwerur"d1"enm. aDrikeieTrtradjieekStotariretbdeedrinegrustnegn, transienten St¨oße(grau gestrichelt) n¨ahert sich etwa ab dem 6. Schritt der station¨aren, stabilen periodischen L¨osung. Die periodische Trajektorie, die 71

KAPITEL 6. VERHALTEN DES STOSSBOHRERS V [m/s] [ rad] 3.9496 3.9495 3.9494 1.5 1.498 3.9493 1.496 3.9492 3.9491 3.949 1.494 1.492 3.9489 0.61200 0.61205 1 0.61210 1.49 0.61200 0.61205 1 Bild 6.16: Bifurkationsdiagramm fu¨r Parameter 1 0.61210 zwischen den zwei gelb hinterlegten Fixpunkten hin- und herspringt, ist tu¨rkis gef¨arbt. Weiteres Verringern von h fu¨hrt zu weiteren Periodenverdoppelungen. Bild 6.14, rechte H¨alfte, zeigt die transiente Trajektorie hin zur stabilen 4-periodischen L¨osung. Mit weiterer schrittweiser Verkleinerung des Flugabstands ergibt sich eine Kaskade immer ¨ofter auftretender Periodenverdoppelungen, die schließlich fu¨r etwa h< 4 . 9994 mm in Chaos mu¨nden. Die Darstellungsweise in den Bifurkationsdiagrammen in Bild 6.13 bew¨ahrt sich, um den Grad der Periodizit¨at der L¨osungen zu erkennen. In Bild 6.15 ist die Ortskurve der periodischen Fixpunkte mit dem Laufparameter h gezeigt, die der Darstellung im Bifurkationsdiagramm entspricht. Farben kennzeichnen den Grad der Periodizit¨at von 1-periodisch(blau) bis 16-periodisch(lila). Ablesen l¨asst sich hier leicht, dass die Bewegungsform mit gr¨oßter Abfluggeschwindigkeit V die 1-periodische Bewegung ist. Bifurkationsparameter 1 In einer weiteren Bifurkationsanalyse wollen wir die Periodizit¨atsexistenz beleuchten, wenn eine der Stoßzahlen variiert wird. In Unterkapitel 6.1 sahen wir, dass die Stoßzahlen in der Praxis schwankende Gr¨oßen sind. Dies motiviert das Interesse der Wahl von 1 fu¨r die Parameterstudie. Mit dem folgenden Parametersatz, fu¨r den wir im vorigen Abschnitt Periode-4-Verhalten 72

6.3. EXPERIMENTELLE UNTERSUCHUNGEN AM STOSSBOHRER sahen, beginnen wir die Variation von 1 = 0 . 6121 aus: Oszillator­ Ampl. Frequenz u 0 f Stoßzahl an Oszillator Bohrer 1 2 Flugabstand h 10 µ m 20 . 19 kHz[0.6120 0.61211] 0.7879 4.9995 mm (6.10) Bild 6.16 zeigt das Bifurkationsdiagramm fu¨r 1 , Bild 4.1 hebt mit einer 200-Schritt-Aufl¨osung des rechten Randes die Periodenverdoppelungskaskade hervor. Je h¨oher 1 , desto geringer ist die Systemd¨ampfung. Entsprechendes l¨asst sich aus den Diagrammen ablesen: Bewegungen mit großen Stoßgeschwindigkeiten ergeben sich fu¨r gr¨oßere Stoßzahlen. Jenseits von 1 = 0 . 61211 herrscht jedoch wieder Chaos vor. Der Bereich, in dem die Stoßzahl variieren kann, um Periodizit¨at zu erm¨oglichen, ist weit schmaler als die Streuung der Werte in der Praxis. Dies erkl¨art, warum periodische L¨osungen in der Praxis so gut wie nie beobachtet werden. 6.3 Experimentelle Untersuchungen am Stoßbohrer Im bisherigen Teil dieses Kapitels wurde das dynamische Verhalten der zentralen Komponente des Ultraschall-Stoßbohrers-- der freien Stoßmasse-untersucht. Der Bohrer, der als Gegenprallelement fu¨r die Stoßmasse auf der einen Seite fungiert, wurde festgehalten und durch eine Hammerbahn ersetzt. Der Antrieb des Bohrers-- der Ultraschall-Transducer, auch Oszillator genannt-- wurde mit einem konstanten Abstand zum Gegenprallelement im Pru¨fstand eingespannt. Die Auswirkungen des in Schwingungsrichtung frei beweglichen Oszillators auf die dynamischen Eigenschaften des Bohrsystems werden in diesem Unterkapitel experimentell und im folgenden simulierend untersucht. 6.3.1 Pru¨fstand fu¨r Messungen am Stoßbohrer Wie im bisherigen Pru¨fstand stellt die Fu¨hrung der Stoßmasse eine Schwierigkeit dar, denn sie soll die Kugel in Flugrichtung nicht beeinflussen, quer dazu jedoch begrenzen. Bild 6.2 auf Seite 55 zeigt den Aufbau eines Stoßbohrers, bei dem als Stoßmasse eine runde Scheibe mit zentralem Loch zur Fu¨hrung gew¨ahlt wurde. Ein Fu¨hrungsstab in L¨angsrichtung verhindert Querbewegungen der Scheibe. Die Reibung an seinem Umfang fu¨hrt jedoch zum Verkanten und Taumeln der Scheibe. Eine verbesserte Version der Fu¨hrung wurde von Dipl.-Ing. Christian Potthast aufgebaut, zu sehen auf den Bildern 6.17 und 6.18 [Potthast et al. 2007a]. Als Stoßmasse bewegt sich eine Stahlkugel mit Durchmesser 8 mm in einer Aluminiumhu¨lse mit geringfu¨gig gr¨oßerem Durchmesser. Eine Aussparung in der Hu¨lse erm¨oglicht das optische Erfassen der Kugelposition. W¨ahrend der Messung sitzt der Bohrer auf einem 73

KAPITEL 6. VERHALTEN DES STOSSBOHRERS Bild 6.17: Gesamter Pru¨fstand mit Linearfu¨hrung fu¨r den Stoßbohrer; Pru¨fstandskonzeption und Bild: C. Potthast Transducer Free Mass Drill Stem Bild 6.18: Stoßbohrer aus Pru¨fstand: Oszillator, Kugel und Bohrer mit Fu¨hrungshu¨lse; Bild: C. Potthast Kraftmesser, auf dem er in vertikaler Richtung nicht fixiert ist. Fu¨r die Messdatenaufzeichnung wurde eine Hochgeschwindigkeitskamera verwendet, die mit einer Frequenz von 9000 Bildern pro Sekunde Photographien speicherte. 74

6.4. SIMULATION DES MODELLS 4. ORDNUNG 6.3.2 Messergebnisse am Stoßbohrer Die Messung am Stoßbohrer mit dem oben beschrieben Pru¨fstand wurde von Dipl.-Ing. Christian Potthast mit dem Parametersatz Ampl. u 0 5 µ m Oszillator­ Frequenz Masse fM 20.19 kHz 264.7 g Kugel­ Masse m 2.1 g Anpress­ Kraft F A 3N GravitationsBeschl. g 9.81 m/s 2 (6.11) durchgefu¨hrt. Einen Ausschnitt von 63 ms der Hochgeschwindigkeitskameraauswertung zeigt Bild 6.19. Sichtbar ist eine Starrk¨orperbewegung des Bohrers (gru¨n), der seinerseits St¨oße mit dem Untergrund(im Pru¨fstand war dies ein Kraftsensor) ausfu¨hrt. Position [mm] 4 3.5 3 2.5 2 1.5 1 0.5 0 0.57 0.58 0.59 0.6 Zeit t[s] 0.61 0.62 Bild 6.19: Gemessene Bewegung von Oszillator, Kugel und Bohrer 6.4 Simulation des Modells 4. Ordnung In den Unterkapiteln 6.1 und 6.2 wurden die m¨oglichen Bewegungsabl¨aufe der Stoßmasse(Kugel) im Inneren des Bohrsystems mithilfe eines auf zwei Zustandsvariablen reduzierten Grundmodells studiert. In Unterkapitel 5.2 wurde durch Verfeinern des vorigen Modells die Fesselung des Oszillators ans Inertialsystem aufgehoben: wie die Stoßmasse kann er sich frei bewegen. Auch wurde der Einfluss der Schwerkraft und einer Anpresskraft mit in das Modell aufgenommen. Diesen Aspekten wird in einer Erweiterung des vormals untersuchten Modells 2. Ordnung Rechnung getragen. Im Gegensatz 75

KAPITEL 6. VERHALTEN DES STOSSBOHRERS zu diesem Modell ist es von Ordnung vier und simuliert ebenfalls diskrete Zeitschritte. Zus¨atzlich zu den beiden Zustandsgr¨oßen des ersten Modells (Stoßphase und Abfluggeschwindigkeit der Kugel) beru¨cksichtigt das erweiterte Modell die Geschwindigkeit und Position des beweglichen Oszillators im Zeitpunkt jedes Stoßes zwischen Stoßmasse und Oszillator. Durch eine einfache Zeitreihensimulation ist es wie zuvor hilfreich und nu¨tzlich, sich zun¨achst einen U¨ berblick u¨ber die zu studierende Dynamik zu machen. Hierzu iterieren wir das Modell 4. Ordnung ausgehend von einem festgelegten Anfangszustand. Die Bewegungen von Stoßmasse und Oszillator halten wir in einer Zeitreihendarstellung fest, wie auch schon auf Seite 66 in Bild 6.11 fu¨r das reduzierte Modell 2. Ordnung. Die Gleichungen(5.43) bis (5.49) iterativ angewendet liefern x k , x k und y k jeweils zu den Stoßzeitpunkten t k . Die Positionen des Oszillatorschwerpunktes x ( t ), der Oszillatorspitze u ( t ) und der Stoßmasse y ( t ) in kontinuierlicher Zeit zwischen den Stoßereignissen werden aus den Anfangswertprobleml¨osungen(5.24), (5.26) und (5.34) berechnet. Fu¨r die Simulationen werden der Parametersatz Oszillator­ Ampl. Frequ. u 0 f 15 µ m 20.19 kHz Stoßzahl an Osz. Bohrer 1 2 0.5 0.6 Massen­ verh¨altn. M/m = µ 264. 7 g 2. 0 g = 132 Beschleunigung aus Anpresskraft Gravitat. F A /M g 3N 0. 2647 kg = 11 m/s 2 9.81 m/s 2 sowie als Anfangsbedingungen fu¨r die Zustandsvariablen Startgeschwindigkeit Startpos. Startzeit Stoßmasse Oszillator t 0 37.7 µ s y 0+ -3.33 m/s x +0 0 m/s x 0 1.42 mm benutzt. Die erste halbe Sekunde der mit den angegebenen Parametern und Startbedingungen simulierten Zeitreihe zeigt Bild 6.20. Im unteren Teil der Zeitreihe ist die Entwicklung der Stoßfrequenzen aufgetragen, d. h. u¨ber jede Stoßzeit ist der Kehrwert der Dauer bis zum n¨achsten Stoß aufgetragen. Die gepunktete Linie( · · · ) bei 1.3 kHz kennzeichnet die mittlere Stoßfrequenz, wobei die(extrem kurzen) Flugdauern von Mehrfachst¨oßen hier und auch im Diagramm nicht auftauchen. Das Histogramm in Bild 6.21 zeigt die H¨aufigkeitsverteilung der Frequenzen: die meisten St¨oße folgen mit einer Frequenz von 500 Hz aufeinander. An der Lage der lokalen Maxima der Stoßfrequenzverteilung unter der Zeitreihendarstellung erkennen wir, dass dichte Stoßfolgen auftreten, wenn die Oszillatorlage ein lokales Minimum einnimmt, bedingt durch minimale Flugwege fu¨r die Stoßmasse. 76

6.4. SIMULATION DES MODELLS 4. ORDNUNG 5 Positionen [mm] Bohrer Flugmasse Oszillator 4 3 2 1 Stoßfrequenz [Hz] 0 10 4 10 3 10 2 0 0.1 0.2 0.3 0.4 0.5 Zeit t[s] Bild 6.20: Zeitreihensimulation(oben): Oszillator, Kugel und feststehender Bohrer; Stoßfrequenzen(unten) mit Markierung( · · · ) der mittleren Stoßfrequenz Detaillierteren Einblick in die Flugbahnen der Kugelmasse zeigen die Bilder 6.22 bis 6.24. Bild 6.22 hebt einen Ausschnitt von 70 ms aus der iterierten Zeitreihe hervor, Bilder 6.23 und 6.24 zeigen Detailausschnitte aus Bild 6.22. Lilafarbene Abschnitte der Fluggeraden markieren Mehrfachst¨oße der Kugel am Oszillator. Mehrfachst¨oße treten auf, wenn die Kugel nach ihrem Oszillatorstoß noch weiter in dessen Richtung fliegt oder wenn sie sich zu langsam in Richtung auf den Bohrer bewegt und von der gleichen Schwingphase des Oszillators eingeholt wird. Selten kommt es vor, dass die Kugel mit geringer Geschwindigkeit den Schwingradius des Oszillators verl¨asst, viele Periodendauern sp¨ater vom Oszillator aber ein weiteres Mal getroffen wird. Auch diese M¨oglichkeit muss im Modell beru¨cksichtigt werden. Vereinzelt kommt es zu Mehrfachst¨oßen am Bohrer. Sie mathematisch zu beschreiben fu¨hrt auf das Paradoxon von Zeno, wenn die Kugel unendlich oft auf den Bohrer trifft, bevor sie wieder in Kontakt mit dem Oszillator kommt. 13 Trifft der Oszillator erst nach dieser endlichen Zeit unendlich vieler 13 Eines von Zeno von Eleas vier Paradoxen von Zeit und Bewegung war, dass Zeit unendlich unterteilbar ist. Die Dichotomie besagt: Um das Intervall[0 , 1] zu durchschreiten, muss zun¨achst einmal die erste H¨alfte dieser Distanz, also das Intervall[0 , 1 2 ] traversiert werden. Hierzu ist es jedoch wiederum notwendig, dessen erste H¨alfte, sprich das Intervall[0 , 1 4 ] zu durchmessen, usw. Folglich sind zun¨achst unendlich viele ReiseEinheiten notwendig, um sich u¨berhaupt vom Start zu entfernen. Dies ist in endlicher Zeit nicht m¨oglich-- daher kann der Reisende nie am Ziel ankommen. Aristoteles erkl¨arte, Achilles k¨onne die Schildkr¨ote niemals einholen, die 10m vor ihm starte und 10 77

KAPITEL 6. VERHALTEN DES STOSSBOHRERS 0 1 2 3 4 5 6 7 8 9 10 Stossfrequenzen[kHz] Bild 6.21: Histogramm­Verteilung der Stoßfrequenzen: die typische Stoßfrequenz liegt bei 500 Hz. St¨oße auf die dann ruhende Kugel, so kann er sie gem¨aß dem in Unterkapitel 5.1 aufgestellten Starrk¨orpermodell nie mehr in Bewegung versetzen. Laborversuche zeigen jedoch, dass dies m¨oglich ist. Fu¨r diese Situation versagt das Modell folglich-- es muss auf ein elastostatisches Modell gewechselt werden. Im Rahmen dieser Arbeit wurde in solchen F¨allen-- anstelle einer Anpassung des Modells-- von weiteren Iterationen abgesehen. 6.5 Vergleich von Modellen und Messungen Im Sichtvergleich der Zeitsimulation in Bild 6.22 mit dem Ausschnitt aus der Messaufzeichnung in Bild 6.19 l¨asst sich eine hohe U¨ bereinstimmung mal langsamer sei. Denn nachdem Achilles die ersten 10 m gelaufen war, hatte sich die Schildkr¨ote um 1m weiterbewegt. In der Zeit, die Achilles aber fu¨r diesen Meter gebraucht hatte, war die Schildkr¨ote um 1 10 m voran gekrochen, usw.. Die unendlich vielen Etappen, die Achilles zuru¨cklegen muss, addieren sich jedoch auf nicht mehr als 10 1 + 10 0 + 10 1 + 10 2 + ... = 1 1 10 k = 10+ . 1 1 1 10 Kommen wir auf die frei fliegende Kugel zuru¨ck, lierte: Angenommen, die fu¨r die Kugel Whitrow pralle mit die der gGleeiscchhewDinidcihgokteoimt v ie 0 ("inZwveeirtteikilaulnegr"R) icfohrtmunugvom feststehenden Bohrer ab und beschleunigt, ohne den Oszillator zu treffen, mit der Erdbeschleunigung g wieder zum Bohrer zuru¨ck, so wird diese Flugdauer 2 v 0 g dauern. Fu¨r den n¨achsten Stoß betr¨agt die Abprallgeschwindigkeit nur noch v 0 · 2 , usw. Die Kugel wird unendlich viele sich verkleinernde und dichter folgende St¨oße mit der Prallplatte ausu¨ben, also scheinbar nie zur Ruhe kommen. Allerdings ben¨otigen diese unendlich vielen sich verku¨rzenden St¨oße nur eine Zeitspanne von 2 v 0 g + 2 v 0 2 g + 2 v 0 22 g + ... = 2 v 0 g 1 1 2 , wonach die Kugel bewegungslos auf dem Bohrer liegen muss. 78

Positionen [mm] 6.5. VERGLEICH VON MODELLEN UND MESSUNGEN 4 3.5 3 2.5 2 1.5 1 0.5 0 0.38 0.39 0.4 0.41 0.42 0.43 0.44 Zeit t[s] Bild 6.22: Ausschnitt aus einer Zeitreihensimulation: Oszillator, freie Stoßmasse und feststehender Bohrer zwischen Messung und Modell feststellen. Allerdings muss der Vergleich mit Vorsicht gezogen werden, da die Bedingungen des Gegenpralls am Bohrstab unterschiedlich sind. Das Modell beru¨cksichtigt weder die Wellenausbreitungsvorg¨ange im Bohrstab, noch dessen Starrk¨orperbewegungen. Die Vor- und Nachgeschichte der Bewegungsabl¨aufe von Simulation und Messung gleichen sich nicht in gleichem Maße wie die gezeigten Ausschnitte, sondern verlaufen ihren deterministischen Regeln folgend chaotisch. Die mittlere H¨ohe, in die sich der Oszillator hebt, liegt zwischen Messung und Simulation dicht beieinander. Einen leichten Unterschied gibt es in der Stoßfrequenz. In der Simulation liegt sie-- wie auch durchschnittlich im gezeigten Ausschnitt-- u¨ber der Stoßfrequenz in der Messung. Dies wird auf den Bohrer mit seiner h¨oheren Aufnahmef¨ahigkeit kinetischer Energie und D¨ampfung am Untergrund zuru¨ckzufu¨hren sein. Die Stoßzahlen des Modells wurden fu¨r die Simulation so angepasst, dass die vom Oszillator erreichten mittleren H¨ohen sich in Messung und Simulation gleichen. Messungen der Kontaktkr¨afte des Bohrers in[Potthast et al. 2007a] zeigen, dass die Kraftspitzen und vermutlich die gr¨oßten Bohrraten dann erfolgen, wenn die Oszillatorposition ein lokales Minimum einnimmt. Dieses Verhalten ist in der Simulation insofern wiederzufinden, als dass die Stoßfrequenz in solchen Phasen Maximalwerte erreicht. Die Mechanismen, die diesem Ph¨anomen jeweils zugrunde liegen, k¨onnen aber in Realit¨at und Simulation unterschiedliche sein. Im Laborversuch ist es m¨oglich, wenn auch nur fu¨r sehr kurze Zeitr¨aume, dass Oszillator, Kugel und Bohrer 79

KAPITEL 6. VERHALTEN DES STOSSBOHRERS 0.3 0.25 Positionen [mm] 0.2 0.15 0.1 0.05 0 0.379 0.3795 0.38 Zeit t[s] 0.3805 0.381 Bild 6.23: Detail aus Bild 6.22 bei t = 0 . 38 s: in dieser Vergr¨oßerung ist die harmonische Schwingung des Oszillators u ( t ) sichtbar. Die Flugbahn der freien Masse ist lilafarben gezeichnet, w¨ahrend die Stoßmasse einen Mehrfachstoß am Oszillator erf¨ahrt. aufeinander liegen, also alle gleichzeitig die Position 0 m einnehmen. Somit u¨bertr¨agt der schwerere Oszillator seinen Stoßimpuls gelegentlich direkt an den Bohrstab, was die Stoßkraftmaxima in[Potthast et al. 2007a] erkl¨art. Im Modell hingegen entspricht der Zustand eines Dreierkontaktes dem Fglleuigchabbestdaenudten"dNuilsl"t,, der einem unendlich was wiederum eine kleinen zeitlichen Stoßabstand unendlich hohe Stoßfrequenz symbolisiert. Diese U¨ berlegungen wu¨rden in der Theorie jedoch nur mit einem nicht-schwingenden Oszillator harmonisieren. Die Schwingungen des Oszillators begrenzen die Dauer solch eines Dreierkontaktes und limitieren auch die Stoßwiederholrate maximal auf die Anregefrequenz. Zusammenfassend l¨asst sich sagen, dass zwischen beiden simulierten Modellen und dem realen Verhalten in Laborpru¨fst¨anden eine große qualitative U¨ bereinstimmung festgestellt wurde. Quantitativ bewegen sich die Zust¨ande in Theorie und Praxis in den gleichen Gr¨oßenordnungen. In beiden Modellen sind die Stoßzahlen die einzigen Parameter, die im Labor nicht eingestellt oder abgelesen werden k¨onnen. Eine messtechnische Bestimmung der Stoßzahlen von u¨ber 500 St¨oßen hat eine Streuung der Werte gezeigt. Obwohl im Modell zeitlich konstante Stoßzahlen eingesetzt werden, war es m¨oglich, sie so einzustellen, dass sich die Oszillatorbewegungen in Simulationen und Experimenten in ann¨ahernd gleicher Flugh¨ohe bewegen. Der Grund fu¨r unterschiedliche Stoßfrequenzen der Stoßmasse im Modell 4. Ordnung und im Ver80

6.5. VERGLEICH VON MODELLEN UND MESSUNGEN 0.35 0.3 Positionen [mm] 0.25 0.2 0.15 0.1 0.05 0 0.43 0.4305 0.431 Zeit t[s] 0.4315 0.432 Bild 6.24: Detail aus Bild 6.22 bei t = 0 . 43 ms: durch die hier geringe H¨ohe des Oszillators u¨ber dem Bohrer(Flugabstand fu¨r Kugelmasse nur 0.3 mm) erh¨oht sich die Stoßfrequenz der Kugelmasse auf 3 kHz, wodurch die Fallbewegung des Oszillators gebremst wird. such am Stoßbohrsystem resultieren aus den verschiedenen Bohrerbedingungen: im Labor war dieser beweglich gefu¨hrt und konnte somit durch sprunghaftes A¨ ndern seiner Starrk¨orpergeschwindigkeit im Stoßzeitpunkt mit der Stoßmasse von dieser mehr kinetische Energie abfu¨hren als im Modell, in dem der Bohrer als starre, teilelastische Prallplatte modelliert wurde. Fu¨r die exakte Abbildung des beweglichen Bohrstabes im Modell k¨onnte eine weitere Verfeinerungsstufe gew¨ahlt werden, die aber nicht mehr im Umfang dieser Arbeit lag. 81

.

Kapitel 7 Numerische Analyse mit Hilfe der mengenorientierten Methoden Kapitel 7 behandelt den Schwerpunkt dieser Arbeit, n¨amlich die Anwendung der mengenorientierten numerischen Methoden. Ihr vielversprechender Einsatz wurde bereits durch mehrere Anwendungen der Methoden auf Systeme mit komplizierter Dynamik nachgewiesen. Die zugrundeliegenden Modelle konnten in solchen F¨allen durch explizite Bewegungsgleichungen in wenigen Zeilen formuliert werden und waren glatter und kontinuierlicher Art, also ohne sprunghafte Wechsel der Bewegungsgleichungen. Ihre Anwendbarkeit auf diskrete, also nicht-kontinuierliche Systeme bzw. solche mit Sprungph¨anomenen wurde bislang noch nicht demonstriert. Weiterhin ist die Verwendung bei Modellen mit transzendenten Bewegungsgleichungen nicht bekannt. Im Rahmen der vorliegenden Arbeit wird der Einsatz der mengenorientierten Methoden erstmals auf ein nichtglattes, zeitdiskretes System demonstriert und beschrieben. Insbesondere wird ein System aus der Praxis behandelt. Dies stellt bereits an die mathematische Modellbildung hohe Anforderungen, auf die in Kapitel 5 eingegangen wurde. Kapitel 7 ist in zwei Teile unterteilt. Im ersten Teil wird das Analysevorgehen anhand des Systems 2. Ordnung demonstriert, das in Unterkapitel 5.1 aufgestellt wurde. Dieses bewusst einfach gehaltene und auf zwei Dimensionen beschr¨ankte Modell erm¨oglicht ein eingehendes Verst¨andnis der Funktionsweise der Methodik. Die Ergebnisse der Analyse lassen sich bequem im 2-dimensionalen Zustandsraum auf der Papierebene darstellen, und die Eigenschaften und Vorteile der Methoden werden daran besprochen. Der zweite Teil dieses Kapitels demonstriert, wie sich der Einsatz der Methoden auf h¨oherdimensionale Systeme, wie sie der Realit¨at n¨aher kommen, eignet. Hierzu verwenden wir das Modell 4. Ordnung, das in Unterkapitel 5.2 hergeleitet wurde. Seine h¨ohere Dimensionszahl erschwert die Interpretation der Ergebnisse, die nur unter Verwendung geeigneter Projektionen auf die 2-dimensionale Papierebene visualisiert werden k¨onnen. 83

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN 7.1 Analyse des Stoß-Modells 2. Ordnung Die Algorithmen, welche die mengenorientierten Methoden ausmachen, sind alle im Softwarepaket Gaio implementiert[Dellnitz, Froyland, Junge 2001]. Gaio steht fu¨r Global Analysis of Invariant Objects. Die Software wird in der und AdyrbneaimtsigsrcuhpepSeys"tAemngee"waamndItnesMtitautthfeu¨mraMtikath--emNautimkedreisrchUeniMveartshite¨amt aPtaikderborn entwickelt und stetig erweitert. Die Softwarepakete sind in der Programmiersprache C geschrieben, daher mu¨ssen die zu untersuchenden Modellgleichungen in dieser Sprache kodiert werden. Hierfu¨r steht eine Schnittstelle zur Verfu¨gung, in der die rechte Seite des Differentialgleichungssystems erster s ide" Ordnung definiert wird. Das geschieht in der C-Funktion " r ight h and void rhs(double*x, double*u, double*y). Im Vektor x werden die Zustandsgr¨oßen dem Modell u¨bergeben; der Vektor y wird w¨ahrend der Modellauswertung berechnet und enth¨alt danach den berechneten Zustand einen Zeitschritt d t sp¨ater(bei kontinuierlichen Modellen) bzw. bei diskreten Modellen einen Iterationsschritt sp¨ater. Im Kopf der Modell-Datei lassen sich die Modellparameter definieren und Einstellungen vornehmen: hier wird angegeben, um welche Art der Modellbeschreibung c(dhuisnkgresstyesAtebmbi"ld o u d n e "g)sveosrsscichhrifhtan" m d a e p lt". oder kontinuierliches Differentialglei7.1.1 Der zylindrische Zustandsraum An dieser Stelle ist es n¨otig, den fu¨r die vorliegende Art von Schwingstoßsystemen angebrachten Phasen- oder Zustandsraum einzufu¨hren. Der Zustandsraum fu¨r den betrachteten Fall hat zwei Dimensionen. Dies sind die Stoßzeit t k und die Abfluggeschwindigkeit V k . Man k¨onnte die Zeit herk¨ommlicherweise auf die Abszisse auftragen und die Geschwindigkeit auf die Ordinate. Ein Systemzustand entspr¨ache einem Koordinatenpunkt in diesem Koordinatensystem und wu¨rde durch einen Abbildungsschritt auf einen weiter rechts liegenden Zustand abgebildet werden. Durch diese Darstellungsweise von Zustandspunkten k¨onnte periodisches oder anders strukturelles Verhalten nicht identifiziert werden. Die(Uhr-)Zeit der aufeinanderfolgenden Stoßzust¨ande spielt in der Analyse eine untergeordnete Rolle, n¨amlich nur zur Bestimmung der Stoßfrequenzen(Stoßwiederholraten); vielmehr ist es entscheidend, in welcher Phasenlage der harmonischen Oszillatorschwingung ein Stoß am Oszillator auftritt. Daher ist es naheliegend, die Stoßzeit t k [s] mit k = t k 2 T mit T = 1 /f zun¨achst auf das Winkelmaß k [rad] zu konvertieren und dieses dann mit k = k modulo 2 84

7.1. ANALYSE DES STOSS-MODELLS 2. ORDNUNG auf eine einzige Periodendauer zuru¨ckzumodulieren 1 . Anstelle des kartesischen Phasenraumes R × R benutzen wir folglich den zylindrischen Phasenraum S 1 × R , der zur besseren Anschauung in Bild 7.1 skizziert ist. Die Bewegungsgleichung des Systems 2. Ordnung ist eine Abbildung des Zylinders auf sich selbst. V k / 2 3 / 2 2 k Bild 7.1: Zylinderf¨ormiger Zustandsraum fu¨r Phasen- und Geschwindigkeitszustand, nach unten durch gru¨ne Sinuslinie begrenzt. rechts: abgerollter Zustandsraum. die DUineetnodploiclhokgeisicthdee"rZzyeliitnlidcrhisfioerrutsncgh"rediteesnPdhenasEennrtawuimckelsunergmin¨ongleirchhatlebs nun, einer Periodendauer darzustellen. Alle Phasendiagramme zum Modell 2. Ordnung in diesem Unterkapitel nutzen diese Darstellung, bei der sie, aufgeschnitten an der 0 Ordinate 2  Ordinate 2 , flach abgerollt wiedergegeben sind. Der stets in Rechteckform gezeichnete Zustandsraum ist u¨ber seine rechte Seite hinaus in seine linke Seite hinein(und umgekehrt) fortgesetzt zu lesen. Somit liegt ein Zustandspunkt an der rechten Seite einem Punkt an der linken Seite sehr nahe. Der abgerollte, flach dargestellte Zustandsraum zeigt sich in der Phasendimension zu beiden Seiten begrenzt. Beim Stoßbohrer ist er außerdem in der Geschwindigkeitsdimension begrenzt: die Ab1 Modulofunktion: fu¨r a, b N gilt: a mod b = Rest von a b ; fu¨r den vorliegenden Fall ist die allgemeinere Definition n¨otig: fu¨r a R + , b R + \ 0 gilt: a mod b = a nb mit n = a b , wobei n auf die n¨achst-tiefere ganze Zahl abgerundet werde. In Worten: a modulo b ist der verbleibende Rest, nachdem man b sooft wie m¨oglich( n mal) von a subtrahiert hat, wobei man sich in den positiven reellen Zahlen bewege. 2 Ordinate: senkrechte Koordinatenachse, oft "y-Achse" genannt. 85

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN fluggeschwindigkeiten V k k¨onnen niemals kleiner als die Oszillatorgeschwindigkeit im Stoßmoment sein. Andernfalls mu¨sste die Stoßmasse nach dem Stoß in den Oszillator eindringen, was physikalisch nicht m¨oglich ist. Die Nicht-Eindringbedingung lautet dementsprechend: V k  - u 0  sin t k . (7.1) Diese untere Schranke ist in Bild 7.1 und in den Zustandsdiagrammen als gru¨ne Kurve eingezeichnet. Sie begrenzt die physikalisch m¨oglichen Zust¨ande nach unten hin. Der gesamte Zustandsraum ist also nur noch nach oben in Richtung+ V k unbeschrankt und ¨ahnelt somit einer Litfaßs¨aule 3 mit schr¨ag abgeschnittenem Boden und unbegrenzter H¨ohe. Bild 7.1 zeigt diesen sogeformten Zustandsraum, der von S 1 × R = { ( , V ) | 0  < 2 ,- u 0  sin  < V<} beschrieben wird. Die Darstellungsweise auf dem Zylinder erm¨oglicht das Erkennen von Fixpunkten, mehr-periodischen L¨osungen und seltsamen Attraktoren. 7.1.2 Bestimmung eines relativen globalen Attraktors Unterkapitel 2.2 beschreibt, wie sich der relative globale Attraktor eines dynamischen Systems approximieren l¨asst. Relativ zu einer Untermenge Q des globalen Zustandsraumes n¨ahert er die attraktive Menge an. Diese umfasst jene Systemzust¨ande, die vom System aufgrund seiner periodischen oder chaotischen Dynamik eingenommen werden k¨onnen. Folglich bildet ein Attraktor eine Teilmenge des Zustandsraumes. Startet das System mit Anfangsbedingungen, die weit außerhalb des begrenzten Attraktorgebietes liegen, so entwickelt sich das System innerhalb weniger Iterationsschritte hin zum Attraktor, auf dem alle zuku¨nftigen Zust¨ande verbleiben. Zielsetzung dieser Arbeit ist das Aufzeigen der Analysem¨oglichkeiten durch Bestimmung relativer globaler Attraktoren A Q mit mengenorientierten Methoden. Fu¨r zwei verschiedene Parameters¨atze wird in diesem und den folgenden Abschnitten das Vorgehen anhand des Systems 2. Ordnung erl¨autert, das in Bild 5.1, Seite 35, skizziert und in Gleichung(5.9), Seite 38, modelliert ist. Zun¨achst werden die fu¨nf Parameter wie folgt gesetzt: u 0 1 µ m 2  · 22 kHz 1 0.9 2 0.9 h 20 µ m . (7.2) Der Boxunterteilungs-Algorithmus(siehe Unterkapitel 2.2) zur Ann¨aherung des relativen globalen Attraktors wurde 14 mal angewandt. Das heißt, in jeder Dimension wurde der zu untersuchende Anfangsbereich (Startbox Q :) 7 mal unterteilt. Dies fu¨hrt auf eine Aufl¨osung von 2 14 Boxen bezogen auf die Startbox. Die 5821 Boxen, die den Attraktor u¨berdecken, 3 Litfaßs¨aule: Eine Anschlags¨aule, benannt nach dem Berliner Buchdrucker E. Litfaß. 86

Abfluggeschwindigkeit V [m/s] 7.1. ANALYSE DES STOSS-MODELLS 2. ORDNUNG zeigt Bild 7.2. Innerhalb der untersuchten Startbox nimmt der Attraktor etwa 36% der Zust¨ande ein. H¨atte man die Startbox zu Beginn in 2 14 Boxen unterteilt, w¨aren 64% der Zellabbildungen unn¨otig vorgenommen worden.(Allerdings waren zur Approximation des Attraktors 14 Unterteilungsschritte n¨otig; in jedem Schritt wurden die vorhandenen Unterboxen untersucht. Dies erh¨oht die Anzahl abzubildender Boxen.) Zur Abbildung der Boxen wurde ein Gitter mit 30 × 30 Testpunkten in jede Box gelegt und deren Bildpunkte berechnet. 1 0.8 0.6 0.4 0.2 0 0 0.5 1 1.5 2 Stoßphase [ rad] Bild 7.2: Ein relativer globaler Attraktor A Q , der nach 14 Unterteilungsschritten durch U¨ berdeckung mit 5821 Boxen angen¨ahert wurde. 7.1.3 Berechnung der Aufenthaltswahrscheinlichkeiten Im vorigen Abschnitt wurde der globale Attraktor zum Parametersatz(7.2) durch U¨ berdeckung mit Boxen bestimmt. Wir lesen an ihm ab, welche m¨oglichen Zust¨ande das System in seinen s¨amtlichen, m¨oglichen Bewegungsformen einnehmen kann, nachdem transiente Einschwingvorg¨ange abgeklungen sind. Im vorliegenden Fall wird ersichtlich, dass Abfluggeschwindigkeiten gr¨oßer 1.1 m/s nie auftreten. Nachdem wir wissen, welche Zust¨ande vom System w¨ahrend seiner zeitlichen Evolution prinzipiell eingenommen werden, stellt sich die Frage, ob es Priorit¨aten bezu¨glich bestimmter Zust¨ande, die dem Attraktor angeh¨oren, gibt. Es ist m¨oglich, die Information, welche Zust¨ande wie wahrscheinlich eingenommen werden, zu erhalten. Anders ausgedru¨ckt l¨asst sich berechnen, wie oft sich im Durchschnitt jeder Systemzustand einstellen wird. Die Genauigkeit fu¨r diese Angabe ist die Feinheit der Boxen. Fu¨r jede Box wer87

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN den wir den Wert ihrer berechnete Information "bAezuu¨fgenlitchhadltesrwazehirtslcichheeinnliEchnktweiitc"kleurnhgaltdeenr. Da diese Dynamik unver¨anderlich ist, bildet Maß", siehe Unterkapitel die Aufenthaltswahrscheinlichkeit 2.3. ein "invariantes Wie in Unterkapitel 2.3 beschrieben, wird zur Berechnung der Aufenthaltswahrscheinlichkeiten die sogenannte U¨bergangswahrscheinlichkeitsmatrix P ben¨otigt. Im entsprechenden Abschnitt auf Seite 13 wurde definiert, dass die U¨ bergangswahrscheinlichkeit die Wahrscheinlichkeit angibt, mit der das System innerhalb eines Iterationsschrittes(also von einem Stoß zum n¨achsten) von einem bestimmten Zustandsbereich in einen anderen(oder den selben) Zustandsbereich wechselt. Jede der Boxen umfasst einen Zustandsbereich. Hier wurden 30 2 Testpunkte auf einem rechtwinkligen Gitter in jeder Box gew¨ahlt, um ihr Bild zu approximieren(vgl. schematische Darstellung 2.2, Seite 14). Da jede der 5821 Boxen des relativen globalen Attraktors auf jede andere Box abbilden kann, besitzt die U¨ bergangswahrscheinlichkeitsmatrix die Gr¨oße 5821 × 5821. Die meisten Boxen bilden nur auf wenige Boxen ab, weshalb viele Eintr¨age der Matrix verschwinden. Die Matrix ist daher du¨nnbesetzt(engl.: sparse), wofu¨r es spezielle Eigenwertl¨oser innerhalb der Software Matlab 4 gibt. Diese werden genutzt, um die in Gleichung(2.9) formulierte Eigenwertaufgabe zu l¨osen. Das Ergebnis davon zeigt, dass die hiesige U¨ bergangswahrscheinlichkeitsmatrix zwei Eigenvektoren y 1 und y 2 mit jeweils dem Eigenwert Eins besitzt. Folglich existieren innerhalb des relativen globalen Attraktors A Q fu¨r den einen festen Parametersatz zwei Aufenthaltswahrscheinlichkeitsverteilungen, die jede fu¨r sich fu¨r eine eigene Art der Dynamik stehen. Aus der Theorie der Eigenwerte ist aber bekannt, dass jede Linearkombination der zwei Eigenwerte, y =  y 1 +  y 2 ,, R , wiederum einen gu¨ltigen Eigenvektor mit Eigenwert Eins liefert. Genau zwei voneinander linear unabh¨angige 5 Vektoren dieser unendlich vielen, m¨oglichen neuen Eigenvektoren beschreiben die zwei Bewegungsm¨oglichkeiten. Zwei Bedingungen mu¨ssen fu¨r diese zwei neuen Eigenvektoren gelten, damit ihre Eintr¨age im Kontext zweier konkurrierender Attraktoren physikalisch sinnvolle Aufenthaltswahrscheinlichkeiten enthalten. Mit Hilfe dieser beiden Bedingungen k¨onnen die Linearfaktoren und bestimmt werden. Zur Erinnerung sei erw¨ahnt, dass die Eintr¨age y ¯ 1 k bzw. y ¯ 2 k der Eigenvektoren je eine von zwei m¨oglichen Aufenthaltswahrscheinlichkeiten der Box B k des relativen globalen Attraktors A Q enthalten. Die zwei Bedingungen lauten: 1. Ausschlussbedingung: Liegt der Systemzustand aufgrund einer Bewegung gem¨aß des einen Attraktors mit einer endlichen Wahrscheinlichkeit in einer bestimmten 4 Matlab: Rechenumgebung mit großer Bibliothek numerischer Algorithmen der Firma The MathWorks, Inc. , Boston, zur Programmierung sowie numerischen und symbolischen Rechnung mit komplexwertigen Matrizen. 5 Lineare Unabh¨angigkeit zweier Vektoren bedeutet, dass die Vektorpfeile nicht durch Verl¨angern oder Verku¨rzen ineinander u¨berfu¨hrbar sind. 88

7.1. ANALYSE DES STOSS-MODELLS 2. ORDNUNG Box, so kann der Systemzustand aufgrund einer Bewegung gem¨aß des anderen Attraktors nie, also mit verschwindender Wahrscheinlichkeit, in dieser Box liegen. Einfacher ausgedru¨ckt: wo der eine Attraktor eine von Null verschiedene und positive Zahl enth¨alt, muss der andere Attraktor eine Null enthalten. Gleiches gilt umgekehrt und muss fu¨r jede der N Boxen des relativen globalen Attraktors erfu¨llt sein. Erfu¨llt sein mu¨ssen: y ¯ 1 k = 0  y ¯ 2 k = 0 und y ¯ 2 k = 0  y ¯ 1 k = 0 , k = 1 , 2 ,..., N. (7.3) 2. Summenbedingung: Durch die attraktive Eigenschaft jedes Attraktors befindet sich der Systemzustand(im station¨aren Verhalten) 100%ig auf dem Attraktor. Da die Wahrscheinlichkeit, dass der Zustand irgendwo auf dem Attraktor liegt, gleich der Summe der Einzelwahrscheinlichkeiten aller den Attraktor u¨berdeckenden Boxen ist, gilt immer y ¯ ik = 1 , i = 1 , 2 . k (7.4) Weitere Details zur Wahl geeigneter Linearfaktoren sind [Dellnitz & Junge 1999] zu entnehmen. 7.1.4 Analyseergebnisse aus Aufenthaltswahrscheinlichkeiten Im vorigen Unterkapitel haben wir zwei m¨ogliche Verteilungen der Aufenthaltswahrscheinlichkeiten auf das in Boxen unterteilte Attraktorengebiet berechnet. Die Resultate haben wir in den Vektoren y ¯ 1 und y ¯ 2 abgelegt. Zur Interpretation der Daten werden wir die Wahrscheinlichkeitswerte einer Farbskala zuordnen und die Boxen entsprechend den in y ¯ 1 bzw. in y ¯ 2 abgelegten Wahrscheinlichkeitswerten einf¨arben. Fu¨r den in(7.2) gegebenen Parametersatz sind die Wahrscheinlichkeitsverteilungen in den Bildern 7.3 bzw. 7.4 gezeigt. Bild 7.3 enth¨alt fu¨r beinahe alle Boxen die Nullwahrscheinlichkeit (dunkelblau gef¨arbte Boxen). Die auftretenden Systemzust¨ande beschr¨anken sich gem¨aß der Grafik nur auf einen kleinen Bereich des globalen Attraktors um den Zustand 0.96 m/s und 1 . 24 rad. Nach dem Abklingen instation¨arer Vorg¨ange konzentriert sich der Systemzustand nur noch um einen einzigen Punkt, dem er sich spiralf¨ormig asymptotisch n¨ahert. Hierbei handelt es sich um einen periodischen Punkt: Stoß fu¨r Stoß l¨auft unter immer wiederkehrend gleichen Bedingungen ab, also bei gleicher Phase und mit gleicher Geschwindigkeit. Von den Anfangsbedingungen oder St¨orungen, die auf das System einwirken, h¨angt es ab, ob die periodische L¨osung eingenommen wird oder die, u¨ber welche Bild 7.4 Auskunft gibt. Dort sind all jene Boxen eingef¨arbt, deren Zust¨ande unter der periodischen Bewegung nicht eingenommen werden. Demnach befindet sich der Systemzustand mit unterschiedlicher 89

Abfluggeschwindigkeit V [m/s] KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN Wahrscheinlichkeit auf dem u¨berwiegenden Teil des globalen Attraktors. Im Vergleich der beiden Bilder sehen wir die Erfu¨llung der Ausschlussbedingung 1, Gleichung(7.3): farbige, nicht dunkelblaue Boxen in Bild 7.3 besitzen in Bild 7.4 Nullwahrscheinlichkeit(dunkelblau und umrandet) und umgekehrt. 18 1.2 16 1 14 0.8 12 0.6 10 8 0.4 6 0.2 4 02 0 0.5 1 1.5 2 Stoßphase [ rad] Bild 7.3: Aufenthaltswahrscheinlichkeiten[%] des periodischen Attraktors y ¯ 1 innerhalb des relativen globalen Attraktors A Q . Rot im Hintergrund ist das Einzugsgebiet der periodischen L¨osung dargestellt. Die Ermittlung eines globalen Attraktors und die Berechnung des Aufenthaltswahrscheinlichkeitsfeldes vermittelt neuartige Einblicke in das dynamische Verhalten eines Prozesses. Die komplizierte Dynamik wird in r¨aumlicher und zeitlicher Dimension ausgedehnt beru¨cksichtigt und mit tiefgehender U¨ bersicht visualisiert. Die Begrenzung des globalen Attraktors gibt schnell einen Eindruck der vom System nie eingenommenen Zust¨ande. In Bild 7.4 erkennen wir, dass in einer chaotischen Bewegung nie Abfluggeschwindigkeiten oberhalb von 0.67 m/s zu erwarten sind, und dies ist unabh¨angig von den Anfangsbedingungen. Unter Beru¨cksichtigung der periorudnisgc"h)enerFkeonrtnstetmzuanngwdeeisteZrhuisnt,anddassrsaduemr sAitntrPahkatosernedinimzeunssaimonm(e"nZhy¨alinngdernisdieesGebiet beschreibt. Die Farbkodierung gibt wieder, wie unterschiedlich oft die in Boxen gefassten Bereiche auftreten. Aus diesen Erkenntnissen folgt der Schluss, dass der zweite Attraktor eine deterministisch chaotische Dynamik repr¨asentiert. Innerhalb des farbigen Gebietes springt der Zustand scheinbar regellos. Doch auch im Chaos gibt es Struktur: Zust¨ande hoher Geschwindigkeiten und Phasen um 3/2 treten bevorzugt auf. Auf die Anwendung des Stoßbohrers u¨bertragen stehen hohe Stoßgeschwindigkeiten fu¨r 90

7.1. ANALYSE DES STOSS-MODELLS 2. ORDNUNG große Stoßenergien; die Phase 3/2 entspricht dem Maximum der Oszillatorschwingungsgeschwindigkeit, welches die hohen Abfluggeschwindigkeiten begru¨ndet. Durch Langzeitsimulationen der chaotischen Dynamik w¨are diese Struktur nicht erkennbar gewesen; das h¨aufige Auftreten von Zust¨anden, die den Bohrfortschritt begu¨nstigen, erh¨alt erst in dieser Analyse statistische Quantit¨at. Eine Bemerkung sei zu den Zahlenwerten der Wahrscheinlichkeitswerte fu¨r die einzelnen Boxen gegeben. Die Verteilung der Wahrscheinlichkeiten infolge der chaotischen Dynamik in Bild 7.4 gibt auf der Farbskala kleine Zahlen wieder. Sie liegen zwischen 0 und 0.14%. Diese Zahlen stehen jede fu¨r sich fu¨r eine relative Wahrscheinlichkeit. Sie bezieht sich auf die Gr¨oße (Bgoexn,audeers:toda s s el " te V n o e lu r mf¨aenll"t)ddeerr einzelnen Box, fu¨r Zustand in genau die den sie gilt. Je kleiner Zustandsbereich, die der von der Box umfasst wird. Dementsprechend kleiner wird der angezeigte Wahrscheinlichkeitswert fu¨r diese Box. Durch Erfu¨llung der Summenbedingung(7.4) ergibt die Integration der Wahrscheinlichkeitsverteilung u¨ber den gesamten Zustandsraum stets Eins. Abfluggeschwindigkeit V [m/s] 0.14 1.2 0.12 1 0.1 0.8 0.08 0.6 0.06 0.4 0.04 0.2 0 0.02 0 0 0.5 1 1.5 2 Stoßphase [ rad] Bild 7.4: Aufenthaltswahrscheinlichkeiten[%] des chaotischen Attraktors y ¯ 2 innerhalb des relativen globalen Attraktors A Q . Der Bereich weißen Hintergrunds markiert das Einzugsgebiet dieses chaotischen Attraktors. Die bisherigen Ergebnisse wurden mit dem Parametersatz(7.2) erzielt. Im Folgenden wird die Bestimmung des globalen Attraktors und der Aufenthaltswahrscheinlichkeiten an einem zweiten Parametersatz demonstriert. Die Anregeamplitude u 0 sowie der Flugabstand h werden auf gr¨oßere Werte gesetzt, die in Unterkapitel 6.2 fu¨r Laborversuche und Simulationen einge91

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN stellt wurden: Oszillator­ Ampl. Kreisfrequenz u 0 10 µ m 2  · 20 . 19 kHz Stoßzahl an Osz. Bohrer 1 2 0.6121 0.7879 Flugabstand h 5 mm . (7.5) Abfluggeschwindigkeit V [m/s] 0.9 4 0.8 3 0.7 0.6 2 0.5 1 0.4 0.3 0 0.2 0.1 -1 0 0.5 1 1.5 2 Stoßphase [ rad] Bild 7.5: Boxu¨berdeckung des relativen globalen Attraktors fu¨r den Parametersatz(7.5). Die Farbkodierung zeigt die Aufenthaltswahrscheinlichkeiten [%] des chaotischen Attraktors. Daru¨ber liegt in roter Farbe das Einzugsgebiet des periodischen Attraktors, dessen Fixpunkt mit einem + markiert ist. Bild 7.5 zeigt den durch die Boxu¨berdeckung angen¨aherten, relativen globalen Attraktor. Es wurden 10 Unterteilungsschritte angewandt, wodurch der gezeigte Zustandsraumausschnitt in vertikaler und horizontaler Richtung in 2 5 = 32 Boxen unterteilt ist. Fu¨r die Boxabbildungen wurde jeweils ein Gitter von 20 × 20 Testpunkten in die Boxen gelegt. Der globale Attraktor zeigt, dass sich mit den eingestellten Parametern Geschwindigkeiten u¨ber 4.2 m/s im Langzeitverhalten nie beobachten lassen. Die Verteilung der Aufenthaltswahrscheinlichkeiten wurde wieder farbkodiert, und es zeigt sich eine Dominanz von Zust¨anden mit hoher Kugelabfluggeschwindigkeit bei einer Phase um 3 / 2 . Bei der Suche nach L¨osungen der speziellen Eigenwertgleichung(2.9) y = P y l¨asst sich genau ein Eigenvektor y 1 finden. Er ist in Bild 7.5 durch die Farbwerte dargestellt. Es kann u¨ber Trajektorienberechnung fu¨r den Parametersatz(7.5) gezeigt werden, dass einer der Fixpunkte stabil ist. Demzufolge ist zu erwarten, dass es einen zweiten Eigenvektor y 2 gibt, der der Eigenwertgleichung genu¨gt. 92

7.1. ANALYSE DES STOSS-MODELLS 2. ORDNUNG Zur Erkl¨arung, warum in dem hier gezeigten Fall kein zweiter Eigenvektor gefunden werden kann, mu¨ssen wir das Einzugsgebiet der periodischen L¨osung betrachten. Ein Verfahren zur Berechnung von Einzugsgebieten wird im folgenden Abschnitt erl¨autert. In Bild 7.5 ist das Einzugsgebiet der periodischen L¨osung rot u¨ber den Attraktor gelegt. Nehmen wir an, es existiere ein Eigenvektor y ¯ 2 , der y ¯ 2 = P y ¯ 2 erfu¨llt und dessen Wahrscheinlichkeitsverteilung einer stabilen, periodischen L¨osung entspricht. Nehmen wir weiter an, der Fixpunkt der periodischen L¨osung bef¨ande sich innerhalb der Box Nr. j . In guter N¨aherung w¨are folglich zu erwarten, dass y 2 ,j = 1, denn s¨amtliche(oder die meisten) Testpunkte aus der Box Nr. j werden erwartungsgem¨aß in Box Nr. j hinein abbilden. Anders ausgedru¨ckt erwartet man, dass das Bild der Box Nr. j innerhalb dieser Box selbst liegt. Der Eintrag p j,j auf der Hauptdiagonalen in der Zeile j der U¨ bergangswahrscheinlichkeitsmatrix w¨are dann p j,j = 100%. Voraussetzung fu¨r die getroffenen Annahmen der Selbstabbildung ist, dass die Box mit dem Fixpunkt vollst¨andig innerhalb des Einzugsgebietes der periodischen L¨osung liegt. Der stabile Fixpunkt ist in Bild 7.5 durch ein+ markiert. Das rot dargestellte Einzugsgebiet der periodischen L¨osung besitzt fu¨r den aktuellen Parametersatz eine filigrane Struktur, die sich als du¨nne F¨aden beschreiben l¨asst, die sich spiralf¨ormig um den Zylinderzustandsraum herum winden und dabei wiederholt teilen. Das Einzugsgebiet ist im Bereich des Fixpunktes so du¨nn, dass es einen großen Teil der Box mit dem Fixpunkt nicht abdeckt. Infolgedessen liegt ein Großteil der Bildpunkte aus dieser Box innerhalb des chaotischen Attraktors, weshalb y 2 ,j << 1 gilt und kein Eigenwert y ¯ 2 existiert. Der globale Attraktor muss in einem solchen Fall durch h¨aufigere Boxunterteilungen in wesentlich feinere Boxen unterteilt werden, um die zweite Eigenl¨osung finden zu k¨onnen. Zur globalen Analyse des Stoßsystems 2. Ordnung wurden relative globale Attraktoren berechnet. Fu¨r Parameters¨atze, bei denen innerhalb des berechneten Attraktors zwei L¨osungen-- eine chaotische und eine periodische Bewegung-- existieren, wurde eine Boxu¨berdeckung des Attraktors gezeigt. Durch Berechnen der Aufenthaltswahrscheinlichkeitsverteilungen fu¨r die Boxmenge des Attraktors wurden in einem Fall beiden L¨osungen statistische Informationen hinzugefu¨gt, im anderen Fall wurde gezeigt, dass sehr hohe Boxaufl¨osungen notwendig sein k¨onnen. Aus der statistischen Verteilung der Zust¨ande auf dem chaotischen Attraktor ließ sich eine Grenzgeschwindigkeit ablesen, oberhalb derer das System im station¨aren Verhalten nie existieren wird. Die Geschwindigkeiten der periodischen Punkte sind gr¨oßer als die Grenzgeschwindigkeit der chaotischen Attraktoren. Hieraus resultiert fu¨r die Praxis, dass eine Betriebsart des Bohrsystems mit periodischem Stoßverhalten fu¨r eine große Bohrrate anstrebenswert ist. Die in der Praxis beobachteten Bohrerfolge [Bar-Cohen et al. 2001a] lassen sich aus der Wahrscheinlichkeitsverteilung der Zust¨ande auf dem chaotischen Attraktor erkl¨aren: die Wahrscheinlichkeit fu¨r Zust¨ande von großer Stoßmassengeschwindigkeit(woraus eine hohe Abbaurate folgt) ist wesentlich h¨oher als fu¨r die u¨brigen Bereiche des At93

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN traktors. 7.1.5 Analyseergebnisse aus Einzugsgebieten Im vorigen Abschnitt studierten wir das L¨osungsverhalten zu Parameters¨atzen, bei denen zwei L¨osungen nebeneinander existieren. Offen war, unter welchen Bedingungen das System welche L¨osung einnimmt. Von Interesse ist ferner, welche Wahrscheinlichkeit dem Erreichen der periodischen und welche der chaotischen L¨osung zuzuordnen ist, wenn man das Schwingstoßsystem unter nicht festgelegten oder zuf¨alligen Bedingungen betreibt. Um auf die Fragestellungen einzugehen, wollen wir die Einzugsgebiete 6 der jeweiligen L¨osungen ermitteln. In den vorigen drei Bildern wurden die Einzugsgebiete in roter Farbe in die Zustandsdiagramme mit hineingezeichnet. Verantwortlich dafu¨r, nach welcher der m¨oglichen stabilen L¨osungen sich das Systemverhalten richtet, ist alleine die Anfangsbedingung, bei der das System startet, bzw. die das System als Trajektorie durchl¨auft. Das Einzugsgebiet einer L¨osung beschreibt genau jenes Gebiet von Anfangsbedingungen im Zustandsraum, von denen aus das System in dPiaersaemLe¨otseursnagtz"ezinwgeeizokgoeenx"istwieirredn.dDe eLr ¨oZsuusntgaenndsbraeusimtzte,inl¨aessstSyssitcehmisn, dessen genau zwei zusammenh¨angende Gebiete unterteilen. Einige transiente Trajektorien, die von Startzust¨anden aus dem Einzugsgebiet nahe um den Fixpunkt herum ausgehen, sind in Bild 7.6 zu sehen. Sie zeigen auch die Entstehung des Punktes an der Parameterstelle h = 5 . 00066 mm am rechten Rand im Bifurkationsdiagramm Bild 6.13. Dort war zu sehen, dass bei nur minimal gr¨oßerem Abstand h die periodische L¨osung abrupt verschwindet und ganz der chaotischen Bewegung weicht. Um das Einzugsgebiet etwa der periodischen L¨osung zu ermitteln, k¨onnte man einen großen Teil des gesamten Zustandsraumes in kleine Unterboxen einteilen(diskretisieren) und fu¨r jede der Unterboxen bestimmen, ob sie Teil des Einzugsgebietes ist. Hierzu k¨onnte man eine oder mehrere Trajektorien mit Startbedingungen aus der Unterbox heraus iterieren und jeweils ihren Verlauf untersuchen. N¨ahert sich die Mehrheit der Trajektorien aus einer Unterbox asymptotisch dem Fixpunkt der periodischen L¨osung, so weist man diese Unterbox dem Einzugsgebiet der periodischen L¨osung zu. Die verbleibenden Unterboxen sind folglich dem Einzugsgebiet der chaotischen L¨osung zuzurechnen. Das geschilderte Verfahren ist ein m¨ogliches, jedoch fu¨r hohe Detailgenauigkeit ein sehr rechenaufw¨andiges. Auf der Grundlage des oben beschriebenen Boxunterteilungsalgorithmus der mengenorientierten Methoden m anifold") wurde der Fortsetzungsalgorithmus entwickelt, der innerhalb der Software g G u a m io ("im g lpobleaml en u tniesrttabislet 6 Einzugsgebiet: engl.: basin of attraction oder domain of attraction 94

7.1. ANALYSE DES STOSS-MODELLS 2. ORDNUNG x 10 -10 5 [m/s] Fixpunkt 0 V-V -5 -1-0.5 0 0.5 1 Fixpunkt [ rad] x 10 -8 Bild 7.6: Transiente Trajektorien fu¨r 1-periodische L¨osung von verschiedenen Anfangszust¨anden des 1-periodischen Einzugsgebietes ausgehend [Dellnitz, Froyland, Junge 2001]. Er n¨ahert das Einzugsgebiet als instabile Mannigfaltigkeit per Ru¨ckw¨artsiteration durch eine Boxmenge an. Fu¨r die Ru¨ckw¨artsiteration der Abbildungsvorschrift ist das Aufstellen der Umkehrabbildung T 1 notwendig, was nicht in allen F¨allen oder nur unter großem Aufwand m¨oglich ist. Zum Berechnen der instabilen Mannigfaltigkeit wird der betrachtete Zustandsraum zun¨achst in Zustandsboxen diskretisiert. Sodann fu¨gt man diejenige Box, die den Fixpunkt enth¨alt, zur bislang noch leeren Boxmenge des Einzugsgebietes hinzu. Ausgehend von dieser Box wird die Mannigfaltigkeit durch Fortsetzung bestimmt. Dies geschieht in zwei Schritten, die wiederholt ausgefu¨hrt werden: 1. Ru¨ckw¨artsabbildung derjenigen Boxen, die im letzten Schritt neu zur Boxmenge hinzugekommen sind. 2. Hinzufu¨gen derjenigen Boxen zur Boxmenge, in denen wenigstens ein Bildpunkt eines Testpunktes aus einer der Abbildungen von Schritt 1 liegt und die noch nicht der Boxmenge angeh¨oren. Der Algorithmus stoppt selbst¨andig, sobald keine neuen Boxen mehr hinzugefu¨gt wurden. Bild 7.7 zeigt die auf diese Weise approximierten Einzugsgebiete fu¨r die periodische(rotes Gebiet) und die chaotische L¨osung(weißes Gebiet). Stellt man sich das Einzugsgebiet um den zylinderf¨ormig dargestellten Zustandsraum(Bild 7.1) herum gelegt vor, so f¨allt auf, dass beide Zustandsgebiete aus einem zusammenh¨angenden Gebiet bestehen. Zur leichteren Betrachtung der zylindrischen Fl¨ache des Zustandsraumes in der Papierebene ist das Zustandsgebiet in der Abbildung u¨ber drei Zyklen fortgesetzt dargestellt. Wie bereits die chaotischen Attraktoren des hier untersuchten Systems nach genu¨gend vielen Boxunterteilungen zeigen, besitzt auch das Einzugs95

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN gebiet der periodischen L¨osung Zu¨ge fraktaler Gebilde, wenn man auf die selbst¨ahnliche Fortsetzung schaut: der Beginn der Entstehung des Gebietes Bild 7.7: Einzugsgebiet zu chaotischer(weiß) und zu periodischer(rot) L¨osung liegt im Bereich des Fixpunkts. Das rotgef¨arbte Einzugsgebiet entwickelt sich in positiver Drehrichtung spiralf¨ormig um den Zustandsraumzylinder herum weiter. Fu¨r ansteigende Phase verzweigt es sich fortw¨ahrend in zwei A¨ ste, von denen der obere Ast gr¨oßerer Abfluggeschwindigkeit wenig sp¨ater endet und der untere Ast sich wiederholt verzweigt. Um auf die Frage der Wahrscheinlichkeit des Erreichens der einen oder anderen L¨osung zu antworten, nehmen wir uns das Einzugsgebiet zu Hilfe: die Wahrscheinlichkeit, mit der die periodische L¨osung bei zuf¨alliger Anfangsbedingung erreicht wird, gleicht dem Fl¨achenanteil, den das dunkle 96

7.2. ANALYSE DES STOSS-MODELLS 4. ORDNUNG Einzugsgebiet am gesamten Zustandsraum hat. Da sich dieser in positive V Richtung des Zustandsraumes ins Unendliche ausdehnt, ist ein begrenzter Bereich zu w¨ahlen, innerhalb dessen die fl¨achenm¨aßige Aufteilung der Einzugsgebiete gewertet wird. Fu¨r das in Bild 7.7 dargestellte Einzugsgebiet ist die Wahrscheinlichkeit, von einem rot gef¨arbten Zustand aus zu starten, 7.89%, wohingegen ein Anfangszustand mit komplement¨aren 92.11% im hellen Teil beginnen wird. Die hohe Wahrscheinlichkeit, bei nicht exakt steuerbarem Ausgangszustand in den chaotischen Attraktor hineinzulaufen, liefert eine weitere Erkl¨arung dafu¨r, dass in der Praxis nie regul¨are Bewegungen beobachtet wurden. Die Ermittlung der Einzugsgebiete der einzelnen, koexistierenden Attraktoren ist ein wichtiger Bestandteil einer umfassenden Systemanalyse. Sie liefert ein Maß fu¨r die Wahrscheinlichkeiten, mit der die einzelnen Attraktoren bei einem Iterationsstart mit zuf¨alligen Anfangsbedingungen erreicht werden. Es wurde gezeigt, dass die Einzugsgebietsberechnung fu¨r Schwingstoßsysteme sehr effizient durch die Anwendung eines BoxFortsetzungsalgorithmus durchgefu¨hrt werden kann. 7.2 Analyse des Stoß-Modells 4. Ordnung Dieses Unterkapitel enth¨alt die Untersuchung des Modells 4. Ordnung mit den mengenorientierten numerischen Methoden. Im ersten Abschnitt wird wieder eine Startbox im Zustandsraum gew¨ahlt und ihr relativer globaler Attraktor bestimmt. Der zweite Abschnitt bespricht die Analyseergebnisse, die aus der Bildung der Attraktoren abgeleitet werden k¨onnen. 7.2.1 Approximation des relativen globalen Attraktors Wie fu¨r das zuvor behandelte, einfachere Modell sollen auch anhand des deutlich komplexeren Modells 4. Ordnung die mengenorientierten Methoden demonstriert werden, um Informationen u¨ber die globale Dynamik des Systems zu gewinnen. Wieder wird der in GAIO implementierte rga Algorithmus(Boxunterteilungsalgorithmus) zur Bestimmung relativer globaler Attraktoren auf die in C-Code u¨bersetzte Modelldatei angewendet. Wir verwenden fu¨r s¨amtliche Analysen innerhalb dieses Unterkapitels die selben System- und Prozessparameter, die in Unterkapitel 6.4 benutzt wurden: Oszillator­ Ampl. Frequ. u 0 f 15 µ m 20.19 kHz Stoßzahl an Osz. Bohrer 1 2 0.5 0.6 Massen­ verh¨altn. M/m = µ 264. 7 g 2. 0 g = 132 Beschleunigung aus Anpresskraft Gravitat. F A /M g 3N 0. 2647 kg = 11 m/s 2 9.81 m/s 2 97

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN Bild 7.8: Projektion der Boxu¨berdeckung des relativen globalen Attraktors A Q nach 21 Unterteilungsschritten in den y - x- x Raum mit Farbkodierung der Aufenthaltswahrscheinlichkeiten in[%]. Fu¨r die schrittweisen Boxunterteilungen innerhalb des rga -Algorithmus wurden fu¨r alle in diesem Unterkapitel gezeigten Rechnungen die Unterteilungsrichtungen zyklisch durch alle 4 Zustandsdimensionen gewechselt in der Reihenfolge y  x x  . Zur approximativen Abbildung der 4-dimensionalen Boxen wurden jeweils 3 4 = 81 Testpunkte mit innerhalb der Boxen einschließlich ihres Randes dzuerfalOlspvteirotneil"t M . o D n i t e e s C e a T r e l s o t"punktanzahl zeigte sich als hinreichend groß fu¨r die Ann¨aherung des dynamischen Verhaltens. Gleichzeitig blieb die Rechenzeit in einem akzeptablen Rahmen. Die gew¨ahlte, 4-dimensionale Startbox Q wurde in allen Bildern durch einen roten Rahmen angedeutet. In der Phasendimension u¨berdeckte sie stets den gesamten Zustandsbereich [0 , 2 ]. Zur Erl¨auterung der Analyse eines Systems 4. Ordnung werden auf diesen Seiten exemplarisch einige Ergebnisse dargestellt: Bild 7.8 zeigt die Ausdehnung der Boxmenge nach 21 Unterteilungsschritten. Durch die zyklische Permutation der Unterteilungsrichtungen wurden die Boxen somit 6-mal in -Richtung und 5-mal in den 3 u¨brigen Dimensionen zweigeteilt. In Bild 7.8 ist ersichtlich, dass die Boxmenge den Rand der Startbox Q in negative x -Richtung und in negative y -Richtung tangiert. Es scheint naheliegend zu sein, die Startbox in diesem Beispiel gr¨oßer zu w¨ahlen, um die globale Dynamik zu beschreiben. Jedoch liegt die Ursache fu¨r die Randberu¨hrung in der langsamen Konvergenz dieses Systems. Es gibt mehrere M¨oglichkeiten, um solch einem Systemverhalten zu begegnen. Das Vorgehen, eine wesentlich gr¨oßere Anzahl an Unterteilungsschritten zu durchlaufen, ist aus praktischen Gru¨nden begrenzt. Durch die 98

7.2. ANALYSE DES STOSS-MODELLS 4. ORDNUNG langsame Konvergenz fallen die transienten Boxen nur sehr langsam weg und mu¨ssen sehr oft mitunterteilt werden. Dies fu¨hrt zu sehr hohen Boxzahlen und hohem Speicherbedarf sowie großer Rechenzeit. Alternativ w¨are es denkbar, eine adaptive Boxunterteilungsstrategie anzuwenden: dabei wird das invariante Maß berechnet und Boxen mit einem Maß unter dem Durchschnittswert aller Boxmaße werden nicht unterteilt bzw. Boxen, deren Maß unter einem gew¨ahlten Schwellenwert liegt, werden aus der Menge gel¨oscht. Eine weitere M¨oglichkeit, mit langsamer Konvergenz umzugehen, ist die Folgende. Anstatt fu¨r die Iterationsvorschrift die einfache Abbildungsvorschrift T anzuwenden, kann die Iterationsschrittweite auf einen Wert gr¨oßer Eins gesetzt werden. Fu¨r das betrachtete System wurde die Schrittweite auf 5 gesetzt. Somit wurde die Abbildung T 5 betrachtet, die eine schnellere Konvergenz zum invarianten Maß besitzt als die urspru¨ngliche Abbildung T 1 . Die Bilder 7.11 und 7.12 zeigen eine Attraktoru¨berdeckung nach 22 Unterteilungsschritten fu¨r die Abbildung T 5 desselben Systems mit denselben Parametern. Der folgende Abschnitt erl¨autert dieses Vorgehen. Bild 7.9: Ausblenden transienter Boxen des Attraktors aus Bild 7.8. 7.2.2 Analyseergebnisse und Diskussion Im vorhergehenden Abschnitt wurde die Berechnung von invarianten Mengen, den relativen globalen Attraktoren des Modells 4. Ordnung, erl¨autert. Da das System im 4-dimensionalen Zustandsraum definiert ist, besitzt auch die Gestalt der Attraktoren 4 Dimensionen. Fu¨r die grafische Darstellung der Attraktoren auf einem 2-dimensionalen Medium(Papier, Bildschirm) oder 99

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN 3-dimensional perspektivisch muss der Attraktor zwei- bzw. einmal entlang einer Zustandsgr¨oße projiziert werden. Dies erschwert die Auswertung, da bei der Betrachtung der projizierten Gestalt Informationen vernachl¨assigt werden mu¨ssen. Daher muss u¨berlegt werden, entlang welcher Achse projiziert werden soll. Bild 7.10: 2D-Ansicht der Attraktordarstellung aus Bild 7.9, farbkodierte Aufenthaltswahrscheinlichkeiten in[%]. Zum Erkennen der dynamischen Eigenschaften des Bohrer-Modells ist die zeitliche Information, die Auskunft u¨ber die Phasenlage k der Stoßkontakte gibt, im Vergleich zu den drei anderen Zustandsgr¨oßen entbehrlich. Integriert man die 4-dimensionalen Boxen des Attraktors u¨ber eine Periodendauer, erh¨alt man eine 3-dimensionale Abbildung des Attraktors. Bild 7.8 zeigt das Resultat dieses Vorgehens. Die Boxmenge approximiert einen Bereich im Zustandsraum. Zur Quantifizierung des dynamischen Verhaltens auf der Boxmenge wurde wie fu¨r das System 2. Ordnung die U¨ bergangswahrscheinlichkeitsmatrix berechnet und deren Eigenvektor zum betragsm¨aßig gr¨oßten Eigenwert, der fast Eins ist, numerisch berechnet. Nach einer Normierung, so dass die Summe aller Eigenvektoreintr¨age Eins ergibt, stellen die Eigenvektorelemente die Aufenthaltswahrscheinlichkeiten des Systems fu¨r jede Box dar. In den dargestellten Boxu¨berdeckungen innerhalb dieses Unterkapitels sind die Aufenthaltswahrscheinlichkeitswerte durch eine Farbkodierung den Boxen zugeordnet, wobei die Skalenwerte auf den Farbskalen stets in der Einheit [%] angegeben sind. Durch die logarithmische Skalierung der Wahrscheinlichkeitswerte auf der Farbskala werden die um einige Gr¨oßenordnungen unterschiedlichen Wahrscheinlichkeiten deutlich. Blaue bis gru¨ne Farbt¨one stehen fu¨r sehr kleine Wahrscheinlichkeiten von etwa 10 25 % bis 10 10 %. Da Langzeititerationen zeigen, dass Geschwindigkeitsbetr¨age von Oszillator bzw. Kugel gr¨oßer 3 m/s bzw. 10 m/s nicht zu erwarten sind, zeigt Bild 7.9 100

7.2. ANALYSE DES STOSS-MODELLS 4. ORDNUNG die gleiche Boxu¨berdeckung wie in Bild 7.8, jedoch sind hier nur Boxen mit Wahrscheinlichkeiten u¨ber dem Schwellenwert von 10 10 % dargestellt. Die verbleibenden Boxen stellen eine bessere Ann¨aherung an die invariante Mannigfaltigkeit des Systems, die die Langzeitdynamik beschreibt, dar. In der Seitenansicht von Bild 7.9, die in Bild 7.10 gezeigt ist, l¨asst sich erkennen, dass sich im station¨aren Systemverhalten nie Kugel- oder Oszillatorabfluggeschwindigkeiten gr¨oßer 2 m/s in positive Richtung(sich vom Bohrstab entfernend) ergeben. In Richtung auf den Bohrstab(in negative Koordinatenrichtung) liegen die maximalen Geschwindigkeiten im Bereich um 8 m/s bzw. 2.5 m/s. Aus den maximalen Geschwindigkeiten lassen sich fu¨r das betrachtete Anwendungsbeispiel des Stoßbohrsystems maximale Stoßenergien ableiten und u¨ber das invariante Maß gewichten. Im vorigen Abschnitt wurde auf die Tangierung des Startboxrandes durch die Boxmenge eingegangen. Ausblenden der Boxen, die aufgrund ihres Wahrscheinlichkeitsmaßes der transienten Menge zugeordnet werden k¨onnen, war eine Maßnahme, um das station¨are Verhalten der langsam konvergenten Menge anzun¨ahern. Die zweite, hier umgesetzte Maßnahme ist die Untersuchung einer mehrfach Iterierten der einfachen Abbildung. Bei fu¨nffacher Iteration der Testpunkte bildet die in den Bildern 7.11 und 7.12 dargestellte Boxu¨berdeckung eine Ann¨aherung an die invariante Menge. Die Boxmenge in Bild 7.11 zeigt einen sehr schmalen Bereich um große, negative Geschwindigkeiten herum. Es l¨asst sich daraus ablesen, dass hohe Geschwindigkeiten sowohl von der Kugel als auch vom Oszillatorschwerpunkt nur in direkter N¨ahe zum Bohrstab(bei minimaler Position x ) auftreten k¨onnen. Die gr¨oßten auftretenden Abst¨ande x von Kugel und Oszillator vom Bohrstab im Stoßaugenblick liegen im Bereich von 2 cm. Kleinere Flugh¨ohen sind jedoch um Gr¨oßenordnungen wahrscheinlicher. Die Ursache der abgeflachten Seite der Attraktormenge in Richtung - x liegt in der Bedingung x k - u 0 cos t k > 0, nach der die Kugel und der Oszillator nicht in den Bohrstab eindringen k¨onnen. Diese Bedingung stellt eine Grenzebene 1 in diesem projizierten Zustandsraum dar: 1 = { ( y , x, x  ) | x = u 0 }. Eine weitere, schr¨ag liegende Ebene 2 begrenzt die Boxmenge in Richtung + y  und - x  durch die Nichteindringbedingung von Kugel und Oszillator direkt nach einem Stoß: 2 = { ( y , x, x  ) | x - u 0  = y }. Bild 7.12 zeigt die Projektion derselben Boxmenge wie in Bild 7.11 in den y x -Raum. Im Unterschied zum Ergebnis der chaotischen Bewegung des Systems 2. Ordnung l¨asst sich hier ablesen, dass keine Phasenlage ein pr¨adestinierter Stoßzustand ist. In dieser Ansicht l¨asst sich wieder die Information gewinnen, dass Wahrscheinlichkeiten fu¨r das Auftreten kleinerer Geschwindigkeiten von Kugel und Oszillator hoch sind. Die Boxmenge in dieser Projektion wird wegen der Nichteindringbedingung nach unten durch 101

KAPITEL 7. ANALYSE MIT MENGENORIENTIERTEN METHODEN Bild 7.11: Projektion des relativen globalen Attraktors A Q der fu¨nffach iterierten Abbildung T 5 in den y - x- x Raum. Blau-gru¨ne Einf¨arbungen der Boxen entsprechen kleiner, rot-braune Farben entsprechen großer Aufenthaltswahrscheinlichkeit des Systems in der jeweiligen Box. die Fl¨ache 3 = { ( y , x, x  ) | x  + u 0  sin t = y } begrenzt, was die sinuswellenf¨ormige Gestalt der Boxmenge begru¨ndet. Zust¨ande jenseits der Grenzebenen 1 , 2 oder 3 k¨onnen physikalisch nicht auftreten. Die Berechnung von relativen globalen Attraktoren des Systems 4. Ordnung erm¨oglicht die Identifizierung struktureller Systemeigenschaften des Langzeitverhaltens. Statistische Eigenschaften k¨onnen als invariantes Maß auch fu¨r h¨oherdimensionale Systeme aus Eigenvektoren berechnet werden. Sie stellen eine weitere Information zur Bewertung des Attraktorgebietes auf der Boxmenge dar. In Kombination mit adaptiven Algorithmen k¨onnen Verfahren angewendet werden, um die Approximation invarianter Mengen langsam konvergenter Systeme zu beschleunigen. 102

7.2. ANALYSE DES STOSS-MODELLS 4. ORDNUNG Bild 7.12: Projektion des relativen globalen Attraktors A Q der fu¨nffach iterierten Abbildung T 5 in den y x -Raum. Blaue Boxen entsprechen einem kleinen, rot-braune Boxen einem großen Wahrscheinlichkeitsmaß. 103

.

Kapitel 8 Zusammenfassung, Diskussion und Ausblick 8.1 Zusammenfassung In der vorliegenden Arbeit stand die Analyse nichtlinearer dynamischer Systeme-- im speziellen die Beschreibung und Analyse nichtglatter Systemdynamik-- im Mittelpunkt. Fu¨r diese Aufgabe, die fu¨r Ingenieure aller Fachrichtungen ein wichtiges Thema darstellt, sollte der Einsatz einer neuartigen uMnedthdoedmenoknlsatsrsieer,tdwerer"dmenen. gDenieosreieMnteiethrtoednennuhmaettreisnchseinchMbeetrheoitdsevn"ie,lfgaecphru¨imft Einsatz auf nichtlineare Systeme bew¨ahrt. Fu¨r die Klasse der nichtglatten dynamischen Systeme war in dieser Arbeit zu untersuchen, inwiefern sie sich fu¨r deren Analyse eignen. Exemplarisch fu¨r eine Vielzahl von Schwingstoßsystemen wurde eine vor wenigen Jahren entstandene Entwicklung eines Stoßbohrsystems fu¨r spr¨odharte Stoffe wie Gestein herangezogen. Dieses System zeichnet sich durch seinen einfachen Aufbau aus und wird dennoch von einer hochkomplexen Dynamik beherrscht. Weiterhin motivierte die Erlangung von tiefergehendem Verst¨andnis des dynamischen Verhaltens des Bohrprozesses die Untersuchung dieses Systems, welches in der Literatur bisher nur in wenigen Arbeiten behandelt wurde. Zur mathematischen Beschreibung des Systemverhaltens waren geeignete Modelle zu bilden. Zur Validierung der theoretischen Simulationen waren Pru¨fst¨ande mit Laborprototypen aufzubauen und mittels geeigneter Messverfahren Zeitreihen der Stoßmassendynamik aufzunehmen. Zentrale Komponente dieses Systems ist eine Stoßmasse, die mit Frequenzen im Schallbereich zwischen einem Bohrstab und einem UltraschallOszillator Flugbewegungen ausfu¨hrt. Die zum Gesteinsbohren n¨otige Stoßenergie wird von der Stoßmasse in ihrer Intensit¨at und Frequenz transformiert. Im Zentrum der Untersuchung stand die Bewegungsanalyse dieser Stoßmasse. Fu¨r das Erreichen der Ziele wurde das Stoßbohrsystem in zwei Schritten untersucht: zun¨achst wurde die Dynamik des Systems in ihrer Komplexit¨at reduziert, indem die Annahme getroffen wurde, dass nur die Stoß105

KAPITEL 8. ZUSAMMENFASSUNG, DISKUSSION UND AUSBLICK masse Bewegungen ausfu¨hren kann. Das reduzierte System wurde aufgrund der Stoßkontakte als zeitdiskretes dynamisches System mathematisch modelliert. Differenzengleichungen 2. Ordnung wurden hergeleitet und ein spezieller Algorithmus zur L¨osung der impliziten, transzendenten Modellgleichungen aufgestellt. In einem zweiten Schritt wurde das Modell verfeinert, um die Dynamik des beweglichen Ultraschalloszillators mit einzubeziehen. Dies fu¨hrte auf ein verfeinertes Differenzengleichungssytem 4. Ordnung. Das Verhalten beider Systeme wurde sowohl in der Praxis als auch in der Theorie untersucht. Ein Laborpru¨fstand wurde fu¨r das reduzierte System entwickelt und aufgebaut. Mit verschiedenen beru¨hrungslosen Messverfahren wurden Zeitreihen des realen Verhaltens aufgenommen. Techniken aus dem Gebiet der nichtlinearen Dynamik wurden angewandt sowie Fixpunkte und periodische L¨osungen bestimmt, in Zeitreihen wurden periodische mit chaotischen L¨osungen verglichen, Bifurkationsanalysen stabiler periodischer L¨osungen wurden fu¨r verschiedene Systemparameter durchgefu¨hrt, und ein Einzugsgebiet der periodischen L¨osung wurde berechnet. Best¨atigt und erg¨anzt durch diese Grunduntersuchungen wurden mithilfe der mengenorientierten numerischen Methoden relative globale Attraktoren und Aufenthaltswahrscheinlichkeitsverteilungen auf diesen Attraktoren bestimmt. Aufbauend auf einem Boxunterteilungsalgorithmus, auf dem die mengenorientierten Methoden basieren, wurde ein Algorithmus zur Berechnung von Einzugsgebieten der periodischen L¨osungen implementiert, mit dem instabile Mannigfaltigkeiten der Umkehrabbildung des Systems iteriert wurden. 8.2 Diskussion der Ergebnisse Die Anwendung der mengenorientierten numerischen Methoden hat verdeutlicht, dass diese Techniken zur Analyse nichtglatter Systeme der Stoßdynamik eine thoden" wichtige und sinnvolle Erg¨anzung bilden. Es hat sich gezeigt, dass ddieer mbeisnhgeernigoernien"tkilearstseinschAenns¨aMtzeedafu¨r keiner Anpassung bedurften. Allerdings ist eine besondere, zeitdiskrete Modellierungsweise, die dem nichtglatten Systemcharakter gerecht wird, n¨otig. Sofern die Beschreibung der Bewegungsabl¨aufe mit Differenzengleichungssystemen m¨oglich ist, kann der mengenorientierte Ansatz angewendet werden, ohne Anpassungen fu¨r Stoßsysteme vornehmen zu mu¨ssen. Gezeigt wurde die Bestimmung relativer Attraktoren, die ein komplettes Bild der gesamten Systemdynamik liefern, ohne dass eine Vorkenntnis u¨ber das Systemverhalten n¨otig ist. Der Nachweis von Aufenthaltswahrscheinlichkeiten u¨ber die Attraktorenzust¨ande war fu¨r nichtglatte Systeme m¨oglich. Fixpunkte, die durch Attraktoren lokalisiert werden konnten, wurden fu¨r eine Robustheitsund Parameterstudie einer Verzweigungsanalyse unterzogen. Dabei hat sich gezeigt, dass periodische L¨osungen nur auf einem sehr schmalen Parameterintervall existieren und durch Periodenverdoppelungsbifurkationen ins Chaos u¨bergehen. Besonderen Wert zeigen die Methoden beim Ermitteln des station¨aren Langzeitverhaltens. Die exponentielle Divergenz benachbarter L¨osungen 106

8.2. DISKUSSION DER ERGEBNISSE erschwert aufgrund der begrenzten Rechenpr¨azision den iterativen Einsatz von Integrationsmethoden. Hier liefern die mengenorientierten Methoden einen wichten Beitrag, da sie das Langzeitverhalten ausschließlich durch Ein-Schritt-Integrationen(im Fall langsamer Konvergenz diskreter Systeme kann die Iterationsschrittweite geringfu¨gig erh¨oht werden) ermitteln. Fu¨r das konkrete Beispiel eines Ultraschall-Stoßbohrers mit eindimensionaler Bewegung einer Stoßmasse und doppelseitigem Stoßkontakt wurde eine Modellierungstechnik dargestellt, die allein auf Differenzengleichungen beruht. Die Technik ist auch anwendbar auf ¨ahnliche Systeme, in denen es sinnvoll ist, von Stoß zu Stoß zu iterieren. Beispielsweise resultierte aus der Analyse, dass in der station¨aren, chaotischen Bewegung nur Geschwindigkeiten unterhalb eines Grenzbereichs erreicht wurden. Die Betrachtung der statistischen Ergebnisse offenbarte, dass Zust¨ande hoher Geschwindigkeit im System 2. Ordnung dominieren und tr¨agt somit zur Erkl¨arung des Funktionsprinzips des Bohrsystems bei. Periodische L¨osungen konnten nur am reduzierten Modell 2. Ordnung gezeigt werden. Sie besaßen fu¨r alle untersuchten Parameterkombinationen eine gr¨oßere Energieu¨bertragungsleistung als chaotische L¨osungen, was aus den h¨oheren Abfluggeschwindigkeiten der Stoßmasse bei periodischem Verhalten hervorgeht. In der Analyse dieses Bohrsystems hat sich gezeigt, dass eine große Sensibilit¨at bezu¨glich der Systemparameter besteht. Einige Parameterbereiche konnten gefunden werden, innerhalb derer verschiedene L¨osungen nebeneinander koexistieren k¨onnen. So wurden mehrfach-periodische L¨osungen neben chaotischen nachgewiesen. Eine große Streuung der Stoßzahlen, die in Versuchen beobachtet wurde, trug dazu bei, dass periodische L¨osungen in der Praxis nicht auftreten konnten. Der Grund hierfu¨r liegt in der starken Abh¨angigkeit des L¨osungstypus(irregul¨ar/regul¨ar) von den Stoßzahlen. Weiterhin zeigte die Berechnung von Einzugsgebieten, dass bei zuf¨alligen Anfangsbedingungen die Wahrscheinlichkeit fu¨r chaotische Bewegung sehr groß ist. Die Auswertung experimenteller Daten best¨atigte eine Hypothese, die der Modellbildung zugrunde lag. Es wurde angenommen, dass die Resonanzl¨angsschwingungen des Oszillators vor jedem Stoß ihre maximale Resonanzamplitude besitzen. Die beim Stoß u¨bertragene Schwingungsenergie war in einem Bruchteil der Stoßmassenflugzeit wieder der Schwingungsform zugefu¨hrt worden, so dass sich die station¨are Schwingung wieder vor dem n¨achsten Stoß einstellen konnte. Die Analyse des verfeinerten Systems 4. Ordnung resultierte in einer begrenzten Boxmenge, aus der die transienten Boxen durch Berechnung eines invarianten Maßes herausgefiltert wurden. So konnten die m¨oglichen Systemzust¨ande unter der chaotisch-deterministischen Langzeitdynamik als station¨are, invariante Menge extrahiert werden. 107

KAPITEL 8. ZUSAMMENFASSUNG, DISKUSSION UND AUSBLICK 8.3 Ausblick Die mengenorientierten numerischen Methoden liefern einen wichtigen Beitrag nicht nur fu¨r die Systemanalyse nichtlinearer Systeme, sondern k¨onnen auch das Systemverst¨andnis nichtglatter Systeme verst¨arken. In der vorliegenden Arbeit wurde ein reduziertes System in 2 Dimensionen und ein verfeinertes System in 4 Dimensionen analysiert. Das verfeinerte System beanspruchte deutlich mehr Rechenleistung, was darauf zuru¨ckzufu¨hren ist, dass die beobachteten Mannigfaltigkeiten die volle Systemdimension beibehalten. Fu¨r die vorliegenden Untersuchungen wurden zur U¨ berwindung der langsamen Konvergenz des Systems dessen Mehrfachiterierte untersucht. Zur Anwendung der mengenorientierten Methoden wird es hilfreich sein, das Konvergenzverhalten realer Systeme- insbesondere langsam konvergenter Systeme- unter verschiedenen Kombinationen von Parametern wie Anzahl der Testpunkte, Boxunterteilungstiefen oder Iterationsschrittweiten exemplarisch zu betrachten. Zwei mechanische Ersatzsysteme wurden aus dem UltraschallBohrsystem zur Untersuchung abgeleitet. Zeitreihensimulationen und Poincar´e-Portraits beider Modelle als auch Laboruntersuchungen an realen Systemen wiesen im Einklang mit mengenorientierten Analysen eine starke Dominanz von chaotischem Verhalten nach. Im Zuge einer Optimierung der Bohrfortschrittsgeschwindigkeit w¨are es denkbar, eine Regelung zu entwerfen, die das System in periodischen Zustand fu¨hrt. Auf die Nichtvorhersagbarkeit der Stoßzahlen und die hohe St¨oranf¨alligkeit der periodischen L¨osungen w¨are hierbei zu achten. 108

Literaturverzeichnis [Abraham& Shaw 2000] Ralph H. Abraham, Christopher D. Shaw, 2000: Dynamics-- the geometry of behavior. Part I, Periodic Behavior-4th ed., ISBN 0-942344-18-9, Aerial Press. 21 [Babitsky 1998] Vladimir I. Babitsky(Babickij), 1998: Theory of VibroImpact Systems and Applications. Springer-Verlag, Berlin, Heidelberg. 28 [Badescu et al. 2005] Badescu, M., Bao, X., Bar-Cohen, Y., Chang, Z., Sherrit, S.: Integrated Modeling of the Ultrasonic/Sonic Drill/Corer-- Procedure and Analysis Results. In Proc. SPIE Smart Struct. Conf., San Diego, CA., SPIE Vol. 5764-37, Marc 7-10, 2005. 5, 29 [Bar-Cohen et al. 2001a] Yoseph Bar-Cohen, Stewart Sherrit, Benjamin P. Dolgin, Xiaqi Bao, Zensheu Chang, Dharmendra S. Pal, Ron Krahe, Jason Kroh, Shu Du and Thomas Peterson, 2001: Ultrasonic/sonic drilling/coring(USDC) for planetary applications. In Proceedings of SPIE's 8th Annual International Symposium on Smart Structures and Materials, March 5-8, 2001, Newport. 18, 67, 93 [Bar-Cohen et al. 2001b] Bar-Cohen, Y., Sherrit, S., Dolgin, B., Bridges, N., Bao, X., Chang, Z., Yen, A., Saunders. R.: Ultrasonic/Sonic Driller/Corer(USDC) as a sampler for planetary exploration. In IEEE Aerospace Conference on the topic of Missions, Systems and Instruments for In Situ Sensing(Session 2.05), March 10-17, 2001 Big Sky, Montana. 18, 29 [Brockhaus 1989] Brockhaus, Naturwissenschaften und Technik, Band 3. Brockhaus GmbH, Wiesbaden, 1989. 5 [Brogliato 1996] Bernard Brogliato, 1996: Nonsmooth Impact Mechanics. Models, Dynamics and Control. Springer-Verlag London Limited, 1996. 28, 30, 36 [Dellnitz& Hohmann 1997] Dellnitz, M. and Hohmann, A., 1997: A subdivision algorithm for the computation of unstable manifolds and global attractors. In Numerische Mathematik Vol. 75, No. 3, Seiten 293-317. 8 [Dellnitz& Junge 1999] Michael Dellnitz, Oliver Junge, 1999: On the approximation of complicated dynamical behavior. SIAM Journal on Nu109

LITERATURVERZEICHNIS merical Analysis, Vol. 36, No. 2, Seiten 491-515, Society for Industrial and Applied Mathematics. 89 [Dellnitz, Froyland, Junge 2001] Dellnitz, M., Froyland, G. and Junge, O., 2001: The algorithms behind GAIO-- Set oriented numerical methods for dynamical systems. In B. Fiedler: Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems. Seiten 145­174, SpringerVerlag, Berlin. 8, 10, 12, 84, 95 [Dellnitz& Junge 2002] Michael Dellnitz and Oliver Junge, 2002: Set oriented numerical methods for dynamical systems. In B. Fiedler, G. Iooss and N. Kopell(eds.): Handbook of Dynamical Systems III: Towards Applications. World Scientific. 8, 9, 10, 17 [Goldschmidt 2003] Stefan Goldschmidt, 2003: Anwendung mengenorientierter numerischer Methoden zur Analyse nichtlinearer dynamischer Systeme am Beispiel der Spurfu¨hrungsdynamik von Schienenfahrzeugen. In J¨org Wallaschek(Hrsg.): Heinz Nixdorf InstitutVerlagsschriftenreihe, Band 112. Dissertation Universit¨at Paderborn. 8 [Guckenheimer& Holmes 1983] John Guckenheimer and Philip J. Holmes, 1983: Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer-Verlag, New York. 29, 36, 38 [Holmes 1982] Philip J. Holmes, 1982: The dynamics of repeated impacts with a sinusoidally vibrating table. In Journal of Sound and Vibration, Vol. 84, No. 2, Seiten 173-189, Academic Press Inc., London. 29, 36, 38 [Hsu 1980] Chieh-Su Hsu, 1980: A theory of cell-to-cell mapping dynamical systems. In Journal of Applied Mechanics, Vol. 47: Seiten 931-939. ASME, Transactions; American Society of Mechanical Engineers, Winter Annual Meeting, Chicago, Ill., Nov. 16-21, 1980. 4, 26 [Hsu 1981] Chieh-Su Hsu, 1981: A generalized theory of cell-to-cell mapping for nonlinear dynamical systems. In Journal of Applied Mechanics, Vol. 48: Seiten 634-642. 26 [Hsu 1987] Chieh-Su Hsu, 1987: Cell-to-Cell Mapping. A Method of Global Analysis for Nonlinear Systems In Applied Mathematical Sciences, Vol. 64, Springer-Verlag, 1987. 8, 26 [Hsu 1992] Chieh-Su Hsu, 1992: Global analysis by cell mapping. In International Journal of Bifurcation and Chaos, Vol. 2, No. 4, Seiten 727-771, World Scientific Publishing Company. 26 [Kraker, van der Spek und van Campen 1999] Bram de Kraker, Jeroen A. W. van der Spek and Dick H. van Campen, 1999: Extensions of Cell Mapping for Discontinuous Systems. In Applied Nonlinear Dynamics and Chaos of Mechanical Systems with Discontinuities, Chapter 4, Seiten 61­102, Editors: Marian Wiercigroch, Bram de Kraker; In World 110

LITERATURVERZEICHNIS Scientific Series on Nonlinear Science, Series A, Vol. 28, Series Editor: Leon O. Chua, 2000. 27 [Leine et al. 2000] Remco I. Leine, D. H. van Campen and B. L. van de Vrande, 2000: Bifurcations in Nonlinear Discontinuous Systems. In: Nonlinear Dynamics, Vol. 23, Seiten 105-164, Kluwer Academic Publishers, Niederlande, 2000. 28 [Leine& Nijmeijer 2004] Remco I. Leine and Nijmeijer, H., 2004: Dynamics and Bifurcations in Non-Smooth Mechanical Systems. Lecture Notes in Applied and Computational Mechanics, Vol. 18, Berlin, Heidelberg, New-York, Springer-Verlag. 28 [Lorenz 1963] Edward N. Lorenz, 1963: Deterministic Nonperiodic Flow. In Journal of the Atmospheric Sciences, Vol. 20, No. 2, 130-141, M¨arz 1963. Ref. aus: Wikipedia, Stand 29. Juni 2007. 21 [Mello, Tufillaro 1985] T. M. Mello and N. B. Tufillaro, 1985: Strange attractors of a bouncing ball. In Am. J. Phys. Vol. 55, No. 4, American Association of Physics Teachers, 1987. 29, 36 [Nayfeh 1995] Ali Hasan Nayfeh and Balakumar Balachandran, 1995: Applied nonlinear dynamics. analytical, computational and experimental methods. Wiley Series in Nonlinear Science. John Wiley& Sons, Inc., New York. 23 [Neumann 2002] Nicolai Neumann, 2002: Nichtlineare Effekte bei L¨angsschwingungen axial polarisierter piezokeramischer St¨abe: Experimentelle Untersuchungen und Parameteridentifikation. Diplomarbeit, Institut fu¨r Mechanik der Technischen Universit¨at Darmstadt, 2002. 2 [Neumann& Sattel 2007] Nicolai Neumann und Thomas R. Sattel, 2007: Set-oriented Numerical Analysis of a Vibro-Impact Drilling System with Several Contact Interfaces. In Journal of Sound and Vibration, Special Issue on Vibro-Impact Systems, Elsevier Science. 36 [Neumann et al. 2007] Nicolai Neumann, Thomas Sattel und J¨org Wallaschek, 2007: On set-oriented numerical methods for global analysis of non-smooth mechanical systems. In Journal of Vibration and Control, Special Issue on Mathematical Methods in Engineering. Hrsg.: Ali H. Nayfeh; Guest Editors: Kenan Tas, J. Tenreiro Machado, Dumitru Baleanu. Verlag Sage Publications, 2007. 5, 28, 29 [Norris 1997] J. R. Norris, 1997: Markov Chains. Cambridge University Press, Cambridge, 1997. 15 [Pavlovskaia& Wiercigroch 2003] Ekaterina Pavlovskaia, Marian Wiercigroch, 2003: Analytical drift reconstruction for visco-elastic impact oscillators operating in periodic and chaotic regimes. In Chaos, Solitons & Fractals Vol. 19, Seiten 151-161, Elsevier Ltd., 2003. 28 111

LITERATURVERZEICHNIS [Potthast et al. 2007a] Christian Potthast, Jens Twiefel and J¨org Wallaschek, 2007: Modelling Approaches for an Ultrasonic Percussion Drill. In Journal of Sound and Vibration, Vibro-Impact Systems, Elsevier Science, 2007. 5, 29, 73, 79, 80 [Potthast et al. 2007b] Christian Potthast, Rocco Eisseler, Uwe Heisel, Detlef Klotz, J¨org Wallaschek, 2007: Piezoelectric Actuator Design for Ultrasonically Assisted Deep Hole Drilling. In Journal of Electroceramics, Springer Netherlands, 2007. 29 [Seifried, Schiehlen, Eberhard 2004] R. Seifried, W. Schiehlen, P. Eberhard, 2004: Numerical and experimental evaluation of the coefficient of restitution for repeated impacts. In International Journal of Impact Engineering, Vol. 32, Seiten 508-­524, Elsevier Ltd., 2005. 63 [Strogatz 1994] Steven H. Strogatz, 1994: Nonlinear Dynamics and Chaos. Perseus Books Publishing, LLC, Cambridge, Massachusetts. 67 [Thompson& Stewart 1986] Michael Thompson, Bruce Stewart, 1986: Nonlinear Dynamics and Chaos. John Wiley and Sons Ltd., Chichester, New York, Brisbane, Toronto, Singapore. 22, 29 [Tongue 1987] Benson H. Tongue and K. Gu, 1987: A higher order method of interpolated cell mapping. In Journal of Sound and Vibration, Vol. 125(1), Seiten 169-179. 1988. 26, 27 [Tufillaro, Mello, Choi, Albano 1986] N. B. Tufillaro, T. M. Mello, Y. M. Choi and A. M. Albano: Period doubling boundaries of a bouncing ball. In Journal Physique, Vol. 47, Seiten 1477-1482. 29 [Van der Pol 1920] Balthasar van der Pol, 1920: A theory of the amplitude of free and forced triode vibrations. Radio Review, 1, 701-710, 754-762. Ref. aus: Takashi Kanamaru, Jan. 2007: Van der Pol Oscillator. Scholarpedia, p.7075. 21 [White& Tongue 1994] M. T. White and Benson H. Tongue, 1994: Application of interpolated cell mapping to an analysis of the Lorenz equations. In Journal of Sound and Vibration, Vol. 188(2), Seiten 209-226. 1995. 26 [Wiercigroch, Wojewoda, Krivtsov 2000] M. Wiercigroch, J. Wojewoda, A. M. Krivtsov, 2000: Dynamics of ultrasonic percussive drilling of hard rocks. In Journal of Sound and Vibration, Vol. 280, Seiten 739-757, Elsevier Ltd., 2004. 29 [Wikipedia 2007] Freie Enzyklop¨adie Wikipedia, 2007: Joseph-Louis Lagrange. Internetreferenz: http://de.wikipedia.org/wiki/JosephLouis Lagrange. Wikimedia Foundation Inc., Stand 30. Juli 2007. 22 112