Finite-Elemente-Verfahren - Wissenschaftliches Rechnen (i788)

Finite-Elemente-Verfahren (kurz: FEM für finite elemente methods) basieren auf der schwachen Formulierung von PDE. Hauptmerkmal ist, dass für die Diskretisierung Basisfunktionen $v_k:B\to\bbR$ mit sehr kleinem Träger

\supp v_k:=\{x\in B:\;v_k(x)\neq 0\}

(1)

verwendet werden. Dadurch entsteht ein lineares Gleichungssystem mit dünn besetzter Systemmatrix, welches effizient gelöst werden kann. Für die Diskretisierung werden keine regelmäßigen Gitter von Stützstellen benötigt, sodass auch PDE auf komplexen Gebieten $B$ mittels FEM gelöst werden können.

11.1Ausgangspunkt¶

Gegeben sei die schwache Formulierung einer PDE in der Form

a(u,v)=b(v)\quad\text{für alle }v\in V,

(2)

wobei $a$ , $b$ , $V$ gegeben sind und $u$ gesucht ist (vgl. Abstrakte Form schwacher Formulierungen).

Weiter sei $V_n\subseteq V$ ein endlichdimensionaler Teilraum von $V$ und $\{v_1,\ldots,v_n\}$ sei eine Basis in $V_n$ , d.h. jede Funktion in $V_n$ kann als Linearkombination der Funktionen $v_1,\ldots,v_n$ dargestellt werden. Durch Einschränkung der schwachen Formulierung auf $V_n$ erhalten wir das diskretisierte Problem

A\,\underline{u}=\underline{b}

(3)

mit der rechten Seite

\underline{b}:={\begin{bmatrix}b(v_1)\\\vdots\\b(v_n)\end{bmatrix}}

(4)

und der Systemmatrix (hier auch Steifigkeitsmatrix genannt)

A:={\begin{bmatrix} a(v_1,v_1)&\cdots&a(v_n,v_1)\\ \vdots&&\vdots\\ a(v_1,v_n)&\cdots&a(v_n,v_n)\\ \end{bmatrix}}

(5)

(vgl. Diskretisierung der abstrakten Formulierung). Aus dem Vektor $\underline{u}\in\bbR^n$ erhalten wir mittels

u:=\underline{u}_1\,v_1+\cdots+\underline{u}_n\,v_n

(6)

die gesuchte Lösung der (schwachen Formulierung der) PDE.

Einziger noch zu klärender Punkt ist die konkrete Wahl der Basisfunktionen $v_k$ .

11.2Wahl der Basis¶

Die Wahl der Basis $\{v_1,\ldots,v_n\}$ erfolgt in zwei Schritten:

Zerlegung des Gebiets $B$ in disjunkte gleichartige Teilmengen $T_1,\ldots,T_N\subseteq B$ (die sogenannten finite Elemente). Dabei ist mit “gleichartig” gemeint, dass zum Beispiel alle Teilmengen Dreiecke sind oder alle sind Vierecke oder alle sind Tetraeder (für $B\subseteq\bbR^3$ ) usw.
Wahl der $v_k$ so, dass die Träger $\supp v_k$ jeweils der Vereinigung von nur wenigen Teilgebieten $T_l$ entsprechen.

Die Zerlegung in Teilgebiete erfolgt stets so, dass sich nur Ecken mit Ecken berühren, nie Ecken mit Kanten oder Flächen. Die Ecken werden auch Knoten genannt und im Folgenden mit $P_1,\ldots,P_m$ bezeichnet.

Sind die $T_l$ Dreiecke, so werden die $v_k$ gern als stückweise lineare Funktionen gewählt: Zu jedem Knoten $P_k$ wird genau eine Basisfunktion gewählt. Diese ist 1 an dem Knoten und Null an allen anderen Knoten (IDVID 1120).

Analog kann bei Tetraedern in $\bbR^3$ vorgegangen werden. Neben stückweise linearen sind auch stückweise quadratische usw. Basisfunktionen üblich. Da Funktionen höherer Ordnung mehr Freiheitsgrade besitzen, werden je nach Bedarf zusätzliche Knoten auf den Berührungskanten oder -flächen der Teilgebiete $T_l$ eingefügt.

11.3Gebietszerlegung¶

Die konkrete Wahl der Zerlegung in Teilgebiete unterliegt diversen Einflussfaktoren:

erforderliche Genauigkeit der Berechnungen (ggf. höhere Knotendichte an “interessanten” Stellen),
exakte Darstellung von Materialgrenzen (diese sollten immer auf Kanten/Flächen der Zerlegung liegen),
exakte Darstellung der Grenzen zwischen verschiedenen Typen von Randbedingungen,
exakte Darstellung von Unstetigkeitsstellen in den Randbedingungen,
aus theoretischen Ergebnissen bekannte Einflüsse auf den Diskretisierungsfehler (z.B. dichtere Knoten in der Nähe konkaver Gebietsgrenzen),
höhere Knotendichte an krumlinigen bzw. krumflächigen Gebietsgrenzen (hinreichend genaue Approximation durch Polygon bzw. Polyeder).

