Modelowanie balistyki z Cx(M) i ρ(h) w Pythonie

Numeryczne modelowanie balistyki zewnętrznej ze zmienną gęstością powietrza

W dalszej analizie zadań balistyki zewnętrznej rozpatruje się przypadki zależności współczynnika oporu czołowego Cx od prędkości oraz gęstości powietrza od wysokości. Czwarty i piąty przypadek wymagają metod numerycznych, ponieważ rozwiązanie analityczne jest niemożliwe. Symulowany jest lot pocisku z parametrami początkowymi: prędkość 600 m/s, kąt 55°, masa 50 kg.

Prawo oporu powietrza z 1943 roku określa Cx jako funkcję liczby Macha M = v/a. Gęstość ρ(h) aproksymowana jest zależnością wykładniczą ρ(h) = ρ₀ × exp(-h/H), gdzie H = 8000 m — wysokość skali.

Zależność Cx od liczby Macha

Współczynnik Cx to funkcja odcinkami:

• 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)²

Siła oporu F = 0.5 ρ(h) v² S Cx (v/|v|), gdzie S — pole przekroju poprzecznego (0.01 m²). Równania ruchu:

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

g = 9.81 m/s², a = 343 m/s.

Implementacja całkowania numerycznego

Do rozwiązania układu równań różniczkowych zwyczajnych stosuje się metodę Rungego-Kutty 4. rzędu (RK4). Krok całkowania dt = 0.01 s, czas końcowy t_end = 98.4 s.

Pełny kod symulacji w Pythonie:

import numpy as np
import matplotlib.pyplot as plt

# Stałe
g = 9.81  # m/s²
m = 50.0  # masa ciała, kg
S = 0.01  # pole przekroju poprzecznego, m²
rho0 = 1.225  # gęstość powietrza na poziomie morza, kg/m³
H = 8000  # wysokość skali, m
a = 343  # prędkość dźwięku, m/s
v0 = 600  # prędkość początkowa, m/s
theta = np.radians(55)  # kąt początkowy, rad
dt = 0.01  # krok całkowania, s
t_end = 98.4  # czas końcowy, s

# Warunki początkowe
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)

# Tablice do przechowywania wyników
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)


# Funkcja do obliczania 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


# Funkcja gęstości powietrza
def rho(h):
    return rho0 * np.exp(-h / H)


# Układ równań różniczkowych zwyczajnych
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])


# Metoda Rungego-Kutty 4. rzędu
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


# Całkowanie
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)

# Wizualizacja
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Trajektoria')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('Trajektoria lotu ciała')
plt.grid(True)
plt.legend()
plt.show()

# Wyświetlenie wyników
print(f"Zasięg lotu: {X[-1]:.2f} m")
print(f"Maksymalna wysokość: {np.max(Y):.2f} m")

Wyniki symulacji

Zasięg lotu: 29 km. Maksymalna wysokość: 11.5 km. Trajektoria uwzględnia przejście przez prędkości naddźwiękowe i rozrzedzenie powietrza na wysokościach.

Kluczowe założenia:

• Model wykładniczy gęstości (uproszczony, rzeczywista atmosfera wymaga dokładniejszych profili).
• Opór kwadratowy bez uwzględnienia innych efektów (wirowanie, wiatr).

| Parametr | Wartość |

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

| Zasięg | 29 000 m |

| Maks. wysokość | 11 500 m |

| Czas lotu | 98.4 s |

| Prędkość pocz. | 600 m/s |

Co jest ważne

• Metody numeryczne są niezbędne dla realistycznych modeli ze zmiennymi Cx(M) i ρ(h).
• RK4 zapewnia dokładność przy dt=0.01 s dla zadań balistycznych.
• Prawo z 1943 r. dokładnie modeluje reżim transsoniczny (M≈0.8–1.0).
• Wykładnicza ρ(h) daje rozsądną aproksymację do 15 km.

— Editorial Team