Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… ·...

24
Inhaltsverzeichnis 1 Einf¨ uhrung ...................................................... 1 1.1 Allgemeine Bemerkungen ....................................... 1 1.2 Einleitende Beispiele ........................................... 3 1.3 ¨ Uberblick ..................................................... 10 2 Modellierung und Simulationen mit finiten Differenzenverfahren . 13 2.1 Modellgleichungen elektrostatischer Probleme ..................... 14 2.2 Das eindimensionale elektrostatische Problem ..................... 16 2.3 Das zweidimensionale elektrostatische Problem .................... 20 2.3.1 Diskrete Beschreibung der Geometrie ...................... 20 2.3.2 Ersetzen der partiellen Ableitungen durch finite Differenzen . . . 22 2.3.3 Aufstellen des linearen Gleichungssystems .................. 24 2.3.4 osen des linearen Gleichungssystems durch geeignete Methoden 25 2.3.5 Lineare Interpolation .................................... 28 2.4 Verallgemeinerung ............................................. 31 2.5 Aufgaben zur finiten Differenzenmethode ......................... 33 3 Randangepasste Gitter .......................................... 35 3.1 Beschreibung anwendungsrelevanter Gebiete ...................... 35 3.2 Erzeugung von randangepassten Gittern .......................... 37 3.3 osen der Poisson-Gleichung auf randangepassten Gittern .......... 43 3.4 Aufgaben zu randangepassten Gittern ............................ 49 4 Finite-Elemente-Methode f¨ ur eindimensionale Probleme ......... 51 4.1 Variationsproblem statt Differenzialgleichung ..................... 51 4.2 Minimierung des Energiefunktionals ............................. 55 4.3 Beispiele ..................................................... 59 4.4 Aufgaben zur Finiten-Elemente-Methode (1D) .................... 64 5 Finite-Elemente-Methode bei elliptischen Randwertproblemen ... 67 5.1 Triangulierung mit linearen Basisfunktionen ...................... 70 5.2 Visualisierung der Finiten-Elemente-Methode ..................... 76 5.3 Triangulierung mit linearen Elementfunktionen .................... 79 5.4 Rechteckzerlegung mit bilinearen Elementfunktionen ............... 82 5.5 Triangulierung mit quadratischen Elementfunktionen .............. 85 5.6 Aufgaben zur Finiten-Elemente-Methode (2D) .................... 94

Transcript of Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… ·...

Page 1: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

Inhaltsverzeichnis

1 Einfuhrung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Allgemeine Bemerkungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Einleitende Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Uberblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Modellierung und Simulationen mit finiten Differenzenverfahren . 132.1 Modellgleichungen elektrostatischer Probleme . . . . . . . . . . . . . . . . . . . . . 142.2 Das eindimensionale elektrostatische Problem . . . . . . . . . . . . . . . . . . . . . 162.3 Das zweidimensionale elektrostatische Problem . . . . . . . . . . . . . . . . . . . . 20

2.3.1 Diskrete Beschreibung der Geometrie . . . . . . . . . . . . . . . . . . . . . . 202.3.2 Ersetzen der partiellen Ableitungen durch finite Differenzen . . . 222.3.3 Aufstellen des linearen Gleichungssystems . . . . . . . . . . . . . . . . . . 242.3.4 Losen des linearen Gleichungssystems durch geeignete Methoden 252.3.5 Lineare Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4 Verallgemeinerung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.5 Aufgaben zur finiten Differenzenmethode . . . . . . . . . . . . . . . . . . . . . . . . . 33

3 Randangepasste Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.1 Beschreibung anwendungsrelevanter Gebiete . . . . . . . . . . . . . . . . . . . . . . 353.2 Erzeugung von randangepassten Gittern . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3 Losen der Poisson-Gleichung auf randangepassten Gittern . . . . . . . . . . 433.4 Aufgaben zu randangepassten Gittern . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4 Finite-Elemente-Methode fur eindimensionale Probleme . . . . . . . . . 514.1 Variationsproblem statt Differenzialgleichung . . . . . . . . . . . . . . . . . . . . . 514.2 Minimierung des Energiefunktionals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.3 Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.4 Aufgaben zur Finiten-Elemente-Methode (1D) . . . . . . . . . . . . . . . . . . . . 64

5 Finite-Elemente-Methode bei elliptischen Randwertproblemen . . . 675.1 Triangulierung mit linearen Basisfunktionen . . . . . . . . . . . . . . . . . . . . . . 705.2 Visualisierung der Finiten-Elemente-Methode . . . . . . . . . . . . . . . . . . . . . 765.3 Triangulierung mit linearen Elementfunktionen . . . . . . . . . . . . . . . . . . . . 795.4 Rechteckzerlegung mit bilinearen Elementfunktionen . . . . . . . . . . . . . . . 825.5 Triangulierung mit quadratischen Elementfunktionen . . . . . . . . . . . . . . 855.6 Aufgaben zur Finiten-Elemente-Methode (2D) . . . . . . . . . . . . . . . . . . . . 94

Page 2: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

viii Inhaltsverzeichnis

6 Einfuhrung in ANSYS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 976.1 Die Benutzeroberflache von ANSYS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 986.2 Elektrostatische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.3 Thermische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1046.4 Mechanische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1136.5 Magnetische Simulation: Stromdurchflossener Leiter . . . . . . . . . . . . . . . 1226.6 Aufgaben zu ANSYS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

7 ANSYS-Simulationen - Projektarbeiten . . . . . . . . . . . . . . . . . . . . . . . . . . . 1357.1 Kraftebestimmung bei Schraubschlussel und Schrauben . . . . . . . . . . . . 1367.2 Modalanalyse eines Ultraschallgebers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1377.3 Kapazitives System zur Fullstandsmessung bei Hubschraubern . . . . . . 1407.4 Simulation eines Beschleunigungsmess-Systems . . . . . . . . . . . . . . . . . . . . 1437.5 Optimierung des Temperaturprofils eines SnO2-Sensors . . . . . . . . . . . . 1467.6 Optimierung einer Fingerspule fur die Kernspintomographie . . . . . . . . 1487.7 Magnetfeldberechnung bei Planarspulen . . . . . . . . . . . . . . . . . . . . . . . . . . 1517.8 Ausbreitung elektromagnetischer Strahlung . . . . . . . . . . . . . . . . . . . . . . . 153

Anhang

A Losen von großen linearen Gleichungssystemen . . . . . . . . . . . . . . . . . . . 157A.1 Direkte Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158A.2 Klassische iterative Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163A.3 Das Verfahren der konjugierten Gradienten . . . . . . . . . . . . . . . . . . . . . . . 172A.4 Aufgaben zum Losen von großen LGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177

B Numerisches Differenzieren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179B.1 Differenzenformeln fur die erste Ableitung . . . . . . . . . . . . . . . . . . . . . . . . 179B.2 Differenzenformeln fur die zweite Ableitung . . . . . . . . . . . . . . . . . . . . . . . 184B.3 Differenzenformeln fur die n-te Ableitung . . . . . . . . . . . . . . . . . . . . . . . . . 185B.4 Aufgaben zum numerischen Differenzieren . . . . . . . . . . . . . . . . . . . . . . . . 186

C Logfiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187C.1 Elektrostatische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187C.2 Statische, thermische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188C.3 Transiente, thermische Simulation: Ein-Last-Simulation . . . . . . . . . . . . 189C.4 Transiente, thermische Simulation: Mehr-Lasten-Simulation . . . . . . . . . 190C.5 Mechanische Simulation (statisch) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191C.6 Mechanische Simulation (Modalanalyse) . . . . . . . . . . . . . . . . . . . . . . . . . . 193C.7 Gleichstrom Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195C.8 Wechselstrom Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196

Literaturverzeichnis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199

Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201

ANSYS-Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203

Homepage zum Buch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205

Page 3: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.Modellierung und Simulationen mit finitenDifferenzenverfahren

Wir werden in diesem Kapitel die Methode der finiten Differenzen einfuhren, um be-liebige partielle Differenzialgleichungen bezuglich den Ortskoordinaten numerisch zulosen. Wir beschranken uns bei der Beschreibung der Methoden auf zweidimensio-nale elektrostatische Probleme, da man fur diese Probleme die Moglichkeit hat, siemesstechnisch zu erfassen. Man ist somit in der Lage, die numerischen Losungen mitexperimentellen Messungen zu vergleichen. Die Methode der finiten Differenzen istaber nicht auf die zweidimensionalen elektrostatischen Probleme beschrankt, sondernsie ist auf den dreidimensionalen Fall und auf andere partielle Differenzialgleichungendirekt ubertragbar.

Das von einem System geladener Elektroden erzeugte elektrische Feld ist qualitativdas gleiche, unabhangig davon, ob es sich in Luft oder in einem Elektrolyten befin-det. In beiden Fallen treten keine feldverzerrenden Raumladungen auf. Diese Tatsacheliefert eine bequeme Methode zur experimentellen Ermittlung des Verlaufs solcher Fel-der: der elektrostatische Trog : Man kopiert das zu untersuchende Elektrodensystem,bringt es in einen Elektrolyten, erzeugt an ihm dem Original entsprechende Span-nungsverhaltnisse und misst die Spannungen im Innern z.B. mit einer Drahtsonde.

Genau genommen ermittelt man Aquipotenziallinien, d.h. Linien mit gleichem elek-trischen Potenzial. Die elektrischen Feldlinien stehen auf diesen Linien senkrecht; daselektrische Feld ist durch den lokalen Gradienten der Potenzialverteilung gegeben. Ma-thematisch formuliert ist das Potenzial Φ eine Funktion von den zwei Ortsvariablenx und y, Φ(x, y), und die elektrische Feldstarke

−→E = −grad(Φ) = −

(∂xΦ(x, y)∂yΦ(x, y)

).

In Abbildung 2.1 sind die Potenzialverlaufe fur ein Drei-Elektroden-System mit Spalt-blende angegeben. Beide Abbildungen sind experimentell ermittelte Daten, die imRahmen eines physikalischen Praktikums an der Hochschule Karlsruhe ermittelt wur-den. In Abb. 2.1 (b) mit und in Abb. 2.1 (a) ohne Durchgriff. Die Potenzialwerte derAquipotenziallinien sind die im Bild unten angegebenen Werte.

Die linke Elektrode liegt in beiden Fallen auf 10 V , die Zwischenelektrode auf 0 V .Der einzige Unterschied zwischen den beiden Konfigurationen besteht im Potenzial-wert der rechten Elektrode: In der linken Anordnung betragt das Potenzial 4 V ; inder rechten 2 V .

Page 4: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

14 2 Modellierung und Simulationen mit finiten Differenzenverfahren

(a) mit Separatrix (b) mit Durchgriff

Abb. 2.1. Potenzialverteilung im elektrostatischen Trog.

Vergleicht man die beiden Abbildungen, erkennt man sehr gut, dass die Potenzial-verteilung im Innern des Gebietes nicht nur quantitativ, sondern auch qualitativ un-terschiedlich ist, je nachdem welchen Wert die Spannung auf der rechten Elektrodehat. Dies ist charakteristisch fur alle elektrostatischen Probleme: Die Rander undRandbedingungen (sowohl die Form als auch die Werte) bestimmen denVerlauf der Losung im Innern des Gebiets. Es mussen folglich auch bei derSimulation immer die Randbedingungen vollstandig spezifiziert werden.

Im Folgenden werden wir numerische Verfahren kennenlernen, um solche elektrostati-schen Konfigurationen zu berechnen. Damit werden wir dann in der Lage sein, nume-rische Experimente durchzufuhren und die Potenzialverteilung bzw. das elektrischeFeld zu gegebener Anordnung zu bestimmen.

Zunachst stellen wir die Modellgleichungen fur elektrostatische Probleme auf. DieseModellgleichungen sind fur alle elektrostatischen Probleme universell gultig. Aller-dings lassen sie sich in der Regel nur naherungsweise losen. Wir fuhren die finiteDifferenzenmethode fur das eindimensionale Problem ein und vergleichen fur diesesProblem die numerische Losung mit der exakten. Anschließend verallgemeinern wirdie Methode auf den zweidimensionalen Fall.

2.1 Modellgleichungen elektrostatischer Probleme

Die Beschreibung elektrostatischer Probleme erfolgt durch drei physikalische Großen:die Ladungsdichte ρ, das elektrostatische Potenzial Φ und das elektrische Feld

−→E . Da

wir uns im Folgenden auf eine zweidimensionale Diskussion beschranken, hangen diese

Page 5: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.1 Modellgleichungen elektrostatischer Probleme 15

Großen dann nur von den beiden Ortsvariablen x und y ab. Die drei Beschreibungs-großen werden durch zwei universell geltende Gleichungen verknupft:

→ Das von einem System geladener Elektroden erzeugte elektrische Feld wird durchden Gradienten des elektrostatischen Potenzials Φ(x, y) beschrieben:

−→E (x, y) = − grad Φ(x, y) = −

(∂xΦ(x, y)∂yΦ(x, y)

). (2.1)

Das elektrische Feld ist also der negative Gradient des Potenzials.

→ Der Gaußsche Satz verknupft das elektrische Feld mit der im System vorhandenenLadungsdichte ρ(x, y)

∇−→E (x, y) =

1ερ(x, y). (2.2)

Dieser Satz besagt, dass die Quellen des elektrischen Feldes die Ladungsdichtensind. ε ist dabei die Dielektrizitatskonstante. Fur Vakuum gilt ε = ε0 = 8.8·1012 F

m .

Setzt man Gleichung (2.1) in (2.2) ein

1ερ(x, y) = ∇

−→E (x, y) = ∂xE1(x, y) + ∂yE2(x, y)

= ∂x(−∂xΦ(x, y)) + ∂y(−∂yΦ(x, y))= −∂2

xΦ(x, y)− ∂2yΦ(x, y)

erhalt man die Poisson-Gleichung

Φxx(x, y) + Φyy(x, y) = −1ερ(x, y). (2.3)

Fur den ladungsfreien Raum (ρ = 0) gilt die Laplace-Gleichung

Φxx(x, y) + Φyy(x, y) = 0. (2.4)

Oftmals kurzt man die linke Seite der Differenzialgleichung durch den Laplace-Operator 4Φ(x, y) := Φxx(x, y) + Φyy(x, y) ab, so dass man die Laplace-Gleichung inder Form 4Φ(x, y) = 0 schreibt.

