Finite-Elemente-Methode

Visualisierung einer FEM-Simulation der Verformung eines Autos bei asymmetrischem Frontalaufprall

Darstellung der Wärmeverteilung in einem Pumpengehäuse mit Hilfe der Wärmeleitungsgleichung. Die „finiten Elemente“ sind mit den Elementkanten als schwarze Linien zu sehen.

Die Finite-Elemente-Methode (FEM), auch „Methode der finiten Elemente“ genannt, ist ein allgemeines, bei unterschiedlichen physikalischen Aufgabenstellungen angewendetes numerisches Verfahren. Am bekanntesten ist die Anwendung der FEM bei der Festigkeits- und Verformungsuntersuchung von Festkörpern mit geometrisch komplexer Form, weil sich hier der Gebrauch der klassischen Methoden (z.B. die Balkentheorie) als zu aufwändig oder nicht möglich erweist. Logisch basiert die FEM auf dem numerischen Lösen eines komplexen Systems aus Differentialgleichungen.

Das Berechnungsgebiet (z.B. der Festkörper) wird in endlich viele Teilgebiete (z.B. Teilkörper) einfacher Form aufgeteilt, z.B. in viele kleine Quader oder Tetraeder. Sie sind die „finiten Elemente“. Ihr physikalisches Verhalten kann aufgrund ihrer einfachen Geometrie mit bekannten Ansatzfunktionen gut berechnet werden. Das physikalische Verhalten des Gesamtkörpers wird dadurch nachgebildet, wie diese Elemente auf die Kräfte, Lasten und Randbedingungen reagieren und wie sich Lasten und Reaktionen beim Übergang von einem Element ins benachbarte fortpflanzen durch ganz bestimmte problemabhängige Stetigkeitsbedingungen, die die Ansatzfunktionen erfüllen müssen.

Die Ansatzfunktionen enthalten Parameter, die in der Regel eine physikalische Bedeutung besitzen, wie z.B. die Verschiebung eines bestimmten Punkts im Bauteil zu einem bestimmten Zeitpunkt. Die Suche nach der Bewegungsfunktion ist auf diese Weise auf die Suche nach den Werten der Parameter der Funktionen zurückgeführt. Indem immer mehr Parameter (z.B. immer mehr, kleinere Elemente) oder immer höherwertige Ansatzfunktionen benutzt werden, kann die Genauigkeit der Näherungslösung verbessert werden.

Die Entwicklung der FEM war in wesentlichen Etappen nur mittels der Entwicklung leistungsfähiger Computer möglich, da sie erhebliche Rechenleistung benötigt. Daher wurde diese Methode von vornherein computergerecht formuliert. Sie brachte einen wesentlichen Fortschritt bei der Behandlung von Berechnungsgebieten beliebiger Form.

Einführung

Generiertes Strukturmodell eines Hubkolbens (Dieselmotor) als Komponente eines Verbrennungsmotors zum Zwecke der Spannungsanalyse

Mit der FE-Methode können Probleme aus verschiedenen physikalischen Disziplinen berechnet werden, da es sich grundsätzlich um ein numerisches Verfahren zur Lösung von Differentialgleichungen handelt. Zunächst wird das Berechnungsgebiet („Bauteil“) in eine große Anzahl von Elementen unterteilt – ausreichend fein. Diese Elemente sind endlich klein (finit), ihre tatsächliche Größe bleibt jedoch mathematisch relevant – sie sind nicht „unendlich klein“ (infinit). Das Aufteilen des Gebiets/Bauteils in eine bestimmte Anzahl Elemente finiter Größe, die sich mit einer endlichen Zahl von Parametern beschreiben lassen, gab der Methode den Namen „Finite-Elemente-Methode“.

Für diese Elemente gibt es Ansatzfunktionen (z.B. lokale Ritz-Ansätze je Element), die beschreiben, wie ein Element auf äußere Einflüsse und Randbedingungen reagiert. Setzt man diese Ansatzfunktionen für alle Elemente in die zu lösenden Differentialgleichungen ein, die die physikalischen Gesetze beschreiben, erhält man zusammen mit den Anfangs-, Rand- und Übergangsbedingungen ein meist sehr großes Gleichungssystem. Es (zumindest näherungsweise) zu lösen ist die Aufgabe des FE-Gleichungslösers. Die Größe des zu lösenden Gleichungssystems hängt maßgeblich von der Anzahl der finiten Elemente ab. Seine Näherungslösung stellt letztlich die numerische Lösung der betrachteten Differentialgleichung dar; wenn für alle Elemente gelöst ist, wie sie sich unter den Lasten verhalten, so hat sich dadurch auch die Reaktion des Gesamtbauteils ergeben.

Finit, Infinit

Mathematisch bleibt die Größe jedes Elements relevant und muss auch in seine Berechnung einfließen, es ist nur 'finit' klein. Bei 'infinit' kleinen Elementen wäre ihre Größe vernachlässigbar und würde in den Gleichungen nicht mehr berücksichtigt. Diesbezüglich bleibt die Elementgröße also relevant.

Eine „ausreichend feine“ Aufteilung des Bauteils in Elemente liegt vor, wenn eine weitere Verfeinerung keinen signifikanten Einfluss auf das Rechenergebnis mehr hat. D. h. das Gesamtergebnis wird diesbezüglich unabhängig von der Elementgröße, die (aus dieser Sichtweise) dann nicht mehr relevant ist. Hat die Elementgröße noch nennenswerten Einfluss auf das Gesamtergebnis, dann gilt i. A. die Vernetzung als nicht fein genug.

Geschichte

