Математическое моделирование частной ограниченной задачи трех тел с учётом второй зональной гармоники в геоцентрической экваториальной системе координат Текст научной статьи по специальности «Механика и машиностроение»

УДК 517.9:521+523.3:629.7

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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЧАСТНОЙ ОГРАНИЧЕННОЙ ЗАДАЧИ ТРЕХ ТЕЛ С УЧЁТОМ ВТОРОЙ ЗОНАЛЬНОЙ ГАРМОНИКИ В ГЕОЦЕНТРИЧЕСКОЙ ЭКВАТОРИАЛЬНОЙ СИСТЕМЕ КООРДИНАТ

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

Математическое моделирование, баллистический анализ, точка либрации, система «Земля - Луна», ограниченная задача трех тел, вторая зональная гармоника, космический аппарат (КА)

Yu.A. Okishev, Yu.V. Klinaev

MATHEMATICAL MODELING OF A RESTRICTED THREE-BODY SUBPROBLEM WITH ACCOUNT TO THE SECOND ZONAL HARMONIC IN THE GEOCENTRIC EQUATORIAL COORDINATE SYSTEM

A disturbing acceleration from the second zonal harmonics caused by the noncentrality of the Earth's gravitational field is considered as part of mathematical modeling of the ballistic flight from low-Earth orbit to the libration point L1 of the «Earth -Moon» as a solution to the restricted problem of three bodies. The obtained results are updated and are close to the real flight, the value of the total impulse velocity and specifications starting point. The resulting model can predict the required total impulse velocity and, consequently, the minimum needs of the fuel mass at any contact time with the libration point.

Mathematical modeling, ballistic analysis, libration point, the system «Earth-Moon», restricted three-body, second zonal harmonic, spacecraft

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

