Интегрирование пассивного участка траектории на основе сочетания метода KS-регуляризации с методом Энке Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКАЯ ФИЗИКА, МАТЕМАТИЧЕСКОЕ

МОДЕЛИРОВАНИЕ

MATHEMATICAL PHYSICS. MATHEMATICAL MODELING

УДК 629.78

DOI: 10.18413/2075-4639-2018-50-1-88-96

ИНТЕГРИРОВАНИЕ ПАССИВНОГО УЧАСТКА ТРАЕКТОРИИ НА ОСНОВЕ СОЧЕТАНИЯ МЕТОДА KS-РЕГУЛЯРИЗАЦИИ С МЕТОДОМ ЭНКЕ

INACTIVE LEG INTEGRATION ON BASIS OF KS-REGULARIZATION IN COMBINATION WITH ENKE'S METHOD

1А.А. Ляпин, 2Б.М. Тюлькин A.A. Lyapin, B.M. Tyulkin

1 Челябинский государственный университет, Миасский филиал, Россия, 456331, Челябинская обл., г. Миасс, пос. Строителей, ул. Керченская, 1. 'Chelyabinsk State University, Miass Branch, Russia, 456313, Chelybinskaia region,

Miass, 1,Kerchenskaya street

2Челябинский государственный университет, Миасский филиал, Россия, 456313, Челябинская обл., г. Миасс, пос. Строителей, ул. Керченская, 1.

2Chelyabinsk State University, Miass Branch, Russia, 456313, Chelybinskaia region,

Miass, 1,Kerchenskaya street

E-mail: lyapin-sasha@mail.ru, tulkin74@rambler.ru

Аннотация

В статье рассмотрен способ представления уравнения движения летательного аппарата на пассивном участке траектории с использованием преобразования Кустаанхеймо-Штифеля в сочетании с методом интегрирования уравнений в отклонениях от аналитически вычисляемой опорной орбиты, предложенным И.Ф.Энке. Приведена сравнительная характеристика исследуемого метода с классическими методами представления уравнений движения. Результаты вычислительных экспериментов показывают, что алгоритм Нистрема интегрирования уравнений движения, представленных в форме Энке в KS-пространстве, существенно эффективнее алгоритмов, основанных на классическом представлении уравнений движения.

Abstract

The article suggests the representation's method of the re-entry vehicle's equations motion at inactive leg with using Kustaanheimo-Stiefel's conversion in combination with Enke's method. Enke's method bases on the integral method in deviations from analytically-calculation orbit reference. I suggest comparative characteristic of researching method with classic methods of motion equations representation. Outcome of calculation-experiments are showed that Nistrem's integral algorithm of equatinos motion in the Enke's form in KS-dimensional space more effective classic methods of motion equations representation.

Ключевые слова: численные методы, пассивный участок, метод Энке, преобразование Кустаанхеймо-Штифеля.

Keywords: calculus of approximations, inactive leg, Enke's method, Kustaanheimo-Stiefel's conversion.

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

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

В текущей работе рассматривается задача интегрирования пассивного участка траектории движения ЛА. Для решения поставленных вопросов предлагается использовать преобразование Кустаанхеймо-Штифеля в сочетании с методом интегрирования уравнений в отклонениях от аналитически вычисляемой опорной орбиты, предложенным И.Ф.Энке.

Уравнение движения ЛА в центральном гравитационном поле Земли имеет следующий вид [1,3]:

- /М^

?х = - — ?х + <Э (1)

г/

В уравнении (1) силу @, действующую на единицу массы ЛА, можно рассматривать как возмущающую силу, по отношению к невозмущенному кеплеровскому движению ЛА в центральном гравитационном поле Земли. Основные свойства этой силы, с точки зрения нашей задачи, состоят в том, что она имеет потенциал, который в инерциальном пространстве определяется только положением ЛА и не зависит от времени и скорости полета.

Уравнение (1) теоретически имеет две особенности: сингулярность и отсутствие устойчивости по Ляпунову. Исследования, выполненные в 1970-1980-е годы [2], показали, что при прогнозе движения искусственных спутников Земли на точность численного решения уравнения (1) оказывают влияние неустойчивость уравнения и неравномерность изменения правой части по траектории полета, а также неустойчивость метода интегрирования, которые препятствуют увеличению шага интегрирования. В этот же период были предложены методы, устраняющие указанные особенности. Одним из таких методов является метод регуляризации уравнения, предложенный Кустаанхеймо и Штифелем [2, 4].

