A&A Title Image Startseite | Programm: wechselwirkende Galaxien | Algorithmen, Lektione, Analysen | Unsere Vorläufer | Kontakt, Datenschutz

Simulation wechselwirkender Galaxien

verfasst von Helmut Jahns (erschienen im VdS-Journal 12)

Einleitung

Lange Zeit wurde angenommen, daß Kollisionen zwischen Galaxien ähnlich wie Zusammenstöße von Sternen eher ausgesprochen unwahrscheinliche Ereignisse sind. Erst in der zweiten Hälfte des vergangenen Jahrhunderts ergaben Simulationen auf Großrechnern, daß unregelmäßig geformte Galaxien ihre Deformation größtenteils den Störungen verdanken, die benachbarte Galaxien in Form von Gezeitenkräften bei Passagen oder Kollisionen einst auf sie ausgeübt. Während für die ursprünglichen Arbeiten noch Supercomputer herangezogen werden mußten, sind derartige numerische Berechnungen inzwischen in Reichweite des heimischen PC?s gelangt. Dies ermöglichte die Realisierung eines einfachen kleinen Programms mit dem Namen Stellardynamik, welches verschiedene Erscheinungsformen wechselwirkender Galaxien veranschaulicht. Die Grundlagen und Vorgehensweise bei der Verwirklichung dieses Programmierprojekts möchte ich im folgenden schildern.

Das Modell

Die Wechselwirkungen zwischen Galaxien werden auf rein gravitative Kräfte aller Sterne untereinander zurückgeführt, die die in ihnen enthaltenen Sterne aufeinander ausüben. Dabei bleiben gasdynamische Effekte, die z.B. zur verstärkten Sternentstehung führen können, unberücksichtigt. Für alle Einzelsterne der Galaxien kennt man zu einem gegebenen Zeitpunkt ihre Positionen und Geschwindigkeiten. Damit lassen sich für jeden Einzelstern die Beschleunigungen, die von allen übrigen Sternen auf ihn ausgeübt werden, berechnen und aufsummieren (strenge n-Körper-Rechnung). Mittels der Gleichungen (hier am Beispiel der x-Koordinate)

rx, neu := rx + vx Δ t und
vx, neu := vx + ax Δ t

(r: Ortskoordinate, v: Geschwindigkeit, a: Beschleunigung) können für alle Sterne neue Koordinaten berechnet werden, die denen einer Momentaufnahme nach Verstreichen einer Zeitspanne Δ t, auch Schrittweite genannt, entsprechen. Werden diese Einzelschritte wie in dieser Simulation ständig wiederholt, so bezeichnet man dieses Vorgehen als iterativ.

Galaxien bestehen jedoch aus einigen zehn bis hundert Milliarden Sternen. Sie alle bei der Berechnung erfassen zu wollen ist derzeit unrealistisch und auch gar nicht zweckmäßig. Anstelle dessen wird angenommen, daß eine Galaxie aus einer erheblich kleineren aber hinreichenden Anzahl an Sternen, im folgenden Modellsterne genannt, besteht. Sie unterscheiden sich von reellen Stern darin, daß ihnen eine entsprechend höhere Masse zugeschrieben wird. Die Gesamtmasse einer Modellgalaxie entspricht somit der Masse einer realen Galaxie.

Die Modellgalaxien werden als zweidimensionale Scheibengalaxien angenommen. In der derzeitigen Programmversion werden die innere Struktur (Spiralarme, Halo, Bulge) vernachlässigt. Die radiale Massenverteilung, also die vom Galaxienzentrum nach außen gehende Massendichte, wurde der einer Spiralgalaxie nachempfunden, d.h. einem Exponentialgesetz der Form e-ar (r: Abstand vom Galaxienzentrum) zugrunde gelegt.

Algorithmus

Eines der großen Probleme solcher Simulationen ist die Teilchenzahl: Die Rechenzeit bei der strengen n-Körper-Rechnung wächst proportional mit dem Quadrat der Teilchenzahl. Glücklicherweise gibt es Verfahren, mit deren Hilfe es möglich ist, das Ansteigen des Rechenaufwands zu begrenzen. Zum einen kann man als Algorithmus die eingeschränkte n-Körper-Rechnung wählen: anstatt sämtliche Wechselwirkungen der Sterne untereinander zu berücksichtigen werden bei diesem Verfahren nur Wechselwirkungen der Sterne mit den Galaxienzentren betrachtet. Der Rechenaufwand steigt hierbei nur proportional zur Teilchenzahl selbst. Dieses Verfahren ermöglicht einen sehr schnellen Code und wird von den meisten der bisherigen Tools eingesetzt. Es hat jedoch einen entscheidenden Nachteil: Gezeiteneffekte lassen sich damit zwar erzeugen, doch Verschmelzungen von Galaxien können nicht modelliert werden. Wenn z.B. zwei Galaxien kollidieren, so sind ihre Galaxienzentren die einzigen kraftausübenden Körper im Ensemble, so daß sie sich wie beim Zwei-Körper-Problem auf Kepler-Bahnen umeinander bewegen und niemals miteinander verschmelzen.

Ziel ist es, die Exaktheit der strengen n-Körper-Rechnung mit dem schnellen Algorithmus des eingeschränkten Problems zu kombinieren. Ein sehr guter Kompromiß ist der Treecode von Barnes und Hut, welcher in [1] beschrieben wurde und im folgenden kurz skizziert werden soll: Grundgedanke ist, daß für die Berechnung der Beschleunigung auf einen Körper entferntere Masseobjekte nur kleinere Beiträge leisten und somit zu größeren Punktmassen, sogenannten Pseudomassen, zusammengefaßt werden können, wobei der dabei entstehende Fehler klein bleibt.

