Ballistik-Modellierung mit Cx(M) und ρ(h) in Python

Numerische Modellierung der äußeren Ballistik mit variabler Luftdichte

In Fortsetzung der Analyse von Problemen der äußeren Ballistik werden Fälle untersucht, bei denen der Luftwiderstandsbeiwert Cx von der Geschwindigkeit und die Luftdichte von der Höhe abhängen. Der vierte und fünfte Fall erfordern numerische Methoden, da eine analytische Lösung unmöglich ist. Der Flug eines Projektils wird mit den Anfangsparametern simuliert: Geschwindigkeit 600 m/s, Winkel 55°, Masse 50 kg.

Das Luftwiderstandsgesetz von 1943 definiert Cx als Funktion der Mach-Zahl M = v/a. Die Dichte ρ(h) wird durch eine exponentielle Abhängigkeit angenähert: ρ(h) = ρ₀ × exp(-h/H), wobei H = 8000 m die Skalenhöhe ist.

Abhängigkeit von Cx von der Mach-Zahl

Der Koeffizient Cx ist eine stückweise Funktion:

• M ≤ 0,73: Cx = 0,157
• 0,73 < M < 0,82: Cx = 0,033M + 0,133
• 0,82 ≤ M < 0,91: Cx = 0,161 + 3,9(M - 0,823)²
• 0,91 ≤ M < 1,00: Cx = 1,5M - 1,176
• M ≥ 1,00: Cx = 0,384 - 1,6(M - 1,176)²

Die Widerstandskraft F = 0,5 ρ(h) v² S Cx (v/|v|), wobei S die Querschnittsfläche (0,01 m²) ist. Bewegungsgleichungen:

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

g = 9,81 m/s², a = 343 m/s.

Implementierung der numerischen Integration

Zur Lösung des Systems von ODEs wird die Runge-Kutta-Methode 4. Ordnung (RK4) angewendet. Integrationsschritt dt = 0,01 s, Endzeit t_end = 98,4 s.

Vollständiger Simulationscode in Python:

import numpy as np
import matplotlib.pyplot as plt

# Konstanten
g = 9.81  # m/s²
m = 50.0  # Körpermasse, kg
S = 0.01  # Querschnittsfläche, m²
rho0 = 1.225  # Luftdichte auf Meereshöhe, kg/m³
H = 8000  # Skalenhöhe, m
a = 343  # Schallgeschwindigkeit, m/s
v0 = 600  # Anfangsgeschwindigkeit, m/s
theta = np.radians(55)  # Anfangswinkel, rad
dt = 0.01  # Integrationsschritt, s
t_end = 98.4  # Endzeit, s

# Anfangsbedingungen
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)

# Arrays zur Speicherung der Ergebnisse
t = np.arange(0, t_end, dt)
X, Y = np.zeros_like(t), np.zeros_like(t)
VX, VY = np.zeros_like(t), np.zeros_like(t)


# Funktion zur Berechnung von cx(M)
def cx(M):
    if M <= 0.73:
        return 0.157
    elif M < 0.82:
        return 0.033 * M + 0.133
    elif M < 0.91:
        return 0.161 + 3.9 * (M - 0.823) ** 2
    elif M < 1.00:
        return 1.5 * M - 1.176
    else:
        return 0.384 - 1.6 * (M - 1.176) ** 2


# Luftdichtefunktion
def rho(h):
    return rho0 * np.exp(-h / H)


# System von ODEs
def f(state, t):
    x, y, vx, vy = state
    v = np.sqrt(vx ** 2 + vy ** 2)
    M = v / a
    rho_h = rho(y)
    cx_val = cx(M)

    Fx = 0.5 * rho_h * v ** 2 * S * cx_val * vx / v
    Fy = 0.5 * rho_h * v ** 2 * S * cx_val * vy / v

    dxdt = vx
    dydt = vy
    dvxdt = -Fx / m
    dvydt = -g - Fy / m

    return np.array([dxdt, dydt, dvxdt, dvydt])


# Runge-Kutta-Methode 4. Ordnung
def rk4_step(state, t, dt):
    k1 = dt * f(state, t)
    k2 = dt * f(state + 0.5 * k1, t + 0.5 * dt)
    k3 = dt * f(state + 0.5 * k2, t + 0.5 * dt)
    k4 = dt * f(state + k3, t + dt)
    return state + (k1 + 2 * k2 + 2 * k3 + k4) / 6


# Integration
state = np.array([x, y, vx, vy])
for i, t_val in enumerate(t):
    X[i], Y[i], VX[i], VY[i] = state
    state = rk4_step(state, t_val, dt)

# Visualisierung
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Flugbahn')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('Flugbahn des Körpers')
plt.grid(True)
plt.legend()
plt.show()

# Ausgabe der Ergebnisse
print(f"Flugweite: {X[-1]:.2f} m")
print(f"Maximale Höhe: {np.max(Y):.2f} m")

Simulationsergebnisse

Flugweite: 29 km. Maximale Höhe: 11,5 km. Die Flugbahn berücksichtigt den Übergang durch Überschallgeschwindigkeiten und die Luftverdünnung in der Höhe.

Wichtige Annahmen:

• Exponentielles Dichtemodell (vereinfacht; die reale Atmosphäre erfordert genauere Profile).
• Quadratischer Luftwiderstand ohne Berücksichtigung anderer Effekte (Drall, Wind).

| Parameter | Wert |

|-----------|-------|

| Flugweite | 29.000 m |

| Max. Höhe | 11.500 m |

| Flugzeit | 98,4 s |

| Anfangsgeschwindigkeit | 600 m/s |

Wichtige Erkenntnisse

• Numerische Methoden sind unerlässlich für realistische Modelle mit variablem Cx(M) und ρ(h).
• RK4 gewährleistet Genauigkeit mit dt=0,01 s für ballistische Aufgaben.
• Das Gesetz von 1943 modelliert den transsonischen Bereich (M≈0,8–1,0) genau.
• Exponentielle ρ(h) bietet eine sinnvolle Näherung bis 15 km.

— Editorial Team