Исследование плоского движения относительно центра масс космического аппарата под действием гравитационного и аэродинамического моментов при снижении с круговых орбит Текст научной статьи по специальности «Физика»

УДК 629.78.015

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

© 2010 Л. В. Глухова, И. А. Тимбай

Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет)

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

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

Рассматривается плоское движение относительно центра масс осесимметричного космического аппарата (КА) при снижении с низких круговых орбит. Исследуются случаи, при которых в процессе снижения происходит изменение характера движения: вращательное движение переходит в колебательное, колебательное движение переходит во вращательное, колебательное движение “скачкообразно” переходит в колебательное движение с другими амплитудными характеристиками. В [1] методом фазовой плоскости исследованы случаи плоского углового движения КА на низких круговых орбитах под действием гравитационного момента и восстанавливающего аэродинамического момента, имеющего синусоидальную зависимость от угла атаки, что характерно для КА сферической формы или формы тонкого конуса, при этом изменением плотности атмосферы в процессе движения пренебрегается. В [2] рассмотрены переходные режимы движения КА относительно центра масс на начальном участке траектории спуска в атмо-

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

Плоское угловое движение осесимметричного КА при снижении с низких круговых орбит под действием гравитационного момента и восстанавливающего аэродинамического момента, имеющего вид бигармонического ряда по углу атаки, описывается уравнением с медленно меняющимися параметрами вида а + а( 1)$>\па+ (Ь( г) + с)8т2а = 0, (1)

где

а - угол атаки;

a(г), Ь(г) - отнесённые к поперечному моменту инерции КА коэффициенты разложения восстанавливающего аэродинамического момента в ряд Фурье;

£ - медленно меняющийся параметр, переменность которого связана с медленным изменением плотности атмосферы в процессе снижения;

с = 3(А - С)^2/(2А) - коэффициент, обусловленный действием гравитационного момента;

А , С - поперечный и продольный моменты инерции КА;

ю = л]^/ Я - угловая скорость движения центра масс КА на круговой орбите;

/и - гравитационная постоянная планеты; Я - расстояние от притягивающего центра до КА.

Коэффициенты а(г) и Ь(г) уравнения (1), если зависимость плотности атмосферы от высоты аппроксимировать экспонентой и пренебречь изменением скорости центра масс КА при снижении, могут быть представлены в виде [2]

а = а0 г, Ь = Ь0 г, (2)

а0 = ааБ1У 2р0 /(2А),

-.а сп/2

b0 = a SlV2p0 /(2A), z = ехр(-Л( H - H 0))

где

a00 и aa - постоянные коэффициенты;

S - характерная площадь; l - характерный размер КА;

V - скорость полёта;

H - высота полёта;

H0, р0 - высота полёта и плотность атмосферы в начальный момент времени t0 соответственно;

X - логарифмический градиент плотности по высоте.

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

Считая, что угол снижения КА, движущегося по почти круговой орбите, мал и меняется медленно, изменение высоты по времени можно определить по приближённой формуле [3]:

H =-

CxV 3S

mg

p,

(З)

где

g = gn (Rn / R )2,

gn - ускорение свободного падения на поверхности планеты,

R - радиус планеты,

n

Cx - коэффициент силы лобового сопротивления, m - масса КА, р - плотность атмосферы.

Учитывая, что p = p0exp(-X(H - H0)),

пренебрегая изменением скорости центра масс и изменением расстояния от притягивающего центра до КА, из (3), полагая t0 = 0 , получим следующую зависимость высоты полёта КА от времени

H = H0 + ln(1 -Xfit)l X, (4)

где в = CxV 3Spo l(mg).

Для выяснения общих свойств движения КА, описываемого системой (1), воспользуемся методом фазовой плоскости. Интеграл энергии системы в случае невозмущённого движения

(z = const, a = const, b = const)

имеет вид

2

а /2 - a cosa-(b + c)cos a = h.

(5)

Тип движения системы определяется соотношением величин а, Ь, с, к. Имеют место три вида фазовых портретов.

