![](https://habrastorage.org/files/a4c/06c/f0d/a4c06cf0dd524789a0131cf391f03be9.png)
Перевод поста Джеффри Брайанта (Jeffrey Bryant), Пако Джейна (Paco Jain) и Майкл Тротта (Michael Trott) "Hidden Figures: Modern Approaches to Orbit and Reentry Calculations".
Код, приведенный в статье, можно скачать здесь.
Выражаю огромную благодарность Полине Сологуб за помощь в переводе и подготовке публикации
Содержание
— Размещение спутника в определенном месте
— Константы и первичная обработка
— Вычисления
— Построение графика
— Как рассчитываются орбиты сегодня
Моделирование возвращаемого спутника
Вышедший недавно в кинотеатрах фильм Скрытые фигуры получил хорошие отзывы. Действие разворачивается в важный период истории США; в нем затрагивается также ряд тем вроде гражданских прав и космической гонки. В центре повествования — история Кэтрин Джонсон и ее коллег (Дороти Воган и Мэри Джексон) из NASA в период развертывания программы Меркурий и ранних исследований пилотируемых космических полетов. Внимание также акцентируется на драматической борьбе за гражданские права афро-американских женщин в NASA, происходившей в то время. Компьютеры в то время едва появились, так что способность Джонсон и ее коллег решать сложные навигационные задачи орбитальной механики без использования компьютера обеспечили важную проверку ранних компьютерных результатов.
![](https://habrastorage.org/getpro/habr/post_images/d86/741/6f7/d867416f77b14b01668cbe99343a66b0.jpg)
Я остановлюсь на двух аспектах ее научной работы, упомянутых в фильме: вычислениях орбиты и расчетах, связанных с вхождением в атмосферу. Для орбитальных вычислений я сначала сделал ровно то же, что и Джонсон, а затем применил более современный прямой подход с использованием инструментов Wolfram Language. В фильме упоминается о решении дифференциальных уравнений методом Эйлера; я же буду сравнивать этот метод с более современным и вычислю возвратную траекторию с помощью данных модели атмосферы, полученных непосредственно из Wolfram Language).
В фильме не освещаются детали решения математических задач командой Джонсон, однако того, что вы увидите, хватит, чтобы прочувствовать и сравнить подход, отображенный в фильме, с тем, который существует на данный момент.
Размещение спутника в определенном месте
Одна из самых ранних работ, в которой Джонсон была соавтором ("Определение азимутального угла в точке выгорания для размещения спутника над выбранной точкой") имеет дело с проблемой проверки, может ли спутник быть размещен над определенной точкой Земли после определенного количества витков вокруг орбиты при определенной стартовой позиции (например, мыс Канаверал, штат Флорида) и траектории вращения. Подход, который команда Джонсон использовала для определения азимутального угла (угол, образованный вектором скорости космического аппарата в момент отключения двигателя, с фиксированным направлением, скажем, на север) для запуска ракеты основывался на других орбитальных параметрах. Это важный шаг для обеспечения правильного местонахождения астронавта для входа в атмосферу на Землю.
![](https://habrastorage.org/getpro/habr/post_images/b91/20f/518/b9120f518a0f40715ff057758c121e39.jpg)
Константы и первичная обработка
В статье Джонсон определяет ряд констант и входных параметров, необходимых для решения задачи вручную. Отступлю немного, чтобы объяснить термин «выгорание», относящийся к остановке ракетного двигателя. После «выгорания» параметры орбиты как бы «замораживаются», и космический аппарат движется только под действием силы тяжести Земли (в соответствии с законами Ньютона). В этом разделе я следую тем определениям единиц, которые были приняты в статье.
![](https://habrastorage.org/getpro/habr/post_images/b5f/cf8/5ef/b5fcf85ef2b5022e4dcccf476d643672.png)
Некоторые функции для удобства определены, чтобы взаимодействовать с углами в градусах вместо радиан. Это позволяет оптимизировать работу при вычислении углов:
![](https://habrastorage.org/getpro/habr/post_images/f1c/644/e2b/f1c644e2b233878fccccd599d3a71847.png)
![](https://habrastorage.org/getpro/habr/post_images/fbd/bb5/81b/fbdbb581ba610e3c842ee97099797f4a.png)
Джонсон описывает несколько других производных параметров, хотя интересно отметить, что иногда она принимает определенные значения для них, а не использует значения, полученные с помощью ее формул. Принятые ею значения часто были близки к значениям, полученным по формулам. Для простоты здесь используются значения из формул.
Semilatus rectum эллиптической орбиты:
![](https://habrastorage.org/getpro/habr/post_images/7ad/ee1/49b/7adee149b93321a66b7032ea3bcc4dbc.png)
Угол в плоскости орбиты между перигеем и точкой выгорания:
![](https://habrastorage.org/getpro/habr/post_images/fdb/67e/1ac/fdb67e1ac534a5a8ea32ed3b99a5b766.png)
Эксцентриситет орбиты:
![](https://habrastorage.org/getpro/habr/post_images/87c/432/bf3/87c432bf3b3fec7a12897b64e5df5dce.png)
Период орбиты:
![](https://habrastorage.org/getpro/habr/post_images/7ac/3fb/cf7/7ac3fbcf721d3f3c0f528bc07c7d7fb0.png)
Эксцентрическая аномалия:
![](https://habrastorage.org/getpro/habr/post_images/c30/ba6/5ab/c30ba65abf6c07d85fac06912e21b589.png)
Для описания следующего параметра проще будет процитировать оригинальную статью: "требование о том, чтобы спутник с точкой выгорания ?1, ?1 проходил через выбранную позицию ?2, ?2 после завершения n прохождений по орбите эквивалентно требованию, в соответствии с которым в ходе первого оборота спутник проходит через эквивалентное место с широтой ?2 так же, как и при выбранном положении, но с долготой ?2, смещенной в восточном направлении от ?2 на величину, достаточную для компенсации вращения Земли в течение n полных вращений, то есть полярного часового угла n?ЕТ. Долгота этой эквивалентной позиции, таким образом, определяется соотношением":
![](https://habrastorage.org/getpro/habr/post_images/a1f/745/b63/a1f745b63e65230247bc346c6859220a.png)
Время от перигея для угла ?:
![](https://habrastorage.org/getpro/habr/post_images/2b9/b9a/6a9/2b9b9a6a96b961ee711b277e8c0b9e99.png)
Вычисления
Часть окончательного решения заключается в определении значения для промежуточных параметров ??1-2е и ?2е. Это может быть сделано несколькими способами. Во-первых, я могу использовать ContourPlot для получения графического решения с помощью уравнений 19 и 20 из статьи:
![](https://habrastorage.org/getpro/habr/post_images/dea/217/61f/dea21761fc0b3d1b3422c0d756d7bdb7.png)
FindRoot можно использовать для нахождения численных решений:
![](https://habrastorage.org/getpro/habr/post_images/8bc/5ef/053/8bc5ef05321437daa502016f8d5b50c7.png)
![](https://habrastorage.org/getpro/habr/post_images/91a/f4e/353/91af4e35327f13131cd882c89fb44a89.png)
Конечно, у Джонсон не было доступа к функциям ContourPlot или FindRoot, так что в ее статье описывается итерационный метод. Я перевел методику, описанную в статье, на язык Wolfram, а также раскрыл ряд других параметров с помощью ее итерационного метода. Поскольку базовые выкладки рассчитаны для сферической формы Земли, поправки на сплющенности включены:
![](https://habrastorage.org/getpro/habr/post_images/9a9/8c1/7ed/9a98c17ed8661e348ae82550c3da681b.png)
Графическое изображение значения ?2е для различных итераций демонстрирует быструю сходимость:
![](https://habrastorage.org/getpro/habr/post_images/d84/58c/b44/d8458cb442f581dc2e4e03c038706488.png)
![](https://habrastorage.org/getpro/habr/post_images/b37/df8/489/b37df848934362ccd0da6a792b31a22d.png)
![](https://habrastorage.org/getpro/habr/post_images/770/089/949/77008994948c14f77b5272bd21bbf316.png)
![](https://habrastorage.org/getpro/habr/post_images/74f/d13/648/74fd1364871f44dbc59f38d3a7302021.png)
Я могу преобразовать метод в команде FindRoot следующим образом:
![](https://habrastorage.org/getpro/habr/post_images/e7d/f62/cb6/e7df62cb62a3dce34bfee68a59830319.png)
![](https://habrastorage.org/getpro/habr/post_images/227/25e/500/22725e500d2736f396d7402905157fab.png)
![](https://habrastorage.org/getpro/habr/post_images/ac9/363/a7d/ac9363a7d1726716e77d56450b25bb17.png)
Интересно, что даже итерационные шаги поиска корня этой более сложной системы сходятся довольно быстро:
![](https://habrastorage.org/getpro/habr/post_images/196/8e3/47b/1968e347bdc925a3e8d1eb408e792b77.png)
![](https://habrastorage.org/getpro/habr/post_images/b5b/c05/bff/b5bc05bff18edc80652671b76bec7c5f.png)
Построение графика
Когда орбитальные параметры определены, можно визуализировать решение. Во-первых, из предыдущих решений должны быть извлечены некоторые параметры:
![](https://habrastorage.org/getpro/habr/post_images/edc/3a8/9a8/edc3a89a89f81f22ccbb0c5ac6b634ae.png)
Далее должны быть получены широта и долгота спутника в зависимости от азимутального угла:
![](https://habrastorage.org/getpro/habr/post_images/980/fca/b23/980fcab23e937ad15cc18180c7ddb629.png)
![](https://habrastorage.org/getpro/habr/post_images/8cb/a85/b8e/8cba85b8ee7d482fb32be6249b62e492.png)
?s и ?s — широта и долгота как функция ?s:
![](https://habrastorage.org/getpro/habr/post_images/d51/d77/a36/d51d77a363bded21d402e5427b4f483a.png)
![](https://habrastorage.org/getpro/habr/post_images/9d5/2bd/dd1/9d52bddd1bf4460ae0419dbcb2146e85.png)
Наземный трек спутника может быть построен путем создания таблицы точек:
![](https://habrastorage.org/getpro/habr/post_images/adb/4ce/3b8/adb4ce3b8655f639ce2abcc28759b115.png)
В статье Джонсон представлена схема орбитального решения, включая маркеры, отмечающие выгорание, а также выбранные и эквивалентные положения. Подобную схему легко воспроизвести:
![](https://habrastorage.org/getpro/habr/post_images/e8d/58b/345/e8d58b3452f1b235ab8ced18a27376bf.png)
![](https://habrastorage.org/getpro/habr/post_images/538/898/11a/53889811a19abe965ed6755804ed46ae.png)
![](https://habrastorage.org/getpro/habr/post_images/095/3c7/3cd/0953c73cdbcaaf8357a62118d0f43f1a.png)
Вот ее оригинальная схема для сравнения:
![](https://habrastorage.org/getpro/habr/post_images/e34/1db/8b4/e341db8b46ffbd99fda0779cb3f1dac0.png)
Более визуально полезная версия может быть построена с помощью функции GeoGraphics, которая преобразует геоцентрические координаты в геодезические:
![](https://habrastorage.org/getpro/habr/post_images/cb1/412/61e/cb141261ef7578e1995aebf4ae687559.png)
![](https://habrastorage.org/getpro/habr/post_images/ef9/fc3/954/ef9fc3954904b7cc49454e89dd184072.png)
![](https://habrastorage.org/getpro/habr/post_images/7b8/f8d/325/7b8f8d325aa2aeb122112670fa836f9f.png)
Как рассчитываются орбиты сегодня
Сегодня практически у каждого из нас есть доступ к вычислительным ресурсам гораздо более мощным, чем те, что были у НАСА в 1960-е годы. Сейчас, используя только настольный компьютер и Wolfram Language, вы можете легко найти прямые численные решения задач орбитальной механики вроде тех, которыми занимались Кэтрин Джонсон и ее команда.
Чтобы вычислить азимутальный угол ? с помощью более современных методов, давайте установим параметры для простой круговой орбиты, берущей свое начало после выгорания над Флоридой; предположим также, что Земля сферически симметрична (я не буду пытаться искать точное соответствие с данными из статьи Джонсон и переопределю некоторые величины, используя современную систему единиц СИ). Начиная с той же высоты околоземной орбиты, используемой Джонсон, и с помощью сферической тригонометрии несложно вывести начальные условия для нашей орбиты:
![](https://habrastorage.org/getpro/habr/post_images/45b/c6a/a3a/45bc6aa3a7565d86d1fa153f414cb282.png)
![](https://habrastorage.org/getpro/habr/post_images/b79/e32/30d/b79e3230d84bc368acd014ab9ffddc6c.png)
![](https://habrastorage.org/getpro/habr/post_images/186/b29/5c4/186b295c47fedb40e21053bab3e95a9c.png)
Соответствующие физические параметры могут быть получены непосредственно с помощью Wolfram Language:
![](https://habrastorage.org/getpro/habr/post_images/cca/87b/1e2/cca87b1e2c8cbc4fa4f23cb41fa3a6c3.png)
Далее я получил дифференциальное уравнение движения нашего космического корабля с учетом гравитационного поля Земли. Есть несколько способов моделирования гравитационного потенциала вблизи Земли. Предположим, что планета сферически симметрична и используем декартову систему координат:
![](https://habrastorage.org/getpro/habr/post_images/780/f26/3bb/780f263bb27456923c080b0d160584c7.png)
Кроме того, вы можете использовать более реалистичную модель гравитации Земли, где в качестве формы планеты берется сплюснутый эллипсоид вращения. Точная форма потенциала такого эллипсоида (при условии постоянной плотности массы эллипсоидальных оболочек) доступна с помощью EntityValue:
![](https://habrastorage.org/getpro/habr/post_images/b45/226/a9f/b45226a9f44b303fdf43a3b3bf65398b.png)
Для общего однородного трехосного эллипсоида потенциал содержит кусочные функции:
![](https://habrastorage.org/getpro/habr/post_images/d9a/943/71c/d9a94371c838ff0d3b372a912312200a.png)
![](https://habrastorage.org/getpro/habr/post_images/a4a/0b3/713/a4a0b37139459ddcc02f5cbd4c3aec9d.png)
Здесь ? является наибольшим корнем выражения х2/(a2 + ?) + у2/(b2 + k) + z2/(c2 + ?) = 1. В случае, если мы имеем дело со сплющенным эллипсоидом, предыдущая формула может быть упрощена до элементарных функций…
![](https://habrastorage.org/getpro/habr/post_images/58a/71e/682/58a71e682a1baee8320b4b7d739616dc.png)
![](https://habrastorage.org/getpro/habr/post_images/503/5ca/658/5035ca65808d50e10a0d1ab4bcc3fb7e.png)
… где ?=((2 z2 (a2-c2+x2+y2)+(-a2+c2+x2+y2)2+z4)1/2-a2-c2+x2+y2+z2)/2.
Более простая форма, которую я буду использовать здесь, задается так называемой общепринятой формулой ускорения силы тяжести (IGF). В данной формуле учитываются отличия от сферически симметричного потенциала до второго порядка в сферических гармониках, и даются численно неотличимые от точного потенциала, на который ссылались ранее, результаты. В терминах четырех измеренных геодезических параметров потенциал IGF может быть определен следующим образом:
![](https://habrastorage.org/getpro/habr/post_images/b06/0df/5a3/b060df5a3ea4e308ff9559555fa42ff5.png)
Я мог бы легко использовать даже еще более точные значения для силы тяжести с помощью функции GeogravityModelData. В исходном положении потенциал IGF только на 0,06% отклоняется от приближения высокого порядка:
![](https://habrastorage.org/getpro/habr/post_images/ff7/3f7/b5c/ff73f7b5c019312168357aa065282954.png)
![](https://habrastorage.org/getpro/habr/post_images/71e/c28/0a9/71ec280a93767e11e76607e4fd40bc96.png)
С этими функциональными формами потенциала нахождение орбитального пути требует вычисления градиента потенциала, чтобы получить вектор гравитационного поля, с применением затем третьего закона Ньютона. Я получаю орбитальные уравнения движения для двух моделей гравитации:
![](https://habrastorage.org/getpro/habr/post_images/416/9d7/c17/4169d7c1726200f1586ba071c7a59fb0.png)
![](https://habrastorage.org/getpro/habr/post_images/aca/75a/c66/aca75ac662d4c46a33e7f26a22092927.png)
Теперь я готов использовать силу NDSolve для вычисления орбитальных траекторий. Однако перед тем, как сделать это, будет неплохо отобразить орбитальный путь в виде кривой в трехмерном пространстве. Чтобы встроить эти кривые в контекст, я построил их по текстурной карте поверхности Земли и спроецировал на сферу. Здесь я построил нужные графические объекты:
![](https://habrastorage.org/getpro/habr/post_images/1d2/586/96f/1d258696f1b5e135134c4c665f13f5a3.png)
![](https://habrastorage.org/getpro/habr/post_images/16e/f76/6d2/16ef766d2859570eefc8c2213b8f3c4b.png)
![](https://habrastorage.org/getpro/habr/post_images/f04/e24/0f8/f04e240f8ad7ac322e623c1e34de85d4.png)
![](https://habrastorage.org/getpro/habr/post_images/913/5e3/4de/9135e34de77c947017f885a3b9438a85.png)
В то время, как орбитальный путь, вычисляемый в инерциальной системе отсчета, образует периодическую замкнутую кривую, получается, что космический корабль будет проходить через различные точки на поверхности Земли во время каждого последующего оборота. Я могу визуализировать этот эффект, добавив дополнительное время вращения к решениям, которые я получаю от функции NDSolve. Примем, что число орбитальных периодов будет три (по аналогии с полетом Джона Гленна), я использовал функцию Manipulate чтобы увидеть, как орбитальный путь зависит от азимутального угла ?, подобно исследованию в статье Джонсон. Я построю один график для сферической Земли (белый цвет) и другой, зеленый, чтобы продемонстрировать эффект сплюснутости (обратите внимание, что расхождение увеличивается с каждым витком):
![](https://habrastorage.org/getpro/habr/post_images/36a/104/120/36a1041208bd2ebc2c2066bd6d1c48b7.png)
![](https://habrastorage.org/getpro/habr/post_images/f87/328/f61/f87328f611eee374127e571dafce3578.png)
В документе, который прилагается к этой записи, вы сможете увидеть функцию Manipulate в действии; обратите внимание на скорость, с которой находится каждое новое решение. Хочется надеяться, что Кэтрин Джонсон и ее коллеги в НАСА были бы впечатлены!
Теперь, изменяя угол ? в точке выгорания, легко вычислить положение космического аппарата после, скажем, трех оборотов:
![](https://habrastorage.org/getpro/habr/post_images/a39/7a2/c22/a397a2c224e2db1e920be87a102ddc86.png)
![](https://habrastorage.org/getpro/habr/post_images/05a/7cf/a48/05a7cfa481256f182fb893c1b05f6ce5.png)
![](https://habrastorage.org/getpro/habr/post_images/ac4/ee5/c83/ac4ee5c839c4269579231eb9b6dd64d5.png)
Моделирование возвращаемого спутника
В фильме также в связи с фазой возвращения упоминается метод Эйлера. После решения исходной задачи нахождения азимутального угла, как это было сделано в предыдущих разделах, пришло время вернуться на Землю. Для замедления вращения происходит запуск ракет, и затем, при переходе из безвоздушной среды космоса в атмосферу, происходит сложный комплекс событий. Для благополучного возвращения космонавта на Землю долны быть приняты во внимание такие факторы, как изменение плотности атмосферы, быстрое торможение и фрикционный нагрев. Необходимо решить проблемы высоты, скорости и ускорения как функции времени. Этот набор проблем может быть решен с помощью метода Эйлера, как это было сделано Кэтрин Джонсон, или с помощью функциональных возможностей решения дифференциальных уравнений в Wolfram Language.
Для простых дифференциальных уравнений можно получить подробное пошаговое решение с помощью указанного метода квадратур. Эквивалентом известной формулы Ньютона F = mа для зависящей от времени массы m(t) является так называемое уравнение для идеальной ракеты (в одном измерении)…
![](https://habrastorage.org/getpro/habr/post_images/ca2/cce/7b8/ca2cce7b8c3bea0d4fa1d38d660bde2a.png)
… где m(t) — масса ракеты, vе скорость выхлопа ракетного двигателя, и m'p(t) — производная массы топлива от времени. При постоянной m'p(t) структура уравнения оказывается относительно простой и легко разрешимой в замкнутой форме:
![](https://habrastorage.org/getpro/habr/post_images/e56/9bf/a59/e569bfa5948fd73a068d7e6c5ccf9cab.png)
![](https://habrastorage.org/getpro/habr/post_images/44d/f3b/a1f/44df3ba1f91012457b2d5f58c855374e.png)
Имея начальные и конечные условия для массы, я получаю знаменитое уравнение движения ракеты ( К.Э.Циолковский, 1903):
![](https://habrastorage.org/getpro/habr/post_images/fd2/f8b/32f/fd2f8b32f5d08488a121819f2c9c37af.png)
![](https://habrastorage.org/getpro/habr/post_images/272/0ef/5bb/2720ef5bb8e097023d0b59960756f6cf.png)
![](https://habrastorage.org/getpro/habr/post_images/dd3/95b/dce/dd395bdce9dd08c1da7415bd36e8df8e.png)
![](https://habrastorage.org/getpro/habr/post_images/b7a/e1e/b38/b7ae1eb3804423f25ff356ad7067a893.png)
![](https://habrastorage.org/getpro/habr/post_images/bff/218/1dc/bff2181dc86ead34f0d8f0c65641162d.png)
![](https://habrastorage.org/getpro/habr/post_images/fa1/daa/27c/fa1daa27c4c3d93f192a0aeed2b581c3.png)
![](https://habrastorage.org/getpro/habr/post_images/be9/ba1/84b/be9ba184b9e3c89dba96565b6d45c931.png)
![](https://habrastorage.org/getpro/habr/post_images/db1/ff5/71b/db1ff571b6be1d1b66037a13f2fd5c50.png)
Детали решения этого уравнения с конкретными значениями параметров я могу взять из Wolfram|Alpha. Вот эти детали вместе со сравнением с точным решением, а также с другими методами численного интегрирования:
![](https://habrastorage.org/getpro/habr/post_images/70f/5ea/929/70f5ea9298c972a4ee4c7d4d6bb0afb6.png)
![](https://habrastorage.org/getpro/habr/post_images/615/4bc/9a0/6154bc9a03298c2501c2da01be0d3fbc.png)
Следуя сюжету фильма, я реализую теперь минималистическую модель процесса повторного входа. Начнем с определения параметров, которые имитируют полет Гленна:
![](https://habrastorage.org/getpro/habr/post_images/3a1/7fe/e4b/3a17fee4bc485db9f549972714f487d2.png)
![](https://habrastorage.org/getpro/habr/post_images/b5f/d34/07c/b5fd3407cd6328b71829ae47426198e2.png)
Я полагаю, что в процессе торможения используется 1% от тяги двигателя первой ступени, и длится он, скажем, 60 секунд. Уравнение движения:
![](https://habrastorage.org/getpro/habr/post_images/4b8/3cf/44c/4b83cf44ceea205c6454d4125cf5fef8.png)
Здесь Fgrav — сила тяжести, Fexhaust(t) — зависящая от времени сила двигателя, и сила трения Ffriction(x(t),v(t)). Последняя зависит как от плотности воздуха в положении x(t), так и от закона трения для v(t).
Для плотности воздуха, которая зависит от высоты, я могу использовать функцию StandardAtmosphereData. Я учитываю ее также из-за парашюта, который открывается на высоте около 8,5 км над землей:
![](https://habrastorage.org/getpro/habr/post_images/004/c53/4f4/004c534f42e503f95207ce3a126df59f.png)
Это дает следующий набор связанных нелинейных дифференциальных уравнений, которые необходимо решить. Функция WhenEvent[...] указывает на конец построения решения ДУ, когда капсула достигает поверхности Земли. Я использую векторные значения положения и скорости переменных X и V:
![](https://habrastorage.org/getpro/habr/post_images/5ca/e8b/045/5cae8b045306035bc4fc3804b7879d7e.png)
С учетом этих определений для веса, выхлопных газов и силы трения воздуха точки зрения…
![](https://habrastorage.org/getpro/habr/post_images/e26/bec/a5d/e26beca5dd8e29ad41a68df8e8ffbba8.png)
… результирующая сила может быть найдена с помощью:
![](https://habrastorage.org/getpro/habr/post_images/a9d/f00/439/a9df004398d767a6ebca5645b8226d84.png)
![](https://habrastorage.org/getpro/habr/post_images/fbb/e44/195/fbbe441951aec356e6f82161cfc43eca.png)
В этой модели я пренебрег вращением Земли, внутренними углами вращения капсулы, активными изменениями угла полета, сверхзвуковыми эффектами силы трения и многим другим. Уравнения, решаемые Кэтрин Джонсон:
![](https://habrastorage.org/getpro/habr/post_images/6f5/0f2/3ec/6f50f23ec50560d15e374b1b33a839db.png)
![](https://habrastorage.org/getpro/habr/post_images/9c9/bd7/6ba/9c9bd76ba85adbcca58d77a8e0d85040.png)
Эту систему уравнений просто решить численно, если дополнить ее начальным положением и скоростью. Для этого нужно всего лишь обратиться к функции NDSolve. Мне не придется беспокоиться об используемом методе, управлении размером шага, контроле ошибок и так далее, потому что Wolfram Language автоматически выбирает те значения, которые гарантируют значимые результаты:
![](https://habrastorage.org/getpro/habr/post_images/6f9/bbc/e73/6f9bbce731f9d0e90ebb38f7b8993e3c.png)
![](https://habrastorage.org/getpro/habr/post_images/b92/158/238/b921582382f6086eeea5bf0fe0012e17.png)
Здесь представлен график высоты, скорости и ускорения в зависимости от времени:
![](https://habrastorage.org/getpro/habr/post_images/9e4/636/4a8/9e46364a89dc85aef9eec96909409e7b.png)
![](https://habrastorage.org/getpro/habr/post_images/d65/4ed/c42/d654edc42428c59892d82e829e5ec83a.png)
![](https://habrastorage.org/getpro/habr/post_images/63b/933/4b4/63b9334b46937bbdfb6e01fb4a592a27.png)
Кривая зависимости от высоты вместо времени демонстрирует, что экспоненциальное возрастание плотности воздуха отвечает за сильное торможение. Это не связано с парашютированием, которое происходит на относительно малой высоте. Пик замедления происходит на очень большой высоте, когда капсула быстро переходит из вакуума к атмосферной среде:
![](https://habrastorage.org/getpro/habr/post_images/8b2/86d/379/8b286d379e5a8cdd2b51255527ade5bf.png)
![](https://habrastorage.org/getpro/habr/post_images/118/1fe/0d1/1181fe0d1816ab6aee50f7cbc757f906.png)
А вот график вертикальной и тангенциальной скоростей капсулы в процессе повторного входа:
![](https://habrastorage.org/getpro/habr/post_images/8cb/55a/db5/8cb55adb5310a09ac7f15ffe25212bba.png)
![](https://habrastorage.org/getpro/habr/post_images/ec0/9af/a48/ec09afa4831c6e64a235b5f774b73011.png)
Теперь я повторю численное решение методом Эйлера с фиксированным шагом:
![](https://habrastorage.org/getpro/habr/post_images/5a7/06e/c98/5a706ec98a04ef88adb1d01e9270f323.png)
![](https://habrastorage.org/getpro/habr/post_images/b17/275/f33/b17275f331e6ca71adf4928e626fba18.png)
Качественно это решение выглядит так же, как и предыдущее:
![](https://habrastorage.org/getpro/habr/post_images/fbc/797/868/fbc797868f34dc79db35d445c05070f5.png)
![](https://habrastorage.org/getpro/habr/post_images/29c/033/97d/29c03397dccfb626ee1ea8b5e6d36ac3.png)
Для используемого размера шага интегрирования по времени накопленная ошибка составляет порядка нескольких процентов. Меньшие размеры шага будут уменьшать ошибку:
![](https://habrastorage.org/getpro/habr/post_images/a0e/1fe/aea/a0e1feaea80e93b2c71da81461c4255e.png)
![](https://habrastorage.org/getpro/habr/post_images/19b/865/b6e/19b865b6e661c0942b091e1eb9c4600e.png)
Обратите внимание, что время посадки, предсказываемое методом Эйлера, отклоняется всего на 0,11% по сравнению с предыдущим временем (для сравнения: если бы я решал уравнение двумя современными методами, — скажем, «BDF» или «Adams», то ошибка была бы меньше на несколько порядков).
В процессе повторного входа генерируется большое количество тепла. Здесь нужен теплозащитный экран. На какой высоте производится наибольшее количество тепла на единицу площади Q? Не вдаваясь в подробности, я могу, чисто из соображений размерности, применить гипотезу
![](https://habrastorage.org/getpro/habr/post_images/0b9/0f4/c0a/0b90f4c0a7bdc010226db04b5cdc89c6.png)
![](https://habrastorage.org/getpro/habr/post_images/e82/67f/231/e8267f231d2f2120926e4e100d47f66b.png)
![](https://habrastorage.org/getpro/habr/post_images/498/326/0fb/4983260fb5ca2123ee1663958ebe7937.png)
![](https://habrastorage.org/getpro/habr/post_images/18c/98f/62a/18c98f62a788be5d2465ecd88c5dc2e7.png)
![](https://habrastorage.org/getpro/habr/post_images/047/b7f/f0d/047b7ff0d58a88e5817593c31b13f5b3.png)
Много чего интересного можно было бы вычислить (Hicks 2009), но точно так же, как фильм должен был закончиться, так и я теперь должен завершить мою статью. Надеюсь, вы простите меня за мои слова: с Wolfram Language вы можете сделать все, что захотите.
Поделиться с друзьями
Комментарии (7)
jaiprakash
19.04.2017 08:42«Burnout» в данном контексте соответствует термину «про?жиг», имхо.
OsipovRoman
19.04.2017 09:15Тут ничего не прожигается. Двигатель прекращает работу в результате выгорания топлива.
jaiprakash
19.04.2017 12:03+1Ну если момент выключения, то «отсечка».
Буквальное соответствие выгорания и отсечки выполняется только для ТТРД, для ЖРД отсечка выполняется по команде, топливо при этом остаётся. А так невольно получается привязка именно к ТТРД. Сейчас и для них пишут «Buster Engine Cut-Off».
Извините за буквоедство, но глаз режет когда вместо русского термина стоит калька с устаревшего английского.
Mysterious
Вот это я понимаю статья для Хабра!
Большое спасибо вам за ваш труд!
mephistopheies
в первую очередь спасибо авторам, а уже затем переводчику, т.к. в этом блоге публикуются только переводы чужих трудов (ну или почти)
OsipovRoman
Мы переводим только посты авторов из нашей компании. Это наши посты из нашего официального блога и других наших англоязычных ресурсов. Так что это не «чужие труды». Также в блоге есть, пока немного, постов, написанных авторами из России.
OsipovRoman
Спасибо, что читаете наш блог! В скором времени выйдет еще много интересных постов.