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

Там очень подробно описано как, зная скорость объекта, можно рассчитать траектории движения различных объектов в пространстве.

 https://habr.com/ru/articles/878168/

В этот раз мы займемся исследованием траектории движения космических ракет.  Сравним формулу Циолковского с законом Ньютона и рассчитаем отправку груза на орбиту земли одноступенчатой ракетой, и двухступенчатой. И все это – в рамках курсов школьной физики и математики с помощью структурного моделирования.

Решим задачу с одноступенчатой ракетой двумя способами:

1) По формуле Циолковского, которая позволяет рассчитать скорость ракеты в любой момент времени:

где:

= 10 000 кг – полезная нагрузка;

 = 200 000 кг – начальная масса топлива;

 = 5 000 м/с – скорость уходящих газов;

= 5 000 кг/с – расход топлива;

– тяга двигателя ракеты.

Те, кто читал предыдущую статью, уже, наверное, догадались, что раз у нас есть формула скорости, то мы поставим интегратор, который нам все проинтегрирует, и мы получим высоту полета. Все примерно так, но по-другому.

Вместо типовых блоков для создания модели мы будем использовать специальный блок –  «Язык программирования». Этот блок работает почти так же, как все остальные: он забирает данные из линий связи на входе, выполняет вычисления и складывает результат в линии связи на выходе.  

Только вычисления задаются в виде скрипта на встроенном языке программирования.