Die gute Nachricht uber die Poisson- bzw. Laplace-Gleichung ist, dass alle elektro-statischen Probleme durch sie beschrieben werden. Die schlechte Nachricht dabei istallerdings, dass man sie in der Regel fur die wenigsten technischen Probleme mathe-matisch exakt losen kann. Jedes Programm, das in der Lage ist diese Differenzial-gleichung naherungsweise zu losen, kann daher beliebige elektrostatische Problemenaherungsweise bestimmen.

Page 6: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

16 2 Modellierung und Simulationen mit finiten Differenzenverfahren

Zusammenfassung: Elektrostatische Modellgleichungen.

Jedes zweidimensionale elektrostatische Problem lasst sich bei gegebener La-dungsdichte ρ durch die Poisson-Gleichung

Φxx(x, y) + Φyy(x, y) = −1ερ(x, y)

beschreiben. Fur den ladungsfreien Raum gilt die Laplace-Gleichung

Φxx(x, y) + Φyy(x, y) = 0.

Im Folgenden werden wir die finite Differenzenmethode anwenden, um eine Nahe-rungslosung fur diese partielle Differenzialgleichung auf achsenparallelen Gittern zubestimmen. Bevor wir jedoch das Problem im Zweidimensionalen losen, werden wirdie Differenzenmethode fur den eindimensionalen Fall einfuhren:

2.2 Das eindimensionale elektrostatische Problem

Gesucht ist die Potenzialverteilung in einem ebenen Plattenkondensator siehe Ab-bildung 2.2 mit Spaltabstand d, Kathodenpotenzial ΦK und Anodenpotenzial ΦA.

Abb. 2.2. Ebener Plattenkondensator.

Fur den ebenen Plattenkondensator ist das Potenzial nur eine Funktion von x, d.h.Φ = Φ(x), so dass die Laplace-Gleichung (2.4) sich reduziert zu

Φ′′(x) = 0.

Durch zweimalige Integration erhalten wir in diesem einfachsten Fall den exaktenVerlauf des Potenzials

Φ(x) = ΦK(x

d) + ΦA(1− x

d),

wobei die Randbedingungen Φ(0) = ΦA und Φ(d) = ΦK berucksichtigt sind. Das Po-tenzial nimmt linear vom Anoden- zum Kathodenpotenzial ab.

Page 7: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.2 Das eindimensionale elektrostatische Problem 17

Wir werden dieses eindimensionale Problem nun numerisch losen. Im Gegensatz zurexakten Losung konnen wir mit dem Rechner das Potenzial nicht an allen Punktenim Raum berechnen. Dies wurde unendlich viel Speicher und Rechenzeit benotigen.Da wir die Potenzialverteilung also nicht kontinuierlich, sondern nur an bestimmten,diskreten Punkten berechnen konnen, fuhren wir diskrete Gitterpunkte auf der x-Achse ein und berechnen die gesuchte Funktion Φ(x) nur an diesen Gitterpunkten.

Abb. 2.3. Diskrete Gitterpunkte.

Fur das obige eindimensionale Beispiel fuhren wir insgesamt 7 Gitterpunkte x0, . . . , x6

ein und bestimmen das unbekannte Potenzial auf den funf Gitterpunkten x1, . . . , x5.Das Potenzial ist am Rand bei x0 = 0 und x6 = d durch die Randbedingung ΦA undΦK vorgegeben. Die gesuchte Funktion Φ ist in jedem inneren Punkt durch die Bedin-gung Φ′′(x) = 0 bestimmt. Denn die zweidimensionale Laplace-Gleichung reduziertsich im Falle von nur einer Variablen x zu dieser gewohnlichen Differenzialgleichung.Also gilt insbesondere an den funf Gitterpunkten x1, . . . , x5

Φ′′(xi) = 0 i = 1, . . . , 5.

Da wir die Funktion nur an den diskreten Gitterpunkten zur Verfugung haben, erset-zen wir die Ableitung Φ′′(xi) durch den zentralen Differenzenquotienten

Φ′′(xi) ∼Φ(xi−1)− 2Φ(xi) + Φ(xi+1)

h2i = 1, . . . , 5

mit h = xi+1 − xi. Mit der Notation Φi = Φ(xi) (vgl. Abbildung 2.4) folgt

Φ′′(xi) ∼Φi−1 − 2Φi + Φi+1

h2.

Die Herleitung der Differenzenformeln fur die Ableitungen auch hoherer Ordnung sindim Anhang B separat in einem eigenstandigen Kapitel beschrieben.

Abb. 2.4. Numerische Approximation.

Page 8: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

18 2 Modellierung und Simulationen mit finiten Differenzenverfahren

Wir ersetzen die kontinuierliche Ableitung der Laplace-Gleichung Φ′′(x) = 0 an jedeminneren Gitterpunkt xi durch die numerische Approximation und erhalten

x1: Φ0 − 2Φ1 + Φ2 = 0x2: Φ1 − 2Φ2 + Φ3 = 0x3: Φ2 − 2Φ3 + Φ4 = 0x4: Φ3 − 2Φ4 + Φ5 = 0x5: Φ4 − 2Φ5 + Φ6 = 0

wobei Φ0 und Φ6 die vorgegebenen Randbedingungen sind: Φ0 = ΦA = 1V , Φ6 =ΦK = 0V .

− 2Φ1 + Φ2 = −1VΦ1 − 2Φ2 + Φ3 = 0

Φ2 − 2Φ3 + Φ4 = 0Φ3 − 2Φ4 + Φ5 = 0

Φ4 − 2Φ5 = 0.

Dies ist ein lineares inhomogenes Gleichungssystem fur die funf Unbekannten Φ1, Φ2,Φ3, Φ4 und Φ5:

−2 1 0 0 01 −2 1 0 00 1 −2 1 00 0 1 −2 10 0 0 1 −2

∣∣∣∣∣∣∣∣∣∣−1

0000

−2 1 0 0 0

0 − 32 1 0 0

0 0 − 43 1 0

0 0 0 − 54 1

0 0 0 0 − 65

∣∣∣∣∣∣∣∣∣∣−1− 1

2− 1

3− 1

4− 1

5

.

Das linke Koeffizientenschema basiert auf einer Tridiagonalmatrix: Nur die Elemen-te der Hauptdiagonalen und den beiden Nebendiagonalen sind von Null verschieden.Durch Losen des Tridiagonalsystems mit dem Thomas-Algorithmus (siehe AnhangA: Losen von großen linearen Gleichungssystemen) erhalt man das rechte Koeffizien-tenschema und durch Ruckwartsauflosen die Losung:

Φ5 =16, Φ4 =

13, Φ3 =

12, Φ2 =

23, Φ1 =

56.

Dies ist die gesuchte Losung an den diskreten Punkten x1, · · · , x5. Fur die graphischeDarstellung verbinden wir die Punkte geradlinig, so dass wir eine stetige, stuckweiselineare Funktion als numerische Losung erhalten:

Abb. 2.5. Graphische Darstellung des Ergebnisses.

Page 9: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.2 Das eindimensionale elektrostatische Problem 19

Diskussion: Im Fall des ebenen Plattenkondensators erhalten wir durch das numeri-sche Verfahren die exakte Losung! Dies liegt daran, dass die Losung des Problems einelineare Funktion ist und das Differenzenverfahren lineare Funktionen exakt differen-ziert. In der Regel wird dies nicht der Fall sein: Denn das numerische Verfahren liefertnur eine Approximation an die exakte Losung! Diese numerische Losung wird umsobesser (genauer) je mehr Gitterpunkte verwendet werden, wenn man Rundungsfehlervernachlassigt. Berucksichtigt man allerdings Rundungsfehler, dann wachst der Ge-samtfehler fur zu kleine Gitterabstande wieder stark an (siehe Anhang B).

