Exercice: simulation de mécanique orbitale 2D (python)

Juste un petit avertissement au préalable: je n'ai jamais étudié l'astronomie ni aucune science exacte d'ailleurs (pas même l'informatique), alors j'essaie de combler cette lacune par l'auto-éducation. L'astronomie est l'un des domaines qui a retenu mon attention et mon idée de l'auto-éducation est à la tête de l'approche appliquée. Donc, droit au but - c'est un modèle de simulation orbital sur lequel je travaille avec désinvolture lorsque j'ai le temps / l'humeur. Mon objectif principal est de créer un système solaire complet en mouvement et de pouvoir planifier des lancements de vaisseaux spatiaux sur d'autres planètes.

Vous êtes tous libres de reprendre ce projet à tout moment et de vous amuser à expérimenter!

mise à jour!!! (10 novembre)

• la vitesse est maintenant deltaV appropriée et le mouvement supplémentaire calcule maintenant le vecteur somme de la vitesse
• vous pouvez placer autant d'objets statiques que vous le souhaitez, à chaque fois que l'objet unité en mouvement vérifie les vecteurs de gravité de toutes les sources (et vérifie les collisions)
• considérablement amélioré les performances des calculs
• un correctif pour tenir compte du mod interactif dans matplotlib. Il semble que ce soit l'option par défaut uniquement pour ipython. Python3 normal requiert explicitement cette instruction.

Fondamentalement, il est désormais possible de "lancer" un vaisseau spatial à partir de la surface de la Terre et de tracer une mission vers la Lune en effectuant des corrections de vecteur deltaV via giveMotion (). Next in line essaie de mettre en œuvre une variable de temps globale pour permettre un mouvement simultané, par exemple la Lune en orbite autour de la Terre pendant que le vaisseau spatial essaie une manœuvre d'assistance à la gravité.

Les commentaires et suggestions d'amélioration sont toujours les bienvenus!

Fait en Python3 avec la bibliothèque matplotlib

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)

Comment ça fonctionne

Tout se résume à deux choses:

1. Création d'objet comme Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)avec des paramètres de position sur la grille (1 unité de grille est 1000 km par défaut mais cela peut aussi être modifié), rayon en unités de grille et masse en kg.
2. Donner à l'objet un certain deltaV tel Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)qu'il doit évidemment Craft = Object(...)être créé en premier lieu comme mentionné au point précédent. Les paramètres ici sont deltaVen m / s (notez que pour l'instant l'accélération est instantanée), motionDirectionest la direction du deltaV en degrés (à partir de la position actuelle, imaginez un cercle de 360 ​​degrés autour de l'objet, donc la direction est un point sur ce cercle) et enfin le paramètre timeest le nombre de secondes après la trajectoire de poussée deltaV de l'objet sera surveillée. Après giveMotion()arrêt de début de la dernière position précédente giveMotion().

Des questions:

1. Est-ce un algorithme valide pour calculer les orbites?
2. Quelles sont les améliorations évidentes à apporter?
3. J'ai envisagé la variable "timeScale" qui optimisera les calculs, car il pourrait ne pas être nécessaire de recalculer les vecteurs et les positions pour chaque seconde. Avez-vous des réflexions sur la façon de le mettre en œuvre ou est-ce généralement une bonne idée? (perte de précision vs performances améliorées)

Fondamentalement, mon objectif est de lancer une discussion sur le sujet et de voir où il mène. Et, si possible, apprenez (ou mieux - enseignez) quelque chose de nouveau et d'intéressant.

N'hésitez pas à expérimenter!

Essayez d'utiliser:

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)

Avec deux brûlures - une prograde en orbite terrestre et une rétrograde en orbite lunaire, j'ai atteint une orbite lunaire stable. Ces valeurs sont-elles proches des valeurs théoriquement attendues?

Exercice suggéré: essayez-le en 3 brûlures - orbite terrestre stable depuis la surface de la terre, brûlure prograde pour atteindre la lune, brûlure rétrograde pour stabiliser l'orbite autour de la lune. Essayez ensuite de minimiser deltaV.

Remarque: je prévois de mettre à jour le code avec des commentaires détaillés pour ceux qui ne connaissent pas la syntaxe python3.

$m_1, m_2$

$r_{21}$$F_{12}=-F_{21}$$r_{12}=-r_{21}$$(x_1,y_1)$$(x_2,y_2)$

et

Avec les positions et vitesses initiales, ce système d'équations différentielles ordinaires (ODE) comprend un problème de valeur initiale. L'approche habituelle consiste à écrire ceci comme un système de premier ordre de 8 équations et à appliquer une méthode de Runge-Kutta ou à plusieurs étapes pour les résoudre.

Si vous appliquez quelque chose de simple comme Euler avant ou Euler arrière, vous verrez la Terre en spirale vers l'infini ou vers le soleil, respectivement, mais c'est un effet des erreurs numériques. Si vous utilisez une méthode plus précise, comme la méthode Runge-Kutta classique du 4e ordre, vous constaterez qu'elle reste proche d'une véritable orbite pendant un certain temps mais finit par s'éteindre à l'infini. La bonne approche consiste à utiliser une méthode symplectique, qui maintiendra la Terre sur l'orbite correcte - bien que sa phase soit toujours désactivée en raison d'erreurs numériques.

Pour le problème à 2 corps, il est possible de dériver un système plus simple en basant votre système de coordonnées autour du centre de masse. Mais je pense que la formulation ci-dessus est plus claire si c'est nouveau pour vous.