Космические исследования, 2019, T. 57, № 4, стр. 308-320

ВВЕДЕНИЕ

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

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

В данной работе представлен способ снижения дозы радиации, полученной бортовыми системами при перелете ММБ с низкой круговой орбиты на ГСО с помощью ядерной ЭРДУ. Данный способ связан со специальным выбором управления движением центра масс ММБ, которое позволяет снизить суммарное количество заряженных частиц РПЗ, поглощаемых бортовыми системами за время перелета. Для перелета с низкой круговой орбиты на ГСО существенная часть накопленной дозы является дозой радиации, полученной от поглощения заряженных частиц РПЗ. Известно [3], что распределение заряженных частиц РПЗ в околоземном пространстве неоднородно и, следовательно, поглощенная доза будет зависеть не только от продолжительности перелета, но и от того, где внутри РПЗ проходит траектория перелета. Меняя расположение траектории относительно поясов РПЗ можно снизить дозовую нагрузку на бортовые системы ММБ.

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

В работе рассматривались перелеты ММБ с модулем полезной нагрузки (МПН) с низкой круговой орбиты с ненулевым наклонением на ГСО с помощью ЭРДУ и ядерной энергетической установки (ЯЭУ). Рассматривался вариант ММБ, описанный в работе [1], c электрической мощностью ЯЭУ равной 1070 кВт. Все используемые далее значения проектных параметров ММБ соответствуют значениям из работы [1], которые согласно этой работе являются оптимальными для данного варианта ММБ. Предполагалось, что ММБ и МПН выводятся на начальную орбиту отдельно, по двухпусковой схеме. Первым пуском выводится ММБ, который с помощью ХРБ переводится на низкую круговую орбиту с ненулевым наклонением, которую будем далее называть радиационно безопасной орбитой (РБО). Вторым пуском выводится МПН. В составе МПН также есть блок расходных компонентов, который содержит рабочее тело для ЭРДУ ММБ и последующего перелета на ГСО. МПН также с помощью ХРБ переводится на РБО и совершает стыковку к ММБ. Далее связка ММБ + МПН с помощью ЭРДУ и бортового ядерного источника энергии совершает перелет на ГСО.

Будем считать, что ММБ и МПН в совместной фазе полета движутся только под действием силы притяжения Земли и силы тяги ЭРДУ. Потенциал притяжения Земли будем считать центральным ньютоновским. Всеми остальными возмущающими силами пренебрежем. Будем считать, что ЭРДУ всегда работает в режиме максимальной тяги, и тяга в течение всего перелета постоянна и равна $P = 27$ Н. Электрическая мощность, необходимая для работы ЭРДУ, генерируется бортовым ядерным источником энергии. Удельный импульс ЭРДУ также будем считать постоянным в течение всего перелета и равным ${{I}_{{sp}}} = 4762$ с. Кроме того, будем считать, что на возможную ориентацию вектора тяги ЭРДУ не наложено ограничений, и при этом соответствующие развороты связки ММБ + МПН будут происходить без затрат рабочего тела.

Пусть в начальный момент связка ММБ и МПН имеет массу ${{m}_{0}} = 40$ т и находится на РБО, высота которой равна ${{h}_{{{\text{Р Б О }}}}} = 800$ км, а наклонение – ${{i}_{0}} = {{51.6}^{^\circ }}.$ Конечной орбитой будем считать ГСО, параметры которой примем следующими: высота равна ${{h}_{T}} = 35793$ км, эксцентриситет и наклонение равны нулю. Угловое положение на начальной орбите фиксировано, а на конечной – произвольно. Для оценки радиационной нагрузки на бортовые системы ММБ будем рассчитывать дозу космического ионизирующего излучения, поглощенную некоторым чувствительным элементом (ЧЭ), который расположен в приборно-агрегатном отсеке ММБ. Для приближенного расчета поглощенной ЧЭ дозы будем считать, что ЧЭ находится за бесконечной алюминиевой стенкой толщиной $\lambda = 4$ мм. По причинам, описанным в разделе целевой функционал, будем считать, что ЧЭ поглощает только протоны и электроны РПЗ. Будем искать такие программы управления вектором тяги ЭРДУ, что соответствующие траектории многовиткового перелета с начальной орбиты на конечную имели бы наименьшую дозу радиации, поглощенную ЧЭ за время перелета.

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

