Моделирование движения космического аппарата по гало-орбите вокруг точки L2 системы Солнце-Земля
Николаева Ю.А., Аксенов С.А.
Московский институт электроники и математики НИУВШЭ Nikolaeva.juli(cp,gmail.сот, aksenov.s.a@gmail.com
Аннотация. Данная работа посвящена разработке стратегии удержания космического аппарата на гало-орбите в окрестности точки L2 системы Солнце-Земля, ее программной реализации и ее анализу. Разработанная стратегия применима не только для гало-орбит, но и для любых ограниченных орбит в окрестности коллинеарных точек либрации.
Ключевые слова: гало-орбита, точка либрации, система Солнце-Земля.
1 Введение
Точками либрации в ограниченной задаче трех тел называют точки равновесия тела с меньшей массой относительно массивных тел. В исследуемой задаче телом с меньшей массой является космический аппарат (КА), телами с большей массой - Солнце и Земля. Существуют 5 точек либрации. Точки либрации Ll, L2 и L3 лежат на прямой, соединяющей массивные тела, L4, L5 - на орбите одного из массивных тел. В данной задаче точка L2 расположена за Землей.
Гало-орбитами называется класс периодических орбит вокруг коллинеарных точек либрации (точек, лежащих на прямой). Они образуются при совпадении периодов обращения КА в плоскости эклиптики и плоскости, перпендикулярной ей.
Впервые гало орбиты были предложены Робертом Фаркуаром в его диссертации [Farquhar, 1970]. В рамках программы Apollo он предложил использовать дополнительный аппарат, помещенный на гало-орбиту вокруг точки L2 системы Земля-Луна, для связи с луноходом, расположенным на обратной стороне Луны. Несмотря на то что данная идея не была реализована, впоследствии гало-орбиты широко использовались в космических миссиях. Наиболее известными из них являются:
■ ISEE-3, АСЕ - помещены на орбиту вокруг точки L1 системы Солнце-Земля;
■ Herschel, Plank - помещены на орбиту вокруг точки L2 системы Солнце-Земля.
Движение аппарата по гало-орбите является неустойчивым, т.е. малое отклонение от начальных параметров приводит к отклонению КА от
орбиты и в итоге к выходу из окрестности точки либрации. Для поддержания движения КА на орбите используются т.н. импульсы коррекции (корректирующие импульсы), призванные нивелировать неустойчивую компоненту движения. Комплекс мер, принимаемых для удержания КА на орбите, называется стратегией удержания КА на гало-орбите. Под стратегией обычно понимаются периодичность, место совершения, величина и направление импульсов коррекции, а также непосредственно методика их расчета. Существуют несколько концепций стратегий удержания КА на орбите, но они, как правило, требуют существенных вычислений или сложны в разработке [Simo et al., 1987, Dunham et al., 2001, Gomez et al.]
В реальности существуют технические ограничения точности определения положения скорости аппарата, а также ограничения на минимальную величину импульса и точность его исполнения [Николаева и др., 2014]. В данном исследовании указанные ограничения не учитываются, т.е. предполагается, что параметры аппарата известны точно и нет ограничений на величину корректирующего импульса.
Задачей данного исследования является разработка и программная реализация стратегии удержания КА на гало-орбите и ее анализ. Полученная в результате исследования стратегия и ее программная реализация применимы для любых ограниченных орбит в окрестности точки L2 системы Солнце-Земля.
2 Математическая модель
Для описания движения КА по гало-орбите введем систему координат, связанную с точкой L2. Центр системы координат расположен в точке L2, ось X направлена вдоль прямой, соединяющей Солнце и Землю, от Солнца к Земле, ось Z направлена в северный полюс эклиптики, ось Y дополняет систему до правой тройки. В данной системе координат гало-орбитой является орбита, периоды обращения КА по которой по осям Z и Y совпадают.
В некоторой инерциальной системе координат уравнения движения КА могут быть представлены в виде:
R=GEf=imij^, (1)
где п - количество притягивающих центров, G - гравитационная постоянная, пц. - масса i-ro притягивающего центра, R - радиус вектор КА, Rj - радиус-вектора i-ro притягивающего центра.
В описанной выше системе координат при переходе к ограниченной задаче трех тел (п = 2) уравнения (1) преобразуются к виду [Gomez et al., 2002; Cobos et al., 2002; Ильин и др., 2012]:
iji - 2y - (2c2 + l)x =
У + 2k- (1 - c2)y = ay, (2)
2 + €2Z =
где с2 - параметр, зависящии от масс тел, ах, ау, а^ - возмущающие ускорения, которые являются функциями координат КА и эксцентриситета орбиты.
Уравнения (2) можно линеаризовать в окрестности точки Ь2. Тогда они принимают следующий вид [НесЫег & а1., 2002, Ильин и др., 2012]:
(х 2у (1 2с2)х = 0, у + 2х - (1 - с2)у = 0, (3)
2 + С22 = О
Система уравнений (3) имеет следующее решение:
х(г) = А1еЛ1: + А2е-М + Ахусоз(а»1 + ф^у)
¡уф ='к1А1е?Л- к2 А2 е~м - к2Аху нт(аЛ + «р^) (4)
=
ГДе А = 4-2-' « = 4-2-' У = к1 = <Р*У
- фаза колебаний в плоскости ХУ, фЕ - фаза колебаний по оси г. Коэффициенты Аь А2, Аху, Ах, и фазы ф~ зависят от начальных условий.
Решения х(1) и у(1) являются линейной комбинацией трех компонент: ограниченной, экспоненциально возрастающей и экспоненциально убывающей. Далее будем предполагать, что решение системы (2) имеет схожую структуру, т.е. состоит из ограниченной, возрастающей и убывающей компонент:
Аху^Ит
О-Фху)
[уСО = к1А1||лпс (0 + к1А2фйес (г) + к2Ахудрцт (г, фху) (5)
г(1) = А^цщй.Фг)
где А^Сй, А^-ДЛ) - возрастающие компоненты, А2^1асГг), А2фй&с(0 -убывающие компоненты, Аху\|1Н
ПИ.1' Ч'ХУ..' 5 М1Г.1
(ЬФг) -
ограниченные компоненты.
Коэффициенты А1? А2 могут быть как положительными, так и отрицательными, поэтому далее по тексту под «возрастающей» компонентой будет подразумеваться возрастающая по модулю компонента, аналогично для «убывающей».
3 Алгоритм подбора начальной скорости и величины корректирующего импульса
Для длительного удержания КА на гало-орбите необходимо, чтобы коэффициент перед возрастающей компонентой равнялся нулю. В данном исследовании производилось численное интегрирование уравнений движения, поэтому подбор начальных условий, при которых коэффициент при возрастающей компоненте обращается в ноль, производился алгоритмически.
Кроме того, численное моделирование не может обеспечить бесконечной точности, поэтому для расчета номинальной гало-орбиты
необходимо совершать корректирующие импульсы. Подбор величин корректирующих импульсов также производился алгоритмически, а в основе методики подбора величины лежит та же идея.
Для моделирования движения КА по гало-орбите использовался пакет GMAT (General Mission Analysis Tool), позволяющий производить интегрирование уравнений движения КА в реальной модели сил с помощью различных численных методов. В GMAT был создан сценарий, в котором описываются начальные параметры КА, подбирается его начальная скорость, величины маневров, а также производится численное интегрирование уравнений движения. В данном исследовании использовался метод Рунге-Кутта 8-9 порядка.
3.1 Алгоритм подбора начальной скорости
Предположим, что КА в начальный момент времени находится в плоскости XZ (Y=0) и двигается ортогонально ей (Vx=Vz=0). Координаты X и Z обеспечивают гало-орбиту. Производится подбор скорости Vyo, обеспечивающей максимальное время нахождения КА на орбите. Подбор скорости осуществляется итерационно и останавливается, когда достигается наибольшая точность.
Из решения x(t) (4) видно, что если значение коэффициента возрастающей компоненты отрицательно, то аппарат отклоняется в сторону отрицательного направления оси, и наоборот. При подборе скорости Vyo данный факт был использован следующим образом. Введем в пространстве виртуальные плоскости Х=Хтах и Х=Хтш таким образом, чтобы гало-орбита заведомо лежала между ними. Затем проинтегрируем уравнения движения до пересечения с этими плоскостями. Т.о. получаем функцию конечной координаты X от начальной скорости Vy. Очевидно, данная функция терпит разрыв при Vy= Vyo. Поиск данного значения производится методом деления отрезка пополам. Иллюстрация метода представлена на рис. 1.
х = х,
X
X = X,
Моделирование движения космического аппарата
по гало-орбите вокруг точки L2 системы Солнце-Земля_
Рис. 1. Подбор начальной скорости КА
Используемый метод позволяет достичь любой точности, вплоть до машинной, которая составляет 10~16. Для гало-орбиты с начальными параметрами X - -277548 км, Y = 0 км, Z - 200000 км, Vx = Vz = 0 км/с найденная скорость Vyo составляет 0.372794445417389 км/с и обеспечивает нахождение аппарата на орбите на протяжении 700 суток (4 оборота).
3.2 Алгоритм подбора величины корректирующего импульса
Алгоритм подбора величины импульса коррекции в целом совпадает с алгоритмом подбора начальной скорости. Основным его отличием является то, что при подборе величины импульса учитывается его направление.
В начальный момент времени задано начальное значение импульса SK, а также направление маневра DirX и DirY (направляющие sin и cos орта вектора коррекции, положенного в плоскость XY). На текущем шаге происходит исполнение импульса, компоненты которого представлены как SK*DirX и SK*DirY. После этого производится интегрирование уравнений движения КА и изменение значения SK.
4 Зависимость затрат характеристического импульса от места и направления совершения маневра
Эффективность импульса коррекции зависит от места и направления его совершения. Введем следующие параметры (см. рис. 2):
■ Угол а - определяет место исполнения импульса;
■ Угол Р - определяет направления исполнения импульса.
Рис. 2. Углы аир
Углу а = 01° соответствует наиболее удаленная от Земли точка орбиты, углам а = ±180® соответствует ближайшая к Земле точка орбиты. Увеличение угла происходит против часовой стрелки.
Важно отметить, что угол Р, по сути, задает ось, вдоль которой может совершаться импульс (как в положительном, так и в отрицательном направлении).
Импульс совершается один раз в оборот в точке, определенной углом
4.1 Нахождение наиболее эффективного места совершения маневра
Для нахождения наиболее эффективного места совершения маневра были исследованы 3 гало-орбиты со следующими начальными координатами:
■ X = -277548 км, У = 0 км, Ъ = 200000 км;
■ X = -373454 км, У = 0 км, Ъ = 400000 км;
■ X = -566256 км, У = 0 км, Ъ = 600000 км.
Проекции движения КА по этим орбитам представлены на рис.3-5.
Рис. 3. Проекция движения КА по гало-орбитам с различными амплитудами на
плоскость ХУ
хЮ9
Рис. 4. Проекция движения КА по гало-орбитам с различными амплитудами на
плоскость ХЪ
т = 2DOODO км г = ДОСООО км г ~ БООООО км
х 10
Рис. 5. Проекция движения КА по гало-орбитам с различными амплитудами на
плоскость XZ
1 ООЕ 09
"i -i
? 1 00Е-10
| 1 00Е-11
: I.OOE-I:
h
2
ШЕ-Ü
[ Ж Л Л у ¿.""Л \/1 Л. \
- 2 = 200000 —г ■ 400000 ■.....г = 600000
i Vi \\ \\ \ & \
{ ^T-'í >1 * * У < Г V* / ; ^-'■•rL-ST /
-:оо
150
:оо
-150 -100 -50 0 50 100
Угол, градусы
Рис. 6. Зависимость математического ожидания значения импульса коррекции от места исполнения импульса для различных амплитуд
На рис. 6 представлена зависимость математического ожидания значения импульса от места исполнения маневра (угол а) при направлении исполнения маневра Р=0°. Для каждого угла а было смоделировано движение КА по гало-орбите на протяжении 350 оборотов.
Из рис.6 видно, что в целом зависимость средней величины импульса от места совершения импульса сохраняет свой характер для различных Z, однако говорить об универсальном месте совершения маневра нельзя.
4.2 Определение наиболее эффективного направления исполнения корректирующего импульса
Для определения наиболее эффективного направления исполнения импульса коррекции для каждого угла а были рассчитаны траектории с заданным фиксированным направлением исполнения импульса (3: -170е < р < ieoe.
На рис. 7 представлены типичные графики зависимости максимальной величины импульса коррекции от направления исполнения импульса для орбиты с начальной координатой Z = 200000 км.
1.00Е-07
1Г.
S
а Я
I t
5s <pd s ^
ú 5 S
S
l.OOE-OS
1.00E-09
1.00E-10
-♦-alpha-140 -B-aIpha-150 -A-alpha-160 alpha-170
«g ....................Mi
-80 -30 20 70
Угол p, градусы
120
170
Рис. 7. Зависимость максимального значения импульса коррекции от направления
исполнения импульса
Т.к. угол Р, как было сказано выше, задает ось, вдоль которой совершается импульс коррекции, данные зависимости являются я-перио-дическими. Это дает возможность рассчитывать зависимость максимальной величины импульса от направления совершения импульса только для одного отрезка по р (в данном случае - отрезка, содержащего 0°).
Данные графики имеют ярко выраженные минимумы. Эти минимумы соответствуют направлению неустойчивости орбиты. Кроме того, из рис. 7
видно, что при приближении к некоторый значениям (3 максимальное значение импульса резко возрастает. Это объясняется тем, что при приближении к направлению устойчивости возможности управления орбитой ограничены и обеспечить необходимые для сбора статистики 350 оборотов не представляется возможным.
5 Заключение
В данном исследовании были разработаны алгоритмы подбора начальной скорости КА и величины импульса коррекции. Данные алгоритмы, реализованные в виде сценариев GMAT, были использованы для исследования наиболее эффективного места и направления исполнения импульса. На основе полученных данных можно утверждать, что наиболее эффективное направление исполнения импульса зависит от места его исполнения и лежит от —10 + жк до 60 + жк градусов, к Е Ж, в зависимости от амплитуды орбиты и места исполнения импульса.
Список литературы
[Ильин и др., 1998] Ильин И.С. и др. Баллистическое проектирование траекторий перелёта с орбиты искусственного спутника Земли на гало-орбиту в окрестности точки L2 системы Солнце-Земля // Препринты ИПМ им. М.В.Келдыша. 2014. № 6. 32 с. [Ильин и др., 2012] Ильин И.С., Сазонов В.В., Тучин А.Г. Построение ограниченных орбит в окрестности точки либрации L2 системы Солнце - Земля // Препринты ИПМ им. М.В.Келдыша. 2012. № 65. 28 с.
[Николаева и др., 2014] Николаева Ю.А., Аксенов С.А., Исследование неустойчивости гало-орбиты вокруг точки L2 системы Солнце-Земля // 13-я Международная конференция «Авиация и космонавтика - 2014». 17-21 ноября 2014 года. Москва. Тезисы. - СПб.: Мастерская печати, 2014 г. - 649-651 с.
[Cobos et al., 2002] Cobos, J., Masdemont, J.J., Astrodynamical Applications of Invariant Manifolds Associated with Collinear Lissajous Libration Orbits // 7th International Conference on Libration Point Orbits and Applications, Parador d'Aiguablava, Girona, Spain, 2002.
[Dunham et al., 2001] 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
[Farquhar, 1970] Farquhar, R.W., The Control and Use of Libration-Point Satellites // NASA TRR-346,1970.
[Gomez et al., 2002] Gomez, G., Marcote, 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.
[Gomez et al.] Gomez, G., Howell, K., Masdemont, J., Simo, C., Station-keeping Strategies for Translunar Libration Points Orbits // AAS 98-168
[Hechler et al., 2002] Hechler, M., Cobos, J., Herschel, Planck and Gaia Orbit Design // 7th International Conference on Libration Point Orbits and Applications, Parador d'Aiguablava, Girona, Spain, 2002.
[Simo et al., 1987] 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