Wie erreicht man dies in der Praxis? Die Pseudo- und Einzelmassen müssen hierarchisch gegliedert werden, so daß entschieden werden kann, ob die kraftausübenden Massen hinreichend weit entfernt sind, daß Pseudomassen verwendet werden können oder ob man eine feinere Massenverteilung wählen muß.


Abb. 1 Einteilung der quaderförmigen Zellen in Unterzellen. Nähere Erläuterung siehe Text.

Zu diesem Zweck stelle man sich eine Zelle in Form eines Quaders vor, die das gesamte Ensemble beinhaltet (s. Abb. 1). Diese Zelle wird gleichmäßig in acht Unterzellen unterteilt (die Kantenlänge der Unterzellen ist genau halb so groß wie die der übergeordneten), welche ihrerseits eine Zusammenfassung der in ihnen enthaltenen Punktmassen repräsentieren. Diese Unterteilungen werden sukzessive fortgesetzt, bis in einer Zelle nur noch eine Einzelmasse enthalten ist. Innerhalb jeder Zelle befindet sich eine Pseudomasse, die sämtliche in der Zelle enthaltenen Einzelmassen repräsentiert. Die Position der Pseudomasse wird durch Bildung des Massenschwerpunktes aus allen enthaltenen Einzelmassen berechnet, während ihre Masse einfach die Summe der Einzelmassen ist. Insbesondere ist jede Einzelmasse auf allen übergeordneten Hierarchieebenen vertreten.

In der Programmierung wird dieser hierarchische Aufbau am besten durch eine Baumstruktur realisiert (vgl. Abb 2). Jeder Knoten innerhalb des Baumes repräsentiert einen Datensatz zu einer Pseudo- oder Einzelmasse sowie bis zu acht Verweise auf Speicherbereiche, in denen jeweils untergeordnete Datensätze abgelegt sind. Solche Baumstrukturen zählen zu den Standardalgorithmen in der Programmierung [3] und haben den Vorteil, daß sie dynamisch sind. Es lassen sich auf die Art beliebig viele massebehaftete Objekte anordnen und verwalten.


Abb. 2 Baumartige Datenstruktur, wie sie in dem Programm verwendet wird (s. Text). Aus Gründen der Übersichtlichkeit wurden in diesem Diagramm zu jedem Knoten nur zwei statt der bis zu acht Nachfolger dargestellt.

Ein sinnvolles Kriterium, ob mit der Gravitationskraft der momentan betrachteten Pseudomasse gerechnet werden kann, ist das Verhältnis der Kantenlänge l der Zelle zum Abstand d des Zellmittelpunktes von derjenigen Einzelmasse, deren Beschleunigung ausgerechnet werden soll. Zu diesem Zweck geht man für alle Einzelmassen rekursiv durch die Baumstruktur und berechnet den Abstand d des Zellmittelpunktes. Ist das Verhältnis von d zu l größer als ein vorgegebener Wert Δ, so kann mit der aktuellen Pseudomasse gearbeitet werden, ansonsten wird eine Ebene tiefer verzweigt.

In der Abb. 3 ist die Ausführungszeit einer Simulation mit dem Treecode derjenigen der strengen n-Körper-Rechnung gegenübergestellt. Es zeigt sich hier sehr deutlich, daß der Performancevorteil des Treecode um so größer ist, je mehr Punktmassen an der Rechnung beteiligt sind. Mit dem Treecode benötigt eine Simulation mit 1000 Punktmassen auf einem 1,46-GHz-AMD-Rechner nur wenige Minuten.
Abb. 3 Rechenzeiten für die strenge n-Körper-Rechnung (hinten) und den Treecode (vorn) für verschiedene Partikelzahlen.

Abb. 4 Sequenz einer Wechselwirkung zweier Galaxien.

Implementation

Ein solches Projekt beginnt man nicht damit, daß man die Idee dazu bekommt und sofort mit dem Codieren beginnt. Eine solche Vorgehensweise führt unweigerlich zu Code, der sich nur schwer warten und erweitern läßt. Ich habe mir vielmehr zu Anfang detailliert überlegt, welche Funktionen die Software beherrschen muß. Erst als der Funktionsumfang festgelegt war, wurde das Klassenkonzept erstellt. Ich habe mich für die Entwicklung mit dem Borland C++ Builder 5.0 Standard entschieden, also auf eine objektorientierte Programmierweise. Besonderes Augenmerk habe ich auf die Trennung von Klassen der betriebssystemabhängigen Benutzeroberflächen und der in reinem ANSI-C++ geschriebenen Programmlogik gelegt. Dies erleichtert die Portierbarkeit des Programms auf andere Plattformen erheblich.

Programmbeschreibung

Das Programm wurde unter dem Gesichtspunkt entwickelt, dem Anwender eine komfortable und zweckmäßige Benutzerschnittstelle zu bieten. Z.B. kann das zu untersuchende Szenario aus Modellgalaxien vom Anwender frei nach seinen Vorstellungen konfiguriert und gespeichert werden. Da das Programmenü während der Berechnung bedienbar bleibt (Multithreading), kann die Simulation jederzeit unterbrochen oder der sichtbare Ausschnitt gezoomt werden. Die Steuerelemente sind genauso wie die Hilfe wahlweise auf Deutsch oder Englisch.

Die Simulation wird fortwährend als Animation aufgezeichnet. Nach der Beendigung der Berechnung kann der Benutzer in den Animationsmodus wechseln und sich das Verhalten seines Galaxienensembles im Zeitraffertempo betrachten oder sich einen besonders interessanten Zeitabschnitt als Standbild auswählen.

Stellardynamik läuft ab Windows 95/NT 4.0 und kann hier (Freeware) heruntergeladen werden.

Literatur