Перевод поста Джеффри Брайанта (Jeffrey Bryant), Пако Джейна (Paco Jain) и Майкл Тротта (Michael Trott) "Hidden Figures: Modern Approaches to Orbit and Reentry Calculations".
Код, приведенный в статье, можно скачать здесь.
Выражаю огромную благодарность Полине Сологуб за помощь в переводе и подготовке публикации

---

Содержание

— Размещение спутника в определенном месте
— Константы и первичная обработка
— Вычисления
— Построение графика
— Как рассчитываются орбиты сегодня
Моделирование возвращаемого спутника

---

Вышедший недавно в кинотеатрах фильм Скрытые фигуры получил хорошие отзывы. Действие разворачивается в важный период истории США; в нем затрагивается также ряд тем вроде гражданских прав и космической гонки. В центре повествования — история Кэтрин Джонсон и ее коллег (Дороти Воган и Мэри Джексон) из NASA в период развертывания программы Меркурий и ранних исследований пилотируемых космических полетов. Внимание также акцентируется на драматической борьбе за гражданские права афро-американских женщин в NASA, происходившей в то время. Компьютеры в то время едва появились, так что способность Джонсон и ее коллег решать сложные навигационные задачи орбитальной механики без использования компьютера обеспечили важную проверку ранних компьютерных результатов.



Я остановлюсь на двух аспектах ее научной работы, упомянутых в фильме: вычислениях орбиты и расчетах, связанных с вхождением в атмосферу. Для орбитальных вычислений я сначала сделал ровно то же, что и Джонсон, а затем применил более современный прямой подход с использованием инструментов Wolfram Language. В фильме упоминается о решении дифференциальных уравнений методом Эйлера; я же буду сравнивать этот метод с более современным и вычислю возвратную траекторию с помощью данных модели атмосферы, полученных непосредственно из Wolfram Language).

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

Размещение спутника в определенном месте

Одна из самых ранних работ, в которой Джонсон была соавтором ("Определение азимутального угла в точке выгорания для размещения спутника над выбранной точкой") имеет дело с проблемой проверки, может ли спутник быть размещен над определенной точкой Земли после определенного количества витков вокруг орбиты при определенной стартовой позиции (например, мыс Канаверал, штат Флорида) и траектории вращения. Подход, который команда Джонсон использовала для определения азимутального угла (угол, образованный вектором скорости космического аппарата в момент отключения двигателя, с фиксированным направлением, скажем, на север) для запуска ракеты основывался на других орбитальных параметрах. Это важный шаг для обеспечения правильного местонахождения астронавта для входа в атмосферу на Землю.

Константы и первичная обработка

В статье Джонсон определяет ряд констант и входных параметров, необходимых для решения задачи вручную. Отступлю немного, чтобы объяснить термин «выгорание», относящийся к остановке ракетного двигателя. После «выгорания» параметры орбиты как бы «замораживаются», и космический аппарат движется только под действием силы тяжести Земли (в соответствии с законами Ньютона). В этом разделе я следую тем определениям единиц, которые были приняты в статье.



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





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

Semilatus rectum эллиптической орбиты:



Угол в плоскости орбиты между перигеем и точкой выгорания:



Эксцентриситет орбиты:



Период орбиты:



Эксцентрическая аномалия:



Для описания следующего параметра проще будет процитировать оригинальную статью: "требование о том, чтобы спутник с точкой выгорания ?1, ?1 проходил через выбранную позицию ?2, ?2 после завершения n прохождений по орбите эквивалентно требованию, в соответствии с которым в ходе первого оборота спутник проходит через эквивалентное место с широтой ?2 так же, как и при выбранном положении, но с долготой ?2, смещенной в восточном направлении от ?2 на величину, достаточную для компенсации вращения Земли в течение n полных вращений, то есть полярного часового угла n?_ЕТ. Долгота этой эквивалентной позиции, таким образом, определяется соотношением":



Время от перигея для угла ?:

Вычисления

Часть окончательного решения заключается в определении значения для промежуточных параметров ??_1-2е и ?_2е. Это может быть сделано несколькими способами. Во-первых, я могу использовать ContourPlot для получения графического решения с помощью уравнений 19 и 20 из статьи:



FindRoot можно использовать для нахождения численных решений:





Конечно, у Джонсон не было доступа к функциям ContourPlot или FindRoot, так что в ее статье описывается итерационный метод. Я перевел методику, описанную в статье, на язык Wolfram, а также раскрыл ряд других параметров с помощью ее итерационного метода. Поскольку базовые выкладки рассчитаны для сферической формы Земли, поправки на сплющенности включены:



Графическое изображение значения ?2е для различных итераций демонстрирует быструю сходимость:









Я могу преобразовать метод в команде FindRoot следующим образом:







Интересно, что даже итерационные шаги поиска корня этой более сложной системы сходятся довольно быстро:

Построение графика

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



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





?s и ?s — широта и долгота как функция ?s:





Наземный трек спутника может быть построен путем создания таблицы точек:



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







Вот ее оригинальная схема для сравнения:



Более визуально полезная версия может быть построена с помощью функции GeoGraphics, которая преобразует геоцентрические координаты в геодезические:

Как рассчитываются орбиты сегодня