Fuhren wir 5 innere Gitterpunkte ein, erhalten wir 5 Unbekannte und eine 5 × 5-Matrix. Wahlen wir 50 oder 100 Gitterpunkte, dann bekommen wir 50 bzw. 100Unbekannte und mussen ein 50× 50- bzw. ein 100× 100-System losen. Die Strukturder Matrix bleibt dabei aber erhalten, da fur jeden Gitterpunkt immer nur die Infor-mation am betrachteten Punkt sowie der direkt benachbarten Punkte einfließt! �

Zusammenfassung: Finite Differenzenmethode.

Der in diesem Abschnitt aufgezeigte Weg ist grundlegend fur alle finite Differen-zenverfahren:

(1) Ersetzen der Geometrie durch diskrete Punkte: Die Geometrie wirddurch die gewahlten Gitterpunkte reprasentiert.

(2) Ersetzen der Ableitungen durch finite Differenzen: Nicht die eigent-liche Modellgleichung wird gelost, sondern die diskretisierte Differenzenglei-chung.

(3) Aufstellen des zugehorigen LGS: Das Gleichungssystem enthalt genau soviele Gleichungen wie Unbekannte (= Anzahl der Gitterpunkte). Fur lineareDifferenzialgleichungen sind die Differenzengleichungen lineare Gleichungs-systeme.

(4) Losen des LGS mit einer geeigneten Methode: In dem diskutiertenBeispiel wird der Thomas-Algorithmus gewahlt.

(5) Lineare Interpolation zwischen den diskreten Punkten: Durch dasLosen der Differenzengleichungen erhalt man die Losung nur an den Git-terpunkten. Man hat keine weitere zusatzliche Information dazwischen. Alseinfachste Methode wahlt man die lineare Verbindung der Losungswerte.

4! Achtung: Die Differenzenmethode besitzt zwei systematische Fehlerquellen:

©1 Die Geometrie wird nicht kontinuierlich beschrieben, sondern durch endlich viele,diskrete Gitterpunkte.

©2 Statt der physikalischen Modellgleichung (Differenzialgleichung) lost man die dis-kretisierte Differenzialgleichung.

Page 10: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

20 2 Modellierung und Simulationen mit finiten Differenzenverfahren

2.3 Das zweidimensionale elektrostatische Problem

Im Folgenden diskutieren wir ein zweidimensionales elektrostatisches Problem amBeispiel des in Abbildung 2.6 angegebenen Drei-Elektroden-Systems mit Spaltblende:

x [cm]

y [cm]

0 5 20

10V2V

6

14

0V

0V

Abb. 2.6. Drei-Elektroden-System mit Spaltblende.

Da in diesem System das Potenzial Φ eine Funktion von x und y ist, muss man ei-ne zweidimensionale Beschreibung wahlen. Wie beim eindimensionalen Fall kann diegesuchte Potenzialverteilung Φ(x, y) nicht als kontinuierliche Funktion im ganzen Be-rechnungsgebiet bestimmt werden, sondern an gewissen (=diskreten) Punkten. Die imAbschnitt 2.2 eingefuhrte Differenzenmethode wenden wir nun auf das zweidimensio-nale Randwertproblem an. Wir folgen der Zusammenfassung fur den eindimensionalenFall und fuhren ein zweidimensionales Berechnungsgitter ein, um Φ(x, y) an den dis-kreten Gitterpunkten zu bestimmen.

2.3.1 Diskrete Beschreibung der Geometrie

Als Rander der Geometrie werden im obigen Problem nicht nur die außeren RanderΦ(x = 0, y) = Φk und Φ(x = d, y) = ΦA berucksichtigt, sondern zusatzlich auch dieinneren Rander Φ(x = 5, 0 ≤ y ≤ 6) = 0V und Φ(x = 5, 14 ≤ y ≤ 20) = 0V . Damitmuss das Gitter so gewahlt werden, dass die Spaltblenden durch Gitterpunkte darge-stellt werden. Die Linien bzw. die Punkte in Abbildung 2.7(a) und Abbildung 2.7(b)dienen als ein solches diskretes Gitter.

Wir definieren die Gitterlinien parallel zu den Koordinatenachsen. Dazu wahlen wireine gleichmaßige Unterteilung sowohl in x- als auch in y-Richtung von 40 Intervallen.Damit erhalt man 41 Gitterpunkte jeweils in x- und 41 Gitterpunkte in y-Richtung(d.h. insgesamt 41× 41 Gitterpunkte).

∆x =d

40=

200mm

40= 5mm; ∆y =

200mm

40= 5mm.

Page 11: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.3 Das zweidimensionale elektrostatische Problem 21

Die Gitterpunkte werden gemaß der Anordnung einer Matrix als Array gespeichert:x[1...41,1...41], y[1...41,1...41]. In unserem Fall sind die Koordinaten der Gitterpunkte

x[i,j] = (i− 1)∆x und y[i,j] = (j − 1)∆y.

Insbesondere wird durch diese Wahl der Gitterpunkte die Spaltblende durch Gitter-punkte reprasentiert!

x [mm]

y [mm]

0,0 50,0 200,0

60,0

140,0

200,0

x [mm]

y [mm]

0,0 50,0 200,0

60,0

140,0

200,0

Abb. 2.7. Gitter fur das Drei-Elektroden-System.

4! Achtung: Ein aquidistantes 42 × 42-Gitter ist fur das Drei-Elektroden-Systemnicht moglich, um die Geometrie zu reprasentieren. Denn in diesem Fall wurden dieinneren Elektroden nicht durch Gitterpunkte erfasst werden; sie waren somit in der Si-mulation nicht vorhanden und man wurde den ebenen Plattenkondensator ohne Spalt-blende simulieren. Entweder muss man auf die Forderung eines aquidistanten Gittersmit konstanten Maschenweiten ∆x, ∆y verzichten und ein achsenparalleles Gittermit unterschiedlichen ∆xi und ∆yj einfuhren. Oder man musste zu einer 81 × 81-Unterteilung ubergehen.

Beispiel 2.1 (Notation). Ein Gitter mit 6 × 5 Gitterpunkten zusammen mit derNotation der Gitterpunkte ist in Abbildung 2.8 angegeben:

(4/3)

Physikalische Gitterpunkte

1/3 2/3 3/3 4/3 5/3 6/3

1/2 2/2 3/2 4/2 5/2 6/2

1/5 2/5 3/5 4/5 5/5 6/5

1/4 2/4 3/4 4/4 5/4 6/4

Nummerierung

1/1 2/1 3/1 4/1 5/1 6/1

Abb. 2.8. Gitterpunkte und Nummerierung.

Page 12: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

22 2 Modellierung und Simulationen mit finiten Differenzenverfahren

Links sind die Gitterpunkte dargestellt, rechts die Nummerierung dieser Gitterpunkte.Die Nummerierung beginnt links unten mit (1/1). Der erste Index gibt die Positionin einer horizontalen Linie (Spaltennummer) und der zweite Index die Liniennummeran.

2.3.2 Ersetzen der partiellen Ableitungen durch finite Differenzen

Nachdem die Geometrie durch ein diskretes Gitter erfasst ist, muss die Modellglei-chung (= Laplace-Gleichung)

