В предыдущей части мы запустили двухступенчатую ракет в космос. Вторая ступень достигла космической скорости по формуле Циолковского и согласно законам Ньютона. Это, конечно, хорошо и правильно, но не совсем. Точнее не совсем правильно. В наших расчетах мы запускали ракету в белый свет, как в копеечку, вертикально вверх. В этом случае первая ступень улетает в открытый космос по инерции и летит, черт знает куда (а черт – потому что бога нет, Гагарин, когда летал, не видел).
Реальные ракеты выходят на орбиту по-другому, не вертикально вверх. После старта ракета начинает отклонятся программой управления с тем, чтобы при выходе на орбиту, она имела направление движения параллельно земле (по-грамотному это называется угол тангажа). Давайте сделаем модель, которая будет это учитывать. Если использовать методы структурного моделирования, это будет сделать не сильно сложнее, чем модель артиллерийского снаряда, который мы перехватывали в задаче про волка и зайца.
Методы структурного моделирования позволят нам создать набор компонентов, из которых, как из кубиков лего, можно собирать одну-, двух- и трехступенчатые ракеты.
А для того, чтобы наша ракета была не абстрактная, возьмём данные по ракете «СОЮЗ», к тому же на хабре уже есть решение этой задачи. Больше спасибо автору, что уже собрал все необходимые данные. https://habr.com/ru/articles/649961/
Тем, кто первый раз пытается создать структурную модель, и кому покажутся сложными физическая модель сферического коня в вакууме или численное интегрирование обыкновенных дифференциальных уравнений, я рекомендую почитать статью про противоракетную оборону Израиля, где все это объясняется на основе знаний математики 4 класса. https://habr.com/ru/articles/878168/
Поскольку мы создаем структурную модель, воспользуемся методикой создания структурных моделей, которая описана здесь https://habr.com/ru/articles/682662/ В том примере мы разбирали модель обогрева комнаты, сейчас по той же самой методике создадим модель ракеты «Союз».
Повторение – мать ученья, поэтому повторим основные стадии создания моделей.
Подготовительная стадия. Данная стадия осуществляется вне среды моделирования и включает следующие шаги:
1. Определение цели моделирования и требований к модели.
2. Разработка структуры модели.
3. Определение упрощения и приближения физической модели и формирование математических уравнений.
4. Сбор данных по параметрам модели, если возможно собрать реальные данные измерения для валидации результата моделирования.
Стадия создания модели. Это стадия разработки модели непосредственно в среде моделирования, после того как подготовительная работа выполнена. Стадия включает следующие шаги:
5. Создание и отладка модели, компонентов модели.
6. Интеграция компонентов в комплексную модель и ее отладка.
7. При наличии возможности сравнение с реальными данными (валидация модели).
8. Корректировка модели по экспериментальным данным.
9. Если модель не отвечает поставленным целям, то возврат на пункт 1, 2 или 3.
Приступим:
1 Цель моделирования – получить удовольствие от повторения расчета траектории реальной ракеты Союз.
2 Разработка структуры модели.
Модель будет состоять из:
1) Ступеней ракеты (чтобы можно было собирать многоступенчатые схемы)
2) Программы полета (изменение угла ориентации, включение ступеней и расстыковка ступеней)
3) Блока расчета гравитации в зависимости от положения ракеты в пространстве.
4) Блока расчета сопротивления в атмосфере.
3. Физические упрощения
Модель представляется в виде точки, учитываются действия сил сопротивления воздуха в атмосфере, кривизны земли, зависимости тяги двигателя от высоты. Далее цитата:
«За физическую модель космической ракеты принимается материальная точка, то есть физические размеры ракеты в расчётах не участвуют. Форма планеты Земля принимается за шар, место старта - точка на поверхности этого шара. Прямоугольная система координат привязана к плоскости, касательной шару в точке запуска с началом координат в точке запуска (далее - плоскость запуска). Также принимается, что траектория полёта есть плоская кривая, поэтому для расчётов будут использованы только две координаты: горизонтальная координата X - дальность, вертикальная координата Y - высота полёта относительно плоскости запуска.»
4. Сбор данных
Данные мы возьмем из статьи про «Союз», тем более там уже все сгруппировано и даже данные телеметрии есть.
Поздравляю! Подготовительная стадия выполнена! Теперь можно начинать создавать компоненты модели. Переходим к пункту 5.
Начнем создание многоступенчатой ракеты» Союз» с создания модели одной ступени.
Модель ступени ракеты
Ступень ракеты будет рассчитывать силу (тягу двигателя) и изменение массы, чтобы можно было собирать многоступенчатую ракету, как из кубиков лего.
Поставим на схему блок «Субмоде��ь» и вызовем контекстное меню, нажав правой кнопкой мыши:

