Движение спускаемого аппарата в атмосфере

Разработать программу интегрирования уравнений движения спускаемого аппарата (СА) в атмосфере (баллистический спуск). Процесс интегрирования должен завершаться при достижении спускаемым аппаратом нулевой высоты.

Уравнения движения записать, используя декартовы координаты СА относительно геоцентрической системы координат \(OX_oY_o\): \[\boldsymbol{r}= \begin{bmatrix} x \\ y \end{bmatrix}.\]

В начальный момент времени СА находится в точке с координатами \((0,y_0)\) и движется со скоростью \(V_0\).

Движение СА происходит в центральном гравитационном поле Земли. На СА действуют сила притяжения \(\boldsymbol G\), определяемая выражением: \[\boldsymbol G = - \frac{\mu m}{r^3} \boldsymbol{r},\]

где \(\mu\) – гравитационный параметр Земли; и сила лобового аэродинамического сопротивления \(\boldsymbol{F}_d\): \[\boldsymbol{F}_d = - \boldsymbol{e}_v C_d q S_m,\]

где \(\boldsymbol{e}_v\) – единичный вектор, определяющий направление скорости СА \[\boldsymbol{e}_v = \frac{\boldsymbol{V}}{|\boldsymbol{V}|}, \quad \boldsymbol{V} = \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix},\]

где \(C_d\) – коэффициент лобового сопротивления, \(S_m\) – характерная площадь (площадь миделя), \(q\) – скоростной напор, зависящий от плотности воздуха и скорости движения СА: \[q = \frac{\rho V^2}{2}\]

При определении скорости движения спускаемого аппарата в атмосфере вращение Земли не учитывается.

ПараметрЗначение
Масса СА, кг3000
\(C_d\)1.6
\(S_m\), м\(^2\)15
\(R_e\), м6371000
\(y_0\), м6461000
\(V_0\), м/c7000
\(\mu\), м\(^3/c^{-2}\)398600.4415\(\cdot 10^9\)
  1. Проинтегрируйте уравнение движения спускаемого аппарата до достижения нулевой высоты.
  2. Постройте график траектории движения спускаемого аппарата.
  3. Постройте зависимость скорости движения спускаемого аппарата от времени.
  4. Постройте зависимость скорости движения спускаемого аппарата от высоты.
  5. Постройте зависимость высоты спускаемого аппарата от времени.
  6. Постройте зависимость скоростного напора от высоты.
  7. Найдите высоту, на которой скоростной напор достигает максимального значения.
  8. Постройте график изменения перегрузки.

Методические указания

  • Для определения плотности воздуха используйте функцию интерполирования по таблице плотности воздуха по стандартной атмосфере в диапазоне высот от -1000 до 100000 м над уровнем моря. Таблица должна загружаться из текстового файла (см. функции dlmread, interp1) до запуска процесса интегрирования.

  • Для определения момента окончания процесса интегрирования используйте функцию Event Locator, передаваемую в функцию численного интегрирования:

function [value,isterminal,direction] = h0EventsFcn(t,q)
  h = ... % вычисление высоты  
  value      =  h; % Отслеживаемая величина (высота)
  isterminal =  1; % Остановить процесс интегрирования
  direction  = -1; % Процесс интегрирования останавливается при уменьшении высоты dh/dt<0, h=0
end

© 2022. All rights reserved.

Powered by Hydejack v9.1.6