Φxx(x, y) + Φyy(x, y) = 0

auf diesen Gitterpunkten gelost werden.

Aber gerade im zweidimensionalen Fall kann die partielle Differenzialgleichung nichtexakt gelost werden. Man muss daher auf Naherungsverfahren zuruckgreifen. Folgenwir der Zusammenfassung des eindimensionalen Problems, ersetzen wir die konti-nuierlichen Ableitungen in der Laplace-Gleichung durch finite Differenzen, um zurdiskreten Laplace-Gleichung zu kommen.

Betrachten wir zunachst einen beliebigen Gitterpunkt (i, j) im Berechnungsgebiet.Dieser Gitterpunkt besitzt die Koordinaten (xij , yij). Da die Laplace-Gleichung injedem inneren Punkt des Gebietes gilt, ist sie auch im Punkt (xij , yij) gultig. Andiesem Punkt ersetzen wir die kontinuierlichen Ableitungen Φxx und Φyy durch finiteDifferenzen:

Φxx(xij , yij) ∼Φ(xi−1,j , yi−1,j)− 2Φ(xij , yij) + Φ(xi+1,j , yi+1,j)

(∆x)2

und

Φyy(xij , yij) ∼Φ(xi,j−1, yi,j−1)− 2Φ(xij , yij) + Φ(xi,j+1, yi,j+1)

(∆y)2,

wobei die Gitterpunkte folgende geometrische Anordnung besitzen.

Abb. 2.9. Geometrische Anordnung des Differenzenoperators.

Wir setzen zur Abkurzung Φij = Φ(xij , yij). Dann ist Φij das Potenzial im Punk-te (xij , yij) und wir konnen die diskrete Laplace-Gleichung an der Stelle (xij , yij)schreiben als

Page 13: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.3 Das zweidimensionale elektrostatische Problem 23

Φxx(xij , yij) + Φyy(xij , yij) ∼Φi−1,j − 2Φij + Φi+1,j

(∆x)2+

Φi,j−1 − 2Φij + Φi,j+1

(∆y)2= 0

bzw.

1(∆x)2

Φi−1,j − 2(1

(∆x)2+

1(∆y)2

)Φij +1

(∆x)2Φi+1,j

+1

(∆y)2Φi,j−1 +

1(∆y)2

Φi,j+1 = 0.

(diskrete Laplace-Gleichung)

Diese Gleichung gilt fur alle Punkte, die keine Randpunkte sind. Randpunkte sindaber nicht nur alle Punkte des außeren Randes, sondern auch alle Punkte im Gebiets-inneren, die Elektroden reprasentieren. Diese Randpunkte mussen getrennt behandeltwerden. Im elektrostatischen Problem sind zwei unterschiedliche Typen von Randbe-dingungen zu betrachten:

Dirichlet-Randbedingung: Gilt fur Punkte auf Elektroden; das Potenzial ist durcheine außere Spannung vorgegeben. Im Falle des Drei-Elektroden-Systems sind diesdie Punkte, die entweder auf Anodenspannung ΦA = 10V , auf KathodenspannungΦK = 2V oder auf Zwischenspannung ΦZ = 0V liegen.

Neumann-Randbedingung: Gilt fur Randlinien, bei denen die Ableitung des Po-tenzials vorgegeben ist. Im elektrostatischen Fall ist diese Ableitung senkrecht zumIsolator gleich Null. Dann nennt man die Randbedingung auch Symmetrie-Bedingung.Bei dem Drei-Elektroden-System nehmen wir an, dass die Elektroden durch ein di-elektrisches Material voneinander isoliert sind. Die physikalische Bedingung fur einenelektrischen Isolator ist, dass sich das elektrische Feld nur entlang aber nicht senkrechtzum Isolator ausbilden kann.

Am oberen und unteren Rand besitzt das elektrische Feld also keine Komponente iny-Richtung, d.h. die Potenziallinien verlaufen parallel zur y-Achse:

∂Φ

∂y= 0.

Diskretisiert bedeutet dies fur den unteren Rand

Φy(x, y = 0) ∼ Φi,2 − Φi,1

∆y= 0

und fur den oberen Rand

Φy(x, y = 20) ∼ Φi,jmax− Φi,jmax−1

∆y= 0.

Page 14: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

24 2 Modellierung und Simulationen mit finiten Differenzenverfahren

Nimmt man alle Feld- und Randpunkte mit ihren Bestimmungsgleichungen zusam-men, erhalten wir 41× 41 lineare Gleichungen fur die 41× 41 Gitterpunkte. Dies istwieder ein lineares Gleichungssystem fur die Potenzialwerte Φij an den Gitterpunkten(xij , yij). Die Losungen der Gleichungen sind dann die Potenzialwerte Φij .

2.3.3 Aufstellen des linearen Gleichungssystems

Um die Struktur des Gleichungssystems besser zu erkennen, wahlen wir im Folgendenein 6× 5-Gitter und zahlen die Gitterpunkte von links unten beginnend durch:

Abb. 2.10. Vereinfachtes Gitter.

Wir haben drei unterschiedliche Arten von Gitterpunkten

Dirichlet-Punkte 1, 7, 13, 19, 25 Potenzial ΦK = 2V6, 12, 18, 24, 30 Potenzial ΦA = 10V3, 9; 21, 27 Potenzial ΦA = 0V

Neumann-Punkte 2, 4, 5; 26, 28, 29Feld-Punkte Rest

An den 6× 5 = 30 Gitterpunkten stellen wir das LGS fur das Potenzial Φij auf. DerEinfachheit wegen setzen wir ∆x = ∆y = 1. Dann ist der diskrete Laplace-Operatoran der Stelle (i, j) gegeben durch: Der Wert des Potenzials an der unteren Stelle(i, j−1) + Potenzial an der oberen Stelle (i, j +1) + Potenzialwert rechts (i+1, j) +links (i− 1, j) minus 4 mal den Potenzialwert an der Stelle (i, j). Dieser Sachverhaltwird kurz durch den 5-Sterne-Operator

1

1 −4 1

1

Page 15: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.3 Das zweidimensionale elektrostatische Problem 25

geometrisch beschrieben. Die Bilanz uber alle Gitterpunkte ergibt1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

1 1 φK

2 -1 1 03 1 φZ

4 -1 1 05 -1 1 06 1 φA

7 1 φK

8 1 1 -4 1 1 09 1 φZ

10 1 1 -4 1 1 011 1 1 -4 1 1 012 1 φA

13 1 φK

14 1 1 -4 1 1 015 1 1 -4 1 1 016 1 1 -4 1 1 017 1 1 -4 1 1 018 1 φA

19 1 φK

20 1 1 -4 1 1 021 1 φZ

22 1 1 -4 1 1 023 1 1 -4 1 1 024 1 φA

25 1 φK

26 -1 1 027 1 φZ

28 -1 1 029 -1 1 030 1 φA

Dies ist das lineare Gleichungssystem, das im vierten Schritt durch geeignete Metho-den gelost wird.

2.3.4 Losen des linearen Gleichungssystems durch geeignete Methoden