В появившемся окне добавим параметры, которые будут определять свойства ступени ракеты. У меня их 3: номинальная тяга, номинальный расход топлива, стартовая масса. См. рисунок 2. Эти параметры я и занес в качестве начальных.

Для каждой ступени будут свои параметры. Но для отладки мы можем рассматривать всю ракету на первой стадии расчета до отделения первой ступени, как одну большую ступень. Это нам позволит сравнить результаты с полученными ранее и телеметрией.
Рассмотрим движение ракеты в произвольный момент времени после старта. Поскольку в отличие от первого варианта мы стартуем не вертикально, нам придется рассматривать движение ракеты в координатах Х Y. За центр координат мы берем точку старта.
Тогда текущее положение относительно точки старта будет Xc Yc (см. рис. 3)

Зная текущее положение и угол наклона ракеты в каждый момент времени, мы можем рассчитать проекции силы тяги и силы тяжести на оси X и Y.
Сначала разберемся с силой тяги:
Для массы ракеты уравнения еще проще: зная расход топлива в секунду, мы используем блок «интегратор», в начальных условиях которого указываем начальную массу.
Войдем внутрь субмодели и соберем схему как на рисунке 4.

Небольшое пояснение для расчета массы. Все ровно так же, как со скоростью. Если мы подаем на вход интегратора скорость движения в м/с, то на выходе будет пройдённый путь. Если мы подаем на вход интегратора скорость расхода топлива в кг/с, то на выходе будет масса потраченного топлива. В нашем случае масса уменьшается, поэтому мы в свойствах источника указываем параметр с отрицательным знаком, в столбце формула. (см. рис. 5)

В блоке «интегратор» в качестве начальных условий мы указываем стартовую массу, таким образом во время расчета масса будет уменьшатся. (см. рис. 6)

Даже такой простой блок уже можно протестировать в автономном режиме.
Выйдем из блока и подключим к нему программу управления. В нашем случае на первом этапе – это просто программа изменения угла. Используем блок из закладки Источники «Кусочно-линейная», в котором зададим угол поворота 0 в начале расчёта и 61 град при 118 секунде, как указано условиях для первой стадии полета. Чтобы внутри блока считался синус и косинус, необходимо пересчитать градусы в радианты. Для этого поставим блок усилитель, в параметрах которого коэффициент усиления pi/180. Подключим графики и получим схему как на рисунке 7.

Eсли запустить на расчет, силы нам покажут дугу, а вот изменение массы будет линейным согласно заданной скорости расхода, и ее уже можно сравнить с исходными данными и убедится, что массу мы теряем правильно.
У нас в расчете получилась общая масса на 118 секунде – 113 816 кг, что в принципе соответствует расчету уменьшения массы при заданном постоянном расходе 313 000 – 1688 *118 = 113 816 кг.
Один компонент мы, похоже, отладили!
Модель действия сил тяжести(гравитации)
На ступень, кроме силы тяги, действует сила тяжести, так же как она действовала на артиллерийский снаряд в задаче про волка и зайца. Нам нужно рассчитать действие силы тяжести (гравитации). В отличие от предыдущей задачи, где старт идет вертикально вверх, здесь идет поворот в процессе полета, а значит направление силы по отношению �� направлению полета тоже меняется. Считаем, что сила тяжести направлена к центру земли, тогда при повороте ракеты и ее смещении по оси Х, будет меняться угол наклона вектора силы тяжести относительно вертикальной оси Y - b. (см. рис. 8 )

Зная текущие координаты можно рассчитать
и
где: – расстояние от цента земли до ракеты
- радиус земли
Из рисунка 7 следует также, что высота ракеты над поверхностью земли h во время полета может быть рассчитана как:
Добавим на схему блок «язык программирования» и запишем эти выражения в виде скрипта (см. рис. 9)

Данный блок получает текущую позицию ракеты и по ней рассчитывает составляющие ускорения притяжения по осям X и Y, и высоту ракеты над уровнем моря.
6. Интеграция компонентов модели
Это пункт 6 из методики структурного моделирования.
Теперь, когда у нас есть расчет силы тяги и силы тяжести, мы можем собрать уравнения Ньютона. У меня получилась схема, как на следующем рисунке:

Разберёмся, как работает эта схема. В модели первой ступени мы рассчитываем силу, разложенную по осям X и Y и массу ступени c учетом сжигания топлива. Далее мы размножаем массу, поскольку нужно считать два ускорения по двум осям. Поделив силу тяги на массу, мы получаем ускорение по оси Х и Y от работы двигателя. К этому ускорению нужно добавить ускорение от силы тяжести. После этого мы ставим интегратор и на выходе из него получаем скорость движения по осям Х и Y, а еще раз интегрируя, получаем положение. Это положение мы передаем в блок гравитации, где рассчитывается направления ускорения свободного падения и высота.
7. Сравнение с реальными данными (валидация)
Выполняем пункт 7 методики.
Запускаем на расчет и получаем немного не то, что у нашего эталонного расчета и показания телеметрии. У меня получилось 62,375 км. И внимательный читатель, конечно, догадался: в нашей модели мы используем максимальную силу тяги в вакууме – 4997904.
Если поставить в модель ступени силу тяжести на старте из условий 4147360, то высота в конце первой стадии будет – 40,10 км. Автор исходного расчета предлагает сделать линейное изменение тяги по времени, когда за 118 секунд сила тяги увеличивается с минимальной до максимальной.
8. Корректировка модели ступени по реальным данным
Переходим к пункту 8 методики структурного моделирования. Уточним модели для учета изменения силы тяги на земле и в вакууме. Это можно реализовать, добавив в модель к номинальной тяге линейный множитель, зависящий от времени. И тогда за 118 секунд тяга увеличится с минимальной до максимальной (см. рис. 11)

Еще одна верификация и валидация.
После доработки наш расчет полностью совпадает с расчетом, выполненным в исходной статье, высота первой стадии полета у меня получилась 47,996 км »48 км. Для проверки скорости добавим в модель структуру расчета абсолютной скорости, которая получается как корень квадратный из суммы квадратов скоростей по х и у. Схема расчета у меня получилась как на рисунке 12.

Скорость получилась равной 6699 км/ч (чуть меньше, чем в эталонном расчете).
Для расчета высоты над уровнем земли и расстояния, измеренного по земле исходя из точки старта, можно воспользоваться теми же соотношениями, которые мы использовали для расчета гравитационного ускорения (см. рис. 8).
Синус угол отклонения силы тяжести направленной к центру:
Где Rr – расстояние от центра земли до ракеты:
= 6.365х106 - радиус земли
Из рисунка 7 следует также, что высота ракеты над поверхностью земли h во время полета может быть рассчитана как:
Расстояние по земле – это длинна дуги, которая рассчитывается как , если угол
выражен в радианах.
Зная значения для синуса, мы можем использовать блок арксинус, и тогда схема расчета положения ракеты относительно поверхности будет выглядеть, как показано на рисунке 13:

Общая модель представлена на рисунке 14:

Теперь мы можем проверить расстояние от точки старта на момент отделения первой ступени 51, 845 – на 155 метров ближе, чем эталонный расчет.
Сравнение с эталонными данными показывает, что модель работает верно.
Продолжаем корректировку модели
Но эталонный расчет не учитывает силу трения воздуха. Давайте добавим ее в расчет. Сила трения рассчитывается по формуле:
Где – плотность воздуха,
– коэффициент сопротивления для конической поверхности;
- площадь поперечного сечения ракеты,
- скорость ракеты.
Для корректного расчета нам необходимо также разложить силу трения на составляющие по оси Х и Y, а значит нам нужно передать в блок расчета угол поворота ракеты. В итоге у меня получается скрипт расчета, представленный на рисунке 15

