Modelowanie balistyki zewnętrznej: analiza i Python

Analiza i symulacja balistyki zewnętrznej: trzy kluczowe przypadki

Zadanie balistyki zewnętrznej opisuje ruch ciała w płaszczyźnie pod wpływem siły grawitacji, rzucanego pod kątem do poziomu z początkową prędkością. Rozważa się trzy przypadki: brak oporu powietrza, liniowy opór zgodnie z prawem Stokesa oraz kwadratowy opór zgodnie z prawem Newtona. Rozwiązania analityczne uzupełnia się metodami numerycznymi, aby uzyskać realistyczne symulacje.

W pierwszym przypadku bez oporu trajektoria ma kształt paraboli. Dla v₀=20 m/s i α=30° zasięg i wysokość maksymalna wyznacza się prostymi wzorami:

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

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

Liniowy opór: prawo Stokesa

Przy liniowym oporze siła jest proporcjonalna do prędkości: F = -k v. Równania różniczkowe:

m dv_x/dt = -k v_x

m dv_y/dt = -m g - k v_y

Rozwiązanie daje wykładnicze zależności prędkości:

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

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

Całkowanie po współrzędnych pozwala wyznaczyć y(x). Dla v₀=70 m/s, α=55°, k=0,04 kg/s, m=0,5 kg trajektoria jest asymetryczna i znacznie krótsza niż w przypadku idealnym.

Kwadratowy opór: prawo Newtona

Na dużych prędkościach dominuje opór kwadratowy: F = -k v². Układ równań różniczkowych:

dv_x/dt = -k v v_x

dv_y/dt = -g - k v v_y

gdzie v = √(v_x² + v_y²). Leonard Euler rozwiązał ten układ analitycznie w XVIII wieku. Współczesny podejście to całkowanie numeryczne metodą Rungego-Kutty rzę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='Trajektoria')
        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='Upadek')
        plt.xlabel('Odległość, m', fontsize=12)
        plt.ylabel('Wysokość, m', fontsize=12)
        plt.title(f'Trajektoria pocisku\nV₀ = {v0} m/s, kąt = {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("=== WYNIKI SYMULACJI ===")
        print(f"Zasięg lotu: {x[-1]:.2f} m")
        print(f"Wysokość maksymalna: {max(y):.2f} m")
        print(f"Czas lotu: {t[-1]:.2f} s")

# Przykład: v0=250 m/s, kąt=60°, d=0.2 m, m=20 kg

Dla v₀=250 m/s, α=60°, średnica 20 cm, masa 20 kg, Cd=0.3 współczynnik k oblicza się jako k = (Cd ρ A)/(2m), gdzie A to pole przekroju. Funkcja solve_ivp z metodą RK45 zapewnia dokładność aż do momentu upadku (y=0).

Parametry i wyniki symulacji

Kluczowe parametry wejściowe:

• Początkowa prędkość v₀: 20–250 m/s
• Kąt rzutu α: 30–60°
• Masa m: 0,5–20 kg
• Średnica: 10–20 cm
• Cd: 0,3 (typowa wartość)
• k: obliczany automatycznie

Wyniki dla przypadku kwadratowego oporu:

• Zasięg lotu: ~10–20 km (zależy od parametrów)
• Wysokość maksymalna: 5–15 km
• Czas lotu: 50–100 s
• Prędkość końcowa: zbliżona do prędkości granicznej

Porównanie przypadków pokazuje skrócenie zasięgu o 30–50% przy uwzględnieniu oporu powietrza.

Co warto wiedzieć

• Bez oporu — trajektoria paraboliczna, maksymalny zasięg przy α=45°.
• Prawo Stokesa stosuje się do niskich prędkości (<30 m/s), daje wykładnicze tłumienie.
• Prawo Newtona — standard w artylerii, wymaga rozwiązania numerycznego równań różniczkowych.
• Kod w Pythonie z scipy.integrate jest uniwersalny, łatwo dostosowalny do rzeczywistych wartości Cd i ρ.

— Editorial Team