Das durch die Diskretisierung der Laplace-Gleichung aufgestellte lineare Gleichungs-system kann man bis auf Rundungsfehler z.B. durch den Gauß-Algorithmus losen.Allerdings ist zu beachten, dass es bei hoherer Gitterauflosung sehr viele Unbekannteenthalt und zu sehr großen Systemen fuhrt. Fur ein 10× 10-Gitter hat man 100 Un-bekannte, was zu einer 100× 100-Matrix fuhrt. Fur diese Matrix benotigt man dannschon 10000 Speichereinheiten. Die Nachteile des Gauß-Algorithmus liegen somit aufder Hand:

– Rundungsfehler sind dominant, da man bei jedem Eliminationsschritt die Koeffizi-enten der verbleibenden Matrix neu berechnet.

– Hohe Rechenzeiten, da die Rechenzeit der Elimination proportional zu n! ist.

– Der Algorithmus verandert die Struktur der Matrix, da er in jedem Eliminations-schritt die Koeffizienten neu berechnet.

– Hoher Speicherbedarf; obwohl im ursprunglichen linearen Gleichungssystem vieleder Koeffizienten zunachst Null sind, werden sie im Verlauf der Verfahrens geandertund durch Zahlen ungleich Null ersetzt. Damit muss im Wesentlichen die gesamteMatrix abgespeichert werden.

Page 16: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

26 2 Modellierung und Simulationen mit finiten Differenzenverfahren

In der Praxis werden deshalb die linearen Gleichungssysteme meist mit iterativenMethoden gelost. Drei der klassischen Iterationsverfahren werden wir in diesem Ab-schnitt skizzieren. Es stellt sich dabei heraus, dass sie einfacher zu handhaben sindals der Gauß-Algorithmus. Die iterativen Methoden haben aber den Nachteil, dasssie nicht fur alle linearen Gleichungssysteme konvergieren, sondern nur fur bestimmtefur die Simulation aber relevanten! Durch diese Einschrankung werden diese Metho-den in der Elementarmathematik in der Regel nicht behandelt. Einen Uberblick ubergangige Verfahren zum Losen von großen linearen Gleichungssystemen findet man imAnhang A.

Wir betrachten im Folgenden ein lineares Gleichungssystem der Form

AΦ = b

mit der (n× n)-Matrix A, der rechten Seite b und dem gesuchten Losungsvektor Φ:

A = (aij)i=1,..,n;j=1,..,n; b = (b1, .., bn)t; Φ = (Φ1, .., Φn)t.

Das obige System besteht aus n Zeilen der Form

n∑j=1

aijΦj = bi (i = 1, ..., n).

2.3.4.1 Allgemeines Iterationsverfahren. Wir losen die i−te Gleichungen nachΦi auf (d.h. Gleichung 1 nach Φ1, Gleichung 2 nach Φ2, . . . , Gleichung n nach Φn)

Φi =1aii

(bi −i−1∑j=1

aijΦj −n∑

j=i+1

aijΦj) (i = 1, . . . , n)

und behandeln dieses System iterativ:

Wir geben uns eine Anfangsschatzung fur die Losung (Φ01, . . . , Φ

0n) des linearen Glei-

chungssystems vor. Im (m+1)-ten Schritt bestimmen wir eine Losung uber die Wertedes m-ten Schritts durch die allgemeine Vorschrift

Φ(m+1)i =

1aii

(bi −i−1∑j=1

aijΦ(m)j −

n∑j=i+1

aijΦ(m)j ) (i = 1, . . . , n)

(Jacobi-Verfahren)

Das Jacobi-Verfahren ist im Worksheet Jacobi sowohl direkt als auch in Form einerProzedur programmiert.

Page 17: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.3 Das zweidimensionale elektrostatische Problem 27

2.3.4.2 Gauß-Seidel-Verfahren. Betrachtet man das Jacobi-Verfahren, so beginntman bei der (m + 1)-ten Iteration mit der 1. Gleichung, um Φ

(m+1)1 durch die “alten“

Werte Φ(m)1 , . . . , Φ

(m)n zu berechnen. Anschließend wahlt man Gleichung 2 und berech-

net Φ(m+1)2 aus den “alten“ Daten Φ

(m)1 , . . . , Φ

(m)n . Dabei kann man aber im Prinzip

ausnutzen, dass Φ(m+1)1 schon berechnet wurde. Man erhalt also ein verbessertes Ver-

fahren, indem man generell bei der Berechnung von Φ(m+1)i die schon aktualisierten

Werte Φ(m+1)1 , . . . , Φ

(m+1)i−1 berucksichtigt:

Φ(m+1)i =

1aii

(bi −i−1∑j=1

aijΦ(m+1)j −

n∑j=i+1

aijΦ(m)j ) i = 1, . . . , n

(Gauß-Seidel-Verfahren)

Das Gauß-Seidel-Verfahren ist im Worksheet Gauß-Seidel sowohl direkt als auch inForm einer Prozedur programmiert.

2.3.4.3 SOR-Verfahren. Das SOR-Verfahren (Successive Overrelaxation) beruhtauf der Tatsache, dass die Konvergenz des Gauß-Seidel-Verfahrens erhoht werdenkann, wenn man eine Linearkombination des alten Wertes Φ

(m)i und des aktualisierten

Wertes Φ(m+1)i wahlt

Φ(neu)i = w · Φ(m+1)

i + (1− w) · Φ(m)i

(SOR-Verfahren)

Der Relaxationsparameter w liegt ublicherweise im Bereich 1 ≤ w ≤ 2. Man kannzeigen, dass 0 < w < 2 sein muss, da sonst keine Konvergenz erfolgt. Außerdemkann gezeigt werden, dass der optimale Relaxationsparameter gegeben ist durchwopt = 2

1+√

1−ρ2wenn ρ der großte Eigenwert der Jacobi-Matrix ist. Das SOR-

Verfahren ist im Worksheet SOR sowohl direkt als auch in Form einer Prozedurprogrammiert.

2.3.4.4 Abbruchkriterium. All diesen Iterationsverfahren ist gemeinsam, dass dieIteration abgebrochen werden muss. Um ein Abbruchkriterium zu erhalten, bestimmtman nach jeder Iteration den Fehler.

Da das lineare Gleichungssystem nicht exakt gelost wird, berechnet man in jedemIterationsschritt nur eine Naherung Φ(m+1). Damit ist A Φ(m+1) 6= b bzw.

A Φ(m+1) − b 6= 0.

Die Einzelfehler nach der (m + 1)-ten Iteration sind somit gegeben durch

Page 18: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

28 2 Modellierung und Simulationen mit finiten Differenzenverfahren

r(m+1)i =

n∑j=1

aijΦ(m+1)j − bi (i = 1, . . . , n)

bzw. der Gesamtfehler durch

R(m+1) = maxi=1,...,n

|r(m+1)i | (Residuum)

Als Abbruchkriterium fordert man in der Regel, dass sowohl

(1) das Residuum kleiner einer gewissen Vorgabe ist, R(m+1) < δ1, als auch dass

(2) die Differenz der Werte zweier aufeinander folgender Iterationen kleiner einer vor-gegebenen Schranke ist

maxi=1,...,n

|Φ(m+1)i − Φ

(m)i | < δ2.

2.3.5 Lineare Interpolation

Durch Losen des linearen Gleichungs-

10.4 9.6 9.4

10.4

10

10

9.4 8.2

8.4

9

Abb. 2.11. Konstruktion der 9.5-Volt-Linie.

