Übung: 2D-Orbitalmechanik-Simulation (Python)

Question

Übung: 2D-Orbitalmechanik-Simulation (Python)

Orbit
Kosmos
Schwere
n-Körper-Simulationen

Zustandsraum

Nur ein kleiner Disclaimer vorweg: Ich habe nie Astronomie oder irgendwelche exakten Wissenschaften (nicht einmal IT) studiert, also versuche ich, diese Lücke durch Selbststudium zu füllen. Astronomie ist einer der Bereiche, die meine Aufmerksamkeit erregt haben, und meine Vorstellung von Selbstbildung ist ein direkt angewandter Ansatz. Also, gleich auf den Punkt - dies ist ein Orbital-Simulationsmodell, an dem ich beiläufig arbeite, wenn ich Zeit / Laune habe. Mein Hauptziel ist es, ein komplettes Sonnensystem in Bewegung zu setzen und Raumschiffstarts zu anderen Planeten zu planen.

Es steht Ihnen allen frei, dieses Projekt jederzeit wieder aufzunehmen und viel Spaß beim Experimentieren!

aktualisieren!!! (November 10)

Die Geschwindigkeit ist jetzt das richtige deltaV und bei zusätzlicher Bewegung wird jetzt der Summenvektor der Geschwindigkeit berechnet
Sie können so viele statische Objekte platzieren, wie Sie möchten, auf jeder Zeiteinheit Objekt in Bewegung prüft auf Gravitationsvektoren aus allen Quellen (und prüft auf Kollision)
die Leistung von Berechnungen erheblich verbessert
Ein Fix zur Berücksichtigung interaktiver Mods in Matplotlib. Sieht so aus, als ob dies die Standardoption nur für Ipython ist. Reguläres python3 erfordert diese Anweisung explizit.

Grundsätzlich ist es jetzt möglich, ein Raumschiff von der Erdoberfläche zu "starten" und eine Mission zum Mond zu planen, indem deltaV-Vektorkorrekturen über giveMotion() vorgenommen werden. Als nächstes wird versucht, eine globale Zeitvariable zu implementieren, um eine gleichzeitige Bewegung zu ermöglichen, z. B. der Mond umkreist die Erde, während das Raumschiff ein Gravitationsunterstützungsmanöver ausprobiert.

Kommentare und Verbesserungsvorschläge sind jederzeit willkommen!

Fertig in Python3 mit Matplotlib-Bibliothek

import matplotlib.pyplot as plt
import math
plt.ion()

G = 6.673e-11  # gravity constant
gridArea = [0, 200, 0, 200]  # margins of the coordinate grid
gridScale = 1000000  # 1 unit of grid equals 1000000m or 1000km

plt.clf()  # clear plot area
plt.axis(gridArea)  # create new coordinate grid
plt.grid(b="on")  # place grid

