n ò ªÐp|É ¬Òçs¨ü¡ Ä{ h e÷ ¢) Ì óc d¿ Ø 1.pdf · elastostatik und -dynamik in...
Post on 06-Feb-2018
219 Views
Preview:
TRANSCRIPT
FINITE-ELEMENT-METHODE
Teil I
(FEM-I)
U. GABBERT
Lehrstuhl für Numerische MechanikInstitut für Mechanik
Nur zum Gebrauch an der Otto-von-Guericke-Universität Magdeburg bestimmt !
Vorbemerkungen III
VorbemerkungenDas Ziel des Teils I der Vorlesung zur FINITE-ELEMENT-METHODE besteht darin,wichtige theoretische Grundlagen dieser Methode zur computergestützten Bauteilbe-rechnung und deren qualifizierter Anwendung im Ingenieurwesen zu vermitteln, wobeider Schwerpunkt auf statischen und dynamischen Problemen der Strukturmechanikliegt. Einleitend wird ein kurzer Überblick über den Stand der Entwicklung und derAnwendung der FEM gegeben. Kerngedanken der FEM werden an Hand einfacherModellprobleme (Stab, Balken) dargestellt. Nach einer kurzen Einführung in die Ener-giemethoden der Mechanik werden die Grundgleichungen der FEM für Probleme derElastostatik und -dynamik in allgemeiner Form abgeleitet. Es werden finite Elementefür wichtige Modellklassen der Mechanik (2D, 3D, Schalen), das isoparametrischeElementkonzept und die Substruktur-Superelement-Technik ausführlich behandelt. EinKapitel zur Fehlerschätzung und Netzadaption, und ein Kapitel zu Fragen der Modell-bildung beschließen den ersten Teil der Vorlesung.
Es werden weiterhin auch Grundkenntnisse über den prinzipiellen Aufbau, dieStruktur und die programmtechnische Realisierung von kommerzieller Berechnungs-software vermittelt. Im Rahmen eines individuellen Semesterbelegs werden dieseKenntnisse vertieft. Der Beleg umfaßt sowohl die Erarbeitung der theoretischenGrundlagen als auch die Erstellen eines kompletten Computerprogramms, einschließ-lich Testung, Dokumentationen und Demonstration. Als Programmiersprache werdenFortran’99 oder C empfohlen.
Praktische Erfahrungen in der Anwendung der FEM werden in einem Rechner-praktikum an Hand von Testbeispielen - überwiegend statische und dynamische Pro-bleme für unterschiedliche Modellklassen der Mechanik - vertieft. Dazu wird daskommerzielle Programmsystems COSAR1 eingesetzt wird. Für individuelle Übungenstehen im Rechnerlabor des Institutes für Mechanik auch andere kommerzielle Pro-grammsysteme zur Verfügung
Der Vorlesungsinhalt beschränkt sich im Teil FEM I auf linear-elastische Pro-bleme. Erweiterungen auf allgemeine Feldprobleme, die Berechnung von Mehrfeld-problemen (für intelligente oder adaptive Struktursysteme) und nichtlineare Aufgabensind Gegenstand des Teils II der Vorlesung.
Umfang und Voraussetzungen
Der Vorlesungsumfang beträgt 4 SWS (2V, 2Ü&P). Die Vorlesung gehört zumPflichtprogramm der Studienrichtung Angewandte Mechanik im Studiengang Maschi-nenbau. Sie ist offen für alle Studienrichtungen, wobei mindestens der erfolgreicheAbschluß des Grundstudiums (Maschinenbau, Mechatronik oder Elektrotechnik) er-forderlich ist. Weiterhin sind Grundkenntnisse auf den Gebieten Elastizitätstheorie,Variationsrechnung, Computeranwendungen in der Mechanik (Programmierung),Computer-Numerik und Flächentragwerke erwünscht. Fehlende Grundlagen sind ge-gebenenfalls individuell nachzuholen, wobei dafür gesonderte Vorlesungsskripte bzw.Literaturempfehlungen zur Verfügung gestellt werden können.
1 COSAR ist ein Warenzeichen der Forschungsgesellschaft für Technische Mechanik FEMCOS mbH, siehe auchhttp://femcos.de
IV Literaturempfehlungen
Leistungsbewertung
Die Leistungsbewertung erfolgt auf der Grundlage einer mündlichen Prüfung (45min.), die die beiden Fächer MNMM-I2 und FEM-I umfasst, wobei der Schwerpunktder Prüfung auf dem Gebiet FEM I liegt. Als Zulassungsvoraussetzung für die Prüfunggilt der erfolgreiche Abschluß des Faches MNMM-I, die Anerkennung des Semester-beleges im Fach FEM-I (wird benotet) und das Testat über die erfolgreiche Teilnahmeam FEM-Rechnerpraktikum. Die Note des Semesterbeleges in FEM-I und die Note imFach MNMM-I bilden die Vornote für die Prüfung, die zu 50% in die Abschlussnote(MNMM-I+FEM-I) eingeht. 3
Literaturempfehlungen[1] Zienkiewicz, O.C.: Methode der finiten Element. VEB Fachbuchverlag Leipzig. 1983; Carl Han-
ser Verlag 1984 (englisches Original: The Finite Element Method, McGraw Hill 1977).[2] Bathe, K.-J.: Finite-Element-Methode. Springer 1990 (englisches Original: Finite Element Proce-
dures in Engineering Analysis. Prentice Hall, Inc. 1982).[3] Szabo, B. A.; Babuška, I.: Finite-Element-Analysis. John Wiley & Sons 1991.
Eine Auswahl deutschsprachiger Bücher[4] Schwarz, H. R.: Methode der finiten Elemente. Teubner Verlag, Stuttgart 1980.[5] Schwarz, H. R.: FORTRAN-Programme zur Methode der finiten Elemente. Teubner Verlag, Stuttgart
1980.[6] Hinton, Owen: FE-Programme für Platten und Schalen. Springer Verlag.[7] Argyris, Mlejnek: Die Methode der FEM, Bd. I, II, III. Fridr. Vieweg Verlag, Wiesbaden, 1986, 1987,
1988.[8] Autorenkollektiv: Die Methode der finiten Elemente in der Festkörpermechanik. VEB Fachbuchver-
lag, Leipzig 1982, Carl Hanser Verlag, 1982.[9] Kämmel, G., Franeck, H.-J., Recke, H. G.: Einführung in die Methode der finiten Elemente. VEB
Fachbuchverlag, Leipzig 1988.[10] Goering, H., Roos, H.-G., Tobiska, L.: Finite-Element-Methode. Akademie Verlag, Berlin 1993.[11] Link, M.: Finite Elemente in der Statik und Dynamik, Teubner-Verlag, Stuttgart, 1989.[12] Chung, T. J.: Finite Elemente in der Strömungsmechanik. VEB Fachbuchverlag, Leipzig 1982, Carl
Hanser Verlag, München 1982.
Fachzeitschriften (Auswahl)[1] Technische Mechanik[2] Int. Journal for Numerical Methods in Engineering[3] Computers and Structures[4] Structural Engineering and Mechanics[5] Computer Methods in Applied Mechanics and Engineering[6] Computational Mechnics Adavances[7] International Journal of Solids and Structures[8] Engineering with Computers[9] Finite Elements in Analysis and Design[10] Engineering Computations[11] Communications in Numerical Methods in Engineering[12] Archives of Computational Methods in Engineering[13] Finite Elements in Analysis and Design[14] Engineering Computations – International Journal for Computer Aided Engineering and Software
2 Mathematische und numerische Methoden der Mechanik – Teil I3 Auf Wunsch wird für jedes der beiden Fächer – MNMM-I und FEM-I auch eine separate Prüfungsbescheini-gung ausgestellt.
Inhaltsverzeichnis V
InhaltsverzeichnisSeite
1 Elementare Einführung in die FEM 11.1 Anwendung der FEM 11.2 Einführung in die Grundlagen der FEM 101.3 Kommerzielle FEM-Softwareprodukte 151.4 Ein Einführungsbeispiel zur Statik von Stäben 171.5 Ein Einführungsbeispiel zur Dynamik von Stäben 271.6 Ein Einführungsbeispiel zur eindimensionalen Wärmeleitung 331.7 Kleiner historischer Rückblick in die frühen Anfänge der FEM 41
2 Grundlagen der FEM für die Elastomechanik 432.1 Differentialgleichungen der linearen Elastostatik 432.2 Variationsformulierungen 482.3 Die Verschiebungsgrößenmethode 522.4 Steifigkeits-, Massen- und Belastungsmatrizen 582.5 Das Gesamtsystem 62
3 Simplex-Elemente 683.1 Das 2-Knoten Stabelement 683.2 Das 2-Knoten Balkenelement (Bernoulli-Balken) 723.3 Das 2-Knoten Rotationsscheibenelement 803.4 Das 2-Knoten Rotationsplattenelement 833.5 Das 3-Knoten Dreieckselement für Scheibenberechnungen 873.6 Das 4-Knoten Rechteckelement für Scheibenberechnungen 923.7 Das 4-Knoten Tetraederelement 963.8 Das 8-Knoten Hexaederelement 100
4 Interpolationsfunktionen als Ansatzfunktionen 1024.1 Lagrangesche Polynome für C0 stetige Elemente 1024.2 Hermitesche Polynome für C1 stetige Elemente 1044.3 Legendresche Polynome als Basis für p-Elemente 1064.4 Konstruktion von Ansatzfunktionen für finite Elemente 107
5 Das isoparametrische Elementkonzept 1115.1 Iso-, sub- und superparametrische finite Elemente 1115.2 Isoparametrische finite Elemente 113
6 Substruktur-Superelementtechnik 1176.1 Motivation und Zielstellung 1176.2 Die statische Kondensation 1176.3 Anwendung auf dynamische Probleme - Gyan-Reduktion 120
7 Modellbildung und Fehleranalyse 1237.1 Einführung 1237.2 Fehlerquellen bei FEM-Analysen 1247.3 Mathematische Grundlagen der Fehleranalyse in der FEM 128
Inhaltsverzeichnis VI
7.4 A priori Fehleranalyse 1317.5 A posteriori Fehleranalyse 1367.6 Adaptive Finite-Element-Methoden 1457.7 Zusammenfassung 146
Achtung: Ab hier wird die Vorlesung weiter aktualisiert!
8 Übersicht über wichtige Elementfamilien 1478.1 C0-stetige 2D- und 3D-Finite Elemente8.2 Platten- und Schalenelemente8.2.1 Übersicht über Elementfamilien und Ansatzfunktionen8.2.2 Die SemiLoof Schalenelemente8.3 Elemente für Faserverbundstrukturen (CLT und Erweiterungen)8.4 Koppelung unterschiedlicher Elementtypen mittels Penalty-Technik8.5 Spezial- und Sonderelemente
1.1. Anwendung der FEM 1
> Konstruktion (CAD)
> Modell
Berechnung (FEM)
Bewertung
Bild 1.1-1 Zusammenspiel von Konstruktion und Berechnung
1. Elementare Einführung in die FEM1.1 Anwendung der FEMBedeutung und Nutzen der FEMDie FEM ist für den Berechnungsingenieur inzwischen zu einem unverzichtbarenAuslegungs- und Simulationswerkzeug in nahezu allen Bereichen des Ingenieurwesensgeworden. Die Ursachen dafür liegen zum einen in den enormen Leistungssteigerun-gen auf dem Gebiet der Computertechnik und den Fortschritten, die vor allem in dercomputerorientierten Mechanik, der numerischen Mathematik sowie der Informatikund Softwareentwicklung erreicht wurden.
Die Anwendung einer FEM-Software stellt auf den ersten Blick scheinbar nurgeringe Anforderungen an die Bedienung und Nutzung (z.B. vergleichbar mit CAD-Systemen). Es lassen sich offensichtlich problemlos auch extrem komplexe und kom-plizierte Probleme berechnen ohne daß der Nutzer weitergehende Kenntnisse von dender Berechnung zu Grunde liegenden Theorie zu haben braucht. So kann man auch inrenommierten Unternehmen häufig fachlich unzureichend qualifizierten Ingenieurenerleben, die zum Beispiel physikalisch und geometrisch nichtlineare Berechnungen fürsicherheitsrelevante Bauteile durchführen ohne sich im mindesten über die den Be-rechnungen zugrundeliegende Modellannahmen oder die verwendeten numerischenLösungsmethoden im Klaren zu sein. Die Ergebnisse solchen Handels werden leidernur manchmal bei spektakuläre Schadensfälle sichtbar. Die Lösung derartiger Proble-me, die teilweise bis vor kurzem noch gar nicht möglich war, erfordert im Gegenteilein fundiertes theoretisches Wissen und gründliche Kenntnisse auf den unterschiedli-chen Fachgebieten, die zur Erzielung eines brauchbaren Ergebnisses zusammenspielenmüssen. Dazu kommt, daß sich derAnwendungsbereich der Finite-Element-Methode ständig erweiterthat und dadurch heute Lösungenund neue Erkenntnisse in nahezuallen Bereichen des Ingenier-wesens, der Physik, der Medizin,der Biologie, der Umwelttechnikusw. ermöglicht werden.
Der Schwerpunkt der An-wendungen liegen aber nach wievor auf dem Gebiet der Struktur-mechanik mit den EinsatzfeldernMaschinenbau, Bauwesen, Fahr-zeugbau, Luft- und Raumfahrt.Durch leistungsfähige Schnittstellen zu CAD-Systemen gelingt es zunehmend besser,die bereits vorhandenen geometrischen Daten zu nutzen und daraus über leistungsfähi-ge Pre-Prozessoren (automatische Netzgeneratoren) geeignete FEM-Berechnungsmodelle zu generieren (siehe Bild 1.1-1). Damit ist dieser arbeitsintensiveAnteil der FEM heute sehr stark vereinfacht worden. Entscheidend für die Qualität ei-ner Berechnung ist aber immer noch die Modellbildung!
2 1. Elementare Einführung in die FEM
Einige Beispiele für die Anwendung der FEMDie folgenden Beispiele zeigen eine kleine Auswahl von Anwendungen der FEM, dieaus Diplomarbeiten von Studenten der Studienrichtung Angewande Mechanik ent-nommen worden sind.
Im Bild 1.1-2 ist die Crash-Simulation eines PKW gezeigt. In der entsprechen-den Diplomarbeit von Jens Grabau1 ging es um die Frage, ob trotz der größeren Dickevon Bauteilen aus Aluminium im Vergleich zu Stahl noch die Nutzung von Schalen-modellen bei der Crash-Simulation zulässig ist oder ob zu 3D-Modellen übergegangenwerden muß. Systematische Untersuchungen, die Jens Grabau im Rahmen seiner Di-plomarbeit durchgeführt hat, haben ergeben, daß auch bei dickeren Aluminiumbautei-len mit Schalenmodellen eine ausreichende Ergebnisqualität bei Crash-Simulationenerreicht wird, so daß auf die wesentlich teureren 3D-Berechnungen (Modellaufberei-tung, Speicherplatzanforderungen, Rechenzeit) verzichtet werden kann.
Finite - Elemente - Modell desGesamtfahrzeugs
Verformte Fahrzeugstruktur 90 msnach dem Aufprall
Längsträger mit Teilen des Vorderwagens für t = 0 und 25 ms_________________________________________________ Quelle: VW AG aus ATZ 90 (1988) 11______
Bild 1.1-2 Crash-Simulation eines Pkw mit LS-DYNA3D
1 Grabau, Jens: Rechnerische Analyse dickwandiger Karosseriebauteile unter Crashbelastungen. Diplomarbeit,Uni Magdeburg, Institut für Mechanik, 1993. (Aufgabenstellung Mercedes-Benz AG Sindelfingen).
1.1. Anwendung der FEM 3
Torsten Mann2 hat in seiner Diplomarbeit Untersuchungen zum Stabilitätsver-halten von Airbus-Rumpfschalen durchgeführt. Dabei ging es auch um die Bewertungdes kommerziellen Finite-Element-Programmsystems ABAQUS hinsichtlich seinerLeistungsfähigkeit bei nichtlinearen Beulrechnungen unter Einbeziehung des Nach-beulverhaltens. Das Bild 1.1-3 zeigt den Aufbau einer solchen Rumpfschale, und dasBild 1.1-4 vermittelt einen Eindruck von dem Aufbau des untersuchten Rumpfseg-mentes.
Bild 1.1-3 Prinzipieller Aufbau eines Airbus Bild 1.1-4 Modell des Segmentes Rumpfsegmentes
Das Bild 1.1-5 zeigt das für die Berechnung verwendete Schalenmodell des Segmentesder Rumpfschale.
Bild 1.1-5 FEM Modell der Rumpfsegmentes Bild 1.1-6 Muster der ersten Beuleigenform
Die nahezu starren Enden der Schale (siehe Bild 1.1-5) wurden in das Modelleingeführt, um definierte Lastsituationen in die Schale einleiten zu können. Sie
2 Mann, Torsten: Berechnung des Verformungsverhaltens einer Airbus-Rumpfstruktur-Sektion unter Einbezie-hung geometrischer Nichtlinearitäten. Diplomarbeit, Uni Magdeburg, Institut für Mechanik, 1999. (Aufgaben-stellung DLR Braunschweig, Dr. Zimmermann)
4 1. Elementare Einführung in die FEM
Bild 1.1-7 Beuleigenform einer Teil- schale eines Airbus A340
entsprechen beispielsweise auch den Lagerbedingungen, die bei Beulexperimenten ineiner Versuchsanlage realisiert werden können.3
Herr Mann hat in seiner Diplomarbeitzunächst eine Reihe von einfacherenVoruntersuchungen (Beulen von einfacherenSchalenmodellen) durchgeführt, um Aussagenüber die Genauigkeit und die Qualität derBerechnungsergebnisse zu erhalten. DieseErgenisse konnten teilweise noch mit“exakten” Lösungen oder mit anderennumerischen Methoden verglichen undabgesichert werden.4 Ein Ergebnis derBeuluntersuchungen an dem komplettenRumpfschalenmodell zeigt Bild 1.1-6. Füreine einfachere Teilschale ist dieBeuleigenform im Bild 1.1-7 dargestellt. DieUntersuchungen zum Nachbeulverhalten deskompletten Rumpfschalensegmentes, wie esin Bild 1.1-5 dargestellt ist, führten teilweisezum Abbruch der Rechnungen oder zu völligunbrauchbaren Ergebnissen. Trotz Einsatz
aller im FEM-System angeboteten Möglichkeiten zur Sicherung der Konvergenz, undauch durch die Nutzung alternativer numerischer Methoden gelang es demDiplomanden nicht, eine brauchbare Lösung zu gewinnen. Weder die bei der DLRvorhandenen Fachleute noch die Hotline des Softwareanbieters konnten das Problemlösen.5
Heikko Rädiger6 hat sich in seiner Diplomarbeit mit der Beurteilung derBetriebsfestigkeit von Motor-Lagerstühlen befaßt. Um die Eingangsdaten für dieBerechnungen zur Betriebsfestigkeit (diese erfolgte mit dem ProgrammsystemFEMFAT7) zu gewinnen, wurde eine Vielzahl von FEM-Berechnungen mit ABAQUSdurchgeführt, um die Spannungen als erforderlich Basisdaten zu gewinnen. DasGesamtmodell dieser Berechnung bestehend aus Zylinderkopf, Kopfdichtung,Kurbelgehäuse, Ölwanne, Zylinderbuchsen und weiteren Anbauteilen enthält etwa170000 finite Elemente (überwiegend 3D Hexaeder- und Pentaederelement) mit etwa330000 Knotenpunkten (ca. 980000 Freiheitsgrade). Es wurden in der Arbeit mitdiesem Modell - das Bild 1.1-8 zeigt nur einen Teil dieses Modells - für 36 Stellungender Kurbelwelle (von 12 Grad bis 702 Grad) komplette Rechnungen durchgeführt und 3 Das Deutsche Forschungszentrum für Luft- und Raumfahrt (DLR) Braunschweig verfügt beispielsweise übereine entsprechende Versuchseinrichtung.4 Ein interessantes Ergebnis ergab sich bei der Nachrechnung eines dieser Testfälle mit COAR. Das dafür inCOSAR genutzte SemiLoof Schalenelement zeigte eine deutlich höhere Genauigkeit als das in ABAQUS im-plementierte Schalenelement, d.h. mit einer sehr groben Vernetzung wurden mit dem SemiLoof schon wesent-lich genauere Ergebnisse erzielt als mit dem feinsten getesteten Netz der ABAQUS Rechnung.5 Es sollte hier angemerkt werden, daß ABAQUS auf dem Gebiet der nichtlinearen Berechnugen zu den weltweitführenden FEM-Systemen zählt!6 Rädiger, Heiko: Festigkeitsbeurteilung von Motor-Lagerstühlen auf der Basis der Finite Element Methode undeinfacher Prüfstandsversuche. (Aufgabenstellung von Steyr-Daimler-Puch AG, Dr. W. Steiner) Diplomarbeit,Uni Magdeburg, Institut für Mechanik, 19997 FEMFAT - Finite Element Method Fatigue, eine spezielle Software für Betriebsfestigkeitsberechnung auf derGrundlage von FEM Berechnungsergebnissen, die mit belieben anderen FE-Systemen ermittelt werden können.
1.1. Anwendung der FEM 5
Bild 1.1-8 Teil des Lagerstuhles eines Dreizylinder-Motors
aus den Ergebnisse mit entsprechenden Last-Zeitfunktionen für einen Zeitraum von400 Stunden Vollastbetrieb des Motors die Datensätze (Spannungs-Zeitverläufe) fürdie Betriebsfestigkeitsrechnungen generiert. Das Bild 1.1-9 zeigt als ein Ergebnisdieser Berechnungen die Verformung des Lagerstuhl als Farbflächenbild.Die Betriebsfestigkeitsrech-nungen geben schließlich einenÜberblick über die vorhande-nen Sicherheiten gegen Dauer-bruch. Deutlich lassen sichbereits in einem sehr frühenStadium des Entwurfes dieje-nigen Bereiche identifizieren,die schädigungsgefährdet sind.Die Computersimulationenführen in Verbindung mitexperimentellen Untersuchun-gen im Labor zu einer deut-lichen Verringerung der Ent-wicklungszeiten.
Das Verformungs- und Tragverhalten moderner Pkw- und LKW-Reifen wirdwesentlich durch den Gürtelverband und die Karkasse geprägt. Diese Reifenbauteilebestehen aus mehreren Lagen Gummi, in die Festigkeitsträger aus Stahl oder textilenMaterialien eingebettet sind. Für die statischen und dynamischen Berechnungen, diehochgradig nichtlinear sind (große Verformungen und nichtlineares Material-
Bild 1.1-9 Verformung des Lagerstuhl beim Zündlastfall 23 dargestellt als Farbflächenbild
6 1. Elementare Einführung in die FEM
verhalten), werden meist firmeninterne FEM-Softwareentwicklungen eingesetzt, diedie hausinterne Kompetenz der großen Reifenkonzerne auf dem Fachgebietverkörpern.
In seiner Diplomarbeit bei Continental hat Andreas Härtwig8 Untersuchungenzur Modellierung durchgeführt. So hat er hat dabei zum Beispiel den Einfluß derverwendeten Mischungsregeln (nach Chamis und Halpin-Tsai) und des Rebar-Konzeptes (das ist. die Berücksichtigung der Stahldrähte zur Bewehrung von Reifensin Form von diskreten finiten Elementen) auf die Ergebnisse untersucht. Ebenso wurdeder Einfluß der Vernetzung, des Lagenaufbaus und der Modellierung der Gürtelkanteund anderes untersucht und bewertet. Das Bild 1.1-10 zeigt den Querschnitt durcheinen der untersuchten PKW-Reifen, und im Bild 1.1-11 ist eines der für dieRechnungen benutzen FE-Modelle abgebildet.
Die Berechnung von ausreichend genauen Spannungen im Bereich von Kerbenist bei vielen Bauteilen entscheidend, weil diese Spannung meist für dieDimensionierung und Lebensdauerabschätzung des Bauteil relevant sind. Durchunzureichende Vernetzungen9 werden die Spannungen im Kerbengrund teilweiseextrem a) überschätzt, wenn die Vernetzung eigentlich ausreichend fein ist abereinspringende Ecken aufweist, die im Berechnungsmodell zu Spannungssingularitätenführen oder b) unterschätzt, wenn die Vernetzung zu grob ist. Besonders dieUnterschätzung der tatsächlichen Spannungen kann zu fatalen Folgen führen, wie zum
8 Härtwig, Andreas: Modellierung der mechanischen Eigenschaften des Gürtelpaketes von Reifen. DiplomarbeitUni Magdeburg, Institut für Mechanik, 1995. (Aufgabenstellung Continental AG Hannover, Dr. Fornefeld).9 Die Vernetzung hängt von der Qualität der verwendeten finiten Elemente und der Lösung (Regularität derLösung) ab, wobei das Problem im Bereich von Kerben der dort auftretende meist große Spannungsgradient ist.
Bild 1.1-10 Querschnitt durch einen Bild 1.1-11 FE-Netz des Reifens PKW Reifen
1.1. Anwendung der FEM 7
Bild 1.1-12 PKW Kurbelwelle
Beispiel der Katastrophen mit der norwegischen Erdölplattform, die durch einenderartigen Fehler zusam-mengebrochen und im Meer versunken ist.
Die Vielzahl vonFEM-Berechnungen anBauteilen mit Kerben, diebei der Entwicklung vonPKW’s durchgeführtwerden, haben VWbewogen, dieses Themaeinmal im Rahmen einerDiplomarbeit gründlichuntersuchen zu lassen.Dabei sollte auch der Ein-fluß, der sich aus der Ver-wendung unterschiedlicherBerechnungssoftware erge-ben könnte, in den Be-
trachung einbezogen wer-den. Als charakteristisches Bauteil wurde sich für eineKurbelwelle entschieden (siehe Bild 1.1-12), wobei als einfachere Referenzlösungauch noch eine biege- und torsionsbelastete Welle mit Absatz untersucht werden sollte(siehe Bild 1.1-15).
Diese Aufgabe wurde von Sven Müller10 im Rahmen seiner Diplomarbeiterfolgreich gelöst, wobei in die Untersuchungen die FEM-Systeme NASTRAN,COSMOS, MECHANICA und PROCISION/EAM sowie das BEM-System BEASYeinbezogen wurden.
Ein Segment der Kurbelwelle zeigt Bild 1.1-13, die berechneten Hauptspannungensind im Bild 1.1-14 dargestellt. Deutlich erkennbar sind die Spannungsextrema im
10 Müller, Sven: Bewertung unterschiedlicher FEM- und BEM-Techniken an Hand einer Kurbelwellenberech-nung. Diplomarbeit, Uni Magdeburg, Institut für Mechanik, 1997. (Aufgabenstellung Volkswagen AG Wolfs-burg, Dr. Stamerjohanns).
Bild 1.1-13 Ausschnitt aus derKurbelwelle Bild 1.1-14 Spannungen in Kurbelwelle
8 1. Elementare Einführung in die FEM
Kerbbereich. Ein entsprechendes Ergebnis für die einfache Welle ist im Bild 1.1-15zu sehen.
Mit den genannten Beispielen aus Diplomarbeiten konnte nur ein kleinerEindruck von den enormen Möglichkeiten der Finite-Element-Methode vermitteltwerden, wobei dabei auch bereits deutlich werden sollte, das die Anwendung nichtimmer ganz unproblematisch ist.
In 62 weiteren am Lehrstuhl für Numerische Mechanik seit 1992 betreutenDiplom- und Studienarbeiten wurden viele interessante Aufgaben und Probleme unterNutzung der FEM gelöst, wobei neben der reinen Anwendung kommerziellerProgrammsystem vielfach auch die Erweiterung und Ergänzung vorhandenerBerechnungssoftware durch Neu- und Weiterentwicklungen, die Gestaltung vonSchnittstellen zum Datentransfer, die Realisierung spezieller Auswerteprogramme undähnliches Bestandteile der Aufgabenstellungen waren. Die Anwendungen reichtendabei von der Auslegung eines Kometenlanders, Untersuchungen an Radiosatteliten,Flugzeugflügeln und PKW Teilen, die Berechnung von Hüftgelenkprothesen und vonKniegelenken, die Simulation von flexiblen von Endoskopen für die Gehirnchirurgie,die Simulation des Langzeitverhaltens von Endlagerstätten im Salzgestein, dieAuslegung von geregelten (smarten oder adaptiven) Strukturen, das Bruchverhaltenvon piezoelektrischen Materialien bis hin zu Methoden und Softwareentwicklungenfür die Nutzung von Parallelrechnern.
Einen Überblick über den Leistungsumfang unseres FEM-Systems COSARkönnen Sie der Internetseite http://www.femcos.de entnehmen.
In einer Vielzahl von Anwendungsgebieten wird die FEM heute erfolgreicheingesetzt. Einen nur unvollkommenden Eindruck von der Vielfalt dieserAnwendungen vermittelt die folgende Zusammenstellung:
Bild 1.1-15 Hauptspannungen in einer auf Biegung belasteten Welle
1.1. Anwendung der FEM 9
• Berechnung von Bauwerken aller Art (Brücker, Hochhäuser, Fundamente, Stau-dämme, Eisenbahntunnel, Kühltürme, Atomkraftwerke mit Simulation der Folgenvon Flugzeugabstürzen und Erdbeben, Auslegung von Brennstofftransportbehäl-tern, und vieles andere mehr.
• Auslegung von Hörschall-Wandlern in Hörgeräten, Simulation des Verhaltens vonKnochen und Gelenken aller Art, Simulationen in der Orthopädie, Kiefer- und Ge-sichtschirurgie, Neurochirurgie, Auslegung von Hüftgelenkimplantaten, Zahnspan-gen, intraokularen Linsen, optimalen Auslegung von Herzschrittmachern, Simula-tion elektrischer Felder im menschlichen Herzen, und viele weitere Anwendungenim Bereich der Medizin und Medizintechnik.
• Entwicklung des optimalen Designs von Turnschuhsohlen, Tennisschlägern, Surf-brettern, Skiern, Golfschlägern, Segel-, Ruder und Paddelbooten, Turngeräten,Bobschlitten und vielen weiteren Sportgeräten für den Spitzensport.
• Thermische Berechnungen an Bügeleisen, dynamische Untersuchungen anWaschmaschinen, Optimierung von Brillengestellen, Auslegung von Korkenzie-hern, Berechnungen von Waschbecken und viele weitere Anwendungen für Haus-haltsgeräte.
• In der Geophysik benutzt man die FEM zu Simulation des Langzeitverhaltes vonEndlagerstätten, untersucht das Auseinanderdriften der Kontinentalschollen, Er-mittelt die Ausbreitung von Erdbebenwellen und ähnliches mehr..
• Anwendung findet die FEM bei Untersuchungen zum Baumwachstums, zum Ver-halten von Getreide, Kartoffeln, Äpfeln,…,. in Speichern, und viele weitere An-wendungen in der Biologie.
10 1. Elementare Einführung in die FEM
1.2. Einführung in die Grundlagen der FEM
Grundidee der MethodeWährend sich bei den klassischen direkten Verfahren der Variationsrechnung, wiedem Verfahren von Ritz oder dem Galerkin-Verfahren, die Ansatzfunktionen über dasgesamte Gebiet (z.B. eine Struktur) erstrecken, werden bei der FEM Ansatzfunktionengewählt, die jeweils nur in einem Teilbereich, einem sogenannten finiten Element, un-gleich Null sind. In Bild 1.2-1 sind die einzelnen Schritte bei der Anwendung der FEMübersichtlich dargestellt.11
Bild 1.2-1 Schritte bei der Anwendung der FEM
Nach der Modellbildung wird zunächst die Geometrie der betrachteten Struktur in ein-zelne Teilgebiete (finite Elemente) zerlegt. Der in Bild 1.2-1 behandelte Kragbalkenwird hier beispielsweise als ebenes Scheibenproblem betrachtet (Modellbildung dis-kutieren). Man zerlegt die Mittelfläche des Kragbalkens in finite Elemente. Dies führtzum sogenannten FE-Netz. Zur Vernetzung des Problems stehen verschiedene finiteElementtypen zur Verfügung. Im vorliegenden Fall sind dazu Dreiecks- und Viereck- 11 Als eine simple Analogie soll nachfolgend die Berechnung der Kreisfläche π=
2rA dienen. Betrachtet man
einen Kreis mit dem Radius r = 1, so ist die exakte Lösung A = π. Ersetzt man den diesen Kreis durch ein n-Eck,läßt sich der Kreisinhalt lediglich mit Kenntnis der Formel für den Flächeninhalt eines Dreiecks approximieren.Durch immer feinere Unterteilung (Vernetzung) in Dreiecke kann der Flächeninhalt beliebig genau angenähertwerden.
1
n = 4 n = 8A = π
A = 2,828A = 2Weitere Unterteilungen liefern: n=64: A=3.136; n=128: A=3.1400; n=256: A=3.1412; n=1000: A=3.1415.
FApproximationder Geometriedurch finiteElemente
F
Element-typen
Elementsteifigkeitsmatrizen
Überlagerung der Matrizen
entsprechend der Element-
knotenzuordnung
Geometrische Randbedingungenund Belastungen
Lösung des Gleichungssystems
> Deformationen
ErmittlungderSpannungen
1.2. Einführung in die Grundlagen der FEM 11
selemente verwendet worden. Das Steifigkeitsverhalten der Gesamtstruktur wird durchdie Summe der Steifigkeitsmatrizen der einzelnen finiten Elemente beschrieben. Diesewerden additiv zur Gesamtsteifigkeitsmatrix überlagert. In Verbindung mit den geo-metrischen Randbedingungen und den Belastungen ergibt sich dann aus der Ge-samtsteifigkeitsmatrix ein lineares Gleichungssystem. Die zu berechnenden unbe-kannten Größen in diesem Gleichungssystem sind die Verformungen an den soge-nannten Knotenpunkten. Aus den Verformungen lassen sich anschließend in jedemeinzelnen finiten Element unter Verwendung der Ansatzfunktionen Spannungen be-rechnen.
Bild 1.2-2 Ein Auswahl finite Elemente des Programmsystems COSAR
Einige Standardelemente für die StrukturmechanikVon großer Bedeutung bei der Anwendung der FEM ist die Wahl der bei der Vernet-zung zu verwenden Elementtypen. Grundsätzlich unterscheidet man in der Struktur-mechanik folgende Modelle:
• Stabmodelle• Balkenmodelle• Scheibenmodelle (2D Modelle)• Plattenmodelle• Schalenmodelle (2.5D Modelle)• rotationssymmetrische und axialsymmetrische Modelle• Volumenmodelle (3D Modelle)• Spezialmodelle (z.B. für Rohleitungen, Stabschalen u.ä.)
Die verschiedenen Elemente, die für die Berechnung der oben genannten Strukturmo-delle erforderlich sind, unterscheiden sich in der Geometrie, der Anzahl und der Artder Freiheitsgrade und den Ansatzfunktionen. Das Bild 1.2-2 zeigt eine Auswahl vonfiniten Elementen, die im Programmsystem COSAR zur Verfügung stehen.
Mathematische FormulierungDie analytische Lösung einer Differentialgleichung, die ein Berechnungsproblem be-schreibt, ist meist nur für Sonderfälle möglich. Der beiden wesentlichen Weg für dienumerische Lösung eines elastostatischen Problems - die Nutzung der FDM oder dieNutzung der FEM - sind in Bild 1.2-3 dargestellt. Der linke Zweig der Darstellung
Hexaeder1 2 3
4567
8
9 10
1112
1314
15
16
17181920
uuu
1
2
3
Φ
12 11
10
98
7
6
54
32
1
Chisel
Rod
uuu
1
2
3
Φ
ξ1
1
2
Tetraeder
1
2 3 4
5
67 8
910
ξ1
ξ2
2 3
4
57
Quadrilateral
6
1
8uu
1
2
Φ
Anvil
12
34
56
78
910
11 1213
14
15161718
6
Wedge
1
2 3 4
57
89
10
111213
ξ1
ξ2
5
6
1 2 3
4
Triangle
Pentaeder
1
2
3
4
5
6 7
8
9
10
11
12
13
14
15
ξ3
ξ2
ξ1
Layered shell elements (semi-loof elements)
ξ3
ξ2
ξ1
uuu
1
2
3
Φ
ξ2
ξ1
uuu
1
2
3
Φ
ξ3 ξ2
2
ξ1
Layered elements
12 1. Elementare Einführung in die FEM
zeigt den Weg, der von der Lösung der Differentialgleichung11 ausgeht und über eineDifferenzenapproximation auf die Lösung eines Gleichungssystems führt. Bei der Fi-nite-Elemente-Methode wird meist vom Prinzip der virtuellen Verschiebungen (Ver-rückungen) ausgegangen. Die Einführung von bereichsweise definierten Näherungsan-sätzen für die unbekannten Verschiebungen führt auf die Methode der finiten Ele-mente, die zur Berechnung der Unbekannten ebenfalls die Lösung eines Gleichungssy-stems erfordert.
Elastostatisches Problem
Differentialgleichungen Extremalaussagen
Grundgl. der Elastizitätstheorie Prinzip der virtuellen Verrückungen
Gleichgew.-Bedingung
Verträglichkeitsbed.
Stoffgesetz
DGL-System
Numerische Lösung
Finite Differenzen
Stoffgesetz
- Gleichgewichtsbed. werden erfüllt- Verträglichkeitsbed. nicht a priori
erfüllt
Numerische Lösung
Finite-Elemente
Funktional
Bild 1.2-3. Die Näherungsverfahren FDM und FEM für elastostatische Probleme
In der Tabelle 1.2-1 sind die FEM und die FDM gegenübergestellt, wobei hier auchnoch die Randelementmethode (BEM) als weiteres wichtiges Verfahren aufgeführt ist.
Variationsformulierungen⇓
Unterteilung des Gebietes infinite Elemente
⇓Ansatzfunktionen in den Ele-
menten⇓
FEM
e
Differentialgleichungen⇓
Einführung eines Gitternetzesmit Gitterpunkten
⇓Differenzenquotienten an den
Gitterpunkte⇓
FDM
Integralgleichungen⇓
Unterteilung des Gebietsran-des in Randelement
⇓Ansatzfunktionen in den Ran-
delementen⇓
BEM
Randelement
Tabelle 1.2-1
11 Formuliert man das Problem in den Verschiebungen, so erhält man die Navier'schen Gleichungen, welche ausdrei partiellen Differentialgleichungen für die drei Verschiebungen im Raum bestehen. Führt man die Spannun-gen als Unbekannte ein, so erhält man die Beltrami-Gleichungen. Eine analytische Lösung der Differentialglei-chungen ist nur für Sonderfälle möglich.
1.2. Einführung in die Grundlagen der FEM 13
Nachfolgend soll dieser Sachverhalt für die FEM noch etwas genauer dargestellt wer-den. Dazu betrachten wir ein beliebiges 2D Gebiet (siehe Bild 1.2-4).
q
y,v
x,u
Bild 1.2-4 2D-Gebiet mit einem Netz aus finiten Dreieckselementen
In den Dreiecken (siehe Bild 1.2-4) soll der Verschiebungszustand durch einen linea-ren Ansatz angenähert werden. In jedem Element gilt also für die Verschiebungen
( )( ) yaxaay,x v
yaxaay,x u
654
321
++=++=
(1.2.-1)
Eine bessere Näherung für den Verschiebungszustand wäre zum Beispiel die Nutzungeines Polynoms höherer Ordnung, das für u(x,y) folgendermaßen geschrieben werdenkann
( ) [ ]
⋅=
n
2
1
22
a.
aa
... y xy y x x 1y,xu (1.2-2)
Die Größen a1 bis an sind die unbekannten Koeffizienten des Ansatzes. Zur Ermittlungdieser Unbekannten wird hier Prinzip vom Minimum des elastischen Potentials be-nutzt, das folgendermaßen formuliert werden kann:
Unter allen kinematisch zulässigen Verschiebungen nimmt das elastische Potential πfür die exakte Lösung einen minimalen Wert an.
Die Unbekannten werden also so bestimmt, daß das elastische Potential, das durch dieeingeführten Formfunktionen und Unbekannten ausgedrückt worden ist, einen mini-malen Wert annimmt. Für eine Stab (Elastizitätsmodul E, Fläche A, Länge l), derdurch eine Volumenlast p belastet ist, lautet das elastische Potential beispielsweise:
∫ ∫−
=π
l l
dxpuAdxxdud
AE21 2
(1.2-3)
Nach Einsetzen der Ansatzfunktionen und Integration wird ( )n21 a,...,a,a π minimiert,d.h. nach jeder einzelnen Unbekannten ai differenziert und das Ergebnis Null gesetzt(Extremwertaufgabe). Es ergibt sich das Gleichungssystem
14 1. Elementare Einführung in die FEM
0ai
=∂∂π
i=1,2,...n (1.2-4).
ZusammenfassungGegeben ist eine Differentialgleichung bzw. ein Differentialgleichungssystem für eineFunktion ( )z,y,xu . Die Finite-Elemente-Methode ist ein numerisches Verfahren, dasdiese Differentialgleichung, die eine kontinuierlichen Lösung besitzt, in ein algebrai-sches Gleichungssystem mit endlich vielen Unbekannten (diskrete Lösung) überführt.Die Differentialgleichung ist physikalisch betrachtet eine Gleichgewichtsaussage. An-stelle der Differentialgleichung wird eine duale Variationsaufgabe gelöst. Sie führt aufein algebraisches Gleichungssystem. In unserem Beispiel ist die Minimalforderung fürdie potentielle Energie des elastischen Systems äquivalent mit der Gleichgewichtsaus-sage der Differentialgleichung. Das Variationsprinzip ist eine schwache Form desGleichgewichts, d.h. bei Nutzung zulässiger Ansatzfunktionen wird das Gleichgewichtnur im integralen Mittel erfüllt, nicht aber lokal. Die gesuchte unbekannte Funktionwird durch bereichsweise definierte Ansatzfunktionen approximiert (Bild 1.2-5), diebestimmte Stetigkeits- und Kontinuitätsforderungen erfüllen müssen, sich aus der Va-riationsformulierung ergeben.1
Bild 1.2-5 Approximation einer Funktion in einem zweidimensionalem Gebiet
Zur Steigerung der Genauigkeit stehen mehrere Möglichkeiten zur Verfügung, z.B. dieh-Methode, die p-Methode oder die r-Methode. Unter dem Begriff h-Methode verstehtman die Steigerung der Anzahl der Elemente im Vernetzungsgebiet (h - ist ein Maß fürdie Elementgröße). Unter der p-Methode versteht man die Steigerung des Polynom-grades der Ansatzfunktionen im Element (p - ist ein Maß für den Polynomgrad). Dier-Methode schließlich hat die Verdichtung des Netzes bei gleichbleibender Anzahl vonElementen in Bereichen hoher Spannungsgradienten zum Ziel. Die einzelnen Metho-den können auch kombiniert werden. Entscheidend für die Auswahl einer Methode istderen Einfluß auf das Konvergenzverhalten der FE-Lösung. Durch die Verfügbarkeitvon Fehlerschätzern und unter Ausnutzung von a-priori Wissen läßt sich z.B. das FE-System so steuern, daß das Netz in Gebieten mit großen Fehlern automatisch verfei-nert wird. 1 Siehe dazu die Kapitel zur Variationsrechnung und zu den Direkten Verfahren der Variationsrechnung derVorlesung Mathematische und numerische Methoden der Mechanik – Teil I.
1.3. Kommerzielle FEM-Softwareprodukte 15
1.3. Kommerzielle FEM-Softwareprodukte
Struktur von FEM-SoftwareEin FEM-System besteht meist aus den folgenden drei Bestandteilen:
• Preprozessor (Dateneingabe)• FEM-Analyse (Berechnung)• Postprozessor (Ergebnisausgabe)
In der Regel gibt es für diese drei Teile separate kommerzielle Softwareprodukte, dieüber Datenschnittstellen miteinander kommunizieren. Üblicherweise kann der Postpro-zessor ebenfalls über Datenschnittstellen Geometriedaten aus CAD-Systemen über-nehmen.
Im Preprocessing (interaktiv) erfolgt im wesentlichen• die Eingabe der Geometriedaten des Modells (häufig über CAD-Schnittstellen),• die Vernetzung des Modells (Elementtypen, Anzahl der Elemente),• die Eingabe von Materialdaten (z.B. Dichte, E-Modul, Querkontraktionszahl),• die Beschreibung der Randbedingungen (z.B. Nullverschiebungen),• die Eingabe der Belastungen (z.B. Einzelkräfte, Linienlasten, Flächenlasten)• die Kontrolle der Eingabedaten.
Im FEM-Analyseteil erfolgt• die Berechnung der Elementsteifigkeitsmatrizen und Elementlastvektoren, • der Aufbau der Struktursteifigkeitsmatrix,• der Aufbau der Lastvektoren,• die Einarbeitung der Randbedingungen,• die Lösung des Gleichungssystems (Berechnung der Knotenverschiebungen),• die Rücksortierung der Elementverschiebungen,• die Berechnung der Dehnungen in den Elementen (falls gewünscht),• die Berechnung der Spannungen in den Elementen.
PostprocessingNach der Berechnung werden die Ergebnisse im Postprocessing ausgewertet. Dabeisteht der Nutzer in der Regel im Dialog mit dem Rechner und entscheidet über dieForm der Ausgabe bzw. über eine Verwertung der Ergebnisse für eine erneute Berech-nung (z.B. durch Netzmodifikation). Die Ausgabe der Ergebnisse kann sowohl gra-fisch als auch in Form von Listen erfolgen. Hierbei ist meist sowohl die Ausgabe vonErgebnisdaten als auch von Modelldaten möglich.
16 1. Elementare Einführung in die FEM
Eine Auswahl von FEM-Software:
• General-Purpose Programmsysteme:ABAQUS, ADINA, ANSYS, COSMOS, COSAR, MARC, MSC/NASTRAN, NISA, PERMAS,ALGOT, ANTRAS, ASAS, ASKA, BERSAFE, FEMAS, GFS-FEM, GIFTS, LUSAS, Mosaic, SAP80, STARDYNE, STRUDL, SYSTUS, UAI/NASTRAN
• Pre- und Postprocessoren (teilweise auch mit eigener FEM-Analyse)PDA/PATRAN, SDRC/SUPERTAB (CAEDS), FEMFAM, FEMAP
• CAD-Systeme mit integriertem FEM-Lösungsmodul (meist lineare Statik)I-DEAS, BRAVO, CADDS-STRESSLAB, EUCLID, INTERGRAPH
• SpezialprogrammeRohrleitungsberechnungen:
KEDRU, NUPIPE, PIPESTRESS, ROHR2Crashsimulationen, Metallumformung:
LS-DYNA3D, LARSTRAN, MSC/DYNA, PAM-CRASH, RADIOSSAkustik:
SYSNOISE, BEMAPMagnetfeld:
ANSYS, FLUX, MAGNET, MAXWELL, MSC/EMAS, PROFI, TOSCA
Strömung:FIDAP, FIRE, FLOTRAN, MSC/PISCES
Spritzgußsimulation:CADMOULD, MEFISTO, MOLDFLOW
1.4. Ein Einführungsbeispiel zur Statik von Stäbe 17
1.4. Ein Einführungsbeispiel zur Statik von Stäben1
Am Beispiel von Stabberechnungen wird nachfolgend die Methode der finiten Ele-mente erläutert. Die Differentialgleichungen der Stabstatik sind nachfolgend zur Erin-nerung aufgeführt.
Zusammenhang zwischen Dehnungen und Verschiebungen:
xu
ll
∂∂
=∆
=ε (1.4-1)
Gleichgewichtsbedingungen:
0px
=−
∂σ∂
(1.4-2)
Materialgesetz (Hookesches Geset)z:
ε=σ E . (1.4-3)
Aus diesen drei Gleichungen folgt die Naviersche Differentialgleichung des Zugstab-problems:
.0pxu
Ex
=+
∂∂
∂∂
(1.4-4)
Das vorliegende Problem ist dual zur Variationsformulierung2 in Form des Prinzipsvom Minimum des elastischen Potentials:
∑∫∫=
→−−ε=πn
1LLL
VV
T !MinFudVpudVEe21
, (1.4-5)
das durch die Verschiebungen (1.4-1) ausgedrückt lautet:
( ) ( )!MinFuxdpuAxd
xu
AE21 n
1LLL
l
2
l
→−−
∂∂
=π ∑∫∫=
(1.4-6)
Als finites Element wird ein Stabelement mit zwei Knoten gewählt (Bild 1.4-1).
l =2le
1 2x
x = - l x = + l
u , F u , F1 1 2 2
Bild 1.4-1 Finites Stabelement
Die Verschiebungen im Element werden durch folgende Ansatzfunktion approximiert3
1 Diese sehr ausführliche Darstellung dient zur Einführung in die FEM für Studenten mit geringeren Eingangs-kenntnisse und dient dem besseren Verständnis der allgemeinen Ableitung der Grundgleichungen der FEM inMatrizenschreibweise, wie sie im Abschnitt 2 erfolgt.2 Siehe Vorlesung Mathematische und Numerische Methoden der Mechanik – Teil I, Abschnitt Variationsrech-nung .
18 1. Elementare Einführung in die FEM
( ) x aaxu 21 += . (1.4-7)
In diesem Polynomansatz lassen sich die Unbekannten 1a und 2a durch mechanischdeutbare Größen - die Verschiebungen der Knotenpunkte 1 und 2 - ersetzen. Die Ver-schiebungen am linken und rechten Rand des Elements sind
( ) l l 211 aauxu −==−= (1.4-8a) ( ) l l 212 aauxu +==+= (1.4-8b)
Daraus folgt
( )
( )122
211
uul2
1a
uu21
a
−=
+=
(1.4-9)
Damit ergibt sich für den Verschiebungsansatz folgender Ausdruck:
2
)x(2
1
)x(1
u
N
)lx
1(21
u
N
)lx
1(21
)x(u4342143421
++−= (1.4-10)
N1(x) und N2(x) werden als Formfunktionen (engl. shape functions) bezeichnet. DieAbleitung der Verschiebung im Element ergibt
( ) ( )21e
212122
11 uu
1uu
l21
u21
u21
u x
Nu
x N
xu
+−=+−=+−=∂
∂+
∂∂
=∂∂
lll
(1.4-11)
Einsetzen des Verschiebungsansatzes in das Variationsproblem liefert für ein finitesElement
2211e
2
eee FuFudxpuAdxxu
EA21
−−−
∂∂
=π ∫∫−−
l
l
l
l
(1.4-12)
( ) ( )
( ) ( )
( ) ( ) 22112211e2221
21
e
ee
22112211e2221
21e2
eee
22112211e2221
212eee
FuFudxpuNuNAuuu2uEA
21
FuFudxpuNuNAuuu2u1
EA21
FuFudxpuNuNAdxuuu2u41
EA21
−−+−+−=
−−+−+−=
−−+−+−=π
∫
∫
∫ ∫
−
−
− −
l
l
l
l
l
l
l
l
l
ll
l
Aus den partiellen Ableitungen des elastischen Potentials bezüglich der freien Para-meter – d.h. den Knotenverschiebungen - ergibt sich folgendes Gleichungssystem:
3 Warum ist diese Ansatzfunktion ausreichend?
1.4. Ein Einführungsbeispiel zur Statik von Stäbe 19
( )
( ) 0FdxpNAuuEA
u
0FdxpNAuuEA
u
)(22e21
e
ee
2
e
)(11e21
e
ee
1
e
e
e
=−−+−=∂
π∂
=−−−=∂
π∂
∫
∫
l
l
l
l
(1.4-13)
In Matrizenschreibweise hat das Gleichungssystem folgende Form:
dxpNN
AFF
uu
1111EA
)( 2
1e
2
1
2
1
e
ee
e
l l
∫
+
=
−
− (1.4-14)
Mit den Abkürzungen Ke - Elementsteifigkeitsmatrix
−
−α=
−
−=
1111
1111EA
ee
eee l
K (1.4-15)
und fe - Elementkraftvektor)e(
p)e(
e FFf += (1.4-16)
mit
= )e(
2
)e(1)e(
FF
F )e(p
)e(e FFf += (1.4-17)
für die Knoteneinzelkräfte und mit
∫
=
e
dxpNN
A2
1l
)e(p
l
F (1.4-18)
für die Volumenlasten kann man Gl. (1.4-14) schreiben als
eee fuK =⋅ . (1.4-19)
Dieser Ausdruck wird als Elementsteifigkeitsbeziehung bezeichnet.
Falls beispielsweise der Kraftvektor infolge Schwerkraft in x – Richtung ermitteltwerden soll, ergibt sich mit
xgp ⋅ρ= , (1.4-20)
aus Gl. (1.4-18)
ρ=
ρ=
+
−⋅⋅ρ= ∫
−
11
gA21
22
21
gAdxx1
x1
21
gA
xee
xexep
l
ll
l
l l
l
F
(1.4-21)
In diesem Fall teilt sich also die Gesamtkraft je zur Hälfte auf die Knoten auf.
20 1. Elementare Einführung in die FEM
Wir untersuchen jetzt das Zusammenwirken der Einzelelemente im Gesamtsystem.Das gesamte elastische Potential berechnet sich aus der Summe der n Elementbeiträgezu
∑=
π=πn
1eegesamt (1.4-22)
Die Minimierung des elastischen Gesamtpotentials liefert
∑= ∂
π∂=
∂
π∂ n
1e i
e
i
gesamt
u
u
(1.4-23)
Setzt man hier die Elementbeiträge entsprechend Gleichung 1.4-13 ein und nimmt an,daß eine aufsteigende Nummerierung der Knotenpunkte vorliegt, dann ergibt sich fürn=3 das Gleichungssystem (1.4-27).4 Das gleiche Ergebnis läßt sich auch über dasKräftegleichgewicht an jedem Knotenpunkt herleiten (siehe Bild 1.4-2).
k-1 k k+1e e+1
f k - äußere Kraft am Punkt k
e e+1
f k
f f f f1
e2e
1 2e+1 e+1
Bild 1.4-2 Kraftgleichgewicht am Knoten k
Das Kraftgleichgewicht am Knoten k liefert:
( ) ( )1k1ek1ee1ke
1kk1ek1ke
1e1
e2k
uu)(u
uuuufff
+++−
++−
+
α−α+α+α−=
−α++−α=+=
(1.4-24)
In Matrizenschreibweise lautet diese Gleichung
[ ] k
1k
k
1k
1e1eee fuu
u =
α−α+αα−
+
−
++ (1.4-25)
Dies führt für das Beispiel mit n=3 Stäben, die durch die drei Elementsteifigkeitsma-trizen
=
⋅
−
−⋅α
2
1
2
11 f
fuu
1111
, (1.4-26a)
4 Beweisen Sie diese Behauptung.
1.4. Ein Einführungsbeispiel zur Statik von Stäbe 21
=
⋅
−
−⋅α
3
2
3
22 f
fuu
1111
, (1.4-26b)
=
⋅
−
−⋅α
4
3
4
33 f
fuu
1111
, (1.4-26c)
beschrieben werden, auf das Gesamtgleichungssystem
=
⋅
αα−α−α+αα−
α−α+αα−α−α
4
3
2
1
4
3
2
1
33
3322
2211
11
ffff
uuuu
000
000
. (1.4-27)
Die Knotennumerierung hat einen entscheidenden Einfluß auf die Struktur der Matrix.Das Ziel muß es sein, eine minimale Bandweite der Matrix zu erhalten. Die in Bild1.4-3 gezeigte nicht fortlaufende Knotennumerierung führt stattdessen auf die Gl.(1.4-28).
1 2 34
f f f f1 2 3 41 2 3
Bild 1.4-3 Beispiel einer nicht fortlaufenden Knotennumerierung
=
⋅
α+αα−α−αα−
α−α−α+αα−α
4
3
2
1
4
3
2
1
2121
33
2332
11
ffff
uuuu
000
000
. (1.4-28)
Es ist zu erkennen, daß die Bandstruktur der Matrix hier verlorengegangen ist.
In Bild 1.4-4 ist ibw das Maß für die halbe Bandweite einer symmetrischen Matrix.Diese kann aus folgender Formel bestimmt werden:5
maxedof dnibw ⋅= , (1.4-29)
wobei max ed die maximale Knotennummerndifferenz an einem Element ist, und ndof istdie Anzahl der Unbekannten (Freiheitsgrade) pro Knoten. Die "Rechenzeit" bei An-wendung von direkten Verfahren zur Gleichungslösung läßt sich näherungsweise nachfolgender Formel bestimmen
( )2cpu ibwnct ⋅⋅≅ . (1.4-30)
Dabei bezeichnet hier n die Gesamtzahl der Gleichungen, und c ist eine vom Prozessorder Rechenanlage abhängige Konstante. Aus Gleichung (1.4-30) wird der Einfluß der
5 Zeigen Sie die Richtigkeit dieser Formel.
22 1. Elementare Einführung in die FEM
Bandweite auf die Rechenzeit ersichtlich. Es ist deshalb bei der Anwendung direkterVerfahren zur Gleichungslösung6 unbedingt erforderlich, bei der Vernetzung auf dieEntstehung einer möglichst kleinen Knotennummerndifferenz zu achten, um die Re-chenzeit gering zu halten7.
Einarbeitung der RandbedingungenEin weiterer wichtiger Punkt bei der FEM-Analyse ist die Einarbeitung der geometri-schen Randbedingungen. Ohne die Berücksichtigung von Verschiebungsrandbedin-gungen ist das Gleichungssystem nicht eindeutig lösbar – es ist singulär. Ein Bauteilmuß also so gelagert sein, dass mindestens die Starrkörperbewegungen unterdrücktwerden!Im folgenden wird die Vorgehensweise zur Einarbeitung verschiedenen Randbedin-gungen aufgezeigt.
Verschiebung der Größe 0u i = am Knoten i ist vorgeschrieben (Nullverschiebung)a) Die i-te Zeile und Spalte werden aus dem Gleichungssystem gestrichen.b) Setzen einer großen Zahlen p
ii10k auf die i-te Hauptdiagonale:8
u
10k
ip
ii
=
⋅
MMMM
MM
M
OO
O
(1.4-31)
Verschiebung der Größe uu i = ist vorgeschriebena) Die i-te Spalte wird mit u multipliziert und von der rechten Seite subtrahieren,
dann werden die i-te Zeile und Spalte gestrichen:
6 Siehe Vorlesung Computer-Numerik.7 Dazu stehen in FEM-Programmen meist Verfahren zur Verfügung, die intern eine automatische Umnumerie-rung mit dem Ziel vornehmen, eine vorhandene Bandweite deutlich zu verringern. Ein einfaches Verfahren, dasgarantiert das Minimumum liefert, gibt es leider nicht (es sei denn man prüft alle möglichen n! Permutationeneinzeln durch).8 Überlegen Sie, wie groß p gewählt werden muß, um eine bestimmte Genauigkeit zu erreichen. Welche Proble-me sind bei dieser Methode zu erwarten?
symm.
0
n
ibw
Bild 1.4-4 Definition der Bandweite ibw
1.4. Ein Einführungsbeispiel zur Statik von Stäbe 23
−
=
ni
i2
i1
k
kk
u M
MMMM
MMMM
OLO
OLO
(1.4-32)
b) Setzen eines großen Zahlenwertes pii10k auf die i-te Hauptdiagonale sowie
Ersetzen der rechten Seite durch piii 10kuf = :
=
MM
M
MM
M
OLO
LOp
iiip
ii 10kuu10k (1.4-33)
Sind die Randbedingungen eingearbeitet worden, kann das Gleichungssystem gelöstwerden. Dazu stehen mehrere Möglichkeiten von direkte oder iterative Verfahren zurVerfügung, wie z.B. das Cholesky-Verfahren, das Gauß-Seidel Verfahren, CG-Verfahren usw. (siehe dazu die Vorlesung Computer-Numerik).
Nach der Lösung des Gleichungssystems sind die Verschiebungen an den Knoten be-kannt, und es lassen sich elementweise die Spannungen berechnen. Im folgenden istdie Ableitung der Spannungen aus den Verschiebungen für ein eindimensionales Sta-belement dargestellt (vergleiche Bild 1.4-1). Aus
xu
EE∂∂
⋅=ε⋅=σ (1.4-31)
folgt mit den Ableitungen der Ansatzfunktionen Gl. (1.4-11) für die Spannung imElement
( )21e
e uuE
+−=σl
, (1.4-32)
die hier elementweise konstant sind.
Ein BeispielAbschließend soll noch ein vollständiges Beispiel vorgerechnet werden (siehe Bild1.4-7). Für die Elementsteifigkeitsmatrix des e-ten Elementes gilt
( ) 1111
ee
−
−⋅α=K ,
e
eee l
EAmit
⋅=α . Mit den Eingabegrößen erhält man für αe:
mm/N100,2,mm/N100,2,mm/N106,1 43
42
41 ⋅=α⋅=α⋅=α .
24 1. Elementare Einführung in die FEM
Für die Gesamtsteifigkeitsbeziehung ergibt sich dann
−
=
⋅
−+−+
−
⋅
00
26000
uuuu
0,2.symm0,20,20,2
00,20,26,1006,16,1
10
4
3
2
1
4 .
Die Lagerbedingungen lauten 0u,0u 41 == . Das Streichen der Zeile/Spalte 1, 4 liefert
−
=
−
⋅=
⋅
−
− −
26,00
26000
10uu
0,40,20,26,3 4
3
2
Für die Verschiebungen ergibt sich daraus
u2 = - 0,1 mm,u3 = - 0,05 mm,
Für die Elementspannungen erhält man schließlich
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( ) 25
433
33
25
322
22
25
211
11
mm/N100005,0100
100,2uu
lE
mm/N5005,01,0200
100,2uu
lE
mm/N401,00500
100,2uu
lE
=++⋅⋅
=+−⋅=σ
=−+⋅⋅
=+−⋅=σ
−=−−⋅⋅
=+−⋅=σ
Beliebige Lage eines Elementes in der EbeneDas finite 2-Knotenelement ist natürlich auch für die Berechnung von ebenen undräumlichen Fachwerksystemen geeignet. Dazu braucht das Element lediglich auf eine
x
F
1
4
3
2
1
2
3
l
l
l
2
1
3
Bild 1.4-7 Stabbeispiel
N 2600F
mm 100l mm 10A
mm 200l mm 20A
mm 500l mm 40A
3,0mm/N 100,2E
12
1
12
2
12
1
25
=
==
==
==
=µ⋅=
1.4. Ein Einführungsbeispiel zur Statik von Stäbe 25
beliebige Lage im Raum transformiert zu werden. Für eine Transformation in der Ebe-ne (vergleiche Bild 1.4-8) ergibt sich:
Bild 1.4-8 Transformation in der Ebene
Die Transformation vom globalen in das lokalen Koordinatensystem lautet
α⋅+α⋅=α⋅+α⋅=
sinvcosuu
sinvcosuu
222
111 (1.4-33)
In Matrizenschreibweise lautet dieses Gleichungssystemer
.
vuvu
sincos0000sincos
uu
2
2
1
1
2
1 Tvv =
⋅
αα
αα=
= (1.4-34)
Durch die T-Matrix kann die lokale Steifigkeitsbeziehung auf das globale Koordina-tensystem transformiert werden. Man erhält
f=vTK
f=vK
⋅⋅
⋅(1.4-35)
Multiplikation von links mit TT liefert eine symmetrische Steifigkeitsmatrix
ffTvK
TKT =TT ⋅=⋅⋅⋅ 43421 , (1.4-37)
mit
αα⋅αα
α−α⋅α−αα⋅α−α−α⋅αα
2
2
22
22
sin.symmcossincos
sincossinsincossincoscossincos
= K
Mit dieser transformierten Steifigkeitsmatrix lassen sich beliebige ebene Fachwerksy-steme berechnen.
26 1. Elementare Einführung in die FEM
Ein Beispiel
1 3
2
F
l l
1 2x
yα
−αα = 45°
gegeben:
l , F, EA
Bild 1.4-8 Ebenes Fachwerksystem
Für das Beispiel in Bild 1.4-8 ergeben sich mit cos 45°= sin 45° = 221 ⋅ für die bei-
den Elemente 1 und 2 die folgenden Elementsteifigkeitsbeziehungen:
⋅
−−−−
⋅=
2
2
1
1
11
vuvu
1.symm111111111
2AEl
vK
⋅
−−
−−
⋅=
3
3
2
2
2
vuvu
1.symm11111
1111
2AEl2vK
Nach Einarbeitung der Randbedingungen folgt daraus das Gleichungssystem
F
0vu
1001
lAE
2
2
−
=
⋅
⋅
⋅,
dessen Lösung
EAF
v,0u 22l
−==
lautet. Zur Spannungsermittlung muß eine Rücktransformation auf das lokale Stabsy-stem erfolgen.
( )
( ) 0uEF
2v2u2u
11
21
221
2211
2
=⋅
−=⋅+=l
und man erhält
( ) ( )( )AF
221
EAF
22E
uuE 1
11
2 −=−=−=σl
ll
1.5 Einführungsbeispiel - Dynamik von Stäben 27
1.5. Ein Einführungsbeispiel zur Dynamik von StäbenIn der Elastodynamik verwenden wir als Variationsprinzip das Hamiltonsche Prinzip
∫ ⇒1
2
t
t
Ldt är station (1.5-1)
mitπ−= TL (1.5-2)
und
∫∫ ρ=
ρ=
)x(
2
V
2
xduA21
Vdtdud
21
T & (1.5-3)
Die Knotenverschiebungen, die in die Approximation des Verschiebungszustandeseingehen, vergleiche Gleichung (1.4-10), sind jetzt zeitlich veränderlich, und es ergibtsich
( ) ( ) ( ) ( ) ( )tuxNtuxNt,xu 2211 += (1.5-4)
2211 uNuNtdud
u &&& +==
(1.5-5)
Als Alternative wäre auch die Verwendung von Ansatzfunktionen im Raum-Zeit Be-reich denkbar, d.h. NI(x,t), was im Fall des hier betrachteten geometrisch eindimensio-nalem Stabelementes zu einem „ebenen“ Raum-Zeit-Element für die Approximationder Verschiebungen in Raum und Zeit führen würde.1 Diese Raum-Zeit-Diskretisierung führt aber auf eine extreme Vergrößerung der Anzahl der Unbekanntenund wird daher nur selten genutzt. Wir gehen deshalb nachfolgend von dem Nähe-rungsansatz (1.5-4) aus und setzen diesen in den Ausdruck für die kinetische Energie(1.5-3) ein. Es ergibt sich
( )∫ ++ρ=e
xdu)x(Nuu)x(N)x(N2u)x(NA21
T 22
222121
21
21eee
l
&&&& (1.5-6)
Die Auswertung der Integrale in (1.5-6) liefert
l llll
ll
l
l
l
l
l
l
l
l
l
32
32
241
x311
x1
x41
xdxx
2141
xdx
121
xdN
32
2
2221
=
+=
+−=
+−=
−=
−
−−−∫∫∫
lll l
l
l31
32
241
xdx
141
xdNN2l
l21 =
−=
−= ∫∫
−−
1 Der Verschiebungsansatz für ein solches Element mit linearem Ansatz im Zeitbereich (Elementlänge in derRaumdimension ist 2a und in der Zeitdimension 2∆t) lautet:
( ) ( ) ( ) ( ) ( ) 44332211 ut,xNut,xNut,xNut,xNt,xu ⋅+⋅+⋅+⋅=mit
∆−⋅
±=
tt
121
ax
121
N 2/1 ,
∆+⋅
±=
tt
121
ax
121
N 4/3 .
28 1. Elementare Einführung in die FEM
l l
l32
xdN 22∫
−
=
Damit ergibt sich für die kinetische Energie eines Elementes
++ρ= 2
22121eee u
32
uu32
u32
A21
T &&&& lll
( )2221
21ee uuuuA
31 &&&& ++ρ= l (1.5-7)
Die Bewegungsdifferentialgleichungen ergeben sich aus den Lagrangeschen Glei-chungen
0
uL
uL
t i
e
i
e =∂∂
−
∂∂
∂∂
& (1.5-8)
Mit( ) ( )
( ) eTeee
Teie
ieiee
21
uT
uuTL
fvvKv +−=
π−=
&
& (1.5-9)
und
( ) ( )
( ) ( )21ee21ee2
e
21ee21ee1
e
u2uA31
u2ut
A31
uT
t
uu2A31
uu2t
A31
uT
t
&&&&&&&
&&&&&&&
+ρ=+∂∂
ρ=
∂∂
∂∂
+ρ=+∂∂
ρ=
∂∂
∂∂
l
l
l
l
,
wobei 2 l = le ist, folgt aus (1.5-8)
eeeee fvKvM =⋅+⋅ && (1.5-10)
Dabei ist
ρ=
2112
A61
eeee lM (1.5-11)
die Massenmatrix des 2-Knoten-Stabelementes, und es gilt weiterhin
=
=
2
1e
2
1e u
uund
uu
vv&&&&
&& .
Im Unterschied zur Statik ist bei dynamischen Problemen ein Differentialgleichungs-system zu lösen. Das Gesamtsystem baut sich analog zur Statik auf (Einspeichern derElementmassen- und -steifigkeitsmatrizen in das Gesamtsystem), und es ergibt sich
fKvvM =+&& (1.5-12)
Falls noch eine geschwindigkeitsproportionale Dämpfung berücksichtigt wird, ergibtsich
fKvvCvM =++ &&& (1.5-13)
1.5 Einführungsbeispiel - Dynamik von Stäben 29
Für die Lösung dieses gewöhnlichen Differentialgleichungssystems (1.5-13/14) wer-den bevorzugt Zeitintegrationsverfahren eingesetzt (z.B. Zentrale Differenzen Metho-de, Newmark Verfahren oder Wilson Verfahren, siehe dazu die Vorlesung Computer-Numerik). Auch Lösungen im Modalraum oder im Frequenzbereich sind möglicheAlternativen.Von großer praktischer Bedeutung ist die Berechnung der freien, ungedämpftenSchwingungen, d.h. die Ermittlung der Eigenfrequenzen und –formen, die sich aus derLösung von (1.5-12) mit 0=f ergeben. Mit dem Separationsansatz ( ) tcost vv 0 ω=
und der 2. Ableitung dieses Ansatzes ( ) tcost 2 vv 0 ωω−=&& ergibt sich aus (1.5-12) dasfolgende Matrizeneigenwertproblem zur Bestimmung der Eigenfrequenzen ω
( ) 002 =⋅⋅ω− vMK (1.5-14)
Ein BeispielAls Beispiel wollen wir die Schwingungen eines einseitig gefesselten Stabes (sieheBild 1.5-1) betrachten.
x,u
lsBild 1.5-1 Freie Längsschwingung eines Stabes
nchwingungeStablängssenzen der Eigenfrequ:.ges
E,l,,A:.geg s
ρ
Die exakte Lösung ergibt sich aus der Lösung der Dgl. uAuEA &&ρ=′′ .2
2 Mit dem Separationsansatz ( ) ( ) ( )tTxUt,xu = sowie TUu &&&& = und TUu ′′=′′ liefert dieser Ansatz
TUATUEA &&ρ=′′ , und daraus folgt 2
TT
UUE
ω−==′′
⋅ρ
&&. Daraus ergeben sich die beiden gewöhnlichen Dgln.
0TT 2 =ω+&& und 0UE
U 2 =ρ
ω+′′ mit den mit den Lösungen tsinBtcosAT ω+ω= und
( ) ωρ
+ωρ
=E
sinDxE
cosCxU . Aus der zweiten Gleichungen lassen sich die Eigenwerte berechnen. Die
Randbedingungen lauten
a) ( ) 0t,0xu == , woraus 0C = folgt und
b) ( ) 0t,xF sN == l , woraus die Eigenwertgleichung 0E
cosE
EADxdud
EA
0
s
s
=ωρ
ωρ
=43421
l
l
folgt. Die Lösung von 0E
cos s =ωρ
l liefert die Eigenwerte .....;23
;2E s l ππ
=ωρ
, Damit ergibt sich für die
ersten beiden Eigenfrequenzen
ρ=
ρπ
=ωρ
=ρ
π=ω
E7124.4E23
,E5708.1E
2 ss2
ss1 llll
.
30 1. Elementare Einführung in die FEM
Nachfolgend werden wir die Eigenfrequenzen mittels der FEM bestimmt. Zunächstwird nur ein finites Element benutzt (Bild 1.5-3).
1 2
Bild 1.5-2 Approximation des Stabes mit einem finiten Element
Nach dem Aufstellen der Elementsteifigkeits- und Massenmatrix, sowie der Einarbei-tung der Randbedingungen ergibt sich folgende Gleichung:
0u2A61EA
2s2
s
=
ρω− l
l (1.5-15)
mit der Lösung
ρ=
ρ=ω
E73.1E3
ss1 ll
. (1.5-16)
Der prozentuale Fehler ergibt sich aus
%10%100exakt 1
exakt 1FEM 1 =⋅ω
ω−ω=ε .
Die Lösung mit zwei finiten Elementen (siehe Bild 1.5-3) ergibt:
1 2 31 2
l ls s0,5 0,5
Bild 1.5-3 Stab mit zwei finiten Elementen
0uuu
210141012
21
A61
110121
011
21EA
3
2
1
s2
s
=
ρω−
−−−
−
ll
(1.5-17)
Um die Randbedingung 0u1 = zu berücksichtigen, werden die 1. Zeile und die 1.Spalte gestrichen. Daraus resultiert folgendes Eigenwertproblem:
0uu
2114
E41
61
1112
3
22s
2 =
ρω−
−
−l .
Mit der Abkürzung E24
1 2s
ρ=β l läßt sich das Eigenwertproblem folgendermaßen
schreiben:
1.5 Einführungsbeispiel - Dynamik von Stäben 31
0uu
2114
1112
3
22 =
⋅
⋅β⋅ω−
−
− .
Die Eigenwerte ergeben sich aus
0211
14222
22
=ω⋅β⋅−ω⋅β−−ω⋅β−−ω⋅β⋅−
.
Die charakteristische Gleichung lautet:
. 0ß1
71
ß1
710
224 =+ω−ω
Mit 24 a=ω erhält man daraus
ß7
235a 2/1
±= .
Die Eigenfrequenz ergibt sich aus
2/11/2 a =ω .
Die ersten beiden Eigenwerte ergeben sich zu
ρ=ω
E6114.1
s1 l
, der Fehler beträgt % 6,2 1 =ε .
l ρ
=ωE6293.5
s2 , der Fehler beträgt % 4,19 2 =ε .
Die Eigenformen für die Vernetzung mit zwei Elementen lassen sich folgendermaßenberechnen. Aus dem Eigenwertproblem
0uu
2114
1112
3
22 =
βω−
−
−
folgt mit 1u 3 = die Lösung
( ) 01u41u2 22
2 =+βω−−
und daraus
2
2
2 421
uβω−
βω+=
32 1. Elementare Einführung in die FEM
Damit ergibt sich für die 1. Eigenform mit l ρ
=ωE59666.2
2s
1 und
1082.0E
l59666.2
E241
2s
2s
21 =
ρρ
=βω l die Lösung
7071.0)1082.0(42
1082.01u 2 =
−+
=
Analog ergibt sich für die 2. Eigenform mit l ρ
=ωE689.31
2s
2 und 32.122 =βω
die Lösung
707.0)32.1(42
32.11u 2 −=
−+
=
Anmerkungen zur Lösung von Eigenwertproblemen• Die Eigenwerte nähern sich stets von "oben " der exakten Lösung an,
d.h. exaktFEM ω>ω .• Die höheren Eigenwerte haben eine geringere Genauigkeit, d.h. die Genauigkeit
nimmt mit steigendem Eigenwert ab.• Die Anzahl der Eigenwerte ist gleich der Anzahl der Unbekannten. Die höheren
Eigenwerte werden durch die FE-Unterteilung bestimmt und sind daher nichtsinnvoll. Die exakte Lösung hat unendlich viele Eigenwerte.
• Die Eigenform (Eigenschwingform) kann nur bis auf einen konstanten Faktorbestimmt werden, d.h. dieser Wert gibt nur das Verhältniswert zwischen deneinzelnen Amplituden an. Für die Darstellung der Eigenschwingform wird inder Regel auf die maximale Amplitude normiert.
1.6. Ein Einführungsbeispiel zur eindimensionalen Wärmeleitung 33
1.6. Ein Einführungsbeispiel zur eindimensionalen Wärmeleitung
Die Lösung der eindimensionalen Wärmeleitungsgleichung wird nachfolgend als Bei-spiel für ein anderes Feldproblem – ein Temperaturfeld – betrachtet. Gleichzeitig wirdan diesem Beispiel die numerische Lösung eines zeitabhängigen Vorganges durch eineinfaches Zeitintegrationsverfahren (zentrale Differenzenmethode) demonstriert.Aus der Fourierschen Gleichung folgt für den eindimensionalen Fall die Differential-gleichung der Wärmeleitung
tT
cQxT2
2
∂∂
ρ=+∂∂
λ in V (1.6-1)
Die thermischen Randbedingungen lauten
TT = auf OT (1.6-2a)
0q xT
=+∂∂
⋅λ auf Oq (1.6-2b)
( ) 0TT xT
=−⋅α+∂∂
⋅λ ∞ auf Oα (1.6-2c)
Es bedeuten:
Als finites Element wird wieder ein eindimensionales Element mit zwei Knoten ver-wendet (siehe Bild 1.6-1).
l l
x
l e Bild 1.6-1 Finites Element
Alle Ableitung erfolgen analog zum Stabelement (siehe Abschnitt 1.4). Der lineareAnsatz für die Temperaturverteilung im Element lautet
in Ktemperatur-UmgebungsT Kmm
Win rgangszahl- Wärmeübe
mmWdichte in Wärmestromq
mmWit in ergiebigkeWärmequellQ -
kg KJazität in . Wärmekapc - spezif
mmkg in - Dichte
KmmWs
in c
mm KWn ähigkeit iWärmeleitf-
in KTemperatur-T
2
2
3
3
3
∞
α
−
ρ
µ=ρ
λ
34 1. Elementare Einführung in die FEM
( ) ( ) ( ) 2211 TxNTxNxT += , (1.6-3)
mit den Formfunktionen
+=
−=
lx
121
N
lx
121
N
2
1
(1.6-4)
Dieser Näherungsansatz entspricht einer elementweisen Approximation der Tempera-turverteilung durch einen linearen Ansatz (d.h., ( ) xaaxT 21 += ).Für die Ableitung der FEM-Gleichungen wird als Ausgangsbasis keine Variationsfor-mulierung benutzt, sondern die "Methode der gewichteten Residuen" verwendet (sieheVorlesungsskript Mathematische und Numerische Methoden der Mechanik - Teil I).
0WdVRV
.Dgl =∫ , (1.6-5)
wobei RDgl. das Residuum der Differentialgleichung der Wärmeleitung (1.6-1) be-zeichnet. Wenn man für die Gewichtsfunktionen W jeweils die Ansatzfunktionen ein-setzt, so entspricht dies dem Galerkinschen Verfahren.
0xdNtT
cQx
TA
0xdNtT
cQx
TA
2)x(
2
2
1)x(
2
2
=
∂∂
ρ−+∂∂
λ
=
∂∂
ρ−+∂∂
λ
∫
∫
(1.6-6)
Mit A wird die Elementfläche bezeichnet, und es wird für T(x) der Ansatz (1.6-3) ein-gesetzt. Vom Zusammenhang zwischen dem Ritzschen Verfahren (bzw. dem Galer-kinschen Verfahren) und der FEM wissen wir, daß die Ansatzfunktionen bis zur (n-1)-ten Ableitung stetig sein müssen (n bezeichnet die höchste unter dem Integral vor-kommende Ableitung), damit die Konvergenz gesichert ist. Bei stückweise (element-weise) definierten Ansatzfunktionen in der FEM bedeutet dies, daß die Bedingungauch zwischen benachbarten finiten Elementen erfüllt sein muß. Die Anwendung derGleichung (1.6-6) bedingt, daß die Ansatzfunktionen stetige erste Ableitungen habenmüssen und daher der Ansatz (1.6-3) eigentlich nicht geeignet ist. Durch partielle Inte-gration läßt sich die Ordnung der Ableitung jedoch reduzieren. Es ergibt sich
xdx
NxT
NxT
xdNx
T
)x(
I
RandII
)x(2
2
∫∫ ∂∂
∂∂
λ−∂∂
λ=
∂∂
λ (1.6-7)
Der Randterm kann durch die Randbedingungen ausgedrückt werden und man erhält:
( )α
∞−⋅α−−=∂∂
λOIOI
RandI NTTqNN
xT
q
(1.6-8)
Damit lauten die Galerkinschen Gleichungen
1.6. Ein Einführungsbeispiel zur eindimensionalen Wärmeleitung 35
( )
0xdN t TcxdQN
NTTqNxdx
NxT
I)x()x(
I
OIOI)x(
Iq
=∂∂ρ−+
+−α−−∂
∂
∂∂
λ−
∫∫
∫ α∞
(1.6-9)
mit I=1,2. Hier gehen nur noch die ersten Ableitungen des Temperaturfeldes ein, sodaß der lineare Ansatz (1.6-4) ausreichend ist. Die Ableitung des Ansatzes liefert
2122
11 T
21
T21
Tx
NT
xN
xT
ll
+−=∂
∂+
∂∂
=∂∂
. (1.6-10)
Somit ergibt sich für I =1 die Gleichung
( ) 0xdNTNTNcxdQNNT
NT21
T21
qNxd21
T21
T21
122111O1
O121O121
q
=+ρ−+α+
+−α−−
−
+−λ−
∫∫
∫
−−∞
−
α
α
ll
llll
l
l
l
l
l
&&
Eine analoge Gleichung ergibt sich auch für I=2. Auswerten der Integrale und Zu-sammenfassen liefert
( ) fTCTKK =++ α& (1.6-11)
mit der Wärmeleitfähigkeitsmatrix
−
−λ=
1111
elK (1.6-12)
und der Wärmeübergangsmatrix
α
α=α
O2
1
N00N
K , (1.6-13)
die jedoch nur für solche Knoten auftritt, die zu einem Rand mit der Randbedingung(1.6-2c) gehören. Mit C wird die Wärmekapazitätsmatrix
ρ=
2112
c61
elC (1.6-14)
bezeichnet. Als Vektor der rechten Seite ergibt sich
+
α+
−=
α
∞ 11
Q21
NN
TNN
q eO2
1
O2
1
q
lf , (1.6-15)
wobei hier Q im Element als konstant angenommen wurde.Nachfolgend wird die Anwendung an einem einfachen Beispiel, der instationärenWärmeleitung in einer Wand (Bild 1.6-2), demo nstriert.
36 1. Elementare Einführung in die FEM
Instationäre Wärmeleitung in einer Wand
A
l = 1m
Wand
Bild 1.6-2 Wärmeleitung in einer Wand
Die Randbedingungen seien
( )( ) C00t,m5.0xT
C1000t,xTT0
°=≥±=°=<=
Das entspricht der Situation, daß eine auf 100 °C erwärmte Wand plötzlich auf 0 °Cabgekühlt wird. Es soll die zeitliche Änderung der Temperatur am Punkt A ermitteltwerden. Die exakte Lösung lautet (siehe Lehrbuch von Tautz):
+
π
−
π
π=
π⋅
−
π
−L
llx3
cose31x
cose4
TT0
2
0
2
F2
3F
20 (1.6-16)
mit 20t4
cF
lρλ
= .
Die FEM-Lösung wird mit nur zwei finiten Elementen unter Ausnutzung der Symme-trie ermittelt (siehe Bild 1.6-3).
A12 3
1 2
l le e
el = l /4
Bild 1.6-3 Berechnungsmodell
Die Wärmeleitfähigkeitsmatrix berechnet sich aus den beiden Elementbeiträgen zu
KmmsW
10105c
KmmW
1016,0
39
3
−
−
⋅=ρ
⋅=λ
1.6. Ein Einführungsbeispiel zur eindimensionalen Wärmeleitung 37
−−−
−λ
=110121
0114
ges lK
Die Einarbeitung der Randbedingung ( ) 00tT1 =≥ ergibt
−
−λ=
11124
ges lK .
Analog erhält man aus den Elementbeiträgen die Wärmekapazitätsmatrix
ρ=
210141012
41
c61
ges lC ,
und nach Einarbeitung der Randbedingung ergibt sich
ρ=
2114
41
c61
ges lC .
Somit erhält man ein gewöhnliches Dgl.-System 1. Ordnung
0gesges =+ TCTK & , (1.6-17)
das durch eine Zeitintegration gelöst werden kann. Dazu führen wir für T& eine einfa-che Differenzenapproximation ein (siehe Bild 1.6-.4).
t
T
t it i + 0,5 ∆ t t i +1
T
T
T
i
i +1
i +0,5
∆ t
Bild 1.6-4 Differenzenapproximation der Temperaturableitung
Die zentrale Differenzenformel lautet:
2TT
T 1iii 2
1+
+
+= (1.6-18a)
tTT
T i1ii 2
1
∆−
= ++
& (1.6-18b)
38 1. Elementare Einführung in die FEM
Damit ergibt sich aus (1.6-17)
( ) ( ) 0t
121
i1iges1iiges =−∆
++ ++ TTCTTK (1.6-19)
und nach Umformung das Lösungsschema
i1igesges t21
t21
TKCTKC
∆−=
∆+ + , (1.6-20)
mit ,...2,1i = wird die Lösung an den Zeitschritten bezeichnet. Aus den Randbedin-gungen folgt der Startvektor
=
100100
0T .
Mit den Zahlenwerten für unser Beispiel und der Schrittweite s1,0t =∆ ergibt sich
i1
1i
i1i
BTAT
BTAT−
+
+
=
=
mit
=
=∆
−
989644,0012409,00062045,0989644,0
s1,0t
1
BA
Damit errechnet man folgende Lösungen:
Zeit TA (t)[ s ] exakte Lösung FE-Lösung
(2 Elemente)0 100,00 100,0010 99,16 100,0420 91,43 87,5230 80,36 75,0240 69,58 64,08M M MM M M
100 28,30 24,8
Deutlich erkennbar ist ein für Wärmeschockprobleme typisches Oszillieren der Lö-sung, das seine Ursache in der Ungenauigkeit des numerischen Lösungsverfahrens hat(siehe Bild 1.6-5).
1.6. Ein Einführungsbeispiel zur eindimensionalen Wärmeleitung 39
0
20
40
60
80
100
120
0 5 10 15 20 30 40 50 60 70 80 90 100
t
TA
Bild 1.4.3-5 Zeitlicher Verlauf der Temperatur am Punkt A der Wand
SchrittweitenabschätzungEin Problem bei der Anwendung von Zeitintegrationsschemen ist die Schrittweitenab-schätzung, so daß eine ausreichende Genauigkeit der Lösung gesichert ist (siehe dazuauch die Vorlesung Computer-Numerik). Es läßt sich zeigen, daß die Schrittweite beider Zentralen Differenzenmethode die Bedingung
ω≤∆
2t (1.6-21)
erfüllen muß, um eine stabile Lösung zu sichern. Dabei ist ω der größte Eigenwertvon ( )KC 1− , d.h.
( )KC 1max
−ω=ω . (1.6-22)
Bei Systemen mit vielen Freiheitsgraden (mehrere Millionen sind heute keine Selten-heit mehr) ist die Ermittlung dieses Eigenwertes extrem aufwendig. Man kann zeigen,daß dieser Wert relativ einfach auf der Basis von Elementinformationen, d.h. mit eKund eC , abschätzt werden kann. Es gilt
( )( )[ ]( )[ ]emine
emaxe
min
max1max min
max
)()(
C
K
CK
KCω
ω≤
ωω
=ω − (1.6-23)
Für unseren eindimensionalen Fall ergibt sich aus den Elementgleichungen
( ) le
emax 2λ
=ω K (1.6-24a)
( ) eemin c61
l ρ=ω C (1.6-24b)
Damit liefert (1.6-23)
( ) 2e
1max c
12lρ
λ≤ω − KC , (1.6-25)
40 1. Elementare Einführung in die FEM
woraus sich für die Schrittweite folgende Abschätzung ergibt:
λρ
=
ρλ
≤∆6c
c21
12
1t
2e
2e
l
l
(1.6-26)
Aus dieser Gleichung wird ersichtlich, daß die Elementlänge le sowie die Material-kenngrößen λ und cρ Einfluß auf die Schrittweite haben.
1.7. Kleiner historischer Rückblick auf die frühen Anfänge der FEM 41
1.7. Kleiner historischer Rückblick in die frühen Anfänge der FEMAls Beginn der Finite-Element-Entwicklung wird allgemein das Jahr 1956 angesehen,in dem Turner, Clough, Martin und Topp eine Arbeit veröffentlichten, in der für denebenen Spannungszustand das Dreieckselement mit 3 Knoten und 6 Freiheitsgraden (2je Knoten) abgeleitet, rechentechnisch realisiert wurde an an Beispielen getestet wurde(nach einem Bericht von Clough wurde diese Arbeit bereits 1953 fertiggestellt aberdurch Turner erst ca. 3 danach zur Veröffentlichung eingereicht). Eine Reihe vonNutzrechnungen an Deltaflügeln von Kampfflugzeigen bewiesen dieLeistungsfähigkeit des Verfahrens, auch wenn die Zahl der Unbekannten nur etwa 50betrug.
Die Entwicklung der Matrizenmethoden der Stabstatik (Langefors, Argyris) undder Rechentechnik waren die wesentlichen Voraussetzungen für diese Pionierleistungder Ingenieure. Das Verfahren war tatsächlich mehr ingenieurmäßig intuitiv alsmathematisch begründet. Das wird deutlich, wenn man die Ableitung derElementsteifigkeitsmatrix des Dreieckselementes betrachtet. Dazu wurdeangenommen, daß im Dreieck ein konstanter Verzerrungszustand vorliegt:
( )ε σ µσ∂∂x x ya
Eux
= = − =1
( )ε σ µσ∂∂y y xb
Euy
= = − =1 (1.7-1)
γ∂∂
∂∂xy xyc
Guy
ux
= = = +1
τ
Durch Integration lassen sich daraus die Verschiebungen ermitteln:
u ax Ay B= + + (1.7-2)
( )v by c A x C= + − +
A,B,C sind Integrationskonstanten (Starrkörperbewegung). Die 6 Konstanten in (1.7-2) lassen sich durch die Knotenverschiebungen ausdrücken, indem dieKnotenkoordinaten eingesetzt werden. Setzt man das so gefundene Ergebnis in (1.7-1)ein, erhält man die Spannungen als Funktion der Knotenverschiebungen. DerZusammenhang zwischen den Spannungen und den Knotenkräften wird aus einfachenGleichgewichtsbetrachtungen am Dreieckselementen erhalten. Damit wurden dieKnotenkräfte näherungsweise durch die Knotenverschiebungen ausdrücken und manerhielt so die bekannte Matrizengleichung für ein finites Element:
fKv = (1.7-3)
Für diesen Fall ist das Ergebnis identisch mit der Lösung, die sich über dieMinimierung des elastischen Potentials ergibt. Verbesserte Elemente sind auf demoben skizzierten Weg nicht ohne weiteres abzuleiten. Die Anwendung des Prinzipsvom Minimum des elastischen Potentials und die Erkenntnis des Zusammenhanges mitdem Ritzschen Verfahren ermöglichte eine mathematische Fundierung und gezielteWeiterentwicklung der Methode der finiten Elemente.Bemerkenswert ist allerdings, daß diese Zusammenhänge längst erkannt und vonCourant in seiner hochinteressanten Arbeit über " Variationsmethoden zur Lösung vonstatischen und dynamischen Problemen" bereits 1943 dargelegt wurden. In der Arbeit
42 1. Elementare Einführung in die FEM
würdigt Courant ausdrücklich die Verdienste von W. Ritz und widmet der Anwendungdes Ritzschen Verfahrens breiten Raum. Courant richtet in der Veröffentlichungbesonderes Augenmerk auf die Möglichkeiten der praktischen Nutzung desVerfahrens. Er beginnt die Arbeit mit dem Zitat eines Ausspruchs von Henri Poincare,das wegen seiner Aktualität hier wiedergegeben werden soll:
H. POINCARE: Die "Lösung eines mathematischen Problems" ist eine Phrase mitunbestimmter Bedeutung. Reine Mathematiker sind irgendwann damit zufrieden, wennsie zeigen können, daß die Nichtexistenz einer Lösung einen logischen Widerspruchenthält, während die Ingenieure ein numerisches Resultat als einziges vernünftigesErgebnis ansehen. Solch einseitigen Orientierungen scheinen mehr menschliche alsobjektive Grenzen zu reflektieren. Die Mathematik ist in sich selbst eine unteilbareorganische Einheit, bestehend aus theoretischer Betrachtung und aktiver Anwendung.
Courant verweist dann auf Gauss und Thompson, die dafür echte Beispiele waren.Courant selbst wendete das Ritzsche Verfahren zur Lösung eines Torsionsproblemsan, wobei er das betrachtete Gebiet trianguliert und für jedes Gebiet einenNäherungsansatz einführt. Er betrachtet sogar verschiedene Ansatzfunktionen undverschieden feine Versetzungen. Er weist abschließend darauf hin, daß er in einergesonderten Publikation die Anwendung dieser Methode auf Platten und andereProbleme, die auch höhere Ableitungen enthalten, darlegen will. Es ist bedauerlich,daß diese Arbeit so wenig Beachtung gefunden hat. Die Zeit der Veröffentlichung magein Grund dafür sein. Die Kenntnis der Arbeit von Courant hätte sicher eineBeschleunigung der Entwicklung der FEM ermöglicht, Fehlentwicklungen hättenvermieden werden können.
top related