Finite-Elemente-Methode

Die Finite-Elemente-Methode (FEM), auch Methode der finiten Elemente und Finite Element Analysen (FEA) 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 

Mit der FEM können Probleme aus verschiedenen physikalischen Disziplinen (statische, dynamische und nichtlinear thermo-physikalische Probleme) 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 schwache Formulierung des zu lösenden Differentialgleichungsproblems ein, das die physikalischen Gesetze beschreibt, 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 Lösung stellt letztlich die numerische Lösung des betrachteten Differentialgleichungsproblems dar.

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.

Die Geschichte der Finite-Elemente-Methode erschließt sich aus den Forschungen und Veröffentlichungen der folgenden Autoren (Auswahl):
 Karl Heinrich Schellbach: Variationsrechnung; Lösung eines Minimalflächenproblems (1851/52)
 Ernst Gustav Kirsch: Die Fundamentalgleichungen der Theorie der Elastizität fester Körper, hergeleitet aus der Betrachtung eines Systems von Punkten, welche durch elastische Streben verbunden sind (1868)
 John William Strutt, 3. Baron Rayleigh (1842–1919): On the theory of resonance. 1870
 Walter Ritz (1878–1909): neue Methode zur Lösung von Variationsproblemen, Ritz’sches Verfahren (1908/09)
 Boris G. Galerkin (1871–1945): Verfahren der gewichteten Residuen (1915)
 Erich Trefftz (1926): lokal begrenzte Ansatzfunktionen; Gegenstück zum Ritz’schen Verfahren
 Hans Ebner (1929): Schubblech als ebenes Element im Flugzeugbau
 Alexander Hrennikoff (1896–1984): Stabmodelle, Ersetzen von Scheiben durch Fachwerke, Platten durch Trägerroste 1940/41
 Richard Courant (1888–1972): Variational methods for the solution of problems of equilibrium and vibration(s). 1943 (Ansatzfunktionen mit lokalem Träger, elementweise Ansätze für Schwingungsprobleme)
 William Prager (1903–1980), John Lighton Synge (1897–1995): Approximation in Elasticity based on the concept of function space. 1947
 John Argyris (1913–2004): Kraft- und Verschiebungsmethode für Stabtragwerke, Matrizenformulierung (1954/55)
 M. J. Turner, Ray W. Clough, H. C. Martin, L. J. Topp: Stiffness and deflection analysis of complex structures. 1956 (erste Strukturberechnung von Flugzeugflügeln bei Boeing, erste Anwendung der FEM mit Computerprogramm, erste Anwendung von Flächenelementen)
 Ray W. Clough (1920–2016): The finite element method in plane stress analysis. 1960 (wahrscheinlich erste Verwendung des Begriffs Finite Elemente)
 Spierig (1963): Entwicklung von Dreieckelementen, Übertragung auf Schalen
 Olgierd Cecil Zienkiewicz (1921–2009), Pionier der FEM und erstes Standardwerk (Lehrbuch): The Finite Element Method in Structural and Continuum Mechanics, 1967 (mit Y. K. Cheung)
 Alfred Zimmer (* 1920) und Peter Groth (* 1938), Pioniere der FEM, erstes deutsches FEM-Lehrbuch: Elementmethode der Elastostatik, 1969 Oldenbourg Verlag München, Wien
 Olga Alexandrowna Ladyschenskaja (1922–2004), Ivo Babuška (* 1926) und Franco Brezzi (* 1945) – Ladyschenskaja-Babuška-Brezzi-Bedingung für die Stabilität eines gemischten Finite-Elemente-Problems mit Sattelpunktstruktur
 Ivo Babuška (1926–2023) – adaptive Finite-Elemente-Algorithmen

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.

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 Diskretisierungsverfahren. Das bedeutet, dass für ein unendlichdimensionales Problem eine endlichdimensionale
Näherung
berechnet wird.

Hierzu wird das Grundgebiet 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 endlichdimensionale Approximationen eingeführt, welche
die unbekannten Knotengrößen als Parameter enthalten. Die lokalen Approximationen werden dann in die schwache Formulierung des Randwertproblems eingeführt. Die dabei entstehenden Elementintegrale
werden i.a. 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.  Auf diese Weise werden Randwertprobleme für lineare partielle Differentialgleichungen in
ein lineares Gleichungssystem  ü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.

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.

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 

