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

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

Т о м V 1 9 7 4 №5

УДК 521.31

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

М. Ю. Беляев, В. П. Семенко

Рассмотрены различные способы описания возмущенного движения для прогнозирования движения центра масс спутника с помощью ЭЦВМ. Проведен сравнительный анализ известных систем уравнений [2, 3] и систем, полученных с помощью вращающейся системы координат [4, 5]. Показано, что при одинаковых затратах машинного времени на интегрирование точность прогнозирования с помощью новых систем уравнений для большинства практически важных случаев выше по сравнению с другими рассмотренными системами. Приведены результаты численных расчетов на ЭЦВМ.

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

Возможны два подхода к решению этого вопроса:

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

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

Оставляя в стороне вопросы, связанные с первым способом повышения скорости расчета движения КА, рассматриваемые, например, в [ 1 ], займемся вторым способом.

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

В настоящее время при решении траекторных задач наиболее широко применяются уравнения движения спутника в декартовой прямоугольной системе координат и уравнения в оскулирующих переменных [2, 3). Наиболее простой является система уравнений в декартовых координатах. Однако для обеспечения сравнительно высокой точности расчета движения спутника с помощью этой системы необходимо выбирать небольшой шаг интегрирования. Последнее обстоятельство ограничивает возможность существенного повышения быстродействия расчета движения спутника в этой системе координат.

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

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

Известны также уравнения движения в цилиндрических координатах [2]. Уравнения движения в цилиндрических координатах используются реже, чем УОЭ или уравнения в декартовых координатах. Однако в [2] отмечается, что эти уравнения являются наиболее выгодными при расчете орбит с малыми эксцентриситетами.

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

В настоящей работе используем форму уравнений движения спутника в ВСК, отличную от полученной в [4]. По аналогии с методом оскулирующих элементов расстояние от начала координат до спутника будем определять по формуле

Я = -, --Р-.----г. (1)

I -г- е cos (и — ши)

где р, е и —параметр, эксцентриситет и аргумент перицентра эллипса в ВСК соответственно; и — аргумент широты.

Отметим, что вследствие специального выбора вращения ВСК скорости спутника в ней и в инерциальной системе координат не равны и, следовательно, эллипс в ВСК будет отличен от соответствующего оскулирующего эллипса.

Параметр р вследствие специфики вращения ВСК постоянен. Вывод уравнений для е и ш„ аналогичен выводу уравнений для соответствующих элементов в системе УОЭ [3]. Поэтому, опуская преобразования, сразу выпишем дифференциальные уравнения для е и ши:

= Jl(S + *)sin (и - „>„);

-^f- ~ — Т~ (S "+" X) cos (и-,ы) + at Le

(2)

где 5—радиальная составляющая возмущающей силы, действующей на спутник;

X = —(XjL —j— Xj — 2 Lkz)\ \Xy, Xz — параметры, определяющие вращение ВСК [4, 5];

R3 *

ф = и — v; V — угол между осью Ох ВСК и радиус-вектором спутника; L — кинетический момент спутника в ВСК.

Использование уравнений (2) вместе с соотношением (1) позволяет исключить дифференциальное уравнение второго порядка для R из системы,полученной в работе |4]. Эта замена оказывается весьма важной при числениом интегрировании уравнений, так как переменные е и ш„ меняются намного медленнее по сравнению с переменной R, что предпочтительнее с точки зрения увеличения шага интегрирования, но в этом случае уравнение для <•>„ имеет особенность при е = 0, которую устраним введением следующих переменных:

■г] = е sin ш„; % = е cos wu. Аналогичную замену переменных обычно проводят и в УОЭ. Дифференциальные уравнения для т) и £ можно получить с помощью (2). Объединяя их с остальными уравнениями ВСК, получим следующую систему:

Í = — -j- (S + X) cos u — £ cos ö -f ^ j;

(3)

R2 J

'xy = FnzR\ К = Fnxy R\ Ь = -^sin и;

' L — А» ■ 9 = —--—cosu; u—-£ — 9 cos it,

где в, ? — углы, определяющие положение плоскости хОу ВСК в инерциальной

системе координат Ох0 y^z^ Fnxy, Fnz — проекции силы, ортогональной R, на плоскость хОу и ось Ог соответственно.

Для сопоставления системы уравнений (3) с другими известными системами уравнений были проведены сравнительные численные расчеты на ЭЦВМ. Рассма-

тривались следующие системы дифференциальных уравнений движения спутника:

А — уравнения движения КА в декартовой прямоугольной системе координат [2];

Б —система уравнений в оскулирующих элементах (УОЭ) [2]. Использовалась система УОЭ, в которой уравнения для аргумента перицентра и эксцентриситета заменяются уравнениями для компонентов вектора Лапласа ^ Эта замена, как известно, устраняет трудности, возникающие при расчете орбит, близких к круговым;