’ ( АТЛ V

1 - exp-

AVS

y у

(1)

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

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

— силы тяготения Луны, Солнца и других небесных тел;

— аэродинамические силы, которые чаще всего учитываются при исследовании движения низких спутников Земли;

— гравитационные силы, связанные с «неправильностью» фигуры Земли или, как еще принято называть, нецентральностью гравитационного поля Земли;

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

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

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

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

2. Исходя из энергетики перелета и необходимости поворота плоскости орбиты, рассматриваемой эпохой выбран 2025 год, т.к. в этот год наклонение орбиты Луны максимально и составляет

28,443°. Таким образом, угол поворота орбиты КА ^ минимален, и, исходя из (2), второй импульс скорости также будет минимальным

А V = ^V2 + vKali - 2 -Vu • V^ • cos А , (2)

где VL1 - вектор скорости точки либрации L1, VKА - вектор скорости КА в точке либрации L1, Ai -

разница между наклонениями перелетной орбиты и орбиты точки либрации L1.

3. Существует два способа попадания в точку либрации - при перелете из окрестности восходящего или нисходящего узла базовой орбиты.

4. Оптимальной датой попадания в точку либрации L1 является 12 часов дня 13 апреля 2025 года, при этом время перелета составляет 4,25 дня. При этом оптимальная траектория реализуется при перелете из восходящего узла орбиты.

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

Очевидно, что в нашем случае аэродинамические силы слишком незначительны, т.к. минимальная высота в рассматриваемом случае равна 300 км, где атмосферное влияние практически отсутствует. Рассмотрим влияние нецентральности гравитационного поля Земли.

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

d2r . mMз

m л? = f—, (3)

где m - масса КА, r - радиус-вектор КА относительно центра Земли, t - время, f - гравитационная постоянная, Mз - масса Земли.

Реальная Земля существенно «сложнее». Прежде всего, на гравитационное ускорение влияет полюсное сжатие Земли. С определенной степенью точности можно считать, что сечение Земли плоскостью, проходящей через её ось, имеет форму эллипса с большой полуосью, равной экваториальному радиусу Язэ = 6378,245км, и малой полуосью, равной полярному радиусу

Rn = 6356,863км.

Коэффициент сжатия Земли рассчитывается как

R — R 1

а = -^-----^ ---------= 0,0033523 (4)

R„ 298,3

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

U = §(3sin2 р-1) (5)

r 3r

где m - гравитационный параметр Земли, r - расстояние от центра Земли до точки, потенциал в которой рассматривается, р - широта точки, § = 66.07 • 103 км2 - константа.

Первое слагаемое справа в равенстве (5) - ньютоновский потенциал. Второе слагаемое учитывает сжатие Земли по полюсам.

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

Для записи возмущающего ускорения, вызванного нецентральностью гравитационного поля Земли, на оси орбитальной системы координат воспользуемся соотношениями, вытекающими из рис. 1.

На этом рисунке положение спутника характеризуется точкой А; плоскость орбиты спутника есть EAF; меридиональная плоскость, проходящая через точку А, - BCD. Угол между местным меридианом и трансверсальным направлением в точку А обозначен у, угол ЕСА - и (аргумент широты), угол ABC - ф - широта точки.

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

S = т § • (3 sin2 р-1) = m § • (3-sin2 /'• sin2 и -1) r r

LI • § LI • §

T =-----3—~ • (sin 2p • cos y) = ——— (sin2 i • sin 2u) (6)

r r

W = L § • (sin 2p • sin y) = L § • (sin 2i • sin u) r r

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

sinp = sinu • sini, cos y = ctgu • tgp, (7)

где 1 - наклонение орбиты, u - аргумент широты.

В уравнениях (6): S - радиальная, T - трансверсальная, W - нормальная составляющая возмущающего ускорения от нецентральности гравитационного поля Земли.

Для того чтобы использовать алгоритмы анализа, описанные в [2-6], необходимо выразить возмущающее ускорение от нецентральности гравитационного пол Земли в геоцентрических координатах x, y, z.

Для этого существует два способа решения.

Первый способ - выражение влияния нецентральности гравитационного поля с учетом изменения оскулирующих элементов орбиты от времени: p(t), e(t), W(t), i(t), w(t), tp(t) - фокальный пара-

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

Используем матрицу перехода от геоцентрической системы координат (x,y,z) к орбитальной (n,r,b) A(8), где г - радиаль, n - траснсверсаль, b - нормаль.

- sin u cos W - cos u sin W cos i - sin u sin W + cos u cos W cos i cos u sin i cos u cos W - sin u sin W cos i cos u sin W + sin u cos W cos i

- sin W sin i cos W sin i

A=

sin u sin і - cos і

(8)

где О - долгота восходящего узла орбиты.

Матричная форма перехода имеет вид

С пЛ

С x Л

=A

у

V z У

(9)

Таблица 1

Элементы матрицы перехода между геоцентрической и орбитальной системой координат

Система координат Орбитальная система координат

N r b

Геоцентрическая орбитальная x - sin u cosW- cosu sin W cosi cos u cos W - sin u sin W cos і - sin W sin і

У - sin u sin W+cosu cosW cosi cos u sin W + sin u cos W cos і cos W sin і

Z cos u sin і sin u sin і - cos і

Используя матрицу перехода, возмущающие ускорения от нецентральности гравитационного поля Земли на оси x, y, z примут вид

ФсжХ = T ■ (sin i cos W - cos u sin W cos i) + S ■ (cos u cos W - sin u sin W cos i) + W ■ (- sin W sin i)

Фсж Y = T ■ (- sin u sin W + cos u cos W cos i) + S ■ (cos u sin W + sin u cos W cos i) + W ■ (cos W sin i) ,(10)

Фсж Z = T cos u sin i + S sin u sin i + W (- cos i)

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

Второй способ - учёт производных возмущающего ускорения нецентральности гравитационного поля Земли по осям геоцентрической системы координат (x, y, z).

Ф

Ф

Ф

cж Z

m3 d rx f 5rz2

r5 V r2

m3 • d- Гу f 5rz2

r5 V r 2

m3 • d rz ' 5г'^

r5 V r2

\

-1

У

Л

(11)

-3

В геоцентрической системе координат система уравнений движения с учетом влияния нецентральности гравитационного поля Земли примет вид

<ІУ г и, и -(г — г ) и, •г

х х г^з г*л V х лх / г*л л

dt

V

dt

ddK

dt

r3 (r - Гл )3

rv Мл -(rv - глу ) Мл •rr

+ Ф

cжX

у лу/ Г л лу + ф

r3 (r - Гл )3 Гл3

r М m • (r - r ) М • r

z г* з г* л \ z л / г* л лі

cж Y ,

(r - rz )3

+ ф

cж Z

г

1

3

3

r

r

z

где г (гх, гу, тг)- радиус вектор КА относительно Земли, г л (глх, глу, глг)- радиус-вектор Луны относительно Земли, т - гравитационный параметр Луны.

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

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

Полученные результаты представлены в табл. 2, 3.

Таблица 2

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

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

110 688,535 4,15 -40,874 145,566 3783,3

112 721,418 3,947 -21,483 156,762 3813,2

114 675,625 3,913 -5,81 170,943 3761

116 699,267 3,697 8,398 -170,305 3784,7

Таблица 3

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

Дата попадания с начала эпохи, дней Д м/ с Т °р, дней а,° и ° ДУт . т , м/с

102 638,211 4,213 8,297 10,538 3737,67

102,5 637,653 4,233 11,414 13,962 3737

103 637,59 4,25 14,648 17,285 3736,9

103,5 638,142 4,263 18,037 20,479 3737,3

104 639,45 4,277 21,621 23,506 3738,5

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

3739

С учетом нецентарльности гравитацонного поля Земли

Без учета нецентральности гравитацонного поля Земли

3733

102

102,5

103

103,5

104

Дата попадания в точку либрации

Рис. 2. Сравнительный график зависимости суммарного импульса скорости от даты попадания КА в точку либрации 1_1 системы «Земля - Луна» при перелете из восходящего узла орбиты

На первый взгляд, значения суммарного импульса скорости изменились несущественно, но при расчете по формуле Циолковского (1) необходимо потратить дополнительно около 4 кг топлива. Следует отметить, что стоимость вывода 1 кг полезного груза только на геостационарную орбиту ракетой-носителем «Союз» составляет около $ 25 000. Таким образом, учет нецентральности гравитационного поля Земли позволяет уточнить массу полезного груза и разумно потратить финансовые средства. Также следует учитывать, что коллинеарная точка либрации L1 системы «Земля - Луна» -неустойчивая. Соответственно, необходимо учитывать влияние нецентральности гравитационного поля Земли для минимизации количества топлива для дальнейшей стабилизации КА в точке либрации L1 системы «Земля - Луна» и, как следствие, увеличение массы полезной нагрузки КА.

Полученная математическая модель позволяет спрогнозировать потребный суммарный импульс скорости и, как следствие, потребную массу для перелета топлива в любой момент времени попадания в точку либрации L1 системы «Земля - Луна». Также можно найти оптимальные даты старты, когда энергетические затраты баллистического перелета минимальны.

Авторы выражают благодарность д.т.н., профессору Константинову Михаилу Сергеевичу (Национальный исследовательский университет «Московский авиационный институт») за ценные советы, конструктивные замечания и рекомендации по учету нецентральности гравитационного поля Земли в поставленной задаче численного моделирования баллистического перелёта КА в точку либрации L1 системы «Земля-Луна».

ЛИТЕРАТУРА

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

2. Окишев Ю.А. Математическая модель баллистического анализа перелета космического аппарата в точку либрации системы «Земля - Луна» / Ю.А. Окишев, Ю.В. Клинаев // Математические методы в технике и технологиях - ММТТ-25:сб. тр. XXV Междунар. науч. конф.: в 10 т. Т. 6. Саратов: Сарат. гос. техн. ун-т, 2012. С. 71-74.

3. Окишев Ю.А. Разработка математической модели баллистического анализа перелёта космического аппарата с низкой околоземной орбиты в точку либрации L1 системы «Земля - Луна», как решение частной ограниченной задачи трёх тел с учётом прецессии орбиты Луны / Ю.А. Окишев, Ю.В. Клинаев // Вестник СГТУ. 2012. № 4 (68). C. 61-68.

4. Окишев Ю.А. Анализ влияния положения Луны на выбор схемы космического перелета в зависимости от даты попадания космического аппарата в точку L1 системы «Земля - Луна» / Ю.А. Окишев, Ю.В. Клинаев // Информационные технологии, системы автоматизированного проектирования и автоматизация: сб. науч. тр. IV Всерос. науч.-техн. конф. с междунар. участием. Саратов: Сарат. гос. техн. ун-т, 2012. С. 114-118.

5. Окишев Ю.А. Алгоритм поиска оптимальной траектории баллистического перелета с низкой околоземной орбиты в точку либрации L1 системы «Земля - Луна» / Ю.А. Окишев, Ю.В. Клинаев // Прикладные аспекты исследований в радиофизике и др.: сб. науч. ст. СГУ и СГТУ. Саратов: ООО Изд. Центр «Рата», 2013. С. 50-61.

6. Окишев Ю.А. Основные подходы к численному моделированию частной задачи трех тел для баллистического анализа перелета космического аппарата с низкой околоземной орбиты в точку L1 системы «Земля - Луна» / Ю.А. Окишев, Ю.В. Клинаев // Вестник СГТУ. 2013. № 1 (69). C. 44-49.

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

аспирант кафедры «Техническая физика Postgraduate

и информационные технологии» Department of Technical Physics

Энгельсского технологического института and Information Technologies,

(филиала) Саратовского государственного Engels Institute of Technology (Branch)

технического университета имени Гагарина Ю A. Yuri Gagarin State Technical University of Saratov

Клинаев Юрий Васильевич - Yuri V. Klinaev -

доктор физико-математических наук, Dr. Sc., Professor

профессор кафедры «Техническая физика Department of Technical Physics

и информационные технологии» and Information Technologies,

Энгельсского технологического института Engels Institute of Technology (Branch)

(филиала) Саратовского государственного Yuri Gagarin State Technical University of Saratov

технического университета имени Гагарина Ю.А.

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