Численное моделирование и компьютерный анализ устойчивости орбитального движения космических аппаратов

Буланов Сергей Георгиевич; Джанунц Гарик Апетович

Раздел III. Информатика

УДК 517.91: 518.1 ББК 22.193

С.Г. Буланов, Г.А. Джанунц

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

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

Ключевые слова: компьютерное моделирование устойчивости; решение обыкновенных дифференциальных уравнений; кусочно-полиномиальная аппроксимация.

S.G. Bulanov, G.A. Dzhanunts

NUMERICAL MODELING AND COMPUTER STABILITY ANALYSIS OF THE ORBITAL MOTION OF THE SPACECRAFT

Abstract. The results of numerical modeling and stability analysis of the orbital motion of the spacecraft in the hypothesis of an amendment to the Law of universal gravitation are presented. Results of the research presented in the form of graphic illustrations of solutions of the system, and the trajectory of the dynamic changes in the nature of management.

Key words: computer modeling of stability; solution of ordinary differential equations; piecewise polynomial approximation.

К настоящему времени механика полета космических летательных аппаратов (КЛА) с малой тягой превратилась в самостоятельный раздел механики космического полета. Система уравнений движения КЛА в плоскости орбиты имеет следующий вид [2]:

'г ( t ) = Vr,

h2

Vr ( t ) = V© r4--r -2 + Ur,

P (1)

V© ( t ) = -VrV© r-1 + U©,

где r , 6 - полярные координаты, Vr, V© - радиальная и трансверсальная составляющие скорости; Ur. U© - составляющие вектора тяги, a, b - большая и малая полуоси эллипса, p - фокальный параметр, e - эксцентриситет эллипса. Для данных параметров имеют место следующие соотношения:

2 , Г. 2 , 2nab h1 , a - b

p = (1 -e2)a, b = J1 -e2 , h = —, — = GM, e =- .

T p a + b

Управляющие воздействия считаются малыми по сравнению с минимальным значением си-

I—2-Г h2

лы тяготения, т.е. JUr + U0 <<-— . При Ur = U0 = 0 и отрицательной полной энергии E < 0

pr

г max

система (1) описывает движение КЛА по эллиптической орбите под действием известного закона всемирного тяготения Ньютона [2].

Законы Кеплера можно записать в виде следующих инвариантов движения:

га 1 = r (1 + e cos 0) - p = 0, (2)

ю 2 = r2 0(t)-h = 0. (3)

Синтез законов управления Ur и U© осуществляется на основе метода АКАР [2]. В качестве инвариантных многообразий выбираются со j(t) = 0 и ra 2(t) = 0. В полярных координатах функция сс j(t) с учетом инварианта ra j из (2) имеет вид

сС j (t) = pr -1 r (t) - er © (t) sin © = pVrr - - eV© sin © . (4)

Имеют место следующие инвариантные соотношения:

Со j(t) + Фсс j(t) = 0, (5)

сс 2(t) + Фю 2(t) = 0. (6)

Совместное решение (5), (6) с учетом (3), (4) влечет следующие законы управления:

h2 (ю2 + h)2 ф/ , ч

Ur =—2--2---(reo j( t) + ею 2sin ©), (7)

pr pr p

ю 2

U©=--2 Ф . (8)

r

На пересечении инвариантных многообразий сс j (t) = ю 2 = 0 законы (7), (8) «обнуляются» и, следовательно, КЛА будет устойчиво двигаться вдоль заданной орбиты - предельного цикла.

При этом важно, что эти законы на пересечении инвариантных многообразий сс j (t) = 0 и

с 2 (t) = 0 редуцируются в известный закон всемирного тяготения Ньютона. Следовательно, законы Ur из (7) и U© из (8) можно назвать динамическими законами гравитационного взаимодействия двух тел.

h

В качестве функции Ф (t) может быть выбрана [3] функция Ф (t) = —г. Результаты иссле-

r

дования представлены для случая сс j( t) Ф 0, ra 2( t) Ф 0.

Для управления движением летательных аппаратов необходимо в режиме реального времени находить приближенное решение системы (1) с требуемой точностью и выполнять анализ устойчивости. Для приближенного решения задачи (1) целесообразно использовать метод кусочно-полиномиального решения задачи Коши для дифференциальных уравнений (ДУ) вследствие малой временной сложности, высокой точности и непрерывного характера приближения [1, 5]. Ниже для ясности кусочно-полиномиальный метод кратко описан для решения ДУ первого порядка.

Пусть требуется найти приближенное решение задачи Коши

^ = f (х, y), y(x0) = y0. (9)

dx

на произвольно фиксированном отрезке [ a, Ъ ]. Предполагается, что на [ a , Ъ ] выполнены все условия существования и единственности решения задачи (9), а также условия сходимости метода Эйлера. Задается система подынтервалов равной длины

p-1

[a, Ъ] = U [х, х+1], P = 2k, k = {0,1,...}.

i=0

Последовательно на каждом подынтервале, начиная с [ х0, х 1], строится интерполяционный полином Ньютона для функции правой части (9), который переводится в форму полинома с числовыми коэффициентами по алгоритму, отличному от формул Виета [5]. В качестве значений зависимой переменной первоначально используются разностные приближения, полученные по методу Эйлера

