Модель движения летающего стенда возвращаемой ступени РН в SimInTech
Модель иллюстрирует движение летающего стенда для отработки мягкой посадки первой ступени многоразовой ракеты-носителя. При помощи подобного летающего стенда компания SpaceX начиная с 2012 года отрабатывала технологию посадки возвращаемой ступени РН Falcon-9.
Схема системы
Схема системы показана на рисунке 1. Движение стенда (далее ступень) рассматривается относительно инерциальной системы координат \(O X_o Y_o\), связанной с землей. Положение ступени определяется координатами центра масс \(x\) (высота), \(y\) и углом наклона ступени \(\theta\) к вертикали.
Рисунок 1 - Схема системы
Движение ступени происходит под действием силы тяжести \(G\), приложенной в центре масс, силы тяги двигателя \(P\), приложенной в точке \(p\), и аэродинамической силы \(F_a\), которая приложена в точке \(a\) (центр давления). Аэродинамическая сила представлена в виде суммы продольной \(X_a\) и поперечной \(Y_a\) составляющих. Направление силы тяги относительно продольной оси ступени определяется углом поворота сопла двигателя \(\gamma_p\).
Уравнения движения
Уравнения движения ступени построены при следующих допущениях:
- Рассматривается действие на ступень только поперечной аэродинамической силы \(Y_a\), которая будет определяться горизонтальной скоростью ступени и скоростью ветра.
- Положение точки приложения аэродинамической силы не изменяется относительно корпуса ступени при развороте ступени: \(ba = x_a\) const. Также не изменяется в процессе движения момент инерции ступени и положение её центра масс относительно корпуса (\(bC = = x_c\) const).
- Полет происходит на небольшой высоте, гравитационное ускорение остается постоянным.
- Не учитывается действие на ступень газодинамических сил от струй двигателя, отражающихся от поверхности земли, при движении на небольшой высоте.
- Динамические свойства системы управления тягой двигателя и вектором тяги не учитываются.
Вертикальное движение ступени определяется тягой двигателя и силой притяжения Земли. Уравнение вертикального движения имеет вид: \[m \ddot{x} = P \cos(\theta + \gamma_p) - G + Y_a \sin \theta\]
где \(G = m g\) – вес ступени, \(g\) – ускорение свободного падения, \(m\) – масса ступени, \(P\) – тяга двигателя, \(Y_a\) – аэродинамическая сила: \[Y_a = S_m C_y \frac{\rho (\dot y - w)^2}{2} \text{sgn}(\dot y - w)\]
где \(S_m\) – характерная площадь ступени (площадь Миделя), к которой отнесен аэродинамический коэффициент поперечной силы \(C_y\), \(w\) – скорость ветра, \(\dot y\) – горизонтальная скорость ступени, \(\rho\) – плотность воздуха.
Горизонтальное движение ступени описывается уравнением: \[m \ddot{y} = P \sin(\theta + \gamma_p) - Y_a \cos \theta\]
Уравнение движения ступени вокруг центра масс: \[J_z \ddot{\theta} = P (x_c-x_p) \sin \gamma_p - Y_a (x_a-x_c)\]
где \(x_a\) – координата x точки приложения аэродинамической силы в базовой системе координат \(b x_b y_b\) ступени, \(x_p\) – координата x точки приложения силы тяги ступени в базовой системе координат, \(x_c\) – координата x центра масс ступени в базовой системе координат, \(J_z\) – поперечный момент инерции ступени.
Уравнения движения ступени интегрируются совместно с дифференциальным уравнением изменения массы при расходе топлива: \[\frac{dm}{dt} = - \frac{P}{I_s}\]
где \(I_s\) – удельный импульс двигателя ступени [м/c].
Модель SimInTech
Структура модели
Модель SimInTech состоит из четырех блоков (субмоделей):
- блок динамической модели ступени – “Попрыгунчик”;
- блок управления тягой двигателя;
- блок управления углом поворота сопла;
- блок, задающий программное изменение высоты;
- блок, задающий программное изменение бокового смещения ступени.
Рисунок 2 - Общий вид модели
Субмодель динамической модели ступени
Динамическая модель ступени создана на основе стандартного блока “Субмодель” с дополнительными свойствами (момент инерции, аэродинамический коэффициент, …). В созданной субмодели используется блок “Язык программирования”, в котором записаны дифференциальные уравнения движения ступени, используя встроенный язык программирования SimInTech. Интегрирование уравнений производится при помощи стандартных блоков типа “Интегратор”. На вход субмодели подаются сила тяги двигателя \(P\), угол поворота сопла \(\gamma_p\) и скорость горизонтального ветра \(W\). Выходы субмодели: кинематические параметры ступени (положение и скорость центра масс, угол поворота ступени, её угловя скорость).
Рисунок 3 - Структура субмодели “Попрыгунчик”
Для построения графиков и упрощения экспорта результатов моделирования в текстовый файл (таблица значений кинематических параметров) в проекте используются сигналы.
Рисунок 4 - Сигналы проекта
Значения в базу сигналов записываются с линий связи субмодели “Попрыгунчик”. Сигналы связываются со значениями с линий связи субмодели при помощи редактора связей
Рисунок 5 - Редактор связей
Код блока “Язык программирования” – дифференциальные уравнения:
// Плоская модель движения возвращаемой ступени
// для отработки посадки
// "Попрыгунчик"
// m - масса
// P - тяга
// w - горизонтальная скорость ветра
// gamma_p - угол поворота сопла
// dx, dy - вертикальная и горизонтальная скорость
// theta - угол наклона ступени
// dtheta - угловая скорость ступени
// rho - плотность воздуха
input m, P, w, gamma_P, dx, dy, theta, dtheta, rho;
output d2x, d2y, d2theta, Fa;
// Аэродинамическая сила
Fa = 0.5*rho*(w-dy)^2*Cy*Sm*sign(w-dy);
// Плечо аэродинамической силы
ha = xa-xc;
// Плечо силы тяги
hp = xc-xp;
// Уравнения
d2x = ( p*cos(theta+gamma_p)-m*g + Fa*sin(theta) )/m;
d2y = ( p*sin(theta+gamma_p)-Fa*cos(theta) )/m;
d2theta = (-p*sin(gamma_p)*hp-Fa*ha )/Jz;
Параметр \(g\) – ускорение свободного падения, объявлена в блоке инициализации модели. Другие параметры, используемые в скрипте блока (Sm, Cy, xc, xp, xa) определены в свойствах субмодели “Попрыгунчик”:
Рисунок 6 - Параметры субмодели “Попрыгунчик”
Для того, чтобы добавить эти свойства к субмодели необходимо перевести SimInTech в режим разработчика (Меню “Вид” - “Режим разработчика”), в общих свойствах субмодели указать новое наименование типа блока, теперь это не “стандартная” субмодель, а блок типа “TJumpStage”. Процедура создания блоков на основе субмоделей описана в справочной системе SimInTech.
Блок управления тягой
Учитывая, что ступень движется с небольшой вертикальной скоростью, тяга двигателя будет близка к весу ступени \(G\). Используемый в блоке управления тяги ПД-регулятор вычисляет “добавку” к этому весу. В блоке управления тягой сравнивается фактическая высота ступени \(x\) с программным значением высоты \(x_p\) и на основе этой разницы (ошибки) формируется значение тяги: \[P = G + k_1 e_x + k_2 \dot e_x\]
где \(e_x = x_p - x\).
Рисунок 7 - Структура блока управления тягой двигателя
Блок управления углом поворота сопла
Горизонтальное движение ступени производится за счет наклона её корпуса. Для наклона корпуса используется момент относительно центра масс, создаваемый тягой двигателя при ненулевом угле поворота сопла \(\gamma_p\). В блоке управления углом поворота сопла при возникновении разницы между программным значением горизонтальной координаты ступени \(y_p\) и фактическим значением \(y\) определяется угол наклона ступени при помощи ПД-регулятора: \[\theta_{PD} = k_3 e_y - k_4 \dot e_y, \quad e_y = y_p-y\]
При этом максимальный угол поворота ступени ограничивается некоторым предельным значением \(\theta_{max}\): \[\theta_p = \begin{cases} \theta_{PD} & |\theta_{PD}| \leq \theta_{max} \\ \theta_{max} & |\theta_{PD}| > \theta_{max} \end{cases}\]
Полученное значение угла поворота ступени используется как программное значение \(\theta_p\), которое на следующем этапе сравнивается с текущим значением угла \(\theta\). Разница между программным и фактическим значением поступает на вход ПИД-регулятора, который формирует сигнал угла поворота сопла двигателя. Этот угол также ограничен некоторым предельным значением: \[\gamma_{PID} = k_5 e_{\theta} + k_6 \dot e_{\theta} + k_7 \int_0^t \dot e_{\theta} d \tau, \quad e_{\theta} = \theta_p-\theta\] \[\gamma_p = \begin{cases} \gamma_{PID} & |\gamma_{PID}| \leq \gamma_{max} \\ \gamma_{max} & |\gamma_{PID}| > \gamma_{max} \end{cases}\]
На вход блока управления углом поворота сопла также поступает горизонтальная скорость блока. При превышении горизонтальной скорости некоторого предельного значения программное значение угла наклона ступени становится равным нулю.
Рисунок 8 - Структура блока управления углом поворота сопла
Блок, задающий программное изменение высоты
Программа управления высотой имеет вид: \[x_p = \frac{h_{max}}{2} \left[ 1 + \sin \left( \frac{2\pi}{t_k}t+\frac{3 \pi}{2} \right) \right]\]
где \(t_k\) – конечное время, \(h_{max}\) – максимальная высота.
Рисунок 9 - Структура блока, задающего программное значение высоты
Блок, задающий программное изменение бокового смещения
Программа управления боковым перемещение ступени определяется выражением: \[y_p = \begin{cases} 0, & t<t_1 \\ \frac{y_k}{2} \left[ \sin( \frac{t-t_1}{t_2-t_1}\pi-\frac{\pi}{2})+1 \right], & t_1 \leq t \leq t_2 \\ 0, & t>t_2 \end{cases}\]
где \(t_1 = \Delta t_v\) – продолжительность вертикального участка при взлете ступени, \(t_2\) - t_{k} - \Delta t_v \(начало вертикального участка спуска ступени,\) t_k \(-- конечное время (продолжительность полета),\) y_k $$ – расстояние от точки взлета до точки посадки.
Это программа управления задается в блоке “Язык программирования”:
// Горизонтальное перемещение
output y_p;
// Начало участка вертикального спуска
t2 = endtime-t_vert;
// Горизонтальное перемещение
yend = y_end;
// Вертикальный участок
if time<=t_vert then y_p = 0;
// Участок горизонтального перемещения
if time>t_vert and time<=t2 then y_p = (sin((time-t_vert)*pi/(t2-t_vert)-pi/2)*0.5+0.5)*yend;
if time>t2 then y_p = yend;
Скрипт модели
В скрипте модели объявляются переменные, используемые в коде и свойствах блоков.
initialization
// Максимальный угол поворота сопла
var gamma_max = 4*pi/180;
// Максимальный угол наклона ступени
var theta_max = 4*pi/180;
// Продолжительность вертикального участка
var t_vert = 10;
// Максимальная горизонтальная скорость ступени
var v_hor_max = 10;
// Горизонтальное перемещение ступени
var y_end = 150;
// Максимальная высота подъема ступени
var h_max = 50;
end
Графики
Для того, чтобы не загромождать модель блоками типа “Временной график” для построения графиков используется менеджер данных:
Рисунок 10 - Менеджер данных
ЗD Визуализация
Для визуализации движения блока средствами SimInTech используются блоки из раздела Визуализация 3D – блоки типа “Просмотрщик” и “Объект”. Один из блоков “Объект” отображает ступень, второй – поверхность неподвижной земли.
Рисунок 11 - Визуализация 3D - Просмотрщик
Свойства блока “Объект”, используемого для отображения ступени показаны на рисунке 12. Положение объекта задается выражениями, записанными при помощи сигналов x, y и theta.
Рисунок 12 - Свойства блока “Объект” для отображения ступени
Результаты
Далее приведены результаты моделирования движения ступени. В соответствии с программой полета ступень поднимается на максимальную высоту 50 метров. Посадка ступени производится на расстоянии 150 метров от места взлета. На рисунке 12 показаны графики изменения высоты и горизонтального перемещения ступени.
Рисунок 13 - Изменение высоты и горизонтального перемещения ступени
На рисунке 13 показаны графики изменения вертикальной и горизонтальной скорости ступени.
Рисунок 14 - Изменение вертикальной и горизонтальной скорости ступени
На рисунке 15 показаны графики изменения угла поворота ступени и угла поворота сопла двигателя. На рисунке можно увидеть реакцию ступени на воздействие горизонтального ветра.
Рисунок 15 - Угол наклона ступени и угол поворота сопла
Видеоиллюстрация результата моделирования движения ступени (Blender 3D):
