10 Numerische Feldberechnung
Elektrische, magnetische und elektromagnetische Felder technischer Komponenten und Gedite besitzen haufig nur geringe bzw. keine Symmetrie und entziehen sich meist einer mit vertretbarem Aufwand durchfiihrbaren analytischen Berechnung. RechnergestUtzte numerische Verfahren sind leistungsfahiger, doch verlangen ihre Komplexitat und die Vielfalt der problemspezifischen Algorithmen einen hohen Einarbeitungsaufwand. Dieses Kapitel versucht, den Leser am Beispiel der numerischen Berechnung elektrischer und elektromagnetischer Felder intuitiv in die Begriffswelt der numerischen Feldberechnung einzufUhren und ihm den Zugang zu den verschiedenen Verfahren zu erleichtem.
Von den heute benutzten Rechnerprogrammen werden folgenden Verfahren bzw. Methoden vorgestellt
- Finite-Elemente-Methoden, - Finite-Differenzen-Verfahren, - Ersatzladungsverfahren, - Boundary-Element-Methoden, - Momenten-Methoden, - Monte-Carlo-Verfahren.
Die meisten Verfahren wandeln die zugrunde liegende Funktionalgleichung (Gleichung, bei der auf mindestens einer Seite mindestens ein Funktional,
A. J. Schwab et al., Begriffswelt der Feldtheorie© Springer-Verlag Berlin Heidelberg 1998
194 10 Numerische Feldberechnung
das heiBt, eine "Funktion einer Funktion" steht) in eine Matrixgleichung bzw. ein lineares Gleichungssystem urn, das auf einem Rechner gelOst wird. Als Funktionalgleichung kommen in der Elektrostatik am haufigsten die Laplace- oder Poisson-Gleichung sowie die Partikularlosung letzterer, das sogenannte Poisson-Integral in Frage (s. 4.2.3),
~<p(r) = -p(r) e
Laplace-Gleichung Poisson-Gleichung Poisson-Integral (10-1)
Bei den meisten technischen Fragestellungen handelt es sich urn Randwertprobleme, die durch die Laplace-Gleichung beschrieben werden, wobei die "Kunst" im Einbringen der Randbedingungen besteht. Haufig erfolgt die Beriicksichtigung der Randbedingungen auch dadurch, daB man das Randwertproblem in ein Poisson-Problem umwandelt (s. 2.1), mit anderen Worten, sich als Ursache des Feldes doch eine bestimmte Linien- oder Flachenladungsdichte p(r), diesmal jedoch auBerhalb der Berandung vorstellt (z.B. Ersatzquellenverfahren mit Ladungen und Stromen innerhalb von Leitern). Das Feldproblem wird dann durch eine Poisson-Gleichung beschrieben, deren homogene Gleichung (Laplace-Gleichung) die Triviallosung Null besitzt.
Die folgenden Abschnitte stellen verschiedene Algorithmen vor, die beispielhaft die Konversion von Funktionalgleichungen vom Typ (10-1) in lineare Gleichungssysteme demonstrieren.
10.1 Finite-Elemente-Methode
Bei der Finite-Elemente-Methode unterscheidet man innerhalb des Verfahrens je nach Art der Ermittlung sogenannter Elementgleichungen grob zwischen Direkter Methode (Kollokationsmethode), Variationsrechnung (Raleigh-Ritz), Gewichteten Residuen (Galerkin et al) und Energiebalance. Fur die Berechnung von Potentialfeldern kommt die Variationsrechnung und die Methode der gewichteten Residuen zur Anwendung, wobei wir den Schwerpunkt auf die Variationsrechnung legen.
10.1 Finite-Elemente-Methode 195
Wahrend in der gewohnlichen Extremwertrechnung gefragt wird
- fur welche Werte Xv eine Funktion f(x) ein Maximum oder Minimum aufweist,
stellt sich in der Variationsrechnung die Frage
- fur welche Funktionen f(x) ein Funktional X(f(x») ein "Maximum oder Minimum" aufweist.
In beiden Fallen erhalt man die Losung durch Differenzieren und Nullsetzen des Differentialquotienten.
Ein von der Physik nahegelegtes Funktional des elektrostatischen Felds ist die im Feldraum gespeicherte elektrische Energie. Bekanntlich verteilen sich auf Elektroden aufgebrachte Ladungen derart, daB sich ein Minimum der potentiellen elektrischen Energie einstellt (ahnlich wie Wasser bergab flieBt, bis sich ein ebener Wasserspiegel eingestellt hat).
a) b)
Bild 10.1: Volumenelement!:N bzw. M!1z eines elektrostatischen Felds, a) dreidimensionales Feld, b) zweidimensionales Feld, z-Achse senkrecht zur Zeichenebene.
Die in einem Volumenelement !lV, s. Bild 10.la, gespeicherte Energie berechnet sich mit der Energiedichte
(10-2)
gemaB Abschnitt 6.3.4 zu
(10-3)
196 10 Numerische Feldberechnung
bzw. mit E(r) = -gradcp(r) = -dcp/dr zu
WIN = .!.Ef (dCP)2 dV = .!.Ef ((dCP)2 + (dCP)2 + (dCP)2JdV 2 /J"V dr 2 /J"V dx dy dz
(10-4)
Zur Vereinfachung beschdinken wir uns auf ein zweidimensionales, yon der z-Koordinate unabhangiges Feld, Bild 10.lb, und erhalten
(10-5)
bzw. bezogen auf die Langeneinheit
W/J"v = .!.Ef ((dCP)2 + (dCP)2JdA !lz 2MdX dy
(10-6)
In (10-6) bezeichnen wir W/J"Y/!:J.z als Funktional X M =f(CPM(X,y» eines Elements !lA, die Summe aller Elementfunktionale "i.W/J"v/!lz als Funktional X = f(cp(x,y») des gesamten Feldgebiets.
Nach diesen yorbereitenden Ausfiihrungen konnen wir uns dem eigentlichen Finite-Elemente-Verfahren zuwenden, das in folgenden Schritten ablauft
1. Diskretisierung des Feldgebiets in einzelne Elemente,
2. Approximation des Potentials CPM(x,y) innerhalb eines Elements,
3. Ermittlung der Elementgleichungen (~Elementmatrizen),
4. Ermittlung der Systemgleichung (~Systemmatrix) (Assemblierung der Elemente zum System),
5. Berucksichtigung der Randbedingungen und Auflosurtg des linearen Gleichungssystems.
10.1 Finite-Elemente-Methode 197
1m folgenden werden die einzelnen Schritte naher betrachtet.
1. Diskretisierung:
Zunachst unterteilt bzw. approximiert man das Feldgebiet durch n Elemente variabler GroBe und Form, z.B. Dreiecke im zweidimensionalen, Tetraeder im dreidimensionalen Fall, Bild 10.2.
Bild 10.2: Diskretisierung eines zweidimensionalen Feldgebiets durch Dreieckselemente.
Urn numerische Probleme zu vermeiden, sollten die Elemente keine extrem stumpfen oder spitzen Winkel enthaIten, m.a.W. moglichst gleichseitig sein. Eine Diskretisierung mit krummlinigen Elementen ist ebenfalls moglich. Sie erlauben bei gleicher Ergebnisgenauigkeit mit kleinerer Elementzahl auszukommen und verringern damit den Eingabeaufwand.
Die gewohnliche Finite-Elemente-Methode eignet sich nur fur geschlossene Feldgebiete. Abhilfe schaffen moderne Verfahren, deren Elementmatrix einen offenen Rand zu beriicksichtigen gestatten (z.B. sog. Hybridverfahren).
2. Approximationsfunktion tpM(x,Y) innerhalb eines Elements:
Die Approximation der unbekannten Potentialfunktion <PM (x,y) innerhalb eines Elements erfolgt im einfachsten Fall durch einen linearen Ansatz,
(10-7)
in dem C1,C2,C3 a priori unbekannte Koeffizienten sind. Ansatze hoherer Ordnung sind moglich bzw. werden z.B. bei isoparametrischen krummlini-
198 10 Numerische Feldberechnung
gen Elementen zwingend erforderlich. Der Einfachheit halber beschranken wir uns im folgenden auf einen linearen Ansatz.
Die drei Koeffizienten C.1,C2,C3 Iassen sich nach Einsetzen der drei Wertepaare Xi Yj, Xj Yj, Xk Y k, auch durch die Knotenpotentiale <Pj, <Pj, <Pk des Elements (sog. generalisierte Koordinaten) ausdriicken,
<Pi = <P(xi, Yi) = C1 + C 2xi + C3Yi ) 3 Gleichungen fur <Pj = <p(Xj,Yj) =C I +C 2Xj +C3Yj 3 Unbekannte
<Pk = <P(xk,Yk) =C I +C 2Xk +C 3Yk CI ,C2,C3 (10-8)
Nach Aufiosen obigen Gleichungssystems erhalt man die Koeffizienten C I ,
C2, C3 als Funktion der Knotenpotentiale und -koordinaten,
(10-9)
Setzt man die Koeffizienten in den urspriinglichen Ansatz (10-7) ein, erhalt man nach Umformung eine zweite Version der Approximationsfunktion, in der das Potential <PM(x,y) durch Interpolation aus jeweils drei Knotenpotentialen <Pi' <Pj' <Pk gewonnen wird,
<PM (x,y) = Ni(x, y) <Pi + Nj(x,y)<pj + Nk(X,Y)<Pk = I. Ny<py y
(v = i,j,k) . (10-tO)
Die Funktionen Nj, Nj, Nk werden Interpolations- bzw. Formfunktionen genannt und hangen von der Form und Lage eines Elements abo Sie sind eine Funktion der Knotenkoordinaten des jeweiligen Elements und demnach von Element zu Element verschieden. Fur einen linearen Ansatz berechnen sich die Nv nach folgender Beziehung
(10-11)
wobei die Indizes zyklisch vertauscht werden.
10.1 Finite-Elemente-Methode 199
BezUglich des Begriffs generalisierte Koordinaten sei daran erinnert, daB man zur eindeutigen Bestimmung der Lage eines Systems von n Massepunkten 3n Koordinaten - fUr jeden Freiheitsgrad eines Massepunkts eine - benatigt. Auch jedes beliebigt:, nicht mechanische System mit K Freiheitsgraden erfordert K GraBen zur eindeutigen Beschreibung seines Zustands. Diese GraBen mUssen nicht notwendigerweise Koordinaten im elementargeometrischen Sinn sein, wie beispielsweise die kartesischen Koordinaten x, y, Z, sondem kannen sich auch als ZustandsgraBen, z.B. Knotenspannungen (bzw. Knotenpotentiale) eines elektrischen Netzwerks usw. oder gar GraBen ohne unmittelbare physikalische Bedeutung manifestieren. Man bezeichnet sie dann oberbegrifflich als verallgemeinerte (generalisierte) Koordinaten, ihre ersten Ableitungen als verallgemeinerte Geschwindigkeiten.
3. Ermittlung der Elementgleichungen bzw. Elementmatrix:
An dieser Stelle unterscheiden sich die eingangs erwahnten Finite-ElementeVersionen durch die Art und Weise, wie man zu den Elementgleichungen gelangt. FUr die vorliegende Aufgabe bedienen wir uns der oben bereits erwahnten Variationsrechnung. Zunachst differenzieren wir den Approximationsansatz <PM (x,y) partiell nach x und y,
(10-12)
und setzen die so erhaltenen Ausdriicke in das Funktional
(10-13)
eines Elements ~A ein. Auf diese Weise erhalten wir das Elementfunktional als eine Funktion der Knotenpotentiale eines Elements "ijk",
(10-14)
200 10 Numerische Feldberechnung
Wie bei der gewohnlichen Extremwertrechnung bilden wir jetzt den Differentialquotienten
dXM, ! 0 bzw. partiell: d{<pM,} ,
(10-15)
Die partielle Differentiation ergibt drei Bestimmungsgleichungen fur die drei unbekannten Potentiale <Ph <Pj, <Pk bzw. ein lineares Gleichungssystem
(10-16)
Fili sich alleine betrachtet besitzt dieses Gleichungssystem nur triviale Losungen <Pi = 0; <Pj = 0; <Pk = O. Eine nichttriviale Losung existiert nur dank der Kopplung aller Gleichungssysteme durch aneinandergrenzende Elemente, insbesondere der an den Rand angrenzenden Elemente.
Die Koeffizienten der Matrix sind Funktionen der Knotenkoordinaten und der PermittivWit £r der Elemente, d.h.
Bei einheitlichem Dielektrikum kann £r vor die Matrix gezogen werden.
Aufgrund ihrer Dimension wird die Matrix Permittivitiitsmatrix bzw. im 3-dimensionalen Fall Kapazitiitsmatrix genannt (im letzteren Fall bleibt az auf der rechten Gleichungsseite); sie entspricht der Steifigkeitsmatrix in der Mechanik. Jedes Element besitzt seine eigene Permittivitatsmatrix, in der jede Zeile einen Elementlmoten und seine VerknUpfung mit den beiden anderen Elementen beschreibt.
4. Ermittlung der Systemgleichung bzw. Systemmatrix:
In diesem Schritt erfolgt die Assemblierung der n Elementmatrizen zur vollstiindigen Systemmatrix. Ausgehend von der Tatsache, daB die gesamte Feld-
10.1 Finite-Elemente-Methode 201
energie gleich der Summe aller Elementenergien ist, bestimmt man das Systemfunktional X = f(<p(x,y)) aus der Summe aller Elementfunktionale XM = f( <PM (x,y)),
(10-17)
Die partielle Differentiation des Systemfunktionals liefert nach Nullsetzen,
ax ~o a{<p}
(10-18)
n Bestimmungsgleichungen fUr n Knotenpotentiale <Pl ... <Pn bzw. das line are G leichungssystem
[P]{<p} = 0 (10-19)
Die Koeffizienten der Systemmatrix erhalt man durch Summation zusammengehoriger Koeffizienten der Elementmatrizen. Jede Zeile der Systemmatrix beschreibt einen Elementknoten und seine VerknUpfung mit den Elementknoten aller anderen anstoBenden Elemente
(10-20)
Betrachten wir beispielsweise Knoten 10 in Bild 10.3, so ergibt sich ein Diagonalelement der Systemmatrix, z.B. PlO,l2' zu
PlO,lO = PlO,lO(l) + PlO,lO(2) + PlO,lO(3) + ... + PlO,lO(6) ,
ein Nichtdiagonalelement der Systemmatrix, z.B. PlO,l2 zu
P lO,l2 = PlO,l2(2) + PlO,l2(3)
202 10 Numerische Feldberechnung
Element 3 Element 2
[PlO,10
Pl1,lO
P12,10
PlO,ll
Pl1,l1
P12,11
PlO,12j {«>1O} Pll,12 «>11 = 0
P12,12 «>12
[PlO,10
P12,10
P13,10
PlO,12
P12,12
P13,12
PlO,13j {«>1O} P12,13 «>12 =0
P13,13 «>13
rp1,1,,_ '"
- "'---"" ---" "
<PI
o
~'---------"---------~,----, , , ',P P P p', P P P',
10,7 10,8 10,9 10,10 10,11 10,12 10,13 ,
------~----------~--------~ " "
<P1O = 0
" " " " " " o
Pn,1 " ~n " ..
Bild 10.3: Assemblierung der Systemmatrix aus Elementmatrizen, gezeigt am Beispiel des Knotens 10. Die Elementgleichungssysteme besitzen eine nichtriviale Losung nur dank ihrer Kopplung untereinander, z.B. durch das gleichzeitige Vorkommen der Potentiale <P1O und <P12 in den Gleichungssystemen der Elemente 3 und 2.
10.1 Finite-Elemente-Methode 203
Da nicht jeder Knoten mit jedem anderen Knoten verbunden ist, wird die Matrix nur schwach besetzt, z.B. tiber 99% Nullelemente. Durch geschickte Numerierung im Rahmen der Behandlung schwach besetzter Matrizen (engl.: sparse matrix tec/1niques) HiBt sich erreiche~, daB die wenigen von Null verschiedenen Elemente in einem symmetrisch zur Diagonalen verteilten schmalen Band liegen (Bandmatrix), was eine kompakte Speicherung ermoglicht.
5. Beriicksichtigung der Randbedingungen:
Die im vorigen Schritt aufgestellte Systemmatrix [P] ist wie die Knotenadmittanzmatrix [Y] der Netzwerktheorie (z.B. bei der LeistungsfluBrechnung in Elektroenergiesystemen) symmetrisch und singular. Mit anderen Worten, ohne Einfiihrung von Randbedingungen besitzt das Gleichungssystem nur die Triviallosung Null. Die Berticksichtigung der Randbedingungen bei der numerischen Losung partieller Differentialgleichungen entspricht der Festlegung der beim analytischen Losen eines unbestimmten Integrals bzw. einer Differentialgleichung auftretenden Integrationskonstanten (vgl. 4.4.1 und 9.5).
Da man die Potentiale der Randknoten kennt - im einfachsten Fall 0% und 100% - konnen die Zeilen der Randknoten gestrichen werden. Mit dem Verschwinden der Randpotentiale im Knotenpotentialvektor wtirden aber auch gleichzeitig aIle Produkte, die Randpotentiale als Faktoren aufweisen, Pmn<!>Rand, verschwinden. Urn dies zu vermeiden, bringt man diese Produkte auf die rechte Gleichungsseite - auf diese Weise wird auch der beim Streichen der Zeilen mit Randpotentialen verlorengegangene quadratische CharaIder der Matrix wieder hergestellt - und erhalt folgendes Gleichungssystem
[:P]{ ip} = [P]Rand {<!> }Rand (10-21)
dessen rechte Seite bekannt ist.
Die Auflosung dieses Gleichungssystems nach dem Knotenpotentialvektor {ip} liefert die gesuchten Potentialwerte in den Elementknoten im Feldraum,
(10-22)
204 10 Numerische Feldberechnung
In praxi wird das Gleichungssystem nicht durch Inversion der Koeffizientenmatrix, sondern durch den Gauf!, Algorithmus bzw. iterativ gelost.
Nach Erhalt aller Knotenpotentiale konnen die Potentialfunktionen CPM(x,y) und die FeldsHirke E = -grad cP in den Elementen ermittelt werden. Hier ist ausdriicklich zu betonen, daE die Aufiosung des Gleichungssystems lediglich diskrete Werte der Potentiale bzw. Feldstarken liefert. Die Darstellung des Ergebnisses in Form von Feldbildern usw. erfordert zusatzliche Zeichenprogramme, die z.B. nach Ermittlung von Punkten gleichen Potentials (Stiitzstellen) diese durch kontinuierliche Linien verbinden, beispielsweise mit Splinefunktionen. Letztere sind stiickweise definierte Funktionen, die an den StoEstellen nicht nur gleiche Funktionswerte, sondern auch gleiche erste Ableitungen besitzen und daher eine knickfreie Interpolation ermoglichen.
1m Fall des hier gewahlten einfachen linearen Ansatzes ergibt sich die Feldstarke innerhalb eines Elements als konstanter Wert. Ansatze hoherer Ordnung ermoglichen deutlich besser annahernde Ergebnisse fur den erwarteten Feldverlauf, Bild 10.4. Alternativ kann man auch von Anfang an feiner diskretisieren.
Bild 10.4: Aquipotentialliniendiagramm bei linearem (links) und quadratischem Ansatz (rechts) der Approximationsfunktion.
10.1 Finite-Elemente-Methode 205
SchlieBlich sei erwahnt, daB man mit Mehrfeldfunktionalen gleichzeitig mehrere Knotenvariable, z.B. Potential und Feldstarke in den Knoten, berechnen kann.
Sollen an Stelle reiner Potentialfelder nichtstationiire elektromagnetische Felder berechnet werden, ist obiger Ansatz der Minimierung der elektrostatischen Energie auf alle Energiebeitrage des von den Randem eingeschlossenen Volumens auszudehnen, d.h.
- gespeicherte elektrische Energie:
(10-23)
- gespeicherte magnetische Energie:
(10-24)
- tiber eingepragte Leitungsstrome zugefiihrte Energie:
(10-25)
- Energie einer Raumladungsdichte p(r):
w = r pcp dV P Jv 2 '
(10-26)
- Wirbelstrom-Verlustenergie:
(10-27)
206 10 Numerische Feldberechnung
- Zusatzterm (Komplementare Energie):
(10-28)
Wie beim elektrostatischen Feld werden auch hier die Feldstiirken durch Potentiale ersetzt, wobei neben dem Skalarpotential <p zusatzlich das magnetische Vektorpotential A (s. a. 5.3) einzufiihren ist. AIle Gleichungen - bis auf den ersten Teil der Gleichungen (10-23), (10-24) und (10-25) - besitzen die allgemeine, auch fUr nichtlineare Probleme giiltige Form. Wie im Kapitel 6 ausfiihrlich erlautert wurde stellt sich bei den Potentialen <p und A im elektromagnetischen Feld die Frage der Eichung (s.a. A6). Da hieriiber in obigen Gleichungen keine Aussagen gemacht werden fiihrt man den Zusatzterm (10-28) ein. Er reprasentiert eine von div A herriihrende Zusatzenergie.
AnschlieBend bildet man aus der Summe aller Energiebeitrage das elektromagnetische Energiefunktional
differenziert dieses nach <p und A, und setzt gleich Null,
of =0 o(<p,A)
(10-29)
Auf ahnliche Weise, wie weiter oben ausfUhrlich fUr das elektrostatische Feld gezeigt, folgt schlieBlich die vollstandige Matrizengleichung:
[M(e)] [~~ ]+[ L(O)][!~ ]+[ K(~)]r:] =~] (10-30)
10.1 Finite-Elemente-Methode
In 10-30 bedeuten M = Permittivitatsmatrix, L = Leitfahigkeitsmatrix, K = Permeabilitatsmatrix.
207
Die Losung des zugehorigen linearen Gleichungssystems mittels problemspezifischem Lasungsverfahren fuhrt auf die gesuchten GraBen <p und A, aus denen sich wiederum die elektrische Feldstarke E und die magnetische Feldsmrke H tiber
bestimmen lassen.
aA E = -grad <p - -
at und H=~rotA
Jl
Offensichtlich mussen bei zeitlich veranderlichen Feldem neben den bekannten Randbedingungen der quasistatischen Fragestellung auch Anfangswerte berucksichtigt werden, das heiBt es handelt sich hier urn kombinierte AnfangswertiRandwertprobleme.
Nach Streichen aller zeitabhangigen Terme geht (10-30) in die Gleichungen fur das elektrostatische und magnetostatische Feld sowie das GleichstromStromungsfeld tiber.
Methode der gewichteten Residuen:
AbschlieBend sei kurz die Vorgehensweise bei der Erstellung der Elementgleichungen nach der Methode der "Gewichteten Residuen" angedeutet.
Die Approximationsfunktion (10-10) erfiHlt die Laplace-Gleichung Ll<p=O in einem Element M nur bis auf eine endliche Abweichung, die man als Residuum R bezeichnet
(10-31)
Grundsatzlich wiinschte man sich, daB das Residuum tiberall im Element verschwindet, was jedoch wegen der endlichen Anzahl der Koeffizienten C1,
C2, C3 im allgemeinen nicht moglich sein kann. Man wird bei einer Nahe-
208 10 Numerische Feldberechnung
rungslOsung aber auch zufrieden sein, wenn das Residuum an moglichst vielen Stellen verschwindet, bzw. tiber ein Element gemittelt, minimal wird,
(10-32)
Die Multiplikation des Residuums mit einer geeigneten Gewichtsfunktion W erlaubt sogar, den gemittelten Fehler exakt zu Null zu machen, das heiBt
(10-33)
Man erhiilt so eine Bestimmungsgleichung fUr die einzelnen Elementgleichungen, die anschlieBend noch zur System-Matrixgleichung assembliert werden, wie das oben fUr die Variationsrechnung gezeigt wurde.
Abhiingig von der Wahl der Gewichtsfunktion erhaIt dann die "Methode der gewichteten Residuen" besondere Namen. Beim Galerkin-Verfahren wiihlt man als Gewichtsfunktionen die Formfunktionen Nv , das heiBt
(10-34)
die letztlich auf das gleiche Gleichungssystems fUhren wie die oben beschriebene Variationsrechnung.
Bei der Kollokationsmethode bzw. Point Matching Methode wiihlt man als Gewichtsfunktion den Dirac-Impuls, z.B. beim Ersatzladungsverfahren, (vgl. a. 4.1).
(10-35)
Beim Verfahren der "Minimalen Fehlerquadrate" multipliziert man das Residuum schlieBlich mit sich selbst.
10.2 Differenzenverfahren 209
10.2 Differenzenverfahren
Das Differenzenverfahren erlaubt die naherungsweise Losung der Laplaceschen Potentialgleichung durch Uberfiihren der Differentialquotienten in Differenzenquotienten. Es eignet sich fUr zwei- und dreidimensionale, iiberwiegend geschlossene Gebiete. Die Vorgehensweise wird hier am zweidimensionalen Beispiel gezeigt. 1m einfachsten Fall wird der Feldraum mit einem aquidistanten Punktgitter iiberzogen, Bild 10.5a.
® ® ® CD
(7)
/ "" ® ®
®
®, hwO thN CD ~ \
hs{ '----.r--'
hE
@
a) b)
Bild 10.5: Zur Erliiuterung des Differenzenverfahrens. a) Diskretisierung des Feldgebiets durch ein Punktgitter, b) Bezeichnung der Gitterabstande fiir h"* const.
Approximiert man die unbekannte Potentialfunktion <p(x,y) im Knoten 0 durch eine nach der zweiten Ableitung abgebrochene Taylorreihe
(10-36)
210 10 Numerische Feldberechnung
so ergeben sich die Potentiale der Nachbarknoten PI bis P4 nach Einsetzen ihrer Koordinaten zu
CPI = cp(xo +h,yo) = CPo + h acp + 0 + h2 a2cp
0 -- + ax 2 ax2
CP2 = cp(xo, Yo +h) = CPo + 0 + h acp ay
+ 0
CP3 = cp(xo -h,Yo) = CPo h acp + 0 + h2 a2cp 2 ax2 ax
CP4 = cp(xo,Yo -h) = CPo - 0
mith=x-xo=y-yo .
acp - h- + 0
ay
Nach Addition aller vier Gleichungen erhiilt man
+ 0
+ 0
+ 0
+ 0
+ h 2 a2cp 2 ay2
+ 0
h2 a2cp + --
2 ay2
(10-37)
Der Term in der Klammer nimmt definitionsgemiiB wegen L\cp = 0 den Wert Null an, so daB man das Potential im Punkt "0" durch die vier Potentiale der Nachbarpunkte darstellen kann,
(10-38)
An den Elektrodenriindern liiBt sich die Aquidistanz des Punktgitters meist nicht aufrechterhalten, d.h. h":F- const. Man arbeitet daher mit der allgemein giiltigen Randgebietsformel, die man aus obigen Gleichungen fUr h":F-const und Ersetzen der Differentialquotienten durch Differenzenquotienten
L\x = hw bzw. hE und L\y = hs bzw. hN {lO-39)
erhiilt,
10.2 Differenzenverfahren 211
~ + ~ + ~ + ~ _ hE {hE +hW) hN{hN +hs) hw{hw +hE) hs{hs +hN)
~o- 1 1 --+--
'-_________ h..:..:W_h-=.E_h....:N-'-h-=-s _______ --" . (10-40)
Der Anschaulichkeit wegen sind hier die Distanzen zu den N achbarpunkten mit Indizes der vier Himmelsrichtungen in englischer Notation bezeichnet (zur Vermeidung einer Verwechslung 0 fur Ost mit 0 fUr Null, vgl. Bild 10.5b). Ausgehend von den Randpotentialen konnen nun die Potentiale der Gitterlmoten im Feldraum Punkt fUr Punkt ermittelt werden.
Die Gleichungen aller Gitterpunkte bilden ein lineares Gleichungssystem, das durch geeignete Umformung der Randgebietsformel erhalten wird. Man erweitert zunachst mit dem Nenner der rechten Seite und bringt alle Terme auf eine Gleichungsseite,
bzw.
+ + +
Fur die Gitterpunkte 0 bis n erhalt man somit folgende Gleichungszeilen
0: hO~ ~o + hOE ~1 + hON ~2 + how ~3 + hos ~4 . = 0
1: hlO~O+hll~l' +h15~5+h16~6+h17~7 =0
In Matrixschreibweise
r~oo . ~onl {~o} = 0
h no h nn ~n
212 10 Numerische Feldberechnung
bzw. I [h]{q>} = 0 (10-41)
Da die Randgebietsformel jeweils nur 5 Gitterpunkte miteinander verkntipft, ist die Matrix nur schwach besetzt, z.B. tiber 99% Nullelemente. Durch geschickte Numerierung im Rahmen der Behandlung schwach besetzter Matrizen (engl.: sparse matrix techniques) HiRt sich erreichen, daR die von Null verschiedenen Elemente in einem schmalen Band urn die Hauptdiagonale liegen (Bandmatrix), was eine kompakte Speicherung ermaglicht.
Wie beim Finite-Elemente-Verfahren ist die Matrix zunachst singular, das Gleichungssystem wird erst nach Einfuhren der Randbedingungen eindeutig lasbar (vgl. 10.1). Man streicht zunachst wieder die Zeilen der Randknoten, deren Potential ja bekannt ist, und wahrt den quadratischen Charakter der Matrix, indem man die Produkte hmnq>Rand auf die rechte Gleichungsseite bringt. Dies ftihrt auf das lineare Gleichungssystem
[h] {ip} = [h ]Rand {q> }Rand (10-42)
mit bekannter rechter Seite. Die Aufiasung dieses Gleichungssystems nach dem Vektor {ip} liefert die gesuchten Potentiale der diskreten Gitterpunkte. Potentiale in den von Gitterpunkten eingeschlossenen Gebieten mtissen durch Interpolation ermittelt werden.
Die Erweiterung auf dreidimensionale Probleme ist formal einfach, indem man ausgehend von der Taylorreihe fur drei unabhangige Variable die sogenannte Wiirfelformel herleitet,
(10-43)
Ahnlich wie sich die Finite-Elemente-Methode auf die Behandlung elektromagnetischer Wellenprobleme erweitern lieR (s. 10.1) erlaubt auch das Differenzenverfahren eine Erweiterung seiner Funktionalitat auf Wellenprobleme (Anfangswert-Randwertprobleme), sog. Finite-Differenzen Zeitbereichsmethode (engl.: FD-TD Method - Finite-Difference Time-Domain Method). Insbesondere eignet es sich zur Lasung elektromagnetischer Pro-
10.2 Differenzenverfahren 213
bleme unmittelbar im Zeitbereich (Streuung transienter elektromagnetischer Felder etc.), was im folgenden noch kurz angedeutet werden solI.
Die Maxwellschen Gleichungen in Differentialform,
rotE aB at
an rotH =JL +
at (10-44)
fuhren, beispielsweise fur ein Kartesisches Koordinatensystem, auf folgendes System skalarer Differentialgleichungen,
a) _ aB x = aEz _ aEy d) aD x = aHz _ aHy _ J at ay az at ay az x
b) _ aBy = aE x _ aEz e) aDy = aH x _ aHz _ J at az ax at az ax y
c) _ aB z = aEy _ aE x f) aDz = aHy _ aH x _ J (10-45 a-f) at ax ay at ax ay z
Diskretisieren wir den Feldraum durch ein kubisches Gitter mit den Gitterabstanden ~,!l.y, !l.z, dann wird ein bestimmter Raumpunkt mit den Laufvariablen i,j,k gekennzeichnet durch P(i~, j!l.y, k!l.z). Da wir zeitlich veranderliche GraEen betrachten und auch die Zeit diskretisieren, werden die FeldgraEen, z.B. Bz , beschrieben durch
(10-46)
Mit dieser Vereinbarung laEt sich z.B. die Differentialgleichung (10-45a) als Differenzengleichung darstellen,
B n+ 112 (i J' +:.. k +:..) - B n -112 (i J' + ~ k +:..) x '2' 2 x '2' 2
!l.t
_ EV(i,j +~,k +l)~EV(i,j +~,k) Er(i,j +l,k +~)-Er(i,j,k +~) !l.z !l.y
(10-47)
214 10 Numerische Feldberechnung
Beginnend mit t=O (n=O) konnen B bzw. H und E alternierend nach jeweils der halben Zeitschrittweite unter Verwendung zweier urn die halbe Schrittweite gegeneinander versetzten kubischen Gittern rekursiv berechnet werden. Bei diesem Algorithmus findet daher keine Matrixinversion statt, es werden lediglich Matrizen mit Vektoren aktualisiert die den Feldzustand zum jeweiligen Zeitpunkt darstellen.
Die Randbedingungen werden dadurch berlicksichtigt, daB bei Annahme eines perfekten Leitermaterials die Tangentialkomponenten der elektrischen Feldstarke und die Normalkomponenten der magnetischen Feldstarke im Gleichungssystem (10-45) verschwinden.
Flir die Diskretisierung gelten folgende Gesichtspunkte:
- Die Abstande Ax, !J.y, !J.z mlissen "elektrisch kurz " , das heiBt klein gegen die Wellenlange (vgl. Kapitel 8) und auBerdem klein gegen die Abmessungen des streuenden Objekts sein.
- Der "Abstand" ~!J.x2 + !J.y2 + !J.Z2 muB groBer sein als der wahrend eines Zeitschritts zuruckgelegte Weg !J.s = c !J.t.
Da man nicht den gesamten Raum diskretisieren kann, wird in akzeptabler Entfernung yom Streuobjekt (z.B. lOAx) eine klinstliche Berandung eingefligt (engl.: lattice truncation), deren Randbedingungen so gewahlt werden, daB der unendlich ausgedehnte Feldraum optimal simuliert wird. Haufig stehen mehrere Moglichkeiten zur Wahl, zum Beispiel Strahlungsbedingung, Absorptionsbedingung oder Lattice-Truncation Bedingung, deren ausflihrliche Erlauterungen jedoch der Spezialliteratur vorbehalten bleiben muE.
10.3 Ersatzladungsverfahren
Beim Ersatzladungsverfahren wird eine Konfiguration von Ersatzladungen (engl.: simulation charges) ermittelt, deren PotenHalfunktion <rE(r) die tatsachliche Potentialfunktion <r(r) in der Umgebung einer an Spannung liegenden technischen Elektrode (Elektrodenpotential <rEi) approximiert. Urn die gesuchte Konfiguration zu find en, positioniert man im Innern der Elektrode an yom Anwender auszuwahlenden Orten n unbekannte Punktladungen Ql . . . Qn' Bild 10.6.
10.3 Ersatzladungsverfahren
Bild 10.6: Aufgeschnittene toroidfi::irmige Ringelektrode (z.B. Steuerelektrode eines StolSspannungsteilers), Positionierung von 6 Ersatzladungen Ql".Q6 auf der zentralen Achse eines Toroids.
215
Mit Hilfe des Superpositionsprinzips berechnet man nun das Potential <PI in einem Punkt PI auf der Elektrodenkontur (s. 4.3),
(10-48)
oder mit Potentialkoeffizienten abgekurzt dargestellt,
6
<PI = L <Pli = Pll QI + Pl2 Q2 + .... + Pl6 Q6 (10-49) i=l
Waren die sechs Ladungen QI ... Q6 bekannt, konnte man aus der Geometrie das Potential <PI berechnen. Bei technischen Feldberechnungsproblemen sind aber nicht die Ladungen bekannt, sondern die Potentiale der Rander (Randwertproblem). Deshalb lautet beim Ersatzladungsverfahren die Fragestellung gerade umgekehrt: Bestimme zu einem bekannten Potential <PI = <PEl sechs Ersatzladungen QI'" Q6, deren Teilpotentiale in ihrer Gesamtheit mit dem vorgegebenen Potential in PI, <PI = <PEl ubereinstimmt. Da wir bei dieser Fragestellung bis jetzt nur eine Gleichung fur sechs Unbekannte haben, schreiben wir obige Gleichung flir flinf weitere Konturpunkte P2",P6 mit einheitlichem bekanntem Potential <PEl an und erhalten folgendes Glei-chungssystem
PllQI +P12Q2 +P13Q3 +P14Q4 +PISQS +P16Q6 = <PEl
P21QI +P22QZ +P23Q3 +P24Q4 +P2SQS +P26Q6 = <PEl
P3lQI +P32Q2 +P33Q3 +P34Q4 +P3sQs +P36Q6 = <PEl
P4lQI +P42Q2 +P43Q3 +P44Q4 +P4SQS +P46Q6 = <PEl
PSIQI +PS2Q2 +PS3Q3 +PS4Q4 +PssQs +PS6Q6 = <PEl
P6lQI +P62Q2 +P63Q3 +P64Q4 +P6SQS +P66Q6 = <PEl
216 10 Numerische Feldberechnung
Bei n angenommenen Ersatzladungen erstellt man n Gleichungen fUr n Konturpunkte.
In Matrixschreibweise
[ ~11 Pnl
~lnl {~l} = {<P~I} Pnn Qn <PEl
bzw.
(10-50)
Die Auflosung des Gleichungssystems liefert die gesuchten Werte der Ladungen Ql ... Qn,
(10-51)
Damit ist die erste NaherungslOsung gefunden, da sich mit den gefundenen Ladungen Ql ... Q6 das Potential <P(r) an jedem Ort des Feldraums durch Superposition naherungsweise berechnen laBt (Exakte Losungen erhillt man nur fur die anfanglich gewahlten Konturpunkte der Elektrode).
Eine Abschatzung der Genauigkeit der NaherungslOsung ist durch die Berechnung der Potentiale weiterer Punkte auf der Elektrodenkontur, sog. Kontrollpunkte, moglich. Die in den Kontrollpunkten berechneten Potentiale weich en in der Regel yom Sollwert <PE mehr oder weniger ab und erlauben damit eine Beurteilung der ZweckmaBigkeit der anfanglich gewahlten Anzahl, Art und Positionierung der Ersatzladungen. In einem zweiten Durchlauf mit mehr Ladungen, Ersatz diskreter Punktladungen durch Linien- und Ringladungen sowie veranderter Positionierung wird die NaherungslOsung verbessert. pie Anzahl der erforderlichen Wiederholungen des Verfahrens hangt von der Erfahrung des Anwenders abo
Zur Wahrung der Ladungsbilanz I.Q+ + I.Q_ = 0 wird beim Ersatzladungsverfahren meist mit einer Spiegelebene gearbeitet. Damit ergeben sich doppelt so viele Potentialkoeffizienten, die aber jeweils paarweise zusammengefaBt werden konnen, so daB die Ordnung der Matrix erhalten bleibt. Die
10.3 Ersatzladungsverfahren 217
Matrix [P] ist voll besetzt, aber vergleichsweise klein, da sie nur von der Anzahl der Ersatzladungen abhangt. Die Berechnung offener Feldgebiete, z.B. die in Bild 10.1 gezeigte Funkenstrecke oder Freileitungsisolatoren, bereiten beim Ersatzladungsverfahren keine Schwierigkeiten (Vergleiche dagegen Finite-Elemente- oder Finite-Differenzen-Verfahren).
Zweidimensionale gestreckte Probleme legen die Verwendung gerader Linienladungsdichten, dreidimensionale rotationssymmetrische Probleme die Verwendung ringfOrmiger Linienladungsdichten (Ringladungen) nahe (vgl. 2.3). Die bei Verwendung von Ladungsdichten erforderliche Integration kann analytisch oder numerisch erfolgen, erstere ist zur Erzielung kurzer Rechenzeiten vorzuziehen. Flir nichtrotationssymmetrische Probleme mlissen die Ringladungen segmentiert oder problemspezifisch stetig variiert werden. Geschichtete Dielektrika lassen sich durch zusatzliche beidseitige Ladungsbelegungen an den Grenzflachen berlicksichtigen.
Das Ersatzladungsverfahren hat wegen seiner didaktischen Transparenz, und seiner einfachen Implementierbarkeit in der Hochspannungstechnik langere Zeit einen beachtlichen Stellenwert besessen, wird aber derzeit zunehmend durch die in der Mechanik schon lange bewahrte Boundary-Element-Methode, die im folgenden beschrieben wird, abgelOst.
10.4 Boundary-Element-Methode
Aufbauend auf dem Verstandnis des Ersatzladungsverfahrens wird im folgenden kurz die Boundary-Element-Methode vorgestellt, die gegenliber dem Finite-Elemente Verfahren und dem Differenzenverfahren den groBen Vorzug bietet, daB nur die Rander (Elektrodenkonturen etc.) diskretisiert werden, nicht der ganze Raum zwischen ihnen.
Zunachst werden die Elektrodenoberflachen in eine Vielzahl von Oberflachenelementen (engl.: boundary elements, patches) diskretisiert, deren Form, GroBe und Ansatz sich nach den geometrischen Gegebenheiten der Aufgabenstellung richten. Ahnlich wie beim Ersatzladungsverfahren wird dann das Randwertproblem durch Annahme der Existenz individueller Flachenladungsdichten p auf den patches in ein Newton-Problem umgewandelt. Jedes Oberflachenelement i bewirkt dann an einem Ort r im Feldraum einen Potentialbeitrag.
218 10 Numerische Feldberechnung
(10-52)
Unter Beriicksichtigung cler Beitdige der Ladungsdichten Pi aller N Oberfliichenelemente erhalt man dann fur das Potential am Ort r
(10-53)
Schreibt man ahnHch wie beim Ersatzladungsverfahren Gleichung (10-53) fur die bekannten Potentiale CPi(ri) = CPEl in der Mitte eines jeden "boundary elements" an und fuhrt zur Vereinfachung der Schreibweise wieder Potentialkoeffizienten ein, erhalt man ein lineares Gleichungssystem
[P][Q] = {CPEl} • (10-54)
Man lOst nach den unbekannten Ladungen auf
(10-55)
und kann dann wie beim Ersatzladungsverfahren cp(r) an jedem Ort des Raumes im Rahmen des "postprocessing" der Ergebnisse ermitteln.
Der wesentliche Unterschied zum Ersatzladungsverfahren besteht darin, daB bei der Boundary-Element-Methode die vorgegebene Geometrie diskretisiert werden kann und die Raumladungsdichten direkt auf dem Rand sitzen dUrfen. Dariiberhinaus wird die Ladungsdichte auf einem OberfUichenelement nicht konstant angenommen sondern Ublicherweise durch eine Approximationsfunktion angenahert. Beim Ersatzladungsverfahren muBten wir die Ersatzladungen auf einer eigens erzeugten geometrischen Raumflache bzw. Kurve hinter dem Rand anordnen, urn zu vermeiden, daB das Elektrodenpotential am Ort einer auf der Oberflache sitzenden Punktladung unendlich groB wUrde. Dieses Problem existiert bei der BEM grundsatzlich auch, laBt sich jedoch durch geeignete singulare Integrationsroutinen lOsen. Ferner erreicht man durch die Verwendung von Ladungsdichten gegenUber Punktladungen den Vorzug, daB sowohl die Ladungsverteilung Uber der Elektrode
10.4 Boundary-Element-Methode 219
als auch das Potentiallangs der Elektrodenkontur einen wesentlich glatteren Verlauf besitzt und damit auch eine hohere Genauigkeit ermoglicht.
10.5 Momenten - Methode
Die Momentenmethode in der hier vorgestellten Form zeigt eine weitere Moglichkeit auf, zeitlich veranderliche, das heiEt, elektrodynamische Felder zu berechnen. Wie bei der Finite-Elemente-Methode fur nichtstationare Felder (s. SchluE von 10.1) berechnen wir auch hier zunachst das elektrische Skalarpotential q>(r) und das magnetische Vektorpotential A(r), wobei wir die ganze Rechnung im Frequenzbereich vornehmen wollen, das heiEt, wir berechnen die komplexen Amplituden ~(r) und A(r) (s. AS). Die komplexen Amplituden der Feldstarken ergeben sich dann zu
~ = -grad q> - jroA und 1
H = -rotA - J.l -
(10-56)
Die Potentialfunktionen q>(r) und A(r) erhalten wir unter Annahme der Lorentz-Eichung durch L6sen der Helmholtz-Gleichungen (s. 6.3.2),
und
p(r) L'1q>(r) + kgq>(r) =-~
- - {O (10-57)
(10-58)
Ersetzen wir noch die Ladungsdichte p(r) gemaE dem Kontinuitatsgesetz in komplexer Schreibweise (s. 3.4, 3.5 und AS),
(10-59)
durch die Stromdichte !(r), erhalten wir fur Gleichung (10-59)
(10-60)
220 10 Numerische Feldberechnung
Die Losungen der Helmholtz-Gleichungen (l0-57) und (10-58) sind bekanntlich die retardierten Potentiale, die im Frequenzbereich folgende Form annehmen (s. 6.3.2):
(10-61)
und
(10-62)
Die mit q indizierten Variablen beziehen sich dabei auf den Ort der Anregungsfunktion J(rq ).
Ahnlich wie beim Ersatzladungsverfahren die GroBe der Ersatzladungen, und beim Boundary-Element-Verfahren die GroBe der Ladungsdichten aus den Randbedingungen errechnet wurden, wird hier die unbekannte Stromdichte !(rq ) aus den Randbedingungen errechnet. Beschdinkt man sich auf dunne Leiter, kann man statt der Stromdichte mit den Leiterstromen rechnen, weswegen das Verfahren dann auch Ersatzstromverfahren genannt wird.
Wir veranschaulichen das weitere Vorgehen an dunnen stromfiihrenden Leitern, beispielsweise Drahtantennen oder Zuleitungen zu StoBspannungsmeBkreisen der Hochspannungspruftechnik, Bild 10.7.
Bild 10.7.: Ermittlung des durch Strahlungskopplung erzeugten Anteils des MeBsignals eines StoBspannungsteilers.
10.5 Momenten-Methode 221
Die Zuleitung wird in inkrementale Drahtelemente As diskretisiert, die Strome in den Elementen werden durch Approximationsfunktionen I(sq) angenahert (vgl. patches der Boundary-ELement-Methode). Die Vereinfachung des Rechnens mit StroIl1-en verdanken wir dem linienfOrmigen Verlauf des Strompfades, fUr den die Volumenintegrale (10-61) und (10-62) in einfache Linienintegrale iibergehen und div!(rq} sich auf d!(sq)/ dSq reduzieren laRt.
Mit dieser Vereinbarung berechnen sich die Potentiale zu
(10-63)
und
(10-64)
Die elektrische Feldstarke in Drahtrichtung berechnet sich dann mit
~(r) = -grad~(r} - jroA(r) (10-65)
zu
E(s) = - ---- e ds 1 lLq 1 a!(Sq) a (e-ikolr-rql) - 41tE 0 jro aSq as Ir - rql s q
(10-66)
und die Spannung langs eines inkrementalen Drahtelements k der Lange As naherungsweise zu
(10-67)
Nahert man auch noch das Integral in Gl. (10-66) durch eine Summe an, erhalt man naherungsweise
222 10 Numerische Feldberechnung
n
!Ik = L ?;k)i = ?;k1 II + ?;k2 !z + .... + ?;kn In i=l
(10-68)
worin 11' 12 , 13 ,,, den I(sq) der einzelnen Strompfadelemente und die restlichen Faktoren zusammengefaBt den Impedanzen Zkl' Zk2 , Zk3 ... entsprechen. SinngemaB laBt sich fur jedes Element eine ahnliche Gleichung anschreiben,
!11 f;11 h + ?;12 h + ... + ?;In In
!12 ?;21 II + ?;22 h + ... + ?;2n In
!13 = ~1 h + ?;32 h + ... + ?;3n In
!In = ?;n1 h + ?;n2 h + ... + ?;nn In (10-69)
Als Matrizengleichung geschrieben,
?;ll ?;12 ?;13 - - - ?;In II !11 ?;21 ?;22 ?;23 - - - ?;2n h !I2 ?;31 ?;32 f;33 - ?;3n h 1IJ
- - - = - - -- - -
?;n1 ?;n2 ?;n3 ?;nn In !In (10-70)
bzw. [?;][I] = [!1] (10-71)
Da bei einem ideal leitenden Leiter die Spannung langs eines Elements den Wert Null besitzt, nehmen im vorliegenden Beispiel aIle Elemente des Spannungsvektors (Randbedingungen) den Wert Null an bis auf das Element, das die Speisespannungsquelle (Testfunktionsgenerator) reprasentiert. Das Gleichungssystem laBt sich damit nach den unbekannten Stromen der einzelnen
10.5 Momenten-Methode 223
Drahtelemente auflosen, woraus dann mit (10-61) und (10-62) die Potentiale ~ und A und schlieBlich die Felder ;g und H folgen.
Flachenhafte Gebilde, z.~. Schirmgehause, werden entweder durch eine endliche Anzahl zu einander senkrechter linienformiger Leiter nachgebildet oder durch Flachenelemente mit den auf Ihnen herrschenden Flachenstromdichten (Strombelage).
SchlieBlich konnen die im Spannungsvektor auftretenden Randbedingungen noch von externen Feldern herriihren, z.B. bei der Behandlung von Streuproblemen (engl.: scattering). Eine auf ein Leitergebilde einfallende elektromagnetische Welle EE und HE influenziert und induziert dort Quellenspannungen- und Quellenstrome, so daB die Randbedingungen aller Elemente im Spannungsvektor von Null verschieden sind. Die influenzierten und induzierten Quellen sind ihrerseits wieder Ursache einer reflektierten elektromagnetischen Welle ER und HR, die sich mit der einfallenden Welle im gesamten Raum zum gesuchten Nettofeld tiberlagert.
Bei eingehender Befassung mit den verschiedenen in den vorangegangenen Abschnitten beschriebenen Algorithmen und intensivem, systematischen Studium des umfangreichen Schrifttums stellt man fest, daB sich aus Sicht der linearen Analysis bzw. Funktionsanalysis aUe vorgestellten Verfahren generisch als Momentenverfahren darstellen lassen. Jedes Verfahren laBt sich formal auf ein Integral yom Typ
J xn f(x)dx (10-72)
zurtickfUhren, das in der Mathematik als n-tes Moment der Funktion f(x) bezeichnet wird. (Beispielsweise liefert das Integral fUr n=1 das mechanische Moment einer Masseverteilung f(x». Diese Zusammenhange lagen wiihrend der historischen Entwicklung der einzelnen Verfahren nicht offen zu Tage und haben sich erst tiber mehrere Jahrzehnte hinweg herauskristallisiert. Dies erklart die Vielfalt der unterschiedlichen intuitiven Vorgehensweisen und die oft nicht eindeutige Namensgebung der Verfahren.
Selbstverstandlich laBt sich die hier in Anlehnung an das Ersatzladungsverfahren vorgest~llte "Momenten-Methode" auch unter Verwendung der im Abschnitt 10.1 "Finite-Elemente"-Verfahren eingefUhrten "Variationsrechnung" und auch der "Gewichteten Residuen" herleiten. Der Anschaulichkeit
224 10 Numerische Feldberechnung
halber wurde jedoch dem "Ersatzstromverfahren" der Vorzug gegeben. Die Alternativen findet der Leser im umfangreichen Spezialschrifttum.
10.6 Monte - Carlo - Methode
Die Monte - Carlo - Methode beruht auf dem Mittelwertsatz der Potentialtheorie, der das Potential im Zentrum einer Kugel als Mittelwert der Potentiale auf der KugeloberfUiche darstellt, Bild 10.B.
Bild 10.8: Mittelwertsatz der Potentialtheorie.
Unterteilt man die KugeloberfUiche in n kleine TeilfUichen IlA gleicher GroEe - mit anderen Worten 4n:r02 = nM-, liiEt sich das Integral durch eine Summe anniihern,
bzw.
1 n
<p(r)=-L<PjM nM.
1=1
1 n <p(r) = - ~ <p.
n~ 1 j=1
(10-73)
(10-74)
10.6 Monte-Carlo-Methode 225
Die Monte-Carlo-Methode der Potentialberechnung stellt nun eine statistische Interpretation dieser Beziehung dar. Beispielsweise uberzieht man bei einem 2-dimensionalen Problem das Feldgebiet mit einem aquidistanten Punktgitter, Bild 10.9.
<P-100%
~""""";:""""""""",,,,~ 1 --- 3
<I>=?
2
~/////////////////-)
<p=0%
Bild 10.9: Zur Veranschaulichung der Monte-Carlo-Methode der Potentialfeldberechnung. Das Punktgitter hat nicht die Aufgabe der Diskretisierung des Feldraums wie beim Differenzenverfahren. Es dient hier lediglich der Festlegung der Wege beim Fortbewegen vom Startpunkt.
AnschlieBend bewegt man sich ausgehend von einem Feldpunkt, des sen Potential gesucht ist, zu einem moglichen Nachbarpunkt weiter, wobei die Richtung in jedem Gitterpunkt erneut von einem Zufallsgenerator mit den Zahlen 1, 2, 3,4 festgelegt wird. Bei N Versuchen gelangt man NlOOoffimal zum Potential <1>100% und Nooffimal zum Potential <1>0%' Fur eine angemessene Zahl von N = NlOO% + NO% Versuchen berechnet sich dann das Potential des Startpunkts zu
(10-75)
Der Vorteil der Monte-Carlo-Methode liegt in ihrer Einfachheit und der wenig aufwendigen Programmiertechnik, ihr Nachteil im groBen Rechenzeitbedarf, da ja nach N Versuchen jeweils nur das Potential eines einzigen Feldpunkts erhalten wird. Eine Steigerung ihrer Effizienz erlaubt die
226 10 Numerische Feldberechnung
Floating-Random- Walk-Methode mit veranderlicher Schrittweite, auf die hier jedoch nicht weiter eingegangen werden soll. Bislang hat die MonteCarlo-Methode in der numerischen Feldberechnung nur vergleichsweise geringe Verbreitung gefunden.
10.7 Allgemeine Bemerkungen zur numerischen Feldberechnung
Grundsatzlich gliedert sich die numerische Feldberechnung in drei Blocke, Bild 10.10.
Eingabe Feldberechnung Ausgabe
Geometriedaten ELV Verfahrenspez. FLV Aufstellen <j)v' Ey, Bv Daten BEM und lOsen
FDM eines Linearen
FEM Gleichungs- Feldlinien, Eingabe- MOM systems Emax, Datensatz MCM Optimieren etc.
Bild 10.10: Grobstruktur der numerischen Feldberechnung. ELV: Ersatzladungsverfahren; FLV: Flachenladungsverfahren; BEM: Boundary-Element-Methode; FDM: Finite-Differenzen-Methode; FEM: Finite-Elemente-Methode; MOM: Momenten -Methode; M CM: Monte-Carlo-Methode.
Abhangig yom Verfahren sind die Schnittstellen variabel, d.h. Teilaufgaben konnen auch dem jeweils benachbarten Block zugeordnet sein,
- Block 1 (engl.: pre-processing) beinhaltet die effiziente Aufbereitung der physikalischen Aufgabenstellung in eine fur den Rechner akzeptable Form, m.a.W. die "anwenderfreundliche" Erstellung des Datensatzes (Problemdefinition).
10.7 Allgemeine Bemerkungen zur numerischen Feldberechnung 227
- Block 2 (eng!.: processing) befaBt sich mit dem Erstellen eines linearen Gleichungssystems, bzw. der verfahrensspezifischen Berechnung der Elemente seiner Koeffizientenmatrix und dessen effizienter Losung.
- Block 3 (eng!.: post-processing) lei stet die anwenderfreundliche Auswertung und Darstellung des Rechenergebnisses in Form von Aquipotential- und FeldsUirkelinien, der Darstellung von Potential- und FeldsHirkeverteilungen langs vorgegebener Raumkurven usw.
Von der gesamten Problematik wurde im vorangegangenen nur der mittlere Block einfUhrend behandelt. Allein dies em Block ware bei weiterer Befassung noch viel hinzuzufiigen. Beispielsweise, wie sich bei den verschiedenen Verfahren Raumladungen, Elektroden auf freiem Potential und verschiedene Dielektrika beriicksichtigen lassen, wie Speicherplatzbedarf verringert, Rechengeschwindigkeit erhoht oder Rechengenauigkeit verbessert werden konnen. Ferner, welche Moglichkeiten zur Abschatzung der Genauigkeit der NaherungslOsungen bestehen, wann und warum bestimmte Problemstellung en zu numerischen Schwierigkeiten fiihren. SchlieBlich konnen verschiedene Verfahren kombiniert werden, beispielsweise in dem man sich in einem abgeschlossenen Gebiet der Vorteile der Finite-Elemente-Methode, auBerhalb der Vorteile der Boundary-Element-Methode bedient.
Nicht weniger komplex sind der erste und dritte Block fur die Kommunikation zwischen Anwender und mittlerem Block (Benutzerschnittstelle). Wegen der Vielfalt der hier in Frage kommenden Moglichkeiten muB auf das fachspezifische Schrifttum bzw. auf die Benutzeranleitungen der verschiedenen Programmpakete verwiesen werden.
In der numerischen Feldberechnung ist noch vieles im FluB, kiinftige Entwicklungen werden insbesondere einer Verbesserung der Ein- und Ausgabeprozeduren gewidmet sein sowie Methoden der Kiinstlichen Intelligenz (KI) beriicksichtigen.