Размещено на http://www.allbest.ru/

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ

«ВЫСШАЯ ШКОЛА ЭКОНОМИКИ»

Московский институт электроники и математики

Национального исследовательского университета

«Высшая школа экономики»

Департамент прикладной математики

Расчет траекторий перелета с низкой околоземной орбиты в окрестности коллинеарных точек либрации системы солнце-земля

Выпускная квалификационная работа - МАГИСТЕРСКАЯ ДИССЕРТАЦИЯ

по направлению 01.04.02 Прикладная математика и информатика

студента группы МА-22 (образовательная программа «Математические методы естествознания и компьютерные технологии»)

Федоренко Юлия Викторовна

Рецензент

Начальник отдела баллистического и навигационного проектирования НПО им. С.А.Лавочкина

А.В. Погодин

Научный руководитель

к.т.н., доцент департамента прикладной математики МИЭМ НИУ ВШЭ

С.А. Аксенов

Москва 2015



АННОТАЦИЯ

Работа посвящена исследованию возможностей одноимпульсного перехода с низкой околоземной орбиты высотой 500 км на квазипериодические орбиты вокруг точки либрации L₂ системы Солнце-Земля. Разработана и реализована методика расчета значения импульса, позволяющего перевести аппарат с низкой околоземной орбиты на орбиту вокруг точки либрации. Исследована взаимосвязь между параметрами перигея отлетной траектории и характеристиками орбиты вокруг точки либрации, на которую он осуществляется. Разработан набор компьютерных инструментов, позволяющих осуществлять поиск времени старта и расчет траекторий перелета на орбиту с заданными характеристиками.



ABSTRACT

The work is devoted research the single-burn transition from low-Earth orbit altitude of 500 km on quasi-periodic orbits around L₂ libration point Sun-Earth system. The aim is to develop and implement a method of calculation the value of maneuver, allowing to transfer the spacecraft from low-Earth orbit to the vicinity of the libration point. The relationship between the parameters of the perigee departure trajectories and the characteristics of the orbit around the libration point, at which the transition occurs, was investigated. A set of computer tools to analyze the start time and calculate the transfer into orbit with given characteristics have been developed.

Supervisor: ______________



Оглавление

Введение

1. Цель работы

2. Математическая модель и инструментарий расчета

2.1 Математическая модель.

2.2 Инструментарий и алгоритмы

2.2.1 Алгоритм подбора начальной скорости КА

2.2.2 Алгоритм подбора величины корректирующего импульса и моделирование отклонений от номинальных значений параметров

3. Расчет и анализ траектории перелета на ограниченную орбиту вокруг точки либрации L2 системы Солнце-Земля

3.1.1 Типы ограниченных орбит вокруг точки L2 системы Солнце-Земля

3.2 Взаимосвязь характеристик отлетного вектора и амплитуд орбиты вокруг точки либрации L2 системы Солнце-Земля

3.3 Взаимосвязь времени старта с возможными характеристиками отлетного вектора на низкой околоземной орбите

3.4 Алгоритм расчета траектории перелета на ограниченную орбиту с заданными характеристиками

3.5 Расчет траектории выхода космического аппарата на гало-орбиту с заданной амплитудой

Заключение

Библиографический список



Введение

В системах двух массивных тел, вращающихся вокруг общего барицентра с постоянной угловой скоростью, существует пять точек, будучи помещенным в которые, тело с бесконечно малой массой будет находиться в состоянии относительного равновесия [1]. Такие точки называются точками либрации или точками Лагранжа и представляют собой частный случай решения ограниченной задачи трех тел [2][3].

Каждая из точек Лагранжа лежит в плоскости орбит массивных тел. Масса одного из крупных тел должна быть больше массы другого. Точки либрации принято обозначать заглавной латинской буквой L с числовым индексом от 1 до 5. Схематичное изображение расположения точек либрации в ограниченной задаче трех тел представлено на Рис.1.

Точки L₁, L₂ и L₃ расположены на оси, соединяющей массивные тела, и называются коллинеарными: L₁ находится между двумя телами системы, ближе к менее массивному телу, L₂ -- снаружи, за менее массивным, и L₃ -- за более массивным. Эти три точки Лагранжа являются неустойчивыми, таким образом, объект, помещенный в одну из них, со временем неизбежно удалится от них. Неустойчивость коллинеарных точек либрации круговой ограниченной задачи трех тел следует из теоремы Ляпунова об устойчивости по первому приближению [3].

Точки L₄ и L₅ называются «троянскими» или же треугольными и являются устойчивыми. Их стабильность обуславливается тем, что расстояния от любой из этих двух точек до массивных тел одинаковы, а значит, и силы притяжения соотносятся в той же пропорции, что и их массы, и, следовательно, результирующая сила направлена к барицентру системы, а также результирующее ускорение связано с расстоянием до центра масс в той же пропорции, что и для двух массивных тел. Поскольку барицентр является одновременно и центром вращения, результирующая сила равна той, что нужна для удержания тела в точке Лагранжа в орбитальном равновесии с остальной системой [4], [5], [6].

Рис.1. Схема пяти точек либрации в системе двух тел (Солнце и Земля). Точки L3, L4, L5 показаны на самой орбите, хотя фактически они будут находиться немного за ней

Впервые о возможности использования точек либрации при разработке космических предположил Артур С. Кларк в 1950 году, в 1966 году Роберт Фаркуар начал проектирование первой космической миссии к одной из таких точек. В 1968 году Фаркуаром и Кэмелом было обнаружено, что увеличение амплитуды осцилляций в орбитальной плоскости может привести к образованию особого класса орбит - гало-орбит [3]. Такие орбиты образуются вокруг коллинеарных точек либрации (L1, L2, L3) при совпадении периодов обращения космического аппарата вокруг точки либрации в плоскости эклиптики и в плоскости, перпендикулярной ей. Периоды орбиты определяются ее амплитудами.

