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

УДК 517.9:521+523.3:629.7

Ю.А. Окишев, Ю.В. Клинаев

РАЗРАБОТКА МАТЕМАТИЧЕСКОЙ МОДЕЛИ БАЛЛИСТИЧЕСКОГО АНАЛИЗА ПЕРЕЛЕТА КОСМИЧЕСКОГО АППАРАТА С НИЗКОЙ ОКОЛОЗЕМНОЙ ОРБИТЫ В ТОЧКУ ЛИБРАЦИИ Ь1 СИСТЕМЫ «ЗЕМЛЯ-ЛУНА» КАК РЕШЕНИЕ ЧАСТНОЙ ОГРАНИЧЕННОЙ ЗАДАЧИ ТРЕХ ТЕЛ С УЧЕТОМ ПРЕЦЕССИИ ОРБИТЫ ЛУНЫ

Разработана численно-аналитическая модель баллистического анализа доставки космического аппарата с химическим ракетным двигателем с низкой околоземной орбиты в точку либрации Ь1 системы «Земля-Луна». Найдены опти-

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

Математическое моделирование, баллистический анализ, точка либрации, Ь1, система «Земля-Луна»

Yu.A. Okishev, Yu.V. Klinaev

DEVELOPMENT OF A MATHEMATICAL MODEL FOR THE BALLISTIC FLIGHT ANALYSIS

OF A SPACECRAFT FROM LOW-EARTH ORBIT TO THE LIBRATION POINT L1 OF THE «EARTH-MOON» AS THE DECISION OF A PRIVATE LIMITED THREE-BODY PROBLEM WITH CONSIDERATION FOR THE MOON ORBIT

A numerical-analytical model has been developed for the ballistic analysis made for the delivery of a spacecraft having a chemical rocket engine from low-Earth orbit in the L1 libration point of the «Earth-Moon». The optimal launch date, flight time, and the needed speed pulses have been determined. To perform the mathematical modeling CAD system and computer algebra MathCAD have been used.

Mathematic modeling, ballistic analysis, libration point, L1, the system «Earth-Moon»

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

Работа посвящена разработке математической модели и алгоритма поиска оптимальной баллистической траектории космического аппарата (КА) с химическим ракетным двигателем (ХРД) при перелете с низкой круговой околоземной орбиты (ее высоту принимаем равной 300 км) в точку либрации L1 системы «Земля - Луна».

В качестве критерия оптимизации предлагается рассматривать суммарный импульс. Задача проектно-баллистического анализа сводится к поиску минимального суммарного импульса скорости AVZ и, как следствие, по формуле Циолковского (1), минимально потребного топлива для перелета. Такой подход является общепризнанным [1] для проведения баллистического анализа.

где т т - масса топлива КА, т0 - начальная масса КА, 1у — удельный импульс тяги, который задается

двигательной установкой КА.

Будем рассматривать схему перелета с двумя включениями ХРД между некомпланарными орбитами, где первый импульс скорости АУ1 реализует переход КА на перелетный эллипс, лежащий в плоскости орбиты базовой орбиты, а второй импульс — скорости АV2 (2) в апогее перелетного эллипса реализует поворот плоскости орбиты на требуемый угол Аi и переход на орбиту точки либрации Ь1 (рис. 1). Для предварительного анализа используется методическая идея импульсной аппроксимации активных участков полета. Исходя из физических свойств коллинеарной точки либрации Ь1 очевидно, что точка Ь1 принадлежит радиус-вектору и плоскости орбиты Луны:

где VL1 — вектор скорости точки либрации Ь1, VКА — вектор скорости КА в точке либрации Ь1, Аг — разница между наклонениями орбит.

(1)

(2)

По теореме косинусов получаем:

(3)

Рис. 1. Схема двухимпульсного перелета

Для проведения анализа полагаем, что начало системы координат расположено в центре Земли, плоскость х-у совпадает с плоскостью эклиптики, ось х направлена в точку весеннего равноденствия. Ось ъ направлена в северный полюс Мира, ось у дополняет до правой тройки. Можно считать, что угол между плоскостью эклиптики и плоскостью земного экватора постоянен и равен 23,4354°.

Наклонение низкой околоземной (базовой) орбиты примем за 51,6° (старт с космодрома «Байконур»).

Из-за прецессии орбиты Луны ее наклонение к плоскости экватора Земли меняется с периодичностью в 18,6 лет. Исходя из энергетики перелета (3), выберем эпоху, когда наклонение орбиты Луны максимально, чтобы значение разности наклонений базовой и Лунной орбиты было минимальным.

За дату рассматриваемой эпохи выбираем 1 января, 12 часов дня каждого рассматриваемого года. Найдем наклонение орбиты Луны по формуле (4):

' , (4)

cos i

Луны

а

а = ГЛуны X V луны

а

■X ■Vy - у ■Vx,

(5)

(6)

где а - вектор интеграла площадей орбиты Луны, аг - проекция вектора - интеграла площадей орбиты Луны на ось z, x и у - проекции радиус-вектора Луны на оси x и y соответственно, Vx и Vy -

проекции вектора скорости Луны на оси x и y соответственно.

Радиус-вектор Луны и значение скорости, а также их проекции по осям возьмем из планетария DE-405, который разработан в JPL (Jet Propulsion Laboratory). Выберем эпоху, когда наклонение орбиты Луны максимально. Как видно из табл. 1, наклонение орбиты Луны максимально и составляет 28,443° в 2025 году, при этом будем считать, что в выбранную эпоху наклонение орбиты не изменяется.

Таблица 1

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

Год i о Луны Год i о Луны Год i о Луны

2012 22,513 2018 20,075 2024 28,195

2013 20,881 2019 21,568 2025 28,443

2014 19,526 2020 23,253 2026 28,258

2015 18,633 2021 24,894 2027 27,638

2016 18,396 2022 26,327 2028 26,584

2017 18,959 2023 27,458 2029 25,174

Так как в данной задаче влияние гравитации Луны на КА существенно, то использовать решения ограниченной задачи двух тел (7) некорректно. Будем использовать решение ограниченной задачи трех тел (8), в которой движение КА массы т, рассматривается в системе двух гравитирующих масс М1 и М 2 (Земли и Луны, соответственно) [2].

Н2 г тМ

= /тМ^, (7)

йг г

где т - масса КА, г - радиус-вектор КА относительно центра Земли, г - время, / - гравитационная постоянная, Мз - масса Земли.

<Н2 Я „М. _ М

_2

A2 " R3 J R23

= /-f^ -R) + f-f(R2 + R) , (8)

Я1 = а----М-------, (9)

1 М 2 + М1

И2 = а------М--------, (10)

2 М 2 + М1

где М1 - масса Земли, М2 - масса Луны, Я - радиус-вектор КА относительно общего барицентра (БЦ), Я1 - расстояние от БЦ до центра Земли, Я2 - расстояние от общего БЦ до центра Луны, а -большая полуось орбиты Луны.

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

Н-Я- = ш2 Я + / М2 2 - / М1 , (11)

йг2 (Я2 - Я)2 (Я1 + Я)2

п

Ш = 2- — , (12)

Т

T = 2-к-

a a

(13)

/■ (М1 + М 2)

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

Приравнивая правую часть уравнения (11) в произвольный момент времени к нулю, найдём расстояние от барицентра системы «Земля-Луна» до точки либрации Я . По уравнению (14) найдем расстояние от центра Земли до точки либрации Ь1 системы «Земля - Луна» ЯЬ1.

яи = Я + Я1, (14)

Для дальнейших расчетов введем коэффициент Х11(15), который характеризует относительное расстояние от центра Земли до точки либрации Ь1. В дальнейшем, умножая на этот коэффициент расстояние от Земли до Луны в любой момент времени можем определить расстояние до точки либрации Ь1. При этом скорость точки либрации УЬ1 (16) найдем из подобного треугольника (рис. 2)

= Я

‘-и _

VL1 _ VMoon ' „ , (16)

R1 + R2

R

' L1 y Moon D D R1 + R2

где RL1 - расстояние от центра Земли до L1, VMoon - скорость Луны относительно Земли.

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

Существует 2 типа решения поставленной задачи: при перелете из восходящего узла орбиты (рис. 3) и из нисходящего узла орбиты (рис. 4).

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

Рис. 2. Упрощенная схема перелета

Рис. 3. Траектория орбиты при перелете из восходящего узла орбиты, где И - точка либрации, 0 - положение КА в момент старта

Рис. 4. Траектория орбиты при перелете из нисходящего узла орбиты, где И - точка либрации, 0 - положение КА в момент старта

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

Решим ограниченную задачу трех тел, применив адаптивный метод Рунге-Кутты 6-го порядка (Кутты-Мерсона) с переменным шагом интегрирования для выбранных долгот восходящего узла и аргумента широты. По проекциям траектории КА на плоскости ХУ и XZ обнаруживаем, что в точку либрации и1 системы «Земля-Луна» не попадаем.

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

65

краевую задачу в среде МаШСАБ [2], найдем долготу восходящего узла, аргумент широты и первый импульс скорости для выбранной даты подлета в точку либрации. При проекции на плоскости ХУ и XZ «промаха» не обнаруживаем. Зная начальные параметры орбиты и первый импульс скорости, мы можем определить второй необходимый импульс скорости.

Решая оптимизационную задачу в среде МаШСАБ [2], найдем такое оптимальное время перелета для выбранной даты попадания в точку либрации, когда импульс скорости минимальный (рис. 5).

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

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

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

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

Таким образом, в рамках рассмотрения апреля в качестве выбранной эпохи, можем сделать вывод, что оптимальной датой попадания в точку либрации Ь1 системы «Земля-Луна» является 103-й день (рис. 6), что соответствует 12 часам дня 13 апреля 2025 года. При этом оптимальное время перелета составляет 4,25 суток (табл. 2) при старте из окрестности восходящего узла базовой орбиты.

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

Таблица 2

Значения второго импульса скорости, оптимального времени перелета, долготы восходящего узла, аргумента широты и суммарного импульса скорости в зависимости от даты подлета к точке либрации L1 системы

«Земля-Луна» для старта из нисходящего узла орбиты

Дата подлета с начала эпохи, дней ау2 , м/с Тарг, дней ^ ,° и ,° /с м/ <

90 713.937 3.523 20.819 -156.875 3798

92 728.963 3.523 42.923 -144.418 3816

94 754.329 3.643 72.787 -141.974 3845

96 790.256 3.81 107.042 -148.97 3883

98 816.785 3.973 140.376 -161.394 3912

100 824.073 4.11 171.446 -175.933 3920

102 810.873 4.21 201.494 169.401-360 3908

104 779.330 4.273 232.07 156.135-360 3876

106 737.161 4.293 263.596 146.128-360 3834

108 700.458 4.257 293.872 141.938-360 3796

110 688.550 4.15 -40.905+360 145.616-360 3782

112 721.435 3.947 -21.519+360 156.809-360 3811

114 675.648 3.913 -5.85+360 170.982-360 3758

116 699.287 3.697 8.352+360 -170.279 3782

118 729.804 3.493 26.369+360 -152.726 3813

120 749.546 3.483 51.451+360 -142.708 3836

Таблица 3

Значения второго импульса скорости, оптимального времени перелета, долготы восходящего узла, аргумента широты и суммарного импульса скорости в зависимости от даты подлета к точке либрации L1 системы

«Земля-Луна» для старта из восходящего узла орбиты

Дата подлета с начала эпохи, дней АУ2, м/с Тар, , дней ^ ,° и ,° АУЕ, м/с

90 873.741 3.527 -129.531 -23.199 3958

92 802.393 3.523 -89.992 -34.544 3889

94 726.620 3.643 -56.42 -36.751 3817

96 679.608 3.81 -32.793 -29.919 3773

98 655.950 3.973 -16.509 -17.831 3751

100 644.094 4,11 -3.72 -3.705 3740

101 640.510 4.167 2.231 3.5 30737

102 638.228 4.213 8.25 10.563 3735

102.5 637.669 4.233 11.366 13.985 3735

103 637.607 4.25 14.6 17.305 3735

103.5 638.159 4.263 17.988 20.498 3735

104 639.467 4.277 21.57 23.523 3737

106 655.593 4.297 38.756 33.316 3753

108 697.085 4,26 62.599 37.308 3793

110 769.007 4.153 -266.269+360 33.423 3863

112 870.575 3.957 -230.792+360 22.170 3961

114 894.654 3.803 -194.74+360 5.29 3971

116 868.553 3.693 -158.396+360 -11.709 3977

118 791.648 3.497 -118.296+360 -27.129 3952

120 732.629 3.483 -78.794+360 -36.203 3903

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

1. Оптимальной датой попадания в точку либрации Ь1 является 13 апреля 2025 года, при этом время перелета составляет 4,25 дня.

2. Суммарный импульс скорости равен 3735 м/с, где первый импульс скорости составляет 3097 м/с, а второй импульс скорости - 637,607 м/с.

3. Оптимальная траектория реализуется при перелете из восходящего узла орбиты.

Таким образом, разработана математическая модель и алгоритм поиска оптимальной баллистической траектории КА с ХРД при перелете с низкой круговой околоземной орбиты в точку либрации Ь1 системы «Земля - Луна» для любой эпохи.

ЛИТЕРАТУРА

1. Константинов М.С. Механика космического полета: учеб. для втузов / М.С. Константинов, Е.Ф. Каменков, Б.П. Перелыгин, В.К. Безвербый, В.П. Мишин; под ред. В.П. Мишина. М.: Машиностроение, 1989. 407 с.

2. Себехей В. Теория орбит: ограниченная задача трех тел; пер. с англ. / В. Себехей; под ред. Г.Н. Дубошина. М.: Наука. Гл. ред. физ.-мат. лит. 1982. 656 с.

3. Ракитин В.И. Руководство по методам вычислений и приложений МАТИСАБ: учеб. пособие / В.И. Ракитин. М.: Физматлит, 2005. 264 с.

Окишев Юрий Александрович -

аспирант кафедры «Техническая физика и информационные технологии»

Энгельсского технологического института (филиала) Саратовского государственного технического университета имени Гагарина Ю.А.

Клинаев Юрий Васильевич -

доктор физико-математических наук, профессор кафедры «Техническая физика и информационные технологии»

Энгельсского технологического института (филиала) Саратовского государственного технического университета имени Гагарина Ю.А.

Статья поступила в редакцию 17.09.12, принята к опубликованию

Yuri A. Okishev -

Postgraduate

of the Department Technical Physics

and Information Technology

Engels Institute of Technology

Part of Gagarin Saratov State Technical University

Yuri V. Klinaev -

Dr. Sc., Professor

of the Department Technical Physics

and Information Technology

Engels Institute of Technology

Part of Gagarin Saratov State Technical University