Der Einsatz der FEM in der Praxis begann in den 1950er Jahren bei einer Strukturberechnung von Flugzeugflügeln in der Luft- und Raumfahrtindustrie (Turner, Clough 1956) und sehr bald auch im Fahrzeugbau. Die Methode basiert hier auf den Arbeiten bei der Daimler AG in Stuttgart, die das selbst entwickelte FEM-Programm ESEM (Elastostatik-Element-Methode) einsetzte, lange bevor die computerunterstützte Konstruktion (CAD) Anfang der 1980er Jahre ihren Einzug hielt. Der Ausdruck Finite-Elemente-Methode wurde erstmals 1960 von R. W. Clough vorgeschlagen und wird seit den 1970er Jahren überall verwendet. Die gängigste deutschsprachige Bezeichnung für industrielle Anwender ist Berechnungsingenieur.

Anwendung

Die erste Anwendung der FEM war die lineare Behandlung von Festkörpern und Strukturen in Form der Verschiebungsmethode und davon ausgehend hat die FEM ihre Anstöße erhalten. Die Bezeichnung „Finite Elemente“ wurde erst etwas später benutzt. Im weiteren Verlauf der Forschung wurde die Finite-Elemente-Methode immer weiter verallgemeinert und kann nunmehr in vielen physikalischen Problemstellungen, u.a. in verschiedenen gekoppelten Feldberechnungen, Wettervorhersagen oder bei technischen Fragestellungen in den Branchen Fahrzeugbau, Medizintechnik, Luft- und Raumfahrttechnik, Maschinenbau oder Konsumgüter in den Ingenieurwissenschaften verwendet werden. Ein Haupteinsatzgebiet der Methode ist die Produktentwicklung, wobei unter anderem mechanische Festigkeitsberechnungen einzelner Komponenten oder beispielsweise komplette Fahrwerks- und Karosseriestrukturen berechnet werden, um aufwändige Crashtests zu sparen.

Vorgehen einer linearen mechanischen Berechnung (exemplarisch)

Programme, welche die Finite-Elemente-Methode verwenden, arbeiten nach dem EVA-Prinzip: Der Anwender erstellt in einem CAD-Programm eine (Bauteil-)Geometrie. Anschließend gibt er im sogenannten FE-Präprozessor weitere Eingaben vor. Ein FEM-Gleichungslöser führt die eigentliche Rechnung durch, und der Benutzer erhält die berechneten Ergebnisse, welche er dann im sogenannten FE-Postprozessor in Form grafischer Anzeigen betrachten kann. Oft sind Prä- und Postprozessor in einem Programm kombiniert oder sogar Bestandteil des CAD-Programms.

Prozesskette lineare Festigkeitsrechnung

Eingabe: Präprozessor

Im CAD-Programm wird das Bauteil konstruiert und mittels einer Direktschnittstelle oder mit einem neutralen Austauschformat wie STEP in den FE-Präprozessor übertragen. Durch die Anwahl von Netzparametern wie Elementgröße und Elementart (z.B. Tetraeder, Hexaeder bei 3D) im Vernetzungsmodul werden mit Hilfe eines Vernetzungsalgorithmus die Finiten-Elemente erzeugt. Für die mechanische Festigkeitsanalyse ist das Materialverhalten einzugeben, das maßgeblich angibt, welche Reaktionen das Bauteil auf äußere Belastungen einnimmt (z.B. Verformung). Je nach Werkstoff ist der Zusammenhang zwischen Spannung und Dehnung unterschiedlich und es ergeben sich verschiedene Verformungen. Wenn dieser Zusammenhang linear ist, werden für die FE-Berechnung lediglich der Elastizitätsmodul und die Poissonzahl benötigt, sonst sind weitere Werkstoffkennwerte und Eingaben im Präprozessor nötig. Weitere Randbedingungen sind zum Beispiel einwirkende Belastungen auf das Bauteil (Kräfte, Druck, Temperatur etc.). Um eine möglichst realitätsnahe Abbildung zu erhalten, werden schließlich die homogenen (Einspannungen) und die inhomogenen Randbedingungen (Verschiebungen) sowie alle zu berücksichtigenden Lasten auf das Modell angegeben.

Verarbeitung: Gleichungslöser

Je nach Programm kommt nun ein separater (eigenständiges Programm) oder ein integrierter Gleichungslöser zum Einsatz. Er berechnet die Simulation, wie sich die Lasten, Kräfte und Randbedingungen auf die Einzelelemente des Bauteils auswirken, und wie sich die Kräfte sowie die Auswirkungen im Bauteil fortpflanzen und auf benachbarte Elemente auswirken. Die Berechnung liefert zunächst eine grobe Näherungslösung. Weitere Iterationen verbessern die Näherung. Meist werden so viele Iterationen berechnet, bis sich nur noch geringste Änderungen ergeben – dann hat die Näherung „konvergiert“ und ist das Ergebnis der Simulation.

Ausgabe: Postprozessor

Im Falle der mechanischen Festigkeitsberechnung erhält der Benutzer als Ergebnis des FEM-Gleichungslösers insbesondere Spannungs-, Deformations- und Dehnungswerte. Diese kann der Postprozessor zum Beispiel in einem Falschfarbenbild darstellen. Die Vergleichsspannungswerte werden beispielsweise zum Festigkeitsnachweis eines Bauteils verwendet.

Allgemeine Funktionsweise

Diskretisierung

Die Finite-Elemente-Methode ist ein diskretes Verfahren, d.h. die Lösung wird auf einer diskreten Untermenge des Grundgebietes berechnet. Hierzu wird dieses in einfache Teilgebiete, die so genannten finiten Elemente zerlegt (Vernetzung, Vermeshen). Die Bezeichnung „finit“ hebt den Unterschied zur analytischen Betrachtung auf infinitesimalen Elementen hervor. Die Ecken der finiten Elemente heißen Knoten. Diese Knoten bilden die diskrete Untermenge für das numerische Verfahren. Auf den Elementen werden Approximationsfunktionen eingeführt, welche die unbekannten Knotengrößen als Parameter enthalten. Die lokalen Approximationen werden in die schwache Formulierung des Randwertproblems eingeführt. Die dabei entstehenden Elementintegrale werden mit numerischer Quadratur berechnet. Dabei werden die Approximationsansätze „herausintegriert“, so dass auf den Elementen nach der Integration nur noch die Knotenwerte als Unbekannte verbleiben. Durch Kontinuitätsforderungen an den Elementgrenzen werden die Elementgleichungen assembliert. Auf diese Weise werden Randwertprobleme für lineare partielle Differentialgleichungen in ein lineares Gleichungssystem mit symmetrischen Systemmatrizen überführt. Für nichtlineare Differentialgleichungen verläuft der Algorithmus analog mit dem Unterschied, dass die nichtlinearen Abhängigkeiten mit geeigneten Methoden (z.B. Newton-Verfahren) iterativ linearisiert werden und das lineare Gleichungssystem in jedem Teilschritt für inkrementelle Größen aufgestellt wird.