systems z.B. mit iterativen Verfahrenerhalt man die gesuchte Losung nahe-rungsweise an den Gitterpunkten. D.h.die Losung wird reprasentiert durch Wer-te Φi,j an den Gitterpunkten (xi,j , yi,j).Zur Interpretation dieser numerischenDaten fuhrt man Aquipotenziallinien ein.Dies sind Linien, auf denen sich derWert des Potenzials nicht andert. DasProblem dabei ist, wie man aus den be-rechneten Daten auf den Gitterpunktenzu diesen Aquipotenziallinien kommt.

Zur Klarung dieser Frage betrachten wirAbbildung 2.11.

Vorgegeben sind die Potenzialwerte anden Gitterpunkten. Gesucht ist z.B. die 9.5-Volt-Linie. Um diese Linie zu finden,wahlt man sich einen Eckpunkt im Gitter (hier links oben mit 10 V) und betrachtetdie beiden angrenzenden Zell-Ecken (10 V nach unten; 9 V nach rechts). Damit istklar, dass die 9.5-Volt-Linie durch die obere Gitterlinie geht. Lineare Interpolationbesagt, dass sie genau durch die Mitte geht. D.h. die 9.5-Volt-Linie startet an deroberen Zell-Linie in der Mitte.

Page 19: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.3 Das zweidimensionale elektrostatische Problem 29

Anschließend werden alle weiteren Linien dieser Zelle untersucht und entschieden, obdie 9.5-Volt-Linie durch eine dieser Linien geht. In unserem Fall ist dies die untere Zell-begrenzung. Lineare Interpolation legt auf dieser unteren Linie den Wert 9.5 V fest.Auf die gleiche Weise wird nun die zweite Zelle durchsucht und die Schnittpunkte der9.5-Volt-Linie mit den Zellbegrenzungen bestimmt. Somit kommt die 9.5-Volt-Linievon der obersten Zelle durch die zweite Zelle bis hin zur untersten Zelle. Dort wirdfestgestellt, dass die Linie in die rechte Zelle geht usw.

Analog zum oben beschriebenen Verfahren bestimmt man auch die 10-Volt-, 9-Volt-bzw. 8.5-Volt-Linie usw. In Abbildung 2.12 (a) sind diese Linien eingezeichnet. Kom-merziell erhaltliche Programme farben die Losungsdarstellungen zusatzlich farblichein, wie dies in Abbildung 2.12 (b) in Grautonen zu sehen ist. Dies liefert aber kei-ne zusatzlichen Informationen, die nicht schon durch die Aquipotenziallinien gegebensind. Lediglich die Bereiche zwischen zwei Aquipotenziallinien werden mit einer Farbegefullt. Die eigentliche Information ist durch die Bereichsgrenzen also Aquipotenzial-linien gegeben.

(a) Aquipotenziallinien. (b) Farbliche Darstellung.

Abb. 2.12. Darstellung der Losung.

Beispiel 2.2 (Mit Maple). In Abbildung 2.13 ist die Losung fur das Drei-Elektroden-System auf einem 41 × 41-Gitter dargestellt. In den beiden Diagrammensind die Aquipotenziallinien fur die Randbedingungen (ΦA = 10V, ΦK = 4V undΦZ = 0V ) bzw. (ΦA = 10V, ΦK = 2V und ΦZ = 0V ) dargestellt. Man beachte, dassder einzige Unterschied im Kathodenpotenzial liegt: In Abbildung 2.13 (a) ist dasKathodenpotenzial ΦK = 4V und in 2.13 (b) ΦK = 2V .

Page 20: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

30 2 Modellierung und Simulationen mit finiten Differenzenverfahren

(a) Aquipotenziallinien fur ΦK = 4V . (b) Aquipotenziallinien fur ΦK = 2V .

Abb. 2.13. Numerische Losung fur das Drei-Elektroden-System.

Man erkennt im Vergleich der beiden Darstellungen sehr schon, dass eine kleine An-derung der Randbedingungen (der einzige Unterschied liegt in ΦK) den qualitativenVerlauf der Losung grundlegend andert. Im zweiten Fall erhalt man einen sog. Durch-bruch der Potenziallinien, der im ersten nicht moglich ist.

Zusammenfassung: Finite Differenzenmethode.

Wie im eindimensionalen Fall besteht die finite Differenzenmethode auch imzweidimensionalen Fall aus funf Schritten:

(1) Ersetzen der Geometrie durch diskrete Punkte: Die Geometrie wirddurch ein zweidimensionales strukturiertes Berechnungsgitter reprasentiert.

(2) Ersetzen der partiellen Ableitungen durch finite Differenzen: Diepartiellen Ableitungen der Differenzialgleichung werden durch finite Diffe-renzen approximiert. So kommt man bei linearen partiellen Differenzialglei-chungen zu linearen Differenzengleichungen.

(3) Aufstellen des zugehorigen LGS: Das lineare Gleichungssystem enthaltgenau so viele Gleichungen wie Unbekannte (= Anzahl der Gitterpunkte).

(4) Losen des LGS mit einer geeigneten Methode: In der Regel werden diedurch finite Differenzenverfahren entstehenden linearen Gleichungen durchiterative Methoden gelost.

(5) Lineare Interpolation zwischen den diskreten Punkten: Durch dasnumerische Losen der Differenzengleichungen erhalt man die Losung an denGitterpunkten. Durch lineare Interpolation erhalt man hieraus die Aquipo-tenziallinien.

Page 21: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.4 Verallgemeinerung 31

2.4 Verallgemeinerung

Durch die Differenzenmethode ist man prinzipiell in der Lage partielle Differenzial-gleichungen in den Ortskoordinaten x, y, und z auf achsenparallelen Gittern zu diskre-tisieren. Fur lineare Differenzialgleichungen fuhrt die Diskretisierung auf ein linearesGleichungssystem, das man in der Regel mit iterativen Methoden lost.

Beispiel 2.3. Am Beispiel der partiellen Differenzialgleichung

4 Φx + Φ + 3 Φxx + 4 Φyy = 10

werden wir aufzeigen, dass die Differenzenmethode auch auf eine solche Differenzial-gleichung anwendbar ist.

Zunachst stellen wir fest, dass die Diskretisierung der Differenzialgleichung unabhan-gig von der Geometrie ist, denn die Geometrie wird reprasentiert durch das diskrete,strukturierte Gitter. Dies bedeutet, dass jeder innere Gitterpunkt einen unteren, obe-ren, linken und rechten Nachbarpunkt besitzt. Wir konnen bei der Diskretisierungalso auf diese Gitterpunkte zuruckgreifen.

Wir erstellen fur die Ableitungen

4 Φx + Φ + 3 Φxx + 4 Φyy

im Punkt (i, j) den zugehorigen finiten Differenzenausdruck und geben fur ∆x =∆y = 1 die geometrische Anordnung des Differenzenoperator an.

4 Φx(i, j) ∼ 1∆x

(2Φi+1,j − 2Φi−1,j)

3Φxx(i, j) ∼ 1∆x2

(3Φi+1,j − 6Φi,j + 3Φi−1,j)

4Φyy(i, j) ∼ 1∆y2

(4Φi,j+1 − 8Φi,j + 4Φi,j−1).

Setzen wir diese Naherungen fur die Ableitungen in die Differenzialgleichung ein,erhalten wir im Punkt (i, j)

4 Φx + Φ + 3 Φxx + 4 Φyy ∼1

∆x(2Φi+1,j − 2Φi−1,j)

+ Φi,j

+1

