Optimierung der Verkn upfung uberlappender sph arischer ...Optimierung der Verkn upfung...

Optimierung der Verkn�upfung

�uberlappender sph�arischer Gitter f�ur

numerische Simulationen

Diplomarbeit

von

Marc Reuting

Institut f�ur Theoretische Physik, Lehrstuhl IVWeltraum- und AstrophysikRuhr-Universit�at Bochum

Oktober 2008

Inhaltsverzeichnis

1 Einleitung 11.1 Grundlagen des Yin-Yang-Gitters . . . . . . . . . . . . . . . . . . . . . . . 31.2 Ziel und Aufbau der Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Konstruktion des Yin-Yang-Gitters 72.1 Vor�uberlegungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 Die Transformation zwischen den Teilgittern Yin und Yang . . . . . . . . . 9

2.2.1 Konstruktion der Einheitsvektoren ri, #i, 'i . . . . . . . . . . . . . 92.2.2 Der Zusammenhang zwischen #1, '1 und #2, '2 . . . . . . . . . . . 102.2.3 R�ucktransformation . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Die Gitter in der #-'-Ebene . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3.1 Das Teilgitter Yin in der #-'-Ebene . . . . . . . . . . . . . . . . . . 122.3.2 Das Teilgitter Yang in der #-'-Ebene . . . . . . . . . . . . . . . . . 122.3.3 Das Yin-Yang-Gitter in der #-'-Ebene . . . . . . . . . . . . . . . . 13

3 Die MHD-Gleichungen 153.1 Beschreibung von Plasmen . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.2 Das Vielteilchenbild . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.3 Das Fluidbild . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.4 Die Momente der Vlasov-Gleichung . . . . . . . . . . . . . . . . . . . . . . 173.5 Die Maxwell-Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4 Grundlagen des numerischen Modells 214.1 Das Integrationsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.2 Das Rechengebiet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

5 Interpolation 255.1 Interpolationsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5.1.1 Lineare Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . 255.1.2 Interpolationspolynome . . . . . . . . . . . . . . . . . . . . . . . . . 255.1.3 Interpolation durch kubische Splines . . . . . . . . . . . . . . . . . 27

5.2 Interpolation auf dem Yin-Yang-Gitter . . . . . . . . . . . . . . . . . . . . 285.2.1 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 285.2.2 Implementation der Spline-Interpolation . . . . . . . . . . . . . . . 305.2.3 Implementation der Newton-Interpolation . . . . . . . . . . . . . . 325.2.4 Vergleich der Interpolationsmethoden . . . . . . . . . . . . . . . . . 34

i

Inhaltsverzeichnis

5.3 Test der Kopplung f�ur den Fall realistischer Simulationen . . . . . . . . . . 35

6 (M)HD-Gleichgewichte 396.1 Gleichgewicht ohne Magnetfeld und mit nur radialer Geschwindigkeit . . . 39

6.1.1 Annahmen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 396.1.2 L�osung des Gleichgewichts f�ur sph�arische Symmetrie . . . . . . . . 406.1.3 St�orung des Gleichgewichtes . . . . . . . . . . . . . . . . . . . . . . 41

6.2 Beispiel f�ur einen axialsymmetrischen Sonnenwind . . . . . . . . . . . . . . 466.2.1 Annahmen f�ur die Gleichgewichtsbedingung . . . . . . . . . . . . . 466.2.2 L�osung der Gleichgewichtsbedingung . . . . . . . . . . . . . . . . . 486.2.3 Das Gleichgewicht in den Koordinaten von Yin und Yang . . . . . . 56

7 Zusammenfassung und Ausblick 61

A Anhang 63A.1 Weitere M�oglichkeit zur Konstruktion von Gleichung (3.17) . . . . . . . . . 63A.2 Entwicklung von Gleichung (3.19) . . . . . . . . . . . . . . . . . . . . . . . 64A.3 Konstruktion von Gleichung (3.25) . . . . . . . . . . . . . . . . . . . . . . 65

Literaturverzeichnis 67

ii

1 Einleitung

Die Neugierde l�asst die Menschheit seit jeher danach streben, Naturph�anomene und ihreEntstehung zu hinterfragen. Die Entwicklung von Technologie f�uhrt h�au�g dazu, dass demMenschen der Alltag erleichtert wird. Vor allem die Physik ist f�ur das Verst�andnis vonZusammenh�angen in Natur und Technik ein wichtiges Teilgebiet der Naturwissenschaften.Je l�anger nach diesen Zusammenh�angen geforscht wird, desto komplexer werden die zuuntersuchenden Ph�anomene, und umso schwieriger wird es, sie bis in das kleinste Detailzu verstehen.Zu diesem Zweck werden inzwischen immer h�au�ger Computer eingesetzt. Mit ihrer Hilfelassen sich solche komplexen Ph�anomene numerisch simulieren. Diese Numerischen Simu-lationen ben�otigen allerdings eine Basis, auf welcher die Berechnungen angestellt werdenk�onnen. Es m�ussen zuerst die f�ur das zu untersuchende Problem relevanten physikalischenGr�o�en und Parameter de�niert werden.Einige L�osungsverfahren benutzen daf�ur ein an die Geometrie des Problems angepasstes

"Objekt\, d. h. man deckt den betrachteten Raum mit einer Art Netz ab, auf welchem imbestimmten Abstand voneinander Punkte ausgew�ahlt werden und auf denen die Gr�o�endann de�niert werden k�onnen. Das so entstandene numerische Gitter besteht aus vielenGitterpunkten, man beschreibt diese durch die Verwendung von Koordinatensystemen.Der Abstand der Gitterpunkte, d. h. die Au �osung des Gitters, und dessen geometrischeForm muss dem zu untersuchenden Problem angepasst sein. Um ein Ph�anomen hinrei-chend gut beschreiben zu k�onnen, muss der Abstand der Gitterpunkte sehr viel kleinersein als die L�angenskalen der betrachteten Objekte oder die �Anderungen der untersuchtenGr�o�en im System.Es ist dabei nicht immer sinnvoll, ein standardisiertes, gleichm�a�iges Gitter zur Beschrei-bung zu verwenden. Wenn die betrachteten Objekte der vorliegenden Problemstellungsich in ihrer Geometrie stark unterscheiden, kann es sinnvoll sein, verschiedene Gitter zuverwenden, d. h. das gesamte Gitter aus verschiedenen Teilgittern zusammenzusetzen.In Abbildung 1.1 erkennt man ein Gitter um den Querschnitt der Trag �ache eines Flug-zeugs mit zwei verstellbaren

"Klappen\ vor und hinter der eigentlichen Trag �ache. Der

Umriss der Trag �ache wird mit einem anderen Gitter abgedeckt als der der beiden Klap-pen, da diese sich e�ektiver mit einem anderen Koordinatensystem beschreiben lassen.Des Weiteren ist zu erkennen, dass die verwendeten Gitter einen �Uberlapp bilden. Diein dieser Abbildung verwendeten Gitter k�onnten beispielsweise dazu genutzt werden, dieLuftstr�omungen um die Trag �ache herum zu simulieren, um durch die gew�ahlte Form derTrag �ache den durch sie hervorgerufenen Auftrieb zu optimieren. Betrachtet man Abbil-dung 1.1, ist nachzuvollziehen, dass es in diesem Fall sinnvoll ist, das geometrische Objektmit mehreren Teilgittern abzudecken, also ein zusammengesetztes, �uberlappendes Gitterzu verwenden.Abbildung 1.2 zeigt ein Beispiel f�ur ein zusammengesetztes Gitter, welches im Ver-

1

1 Einleitung

Abbildung 1.1: Darstellung eines Gitters angepasst an den Querschnitt der Trag �ache einesFlugzeugs (Chesshire & Henshaw 1990)

Abbildung 1.2: Darstellung eines zusammengesetzten, nicht �uberlappenden sph�arischen Gitters("cubed sphere\, vgl. Koldoba et al. (2002))

gleich zu dem Gitter in Abbildung 1.1 lediglich ein einzelnes geometrisches Objekt -n�amlich eine Kugel - abdecken soll. Bei der Betrachtung von Abbildung 1.2 ist nichtsofort ersichtlich, warum man eine Kugel mit Hilfe von mehreren Teilgittern abdeckensollte. Eigentlich k�onnte man vermuten, dass sich Kugeln hinreichend gut mit einemnicht-zusammengesetzten Gitter in sph�arischen Polarkoordinaten abdecken lassen. DieBeschreibung mit Hilfe dieser Koordinatensysteme st�o�t aber schon an seine Grenzen,sobald E�ekte an den Polen der Kugel untersucht werden sollen. Der Grund daf�ur liegtin dem divergierenden Skalenfaktor h' der Di�erentialoperatoren in sph�arischen Polar-koordinaten (vgl. Kageyama (2005)). Eine e�ziente L�osung dieses Problems versprichtdie Verwendung des so genannten Yin-Yang-Gitters (vgl. Kageyama & Sato (2004)), beiwelchem die Kugel durch zwei unterschiedliche, gegeneinander verdrehte, strukturglei-

2

1.1 Grundlagen des Yin-Yang-Gitters

che Teilgitter abgedeckt wird. In dieser Arbeit konzentrieren wir uns auf dieses Gitter,mit welchem Objekte in sph�arischer Geometrie beschrieben werden k�onnen, da wir ander Beschreibung solcher Objekte interessiert sind. Viele Objekte der Weltraumphysik,wie beispielsweise die Sonne, lassen sich durch ihre sph�arische Geometrie sinnvoll durchsph�arische Koordinaten beschreiben. Ist dabei von

"Kugeln\ die Rede, so ist damit die

Vollkugel mit ihren drei Dimensionen in sph�arischen Polarkoordinaten (r,#,') mit demAbstand r vom Ursprung und den Winkeln #, welcher den Breitenwinkel auf der Kugelangibt, und ', welcher den L�angenwinkel angibt, gemeint. Mit

"Sph�are\ ist lediglich die

Kugelober �ache gemeint. Ein Gitter, welches durch sph�arische Polarkoordinaten beschrie-ben wird, ist in Abbildung 1.3 dargestellt.

Abbildung 1.3: Beispiel eines nicht-zusammengesetzten Gitters, welches die gesamte Kugelober- �ache abdeckt

Im Folgenden werden wir f�ur dieses Gitter den Begri�"herk�ommliches Gitter\ verwenden.

Soll nun die Sph�are durch das wie in Abbildung 1.3 dargestellte herk�ommliche Gitter ab-gedeckt werden, ergibt sich das angesprochene Problem des divergierenden Skalenfaktorsh' an den Polen (vgl. Kapitel 2 und Kageyama (2005)).

1.1 Grundlagen des Yin-Yang-Gitters

Die Idee ist, dass man die Kugelober �ache in zwei identische, aber gegeneinander verdrehteTeilgitter zerlegt, deren Ausdehnung so gew�ahlt ist, dass sie zusammen die gesamte Kuge-lober �ache abdecken und in ihren Randbereichen einen �Uberlapp bilden. Diese Teilgitterwerden Yin und Yang genannt - zwei Begri�e, die urspr�unglich aus der asiatischen Philo-sophie bekannt sind und hier als Namen f�ur die beiden Teilgitter verwendet werden. Daserste Gitter ist das nicht gedrehte Teilgitter Yin. Es ist dadurch gekennzeichnet, dass esdurch dasselbe Koordinatensystem repr�asentiert wird wie das herk�ommliche Gitter, ledig-lich der Bereich der beiden Winkel # und ' wird eingeschr�ankt (vgl. Gleichung (2.4)). Umdas Teilgitter Yang zu erhalten, dreht man das Teilgitter Yin auf die Weise (vgl. Kapitel2), dass es nach dem Zusammenf�ugen beider Teilgitter die komplette Sph�are abdeckt. InAbbildung 1.4 sind beide Teilgitter einzeln dargestellt.

3

1 Einleitung

Abbildung 1.4: a) Das Teilgitter Yin mit seinen sph�arischen Einheitsvektoren (links), b) DasTeilgitter Yang mit den zugeh�origen Einheitvektoren (rechts)

In Abbildung 1.4 a) ist das Teilgitter Yin mit den sph�arischen Einheitsvektoren (r,#,')und in Abbildung 1.4 b) das Teilgitter Yang mit seinen sph�arischen Einheitsvektoren undf�ur den Fall r = 1 dargestellt. F�ugt man beide Teilgitter zusammen, ergibt sich das inAbbildung 1.5 dargestellte zusammengesetzte Yin-Yang-Gitter:

Abbildung 1.5: Darstellung des zusammengesetzten Yin-Yang-Gitters

4

1.2 Ziel und Aufbau der Arbeit

1.2 Ziel und Aufbau der Arbeit

Das prim�are Ziel dieser Arbeit liegt darin, durch eine e�ektive Kopplung der beidenvon Kageyama & Sato (2004) eingef�uhrten Teilgitter deren N�utzlichkeit zu testen. DasProblem des verwendeten Yin-Yang-Gitters besteht darin, dass sich innerhalb des Re-chengebietes R�ander be�nden (vgl. Abbildung 1.5), welche physikalisch nicht existieren.Deshalb ist es notwendig, die beiden Teilgitter aneinander zu koppeln, w�ahrend der Vor-teil in der Beschreibung von E�ekten an den Polen liegen sollte (keine divergierendenSkalenfaktoren). Man muss daf�ur sorgen, dass sich die vorgegebenen Startwerte in dem�Uberlappungsbereich der beiden Teilgitter mit der Zeit nicht unterschiedlich entwickeln,die vorgegebenen Funktionswerte auf beiden Teilgittern d�urfen mit der Zeit nicht

"aus-

einander laufen\.Die Arbeit ist wie folgt aufgebaut. In Kapitel 2 wird die Konstruktion des Yin-Yang-Gitters erl�autert und die Form der Einheitsvektoren in Koordinaten von Yin und Yangberechnet. In Kapitel 3 werden die f�ur die Tests f�ur die Transformationen und die Kopp-lung verwendeten Gleichungen vorgestellt. Im Kapitel 4 werden kurz das numerische Mo-dell und das Rechengebiet umrissen. In Kapitel 5 wird dargestellt, wie man die beidenTeilgitter e�zient aneinander koppeln kann. Schlie�lich werden in Kapitel 6 die Testserl�autert und ihre Ergebnisse diskutiert. Kapitel 7 fasst die erhaltenen Ergebnisse nocheinmal zusammen.

5

1 Einleitung

6

2 Konstruktion des Yin-Yang-Gitters

2.1 Vor�uberlegungen

In Kapitel 1 wurde bereits erw�ahnt, dass bei Verwendung von sph�arischen Polarkoordina-ten der Skalenfaktor h' an den Polen der Kugel divergiert und somit numerische Modelle,die auf der Beschreibung durch solche Koordinatensysteme basieren, an ihre Grenzensto�en, sobald E�ekte an den Polen simuliert werden sollen. Um zu verstehen, wie mandieses Problem durch Verwendung des Yin-Yang-Gitters vermeiden kann, betrachten wirzun�achst sph�arische Polarkoordinaten und deren Di�erentialoperatoren.Der Zusammenhang zwischen Cartesischen Koordinaten (x,y,z) und sph�arischen Polar-koordinaten (r,#,') wird durch folgende Gleichung repr�asentiert (vgl. Bronstein et al.(2000)): 0@ x

yz

1A =

0@ r sin# cos'r sin# sin'r cos#

1A ; (2.1)

wobei der Bereich f�ur den Winkel # 2 [0; �] und f�ur den Winkel ' 2 [0; 2�] ist.

Der vektorielle Di�erentialoperator ~r hat in Cartesischen Koordinaten folgende Gestalt:

~r =

�@

@x;@

@y;@

@z

�(2.2)

Die Komponenten des Operators enthalten in Cartesischen Koordinaten (x,y,z) die Ab-leitungen in die Koordiantenrichtungen. Nach Transformation in sph�arische Polarkoordi-naten ergibt sich f�ur den Operator ~r (siehe z. B. Huba (2006)) :

~r = (rr;r#;r') =

�@

@r;1

r

@

@#;

1

r sin#

@

@'

�(2.3)

An Gleichung (2.3) ist ersichtlich, dass der Skalenfaktor in '-Richtung (h') an den beidenPolen (Nordpol: # = 0, S�udpol: # = �) divergiert:

h'(# = 0; # = �)!1;

Im Folgenden wird nun diskutiert, wie die Verwendung des Yin-Yang-Gitters dieses Pro-blem l�osen soll. Der Index 1 kennzeichnet im Folgenden das Teilgitter Yin und der Index2 das Teilgitter Yang. Das Teilgitter Yin ist - wie bereits in Kapitel 1 erl�autert - durchfolgende Winkelbereiche gekennzeichnet:

#1=2 2��

4;3�

4

�und '1=2 2

��

4;7�

4

�(2.4)

7


Das Teilgitter Yang erh�alt man durch Drehung des Teilgitters Yin in Cartesichen Koor-dinaten (x, y, z) zuerst um 'D = 180� um die z-Achse, und anschlie�end um #D = 90�

um die x-Achse.Mithilfe von Drehmatrizen l�asst sich die Transformation von Yin nach Yang folgender-ma�en darstellen:

0@ x2y2z2

1A =

0@1 0 00 cos'D sin'D0 � sin'D cos'D

1A0@ cos#D sin#D 0� sin#D cos#D 0

0 0 1

1A0@ x1y1z1

1A=

0@1 0 00 0 10 �1 0

1A0@�1 0 00 �1 00 0 1

1A0@ x1y1z1

1A =

0@ �x1z1y1

1AZusammenfassend gilt also:

0@ x2y2z2

1A =

0@ �x1z1y1

1A (2.5)

Gleichung (2.5) liefert einen Zusammenhang zwischen den Koordinaten auf Yin und Yang.Nutzt man nun den Zusammenhang zwischen Cartesischen Koordinaten und Kugelkoordi-naten aus Gleichung (2.1), so ergibt sich folgender Zusammenhang zwischen den Winkeln#1, '1 in Yin und den Winkeln #2, '2 in Yang:

sin#2 cos'2 = � sin#1 cos'1 (2.6a)

sin#2 sin'2 = cos#1 (2.6b)

cos#2 = sin#1 sin'1 ; (2.6c)

so dass man insgesamt das �uberlappende sph�arische Yin-Yang-Gitter erh�alt (vgl. Abbil-dung 1.5).Die Gr�o�e des �Uberlapps l�asst sich berechnen, indem man die Fl�ache jedes der beidenGitter berechnet und mit der Kugelober �ache vergleicht:

A =

Z 3=4�

�=4

Z 7=4�

1=4�

sin# d# d' � 2; 21� (2.7)

Der �Uberlapp f�ur das Yin-Yang-Gitter ergibt sich demnach zu: A�Uberlapp � 0; 42�

8

2.2 Die Transformation zwischen den Teilgittern Yin und Yang


2.2.1 Konstruktion der Einheitsvektoren ri, #i, 'i

Die f�ur dieses Kapitel verwendeten Identit�aten wurden dem Buch von Neutsch (1995)

entnommen. Den Zusammenhang zwischen den sph�arischen Einheitsvektoren r, # und 'und den Cartesischen Einheitsvektoren x, y sowie z liefern die folgenden Gleichungen:

r = sin# cos' x+ sin# sin' y + cos# z (2.8a)

# = cos# cos' x+ cos# sin' y � sin# z (2.8b)

' = � sin' x+ cos' y (2.8c)

Andererseits lassen sich die Cartesischen Einheitsvektoren durch die sph�arischen mitHilfe der folgenden Formeln berechnen:

x = sin# cos' r + cos# cos' #� sin'' (2.9a)

y = sin# sin' r + cos# sin' #+ cos'' (2.9b)

z = cos# r � sin# # (2.9c)

Wie wir anhand von Gleichung (2.5) gesehen haben, h�angen die Gitter Yin und Yangin Cartesischen Koordinaten auf folgende Weise zusammen:0@ x2

y2z2

1A =

0@ �x1z1y1

1ANun werden wir mit Hilfe der Gleichungen (2.5) und (2.8a) - (2.8c) einen Zusammenhangzwischen den sph�arischen Einheitsvektoren auf Yin und den Einheitsvektoren auf Yangherstellen:

r2 = sin#2 cos'2 x2 + sin#2 sin'2 y2 + cos#2 z2 (2.10a)

#2 = cos#2 cos'2 x2 + cos#2 sin'2 y2 � sin#2 z2 (2.10b)

'2 = � sin'2 x2 + cos'2 y2 (2.10c)

Setzt man nun Gleichung (2.5) in Gleichungen (2.10a) - (2.10c) ein, ergibt sich:

r2 = � sin#2 cos'2 x1 + sin#2 sin'2 z1 + cos#2 y1 (2.11a)

#2 = � cos#2 cos'2 x1 + cos#2 sin'2 z1 � sin#2 y1 (2.11b)

'2 = sin'2 x1 + cos'2 z1 (2.11c)

9


Benutzt man nun die Gleichungen (2.9a) - (2.9c), so ergeben sich durch Einsetzen inGleichungen (2.11a) - (2.11c) die folgenden Beziehungen:

r2 = (� sin#2 cos'2 sin#1 cos'1 + sin#2 sin'2 cos#1 + cos#2 sin#1 sin'1) r1+(� sin#2 cos'2 cos#1 cos'1 � sin#2 sin'2 sin#1 + cos#2 cos#1 sin'1) #1+(sin#2 cos'2 sin'1 + cos#2 cos'1) '1

(2.12a)

#2 = (� cos#2 cos'2 sin#1 cos'1 + cos#2 sin'2 cos#1 � sin#2 sin#1 sin'1) r1+(� cos#2 cos'2 cos#1 cos'1 � cos#2 sin'2 sin#1 � sin#2 cos#1 sin'1) #1+(cos#2 cos'2 sin'1 � sin#2 cos'1) '1

(2.12b)'2 = (sin'2 sin#1 cos'2 + cos'2 cos#1) r1

+(sin'2 cos#1 cos'1 � cos'2 sin#1) #1+(� sin'2 sin'1) '1

(2.12c)

Mit den Gleichungen (2.12a) - (2.12c) ist ein Zusammenhang zwischen den sph�arischen

Einheitsvektoren auf Yin (r1,#1,'1) und den Einheitsvektoren auf Yang (r2,#2,'2) mitHilfe der Winkel #1 und '1 auf dem Gitter Yin, sowie der Winkel #2 und '2 auf Yanghergestellt.Nun sollen die Einheitsvektoren auf einem der beiden Teilgitter berechnet wer-den, wenn nur die eines der beiden Winkelpaare (#1; '1) oder (#2; '2) und die zugeh�origenEinheitsvektoren auf dem jeweils anderen Gitter bekannt sind.Ersetzt man nun mit Hilfe der Identit�aten (2.6a) - (2.6c) alle Terme, die die Winkel #1und '1 enthalten, so erh�alt man Beziehungen mit denen sich die Einheitsvektoren aufYang nur durch Gr�o�en auf Yin berechnen lassen. Es ergeben sich folgende Beziehungenf�ur die Einheitsvektoren auf Yang in Abh�angigkeit von Gr�o�en auf Yin:

r2 = r1 (2.13)

#2 = 1p1�sin2 #1 sin2 '1

h� cos#1 sin'1 #1 � cos'1 '1

i(2.14)

'2 = 1p1�sin2 #1 sin2 '1

hcos'1 #1 � cos#1 sin'1 '1

i(2.15)

2.2.2 Der Zusammenhang zwischen #1, '1 und #2, '2

Abschlie�end in diesem Kapitel sollen die Winkel eines Teilgitters nur aus denen des je-weils anderen berechnet werden.Gleichung (2.6c) hat uns bereits eine Gleichung f�ur # geliefert:

cos#2 = sin#1 sin'1 (2.16)

sie lautet also:

10


#2 = arccos (sin#1 sin'1) (2.17)

Teilt man die beiden Identit�aten (2.6b) und (2.6c) durcheinander, l�asst sich folgenderZusammenhang f�ur ' ableiten:

'2 = arctan

�cos#1

� sin#1 cos'1

�(2.18)

Dabei ist das Argument der Funktion nicht als Quotient der beide Zahlen zu verstehen,sondern als:

'2 = atan2 [cos#1;� sin#1 cos'1]

Wobei die Funktion atan2[a; b] eine Funktion aus FORTRAN darstellt, welche zwei Ar-gumente a und b besitzt. Die Vorzeichen vor den Argumenten sind wie die Phase einerkomplexen Zahl aufzufassen.

2.2.3 R�ucktransformation

Nun sind sowohl die sph�arischen Einheitsvektoren auf Yang in Abh�angigkeit von denenauf Yin bekannt, als auch die beiden Winkel #2 und '2. Die Berechnung der R�ucktrans-formation erfolgt v�ollig analog, wie man schon an der Symmetrie der Gleichungen (2.6a)- (2.6a) erkennen kann. Die R�ucktransformationen lauten:

r1 = r2 (2.19)

#1 = 1p1�sin2 #2 sin2 '2

h� cos#2 sin'2 #2 � cos'1 '2

i(2.20)

'1 = 1p1�sin2 #2 sin2 '2

hcos'2 #2 � cos#2 sin'2 '2

i(2.21)

#1 = arccos (sin#2 sin'2) (2.22)

'1 = arctan

�cos#2

� sin#2 cos'2

�(2.23)

11


2.3 Die Gitter in der #-'-Ebene

2.3.1 Das Teilgitter Yin in der #-'-Ebene

Stellt man das Teilgitter Yin in die #-'-Ebene projiziert dar, handelt es sich in den Ko-ordinaten von Yin per constructionem um ein Rechteck (Abbildung 2.1).Man erkennt in Abbildung 2.1 die schon in Abschnitt 1.1 angesprochenen Grenzen des

Abbildung 2.1: Das Gitter Yin in die #-'-Ebene projiziert

Teilgitters Yin. In Abbildung 2.1 wird der Winkel # von oben nach unten und der Winkel' von links nach rechts gez�ahlt.

2.3.2 Das Teilgitter Yang in der #-'-Ebene

Nachdem wir das Teilgitter Yin auf einfache Weise in seinen Koordinaten dargestellt ha-ben, soll nun das Teilgitter Yang ebenfalls in den Koordinaten von Yin in die #-'-Ebeneprojiziert dargestellt werden. Dazu m�ussen die Winkel #2 und '2 des Teilgitters Yang mitHilfe der Beziehungen (2.17) und (2.18) in die Koordinaten von Yin transformiert werden.Es ergibt sich das in Abbildung 2.2 dargestellte Gebiet.Zus�atzlich eingezeichnet sind hier die Intervallgrenzen der Winkel # und ' des in rot dar-gestellten Teilgitters Yang. Weiterhin ist zu beachten, dass die Winkel in der Darstellungvon Yin nicht mehr von oben nach unten, sondern wie in Abbildung 2.2 eingezeichnetgez�ahlt werden. Die Pole der Kugel be�nden sich aber wie gew�ohnlich bei # = 0 und# = �

12

2.3 Die Gitter in der #-'-Ebene

Abbildung 2.2: Das Gitter Yang in die #-'-Ebene projiziert

2.3.3 Das Yin-Yang-Gitter in der #-'-Ebene

Zeichnet man nun beide Gitter zusammen, so ergibt sich Abbildung 2.3.In dieser Abbildung ist erneut das Teilgitter Yin blau und das Teilgitter Yang rot ge-

Abbildung 2.3: Die Gitter Yin und Yang in die #-'-Ebene projiziert

kennzeichnet.Man erkennt in Abbildung 2.3, dass es keine Stelle in der Ebene gibt, die nicht von min-

13


destens einem der beiden Teilgitter abgedeckt wird. An mehreren Stellen jeweils am Randder beiden Gitter kann man einen �Uberlapp erkennen. An den Stellen, an welchen kein�Uberlapp entsteht, grenzen die Teilgitter direkt aneinander.

14

3 Die MHD-Gleichungen

In Kapitel 1 und 2 wurde das Yin-Yang-Gitter vorgestellt. Wie bereits in Kapitel 1 an-gesprochen, soll dieses insbesondere in Bezug auf die Kopplung der beiden Teilgitter aufdie N�utzlichkeit f�ur numerische Simulationen getestet werden, deren Anwendungen imBereich der Weltraumphysik liegen. Die hier untersuchten Ph�anomene aus der Weltraum-phsik werden durch Magnetohydrodynamische (MHD) Gleichungen beschrieben. DieseGleichungen werden im Folgenden kurz erl�autert. Wir folgen an dieser Stelle zumeist derDarstellung von Chen (1974).

3.1 Beschreibung von Plasmen

"Unter einem Plasma versteht man ein quasineutrales Gas aus geladenen und ungeladenenTeilchen, welches kollektives Verhalten zeigt.\(Chen 1974)

F�ur die Beschreibung von Plasmen wird hier der Ansatz des Fluidbildes verwendet. Umdiesen Ansatz verwenden zu k�onnen, muss das Plasma bestimmte Bedingungen erf�ullen:es m�ussen die internen L�angenskalen des Plasmas (z. B. Gyrationsradien und vor allem dieDebyel�ange) sehr klein sein gegen�uber den L�angenskalen der makroskopischen �Anderungen(wie beispielsweise die �Anderung der Teilchendichte oder des Magnetfeldes). Man vollziehtan dieser Stelle den �Ubergang vom Einzelteilchen-Ensemble zum kontinuierlichen Fluid.Um diesen �Ubergang durchf�uhren zu k�onnen, geht man zun�achst �uber zum so genanntenVielteilchenbild.

3.2 Das Vielteilchenbild

Nehmen wir an, es bef�anden sich N Teilchen im System. Im Vielteilchenbild ist man nichtam genauen Verhalten aller N Teilchen interessiert. Das Verhalten aller N Teilchen zuuntersuchen w�are aufgrund der gro�en Anzahl der Teilchen im System weder m�oglich nochsinnvoll. Das Ziel ist es nicht, das Verhalten aller Teilchen im Detail zu kennen, sonderndas Verhalten des Plasmas in seiner Gesamtheit. Zu diesem Zweck f�uhrt man die Funktionef ein, welche die Position der N Teilchen im Phasenraum beschreibt. Die Koordinatendes Phasenraums sind der Ort ~x und die Geschwindigkeit ~u der Teilchen des Plasmas:

ef( ~x1; : : : ; ~xN ; ~u1; : : : ; ~uN ; t)Diese Funktion ist zu verstehen als �Uberlagerung der L�osung der Bewegungsgleichungenvon N Teilchen. Sie ist also abh�angig von den Ortsvektoren ~xi und den Geschwindig-keitsvektoren ~ui aller N Teilchen und der Zeit t, also von 6N + 1 Gr�o�en. Diese hohe

15


Dimensionalit�at des Problems l�asst sich allerdings stark reduzieren, wenn man die Ver-teilungsfunktion f einf�uhrt, welche die Wahrscheinlichkeit daf�ur angibt, ein Teilchen aneinem bestimmten Ort ~x bei einer Geschwindigkeit ~u vorzu�nden:

f(~x; ~u; t)

Diese Verteilungsfunktion h�angt nur noch ab von 7 unabh�angigen Variablen. Weiter redu-zieren l�asst sich die Zahl der Variablen, wenn die Verteilungsfunktion �uber den Geschwin-digkeitsraum integriert wird:

n(~x; t) =

Zf(~x; ~u; t) d3u (3.1)

Man erh�alt dann die Teilchenzahldichte n des Plasmas. Der �Ubergang zu makroskopischenGr�o�en, wie der Teilchenzahldichte n, der Str�omungsgeschwindigkeit ~v des Plasmas oderdem Druck P , l�asst sich dann verstehen als �Ubergang zum Fluidbild. Die Str�omungsge-schwindigkeit ~v des Plasmas ist gegeben durch:

~v =

Z~uf(~x; ~u; t) d3u (3.2)

3.3 Das Fluidbild

Man ist nun an der zeitlichen Entwicklung dieser Verteilungsfunktion interessiert. UnterVernachl�assigung von St�o�en der Teilchen untereinander ergibt sich aus dem Liouville-Theorem

df

dt= 0; (3.3)

welches besagt, dass in diesem Fall die totale Zeitableitung der Verteilungsfunktion ver-schwindet. Schreibt man die totale zeitliche Ableitung als Summe der partiellen Ableitun-gen, so ergibt sich:

df

dt=@f

@t+@f

@~x� @~x@t

+@f

@~u� @~u@t

= 0 (3.4)

Mit der De�nition der Lorentz-Kraft ~FL

~FL = q�~E + ~v � ~B

�; (3.5)

wobei q die Ladung des Teilchens ~B das Magnetfeld und ~E das elektrische Feld darstellen,ergibt sich schlie�lich die so genannte Vlasov-Gleichung:

df

dt=@f

@t+ ~u � @f

@~x+

q

m

�~E + ~v � ~B

�� @f@~u

= 0 (3.6)

Betrachtet man nun Momente dieser Vlasov-Gleichung ergeben sich daraus die MHD-Gleichungen. Das nullte Moment (M0) der Verteilungsfunktion bzw. einer Gleichung ist

16

3.4 Die Momente der Vlasov-Gleichung

dabei de�niert als:

M0(~x; t) =

Zf(~x; ~u; t) d3u ; (3.7)

Es exitieren allerdings noch weitere Momente:

~M1(~x; t) =

Z~uf(~x; ~u; t) d3u ; (3.8)

wobei anzumerken ist, dass ~M1 ein Vektor ist.

M2 =

Z[~u~u] f(~x; ~u; t) d3u (3.9)

Das Moment M2 ist ein Tensor und [~u~u] ein Tensorprodukt. Weitere Momente k�onnendann analog berechnet werden.

3.4 Die Momente der Vlasov-Gleichung

Die MHD-Gleichungen (siehe z. B. Chen (1974)) erh�alt man aus der Addition der erstendrei Momentengleichungen f�ur Elektronen und Ionen .

0-tes Moment

Aus dem 0-ten Moment ergibt sich die Bilanzgleichung f�ur die Massendichte � (oder Kon-tinuit�atsgleichung):

@�

@t+ ~r � (�~v) = 0 (3.10)

1-tes Moment

Aus dem 1-ten Moment folgt die Bilanzgleichung f�ur die Impulsdichte �~v:

@(�~v)

@t+ ~r � (�[~v~v])� ~J � ~B + ~rP = 0 (3.11)

Gleichung (3.11) ist keine Momentengleichung, sondern ergibt sich erst aus der Summevon Elektronen- und Ionengleichung. Das Produkt [~v~v] ist erneut ein Tensorprodukt; auch

der Druck P ist im allgemeinen ein Tensor, ~B ist das �au�ere und vom Plasma selbst in-duzierte Magnetfeld, ~J stellt die Stromdichte dar.

17


2-tes Moment

Aus dem 2-ten Moment ergibt sich die Bilanzgleichung f�ur die innere (thermische) Ener-giedichte ":

@(�")

@t+ ~r � (�"~v) + ~r � ~q + P (~r � ~v) = 0 (3.12)

In Gleichung (3.12) stellt der zweite Term die Konvektion dar, der dritte Term die W�arme-leitung und der vierte Term die adiabatische Energie�anderung. F�ur ein idales, polytro-pisches Gas gilt folgender Zusammenhang zwischen der inneren Energiedichte " und derTemperatur T (vgl. Priest (1982)):

" = cV T ; (3.13)

wobei cV die spezi�sche W�armekapazit�at bei konstantem Volumen ist.

cV =1

� 1

kBm

; (3.14)

mit der Boltzmannkonstanten kB. Der Adiabaten-Exponent ist de�niert als das Verh�alt-nis der spezi�schen W�armekapazit�aten bei konstantem Druck und konstantem Volumen,weiterhin h�angt er auch auf folgende Weise von der Zahl nf der Freiheitsgrade ab:

:=Cp

CV=nf + 2

nf

Zwei wichtige F�alle sind die eines monoatomaren Gases mit nf = 3, also = 5=3 undeines isothermen Gases mit = 1. Nutzt man nun das ideale Gasgesetz

P =�

mkBT (3.15)

folgt daraus

�" =P

� 1: (3.16)

Vernachl�assigt man den W�armestrom ~q erh�alt man die im Simulationscode verwendeteForm dieser Gleichung (3.17):

@P

@t+ ~r � (P~v) + ( � 1)P (~r � ~v) = 0 (3.17)

f�ur den Fall, dass die Entropie erhalten ist, kann Gleichung (3.17) auch aus dieser abge-leitet werden, dies wird im Anhang gezeigt. Da Gleichung (3.17) aber nicht konservativist, nimmt man noch eine Ersetzung vor (vgl. Otto (1987)):

P = u (3.18)

18

3.5 Die Maxwell-Gleichungen

Mit dieser Ersetzung l�asst sich Gleichung (3.17) schreiben als:

@u

@t+ ~r � (u~v) = 0 (3.19)

Eine genauere Darstellung dieser Ersetzung ist erneut im Anhang nachzulesen.

3.5 Die Maxwell-Gleichungen

Neben den Gleichungen (3.10), (3.11) sowie (3.19) m�ussen die elektromagetischen Felder~E und ~B zus�atzlich noch die vier Maxwell-Gleichungen erf�ullen

~r� ~E = �@~B

@t(3.20)

~r� ~B = �0

~J + "0

@ ~E

@t

!(3.21)

~r � ~E = �="0 (3.22)

~r � ~B = 0; (3.23)

wobei ~J die Stromdichte und � die Ladungsdichte ist. Gleichung (3.22) wird in der MHDin der Regel vernachl�assigt, da es sich bei einem Plasma um ein quasineutrales Gas han-delt und somit � � 0 gilt.Der Term �0"0@ ~E=@t in Gleichung (3.21) gibt den Verschiebungstrom an und kann hiervernachl�assigt werden, da die Plasmageschwindigkeit sehr klein ist gegen�uber der Vakuum-Lichtgeschwindigkeit c = 1=

p�0"0 = 3 � 108m=s (siehe z.B. Priest (1982)). Dann ergibt

sich f�ur Gleichung (3.21):

~J =1

�0(~r� ~B)

Zus�atzlich wird das Ohmsche Gesetz genutzt, um das elektrische Feld ~E aus den Glei-chungen zu eliminieren (siehe beispielsweise Priest (1982)). Es lautet

~E + ~v � ~B = � ~J; (3.24)

wobei � die Resistivit�at des Plasmas ist. Mit den Maxwell-Gleichungen (3.20) und (3.23)kombiniert, ergibt sich die Induktionsgleichung:

@ ~B

@t=�~r�

�~v � ~B

��+

�

�0� ~B � ~r� � ~J (3.25)

Wie man Gleichung (3.25) erh�alt, ist erneut im Anhang nachzulesen.Abschlie�end werden alle Gleichungen, welche f�ur die Simulation verwendet werden auf-gef�uhrt in ihrer endg�ultigen Form aufgef�uhrt. Diese wird ausf�uhrlich bei Laukert (2008)

19


erl�autert.

@�

@t+ ~r � (�~v) = 0 (3.26)

@(�~v)

@t+ ~r � (�[~v~v])� ~J � ~B + ~ru = 0 (3.27)

@ ~B

@t+ ~r� (~v � ~B) +

�

�0� ~B � ~r� � ~J = 0 (3.28)

@u

@t+ ~r � (u~v) = 0 ; (3.29)

wobei Gleichung (3.27) im Falle einer nicht-verschwindenden Gravitationsbeschleunigungg noch um einen Term erweitert werden muss:

@(�~v)

@t+ ~r � (�[~v~v])� ~J � ~B + ~ru = �gr (3.30)

Der Vektor r ist erneut der Einheitsvektor in r-Richtung in sph�arischen Polarkoordinaten.Zus�atzlich ist anzumerken, dass die im Code verwendeten Gleichungen auf typische Gr�o�ennormiert sind (vgl. Laukert (2008)).

20

4 Grundlagen des numerischen Modells

4.1 Das Integrationsverfahren

In Kapitel 3 wurden die MHD-Gleichungen erl�autert. Nun soll beschrieben werden, wiediese in dem von Laukert (2008) entwickelten Code gel�ost werden. Als Integrationsver-fahren f�ur die hyperbolischen Gleichungen (3.26), (3.27) und (3.29) verwenden wir dasvon Laukert (2008) detailliert beschriebene Leapfrog-Verfahren. Betrachten wir z. B. dieeindimensionale, hyperbolische Di�erentialgleichung:

@y

@t+ c

@y

@x= 0 ; (4.1)

mit der gesuchten Funktion y = f(x; t) abh�angig vom Ort x in einer Dimension und derZeit t sowie einer beliebigen Konstanten c. Die Funktionswerte y = f(x; t) sind au�erzum Zeitpunkt t = 0 nicht in kontinuierlicher Weise bekannt, sondern liegen bei einernumerischen Simulation nur in diskreter Form vor, d. h. Funktionswerte sind nur auf dendiskreten Gitterpunkten bekannt. Die Diskretisierung entlang x und t erfolge �aquidistant,xi sei der i-te Gitterpunkt und tn der n-te Zeitpunkt. Die Gr�o�e �x = xi+1 � xi gibtdann den Gitterabstand an und �t = tn+1� tn den Zeitschritt. Das Gitter wird zun�achstin nx Gitterpunkte, d. h. nx � 1 Gitterabst�ande eingeteilt, mit einem ersten Gitterpunktxmin = x(1) und einem letzten Gitterpunkt xmax = x(nx), so dass i 2 [1; nx] gilt. F�ur deni-ten Gitterpunkt folgt dann:

xi = x(i) = xmin + (i� 1)�x (4.2)

Aus x(nx) = xmax = xmin + (nx � 1)�x folgt f�ur den Gitterabstand �x:

�x =xmax � xmin

nx � 1(4.3)

Also insgesamt:

xi =xmax � xmin

nx � 1(i� 1) + xmin (4.4)

Diskretisiert man Gleichung (4.1) nach dem Leapfrog-Verfahren, so erh�alt man:

yn+1i = yn�1i � c

�t

�x(yni+1 � yni�1) (4.5)

Das Leapfrog-Verfahren ben�otigt - um den Funktionswert yn+1i aus einer Gleichung der

Form (4.5) zu einem neuen Zeitpunkt tn+1 an einer Stelle xi berechnen zu k�onnen - dieFunktionswerte yni+1 und y

ni�1 an den Nachbargitterpunkten xi+1 und xi�1. Auf den Rand-

punkten kann also nicht integriert werden, da der Funktionswert y auf einem der beiden

21


Nachbarpunkte unbekannt ist. Man f�ugt deswegen an beiden R�andern einen weiterenGitterpunkt - den sogenannten numerischen Rand - an.

4.2 Das Rechengebiet

Das Rechengebiet soll nun um zwei Gitterpunkten erweitert werden, d. h. aus nx Gitter-punkten werden n0x = nx + 2 Gitterpunkte. Das physikalische Rechengebiet durchl�auftjetzt i 2 [2; n0x � 1]. Die neuen Grenzen verschieben sich mit x0min = xmin � �x undx0max = xmax +�x. Daraus folgt f�ur den i-ten Gitterpunkt x0i:

x0i = x0(i) = x0min + (i� 1)�x = (xmin ��x) + (i� 1)�x (4.6)

(4.7)

Der Gitterabstand �x ergibt sich:

�x =xmax � xmin

nx � 1=xmax � xmin

n0x � 3(4.8)

Und damit f�ur den i-ten Gitterpunkt xi:

xi =xmax � xmin

n0x � 3(i� 2) + xmin (4.9)

In Abbildung 4.1 sind zun�achst die beiden Teilgitter Yin und Yang einzeln mit ihremnumerischen Rand dargestellt.Das Zusammenf�ugen beider Teilgitter ergibt Abbildung 4.2.

Abbildung 4.1: a) Das Teilgitter Yin in der #-'-Ebene mit zus�atzlichem numerischen Rand(gr�un, links), b) Das Teilgitter Yang in der #-'-Ebene mit zus�atzlichem nume-rischen Rand (gelb, rechts)

22

4.2 Das Rechengebiet

Abbildung 4.2: Das Yin-Yang-Gitter in der #-'-Ebene mit zus�atzlichem numerischen Rand, Yinund Yang sind erneut in blau und rot dargestellt, der numerische Rand von Yingr�un und der numerische Rand von Yang gelb

Durch das Anf�ugen dieses numerischen Randes wird der �Uberlapp an den Stellen gr�o�er,an denen auch vorher schon ein �Uberlapp war und an den Stellen, an welchen die Teilgit-ter genau aneinander passten, entsteht nun erst ein �Uberlapp.Bei bekannten Funktionswerten auf dem numerischen Rand lassen sich die Funktionswerteauf einem Gitterpunkt mit dem Leapfrog-Verfahren berechnen. Diese Werte sind allerdingsunbekannt. Der Vorteil bei der Verwendung des Yin-Yang-Gitters liegt nun darin, dassdie fehlenden Funktionswerte durch Interpolation aus den Werten auf dem jeweils ande-ren Teilgitter berechnet werden k�onnen und nicht extrapoliert werden m�ussen. Denn diebeiden Teilgitter decken zusammen die ganze Sph�are ab und es entstehen keine Gebiete,die nicht von einem der beiden Teilgitter abgedeckt werden. Eine ausf�uhrliche Darstellungdieser Interpolation erfolgt im folgenden Kapitel 5. Wie man die Funktionswerte auf demnumerischen Rand berechnet, wird im folgenden Kapitel 5 erl�autert.

23


24

5 Interpolation

5.1 Interpolationsverfahren

Die Problemstellung bei der Interpolation ist, dass man eine reellwertige Funktion y =f(x) abh�angig von einer reellen Variablen x nur an diskreten Stellen - St�utzstellen genannt- kennt. Bekannt seien (n+1) diskrete, paarweise verschiedene St�utzstellen x0 < � � � < xnund die zugeh�origen St�utzwerte y0; : : : ; yn. Die Aufgabe ist es nun die Funktion zwischenden St�utzstellen geeignet zu approximieren.Es werden nun exemplarisch einige Interpolationsverfahren erl�autert, dabei wird zumeistder Darstellung von Schwarz (1986) gefolgt.

5.1.1 Lineare Interpolation

Bei diesem Interpolationsverfahren werden zwischen den St�utzstellen jeweils st�uckwei-se Geraden aneinandergesetzt. Man l�ost dieses Interpolationsaufgabe also

"lokal\, indem

man den Funktionswert an der unbekannten Stelle durch Verbinden der beiden Funkti-onswerte an den benachbarten St�utzstellen erh�alt. Die lineare Interpolationsformel ergibtsich zu:

y =x� xi

xi+1 � xiyi+1 +

xi+1 � x

xi+1 � xiyi (5.1)

Hier steht y f�ur den unbekannten Funktionswert an der Stelle x. Die Punkte xi+1 und xisind die zur Stelle x benachbarten St�utzstellen mit den Funktionswerten yi+1 bzw. yi miti 2 [0; : : : ; n].

5.1.2 Interpolationspolynome

Lagrange-Interpolation

Es besteht aber auch die M�oglichkeit die Interpolationsaufgabe nicht"lokal\, sondern

"global\, also auf dem gesamten betrachteten Intervall zu l�osen. Gesucht wird ein Poly-nom n-ten Grades

Pn(x) = a0 + a1x+ � � �+ anxn ; (5.2)

welches die Interpolationsbedingungen

Pn(xi) = yi (i = 0; 1; : : : ; n) (5.3)

erf�ullt. Dieses Polynom existiert und ist eindeutig, siehe hierzu Schwarz (1986).Die Berechnung dieses Interpolationspolynoms ist jedoch nicht eindeutig. Die folgende

25

5 Interpolation

Darstellung nennt man die Lagrangesche Interpolationsformel, nach dem italienischenMathematiker und Astronom Joseph-Louis Lagrange (1736-1813, geboren als GiuseppeLudovico Lagrangia):

Pn(x) =nXi=0

yi(x� x0) : : : (x� xi�1)(x� xi+1) : : : (x� xn)

(xi � x0) : : : (xi � xi�1)(xi � xi+1) : : : (xi � xn)(5.4)

Man kann erkennen, dass mit der Anzahl der verwendeten St�utzstellen der Grad des Po-lynoms ansteigt. Gleichung (5.4) liefert also eine explizite Darstellung des Polynoms. F�urdie Interpolation nach Lagrange formt man dieses Polynom aber noch algebraisch um, um

den Rechenaufwand zu verk�urzen. Dazu werden die sogenannten St�utzkoe�zienten �(k)i

de�niert, die rekursiv berechnet werden k�onnen (siehe dazu Schwarz (1986)):

Start: �(0)0 = 1

Rekursion: f�ur k = 1; 2; : : : ; n :

f�ur i = 0; 1; : : : ; k � 1 :

�(k)i =

�(k�1)i

xi � xk

�(k)k = �

k�1Xi=0

�(k)i (5.5)

Aus den �(n)i = �i k�onnen die Gr�o�en �i berechnet werden:

�i :=�i

x� xi(i = 0; : : : ; n) (5.6)

Schlie�lich erh�alt man damit die Formel der Lagrange-Interpolation zur Berechnung deszu interpolierenden Wertes y = Pn(x):

Pn(x) =

Pni=0 �iyiPni=0 �i

(5.7)

Newton-Interpolaton

Eine weitere Form der Darstellung f�ur das Polynom Pn(x) ist die Newtonsche Interpola-tionsformel, nach dem englischen Physiker, Mathematiker und Astronom Isaac Newton(1643-1727):

Pn(x) = c0 + c1(x� x0) + c2(x� x0)(x� x1) + c3(x� x0)(x� x1)(x� x2) + : : :

+cn(x� x0)(x� x1) : : : (x� xn) (5.8)

Die unbekannten Koe�zienten c0; c1; : : : ; cn lassen sich sukzessive aus den Interpolations-bedingungen berechnen:

26

5.1 Interpolationsverfahren

Pn(x0) = c0 = y0Pn(x1) = c0 + c1(x1 � x0) = y1Pn(x2) = c0 + c1(x2 � x0) + c2(x2 � x0)(x2 � x1) = y2usw.

(5.9)

5.1.3 Interpolation durch kubische Splines

Ein weiteres Verfahren, welches die Interpolationsaufgabe als eine Art"Hybridverfahren\,

also weder"global\ noch

"lokal\ l�ost, ist die Interpolation durch kubische Splines. Hierbei

werden wieder - wie bei der linearen Interpolation - st�uckweise Funktionen zwischen denSt�utzstellen aneinandergesetzt, allerdings nun so, dass die Gesamtl�osung im Gegensatzzur linearen Interpolation stetig di�erenzierbar bleibt (vgl. Schwarz (1986)).Gesucht wird nun eine mindestens einmal stetig di�erenzierbare Interpolationsfunktions(x). Dazu geht man von dem Modell aus, dass durch die gegebenen St�utzpunkte eineeine d�unne, homogene Latte - woher auch der Name f�ur die Methode kommt,

"spline\ ist

der englische Begri� f�ur"Latte\ - gelegt sei, die in den St�utzpunkten gelenkig gelagert sei

und dort keinen �au�eren Kr�aften unterliege. Dann soll die Biegelinie der Latte die L�osungs(x) der Interpolationsaufgabe sein.Diese L�osung wird bestimmt, indem man nach Extremalprinzipien die Deformations-energie der Latte durch ihre angenommene Form minimiert (vgl. Schwarz (1986)). AlsAnsatz w�ahlt man zwischen jeweils zwei St�utzstellen ein kubisches Polynom mit dennoch zu bestimmenden Gr�o�en ai; bi; ci; di, man erh�alt daraus schlie�lich n� 1 Polynome(i = 1; : : : ; n) bei der Verwendung von n St�utzstellen

si(x) = ai(x� xi)3 + bi(x� xi)

2 + ci(x� xi) + di (5.10)

Diese k�onnen mit Hilfe folgender Beziehungen berechnet werden (siehe Schwarz (1986))

ai =1

6hi(y00i+1 � y00i ) (5.11)

bi =1

2y00i (5.12)

ci =1

hi(yi+ 1� yi)� 1

6hi(y

00

i+1 + 2y00i ) (5.13)

di = yi ; (5.14)

wobei

hi = (xi+1 � xi) (5.15)

und

yi = diyi+1 = aih

3i + 2bih

2i + cihi + di

y00i = 2biy00i+1 = 6aihi + 2bi :

27

5 Interpolation

Mit diesen Gleichungen ist das Vorgehen bei der kubischen Spline-Interpolation vorgege-ben und die gesuchten Polynome k�onnen sukzessive berechnet werden. Die Vorteile diesesweder

"lokalen\ noch

"globalen\ Verfahrens gegen�uber den Verfahren, die

"lokal\ - wie die

lineare Interpolation - sind, dass die L�osungsfunktion der Interpolationsaufgabe auch anden St�utzstellen di�erenzierbar bleibt und gegen�uber den

"globalen\ Verfahren, die mit

dem Interpolationspolynom arbeiten ist, dass die L�osungsfunktion auf jedem Intervall einPolynom dritten Grades unagh�angig von der Anzahl der verwendeten St�utzstellen bleibtund somit nicht zu Oszillationen bei der Verwendung vieler St�utzstellen neigt.

5.2 Interpolation auf dem Yin-Yang-Gitter

5.2.1 Randbedingungen

An dieser Stelle sei noch einmal an das Ende von Kapitel 4 erinnert. Dort wurde erl�autert,dass das Leapfrog-Verfahren die Funktionswerte an den Nachbargitterpunkten ben�otigt.Die Randpunkte der jeweiligen Teilgitter besitzten aber nicht in beide Richtungen Nach-barpunkte, was zu der Einf�uhrung der numerischen Randes f�uhrt. Wie man die Funkti-onswerte auf dem numerischen Rand per Interpolation aus dem jeweils anderen Gitterberechnen kann, soll im Folgenden erl�autert werden.Um die genaue Vorgehensweise bei der Interpolation im einzelnen zu verstehen, betrach-ten wir erneut Abbildung 2.3. In dieser Abbildung wird das Gitter in der #-'-Ebene mitdem numerischen Rand gezeigt.Die beiden Teilgitter Yin und Yang wurden so gew�ahlt, dass die Gitterpunkte �aquidistantverteilt sind. Es wird an dieser Stelle exemplarisch lediglich ein Wert auf dem #-Randberechnet, die Berechnung der '-Werte erfolgt analog.

Zun�achst wird berechnet, wo der Randpunkt in dem jeweils anderen Teilgitter liegt. Dazumuss beispielweise der Winkel #2 des entsprechenden Gitterpunktes auf dem numerischenRand im Teilgitter Yang aus den bekannten Winkeln #1 und '1 im Teilgitter Yin berech-net werden. Die Berechnung erfolgt anhand der Transformationsformel (2.6c):

#2 = arccos (sin#1 sin'1)

Aus diesem Winkel wird der Gitterindex i berechnet. Dazu benutzen wir Gleichung (4.9),mit x! #:

#i =#max � #min

n# � 3(i� 2) + #min (5.16)

Dabei ist allerdings zu beachten, dass dieser berechnete Gitterindex durch die Umrechnungnun nicht mehr notwendigerweise ganzzahlig ist, deswegen ersetzen wir i durch r# 2 IR(vgl. Abbildung 5.1).

r# =#2 � #min#max � #min

(n# � 3) + 2

28


Die Gr�o�en #min = �4bzw. #max =

3�4sind die vorher schon erw�ahnten Intervallgrenzen

f�ur den Winkel #, n# ist die Anzahl der Gitterpunkte in #-Richtung.Die ersten beiden Nachbargitterpunkte �ndet man, indem man als ersten zu verwen-

Abbildung 5.1: Der Rand der beiden Gitter

denden Gitterpunkt i#2 den Punkt nimmt, der dem ganzzahligen Anteil der berechnetenStelle entspricht:

i#2 = [r#] (5.17)

Die erste St�utzstelle nennen wir i#2 , da f�ur die Interpolation nach Newton und die Splin-einterpolation insgesamt vier St�utzstellen verwendet werden sollen. Den n�achstgr�o�erenNachbargitterpunkt i#3 erh�alt man durch:

i#3 = i#2 + 1 (5.18)

Damit h�atte man die beiden direkten Nachbargitterpunkte links und rechts der interpolie-renden Stelle gefunden, auf welchen man die zugeh�origen Funktionswerte kennt und mankann bereits linear interpolieren. Genauer wird die Interpolation allerdings, wenn man aufbeiden Seiten jeweils einen Gitterpunkt hinzuf�ugt, welche als dritte (i#1) und vierte (i#4)dienen. Diese ergeben sich aus:

i#1 = i#2 � 1 (5.19)

i#4 = i#3 + 1 (5.20)

F�ur den Fall, dass einer der beiden �au�eren Gitterpunkte i#1 oder i#4 auf einen nume-rischen Randpunkt f�allt, was nur an wenigen Stellen geschehen kann, an welchen beideTeilgitter genau aneinander grenzen (siehe Abbildung 2.3), werden alle Punkte verscho-ben (i#k ! i#k + 1, bzw. i#k ! i#k � 1).Es l�asst sich also insgesamt festhalten, dass f�ur die Interpolation in der Simulation ent-weder zwei (lineare Interpolation) oder vier St�utzstellen verwendet werden. Die zu inter-polierende Stelle liegt aufgrund des Verfahrens immer zwischen den St�utzstellen 2 und3, es sei denn die Stellen m�ussen um �1 verschoben werden, weil sie genau auf einemRandpunkt liegen. Prinzipiell lassen sich nun alle Funktionswerte durch Interpolation derFunktionswerte auf dem jeweils anderen Gitter berechnen. Nutzt man die �Aquidistanz des

29

5 Interpolation

Gitters aus, so lassen sich die L�osungen der Interpolationsaufgabe explizit ausrechnen. Esergeben sich dann Polynome in den bekannten Di�erenzen der zu interpolierenden Stel-le von den St�utzstellen. Die nach der Spline-Interpolation und der Newton-Interpolationerhaltenen L�osungen werden im Folgenden dargestellt.

5.2.2 Implementation der Spline-Interpolation

Bei der Spline-Interpolation mit Hilfe von vier St�utzstellen ergeben sich drei Polynomeentsprechend der Intervalle, diese drei Polynome werden s1, s2 und s3 genannt (vgl. Ab-bildung 5.2). Berechnet man die Vorfaktoren ai, bi, ci und di aus Gleichung (5.10) ergibtsich f�ur die Vorfaktoren zu dem Polynom s1(x), welches die L�osung zwischen der erstenund zweiten St�utzstelle repr�asentiert:

a1 = �4(f(i#2)� f(i#1))

15h3+f(i#3)� f(i#2)

3h3� f(i#4)� f(i#3)

15h3(5.21)

b1 = 0 (5.22)

c1 =19(f(i#2)� f(i#1))

15h� f(i#3)� f(i#2)

3h+f(i#4)� f(i#3)

5h(5.23)

d1 = f(i#1) ; (5.24)

wobei hi = xi+1 � xi = h im Fall �aquidistanter St�utzstellen. Die Vorfaktoren zu demPolynom s2(x) lauten:

a2 =f(i#2)� f(i#1)

5h3� 2(f(i#3)� f(i#2))

3h3+f(i#4)� f(i#3)

5h3(5.25)

b2 = �f(i#2)� f(i#1)

5h2+f(i#3)� f(i#2)

h2� 4(f(i#4)� f(i#3))

5h2(5.26)

c2 =3(f(i#2)� f(i#1))

5h+5(f(i#3)� f(i#2))

3h+7(f(i#4)� f(i#3))

15h(5.27)

d2 = f(i#2) (5.28)

Die Vorfaktoren zu dem Polynom s3(x) ergeben sich zu:

a3 = �f(i#2)� f(i#1)

15h3+f(i#3)� f(i#2)

3h3� 4(f(i#4)� f(i#3))

15h3(5.29)

b3 = �f(i#2)� f(i#1)

5h2+f(i#3)� f(i#2)

h2� 4(f(i#4)� f(i#3))

5h2(5.30)

c3 = �2(f(i#2)� f(i#1))

15h+2(f(i#3)� f(i#2))

3h+7(f(i#4)� f(i#3))

15h(5.31)

d3 = f(i#3) (5.32)

Um Rechenzeit zu sparen, bringen wir die Polynome auf die Form

si(#2) = pi1(�)f(i#1) + pi2(�)f(i#2) + pi3(�)f(i#3) + pi4(�)f(i#4) (5.33)

30


mit nur vom Gitter abh�angigen Polynomen pji die an der Stelle

� =#2 � #(i#2)

#(i#3)� #(i#2)(5.34)

ausgewertet werden, so dass die pji nur zu Beginn der Simulation werden m�ussen. DieVorfaktoren pji ergeben:

p11 =4

15�3 � 19

15� + 1

p12 = �3

5�3 +

8

5�

p13 =2

5�3 � 2

5�

p14 = � 1

15�3 +

1

15� (5.35)

p21 = �1

3�3 +

4

5�2 � 7

15�

p22 = �3 � 9

5�2 � 1

15� + 1

p23 = ��3 + 6

5�2 +

4

5�

p24 =1

3�3 � 1

5�2 � 2

15� (5.36)

p31 =1

15�3 � 1

5�2 +

2

15�

p32 = �2

5�3 +

6

5�2 � 4

5�

p33 =3

5�3 � 9

5�2 +

1

5� + 1

p34 = � 4

15�3 +

4

5�2 +

7

15� (5.37)

Von der Lage der zu interpolierenden Stelle #2 h�angt ab, welches der Polynome s1; : : : ; s3zu w�ahlen ist: Aufgrund der Wahl der St�utzstellen liegt die zu interpolierende Stelle #im Normalfall zwischen den St�utzstellen i#2 und i#3 , woraus folgt, dass das Polynom s2zu w�ahlen ist (vgl. Abbildung 5.2). Es sei denn, einer der beiden St�utzstellen i#1 oder i#4liegt auf einem numerischen Randpunkt, dann folgt daraus eine Verschiebung der Stellen(i#k ! i#k � 1, bzw. i#k ! i#k + 1). Werden beispielsweise die vier Stellen i#k ! i#k + 1

31

5 Interpolation

verschoben, so hat das eine �Anderung der Wahl des zu verwendenden Polynoms zufolge.In diesem Fall w�are das Polynom s1 zur Interpolation zu w�ahlen. Der Wert von � �andertsich nicht (Gleichung (5.39)), da das Gitter �aquidistant ist. Die L�osungen im Fall der

Abbildung 5.2: Wahl des richtigen Polynoms si

Verwendung der Spline-Interpolation ergeben sich aus den L�osungen si aus Gleichungen(5.33) mit den angegebenen Vorfaktoren pji.

5.2.3 Implementation der Newton-Interpolation

Wird statt der Splineinterpolation die Interpolation nach Newton verwendet, ergibt sichzun�achst das folgende Polynom:

Pn(r#) = f(i#1)

+(r# � i#1)

�f(i#1)� f(i#2)

h+ (r# � i#2)

�f(i#3)� 2f(i#2) + f(i#1)

2h2

��(r# � i#1)

�(r# � i#2)

�(r# � i#3)(�f(i#4) + 3f(i#3)� 3f(i#2) + f(i#1))

6h3

��Nutzt man nun die �Aquidistanz des Gitter und bezieht alle verwendeten Gitterpunkte i#auf die St�utzstelle i#2

r# � i#2 = �h (5.38)

r# � i#1 = i#2 + �h� (i#2 � h) = h(� + 1)

r# � i#3 = i#2 + �h� (i#2 + h) = h(� � 1) ;

Die Wahl der Gr�o�e � im Fall der Newton-Interpolation stimmt mit dem � der Spline-Interpolation �uberein, denn es gilt:

� =#� #(i#2)

#(i#3)� #(i#2)=

r# � i#2i#3 � i#2

=r# � i#2

h(5.39)

32


Daraus ergeben sich erneut eine Form wie in (5.33) f�ur das die L�osung Pn:

Pn = p1(�)f(i#1) + p2(�)f(i#2)(�) + p3(�)f(i#3) + p4(�)f(i#4) (5.40)

Die vier Vorfaktoren ergeben sich zu:

p1 = �1

6�3 +

1

2�2 � 1

3�

p2 =1

2�3 � �2 � 1

2� + 1

p3 = �1

2�3 +

1

2�2 + �

p4 =1

6�3 � 1

6�

Aufgrund von Gleichung (5.39) ergibt sich f�ur den Fall der Verwendung der Newton-Interpolation ohne Verschiebung:

0 � � � 1 (5.41)

Wird allerdings eine Verschiebung vorgenommen, gilt:

�1 � � � 0 (5.42)

bzw. 1 � � � 2 (5.43)

In Abbildung 5.3 sind die vier Polynome p1, p2, p3 und p4 in Abh�angigkeit von �, alsoabh�angig von der Lage der Neustelle dargestellt.

Abbildung 5.3: Die Polynome p1 (rot), p2 (blau), p3 (gr�un), p4 (t�urkis) in Abh�angigkeit von �,die waagerechte gestrichelte Linie markiert den Funktionswert p = 1

Man erkennt in Abbildung 5.3, dass f�ur den Vorfaktor von f(i#1) = 1 - also p1 (roter

33

5 Interpolation

Graph) - gilt, f�ur den Fall � = �1, was gleichbedeutend mit der Tatsache #(r#) = #(i#1)ist. F�ur alle anderen Vorfaktoren gilt in diesem Fall pi = 0; i = 2; 3; 4. Entsprechen-des gilt f�ur die weiteren drei F�alle: #(r#) = #(i#2), #(r#) = #(i#3), #(r#) = #(i#4). F�urden Fall 0 < � < 1 �uberwiegen beispielsweise die Beitr�age der Vorfaktoren p2 und p3gegen�uber denen von p1 und p4. Je nach Gr�o�e von � unterscheiden die Beitr�age derVorfaktoren pi. Die im Code vorkommende L�osung der Interpolationsaufgabe bei Ver-wendung der Newton-Interpolation ergibt sich demnach aus Gleichung (5.40). Die beidenhier vorliegenden Koordinaten, also die #- sowie die '-Richtung werden als unabh�angigvoneinander angenommen und somit v�ollig analog behandelt. Die zweidimensionale In-terpolation ergibt sich einfach als Produkt der beiden L�osungen. Die zweidimensionaleInterpolationsaufgabe wird also gel�ost, indem wir die zweidimensionale Interpolation alszwei mal eindimensionale Interpolation aufassen.

5.2.4 Vergleich der Interpolationsmethoden

In diesem Abschnitt sollen die im letzten Abschnitt beschriebenen, umgeformten Verfah-ren - die Splineinterpolation und die Interpolation nach Newton - bei der Verwendungvon vier St�utzstellen untersucht und mit der linearen Interpolation verglichen werden.Als eine erste Testfunktion soll die Sinus-Funktion f(x) = sinx (Abbildung 5.4) dienen.

Abbildung 5.4: a) Darstellung der Sinus-Funktion f(x) = sinx auf dem Intervall [0; 0:01] (links),b) Darstellung der logarithmierten Di�erenz der interpolierten Funktionswer-te zu den bekannten Funktionswerten von f(x) = sinx (durchgezogene Li-nie: Newton-Interpolation, gepunktet: Spline, gestrichelt: lineare Interpolation,rechts)

In Abbildung 5.4 ist die logarithmierte Di�erenz der interpolierten Funktionswerte zu denbekannten Funktionswerten von f(x) = sinx dargestellt (durchgezogene Linie: Newton-Interpolation, gepunktet: Spline, gestrichelt: lineare Interpolation). Es ist zu erkennen,dass die maximale Di�erenz bei Verwendung der linearen Interpolation mit exp(�11) �2 � 10�5 (ergibt sich aus der logarithmierten Darstellung) am gr�o�ten ist. Als zweitgr�o�te

34

5.3 Test der Kopplung f�ur den Fall realistischer Simulationen

Abweichung ergibt sich die nach der Splineinterpolation mit etwa exp(�12) � 6 �10�6. DieDi�erenz bei Verwendung der Interpolation nach Newton ist mit etwa exp(�20) � 2 �10�9um einiges geringer.

Als n�achste Testfunktion dient die in Abbildung 5.5 dargestellte Exponentialfunktionf(x) = ex.Die Wahl der Linienarten ist dieselbe wie in Abbildung 5.4. In Abbildung 5.5 b) ist zu

Abbildung 5.5: a) und b) Die Darstellung ist analog zu der in Abbildung 5.4, allerdings mit derTestfunktion f(x) = ex

erkennen, dass die maximale Di�erenz erneut bei Verwendung der linearen Interpolationmit etwa exp(�9) � 1 � 10�4 erneut am gr�o�ten ist. Als zweitgr�o�te Abweichung ergibtsich die nach der Splineinterpolation mit etwa exp(�10) � 5 � 10�5. Die Di�erenz beiVerwendung der Interpolation nach Newton ist erneut mit etwa exp(�16) � 1 � 10�7 umeiniges geringer.

Die letzte Testfunktion ist f(x) = x2ex (Abbildung 5.6):Die Linienarten wurden gew�ahlt wie in Abbildung 5.4. Auch hier best�atigt sich das Er-gebnis der anderen Tests, allerdings ist die Abweichung mit etwa exp(�8) � 3 � 10�4insgesamt gr�o�er und bei allen Verfahren von �ahnlicher Gr�o�enordnung. Und man kannzusammenfassend festhalten, dass bei den hier verwendeten Testfuntionen f�ur die Inter-polationsverfahren angepasst an die Bedingungen auf dem Gitter die Interpolation nachNewton die geringsten Abweichungen aufweist, w�ahrend die lineare Interpolation am un-genauesten ist. Ingesamt bleiben die Abweichungen aber mit maximal 3% (enspricht demDi�erenzen bezogen auf die ensprechenden Funktionswerte) vernachl�assigbar.


Simuliert wurde hier die Ausbreitung von Alfv�en �ugeln, welche auf�uhrlich bei Laukert(2008) diskutiert werden. Alfv�en �ugel entstehen, wenn ein Plasma senkrecht zu einemMagnetfeld str�omt und dabei einen leitenden K�orper umstr�omen muss.

35

5 Interpolation

Abbildung 5.6: a) und b) Die Darstellung ist analog zu der in Abbildung 5.4, hier mit derTestfunktion f(x) = x2ex

Ein Beispiel zum Vergleich der linearen Interpolaltion, welche zwei St�utzstellen verwendet,zur Newtoninterpolation bei Verwendung von vier St�utzstellen mit wirklichen Simulati-onsrechnungen liefern Abbildungen 5.7 und 5.8. Dargestellt ist hier das Dichtepro�l beifestem Radius r der kompletten #-'-Ebene. Zu sehen ist das vollst�andige sph�arische Re-chengebiet, abgedeckt durch beide Teilgitter Yin und Yang. Wenn eine hinreichend guteKopplung beider Teilgitter durch die Verwendung von Interpolationsverfahren gew�ahr-leistet ist, sollte der Rand zwischen beiden Teilgittern, welcher dadurch entstanden ist,dass ein zusammengesetztes Gitter verwendet wurde, nicht sichtbar sein.In Abbildung 5.7 kann man deutlich diesen Rand erkennen, was darauf hinweist, dass

Abbildung 5.7: Dichte � zum Zeitpunkt mit Hilfe linearer Interpolation

36


eine zufriedenstellende Kopplung der Teilgitter nicht gew�ahrleistet ist. Es ist noch anzu-merken, dass in den Abbildungen n = �

mgilt, mit m = 1 folgt n = �.

Wird allerdings mit Hilfe der Newton-Intepolation interpoliert, ergibt sich daraus eine

Abbildung 5.8: Dichte � in der #-'-Ebene mit Hilfe der Interpolation nach Newton

erhebliche Verbesserung der Kopplung der beiden Teilgitter (vgl. Abbildung 5.8).In Abbildung 5.9 ist die Dichte � in Cartesischer Darstellung und zu einem sp�ateren Zeit-punkt zu erkennen.In der hier dargestellten x-z-Ebene bei y = 0 ist der gr�o�ere Teil des Rechengebietes vom

Abbildung 5.9: Dichte � in Cartesischer Darstellung, b) Dichte � ebenfalls in Cartesischer Dar-stellung (gegl�attet)

Teilgitter Yin abgedeckt. Der kleinere Teil, welcher sich im Bereich x < 0 be�ndet und

37

5 Interpolation

dort ein Dreieck bildet, wird von dem Teilgitter Yang abgedeckt. In Abbildung 5.9 a) istdieser Bereich deutlich sichtbar. Die Dichte ist dort etwas erh�oht im Vergleich zu ihrerUmgebung, zu erkennen ist dies an der st�arkeren Rotf�arbung dieses Bereichs. Man kannin Abbildung 5.9 a) also erneut die Teilgitter Yin und Yang erkennen.Verbessern l�asst sich die Kopplung der Teilgitter, in dem man die Funktionswerte derphysikalischen Gr�ossen, wie z. B. die Dichte �, auf dem gesamten Gitter gl�attet.Hierzu werden zuerst die Werte auf beiden Teilgittern auf das bei Laukert (2008) be-schriebene

"sph�arisches Gitter\ transformiert. Anschlie�end werden die Funktionswerte

auf diesem sph�arischen Gitter gegl�attet, indem jeder Funktionswert aus dem Funktions-wert an dieser Stelle und den Funktionswerten auf den Nachbargitterpunkten gemitteltwird. Abschlie�end werden alle Funktionswerte wieder zur�ucktransformiert. Das Ergebnisnach dieser Gl�attung ist in Abbildung 5.9 b) dargestellt. Man erkennt, dass die Kopplungdurch die beschriebene Gl�attung erheblich besser funktioniert.

38

6 (M)HD-Gleichgewichte

Zum Test der korrekten Funktionsweise der Kopplung der beiden Teilgitter des Yin-Yang-Gitters sollen analytisch MHD-Gleichgewichte berechnet und anschlie�end in den Codeimplementiert werden, um zu testen, ob das Gleichgewischt stabil bleibt. Wenn ein Sys-tem im Gleichgewicht ist, muss f�ur eine beliebige physikalische Gr�o�e A gelten:

@A

@t= 0

6.1 Gleichgewicht ohne Magnetfeld und mit nur radialerGeschwindigkeit

6.1.1 Annahmen

Bevor man das Gleichungssystem - bestehend aus Gleichungen (3.26), (3.30), (3.28) - ana-lytisch l�osen kann, sind noch weitere Annahmen notwendig. Nehmen wir zur Konstruktiondes ersten Gleichgewichts an, das kein Magnetfeld ~B existiere und die Geschwindigkeit ~vnur eine radiale Komponente vr = v0

pr besitze. Weiterhin sei die Temperatur konstant,

also P = T0� (vgl. Gleichung (3.15)), es handelt sich also um ein isothermes Gas mit = 1. Zusammengefasst bedeutet dies:

~B = ~0

vr = v0pr

v# = 0

v# = 0

T = T0 ) P = T0� = u

Aufgrund dieser Annahmen vereinfacht sich das Gleichungssystem und reduziert sich auflediglich zwei:

~r � (�~v) = 0 (6.1)

~r � (� [~v~v]) + ~rP + g�r = 0 ; (6.2)

wobei r erneut f�ur den radialen Einheitsvektor in sph�arischen Polarkoordinaten steht.Im Folgenden noch zu bestimmen sind die Dichte � und die Gravitationsbeschleunigungg, die hier nicht vorgegeben, sondern so berechnet werden, dass sich eine analytischesonnenwind�ahnliche L�osung ergibt.

39


6.1.2 L�osung des Gleichgewichts f�ur sph�arische Symmetrie

In sph�arischer Symmetrie und unter der Annahme, dass sowohl die Geschwindigkeit ~v nureine r-Komponente vr besitzt als auch die Dichte � nur eine Funktion von r ist (� = �(r)),haben Gleichungen (6.1) und (6.2) folgende Gestalt:

1

r2@

@r

�r2�vr

�= 0 (6.3)

1

r2@

@r

�r2�v2r

�+ T0

@�

@r+ g� = 0 (6.4)

Aus Gleichung (6.3) ergibt sich mit einer beliebigen Konstanten K direkt:

r2�v0pr = K (6.5)

�(r) =K

v0r5=2� �0

r5=2(6.6)

Mit der gew�ahlten Form f�ur die Geschwindigkeit vr und der gem�a� (6.6) berechnetenDichte � wird Gleichung (6.3) erf�ullt. Einsetzen in Gleichung (6.4) ergibt:

1

r2@

@r

0BB@r2 �0r5=2

v20r| {z }=�0v20r

1=2

1CCA� 5

2T0

�0r7=2

+ g�0r5=2

= 0 (6.7)

1

2r2�0v

20r�1=2 � 5

2T0

�0r7=2

+ g�0r5=2

= 0 (6.8)

1

2v20 �

5

2T0r

�1 + g = 0 (6.9)

Daraus ergibt sich f�ur die Gravitationsbeschleunigung g(r):

g(r) =5T0r

�1 � v202

(6.10)

Analytisch ergibt sich also mit der Gravitationbeschleunigung aus Gleichung (6.10) einGleichgewicht. Dieses wurde in den Code implementiert, um zu pr�ufen, ob das Gleichge-wicht stabil bleibt. In Abbildung 6.1 sind die Startl�osungen der Geschwindigkeitskompo-nente vr sowie der Dichte � dargestellt.

Dargestellt sind die Di�erenzen von � und vr zwischen den Werten bei t = 15 und derStartl�osung

Es ist noch anzumerken, dass in den Abbildungen erneut n = �m

gilt, mit m = 1 folgtn = �. Man erkennt in Abbildungen 6.1 und 6.2 wie sich die Geschwindigkeitskomponentevr und die Dichte � mit dem Radius r �andern. Insgesamt sind die Di�erenzen mit ma-ximal etwa 10�8 im Bereich der numerischen Zahlendarstellung, so dass sich festhaltenl�asst, dass das analytische Gleichgewicht auch numerisch bestehen bleibt.

40

6.1 Gleichgewicht ohne Magnetfeld und mit nur radialer Geschwindigkeit

Abbildung 6.1: a) Die Geschwindigkeitskomponente vr als Funktion vom Radius r (links), b)Die Dichte � als Funktion vom Radius r (rechts)

Abbildung 6.2: a) Di�erenz der Geschwindigkeitskomponete vr zwischen dem Zeitpunkt t = 15und der Startl�osung , b) Die Dichte � in der Darstellung von Teil a)

6.1.3 St�orung des Gleichgewichtes

Um die f�ur die Simulationen verwendeten Transformationen zwischen den Teilgittern zutesten, sollen die Dichte � und die Geschwindigkeitskomponente vr des in Abschnitt 6.1beschriebenen Gleichgewichtes nun gest�ort werden. Au�erdem soll die Propagation derSt�orung beobachtet werden. Die Richtigkeit der Transformationen und die Qualit�at derKopplung soll daran gepr�uft werden, ob die St�orungen - zum Startpunkt an verschiedeneStellen des Gitters platziert - im zeitlichen Verlauf der Simulation vergleichbare Ergeb-nisse liefern. Zuerst wird die St�orung in die �Aquatorebene platziert, welche sich genaumittig in dem Teilgitter Yin be�ndet. Als weitere Stellen werden Randpunkte zwischenden beiden Teilgittern und die Pole dienen. Sollte die Simulation an den verschiedenenStellen vergleichbare Ergebnisse liefern, ist der Test der Transformationen erfolgreich, da

41


die Entwicklung der St�orung aufgrund der sph�arischen Symmetrie nicht von der Stelleabh�angen sollte, wo sie zum Startpunkt platziert wurde. Als letztes soll die St�orung nochan einem der Pole platziert werden, da die Pole die Bereiche f�ur sph�arische Simulationensind, welche man mit Hilfe von den in Kapitel 1 angesprochenen herk�ommlichen Gitternnicht untersuchen kann. Im Folgenden soll die Form der gew�ahlten St�orung beschriebenwerden.Die St�orung bei der Wahl von sph�arischen Polarkoordinaten ist in alle drei Raumrichtun-gen gleich, deswegen beschr�anken wir uns bei der Beschreibung auf die #-Richtung, dieweiteren zwei Raumrichtungen erfolgen analog. Wir beschr�anken uns bei der Beschreibungauf die St�orung der Dichte �, da die der Geschwindigkeitskomponente vr analog gest�ortwird. Die St�orung der Dichte �s ist wie folgt gew�ahlt worden:

�s =

(sin2

��2#�(#m�5�#)

5�#

�; f�ur#m � 5�# � # � #m + 5�#

0 sonst(6.11)

wobei �# der Gitterabstand und #m der Ort sei, an welchem die St�orung ihren Maximal-wert besitzt. Die Breite dieser St�orung soll in alle drei Raumrichtungen gleich gew�ahltwerden.Durch die Wahl dieser Form der St�orung wird erreicht, dass die St�orung �s an der Stelle#m = 1 ihren Maximalwert besitzt und zu jeder Seite ab einem gewissen Abstand (nach5�#) auf Null abgefallen ist. Die Amplitude von �s wurde so gew�ahlt, dass die St�orunggro� genug ist, um sie in den Abbildungen erkennen zu k�onnen und von vr in der Weise,dass die Geschwindigkeit der St�orung gro� genug ist, dass diese nicht zu schnell zerf�alltund weiter nach au�en propagiert, um eine Entwicklung beobachten zu k�onnen. Es hatsich als sinnvoll herausgestellt, in einem Gebiet (11 � r � 21, in normierten Einheiten) zusimulieren und rm = 15 zu w�ahlen, da f�ur gr�o�er werdende Radien r der Dichtegradientkleiner wird (�(r) / r�5=2). In Abbildung 6.3 ist die Entwicklung der Dichte � am Beispielvon �s zu erkennen.In dieser Abbildung ist zun�achst der Verlauf der homogenen Hintergrunddichte in Ab-h�angigkeit von Radius r dargestellt. �s bewegt sich mit der Zeit nach au�en und wirddabei kleiner. Zum Zeitpunkt t = 10 ist �s zu klein geworden, um diese noch in der Ab-bildung 6.3 zu erkennen. Abbildung 6.4 zeigt vier mal die Dichte zum Zeitpunkt t = 10an den beschriebenen verschiedenen Orten. Alle Abbildungen zeigen prinzipiell die glei-chen Ergebnisse, so dass sich festhalten l�asst, dass eine zufrieden stellende Kopplung derbeiden Teilgitter gew�ahrleistet ist. Es bleibt allerdings anzumerken, dass unterschiedli-che Amplituden an den betrachteten Stellen gew�ahlt werden mussten, so dass dort nochVerbesserungspotential bleibt. Trotz dieser Anmerkung ist festzustellen, dass der Testgelungen ist, da prinzipiell die Entwicklung der St�orung an allen Stellen gleich ist, derUnterschied bei den Ergebnissen liegt haupts�achlich in der Amplitude.Abbildung 6.5 ist analog zu Abbildung 6.4, es wird allerdings die radiale Geschwindig-keitskomponente vr gezeigt.Die dunklen Bereiche zwischen der St�orung und dem Mittelpunkt der Kugel in Abbil-dung 6.5 bedeuten, dass dort die Geschwindigkeit kleiner wird, was mit einer leichtenDichteerh�ohung konsistent ist In Abbildung 6.6 erkennt man den Verlauf der Dichte zueinem sp�ateren Zeitpunkt.

42


Abbildung 6.3: Cartesische Darstellung der Dichte � mit der St�orung �s am Nordpol zu ver-schiedenen Zeiten (t=2,5 (oben links), t=5 (unten links), t=7,5 (oben rechts),t=10 (unten rechts))

Abbildung 6.4: Darstellung von � (�s in der �Aquatorebene (oben links), am Rand # = �=4(unten links), am Rand # = 3�=4 (oben rechts), am Nordpol (unten rechts))

43


Abbildung 6.5: Darstellung analog zu Abbildung 6.4 hier mit vr

Es wird deutlich, dass �s kleiner geworden ist im Vergleich zu Abbildung 6.4 und sich

Abbildung 6.6: Vgl. Abbildung 6.4 zu einem sp�ateren Zeitpunkt

weiter zum Rand des Simulationsgebietes bewegt hat. Abbildung 6.7 stellt erneut vr dar

44


(vgl. 6.5). Auch hier ist zu erkennen, dass die St�orung sich weiter in Richtung des Randes

Abbildung 6.7: Vgl. Abbildung 6.5 zu einem sp�ateren Zeitpunkt

bewegt hat. Die dunklen Bereiche sind verschwunden, da auch �s klein geworden ist. Ins-gesamt stellt man fest, dass trotz der Modi�kationen der Amplitude, die Tests bzgl. derKopplung als erfolgreich bezeichnet werden k�onnen.

45


6.2 Beispiel f�ur einen axialsymmetrischen Sonnenwind

In diesem Abschnitt soll als weiteres Beispiel ein MHD-Gleichgewicht analytisch unter-sucht werden, welchem die Annahme eines globalen - zum �Aquator symmetrischen, teil-weise o�enen - Magnetfeldes ~B zugrunde liegt, welches der Arbeit von Tsinganos & Low(1989) entnommen ist. Diese Kon�guration repr�asentiert die Situation eines Sonnenwin-des, welcher in den magnetisch o�enen Regionen zentriert um die Pole entweicht, dasPlasma bleibt in einem magnetisch geschlossenem G�urtel um den �Aquator eingeschlos-sen. Die Geschwindigkeit ~v wird dabei als parallel zum Magnetfeld ~B angenommen. Diezu erf�ullenden MHD Gleichungen f�ur ein Gleichgewicht sind (vgl. Gleichungen (3.26) -(3.30)):

~r � (�~v) = 0 (6.12)

1

�0(~r� ~B)� ~B � ~rP � �GM

r2r = �(~v � ~r)~v (6.13)

~r� (~v � ~B) = 0 (6.14)

1

� 1(~v � ~r)

�P

�

�+ P (~v � ~r)1

�= � (6.15)

~r � ~B = 0 (6.16)

G (g = GMr2) in Gleichung (6.13) stellt die Gravitationsbkonstante dar. M ist die Sonnen-

masse.

6.2.1 Annahmen f�ur die Gleichgewichtsbedingung

Folgende Gr�o�en erf�ullen Gleichungen (6.12) - (6.16) (siehe Tsinganos & Low (1989):

~B(r) =1

r sin#

�1

r

@A

@#r � @A

@r#

�; mitA = D(r) sin2 # (6.17)

P (r) = P0(r) + P1(r) sin2 # (6.18)

�(r) =��0(r) + �1(r) sin

2 #�cos# (6.19)

W (r) =�

�(r)(6.20)

~v(r) = W (r) ~B(r) ; (6.21)

wobei die f�unf Gr�o�en P0, P1,W , �0, �1 Funktionen des Radius r sind und � ein konstanterParameter ist. Die Gr�o�e � beschreibt ein externe Energiequelle (Heizung bzw. K�uhlung).Die Auswertung von Gleichung (6.13) f�uhrt zun�achst auf zwei voneinander unabh�angigeGleichungen f�ur die beiden Koordinaten r und #. Es l�asst sich allerdings eine der beidenin zwei Gleichungen separieren, da sie sowohl aus von # abh�angigen als auch von # un-abh�angigen Summanden besteht, die unabh�angig voneinander verschwinden m�ussen. Esergeben sich also drei Gleichungen f�ur die vier Unbekannten P0, P1,W und D, so dass sich

46


P0, P1, W abh�angig von D ausdr�ucken lassen. Gleichung (6.15) l�asst sich dazu nutzen,die Energiequelle � festzulegen.

dP0

dr= �4�

r2D

d

dr

�WD

r2

�� GM

Wr2(6.22)

dP1

dr=

4�

r2D

d

dr

�WD

r2

�� 2�WD

r4dD

dr+

�

r3W

�dD

dr

�2

� 1

�0r2L(D)

dD

dr(6.23)

P1 =�

r2D

d

dr

�W

dD

dr

��

2r2W

�dD

dr

�2

� 1

�0r2L(D)D (6.24)

�0 =2WD

�r2

�1

� 1

d

dr(P0W ) + P0

dW

dr

�(6.25)

�1 =2WD

�r2

�1

� 1

d

dr(P1W ) + P1

dW

dr

�� 2W 2P1

�( � 1)r2dD

dr; (6.26)

wobei der Operator L eingef�uhrt wurde:

L =d2

dr2� 2

r2; (6.27)

welcher auf folgender Weise mit der elektrischen Stromdichte ~J zusammenh�angt:

~J = � 1

�0rL(D) sin# ' (6.28)

In dieser Formulierung bleibt die Funktion D eine frei w�ahlbare Funktion, die die Formdes Magnetfeldes ~B bestimmt. Setzt man die Funktion D zu

D = D4r4 +D2r

2 +D1

r(6.29)

im Bereich r < r� und es gilt

D = D� (6.30)

f�ur r > r�, wobei D1, D2, D4, D� und r� Konstanten sind und r� steht die so genann-te

"source surface\ bezeichnet (vgl. Priest (1982)), es ist der Abstand r von der Sonne,

ab welchem das Magnetfeld nur noch eine radiale Komponente besitzt. Setzt man nunGleichung (6.29) in Gleichung (6.28), so ergibt sich, dass die Stromdichte ~J nur von D4

abh�angt:

~J = � 1

�0r

�d2

dr2� 2

r2

��D4r

4 +D2r2 +

D1

r

�sin#'

= �10D4

�0r sin#' (6.31)

47


Einsetzen von Gleichung (6.29) in Gleichung (6.17) liefert die Form des Magnetfeldes imBereich r < r�:

~B = ~B1 + ~B2 + ~B4 (6.32)

~B1 =D1

r3

�2 cos# r + sin# #

�(6.33)

~B2 = 2D2

�cos# r � sin# #

�(6.34)

~B4 = 2D4r2�cos# r � 2 sin# #

�(6.35)

Das Magnetfeld ist eine Superposition eines Dipol-Feldes ~B1, eines homogenen Feldes ~B2

in z-Richtung (vgl. Abbildung 6.8) und eines nicht-stromfreien (vgl. Gleichung (6.31)) Ma-

gnetfeldes ~B4 , deren Amplituden durch die Parameter D1, D2, D4 gegeben sind. W�ahltman nun D so, dass die Ableitung der Funktion D an der Grenze r = r� verschwindet

4D4r5�+ 2D2r

3��D1 = 0 ; (6.36)

so ist das Magnetfeld an der Stelle r = r� stetig di�erenzierbar (vgl. Gleichung (6.17)).F�ur r > r� gilt f�ur das Magnetfeld nach Gleichung (6.30):

~B =2D� cos#

r2r ; (6.37)

wobei der Parameter D� durch

D� = D(r�) (6.38)

gegeben ist. Mit der hier dargestellten Form f�ur D k�onnen nun die �ubrigen radialenFunktionen P0, P1, W , �0, �1 anhand folgender Gleichungen ((6.22)-(6.26)) berechnetwerden.

6.2.2 L�osung der Gleichgewichtsbedingung

Bestimmung der Funktion W

Aus beiden Gleichungen (6.24) und (6.23) l�asst sich P1 eliminieren, kombiniert man an-schlie�end beide Gleichungen, erh�alt man eine Gleichung f�ur W . Diese Gleichung kannanalytisch integriert werden (Tsinganos & Low 1989). F�ur r < r� erh�alt man f�ur W :

W = W0

Z r

0

�5=2 exp[�20D4I(�)]

(D4�5 +D2�3 +D1)1=2d� +W1 ; (6.39)

mit

I(�) =

Z �

0

� 04 d� 0

4D4� 05 + 2D2� 03 �D1; (6.40)

48


wobei W0 und W1 Integrationskonstanten sind, � ist in beiden Gleichungen die Integrati-onsvariable. F�ur r > r� ergibt sich:

W =1

�

�Cr2 +

1

�0

�; (6.41)

dabei ist C eine weitere Integrationskonstante. Um nun die Gr�o�en aus Gleichungen (6.22)bis (6.26) f�ur ein Beispiel explizit berechnen zu k�onnen, muss zun�achst das Integral f�urW (Gleichung (6.39)) gel�ost werden.Mit der Annahme D4 = 0, d. h. f�ur ein stromfreies Magnetfeld, l�asst sich das Integral f�urW analytisch l�osen. Es vereinfacht sich zu:

W = W0

Z r

0

�5=2

(D2�3 +D1)1=2d� +W1 (6.42)

Mit der Substitution � = ry ergibt sich:

� = ry ) d� = r dy (6.43)

W = W0

Z r

0

�5=2

(D2�3 +D1)1=2d� +W1 (6.44)

= W0

Z 1

0

(ry)5=2r

(D2(ry)3 +D1)1=2dy +W1 (6.45)

= W0r7=2| {z }

:=W 0

0

Z 1

0

y5=2

(D2r3| {z }

:=D0

2

y3 +D1)1=2dy +W1 (6.46)

= W 0

0

Z 1

0

y5=2

(D0

2y3 +D1)1=2

dy +W1 (6.47)

Eine weitere Substitution y = z1=3 ergibt:

y = z1=3 ) dy =1

3z�2=3 dz (6.48)

W = W 0

0

Z 1

0

y5=2

(D0

2y3 +D1)1=2

dy +W1 (6.49)

=1

3W 0

0

Z 1

0

�z1=3

�5=2z�2=3

(D0

2z +D1)1=2dz +W1 (6.50)

=1

3W 0

0

Z 1

0

z1=6

(D0

2z +D1)1=2dz +W1 (6.51)

=1

3W 0

0

1

D1=21

Z 1

0

z1=6�1 +

D0

2

D1z�1=2 dz +W1 (6.52)

49


Mit der Substitution z = D1

D0

2

v folgt:

W =1

3W 0

0

1

D1=21

Z 1

0

z1=6�1 +

D0

2

D1z�1=2 dz +W1 ; z =

D1

D0

2

v ) dz =D1

D0

2

dv (6.53)

=1

3W 0

0

1

D1=21

Z D0

2

D1

0

�D1

D0

2

�1=6v1=6D1

D0

2

(1 + v)1=2dv +W1 (6.54)

=1

3W 0

0

D2=31

D07=62

Z D0

2

D1

0

v1=6

(1 + v)1=2dv +W1 (6.55)

F�ur die Integrationsgrenze giltD0

2

D1< 1

2D2r3�| {z }

:=D00

2

�D1 = 0) 2D00

2 = D1 ) 2D00

2

D1= 1 ) D00

2

D1< 1 (6.56)

D0

2 = D2r3 = D2r �3

�r

r�

�3

= D00

2

�r

r�

�3

) D0

2

D1=D00

2

D1

�r

r�

�3

< 1 ; (6.57)

da die L�osung f�urW f�ur den Fall r < r� berechnet wird. Die Form der hypergeometrischenReihe lautet (siehe Storch & Wiebe (2003) 13.C.13 Beispiel):

1Xn=0

(a)nxn

n!=

1

(1� x)a; f�ur alle jxj < 1 ; (6.58)

mit dem Pochhammer-Symbol (a)n:

(a)n = a � (a+ 1) � (a+ 2) � : : : (a+ n� 1) (6.59)

An dieser Stelle kann der Integrand in eine hypergeometrischen Reihe �uberf�uhrt werden(vgl. Gleichung (6.58), und es folgt:

50


W =1

3W 0

0

D2=31

D07=62

Z D0

2

D1

0

v1=6

(1 + v)1=2dv +W1 (6.60)

=1

3W 0

0

D2=31

D07=62

Z D0

2

D1

0

v1=61Xn=0

�1

2

�n

(�v)nn!

dv +W1 (6.61)

=1

3W 0

0

D2=31

D07=62

1Xn=0

�1

2

�n

(�1)nn!

Z D0

2

D1

0

v1=6+n dv +W1 (6.62)

=1

3W 0

0

D2=31

D07=62

1Xn=0

�1

2

�n

(�1)nn!

�vn+7=6

n+ 7=6

�D0

2

D1

0

+W1 (6.63)

=1

3W 0

0

D2=31

D07=62

1Xn=0

�1

2

�n

(�1)nn!

�D0

2

D1

�n+7=6

n+ 7=6+W1 (6.64)

=1

3W 0

0

1

D1=21

1Xn=0

�1

2

�n

1

n!

��D0

2

D1

�nn+ 7=6

+W1 (6.65)

Da fogender Zusammenhang gilt (vgl. beispielsweise Abramovitz & Stegun (1974))

(z)n =�(z + n)

�(z); (6.66)

wobei �(z) die Gamma-Funktion ist, ergibt sich mit

�(z + 1) = z�(z) ; (6.67)

folgende Gleichung:

�76

�n�

136

�n

=�(7

6+ n)

�(76)

�(136)

�(136+ n)

(6.68)

=�(7

6+ n)

�(76)

76�(7

6)

(7=6 + n)�(76+ n)

(6.69)

=76

7=6 + n(6.70)

Daher folgt:

1

7=6 + n=

�76

�n�

136

�n

� 67; (6.71)

51


Aus dieser Identit�at folgt schlie�lich Gleichung (6.71). Damit folgt also insgesamt f�ur W :

W =2

7W 0

0

1

D1=21

1Xn=0

�12

�n

�76

�n�

136

�n

� 1n!

��D2r

3

D1

�n

+W1 (6.72)

Mit der verallgemeinerten hypergeometrischen Reihe (vgl. Storch & Wiebe (2003) 13.C.13Beispiel)

1Xn=0

(a1)n(a2)n : : : (aq)n(b1)n(b2)n : : : (bp)n

� xn

n!=: qFp(a1; : : : ; aq; b1; : : : ; bp;x) (6.73)

ergibt sich f�ur W :

W =2

7

W0

D1=21

r7=2 2F1

�1

2;7

6;13

6;�D2r

3

D1

�+W1 (6.74)

Um nun die weiteren radialen Funktionen P0, P1, �0, �1 berechnen zu k�onnen, betrachtenwir zwei F�alle:

� limD1!0 2F1

�12; 76; 136;�D2r3

D1

�� limD2!0 2F1

�12; 76; 136;�D2r3

D1

�Der Fall D2 ! 0 ergibt sich direkt aus Gleichung (6.74)

W =7

2

W0

D1=21

r7=2 +W1 ; (6.75)

da von 2F1

�12; 76; 136; 0�nur der Term zu n = 0 �ubrig bleibt und dieser 1 liefert (vgl. Glei-

chung (6.73)). Der erste der beiden F�alle ist schwieriger auszuwerten. Dazu benutzen wirfolgende Identit�at (vgl. Abramovitz & Stegun (1974), Formel (15.3.6)):

2F1(a; b; c; z) =�(c)�(b� a)

�(b)�(c� a)(�z)�a 2F1

�a; 1� c+ a; 1� b+ a;

1

z

�+�(c)�(a� b)

�(a)�(c� b)(�z)�b 2F1

�b; 1� c+ b; 1� a+ b;

1

z

�(6.76)

Im hier vorliegenden Fall ergibt sich:

2F1

�1

2;7

6;13

6;�D2

D1r3�

=�(13

6)�(2

3)

�(76)�(5

3)

�D2

D1r3��

1

2

+�(13

6)�(�2

3)

�(12)�(1)

�D2

D1r3��

7

6

(6.77)

52


Daraus folgt:

2F1

�1

2;7

6;13

6;�D2

D1r3�

=�(13

6)�(2

3)

�(76)�(5

3)| {z }

= 7

4

�D2

D1

��

1

2

r�3=2

+�(13

6)�(�2

3)

�(12)�(1)

�D2

D1

��

7

6

r�7=2 (6.78)

Somit ergibt sich f�ur den Fall D1 ! 0 f�ur W durch Einsetzen in Gleichung (6.74):

W =2

7

W0

D1=21

r7=27

4

�D1

D2

�1=2

r�3=2 (6.79)

+2

7

W0D2=31

D7=62

�(136)�(�2

3)

�(12)�(1)| {z }

!0

+W1 (6.80)

=1

2

W0

D1=22

r2 +W1 (6.81)

Die Gleichungen (6.75) und (6.79) ergeben sich auch, wenn man direkt in Gleichung (6.42)D1 bzw. D2 zu Null setzt, was das erhaltene Ergebnis ohne Vereinfachung best�atigt.

Bestimmung der restlichen radialen Funktionen aus dem W

Im Folgenden sollen die beiden F�alle D1 = D4 = 0 (homogenes Magnetfeld) und D2 =D4 = 0, (Dipol-Feld).Mit Hilfe der beschriebenen F�alle sollen nun aus den daraus entstehenden zwei unter-schiedlichen Funktionen W die �ubrigen Gr�o�en nach Gleichungen (6.22) bis (6.26) be-stimmt werden.

Berechnung f�ur D1 = D4 = 0

D = D2r2 (6.82)

W =1

2

W0

D1=22

r2 +W1 (6.83)

F�ur die weitere Berechnung soll der Einfachheit halber W1 = 0 gelten, d. h.:

W =1

2

W0

D1=22

r2 (6.84)

53


Durch Einsetzen von W und D in Gleichung (6.24) �ndet man f�ur den Druck P1:

P1 = 2�D3=22 W0r

2 (6.85)

Setzt man nun W und D in Gleichung (6.22) ein, ergibt sich:

dP0

dr= �4�D3=2

2 W0r � 2�GMD1=22

W0r4(6.86)

Integration �uber r ergibt f�ur P0:

P0 =

Z �4�D3=2

2 W0r0 � 2�GMD

1=22

W0r04

!dr0 (6.87)

= �2�D3=22 W0r

2 +2�GMD

1=22

3W0r3+ P0 ; (6.88)

wobei P0 eine Integrationskonstante ist. Einsetzen von D, W und P0 in Gleichung (6.25)ergibt nach l�angerer, aber einfacher Rechnung f�ur den ersten Term der Energiequelle �0:

�0 = �2 + 1

� 1D

3=22 W 3

0 r5 +

1

3

2 � 3

� 1GMD

1=22 W0 (6.89)

Schlie�lich ergibt f�ur den zweiten Term der Energiequelle �1:

�1 = 2

� 1D

3=22 W 3

0 r5 (6.90)

F�ur den Fall = 1 muss gelten �0 = �1 = 0 gelten, da dann ein isothermes Gas betrachtetwird. Zusammengefasst ergibt sich daher f�ur alle betrachteten physikalischen Gr�o�en f�urden Fall D1 = 0 und D4 = 0:

~B = 2D2


�~v = W0D

1=22 r2


�P = �2�D3=2

2 W0r2 +

2�GMD1=22

3W0r3+ 2�D

3=22 W0r

2 sin2 #+ P0

= 2�D3=22 W0

�GM

3W 20D2r3

� r2 cos2 #

�� =

��2 + 1

� 1D

3=22 W 3

0 r5 +

1

3

2 � 3

� 1GMD

1=22 W0 + 2

� 1D

3=22 W 3

0 r5 sin2 #

�cos#

=2

� 1cos#D

3=22 W 3

0 r5

�1

6(2 � 3)

GM

D2W 20

� cos2 #� 1

�� =

�

W= 2�

D1=22

W0r�2

54


Berechnung f�ur D2 = D2 = 0

D =D1

r(6.91)

W =2

7

W0

D1=21

r7=2 +W1 (6.92)

Erneut wird W1 = 0 angenommen, d. h.:

W =2

7

W0

D1=21

r7=2 (6.93)

Durch Einsetzen von W und D in Gleichung (6.24) �ndet man f�ur P1:

P1 = �4

7�D

3=21 W0r

�5=2 (6.94)

Setzt man nun W und D in Gleichung (6.22) ein, ergibt sich:

dP0

dr= �4

7�D

3=21 W0r

�7=2 � 7�GMD1=21

2W0r11=2(6.95)

Integration �uber r ergibt f�ur P0:

P0 =

Z �4

7�D

3=21 W0r

0�7=2 � 7�GMD1=21

2W0r011=2

!dr0 (6.96)

=8

35�D

3=21 W0r

�5=2 +7�GMD

1=21

9W0r9=2+ P 0

0 ; (6.97)

wobei P 00 eine weitere Integrationskonstante ist. Nach Einsetzen von D, W und P0 in

Gleichung (6.25) erh�alt man f�ur �0:

�0 =1

1715

224 � 160

� 1D

3=21 W 3

0 r1=2 +

1

63

28 � 36

� 1GMD

1=21 W0r

�3=2 (6.98)

Schlie�lich ergibt sich f�ur �1:

�1 = � 1

343

112 � 48

� 1D

3=21 W 3

0 r1=2 (6.99)

Zusammengefasst ergibt sich daher f�ur alle betrachteten physikalischen Gr�o�en f�ur denFall D1 = 0 und D4 = 0:

55


~B = D1

�2 cos#

r3r +

sin#

r3#

�~v =

7

2W0D

1=21 r1=2

�2 cos# r + sin# #

�P =

8

35�D

3=21 W0r

�5=2 +7�GMD

1=21

9W0r9=2� 4

7�D

3=21 W0r

�5=2 sin2 #

= �D1=21

�8

35D1W0r

�5=2 � 7GM

9W0r9=2� 4

7D1W0r

�5=2 sin2 #

�� =

1

1715

224 � 160

� 1D

3=21 W 3

0 r1=2 cos#+

1

63

28 � 36

� 1GMD

1=21 W0r

�3=2 cos#

� 1

343

112 � 48

� 1D

3=21 W 3

0 r1=2 sin2 # cos#

� =�

W=

7

2�D

1=21

W0r�7=2

6.2.3 Das Gleichgewicht in den Koordinaten von Yin und Yang

Nun m�ussen alle Gr�o�en von sph�arischen Polarkoordinaten in Koordinaten auf Yin undYang transformiert werden.In diesem Abschnitt steht der �Ubersicht halber der Index i f�ur eine Gr�o�e auf Yin undder Index a f�ur eine Gr�o�e auf Yang.Die Gr�o�en P , �, � und W in Koordinaten von Yin lauten:

Pi = P0(r) + P1(r) sin2 #i (6.100)

�i = �i(r) (6.101)

�i =��0(r) + �1(r) sin

2 #i�cos#i (6.102)

Wi =�

�i(6.103)

Das Magnetfeld ~B und die Geschwindigkeit ~v f�ur den Bereich r < r� in den Koordinatenvon Yin lauten:

~Bi = ~B1i + ~B2i + ~B4i (6.104)

~B1i =D1

r3

�2 cos#i ri + sin#i #i

�(6.105)

~B2i = 2D2

�cos#i ri � sin#i #i

�(6.106)

~B4i = 2D4r2�cos#i ri � 2 sin#i #i

�(6.107)

~vi = W ~Bi (6.108)

56


F�ur den Bereich r > r� lauten das Magnetfeld ~B und die Geschwindigkeit ~v in den Koor-dinaten von Yin (vgl. 6.41):

~Bi =2D� cos#i

r2ri (6.109)

~vi =1

�

�Cr2 +

1

�0

�~Bi (6.110)

Durch Einsetzen der Gleichungen (2.13), (2.14), (2.15) sowie Gleichungen (2.17), (2.18)in alle berechneten Gr�o�en auf Yin liefert die Gr�o�en auf Yang:

Pa = P0(r) + P1(r)(1� sin2 #a sin2 'a) (6.111)

�a = �a(r) (6.112)

�a =��0(r) + �1(r)(1� sin2 #a sin

2 'a)�sin#a sin'a (6.113)

Wa =�

�a(6.114)

Das Magnetfeld ~B und die Geschwindigkeit ~v f�ur den Bereich r < r� in den Koordinatenvon Yang lauten:

~Ba = ~B1a + ~B2a + ~B4a (6.115)

~B1a = D1

�2 sin#a sin'a

r3ra � cos#a sin'a

r3#a � cos'a

r3'a

�(6.116)

~B2a = 2D2

�sin#a sin'a ra + cos#a sin'a #a + cos'a 'a

�(6.117)

~B4a = 2D4r2�sin#a sin'a ra + 2 cos#a sin'a #a + 2 cos'a 'a

�(6.118)

~va = W ~Ba (6.119)

F�ur den Bereich r > r� lauten das Magnetfeld ~B und die Geschwindigkeit ~v in den Koor-dinaten von Yang:

~Ba =2D� sin#a sin'a

r2ra (6.120)

~va = W ~Ba ;mitW =1

�

�Cr2 +

1

�0

�~Ba (6.121)

Das Magnetfeld ~B soll zum Test der Transformation in den Code implementiert werdenum das Magnetfeld in Koordinaten von Yin und Yang zu zeichnen. Zun�achst werden dieeinzelnen Bestandteile ~B1, ~B2, ~B4 dargestellt. Zu diesem Zweck darf dann lediglich derBereich innerhalb der

"source surface\ betrachtet werden, da Bedingung (6.36) nicht mehr

erf�ullt werden kann. Abbildung 6.8 stellt ~B1, ~B2 und ~B4 dar.

57


Abbildung 6.8: a) Das Dipol-Magnetfeld ~B1 (links), b) Das homogene Magnetfeld ~B2 (Mitte),c): Das nicht-stromfreie Magnetfeld ~B4 (rechts)

Die gestrichelte Linie gibt den Rand des Rechengebietes an, d. h. simuliert wurde bis zueinem Radius von r = 6. Es ist zu erkennen, dass das Magnetfeld ~B1 die Form eines Dipol-feldes hat, ~B2 ist ein homogenes Magnetfeld und zeigt in z-Richtung. Ist das Magnetfeld~B4 ungleich Null, ist das betrachtete System nicht mehr stromfrei (vgl. auch Gleichung(6.31)).In Abbildung 6.9 wurden die Parameter D1, D2, D4 sowie r� folgenderma�en gew�ahlt(siehe Tsinganos & Low (1989)):

D1 =4

5r3�; D2 = 1 ; D4 = � 3

10

1

r2�

; mit r� = 4 (6.122)

Daraus ergibt sich f�ur D(r�) = D�:

D� = D4r4�+D2r

2�+D1

r�=

3

2r2�

(6.123)

In Abbildung 6.9 ist das in Yin- und Yang-Koordinaten transformierte Magnetfeld ~Bdargestellt, welches sich als �Uberlagerung (gem�a� Gleichung (6.122)) der in Abbildung6.8 dargestellten Summanden ergibt.Vergleicht man Abbildung 6.9 mit der entsprechenden Abbildung aus Tsinganos & Low(1989), welches mit Hilfe sph�arischer Polarkoordinaten berechnet wurde, ist zu erkennen,

dass die Transformation des Magnetfeldes ~B aus Abschnitt 6.2.3 in Koordinaten von Yinund Yang die richtigen Ergebnisse liefert, da beide Abbildungen das selbe Resultat zeigen(vgl. Abbildung 6.10).

58


Abbildung 6.9: Das Magnetfeld ~B aus Gleichungen (6.104) bzw. (6.109) und (6.115) bzw. (6.120)

Abbildung 6.10: Das Magnetfeld ~B in sph�arischen Polarkoordinaten aus Tsinganos & Low (1989)

Die Transformationsformeln aus Kapitel 2.2 liefern die richtigen Ergebnisse.

59


60

7 Zusammenfassung und Ausblick

Der Schwerpunkt dieser Arbeit lag in der Untersuchung der N�utzlichkeit des von Ka-geyama & Sato (2004) eingef�uhrten �uberlappenden sph�arischen Yin-Yang-Gitters, sowieinsbesondere die Analyse einer e�ektiven Kopplung der beiden Teilgitter.Numerische Simulationen zu physikalischen Prozessen in sph�arischer Geometrie werdenklassisch mit Hilfe von sph�arischen Gittern realisiert. Daf�ur bietet es sich an, das in Kapitel1 beschriebene herk�ommliche Gitter zu verwenden. Der Umgang mit diesem klassischenGittersystem ben�otigt eine genaue Kenntnis der zu untersuchenden Gr�o�en in sph�arischenPolarkoordinaten. Dieses Koordinatensystem hat den Vorteil, dass der Umgang damit re-lativ vertraut ist und die Transformationsformeln in dieses System bekannt sind. EinNachteil dieses Systems ist allerdings, dass keine E�ekte in den polaren Regionen unter-sucht werden k�onnen, da es zu Divergenzen der Skalenfaktoren von Di�erentialoperatorenin diesem System kommt (vgl. Kageyama (2005)). Das von Kageyama & Sato (2004) vor-gestellte zusammengesetzte �uberlappende Yin-Yang-Gitter deutet auf eine L�osung diesesProblems durch die Verwendung zweier identischer, aber gegeneinander verdrehter Teil-gitter hin, welches in Kapitel 1 eingef�uhrt und in Kapitel 2 n�aher erl�autert wurde. DesWeiteren werden in diesem Kapitel die Transformationsformeln berechnet und die Dar-stellung der Teilgitter in der #-'-Ebene gezeigt. Der Vorteil der sinnvollen Beschreibungvon E�ekten in den polaren Regionen bringt allerdings den Nachteil mit sich, dass diebeiden Teilgitter aneinander gekoppelt werden m�ussen, denn durch die Zusammensetzungzweier Teilgitter entstehen nicht-physikalische R�ander im Rechengebiet. An ihren R�andern�uberlappen die beiden Teilgitter und es muss daf�ur gesorgt werden, dass sich in diesem Be-reich die Werte der betrachteten Gr�o�en auf den beiden Teilgittern nicht unterschiedlichentwickeln. Ein weiterer Vorteil dieses Gitters liegt darin, dass die beiden Teilgitter einen�Uberlapp bilden und man dadurch die Kopplung durch die in Kapitel 5 beschriebene unddiskutierte Interpolation gew�ahrleisten kann. Die zahlreichen Tests haben gezeigt, dassdie Interpolation nach Newton die genauesten Resultate liefert und deswegen auch als In-terpolationsmethode zu w�ahlen ist. Die von Laukert (2008) durchgef�uhrten Simulationenverwendeten das Yin-Yang-Gitter bereits, allerdings w�aren die Tests mit der optimiertenNewton-Interpolation zu wiederholen.Die korrekte Funktionsweise der Kopplung der Teilgitter wurde dann mit Hilfe der in Ka-pitel 6 beschriebenen Simulationsrechnungen getestet. Zun�achst stellte sich in Abschnitt6.2.3 heraus, dass die in Abschnitt 2.2 berechneten Transformationsformeln die richtigenResultate liefern. In Abschnitt 6.1 wurde analytisch ein Gleichgewicht berechnet und an-schlie�end in den Code implementiert, um festzustellen, ob das berechnete Gleichgewichtstabil bleibt.Nachdem sich ergab, dass dies der Fall ist, wurde dieses Gleichgewicht in Abschnitt 6.1.3gest�ort, um durch Platzieren der St�orung an verschiedenen Stellen - insbesondere in dem�Uberlappungsbereich der Teilgitter und an einem der beiden Pole der Kugel - zu testen,

61

7 Zusammenfassung und Ausblick

ob sich die St�orung wie theoretisch vorhergesagt verh�alt und ob sich zeigt, dass die Plat-zierung der St�orung in sph�arischer Symmetrie keinen Ein uss auf deren Verhalten hat.Dabei hat sich herausgestellt, dass die St�orung der Dichte � und der Geschwindigkeits-komponente vr tats�achlich zu prinzipiell gleichen Resultaten f�uhrt. Die lokale St�orung derDichte � bewegt sich mit der Zeit zum �au�eren Rand des Rechengebietes und wird dabeibreiter und acher, l�auft also wie erwartet auseinander.Dieses Verhalten zeigt sich wie vorhergesagt unabh�angig von der Position der St�orung.Zuerst hatten wir die St�orung auf die �Aquatorebene und anschlie�end an verschiedeneStellen auf dem �Uberlappungsbereich platziert, um zu testen, ob Rande�ekte zwischenden Teilgittern einen Ein uss auf das Ergenis haben. Es stellte sich heraus, dass dasSt�orungsverhalten unabh�angig von Rande�ekten ist.Im Anschluss daran wurde die St�orung am Nordpol der Kugel platziert, eine Region diemit einem herk�ommlichen Gitter nicht untersucht werden kann. Auch hier bleib die Formder St�orung erhalten, was zeigt, dass das Yin-Yang-Gitter prinzipiell - sogar in den pola-ren Regionen - gute Ergebnisse liefert und eine zufriedenstellende Kopplung der Teilgittergew�ahrleistet wird.An dieser Stelle zeigte sich allerdings auch noch ein Schwachpunkt dieses Gittermo-dells. Die St�orung schien sich zwar unabh�angig von der gew�ahlten Position zu verhalten.Die Amplitude der St�orung musste jedoch abh�angig von der Lage nachgeregelt werden.Vor allem die Amplitude der St�orung am Nordpol der Kugel nahm schneller ab als diein der �Aquatorebene. An dieser Stelle besteht daher noch Handlungsbedarf. Trotz derSchw�achen, die in Zukunft noch n�aher zu untersuchen sind, k�onnen die durchgef�uhrtenTests insgesamt als gelungen eingestuft werden.Somit gelangt man insgesamt zu dem Fazit, dass das Yin-Yang-Gitter eine gute Alternati-ve zu einem herk�ommlichen Gitter (vgl. Kapitel 2) darstellt, vor allem wenn physikalischeProzesse in polaren Regionen untersucht werden sollen. Auch wenn an einigen Stellennoch Schwachpunkte auftreten, die noch zu beheben sind.

62

Anhang A

Anhang

Die Vektoridentit�aten sind Huba (2006) entnommen.

A.1 Weitere M�oglichkeit zur Konstruktion von Gleichung (3.17)

Wenn die Entropie erhalten bleibt, kann Gleichung (3.17) auch aus deren Erhaltung ab-geleitet werden:Die Entropie ist ein Ma� f�ur die Unordnung in einem System und sie ist de�niert als (siehez. B. Rebhan (2005):

S = �cv ln�P

P0

��0�

� �+ eS0 (A.1)

Formt man die Entropie um, indem man die konstanten Terme P0 und �0 mit in die Kon-stante eS0 schreibt, ergibt sich:

S = �cv ln�P

�

�+ S0 (A.2)

Man erh�alt damit f�ur Gleichung (3.17):

dS

dt=@S

@t+ (~v � ~r)S = 0 (A.3)

Setzt man nun die De�nition f�ur S aus Gleichung (A.2) in Gleichung (A.3) ein, so ergibtsich f�ur den ersten Summanden:

@S

@t= �cv

�1

P

@P

@t�

1

�

@�

@t

�(A.4)

Der zweite Summand aus Gleichung (A.3) ergibt:

(~v � ~r)S = �cv~v ��1

P~rP �

1

�~r��

(A.5)

63

Anhang A Anhang

Durch Addition der beiden Terme und weitere Umformungen ergibt sich:

1

P

�@P

@t+ ~v � (~r � P )

�=

1

�

�@�

@t+ ~v � (~r�)

�

= 1

�

0BB@ @�

@t+ ~r � (�~v)| {z }

=0 (Gleichung (3.10))

��(~r � ~v)

1CCAWeiteres Umformen ergibt:

@P

@t+ ~v � (~rP ) + P (~r � ~v) = 0 (A.6)

@P

@t+ ~r � (~vP )� P (~r � ~v) + P (~r � ~v) = 0

Daraus ergibt sich schlie�lich Gleichung (3.17):

@P

@t+ ~r � (~vP ) + ( � 1)P (~r � ~v) = 0 (A.7)

A.2 Entwicklung von Gleichung (3.19)

Der Ausgangspunkt ist Gleichung (A.6):

@P

@t+ ~v � (~rP ) + P (~r � ~v) = 0

Der Druck P soll nun ersetzt werden �uber die skalare Gr�o�e u:

P = u

Eingesetzt ergibt:

u �1@u

@t+ u �1~v � (~ru) + u (~r � ~v) = 0

Zusammenfassen ergibt:

@u

@t+ ~v � (~ru) + u(~r � ~v) = 0

Man erh�alt schlie�lich Gleichung (3.19):

@u

@t+ ~r � (u~v) = 0 (A.8)

Gleichung (A.8) ist konservativ und wird in dieser Form im Code verwendet.

64

A.3 Konstruktion von Gleichung (3.25)

A.3 Konstruktion von Gleichung (3.25)

Ausgangspunkt hier ist Gleichung (3.24):

~E + ~v � ~B = � ~J

Eingesetzt in Gleichung (3.20):

@ ~B

@t= �(~r� ~E)

= ��~r�

�� ~J �

�~v � ~B

��Unter Ausnutzung von

~J =1

�0(~r� ~B)

ergibt sich:

@ ~B

@t= �

�~r�

�� ~J �

�~v � ~B

��= �

�~r� � ~J � �

�0

�~r�

�~r� ~B

�� ~r�

�~v � ~B

��= �

�~r� � ~J

��

�0

~r ~r � ~B| {z }=0

!!+

�

�0

�4 ~B

�+ ~r�

�~v � ~B

�(A.9)

Somit ergibt sich schlie�lich Gleichung (3.25):

@ ~B

@t=�~r�

�~v � ~B

��+

�

�04 ~B � ~r� � ~J

65

Anhang A Anhang

66

Literaturverzeichnis

Abramovitz, M. & Stegun, I. A. Handbook of Mathematical Functions. Dover Publications,New York, 1974

Bronstein, I. N., Semendjajew, K. A., Musiol, G. & M�uhlig, H. Taschenbuch der Mathe-

matik. Verlag Harry Deutsch, Thun und Frankfurt am Main, 2000

Chen, F. F. Introduction to Plasma Physics and Controlled Fusion. Springer Science +Business Media, LLC, New York, 1974

Chesshire, G. & Henshaw, W. D. Composite overlapping meshes for the solution of partialdi�erential equations. Journal of Computational Physics, 90, 1, 1990

Huba, J. D. NRL Plasma Formulary. Beam Physics Branch, Plasma Physics Division,Naval Research Laboratory, Washington, DC 20375, 2006

Kageyama, A. Dissection of a Sphere and Yin-Yang Grids. Journal of the Simulator, 3,20, 2005

Kageyama, A. & Sato, T. \Yin-Yang grid": An overset grid in spherical geometry. Geo-chemistry, Geophysics, Geosystems, 5, Q09005, 2004

Koldoba, A. V., Romanova, M. M., Ustyugova, G. V. & Lovelace, R. V. E. Three-dimensional Magnetohydrodynamic Simulations of Accretion to an Inclined Rotator:The \Cubed Sphere" Method. Astrophys. J. Lett., 576, L53, 2002

Laukert, M. MHD-Simulationen auf zusammengesetzten Gittern in sph�arischer Geome-

trie. Diplomarbeit, Ruhr-Universit�at Bochum, 2008

Neutsch, W. Koordinaten: Theorie und Anwendungen. Spektrum Akademischer VerlagGmbH, Heidelberg, 1995

Otto, A. Zur linearen und nichtlinearen Analyse resistiver Instabilit�a tsprozesse in schwach

zweidimensionalen Gleichgewichten. Dissertation, Ruhr-Universit�at Bochum, 1987

Priest, E. Solar Magnetohydrodynamics. D. Reidel Publishing Company, Dordrecht, 1982

Rebhan, E. Theoretische Physik 2. Elsevier GmbH, Spektrum Akademischer Verlag,M�unchen, 2005

Schwarz, H. R. Numerische Mathematik. B.G. Teubner, Stuttgart, 1986

67

Literaturverzeichnis

Storch, U. &Wiebe, H. Lehrbuch der Mathematik, Band 1, Analysis einer Ver�anderlichen.Spektrum Akademischer Verlag GmbH, Heidelberg, 2003

Tsinganos, K. & Low, B. C. Steady hydromagnetic ows in open magnetic �elds. II -Global ows with static zones. Astrophys. J., 342, 1028, 1989

68

Danksagung

An dieser Stelle m�ochte ich mich herzlichst bei allen bedanken, die mich bei der Erstellungdieser Arbeit unterst�utzt haben.In erster Linie bedanke ich mich bei PD Dr. Horst Fichtner, der mir durch die Themenstel-lung erst erm�oglicht hat, diese Arbeit zu schreiben. Er hat mir stets mit seinem Fachwissenkompetent zur Seite gestanden und war jederzeit o�en f�ur Gespr�ache und R�uckfragen.Vor allem m�ochte ich mich auch f�ur die geduldige und hilfsbereite Betreuung bei Dr. An-dreas Kopp bedanken, der mir nicht nur bei numerischen Problemen viele wertvolle Tipsgegeben hat. Aufgrund der gro�en Entfernung hat dies zahlreiche Wochenenden und denAustausch diverser E-Mails in Anspruch genommen.Bedanken m�ochte ich mich auch bei Christian R�ocken f�ur seine Freundschaft, st�andigeHilfsbereitschaft und die vielen konstruktiven Fachgespr�ache und Hinweise bei der Erstel-lung dieser Arbeit.Bei Prof. Dr. Reinhard Schlickeiser m�ochte ich mich f�ur die Aufnahme in einen kollegialen,von Freundlichkeit und Hilfsbereitschaft gepr�agten Lehrstuhl bedanken.Weiterhin bedanke ich mich bei Herrn Bernd Neubacher f�ur die schnelle Hilfe bei Computer-Problemen aller Art.Ein gro�er Dank geht an Hanna f�ur ihre Geduld und Unterst�utzung und ihr stets o�enesOhr. Sie war mir w�ahrend der ganzen Phase der Arbeit und dar�uber hinaus immer eingro�er R�uckhalt.Ein gro�er Dank geht nat�urlich auch an meine Eltern, die es mir �uberhaupt erst erm�oglichthaben, das Physikstudium sorgenfrei zu absolvieren und immer an mich geglaubt habenund mir immer wieder hilfsbereit zur Seite standen.

69

Optimierung der Verkn upfung uberlappender sph arischer ...Optimierung der Verkn upfung...

Documents

Transcript of Optimierung der Verkn upfung uberlappender sph arischer ...Optimierung der Verkn upfung...