Введем орбитальную систему координат $OXYZ$ следующим образом. Начало координат $O$ расположено в центре масс связки ММБ и МПН. Радиальная ось $OY$ направлена из центра Земли в начало координат O, бинормальная ось $OZ$ сонаправлена с вектором орбитального кинетического момента ММБ. Трансверсальная ось $OX$ лежит в плоскости оскулирующей орбиты, перпендикулярна оси $OY$ и направлена в сторону движения ММБ по орбите. Будем задавать направление вектора реактивного ускорения, создаваемого ЭРДУ, с помощью угла тангажа и угла рысканья. Углом тангажа $\vartheta $ будем называть угол между проекцией вектора тяги на плоскость $OXY$ оскулирующей орбиты ММБ и трансверсальным направлением $OX.$ Углом рысканья $\psi $ будем называть угол между вектором тяги и плоскостью $OXY$ оскулирующей орбиты ММБ. Тогда проекции реактивного ускорения ММБ на оси орбитальной системы координат будут иметь вид:

где $\delta $ – функция включения ЭРДУ (при $\delta = 0$ двигательная установка выключена, при $\delta = 1$ – включена), $P$ – максимальная величина реактивной тяги ЭРДУ, $m$ – текущая масса связки ММБ и МПН.

Чтобы избежать особенностей в окрестности нулевых эксцентриситета и наклонения будем использовать для описания движения ММБ на орбите уравнения в равноденственных элементах:

где равноденственные элементы $h,{{e}_{x}},{{e}_{y}},{{i}_{x}},{{i}_{y}},F$ выражаются через фокальный параметр $p$ оскулирующей орбиты, эксцентриситет e, аргумент перицентра ω, наклонение i, долготу восходящего узла Ω, истинную аномалию $\nu $ следующим образом: $h = \sqrt {p{\text{/}}\mu } ,$ ${{e}_{x}} = e\cos \left( {\Omega + \omega } \right),$ ${{e}_{y}} = e\cos \left( {\Omega + \omega } \right),$ ${{i}_{x}} = \operatorname{tg} \left( {i{\text{/}}2} \right)\cos \Omega ,$ ${{i}_{y}} = \operatorname{tg} \left( {i{\text{/}}2} \right)\sin \Omega ,$ $F = \nu + \omega + \Omega .$ В уравнениях также используются обозначения: $\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},$ $w = {{I}_{{sp}}}g$ – скорость истечения рабочего тела ЭРДУ, $g = 9.80665$ м/с² – ускорение свободного падения, $\mu $ – гравитационный параметр Земли.

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

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

