Im Rahmen des Projektlabor Mathesis, welches im MINTgrün Orientierungsstudium angesiedelt ist, wurde dieses Projekt umgesetzt. Wie bereits der Name andeutet, ist das Ziel des Projekts eine Simulation von Partikeln im dreidimensionalen Raum. Besonders auf die drei Kerneigenschaften Gravitation, Reibung und Kollision wird hierbei wertgelegt.
Während der Arbeit an unserem Projekt hielten wir unsere Fortschritte in diesem Protokoll fest. Der überwiegende Teil des Projekts entstand beim gemeinsamen Arbeiten in der Gruppe.

Unser Programm basiert auf der Bibliothek vpython, welche die graphische Umsetzung ermöglicht. Im Anfang unseres Programms werden die Partikel zweimal erzeugt. Einmal als die Liste „partikel“, welche die einzelnen Partikel als Dictionaries enthält, und einmal als die Liste „kugeln“, welche vpython Objekte enthält, die graphisch dargestellt werden können.

Nach dem definieren einiger Funktionen, auf welche im Detail weiter unten eingegangen wird, findet sich die Hauptschleife in Form eines while Loops, eine Schleife, welche theoretisch unendlich laufen kann. Die Grundidee hierbei ist, dass jede Iteration in dieser Hauptschleife einen Frame der Simulation bildet. Die Dictionaries aus „partikel“ enthalten sieben Einträge. Hierbei stellen drei die Position im dreidimensionalen Raum dar und- drei weitere Einträge den Geschwindigkeitsvektor. Der letzte Eintrag ist für die Masse, welche zum Beispiel bei der Gravitationsfunktion benötigt wird. In jeder Iteration wird dann die Position eines jeden Partikels aktualisiert, indem auf die alte Position $\vec{p}_{alt}$ der Geschwindigkeitsvektor $\vec{v}$ hinzuaddiert wird.
Es gilt also für die neue Position $\vec{p}_{neu}=\vec{p}_{alt}+\vec{v}$. Auf diese Weise würde zunächst eine gradlinige Bewegung für alle Partikel entstehen. Die Simulation wird aber dadurch interessant, dass in der Hauptschleife vorneweg noch die anderen Funktionen die Geschwindigkeitsvektoren verändern, teilweise entfallen sogar ganze Partikel. Im folgendem werden diese Kerneigenschaften und die damit verbundenen Funktionen erläutert, welche den eigentliche Grund der Simulation darstellen.

Wie oben erklärt besteht unser Algorithmus aus drei Hauptfunktionen, diese werden im folgendem genauer beschrieben.

Die Gravitation war die erste und bleibt nach wie vor die wichtigste programmierte Eigenschaft in unserem Code.

$F=m\cdot a$

Als Gravitation bezeichnet man die Anziehungskraft zwischen Massen. Jede Masse im Universum zieht sich gegenseitig an, $F_G=\gamma\cdot\frac{m_1\cdot m_2}{r^2}$ dabei ist die Stärke der Anziehungskraft abhängig von Masse und Abstand der Massekörper. Die Abhängigkeit vom Abstand zwischen den Körpern führt dazu, dass bei sehr geringen Abständen (limes r gegen null) die Gravitationskraft sehr groß wird. Da Körper ihre Masse nicht in einen Punkt konzentrieren, dreidimensionale Objekte sind, kommt es zu Kollisionen, bevor die Stärke der Gravitation sozusagen ins unendliche schießt. [Zur Kollision in unserer Simulation folgt später ein Abschnitt] Die Gravitationskonstante $\gamma$ beträgt im Universum ungefähr $6,67\cdot10^{-11}\frac{m^3}{kg\cdot s^2}$, weshalb die Gravitation zwischen kleinen Massen üblicherweise nicht zu einer spürbaren Beschleunigung führt. In unserer Simulation haben wir größere Konstanten gewählt, und somit die Stärke der Anziehungskraft deutlich erhöht.