Преобразование Кустаанхеймо-Штифеля (^-преобразование) [2] позволяет регу-ляризировать уравнение (1), а также его решение. КБ-преобразование задает отображение физического пространства (х) в четырехмерное ^-пространство (—) и совместно с временным преобразованием ^ = |х|^ = гйБ, осуществляющим переход к новой независимой переменной ^ , называемой фиктивным временем, представляется в виде

х = ¿(—)—, где ¿(—) - матрица КБ-преобразования,

¿(—) =

—2 - Й3 Й4

Й2 —х - Й4 - Й3

—3 —4 —1 —2

—4 - Й3 Й2 ——1

(2)

(3)

Дифференциальные уравнения для и, I, к в функции от ^ имеют вид

Л 1 д ^ _

4 ЗУ (4)

Л'=-2|—|2—V = |—|2. д

Результаты вычислительных экспериментов по исследованию эффективности алгоритмов численного прогнозирования движения искусственных спутников земли, приведенные в [5], показывают, что алгоритмы, основанные на представлении уравнений движения в КБ-переменных, существенно превосходят алгоритмы, основанные на классическом представлении уравнений движения в прямоугольных координатах, как по точности, так и по времени расчета.

Рассматриваемый в настоящей работе класс орбит (траекторий свободного полета ЛА) отличается от орбит, приведенных в [5]. Большая часть этих орбит незамкнута, т.е. они пересекаются с поверхностью Земли. Эксцентриситеты орбит е > 0.2 , т.е. они относятся к классу средне- и высокоэллиптических орбит. Высота апогея лежит в интервале [800,7000] км. Интервал интегрирования участка свободного полета не превышает 4500 с.

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

1. Уравнения движения в КБ-переменных имеют порядок 9 вместо 6 в прямоугольных координатах.

2. В процессе численного интегрирования системы (4) на каждом шаге требуется проведение КБ-преобразования для определения вектора положения х в физическом пространстве, чтобы вычислить частные производные потенциала У24 по параметрам и1, и2, и3, и4 и расстояние до граничной поверхности интегрирования.

При решении системы уравнений (4) будем использовать метод численного интегрирования уравнений «в отклонениях от аналитически вычисляемой опорной орбиты», который также называется интегрированием «по оскулирующей орбите». Такой подход для повышения точности прогноза орбит впервые реализовал Иоганн Франц Энке (1857) [5].

Для вывода дифференциальных уравнений метода Энке в КБ-переменных вводится новая независимая переменная Е = - обобщенная эксцентрическая аномалия, для которой выполняются следующие соотношения [2]:

й й

йБ йЕ' (5)

й2 _ й2 йш й

-= 4ш2--+ 4ш--.

йЕ2 йЕ йЕ

Здесь, ш = - частота системы.

Вместо физического времени I в системе (5) предлагается [2] использовать временной элемент т, связанный с физическим временем соотношением

т = Ь+--- и •и =1—и • и . (6)

2ш2 П у '

При отсутствии возмущающего потенциала система уравнений принимает

вид

й2Щ 1_>

-¿Ё2Г + 1ик = 0, (7)

йтк д йшк

, йЕ йЕ

Решение системы (7) в КБ-пространстве представляется в виде

'-» -» Е Е Р

ик = аксо5- + РцБЫ-, шк = ш0, тк=—^Е + т0 2 2

Е , Е йик (8)

<*к = С05-^ик - 2зЫ-

2 2 йЬ

-> Е ^ Е йик

Рк = ып-• ик + 2С08-• —

Здесь ¿^ и / - векторные константы, определяющие положение перицентра и орбитальной плоскости в ^-пространстве. Они определяются из начальных данных.

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

d<5r

d25U 1 ^ 1 Г /

--з) + 2Г^24

1 , ra^24> '-u +

4 du

\ 1 dw dU

/ w d^ dK24 2 dût

+

16w3

ut

ut

dût w2 d^ d^