Первое практическое использование точек либрации было осуществлено Робертом Фаркуаром в рамках «Международной программы по исследованию Солнца-Земли» ISEE (InternationalSun/EarthExplorer) [4][7]. Миссия была успешно осуществлена в 1978 году. Космический аппарат ISEE-3 был размещен вблизи точки либрации L₁ системы Солнце-Земля.

Помимо миссии ISEE-3, точки либрации системы Солнце-Земля нашли себе применение в ряде известных космических миссий:

· ISEE-3 (NASA) к точке либрации L₁ в 1978 году для изучения частиц солнечного ветра и космических солнечных лучей;

· WIND (NASA) к точке либрации L₁ в 1994 для изучения частиц солнечного ветра и магнитосферы Земли [6];

· SOHO (ESA-NASA) к точке либрации L₁ в 1996 для изучения солнечной динамики и глубинных слоях Солнца [6];

· ACE (NASA) к точке либрации L₁ в 1997 для изучения частиц солнечного ветра и межзвездной среды;

· MAP (NASA) к точке либрации L₂ в 2001 для изучения реликтового излучения;

· GENESIS (NASA) к точкам либрации L₁, L₂ Солнце-Земля в 2001 для изучения частиц солнечного ветра;

· WSO к точке либрации L₂ Солнце-Земля в 2006 для изучения электромагнитного излучения Вселенной в ультрафиолетовом диапазоне;

· FIRST/HERSCHEL (ESA) к точке либрации L₂ в 2007 для изучения инфракрасного излучения в космосе;

· PLANK (ESA) к точке либрации L₂ в 2007 для изучения реликтового излучения;

· TRIANA (NASA) к точке либрации L₁ в 2008 в качестве космической обсерватории земли;

· GAIA (ESA) к точке либрации L₂ в 2010-2012 для составления карты распределения звёзд нашей Галактики;

· NGST/JWST (NASA) к точке либрации L₂ в 2011 в качестве космического телескопа для изучения космического пространства;

· Constellation X (NASA) к точке либрации L₂ в 2013 для изучения космических объектов по их рентгеновскому излучению;

· DARWIN (ESA) к точке либрации L₂ в 2014 для наблюдения экзопланет и поиска жизни на них;

· TPF (NASA) к точке либрации L₂ в 2015 для исследования окрестности далеких звезд в поисках экзопланет, схожих с планетами земной группы;

· SAFIR (NASA) к точке либрации L₂ в 2015 в качестве инфракрасного телескопа;

Неустойчивость орбит вокруг коллинеарных точек либрации приводит к тому, что малые отклонения от номинального вектора состояния космического аппарата (КА) приводят к сходу аппарата с орбиты. Таким образом, длительное нахождение КА на орбите вокруг точки либрации требует периодических коррекций его движения. Стратегии удержания КА в окрестности точки либрации отличаются частотой, применения корректирующих импульсов и методами расчета их значений [3][4].

Данная работа посвящена исследованию возможностей одноимпульсного перехода с низкой околоземной орбиты высотой 500 км на квазипериодические орбиты вокруг точки либрации L₂ системы Солнце-Земля.

Результаты данного исследования были получены по заказу научно-производственного объединения имени Лавочкина в рамках реализации мероприятий по проектированию миссий «Спектр-РГ», которую предполагается запустить в 2017, и «Спектр-М», запуск которой планируется на 2019 год [8].

Концепция проекта «Спектр-РГ» была сформирована ещё в 1987 году совместно учёными СССР, Финляндии, ГДР, Дании, Италии и Великобритании, однако возможности для осуществления данной миссии появились лишь в 2002 году. Тогда же и возобновилась разработка проекта. С помощью космической обсерватории «Спектр-РГ» планируется производить обзор космического пространства в рентгеновском и гамма-диапазоне электромагнитного спектра, а также совершать поиск скоплений галактик, изучение чёрных дыр, нейтронных звёзд, вспышек сверхновых и галактических ядер. Аппарат с массой 2385 кг (масса топлива 370 кг) будет запущен в точку либрации L₂, где уравновешивается тяготение Солнца и Земли, и станет первым российским аппаратом в окрестности этой точки.

Обсерватория «Спектр-М» ? это следующая российская миссия, в рамках которой планируется разместить КА в окрестности точки либрации Солнце-Земля. Миссия будет направлена на исследования удаленных объектов Вселенной в миллиметровом и инфракрасном диапазонах.

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

1. Цель работы

Целью работы является разработка инструментария, позволяющего рассчитывать траектории перелета от Земли в окрестность точки либрации L₂ системы Солнце-Земля, обеспечивающей переход без приложения корректирующих импульсов на ограниченную орбиту с заданными параметрами вокруг этой точки либрации.

Задачи работы:

· Разработка и реализация алгоритма вычисления значений переходного импульса, позволяющего вывести аппарат с заданным вектором состояния на ограниченную орбиту вокруг точки L₂ системы Солнце-Земля, путем изменения модуля его скорости.

· Построение численной аппроксимации взаимосвязи характеристик отлетного вектора КА с амплитудами орбиты вокруг точки либрации L₂ системы Солнце-Земля.

· Реализация алгоритма расчета времени даты старта, обеспечивающего заданный вектор отлета с низкой околоземной орбиты.

· Исследование взаимосвязи параметров ограниченной орбиты и характеристик отлетного вектора.

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

