Ejercicio: simulación de mecánica orbital 2D (python)

12

Solo un pequeño descargo de responsabilidad de antemano: nunca he estudiado astronomía o ciencias exactas para ese asunto (ni siquiera TI), así que estoy tratando de llenar este vacío mediante la autoeducación. La astronomía es una de las áreas que ha captado mi atención y mi idea de la autoeducación es de enfoque aplicado. Entonces, directamente al grano: este es un modelo de simulación orbital en el que estoy trabajando casualmente cuando tengo tiempo / estado de ánimo. Mi objetivo principal es crear un sistema solar completo en movimiento y la capacidad de planificar el lanzamiento de naves espaciales a otros planetas.

¡Todos ustedes son libres de elegir este proyecto en cualquier momento y divertirse experimentando!

¡¡¡actualizar!!! (10 de noviembre)

  • la velocidad ahora es deltaV adecuada y al dar un movimiento adicional ahora se calcula el vector suma de velocidad
  • puede colocar tantos objetos estáticos como desee, en cada objeto de unidad en movimiento, verifique los vectores de gravedad de todas las fuentes (y verifique la colisión)
  • mejoró enormemente el rendimiento de los cálculos
  • Una solución para tener en cuenta el mod interactivo en matplotlib. Parece que esta es la opción predeterminada solo para ipython. Python3 regular requiere esa declaración explícitamente.

Básicamente, ahora es posible "lanzar" una nave espacial desde la superficie de la Tierra y planear una misión a la Luna haciendo correcciones vectoriales deltaV a través de giveMotion (). El siguiente en la línea está tratando de implementar una variable de tiempo global para permitir el movimiento simultáneo, por ejemplo, la Luna orbita la Tierra mientras la nave espacial prueba una maniobra de asistencia por gravedad.

¡Los comentarios y sugerencias de mejoras son siempre bienvenidos!

Hecho en Python3 con la biblioteca 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)

Cómo funciona

Todo se reduce a dos cosas:

  1. Crear objetos como Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)con parámetros de posición en la cuadrícula (1 unidad de cuadrícula es 1000 km por defecto, pero esto también se puede cambiar), radio en unidades de cuadrícula y masa en kg.
  2. Darle al objeto algo deltaV, como Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)obviamente requiere Craft = Object(...)que se cree en primer lugar como se menciona en el punto anterior. Los parámetros aquí están deltaVen m / s (tenga en cuenta que por ahora la aceleración es instantánea), motionDirectiones la dirección del deltaV en grados (desde la posición actual imagine un círculo de 360 ​​grados alrededor del objeto, entonces la dirección es un punto en ese círculo) y finalmente el parámetro timees cuántos segundos después de la trayectoria de empuje deltaV del objeto será monitoreado. giveMotion()Comienzo posterior desde la última posición de la anterior giveMotion().

Preguntas:

  1. ¿Es este un algoritmo válido para calcular órbitas?
  2. ¿Cuáles son las mejoras obvias que se harán?
  3. He estado considerando la variable "timeScale" que optimizará los cálculos, ya que podría no ser necesario volver a calcular los vectores y las posiciones por cada segundo. ¿Alguna idea sobre cómo debería implementarse o es generalmente una buena idea? (pérdida de precisión versus rendimiento mejorado)

Básicamente, mi objetivo es comenzar una discusión sobre el tema y ver a dónde conduce. Y, si es posible, aprenda (o incluso mejor, enseñe) algo nuevo e interesante.

¡Siéntete libre de experimentar!

Intenta usar:

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)

Con dos quemaduras: un programa en órbita terrestre y otro retrógrado en órbita lunar, logré una órbita lunar estable. ¿Están estos cerca de los valores teóricamente esperados?

Ejercicio sugerido: Pruébelo en 3 quemaduras: órbita estable de la Tierra desde la superficie de la Tierra, programe quemaduras para llegar a la Luna, quemadura retrógrada para estabilizar la órbita alrededor de la Luna. Luego intente minimizar deltaV.

Nota: planeo actualizar el código con comentarios extensos para aquellos que no están familiarizados con la sintaxis de python3.

