Steffen Roth Baustatik und Baudynamik Algorithmen zur nichtlinearen Stabilitätsanalyse dünnwandiger Strukturen Algorithmen zur nichtlinearen Stabilitätsanalyse dünnwandiger Strukturen von Steffen Roth Bericht Nr. 71 Institut für Baustatik und Baudynamik der Universität Stuttgart Professor Dr.-Ing. habil. M. Bischoff 2020 © Steffen Roth Berichte können bezogen werden über: Institut für Baustatik und Baudynamik Universität Stuttgart Pfaffenwaldring 7 70550 Stuttgart Tel.: 0711 - 685 66123 Fax: 0711 - 685 66130 E-Mail: sekretariat@ibb.uni-stuttgart.de http://www.ibb.uni-stuttgart.de/ Alle Rechte, insbesondere das der Übersetzung in andere Sprachen, vorbehalten. Ohne Geneh- migung des Autors ist es nicht gestattet, diesen Bericht ganz oder teilweise auf photomechani- schem, elektronischem oder sonstigem Wege zu kommerziellen Zwecken zu vervielfältigen. D 93 - Dissertation an der Universität Stuttgart ISBN 978-3-00-067327-6 Algorithmen zur nichtlinearen Stabilitätsanalyse dünnwandiger Strukturen Von der Fakultät Bau- und Umweltingenieurwissenschaften der Universität Stuttgart zur Erlangung der Würde eines Doktor-Ingenieurs (Dr.-Ing.) genehmigte Abhandlung vorgelegt von Steffen Roth aus Bad Kissingen Hauptberichter: Prof. Dr.-Ing. habil. Manfred Bischoff, Stuttgart Mitberichter: Prof. Dr.-Ing. habil. Werner Wagner, Karlsruhe Tag der mündlichen Prüfung: 22. Juli 2020 Institut für Baustatik und Baudynamik der Universität Stuttgart 2020 Kurzfassung Kurzfassung Die wissenschaftliche Auseinandersetzung mit der Stabilität von Tragwerken reicht in der Geschichte weit zurück. Leonhard Euler befasste sich bereits im Jahr 1744 mit dem Stabilitätsversagen druckbelasteter Stäbe. Bis heute ist die Frage nach der Stabilität in vielen Bereichen, vom Bauwesen bis zur Luft- und Raumfahrt, relevant für die Bemes- sung und Konstruktion einzelner Bauteile sowie gesamter Tragstrukturen. Die Analyse eines Tragwerks hinsichtlich seiner Stabilität kann dabei in drei Abschnitte untergliedert werden: erstens die Untersuchung des Systems vor dem Stabilitätsversagen, zweitens die Eigenschaften des Systems an einem kritischen Punkt sowie drittens die Untersuchung des Verhaltens im Nachbeulbereich. Für eine geometrisch nichtlineare Analyse des Vor- und Nachbeulbereichs existiert eine Vielzahl inkrementell-iterativer Lösungsverfahren, die seit den 1970er Jahren entwickelt wurden. Sie ermöglichen die zuverlässige Berechnung nichtlinearer Gleichgewichtspfade einer Struktur, auch wenn diese in Last oder Verschiebung rückläufig sind. In den ver- gangen Jahrzehnten wurden diese Methoden stetig weiterentwickelt und verbessert, etwa durch die Erweiterung auf mehrere Parameter. So lassen sich für Systeme mit mehreren Veränderlichen die resultierenden Gleichgewichtspfade zu Gleichgewichtsflächen erwei- tern. Mithilfe begleitender Maßnahmen kann zudem an jedem Gleichgewichtspunkt eine Aussage über die Stabilität des Systems getroffen werden. Der Fokus der vorliegenden Arbeit liegt auf dem zweiten der oben genannten Punkte. Es werden Methoden präsentiert, die eine direkte Berechnung von Durchschlags- und Verzweigungspunkten für Systeme ermöglichen, die durch eine von außen aufgebrachte Knotenverschiebung belastet werden. Hierfür wird die für Kraftlastfälle bereits verfüg- bare Methode der erweiterten Systeme, die in der Literatur meistens mit dem englischen Begriff Extended Systems bezeichnet wird, modifiziert und ergänzt. Für die direkte Be- rechnung von Durchschlagspunkten wird die Residuumsgleichung für das Gleichgewicht um eine Stabilitätsbedingung erweitert. Diese ist ausschließlich an kritischen Punkten erfüllt und lässt sich aus dem Indifferenzkriterium nach Euler herleiten. Eine Neben- bedingung, die eine triviale Lösung vermeidet und die Anzahl der Gleichungen und die Anzahl der Unbekannten ausgleicht, ergibt sich aus der Normierung der Länge des Ei- genvektors. Für die direkte Ermittlung von Verzweigungspunkten muss zusätzlich das Kriterium zur Unterscheidung von Durchschlags- und Verzweigungspunkten für den Fall inhomogener Dirichlet-Randbedingungen neu hergeleitet werden. Nach einer konsisten- ten Linearisierung des nichtlinearen Gleichungssystems lässt sich die Lösung mittels inkrementell-iterativer Methoden ermitteln. Mit den präsentierten Methoden wird das Anfangsnachbeulverhalten für Systeme mit inhomogenen Dirichlet-Randbedingungen numerisch analysiert, wodurch die bisherige Kategorisierung nach Koiter (1967) er- weitert wird. Zur Linearisierung des nichtlinearen Gleichungssystems ist die Ermittlung einer Rich- tungsableitung der Steifigkeitsmatrix in Richtung des Beulvektors erforderlich. In der i Kurzfassung bisherigen Literatur wird die Ableitung mittels Differenzenquotienten angenähert. Die- se Methoden benötigen eine geringe Rechenzeit, ihre Eigenschaften hängen jedoch von einem Parameter, der Größe des Differenzenschritts, ab. Wird dieser zu groß oder zu klein gewählt, ist der resultierende Fehler zu groß und die iterative Berechnung konver- giert nur langsam oder gar nicht zur gesuchten Lösung. In dieser Arbeit werden erstmals Methoden zur Ermittlung totaler und partieller Ableitungen, die auf hyperkomplexen Zahlen basieren, mit der Methode der Extended Systems kombiniert. Diese Methoden sind unabhängig von der Parameterwahl und ermöglichen eine exakte Berechnung erster und zweiter Ableitungen. In der vorliegenden Arbeit werden die beschriebenen Metho- den auf die Ermittlung der Richtungsableitung übertragen. Bei der Untersuchung der Tragfähigkeit einer stabilitätsgefährdeten Struktur haben Im- perfektionen einen großen Einfluss auf die Traglast. Bereits kleinste Abweichungen in der Geometrie eines Tragwerks oder Eigenspannungen im Material können, besonders bei schlanken Tragwerken, zu einer deutlichen Abminderung der Traglast führen. Um den Einfluss geometrischer Imperfektionen genauer untersuchen zu können, werden in dieser Arbeit zwei Methoden präsentiert, die sich in ihrem Umgang mit Imperfektionen grund- legend voneinander unterscheiden. Bei der ersten Methode ist die Imperfektionsform eine gegebene Größe. Durch eine effiziente, sich wiederholende Berechnung kritischer Punkte für eine Vielzahl an Imperfektionsformen oder -amplituden lassen sich damit kritische Pfade ermitteln. Dies geschieht in der vorliegenden Arbeit durch eine geeignete Wahl der Prädiktoren für die unterschiedlichen Imperfektionsformen bzw. über eine zusätzli- che Erweiterung der Gleichungssysteme zur direkten Berechnung kritischer Punkte um eine Zusatzgleichung zur Pfadverfolgung. Bei der zweiten Methode werden die Imper- fektionen als Unbekannte in das Gleichungssystem eingebracht. Durch eine zusätzliche Bedingung, die sich aus der Variation des Potentials nach den Imperfektionen ergibt, sowie das Einbinden der Nebenbedingung über eine Straffunktion lassen sich mit dieser Methode ungünstigste Imperfektionsformen finden. Diese Methode ist eine Weiterent- wicklung der von Deml und Wunderlich (1997) vorgschlagenen Vorgehensweise. Abschließend werden die hier vorgestellten Methoden und Algorithmen an einer Aus- wahl an numerischen Experimenten demonstriert. Anhand der Beispiele wird u. a. ge- zeigt, dass die Richtungsableitung für eine Vielzahl an Elementformulierungen effizient ermittelt werden kann. Hierbei kommen vom einfachen ebenen Stabelement über isogeo- metrische Schalenelemente bis hin zu Volumenschalenelementen, die durch eine Erwei- terung des Verzerrungsansatzes künstliche Versteifungseffekte reduzieren, zum Einsatz. Zudem wird die Anwendbarkeit der in der vorliegenden Arbeit entwickelten Methoden auf nichttriviale Systeme im Rahmen der Finite-Elemente-Methode mit vielen Freiheits- graden nachgewiesen. ii Abstract Abstract The scientific debate on the stability of structures goes back a long way in history. Leonhard Euler dealt with the stability failure of pressure-loaded bars as early as 1744. To this day, the question of stability is relevant for the design and construction of individual components and entire supporting structures in many areas of engineering, from civil engineering to aerospace. The analysis of a structure regarding its stability can be divided into three sections: first, the investigation of the system prior to stability failure, second, the system properties directly at a critical point and third, the investigation of the postbuckling behavior. For a geometrically non-linear analysis of the pre- and postbuckling behavior, a variety of incremental-iterative solution methods exist, which have been developed since the 1970s. They allow the reliable calculation of non-linear equilibrium paths with snap- backs ans snap-throughs in load and displacement. In the past decades, these methods have been continuously developed and improved, for example by extending them to multiple parameters. Thus, the resulting equilibrium areas can be determined for systems with several variables. Furthermore, conclusions can be drawn about the stability of the system at each equilibrium point with the help of accompanying measures. The focus of the present thesis is laid on the second of the above-mentioned issues. Methods are presented that allow a direct calculation of limit and bifurcation points for systems loaded by an externally applied nodal displacement. For this purpose, the method of extended systems, which already exists for force load cases, is modified and supplemented. For the direct calculation of limit points, the residual equation for equi- librium is extended by a stability condition. This condition is only fulfilled at critical points and can be derived from the indifference criterion according to Euler. A constraint that avoids a trivial solution and balances the number of equations and the number of unknowns is obtained by normalizing the length of the eigenvector. For the direct de- termination of bifurcation points, the criterion for distinguishing limit and bifurcation points for the case of non-homogeneous Dirichlet boundary conditions must also be de- rived accordingly. After consistent linearization of the non-linear system of equations, the solution can be determined by incremental-iterative procedures. With the presented methods the initial post buckling behaviour according to Koiter (1967) is analyzed nu- merically for systems including inhomogeneous Dirichlet boundary conditions, whereby the previously described categorization can be extended. Due to the linearization of the non-linear system of equations, it is necessary to deter- mine a directional derivative of the stiffness matrix into the direction of the buckling vector. In the existing literature, the derivation is approximated with difference quoti- ents. These methods require only a short computing time, but depend on a parameter to be selected. If this parameter is chosen too large or too small, the resulting error is too large, respectively and the iterative calculation requires more iterations or diver- ges within a finite number of iterations. In this thesis, methods for determining total iii Abstract and partial derivatives are combined with the method of extended systems based on hyper-complex numbers for the first time. These methods are independent of parameter choice and allow an exact calculation of first and second derivatives. Moreover, in the present thesis, the described methods are applied to the determination of directional derivatives. When investigating the load-bearing capacity of a structure at risk of stability, imper- fections have a major influence on the critical load. Even the smallest deviations in the geometry of a structure or residual stresses in the material can lead to a significant re- duction in the load-bearing capacity, especially in the case of slender structures. In order to investigate the influence of geometric imperfections more precisely, two methods are presented in this thesis, which differ fundamentally in their handling of imperfections. In the first method, the shape of the imperfection is prescribed. By an efficient, repetitive calculation of critical points for a large number of imperfection shapes or amplitudes, critical paths, also called fold lines, can be determined. This is done in the present work by a suitable predictor choice for the different imperfection modes or by an additional extension of the systems of equations for the direct calculation of critical points by an additional equation for path tracing. In the second metho, the imperfections are in- troduced as unknowns into the system of equations. By using an additional condition resulting from the variation of the potential energy with respect to the imperfections as well as by including the constraint via a barrier function, this method can be used to find unfavorable imperfection shapes. Finally, a selection of the methods and algorithms presented in this thesis are demonstra- ted in numerical experiments. Among other things, the examples show that the directio- nal derivative can be efficiently determined for a large number of element formulations, from simple truss elements to isogeometric mid-surface shell elements to solid shell ele- ments, which reduce locking effects by an extension of the strain ansatz. In addition, the applicability of the methods developed in this thesis to non-trivial systems within the framework of FEM with many degrees of freedom is demonstrated. iv Danksagung Danksagung Die vorliegende Arbeit entstand während meiner Zeit als wissenschaftlicher Mitarbeiter am Institut für Baustatik und Baudynamik der Universität Stuttgart. Ein ganz besonderer Dank gebührt meinem Doktorvater Prof. Dr.-Ing. habil. Manfred Bischoff. Lieber Manfred, vielen Dank, dass du dich gemeinsam mit mir auf das Ex- periment eingelassen hast. Die vielen fachlichen Diskussionen mit dir haben regelmäßig neue spannende Fragen aufgeworfen und meine Forschungsarbeit immens bereichert. Auch der Einsatz für dein Institut und deine Mitarbeiter ist vorbildhaft und hat mich, besonders in den schwierigen Phasen immer aufs Neue motiviert. Daneben sorgst du mit deinem Humor und deinem weit über die Baustatik hinausreichendem Interesse für eine unvergleichbare Stimmung am Institut. Des Weiteren bedanke ich mich bei Herrn Prof. Dr.-Ing. habil. Werner Wagner für das Interesse an meiner Arbeit, für die Übernahme des Mitberichtes sowie das schnelle Anfertigen des Gutachtens. Bei Herrn Prof. Dr.-Ing. Dr.-Ing. E.h. Dr. h.c. Ekkehard Ramm bedanke ich mich für die fruchtbaren und erkenntnisreichen Diskussionen rund um mein Forschungsthema und darüber hinaus. Ein besonderer Dank geht an alle meine ehemaligen Kolleg*innen am Institut für Bausta- tik und Baudynamik. Ihr habt durch den fachlichen Austausch und die zahlreichen Dis- kussionen, aber auch durch die vielen erfreulichen Abende, zum Erfolg dieser Arbeit mit beigetragen. Besonders hervorheben möchte ich an dieser Stelle Dr.-Ing. Malte von Scheven, der mir und meinem PC rund um die Uhr aus der Patsche geholfen hat sowie Florian Geiger, der über Jahre das Büro mit mir teilte und somit mein erster Ansprech- partner bei allen wissenschaftlichen und privaten Fragen war. Für die gewissenhafte Korrektur und die konstruktive Kritik an meiner schriftlichen Arbeit bedanke ich mich bei allen Kolleg*innen und Freund*innen, die mir dabei geholfen haben. Zu guter Letzt bedanke ich mich bei all meinen Freunden, meiner Familie und mei- ner Freundin Florina von ganzem Herzen für all die Liebe und Unterstützung während meines Studiums und der Promotion. Stuttgart, im Oktober 2020 Steffen Roth „Ich habe fertig.“ – Giovanni Trapattoni– v vi Inhaltsverzeichnis Abbildungsverzeichnis xi Tabellenverzeichnis xv Abkürzungen und Bezeichnungen xvii 1 Einleitung 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Stand der Technik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Zielsetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Gliederung der Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Grundlagen der Kontinuumsmechanik 9 2.1 Nichtlineare Festkörpermechanik . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 Bewegung, Kinematik und Verzerrung . . . . . . . . . . . . . . . 10 2.1.2 Spannungsmaße und Konstitutivgleichungen . . . . . . . . . . . . 11 2.1.3 Gleichgewicht . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.4 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Elastostatisches Randwertproblem . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Lösungsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Diskretisierung im Raum . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Linearisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3 Phänomene der geometrisch nichtlinearen Strukturmechanik 21 3.1 Nichtlineare Effekte in der Strukturmechanik . . . . . . . . . . . . . . . . 22 3.2 Statische Gleichgewichtspfade und Gleichgewichtsflächen . . . . . . . . . 24 3.3 Stabilität statischer Gleichgewichtszustände . . . . . . . . . . . . . . . . 28 3.4 Kritische Punkte und kritische Pfade . . . . . . . . . . . . . . . . . . . . 29 3.4.1 Durchschlagspunke bei Kraftlastfällen . . . . . . . . . . . . . . . . 30 vii Inhaltsverzeichnis 3.4.2 Durchschlagspunkte bei Verschiebungslastfällen . . . . . . . . . . 31 3.4.3 Verzweigungspunke . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4.4 Kritische Pfade bei Mehrparametersystemen . . . . . . . . . . . . 32 3.5 Nachbeulverhalten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.5.1 Kraftlastfälle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.5.2 Verschiebungslastfälle . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.5.3 Mehrparametersysteme . . . . . . . . . . . . . . . . . . . . . . . . 39 3.6 Imperfektionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4 Geometrisch nichtlineare Strukturanalyse 43 4.1 Pfadverfolgungsmethoden . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1.1 Prädiktorschritt . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1.2 Konsistente Linearisierung . . . . . . . . . . . . . . . . . . . . . . 51 4.1.3 Korrektoriteration . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.1.4 Bogenlängenkontrolle für inhomogene Dirichlet-Randbedingungen 53 4.2 Auffinden kritischer Punkte . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3 Direkte Berechnung kritischer Punkte für Kraftlastfälle . . . . . . . . . . 55 4.3.1 Unterscheidung von Durchschlags- und Verzweigungspunkten für Kraftlastfälle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3.2 Durchschlagspunkte . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3.3 Verzweigungspunkte . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.4 Direkte Berechnung kritischer Punkte für Verschiebungslastfälle . . . . . 62 4.4.1 Durchschlagspunkte . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.4.2 Verzweigungspunkte . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.5 Wahl der Startwerte bei Extended Systems . . . . . . . . . . . . . . . . . 67 5 Geometrisch nichtlineare Imperfektionsanalyse 69 5.1 Berechnung kritischer Pfade imperfekter Systeme . . . . . . . . . . . . . 70 5.1.1 Wiederberechnung kritischer Punkte an imperfekten Systemen . . 70 5.1.2 Automatisierte Berechnung kritischer Pfade . . . . . . . . . . . . 72 5.2 Berechnung ungünstiger Imperfektionsformen . . . . . . . . . . . . . . . 74 6 Numerische Ableitungen 79 6.1 Klassische Differenzenverfahren . . . . . . . . . . . . . . . . . . . . . . . 80 6.2 Ableiten mit der Complex-step-Methode . . . . . . . . . . . . . . . . . . 84 6.3 Ableiten mit hyperdualen Zahlen . . . . . . . . . . . . . . . . . . . . . . 87 6.4 Rechnen mit hyperkomplexen Zahlen . . . . . . . . . . . . . . . . . . . . 89 6.5 Richtungsableitung der Steifigkeitsmatrix . . . . . . . . . . . . . . . . . . 91 6.6 Automatisches Differenzieren . . . . . . . . . . . . . . . . . . . . . . . . . 94 viii Inhaltsverzeichnis 7 Numerische Experimente 99 7.1 Berechnung kritischer Punkte . . . . . . . . . . . . . . . . . . . . . . . . 99 7.2 Berechnung kritischer Pfade . . . . . . . . . . . . . . . . . . . . . . . . . 110 8 Zusammenfassung und Ausblick 113 8.1 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 8.2 Ausblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 A Anhang 117 A.1 Analyse eines Dreigelenktragwerks mit variabler Höhe . . . . . . . . . . . 117 A.2 Symmetrische entfestigende Feder . . . . . . . . . . . . . . . . . . . . . . 119 Literaturverzeichnis 123 ix Inhaltsverzeichnis x Abbildungsverzeichnis 2.1 Kinematik eines Festkörpers B in der Referenz- und der Momentankonfi- guration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.1 Übersicht über nichtlineare Effekte in der Strukturmechanik, vgl. Wag- ner (2017) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Modell des flachen, symmetrischen Dreigelenktragwerks (links) und der dazugehörige Gleichgewichtspfad (rechts) . . . . . . . . . . . . . . . . . . 25 3.3 Modell des flachen, asymmetrischen Dreigelenktragwerks (links oben) und der dazugehörige Gleichgewichtspfad für die Freiheitsgrade D1 und D2 (rechts oben & links unten), räumliche Darstellung mit Projektionen (rechts unten) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.4 Modell des hohen, symmetrischen Dreigelenktragwerks (links) und die dazugehörige Gleichgewichtskurve (rechts) mit Primär- und Sekundärpfad 27 3.5 Modell eines Dreigelenktragwerks mit variabler Höhe (links) und die da- zugehörige Gleichgewichtsfläche (rechts) . . . . . . . . . . . . . . . . . . 28 3.6 Stabilität eines Gleichgewichtszustand am Beispiel einer Kugel. Stabiler Gleichgewichtszustand (links), indifferenter Gleichgewichtszustand (mit- tig) und instabiler Gleichgewichtszustand (rechts) . . . . . . . . . . . . . 29 3.7 Durchschlagspunkte für Kraftlastfall am Beispiel eines Dreigelenktrag- werks (qualitativ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.8 Durchschlagspunkte für Verschiebungslastfall am Beispiel eines modifi- zierten Dreigelenktragwerks (qualitativ) . . . . . . . . . . . . . . . . . . . 31 3.9 Verzweigungspunkt am Beispiel eines Dreigelenktragwerks (qualitativ) . . 32 3.10 Modell eines Dreigelenktragwerks mit variabler Höhe (links) und die da- zugehörige Gleichgewichtsfläche (rechts) mit kritischen Pfaden . . . . . . 33 3.11 Klassifikation kritischer Punkte für Kraftlastfälle, angelehnt an Koiter (1967) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 xi Abbildungsverzeichnis 3.12 Stabilität des Nachbeulbereichs bei Durchschlagspunkt in Anlehnung an Koiter (1967) unter Verschiebungslast . . . . . . . . . . . . . . . . . . . 36 3.13 Stabiler symmetrischer Verzweigungspunkt in Anlehnung an Koiter (1967) unter Verschiebungslast . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.14 Asymmetrischer Verzweigungspunkt in Anlehnung an Koiter (1967) un- ter Verschiebungslast . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.15 Instabiler symmetrischer Verzweigungspunkt in Anlehnung an Koiter (1967) unter Verschiebungslast . . . . . . . . . . . . . . . . . . . . . . . . 38 3.16 Modell eines Dreigelenktragwerks mit variabler Höhe (links) und die da- zugehörige Gleichgewichtsfläche mit Stabilitätseigenschaften (rechts) . . . 39 3.17 Auswirkungen geometrischer Imperfektionen für Durchschlagsprobleme am Beispiel eines modifizierten Dreigelenktragwerks . . . . . . . . . . . . 40 3.18 Auswirkungen geometrischer Imperfektionen für Verzweigungsprobleme am Beispiel eines druckbelasteten Pendelstabs mit horizontaler Feder . . 41 4.1 Vergleich von Last- und Verschiebungskontrolle sowie deren Grenzen in Anlehnung an Bischoff u. a. (2016) . . . . . . . . . . . . . . . . . . . . 47 4.2 Vergleich der Bogenlängenkontrollen nach Riks & Wempner (links) und Ramm (rechts) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.3 Flussdiagramm zur Berechnung kritischer Punkte bei inhomogenen Dirichlet- Randbedingungen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.1 Exemplarische Darstellung eines Gleichgewichtspfads bei steigender Im- perfektionsamplitude (links), direkte Berechnung der kritischen Punkte (rechts), in Anlehnung an Reitinger (1994) . . . . . . . . . . . . . . . . 71 5.2 Modifiziertes Dreigelenktragwerk mit variabler Imperfektion . . . . . . . 71 5.3 Gleichgewichtsfläche des modifizierten Dreigelenktragwerks mit Gleichge- wichtspfad des perfekten Systems (pink) und kritischem Pfad (blau) . . . 72 5.4 Imperfektes Dreigelenktragwerk . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 Kritischer Pfad des asymmetrischen Dreigelenktragwerks . . . . . . . . . 78 6.1 Übersicht über die in dieser Arbeit vorgestellten Methoden zur numeri- schen Ableitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.2 Steigung der Sekante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.3 Vergleich des relativen Fehlers von Vorwärtsdifferenzen- und zentralem Differenzenquotient für die erste Ableitung von Gl. (6.9) . . . . . . . . . 83 6.4 Vergleich des relativen Fehlers von Vorwärtsdifferenzen- und zentralem Differenzenquotient sowie CSDA für die erste Ableitung von Gl. (6.9) . . 85 6.5 Vergleich des relativen Fehlers von Vorwärtsdifferenzen- und zentralem Differenzenquotient, CSDA und CSDA mit dualen Zahlen für die erste Ableitung von Gl. (6.9) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 xii Abbildungsverzeichnis 6.6 Vergleich des relativen Fehlers von Vorwärtsdifferenzen- und zentralem Differenzenquotient sowie CSDA für die zweite Ableitung von Gl. (6.9) . 87 6.7 Vergleich des relativen Fehlers von Vorwärtsdifferenzen- und zentralem Differenzenquotient, CSDA und hyperdualen Zahlen für die zweite Ablei- tung von Gl. (6.9) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.8 Ebene Scheibe unter Verschiebungslast, Modelldaten und Randbedingun- gen (links), Diagramm mit Vergleich der relativen Rechenzeiten (rechts) . 93 7.1 Statisches System des L-förmigen Rahmentragwerks in Anlehnung an Lee u. a. (1968) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 7.2 Details zur Modellierung des Rahmentragwerks mit vier linearen Q1- Scheibenelementen über die Dicke t . . . . . . . . . . . . . . . . . . . . . 100 7.3 Details zur Modellierung des Rahmentragwerks mit zwei Volumenscha- lenelement über die Dicke t und einem über die Tiefe w . . . . . . . . . . 101 7.4 Gleichgewichtspfade für x-Verschiebung mit Durchschlagspunkten bei un- terschiedlicher Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . 102 7.5 Gleichgewichtspfade für y-Verschiebung mit Durchschlagspunkten bei un- terschiedlicher Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . 102 7.6 Detail der Gleichgewichtspfade für x-Verschiebung mit Durchschlagspunk- ten bei unterschiedlicher Diskretisierung . . . . . . . . . . . . . . . . . . 103 7.7 System der ebenen Scheibe mit Randbedingungen für eine Diskretisierung mit 200 Elementen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.8 Primärpfad (PP) und Sekundärpfade (SP+ und SP–) für die Mittelpunkt- verschiebung des ebenen Scheibentragwerks . . . . . . . . . . . . . . . . . 105 7.9 Konturplots der Verschiebungen und der Eigenvektoren für verschiedene Netzfeinheiten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 7.10 System der gekrümmten Schale mit Randbedingungen . . . . . . . . . . . 107 7.11 Detail des Gleichgewichtspfads der gekrümmten Schale . . . . . . . . . . 108 7.12 Verschiebungen im Nachbeulbereich (links) und Eigenvektor am kriti- schen Punkt (rechts) der gekrümmten Schale bei grober Diskretisierung, überhöhte Darstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.13 Verschiebungen im Nachbeulbereich (links) und Eigenvektor am kriti- schen Punkt (rechts) der gekrümmten Schale bei feiner Diskretisierung, überhöhte Darstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.14 Grundriss und Ansicht der Gitterschale nach Forde (1986) . . . . . . . 110 7.15 Räumliche Ansicht des unverformten Systems (links) und des Eigenvek- tors (rechts) in 50-facher Überhöhung . . . . . . . . . . . . . . . . . . . . 110 7.16 Gleichgewichtsfläche der Gitterschale (grau) mit dem Gleichgewichtspfad des perfekten Systems (pink) und dem kritischem Pfad (blau) . . . . . . 111 xiii Abbildungsverzeichnis A.1 Dreigelenktragwerk mit variabler Höhe . . . . . . . . . . . . . . . . . . . 117 A.2 Gleichgewichtspfade für Dregelenktragwerke unterschiedlicher Höhe mit Durchschlagspunkten (LP) und Verzweigungspunkten (BP) . . . . . . . . 119 A.3 Ungedehnte (links) und gedehnte Feder (rechts) . . . . . . . . . . . . . . 119 A.4 Federkraft (links) und Federsteifigkeit (rechts) . . . . . . . . . . . . . . . 120 xiv Tabellenverzeichnis 7.1 Übersicht über die unterschiedlichen Diskretisierungen . . . . . . . . . . 101 7.2 Vergleich der Rechenzeiten für verschiedene Ableitungsmethoden . . . . . 103 7.3 Konvergenzverhalten für unterschiedlich große Werte des Intervalls h bei Ermittlung der Richtungsableitung mit CSDA und VD . . . . . . . . . . 104 7.4 Ergebnisse für λkrit in Abhängigkeit der Netzfeinheit für ebene Scheibe . 106 7.5 Ergebnisse für λkrit in Abhängigkeit der Netzfeinheit für gekrümmte Schale108 xv Tabellenverzeichnis xvi Abkürzungen und Bezeichnungen Abkürzungen AD . . . . . . . . . . automatisches Differenzieren B . . . . . . . . . . . . Verzweigungspunkt (bifurcation point) CAD . . . . . . . . . computer-aided design CSDA . . . . . . . Complex-step Derivative Approximation CSDAR . . . . . . Richtungsableitung mit CSDA ES . . . . . . . . . . . Extended Systems FEM . . . . . . . . . Finite-Elemente-Methode L . . . . . . . . . . . . Durchschlagspunkt (limit point) UMFPACK . . Unsymmetric MultiFrontal PACKage Mathematische Notationen (•)−1 . . . . . . . . . Inverse von (•) (•)⊺ . . . . . . . . . . Transponierte von (•) |(•)| . . . . . . . . . . Betrag einer skalaren Größe (•) ‖(•)‖2 . . . . . . . Euklid’sche Norm eines Vektors (•) (•)−⊺ . . . . . . . . . Transponierte der Inversen von (•) (• · •) . . . . . . . . Skalarprodukt zwischen Vektoren bzw. zwischen Tensoren (inneres Produkt) (• ⊗ •) . . . . . . . Dyadisches Produkt δ(•) . . . . . . . . . . virtuelle Größe det(•) . . . . . . . . Determinante von (•) div(•) . . . . . . . . Divergenz von (•) bezogen auf die Momentankonfiguration DIV(•) . . . . . . Divergenz von (•) bezogen auf die Referenzkonfiguration xvii Abkürzungen und Bezeichnungen lin(•) . . . . . . . . Linearisierung von (•) d(•) d(•) . . . . . . . . . . totale Ableitung ∂(•) ∂(•) . . . . . . . . . . partielle Ableitung tr(•) . . . . . . . . . Spur von (•) ⋃ . . . . . . . . . . . . Assemblierungsoperator f ′ v(x) . . . . . . . . . Richtungsableitung der Funktion f (x) in Richtung des Vektors v Lateinische Buchstaben b . . . . . . . . . . . . Volumenlast in der Momentankonfiguration B0 . . . . . . . . . . . Volumenlast in der Referenzkonfiguration B . . . . . . . . . . . . B-Operator B . . . . . . . . . . . . Materieller Körper c . . . . . . . . . . . . . allgemeine zu inkrementierende Größe der Zusatzgleichung f ĉ . . . . . . . . . . . . . allgemeine vorgegebene Schrittweite in Zusatzgleichung f C . . . . . . . . . . . . Materialmatrix CCG . . . . . . . . . rechter Cauchy-Green-Tensor C . . . . . . . . . . . . vierstufiger Materialtensor D(•) . . . . . . . . . inkrementelle Änderung der Größe (•) Di . . . . . . . . . . . einzelne Einträge aus D dX, dx . . . . . . differenzielles Linienelement in der Referenz- bzw. Momentankonfi- guration d . . . . . . . . . . . . Vektor der diskreten Knotenverschiebungen auf Elementebene D . . . . . . . . . . . . globaler Vektor der diskreten Knotenverschiebungen D̂ . . . . . . . . . . . . globaler Vektor der vorgeschriebenen Knotenverschiebungen D̃ . . . . . . . . . . . . globaler Vektor der diskreten Knotenverschiebungen in einem Nach- barzustand e . . . . . . . . . . . . . Element der Finite-Elemente-Diskretisierung e1, e2, e3 . . . . Normierte Basisvektoren des kartesischen Koordinatensystems E . . . . . . . . . . . . Elastizitätsmodul E . . . . . . . . . . . . Green-Lagrange’scher Verzerrungstensor f . . . . . . . . . . . . . Zusatzgleichung F0 . . . . . . . . . . . Einzellast fint, fext . . . . . . Vektoren der inneren und äußeren Kräfte auf Elementebene F . . . . . . . . . . . . Deformationsgradient xviii Abkürzungen und Bezeichnungen Fint, Fext . . . . globale Vektoren der inneren und äußeren Kräfte FDir . . . . . . . . . . Ableitung der internen Kräfte nach der vorgegebenen Verschiebungs- änderung G . . . . . . . . . . . . Schubmodul h . . . . . . . . . . . . . Größe des Intervalls bei numerischen Ableitungen i . . . . . . . . . . . . . Iterationsschritt, in Kapitel 6 der Imaginärteil einer komplexen Zahl I . . . . . . . . . . . . . zweistufiger Identitätstensor bzw. Einheitsmatrix Isym . . . . . . . . . . symmetrischer vierstufiger Identitätstensor ℑ(•) . . . . . . . . . Imaginärteil von (•) ℑǫ1 (•) . . . . . . . erster Imaginärteil der hyperdualen Zahl (•) ℑǫ2 (•) . . . . . . . zweiter Imaginärteil der hyperdualen Zahl (•) ℑǫ1ǫ2 (•) . . . . . . dritter Imaginärteil der hyperdualen Zahl (•) J . . . . . . . . . . . . Jacobi-Determinante kT . . . . . . . . . . . Tangentensteifigkeitsmatrix auf Elementebene Ke . . . . . . . . . . . globale elastische Steifigkeitsmatrix Kg . . . . . . . . . . . globale geometrische Steifigkeitsmatrix KT . . . . . . . . . . globale Tangentensteifigkeitsmatrix Ku . . . . . . . . . . . globale Anfangsverschiebungssteifigkeitsmatrix KDir . . . . . . . . . Ableitung der globalen Tangentensteifigkeitsmatrix KT nach der vor- gegebenen Verschiebungsänderung KΦ . . . . . . . . . . Richtungsableitung der globalen Tangentensteifigkeitsmatrix KT in Richtung von Φ Ktot . . . . . . . . . Gesamtsteifigkeitsmatrix L . . . . . . . . . . . . untere Dreiecksmatrix m . . . . . . . . . . . . inkrementeller Schritt n . . . . . . . . . . . . Anzahl der Freiheitsgrade des Systems nele . . . . . . . . . . Anzahl der finiten Elemente N, n . . . . . . . . . Normale in der Referenz- bzw. Momentankonfiguration Ne . . . . . . . . . . . Matrix der Formfunktionen des Elements e Ni . . . . . . . . . . . Formfunktion P . . . . . . . . . . . . Materieller Punkt P . . . . . . . . . . . . erster Piola-Kirchhoff-Spannungstensor R . . . . . . . . . . . . Residuum der Ungleichgewichtskräfte Rtot . . . . . . . . . . Residuum des erweiterten Gleichungssystems R 3 . . . . . . . . . . . euklidischer Raum xix Abkürzungen und Bezeichnungen R(•) . . . . . . . . . Imaginärteil von (•) s . . . . . . . . . . . . . Bogenlänge, in Abschnitt 4.3.1 und 4.4.2 auch Bahnparameter ŝ . . . . . . . . . . . . . vorgeschriebene Bogenlänge S . . . . . . . . . . . . zweiter Piola-Kirchhoff-Spannungstensor SP . . . . . . . . . . . Steifigkeitsparameter t . . . . . . . . . . . . . Zeit bzw. Pseudozeit t . . . . . . . . . . . . . Cauchy-Spannungsvektor t0 . . . . . . . . . . . . Anfangszeitpunkt des betrachteten Zeitintervalls T̂, t̂ . . . . . . . . . . vorgeschriebene Spannung in der Referenz- bzw. Momentankonfigu- ration u . . . . . . . . . . . . Verschiebungsvektor eines materiellen Punkts û . . . . . . . . . . . . vorgeschriebene Verschiebung U . . . . . . . . . . . . obere Dreiecksmatrix Ūii . . . . . . . . . . . Diagonaleinträge der oberen Dreiecksmatrix U v . . . . . . . . . . . . Ableitungsrichtung der Richtungsableitung W . . . . . . . . . . . Verzerrungsenergiedichte x . . . . . . . . . . . . Ortsvektor eines materiellen Punkts in der Momentankonfiguration X . . . . . . . . . . . . Ortsvektor eines materiellen Punkts in der Referenzkonfiguration Griechische Buchstaben α . . . . . . . . . . . . Variable zur Regularisierung des erweiterten Gleichungssystems bei der Berechnung von Verzweigungspunkten β . . . . . . . . . . . . Skalierungsfaktor Γ, γ . . . . . . . . . . Rand des Gebiets Ω0, Ωt in der Referenz- bzw. Momentankonfigura- tion ΓD, γD . . . . . . . Oberfläche mit Dirichlet-Randbedingung in der Referenz- bzw. Mo- mentankonfiguration ΓN, γN . . . . . . . Oberfläche mit Neumann-Randbedingung in der Referenz- bzw. Mo- mentankonfiguration ∆(•) . . . . . . . . . Änderung der Größe (•) zwischen zwei Iterationsschritten δΠ . . . . . . . . . . . virtuelle Arbeit ∆t . . . . . . . . . . . Tangentenvektor in Zusatzgleichung nach Riks/Wempner ∆u . . . . . . . . . . iterative Änderung des Normalenvektors in Zusatzgleichung nach Riks/Wempner ε, ε+, ε− . . . . . geometrische Imperfektionen xx Abkürzungen und Bezeichnungen ǫ . . . . . . . . . . . . . Imaginäre Einheit einer dualen Zahl ǫ1, ǫ2, ǫ1ǫ2 . . . Imaginäre Einheiten einer hyperdualen Zahl λ, λi . . . . . . . . . Steigerungsfaktor λL . . . . . . . . . . . 1. Lamé-Konstante µL . . . . . . . . . . . 2. Lamé-Konstante ν . . . . . . . . . . . . . Querdehn- bzw. Poissonzahl Π . . . . . . . . . . . . potenzielle Energie ρ0, ρ . . . . . . . . . Dichte in der Referenz- bzw. Momentankonfiguration σ . . . . . . . . . . . . Cauchy-Spannungstensor Φ . . . . . . . . . . . . Eigenvektor χ . . . . . . . . . . . . Abbildungsfunktion zwischen Referenz- und Momentankonfiguration Ψ . . . . . . . . . . . . Korrekturfaktor ω . . . . . . . . . . . . Eigenwert Ω0, Ωt . . . . . . . Gebiet des materiellen Körpers B in der Referenz- bzw. Momentan- konfiguration Ωe . . . . . . . . . . . Gebiet eines finiten Elements Kopf- und Fußzeiger (•)i . . . . . . . . . . Größe (•) im Iterationsschritt i (•)e . . . . . . . . . . Größe (•) gehört zum finiten Element e (•)krit . . . . . . . . Größe (•) an einem kritischen Punkt (•)tot . . . . . . . . . Größe (•) mit mehr als n Einträgen xxi 1 Einleitung „Es bewegt sich alles, Stillstand gibt es nicht“ schrieb der Künstler Jean Tinguely in seinem Manifest „Für Statik“ im Jahr 1959. In einer Welt ohne Stillstand kann es somit auch keine Stabilität geben, leitet sich der Begriff der Stabilität doch vom Lateinischen stabilis ab, dass mit fest stehend übersetzt werden kann. Anders als in dieser radikalen Verneinung der Stabilität im Sinne statischer Gleichgewichtszustände, die aus dem Zitat Tinguelys abgeleitet werden kann, ist die Existenz stabiler Zustände, aber auch deren Grenzen, in vielen wissenschaftlichen Disziplinen von großer Bedeutung. In der Kernphy- sik etwa, in der die Stabilität eines Atomkerns gegenüber radioaktivem Zerfall untersucht wird. Auch in den Wirtschaftswissenschaften oder der Thermodynamik werden Systeme auf ihre Stabilität und deren Grenzen hin analysiert. Dabei stellt sich immer die Frage nach der maximalen Größe einer auf das System einwirkenden Störung, die gerade noch zu keiner signifikanten Abweichung in der Systemantwort führt (vgl. Schmitz (2016)). In der vorliegenden Arbeit geht es um die Frage nach der Stabilität in der geometrisch nichtlinearen Strukturmechanik. Neben dem Verhalten der Struktur im stabilen und instabilen Zustand ist vorrangig die Frage nach der Berechnung stabilitätskritischer Zu- stände von Interesse. Im folgenden Abschnitt wird zunächst kurz die Motivation erläutert, die zu dieser Arbeit geführt hat. Es folgen ein kurzer Überblick über den Stand der Technik, die Ziele der Arbeit sowie deren Gliederung. 1.1 Motivation Die Frage nach der Stabilität einer dünnwandigen, druckbelasteten Struktur ist in vie- len Bereichen der Ingenieurwissenschaften von großem Interesse. Dabei wird der Begriff 1 1 Einleitung häufig im Zusammenhang mit dem Versagen eines Bauteils oder eines gesamten Trag- werks in Verbindung gebracht. Die Versagensszenarien des Knickens, Biegedrillknickens oder Beulens sowie die rechnerischen und konstruktiven Vorschriften, um eben dies zu vermeiden, sind für die jeweiligen Strukturen und Werkstoffe in den entsprechenden Normen geregelt. Daneben gibt es seit einigen Jahren immer mehr Forschungsarbeiten mit dem Ziel, Sta- bilitätsversagen gezielt zu nutzen. Beispiele hierfür sind z. B. die Untersuchung und Entwicklung auxetischer Metamaterialien (Lazarus und Reis (2015)) oder die ge- zielte Veränderung der Topografie und damit der aerodynamischen Eigenschaften von Tragflächen (Terwagne u. a. (2014)). In der Arbeit von Reis (2015) werden diese verschiedenen Sichtweisen mit Buckliphobia und Buckliphilia bezeichnet. Für beide Fälle ist zunächst das nichtlineare Verhalten der Struktur, also der nichtlinea- re Zusammenhang zwischen der Kraft, die auf ein Tragwerk wirkt und der deformierten Geometrie von Interesse. Da aufgrund der Nichtlinearität des Problems die Lösung nicht direkt ermittelt werden kann, kommen hierfür inkrementell-iterative Methoden zum Ein- satz. Den Übergang zwischen stabilem und instabilem Gleichgewicht bilden sogenannte kritische Punkte, an denen die Steifigkeitsmatrix singulär wird. Man unterscheidet dabei zwischen Durchschlags- und Verzweigungspunkten, wobei Verzweigungspunkte nochmals untergliedert werden können in symmetrisch stabile, symmetrisch instabile sowie asym- metrische Verzweigungspunkte. Sowohl für die Vermeidung von Schäden an Tragwerken als auch für die gezielte Manipulation gewünschter Struktureigenschaften ist eine exakte Berechnung dieser kritischen Punkte hilfreich. Imperfektionen können die Traglast dünnwandiger Strukturen stark reduzieren. Hierzu existiert bereits eine Vielzahl an Forschungsarbeiten. Die Suche nach der ungünstigsten Imperfektionsform gleicht dabei häufig der Suche nach der Nadel im Heuhaufen. Ob für reale Tragwerke die eine ungünstigste Imperfektionsform überhaupt existiert, wird in Schneider u. a. (2005) zudem hinterfragt. Aus akademischer Sicht ist dennoch von Interesse, Imperfektionsformen systematisch zu untersuchen und die zumindest theore- tisch ungünstigste Imperfektionsform zu finden. 1.2 Stand der Technik Bei der Erstellung des folgenden Abschnitts wurde auf die Arbeiten von Wagner (1991) und Reitinger (1994) zurückgegriffen. Die ersten Stabilitätsuntersuchungen gehen in der Zeit weit zurück. Euler untersuchte bereits 1744 das Ausknicken druckbelasteter Stäbe. Ausführliche Untersuchungen des 2 1.2 Stand der Technik Stabilitätsverhaltens von Schalen erfolgen vor etwa einem Jahrhundert in den Arbeiten von Zoelly (1915) (Kugelschalen) sowie Lorenz (1908), Timoschenko (1910) und Flügge (1932) (Kreiszylinderschalen). Die hierbei rechnerisch ermittelten Grenzlasten wichen jedoch weit von den experimentell ermittelten Lasten ab. In den Arbeiten von Von Kármán und Tsien (1939) und Von Kármán und Tsien (1941) wurde aus diesem Grund die Untersuchung auf den Nachbeulbereich erweitert. Eine asymptotisch exakte Näherung des Nachbeulbereichs für kontinuierliche elasti- sche Strukturen, auch Anfangsnachbeulanalyse genannt, lieferte die Arbeit von Koiter (1967), die bereits im Jahr 1945 in niederländischer Sprache erschien. Durch eine Reihen- entwicklung des Belastungs- und Verschiebungszustands im unmittelbaren Nachbeulbe- reich lässt sich der Sekundärpfad in der Umgebung des Verzweigungspunkts ermitteln. Das Verhalten des Sekundärpfads im Nachbeulbereich ermöglicht schließlich eine Ka- tegorisierung der verschiedenen Verzweigungspunkte. Etwa zwanzig Jahre nach Koiter erfolgte eine Weiterentwicklung der Anfangsnachbeulanalyse u. a. durch Arbeiten von Sewell (1965), Thompson (1969), Budiansky und Hutchinson (1972), Thomp- son und Hunt (1973b) sowie Budiansky (1974). Fast alle bisher aufgeführten Arbeiten beschäftigen sich mit der Untersuchung des nicht- linearen Systemverhaltens und der Stabilität kontinuierlicher Systeme. Mit dem Auf- kommen des Computers in den 60er Jahren des vergangenen Jahrhunderts verlagern sich die Untersuchungen hin zur numerischen Analyse diskreter Systeme. Die vermut- lich erste Arbeit, die sich aus numerischer Sicht mit dem Verzweigen befasst, stammt von Thurston (1969) und lässt sich als numerisches Pendant zur Koiter’schen An- fangsnachbeulanalyse interpretieren. Mit der immer schneller voranschreitenden Entwicklung der Finite-Elemente-Methode in den 1960er und 1970er Jahren verschob sich der Forschungsschwerpunkt immer weiter hin zu numerischen Stabilitätsuntersuchungen. Dies liegt zum einen an der Entwicklung neuer Elementformulierungen für Stäbe, Scheiben und Schalen, mit denen sich schlanke Tragwerke untersuchen lassen. Es liegt zum anderen an der Entwicklung verschiedener Methoden zur Pfadverfolgung bei nichtlinearem Strukturverhalten. Diese Entwicklung fand sowohl in den Ingenieurwissenschaften wie auch in der numerischen Mathema- tik statt. Wichtige Arbeiten in der Ingenieurliteratur stammen dabei von Wempner (1971), Riks (1972), Riks (1979), Crisfield (1981), Ramm (1981) sowie Schweizer- hof und Wriggers (1986). Im Bereich der Mathematik sind hier die Arbeiten von Keller (1977), Rheinboldt (1981) oder Fried (1984) zu nennen. Eine Erweiterung der Problemstellung auf mehrere Kontrollparameter, in der Literatur auch als generali- sierte Pfadverfolgung bezeichnet, findet sich u. a. in Spence und Jepson (1984) sowie Eriksson (1998). 3 1 Einleitung Zur Berechnung kritischer Punkte wurden zunächst Näherungsverfahren entwickelt. Da- bei werden begleitend zu einer Pfadverfolgung z. B. die Determinante der Steifigkeits- matrix oder deren Eigenwerte analysiert. Wird entlang eines Gleichgewichtspfads ein kritischer Punkt überschritten, kann man sich diesem über Bisektionsmethoden annä- hern. Hier sei u. a. auf die Arbeiten von Brendel (1979), Brendel und Ramm (1982) und Wagner und Wriggers (1988) verwiesen. Die direkte Berechnung kritischer Punkte findet sich zuerst in der mathematischen Li- teratur. Durch eine Erweiterung der Gleichgewichtsbedingung um eine Stabilitätsbedin- gung, die nur an einem kritischen Punkt erfüllt ist, ließen sich erstmals Durchschlags- und Verzweigungspunkte direkt berechnen. Wichtige Arbeiten sind hierbei Keener und Keller (1973), Seydel (1977), Moore (1980) sowie Spence und Jepson (1984). Die genannten Arbeiten beschränken sich dabei jedoch auf Systeme mit nur wenigen Freiheitsgraden. Eine Übertragung dieser Methoden auf Anwendungen im Ingenieurwe- sen, bei denen die Zahl der Freiheitsgrade deutlich größer ist, gelang durch die Arbeit von Wriggers u. a. (1988) und Wriggers und Simo (1990). Auch hier lässt sich, wie bereits bei der Pfadverfolgung, die Problemstellung auf Systeme mit mehr als einem Lastparameter erweitern. Untersuchungen hierzu finden sich beispielsweise in Eriks- son (1994), Eriksson (1997), Rezaiee-Pajand und Moghaddasie (2014), Cox u. a. (2018) sowie Groh u. a. (2018). Wie bereits erwähnt überschätzen rechnerisch ermittelte Grenzlasten die bei realen Ex- perimenten ermittelten Grenzlasten deutlich. Dies lässt sich durch Imperfektionen in der Geometrie oder dem Material erklären, die in einem Rechenmodell mit perfekter Geo- metrie und homogenem Materialverhalten nicht berücksichtigt werden. Im Bauwesen werden daher häufig eigenformaffine Vorverformungen berücksichtigt. Um die ungüns- tigste Imperfektionsform rechnerisch ermitteln zu können, wurden die Imperfektionen in der Arbeit von Deml (1997) und Deml und Wunderlich (1997) als zusätzliche Unbekannte eingeführt. Einen anderen Ansatz, der sich stärker an realen Tragwerken orientiert und auf einem stochastischen Ansatz für Imperfektionsformen beruht, verfol- gen die Arbeiten von Lauterbach u. a. (2018) und Fina u. a. (2020). Wie bereits bei den Pfadverfolgungsmethoden wurden auch die Methoden zur exak- ten Berechnung erster und zweiter Ableitungen auf unterschiedliche Weise sowohl in den Ingenieurwissenschaften als auch in der Mathematik entwickelt. In der numerischen Mathematik werden die hier vorgestellten Methoden unter dem Bergriff automatic diffe- rentiation zusammengefasst. Sie beruhen auf einer Zerlegung der abzuleitenden Funkti- on in grundlegende mathematische Operationen und Funktionen, die anschließend unter Berücksichtigung der Kettenregel abgeleitet werden können. Eine erste Beschreibung der Methode findet sich in der Arbeit von Wengert (1964), die in Rall (1983) und Rall (1986) weiterentwickelt wurde. Eine sehr gute Zusammenfassung findet sich dazu 4 1.3 Zielsetzung in Griewank und Walther (2008). Die in den Ingenieurwissenschaften entwickelten Methoden lassen sich durch die Zuhilfenahme hyperkomplexer Zahlen aus den klas- sischen Differenzenquotienten herleiten. Erste Arbeiten finden sich bereits in Lyness und Moler (1967) und Lyness (1967). Sie bilden die Grundlage für die Algorithmen in Squire und Trapp (1998) und Martins u. a. (2003) zur Berechnung erster Ab- leitungen und in Fike und Alonso (2011) zur exakten Berechnung erster und zweiter Ableitungen. 1.3 Zielsetzung Das Ziel der vorliegenden Arbeit ist die Entwicklung zuverlässiger Algorithmen zur direk- ten Berechnung kritischer Punkte für statische Systeme, die durch inhomogene Dirichlet- Randbedingungen belastet sind. Weiterhin sollen diese Algorithmen zur Untersuchung der Imperfektionsempfindlichkeit schlanker Tragwerke erweitert werden. Zur Berücksichtigung der Verschiebungslastfälle bei der direkten Berechnung von Durch- schlags- und Verzweigungspunkten werden die auch als Extended Systems bezeichneten Methoden, die in Wriggers u. a. (1988) erstmals im Zusammenhang mit der Finite- Elemente-Methode vorgestellt wurden, angepasst. Dafür müssen alle notwendigen Bedin- gungen für Verschiebungslasten neu hergeleitet werden. Die neuen Bedingungen lassen sich anschließend in einem erweiterten nichtlinearen Gleichungssystem zusammenfas- sen. Bevor das nichtlineare Gleichungssystem gelöst werden kann, muss dieses zunächst kon- sistent linearisiert werden. Im anschließenden iterativen Lösungsprozess muss die Rich- tungsableitung der Steifigkeitsmatrix in Richtung des zum Nulleigenwert gehörenden Eigenvektors bestimmt werden. Zur Berechnung der Richtungsableitung werden in die- sem Zusammenhang erstmals numerische Ableitungsalgorithmen angewandt, die auf hy- perkomplexen Zahlenformaten beruhen. Die Algorithmen zur numerischen Berechnung erster und zweiter Ableitungen nach Martins u. a. (2003) und Fike und Alonso (2011) werden hierfür auf die Berechnung der Richtungsableitung übertragen. Ziel ist dabei die zuverlässige Berechnung der Richtungsableitung, die mittels komplexzahliger Ableitungsmethoden praktisch fehlerfrei und parameterunabhängig möglich ist. Aufbauend auf den bisher entwickelten Methoden soll in einem nächsten Schritt die Im- perfektionsempfindlichkeit schlanker Strukturen untersucht werden. Hierfür müssen die Algorithmen erneut erweitert werden. Ziel ist dabei zum einen eine effiziente Berech- nung kritischer Punkte für eine Vielzahl imperfekter Ausgangssysteme. Zum anderen soll es eine Modifikation der in Deml (1997) vorgestellten Methoden zur Bestimmung der ungünstigsten Imperfektionsformen ermöglichen, diese auch dann zu finden, wenn 5 1 Einleitung bei zunehmender Imperfektionsamplitude die kritische Last nach einem anfänglichen Abfall wieder ansteigt. 1.4 Gliederung der Arbeit Kapitel 2 liefert einen Überblick über die in dieser Arbeit notwendigen Grundglei- chungen der nichtlinearen Festkörpermechanik. Sie ermöglichen die Formulierung des elastostatischen Randwertproblems in der starken Form, das anschließend mit der Me- thode der gewichteten Residuen in die schwache Form überführt wird. Diese bildet die Grundlage zur Lösung mit der Finite-Elemente-Methode. In Kapitel 3 wird nach einem allgemeinen Überblick über nichtlineare Effekte in der Strukturmechanik der Fokus auf geometrisch nichtlineare Phänomene gelegt. Neben den Eigenschaften statischer Gleichgewichtspfade und deren Erweiterung zur Gleichgewichts- fläche wird auf den Begriff der Stabilität in dieser Arbeit sowie die Definition kritischer Punkte genauer eingegangen. Des Weiteren wird die Kategorisierung des Anfangsnach- beulverhaltens nach Koiter (1967) auf Verschiebungslastfälle erweitert. Abschließend wird der Einfluss von Imperfektionen auf die Stabilität schlanker Strukturen beschrie- ben. Die Kapitel 4 und 5 befassen sich anschließend mit der Analyse der in Kapitel 3 be- schriebenen Phänomene. Kapitel 4 gibt hierfür zuerst einen Überblick über Methoden zur Verfolgung stati- scher Gleichgewichtspfade, bevor verschiedene Methoden zur Analyse der Stabilität von Gleichgewichtspunkten verglichen werden. Im Anschluss werden erweiterte Gleichungs- systeme zur direkten Berechnung kritischer Punkte für Kraftlastfälle vorgestellt. Darauf aufbauend werden die erweiterten Gleichungssysteme zur direkten Berechnung kritischer Punkte für Verschiebungslastfälle hergeleitet. Kapitel 5 beschäftigt sich mit Möglichkeiten zur systematischen Untersuchung der Imperfektionsempfindlichkeit schlanker Tragwerke. Hierfür wird neben effizienten Algo- rithmen zur wiederholten Berechnung kritischer Punkte für eine Vielzahl imperfekter Systeme auch ein Algorithmus zur Suche nach der ungünstigsten Imperfektionsform vorgestellt, der auf der Arbeit von Deml (1997) beruht. Kapitel 6 befasst sich mit der numerischen Berechnung von partiellen und totalen Ableitungen analytischer Funktionen. Hierfür werden zunächst die klassischen Diffe- renzenquotienten erläutert und die dabei auftretenden Fehler untersucht. Anschließend werden neuere Methoden, die auf den Arbeiten von Martins u. a. (2003) und Fike 6 1.4 Gliederung der Arbeit und Alonso (2011) basieren, vorgestellt. Diese Methoden basieren auf hyperkomple- xen Zahlen und ermöglichen eine fehlerfreie Berechnung von Ableitungen. Abschließend werden die präsentierten Methoden zur Ermittlung der Richtungsableitung der Steifig- keitsmatrix in Richtung des Beulvektors erweitert. In Kapitel 7 wird eine Auswahl der neuen Methoden, die in dieser Arbeit vorgestellt werden, an numerischen Beispielen demonstriert. Hierbei werden die Vorteile und Be- sonderheiten aufgezeigt und mit bisherigen Verfahren oder Beispielen aus der Literatur verglichen. Kapitel 8 fasst alle Neuerungen und Erkenntnisse die sich aus der Arbeit geben, noch einmal zusammen. Abschließend wird an ausgewählten Beispielen aufgezeigt, wie die beschriebenen Methoden weiterentwickelt oder außerhalb des bisherigen Anwendungs- gebiets eingesetzt werden können. 7 1 Einleitung 8 2 Grundlagen der Kontinuumsmechanik Die Kontinuumsmechanik ist die Grundlage strukturmechanischer Problemstellungen. Im ersten Abschnitt 2.1 wird daher der Zusammenhang zwischen Deformationen, Ver- zerrungen und Spannungen aufgezeigt. Diese Beziehungen führen – zusammen mit den statischen und geometrischen Randbedingungen – auf das elastostatische Randwert- problem, das in Abschnitt 2.2 von der starken in die schwache Form überführt wird. In Abschnitt 2.3 wird auf die grundlegende Vorgehensweise bei der Methode der fini- ten Elemente eingegangen. Hierbei wird die schwache Form des kontinuierlichen Pro- blems zunächst in eine diskrete Form überführt. Auch das diskretisierte Problem ist aufgrund nichtlinearer Zusammenhänge noch nicht direkt lösbar, weshalb es zunächst linearisiert und anschließend mit inkrementell-iterativen Methoden gelöst wird. Das fol- gende Kapitel hat jedoch keineswegs den Anspruch, einen umfassenden Überblick über die Grundlagen der Kontinuumsmechanik zu bieten. 2.1 Nichtlineare Festkörpermechanik Die Festkörpermechanik ist das Teilgebiet der Kontinuumsmechanik, das sich auf die Be- schreibung der Bewegung deformierbarer Körper konzentriert. Die weiteren Teilgebiete der Kontinuumsmechanik – die Strömungslehre und die Gastheorie – spielen für diese Arbeit keine Rolle und werden daher nicht näher betrachtet. Neben einer komprimierten Übersicht über die grundlegenden kontinuumsmechanischen Gleichungen soll das folgen- de Kapitel vor allem eine Einführung in die in dieser Arbeit verwendete Notation bieten. Für ein tiefer gehendes Studium der Kontinuumsmechanik sei auf die Lehrbücher von Holzapfel (2000), Bonet und Wood (1999), Altenbach und Altenbach (1994) und Ibrahimbegovic (2009) verwiesen, an denen sich dieser Abschnitt orientiert. 9 2 Grundlagen der Kontinuumsmechanik e3 e2 e1 χ X(t0) x(t) u(t) des Körpers B Ω0 Ωt Γ γ P P dX dx zum Zeitpunkt t0 Referenzkonfiguration des Körpers B zum Zeitpunkt t Momentankonfiguration Abbildung 2.1: Kinematik eines Festkörpers B in der Referenz- und der Momentankon- figuration 2.1.1 Bewegung, Kinematik und Verzerrung Die Abbildung 2.1 zeigt einen Körper B im euklidischen Raum R 3 zum Zeitpunkt t0 in seiner Referenzkonfiguration – die in dieser Arbeit immer der Ausgangskonfiguration entspricht – sowie zum Zeitpunkt t in der aktuellen bzw. Momentankonfiguration. Das vom Körper B eingenommene Gebiet zum Zeitpunkt t0 wird mit Ω0 bezeichnet, zum Zeitpunkt t mit Ωt . Alle weiteren Größen werden in der Referenzkonfiguration groß und in der Momentankonfiguration klein geschrieben. Der Rand des Gebietes (Γ bzw. γ) wird, je nach äußerer Beanspruchung, weiter unterteilt in den Dirichlet-Rand (ΓD bzw. γD) und den Neumann-Rand (ΓN bzw. γN). Die Lage eines materiellen Punkts P im Gebiet Ω0 ist in seinem Ausgangszustand durch den Ortsvektor X beschrieben. Mithilfe der Abbildungsfunktion χ lässt sich die Lage des Punkts P in der Momentankonfiguration beschreiben als x = χ (X,t) . (2.1) Damit ergibt sich der Verschiebungsvektor u des materiellen Punkts P aus der Differenz der beiden Ortsvektoren zu u = x − X. (2.2) 10 2.1 Nichtlineare Festkörpermechanik Zur Beschreibung der Deformationen wird der materielle Deformationsgradient F ein- geführt. Er ergibt sich aus der partiellen Ableitung des aktuellen Ortes x nach dem Ortsvektor in der Referenzkonfiguration X: F = ∂x ∂X = ∂χ (X,t) ∂X . (2.3) Er lässt sich auch geometrisch als Abbildung eines differenziellen Linienelements der Referenzkonfiguration dX auf ein differenzielles Linienelement der Momentankonfigu- ration dx interpretieren: dx = F dX. (2.4) Der Deformationsgradient liefert eine vollständige Beschreibung der Deformation eines Körpers. Er ist jedoch zur Beschreibung der Verzerrungen ungeeignet, da er auch Starr- körperbewegungen enthält. Dies führt bei Verwendung von F als Verzerrungsmaß zu unphysikalischen Spannungen bei reinen Translationen oder Rotationen. Um dies zu vermeiden, wird in der vorliegenden Arbeit der Green-Lagrange’sche Verzerrungsten- sor E als Verzerrungsmaß gewählt. Dieser ist besonders für große Verschiebungen bei gleichzeitig kleinen Verzerrungen geeignet: E = 1 2 ( F ⊺ F − I ) = 1 2 (CCG − I) . (2.5) Der Verzerrungstensor E ist ein richtungsunabhängiger und symmetrischer Tensor zwei- ter Stufe, der quadratisch in den Verschiebungen u ist. Er lässt sich sowohl aus dem Deformationsgradienten als auch aus dem rechten Cauchy-Green-Tensor CCG und dem Identitätstensor I ermitteln. Tiefergehende Ausführungen zur Ermittlung von CCG fin- den sich u. a. in Marsden und Hughes (1994). 2.1.2 Spannungsmaße und Konstitutivgleichungen Wirken äußere Lasten auf einen Körper B, werden in seinem Inneren Spannungen er- zeugt. Zur Beschreibung dieser Spannungen können verschiedene Spannungsmaße her- angezogen werden. Die Cauchy-Spannungen beziehen die aktuelle Kraft auf die echte, verformte Fläche in der Momentankonfiguration. Sie werden daher häufig auch als wahre Spannungen bezeichnet. Der Cauchy-Spannungsvektor t ist nach dem Cauchy-Theorem über den Cauchy-Spannungstensor σ und die Einheitsnormale der Schnittfläche n bzgl. 11 2 Grundlagen der Kontinuumsmechanik der aktuellen Konfiguration definiert zu t = σn mit σ = σ ⊺ . (2.6) Zur Beschreibung der Spannungen in der Referenzkonfiguration wird ein weiteres Span- nungsmaß benötigt, die 2. Piola-Kirchhoff-Spannungen S, auch PK2-Spannungen ge- nannt. Sie lassen sich mithilfe des materiellen Deformationsgradienten F und dessen Determinante – häufig auch als Jacobi-Determinante J bezeichnet – als Funktion der Cauchy-Spannungen darstellen: S = JF−1σF−⊺ . (2.7) Der 2. Piola-Kirchhoff-Spannungstensor stellt das zu den Green-Lagrange’schen Verzer- rungen (Gleichung (2.5)) energetisch konjugierte Spannungsmaß dar. Der Zusammenhang zwischen Spannungen und Verzerrungen wird in den Konstitutiv- oder Werkstoffgleichungen beschrieben. Das in dieser Arbeit verwendete hyperelastische Saint-Venant-Kirchhoff-Materialmodell stellt eine Erweiterung der linearen Elastizität dar, bei der die Green-Lagrange’schen Verzerrungen E sowie die PK2-Spannungen S die linearen Verzerrungen und Spannungen ersetzen. Hyperelastische Materialimodelle sind nicht pfadabhängig. Die Spannungen lassen sich über deren Potential, in diesem Fall auch als Verzerrungsenergiedichte W (E) bezeichnet, ausdrücken: S = ∂W (E) ∂E . (2.8) Für das hier verwendete Saint-Venant-Kirchhoff-Materialmodell berechnet sich die Ver- zerrungsenergiedichte wie folgt: W (E) = 1 2 λL ( tr(E) )2 + µL tr(E : E). (2.9) Die Ableitung nach den Verzerrungen liefert schließlich die PK2-Spannungen S = λL tr(E)I + 2µLE (2.10) mit den beiden Lamé-Konstanten λL und µL, die sich aus dem Elastizitätsmodul E und der Querdehnzahl ν berechnen lassen: λL = Eν (1 + ν)(1 − 2ν) , µL = G = E 2(1 + ν) . (2.11) 12 2.1 Nichtlineare Festkörpermechanik Die zweite Lamé-Konstante µL entspricht dabei dem Schubmodul G. Der Zusammenhang zwischen Spannungen und Verzerrungen lässt sich auch über den vierstufigen Materialtensor C ausdrücken, der bei St.-Venant-Kirchhoff-Material kon- stant ist: S = C : E (2.12) mit C = ∂S(E) ∂E = λL(I ⊗ I) + 2µLIsym (2.13) Hierbei steht Isym für den vierstufigen, symmetrischen Identitätstensor. Im allgemeinen Fall besitzt der Materialtensor C 81 unabhängige Einträge. Aufgrund der Symmetrie von S und E lässt sich die Zahl der unabhängigen Einträge auf 21 redu- zieren. Wird wie in der vorliegenden Arbeit ausschließlich isotropes Materialverhalten angenommen, reduziert sich die Zahl der unabhängigen Einträge weiter von 21 auf die zwei Lamé-Konstanten. Weitere Modelle zur Beschreibung des Materialverhaltens sowie ausführlichere Erläute- rung finden sich unter anderem in Simo und Hughes (1998), Holzapfel (2000) und Belytschko u. a. (2014). 2.1.3 Gleichgewicht Die kinematischen Gleichungen stellen den Zusammenhang zwischen den Verschiebun- gen und den Verzerrungen her, die Konstitutivgleichungen die Verbindung von Verzer- rungen und Spannungen. Die nun folgenden Gleichgewichtsbedingungen verknüpfen die äußeren Lasten mit den Spannungen im Inneren des Körpers B. Für statische Problemstellungen kann das Gleichgewicht allein aus der Impulsbilanz hergeleitet werden. Beschleunigungen werden hierbei vernachlässigt. div σ + ρb = 0 (2.14) Der Ausdruck ρb entspricht dabei den auf den Körper wirkenden Volumenkräften, div σ beschreibt die internen Kräfte in der Momentankonfiguration. In der Referenz- konfiguration lassen sich die internen Kräfte aus der Divergenz des 1. Piola-Kirchhoff- Spannungstensors P berechnen, welcher sich durch Multiplikation mit dem inversen De- formationsgradienten in den 2. Piola-Kirchhoff-Spannungstensors S = F−1P überführen 13 2 Grundlagen der Kontinuumsmechanik lässt. Somit lautet die Gleichgewichtsbedingung in der Referenzkonfiguration Div(FS) + ρ0B0 = 0. (2.15) 2.1.4 Randbedingungen Neben den drei im Gebiet Ω0 aufgestellten Bedingungen für die Kinematik (Gl. (2.5)), die Werkstoffbeziehung (Gl. (2.12)) und das Gleichgewicht (Gl. (2.15)) müssen nun noch die Bedingungen für den Rand Γ bestimmt werden. Hierbei wird zwischen den statischen (Spannungs- bzw. Neumann-) Randbedingungen und kinematischen (Verschiebungs- bzw. Dirichlet-) Randbedingungen unterschieden. Auf dem Neumannrand ΓN muss der Spannungsvektor T mit dem aufgebrachten Spannungsvektor T̂ im Gleichgewicht ste- hen. Über das Cauchy-Theorem in der Referenzkonfiguration und die Überführung des 2. in den 1. Piola-Kirchhoff-Tensor ergibt sich die Neumann-Randbedingung zu FSN = T̂, (2.16) mit dem Normalenvektor N auf dem Neumannrand ΓN. Die Verschiebungen u am Dirichlet-Rand ΓD müssen den vorgeschriebenen Verschiebun- gen û entsprechen: u = û. (2.17) Sind die Verschiebungen am Dirichlet-Rand gleich null, spricht man auch von homoge- nen Dirichlet-Randbedingungen. Dies entspricht der Situation an einer Lagerung. Bei inhomogenen Dirichlet-Randbedingungen wird eine vorgeschriebene Verschiebung un- gleich null auf den Rand aufgebracht. Inhomogene Dirichlet-Randbedingungen werden im Ingenieurwesen häufig als Verschiebungslastfall bezeichnet. 2.2 Elastostatisches Randwertproblem Die Gleichungen (2.16) und (2.17) komplettieren die Grundgleichungen des elastosta- tischen Randwertproblems (siehe (2.18)). Sie sind an jedem Materialpunkt erfüllt und liegen somit in der starken Form vor. Sie sind hier noch einmal übersichtlich für die Referenzkonfiguration zusammengefasst: 14 2.2 Elastostatisches Randwertproblem Kinematik: E = 1 2 (F ⊺ F − I) in Ω0 Werkstoff: S = C: E in Ω0 Gleichgewicht: Div(FS) + ρ0B0 in Ω0 Neumann-Randbedingungen: FSN = T̂ auf ΓN Dirichlet-Randbedingungen: u = û auf ΓD (2.18) In der starken Form lassen sich die Differentialbeziehungen nur für wenige Sonderfäl- le analytisch lösen. Wird die starke Form jedoch in eine integrale Form überführt – die sogenannte schwache Form – kann eine approximierte Lösung mit geeigneten Nähe- rungsverfahren, etwa der Finite-Elemente-Methode, bestimmt werden. Zur Überführung der starken Form in die schwache Form kommt die Methode der gewichteten Residuen zum Einsatz. Dabei wird zunächst eine Testfunktion δu eingeführt, die die homogenen Verschiebungsrandbedingungen erfüllen muss. Multipliziert man die Gleichgewichtsbedi- ungung (2.15) und die Kraftrandbedingungen (2.16) mit δu und integriert den Ausdruck über das Gebiet Ω0 bzw. den Rand ΓN, erhält man den Ausdruck ∫ Ω0 (Div(FS) + ρ0B0) · δu dΩ0 + ∫ ΓN (T̂ − FSN) · δu dΓN = 0 (2.19) und damit einen Ausdruck für die virtuelle Arbeit. Stellt δu ein virtuelles Verschiebungsfeld dar, ergibt sich aus der Methode der gewich- teten Residuen das Prinzip der virtuellen Verschiebungen (PvV). Mehrmalige Umfor- mungen, die z.B. in Belytschko u. a. (2014) ausführlich beschrieben sind, führen auf folgenden Ausdruck für das PvV: δΠPvV(u) = ∫ Ω0 δE: S dΩ0 ︸ ︷︷ ︸ −δΠint − ∫ Ω0 δu · B0 dΩ0 − ∫ ΓN δu · T̂ dΓN ︸ ︷︷ ︸ δΠext = 0 (2.20) Für Problemstellungen, die ein Potential besitzen, kann das PvV auch mittels des Prin- zips vom Minimum der potentiellen Energie (PMPE) hergeleitet werden. Die Variation des Potentials 15 2 Grundlagen der Kontinuumsmechanik Π(u) = ∫ Ω0 Wint(E) dΩ0 − ∫ Ω0 ρ0b · u dΩ0 − ∫ ΓN T̂ · u dΓN (2.21) führt nach mehreren Umformungen ebenfalls auf den Ausdruck in Gl. (2.20). Neben dem hier vorgestellten Einfeldfunktional, bei dem das Verschiebungsfeld u als einzige Primärvariable auftritt, wird in dieser Arbeit auch eine Elementformulierung verwendet, die auf einem Mehrfeldfunktional basiert. Für das achtknotige Volumenscha- lenelement Shell 10 wird das Funktional nach Simo und Rifai (1990) verwendet, das ei- ne modifizierte Variante des Prinzips von Hu-Washizu (Hu (1955) und Washizu (1955)) darstellt. Hierauf näher einzugehen, würde jedoch den Rahmen dieses Kapitels spren- gen. Daher sei an dieser Stelle auf Irslinger (2013) und Bischoff (1999) verwiesen. Alle weiteren in der Arbeit verwendeten Elemente basieren auf einem Einfeldfunktio- nal. Aus diesem Grund wird auf weitere Mehrfeldfunktionale nicht tiefer eingegangen. Ausführliche Herleitungen der Prinzipien finden sich für den linearen Fall z. B. in von Scheven und Bischoff (2016) und für den nichtlinearen Fall u. a. in Belytschko u. a. (2014). 2.3 Lösungsverfahren Eine analytische Lösung des elastostatischen Randwertproblems lässt sich nur für weni- ge Sonderfälle finden, weshalb in der Regel numerische Näherungsverfahren zum Einsatz kommen. In der Strukturmechanik ist dies häufig die Finite-Elemente-Methode (FEM), die das kontinuierliche Problem in eine diskrete Form überführt. Dazu wird das konti- nuierliche Primärvariablenfeld durch eine Näherung ersetzt, die mithilfe einer endlichen Anzahl n unbekannter Knotenwerte gebildet wird. 2.3.1 Diskretisierung im Raum Bei der FEM wird zuerst das kontinuierliche Gebiet Ω in eine endliche Zahl finiter Elemente unterteilt, die an Knoten miteinander verbunden sind, an denen wiederum die Freiheitsgrade definiert sind: Ω ≈ Ωh = nele⋃ e=1 Ωe. (2.22) 16 2.3 Lösungsverfahren Hierbei steht der obere Index h für diskretisierte Größen, der Index e für Größen eines Elements e, nele bezeichnet die Anzahl der Elemente, ⋃ ist der Assemblierungsopera- tor. Um Knotenwerte im Element interpolieren zu können, werden Formfunktionen Ni defi- niert, die in der Matrix der Formfunktionen Ne elementweise zusammengefasst werden. Ein Ansatz für die diskretisierte Geometrie lässt sich somit durch den Ortsvektor der Knoten Xi und den Formfunktionen aufstellen. Aus mathematischer Sicht wird dabei ein System partieller Differentialgleichungen in ein System algebraischer Gleichungen überführt. Neben der Geometrie werden auch das kontinuierliche Verschiebungsfeld u sowie die virtuellen Verschiebungen δu durch die Formfunktionen Ne und die diskreten Vektoren der Element-Knotenverschiebungen d innerhalb eines Elements approximiert: uh e = Ned, (2.23) δuh e = Neδd. (2.24) Werden für die Approximation der Verschiebungen die gleichen Formfunktionen gewählt wie für die Approximation der Geometrie, spricht man vom isoparametrischen Konzept. Wenn für den Ansatz der virtuellen Verschiebungen die gleichen Formfunktionen ver- wendet werden wie für physikalische Verschiebungen, spricht man auch vom Bubnov- Galerkin-Verfahren. Beide Prinzipien sind Grundlage vieler verschiebungsbasierter Ele- mentformulierungen, wie sie auch in dieser Arbeit verwendet werden. Wenn für die Beschreibung der Geometrie und der Verschiebungen die aus der CAD-Modellierung be- kannten Spline-Funktionen verwendet werden, spricht man zudem von isogeometrischen finiten Elementen. Um die Lesbarkeit der Gleichungen in diesem Abschnitt zu erleichtern, wird in der Folge bei diskreten Größen auf den Kopfzeiger (•)h verzichtet. Das PvV aus Gleichung (2.20) lässt sich damit in diskretisierter Form wie folgt schreiben: δΠPvV(D) = −δΠint PvV(D) − δΠext PvV(D) = δD ⊺ ( Fint(D) − Fext(D) ) = 0. (2.25) Hierbei ist D der globale Vektor der Knotenverschiebungen, Fint und Fext sind die Vek- toren der sog. inneren und äußeren Kräfte. Diese globalen Größen lassen sich wieder aus den Elementgrößen assemblieren: Fint = nele⋃ e=1 f e int, Fext = nele⋃ e=1 f e ext, (2.26) 17 2 Grundlagen der Kontinuumsmechanik mit f e int = ∫ Ωe B ⊺ S dΩe. (2.27) Das Gleichgewicht der inneren und äußeren Arbeiten kann nach Gleichung (2.25) auch in das Kräftegleichgewicht Fint(D) − Fext(D) = R(D) = 0 (2.28) überführt werden, da für die virtuellen Knotenverschiebungen δD eine beliebige Test- funktion gewählt werden kann, die lediglich die Verschiebungsrandbedingungen erfüllen muss. 2.3.2 Linearisierung Aufgrund der nichtlinearen Beziehung zwischen den inneren Kräften Fint und den Kno- tenverschiebungen D ist das Kräftegleichgewicht nicht direkt lösbar, weshalb zur Lösung des nichtlinearen Gleichungssystems ein iteratives Verfahren notwendig ist. Eine Mög- lichkeit hierfür ist das Newton-Raphson-Verfahren, das in der Nähe der Lösung quadra- tische Konvergenz aufweist. Hierfür muss zunächst Gleichung (2.28) linearisiert werden. Dies entspricht der Entwicklung einer Taylorreihe bis zum linearen Glied um die Lösung des i-ten Iterationsschritts: lin R(Di+1) = R(Di) + ∂R(D) ∂D ∣ ∣ ∣ ∣ ∣ Di ∆Di+1 (2.29) mit ∆Di+1 = Di+1 − Di . (2.30) In dieser Arbeit werden keine mitgehenden Lasten berücksichtigt, daher entsteht keine Abhängigkeit der äußeren Kräfte Fext von den Knotenverschiebungen. Eine anschauliche Darstellung verschiebungsabhängiger Drucklasten findet sich in Schweizerhof und Ramm (1984). Gleichung (2.28) vereinfacht sich damit zu Fint(D) − Fext = R(D) = 0 (2.31) 18 2.3 Lösungsverfahren Somit lässt sich die partielle Ableitung in Gleichung (2.29) zusammen mit Gleichung (2.20) umformulieren zu ∂R(D) ∂D = ∂Fint(D) ∂D = nele⋃ e ∫ Ωe ( B ⊺ CB + B,dS ) dΩe ︸ ︷︷ ︸ kT ︸ ︷︷ ︸ KT , (2.32) wobei (•),d eine vereinfachte Schreibweise für die partielle Ableitung auf Elementebene nach den Elementverschiebungen d ist. Der nichtlineare B-Operator B setzt die Ver- zerrungen im Element mit den Knotenverschiebungen in Beziehung, C steht für die Materialmatrix. Sie entspricht dem Materialtensor vierter Stufe C in Voigt’scher Nota- tion. Die daraus resultierende globale Tangentensteifigkeitsmatrix KT, die aus den auf Ele- mentebene berechneten Elementsteifigkeitsmatrizen assembliert wird, lässt sich weiter in drei Teile zerlegen: • die elastische Steifigkeitsmatrix Ke, • die Anfangsverschiebungssteifigkeitsmatrix Ku und • die geometrische Steifigkeitsmatrix Kg. Die geometrische Steifigkeitsmatrix Kg entspricht dabei dem zweiten Summanden im Integral in Gleichung (2.32). Sie beinhaltet die Steifigkeit, die durch Änderung der inne- ren Kräfte im Element entsteht, weshalb sie im Englischen auch als initial stress stiffness matrix bezeichnet wird. Der erste Summand in Gleichung (2.32) entspricht der Matrix Keu. Die Anfangsverschiebungssteifigkeitsmatrix Ku, die die Änderung der elastischen Steifigkeit aufgrund der Änderung der Geometrie beschreibt, lässt sich durch Ku = Keu − Ke (2.33) berechnen. Die elastische Steifigkeitsmatrix Ke beschreibt die Steifigkeit des unverform- ten Ausgangssystems und ist somit unabhängig von den Verschiebungen. Weitere In- formationen hierzu finden sich u. a. in Bischoff u. a. (2016), de Borst u. a. (2012) oder Belytschko u. a. (2014). Wie bereits erwähnt werden in dieser Arbeit keine zeitabhängigen, dynamischen Effekte berücksichtigt. Aus diesem Grund ist eine Diskretisierung in der Zeit nicht erforderlich. 19 2 Grundlagen der Kontinuumsmechanik 20 3 Phänomene der geometrisch nichtlinearen Strukturmechanik Bei strukturmechanischen Analysen werden häufig lineare Zusammenhänge zwischen Spannungen, Verzerrungen und Verschiebungen angenommen. Diese linearen Theori- en stellen häufig eine ausreichend genaue Näherung dar. So sind etwa im Bauwesen die Anforderungen an die Gebrauchstauglichkeit so streng, dass die Deformationen der einzelnen Bauteile oftmals sehr klein sind. Es genügt in diesen Fällen daher das Gleich- gewicht am unverformten System zu bilden. Doch es gibt im Bauwesen auch Tragwerke (z.B. Seilnetze), für die die Strukturantwort des mechanischen Modells nach linearer Theorie das reale Verhalten des Tragwerks nicht ausreichend genau abbilden kann. In diesen Fällen muss mit nichtlinearen Theorien gearbeitet werden. Eine Übersicht über die verschiedenen Ursachen nichtlinearen Strukturverhaltens fin- det sich in Abschnitt 3.1. Hierin wird näher auf die Hintergründe und die Unterschiede zwischen geometrischen und materiellen Nichtlinearitäten sowie Nichtlinearitäten durch Randbedingungen eingegangen. Abschnitt 3.2 beschäftigt sich zunächst mit den Möglich- keiten der Darstellung des nichtlinearen Zusammenhangs zwischen aufgebrachter Last und resultierender Verschiebung. Hierfür werden ebene und räumliche Darstellungen der primären und sekundären Gleichgewichtspfade für Systeme mit einem Lastparameter vorgestellt. In einem nächsten Schritt wird für den Fall mehrerer veränderlicher Pa- rameter die Darstellung der Last-Verschiebungs-Beziehung über Gleichgewichtsflächen beschrieben. In Abschnitt 3.3 wird zunächst ein kurzer Überblick über verschiedene Definitionen von Stabilität gegeben. Anschließend wird der in der Arbeit verwendete Begriff der Stabilität statischer Gleichgewichtszustände in der Strukturmechanik genau- er erläutert. In Abschnitt 3.4 wird auf die kritischen Punkte eingegangen. Es wird das unterschiedliche Verhalten von Durchschlags- und Verzweigungspunkten aufgezeigt, zu- dem werden auch hier die Besonderheiten bei Mehrparametersystemen vorgestellt. Im Anschluss folgt in Abschnitt 3.5 die Frage nach dem Verhalten unmittelbar nach dem 21 3 Phänomene der geometrisch nichtlinearen Strukturmechanik Erreichen eines kritischen Punkts. Neben der Einteilung nach Koiter (1967) wird in dieser Arbeit erstmals auch auf die Unterscheidung kritischer Punkte für Systeme mit inhomogenen Dirichlet-Randbedingungen in Bezug auf das Nachbeulverhalten eingegan- gen. Abschließend befasst sich Abschnitt 3.6 mit dem Einfluss von Imperfektionen auf die Stabilität von Tragwerken, auf die Eigenschaften kritischer Punkte sowie auf das Nachbeulverhalten. Die Ergebnisse aller gezeigten Beispiele in diesem Kapitel sind rein qualitativ. Sie die- nen lediglich der Beschreibung und Veranschaulichung ausgewählter Effekte. Wenn nicht anders angegeben sind bei allen Beispielen die Querschnittsfläche sowie der Elastizitäts- modul konstant. 3.1 Nichtlineare Effekte in der Strukturmechanik Die in der Festkörpermechanik auftretenden Nichtlinearitäten lassen sich nach ihrem Ursprung in drei Gruppen einteilen: Nichtlinearität geometrisch materiell durch Rand- bedingungen große Deformationen große Verzerrungen Plastizität Viskosität Schädigung geometrisch statisch Abbildung 3.1: Übersicht über nichtlineare Effekte in der Strukturmechanik, vgl. Wag- ner (2017) Von geometrischer Nichtlinearität spricht man, wenn der Zusammenhang zwischen den Verzerrungen und den resultierenden Verschiebungen nicht mehr linear ist. Diese Nicht- linearität entsteht durch Änderungen der Steifigkeit aufgrund der Geometrieänderung, die, anders als bei einer linearen Berechnung, nicht mehr vernachlässigbar klein ist. Geometrisch nichtlineare Effekte treten also bei großen Deformationen auf. Die großen Deformationen führen dazu, dass das Gleichgewicht nicht mehr am unverformten Sys- tem gebildet werden kann. Die Gleichgewichtsbedingung gilt also für den deformierten Zustand, dieser ist zu Beginn der Berechnung jedoch unbekannt. Die Lösung kann somit 22 3.1 Nichtlineare Effekte in der Strukturmechanik nicht in einem Schritt ermittelt werden, sondern muss durch einen iterativen Prozess be- stimmt werden. Hierbei wird die Gleichgewichtsbedingung in jedem Iteratationsschritt am deformierten System neu aufgestellt. Bei großen Deformationen muss man zwischen großen Verschiebungen bzw. Rotatio- nen bei kleinen Verzerrungen und bei großen Verzerrungen unterscheiden. Diese Arbeit konzentriert sich auf Fälle, bei denen die Verschiebungen und Rotationen groß, die Ver- zerrungen jedoch klein sind. Beispiele hierfür sind z. B. Schalen- oder Seilkonstruktionen. Große Verschiebungen und Rotationen in Kombination mit großen Verzerrungen, wie sie etwa bei Umformprozessen auftreten, werden hier nicht näher analysiert. Von materieller Nichtlinearität ist die Rede bei einem nichtlinearen Zusammenhang zwischen den Spannungen und den Verzerrungen. Die Nichtlinearität kann wiederum verschiedene Ursachen haben, weshalb zwischen den drei folgenden Effekten unterschie- den wird. • Gehen nach Be- und Entlastungsvorgang die Deformationen nicht vollständig zu- rück, spricht man von Plastizität. Das Materialverhalten wird auch als pfadabhän- gig bezeichnet. • Von Schädigung spricht man, wenn das Materialgesetz Schädigungen auf Mikro- ebene berücksichtigt. • Ist das Materialgesetz zudem zeit- bzw. geschwindigkeitsabhängig, spricht man von Viskosität. Ein bekannter Effekt ist hierbei das Kriechen. Es existieren daneben verschiedene Mischformen wie die nichtlineare Viskoplastizität. Einen ausführlichen Überblick über verschiedene lineare und nichtlineare Materialmo- delle und deren Implementierung findet man in Belytschko u. a. (2014). Als Drittes bleiben noch die Nichtlinearitäten, die durch Randbedingungen hervorgerufen werden. Ein Beispiel für Nichtlinearitäten durch statische Randbedingungen sind mitge- hende Lasten wie Gas- oder Flüssigkeitsdruck. Diese Lasten stehen immer senkrecht auf der Oberfläche und ändern somit in jedem Iterationsschritt ihre Richtung. Ein Beispiel für Nichtlinearitäten durch Verschiebungsrandbedingungen sind Kontaktprobleme. Auch hier sei wieder auf die Literatur verwiesen. Mitgehende Lasten werden u. a. in Wrig- gers (2001) behandelt. Eine ausführliche Auseinandersetzung mit Kontaktproblemen findet sich u. a. in Wriggers (2007). Die geometrische und die materielle Nichtlinearität zeigen sich in den ersten beiden der drei Feldgleichungen aus Gl. (2.18). Der dritte Punkt, die nichtlinearen Randbedingun- gen, entsprechen den letzten beiden Gleichungen aus Gl. (2.18). Lediglich die dritte Feldgleichung, die Gleichgewichtsbedingung, ist für den linearen und den nichtlinearen Fall identisch. In der Literatur ist vereinzelt die Rede von statischer Nichtlinearität, wenn 23 3 Phänomene der geometrisch nichtlinearen Strukturmechanik das Gleichgewicht am verformten System gebildet wird. Da die entsprechende Feldglei- chung jedoch unverändert bleibt, erscheint diese Bezeichnung unpassend und wird hier nicht weiter verwendet. Im Rahmen dieser Arbeit soll der Fokus auf geometrisch nichtlineare Effekte gelegt werden. Materielle Nichtlinearitäten sowie Nichtlinearitäten durch Randbedingungen werden nur der Vollständigkeit halber kurz erläutert. 3.2 Statische Gleichgewichtspfade und Gleichgewichtsflächen Wird der Zusammenhang zwischen der bei einem realen Versuch auf einen Probe- körper aufgebrachten Last und der daraus resultierenden Verschiebung für ein reales Experiment grafisch aufgetragen, spricht man von einer Last-Verschiebungskurve. Bei numerischen Experimenten lässt sich dieser Zusammenhang ebenfalls darstellen, man spricht dann von einem Gleichgewichtspfad. Diese beiden Kurven können identisch sein, wenn das Experiment quasi-statisch – also ohne Berücksichtigung dynamischer Effekte – durchgeführt wird. Gleichgewichtspfade in der Ebene Die Eigenschaften und unterschiedlichen Arten der Darstellung von Gleichgewichtspfa- den und -flächen sollen zunächst anhand einiger einfacher Beispiele gezeigt werden. Hier- für wird zunächst das in Abbildung 3.2 dargestellte flache, symmetrische und durch eine Einzellast F belastete Dreigelenktragwerk untersucht. Das Verhalten des Tragwerks lässt sich dabei mit nur einem Freiheitsgrad D beschreiben. Dies liegt zum einen darin be- gründet, dass sowohl das System als auch die Last symmetrisch sind. Zum anderen ist die Geometrie so gewählt, dass kein Verzweigen auftritt. Für welche Geometrien dies zutrifft, wird im Anhang A.1 analysiert und beschrieben. Im rechten Teil von Abbildung 3.2 ist der Gleichgewichtspfad dargestellt. Er entspricht der Last-Verschiebungskurve für ein verschiebungsgesteuertes Experiment. An jedem Punkt entlang dieser Kurve befindet sich das System im Gleichgewicht. Beginnt der Pfad wie hier im Ursprung, spricht man vom Primärpfad. Bei Tragwerken, die durch eine äußere Kraft F belastet sind, wird hierbei eine ausgewähl- te Verschiebung D auf der Abszissenachse angetragen, die aufgebrachte Last F = λF0 24 3.2 Statische Gleichgewichtspfade und Gleichgewichtsflächen a a a /2 D F = λF0 λF0 D Abbildung 3.2: Modell des flachen, symmetrischen Dreigelenktragwerks (links) und der dazugehörige Gleichgewichtspfad (rechts) oder der Lasterhöhungsfaktor λ auf der Ordinatenachse. Wirkt dagegegen eine Verschie- bungslast, befindet sich an jedem belasteten Knoten ein zusätzliches Lager. Auf der Abs- zissenachse wird die an den zusätzlichen Lagern aufgebrachte Verschiebung D̂ = λD̂0 angetragen, auf der Ordinatenachse die resultierende Lagerkraft F , die im Nachgang aus den Verschiebungen rückgerechnet werden kann. Werden mehrere Knoten verscho- ben, bietet es sich in vielen Fällen an, eine Resultierende zu ermitteln und für F diese Resultierende durch die Anzahl der verschobenen Knoten zu dividieren. Da bei diesem Beispiel nur ein veränderlicher Parameter λ und eine Verschiebung exis- tieren, liegt die Gleichgewichtskurve in einer Ebene. Gleichgewichtspfade im Raum Im nächsten Beispiel betrachten wir das in Abbildung 3.3 dargestellte asymmetrische Dreigelenktragwerk. Durch die Asymmetrie entstehen in diesem Beispiel nicht nur ver- tikale, sondern auch horizontale Verschiebungen. Es handelt sich um ein System mit einem Lastparameter und zwei Freiheitsgraden, der Gleichgewichtspfad liegt also im dreidimensionalen Raum. In der räumlichen Darstellung in Abbildung 3.3 erkennt man, dass es sich bei den beiden ebenen Gleichgewichtspfaden für die Verschiebungen D1 und D2 um Projektionen des Gleichgewichtspfads in zweidimensionale Unterräume handelt. Für Systemen mit mehr als zwei Freiheitsgraden ist es immer erforderlich, Projektionen des n + 1-dimensionalen Pfads in den zwei-, seltener auch den dreidimensionalen Unter- raum zu betrachten. Hierbei steht n für die Zahl der Freiheitsgrade des Systems. Je nach Wahl des Unterraums, in den die Kurve projiziert wird, können jedoch Eigenschaften des Pfads und somit Verhaltensweisen des Tragwerks nicht sofort ersichtlich sein. So können etwa plötzliche Richtungsänderungen der Kurve nur in einigen der Projektionen auftreten, während der Pfad in anderen Projektionen keine Änderung seiner Richtung 25 3 Phänomene der geometrisch nichtlinearen Strukturmechanik a + h a − h a /2 D1 λF0 D2 λF0 D1 λF0 D2 λF0 D2 D1 Abbildung 3.3: Modell des flachen, asymmetrischen Dreigelenktragwerks (links oben) und der dazugehörige Gleichgewichtspfad für die Freiheitsgrade D1 und D2 (rechts oben & links unten), räumliche Darstellung mit Projektionen (rechts unten) aufweist. An dieser Stelle sei auf die Arbeit von Pohl (2014) verwiesen, in der Metho- den zur effizienten Berechnung komplexer Gleichgewichtspfade vorgestellt werden, die sich diese Eigenschaften zunutze machen. Das dritte untersuchte Beispiel ist ein hohes Dreigelenktragwerk, wie es in Abbildung 3.4 zu sehen ist. Die Abmessungen wurden so gewählt, dass ein sekundärer Gleichgewichts- pfad existiert (weitere Informationen hierzu siehe Anhang A.1). Ausgehend vom Punkt, an dem die beiden Pfade verzweigen, kann sich das Tragwerk sowohl in D1 als auch in D2-Richtung verformen und befindet sich dabei im Gleichgewicht. Für komplexe Tragwerke mit mehreren hundert oder tausend Freiheitsgraden existiert im Allgemeinen eine Vielzahl von Gleichgewichtspfaden. Diese können, wie hier, den Primär- oder einen anderen Gleichgewichtspfad schneiden. Es existieren aber auch iso- lierte Pfade, die keinerlei Verbindung zu anderen Pfaden besitzen. 26 3.2 Statische Gleichgewichtspfade und Gleichgewichtsflächen a a 5a /2 D1 λF0 D2 D1 λF0 D2 Abbildung 3.4: Modell des hohen, symmetrischen Dreigelenktragwerks (links) und die dazugehörige Gleichgewichtskurve (rechts) mit Primär- und Sekundär- pfad Gleichgewichtsflächen Bei allen bisherigen Beispielen wurde immer von einem Lastparameter λ ausgegangen. Im folgenden Beispiel, einem zunächst flachen Dreigelenktragwerk, wird ein zweiter Pa- rameter eingeführt. Diese zweite zu parametrisierende Größe kann dabei frei gewählt werden. Es kann sich dabei um eine zweite Last λ2F0 handeln, aber auch um eine mate- rielle oder geometrische Größe. So lassen sich beispielhaft der E-Modul in ausgewählten Elementen oder die ursprüngliche Lage einzelner Knoten parametrisieren. Im folgenden Beispiel wird die Höhe h0 + λ2h des Dreigelenktragwerks variiert. Durch den zusätzlichen Parameter λ2 ergibt sich die rechts in Abbildung 3.5 dargestellte Gleichgewichtsfläche. Sie lässt sich als eine kontinuierliche Aneinanderreihung der ein- zelnen Gleichgewichtspfade für die variierende Höhe verstehen. Die dargestellte Fläche ist jedoch nur eine Projektion der vierdimensionalen Hyperfläche (λ1, λ2, D1, D2) in den dreidimensionalen Unterraum (λ1, λ2, D1). Auf dieser Fläche können wiederum Pfa- de mit ausgewählten Eigenschaften gefunden werden. Auf der hier abgebildeten Fläche kann etwa der Pfad der kritischen Punkte, im Englischen auch als Fold Line bezeichnet, verfolgt werden. Hierzu sei auf das Kapitel 5 verwiesen, in dem genau diese Pfade näher untersucht werden. Auch bei einer Fläche, die durch einen zweiten Lastparameter entsteht, kann nach Pfa- den auf der Fläche gesucht werden, die wiederum spezielle Eigenschaften besitzen. So lassen sich etwa Pfade finden, auf denen die beim Ablaufen zu leistende Arbeit minimal ist. Andere Pfade verbinden zwei Gleichgewichtspunkte, für die jeweils eine gewünschte Geometrie vorgegeben ist, ohne dabei instabil zu werden. 27 3 Phänomene der geometrisch nichtlinearen Strukturmechanik a a h 0 λ 2 h D1 λ1F0 D2 D1 λ2 λ1F0 Abbildung 3.5: Modell eines Dreigelenktragwerks mit variabler Höhe (links) und die dazugehörige Gleichgewichtsfläche (rechts) 3.3 Stabilität statischer Gleichgewichtszustände Leonhard Euler analysiert 1744 erstmals das Knicken schlanker, gerader, druckbean- spruchter Stäbe. Noch heute sind die vier Euler’schen Knickfälle Bestandteil der Lehre für angehende Ingenieure. Im Jahr 1788 befasst sich Joseph-Louis Lagrange in „Mécha- nique Analytique“ ebenfalls mit der Frage nach Stabilität. Anders als Euler, der sich mit der Stabilität statischer Systeme auseinandersetzt, befasst sich Lagrange mit der Sta- bilität dynamischer Systeme. Die Definition von stabilem Gleichgewicht geht dabei auf das Minimum der potentiellen Energie zurück. Im Jahr 1892 liefert Alexander Michailo- witsch Liapunov in seiner Dissertation eine Stabilitätstheorie nichtlinearer dynamischer Systeme. Danach ist ein System bzw. ein Prozess dann stabil, wenn die Lösung eines gestörten Problems nicht über alle Grenzen hinweg von der Lösung des ungestörten Systems abweicht. Weitere Ausführungen zu den genannten Theorien finden sich u. a. in Loria und Panteley (2017), Timoschenko (1936), Kollbrunner und Meister (1955), Bürgermeister und Steup (1966) sowie Thompson und Hunt (1984). Die vorliegende Arbeit befasst sich mit der Frage der Stabilität von Gleichgewichts- zuständen. Auch hierauf kann die Definition nach Liapunov übertragen werden. Ein System befindet sich demnach in einem stabilen Gleichgewichtszustand, wenn es nach einer hinreichend kleinen Störung wieder in seine ursprüngliche Gleichgewichtslage zu- rückkehrt. Ist ein System durch konservative Kräfte belastet, lässt sich das Verhalten gut am Beispiel der Kugel in Abbildung 3.6 beschreiben. Wird die Kugel im linken Bild ausgelenkt, kehrt sie wieder in den Ausgangszustand zurück. Hierbei lässt sich die Kurve auch als potentielle Energie betrachten. Weist die potentielle Energie ein lokales Minimum auf, befindet sich das System in einem stabilen Gleichgewicht. 28 3.4 Kritische Punkte und kritische Pfade Abbildung 3.6: Stabilität eines Gleichgewichtszustand am Beispiel einer Kugel. Stabiler Gleichgewichtszustand (links), indifferenter Gleichgewichtszustand (mit- tig) und instabiler Gleichgewichtszustand (rechts) Ein Gleichgewichtszustand ist hingegen instabil, wenn bereits eine kleine Störung dazu führt, dass sich das System unaufhaltsam von der ursprünglichen Gleichgewichtslage entfernt. Dies lässt sich an der rechten Darstellung in 3.6 erkennen. Lenkt man die Kugel nur minimal aus, rollt diese nicht wieder zurück in ihre ursprüngliche Ausgangsposition. Interpretiert man die Kurve wieder als potentielle Energie, so weist sie in diesem Fall ein lokales Maximum auf. Der Fall des indifferenten Gleichgewichts ist im mittleren Bild zu sehen und bildet den Übergang zwischen stabilem und instabilem Gleichgewicht. Nach einer kleinen Auslen- kung der Kugel kehrt diese nicht wieder in die Ausgangslage zurück. Die Auslenkung wächst allerdings auch nicht über alle Grenzen hinweg an. 3.4 Kritische Punkte und kritische Pfade Von besonderem Interesse sind die Punkte auf dem Gleichgewichtspfad oder der Gleich- gewichtsfläche, an denen die Grenze des stabilen Gleichgewichts erreicht wird. Diese kritischen Punkte zeichnen sich – frei nach Euler – dadurch aus, dass für die gegebene Belastung neben dem aktuellen Gleichgewichtszustand mindestens ein weiterer, unmit- telbar benachbarter Gleichgewichtszustand existiert. In Anlehnung an Mang und Hofstetter (2000) lassen sich kritische Punkte nach ihrer Versagensart im Stabilitätsfall in drei Arten unterteilen: Durchschlagen, Verzwei- gen des Gleichgewichts und das Erreichen der Traglast. Da es sich im dritten Fall, dem Erreichen der Traglast, definitionsgemäß um ein materiell nichtlineares Phänomen han- delt, konzentriert sich diese Arbeit auf die Untersuchung von Durchschlagspunkten, im Englischen als limit points bezeichnet, und Verzweigungspunkten, im Englischen bifur- cation point genannt. Wenn im Folgenden von kritischen Punkten die Rede ist, bezieht sich dies ausschließlich auf Durchschlags- und Verzweigungspunkte. Kritische Punkte kennzeichnen sich nach Euler also durch des Vorhandensein eines un- mittelbar benachbarten Gleichgewichtszustands bei gleichbleibender Last. Dies bedeu- 29 3 Phänomene der geometrisch nichtlinearen Strukturmechanik λF0 Durchschlagspunkt D F = λF0 Abbildung 3.7: Durchschlagspunkte für Kraftlastfall am Beispiel eines Dreigelenktrag- werks (qualitativ) tet, dass an einem kritischen Punkt die Steifigkeit in mindestens eine Richtung gleich null sein muss. Dies hat zur Folge, dass die Tangentensteifigkeitsmatrix KT an einem kritischen Punkt singulär wird und nicht mehr positiv definit ist. Das unterschiedliche Verhalten beim Erreichen bzw. Überschreiten von kritischen Punk- ten lässt sich in den Abbildungen 3.7, 3.8 und 3.9 erkennen. Die Art des kritischen Punkts hängt dabei von der Richtung ab, in die der Eigenvektor am kritischen Punkt orientiert ist. Hierauf wird im Abschnitt 4.3.3 näher eingegangen. 3.4.1 Durchschlagspunke bei Kraftlastfällen In Abbildung 3.7 ist ein flaches Dreigelenktragwerk mit dem dazugehörigen qualitativen Gleichgewichtspfad dargestellt. Hierbei ist die aufgebrachte Last λF0 am mittleren Kno- ten über die resultierende vertikale Verschiebung des Knotens aufgetragen. Neben der in schwarz dargestellten Ausgangsgeometrie sind mehrere deformierte Zustände dargestellt. Die zu den Deformationszuständen gehörenden Punkte auf der Gleichgewichtskurve sind in den jeweils passenden Farben gekennzeichnet. Am ersten Durchschlagspunkt, hier in blau dargestellt, weist der Gleichgewichtspfad eine horizontale Tangente auf. Diese ho- rizontale Tangente lässt sich auch anschaulich interpretieren: Das System besitzt keine Steifigkeit in Richtung der aufgebrachten Kraftlast λF0. Ohne λ weiter zu steigern lässt sich das System also in einen infinitesimal benachbarten Gleichgewichtszustand führen. Dies entspricht der Definition des indifferenten Gleichgewichts aus Abschnitt 3.3. Bei einem lastkontrollierten Experiment ist das Tragwerk nach dem Erreichen des Durch- schlgspunkts instabil und schlägt dynamisch durch. Nach dem Ausschwingen kommt es auf einem stabilen Ast des Gleichgewichtspfads in der grau gezeichneten Konfiguration zur Ruhe. Wird die Last im Anschluss wieder zurückgefahren, ist das System bis zum Erreichen des pink gekennzeichneten Durchschlagspunkts stabil und schlägt von dort in eine hier nicht dargestelle Gleichgewichtslage zurück. 30 3.4 Kritische Punkte und kritische Pfade Durchschlagspunkt F λD0 D = λD0 D1 F D1 F D1 λD0 Abbildung 3.8: Durchschlagspunkte für Verschiebungslastfall am Beispiel eines modifi- zierten Dreigelenktragwerks (qualitativ) 3.4.2 Durchschlagspunkte bei Verschiebungslastfällen Abbildung 3.8 zeigt die Modifikation eines Dreigelenktragwerks. Es ist am rechten Kno- ten mit einer inhomogenen Verschiebungsrandbedingung belastet. Die dazugehörigen Gleichgewichtspfade sind daneben und darunter abgebildet. Sie zeigen die resultierende vertikale Lagerkraft F am rechten Knoten, angetragen über die aufgebrachte Verschie- bung λD0 bzw. die vertikale Verschiebung D1 am freien Knoten. Anders als bei Kraft- lastfällen, bei denen die aufgebrachte Last λF0 über die resultierende Verschiebung D angetragen wird, tritt der erste Durchschlagspunkt folglich auch nicht an der ersten horizontalen, sondern an der ersten vertikalen Tangente des Gleichgewichtspfads auf. Ähnlich wie die horizontale Tangente bei einem Kraftlastfall lässt sich auch die vertikale Tangente bei einem Verschiebungslastfall anschaulich interpretieren: am Durchschlags- punkt besitzt das System keine Flexibilität in Richtung der aufgebrachten Verschiebung D. Es existiert also in der unmittelbaren Umgebung für die gleiche aufgebrachte Ver- schiebung D ein weiterer Gleichgewichtszustand. Bei einem verschiebungsgesteuerten Experiment wird das System mit Erreichen des Durchschlagspunkts instabil. Auch hier schlägt das System dynamisch durch und kommt auf einem stabilen Ast des Gleichgewichtspfads zur Ruhe. 3.4.3 Verzweigungspunke In Abbildung 3.9 ist ein hohes Dreigelenktragwerk mit dazugehörigem Gleichgewichts- pfad dargestellt. Das Dreigelenktragwerk wird durch eine äußere Kraftlast F = λF0 31 3 Phänomene der geometrisch nichtlinearen Strukturmechanik λF0 D1 D2 λF0 Verzweigungspunkt Abbildung 3.9: Verzweigungspunkt am Beispiel eines Dreigelenktragwerks (qualitativ) belastet. Zum besseren Verständnis ist der Gleichgewichtspfad im Raum gegeben. Hier- bei stehen D1 für die vertikale und D2 für die horizontale Verschiebung des belasteten Knotens. Anders als bei dem niedrigen Dreigelenktragwerk in Abbildung 3.7 wird hier der Primär- pfad vor dem ersten Durchschlagspunkt – also vor dem ersten lokalen Maximum – von einem Sekundärpfad gekreuzt. Man spricht daher von einem Verzweigungspunkt. Auch an einem Verzweigungspunkt wird die Steifigkeitsmatrix des diskretisierten Systems sin- gulär. Anders als bei Durchschlagspunkten verschwindet bei Verzweigungspunkten die Steifigkeit jedoch nicht in Richtung der aufgebrachten Last, sondern orthogonal dazu. Ist das System perfekt, wird sich das Tragwerk entsprechend dem Primärpfad weiter verformen. Minimale Störungen der Gleichgewichtslage reichen jedoch aus, um vom Pri- märpfad in den Sekundärpfad (schwarz dargestellt) zu wechseln. Anders als bei echten Durchschlagsproblemen, bei denen sich der kritische Punkt durch die starke Änderung der Steigung der Gleichgewichtskurve ankündigt, treten Verzwei- gungsprobleme ohne Vorwarnung auf. Typische Beispiele hierfür aus der Praxis sind das Beulen von Zylindern, wie z. B. Silos, oder das Ausknicken druckbelasteter Stützen. 3.4.4 Kritische Pfade bei Mehrparametersystemen Wie schon in Abschnitt 3.2 gezeigt, lässt sich für Systeme mit mehreren Parametern die Gleichgewichtsbeziehung als Fläche darstellen (vgl. Abb. 3.5). So wie aus einem Gleich- gewichtspfad im zweidimensionalen Raum bei Hinzunahme eines weiteren Parameters (in diesem Fall die Höhe des Dreigelenktragwerks) eine Fläche im dreidimensionalen Raum entsteht, so wandelt sich der kritische Punkt zum kritischen Pfad. Wir betrachten erneut das Dreigelenktragwerk mit variabler Höhe. In Abbildung 3.10 ist die schon bekannte Gleichgewichtsfläche dargestellt, bei der der Lasterhöhungsfaktor λ1 32 3.5 Nachbeulverhalten a a h 0 λ 2 h D1 λ1F0 D2 D1 λ2 λ1F0 VerzweigungspunkteDurchschlagspunkte Abbildung 3.10: Modell eines Dreigelenktragwerks mit variabler Höhe (links) und die dazugehörige Gleichgewichtsfläche (rechts) mit kritischen Pfaden über die vertikale Ve