v4.6 – Implementierung

 

4.6.1 Allgemeine Vorgehensweise

Zunächst muss in geeigneter Form das Problem beschrieben werden, nämlich die Differentialgleichung, die Randbedingungen und das Gebiet. Ebenfalls zum Preprocessing gehört die Vernetzung des Gebiets. Je allgemeiner die zugelassenen Gebiete sind, desto höher sind die Anforderungen an den Vernetzungsalgorithmus.

Auch wenn man bei der Aufstellung des Gleichungssystems \left( {14} \right) an knotenweise definierte Funktionen denkt, wird die Steifigkeitsmatrix und die rechte Seite elementweise erzeugt und assembliert. Die Integration wird dabei näherungsweise mit Quadraturformeln ausgeführt. Der dabei entstehende Fehler kann mit Hilfe des 2. Lemmas von Strang abgeschätzt werden, wird jedoch in den theoretischen Betrachtungen meist vernachlässigt.

Ein Vorteil der Finite-Elemente-Methode besteht darin, dass die Steifigkeitsmatrix nur dünn besetzt ist (sparse). Für deren Speicherung stehen verschiedene Techniken zu Verfügung. Verbreitet ist die Kompaktliste zeilenweise (compressed row storage), wobei nur die Nichtnullelemente der Matrix gespeichert werden.

Das entstehende Gleichungssystem kann mit einem direkten Verfahren gelöst werden, das die dünne Besetztheit geschickt ausnutzt. Bei großer Dimension des Gleichungssystems sind jedoch iterative Verfahren wie das Verfahren der konjugierten Gradienten, effektiver, vor allem wenn sie mit einer guten Vorkonditionierung kombiniert werden. Die besten Vorkonditionierungen basieren auf einer Netzhierarchie, die die zu implementierende Software entsprechend komplex werden lässt.

Schließlich gehört zu jeder Implementierung ein Postprocessing, das die graphische Darstellung der Ergebnisse, die Berechnung mittelbarer Größen (Ableitungen oder Funktionalwerte der Lösung) und eine a posteriori Fehlerbewertung des Ergebnisses beinhaltet. An dieser Stelle entsteht im Allgemeinen eine Schleife, wobei man aus den a posteriori Fehlerinformationen einen neuen, besser geeigneten Finite-Elemente-Raum konstruiert und wieder in die Phase der Assemblierung und Lösung des Gleichungssystems eintritt.

Auf Grund der Spezifik unterschiedlicher Typen partieller Differentialgleichungen und der Vielfalt möglicher Diskretisierungen und anderer Einzelbausteine (Speichertechniken, Löser, etc.) gibt es keine universellen Softwarelösungen.

4.6.2 Datenstrukturen

Ein Finite-Elemente-Netz besteht im 2D-Fall aus Knoten V (vertices), Kanten E (edges) und Elementen T. Im 3D-Fall kommen noch Seiten F (faces) hinzu.

In jedem Fall gibt es ein Feld V2X\left( {d,{N_V}} \right), in dem zu jedem Knoten die Koordinaten gespeichert sind.

Im einfachen Fall (reduzierte Datenstruktur) gibt es nur noch ein Feld T2V, in dem zu jedem Element die Knotennummern gespeichert sind.

Im allgemeinen Fall (vollständige Datenstruktur) gibt es folgende Felder

2D: T2E,\:\:E2V
3D: T2F,\:\:F2E,\:\:E2V

Weitere Felder können nützlich sein:

  • Liste der Randseiten
  • Liste der Nachbarelemente
  • Vater-Sohn-Beziehungen in hierarchischen Netzen
  • Materialbeschreibung

Beispiel: Dreiecksnetz für \Omega = {\left( {0,1} \right)^2}

num-406-dreiecksnetz-datenstruktur-beispiel

Listen:

num-406-knotenliste-kantenliste-relation-beispiel

T2V entspricht einer Zuordnung von lokalen Knotennummern zu globalen Knotennummern. Diese ist nicht eindeutig.

4.6.3 Aufbau des Finite-Elemente-Gleichungssystems

Zu lösen ist das lineare Gleichungssystem {K_h}{u_h} = {f_h} mit

{K_h} = \left[ {a\left( {{\varphi _j},{\varphi _i}} \right)} \right]_{i,j = 1}^n,\quad \quad {f_h} = \left[ {\left( {f,{\varphi _i}} \right)} \right]_{i = 1}^n