Beispiel für die Anwendung eines adaptiven Gitters zur Berechnung der Luftströmung um einen Flugzeugflügel.

Bei gewissen Aufgabenstellungen ist die Unterteilung in Elemente durch das Problem bereits weitgehend vorgegeben, zum Beispiel bei räumlichen Fachwerken, bei denen die einzelnen Stäbe die Elemente der Konstruktion bilden. Das gilt auch bei Rahmenkonstruktionen, wo die einzelnen Balken oder unterteilte Balkenstücke die Elemente der Aufgabe darstellen. Bei zweidimensionalen Problemen wird das Grundgebiet in Dreiecke oder Vierecke eingeteilt. Selbst wenn nur geradlinige Elemente verwendet werden, erreicht man mit einer entsprechend feinen Diskretisierung eine recht gute Approximation (Annäherung) des Grundgebietes. Krummlinige Elemente erhöhen die Güte der Annäherung. Jedenfalls erlaubt diese Diskretisierung eine flexible und auch dem Problem angepasste Erfassung des Grundgebietes. Allerdings muss darauf geachtet werden, dass sehr spitze oder überstumpfe Winkel an den Element-Eckknoten vermieden werden, um numerische Schwierigkeiten auszuschließen. Dann wird das gegebene Gebiet durch die Fläche der approximierenden Elemente ersetzt. Mit dem Patch-Test kann man später überprüfen, ob das gut gelungen ist.

Räumliche Probleme werden mit einer Unterteilung des dreidimensionalen Gebietes in Tetraederelemente, Quaderelemente oder andere dem Problem angepasste, möglicherweise auch krummflächig berandete Elemente, dies sind i.d.R. Serendipity- oder Lagrange-Elemente, bearbeitet.

Die Feinheit der Unterteilung, d.h. die Dichte des Netzes, hat maßgeblichen Einfluss auf die Genauigkeit der Resultate der Näherungsrechnung. Da gleichzeitig der Rechenaufwand bei der Verwendung feinerer und dichterer Netze steigt, gilt es, möglichst intelligente Vernetzungslösungen zu entwickeln.

Element-Ansatz

In jedem der Elemente wird für die gesuchte Funktion, bzw. allgemeiner für die das Problem beschreibenden Funktionen, ein problemgerechter Ansatz gewählt. Im Besonderen eignen sich dazu ganze rationale Funktionen in den unabhängigen Raumkoordinaten. Für eindimensionale Elemente (Stäbe, Balken) kommen Polynome ersten, zweiten, dritten und gelegentlich sogar höheren Grades in Frage. Bei zweidimensionalen Problemen finden lineare, quadratische oder höhergradige Polynome Verwendung. Die Art des Ansatzes hängt dabei einerseits von der Form des Elementes ab, und andererseits kann auch das zu behandelnde Problem den zu wählenden Ansatz beeinflussen. Denn die Ansatzfunktionen müssen beim Übergang von einem Element ins benachbarte ganz bestimmte problemabhängige Stetigkeitsbedingungen erfüllen. Die Stetigkeitsanforderungen sind häufig aus physikalischen Gründen offensichtlich und aus mathematischen Gründen auch erforderlich. Zum Beispiel muss die Verschiebung eines zusammenhängenden Körpers in einer Richtung beim Übergang von einem Element zum anderen stetig sein, um die Kontinuität des Materials zu gewährleisten. Im Fall der Balken- oder Plattenbiegung sind die Stetigkeitsanforderungen höher, da dort aus analogen physikalischen Gründen sogar die Stetigkeit der ersten Ableitung bzw. der beiden ersten partiellen Ableitungen gefordert werden muss. Elemente mit Ansatzfunktionen, welche den Stetigkeitsbedingungen genügen, heißen konform.

Um nun die Stetigkeitsanforderungen tatsächlich zu erfüllen, muss der Funktionsverlauf im Element durch Funktionswerte und auch durch Werte von (partiellen) Ableitungen (den Knotenpunktverschiebungen) in bestimmten Punkten des Elementes, den Knotenpunkten, ausgedrückt werden. Die in den Knotenpunkten benutzten Funktionswerte und Werte von Ableitungen nennt man die Knotenvariablen des Elements. Mit Hilfe dieser Knotenvariablen stellt sich die Ansatzfunktion als Linearkombination von sogenannten Formfunktionen mit den Knotenvariablen als Koeffizienten dar.

Es ist zweckmäßig, für die Knotenpunktkoordinaten neben einem elementbezogenen lokalen ein globales Koordinatensystem zu verwenden. Beide werden durch Transformationsfunktionen miteinander verknüpft. Werden für diese Transformation dieselben Formfunktionen wie für den Verformungsansatz benutzt, so sind es isoparametrische Elemente, bei Funktionen niedrigeren bzw. höheren Grades sub- bzw. superparametrische Elemente.

Randbedingungen

Problemstellung	Dirichlet-Randbedingung/Funktionswert	Neumann-Randbedingung
statisches Problem	Auflagerbedingung/Verschiebung	Kraft
Sickerströmung	Standrohrspiegelhöhe>	Quelle oder Senke
Wärmeleitung	Temperatur	Wärmestrom bzw. Wärmestromdichte
elektrischer Strom	elektrische Spannung	Stromstärke
Elektrostatik	elektrische Spannung	elektrische Ladung
Magnetostatik	magnetisches Potenzial	magnetischer Fluss

