Unterabschnitte

• Elektrostatisches Feld
  • Die elektrischen Ladungen
  • Schwebende Randbedingungen
  • Die Feldenergie
  
  
• Leitungsströme
  • Neumann-Bedingungen
  • Die Flächenstromdichte
  • Für schwebende Randbedingungen
  • Die Verlustleistung
  
  
• Elektro-quasistatische Feldberechnung
  • Die Zeitdiskretisierung
  
  
• Wärmeleitung stationär
• Wärmeleitung transient

4.5 Assemblierung

Nachdem nun die einzelnen Integralterme der Gleichungen des Galerkin Ansatzes elementweise bestimmt wurden, sollen nun die Beiträge der Elemente aufsummiert werden, um das Integral über das gesamte Gebiet $\Omega$ sowie den Rand $\Gamma$ zu bilden. Daraus kann ein Gleichungssystem für die unbekannten $u_i$ aufgestellt werden.

Elektrostatisches Feld

Das Gleichungssystem (4.28) kann in Matrixschreibweise angegeben werden

$$[A]\{\varphi\}=\{b\}\;.$$ (4.64)

Dazu berechnet man zunächst die einzelnen Elementmatrizen mittels (4.53) und (4.56), wobei für den Tensor $\makebox{\boldmath $\underline m$}$ die Permittivität $\makebox{\boldmath $\underline\varepsilon$}$ eingesetzt wird. Die Matrix $[A]$ wird dann durch Summierung über alle Elementmatrizen gebildet

$$[A]=\sum_{m=1}^M [A^{e_m}_L]\;.$$ (4.65)

Diese Gleichung ist allerdings nur symbolisch zu verstehen. Bei dieser Summierung wird nämlich von der lokalen Knotennummerierung auf die globale übergegangen. Hat beispielsweise der erste Knoten eines Elements den globalen Index 17 und der zweite den Index 62, so bedeutet dies, dass das Element (1,2) der Matrix $[A^{e_m}_L]$ an der Position (17,62) in der Matrix $[A]$ addiert wird. Dieser Vorgang wird als Assemblierung bezeichnet.

Gleichermaßen geht man bei der Summierung der Randterme von lokaler auf globale Knotennummerierung über

$$\{b\}=\sum_{m=1}^{M_B} [A^{e_m}_B] \{\sigma\}\;.$$ (4.66)

Die Spaltenmatrix $\{b\}$ ergibt sich als Summe der Beiträge der Elemente am Rand des Simulationsgebietes, kann also nur für Randknoten einen Wert ungleich Null haben. Genauer gesagt, tritt eine Flächenladung $\sigma$ nur an den Elektroden auf--das sind jene Knoten, die eine Dirichlet-Randbedingung haben. Deshalb können auch nur an diesen Stellen des $\{b\}$-Vektors Einträge ungleich Null vorhanden sein. Das Gleichungssystem (4.64) kann auch detailliert angeschrieben werden, wobei die obere Hälfte des Gleichungssystems den inneren Knoten und die untere Hälfte den Randknoten auf $\Gamma _1$ zugeordnet ist (die unbekannten Werte sind unterstrichen):

$$\left[\begin{array}{ccc\vert ccc} a_{11} &\cdots&a_{1,N_A} &a_{1,... ...line \underline b_{N_A+1} \\ \vdots \\ \underline b_N \\ \end{array}\right\}$$ (4.67)

oder als Teilmatrizen dargestellt

$$\begin{bmatrix}[A_A]\phantom{{}^T} & [A_B] \\ {} [A_B]^T & [A_D]... ...ay}\right\}=\left\{\begin{array}{c} \{0\} \\ \{ b_D \} \\ \end{array}\right\}$$ (4.68)

Um nun die unbekannten Potenzialwerte  $\{\varphi_A\}$ im Inneren des Leiters zu erhalten, muss man das folgende Gleichungssystem lösen:

$$[A_A]\cdot\{\varphi_A\}=-[A_B]\cdot\{\varphi_D\}$$ (4.69)

Die elektrischen Ladungen

auf den Elektroden werden berechnet, indem man die Spaltenmatrix $\{b_D\}$ mit der folgenden Gleichung ermittelt:

$$\{b_D\}=[A_B]^T\cdot\{\varphi_A\}+[A_D]\cdot\{\varphi_D\}$$ (4.70)

und dann das Gleichungssystem

$$[B]\{\sigma\}=\{b\}$$ (4.71)

löst. Die Matrix $[B]$ setzt sich aus den Elementmatrizen der Randelemente, die an den Oberflächen der Elektroden liegen, zusammen:

$$[B]=\sum_{m=1}^{M_B} [A^{e_m}_B] \,.$$ (4.72)

Bei der Summierung muss wieder von der lokalen Knotennummerierung der Elementmatrizen $[A^{e_m}_B]$ auf die globale Nummerierung übergegangen werden. Nummeriert man die Gitterknoten in einer Reihenfolge, sodass Knoten auf gleichen Elektroden durchlaufende Indizes haben, dann ergibt sich für die Matrix $[B]$ eine Blockstruktur.

$$\begin{bmatrix}[0]& [0] & [0] & [0] \\ {} [0] & [B_1] & [0] & [0... ...\{\begin{array}{c} \{0\} \\ \{b_1\} \\ \vdots \\ \{b_n\} \end{array}\right\}$$ (4.73)

Jeder der $n$ Blöcke kann einer Elektrode zugeordnet werden, und das Gleichungssystem kann nach Elektroden separiert gelöst werden.

Schwebende Randbedingungen

erfordern eine spezielle Behandlung, da hier weder das Potenzial der Elektrode noch die Verteilung der Ladung bzw. die entsprechenden Einträge im $\{b\}$-Vektor bekannt sind. Im Folgenden wird angenommen, die Knoten mit den globalen Indizes $k\ldots l$ seien Randknoten einer schwebenden Elektrode. Dann folgt aus (2.20) unmittelbar

$$\varphi_k=\varphi_i=\varphi_l \quad k<i<l.$$ (4.74)

Aus (4.61) sowie (4.48) lässt sich zeigen

$$\oint_{\Gamma _f}\tilde\sigma\,\textrm{d}A = \sum_{j=k}^l\sigma_i... ...sigma_i\sum_{i=k}^l\,\oint_{\Gamma _f}N_i N_j\,\textrm{d}A= \sum_{i=k}^{l} b_i.$$ (4.75)

Aus (2.21) resultiert somit die Forderung

$$\sum_{i=k}^{l} b_i=0.$$ (4.76)

Bei der Lösung des Gleichungssystems müssen nun (4.74) und (4.76) als zusätzliche Gleichungen in das System aufgenommen werden, oder man nutzt (4.74) um die Gleichungen $k$ bis $l$ zu eliminieren. Dazu ersetzt man in der Matrix $[A]$ die Zeilen $k$ bis $l$ durch eine einzelne Zeile mit den Spaltensummen der eliminierten Zeilen. Dasselbe wird natürlich auch für die Spalten von $[A]$ durchgeführt. Beim $\{b\}$-Vektor werden ebenfalls die Zeilen $k$ bis $l$ entfernt und durch eine einzelne Null ersetzt.

Die Feldenergie

kann mittels

$$W=\int_{\Omega}\!\mathchoice{\mbox{\boldmath$\displaystyle E$}} {... ...la N_i)\makebox{\boldmath$\underline\varepsilon$}(\nabla N_j)\,\textrm{d}\Omega$$ (4.77)

berechnet werden und ergibt

$$W=\frac12\{\varphi\}^T[A]\{\varphi\}=\frac12\{\varphi\}^T\{b\} =\frac12\{\varphi_D\}^T\{b_D\}\;.$$ (4.78)

Nimmt man wieder an, dass die Knoten, die zu einer Elektrode gehören, fortlaufend nummeriert sind so kann obige Gleichung in Teilsummen, die den einzelnen Elektroden zugeordnet sind, gegliedert werden

$$W=\varphi_1\sum_{j=K_1}^{K_2-1}b_j+\varphi_2\sum_{j=K_2}^{K_3-1}b_j\ldots$$ (4.79)

Die Potenzialwerte $\varphi_1, \varphi_2, \ldots$ sind auf den Elektroden konstant und wurden deshalb vor die Summe gehoben. Die Summierung über die $b_j$ entspricht gemäß (4.75) der Ladung auf der jeweiligen Elektrode. Die Energie des Feldes entspricht deshalb der Summe der Produkte aus Ladung und Potenzial der einzelnen Elektroden. Für die Kapazitätsberechnung mit Finiten Elementen bedeutet dies, dass die Methode der Ladungsintegration und die Energiemethode völlig äquivalent sind.

Leitungsströme

Die Berechnung der Potenzialverteilung und Stromdichte in elektrischen Leitern kann analog zur Berechnung des Elektrischen Feldes in dielektrischen Materialien erfolgen. Man erhält wieder ein Gleichungssystem der Form

$$[A]\{\varphi\}=\{b\}$$ (4.80)

mit

$$[A]=\sum_{m=1}^M [A^{e_m}_L] \quad\mathrm{und}$$ (4.81)

$${} \{b\}=\sum_{m=1}^{M_B} [A^{e_m}_B] \{J\}\,.$$ (4.82)

Die Berechnung der Elementmatrizen $[A^{e_m}_L]$ erfolgt mit (4.56) und (4.54), wobei für den Materialfaktor $\makebox{\boldmath $\underline m$}$ die elektrische Leitfähigkeit $\gamma$ eingesetzt wird. In (4.81) und (4.82) wird wieder bei der Summenbildung von lokalen auf globale Knotennummerierung übergewechselt. Die Matrix $[A]$ kann äquivalent zu (4.68) in die Teilmatrizen $[A_A]$, $[A_B]$, $[A_B]^T$ und $[A_D]$ aufgeteilt werden. Auch die Spaltenmatrix $\{b\}$ wird aufgeteilt in

$$\{b\}=\left\{\begin{matrix}\{b_A\}\\ \{b_D\}\end{matrix}\right\}$$ (4.83)

Neumann-Bedingungen

zeichnen sich dadurch aus, dass in der Spaltenmatrix $\{b_A\}$, die den entsprechenden Randknoten zuordenbare Zeilen, von Null verschiedene Einträge haben.

Die unbekannten Potenzialwerte $\{\varphi_A\}$ erhält man nun durch Lösen des Gleichungssystems

$$[A_A]\{\varphi_A\}=\{b_A\}-[A_B]\{\varphi_D\}.$$ (4.84)

Die Flächenstromdichte

$J$ auf Dirichlet'schen Rändern kann mit (4.70) bis (4.73) berechnet werden, indem formal $\sigma$ durch $J$ ersetzt wird.

Für schwebende Randbedingungen

gilt für das Potenzial (4.74) und für den Gesamtstrom in Analogie zu (4.76)

$$\sum_{i=k}^{l} b_i=I_f\,.$$ (4.85)

Diese Bedingungen können, wie zuvor bei der Berechnung des elektrostatischen Feldes, in das Gleichungssystem eingebaut werden, indem man die entsprechenden Zeilen und Spalten zusammenfasst.

Die Verlustleistung

ergibt sich als

$$W=\{\varphi\}^T[A]\{\varphi\}=\{\varphi\}^T\{b\}$$ (4.86)

Elektro-quasistatische Feldberechnung

Gleichung 4.37 kann in Matrixschreibweise in folgendermaßen angeschrieben werden:

$$[M_1]\{\dot\varphi\}=-[M_2]\{\varphi\}+[M_3]\{J\}.$$ (4.87)

Die Zeitdiskretisierung

wird vorgenommen, indem man eine Integration über einen Zeitschritt $\Delta t=t_n-t_{n-1}$ durchführt. (Zwecks Abkürzung wird für das Potenzial $\varphi(t_n)$ zum Zeitschritt $t_n$ einfach $\varphi_n$ geschrieben.)

$$[M_1]\bigl(\{\varphi_n\}-\{\varphi_{n-1}\}\bigr)= -[M_2]\int_{t_{... ...{t_n}\!\{\varphi\}\,\textrm{d}t + [M_3]\int_{t_{n-1}}^{t_n}\!\{J\}\,\textrm{d}t$$ (4.88)

Zur Lösung dieser Integrale sollen nun zwei verschiedene Näherungsverfahren betrachtet werden, nämlich das Rückwärts-Euler-Verfahren, das sich durch numerische Stabilität auszeichnet und die Trapezmethode, (oft auch als Crank-Nicolson-Verfahren bezeichnet) mit der im Allgemeinen eine höhere Genauigkeit erreicht werden kann. Beim Rückwärts-Euler-Verfahren wird das Integral über einen Zeitschritt durch die Fläche eines Rechtecks angenähert

$$\int_{t_{n-1}}^{t_n}\!\{\varphi\}\,\textrm{d}t \approx \Delta t\{\varphi_n\},$$ (4.89)

bei der Trapezmethode hingegen durch die Fläche eines Trapezes

$$\int_{t_{n-1}}^{t_n}\!\{\varphi\}\,\textrm{d}t \approx \Delta t\frac12 \bigl(\{\varphi_n\}+\{\varphi_{n-1}\}\bigr).$$ (4.90)

Setzt man nun eine der beiden Näherungsformeln in (4.88) ein, erhält man wiederum ein Gleichungssystem der Form

$$[A]\{\varphi_n\}=\{b\}$$ (4.91)

mit

$$\{b\}=[\bar A]\{\varphi_{n-1}\}+[B]\int_{t_{n-1}}^{t_n}\!\{J\}\,\textrm{d}t$$ (4.92)

mit dem man das Potenzial zum aktuellen Zeitpunkt $t_n$ ausgehend vom Potenzial zum vorhergegangenem Zeitpunkt $t_{n-1}$ berechnen kann. Dabei beginnt man zum Zeitpunkt $t_0$ mit dem Potenzial $\{\varphi_0\}$ laut Anfangsbedingung (2.41).

Für die Diskretisierung mittels Rückwärts-Euler-Verfahren gilt

$$[A] = [M_1]+\Delta t[M_2]= \sum_{m=1}^M \bigl([A^{e_m}_{L,\makebox{\b... ...erline\varepsilon$}}]+\Delta t [A^{e_m}_{L,\gamma}] \bigr) \qquad\mathrm{sowie}$$ (4.93)

$${} [\bar A] = [M_1]\phantom{+\Delta t[M_2]}\;= \sum_{m=1}^M [A^{e_m}_{L,\makebox{\boldmath$\underline\varepsilon$}}]\;.$$ (4.94)

Bei der Trapezmethode gilt

$$[A] = [M_1]+\frac{\Delta t}2[M_2] = \sum_{m=1}^M \bigl([A^{e_m}_{L,\m... ...epsilon$}}]-\frac{\Delta t}{2} [A^{e_m}_{L,\gamma}] \bigr) \qquad\mathrm{sowie}$$ (4.95)

$${} [\bar A] = [M_1]-\frac{\Delta t}2[M_2] = \sum_{m=1}^M \bigl([A^{e_m}_{L,\m... ...dmath$\underline\varepsilon$}}]+\frac{\Delta t}{2} [A^{e_m}_{L,\gamma}] \bigr).$$ (4.96)

In beiden Fällen erhält man die Matrix $[B]$ durch Summierung über alle Randelemente mit inhomogenen Neumann-Bedingungen

$$[B]=[M_3]=\sum_{m=1}^M [A^{e_m}_B].$$ (4.97)

Die Summenformeln in (4.93) bis (4.97) sind hier wieder nur symbolisch zu verstehen, da während der Summierung von der lokalen zur globalen Knotennummerierung gewechselt wird.

Die Matrix $[A]$ kann wieder, wie in (4.68), in die Teilmatrizen $[A_A]$, $[A_B]$, $[A_B]^T$ und $[A_D]$ aufgeteilt werden, ebenso die Spaltenmatrix $\{b\}$, wie in (4.83). Damit kann man die unbekannten Potenzialwerte $\{\varphi_A\}$ wie in (4.84) berechnen.

Wärmeleitung stationär

Das Gleichungssystem zur Berechnung der stationären Temperaturverteilung (4.39) kann in Matrixschreibweise als

$$[A]\{T\}=\{b\}$$ (4.98)

angeschrieben werden mit

$$[A]=\sum_{m=1}^M [A^{e_m}_{L,\gamma_T}]$$ (4.99)

und

$$\{b\}=\sum_{m=1}^M [A^{e_m}_P] \{p\} +\sum_{m=1}^M [A^{e_m}_B]\{\Theta\}\,.$$ (4.100)

Beim Bilden der Summen wird hier wieder von lokaler auf globale Knotennummerierung übergegangen.

Wählt man die Reihenfolge der Einträge der Spaltenmatrix der Temperatur wieder so, dass Knoten mit Dirichlet-Bedingungen am Schluss aufscheinen, kann man die Gleichung wie in (4.68) aufspalten und man erhält die unbekannten Temperaturwerte $\{T_A\}$ als Lösung des folgenden Gleichungssystems:

$$[A_A]\{T_A\}=\{b_A\}-[A_B]\{T_D\} \,.$$ (4.101)

Wärmeleitung transient

Gleichung 4.40 hat die Form

$$[M_1]\{\dot T\}=-[M_2]\{T\}+[M_3]\{p\}+[M_4]\{\Theta\}$$ (4.102)

Die Zeitdiskretisierung wird äquivalent zu (4.88)-(4.90) mit dem Rückwärts-Euler- oder Trapezverfahren vorgenommen. Man erhält wiederum ein Gleichungssystem der Form

$$[A]\{T_n\}=\{b\}$$ (4.103)

mit

$$\{b\}=[\bar A]\{T_{n-1}\} +\sum_{m=1}^M [A^{e_m}_P]\int_{t_n}^{t_... ...{d}t +\sum_{m=1}^M [A^{e_m}_B]\int_{t_n}^{t_{n-1}}\!\{\theta\}\,\textrm{d}t \,.$$ (4.104)

Beim Rückwärts-Euler-Verfahren erhält man die Matrizen $[A]$ und $[\bar A]$ mit

$$[A] =[M_1]+\Delta t[M_2] = \ \sum_{m=1}^M \bigl(c_p\rho_m[A^{e_m}_P] +\Delta t[A^{e_j}_{L,\gamma_T}]\bigr)$$ (4.105)

$${} [\bar A] =[M_1]\phantom{\Delta t[M_2] }= \sum_{m=1}^M c_p\rho_m[A^{e_m}_P]$$ (4.106)

und beim Trapezverfahren mit

$$[A] = [M_1]+\frac{\Delta t}2[M_2] = \sum_{m=1}^M \bigl(c_p\rho_m[A^{e_m}_P] +\frac12\Delta t[A^{e_m}_{L,\gamma_T}]\bigr)$$ (4.107)

$${} [\bar A] = [M_1]-\frac{\Delta t}2[M_2] = \sum_{m=1}^M \bigl(c_p\rho_m[A^{e_m}_P] -\frac12\Delta t[A^{e_m}_{L,\gamma_T}]\bigr)\,.$$ (4.108)

Im Fall, dass der thermische Leitwert $\gamma_T$ oder die Wärmekapazität $c_p$ einzelner Materialien nicht konstant ist, muss das natürlich bei der Assemblierung der Matrizen $[A]$ und $[\bar A]$ berücksichtigt werden. Und zwar müssen für $[A]$ die Werte zum Zeitpunkt $t_n$ und für $[\bar A]$ zum Zeitpunkt $t_{n-1}$ verwendet werden.

Da die elektrischen Zeitkonstanten in der Regel sehr viel kürzer sind als die thermischen, wird für die Berechnung des elektrischen Systems eine stationäre Stromverteilung angenommen. Unter der Voraussetzung, dass die elektrischen Parameter keine Temperaturabhängigkeit aufweisen, ist deshalb auch die Verlustleistungsdichte $\{p\}$ konstant und das Integral aus (4.104) vereinfacht sich zu

$$\int_{t_n}^{t_{n-1}}\!\{p\}\,\textrm{d}t=\Delta t\{p\} \,.$$ (4.109)

Ist hingegen etwa die elektrische Leitfähigkeit von der Temperatur abhängig, muss zur Bestimmung der Verlustleistungsdichte das elektrische System für jeden Zeitschritt neu berechnet werden und das Integral aus (4.104) kann mittels der Trapezformel genähert werden

$$\int_{t_n}^{t_{n-1}}\!\{p\}\,\textrm{d}t=\Delta t\frac{\{p_{n-1}\}+\{p_n\}}{2} \,.$$ (4.110)

Gleichung 4.103 kann dann wieder, wie in (4.68), aufgespalten werden und man erhält die unbekannten Temperaturwerte $\{T_A\}$ als Lösung des Gleichungssystems (4.101).