Для расчета потоков поглощенных ЧЭ электронов и протонов использовались модели AE8/AP8 MIN, AE8/AP8 MAX [5–7] и AE9/AP9 [8]. Расчеты проводились с помощью программного комплекса IRENE (https://www.vdl.afrl.af.mil/programs/ae9ap9/ downloads.php). Согласно рекомендациям работы [9] для расчета магнитных координат при использовании моделей AE8 MIN, AE8 MAX, AP8 MIN в качестве моделей магнитного поля внутренних источников Земли использовалась модель Jensen и Cain [10]. Для модели AP8 MAX использовалась экстраполированная на эпоху 1970 модель магнитного поля внутренних источников GSFC 12/66 [11]. Магнитное поле внешних токов магнитосферы Земли при использовании всех моделей типа AE8/AP8 не учитывалось. При расчете потоков частиц РПЗ при помощи модели AE9/AP9 для вычисления магнитных координат использовалась связка модели внутренних источников IGRF 2015 с моделью внешних токов Olson-Pfitzer Quiet (OPQ77). На основе рассчитанных потоков протонов и электронов РПЗ расчет поглощенной ЧЭ дозы проводился с использованием программы Shieldose2 [12]. Данная программа рассчитывает дозу, поглощенную чувствительным элементом за алюминиевой защитой некоторой толщины и геометрии от всенаправленного потока электронов и протонов произвольного спектра. В данной работе все дозы рассчитывались, предполагая, что чувствительный элемент ММБ состоит из кремния и находится за бесконечной алюминиевой стенкой толщиной $\lambda = 4$ мм (или примерно 1 г/см²).

Суммарную дозу, поглощенную ЧЭ за весь перелет, можно записать следующим образом:

где $T$ – время перелета ММБ с начальной орбиты на ГСО, ${{N}_{{D1}}}\left( {t,{\mathbf{r}}(t)} \right)$ – мощность дозы, поглощаемая ЧЭ в зависимости от текущего времени $t$ и положения ${\mathbf{r}}(t)$ ММБ в околоземном пространстве. Для применения принципа максимума Понтрягина к системе с целевым функционалом (8) необходимо к уравнениям движения добавить следующее уравнение с начальным условием ${{D}_{1}}(0) = 0$. При этом функция ${{N}_{{D1}}}$ рассчитывается численно на каждом шаге интегрирования с помощью программного комплекса IRENE. В таком случае численное интегрирование траекторий сильно осложняется наличием у функции ${{N}_{{D1}}}$ разрывов частных производных и областей с большими градиентами, что сильно увеличивает количество шагов численного интегрирования и повышает жесткость рассматриваемой системы дифференциальных уравнений. Как будет показано в разделе численные результаты, далее можно рассматривать упрощенную зависимость мощности дозы от времени и положения ММБ в околоземном пространстве. Будем считать, что мощность дозы зависит только от расстояния до центра Земли r, наклонения $i$ и истинной долготы F, и не зависит от текущего времени. Для снижения жесткости системы уравнений движения МКБ рассмотрим осредненную по времени за ${{n}_{a}}$ витков на круговой орбите радиусом $r$ и наклонением $i$ функцию мощности дозы где ${{T}_{r}}$ – период кеплеровского движения ММБ по круговой орбите с радиусом r. При осреднении (9) долгота восходящего узла и аргумент перигея рассматриваемых орбит считались нулевыми. Чтобы избежать разрывов производных функции мощности дозы, построим гладкую аппроксимацию функции ${{\tilde {N}}_{D}}$ на некоторой сетке на плоскости $\left( {i,r} \right)$ и обозначим эту функцию ${{N}_{D}}.$ Для возможности применения метода численного интегрирования высокого порядка была построена функция ${{N}_{D}},$ обладающая непрерывными частными производными вплоть до 10 порядка включительно. Для этого был произведен расчет функции ${{\tilde {N}}_{D}}(r,i)$ согласно (9) на некоторой опорной сетке из $20 \times 25 = 500$ круговых орбит при числе витков осреднения ${{n}_{a}} = 10.$ Затем по полученным значениям для сглаживания шумовой составляющей был построен кубический сглаживающий сплайн с помощью функции spaps в среде MATLAB (https:// www.mathworks.com/help/curvefit/spaps.html). Далее на некоторой другой сетке были вычислены значения сглаживающего сплайна, и был построен обыкновенный сплайн 11-го порядка при помощи функции spapi в среде MATLAB (https://www. mathworks.com/help/curvefit/spapi.html). Именно этот сплайн был использован при численном интегрировании уравнений движения ММБ. После всех преобразований будем считать, что целевым функционалом задачи является суммарная поглощенная ЧЭ доза, которая рассчитывается как

Тогда систему уравнений движения ММБ нужно дополнить уравнением

с начальным условием $D(0) = 0,$ где $r = \mu {{{{h}^{2}}} \mathord{\left/ {\vphantom {{{{h}^{2}}} \xi }} \right. \kern-0em} \xi }$ и $i = \arccos \left( {{{\left( {2 - \tilde {\varphi }} \right)} \mathord{\left/ {\vphantom {{\left( {2 - \tilde {\varphi }} \right)} {\tilde {\varphi }}}} \right. \kern-0em} {\tilde {\varphi }}}} \right).$ Для оценки влияния различных моделей потоков частиц РПЗ на возможность снижения дозы радиации было построено указанным выше способом четыре функции ${{N}_{D}}.$ Три из них являются средними по времени и состоянию магнитосферы мощностями дозы для моделей AE8/AP8 MIN, AE8/AP8 MAX и AE9/AP9. Последняя же функция ${{N}_{D}}$ содержит помимо информации о средних потоках частиц РПЗ еще и информацию о степени неопределенности потоков в областях околоземного пространства, где отсутствуют измерения. При расчете этой функции величины потоков частиц РПЗ в модели AE9/AP9 рассматривались как случайные поля с малой дисперсией в областях, где присутствуют измерения с различных спутников, и с большой дисперсией там, где измерений нет. Четвертая функция ${{N}_{D}}$ является аппроксимацией 95-го перцентиля мощности дозы, рассчитанного по формуле (9) по 40 реализациям случайного поля потоков частиц РПЗ модели AE9/AP9 для каждой из 500 орбит опорной сетки. Все полученные функции мощности дозы представлены на рис. 1. Значения мощности дозы для каждой полученной функции ${{N}_{D}}$ показаны оттенками серого. Подробно четвертая функция ${{N}_{D}}$ показана в виде трехмерной поверхности на рис. 2, где отмечены области, соответствующие внутреннему и внешнему радиационным поясам Земли.

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

Будем решать задачу оптимального быстродействия [13] с дополнительным уравнением (11) для накопленной дозы $D$ и дополнительным условием на правом конце:

Рассмотрим серию задач оптимального быстродействия с условием (12), где в каждой задаче ${{D}_{f}}$ постоянно. Последовательно снижая значение ${{D}_{f}}$ от задачи к задаче, начиная от величины поглощенной дозы на траектории оптимального быстродействия ${{D}_{{\min T}}}$, и решая получающиеся краевые задачи, можно найти минимальное значение ${{D}_{f}} = {{D}_{{\min }}}$, при котором еще существует решение краевой задачи. Таким образом, можно найти траекторию перелета с наименьшим конечным значением поглощенной дозы.

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

где введены обозначения и p_h, p_ex, p_ey, p_ix, p_iy, p_F, p_m, p_D – сопряженные переменные к h, e_x, e_y, i_x, i_y, F, $m$ и D, соответственно. Для рассматриваемой задачи оптимального быстродействия управление, удовлетворяющее условию максимума функции Понтрягина, имеет следующий вид:

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

где $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}} ,$ ${{H}_{F}} = {{p}_{F}}{{{{\xi }^{2}}} \mathord{\left/ {\vphantom {{{{\xi }^{2}}} {{{h}^{3}}}}} \right. \kern-0em} {{{h}^{3}}}},$ ${{H}_{D}} = {{p}_{D}}{{N}_{D}}.$ Тогда, если ввести векторные обозначения ${\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}},$ уравнения оптимального движения запишутся в следующей форме:

Из (21) видно, что переменная, сопряженная к поглощенной дозе, постоянна в течение всего перелета ${{p}_{D}} = {\text{const}}$ и имеет смысл множителя Лагранжа в поставленной задаче оптимизации.

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

Также упростится выражение для гамильтониана, который примет вид

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

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

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

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

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

где граничные равноденственные элементы выражаются через высоту ГСО, высоту и наклонение РБО следующим образом: где ${{R}_{{\text{E}}}} = 6371$ км – средний радиус Земли. В последней строчке (28) под величиной ${{\bar {H}}_{{{\text{opt}}}}}(T)$ имеется в виду выражение для гамильтониана (26), осредненное по схеме (27) на момент конца перелета. Условия из второй, третьей и четвертой строчек (28) содержат значения равноденственных элементов, дозы и гамильтониана на правом конце. Для их расчета необходимо проинтегрировать систему уравнений оптимального движения (22), (24), (25) на промежутке $t \in \left[ {0;T} \right]$ с начальными значениями из первой и второй строчек (28) и применением схемы осреднения (27). Для этого необходимо знать время перелета $T$ и недостающие начальные значения сопряженных переменных ${\mathbf{p}}(0)$ и ${{p}_{D}},$ которые запишем в виде вектора неизвестных параметров

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

