Unterabschnitte

• Koordinatentransformation
• Ansatzfunktionen
• Elementintegral für den Laplace-Term
• Elementintegral für den Poisson-Term
• Elementintegral am Rand

4.4 Elementformulierung

Für die Berechnung der Näherungslösung müssen Volumsintegrale über das gesamte Simulationsgebiet $\Omega$ und über den Rand $\Gamma$ gelöst werden. Die Volumsintegrale über $\Omega$ können als Summe der Anteile der einzelnen Elemente angeschrieben werden.

$$I=\sum_{i=1}^M I^{e_i}$$ (4.41)

Als Gitterelemente sollen hier ausschließlich Tetraeder verwendet werden, außer natürlich für die Beschreibung des Randes des Simulationsgebietes (zur Vorgabe von Randbedingungen) -- hier werden Dreiecke verwendet, die konform zu den an den Rand grenzenden Tetraederflächen sind. Für zweidimensionale Simulationen werden stattdessen Dreiecke und Linienelemente verwendet. Die ausführliche Behandlung dieser Elemente soll jedoch im Rahmen dieser Arbeit nicht erfolgen, da sie völlig analog zum dreidimensionalen Fall durchgeführt werden kann.

Prinzipiell ist man bei der Finiten Elemente Methode nicht auf Tetraeder bzw. Dreiecke beschränkt. Ebenso können beispielsweise Hexaeder, Prismen oder Oktaeder verwendet werden. Auch Gitter mit gemischten Elementtypen sind möglich.

Koordinatentransformation

Da die Integration über ein Element für Tetraeder und Dreiecke meist analytisch durchführbar ist, soll sie im lokalen Elementkoordinatensystem $(\xi ,\eta ,\zeta )$ erfolgen. Der Integrand muss dafür vom $(x,y,z)$- ins $(\xi ,\eta ,\zeta )$-Koordinatensystem transformiert werden (siehe Abbildungen 4.1 und 4.2).

Abbildung 4.1: Allgemeines Tetraederelement im $(x,y,z)$-Koordinatensystem

Abbildung 4.2: Tetraederelement im normierten $(\xi ,\eta ,\zeta )$-Koordinatensystem

Mit folgender Abbildungsvorschrift kann von lokalen Elementkoordinaten ins globale Koordinatensystem umgerechnet werden

$$\begin{array}{r@{}r@{}r@{}r@{}r} x=& x_0 +& (x_1-x_0)\xi +& (... ... +& (z_1-z_0)\xi +& (z_2-z_0)\eta +& (z_3-z_0)\zeta \end{array}$$ (4.42)

bzw. in Matrixschreibweise

$$\{r\}=\{r_0\} + [J]\cdot\{\delta\}$$ (4.43)

oder für die Umrechnung von globalen in lokale Koordinaten

$$\{\delta\}=[J]^{-1}\big(\{r\}-\{r_0\}\big)$$ (4.44)

mit der sogenannten Jakobi-Matrix

$$[J]=\left(\begin{array}{ccc} x_1-x_0 & x_2-x_0 & x_3-x_0 \\ y_1-y_0 & y_2-y_0 & y_3-y_0 \\ z_1-z_0 & z_2-z_0 & z_3-z_0 \end{array}\right)$$ (4.45)

sowie $\{r\}=(x,y,z)^T$ und $\{\delta\}=(\xi,\eta,\zeta)^T$.

Ansatzfunktionen

Die Näherungslösung $\tilde u$ wird gemäß (4.6) als Linearkombination der Ansatzfunktionen $N_j$ dargestellt. Für eine Behandlung mit der Methode der gewichteten Residuen müssen die Ansatzfunktionen am Übergang von einem Element zum Nachbarelement zumindest stetig sein und innerhalb eines Elements zumindest einmal stetig differenzierbar sein. Die Ansatzfunktionen $N^e_k(\xi,\eta,\zeta)$ werden nun lokal auf dem Einheitstetraeder definiert. Dabei ist anzumerken, dass sich die globale Ansatzfunktion $N_j(x,y,z)$ aus den lokalen Ansatzfunktionen der am Knoten $j$ angrenzenden Tetraederelemente zusammensetzt. Die einfachste Art von Ansatzfunktionen sind lineare Funktionen, die für den Einheitstetraeder wie folgt gegeben sind:

$$\begin{split}N^e_{L0}&=1-\xi-\eta-\zeta \\ N^e_{L1}&=\xi\\ N^e_{L2}&=\eta\\ N^e_{L3}&=\zeta \end{split}$$ (4.46)

Um eine beliebig quadratische Funktion im dreidimensionalen Raum anzugeben benötigt man zehn Parameter, daher benötigen Tetraederelemente mit quadratischem Ansatz auch zehn Stützpunkte. Man wählt deshalb außer den vier Eckpunkten des Tetraeders zusätzlich noch sechs Punkte, die jeweils in der Mitte der Kanten liegen.

Abbildung 4.3: Knotennummerierung bei einem Tetraeder mit quadratischen Ansatzfunktionen

Die quadratischen Ansatzfunktionen im Einheitstetraeder lassen sich durch die linearen ausdrücken:

$$\begin{split}N^e_{Q0}&=N^e_{L0}(2N^e_{L0}-1)\\ N^e_{Q1}&=N^e... ...Q8}&=N^e_{L1}N^e_{L3}\\ N^e_{Q9}&=N^e_{L2}N^e_{L3} \end{split}$$ (4.47)

