Numerické modelování vnější balistiky s proměnnou hustotou vzduchu
V pokračování analýzy úloh vnější balistiky se zkoumají případy se závislostí koeficientu čelního odporu Cx na rychlosti a hustoty vzduchu na výšce. Čtvrtý a pátý případ vyžadují numerické metody, protože analytické řešení není možné. Modeluje se let střely s počátečními parametry: rychlost 600 m/s, úhel 55°, hmotnost 50 kg.
Zákon odporu vzduchu z roku 1943 definuje Cx jako funkci Machova čísla M = v/a. Hustota ρ(h) je aproximována exponenciální závislostí ρ(h) = ρ₀ × exp(-h/H), kde H = 8000 m je měřítková výška.
Závislost Cx na Machově čísle
Koeficient Cx je po částech definovaná funkce:
- 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)²
Síla odporu F = 0.5 ρ(h) v² S Cx (v/|v|), kde S je plocha průřezu (0.01 m²). Pohybové rovnice:
ẋ = vx
ẏ = vy
v̇x = -Fx/m
v̇y = -g - Fy/m
g = 9.81 m/s², a = 343 m/s.
Implementace numerické integrace
K řešení soustavy obyčejných diferenciálních rovnic se používá metoda Runge-Kutta 4. řádu (RK4). Krok integrace dt = 0.01 s, konečný čas t_end = 98.4 s.
Úplný kód simulace v Pythonu:
import numpy as np
import matplotlib.pyplot as plt
# Konstanty
g = 9.81 # m/s²
m = 50.0 # hmotnost tělesa, kg
S = 0.01 # plocha průřezu, m²
rho0 = 1.225 # hustota vzduchu na úrovni moře, kg/m³
H = 8000 # měřítková výška, m
a = 343 # rychlost zvuku, m/s
v0 = 600 # počáteční rychlost, m/s
theta = np.radians(55) # počáteční úhel, rad
dt = 0.01 # krok integrace, s
t_end = 98.4 # konečný čas, s
# Počáteční podmínky
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)
# Pole pro ukládání výsledků
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)
# Funkce pro výpočet 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
# Funkce hustoty vzduchu
def rho(h):
return rho0 * np.exp(-h / H)
# Soustava ODR
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 Runge-Kutta 4. řá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
# Integrace
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)
# Vizualizace
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Trajektorie')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('Trajektorie letu tělesa')
plt.grid(True)
plt.legend()
plt.show()
# Výsledky
print(f"Dosah letu: {X[-1]:.2f} m")
print(f"Maximální výška: {np.max(Y):.2f} m")Výsledky simulace
Dosah letu: 29 km. Maximální výška: 11.5 km. Trajektorie zohledňuje přechod přes nadzvukové rychlosti a zředění vzduchu ve výškách.
Klíčová zjednodušení:
- Exponenciální model hustoty (zjednodušený, reálná atmosféra vyžaduje přesnější profily).
- Kvadratický odpor bez zohlednění dalších efektů (rotace, vítr).
| Parametr | Hodnota |
|----------|----------|
| Dosah | 29 000 m |
| Max. výška | 11 500 m |
| Doba letu | 98.4 s |
| Poč. rychlost | 600 m/s |
Co je důležité
- Numerické metody jsou nezbytné pro realistické modely s proměnnými Cx(M) a ρ(h).
- RK4 zajišťuje přesnost při dt=0.01 s pro balistické úlohy.
- Zákon z roku 1943 přesně modeluje transsonický režim (M≈0.8–1.0).
- Exponenciální ρ(h) poskytuje rozumnou aproximaci do 15 km.
— Editorial Team
Zatím žádné komentáře.