Космические исследования, 2019, T. 57, № 5, стр. 386-400

ВВЕДЕНИЕ

Одним из вариантов выведения КА на геостационарную орбиту (ГСО) является комбинированная схема выведения с помощью разгонного блока (РБ) и электроракетной двигательной установки (ЭРДУ). Данная схема обладает преимуществами и недостатками как схемы выведения только с помощью РБ, так и полностью электроракетной схемы. Анализу различных аспектов комбинированной схемы посвящено множество работ [1–5]. Одной из проблем, возникающих во время электроракетной фазы выведения КА на ГСО, является накопление радиационных повреждений в фотоэлектрических преобразователях (ФЭП) СБ при прохождении КА радиационных поясов Земли (РПЗ). Радиационные повреждения ФЭП приводят к существенной деградации его основных характеристик, в том числе и к снижению выходной электрической мощности, необходимой для функционирования бортовых систем, включая ЭРДУ во время выведения. Ситуация осложняется неблагоприятными радиационными условиями на ГСО и длительными сроками активного существования (САС) на целевой орбите, характерными для связных КА, являющихся большинством среди выводимых на ГСО аппаратов. Кроме того, наблюдается эффект “лавинной” деградации характеристик ФЭП [6, 7], когда после некоторого количества повреждений наблюдается экспоненциальное падание характеристик ФЭП. Все эти факторы приводят к необходимости учета радиационной деградации электрической мощности СБ при проектировании КА. Для парирования падения мощности СБ и получения гарантированной выходной мощности на конец САС необходимо при проектировании закладывать избыточную мощность СБ, что приводит к увеличению площади СБ, сложности их конструкции и повышению стоимости КА.

В данной работе рассматривается способ снижения радиационной нагрузки на ФЭП и бортовые системы в целом, связанный с изменением формы траектории выведения. Изменение формы траектории перелета может привести к снижению суммарного количества заряженных частиц (флюенса), поглощенных ФЭП, благодаря тому, что распределение потоков заряженных частиц РПЗ в околоземном пространстве неоднородно [7]. Следовательно, за одно и тоже время пребывания в РПЗ на поверхность ФЭП могут упасть различное количество заряженных частиц в зависимости от того, где внутри РПЗ находится КА.

Данная работа является продолжением работ [8, 9], где рассматривалась оптимизация траекторий комбинированного выведения КА с СБ на основе кремниевых ФЭП. Предполагается, что в данной работе КА оснащен СБ с перспективными трехкаскадными фотоэлектрическими преобразователями ФЭП, обладающими на данный момент самым высоком коэффициентом полезного действия [10].

Помимо задач, связанных с комбинированным выведением, проблема большой радиационной нагрузки может возникнуть при электроракетном выведении полезного груза на ГСО с помощью многоразового межорбитального буксира с ядерной энергетической установкой. В этой задаче перелет начинается с круговой орбиты высотой 800 км, что приводит к еще большей радиационной нагрузке на бортовые системы буксира и полезной нагрузки. Задача снижения поглощенной в РПЗ дозы радиации с помощью изменения формы траектории перелета межорбитального буксира рассмотрена в работе [11].

ПОСТАНОВКА ЗАДАЧИ

Рассматривается следующая схема выведения КА на ГСО. На первом этапе ракета-носитель типа “Союз” выводит КА с разгонным блоком, аналогичным по характеристикам РБ “Фрегат”, на низкую опорную орбиту. Затем РБ выполняет доставку КА на промежуточную орбиту и отделяется от КА. На заключительном этапе КА с помощью ЭРДУ переводится по многовитковой траектории на ГСО.

Будем считать, что КА на этапе выведения с помощью ЭРДУ движется только под действием силы притяжения Земли и силы тяги ЭРДУ. В разложении гравитационного потенциала притяжения Земли будем учитывать только главный член и вторую зональную гармонику. Всеми остальными возмущающими силами пренебрежем. Электрическая мощность, необходимая для работы ЭРДУ, генерируется только бортовыми СБ на основе трехкаскадных ФЭП. Предполагается, что суммарная электрическая мощность СБ выбрана с избытком так, что на конец САС КА (15 лет без учета выведения) генерируемая мощность будет достаточной для функционирования целевой аппаратуры, служебных систем и выполнения КА целевой задачи в целом. Также будем считать, что максимальная тяга ЭРДУ $P$ и удельный импульс ${{I}_{{sp}}}$ в течение всего этапа электроракетного выведения постоянны. Пусть на борту рассматриваемого КА в качестве маршевой ЭРДУ установлены два двигателя СПД-140Д [12] суммарной максимальной тягой $P = 0.58$ Н. Удельный импульс ЭРДУ принимаем равным ${{I}_{{sp}}} = 1780$ с. Кроме того, примем допущение, что на возможную ориентацию вектора тяги ЭРДУ не наложено ограничений, и при этом соответствующие развороты КА будут считаться не требующими затрат рабочего тела. Наконец, будем считать, что деградация СБ происходит только под воздействием на ФЭП протонов и электронов РПЗ, воздействием других источников космической радиации пренебрежем.

Пусть на начало этапа выведения с помощью ЭРДУ КА находится в перицентре промежуточной орбиты и имеет некоторую массу ${{m}_{0}},$ высота перицентра промежуточной орбиты равна ${{h}_{{\pi ,0}}},$ высота апоцентра ${{h}_{{\alpha ,0}}},$ наклонение – ${{i}_{0}} > 0.$ Аргумент перигея и долгота восходящего узла промежуточной орбиты равны ${{\omega }_{{\pi ,0}}}$ и ${{\Omega }_{0}}$, соответственно. Целевой орбитой КА будем считать ГСО: высота ${{h}_{T}} = 35{\kern 1pt} 793$ км, нулевые эксцентриситет и наклонение ${{i}_{T}} = 0.$ Угловое положение КА на конечной орбите не фиксируется. Для заданных фиксированных значений ${{h}_{{\pi ,0}}},$ ${{h}_{{\alpha ,0}}},$ ${{i}_{0}},$ ${{\omega }_{{\pi ,0}}},$ ${{\Omega }_{0}}$ и ${{m}_{0}}$ найдем такие траектории перелета c помощью ЭРДУ, чтобы падение максимальной выходной электрической мощности СБ КА вследствие влияния космической радиации за время перелета с помощью ЭРДУ было наименьшим. Поскольку в течение САС мощность СБ только уменьшается и при проектировании КА необходимо гарантировать достаточный уровень энергоприхода на конец САС, то в целевой функционал для большей наглядности будем включать деградацию СБ за время нахождения КА на ГСО вплоть до конца САС. Тогда рассматриваемую задачу можно записать следующим образом $\Delta N = {{N}_{0}} - {{N}_{{{\text{С А С }}}}} \to \min ,$ где ${{N}_{0}}$ – мощность СБ на начало этапа выведения с помощью ЭРДУ, а ${{N}_{{{\text{С А С }}}}}$ – мощность СБ на конец САС. При фиксированной ${{N}_{{{\text{С А С }}}}}$ эта задача эквивалентна задаче максимизации относительной мощности СБ на конец САС

Таким образом, далее будем рассматривать задачу поиска траекторий электроракетного перелета с максимальным значением относительной мощности ${{\beta }_{{{\text{С А С }}}}}$ СБ на конец САС для каждого выбранного набора значений ${{h}_{{\pi ,0}}},$ ${{h}_{{\alpha ,0}}},$ ${{i}_{0}},$ ${{\omega }_{{\pi ,0}}},$ ${{\Omega }_{0}}$ и ${{m}_{0}}.$

УРАВНЕНИЯ ДВИЖЕНИЯ

Введем орбитальную систему координат $OXYZ$ следующим образом. Начало координат $O$ расположено в центре масс КА. Радиальная ось $OX$ направлена из начала координат $O$ в направлении от центра Земли, бинормальная ось $OZ$ сонаправлена с вектором орбитального кинетического момента КА. Трансверсальная ось $OY$ лежит в плоскости оскулирующей орбиты, перпендикулярна оси $OX$ и направлена в сторону движения КА по орбите. Чтобы избежать особенностей в окрестности нулевых эксцентриситета и наклонения будем использовать для описания движения КА на орбите уравнения в равноденственных элементах:

где равноденственные элементы $h,$ ${{e}_{x}},$ ${{e}_{y}},$ ${{i}_{x}},$ ${{i}_{y}},$ $F$ выражаются через фокальный параметр $p$ оскулирующей орбиты, эксцентриситет $e,$ аргумент перицентра $\omega ,$ наклонение $i,$ долготу восходящего узла $\Omega ,$ истинную аномалию $\nu $ следующим образом: $h = \sqrt {{p \mathord{\left/ {\vphantom {p \mu }} \right. \kern-0em} \mu }} ,$ ${{e}_{x}} = e\cos \left( {\Omega + \omega } \right),$ ${{e}_{y}} = e\,{\text{sin}}\left( {\Omega + \omega } \right),$ ${{i}_{x}} = \operatorname{tg} \left( {{i \mathord{\left/ {\vphantom {i 2}} \right. \kern-0em} 2}} \right)\cos \Omega ,$ ${{i}_{y}} = \operatorname{tg} \left( {{i \mathord{\left/ {\vphantom {i 2}} \right. \kern-0em} 2}} \right)\sin \Omega ,$ $F = \nu + \omega + \Omega ,$ ${{a}_{X}},$ ${{a}_{Y}},$ ${{a}_{Z}}$ – компоненты возмущающего ускорения, действующего на КА на орбите, в проекции на оси орбитальной системы координат, $m$ – текущая масса КА, $w = {{I}_{{sp}}}g$ – скорость истечения рабочего тела ЭРДУ, $g = 9.80665$ м/с² – ускорение свободного падения, $\mu = {\text{3}}{\text{.98600436}} \times {\text{1}}{{{\text{0}}}^{{14}}}$ м³/с² – гравитационный параметр Земли. В уравнениях движения также используются обозначения: $\xi = 1 + {{e}_{x}}\cos F + {{e}_{y}}\sin F,$ $\eta = {{i}_{x}}\sin F - {{i}_{y}}\cos F$ и $\tilde {\varphi } = 1 + i_{x}^{2} + i_{y}^{2}.$

Будем задавать направление вектора реактивного ускорения, создаваемого ЭРДУ, с помощью угла тангажа и угла рысканья. Углом тангажа $\vartheta $ будем называть угол между проекцией вектора тяги на плоскость $OXY$ оскулирующей орбиты КА и трансверсальным направлением $OY.$ Углом рысканья $\psi $ будем называть угол между вектором тяги и плоскостью $OXY$ оскулирующей орбиты КА. Тогда проекции возмущающего ускорения от работы ЭРДУ и второй зональной гармоники разложения гравитационного потенциала Земли на оси орбитальной системы координат будут иметь следующий вид:

где $\delta $ – функция включения ЭРДУ (при $\delta = 0$ двигательная установка выключена, при $\delta = 1$ – включена) и введены следующие обозначения а также ${{J}_{2}} = {\text{1}}{\text{.08263}} \times {\text{1}}{{{\text{0}}}^{{ - 3}}}$ – коэффициент второй зональной гармоники и ${{R}_{{\text{э }}}} = {\text{6378}}{\text{.14}}$ км – экваториальный радиус Земли.

ЦЕЛЕВОЙ ФУНКЦИОНАЛ

Максимальная выходная электрическая мощность СБ КА на орбите зависит от суммарного количества заряженных частиц (флюенса) космической радиации, упавших на единицу площади чувствительной поверхности ФЭП СБ за время от начала перелета до текущего момента. Для расчета флюенсов электронов и протонов РПЗ различных энергий в конкретной точке околоземного пространства будем использовать модели AE8 MAX и AP8 MAX [13–15], реализованные в программном комплексе IRENE (https://www.vdl.afrl. af.mil/programs/ae9ap9/downloads.php). Согласно рекомендациям работы [16] для расчета магнитных координат при использовании модели AE8 MAX в качестве моделей магнитного поля внутренних источников Земли использовалась модель Jensen и Cain [17]. Для модели AP8 MAX использовалась экстраполированная на эпоху 1970 модель магнитного поля внутренних источников GSFC 12/66 [18]. Магнитное поле внешних токов магнитосферы Земли не учитывалось.

Падение максимальной электрической мощности СБ рассчитывалось с помощью метода эквивалентного потока [6]. Его суть состоит в переводе всенаправленного потока протонов и электронов через набор коэффициентов в моноэнергетический нормально падающий поток электронов с энергией 1 МэВ, который создаст то же количество однотипных повреждений ФЭП в единицу времени, что и исходные потоки. Таким образом, уровень повреждений, нанесенных космической радиацией ФЭП, можно измерять в эквивалентном флюенсе нормально падающих на ФЭП электронов с энергией 1 МэВ. Зная величину эквивалентного флюенса электронов с энергией 1 МэВ можно по экспериментальным данным оценить величину деградации различных электрических характеристик ФЭП, в том числе и снижение максимальной электрической мощности ФЭП. Для расчета эквивалентного флюенса по спектрам падающих на ФЭП электронов и протонов РПЗ использовалась программа EQFRUX [6]. Переводные коэффициенты соответствовали трехкаскадному ФЭП производства компании Spectrolab (http://www.spectrolab.com/DataSheets/ TJCell/tj.pdf) и были взяты из программы General EQFlux (https://opensource.gsfc.nasa.gov/projects/ eqflux/index.php). При расчете падения мощности СБ учитывалось, что фоточувствительные элементы ФЭП находятся за слоем защитного стекла, которое частично ослабляет поток энергичных частиц РПЗ падающих на фотоэлектронные преобразователи СБ. Толщина защитного покрытия влияет на переводные коэффициенты для расчета эквивалентного флюенса. В расчетах было принято, что ФЭП находятся за слоем кварцевого стекла толщиной ${{\delta }_{п }} = 150$ мкм плотностью 2.2 г/см³ (как в [6]). Предполагалось, что задняя поверхность ФЭП закрыта элементами конструкции СБ с такими же ослабляющими свойствами как и у кварцевого стекла толщиной ${{\delta }_{{\text{з }}}} = 500$ мкм. Согласно [6] для учета повреждений ФЭП, нанесенных потоком частиц, падающих как с передней, так и с задней поверхности СБ, рассчитывались отдельно эквивалентные флюенсы для частиц, падающих на переднюю и заднюю поверхности СБ, а затем складывались.

Интегральный эквивалентный флюенс электронов 1 МэВ накопленный единицей площади ФЭП от начала перелета и до момента времени $t$ можно рассчитать следующим образом:

где ${{\bar {U}}_{n}}\left( {t,{\mathbf{r}}\left( t \right)} \right)$ – эквивалентный поток электронов 1 МэВ в зависимости от текущего времени $t$ и положения ${\mathbf{r}}\left( t \right)$ КА в околоземном пространстве. Функция $\,{{\bar {U}}_{n}}\left( {t,{\mathbf{r}}\left( t \right)} \right)$ рассчитывается численно с помощью программного комплекса IRENE как описано выше. Согласно [6, 7] относительную мощность СБ на конец САС можно записать в виде следующей полуэмпирической формулы: где $C = 0.1037$ и Φ_X$ = {\text{2}}{\text{.432}} \times {{10}^{{14}}}$ см^–2 – константы аппроксимации экспериментальных данных о деградации ФЭП из документации фирмы производителя ФЭП, ${{\Phi }_{{{\text{Г С О }}}}} = {\text{3}}{\text{.54}} \times {\text{1}}{{{\text{0}}}^{{14}}}$ см^–2 – эквивалентный флюенс электронов 1 МэВ накопленный ФЭП за время после выведения на ГСО и до конца САС, рассчитанный с помощью программного комплекса IRENE за период с 1.I.2021 года по 1.I.2036 года, $T$ – время выведения КА на ГСО с помощью ЭРДУ. Преимущество описания деградации мощности ФЭП именно формулой (9) состоит в том, что эта формула имеет под собой физические основы, и также данная формула содержит всего два неизвестных коэффициента определяемых из экспериментальных данных. Зависимость ${{\beta }_{{{\text{С А С }}}}}$ от суммарного эквивалентного флюенса 1 МэВ электронов ${{\Phi }_{n}}(T) + {{\Phi }_{{{\text{Г С О }}}}}$ показана на рис. 1. Далее введем следующие функции: Тогда

Для применения принципа максимума Понтрягина к системе с целевым функционалом (9) н-еобходимо к уравнениям движения КА добавить уравнение (10) с начальным условием ${{\beta }_{n}}(0) = 1 - С \ln \left( {1 + {{{{\Phi }_{{{\text{Г С О }}}}}} \mathord{\left/ {\vphantom {{{{\Phi }_{{{\text{Г С О }}}}}} {{{\Phi }_{X}}}}} \right. \kern-0em} {{{\Phi }_{X}}}}} \right).$

Но в связи с особенностями реализации расчета спектров частиц РПЗ в программном комплексе IRENE и резким нарастанием величин потоков частиц РПЗ на нижней границе внутреннего радиационного пояса Земли функция ${{U}_{n}}$ имеет разрывы частных производных и области с большими градиентами. Этот факт сильно осложняет численное интегрирование траекторий КА и повышает жесткость рассматриваемой системы дифференциальных уравнений. Аналогично работам [9, 11] для снижения жесткости полученной системы воспользуемся методом осреднения эквивалентного потока. Пусть осредненный эквивалентный флюенс $\Phi $ и осредненная относительная мощность СБ $\beta $ заданы следующим образом:

где функция осредненного эквивалентного потока $U\left( {r,i} \right)$ зависит только от радиуса $r(t) = \mu {{{{h}^{2}}(t)} \mathord{\left/ {\vphantom {{{{h}^{2}}(t)} {\xi (t)}}} \right. \kern-0em} {\xi (t)}}$ и наклонения $i(t) = \arccos \left( {{{\left( {2 - \tilde {\varphi }(t)} \right)} \mathord{\left/ {\vphantom {{\left( {2 - \tilde {\varphi }(t)} \right)} {\tilde {\varphi }(t)}}} \right. \kern-0em} {\tilde {\varphi }(t)}}} \right)$ оскулирующей орбиты КА. Расчет функции осредненного эквивалентного потока производится по следующей формуле где ${{T}_{r}}$ – период кеплеровского движения КА по круговой орбите с радиусом r, ${{n}_{a}}$ – число витков осреднения, ${{t}_{0}}$ – фиксированная дата, соответствующая 00.00.00 1.I.2020 года, ${{{\mathbf{r}}}_{d}}(t)$ – вектор-функция, описывающая положения КА в околоземном пространстве при движении в течении ${{n}_{a}}$ витков по круговой орбите радиусом r, наклонением $i$ и нулевой долготой восходящего узла.

Чтобы избежать разрывов частных производных функции U, гладкость которых необходима для численного интегрирования систем дифференциальных уравнений методами высокого порядка, аналогично работам [9, 11] построим гладкую аппроксимацию функции $U$ с помощью сплайнов 11-го порядка. Для этого был произведен расчет функции $U$ согласно формуле (11) на опорной сетке из $30 \times 40 = 1200$ круговых орбит в области от 0 до 90 градусов по наклонению и от 7071 до 63 371 км по радиусу. Число витков осреднения выбиралось таким образом, чтобы оно было целым и время движения по круговой орбите осреднения составляло не менее 3 суток. Кроме того, дополнительно ко всем операциям по построению гладкого сплайна $U$ из работ [9, 11] в данной работе применено дополнительное преобразование входных и выходных переменных следующего вида

где ${{R}_{{\text{E}}}} = {\text{6371}}{\text{.2}}$ км – средний радиус Земли, ${{U}_{s}},$ ${{r}_{b}},$ ${{r}_{s}},$ ${{i}_{s}}$ – константы преобразования, которые выбирались таким образом, чтобы вписать получаемую поверхность $W = W\left( {x,y} \right)$ в единичный куб ${{\left[ {0;1} \right]}^{3}}.$ Это преобразование позволило сделать более плавным изменение эквивалентного потока на границах внутреннего РПЗ и уменьшить погрешность аппроксимации значений осредненного эквивалентного потока, рассчитанных с помощью программного комплекса IRENE. Получившийся в итоге сплайн $U(r,i)\,$ представлен на рис. 2. Видно, что вклад внешнего электронного радиационного пояса в суммарный эквивалентный поток на несколько порядков меньше, чем вклад протонов внутреннего радиационного пояса. Таким образом, задача пересечения РПЗ с наименьшими потерями мощности СБ сводится к задаче оптимального пересечения области внутреннего радиационного пояса.

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

с начальным условием β(0) = 1 – C × $ \times \,\ln \left( {1 + {{{{\Phi }_{{{\text{Г С О }}}}}} \mathord{\left/ {\vphantom {{{{\Phi }_{{{\text{Г С О }}}}}} {{{\Phi }_{X}}}}} \right. \kern-0em} {{{\Phi }_{X}}}}} \right).$ В качестве целевого функционала в данной работе будем рассматривать значение неосредненной относительной мощности СБ на конец САС ${{\beta }_{n}}(T).$ Для упрощения задачи согласно результатам работ [9, 11] допустимо рассматривать в качестве целевого функционала значение осредненной относительной мощности СБ на конец САС $\beta (T).$

ОПТИМАЛЬНОЕ УПРАВЛЕНИЕ

Задачу на максимизацию относительной мощности СБ на конец САС можно формализовать как задачу оптимального быстродействия [19] с дополнительной фазовой переменной $\beta (t),$ удовлетворяющей уравнению (12), и дополнительным ограничением на правом конце следующего вида

Рассмотрим серию задач оптимального быстродействия с ограничением (13), где в каждой задаче ${{\beta }_{T}}$ постоянно. Постепенно увеличивая ${{\beta }_{T}}$ и решая получаемые краевые задачи можно найти максимальное значение $\beta (T),$ при котором еще существует решение краевой задачи. Таким образом, можно найти траекторию перелета с наибольшим возможным конечным значением относительной мощности СБ на конец САС.

Аналогично работе [19] воспользуемся принципом максимума Понтрягина для получения оптимального в смысле времени перелета управления системой (1)–(7), (12). Функция Понтрягина для данной системы будет иметь вид

где введены обозначения и ${{p}_{h}},$ ${{p}_{{ex}}},$ ${{p}_{{ey}}},$ ${{p}_{{ix}}},$ ${{p}_{{iy}}},$ ${{p}_{F}},$ ${{p}_{m}},$ ${{p}_{\beta }}$ – сопряженные переменные к $h,$ ${{e}_{x}},$ ${{e}_{y}},$ ${{i}_{x}},$ ${{i}_{y}},$ $F,$ $m$ и $\beta ,$ соответственно. Знак T в верхнем индексе после вектора-строчки здесь и далее будет означать операцию транспонирования. Запись $\left( {{\mathbf{g}},{\mathbf{A}}} \right)$ означает скалярное произведение векторов ${\mathbf{g}}$ и ${\mathbf{A}}$. Для данной задачи управление, удовлетворяющее условию максимума функции Понтрягина, имеет следующий вид:

В таком случае можно исключить из рассмотрения уравнения для массы КА $m$ и для сопряженной к ней переменной ${{p}_{m}},$ используя в уравнениях движения явную зависимость массы от времени $m = {{m}_{0}} - \left( {P{\text{/}}w} \right)t.$ Подставим выражения для оптимального управления (15), (16) в функцию Понтрягина (14) и получим выражение для гамильтониана

где $k = {h \mathord{\left/ {\vphantom {h {\left( {m\xi } \right)}}} \right. \kern-0em} {\left( {m\xi } \right)}},$ $A = \sqrt {A_{\tau }^{2} + A_{r}^{2} + A_{n}^{2}} .$ Тогда, если ввести векторные обозначения ${\mathbf{x}} = {{\left( {h,{{e}_{x}},{{e}_{y}},{{i}_{x}},{{i}_{y}}} \right)}^{T}}$ и ${\mathbf{p}} = {{\left( {{{p}_{h}},{{p}_{{ex}}},{{p}_{{ey}}},{{p}_{{ix}}},{{p}_{{iy}}}} \right)}^{T}}$ уравнения оптимального движения запишутся в следующей форме:

Как будет показано ниже к системе уравнений (18)–(23) применяется метод осреднения [19]. Для осредненной системы уравнений гамильтониан не будет зависеть от истинной долготы F, поэтому ${{p}_{F}} = {\text{const}}{\text{.}}$ А поскольку из условий трансверсальности на правом конце ${{p}_{F}}(T) = 0,$ то ${{p}_{F}} \equiv 0$ на всей траектории перелета КА. В таком случае уравнение для ${{p}_{F}}$ выпадает из рассмотрения. Если заметить, что уравнения (20) и (23) имеют общий множитель то несложно показать, что уравнение (23) имеет явное решение

Подставляя выражение (24) в остальные уравнения оптимального движения (18)–(23) и гамильтониан (17), можно исключить из рассмотрения уравнение (23) и получить систему уравнений оптимального движения КА следующего вида

где ${{С }_{0}} = С \exp \left( {{{ - \left( {1 - \beta (0)} \right)} \mathord{\left/ {\vphantom {{ - \left( {1 - \beta (0)} \right)} C}} \right. \kern-0em} C}} \right) = {\text{const}}{\text{.}}$ Также упростится выражение для гамильтониана, который примет вид

ОСРЕДНЕНИЕ ДВИЖЕНИЯ

Следует отметить, что система уравнений оптимального движения (25)–(28) содержит одну быструю переменную F, а все остальные переменные – медленные. Поэтому для повышения численной устойчивости и снижения вычислительной сложности краевой задачи принципа максимума к системе (25)–(28) можно применить метод осреднения [19]. Осреднение дифференциальных уравнений оптимального движения по времени за один виток проводится по следующей формуле

где $t$ – текущее время, ${{T}_{c}}$ – оскулирующий период орбиты КА, ${{{\mathbf{f}}}_{e}}\left( {{\mathbf{y}},F,t} \right)$ – правые части уравнений оптимального движения, ${\mathbf{y}} = {{\left( {{{{\mathbf{x}}}^{T}},\beta ,{{{\mathbf{p}}}^{T}}} \right)}^{T}},$ n = μ^–1 ×$ \times \,\,{{h}^{{ - 3}}}{{\left( {1 - e_{x}^{2} - e_{y}^{2}} \right)}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}$ – среднее движение КА и ${{ds} \mathord{\left/ {\vphantom {{ds} {dF}}} \right. \kern-0em} {dF}} \equiv {{dt} \mathord{\left/ {\vphantom {{dt} {dF}}} \right. \kern-0em} {dF}} = {{{{h}^{3}}} \mathord{\left/ {\vphantom {{{{h}^{3}}} {{{\xi }^{2}}.}}} \right. \kern-0em} {{{\xi }^{2}}.}}$ При работе с осредненными уравнениями оптимального движения уравнение для истинной долготы (26) исключалось из системы, поскольку правые части медленных уравнений после осреднения не зависят от истинной долготы F. И, напротив, при необходимости интегрировать неосредненные уравнения движения уравнение для истинной долготы (26) вновь добавлялось в систему.

КРАЕВАЯ ЗАДАЧА

Для удовлетворения условий принципа максимума необходимо решить краевую задачу для системы осредненных уравнений оптимального движения (25), (27), (28) со следующими краевыми условиями на левом и правом концах:

где граничные равноденственные элементы выражаются через высоту ГСО и параметры промежуточной орбиты следующим образом:

В четвертой строчке (31) под величиной ${{\bar {H}}_{{opt}}}(T)$ имеется в виду выражение для гамильтониана (29), осредненное по схеме (30) на момент конца перелета с помощью ЭРДУ. Условия из третьей и четвертой строчек (31) содержат значения равноденственных элементов, относительной мощности СБ и гамильтониана на правом конце. Для их расчета необходимо проинтегрировать систему уравнений оптимального движения (25), (27), (28) на промежутке $t \in \left[ {0;T} \right]$ с начальными значениями из первой и второй строчек (31) и применением схемы осреднения (30). Для этого необходимо знать время перелета $T$ и недостающие начальные значения сопряженных переменных ${\mathbf{p}}(0)$ и ${{p}_{\beta }}(0),$ которые запишем в виде вектора неизвестных параметров

Таким образом, для решения краевой задачи принципа максимума необходимо подобрать такие значения ${\mathbf{z}},$ что при интегрировании осредненной системы уравнений оптимального движения КА удовлетворились условия (31). Тогда краевая задача принципа максимума сводится к задаче решения нелинейной системы уравнений следующего вида: с неизвестными параметрами ${\mathbf{z}}$ из (32).

РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ

Для решения нелинейной системы уравнений (33), возникающей из принципа максимума, в данной работе применяется метод продолжения решения по параметру. Используемый вариант этого метода подробно описан в работах [19–21]. Согласно методу продолжения по параметру решение нелинейной системы уравнений (33) сводится к решению задачи Коши для системы обыкновенных дифференциальных уравнений вида

с начальным условием ${\mathbf{z}}\left( 0 \right) = {{{\mathbf{z}}}_{0}},$ где ${\mathbf{f}}\left( {{{{\mathbf{z}}}_{0}}} \right) = {\mathbf{b}},$ а $\tau $ – параметр продолжения. Решая численно каким-либо методом систему (34) на интервале $\tau \in \left[ {0;1} \right]$ получим решение исходной системы уравнений (33) как ${\mathbf{z}}* = {\mathbf{z}}(1).$

Для непрерывного перехода от задачи оптимального быстродействия, которая решается относительно легко, к задаче с заданной конечной мощностью СБ будем задавать падение мощности на конец САС в зависимости от параметра продолжения $\tau $ следующим образом:

где ${{\beta }_{{\min T}}}$ – падение мощности на конец САС на траектории оптимального быстродействия, а ${{\beta }_{{\max }}}$ – некоторое значение, такое что ${{\beta }_{{\min T}}} < {{\beta }_{{\max }}} < \beta (0) < 1.$ Тогда выражение для производной от ${\mathbf{f}}$ по параметру продолжения можно выписать в явном виде Следует отметить, что при ${{p}_{\beta }}(0) = 0$ уравнения оптимального движения КА (25), (27), (28) и гамильтониан (29) совпадают с уравнениями оптимального движения КА и гамильтонианом для осредненной задачи оптимального быстродействия. Поэтому в качестве начального приближения ${{{\mathbf{z}}}_{0}}$ будем использовать следующие значения где ${\mathbf{p}}{\text{*}}(0)$ и $T{\text{*}}$ получены из решения задачи оптимального быстродействия для перелета КА между промежуточной орбитой и ГСО.

Для вычисления невязок на правом конце траектории перелета КА, входящих в (33), при численном интегрировании осредненной системы оптимального движения использовался метод Дорманда–Принса 7(8) порядка [22]. Осреднение правых частей уравнений оптимального движения проводилось согласно (30) с помощью метода трапеций на равномерной сетке из 128 точек. Для вычисления частных производных от функции ${\mathbf{f}}$ по каждому из неизвестных параметров в (34) использовался метод конечных разностей первого порядка. При этом возникает численная ошибка [22], связанная с тем, что при последовательном интегрировании системы дифференциальных уравнений оптимального движения КА с близкими начальными условиями интегратор проходит при разных интегрированиях по разным значениям независимой переменной. Для исключения этой ошибки интегрировались одновременно шесть копий системы уравнений (25), (27), (28) с различными начальными условиями. Разрешение относительно производной ${{d{\mathbf{z}}} \mathord{\left/ {\vphantom {{d{\mathbf{z}}} {d\tau }}} \right. \kern-0em} {d\tau }}$ в (34) производилось путем численного решения системы линейных уравнений методом гауссового исключения.

Поскольку наибольшая эффективность использования метода продолжения решения по параметру достигается при использовании методов численного интегрирования высокого порядка [19], в данной работе для интегрирования (34) также использовался метод Дорманда–Принса 7(8) порядка. Однако для еще большего увеличения длины шага по параметру продолжения и для снижения количества отброшенных шагов в данной работе после каждого принятого шага запускалась процедура коррекции полученного значения вектора неизвестных параметров ${{{\mathbf{z}}}_{1}}.$ Процедура коррекции состояла в решении при фиксированном $\tau {\text{*}}$ следующего уравнения

относительно вектора ${\mathbf{z}}\left( {\tau {\text{*}}} \right).$ В отличии от работ [9, 11], использующих метод Ньютона для коррекции, решение уравнения (37) проводилось методом Левенберга–Марквардта с начальным приближением ${\mathbf{z}}\left( {\tau {\text{*}}} \right) = {{{\mathbf{z}}}_{1}}.$ Применение метода Левенберга–Марквардта при больших невязках возникающих при интегрировании (34) позволяет за меньшее число итераций сойтись к решению (37) и повысить точность решения краевой задачи принципа максимума.

Для проверки насколько отличается полученное управление от оптимального вследствие применения метода осреднения после каждого принятого шага и проведенной процедуры коррекции проводилось интегрирование неосредненных уравнений движения (25)–(28). При этом начальные значения сопряженных переменных и время перелета брались из ${\mathbf{z}}\left( {\tau {\text{*}}} \right),$ сопряженная переменная к истинной долготе полагалась равной нулю в течение всего перелета ${{p}_{F}}(t) \equiv 0,$ а начальное значение истинной долготы полагалось нулевым $F(0) \equiv 0.$ Кроме того, после получения неосредненных траекторий на этих траекториях рассчитывались с помощью программного комплекса IRENE неосредненный эквивалентный флюенс на конец перелета на ГСО (8) и неосредненная относительная мощность СБ на конец САС (9). Численное интегрирование неосредненных траекторий оптимального движения также проводилось методом Дорманда–Принса 7(8) порядка.

ОПТИМИЗАЦИЯ АРГУМЕНТА ПЕРИГЕЯ ПРОМЕЖУТОЧНОЙ ОРБИТЫ

Рассмотрим несколько промежуточных орбит, параметры которых указаны в табл. 1. Эта таблица определяет энергетические параметры промежуточных орбит, но не фиксирует для каждой орбиты долготу восходящего узла ${{\Omega }_{0}}$ и аргумент перигея ${{\omega }_{{\pi ,0}}}.$ В задаче оптимального быстродействия (в центральном поле притяжения Земли) для перелета между некомпланарными эллиптической и круговой орбитами оптимальным значением аргумента перигея промежуточной орбиты является ${{\omega }_{{\pi ,0}}} = 0,$ а значение ${{\Omega }_{0}}$ на время выведения не влияет. Для осредненного эквивалентного потока в силу симметрии относительно оси вращения Земли значение ${{\Omega }_{0}}$ не будет влиять на относительную мощность на конец САС. Поэтому далее будем рассматривать задачу оптимального быстродействия с учетом второй зональной гармоники разложения гравитационного потенциала Земли для каждой промежуточной орбиты из таблицы при ${{\Omega }_{0}} = 0$ и различных значениях ${{\omega }_{{\pi ,0}}}.$

Результаты расчета осредненной и неосредненной относительной мощности СБ на конец САС в зависимости от ${{\omega }_{{\pi ,0}}}$ приведены на рис. 3. Расчеты проведены для перелета на ГСО с промежуточной орбиты №1 из табл. 1. Виден периодический характер зависимости осредненной относительной мощности от аргумента перигея промежуточной орбиты. Период этой зависимости составляет 180 градусов. На рис. 4 изображена зависимость осредненной относительной мощности СБ на конец САС от аргумента перигея для траекторий перелета на ГСО со всех промежуточных орбит из табл. 1. На рис. 5 отложены значения осредненной относительной мощности СБ на конец САС и потребной характеристической скорости для выведения КА по полученным траекториям на ГСО при различных значениях ${{\omega }_{{\pi ,0}}}.$ Точки, соответствующие одной строчке из табл. 1, соединены линиями. Из рис. 5 видно, что в рамках осредненной модели эквивалентного потока существует при одних и тех же затратах характеристической скорости два значения ${{\omega }_{{\pi ,0}}},$ одному из которых соответствует траектория оптимального быстродействия с большим значением $\beta \left( T \right).$

ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ РЕШЕНИЯ КРАЕВОЙ ЗАДАЧИ

Для каждой орбиты из табл. 1 выбиралось значение аргумента перигея промежуточной орбиты с минимальной и максимальной относительной мощностью СБ на конец САС на траекториях оптимального быстродействия и решалась краевая задача, сведенная выше к задаче Коши (34). Величина ${{\beta }_{{\max }}}$ во всех полученных краевых задачах полагалась равной ${{\beta }_{{\max }}} = 0.906,$ что меньше $\beta (0) = {\text{0}}{\text{.90679615}}{\text{.}}$ При некотором значении параметра продолжения $\tau = {{\tau }_{f}} < 1$ наблюдалось существенное увеличение невязок на правом конце. При этом корректор после каждого шага не мог свести невязки на правом конце к нулю. В этот момент интегрирование по параметру продолжения прерывалось. Следует однако отметить, что подобное прерывание интегрирования системы (34) не гарантирует, что нарастание невязок на правом конце произошло из-за отсутствия решения краевой задачи при $\tau > {{\tau }_{f}}.$ Но тем не менее, полученные траектории позволяют оценить величину максимального возможной мощности СБ на конец САС.

Полученные в ходе решения (34) значения ${\mathbf{z}}(\tau ),$ использовались для интегрирования неосредненных уравнений оптимального движения (25)–(28) и расчета неосредненной относительной мощности СБ на конец САС с помощью программного комплекса IRENE. На рис. 6 представлены результаты расчета осредненной и неосредненной относительной мощности СБ на конец САС на полученных неосредненных траекториях. Все расчеты, представленные на рис. 6, соответствуют выбору ${{\omega }_{{\pi ,0}}}$ с минимальной относительной мощностью СБ на конец САС на траекториях оптимального быстродействия. Осредненная относительная мощность нарисована на всех графиках рис. 6 светлыми линиями, неосредненная – темными. В заголовке каждого графика на рис. 6 указано наклонение и выбранное значение аргумента перигея промежуточной орбиты. На абсциссах графиков отложена величина дополнительных затрат характеристической скорости на перелет с помощью ЭРДУ по сравнению с траекторией оптимального быстродействия для выбранных параметров промежуточной орбиты. Все графики показывают монотонное увеличение относительной мощности с ростом затрат характеристической скорости. Максимальное увеличение мощности СБ на конец САС в этой группе расчетов составляет 0.6% от мощности СБ на начальной орбите. При этом дополнительные затраты характеристической скорости не превышают 400 м/с.

На рис. 7 представлены результаты расчета осредненной и неосредненной относительной мощности СБ на конец САС на полученных неосредненных траекториях, соответствующих ${{\omega }_{{\pi ,0}}}$ с максимальной относительной мощностью СБ на конец САС на траекториях оптимального быстродействия. Все обозначения аналогичны рис. 6. Стоит только отметить, по абсциссе отложены дополнительные затраты характеристической скорости по сравнению с траекторией оптимального быстродействия с соответствующей промежуточной орбиты, но с аргументом перигея соответствующим минимальной относительной мощностью СБ. Такой выбор отсчета дополнительных затрат характеристической скорости позволяет проиллюстрировать результаты оптимизации относительной мощности одновременно по ${{\omega }_{{\pi ,0}}}$ и по ${\mathbf{z}}.$ Следует заметить, что зависимости, представленные на рис. 7 имеют немонотонный характер: относительная мощность растет, а дополнительные затраты характеристической скорости сначала убывают, а потом нарастают. Такой характер изменения характеристической скорости может быть связан с тем, что начальная траектория, с которой запускается метод продолжения по параметру, не является оптимальной по быстродействию в смысле выбора аргумента перигея промежуточной орбиты. Поэтому сначала при оптимизации увеличивается относительная мощность за счет снижения времени перелета и, как следствие, эквивалентного флюенса. Затем в какой-то момент эти возможности исчерпываются и дальнейшее увеличение мощности СБ на конец САС получается только за счет увеличения времени перелета. Максимальное увеличение мощности СБ на конец САС в этой группе расчетов составляет 0.57% от мощности СБ на начальной орбите. При этом дополнительные затраты характеристической скорости не превышают 1258 м/с.

По результатам оптимизации аргумента перигея промежуточной орбиты и применения принципа максимума удалось в зависимости от промежуточной орбиты увеличить мощность СБ на конец САС на 0.16–0.66% от мощности СБ на промежуточной орбите. При этом дополнительные затраты характеристической скорости составили от 13 до 1087 м/с.

На рис. 8 показаны зависимости необходимой толщины защитного стекла от длительности перелета при различных фиксированных на конец САС относительных мощностях СБ. Значения фиксированной на конец САС относительной мощности указаны на изолиниях. График получен путем варьирования толщины лицевого защитного стекла СБ на траекториях полученных для случая перелета с промежуточной орбиты с наклонением ${{i}_{0}}$ = 46° и аргументом перигея ${{\omega }_{{\pi ,0}}}$ = 180°.

На рис. 9 показано как менялась траектория перелета с промежуточной орбиты с наклонением ${{i}_{0}}$ = 46° и аргументом перигея ${{\omega }_{{\pi ,0}}}$ = 180° по мере увеличения параметра продолжения и, соответственно, увеличения относительной мощности СБ на конец САС. Траектории при каждом значении параметра продолжения представлены в осях “радиус–наклонение”. Две пунктирные линии на каждом графике показывают зависимости высоты перигея и апогея от наклонения на траектории оптимального быстродействия. Две непрерывные линии соответствуют зависимости высоты перигея и апогея от наклонения на траектории, полученной при данном значении параметра продолжения. Оттенками серого на каждом графике показана величина эквивалентного потока электронов 1 МэВ, падающего на СБ. На графиках видна тенденция к вынесению перигейных участков траектории выше границы внутреннего РПЗ, а также к увеличению максимальной высоты апогея. Также на рис. 9 можно заметить, что по мере увеличения конечной относительной мощности СБ изменение наклонения орбиты все меньше происходит в области пересечения траекторией КА внутреннего РПЗ и все больше – вне этой области. Этот факт частично подтверждает правильность расчетов, так как видно, что с уменьшением наклонения эквивалентный поток внутри радиационных поясов растет, и поэтому выгоднее в плане радиационной нагрузки пересекать РПЗ (другими словами набирать высоту) в области больших наклонений.

ЗАКЛЮЧЕНИЕ

В работе рассмотрены возможности снижения радиационной деградации мощности СБ КА на конец САС за счет изменения формы траектории выведения КА с различных промежуточных орбит на ГСО с использованием ЭРДУ и СБ на основе трехкаскадных ФЭП. Математический аппарат построения оптимальных по быстродействию траекторий с фиксированной относительной мощностью СБ на конец САС позволил увеличить мощность СБ на конец САС на 0.16–0.66% от мощности СБ на начало перелета в зависимости от параметров промежуточных орбит. Дополнительные затраты характеристической скорости увеличились относительно траекторий оптимального быстродействия на 13–1087 м/с. При оптимизации учтена вторая зональная гармоника в разложении гравитационного потенциала притяжения Земли. Относительно небольшое увеличение конечной мощности СБ обуславливается большей радиационной стойкостью трехкаскадных ФЭП по сравнению с кремниевыми и тем фактом, что большую часть времени перелета КА находится в области с малым эквивалентным потоком, поскольку существенное влияние на деградацию трехкаскадных ФЭП имеет только внутренний РПЗ. В работе рассмотрена постановка задачи, где к уравнениям движения КА добавляется уравнение для учета радиационной деградации мощности СБ на конец САС. Применен принцип максимума к полученной задаче, получено в явном виде решение сопряженного уравнения к относительной мощности СБ на конец САС.

Работа выполнена за счет средств гранта Российского научного фонда (соглашение № 16-19-10429).

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