Если кликнуть по блоку два раза, то вместо свойства, как у других блоков, появляется окно редактирования, где мы можем написать программу на языке SimInTech, так как мы это делали для общего скрипта программы, где заносили константы. (см. https://habr.com/ru/articles/878168/ )

Все достаточно просто. Сначала запишем параметры ракеты в список констант, разделенных запятыми, а дальше запишем формулу Циолковского. На выход отправим переменную скорость V

Во время расчета на каждом шаге моделирования этот блок будет выдавать значение скорости с учетом времени. Глобальная переменная time на каждом шаге отображает текущее модельное время.

Переменные, которые мы хотим выдать из блока, перечисляются после ключевого слова output. Если сохранить модель, нажав на зеленую галочку вверху, то в блоке исчезнет вход, а выход будет подписан как v.

Подключим к выходу нашей модели блок-интегратор. Как мы уже знаем, если на входе в интегратор – скорость, то на выходе будет путь, а в нашем случае высота.

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

Перейдем к настройке параметров расчета и установим конечное время расчета – 40 секунд. Обратим внимание на настройки шага интегрирования hmin, hmax = 0.001 и метод интегрирования (Эйлера).

Запускаем на расчет и видим, что у нас получилось в конце расчета:

скорость – 7.655 км/с; высота  - 111,7 км

Функционально-блочная диаграмма или язык программирования?

Вместо блока языка программирования мы могли бы набрать формулу Циолковского в виде блок-схемы, упаковав его в субмодель. У меня получилось примерно, как на следующем рисунке:

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

Забавно, но у меня при создании этой схемы получилось сократить вычисления. Внимательно посмотрев на схему, вы не найдете здесь вычисления силы тяги

А все потому, что в формуле Циолковского она делится на и ее можно сократить, и формула в этом случае приберет вид:

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

На численный результат вычислений это не влияет, цифры получаются те же самые. Но мысократили вычисления на одну операцию, а с учетом того, что мы их выполняем на каждом шаге интегрирования (в нашем случае 40 000 раз за один прогон), мы снизили электропотребление процессора и уменьшили наш углеродные след. А также отсрочили тепловую смерть вселенной, и это хорошо!

Получается, что хотя нарисовать диаграмму дольше, чем написать скрипт, оптимизировать и отлаживать диаграмму легче. Поэтому ответ на вопрос, что лучше диаграмма или скрипт, такой же, как и на вопрос: если жизнь на марсе, нет ли жизни на марсе? Это науке неизвестно!

Интегрировать или не интегрировать, вот в чем вопрос.

Вернёмся от красивой диаграммы к нашему волшебному блоку «Язык программирования». Мы не зря назвали наш блок волшебным, он может не только рассчитывать по заданному скрипту значения скорости ракеты в зависимости от времени, но может и сам считать интеграл, а значит может считать траекторию полета.

Скопируем блок и поставим на схему копию блока, которую будем модифицировать. Нам достаточно внести небольшие изменения в программу на языке программирования. Поскольку скорость у нас — это производная от пути, то мы можем так и записать выражение для производной пути: s’ = v. 

Также необходимо указать начальное состояние интегрируемой переменной, s = 0 указываем после ключевого слова init

После этого мы можем указать переменную s  в выходных переменных. Полный текст будет выглядеть, как на рисунке:

Теперь из блока у нас выходят две линии связи: по одной идет скорость, по второй высота. Сравним результаты расчета. Как бы это ни было удивительно, но они полностью совпадают! Волшебный блок работает.

Проверим Циолковского Ньютоном

Формула Циолковского – это конечно хорошо и здорово. Но, как сказал бог Аврааму, доверяй, но проверяй! 

Давайте посчитаем высоту полета и скорость этой же ракеты, не по готовой формуле Циолковского, а на основе физических законов, а именно закона Ньютона.

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

Закон Ньютона:

 Сила, действующая на тело равна массе тела, умноженной на ускорение его движения:

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

В первом приближении не будем учитывать сопротивление воздуха (в формуле Циолковского сопротивления воздуха тоже нет). Тогда на ракету у нас действуют всего две силы: отрицательная сила тяжести    и сила тяги : 

Тогда ускорение ракеты можно посчитать, зная массу и сумму приложенных сил:

где: - суммарная сила, действующая на ракету.

Ранее мы по скорости в блоке посчитали путь, теперь по ускорению посчитаем сначала скорость, а потом путь. 

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

В модели ракеты получаются три меняющиеся переменные:

- высота (s)

- скорость (v)

- масса (M)

Для каждой из них у нас есть формула расчета производной. Для все трех переменных нужно задать начальное состояние для интегрирования. Список констант, описывающих ракету, не поменяется, а вот текст программы будет выглядеть, как показано на рисунке:

Чтобы сравнить результаты расчета по формуле Циолковского, добавим на схему сравнивающий блок и выведем результат на график:

Запустим на расчет и, если вы все сделали правильно, мы увидим, что у нас результат расчета отличается.

  Скорость по формуле Циолковского - 7 654.7896 м/с, скорость по уравнениям Ньютона - 7 654.5896 м/с, погрешность 0.2 м/с.

Вот это поворот! И кто тут прав: Ньютон или Циолковский? Правильный ответ – оба!  И Циолковский прав, и Ньютон. А виноват в том, что ответы разные, другой великий ученый – старина Эйлер, наш российский ученый, похороненный в Санкт-Петербурге. 

Дело в том, что для расчета интеграла мы использовали метод Эйлера, описание которого подробно приведено в предыдущей статье. Изобретенный 250 лет назад, он по умолчанию содержит в себе отклонения, когда плавная непрерывная функция считается по шагам. И поэтому скорость, рассчитанная по формуле Циолковского, более точная, чем полученная интегрированием по методу Эйлера.

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

Результат: отклонение также уменьшилось в 10 раз, было 0.2 стало 0.02. Но скорость вычисления значительно замедлилась. Мы стали считать только в 1.5 быстрее реального времени, а это очень медленно. Скорость можно посмотреть в нижней части схемного окна. 

Кажется логичным, что чем больше шагов по времени, тем больше вычислений требует метод Эйлера. Можно сократить вычисления, уменьшить шаг, тогда у нас будет меньше вычислений. Например, если задать шаг интегрирования 0.1, то вычисление будет мгновенным, нажми на кнопку – получишь результат. Скорость расчета в 1300 раз быстрее реального времени. Но погрешность расчета уже 20 м/с.

Получается, как в анекдоте:

- С какой скоростью вы можете печатать?

- 20 слов в минуту.

- А можете 50 слов в минуту?

- Я могу и 100 слов минуту, но такая фигня получается. 

Казалось бы, некуда деваться: либо быстро, либо точно, либо одно из двух. Но, на наше счастье, еще Ломоносов отмечал:

 «…что может собственных Платонов

 И быстрых разумом Невтонов 

Российская земля рождать.” 

Поэтому не Эйлером единым ограничивается вклад Российских и Советских ученых в методы расчета динамических систем.  

Представляю вам Леонида Марковича Скворцова. Его методы используются в SimInTech и, если раскрыть свойства «Методы интегрирования», то мы увидим целую пачку методов с непонятными аббревиатурами SRK2m1 или ARK32. Так вот все эти методы разработаны именно Леонидом Марковичем.

И нельзя сказать, что он секретный физик. Специалисты в области численного решения дифференциальных уравнений его неплохо знают. Вот, например, статья в цифровой библиотеке, на сайте Гарварда от 2022 года. «Диагонально-неявные методы Рунге-Кутты третьего и четвертого порядка для жестких и дифференциально-алгебраических задач.»  

https://ui.adsabs.harvard.edu/abs/2022CMMPh..62..766S/abstract

А вот свежие статьи уже зарубежных авторов с цитатами из этой статьи

https://ui.adsabs.harvard.edu/abs/2022CMMPh..62..766S/citations

Ну, и максимально близкий к теме статьи доклад космического агентства NASA  DiagonallyImplicit Runge-Kutta Methods for Ordinary Differential Equations. A Review.

В этом докладе  -  5 упоминаний фамилии Skvortsov по тексту.

https://t.me/Tech_Petuhoff/655

Для чего нужны эти методы? Они обеспечивают возможность менять шаг интегрирования в процессе расчета, чтобы можно было ускорить расчет, не теряя точности. Увеличивая шаг интегрирования, когда скорости процессов маленькие или уменьшать, когда скорости большие.

И пусть вас не вводит в заблуждение имя Рунге-Кутта в названии. Так называют целый класс методов, а вот конкретная реализации в этом классе — это метод Скворцова.

Давайте выберем один из методов Скворцова, например (Адаптивный 3).

Зададим минимальный шаг – 0.0001 и максимальный шаг – 0.1 

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

Запустим на расчет и восхитимся результатами, полученными с помощью советских и российских математиков.

Скорость расчета в 1200 раз быстрее реального времени, а ошибка - 4.6x10^-5

 Сравните с методом Эйлера. Счет, как говорится на табло! Именно поэтому методы Скворцова используются в NASA и в SimInTech. Так что наш вклад в мировую космонавтику – не только «Юра мы все просрали!»

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

Воспользуемся той же формулой, которую мы использовали в статье для снаряда:

Где: v - скорость тела, k - коэффициент сопротивления воздуха.

Где: S - площадь лобового сопротивления снаряда; - плотность воздуха;

Плотность воздуха в зависимости от высоты может рассчитываться специальным блоком в SimInTech.

 Чтобы забрать плотность в модель ракеты, мы должны прописать служебное слово input и дать имя переменной, которую мы будем забирать из блока расчета атмосферы. После этого достаточно добавить недостающие параметры для расчета силы трения воздуха и добавить силу в уравнения:

Для корректной работы схемы требуется также задержка на шаг интегрирования перед блоком расчета свойств атмосферы. Поскольку для расчета высоты нужна плотность на данной высоте, а для расчета плотности нужна высота, это называется алгебраическая петля. Задержка на шаг позволяет исправить это противоречие. В ней мы задаем начальную высоту, а далее на каждом следующем шаге для расчета полётности мы используем значение высоты с предыдущего шага, сохраненное в задержке.  В итоге у нас получается вот такая схема:

Запускаем на расчет и получаем следующие результаты:

Отклонение, связанное с трением, у нас растет с ростом скорости. Чем больше скорость, тем сильнее трение, и тем сильнее отличается скорость по формуле Циолковского от нашего физического расчета. По мере увеличения высоты плотность воздуха падает, примерно на 35 секунде (высота 80 км) трение уже практически не снижает скорость, и отклонение остается постоянным (у нас это примерно 172 м/с).

По высоте отличие составляет 9 км, по идеальной формуле Циолковского высота 111.7 км, по физическому расчету – 107.3 км

Но главный результат нашего расчета то, что ракета не вышла на первую космическую скорость, это какая-то ФАУ-2 получается, а не космическая ракета. Мы мирные люди и нам нужен мирный космос. Поэтому добавим вторую ступень!

Вторая ступень

Двухступенчатая ракета – это, по своей сути, две ракеты, поставленные одна на другую. Сначала первая ракета (первая ступень) разгоняется, насколько ей хватает топлива, потом запускается вторая ракета (вторая ступень), а первая благополучно падает на землю (или садится, если это ракета Илона нашего Маска).

Поскольку мы находимся в среде динамического моделирования, то созданную модель можно просто скопировать и получить сразу две ракеты. А если их поставить одну на другую, мы и получим двухступенчатую ракету.  

Конечно, чтобы одна из них была второй ступенью, придется немного модифицировать программу в блоке.  

Тут у нас появляется программа полета. Мы должны иметь возможность включать и выключать двигатель. Для этого добавим новую входную переменную enginecontrol, добавим ее после служебного слова input

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

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

Теперь у нас расход топлива и сила тяги появляются только тогда, когда управляющее воздействие больше 0.

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

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

В текст программы внесем следующие изменения:

1)        В параметрах ракеты зададим данные по второй ступени

