Исследование зависимости формы ограниченной орбиты КА от начального вектора состояния в окрестности точки либрации Ь2 системы Солнце-Земля_
Исследование зависимости формы ограниченной орбиты КА от начального вектора состояния в окрестности точки либрации Ьд системы Солнце-Земля
Бобер С.А., Аксенов С.А., Николаева Ю.А.
Московский институт электроники и математики НИУВШЭ $1а$. ЬоЬег@^аИ. сот, ак$епо\.8. а@£тай. сот, nikolaeva.juli@gm.ail. сот
Аннотация. В работе описана математическая модель движения КА в окрестности точки либрации Ьг системы Солнце-Земля, методика подбора начальных условий и величин корректирующих импульсов для минимизации неустойчивой компоненты движения. Рассчитаны карты амплитуд и разбросов ограниченных орбит, дана их классификация.
Ключевые слова: точки либрации, орбита, космический аппарат, гало, квазигало, орбиты Лиссажу, орбиты Ляпунова, резонанс, Рунге-Кугта, точка \л, ОМАТ
1 Введение
Точки либрации, или точки Лагранжа, - это частные решения ограниченной задачи трех тел, в которых третье тело (например, космический аппарат1) с пренебрежимо малой массой будет находиться в состоянии равновесия относительно двух массивных тел. Всего таких точек пять: 1л,...,Ь5, три из которых 1л-Ьз являются коллинеарными, т.е. лежат на прямой, соединяющей два массивных тела (Рис. 1).
Движение по орбите вокруг коллинеарных точек либрации является
неустойчивым, т.е. малые отклонения от номинальных значении параметров приводят к выходу аппарата из окрестности точки либрации.
1 Далее термин «космический аппарат» будет сокращенно обозначаться «КА».
Одной из важнейших задач является разработка стратегии удержания аппарата на орбите. Под стратегией понимается периодичность исполнения корректирующих маневров, направление их исполнения и т.д., а также непосредственно методика расчета их величины.
Точки либрации Li и L2 являются самыми востребованными с научной точки зрения, т.к. из их окрестностей возможны и удобны как наблюдения за Землей и Солнцем, так и изучение глубокого космоса. С 1978 по 2015 запущено и планируется запустить на орбиты в окрестностях этих точек [Аксенов С.А. и др., 2013], [Николаева Ю.А. и др., 2013], [Федоренко Ю.В. и др., 2013] 16 космических аппаратов (NASA, ESA), среди которых:
■ ISEE-3 (1978, LI) - первый КА, направленный к точке либрации,
■ WIND (1994, LI) - изучение солнечного ветра и магнитосферы Земли,
■ HERSCHEL (2007, L2) - изучение звезд в инфракрасном диапазоне,
■ PLANK (2007, L2) - изучение реликтового излучения.
НПО им. Лавочкина также разрабатывает два КА предназначенных для запуска в окрестность L2 Солнце-Земля: Спектр-М (в соответствии с федеральной космической программой России на 2006-2015 годы) и Спектр-РГ - совместный с Германией проект [Шматов и др., 2013]. Спектр-М предназначен для астрофизических исследований в миллиметровом, субмиллиметровом и инфракрасном диапазонах, Спектр-РГ - изучения Вселенной в гамма- и рентгеновском жестком диапазоне энергий.
2 Математическая модель
Солнце
с
х
Рис. 2. Система координат
Исследование зависимости формы ограниченной орбиты КА от начального вектора состояния в окрестности точки либрации L2 системы Солнце-Земля_
Движение КА рассматривается во вращающейся системе координат с фиксированным направлением Солнце-Земля2. Центр системы координат находится в точке L2. Ось X направлена от Солнца к Земле, вдоль соединяющей их прямой. Ось Z направлена к северному полюсу эклиптики, ось Y дополняет систему до правой тройки. Указанная система координат представлена на Рис. 2.
Уравнения движения КА в системе N массивных тел могут быть в некоторой инерциальной системе координат представлены в виде:
(D
где п - количество притягивающих центров, G - гравитационная постоянная, R - радиус-вектор КА, m; - масса i-ro тела, Ri - радиус-вектор i-ro тела.
В ограниченной задаче трех тел (N = 2) при переходе к вращающейся системе координат, описанной выше, уравнения (1) могут быть преобразованы к следующему виду [Эльясберг и др., 1985], [Gomez et al., 2002]:
'х - 2у - (2,с2 + 1)х = а1( у - 2х - (1 - с2)у = ау, (2)
Z + C-rZ = Щ,
где с2 - параметр, зависящий от масс тел, ах,ау, ая - возмущающие ускорения, которые являются функциями координат КА и эксцентриситета орбиты.
После линеаризации в окрестности точки L2 система (2) может быть приведена к виду [Farquhar, 1970], [Gomez et al., 2002]:
fx - 2y - (1 + 2c2)x = 0. У — 2x — (1 — c2 )y = 0, (3)
z + c2z = 0
Данная система имеет следующее решение:
x(t) = A_e?t + Aze_At + Asy cosicjt + <pxy) y ft) = k.A.eÀ: - k2A2e~At - k2Asy sinlojt + (4)
z(t) =
где A = J-3L-to= j-3L-, y = ^2, = —, k2 -
2 ^ 2 ' 1 23L ' г 2ffl '
4>щ ~ фаза, колебаний в плоскости XY, фЕ - фаза колебаний по оси Z. Коэффициенты А±, А2, Аху, As и фазы ф.1у, ф2 зависят от начальных условий.
Решения x(t) и y(t) являются линейной комбинацией трех компонент -ограниченной (sm и cos), экспоненциально возрастающей и
2 Далее, если не сказано иное, координаты, скорости, уравнения и др. будут рассматриваться в этой системе координат
убывающей к нулю частей. Далее будем предполагать, что решение
системы (1), записанное во вращающейся системе координат с центром в точке L2, имеет схожую структуру, т.е. состоит из возрастающей, убывающей и ограниченной компонент.
Для длительного удержания КА в окрестности точки либрации необходимо, чтобы возрастающие компоненты системы (5) были минимальны по модулю, т.е. близки к нулю.
В этой работе анализ движения КА строится на численном интегрировании уравнений движения. Подбор начальных условий, равно как и величин корректирующих маневров, обеспечивающих минимизацию возрастающих компонент, происходит алгоритмически. Неточность знания вектора состояния КА и исполнения корректирующих маневров в данной работе во внимание не принимаются.
3 Инструментарий
Для моделирования движения КА в окрестности точки Lг было использовано программное обеспечение GMAT (General Mission Analysis Tool), разрабатываемое NASA и находящееся в открытом доступе вместе с исходными кодами. В данном исследовании интегрирование производилось методом Рунге-Кутта 8-9 порядка в реалистичной модели сил.
В GMAT были созданы сценарии для подбора начальной скорости аппарата, величины маневров, интегрирования уравнений движения КА для вычисления характеристик орбит.
В связи с большим количеством расчетов для составления карт амплитуд (10000 точек) и с тем, что GMAT не позволяет распараллеливать расчеты, была написана программа, позволяющая разделять один общий сценарий на несколько частей. Это позволило задействовать 40 процессорных ядер и сократить время расчета карт с 200 суток до недели.
4 Алгоритм подбора начальной скорости КА
Начальное положение КА ограничено плоскостью XZ, т.е. (Х0, 0, Z0). От конкретного положения в этой плоскости зависит тип орбиты и ее характеристики. Вектор начальной скорости КА выбирается ортогональным плоскости XZ, т.е. (0, Vy, 0).
Таким образом, при заданных координатах КА (Х0, Z0) требуется подобрать компоненту скорости Vy, минимизирующую возрастающие компоненты движения с вычислительной точностью GMAT.
Алгоритм заключается в подборе таких начальных условий для системы уравнений (2), чтобы в ее решении (4) коэффициент возрастающей компоненты равнялся нулю. Данная задача решается итерационно. Из уравнений (4, 5) следует, что если значение коэффициента Аг меньше нуля, то координата х КА экспоненциально
Исследование зависимости формы ограниченной орбиты КА от начального вектора
состояния в окрестности точки либрации Ь2 системы Солнце-Земля_
убывает в сторону отрицательных значений; в противном случае, координата х экспоненциально возрастает в сторону положительных значений. Такое поведение КА в окрестности точки либрации было использовано для подбора компоненты скорости Уу следующим образом. Были выбраны две плоскости X = Хшах и X = Хтшп таким образом, чтобы ограниченная орбита находилась между этими плоскостями. С заданной начальной скоростью производилось численное интегрирование уравнений движения КА в ОМАТ до достижения космическим аппаратом любой из указанных плоскостей. Т. о. была определена функция Б(Уу) - конечная координата X КА, начинающего движение со скоростью Уу, при достижении им одной из плоскостей (Рис. 3).
Очевидно, что данная функция терпит разрыв при Уу = УуО, соответствующей нулевому значению возрастающей компоненты движения. Задача вычисления значения УуО решается методом деления отрезка пополам, который позволяет достигнуть любой наперед заданной точности вплоть до машинного ограничения 10"16.
Рис. 3. Функция Б(Уу) терпит разрыв в точке УуО
Скорость, подобранная с максимальной точностью, позволяет КА находиться на ограниченной орбите на протяжении 4 оборотов (700 суток). Для численного расчета ограниченной орбиты на более долгий период времени необходимы периодические математические коррекции скорости аппарата, компенсирующие возрастающую компоненту движения.
Алгоритм подбора величины вектора корректирующего импульса был построен исходя из аналогичных соображений с единственной разницей -вектор импульса коррекции располагается в плоскости ХУ под определенным заранее углом. Оптимальный выбор данного угла - тема отдельной работы и здесь она не рассматривается.
§ Карты амплитуд ограниченных орбит
На Рис. 4 изображены основные параметры ограниченной орбиты: амплитуды в положительном и отрицательном направлении соответствующих осей: Ах+, Ах-, Аг-, Ау = Ау+ = Ау-; а также
разбросы ААг+, ААг-, ААу.
х 10
Рис. 4. Амплитуды ограниченной орбиты (на примере квазигало)
В данной работе были проведены расчеты карт амплитуд и разбросов орбит для следующих наборов начальных положений КА:
Х0: от -1000000 км до 0 км с шагом в 10000 км, от 0 км до 1000000 км с шагом в 10000 км.
Таким образом, были рассчитаны карты с разрешением 100x100 точек, в каждой из которых было рассчитано 30 оборотов орбиты (15 лет) с подбором начальной скорости, расчётом и применением корректирующих импульсов каждые 2 оборота под углом 45° в плоскости ХУ для удержания КА в окрестности точки либрации.
состояния в окрестности точки либрации Ь2 системы Солнце-Земля
Рис. 5. Карты амплитуд ограниченных орбит Аг+, Ах+, Ау
На Рис. 5 изображены карты, на которых каждой начальной точке траектории КА на плоскости ХЯ ставится в соответствие одна из указанных амплитуд. На карте можно наблюдать разрыв,
происхождение которого объяснено далее.
6 Классификация ограниченных орбит
Анализ амплитудных карт и форм орбит выявил следующие классы.
Гало-орбитами называется класс периодических орбит вокруг коллинеарных точек либрации (Ьх-Ьз). Они образуются при совпадении частот обращения КА в плоскости эклиптики и в плоскости, перпендикулярной ей. В реальной системе сил не существует «чистых» гало-орбит, т.к. из-за возмущающих воздействий от массивных тел Солнечной системы они вырождаются в квазипериодические орбиты, обладающие разбросами ААг+, ААг-, ААу. В этом случае гало-орбигами, среди всех квазигало-орбит, будем называть таковые с наименьшими разбросами ДАг-, ААу.
Орбиты Лиссажу - непериодические орбиты, частоты обращений которых в указанных плоскостях не кратны.
Вертикальные орбиты Ляпунова - класс периодических орбит, проходящих через точку либрации. Кратность частот - 2.
На Рис. 6, 8 наглядно изображено распределение классов орбит по начальным точкам орбит в плоскости XZ.
-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 О
Х,кт х105
Рис. 6. Классы орбит на карте амплитуд Аг+
Разрыв функции Аг+ (ХО, Z0) связан с переходом от квазигало-орбит к орбитам Лиссажу, механизм этого перехода проиллюстрирован на Рис.7. На рисунке показаны орбиты, отвечающие различным начальным положениям аппарата, лежащим на прямой /0=400000 км, при этом начальная координата Х0 изменяется от значения, соответствующего гало-орбите, лежащего слева от разрыва Аг+ (Х0, до значений, лежащих справа от него. Можно видеть, что при увеличении координаты Х0 увеличивается разбросы орбиты ДАг+ и ААг-, при этом орбита остается асимметричной относительно эклиптики. Разрыв функции Аг+ (Х0, Х0>) отвечает точке, в которой значения ААг+ и ДАг- достигают соответствующих значений Аг+ и Аг-, при этом появляются точки, в которых стоимость перехода от северной квазигало-орбиты к южной становится нулевой и орбита начинает "отражаться" от эклиптики, это приводит к тому, что А1+ и Аг- становятся равны между собой, и Аг+ терпит разрыв. Комбинации (Х0, ХО), лежащие справа от разрыва приводят к возникновению орбит Лиссажу.
Исследование зависимости формы ограниченной орбиты КА от начального вектора состояния в окрестности точки либрации Ь2 системы Солнце-Земля_
О.ОЕ+ОО
-1.5Е+06 -1.3Е+06 -1.1Е+06 -9.0Е+05 -7.0Е+05 -5.0Е+05 -З.ОЕ+ОБ -1.0Е+05 1.0Е+05 З.ОЕ+ОБ
X, КМ
Рис. 8. Резонансные орбиты
На Рис. 8 приведена карта начальных условий для получения резонансных орбит. Гало орбиты отвечают случаю сау/щ = 1, начальные
-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0
X, кгп х 105
Рис. 7. Механизм перехода от квазигало-орбит к орбитам Лиссажу
3.0Е+О6
2.5Е+06
2.0Е+06
* 1.5Е+06
N
1.0Е+06
5.0Е+05
условия, приводящие к гало-орбитам показаны красными линиями. Орбиты Ляпунова отвечают ситуации, когда соу/щ = 2. Траектория в этом случае формирует восьмерку в проекции на плоскость YZ, начальные условия, приводящие к таким орбитам, показаны синей линией. Определенные сочетания начальных условий (зеленые линии) могут приводить к сложным резонансным квазигармоническим движениям, образующим многолепестковые структуры, примеры которых приведены на Рис. 8.
Список литературы
[Аксенов С.А. и др., 2013] Аксенов С.А., Ефремова Е.В., Данхэм Д.У., Компьютерное моделирование миссии к точке либрации L2 системы Земля-Луна // Инновационные информационные технологии, 2013, Т. 2, № 2, С. 545-549.
[Николаева Ю.А. и др., 2013] Николаева Ю.А., Аксенов С.А., Данхэм Д.У., Расчет окон запуска космического аппарата для траектории земля - точка L2 системы Земля-Луна // Инновационные информационные технологии, 2013, Т. 2, № 2, С. 567-573.
[Федоренко Ю.В. и др., 2013] Федоренко Ю.В., Аксёнов С.А., Данхэм Д., Оценка времени видимости космического аппарата при движении вокруг точки либрации L2 системы Земля-Луна // Инновационные информационные технологии, 2013, Т. 2, № 2, С. 573-577.
[Шматов и др., 2013] Шматов С.И., Мордвинкин А.С., Возмущающее воздействие солнечного излучения на космический аппарат «Спектр-РГ» на рабочей орбите // Вестник НПО им. Лавочкина, 2013, №5, с. 21.
[Эльясберг и др., 1985] Исследование движения космического аппарата в окрестности коллинеарного центра либрации / Эльясберг П.Е., Синицын В.М., Тихомирова Т.А. // Институт космических исследований (ИКИ) АН СССР, 1985.
[Farquhar, 1970] Farquhar, R.W., The Control and Use of Libration-Point Satellites // NASA TR R-346,1970.
[Gomez et al., 2002] Gomez, G., Marcóte, M., Masdemont, J.J., Trajectory Correction Maneuvers to Libration Point Orbits // 7th International Conference on Libration Point Orbits and Applications, Parador d'Aiguablava, Girona, Spain, 2002.