Diplomarbeit CalculiX Hiller

Vorgelegt von: Jörg Hiller (183240) Externer Betreuer: Dr.-Ing. Jan Reger Hochschul Betreuer: Prof. Dr.-Ing. Michael Wahle

Evaluierung von Open-Source Finite-Elemente Solvern für den Einsatz in Optimierungsprozessen

Evaluierung von Open-Source FiniteElemente Solvern für den Einsatz in Optimierungsprozessen

Sperrvermerk

Sperrvermerk Die vorliegende Diplomarbeit beinhaltet vertrauliche Informationen der Firma P+Z Engineering GmbH. Die Weitergabe des Inhaltes der Arbeit im Gesamten oder in Auszügen ist grundsätzlich untersagt. Es dürfen keine Kopien oder Abschriften angefertigt werden. Ausnahmen bedürfen der schriftlichen Genehmigung der Firma P+Z Engineering GmbH. Es wird jedoch ausdrücklich darauf hingewiesen, dass die im Anhang dargestellten Beispiele (Test1 – Test3 mit den dazugehörigen Ergebnissen und Abbildungen im Hauptteil) vom Betreuer, Herrn Prof. Wahle, zu Lehrzwecken genutzt werden dürfen.

München, Apri 2008

Dipl.- Ing. Jan Reger P+Z Engineering GmbH Anton-Ditt-Bogen 3 D - 80939 München E-mail: [email protected]

Fachhochschule Aachen

Jörg Hiller, Mat.-Nr. 183240

I

Erklärung

Erklärung Hiermit erkläre ich, Jörg Hiller, geboren am 22.02.1978, in Düren, die vorliegende Diplomarbeit selbständig angefertigt zu haben. Es wurden keine anderen, als die angegebenen Quellen und Hilfsmittel benutzt.



II

Kurzfassung

Kurzfassung

In dieser Arbeit werden Open-Source Finite-Elemente-Programme auf ihre Einsatzfähigkeit in Optimierungsprozessen untersucht.

Um eine Optimierung schnell durchzuführen, wird eine große Anzahl an Simulationen benötigt. Häufig stehen jedoch nicht genügend Lizenzen zu Verfügung. Hier bietet sich der Einsatz von frei verfügbaren Programmen an.

In

Rahmen

der

Diplomarbeit

werden

verschiedene

Programme

aus

unterschiedlichen Anwendungsgebieten kurz vorgestellt. Geeignete Programme werden anhand von Standard-Beispielsimulationen auf ihre Ergebnisqualität untersucht und mit analytischen Lösungen und den Ergebnissen der kommerziellen Programme NASTRAN und Abaqus verglichen. Anschließend wird die maximale Modellgröße auf verschiedenen Computersystemen bestimmt und die benötigte Rechenzeit sowie der Speicherverbrauch getestet. Zum Ende wird die Software in einen Optimierungsprozess integriert.

Als Ergebnis der Arbeit wurde eine Software gefunden, die sich ohne Probleme in schon bestehende Prozesse integrieren lässt. Das Programm ist kompatibel zu kommerziellen

Preprocessoren

und

hat

ein

einfaches

Inputformat.

Die

Lösungsqualität ist gut und größere Modelle können in einer angemessenen Zeit berechnet werden. Ferner lässt sich die Software in Optimierungsprozesse integrieren.



III

Inhaltsverzeichnis

Inhaltsverzeichnis Abbildungsverzeichnis .......................................................................................................... 2 Tabellenverzeichnis .............................................................................................................. 3 1 Einleitung..................................................................................................................... 4 2 Anforderungen an Finite-Elemente Programme............................................................. 5 2.1 Kriterien zur Auswahl der Open-Source Programme ................................................... 6 2.2 Übersicht Open-Source FE-Programme ..................................................................... 7 2.2.1 CalculiX ...................................................................................................... 10 2.2.2 Open Foam................................................................................................... 11 2.2.3 Impact ......................................................................................................... 12 2.2.4 MBDyn........................................................................................................ 14 2.2.5 OOFEM....................................................................................................... 15 2.3 Besonderheiten von CalculiX/OOFEM.................................................................... 17 2.3.1 CalculiX ...................................................................................................... 17 2.3.2 OOFEM....................................................................................................... 20 3 Evaluation der FE-Programme.................................................................................... 22 3.1 Durchbiegung Blattfeder........................................................................................ 23 3.1.1 Fazit ............................................................................................................ 26 3.2 Biegeeigenfrequenz Balken.................................................................................... 28 3.2.1 Fazit ............................................................................................................ 31 3.3 Eigenfrequenzen Platte .......................................................................................... 32 3.3.1 Fazit ............................................................................................................ 35 3.4 Test Modellgröße.................................................................................................. 36 3.4.1 Fazit ............................................................................................................ 38 3.5 Zusammenfassung ................................................................................................ 40 4 Optimierung ............................................................................................................... 42 4.1 Einführung Sensitivitätsanalyse / Robustheitsanalyse ................................................ 43 4.1.1 Optimierungsprozess Finite-Elemente Querlenker............................................. 44 4.2 Fazit der Optimierung............................................................................................ 46 5 Steifigkeitsberechnung Motorhaube ............................................................................ 49 5.1 Kopfaufpralldefinition ........................................................................................... 50 5.1.1 Modellbeschreibung FE-Motorhaube............................................................... 50 5.2 Diskussion der Simulationsergebnisse ..................................................................... 51 6 Diskussion der Ergebnisse ........................................................................................... 54 A1 Open-Source FE-Programme .................................................................................. 57 A2 Inputdecks.............................................................................................................. 60 A2.1 Test 1 .................................................................................................................. 60 A2.2 Test 2 .................................................................................................................. 63 A2.3 Test 3 .................................................................................................................. 65 A3 Analytische Herleitung ............................................................................................ 67 A3.1 Test 1. Durchbiegung Blattfeder ............................................................................. 68 A3.2 Test 2. Eigenfrequenz Biegebalken ......................................................................... 72 A3.3 Test 3. Eigenfrequenzen fest eingespannte Platte ...................................................... 79 A4 Literaturverzeichnis................................................................................................ 84



1

Abbildungsverzeichnis