Die Vernetzung (d.h. die Wahl der Zerlegung) kann manuell oder durch Algorithmen erfolgen. Ein verbreiteter Algorithmus ist die Delaunay-Triangulation.

Gebietszerlegungen können besser oder schlechter für FEM geeignet sein als andere. Insbesondere sollte die Form der Teilmengen hinreichend einheitlich sein. Diese Forderung kann wie folgt präzisiert werden:

Neben Algorithmen zur Erzeugung von Vernetzungen existieren auch Algorithmen zur Verbesserung von Netzen im Sinne der Isotropie (Netzglättung).

Auch können gegebene Netze weiter verfeinert werden durch weiteres Unterteilen der vorhandenen Teilgebiete (Netzverfeinerung). Insbesondere die adaptive Netzverfeinerung ist praktisch von Bedeutung: In Regionen, in denen eine vorherige Netzverfeinerung zu einer deutlichen Änderung der Lösung geführt hat, wird weiter verfeinert. In Regionen ohne nennenswerte Änderung ist keine weitere Verfeinerung nötig. Durch dieses Vorgehen können Rechenzeit und Speicherplatz gespart werden bei gleichzeitig kleinem Diskretisierungsfehler.

11.4Aufstellen der Steifigkeitsmatrix¶

Ausgehend von einer gegebenen Zerlegung des Gebiets $B$ in Teilgebiete $T_1,\ldots,T_N$ soll die Steifigkeitsmatrix $A$ berechnet werden. Wir demonstrieren das Vorgehen anhand einer Dreiecksvernetzung im $\bbR^2$ und stückweise linearen Basisfunktionen. Für andere Vernetzungen und Basisfunktionen ist das Vorgehen völlig analog.

11.4.1Globale und lokale Knotennummerierung¶

Die Eckpunkte der Dreiecke seien $P_1,\ldots,P_m$ (Knoten). Punkte, an denen zwei oder mehr Dreiecke sich berühren, werden nur einfach gezählt. Die Nummerierung der Knoten mit $1,\ldots,m$ wird als globale Nummerierung bezeichnet.

Zu jedem Dreieck $T_l$ werden die drei Eckpunkte zusätzlich noch mit $P^l_1,P^l_2,P^l_3$ bezeichnet. Dies ist die lokale Knotennummerierung.

Im Speicher werden einerseits die Koordinaten der Knoten $P_1,\ldots,P_m$ zusammen mit der jeweiligen globalen Knotenummer abgelegt; andererseits die Elementzusammenhangstabelle, welche für jedes Teilgebiet $T_l$ (Zeilen) folgende Werte (Spalten) enthält:

Teilgebietnummer $l$ ,
globale Nummer des Knotens $P^l_1$ ,
globale Nummer des Knotens $P^l_2$ ,
globale Nummer des Knotens $P^l_3$ ,
Parameter in der PDE für dieses Teilgebiet (Material, Leitfähigkeit,...),

Durch Knotenkoordinaten und Elementzusammenhangstabelle ist die Vernetzung vollständig beschrieben (IDVID 1130).

11.4.2Referenzelemente¶

Für die Berechnung von $a(v_k,v_l)$ müssen Integrale über den Teilgebieten $T_l$ berechnet werden. Um dies effizient durchführen zu können, werden die auftretenden Integrale so transformiert, dass der Integrationsbereich immer gleich ausssieht (Referenzdreieck oder allgemeiner: Referenzelement).

Als Referenzdreieck wird üblicherweise das Dreieck mit den Eckpunkten $(0,0)$ , $(1,0)$ , $(0,1)$ verwendet. Liegt dieses in einem $\xi_1$ - $\xi_2$ -Koordinatensystem, ist die Transformation auf das Teilgebiet $T_l$ gegeben durch

{\begin{bmatrix}x_1\\x_2\end{bmatrix}}=P^l_1+\xi_1\,(P^l_2-P^l_1)+\xi_2\,(P^l_3-P^l_1).

(8)

11.4.3Formfunktionen über dem Referenzelement¶

Jede Basisfunktion $v_k$ soll an genau einem Knoten den Wert 1 annehmen, an allen anderen den Wert 0. Auf einem einzelnen Teilgebiet $T_l$ betrachtet sollen die Basisfunktionen linear verlaufen. Für dreieckige Teilgebiete gibt es somit nur vier mögliche Funktionsverläufe:

überall 0,
an einem der drei Eckpunkte 1, sonst 0.

Auf dem Referenzdreieck setzen wir:

$\varphi_1(\xi_1,\xi_2):=1-\xi_1-\xi_2\qquad$ (Wert 1 bei $(0,0)$ ),
$\varphi_2(\xi_1,\xi_2):=\xi_1\qquad$ (Wert 1 bei $(1,0)$ ),
$\varphi_3(\xi_1,\xi_2):=\xi_2\qquad$ (Wert 1 bei $(0,1)$ ).

11.4.4Basisfunktionen¶