Die häufigsten Randbedingungen für partielle Differentialgleichungen sind Randbedingungen erster Art (Dirichlet-Bedingungen), zweiter Art (Neumann-Bedingungen) und dritter Art (Robin-Bedingungen). Bei Dirichlet-Bedingungen werden am Rand Funktionswerte vorgeschrieben, bei Neumann-Bedingungen Normalableitungen und bei Robin-Bedingungen Linearkombinationen davon.
 Je nach Art des physikalischen Problems kann es sich um verschiedene physikalische Größen handeln, wie in der Tabelle beispielhaft dargestellt.

Weitere Arten von Randbedingungen 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.

Nur bei Dirichlet-Bedingungen ist die Diskretisierung einfach, hier werden in den Knoten auf dem Rand die gegebenen Randwerte einfach übernommen. Bei homogenen Neumann-Bedingungen hat man gar nichts zu tun, dies sind sogenannte natürliche Randbedingungen. Bei inhomogenen Randbedingungen zweiter Art oder
Randbedingungen dritter Art hingegen beeinflussen die Randbedingungen die schwachen Formulierungen des Problems!

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:

 und  sind Massen-, Dämpfungs- und Steifigkeitsmatrix des Systems;  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  auf ein System mit Dämpfung,
dann wird sich nach ausreichend langer Zeit auch ein konstanter Verschiebungsvektor  einstellen.
Für  ist dann  , und das Gleichungssystem vereinfacht sich zu  mit der Lösung 

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  ausgesetzt,
dann kann es erforderlich sein, die Eigenfrequenzen des Systems zu ermitteln, um Resonanzen im Betrieb zu vermeiden.

Eigenfrequenzen sind alle Frequenzen , für die ein Verschiebungsvektor
 eine Lösung des unbelasteten () und ungedämpften
() Gleichungssystems darstellt. Für den Geschwindigkeits- und Beschleunigungsvektor gilt dann

und das Gleichungssystem lautet damit

Um die Eigenfrequenzen  und die dazugehörigen Eigenformen  zu berechnen, muss der implizite Solver also das Eigenwertproblem

lösen.

Explizite FEM-Solver

Explizite Solver berechnen die Verschiebungsvektoren  zu bestimmten diskreten Zeitpunkten 
innerhalb eines vorgegebenen Zeitintervalls. Knotengeschwindigkeiten und -beschleunigungen werden durch Differenzenquotienten
aus den Verschiebungen  zu aufeinanderfolgenden Zeitpunkten
 angenähert. Mit konstanter Zeitschrittweite  gilt

hat das diskretisierte Gleichungssystem die Form

Durch Auflösen dieser Gleichung erhält man eine Beziehung, mit der der Verschiebungsvektor 
aus den vorher berechneten Vektoren  und  ermittelt werden kann:

Die Berechnung der Inversen  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.

Der mathematische Zugang zur Finite-Elemente-Methode

Die schwache Formulierung einer Randwertaufgabe und das Galerkin-Verfahren 

Betrachtet wird  eine elliptische Randwertaufgabe zweiter Ordnung, als Beispiel die erste Randwertaufgabe für eine Poisson-Gleichung in einem zweidimensionalen Gebiet  mit Rand :

Besitzt  Ecken, so sagt die Lösungstheorie elliptischer Gleichungen entgegen der Erwartung, dass dieses Problem gar keine zweimal differenzierbare Lösung besitzt! Deshalb ist die Nutzung einer schwachen Formulierung wichtig, extrem wichtig auch hinsichtlich von Näherungsverfahren zur Berechnung von Lösungen wie der Finiten-Element-Methode (FEM), die nur stetige, aber keine differenzierbaren Ansatzfunktionen nutzen (nur für Probleme höherer, z. B. vierter Ordnung, werden hin und wieder differenzierbare Ansatzfunktionen verwendet).

Zur Herleitung der schwachen Formulierung multipliziert man die gegebene Gleichung mit einer Funktion , integriert über  und beseitigt dann die zweiten Ableitungen durch partielle Integration (bzw. der Anwendung eines Integralsatzes). Das ergibt

