Исследование оптимальных пространственных траекторий перелета с малой тягой в системе Земля — Луна Текст научной статьи по специальности «Механика и машиностроение»

УДК 629.787

EDN КОШКС

Исследование оптимальных пространственных траекторий перелета с малой тягой в системе Земля — Луна

© Е.Ю. Кувшинова АО ГНЦ «Центр Келдыша», Москва, 125438, Российская Федерация

Рассмотрен пространственный перелет космического аппарата с электроракетной двигательной установкой между низкими окололунной и околоземной орбитами за минимальное время. Цель статьи — оценка влияния на затраты характеристической скорости, необходимой на перелет между окололунной и околоземной орбитами, таких баллистических параметров траектории, как угловая дальность перелета (истинная долгота) и положение плоскости орбит у Земли и у Луны. Направление вектора тяги электроракетной двигательной установки определялось с использованием принципа максимума Понтрягина. Оптимизация траекторий перелета в системе Земля — Луна в части определения направления вектора тяги электроракетной двигательной установки проводилась в рамках единой постановки («сквозной» расчет траектории) с учетом гравитации Земли и Луны (без использования метода сфер действия), положение которых определялось согласно эфемеридной модели ЕРМ 2008.

Ключевые слова: космический аппарат, оптимальные траектории перелета, электроракетная двигательная установка, малая тяга, перелет в системе Земля — Луна

Введение. Численные решения задач перелета космических аппаратов (КА) между орбитами искусственного спутника Земли и Луны приведены в литературе в ограниченном количестве. При этом вопросы влияния баллистических параметров траектории (в частности, угловой дальности перелета, положения плоскости орбит у Земли и у Луны) на оценку затрат характеристической скорости не получили достаточного освещения в имеющихся материалах. В силу этих причин возникла потребность в проведении исследований. С применением разработанной методики расчета оптимальных пространственных траекторий перелета КА с электроракетной двигательной установкой (ЭРДУ) получены результаты оценки затрат характеристической скорости, требующейся на перелет между окололунной и околоземной орбитами.

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

Работы с указанием используемого метода, особенностей решения и результатов оценки затрат характеристической скорости (Кхар), необходимой на перелет

Работа Особенности метода решения ^хар, м/с

[1, 2] Низкие окололунная и околоземная орбиты: • метод сфер действия; • оптимальный закон управления (принцип максимума) с пассивными участками на этапе согласования 79551

[3-5] • метод сфер действия; • оптимизация траектории между сферами действия Земли и Луны (задача трех тел; комбинированный метод — прямой и непрямой) с пассивными участками 78102

[6] • метод сфер действия; • оптимизация траектории на втором участке (прямой метод); • задача четырех тел; • с пассивными участками 99023

[7] • квазиоптимальное управление с обратной связью (метод КОУСОС [8]) с пассивными участками; • задача четырех тел (с использованием промежуточной орбиты в районе точки либрации Ь\); • для определения координат Луны и Солнца используются 1РЬ-эфемериды ББ403 97834

[9] Высокие окололунная и околоземная орбиты:* • оптимальный закон управления (принцип максимума) с параллельным использованием для расчета траекторий 400 процессоров; • задача трех тел; • метод разработан для трехмерной задачи, расчетный пример — для плоской задачи > 20005

[10] • оптимальный закон управления (принцип максимума) с решением задачи компланарного перелета; • задача четырех тел 312031806

[11] • оптимальный закон управления (принцип максимума) при условии прохождения через точку либрации; • задача четырех тел (эфемеридная модель движения планет) ~2970/ 285629617

Примечание. Расчеты проводились при следующих исходных данных: 1 — околоземная орбита — круговая высотой 800 км, окололунная орбита — круговая высотой 100 км, начальное ускорение — 7,75-Ю-4 м/с2, удельный импульс тяги ЭРДУ — 58 840 м/с; 2 — околоземная орбита — круговая высотой 407 км и наклонением 0, окололунная орбита — полярная круговая высотой 100 км [5]; 3 — околоземная орбита — круговая высотой 400 км и наклонением 5,2°, окололунная орбита — круговая полярная высотой 100 км; 4 — околоземная орбита — высотой ~500 км и наклонением 51,5°, окололунная орбита — полярная высотой 300 км; 5 — околоземная орбита — круговая высотой 125 000 км, окололунная орбита — круговая высотой 50 000 км; 6 — околоземная орбита — круговая высотой 167 км, окололунная орбита — круговая высотой 5000 км; 7 — околоземная орбита — эллиптическая с высотой перигея, равной 4500 км, высотой апогея — 50 000 км, наклонением — 28°, окололунная орбита — круговая высотой 5000 км и наклонением 35° / эллиптическая орбита — высотой периселения 6000 км и высотой апоселения 10 000 км [11]. * За исключением работы [10], где высота круговой околоземной орбиты равна 167 км.

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

Цель работы — оценка влияния на затраты характеристической скорости, необходимой на перелет между низкими окололунной и околоземной орбитами, различных баллистических параметров траектории (угловой дальности, положения плоскостей орбит у Земли и у Луны), полученные на основе используемой методики расчета пространственных траекторий перелета КА с ЭРДУ между окололунной и околоземной орбитами. Следует отметить, что методика усовершенствована и является продолжением исследований, приведенных в [12].

Основные методические положения. При исследовании оптимальных пространственных траекторий перелета с малой тягой в системе Земля — Луна были приняты шесть методических положений.

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

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

Проведение «сквозного» расчета траектории перелета. Отсутствие в методике упрощающих положений, предусматривающих разбиение траектории на сферы действия, обеспечило повышение точности результатов оценки затрат характеристической скорости. «Сквозной» расчет траектории перелета между окололунной и околоземной орбитами позволяет исключить влияние на оценку затрат характеристической скорости методических погрешностей, связанных с разбиением траектории перелета на участки с разным составом внешних воздействий, влияющих на КА, и, следовательно, с определением оптимальных параметров орбиты перехода между участками траектории перелета.

Использование равноденственных элементов. Для интегрирования уравнений движения применялись равноденственные элементы {h, ex, ey, ix, i , F} [1], связанные с традиционно используемыми оскулирующими элементами следующими соотношениями: i Р

h = .—, ex = e cos(Q + ro), e = e sin(Q + ro), \ д

ix = tg-2- • cos Q, iy = tg-2 • sin Q, F = Q + ro +

где p — фокальный параметр орбиты; д — гравитационный параметр; e — эксцентриситет орбиты; Q — долгота восходящего узла; ю — аргумент перицентра; i — наклонение орбиты; F — угловая дальность (истинная долгота); $ — истинная аномалия.

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

Применение двух систем орбитальных элементов для расчета траектории перелета. Для того чтобы исключить возможное вырождение уравнений движения в равноденственных элементах (вблизи другого притягивающего центра при расчете траектории только в одной из двух систем орбитальных элементов — селеноцентрической или геоцентрической), расчет траектории проводился сначала в селеноцентрической системе орбитальных элементов, а затем — в геоцентрической (переход из селеноцентрической в геоцентрическую систему орбитальных элементов осуществлялся тогда, когда эксцентриситет орбиты относительно Луны становился больше 2). Для перевода сопряженных переменных из селеноцентрической в геоцентрическую систему координат были использованы канонические преобразования, представленные в работе [13].

