albertus
Registriertes Mitglied
20.3 Der finale 3D-Solver in Python (Inklusive statistischer Streuung)
Das folgende Skript führt eine Schleife von 50 autonomen Wiedereintrittsflügen durch. Jedes Mal werden die aerodynamischen Parameter leicht variiert, um die reale Streuung am Boden sichtbar zu machen.
Python:
import numpy as np
import matplotlib.pyplot as plt
# Physikalische Basiskonstanten
G = 6.67430e-11
M_EARTH = 5.972e24
R_EARTH = 6371000
RHO_0 = 1.225
H_SCALE = 7500
OMEGA_EARTH = 7.292115e-5
LATITUDE = np.radians(46.0) # Zielgebiet Kasachstan
def run_single_entry(rho_noise, cd_noise):
""" Simuliert einen einzelnen 3D-Wiedereintritt mit Rauschen """
h = 100000.0 # Start bei 100 km Hoehe
# 3D-Geschwindigkeitsvektoren (m/s)
v_z = 200.0 # Initiale Sinkgeschwindigkeit
v_forward = 7500.0 # Orbitale Vorwaertsgeschwindigkeit in der Flugbahn
v_lateral = 0.0 # Seitliche Drift (Start bei 0)
# Positionen (m)
x_forward = 0.0
y_lateral = 0.0
m = 2400.0 # Kapselmasse
dt = 0.1 # Zeitschritt
while h > 0 and v_z >= 0:
# Dynamische Atmosphaerendichte mit stochastischem Rauschen
rho = RHO_0 * np.exp(-h / H_SCALE) * rho_noise
g_h = (G * M_EARTH) / (R_EARTH + h)**2
# Strukturumschaltung der Flaechen (Kapsel -> Bremsschirm -> Hauptschirm)
if h > 9500:
A_eff, cd_base = 4.9, 1.3
elif 7500 < h <= 9500:
A_eff, cd_base = 24.0, 1.4
else:
A_eff, cd_base = 1000.0, 1.3
# Widerstandsbeiwert mit stochastischem Rauschen beaufschlagen
cd_eff = cd_base * cd_noise
# 1. Vertikale Achse (Z)
drag_z = 0.5 * rho * v_z**2 * cd_eff * A_eff
a_z = g_h - (drag_z / m)
v_z += a_z * dt
h -= v_z * dt
# 2. Orbitale Vorwaertsachse (X)
drag_forward = 0.5 * rho * v_forward**2 * cd_eff * A_eff
a_forward = - (drag_forward / m)
v_forward += a_forward * dt
x_forward += v_forward * dt
# 3. Laterale Coriolis-Achse (Y)
a_coriolis = 2.0 * OMEGA_EARTH * v_z * np.cos(LATITUDE)
drag_lateral = 0.5 * rho * v_lateral**2 * cd_eff * A_eff * np.sign(v_lateral)
a_lateral = a_coriolis - (drag_lateral / m)
v_lateral += a_lateral * dt
y_lateral += v_lateral * dt
return x_forward / 1000.0, y_lateral / 1000.0 # Rueckgabe in km
# --- MONTE-CARLO-SCHLEIFE (50 Durchlaeufe) ---
np.random.seed(42) # Fuer reproduzierbare Ergebnisse
num_simulations = 50
landings_x = []
landings_y = []
for _ in range(num_simulations):
# Zufaellige Variationen der Natur (Normalverteilung)
rho_variation = np.random.normal(1.0, 0.03) # 3% Standardabweichung Luftdichte
cd_variation = np.random.normal(1.0, 0.015) # 1.5% Standardabweichung Strömungswiderstand
x_land, y_land = run_single_entry(rho_variation, cd_variation)
landings_x.append(x_land)
landings_y.append(y_land)
# Plot der 3D-Landeellipse
plt.figure(figsize=(10, 6))
plt.scatter(landings_x, landings_y, color='crimson', marker='o', alpha=0.7, label='Simulierte Landepunkte')
# Berechne Mittelpunkt und Dimensionen fuer die visuelle Ellipse
mean_x, mean_y = np.mean(landings_x), np.mean(landings_y)
plt.axhline(mean_y, color='gray', linestyle='--', alpha=0.5)
plt.axvline(mean_x, color='gray', linestyle='--', alpha=0.5)
plt.title('Monte-Carlo-Simulation: Reale 3D-Landeellipse der Sojus-Kapsel')
plt.xlabel('Flugstrecke in der Orbitalebene (km)')
plt.ylabel('Laterale Abweichung durch Coriolis-Drift (km)')
plt.grid(True)
plt.legend()
plt.show()
20.4 Mathematisch-physikalische Auswertung des Plots (Wo steckt das 3D?)
Auf den ersten Blick sieht das Diagramm wie eine rein zweidimensionale, flache Punktwolke aus. Als Ingenieur stellt man sich sofort die Frage: Wo sehe ich hier die dreidimensionale Dynamik?Die Antwort ist ein klassischer Fall von Projektion: Der Plot zeigt uns die Draufsicht (Vogelperspektive) der kasachischen Steppe am Aufschlagpunkt (h = 0). Die dritte Dimension – die vertikale Achse (Z) – ist in der Grafik also implizit dadurch enthalten, dass jeder rote Punkt das Ende einer individuellen, dreidimensionalen Trajektorie markiert, die sich aus 100 km Höhe schräg nach unten durch den Raum gefressen hat.
Die 3D-Dynamik offenbart sich unmissverständlich in der strikten Korrelation der Punkte, die eine saubere Diagonale von links oben nach rechts unten bilden:
- Ein Punkt ganz links oben (z. B. bei 746 km Flugstrecke): Hier war die stochastisch simulierte Atmosphäre dichter als normal. Die Kapsel wurde in der orbitalen Flugbahn (X) stark abgebremst und fällt früher zu Boden. Weil sie aber in der Vertikalen (Z) so massiv eingebremst wurde und dadurch viel länger am Fallschirm nach unten schwebte, hatte die Erddrehung deutlich mehr Zeit, die Kapsel über die Corioliskraft seitlich (Y) abzulenken. Daher der hohe Wert auf der Y-Achse.
- Ein Punkt ganz rechts unten (z. B. bei 752 km Flugstrecke): Hier war die Atmosphäre dünner. Die Kapsel erlebte weniger Widerstand, schoss in der Orbitalebene weiter nach vorne, durchmaß die Vertikale schneller und gab der Coriolis-Drift kaum Zeit, sie seitlich zu versetzen.