· Осуществление тестовых расчетов траекторий перелета КА на ограниченные орбиты вокруг точки либрации L₂ системы Солнце-Земля, обладающие заданными характеристиками.

2. Математическая модель и инструментарий расчета

2.1 Математическая модель

В данной работе, для описания движения КА, была использована вращающаяся система координат с фиксированным направлением Солнце-Земля, центр которой находится в точке L₂. Ось X направлена вдоль прямой, соединяющей Солнце и Землю, в направлении от Солнца к Земле, ось Z направлена к северному полюсу эклиптики, ось Y дополняет систему до правой тройки. Схематичное изображение системы координат представлено на Рис. 2.1.

Рис. 2.1. Система координат

В некоторой инерциальной системе координат уравнения движения КА для системы N-массивных тел могут быть представлены в виде:

,(2.12.2@ \* MERGEFORMAT )

где - гравитационная постоянная, - количествo притягивающих центрoв, R - радиус-вектор КА, - мaсса i-гo тела, - радиус-вектор i-го телa.

В случае ограниченной задачи трех тел, где количество массивных тел , при переходе вращающейся системе координат, уравнения (2.1), представленные выше, могут быть приведены к следующему виду [7], [9], [10]:

(2.2)

где -зависящий от масс тел параметр, - возмущающие ускорения, являющиеся функциями координат КА.

В результате линеаризации системa (2.2) может быть приведенa к виду [1], [6], [9], [11]:

(2.3)

Система (2.3) имеет решение в виде:

(2.4)

где , , , ,

, - фаза колебаний в плоскости XY, - фаза колебаний по оси Z. Коэффициенты , , , и фазы , непрерывно зависят от начальных условий.

Рeшeния и составляют линейную комбинацию трех характеристик: возрастающей по модулю (неустойчивой) A₁e^лt, убывающей по модулю к нулю (устойчивой) A₂e^-лt частей и ограниченной. Поиск вектора состояния, соответствующего определенной ограниченной орбите, сводится к поиску начальных условий, где значение коэффициента A₁ будет равным нулю.

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

(2.5)

гдe , , - непрерывные функции своих аргументов, ограниченные по t, , - непрерывные функции своих аргументов, по убывающие по модулю при возрастании t и , - нeпрерывные функции своих аргументов, возрастающие по модулю при возрастании . Причeм , возрастают по первому аргументу и .

Соответствующие нулевому коэффициенту решения образуют устойчивое многообразие, характерное тем, что при они не покидают ограниченной окрестности точки L₂. Аналогично, соответствующие нулевому коэффициенту решения образуют неустойчивое многообразие, характерное тем, что при они не покидают ограниченной окрестности точки L₂ [12][13][14].

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

2.2 Инструментарий и алгоритмы

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

Для моделирования описания движения космического аппарата на орбите и моделирования его траектории использовался программный пакет GMAT (General Mission Analysis Tool). GMAT - это открытый программный комплекс, предназначенный для анализа, оптимизации и моделирования траекторий космических аппаратов, разработанный группой космических корпораций при участии NASA. Программа находится в открытом доступе и с открытыми исходниками, оснащена скриптовым языком, позволяющим производить численное интегрирование уравнения движения космического аппарата в реалистичной модели сил с помощью методов Рунге-Кутта различных порядков и другие. В данной работе был использован метод Рунге-Кутта 8-9 порядка. С помощью языка сценариев в GMAT были описаны алгоритмы, которые производили подбор начальной скорости аппарата, величины маневров, а также интегрирование уравнений движения космического аппарата по орбите. Тексты сценариев представлены в приложении.

2.2.1 Алгоритм подбора начальной скорости КА

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

Алгоритм заключается в подборе таких начальных условий для системы уравнений (2.1), которые обеспечивают равенство нулю коэффициента в ее решении (2.5). Данная задача решается итерационно. На первом шаге необходимо выбрать две плоскости таким образом, чтобы ограниченная орбита, на которой должен находиться аппарат, располагалась между ними. Для заданного вектора состояния численное интегрирование орбиты производится до момента пересечения с одной из плоскостей x_max или x_min. Конечная координата x_c является функцией начального вектора состояния КА, причем x_c = x_max, если и x_c = x_min, если . На множестве векторов состояния, обеспечивающих ограниченную орбиту, функция x_c терпит разрыв. Таким образом, задача поиска начальных условий, приводящих к ограниченной орбите, сводится к задаче поиска точки разрыва функции x_c определяемой численно. Из уравнений (2.4) следует, что если значение коэффициента меньше нуля, то аппарат отклоняется от орбиты в сторону отрицательных значений координаты X; в случае, когда коэффициент больше нуля, аппарат отклоняется в сторону положительных значений X.

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

Рис 2.1. Визуализация работы алгоритма подбора начальной скорости КА.

На Рис. 2.1 приведен пример поиска начальной скорости V_y для КА, находящегося в плоскости XZ и двигающегося ортогонально ей:

X = -277548 км,

Y = 0,

Z = 200000 км,

V_x = V_z = 0.



Рис. 2.2. Блок-схема алгоритма подбора начальной скорости и корректирующего импульса.

Данным методом была рассчитана начальная скорость V_y, равная -0.372794445417389 км/с. Скорость, подобранная с максимальной точностью, позволяет рассчитать в реалистичной системе сил орбиту, остающуюся в пределах заданной окрестности точки либрации L₂ системы Солнце-Земля до 900 суток (порядка 4 полных оборотов). На Рис.2.2 показана блок-схема алгоритма подбора начальной скорости.

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

