Vielteilchensimulation mit anziehenden Kräften zwischen den Teilchen

Dieses Projekt erzeugt eine Simulation von Teilchenbahnen auf Grund der Wechselwirkungen(hier Gravitationskräfte) die zwischen den einzelnen Teilchen wirken. Zusätzlich wird die Stabilität der Bahnen untersucht und mit Kegelschnitten verglichen.

Physikalische Grundlagen(ident mit Projekt: "Simulation unseres Sonnensystems")

Damit die Teilchenbahnen berechnet werden können, muss folgendes Differentialgleichungssystem gelöst werden.

Anfangsbedingungen

• Startposition von n Planeten ([math]\vec r_1,\;\vec r_2,\;\cdots\;\vec r_n[/math])
• Startgeschwindigkeiten von n Planeten ([math]\vec v_1,\;\vec v_2,\;\cdots\;\vec v_n[/math])

Differentialgleichungen

Geschwindigkeiten
[math]\vec \dot{r_1}\;=\;\vec v_1[/math]
[math]\vec \dot{r_2}\;=\;\vec v_2[/math]
[math]\cdots\;[/math]
[math]\vec \dot{r_n}\;=\;\vec v_n[/math]

Beschleunigungen
[math]\vec \dot{v_1}\;=\;\vec a_1[/math]
[math]\vec \dot{v_2}\;=\;\vec a_2[/math]
[math]\cdots\;[/math]
[math]\vec \dot{v_n}\;=\;\vec a_n[/math]

Die Teilbeschleunigungen, die ein Planet durch die Gravitationsfelder aller anderen Planeten erfährt, werden mit Hilfe des NEWTONschen Gravitationsgesetzes, das die Anziehungskraft zweier Massepunkte bzw. zweier homogener Kugeln, wie es Planeten in diesem Modell näherungsweise sind, berechnet:

[math]\vec F_n\;=\;-\gamma\;\frac{m_nm_2}{r^2}\frac{\vec r}{\begin{Vmatrix} r \end{Vmatrix}} [/math]    und    [math]\vec F_n\;=\;m_n\vec a_n [/math] [math]\Longrightarrow \vec a_n\;=\;-\gamma\;\frac{m_2}{r^2}\frac{\vec r}{\begin{Vmatrix} r \end{Vmatrix}} [/math]

[math]\vec F_n\;[/math]	[math]\cdots\;[/math]	[math]\mbox{Gravitationskraft, die auf Planeten }n\mbox{ wirkt}\;[/math]
[math]m_n\;[/math]	[math]\cdots\;[/math]	[math]\mbox{Masse des Planeten }n\;[/math]
[math]m_2\;[/math]	[math]\cdots\;[/math]	[math]\mbox{Masse des weiteren Planeten zu dem die Gravitationskraft gerichtet ist}\;[/math]
[math]r\;[/math]	[math]\cdots\;[/math]	[math]\mbox{Abstand der beiden Planeten}\;[/math]
[math]\vec a_n\;[/math]	[math]\cdots\;[/math]	[math]\mbox{Gravitationsbeschleunigung des Planeten }n\;[/math]
[math]\gamma\;=\;6.67\;10^{-11}\;[/math]	[math]\cdots\;[/math]	[math]\mbox{Gravitationskonstante}\;[/math]

Um zur Gesamtbeschleunigung zu gelangen, werden alle Teilbeschleunigungen vektoriell addiert: [math]\vec a_{ges}\;=\;\vec a_1\;+\;\vec a_2\;+\;\cdots\;+\;\vec a_n[/math]

Programmtechnische Hintergründe

Ein wesentlicher Part dieses Projektes ist die Betrachtung der Stabilität der Bahnen. Hier werden die Ergebnisse der Lösung des DG-Systems nach der manuell programmierten simplen Eulermethode und dem matlabinternen Ode23t-Solver mit einem Kegelschnitt verglichen. Ode23t wird dann auch in einem weiteren Programm zur Berechnung der Teilchenbahnen unter Berücksichtigung aller Gravitations-Wechselwirkungen verwendet.

Beschreibung der Files

planetensystem_stabilitaetstest.m nimmt als Default die Mondbahn und zeigt die Stabilität(anhand Vergleich mit Kegelschnitt gefittet aus der ersten Umrundung) der Bahn. Oder es kann aus den 2 Fällen: 'Ellipse' und 'Hyperbel' gewählt werden. Zusätzlich werden einige dinge wie Form der Bahn(entweder Ellipse oder Hyperbel(Kreis bzw Gerade sind Grenzfälle dazu)), falls vorhanden, Umlaufzeit, kleinste- mittlere grösste Annäherung/Geschwindigkeit ausgegeben.

planetensystem_nteilchenwechselwirkung.m ist die Realisierung der Simulation mit wählbarer Anzahl an Himmelskörpern. Entfernung von der Sonne in etwa im Intervall der innersten vier Planeten. Programmtechnisch kann man sich entweder entscheiden die Darstellung gleichzeitig mit der Berechnung auszuführen, oder eben zurerst die ganze Berechnung durchzuführen und am Ende die Simulation zu starten. Beim starten des Programms kann man diese Wahl treffen.

intacc.m Funktion für Odesolver im Programm planetensystem_stabilitaetstest.m

intacc_nteilchen.m Funktion für Odesolver im Programm planetensystem_nteilchenwechselwirkung.m

Beschreibung einiger Ergebnisse

Das Phasendiagramm(x in Abhängigkeit von vx) zeigt im Fall der Lösung mit Ode23t eine geschlossene Kurve die bei der periodischen Bewegung(Ellipsenbahn) auftreten sollte.

>>>Bild, planetensystem_stabilitaetstest.m; 'ellipse'; Figure No. 2<<<

Hingegen zeigt das Phasendiagramm im Fall der Lösung mit der Eulermethode keine geschlossene Kurve. (die Schittweite wurde demonstrativ zu gross gewählt, jedoch würde eine kleinere Schrittweite bei der Eulermethode ohnehin verhältnissmäßig lange Rechenzeit erfordern

>>>Bild, planetensystem_stabilitaetstest.m; 'ellipse'; Figure No. 3<<<

Hier der Vergleich mit einem Kegelschnitt,

Im Fall der Ellipse:

>>>Bild, planetensystem_stabilitaetstest.m; 'ellipse'; Figure No. 1<<<

Im Fall einer Hyperbel:

>>>Bild, planetensystem_stabilitaetstest.m; 'hyperbel'; Figure No. 1<<<

Hier sieht man das die Bahn mit Ode23t über den gewählten Zeitraum(10 Umrundungen) sehr schön auf dem Kegelschnitt bleibt, mit Eulermethode jedoch mehr oder weniger davon abweicht. (blaue Kreuze Ode23t, schwarze Kreuze Euler)