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

Электронный журнал «Труды МАИ». Выпуск № 68

www.mai.ru/science/trudy/

УДК 629.78

Методика определения оптимальной траектории перелета с малой тягой между околоземной и окололунной орбитами

Кувшинова Е. Ю.

Исследовательский центр имени М. В. Келдыша, Онежская ул., д. 8, Москва, 125438,

Россия E-mail: kerc@elnet.msk.ru

Аннотация

В работе рассматривается задача перелёта космического аппарата с электроракетной двигательной установкой между околоземной и окололунной орбитами за минимальное время.

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

Одной из особенностей представленной методики является расчёт оптимальной траектории перелёта без разбиения на участки, соответствующие сферам действия, и связанных с этим методическими погрешностями. Другие особенности заключаются в использовании точных эфемерид Луны EPM 2008 и в постановке краевой задачи с учётом неединственности ее решения.

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

Введение

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

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

При определении оптимальных законов управления вектором тяги необходимо учитывать гравитационные воздействия двух притягивающих центров в рамках задачи трех тел. Традиционно поиск траекторий перелета в системе Земля-Луна проводится с разбиением траектории на отдельные участки, обычно на сферы действия. Такой подход вносит методическую погрешность, связанную с разбиением траектории на участки и с дальнейшей их стыковкой. Для построения оптимальной траектории точку стыковки участков траектории необходимо оптимизировать, что потребует гораздо большего объёма расчётов, чем расчёт траектории напрямую без разбиения траектории на участки.

Оценки значения характеристической скорости на перелёт с малой тягой между околоземной и окололунной орбитами у разных авторов [1-5] различны и находятся в широком диапазоне от 7,5 до 9,8 км/с. Различия в характеристической скорости могут быть обусловлены использованием приближенных методов расчета траекторий перелета, наличием большого числа экстремумов рассматриваемой задачи, зависимостью характеристической скорости от исходных данных, в частности, даты старта, параметров начальной и целевой орбит.

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

Траектория перелета определялась с одновременным учетом гравитационных полей Земли и Луны (ограниченная задача трех тел) и без разбиения на участки. Направление вектора тяги определялось по закону, полученному с использованием принципа максимума Понтрягина [6-7]. Тяга и удельный импульс тяги предполагались постоянными на всей траектории. Рассматривался непрерывный режим работы электроракетных двигателей. При расчётах действием Солнца и других небесных тел, нецентральностью гравитационных

полей Земли и Луны, наличием атмосферы Земли и солнечного давления пренебрегалось. Эфемериды Луны определялись согласно модели EPM 2008 [8], разработанной «Институтом прикладной астрономии РАН».

Математическая модель движения космического аппарата

При выборе равноденственных элементов для уравнений движения в рамках задачи 2-х тел реализуется хорошая сходимость краевых задач, а решённая задача дает хорошее начальное приближение сопряженных переменных при варьировании исходных данных (ускорения, удельного импульса). Эти свойства системы равноденственных элементов обусловили их использование в данной работе в рамках ограниченной задачи трех тел. Другими положительными свойствами от использования равноденственных элементов является то, что они позволяют:

- исключить вырождение уравнений движения при нулевых значениях эксцентриситета и наклонения;

- исследовать влияние угловой дальности перелета на набор характеристической скорости;

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

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

Принятые для расчёта равноденственные элементы F} [1, 9] связаны

следующими соотношениями с оскулирующими элементами:

Ь = £

ех = е • 008(О + ю)

еу = е • 8т(0 + ю) (1)

¿х = • с°8(Д)

1у = • вт^

Б=0+ю+Э

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

Использование эфемерид Луны, явно зависящих от времени, и расчёт траектории со сменой начала отсчёта ведут к нестандартной записи ускорений, действующих на КА. Ниже приводится анализ вида ускорений. Рассматривалась система векторов (рисунок 1), определяющая положение трёх тел (Земля-Луна-КА) в инерциальной системе координат ОиХиУи2и. На рисунке 1 вектора с началом в центре инерциальной системы координат Ои обозначены буквой р, вектора с началом в центре Земли обозначены буквой Я, вектора с началом в центре Луны обозначены буквой г, вектора с окончанием в центре Земли обозначены индексом З, вектора с окончанием в центре Луны обозначены индексом Л, вектора с окончанием в центре масс КА - без индекса.