Nachdem ein gegebenes Problem diskretisiert ist und die Elementmatrizen aufgestellt sind, führt man vorgegebene Randbedingungen ein. Ein typisches FE-Problem kann zwei Arten von Randbedingungen haben: Dirichlet-Randbedingungen und Neumann-Randbedingungen. Sie gelten (wirken) immer an den Knotenpunkten.

Eine Dirichlet-Randbedingung gibt einen Funktionswert direkt vor und eine Neumann-Randbedingung gibt eine Ableitung eines Funktionswertes vor. Ist eine Dirichlet-Randbedingung vorgegeben, bedeutet dies, dass das Problem einen Freiheitsgrad weniger bekommt und die zugehörige Zeile und Spalte in der Gesamtmatrix gestrichen wird. Ist die Dirichlet-Randbedingung ungleich Null, so wird der Wert entsprechend seinem Vorfaktor der Linearform („rechten Seite“) hinzugefügt. Je nach Art des physikalischen Problems kann es sich um verschiedene physikalische Größen handeln, wie in der Tabelle beispielhaft dargestellt. Die Neumann-Randbedingungen haben des Weiteren einen Anteil an der Linearform („rechte Seite“).

Eine weitere Variante sind periodische Randbedingungen, bei denen die Werte an einem Rand als Daten für einen anderen Rand genommen werden und so ein periodisch unendlich fortgesetztes Gebiet simuliert wird. Für rotationssymmetrische Probleme werden sogenannte zyklische Randbedingungen definiert.

Grundgleichungen der Verschiebungsmethode

Die Verschiebungsmethode ist die Standardformulierung der Finite-Elemente-Methode, bei der die Verschiebungen die primären Unbekannten sind, die die Translation, Rotation und Verformung eines Festkörpers beschreiben. Die Verschiebungsmethode ist in allen gängigen Finite-Elemente-Programmen verfügbar, mit denen Probleme der Festkörpermechanik berechnet werden können. Für die Lösung von Festkörper-Problemen liegen mehrere Grundgleichungen vor.

Prinzip von d'Alembert in der Lagrangeschen Fassung

Eine der Verschiebungsmethode zugrunde liegende Gleichung, mit der allgemeine Probleme der Festkörpermechanik behandelt werden können, ist das Prinzip von d’Alembert, wie es die Kontinuumsmechanik in der Lagrangeschen Beschreibung formuliert. Mit diesem Prinzip können sowohl lineare Probleme, wie die Frage nach Eigenschwingungen, als auch hoch nichtlineare Probleme, wie Crashtests, analysiert werden. Hier wird die Methode der gewichteten Residuen nach Galerkin, auch Galerkin-Methode oder Galerkin-Ansatz genannt, verwendet.

Prinzip vom Minimum der potenziellen Energie

In konservativen Systemen können bei einem statischen Problem die Knotenpunktverschiebungen aus der Bedingung ermittelt werden, dass im gesuchten Gleichgewichtszustand die potenzielle Energie ein Minimum hat. Mit dem Prinzip vom Minimum der potenziellen Energie können die Steifigkeitsgleichungen finiter Elemente direkt bestimmt werden. Die potenzielle Energie einer Konstruktion ist die Summe aus der inneren Verzerrungsenergie (der elastischen Formänderungsenergie) und dem Potenzial der aufgebrachten Lasten (der von äußeren Kräften geleisteten Arbeit).

Bogenlängenverfahren (arc-length method)

Das Bogenlängenverfahren ist eine Methode, bei der man kraftgesteuert bis über das Maximum der Traglast hinaus rechnen kann. Die Notwendigkeit von kraftgesteuerten Methoden liegt darin, dass man im Gegensatz zu verschiebungsgesteuerten Methoden mehrere Lasten direkt proportional steigern kann. Beim Bogenlängenverfahren wird die Last wie vorgegeben gesteigert; würde diese Belastungssteigerung zu einer zu großen Deformation führen, so wird die Last mit einem Faktor kleiner als 1 multipliziert, nach Erreichen der Traglast sogar mit negativen Werten.

Stochastische Finite-Elemente-Methode

Bei der Variante der stochastischen Finiten-Elemente-Methode (SFEM) werden Eingangsgrößen des Modells, welche mit einer Unsicherheit behaftet sind, zum Beispiel Materialfestigkeiten oder Belastungen, durch stochastische Größen modelliert. Dies kann mithilfe gewöhnlicher Zufallsvariablen erreicht werden. Oft werden auch Zufallsfelder verwendet, wobei es sich um zufällig variierende, stetige mathematische Funktionen handelt. Eine geläufige Berechnungsmethode ist dabei die Monte-Carlo-Simulation. Dabei wird die FE-Berechnung für viele zufällige Realisierungen (samples) der Eingangsgrößen wiederholt, bis man einen gewissen, im Vorfeld definierten, stochastischen Fehler unterschreitet. Anschließend werden aus allen Ergebnissen die Momente, also Mittelwert und Varianz, berechnet. Je nach Streuung der Eingangsgrößen sind oftmals sehr viele Wiederholungen der FE-Berechnung nötig, was viel Rechenzeit in Anspruch nehmen kann.

Implizite und explizite FE-Löser

Strukturmechanische FEM-Systeme werden durch lineare Gleichungssysteme 2. Ordnung dargestellt:

$M{\ddot {u}}(t)+D{\dot {u}}(t)+Ku(t)=P(t)$

$M,D$ und sind Massen-, Dämpfungs- und Steifigkeitsmatrix des Systems; P(t) ist der Vektor der externen Kräfte, die auf das Modell wirken. ist der Vektor der Freiheitsgrade.

Oft bestehen komplexe Bauteilmodelle aus mehreren Millionen Knoten, und jeder Knoten kann bis zu 6 Freiheitsgrade besitzen. Somit müssen FEM-Solver (Gleichungssystemlöser) gewisse Anforderungen in Bezug auf effektives Speichermanagement und ggf. Nutzung mehrerer CPUs erfüllen. Es gibt zwei grundsätzlich verschiedene Arten von FEM-Solvern: implizite und explizite.