В — система уравнений движения в цилиндрических координатах [2];

Г — система (3).

С целью выяснения возможностей метода оскулирующих переменных и метода ВСК исследовались различные модификации систем УОЭ и систем уравнений движения, записанных с помощью ВСК. Из рассмотренных систем УОЭ для дальнейшего сравнительного анализа была выбрана система Б, поскольку она является в рассматриваемом случае достаточно универсальной, так как позволяет рассчитывать большой класс орбит со сравнительно высокой точностью. Аналогичным образом из систем, записанных с помощью ВСК, была выбрана система Г.

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

За меру точности определения положения спутника бралась величина

д/г = 1я.

(4)

В качестве эталонной принималась орбита, рассчитанная с небольшим шагом интегрирования.

Интегрирование всех систем проводилось методом Рунге—Кутта с постоянным шагом. Величина шага для каждой системы выбиралась в зависимости от вида правых частей дифференциальных уравнений таким образом, чтобы времена интегрирования'на одинаковых интервалах полета спутника были равны.

В качестве возмущения рассматривались вторая гармоника в разложении геопотенциала в ряд по сферическим функциям и сопротивление атмосферы. Выражения для составляющих возмущающей силы, действующей на спутник в указанных случаях, имеются, например, в [3]. Значение баллистического коэффициента спутника бралось равным 0,003 м2/кг.

На точность прогнозирования движения спутника влияют интервал и шаг интегрирования.

Увеличение шага интегрирования оказывает различное влияние на точность интегрирования рассматриваемых систем на одинаковом интервале. Кроме того» увеличение интервала может существенным образом изменить картину сравнения. Так, например, использование системы В для численного прогнозирования на интервале, не превышающем 3—5 суток, дает сравнительно высокую точность. Однако при увеличении интервала интегрирования до 8—10 и более суток точность расчета с помощью системы В резко падает. Это связано с тем, что одна из интегрируемых функций — величина выхода спутника из базовой плоскости, проходящей через начальные векторы скорости и положения спутника, становится, вследствие эволюции орбиты, довольно большой, что приводит к резкому понижению точности при интегрировании системы В со сравнительно большим шагом. В работе для заданного интервала полета определялось максимальное

значение ДЛ? в зависимости от времени интегрирования каждой системы (которое, очевидно, зависит от шага интегрирования).

Сравнительные расчеты проводились на ЭЦВМ М-220 для орбит со следующими начальными значениями эксцентриситета, наклонения и параметра, представленными в табл. 1.

На фиг. 1—4 показаны максимальные значения Д/? для трех суток полета спутника в зависимости от времени, затрачиваемого на интегрирование уравнений движения в течение этого интервала с помощью каждой из систем.

Таблица 1

№ орбиты е 1 р, км-

1 0,0018 51,51° 6634

11 0,2 51,51° 7900

111 0,05 63,43° 6996

IV 0,0018 Г 6634

л Я/* Орби та Е

Фиг. 2

На основе результатов численных расчетов можно сделать некоторые выводы. Системы уравнений А и Б в рассмотренном диапазоне интегрирования характеризуются наименьшей точностью (из анализируемых четырех систем). Из рассмотрения фиг. 1—4 видно, что эти две системы дают меньшие точности по сравнению с системой В для всех рассмотренных орбит. Для всех орбит, кроме орбиты II, эти системы также значительно хуже системы Г. Таким образом, системы уравнений А и Б менее всего подходят для численного прогнозирования движения спутника. Любопытно отметить, что именно эти системы (особенно система А), по-видимому, наиболее широко используются при численном прогнозировании орбитального движения спутников [2].

В работе [2] отмечается, что для орбит, имеющих эксцентриситет использование уравнений движения в цилиндрических координатах (система В) позволяет существенно повысить быстродействие расчетов. Из фиг. 1—4 видно, что в случае использования системы В выигрыш во времени интегрирования по сравнению с системами А и Б достигается для всех рассмотренных случаев. Однако в целом этот выигрыш, особенно по сравнению с системой Б, невелик. Кроме того, как показали соответствующие расчеты, при увеличении интервала прогнозирования точность расчета движения спутника с помощью этой системы, в силу указанной выше причины, резко падает. Системы же А, Б, а также Г лишеиы этого недостатка.

Наиболее предпочтительной из всех здесь рассмотренных систем для расчета орбит с малыми и умеренными эксцентриситетами оказывается система уравнений Г (см. фиг. 1, 3, 4).

Таким образом, можно сделать следующий вывод.

Фиг. 3

Фиг. 4

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

Теперь остановимся на расчете орбит, имеющих большой эксцентриситет. Покажем, что простой переход к новой независимой переменной в системе (3) позволяет повысить точность численного прогнозирования и для этого класса орбит. В дифференциальных уравнениях движения (3) перейдем к новой независимой переменной ак—невозмущенному аргументу широты спутника в ВСК. Используя уравнение

= А., (5>

- Л

получим:

di} dak

di duk

-f-(S + X) cos« + P

77 (S + .Y) sinu + ii^j-

„ rfcp * du ь

cos i

rf®

dX rV

= Л: ff;

dak

d\z dl=

_= 7 —_ sin u:

duk R*

¿ф. = _. 7__

duk R2 sin $

cos «;

du

А - Аг ff*

d<? dub

cos ft ,

(6)

где -г = Rl|L .

Переход к независимой переменной ик ставит реальному движению спутника в соответствие кеплерово движение. Эта замена позволяет исключить время из системы уравнений и определять его с помощью уравнения Кеплера [3]. Отметим, что аналогичный переход к новой независимой переменной может быть осуществлен и для системы УОЭ. Однако переход к новой переменной в ВСК в силу ее особенностей дает некоторые преимущества. Специальные проверки показали, что основная погрешность численного прогнозирования с помощью системы (3) связана с интегрированием последнего уравнения этой системы (или аналогичного уравнения для быстроменяющейся переменной в случае использования УОЭ [6]). Ясно, что возможность увеличения шага при численном интегрировании связана с поведением интегрируемой функции В случае использования УОЭ поведение этой функции для какого-либо заданного движения спутника определяется единственным образом. При использовании же неинерциальной ВСК в силу переопределенности этой системы [4, 5] выбором начальной ее угловой скорости можно получать различные траектории спутника в ней, т. е. тем или иным образом влиять на поведение интегрируемых функций.

Основная погрешность прогноза с помощью системы (6) также обусловлена интегрированием последнего уравнения этой системы. Интегрируемая функция и является суммой быстрой переменной V и медленно меняющейся Ол Дифференциальное уравнение для переменной V записываем в виде

dah

R

где R/г — радиус кеплеровой траектории спутника.

Отсюда ясно, что желательно иметь такое движение спутника в ВСК, при котором Rk близко к R. В [5| показано, что кеплеров эллипс, получающийся с помощью ВСК, оказывается существенно ближе к реальной траектории спутника, чем кеплеров эллипс, полученный обычным способом в неподвижной системе координат. Для оценки возможности использования системы (6) для численного прогнозирования были проведены сравнительные расчеты на ЭЦВМ

следующих вариантов орбит, данных

Таблица 2

№ орбиты е 1

V 0,5 9 900 51°, 51

VI 0,7 11 220 63°, 43

VII 0,85 12210 1 =

в табл. 2.

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

йК^

1200-

Орбиты Т,ЛГ, Ш

1000

600

500

т

200

д~т

О 40 80 ¡20 160 200 2М 200 320 с, с

Фиг. 5

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

Результаты сравнительных расчетов приведены на фиг. 5. Цифра у каждой кривой на этой фигуре соответствует номеру варианта, а буква — используемой системе, причем система (6)обозначена буквой Д. График Б-У1 на фигуре отсутствует, поскольку погрешность расчета орбиты VI с помощью системы Б слишком велика.

Из сравнения приведенных результатов следует, что точность расчета движения спутника с помощью системы Д оказывается значительно выше, чем точность расчета с помощью системы Б. Поэтому для расчета движения спутника по сильно вытянутым орбитам целесообразно использовать систему (6). При расчете орбит с малыми и умеренными эксцентриситетами использование системы Д нецелесообразно, так как в этом случае она уступает, например, системам Б и Г.

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

1. Карпов И. И., Платонов А. К. Ускорение численного интегрирования дифференциальных уравнений в небесной механике. Космические исследования, т. X, вып. 6, 1972.

2. Основы теории полета космических аппаратов. Под ред. Г. С. Нариманова и М. К. Тихонравова. М., „Машиностроение", 1972.

3. Эльясберг П. Е. Введение в теорию полета искусственных спутников Земли. М., „Наука", 1965.

4. С е м е н к о В. П. О д'вижении спутника в неинерциальной системе координат. Труды вторых чтений, посвященных разработке научного наследия Ф. А. Цандера. М., 1973.

5. С е м е н к о В. П. Исследование уравнений движения спутника в неинерциальной системе координат. „Ученые записки ЦАГИ", т. V, № 4, 1974.

6. Бэттин Р.. Г. Наведение в космосе. Пер. с англ. М., „Машиностроение", 1966.

ЛИТЕРАТУРА

Рукопись поступила 2111Х 1973 г.