Сегодня практически у каждого из нас есть доступ к вычислительным ресурсам гораздо более мощным, чем те, что были у НАСА в 1960-е годы. Сейчас, используя только настольный компьютер и Wolfram Language, вы можете легко найти прямые численные решения задач орбитальной механики вроде тех, которыми занимались Кэтрин Джонсон и ее команда.

Чтобы вычислить азимутальный угол ? с помощью более современных методов, давайте установим параметры для простой круговой орбиты, берущей свое начало после выгорания над Флоридой; предположим также, что Земля сферически симметрична (я не буду пытаться искать точное соответствие с данными из статьи Джонсон и переопределю некоторые величины, используя современную систему единиц СИ). Начиная с той же высоты околоземной орбиты, используемой Джонсон, и с помощью сферической тригонометрии несложно вывести начальные условия для нашей орбиты:







Соответствующие физические параметры могут быть получены непосредственно с помощью Wolfram Language:



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



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



Для общего однородного трехосного эллипсоида потенциал содержит кусочные функции:





Здесь ? является наибольшим корнем выражения х²/(a² + ?) + у²/(b² + k) + z²/(c² + ?) = 1. В случае, если мы имеем дело со сплющенным эллипсоидом, предыдущая формула может быть упрощена до элементарных функций…





… где ?=((2 z² (a²-c²+x²+y²)+(-a²+c²+x²+y²)²+z⁴)^1/2-a²-c²+x²+y²+z²)/2.

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



Я мог бы легко использовать даже еще более точные значения для силы тяжести с помощью функции GeogravityModelData. В исходном положении потенциал IGF только на 0,06% отклоняется от приближения высокого порядка:





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





Теперь я готов использовать силу NDSolve для вычисления орбитальных траекторий. Однако перед тем, как сделать это, будет неплохо отобразить орбитальный путь в виде кривой в трехмерном пространстве. Чтобы встроить эти кривые в контекст, я построил их по текстурной карте поверхности Земли и спроецировал на сферу. Здесь я построил нужные графические объекты:









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





В документе, который прилагается к этой записи, вы сможете увидеть функцию Manipulate в действии; обратите внимание на скорость, с которой находится каждое новое решение. Хочется надеяться, что Кэтрин Джонсон и ее коллеги в НАСА были бы впечатлены!

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

Моделирование возвращаемого спутника

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

Для простых дифференциальных уравнений можно получить подробное пошаговое решение с помощью указанного метода квадратур. Эквивалентом известной формулы Ньютона F = mа для зависящей от времени массы m(t) является так называемое уравнение для идеальной ракеты (в одном измерении)…



… где m(t) — масса ракеты, v_е скорость выхлопа ракетного двигателя, и m^'_p(t) — производная массы топлива от времени. При постоянной m^'_p(t) структура уравнения оказывается относительно простой и легко разрешимой в замкнутой форме:





Имея начальные и конечные условия для массы, я получаю знаменитое уравнение движения ракеты ( К.Э.Циолковский, 1903):

















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





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





Я полагаю, что в процессе торможения используется 1% от тяги двигателя первой ступени, и длится он, скажем, 60 секунд. Уравнение движения:



Здесь F_grav — сила тяжести, F_exhaust(t) — зависящая от времени сила двигателя, и сила трения F_friction(x(t),v(t)). Последняя зависит как от плотности воздуха в положении x(t), так и от закона трения для v(t).

Для плотности воздуха, которая зависит от высоты, я могу использовать функцию StandardAtmosphereData. Я учитываю ее также из-за парашюта, который открывается на высоте около 8,5 км над землей:



Это дает следующий набор связанных нелинейных дифференциальных уравнений, которые необходимо решить. Функция WhenEvent[...] указывает на конец построения решения ДУ, когда капсула достигает поверхности Земли. Я использую векторные значения положения и скорости переменных X и V:



С учетом этих определений для веса, выхлопных газов и силы трения воздуха точки зрения…



… результирующая сила может быть найдена с помощью:





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





Эту систему уравнений просто решить численно, если дополнить ее начальным положением и скоростью. Для этого нужно всего лишь обратиться к функции NDSolve. Мне не придется беспокоиться об используемом методе, управлении размером шага, контроле ошибок и так далее, потому что Wolfram Language автоматически выбирает те значения, которые гарантируют значимые результаты:





Здесь представлен график высоты, скорости и ускорения в зависимости от времени:







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





А вот график вертикальной и тангенциальной скоростей капсулы в процессе повторного входа:





Теперь я повторю численное решение методом Эйлера с фиксированным шагом:





Качественно это решение выглядит так же, как и предыдущее:





Для используемого размера шага интегрирования по времени накопленная ошибка составляет порядка нескольких процентов. Меньшие размеры шага будут уменьшать ошибку:





Обратите внимание, что время посадки, предсказываемое методом Эйлера, отклоняется всего на 0,11% по сравнению с предыдущим временем (для сравнения: если бы я решал уравнение двумя современными методами, — скажем, «BDF» или «Adams», то ошибка была бы меньше на несколько порядков).

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









Много чего интересного можно было бы вычислить (Hicks 2009), но точно так же, как фильм должен был закончиться, так и я теперь должен завершить мою статью. Надеюсь, вы простите меня за мои слова: с Wolfram Language вы можете сделать все, что захотите.