Implizite FEM-Solver gehen von bestimmten Annahmen aus, unter denen der berechnete Lösungsvektor gültig ist. Wirkt z.B. eine zeitlich unveränderliche Last $P(t)=P=const.$ auf ein System mit Dämpfung, dann wird sich nach ausreichend langer Zeit auch ein konstanter Verschiebungsvektor $u(t)=u=const.$ einstellen. Für $t\to \infty$ ist dann ${\ddot {u}}(t)={\dot {u}}(t)=0$ , und das Gleichungssystem vereinfacht sich zu $Ku=P$ mit der Lösung $u=K^{-1}P$

Für einen gegebenen Lastvektor kann der Verschiebungsvektor mit Hilfe des Gauß-Algorithmus oder durch QR-Zerlegung von berechnet werden.

Ist ein mechanisches System einer harmonischen Anregung $P(t)={\hat {P}}\sin(\Omega t)$ ausgesetzt, dann kann es erforderlich sein, die Eigenfrequenzen des Systems zu ermitteln, um Resonanzen im Betrieb zu vermeiden.

Eigenfrequenzen sind alle Frequenzen $\omega _{i}$ , für die ein Verschiebungsvektor $u(t)={\hat {u}}_{i}\sin(\omega _{i}t)$ eine Lösung des unbelasteten ( $P(t)=0$ ) und ungedämpften ( D=0 ) Gleichungssystems darstellt. Für den Geschwindigkeits- und Beschleunigungsvektor gilt dann

${\dot {u}}_{i}={\hat {u}}_{i}\ \omega _{i}\cos(\omega _{i}t),\ {\ddot {u}}_{i}=-{\hat {u}}_{i}\omega _{i}^{2}\sin(\omega _{i}t)$

und das Gleichungssystem lautet damit

$(-M\omega _{i}^{2}+K){\hat {u}}_{i}\sin(\omega _{i}t)=0$

Um die Eigenfrequenzen $\omega _{i}$ und die dazugehörigen Eigenformen ${\hat {u}}_{i}$ zu berechnen, muss der implizite Solver also das Eigenwertproblem

$-M\omega _{i}^{2}+K=0$

lösen.

Explizite FEM-Solver

Explizite Solver berechnen die Verschiebungsvektoren $u_{i}$ zu bestimmten diskreten Zeitpunkten $t_{i}$ innerhalb eines vorgegebenen Zeitintervalls. Knotengeschwindigkeiten und -beschleunigungen werden durch Differenzenquotienten aus den Verschiebungen $u_{i+1},u_{i},u_{i-1}$ zu aufeinanderfolgenden Zeitpunkten $t_{i+1},t_{i},t_{i-1}$ angenähert. Mit konstanter Zeitschrittweite $t_{i+1}-t_{i}=\Delta \ t$ gilt

${\dot {u}}={\frac {u_{i}-u_{i-1}}{\Delta \ t}}$

${\ddot {u}}={\frac {{\frac {u_{i+1}-u_{i}}{\Delta \ t}}-{\frac {u_{i}-u_{i-1}}{\Delta \ t}}}{\Delta \ t}}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{(\Delta t)^{2}}}$

hat das diskretisierte Gleichungssystem die Form

$M{\frac {u_{i+1}-2u_{i}+u_{i-1}}{(\Delta t)^{2}}}+D{\frac {u_{i}-u_{i-1}}{\Delta \ t}}+Ku_{i}=P_{i}$

Durch Auflösen dieser Gleichung erhält man eine Beziehung, mit der der Verschiebungsvektor $u_{{i+1}}$ aus den vorher berechneten Vektoren $u_{i}$ und $u_{i-1}$ ermittelt werden kann:

$u_{i+1}=M^{-1}((\Delta t)^{2}(P_{i}-Ku_{i})-D\Delta \ t(u_{i}-u_{i-1}))+2u_{i}-u_{i-1}$

Die Berechnung der Inversen $M^{{-1}}$ wird in der Praxis nicht durchgeführt, da explizite Solver in der Regel als Diagonalmatrix annehmen und daher jede Zeile des Gleichungssystems nur durch den Diagonaleintrag in der entsprechenden Zeile von geteilt werden muss.

Explizite Solver werden u.A. im Fahrzeugbau für die Berechnung von Crash-Lastfällen verwendet.

Der Vorteil direkter Gleichungslöser nach dem Gauß-Verfahren liegt für die praktische Anwendung in der numerischen Stabilität und dem Erhalt eines exakten Ergebnisses. Nachteilig sind die schlechte Konditionierung der üblicherweise dünn besetzten Steifigkeitsmatrizen und der hohe Speicherbedarf, wie oben erwähnt. Iterative Gleichungslöser sind unempfindlicher bei schlechter Konditionierung und benötigen weniger Speicher, wenn die Nicht-Null-Elemente-Speicherung verwendet wird. Allerdings verwenden iterative Solver ein Abbruchkriterium für die Berechnung der Ergebnisse. Wenn dieses erreicht wird, bevor eine annähernd exakte Lösung gefunden wurde, kann das Ergebnis, beispielsweise ein Spannungsverlauf, leicht fehlinterpretiert werden.

In manchen Implementierungen werden für die häufig auftretenden dünn besetzten Matrizen lediglich die Positionen und Werte der Einträge, die von Null abweichen, gespeichert. Damit kann man die Gleichungssysteme weiterhin direkt lösen, spart aber erheblich Speicherplatz.

Isogeometrische Analyse

Damit bei einer (kleineren) Änderung der (CAD-)Bauteil-Geometrie nicht aufwendig neu in finite Elemente unterteilt werden muss, können manche Programme ein bereits vorhandenes FE-Netz einer (sehr ähnlichen, neuen) CAD-Geometrie anpassen, was meist deutlich weniger Rechenzeit benötigt.