2.2.2 Алгоритм подбора величины корректирующего импульса и моделирование отклонений от номинальных значений параметров

Блок-схема, представленная на Рис.2.2 в предыдущем разделе, подходит для описания алгоритма подбора величины импульса. Как было замечено ранее, чтобы нивелировать неустойчивую (возрастающую) компоненту движения, необходимо производить коррекции движения КА в окрестности точки либрации. В рамках данной работы был создан алгоритм, который аналогичен алгоритму подбора начальной скорости, описанному выше, но при подборе величины импульса учитывается направление его исполнения, изменение скорости по которому приводит к максимальному изменению неустойчивой компоненты.

Задано начальное значение величины импульса коррекции SK, направляющие и орта вектора коррекции, положенного в плоскость XY (направление импульса DirX и DirY). На текущем шаге происходит исполнение импульса, где его компоненты в разложении по осям рассчитываются как SK*DirX и SK*DirY. Затем производится интегрирование уравнений движения КА, после чего значение SK изменяется.

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

3. Расчет и анализ траектории перелета на ограниченную орбиту вокруг точки либрации L2 системы Солнце-Земля

3.1.1 Типы ограниченных орбит вокруг точки L2 системы Солнце-Земля

Анализируя решение (2.4) линеаризованной системы (2.3), можно заключить, что амплитуды орбиты по осям X и Y зависят друг от друга линейно, а амплитуда по Z является независимой, при этом колебания по X и по Y происходят с одной частотой, а колебания по Z - с другой. Однако, это решение справедливо лишь в малой окрестности точки либрации для ограниченной задачи трех тел. В реальности, амплитуды и частоты колебаний координат КА вокруг точки связаны между собой сложным нелинейным законом [15].

Различные комбинации амплитуд и частот колебаний КА вокруг точки либрации приводят к орбитам различных типов. Поскольку амплитуды и частоты орбиты определяются начальными условиями, необходимо уметь сопоставить исходный вектор состояния КА с типом и характеристиками орбиты, которую он порождает. С целью исследования возможных комбинаций амплитуд ограниченных орбит вокруг точки либрации L₂ и возможных типов таких орбит, в работе осуществлен расчет траекторий, начальный вектор состояния которых обладает следующими характеристиками: , . Произведен анализ зависимости типа орбиты и комбинаций ее амплитуд от начальных координат аппарата и .

Ограниченные орбиты, периоды которых по осям X, Y и Z равны времени одного оборота вокруг точки либрации, называются гало-орбитами. При движении по такой орбите аппарат описывает замкнутую кривую, симметричную относительно плоскости XZ. Амплитуды движения в положительном (северном) и отрицательном (южном) направлении оси Z не равны друг другу [16]. В зависимости от того, в каком направлении отклонение от эклиптики является наибольшим, орбита называется северной или южной. При этом для каждой северной гало-орбиты существует симметричная ей южная. На Рис. 3.1 представлены проекции северной и южной гало-орбит на плоскости XZ, YZ и XY.

Рис. 3.1. Примеры северной (красный цвет) и южной (синий цвет) гало-орбит

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

В зависимости от характера рассогласования колебаний КА вокруг точки либрации по различным осям, будем разделять орбиты на два типа: квазигало-орбиты и орбиты Лиссажу [17][18].

Рис. 3.2. Проекции квазигало-орбиты на плоскости XY, XZ, YZ.

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

· Фигура, заметаемая орбитой несимметрична относительно плоскости эклиптики () [19][20].

· Проекция траектории на плоскость YZ никогда не пересекает некоторой окрестности точки L₂, при этом если то ; если то [21][22].

Значения , и характеризуют разброс проекции орбиты на плоскость XZ вокруг некоторой средней линии. Данные характеристики крайне важны на этапе выбора орбиты вокруг точки L₂ системы Солнце-Земля, поскольку обуславливают возможность попадания КА в полутень Земли. Пример квазигало-орбиты с отмеченными характеристиками , , , , , , и приведен на Рис. 3.2.

Пример орбиты Лиссажу приведен на Рис. 3.3. Из рисунка видно, в отличие от квазигало-орбит, габариты орбит Лиссажу симметричны относительно эклиптики и относительно плоскости XZ. При этом амплитуды по X в положительном и отрицательном направлениях не совпадают. Таким образом, для орбит Лиссажу справедливо следующее: , и .

Начальные условия, тип и характеристики орбиты будут зависеть от исходных координат и Соответствующих для КА, находящемуся в плоскости XZ и движущемуся ортогонально ей, н. Рассчитанные карты характеристик , , и , описывающие их зависимость от начальных координат КА и приведены на Рис. 3.4. На Рис. 3.5 приведены карты характеристик , , и .

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



Рис. 3.3. Проекции орбиты Лиссажу на плоскости XY, XZ, YZ.

На Рис. 3.4. хорошо виден разрыв функции Az+ (, ), который связан с переходом от квазигало-орбит к орбитам Лиссажу. Механизм этого перехода проиллюстрирован на Рис. 3.6. На рисунке показаны орбиты, отвечающие различным начальным положениям аппарата, лежащим на прямой км, при этом начальная координата изменяется от значения, соответствующего гало-орбите, лежащего слева от разрыва Az+ (, ) до значений, лежащих справа от него.

Рис. 3.4. Карты характеристик Ax-, Ax+, Az-, Az+

Можно заметить, что при увеличении координаты увеличиваются разбросы орбиты ?Az+ и ?Az-, при этом орбита остается асимметричной относительно эклиптики. Разрыв функции Az+(, ) отвечает точкам, в которых значения ?Az+ и ?Az- достигают соответствующих значений Az+ и Az-. Комбинации (, ), лежащие справа от разрыва приводят к возникновению орбит Лиссажу, а слева - к возникновению квазигало-орбит.