Nun muss man sich fragen, für welche  diese Gleichung Sinn macht. Wichtig ist offenbar, dass  und  und deren Ableitungen quadratisch integrierbar sind, denn dann existieren alle vorkommenden Integrale. Präzise sagt man:  und  seien aus dem Sobolev-Raum . Die Menge der Funktionen aus dem , die zudem auf dem Rand von  gleich Null sind, nennt man .

Damit ergibt sich folgende schwache Formulierung:  genügt

Solch eine Gleichung nennt man auch Variationsgleichung. Fundamental ist, dass sich bei anderen Randbedingungen die schwache Formulierung ändert.
Hat man etwa die Randbedingung 2. Art  auf  , so entsteht mit  die schwache Formulierung

Ist aber die Randbedingung 3. Art  gegeben, so bleibt es bei , aber es entsteht

Das bedeutet später für die Finite-Elemente-Methode: Randbedingungen 1. Art sind unmittelbar zu berücksichtigen, Randbedingungen 2. oder 3. Art sind natürliche Randbedingungen und werden über die adäquate schwache Formulierung indirekt eingebaut.

Jede dieser drei schwachen Formulierungen bzw. Variationsgleichungen besitzt folgende allgemeine Form: Gesucht ist , einem Teilraum vom Raum 
mit

 für alle .

Hierbei ist  eine Linearform,  nennt man Bilinearform.

Da Randwertaufgaben für partielle Differentialgleichungen nur in Spezialfällen exakt lösbar sind, sucht man nun Näherungsverfahren. Ein wichtiges solches Verfahren ist das Galerkin-Verfahren für Variationsgleichungen.

Dazu wählt man einen endlichdimensionalen Teilraum  von  und nennt  Galerkin-Näherung, wenn gilt

 für alle .

Dies ist äquivalent zu einem linearen Gleichungssystem mit endlich vielen Unbekannten, und damit ist die Näherungslösung  berechenbar. Dies sieht man folgendermaßen: Ein linearer endlichdimensionaler Raum besitzt eine endliche Basis, diese Basis von  sei  für
. Dann gibt es Konstanten  mit

Setzt man dies in die Diskretisierung ein sowie  für , so erhält man das lineare Gleichungssystem

 für  

zur Berechnung der unbekannten Koefficienten  der Näherungslösung. Aus den gegebenen Daten und den gewählten Basisfunktionen des endlichdimensionalen Raumes sind die Elemente der Koeffizientenmatrix und die ,rechte' Seite des Gleichungssystem berechenbar, gegebenenfalls durch numerische Berechnung der auftretenden Integrale. Im Fall einer symmetrischen Bilinearform ist die Koeffizientenmatrix symmetrisch. Unter gewissen Voraussetzungen an die Bilinearform besitzen sowohl das Variationsproblem als auch das diskrete Problem (das erzeugte Gleichungssystem) eine eindeutige Lösung.

Hingewiesen sei noch darauf, dass im Fall einer symmetrischen Bilinearform man das Variationsproblem

mit dem Ritz-Verfahren

näherungsweise lösen kann, die notwendige Optimalitätsbedingung für dieses Optimierungsproblem aber wieder zu der Galerkin-Formulierung
 mit

 für alle 

führt. Während im symmetrischen Fall also egal ist, ob man mit dem Ritz- oder dem Galerkin-Verfahren diskretisiert, steht im nichsymmetrischen Fall nur der Galerkin-Zugang zur Verfügung.

Entscheidend für die Realisierung einer Galerkin-Diskretisierung  (oder auch Ritz) ist die Wahl des endlichdimensionalen Teilraumes  von
. Bei der spektralen Galerkin-Methode wählt man Polynome und speziell orthogonale Polynome als Basisfunktionen. Bei der Finite-Elemente Methode
sind die Basisfunktionen Splines, d. h., stückweise polynomiale Funktionen. Da man die finiten Elemente lokal definiert, muss man sorgfältig darauf achten, ob der sich global ergebene Finite-Elemente-Raum  die Eigenschaft  besitzt. Wenn ja, so heißt die resultierende Finite-Elemente-Methode konform.

Finite Elemente und Finite-Elemente-Räume 

Ausgangspunkt ist die Zerlegung des Grundgebietes  in endlich viele einfache Teilgebiete. Das sind im zweidimensionalen Fall Dreiecke und Vierecke, im dreidimensionalen Tetraeder und Parallelepipede. Der Einfachheit halber sei hier  zweidimensional und ein polygonales Gebiet. Zerlegt wird in Dreiecke.