y¡p = y(p-1) + h•f (X (py(p-ц^ p = 1 n .

На i -м подынтервале полиномиальное приближение функции правой части (9) имеет вид:

(х) = ai 0 +Z

( х - х ^

i 0

"i t

t=1

h

(10)

Полином (10) приближает производную искомого решения задачи (9). Приближение самого решения строится как первообразная от (10) при условии подстановки в качестве постоянного значения решения на левой границе текущего подынтервала. Для подынтервала с номером i = 0 начальное значение известно из условия задачи Коши - у00 = у 0, для подынтервалов [х1, х1+1],

i = 1,...,2к -1, оно принимается равным значению приближенного решения в конце предыдущего подынтервала - у10 = у-1 п.

t

Таким образом, искомое приближение решения задачи на i -м подынтервале имеет вид

(х) = У 0 + I V,» (х) ^,

или,

zi (X) = yt 0 + h Ё

t + 1

X - X,.

(11)

Вычисление значений полинома (11) производится по схеме Горнера при t =

( ((

Z (t) = yt 0 + h

ai {n-1)

V v

n +1

t+

ai {n-2)

n - 1

t +... + a,,

t.

(12)

Значения у р = zi (х1р), р = 1, п , вычисленные из (12) оказываются более точными приближениями решения, чем получаемые непосредственно с помощью метода Эйлера. Эти значения принимаются за новые приближения в узлах интерполяции для последующего на их основе построения полинома (10). Данный рекуррентный процесс позволяет существенно уточнить полученные приближения.

Аналогичное приближение строится на следующем подынтервале и т.д., до исчерпания отрезка [ а, Ъ ].

Степень полинома и количество подынтервалов 2к варьируются в априори заданном диапазоне значений, в результате фиксируются оптимальный набор значений для конкретной задачи.

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

Ш—18

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

методов [1, 5].

Моделирование устойчивости решения системы (1) выполняется ниже на основе критерия устойчивости для систем вида:

dY

d" = F(t, Y),

(13)

Y (to) = Y.,

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

представленного в [4].

При анализе устойчивости решения задачи (13) предполагается, что существует 8 0 > 0, такое, что для всех невозмущенных и возмущенных решений системы (13), начальные условия которых удовлетворяют условию ||y0 - Y0 || <8 0, выполнены все условия существования и единственности на полуоси [t0, да). Дополнительно предполагается, что k -я компонента функции F (t, Y), к = 1,..., n , удовлетворяет аналогу условия Липшица с единой константой на произвольном промежутке полуоси:

\fk (t, Y) - fk (t, Y)| < L\yk - , L = const, V t е[ T ]aVT е( да). В рассматриваемых условиях для устойчивости решения задачи (13) необходимо и достаточно существование А1, 0 <Д1 <8 0, такого, что для всех решений Y = Y(t), Y (t0) = Y0 при ог-

раничении 0 < Y0 - Y0 < Д j выполняется неравенство

~k (t) - Ук (t)

~k 0 - Ук 0

<с , с = const, Vtе [t0, да), k = 1,...,n .

(14)

Для асимптотической устойчивости в тех же условиях необходимо и достаточно, чтобы выполнялось (14) и существовало Д2 <Д^ такое, что неравенство 0 < - Y0 <Д 2 влечет

lim

t ^да

~ (t) -Ук (t)

= 0, к = 1,...,п .

Ук 0 — Ук 0

Для большей достоверности компьютерного анализа устойчивости приближенные значения ук (?) и ук (?) в (14) целесообразно вычислять на основе описанного выше кусочно-полиномиального метода.

X

i0

t + 1

a

t

h

t=0

X — X

i 0

h

a

n

Моделирование решения системы (1) выполняется на промежутке [0; 2 106] при начальных условиях

г(0) = г,, у(0) =у, 0(0) = 0, (0) = 0. (15)

При этом для констант из (1) задаются значения [3]:

р ек

р = 36000км, е = 0,5, к = 145106 км /ссе, г0 =-км, У0 = — км/сек.

1 + е р

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

Для управления Ф (t) = -— система (1) преобразуется в систему

г

' * = V

аг г'

с1Уг V®2 кУг ек2 81И(0)

аг

а0 = у аг г

г2 р

(16)

аУ0 = уу к¥0 +

аг г г2 г3

и -У0-

г 2

рг

2

Р

к¥г ек2 зш(0)

~ + 2 г г р

и& = -

кУ 0 к2 г2 г

(17)

Первоначально численное решение системы (16) выполняется с шагом интегрирования равным 1 сек. Это делается для того, чтобы получать результаты за приемлемое время. При этом теоретически может неконтролируемо расти погрешность решения, однако результат совпадает с решением масштабированной системы, которое приводится в дальнейшем.

Результаты моделирования решения иллюстрируют рис. 1-7.

Рис. 1. Изменение радиуса г ( ?) из (16) при начальных условиях (15)

Рис. 2. Изменение скорости Уг (?) из (16) при начальных условиях (15)

Рис. 3. Изменение 0 ( ?) из (16) при начальных условиях Рис. 4. Изменение скорости У0 ( ?) из (16) при начальных

2

г

(15)

условиях (15)

Рис. 6. Изменение иг (?) из (17) при начальных условиях (15)

Рис. 7. Изменение и0 (?) из (17) при начальных условиях (15)

При е = 0 результаты моделирования иллюстрируют рис. 8-14.

"I

0.5 1.0 1.5 2.0

Рис. 8. Изменение радиуса г ( ?) из (16) при начальных Рис. 9. Изменение скорости Уг ( ?) из (16) при начальных

условиях (15)

Рис. 10. Изменение 0( ?) из (16) Рис. 11. Изменение скорости У0 (?) из (16)

при начальных условиях (15) при начальных условиях (15)

Рис. 12. Фазовый портрет траектории КЛА, описываемого системой (16) при начальных условиях (15)

Рис. 13. Изменение иг (?) из (17) Рис. 14. Изменение и0 (?) из (17)

при начальных условиях (15) при начальных условиях (15)

Анализ устойчивости системы (16) выполняется на основе критерия устойчивости (14). Согласно эксперименту [3] для достоверности анализа устойчивости на основе (14) шаг численного интегрирования необходимо выбирать порядка 10—8 -10—6 или меньше. Это требует значительных временных затрат при моделировании на промежутке [0; 2 • 106 ]. Выбор величины шага порядка 0.01 -1 позволяет выполнить анализ устойчивости и найти приближенное решение системы при меньших временных затратах, но при этом может неограниченно расти погрешность. С целью сведения погрешности разностного решения к минимуму система (16) масштабируется. Вводятся новые переменные:

R1 =—, V =

.©1 =

©

©I

V,

t

А = —, 1 т

(18)

где г V

та^' г тах ^6

К.

- максимальные значения переменных г, Уг, ©, V© из (16),

Т = 2-106.

Выражая из (18) переменные г , Уг, © , V© и подставляя их в (16), получим следующую систему:

тах 2

^ = V, тах Т V 1'

й©1 V

V2 ТУ

у© тах у2

г V

тах г тах

TУ

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

У©2max Т У2 кТУх ек1 Т б1П (©тах©1 ) Т

PУг т

©тах гтах Л1

V У2

гтах Л1

г2ах Л2

к Т

pУ г2 Л,2

г г тах тах 1

hTУ1

г2 Л2 г3 V Л

тах 1 тах © тах 1

(19)

Результаты численного анализа устойчивости системы (19) при соответственно преобразованных начальных условиях

г V © V

г V © Уr

(20)

представлены в таблице 1.

Таблица 1

Результаты анализа устойчивости системы (19) при начальных условиях (20)

г

Начальное возмущение Приближение к тах 1< к < я ~к ( А ) - Ук ( А ) Ук 0 — Ук 0 по шагам разностной схемы

5 = 5 1 0.1 0.2 0.3 0.4 0.5

погта 27.09 73.63 119.76 165.86 211.75

1 0.6 0.7 0.8 0.9 1

погта 257.60 303.53 349.37 395.13 440.99

5 = 2 1 0.1 0.2 0.3 0.4 0.5

погта 1.04 1.04 1.04 1.04 1.04

1 0.6 0.7 0.8 0.9 1

погта 1.04 1.04 1.04 1.04 1.04

5 = 1,5 1 0.1 0.2 0.3 0.4 0.5

погта 1.05 1.05 1.06 1.05 1.07

1 0.6 0.7 0.8 0.9 1

погта 1.05 1.05 1.07 1.05 1.06

5= 1,2 1 0.1 0.2 0.3 0.4 0.5

погта 1.47 1.03 1.03 1.04 1.03

1 0.6 0.7 0.8 0.9 1

погта 1.03 1.04 1.03 1.03 1.04

t 0.1 0.2 0.3 0.4 0.5

погта 1.04 1.34 1.18 1.19 1.34

5 = 1 t 0.6 0.7 0.8 0.9 1

погта 1.04 1.22 1.34 1.04 1.21

t 0.1 0.2 0.3 0.4 0.5

5 = 0,1 погта 11.62 10.16 6.39 14.86 10.31

t 0.6 0.7 0.8 0.9 1

погта 8.90 15.33 10.47 10.41 15.36

t 0.1 0.2 0.3 0.4 0.5

5 = 0,001 погта 263.02 209.79 179.69 156.55 195.66

t 0.6 0.7 0.8 0.9 1

погта 178.84 208.00 203.24 181.31 206.78

При возмущении из диапазона 5 < 5 решение системы (16) устойчиво, но не асимптотически.

Таблица 2

Результаты анализа устойчивости системы (19) при начальных условиях (20) ( е = 0 )

Начальное возмущение Приближение к тах 1< к < п ~к (t) - Ук (t) Ук 0 — Ук 0 по шагам разностной схемы

5 = 5 t 0.1 0.2 0.3 0.4 0.5

погта 1.03 1.02 1.02 1.02 1.02

t 0.6 0.7 0.8 0.9 1

погта 1.02 1.02 1.02 1.02 1.02

5 = 2 t 0.1 0.2 0.3 0.4 0.5