Рис. 3.5. Карты характеристик Ay, ?Ay, ?Az-, ?Az+

Анализируя Рис. 3.5 можно видеть, что существуют множество комбинаций , соответствующее «оврагам» приведенных на нем функций функции , , и . Это множество отвечает периодическим решениям задачи трех тел: гало-орбитам и вертикальным орбитам Ляпунова.



Рис. 3.6. Механизм перехода от квазигало-орбит к орбитам Лиссажу

Значения начальных координат КА, приводящих к гало-орбитам, приведены в табл. 3.1, где приведены начальные значения (X₀, Z₀) для гало-орбит, соответствующие им значения начальной скорости V_y0, амплитуда по Y - A_y , амплитуды A_z-, A_x+ (для гало-орбит A_z+ = Z₀, A_x- = -X₀). Представленные в таблице значения скоростей даны для оценки порядка величин.

Табл. 3.1 Начальные условия, приводящие к гало-орбитам

Z0 = Az+	X0 = -Ax-	Vy0	Az-	Ay	Ax+
700000	-750584	0.590682	-1161363	1123996	27247.93
650000	-644194	0.531138	-1019016	1055838	69348.51
600000	-566256	0.486786	-902519.1	998711.8	97547.11
550000	-504361	0.450970	-801092.7	948178.6	118026.4
500000	-453098	0.420705	-709166.7	902584.7	133514.3
450000	-410038	0.394814	-624317.6	861478.1	145426.4
400000	-373454	0.372361	-544603.3	824361.1	154621.5
350000	-342355	0.352876	-468881.8	791056.5	161665.1
300000	-316304	0.336297	-396527.0	761937.2	167201.8
250000	-294868	0.322479	-326850.3	737146.2	171903.2
200000	-277549	0.311127	-259115.5	716438.1	174948.9
150000	-264191	0.302247	-192948.4	700020.8	176782.2
100000	-254875	0.296057	-128035.5	688465.3	178864.3
50000	-249171	0.292184	-63817.81	681167.7	179489.0

3.2 Взаимосвязь характеристик отлетного вектора и амплитуд орбиты вокруг точки либрации L2 системы Солнце-Земля

Известно, что перелет на орбиту вокруг точки либрации L₂ системы Солнце-Земля может быть осуществлен совершением одного импульса на низкой околоземной орбите [2], [12], [15], [23]. Фактически, данный перелет осуществляется по орбите, принадлежащей устойчивому многообразию решений задачи трех тел, то есть обладающей нулевой неустойчивой составляющей.

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

В момент отлета с рассматриваемой околоземной орбиты, аппарат находится на сфере, описывающей Землю на высоте 500 километров, а вектор скорости направлен по касательной к этой сфере. При осуществлении импульса перелета в направлении движения КА, существует единственное значение , позволяющее перевести аппарат на орбиту, принадлежащую стабильному многообразию. Таким образом, положение и направление движения КА в момент старта с околоземной орбиты однозначно определяет характеристики ограниченной орбиты вокруг точки L₂, на которую он попадет. Целью настоящей главы является исследование влияния параметров перигея отлета на характеристики орбиты вокруг точки либрации при старте с круговой орбиты высотой 500 км.

Перигей отлета можно определить тремя параметрами, два из которых - это угловые координаты, описывающие положение аппарата (RA и DEC), а третий - угол, определяющий направление его движения (AZI). Иллюстрация этих параметров приведена на Рис. 3.7.

Рис. 3.7. Характеристики, описывающие положение и направление отлетного вектора КА с околоземной орбиты

Во вращающейся системе координат с фиксированной осью Солнце-Земля и центром в Земле (ось X направлена от Земли к точке , Z - на северный полюс эклиптики), угол RA - это угол между проекцией радиус-вектора КА на эклиптику и осью X, DEC - угол между радиус-вектором КА и плоскостью эклиптики, AZI - угол между вектором, направленным по меридиану в сторону Северного Полюса, и вектором скорости КА.

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

Характеристики разброса орбит, для которых были построены карты, наглядно представлены на Рис. 3.2. A_z+ - амплитуда в положительном направлении оси Z, A_z- - в отрицательном направлении оси Z, A_y - амплитуда по оси Y, A_x+ - амплитуда в положительном направлении оси X, A_x- - амплитуда в отрицательном направлении оси X. ?A_z+, ?A_z-, ?A_y ? разбросы орбиты. На Рис. представлены зависимости характеристик и от параметров DEC и AZI при RA, меняющимся от 170° до 200° с шагом в 10 градусов.

Так, например, для при RA = 200° можно увидеть, что наиболее близкие к гало орбите траектории можно получить при параметрах перигея отлета порядка DEC = 20° и AZI = 60°, при этом полученной орбиты составит порядка 350000 км.

На рис. представлены зависимости характеристик Ay+ от параметров DEC и AZI при RA, меняющимся от 170° до 200° с шагом в 10 градусов.



Рис. 3.8. Суперпозиция цветовой карты значений A_z+ и контурной карты разброса ДA_z+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 170°

космический аппарат орбита солнце

Рис. 3.9. Суперпозиция цветовой карты значений A_z+ и контурной карты разброса Д A_z+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 180°



Рис. 3.10. Суперпозиция цветовой карты значений Az+ и контурной карты разброса Д Az+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 190°

Рис. 3.11. Суперпозиция цветовой карты значений A_z+ и контурной карты разброса Д A_z+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 200°



Рис. 3.12. Карта значений Ay+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 170°

