Посадка на Луну
Моделирование процесса мягкой посадки на Луну космического аппарата в среде Python
Алгоритм посадки
Предположим, что КА находится на опорной круговой орбите вокруг Луны высотой \(h_0\). Посадка КА выполняется в четыре этапа:
- двухимпульсный переход с опорной круговой орбиты на круговую посадочную орбиту или орбиту ожидания с которой будет выполнятся спуск;
- спуск (гравитационный разворот) с посадочной орбиты до некоторой небольшой высоты начала вертикального участка мягкой посадки (~100 м) – третье включение двигателя;
- вертикальное движение КА c выключенным двигателем (свободное падение);
- гашение вертикальной скорости и касание поверхности (четвёртое включение двигателя).
На первом этапе выполняется переход с начальной опорной круговой орбиты \(h_0\) на посадочную круговую орбиту \(h_1\), с которой будет выполнятся переход на траекторию спуска. Переход на посадочную орбиту выполняется по гомановской траектории.
Второй этап (от \(t_1\) до \(t_2\) ) – это этап спуска с посадочной круговой орбиты до достижения заданной высоты \(h_2\) и скорости \(v_2\). На этом этапе двигатель КА создает постоянную тягу \(F\), направленную против вектора скорости. При такой схеме торможения обеспечивается расход топлива близкий к минимальному. Эта схема посадки использовалось КА Сервейер.
Двигатель выключается на заданной высоте \(h_2\) при достижении скорости снижения значения \(v_2\) (например, \(h_2 = 100\) м и \(v_2 = -1\) м/с соответственно). Для выполнения этого условия должна быть определена высота круговой орбиты \(h_1\), с которой начинается спуск. Это высота может быть определена методом “пристрелки” путем интегрирования дифференциальных уравнений движения КА на втором этапе при различных значениях начальной высоты \(h_1\) до достижения условия \[h(t_2) = h_2, \quad v(t_2) = v_2\]
После достижения высоты \(h_2\) двигатель КА выключается. Далее в течение некоторого времени (третий этап) КА свободно падает до высоты запуска двигателя \(h_3\). Работа двигателя на четвертом участке обеспечивает касание поверхности Луны с нулевой вертикальной скоростью.
Время запуска двигателя для мягкой посадки может быть определено из уравнений вертикального движения КА в однородном гравитационном поле как материальной точки переменной массы \[v_2 - g(t_{3e}+t_{23}) - I_s \ln \frac{m_2-bt_{3e}}{m_2} = 0\] \[h_2+v_2(t_{3e}+t_{23})-\frac{g}{2}(t_{3e}+t_{23})^2-I_s\left\{\left[t_{3e}-\frac{m_2}{b}\right]\left[\ln(1-\frac{b}{m_2}t_{3e})-1\right]-\frac{m_2}{b}\right\}=0\]
где \(t_{23} = t_3 - t_2\) – продолжительность пассивного вертикального участка, \(t_{3e} = t_e - t_3\) – продолжительность работы двигателя.
Уравнения движения
Уравнение движения посадочного модуля, как материальной точки переменной массы в центральном гравитационном поле Луны имеет вид: \[\ddot{\vec{r}} = - \mu \frac{\vec r}{r^3} - \vec{e}_v \frac{F}{m}, \quad \dot m = - \frac{F}{I_s}\]
где \(\mu\) – гравитационный параметр Луны, \(F\) – сила тяги двигателя, \(I_s\) – удельный импульс, \(\vec r\) – радиус-вектор посадочного модуля относительно центра Луны, \(\vec{e}_v\) вектор направления скорости посадочного модуля: \[\vec{e}_v = \frac{\dot{\vec r}}{|\dot{\vec r}|}.\]
Проинтегрируем эти уравнения численно в среде Python.
Модель посадки в среде Python
Подключим необходимые библиотеки
# Численные методы
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
# Класс для хранения данных
from dataclasses import dataclass
from typing import Callable
# Графики
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.pylab as pylab
Параметры рассматриваемой системы будем хранить в объекте типа dataclass.
@dataclass
class Data:
# Гравитационный параметр Луны
mu : float = 4890.0e9;
# Радиус Луны
Rm : float = 0.5*3476e3;
# Начальная масса
m : float = 1500
# Зависимость тяги двигателя
# посадочного модуля от времени
F : Callable[[float], float] = lambda t: 0
# Максимальная тяга
Fmax : float = 10e3
# Удельный импульс
Isp : float = 3500
Создаем экземпляр класса Data (объект)
p = Data()
Объявим функцию правых частей дифференциальных уравнений движения посадочного модуля. Первый аргумент функции – время, второй – вектор состояния системы на момент времени \(t\), третий аргумент – структура с данными системы. Вектор состояния системы это массив из 7 элементов: \[q = [x, y, z, v_x, v_y, v_z, m_f]\]
где \(x, y, z\) – координаты центра масс модуля в луноцентрической системе координат OXYZ, \(v_x, v_y, v_z\) – проекции скорости посадочного модуля в системе координат OXYZ, \(m_f\) – масса израсходованного топлива.
def dqdt(t, q, p):
r = q[0:3]
v = q[3:6]
# Израсходовано топлива
mf = q[6]
# Модуль скорости
vnorm = np.linalg.norm(v)
# Направление вектора скорости
ev = v/vnorm
# Текущая масса модуля
mass = p.m - mf
# Ускорение
a = - p.mu*r/(np.linalg.norm(r)**3) - ev*p.F(t)/mass
# Секундный расход
dm = np.abs(p.F(t))/p.Isp
return np.hstack([v,a,dm])
Объявим функцию, контролирующую достижение заданной радиальной скорости \[v_r = \frac{\vec r}{|\vec r|} \cdot \dot{\vec r}\]
Эта функция будет использоваться для остановки процесса интегрирования при достижении заданной радиальной (вертикальной) скорости vk. Аргументы функции: время t, вектор состояния на момент t, параметры системы p и контролируемое значение скорости vk:
def event_hv(t, q, p, vk):
r = np.linalg.norm(q[0:3])
h = r - p.Rm
er = q[0:3]/r
vr = np.dot(q[3:6],er);
return vr-vk
Функция будет возвращать ноль при равенстве вертикальной скорости заданному значению.
Объявим функцию, которая определяет высоту, при которой достигается заданная вертикальная скорость в конце второго этапа движения (торможения) \(v_2\) при движении с начальной круговой орбиты высотой \(h_1\). В функции производится интегрирование уравнений движения посадочного модуля до достижения условия \(v_r = v_2\):
def get_he_v_eq_0(h1, p, v2):
# Начальные условия
R0 = p.Rm+h1;
V0 = np.sqrt(p.mu/R0)
# Вектор начальных условий
q0 = np.array([R0,0,0, 0,V0,0, 0.0])
# Уточняем функцию "событий"
ev = lambda t,q: event_hv(t,q,p,v2)
# Останавливаем процесс интегрирования
ev.terminal = True
# если функция event_hv (ее возвращаемое значение --
# разность скоростей) пересекает ноль при возрастании
ev.direction = 1
# Запускаем процесс численного интегрирования
sol = solve_ivp(lambda t,q: dqdt(t,q,p), [0,500], q0, method='RK45', events = [ev,], rtol = 1e-8)
# Определяем и возвращаем конечную высоту:
# высота на которой скорость стала равна vk
he = np.sqrt(np.sum(sol.y[0:3]**2,axis=0))-p.Rm
return he[-1]
Начальные условия интегрирования соответствуют движению КА по круговой орбите: \[\mathbf r(0) = [ R_m + h_1,0,0]^T, \quad \mathbf v(0) = [0,\sqrt{\mu/(R_m+h_1)},0]^T\]
Задавая различные значения начальной высоты h1, мы можем получить значения высоты, при которой достигается заданная вертикальная скорость v2.
“Включаем” постоянную тягу
p.F = lambda t: p.Fmax
При помощи функции scipy.optimize.fsolve находим высоту посадочной круговой орбиты: включив двигатель на этой высоте и поддерживая направление силы тяги против вектора скорости, посадочный модуль достигнет высоты \(h_2\) = 100 м с вертикальной скоростью \(v_2 = -1\) м/с.
h2 = +100.0
v2 = -1.0
h1 = fsolve(lambda x: get_he_v_eq_0(x[0], p, v2)-h2, [10000.0])[0]
Функция scipy.optimize.fsolve численным методом находит решение нелинейного уравнения \[f(\vec x) = 0\]
В рассматриваемом случае искомым аргументом этой функции является высота посадочной орбиты \(h_1\). Получим численное решение для найденной начальной высоты \(h_1\):
R0 = p.Rm+h1
V0 = np.sqrt(p.mu/R0)
q0 = np.array([R0,0,0, 0,V0,0.0, 0.0])
ev = lambda t,q: event_hv(t,q,p,-1.0)
ev.direction = 1
ev.terminal = True
sol = solve_ivp(lambda t,q: dqdt(t,q,p), [0,500], q0, method='RK45', events = [ev,], rtol = 1e-8)
Построим графики изменения высоты и вертикальной скорости посадочного модуля
# Модуль радиус-вектора
r = np.sqrt(np.sum(sol.y[0:3]**2,0))
# Высота
height = r-p.Rm
# Модуль скорости
velocity = np.sqrt(np.sum(sol.y[3:6]**2,0))
# Единичный вектор местной вертикали
er = sol.y[0:3].transpose()/np.repeat(np.sqrt(np.sum(sol.y[0:3].transpose()**2,axis=1)),3).reshape(-1,3)
# Вертикальная скорость
velocity_r = np.sum(er*sol.y[3:6].transpose(),axis=1)
plt.figure(figsize=[14,4])
plt.subplot(1,2,1)
plt.plot(sol.t,height)
plt.xlabel('t, c')
plt.ylabel('h, м')
plt.grid(ls=':')
plt.title('Высота');
plt.subplot(1,2,2)
plt.plot(sol.t,velocity_r)
plt.xlabel('t, c')
plt.ylabel('$v_r$, м/c')
plt.grid(ls=':')
plt.title('Вертикальная скорость');
Определим высоту начала пассивного участка, зададим скорость посадки (0) и вычислим секундны расход топлива:
# Высота начала вертикального участка
h2 = (np.sqrt(np.sum(sol.y[0:3]**2,0))-p.Rm)[-1];
# Скорость посадки
Ve = 0;
# Секундный расход топлива
b = abs(p.Fmax)/p.Isp;
# Ускорение свободного падения (постоянное)
g = p.mu/p.Rm**2;
Зададим массу посадочного модуля на момент \(t_2\)
m2 = p.m-sol.y[6][-1]
Решим систему уравнений для определения продолжительности свободного падения (\(t_3 - t_2\)) и продолжительности работы двигателя (\(t_e - t_3\))
f1 = lambda tp,ta: v2-g*(ta+tp)-p.Isp*np.log((m2-b*ta)/m2)-Ve
f2 = lambda tp,ta: h2+v2*(ta+tp)-g/2*(ta+tp)**2-p.Isp*((ta-m2/b)*(np.log(1-b/m2*ta)-1)-m2/b);
f = lambda x: (f1(x[0],x[1]), f2(x[0],x[1]));
t23e = fsolve(f,[1,1])
print('Начальная высота {:3.1f} м'.format(h2))
print('Продолжительность свободного падения {:3.1f} c'.format(t23e[0]))
print('Продолжительность торможения {:3.1f} c'.format(t23e[1]))
Начальная высота 100.0 м
Продолжительность свободного падения 9.7 c
Продолжительность торможения 1.8 c
С начальными условиями, соответствующими моменту \(t_2\)
t0 = sol.t[-1];
q0 = sol.y[:,-1];
Проинтегрируем уравнения движение посадочного модуля на интервале \(t_e - t_2\) с учетом того, что двигатель запускается через \(t_3 - t_2\) после начала движения
# Конечное время
te = t23e[0]+t23e[1];
# Закон изменения тяги
p.F = lambda t: (t>t0+t23e[0])*p.Fmax
# Интегрируем...
sol2 = solve_ivp(lambda t,q: dqdt(t,q,p), [t0,t0+te], q0, method='BDF', rtol = 1e-8)
Построим графики изменения высоты и скорости посадочного модуля на третьем и четвертом этапах движения:
# Модуль радиус-вектора
r = np.sqrt(np.sum(sol2.y[0:3]**2,0))
# Высота
h = r-p.Rm
# Модуль скорости
v = np.sqrt(np.sum(sol2.y[3:6]**2,0))
# Единичный вектор местной вертикали
er = sol2.y[0:3].transpose()/np.repeat(r,3).reshape(-1,3)
# Радиальная скорость
vr = np.sum(er*sol2.y[3:6].transpose(),axis=1)
plt.figure(figsize=[14,5])
plt.subplot(1,2,1)
plt.plot(sol2.t,h,lw=3)
plt.xlabel('t, c')
plt.ylabel('h, м')
plt.grid(True, lw=0.5, zorder=0,ls=':')
plt.title('Высота');
plt.subplot(1,2,2)
plt.plot(sol2.t,vr,lw=3)
plt.xlabel('t, с')
plt.ylabel('$v_r$, м/с')
plt.grid(True, lw=0.5, zorder=0,ls=':')
plt.title('Вертикальная скорость');
В свободном падении с высоты 100 м с начальной скоростью минус 1 м/с посадочный модуль до запуска двигателя достигает скорости снижения более 15 м/с, а затем в течение 1,8 секунд гасит вертикальную скорость до момента касания поверхности Луны. Суммарный расход топлива без учета затрат на гомановский переход с опорной орбиты на посадочную орбиту составляет около 584 кг:
sol2.y[6,-1]
583.89