КА

Земля

Луна

Уи

Хи

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

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

Р =

Я

з я Я3

^ , -возм. л 3 + а КА г3

_ .. . хз , -возм. рл = дз + ал ,

гз

— 11 . Я л , -возм. рз = Дл + аз :

Я3

(2)

(3)

где дз=3,986 1014 м3/с2 - гравитационный параметр Земли, дл=0,049 1014 м3/с2 -

гравитационный параметр Луны, a КА ', a л ', a3 ' - ускорения от планет (кроме Земли и

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

Из треугольников Ои - Земля - КА и Ои - Луна - КА следует:

Р = Рл + r = Рз + Ял + (5)

Р = Рз + R (6)

После двукратного дифференцирования соотношений (5) и (6) получено:

Р = Рз + R л + Г , (7)

Р = Рз + R (8)

С целью исключения векторов р и рз с началом в инерциальной системе координат Ои соотношения (2 - 4) были подставлены в (7 - 8): R r —возм. _ ,, R л . —возм.

- д з • ^у - д л -^+а КА = д л • -у+а возм + я л + г, (9)

Я3 г3 ял

II Я II г + возм. _.. Я л, -возм. пт

-дз '^-дл ^ + аКА = дл +а з + Я (10) Я г Ял

Искомые гравитационные ускорения г и Я имеют вид:

7-11 Я II ^и ^Т I оВозм. ^возм. /Пч

г = -дз •-у-дл •"у-дл ^-у--Ял + аКА -аз , (11) Я3 г3 Я3л

Я = -д з Л - д л •^ - д л •Ял+аКААм. - авозм. (12)

Я г Ял

Расстояние между КА и Землёй по астрономическим меркам пренебрежимо мало, следовательно, ускорения, действующие на КА и Землю авд^ и авозм., приблизительно равны и их разность равна нулю. С учётом этого соотношения (11 - 12) примут вид:

.. Я г Я -и- Я + г г Я ••

г = -дз •^Г -дл ^-у -дл ^ - Ял =-дз--л-У -дл ^-у -дл - Ял, (13)

Я3 г3 Я л (Я л + г)3 г3 Я л

^ Я г Я л Я Я - Я л Я л

Я = -дз ^-дл •-у-дл —л = -дз -дл---дл —л (14)

Я3 г3 Я3 Я3 (Я - Ял)3 Я3

При составлении уравнений движения в оскулирующих элементах составляющие

„ г R „

ускорений -д л • — и - дз ——, соответствующие ускорению в ограниченной задаче двух Г3 Я3

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

а - и Ял + г дл •Ял Я (15)

ал --дз---о--Я л, (15)

(Я л + г)3 ял

аз -Ил-"ф (16)

л (Я-Ял)3 я3л

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

вторые производные Луны Я л .

Уравнения для фазовых переменных в безразмерных равноденственных элементах [1, 9] имеют вид:

ёЬ 2 т - Ь •Т•к

- Ь • (Т • ((1 + к) • соб Б + к ех ) + 8 • Б1п Б + W • еу -л- к)

ёт

ёех

ёт ёеу

—— - Ь • (Т • ((1 + к) • б1п Б + к • еу ) - 8 • соб Б - W • ех •л к) ёт у

^--!• Ь • W •к^ собБ (17)

ёт 2

¿1у 1

—— ----Ь • W • к • ф^ б1п Б

ёт 2 у

ёБ 1

ёт Ь3 • 2 ёт _ 2 • N •л ё1

* -1

ёт

- •Ь • W • л^к

к2

ёт тэрду2

т

где h, ex, ey, ix, iy, F, m - фазовые переменные, N - мощность, подаваемая на электроракетную двигательную установку, лт - тяговый КПД ЭРДУ, Ьэрду - удельный импульс тяги ЭРДУ,

к =_1_, Л = ix ' sinF-iy • cosF, ф = 1 + ix2 + iy2, T, S, W - проекции вектора

1 + ex • cosF + ey • sinF

ускорения a от тяги электроракетной двигательной установки и возмущающих ускорений на оси орбитальной системы координат (трансверсальное, радиальное и нормальное направления). Выражение для ускорения a с учетом (15, 16) принимает следующий вид:

- 2• N'Лт ^ л -a =-— • R0 + A • a¿,

1эрду • m

где aj = {aз ,aл} - возмущающие гравитационные ускорения, действующие на КА, в

геоцентрической системе координат или в селеноцентрической системе координат,

'cos у cos аЛ Ro = sin у cos а vsin а

а и у - углы, рассматриваемые в качестве параметров управления,

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

А

2 2 2 - sinF + - • iy • (cosF • ix + sinF • iy) cosF - —• ix • (cosF • ix + sinF • iy) — • (cosF • ix + sinF • iy)

А =

cosF + y iy • (sinF • ix - cosF • iy) sinF - —• ix • (sinF • ix - cosF • iy) — • (sinF • ix - cosF • iy)

ф y v x ф x v x ф

ix -iy

, • 2 • 2

2 . 2 . 1 -ix -iy

- — i — i iy л ix

ф 1 ф л ф

Безразмерные переменные связаны с соответствующими размерными переменными (размерные величины отмечаются звёздочкой) следующими соотношениями:

t r V m

t =-, r =-, V =-, m =-

t хар гхар "хар m хар

где 1 хар - гха^^|-хар— характерное время, гх&р=384 400 000 м - характерное расстояние

Гд

равное большой полуоси орбиты Луны, Ухар ---характерная скорость,

\ гхар

где д - гравитационным параметр основного притягивающего центра, тхар - характерная масса.

Согласно принципу максимума Понтрягина гамильтониан Н для задачи

максимального быстродействия имеет следующий вид: 8 2

H = • fi = Yh • h2 • T+ • h• (T • ((1 + к)• cos(F) + к-ex ) + S• sin(F) + W• ey -л-к) + i=1 x

+ Yey • h• (T• ((1 + к)• sin(F) + к• ey )-S• cos(F)-W• ex ^к) + -Yix • 1 • h • W • к • ф • cos(F)-Yiy • 1 • h • W • к • ф • sin(F) +

+ • (^r^T - h • W • л • к) - m •Ym +Yt, (18)

h3 •к2

где vh • Ve Ve Vi Vi VF Vm Vt - переменные, сопряжённые переменным h, ex,

' x' y y' ' '

ey, ix, iy, F, m, t

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

A

cosу^cosа = , х , (19)

Va х2 + Ar2 + An2

A

sin у^ cos а = , r , (20)

"\¡A х2 + Ar2 + An2

A

sin а= n , (21)

Л/А х2 + Ar2 + An2

где Ах = h •^Vh + h• (1 + к)• cos(F)+ hex>yex + h• (1 + к)• sin(F)+ hey>ye ,

Ar = h• sin(F)•Ve* -h• cos(F)•yey , An =-h^куf + h•Л•к•ey•V^ -h•л•к•ex•V^

-1 • h • к • ф • cos(F) • Vix -1 • h • к • ф • sin(F) • Viy •

Уравнения для сопряжённых переменных имеют вид:

^ = _( _ 3 • + т ^ + Ах ах ( к • ь4 аь 1

ат _ 5АГ

— + Б--L + Аг

аь аь г

аБ аАп л

— + W--— + А

аь аь

аш п а^

^ ех -( - 2 • V р ак + т • аА х + А х ат + Б • аАг + Аг аБ + аАп + Ап аш

ах к3 • ь3 аех аех аех аех аех аех аех

^ еу -( - 2 • V Е ак + т • аА х + А х ат + Б • аАг + Аг аБ + аАп + Ап аш

ах к3 • ь3 аеу аеу аеу аеу аеу аеу аеу

¿V\ аА

^ =-(Т- ^ + А

ах

ах

¿V

а1х

У =-(Т + Ах

а, 1

ат _ аАг

х -+ Б--L + Аг

аи аи

= -( - 2 •

У

V,

'X

ат

аГ

с аАг л

+ Б--L + Аг

у а1у

аБ аАп л

--+ Ш--— + А

а1х а1х

+ ш + Аг

а1у

аш

п • ^ )

а1х

аш

а1у

(22)

к3 • ь

-ак+т ^+а х 3 аБ аБ 1

ат _ аАг

— + Б--L + Аг

аБ аБ г

аБ аАп

— + Ш--- + Ап

аБ аБ г

аш

Ж

+А

т /л ат аБ

—т = -(А х--+ Аг--

ах ат ат

¿V х ,к ат аБ .

—— = -(А х--+ Аг--+ Ап

I ах х ах г ах г

ат

аш

)

ах

)

Ч

астные произв одные,

входящие в (22), не приводятся ввиду их громоздкости.

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

Координаты Луны явным образом зависят от времени. Для обеспечения автономности в систему (17) введено последнее уравнение. Так как в автономных системах гамильтониан постоянен, свойство автономности оказывается полезным при контроле соответствия системы уравнений для сопряжённых переменных (22) системе уравнений движения (17).

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

1А ах

16-е уравнение - — = 1.

ах

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

Использование принципа максимума Понтрягина приводит к двухточечной краевой задаче с краевыми условиями, приводимыми далее. Решение краевой задачи сводится к задаче математического программирования, которая решалась совместно методами Ньютона, градиентного спуска и координатного спуска [13 - 14].

)

)

)

Верификация программного обеспечения

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

При составлении программы расчёта оптимальных траекторий проверялось соответствие сопряженных уравнений уравнениям движения. Для этого использовалось свойство постоянства гамильтониана для автономных систем. В результате интегрирования уравнений движения для всех проведенных расчетов изменение гамильтониана не превышало 10-6 % вне зависимости от исходных данных и точности интегрирования. Изменение гамильтониана такого порядка свидетельствует о соответствии сопряжённых уравнений уравнениям движения.

При моделировании траектории последовательно в равноденственных селеноцентрических элементах, а затем в геоцентрических равноденственных элементах необходимо чтобы она оставалась неизменной независимо от выбранных элементов или момента перехода из одной системы в другую. Для подтверждения этого была составлена программа интегрирования уравнений движения параллельно в селеноцентрических равноденственных элементах и в геоцентрических равноденственных элементах в течение всей траектории. Параллельное интегрирование необходимо для того, чтобы переменные, определяющие траекторию и в селеноцентрических, и в геоцентрических равноденственных элементах определялась в одни и те же моменты времени. Для тестовых расчётов были рассмотрены траектории, полученные при параллельном интегрировании уравнений движения в геоцентрических и селеноцентрических равноденственных элементах, в течение 300-суточного пассивного движения по круговой ОИСЛ с начальной высотой 7 000 (10 000, 30 000) км.

При интегрировании уравнений движения существует параметр, обозначаемый То1 [15], которым ограничивается оценка ошибки интегрирования на шаге, определяемая методом интегрирования. Этот параметр в методе Рунге-Кутта-Фельберга определяет точность интегрирования системы дифференциальных уравнений. Увеличение параметра То1 приводит к снижению точности численного решения и к увеличению быстроты счёта. Соответственно уменьшение параметра То1 приводит к увеличению точности, к уменьшению быстроты счета (за счёт большего числа шагов) и ограничено точностью машинного представления числа. Значения расхождений декартовых координат, полученных при параллельном интегрировании траекторий, с учётом изменения параметра То1 приведены в

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

Таблица 1

Результаты параллельного интегрирования траекторий в геоцентрической и

селеноцентрической системах координат.

Расхождения по декартовым координатам

Параметр То1 Ах, км Ау, км А2, км Число шагов интегрирования

10-10 6,95 10-2 1,88 10-1 9.66 10-3 233 732

10-11 2,04 10-3 5,43 10-3 2,75 10-4 329 836

10-12 8,45 10-3 2,22 10-2 1.02 10-3 566 766

10-13 2,00 10-2 5,30 10-2 2.52 10-3 1 183 029

Начальные параметры орбиты, по которой происходит движение буксира: аргумент широты и:=0°, долгота восходящего узла й=0°, наклонение орбиты 1=0°, эксцентриситет е=0, радиусы перицентра и апоцентра г„=га=8 738 км

10 1,22 10-3 2,84 10-3 1,21 10-3 110 018

10-11 2,08 10-2 4,65 10-2 1,96 10-2 152 436

10-12 5,69 10-4 1,03 10-3 5,73 10-4 256 818

10-13 1,92 10-2 4,26 10-2 1,79 10-2 523 149

Начальные параметры орбиты: и:=0°, й=0°, 1=0°, е=0, гл=га=11 738 км

10 2,18 10-2 2,37 10-2 6,07 10-3 7 083

10-11 7,33 10-3 1,58 10-2 5,39 10-3 9 339

10-12 1,19 10-2 1,89 10-2 5,84 10-3 12 763

10-13 1,98 10-2 2,34 10-2 6,35 10-3 20 750

Начальные параметры орбиты: и:=0°, й=0°, 1=0°, е=0, гл=га=31 738 км

Как видно из таблицы 1, при уменьшении параметра То1 до определённого предела расхождение по координатам (характеризующее ошибку интегрирования) уменьшается, что является характерной особенностью при численном интегрировании дифференциальных уравнений. При моделировании пассивного движения КА по окололунным орбитам с высотами 7 000, 10 000 и 30 000 км в геоцентрической и селеноцентрической системах координат расхождения по декартовым координатам составили от 6 до 190 м. Приведённые выше расхождения по декартовым координатам показывают малость ошибок интегрирования. Это даёт основание полагать, что эти ошибки незначительным образом скажутся на величине определяемой характеристической скорости.

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

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

Т, сутки

Рисунок 2 - Зависимость радиуса апоцентра, перицентра и наклонения орбиты относительно Земли от времени перелёта

Расхождения значений декартовых координат в конце траектории, через 30 суток пассивного движения, составило от 100 до 500 м. Основную долю в эти величины вносят погрешности декартовой системы координат, о чём свидетельствует количество шагов интегрирования для декартовой системы координат равное 10 171, для равноденственных элементов - 2 619. Незначительность ошибок интегрирования указывает на малость их влияния на величину характеристической скорости перелёта.

Методика определения оптимальной траектории

Как показано в работе [16] при перелёте на геостационарную орбиту существуют экстремали, полностью удовлетворяющие условиям принципа максимума и различающиеся угловой дальностью (истинной долготой) перелёта. При перелёте между окололунной и околоземной орбитами возможно существование подобных экстремалей, но различающихся угловой дальностью перелёта, как у Земли, так и у Луны. В работе вводится аналог угловой дальности Б, учитывающий целое число реализуемых витков относительно Луны и число витков у Земли вместе с дробной частью:

Б = [Бл ] • 2п + Бз,

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

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

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

(

1

ех

еУ ■х

¡У Б

Л

у 1=1

0

х

еУ

¡0

¡0 Б0 ,

п

а

Г к ^

Гп

Г к

к

г=г

(23)

к

Если начальное число витков у Луны известно, то ставится задача определения траектории с заданным числом витков у Луны и реализующая "раскрутку" у Луны и "скрутку" у Земли. Граничные условия для такой задачи могут быть следующие:

( \

1

ех

еУ ■х

¡У Б

у г=г

11° ^

е0 ( Гп

еУ ■0 , Га

■

¡0 V Бл

Б0 ,

л

Г, к ^

г=г

к

п

Г к

Га к

, Б к ,

V у

(24)

0

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

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

При реализации такого подхода целесообразно для обеспечения сходимости краевой задачи принять начальную высоту околоземной орбиты ~250 000 км. Когда приближенное решение задачи перелёта на высокую ОИСЗ будет получено, снижать высоту околоземной орбиты шагами ~50 000 км до заданной высоты 800 км.

2. После получения приближенного решения краевую задачу целесообразно решать точно. Граничные условия имеют вид:

(

И

ех

еУ 1х

Ч

Б

Л

У 1=1

0

И0 ^ ( И

е0 ех

е0 еУ

10 , 1х

10 1У

г0, V ~

Л

Г Ьк ^

'X

1=\

1х •к

Vг У

(25)

к

где Ьк, еХ, ек, £,

к ~к

еу, Г - зафиксированы.

Возможна постановка краевой задачи перелёта при свободной переменной F(tk), тогда сопряжённая переменная ^г(^к) должна быть равна нулю (то есть выполняться условия трансверсальности).

В результате решении краевых задач, интегрируя уравнения в равноденственных элементах, дробная часть числа витков у Луны в момент перехода из

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

к

к

е

У

к

У

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

С помощью разработанного программно-методического обеспечения в рамках ограниченной задачи трёх тел решена краевая задача пространственного перелета с ОИСЛ высотой 100 км и наклонением плоскости орбиты перпендикулярной опорной плоскости на ОИСЗ высотой 800 км и наклонением 51,6°. За опорную плоскость принята плоскость земного экватора. Долгота восходящего узла стартовой окололунной орбиты принята равной нулю, целевой околоземной орбиты - 286°. Расчёт был проведён на дату старта с

«-» 3

окололунной орбиты - 10 января 2020 года. Начальное ускорение принято равным ~1,6 102 ^ ^ м/с , удельный импульс тяги электроракетных двигателей - 3000 с. Проекции траектории

межорбитального перелёта ОИСЛ-ОИСЗ приведены на рисунке 3.

, км

У.

-600

400 000 -л

Хз, к»

000 -400 000 -200 000 0 200 000 400 000 60

-400 000

-600

200 000 п

Тг4 ч

\

я § \ Хз, км

000 -400 000 -200 000 0 200 000 400 000 60

V

-200 000

а)

б)

У™

-50

14 000 П

> \Xл,

000 -40 000 -30 000 -20 000 -10 000 0 10 0 00 20 0

-6 000

-50

40 000

\

\ у,

000 -40 000 -30 000 -20 000 -10 000 0 И^/10 1 00 20 0

У

-30 000

в) г)

Рисунок 3 - Проекции траектории межорбитального перелёта ОИСЛ - ОИСЗ на плоскость

экватора Земли а, б) в геоцентрической, в, г) в селеноцентрической системах координат

км

Точность решения краевой задачи составляет несколько метров от заданного значения высоты круговой орбиты 8ОО км.

Значение характеристической скорости, необходимой на перелёт ОИСЛ-ОИСЗ, составило 8665,5 м/с при продолжительности перелёта 51 сутки. Вопросы, связанные с влиянием на значение характеристической скорости различных баллистических параметров (даты старта, параметров начальной и целевой орбит), требуют дальнейших исследований.

Заключение

1. Разработано программно-методическое обеспечение для расчёта оптимальных траекторий перелёта с малой тягой между околоземной и окололунной орбитами в рамках ограниченной задачи трёх тел без разбиения траектории на участки.

2. Проведёны верификационные расчеты, которые позволяют сделать вывод о работоспособности разработанного программно-методического обеспечения,

3 2

3. При дате старта 1О января 2О2О года, начальном ускорении 1,61О- м/с , фиксированном числе витков у Луны (34 витка) и у Земли(143 витка) получена оценка затрат характеристической скорости на перелёт ОИСЛ-ОИСЗ, которая составила 8665,5 м/с.

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

1. Гришин С. Д., Захаров Ю. А., Оделевский В. К. Проектирование космических аппаратов с двигателями малой тяги. М.: Машиностроение, 199О - 224 с.

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

3. Бэттин Р. Наведение в космосе. М.: Машиностроение, 1966 - 448 с.

4. 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. 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2007

5. Kluever C. A., Pierson B. L. Optimal Earth-Moon Trajectories Using Nuclear Electric Propulsion. Journal of Guidance, Control, and Dynamics Vol. 20, No. 2, March-April 1997

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

7. Брайсон А., Хо Ю-Ши Прикладная теория оптимального управления. М.: Мир, 1972 - 544 c.

8. Питьева Е. В. Эфемериды EPM2008. Институт прикладной астрономии РАН ftp://quasar.ipa.nw.ru/incoming/EPM2008/

9. Петухов В. Г. Квазиоптимальное управление с обратной связью для многовиткового перелета с малой тягой между некомпланарными эллиптической и круговой орбитами. Космические исследования, 2011, том 49, № 2 , с. 128-137

10. Fehlberg E. Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control. NASA TR R-287, 1968 - 88 p.

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

12. Butcher J. C. Numerical Methods for Ordinary Differential Equations. England, John Wiley&Sons, 2003 - 425 с.

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

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

15. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980 - 279 с.

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