Рис. 3.13. Карта значений Ay+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 180°

Рис. 3.14. Карта значений Ay+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 190°

Рис. 3.15. Карта значений Ay+ в зависимости от характеристик отлетного вектора DEC и AZI для RA = 200°

На рис. представлены зависимости значений импульса перехода от параметров DEC и AZI при RA, меняющимся от 170° до 200° с шагом в 10 градусов.

Рис. 3.16. Карта значений импульса перехода в зависимости от характеристик отлетного вектора DEC и AZI для RA = 170°

Рис. 3.17. Карта значений импульса перехода в зависимости от характеристик отлетного вектора DEC и AZI для RA = 180°



Рис. 3.18. Карта значений импульса перехода в зависимости от характеристик отлетного вектора DEC и AZI для RA = 190°

Рис. 3.19. Карта значений импульса перехода в зависимости от характеристик отлетного вектора DEC и AZI для RA = 200°

Приведенные зависимости свидетельствуют о том, что возрастание широты DEC приводит к значительному увеличению амплитуд орбиты вокруг точки либрации и значения импульса, необходимого для перехода на эту орбиту. Влияние направления движения КА (AZI) на характеристики орбиты вокруг точки либрации уменьшается с увеличением DEC, это проявляется в том, что линии уровня на Рис. 3.8 - 3.15 становятся практически вертикальными при больших значениях DEC.

3.3 Взаимосвязь времени старта с возможными характеристиками отлетного вектора на низкой околоземной орбите

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

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

Расчет этого множества выполнен при помощи сценария GMAT, который позволяет вывести в текстовый файл вектора (DEC, AZI), полученные при интегрировании уравнений движения КА до заданного RA, начиная с заданного вектора состояния в гринвичской системе, с заданным шагом по начальному времени.



Рис. 3.20. Множество отлетных векторов на плоскости DEC-AZI в течение одного дня (1 февраля)

Для анализа связи времени старта с возможными характеристиками перигея отлета произведен расчет множества точек RA-DEC-AZI, соответствующих, следующему вектору старта в гринвичской системе координат:

X = -3214 км,

Y = 2880 км,

Z = 5355 км,

V_X = -4.373 км/с ,

V_Y = -5.822 км/с ,

V_Z = 0.506 км/с.

Моделирование осуществлялось для даты старта, которая изменялась в диапазоне одного года с шагом 0.01 суток. Интегрирование велось до значения RA = 160°. Результаты моделирования для имеющегося начального вектора представлены на Рис. 3.20 - 3.23.

Для фиксированного значения RA, точка, соответствующая моменту старта, описывает на плоскости DEC-AZI кривую, проиллюстрированную на Рис. 3.20. При изменении даты старта на один день, данная кривая трансформируется, ее эволюция в течение месяца проиллюстрирована на Рис. 3.21, в течение двух месяцев - на Рис. 3.22.

На Рис. 3.23 представлены кривые, соответствующие множеству отлетных векторов на плоскости DEC-AZI лишь для первого дня января, апреля, июля и октября. В течение года данная кривая заметает множество возможных отлетных векторов во вращающейся системе координат с фиксированной осью Солнце-Земля, показанное серым цветом на Рис. 3.24.

Рис. 3.21. Множество отлетных векторов на плоскости DEC-AZI в течение одного месяца (февраль). Кривая, соответствующая первому дню, отмечена красным.



Рис. 3.22. Множество отлетных векторов на плоскости DEC-AZI с февраля по апрель. Кривая, соответствующая первому дню отмечена красным, первому месяцу - синим.

Рис. 3.23. Множество отлетных векторов на плоскости DEC-AZI в течение одного года. Кривые, соответствующие соседним месяцам, отмечены чередующимися цветами.



Рис. 3.24. Множество отлетных векторов на плоскости DEC-AZI в течение одного года. Отображены только первые дни каждого четвертого месяца.

Рис. 3.25. Множество отлетных векторов на плоскости DEC-AZI в течение одного года. Отображены только первые дни каждого месяца.

3.4 Алгоритм расчета траектории перелета на ограниченную орбиту с заданными характеристиками

При старте с круговой орбиты заданной высоты, вывод на которую осуществляется с определенного космодрома, диапазон характеристик орбиты вокруг точки L₂, на которую возможен одноимпульсный перелет, ограничен [24]. Использование заранее рассчитанных карт соответствия характеристик орбиты вокруг L₂ параметрам перигея отлета (RA, DEC, AZI) позволяет проверить осуществимость перелета на заданную орбиту вокруг точки либрации при старте с заданного космодрома в заданном временном окне.

Выбор времени старта и момента осуществления переходного импульса, переводящего КА с низкой околоземной орбиты на требуемую орбиту вокруг точки либрации, может быть осуществлен методом простого перебора. При наличии предварительно рассчитанной карты, связывающей параметры RA, DEC и AZI с характеристиками орбиты вокруг L₂ данная задача решается на современном ПК за достаточно короткое время.

Было разработано программное обеспечение, позволяющее решать данную задачу. Предполагается, что старт осуществляется с круговой орбиты, высотой 500 км. Пользователь задает вектор выведения на круговую орбиту в Гринвичской системе координат, диапазон, в котором осуществляется поиск времени старта и ограничения на характеристики требуемой орбиты вокруг точки L₂. Программа возвращает набор возможных моментов старта, каждому из которых соответствует момент осуществления переходного импульса, значение и характеристики получаемой орбиты вокруг L₂.

3.5 Расчет траектории выхода космического аппарата на гало-орбиту с заданной амплитудой