Bei den Berechnungen stoßen wir, sobald mehr als zwei Massekörper im Spiel sind, auf das sogenannte Dreikörperproblem ¹⁾ : für drei und mehr Körper ist es bisher (bis auf Spezialfälle) nicht gelungen, exakte analytische Lösungen zu finden und Bahnkurven genau vorherzusagen, die Differentialgleichungen beinhalten zu viele unbekannte Variablen. Numerische Berechnungen sind dennoch möglich, wobei für jedes Paar von Massekörpern in jedem Zeitschritt die Kraft bestimmt -und aus allen Kräften für die einzelnen Partikel die resultierende Beschleunigung aufsummiert wird. Mit dieser neuen Beschleunigung ergibt sich eine aktuelle Momentangeschwindigkeit für jedes Partikel, und die Position im nächsten Zeitschritt wird für die erneute Berechnung verwendet. Mithilfe dieser Methode konnten wir in unserer Simulation die Wechselwirkung vieler Massen aufeinander darstellen.

In jedem Frame, also Zeitschritt, unserer Simulation, werden die einzelnen Partikel nacheinander betrachtet. Dieses Schema zeigt grob die durchlaufenen Schleifen:

Zunächst werden die Abstände zu den anderen Partikeln berechnet, und anschließend entweder eines der beiden Kollisions-szenarien (Fusion/Explosion), oder die Gravitations-beschleunigung berechnet. Für die Berechnung der Beschleunigung erstellten wir die Funktion „Gravitation“. Als Input nimmt sie die Positionen zweier Partikel (diese sind als x,y,z Koordinate in der Partikelliste gegeben) und die Stärke der Gravitationskonstante, und als Output erhalten wir die Beschleunigung, welche auf Partikel i wirkt. In einem Zeitschritt von einer Sekunde entspricht diese Beschleunigung dem Geschwindigkeitsvektor, welcher durch diese Beschleunigung an dem Partikel ansetzt, da $a=\frac{v}{s}$ gilt. Diese Geschwindigkeitsvektoren werden alle auf die vorherige Geschwindigkeit von Partikel i hinzu addiert, und es ergibt sich die neue Momentangeschwindigkeit.

Aus physikalischer Sicht steht jeder Frame für eine Sekunde, wobei die Simulation (aus ästhetischen und praktischen Gründen) deutlich schneller abgespielt wird.

Abbildung 1: Die Umsetzung der Gravitation und der Reibung

Als einen weiteren Bestandteil unserer Simulation haben wir Reibung eingebaut, welche die Bewegung dämpft, da kinetische Energie umgewandelt wird. Es gibt verschiedene Arten von Reibung, Abhängig von beispielsweise der größe der Angriffsfläche des Objekts und der Beschaffenheit der Umgebung durch die es sich bewegt. Wir haben uns dafür entschieden, abhängig von Geschwindigkeit und Masse der Partikel, einen geringen Geschwindigkeitsverlust zu kalkulieren.

Wenn sich also einige Partikel in einem Kräftegleichgewicht auf stabilen Bahnen bewegen, führt die Reibung in der Simulation dazu, dass sie stetig kleine Anteile ihrer Bewegungsenergie verlieren, und die Bahnen nicht endlos als stabile Bewegung gegeben bleiben.

Die Kollision findet genau dann statt, wenn sich zwei Partikel zu nahe kommen, hierbei wird zwischen der Fusion bei langsamen -und Explosion bei schnellen Zusammenstößen unterschieden. Die besondere Herausforderung hierbei bestand darin, dass sich die eigentliche Partikelanzahl verändert. Ein Problem ähnlich wie bei der Gravitation besteht darin, dass auch hier zwischen je zwei Partikeln zunächst der Abstand ermittelt werden muss, wodurch eine hohe Rechenzeit entsteht. Die beiden folgenden Fälle treten nur ein, wenn der Abstand zwischen zwei Partikeln klein genug ist.

Anschaulich gesprochen, verschmelzen bei der Fusion zwei Partikel zu einem. Im Programm wird hierbei zunächst ein Partikel übernommen. Seine Position bleibt gleich, da sich beide Partikel relativ nah sein sollen. Als nächstes wird der neue Richtungsvektor $\vec{v}_{neu}$ ermittelt. Da hier der Impulserhaltungssatz wichtig ist, verwenden wir folgende Formel:
$\vec{v}_{neu}\cdot(m_1+m_2)=m_1\cdot\vec{v}_1+m_2\cdot\vec{v}_2.$
Hier stehen $m_1$ und $m_2$ für die Massen der zwei Partikel und $\vec{v}_1$ und $\vec{v}_2$ für die entsprechenden alten Geschwindigkeitsvektoren. Zum Schluss wird als neue Masse die Summe der beiden alten verwendet und dass übergebliebene Partikel entfernt und die Fusion ist abgeschlossen.

