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

_______УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Т о м II 197 1

№ 5

УДК 521.4:629.76+629.78

О СБЛИЖЕНИИ КОСМИЧЕСКИХ АППАРАТОВ ПО МЕТОДУ СВОБОДНЫХ ТРАЕКТОРИЙ

Ю. П. Яблонько

Рассматривается задача сближения двух космических аппаратов, находящихся на близких эллиптических орбитах, при использовании метода свободных траекторий. Предложены два метода приближенного решения уравнений относительного движения с переменными коэффициентами: асимптотического интегрирования и рекуррентных соотношений. Определены условия сходимости решений. На основании полученных решений найдены аналитические соотношения, используемые для определения характеристических скоростей, необходимых для реализации сближения.

Одним из перспективных методов управления при сближении космических аппаратов(КА) является метод импульсного корректирования кеплеровой траектории КА.

Частным случаем этого метода является двухимпульсный переход, называемый иногда методом свободных траекторий. При его реализации корректирующие импульсы могут быть определены либо интегрированием уравнений относительного движения [1], [2], либо варьированием интегралов уравнений движения [3].

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

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

Метод асимптотического интегрирования, аналогичный методу малого параметра Пуанкаре [4], [5], позволяет получить решение при малых значениях эксцентриситета опорной орбиты для траекторий дальностью порядка одного оборота. Метод рекуррентных соотношений [6], основанный на отыскании решения в виде ряда по степеням независимой переменной, применим при любых эксцентриситетах е<1, но для траекторий с ограниченным временем сближения.

Относительное движение двух КА рассматривается в подвижной правой ортогональной системе координат xyz, начало которой совпадает с центром масс одного из них, двигающегося по известной опорной орбите. При этом ось у направлена вдоль радиус-вектора аппарата от центра притяжения, ось х перпендикулярна оси у в плоскости опорной орбиты. Уравнения относительного движения в случае линеаризованного гравитационного поля имеют вид [1]—[3]

х — с2 erf cos ft* — 2erf у + 2c3 erf sin fty = 0;

у — с1 -г)3 (3 -f e cos ft)_y -f- 2crf x — 2c2 erf sin & x = 0; (1)

z -|- с2 rf z = 0;

здесь с = ^l2/p3l2\ = 1+4? cos jj. — гравитационная постоянная;

e, ft, p — эксцентриситет, истинная аномалия и параметр опорной орбиты.

Метод асимптотического интегрирования. Вводя переменные Х\ = X, Х2 = х, х3 =у, х4 =у, (2)

представим первые два уравнения системы (1) в нормальном виде:

dX = U (t)X,

где

U{t):

0

С2 erf COS ft 0

2с2 erf sin 8

dt

1

0

0

-2crf Х =

о

—2с2 erf sin ft 0

с2 (3 ~f~ е cos ft) rf

X,

x2

Xg X.

0

2crf

1

0

(3)

(4)

(5)

Третье уравнение системы (1) является независимым, и его решение может быть найдено отдельно.

Следуя [4], [5], введем вместо (3) систему

ах

dt

= U(x)X;

зде сь

