Numerische Modellierung/Ablaufschema

Modellierungsansatz: Der Modellierung liegen die in Abschnitt 4.6 getroffenen Annahmen zugrunde. Die Implementierung erfolgt in 3D und als Referenz dienen eine Reihe von Arbeiten, in deren Rahmen vergleichbare Codes entwickelt wurden [10; 121; 125–127]. Bei dem gewählten Ansatz werden die Elektronen als Kontinu-um simuliert, während die Ionen als Teilchen mit Hilfe der Particle-in-Cell (PIC)-Methode [110; 128–131] nachgebildet werden. Der Grundgedanke letzterer (PIC)-Methode ist es, nicht die Kraft von jedem Teilchen auf jedes Teilchen simulieren zu müssen.

Stattdessen wird, wie in Abb. 7.2 ersichtlich und im Anhang in Abschnitt B.2 detail-liert, das Simulationsgebiet in Voxel unterteilt. Das Potential sowie die Elektronen-und Ionendichten werden auf den Eckpunkten (Stützpunkte des Feldes) berech-net. Durch dieses Potential lässt sich die „Weitdistanzwechselwirkung“ der Teilchen, unter Vernachlässigung der „Kurzdistanzwechselwirkungen“ innerhalb einer Zelle, berechnen. Dieses Vorgehen erlaubt es die Komplexität (im Sinne der Informatik) als Funktion der Teilchenzahl zu reduzieren. Teilcheninteraktionen (Reflexion oder Löschung¹) können mit der Wand jedes Voxels stattfinden.

1 Löschung heißt hierbei, dass das Teilchen mit seinen Orts- und Geschwindigkeitsdaten aus dem Datenspeichern entfernt wird. Dies kann z. B. notwendig sein, wenn ein Teilchen das Simulationsgebiet verlässt.

96

7.4 Simulation der Ionenoptik

DIN A4 210 mm x 297 mm

Seite ist 15 cm breit. Jeder Block ist 4,6 mm Breit.

So ergibt sich links und rechts 0,2 mm Abstand vom Rand.

ja nein

Abbildung 7.9: Ablaufschema für die Simulation der Ionenoptik

Ablaufschema: In Abb 7.9 ist das Ablauf-schema zur Simulation der Ionenoptik dar-gestellt. Es folgt die Erklärung der einzelnen Punkte entsprechend der Reihenfolge des Ab-laufschemas, wobei die Randbedingungen in Abschnitt 7.4.2 gesondert betrachtet werden.

Elektrostatisches Potential berechnen:

In dem verwendeten Modell wechselwirken die Elektronen nur über ihr elektrostati-sches Fernfeld mit den Ionen (keine Coulomb-Stöße). Dementsprechend ist es möglich, die Elektronen(dichte) indirekt über das berech-nete Potential zu berücksichtigen. Hierzu wird die Ladungsdichte ρ in der Poisson-Gleichung (4.39) mit ρ = ρ_i −ρ_e in die La-dungsdichte der Ionenρ_i und der Elektronen ρ_e separiert. Für letztere gilt ρ_e = n_ee, wo-bein_emit der Boltzmann-Beziehung 4.44 be-schrieben wird. Es folgt die nichtlineare par-tielle Differentialgleichung

∆φ_(~r) =−ρ_i−en_ese

e(^φ(~r)−φs)

kBTe

ε0

. (7.20) Zur Lösung wird, wie in Referenz [132] vor-geschlagen, die Finite-Differenzen-Methode in Kombination mit dem mehrdimensiona-len Newtonverfahren verwendet. So wird für eine gegebene Raumladungsdichte der Io-nen das Potential berechnet, welches intrin-sisch die Raumladungsdichte der Elektronen beinhaltet. Eine Beschreibung der Berech-nung erfolgt in Abschnitt B.4 bzw. Unterab-schnitt B.4.4. In der ersten Iteration des Ab-laufschemas ist die Raumladungsdichte der Ionen im gesamten Raum gleich 0.

Elektrostatisches Feld berechnen:

Nachdem das Potential auf den Stützpunk-ten berechnet ist, wird das elektrostatische Feld an diesen Punkten durch Anwendung

der Finite-Differenzen-Methode auf Gl. (4.38) ermittelt. Hierbei muss insbe-sondere der Umgang mit Dirichlet- und Neumann-Randbedingungen betrachtet werden. Eine detailliertere Beschreibung erfolgt im Anhang in Abschnitt B.4 bzw.

Unterabschnitt B.4.6.

97

Kapitel 7 Extraktionssystem

Ionentrajektorien berechnen: Bei der Initialisierung werden N Startteilchen für den Zeitschritt t = 0 festgelegt. Deren Verteilung in y- und z-Richtung sowie die Position in x-Richtung und deren Geschwindigkeit werden in Abschnitt 7.4.2 be-schrieben. Für jedes Startteilchen wird eine Trajektorie berechnet. Diese besteht aus mehreren Teilchen mit dem zeitlichen Abstand ∆t, wobei jedes Teilchen die Positi-on und Geschwindigkeit zum jeweiligen Zeitpunkt beinhaltet. Ein solches Teilchen repräsentiert eine bestimmte Anzahl an Ionen und wird dementsprechend auch als Makro-Teilchen bezeichnet. Masse und Ladung des Makroteilchens sind frei einstell-bar. Die einzelnen Teilchen einer Trajektorie sind sequenziell zu berechnen, jedoch erfolgt die Berechnung der Trajektorien auf mehreren CPU-Kernen getrennt. Der Ablauf jeder Trajektorienberechnung erfolgt ausgehend vom jeweiligen Startteilchen.

Dieses wird so lange um den Zeitschritt ∆t weiterbewegt, bis es bei der Kollision mit einem Extraktionsgitter oder mit der Weltraumgrenze gelöscht wird und somit die Trajektorienberechnung endet. Bei der Löschung werden Ladung und Impuls des Teilchens auf dem entsprechenden Objekt¹ deponiert. Die Weiterbewegung um einen Zeitschritt, inklusive der Behandlung der hierfür notwendigen Randbedingungen, ist in Abschnitt B.6.2 beschrieben.

Ionenladungsdichte berechnen: Nachdem alle Trajektorien simuliert wurden, liegt eine Anzahl an Teilchen im Raum vor. Der Prozess, um die an annähernd kon-tinuierlichen Positionen vorliegenden Punktladungen den Stützpunkten des diskret platzierten Berechnungsgitters zuzuordnen, wird als Scatter bezeichnet. Hierbei ist wichtig, das gleiche Schema wie beim Gather-Verfahren (Bestandteil der Teilchen-bewegung, siehe Abschnitt B.6.2) zu verwenden. Beim Scatter-Verfahren wird eine Schleife durchlaufen, in der für jedes einzelne Teilchen die Gewichtungsfaktoren wie bei der trilinearen Interpolation (siehe Abschnitt B.6.1) berechnet werden. Die La-dung eines jeden einzelnen Teilchens wird entsprechend der Gewichtungsfaktoren auf die 8 benachbarten Stützpunkte aufgeteilt. Der Stützpunkt der dem Teilchen am dichtesten ist bekommt den größten Anteil der Ladung zugeordnet.

Unterrelaxation: Bei der Berechnung der Ionenladungsdichte wird die Ladungs-dichte der Trajektorien zuerst in einem temporären Array ρ^Tr_i gespeichert. Dies ist notwendig, da die zu Beginn berechneten Trajektorien, gerade in der Nähe des Per-veanzlimits, also bei hohen Ionendichten als Randbedingung, sehr stark von der realen Lösung abweichen können. Würde die Ladung der berechneten Trajektorien in der nächsten Iteration verwendet, entstünden ggf. Zwischenergebnisse, auf deren Basis die Simulation nicht konvergiert. Deswegen wird ein Relaxationsparameter ω zur Unterrelaxation gemäß

ρ_i =ρ_i+ρ^Tr_i −ρ_iω (7.21) verwendet, wobeiρi das Array der tatsächlich verwendeten Ionenladungsdichte ist.

Er hat einen Wert zwischen 0 und 1 und sorgt dafür, dass die Ladungsdichte im Raum in jeder Iteration nur etwas korrigiert wird.

Abbruchkriterium: Für das Abbruchkriterium wird üblicherweise aus einem Systemzustand, wie z. B. der Ladungsdichte im Raum, ein skalarer Wert berechnet.

1 Dies sind zum einen die jeweiligen Gitter, aber auch die Randfläche zum Weltraum hin.

98

7.4 Simulation der Ionenoptik

Da das exakte Ergebnis nicht bekannt ist, wird nicht der Fehler, sondern die Änderung des Wertes zwischen den letzten beiden Iterationen betrachtet. Ist diese klein genug, wird von einer Konvergenz der Lösung mit geringer Abweichung von dem exakten Ergebnis ausgegangen und die Simulation beendet. Experimentell lässt sich ermitteln, dass selbst bei annähernd gleich bleibenden Gitterströmen und Raumladungen (Mittelwert über das gesamte Gebiet) noch Änderungen der Trajektorien bzw. der Ladungsdichte am Ende des Simulationsgebietes auftreten.

Deswegen wird als Abbruchkriterium bewertet, wie stark sich die Teilchenpositionen innerhalb einer Iteration geändert haben. Als Vergleichsgröße wird die mittlere quadratische PositionsabweichungR verwendet:

R=

v u u t

N−1

X

n=0 M−1

X

m=0

~r_[n,m,i]−~r_[n,m,i+1]² (7.22) Hierbei ist~r_[n,m,i] das m-te Teilchen der Trajektorie n in Iteration i. Zur Verringe-rung der erforderlichen Rechenzeit wird nur eine reduzierte Anzahl an Trajektorien verglichen.