Abbildungsverzeichnis Abb. 2-1: Postprocessor CalculiX GraphiX [http://www.calculix.de/]................................... 10 Abb. 2-2: Pre-/Postprocessing Open Foam [http://www.opencfd.co.uk/openfoam/doc/] ..... 12 Abb. 2-3: Pre-/Postprocessor Impact [http://impact.sourceforge.net/]................................. 13 Abb. 2-4: Pre-/Postprocessing Blender [http://www.osengineer.com/blendedmbdynusersguide.htm].............................................. 14 Abb. 2-5: Postprocessing OOFEM [http://www.oofem.org/]................................................ 15 Abb. 2-6: Expansion Plattenelemente ................................................................................ 18 Abb. 2-7: Expansion Balkenelement .................................................................................. 19 Abb. 3-1: FE-Modell Test 1 ................................................................................................ 23 Abb. 3-2: Einfluss Elementdiskretisierung auf Verschiebung (Bild: CalculiX-GraphiX) ....... 24 Abb. 3-3: Vergleich analytische Lösung / Open-Source ..................................................... 25 Abb. 3-4: Varianten NASTRAN (Schubweich) / Open-Source............................................ 26 Abb. 3-5: Vergleich aller FE-Ergebnisse mit der analytischen Lösung. .............................. 27 Abb. 3-6: FE-Model Test 2................................................................................................. 28 Abb. 3-7: Eigenformen CalculiX (Bild: CalculiX-Graphix) ................................................... 29 Abb. 3-8: Vergleich analytische Lösung (nach Timoshenko und Bernoulli) / CalculiX......... 30 Abb. 3-9: Vergleich NASTRAN / CalculiX........................................................................... 30 Abb. 3-10: FE-Modell Test 3 .............................................................................................. 32 Abb. 3-11: Eigenformen CalculiX (Bild:CalculiX-CrunchiX) ................................................ 33 Abb. 3-12: Eigenformen NASTRAN (Bild: Patran).............................................................. 33 Abb. 3-13: Vergleich Eigenfrequenzen Hexaeder- / Schalenmodell ................................... 34 Abb. 3-14: Speicherbedarf / Modellgröße (64Bit-System) .................................................. 37 Abb. 3-15: Rechenzeit: Solver / Modellgröße..................................................................... 38 Abb. 4-1: "Robust Design" ................................................................................................. 43 Abb. 4-2: Querlenkermodell ............................................................................................... 44 Abb. 4-3: Ausschnitt Inputdeck CalculiX............................................................................. 45 Abb. 4-4: Dakotainput ........................................................................................................ 46 Abb. 4-5: Änderung E-Modul durch DAKOTA .................................................................... 47 Abb. 4-6: Scatterplot .......................................................................................................... 48 Abb. 5-1: Kopfaufprall [www.euroncap.com] ...................................................................... 50 Abb. 5-2: FE-Modell Motorhaube ....................................................................................... 51 Abb. 5-3: Verschiebung Motorhaube.................................................................................. 52 Abb. A3-1: Blattfeder.......................................................................................................... 68 Abb. A3-2: Bernoulli Annahme........................................................................................... 69 Abb. A3-3: Zusammenhang Verdrehung und Winkeländerung .......................................... 70 Abb. A3-4: Gesamtverschiebung ....................................................................................... 71 Abb. A3-5: Gedrungener Balken ........................................................................................ 72 Abb. A3-6: Verdrehung Balken .......................................................................................... 72 Abb. A3-7: Balkenelement ................................................................................................. 72 Abb. A3-8: Bewegungsgleichung Balken ........................................................................... 74 Abb. A3-9: Vergleich der Balkentheorien [GRA91]............................................................. 76 Abb. A3-10: Platte.............................................................................................................. 79 Abb. A3-11: Bewegungsgleichung Platte ........................................................................... 79 Abb. A3-12: Eigenformen einer quadratischen Platte mit Seitenlänge a und Biegesteifigkeit D[LEI61]............................................................................................................................. 80 Abb. A3-13: Schub- und Rotationsträgheitseinfluss auf die Wellengeschwindigkeit [MIN61] .......................................................................................................................................... 81



2

Tabellenverzeichnis

Tabellenverzeichnis Tabelle 2-1. Übersicht Open-Source Programme ................................................................ 9 Tabelle 3-1: Modellgrößen des Benchmarks...................................................................... 36 Tabelle A3-1: Charakteristische Gleichungen für verschiedene Randbedingungen [GRO99] .......................................................................................................................................... 74 Tabelle A3-2: Vergleich der Eigenfrequenzen nach Quellen .............................................. 77 Tabelle A3-3: Eigenfrequenzen.......................................................................................... 80 Tabelle A3-4: Einsatzgebiete Plattentheorie ...................................................................... 82



3

1 Einleitung

1 Einleitung Die Finite-Elemente Methode (FEM) ist derzeit das wichtigste Werkzeug im Ingenieurbereich. Erstmals wurde die FEM im Bereich der Struktur- und Strömungsmechanik 1956 von Boeing zur Untersuchung gepfeilter Flugzeugflügel eingesetzt. Seitdem hat es eine rasante Entwicklung in den Simulationsverfahren gegeben, die auf die steigenden Rechnerkapazitäten bei sinkenden Kosten, den Einzug

computergestützter

Produktentwicklung

und

das

rechnergestützte

Konstruieren beruht. Die Bedeutung der Finite-Elemente Methode ergibt sich aus dem großen Einsatzspektrum. So können schon in der Entwurfsphase verschiedene Konzepte auf ihre Tauglichkeit untersucht und die unterschiedlichsten Einflüsse sondiert werden. Das Computer Aided Engineering (CAE) trägt so dazu bei, die Entwicklungskosten und –zeit zu senken, die Anzahl von teuren Versuchsreihen zu reduzieren und das Produkt optimal auszulegen. Voraussetzung für die Erfüllung dieser Ziele ist dabei eine leistungsfähige Soft- und Hardware, hinreichende Einarbeitungszeit in die verwendeten Softwareprodukte sowie eine kritische Beurteilung der berechneten Ergebnisse. Um alle Vorteile des virtuellen Prototyping ausnutzen zu können, ist ein großes Softwareportfolio notwendig, wodurch zum Teil enorme Lizenzkosten anfallen. Dies ist insbesondere bei Optimierungen der Fall, wo nach Möglichkeit viele Simulationen parallel laufen, um ein schnelles und aussagekräftiges

Ergebnis

zu

erhalten.

Um

die

Lizenzkosten

bei

Optimierungsprozessen niedrig zu halten bietet sich der Einsatz von Open-SourceSoftware an. Hierbei muss beachtet werden, dass der Wegfall der Lizenzkosten durch aufwendiges Einarbeiten in das neue Softwareprodukt oder durch umständliches Integrieren in schon bestehende Prozesse relativiert werden kann. Somit müssen bei der Verwendung von Open-Source-Softwareprogrammen die Faktoren Einarbeitung, Integration, Ergebnisqualität und Kompatibilität bedacht werden. In dieser Arbeit wird einleitend ein Überblick über den Markt der Open-SourceFinite-Elemente Programme gegeben und geeignete Pakete weiter untersucht. Dazu werden einige Testbeispiele simuliert und die Ergebnisqualität anhand Berechungen mit kommerziellen Programmen und einem analytischen Ergebnis bewertet. Abschließend wird die Integration in einen Optimierungsprozess getestet. Fachhochschule Aachen


4

2 Anforderungen an Finite-Elemente Programme

2 Anforderungen an Finite-Elemente Programme In der Automobilindustrie ist heute das Simulieren von Gesamtfahrzeugen oder einzelnen Komponenten integrativer Bestandteil des Entwicklungsprozesses. Frühzeitig werden neue Entwicklungen in der Softwaretechnologie angewendet und vorhandene Programme durch neue Anforderungen an den Hersteller erweitert. So umfasst das Einsatzspektrum an Simulationen zum Beispiel: •

Fahreigenschaftsimulationen.

•

Numerische Strömungssimulationen zur Bewertung des Klimakomforts.

•

Festigkeitsberechnungen, um das Bauteilversagen vorherzusagen.

•

Kontaktanalysen (Crash) für Unfallsimulationen.

•

Akustikuntersuchungen zur Geräuschminderung.

•

Dynamiksimulationen, um das Schwingungsverhalten zu bewerten.

Für jede der oben aufgeführten Simulationen kommen jeweils spezielle Programme aus verschiedenen Bereichen der nummerischen Simulation zum Einsatz. Der Großteil an Berechnungen mittels der Finite-Elemente-Theorie befasst sich jedoch mit dem Gebiet der Statik und Dynamik. Jeder Ingenieur lernt an der Hochschule die analytischen Grundlagen, um Festigkeitsprobleme und das dynamische Verhalten von Körpern zu berechnen. Um den Studenten einen kostenlosen Zugang zur Finite-Elemente-Theorie und diese durch Übungen zu verfestigen, werden häufig an den Universitäten eigene Programme entwickelt. Daher ist der Umfang an Software, mit der sich Probleme aus diesen Bereichen berechnen lassen, groß.



5


2.1

Kriterien zur Auswahl der Open-Source Programme

Bei der Recherche zu Open-Source-Programmen wurden weit über 100 Softwarepakete gefunden (siehe Anhang A1). Grundsätzlich lassen sich die Programme in drei Kategorien einteilen. Zum einen sind das die Softwarepakete aus dem Umfeld einer Universität. Hier liegt der Schwerpunkt in der Simulation von kleinen Problemen, die sich auch mit einer Handrechnung lösen lassen, um den Studenten den Einstieg in die Finite-ElementeTheorie zu ermöglichen. Die Benutzerfreundlichkeit bei der Erstellung des Inputs ist dabei nicht immer gegeben. Auch sind die Darstellungsmöglichkeiten der Ergebnisse beschränkt und es gibt keine Importschnittstellen zu anderen FEProgrammen. Die zweite Kategorie sind die „semi-professionellen“ Open-Source Projekte. Hierbei handelt es sich um Programme, die schon länger existieren und durch Webseiten einer breiten Öffentlichkeit zugänglich gemacht wurden. Dies hat zur Folge, dass immer mehr Personen an der Entwicklung beteiligt sind und eine breite Benutzergemeinde durch Anwenden des Programms zur Fehlerbeseitigung beitragen. Diese Programme besitzen meist Importschnittstellen zu anderen Programmen und ermöglichen eine grafische Darstellung der Ergebnisse. Der Code lässt meist das Parallelisieren der Simulation zu und verfügt über effektiv arbeitende Gleichungslöser, um auch große Probleme zu berechnen. Als letzte Kategorie gibt es die „low-cost“ Programme. Diese Software bietet einen mit kommerziellen Programmen vergleichbaren Funktionsumfang und liefert meist ausgereifte Pre- und Postprocessoren mit. Für die Verwendung der Software und den Support fallen geringe Lizenzkosten an. Bei den für diese Diplomarbeit untersuchten Programmen handelt es sich ausschließlich um frei verfügbare (Open-Source) Produkte. Zur Auswahl geeigneter Programme werden folgende Kriterien formuliert, die in Teilen erfüllt werden sollten: •

Importfunktion für CAD-Formate, wie z.B.: IGES, STEP

•

Importfilter für Finite-Elemente Netze aus anderen Programmen, wie NASTRAN,

Abaqus,

ANSYS.

Hierdurch

besteht

die

Möglichkeit,

kommerzielle Preprocessoren für die Vernetzung zu nutzen. Vorteilhaft an



6


diesem Vorgehen ist, dass diese Programme mehr Funktionen für eine Vernetzung zur Verfügung stellen. •

Einfaches und klar definiertes Inputformat, um die Einarbeitung zu erleichtern und Fehler zu vermeiden.

•

Eigener Preprocessor oder kompatibel zu kommerziellen Programmen.

•

Eine umfangreiche Elementbibliothek. Mindestens sollten hier Schalen- und Volumenelemente vorhanden sein, die für eine Simulation im Raum (6 Freiheitsgrade) geeignet sind.

•

Linear-elastisches Werkstoffverhalten.

•

Fähigkeit stationäre, statische, transiente und dynamische Probleme zu simulieren.

•

Eigener Postprocessor oder kompatibel zu kommerziellen Programmen.

•

Multiprozessorfähigkeit oder die Möglichkeit, Modelle parallel zu simulieren, um die teilweise sehr großen Modelle berechnen zu können.

2.2

Übersicht Open-Source FE-Programme

Die im vorherigen Kapitel formulierten Anforderungen konnten nur wenige Programme erfüllen. Viele Programme sind nur für die Berechnung von kleineren Problemen geeignet, was daran liegt, dass der Hauptanteil der Programme in dem Umfeld

einer

Universität

entwickelt

wurde.

Somit

waren

meistens

die

Grundvoraussetzungen wie Preprocessorfähigkeit, d.h. die Möglichkeit aus CADDaten Finite-Elemente Netze zu erzeugen oder schon vorhandene Netze aus anderen Formaten zu importieren, nicht gegeben. Oder Simulationen können nur in der Ebene und nicht im dreidimensionalen Raum durchgeführt werden. Aus den gefundenen Programmen werden nur zwei als aussichtsreiche Kandidaten weiter untersucht. Das Programm OOFEM, welches an der Prager Universität für Maschinenbau entwickelt wird und CalculiX, das von Mitarbeitern bei MTU in München entwickelt wird. In Tabelle 2-1 findet sich eine Funktionsübersicht der Programme. Hier sind auch Open-Source Programme aus anderen Anwendungsgebieten aufgeführt, die im Rahmen der Recherche gefunden wurden. Diese zeichnen sich durch einen sehr fortgeschrittenen Entwicklungsstand aus oder sind die einzigen in ihrem Fachhochschule Aachen


7


Einsatzbereich. Eine kurze Beschreibung der aufgeführten Programme findet sich in den nachfolgenden Kapiteln. Als sehr gutes Beispiel einer Open-Source Software muss das Programm Code_Aster kurz erwähnt werden. Code-Aster wird von Electricitè de France (EDF), dem französischen nationalen Energieversorgungsunternehmen, seit über 10 Jahren entwickelt und seit 2001 unter der GNU-Lizenz samt dem Quellcode zur Verfügung gestellt wird. Das Programm kann unter http://www.code-aster.org/ herunter

geladen

werden.

Es

bietet

den

vollen

Funktionsumfang

einer

kommerziellen Software, wie z.B. MD NASTRAN, an und ist zu vielen Pre- und Postprocessoren kompatibel. Ein gravierender Nachteil ist die französische Dokumentation. Die notwendigen französischen Begriffe für das Inputformat und den mitgelieferten Pre-/Postprocesor Salome lassen sich sicherlich erlernen, machen einen Einstieg aber nicht gerade einfach. Zur Zeit gibt es große Anstrengungen, um Code_Aster einer englischsprachigen Benutzergemeinde zugänglich zu machen. Hier ist das Projekt CAElinux [http://caelinux.com/CMS/] sehr weit fortgeschritten und bietet einen ersten Einstieg in Code_Aster. Aufgrund der Sprachbarriere und des großen Funktionsumfangs wird Code_Aster nicht weiter in der Diplomarbeit untersucht.



8

www.calculix.de

9. August 2007

Homepage

Letztes Update

Impact



Kompatible VDA, STL CAD-Formate

Kompatible Inputformate

--

NASTRAN, Abaqus, Fluent, STAR-CD Ansys, Code_Aster, Duns, ISAAC, OpenFOAM

(Paraview, Ensight)

--

NASTRAN, UNV (Salome), Gmsh Version 2.

(GiD)

Postprocessor

√

-(GiD)

--

Nichtlinear Elastisch, Linear Elastisch

Grafische Postprocessing Pre- und BenutzerPostprocessing oberfläche Preprocessor (Kommandozeile) √

Linear Elastisch, Hyperelastisch, Materialmodelle Viskoelastisch, Elastoplatisch

1-D, 2-D, 3-D, Spezialelemente

C++ Toolbox, Finite-Volumen Methode

1-D, 2-D, 3-D, Elementtypen Spezialelemente

11. Juli 2007

1.4.1 0.7.5 StrömungsCrashsimulation, mechanik, Metallumformung, Wärmeströmungen, Topologieoptimierung Elektrostatik, Verbrennungsvorgänge

3. August 2007

www.opencfd.co.uk/ http://impact.sourceforge.net/

Open Foam

1.7 Statik, Dynamik, Nichtlineare Probleme, Kontakt, Wärmeströmungen, Einsatzgebiet Inkompressible Gase und Flüssigkeiten

Version

Calculix

Software

--

--

(Blender, MBDyn sim suite) (Blender, GnuPlot, MBDyn sim suite)

Pre- und Postprocessing

(Abbildung der flexiblen Körper) Linear Elastisch, Viskoelastisch, Elastoplatisch

Starre und flexible Körper

1.3.3 Beta Mehr Körper Simulation

26. Januar 2008

www.aero.polimi.it/~mbdyn

MBDyn

--

--

(t3d, Targe2) OOFEG, (Paraview)

Postprocessing

Linear Elastisch, Plastisch, Spezielle

1-D, 2-D, (3-D)

www.oofem.org 4. November 2006 1.7 Statik (2-D), Dynamik (2-D)

OOFEM


Tabelle 2-1. Übersicht Open-Source Programme

9


2.2.1 CalculiX Die Simulationsumgebung CalculiX kann statische, dynamische, nichtlineare Probleme sowie Wärmeleitung, inkompressible Flüssigkeiten und Kontaktprobleme berechnen. In dem Softwarepaket sind der Solver CalculiX CrunchiX und der Pre/Postprocessor

CalculiX

GraphiX

enthalten

(Abb.

2-1).

Zur

Lösung

der

Gleichungssysteme stehen zwei iterative und ein direktes Verfahren zur Auswahl. Falls die Programmexecutables aus den Sourcedateien kompiliert werden, bestehen zusätzlich die Optionen einen „out-of-core“ Solver einzubinden, Dual-Core CPUs anzusteuern und das Message Passing Interface (MPI) für Clustersysteme zu implementieren. Der Simulationsinput gleicht dem Abaqusformat, wodurch auch kommerzielle Preprocessoren für die Erstellung des Inputdecks verwendet werden können. Der mitgelieferte Preprocessor hat keine grafische Benutzeroberfläche. Die Vernetzung von Strukturen erfolgt daher innerhalb der Linuxshell. Einige Eingaben werden in dem Grafikfenster des Postprocessor durchgeführt und der Fortschritt kann hier verfolgt werden. Für komplexere Strukturen kann daher nur das automatische

Vernetzen

mit

Tetraederelementen

angewandt

werden.

Die

Elementbiliothek ist sehr umfangreich und beinhaltet sogar Sonderelementen wie z.B.: Gapelemente oder Multi Point Constraint Equations (MPC).

Abb. 2-1: Postprocessor CalculiX GraphiX [http://www.calculix.de/]

Die

Materialdatenbank

umfasst

linear-elastisches,

hyperelastisches,

viskoelastisches und elastoplatisches Materialverhalten. Inputdecks aus den Fachhochschule Aachen


10


Simulationsprogrammen NASTRAN, Abaqus, Ansys, Code_Aster, Duns, ISAAC und OpenFOAM können importiert werden. Für die Darstellung der Ergebnisse durch CaluliX GraphiX stehen Fringeplots und Vektorplots oder „time histoy plots“ zur Verfügung. Animierte Strukturen werden in Form von Gif-Animationen exportiert. Auf der Homepage von CalculiX lassen sich Tutorials, Beispiele und die Dokumentationen zu CalculiX CrunchiX und CalculiX GraphiX herunterladen. Für Fragen, die nicht mit dem Manual geklärt werden, gibt es ein sehr aktives Forum.

2.2.2 Open Foam Open Foam ist ein Programm zur Simulation von kompressiblen / inkompressiblen Fluiden, Wärmekonvektionsströmungen und Elektrodynamik. Das Programm ist eine Softwaresammlung bestehend aus verschiedenen C++ Modulen. So besteht die Möglichkeit, die zur Verfügung stehenden Solver, Pre-/Postprocessoren, Datenvisualisierungsmodule

und

Netz-/Inputmanipulierungsprogramme

einzeln

anzusteuern oder aber die Bibliotheken im vordefinierten Prozess anzuwenden. Die Aufspaltung in einzelne Module ermöglicht dem Anwender eine Anpassung an seine speziellen Bedürfnisse. In Kombination mit der gängigen Programmiersprache C++ lassen sich so leicht neue Prozesse implementieren und das Programm erweitern, ohne dass die gesamte Softwarestruktur bearbeitet werden muss. Für den Import von vorhandenen Inputdateien stehen Schnittstellen zu den kommerziellen Programmen Fluent und Star-CD zur Verfügung. Dadurch wird das Erstellen neuer Simulationsaufgaben erheblich vereinfacht und das Einarbeiten in den Open Foam Input durch Vergleichen der verschiedenen Formate erleichtert. Finite-Elemente Netze können aus den Programmen ANSYS, CFX, Ideas, gmsh und Netgen importiert oder durch das implementierte Vernetzungsmodul erstellt und bearbeitet werden. Dafür wird die Struktur zuerst in einer ASCII-Datei definiert und anschließend automatisch vernetzt. Hier ist es von Vorteil, dass das Netz in einem kommerziellen Programm mit grafischer Benutzeroberfläche erstellt werden kann und hinterher in einem importfähigen Format eingelesen wird. Zum Manipulieren stehen

innerhalb

von

Open

Foam

grundsätzliche

Befehle

wie

z.B.

Zusammenführen, Anbinden, Überprüfen, Verfeinern, Trennen usw. zur Verfügung.



11


Abb. 2-2: Pre-/Postprocessing Open Foam [http://www.opencfd.co.uk/openfoam/doc/]

Um einen Simulationslauf zu definieren bzw. zu editieren wird auf den eigenen grafischen Pre-/Postprocessor zurückgegriffen (Abb. 2-2). Hier können alle Aufgaben von Vernetzung bis hin zum grafischen Darstellen der Ergebnisse grafisch abgearbeitet werden. Für das Postprocessing bietet Open Foam ein eigenes Tool an oder aber die Ergebnisse werden in VTK-Format heraus geschrieben

und

in

dem

Open-Source

Programm

Paraview

[http://www.paraview.org] bzw. dem kommerziellen Tool Ensight dargestellt. Für Fragen oder einen ersten Einstieg gibt es ein aktives Diskussionsforum sowie Tutorials und eine sehr ausführliche Dokumentation. 2.2.3 Impact Mit der Software Impact können Crashanalysen und Umformprozesse simuliert werden. Impact ist in JAVA geschrieben wodurch sehr leicht Simulationen auf verschiedenen

Computerarchitekturen

durchführen

lassen,

da

JAVA

eine

plattformunabhängige Programmiersprache ist. Das Programm befindet sich jedoch noch in einer frühen Entwicklungsphase. Somit sind die Kontaktdefinitionen, Materialmodelle sowie die Elementbibliothek noch nicht sehr fortgeschritten. Impact liefert

aber

bereits

eine

komplette

Pre-/Postprocessingumgebung

für

die

Vernetzung, Simulationsdefintion und die grafische Darstellung der Ergebnisse Fachhochschule Aachen


12


(Abb. 2-3). Weiterhin wird der Pre-/Postprocessor GiD und die Open-Source Vernetzungssoftware Gmsh unterstützt. Für das Importieren von vorhandenen Strukturen bietet Impact eine rudimentäre Schnittstelle zu NASTRAN, wobei Randbedingungen und Materialdefinitionen nicht mit übersetzt werden. In Zukunft wird hier aber mehr unterstützt, um eine große Zielgruppe durch das weit verbreitete Format anzusprechen bzw. größtmögliche Kompatibilität zu gewährleisten. Hierfür werden dann auch in einer fortgeschrittenen Version die Inputformate der kommerziellen Crashprogramme Radioss und PAM-CRASH importiert werden können.

Abb. 2-3: Pre-/Postprocessor Impact [http://impact.sourceforge.net/]

Als Zusatzfeature kann mit Impact eine Topologieoptimierung durchgeführt werden. Im ersten offiziellen stabilen Release soll dann auch noch eine Parametrisierung des

Modells

bei

der

Topologieoptimierung

möglich

sein.

Die

gesamte

Dokumentation und die Anwenderbeispiele befinden sich ebenfalls noch in einer Betaphase.

Aus

der

Mailingliste

lässt

sich

eine

sehr

aktive

Entwickler/Anwendergemeinde ableiten und Fragen werden relativ schnell und kompetent beantwortet.



13


2.2.4 MBDyn MBDyn steht für MultiBody Dynamics Software. Es handelt sich hierbei um eine Software zum Simulieren von Mehrkörpersystemen. Für das grafische Pre- und Postprocessing

stehen

verschiedene

Tools

anderer

Softwareprojekte

zur

Verfügung. Da es sich (noch) um einen Forschungscode handelt, haben diese Projekte die Entwicklung der grafischen Benutzeroberfläche übernommen (Abb. 2-4), während sich die Entwicklergruppe von MBDyn weiterhin um die Verbesserung des

Codes

kümmert.

Für

nähere

Informationen

über

die

Pre-

/Postprocessorprogramme siehe: •

MBDyn sim suite Projekt [http://mbdynsimsuite.sourceforge.net/]

•

Blender [http://www.osengineer.com/blendedmbdynintro.htm]

Neben den oben genannten Möglichkeiten steht noch eine Exportschnittstelle für die Ergebnisanimation zu dem kommerziellen Programm ADAMS/View zur Verfügung.

Abb. 2-4: Pre-/Postprocessing Blender [http://www.osengineer.com/blendedmbdynusersguide.htm]

Um die Ergebnisqualität zu verbessern, können flexible Körper mit Hilfe von NASTRAN generiert werden. Als Hilfe stehen neben den verschiedenen Mailinglisten eine ausführliche Dokumentation sowie Tutorials und Beispiele zur Verfügung.



14


2.2.5 OOFEM OOFEM steht für Object Oriented Finite Element Solver. Mit OOFEM können dynamische,

statische

und

strömungsmechanische

Berechungen

sowie

Wärmekonvektionsströmungen simuliert werden. OOFEM kann auf mehreren Rechnern parallel ausgeführt werden. Hierdurch besteht die Möglichkeit, auch größere Strukturen zu simulieren. Die Elementbibliothek beinhaltet alle gängigen Elementtypen

für

strukturmechanische

Berechnungen,

wobei

aber

das

Haupteinsatzgebiet der Elemente im zweidimensionalen Raum liegt. Für die räumliche Abbildung von Strukturen stehen Tetraederelemente zur Verfügung. In der

Materialdatenbank

finden

sich

elastische,

plastische

und

spezielle

Materialmodelle, wie z.B. zur Simulation der Rissfortbildung in einem Körper.

Abb. 2-5: Postprocessing OOFEM [http://www.oofem.org/]

Das Erstellen eines Finite-Elemente Netzes für OOFEM muss mit den externen Programmen t3d und Targe2 erfolgen. Beide Vernetzer fallen unter die Rubrik „low cost“-Programme und bieten für den nicht kommerziellen Einsatz kostenlose Versionen auf Anfrage an. Für das Postprocessing wird das Programm OOFEG mitgeliefert (Abb. 2-5). Hiermit lassen sich Contourplots, Animationen und Farbplots erstellen. Die Ergebnisse können mittels externer Programme, wie z.B. GNUPlot Fachhochschule Aachen


15


oder auch Excel, zu X-Y Diagrammen weiter verarbeitet werden. Als alternative Möglichkeit bietet sich das Exportieren der Ergebnisse in das VTK-Format an. Die Darstellung erfolgt dann mit der Open-Source Software Paraview. Hilfe findet man im Manual oder dem Userforum. Fragen im Forum werden aufgrund der kleinen Benutzer-gemeinde

im


Moment

überwiegend

vom

Entwickler


beantwortet.

16


2.3

Besonderheiten von CalculiX/OOFEM

Bei der Durchführung der Beispielsimulationen zur Klärung der Einsatzgrenzen und der

Ergebnisqualität

der

Programme

CalculiX

und

OOFEM

sind

einige

Besonderheiten aufgefallen. So ist zum Beispiel die in der Dokumentation von OOFEM beschriebene Elementbibliothek nicht so umfangreich bzw. einsatzfähig und der Postprocessor konnte selbst nach ausführlichem Kontakt mit dem Programmierern nicht zum Laufen gebracht werden. CalculiX verwendet eine spezielle Vorgehensweise, um Schalen- und Balkenelemente abzubilden. Schalenund

Balkenelemente

werden

zu

Durchbiegung richtig abzubilden,

Hexaederelementen können

expandiert.

nur Elemente mit

Um

die

quadratischer

Ansatzfunktion (Elementen mit Mittelknoten auf den Elementkanten) verwendet werden. Durch dieses Vorgehen wird die Modellgröße künstlich vergrößert. Für Strukturen mit Schalenelementen erhöht sich die Modellgröße um den Faktor 2,5 und für Balkenstrukturen sogar um den Faktor 6,7. Im folgenden Abschnitten werden einige Features und besondere Eigenheiten der Open-Source Programme aufgeführt.

2.3.1 CalculiX CalculiX verwendet bei der Abbildung von Schalen- und Balkenelemente eine besondere Strategie. Für Schalenelemente wird derselbe Ansatz wie in Abaqus (SC8R Continuum Shell Element) für dicke Schalen angesetzt. Zweidimensionale Elemente werden entlang der Schalenelementnormalen zu dreidimensionalen Hexaederelementen bzw. dreieckige Schalen zu Keilelementen expandiert. Dabei definiert das Originalschalenelement die Mittelfläche, von der aus jeweils die Hälfte der Schalendicke in die positive und die negative Normalenachse extrudiert wird (Abb. 2-6). Die Knotenanzahl für das Schalenelement erhöht sich hierdurch von acht Knoten auf zwanzig und damit um den Faktor 2,5. Somit wird die noch rechenfähige Modellgröße bei der Verwendung von Schalen- und Balkenelementen deutlich reduziert. Hier wird sehr schnell die Leistungsgrenze eines Computers erreicht und für die Simulation wird ein iteratives Verfahren notwendig. Eine Übersicht der in CalculiX enthaltenen Gleichungslöser und deren Einsatzgrenzen findet sich in Kapitel 3.4 Test Modellgröße. Fachhochschule Aachen


17


Abb. 2-6: Expansion Plattenelemente

Um die Verbindung der neu generierten Volumenelemente zu gewährleisten, werden sogenannte „knot“ eingeführt. Die neuen Knoten sind in einem „knot“ nummerisch zusammengefasst und es wird sicher gestellt, dass bei Belastung der rotatorischen Freiheitsgrade diese auf die translatorischen der Hexaederelemente umgerechnet werden. „Knots“ werden generiert, wenn •

benachbarte

Elemente

eine

unterschiedliche

Elementnormalenrichtung

haben, •

verschiedene Elementtypen miteinander verbunden sind,

•

die Dicke variiert,

•

eine Randbedingung auf den rotatorischen Freiheitsgraden vorliegt oder

•

ein Biege- oder Torsionsmoment angreift.

Wenn einer der oben aufgeführten Punkte erfüllt ist, wird bei Balkenelementen, wie für Schalenelemente, ein „knot“ generiert. Die Knotenanzahl erhöht sich dabei um den Faktor 6,7 für Balkenelemente und um den Faktor 2,5 bei Schalenelementen durch die Expansion (Abb. 2-7). Bei Belastung eines Schalenelementknotens mit einer Punkkraft wird die Last auf den Referenzknoten des „knots“ aufgebracht. Falls kein „knot“ vorhanden ist und ein Mittelknoten belastet wird, greift jeweils die Hälfte der Kraft auf die expandierten Knoten an. Für die äußeren Elementknoten wird die Fachhochschule Aachen


18


Kraft zu 1/6 - 2/3 - 1/6 aufgeteilt. Bei Belastung eines Balkenknotens mit einer Punktlast wird diese zu 1/4 - 1/4 - 1/4 - 1/4 für Elementmittenknoten bzw. zu (-1/12) - (1/3) - (-1/12) - (1/3) - (-1/12) - (1/3) - (-1/12) - (1/3) für Endknoten aufgeteilt. Durch die Generierung von Hexaederelementen für Balkenstrukturen (Abb. 2-7) können nur rechteckige und elliptische bzw. als Sonderfall der Ellipse kreisförmige Balkenquerschnitte verwendet werden.

Abb. 2-7: Expansion Balkenelement

Bei Modellen mit Schalen- oder Balkenelementen, wird wie oben erläutert, die zu berechnende Modellgröße deutlich erhöht. Hier kann schnell die zu Verfügung stehende

Rechnerkapazität

ausgeschöpft

sein.

Wenn

CalculiX

aus

dem

Sourcecode selber kompiliert wird, kann ein „out-of-core“ Solver implementiert werden. Vorteil hierbei ist, dass, falls die Rechnerkapazitäten nicht ausreichen, die schon berechneten Gleichungen auf die Festplatte ausgelagert werden. Nachteilig an dem Vorgehen ist die langsame Rechengeschwindigkeit, die durch den Schreibund Lesevorgang auf die Festplatte hervorgerufen wird. Das Implementierten des „out-of-core“ Solvers bzw. das Kompilieren von CaluliX aus den Sourcedateien erfordert jedoch einige Linux- und Programmierkenntnisse. Als weitere Möglichkeit zur

Simulation

großer

Computerclustersystemen

Modelle installiert

zu

simulieren

werden.

Auch

kann

CalculiX

auf

hierfür

müssen

die

Programmexecutables selber aus den Sourcedateien erstellt werden, damit das benötigte Message Passing Interface (MPI) zur Verfügung steht. Als letzte Option für eine effiziente Simulation von großen Modellen kann das Multiphreading aktiviert werden, wodurch bei mehr Prozessorkernsystemen beide CPUs verwendet werden. Da heutige Computersysteme immer häufiger mit Dual-Core CPUs (zwei Fachhochschule Aachen


19


Prozessorkerne) ausgeliefert werden, bietet sich diese Option an und sollte immer mit kompiliert werden. CalculiX verwendet für die Definition des Simulationsinputs ein dem kommerziellen Programm Abaqus ähnliches Format. Für Benutzer von Abaqus entsteht so keine Einarbeitungszeit. Vorteilhaft an dem Vorgehen ist, dass jeder kommerzielle Preprocessor für die Erstellung von CalculiX-Inputdateien benutzt werden kann. Kommerzielle

Programme

bieten

meist

eine

benutzerfreundliche

grafische

Oberfläche, so dass die Erstellung der Inputdateien einfacher fällt. Auch bieten kommerzielle Programme deutlich mehr Möglichkeiten an, um komplizierte Strukturen zu vernetzen. Somit lässt sich CaluliX in schon bestehende Prozesse einfach implementieren, ohne dass Ausfallzeiten oder Komplikationen für neue Anwender anfallen. Anzumerken ist jedoch, dass das verwendete Format einem älteren Abaqusstand entspricht, so dass für den Anfang immer im CalculiX-Manual nach der genauen Definition der Einträge geschaut werden muss.

2.3.2 OOFEM Bei OOFEM hat sich der sehr kompliziert zu erstellende Input als ein großes Hindernis herausgestellt. So muss zum Beispiel bei der Knotenpunktdefinition für jeden Knoten die Lasteinleitung und die Randbedingung angegeben werden (siehe Kapitel A2.1). In Kombination mit der sehr technisch gehaltenen Dokumentation unterlaufen so leicht Fehler im Inputdeck. Falls ein Fehler auftritt, erhält der Benutzer kein Feedback vom Programm über die Art des Problems, was das Beseitigen

des

Fehlers

erschwert.

Das

zu

OOFEM

kompatible

Vernetzungsprogramm targe2 konnte nach einer ausführlichen Internetrecherche nicht gefunden werden. Auf Nachfrage beim Entwickler von OOFEM hat sich herausgestellt, dass die Entwicklung von targe2 eingestellt wurde. Der alternative Preprocessor t3d besitzt keine Importschnittstelle für CAD-Formate, daher müssen Geometriedaten neu in t3d erstellt werden. Finite-Elemente Netze aus anderen Programmen können ebenfalls nicht verwendet werden. Eine Bewertung über die Benutzerfreundlichkeit des Preprocessors t3d kann nicht erfolgen, da die erforderlichen ASCII-Inputdateien für die durchgeführten Tests von Hand erzeugt wurden. Eine genaue Beschreibung des Programms findet sich auf der Website Fachhochschule Aachen


20


http://mech.fsv.cvut.cz/~dr/t3d.html. Im Verlauf der Arbeit hat sich bei der Durchführung der Beispieltests gezeigt, dass die Elementbibliothek weniger umfangreich als angegeben ist. Es wurde festgestellt, dass der Haupteinsatzbereich von OOFEM im zweidimensionalen Raum liegt, da hierfür die meisten Elemente implementiert sind. Zur Simulation von dreidimensionenaler Probleme stehen nur Tetraederelemente zur Verfügung und zum Berechnen von Schalen und Platten können nur Dreieckselemente verwendet werden. Bei beiden Elementtypen muss daher immer die schlechte Ergebnisqualität durch zu grobe Vernetzung beachtet werden. Für die Darstellung der Ergebnisse liefert OOFEM den Postprocessor OOFEG mit. Auf den verwendeten Computersystemen mit RedHat Linuxdistribution stellten sich einige Fehler ein. Die simulierten Finite-Elemente Netze konnten zwar dargestellt werden, aber die farbliche Kodierung der Ergebnisse funktionierte nicht. Auch ist der Export in das VTK-Format, um den Simulationslauf mittels Paraview darzustellen, noch nicht stabil. Auf Anfrage bei den Entwicklern hat sich herausgestellt, dass OOFEG im Moment komplett umgeschrieben wird und daher keine fehlerfrei lauffähige Version vorhanden ist. Die Probleme mit dem VTK-Format liegen an dem frühen Entwicklungsstadium der Schnittstelle. Abschließend sei darauf hingewiesen, dass OOFEM von der „Czech Technical University“ in Prag entwickelt wird. Der Grund für die Programmierung von OOFEM ist, den Studenten der technischen Fakultät zu einem kostenlosen Zugang zur Finite-Elemente-Theorie

zu

verhelfen.

Hiermit

soll

ermöglicht

werden

die

Vorlesungsinhalte durch praktische Anwendung zu vertiefen. OOFEM ist ein weiteres FE-Programm, das aus dem universitären Umfeld kommt und wurde nur deshalb einer tiefer gehenden Untersuchung unterzogen, um die Einsatzgrenzen zu deutlich fortgeschritteneren Programmen aufzuführen. Bei all der Kritik und den Einschränkungen sollte dies immer berücksichtig werden und es ändert nichts an der Tatsache, dass es sich hier um ein voll einsatzfähiges Finite-Elemente Paket handelt.



21

3 Evaluation der FE-Programme

3 Evaluation der FE-Programme Um einen Überblick über die Ergebnisqualität und Benutzerfreundlichkeit der OpenSource

Finite-Elemente

Programme

zu

erhalten,

werden

drei

Tests

mit

Standardsimulationsarten durchgeführt. Eine statische Simulation, bei der die Enddurchbiegung einer dünnen Schalenstruktur berechnet wird. Der zweite Test simuliert die ersten drei Biegeeigenfrequenzen eines Balkens. Im letzten Test sollen die ersten sechs Eigenfrequenzen einer quadratischen Plattenstruktur berechnet. Um einen Überblick über die Lösungsqualität der verwendeten Elementdefinition zu erhalten, werden die Finite-Elemente Netz der ersten beiden Tests iterativ verfeinert. Die Ergebnisse der Simulationen werden anschließend mit den Lösungen der kommerziellen Programme Abaqus und NASTRAN sowie einer analytischen Lösung verglichen. Bei der Eigenschwingungssimulation der quadratischen Platte wird

auf

eine

Verfeinerung

des

Netzes

verzichtet.

Hier

wird

die

Simulationsmodelle

die

Benutzerfreundlichkeit des Postprocessors CalculiX-CrunchiX untersucht. Zum

Schluss

werden

durch

verschieden

große

Einsatzgrenzen von CalculiX bestimmt. Die analytische Herleitungen zu den Tests findet sich im Anhang 3.



22


3.1

Durchbiegung Blattfeder

l 100mm Simuliert wird eine dünne schlanke Struktur, Länge zu Dicke ( = ) und einer t 2mm

Breite von 10mm , mit einseitiger fester Einspannung. Am gegenüberliegenden Ende wird die Struktur mit einer Einheitslast belastet (Abb. 3-1). Ausgewertet wird die maximale

Verschiebung

am

Lastangriffspunkt.

Durchgeführt

werden

neun

Rechnungen mit unterschiedlicher Elementdiskretisierung. Dabei ändert sich die Verteilung der Elemente im Verhältnis über der Länge und Breite.

Abb. 3-1: FE-Modell Test 1

Aufgrund der Eigenheiten bei der Schalendefinition von CalculiX (Expansion zu Hexaederelement mit 20 Knoten, siehe Kapitel 2.3.), verhält sich die Struktur bei einer geringen Elementanzahl über der Länge zu steif. Es treten deutliche LockingEffekte auf, die sich in einer zu kleinen Verschiebung des Lastangriffspunktes zeigen. In der Ergebnisauswertung (Abb. 3-3/Abb. 3-4) ist dieser Effekt deutlich zu beobachten. Das relativ große Locking, das sich bei der anfänglichen geringen Elementanzahl einstellt, lässt sich durch das schlechte Seitenverhältnis (engl.: „aspect ratio“) bei der Berechnung dünner Strukturen durch Volumenelemente erklären. Mit feiner werdender Diskretisierung über der Länge sinkt das Verhältnis von

h 5 Elementlänge h 20 = auf = und es stellt sich ein brauchbares Ergebnis ein. Elementdicke t 2 t 2 Fachhochschule Aachen


23


Die farbliche Auswertung der Verschiebung (Abb. 3-2) zeigt wie sich die Durchbiegung

durch

Steigerung

der

Elementanzahl

über der

Länge

zur

Einspannung hin verschiebt, d.h. die Struktur wird weicher.

Abb. 3-2: Einfluss Elementdiskretisierung auf Verschiebung (Bild: CalculiX-GraphiX)

Bei der Auswertung der Ergebnisse wird deutlich, dass das Lösungsverhalten maßgeblich von der Elementdiskretisierung über der Länge abhängt (Abb. 3-3). Mehr Elemente über die Breite verbessern die Ergebnisse kaum. Allgemein bleibt die Differenz, gegenüber dem analytischen Ergebnis zwischen ca. 4 - 6% für die grobe Diskretisierung und verringert sich auf 1,5% für die Feinste. Im Vergleich mit der NASTRAN-Lösung, sowohl mit als auch ohne Mittelknoten auf den Elementkanten, wird kein großer Unterschied zu den Ergebnissen der analytischen Lösung festgestellt (Abb. 3-4). Die Simulation durch NASTRAN liefert bei schon fünf Elementen über der Länge eine gute Lösung, die bei Erhöhung der Elementanzahl weiterhin relativ konstant bleibt. Eine Besonderheit von CalculiX ist, dass bei Simulationen von Schalen- und Balkenelementen die Berechnung durch einen iterativen Solver erfolgt. Der Gebrauch einer solchen Lösungsstrategie ist für diese Berechnung durchaus fragwürdig. Durch die kleine Krafteinleitung kommt es nur zu sehr geringen Verschiebungen im System. Hierdurch wird das Konvergenzkriterium schnell erfüllt



24


und die Rechnung wird zu frühzeitig abgebrochen, obwohl evtl. noch nicht die richtige Lösung erreicht wurde.

Differenz [%]

0 -2 -4 -6 -8 -10 -12 5x1

10x1

20x1

5x2

10x2

20x2

Diskretisierung (Länge x Breite )

5x3 CalculiX

10x3 20x3 OOFem

Abb. 3-3: Vergleich analytische Lösung / Open-Source

Die Qualität der OOFEM Simulation liegt deutlich unter der von CalculiX. Im Vergleich mit dem analytischen Ergebnis werden wie bei CalculiX zu geringe Werte berechnet (Abb. 3-3). Eine Verbesserung wird auch nicht durch eine Erhöhung der Elementanzahl erreicht. Der Grund hierfür liegt darin, dass OOFEM für die Simulation von dünnen Strukturen nur lineare dreieckige Plattenelemente ohne Knoten auf dem Elementrand anbietet. Die hier verwendete Elementanzahl ist eindeutig zu gering, um ein akzeptables Ergebnis zu erhalten. Es wäre ein deutlich feineres FiniteElemente Netz nötig, um den Fehler durch die zu steifen Biegeeigenschaften der linearen Dreieckselemente zu verringern. Nastran (lin) - OpenSource 0 Differenz [%]

-2 -4 -6 -8 -10 -12 5x1

10x1


20x1

5x2

10x2

20x2

5x3


10x3

20x3

25


Nastran (quad) - OpenSource

Differenz [%]

0 -2 -4 -6 -8 -10 -12 5x1

10x1

20x1

5x2

10x2

20x2

Diskretisierung (Länge x Breite)

5x3 CalculiX

10x3

20x3

OOFem

Abb. 3-4: Varianten NASTRAN (Schubweich) / Open-Source

3.1.1 Fazit Allgemein ist die Ergebnisqualität von CalculiX bei einer sinnvollen Diskretisierung zufrieden stellend. Die Elementabbildung von Schalenelementen, innerhalb von CalculiX, durch Hexaederelemente führt zu einem steifen Strukturverhalten (Abb. 3-5). Als kritischer Parameter, von dem das Konvergenzverhalten der Lösung abhängt,

hat

sich

das

Seitenverhältnis

der

Volumenelemente

größte Elementseitenlänge l = (engl. aspect ratio) heraus gestellt. Je feiner die Elementdicke t

Diskretisierung ist, desto kleiner wird das Seitenverhältnis. Die Elemente werden zunehmend weicher und die Lösung streb hin zum analytischen Ergebnis. Die restliche Abweichung des Ergebnisses ergibt sich durch die iterative Lösungsstrategie. Durch die Schalenelemente und die Randbedingung auf den rotatorischen Freiheitsgraden (siehe Kapitel 2.3) erfolgt die Berechnung iterativ. Die kleinen Verschiebungen im System führen zu einem schnellen erreichen des Konvergenzkriteriums und die Rechnung wird zu früh beendet. OOFEM ist für den professionellen Einsatz nicht geeignet. Das Fehlen von Rechteckplattenelementen bzw. Elementen mit Knoten auf den Elementkanten erweist sich doch als gravierender Nachteil, neben den schon in Kapitel 2.3.2 aufgeführten Einschränkungen. Es zeigt sich hier das deutlich schlechtere Verhalten von Dreiecksplattenelementen. Auch konnte der integrierte Postprocessor nicht zum Laufen gebracht werden, wodurch keine grafische Auswertung erfolgen konnte. Als Alternative bietet sich das Fachhochschule Aachen


26


Exportieren der Ergebnisse in das VTK-Format an. Aber diese Möglichkeit funktionierte nicht direkt und kann nur mit einigem Programmieraufwand erneut funktionsfähig gemacht werden.

Differenz [%]

0 -0,5 -1 -1,5 -2 20x1

20x2

20x3

Diskretisierung (LxB)

CalculiX Abaqus S8R

Nastran linear Abaqus S8R5

Nastran quadratisch

Abb. 3-5: Vergleich aller FE-Ergebnisse mit der analytischen Lösung.

In Abb. 3-5 sind die Ergebnisse der kommerziellen Programme für das Problem dargestellt. Auch hier lassen sich kleine Unterschiede zum analytischen Wert beobachten. Jedoch ist die Ergebnisqualität schon bei der groben Diskretisierung deutlich besser als bei den Open-Source Programmen.



27


3.2

Biegeeigenfrequenz Balken

Simuliert wird ein kurzer gedrungener Balken mit Rechteckquerschnitt, mit Hilfe von Balkenelementen. Der Balken ist einfach gelagert, was bedeutet, dass keine Zwangsbedingungen auf die rotatorischen Freiheitsgrade (Einspannung, Fixierung) angewendet werden. Eine Verschiebung in x-Richtung ist am Lager A möglich (Abb. 3-6). Die anderen translatorischen Freiheitsgrade der Lagerung sind fixiert. Berechnet werden die ersten drei Biegeeigenfrequenzen, bei unterschiedlichen Diskretisierungsgrad des Balkens.

Abb. 3-6: FE-Model Test 2

Die Elementanzahl der verschiedenen Eigenfrequenzsimulationen beträgt 3, 5, 6 und

30

Balkenelemente.

Dabei

ändert

sich

die

Elementlänge

von

l 30mm Gesamtlänge des Balkens l 30mm = = auf = . Anzahl Elemente n 3 n 30

Die Eigenformen der Biegemoden werden von CalculiX richtig berechnet. Die Darstellung durch den Postprocessor CalculiX-CrunchiX funktioniert problemlos (Abb. 3-7). Auffällig ist, dass der dritte Eigenmode mit nur fünf Balkenelementen richtig dargestellt wird. Dies erklärt sich dadurch, dass die Balkenelemente zu Fachhochschule Aachen


28


Hexaederelementen mit Mittelknoten (siehe Kapitel 2.3.) expandiert werden. Hierdurch wird eine Verformung der Balkenelementmitte möglich, während bei der Simulation mit NASTRAN und Abaqus lineare Balkenelemente verwendet werden.

I Mode / 5 Elm

II Mode / 5 Elm

III Mode / 5 Elm

I Mode / 30 Elm

II Mode / 30 Elm

III Mode / 30 Elm

Abb. 3-7: Eigenformen CalculiX (Bild: CalculiX-Graphix)

Im Vergleich mit der analytischen Lösung nach Bernoulli liefert CalculiX bei drei Elementen ein gutes Ergebnis für die ersten drei Eigenfrequenzen. Mit steigender Elementanzahl und höheren Eigenmoden wird die FE-Struktur zunehmend weicher (Abb. 3-8) und die FE-Lösung divergiert zunehmend zu niedrigen Eigenfrequenzen. Für die Eigenfrequenzen nach der Timoshenkotheorie fallen die Ergebnisse für die kleinen Diskretisierungsgrade bei steigendem Eigenmode schlechter aus. Eine Erhöhung der Elementanzahl verbessert die Werte für die Eigenfrequenzen und das Ergebnis konvergiert hin zu den analytischen Werten (Abb. 3-8). Der Vergleich mit NASTRAN ergibt ein anderes Bild der Lösungsqualität von CalculiX. Hier fällt auf, dass für drei Balkenelemente kein Ergebnis für den dritten Eigenmode vorliegt. NASTRAN liefert hierfür keinen Wert, was daran liegt, dass die Eigenform mit drei linearen Balkenelementen nicht abgebildet werden kann, während bei CalculiX durch die Expansion zu Hexaederelementen mit Mittelknoten eine Deformation in der Balkenmitte zulässig ist. Die Ergebnisse der beiden Simulation differieren in den Schranken von ±8% im Vergleich mit der NASTRANLösung und korrelieren deutlich besser als im Vergleich mit den analytischen Lösungen.



29


25 20 Differenz [%]

15 10 5 0 -5 -10 -15 -20 3

5

Mode 1 (Bernoulli) Mode 1 (Timoshenko)

6 Elementanzahl Mode 2 (Bernoulli) Mode 2 (Timoshenko)

30 Mode 3 (Bernoulli) Mode 3 (Timoshenko)

Abb. 3-8: Vergleich analytische Lösung (nach Timoshenko und Bernoulli) / CalculiX

Im Vergleich mit den NASTRAN-Bernoulli-Lösungen wird das Strukturverhalten zu weich

wiedergegeben.

Mit

zunehmender

Elementanzahl

und

Eigenmode

divergieren die Ergebnisse schnell (Abb. 3-9). Ein ähnliches Ergebnis lässt sich für die

NASTRAN-Timoshenko-Lösungen

feststellen,

wobei

jedoch

die

Eigenfrequenzen anfänglich zu hoch berechnet.

Nastran

8 6 Differenz [%]

4 2 0 -2 -4 -6 -8 3

5

Mode 1 (Bernoulli) Mode 1 (Timoshenko)

6 Elementanzahl Mode 2 (Bernoulli) Mode 2 (Timoshenko)

30 Mode 3 (Bernoulli) Mode 3 (Timoshenko)

Abb. 3-9: Vergleich NASTRAN / CalculiX



30


3.2.1 Fazit CalculiX berechnet die Eigenfrequenzen mit relativ guter Übereinstimmung zu den Werten von NASTRAN (Abb. 3-9). Es lässt sich feststellen, dass CalculiX immer niedrigere Werte im Vergleich mit NASTRAN liefert. Der verwendete Ansatz von CalculiX, Balkenelemente durch Hexaederelemente abzubilden, kann durchaus verwendet

werden

und

zeigt

bei

vernünftiger

Diskretisierung

eine

gute

Ergebnisqualität. Mit zunehmender Elementanzahl divergiert das Ergebnis jedoch deutlich zu niedrigeren Eigenfrequenzen, was darauf schließen lässt, dass das Strukturverhalten bei Simulationen mit Balkenelementen zu weich wiedergegeben wird. Im Vergleich mit der analytischen Lösung nach der Bernoulli-Theorie werden die Ergebnisse von CalculiX zu niedrig berechnet. Bezüglich der TimoshenkoTheorie sind die Werte anfangs jedoch zu hoch, während bei steigender Elementanzahl konvergiert das Ergebnis zu dieser analytischen Lösung konvergiert.



31


3.3

Eigenfrequenzen Platte

In der hier durchgeführten Simulation werden die ersten sechs Eigenfrequenzen einer dünnen quadratischen Platte berechnet (Abb. 3-10). Die Struktur ist an den Rändern

fest

eingespannt,

so

dass

in

den

sechs

Freiheitsgraden

der

Plattenelementrandknoten keine Verschiebung stattfinden kann.

Abb. 3-10: FE-Modell Test 3

Die Eigenformen der verschiedenen Eigenmoden werden von CalculiX-GraphiX (Postprocessor siehe Kapitel 2.3) richtig dargestellt. Hierfür wird die Möglichkeit Farbplots (engl. fringe plot) zu erstellen genutzt, wodurch die unterschiedlichen Deformationsgrade der einzelnen Elemente farblich getrennt aufgetragen werden. Im direkten Vergleich der Plots mit NASTRAN (diese wurden mit dem Postprocessor PATRAN erstellt) lässt sich jedoch eine kleine Drehung bei der Plattendeformation der Modes 2 und 6 um die Symmetrieachsen erkennen. Diese Abweichung ist kein Fehler in der grafischen Darstellung, sondern kann leicht durch nummerische Ungenauigkeiten bei der Berechnung (Nachkommastellenanzahl, 64/32 Bit Betriebsystem) oder Vernetzung (Knotenkoordinaten nicht exakt genug) entstehen. Fachhochschule Aachen


32


Abb. 3-11: Eigenformen CalculiX (Bild:CalculiX-CrunchiX)

Grundsätzlich wird das Schwingungsverhalten der Platte in guter Korrelation zu den Farbplots von NASTRAN dargestellt (vergleiche Abb. 3-11 / Abb. 3-12). Ein Nachteil bei der farblichen Darstellung von Ergebnissen mittels CalculiX-GraphiX ist, dass sich die Farbskalierung nicht verändern lässt. Es werden immer die in CalculiXGraphiX voreingestellten Farben für die Darstellung benutzt. Nur die dazugehörige Werteskalierung kann geändert werden. Generell ist der Postprocessor sehr intuitiv zu bedienen, und nach kurzer Einarbeitung kann der volle Funktionsumfang genutzt werden (Siehe Kapitel 2.3).

Abb. 3-12: Eigenformen NASTRAN (Bild: Patran)

Da durch die Expansion der CalculiX-Schalenelemente zu Volumenelementen immer ein Hexaederelementnetz berechnet wird und in den vorhergehenden Test eine

teilweise

gute


Übereinstimmung

der

Ergebnisse

mit


Schalen-

und 33


Balkenelementen festgestellt wurde, wird hier das Referenzergebnis durch eine Volumenstruktur ermittelt, um einen Vergleich mit der „speziellen Schalen/Hexaederelementdiskretisierung“ zu ermöglichen. Für diesen Vergleich wird die Platte mit derselben Elementanzahl von 100x100 Elementen über der Länge und Breite (der Schalenstrukturen) mit Hexaederelementen diskretisiert. Über der Dicke wird mit drei Elementen vernetzt. Die Berechnung der Eigenfrequenzen der Volumenstruktur erfolgt mittels Abaqus. Alle anderen Ergebnisse der kommerziellen Programme beziehen sich auf die Simulation mit Schalenelementen. Die mittels CalculiX berechneten Eigenfrequenzen weisen eine konstante Differenz zu der Lösung der Hexaederstruktur auf. Es lässt sich keine Änderung im Ergebnis feststellen. Da hier nur ein Finite-Elemente-Netz berechnet wird, kann keine Aussage über das Lösungsverhalten (Konvergenzverhalten, Lockingphänomene) gemacht werden. Das im zweiten Test festgestellte zu weiche Verhalten durch die Hexaederelementabbildung lässt sich hier wieder beobachten. Mit aber nur ca. -1,7% Differenz zur Volumenstruktur und 1,5% zu den NASTRAN-Lösungen wird hier eine sehr gute Ergebnisqualität erreicht. Hexaedermodel 1,5 1

Differenz [%]

0,5 0 -0,5 -1 -1,5 -2 1

2

3

4

5

6

Mode

Calculix Nastran quadratisch

Abaqus S8R5 Analytisch

Nastran linear

Abb. 3-13: Vergleich Eigenfrequenzen Hexaeder- / Schalenmodell

In Abb. 3-13 ist neben den Finite-Elemente-Ergebnissen eine analytische Lösung der Eigenfrequenzen für eine fest eingespannte Platte dargestellt (nach Young Fachhochschule Aachen


34


[LEI61]). Diese wurde im Rahmen einer Untersuchung der Nasa [LEI69] zum Schwingungsverhalten von Schalen- und Plattenstrukturen veröffentlicht. Im Anhang sind die Berechnungsformel und die Eigenformen aus der Untersuchung aufgeführt. Auf eine Bewertung der Eigenfrequenzwerte der analytischen Lösung bzw. der Finite-Elemente Simulation wird hier verzichtet. Eine Diskussion zum Thema ist im Anhang zu finden.

3.3.1 Fazit Das Ergebnis der Eigenfrequenzsimulation mittels CalculiX ist durchaus zufrieden stellend. Für die hier verwendete feine Vernetzung der Platte ergibt sich eine gute Lösungsqualität. Auch hier lässt sich ein kleiner Unterschied zu den Werten der kommerziellen

Programmen

feststellen,

der

durch

die

spezielle

Schalenelementabbildung hervorgerufen wird. Durch dieses Vorgehen wird das Strukturverhalten zu weich bei dynamischen Simulationen abgebildet. Bei weiterer Verfeinerung

des

Netzes

ist

davon

auszugehen,

dass,

wie

in

der

Eigenfrequenzberechnung der Balkenstruktur festgestellt wurde, eine weitere Steigerung der Genauigkeit erfolgt. Grundsätzlich ist jedoch eine gute Näherung zur Referenzlösung mit einer konstanten Differenz von ca. 1,8% gefunden worden. Das Auswerten

der

Ergebnisse

durch

CalculiX-GraphiX

erfolgt

nach

kurzer

Einarbeitungszeit ohne weitere Probleme. Die Oberfläche des Postprocessors ist sehr einfach gehalten und intuitiv zu bedienen. Alle wichtigen Funktionen werden mittels der rechten Maustaste aktiviert. Negativ aufgefallen ist nur, dass sich die Farben für die Darstellung der Ergebnisse nicht einstellen lassen. Hierfür werden immer Voreinstellungen verwendet und nur die dazugehörigen Skalierungswerte können geändert werden.



35


3.4

Test Modellgröße

Die Einsatzfähigkeit eines Finite-Elemente Programms wird hauptsächlich durch die Größe der rechenfähigen Finite-Elemente Netze bestimmt. Bei Simulationen in der Automobilindustrie werden Strukturen mit mehreren Millionen Elementen berechnet. In Optimierungsprozessen sind meist nur Teilssysteme mit einigen tausend Elementen zu simulieren, wobei es hier auf die Schnelligkeit der Lösung ankommt, um den Optimierungsprozess zeitlich nicht unnötig zu verlängern. Um die maximal mögliche

Modellgröße

der

simulationsfähigen

Finite-Elemente-Strukturen

zu

bestimmen, wird ein Quader mit Hexaederelementen vernetzt. Die Anzahl der Elemente wird für die Simulation sukzessive erhöht. Für die durchgeführten Rechnungen ergibt sich folgende Element-/Knotenanzahl (Tabelle 3-1).

Simulation Elementanzahl Knotenanzahl

I 8000 9261

II 64000 68921

III 216000 226981

IV 512000 531441

V 1020100 1050804

Tabelle 3-1: Modellgrößen des Benchmarks

Die

Volumenstruktur

ist

an

der

unteren

Fläche

fest

eingespannt.

Eine

Einheitsflächenlast wird auf der gegenüberliegenden Seite aufgebracht. Zur Durchführung der Simulationen werden ein 32Bit und ein 64 Bit System verwendet, die jeweils mit 2 Prozessoren ausgestattet sind. Der 32Bit Computer kann bei einer Simulationen pro CPU maximal 2GB Arbeitsspeicher (RAM) adressieren. Für das 64Bit System können höchstens 12GB Arbeitspeicher für die Rechnung allokiert werden, also pro CPU 6Gb RAM. CalculiX bietet für die Berechnung von Finite-Elemente Modellen drei verschiedene Lösungsverfahren an. Einen direkten Gleichungslöser und zwei iterative Verfahren. Bei

kleinen

Modellen

und

genügend

Arbeitsspeicher

kann

der

direkte

Lösungsalgorithmus „SPOOLES“ verwendet werden. Bei dem direkten Verfahren werden

die

Systemmatrizen

komplett

in

den

Arbeitsspeicher geschrieben

(gespeichert) und gelöst. Hierdurch wird die maximale Modellgröße durch den zur Verfügung

stehenden

Arbeitsspeicher

beschränkt.

Die

iterativen

Verfahren

benötigen deutlich weniger Arbeitsspeicher, da nur die aktuell zu berechnenden Matrixeinträge und Lösungen der vorher durchgeführten Operationen gespeichert werden. Vorzugsweise werden diese Verfahren bei Simulationen mit vielen Fachhochschule Aachen


36


Elementen oder bei nichtlinearen Problemstellungen verwendet. Hierfür bietet CaculiX die iterativen Verfahren „iterative Scaling“ und „iterative Cholesky“ an. In Abb. 3-14 ist der Vorteil der iterativen Verfahren in Bezug auf den benötigten Arbeitsspeicher deutlich zu erkennen. Der Speicherverbrauch der iterativen Solver ist bei dem hier simulierten Modell um den Faktor 4 kleiner als bei der Berechung mit dem direkten Solver „SPOOLES“. Ab einer Problemgröße von ca. 100 000 Elementen stellt sich auch ein Unterschied bei den iterativen Verfahren ein. 14 12

Speicher [Gb]

10 8 6 4 2 0 10000

100000

1000000

Elementanzahl SPOOLES

Iterative Scaling

Iterative Cholesky

Abb. 3-14: Speicherbedarf / Modellgröße (64Bit-System)

Hierfür ist die vorgeschaltete Konditionierungsstrategie (engl.: Preconditioner) verantwortlich. Das „iterative Scaling“ Verfahren bearbeitet nur die Diagonaleinträge der Matrix, wodurch weniger gespeichert werden muss. Der „iterative Cholesky“ Algorithmus bearbeitet die gesamten Matrixeinträge. Dies führt zu einem höheren Speicherverbrauch. Durch das Verbessern der gesamten Matrixstruktur arbeitet das „iteratve Cholesky“ Verfahren bei größeren Finite-Element Netzen effektiver, da durch die bessere Struktur der Matrix schneller Konvergenz erreicht wird (Abb. 3-15). Im Geschwindigkeitsvergleich arbeitet das direkte Verfahren weniger effizient als die iterativen Verfahren. Die Simulationszeit steigt doch deutlich mit zunehmender Modellgröße an. Im Vergleich ist die Rechenzeit bei nur 100000 Elementen circa um den Faktor 1000 langsamer.



37


4000

Rechenzeit [sec]

3000

2000

1000

0 10000

100000

1000000

Volumenelemente iter.Cholesky 32 Bit iter.Cholesky 64 Bit

iter.Scaling 32 Bit iter.Scaling 64 Bit

SPOOLES 32 Bit SPOOLES 64 Bit

Abb. 3-15: Rechenzeit: Solver / Modellgröße

Auch die iterativen Verfahren unterscheiden sich in der benötigten Rechenzeit, wie in Abbildung 3-15 zu beobachten ist. Bei großen Modellen ist die Rechenzeit beim „iterative Cholesky-Verfahren“ im Gegensatz zum iterativen Scaling-Verfahren geringer. Dieser Geschwindigkeitsvorteil wird durch die bessere Matrixstruktur hervor gerufen. Zwar braucht das Cholesky-Verfahren länger als das ScalingVerfahren, um die Struktur der Systemmatrizen zu verbessern, der anfängliche hohe Zeitaufwand wird jedoch dann durch das einfachere Berechnen der Gleichungen gut gemacht.

3.4.1 Fazit CalculiX bietet für die Berechnung der Systemmatrizen eine ausreichende Auswahl an Gleichungslösern. Mit dem direkten Verfahren lassen sich bei einem 32Bit System Modelle mit bis zu ca. 70000 Elementen berechnen. Bei 64Bit Systemen und 12GB Arbeitsspeicher können Simulationen mit ca. 200000 Elementen durchgeführt werden. Die Rechenzeit ist im Gegensatz zu den iterativen Verfahren deutlich

länger.

Dafür


sind

die

iterativen

Verfahren


effektiver

in

der 38


Speichernutzung. Bei vorsichtiger Abschätzung können hier Modelle mit bis zu 2000000 Elementen simuliert werden. Hierfür ist jedoch ein 64Bit Computersystem mit 12 GB Arbeitsspeicher nötig. Ein 32Bit System kann Modelle mit circa 200000 Elementen

berechnen.

Das

„iterative

Scaling“

Verfahren

hat

einen

Geschwindigkeitsvorteil gegenüber dem „iterative Cholesky“ Verfahren, welcher sich jedoch bei größeren Modellen schnell relativiert. In der Speichernutzung ist das „iterative Scaling“ Verfahren geringfügig effektiver. Grundsätzlich handelt es sich bei den hier angegebenen Modellgrößen, für die eine Berechnung noch möglich ist, um Schätzwerte. Bei dem hier verwendeten Testmodell handelt es sich um ein einfaches Hexaederelementmodell. Für komplexeren Modelle, z.B. Simulationen mit Zwangsbedingungen oder nichtlineare Simulationen, kann sich die maximale Elementgrenze verringern. Weiterhin kann sich die Berechung mittels des direkten Verfahrens, bedingt durch das schlechte Konvergenzverhalten von komplexen Strukturen, als vorteilhafter herausstellen und schneller zu einem Ergebnis führen. Nicht mehr getestet wurden die MPI (Message Passing Interface) Schnittstelle und ein „out-of-core“ Solver. Beide konnten aufgrund des erhöhten Zeitaufwandes beim Kompilieren von CalculiX nicht mehr in das Programm implementiert werden. Die MPI-Schnittstelle ermöglicht eine Simulation auf einem Clustersystem bzw. auf mehreren Computern parallel. Hierzu wird das Modell in kleinere Teile gesplittet und dann jeweils auf einem Rechner gelöst. Der „out-of-core“ Solver schreibt Teile der Systemmatrizen

auf

die

Festplatte,

wenn

der

zur

Verfügung

stehende

Arbeitsspeicher nicht mehr ausreicht. Grundsätzlich wird durch diese beiden Features die maximale Modellgröße erhöht.



39


3.5

Zusammenfassung

Durch die verschiedenen Simulationen sind einige Besonderheiten durch die speziellen Schalen- und Balkenelementdefinition innerhalb von CalculiX aufgefallen. Grundsätzlich können nur quadratische Schalen- und Balkenelemente verwendet werden, d.h. Elemente mit Knoten auf den Elementkanten, da die Elemente, wie in Kapitel 2.3 beschrieben zu Hexaederelementen mit 20 Knoten (Abb. 2-6/Abb. 2-7) expandiert werden. Durch kleine Veränderungen in Abaqusinputdateien lassen sich leicht Inputdecks für CalculiX erstellen. Hierbei sollte jedoch im Kapitel „Input Deck Format“ des CalculiX-Manuals nach Unterschieden und Einschränkungen geschaut werden. Bei statischen Simulationen verhalten sich die Schalenelemente bei zu groben Finite-Elemente Netzen zu steif. Durch Steigerung der Elementanzahl und das dadurch verbesserte Seitenverhältnis der Hexaederelemente konvergiert die Lösung hin zu den Ergebnissen von NASTRAN bzw. der analytischen Lösung. Die Eigenfrequenzsimulation für schubsteife Balkenelemente (Bernoulli-Theorie) mit CalculiX ergab zu niedrige Frequenzen im Vergleich mit der NASTRAN-Lösung und den

analytischen Werten.

Mit

steigendem

Vernetzungsgrad

und

höheren

Eigenmoden verschlechtern sich die Ergebnisse zunehmend. Für

schubweiche

Balkenelemente

(Timoshenko-Theorie)

werden

die

Eigenfrequenzen anfänglich zu hoch (zu große Steifigkeit) berechnet. Durch Verfeinerung des Netzes konvergiert das Ergebnis hin zu den analytisch berechneten

Werten.

Im

Vergleich

mit

NASTRAN

werden

jedoch

die

Eigenfrequenzen mit steigender Netzfeinheit zu niedrig berechnet. Bei der dynamischen Simulation mittels Schalenelemente wird eine zu weiche Steifigkeit

festgestellt.

Hierfür

wurde

als

Referenzlösung

eine

Platte

mit

Hexaederelementen vernetzt, um einen Vergleich der Lösungsqualität durch die speziellen Schalen-/Balkenelemente zu erhalten. Die berechneten Eigenfrequenzen der Platte durch CalculiX liegen leicht unter denen von NASTRAN und der Referenzlösung der Hexaederelementstruktur. Das Auswerten der Ergebnisse von CalculiX durch den mitgelieferten Pre/Postprocessor CalculiX-GraphiX erfolgt nach kurzer Einarbeitungszeit und unter Verwendung des gelungenen Manuals intuitiv. Der Funktionsumfang umfasst alle Fachhochschule Aachen


40


wichtigen Darstellungsmöglichkeiten von kommerziellen Programmen, wie z.B. Farbplots,

Animationen

oder

Erstellen

von

Bilddateien,

um

die

Simulationsergebnisse grafisch auszuwerten. CalculiX-GraphiX bietet ebenfalls grundlegende Funktionen für die Vernetzung von CAD-Daten, Konvertieren von Simulationsinputdecks und Erstellen von simplen Finite-Elemente-Netzen. Da verschiedene Simulationsinputdateien für die Tests benötigt werden, sind alle Strukturen mit dem kommerziellen Preprocessorprogramm ANSA vernetzt. Daher lässt sich keine Aussage über die Qualität der Vernetzung mittels CalculiX-GraphiX bzw. der Benutzerfreundlichkeit machen. Der Vorteil durch dieses Vorgehen liegt darin, nur ein Finite-Elemente-Netz zu erstellen und dieses in die verschiedenen Formate von NASTRAN, Abaqus und CalculiX zu überführen. CalculiX bietet eine gute Grundausstattung an direkten und iterativen Solvern. Die maximale Modellgröße für das direkte Verfahren liegt bei ca. 70 000 Elementen für das 32Bit System bzw. 200 000 Elementen für das 64Bit Computersystem. Mit einem 32Bit Computer können die iterativen Lösungsalgorithmen je nach Verfahren bis ca. 500 000 Elemente simulieren und bei dem hier verwendeten 64Bit Rechner bis ca. 2 000 000 Elementen, wobei nur bis 1 000 000 getestet wurde und die hier angegebenen Grenze eine vorsichtige Hochrechnung ist. Die durchgeführte Simulation mittels OOFEM weist sehr schlechte Werte auf. Im Vergleich zu CalculiX, NASTRAN und dem analytischen Ergebnis ist die berechnete Verschiebung deutlich zu klein und selbst bei Verfeinerung des Netzes wird keine Verbesserung

erreicht.

Hier

besteht

noch

Entwicklungsbedarf

bei

der

Elementbibliothek. Diese ist deutlich zu klein und bietet nur lineare dreieckige Schalenelemente für Strukturberechnungen. Das implementierte Element kann nicht in dynamischen Simulationen verwendet werden. Der Input der Simulationsdateien ist sehr kryptisch und es lassen sich vorhandene Finite-Elemente Netze nicht in OOFEM-Dateien konvertieren.



41

4 Optimierung

4 Optimierung Neben der Finite-Elemente-Methode hat sich die Optimierung als weiteres Werkzeug zur Lösung technischer Probleme etabliert. Soll zum Beispiel ein Bauteil möglichst leicht werden, aber gleichzeitig die Steifigkeit erhöht werden, nutzt man eine Optimierungsschleife, um diesen Zielkonflikt schnell und effizient zu lösen. Hierfür

bieten

die

Optimierungssoftwareprogramme

eine

Vielzahl

von

Optimierungsarten und –algorithmen an, um eine wirkungsvolle Durchführung, bei gegebener Rechnerkapazität und Problemstellung zu gewährleisten. Einen einführenden Überblick über das Thema der Optimierung bietet [BAI94]. Frei

verfügbare

Finite-Elemente

Programme

sind

für

den

Einsatz

in

Optimierungsanalysen sehr interessant. Hier ist es oft notwendig, viele Rechnungen parallel laufen zu lassen, um schnell zu einem Ergebnis zu gelangen. Dabei fallen schnell immense Lizenzkosten an, die durch den Einsatz von Open-SourceSoftware minimiert werden können. Um zu testen, ob sich CalculiX in einen Optimierungsprozess implementieren lässt, wird eine Sensitivitätsanalyse durchgeführt. Hierfür wird ein Querträger aus einem aktuellen Fahrzeugmodell untersucht. Die Optimierung wird von dem Open-SourceSoftwarepaket DAKOTA gesteuert.



42

4 Optimierung

4.1

Einführung Sensitivitätsanalyse / Robustheitsanalyse

Die Simulation von Bauteilen bzw. des Bauteilverhaltens mittels der FiniteElemente-Methode oder die Optimierung eines Designentwurfs wird meistens unter der Annahme deterministischer Eingangsgrößen durchgeführt. Dabei wird außer Acht gelassen, dass aufgrund verschiedener Einflüsse diese gar nicht genau eingestellt

werden

können.

Die

Einflüsse,

die

zur

Streuungen

in

den

Eingangsparametern führen, ergeben sich zum Beispiel aus Fertigungsprozessen (Material, Geometrie) oder durch äußere Parameter wie Temperatur und Belastungsschwankungen. Die Streuungen in den Eingangsparametern eines Bauteildesigns führen dazu, dass sich das Bauteilverhalten nicht mehr genau beschreiben lässt und somit die Zielgrößen bei der Bewertung des Systems erheblich streuen können. Halten sich die Streuungen in gewissen (tolerierbaren) Grenzen

bzw.

lässt

sich

das

Bauteilverhalten,

trotz

Variation

in

den

Eingangsgrößen, annähernd genau vorhersagen, spricht man von einem robusten Design (Abb. 4-1). Das bedeutet, dass kleine Änderungen im Input zu kleinen Änderungen in der Antwort führen.

Abb. 4-1: "Robust Design"

Robustheit beschreibt somit die Sensitivität der Zielgrößen gegenüber Variabilität in den Eingangsparametern. Somit lässt sich die Robustheitsoptimierung als die Suche nach dem optimalen Punkt bei gleichzeitig minimaler Streuung der Zielgröße Fachhochschule Aachen


43

4 Optimierung

beschreiben. Bei ungeschickter Wahl der Eingangsvariablen kann sich der Rechenaufwand für die Optimierung auf eine große Anzahl von Varianten vergrößern. Um die Anzahl der Eingangsvariablen zu minimieren, führt man im Vorfeld eine Sensitivitätsanalyse durch. Hierbei werden die Eingangsparameter auf die reduziert, die besonders großen Einfluss auf das Bauteilverhalten haben.

4.1.1 Optimierungsprozess Finite-Elemente Querlenker Zur Durchführung der Sensitivitätsanalyse wird ein Querlenker aus einem fertig vernetzten Fahrzeugmodell genutzt. Das Modell (Abb. 4-2) besteht aus 24110 Tetraederelementen mit 10 Knoten pro Element bzw. 41399 Knoten im gesamten Modell.

Abb. 4-2: Querlenkermodell

Als Material wird ein isotroper Stahl mit elastoplastischem Verhalten gewählt. Die Randbedingungen sind im Punkt •

A: Translation z fixiert (rote Markierung)

•

B: Translation y,z fixiert, Masterknoten der RBE-Spinne

•

C: Translation x,y,z fixiert, Masterknoten der RBE-Spinne Fachhochschule Aachen


44

4 Optimierung

Belastet wird die Struktur durch eine Kraft von 16000 N in X – Richtung im Punkt A. Im rot markierten Auswertungsbereich wird die plastische Dehnung der Elemente ausgegeben. Für die Steuerung der Optimierung und Definition der Variablen muss der FiniteElemente-Input des Querlenkers modifiziert werden. Die Optimierungsvariablen werden in der ASCII-Datei des Inputs durch geschweifte Klammern deklariert (siehe Abb. 4-3). Die so markierten Einträge werden dann während der Optimierung durch DAKOTA

mit

den

entsprechenden

Optimierungsparameter

ausgetauscht.

Abb. 4-3: Ausschnitt Inputdeck CalculiX

Als Parameter wird die Position des Kraftangriffpunktes in den drei Raumrichtungen variiert. Weiterhin wird das Materialverhalten durch den Elastizitätsmodul sowie der Spannungs-Dehnungskurvendefinition des „Hardening“-Eintrages gestreut. Die Grenzen der Optimierungsvariablen werden im Steuerungsfile von DAKOTA definiert. Hierfür müssen die entsprechenden Einträge in der ASCII-Datei eingefügt werden

(Abb.

4-4).

Die

diskreten

Werte

für

den

Spannungsoffset

(Parallelverschiebung der Kurven), die Variation des E-Moduls und die Änderung der Krafteinleitungsposition werden dann durch statistische Algorithmen des Optimierers generiert und in die CalculiX-Inputdatei neu eingefügt. Fachhochschule Aachen


45

4 Optimierung

Abb. 4-4: Dakotainput

Zum Schluss werden noch die Optimierungsstrategie und die Art des Outputs der Optimierung bestimmt. Abschließend werden die Dateien des Finite-Elemente Inputs für CalculiX und der Dakota-Steuerungsfile

in

ein

Queingsystem

importiert.

Das

Queingsystem

koordiniert die einzelnen Simulationen auf den verschiedenen Rechnern. Hierfür wurden vorher Computersysteme reserviert, auf denen CalculiX installiert wurde und ausschließlich die Optimierungsrechnungen laufen.

4.2

Fazit der Optimierung

Ziel dieser Optimierung ist es zu zeigen, dass sich CalculiX in eine Optimierung mittels DAKOTA integrieren lässt und die berechneten Ergebnisse sich weiter verarbeiten bzw. auswerten lassen. Die Integration von CalculiX in den Prozess und die Auswertung der berechneten Ergebnisse mittels Excel und Knut verlief ohne Schwierigkeiten. Auf eine explizite Ergebnissauswertung, dass heißt eine Bewertung der Sensitivitäten wird an dieser Stelle verzichtet. Für die Optimierung wurden insgesamt 2000 Simulationen auf bis zu 100 Computersystemen gleichzeitig durchgeführt. Die Optimierung dauerte 6,5 Stunden, Fachhochschule Aachen


46

4 Optimierung

wobei im Schnitt eine Rechnung ca. 5 Minuten benötigt. Die Installation von CalculiX auf den verschiedenen Computern verlief ohne Schwierigkeiten. Hierbei ist es von Vorteil gewesen, dass CalculiX auf der verwendeten Linuxdistribution neu kompiliert wurde und so alle benötigten Programme und Pfaddefinitionen auf allen Rechner

identisch

waren.

Die

Verwaltung

der

Rechnungen

durch

das

Queingsystem verlief dabei problemlos und es wurde kein Computerabsturz beobachtet. Wie die meisten Finite-Elemente-Programme verwendet auch CalculiX für die Definition des Inputs eine ASCII-Datei. Hierdurch ist es dem Optimierungsprogramm DAKOTA möglich, änderung im Input vorzunehmen bzw. die Optimierungsvariablen einzufügen. Die vorgenommenen Änderungen im Input lassen sich z.B. mit Excel visualisieren. In Abbildung 4-5 ist beispielhaft die durch DAKOTA vorgenommene Variation im Elastizitätsmodul in der Inputdatei aufgetragen. Es lässt sich hier sehr gut die Streuung in den einzelnen Samples beobachten. Variation E-Modul

E-Modul [N/mm²]

213000 212000 211000 210000 209000 208000 207000 Samples Abb. 4-5: Änderung E-Modul durch DAKOTA

Die Ergebnisse der einzelnen Simulationen der Optimierung können im ASCIIFormat heraus geschrieben werden und anschließend mit Programmen, wie z.B. Excel oder der speziell für DAKOTA entwickelten Tool „Knut“, ausgewertet werden. Standardmäßig wird für das Postprocessing einer DAKOTA Optimierung bei P+Z Knut verwendet. Knut ist eine eigene Entwicklung, die für die mathematische Manipulation und das Visualisieren der Ergebnisse auf Matlab zurückgreift. Hierfür müssen die Ergebnisse in einem speziellen zeilenorientierten ASCII Format vorliegen, um ausgewertet werden zu können. Die CalculiX-Ergebnisse lassen sich Fachhochschule Aachen


47

4 Optimierung

problemlos in das spezielle Format importieren bzw. zu einem Gesamtergebnisfile umschreiben.

Ein

typisches

Beispiel

für

einen

Auswerteplot

einer

Sensitivitätsanalyse ist der Abbildung 4-6 zu entnehmen.

max. Dehnung

Scatterplot

z-Position

Abb. 4-6: Scatterplot

In dem sogenannten Scatterplot lassen sich Abhängigkeiten von Variablen identifizieren. Hier ist zu erkennen, dass die maximale Dehnung von der Lastangriffsposition in der Z-Richtung abhängt. Wenn zwei Parameter unabhängig voneinander sind oder nur eine schwache Abhängigkeit zueinander aufweisen, wäre eine Punktewolke ohne Struktur oder eine willkürliche Verteilung der Punkte im Plot aufgetragen. Eine genaue Beschreibung der Inputparameter von DAKOTA und ein kompletter Überblick der Optimierungsalgorithmen würden den Umfang der Arbeit übersteigen. Verwiesen sei hierfür auf das Manual, das sich im Internet unter der Adresse http://www.cs.sandia.gov/DAKOTA befindet.



48

5 Steifigkeitsberechnung Motorhaube

5 Steifigkeitsberechnung Motorhaube Für einen abschließenden Test wird eine relativ große komplexe Schalenstruktur gewählt. Hierbei ist es von Interesse, ob Simulationen mit den meist sehr detaillierten Modellen aus der Automobilindustrie durchgeführt werden können. Um eine Aussage über die Qualität der Berechnung zu machen, wird dieselbe Struktur in Abaqus simuliert. Das Ergebnis wird hinsichtlich der maximalen Durchbiegung aufgrund eines Kopfaufpralls und der Rechenzeit verglichen. Bei der hier simulierten Motorhaube handelt es sich um ein vernetztes Fahrzeugmodell eines Ford Taurus, Baujahr 2001. Das Finite-Elemente-Modell ist ein

frei

verfügbares

Crashmodell,

was

im

Internet

(http://www.ncac.gwu.edu/vml/models.html) herunter geladen werden kann. Das amerikanische „FHWA/NHTSA National Crash Analysis Center“ in Ashburn, Virigina stellt verschiedene Fahrzeugmodelle und Crashbarrieren zum Download bereit. Die Modelle sind im LS-Dyna-Inputformat definiert. Durch einen kommerziellen Preprocessors können die Daten in das Abaqus-Format und somit auch in das CalculiX-Format konvertiert werden.



49


5.1

Kopfaufpralldefinition

Für die Simulation wird ein Kopfaufpralltest nach EEVC WG-17 bzw. EURO-NCAP durchgeführt. Der Kopfaufpralltest nach EEVC WG-17 fordert für einen DummyKopf eines erwachsenen Menschen einen resultierenden HIC-Wert von 1000 oder niedriger. Das Kopfmodell wird in einem Winkel von 65° mit 40km/h auf verschiedene Bereiche der Motorhaube geschossen (Abb. 5-1). Eine genaue Beschreibung der Testprozedur findet sich im Internet unter der Adresse http://www.EuroNCAP.com. Das Gesetz der Europäischen Union für die Berechung des

HIC-Wertes

und

weitere

Einzelheiten

der

EEVC

WG-17

können

http://ec.europa.eu/enterprise/automotive/directives/vehicles/dir2003_102_ce.htm entnommen werden. Der HIC-Wert errechnet sich aus der Kopfbeschleunigung, die über ein Zeitfenster von 3ms erreicht wird, nach der Formel 2 ,5

 t2 1  a HIC = (t 2 − t1 )  ∫ a r dt •  , mit a r = . t 2 − t1  g  t1 Die Aufschlagbereiche werden für die verschiedenen Fahrzeugtypen bzw. Formen individuell bestimmt, wobei die resultierende Kopfbeschleunigung bei modernen Fahrzeugstrukturen im Mittel bei ca. 11

m liegt, um den HIC-Wert von 1000 zu s2

unterschreiten.

Abb. 5-1: Kopfaufprall [www.euroncap.com]

5.1.1 Modellbeschreibung FE-Motorhaube Für die Simulation wird nur eine Hälfte der Motorhaube berechnet. Das Motorhaubenmodell hat 12250 Schalenelemente (Abb. 5-2) und es erfolgt eine statische Berechung. Ausgewertet wird die maximale Durchbiegung aufgrund einer Flächenlast.

Als

Randbedingungen


gilt

für

das

Gelenk


und

das 50


Motorhaubenschloss, dass in allen Richtungen keine Verschiebung möglich ist. Für die Symmetrieebene gilt keine Verschiebung in y und keine Rotation um y, z. Die Lasteinleitung erfolgt kreisflächig an einer sehr steifen Stelle der Motorhaube. Die Fläche

der

Lasteinleitung

ergibt

sich

aus

dem

Kopfdurchmesser

eines

Erwachsenenkopfes von 165mm. Mit einer mittleren Kopfbeschleunigung von 11

m s2

und einen Gewicht von 4,8 kg für den Erwachsenenkopf ergibt sich eine resultierende mittlere Kraft von F=52,8 N. Um die zu erwartenden hohen Verformungen mit plastischer Deformation abbilden zu können, wird ein Stahlmaterial mit elastoplastischem Verhalten gewählt. Hierdurch wird eine nichtlineare Simulation erforderlich.

Abb. 5-2: FE-Modell Motorhaube

5.2

Diskussion der Simulationsergebnisse

Die Simulation wird in der Automobilindustrie durch eine Kontaktsimulation durchgeführt. Die Kontaktdefinitionen von CalculiX befinden sich jedoch noch in einem Prereleasestadium. Um Probleme bei der Simulation zu vermeiden, wurde hier eine vereinfachte statische Simulation durchgeführt. Fachhochschule Aachen


51


CalculiX berechnet für die maximale Verschiebung in Z-Richtung im Lastbereich einen Wert von 53mm auf einem 32Bit System. Dasselbe Modell auf einem 64Bit System ergibt ebenfalls einen Wert von 53mm. Hier kann sehr schön der Unterschied, der durch die maximal mitgeführten Nachkommastellen bei 64Bit Computersystemen

entsteht,

beobachten

werden.

Die

Anzahl

der

Nachkommastellen die für eine Berechnung verwendet werden können, ist bei 64Bit um den Faktor 1010 größer. Daher wird das Konvergenzkriterium aufgrund der höheren Genauigkeit später erreicht. Eine Kontrollrechnung mit Abaqus ergibt bei gleicher Modellierung eine maximale Verschiebung von 61mm.

Abb. 5-3: Verschiebung Motorhaube

Bei realen Motorhaubenstrukturen werden heute im Mittel circa 60mm durch den Aufschlag des erwachsenen Kopfes erreicht. Somit würde das mit CalculiX ermittelte Ergebnis zu einer falschen Bewertung der Motorhaube führen. Das zu steife Verhalten der Struktur lässt sich, wie in Kapitel 3.1 erwähnt, auf das schlechte Seitenverhältnis der Elemente zurückführen. Hier kann eine weitere Netzverfeinerung das Ergebnis positiv beeinflussen. Weiterhin werden, wie in Kapitel 2.3.1 aufgeführt, durch die gebogene Motorhaubenstruktur und das Festhalten der rotatorischen Freiheitsgrade an der Symmetrieebene „knots“, also Zwangsbedingungen eingeführt. Diese behindern die Struktur an der freien Durchbiegung. Fachhochschule Aachen


52


CalculiX benötigt für die Simulation ca. 30min., während Abaqus in ca. 2min. das Ergebnis berechnete. Der Rechenzeitunterschied lässt sich auf die zusätzliche Einführung der Zwangsbedingungen und die Expansion der Elemente zu Hexaedern zurückführen. Auch muss ein deutlich größeres Modell berechnet werden, da sich die Knotenanzahl durch die neuen Knoten der Hexaederelemente vergrößert. Die grafische Darstellung durch CalculiX-GraphiX lässt eine gute Korrelation

zu

dem

Abaqusergebnis

erkennen

(Abb.

5-3:

Verschiebung

Motorhaube). Hier ist leider, wie schon erwähnt, von Nachteil, dass sich die Farbverläufe in CalculiX-GraphiX nicht einstellen lassen. Allgemein ist das Ergebnis des Tests durchaus zufrieden stellend. Wenn bei komplexeren Schalenstrukturen bedacht wird, dass sich die Modellgröße erhöht und auf eine feine Diskretisierung geachtet wird, können die Simulationsergebnisse mit kommerziellen Programmen mithalten. Hier begrenzt jedoch der Arbeitsspeicher die maximale Modellgröße.



53

6 Diskussion der Ergebnisse

6 Diskussion der Ergebnisse Das Programm OOFEM ist für einen Einsatz in Optimierungsprozessen bzw. zur Berechnung komplexer Strukturen nicht geeignet. Es bietet nur eine beschränkte Elementbibliothek mit linearen Elementtypen an, also jeweils kein Knoten auf der Elementkantenmitte. Das Fehlen eines lauffähigen Pre-/Postprocessors und die sehr komplexe Inputdefinition erschweren ein Einarbeiten deutlich. Für den ersten Test, der statischen Balkenverformung, errechnet OOFEM viel zu niedrige Werte. Durch die lineare Elementabbildung sind hier deutlich mehr Elemente

erforderlich,

um

ein

zufriedenstellendes

Ergebnis

zu

erhalten.

Dynamische Berechnungen (Test 2/3) sind im jetzigen Entwicklungsstand und mit Balken- und Schalenelementen nicht möglich.

Das Softwarepaket CalculiX, mit dem Solver CalculiX-CrunchiX und dem Pre/Postprocessor CalculiX-GraphiX, lässt sich ohne Schwierigkeiten in Optimierungsprozesse einbinden. Auch können komplexe Strukturen, nichtlineare und Kontakt Berechnungen durchgeführt werden. CalculiX verwendet, wie in Kapitel 2.3.1 beschrieben, eine eigene Strategie, um Berechnungen mit Schalen- und Balkenelementen durchzuführen. Die statische Simulation ergibt ein gutes Ergebnis für die statische Balkenverformung und somit für die Steifigkeit der Schalenelemente. Für das feinste Modell beträgt die Differenz 1,8%. Der Vergleich mit der Lösung von NASTRAN ergibt eine Differenz von ca. 1,5%. Allgemein korreliert das Ergebnis jedoch gut. Als kritischer Parameter für die Ergebnisqualität konnte die Elementkantenlänge in Balkenlängsrichtung identifiziert werden. Eine Verfeinerung des Netzes über der Breite beeinflusst das Ergebnis nicht.

Durch den Test der Eigenfrequenzen einer Balkenstruktur wird festgestellt, dass durch die Abbildung der Elemente als Hexaeder der Schubeinfluss mit berücksichtigt wird. Also werden die Frequenzen in Analogie zur TimoshenkoTheorie berechnet. Im Vergleich mit der analytischen Lösung nach Bernoulli und der NASTRAN-Simulation mit schubsteifen Balkenelementen ergeben sich zu niedrige



54


Werte. Die Differenz steigt mit der Anzahl an Balkenelemente und höherem Eigenmode deutlich und das Strukturverhalten wird zu weich abgebildet. Für die schubweichen Ergebnisse korrelieren die Lösungen besser. Bei 30 Elementen beträgt die Differenz zur analytischen Lösung kaum 0,5%. Im Vergleich mit NASTRAN werden die Frequenzen jedoch mit steigendem Eigenmode und höherer Elementanzahl zu niedrig berechnet. Ein Vergleich der von NASTRAN berechneten Ergebnisse zeigte auch eine Differenz zu den analytischen Werten.

Die

Eigenfrequenzsimulation

der

quadratischen

Platte

ergibt

eine

gute

Übereinstimmung zu den Werten der Lösung von NASTRAN. Hier liegen die Werte von CalculiX wie bei den vorherigen Tests um die 1,5% unter denen der Referenzlösung von NASTRAN. Eine weitere Steigerung der Lösungsqualität lässt sich durch die Verfeinerung des FE-Netzes erreichen.

Der Test mit einer realen Motorhaubenstruktur (statische Analyse) ergibt ein anderes Bild der Lösungsqualität von CalculiX. Bei der komplexen Struktur wird die Durchbiegung im Gegensatz zu der Referenzlösung von Abaqus zu gering berechnet. Auch ist die Rechenzeit um ein Vielfaches länger. Durch die Expansion wird ein großer Teil der Rechenzeit für den Aufbau der Systemmatrizen benötigt. Das abweichende Ergebnis ist durch die schlechte Elementqualität zu erklären. Für die Berechnung ist ein feineres Netz erforderlich, um das schlechte Verhalten der Hexaederelemente bei komplexen Strukturen zu kompensieren.

Grundsätzlich sind die Ergebnisse aller Tests zufrieden stellend im Hinblick auf die Ergebnisqualität. Im Vergleich mit NASTRAN und der analytischen Lösung lässt sich bei einer der Problemstellung angemessenen Vernetzung eine gute Korrelation erreichen. Die Untersuchung der in CalculiX implementierten Solver ergab maximale Modellgrößen, die noch zu berechnen sind, von ca. 70000 Elementen für den direkten Solver und ca. 200000 Elementen für die iterativen Verfahren bei einem 32 Bit Computersystem. Auf 64 Bit Systemen verschieben sich die Grenzen hin zu ca. 200000 Elementen für den direkten Gleichungslöser und ca. 2000000 Elementen für die iterativen Solver. Die hier ermittelten Modellgrößen sind für eine Open-Source Fachhochschule Aachen


55


Software sehr gut. Der Haupteinflussfaktor für die Modellgrößen ist der verfügbare Arbeitsspeicher. Durch das Implementieren eines „out-of-core“ Solvers oder der MPI-Schnittstelle in CalculiX lassen sich diese Grenzen weiter nach oben verschieben. Problematisch hat sich die Expansion der Elemente zu Hexaedern heraus gestellt. Hierdurch ist das tatsächlich berechnete Modell deutlich größer als die des Schalenmodells. Das Implementieren von CalculiX in die Optimierungsschleife verlief ohne Schwierigkeiten. Die Lösungsdauer für die Optimierung ist mit ca. 6,5 Stunden zufriedenstellend. Das relativ große FE-Modell wurde im Schnitt in ca. 5 Minuten gelöst. Die Ergebnisse konnten mittels Excel und Knut weiter verarbeitet werden.

Abschließend kann gesagt werden, dass die Open-Source Software CalculiX für einen Einsatz in Optimierungsprozessen geeignet ist. Durch das zu Abaqus kompatible Inputformat erfolgt eine Einarbeitung relativ schnell. Ein weiterer Vorteil ist, dass kommerzielle Preprocessoren für die Vernetzung verwendet werden können. Die Ergebnisqualität von CalculiX ist zufriedenstellend. Bei komplexen Schalenelementstrukturen verhalten sich jedoch die Elemente von CalculiX zu steif im Vergleich zu „richtigen“ Schalenelementen. Der Fehler kann jedoch durch eine feinere Vernetzung korrigiert werden. Mit den implementierten Gleichungslösern können auf einem 64Bit Computer Modelle mit bis zu 2000000 Elementen simuliert werden.

Der

mitgelieferte

Postprocessor

bietet

grundlegende

Auswertemöglichkeiten und ist intuitiv zu bedienen. Das Integrieren in einen Optimierungsprozess mit DAKOTA verläuft aufgrund des im ASCII-Format definierten Inputs ohne Probleme.



56

A1 Anhang

A1 Open-Source FE-Programme (Stand Juni 2007) ALADDIN:

DOLFIN:

http://www.isr.umd.edu/~austin/aladdin.htm

http://www.nada.kth.se/~jhoffman/software/i

l

ndex.html

AxisVM: http://www.axisvm.com

DUNS: http://sourceforge.net/projects/duns

ADVENTURE:

DIANA:

http://adventure.q.t.u-tokyo.ac.jp/software

http://www2.tnodiana.com/gs/rdas/papp.asp?

AFEMS:

cmd=Applications

http://www.femengineering.com/index2.ht

DOUG:

ml

http://www.maths.bath.ac.uk/~parsoft/doug

ALBERTA: http://www.albertafem.de

DOLFIN: http://www.fenics.org/projects

Ansoft: http://www.ansoft.com

EMAP:

AnyBody: http://www.anybodytech.com

http://www.emclab.umr.edu/emap.html

BASIN:

EXPDE:

mailto:[email protected]

http://www10.informatik.unierlangen.de/~pfl

CAPCAST: http://www.ekkinc.com

aum/expde/public_html

Channelflow:

EULER:

http://www.cns.gatech.edu/channelflow

http://www.iee.cas.cz/staff/solin/euler

CLAWPACK:

ELFEN: http://rsazure.swan.ac.uk/elfen.htm

http://www.amath.-washington.edu/~claw

FeatFlow: http://www.featflow.de

COSMOS: http://www.cosmosm.com

FeaTure:

CURLY3D: http://curly3d.sourceforge.net

http://www.tm.ctw.utwente.nl/onderzoek/Di

COSHELL:

ekA/Feature

http://www.digitaladdis.com/sk/programs.h

FEC: ftp://ftp.math.uh.edu/pub/Math

tm

FELIB:

CodeAster: http://www.codeaster.org

http://www.softeng.cse.-clrc.ac.uk/felib4

COMSOL: http://www.femlab.com

FemObject:

DANSOFT: mailto:[email protected]

http://www.zace.com/femobj.htm

Deal II: http://www.dealii.org

FieldMagic: http://www.intactsolutions.com



57

A1 Anhang

FEMLAB:

Haplo: http://haplo.info

http://homepage.usask.ca/~ijm451/finite/fe

HMD: http://www.heldeneng.com

_resources/node164.html

Impact: http://impact.sourceforge.net

FE/Pipe: http://www.paulin.com

ISAAC: http://isaac-cfd.sourceforge.net

FELyX: http://felyx.sourceforge.net

Infolytica: http://www.infolytica.com

FreeFem: http://www.freefem.org/ff++

IMTEK:

FEMLIB:

http://www.imtek.unifreiburg.de/simulation/


mathematica/IMSweb


KASKADE:

Femlisp: http://www.femlisp.org

http://homepage.usask.ca/~ijm451/finite/fe_r

FELIPE: http://www.felipe.co.uk

esources/node181.html

FEMM: http://femm.fostermiller.com

KeyFE2:

FEMOctave:

http://users.skynet.be/keyFE2/keyFE2.html


KFEM: http://kfem.sourceforge.net


LUGR:

FEMSET:

http://homepages.cwi.nl/~gollum/LUGR/ind

http://www.hawhamburg.de/rzbt/dnksoft/ca

ex.html

mmpus/femset.html

LISA: http://www.lisa-fet.com/index640.htm

Free Finite Element Package:

LUSAS: http://www.lusas.com

http://www.free-finite-

MODFE:

elementpackage.smial.de

http://water.usgs.gov/software/modfe.html

FEAP:

MOUSE:

http://www.ce.berkeley.edu/~rlt/feap

http://www.vug.uniduisburg.de/MOUSE

Felt: http://felt.sourceforge.net

MOLDFLOW:

FreeFEM, FreeFEM3D:

http://www.moldflow.com/stp

http://www.freefem.org

MiniFem:

Geocrack:

http://www.rwthaachen.de/mmw/minifem/m

http://ww2.mne.ksu.edu/~geocrack

inifem.html

Gerris: http://gfs.sourceforge.net

MEANS: http://www.femcad.de

GeoFEM: http://geofem.tokyo.rist.or.jp

MGGHAT:

GEODE: http://www.geodesol.com

http://homepage.usask.ca/~ijm451/finite/fe_r

Getfem++:

esources/node190.html

http://www-gmm.insatoulouse.fr/getfem Fachhochschule Aachen


58

A1 Anhang

NASTRAN (NASA):

PHOENICS: http://www.cham.co.uk

http://www.openchannelsoftware.org/projec

PSU: http://www.isd.uni-stuttgart.-de/psu

ts/NASTRAN

Rheolef:

NaSt3DGP: http://wissrech.iam.uni-

http://www.sai.msu.su/sal/B/2/RHEOLEF.ht

bonn.de/research/projects/Na-St3DGP

ml

NIKE3D:

STAAD: http://www.reiusa.com

http://wwweng.llnl.gov/mdg/mdg_codes_ni

Solvia: http://www.solvia.se

ke3d.html

Straus 7: http://www.straus7.com

NIST:

STARS:


http://www.openchannelfoundation.org/order


s/index.php?group_id=191

NLFET: http://nlfet.sourceforge.net

STRAP: http://www.atirsoft.com

OFELI: http://ofeli.sourceforge.net

SLFFEA:

OpenSees: http://opensees.berkeley.edu

http://www.geocities.com/Athens/2099/slffe

Open Foam:

a.html

http://www.opencfd.co.uk/openfoam

TOCHNOG: http://tochnog.sourceforge.net

Open Flower:

VAMUCH:

http://openflower.sourceforge.net

http://www.mae.usu.edu/faculty/wenbin/ht_d

OpenFem:

ocs/vamuch.html

http://wwwrocq.inria.fr/OpenFEM

VAPAS:

ParaFEM:

http://www.mae.usu.edu/faculty/wenbin/ht_d

http://www.sve.man.ac.uk/Research/AtoZ/

ocs/vapas.html

ParaFEM

Warp3D:

PARSYS:

http://cern49.cee.uiuc.edu/cfm/warp3d.html

http://www.iee.cas.cz/staff/solin/parsys

Zebulon:

PHYSICA:

http://www.nwnumerics.com/Zebulon

http://www.sofistik.de/loesungen/statik-

Z88: http://www.z88.org

fem/multiphysics-cfd PolyDE: http://www.tuharburg.de/mst/english/polyd e/PolyDE.html



59

A2 Anhang

A2 Inputdecks Für die einzelnen Testbeispiele ist im folgenden Abschnitt beispielhaft jeweils ein Inputdeck aufgeführt.

A2.1 Test 1 Calculix *HEADING ** Test 1 20x3 Elemente ** ** Knotenkoordinaten ** *NODE,NSET=ALL 1,

0.,

100.,

0.

……. 653,

6.6666675,

95.,

0.

654,

8.333334,

95.,

0.

655,

6.6666675,

97.5,

0.

** ** Elementdefinition, quadratisches Plattenelement mit 8-Knoten mit voller Integration ** *ELEMENT, TYPE=S8, ELSET=PSHELL 1,

17,

464,

523,

512,

483,

522,

521,

515

2,

512,

523,

527,

511,

521,

525,

524,

514

……. ** ** Zuweisen der Eigenschaften der Plattenelemente ** *SHELL SECTION, ELSET=PSHELL, MATERIAL=STEEL 2., ** ** Materialparameter ** *MATERIAL, NAME=STEEL *DENSITY 7.85e-09,



60

A2 Anhang

*ELASTIC 210000., 0.3 ** ** Definition der Rechnung, Statik ** *STEP *STATIC ** ** Lastendefinition, Einheitskraft in neg. Z-Richtung ** *CLOAD 17,

3,

-1.

** ** Randbedingung, fest Eingespannt in allen Raumrichtungen ** *BOUNDARY, TYPE=DISPLACEMENT 1,

1,

6,

0.

516,

1,

6,

0.

517,

1,

6,

0.

518,

1,

6,

0.

519,

1,

6,

0.

520,

1,

6,

0.

3,

1,

6,

0.

** ** Definition des Outputs ** ** Verschiebung und Reaktionskräfte der Knoten *NODE PRINT, NSET=ALL U, RF ** ** Spannungen und Dehnungen für die Elemente *ELEMENT PRINT, ELSET=PSHELL S,E *N FILE U, RF *EL FILE S,E *END STEP



61

A2 Anhang

OOFem # Name Outputfile test1_20x3.out # Testbeschreibung job description: small beam test # Analyseart, Statik LinearStatic nsteps 1 # OOFem spezifische Parameter (siehe Dokumentation) domain 3dShell OutputManager tstep_all dofman_all element_all # Knotenanzahl|Elementanzahl|Querschnitt|Material|Randbedingungen| ndofman 84 nelem 120 ncrossSect 1 nmat 1 nbc 2 # Knotendefinition mit Lasteinleitung und Randbedingungen nodes

1 coords 3

100.

0. 0.0 bc 6 1 1 0 1 0 1 load 1 2

nodes

2 coords 3

100.

3.33333 0.0 bc 6 1 1 0 1 0 1

nodes

3 coords 3

100.

6.66667 0.0 bc 6 1 1 0 1 0 1

……. # Elementdefinition (Dreiecksplatte) mit Querschnitts- ,Material- und Intergrationspunktdefinition rershell

1 nodes 3

1

2

6 crossSect 1 mat 1 NIP 1

rershell

2 nodes 3

6

5

1

crossSect 1 mat 1 NIP 1

……. # Querschnitsparameter SimpleCS 1 thick 2.0 # Materialdaten IsoLE 1 d 7.85e-9 E 2.1e5 n 0.3 talpha 0.0 # Randbedingungen mit Lastaufbringung BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0 NodalLoad 2 loadTimeFunction 1 Components 6 0.0 0.0 1.0 0.0 0.0 0.0 ConstantFunction 1 f(t) 1.0



62

A2 Anhang

A2.2 Test 2 Calculix *HEADING ** *NODE 1,

0.,

0.

2,

0.5,

0.

……. ** Elementdefinition, Balken mit Zwischenknoten *ELEMENT, TYPE=B32 , ELSET=BEAM 1,

1,

2,

3

2,

3,

4,

5

……. ** ** Balkeneigenschaften *BEAM SECTION, ELSET=BEAM, SECTION=RECT, MATERIAL=STEEL 3., 5. 0.,0.,-1. ** *MATERIAL, NAME=STEEL *DENSITY 7.8e-09, *ELASTIC 210000., 0.3 ** ** Eigenfrequenzanalyse, ersten 20 Frequenzen warden berechnet ** *STEP *FREQUENCY 20, ** ** ** Output wird für jede Frequenz rausgeschrieben **EL PRINT,FREQUENCY=1 **NODE PRINT, FREQUENCY=1 *NODE FILE U



63

A2 Anhang

** Zur Darstellung des Querschnittes wird der Output in 3-Dimensionaler Form angefordert *EL FILE,OUTPUT=3D S *BOUNDARY 1,2,3 61,1,3 *END STEP

OOFem Eine Berechnung des Problems mit OOFem ist nicht möglich.



64

A2 Anhang

A2.3 Test 3 Calculix ** *HEADING

** ** ** *NODE 1,

0.,

100.,

28,

100.,

100.,

55,

100.,

0.,

82,

0.,

0.,

0. 0. 0. 0.

….. ** ** Elementdefinition, quadratisches Plattenelement mit 8-Knoten mit reduzierter Integration ** *ELEMENT, TYPE=S8R, ELSET= PSHELL 1,

1443,

1135,

82,

1238,

1442,

1237,

1339,

1441

2,

1134,

1135,

1443,

1448,

1236,

1442,

1446,

1444

3,

1238,

1239,

1449,

1443,

1340,

1445,

1447,

1441

……. ** ** ** SECTION DATA ** *SHELL SECTION, ELSET=PSHELL, MATERIAL=STEEL 2., ** ** *MATERIAL, NAME=STEEL *DENSITY 2.7E-9, *ELASTIC, TYPE=ISOTROPIC 70000.,

0.3

** ** Setdefinition um Randbedingungen komfortabler aufzubringen



65

A2 Anhang

** *NSET, NSET=BC 1,

28,

55,

82,

629,

630,

631,

632,

633,

634,

……. ** ** *STEP *FREQUENCY 20 ** ** BOUNDARY ** *BOUNDARY, TYPE=DISPLACEMENT BC,

1,

6,

0.

*NODE FILE U, *END STEP

OOFem Eine Berechnung ist mit der vorliegenden Version nicht möglich.



66

A3 Anhang

A3 Analytische Herleitung An dieser Stelle erfolgt eine Herleitung analytischer Ergebnisse für die Testbeispiele. Für das erste Testbeispiel wird die Durchbiegung infolge der Belastung einer exzentrisch angreifenden Kraft berechnet. Die Gesamtverschiebung des flachen langen Kragträgers ergibt sich aus Biege-, Schub- und Torsionseinfluss.

Für den zweiten Test werden die ersten drei Biegeeigenfrequenzen eines kurzen dicken

Balkens

ermittelt.

Dabei

findet

eine

Erweiterung

der

klassischen

Bernoullitheorie um die Wirkung Schub- und Rotationsträgheitseinfluss nach Timoshenko statt. Es werden die Grenzen der einzelnen Theorien gezeigt und weitere aus der Literatur entnommene Formeln für die Bestimmung der Biegeeigenfrequenzen aufgeführt.

Das dritte Beispiel behandelt die Eigenschwingungen von quadratischen, fest eingespannten Platten. Für diesen Fall werden die Ergebnisse aus der Literaturrecherche aufgeführt, die Bewegungsgleichung für die Platte in Analogie zur

Timoshenko-Theorie

angegeben

und

die

Grenzen

der

klassischen

Plattentheorie für den Einsatz in der Dynamik gezeigt.



67

A3 Anhang

A3.1 Test 1. Durchbiegung Blattfeder Im Test 1 wird die Durchbiegung eines schlanken Balkens (t << b) aufgrund einer exzentrisch angreifenden Kraft berechnet (Abb. A3-1). Dabei wird die Absenkung im Lastangriffpunkt betrachtet, die sich aus den Anteilen Biegung, Schub sowie Torsion ergibt.

x E,G

F t

l

b

Abb. A3-1: Blattfeder

Den Biegeanteil erhält man aus den Gleichgewichtsbedingungen für den Balken Q′ = − q , M ′ = Q ,

dem Elastizitätsgesetz für das Biegemoment M = EI ψ′

und bei Annahme, dass die Schubsteifigkeit sehr groß ist ( χGA → ∞ ), wird aus dem Elastizitätsgesetz für die Querkraft Q = χGA( w′ + ψ) w′ + ψ = 0 w′′ = −ψ′ .

durch Eliminieren von ψ′ die Differentialgleichung der Biegelinie zu w′′ = −

M . EI

Für das Moment gilt M = − F (l − x) Einsetzen und Integrieren führt zu EIw′′ = F (− x + l )  − x2  EIw′ = F  + lx  + C1  2   x 3 lx 2  EIw = F  − +  + C1 x + C2 2   6 Fachhochschule Aachen


68

A3 Anhang

Mit den Randbedingungen w′(0) = 0 , w(0) = 0 folgt für die Integrationskonstanten C1 = 0 , C2 = 0 . Damit wird der Biegeverlauf zu

Fl 3  x3 x2  − + 3  . l2  6 EI  l 3

w( x) =

Die maximale Durchbiegung tritt an der Krafteinleitungsstelle x = l auf: wmax = f B =

Fl 3 . 3EI

Bei der bis hier gemachten Herleitung wird davon ausgegangen, dass der Balken schubstarr ist, was bedeutet, dass der Querkrafteinfluss keine Winkeländerung hervorruft und Balkenquerschnitte, die vor der Deformation senkrecht zur neutralen Faser standen, auch nach der Deformation senkrecht zu ihr stehen bleiben (Abb. A3-2).

z

x

w

w′

−ψ

Abb. A3-2: Bernoulli Annahme

Als nächstes wird der Schubeinfluss auf den Biegebalken betrachtet. Das Elastizitätsgesetz für die Querkraft lautet w′ + ψ =

Q GAS

Mit Einführung eines Schubkorrekturfaktors für einen Rechteckquerschnitt χ =

5 , 6

AS = χA und wB′ = −ψ , ergibt sich für die Gesamtneigung (Indizies B: Biegung und

S: Schub) w′ = w′B + wS′

mit wS′ =

Q . Durch Integration und unter Beachtung der Randbedingungen GAS

wS (0) = 0 wird nun der Biegeverlauf des Balkens aufgrund des Schubeinflusses zu Fachhochschule Aachen


69

A3 Anhang

wS =

Q x. GAS

Die größte Durchbiegung des Balken tritt an der Kraftangriffstelle x = l auf und mit Q = F wird die Gesamtdurchbiegung aufgrund von Biege- und Schubeinfluss zu

Fl 3 Fl w= + . 3EI GAS Für den Torsionseinfluss auf die Balkenbiegung werden folgende Annahmen getroffen [GRO99]: •

Querschnitte verdrehen sich unter Torsionseinfluss als Ganzes, d.h. Punkte des Querschnitts, die vor der Verformung auf einer Geraden liegen, sind auch nach der Verformung auf einer Geraden.

•

Ebene Querschnitte bleiben eben. Es tritt keine Verformung aus der Ebene heraus auf (keine Verwölbung). γ

dx

dυ

b

Abb. A3-3: Zusammenhang Verdrehung und Winkeländerung

Bei kleinen Verformungen gilt für die Verdrehung d υ und die Winkeländerung γ der folgende Zusammenhang (Abb. A3-3) b b dυ d υ = γdx → γ = . 2 2 dx

Mit dem Elastizitätsgesetz τ = G γ folgt dann τ=G

b dυ b = G υ′ . 2 dx 2

Für das Torsionsmoment gilt b M T = ∫ τdA 2 M T = Gυ′∫

b2 dA 4

M T = Gυ′I P .



70

A3 Anhang

Dabei wird die neu eingeführte Größe I P als das polare Flächenträgheitsmoment bezeichnet, das im weiteren Verlauf als Torsionsträgheitsmoment IT bezeichnet wird. Durch die exzentrisch angreifende Kraft wirkt ein konstantes Torsionsmoment MT = F

b . Mittels Integration über die Balkenlänge von 2

υ′ =

MT GIT

ergibt sich die Verdrehung am Balkenende aufgrund des Torsionsmomentes zu υ=

MT l b und die Verschiebung zu wT = υ . GIT 2

Die resultierende Gesamtverschiebung (Abb. A3-4) setzt sich aus einem Biege-, Schub- und Torsionsanteil zusammen. F wB + wS

wT Abb. A3-4: Gesamtverschiebung

Somit folgt für die Gesamtverschiebung Fl 3 Fl b 2 Fl w= + + 3EI GAS 4 GIT

und mit I=

bt 3 1 5 , IT = bt 3 , AS = bt 12 3 6

wird die resultierende Verschiebung zu

w=


Fl  4l 2 6t 2 3b 2  + +  . bt 3  E 5G 4G 


71

A3 Anhang

A3.2 Test 2. Eigenfrequenz Biegebalken In Test 2 werden die ersten drei Biegeeigenfrequenzen eines kurzen dicken Balkens (Abb. A3-5)

Balkenlänge l 30 = ermittelt. Dabei ist der Balken statisch bestimmt Balkentiefe t 3

gelagert und der Einfluss der Rotationsträgheit sowie des Schubes werden berücksichtigt (Timoshenkotheorie).

w ( x, t )

z

x

Abb. A3-5: Gedrungener Balken

Schwingt der Balken um seine y-Achse, so erfolgt eine Verschiebung w ( x, t ) in zRichtung (Abb. A3-5) und eine Drehung ψ ( x, t ) um die y-Achse (Abb. A3-6).

ψ ( x, t ) Abb. A3-6: Verdrehung Balken

Folglich lassen sich für das Balkenelement der Schwerpunktsatz in z-Richtung und der Drallsatz durch den Schwerpunkt (parallel zur y-Achse) formulieren (Abb. A3-7). ρAdx ρI

∂2w ∂Q = −Q + Q + dx 2 ∂t ∂x

∂ 2ψ ∂M dx = − M + M + dx − Qdx 2 ∂t ∂x M+

M

z

Q dx

∂M dx ∂x

∂Q Q+ dx ∂x

x

Abb. A3-7: Balkenelement



72

A3 Anhang

Weiterhin wird das Elastizitätsgesetz für das Biegemoment M = EI ψ′

und für die Querkraft

Q = GAS ( w′ + ψ ) benötigt. Durch die vier Differentialgleichungen werden der Biegeeinfluss, der Einfluss der Schubdeformation sowie die Drehträgheit beschrieben. Um

die

Bewegungsgleichung

des

Balkens

zu

erhalten,

wird

das

Momentengleichgewicht in der Form ρI

∂2 ∂2 ′ w + ψ − ρ I w′ − M ′ + Q = 0 ( ) ∂t 2 ∂t 2

geschrieben. Durch Einsetzen des Elastizitätsgesetzes für das Biegemoment in der Form M = EI ( w′ + ψ )′ − EIw′′ und Differenzieren nach x wird die Gleichung für das Momentengleichgewicht zu ρI

∂2 ∂2 ′ ′ w + ψ − ρ I w′′ − EI ( w′ + ψ )′′′ + EIw IV + Q′ = 0 . ) 2 ( 2 ∂t ∂t

Mit Q′ = ρA

∂2 Q w und ( w′ + ψ ) = 2 ∂t GAS

wird die Differentialgleichung zu

 EA  IA && && − ρI 1 + &&′′ + ρ2 && = 0 . EIw IV + ρAw w w GAS  GAS  Dies ist die Bewegungsgleichung einer freien Schwingung nach der Timoshenko Balkentheorie [TIM55]. Mit dem Ansatz für harmonische Schwingungen

w ( x, t ) = W ( x ) cos ( ωt − α ) ergibt sich die Differentialgleichung

EI

d 4W EA  d 2W 2 2 2 IA − ρ A ω W − ρ I ω 1 + ω4W = 0 .   2 +ρ 4 dx GAS  GAS  dx



73

A3 Anhang

Für den Fall des beidseitig frei drehbaren Balkens und unter Einbeziehung der Randbedingungen gilt folgende Lösungsansatz Wk ( x ) = Bk sin

k πx , k = 1, 2.... . l

Durch Einsetzen der Lösung in die Differentialgleichung werden die Eigenwerte bestimmt und es ergibt sich  EA  kπ  EI   − ρAω2 − ρI ω2 1 +  l   GAS 4

Bernoulli

 kπ  2 AI ω4 = 0   + ρ l GA   S  2

Timoshenko

(Rotationsträgheit & Schubeinfluss)

Abb. A3-8: Bewegungsgleichung Balken

oder mit a 4 = ω2

ρA I und i 2 = folgt EI A 4

2

EA   k π  EA  kπ  4 4 2 − a − a i 1 + + a8i 4 =0    l   GAS    GAS   l  In Tabelle A3-1 sind weitere charakteristische Gleichungen für verschiedene Lagerungen aufgeführt. κk l

Lagerung

charakt. Gleichung

gelenkig-gelenkig

sin κl = 0

eingespannt-gelenkig

tan κl − tanh κl = 0

≈ ( 4k + 1)

π 4

cosh κl cos κl + 1 = 0

≈ ( 2k − 1)

π 2

≈ ( 2k + 1)

π 2

eingespannt-frei eingespannt-eingespannt

cosh κl cos κl − 1 = 0

frei-frei

cosh κl cos κl − 1 = 0

κπ

Tabelle A3-1: Charakteristische Gleichungen für verschiedene Randbedingungen [GRO99]

Wie in Abb. 2-1 gezeigt, lassen sich die einzelnen Terme der Bewegungsgleichung den verschiedenen Theorien zuordnen. Werden nur die ersten beiden Terme betrachtet, so erhält man den Euler-Bernoulli-Balken mit Fachhochschule Aachen


74

A3 Anhang

ωk = k 2 π2 Der

zweite

Term

beschreibt

die

EI . ρAl 4 Erweiterungen

um

Schub-

und

Rotationsträgheitseinfluss, welche erstmals von Stepan Prokofyevich Timoshenko [MIN51] eingeführt wurden. Dabei berücksichtigt der erste Teil die Rotationsträgheit und der Zweite den Schubeinfluss. Durch Hinzufügen des einen oder anderen Terms lässt sich so die Bewegungsgleichung herleiten, die nur die Rotationsträgheit oder den Schubeinfluss mit einbezieht. Der letzte Term, der die Rotationsträgheit und Schubeinfluss koppelt, ist bei praktischen Anwendungsfällen, bei denen davon ausgegangen wird, dass

kiπ << 1 , von sekundärer Bedeutung (es ergibt sich auch l

keinerlei technische Aussage, die aus diesem Term gewonnen werden könnte). Um die relative Größe des letzten Terms im Vergleich zum „Timoshenko-Term“ zu 4

 kπ  bilden, wird dieser mittels a =   zu  l  4

EA  kπ  ai ≡ a 4i 2   GAS  l  8 4

2

 kiπ 2 EA      l  GAS 

umgeformt. Nun wird erkennbar, dass 2 2 2 EA   kπ   kiπ  EA  4 2  kπ   a i      << a i   1 + ,  l   l  GÁS   l   GÁS  4 2

da

kiπ << 1 ist. Der letzte Teil wird für die Bestimmung der Eigenfrequenzen nicht l

weiter berücksichtigt. Die Eigenkreisfrequenzen nach der Timoshenko-Theorie ergeben sich somit zu

 EA  kπ   ωk =   1 + i 2 1 +  l    GAS 2


2  kπ        l  

− 12

EI . ρA


75

A3 Anhang

In Abb. A3-9 findet ein Vergleich der verschiedenen Balkentheorien statt. Es wird deutlich, dass ab einer bestimmten Wellenlänge (Wellenlänge und geschwindigkeit sind dimensionslos angegeben.) bzw. einer entsprechenden Frequenz, der Einfluss des Schubes sowie die Rotationsträgheit nicht mehr vernachlässigt werden dürfen. Es ist bemerkenswert wie gut die TimoshenkoTheorie mit dem exakten Ergebnis übereinstimmt. Im Vergleich dazu weist die Bernoulli- und Rayleightheorie (nur Rotationsträgheit) zu hohe Werte auf und die Struktur wird zu steif bewertet.

Abb. A3-9: Vergleich der Balkentheorien [GRA91]



76

A3 Anhang

Durch eine Literaturrecherche wurden weitere Formeln für die Timoshenko-Theorie recherchiert. Die Differenzen dieser Formeln lassen sich dadurch erklären, dass verschiedene Näherungen oder Potenzreihenentwicklungen benutzt wurden. Anzumerken

ist

dabei,

dass

alle

Autoren

den

letzten

Term

in

der

Differenzialgleichung wegfallen lassen und dass durch die unterschiedlichen Näherungen die Ergebnisse sich teilweise stark voneinander unterscheiden. Tabelle A3-2 listet die Ergebnisse der ersten drei Biegeeigenfrequenzen für die verschiedenen Balkentheorien auf. Deutlich erkennbar ist, wie sich der Schub- bzw. Rotationsträgheitseinfluss bei höheren Frequenzen auf das Ergebnis auswirkt.

Quelle [TIM55] Bernoulli [TIM55] 1. Bernoulli+Rotationsträgheit [NUK99] 2. Bernoulli+Schubeinfluss 3. [TIM55] Diplomarbeit 4. [GRO99] 5. [CLO75]

Erste Zweite Dritte Biegeeigen- Biegeeigen- Biegeeigenfrequenz [Hz] frequenz [Hz] frequenz [Hz] 7817,769 31271,075 70359,919 7785,619

30756,686

67755,826

7743,608

30132,725

64958,676

7685,314 7688,588 7619,086 7718,427

29151,794 29345,472 28092,154 29681,614

59631,059 61592,089 54266,629 62313,274

Tabelle A3-2: Vergleich der Eigenfrequenzen nach Quellen

2

1.

 kπ  ωk =    l 

2 EI  1  k πi   1 −    ρA  2  l  

2.

 kπ  ωk =    l 

2

2  EI    k πi   EA     1 +   ρA    l   GAS

2

3.

 kπ  ωk =    l 

EI  1  k πi  1 −   ρA  2  l 

4.

ωk = k π

5.

ωk = k 2 π2

2


2

2

   

 EA 1 +  GAS

−

1 2

   

2   EA   k πi   1 − 1 +      GAS   l   2 EI  1  kiπ   EA   1−    1 + 4  ρAl  2  l   GAS  

EI ρAl 4


77

A3 Anhang

Die verschiedenen Werte in Tabelle A3-2 ergeben sich hauptsächlich durch numerische Ungenauigkeiten. Bei einer Vergleichsrechnung der durch EXCEL und mittels Matlab ermittelten Werte, wurden weitere Abweichungen beobachtet. Auch ist zu beobachten, dass verschiedene Werte für die Formeln 3 und 5 errechnet werden, obwohl die Formeln sich nur durch die Position der Länge (In der Wurzel oder vor die Wurzel gezogen.) unterscheiden.



78

A3 Anhang

A3.3 Test 3. Eigenfrequenzen fest eingespannte Platte Beim dritten Testbeispiel werden die ersten sechs Eigenfrequenzen für eine fest eingespannte homogene quadratische Platte (Abb. A3-10) berechnet. Für eine Herleitung der Schwingungsgleichung für eine Platte, die die klassische Plattengleichung um Schub- und Rotationsträgheit erweitert, sei an dieser Stelle auf [MIN51] verwiesen. Alle h

Ränder

fest

eingespannt.

a=b

a b Abb. A3-10: Platte

Die Differentialgleichung der Platte lautet in Analogie zur Schwingungsgleichung nach Timoshenko für den Balken B∆ 2 w − ∆

ρh 3 ∂ 2 ρ ∂2 ρ2 h3 ∂ 4 ∂2 w − B ∆ w + w + ρ h w=0 12 ∂t 2 G′ ∂t 2 12G′ ∂t 4 ∂t 2

mit B=

Eh3 Biegesteifigkeit und G ′ = κ 2G 2 12 (1 − υ )

korrigierter Schubmodul mit κ 2 ≈ 0.3 × υ + 0.76 . Hierbei lassen sich die einzelnen Terme für Schub-, Biege- und Rotationsträgheitseinfluss identifizieren (Abb. A3-11):

ρh 3 ∂ 2 ρ ∂2 ρ2h3 ∂ 4 ∂2 B∆ w − ∆ w − B∆ w+ w + ρh 2 w = 0 12 ∂t 2 G ′ ∂t 2 12 G ′ ∂ t 4 ∂t 2

Biegeeinfluss

Rotationsträgheit

Schubeinfluss

(Rotationsträgheit & Schubeinfluss)

Abb. A3-11: Bewegungsgleichung Platte

mit ∆ =

∂2 : Laplace − Operator. ∂x 2 Fachhochschule Aachen


79

A3 Anhang

Der Fall der quadratischen Platte ist ein Spezialfall der rechteckigen Platte. Die Lösung

der

Bewegungsgleichung

mit

der

Ansatzfunktion,

die

auch

die

Randbedingungen mitberücksichtigt, würde den Umfang dieser Arbeit weit überschreiten. Meist wird auch in der Fachliteratur nur der Fall der rechteckigen Platte bzw. Membran mit Auflagerrandbedingungen (an den Rändern frei drehbar) behandelt. Einen umfangreichen Überblick zum Thema Schwingungen von Platten gibt Leissa [LEI61] in seiner Arbeit „Vibration of Plates (NASA SP-160)“, die im Auftrag der NASA erstellt wurde. Dabei werden verschiedene Randbedingungen sowie Formen von Platten ausführlich behandelt. Für dem hier vorliegenden Fall ergeben sich aus der Veröffentlichung für die ersten sechs Eigenfrequenzen nach Young [LEI61], mit D =

Eh 3 und ρ = ρ × h , 12(1 − υ 2 )

folgende Werte:

Eigenmode 1 Young

2

3

4

5

6

1793,063 3657,370 5394,134 6558,454 6588,845 8227,961 Tabelle A3-3: Eigenfrequenzen

Die Eigenformen einer fest eingespannten quadratischen Platte können Abb. A3-12 entnommen werden.

Abb. A3-12: Eigenformen einer quadratischen Platte mit Seitenlänge a und Biegesteifigkeit D[LEI61]

Leider wird weder der Schubeinfluss noch der Einfluss der Rotationsträgheit in den oben genannten Formeln berücksichtigt. Die in Abb. A3-11 betrachtete Formel beinhaltet alle Einflüsse, die auf eine Platte wirken und ist aus der Arbeit von Mindlin „Influence of Rotatory Inertia and Shear on Flexural Motions of Isotropic, Elastic Plates” [MIN51] hergeleitet. Welche Auswirkungen diese Terme auf die Eigenfrequenzen haben, geht aus Abb. A3-13 Fachhochschule Aachen


80

A3 Anhang

hervor. Es lässt sich erkennen, dass Oberhalb eines Verhältnisses von Plattendicke zu

Wellenlänge,

l i

λ= =

Wellenlänge , i − te Schwingung

von

h ≈ 0,125 λ

die

klassische

Plattentheorie nicht mehr gilt. In unserem Beispiel würde bei einer Plattendicke von 2mm eine charakteristische Wellenlänge von λ ≈ 16mm vorliegen (Abb. A3-14). Damit ergibt sich eine Grenzfrequenz, für die die Rotationsträgheit und der Schubeinfluss nicht mehr vernachlässig werden darf, von f ≈ 760 Hz .

Abb. A3-13: Schub- und Rotationsträgheitseinfluss auf die Wellengeschwindigkeit [MIN61]



81

A3 Anhang

Obwohl nach der allgemeinen Definition hier eine dünne Platte mit dem Verhältnis h h 2mm = = a b 100mm

vorliegt (siehe Tabelle A3-4), zeigen doch die vorherigen

Annahmen, dass diese Definition nicht für die Bestimmung der Eigenfrequenzen zutrifft (Dynamik). Die wesentlichen Vereinfachungen der Plattentheorie sind: •

linear elastisches Materialverhalten (gilt nur für kleine Verformungen)

•

konstante Plattendicke h

•

die Spannungen normal zur Schalenmittelfläche sind vernachlässigbar ( σ Z = 0, ε Z = 0 ).

•

dicke Platte:

Reissner-Mindlin (Timoshenko-Theorie bei Balken): Die Normalen zur

Mittelfläche sind vor der Verformung gerade, nach der Verformung auch gerade, müssen jedoch nicht senkrecht zur Mittelfläche stehen. •

dünne Platte:

Kirchhoff-Theorie (Bernoulli-Theorie bei Balken): Die Normalen zur

Mittelfläche stehen auch nach der Verformung noch senkrecht auf der Mittelfläche

h h ; a b

Plattentheorie

Dicke Platte

Dünne Platte

1 1 L 5 10

1 1 L 10 50

Reissner-Mindlin

Kirchhoff

(Schubweich)

(Schubstarr)

Tabelle A3-4: Einsatzgebiete Plattentheorie

Besonders bei der Betrachtung von höheren Frequenzen, wo in Analogie zur Balkentheorie der Schubeinfluss mit höheren Frequenzen immer weiter dominiert, darf die Theorie nach Kirchhoff nicht benutzt werden. Dieser Theorie folgend, würde sich die Struktur zu steif verhalten, wodurch es zu höheren Werten der Eigenfrequenzen kommt.



82

A3 Anhang

Abb. A3-14: Charakteristische Wellenlänge ( λ ) über Frequenz ( Eisenplatten [CRE67]



f ) für

83

A4 Anhang

A4 Literaturverzeichnis [TIM55]

Timoshenko, Stephan P.: Vibration Problems in Engineering, New York 1955

[CLO75]

Clough, Ray W.; Penzien, Joseph : Dynamics of Structures, Tokio 1975

[GRO99]

Gross, Dietmar; Hauger, Werner; Schnell, Walter; Wriggers, Peter: Technische Mechanik Bnd. 4, Berlin 1999

[NUK99]

Nukala, Phani K.: Implementation of Rotary Inertia Effect on the Free Vibration Response of Beams, Lawrence Livermore National Laboratory 1999

[LEI69]

Leissa, Arthur W.: Vibration of Plates (NASA SP-160), Washington D.C. 1969

[MIN51]

Mindlin, Raymond D., Influence of Rotatory Inertia and Shear on Flexural Motions of Isotropic, Elastic Plates, New York 1951

[CRE67]

Cremer, Lothar; Heckl, Manfred : Körperschall, Berlin 1967

[GRA91]

Graff, Karl F.: Wave Motion in Elastic Solids, New York 1991

[Cal1.6]

Dokumentation CalculiX, http://www.dhondt.de/

[B2000]

Dokumentation B2000, http://www.smr.ch/local/doc/B2000/html/index.html

[OOFem]

Dokumentation Object Oriented Finite Element Solver, http://www.oofem.org/documentation/manual.html

[DAK07]

Homepage Dakota (Stand Julie 2007), http://www.cs.sandia.gov/DAKOTA/software.html,

Sandia

National

Labroatories [BAI94]

Baier, Horst, Seeßelberg, Christoph, Specht, Bernhard : Optimierung in der Strukturmechanik, Wiesbaden 1994



84

Diplomarbeit CalculiX Hiller

Recommend Documents