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.

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.

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.

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:

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

${\displaystyle M,D}$ und ${\displaystyle K}$ sind Massen-, Dämpfungs- und Steifigkeitsmatrix des Systems; ${\displaystyle P(t)}$ ist der Vektor der externen Kräfte, die auf das Modell wirken. ${\displaystyle u}$ 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 ${\displaystyle u}$ gültig ist. Wirkt z. B. eine zeitlich unveränderliche Last ${\displaystyle P(t)=P=const.}$ auf ein System mit Dämpfung, dann wird sich nach ausreichend langer Zeit auch ein konstanter Verschiebungsvektor ${\displaystyle u(t)=u=const.}$ einstellen. Für ${\displaystyle t\to \infty }$ ist dann ${\displaystyle {\ddot {u}}(t)={\dot {u}}(t)=0}$ , und das Gleichungssystem vereinfacht sich zu ${\displaystyle Ku=P}$ mit der Lösung ${\displaystyle u=K^{-1}P}$

Für einen gegebenen Lastvektor ${\displaystyle P}$ kann der Verschiebungsvektor ${\displaystyle u}$ mit Hilfe des Gauß-Algorithmus oder durch QR-Zerlegung von ${\displaystyle K}$ berechnet werden.

Ist ein mechanisches System einer harmonischen Anregung ${\displaystyle 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 ${\displaystyle \omega _{i}}$, für die ein Verschiebungsvektor ${\displaystyle u(t)={\hat {u}}_{i}\sin(\omega _{i}t)}$ eine Lösung des unbelasteten (${\displaystyle P(t)=0}$) und ungedämpften (${\displaystyle D=0}$) Gleichungssystems darstellt. Für den Geschwindigkeits- und Beschleunigungsvektor gilt dann

${\displaystyle {\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

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

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

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

lösen.

Explizite FEM-Solver

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

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

${\displaystyle {\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

${\displaystyle 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 ${\displaystyle u_{i+1}}$ aus den vorher berechneten Vektoren ${\displaystyle u_{i}}$ und ${\displaystyle u_{i-1}}$ ermittelt werden kann:

${\displaystyle 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 ${\displaystyle M^{-1}}$ wird in der Praxis nicht durchgeführt, da explizite Solver in der Regel ${\displaystyle M}$ als Diagonalmatrix annehmen und daher jede Zeile des Gleichungssystems nur durch den Diagonaleintrag in der entsprechenden Zeile von ${\displaystyle M}$ 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.

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 ${\displaystyle \Omega }$ mit Rand ${\displaystyle \partial \Omega }$:

${\displaystyle -\triangle u+cu=f\quad {\rm {in}}\,\,\Omega ,\,\,u=0\,\,{\rm {auf}}\,\,\partial \Omega =\Gamma .}$

Besitzt ${\displaystyle \Omega }$ 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 ${\displaystyle v}$, integriert über ${\displaystyle \Omega }$ und beseitigt dann die zweiten Ableitungen durch partielle Integration (bzw. der Anwendung eines Integralsatzes). Das ergibt

${\displaystyle \int _{\Omega }\nabla u\nabla v-\int _{\Gamma }{\frac {\partial u}{\partial n}}v+\int _{\Omega }cuv=\int _{\Omega }fv.}$

Nun muss man sich fragen, für welche ${\displaystyle u,v}$ diese Gleichung Sinn macht. Wichtig ist offenbar, dass ${\displaystyle u}$ und ${\displaystyle v}$ und deren Ableitungen quadratisch integrierbar sind, denn dann existieren alle vorkommenden Integrale. Präzise sagt man: ${\displaystyle u}$ und ${\displaystyle v}$ seien aus dem Sobolev-Raum ${\displaystyle H^{1}(\Omega )}$. Die Menge der Funktionen aus dem ${\displaystyle H^{1}(\Omega )}$, die zudem auf dem Rand von ${\displaystyle \Omega }$ gleich Null sind, nennt man ${\displaystyle H_{0}^{1}(\Omega )}$.

Damit ergibt sich folgende schwache Formulierung: ${\displaystyle u\in H_{0}^{1}(\Omega )}$ genügt

${\displaystyle \int _{\Omega }\nabla u\nabla v+\int _{\Omega }cuv=\int _{\Omega }fv\quad \forall v\in H_{0}^{1}(\Omega ).}$

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 ${\displaystyle {\frac {\partial u}{\partial n}}=0}$ auf ${\displaystyle \Gamma }$, so entsteht mit ${\displaystyle u\in H^{1}(\Omega )}$ die schwache Formulierung

${\displaystyle \int _{\Omega }\nabla u\nabla v+\int _{\Omega }cuv=\int _{\Omega }fv\quad \forall v\in H^{1}(\Omega ).}$

Ist aber die Randbedingung 3. Art ${\displaystyle {\frac {\partial u}{\partial n}}+pu=0}$ gegeben, so bleibt es bei ${\displaystyle u\in H^{1}(\Omega )}$, aber es entsteht

${\displaystyle \int _{\Omega }\nabla u\nabla v+\int _{\Gamma }puv+\int _{\Omega }cuv=\int _{\Omega }fv\quad \forall v\in H^{1}(\Omega ).}$

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 ${\displaystyle u\in V}$, einem Teilraum vom Raum ${\displaystyle H^{1}(\Omega )}$ mit

${\displaystyle a(u,v)=(f,v)}$ für alle ${\displaystyle v\in V}$.

Hierbei ist ${\displaystyle (f,v)=\int _{\Omega }fv}$ eine Linearform, ${\displaystyle a(\cdot ,\cdot )}$ 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 ${\displaystyle V_{h}\subset V}$ von ${\displaystyle V}$ und nennt ${\displaystyle u_{h}\in V_{h}}$ Galerkin-Näherung, wenn gilt

${\displaystyle a(u_{h},v_{h})=(f,v_{h})}$ für alle ${\displaystyle v_{h}\in V_{h}}$.

Dies ist äquivalent zu einem linearen Gleichungssystem mit endlich vielen Unbekannten, und damit ist die Näherungslösung ${\displaystyle u_{h}}$ berechenbar. Dies sieht man folgendermaßen: Ein linearer endlichdimensionaler Raum besitzt eine endliche Basis, diese Basis von ${\displaystyle V_{h}}$ sei ${\displaystyle {\phi _{i}}}$ für ${\displaystyle i=1,\cdots ,N}$. Dann gibt es Konstanten ${\displaystyle c_{j}}$ mit

${\displaystyle u_{h}=\sum _{j=1}^{N}c_{j}\phi _{j}\,\,.}$

Setzt man dies in die Diskretisierung ein sowie ${\displaystyle v_{h}:=\phi _{i}}$ für ${\displaystyle i=1,\cdots ,N}$, so erhält man das lineare Gleichungssystem

${\displaystyle \sum _{j=1}^{N}a(\phi _{j},\phi _{i})c_{j}=(f,\phi _{i})}$ für ${\displaystyle i=1,\cdots ,N}$

zur Berechnung der unbekannten Koefficienten ${\displaystyle c_{j}}$ 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

${\displaystyle \min _{v\in V}[{\frac {1}{2}}a(v,v)-(f,v)]}$

mit dem Ritz-Verfahren

${\displaystyle \min _{v_{h}\in V_{h}}[{\frac {1}{2}}a(v,v)-(f,v)]}$

näherungsweise lösen kann, die notwendige Optimalitätsbedingung für dieses Optimierungsproblem aber wieder zu der Galerkin-Formulierung ${\displaystyle u_{h}\in V_{h}}$ mit

${\displaystyle a(u_{h},v_{h})=(f,v_{h})}$ für alle ${\displaystyle v_{h}\in V_{h}}$

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 ${\displaystyle V_{h}}$ von ${\displaystyle V}$. 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 ${\displaystyle V_{h}}$ die Eigenschaft ${\displaystyle V_{h}\subset V}$ 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 ${\displaystyle \Omega }$ in endlich viele einfache Teilgebiete. Das sind im zweidimensionalen Fall Dreiecke und Vierecke, im dreidimensionalen Tetraeder und Parallelepipede. Der Einfachheit halber sei hier ${\displaystyle \Omega }$ zweidimensional und ein polygonales Gebiet. Zerlegt wird in Dreiecke.

Ein Dreieck bzw. Element der Zerlegung sei ${\displaystyle K_{i}}$, einschließlich des Randes ${\displaystyle {\overline {K}}_{i}}$. Die Zerlegung sei grundsätzlich zulässig, d. h., für ${\displaystyle i\not =j}$ sei ${\displaystyle {\overline {K}}_{i}\cap {\overline {K}}_{j}}$ 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 ${\displaystyle K}$, einer linearen Menge ${\displaystyle P}$ von Polynomen auf ${\displaystyle K}$ und einer Menge von Vorgaben (oder linearen Funktionalen), die ein Polynom aus ${\displaystyle P}$ eindeutig festlegen. Die Anzahl dieser Vorgaben ist also gleich der Dimension von ${\displaystyle P}$. Diese Vorgaben sind oft Werte in Punkten von ${\displaystyle {\overline {K}}}$ (nicht nur in den Knoten bzw. Eckpunkten), aber auch Integrale über Kanten von ${\displaystyle K}$ oder Integrale über ${\displaystyle K}$. 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 ${\displaystyle V_{h}}$ nun einfach durch stückweises Zusammensetzen. Für die Galerkin-Diskretisierung von Randwertaufgaben zweiter Ordnung ist nun extrem wichtig, wann ${\displaystyle H^{1}}$-Konformität vorliegt. Zentral ist folgende Aussage:

Ist ${\displaystyle v_{h}\in V_{h}}$ stückweise polynomial und global stetig, so gilt ${\displaystyle v_{h}\in H^{1}(\Omega )}$.

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 ${\displaystyle v_{h}\in H_{0}^{1}(\Omega )}$, so ist zusätzlich dafür zu sorgen, dass ${\displaystyle v_{h}}$ 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 ${\displaystyle K}$ ein Dreieck ist und ${\displaystyle P}$ die Menge aller Polynome ersten Grades auf ${\displaystyle K}$, also die Menge aller Funktionen der Form ${\displaystyle \alpha x+\beta y+\gamma }$. Die Werte in drei nicht auf einer Geraden liegenden Punkte bestimmen eine Funktion aus ${\displaystyle P}$ 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 ${\displaystyle P_{1}}$-Element, das zweite ist das nichtkonforme ${\displaystyle P_{1}}$-Element.

Es sei nun ${\displaystyle P}$ die Menge aller Polynome vom Grad ${\displaystyle k}$. Die Dimension von ${\displaystyle P}$ ist ${\displaystyle (k+2)(k+1)/2}$. Man benötigt also Werte in ${\displaystyle (k+2)(k+1)/2}$ Punkten, um das konforme ${\displaystyle P_{k}}$-Element definieren. Dazu und für die Beschreibung von lokalen Basisfunktionen von ${\displaystyle P_{k}}$ sind baryzentrische Koordinaten nützlich.

Ist ${\displaystyle K}$ ein Dreieck mit den Ecken ${\displaystyle p^{1},p^{2},p^{3}}$, dann sind die baryzentrischen Koordinaten ${\displaystyle \lambda _{1},\lambda _{2},\lambda _{3}}$ eindeutig durch die beiden Gleichungen ${\displaystyle \lambda _{1}+\lambda _{2}+\lambda _{3}=1}$ und

${\displaystyle x=\sum _{i=1}^{3}\lambda _{i}p^{i}}$

den Punkten ${\displaystyle x\in {\overline {K}}}$ zugeordnet. Die Ecke ${\displaystyle p^{1}}$ z. B. besitzt die baryzentrischen Koordinaten ${\displaystyle (1,0,0)}$, der Schwerpunkt des Dreiecks die Koordinaten ${\displaystyle (1/3,1/3,1/3)}$.

Für das konforme ${\displaystyle P_{1}}$-Element sind nun lokale Basisfunktionen leicht angebbar, diese sind ${\displaystyle \lambda _{1},\lambda _{2},\lambda _{3}}$ und es gilt mit der Kronecker-Delta Funktion ${\displaystyle \delta }$

${\displaystyle \lambda _{i}(p^{j})=\delta _{ij}}$.

Als nächstes wird die Frage beantwortet, in welchen Punkten eines Dreiecks man die Funktionswerte für das ${\displaystyle P_{2}}$-Element und das ${\displaystyle P_{3}}$-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 ${\displaystyle (1/2,0,0)}$, ${\displaystyle (0,1/2,0)}$ und ${\displaystyle (0,0,1/2)}$. Die Menge der kubischen Polynome besitzt die Dimension 10. Die Menge der Vorgabepunkte für die Funktionswerte für das ${\displaystyle P_{3}}$-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 ${\displaystyle P_{2}}$-Element sind

${\displaystyle \lambda _{1}(2\lambda _{1}-1),\lambda _{2}(2\lambda _{2}-1),\lambda _{3}(2\lambda _{3}-1),\lambda _{1}\lambda _{2},\lambda _{1}\lambda _{3},\lambda _{2}\lambda _{3}.}$

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 ${\displaystyle P_{3}}$-Elemente.

Sei ${\displaystyle q}$ einer der Vorgabepunkte bei einem ${\displaystyle P_{k}}$- Finite-Elemente-Raum. Dann besitzt eine globale Ansatzfunktion die folgende Eigenschaft: sie ist gleich Eins im Punkt ${\displaystyle q}$ 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 ${\displaystyle P_{k}}$-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 ${\displaystyle K}$ 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 ${\displaystyle c=0}$ und ihre FEM-Diskretisierung mit linearen finiten Elementen. Es sei

${\displaystyle u_{h}=\sum _{j=1}^{N}c_{j}\phi _{j}}$

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 ${\displaystyle p^{r}}$ der ${\displaystyle r}$-te der ${\displaystyle N}$ inneren Eckpunkt der Triangulation, so gilt ${\displaystyle c_{r}=u_{h}(p^{r})}$.

Die unbekannten ${\displaystyle c_{j}}$ berechnet man aus dem Gleichungssystem

${\displaystyle A_{h}c=f^{h}.}$

Dabei ist ${\displaystyle A_{h}}$ die sogenannte Steifigkeitsmatrix mit den Elementen

${\displaystyle a_{ik}=\int _{\Omega }\nabla \phi _{i}\nabla \phi _{k}}$

und der rechten Seite mit ${\displaystyle f^{h}=(f_{i}^{h})}$ und ${\displaystyle f_{i}^{h}=\int _{\Omega }f\phi _{i}}$.

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

Dazu erklärt man zu einem Dreieck bzw. Element ${\displaystyle K_{j}}$ gehörende Elementsteifigkeitsmatrizen durch

${\displaystyle A_{h}^{j}=(a_{ik}^{j})_{i,k\in I_{j}}}$

mit

${\displaystyle a_{ik}^{j}=\int _{K_{j}}\nabla \phi _{i}\nabla \phi _{k}\quad {\rm {und}}\,\,I_{j}=\{i:supp\,\phi _{i}\cap K_{j}\not =\emptyset \}.}$

Analog wird durch

${\displaystyle f^{j}=(f_{i}^{j})_{i\in I_{j}}\quad {\rm {mit}}\,\,f_{i}^{j}=\int _{K_{j}}f\,\phi _{i}.}$

eine elementweise rechte Seite erklärt. Durch Summation über alle ${\displaystyle K_{j}}$ 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 ${\displaystyle K}$ ein Element mit den Ecken ${\displaystyle p^{i},p^{k},p^{l}}$ und z. B. ${\displaystyle p^{i}=(x_{1},y_{1})}$. Dann kann man ${\displaystyle K}$ durch die Transformation

${\displaystyle x=x_{1}+(x_{2}-x_{1})\xi +(x_{3}-x_{1})\eta ,\quad y=y_{1}+(y_{2}-y_{1})\xi +(y_{3}-y_{1})\eta }$

auf das Referenzelement mit den Ecken ${\displaystyle (0,0),(1,0),(0,1)}$ in der ${\displaystyle (\xi ,\eta )}$-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 ${\displaystyle K}$ in der ${\displaystyle (\xi ,\eta )}$-Ebene die Basisfunktionen

${\displaystyle 1-\xi -\eta ,\,\,\xi \,\,{\rm {und}}\,\,\eta .}$

Die Berechnung deren Ableitung nach ${\displaystyle \xi }$ und ${\displaystyle \eta }$ ist trivial. Die Funktionaldeterminante ${\displaystyle D}$ der Transformation ist

${\displaystyle D=(x_{2}-x_{1})(y_{3}-y_{1})-(x_{3}-x_{1})(y_{2}-y_{1}),}$

der Betrag von ${\displaystyle D}$ ist gleich dem doppelten Flächeninhalt von ${\displaystyle K}$. Letztlich muss man noch die Ableitungen nach ${\displaystyle x}$ und ${\displaystyle y}$ in Ableitungen nach ${\displaystyle \xi }$ und ${\displaystyle \eta }$ umrechnen. Es gilt

${\displaystyle D\,{\frac {\partial }{\partial x}}=(y_{3}-y_{1}){\frac {\partial }{\partial \xi }}+(y_{1}-y_{2}){\frac {\partial }{\partial \eta }}\quad {\rm {und}}\,\,D\,{\frac {\partial }{\partial y}}=(x_{1}-x_{3}){\frac {\partial }{\partial \xi }}+(x_{1}-x_{2}){\frac {\partial }{\partial \eta }}.}$

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

${\displaystyle a_{11}={\frac {1}{2|D|}}[(y_{2}-y_{3})^{2}+(x_{3}-x_{2})^{2}],\,\,a_{12}={\frac {1}{2|D|}}[(y_{2}-y_{3})(y_{3}-y_{1})+(x_{3}-x_{2})(x_{1}-x_{3})]\,\,a_{13}={\frac {1}{2|D|}}[(y_{2}-y_{3})(y_{1}-y_{2})+(x_{3}-x_{2})(x_{2}-x_{1})].}$

Zurückkehrend zu der Bezeichnung ${\displaystyle K_{j}}$ mit den Ecken ${\displaystyle p^{i}=(x_{i},y_{i})}$, ${\displaystyle p^{k}=(x_{k},y_{k})}$ und ${\displaystyle p^{l}=(x_{l},y_{l})}$ ergibt sich für die Elementsteifigkeitsmatrix mit ${\displaystyle i\not =k}$

${\displaystyle a_{ik}^{j}={\frac {1}{2|D_{j}|}}[(y_{l}-y_{k})(y_{i}-y_{l})+(x_{l}-x_{k})(x_{i}-x_{l})]}$

und für die Diagonalelemente

${\displaystyle a_{ii}^{j}={\frac {1}{2|D_{j}|}}[(y_{k}-y_{l})^{2}+(x_{k}-x_{l})^{2}].}$

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.