Движение спускаемого аппарата в атмосфере
Разработать программу интегрирования уравнений движения спускаемого аппарата (СА) в атмосфере (баллистический спуск). Процесс интегрирования должен завершаться при достижении спускаемым аппаратом нулевой высоты.
Уравнения движения записать, используя декартовы координаты СА относительно геоцентрической системы координат \(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\), м/c | 7000 |
| \(\mu\), м\(^3/c^{-2}\) | 398600.4415\(\cdot 10^9\) |
- Проинтегрируйте уравнение движения спускаемого аппарата до достижения нулевой высоты.
- Постройте график траектории движения спускаемого аппарата.
- Постройте зависимость скорости движения спускаемого аппарата от времени.
- Постройте зависимость скорости движения спускаемого аппарата от высоты.
- Постройте зависимость высоты спускаемого аппарата от времени.
- Постройте зависимость скоростного напора от высоты.
- Найдите высоту, на которой скоростной напор достигает максимального значения.
- Постройте график изменения перегрузки.
Методические указания
Для определения плотности воздуха используйте функцию интерполирования по таблице плотности воздуха по стандартной атмосфере в диапазоне высот от -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