Die Idee hierbei war, dass zwei schnell aufeinander prallende Partikel nicht fusionieren, sondern eher auseinander gesprengt werden. Hierbei wird eine Anzahl an Splitterpartikeln erzeugt, so dass später keine kinetische Energie verloren geht. In der Variabel „v“ wird hierfür zunächst die Geschwindigkeit der einzelnen Splitterpartikel, wir haben angenommen, dass alle neuen Partikel die gleiche Geschwindigkeit und Masse besitzen. Diese neuen Splitterpartikel werden nun erzeugt, zunächst nur eine Hälfte an dem Ort der Kollision mit einem zufälligem Richtungsvektor und dann die andere Hälfte, so dass je zwei Splitterpartikel entgegengesetzte Richtungsvektoren besitzen. Die Splitterpartikelmasse ist hierbei die Anzahl der Splitterpartikel geteilt durch die Summe der Massen der beiden alten Partikel. Zu guter Letzt werden alle Richtungsvektoren so skaliert, dass die Norm der Richtungsvektoren gleich zu „v“ ist. Dann können die alten Partikel entfernt und die neuen Partikel der Simulation hinzugefügt werden.

Um zu entscheiden, welche der beiden Vorgänge bei einer Kollision passieren, nehmen wir zunächst die Norm der Differenz der Richtungsvektoren der beiden alten Partikel. Wenn nun dieses Ergebnis kleiner als eine Konstante namens „threshold“ ist, entsteht eine Fusion, sonst eine Explosion. Es ist schwierig diesen Grenzwert so zu wählen, dass weder alles Explodiert, noch alles in sich zusammenfällt.

Die erste Herausforderung dieses Projektes war eine passende Bibliothek zu finden, um die graphische Umsetzung zu ermöglichen. Anfangs wollten wir matplotlib nutzen, aber schnell war klar, dass diese Bibliothek nicht unseren Vorstellungen entsprach. Alles funktionierte nur sehr langsam und war zuzüglich sehr umständlich. Wir entschieden uns für die Bibliothek vpyhthon, welche uns auf eine einfache Weise ermöglichte, Körper in einer dreidimensionalen Welt zu platzieren. In dieser Welt kann man sich sogar bewegen, die Simulation konnten wir uns so von allen Seiten anschauen und waren nicht auf einen Blickwinkel beschränkt. Als nächstes mussten wir eine sinnvolle Struktur finden, um die Partikel im Code zu repräsentieren. Wir entschieden uns für die Dictionaries, da mit diesem Datentyp verständlich auf Eigenschaften der Partikel zugegriffen werden kann. Wir begannen die Simulation, indem wir zunächst die Hauptschleife sehr einfach aufbauten, der einzige Vorgang war hierbei die Aktualisierung der Position durch das Hinzurechnen der Geschwindigkeitsvektoren bei allen Partikeln. Anschließend setzten wir die Gravitation und die Reibung um. Bei beiden Eigenschaften fingen wir mit physikalischen Vorüberlegungen an und erstellten hinterher die entsprechenden Funktionen mit der passenden Implementierung innerhalb der Hauptschleife. Schwierigkeiten bereitete uns auch die Kollision. Zunächst suchten wir eine Möglichkeit Objekte aus unserer Simulation verschwinden zu lassen, vpython bot entsprechende Möglichkeiten an.

Abbildung 2: Hier wurden die Partikel noch nicht entfernt

Bei der eigentlichen Programmierung bestand die Herausforderung darin, mit entfernten Objekten umzugehen. Anders als bei zum Beispiel der Gravitation konnten wir nun nicht mehr zwei for-Schleifen in einander verschachteln, um mögliche Kollisionen festzustellen und dann auch auszuführen, da sobald Partikel aus der Kollision entfernt wurden, die Schleifen frühzeitig abgebrochen werden konnten. Wir nutzten zwei while-Schleifen um diese Problematik zu umgehen. Die Fusion ließ sich recht problemlos umsetzten, aber die Explosion bereitete uns Schwierigkeiten und bot eine Menge Fehlerquellen auf.

Abbildung 3: Das Ergebnis eines Fehlers auf dem Weg zur Explosion