с неизвестными параметрами ${\mathbf{z}}$ из (29).

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

Для решения нелинейной системы уравнений (30), возникающей из принципа максимума, в данной работе применяется метод продолжения решения по параметру. Описание данного метода и примеры его использования приведены в работах [13–16]. Далее в этом разделе вкратце описана суть данного метода.

Пусть задано некоторое начальное приближение ${{{\mathbf{z}}}_{0}}$ решения системы уравнений (30), и пусть ${\mathbf{f}}\left( {{{{\mathbf{z}}}_{0}}} \right) = {\mathbf{b}}.$ Тогда рассмотрим параметрическое семейство систем нелинейных уравнений, зависящее от скалярного параметра $\tau $ и имеющее вид:

Начальное приближение ${{{\mathbf{z}}}_{0}}$ является решением системы уравнений (31) при $\tau = 0$ по построению: ${\mathbf{F}}\left( {{{{\mathbf{z}}}_{0}},0} \right) = 0.$ Заметим, что исходная система уравнений (30) содержится в семействе (31) при $\tau = 1{\text{:}}$ ${\mathbf{F}}\left( {{\mathbf{z}},1} \right) = {\mathbf{f}}\left( {\mathbf{z}} \right) = 0.$ Для нахождения семейства решений (31) ${\mathbf{z}} = {\mathbf{z}}\left( \tau \right),$ зависящего от параметра τ, продифференцируем по $\tau $ систему уравнений (31) и, разрешив полученную систему относительно производной по параметру, получим систему обыкновенных дифференциальных уравнений вида

с начальными условиями

Решая полученную систему дифференциальных уравнений каким-либо методом при $\tau \in \left[ {0;1} \right],$ получаем решение исходной системы уравнений (30) как ${\mathbf{z}}* = {\mathbf{z}}(1).$ Таким образом, решение исходной системы нелинейных уравнений сводится к решению задачи Коши для системы обыкновенных дифференциальных уравнений (32), (33). В качестве начального приближения ${{{\mathbf{z}}}_{0}}$ для решения краевой задачи (30) в данной работе использовались значения ${\mathbf{p}}(0),$ $T$ из решения задачи оптимального быстродействия для перелета ММБ между РБО и ГСО. При этом начальное значение${{p}_{D}}$ полагалось равным нулю.

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

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

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

ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Практический интерес имеет только снижение конечной дозы ниже ее значения на траектории оптимального быстродействия, то есть когда ${{D}_{f}} < {{D}_{{\min T}}}.$ Но величина максимального снижения поглощенной дозы заранее неизвестна, поэтому проводилась следующая процедура. Сначала требуемое значение поглощенной дозы полагалось ${{D}_{f}} = 0$ и запускалось интегрирование по параметру продолжения системы (32), (33). При некотором значении параметра продолжения $\tau = {{\tau }_{f}} < 1$ наблюдалось существенное увеличение невязок на правом конце. При этом часто либо минимальная высота перицентра на траектории опускалась ниже 300 км, либо максимальная высота апоцентра поднимались выше орбиты Луны. Матрица в правой части (32) при этом оказывалась необратимой. В этот момент интегрирование по параметру продолжения прерывалось. Далее если шаги по $\tau $ были слишком большие и была необходимость получить более подробное решение ${\mathbf{z}}\left( \tau \right),$ $\tau \in \left[ {0;{{\tau }_{f}}} \right],$ то интегрирование системы (32), (33) перезапускалось со значением ${{D}_{f}},$ равным значению осредненной дозы $D(T)$ при $\tau = {{\tau }_{f}}.$ При этом интегрирование проводилось в пределах $\tau \in \left[ {0;{{\tau }_{f}}} \right],$ и выбирался более мелкий шаг по параметру продолжения. Следует отметить, что траектории при $\tau = {{\tau }_{f}}$ полученные таким образом из-за ошибок численных методов необязательно имеют наименьшую в принципе возможную величину конечной поглощенной дозы. Но эти траектории позволяют получить оценку максимального возможного снижения поглощенной дозы радиации на конец перелета.