1. IЬ + с 1< 0.51 а I. Фазовый портрет аналогичен фазовому портрету колебательных систем маятникового типа и для случая а > 0 изображён на рис. 1а. При а < 0 фазовая картина сдвигается по оси а на величину п. Вращательному движению соответствует условие

к > I а I -(Ь + с), колебательному -к < I а I -(Ь + с). На границе перехода из

области вращения в область колебания к = I а I -(Ь + с).

2. (Ь + с) > 0.51 а I, (Ь + с) > 0. На фазовом портрете появляются дополнительные особые точки типа «седло», соответствующие значениям угла атаки

а* = ± агеео8(-0.5а /(Ь + с)) + 2пп (п = ±1, ± 2,...), и имеют место три области движения -вращательная А3 и две колебательные А1 и А2 , разделённые сепаратрисами (рис.1б). Вращательному движению соответствует условие к > а2 /(4(Ь + с)), колебательному относительно а= 0 или а = п - к < а 2/(4(Ь + с)). На границе перехода к = а2 /(4(Ь + с)).

3. IЬ + с 1> 0.51 а I, (Ь + с) < 0. Фазовый портрет для случая а > 0 изображён на рис.1 в (при а < 0 картина сдвинута по оси а на величину п). Здесь значениям угла атаки

а* = ± агеео8(-0.5а /(Ь + с)) + 2пп (п = ±1, ± 2,...) соответствуют особые точки типа «центр». На фазовом портрете располагаются четыре характерные области движения - вращательная и три колебательные. Вращательному движению соответствует условие к > I а I -(Ь + с), на границе перехода между областью вращения и областью колебания относительно а=0 -

к = I а I -(Ь + с). Колебательному движению относительно а = 0 соответствует условие

I а I -(Ь + с) > к >-1 а I -(Ь + с). Границе перехода между колебательной областью относительно положения равновесия по углу атаки а = 0 и одной из двух колебательных областей относительно положений равновесия а*=± агееов(-0.5а /(Ь + с)) соответствует условие к = -1 а I -(Ь + с).

а

^

^—

б) (Ь + с) > 0.51 а I, (Ь + с) > 0

а

в) IЬ + с !> 0.51 а I, (Ь + с) < 0

Рис. 1. Фазовые портреты

В возмущённом движении с изменением величин коэффициентов а и Ь , обусловленным ростом плотности атмосферы в процессе снижения, происходит эволюция фазовых траекторий, в результате которой они могут пересекать сепаратрисы, попадая в различные области фазового портрета. Это явление сопровождается качественным изменением характера движения: вращательное движение может переходить в колебательное, колебательное движение может переходить во вращательное, колебательное движение может “скачкообразно” переходить в колебательное же движение, но с другими амплитудными характеристиками. Следует отметить, что если не учитывать влияние гравитационного момента и рассматривать движение только под действием аэродинамического момента, то эволюция фазовых траекторий происходит в рамках одного из трёх указанных фазовых портретов (коэффициенты а и Ь изменяются пропорционально), при этом имеют место однократные проходы через сепаратрисы фазовой точки из внешних областей во внутренние [2]. При совместном действии

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

0.5 I a I „

тов -------. Поскольку с ростом плотно-

I b + c I

сти атмосферы коэффициенты a и b растут, а коэффициент c остаётся практически неизменным (полагаем c = const ), то указанное соотношение может меняться, и в момент выполнения равенства

0.5 I a I ,

-------= 1 происходит смена фазового

Ib+cI

портрета. В процессе снижения КА возможны следующие случаи:

- вид фазового портрета не меняется;

- происходит смена фазового портрета, изображённого на рис. 1а, на фазовый портрет, изображённый на рис. 1 б, или наоборот;

- происходит смена фазового портрета, изображённого на рис. 1а, на фазовый портрет, изображённый на рис. 1 в, или наоборот;

- происходит смена фазового портрета, изображённого на рис. 1б, на фазовый портрет, изображённый на рис. 1 а, затем происходит смена на фазовый портрет, изображённый на рис. 1 в, или наоборот.

Высота полёта, соответствующая моменту смены фазового портрета, изображённого на рис. 1 а, на фазовый портрет, изображённый на рис. 1 б, или наоборот, учитывая (2), определяется по формуле