In Abb. 4.4 sind die Ansatzfunktionen auf Einheitstetraedern in Form von Isoflächen dargestellt. Dabei entspricht die Farbe Blau einem Wert von Null und Rot einem Wert von eins. Hier fällt auf, dass die Funktionen $N^e_{Q0}$ bis $N^e_{Q3}$ auch negativ werden (Bereich zwischen den beiden dunkelblauen Isoflächen).

Abbildung 4.4: Isoflächendarstellung der quadratischen Ansatzfunktionen $N^e_{Q0}\ldots N^e_{Q9}$ (in dieser Reihenfolge) auf einem Tetraederelement. Jede der Ansatzfunktionen hat auf genau einem Knoten einen Wert von 1 und 0 auf allen anderen.

Zu jedem Elementknoten $k$ gibt es genau eine Ansatzfunktion $N^e_k$, die an diesem Knoten gleich 1 ist und 0 auf allen anderen. Für die Summe aller Ansatzfunktionen gilt

$$\sum_{k=0}^{nl-1}N^e_k(\xi,\eta,\zeta)=1$$ (4.48)

Im Folgenden sollen die Integrale über $\Omega$ und $\Gamma$ auf einem einzelnen Gitterelement gelöst werden, indem das Integral in das lokale Elementkoordinatensystem transformiert wird und dort nach einer analytischen Lösung gesucht wird.

Elementintegral für den Laplace-Term

Für ein einzelnes Element $E$ sieht das Integral für den Term

$$I^e_L=\int_E(\nabla N_i)\makebox{\boldmath$\underline m$}(\nabla N_j)\,\textrm{d}x\,\textrm{d}y\,\textrm{d}z$$ (4.49)

folgendermaßen aus:

$$\begin{split}I_L^e &=\int_E\Big[ m_{11}\frac{\partial{N^e_k}}... ...\Big] \,\textrm{d}x \,\textrm{d}y \,\textrm{d}z \\ \end{split}$$ (4.50)

Das Integral kann nun auf das lokale $(\xi ,\eta ,\zeta )$-Koordinatensystem im Element transformiert werden.