2)        Для учета скорости и положения будем учитывать набранные первой ступенью скорость и высоту.

3)        Пока двигатель не включен (enginecontrol = 0), ускорение и скорость также равны 0.

4)        Как только двигатель включается (enginecontrol = 1), начинается расчет движения второй ступени. 

5)        В качестве результата отдаются суммы скоростей и высоты первой ступени и второй. 

Теперь, когда вторая ступень содержит в себе все необходимое для старта, займемся программой полета. 

У нас сейчас все очень просто: первые 32 секунды работает двигатель первой ступени, затем он выключается и включается двигатель второй ступени. Такую программу можно реализовать с помощью блока «ступенька», который выдает 1 в течение 32 первых секунд и 0 после.

Для первой ступени мы используем это сигнал напрямую, для второй мы его инвертируем. Тогда во второй ступени первые 32 секунды полета будет 0, потом 1. Далее нам нужно запомнить высоту и скорость, на которой происходит переключение режимов, и передать их из первой ступени во вторую. Для этого у нас есть блок запоминания сигнала. Управлять запоминанием мы будем той же ступенькой, что и двигателем. Пока сигнал 1, передаются высота и скорость, как только сигнал 0, набранные скорость и высота запоминаются. И используются во второй ступени.  Схема расчета достаточна проста:

Если запустить на расчет, то мы увидим, что вторая ступень достигает первой космической скорости с запасом и выводит спутник на орбиту!

Так среда структурного моделирования позволяет моделировать полет ракеты даже школьнику.

 И еще, благодаря структурному моделированию, у нас теперь есть новый блок – модель ступени, и отдельная программа включения и выключения двигателей. 

Модели из статьи можно взять здесь...

Можно попробовать посчитать что-то приближенное к реальности, например ракету Союз но это тема следующей статьи. 

Подписывайтесь на канал Технолог Петухов https://t.me/Tech_Petuhoff/680 там без цензуры. (Кормящим матерям и беременным женьщинам обоего пола минздрав не рекомендует)