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

УДК 629.78

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

© 2003 Е.А. Кеньшов, И.А. Тимбай Самарский государственный аэрокосмический университет

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

Постановка задачи

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

указанной асимметрии может быть представлен в виде четного ряда Фурье по углу атаки с постоянным членом (средним значением момента) и первым косинусоидальным членом. Тогда, добавляя косинусоидальную компоненту к основному восстанавливающему аэродинамическому моменту, описываемому нечетным рядом Фурье с двумя первыми гармониками [2], считая среднее значение момента асимметрии и значение демпфирующего момента малыми величинами, учитывая, что коэффициенты моментов медленно меняются во времени из-за медленного изменения плотности атмосферы в процессе спуска, движение КА можно описать следующим уравнением [2, 4]

а + a(z)sina + b( z) sin 2a +

+ c( z )cosa = e<7( z)a + £%( z ^ Q)

где a - угол атаки; a(z), b(z), c(z) - коэф-

фициенты моментной характеристики; a(z)

и ^(z) - усредненные по углу атаки коэффи -циенты демпфирования и асимметрии; z -медленно меняющийся параметр, переменность которого связана с медленным изменением плотности атмосферы в процессе спус-

ка; е - малый параметр, характеризующий малость коэффициентов демпфирования и асимметрии.

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

a(z) = a0 z, b(z) = Ь0 z , c(z) = с0 z, ъ(,) = о 0z, ^(,) = £ 0z, а„ = -т£ш0;р(^)/(21),

Ь„ =-mtSlv02p(to)/(2/), с0 =-АуС^П'02р(10)/(2/),

= (С + Дт, )!рЛ )/(),

<3, = (КI /1 - С“ 0 / m)SV„p(t0) / 2,

z = exp(^(t -10)), e = rk¥o\sin 0o |,

где ma,mb - коэффициенты основного восстанавливающего аэродинамического момента; m“ - усредненный по углу атаки коэффи -циент аэродинамического момента демпфи-

С a

>a0 - усредненная производная по углу атаки коэффициента аэродинамической подъемной силы; Ст0, Ст1 - первые два члена разложения коэффициента тангенциальной аэродинамической силы в четный ряд Фурье;

Amz - коэффициент аэродинамического момента асимметрии; Ay = Ay/l - безразмерное смещение центра масс КА от оси его геометрической симметрии; l - характерный размер; S - характерная площадь; / - поперечный момент инерции; m - масса аппарата; V0 - скорость полета; 0О - угол наклона траектории; p(to) - плотность атмосферы в начальный момент времени; X - логарифмический градиент плотности атмосферы по высоте.

Невозмущенное движение КА

В случае невозмущенного движения ( £ = 0, z = const) уравнение движения КА запишется в виде

а + а£ша + bsin 2а + ссо& а = 0. (2)

Следует отметить, что уравнение движения физического маятника по углу его отклонения от вертикали, когда он установлен на вращающейся платформе и точка его подвеса смещена относительно оси вращения платформы [5], имеет вид аналогичный (2).

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

й2

—— асо^а-Ьсо52 а + с&іпа = Ъ. (3)

Фазовым пространством рассматриваемой системы является цилиндр с координатами а, а . Уравнение (3) определяет связь

между а и а, тем самым, является уравнением фазовых траекторий. Из (3) следует, что

a = ±^2( f( a ) + h)

где

f( a) = acos a + bcos2 a-csin a. (4)

Экстремальные значения функции f(a) соответствуют состояниям равновесия уравнения (2), т.е. особым точкам на фазовом цилиндре a, a.

Из выражения (4) функции f(a) следует, что в зависимости от значений параметров a, b, c могут существовать или две, или четыре особые точки. Бифуркационное соот-

ab

ношение параметров Л = ^ и z ^, разделяющее эти случаи, находится из условия слияния двух особых точек, что осуществляется при одновременном выполнении двух условий:

f'(a) = 0 и f"(a) = 0 .

Разрешая эти соотношения относительно П и Z , получаем параметрическое представление бифуркационной кривой П = П(Z) : П = ctg3a, Z = -(2sin3 a)~l. (5)

Согласно (5) на плоскости п, £ протекают четыре ветви кривой п = П( С), соответствующие изменению а в интервале [- п, п]. Данная кривая, изображенная на рис.1, разбивает плоскость на три области. При параметрах П, С , находящихся внутри области, ограниченной бифуркационной кривой, на фазовом цилиндре имеются две особые точки, а в других областях - четыре особые точки. Для параметров, принимающих значения £ = ±0,5, п = 0 происходит слияние трех особых точек: при а = 0 для £ = +0,5 , п = 0 и

при а = ±п для £ = -0,5, п = 0.

На рис.2 показан фазовый портрет системы, когда имеются четыре особые точки, на рис.3 - две особые точки.

Как видно, космический аппарат в зависимости от начальных условий (И) и величин а, Ь, с может совершать как вращательные, так и колебательные движения, области которых разделены сепаратрисами.

Найдем решение для угла атаки в невозмущенном вращательном движении. Вводя

= а

новую переменную и = ^ —, уравнение (3) можно записать в виде

и2 = — [И - а + Ь ] 4 -—си 3 + — [И - Ь ] 2 -

а

-—си + — \к + а + ЬІ = О(и) 2 4

. (6)

Разделяя переменные в уравнении (6), получим

а

Рис.2. Фазовый портрет: а = 3, Ь = 4, с = 1 (п = 3, £ = 4)

йі = ±-

йи

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

і — і0

и

=1-

йи

(7)

,4см ■

Интеграл (7) относится к классу эллиптических интегралов и приводится к нормальному эллиптическому интегралу Лежандра первого рода [6]. Результат интегрирования зависит от типа корней полинома четвертой степени 0(и). В случае вращения все четыре корня полинома комплексные. Введем следующее правило нумерации этих корней:

Ь— ± /с1, Ь2 ± гс2, Ь— > Ь2, с > 0, с2 > 0.

Приведение интеграла (7) к нормальному эллиптическому интегралу первого рода осуществляется посредством преобразования

П

52 особы е точки

3-

4 особые точки 1 4 особые точки

-4 -2^>’0_ -32 особы \ 2 е точки 4

Рис.1. Бифуркационная кривая

Рис.3. Фазовый портрет: а = — 3, Ь = 2, с = —1 (п = 3, С = —2)

и = и(ф), отображающего интервал интегрирования [0, и] в соответствующий интервал действительного аргумента [ф0, ф]. Вид преобразования и = и(ф ) зависит от типа и сочетания корней, а также от знака старшего коэффициента полинома 0(и) [6]. Для рассматриваемого случая вращения все четыре корня полинома комплексные и преобразование может быть представлено в виде

u = b1 + c1tg (ф + ф*),

где ф, -

03 +04

034 = arctg

(В)

c1 ± c2

b1 - b2

- острые углы.

2 ’ ~3,4

Полагая, что в момент времени і = і0, а = а 0 , при этом

a0

#-2-- b1

ф0 = arctg-----------------2-ф,

(9)

Тогда, учитывая (8), интеграл (7) преобразуется к следующему виду:

*- *0 = —г[^( ф,к) - ¥( Ф0,к)], (10)

где

F( ф,k) = Ь=

o -у1 -

dф

- эллиптический

k sin ф

эллиптического интеграла,

8 =

h - a + b|

2yj cos 05

Vc1cT, =

0? cos 03 1

_ ___________0

Icos05 2 , 4 cos04 , 2 5

а = 2аrcfg[b1 + с^(ат((* - *0) + т0 ,к)+ф„ )]. (12)

Следует отметить, что формула (12) задает непрерывную функцию на отрезке