t/(T) = £e*i7*(0 = t/e + ei;1(0 +

(6)

(7)

*=о

т = е£, е — малый параметр.

Будем искать решение с учетом двух членов разложения, что дает достаточную точность при е = 1, так как в этом случае системы (3) и (6) совпадают. Поэтому в дальнейшем

£/(*) = £/(*)= */0+*/,(*). (8)

Следует подчеркнуть, что истинным малым параметром системы (3) является эксцентриситет е опорной орбиты, а не введенный параметр е.

Матрица 6/0 получается из 1Щ) при « = 0:

и& —

0 1 0 0

0 0 0 2 с

0 0 0 1

0 -2с Зс2 0

Разобьем собственные числа матрицы и0 на р групп , 4а’ гДе

Р

о=1,.. „р; 2^а = /1 (л = 4 для матрицы и0)-

о— 1

Предположим, что выполняется условие

1^')-^*)1>с>0 (»**; /=1,.. ,Й0;У=1.......**>•

Определив собственные числа матрицы (9),

^1,2 = ^3,4 = зЬ с • *’>

разобьем их на две группы (р = 2). Представим матрицу Последующим образом:

р

ио^'ЩК'А'М,, (П>

0=1

где /Са.'Лв> Л1а —матрицы типа пХ^а, £яХ*а и соответствено. Будем искать

решение в виде

*(, = £вУ. = (*„ + К111)1'в (а = 1,2), (12)

где Уа — вектор-функция, являющаяся решением уравнений ЛУа

-2Г “ А» К. = (Д. + Л?1) У. (в=1-2>- <13>

Как показано в работе [5], матрицы Ка можно получить из вспомогательных матриц вида

Р *5

М<Л>)-ПП<У®-Х)')£*>. (14)

5=1 у—1

где Еп — единичная матрица п-го порядка. В рассматриваемом случае

*1 ((У0)^(и0~Х3Е)(и0-Х4Е);

А (и0) = (и0-к Е)(Ц0 -\,Е).

В качестве матриц /(а можно взять любые матрицы, составленные из линейно независимых столбцов матриц Да(Уо). Подставив (10) в (15), найдем

О

3 с О 2

О О

>.)

(15)

м = к~1 =

мг м»

Л = миоК =

А 0

0 Л2

0 2

2с 0

1 0 ’

0 —с

0 1 с 2

1 0 0

0 2 с -3

0 0 0 —

О О о

3 С 0 | 1

0 1 0 — с 1

| с 0

(16)

Л!» = - М,йМ (» = 1,2),

где = (/)/с.;

так как Ка и Ма определяются из (16), то

ЛК,

Л-=о и д1°1 = -£/,(0^; | К[ц = /со'а11, \

(18)

где

Полагая (?^У = 0, остальные члены найдем из уравнений А, <2$ = К + м, О™ (5 ф а; а = 1,2).

(19)

Запишем в развернутом виде дифференциальные уравнения вспомогательной системы (13):

^- = <а1 + л(Ч)Г,-14“’1'!|п8

— сегр СОБ 9 с [3 + 4т)8 (г;е соз 0 -(- Зт) — 3)] 4сег,3 51п в

ау2 /Л I *ИК V II —4се7)38*П» с[3+ 4*1*(1]гСО$& — 1)]

1Г = (Л2 + Л2 ) У2 =_________________. ' _

Ух

У2

'* ‘ ' 1 ' " || — сгр (у\е соя $ + Зк) — 4)

Систему (20) решаем методом Пикара

• 4сет)3 81П 0

(20)

где <?! =

с\ 1 — || «3 и с2 =

с21! II «4

( *

К] = Есх + | Л! С! <Н, У 2 = Есч + ^ Л2 с2 ЛЬ, о о

■ произвольные постоянные.

(21)

У =

Ух *11 *12 0

У2 *21 *22

Уз У4 0 *33 *34 *43 *44

Вс,

(22)

а / в

где Ьп = Ьп — 1 + 2е (е эШ2 *> — 2 сое й)9о; Ьп = — е (^п 8- + -у + вШ 2»д ; Ь21 = 116е з1п & -+- 2е2 2^) '

' ------Чз7Гагс‘8(

Зе 81 п 0 (1 —ег)(1 -(-ЙСОБ^

(1 —е^у

1 — е »

уг^7* ‘8 2

*зз = *44 = 1 — 2е (е ®!Г12 8—2 соб9)^; Ьи — ^4е вШ 0 + 2е2 у вт 20

/ 1 — е »

1 \ УТ^ё* ^ 2

Зе 81*1 д

~ (1 — е2) (1 + е сое») + (1 — еа)3/2 агс,8

*43 = |& — Зе вш а — е (в!п » -у- + -|- вт 2»^а; В = || &у ||; с

Решение системы (3) имеет вид

X = KY =Ш=\\аи\\-\\Ьи\\

■ с (i,J = 1.......4),

где

ап = — 28ет,8 sin 0; я13 = — 14etj8sin 9; я21 = с [3 + 4тг)2 (3-г) -f iet] cos & — 3)]; а2з = 2с [4 -J- Зт)2 (f\e cos & — 1)];

«3i = 2[l+r)2(4erj cos & -(- Зт] — 3)]; а33 = 5 + 4г)2 (к)е cos & — 1); й41 = 14сег,з sin 9;

Окончательно решение запишем в виде

* = 1ИуИ<7-[| (Л У=1.........4),

я12 = 1 -f 4ет)3 cos 8; аи = 2 [5 + т)2 (iei) cos & Зт)

а22 = 4cetf sin 0; а24 = Qcetf sin 9; а3 2 = 2ет)3 sin 9; яз4 = 4ет|3 sin ft; д42 == — 2cer,3 cos

7)];

a43 — aU— C-

(24)

(25)

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

/</=1.....4).

(26)

где Д — определитель, полученный из матрицы К при 9 = Д/— определитель, полученный из Д заменой /-го столбца столбцом начальных условий (лго, х0, Уо,Уо)-

Аналогичным методом решается третье уравнение системы (1). В этом случае

