ballistic_rocket

Расчёт точки падения головной части.

Дано

• геодезические координаты точки старта (долгота, широта, высота) баллистической ракеты (БР);
• азимут пуска;
• начальная вертикальная скорость выброса БР из шахты;
• временная программа управления по углу тангажа для 1, 2 и 3 ступеней БР;
• временная таблица изменения тяги двигательных установок для 1, 2 и 3 ступеней БР;
• временная таблица изменения массы 1, 2 и 3 ступеней БР;
• площадь среза сопла 1, 2 и 3 ступеней БР;
• площадь Миделя 1, 2 и 3 ступеней БР;
• масса БГ;
• площадь Миделя БГ;
• баллистический коэффициент БГ.

Задачи

• моделирование движения БР на основе системы дифференциальных уравнений в начальной стартовой системе координат (НССК) на активном атмосферном и внеатмосферном участках траектории;
• моделирование движения БГ на основе системы дифференциальных уравнений в НССК на пассивном внеатмосферном участке траектории;
• моделирование движения БГ на основе системы дифференциальных уравнений в НССК на пассивном атмосферном участке траектории;
• расчет трассы полета БР (геодезической широты и долготы) от момента старта до точки падения и ее отображение на фоне карты Земли.

Алгоритм

Для выполнения поставленных задач задание было разбито на следующие этапы:

• Моделирование движения БР на активном атмосферном участке
• Моделирование движения БР на активном внеатмосферном участке
• Моделирование движения БГ на пассивном внеатмосферном участке
• Моделирование движения БГ на пассивном атмосферном участке
• Визуализация

Активный атмосферный участок

Для моделирования БР будет использован метод Рунге-Кутты 4 порядка с шагом в 3 микросекунды. Будут использованы следующие ДУ в ГОСК:
  $\dot{r}=v$
  $\dot{v}=\frac{1}{m}(P-F_D-F_g)$
Где:
$P$ – Мощность двигателя в текущий момент времени
$F_D = \frac{1}{2} Cd A v^2$ – Сила лобового сопротивления
$F_g = \frac{G M_E m r}{|r|^3}$ – Сила гравитационного притяжения
$Cd$ – Коэффициент сопротивления текущей ступени при текущей скорости
$A$ – Площадь Миделя текущей ступени

Масса и мощность двигателя дана в таблице. Для промежуточных значений применяется линейная интерполяция с разрывами в точках отделения ступеней. Направление мощности определяется по временной программе управления по углу тангажа.

Для расчёта Cd используется таблица коэффициентов при разных высотах и числах Маха, которая интерполируется для промежуточных значений. Скорость звука, необходимая для расчёта числа Маха, вычисляется согласно ГОСТ 4401-81. Плотность воздуха также вычисляется согласно данному ГОСТ.

Активный внеатмосферный участок

Единственное отличие от атмосферного участка заключается в том, что в ДУ для скорости отсутствует сила лобового сопротивления.

Пассивный внеатмосферный участок

Двигатели прекращают свою работу, остаётся только боеголовка. Вне атмосферы на неё действует только сила притяжения Земли. Дифференциальные уравнения будут выглядеть следующим образом:
  $\dot{r}=v$
  $\dot{v}=-\frac{1}{m} F_g$

Пассивный атмосферный участок

Боеголовка входит в атмосферу, на неё начинает сила лобового сопротивления.
  $\dot{r}=v$
  $\dot{v}=\frac{1}{m}(-F_D-F_g)$

При столкновении боеголовки с Землей работа программы прекращается.

Визуализация

Для визуализации был написан скрипт на Python, строящий трассу и траекторию полёта по выводу программы.

Результаты, примеры работы программы

Программа была запущена для входных данных, указанных в файле datafiles/launch_ini.txt и нескольких файлов, описывающих изменение тяги, массы и угла тангажа, а также включающих в себя коэффициенты сопротивления для разных ступеней в зависимости от высоты и числа Маха.
Результаты моделирования:

• Время полёта: 2131.11 с.
• Точка старта: [28.889 283.902 100] - [B L H]
• Точка падения: [-1.12955 -9.23319 -2.25659] - [B L H]

Трасса полёта:

График изменения высоты относительно дальности полёта