МЕТОД РАСЧЕТА ВОЗМУЩЕННОЙ ТРАЕКТОРИИ ДВУХИМПУЛЬСНОГО ПЕРЕЛЕТА МЕЖДУ ГАЛО-ОРБИТОЙ В ОКРЕСТНОСТИ ТОЧКИ L2 СИСТЕМЫ СОЛНЦЕ - ЗЕМЛЯ И ОКОЛОЛУННОЙ ОРБИТОЙ Текст научной статьи по специальности «Физика»

УДК 629.78

БОТ: 10.18698/2308-6033-2020-11-2033

Метод расчета возмущенной траектории двухимпульсного перелета между гало-орбитой в окрестности точки Ь2 системы Солнце — Земля и окололунной орбитой

© Чжоу Жуи Московский авиационный институт, Москва, 125993, Россия

Предложен новый метод решения задачи расчета возмущенной траектории двухимпульсного перелета между окололунной орбитой и гало-орбитой в окрестности точки Ь2 системы Солнце — Земля. В отличие от традиционных численных методов, этот метод имеет лучшую сходимость. На всех участках траектории учитываются ускорения от сил притяжения Земли, Луны и Солнца как точечных масс и ускорение от второй зональной гармоники геопотенциала. Расчет траектории перелета сводится к решению двухточечной краевой задачи для системы обыкновенных дифференциальных уравнений. Разработанный метод основан на методе продолжения по параметру и не требует выбора начального приближения для решения краевой задачи. Приведены примеры и результаты анализа, выполненного на основе предложенного метода.

Ключевые слова: космический аппарат, двухимпульсный перелет, задача четырех тел, круговая ограниченная задача трех тел, гало-орбита, окололунная орбита, метод продолжения, проектно-баллистический анализ

Введение. Рассматривается задача расчета возмущенной траектории двухимпульсного перелета космического аппарата (КА) между окололунной орбитой и гало-орбитой в окрестности точки Ь2 системы Солнце — Земля. Это задача не новая, к ней неоднократно обращались в связи с подготовкой и реализацией космических проектов с размещением КА в окрестности коллинеарных точек либрации, например, в работах [1, 2]. Для расчета таких траекторий необходимо рассмотреть сложную динамическую модель задачи четырех тел с учетом сил притяжения Земли, Луны и Солнца на всех участках перелета, что приводит к необходимости решения достаточно сложной краевой задачи для системы обыкновенных дифференциальных уравнений.

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

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

Цель работы — применить к решению краевой задачи перелета между гало-орбитой вокруг точки либрации Ь2 системы Солнце — Земля и окололунной орбитой устойчивый метод продолжения по параметру [3-8], хорошо зарекомендовавший себя при решении задач по оптимизации траекторий космических аппаратов с малой тягой.

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

Схема перелета. Рассматриваются двухимпульсные перелеты между гало-орбитой вокруг точки либрации Ь2 системы Солнце — Земля и низкой окололунной круговой орбитой. Первый импульс скорости выдается в некоторой выбираемой точке гало-орбиты, а второй — в некоторой выбираемой точке заданной окололунной орбиты. В рамках данной математической модели учитываются притяжение Земли, Луны и Солнца как точечных масс и вторая зональная гармоника геопотенциала.

Исходными данными для вычисления траектории являются время выдачи первого импульса Ь, амплитуда А гало-орбиты, фазовый угол ф, определяющий положение КА на гало-орбите, время Т перелета (промежуток времени между моментами выдачи первого и второго импульсов скорости), высота кк и наклонение ¡к конечной круговой окололунной орбиты.

Для решения задачи требуется вычислить векторы двух импульсов скорости: первого импульса скорости А VI, обеспечивающего выведение КА на орбиту перелета к Луне, и второго импульса скорости АV2, обеспечивающего выведение КА с пролетной траектории на конечную окололунную орбиту.

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

Пусть три материальные точки (небесные тела) П1, П2 и П3 движутся под действием сил взаимного притяжения. В безразмерных

единицах массы точек П1 и П2 равны т1 и т2 соответственно; т1 >> т3; т1 + т2 = 1; расстояние между П1 и П2 составляет Ь = 1; угловая скорость системы ю = 1; постоянная тяготения О = 1; т2 = ц.