Указанная выше процедура проводилась для всех четырех функций мощности дозы ${{N}_{D}}.$ Далее полученные неизвестные параметры осредненной краевой задачи ${\mathbf{z}}\left( \tau \right)$ использовались для интегрирования неосредненных уравнений оптимального движения (22)–(25) и расчета неосредненной дозы ${{D}_{1}}.$ Зависимость осредненной (10) и неосредненной дозы (8) от времени перелета на полученных траекториях в случае модели AE8/AP8 MAX показана на рис. 3 слева. Каждой точке на графике соответствует траектория при определенном значении параметра продолжения. Видно, что осредненная и неосредненная дозы достаточно хорошо совпадают. Ошибка аппроксимации конечной неосредненной дозы осредненной дозой на полученных траекториях в случае моделей AE8/AP8 MIN, AE8/AP8 MAX и среднего значения AE9/AP9 не превышает 2.2%. Поскольку последняя функция ${{N}_{D}}$ является аппроксимацией 95-го перцентиля мощности дозы, рассчитанного по 40 реализациям случайного поля, то значения этой функции и интеграл от нее по траектории будет зависеть от конкретных реализаций случайного поля потоков частиц РПЗ. Если в набор реализаций попали случаи с низкой мощностью дозы вдоль траектории, то неосредненная доза может существенно отличаться от осредненной дозы. Но все же в большинстве случаев неосредненная доза не будет превосходить осредненную. Поэтому осредненную дозу для четвертой функции ${{N}_{D}}$ можно рассматривать как оценку накопленной дозы в худшем случае. На рис. 3 справа приведена зависимость осредненной и неосредненной дозы от времени перелета на полученных траекториях для четвертой функции ${{N}_{D}}.$ Неосредненная доза рассчитывалась четыре раза с различными начальными значениями генератора случайных чисел модели потоков РПЗ. Видно, что все четыре расчета неосредненной дозы отличаются примерно на константу от осредненной дозы. Но при этом следует отметить, что осредненная доза для всех моделей мощности дозы является хорошим индикатором относительной дозовой нагрузки: если осредненная доза от траектории к траектории снижается, то и, скорее всего, снизится неосредненная доза. Поэтому принятое упрощение о зависимости мощности дозы только от радиуса и наклонения оскулирующей орбиты вполне допустимо в данной работе.

Для сравнения возможностей снижения конечной дозы для каждой функции мощности дозы был проведен расчет отношения неосредненной дозы на полученных траекториях к неосредненной дозе на траектории оптимального быстродействия. Результаты расчетов этого отношения в зависимости от времени перелета показаны на рис. 4 и в зависимости от разности затрат характеристической скорости на полученных траекториях и на траектории оптимального быстродействия – на рис. 5. На этих рисунках видно, что предложенная выше методика снижения дозовой нагрузки на траектории выведения ММБ с помощью специального управления работоспособна и позволяет получить траектории с дозой, меньше чем на траектории оптимального быстродействия. Максимально удается снизить конечную дозу радиации на 25–38% от дозы на траектории оптимального быстродействия в зависимости от модели потоков РПЗ. При этом время выведения увеличивается на 4–7% от оптимального времени выведения, а дополнительные затраты характеристической скорости увеличиваются относительно траектории оптимального быстродействия на 320–560 м/с.