Для этого достаточно знать плотность воздуха и скорость движения ракеты. Плотность у нас рассчитывает специальный блок «атмосфера». Собираем все это в одну субмодель, у меня получается схема как на рисунке 16:

А общая модель с учетом сопротивления воздуха представлена на рисунке 17. Для того чтобы не плодить линий связи, высоту над поверхностью я передаю в блок в память с именем h, чтобы потом использовать в любом месте. И используем это значение, передав его в блок расчета сопротивления.

Для расчета перегрузки мы используем блок аналогичный расчета скорости, только в место пересчета м/с в км/ч введем коэффициент 1/g (см. рис. 18).

Результаты моделирования траектории полета представлены на рисунке 19

Валидация модели
Сравнение результатов моделирования с эталонными данными и показаниями телеметрии
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 45 | 48 | 45,296 |
Дальность, км | 48 | 52 | 50,232 |
Скорость, км/ч | 6312 | 6702 | 6478 |
Перегрузка, g | 4 | 4 | 4,06 |
Как видим, результат с учетом сопротивления воздуха получается ближе к данным телеметрии. Ура! Космос будет наш!
Сделаем еще одно уточнение. Мне не нравится, что в текущей модели эффективность двигателя снижается линейно по времени. В реальности, эффективность двигателя зависит от давления атмосферы, которое в свою очередь нелинейно зависит от высоты. Самое интересное, что в нашем блоке расчета атмосферы уже есть параметр давления на разной высоте. Остается только рассчитать коэффициент пропорциональности, который на нулевой высоте будет выдавать минимальную тягу из данных, а в вакууме –максимальную.
Тогда наша модель станет универсальной. Вместо 118 секунд увеличения тяги, которое задается в исходной модели, мы можем разгонятся как быстрее, так и медленнее, но тяга будет зависеть только от высоты ступени.
Добавляем в модель ступени вход высота и меняем линейное увеличение тяги на зависимость от высоты. Схема получается примерно, как на рисунке 20.

Для пересчета давления у меня получилась зависимость как на рисунке 20.
На графике рисунка 21. показано отличие коэффициентов линейного увеличения тяги двигателя в зависимости от высоты.


Повторим расчет и сравнение с эталонными и показателем телеметрии
Результаты расчета первой стадии полета:
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 45 | 48 | 46,206 |
Дальность, км | 48 | 52 | 51,342 |
Скорость, км/ч | 6312 | 6702 | 6648 |
Перегрузка, g | 4 | 4 | 4,06 |
Видно, что наша модель ближе к данным телеметрии, чем и эталонная, что не может не радовать.
Теперь, когда мы проверили и отладили модель одной ступени на первом участке полета, можно использовать эту же модель для исследования следующего участка полета.
Необходимо отметить, что у нас вторая ступень работает одновременно с первой и продолжает работать, когда первая уже отстыковывается. Поэтому нам необходимо создать программу управления, в которой отдельно указывается время включения двигателя, время отключения двигателя и время отделения ступени, и они будут разными для разных ступеней.
Соберем программу управления в отдельный блок субмодели, где и будем задавать время включения и выключения двигателей и время отделения ступеней.
Для первой и второй ступени достаточно указать время отключения и время отделения ступени. Для этого можно использовать блок ступенька, которая равна 1 на старте и становится равна 0 в момент отключения или отстыковки. Потому что отделение происходит одновременно с выключением двигателя.

Для того что бы отрабатывать команды, необходимо доработать модель ступени, которая будет принимать команды на включение и выключение двигателя и отстыковку ступени. У меня получается схема как на рисунке 24.

Два ключа, которые замкнуты при значении 1, обеспечивают передачу значения по линии связи, при 0 они разомкнуты и по линии предаётся 0. В этом случае при отключении двигателя скорость изменения массы не передается в интегратор, и масса остается постоянной. Второй ключ при размыкании обнуляет значение массы, предаваемой из модели в общую модель ракеты. Такую модель можно копировать и подключать в общую модель. (см. Рисунок 25)