В этом случае уравнения движения точки пренебрежимо малой массы П3 в модели круговой ограниченной задачи трех тел в синодической (вращающейся) системе координат имеет вид

'х - 2 у = -их; у + 2 х = -иу;

где

, = -и,

1 ^ 1 , 2 2. 1 . и =-------г (х2 + у2) --|д(1 -ц);

1 = 7 (х + |д)2 + у2 + ,2; г2 = -у/(х -1 + д)2 + у2 + ,2.

С помощью приближенного решения третьего порядка Ричардсона для ограниченной круговой задачи трех тел [9] можно построить математическую модель начальной гало-орбиты во вращающейся (синодической) системе координат. В этой модели координаты и компоненты скорости КА на начальной гало-орбите можно вычислить по следующим соотношениям:

х = а21А2 + а22А, - Ах сов ф + (а23А1 - а24А

|соб 2ф-

+ (аз1А3 - аз2АхА, )3ф;

у = кАх ми ф + (¿21 А2 - ¿22А,2)) 2ф + (¿31 Ах3 - ¿32АХА,)) 3ф; , = 8А, соб ф + 5^21 АХА, (соб 2ф - 3) - 5 (Аъх - йЪ2АХА,)) 3ф; (1) Ах ф - 2 (3Ах2 - а24А,2 ) ) 2ф - 3 (А3Х - а32АхА,2 ) ) 3ф]; кАх сов ф + 2 (¿21А2Х - ¿22А,2)сов 2ф + 3 (¿31А3Х - ¿Ъ1 АхА,2) сов 3ф А ф - 2^21 АхА, ми 2ф + 3 (А;! - й32АхА,2)) 3ф

Ух = ША

Уу =шй

V, = ША5

где Ах, А, — амплитуды гало-орбиты в направлении осей X и 2 соответственно; а21, а22, а2э, а24, аэ1, аэ2, ¿21, ¿22, ¿31, ¿32, йл, йц, йЪ2, ша, — коэффициенты, являющиеся константами; ф — фазовый угол, опре-

деляющий положение КА на гало-орбите; 5 — параметр направления гало-орбиты, 5 = 1 или -1 (значение «1» соответствует северной гало-орбите, значение «-1» — южной).

Для того чтобы использовать эти решения в рамках модели возмущенных эфемерид, необходимо перейти из вращающейся системы координат, связанной с барицентром системы или точкой либрации, в фиксированную геоцентрическую экваториальную систему координат ЕМЕ2000, отнесенную к среднему экватору и среднему равноденствию в фундаментальную эпоху 12000 — 12 часов 1 января 2000 г. по шкале земного времени (ТТ).

Для этого записываем единичные векторы осей мгновенной синодической системы координат в геоцентрической экваториальной системе координат ЕМЕ2000 с использованием векторов мгновенного положения и мгновенной скорости Солнца, полученных с помощью эфемеридного обеспечения БЕ405 [10], следующим образом:

е х =■

ег =

^ х у 5 ^ х у

е у = е х х е z,

Ш =■

15 х у 5

Хс

где х5 — вектор мгновенного геоцентрического положения Солнца; у5 — вектор мгновенной геоцентрической скорости Солнца; ю — мгновенная угловая скорость движения Солнца относительно Земли.

После этого векторы положения х и скорости у КА относительно Земли в инерциальной геоцентрической экваториальной системе координат ЕМЕ2000 могут быть вычислены по безразмерным векторам положения х8;поа и скорости у8;поа в геоцентрической вращающейся системе координат в следующем виде:

Г х >

V у у