∆x2(3Φi+1,j − 6Φi,j + 3Φi−1,j)

+1

∆y2(4Φi,j+1 − 8Φi,j + 4Φi,j−1) = 10.

Page 22: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

32 2 Modellierung und Simulationen mit finiten Differenzenverfahren

Zur ubersichtlicheren Schreibweise setzen wir fur ∆x = ∆y = 1. Dann erhalten wirim Punkt (i, j)

4 Φi,j−1 + 1 Φi−1,j − 13 Φi,j + 5 Φi+1,j + 4 Φi,j+1 = 10.

Dieser Sachverhalt wird wieder kurz durch einen 5-Sterne-Operator

4

1 −13 5

4

dargestellt.

Soll die partielle Differenzialgleichung fur das Gebiet aus Abbildung 2.10 gelost wer-den, so erhalten wir fur die angegebenen Randbedingungen das folgende lineare Glei-chungssystem fur die Losung auf den Gitterpunkten:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 301 1 φK

2 -1 1 03 1 φZ

4 -1 1 05 -1 1 06 1 φA

7 1 φK

8 4 1 -13 5 4 109 1 φZ

10 4 1 -13 5 4 1011 4 1 -13 5 4 1012 1 φA

13 1 φK

14 4 1 -13 5 4 1015 4 1 -13 5 4 1016 4 1 -13 5 4 1017 4 1 -13 5 4 1018 1 φA

19 1 φK

20 4 1 -13 5 4 1021 1 φZ

22 4 1 -13 5 4 1023 4 1 -13 5 4 1024 1 φA

25 1 φK

26 -1 1 027 1 φZ

28 -1 1 029 -1 1 030 1 φA

Dieses lineare Gleichungssystem wird im vierten Schritt durch iterative Methodengelost. Anschließend wird im funften Schritt durch lineare Interpolation die Losunggraphisch in Form von Aquipotenziallinien dargestellt.

Page 23: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

2.5 Aufgaben zur finiten Differenzenmethode 33

2.5 Aufgaben zur finiten Differenzenmethode

2.1 Gegeben ist die Differenzialgleichung

y′′(x) + y(x) = 0,

die auf dem Intervall [0, 5] mit den Randwerten y(0) = 1 und y(5) = 0 numerisch mitder finiten Differenzenmethode gelost werden soll.

a) Unterteilen Sie das Intervall [0, 5] in 5 Teilintervalle und stellen Sie das lineareGleichungssystem fur die Unbekannten yi = y(i) fur i = 1, . . . , 4 auf.

b) Losen Sie das LGS und skizzieren Sie die Losung.

2.2 Gegeben ist die Differenzialgleichung

y′′(x) + y′(x) = 0,

die auf dem Intervall [0, 5] mit den Randwerten y(0) = 1 und y(5) = 0 numerisch mitder finiten Differenzenmethode gelost werden soll.

a) Unterteilen Sie das Intervall [0, 5] in 5 Teilintervalle und stellen Sie das lineareGleichungssystem fur die Unbekannten yi = y(i) fur i = 1, . . . , 4 auf.

b) Losen Sie das LGS und skizzieren Sie die Losung.

2.3 Gegeben ist das lineare Gleichungssystem Ay = b mit

A =

−2 1 0 0

1 −2 1 00 1 −2 10 0 1 −2

und b =

1000

.

a) Losen Sie das LGS iterativ mit der Jacobi-Methode, indem Sie vier Iterationsschrit-te ausfuhren. Der Startvektor sei y(0) = (0, 0, 0, 0)t.

b) Losen Sie das LGS, indem Sie einen Iterationsschritt mit dem Gauß-Seidel-Verfahrendurchfuhren.

2.4 Gegeben ist das Potenzialproblem

Φxx(x, y) + Φyy(x, y) = 0,

welches fur die nachfolgende Geometrie (vgl. Abb. 2.14 (a)) numerisch mit der finitenDifferenzenmethode gelost werden soll.

a) Spezifizieren Sie fur jeden Gitterpunkt den Typ (Dirichlet-, Neumann- oder Feld-punkt.

b) Erstellen Sie das zugehorige LGS zum Losen der Potenzialgleichung unter der An-nahme, dass ∆x = ∆y = 1. (vgl. Abb. 2.14 (b))

c) Skizzieren Sie eine mogliche Losung, indem Sie vier Aquipotenziallinien in Abb.2.14 (a) qualitativ einzeichnen.

Page 24: Modellbildung und Simulation - home.hs-karlsruhe.deweth0002/buecher/simulation/pdf/kapitel2… · viii Inhaltsverzeichnis 6 Einf¨uhrung in ANSYS ..... 97 6.1 Die Benutzeroberfl¨ache

34 2 Modellierung und Simulationen mit finiten Differenzenverfahren

Abb. 2.14. Potenzialproblem.

2.5 Gegeben ist die partielle Differenzialgleichung

Φxx(x, y)− 4Φx(x, y) + 2 Φy(x, y) + Φyy(x, y) = 16,

die fur das Gebiet in Abb. 2.14 numerisch mit der finiten Differenzenmethode gelostwerden soll. Erstellen Sie das zugehorige LGS zum Losen der Potenzialgleichung unterder Annahme, dass ∆x = ∆y = 1.

2.6 Gegeben ist ein achsenparalleles aquidistantes Gitter mit den Abstanden ∆x = ∆y = 1.Erstellen Sie finite Differenzenausdrucke fur die Ableitungen

fx, fy, fxy, fxxyy,

wenn die Funktion fij nur an den Gitterpunkten (i, j) gegeben ist. Skizzieren Sie imGitter, welche Punkte fur die Berechnung der jeweiligen Ableitung mit einbezogenwerden mussen. Wie sind die zugehorigen Gewichte?

*2.7 Losen Sie folgende Problemstellung, indem Sie das zugehorige Maple-Worksheetverwenden:

a) Fuhren Sie fur ein Zwei-Elektroden-System ohne Zwischenelektrode ein (6 × 5)-Gitter ein (siehe Abb. 2.10) und erstellen Sie mit der finiten Differenzenmethodedas zugehorige lineare Gleichungssystem (ΦA = 10V ; ΦK = 2V ).

b) Losen Sie das unter (a) erhaltene lineare Gleichungssystem mit dem Jacobi-Verfahren, dem Gauß-Seidel-Verfahren und dem SOR-Verfahren (w = 1.3).

c) Losen Sie das lineare Gleichungssystem fur das Drei-Elektroden-System mit Zwi-schenelektrode (ΦA = 10V, ΦK = 2V, ΦZ = 0V ) mit dem Jacobi-Verfahren, demGauß-Seidel-Verfahren und dem SOR-Verfahren.

*2.8 Gegeben ist die Poisson-Gleichung im Innern des Einheitsquadrates I2

−uxx − uyy = −x2 − y2 + x + y

u(x, y) = 0 auf dem Rand von I2.

Die exakte Losung ist gegeben durch u(x, y) = 12

x(x− 1)y(y − 1).

a) Losen Sie die partielle Differenzialgleichung mit der Differenzenmethode fur dieWerte ∆x = ∆y = 1

10, ∆x = ∆y = 1

20, ∆x = ∆y = 1

40, ∆x = ∆y = 1

100.

b) Vergleichen Sie die numerische Losung mit der exakten Losung. Welche Aussagenuber die Genauigkeit kann man in Abhangigkeit der Schrittweite machen?