H а^б = H 0 - ln(c /(0.5ao - bo))/ Л. (б)

Высота полёта, соответствующая моменту смены фазового портрета, изображённого на рис. 1 а, на фазовый портрет, изображённый на рис. 1 в, или наоборот, учитывая (2), определяется по формуле

H а^в = H 0 - ln(c /( °.5ao - b0))/Л. (7)

Для описания движения системы с медленно меняющимися параметрами (1) будем использовать интеграл действия, записанный в форме

amax

I = J a da, (8)

°min

где amin, amax - амплитудные значения угла атаки (при вращении (^тт = -п,

amax =п; при колебаниях amin = ^max),

величина а определяется из (5).

Амплитудные значения угла атаки определяются из решения уравнения

h = -a COS amax - (b + C) COS2 amax .

Для системы (1) равенство I = const справедливо для большинства начальных условий с точностью O(elne) на временах порядка Не [4], где £ - малый параметр, характеризующий скорость изменения параметра г . Исключительное множество начальных условий, для которых эта оценка несправедлива, имеют меру O(£n), где n > l - любое наперёд заданное число. Режимы движения, соответствующие данным начальным условиям (режимы зависания аппарата в окрестности неустойчивого равновесия) в статье не рассматриваются.

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

по начальным условиям.

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

Предполагается, что космическим аппаратом за время движения от t = 0 до границы перехода совершено несколько оборотов или колебаний.

ЗЗ

Поскольку после перехода через сепаратрису фазовая точка может попадать в различные колебательные области, возникает задача выбора области продолжения движения. Внутренние области движения А1 и А2 отделяются от внешней

области А3 сепаратрисами 11 и 12, соответственно (рис. 1б, 1в). Вероятность Рі, і = 1,2 захвата в одну из колебательных областей А1 или А2 вычисляется по формулам [4]:

Р1_ 011

P2 ®l2

rdF

0 =-firf^dt (i = 1,2,3),

i і dz

F = H (a ,a, z) - H (0,a*, z),

(9)

(10)

где

H = 2 / 2 - a cos a — (b + c)cos2 a - гамильтониан;

fz = г ;

a = 0, a = a* - координаты седловой особой точки на фазовом портрете (E = 0 в седловой особой точке и на сепаратрисах, E > 0 в A3, E < 0 в A12).

Интегралы (10) вычисляются вдоль сепаратрис l1 , l2 и l3 , параметризованных временем t невозмущённого движения по ним. Физический смысл величины 0{ -

скорость приближения площади, ометае-мой фазовой траекторией, к площади, ограниченной сепаратрисой. При 0Ь > 0

фазовая траектория «погружается» из внешней области во внутреннюю, при 0Ь < 0 фазовая траектория «выталкивается» из внутренней области во внешнюю. Определить функции 0Ь можно, используя также выражения для интеграла действия, вычисленного вдоль сепаратрис, по формуле [4]

dIi

0, = f—L .

li z dz

(11)

рис, получим аналогично [2], учитывая при этом обусловленный действием гравитационного момента коэффициент с.

Для фазовых портретов, изображённых на рис. 1а и 1в:

= 4Л/2(Ь0 г» + с) [^ 1 + )],

при (b0 z* + c) > 0,

(12)

I3 — 4^/_2(b0 z» + c) [л]и, +1 + u»ln((1+>/u» +1)/ )]

при (b0 z* + c) < 0,

(13)

Ii1 = Ii2 = U -2(b0z*+c) (%A - u - ut in((1+V1 - u,) /)) .,(1 4)

где

u, =10.5a0z* /(b0z* + c) I, z* - значение параметра z в момент перехода через сепаратрису.

Для фазового портрета, изображённого на рис. 1б:

Ik = 2j2(b0 z* + c) (sin a, - a,cos a,), (15)

I2 = 2^2(b0 z* + c) (sin a, + (п-а, )cos a,) (16)

где a* = arccos(-0.5a0 z* / (b0 z* + c)).

Значение интеграла действия I0 определяется по начальным условиям и может быть вычислено по формуле (8), а также по аналитическим выражениям для интеграла действия, приведённым в [2].

