CC BY
70УДК 531.01:629.78
БОТ 10.18698/2308-6033-2017-6-1666
Алгоритм нахождения траектории перелета между двумя эллиптическими орбитами
© Н.В. Островский
Вятский государственный университет, г. Киров, 610000, Россия
Впервые представлен алгоритм расчета перелетной эллиптической орбиты между двумя эллиптическими орбитами. Решены задачи построения эллипса по двум радиус-векторам, выходящим из одного фокуса, расчета импульсов скоростей разгона и торможения, расчета ускорения, необходимого для изменения плоскости орбиты. Актуальность проделанной работы связана с тем, что в известных публикациях подобные задачи или не были рассмотрены, или отсутствуют конкретные алгоритмы расчетов. Расчет длины перелетной орбиты выполнен методом численного интегрирования. Движение по эллиптической орбите рассмотрено как суперпозиция кругового и радиального движений. Показано, что с увеличением длины большой полуоси перелетного эллипса увеличивается его эксцентриситет, что приводит с росту радиальной скорости космического аппарата на перелетной орбите, вследствие чего происходит увеличение скоростей разгона и торможения.
Ключевые слова: перелетная орбита, эллиптическая орбита, траектория полета к Марсу, численное интегрирование
Введение. Ключевая задача при проектировании межпланетных космических полетов — построение перелетной орбиты, соединяющей исходную и конечную орбиты планет. В качестве основных типов перелетных орбит рассматриваются гиперболические, эллиптические, гомановские и спиральные [1].
Для вывода космических аппаратов (КА) на гиперболические орбиты необходимы большие затраты энергии: вывод на них КА сопряжен с использованием ускорения до 10^. В качестве альтернативного метода для выхода на спиральную орбиту может быть использовано ускорение порядка 0,001^. Однако в таком случае перелет занимает много времени.
Наилучшей совокупностью ускорения, необходимого для вывода КА на перелетную орбиту, и времени перелета обладают эллиптические орбиты.
Рассмотрим частный случай эллиптических орбит — гоманов-ские орбиты, когда эллипс является касательным к двум окружностям (рис. 1). Нахождение эллипса гомановской орбиты не вызывает проблем, поскольку длина большой оси будет равна сумме радиусов круговых орбит, а длину малой оси можно выбрать, исходя из требуемого времени перелета и приемлемой скорости разгона КА.
Рис. 1. Гомановская орбита для перелета с одной круговой орбиты на другую компланарную круговую орбиту [1]:
1 — начальная орбита; 2 — конечная орбита; 3 — гомановская орбита для перелета
Более сложная задача заключается в построении эллиптической орбиты для перелета между двумя эллиптическими орбитами. Алгоритм нахождения такой орбиты в общедоступных средствах информации отсутствует [1, 2].
Главная цель данной статьи — построить алгоритм вычисления параметров эллиптической орбиты в зависимости от расположения точек отлета и прилета на исходной и конечной орбитах. Найдя геометрические параметры орбиты, можно вычислить свойственные ей параметры движения КА, а на их основе — скорости разгона и торможения.
Методы решения задач. Построение перелетной эллиптической орбиты сводится к нахождению формы эллипса, одна точка (точка В) которого находится на 1-й эллиптической орбите, а другая (точка С) — на 2-й (рис. 2). Причем все три эллипса имеют общий фокус А.
Решение данной задачи было найдено с использованием известных свойств эллипса и треугольника (теорема косинусов).
Рис. 2. Перелетная орбита между двумя эллиптическими орбитами
Скорость движения тела по круговой орбите легко найти из присущего этому виду движения равенства сил, воздействующих на тело, — силы тяготения и центробежной силы:
= mv^ r2 r
где Д — гравитационный параметр центрального тела; m — масса движущегося тела; r — длина радиус-вектора движущегося тела относительно центрального; vc — скорость кругового движения. Отсюда следует, что скорость кругового движения
Vc . (2)
На эллиптической орбите равенство, описываемое уравнением (1), соблюдается только в двух точках орбиты, когда угол между радиус-вектором тела относительно фокуса и большой осью эллипса равен 90°. Этот радиус называется фокальным параметром эллипса р, его величину определяют из отношения
p = b2/a, (3)
где b — длина малой полуоси; a — длина большой полуоси.
Вычислить круговую скорость тела в любой точке орбиты можно, используя закон сохранения момента количества движения (углового момента) [3]:
K = mvcr = const, (4)
где K — угловой элемент.
Это свойство нашло свое выражение во втором законе Кеплера, согласно которому радиус-вектор планеты описывает площади, прямо пропорциональные промежуткам времени [4].
Длина радиус-вектора точки на эллиптической кривой относительно фокуса эллипса в зависимости от величины угла (обозначим его 0 ), отсчитываемого от перицентра эллипса, может быть вычислена по уравнению [5]:
r (5)
1 + e cos 0
где е — эксцентриситет эллипса.
Зная скорость тела в любой точке, можно найти время движения тела по дуге ВС, вычислив интеграл:
0Cr 2 (0)d0
tBC = K • (6)
0b K
2
Однако данный интеграл имеет решение, только если e > 1, в то время как для эллипса e < 1. Тогда tBC можно найти с помощью численного интегрирования на электронно-вычислительной машине (ЭВМ). В представленных ниже результатах расчетов интегрирование было выполнено с шагом 0,0001 от величины угла BAC.
Более полное описание движения тела по эллиптической орбите позволяет сформировать модель орбитального движения небесных тел, в которой эллиптическое движение рассматривается как суперпозиция радиального и кругового движений (см., например, [6]).
На небесное тело, движущееся по эллиптической орбите, действует радиальное ускорение aR, величина которого равна векторной сумме ускорения силы тяжести Qq и центробежного ускорения Qq :
aR = aQ - QC . (7)
Под действием радиального ускорения за интервал времени At тело приобретает радиальную скорость:
VR (t) = VR (0) + QR (t)At, (8)
а движение этого тела вдоль радиус-вектора ведет к изменению его длины
r (t) = r (0) - Vr (t)At . (9)
Величина радиальной скорости имеет положительное значение при движении тела к центру тяготения и отрицательное — при движении от центра.
Величина круговой скорости связана с длиной радиус-вектора через К:
ус (О = К/го(0г(0. (10)
Данная система уравнений не имеет аналитического решения, поэтому был использован метод численного интегрирования. Расчеты проведены с шагом & = 60 с. Эфемериды Земли и Марса были найдены по программе Р1аперЬ 4.2 [7].
Алгоритм нахождения параметров перелетного эллипса. Эллипс однозначно можно описать любыми двумя из перечисленных параметров [5]:
• длина большой полуоси а;
• длина малой полуоси Ь;
• половина фокусного расстояния с;
• эксцентриситет эллипса е.
Данные параметры связаны между собой следующими уравнениями:
а2 = Ь2 +с2; (11)
Итак, нам известны длины АВ и АС и угол ВАС. Через точки В и С можно провести целую группу эллипсов, характеризующихся определенной взаимосвязью параметров. Иными словами, можно, например, найти е как функцию от а (рис. 3).
с
Е
Рис. 3. Чертеж для нахождения параметров эллипса
Выделим 2-й фокус эллипса Б. Теперь, исходя из свойств эллипса и треугольника, можно составить систему уравнений:
BC = yjAB2 + AC2 - 2AB • AC • cosZBAC ; (13)
DB = 2a - AB ; (14)
DC = 2a - AC ; (15)
BCD BC2 + DC2 - BD2 , (16)
cos ZBCD =-; (16)
2BC • DC
.BC. BC2 + AC2 - AB2 ,
cos ZBCA =-; (17)
2BC•AC
ZA CD = ZBCD - ZBCA ; (18)
2c = AD = VAC2 + DC2 - 2AC • DC • cos ZACD ; (19)
AB2+DB2 -AD2 ,
cos ZBAD =-; (20)
2 AB • DB
ZEAB = %-ZBAD ; (21)
ZEAC = ZEAB ±ZBAC . (22)
Выбор знака в уравнении (22) зависит от направления отсчета угла BAC.
Расчет параметров перелетного эллипса. Для иллюстрации использования данного алгоритма возьмем, например, одну точку на орбите Земли с долготой 46,673° и радиусом АВ = 1,4818-1011 м, а вторую — на орбите Марса с долготой 255,115° и радиусом АС = 2,2274 • 1011 м. В результате угол ВАС составит 208,442°. Результаты расчета эксцентриситета в зависимости от заданной величины полуоси орбиты представлены в табл. 1.
Таблица 1
Результаты расчетов эллиптических перелетных орбит
a, м • 1011 e p, м • 1011 vC, м/с L, м • 1011 ДГ, сут
1,9 0,22026 1,8078 27 094 6,9960 326,23
2,0 0,27205 1,8520 26 769 7,6437 386,56
2,3 0,40330 1,9259 26 250 9,0509 542,19
2,6 0,49235 1,9697 25 957 10,1830 693,98
3,0 0,57482 2,0088 25 703 11,4490 902,95
3,5 0,64533 2,0415 25 496 12,9500 1179,20
4,0 0,69562 2,0645 25 354 14,2640 1473,10
4,7 0,74560 2,0872 25 216 15,9390 1913,90
5,7 0,79373 2,1090 25 085 18,0900 2599,90
По данным табл. 1, с увеличением длины большой полуоси эллипса перелетной орбиты возрастает его эксцентриситет, а значит, и длина дуги ВС и время перелета.
С 1964 г. к Марсу было направлено около 30 КА и 23 из них достигли Красной планеты, но публикации, детально описывающие траектории этих полетов, отсутствуют. В работе [8] приведены расчеты перелета Земля — Марс, выполненные в приближении круговых орбит этих планет. В этом случае перелетная орбита должна быть близка к гомановской. Для тех же исходной и конечной точек вычисленное время перелета составляет 306 сут.
Разгон и торможение КА. Для того чтобы перевести КА с текущей орбиты небесного тела, с которым он был связан силой тяготения, на перелетную орбиту, ему нужно придать скорость, соответствующую данной точке перелетной орбиты. Как было показано выше, движение тела на эллиптической орбите характеризуется двумя скоростями: круговой и радиальной. Для точки на орбите Земли с долготой 46,673° и радиусом 1,4818-1011 м круговая скорость ус е = 30 053 м/с, а радиальная — уке = 417 м/с. Суммарная величина скорости разгона вычислена из суммы квадратов радиальной и круговой скорости. Результаты расчетов для рассматриваемого случая, а также величины радиальной и круговой скорости для перелетной орбиты в зависимости от длины большой полуоси эллипса и соответствующие им скорости разгона приведены в табл. 2.
Таблица 2
Скорости разгона для различных эллиптических перелетных орбит
a, м • 1011 vс, м/с vR , м/с Л^, м/с ЛvR, м/с Л^, м/с
1,9 33 055 -294,8 3002 -712 3085
2,0 33 456 -2885,0 3403 -3302 4742
2,3 34 117 -7085,0 4064 -7502 8532
3,0 34 842 -11 609,0 4790 -12 026 12 397
Как видим, с увеличением эксцентриситета перелетной орбиты резко возрастает величина радиальной скорости, что приводит к значительному увеличению суммарной скорости разгона Дуе . В работе [8] в результате расчетов скорость разгона составила 3628 м/с, при этом радиальная составляющая отсутствует.
В табл. 3 приведены скорости торможения КА для перехода с перелетной орбиты на орбиту Марса в точке с долготой 255,115° и радиусом 2,2274-Ю11 м, для которой ус м = 24 577 и ук М = 2235 м/с.
В рассматриваемом случае для перехода с перелетной орбиты на орбиту Марса необходимо увеличить круговую скорость КА. И так же, как и в предыдущем случае, увеличение эксцентриситета ведет к значительному росту скорости торможения.
Скорость торможения, найденная в работе [8], составляет всего 718 м/с, что свидетельствует о том, что допущения, при которых орбита Марса рассматривается как круговая, ведут к существенному искажению результатов.
Таблица 3
Скорости торможения для различных эллиптических перелетных орбит
a, м-1011 vC, м/с vR, м/с AvC, м/с AvR, м/с Avz, м/с
1,9 22 007 3101 2570 -866 2712
2,0 22 222 5710 2355 -3475 4198
2,3 22 684 9976 1893 -7741 7969
3,0 23 198 14 555 1379 -12 320 12 397
Главный вывод, который необходимо сделать по данным табл. 2 и 3, состоит в том, что с увеличением эксцентриситета перелетной орбиты значительно увеличивается радиальная скорость как в точке отлета, так и в точке прилета, так что целесообразно использовать эллиптическую орбиту с минимально возможным эксцентриситетом.
Для перехода КА на планетоцентрическую орбиту нет необходимости нацеливать КА точно на планету и доводить его кинетические параметры до параметров планеты.
Изменение плоскости орбиты. Рассмотренный выше алгоритм был найден для случая, когда исходная и конечная орбиты КА лежат в одной плоскости. Для изменения плоскости орбиты необходимо придать КА дополнительное ускорение, перпендикулярное плоскости орбиты [1]. При этом величина изменения вектора скорости может быть вычислена по уравнению:
Av = 2v sin (Ф/2) , (23)
где ф — угол, на который нужно повернуть плоскость орбиты (рис. 4).
Рис. 4. Изменение направления движения КА:
Vl — начальное направление вектора скорости; У2 — конечное направление после поворота на угол ф
Ускорение, нормальное к плоскости орбиты, должно менять не величину вектора круговой скорости, а только его направление. Уравнение (23) справедливо для равностороннего треугольника. Конечно, прямоугольный треугольник не может быть равносторонним, но если Ду ^ 0, то два прилегающих угла можно принять прибли-
зительно равными 90°. Фактически траектория, по которой КА совершает движение между плоскостями, будет представлять дугу.
Примем ф = 0,0134 рад (0,768°), ус = 28 000 м/с. Тогда, согласно
уравнению (23), Ау = 375,2 м/с. Длину дуги можно вычислить,
умножив ф на средний радиус орбиты гм;а = 2 -1011 м, в результате
вычислений получим 2,68 109 м. Разделив длину дуги на скорость, получим время перехода между плоскостями — 82,67 сут. Разделив скорость на интервал времени, получим ускорение — 5,25 10-5 м/с2.
Нормальное ускорение аN можно также найти по известной теореме об изменении момента количества движения: производная по времени от кинетического момента механической системы относительно некоторого неподвижного центра равна геометрической сумме моментов всех внешних сил, действующих на систему, относительно того же центра (рис. 5) [3]:
ёК ё ф ё ф ( Л
_ = к_ = Ег ^ % = , (24)
где ^ — сила, действующая на тело и перпендикулярная плоскости орбиты.
Рис. 5. Изменение плоскости орбиты
Приняв Аф = 0,0134 рад, At = 82,67 сут (7,14106 с), ус = 28 000 м/с, получим ам = 5,25 • 10"5 м/с2.
Известно, что прецессионное движение гироскопа является безынерционным, т. е. оно прекращается при исчезновении внешней силы, его вызвавшей [9]. Поворот плоскости орбиты небесного тела при сохранении его углового момента, логично также рассматривать как безынерционное движение.
Заключение. В работе впервые представлен алгоритм нахождения параметров перелетной эллиптической орбиты между двумя эллиптическими орбитами. Показано, что с увеличением длины большой полуоси перелетного эллипса увеличивается его эксцентриситет, что приводит к росту радиальной скорости КА на перелетной орбите, и, как следствие, увеличению скоростей разгона и торможения.
ЛИТЕРАТУРА
[1] Larson W.J., Wertz J.R., ed. Space Mission Analysis and Design. Dordrecbt, Boston, London, Microcosm Press, Kluwer Academic Publishers, 2005, 504 p.
[2] Хартов В.В., Ефанов В.В., ред. Проектирование автоматических космических аппаратов для фундаментальных научных исследований. В 3 т. Т. 2. Москва, Изд-во МАИ-Принт, 2014, 544 с.
[3] Митишов Е.А., Берестова С.А. Теоретическая механика: статика, кинематика, динамика. Москва, Ижевск, Институт компьютерных исследований, 2005, 17б с.
[4] Драчев М.М., Демин В.Г., Климишин И.А., Чурагин В.М. Астрономия. Москва, Просвещение, 19S3, 3S4 с.
[5] Александров А.Д., Нецветаев Н.Ю. Геометрия. Москва, Наука, Гл. ред. физ.-мат. лит., 1990, б72 с.
[6] Ostrovskiy N.V. Modeling of the celestial body transition from heliocentric orbit to planet-centric. Reports of International astronomical congress "Astrokazan-2011 ", Kazan, August 22-30. Kazan, Kazan Federal University, 2011, pp. 1SS-190.
[7] Chapront J., Francou G. Ephemerides of planets between 1900 and 2100 (199S update). Bureau des Longitudes, Group: Dynamics of Solar System (199б). URL: http://cdsarc.u-strasbg.fr/viz-bin/Cat7VI/S7
[S] Симонов А.В. Использование орбиты ожидания у Марса для полета к Главному поясу астероидов. Вестник ФГУП «НПО им. С.А. Лавочкина», 2010, № 1, с. 14-23.
[9] Каримов И. Лекция 11. Гироскопы. Теоретическая механика: электрон. учеб. курс для студ. оч. и заоч. форм обуч. URL: http ://www.teoretmeh.ru/dinamika9. htm
Статья поступила в редакцию 17.04.2017
Ссылку на эту статью просим оформлять следующим образом:
Островский Н. В. Алгоритм нахождения траектории перелета между двумя эллиптическими орбитами. Инженерный журнал: наука и инновации, 2017, вып. б. http://dx.doi.org/10.18698/2308-6033-2017-6-1666
Статья подготовлена по материалам доклада, представленного на XLI Академических чтениях по космонавтике, посвященных памяти академика С.П. Королёва и других выдающихся
отечественных ученых — пионеров освоения космического пространства.
Москва, МГТУ им. Н.Э. Баумана, 24-27января 2017г.
Островский Николай Владимирович — канд. техн. наук, доцент Вятского государственного университета. Автор более 30 работ в области небесной механики и теории гравитации. e-mail: onv1@yandex.ru
Algorithm for identifying the transfer trajectory between two elliptic orbits
© N.V. Ostrovskiy Vyatka State University, Kirov, 610000, Russia
This article for the first time ever presents an algorithm for the transfer elliptic orbit between two elliptic orbits. We have solved the following problems: constructing the ellipse in two radius-vectors emerging from the same focal point, calculating acceleration and deceleration velocities impulses, computing the deceleration needed for the orbit plane change. This study is of great importance because the other well-known publications either do not consider such problems or do not suggest particular computation algorithms. We have calculated the transfer orbit length by the numerical quadrature method. The motion along the elliptic orbit is viewed as the superposition of the circular and radial motions. It is shown that with the increase of the transfer ellipse semi-major axis length its eccentricity also increases. It results in the increase of the spacecraft radial velocity in the transfer orbit whereupon the acceleration and deceleration velocities also increase.
Keywords: transfer orbit, elliptic orbit, flight trajectory to Mars, numerical quadrature method
REDERENCES
[1] Larson W.J., Wertz J.R., ed. Space Mission Analysis and Design. Torrance, California, Microcosm Press Publ., Dordrecht, Boston, London, Kluwer Academic Publ., 2005, 504 p.
[2] Khartov V.V., Yefanov V.V., ed. Proyektirovaniye avtomaticheskikh kosmich-eskikh apparatov dlya fundamentalnykh nauchnykh issledovaniy. V 3 tom. Tom 2 [Designing unmanned spacecrafts for fundamental scientific research. In 3 vol., vol. 2]. Moscow, MAI-Print Publ., 2014, 544 p.
[3] Mitishov Ye.A., Berestova S.A. Teoreticheskaya mekhanika: statika, kinema-tika, dinamika [Theoretical mechanics: statics, kinematics, dynamics]. Moscow, Izhevsk, Institute of Computer Science Publ., 2005, 176 p.
[4] Drachev M.M., Demin V.G., Klimishin I.A., Churagin V.M. Astronomiya [Astronomy]. Moscow, Prosveshcheniye Publ., 1983, 384 p.
[5] Aleksandrov A.D., Netsvetayev N.Yu. Geometriya [Geometry]. Moscow, Nauka Publ., 1990, 672 p.
[6] Ostrovskiy N.V. Modeling of the celestial body transition from heliocentric orbit to planet-centric. Proc. of International astronomical congress "Astrokazan-2011", Kazan, August 22-30. Kazan, Kazan Federal University, 2011, pp. 188-190.
[7] Chapront J., Francou G. Ephemerides of planets between 1900 and 2100 (1998 update). Bureau des Longitudes, Group: Dynamics of Solar System (1996). Available at: http://cdsarc.u-strasbg.fr/viz-bin/Cat7VI/87
[8] Simonov A.V. Vestnik NPO im. S.A. Lavochkina — Bulletin of Lavochkin Science and Production Association, 2010, no. 1, pp. 14-23.
[9] Karimov I. Lektsiya 11. Giroskopy. Teoreticheskaya mekhanika [Theoretical mechanics. Lecture 11. Gyroscopes]. Available at:
http ://www.teoretmeh.ru/dinamika9.htm
Ostrovskiy N.V., Cand. Sc. (Eng.), Assoc. Professor of Vyatka State University. Author of over 30 works on celestial mechanics and gravitation theory. e-mail: onv1@yandex.ru
