ДОКЛАДЫ АКАДЕМИИ НАУК, 2012, том 445, № 1, с. 32-36
= МЕХАНИКА =
УДК 521.1 + 629.78
ОПТИМАЛЬНЫЕ ТРАЕКТОРИИ ПЕРЕЛЕТА КОСМИЧЕСКОГО АППАРАТА С МАЛОЙ ЭЛЕКТРОРЕАКТИВНОЙ ТЯГОЙ
К АСТЕРОИДУ АПОФИС © 2012 г. В. В. Ивашкин, И. В. Крылов
Представлено академиком Т.М. Энеевым 06.03.2012 г. Поступило 14.03.2012 г.
В настоящее время весьма актуальна задача организации исследовательской экспедиции к сближающемуся с Землей астероиду Апофис. 13 апреля 2029 г. Апофис пролетит вблизи нашей планеты на номинальном расстоянии ~38 тыс. км, вследствие чего подвергнется сильному гравитационному воздействию, которое значительно изменит его орбиту. В трубке рассеивания, построенной вокруг новой номинальной орбиты астероида, есть траектории, оканчивающиеся столкновением с Землей в 2036 г. [1]. Поэтому важно уточнить имеющуюся у нас информацию о движении и особенностях физического строения Апофиса, например, заблаговременно послав к нему космический аппарат (КА). В данной работе рассматривается вопрос определения энергетически оптимальной траектории полета к астероиду КА с малой электрореактивной тягой.
Анализируется схема экспедиции, которая включает в себя три основных этапа: геоцентрический этап разгона КА с помощью двигателя большой тяги, гелиоцентрический этап, на котором осуществляется перелет КА к астероиду с выравниванием их координат и скоростей при помощи двигателя малой тяги, и, наконец, этап торможения КА с выходом на орбиту искусственного спутника астероида.
На первом этапе (геоцентрического движения) предполагается использовать ракету-носитель "Союз-ФГ", которая обеспечит доставку КА массой т1 = 7130 кг на промежуточную орбиту высотой Н1 = 200 км. Последующий разгон аппарата может быть осуществлен с помощью блока "Фрегат" [2]. При этом считается, что набор необходимой геоцентрической энергии обеспечивается несколькими активными участками и гравитационные потери малы. Тогда масса аппарата
Институт прикладной математики им. М.В. Келдыша Российской Академии наук, Москва Московский государственный технический университет им. Н.Э. Баумана
т ) в момент начала гелиоцентрического движения (в импульсном приближении) определяется так:
m
(t0} = Щ eXP [W)
- m„
(1)
где Д V =
vl + ^
Hi
Rj
— величина разгонного
импульса, Wф — скорость реактивной струи "Фрегата" (с удельной тягой 326 с), Шф = 970 кг — "сухая" масса "Фрегата", V» — гиперболический избыток скорости КА на бесконечности, ц З — гравитационный параметр Земли, Rj — радиус промежуточной орбиты. При V» = 0 масса m (t0) = 1630 кг.
На втором этапе, гелиоцентрического движения при помощи двигателей малой тяги, рассматривается задача максимизации конечной массы КА. Уравнения движения КА в прямоугольной гелиоцентрической эклиптической системе координат OXYZ и имеют вид
Г = V,
(2)
V = g (r) + u,
где г и v — гелиоцентрические радиус-вектор и
r
вектор скорости КА, g(r) = -ц— — ускорение силы
r
тяжести, ц — гравитационный параметр Солнца,
F
r = | r| — расстояние от КА до центра Солнца, u =--
m
управляющее ускорение двигателей малой тяги, F — их тяга. При условии постоянства мощности двигателя масса КА в конце гелиоцентрического движения m(tf) тем больше, чем меньше интеграл
[3]:
J =
1
J u 2dt,
(3)
где и = |и|. Первоначально полагаем, что при разгоне КА на первом участке Ух = 0, и граничные условия второго участка заданы в виде
о
г(?о) = Го, v(to) = vo, (4)
r(tf) = Гf, v(tf) = v f, (5)
где параметры r0, v0 соответствуют орбите Земли, а rf, v f — орбите астероида. Полагая также на первом этапе анализа, что ограничений на ускорение нет (случай идеальной тяги ИТ), сформулируем следующую задачу.
Задача 1. Найти такое управление u = u(opf>, что траектория (2) удовлетворяет условиям (4), (5) и функционал (3) минимален.
Эффективность решения задачи 1 во многом зависит от выбора начального приближения u(0) для u(opt). В данной работе этот выбор осуществляется следующим образом. Первоначально траектория КА строится методом динамического программирования в априорно заданной области изменения параметров [4, 5], что позволяет находить экстремум функции, не делая никаких предположений о ее выпуклости. Поскольку решение задачи в точной постановке ввиду "проклятия размерности" потребовало бы проведения очень большого объема вычислений, здесь принимается допущение о близости орбиты перелета к эклиптике, что позволяет независимо рассматривать движения КА в плоскостях OXY и OXZ [5]. Последующее уточнение полученных результатов осуществляется методом "блуждающей" трубки [6]. Окончательно, искомая траектория строится в рамках полной пространственной модели (1) методом локального варьирования, который требует меньших по сравнению с методом "блуждающей" трубки вычислительных ресурсов, обладая при этом худшей устойчивостью к локальным экстремумам [7]. Затем система (1) линеаризуется относительно полученной векторной функции r(0)(t):
(6)
Г = V,
V = §(г(0)(?)) + С(Г(0)(0)(Г - Г(0)(0) + и,
где в (г) = —— матрица размерности 3 х 3. дг
Минимум функционала (3) для системы (6) при условиях (4), (5) отыскивается на основе принципа максимума Понтрягина [8]. Тогда сопряженные уравнения имеют вид
¥v =
у, = -GT (г(0) (t , а управление u находится из условия
-у, (0)
u
2 '
\T,T
(7)
(8)
При этом у = , уЦ)г
В результате решения линейной краевой задачи (4)—(8) регулярными методами определяется соответствующий вектор у(0) ), а по фор-
lg J, км2/с3 -4.00
-4.25
-4.50
-4.75
-5.00
-5.25
-5.50
-5.75
-6.00
2 3
(t0 - 2020), г
Рис. 1. Зависимости функционала J от момента начала гелиоцентрического участка, полученные на интервале 25.06.2019 г.—27.09.2022 г. (Г0 - 2020 г.) для различных значений продолжительности перелета в сутках.
муле (8) — начальное приближение и(0) для задачи 1. Следует отметить, что изложенный выше подход [5] является по сути обобщением широко известного метода транспортирующих траекторий [9] на случай произвольных значений угловой дальности перелета.
После того, как управление и(0) получено, задача 1 решается в точной постановке на основе использования принципа максимума Понтрягина. Система сопряженных уравнений имеет вид
Vv=-V,, V, = -GT (r)¥v •
(9)
Двухточечная нелинейная краевая задача (2), (4), (5), (8), (9), решается методом продолжения по параметру, уравнение которого
d = -A-1(р)а(¥ (0)(to)),
d т
(10)
где oT(у (to)) = [(v(tf) - vf )T, (r(tf) - rf )T] - вектор
невязок краевых условий, A(y(t0)) =
др(у(о)) _ Зу(^) .
= р (т). Интегрируя (10) по те [0; 1] при начальном условии р(т)| т=0=у (0)(?0), получаем при т = 1 искомый вектор начальных сопряженных пере-
m(tf), кг 1470
1460 -
1450
1440
1430 h
1420
1410
1400
1390
где l — единичный вектор (|l| = 1). Ставится
30 40
Í0-Í1, сут
Рис. 2. Зависимости конечной массы КА от момента начала гелиоцентрического участка для нулевых и оптимальных значений гиперболического избытка скорости Vх для момента ^ 10.03.2020 г.
менных у{opt)(t0) = р(т)| T=i, а также оптимальную
(opt) л
траекторию и управление u для задачи 1.
С помощью описанной выше процедуры на интервале дат старта Tx = [25.06.2019; 27.09.2022] для различных значений продолжительности перелета в сутках А/ g {185, 230, 275, 320, 365} были определены оптимальные траектории для задачи 1. На рис. 1 приведены зависимости функционала (3) от времени отлета с орбиты Земли и величины А/. Анализ полученных данных показывает монотонное уменьшение функционала с ростом А/. Поэтому в дальнейших исследованиях в качестве основного варианта принято А/ = 365 сут. При этом выделяется оптимальная дата отлета /0 = = 23.03.2020 г. Соответствующую оптимальную траекторию (и точку на рис. 1 для нее) далее обозначим символом Тр1. При оценке массовых характеристик КА полагаем, что мощность двигателя малой тяги в струе ^дв = 3.75 кВт, а скорость истечения We = 25 км/с [10]. Тогда для данной траектории конечная масса КА m(tf ) = 1409 кг.
Дальнейшее увеличение массы m(tf ) возможно
за счет сообщения КА при разгоне с промежуточной орбиты у Земли дополнительного импульса AV, что приведет к росту скорости Vœ. В этом случае начальные условия второго, гелиоцентрического, участка определятся так:
r(t о) = Го, v(to) = v 0 + VJ (11)
Задача 2. Найти такие u = u
(opt)
l
l(opt) и
Уо = У^0, что траектория (2) удовлетворяет условиям (11), (5) и функционал ) (конечная масса) максимален.
Решение задачи 2 находим на основе определения т (/0) после разгона у Земли и последующего решения задачи 1. С этой целью для начального момента /0 исследуем временной диапазон Т2 = = [10.03.2020; 19.04.2020], лежащий в окрестности Тр1. Продолжительность перелета А ^ между орбитами Земли и Апофиса считаем равной 365 сут. Для каждой выбранной даты старта из интервала Т2 задаются значения гиперболического избытка скорости из множества Уад е [0; 1.5] км/с, а также допустимые величины компонент вектора Ь Начальная масса КА определяется по (1), после чего решается задача 1 в исходной постановке. Конечная масса КА определяется по формуле
о)
m(tf )
2NдВ + m(t o)J
(12)
Оптимальность текущего вектора l оценивается путем проверки соотношения 1 {ор1) ТТ и{ор1) (0), которое вытекает из условий трансверсальности [10].
На рис. 2 представлено изменение конечной массы КА в окрестности оптимальной траектории Тр1 для Уад = 0 и Уда = У^рР). Как видно, максимум функционала, соответствующий У^ = = 0.95 км/с, смещается в точку 10 = 09.04.2020 г., т.е. наилучший момент старта сдвигается вправо по временной оси примерно на две недели. Эту траекторию, оптимальную для задачи 2, обозначим Тр2. Для нее на втором участке начальная масса т^0) = 1598 кг. Конечная масса КА т(^) = 1461 кг увеличивается на 52 кг, т.е. на ~3.7%.
Далее учтем, что на практике величина тяги двигателей В = |Р| всегда ограничена и, как правило, носит ступенчатый характер (случай кусочно-постоянной тяги — КПТ):
F е {0; Fm
N
F = 2 дв
^ max ^ • We
(13)
Определяя текущую массу аппарата из уравнения
dm ___F_
dt ~ We
(14) = 1),
ре
и полагая и = —, где e — единичный вектор ( т
сформулируем следующую задачу.
Задача 3.
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.