Für jeden Knoten $P_k$ sei $I_k$ die Indexmenge, die gerade die am Knoten anliegenden Teilgebiete beschreibt. Die Basisfunktion $v_k$ ist dann gegeben durch

v_k(x)=\begin{cases} \varphi_{i(k,l)}(x),&\text{für }x\in T_l,\;l\in I_k,\\ 0,&\text{sonst}, \end{cases}

(9)

wobei $i(k,l)\in\{1,2,3\}$ gerade so gewählt wird, dass $P^l_{i(k,l)}=P_k$ gilt.

Sind die Funktionen in $V$ auf einem Teil des Randes Null (z.B. als Folge von Dirichlet-Randbedingungen in der PDE), so werden für die entsprechenden Randknoten keine Basisfunktionen benötigt!

11.4.5Bereichsintegrale¶

Zu einem Knoten $P_k$ einer Dreieckszerlegung enthalte $J_k$ alle Indizes von benachbarten Knoten. Ist die Bilinearform $a$ durch ein Integral über das Gebiet $B$ gegeben, so gilt $a(v_k,v_l)=0$ , falls $l\notin J_k$ . Jede Zeile und jede Spalte der Steifigkeitsmatrix hat also höchstens so viele Nicht-Null-Einträge wie die beiden relevanten Knoten Nachbarn haben (Zeile $l$ hat $|J_k|$ Einträge, Spalte $k$ hat $|J_l|$ Einträge. Insbesondere ist die Steifigkeitsmatrix dünn besetzt, was FEM für dreidimensionale Gebiete oder sehr feine Netze überhaupt erst möglich macht.

Die Berechnung der Bereichsintegrale über $B$ erfolgt praktisch immer numerisch, auch wenn manuelles Rechnen möglich wäre. Dazu werden die Integrale über den Teilgebieten $T_l$ zunächst auf das Referenzdreieck transformiert. Anschließend erfolgt die numerische Integration mittels Standardverfahren.

Sei

\Phi_l(\xi_1,\xi_2):=P^l_1+\xi_1\,(P^l_2-P^l_1)+\xi_2\,(P^l_3-P^l_1)

(10)

die Abbildung von $T_l$ auf das Referenzdreieck. Dann gilt

\int_{T_l}h\diff b=\int_0^1\int_0^{1-\xi_1} h\bigl(\Phi_l(\xi_1,\xi_2)\bigr)\,|\det J_{\Phi_l}(\xi_1,\xi_2)|\diff\xi_2\diff\xi_1

(11)

mit der zu integrierenden Funktion $h$ , wobei

\det J_{\Phi_k}(\xi_1,\xi_2)=[P^l_2-P^l_1]_1\,[P^l_3-P^l_1]_2-[P^l_2-P^l_1]_2\,[P^l_3-P^l_1]_1

(12)

(IDVID 1140).

11.4.6Randbedingungen¶

Eventuell vorhandene Dirichlet-Randbedingungen sind bereits in den Funktionenraum $V$ bzw. dessen endlichdimensionale Entsprechung $V_n$ integriert und wurden bei der Wahl der Basisfunktionen $v_1,\ldots,v_n$ berücksichtigt (keine Basisfunktionen an Randknoten im Bereich von Dirichlet-Randbedingungen).

Randteile, auf denen Neumann-Randbedingungen gegeben sind, führen üblicherweise zu Kurvenintegralen (2D-Gebiet) oder Oberflächenintegralen (3D-Gebiet) in der schwachen Formulierung der PDE. Diese können völlig analog zu den Bereichsintegralen numerisch berechnet werden. Führen dies hier nur für 2D-Gebiete näher aus.

Sei $\Gamma_2\subseteq\partial B$ der Teil des Randes, auf dem ein Kurvenintegral

\int_{\Gamma_2}h\diff s

(13)

mit einer Funktion $h:B\to\bbR$ zu berechnen ist. Das Integral kann in eine Summe von Integralen über die einzelnen Randkanten der Gebietszerlegung geschrieben werden. Ist $K$ eine solche Randkante mit den Endknoten $P_l$ und $P_k$ , so überführt man das Integral

\int_K h\diff s

(14)

in das Integral

\int_0^1 h\bigl(P_l+t\,(P_k-P_l)\bigr)\,\sqrt{[P_k-P_l]_1^2+[P_k-P_l]_2^2}\diff t

(15)

auf der Referenzkante $[0,1]$ (IDVID 1150).

Ist der Integrand $h$ als Linearkombination der Basisfunktionen $v_1,\ldots,v_n$ gegeben, so haben nur die beiden zu den Endknoten gehörenden Basisfunktionen $v_l$ und $v_k$ Einfluss auf das Integral.

11.5Aufstellen der rechten Seite¶

Das Aufstellen der rechten Seite $\underline{b}$ im zu lösenden Gleichungssystem erfolgt analog zur Systemmatrix durch numerische Berechnung der in der Linearform $b$ enthaltenen Integrale über jedem Teilgebiet $T_l$ .