orbit gravity n-body-simulations espacio de Estados
fuente
¡Una muy buena idea para la autoeducación! ¿Sería posible resumir sus fórmulas para aquellos de nosotros que no estamos familiarizados con la sintaxis de Python?
Seguro, creo. Haré comentarios más extensos en el código para aquellos que quieran recogerlo y resumir la lógica general en la pregunta misma.
Statespace
Fuera de mi cabeza: considere usar un vector para la velocidad en lugar de tratar la velocidad y la dirección de manera diferente. Cuando dices "F = m * v", ¿quieres decir "F = m * a"? ¿Estás asumiendo que la Tierra no se mueve porque es mucho más pesada que el asteroide? Considere mirar github.com/barrycarter/bcapps/blob/master/bc-grav-sim.pl
barrycarter
Puedes mover cualquier objeto, incluida la Tierra. Para fines de prueba, incluí solo objeto -> Relación de tierra en el bucle principal. Se puede convertir fácilmente que cada objeto se relaciona con todos los demás objetos que se crean. Y cada objeto puede tener su propio vector de movimiento. Razón por la que no lo hice: cálculos muy lentos incluso para 1 objeto. Espero que las unidades de tiempo de escala ayuden mucho, pero todavía no estoy seguro de cómo hacerlo correctamente.
Statespace
1
OKAY. Un pensamiento: ¿hace la simulación para dos objetos reales (por ejemplo, Tierra / Luna o Tierra / Sol) y compara sus resultados con ssd.jpl.nasa.gov/?horizons para mayor precisión? No será perfecto debido a las perturbaciones de otras fuentes, pero ¿le dará alguna idea de la precisión?
barrycarter

Respuestas:

11

m1,m2

F=ma
a

F21=Gm1m2|r21|3r21

r21F12=−F21r12=−r21(x1,y1)(x2,y2)

r21=(x1−x2y1−y2).

y

|r|=(x1−x2)2+(y1−y2)2.
a=F/m

x1″(t)=Gm2(x2−x1)|r|3y1″(t)=Gm2(y2−y1)|r|3x2″(t)=Gm1(x1−x2)|r|3y2″(t)=Gm1(y1−y2)|r|3.

Junto con las posiciones y velocidades iniciales, este sistema de ecuaciones diferenciales ordinarias (EDO) comprende un problema de valor inicial. El enfoque habitual es escribir esto como un sistema de primer orden de 8 ecuaciones y aplicar un método Runge-Kutta o multipaso para resolverlas.

Si aplica algo simple como Euler hacia adelante o hacia atrás, verá que la Tierra gira en espiral hacia el infinito o hacia el sol, respectivamente, pero eso es un efecto de los errores numéricos. Si usa un método más preciso, como el método clásico Runge-Kutta de 4to orden, encontrará que permanece cerca de una órbita verdadera por un tiempo, pero aún así se va al infinito. El enfoque correcto es utilizar un método simpléctico, que mantendrá a la Tierra en la órbita correcta, aunque su fase seguirá estando apagada debido a errores numéricos.

Para el problema de los 2 cuerpos, es posible derivar un sistema más simple basando su sistema de coordenadas alrededor del centro de masa. Pero creo que la formulación anterior es más clara si esto es nuevo para usted.

David Ketcheson
fuente
Esto tomará algún tiempo para digerir.
Statespace
Aún digiriendo. Demasiadas palabras desconocidas para mí, pero de alguna manera siento que en algún momento llegaré allí. Por ahora, mi propio algoritmo es lo suficientemente bueno para que las cosas simplemente funcionen. Pero cuando conecte el movimiento simultáneo, me veré obligado a acceder a la literatura y leer sobre algoritmos adecuados. Dado que las limitaciones del hardware moderno son mucho más flexibles, puedo permitirme perder el tiempo con ecuaciones simples. Sin embargo, no tengo miedo por mucho tiempo.
Statespace
De hecho, los métodos simplécticos son, con mucho, los más precisos, pero creo que es difícil para alguien sin experiencia en la ciencia implementarlos. En su lugar, puede utilizar el método de Euler muy simple junto con la corrección de Feynman. No creo que necesite algo más complejo que eso para fines de autoeducación.
chrispap