Задача внешней баллистики. Третья часть

2026-04-08 в 12:00, admin, рубрики: баллистика, баллистическое движение, дифференциальные уравнения, кривизна, математика, механика, сопротивление воздуха, траектория, физика, численные методы

В этой статье мы завершим решение задачи внешней баллистики разбором шестого и седьмого случаев. В них мы учтём уменьшение гравитации с высотой, а также кривизну Земли.

Реализуем всё вышенапечатанное на Python:

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

# Константы
g0 = 9.80665  # м/с² — ускорение свободного падения на уровне моря
RE = 6371000  # м — радиус Земли
rho0 = 1.225  # кг/м³ — плотность воздуха на уровне моря
H = 8000  # м — масштабная высота

# Скорость звука как функция высоты (упрощённая модель)
def sound_speed(h):
    if h < 11000:  # Тропосфера
        return 340 - 0.0065 * h
    else:  # Стратосфера и выше
        return 295

# Параметры снаряда
m = 2000  # кг — масса
S = 0.3  # м² — площадь миделева сечения

# Начальная скорость и угол
v0 = 4000  # м/с
theta0 = 50  # градусов

def cx(M):
    """Коэффициент лобового сопротивления по закону 1943 года"""
    if M < 0.8:
        return 0.15 + 0.25 * M
    elif M < 1.2:
        return 0.4 - 0.3 * (M - 0.8)
    elif M<4:
        return 0.28 + 0.05 * (M - 1.2)
    else:
        return 1

def rho(h):
    """Плотность воздуха с экспоненциальным убыванием"""
    return rho0 * np.exp(-h / H)

def g(h):
    """Ускорение свободного падения с высотой"""
    return g0 * (RE / (RE + h))**2

def equations(t, state):
    """
    Система дифференциальных уравнений движения
    state = [r, phi, vr, vphi]
    r — расстояние от центра Земли
    phi — полярный угол
    vr — радиальная скорость
    vphi — азимутальная скорость
    """
    r, phi, vr, vphi = state
    h = r - RE  # высота над поверхностью

    # Проверка выхода за пределы расчёта
    if h < -100:  # снаряд ушёл под землю
        return [0, 0, 0, 0]
    if h > 1e6:  # высота более 1000 км — сопротивление пренебрежимо мало
        # Упрощённые уравнения без сопротивления
        dvr_dt = -g(h) + vphi**2 / r
        dvphi_dt = -vr * vphi / r
        return [vr, vphi / r, dvr_dt, dvphi_dt]

    # Скорость и число Маха
    v = np.sqrt(vr**2 + vphi**2)
    if v == 0:
        M = 0
    else:
        a = sound_speed(h)
        M = v / a

    # Плотность воздуха и коэффициент сопротивления
    rho_h = rho(h)
    cx_val = cx(M)

    # Силы сопротивления
    F_drag = 0.5 * rho_h * v**2 * S * cx_val
    Fr = -F_drag * vr / v  # радиальная компонента
    Fphi = -F_drag * vphi / v  # азимутальная компонента

    # Ускорения
    dvr_dt = Fr / m - g(h) + vphi**2 / r
    dvphi_dt = Fphi / m - vr * vphi / r

    return [vr, vphi / r, dvr_dt, dvphi_dt]

def event_ground(t, state):
    """Событие — достижение поверхности Земли"""
    r = state[0]
    return r - RE  # ноль при r = RE
event_ground.terminal = True  # останавливаем интегрирование
event_ground.direction = -1  # только при уменьшении r

# Начальные условия
r0 = RE  # начальное расстояние от центра Земли
phi0 = 0  # начальный угол
vr0 = v0 * np.sin(np.radians(theta0))  # радиальная скорость
vphi0 = v0 * np.cos(np.radians(theta0))  # азимутальная скорость

initial_state = [r0, phi0, vr0, vphi0]

# Временной интервал (максимальное время)
t_span = (0, 1000)  # 0–1000 секунд (достаточно для большинства случаев)

# Решение системы ОДУ с событием
solution = solve_ivp(
    equations,
    t_span,
    initial_state,
    method='RK45',
    events=event_ground,
    dense_output=True
)

# Извлечение результатов
t = solution.t
r = solution.y[0]
phi = solution.y[1]
vr = solution.y[2]
vphi = solution.y[3]

# Преобразование в декартовы координаты для визуализации
x = np.sqrt(vr**2 + vphi**2)*t * np.cos(np.radians(phi))

# Расчёт дополнительных параметров
h = r - RE
v_total = np.sqrt(vr**2 + vphi**2)
M_history = [v_total[i] / sound_speed(h[i]) if v_total[i] > 0 else 0 for i in range(len(v_total))]

# Поиск точки падения (из события)
if solution.t_events[0].size > 0:
    fall_time = solution.t_events[0][0]
    # Интерполяция для точного определения координат в момент падения
    state_fall = solution.sol(fall_time)
    r_fall = state_fall[0]
    phi_fall = state_fall[1]
    range_m = (r_fall - RE) * phi_fall  # приближённо для малых углов
    max_height = np.max(h)  # максимальная высота
else:
    fall_time = t[-1]
    range_m = x[-1]
    max_height = np.max(h)

# Визуализация
plt.figure(figsize=(14, 10))

# Траектория
plt.subplot(2, 2, 1)
plt.plot(x/1000,h/1000, label='Траектория', linewidth=2)
plt.xlabel('X, км')
plt.ylabel('Y, км')
plt.title('Траектория полёта')
plt.grid(True, alpha=0.7)
plt.axis('equal')
plt.legend()

# Высота во времени
plt.subplot(2, 2, 2)
plt.plot(t, h/1000, label='Высота', linewidth=2, color='orange')
plt.xlabel('Время, с')
plt.ylabel('Высота, км')
plt.title('Высота полёта')
plt.grid(True, alpha=0.3)
plt.legend()

# Скорость во времени
plt.subplot(2, 2, 3)
plt.plot(t, v_total, label='Скорость', linewidth=2, color='green')
plt.xlabel('Время, с')
plt.ylabel('Скорость, м/с')
plt.title('Скорость полёта')
plt.grid(True, alpha=0.3)
plt.legend()

# Число Маха во времени
plt.subplot(2, 2, 4)
plt.plot(t, M_history, label='Число Маха', linewidth=2, color='red')
plt.xlabel('Время, с')
plt.ylabel('Число Маха')
plt.title('Число Маха полёта')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

Результат моделирования можно увидеть ниже:

Таким образом, исходная задача решена, создана максимально реалистичная математическая модель в общем виде, учитывающая очень многие параметры и факторы, влияющие на движение тела в атмосфере Земли, и предложен алгоритм, который её решает за разумное время с большой точностью.

Литература: