КОСМИЧЕСКИЕ ИССЛЕДОВАНИЯ, 2014, том 52, № 2, с. 113-124
УДК 521.1; 629.78
ОПТИМИЗАЦИЯ ТРАЕКТОРИЙ ПЕРЕЛЕТА КОСМИЧЕСКОГО АППАРАТА С БОЛЬШОЙ И МАЛОЙ ТЯГОЙ К АСТЕРОИДУ АПОФИС
© 2014 г. В. В. Ивашкин1, И. В. Крылов2
1 Институт прикладной математики им. М.В. Келдыша РАН, г. Москва 2 Московский государственный технический университет им. Н.Э. Баумана Ivashkin@keldysh.ru, krylov_i_v@mail.ru Поступила в редакцию 19.07.2012 г.
Исследована задача оптимизации перелета КА к астероиду Апофис. Проанализирована схема полета, которая включает в себя геоцентрический этап разгона КА с двигателем большой тяги, гелиоцентрический этап управления при помощи двигателя малой тяги и этап торможения с выходом на орбиту спутника астероида. При этом решается задача оптимального управления для случаев идеальной и кусочно-постоянной малой тяги, определяется оптимальная величина и направление гиперболической скорости КА на "бесконечности" при отлете от Земли. Траектории КА находятся на основе разработанного комплексного метода оптимизации, сочетающего метод динамического программирования на первом этапе анализа и принцип максимума Понтрягина на заключительном этапе, а также метод продолжения по параметру, Получены оценки конечной массы КА, а также полезной массы, которая может быть доставлена к астероиду с помощью ракеты-носителя Союз-Фрегат.
БО1: 10.7868/80023420613060022
ВВЕДЕНИЕ
К числу актуальных задач, стоящих в настоящее время перед космонавтикой, следует отнести организацию экспедиции к сближающемуся с Землей астероиду Апофис. Исследования показали, что 13.IV.2029 г. Апофис пролетит вблизи нашей планеты на номинальном расстоянии ~38 тыс. км, вследствие чего подвергнется сильному гравитационному воздействию, которое существенно изменит его орбиту. Трубка рассеивания, построенная вокруг новой расчетной орбиты астероида, содержит траектории, приводящие к столкновению с Землей, в первую очередь — в 2036 году ([1] и др.). Энергия подобного столкновения Ест оценивается примерно в 500 МТ ТНТ (для Тунгусского события Ест ~ 12 МТ ТНТ). Поэтому важно заблаговременно вывести к Апофису оснащенный специальной исследовательской аппаратурой космический аппарат, который поможет улучшить информацию об элементах движения и особенностях физического строения астероида.
В данной работе рассматривается задача формирования энергетически оптимальных траекторий для такой экспедиции. При этом исследуется схема перелета, включающая в себя три основных этапа. На первом этапе (геоцентрического движения) КА, находящийся на низкой (с высотой Н ~ ~ 200 км) промежуточной околоземной орбите, разгоняется при помощи двигателей большой тя-
ги (например, блока Фрегат [2]) до второй космической или гиперболической скорости. На втором этапе (гелиоцентрического движения) осуществляется выравнивание пространственных координат и скоростей КА и астероида при помощи электрореактивных двигателей малой тяги (ДМТ) [3]. Поскольку ДМТ обеспечивают высокие скорости истечения реактивной струи (более десяти километров в секунду), их применение особенно выгодно в условиях продолжительного гелиоцентрического перелета. Наконец, на третьем этапе (движения в сфере действия астероида) КА гасит свою скорость с параболической до круговой и становится искусственным спутником Апофиса.
1. ЭТАП ГЕОЦЕНТРИЧЕСКОГО ДВИЖЕНИЯ
Будем полагать, что вывод КА на низкую околоземную орбиту в нашем случае осуществляется при помощи ракеты-носителя Союз-ФГ, которая обеспечивает доставку массы т1 « 7130 кг на высоту 200 км, а для разгона КА до второй космической или гиперболической скорости используется блок Фрегат. Для этого блока удельная тяга двигателя ~326 с, скорость истечения Жф ~ 3.198 км/с. Будем также полагать, что набор необходимой геоцентрической энергии производится несколькими активными участками и гравитационные потери малы. Тогда масса аппарата т (?0) в момент
начала гелиоцентрического движения приближенно определится следующим образом:
т (/о) = тт ехр
г \ АУ
V ^Ф У
- тф,
(1)
где ДУ = , № + ^
— величина разгонного
импульса, тф = 970 кг — "сухая" масса Фрегата, У^ - гиперболическая скорость КА на "бесконечности", цЗ — гравитационный параметр Земли, Ят - радиус промежуточной орбиты. Если Уш = 0, то масса КА в момент начала гелиоцентрического участка т (/о) равна 1630 кг.
2. ЭТАП ГЕЛИОЦЕНТРИЧЕСКОГО ДВИЖЕНИЯ
2.1. Задача об идеальной тяге. На втором этапе движение КА рассматривается в гелиоцентрической эклиптической прямоугольной системе координат, ось ОХ которой направлена в точку весеннего равноденствия, ось OZ - на северный полюс эклиптики, ось ОУ дополняет систему до правой. Будем учитывать только силу тяготения Солнца и реактивную тягу ДМТ. Тогда уравнения движения КА имеют вид:
Г = V
V = g (г) + и,
(2)
= -р,
(3)
р = - секундный массовый расход;
р = = ти; и = |и|; Же - скорость истечения частиц в реактивной струе для ДМТ. Полагая, что мощность ДМТ
Nдв = в ^/2
(4)
постоянна на протяжении всего перелета, интегрируем (3) с учетом (4) [3]. Тогда масса в конце гелиоцентрического этапа:
) =
2N д
2^ + т(^
т( о),
(5)
где t0, - заданные значения времени начала и конца второго этапа перелета,
1 =
tf
Iи
Л.
(6)
Из (5) следует, что конечная масса КА тем больше, чем меньше величина интеграла (6). При этом граничные условия для данного этапа задаются в виде:
Г/о) = Го, v(tо) = V о + УД
Г(tf) = Г-, v(tf) = vf,
(7)
где гт = [гх,гу, г1 ] и Vт = \ух, vy, V. ] — гелиоцентрические радиус — вектор и вектор скорости КА (знак "т" здесь и далее означает транспонирование); g(г) = -ц (г/г3) — ускорение силы тяжести;
ц — гравитационный параметр Солнца; г = | г | — расстояние от КА до центра Солнца;
т РТ
и = [их, иу, и.,] = — — управляющее ускорение,
т
создаваемое ДМТ; Е — тяга ДМТ; т - текущая масса КА, определяемая уравнением
где 1 — единичный вектор (|1| = 1), параметры г0, v0 соответствуют радиусу-вектору и вектору скорости Земли в момент /о, а параметры , V^ - радиусу-вектору и вектору скорости астероида в момент tf. Полагая, что никаких ограничений на управляющее ускорение не накладывается, сформулируем следующую задачу об идеальной тяге (ИТ):
Задача 1. Найти такое управление и = и(ор1), что траектория (2) удовлетворяет условиям (7) при У« = 0, а функционал (6) минимален.
В работе рассматриваются три возможных способа решения задачи 1. Первый способ основывается на принципе оптимальности Беллмана (динамическое программирование), второй способ предполагает использование принципа максимума Понтрягина, третий, комплексный способ построен на сочетании первых двух. Далее кратко рассмотрим эти подходы.
2.2. Решение задачи об идеальной тяге методом динамического программирования. Будем пола-
(орЧ
гать, что процедура отыскания управления и методом динамического программирования включает в себя три основных стадии, которые выполняются последовательно. На первой стадии задача 1 решается приближенно. Для этого задаются априорные диапазоны изменения параметров траектории КА. В нашем случае: гх, гу е [-180.0 • 106; 180.0 • • 106] км, ^ е [-10.0 • 106; 10.0 • 106] км, V,, vy е [-40.0; 40.0] км/с и V. е [-5.0; 5.0] км/с. Принимается также допущение о близости орбиты КА к плоскости эклиптики. Тогда, система (2) преобразуются к виду:
ук = ук+2; ук+2 =
.1 — .2; 12
у2 + у22) н_
■у к + ч; (8)
л/(112 + ¿2уШ
.1 + щ;
(9)
о
где к = 1, 2; УТ = [[,У2,Уз,У4] = [г„Гу,Vх,V.],
гТ = [¿1, ¿2 ] = ]; dy.it) = Ух2(/> + У22(0; функции у1 (t> и у2^> — получаются в результате интегрирования (8) при задании функций ик (0. Далее минимизируется функционал:
Jy =
j
J М + u2)dt
min
(10)
для плоского движения (8). Затем полученные зависимости у1 (^ и у2^> подставляются в (9), и минимизируется функционал
Jz =
j
J u3
dt
min
(11)
N-1
Jy = X Jy(y(t'y(ti+1)) ^ min,
i=0 N-1
Jz = X Jz(z(ti), z(ti+1)) ^ min.
(12)
(13)
i=0
ют вид:
Jy (y (h), y(ti+i)) =
C2 + C22 - (CA + C2C4)h + (C3 + C4)h2
Jz (z(ti), z(t i+i)) =
C 2h2"
C52 - C5C6h + C6h-
h (14) 4,
(15)
где постоянные интегрирования Сь C2, .., C6 в (14) и (15) вычисляются по формулам:
Ci = Zi(yi(ti), Уз&), yi(ti+i), y3(ti+i), g*),
C2 = Zl(y2(ti), y4 (ti), y2(ti+l), y4 (ti+l), gy ),
C3 = z 2Ш), ys(ti), yi(ti+i), ys(ti+i)), (16) C4 = Z l(yi(tt), y4(ti), y2(ti+l), y4(ti+l)), C5 = Zl(zi(ti), z2(ti), zi(ti+i), z2(ti+i), gz), C6 = Z 2(zi(ti), z 2(ti), zi(ti+i), z 2(ti+i)).
Здесь
(17)
для системы (9).
Введем в пространствах у) и г) координатные сетки Qy и Qг, которые характеризуются шагами h, И у и Иг по временной и пространственным координатам, а также элементарные операции, которые каждой паре узлов из Qy и Qz ставят в соответствие оптимальное управление, переводящее у(^> в у(^+1> и г()1 > в г(^+1>. При этом i = [0, 1, 2, ..., N — 1] — индекс узла на временных осях сеток Qy и Qг, tN = tf [4]. Тогда (10) и (11) можно записать следующим образом:
Следуя [4], будем полагать, что гравитационное поле на интервале [7,, постоянно и характери-
т( Г + г/+1 11 + ti зуется осредненным вектором g I , -—1
= [ёх, ёу, ёг. ]. Тогда элементарные операции име-
re \ 12(Пз -П1) - 8n2h - 4n4h _
Zi(ni, П2, Пз, П4, g) = ——— , 2 -— - 2g;
h
у, \ 24(п3 -П1) - 12(n2 +n4)h
Z 2СП1, П2, Пз, П4) = ——— 3 ——.
h
В работе полагается, что сетки Qy и Qz — прямоугольные параллепипеды и содержат по 15 узлов вдоль любого измерения. Решение задачи о минимизации функционала (12) отыскивается следующим образом. Для каждого момента ti+1 на временной оси Qy и каждого текущего узла P(tt+1) из Qy(ti+1) проверяются траектории, проходящие через все узлы сетки Qy(tt) и выбранный узел P(ti+1). Вариант, удовлетворяющий реккурентно-му соотношению Беллмана [5]:
S(y(ti+i)) = min [S(y(tf)) + Jy(y(tf), yfc+i))], i > 0;
y(ti)eQy(ti) (18)
S(y(to)) = 0,
где S(y(ti)), S(y(ti+1)) — значения (12) вдоль оптимальной траектории, считается оптимальным для P(ti+1) и запоминается. Таким образом, для сетки Qy строится последовательность а y перехода из единственного начального узла P(t0) в единственный конечный узел P(tN). Решение задачи о минимизации (13) (последовательность аz)
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.