a\left( {{\varphi _j},{\varphi _i}} \right) = \int\limits_\Omega {\nabla {\varphi _j} \cdot \nabla {\varphi _i}} ,\quad \quad \left( {f,{\varphi _i}} \right) = \int\limits_\Omega {f{\varphi _i}}

Die Berechnung der Bilinearform wurde bereits am Beispiel erläutert. Obwohl hier knotenbezogene Größen eingehen ({\varphi _i}), läuft die Schleife über die Elemente:

\int\limits_\Omega {\nabla {\varphi _j} \cdot \nabla {\varphi _i}} = \sum\limits_{T \in {\mathcal{T}_h}} {\int\limits_T {\nabla {\varphi _j} \cdot \nabla {\varphi _i}} }

Für ein Element T \in {\mathcal{T}_h} ist das Integral nur ungleich 0, wenn {x_i} und {x_j} Knoten von T sind.

Wir betrachten das lineare Dreieckselement {T_\mu } mit den Ecken , und . Der Beitrag von zur Steifigkeitsmatrix hat die Form

{\tilde K_{h,\mu }} = \left( {\begin{array}{*{20}{c}}{} & {} & {} & {} & {} & {} \\{} & {K_{i,i}^\mu } & {K_{i,j}^\mu } & {} & {K_{i,l}^\mu } & {} \\{} & {K_{j,i}^\mu } & {K_{j,j}^\mu } & {} & {K_{j,l}^\mu } & {} \\{} & {} & {} & {} & {} & {} \\{} & {K_{l,i}^\mu } & {K_{l,j}^\mu } & {} & {K_{l,l}^\mu } & {} \\{} & {} & {} & {} & {} & {} \\ \end{array} } \right)

(überall anders stehen Nullen!) Die auf Nicht-Nullelemente reduzierte Matrix

{K_{h,\mu }} = \left( {\begin{array}{*{20}{c}}{K_{i,i}^\mu } & {K_{i,j}^\mu } & {K_{i,l}^\mu } \\{K_{j,i}^\mu } & {K_{j,j}^\mu } & {K_{j,l}^\mu } \\{K_{l,i}^\mu } & {K_{l,j}^\mu } & {K_{l,l}^\mu } \\ \end{array} } \right)\quad \in {\mathbb{R}^{3 \times 3}}

heißt Elementsteifigkeitsmatrix. Die Beziehung zwischen {\tilde K_{h,\mu }} und {K_{h,\mu }} wird durch T2V definiert. Die Addition der Elemente von {K_{h,\mu }} an die entsprechenden Stellen von {K_h} heißt Assemblierung.

Auch für den Lastvektor (rechte Seite) führen wir die Berechnung am Referenzelement durch und transformieren dann auf das echte Element:

\int\limits_T {f\left( x \right){\varphi _i}\left( x \right)dx} = \int\limits_{\hat T} {f\left( {F\left( {\hat x} \right)} \right){{\hat \varphi }_i}\left( {\hat x} \right)\left| {\frac{{\partial \left( {{x_1},{x_2}} \right)}}{{\partial \left( {{{\hat x}_1},{{\hat x}_2}} \right)}}} \right|d\hat x}

\approx \sum\limits_{j = 1}^{{n_{\operatorname{int} }}} {{w_j}f\left( {F\left( {{{\hat x}_j}} \right)} \right){\varphi _i}\left( {{{\hat x}_j}} \right)J\left( {{{\hat x}_j}} \right)}

Die Daten werden also nur an Punkten ausgewertet. Dies wird bei stark oszillierenden Funktionen kritisch (Shannon-Theorem!).

4.6.4 Berechnung der Element-Steifigkeitsmatrix

Die Integrale \int\limits_T {\nabla {\varphi _j} \cdot \nabla {\varphi _i}} berechnet man mit Hilfe von Quadraturformeln, die Quadraturpunkte über das Referenzsystem:

num-406-quadraturpunkte-referenzelement-allgemein

Einbau von homogenen Dirichlet-Randbedingungen: Man assembliert die Matrix und die rechte Seite und streicht dann Zeilen und Spalten, die zu Randknoten gehören.

Eigenschaften der Systemmatrix:

Die Matrix {K_h}

  • hat die Dimension n \approx {h^{-d}}
  • ist dünn besetzt, die Anzahl der Nicht-Nullelemente pro Zeile ist O\left( 1 \right), die der gesamten Matrix O\left( n \right)
  • ist schlecht konditioniert: \kappa \left( {{K_h}} \right) = O\left( {{h^{-2}}} \right) unabhängig von der Raumdimension

4.6.5 Methoden zur Netzgenerierung

4.6.5.1 Manuelle Methode

Eingabe aller Daten per Hand. Das ist praktikabel für Gebiete mit einfacher Geometrie und sehr große Netze.

4.6.5.2 Produktmethode

Geeignet für prismatische oder zylinderförmige Gebiete:

num-406-produktmethode-netzgenerierung-finite-elemente

4.6.5.3 Transformation von Referenznetzen

Nützlich bei topologisch einfachen Gebieten mit gekrümmten Rändern.

num-406-transformation-netzgenerierung-finite-elemente

4.6.5.4 Halbautomatische Methode

Vorgehen:

  • Zerlegung des Gebiets in Teilgebiete
  • Teilgebiete werden mit einer der vorherigen Methoden vernetzt

num-406-halbautomatisch-netzgenerierung-finite-elemente

Achtung: Sichern einer zulässigen Vernetzung!

4.6.5.5 Quadtree- und Octree-Methode

Vorgehen:

  • Start mit 1 Rechteck, das \Omega enthält
  • rekursive Teilung der Elemente, bis der Rand von \Omega hinreichend gut approximiert wird
  • Streichen aller Elemente, die außerhalb von \Omega liegen
  • Verschieben von Knoten und Kanten zur besseren Randapproximation

num-406-quadtree-octree-netzverfeinerung-kante

Man kann hier sehr schnell bestimmen, in welchem Element ein gegebener Punkt liegt.

4.6.5.6 Delauney-Vernetzung

Gegeben: N Punkte {x_i} \in {\mathbb{R}^d}

Vorgehen:

  • Konstruiere Voronoi-Polyeder V\left( {{x_i}} \right) = \left\{ {x \in {\mathbb{R}^d}:\left\| {x-{x_i}} \right\| < \left\| {x-{x_j}} \right\|\:\:\forall j \ne i} \right\} durch Mittelsenkrechten

    num-406-netzgenerierung-voronoi-polyeder

  • Verbinde diejenigen Punkte, deren Voronoi-Polyeder eine gemeinsame Kante oder Fläche besitzen

    num-406-delauney-triangulierung-netzgenerierung

Die Delauney-Triangulierung trianguliert die konvexe Hülle der Punkte {x_i}.

Im Inneren jeden Umkreises zu einem Dreieck liegt kein weiterer Punkt. D.h. die Summe der einer Seite gegenüberliegenden Innenwinkel ist kleiner als oder gleich 180°.

4.6.5.7 Advancing-Front-Methoden

Vorgehen:

  • Vernetze den Gebietsrand
  • Generiere ein Element, das mindestens eine Seite auf der Front besitzt, nach verschiedenen Kriterien, z.B.
    • Lücken schließen
    • spitze Winkel in Front eliminieren
  • Definiere neue Front (Rand des noch nicht vernetzten Gebietes). ist diese leer, gehe zum zweiten Schritt.

4.6.6 Netzverfeinerung

Man kann geeignete Finite-Elemente-Netze auch durch sukzessive Verfeinerung von groben Ausgangsnetzen gewinnen. Die lokale Netzdichte kann gesteuert werden durch

  • a-priori Informationen
  • a-posteriori Informationen (lokale Fehlerindikatoren, die basierend auf der Lösung auf den bisherigen Netz berechnet werden)

Teilung von Elementen:

2D: Viertelung von Dreiecken (rote Verfeinerung) erzeugt kongruente Dreiecke

num-406-netzverfeinerung-rote-viertelung

3D: Achtelung von Tetraedern liefert im Allgemeinen keine kongruenten Tetraeder. Es existieren verschiedene Strategien.

Falls nicht alle Dreiecke geviertelt werden, entsteht im Allgemeinen keine zulässige Zerlegung (hängende Knoten entstehen). Abhilfe: angrenzende Dreiecke werden halbiert (grüne Verfeinerung)

num-406-grune-verfeinerung-netzgenerierung

Um nicht immer kleinere Winkel zu erhalten, werden grüne Dreiecke nicht nochmals halbiert.

num-406-schlechte-richtige-netzverfeinerung-element

Ähnliche Artikel

Kommentar verfassen