Ein Dreieck bzw. Element der Zerlegung sei , einschließlich des Randes . Die Zerlegung sei grundsätzlich zulässig, d. h.,  für  sei  leer oder ein Eckpunkt (Knoten) oder eine gemeinsame Kante der beiden Dreiecke. Ausgeschlossen wird also bei zwei benachbarten Dreiecken, dass eine Ecke eines Dreiecks auf einer Kante des Nachbardreiecks liegt, aber keine Ecke desselben ist.

Ein finites Element besteht nun aus einem Dreieck , einer linearen Menge  von Polynomen auf  und einer Menge von Vorgaben (oder linearen Funktionalen), die ein Polynom aus  eindeutig festlegen. Die Anzahl dieser Vorgaben ist also gleich der Dimension von . Diese Vorgaben sind oft Werte in Punkten von  (nicht nur in den Knoten bzw. Eckpunkten), aber auch Integrale über Kanten von  oder Integrale über . Sind nur Funktionswerte vorgegeben, spricht man von Lagrange-Elementen.
Strebt man (für Probleme vierter Ordnung) global stetig differenzierbare Elemente an, sind Vorgaben auch Werte von Ableitungen in gewissen Punkten. Solche Elemente heißen Hermite-Elemente.

Da ein finites Element also zunächst nur lokal definiert ist, erhält man ein Element des Finiten-Elemente-Raumes  nun einfach durch stückweises Zusammensetzen. Für die Galerkin-Diskretisierung von Randwertaufgaben zweiter Ordnung ist nun extrem wichtig, wann -Konformität vorliegt. Zentral ist folgende Aussage:

Ist  stückweise polynomial und global stetig, so gilt .

Das bedeutet also, dass man bei der zunächst lokalen Konstruktion nur einen stetigen Übergang von Element zu Element sichern muss, die Vorgaben sind dahingehend zu überprüfen. Benötigt man , so ist zusätzlich dafür zu sorgen, dass  auf dem Rand des Gebietes gleich Null ist. Bei Randbedingungen 2. oder 3. Art ist die korrekte schwache Formulierung für das Einarbeiten der Randbedingungen zuständig.

Ein grundlegendes Beispiel sind zwei Elemente, für die  ein Dreieck ist und  die Menge aller Polynome ersten Grades auf , also die Menge aller Funktionen der Form . Die Werte in drei nicht auf einer Geraden liegenden Punkte bestimmen eine Funktion aus  eindeutig. Wählt man die drei Funktionswerte in den Ecken, so erhält man globale Stetigkeit, denn die Werte in den zwei Ecken einer Kante definieren eine lineare Funktion einer Variablen eindeutig. Wählt man jedoch die drei Funktionswerte in den Seitenmitten, so erhält man ein unstetiges Element. Das erste beschriebene Element heißt konformes -Element, das zweite ist das nichtkonforme -Element.

Es sei nun  die Menge aller Polynome vom Grad . Die Dimension von  ist . Man benötigt also Werte in  Punkten, um das konforme -Element definieren. Dazu und für die Beschreibung von lokalen Basisfunktionen von  sind baryzentrische Koordinaten nützlich.

Ist  ein Dreieck mit den Ecken , dann sind die baryzentrischen Koordinaten  eindeutig durch die beiden Gleichungen  und

den Punkten  zugeordnet. Die Ecke  z. B. besitzt die baryzentrischen Koordinaten , der Schwerpunkt des Dreiecks die Koordinaten .

Für das konforme -Element sind nun lokale Basisfunktionen leicht angebbar, diese sind  und es gilt mit der Kronecker-Delta Funktion 

.

Als nächstes wird die Frage beantwortet, in welchen Punkten eines Dreiecks man die Funktionswerte für das -Element und das -Element vorschreibt. Die Menge der quadratischen Funktionen besitzt die Dimension 6. Deshalb wählt man die drei Eckpunkte und die drei Seitenmitten, diese besitzen die baryzentrischen Koordinaten ,  und .
Die Menge der kubischen Polynome besitzt die Dimension 10. Die Menge der Vorgabepunkte für die Funktionswerte für das  -Element  sind die drei Ecken, die zwei Punkte auf jeder Kante, die die Kante in drei gleiche Teile zerlegen und der Schwerpunkt des Dreiecks.