Использование вспомогательных декартовых систем координат. В качестве вспомогательной рассматривалась декартова система координат для формул перехода между селеноцентрической и геоцентрической системами координат, а также для определения возмущающих ускорений. Выбранная геоцентрическая система координат OXYZ является прямоугольной экваториальной с началом координат (точка О) в центре Земли и следующим направлением осей: ось ОХ лежит в плоскости экватора и направлена в точку весеннего равноденствия, ось OZ направлена перпендикулярно плоскости экватора к Северному полюсу мира, ось OY проходит в плоскости экватора перпендикулярно оси ОХ, дополняя систему координат до правой. Селеноцентрическая система координат OxXyZ была получена параллельным сдвигом

геоцентрической системы координат в центр Луны, за опорную плоскость селеноцентрической системы координат принята плоскость экватора Земли.

Определение эфемерид Луны согласно модели ЕРМ 2008, разработанной «Институтом прикладной астрономии РАН» [14].

Математическая модель движения КА. Система уравнений для фазовых переменных в равноденственных элементах имеет вид

^ = И2Т к, ё т

dex

d х

dey

d х

dix

= h(t( (1 + k)cosF + кех ) + SsinF + Weyr\Kj, = h(t( (1 + K)sinF + кеу )-ScosF -Жехг\к),

1

—x = — hwkф cos F, d х 2

diy 1

—— = — hWкф sin F, d х 2

dF_ d х dm d х

dt d х

—^r - hWqK, h3K2

2 N ЛтЛспу

L

ЭРДУ

= 1,

где h, ex, ey, ix, iy, F, m — фазовые переменные (m — масса КА); N — электрическая мощность, подаваемая на ЭРДУ; T, S, W — проекции вектора ускорения а от тяги ЭРДУ и возмущающего ускорения на трансверсальное направление, радиальное направление и нормальное направление к плоскости орбиты соответственно;

к = У(1 + ех cosF + еу sinf); л = ix sinF -iy cosF; ф = 1 + ix2 +iy2;

Лт — тяговый коэффициент полезного действия ЭРДУ; ЛСПУ — коэффициент полезного действия системы преобразования и управления; 1эрдУ — удельный импульс тяги ЭРДУ.

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

2Nлтлспу R

а = ■

LЭРДУ 'm

-Ro + Aai.

2

С соб у соб а ^

Здесь Е0 =

v

Б1п у соб а Бт а

вектор, содержащим углы а и у, рассмат-

у

риваемые в качестве параметров управления (а — угол между вектором тяги и орбитальной плоскостью, у — угол между проекцией вектора тяги на орбитальную плоскость и трансверсальным направлением); А — матрица перехода от декартовой системы координат

к орбитальной системе координат, А =

С а11 а12 а13 "

а21 а22 а23

V а31 а32 а33 у

где

ап =

а12 = собР -

б1пР + Фгу (собР• 1Х + б1пР• гу );

Ф

2

ф

1х (Р • 1х+^ • *у);

а13 =■

ф

(собР • 1х + б1пР • 1у );

а21 =

собР +—(пР• ¡т - СОБР• );

Ф п х у>

а22 = Бт Р - -

Ф

1х • 1х- со8 Р • 1у);

а23 = — (б1пР• 1х - собР• гу); Ф

аз1 = -- гу;

2

Ф

2. Ф

а32 = гх; а33 =

(1 --1

Ф

<2).

а1 = |аз, ал ] — возмущающие гравитационные ускорения, действующие на КА, в геоцентрической и селеноцентрической декартовой системах координат соответственно.

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

Рассматривалась система векторов (рис. 1), определяющая положение трех тел (Земля — Луна — КА) в инерциальной системе координат. На рисунке векторы с началом в центре инерциальной системы координат Ои обозначены р, векторы с началом в центре Земли — Е, векторы

с началом в центре Луны — Г, векторы с окончанием в центре Земли — индексом «З», векторы с окончанием в центре Луны — индексом «Л», векторы с окончанием в центре масс КА — без индекса.

КА

Земля

Луна

Рис. 1. Система векторов для КА, Луны и Земли

Согласно второму закону Ньютона, ускорения, действующие на КА, Луну и Землю, определяются следующим образом:

•• Я Г

^ —возм

Р = _ИЗТ3 _ Ил"Г + аКА ,

Я г

(1)

Р = и -З- + авозм

рЛ - ИЗ 3 + аЛ ,

(2)

Рз = Ил

Я

я3

Л , —возм

" "Г С1о ,

(3)

где аКАзм

—возм