Учитывая постоянство интеграла действия, из решения трансцендентного уравнения Ih (z*) = I0 можно определить

значение параметра z* в момент, соответствующий переходу через сепаратрису. Тогда, учитывая, что z* = exp(-l(H, - H0)),

высота H* , соответствующая моменту перехода из одной области движения в другую, определяется по формуле

H, = H0 -ln(z,)/A. (17)

Учитывая (11), дифференцируя выражения для интеграла действия Il по параметру z , получим следующие аналитические выражения для функций 0L:

для фазовых портретов, изображённых на рис. 1 а и 1 в:

Аналитические выражения для интеграла действия , взятого вдоль сепарат-

( 4Ь

[V и,-1 + и.агеІ£^1/(и,-1))] +

-\]2(Ь0 г* + с)

ч + 4^2^г* + с)каг^^и, -1))

при (Ь0 г* + с) > 0, (18)

( -4Ь

0 = І.г

I [д/Й +1 + М,1п((1^и, + 1)/^/и7)] -

^ 2(Ь0 г* + с)

- 44-2(Ь0 г* + с)у,1п((1^л/иТП)^7и7) при (Ь0 г* + с) < 0, (19)

(20)

-2Ь,

УТ-ЙГ - и* 1п((1 + ^1-и* )/^ЙТ)] +

0, =0, = / л]~2Фог* + с)[

11 12 ■> I 4

^ + 2]-2(Ьо г* + с)V* 1п((1^ 1- и* )/^/й*)

где V» = 0.51 а01 с /(Ь0г» + с)2;

для фазового портрета, изображённого на рис. 1б:

2Ь

г^ІГ

8іп а* -а* сов а* 1 +

л/2(^0 г* + с)

+ 2^2(Ь0 г* + с)а*г* 8Іп а

Л

(21)

2Ь„

г^іпа* + (п-а*)со8 а* ]-

^2(Ь0 г* + с)

- 2^2(Ь0г* + с) (п - а* )у* віп а*

(22)

0.5а0с0

0.25(а0 г*)2

1 -

(Ь0 г* + с)

2 (Ь0 г* + с)2

На рис.2 показан характер изменения плоского движения КА для следую-

щих начальных данных:

-7-2

а0 = 1.6-10-7 с-2, Ь0 = 5.8-10-7 с-2,

с = -1-10 с-, а0 = 0.3 рад, а = 6.9 -10"4 рад/с, Я0 = 300000 м,

Л = 1/ 43000, в = 0.06924 м/с. Поскольку в данном случае

I Ь + с 1> 0.51 а{) I при (Ь + с) < 0

и

I а01 -(Ь0 + с) > ^0 > -1 а01 -(Ь0 + с), то КА в начальный момент совершает колебания относительно а = 0 в рамках фазового портрета, изображённого на рис. 1в. Высоты полёта На^в и На^б, на которых происходит смена фазовых портретов, с учётом (6) и (7) равны 282132 м и 270194 м, соответственно. Таким образом, движение до высоты На^в происходит в рамках фазового портрета, изображённого

на рис. 1в, затем до высоты На^ в рамках фазового портрета, изображённого на рис. 1а, и продолжается в рамках фазового портрета, изображённого на рис. 1 б. Вычисление высот Н*, соответствующих

моментам перехода из одной области движения в другую, в соответствии с (17) даёт следующий результат. При движении в рамках фазового портрета, изображённого на рис. 1в, КА совершает только колебания. При движении в рамках фазового портрета, изображённого на рис. 1а, КА на высоте Н* = 275050 м переходит во вращательное движение. При движении в рамках фазового портрета, изображённого на рис. 1б, КА на высоте Н*2 = 261351 м вновь переходит в колебательное движение. В соответствии с (9), учитывая (21), (22), вероятность попадания в область колебания относительно а = 0 равна 34%, вероятность попадания в область колебания относительно а = п равна 66%. Значения высот перехода, полученные по аналитическим выражениям, с точностью до периода колебаний совпадают с результатами численного интегрирования уравнения (1). Значения вероятностей, полученные путём многократного численного интегрирования дифференциального уравнения (1) при незначительном варьировании начальных условий (значения а задавались равномерно из промежутка [6.8-10-4;7.3-10-4]) соответствуют значениям вероятностей, вычисленных по аналитическим выражениям.