Die lokalen Basisfunktionen für das  -Element sind

Jede der lokalen Basisfunktionen besitzt folgende Eigenschaft: sie ist gleich Eins in genau einem der 6 Vorgabepunkte und gleich Null in den anderen fünf.
Analog konstruiert man eine lokale Basis für  -Elemente.

Sei   einer der Vorgabepunkte bei einem - Finite-Elemente-Raum. Dann besitzt eine globale Ansatzfunktion die folgende Eigenschaft: sie ist gleich Eins im Punkt  und gleich Null in allen anderen Vorgabepunkten. Dies impliziert, dass der Träger einer globalen Ansatzfunktion (der Bereich, in dem die Ansatzfunktion ungleich Null ist) nur aus wenigen Dreiecken der Zerlegung besteht! Damit entstehen in der Koeffizientenmatrix des Gleichungssystems, das mittels der FEM erzeugt wird, viele Nullen, die Matrix ist schwach besetzt. Dies ist ein Vorteil von Splines als Ansatzfunktionen beim Galerkin-Verfahren.

Bei sogenannten affinen Familien von finiten Elementen, und dazu gehören die -Elemente, werden die globalen Ansatzfunktionen zur Erzeugung des diskreten Problems (Gleichungssystem) dadurch nicht benötigt, dass man jedes Element durch eine lineare Abbildung auf ein Referenzelement abbildet (dadurch werden auch die Ansatzfunktionen und Vorgabewerte auf   auf ebensolche auf dem Referenzelement abgebildet). Dies wird im Folgenden für lineare Elemente demonstriert.

Die Erzeugung des diskreten Problems durch elementweises Herangehen und Transformation 

Betrachtet wird exemplarisch das obige Poissonproblem mit  und ihre FEM-Diskretisierung mit linearen finiten Elementen. Es sei

mit den globalen nodalen Ansatzfunktionen, die genau in einem der inneren Eckpunkte der Triangulation gleich Eins sind, in allen anderen Eckpunkten, insbesondere in denen auf dem Rand des Gebietes, gleich Null. Damit erfüllt der Ansatz die homogenen Randbedingungen. Ist  der -te der   inneren Eckpunkt der Triangulation, so gilt .

Die unbekannten   berechnet man aus dem Gleichungssystem

Dabei ist   die sogenannte Steifigkeitsmatrix mit den Elementen

und der rechten Seite mit   und  .

Der entscheidende Trick ist nun, dass man die Integrale über die einzelnen Dreiecke  der Triangulation berechnet und dann über alle Beiträge der Elemente summiert, zudem nicht wirklich über  integriert, sondern durch Transformation die Integration über ein Referenzelement
realisiert.

Dazu erklärt man zu einem Dreieck bzw. Element  gehörende Elementsteifigkeitsmatrizen durch

mit

Analog wird durch

eine elementweise rechte Seite erklärt. Durch Summation über alle  entsteht das final zu lösende Gleichungssystem.

In dem betrachteten Spezialfall kann man die Elementsteifigkeitsmatrizen exakt berechnen, andernfalls ist eine numerische Integration erforderlich.
Es sei  ein Element mit den Ecken  und z. B. . Dann kann man  durch die Transformation

auf das Referenzelement mit den Ecken   in der  -Ebene abbilden. Die Integrale zur Berechnung der Elementsteifigkeitsmatrix werden durch diese Transformation auf dem Referenzelement berechnet. Bei der Transformation werden aus den drei nodalen lokalen Basisfunktionen auf   in der -Ebene die Basisfunktionen

Die Berechnung deren Ableitung nach  und  ist trivial. Die Funktionaldeterminante  der Transformation ist

der Betrag von  ist gleich dem doppelten Flächeninhalt von . Letztlich muss man noch die Ableitungen nach   und  in Ableitungen nach    und   umrechnen. Es gilt

Insgesamt erhält man dann z. B. für die erste Zeile der Elementsteifigkeitsmatrix

Zurückkehrend zu der Bezeichnung  mit den Ecken ,  und 
ergibt sich für die Elementsteifigkeitsmatrix mit 

und für die Diagonalelemente

Zur Lösung des diskreten Problems 

Die Diskretisierung eines linearen Randwertproblems mit der FEM führt auf ein lineares Gleichungssystem großer Dimension.