(Я ~ <7о\ (27)

где

q — & + е sin 9 — ■

т (q — 9о\ , • « . Л

= z°ИГа cos\Г2-) + z^sln{

5 — if ‘ . « =

2

З + Yj3 4 : n=_l~

[t

e sin &

■ e2 I 1 + e cos 0 у і m0 = m (&o); n0 = n (&o); 9o = <7 (’%)•

== arctg xqrr tg-J-)]; f (28)

Изменение величины q($) в зависимости от эксцентриситета опорной орбиты показано на фиг. 1.

Анализ полученных решений и сравнение их с точным решением (фиг. 2) показывают, что асимптотический метод дает хорошую точность для е<0,1 во всем диапазоне изменения истинной аномалии (0-< ft360°). В частном случае прие = 0 приближенное решение совпадает с точным.

Решение методом рекуррентных соотношений. Перейдя к новой независимой переменной 0, вводим новые неизвестные

*1 =*; х2=у;

*3 ==

Xi = sin ft; хъ = cos 8;

1

L6

xR = -

7) (ft) ’

dx

w

dv

W dz Ж

^j0 = xcosft; xu =y cos ft; x12 = x sin ft; Xj3 = y sin ft;

Xr, = sin ft

x8 = sin ft x9 — sin ft

dx

Хц Ж ;

•*•15 = dy ~db’

dz

■*16 = db ’

(29)

Фиг. 2

которые дают возможность преобразовать исходную систему (1) к более удобной для последующих выкладок системе дифференциальных и алгебраических уравнений: dxi dxч___ dx3

dft

dxK

rfft

VJ4’ db x\ + x\ -

15-

= X

16>

1; Xq — 1 &x5 Xq~,

x

10 ■

Xi—XiXn, Xfj—ЛГ4ЛГ15, Xg Х4ХЩ,

~ X§ X}] Xji — X§ X.2> X-i2 “ X^ Xfj AJjg == X4 X2,

dx

14

db

dx

is

= lex^x-, -f ex6x10 + 2xis — 2ex6 xl3;

2ex§x$ Sx6xt -г- ex§ Xi — 2xj4-j-2exg x^2\ = 2ехЁх9 — x6x3.

dx i6

rfft

(32)

Представим неизвестные .*;*(&) (1-<г<16) в виде разложений по степеням 8:

00

*,(»)= (1<г<16), (31)

л-О

где х№ — неизвестные коэффициенты. С учетом (31) выпишем

рекуррентные формулы для определения неизвестных коэффициентов:

(п + 1) 441 = 4,4>; (я + 1) х%>+1 = х^);

(П + 1) 43|1 = х(п,6)’ (« + П *<4>1 = *<*>;

244) 44) = - £, *<4) Ч4АТ - х -ет;

7=0 7 = 0

Х(п] = - е 2 ХТ5’ Х<®-Т; 1' = X ХТ4> Хп-Г

т=0 7=0

^ = £*(«*£!); *(п9) = Х**4*™;

7=0 7=0

*(Г= £*?’*!,",; 4"’-

7=0 7 = 0

412) = 2 Х'4> Л«-Т; Л«3) = Е ХТ4) Л«-7;

7 = 0 7 = 0

<« +1)4'|\ = 2е^ х^ 471т + е 2 ХТ хп-1 + 2хп6) ~2еЕ хт6)

7=0 7 = 0 7=0

(п + 1) х(*\ =2е£ х«> х^ ■ 3 £ х[в) 42>7 +

*[•—0 7=0

+ в 2 х16) 4’Д - 2414> + 2еХ -^6) ^Д;

7=0 7 = 0

(Я + 1)^ = 2^- £*‘вЧ81т.

7 = 0 7=0 ;

где я> 0.

Определим условия, при которых разложения (31) сходятся. Допустим, что имеем неравенства

!*Ю|<а,/Ст (1С/С16); (33)

здесь

(34)

В,— постоянные, определяемые начальными значениями неизвест-

79

ных в момент £ = 0;Хиа(— постоянные, подлежащие определению. Предварительно учтем очевидные соотношения

7=2

ЕЧ*»-т

7=1

я + 2 ^ 3

я (л + 1) ^2

1

(я — 1 ~Ь 2^д_1) ___

(я-}-1)(я + 2) 2 (я—1)

И

п

2 (я — 1 -|- 25п-1) я (я + 1) (я + 2)

X":

-1 + 25л—1 ^ 5 1

Т+2-----<Т;П> '

(35)

где

п 1

*.-2т-

7=1

Анализируя с помощью (33)—(35) каждое уравнение системы (32), получим систему неравенств, дающую достаточные условия справедливости рекуррентных уравнений и включающую постоянные аь В1 (1 16) и X.

Если удовлетворяются неравенства (33) и все полученные неравенства, то выражения для неизвестных (31) сходятся при условии, что сходится сумма

2*7

7*1

■Е-

7 = 1

ХТ

1(7 + 1)

8Т1

1

Это условие реализуется при Л>0, ?>-1. Зная зна-

чения р и е, определив Вк (1 А < 16) из начальных условий согласно (29) и принимая одни коэффициенты а1 (1-<г-<16) достаточно малыми, другие — достаточно большими, чтобы удовлетворялась рассматриваемая система неравенств, найдем значение Х = Хтах, соответствующее всем неравенствам и определяющее интервал значений 8, на котором имеет место сходимость разложений (31). Значению Хтах отвечает минимальная величина независимой переменной 8, так как |8|<+/Хтах.

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

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

Выясним погрешности метода при ограничении несколькими первыми членами ряда.

Остаточный член ряда (31) имеет вид

/?5?(в)| =

«о «о

2 ■'?’8' <в'27Тпт)1

у=л+1 л + 1 III”1/

X7»7 .....л+1 £ 1

а7 к

т(т+ 1)

7=л+1 1 41 7

7 = л + 1

Т (т + 1)'

Но ііт V -т——гг = ——г . Отсюда можно заключить, что ряд ^т+1) п + Х

(31) сходится, так как

<37>

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

В качестве примера рассмотрен вариант, соответствующий следующим начальным условиям: лг0 = 100 «:лг, у0 = г0 = 50 км, Л'0 = — 200 м/сек, у0 — г0 = — 50 м/сек, 80 = 1,1 • м/сек, 80 = 0, е =0,05.

Результаты расчета, представленные на фиг. 4, показали, что максимальная ошибка определения координат по приближенным

Необходимо отметить, что приближенные формулы, получаемые из рекуррентных соотношений, имеют простой вид, хорошо аппроксимируют точное решение для любых значений е<С 1, а сходимость полученных разложений имеет место на достаточном для сближения интервале времени (см. фиг. 3).

Расчет корректирующих импульсов, характеристической скорости и времени сближения. При решении задачи импульсного корректирования траектории сближения КА удобно использовать метод, основанный на определении переходной матрицы, дающей связь между текущими и будущими значениями фазовых координат [2], [7]. Используя решение, полученное методом рекуррентных соотношений, найдем значение переходной матрицы.

6—Ученые записки К» 5

81

Представим уравнения (1) в форме

* ДО = О (*)•*(*)=• О (*)•

*Ц) І *(*)!

где б(£) — матрица коэффициентов системы, X (() — вектор-столбец

фазовых координат, х(Ь), х(() — векторы компонентов положения и скорости соответственно.

Решение (38) может быть представлено в виде

Х(*).

х(%) <3і (х, 0 <?2 (х, Ї) I II *(0

Х(х) = <2(х, І)-Х{І) = <3з(х, 0 <ЭЛХ> *) 1 1 х (Ґ)

(39)

где <3(т, ^ — переходная матрица.

Найдем алгоритм двухимпульсной коррекции при сближении в некоторый фиксированный момент встречи т. Предполагаем, что скорость в момент коррекции изменяется мгновенно: первая коррекция производится в момент t=t0 таким образом, чтобы ликвидировать рассогласование по дальности (при 1 — х = 0), вторая коррекция ликвидирует рассогласование по скорости в момент

сближения (при £ = х, л = 0, х = 0).

Тогда из (39) при t = tй следует, что

X со = Я1 (х, *„) х (*0) + <32 (х, 4) л (*0) = 0, откуда первый корректирующий импульс

где

К=Х*(*о)'

'О* і', І0)(3\(Х> *о) Х(*о)>

(40)

(41)

У0 и Ук— векторы скорости в момент / = (0 до и после проведения коррекции соответственно.

Второй корректирующий импульс находится также из уравнения (39) при условии £ = т:

где

Ун-

\Р2= У,-Ун,

■■ Хн (-С) = [<?4 (х) 5 (х, *0) + <23 (х)] X (*„), Б (X, І)= — <ЗҐ (х, і0) (І! (X, /0),

(42)

(43)

УИ и Ут — векторы скорости в момент t = x соответственно до и после проведения коррекции, причем Ух = 0.

Закон наведения полностью определится уравнениями (40)—(43), позволяющими рассчитать потребные скорости и корректирующие импульсы.

В решении (31) учтем первые четыре члена разложения. Тогда

где Е—единичная матрица 3X3,

1 О, Р О

Q2 (0 —

2

Gj t* t

о

Q<(*) =

1 G^t О

-Gxf 10 О 0 1

•*о + ~2~ @212 Н—G313 G2* + -i-G8f2

Я|(0 = Уо + ^412, + -Q-G513 ; Я2(0 = o4*+i-G5^

20-4-G6^2 — G61

G1=2ey?{b0);

G2 = с2 ет)3 (&0) (cos 90 x0 — 2 sin Ь0у0)\

G3 = 2с1 т] (»0) [2се sin 80 х0 — 2х0 + erf (»0) (3 + е cos &0)у0 — е sin j/0];

G4 = с2 rf (&0) [2 е sin &0 х0 + (3 4- е cos в0)л1;

<35 = 2 с2 7J (&0) [— с^3 (&0) cos &о х0 + e-q (&„) sin &0 л:0 + 2се sin &0_у0 — 2у0 ];

G6 — c2 т)3 (&0) z0.

С помощью формул (44), учитывая (41) и (43), можно получить значения импульсных коррекций скорости

(45)

где

Qi('t)=

/^2 9 Oi X-

<32(Ф

■Gb2

4 T — 2G, 0

2Gj 4 T 0

0 0 4 +G ix2 T

2 (2+G?t3) 2 Gx 0

T

- 2Gj 2(2+G? x2) 0

T

0 0 4 -|- Gi x2

Суммарные энергетические затраты, выраженные через характеристическую скорость, необходимые для сближения, полностью определятся заданным временем т и начальными значениями фазовых координат:

где V,- (2=1, 2, 3) —проекции вектора скорости на оси х, у и г соответственно.

В качестве примера рассмотрен случай сближения аппаратов для эллиптической опорной орбиты с перигеем Гл = 6870 км при Ро = 20 км, К0 = 20 и 100 м/сек, ф0 = 30°, в0 = 150°, где Ро и У0 — начальные значения дальности и скорости, ф0 и ®о ~ углы между линией визирования, вектором скорости и положительным направлением оси х соответственно.

Расчеты показали, что значения характеристической скорости существенно изменяются в зависимости от эксцентриситета опорной траектории и начальных значений истинной аномалии (фиг. 5 и 6).

Анализ данных показывает, что легко могут быть выделены области значений е и &0> при которых достигает минимума.

Найдем оптимальное время,

кую скорость, из условия равенства нулю частной производной Wl по т. Из уравнений (37) можно получить значения компонентов скорости (при < = т) и их частных производных по х, подстановка которых в указанное условие приводит к алгебраическому уравнению девятого порядка относительно *, не разрешимому в общем виде. Искомый корень этого уравнения можно приближенно представить в виде ___________________

т ^ 1 — 12 с2 тг)4 (80) хр 1

т»- 6 '

Ро

где т0 =—™ , р0 и р0 — начальное расстояние между аппаратами

Ро

и скорость сближения вдоль линии визирования. Можно показать, что соотношение (47) имеет место при выполнении условия

~\/~3

1Ро^ 6СТ|2(&0) ‘

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

Полученная формула (47) характеризует оптимальное время сближения в центральном гравитационном поле в зависимости от эксцентриситета опорной орбиты, начального положения аппарата на орбите и значения времени сближения х0, полученного для плоскопараллельного поля. Расчеты позволяют установить практическое совпадение хт с т0 для сравнительно малого времени сближения (до тт< 100 сек), когда центральностью гравитационного поля можно пренебречь.

ЛИТЕРАТУРА

1. Белецкий В. В., Егоров В. А. Межпланетные полеты с двигателями постоянной мощности. .Космические исследования', т. 2, вып. 3, 1964.

2. Ч а р н ы й В. И. Об изохронных производных, Сб. „Искусственные спутники Земли". М., АН СССР, № 16, 1963.

3. Е р м и л о в Ю. А. О расчете импульсов для сближения спутников в центральном ньютоновском поле тяготения. .Космические исследования", т. 6, вып. 6, 1968.

4. Р а п о п о р т И. М. Об устойчивости регулируемых процессов. ДАН СССР, т. 158, № 2, 1964.

5. А б г а р я н К. А. Асимптотическое расщепление уравнений регулируемого объекта и системы регулирования. ДАН СССР, т. 158, № 3, 1964.

6. Stef fens en J. F. On the differential equations of Hill in the theory motion of the moon. „Acta Mathematical, v. 93, No. 3—4, 1955 (I); v. 95, No. 1-2, 1956 (II).

7. S с о n z о P. An astronomical approach to the problem of satellite rendezvous. .Astronautica Acta", v. IX, No. 5—6, 1963.

Рукопись поступила 301XI 1970 г.