(Л Т

—возм

6/о

л , аЗ — возмущающие ускорения от Солнца, планет (кроме Земли и Луны) и других объектов Солнечной системы, действующие на КА, Луну и Землю соответственно.

После двойного дифференцирования соотношений из векторных треугольников «Ои — Земля — КА» и «Ои — Луна — КА» следует

Р = Рз + ял +г, Р = Рз +я.

(4)

(5)

Замена в соотношениях (4) и (5) векторов р и Рз с началом в точке отсчета Ои инерциальной системы координат ОиХиУи1и формулами (1)—(3) приводит к следующим выражениям для гравитационных ускорений Г и Я:

Г = -Цз Я - цЛ -Т- - Цд ЯД - Яд + —КГ - <3м, (6)

Я3 г3 ЯЛ

Я — . . Я ,, г . . ЯЛ I —возм — возм (7)

Я г ЯД

Поскольку расстояние между КА и Землей по астрономическим меркам относительно мало, разностью возмущающих ускорений —КА3м - —з°3м от воздействий Солнца, планет (кроме Земли и Луны) и других объектов Солнечной системы, действующих на КА и Землю, можно пренебречь. Возмущающие воздействия от Солнца и других небесных тел учитывались опосредованно через использование

высокоточных эфемерид Луны с точностью до разности —Козм - —з°зм,

которой пренебрегли.

С учетом этого соотношения (6) и (7) примут вид

,, Я г ЯЛ - /оч

г —-^ЗТГ3-1д — -1дТз -Ял, (8)

Я3 г3 ЯЛ

,, Я г ял ^

Я — -^зт;з-1д —^Л^. (9)

Я3 г3 ЯЛ

Подстановка соотношений (8) в (4) и (9) в (5) приводит к одному уравнению в инерциальной системе координат (относительно барицентра Солнечной системы):

Р Р Я г ял

Р — Р З-Д ——Д^ Я г ЯД

Таким образом, траектории, рассчитанные в системах координат относительно Земли и Луны и относительно барицентра Солнечной системы, должны совпадать.

Ч Яг й

Члены -|3 —3 и -|Л -3 соотношений, приведенных выше, являются основными гравитационными ускорениями, действующими на КА в геоцентрической и селеноцентрической системах координат, поэтому они не входят в выражения для возмущающих гравитационных уравнений. Тогда возмущающие гравитационные ускорения, действующие на КА, в геоцентрической системе координат а3

и в селеноцентрической системе координат —д принимают следующий вид:

r R

аз = -^л — (10)

r Rл

R R

ал -^л^Т -Rл • (11)

R3 ™ Rj|

Использование принципа максимума Понтрягина. Рассматривалась вариационная задача по определению законов управления вектором тяги ЭРДУ при перелете КА между круговыми окололунной и околоземной орбитами за минимальное время (задача быстродействия). Направление вектора тяги определялось по закону, полученному с использованием принципа максимума Понтрягина [15].

Поскольку рассматривалась задача быстродействия, согласно принципу максимума Понтрягина, режим работы электроракетных двигателей должен быть непрерывным (без выключений) [16]. Тяга и удельный импульс тяги ЭРДУ в процессе перелета предполагались постоянными.

В соответствии с принципом максимума Понтрягина вводятся сопряженная вектор-функция ¥ и гамильтониан Н, который для задачи быстродействия (минимального времени перелета) имеет вид

H = £Y if = Y hh2T к + i=1

+ Y eh (t ((1 + к) cos F + кех ) + S sin F + Wey ^к) + + Yeyh(t((1 + к) sinF + кеу ) - ScosF - Wexлк) -

- Yix • 2 hW кф cosF - Yiy • 1 hWкф sinF +

+ Yf ^h^ - hWЛк|- mYm + Y,,

где Y h, Y ex, Y ey, Y ix, Yiy, YF, Y m, Y t — переменные, сопряженные переменным h, ex, ey, ix, i , m, t (t — продолжительность перелета).

Углы ориентации вектора тяги ЭРДУ (а, у), согласно принципу максимума, определяются из следующих соотношений:

А л„

cosy cosa = . , siny cosa =

i , oiii y vuoa — i

Va?+A+лП VA?+A+Л

4

2

sina =

Va?+л?+л2

где

Ат—И2 к¥И + И (1 + к) 008 Р ■х¥е + И кех ¥е + И (1 + к) 81п Р ■х¥е + И кеу ■х¥е ;

Аг — И б1ПР е - ИооъР е ;

Лп — -ИцкФр + Илкеу-Илкех-1 Икфоо8Р-1 Икф81пР.

Ух у 2 х 2 у

Проекции вектора ускорения — от тяги ЭРДУ и возмущающего ускорения в селеноцентрической системе координат на трансвер-сальное направление, радиальное направление и нормальное направление к плоскости орбиты соответственно имеют вид

Т —

2 NЛтЛспу

(

I.

эрдут л/лт2 + Л2Г + А2

3-

х + х

Л

д

Я3

1 д л3 хд ЯД

—11 +

^ у+уЛ уЛ ^

-1з--3--1Д ТЛ " ул —12 +

V Я Яд

I + 1Л

-|з-3Л

3 Я3

1Л Я3 ^Л

Я

13'

где х, у, I и Хд, уд, 1д — координаты КА и Луны в декартовой селеноцентрической системе координат соответственно;

5 —

2 nлтЛспу

(

I

эрду т ^А+А+А

3-

х + х

Л

Л

Я3

1 д 7ГГ хД ЯД

—21 +

^ у+уЛ уЛ .. ^ -13--3--1Л ЯТ - уд —22 +

V Я ЯЛ )

I + 1Л

-|3-

Я3

1 дт^ гд Я

23'

ж —

2 NЛт Лспу А

I

эрду т ^А + аг2 + А

2

п V

-|3-

х + х

Л

Я3

1 д 7ГГ хД ЯД

—31 +

^ у+уЛ уЛ .. ^ |з Я3—|д тг - уд —32 +

V Я ЯД )

-|з

I + 1т

Я3

1Л гд Я

—

33 •

Аналогичным образом проекции Т, 5, Ж были получены в геоцентрической системе координат.

Наличие координат Луны, зависящих от независимой переменной (времени), делает уравнения движения неавтономными, поэтому вводится новая независимая переменная т, связанная с переменной времени

дифференциальным уравнением — = 1. Введение дополнительного диф-

ё т

ференциального уравнения и новой независимой переменной не меняет

Исследование оптимальных пространственных траекторий перелета...

форму сопряженных уравнений, лишь добавляется дифференциальное уравнение для переменной, сопряженной времени:

^ = _Г 4 &Т+а/А+Апм \

& т I & & & )

Система уравнений для сопряженных переменных имеет вид