Bei nicht zu großer Dimension löst man lineare Gleichungssysteme mit einem direkten Verfahren, Standard ist ein Gaußsches Eliminationsverfahren und seine Varianten. Bei extrem großer Dimension greift man zu einem iterativen Verfahren. Zu den sogenannten Krylow-Unterraum-Verfahren gehören das bekannte konjugierte Gradientenverfahren für symmetrische Probleme und GMRES für nichtsymmetrische Systeme, beide Verfahren werden oft mit Vorkonditionierung genutzt. Da bei der Methode der finiten Elemente durch die notwendige Gitterverfeinerung zur Erreichung der gewünschten Genauigkeit ohnehin mehrere Gitter verwendet werden, bieten sich Mehrgitterverfahren zur Lösung der diskreten Probleme an.

Fehlerabschätzung 

Üblich sind Fehlerabschätzungen für die Finite-Element-Methode im Fall elliptischer Randwertaufgaben zweiter Ordnung in der -Norm, dabei werden Methoden der Funktionalanalysis angewandt.

Ein typisches Ergebnis wird nun skizziert. Die verwendete -Norm ist definiert durch

Betrachtet werden lineare finite Elemente, ein polygonales Gebiet  und eine zulässige Zerlegung in Dreiecke. Es sei
 die längste Seite eines Dreiecks  und . Dann erhält man unter gewissen Voraussetzungen an die der elliptischen Randwertaufgabe zugeordneten Bilinearform und für  die Fehlerabschätzung

wenn zudem die Dreieckszerlegung z. B. der Minimalwinkelbedingung genügt. Das bedeutet: der minimale Innenwinkel aller Dreiecke auf der Familie der betrachteten Zerlegungen ist (auch bei Verfeinerung der Zerlegung) nach unten beschränkt, zu spitz dürfen also Dreiecke bei der Minimalwinkelbedingung nicht werden.
Hingewiesen sei darauf, dass für konvexe Gebiete die Voraussetzung  realistisch ist, für Gebiete mit einspringenden Ecken dagegen nicht.

Für -Elemente hofft man auf eine höhere Fehlerordnung, d. h.

 

Dies ist aber nur zu erwarten, wenn die Lösung der gegebenen Randwertaufgabe zusätzliche Glattheit aufweist.

Erweiterungen der Finite-Elemente-Methode 

Gemischte finite Elemente nutzt man vorwiegend bei Problemen mit Nebenbedingungen, z. B. bei der Divergenzfreiheit im Stokes-Problem.
Nichtkonforme finite Elemente spielen z. B. bei Problemen 4. Ordnung wie der Plattengleichung eine Rolle, wenn man dieselben Elemente wie für Gleichungen 2. Ordnung verwenden will. Die diskontinuierliche Galerkin-Methode ist ein sehr variables Werkzeug zur Anwendung unstetiger Elemente.

Stabilisierte Finite-Element-Methoden wie die Stromliniendiffusion-Finite-Element-Methode benötigt man für Probleme mit dominanter Konvektion. Besitzt man Informationen über existierende Grenzschichten (s. auch Singuläre Störungen), kann man zudem Grenzschichtangepasste Gitter einsetzen.

Literatur 
 
 
 
 
 Band 1: Grundlagen: Basiswissen und Arbeitsbeispiele zu FEMAnwendungen. 2007, ISBN 978-3-8169-2685-6.
 Band 2: Strukturdynamik. 2008, ISBN 978-3-8169-2842-3.
 Band 3: Temperaturfelder. 2007, ISBN 978-3-8169-2714-3.
 Band 4: Elektrotechnik. 2009, ISBN 978-3-8169-2841-6.

 
 
 
 
 
 
 René de Borst, Mike Crisfield, Joris Remmers, Clemens Verhoosel: Nichtlineare Finite-Elemente-Analyse von Festkörpern und Strukturen, Wiley-VCH 2014
 Karl-Eugen Kurrer: The History of the Theory of Structures. Searching for Equilibrium, Ernst und Sohn 2018, S. 881–914, ISBN 978-3-433-03229-9

Weblinks 

 Dr.-Ing. Ulrich Simon: Numerische Methoden in der Biomechanik
 
 F.U. Mathiak:  [PDF; 3,86 MB].

Einzelnachweise 

Numerische Mathematik
Computer Aided Engineering