В данной статье приведен один из возможных вариантов математического обеспечения для исследования динамики неуправляемого движения спускаемых аппаратов, созданное на основе обобщения опыта проектирования, отработки и эксплуатации аппаратов типа "Венера" и "Марс" второго и третьего поколений. Скрипт рассчитывающий параметры траектории написан в среде математического моделирования Engee (https://engee.com), ссылка на данный скрипт на портале Engee: https://engee.com/community/ru/catalogs/projects/modelirovanie-posadki-komicheskogo-apparata-na-veneru . Приводятся обобщенные данные, характеризующие динамику спуска в атмосферах планет типа "Венера 9-16", "Вега-1,2" и "Марс-3,6". При данном моделирования рассматривается процесс посадки на Венеру и заданы параметры её атмосферы.
Подключаемые библиотеки
using Plots
using Pkg
using Polynomials
using LinearAlgebra
#Константы характеризующие атмосферу планеты
# высота (метры), температура (по данным), давление (в мм рт.ст.)
data = [
0 462 92.10;
5 424 66.65;
10 385 47.39;
15 348 33.04;
20 306 22.52;
25 264 14.93;
30 222 9.851;
35 180 5.917;
40 143 3.501;
45 110 1.979;
50 75 1.066;
55 27 0.5314;
60 -10 0.2357;
65 -30 0.09765;
70 -43 0.03690;
80 -76 0.004760;
90 -104 0.0003736;
100 -112 0.00002660
]
# Высота в метрах
h = data[:, 1] * 1000
# Температура в Кельвинах
T = data[:, 2] .+ 273
# Атмосферное давление в мм рт. ст.
P = data[:, 3]
# Молярная масса воздуха (г/моль)
M = 28.98
# Универсальная газовая постоянная (Дж/(моль·К))
R = 8.314
# Расчет плотности воздуха
ro = P .* M ./ (R .* T)
x = h
y = ro
Интерполяция параметров атмосферы с помощью МНК
Так как параметры атмосферы сильно зависят от высоты дискретной таблицы с константами недостаточно для расчета, для получения непрерывной функции характеризующей атмосферы будем использовать метод наименьших квадратов шестого порядка
n = length(x);
Sx = sum(x)
Sx2 = sum(x.^2)
Sx3 = sum(x.^3)
Sx4 = sum(x.^4)
Sx5 = sum(x.^5)
Sx6 = sum(x.^6)
Sx7 = sum(x.^7)
Sx8 = sum(x.^8)
Sx9 = sum(x.^9)
Sx10 = sum(x.^10)
Sx11 = sum(x .^11)
Sx12 = sum(x .^12)
Sy = sum(y)
Sxy = sum(x .*y)
Sx2y = sum(x .^2 .*y)
Sx3y = sum(x.^3 .*y)
Sx4y = sum(x.^4 .*y)
Sx5y = sum(x.^5 .*y)
Sx6y = sum(x.^6 .*y)
#F(x) = a0+a1*x+a2*x^2+a3*x^3+a4*x^4+a5*x^5+a6*x^6;
# поиск коэффициентов
A = [n Sx Sx2 Sx3 Sx4 Sx5 Sx6;
Sx Sx2 Sx3 Sx4 Sx5 Sx6 Sx7;
Sx2 Sx3 Sx4 Sx5 Sx6 Sx7 Sx8;
Sx3 Sx4 Sx5 Sx6 Sx7 Sx8 Sx9;
Sx4 Sx5 Sx6 Sx7 Sx8 Sx9 Sx10;
Sx5 Sx6 Sx7 Sx8 Sx9 Sx10 Sx11;
Sx6 Sx7 Sx8 Sx9 Sx10 Sx11 Sx12];
b = [Sy,
Sxy,
Sx2y,
Sx3y,
Sx4y,
Sx5y,
Sx6y ];
Ai = inv(A)
A2 = Ai*b; # коэффициенты полинома
#Функция описывающая параметры атмосферы
function F(X,A2)
return A2[1] + A2[2]*X+ A2[3]*X^2+A2[4]*X^3+A2[5]*X^4+A2[6]*X^5+A2[7]*X^6
end
y2 = zeros(length(x),1)
for i = 1:length(x)
y2[i] = F(x[i],A2)
end
#Построение исходных точек
plot(x,y,title="Плотность атмосферы (замер)")
plot(x,y2,title="Плотность атмосферы (интерполяция)")
По данным графикам видно, что плотность атмосферы снижается с высотой.
Моделирование траектории движения спускаемого аппарата (кусочно)
В рамках данной работы предполагается что траектория спускаемого аппарата состоит из трех участков:
Пассивное торможение. Спускаемый аппарат входит в атмосферу планеты и под воздействием динамического трения у него снижается скорость движения. Данный режим длится до тех пор, пока скорость не позволит раскрыть парашют.
Пассивное торможение с раскрытым парашютом. Наиболее длинный участок траектории, спускаемый аппарат плавно теряет скорость и высоту.
Активное торможение. Завершающий участок траектории на котором происходит включение тормозного двигателя и торможение. При этом парашют остается открытым и отстреливается непосредственно у поверхности.
Для моделирования движения будет использоваться система дифференциальных уравнений которая описывает возмущенное баллистическое движение (с динамическим трением и реактивным движением). Результаты интегрирования данной системы ДУ с тремя различными наборами входных данных будут сопрягаться на основании данных о высота на которой будет находиться спускаемый аппарат. Вектор содержащий высоты на которых будет осуществляться сопряжение выглядит следующим образом [70000 50000], погрешность сопряжение составляет 100м.
Дифференциальное уравнение описывающее движение
function odefcn2(t, x, S, A2, F_torm)
k = F(x[4],A2) # коэффициент лобового сопротивления
m = 10 # масса тела
dxdt = zeros(4)
dxdt[1] = -x[1]*k/m*S # проекция скорости по X
dxdt[2] = -9.87 - k/m*x[2]+F_torm/m # проекция скорости по Y
dxdt[3] = x[1] # проекция координат по x
dxdt[4] = x[2] # проекция координат по y
return dxdt
end
odefcn2 (generic function with 1 method)
Функция в которой реализован метод Эйлера
function euler_step(f, t, x, dt, S,A2,F_torm)
dx = f(t, x, S, A2,F_torm) # вызов функции содержащей дифференциальное уравнение
return x + dt * dx # сумма за 1 шаг
end
Расчет точек первого участка траектории
# Начальные условия
t = 0.0
t_final = 200
dt = 0.05
# в рамках моделирования для решения ДУ время берем с запасом чтобы не потерять точки сопряжения участков траектории
# лишний массив данных будет отсекаться и сопрряжение с участком активного
# торможения будет проходить по высоте
H_sw1 = 70000; # высота раскрытия парашюта
H_sw2 = 50000; # высота начала торможения щитком
h_eps = 100; # погрешность определения высоты
V0 = 11200; # начальная скорость
fi = 20/180*pi; # угол под которым космический аппарат входит в атмосферу
Vx0 = V0*cos(fi); # начальная проекция скорости по Х
Vy0 = -V0*sin(fi); # начальная проекция скорости по Х
X0 = 0;
Y0 = 100000;
S = 3.0 # площадь лобового сопротивления спускаемого аппарата
# Вектор начальных условий
x = [Vx0, Vy0, X0, Y0]
F_torm = 0 # тяга двигателя активного торможения (на данном участке траекто��ии двигатель выкючен)
times = [] # массив времени
N = Int((t_final-t)/dt) # число элементов в массивах времени и параметров движения
VX1 = zeros(N,1) # пустые матрицы содержащие параметры траектории
VY1 = zeros(N,1)
X1 = zeros(N,1)
Y1 = zeros(N,1)
times1 = zeros(N,1)
j = 1 # обнуление счетчика для заполнения выходных массивов
while t <= t_final # цикл содержащий метод Эйера и заполнение матриц содержащих параметры траектории
x = euler_step(odefcn2, t, x, dt, S,A2,F_torm)
VX1[j] = x[1]
VY1[j] = x[2]
X1[j] = x[3]
Y1[j] = x[4]
times1[j] = t
j = j+1
t += dt
end
Расчет ускорений на первом участке траектории
AX1 = zeros(N,1) # пустые матрицы содержащие проекции ускорения и угол направления движения.
AY1 = zeros(N,1)
alfa = zeros(N,1)
for i = 1:(length(VX1)-1)
AX1[i] = (VX1[i+1]-VX1[i])/dt; # расчет проекций ускорения с помощью численного дифференциарования
AY1[i] = (VY1[i+1]-VY1[i])/dt;
alfa[i] = atan((Y1[i+1]-Y1[i])/(X1[i+1]-X1[i])); # направление вектора ускорения
end
Расчет точки окончания первого участка траектории
N_sw1 = 0 # начальное значение счетчика содержащего индекс массива на котором закончился первый участок
for i = 1:(length(VX1)-1)
if (Y1[i]< H_sw1+h_eps) && (Y1[i]> H_sw1-h_eps) # условия по которым определяем окончание первого участка движения с учетом погрешности определения высоты
N_sw1 = i #обновление значения счетчика
Break # выход из цикла при достижении условия окончания участка
end
end
print(N_sw1) # вывод найденного значения счетчика
155
Расчет точек второго участка траектории
Vx0_1 = VX1[N_sw1]; # начальная проекция скорости по Х
Vy0_1 = VY1[N_sw1]; # начальная проекция скорости по Y
X0_1 = X1[N_sw1]; # начальная проекция координаты по Х
Y0_1 = Y1[N_sw1]; # начальная проекция координаты по Y
# начальные условия для расчета второго участка берем из конечной точки первого участка
# вектор начальных условий для второго участка траектории
x = [Vx0_1, Vy0_1, X0_1, Y0_1]
t = 0 # сброс счетчика времени
S = 100; # площадь парашютов
N = Int((t_final-t)/dt) # число точек в массиве содержащем второй участок траектории
VX2 = zeros(N,1) # пустые массивы для хранения параметров траектории
VY2 = zeros(N,1)
X2 = zeros(N,1)
Y2 = zeros(N,1)
times2 = zeros(N,1)
F_torm = 0 # на втором участке траектории тормозной двигатель выключен
j = 1
while t <= t_final
x = euler_step(odefcn2, t, x, dt, S,A2,F_torm)
VX2[j] = x[1] #накопление массивов результатов интегрирования ДУ
VY2[j] = x[2]
X2[j] = x[3]
Y2[j] = x[4]
times2[j] = t
j = j+1
t += dt
end
Проекции ускорения на втором участке траектории
AX2 = zeros(N,1)
AY2 = zeros(N,1)
alfa2 = zeros(N,1)
for i = 1:(length(VX2)-1)
AX2[i] = (VX2[i+1]-VX2[i])/dt; # расчет проекций ускорения для второго участка траектории
AY2[i] = (VY2[i+1]-VY2[i])/dt;
alfa2[i] = atan((Y2[i+1]-Y2[i])/(X2[i+1]-X2[i]));
end
Расчет точки окончания второго участка траектории
N_sw2 = 0
for i = 1:(length(VX2)-1)
if (Y2[i]< H_sw2+h_eps) && (Y2[i]> H_sw2-h_eps) # критерий достижения окончания второго участка
N_sw2 = i
break
end
end
print(N_sw2) # вывод точки окончания второго участка
102
Расчет точек третьего участка траектории
Vx0_2 = VX2[N_sw2]; # начальная проекция скорости по Х
Vy0_2 = VY2[N_sw2]; # начальная проекция скорости по Х
X0_2 = X2[N_sw2];
Y0_2 = Y2[N_sw2];
# вектор начальных условий для третьего участка
x = [Vx0_2, Vy0_2, X0_2, Y0_2]
t = 0 # cброс счетчика времени
S = 100; # площадь лобового сопротивления (парашют всё еще раскрыт)
F_torm = 1300 # включение тормозного двигателя
N = Int((t_final-t)/dt)
VX3 = zeros(N,1) # пустые массивы для накопления результатов интегрирования для третьего учаастка
VY3 = zeros(N,1)
X3 = zeros(N,1)
Y3 = zeros(N,1)
times3 = zeros(N,1)
j = 1
while t <= t_final
x = euler_step(odefcn2, t, x, dt, S,A2,F_torm) # вызов функции метода эйлера
VX3[j] = x[1]
VY3[j] = x[2]
X3[j] = x[3]
Y3[j] = x[4]
times3[j] = t
j = j+1
t += dt
end
Расчет проекций ускорения для третьего участка траектории
AX3 = zeros(N,1) # пустые массивы для накопления проекций вектора ускорения
AY3 = zeros(N,1)
alfa3 = zeros(N,1)
for i = 1:(length(VX3)-1)
AX3[i] = (VX3[i+1]-VX3[i])/dt; # запись в массивы проекций вектора ускорения
AY3[i] = (VY3[i+1]-VY3[i])/dt;
alfa3[i] = atan((Y3[i+1]-Y3[i])/(X3[i+1]-X3[i]));
end
Точка столкновения с поверхностью
N_sw3 = 0
for i = 1:(length(VX3)-1)
if (Y3[i]< 0+h_eps) && (Y3[i]> 0-h_eps) # критерий по которому определяем что поверхность достигнута (с учетом допуска)
N_sw3 = i
break
end
end
print(N_sw3) # вывод точки окончания третьего участка
401
После того как построены участки траектории соберем их в один массив
X_all = vcat(X1[1:N_sw1],X2[1:N_sw2],X3[1:N_sw3]) # проекции координаты Х для всех трех участках
Y_all = vcat(Y1[1:N_sw1],Y2[1:N_sw2],Y3[1:N_sw3])# # проекции координаты Y для всех трех участках
V_all = zeros(length(X_all)-1,1)# пустой массив для хранения вектора скорости
A_all = zeros(length(X_all)-1,1)# пустой массив для хранения вектора ускорения
VX_all = vcat(VX1[1:N_sw1],VX2[1:N_sw2],VX3[1:N_sw3]) # массив проекций скорости Х для всех трех участков
VY_all = vcat(VY1[1:N_sw1],VY2[1:N_sw2],VY3[1:N_sw3])# массив проекций скорости Y для всех трех участков
AX_all = vcat(AX1[1:N_sw1],AX2[1:N_sw2],AX3[1:N_sw3])# массив проекций ускорения Х для всех трех участков
AY_all = vcat(AY1[1:N_sw1],AY2[1:N_sw2],AY3[1:N_sw3])# массив проекций ускорения Y для всех трех участков
for i = 1:(length(X_all)-1)
V_all[i] = sqrt(VX_all[i]^2+VY_all[i]^2); # векторная сумма проекций скорости
A_all[i] = sqrt(AX_all[i]^2+AY_all[i]^2); # векторная сумма проекций ускорения
end
Построение графиков содержащих параметры движения
plot(X_all,Y_all,title="Координаты спускаемого аппарата")
plot(V_all,title="Суммарный вектор скорости спускаемого аппарата")
plot(A_all, title="Суммарный вектор ускорения спускаемого аппарата")
Перспективы дальнейшего использования данного скрипта
При необходимости дальнейшего исследования скрипт может быть доработан, в частности можно добавить параметры в дифференциальное уравнение такие как:
Учет изменения массы при работе тормозного двигателя
Влияние горизонтального и вертикального ветра на спускаемый аппарат
Угловые колебания космического аппарата (а моменты инерции по трем осям для которые также будут меняться при изменении массы космического аппарата)
Кроме того можно варьировать постоянные параметры и начальные условия таким образом чтобы максимально увеличить запасы устойчивости при движении спускаемого аппарата.