(

0 0 0

vшeу -Шех 0 ех еу еху

х

sinod

Л

,|у5 |

(2)

В промежутке времени между двумя импульсами скорости на КА действуют только внешние силы. В рассматриваемой математической модели движения КА учитываются ускорения от сил притяжения Земли, Луны и Солнца как материальных точек и возмущающее ускорение от второй зональной гармоники геопотенциала. Система дифференциальных уравнений движения КА в инерциальной геоцентрической экваториальной системе координат ЕМЕ2000 имеет вид

ёх

йг

• = у;

йу = дЕ

~Г= Т Х + а J 2 + Мм ёг г

хм х

хм х

м

м

+ ^5

^ - х

Хр - х

Г

(3)

где х = (х, у, _)Т — вектор геоцентрического удаления КА; г — время (в шкале всемирного координированного времени ЦТС); у = (ух, уу, у2)Т — вектор геоцентрической скорости КА; дЕ, , — гравитационные параметры Земли, Луны и Солнца соответственно; г = |х| — геоцентрическое удаление КА; х5 — вектор геоцентрического положения Солнца; ал — вектор возмущающего ускорения, обусловленного второй зональной гармоникой геопотенциала,

(( _2 Л Л х

а л 2 =

1 К

5 4 -1

V г у ( 2 Л

5 ^ -1

V г у

у

( _ 2

Л

- 3

V V

у у

(Л2 — коэффициент второй зональной гармоники геопотенциала; Ке — экваториальный радиус Земли); хм — вектор геоцентрического положения Луны; гм = |хм| — геоцентрическое удаление Луны; г5 = |х5| — геоцентрическое удаление Солнца.

В конце траектории перелета КА должен выйти на орбиту вокруг Луны. Для анализа селеноцентрического движения КА используется инерциальная селеноцентрическая селеноэкваториальная система координат ЬМЕ2000. Таким образом, необходимо перейти из инерци-альной экваториальной геоцентрической системы координат ЕМЕ2000 в инерциальную селеноцентрическую селеноэкваториальную систему координат ЬМЕ2000.

Единичные векторы осей системы координат ЬМЕ2000 в системе координат ЕМЕ2000 записывают в виде

е2= (сова соб5 .

Бта соб5 ,

sin5 )

е х =

е_ х е3

Iе_ х ез|

еу = ег х ех,

где ар, 5р — прямое восхождение и склонение вектора в направлении северного полюса Луны в системе координат ЕМЕ2000, которое рассчитывается по методике Международного астрономического союза [11]; е3 = (0, 0, 1) — единичный вектор в направлении оси 2 системы координат ЕМЕ2000.

Матрица перехода от селеноцентрических координат КА в ЬМЕ2000 хьМЕ к геоцентрическим координатам КА в ЕМЕ2000 хЕМЕ имеет вид

М = (е х,е у ,е _ )

хЕМЕ МхЬМЕ>

кЬМЕ

М Т х

ЕМЕ-

Для обеспечения возможности выведения КА с подлетной траектории на круговую орбиту вокруг Луны с заданной высотой и заданным наклонением к экватору Луны 4 без изменения плоскости орбитального движения при выполнении маневра необходимо соблюдение следующих трех условий в момент tk (Д = Ь + Т [3]:

г

т

-(+ \ ) = 0;

хт V т = 0; (4)

т т

С

mz

Ст

- соб 4 = 0.

Здесь гт = |хт| — селеноцентрическое удаление КА; хт = Мт (х - хм) — вектор положения КА в инерциальной селеноцентрической селеноэк-ваториальной системе координат ЬМЕ2000; ут = М (у - ум) — вектор селеноцентрической скорости КА в инерциальной селеноцентрической селеноэкваториальной системе координат ЬМЕ2000; ст = |ст|, где Г \

ст хт Х ут

С

тх

С

ту

С

V /

векторная постоянная площадей селеноцен-

трического движения КА в момент tk.

Для того чтобы решить краевую задачу по расчету двухимпульс-ного перелета, необходимо найти неизвестный вектор первого импульса скорости ДУь изменяющий вектор скорости КА в начальный момент времени ^ таким образом, чтобы после интегрирования дифференциальных уравнений движения КА (3) в конечный момент ^ удовлетворялись конечные условия (4). Начальный фазовый вектор КА в момент ^ до приложения первого импульса скорости определяется в безразмерной синодической системе координат последними тремя уравнениями системы уравнений (1) и переводится в ЕМЕ2000 с помощью преобразования (2).

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

ДУг = |ут (Ч )|-

Мм

]Км + К

где ут (к ) = Мт (у (tk ) - ум (tк )) — селеноцентрическая скорость КА в ЬМЕ2000.

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

приближения рассматривается траектория перелета с гало-орбиты к Луне, лежащая в плоскости орбиты Луны. Для этого требуется выполнение условия

Ит • ГКА = 0 (5)

где Ит = гт1 х гт2 (гт1, гт2 — векторы геоцентрического положения Луны в системы координат ЕМЕ2000 в начальный и конечный момент времени соответственно); гКА = Г(ф) — вектор положения КА на гало-орбите, определяемый первыми тремя уравнениями системы (1) и преобразованием (2).

Решив нелинейное уравнение (5), можно определить фазовый угол ф КА на гало-орбите в момент пересечения плоскости орбиты Луны. Очевидно, что при достаточно большой амплитуде гало-орбиты таких решений будет два — по одному на нисходящей и восходящей ветви гало-орбиты. В дальнейшем будем различать эти решения с помощью параметра п, равного 1 или 2 соответственно.

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

К™ + егКА | = 0 (6)

Гт2

где егт = . . — единичный геоцентрический вектор направления на

|Гт21

Луну в системе координат ЕМЕ2000 в конечный момент времени tk = Ь + Т (^ — начальный момент времени; Т — длительность перелета); егКА = . — единичный геоцентрический вектор направления

|ГКА I

на КА в системе координат ЕМЕ2000 в начальный момент времени

Пусть гт и гКА — векторы положения Луны и КА в инерциальной геоцентрического системе координат ЕМЕ2000, Я т и Яка модули векторов гт и гКА. Очевидно, что егт является почти периодической функцией от tk = Ь + Т, а егКА — почти периодической функцией от При фиксированном значении начального времени ^ существует множество значений длительности Т перелета, удовлетворяющих уравнению (6) и отличающихся друг от друга примерно на величину 2л:/(шт - шД где шт, ш,5 — средние угловые скорости геоцентрического движения Луны и Солнца соответственно.

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

ау1 =±у ^-у

КА'

где V =

, 2^E rm

|гка|( (гКА + Ы)

; ev = h, h = hm xr^; vka — скорость КА h

на гало-орбите непосредственно перед выдачей первого импульса скорости.

Схема выбора начального приближения АРЦ представлена на рис. 1. Обозначения: гт1, гт2 — векторы положения Луны в начальный и конечный момент времени; гКА — вектор положения КА; еА^ — вектор направления импульса скорости; еу — вектор направления скорости КА после выдачи АУ1; е^Кд — вектор направления скорости КА.

гт\

Гало-орбита вокруг точки либрации

h/n rmi • rm2

Тт2

Земля

Орбита Луны

lA^il • едК] /

ГКА hm•Гка = о

' I^KaI ■ еКкд

Рис. 1. Выбор начального приближения А V

Поскольку угол ime между плоскостями эклиптики и орбиты Луны ненулевой и равен примерно 5°, можно вычислить минимальную амплитуду гало-орбиты в направлении оси Z:

AZ min = XL2 tg (*me ) •

Рассматриваемое начальное приближение существует, если амплитуда AZ гало-орбиты в направлении оси Z не меньше AZmin. При AZ > AZmin, уравнение (5) имеет два решения.

Для решения краевой задачи используется метод продолжения по параметру [4, 6, 7]. Суть этого метода — формальная редукция рассматриваемой краевой задачи к задаче Коши. Краевую задачу для системы дифференциальных уравнений (3) с конечными условиями (4) можно представить в виде уравнения для невязок на правом конце траектории:

Метод расчета возмущенной траектории двухимпульсного перелета...

f (z) = 0, (7)

где вектор f состоит из правых частей конечных условий (4); z = AVi = = (Av1x, Av1y, Av1z) — вектор неизвестных параметров краевой задачи.

При некотором начальном приближении для неизвестных параметров краевой задачи z0 вычисляется вектор невязок (7):

f (zo) = b. (8)

Рассмотрим погружение уравнения (7) в однопараметрическое семейство

f (z) = (1 -т) b, (9)

где т — параметр продолжения.

Представим вектор z в виде функции от этого параметра: z = z(t), причем z(0) = z0 из уравнения (8). Потребуем выполнения равенства (9) для любого 0 < т < 1. Очевидно, что при т = 0 уравнение (9) совпадает с формулой (8), а при т = 1 — с уравнением для невязок для искомой краевой задачи (7).

Уравнение (9) фактически представляет собой ньютоновскую го-мотопию между системой уравнений f(z) - b = 0 с известным решением z = z0 и исходной системой уравнений f(z) = 0.

Дифференцируя уравнение (9) по параметру продолжения т и разрешая полученное выражение относительно производной dz/dr, получим формальную редукцию уравнения (7) к задаче Коши:

dz Y1 b

^т" lez J ; (10)

z (0) _ z 0, 0 <т<1.

Очевидно, что, интегрируя выражение (10) по т от 0 до 1, в силу уравнений (8) и (9) можно определить искомый вектор неизвестных параметров краевой задачи (7) в виде z = z(1). Фактически для рассматриваемой краевой задачи выражение (10) реализует продолжение по конечным краевым условиям от начального решения, соответствующего начальному приближению для неизвестных параметров краевой задачи, до искомого решения.

Дифференциальные уравнения (10) будем называть дифференциальными уравнениями метода продолжения. Метод продолжения можно применить, если существует матрица частных производных Sf/Sz, невырожденная на всем интервале продолжения т е [0; 1]. Чтобы улучшить сходимость процесса решения рассматриваемой краевой задачи, целесообразно ввести параметр продолжения т в правые

части дифференциальных уравнений движения и в краевые условия. В этом случае функция невязок Г явно зависит от параметра продолжения т, а дифференциальные уравнения метода продолжения примут вид

— = -{-1 йг г

В рассматриваемой краевой задаче параметр продолжения вводится в краевые условия следующим образом:

^(г, т) = тГ (г) = 0,

а в дифференциальные уравнения движения — так:

йх

йг

йу

йг

■ = у;

Д Е

= - ¿-Е х + т

(

aJ 2 + Мм

Л

хМ х

v iхм х|

М

М

(

+ Д5

хе -х

VIх 5 -х1

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

йГ = 1т(Г (г + /е^л)) + 2 ).

И

Г = 1т(ТТ + /И)) + о(И2). йт И

И = 10

-20

где ек — единичный вектор, у которого к-я компонента равна 1 (к = 1, 2, 3), а остальные компоненты — нулевые.

Численные примеры. Рассмотрим расчет траектории перелета с южной гало-орбиты вокруг точки либрации Ь2 системы Солнце — Земля, имеющей амплитуду 650 000 км, на круговую окололунную орбиту высотой 100 км с наклонением 90° при времени начала перелета г0 = 00:00:00 06.07.2020 и длительности перелета меньше 45 сут.

В этом случае, применяя процедуру, описанную в разделе «Метод решения краевой задачи», получим для п = 1 время перелета Т = 38,0896 сут. На рис. 2 приведена полученная траектория перелета между гало-орбитой вокруг точки либрации Ь2 системы Солнце —

Земля и окололунной орбитой (единица длины — радиус орбиты Луны, начало координат — Земля). На рис. 3 в крупном масштабе показан участок стыковки между переходной и окололунной орбитами.

Рис. 2. Траектория перелета между гало-орбитой вокруг точки либрации Ь2

и окололунной орбитой

Рис. 3. Участок стыковки переходной и окололунной орбит

В таблице приведена зависимость требуемых импульсов скорости от длительности перелета для перелетов в плоскости орбиты Луны при Ь = 00:00:00 06.07.2020, Л2 = 650 000 км, п = 1 и п = 2.

Из таблицы видно, что минимальный суммарный импульс скорости достигается в случае п = 1, Те [25, 45] сут. Поэтому последующие расчеты траектории с минимальной величиной суммарного импульса скорости при разных датах старта, разных значениях амплитуды гало-орбиты и начального фазового угла КА на гало-орбите выполняются для данного случая.

Зависимость импульсов скорости от длительности перелета для перелетов в плоскости орбиты Луны

Т, сут △ У\, км/с △ У2, км/с АУ%, км/с

п = 1

10,71193922 1,7953 1,5831 3,3784

38,06959355 0,3352 0,6740 1,0093

65,39604632 0,5378 0,6812 1,2190

п = 2

16,42446535 1,3100 1,0001 2,3101

43,79707195 0,6822 1,6090 2,2912

т, сут

^ ^ ^ ^ ¿Р

\ 'Ь Ъ л

Дата а

Л*сумма> КМ/МС

Дата б

Рис. 4. Зависимость времени перелета (а) и суммарного импульса скорости (б) от даты старта в июле 2020 г.

На рис. 4 приведены графики зависимости времени перелета и суммарного импульса скорости от даты старта в июле 2020 г. (п = 1, Т е [25, 45] сут).

Из данных графиков видно, что время перелета периодически изменяется с датой старта. При этом чем короче время полета, тем больше суммарный импульс скорости. Минимальный суммарный импульс скорости достигается при ^ = 00:00:00 02.08.2020. Поэтому последующие расчеты выполнялись для данного случая.

Был рассмотрен диапазон изменений амплитуды начальной гало-орбиты от -800 000 до 800 000 км (отрицательная величина амплитуды соответствует северной гало-орбите, а положительная — южной). Зависимость времени перелета и суммарного импульса скорости от амплитуды начальной гало-орбиты приведены на рис. 5.

Северные гало-орбиты

Южные гало-орбиты

-10 -8 -6 -4 -2 0 2 4 Амплитуда, 105 км а

Северные гало-орбиты

Южные гало-орбиты

сумма

-10 -8 -6 -4 -2 0 2 4 Амплитуда, 105 км б

Рис. 5. Зависимость времени перелета (а) и суммарного импульса скорости (б) от амплитуды начальной гало-орбиты

Из графиков видно, что в этих условиях чем больше время полета и чем меньше амплитуда гало-орбиты, тем меньше суммарный импульс скорости. А минимум суммарного импульса скорости достигается при п = 1, ^ = 00:00:00 02.08.2020, Az = -160 000 км, T = 41,4642 сут. Поэтому расчеты оптимального начального положения КА на начальной гало-орбите выполнялись для данного случая.

На рис. 6 показаны графики зависимости суммарного импульса скорости от начального фазового угла КА на начальной гало-орбите. Из этих графиков видно, что оптимальный начальный фазовый угол КА на начальной гало-орбите равен 290,8448°, при этом значение суммарного импульса скорости минимально и составляет 795,6 м/с (из которых AV1 = 128,5 м/с, ЛV2 = 667,1 м/с).

AV, км/с

сумма

260 265 270 275 280 285 290 ф, град

Рис. 6. Зависимость суммарного импульса скорости от начального фазового угла ф космического аппарата на начальной гало-орбите

Заключение. Разработан численный метод для расчетов двух-импульсных траекторий перелета между гало-орбитой вокруг точки либрации L2 системы Солнце — Земля и окололунной орбитой. Показаны эффективность, хорошая сходимость и вычислительная устойчивость данного метода при расчете траекторий для разных дат старта и разных значений длительности перелета.

ЛИТЕРАТУРА

[1] Howell K.C., Kakoi M. Transfers between the Earth-Moon and Sun-Earth systems using manifolds and transit orbits. Acta Astronautica, 2006, vol. 59, pp. 367-380.

[2] Li Mingtao, Zheng Jianhua, Yu Xizheng, Gao Dong. Two impulses transfer trajectory design for Sun-Earth libration point missions. J. Beijing University of Aeronautics and Astronautics, 2009, vol. 35, no. 7, рр. 865-868.

[3] Петухов В.Г., Чжоу Жуи. Расчет возмущенной импульсной траектории перелета между околоземной и окололунной орбитами методом продолжения по параметру. Вестник МАИ, 2019, т. 26, № 2, с. 155-165.

[4] Petukhov V.G. One Numerical Method to Calculate Optimal Power-Limited Trajectories. IEPC-95-221, Moscow, 1995.

[5] Петухов В.Г. Оптимизация многовитковых перелетов между некомпланарными эллиптическими орбитами. Космич. исслед., 2008, т. 42, № 3, с. 260-279.

[6] Петухов В.Г. Оптимизация межпланетных траекторий космических аппаратов с идеально-регулируемым двигателем методом продолжения. Космич. исслед., 2008, т. 46, № 3, с. 224-237.

[7] Петухов В.Г. Метод продолжения для оптимизации межпланетных траекторий с малой тягой. Космич. исслед., 2012, т. 50, № 3, с. 258-270.

[8] Николичев И.А. Оптимизация многовитковых межорбитальных перелетов с двигателями малой тяги. Вестник МАИ, 2013, т. 20, № 5, с. 66-76.

[9] Richardson D.L. Halo Orbit Formulation for the ISEE-3 Mission. J. Guidance and СстШ, 1980, vol. 3, no. 6, pp. 543-548.

[10] Ivanyukhin A.V., Petukhov V.G. Low-Energy Sub-Optimal Low-Thrust Trajectories to Libration Points and Halo-Orbits. Cosmic Research, 2019, vol. 57, no. 5, pp. 378-388.

[11] Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009. Celestial Mechanics and Dynamical Astronomy, 2011, vol. 109, pp. 101-135.

Статья поступила в редакцию 23.03.2020

Ссылку на эту статью просим оформлять следующим образом:

Чжоу Жуи. Метод расчета возмущенной траектории двухимпульсного перелета между гало-орбитой в окрестности точки L2 системы Солнце — Земля и окололунной орбитой. Инженерный журнал: наука и инновации, 2020, вып. 11. http://dx.doi.org/10.18698/2308-6033-2020-11-2033

Статья подготовлена по материалам доклада, представленного на ХЫУАкадемических чтениях по космонавтике, посвященных памяти академика С.П. Королёва и других выдающихся отечественных ученых — пионеров освоения космического пространства «Королёвские чтения — 2020», Москва, 29-31 января 2020 г.

Чжоу Жуи — аспирант Московского авиационного института. Область деятельности и научных интересов: динамика, баллистика, управление движением летательных аппаратов. e-mail: rui.zhou0225@qq.com

Method for calculating the perturbed trajectory of a two-impulse flight between the halo orbit in the vicinity of the L2 point of the Sun — Earth system and the near-lunar orbit

© Zhou Rui Moscow Aviation Institute, Moscow, 125080, Russia

The paper introduces a new method for solving the problem of calculating the perturbed trajectory of a two-impulse flight between a near-lunar orbit and a halo orbit in the vicinity of the L2 point of the Sun — Earth system. Unlike traditional numerical methods, this method has better convergence. Accelerations from the gravitational forces of the Earth, the Moon and the Sun as point masses and acceleration from the second zonal harmonic of the geopotential are taken into account at all sections of the trajectory. The calculation of the flight path is reduced to solving a two-point boundary value problem for a system of ordinary differential equations. The developed method is based on the parameter continuation method and does not require the choice of an initial approximation for solving the boundary value problem. The last section of the paper provides examples and results of the analysis based on this method.

Keywords: spacecraft, two-impulse flight, four-body problem, circular restricted three-body problem, halo orbit, near-lunar orbit, continuation method, design-ballistic analysis

REFERENCES

[1] Howell K.C., Kakoi M. Transfers between the Earth — Moon and Sun — Earth systems using manifolds and transit orbits. Acta Astronautica, July-September 2006, vol. 59, pp. 367-380.

[2] Li Mingtao, Zheng Jianhua, Yu Xizheng, Gao Dong. Two impulses transfer trajectory design for Sun — Earth libration point missions. Journal of Beijing University of Aeronautics and Astronautics, July 2009, vol. 35, no. 7, pp. 865-868.

[3] Petukhov V.G., Zhou R. Vestnik MAI — Airospace MAI Journal, 2019, vol. 26, no. 2, pp. 155-165.

[4] Petukhov V.G. One Numerical Method to Calculate Optimal Power-Limited Trajectories. IEPC-95-221, Moscow, 1995.

[5] Petukhov V.G. Kosmicheskie issledovaniya — Cosmic Research, 2008, vol. 42, no. 3, pp. 260-279.

[6] Petukhov V.G. Kosmicheskie issledovaniya — Cosmic Research, 2008, vol. 46, no. 3, pp. 224-237.

[7] Petukhov V.G. Kosmicheskie issledovaniya — Cosmic Research, 2012, vol. 50, no. 3, pp. 258-270.

[8] Nikolichev I.A. Vestnik MAI — Airospace MAI Journal, 2013, vol. 20, no. 5, pp. 66-76.

[9] Richardson D.L. Halo Orbit Formulation for the ISEE-3 Mission. J. Guidance and ^ntrol, 1980, vol. 3, no. 6, pp. 543-548.

[10] Ivanyukhin A.V., Petukhov V.G. Low-Energy Sub-Optimal Low-Thrust Trajectories to Libration Points and Halo-Orbits. Cosmic Research, 2019, vol. 57, no. 5, pp. 378-388.

[11] Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009. Celestial Mechanics and Dynamical Astronomy, 2011, vol. 109, pp. 101-135.

Zhou Rui, post-graduate student, Moscow Aviation Institute. Research interests: dynamics, astrodynamics, spacecraft motion control. e-mail: rui.zhou0225@qq.com