dSl = -( +t ^^ + s A + A,. ® + W A+An 8W 1

dx l K2h4 8h 4 8h 8h r 8h 8h n 5h I

d 44

d x

d

d x

d

d x v

(

( 2 4f 8K +t 8a. + Ax 8T + S ^ + Ar 8S

l K3h3 8ex 8ex x 8ex 8ex r 8ex

( -2 ТF3 8K + T 8Ax + Ax 8T + S + Ar 8S

l K3h3 8ey 8ey x 8ey 8ey r 8ey

u/ 8An . 8W 1

+ W—n + An-

8ex 8ex

ш 8An . 8W + W—- + An-

8e„ 8e

Л

y

y

f

,8Ax

-8An

8W

d x d 4

T^- + Ax — + S8Ai + Ar — + W^ + A

8ix 8ix 8ix 8ix 8ix 8ix

T8Ax . 8T S 8Ar . 8S W8An 8W T—x + Ax — + S—- + Ar — + W—n + An

\

\

8i,.

8i,.

8iy

8iy

8iy 8iy I

d x d 4

F I . 4F 8K 8Ax 8T 8Ar 8S 8Д F = -l -2—rF--+ T—x + Ax— + S—- + Ar — + W

8W

c3h3 8F

8F

8F 8F

8F

+ An 8F n 8F

I Л 8T 8S . 8W = -| A — + Ar — + An-

d x l 8m 8m 8m

dyt (.8T .8S . 8W

—— = -l Ax — + Ar — + An-

d x l dt dt dt

Для интегрирования систем дифференциальных уравнений движения и сопряженных переменных в селеноцентрической и в геоцентрической системах координат использовался метод Рунге — Кутты — Фельберга 7(8)-го порядка с автоматическим определением шага [17-19].

Использование принципа максимума Понтрягина приводит к двухточечной краевой задаче для системы дифференциальных уравнений.

Решение краевой задачи сводится к задаче минимизации функции (не*

вязки) S , которая решается совместно методами Ньютона (поиск ну*

лей компонентов вектор-функции Sj) и градиентного спуска (минимизируется норма функции невязки:

* - ^ * 2 S = Z (аА) , i=1

неизв с* ^

где п — количество неизвестных краевой задачи; щ — весовой

* *

коэффициент; — компоненты минимизируемой функции Sj,

представляющие собой разницу между требуемыми и получаемыми в конце траектории элементами (зависят от вида граничных условий)) [20-21].

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

При перелете между окололунной и околоземной орбитами также возможно существование экстремалей, различающихся угловой дальностью перелета, как у Луны РД, так и у Земли Р . Угловая дальность Р не может учесть все свойства траектории, связанные с количеством витков у Луны Рд и у Земли Р3 . Поэтому был введен

аналог угловой дальности Р, учитывающий целое число реализованных витков относительно Луны на момент перехода из селеноцентрической системы координат в геоцентрическую и угловую дальность перелета, рассчитываемую относительно Земли, начиная с момента перехода в геоцентрическую систему координат:

Р — [ !д/ 2я] 2л + !3,

где [...] — целая часть числа; РЛ — угловая дальность селеноцентрического участка; Р3 — угловая дальность геоцентрического участка.

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

разные семейства экстремалей.

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

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

конце траектории (у Земли);

- при свободных значениях аналога угловой дальности перелета Р ) и долготы восходящего узла 0( ^ ) на правом конце траектории (у Земли);

- при свободных значениях аналога угловой дальности перелета Р ) и долготы восходящего узла 0( ) на правом конце траектории (у Земли), при свободной долготе восходящего узла 0(0 ) на левом конце траектории (у Луны).

В общем, для всех краевых задач рассматривается пространственная задача перелета КА с начальной низкой круговой окололунной орбиты на конечную низкую круговую околоземную орбиту за минимальное время. Параметры начальной орбиты — высота Н (¿0) =

= 100 км, наклонение /Д (¿0) = 90° (относительно плоскости земного экватора). Параметры конечной орбиты — высота Н (¿к) = 800 км, наклонение /3 (¿к) = 51,6°. Здесь и далее индексы «Л» и «3» обозначают селеноцентрическую и геоцентрическую системы координат соответственно, «0» и «к» — параметры начальной и конечной орбит. Начальная масса КА у Луны тЛ (¿0 ) зафиксирована.

Краевая задача межорбитального перелета при фиксированном аналоге угловой дальности перелета или со свободным аналогом угловой дальности перелета на правом конце траектории (у Земли). Для рассматриваемой краевой задачи приняты следующие дополнения параметров для начальной и конечной орбит соответственно: начальная точка лежит на оси ОХ, связанной с Луной; конечный аналог угловой дальности Р (¿к ) — зафиксирован или свободен (в зависимости от рассматриваемой постановки задачи), конечная долгота восходящего узла & (¿к ) — зафиксирована.

Неизвестными в решаемой краевой задаче являлись начальные значения сопряженных переменных (на левом конце траектории

у Луны) — , , ¥Л, ¥Л, и продолжительность перелета

Продолжительность перелета включена в состав минимизируемой функции и использована нормировка сопряженных переменных (при

условии ¥д (¿0 ) = 1), для того чтобы исключить из числа неизвестных краевой задачи одну сопряженную переменную ¥д. Начальные значения сопряженных переменных ¥т и ¥; могут быть выбраны произвольно, так как они не влияют на закон управления.

Для поиска траекторий перелета использовался следующий двух-стадийный подход.

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

' гЗ (к)- гЗкЛ

* «1 = Га (к )- Гак , (12)

К ь(к)-4 ,

где г, — радиус перицентра; tk — конечный момент времени на правом конце траектории или дата прибытия КА на целевую околоземную орбиту; га — радиус апоцентра; / — наклонение орбиты.

Такая постановка определяет незамкнутую краевую задачу, у которой количество неизвестных переменных превышает количество краевых условий, она позволяет более просто получить траекторию с заданными свойствами («раскрутку» — у Луны и «скрутку» — у Земли).

*

Краевую задачу с указанными выше компонентами невязки ^

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

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

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

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

(

¿2 =

кЗ ( 'к )- - кЗк ]

еЗ ('к )- -еЗ к

еЗ ('к)- -еЗ к

гЗ ('к )- ■ з к гх

гЗ ('к )- ■ З к ~гУ

17 ('к )- 1 к)

(13)

где кЗ

3 к 3 к -З к -З к

"У

IV

У

17 — зафиксированы. При постановке краевой задачи перелета со свободной переменной 17 ('к ) сопряженная переменная ¥ 1 () должна быть равна нулю (т. е. должны выполняться условия трансверсальности) — ¥ 1 (к ) = 0.

В этом случае компоненты невязки будут иметь вид

Г кз ('к)-кзк ^

еЗ (к)-еЗк

% =

!('к )-('к)-

У ('к )-*У ¥ 1 ('к)

З

З

х З к

У

• З к

• З к

(14)

где кк

е

З к З к З к З к

г,

У , ,х , у — зафиксированы.

Краевая задача межорбитального перелета при свободных значениях аналога угловой дальности перелета и долготы восходящего узла на правом конце траектории (у Земли). Для рассматриваемой краевой задачи дополнения параметров для начальной и конечной орбит соответственно следующие: начальная долгота восходящего узла & ('0 ) = 0,01°, конечные аналог угловой дальности

перелета 17 ('к ) и долгота восходящего узла & ('к ) являются свободными переменными. Для того чтобы получить траекторию с заданными свойствами так же, как и в предыдущей постановке задачи, требуется определить начальные значения пяти сопряженных переменных — ¥

1, ¥Л , ¥Л, ¥Л, ¥Л и продолжительности перелета '.

X У X У

Первая стадия решения краевой задачи будет аналогична постановке задачи, приведенной в предыдущем разделе при фиксирован-

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

Вторая стадия., на которой при постановке краевой задачи со свободным аналогом угловой дальности перелета Р ) и свободной долготой восходящего узла ) требуется. чтобы на правом конце траектории выполнялись также условия трансверсальности — ¥р ) = 0.

Долгота восходящего узла в конечный момент времени I = tk может быть задана следующим многообразием:

а = еЗ 2 + е3 2 - еЗк 2 = 0.

•З 2 , З 2 .2 £2 = гх +Ч -

( к \

V 2 J

= 0.

Условия трансверсальности применительно к заданному выше многообразию имеют вид

¥Зх (tk ) = 2^е/

¥Зз (tk ) = 2у1ез

¥ЗХ (tk ) = 2у2/, ¥ З (к ) = ^ 21з

е3¥3 - еЗ¥3 = 0

х ез 3 ех

/З¥З - /З¥З = 0

где у1 и у2 — коэффициенты.

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

( (tk)-hзk Л

,З (+ \ хт/З (+ \ „З (+ \ хт/З

ез З

13

3 2 ,

=

!(tk)¥% (Ч)-еЗ(tk)¥Зх (tk)

Х (tk)¥3 (tk)- 1эу (н)¥Зх (tk)

в ) + е_З2 (tk )

? (tk)+/3 2 (ч)-1 ¥ Р (Ч)

(15)

где ИЗ. 4 — зафиксированы.

Краевая задача межорбитального перелета при свободных значениях аналога угловой дальности перелета и долготы восходящего узла на правом конце траектории (у Земли), свободной долготе восходящего узла на левом конце траектории (у Луны).

Для рассматриваемой краевой задачи дополнения параметров для

начальной и конечной орбит соответственно следующие: начальная точка лежит на оси ОХ, связанной с Луной, начальная долгота восходящего узла 0('0 ) — свободная переменная; конечные аналоги угловой дальности перелета 17 ('к) и долгота восходящего узла 0('з) —

свободные переменные.

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

¥Л, ¥Л, ¥Л, ¥Л, неизвестная продолжительность перелета '

и долгота восходящего узла у Луны 0( '0 ).

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

еЗ ('к )- еЗ Рз ('к )- Рз

гЗ ('з )- гЗ к

гЗ ('к)- гЗк

=

■з ('к)-

к

■З2 ('к) + ¡?2 ('к)-гЗ

(16)

где е'З, рз, глз к, г3 к, ¡З — зафиксированы.

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

По сравнению с предыдущей постановкой краевой задачи, где только аналог угловой дальности перелета 1 ('к ) и долгота восходящего узла 0('з ) на правом конце траектории у Земли полагались свободными, в данной постановке задачи дополнительно считалась свободной еще и долгота восходящего узла 0('0) на левом конце

траектории у Луны. В связи с этим краевая задача дополнялась следующим условием трансверсальности (аналогично условиям трансверсальности при свободной переменной 0('к )):

¡Л ¥ Л - ¡Л ¥ Л = 0.

У

У ¡X

В этом случае условие ед ¥л - еЛ ¥Л = 0 не будет использоваться, так как в начальный момент времени to окололунная орбита является круговой, т. е. е = ед = ед = 0 (условие трансверсальности автоматически выполняется за счет инициализации исходных данных).

Компоненты невязки будут иметь вид

Г Ьз (з)-л

£ (з)¥е^ (3)-е) (Ч)¥\ (3)

43 (з )¥* (3)-1) (3 )¥=? (3)

еХ 2 (з) + е) 2 (3 )

Я 2 (з)+% 2 (з)-4

¥Р (3 )

й* =

Ч (to)¥Л (to)-^Ул (о)¥Л (to)

(17)

где И3 , 13 — зафиксированы.

Таким образом, рассматриваемая методика расчета траекторий перелета КА с ЭРДУ включает в себя совокупность следующих положений:

- рассматривался пространственный перелет между низкими окололунной и околоземной орбитами за минимальное время;

- определение направления вектора тяги ЭРДУ проводилось с использованием принципа максимума Понтрягина и в рамках единой постановки («сквозной» расчет траектории) с учетом гравитации З емли и Луны (без применения метода сфер действия);

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

- вводился аналог угловой дальности Р = [РЛ /2л] 2л + РЗ, позволяющий учесть все свойства траектории, связанные с количеством витков и у Луны — Рл , и у Земли — Р3.

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

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

Исследование влияния аналога угловой дальности перелета при фиксированных значениях долготы восходящего узла на левом и правом концах траектории. Получение оптимальных пространственных траекторий перелета с «раскручивающейся» спиралью у Луны и «скручивающейся» спиралью у Земли характеризуется вычислительными проблемами при низких высотах стартовой окололунной и целевой околоземной орбит, поэтому во многих методиках (например в [9]) для ее преодоления используются высокие начальные и/или конечные орбиты. В выполненном исследовании для получения оптимальных пространственных траекторий перелета между низкими окололунной и околоземной орбитами был реализован подход с применением решения вспомогательных краевых задач. При его использовании для расчета траектории перелета в равноденственных элементах начальный вектор сопряженных переменных можно определить варьированием сопряженных переменных с большим шагом от начального вектора из нулей и единиц следующего вида:

{¥Л = 1, ¥ Л = ¥ Л =¥Л =¥Л = ¥ ^ = 0}.

Таким образом, получается первое приближенное решение вспомогательной краевой задачи с компонентами невязки (12). Далее с теми же компонентами невязки реализуется метод спуска по параметру (по высоте орбиты).

Решение краевых задач пространственного перелета между окололунной и околоземной орбитами проведено со следующими исходными данными: стартовая окололунная орбита высотой 100 км; плоскость стартовой орбиты лежит в плоскости, перпендикулярной плоскости земного экватора; целевая околоземная орбита высотой 800 км и наклонением 51,6°; начальное ускорение принято равным 1,7 10-3 м/с2; удельный импульс тяги ЭРДУ составляет 29 420 м/с; дата старта с окололунной орбиты — 10 января 2020 года (такая дата старта представлена в работе [12] на примере расчета траектории перелета с фиксированным аналогом угловой дальности 17, которая является начальным приближением для настоящего исследования).

Аналог угловой дальности 1 допускает неоднозначность в задании траектории с точностью до перераспределения целого количества витков между геоцентрическим и селеноцентрическим участками. Однако практика расчетов показывает, что аналог угловой дальности 17 обладает свойством незначительно изменять угловую дальность селе-

ноцентрического участка 1Л при варьировании угловой дальности геоцентрического участка 1З. В результате решения вспомогательной краевой задачи с компонентами невязки (12) случайным образом было определено численное значение аналога угловой дальности I7, которое было зафиксировано и использовано для решения краевой задачи с компонентами невязки (13). От полученной экстремали строилась последовательность решений с перебором величины 17, т. е.

строилась зависимость ¥хар (I). Для получения зависимости Ухар (I)

использовался метод спуска по параметру I. При этом в процессе решения угловая дальность селеноцентрического участка 1Л оставалась практически неизменной, а угловая дальность геоцентрического участка I3 изменялась малыми шагами.

Для принятых исходных данных были получены оценки затрат характеристической скорости в зависимости от аналога угловой дальности I7. Зависимости величины затрат характеристической скорости ¥хар и конечные значения сопряженной переменной ^^ от угловой дальности геоцентрического участка 13 при двух фиксированных значениях долготы восходящего узла & (1к) в конечный момент времени у Земли, равных 290° и 0, представлены на рис. 2 и 3 соответственно.

Минимальные значения оценки затрат характеристической скорости, равные ~8220 м/с при значении долготы восходящего узла у Земли & (1к ) = 290° и ~ 8428 м/с — при & (^ ) = 0, достигаются

Гхар, м/с т|х10"3

Рис. 2. Зависимость ¥хар и конечного значения от 13 при долготе восходящего узла у Земли & (tk ) = 290°

^хар. М/С

Ч^хЮ"4

Рис. 3. Зависимость ¥хар и конечного значения ¥р от РЗ при долготе восходящего узла у Земли ^ (¡к ) = 0

при близких значениях угловой дальности геоцентрического участка Р3 , составляющих 52 515° и 52 540°, или ~146 витков (при этом

угловая дальность селеноцентрического участка рл ~12 240°, или 34 витка).

В точках, где оценка затрат характеристической скорости принимает минимальное значение, сопряженная переменная ¥р (¡к ) = 0, т. е. выполняется условие трансверсальности для постановки краевой задачи перелета при свободной переменной Р (¿к ). Зависимости ¥хар (РЗ ), приведенные на рис. 2 и 3, объединяет наличие одной экстремали в отличие от аналогичных зависимостей Ухар (Р) в задаче перелета на

геостационарную орбиту [22].

Таким образом, предположение о наличии только одного минимума в зависимости оценки затрат характеристической скорости от аналога угловой дальности позволяет использовать условие ¥р (д) = 0 и не задавать заранее неизвестную переменную Р (¡к), а определять

как оптимальную.

Исследование влияния долготы восходящего узла на правом конце траектории при фиксированном значении долготы восходящего узла на левом конце траектории у Луны и свободном значении аналога угловой дальности перелета. В решенных ранее краевых задачах долгота восходящего узла на правом конце траектории у Земли Q (¡к) определялась случайным образом из решения краевой задачи с компонентами невязки (13). Определение долготы восходящего

узла & (д) оптимальным образом позволит получить меньшее значение оценки затрат характеристической скорости, необходимой на перелет.

В предыдущем примере исследование влияния аналога угловой дальности на оценку затрат характеристической скорости приведено при двух значениях долготы восходящего узла у Земли (290° и 0). Ниже даны результаты оценки затрат характеристической скорости при изменении долготы восходящего узла у Земли & ) от 0 до 360°

с шагом 1°. Для этого был решен ряд краевых задач с компонентами невязки (14) при (д) = 0, где в качестве достаточно хорошего начального приближения использовалось ранее полученное решение при & ) = 290°, и с исходными данными, приведенными в предыдущем иллюстрирующем примере. Полученная зависимость оценки затрат характеристической скорости от долготы восходящего узла у Земли & (¿к ) при (¿к ) = 0 представлена на рис. 4 (данная зависимость была получена в предыдущих исследованиях и приведена в [23]).

^хар, м/с

Рис. 4. Оценка затрат характеристической скорости V в зависимости от долготы восходящего узла у Земли & (^ ) при ¥I (к ) = 0

Функция ¥хар (&)), показанная на рис. 4, при изменении долготы восходящего узла у Земли от 0 до 360° имеет два минимума и один максимум, которые соответствуют значениям долготы восходящего узла у Земли, равным 61°, 288° и 185°. При оптимальной (обеспечивающей минимальные энергетические затраты) угловой дальности геоцентрического участка I ~ 52 485° (~146 витков) минимальная оценка затрат характеристической скорости, равная 8161 м/с, достигается при

долготе восходящего узла у Земли Q (¡к) = 61°. Максимальное значение оценки затрат характеристической скорости, равное 8971 м/с, достигается при долготе восходящего узла у Земли О(¡к) = 185° (угловая дальность геоцентрического участка РЗ ~ 51 745° (~144 витка)).

Разница между минимальным и максимальным значениями оценки затрат характеристической скорости в зависимости от конечной долготы восходящего узла у Земли составила ~810 м/с (~9 % от максимального значения).

С целью определения оптимальных значений на правом конце траектории не только аналога угловой дальности, но и долготы восходящего узла была решена краевая задача с компонентами невязки (15) (в качестве начального приближения был использован вектор неизвестных переменных, полученный при Q (¡к) = 61°) со свободными

переменными Р (¡к) и Q (¡к). Для решенной краевой задачи оптимальное (соответствующее минимальному значению оценки затрат характеристической скорости) значение долготы восходящего узла у Земли составило 60,7°, оптимальное значение истинной долготы селеноцентрического участка — Рл ~ 12 317° (~34 витка), геоцентрического участка — Рл ~ 52 486° (~146 витков). При этом оценка затрат

характеристической скорости на перелет имела значение 8161 м/с.

*

Искомые значения неизвестных краевой задачи равны соответственно: = -0,0017751174; = -0,0030719190; уЛ = -0,0137890768;

ул = -0,0093318108; уР = 6,83832-10-5, X = 1,2320245878. Точность

решения краевой задачи составила несколько десятков метров от заданного значения высоты круговой орбиты 800 км. Точность по наклонению орбиты в градусах составила 1,0 • 10-4 при заданном значении 51,6°. Продолжительность перелета с окололунной орбиты на околоземную — 48,5 сут.

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

Значения неизвестных краевой задачи приведены с точностью до 10-го знака после запятой. Расчеты показали, что при уменьшении количества знаков после запятой до 8 точность решения краевой задачи снижается. В этом случае точность решения составляет несколько десятков километров от заданного значения высоты круговой орбиты, равного 800 км, по наклонению орбиты в градусах - 1,3-10-1 при заданном значении / = 51,6°.

орбиты, лежащей в плоскости, перпендикулярной плоскости земного экватора, на околоземную орбиту высотой 800 км и наклонением 51,6°; начальное ускорение принято равным ~1,710-3 м/с2, удельный импульс тяги электроракетных двигателей — 29 420 с, дата старта — 22 августа 2038 года (выбор даты обусловлен близостью оскулирую-щих элементов Луны к соответствующим оскулирующим элементам на дату 10 января 2020 года, для которой ранее были получены решения краевых задач).

Краевые задачи решались с компонентами невязки (15) при свободных значениях аналога угловой дальности II (/к) и долготы восходящего узла & (д) на правом конце траектории у Земли. Долгота восходящего узла на левом конце траектории & (/0 ) у Луны варьировалась на интервале от 0 до 360°.

В результате расчетов получена оценка затрат характеристической скорости от долготы восходящего узла у Луны & (^ ) на интервале от 0 до 360° с изменяемым шагом (рис. 5). Краевые задачи решены с высокой точностью, которая составила от 2 км до 10 м от заданного значения высоты круговой орбиты (800 км). Точность по наклонению в градусах была в диапазоне от 7,7-10- до 1,0-10- от заданного значения 51,6°.

^хар.^

Рис. 5. Оценка затрат характеристической скорости V в зависимости от долготы восходящего узла у Луны & ()

На рис. 5 видно, что зависимость Кхар (& (^)) при изменении

долготы восходящего узла у Луны от 0 до 360° имеет два минимума и два максимума, которые соответствуют значениям долготы восходящего узла у Луны, равным 37°, 200° и 115°, 300°. Минимальное

значение оценки затрат характеристической скорости, составляющее 7911 м/с, достигается при следующих оптимальных, т. е. обеспечивающих минимальные энергетические затраты, параметрах: долготе восходящего узла у Луны О (¡0) = 200°, угловой дальности геоцентрического участка РЗ ~ 52 800° (~146 витков), угловой дальности селеноцентрического участка Рл ~ 12 624° (~35 витков).

Максимальное значение оценки затрат характеристической скорости ~ 8122 м/с достигается при долготе восходящего узла у Луны О (0) = 115°, при этом РЗ ~ 52 441° (~146 витков), Рл ~ 12 610°

(~35 витков). Разница между минимальным и максимальным значениями оценки затрат характеристической скорости в зависимости от начальной долготы восходящего узла у луны составила ~210 м/с (~3 % от максимального значения). Зависимость ¥хар (О (¡0)), построенная для принятых исходных данных, в отличие от зависимости Ухар (О(¡к)), имеет два максимума, которые соответствуют значениям долготы восходящего узла у луны 115° и 300°, с разницей в оценке затрат характеристической скорости, равной 34 м/с.

С целью определения оптимальных значений долготы восходящего узла на левом и правом концах траектории, а также аналога угловой дальности на правом конце траектории была решена вспомогательная краевая задача с компонентами невязки (16), а затем основная краевая задача с компонентами невязки (17) со свободными переменными О (¡0 ), О (¡к ), Р (¡к ). Для решенной краевой задачи

были получены оптимальные (соответствующее минимальному значению оценки затрат характеристической скорости) значения долготы восходящего узла у Земли О (¡к) = 11,7°, долготы восходящего

узла у луны О (¡0 ) = 199,1°, угловой дальности селеноцентрического участка Рл ~ 12 623° (~35 витков) и геоцентрического участка РЗ ~ 52 800° (147 витков). При этом оценка затрат характеристической скорости на перелет составила 7911 м/с.

Искомые значения неизвестных краевой задачи равны соответственно = -0,0272170617; = -0,0256438466; у£ =

= 0,0516461101; ул = 0,0179191675; уР = 0,0015570103, X =

= 1,1991805976; О (0 ) = 3,4755572408. Точность решения краевой задачи составила несколько десятков метров от заданного значения высоты круговой орбиты, равной 800 км, точность по наклонению орбиты в градусах — 1,7-10-6 при заданном значении 51,6°, продолжительность перелета с окололунной орбиты на околоземную — 47,2 сут.

Таким образом, разница между минимальным и максимальным значениями оценки затрат характеристической скорости, необходимой на перелет, при варьировании долготы восходящего узла у Земли О (¡к), равная 810 м/с почти в 4 раза превышает разницу между минимальным и максимальным значениями при варьировании долготы восходящего узла у луны О (¡0 ) (210 м/с). Такая разница может быть

обусловлена массой Земли, которая больше в 81,3 раза массы луны.

Заключение. Усовершенствована методика расчета пространственных траекторий перелета КА с ЭРДУ, позволяющая с высокой точностью определять оценку затрат характеристической скорости, необходимой на перелет между низкими окололунной и околоземной орбитами. Особенностью методики является введение вместо угловой дальности Р аналога угловой дальности Р, позволяющего учесть все свойства траектории, связанные с количеством витков у луны Рл и у Земли Р3.

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

При постановке краевых задач межорбитального перелета с фиксированным аналогом угловой дальности перелета на правом конце траектории — у Земли (при фиксированных значениях долготы восходящего узла на левом и правом концах траектории) показано наличие только одного минимума в зависимости оценки затрат характеристической скорости от аналога угловой дальности, что позволило использовать постановку задачи со свободной переменной Р. Так, при значении долготы восходящего узла у Земли, равной 290°, минимальное значение оценки затрат характеристической скорости составило ~8220 м/с.

Получены оценки влияния на затраты характеристической скорости долготы восходящего узла у Земли при свободном аналоге угловой дальности перелета Р (к ) и фиксированном значении долготы восходящего узла на левом конце траектории у луны. Показано наличие двух минимумов и одного максимума в зависимости Ухар (О (¡к ) ), минимальное значение оценки затрат характеристической скорости составило 8161 м/с при значении долготы восходящего узла у Земли, равном 60,7°. Разница между минимальным и максимальным (~8970 м/с) значениями оценки затрат характеристической скорости составила ~810 м/с, или ~9 % от максимального значения.

Получены оценки влияния на затраты характеристической скорости долготы восходящего узла у Луны при свободных аналоге угловой дальности перелета F (tk ) и долготе восходящего узла на правом конце траектории у Земли. Показано наличие двух минимумов и двух максимумов в зависимости Кхар (Q (t0 )). Минимальное значение

оценки затрат характеристической скорости составило 7911 м/с при значениях долготы восходящего узла у Луны, равном 199,1°, и долготы восходящего узла у Земли, равном 11,7°. Разница между минимальным и максимальным (~8122 м/с) значениями оценки затрат характеристической скорости составила ~210 м/c, или ~3 % от максимального значения.

ЛИТЕРАТУРА

[1] Захаров Ю.А. Проектирование межорбитальных космических аппаратов. Москва, Машиностроение, 1984, 176 с.

[2] Гришин С.Д., Захаров Ю.А., Оделевский В.К. Проектирование аппаратов с двигателями малой тяги. Москва, Машиностроение, 1990, 224 с.

[3] Pierson B.L., Kluever C.A. Three-Stage approach to optimal low-thrust Earth-Moon trajectories. Journal of Guidance, Control and Dynamics, 1994, vol. 17, no. 6, pp. 1275-1282.

[4] Kluever C.A., Pierson B.L. Optimal low-thrust three-dimensional Earth-Moon trajectories. Journal of Guidance, Control and Dynamics, 1995, vol. 18, no. 4, pp. 830-837.

[5] Kluever C.A., Pierson B.L. Optimal Earth-Moon trajectories using nuclear electric propulsion. Journal of Guidance, Control and Dynamics, 1997, vol. 20, no. 2, pp. 239-245.

[6] Casaregola C., Geurts K., Pergola P., Biagioni L., Andrenucci M. Mission analysis and architecture definition for a small electric propulsion transfer module to the Moon. In: Proceedings of the 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. AIAA 2007-5232, Cincinnati, USA, July 8-11, 2007.

[7] Ивашкин В.В., Петухов В.Г. Траектории перелета с малой тягой между орбитами спутников Земли и Луны при использовании орбиты захвата Луной. ИПМим. М.В. Келдыша РАН. Препринт № 81. Москва, 2008, 32 с.

[8] Петухов В.Г. Робастное квазиоптимальное управление с обратной связью для перелета с малой тягой между некомпланарными эллиптической и круговой орбитами. Вестник МАИ, 2010, т. 17, № 3, с. 50-58.

[9] Grebow D.J., Ozimek M.T., Howell K.C. Advanced modeling of optimal low-thrust Lunar pole-sitter trajectories. In: Proceedings of the 60th International As-tronautical Congress. IAC-09-C1.5.4, Daejeon, Korea, October 12-16, 2009.

[10] Perez-Palau D., Epenoy R. Indirect optimization of low-thrust Earth-Moon transfers in the Sun-Earth-Moon system. In: Proceedings of the 68th International Astronautical Congress. IAC-17-C1.6.2, Adelaide, Australia, September 25-29, 2017.

[11] Иванюхин А.В., Петухов В.Г., Юн Сон Ук. Траектории перелета к Луне с минимальной тягой. Космические исследования, 2022, т. 60, № 6, с. 517-527. DOI: 10.31857/S002342062205003X

[12] Кувшинова Е.Ю. Методика определения оптимальной траектории перелета с малой тягой между околоземной и окололунной орбитами. Труды МАИ, 2013, вып. 68. URL: http://trudy.mai.ru/published.php?ID=41742

[13] Powers W.F., Tapley B.D. Canonical transformation applications to optimal trajectory analysis. American Institute of Aeronautics and Astronautics Journal, 1969, vol. 7, no. 3, pp. 394-399. https://doi.org/10.2514/3.5119

[14] Питьева Е.В. и др. Эфемериды EPM2008. Институт прикладной астрономии РАН. URL: ftp://ftp.iaaras.ru/pub/epm/EPM2008/ (дата обращения: 01.04.2024).

[15] Понтрягин Л.С., Болтянский В.Г., Гамкрелидзе Р.В., Мищенко Е.Ф. Математическая теория оптимальных процессов. 4-е изд. Москва, Наука, 1983, 392 с.

[16] Petukhov V.G. Minimum-thrust problem and its application to trajectory optimization with thrust switchings. In: Proceedings of the 64th International Astro-nautical Congress. IAC-13-C1.6.2, Beijing, China, September 23-27, 2013.

[17] Fehlberg E. Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta formulas with stepsize control. NASA TR R-287, 1968.

[18] Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. Москва, Мир, 1990, 512 с.

[19] Butcher J.C. Numerical Methods for Ordinary Differential Equations. John Wiley&Sons, Ltd., The Atrium, 2003.

[20] Шаманский Е.В. Методы численного решения краевых задач на ЭЦВМ. Часть 2. Нелинейные краевые задачи и задачи на собственные значения для дифференциальных уравнений. Киев, Наукова думка, 1966, 244 с.

[21] Аоки М. Введение в методы оптимизации. Основы и приложения нелинейного программирования. Москва, Наука, 1977, 344 с.

[22] Ахметшин Р.З. Плоская задача оптимального перелета космического аппарата с малой тягой с высокоэллиптической орбиты на геостационар. Космические исследования, 2004, т. 42, № 3, c. 248-259.

[23] Кувшинова Е.Ю. Применение многоразовых буксиров с ядерной электроракетной двигательной установкой в лунной программе. Космонавтика и ракетостроение, 2017, № 6 (99), c. 75-80.

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

Ссылку на эту статью просим оформлять следующим образом: Кувшинова Е.Ю. Исследование оптимальных пространственных траекторий перелета с малой тягой в системе Земля — Луна. Инженерный журнал: наука и инновации, 2024, вып. 7. EDN KOILKC

Кувшинова Екатерина Юрьевна — канд. техн. наук, ведущий научный сотрудник АО ГНЦ «Центр Келдыша». Область научных интересов: двигательные и энергетические установки космических средств межорбитальной транспортировки, динамика полета космических аппаратов с малой тягой. e-mail: [email protected]

Studying optimal spatial trajectories in a low thrust flight within the Earth-Moon system

© E.Yu. Kuvshinova

State Scientific Center of the Russian Federation "Keldysh Research Center", Moscow, 125438, Russian Federation

The paper considers spatial flight of the electric propulsion system powered spacecraft between the low lunar and near-Earth orbits with a minimum time. The paper objective is to assess the influence of characteristic velocity required for a flight between the lunar and near-Earth orbits, such as the trajectory ballistic parameters, flight angular range (true longitude) and the orbital plane position in vicinity of the Earth and the Moon, on the expenses. Thrust vector direction of the electric rocket propulsion system is determined using the Pontryagin maximum principle. Flight trajectories in the Earth—Moon system were optimized in terms of determining the thrust vector direction of the electric rocket propulsion system within the framework of a single formulation (end-to-end trajectory calculation) taking into account the Earth and the Moon gravity (without using the sphere of action method). Their position was determined according to the EPM 2008 ephemeris model.

Keywords: spacecraft, optimal flight trajectories, electric rocket propulsion system, low thrust, flight in the Earth-Moon system

REFERENCES

[1] Zakharov Yu.A. Proektirovanie mezhorbitalnykh kosmicheskikh apparatov [Design of inter-orbital spacecraft]. Moscow, Mashinostroenie, 1984, 176 p.

[2] Grishin S.D., Zakharov Yu.A., Odelevsky V.K. Proektirovanie apparatov s dvigatelyami maloy tyagi [Design of spacecraft equipped with low-trust engines]. Moscow, Mashinostroenie, 1990, 224 p.

[3] Pierson B.L., Kluever C.A. Three-stage approach to optimal low-thrust Earth— Moon trajectories. Journal of Guidance, Control and Dynamics, 1994, vol. 17, no. 6, pp. 1275-1282.

[4] Kluever C.A., Pierson B.L. Optimal low-thrust three-dimensional Earth—Moon trajectories. Journal of Guidance, Control and Dynamics, 1995, vol. 18, no. 4, pp. 830-837.

[5] Kluever C.A., Pierson B.L. Optimal Earth—Moon trajectories using nuclear electric propulsion. Journal of Guidance, Control and Dynamics, 1997, vol. 20, no. 2, pp. 239-245.

[6] Casaregola C., Geurts K., Pergola P., Biagioni L., Andrenucci M. Mission analysis and architecture definition for a small electric propulsion transfer module to the Moon. In: Proceedings of the 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. AIAA 2007-5232, Cincinnati, USA, July 8-11, 2007.

[7] Ivashkin V.V., Petukhov V.G. Traektorii pereleta s maloy tyagoy mezhdu orbitami sputnikov Zemli i Luny pri ispolzovanii orbity zakhvata Lunoy [Trajectories for flight between Earth and Moon satellite orbits using intermediate orbit with capture by Moon]. IPM im. M.V. Keldysha RAN. Preprint no. 81 [Keldysh Institute of Applied Mathematics, Russian Academy of Sciences. Preprint no. 81]. Moscow, 2008, 32 p.

[8] Petukhov V.G. Robastnoe kvazioptimalnoe upravlenie s obratnoy svyazyu dlya pereleta s maloy tyagoy mezhdu nekomplanarnymi ellipticheskoy i krugovoy

orbitami [Robust suboptimal feedback control for low-thrust transfer between noncoplanar elliptical and circular orbits]. Vestnik MAI — Aerospace MAI Journal,, 2010, vol. 17, no. 3, pp. 50-58.

[9] Grebow D.J., Ozimek M.T., Howell K.C. Advanced modeling of optimal low-thrust Lunar pole-sitter trajectories. In: Proceedings of the 60th International As-tronautical Congress. IAC-09-C1.5.4, Daejeon, Korea, October 12-16, 2009.

[10] Perez-Palau D., Epenoy R. Indirect optimization of low-thrust Earth-Moon transfers in the Sun-Earth-Moon System. In: Proceedings of the 68th International Astronautical Congress. IAC-17-C1.6.2, Adelaide, Australia, September 25-29, 2017.

[11] Ivanyukhin A.V., Petukhov V.G., Yoon S.W. Traektorii pereleta k Lune s mini-malnoy tyagoy [Minimum-thrust flight trajectories to the Moon]. Kosmicheskie issledovaniya — Cosmic Research, 2022, vol. 60, no. 6, pp. 481-490. https://10.1134/S0010952522050033X

[12] Kuvshinova E.Yu. Metodika opredeleniya optimalnoy traektorii pereleta s maloy tyagoy mezhdu okolozemnoy i okololunnoy orbitami [The methodology for determination of optimal flight path with low thrust between the low-Earth and low-Lunar orbits]. Trudy MAI, 2013, iss. 68. Available at: https://www.trudymai.ru/published.php?ID=41742

[13] Powers W.F., Tapley B.D. Canonical transformation applications to optimal trajectory analysis. American Institute of Aeronautics and Astronautics Journal, 1969, vol. 7, no. 3, pp. 394-399. https://doi.org/10.2514Z3.5119

[14] Pityeva E.V. et al. Efemeridy EPM2008 [Ephemeris EPM2008]. Institut pri-kladnoy astronomii RAN [Institute of Applied Astronomy, Russian Academy of Sciences]. https://ftp.iaaras.ru/pub/epm/EPM2008/

[15] Pontryagin L.S., Boltyansky V.G., Gamkrelidze R.V., Mishchenko E.F. Ma-tematicheskaya teoriya optimalnykh protsessov [Mathematical theory of the optimal processes]. 4th edition. Moscow, Nauka Publ., 1983, 392 p.

[16] Petukhov V.G. Minimum-thrust problem and its application to trajectory optimization with thrust switchings. In: Proceedings of the 64th International Astro-nautical Congress. IAC-13-C1.6.2, Beijing, China, September 23-27, 2013.

[17] Fehlberg E. Classical fifth-, sixth-, seventh-, and eighth-order Runge—Kutta formulas with stepsize control. NASA TR R-287, 1968.

[18] Hairer E., Norsett S.P., Wanner G. Solving Ordinary Differential Equations I. Nonstiff Problems. Berlin, Springer-Verlag, 1987 [In Russ.: Khayrer E., Nersett S., Vanner G. Reshenie obyknovennykh differentsialnykh uravneniy. Nezhestkie zadachi. Moscow, Mir Publ., 512 p.].

[19] Butcher J.C. Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, Ltd., The Atrium, 2003.

[20] Shamanskiy E.V. Metody chislennogo resheniya kraevykh zadach na ETsVM. Chast 2. Nelineynye zadachi i zadachi na sobstvennye znacheniya dlya differentsialnykh uravneniy [Methods for numerical solution of boundary value problems on a computer. Part 2. Nonlinear boundary value problems end eigenvalue problems for differential equations]. Kyiv, Naukova Dumka Publ., 1966, 244 p.

[21] Aoki M. Introduction to Optimization Techniques. Fundamentals and Applications of Nonlinear Programming. Macmillan, 1971 [In Russ.: Aoki M. Vvedenie v metody optimizatsii. Osnovy i prilozheniya nelineynogo programmirovaniya. Moscow, Nauka Publ., 1977, 344 p.].

[22] Akhmetshin R.Z. Ploskaya zadacha optimalnogo pereleta kosmicheskogo appa-rata s maloy tyagoy s vysokoellipticheskoy orbity na geostatsionar [Plane trajectories of optimal transfer of spacecraft with low thrust from elliptical orbit to geostationary one]. Kosmicheskie issledovaniya — Cosmic Research, 2004, vol. 42, no. 3, pp. 238-249.

[23] Kuvshinova E.Yu. Primenenie mnogorazovykh buksirov s yadernoy elektro-raketnoy dvigatelnoy ustanovkoy v lunnoy programme [Application of reusable tugs with nuclear electric propulsion system in the lunar program]. Kosmonavti-ka i raketostroenie — Cosmonautics and Rocket Engineering, 2017, no. 6 (99), pp. 75-80.

Kuvshinova E.Yu., Cand. Sc. (Eng.), Leading Researcher, State Scientific Center of the Russian Federation "Keldysh Research Center". Research interests: propulsion and power plants of the interorbital transportation space vehicles, dynamics of the low-thrust spacecraft flight. e-mail: [email protected]