На рис. 6 показано как менялась траектория перелета в случае модели AE8/AP8 MAX по мере увеличения параметра продолжения $\tau $ и, соответственно, снижения конечной дозы радиации. Траектории при каждом значении параметра продолжения представлены в осях “радиус–наклонение”. Две пунктирные линии на каждом графике показывают зависимости высоты перигея и апогея от наклонения на траектории оптимального быстродействия. Две непрерывные линии соответствуют зависимости высоты перигея и апогея от наклонения на траектории, полученной при данном значении параметра продолжения. Оттенками серого на каждом графике показана величина осредненной мощности дозы. Видно, что сначала наблюдается тенденция к увеличению максимальной высоты перигея, чтобы вынести перигейные участки траектории из области высокой мощности дозы. Затем примерно при $\tau \approx 0.19$ высота перигея становится равной высоте апогея и оскулирующая орбита становится круговой во время всего перелета. Далее по мере увеличения параметра продолжения максимальная высота траектории увеличивается, пока примерно при $\tau \approx 0.36$ снова на траектории появляется ненулевой эксцентриситет. Далее по мере увеличения параметра продолжения максимальный эксцентриситет растет вплоть до максимальных полученных значений τ. Также на рис. 6 можно заметить, что по мере снижения конечной дозы радиации изменение наклонения орбиты все меньше происходит в области пересечения траекторией ММБ радиационных поясов Земли и все больше – вне этой области. Этот факт частично подтверждает правильность расчетов, так как видно, что с уменьшением наклонения мощность дозы внутри радиационных поясов растет, и поэтому выгоднее в плане радиационной нагрузки пересекать РПЗ (другими словами набирать высоту) в области больших наклонений.

На рис. 7 изображены проекции двух траекторий на координатные плоскости инерциальной геоэкваториальной системы координат. Траектория слева является траекторией оптимального быстродействия. Справа изображена траектория с наименьшей полученной дозой для модели потоков частиц РПЗ AE8/AP8 MAX. Все проекции на рис. 7 изображены в одном масштабе. В центре каждой из проекций условно изображена Земля. Жирной линией в плоскости $Oxy$ показана геостационарная орбита.

Для траекторий, изображенных на рис. 7, на рис. 8 представлены графики зависимости углов тангажа и рысканья, а также дозы радиации, поглощенной ЧЭ, в зависимости от времени перелета ММБ на ГСО с помощью ЭРДУ. Три графика слева соответствуют траектории оптимального быстродействия, а графики справа относятся к траектории с наименьшей полученной дозой для модели потоков частиц РПЗ AE8/AP8 MAX. На каждом из графиков отмечены серым цветом участки траектории, когда поглощенная ЧЭ доза активно растет. Эти участки примерно соответствуют пересечению ММБ внутреннего и внешнего радиационных поясов. Как видно на графиках угла рысканья для траектории с наименьшей дозой амплитуда угла рысканья снижается внутри серых областей, что приводит к меньшему изменению наклонения внутри РПЗ, наблюдавшемуся на рис. 6.

ЗАКЛЮЧЕНИЕ

В работе были рассмотрены возможности оптимизации поглощенной дозы за счет изменения формы траектории выведения ММБ с низкой круговой орбиты на ГСО с использованием ЭРДУ и ЯЭУ. Математический аппарат построения оптимальных по быстродействию траекторий с фиксированной конечной дозой позволил снизить поглощенную за перелет дозу радиации на 25–38% в зависимости от используемой модели потоков частиц РПЗ. При этом время выведения увеличилось на 4–7% от оптимального времени выведения. Дополнительные затраты характеристической скорости увеличиваются относительно траектории оптимального быстродействия на 320–560 м/с. Для эффективного решения краевых задач принципа максимума, возникающих при решении задач оптимального быстродействия с фиксированной конечной дозой, был предложен метод аппроксимации мощности дозы с помощью осреднения и построения сплайнов 11-го порядка. Полученная аппроксимация является достаточно гладкой для применения в численном интегрировании уравнений движения ММБ интегратором 7(8) порядка и достаточно хорошо отражает характер зависимости накопленной дозы в зависимости от формы траектории выведения ММБ.

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

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