ае [-п; п]. В рассматриваемом случае вращения необходимо продолжить функцию а на соответствующие отрезки, исходя из ее непрерывности.

Возмущенное движение КА

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

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

интеграл первого рода, k2 = sin2 05 - модуль

a0 = -0,02 с-

b0 =-0,005 с"2,

острый угол.

Обращая эллиптический интеграл, из равенства (10) получим

ф = ат[(* - *0) + т0 ,к ], (11)

где т0 = Г(ф0,к) .

Учитывая формулы (8), (11) и то, что а

и = ^ —, имеем следующее решение для угла атаки во вращательном движении:

с0 = —0,005 с"2, = 0,0008 с"2,

а0 = —0,00001 с_1, а0 = 10 град,

а0 = 30 град / с, в= 0,05 с _1).

Будем полагать, что критерий применимости асимптотических методов в задачах спуска КА в атмосфере выполняется. Для плоского движения он имеет вид [4]

XV0 \sin 00

2a 0 в

>> 1. (13)

Тогда для описания вращательного движения системы с медленно меняющимися параметрами (1) воспользуемся законом изменения интеграла действия [4]

c

2

а, рад/с

Рис.4. Переход вращательного движения в колебательное

где

Ig = С0+ 2п^0г ■

= | Ша ,

(14)

(15)

а а определяется из интеграла энергии (3).

Исключая из уравнения (14) время, принимая *0 = 0 и учитывая, что в этом случае , = ехр(в*), получим

С0^ , 2п^0

+

ё, в в После интегрирования имеем

I (,) = А - Ве-С(г-1;.

(16)

Здесь

- 2п^0 С0

В=-

2п^ 0

+1

С0

g 0

С = -^. в ;

ствия в процессе спуска изменяется. В момент времени, когда кривая I (г) пересекает

кривую интеграла действия, взятого на сепаратрисе, разделяющей вращательную и колебательную области, произойдет переход вращательного движения в колебательное.

Вычислим интеграл действия вдоль сепаратрисы. На сепаратрисе многочлен четвертой степени 0(и) имеет два комплексных корня и два совпадающих действительных ( и34 ). Этот кратный корень должен одновременно удовлетворять условиям 0(и34 ) = 0, с>34 ) = 0. Этим условиям также удовлетворяет величина полной энергии И. В результате преобразований интеграла (15), ис-

= а

пользуя подстановку и = ^ —, имеем следующее аналитическое выражение для интеграла действия, вычисленного вдоль сепаратрисы:

/;= , (18)

где

В = |ё1 /1 + ё21 ^ + 13 , (19)

/1 =

г1п

VI. 2 Г~2 2

У + г -д/г -р

у]X + 2вт +ут2у/г2 -р2 7у2 -р2

22 при р < г ,

ц/х + 20т + 2 д/р2 -;

22 при р > г .

.л/Т^

п

----аг^- __________

2 2 2 2 д/р - г

^ 0 = I &(* 0), - значения интеграла действия

в начальный момент времени.

В случае, когда на границе атмосферы коэффициенты а, Ь, с существенно малы по сравнению с угловой скоростью вращения а0 , начальное значение интеграла действия определяется по формуле

^ 0 = 2па 0. (17)

Согласно формуле (16) интеграл дей-

arctg

■<]х + 20т + т2 р^г2 - р2 р^у2 + г

■\]х + 2вт + ут2 р^р2 - г2

I 2 2 I 2 2

1п(г)- 1пУ^р -г + ^У + г

