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

2026-04-05 в 11:01, admin, рубрики: баллистика, баллистическое движение, математика, механика, плотность газа, сопротивление воздуха, физика, численное моделирование, численные методы

В этой статье мы продолжим решать задачу внешней баллистики и разберём четвёртый и пятый случаи. Они более общие.

Мы учтём зависимость коэффициента лобового сопротивления от скорости, а также зависимость плотности воздуха от высоты.

Показано, что в этом случае аналитическое решение невозможно, поэтому нужно использовать численные методы.

Написан алгоритм численного решения данной задачи.

Сейчас напишем на Python код для моделирования полёта снаряда с учётом вышеизложенного:

(начальная скорость 600 м/с, угол 55 градусов, масса снаряда 50 кг)

import numpy as np
import matplotlib.pyplot as plt

# Константы
g = 9.81  # м/с²
m = 50.0  # масса тела, кг
S = 0.01  # площадь поперечного сечения, м²
rho0 = 1.225  # плотность воздуха на уровне моря, кг/м³
H = 8000  # масштабная высота, м
a = 343  # скорость звука, м/с
v0 = 600  # начальная скорость, м/с
theta = np.radians(55)  # начальный угол, рад
dt = 0.01  # шаг интегрирования, с
t_end = 98.4  # конечное время, с

# Начальные условия
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)

# Массивы для хранения результатов
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)


# Функция для вычисления 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


# Функция плотности воздуха
def rho(h):
    return rho0 * np.exp(-h / H)


# Система ОДУ
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])


# Метод Рунге-Кутты 4-го порядка
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


# Интегрирование
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)

# Визуализация
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Траектория')
plt.xlabel('X, м')
plt.ylabel('Y, м')
plt.title('Траектория полёта тела')
plt.grid(True)
plt.legend()
plt.show()

# Вывод результатов
print(f"Дальность полёта: {X[-1]:.2f} м")
print(f"Максимальная высота: {np.max(Y):.2f} м")

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

Видно, то снаряд в таком случае пролетит 29 км и залетит на высоту 11,5 км.

При этом было допущение, что зависимость плотности воздуха от высоты экспоненциальная

(это, конечно не так, это требует написания отдельной статьи).

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

Написан код, который моделирует траекторию полёта тела в этом случае.

Продолжение исследования будет далее, в следующих статьях......

Литература:

Окунев Б. Н. «Решение основной задачи внешней баллистики при квадратичном законе сопротивления воздуха» (1932).
Окунев Б. Н. «Основная задача внешней баллистики и аналитические методы её решения» (1934).
Дмитриевский А. А., Лысенко Л. Н. «Внешняя баллистика» (2005).
Лысенко А. Н. «Внешняя баллистика» (2024).
Шапиро Я. М. «Внешняя баллистика» (1946).
Беляева С. Д. «Внешняя баллистика с примерами и задачами» (2023).