УДК 629.78
ХАОТИЧЕСКИЙ АНАЛИЗ ДВИЖЕНИЯ КОСМИЧЕСКОГО АППАРАТА С МАЛОЙ МАССОВОЙ АСИММЕТРИЕЙ ПРИ СПУСКЕ В АТМОСФЕРЕ
© 2009 В.С. Асланов, А.С. Ледков Самарский государственный аэрокосмический университет Поступила в редакцию 15.01.2009
Методами хаотической динамики исследовано неуправляемое движение космического аппарата с малой массовой асимметрией при спуске в атмосфере. На космический аппарат действует восстанавливаемый момент, описываемый бигармонической зависимостью. С помощью метода Мельникова исследовано поведение системы в окрестности сепаратрисы. Получены аналитические решения движения космического аппарат по гомоклиническим траекториям. Найдены критерии возникновения хаоса и представлены результаты численного моделирования, подтверждающие справедливость полученных решений. Исследовано изменение толщины хаотического слоя в процессе спуска космического аппарата в атмосфере. Даны рекомендации по выбору уровня демпфирования, устраняющего возможность появления хаотических переходов на всей траектории спуска. Ключевые слова: космический аппарат, бигармонический момент, хаос, метод Мельникова, гомокли-ническая траектория, аналитическое решение, хаотический слой, устойчивость, сечения Пуанкаре.
ПОСТАНОВКА ЗАДАЧИ
В настоящее время в большинстве космических программ, включающих этапы доставки груза с орбиты на поверхность планеты, используются затупленные сегментально-коничесие космические аппараты (КА) малого удлинения (Beagle2, Huygens, Phoenix, проект NASA "Mars Sample Return") [1-4]. Такая форма позволяет обеспечить эффективное торможение в разряженных слоях атмосферы и обладает особенностью аэродинамической природы: аппараты такой формы могут иметь три балансировочных положения равновесия по пространственному углу атаки а (рис. 1). Учесть эту особенность позволяет аппроксимация коэффициента восстанавливающего момента би-гармонической зависимостью:
■0-2
0.41-
О -тгЧ тг/2 3it'4 гт
а
Рис. 1. Зависимость коэффициента восстанавливающего момента от времени
Асланов Владимир Степанович, доктор технических наук, профессор, заведующий кафедрой теоретической механики. E-mail: [email protected]. Ледков Александр Сергеевич, аспирант. E-mail: [email protected].
та= а Бта + Ь вт2а. (1)
Поскольку устойчивыми будут два положения равновесия а = 0 и а = ж, в процессе спуска КА его пространственный угол атаки может принимать как малые, так и большие значения. Это делает невозможным применение метода линеаризации. Для исследования движения КА необходим анализ нелинейных дифференциальных уравнений. В частном случае движения осесимметричного КА, когда состояние системы определяется изменением одной быстрой переменной - пространственного угла атаки, получены приближенные аналитические решения нелинейной возмущенной системы[5]. В случае движения КА с малой асимметрией система уравнений становится двухчастотной. При этом создаются предпосылки для возникновения различных видов параметрических резонансов. Нелинейный анализ резонансов строится, как правило, на использовании аналитических решений невозмущенного движения[6]. Эти решения имеют весьма сложный вид, что затрудняет их применение.
В рамках данной работы для анализа нелинейного движения КА применяются методы хаотической динамики, в частности метод Мельникова [7-10]. При хаотическом поведении системы в окрестности сепаратрисы могут происходить переходы из одной области движения в другую. Эти переходы можно рассматривать как проявления резонанса. Для иллюстрации эффективности методов хаотической динамики при решении задачи, рассмотрим один из видов асимметрии -смещение центра масс КА с оси симметрии формы на расстояние 1Т (рис. 2). В этом случае состояние системы определяется двумя быстрыми переменными: пространственным углом атаки а и углом аэродинамического крена р.
Рис. 2. Космический аппарат НЕЛИНЕЙНАЯ СИСТЕМА
Запишем систему дифференциальных уравнений, описывающую движение КА с малой асимметрией (Ф 0 ) относительно его центра масс О при спуске в атмосфере [6]
(R - G cos a)(G - R cosa)
a-тЦ-1 - Ma =
sin a
= s(Pa(a,q), z),
I _ G-Rcosa ф =— R-■
Ix
sin2 a
z = sФz (a, ф, z);
, „ qSa qSL . sФa = D—— zTcTsin ф,
(2)
(3)
(4)
V
Ma= qfma, D =
I
<;L2
zn
I
ya
m„
Здесь R и G - с точностью до множителя проекции вектора кинетического момента на продольную ось и на направление скорости, V - скорость движения центра масс КА, Б - коэффициент демпфирования, q = р(Н V2 /2 - скоростной напор, р(Н) - плотность атмосферы, н - высота полета, S - площадь миделевого сечения КА, Ь - характерный размер КА, 1х, I - продольный и по-
перечный моменты инерции КА, mzn
m* =■
dmz да.
компоненты матрицы производных аэродинамических коэффициентов по безразмерным угловым скоростям, заданных в системе координат, связанной с пространственным углом атаки, агп = агпЬ / V, агп - проекция угловой скорости на ось г , та- коэффициент восстанавливающе-
а
го аэродинамического момента, суа - производная коэффициента аэродинамической подъемной
силы cya (a) по углу атаки, mKA - масса КА, z = {R, G, V, H, в} - вектор медленных переменных, в - угол наклона траектории.
Коэффициент cT разложим в ряд Фурье и ограничимся рассмотрением нескольких первых членов
cT= c0 + c1 cos a + c2 cos 2a . (5)
Будем также считать, что коэффициент D постоянен.
При s = 0 система (1)-(3) приводится к невозмущенной:
(R - G cos a)(G - R cosa)
a-г4-L-ма = о. (6)
sin a
Уравнение (6) имеет аналитическое решение, форма которого зависит вида от корней полинома
f (u) = 2(1 - u 2)(h -
aSL
I
■qu
bSL
qu2) +
+ 2GRu - G2 - R2
(7)
где h - энергия системы. В случае, когда все корни полинома f (u) действительны ( щ , u2, u3, u4 ) и пронумерованы так, что щ > u2, u3 > u4, а u3 < u2 или u4 > u,, решение имеет вид [6]:
u, (u2 - u3) + u3 (u, - u2) cos2 Y
cosa = ^-(-) + ( ) 2-, (8)
(2 - u3) + ( - u2) cos y
где y = am(Pt + т0, k),
fi=J-
qSL(u1 - u3 )(u2 - u4 )
2I
k =
(u1 - u2 )(u3 - u4 )
, w \ , T0 - постоянная интег-« (u, -u3)(u2 -u4) 0
рирования. В случае, когда полином f (u) имеет как действительные u,, u2 (u, > u2 ), так и мнимые u34 ± iv корни,
(u2 + щ£)-(u2 -cosy C0Sa= (1 + ^)-(l-^)cosY , (9)
где a, a =л1 (ui- u34)2,
e=J-
qSLa1a2
k=
2I
1 (u -u^Xu -u34)+V
2aya1
Здесь и всюду далее i = 1,2.
УРАВНЕНИЯ ГОМОКЛИНИЧЕСКИХ ТРАЕКТОРИЙ
Для использовании метода Мельникова [7-10] необходимо найти гомоклинические траекторий а°г- (}) невозмущенной системы (6) (рис. 3).
n
Рис. 3. Фазовый портрет невозмущенной системы
Движению по сепаратрисе соответствует энергия, вычисленная в седловой точке a0
(G- R coso* )2 qSL, , 2 л h = ---—т^—I--(acosa0 + bcos а0).
_ . 2 * 2sin а
I
где
P = .
qSL(v0 - V2)(Vo - Vj)
2I
cn - эллипти-
ческий косинус [11]; получим уравнения гомок-линических траекторий в виде
cosa°(t) = v. +
v0 - v.
1 + ch-2(et + 70) , (12)
v - V,
00 =
f (cos o0(t))
(13)
1 - cos2 o0(t) .
Здесь траектории а01 соответствует значение j = 2, а траектории a02 - j = 1.
ОЦЕНКА УСТОЙЧИВОСТИ
ОБЛАСТЕЙ ДВИЖЕНИЯ
Переменные z медленно меняются в процессе спуска КА. Однако они не являются периодическими по времени величинами, поэтому применение описанного в [9, 10] модифицированного метода Мельникова для систем с медленными параметрами невозможно.
Рассмотрим возмущенную систему, получаемую из (2)-(4) при фиксации некоторого момента времени.
л =
Р = fp +£gp, а = L+sgo;
(R - G cosa)(G - R cosa)
I sin3 a - qSL (a sin a + b sin 2o),
(14)
gp = zT (c0 + c1 cos a + c2 cos 2a) qSL sin p -
- D
qSPa V
Полином f (u) при этом имеет четыре действительных корня, два из которых совпадают u1 = u4 = u0. Перенумеруем корни:
a.) u1 = v0, u2 = v1, u3 = v2, u4 = v0; (10)
b.) u1 = v0, u2 = v2, u3 = v1, u4 = v0. (11) Вариант (a) соответствует траектории a01, а
вариант (b) - траектории a02. Подставляя корни (10) и (11) в решение (8) и учитывая, что cos у = cn(Pt + г0,1) = ch-1 (в + 70) .
■'а I , 8а
где р = 1а - обобщенный импульс. Величины q и V, G и л будем считать постоянными, а в качестве закона изменения угла аэродинамического крена - использовать зависимость р = XI,
где X - осредненная по периоду колебания угла атаки угловая скорость
- T if fR-
(G - R cosa) cosa
Л
sin2 a
J
О '-"'max / T
Ul-fR
= R
I
n ^ f--?
J* j
(G - R cos a) cosa sin2 a
_ 2 (GJ - RJ2)
T '
dt =
da a
(15)
= r cosa da = f
J1 = J "; J 2 =1
J cm Г/ Г/ ' j
da
sin2 a a
a sin2 a '
max
T = 2
J a '
(16)
Для вычисления интегралов (16) необходимо представить функцию (7) в виде
/(и) = (и - Щ )(щ - и2)(и - и3 )(и - и4 )
для случая, когда / (и) имеет только действительные корни;либо в виде
f (u) = (u - u1)(u - u2) ((u - u34)2 +v2)
для случая, когда присутствуют комплексные корни. Воспользовавшись заменой u = cosa и выражениями (8) или (9), перейдем в интегралах (16) от переменной a к переменной Y. Пос-
a
a
а
a
a
ле алгебраических преобразований подынтегральных выражений, (16) сводятся к полным нормальным эллиптическим интегралам Лежан-дра первого К (к) и третьего П (п, к ) рода [11]
=в ( П (П1, к ) + СП (п2, к) + С3 К (к)), Зг =в (С, П (п,, к )-С2 П (п2, к ) + С4 К (к)),
где
K (k) = j
dy
ф - k2 si
sin у
П (n, k) — j-
dy
o (l + и sin2 y)) - k2 sin2 y
В случае, когда все корни (7) действительные
T —
2K (k) ^ = u2 - u3
в ; 1 2(1 - u2)(1 - u3)
C = u2 u3
2(1 + u2)(1 + u3)
C = u3
, . 2
1 - u
C4 =
1
1 - u
2 , 1
n=
(u2 - u1)(1 - u3)
(u1 - u3)(1 - u2)
n2 =
(u2 - u1 )(1 + u3)
(u1 - u3 )(1 + u2) '
Для случая комплексных корней T = 4K (k) с =. ux - u2 X1 - X3
в
2(1 -u1)(1 -u2) x2 - x4
C u 2 x + X3 с 2 X2 X4
2(1 + u1)(1 + u2) x2 + x¿
3 2 _ 2 ,
4 x2 - x4
C4=
2 x2
4 „2 „2 , "1
n=
(x2 - x4 )2
4a1a2(1 - u1)(1 - u2)
n2 =
(x2 + x4 ) =
4a1a2(1 + u1)(1 + u2) , x1 = a1 + a2,
x2 — a a2, x3 — u2a1 + u1a2, x4 — u2a1 u1a2.
Запишем критерий Мельникова [10] для системы (14).
Мг (to) — £ fa (a0, (t))gp (a0i(t), t + to)- fP (a°t (t) )ga(a\ (t), t + to )d t, (17)
где gp (a0, (t), t +10) - вектор возмущений, вычисленный на i -ой гомоклинической траектории. Подставляя в (17) f и g из (14), получим
Mt (t0) — qSLzaO (t) sin ( (t +10)) x x (c0 + c1 cos a° (t) + c2 cos 2a° (t)) -
- DÍSSOtL d,
V
(18)
Здесь ) и а) вычисляются по формулам (12), (13) соответственно.
Если функция Мельникова Mi (¿0) имеет простой корень в точке t0, то при достаточно малом в устойчивые и неустойчивые многообразия пересекаются трансверсально. Если функция Mi ) не имеет простых корней, то устойчивые и неустойчивые многообразия не пересекаются и хаотических режимов, обусловленных пересечением гомоклинических структур не возникает. Другими словами критерий Мельникова позволяет судить об устойчивости областей движения системы (14). В случае, если функция Мельникова имеет простые корни, то в системе присутствует хаос и фазовая траектория может пересечь сепаратрису и покинуть свою область.
ХАОТИЧЕСКИЙ АНАЛИЗ
Рассмотрим движение КА, имеющего следующие массово-геометрические параметры: Ь = 0.35 м, £ = 2.8 м3, 1Х = 1 кг• м2, I = 3 кг • м2, гТ = 0.02 м. Аэродинамические коэффициенты КА зададим с помощью выражений (4) и (5), в которых а = 0.036 , Ъ = -0.178 с0 = 0.069, с, = 1.112, с2 = 0.206 (рис. 1). Построим критерий Мельникова (18) для системы (14), имеющей следующее параметры Я = 0.5 с-1, О = 0.2 с-1, Я = 2 с-1, V = 1000 м/ с, q = 2000 кг/мс2.
При П = 0 М1 и М2 обращаются в ноль в нескольких точках. В этом случае в системе возникает хаос. При П = 0.05 М2 не обращается в ноль, то есть в окрестности второй гомоклинической траектории хаоса не возникает (рис. 4). На рис. 5, 6 показаны сечения Пуанкаре, построенные для фазовых траекторий, начинающихся только в областях А1 и А2. Область А2 устойчива: фазовые траектории, начинающиеся внутри этой области, не покидают ее (рис. 5). Область А1 неустойчива (подвержена хаосу): некоторые фазовые траектории, начинающиеся в этой области, пересекают сепаратрису и попадают в другие области движения (рис. 6). При П = 0.6 ни М1 ни М2 не обращаются в 0 и в окрестности гомоклинических траекторий хаоса не возникает, то есть обе области устойчивы.
2
,Ví
ч 1 ■W, А
\ / \ А А
\ / \ i
V У У У У
м,
12
15 ¡.с
Рис. 4. Критерий Мельникова для областей Л1 и Д , вычисленный при В = 0.05
10
-10
. _ _
1
А?
>>1 у 1 % % я,. ,
1 \ ___ - --
.т'4
з-.ч
Рис. 5. Сечение Пуанкаре для фазовых траекторий, начинающихся в области A2 при D = 0.05
Поскольку хаотический анализ выполнялся для фиксированного момента времени, следует определить точку на траектории спуска, в которой КА наиболее подвержен хаосу. Критерий Мельникова Мг (t0) характеризует толщину хаотического слоя в окрестности сепаратрисы [10]. Изучим, как меняется толщина этого слоя в процессе спуска КА. Пусть рассмотренный ранее КА осуществляет спуск в атмосферу Марса [12] при следующих начальных условиях: а0 =п/6, а0 = 0 с-, V0 = 1000 м/с , H0 =120 км , в0 = -п/12 . Зададим коэффициент демпфирования D = 0.01. В качестве приближенной характеристики толщины хаотического слоя будем рассматривать величину
Лг = max Мг (t0) - min Мг (t0) .
На рис. 7 показано изменение Лг в процессе спуска. Максимальная толщина хаотического слоя наблюдается при максимальном значении скоростного напора (при t = t * рис. 8). Выбирая коэффициент демпфирования D достаточно большим для подавления хаоса при t = t *, мы тем самым устраним его на всей траектории спуска. Значение параметров системы в момент t = t * можно получить, рассмотрев движение центра масс КА при спуске в атмосфере без учета движения относительно центра масс.
Рис. 6. Сечение Пуанкаре для фазовых траекторий, начинающихся в области Л1 при В = 0.05
I..Í i
ы t
/У V -
// д" —-
ion i:n míi i С-.П iw :íiü ;;<: 210 2№
2 к n
Рис. 7. Изменение толщины хаотического слоя в процессе спуска
i20 ico 2-ю
ix
i'mc. о. изменение скоростного напора при неуправляемом спуске КА
ЗАКЛЮЧЕНИЕ
В работе предложен новый, основанный на методах хаотической динамики, подход к исследованию пространственного движения КА с малой асимметрией при спуске в атмосферу. Получены уравнения гомоклинических траекторий нелинейной системы дифференциальных уравнений движения КА. На основании метода Мельникова записан критерий устойчивости движения в окрестности гомоклинических траекторий. Исследовано изменение толщины хаотического слоя в процессе спуска КА в атмосфере и даны рекомендации по выбору уровня демпфирования, устраняющего возможность появления ха-
отических переходов на всей траектории спуска. Несмотря на то, что в работе рассматривался лишь один из видов асимметрии, предложенный подход может быть легко распространен и на другие ее виды.
Работа выполнена при поддержке Российского фонда фундаментальных исследований (0601-00355).
СПИСОК ЛИТЕРАТУРЫ
1. Sancisi-Frey S., Spry J.A., Garry J., Pillinger J.M. Atmospheric entry simulations of Mars lander bioload -experiments in support of Beagle 2 // Research in Microbiology. 2006. Vol. 157. pp. 25-29.
2. Gershman R, Adams M, Mattingly R, Rohatgi N, Corliss J, Dillman R., Fragola J., Minarick J. Planetary protection for Mars sample return // Advances in Space Research. 2004. Vol. 34. pp. 2328-2337.
3. Collinet J., Brenner P., PalermS. Dynamic stability of the HUYGENS probe at Mach 2.5 // Aerospace Science and Technology. 2007. №11 pp. 202-210.
4. Shotwell R. Phoenix—the first Mars Scout mission //
Acta Astronautica. 2005. Vol. 57. pp. 121-134.
5. Асланов В.С, Ледков А. С. Особенности вращательного движения КА при спуске в атмосфере Марса // Космические исследования. 2007. Т 45, № 4 С. 351-357.
6. Асланов В.С. Пространственное движение тела при спуске в атмосфере. М.: Физматлит, 2004. 160 с.
7. Wiggins S, Holmes P. Homoclinic orbits in slowly varying oscillators // SIAM J. Math. Anal. 1987. Vol. 18. № 3. pp. 612-629.
8. Holmes P. and Marsden J. Horseshoes in perturbations of hamiltonian systems with two degrees of freedom // Commun. Math. Phys. 1982. Vol. 82. pp 523-544.
9. Faouzi Lakrad, Moulay Mustapha Charafi Perturbation methods and the Melnikov functions for slowly varying oscillators // Chaos, Solitons and Fractals. 2005. Vol. 25. pp. 675-680.
10. Симиу Э. Хаотические переходы в детерминированных и стохастических системах. Применение метода Мельникова в технике, физике и нейрофизиологии. - М.:ФИЗМАТЛИТ, 2007. - 208с.
11. Янке Е., Эмде Ф, Леш Ф. Специальные функции (формулы, графики, таблицы). М.: Наука, 1968. 344с.
12. Мороз Б.И, Кержанович В.В, Краснопольский В.А Инженерная модель атмосферы Марса для проекта Марс-94 (МА-90) // Космич. исслед. 1991. Т. 29. № 1. С. 3-19.
THE CHAOTIC ANALYSIS OF MOVEMENT OF THE SPACECRAFT WITH SMALL MASS ASYMMETRY AT DESCENT IN AN ATMOSPHERE
© 2009 V.S. Aslanov, A.S. Ledkov
Samara State Aerospace University
The uncontrolled motion of the spacecraft with small mass asymmetry at descent in atmosphere is investigated by the methods of chaotic dynamics. The restoring moment acts on the spacecraft. It approximated by biharmonic expression. With the help of Melnikov method the behaviour of a system in neighbourhood of the separatrix is investigated. Analytical solutions for homoclinic trajectories are received. Criteria of occurrence of chaos are found. The results of the numerical modelling are show validity of received solutions. Change of thickness of a chaotic layer is investigated during descent of the spacecraft at atmosphere. Recommendations are given at the choice of the damping eliminating a capability of occurrence of chaotic transitions on all trajectory of descent.
Key words: spacecraft, biharmonic moment, chaos, Melnikov method, homoclinic trajectory, the analytical solution, a chaotic layer, stability, Poincare sections.
Vladimir Aslanov, Doctor of Technics, Professor, Head of
Theoretical Mechanics Department.
E-mail: [email protected]
Alexander Ledkov, Graduate Student.
E-mail: [email protected].