Посадка на Луну

Моделирование процесса мягкой посадки на Луну космического аппарата в среде Python

Алгоритм посадки

Предположим, что КА находится на опорной круговой орбите вокруг Луны высотой \(h_0\). Посадка КА выполняется в четыре этапа:

  1. двухимпульсный переход с опорной круговой орбиты на круговую посадочную орбиту или орбиту ожидания с которой будет выполнятся спуск;
  2. спуск (гравитационный разворот) с посадочной орбиты до некоторой небольшой высоты начала вертикального участка мягкой посадки (~100 м) – третье включение двигателя;
  3. вертикальное движение КА c выключенным двигателем (свободное падение);
  4. гашение вертикальной скорости и касание поверхности (четвёртое включение двигателя).

На первом этапе выполняется переход с начальной опорной круговой орбиты \(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

© 2023. All rights reserved.

Powered by Hydejack v9.1.6