У меня получилась примерно такая модель, как на рисунке 24. В ней рассчитывается первая ступень и вторая ступень.
Обратимся к данным из исходной статьи:
public StageTwo() {
t = 182.0; // Время работы второй ступени после отделения первой, с
F = 921200.0; // Тяга двигателя второй ступени, Н
M = 96752.0; // Масса ракеты с космическим кораблём после отделения первой ступени, кг
angle = 16.0; // Конечный угол поворота ракеты относительно вертикальной оси за время работы второй ступени, град
k = 310.0; // Расход массы второй ступени, кг/с
calculateTimeValues();
Распределим массу, тягу и расход из одной ступени в две так, чтобы параметры соответствовали исходным данным. В константу для МА добавим массу третьей ступени и космического аппарата, чтобы общая стартовая масса осталась прежней. Теперь можно добавить угол поворота в программе полета, время работы второй ступени и запустить расчет на 300 секунд, чтобы получить траекторию двух первых стадий полета.
У меня получились вот такие красивые графики:

Сравнение с результатом второй стадии полета показывает, что мы опять получили результат более приближенный к данным телеметрии.
Результаты расчета второй стадии полета:
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 154 | 164 | 154,962 |
Дальность, км | 452 | 527 | 503,975 |
Скорость, км/ч | 13 732 | 14 480 | 13 934 |
Перегрузка, g | 2.3 | 2.27 | 2,27 |
Добавляем третью ступень. Поскольку она работает уже за пределами атмосферы, то нам не нужно рассчитывать зависимость тяги от высоты, ее модель будет как на рисунке 27.

Используя данные из исходной статьи, заполним свойства ступени.
public StageThree() {
t = 224.0; // Время работы третьей ступени, с
F = 294000.0; // Тяга двигателя третьей ступени, Н
M = 30693.0; // Масса ракеты с космическим кораблём после отделения второй ступени, кг
angle = 8.0; // Конечный угол поворота ракеты относительно вертикальной оси за время работы третьей ступени, град
k = 84.7; // Расход массы третьей ступени, кг/с
calculateTimeValues(); }
Хьюстон у нас проблемы!
А вот с третьей ступенью у меня ничего не сложилось. Если взять тягу двигателя, указанную в задаче, то ни скорости, ни высоты у меня не получается, и перегрузка при такой тяге составляет всего 2.4 вместо 2.9, скорость 25 200 км/ч вместо 25 600 км/ч, высота 172 км вместо 202 км
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 202 | 200 | 173 |
Дальность, км | 1675 | 1738 | 1637 |
Скорость, км/ч | 26 737 | 27 019 | 25 270 |
Перегрузка, g | 2.9 | 2.61 | 2.4 |
Для меня это странно, поскольку сравнение идет с реальной телеметрией и до этого момента все совпадало. Что делать?
Методика создания структурных моделей на этом месте предлагает два варианта:
8 Корректировка модели по экспериментальным данным.
9 Если модель не отвечает поставленным целям, то возврат на пункт 1, 2 или 3
Вернемся на пункт 3 и поменяем исходные упрощения и приближения, принятые в математической модели.
В качестве гипотезы предположим, что нужно учитывать изменение силы тяжести при удалении от земли. Формула находится достаточно быстро в интернете:
Вот чем хорошо структурное моделирование: можно легко и просто уточнять модель. Достаточно просто поменять один блок в модели, не переделывая всю структуру. Переходим в модель расчета гравитации, где есть и радиус земли, и высота. Просто добавляем формулу зависимости тяжести от высоты в расчет. У меня получается вот такой код:

Теперь у нас сила тяжести уменьшается с высотой. Ракета должна полететь выше. Запускаем на расчет и видим, что наше предположение оправдалось:
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 202 | 200 | 206 |
Дальность, км | 1 675 | 1 738 | 1 631 |
Скорость, км/ч | 26 737 | 27 019 | 25 212 |
Перегрузка, g | 2.9 | 2.61 | 2.4 |
Получается по высоте ближе к правде (телеметрии), но скорость и ускорение ниже, чем у по телеметрии «Союза».
Вообще-то, все это я так долго собирал, чтобы посмотреть, как ракета «Союз» выходит на орбиту. Поэтому добавляю блок расчета поверхности земли (см. рисунок 29) и вывожу его результат на фазовый портрет.

Немного подправив модель, чтобы было поменьше линий связи, общая модель у меня получилась как на рисунке 30.

Логично, что ни хрена не получилось! Скорость меньше первой космической, см. рисунок 31.

Ничего, у Илона нашего Маска тоже ракеты не сразу полетели. Продолжаем исследования.
Я еще раз проверил данные и нашел ошибку в исходной статье:
· Расход массы третьей ступени: k3 = (27755 - 2355)/300 = 84,7 кг/с
Но ступень работает 224 секунды!
А значит. правильная скорость расхода на третьей ступени:
k3 = (27755 - 2355)/224 = 113.4 кг/с!
Модель-то уже готова, нужно поменять одну цифру и посмотреть, что получится. Сначала считаем 524, получаются такие результаты:
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 202 | 200 | 229 |
Дальность, км | 1 675 | 1 738 | 1 659 |
Скорость, км/ч | 26 737 | 27 019 | 29 505 |
Перегрузка, g | 2.9 | 2.61 | 5.5 |
Если в первом варианте была явно недостаточная скорость, то сейчас явный перебор. Явно должны выйти на орбиту. Меняем время расчета на 15000 сек и смотрим за результатом. Результат совсем не похож на правду. Ракета улетает на 2000 км, потом возвращается, цепляет атмосферу и благополучно втыкается в землю в самом начале второго витка! (см. рис. 32)

Оказывается, просто большая скорость нам не гарантирует нужного результата.
Опять корректируем модель по экспериментальным данным.
Возвращаемся обратно, меняем расход топлива как в условиях (хотя это не верно), но увеличиваем тягу, чтобы достичь заданных показателей. Так же пришлось поменять конечный на 92 градусов. И вот при тяге третьей ступени 335000Н, заданные параметры телеметрии практически идеально!
Результаты расчета третьей стадии полета:
| Данные телеметрии | Расчет эталонный | Расчет текущий |
Высота, км | 202 | 200 | 202.8 |
Дальность, км | 1675 | 1738 | 1678 |
Скорость, км/ч | 26 737 | 27 019 | 27 178 |
Перегрузка, g | 2.9 | 2.61 | 2.9 |
Решение задачи, которое совпадает c данными телеметрии на первых 524 секундах полета! Вот оно счастье! Но нет, опять опыт. Ракета с данными, похожими на телеметрию, отказывается выходить на орбиту!
Смотрите: после отделения третьей ступени аппарат поднимается на высоту 300 км, а потом спускается по орбите. При этом он касается атмосферы и начинает активно тормозиться, спускаясь до высоты 60 км.
Несмотря на то, что ему хватит инерции выйти еще раз на высоту 150 км, за атмосферу, но скорость уже ниже первой космической. Аппарат еще раз ныряет в атмосферу, где окончательно тормозится и падает.

Раз мы так высоко забираемся, может уменьшить угол? В исходной статье говорится о повороте на 8 градусов во время полета третей ступени, это вы итоге дает нам 85 в конце разгонного участка. Поскольку мы можем менять программу легко и непринужденно, сделаем угол 100. И добавим чуть тяги третей ступени до 350 кН. И тут случилось чудо! Я почувствовал себя Королевым, ракета вышла на орбиту:

Что характерно, при таком угле третья стадия полета не совпадает с исходными данными. У меня получается, что ракета набирает только 185 км:

Дальше можно бесконечно продолжать экспериментировать с разными ступенями, силами тяги и даже физическими законами.
Например, если убрать зависимость ускорения свободного падения от расстояния, то получается вот такая забавная траектория:

Модель можно взять здесь...
В принципе, можно попробовать осуществить посадку первой ступени и повторить модель ракеты Илона нашего Маска, просто меняя тягу и расход топлива первой ступени. И вообще, созданная нами модель позволяет задать параметры любой ракеты. А все потому что мы используем структурное моделирование.
Подписыватесь на канал технолог Петухов, там еще больше про математическое моделирование и без цензуры.
Видео с деомнстрацие модели: