802 - Simulation einer Pierce-Diode -...

802 – Simulation einer Pierce-DiodeVersuchsprotokoll zum F-Praktikum an der Ruhr-Universität Bochum

Tobias Krähling <[email protected]> 14.11.2009Version 1.0

Inhaltsverzeichnis

1 Einführung 1

2 Pierce-Diode 12.1 Plasmaoszillationsmodell von Tonks und

Langmuir . . . . . . . . . . . . . . . . . . . 22.2 Modifikationen im Pierce-Modell . . . . . 2

2.2.1 Dimensionslose Formulierung . . 32.3 Zu untersuchende Pierce-Dioden-

Modifikation . . . . . . . . . . . . . . . . . 3

3 Particle-in-Cell Simulation (PIC) 43.1 Implementierung . . . . . . . . . . . . . . 43.2 Tridiagonalbandlöser . . . . . . . . . . . . 6

3.2.1 Übertragung auf die Poisson-Gleichung . . . . . . . . . . . . . . 7

4 Programm 7

5 Simulationsergebnisse 95.1 Variation des freien Parameters α bei kon-

stanter Geschwindigkeit . . . . . . . . . . 95.2 Einfluss der Geschwindigkeitsverteilung 12

6 Zusammenfassung 13

Literatur 13

1 Einführung

Im Gegensatz zu interstellaren Plasmen, dieaufgrund ihrer Dimension über weite Bereicheals unbegrenzte Plasmen angesehen werdenkönnen, müssen bei technischen Plasmen dieendlichen Abmessungen berücksichtigt wer-den. Die Begrenzungen des Plasmas führendabei zu Veränderungen der Plasmadynamik.Ein einfaches thermionisches Entladungsmo-dell für begrenzte Plasmasysteme (BoundedPlasma Systems, BPS’s) wurde bereits 1944 vonPierce vorgeschlagen und erweiterte das Plas-maoszillationsmodell von Tonks und Lang-muir (1929) um physikalisch motivierte Rand-bedingungen. Modellierung und Simulationvon BPS stellen dabei einen wichtigen Bei-trag für die quantitative Beschreibung vonsolchen Plasmasystemen dar, deren Ergebnis-se für viele Plasmabereich relevant sind. EineÜbersicht über verschiedene mikroskopischeBPS-Simulationen ist in Kuhn (1994) zu finden.

2 Pierce-Diode

Das Modell der Pierce-Diode ist ein verhältnis-mäßig einfaches – jedoch wichtiges und reprä-sentatives – ’archetypisches’ BPS-Modell. Mitihm konnten erstmals Instabilitätsmechanis-men beschrieben werden, die als Ursache für

– 1 –

802 – Simulation einer Pierce-Diode Tobias Krähling

L

Φ(0) = 0 Φ(L) = 0

El(0) = −0vEl(0) = v0

Elektronenstrom

neutralisierenderIonenhintergrund

Ion = 0

Abbildung 1: Schematischer Aufbau einerPierce-Diode

Plasmaoszillationanregung interpretiert wer-den können. Bei dem ursprünglichen, von Pier-ce (1944) vorgeschlagenen, Modell wird davonausgegangen, das sich ein kalter Elektronen-strahl vor einem neutralisierenden Ionenhin-tergrund von der Kathode zur Anode bewegt(siehe Abbildung 1). Anode und Kathode sinddabei kurzgeschlossen und auf Nullpotential,externe magnet- und elektrische Felder wirkennicht auf das BPS ein.

Bereits dieses, wohl einfachste, geschlosseneModell einer thermionischen Entladung zeigteine vielfältige Dynamik, verursacht durchverschiedene Instabilitäten, die auf komplexe-re Plasmasysteme übertragen werden können.Stationäre Gleichgewichte (Stabilität der ho-mogenen Lösung) können bereits durch eineinfaches hydrodynamisches Modell unter-sucht werden, für Kathodenoszillationen istein kinetisches Modell notwendig1. Weiterhinkönnen zeitlich chaotische Plasmaoszillatio-nen auftreten.

1 Kathodenoszillationen sind Relaxationsoszillationen,die durch die Zunahme der negativen Raumladungvor der Kathode, die sogenannte virtuelle Kathode,entstehen.

2.1 Plasmaoszillationsmodell vonTonks und Langmuir

Im Modell von Tonks und Langmuir (1929)werden die Elektronen als eindimensionaleFlüssigkeit betrachtet. Die Beschreibung er-folgt über die Kontinuitätsgleichung

∂%

∂t+∂∂x

(v%) = 0, (2.1)

der Impulsbilanz

∂v∂t

+12∂v2

∂x=

eme

∂∂x

Φ (2.2)

und der Poisson-Gleichung

∂2Φ

∂x2 = − 1ε0

(%0 + %), (2.3)

aus der sich selbstkonsistent das elektrischePotential Φ ergibt. Dabei bezeichnet, wie auchim Folgenden, %die Elektronendichte, %0 die Io-nendichte, v die Elektronengeschwindigkeit, edie Elementarladung, me die Elektronenmasseund ε0 die Dielektrizitätskonstante.

2.2 Modifikationen im Pierce-Modell

Im Pierce-Modell (Pierce, 1944) ist die entschei-dende Veränderung zum Plasmaoszillations-modell die Einführung von Randbedingungender Form:

Φ(0) = Φ(L) = 0, (2.4a)

%(0) = −%0 und (2.4b)

v(0) = v0 (2.4c)

mit der Systemlänge L und der Kathode beix = 0 (siehe auch Abbildung 1). Diese Rand-bedingungen haben entscheidenden Einflussauf die Dynamik der Pierce-Diode und sindnachfolgend aufgeführt:

– 2 – v1.0 – 14.11.2009


I An der Kathode werden dem System kon-tinuierlich bewegte Ladungsträger (Elek-tronen) hinzugefügt, was einer kontinuier-lichen Energiezuführung entspricht undeinem treibenden Term darstellt (Bedin-gungen für Ladungsdichte und Elektro-nengeschwindigkeit).

I An den beiden Plasmabegrenzungen (Ka-thode und Anode) tritt Teilchenverlustauf, was durch die Bedingungen für dasPotential beschrieben ist. Dieser Teilchen-verlust kann als Dissipationsmechanis-mus des Systems angesehen werden.

Das Modell beschreibt demnach ein Systemmit kontinuierlicher Ladungsträgerquelle und-senke – ein getriebenes, dissipatives System.

2.2.1 Dimensionslose Formulierung

Die Gleichungen (2.1), (2.2) und (2.3) sowiedie Randbedingungen (2.4) lassen sich mittelsder Variablentransformationen

x ↦→ xL, t ↦→ v0t

L, % ↦→ %

%0,

v ↦→ vv0

und Φ ↦→ eΦmv2

0

in die dimensionslose Form∂%

∂t+∂∂x

(v%) = 0 (2.5)

∂v∂t

+12∂v2

∂x=∂Φ∂x

(2.6)

∂2Φ

∂x2 = −α2%0(1 + %) (2.7)

Φ(0) = Φ(L) = 0, %(0) = −1, v(0) = 1 (2.8)

überführen. Für den einzig verbleibende Para-meter α gilt dabei

α BωPLv0

(2.9)

mit der Plasmafrequenz

ωP B

√%0eε0me

(2.10)

Abbildung 2: Anwachsraten der instabilstenModen in Abhängigkeit von α (RUB, 2003)

Die homogene Lösung mit % = −1 und v = 1ergibt sich direkt aus den Gleichungen (2.5)–(2.8). In Abbildung 2 ist für die homogeneLösung die auf die Plasmafrequenz ωP be-zogene Anwachsrate −ℑ(ω/ωP) der jeweilsinstabilsten Mode in Abhängigkeit des frei-en Parameters α dargestellt. Für den Bereichπ < α < 2π tritt die klassische, nicht oszil-latorische, Pierce-Instabilität auf und im Be-reich 2π < α < 3π die oszillatorische Pierce-Instabilität. Stabile Lösungen existieren dage-gen im Bereich 0 < α ≤ π, α = 2π und α ≈ 3π.

2.3 Zu untersuchende Pierce-Dioden-Modifikation

Zu dem im vorherigen Abschnitt aufgeführ-ten klassischen Modell der Pierce-Diode wur-den bereits verschiedene Modifikationen un-tersucht, z. B. die Beschaltung mit einem exter-nen RLC-Kreis (Lawson, 1989a,b; Kuhn, 1994).Im Rahmen des F-Praktikums soll hier diefolgende Modifikation der klassischen Pierce-Diode untersucht werden:

I Auf beiden Seiten der Pierce-Diode wer-den kontinuierlich Elektronen mit Ge-schwindigkeit vL = v bzw. vR = −v er-zeugt – es werden also zwei entgegen-gesetzte Elektronenstrahlen erzeugt und

– 3 – v1.0 – 14.11.2009


beide Begrenzungen sind gleichzeitig Ka-thode und Anode.

I In einem zweiten Schritt soll der Ge-schwindigkeitsbetrag der in die Pierce-Diode eintretenden Elektronen nicht kon-stant, sondern Normalverteilt (Gauß-Verteilung) sein.

Hierdurch werden für die erste Modifika-tion, neben den Randbedingungen in Glei-chung (2.8) die Randbedingungen %(L) = −1und v(L) = −1 eingeführt. Ziel ist es, diesesveränderte Pierce-Modell hinsichtlich seinerStabilität (Bifurkation) numerisch zu untersu-chen.

3 Particle-in-Cell Simulation(PIC)

Zur Simulation der modifizierten Pierce-Diodemüssen die Newton’schen Bewegungsglei-chungen

xi = vi (3.1a)

vi = Fi(x1, . . . ,xN) (3.1b)

für die Elektronenbewegung verwendet wer-den, da das hydrodynamische Modell nichtmehr ausreicht. Eine vollständige mikroskopi-sche Simulation des Systems wäre aufgrundder großen Anzahl an Elektronen nicht durch-zuführen, so dass mehrere Elektronen zusogenannten Makropartikeln zusammenge-fasst werden. Die Möglichkeit, das Systemdurch die Definition der Makropartikel zu si-mulieren, wird jedoch durch einen höherenRauschanteil ’erkauft’.

Bei der PIC-Simulation, die seit nun gut 50Jahren für die Plasma-Simulationen verwen-det wird, werden Partikeltranslationen und-geschwindigkeitsänderungen in einem konti-nuierlichen Phasenraum durchgeführt, wäh-rend Teilchendichten und selbstkonsistente

elektrische Felder auf einem festen Gitter be-rechnet werden (siehe auch Dawson (1983)).Hierfür wird die Dichte jedes einzelnen Par-tikels in einer Gitterzelle auf die umgeben-den Gitterknoten gewichtet und die Poisson-Gleichung gelöst. Anschließend wird die Kraft,die auf das Partikel im Phasenraum wirkt, be-stimmt und die Geschwindigkeit des Teilchensentsprechend geändert.

3.1 Implementierung

Die PIC-Simulation der modifizierten Pierce-Diode wurde in C++ (GNU C/C++ compiler)mit einem graphische QT-Benutzerinterface(QT-Bibliothek, http://qt.nokia.com) imple-mentiert, für die Compilersteuerung wur-de CMake (http://www.cmake.org) verwen-det. Zusätzlich wurde eine eigene Biblio-thek für Datenstrukturen, Portabilitätsdefini-tionen usw. eingesetzt (libsbstd in Version 1.5).Aus der Boost-Bibliothek (http://www.boost.org) wurde der Zufallszahlengenerator (Lag-ged Fibonacci Generator mit p = 23209 undq = 13470) und die Normalverteilung verwen-det. Ausgeführt wurde die Simulation auf ei-nem Intel Core2Duo E8400 unter einem 64bit-Linux.

Aus den Randbedingungen (siehe Abschnitt2.3) folgt, dass periodisch im Zeitabstand Tneue Partikel am linken und rechten Rand mitden Geschwindigkeiten v bzw.−v hinzugefügtwerden. Diese Teilchen tragen die Ladung qc =

−T/2, was aus der Bedingung an % folgt (beimklassischen Pierce-Modell tragen die Teilchendie Ladung qc = −T, da nur an der linken Seiteneue Teilchen hinzugefügt werden).

Die Klasse Particle modelliert ein einzelnes Par-tikel. Über die Methode setF kann die Kraft,die auf das Partikel einwirkt, gesetzt werden.moveV ändert die Geschwindigkeit anhandder auf das Partikel einwirkenden Kraft und

– 4 – v1.0 – 14.11.2009

http://qt.nokia.com

http://www.cmake.org

http://www.boost.org

http://www.boost.org


der Einwirkzeit, moveX analog die Position an-hand der aktuellen Geschwindigkeit und derübergebenen Zeit. FGenerateVelocity ist eineFunktorklasse, die dazu verwendet wird, dieGeschwindigkeit eines neuen Partikels zu ge-nerieren – je nach Simulationsmodi liefert derFunktor entweder eine konstante Geschwin-digkeit oder eine gaußverteilte Geschwindig-keit bei vorgegebener mittleren Geschwindig-keit und Standardabweichung (mittels derKlasse RNGGauss).

In der Klasse PIC ist die Simulation implemen-tiert und hält alle in der Simulation befind-lichen Partikel in einer doppelt verkettetenListe, die über die Methode getList nur-lesbarzugreifbar ist (wird von der GUI-Klasse fürdie Darstellung der Partikel benötigt). Alleeinzelnen Schritte, die für die PIC-Simulationbenötigt werden, sind gekapselt und könnennur in ihrer Gesamtheit über die öffentlicheMethode integrate ausgeführt werden. Diessind, ihrer Reihenfolge nach (mit ∆x wird dieLänge einer Gitterzelle bezeichnet, die sich beiN Gitterknoten über ∆x = (N − 1)−1 ergibt):

1. Verschiebe alle Partikel gemäß der Be-wegungsgleichung (3.1) um einen halbenZeitschritt dt, was mit der Methode pMo-veX(dt) für alle Partikel in der Liste erfolgt.Dies entspricht dem Leapfrog-Verfahren(siehe Abbildung 3):

xn+1/2i = xn−1/2

i + ∆tvni (3.2)

2. Entferne alle Partikel, die sich außerhalbder Plasmabegrenzung befinden, also fürdie gilt: x < 0 ∨ x > 1 (Methode pElimina-teParticles).

3. Ermittlung der Dichte % auf einem re-gelmäßigen Gitter mit Knotenanzahl N(Methode pWeightRho). Diese Wichtungauf die Gitterknoten erfolgt über die Glei-

t

(n − 1)∆t

vn−1i

(n)∆t

vni

(n + 1)∆t

vn+1i

(n + 2)∆t

vn+2i

(n − 12)∆t

xn− 1

2

i

(n + 12)∆t

xn+ 1

2

i

n+12 ,En+ 1

2 ,Fn+12

(n + 32)∆t

xn+ 3

2

i

(n + 52)∆t

xn+ 5

2

i

Abbildung 3: Leapfrog-Verfahren

chungen:

%i = q j(i + 1)∆x − x j

∆x2 (3.3a)

%i+1 = q jx j − i∆x

∆x2 (3.3b)

wobei x j und q j die x-Koordinate und La-dung des aktuellen Partikels bezeichnetund i die Gitterzelle, indem sich das Parti-kel befindet. Die Gesamtdichte an einemKnoten setzt sich aus allen Einzelpartikel-dichten an diesem Knoten additiv zusam-men.

4. Lösung der Poisson-Gleichung (2.7) mitden Randbedingungen auf dem regelmä-ßigen Gitter mittels Tridiagonalbandlöser(Methode pPoisson), näheres in Abschnitt3.2. Das elektrische Feld auf den Gitter-knoten ergibt sich dann zu:

Ei = −Φi+1 −Φi−1

2∆x(3.4)

5. Wichtung der Kraft auf ein Partikel voneinem Gitterknoten auf den Ort des Parti-kels über

F j = − (i + 1)∆x − x j

∆xEi −

x j − i∆x∆x

Ei+1

(3.5)wobei j das konkrete Partikel und i dieGitterzelle bezeichnet, indem sich das Par-tikel befindet (Methode pWeightForce).

– 5 – v1.0 – 14.11.2009


6. Gemäß der Bewegungsgleichung (3.1)werden die Geschwindigkeiten aller Parti-kel um einen Zeitschritt dt erhöht (Metho-de pMoveV), wieder mittels des Leapfrog-Verfahrens (siehe Abbildung 3):

vn+1i = vn

i + ∆tFi({xn+1/2j }) (3.6)

7. Verschiebe alle Partikel um einen halbenZeitschritt wie in 1. – somit wurde dasTeilchen um einen ganzen Zeitschritt dtverschoben.

8. Bei Bedarf werden an beiden Seiten derPierce-Diode (bzw. an der linken Seiteim klassischen Pierce-Modell) jeweils einneues Partikel in die Diode injiziert, fallsdie letzte Injektion länger als das angege-bene Injektionsintervall T zurückliegt.

Iterativer Aufruf der Methode integrate ermög-licht somit, das Verhalten der Partikel in derPierce-Diode in seiner zeitlichen Entwicklungzu beobachten, wobei die Zeitauflösung überdie Zeitschrittwahl dt erfolgt.

3.2 Tridiagonalbandlöser

Zur Lösung der Poisson-Gleichung (2.7) wirdein Tridiagonalbandlöser verwendet, da sichdie Poisson-Gleichung als Gleichungssystemmit tridiagonale Matrix beschreiben lässt. Derverwendete Algorithmus ist dabei aus Presset al. (1992, S. 50ff) und Engeln-Müllges undReutter (1993, S. 87ff) entnommen und entspre-chend adaptiert worden.

Eine Gleichungssystem Ax = a mit tridiagona-ler Matrix hat die Form:

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

a11 a12 0

a21 a22 a23 0

0 a32 a33 a34 0. . .

. . .. . .

. . .

0 an−1n−2 an−1n−1 an−1n

0 ann−1 ann

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠⏟ ⏞ A

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

x0

x1

x2

...

xn−1

xn

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠⏟ ⏞ x

=

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

a0

a1

a2

...

an−1

an

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠⏟ ⏞ a

Ist L eine bidiagonale untere Dreiecksmatrixund R eine normierte bidiagonale obere Drei-ecksmatrix, so lässt sich Ax = a mit der Zer-legung A = LR in ein äquivalentes SystemRx = r überführen. Hierzu sind die Schritte

1. Faktorisierung: A = LR⇒ L,R

2. Vorwärtselimination: a = Lr⇒ r und

3. Rückwärtselimination: Rx = r⇒ x

durchzuführen. Dabei werden die Elementeder Matrizen A, L, R und des Vektors r wiefolgt bezeichnet und als Vektor abgespeichert:

A =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

d0 c0 0

e1 d1 c1 0

0 e2 d2 c2 0. . .

. . .. . .

. . .

0 en−2 dn−2 cn−2

0 en−1 dn−1

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠, r =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝r1

r2

...

rn−1

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

R =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

1 γ0

1 γ1

. . .. . .

1 γn−2

1

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠, L =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝α0

β1 α1

. . .. . .

βn−1 αn−1

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠Koeffizientenvergleich von A mit LR liefert:

∀i : βi = ei.

Nachfolgend wird der Algorithmus ausEngeln-Müllges und Reutter (1993) für diebenötigten drei Schritte aufgeführt, wobei βi

direkt mit ei ersetzt wurde:

– 6 – v1.0 – 14.11.2009


Algorithmus

1. Zerlegung: A = LR1.1 α0 = d0

1.2 γ0 = c0/α0

1.3 Berechne für alle i = 1 . . . n − 1:1.3.1 αi = di − eiγi−1

1.3.2 γi = ci/αi

1.4 αn−1 = dn−1 − en−1γn−2

2. Vorwärtselimination: a = Lr2.1 r0 = a0/d0

2.2 Berechne für alle i = 1 . . . n − 1:2.2.1 ri = (ai − eiri−1)/αi

3. Rückwärtselimination: Rx = r3.1 xn−1 = rn−1

3.2 Berechne für alle i = (n − 2) . . . 0:3.2.1 xi = ri − γixi+1

3.2.1 Übertragung auf die Poisson-Gleichung

Die dimensionslose Poisson-Gleichung (2.7)wird bei der Betrachtung auf einem diskretenGitter mit den Gitterpunkten j (mit %0 = 1 undfür j = 1 . . . n − 2) zu

Φ j−1 − 2Φ j + Φ j+1

(∆x)2 = −α2(1 + % j) (3.7)

und aus den Randbedingungen folgt Φ0 =

Φn−1 = 0. Dies kann als Matrixgleichung derForm

AΦ = −(∆x)2α2(1 + %)⏟ ⏞ =a

(3.8)

aufgefasst werden. Die Matrix A hat dabei dieForm:

A =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

1 0 0

1 −2 1

0 1 −2 1. . .

. . .. . .

. . .

0 1 −2 1

0 0 1

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

Die Schritte Zerlegung und Vorwärtseliminati-on wurden zusammengefasst und sparen hier-durch die Verwendung von mehreren Vekto-ren. Man erhält

γi =

⎧⎪⎪⎨⎪⎪⎩0 ; i = 0ciαi

= 1di−eiγi−1

→ 12−γi−1

; i = 1 . . . n − 2

und

ri =

⎧⎪⎪⎨⎪⎪⎩0 ; i = 0ai−eiri−1αi→ γi(ri−1 − ai) ; i = 1 . . . n − 1

wobei das negative Vorzeichen von di heraus-gezogen und in beiden Gleichungen angepasstwurde, was durch den Pfeil gekennzeichnetist.

Aus der Randbedingung Φ(L) = 0 folgt für dieRückwärtselimination

xi =

⎧⎪⎪⎨⎪⎪⎩0 ; i = n − 1

ri − γixi+1 → ri + γixi+1 ; i = (n − 2) . . . 0

wobei auch hier die Anpassung bezüglich desVorzeichens von di durchgeführt wurde. xi ent-hält nun die Potentiale auf den Gitterpunkten,die den Randbedingungen Φ(0) = Φ(L) = 0genügen. In der Implementierung sind dieMatrixelemente in Vektoren gespeichert undfolgendermaßen benannt:

I γi⇔ TVector<double>::a

I ri⇔ TVector<double>::b

I ai⇔ TVector<double>::f

I xi⇔ TVector<double>::pV_phi

4 Programm

Die graphische Benutzeroberfläche des Simu-lationsprogramms ist in Abbildung 4 wieder-gegeben und bietet die Möglichkeit, den Pha-senraum der Partikel in der Pierce-Diode dar-zustellen und die möglichen Simulationspara-

– 7 – v1.0 – 14.11.2009


Abbildung 4: GUI der Pierce-Dioden Simulation

– 8 – v1.0 – 14.11.2009


meter einzustellen. Auf der linken Seite befin-det sich die graphische Darstellung des Pha-senraums, wobei horizontal die x-Koordinateund vertikal die Geschwindigkeit v aufgetra-gen ist. Die x-Koordinate verläuft in den di-mensionslosen Variablen von x = 0 bis x = 1,bei der Geschwindigkeitsachse ist in der Mit-te v = 0 und die in der Simulationssteue-rung eingestellte Geschwindigkeit erstrecktsich über einviertel der Höhe. Auf der rechtenSeite ist die Steuerung der Simulation mög-lich, deren Einstellungsmöglichkeiten zumTeil selbsterklärend sind, zum anderen zu derBeschreibung im vorherigen Abschnitt korres-pondiert. Im Reiter Simulationsparameter kön-nen zunächst die grundlegenden Parametereingestellt werden, wobei hier über Simulati-onstyp die klassische Pierce-Diode oder diemodifizierte Pierce-Diode mit zwei entgegen-gesetzten Kathoden ausgewählt werden kann,dt gibt die Zeitauflösung an und T das Injek-tionsintervall. Ist als Geschwindigkeitsdistri-bution Gauß ausgewählt, so können über denReiter Gauß-Distribution die Standardabwei-chung σwie auch der Seed für den Zufallszah-lengenerator eingestellt werden (bei auto wirdein automatischer Seed bei jedem Starten derSimulation erzeugt und angezeigt).

Alle Änderungen an den Simulationsparame-tern werden erst beim nächsten Starten einerSimulation aktiviert, Änderungen an Farbenund der Aktualisierungsrate können währendder Simulation durchgeführt werden und er-folgen bei der nächsten Aktualisierung.

5 Simulationsergebnisse

Für alle Simulationen wurden 100 Gitterkno-ten, eine Iterationszeit von dt = 0,001 undeine Injektionszeit von T = 0,005 verwendet.In der ersten Variation wurde dann für dieZweistrom-Piercediode (Modus: diametrical)

bei einer konstanten Injektionsgeschwindig-keit v = 1,0 der freie Parameter α im Bereich0 − 10 variiert und das Verhalten der Pierce-Diode analysiert.

5.1 Variation des freien Parameters αbei konstanter Geschwindigkeit

Um das Verhalten der Pierce-Diode bei Variati-on des freien Parameters α zu bestimmen, wur-de jeweils die Simulation bis zu einer Zeit vont ≃ 1500 − 2000 durchgeführt – die Ergebnissesind in den Tabellen 1 und 2 wiedergegeben.Im Bereich 0 ≤ α . π war das Verhalten stabilund beide entgegengesetzt injizierte Elektro-nenstrahlen beeinflussen sich nur sehr gering-fügig. Dabei wurde die Lösung als stabil ange-sehen, wenn sich beide Elektronenstrahlen vonder Ursprungskonfiguration nicht mehr als ei-ne Strahlbreite veränderten2. Dieses Verhaltenentspricht demjenigen der klassischen Pierce-Diode, eine schematische Darstellung der Pha-senraumkonfiguration ist in Abbildung 5(a)zu finden. Ab α ∼ π ist zunächst eine instabilePhasenraumkonfiguration zu beobachten, dieanschließend in eine stationäre Lösung über-geht (siehe Abbildung 5(b)). Dabei ist die in-stabile Phasenraumkonfiguration nur dadurchgekennzeichnet, dass die Amplitude zunächststetig steigt, bis der stationäre Zustand erreichtist – mit zunehmenden α steigt die Ampli-tude des stationären Zustands immer weiteran, wohingegen die charakteristische Zeit, biszu der der stationäre Zustand erreicht wird,kontinuierlich sinkt. Analog zur klassischenPierce-Diode ist der stationäre Zustand eben-falls durch ein Maximum gekennzeichnet.

Ab α ≥ 6,0785 ist die instabile Phase nicht nurdurch ein Anwachsen der Amplitude bis zum

2 Dies konnte gut an dem oberen Strahl bestimmt wer-den, da sich dieser in der Ursprungskonfigurationdirekt unterhalb der Hilfslinie befindet und als stabilangesehen wurde, wenn sich dieser nicht in einemBereich vollständig oberhalb der Hilfslinie befindet.

– 9 – v1.0 – 14.11.2009


Tabelle 1: Ergebnisse der Simulation bei Va-riation des freien Parameters α und konstanterGeschwindigkeit

α Verhalten

0,0 stabil...

...

3,1764 stabil

3,1765 anwachsend, dann stationär miteinem Maximum

... (Amplitude steigt, Zeit bis zumstationären Zustand sinkt)

6,0784 stationär mit einem Maximum

6,0785 zunächst instabil (Einschlussvon Partikel), dann stationär

......

7,8500 zunächst instabil (Einschlussvon Partikel), dann stationär

7,8501 Bildung einer virtuellen Katho-de rechts, die nach links wan-dert; dann quasistationär.Die Bildung der virtuellen Ka-thode erfolgt bei steigendem αmal von links, mal von rechts,teilweise auch von der Mitte aus-gehend zu beiden Seiten, beson-dere Punkte sind

. 8,7 nachfolgend aufgeführt.

> 8,7 zunächst chaotisch (gebildetevirtuelle Kathoden werden wie-der aufgebrochen, anschließendquasistationär

Tabelle 2: Ergebnisse der Simulation bei Va-riation des freien Parameters α und konstanterGeschwindigkeit (forts.)

besondere Punkte

α Verhalten

7,8612 virtuelle Kathode l→ r

7,8613 instabil→ stationär...

...

7,8616 instabil→ stationär

7,8617 virtuelle Kathode l← r


7,8760 instabil→ stationär...

...







0

v

−v

0 1x→

(a) stabil

0

v

−v

0 1x→

(b) stationär

b

b

b

b

b

b

b

b

b

b

bb

b

b

b

b

b

b

b

b

b

b

b

b

b

bb

b

0

v

−v

0 1x→

(c) instabil

0

v

−v

0 1x→

(d) quasistationär

Abbildung 5: Schematische Darstellung derbeobachteten Phasenraumkonfigurationen ander modifizierten Pierce-Diode

– 10 – v1.0 – 14.11.2009


stationären Zustand gekennzeichnet – viel-mehr werden Partikel innerhalb der Pierce-Diode eingeschlossen, was jeweils durch ei-ne Schlängelbewegung des jeweils hinterenTeils des Partikelstrahls (in Abbildung 5(c)durch die farbigen Bereiche gekennzeichnet)erfolgt. Wenn alle so eingeschlossenen Parti-kel die Pierce-Diode verlassen konnten gehtdiese wieder in den stationären Zustand miteinem Maximum (Abbildung 5(b)) über, wobeidie Amplitude proportional zum Parameterα ist. Bei diesem Verhalten ist kein klares Ver-halten der charakteristischen Zeit, bis zu derder stationäre Zustand erreicht wird, und demParameter α zu beobachten – diese Zeit zeigteher ein chaotisches Verhalten.

Das im vorherigen Abschnitt beschriebeneVerhalten ist bis α = 7,8500 zu beobachten,anschließend verändert sich das Verhaltender modifizierten Pierce-Diode grundlegend.Nach einer Zeit mit chaotischem Verhaltenbildet sich überwiegend an einer der beidenSeiten eine virtuelle Kathode, die sich anschlie-ßend bis zur anderen Seite ausbreitet, so dasssich vor beiden Kathoden jeweils eine virtuelleKathode befindet – teilweise erfolgt die Bil-dung der virtuellen Kathode in der Mitte undbreitet sich zu beiden Seiten aus. Innerhalb derPierce-Diode bildet sich ein Steifen aus, indemsich viele Partikel mit geringer Geschwindig-keit befinden (etwa in einem Intervall ± 1

4 v).Dieser Zustand ist an sich nicht-stationär, än-dert sich jedoch über größere Zeitskalen beiMittlung nicht, weswegen dieser als quasista-tionärer Zustand bezeichnet wird. In Abbil-dung 5(d) ist dieser schematisch als Mittlungüber eine größere Zeitskala dargestellt. Auffäl-lig ist noch, dass das sich ausbildende Bandmit geringer Geschwindigkeit zeitlich oszil-liert und ein oder zwei Knoten bildet. Es konn-te beobachtet werden, dass diese Oszillationdes Bandes zeitlich nicht stabil ist und das imZeitverlauf Oszillationen mit einem als auchmit zwei Knoten beobachtet werden können

– beide virtuellen Kathoden wirken dabei alstreibendes System für diese Schwingungen.

Bemerkenswert sind die Bereich α ≈ 7,8614,α ≈ 7,8763 und α ≈ 7,9758 – in diesen dreiBereichen konnte, nach einem instabilen Zu-stand mit Partikeleinschluss, ein stationärerZustand beobachtet werden, analog zum Be-reich 6,0785 ≤ α ≤ 7,8500 (detaillierte Auf-führung in Tabelle 2). Dabei fällt auf, dassdie Bildung dieser stationären Inseln immerzwischen einem Wechsel derjenigen Seite auf-tritt, wo die erste virtuelle Kathode gebildetwird. Ein Wechsel der ersten virtuellen Ka-thodenseite führt jedoch nicht zwangsläufigzu einem dazwischenliegenden stationärenBereich (z. B. bei α = 7,9791 − 7,9792 undα = 8,2433 − 8,2434).

Ab α & 8,7 zeigt sich vor der Bildung desquasistationären Zustands ein chaotisches Ver-halten. In dieser chaotischen Phase könnenauch bereits gebildete virtuelle Kathoden zer-stört werden. Desweiteren können sich in demBand niedriger Partikelgeschwindigkeit Berei-che bilden (häufig Kreise, teilweise auch El-lipsen), in denen sich keine Partikel befinden– das Band wird quasi an einer Stelle geteiltund fließt um diesen Störbereich herum. Die-ser Störbereich pendelt dann häufig zwischenden beiden virtuellen Kathoden und wird imzeitlichen Verlauf an den virtuellen Kathodenverkleinert und schlussendlich vernichtet (die-se Störbereiche können sich immer bei der Aus-bildung des Bandes beim quasistationären Zu-stand zunächst einstellen). Nach unterschied-lichen Zeiten geht das System dann in denquasistationären Zustand mit der Bandoszilla-tion über.

Zusammenfassend ähnelt das Verhalten dermodifizierten Zwei-Strahl Pierce-Diode demje-nigen der klassischen Pierce-Diode (in demhier betrachteten Bereich 0 ≤ α ≤ 10). Diestabilen Bereiche bei α = 2π und α ≃ 3π der

– 11 – v1.0 – 14.11.2009


klassischen Pierce-Diode sind jedoch nicht vor-handen und als stabile Lösungen kommen nurdiejenigen mit einem Maximum in Betracht.

5.2 Einfluss der Geschwindigkeits-verteilung

In einem weiteren Teil wurde die konstanteGeschwindigkeit, mit der die Partikel in diePierce-Diode injiziert werden, durch eine Ge-schwindigkeitsverteilung mit Normalvertei-lung ersetzt. Motivation für diese Betrachtungist die Möglichkeit, kleine Störungen des Sy-stems – hier durch Variation der Geschwindig-keiten (Variation der Standardabweichung σ) –zu betrachten und den Einfluss auf die im vor-herigen Abschnitt gefundenen stabilen undstationären Bereiche abzuschätzen. Im folgen-den werden die Grenzen für die Verhaltensän-derung ab einer bestimmten Standardabwei-chung jeweils dann akzeptiert, wenn alle fünfdurchgeführten Simulationsläufe bis t ≃ 2000das das selbe Verhalten zeigen wie bei derSimulation mit konstanter Geschwindigkeit.

Zunächst wurde der Bereich α = 3,1764 unter-sucht, der ein stabiles Verhalten zeigte (sieheTabelle 1). Dieser Bereich scheint für kleineStörungen sehr empfindlich zu sein, stabileLösungen konnten nur für σ ≤ 0,005 beobach-tet werden. Für größere σ geht das Systemnach unterschiedlichen Zeiten in den quasista-tionären Zustand über, wobei die virtuellenKathoden wesentlich näher zusammenliegen– liegt das σ nur geringfügig über der Gren-ze von 0,005, so gibt es teilweise auch sta-bile Lösungen3. Zu beobachten ist, dass beigrößeren Störungen, d.h. bei größerem σ, dasSystem schneller in den quasistationären Zu-stand übergeht. Es zeigt sich, das bei kleineremα die Störanfälligkeit stark absinkt. Für α = 3,0

3 Aus diesem Grund wurde entschieden, fünf Simulati-onsläufe für die Grenzenbestimmung zu verwenden.

liegt die Grenze bei σ ≃ 0,030 und für α = 2,5bei σ ≃ 0,125.

Als zweiter Verhaltensübergangsbereich wur-de α ≤ 7,85 untersucht und folgende Ergebnis-se erhalten:

I Für α & 7,65 ist bereits eine Störung mitσ = 0,001 so groß, dass das System teil-weise in den quasistationären Zustandübergeht.

I Für α = 7,60 darf σ = 0,002 nicht über-steigen, um den stationären Zustand zuerreichen.

Es zeigt sich, das im Bereich α . 7,85 das Errei-chen des stationären Zustandes sehr störanfäl-lig ist und bereits kleine Abweichungen der Ge-schwindigkeit zum quasistationären Zustandführen können. Auffällig ist noch, dass derZustandsweg, den dass System geht (stationäroder quasistationär), sehr früh festgelegt wirdund sich dann nicht mehr ändert. Erreicht dasSystem den instabilen Zustand, indem ein Teilder Partikel eingeschlossen werden (siehe Ab-bildung 5(c)), so führt dieser immer zum statio-nären Zustand – andererseits nimmt die Wahr-scheinlichkeit, diesen stationären Zustand zuerreichen, mit steigender Störung ab.

Als letzter Bereich, der hier untersucht werdensoll, ist die stationäre Insel bei α ≈ 7,8763.Hier zeigt sich, dass bereits eine Störung mitσ = 0,001 diese stationäre Insel zerstört undnur noch quasistationäre Lösungen auftreten.Dabei ist zu beobachten, dass die Bildung dervirtuellen Kathode nicht mehr immer von nureiner Seite, sondern mal von links, mal vonrechts erfolgt.

– 12 – v1.0 – 14.11.2009


6 Zusammenfassung

Es zeigt sich, dass die durchgeführten Modi-fikationen an dem klassischen Pierce-Modellzu einigen Verhaltensänderungen führen. Fürα . π ist die Lösung stabil, während fürα . 7,85 nach einer instabilen Phase das Sy-stem in einen stationären Zustand übergeht.Bei weiterer Erhöhung des freien Parametersα ergibt sich ein quasistationärer Zustand, wo-bei drei stationäre Inseln gefunden werdenkonnte.

Weiterhin zeigte sich, indem die konstante Ge-schwindigkeit durch eine normalverteilte Ge-schwindigkeit ersetzt wurde, dass die Grenzensehr störanfällig sind – dies trifft insbesondereauf den stationären Bereich zu, die stationärenInseln verschwinden bei kleinsten Störungenbereits völlig.

Literatur

John M. Dawson. Particle simulation of plasmas. Reviewsof Modern Physics, 55(2):403–447, 1983. DOI: 10.1103/

RevModPhys.55.403.

G. Engeln-Müllges und F. Reutter. Numerik-Algorithmenmit ANSI C-Programmen. BI Wissenschaftsverlag, 1993.ISBN 3-411-16291-0.

S. Kuhn. The Physics of Bounded Plasma Systems(BPS’s): Simulation and Interpretation. Contrib. PlasmaPhys., 34(4):495 – 538, 1994.

William S. Lawson. The Pierce diode with an externalcircuit. I. Oscillations about nonuniform equilibria.Physics of Fluids B: Plasma Physics, 1(7):1483 – 1492,1989a. DOI: 10.1063/1.858925.

William S. Lawson. The Pierce diode with an externalcircuit. II. Chaotic behavior. Physics of Fluids B: PlasmaPhysics, 1(7):1493 – 1501, 1989b. DOI: 10.1063/1.859199.

J. R. Pierce. Limiting Stable Current in Electron Beamsin the Presence of Ions. Journal of Applied Physics, 15:721 – 726, 1944.

William H. Press, Saul A. Teukolsky, William T. Vetterling,und Brian P. Flannery. Numerical Recipes in C. The Artof Scientific Computing. Cambridge University Press,2. Auflage, 1992. ISBN 0-521-43108-5.

F-Praktikum Physik: Simulation einer Pierce-Diode. RUB,2003. Versuchsanleitung zum Praktikum für Fortge-schrittene an der Ruhr-Universität Bochum, Versionvom 29.09.2003.

Lewi Tonks und Irving Langmuir. Oscillations in IonizedGases. Physical Review, 33:195 – 210, 1929.

– 13 – v1.0 – 14.11.2009

802 - Simulation einer Pierce-Diode -...

Documents

Transcript of 802 - Simulation einer Pierce-Diode -...