Модель мягкой посадки многоразовой ракеты-носителя Falcon-9
Простейшая модель взаимодействия посадочной опоры многоразовой первой ступени РН Falcon-9 с поверхностью земли в процессе мягкой посадки.
Схема посадочной опоры
Рассмотрим упрощенную модель взаимодействия посадочной опоры многоразовой ступени ракеты-носителя с поверхностью земли. На рисунке 1 приведена фотография посадочных опор возвращаемой ступени РН Falcon-9. Ступень садится на 4 опоры. Каждая опора состоит из демпфирующего подкоса переменной длины BK и опорного подкоса АК постоянной длины, которые шарнирно закрепляются на корпусе РН.
Рисунок 1 - Посадочные опоры
Схема взаимодействия опоры с поверхностью земли приведена на рисунке 2. На рисунке показан корпус возвращаемого блока АВ, подкос постоянной длины \(AK = l\), демпфирующий подкос \(BK = s\). Со стороны поверхности земли на опору действуют сила реакции опоры \(\mathbf N\) и сила трения \(F\). После касания опорой поверхности земли демпфирующий подкос BK поглощает кинетическую энергию, уменьшаясь в длине.
Рисунок 2 - Силы, действующие на посадочную опору
Предположим, что в демпфирующем подкосе используется простейшее демпфирующее устройство – сотовый демпфер, который сминается при сжатии поглощая энергию. Такой способ поглощения кинетической энергии использовался в посадочных модулях Аполлон (рисунок 3). Предположим, что при смятии демпфера создается постоянная сила сопротивления \(D\).
Рисунок 3 - Сотовый демпфер
Математическая модель
Запишем уравнение вертикального поступательного движения возвращаемой ступени в процессе взаимодействия всех четырех посадочных опор с поверхностью земли, пологая, что все опоры “работают” одинаково, т.е. все опоры создают одинаковую силу сопротивления \(\mathbf D\) при сжатии демпфера и силу реакции опоры (земли) \(\mathbf N\). Уравнения движения запишем в форме уравнений Лагрнажа II-го рода, выбрав в качестве обобщенной координаты длину демпфирующего подкоса \(s = BK\).
Кинематические соотношения
Из теоремы косинусов для треугольника ABK найдём зависимость длины демпфирующего подкоса от высоты шарнира А над поверхностью земли h = AE, который определяет вертикальное положение ступени: \[s^2 = l^2 + a^2 - 2 a l \cos \alpha = l^2 + a^2 + 2 a l \cos (\pi - \alpha)\]
Косинус угла \(\cos (\pi - \alpha)\) определим из треугольника AEK \[\cos (\pi - \alpha) = \frac{h}{l}\]
Тогда: \[s^2 = l^2 + a^2 - 2 a l \cos \alpha = l^2 + a^2 + 2 a h\]
Продифференцировав это уравнение, получим зависимость скорости деформации подкоса \(\dot s\) от изменения вертикальной скорости ступени \(\dot h\): \[2 s \dot s = 2 a \dot h \quad \rightarrow \quad \dot s = \frac{a}{s} \dot h\]
Скорость деформации демпфера телескопического подкоса в \(s/a\) раз меньше вертикальной скорости посадочной ступени. Найденная зависимость \(s\) от \(h\) также позволяет записать уравнение в вариациях: \[s \cdot \delta s = a \cdot \delta h \quad \rightarrow \quad \delta s = \frac{a}{s} \delta h,\]
которое пригодится для определения обобщенных сил при выводе уравнения Лагранжа II-го рода.
Уравнения движения
Уравнение вертикального движения корпуса возвращаемой ступени при взаимодействии четырех опор с поверхностью земли будет иметь вид: \[\frac{d}{dt} \frac{\partial T}{\partial \dot s} - \frac{\partial T}{\partial s} = Q_s\]
Будем считать, что масса посадочных опор мала в сравнении с массой всей ступени. В этом случае кинетическая энергия \(T\) рассматриваемой механической системы будет определяться кинетической энергией корпуса: \[T = \frac{m \dot{h}^2}{2}\]
где \(m\) – масса возвращаемой ступени в момент посадки. Перепишем выражение для кинетической энергии, учитывая что \(\dot{h} = \dot h = \frac{s}{a} \dot s\): \[T = \frac{m \dot{s}^2}{2} \frac{s^2}{a^2}\]
Левая часть уравнения Лагранжа запишется в виде: \[\frac{d}{dt} \frac{\partial T}{\partial \dot s} - \frac{\partial T}{\partial s} = m \frac{d}{dt} \left( \frac{s^2}{a^2} \dot{s} \right) - m \frac{s}{a^2} \dot{s}^2 = \frac{ms}{a^2} \left( \dot s^2 + s \ddot s \right)\]
Обобщенную силу \(Q_s\) определим, записав выражение для элементарной работы сил на элементарном перемещение \(\delta s\): \[\delta A = \delta A_g + 4 \delta A_D + 4 \delta A_F\]
Коэффициент 4 учитывает то, что ступень садится на 4 посадочные опоры и при взаимодействии опор с поверхностью земли все четыре силы демпфера и все четыре силы трения совершают работу. Элементарная работа силы тяжести при вертикальном перемещении корпуса ступени определяется как \[\delta A_g = - m g \delta h = - m g \frac{s}{a} \delta s\]
Знак минус учитывает, что при положительной вариации \(\delta h\) (высота увеличивается) работа совершаемая силой тяжести будет отрицательной. Элементарная работа демпфера при отрицательном элементарном перемещении (сжатии демпфера) будет отрицательной, поэтому \[\delta A_D = D \delta s,\]
где \(D\) - сила сопротивления демпфера при его деформации. Элементарная работа силы трения при элементарном перемещении точки К вдоль поверхности земли \(\delta K\): \[\delta A_F = - F \delta K = - N \mu \frac{\partial }{\partial h} \left( \sqrt{l^2-h^2} \right) \delta h = N \mu \frac{hs}{a \sqrt{l^2-h^2}} \delta s\]
где \(\mu\) – коэффициент трения. Силу реакции опоры \(N\) определим, записав уравнение вертикального движения центра масс возвращаемой ступени под действием внешних сил: \[m \ddot{h} = 4N - m g\]
из которого следует, что \[N = (m \ddot{h} + m g)/4.\]
Вертикальное ускорение \(\ddot h\) определим из полученного ранее кинематического уравнения \[\dot s = \frac{a}{s} \dot h\]
продифференцировав его \[\ddot h = \frac{1}{a}(\ddot s s + \dot{s}^2)\]
Выражение для элементарной работы примет вид: \[\delta A = - m g \frac{s}{a} \delta s + 4 D \delta s + \mu m \frac{hs(\frac{1}{a}(\ddot s s + \dot{s}^2) + g)}{a \sqrt{l^2-h^2}} \delta s\]
Запишем обобщенную силу: \[Q_s = - m g \frac{s}{a} + 4 D + \mu m \frac{hs(\frac{1}{a}(\ddot s s + \dot{s}^2) + g)}{a \sqrt{l^2-h^2}}\]
Уравнение движения примет вид: \[m s \dot s^2 + m s^2 \ddot s = - m g s a + 4 D a^2 + \mu m \frac{hs(\ddot s s + \dot{s}^2 + a g)}{\sqrt{l^2-h^2}}\]
Разрешив это уравнение относительно старшей производной, получим: \[\ddot{s} = \frac{1}{s} \left( \frac{4 a^2 x}{m s \left(x-\mu h\right)}D-a g-\dot s^2 \right)\]
где \(x = \sqrt{l^2-h^2}\).
Стационарное решение этого уравнения при \(\ddot s = \dot s = 0\) \[\frac{4 a^2 x}{m s \left(x-\mu h\right)}D - a g = 0\]
Из стационарного решения можно вывести минимальное усилие, создаваемое демпфером, при котором система будет находится в равновесии: \[D_{min} = ag \frac{m s \left(x-\mu h\right)}{4 a^2 x},\]
которое в случае отсутствия трения опоры о поверхность земли (\(\mu = 0\)) примет вид: \[D_{min} = \frac{mg}{4} \frac{s}{a}\]
Энергия, поглощаемая демпфером
После касания опорами поверхности земли уменьшается кинетическая энергия ступени и её потенциальная энергия в поле силы тяжести за счет перемещения центра масс вниз на \(h_0 - h_k\).
Эта энергия поглощается демпферами посадочных опор: \[\frac{m V_0^2}{2} + m g (h_0-h_k) = 4 D(s_0-s_k)\]
В правой части этого уравнения энергия поглощаемая всеми четырьмя демпферами посадочных опор. Энергия поглощаемая одним демпфером: \[E_D = \frac{1}{4} \left[ \frac{m V_0^2}{2} + m g (h_0-h_k) \right]\]
Учитывая зависимость \(s\) от \(h\), получим \[\frac{m V_0^2}{2} + m g \frac{s_0^2-s_k^2}{2a} = 4 D(s_0-s_k)\]
Выполнив замену \(s_k = s_0 - \Delta s\), запишем квадратное уравнение для определения деформации демпфера \(\Delta s\) при заданной силе сопротивления \(D\): \[-\frac{mg}{2a} (\Delta s)^2 + (\frac{mg}{a}s_0 - 4D)\Delta s + \frac{m V_0^2}{2} = 0\]
и выражение для определения силы смятия демпфера \(D\) при заданной его деформации и начальной вертикальной скорости: \[D = \frac{m g}{4 a}\left[s_0 + \frac{V_0^2 a}{2 g \Delta s} - \frac{\Delta s}{2} \right]\]
На следующем рисунке показаны графики минимального усилия смятия демпфера в зависимости от его максимальной деформации и начальной вертикальной скорости возвращаемого блока.
Энергия поглощаемая демпфером одной посадочной опоры: \[E_D = \frac{m g}{4 a}\left[s_0 + \frac{V_0^2 a}{2 g \Delta s} - \frac{\Delta s}{2} \right] \Delta s\]
Численное интегрирование уравнения движения
Численно проинтегрируем полученное уравнение движения в среде python.
# Именованный кортеж для хранения параметров системы
from collections import namedtuple
# Массивы и матрицы
import numpy as np
# Численное интегрирование дифференциальных уравнений
from scipy.integrate import solve_ivp
# Графики
import matplotlib.pyplot as plt
Функция правых частей системы дифференциальных уравнений:
def dydt(t, y, p):
# Длина демпфера
s = y[0]
# Скорость изменения длины демпфера ds/dt
ds = y[1]
# Высота
h = np.sqrt((s**2 - p.l**2 - p.a**2)/(2*p.a))
# x
x = np.sqrt(p.l**2-h**2)
D = p.D
# dv/dt =
d2s = 4*p.a*p.a*D*x/(p.mass*(x-p.mu*h)*s*s) - (p.a*p.g+ds*ds)/s
return np.array([ds,d2s])
Функция для остановки процесса интегрирования, когда скорость изменения деформации демпфера достигнет нуля:
def event_ds_eq_0(t, y):
# Функция-"детектор", передаваемая в интегратор (параметр events),
# для определения времени достижения нулевой скорости вертикальной скорости
# остановки процесса интегрирования
return y[1]
# функция определяется условие h = 0 при движении "вниз"
event_ds_eq_0.direction = 1
# функция-детектор активна
event_ds_eq_0.terminal = True
Создаем именованный кортеж params с параметрами системы
params = namedtuple("params", "l a mass g D")
# Масса ступени
params.mass = 30000.0
# Расстоянием между шарнирами подкосов
params.a = 3.0
# Длина подкоса постоянной длины
params.l = 6.5
# Ускорение свободного падения
params.g = 9.807
# Коэффициент трения
params.mu = 0.0
# Сила, создаваемая демпфером
params.D = 250000
# Начальные условия:
# начальная вертикальная скорость
dh0 = -2.0
# начальная длина демпфирующего подкоса
s0 = 8.5
# начальная скорость деформации демпфера
ds0 = dh0*params.a/s0
Запуск процесса численного интегрирования:
# Без учета трения
sol = solve_ivp(lambda t,y: dydt(t,y,params), [0, 1], [s0, ds0], method='LSODA', events = [event_ds_eq_0], rtol = 1e-9)
# С учетом трения
params.mu = 0.4
sol_mu01 = solve_ivp(lambda t,y: dydt(t,y,params), [0, 1], [s0, ds0], method='LSODA', events = [event_ds_eq_0], rtol = 1e-9)
Для сравнения получены решения для нулевого трения и для трения \(\mu = 0.4\). Используя полученные решения, вычисляем (таблицы) высоту и вертикальную скорость:
# Определяем таблицу изменения высоты
h = (sol.y[0]**2-params.l**2-params.a**2)*0.5/params.a
h_mu01 = (sol_mu01.y[0]**2-params.l**2-params.a**2)*0.5/params.a
# и вертикальной скорости
dh = sol.y[0]*sol.y[1]/params.a
dh_mu01 = sol_mu01.y[0]*sol_mu01.y[1]/params.a
# и время остановки (достижение нулевой скорости)
tmax = sol.t.max()
Построение графиков изменения длины подкоса BК:
plt.figure(figsize=[11,5])
plt.subplot(121)
plt.plot(sol.t,sol.y[0]-s0);
plt.plot(sol_mu01.t,sol_mu01.y[0]-s0);
plt.xlabel('t, c'); plt.ylabel('$\Delta s$, м');plt.xlim([0,tmax])
plt.grid(True,ls=':')
plt.legend(['$\mu=0.0$','$\mu=0.4$'])
plt.title('Деформация демпфера')
plt.subplot(122)
plt.plot(sol.t,sol.y[1]);
plt.plot(sol_mu01.t,sol_mu01.y[1]);
plt.xlabel('t, c'); plt.ylabel('ds/dt, м/c'); plt.xlim([0,tmax])
plt.grid(True,ls=':')
plt.legend(['$\mu=0.0$','$\mu=0.4$'])
plt.title('Скорость деформации демпфера')
plt.tight_layout()
plt.savefig('Falcon_s_ds.png',dpi=160)
При выбранных параметрах сотовый демпфер поглощает кинетическую энергию на перемещении (деформации) около 0,33 м за 0,9 секунд. Можно оценить поглощенную энергию или суммарную по модулю работу демпферов : \(А = 4 D \cdot \Delta s = 4 \cdot 250 кН \cdot 0.33 м = 83 кДж\).
Построение графиков изменения высоты:
plt.figure(figsize=[14,4])
plt.subplot(121)
plt.plot(sol.t,h);
plt.plot(sol_mu01.t,h_mu01);
plt.xlabel('t, c'); plt.ylabel('h, м');plt.xlim([0,tmax])
plt.grid(ls=':')
plt.legend(['$\mu=0.0$','$\mu=0.4$'])
plt.title('Высота')
plt.subplot(122)
plt.plot(sol.t,dh);
plt.plot(sol_mu01.t,dh_mu01);
plt.xlabel('t, c'); plt.ylabel('dh/dt, м/c');plt.xlim([0,tmax])
plt.grid(ls=':')
plt.legend(['$\mu=0.0$','$\mu=0.4$'])
plt.title('Вертикальная скорость')
plt.savefig('Falcon_h_dh.png',dpi=150)
Задание
Масса возвращаемой ступени в момент посадки равна 30 т, AK=6.5 м, BK = 8 м, \(a\) = 3 м. Ступень касается земли с начальной вертикальной скоростью 3 м/с.
- Найти минимальное необходимое усилие смятия демпфера \(D_{0.4}\) при \(\mu = 0.4\) и деформации демпфера не более 0.5 м.
- Найти минимальное необходимое усилие смятия демпфера \(D_{0.0}\) при \(\mu = 0.0\) и деформации демпфера не более 0.5 м.
- Верифицировать модель, убедившись в равенстве нулю суммы начальной кинетической энергии возвращаемой ступени, работы сили тяжести, работы демпферов и сил трения опор.
