Modelování vnější balistiky: analýza a Python

Analytické a numerické modelování vnější balistiky: tři klíčové případy

Úloha vnější balistiky popisuje rovinný pohyb tělesa působením gravitace, které je vystřeleno pod úhlem k horizontu s počáteční rychlostí. Zvažujeme tři případy: absence odporu vzduchu, lineární odpor podle Stokesova zákona a kvadratický odpor podle Newtonova zákona. Analytická řešení jsou doplněna numerickými metodami pro realistické modelování.

V prvním případě bez odporu je trajektorie parabola. Pro v₀=20 m/s a α=30° se dosah a výška určí jednoduchými vzorci:

x(t) = v₀ cos(α) t

y(t) = v₀ sin(α) t - (g t²)/2

Lineární odpor: Stokesův zákon

Při lineárním odporu je síla úměrná rychlosti: F = -k v. Diferenciální rovnice:

m dv_x/dt = -k v_x

m dv_y/dt = -m g - k v_y

Řešení dává exponenciální závislosti rychlostí:

v_x(t) = v₀ cos(α) e^(-kt/m)

v_y(t) = (v₀ sin(α) + (m g)/k) e^(-kt/m) - (m g)/k

Integrací podle souřadnic lze získat y(x). Pro v₀=70 m/s, α=55°, k=0,04 kg/s, m=0,5 kg je trajektorie asymetrická s menším dosahem než v ideálním případě.

Kvadratický odpor: Newtonův zákon

Při vysokých rychlostech dominuje kvadratický odpor: F = -k v². Soustava diferenciálních rovnic:

dv_x/dt = -k v v_x

dv_y/dt = -g - k v v_y

kde v = √(v_x² + v_y²). Leonard Euler vyřešil tuto soustavu analyticky v 18. století. Moderní přístup spočívá v numerické integraci pomocí Runge-Kutta řádu 4–5.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

class BallisticsSimulator:
    def __init__(self):
        self.g = 9.81
        self.rho = 1.225

    def calculate_k(self, diameter, mass, Cd=0.3):
        A = np.pi * (diameter / 2) ** 2
        k = (Cd * self.rho * A) / (2 * mass)
        return k

    def equations_of_motion(self, t, state, k):
        x, vx, y, vy = state
        v = np.sqrt(vx**2 + vy**2)
        ax = -k * v * vx
        ay = -self.g - k * v * vy
        return [vx, ax, vy, ay]

    def simulate(self, v0, angle_deg, x0=0, y0=0,
              diameter=0.1, mass=10, Cd=0.3,
              t_max=100, dt=0.01):
        k = self.calculate_k(diameter, mass, Cd)
        angle_rad = np.radians(angle_deg)
        vx0 = v0 * np.cos(angle_rad)
        vy0 = v0 * np.sin(angle_rad)
        initial_state = [x0, vx0, y0, vy0]
        t_span = (0, t_max)
        t_eval = np.arange(0, t_max, dt)
        solution = solve_ivp(
            self.equations_of_motion,
            t_span,
            initial_state,
            args=(k,),
            t_eval=t_eval,
            method='RK45'
        )
        t = solution.t
        x = solution.y[0]
        y = solution.y[2]
        fall_index = np.where(y < 0)[0]
        if len(fall_index) > 0:
            fall_index = fall_index[0]
            t = t[:fall_index + 1]
            x = x[:fall_index + 1]
            y = y[:fall_index + 1]
        return t, x, y, k

    def plot_trajectory(self, t, x, y, k, v0, angle):
        plt.figure(figsize=(12, 8))
        plt.plot(x, y, 'b-', linewidth=2, label='Trajektorie')
        plt.scatter(x[0], y[0], color='red', s=100, zorder=5, label='Start')
        plt.scatter(x[-1], y[-1], color='green', s=100, zorder=5, label='Pád')
        plt.xlabel('Vzdálenost, m', fontsize=12)
        plt.ylabel('Výška, m', fontsize=12)
        plt.title(f'Trajektorie střely\nV₀ = {v0} m/s, úhel = {angle}°, k = {k:.6f}',
                  fontsize=14)
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.axis('equal')
        plt.show()

    def print_results(self, t, x, y):
        print("=== VÝSLEDKY MODELOVÁNÍ ===")
        print(f"Dosah letu: {x[-1]:.2f} m")
        print(f"Maximální výška: {max(y):.2f} m")
        print(f"Čas letu: {t[-1]:.2f} s")

# Příklad: v0=250 m/s, úhel=60°, d=0.2 m, m=20 kg

Pro v₀=250 m/s, α=60°, průměr 20 cm, hmotnost 20 kg, Cd=0.3 je koeficient k vypočítán jako k = (Cd ρ A)/(2m), kde A je průřez. Funkce solve_ivp s metodou RK45 zajišťuje přesnost až do okamžiku dopadu (y=0).

Parametry a výsledky modelování

Klíčové vstupní parametry:

• Počáteční rychlost v₀: 20–250 m/s
• Úhel vystřelení α: 30–60°
• Hmotnost m: 0,5–20 kg
• Průměr: 10–20 cm
• Cd: 0,3 (typická hodnota)
• k: vypočítáno automaticky

Výsledky pro kvadratický případ:

• Dosah letu: ~10–20 km (závisí na parametrech)
• Maximální výška: 5–15 km
• Čas letu: 50–100 s
• Konečná rychlost: blíží se termální rychlosti

Srovnání případů ukazuje zkrácení dosahu o 30–50 % při započtení odporu vzduchu.

Co je důležité

• Bez odporu — parabolická trajektorie, maximální dosah při α=45°.
• Stokesův zákon platí pro nízké rychlosti (<30 m/s), dává exponenciální tlumení.
• Newtonův zákon je standard pro artilerii, vyžaduje numerické řešení diferenciálních rovnic.
• Python kód s scipy.integrate je univerzální, snadno přizpůsobitelný reálným hodnotám Cd a ρ.

— Editorial Team