- PVSM.RU - https://www.pvsm.ru -
Перевод поста Джеффри Брайанта (Jeffrey Bryant), Пако Джейна (Paco Jain) и Майкл Тротта (Michael Trott) "Hidden Figures: Modern Approaches to Orbit and Reentry Calculations [1]".
Код, приведенный в статье, можно скачать здесь [2].
Выражаю огромную благодарность Полине Сологуб [3] за помощь в переводе и подготовке публикации
— Размещение спутника в определенном месте [4]
— Константы и первичная обработка [5]
— Вычисления [6]
— Построение графика [7]
— Как рассчитываются орбиты сегодня [8]
Моделирование возвращаемого спутника [9]
Я остановлюсь на двух аспектах ее научной работы, упомянутых в фильме: вычислениях орбиты и расчетах, связанных с вхождением в атмосферу. Для орбитальных вычислений я сначала сделал ровно то же, что и Джонсон, а затем применил более современный прямой подход с использованием инструментов Wolfram Language [14]. В фильме упоминается о решении дифференциальных уравнений методом Эйлера [15]; я же буду сравнивать этот метод с более современным и вычислю возвратную траекторию с помощью данных модели атмосферы, полученных непосредственно из Wolfram Language).
В фильме не освещаются детали решения математических задач командой Джонсон, однако того, что вы увидите, хватит, чтобы прочувствовать и сравнить подход, отображенный в фильме, с тем, который существует на данный момент.
Одна из самых ранних работ, в которой Джонсон была соавтором ("Определение азимутального угла в точке выгорания для размещения спутника над выбранной точкой [16]") имеет дело с проблемой проверки, может ли спутник быть размещен над определенной точкой Земли после определенного количества витков вокруг орбиты при определенной стартовой позиции (например, мыс Канаверал, штат Флорида) и траектории вращения. Подход, который команда Джонсон использовала для определения азимутального угла (угол, образованный вектором скорости космического аппарата в момент отключения двигателя, с фиксированным направлением, скажем, на север) для запуска ракеты основывался на других орбитальных параметрах. Это важный шаг для обеспечения правильного местонахождения астронавта для входа в атмосферу на Землю.
В статье Джонсон определяет ряд констант и входных параметров, необходимых для решения задачи вручную. Отступлю немного, чтобы объяснить термин «выгорание», относящийся к остановке ракетного двигателя. После «выгорания» параметры орбиты как бы «замораживаются», и космический аппарат движется только под действием силы тяжести Земли (в соответствии с законами Ньютона). В этом разделе я следую тем определениям единиц, которые были приняты в статье.
Некоторые функции для удобства определены, чтобы взаимодействовать с углами в градусах вместо радиан. Это позволяет плавно регулировать время вычисления углов:
Джонсон описывает несколько других производных параметров, хотя интересно отметить, что иногда она принимает определенные значения для них, а не использует значения, полученные с помощью ее формул. Принятые ею значения часто были близки к значениям, полученным по формулам. Для простоты здесь используются значения из формул.
Semilatus rectum [17] эллиптической орбиты:
Угол в плоскости орбиты между перигеем и точкой выгорания:
Эксцентриситет орбиты:
Период орбиты:
Эксцентрическая аномалия:
Для описания следующего параметра проще будет процитировать оригинальную статью: "требование о том, чтобы спутник с точкой выгорания φ1, λ1 проходил через выбранную позицию φ2, λ2 после завершения n прохождений по орбите эквивалентно требованию, в соответствии с которым в ходе первого оборота спутник проходит через эквивалентное место с широтой φ2 так же, как и при выбранном положении, но с долготой λ2, смещенной в восточном направлении от λ2 на величину, достаточную для компенсации вращения Земли в течение n полных вращений, то есть полярного часового угла nωЕТ. Долгота этой эквивалентной позиции, таким образом, определяется соотношением":
Время от перигея для угла θ:
Часть окончательного решения заключается в определении значения для промежуточных параметров δλ1-2е и θ2е. Это может быть сделано несколькими способами. Во-первых, я могу использовать ContourPlot [18] для получения графического решения с помощью уравнений 19 и 20 из статьи:
FindRoot [19] можно использовать для нахождения численных решений:
Конечно, у Джонсон не было доступа к функциям ContourPlot или FindRoot, так что в ее статье описывается итерационный метод. Я перевел методику, описанную в статье, на язык Wolfram, а также раскрыл ряд других параметров с помощью ее итерационного метода. Поскольку базовые выкладки рассчитаны для сферической формы Земли, поправки на сплющенности включены:
Графическое изображение значения θ2е для различных итераций демонстрирует быструю сходимость:
Я могу преобразовать метод в команде FindRoot следующим образом:
Интересно, что даже итерационные шаги поиска корня этой более сложной системы сходятся довольно быстро:
Когда орбитальные параметры определены, можно визуализировать решение. Во-первых, из предыдущих решений должны быть извлечены некоторые параметры:
Далее должны быть получены широта и долгота спутника в зависимости от азимутального угла:
φs и λs — широта и долгота как функция θs:
Наземный трек спутника может быть построен путем создания таблицы точек:
В статье Джонсон представлена схема орбитального решения, включая маркеры, отмечающие выгорание, а также выбранные и эквивалентные положения. Подобную схему легко воспроизвести:
Вот ее оригинальная схема для сравнения:
Более визуально полезная версия может быть построена с помощью функции GeoGraphics [20], которая преобразует геоцентрические координаты в геодезические:
Сегодня практически у каждого из нас есть доступ к вычислительным ресурсам гораздо более мощным, чем те, что были у НАСА в 1960-е годы. Сейчас, используя только настольный компьютер и Wolfram Language, вы можете легко найти прямые численные решения задач орбитальной механики вроде тех, которыми занимались Кэтрин Джонсон и ее команда.
Чтобы вычислить азимутальный угол ψ с помощью более современных методов, давайте установим параметры для простой круговой орбиты, берущей свое начало после выгорания над Флоридой; предположим также, что Земля сферически симметрична (я не буду пытаться искать точное соответствие с данными из статьи Джонсон и переопределю некоторые величины, используя современную систему единиц СИ). Начиная с той же высоты околоземной орбиты, используемой Джонсон, и с помощью сферической тригонометрии несложно вывести начальные условия для нашей орбиты:
Соответствующие физические параметры могут быть получены непосредственно с помощью Wolfram Language:
Далее я получил дифференциальное уравнение движения нашего космического корабля с учетом гравитационного поля Земли. Есть несколько способов моделирования гравитационного потенциала вблизи Земли. Предположим, что планета сферически симметрична и используем декартову систему координат:
Кроме того, вы можете использовать более реалистичную модель гравитации Земли, где в качестве формы планеты берется сплюснутый эллипсоид вращения. Точная форма потенциала такого эллипсоида (при условии постоянной плотности массы эллипсоидальных оболочек) доступна с помощью EntityValue [21]:
Для общего однородного трехосного эллипсоида потенциал содержит кусочные функции:
Здесь κ является наибольшим корнем выражения х2/(a2 + κ) + у2/(b2 + k) + z2/(c2 + κ) = 1. В случае, если мы имеем дело со сплющенным эллипсоидом, предыдущая формула может быть упрощена до элементарных функций…
… где κ=((2 z2 (a2-c2+x2+y2)+(-a2+c2+x2+y2)2+z4)1/2-a2-c2+x2+y2+z2)/2.
Более простая форма, которую я буду использовать здесь, задается так называемой общепринятой формулой ускорения силы тяжести (IGF). В данной формуле учитываются отличия от сферически симметричного потенциала до второго порядка в сферических гармониках, и даются численно неотличимые от точного потенциала, на который ссылались ранее, результаты. В терминах четырех измеренных геодезических параметров потенциал IGF может быть определен следующим образом:
Я мог бы легко использовать даже еще более точные значения для силы тяжести с помощью функции GeogravityModelData [22]. В исходном положении потенциал IGF только на 0,06% отклоняется от приближения высокого порядка:
С этими функциональными формами потенциала нахождение орбитального пути требует вычисления градиента потенциала, чтобы получить вектор гравитационного поля, с применением затем третьего закона Ньютона. Я получаю орбитальные уравнения движения для двух моделей гравитации:
Теперь я готов использовать силу NDSolve [23] для вычисления орбитальных траекторий. Однако перед тем, как сделать это, будет неплохо отобразить орбитальный путь в виде кривой в трехмерном пространстве. Чтобы встроить эти кривые в контекст, я построил их по текстурной карте поверхности Земли и спроецировал на сферу. Здесь я построил нужные графические объекты:
В то время, как орбитальный путь, вычисляемый в инерциальной системе отсчета, образует периодическую замкнутую кривую, получается, что космический корабль будет проходить через различные точки на поверхности Земли во время каждого последующего оборота. Я могу визуализировать этот эффект, добавив дополнительное время вращения к решениям, которые я получаю от функции NDSolve. Примем, что число орбитальных периодов будет три (по аналогии с полетом Джона Гленна), я использовал функцию Manipulate [24] чтобы увидеть, как орбитальный путь зависит от азимутального угла ψ, подобно исследованию в статье Джонсон. Я построю один график для сферической Земли (белый цвет) и другой, зеленый, чтобы продемонстрировать эффект сплюснутости (обратите внимание, что расхождение увеличивается с каждым витком):
В документе, который прилагается к этой записи, вы сможете увидеть функцию Manipulate в действии; обратите внимание на скорость, с которой находится каждое новое решение. Хочется надеяться, что Кэтрин Джонсон и ее коллеги в НАСА были бы впечатлены!
Теперь, изменяя угол ψ в точке выгорания, легко вычислить положение космического аппарата после, скажем, трех оборотов:
В фильме также в связи с фазой возвращения упоминается метод Эйлера. После решения исходной задачи нахождения азимутального угла, как это было сделано в предыдущих разделах, пришло время вернуться на Землю. Для замедления вращения происходит запуск ракет, и затем, при переходе из безвоздушной среды космоса в атмосферу, происходит сложный комплекс событий. Для благополучного возвращения космонавта на Землю долны быть приняты во внимание такие факторы, как изменение плотности атмосферы, быстрое торможение и фрикционный нагрев. Необходимо решить проблемы высоты, скорости и ускорения как функции времени. Этот набор проблем может быть решен с помощью метода Эйлера, как это было сделано Кэтрин Джонсон, или с помощью функциональных возможностей решения дифференциальных уравнений в Wolfram Language.
Для простых дифференциальных уравнений можно получить подробное пошаговое решение с помощью указанного метода квадратур. Эквивалентом известной формулы Ньютона F = mа для зависящей от времени массы m(t) является так называемое уравнение для идеальной ракеты (в одном измерении)…
… где m(t) — масса ракеты, vе скорость выхлопа ракетного двигателя, и m'p(t) — производная массы топлива от времени. При постоянной m'p(t) структура уравнения оказывается относительно простой и легко разрешимой в замкнутой форме:
Имея начальные и конечные условия для массы, я получаю знаменитое уравнение движения ракеты ( К.Э.Циолковский, 1903 [25]):
Детали решения этого уравнения с конкретными значениями параметров я могу взять из Wolfram|Alpha [26]. Вот эти детали вместе со сравнением с точным решением, а также с другими методами численного интегрирования:
Следуя сюжету фильма, я реализую теперь минималистическую модель процесса повторного входа. Начнем с определения параметров, которые имитируют полет Гленна:
Я полагаю, что в процессе торможения используется 1% от тяги двигателя первой ступени, и длится он, скажем, 60 секунд. Уравнение движения:
Здесь Fgrav — сила тяжести, Fexhaust(t) — зависящая от времени сила двигателя, и сила трения Ffriction(x(t),v(t)). Последняя зависит как от плотности воздуха в положении x(t), так и от закона трения для v(t).
Для плотности воздуха, которая зависит от высоты, я могу использовать функцию StandardAtmosphereData [27]. Я учитываю ее также из-за парашюта, который открывается на высоте около 8,5 км над землей:
Это дает следующий набор связанных нелинейных дифференциальных уравнений, которые необходимо решить. Функция WhenEvent [28][...] указывает на конец построения решения ДУ, когда капсула достигает поверхности Земли. Я использую векторные значения положения и скорости переменных X и V:
С учетом этих определений для веса, выхлопных газов и силы трения воздуха точки зрения…
… результирующая сила может быть найдена с помощью:
В этой модели я пренебрег вращением Земли, внутренними углами вращения капсулы, активными изменениями угла полета, сверхзвуковыми эффектами силы трения и многим другим. Уравнения, решаемые Кэтрин Джонсон:
Эту систему уравнений просто решить численно, если дополнить ее начальным положением и скоростью. Для этого нужно всего лишь обратиться к функции NDSolve. Мне не придется беспокоиться об используемом методе, управлении размером шага, контроле ошибок и так далее, потому что Wolfram Language автоматически выбирает те значения, которые гарантируют значимые результаты:
Здесь представлен график высоты, скорости и ускорения в зависимости от времени:
Кривая зависимости от высоты вместо времени демонстрирует, что экспоненциальное возрастание плотности воздуха отвечает за сильное торможение. Это не связано с парашютированием, которое происходит на относительно малой высоте. Пик замедления происходит на очень большой высоте, когда капсула быстро переходит из вакуума к атмосферной среде:
А вот график вертикальной и тангенциальной скоростей капсулы в процессе повторного входа:
Теперь я повторю численное решение методом Эйлера с фиксированным шагом:
Качественно это решение выглядит так же, как и предыдущее:
Для используемого размера шага интегрирования по времени накопленная ошибка составляет порядка нескольких процентов. Меньшие размеры шага будут уменьшать ошибку:
Обратите внимание, что время посадки, предсказываемое методом Эйлера, отклоняется всего на 0,11% по сравнению с предыдущим временем (для сравнения: если бы я решал уравнение двумя современными методами, — скажем, «BDF» или «Adams», то ошибка была бы меньше на несколько порядков).
В процессе повторного входа генерируется большое количество тепла. Здесь нужен теплозащитный экран. На какой высоте производится наибольшее количество тепла на единицу площади Q? Не вдаваясь в подробности, я могу, чисто из соображений размерности, применить гипотезу :
Много чего интересного можно было бы вычислить (Hicks 2009 [29]), но точно так же, как фильм должен был закончиться, так и я теперь должен завершить мою статью. Надеюсь, вы простите меня за мои слова: с Wolfram Language вы можете сделать все, что захотите.
Автор: Wolfram Research
Источник [30]
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/programmirovanie/253023
Ссылки в тексте:
[1] Hidden Figures: Modern Approaches to Orbit and Reentry Calculations: http://blog.wolfram.com/2017/02/24/hidden-figures-modern-approaches-to-orbit-and-reentry-calculations/
[2] здесь: http://blog.wolfram.com/data/uploads/2017/02/Hidden-figures-modern-approaches-to-orbit-and-reentry-calculations.cdf
[3] Полине Сологуб: https://vk.com/id65091763
[4] Размещение спутника в определенном месте: https://habrahabr.ru/company/wolfram/blog/NUMBER/#1
[5] Константы и первичная обработка: https://habrahabr.ru/company/wolfram/blog/NUMBER/#2
[6] Вычисления: https://habrahabr.ru/company/wolfram/blog/NUMBER/#3
[7] Построение графика: https://habrahabr.ru/company/wolfram/blog/NUMBER/#4
[8] Как рассчитываются орбиты сегодня: https://habrahabr.ru/company/wolfram/blog/NUMBER/#5
[9] Моделирование возвращаемого спутника: https://habrahabr.ru/company/wolfram/blog/NUMBER/#6
[10] Скрытые фигуры: http://www.wolframalpha.com/input/?i=Hidden+Figures
[11] Кэтрин Джонсон: https://ru.wikipedia.org/wiki/%D0%94%D0%B6%D0%BE%D0%BD%D1%81%D0%BE%D0%BD,_%D0%9A%D1%8D%D1%82%D1%80%D0%B8%D0%BD
[12] Дороти Воган: https://en.wikipedia.org/wiki/Dorothy_Vaughan
[13] Мэри Джексон: https://en.wikipedia.org/wiki/Mary_Jackson_(engineer)
[14] Wolfram Language: http://www.wolfram.com/language/
[15] методом Эйлера: http://mathworld.wolfram.com/EulerForwardMethod.html
[16] Определение азимутального угла в точке выгорания для размещения спутника над выбранной точкой: https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980227091.pdf
[17] Semilatus rectum: http://mathworld.wolfram.com/SemilatusRectum.html
[18] ContourPlot: http://reference.wolfram.com/language/ref/ContourPlot.html
[19] FindRoot: http://reference.wolfram.com/language/ref/FindRoot.html
[20] GeoGraphics: http://reference.wolfram.com/language/ref/GeoGraphics.html
[21] EntityValue: http://reference.wolfram.com/language/ref/EntityValue.html
[22] GeogravityModelData: http://reference.wolfram.com/language/ref/GeogravityModelData.html
[23] NDSolve: http://reference.wolfram.com/language/ref/NDSolve.html
[24] Manipulate: http://reference.wolfram.com/language/ref/Manipulate.html
[25] К.Э.Циолковский, 1903: http://www.tsiolkovsky.ru/index.php?option=com_content&task=view&id=64
[26] Wolfram|Alpha: http://www.wolframalpha.com/
[27] StandardAtmosphereData: http://reference.wolfram.com/language/ref/StandardAtmosphereData.html
[28] WhenEvent: http://reference.wolfram.com/language/ref/WhenEvent.html
[29] Hicks 2009: http://www.dtic.mil/dtic/tr/fulltext/u2/a505342.pdf
[30] Источник: https://habrahabr.ru/post/326756/?utm_source=habrahabr&utm_medium=rss&utm_campaign=best
Нажмите здесь для печати.