Schlussendlich konnten wir diese Punkte aber umsetzen. Man könnte das Projekt noch gut fortführen, indem bespielsweise noch interaktive Funktionen hinzugefügt werden, um Konstanten der Simulation, wie zum Beispiel bei der Gravitation, zu verändern. Auch die Codestruktur kann noch verbessert werden, indem zum Bespiel Vektoren zu einem Objekt zusammengefasst werden. Unser Projekt bot eine Menge Potential für kreative Ideen und wir konnten auch einige darin verwirklichen.

Damit das Programm funktioniert, muss die Bibliothek vpython installiert sein. Es sollte sich beim Ausführen ein Browserfenster mit einem Bild öffnen. Dieses Bild kann vergrößert werden, indem man die Ränder des Bildes nach außen zieht. Daher startet die Simulation erst nach fünf Sekunden. In der Simulation bewegen kann man sich mit (Shift/Strg/Alt)+Linksklick für Bewegung, Rotation und Zoom.

'''
Projekt Partikelsimulation aus dem Wintersemester 2020/21 im Rahmen des Mathesis Projektlabors von
Paulina Marie Fritze,
Berhan Rêzan Isik,
Moritz Merkel und
Amelie Schäfer
'''
import random as rd
import vpython as vp
import time
#Konstanten
height=width=length=3000
n=35
 
partikel=[] #{x,y,(z),dx,dy,(dz),m}
for i in range(n):
    partikel.append({'x': rd.randint(-width,width),'y': rd.randint(-height,height),'z': rd.randint(-length,length),'dx': 6*rd.random()-3,'dy': 6*rd.random()-3,'dz': 6*rd.random()-3,'m': rd.random()*1000})
 
kugeln=[vp.sphere(pos=vp.vector(partikel[i]['x'],partikel[i]['y'],partikel[i]['z']),radius=20) for i in range(n)]
time.sleep(5)
 
def norm(x,y,z):
    return (x**2+y**2+z**2)**0.5
    #alternative Metrik 
    #return abs(x)+abs(y)+abs(z)
    #supremumsnorm
    #return max(abs(x),abs(y),abs(z)
 
    #F=G*m1*m2/(r^2)
def Gravitation(p1, p2, G=70):
#punkt + einheitsrichtungsvektor * beschleunigung
    dx,dy,dz=(p2['x']-p1['x']),(p2['y']-p1['y']),(p2['z']-p1['z'])
    d=norm(dx,dy,dz)
    if d==0:
        return (0,0,0)
    F=G*p1['m']*p2['m']/d**2
    a=F/p1['m']
    dx,dy,dz=a*dx/d,a*dy/d,a*dz/d
    return (dx,dy,dz)
def Reibung(p1,c=-0.001):
    return (c*p1['dx']/p1['m'],c*p1['dy']/p1['m'],c*p1['dz']/p1['m'])
 
#main loop
while True:
    vp.rate(60)
    i=0
    j=0
    gesamtsplitteranzahl=0
    while i < len(partikel)-gesamtsplitteranzahl:
        j=i+1
        while j < len(partikel)-gesamtsplitteranzahl:
            d=norm(partikel[j]['x']-partikel[i]['x'],partikel[j]['y']-partikel[i]['y'],partikel[j]['z']-partikel[i]['z'])
            if d<(kugeln[i].radius+kugeln[j].radius)/2:
                differenz=norm(partikel[i]['dx']-partikel[j]['dx'],partikel[i]['dy']-partikel[j]['dy'],partikel[i]['dz']-partikel[j]['dz'])
                treshold=80
                if differenz<treshold or len(partikel)>150:
                    #Fusion
                    partikel[i]['dx']=(partikel[i]['m']*partikel[i]['dx']+partikel[j]['m']*partikel[j]['dx'])/(partikel[i]['m']+partikel[j]['m'])
                    partikel[i]['dy']=(partikel[i]['m']*partikel[i]['dy']+partikel[j]['m']*partikel[j]['dy'])/(partikel[i]['m']+partikel[j]['m'])
                    partikel[i]['dz']=(partikel[i]['m']*partikel[i]['dz']+partikel[j]['m']*partikel[j]['dz'])/(partikel[i]['m']+partikel[j]['m'])
                    partikel[i]['m']+=partikel[j]['m']
                    partikel.pop(j)
                    kugeln[j].visible=False
                    kugeln.pop(j)
                else:
                    #Explosion
                    splitteranzahl=4*2#gerade zahl
                    masse=(partikel[i]['m']+partikel[j]['m'])/splitteranzahl
                    #splitteranzahl*m*v^2=m1*v1^2+m2*v2^2=c
                    v=((partikel[i]['m']*norm(partikel[i]['dx'],partikel[i]['dy'],partikel[i]['dz'])**2+partikel[j]['m']*norm(partikel[j]['dx'],partikel[j]['dy'],partikel[j]['dz'])**2)/(splitteranzahl*masse))**0.5
                    splitterpartikel=[{'x':0,'y':0,'z':0,'dx':2*rd.random()-1,'dy':2*rd.random()-1,'dz':2*rd.random()-1,'m':0} for k in range(int(splitteranzahl/2))]
                    splitterpartikel=[{'x':partikel[i]['x'],'y':partikel[i]['y'],'z':partikel[i]['z'],'dx':splitterpartikel[k]['dx']*v/norm(splitterpartikel[k]['dx'],splitterpartikel[k]['dy'],splitterpartikel[k]['dz']),                                                                                                                                                                                                                                                                         'dy':2*splitterpartikel[k]['dy']*v/norm(splitterpartikel[k]['dx'],splitterpartikel[k]['dy'],splitterpartikel[k]['dz']),'dz':splitterpartikel[k]['dx']*v/norm(splitterpartikel[k]['dz'],splitterpartikel[k]['dy'],splitterpartikel[k]['dz']),'m':masse} for k in range(int(splitteranzahl/2))]
                    splitterpartikel+=[{'x':partikel[i]['x'],'y':partikel[i]['y'],'z':partikel[i]['z'],'dx':-splitterpartikel[k]['dx'],'dy':-splitterpartikel[k]['dy'],'dz':-splitterpartikel[k]['dx'],'m':masse} for k in range(int(splitteranzahl/2))]
                    for k in range(splitteranzahl):
                        splitterpartikel[k]['dx']+=(partikel[i]['m']*partikel[i]['dx']+partikel[j]['m']*partikel[j]['dx'])/(partikel[i]['m']+partikel[j]['m'])
                        splitterpartikel[k]['dy']+=(partikel[i]['m']*partikel[i]['dy']+partikel[j]['m']*partikel[j]['dy'])/(partikel[i]['m']+partikel[j]['m'])
                        splitterpartikel[k]['dz']+=(partikel[i]['m']*partikel[i]['dz']+partikel[j]['m']*partikel[j]['dz'])/(partikel[i]['m']+partikel[j]['m'])
                    partikel.pop(j)
                    kugeln[j].visible=False
                    kugeln.pop(j)
                    partikel.pop(i)
                    kugeln[i].visible=False
                    kugeln.pop(i)
 
                    gesamtsplitteranzahl+=splitteranzahl
                    partikel+=splitterpartikel
                    kugeln+=[vp.sphere(pos=vp.vector(splitterpartikel[k]['x'],splitterpartikel[k]['y'],splitterpartikel[k]['z']),radius=20) for k in range(splitteranzahl)]
            else:
                j+=1
        i+=1
 
    #Geschwindigkeit aktualisieren
    #Gravitation
    for i in range(len(partikel)):#Hauptpartikel
        for j in range(len(partikel)):
            if i==j:
                continue
            d=Gravitation(partikel[i],partikel[j])#(dx,dy)
            partikel[i]['dx']+=d[0]
            partikel[i]['dy']+=d[1]
            partikel[i]['dz']+=d[2]
 
    for i in range(len(partikel)):
        d=Reibung(partikel[i])
        partikel[i]['dx']+=d[0]
        partikel[i]['dy']+=d[1]
        partikel[i]['dz']+=d[2]
    #Position aktualisieren
 
    for i in range(len(partikel)):
        partikel[i]['x']+=partikel[i]['dx']
        partikel[i]['y']+=partikel[i]['dy']
        partikel[i]['z']+=partikel[i]['dz']
    #vpython
    for i in range(len(partikel)):
        kugeln[i].radius=15*partikel[i]['m']**(1/3)
        kugeln[i].pos=vp.vector(partikel[i]['x'],partikel[i]['y'],partikel[i]['z'])

Zu guter letzt ist hier der Code und die erforderlichen Anforderungen für die Bibliotheken zusammengefasst. partikelsimulation.zip