Моделирование посадки космического аппарата на Венеру
goltsevРасчет параметров траектории спускаемого аппарата при посадке на Венеру
Введение
В данной статье приведен один из возможных вариантов математического обеспечения для исследования динамики неуправляемого движения спускаемых аппаратов, созданное на основе обобщения опыта проектирования, отработки и эксплуатации аппаратов типа "Венера" и "Марс" второго и третьего поколений. Приводятся обобщенные данные, характеризующие динамику спуска в атмосферах планет типа "Венера 9-16", "Вега-1,2" и "Марс-3,6". При данном моделирования рассматривается процесс посадки на Венеру и заданы параметры её атмосферы.
Венера является одной из самых изучаемых планет Солнечной системы. И не удивительно, поскольку Венера- это ближайшая к нам планета Солнечной системы. На протяжении второй половины 20века и начала 21 века были запущены на Венеру очень много космических аппаратов(начиная с Венеры 1). Тем более, что 8 июня 1975 года космическим аппаратом Венера 9 были получены первые черно-белые снимки поверхности Венеры. От того ,как пройдет посадка спускаемого аппарата зависит дальнейшее исследование планеты и получение информации о ней. Посадка бывает активной и пассивной. Активная посадка отличается от пассивной тем, что в пассивной посадке не используются посадочные механизмы, такие как парашют или тормозной механизм. В активной посадке используется парашют, чтобы снизить скорость аппарата.
В настоящее время используется два основных типа траекторий интенсивного торможения: баллистическое уменьшение скорости за счет силы лобового сопротивления аппарата и планирующее - с использованием подъёмной силы аппарата для увеличения времени торможения и, следовательно, для снижения конечной скорости, уменьшения максимальных перегрузок и тепловых потоков. Для организации последующего торможения на втором участке могут быть использованы также две схемы: активное торможение с помощью реактивных двигателей или торможение с помощью парашютных либо других аэродинамических систем. Возможно применение и комбинированной схемы торможения.
Подключаемые библиотеки:
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="Плотность атмосферы (интерполяция)")
Моделирование траектории движения спускаемого аппарата (кусочно)
В рамках данной статьи, предполагается что траектория спускаемого аппарата состоит из трех участков:
-
Пассивное торможение. Спускаемый аппарат входит в атмосферу планеты и под воздействием динамического трения у него снижается скорость движения. Данный режим длится до тех пор пока скорость не позволит раскрыть парашют.
-
Пассивное торможение с раскрытым парашютом. Наиболее длинный участок траектории, спускаемый аппарат плавно теряет скорость и высоту.
-
Активное торможение. Завершающий участок траектории на котором происходит включение тормозного двигателя и торможение.
Для моделирования движения будет использоваться система дифференциальных уравнений
Дифференциальное уравнение описывающее движение
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 # Vx
dxdt[2] = -9.87 - k/m*x[2]+F_torm/m # Vy
dxdt[3] = x[1] # x
dxdt[4] = x[2] # y
return dxdt
end
Функция в которой реализован метод Эйлера
function euler_step(f, t, x, dt, S,A2,F_torm)
dx = f(t, x, S, A2,F_torm)
return x + dt * dx
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)
Расчет точек второго участка траектории
Vx0_1 = VX1[N_sw1]; # начальная проекция скорости по Х
Vy0_1 = VY1[N_sw1]; # начальная проекция скорости по Х
X0_1 = X1[N_sw1];
Y0_1 = Y1[N_sw1];
# начальное состояние
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)
Расчет точек третьего участка траектории
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
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)
После того как построены участки траектории, соберем их в один массив
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])
L_all = zeros(length(X_all)-1,1)
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])
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])
for i = 1:(length(X_all)-1)
L_all[i] = sqrt(X_all[i]^2+Y_all[i]^2);
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="Суммарный вектор ускорения спускаемого аппарата")