Programme

Finite-Elemente-Software und ihre Anwendung ist mittlerweile eine Industrie mit mehreren Milliarden US-Dollar Jahresumsatz.

In der Praxis sind viele verschiedene Standalone-Großprogramme mit ähnlichem Anwendungsspektrum im Einsatz; die Auswahl, welches Programm zum Einsatz kommt, ist nicht nur von der Verwendung, sondern auch von Faktoren wie Verfügbarkeit, Zertifizierungsstandard im Unternehmen oder Lizenzkosten abhängig.
Mit den in kommerziellen CAD-Systemen integrierten Finite-Elemente-Paketen können einfachere (i.d.R. lineare) Problemstellungen berechnet und mithilfe des CAD-Systems anschließend direkt ausgewertet werden. Die einzelnen Schritte, z.B. der Vernetzungs-Prozess (meshing) laufen automatisch im Hintergrund ab.
Da mitunter sehr viel Rechenleistung nötig ist, um die Berechnung durchzuführen, stellen die ersten Unternehmen ihren Nutzern Rechenleistung in Form von Cloud-Diensten zur Verfügung.
Es gibt Prä-/Postprozessoren mit graphischer Oberfläche und getrennten FE-Lösern.
Es gibt Programmframeworks ohne graphische Oberfläche, meist als Präprozessor mit integriertem Gleichungslöser, die per Programmiersprache bedient werden, um beispielsweise mit selbst angefertigten Zusatzroutinen den FE-Solver zu steuern.

Mathematische Herleitung

Das untersuchte Lösungsgebiet wird zunächst in Teilgebiete $G_{1},\dotsc ,G_{m}\subset G$ , die finiten Elemente, eingeteilt:

$G=\bigcup _{{e=1}}^{m}G_{e}$ .

Innerhalb werden für die gesuchte Lösungsfunktion $f\colon G\to \mathbb{R} ^{k}$ nun verschiedene Ansatzfunktionen $\psi _{1},\dotsc ,\psi _{n}\colon G\to \mathbb {R} ^{k}$ definiert, die jeweils nur auf wenigen Elementen ungleich Null sind. Diese Eigenschaft ist der eigentliche Grund für die Bezeichnung „finite“ Elemente.

Durch eine Linearkombination der Ansatzfunktionen werden mögliche Lösungen der numerischen Näherung festgelegt

$f(x)=\sum _{{i=1}}^{n}c_{i}\cdot \psi _{i}(x)\quad ,c_{i}\in \mathbb{R}$ .

Da jede Testfunktion auf den meisten Elementen verschwindet, lässt sich umgekehrt die auf das Element $G_{e}$ eingeschränkte Funktion $f|G_{e}$ , durch die Linearkombination weniger Testfunktionen $\psi _{i}$ beschreiben.

Lassen sich die Differentialgleichungen und die Randbedingungen des Problems als lineare Operatoren bezüglich der Funktionen darstellen, führt dies zu einem linearen Gleichungssystem bzgl. der freien Variablen der Linearkombination $c_{i}$ :

mit

= Lineare Abbildung von $\mathbb {R} ^{n}$ in einen Funktionsraum

= Vektor der Linearkombinationsfaktoren $c_{i}$

= Funktion, die die Differentialgleichung und die Randbedingungen repräsentiert

Um ein endliches lineares Gleichungssystem zu erhalten, wird auch der Wertebereich von über Ansatzfunktionen $\phi _{1},\dotsc ,\phi _{n}$ modelliert. Dann lässt sich über Linearkombinationen der $\phi _{i}$ beschreiben:

$g(x)=\sum _{{i=1}}^{n}b_{i}\cdot \phi _{i}(x)\quad ,b_{i}\in \mathbb{R}$

und man erhält insgesamt das Gleichungssystem

$A\cdot c=b$

mit

= quadratische Matrix mit $A=\phi ^{{-1}}\cdot D$

= Vektor der Linearkombinationsfaktoren $c_{i}$

= Vektor der Linearkombinationsfaktoren $b_{i}$

Die Dimension der Matrix ergibt sich aus der Anzahl der Ansatzfunktionen multipliziert mit den dem physikalischen Modell zugrunde liegenden Freiheitsgraden . Die Dimension der Matrix ist die Anzahl der Gesamtfreiheitsgrade, wobei dem Modell entsprechende Festlegungen bezüglich der Eindeutigkeit des Problems (z.B. im Fall eines elastischen Körpers die Starrkörperverschiebungen) ausgeschlossen werden müssen.

Weil jedes Element nur mit wenigen benachbarten Elementen verbunden ist, sind die meisten Werte der Gesamtmatrix Null, so dass sie „dünnbesetzt“ ist. In den meisten Anwendungsfällen werden die gleichen Funktionen als Ansatzfunktionen $\psi _{i}$ und Testfunktionen $\phi _{i}$ verwendet. In diesem Fall ist die Matrix außerdem symmetrisch zu ihrer Hauptdiagonale.

Ist die Anzahl der Freiheitsgrade nicht allzu groß (bis ca. 500.000), lässt sich dieses Gleichungssystem am effizientesten mittels eines direkten Verfahrens lösen, zum Beispiel mit dem gaußschen Eliminationsverfahren. Hierbei kann die dünnbesetzte Struktur des Gleichungssystems effektiv genutzt werden. Während beim Gauß-Algorithmus der Berechnungsaufwand für Gleichungen ${\mathcal O}(N^{3})$ beträgt, lässt sich der Aufwand durch geschickte Pivotwahl (zum Beispiel Markowitz-Algorithmus oder graphentheoretische Ansätze) aber deutlich reduzieren.

