Движение спускаемого аппарата в атмосфере
Простейшая плоская модель движения СА в атмосфере (Nguyen X. Vinh, Adolf Busemann, Robert D. Culp Hypersonic and planetary entry flight mechanics) The Unevrsity of Michigan Press
Спускаемый аппарат, рассматриваемый как материальная точка с заданной массой, коэффициентом лобового аэродинамического сопротивления \(C_d\) и коэффициентом подъемной силы \(C_l\) начинает движение на круговой орбите высотой 100 км. Уравнения движения СА имеют вид: \[\frac{dV}{dt} = - \frac{C_d S_m q}{m} - g \sin \theta,\] \[V \frac{d \theta}{dt} = \frac{C_l S_m q}{m} - \left(g-\frac{V^2}{r}\right) \cos \theta,\] \[\frac{dr}{dt} = V \sin \theta\]
где \(V\) - скорость СА; \(\vartheta\) – угол наклона траектории (угол между вектором скорости и линией местного горизонта); \(r\) – расстояние от центра Земли до центра масс СА; \(g\) – ускорение свободного падения на высоте \(h\); \(m\) – масса СА; \(S_m\) – характерная площадь (площадь Миделя); \(C_d\) – аэродинамический коэффициент лобового сопротивления; \(C_l\) – аэродинамический коэффициент подъёмной силы; \(q\) – скоростной напор, зависящий от плотности воздуха и скорости движения СА: \[q = \frac{\rho V^2}{2}\]
| Параметр | Значение |
|---|---|
| Масса СА, кг | 1000 |
| \(C_d\) | 1.2 |
| \(C_l\) | 0.0 |
| \(S_m\), м\(^2\) | 5 |
Код на языке Python
Подключение необходимых библиотек
from collections import namedtuple
import scipy
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.pylab as pylab
params = {'legend.fontsize': 14, 'figure.figsize': (10, 7), 'axes.labelsize': 14,
'axes.titlesize':14, 'xtick.labelsize':14,'ytick.labelsize':14}
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.pylab as pylab
params = {'legend.fontsize': 14, 'figure.figsize': (10, 7), 'axes.labelsize': 14,
'axes.titlesize':14, 'xtick.labelsize':14,'ytick.labelsize':14}
pylab.rcParams.update(params)
Модель изменения плотности воздуха в диапазоне от 0 до 100 км с погрешностью 1,5% (В. А. ЯРОШЕВСКИЙ “АППРОКСИМАЦИЯ МОДЕЛИ СТАНДАРТНОЙ АТМОСФЕРЫ” Ученые записки ЦАГИ, т. XL, №3, 2009.)
a = np.array([-6.3759,-7.3012,-1.1817])
b = np.array([-0.4754,-0.0096,-0.0068,-0.0120,0.0042]);
c = np.array([ 0.1803, 0.0872,-0.0153, 0.0145,0 ]);
def rho(h):
# Функция вычисления плотности воздуха
# высота h задаётся в километрах
# результат -- плотность кг/м3
x = h/50.0-1
sa = a[0] + a[1]*x + a[2]*x*x
sbc= sum( (b[i]*np.cos((i+1)*np.pi*x) + c[i]*np.sin((i+1)*np.pi*x) for i in range(5)) )
return np.exp(sa + sbc)
Константы
# Гравитационный параметр Земли
mu = 398600.4415e9
# Радиус Земли
Re = 6371000.0
Функция вычисления ускорения свободного падения м/с\(^2\) на высоте h (в метрах)
def g_acc(h):
#
return mu/(Re+h)**2
Функция правых частей
Функция правых частей дифференциальных уравнений: вычисление правых частей дифференциальных уравнений для момента времени t и вектора состояния y = [r(t), v(t), theta(t)]:
def dydt(t, y, p):
# Радиус-вектор точки в момент времени t
r = y[0]
# Скорость в момент времени t
v = y[1]
# Угол наклона траекториии в момент времени t
theta = y[2]
# Высота (км)
h = (r - Re)*0.001
# Скоростной напор Н/м^2
q = rho(h)*v*v/2
# Ускорение свободного падения
g = g_acc(r - Re)
# dv/dt =
dv = - q*p.CD*p.Sm/p.mass - g*np.sin(theta)
# dtheta/dt =
dtheta= (q*p.CL*p.Sm/p.mass - (g - v*v/r)*np.cos(theta))/v
# dr/dt =
dr = v*np.sin(theta)
# ds/dt
ds = (Re/r)*v*np.cos(theta)
return (dr,dv,dtheta,ds)
Функция для остановки процесса интегрирования при достижении нулевой высоты
def event_h_eq_0(t, y):
# Функция-"детектор", передаваемая в интегратор (параметр events),
# для определения времени достижения нулевой высоты и
# остановки процесса интегрирования
return y[0]-Re
# функция определяется условие h = 0 при движении "вниз"
event_h_eq_0.direction = -1
# функция-детектор активна
event_h_eq_0.terminal = True
# функция возвращает высоту
Параметры системы
params = namedtuple("params", "CD CL mass")
# Масса тела
params.mass = 60000.0
# Аэродинамические коэфиициенты,
# которые в общем случае зависят от числа Маха,
# угла атаки
params.CD = 1.4 # Коэффицент лобового сопротивления
params.CL = 0.0 # Коэффицент подъёмной силы
# Площадь миделя
params.Sm = np.pi*3**2/4
Начальные условия
# Начальные условия
h0 = 120e3 # Начальная высота [м]
v0 = 2800 # Начальная скорость [м/c]
theta0 = np.deg2rad(30.0) # начальный угол наклона траектории [радиан]
Запуск процесса интегрирования
sol = solve_ivp(lambda t,y: dydt(t,y,params), [0, 1000.0], [Re+h0, v0, theta0, 0], method='LSODA', events = [event_h_eq_0], rtol = 1e-9)
Результаты
График изменения высоты от времени
plt.plot(sol.t,(sol.y[0]-Re)*0.001);
plt.xlabel('t, c');plt.ylabel('h, км');
Зависимость скоростного напора от высоты
plt.plot((sol.y[0]-Re)*0.001,0.001*rho((sol.y[0]-Re)*0.001)*sol.y[1]*sol.y[1]/2.0);
plt.xlabel('h, км');plt.ylabel('q, кН/м$^2$');