-лр2( и324 + 1 х + 20и34 +

(и324 + 1 )

■7и324

2

2

2

при р < г ,

2

при р > г

3

т, п ■

У

Х±4 ( У-Х )2 + 402

п2 +1 т2 +1

20

Х + 20п + уп2

Х + 20т + ут2

У

п - и

34

и34 - т

42( т - п )

1 т2 +1

(0 - хи34 + тт -0и34 т),

42( т - п )

2 т2 +1

(0-Хи34 +уп -0и34 п),

у = И* - а0 + Ь0, х = 3ти34 - 4с0и34 + 2И* - 2Ь0,

0 = Уи34 - с0 .

Приравнивая выражение для интеграла действия, описывающего его изменения в процессе спуска (16), выражению интеграла действия, вычисленного вдоль сепаратрисы (18), получим следующее уравнение для определения значения параметра г , при котором происходит переход вращательного движения в колебательное

А - Ве ~С (г-1)= . (20)

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

2п^ 0(г -1) . .

1— • (21)

Тогда приравнивая выражение (21) к выражению для интеграла, вычисленного на сепаратрисе (18), и разрешая полученное уравнение относительно параметра 2, получим

г =

1.0 -

2п^0

4п^0

можно определить значение времени перехода по формуле

* = 1пг/в . (23)

В случае, когда кривые (16) и (18) не пересекаются, переход вращательного движения в колебательное не происходит, возникает длительное вращение КА (плоская авторотация). На рис.5 показан случай такого движения (начальные условия движения:

а0 =-0,02 с“2, Ь0 =-0,005 с“2, с0 =-0,005 с"2, = 0,00087 с"2,

с0 =-0,00001 с_1, а0 = 10 град, а 0 = 30 град / с, в = 0,05 с _1).

Критическую угловую скорость (а 0 )кр,

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

(а0)кр = [а-{а-Р/(2С)}ерг/(4ВС]/(2п), (24)

где Р = АС -VА2С2 -2СВ2 .

В случае, если демпфированием можно пренебречь, выражение для критической угловой скорости принимает вид

(а „)„ = В !в / (16п2 ^). (25)

а, рад/с

(22)

Зная величину г, при которой происходит переход вращательного движения в колебательное, и, учитывая, что г = ехр(в*),

Рис.5. Плоская авторотация

2

2

г

2

в

в

Следует отметить, что в отличие от аналогичных формул для критической угловой скорости, приведенных в [3] в интегральном виде, в формулах (24) и (25) параметр в выражен через элементарные функции (19). Формулы (24), (25) дают практически точный результат, когда справедливо уравнение для изменения интеграла действия (14), что соответствует большим значениям Д (условие (13)).

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Кузмак Г.Е. Динамика неуправляемого движения летательных аппаратов при входе в атмосферу. М.: Наука, 1970.

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

3. Ярошевский В.А. Исследование условий возникновения плоской авторотации кос -мического аппарата на этапе входа в атмосферу // Космические исследования. 2001. Т.39. №1.

4. Ярошевский В.А. Движение неуправляемого тела в атмосфере. М.: Машиностроение, 1978.

5. Бутенин Н.В., Неймарк Ю.И., Фуфаев Н. А. Введение в теорию нелинейных колебаний. М.: Наука, 1987.

STUDY OF TRANSITION OF SPACECRAFT WITH ASYMMETRY FLAT ROTATION TO FLUCTUATING MOTION WHILE ENTERING THE ATMOSPHERE

© 2003 E.A. Kenshov, I.A. Timbai

Samara State Aerospace University

Motion of Spacecraft with little asymmetry, restoring aerodynamic moment of which is described by Fourier series in the angle of attack with the first two sine and the first one-cosine terms, is considered in relation to its center of mass. Transition of Spacecraft flat rotation to fluctuation, conditioned by slow alteration of the coefficients of moment characteristics as well as by presence of little asymmetry and damping coefficients and their slow alteration, is studied. Analytical expression for the action integral, taken along the separatrixes that divide rotation and fluctuation regions of the system phase portrait, has been obtained. Analytical formulas for determination of the times of transition of rotation to fluctuation, and for critical angular velocity of the extra-atmosphere rotation, which causes prolonged body rotation (flat auto-rotation) when exceeded, have been derived.