Für mehr als 500.000 Unbekannte bereitet die schlechte Kondition des Gleichungssystems den direkten Lösern zunehmend Schwierigkeiten, so dass man für große Probleme im Allgemeinen iterative Löser, die schrittweise eine Lösung verbessern, verwendet. Einfache Beispiele dafür sind das Jacobi- und Gauß-Seidel-Verfahren, praktisch werden aber eher Mehrgitterverfahren oder vorkonditionierte Krylow-Unterraum-Verfahren, wie das Verfahren der konjugierten Gradienten oder GMRES, verwendet. Aufgrund der Größe der Gleichungssysteme ist manchmal der Einsatz von Parallelrechnern nötig.

Ist die partielle Differentialgleichung nichtlinear, ist auch das resultierende Gleichungssystem nichtlinear. Ein solches lässt sich in der Regel nur über numerische Näherungsverfahren lösen. Ein Beispiel für ein solches Verfahren ist das Newton-Verfahren, in dem schrittweise ein lineares System gelöst wird.

Es gibt heute eine Vielzahl von kommerziellen Computerprogrammen, die nach der Methode der Finiten Elemente arbeiten.

Schwache Formulierung

Eine elliptische partielle Differentialgleichung lässt sich schwach formulieren, d.h. die Problemstellung kann auf eine Art ausgedrückt werden, die von der Lösung weniger Glattheit verlangt. Dies geschieht wie folgt.

Gegeben sei ein Hilbert-Raum , ein Funktional (Funktion aus dessen Dualraum) $f\in H':=\{g:\;H\rightarrow {\mathbb {R}}|\;g{\text{ ist linear und stetig}}\}$ , sowie eine auf stetige und elliptische Bilinearform $a(\cdot ,\cdot )$ , so heißt $u\in H$ Lösung des Variationsproblems, wenn

$a(u,v)=f(v)\ \forall v\in H$ .

Existenz und Eindeutigkeit der Lösung liefert der Darstellungssatz von Fréchet-Riesz (für den Fall, dass die Bilinearform symmetrisch ist) bzw. das Lemma von Lax-Milgram (allgemeiner Fall).

Wir wissen, dass der Raum $L^{2}(\Omega ):=\{f:\Omega \rightarrow {\mathbb {R}}\ |\ \|f\|_{{L^{2}}}<\infty \}$ ein Hilbert-Raum ist. Ausgehend hiervon kann man die Sobolewräume $H^{s}(\Omega )$ über die sogenannte schwache Ableitung definieren.

Das Problem a(u,v)=f(v) kann man als eine Variante einer partiellen Differentialgleichung auf einem Gebiet $\Omega$ auffassen.

Das Poissonproblem als Beispiel:

$-\Delta u(x)=f(x)\ \forall x\in \Omega$

$u(x)=0\ \forall x\in \partial \Omega$

wobei hier $\Delta$ den Laplace-Operator bezeichnet. Eine Multiplikation mit unendlich oft stetig differenzierbaren Funktionen $\psi \in C_{0}^{{\infty }}(\Omega )$ ergibt nach einer Integration

$\Leftrightarrow -\int _{\Omega }\Delta u\,\psi \,dx=\int _{\Omega }f\,\psi \,dx\ \forall \psi \in C_{0}^{\infty }(\Omega ).$

Eine partielle Integration (Erste Greensche Formel) sowie die Nullrandbedingungen für $\psi$ liefern dann

$\Leftrightarrow \int _{\Omega }\nabla u\,\nabla \psi \,dx=\int _{\Omega }f\,\psi \,dx\ \forall \psi \in C_{0}^{\infty }(\Omega )$

Nun ist $a(u,\psi ):=\int _{\Omega }\nabla u\,\nabla \psi \,dx$ eine elliptische und stetige Bilinearform auf $H^1_0(\Omega):=\overline{C^{\infty}_0(\Omega)}^{\|\cdot\|_{H^1}}$ , sowie die rechte Seite $(f,\psi )_{{L^{2}}}=f(\psi )$ eine stetige Linearform auf $H_{0}^{1}(\Omega )$

Besitzt der betrachtete Funktionenraum/Hilbert-Raum eine endliche Basis, so kann man ein lineares Gleichungssystem aus der Variationsformulierung gewinnen.

Für Funktionenräume entscheidet die Wahl der Basis über die Effizienz des Verfahrens. Gängig sind hierbei die Verwendung von Splines mit Triangulierungen, sowie in bestimmten Fällen die diskrete Fourier-Transformation (Aufspaltung in Sinus und Cosinus).

Aufgrund von Flexibilitätsüberlegungen bezüglich der Geometrie des Gebietes $\Omega$ wird in der Regel folgender Ansatz gewählt.

Man diskretisiert das Gebiet $\Omega$ , indem man es in Dreiecke zerteilt und man benutzt Splines $\lambda _{p}(x)$ , assoziiert mit den Eckpunkten p, um den endlichdimensionalen Funktionenraum auf $\Omega$ aufzuspannen. Die Splines erfüllen an festgelegten Punkten auf den Dreiecken $\lambda _{p}(q)=\delta _{{pq}}$ (wobei δ das Kronecker-Delta ist). Damit kann man dann eine diskrete Funktion $u_{h}(x)$ darstellen durch

$u_{h}(x)=\sum _{p}u_{p}\,\lambda _{p}(x)$

mit $u_{p}$ den Koeffizienten bezüglich der Basisdarstellung. Aufgrund der endlichen Basis muss man nicht mehr gegen alle $\psi \in H_{0}^{1}$ testen, sondern nur noch gegen alle Basisfunktionen, die Variationsformulierung reduziert sich aufgrund der Linearität auf

$a(u_{h},\lambda _{q})=\sum _{p}u_{p}\,a(\lambda _{p},\lambda _{q})=(f,\lambda _{q})\ \forall q$

Also haben wir ein lineares Gleichungssystem zum Lösen gewonnen

$A\cdot u=f$ ,

mit

$A_{{pq}}=a(\lambda _{p},\lambda _{q})$

und

$f_{q}=(f,\lambda _{q})$

Dieses Resultat erhält man mit jeder endlichen Basis des Hilbert-Raumes.

Beispiel

Formale Definition des finiten Elementes (nach Ciarlet)

Ein finites Element ist ein Tripel $E=(T,\Pi ,\Sigma )$ , wobei:

ist ein nichtleeres Gebiet (z.B. Dreiecke, Vierecke, Tetraeder usw.)
$\Pi$ ist ein endlichdimensionaler Raum von Ansatzfunktionen (lineare, quadratische oder kubische Formfunktionen, also Splines; Sinus usw.) $\rightarrow$ Formfunktionen
$\Sigma$ ist eine Menge von linear unabhängigen Funktionalen auf $\Pi$ $\rightarrow$ Knotenvariablen

Es gelte für die Funktionale, dass sie zu Funktionen der Basis assoziiert seien:

$\sigma _{i}\in \Sigma :\sigma _{i}(\pi _{j})=\delta _{{ij}}\ \ j\in \{1,\ldots ,\dim(\Pi )\}$

So gilt für jede Funktion

$v\in \Pi :v(x)=\sum _{i}\sigma _{i}(x)\,\pi _{j}$ .

Für Sinus als Basisfunktion im ${\mathbb (}R)^{1}$ ist dann

$\operatorname {span}\,\left\{\sin(x),\sin(2x),\ldots ,\sin(nx)\right\}=\Pi$

und die Funktionale

$\sigma _{i}(\psi ):=(\sin(ix),\psi )_{{L^{2}}}$ .

Für Splines genügt hingegen die Punktauswertung auf den festgelegten Punkten der Dreiecke: $\sigma _{p}(\psi ):=\psi (p)$ .

S1: Lineare Elemente auf Dreiecken

Der FEM-Raum der stetigen, stückweise lineare Funktionen ist definiert als:

$S^{1}(\Omega ,T):=\left\lbrace u\in C^{0}(\Omega )\;|\;{u|}_{\Delta }(x,y)=a_{\Delta }+b_{\Delta }x+c_{\Delta }y\quad \forall \Delta \in T\right\rbrace$ ,

wobei $\Omega \subset \mathbb {R} ^{2}$ ein Gebiet und $T\subset \Omega$ die Triangulierung des Gebietes mit Dreiecken $\Delta$ ist sowie ${u|}_{\Delta }$ die Einschränkung der stetigen Funktion auf das Dreieck $\Delta \in T$ bezeichnet.

P1: Lineares Referenzelement auf einem Dreieck

Das Referenzelement ${\hat {T}}$ ist definiert als:

${\hat {T}}=\left\lbrace (x,y)\in \mathbb {R} ^{2}|0<x<1\;{\text{and}}\;0<y<1-x\right\rbrace .$

Die linearen Elemente sind Funktionen vom Typ:

$p\colon {\hat {T}}\to \mathbb {R} ,(x,y)\mapsto a+bx+cy\quad {\text{mit}}\;a,b,c\in \mathbb {R} .$

Zur Definition der Funktion reicht es dabei schon aus die Werte an den Eckpunkten ${\hat {E_{1}}}:=(0,0),{\hat {E_{2}}}:=(1,0),{\hat {E_{3}}}:=(0,1)$ vorzugeben. Deshalb kann man alle Funktionen mittels Basisfunktionen ${\hat {\phi }}_{i}$ darstellen:

$p(x,y)=\sum _{i=1}^{3}\alpha _{i}{\hat {\phi }}_{i}(x,y)\quad {\text{mit}}\;\alpha _{i}\in \mathbb {R} .$

Die Basisfunktionen sind gegeben als lineare Funktionen, die jeweils nur an einem Eckpunkt ungleich Null sind:

${\hat {\phi }}_{i}({\hat {E_{j}}})=\delta _{i,j}\quad \forall i,j=1,2,3$

wobei $\delta$ die Kronecker-Delta Funktion ist.

Transformation des Referenzelementes

Zur Verknüpfung des Referenzelementes mit einem beliebigen Dreieck (Eckpunkte: $E_{1},E_{2},E_{3}$ ) nutzt man eine lineare Transformation $F_{T}$ :

${\begin{aligned}F_{T}&\colon {\hat {T}}\to T,\quad F(v)=B_{T}v+b_{T}\\B_{T}&:=\left(E_{2}-E_{1},E_{3}-E_{1}\right),\quad b_{T}:=E_{1}.\end{aligned}}$

Bei vielen Aufgaben bezüglich partieller Differentialgleichungen muss das $L^{2}$ -Skalarprodukt von Basisfunktionen $\phi _{i}$ (definiert auf einem beliebigen Dreieck ) berechnet werden:

$(\nabla \phi _{i},\nabla \phi _{j})_{L^{2}(T)}=\int _{T}\nabla \phi _{i}(x)\cdot \nabla \phi _{j}(x)dx.$

Mithilfe des Transformationssatzes kann man die Integration auf das Referenzelement verschieben:

$\int _{T}\nabla \phi _{i}(x)\cdot \nabla \phi _{j}(x)dx=\int _{\hat {T}}\left(\left(DF_{T}(z)\right)^{-T}\nabla {\hat {\phi }}_{i}(z),\left(DF_{T}(z)\right)^{-T}\nabla {\hat {\phi }}_{j}(z)\right)_{\mathbb {R} ^{2}}|\det DF_{T}(z)|dz.$

Literatur

Martin Mayr/Ulrich Thalhofer: Numerische Lösungsverfahren in der Praxis: FEM-BEM-FDM. Hanser, 1993, ISBN 3-446-17061-8.
D. Braess: Finite Elemente – Theorie, schnelle Löser und Anwendungen in der Elastizitätstheorie. 4. Auflage. Springer, 2007, ISBN 978-3-540-72449-0.
Klaus-Jürgen Bathe: Finite-Elemente-Methoden. 2. Auflage. Springer-Verlag, 2002, ISBN 3-540-66806-3.
René de Borst, Mike Crisfield, Joris Remmers, Clemens Verhoosel: Nichtlineare Finite-Elemente-Analyse von Festkörpern und Strukturen, Wiley-VCH 2014

Basierend auf einem Artikel in:

Wikipedia.de