$$\begin{split}I_L^e =\int_E\Bigg[ m_{11}&\!\left(\frac{\partia... ...[J]\,\textrm{d}\xi\,\textrm{d}\eta\,\textrm{d}\zeta \end{split}$$ (4.51)

Man kann nun den Integrand nach den darin auftretenden Ableitungen im lokalen Koordinatensystem zusammenfassen und erhält

$$\begin{split}I_L^e =\int_E\Big(& g_0 \frac{\partial{N^e_k}}{\... ...ig)\,\textrm{d}\xi\,\textrm{d}\eta\,\textrm{d}\zeta \end{split}$$ (4.52)

Zur Vereinfachung definiert man die Matrix $[K]$ als $[K]=\det[J]\cdot[J]^{-1}$, setzt ihre Elemente in (4.51) ein und erhält somit für die Faktoren $g_0$ bis $g_5$

$$\begin{split}g_0 & = \left( m_{11} K_{11}^2 + m_{22} K_{21}^2... ...3} ( K_{33} K_{11} + K_{31} K_{13} )\right)/\det[J] \end{split}$$ (4.53)

Für elementweise konstantes $\makebox{\boldmath $\underline m$}$ sind auch die sechs Koeffizienten $g_0 \ldots g_5$ auf dem gesamten Gitterelement konstant, sie sind lediglich von der Lage des Elements im globalen Koordinatensystem $(x,y,z)$ und den Komponenten des Tensors  $\makebox{\boldmath $\underline m$}$ abhängig.

Für isotrope Medien gilt $m_{11}= m_{22}= m_{33} = m$ und $m_{12}= m_{23}= m_{13}=0$. Die Koeffizienten lassen sich dadurch vereinfachen:

$$\begin{split}g_0 &= m(K_{11}^2 + K_{21}^2 + K_{31}^2)/\det[J]... ...3}K_{11} + K_{23}K_{21} + K_{33}K_{31})/\det[J] \\ \end{split}$$ (4.54)

Das Integral (4.52) kann nun in sechs Integrale über die einzelnen Terme aufgespalten und die Faktoren $g_0 \ldots g_5$ vor diese Integrale gehoben werden. Die einzelnen Integranden bestehen jetzt nur mehr aus Produkten der Ableitungen der Ansatzfunktionen im lokalen Elementkoordinatensystem und können analytisch berechnet werden.^4.4Diese Integrationen brauchen für alle Gitterelemente nur einmal ausgeführt werden. Für die praktische Anwendung berechnet man sie für alle Paarungen $(k,l)$ von Ansatzfunktionen eines Elements und speichert sie in Matrizen (den sogenannten Elementmatrizen) die wie folgt gebildet werden:

$$\begin{split}[A_{L0}^e]&=\int_E \left\{\frac{\partial{N^e}}{\... ...,\textrm{d}\xi\,\textrm{d}\eta\,\textrm{d}\zeta \\ \end{split}$$ (4.55)

Die sechs Elementmatrizen $A^e_{L0}\ldots A^e_{L5}$ können nun zu einer gesamten Elementmatrix $A^e_L$ aufsummiert werden. Da jede der Elementmatrizen symmetrisch ist, hat natürlich auch die gesamte Elementmatrix Symmetrieeigenschaft.

$$[A^e_L]= g_0[A_{L0}^e] + g_1[A_{L1}^e] + g_2[A_{L2}^e] + g_3[A_{L3}^e] + g_4[A_{L4}^e] + g_5[A_{L5}^e]$$ (4.56)

Dabei sind die Einträge der gesamten Elementmatrix $[A^e_L]$ die Lösungen von (4.49) für alle Kombinationen von lokalen Ansatzfunktionen.

Elementintegral für den Poisson-Term

Das Integral

$$I^e_P=\int_E N_i N_j\,\textrm{d}x \,\textrm{d}y\,\textrm{d}z$$ (4.57)

wird zuerst ins Elementkoordinatensystem $(\xi ,\eta ,\zeta )$ transformiert

$$I_P^e=\int_E N^e_k N^e_l\det[J]\,\textrm{d}\xi\,\textrm{d}\eta\,\textrm{d}\zeta$$ (4.58)

und kann dann analytisch gelöst werden, indem die Jakobi-Determinante vor das Integral gehoben wird.^4.5Man kann nun für alle Paarungen der lokalen Ansatzfunktionen eine entsprechende Elementmatrix aufstellen

$$[A_{P0}^e]=\int_E\{N^e\}\cdot\{N^e\}^T\,\textrm{d}\xi\,\textrm{d}\eta\,\textrm{d}\zeta \,,$$ (4.59)

die ausschließlich von den Ansatzfunktionen abhängig ist und für alle Elemente nur ein einziges Mal berechnet werden muss. Die Matrix

$$[A_P^e]=\det[J]\cdot [A_{P0}^e]$$ (4.60)

enthält dann die Lösungen von (4.57) für alle Kombinationen von lokalen Ansatzfunktionen.

Elementintegral am Rand

Das Integral über die Oberfläche des Gebiets $\Omega$ kann als Summe der Anteile der Randelemente aufgefasst werden.

$$I_B^e=\int_{E_B}\!N_i N_j \,\textrm{d}A$$ (4.61)

Bei einem dreidimensionalen Simulationsbereich mit Tetraedern als Gitterelemente sind die Randelemente Dreiecke. Bei der Koordinatentransformation wird nun ein Punkt auf dem Dreieck mit den globalen Koordinaten $(x,y,z)$ in die normierten Elementkoordinaten $(\xi,\eta)$ transformiert. Dabei ergibt sich folgende Jakobi-Determinante:

$$\det[J]=\vert(\mathchoice{\mbox{\boldmath$\displaystyle r_1$}} {\... ...{\boldmath$\scriptstyle r_0$}} {\mbox{\boldmath$\scriptscriptstyle r_0$}})\vert$$ (4.62)

wobei mit $\mathchoice{\mbox{\boldmath $\displaystyle r_0$}} {\mbox{\boldmath $\textstyle... ...mbox{\boldmath $\scriptstyle r_2$}} {\mbox{\boldmath $\scriptscriptstyle r_2$}}$ die Ortsvektoren der Eckpunkte des Dreiecks bezeichnet werden. Man kann das Integral ebenfalls analytisch lösen^4.6

$$[A_B^e]=\det[J]\int_{E_B}\{N^{e_B}\}\cdot\{N^{e_B}\}^T \,\textrm{d}\xi\,\textrm{d}\eta$$ (4.63)

wobei $\{N^{e_B}\}$ die entsprechenden Ansatzfunktionen der Randelemente sind, die gemäß (4.46) bzw. (4.47) gegeben sind, wobei für die $\zeta$-Koordinate immer 0 einsetzt wird.

Fußnoten

... werden.^4.4

Eine analytische Lösung in dieser Form ist nur dann möglich wenn die Faktoren $g_0 \ldots g_5$ über das gesamte Element konstant sind. Das ist genau dann der Fall, wenn das Element eine konstante Jakobi-Determinante besitzt (ist bei den hier verwendeten Tetraederelementen immer erfüllt) und der Materialfaktor  $\makebox{\boldmath$\underline m$}$ im Element ebenfalls konstant ist. Anderenfalls muss für die Integration eine numerische Approximation verwendet werden.

... wird.^4.5

Voraussetzung dafür ist wieder, dass die Jakobi-Determinante über das gesamte Element konstant ist.

... lösen^4.6

Eine konstante Jakobi-Determinante ist wieder Voraussetzung, ist aber bei Dreieckelementen immer gegeben.