(9)

5U = ut — Ufc, = w — <5t = t — rfc

dw

Г d^24

= 0

Численное интегрирование полученной системы уравнений проводится методом Нистрема 4-го порядка.

Методология оценки эффективности

Эффективность метода Нистрема интегрирования уравнений движения (9) определяется в сравнении со следующими методами [6]:

- методом Рунге-Кутта 4-го порядка в физическом пространстве,

- методом Нистрема 4-го порядка в физическом пространстве. Сравнительная характеристика методов интегрирования уравнений движения ЛА

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

В качестве эталонного решения используются результаты интегрирования методом Рунге-Кутта 4-го порядка точности при шаге по времени 10 сек. Практика показывает, что точность расчета параметров траектории при таком шаге не превышает 1м. Будем считать, что необходимая точность достигнута, если длина вектора отклонения полученного решения от эталонного Аг<6м (рисунок 1).

А г

Рис 1. Отклонение полученного решения от эталонного (rk - полученное решение,

Гэ - эталонное решение) Fig. 1. Deviation of computed solution from reference one (rk - computed solution, Гэ - reference solution)

Результаты расчетов этими методами будем анализировать по следующей схеме:

1. Анализ влияния формы траектории на объем вычислений.

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

2. Анализ вычислительных затрат каждого метода.

2.1. Оценим временные затраты на расчет одного шага каждым из методов, с учетом следующих особенностей. Во-первых, из-за малости величины временных затрат на 1 шаг, программные средства не позволяют получить их точную оценку. Для устранения данной проблемы возьмем «усредненную» оценку времени, полученную следующим образом: зациклим время расчета основного участка траектории 1000 раз (обозначим полученное

т

время Т), тогда «усредненное» время 1 шага Т[ = где К-количество шагов интегри-

рования.

Во-вторых, полученное «усредненное» время будет действительно для той машины, на которой проводились расчеты. Для получения универсальной оценки введем условную единицу. Полученные относительно этой единицы временные затраты будут сохраняться при переходе как на другую машину, так и на другую операционную систему. За условную единицу (у.е.) принимается время, затрачиваемое на расчет одного шага методом Рунге-Кутта 4-го порядка точности системы дифференциальных уравнений движения ЛА в инерциальной СК.

2.2. Сопоставим полученные временные затраты количеству шагов, необходимому для достижения заданной точности, полученному в п.1.

Для реализации вышеописанной схемы использовалось специальное программное средство отработки и тестирования задач расчета и контроля полетного задания [7 - 11].

Результаты

Для проведения сравнительной характеристики составлен список входных элементов различных форм траекторий (всего было проведено более 10000 расчетов). Для наглядного представления результатов вычислений в таблице 1 представлены некоторые характерные траектории. Основным критерием выбора этих начальных данных является различие по эксцентриситету и параметру конического сечения, величины этих параметров оказывают основное влияние на форму и длину траектории (таблица 2).

Введем следующие обозначения: г= [тх,гу,г2 ] - вектор положения, М = [Цс, ~] — вектор скорости, Ьнач - начальное время, *;пол — время полета, Н5 — высота конечной точки траектории над ОЗЭ, е - эксцентриситет, р - параметр конического сечения.

Таблица 1. Начальные данные для контрольных траекторий Table 1. Initial data for test trajectories

Траектория Гх ry rz Vx Vy Vz ^ач

№1 670572 7472058 -1380119 -6536,918 2566,391 301,875 531,64

№2 675869 -1380119 7472058 -4317,732 2857,410 -190,211 531,64

№3 670572 7472058 -1380119 -3317,732 2857,410 -190,211 531,64

№4 675869 -1380119 7472058 -1317,732 1857,410 -190,211 531,64

№5 -670572 7472058 1380119 -3317,732 4527,320 -2190,211 531,64

№6 680572 7359860 -2356412 -5536,918 2566,391 -3301,875 531,64

Таблица 2. Основные параметры движения Table 2. Main moving parameters

Траектория e tncra p

№1 0,273207 3829,1183 6693993,26

№2 0,519602 582,6074 3746899,41

№3 0,775737 1439,4918 1860869,05

№4 0,908264 523,9485 702831,79

№5 0,865294 3082,5653 1468345,668

№6 0,429227 4091,3472 5953743,441

Эталонные решения приведены в таблице 3.

Таблица 3. Эталонные решения Table 3. Test solutions

Траектория Tx ry Tz Vx Vy Vz

№1 -2168325,57 -5843348,55 1266567,24 8284,26 -994,238 -771,691

№2 -1798674,55 415475,73 6105051,88 -3861,36 3131,13 -4759,04

№3 -3372209,13 5359200,63 -627915,6 -1148,57 -6104,73 1181,78

№4 -71814,60 -268393,16 6372048,35 -1476,97 2317,79 -4264,65

№5 -4678431,43 3407362,56 -2665443,09 2749,67 -6666,72 2231,36

№6 -4289589,90 -4701846,40 -154063,95 7210,16 -2002,97 3824,43

Введем следующие обозначения: RK - метод Рунге-Кутта 4-го порядка точности, Nst - метод Нистема в физическом пространстве, ^ - метод Нистема в ^-пространстве. В таблице 4 представлено влияние формы траектории на количество шагов интегрирования для достижения заданной точности.

Таблица 4 Оценка влияния формы траектории на расчет Table 4. Valuation of dependence of trajectory form on calculations

Траектория e t P N(RK) N(Nst) N(KS)

№1 0,273207 3829,1183 6693993,26 51 56 6

№2 0,519602 582,6074 3746899,41 5 5 1

№3 0,775737 1439,4918 1860869,05 13 13 3

№4 0,908264 523,9485 702831,79 4 4 1

№5 0,865294 3082,5653 1468345,67 19 20 5

№6 0,429227 4091,3472 5953743,44 29 32 9

Анализируя таблицу 4 можно сделать следующие выводы:

1) При классическом представлении уравнений движения форма траектории влияет на объем требуемых вычислений. Для траекторий, близких к круговым, у которых величина эксцентриситета близка к нулю, количество шагов интегрирования, требуемое для обеспечения заданной точности расчета, больше, чем для траекторий, у которых величина эксцентриситета близка к единице (например, траектории №1 и №5).

2) При представлении уравнений движения ЛА в KS-переменных форма траектории практически не влияет на объем вычислений, он зависит от продолжительности полета, т.е. от размера траектории.

Исследуем временные затраты различных алгоритмов для достижения заданной точности. В результате расчетов выявлено, что временные затраты на 1 шаг интегрирования методом Нистрема в физическом пространстве Ti=0.66 у.е, а в KS-пространстве Ti=2.02 у.е. Составим таблицу для сравнения общих временных затрат для всех рассмотренных методов (таблица 5).

Таблица 5. Общие временные затраты Table 5. Overall time expenses

Траектория T(RK) (у.е) T(Nst) (у.е) T(KS) (у.е)

№1 51 36,96 12,12

№2 5 3,3 2,02

№3 13 8,58 6,06

№4 4 2,64 2,02

№5 19 13,2 10,1

№6 29 21,12 18,18

Из таблицы следует вывод: несмотря на то что объем вычислений в KS-пространстве на каждом шаге больше, чем при классическом представлении уравнений, общий объем вычислений значительно ниже.

Заключение

Результаты вычислительных экспериментов, показывают, что алгоритм Нистема для уравнений движения представленных в форме Энке в KS-пространстве, существенно эффективнее алгоритмов, основанных на классическом представлении уравнений движения. Во-первых, несмотря на то, что объем вычислений на каждом шаге интегрирования больше, чем у классических методов, общий объем вычислений для получения заданной точности будет меньше. Во-вторых, форма траектории не влияет на точность расчетов, что делает метод Нистрема в KS-пространстве универсальным.

Список литературы References

1. Аппазов Р.Ф., Сытин О.Г. 1987. Методы проектирования траекторий носителей и спутников Земли. Издательство «Наука». Главная редакция физико-математической литературы, 435.

Appazov R. F., Sytin O.G. 1987. Methods for designing the trajectories of carriers and Earth satellites. "Nauka", 435.

2. Штифель Е., Шейфеле Г. 1975. Линейная и регулярная небесная механика. -М.:Наука, 30.

Shtifel E., Sheifel G. 1975. Linear and regular celestial mechanics. Nauka, 304.

3. Эскобал П. 1970. Методы определения орбит, М: Мир, 472.

Escobal P. 1970. Methods for determining orbits. MIR, 472 .

4. Иванюхин А.В. 2014. Оптимизация траектории космического аппарата с идеально регулируемым двигателем в переменных Кустаанхеймо-Штифеля, «Труды МАИ», выпуск 75.

Ivanyukhin A.V. 2014. Optimization of the trajectory of a spacecraft with an ideally adjustable engine in Kustaanheimo-Stiefel variables, "Proceedings of the MAI", 75.

5. Бордовицына Т.В. 1984. Современные численные методы в задачах небесной механики. - М.: Наука Главная редакция физико-математической литературы, 136.

Bordovitsyna T.V. 1984. Modern numerical methods in the problems of celestial mechanics, Moscow, Nauka, 136.

6. Хайрер Э., Нерсетт С., Ваннер Г. 1990. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи - М.: Мир.

Hirer E., Nersett S., Vunner G. 1990. Solution of regular differential equations. Lenient tasks. Moscow: MIR.

7. Ляпин А.А. Методология тестирования ПО РКПЗ. Материалы XX Юбилейной Международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС'2017), 24-31 мая 2017 г., Алушта. - М.: Изд-во МАИ, 2017. 154-156 с. ISBN 978-5-4316-0401-0.

Lyapin A.A. Methodology of CCFT software testing. Materials of XX Anniversary International Conference covers the computing mechanics and state-of-the-art applications systems (CMSAS' 2017), May 24-31th, 2017, Alushta. - M: MAI, 2017. 154-156 p.

8. Ляпин А.А. Разработка программного комплекса для баллистического проектирования: построение структурной схемы/ Наука XXI века: проблемы, поиски, решения: материалы XL научно-практической конференции, посвященной 40-летию Челябинского государственного университета и 20-летию Миасского филиала ЧелГУ (Миасс, 22 апреля 2016г.) / под ред. канд. филол. н., доц. А. Г. Бент; Миасский филиал ФГБОУ ВПО «ЧелГУ». - Миасс: Геотур, 2016. -382 с. - ISBN 978-5-89204-262-8.

Lyapin A.A. Development of a software package for ballistic design: construction of a structural diagram. Science of the XXI century: problems, searches, solutions: materials of the XL scientific-practical conference dedicated to the 40th anniversary of the Chelyabinsk State University and the 20th anniversary of the Miass branch of the CSU (Miass, April 22, 2016), Geotur, 2016. - 382 p.

9. Ляпин А.А. Автоматизация инженерных процессов отработки и тестирования программных модулей расчета и контроля полетного задания. Сборник материалов Всероссийской молодежной научно-практической конференции "Орбита молодежи" и перспективы развития российской космонавтики, 2017, Томск, издательство томского политехнического университета, 95-96 с. ISBN 978-5-4387-0776-9.

Lyapin A.A. Automation of engineering processes for testing of programm modules for calculating and controlling a flight task. Collection of materials of the All-Russian Youth Scientific and Practical Conference "Youth Orbit" and prospects for the development of Russian cosmonautics, 2017, Tomsk, Tomsk Polytechnic University, 95-96 p.

10. Ляпин А.А. 2017. Тестирование программных модулей расчета и контроля полетного задания на основе построения приоритетов исходных данных. Труды МАИ, № 97.

Lyapin A.A. Programming modules calculation and control of the flight task testing based on prioritization of source data. "Proceedings of the MAI", vol. 97, 2017.

11. Программа BTTesting: программа для ЭВМ №2017664117 / А.А. Ляпин; правообладатель: Акционерное Общество «Государственный ракетный центр имени академика В.П. Макеева» (RU). Зарегистрировано 18.12.2017.

Lyapin A.A. BTTesting: programm for IBM. Patent RF 96875, no. 2017664117; decl. 18.12.2017.

12. Сихарулидзе Ю.Г. 1982. Баллистика летательных аппаратов. М.: Наука. Главная редакция физико-математической литературы, 352.

Sikharulidze Y.G. Ballistics of aircraft. - Moscow: Nauka, 1982, 352 p.