а) зависимость угла атаки от высоты

следования. - 1997. - Т.35. - №3. - С. 279286.

3. Балк, М. Б. Элементы динамики космического полёта [Текст] / М. Б. Балк.

- М.: Наука, 1965.

4. Нейштадт, А. И. Об изменении адиабатического инварианта при переходе через сепаратрису [Текст] / А. И. Ней-штадт // Физика плазмы. - 1986. - Т.12. -Вып.8. - С. 992-1001.

References

1. Beletsky, V. V. Artificial satellite motion around its centre of mass [Text] / V. V. Beletsky. - Moscow: Nauka, 1965.

2. Aslanov, V. S. Transient modes of spacecraft angular motion at the upper section of the reentry trajectory [Text] / V. S. Aslanov, I. A. Timbay // Kosmitcheskiye Issledovanya (Space Research). - 1997. -Vol.35. - No.3. - PP. 279-286.

3. Balk, M. B. Elements of space flight dynamics [Text] / M. B. Balk. - Moscow: Nauka, 1965.

4. Neuschtadt, A. I. Variation of adiabatic invariant during the crossing of a sepa-ratrix [Text] / A. I. Neuschtadt // Physica Plasmy (Plasma Physics). - 1986. - Vol.12.

- Issue8. - PP. 992-1001.

RESEARCH OF PLANAR MOTION OF SPACECRAFT AROUND ITS CENTER OF MASS UNDER THE INFLUENCE OF THE GRAVITATIONAL AND THE AERODYNAMIC MOMENTS DURING ITS DESCENT FROM CIRCULAR ORBITS

© 2010 L. V. Glukhova, I. A. Timbay

Samara State Aerospace University named after academician S. P. Korolyov

(National Research University)

The paper deals with the planar motion of an uncontrolled axisymmetric spacecraft around its center of mass under the influence of the gravitational moment and the restoring aerodynamic moment. The restoring aerodynamic moment of the spacecraft is described by an odd Fourier series in the angle of attack with the two first harmonics. The evolution of phase trajectories is studied on the basis of analyzing the action integral. The evolution of the phase trajectories is caused by the descent of the spacecraft from a low circular orbit due to the atmospheric deceleration. The moments of transition between various areas of the phase plane are determined. For the cases of motion, when, intersecting the separatrix, the phase point may fall into various oscillation regions, the formulas for determining the possibility of capture into any region are found.

Planar motion, gravitational moment, aerodynamic moment, phase portrait, action integral, probability.

б) фазовая траектория

Рис. 2. Характер изменения углового движения КА

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

Библиографический список

1. Белецкий, В. В. Движение искусственного спутника относительно центра масс [Текст] / В. В. Белецкий. - М.: Наука, 1965.

2. Асланов, В. С. Переходные режимы углового движения КА на верхнем участке траектории спуска [Текст] / В. С. Асланов, И. А. Тимбай // Космические ис-

Информация об авторах

Тимбай Иван Александрович, д.т.н., профессор, заведующий кафедрой высшей математики. Самарский аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет). С. П. Королёва. Область научных интересов: динамика летательных аппаратов. E-mail: timbai@ssau.ru.

Глухова Лия Валерьевна, аспирант кафедры высшей математики. Самарский аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет). Область научных интересов: динамика летательных аппаратов. E-mail: bazonchik @ gmail.com.

Timbay Ivan Alexandrovitch, doctor of technical sciences, professor, head of the department of higher mathematics. Samara State Aerospace University named after academician

S. P. Korolyov (National Research University). Area of research: aircraft dynamics. E-mail: timbai@ssau.ru.

Glukhova Liya Valeryevna, postgraduate student of higher mathematics department. Samara State Aerospace University named after academician S. P. Korolyov (National Research University). Area of research: aircraft dynamics. E-mail: bazonchik @ gmail.com.