Aufgaben
Lösung zu Aufgabe 8.1
Durch wahlloses Herumprobieren ist es gar nicht so einfach eine Bahnkurve zu skizzieren, die dem Swing-By-Manöver aus der Aufgabe gleicht. Es ist daher hilfreich, sich eine geeignete Strategie zu überlegen:
Wir lassen den Jupiter mit der Geschwindigkeit von $13\,\textrm{km/s}$ aus einem bestimmten Abstand $d_\textrm{Jupiter}$ einfliegen. Die Sonde starten wir in einer gegebenen Entfernung vom Koordinatenursprung und lassen diese mit einer Geschwindigkeit von $9\,\textrm{km/s}$ exakt auf den Koordinatenursprung zufliegen, wobei sie unter einem Winkel $\alpha$ zur negativen $x$-Achse startet. Die Masse der Sonde setzen wir mit $1000\,\textrm{kg}$ an. Da der Jupiter viele Größenordnungen mehr Masse hat, spielt der genaue Zahlenwert keine Rolle.
Als Anfangsabstand der Sonde vom Koordinatenursprung wählen wir einen großen Abstand von $1{,}5\cdot10^{10}\,\textrm{m}$ und als Anflugwinkel $\alpha=60^\circ$. Der einzige Parameter der jetzt noch bestimmt, um welchen Winkel die Sonde aus ihrer ursprünglichen Flugbahn abgelenkt wird, ist die Anfangsposition des Jupiters, die durch den Abstand $d_\textrm{Jupiter}$ gegeben ist. Diese Zahl wählen wir durch Probieren nun so, dass die Flugbahn schön symmetrisch bezüglich der $x$-Achse wird.
Das unten abgedruckte Programm gibt eine Animation der Bahnen der Sonde und des Jupiters aus.
Außerdem wird der Zeitverlauf des Abstandes der beiden Körper und der Betrag der Geschwindigkeit der Sonde dargestellt.
Sie erkennen, dass die Geschwindigkeit der Sonde sich während des Manövers von 9 km/s auf ca. 22 km/s erhöht.
Die maximal mögliche Geschwindigkeit, die die Raumsonde bei einem solchen Manöver erreichen könnte, würde sich ergeben, wenn die Raumsonde exakt in positiver $x$-Richtung anfliegen würde und am Ende wieder exakt in negativer $x$-Richtung wegfliegen würde. Das kann man sich anschaulich so vorstellen, als würde die Sonde wie ein Gummiball vom Jupiter abprallen. In dieser Modellvorstellung kann man die maximale Geschwindigkeit der Sonde leicht ausrechnen: Im Ruhesystem des Jupiters bewegt sich die Sonde mit einer Geschwindigkeit von $22\,\textrm{km/s}$. Da der Jupiter viel schwerer ist als die Sonde, bewegt sich die Sonde nach dem Stoß mit einer Geschwindigkeit von $22\,\textrm{km/s}$ vom Jupiter weg. In dem Bezugssystem, in dem sich der Jupiter mit einer Geschwindigkeit von $13\,\textrm{km/s}$ in negative $x$-Richtung bewegt, ist das also eine Geschwindigkeit von $35\,\textrm{km/s}$.
Sie können bei dem Swingby-Manöver die Endgeschwindigkeit der Sonde erhöhen, indem Sie die Sonde mit einem flacheren Winkel einfliegen lassen und die Anfangsposition des Jupiters entsprechend anpassen. Sie werden feststellen, dass man sich für kleine Anflugwinkel der Endgeschwindigkeit von $35\,\textrm{km/s}$ nähert. Allerdings wird dabei der Minimalabstand zwischen Sonde und Jupiter schnell kleiner, sodass es bei zu kleinem Anflugwinkel zu einer Kollision kommen würde.
Mehrteilchen/Loesungen/swingby.py
"""Simulation eines Swing-by-Manövers."""
import math
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation
import scipy.integrate
# Konstanten: 1 Stunde, 1 Tag und 1 Jahr [s].
stunde = 60 * 60
tag = 24 * stunde
jahr = 365.25 * tag
# Simulationszeit und Zeitschrittweite [s].
t_max = 36 * tag
dt = 1 * stunde
# Masse des Planeten und der Raumsonde [kg].
m_planet = 1.898e27
m_sonde = 1e3
# Anflugwinkel der Raumsonde relativ zur Bewegungsrichtung
# des Planeten.
alpha = math.radians(60)
# Anfangsentfernung der Körper vom Koordinatenursprung.
abstand_sonde = 15e9
abstand_planet = 20.18e9
# Radius des Planeten.
radius_planet = 6.9911e7
# Anfangsgeschwindigkeiten der Körper.
betrag_v0_sonde = 9e3
betrag_v0_planet = 13e3
# Anfangspositionen der Körper [m].
r0_planet = np.array([abstand_planet, 0.0])
r0_sonde = abstand_sonde * np.array([-np.cos(alpha),
-np.sin(alpha)])
# Gravitationskonstante [m³ / (kg * s²)].
G = 6.6743e-11
# Anfangsgeschwindigkeiten der Körper [m/s].
v0_planet = betrag_v0_planet * np.array([-1.0, 0.0])
v0_sonde = betrag_v0_sonde * np.array([np.cos(alpha),
np.sin(alpha)])
def dgl(t, u):
"""Berechne die rechte Seite der Differentialgleichung."""
r1, r2, v1, v2 = np.split(u, 4)
a1 = G * m_sonde / np.linalg.norm(r2 - r1) ** 3 * (r2 - r1)
a2 = G * m_planet / np.linalg.norm(r1 - r2) ** 3 * (r1 - r2)
return np.concatenate([v1, v2, a1, a2])
# Lege den Zustandsvektor zum Zeitpunkt t=0 fest.
u0 = np.concatenate((r0_planet, r0_sonde, v0_planet, v0_sonde))
# Löse die Bewegungsgleichung numerisch.
result = scipy.integrate.solve_ivp(dgl, [0, t_max], u0, rtol=1e-9,
t_eval=np.arange(0, t_max, dt))
t = result.t
r_sonde, r_planet, v_sonde, v_planet = np.split(result.y, 4)
# Berechne den Abstand der Raumsonde vom Planeten.
abstand = np.linalg.norm(r_sonde - r_planet, axis=0)
# Gibt die minimale Entfernung der Raumsonde vom Planeten in
# Vielfachen des Radius des Planeten an.
abstand_min = np.min(abstand)
print(f"Min. Abstand: {abstand_min / radius_planet:.1f} "
"Planetenradien")
# Berechne den Geschwindigkeitsbetrag der Raumsonde.
betrag_v_sonde = np.linalg.norm(v_planet, axis=0)
# Erzeuge eine Figure und eine Axes für die Animation.
fig_anim = plt.figure()
fig_anim.set_tight_layout(True)
ax_anim = fig_anim.add_subplot(1, 1, 1)
ax_anim.set_xlabel('$x$ [m]')
ax_anim.set_ylabel('$y$ [m]')
ax_anim.grid()
# Plotte die Bahnkurve der beiden Körper.
ax_anim.plot(r_sonde[0], r_sonde[1], '-r')
ax_anim.plot(r_planet[0], r_planet[1], '-b')
# Erzeuge eine zweite Figure für die Plots.
fig_plot = plt.figure()
fig_plot.set_tight_layout(True)
# Erzeuge eine Axes für den Abstand als Funktion der Zeit.
ax_abstand = fig_plot.add_subplot(1, 2, 1)
ax_abstand.set_xlabel("$t$ [Tage]")
ax_abstand.set_ylabel("Abstand [m]")
ax_abstand.grid()
ax_abstand.plot(t / tag, abstand)
# Erzeuge eine Axes für die Geschwindigkeit der Raumsonde.
ax_geschwindigkeit = fig_plot.add_subplot(1, 2, 2)
ax_geschwindigkeit.set_xlabel("$t$ [Tage]")
ax_geschwindigkeit.set_ylabel("Geschwindigkeit [km/s]")
ax_geschwindigkeit.grid()
ax_geschwindigkeit.plot(t / tag, betrag_v_sonde / 1e3)
# Erzeuge Punktplots für die Positionen der Himmelskörper.
plot_sonde, = ax_anim.plot([], [], 'o', color='red')
plot_planet, = ax_anim.plot([], [], 'o', color='blue')
def update(n):
"""Aktualisiere die Grafik zum n-ten Zeitschritt."""
plot_sonde.set_data(r_sonde[:, n].reshape(-1, 1))
plot_planet.set_data(r_planet[:, n].reshape(-1, 1))
return plot_sonde, plot_planet
# Erzeuge das Animationsobjekt und starte die Animation.
ani = mpl.animation.FuncAnimation(fig_anim, update, frames=t.size,
interval=40, blit=True)
plt.show()