class Object:
    _instances = []
    def __init__(self, name, position, radius, mass):
        self.name = name
        self.position = position
        self.radius = radius  # in grid values
        self.mass = mass
        self.placeObject()
        self.velocity = 0
        Object._instances.append(self)

    def placeObject(self):
        drawObject = plt.Circle(self.position, radius=self.radius, fill=False, color="black")
        plt.gca().add_patch(drawObject)
        plt.show()

    def giveMotion(self, deltaV, motionDirection, time):
        if self.velocity != 0:
            x_comp = math.sin(math.radians(self.motionDirection))*self.velocity
            y_comp = math.cos(math.radians(self.motionDirection))*self.velocity
            x_comp += math.sin(math.radians(motionDirection))*deltaV
            y_comp += math.cos(math.radians(motionDirection))*deltaV
            self.velocity = math.sqrt((x_comp**2)+(y_comp**2))

            if x_comp > 0 and y_comp > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity))  # update motion direction
            elif x_comp > 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 90
            elif x_comp < 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 270

        else:
            self.velocity = self.velocity + deltaV  # in m/s
            self.motionDirection = motionDirection  # degrees
        self.time = time  # in seconds
        self.vectorUpdate()

    def vectorUpdate(self):
        self.placeObject()
        data = []

        for t in range(self.time):
            motionForce = self.mass * self.velocity  # F = m * v
            x_net = 0
            y_net = 0
            for x in [y for y in Object._instances if y is not self]:
                distance = math.sqrt(((self.position[0]-x.position[0])**2) +
                             (self.position[1]-x.position[1])**2)
                gravityForce = G*(self.mass * x.mass)/((distance*gridScale)**2)

                x_pos = self.position[0] - x.position[0]
                y_pos = self.position[1] - x.position[1]

                if x_pos <= 0 and y_pos > 0:  # calculate degrees depending on the coordinate quadrant
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+90

                elif x_pos > 0 and y_pos >= 0:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))+180

                elif x_pos >= 0 and y_pos < 0:
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+270

                else:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))

                x_gF = gravityForce * math.sin(math.radians(gravityDirection))  # x component of vector
                y_gF = gravityForce * math.cos(math.radians(gravityDirection))  # y component of vector

                x_net += x_gF
                y_net += y_gF

            x_mF = motionForce * math.sin(math.radians(self.motionDirection))
            y_mF = motionForce * math.cos(math.radians(self.motionDirection))
            x_net += x_mF
            y_net += y_mF
            netForce = math.sqrt((x_net**2)+(y_net**2))

            if x_net > 0 and y_net > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce))  # update motion direction
            elif x_net > 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 90
            elif x_net < 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 270

            self.velocity = netForce/self.mass  # update velocity
            traveled = self.velocity/gridScale  # grid distance traveled per 1 sec
            self.position = (self.position[0] + math.sin(math.radians(self.motionDirection))*traveled,
                             self.position[1] + math.cos(math.radians(self.motionDirection))*traveled)  # update pos
            data.append([self.position[0], self.position[1]])

            collision = 0
            for x in [y for y in Object._instances if y is not self]:
                if (self.position[0] - x.position[0])**2 + (self.position[1] - x.position[1])**2 <= x.radius**2:
                    collision = 1
                    break
            if collision != 0:
                print("Collision!")
                break

        plt.plot([x[0] for x in data], [x[1] for x in data])

Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(100.0, 100.0), radius=1.737, mass = 7.347e22)  # position not to real scale
Craft = Object(name="SpaceCraft", position=(49.0, 40.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)
Craft.giveMotion(deltaV=2000.0, motionDirection=90, time=60000)
plt.show(block=True)

Wie es funktioniert

Es läuft alles auf zwei Dinge hinaus:

Erstellen von Objekten wie Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)mit Parametern der Position auf dem Gitter (1 Einheit des Gitters ist standardmäßig 1000 km, aber dies kann auch geändert werden), Radius in Gittereinheiten und Masse in kg.
Dem Objekt etwas deltaV zu geben, wie es Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)offensichtlich erforderlich ist, um es Craft = Object(...)zuerst zu erstellen, wie im vorherigen Punkt erwähnt. Die Parameter hier sind deltaVin m/s (beachten Sie, dass die Beschleunigung momentan momentan ist), motionDirectiondie Richtung von deltaV in Grad (stellen Sie sich von der aktuellen Position aus einen 360-Grad-Kreis um das Objekt vor, also ist die Richtung ein Punkt auf diesem Kreis) und schließlich timeist der Parameter, wie viele Sekunden danach wird die deltaV-Push-Trajektorie des Objekts überwacht. Nachfolgender giveMotion()Start von der letzten Position des vorherigen giveMotion().

Fragen:

Ist dies ein gültiger Algorithmus zur Berechnung von Umlaufbahnen?
Was sind die offensichtlichen Verbesserungen, die vorgenommen werden sollten?
Ich habe über die Variable "timeScale" nachgedacht, die die Berechnungen optimiert, da es möglicherweise nicht erforderlich ist, Vektoren und Positionen für jede Sekunde neu zu berechnen. Irgendwelche Gedanken darüber, wie es implementiert werden sollte oder ist es im Allgemeinen eine gute Idee? (Verlust der Genauigkeit vs. verbesserte Leistung)

Grundsätzlich ist es mein Ziel, eine Diskussion über das Thema zu beginnen und zu sehen, wohin es führt. Und, wenn möglich, etwas Neues und Interessantes lernen (oder noch besser - lehren).

Fühlen Sie sich frei zu experimentieren!

Versuchen Sie es mit:

Earth = Object(name="Earth", position=(50.0, 100.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(434.0, 100.0), radius=1.737, mass = 7.347e22)
Craft = Object(name="SpaceCraft", position=(43.0, 100.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=10575.0, motionDirection=180, time=322000)
Craft.giveMotion(deltaV=400.0, motionDirection=180, time=50000)

Mit zwei Zündungen - einer prograden auf der Erdumlaufbahn und einer rückläufigen auf der Mondumlaufbahn - erreichte ich eine stabile Mondumlaufbahn. Sind diese nahe an den theoretisch erwarteten Werten?

Vorgeschlagene Übung: Probieren Sie es in 3 Durchläufen aus – stabile Erdumlaufbahn von der Erdoberfläche, progrades Durchbrennen, um den Mond zu erreichen, rückläufiges Durchbrennen, um die Umlaufbahn um den Mond zu stabilisieren. Versuchen Sie dann, deltaV zu minimieren.

Hinweis: Ich plane, den Code mit ausführlichen Kommentaren für diejenigen zu aktualisieren, die mit der python3-Syntax nicht vertraut sind.

Benutzer2449

Eine sehr gute Idee zum Selbststudium! Wäre es möglich, Ihre Formeln für diejenigen von uns zusammenzufassen, die mit der Python-Syntax nicht vertraut sind?

Zustandsraum

Sicher, denk ich. Ich werde im Code ausführlichere Kommentare für diejenigen machen, die ihn aufgreifen möchten, und die allgemeine Logik in der Frage selbst zusammenfassen.

Benutzer21

Aus dem Kopf heraus: Erwägen Sie die Verwendung eines Vektors für die Geschwindigkeit, anstatt Geschwindigkeit und Richtung unterschiedlich zu behandeln. Wo Sie "F = m * v" sagen, meinen Sie "F = m * a"? Du gehst davon aus, dass sich die Erde nicht bewegt, weil sie so viel schwerer ist als ein Asteroid? Erwägen Sie einen Blick auf github.com/barrycarter/bcapps/blob/master/bc-grav-sim.pl

Zustandsraum

Sie können jedem Objekt, einschließlich der Erde, Bewegung verleihen. Zu Testzwecken habe ich nur Objekt -> Erdbeziehung in die Hauptschleife aufgenommen. Es lässt sich leicht umrechnen, dass sich jedes Objekt auf alle anderen erstellten Objekte bezieht. Und jedes Objekt kann seinen eigenen Bewegungsvektor haben. Grund, warum ich es nicht getan habe - sehr langsame Berechnungen selbst für 1 Objekt. Ich hoffe, dass die Skalierung von Zeiteinheiten sehr hilfreich sein wird, aber ich bin mir immer noch nicht sicher, wie ich es richtig machen soll.

Benutzer21

OK. Ein Gedanke: Führen Sie die Simulation für zwei reale Objekte durch (z. B. Erde/Mond oder Erde/Sonne) und vergleichen Sie Ihre Ergebnisse auf Genauigkeit mit ssd.jpl.nasa.gov/?horizons ? Es wird aufgrund von Störungen durch andere Quellen nicht perfekt sein, aber es wird Ihnen eine Vorstellung von der Genauigkeit geben?

Zustandsraum

Super, danke. Allein dieser Link war es wert, die Diskussion zu beginnen. Also, erster Schritt - ich werde versuchen, den Code so umzuschreiben, dass Objekte alle anderen Objekte auf Schwerkraftvektoren überprüfen. Übrigens, das dringendste Problem der Zeitskalierung ist immer noch hier. Und es wäre interessant zu prüfen, wie sich die Ergebnisse verbessern oder verschlechtern, wenn der Zeitskalenparameter geändert wird.

Zustandsraum

Oh, und F=m*v ist wirklich das. Die Annahme ist, dass das Objekt, wenn es zuerst in Bewegung versetzt wird, zum absoluten Stillstand kommt und dadurch die Geschwindigkeit gleich der Beschleunigung ist. Wie die sofortige Geschwindigkeit von 0 -> beliebiger Wert, und geht mit dieser Geschwindigkeit unbegrenzt weiter, wenn keine anderen Kräfte auf das Objekt einwirken.

Antworten (1)

Übung: 2D-Orbitalmechanik-Simulation (Python)

Eine sehr gute Idee zum Selbststudium! Wäre es möglich, Ihre Formeln für diejenigen von uns zusammenzufassen, die mit der Python-Syntax nicht vertraut sind?
Sicher, denk ich. Ich werde im Code ausführlichere Kommentare für diejenigen machen, die ihn aufgreifen möchten, und die allgemeine Logik in der Frage selbst zusammenfassen.
Aus dem Kopf heraus: Erwägen Sie die Verwendung eines Vektors für die Geschwindigkeit, anstatt Geschwindigkeit und Richtung unterschiedlich zu behandeln. Wo Sie "F = m * v" sagen, meinen Sie "F = m * a"? Du gehst davon aus, dass sich die Erde nicht bewegt, weil sie so viel schwerer ist als ein Asteroid? Erwägen Sie einen Blick auf github.com/barrycarter/bcapps/blob/master/bc-grav-sim.pl
Sie können jedem Objekt, einschließlich der Erde, Bewegung verleihen. Zu Testzwecken habe ich nur Objekt -> Erdbeziehung in die Hauptschleife aufgenommen. Es lässt sich leicht umrechnen, dass sich jedes Objekt auf alle anderen erstellten Objekte bezieht. Und jedes Objekt kann seinen eigenen Bewegungsvektor haben. Grund, warum ich es nicht getan habe - sehr langsame Berechnungen selbst für 1 Objekt. Ich hoffe, dass die Skalierung von Zeiteinheiten sehr hilfreich sein wird, aber ich bin mir immer noch nicht sicher, wie ich es richtig machen soll.
OK. Ein Gedanke: Führen Sie die Simulation für zwei reale Objekte durch (z. B. Erde/Mond oder Erde/Sonne) und vergleichen Sie Ihre Ergebnisse auf Genauigkeit mit ssd.jpl.nasa.gov/?horizons ? Es wird aufgrund von Störungen durch andere Quellen nicht perfekt sein, aber es wird Ihnen eine Vorstellung von der Genauigkeit geben?
Super, danke. Allein dieser Link war es wert, die Diskussion zu beginnen. Also, erster Schritt - ich werde versuchen, den Code so umzuschreiben, dass Objekte alle anderen Objekte auf Schwerkraftvektoren überprüfen. Übrigens, das dringendste Problem der Zeitskalierung ist immer noch hier. Und es wäre interessant zu prüfen, wie sich die Ergebnisse verbessern oder verschlechtern, wenn der Zeitskalenparameter geändert wird.
Oh, und F=m*v ist wirklich das. Die Annahme ist, dass das Objekt, wenn es zuerst in Bewegung versetzt wird, zum absoluten Stillstand kommt und dadurch die Geschwindigkeit gleich der Beschleunigung ist. Wie die sofortige Geschwindigkeit von 0 -> beliebiger Wert, und geht mit dieser Geschwindigkeit unbegrenzt weiter, wenn keine anderen Kräfte auf das Objekt einwirken.

David Ketcheson · Answer 1

Die beiden beteiligten Körper sollen Massen haben $m_1, m_2$ . Beginnen Sie mit dem zweiten Newtonschen Gesetz

F = m a

$F = ma$ wo

a

$a$ ist Beschleunigung. Die Gravitationskraft auf Körper 2 von Körper 1 ist gegeben durch

F_{21} = \frac{G m_{1} m_{2}}{| r_{21} |^{3}} r_{21}

$F_{21} = \frac{G m_1 m_2}{|r_{21}|^3}r_{21}$

wo $r_{21}$ ist der relative Positionsvektor für die beiden fraglichen Körper. Die Kraft auf Körper 1 von Körper 2 ist natürlich $F_{12}=-F_{21}$ seit $r_{12}=-r_{21}$ . Wenn wir die Positionen mit bezeichnen $(x_1,y_1)$ und $(x_2,y_2)$ , dann

r_{21} = (\begin{matrix} x_{1} - x_{2} \\ j_{1} - j_{2} \end{matrix}) .

$r_{21} = \begin{pmatrix} x_1-x_2 \\ y_1-y_2 \end{pmatrix}.$

und

| r | = \sqrt{(x_{1} - x_{2})^{2} + (j_{1} - j_{2})^{2}} .

$|r| = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2}.$ Dann Newtonsches Gesetz

a = F / m

$a=F/m$ wird

\begin{aligned} x_{1}^{″} (t) & = \frac{G m_{2} (x_{2} - x_{1})}{| r |^{3}} \\ j_{1}^{″} (t) & = \frac{G m_{2} (j_{2} - j_{1})}{| r |^{3}} \\ x_{2}^{″} (t) & = \frac{G m_{1} (x_{1} - x_{2})}{| r |^{3}} \\ j_{2}^{″} (t) & = \frac{G m_{1} (j_{1} - j_{2})}{| r |^{3}} . \end{aligned}

$\begin{align} x_1''(t) & = \frac{G m_2 (x_2-x_1)}{|r|^3} \\ y_1''(t) & = \frac{G m_2 (y_2-y_1)}{|r|^3} \\ x_2''(t) & = \frac{G m_1 (x_1-x_2)}{|r|^3} \\ y_2''(t) & = \frac{G m_1 (y_1-y_2)}{|r|^3}. \end{align}$

Zusammen mit den Anfangspositionen und -geschwindigkeiten umfasst dieses System gewöhnlicher Differentialgleichungen (ODEs) ein Anfangswertproblem. Der übliche Ansatz besteht darin, dies als System erster Ordnung aus 8 Gleichungen zu schreiben und eine Runge-Kutta- oder Mehrschrittmethode anzuwenden, um sie zu lösen.

Wenn Sie etwas Einfaches wie Vorwärts-Euler oder Rückwärts-Euler anwenden, sehen Sie, wie sich die Erde ins Unendliche bzw. in Richtung Sonne windet, aber das ist ein Effekt der numerischen Fehler. Wenn Sie eine genauere Methode verwenden, wie die klassische Runge-Kutta-Methode 4. Ordnung, werden Sie feststellen, dass sie eine Weile nahe an einer wahren Umlaufbahn bleibt, sich aber schließlich immer noch nach innen windet (solange der Zeitschritt klein genug ist). oder geht gegen unendlich (wenn der Zeitschritt zu groß ist). Der richtige Ansatz besteht darin, eine symplektische Methode zu verwenden, die die Erde auf der richtigen Umlaufbahn hält – obwohl ihre Phase aufgrund numerischer Fehler immer noch verschoben sein wird.

Für das 2-Körper-Problem ist es möglich, ein einfacheres System abzuleiten, indem Sie Ihr Koordinatensystem um den Massenmittelpunkt herum aufbauen. Aber ich denke, die obige Formulierung ist klarer, wenn dies neu für Sie ist.

Immer noch verdauen. Zu viele unbekannte Wörter für mich, aber irgendwie habe ich das Gefühl, dass ich es irgendwann schaffen werde. Im Moment ist mein eigener Algorithmus gut genug, damit die Dinge einfach funktionieren. Aber wenn ich die simultane Bewegung einstecke, werde ich gezwungen sein, in die Literatur zu gehen und etwas über geeignete Algorithmen zu lesen. Angesichts der Tatsache, dass die Einschränkungen moderner Hardware viel lockerer sind, kann ich es mir leisten, mit einfachen Gleichungen herumzuspielen. Befürchte aber nicht mehr lange.
Tatsächlich sind die symplektischen Methoden bei weitem die genauesten, aber ich denke, dass es für jemanden ohne wissenschaftlichen Hintergrund schwierig ist, sie umzusetzen. Stattdessen können Sie die sehr einfache Euler-Methode zusammen mit der Feynman-Korrektur verwenden. Ich glaube nicht, dass Sie etwas Komplexeres als das für Selbstbildungszwecke brauchen.

Übung: 2D-Orbitalmechanik-Simulation (Python)

Zustandsraum

Benutzer2449

Zustandsraum

Benutzer21

Zustandsraum

Benutzer21

Zustandsraum

Zustandsraum

Antworten (1)

David Ketcheson

Zustandsraum

Zustandsraum

Chris Papageorgakis

Ist es möglich, ein stabiles 3-Körper-System zu haben, das in einem perfekten Kreis kreist?

Wenn eine Sonde am Rand der Einflusssphäre der Erde kreisen würde, wie langsam würde sie umkreisen?

Können Sie das Größenmuster der Hill-Sphären der Objekte des Sonnensystems erklären?

Warum benötigen wir bei der Berechnung der sechs Kepler-Orbitalparameter sowohl die Exzentrizität als auch die große Halbachse? Sagt dir das eine nicht das andere?

Könnte eine Drohne der Marsatmosphäre entkommen?

Was ist der Unterschied zwischen Einflusssphäre und Hügelsphäre?

Was ist die optimale Fluchtbahn aus der Nähe eines Schwarzen Lochs?

Ist die Erdoberfläche „in der Umlaufbahn“?

Wie kann ich die Umlaufzeiten in einem Doppelsternsystem berechnen?

Könnte ein Mond von der Umlaufbahn eines Planeten in einen anderen überführt werden?