В рамках данной работы производился расчет параметров отлетного вектора при заданных ограничениях на геометрию орбиты. С учетом заданных характеристик Az Az+ находится вектор старта по которому определяется отлетный вектор.

В программном комплексе Matlab были заданы заранее рассчитанные соответствия между отлетным вектором (RA, DEC, AZI) во вращающейся системе координат с фиксированной осью Солнце-Земля и параметрами ограниченной орбиты вокруг точки L₂ (, , ), где - высота зоны потенциального затенения от Земли, попадание в которую для КА нежелательно. По заданному отлетному вектору (RA, DEC, AZI) осуществляется интерполяция значений (, , ) между ближайшими значениями.

Разработанный в GMAT сценарий осуществляет подбор отлетных параметров в заданном диапазоне времени и при заданных ограничениях на сверху и снизу и ограничении на снизу. Далее производится корректировка значения импульса dV, полученного ранее, а затем рассчитывается перелет в окрестность точки L₂ и 30 оборотов КА на ограниченной орбите с применением корректирующих импульсов.

Листинги сценариев и функций Matlab приведены в приложениях 1-4.

Алгоритм работы с указанными инструментами для получения орбиты перелета при заданных ограничениях на геометрию орбиты работает по следующей схеме:

1. В первом сценарии (transfer_prediction.script) задаются:

a. Вектор выхода на околоземную орбиту во вращающейся гринвичской системе координат;

b. Временной отрезок выхода на околоземную орбиту [T0, T1] и шаг dt, с которым будет производиться перебор стартового времени внутри отрезка;

c. Ограничения:

i. - ограничение на A_z+ снизу;

ii. - ограничение на A_z+ сверху;

iii. - ограничение на H_z+ снизу.

2. Сценарий (transfer_prediction.script) запускается на выполнение, в результате работы которого получается текстовый файл, каждая строка которого содержит:

a. время выхода на околоземную орбиту t,

b. отлетный вектор (RA, DEC, AZI) во вращающейся системе координат с фиксированной осью Солнце-Земля,

c. параметры ограниченной орбиты A_z+, H_z+,

d. приближенное значение отлетного импульса dV,

e. отлетное время,

f. отлетный вектор состояния в гринвичской системе координат.

3. Полученное время выхода на околоземную орбиту, отлетная координата RA, приблизительное значение отлетного импульса dV вводятся во второй сценарий (TRANSFER.script), который:

a. уточняет значение импульса,

b. применяет его,

c. интегрирует уравнения движения до достижения окрестности точки L2,

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

e. строит проекции полученной орбиты.

Пример 1. (Az+ в диапазоне 350-400 тыс. км.)

Пример траектории номинальной орбиты перелета был рассчитан с помощью созданных инструментариев GMAT в модели сил, с использованием гравитационной модели, учитывающей влияние всех планет Солнечной систем На Рис. 3.26-3.28 приведены проекции орбиты траектории КА на оси XY, XZ и YZ,.

Поиск данной орбиты осуществлялся следующим образом. Был задан вектор выведения КА на круговую околоземную орбиту в Гринвичской системе координат:

X = -3214 км

Y = 2880 км

Z = 5355 км

V_X = -4.373 км/с

V_Y = -5.822 км/с

V_Z = 0.506 км/с).

Также был задан диапазон дат старта, внутри которого осуществлялся поиск (1-31 марта 2015 года), и требования на характеристики и . Осуществлялся поиск орбит с амплитудой и минимально возможным разбросом .

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

· Время вывода на круговую орбиту: 27082.83041 (01 марта 2015 07:55:12),

· Характеристики отлетного вектора: RA=200°, DEC=22.6°, AZI=67°,

· Величина выхода из эклиптики в южном направлении: =384132 км,

· Время старта с круговой орбиты: 27082.88556 (01 марта 2015 09:14:37),

· Первое приближение величины импульса перелета: = 3.12271 км/с.

Полученные данные использовались для моделирования движения КА до момента осуществления перелетного импульса. В перигее отлета значение импульса корректировалось для получения орбиты, принадлежащей устойчивому многообразию орбит вокруг точки L₂. Для получения номинальной траектории, после достижения окрестности точки L₂, один раз в оборот применялись поправки скорости КА, устраняющие неустойчивую компоненту движения. Интегрирование орбиты осуществлялось в среде GMAT с использованием гравитационной модели, учитывающей влияние всех планет Солнечной системы.



Рис. 3.26. Проекция полученной в результате работы сценария орбиты на плоскость XY

Рис. 3.27. Проекция полученной в результате работы сценария орбиты на плоскость XZ



Рис. 3.28. Проекция полученной в результате работы сценария орбиты на плоскость YZ

Пример 2. (Az+ в диапазоне 750-800 тыс. км.)

На рис. 3.29-3.31 приведены проекции орбиты траектории КА, полученные в результате перелета с низкой околоземной орбиты на квазигало-орбиту вокруг точки L₂ системы Солнце-Земля, построенной с использованием гравитационной модели, учитывающей влияние всех планет Солнечной системы.

Вектор выведения КА на круговую околоземную орбиту в Гринвичской системе координат совпадает с заданным вектором, описанным в предыдущем примере (X = -3214 км; Y = 2880 км; Z = 5355 км; V_X = -4.373 км/с; V_Y = -5.822 км/с;

V_Z = 0.506 км/с), взят тот же диапазон дат старты, но требования на характеристики и изменились: .

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

· Время вывода на круговую орбиту: 27083.35041 (01 марта 2015 20:24:00),

· Характеристики отлетного вектора: RA=170°, DEC=40°, AZI=159°,

· Величина выхода из эклиптики в южном направлении: = 750201 км,

