Описание движения космического аппарата с помощью дифференциальных уравнений Текст научной статьи по специальности «Математика»

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

Федоренко Ю.В., Аксенов С.А.

Национальный исследовательский университет «Высшая школа экономики» jufedorenko@gmail.сот. aksenov.s.a@gmail.сот

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

Ключевые слова: точки либрации, точки Лагранжа, траектория Лиссажу, математическое моделирование, компьютерное моделирование

1 Общие сведения

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

Существует два семейства орбит вокруг точек либрации: гало-орбиты и орбиты Лиссажу. Оба семейства орбит имею как преимущества, так и недостатки. Так, например, если требуется непрерывно наблюдать космический аппарат с Земли, то необходимо использовать гало-орбиты, так как аппарат, находясь на таких орбитах, никогда не скроется в тени Луны. Недостатком гало-орбит являются большие значения амплитуд (47000км и более), что не позволяет качественно осуществлять связь, например, с портативной станции на обратной стороне Луны. В этом случае необходимо использовать движение аппарата по орбитам Лиссажу, амплитуда которых может быть в десятки раз меньше гало-орбиты.

2 Орбиты Лиссажу

Орбиты Лиссажу являются квазипериодическими потому, что период колебания по оси Y не совпадает с периодом колебания по оси Z. Это различие в периодах в результате дает сложное гармоническое движение, описанное кривыми Лиссажу [Cuevas, 2002]. В дополнение, период для каждой оси является функцией от собственной амплитуды.

Бесконечное количество комбинаций периодов на осях Y и Z соответствует такому же обширному спектру орбит Лиссажу, каждая из которых остается центрированной относительно точки либрации в плоскости XY. Орбиты Лиссажу могут претерпевать затмение Луной, что означает затруднение связи с Землей. Подобные промежутки времени без сообщения с Землей могут продолжаться от нескольких часов до суток и даже больше. Также, асинхронный характер колебаний означает, что кривая может не пролечь на том же самом месте несколько сотен дней, что делает проектирование миссий гораздо более сложным.

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

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

Приближенная аналитическая оценка отношения времени затмения аппарата Луной ко времени нахождения на орбите была предложена Р.Фаркуаром в [Hopkins, 2013]

где

Ам = 3 099 км? Ау1 = Аг1 = 0,02 = 7688 км_

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

3 Постановка задачи и аналитическое решение

Если рассматривать движение в декартовой системе координат с центров в точке Ь2, то уравнения движения имеют вид

х — 2(1 + 1?г)у + (1 + уь — 1?гу =

С1 + 1?г)2[(1 + ГГ1)(1 + д) + х]-~(1 + -

1-/4

^ла+гг^а +«?) + х] - %

>Ь Уь

П

т2а| Уь

(1-м+?1г)С1+е)+у£ Сх-х£

+

г( (.1.4

'о

йо

у + 2(1 + %)[(1 + П1)^^ + ^[(1 + ГГ^О + С?) + -г] - + == + (1 + г^)2]у -

1-Е |

у1

2 3

У-Уз | Уз

'О

я.

о.

1)

г 4-

+ ^У = УхЖ - рх(1 + г?г)[(1 + у^ 1)(1 -

+ л ] -

1-м

~ Л»3

2 3

2 ™ .

_ 1.

¡1 № Ц2 зо

- массовый параметр для системы Земля-Луна ( ).

Уь

- отношение расстояний для точки Ь2

У = с1 /И = 0.1678331476 ( )• т

- соотношение между средним движением Луны и Солнца (т=0.0748013263).

Г Г ЛЛ /у*

©' *ф'((I

- расстояния между малым телом и Солнцем, Землей и луной соответственно.

- расстояние между Землей и Солнцем.

П

- расстояние между барицентром Земля-Луна и Солнцем

Уз*

(обозначение компонент " " ).

С£.3

- среднее расстояние между барицентром Земля-Луна и Солнцем.

2

Переменные т , используемые в уравнении (2), связаны с прямым солнечным возмущением. Значения переменных Ч, рх и зависят от

косвенного солнечного возмущения и от эксцентриситета орбиты Луны. Пренебрегая членами более высокого порядка, уравнения движения можно записать в развернутом виде [Hopkins, 2013]

х - 2(1 + vz) у = [(1 + vz)2 + 2(1 + q)~3Bl]x + (уz)y ~ vxvzz - -\cL{ 1 + q)~4[2x2 -(y2 + z2)] + 2Dl(1 + Q)-5 [2x2 - 3(y2 + z2)]x ---El( 1 + q)~6 [8*4 + 3(y4 + z4) - 24x2(y2 +

z2) + 6yV] _ _m2 (b) {[! _ 3 (I)'

x —

у + (1 + vz)x =[(14- vz)2 -(14- q)-3Bl + V*]y -vzx -I- vxz -I- 2vxz + +30,(1 + q)

xy-

5

-DLQ1 + e)"s[4x2 - 2(y2 + z2)]y + f^(l + +a)"6 [4х2 - 3(y2 + z2)]xy +

m2ib -4m

X

- я fe

1-3

\rs

y-

(2)

Z-"= \Ух - (1 + q)~3Bl]z - vx(l + vz)x - vxy -2vxy+3CL(l + -\dlH + qT5№2 -

(y2 + z2)]z + Iel(1 + Qy6 [4x2 - 3(y2 +

где bl = 3.190423657 CL = 2,659334398 DL = 2.5830098И EL = 2.572040953 и

(1 + ^) 1 = 1 4 е соз ф 4 е2 соз 2ф 4 ^ш2 4

т2 соз 4 те--ф) — ^е3 соз ф 4

9 з 19 з те 15 (а\ г

-е соэЗ ф Н—ш со5 2£--т. —). соз < -

8 г 6 * 16 \а'/ >

—-еу2 соз(ф- — 217) -I- —те2 соз2£ 4 ^■т2е соб(2% — ф) 4 4 шее' соз(2^ — ф -

15

ф') Н—шее' — ф 4 ф') Н-

§

4

-(-7) е*4 ф') — — т2£?со5^ + ^

4 \_dT- ^ 12

— ш2есо5(2^ + ф)--—те2соя(2£ — Зф) —

16

3 2 I л. ь I ^ 2

2

— тее'со^(ф 4 ф')

-т2е' созф' 4 -т2е'соБ(2£ — ф') — — -т2е' соз(2<^ 4- ф') 4 т^тее'соб(Ф - ф') —

- 2е со б ф 4 - е2 соб 2ф 4 — т2 шз 2\ 4

2 4

— те соз(2^ — ф)--— соз ф 4 — е3 соз Зф 4

85 з 15 / а \ г 5 2

— т соз 2с--т — I соз Ь — соз( ф —

12 ъ е \а'/ ъ 4 г ^

—2г||) + — те2 соз 2\ 4 — т2е соз(2^ — ф) 4

— тее'соз!^ — ф — ф')--— теег соз(2^ — ф

ф') 4 ^ е' соз(% 4 ф') — ^ т2 е соз ф 4

4 ^гт2е соз(2^ 4 ф) — ^те2 соз(2^ — Зф) —

71

Зт2е'соз ф' 4 4 — т2е' соз(2^— ф') —

— т2е' га (2(4- фО 4 — тее' соз(ф — ф') —

— — тее' соя(ф 4 фО

4

г?х - Зт2у со б % зт(уг — П) при ЭТОМ

ф = ^т2 - ^р-т3^ £ 4 £ — ¿6

(4)

(5)

f = (1 - m)t + Е - £f (7)

(л _L 3 2 9 3 273 VU ■ Г»

п =1 4 --m15--m t 4 £ — íin ,n4

' V 4 32 1.28 / u (8)

ф'' = mt 4- £.' — й/ (9)

c;

г/ = mt + £' 4- 2ef sin<&f 4- sin 2ф* /1 „ч

4 г (Ш)

Í2:

гч f3 2 9 3 273 Л*

= íin — -m--туп--m £

и 44 32 1.28 /

(П)

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

(j^j: = cos f - 2(е sin ф - ef sin фf) sin % —

cósale2 4- el2-y2--(e2 eos 2ф 4- e'2 eos 2 ф') -

4

2ee'[cos((f>- фг) - соз(ф 4 ф')]} "

—sinf ¡~e2 sin 2ф — ~ef2 sin 2ф' 4- —m2 sin Zf 4 (12)

14- 4 £>

—mesin(2$ — ф)--3те' зт ф'} + -y2 cos(2r¡ —

4 J 4?

O

{^j = — sin £ — 2(e sin ф — e' sin фf) cos £ 4-sin^ {e2 4 ef2-y2--(e2 eos 2ф 4- ef2eos 2ф') —

l 4

2ee'[cos(ф — — cos(</> 4 ф')]] ~~ —cos^ j-e2 sin2ф - -e'2 sin 2ф' 4- -^m2 sin 2f -I-

4 S

— mesin(2f — <£■)--Зше' sin — -y2 sin(

4 J 4

f)

(13)

= т(У-П) (14)

г^)3 = 1 + Зе'соБф' + -е,2(1 + Зсоз2 ф1) (15)

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

е - эксцентриситет орбиты Луны (е = 0,054900489^ е - эксцентриситет орбиты Земли = 0,0167217) (а./а ) _ изменение отношений между большой полуосью орбиты Земли и большой полуосью орбиты Луны [Ся/а ) =0.0025093523]. У - тангенс среднего наклона орбиты Луны (У = 0,0900463 06

£,£ _ средняя долгота эпох среднего движения Луны и среднего движения Солнца.

¿6,6/ _ средняя долгота лунного перигея и средняя долгота солнечного перигея.

^о - долгота среднего восходящего узла лунной орбиты. Поиск решений приведенных выше уравнений движения осуществляется с помощью вычислений на ЭВМ. В результате многочисленных алгебраических преобразований, можно получить точные аналитические представления некоторых семейств квазипериодических орбит вокруг точки Ь2.

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

х — ау — Ьх = О

у + ах + су = 0 (16)

М + = 0 где

а = 2^1 + -у2т2) = 2.000034

;п)

b = В, [2 4 Зе2 4 т2 4 — е4 4 —т2е2 - — т* + L V 4 64 1.6

-mVz) 4 + f 1 + 2е2 4 — 4 - е4 + — mV 4

2 } V 2 8 32

— m4 4 -mV2l = 7.436984

32 4 /

j-j, i 3 2 , . 15 4 i 483 2 2 9 4 ,

с = B, [ 1 4 - 4--h — H--mre¿--mr 4

L V 2 2 8 128 32

3 2 гг\ /-1 I i 2 I I 21 4 I 129 2 2 I

-7ТГС *--14 2e¿ 4--— --m¿ez 4

4 J \ 2 8 32

—ra4 + -mV2) = 2,204904

32 4 /

л n (Л . 3 2 , , 483 2 2 9 4 i d = В, 1 + -e 4--h —m^e*--rcr 4

** \ 2 2 128 32

^raV2) 4 ra2 (l 4 V2 - - ¿y2) = 3.219481

Параметры ТП, в, у И е' использовались в вычислениях выше для упрощенного анализа уравнений более высокого порядка, а также для вычисления более точных значений для основных частот. Легко увидеть, что в уравнении (17) решение по оси Z есть простая гармоника и не зависит от решений по X и Y. Движение в плоскости XY является согласованным и характеристическое уравнение может быть представлено в виде:

s4 4 (а2 — b — c)s2 - Ьс — 0 (18)

Это уравнение имеет два действительных корня, равных по величине, но противоположных по знаку. Два оставшихся корня являются чисто мнимыми. Таким образом, движение в плоскости XY имеет экспоненциальный и осциллирующий режимы. Однако, при соблюдении определенных начальных условий, будет реализован только колебательный режим и решение будет ограниченным [Farquhar, 1970]. В этом случае линеаризованное решение принимает вид:

х10 = кАу0 sin(>xy01 4 0j_)

Ую - Ayo cos(>xyü t 4 (19)

z10 = Az0 sin(ft>z0 t 4 в2) где

й)ху0 ~ 1.865485

= 1-794291 (20)

к = Шзв- = £1^ = 0.341763

и четыре константы обозначаются как Ау0, б1и б2 и определяются начальным положением и скоростью аппарата.

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

х — тх1 4- тгх2 4 т3х3

у - ту1 + т2у2 + т3у3

(21)

г = тг^ 4- т2г2 + ш3г3

Асимптотическое разложение здесь произведено по малому параметру т. В данной схеме, первые слагаемые в правых частях уравнений (22) принимаются равными 0(т) , следующие слагаемые принимаются равными 0(т2) и так далее. Константы е/е' также задаются как 0(т), а а/а' приравниваются 0(т2).

Уравнения (22) подставляются в уравнение (3), чтобы получить систему линейных неоднородных дифференциальных уравнений в представлении (xi, yi, zi), i = 1, 2 , 3 ... . Это осуществляется посредствам приравнивания коэффициентов при одинаковых степенях m за исключением постоянных коэффициентов линейных членов, которые описаны в уравнениях (17) и (18),

Чтобы избежать появления секулярных членов в (хЗ , уЗ , z3 ) может быть использован метод Линдшгедта - Пуанкаре [Moulton, 1914]. В этом методе время нормализуется так, что допускается зависимость частоты колебаний от амплитуды. После такой процедуры (xl , yl , zl ) можно записать в виде:

х1 = 0.341763 Ау sinTY

Ух — Ау eos Тг

(22)

zL = Az sin T2j где

+ ®i (23)

Т2 = Mz{AyAz)b + при

шжуО

(úXy —

0)z =

1 + m2ú}xyZ(AyfAz)

(24)

1 + m¿(i>z2(Ay,Az) Дифференциальные уравнения для (х2'У2/г2) имеют вид:

d2x2 [ а \ dy2 Ь _ F-

~2 I I ? — ~

йТг \Wxyo/ йТг wxyü wxyo

2

d2y f a \ dx2 с G

(—)

\Wxy о/

+ — )^Г + 35-У2 =

2

dT2 \ W-rvn í dTi W^n У ¿ W2

i (25)

£Í2Z2 , H2

dT2 2 w|0 где

F2 - 2(3BL + 2)xt eos ф — 2 yx sin ф + 4 (m) drf COS ^ ~ ~ 2 ~ X? ~ Zí)

Gz = 2: ^ sin ф - Q (3BL - 4)yx eos ф - (26)

4(m)drf+ +3C¿Ai>'i

H2 = — 3 j BLz± eos ф + 3 Решая уравнения (26) получаем

= G.5549Ü4Í—UV sin(ф - ГЛ

+ 0.49321 з(—) sin(<¿> Hh Т,)

Vm/ J

—0.09588405 Ау соэ 2Тг + 0.12:8774 А22соъ 2Т2

- 0.26186 А2 --0.205537 Ау

у2 = -1.90554 (£).Ау соз(ф - Гх) +

1.210699 соз(ф + Гх) —

-0.055296 А2 эш21\ - 0.08659705 А*&т2Т2

= 1.052082 (£)А25т(ф + Т2) +

1.856918 (£] Аг - Г2) +

+0.4241194 АуАг сок(Т2 - Гх) -0.1339910 АуАг со5(Т2 + Т^

Силовые функции для уравнения третьего порядка имеют вид

Я, =

2шху2(Ьк - - ш2у0

0342:854 Г-) +

\т/

0.008202 А2 + +0.108050 А2

| Ау si.ii Тх +

6д - -■ 2тху2(аксоху0 + с) -

0.325455 (£)" - 0.058141 А2

хуО

+0.102210 А2

+

Ау соз Тх + •••

Я-> = ш

1/196357 (£)2 + 0.078465 А2-

0.065877 А2 - 2^2

Аг • ■ эш Т? +

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

íTl • 'TI1 fTl *ТТ »TI • 'TI

секулярных членах, содержащих 1i. sm 1i. cos 11 и 1 г sm Эти члены

могут быть устранены путем правильного выбора значений шху2, 0)zi и

произвольной константы Q, которая вводится, чтобы хз содержал члены

СгАу sin Ддадиз приводит к трем алгебраическим уравнениям, которые используются для определения условий существования данных констант:

Ct = 0.09210089 (£) + 0.02905486 А2 + 0.007644849 А\

шкуг = 0.1387811 (£)2 + 0,04349909 Ау - 0.04060812 А\ ^ o)z2 = 0.5981779 (£f - 0.03293845 A¡, + 0,03923249 А\

Конечный результат можно представить в виде

х3 = (£): Ау[—0.122841 sin(2</; - TL) +

0,643204 sin(2</> + Fx)] +

+ (£) A¡[0.198:388:cos </> - 0,387184cos(ф -

2 T2) + 0.335398: cos (ф + +2Г2)] +

(j¿) A2y[0.173731 cqs<¿> + 0.32:5999 cos(<£- 2TJ -

—0,270446cos(</) + 2TJ] +

Q: Ay[-1.10033 sin(<£- TL - 2f) - (30)

-1.189247 sin (ф+Т1 - — 2f)] + AyÁ2[—0.430448sin(2T2 - TL) -—0.0313 0 2 sin ( 2 T2 + TL) ] Hh

[0.027808sin ЗТг] + C\Ay sinT^ + —/1,.[—0,38856sin(Тг - 2:{) Hh +0.455452:sin(Tx + 20]

Уз = [0.608685 соз(2ф-7!) +

1.407026 сог(2ф + Г1)] +

+ (£),Л1[-0.116822зтф - 0.214742:51п(ф -

2Т±) - 0.232503БшСф + +2Т2)] +

(£)Лу[-0.109499 зш ф - 0.144553этО -

27^) - —0.155751 ът(ф + 27^)] +

(£)Лу[2.733376.соэСф - - 20-

-3.848485соз(ф + 7\ - 20] +

ЛуЛ§[—1.191421 соз(2Х2 - 7\ ) -

-0.000165 003(7*2 + т1.)] +

Лу[—0.027574соз ЪТг] + +Ау[-1.743411 соб(Тх -

20 + 0.741825 соз^ + 2{>]

= ^[0.5366:52 зт(2:ф-Г2) +

1.103381 зт(2ф + Г2)] +

+ (¿)АуД2[-0.353754оо5(ф - Т2 ~ 7\) +

0.367360 соэСф + Т: + Г2) + +0.063629 со5(ф + Т2 - 7\) ~ 0.034729 собСф + Т2 - Т^] +

+ (£);Лг[-2.353465 5т(ф - Т2 - 20 "

3.831413 зт(ф + Т2 - 20] + +Л|[0.0176645тЗГ2] + Л^Л2 [-0.86684зт(Г2 -

271) "

—0.044724 5т(Т2 + 271)] + Лг[—1.487917 5т(Т2 - 20 + +0.475507зт(Г2 + 20]

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

4 Список литературы

[Cuevas, 2002] О. Caevas An Overview of Trajectory Design Operations for the Microwave Anisotropy Probe Mission // AIAA Paper 2002-4425, presented at the AIAA/AAS Astrodynamics Specialist Conference, Monterey, California, August 5, 2002.

[Burns, 2013] Burns, J. O., Kring, D. A., Hopkins, J. В., Norris, S., Lazio, T. J., and Kasper, J., "A Lunar L2-Farside Exploration and Science Mission Concept With the Orion Multipurpose Crew Vehicle and a Teleoperated Lander/Rover," Advances in Space Research, 52 (2013) pp. 306-320.

[Hopkins, 2013] Hopkins, J.B. Proposed orbits for human missions to the Earth-Moon L2 region // J. Hopkins, Lockheed Martin Corporation, United States // 64th International Astronautical Congress, Beijing, China, 2013

[Blaquiere, 1966] Blaquiere, A.: 1966, Nonlinear System Analysis, Academic Press, New York.

[Farquhar, 1970] Farquhar, R.W. The Control and Use of Libration-Point Satellites / Robert W. Farquhar // Goddard Space Flight Center Greenbelt, Maryland 20771, 1970.