· Время старта с круговой орбиты: 27083.36153 (01 марта 2015 20:40:01),

· Первое приближение величины импульса перелета: = 3.12708 км/с.

Рис. 3.29. Проекция полученной в результате работы сценария орбиты на плоскость XY



Рис. 3.30. Проекция полученной в результате работы сценария орбиты на плоскость XZ

Рис. 3.31. Проекция полученной в результате работы сценария орбиты на плоскость YZ



Заключение

В работе реализована методика расчета ограниченных орбит вокруг точки либрации L₂ системы Солнце-Земля. Движение КА в окрестности точки либрации представляется как суперпозиция трех составляющих: устойчивой (стремящейся к нолю), неустойчивой (стремящейся к бесконечности), и ограниченной. Были разработаны методики, позволяющие осуществить расчеты, помогающие исследовать особенности орбитального движения в окрестности точки либрации L₂ системы Солнце-Земля.

С помощью осуществленных расчетов была построена и проанализирована зависимость типа и характеристик ограниченных орбит вокруг точки L₂ системы Солнце-Земля от параметров начального вектора состояния КА. Использование результатов этих расчетов позволяет подбирать начальные условия, приводящие к периодическим орбитам, квазигало-орбитам и орбитам Лиссажу, обладающими заданными характеристиками.

Разработан набор компьютерных инструментов, позволяющих осуществлять расчеты траекторий, принадлежащих устойчивому многообразию орбит вокруг точки L₂ системы Солнце-Земля. Данные инструменты применены для анализа возможных траекторий одноимпульсного перелета с низкой околоземной орбиты на ограниченную орбиту вокруг точки либрации. Возможности перелета на такие орбиты ограничены параметрами круговой орбиты, с которой осуществляется перелет. Разработанные инструменты позволяют для заданного в Гринвичской системе координат вектора выведения найти дату старта, обеспечивающую перелет КА на ограниченную орбиту вокруг точки L₂, обладающую заданными характеристиками.



Библиографический список

[1] Farquhar, R.W., The Control and Use of Libration-Point Satellites // NASA TR R-346, 1970.

[2] Dunham, D.W., Roberts, C.E., Stationkeeping Techniques for Libration-Point Satellites // The Journal of the Astronautical Sciences, Vol.49, No. 1, January-March 2001, pp.127-144

[3] Gomez, G., Howell, K., Masdemont, J., Simo, C., Station-keeping Strategies for Translunar Libration Points Orbits // AAS 98-168

[4] Simo, C., Gomez, G., Llibre, J., Martinez, R., Rodriguez, J., On the Optimal Station Keeping Conrol of Halo Orbits // Acta Astrinautica, Vol. 15, No. 6/7, pp. 391-397, 1987

[5] H. Hechler, J. Cobos, Herschel, Planck and Gaia Orbit Design // Proceedings of the Conference “Libration Point Orbits and Applications”, Aiguablava, Spain, 10-14 June 2002.

[6] C.E. Roberts, Long term missions at the Sun-Earth libration point L1: ACE, SOHO, and WIND // AAS/AIAA Astrodynamics Specialist Conference; 31 Jul. - 4 Aug. 2011; Girdwood, AK; United States.

[7] R. W. Farquhar and A. A. Kamel. Quasi-Periodic Orbits About the Translunar Libration Point. Celestial Mechanics, 7(4):458-473, 1973.

[8] Ильин И.С., Сазонов В.В., Тучин А.Г., Построение ограниченных орбит в окрестности точки либрации L2 системы Солнце-Земля // Препринты ИПМ им. М.В.Келдыша, 2012, №65. 28 с. URL: http://library.keldysh.ru/preprint.asp?id=2012-65

[9] C. Ocampo, An Architecture for a Generalized Spacecraft Trajectory Design and Optimization System // Proceedings of the Conference “Libration Point Orbits and Applications”, Aiguablava, Spain, 10-14 June 2002.

[10] J.A. Kechichian, E.T. Campbell, M.F. Werner, E.Y. Robinson, Solar Surveilance Zone Population Strategies with Picosatelliter Using Halo and Distant Retrograde Orbits // Proceedings of the Conference “Libration Point Orbits and Applications”, Aiguablava, Spain, 10-14 June 2002.

[11] G. Gomez, J.J. Masdemont, J.M. Mondelo, Libration Point Orbits: a Survey from the Dynamical Point of View // Proceedings of the Conference “Libration Point Orbits and Applications”, Aiguablava, Spain, 10-14 June 2002.

[12] G. Gomez, A. Jorba, J. Masdemont, C. Simo, Dynamics and Mission Design near Libration Points, Vol. III Advanced Methods for Collinear Points // World Scientific Publishing Co. Ptc. Ltd., 2001

[13] E. Perozzi, S. Ferraz-Mello, Space Manifold Dynamics: Novel Spaceways for Science and Exploration // Springer, 2010

[14] J. Cobos, J. Masdemont, Astrodynamical Applications of Invariant Manifolds Associated with Collinear Lissajous Libration Orbits // Proceedings of the Conference “Libration Point Orbits and Applications”, Aiguablava, Spain, 10-14 June 2002.

[15] D.C.Folta et al., Earth-Moon libration point orbit station-keeping: Theory, Model and operations // Acta Astronautica, 2013, http://dx.doi.org/10.1016/j.actaastro.2013.01.022

[16] C.E. Roberts, The SOHO Mission L1 Halo Orbit Recovery from the Attitude Control Anomalies of 1998 // Proceedings of the Conference “Libration Point Orbits and Applications”, Aiguablava, Spain, 10-14 June 2002.