Исследование режима одноосной солнечной ориентации спутника

( Analysis of the single axis solar orientation of the Earth artificial satellite
Preprint, Inst. Appl. Math., the Russian Academy of Science)

ВВЕДЕНИЕ

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

Работа состоит из двух разделов. В первом разделе вычисляется действующий на солнечный парус момент сил светового давления со стороны Солнца. Рассмотрены симметричный парус с одинаковыми лепестками произвольной формы и парус в виде кусочно-линейной триангулированной поверхности. Полученные формулы применены для расчета момента сил светового давления, действующих на парус, который испытывает тепловые деформации из-за нагрева солнечными лучами.

Второй раздел посвящен численному интегрированию уравнений движения спутника относительно центра масс. Здесь использованы два типа уравнений такого рода, по-разному описывающих вращательное движение спутника, близкое регулярной прецессии Эйлера. Первый тип – достаточно подробные уравнения, включающие динамические уравнения Эйлера вращательного движения твердого тела; второй тип – так называемые эволюционные уравнения, описывающие вековое изменение (эволюцию) наиболее существенных параметров вращательного движения спутника. При расчете правых частей эволюционных уравнений использовано численное усреднение приложенного к спутнику момента сил светового давления. Такой подход делает использованную методику исследования движения универсальной – пригодной для любой модели паруса. В разных разделах использована независимая нумерация формул.

1.     ВЫЧИСЛЕНИЕ ДЕЙСТВУЮЩЕГО НА ПАРУС МОМЕНТА СИЛ

СВЕТОВОГО ДАВЛЕНИЯ

1.1. Общие формулы для случая симметричного паруса. Солнечный парус состоит из восьми одинаковых лепестков и представляет собой осесим­метричную поверхность, которая переходит сама в себя при поворотах вокруг оси симметрии на . Лепестки можно разбить на две одинаковые группы с такими же свойствами симметрии. Каждая группа получается из другой сдвигом вдоль общей оси симметрии и поворотом вокруг этой оси на угол . Вследствие симметрии вычисление приложенного к парусу момента сил светового давления естественно начать с вычисления такого момента, действующего на одну из указанных четверок лепестков.

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

,    ,    ;

,    ,   ;

,  ,  ;

,   ,    ;

причем каждый лепесток в собственной системе координат задается одним и тем же уравнением .

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

,

,   ,   .

Здесь  – световое давление на плоское идеальное зеркало, расположенное перпендикулярно солнечным лучам в точке ,  и  – площадь элемента поверхности лепестка в точке с радиусом-вектором  и нормаль к этому элементу. Направление нормали выбрано так, что . В собственной системе координат лепестка  нормаль имеет вид

,   .

Интегрирование в выписанных интегралах производится по всей поверхности паруса, и каждый интеграл естественным образом разбивается на сумму четырех интегралов, вычисляемых по поверхности отдельных лепестков. Интегралы по поверхности лепестка представляют собой векторы. Они записываются покомпонентно в собственной системе координат лепестка, преобразуются к удобной форме и затем пересчитываются в систему . В этой системе полученные векторные выражения складываются. Несложные, но громоздкие вычисления приводят к довольно простым (вследствие симметрии задачи) соотношениям. Пусть в системе координат

,    ,    .

Тогда

,

,    ,

,   ,   ,

,   ,

,

,   ,   .

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

,   ,

,   , ,   .

1.2. Частные случаи формы паруса. В частных случаях формы лепестка интегралы, выписанные в конце предыдущего пункта, легко вычисляются. Например, это имеет место в случае плоских лепестков [1], т.е. для линейной функции . В этом случае const  и рассматриваемые шесть интегралов выражаются через три интеграла

    ,                                        (1)

имеющих прозрачный геометрический смысл. Здесь – площадь лепестка,  – координаты его геометрического центра масс в системе . При , , , , , ,  получаем лепесток па­руса, рассматривавшегося в [1]. Для него , , , , . Составляющие момента сил светового давления, действующего на четверку таких лепестков, имеют вид

,   ,

где  – орт оси .

Рассмотрим вторую четверку тех же лепестков, но с . Остальные параметры сохраним прежними. Вторая четверка сдвинута относительно первой вдоль оси , и ее можно считать повернутой вокруг этой оси относительно первой четверки на угол  (с математической точки зрения поворот здесь возможен на любой угол; угол  и размеры лепестков выбраны из условия не затенения паруса). Такие две четверки образуют парус, рассматривавшийся в [1]. Составляющие  и  момента сил светового давления, приложенного ко второй четверке лепестков, выражаются указанными выше формулами, в которых следует сделать замену . Складывая моменты сил светового давления, действующие на обе четверки, получим формулу

,                                            (2)

.

Эта формула была получена в [1] другим, более элементарным способом. Из­менив одно принятое соотношение:   , из при­веденных выше общих выражений можно получить применявшуюся в [1] формулу момента сил светового давления для паруса с пропеллерным эффектом.

          Рассмотрим теперь другой частный случай. Пусть в собственной системе координат лепесток симметричен относительно плоскости . Это значит, проекция лепестка на плоскость  симметрична относительно оси  и . Из последнего соотношения следуют равенства

,  ,  .

Предположим для определенности, что проекция лепестка на плоскость  имеет вид: , . Тогда

.

Подынтегральное выражение во внутреннем интеграле повторного интеграла справа – нечетная функция  при любом , поэтому 

.

Отсюда . Аналогично устанавливаются равенства  и . В ре­зультате составляющие момента сил светового давления, действующего на четверку лепестков, имеют вид

,   .

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

          Приведем вычисление  для лепестка паруса, задаваемого функциями

,     ,

где ,  ,  ,  . Эти функции приближенно описывают лепесток в форме равнобедренного треугольника, изогнутого по поверхности прямого кругового цилиндра радиуса , причем  ось симметрии треугольника изгибается по окружности, перпендикулярной образующим цилиндра. Данный лепесток моделирует плоский лепесток, испытывающий тепловую деформацию (см. ниже). Если , то деформация отсутствует и такой лепесток является частным случаем плоского лепестка, рассмотренного выше, в частности, нормаль к лепестку имеет вид . В приводимых ниже выражениях опущены слагаемые, имеющие порядок . Для одного лепестка имеем

,

,   .

Полученные соотношения удобно выразить через величины ,  и , которые использовались выше при рассмотрении плоского лепестка. Чтобы установить соотношение между этими величинами и параметрами данного примера воспользуемся соотношениями (ср. выше)

,

,

.

Отсюда находим

,    ,    .

Воспользовавшись приведенными соотношениями, выражения для ,  и  можно представить следующим образом

,

,    .

Рассмотрим вторую четверку таких же лепестков, только заменим ее параметр  на . Остальные параметры оставим прежними. Вторая четверка сдвинута относительно первой вдоль оси , и ее можно считать повернутой вокруг этой оси относительно первой четверки на угол . Такие две четверки образуют парус из восьми лепестков. Вычисление действующего на этот парус момента сил светового давления сводится к суммированию (с нужными весами) компонент ,  для обеих четверок и приводит к формуле (2), в которой , причем  совпадает с прежним ,

.

          Приведем числовые оценки. На орбите Земли  Н/м. Па­раметры формы исходного плоского лепестка возьмем такие же, как в [1]: , , , м, м, м. В этом случае Нм. Для  рассмотрим два значения  м (соответствует коэффициенту теплового расширения трубы каркаса лепестка  1/град. и перепаду температур на освещенной и теневой сторонах трубы С) и м ( 1/град., перепад температур С). В первом случае  Нм и отношение . Во втором случае  Нм и отношение .

          1.3. Вычисление момента сил светового давления с учетом тепловой деформации паруса. Тепловая деформация лепестков паруса, вызываемая па­дающими на них солнечными лучами, в общем случае не одинакова из-за раз­ных условий их освещенности. Это обстоятельство делает парус несимметрич­ным и требует анализа приемлемости указанных выше формул для расчета момента сил светового давления. Чтобы оценить ошибки, вносимые в эти формулы тепловой деформацией паруса, воспользуемся следующими со­ображениями. Каждый лепесток паруса представляет собой перепонку из тонкой пленки, натянутую между двумя надутыми газом трубами, которые также изготовлены из тонкой пленки. Диаметр каждой трубы – 15 см, длина – более 14 м. Нагрев трубы солнечными лучами приводит к перепаду температур около K на ее освещенной и теневых сторонах. В результате теплового расширения освещенная сторона трубы становится длиннее ее теневой стороны, и труба деформируется.

В первом приближении каждую трубу можно считать изогнутой по дуге окружности, которая лежит в плоскости, содержащей продольную ось неде­формированного паруса . В обоснование этого допущения заметим, что подходящей окружностью можно аппроксимировать любую гладкую кривую и что ось  составляет малый угол с ортом . Радиус аппроксимирующей окружности находится так. Поскольку ось трубы практически перпендикулярна , примем, что все точки наиболее нагретой образующей трубы имеют одинаковую температуру , а все точки наиболее холодной образующей – одинаковую температуру . Примем также, что эти образующие и ось лежат в одной плоскости. Пусть  – искомый радиус,  – диаметр трубы,  – коэффициент линейного теплового расширения пленки. Из геометрического подобия наиболее нагретой и наиболее холодной образующих: а также из закона теплового расширения следует соотношение . Отсюда находим

.

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

Теперь усложним модель, и будем учитывать изгиб каждой трубы инди­видуально. Каждую трубу, точнее, ее ось аппроксимируем ломаной из  звеньев одинаковой длины . Полагаем, что эта ломаная лежит в плоскости, содержащей ось . Пронумеруем звенья ломаной числами 1, 2, …, , начиная от места крепления трубы к центральному телу спутника. Орт звена с номером  и нормаль к этому звену, лежащую в плоскости ломаной и оси , обозначим соответственно  и . Ось  и орты ,  компланарны. Направление нормали выберем из условия . Введем угол  положив

,        .

Как нетрудно видеть,  – угол между нормалью  и проекцией орта  на плоскость, содержащую ось  и орт . Для каждого звена ломаной введем температуры  и  наиболее нагретой и наиболее холодной образующих со­ответствующего куска трубы. Для трубы, изготовленной из пленки ПЭТ-М, эти температуры (выраженные в ) определяются системой уравнений

,

,

,   ,

,   .

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

.

Нет необходимости вычислять все введенные углы в явном виде. Можно воспользоваться формулами

,        .

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

          Вычисление момента сил светового давления, действующего на парус, выполняется так (ср. [2] ). Рассмотрим какой-либо лепесток и две ломаные, мо­делирующие его трубы. Звенья этих ломаных с одинаковыми номерами считаем противоположными сторонами пространственных четырехугольников (вершины четырехугольника в общем случае не лежат в одной плоскости). Другие стороны этих четырехугольников суть прямые, соединяющие вершины обеих ломаных. При этом первая от центрального тела вершина одной ломаной соединяется с аналогичной первой вершиной второй ломаной, вторая вершина первой ломаной соединяется со второй вершиной второй ломаной, третья с третьей и т.д. Каждый такой пространственный четырехугольник диагональю, идущей, если смотреть с положительного конца оси , из правой дальней его вершины в левую ближнюю, разбивается на два треугольника. В результате каждый лепесток паруса разбивается на треугольники. Вычисление момента сил светового давления на парус сводится к суммированию моментов, создаваемых всеми такими треугольниками. Расчетные формулы имеют вид

,                                          (3)

,   ,   .

Здесь  – площадь треугольника с номером ,  – радиус-вектор геометриче­ского центра масс этого треугольника (точки пересечения его медиан) относи­тельно точки ,  – нормаль к плоскости -ого треугольника, составляющая острый угол с осью , индекс  пробегает значения . Напомним,  – число звеньев ломаной, моделирующей трубу.

Результат расчета момента указанным довольно громоздким способом интересно сравнить с результатом расчета по формуле (2). Интересно даже не совпадение или отличие численных значений момента, а структура зависимости . Чтобы изучить эту структуру, запишем формулу (2) в координатной форме в системе . Орт  в этой системе зададим углами  и , положив . Получим

,   ,   .

Отсюда следует

,    .

Таким образом, чтобы проверить, сводится ли новая зависимость  к виду (3), следует проверить, равны ли тождественно нулю две функции и тож­дественна ли некоторой постоянной, отличной от нуля, третья функция. Все эти функции зависят от двух переменных. Проверку сведем к вычислению величин  и ,  , определенных соотношениями

,    ,

,    .

В общем случае эти величины – функции , но если с высокой точностью вы­полняются соотношения  , то можно считать справедливой формулу (2). Если дополнительно окажется, что const, то ситуация будет в точности той, которая рассматривалась в [1].

Для приближенного вычисления величин  и  , воспользу­емся формулой прямоугольников на равномерной сетке. В случае интегрирования периодической функции по периоду эта формула является весьма точной. Расчеты проводились для лепестка, имеющего в недеформируемом состоянии форму равнобедренного треугольника с углом при вершине  и высотой, опущенной на основание, 14.2 м. Диаметр трубы – 15 см,  1/град.,  и . В данном случае две верных значащих цифры при вычислении  можно получить при 10 и более точек в указанной формуле. Результаты расчетов приведены в пяти первых столбцах следующей таблицы.

Таблица 1

Как показывают значения величин  в таблице, формула (2) оказалась достаточно точной. Однако хотя эти величины малы, коэффициент  при изменении  от нуля до  меняется примерно на 1.3%. Таким изменением можно пренебречь, но интересно заметить, что найденные значения  достаточно точно аппроксимируются формулой , где  и  – постоянные коэффициенты. Определение этих коэффициентов методом наименьших квадратов дает  Нм,  Нм. Ошибки аппроксимации  приведены в последнем столбце таблицы.

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

Таблица 2

Приведенные результаты показывают, что в рамках принятой модели влияние тепловой деформации приводит в основном к увеличению модуля ко­эффициента . Новое значение этого коэффициента лежит в пределах, рас­смотренных в [1]. В [1] параметры  и в формуле (2) парметризовались углом : , . Рассматривались несколько значений этого угла из диапазона . При  имеем  Нм. Таким об­разом, учет тепловой деформации паруса не повлиял даже на количественные выводы работы [1]. Заметим, что точность одноосной солнечной ориентации при этом новом значении  увеличивается. Даже если учесть зависимость  от , то структура момента светового давления не изменится – он по-прежнему будет потенциальным. Это значит, методы и выводы работы [1] применимы и в рассматриваемой ситуации.

1.4. Интерполяция момента сил светового давления. Аппроксимация поверхности паруса системой треугольников и расчет действующего на парус момента сил светового давления по формулам (2) универсальны и практичны. Однако в случае паруса сложной формы или использования достаточно де­тальной модели паруса работа компьютерной программы, основанной на этих формулах, требует сравнительно больших затрат процессорного времени. В процессе интегрирования уравнений вращательного движения спутника желательно вообще не обращаться к такой программе. Следует один раз насчитать функцию  на некоторой сетке и хранить насчитанные данные в отдельном файле. Перед интегрированием уравнений движения эти данные следует считать в массив и в процессе интегрирования вычислять по ним  с помощью интерполяционных формул [2, 3].

          Интерполяционные формулы составим исходя из факта, что орт  в про­цессе ориентированного движения спутника не отклоняется от оси  более чем на . Пусть в системе координат  . Введем параметры , ,  и будем рассматривать  как функцию , заданную в квадрате . В плоскости  введем сетку с узлами в точках   и вычислим массив . Интерполяция функции  осуществляется так. Для произвольной точки  выбираются целые числа  и , удовлетворяющие неравенствам , , и находятся , . В качестве интерполяционного значения  принимается [4]

,

,     ,

,

,  ,  ,  ,  .

Функция  непрерывно дифференцируема, причем .

          В данной работе описанный способ интерполяции применялся только в случае, когда учитывалась тепловая деформация паруса. В остальных случаях формулы (3) использовались непосредственно.

1.5. Модели паруса, использованные в расчетах. При численном моделировании режима одноосной ориентации спутника рассматривались четыре варианта задания момента сил светового давления. В вариантах 1 – 3 лепестки паруса считались одинаковыми равнобедренными треугольниками с углом при вершине  и высотой, опущенной на основание,  м. Каждая такая высота лежит в плоскости, содержащей ось  и составляет с ней угол . Эти плоскости образуют между собой углы, кратные . Вершины треугольников, противоположные основанию, лежат на оси . При этом вершины двух соседних треугольников сдвинуты друг относительно друга на 0.5 м, и если треугольники занумеровать последовательно, то вершины всех треугольников с номерами одинаковой четности совпадают. Если основания всех треугольников параллельны плоскости , то такой парус в точности совпадает с парусом идеальной формы, рассматривавшемся в [1]. В данной работе такой вариант модели паруса называется вариантом 3. Полагаем при этом, что центр масс спутника совпадает с вершинами, имеющими меньшее значение координаты .

          Повернем каждый треугольник на одинаковый угол  вокруг высоты, опущенной на основание. Знак этого угла считаем положительным, если для наблюдателя на вершине треугольника поворот выполняется против часовой стрелки. Полученный парус представляет собой парус с пропеллерным эф­фектом из работы [1]. Этот вариант модели паруса назовем вариантом 1.

          Вариант 2 получается из варианта 3, если в парусе идеальной форме повернуть на одинаковые по модулю, но противоположные по знаку углы только два противоположных лепестка одинаковой четности. Для определенности лепестки паруса занумеруем числами 1, 2, …, 8 и будем считать, что вершины лепестков 1, 3, 5, 7 (напомним, они совпадают и лежат на оси ) имеют координату больше той же координаты вершин лепестков 2, 4, 6, 8. Поворачиваются лепестки 1 и 5. Угол поворота основания лепестка 1 обозначим . Тогда угол поворота лепестка 5 будет . Как нетрудно видеть, при  варианты 1 и 2 совпадают с вариантом 3. Вариант 2 рассматривается для иллюстрации того факта, что повороты лепестков паруса вокруг высоты, опущенной на основание, могут не вызывать пропеллерного эффекта.

          Вариант 4 – это по существу вариант 3, подверженный описанной выше тепловой деформации. Детальное описание этого варианта приведено в п. 1.3. Деформация труб рассчитывалась при аппроксимации их ломаными с 40 звеньями (. Для расчета момента светового давления при интегрировании уравнений движения спутника использовался способ интерполяции описанный в п. 1.4. Интерполяция строилась при , погрешность интерполяции не превышала Нм.

2.     ОДНООСНАЯ СОЛНЕЧНАЯ ОРИЕНТАЦИЯ СПУТНИКА

1. Уравнения движения спутника относительно центра масс. Будем использовать две системы таких уравнений. Первая система предназначена для непосредственного численного интегрирования, вторая служит отправной точкой при получении с помощью метода усреднения так называемых эволюционных уравнений, которые описывают вековое изменение наиболее существенных параметров вращательного движения спутника. Расчет правых частей эволюционных уравнений и их интегрирование выполняются численно. Возможность применения этих уравнений обусловлена специальным видом изучаемых движений.

Спутник считаем динамически симметричным твердым телом. Для записи указанных систем его уравнений движения введем три правые декартовы системы координат. Система  образована главными центральными осями инерции спутника. Тензор инерции спутника в этой системе имеет матрицу diag. Оси системы  сохраняют неизменные направления в абсолютном пространстве. Плоскость  параллельна плоскости земного экватора, оси  и  направлены в точку весеннего равноденствия и северный полюс мира. Система  связана с вектором  кинетического момента спутника в его движении относительно центра масс; ось  направлена по , ось  лежит в плоскости .

                Положение системы координат  относительно системы  будем задавать при помощи углов  и . Первый из этих углов образован осями  и , второй – осями  и . Направления отсчета этих углов согласованы с направлениями осей  и  соответственно. Матрицу перехода от системы  к системе  обозначим , где  – косинус угла между осями  и . Элементы этой матрицы имеют вид

,    ,    ,

,    ,     ,

,            ,          

                Положение системы координат  относительно системы  зададим углами ,  и , определив их следующим образом. Система  может быть переведена в систему  тремя последовательными поворотами: 1) на угол  вокруг оси , 2) на угол  вокруг новой оси , 3) на угол  вокруг новой оси , совпадающей с осью . Матрицу перехода от системы  к системе  обозначим , где  – косинус угла между осями  и . Элементы этой матрицы выражаются формулами

,        ,

,                 ,

,     ,

,

Матрицу перехода от системы  к системе  обозначим ,  – косинус угла между осями  и . Введенные матрицы связаны между собой формулой .

Пусть  – вектор , ,  – его компоненты в одной из введенных систем координат. Если эти компоненты относятся к системе , то будем писать . Если же они указаны в системе  или , то соответственно напишем  или  .

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

,      ,

,                           (1)

,        ,

,        ,

,        .

Здесь  – время,  – абсолютная угловая скорость спутника,  – гравитационный параметр Земли,  – геоцентрический радиус-вектор точки , ,  – действующий на спутник момент сил светового давления. Гравитационный момент в выписанных уравнениях задан явными выражениями.

Первые три уравнения (1) представляют собой динамические уравнения Эйлера вращательного движения твердого тела, остальные уравнения этой системы – кинематические уравнения Пуассона для направляющих косинусов осей  и  в системе координат . В силу последнего замечания переменные  и  не являются независимыми – они связаны соотношениями (некоторые из условий ортогональности матрицы  )

,     ,     ,

являющимися интегральными соотношениями системы (1). Эти соотношения должны быть выполнены при задании начальных условий вычисляемых решений.

В процессе интегрирования уравнений (1) величины  рассчитывались по формулам ,  и , движение точки  задавалось формулами задачи двух тел, в которых учитывались вековые изменения элементов орбиты, обусловленные главным членом сжатия Земли. Необходимое для расчета момента  геоцентрическое движение Солнца находилось по формулам [5]. При этом в качестве орта  направления “спутник – Солнце” (см. предыдущий раздел) использовался орт “Земля – Солнце”.

Переходим ко второй системе. Для реализации одноосной солнечной ориентации спутника на длительных интервалах времени предлагается использовать движения с достаточно большим собственным кинетическим моментом [1]. Эти движения близки регулярной прецессии Эйлера осесимметричного твердого тела, и при изучении этих движений внешние моменты можно рассматривать как малые возмущения. Приближенное исследование движения спутника в такой ситуации можно провести методом усреднения. Такое исследование проведено в [1], но там принято, что парус имеет номинальную форму и для момента сил светового давления использованы простые аналитические выражения. В частности, в большинстве рассмотренных в [1] задач момент светового давления оказался потенциальным. Основное внимание в [1] уделено получению приближенных аналитических оценок, характеризующих точность одноосной солнечной ориентации и время ее существования.

В случае произвольной формы солнечного паруса детальное исследование движения спутника можно провести только численными методами. Непосредственное численное интегрирование уравнений (1) позволяет решить эту задачу, но для надежности полученные таким способом результаты следует подтвердить другим методом. В качестве альтернативного метода использовалось численное интегрирование так называемых эволюционных уравнений. Эволюционные уравнения выводятся в рамках той же схемы усреднения по регулярной прецессии Эйлера [6], которая применялась в [1], но записаны они в других переменных, имеют более высокий порядок и дают возможность изучить более широкий класс движений. Использованная в данной работе методика численного расчета правых частей этих уравнений является частным случаем методики [2,3].

Пусть ,  – вычисленный относительно точки  главный момент внешних сил, действующих на спутник. Уравнения движения спутника вокруг центра масс запишем в виде

,     ,     ,

,                (2)

,

.

Вывод первых трех из этих уравнений приведен в [6], уравнения для углов ,  и  можно вывести по схеме, которая использованна в [6] для вывода аналогичных уравнений для углов Эйлера, задающих положение системы координат  относительно системы .

В общем случае использованные в [6] углы Эйлера более удобны для описания регулярной прецессии твердого тела и вывода эволюционных уравнений. Однако при совпадении оси собственного вращении спутника с вектором его кинетического момента (осью прецессии) эти углы вырождаются, и именно такие движения планируется использовать для реализации режима одноосной солнечной ориентации [1]. Указанное вырождение делает невозмож-ным применение углов Эйлера в численных расчетах представляющих интерес движений. Чтобы избежать вырождения, в данной работе введены углы ,  и , хотя это несколько усложнило вывод эволюционных уравнений.

Поясним смысл переменных системы (2). Величины ,  и  задают вектор  в абсолютном пространстве, углы ,  и  характеризуют движение спутника относительно этого вектора. Главный момент  внешних сил полагаем суммой гравитационного момента и момента сил светового давления. Тогда   суть функции , элементов и  введенных матриц перехода, а также положений Солнца и спутника в системе координат . Геоцентрические движения Солнца и спутника будем задавать так же, как в уравнениях (1). В результате придем к функциям , периодически зависящим от угла  с периодом . Получили замкнутую систему уравнений вращательного движения.

2.2. Эволюционные уравнения. Влияние внешних моментов на вращательное движение спутника будем считать малым. Это допущение позволяет воспользоваться для приближенного интегрирования уравнений (2) методом усреднения. Сначала приведем основные формулы невозмущенного движения.

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

,     ,                               (3)

,                                    (4)

в которых  – параметр. Уравнения (3) не содержат  и образуют замкнутую подсистему с первым интегралом const. Замена переменных ,  переводит систему (3) в систему

,       ,

общее решение которой запишем в виде

,       ,       ,       .

Здесь  и  – произвольные постоянные, . Задание переменных  и  не определяет углов  и  единственным образом. Ниже будем рассматривать движения, в которых эти углы близки к нулю. В этом случае указанный выше первый интеграл системы (3) принимает вид .

          Подставив найденное решение системы (3) в уравнение (4), получим

.

Поскольку , из последнего соотношения следует

.

Полученное выражение для  можно представить следующим образом

,

где  – произвольная постоянная,

,

,     ,

.

При  выписанное решение принимает вид ,  и описывает стационарное вращение спутника вокруг оси , совпадающей с осью .

          Регулярную прецессию спутника будем задавать величинами , ,  и . В случае  изменение  описывается уравнением

.                          (5)

При переменные , ,  и  являются медленными, переменные ,  и  – быстрыми. Зависимость компонент внешнего момента   от времени при фиксированных значения остальных переменных выражается квазипериодическими функциями времени, наибольшая частота которых равна орбитальной частоте (среднему движению спутника) . Поскольку для представляющих интерес движений , , время в рассматриваемой задаче также является медленной переменной.

Движение спутника можно рассматривать как медленную эволюцию регулярной прецессии. Дифференциальные уравнения, приближенно описы-вающие в нерезонансном случае вековое изменение (эволюцию) медленных переменных, так называемые эволюционные уравнения [6], получаются усреднением правых частей первых трех уравнений (2) и уравнения (5) по регулярной прецессии Эйлера при фиксированных значениях всех медленных переменных. Усреднение проведем следуя рецепту [7], первоначально предло-женному для усреднения по более сложному движению Эйлера – Пуансо. С помощью указанного выше решения системы (3) тригонометрические функции углов  и  выразим через величины  и , элементы матрицы перехода  выразим через ,  и . В результате компоненты внешнего момента станут функциями величин , , , , ,  и , периодическими по  и  с периодом . Уравнение (5) при этом перейдет в уравнение

.

В нерезонансном случае (при иррациональном отношении частот  и ) эволюционное уравнение для  получается усреднением правой части последнего уравнения по  и

.

Аналогичным образом получаются эволюционные уравнения относительно ,  и . Операцию усреднения, определяемую правой частью последнего уравнения, будем обозначать символом . С использованием этого обозначения эволюционные уравнения можно записать так

 ,     ,     ,           (6)

Формально условие отсутствия резонанса означает, что отношение  должно быть иррациональным числом. На самом деле эволюционные уравнения (6) могут быть неадекватны лишь в тех областях изменения переменных  и , в которых это выражение приближенно равно отношению небольших натуральных чисел.

Поскольку правые части уравнений (6) линейны по моменту , при их вычислении операцию усреднения  можно выполнить над каждым учитываемым внешним моментом по отдельности, а затем полученные результаты сложить. Для некоторых моментов, например, для гравитационного, усреднение  быть проведено аналитически [6], но в общем случае момента сил светового давления приходится применять численные методы. Указанное свойство линейности позволяет выбирать при усреднении каждого момента наиболее удобный способ.

В случае рассматриваемого спутника с солнечным парусом слагаемые правых частей уравнений (6), отвечающие гравитационному моменту, рассчитывались по формулам [6]

,         ,

,       ,

где ; слагаемые, отвечающие моменту сил светового давления, находились численно.

Опишем численное выполнение операции . Оно заключается в вычислении двойного интеграла по квадрату  от функции двух переменных, причем по каждой переменной эта функция – периодическая с периодом . При вычислении таких интегралов в случае достаточно гладкой подынтегральной функции весьма эффективны теоретико-числовые квадратур- ные формулы [8]. Для функции двух переменных оптимальная по точности формула выглядит следующим образом. Пусть функция  – периодическая по  и  с периодом 1,  и  – соседние члены последовательности чисел Фибоначчи 1, 1, 2, 3, 5, …Тогда принимается

,                            (7)

где символ  означает дробную часть числа. В расчетах, результаты которых описаны ниже, результат практически не зависел от  при .

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

Здесь  – орт орбитального кинетического момента спутника (нормаль к плоскости орбиты),  – орт направления “Земля – Солнце”,  км – радиус Земли,  – большая полуось орбиты спутника. Коэффициент  представляет собой отношение длины отрезка орбиты, освещенного Солнцем, к длине всей орбиты. При выводе выражения для  орбита спутника считалась круговой радиуса . В данном случае это допущение оправдано, поскольку эксцентриситет орбиты мал.

Слагаемые правых частей уравнений (6), отвечающие гравитационному моменту и усредненные по орбитальному моменту, рассчитывались по формулам [2]

,

,

где  – эксцентриситет орбиты.

Ниже, говоря о системе (6), усредненной по орбитальному движению, будем называть ее системой (6). Именно система (6) имелась в виду, когда в начале данного раздела упоминались эволюционные уравнения, выводимые из второй системы. Система (6) использовалась только для проверки программы численного интегрирования системы (6).

2.3. Результаты численных расчетов. Исследование движения спутника в режиме одноосной солнечной ориентации сводилось к численному интегрированию уравнений (1), (6) и (6), причем основное внимание уделялось системам (1) и (6).. Начальный момент  вычисляемых решений был отнесен к  декретного московского (зимнего) времени 22.09.2001. Принималось, что спутник в этот момент находился в восходящем узле орбиты. Параметры орбиты были взяты следующие (ср. [1]): большая полуось –  км, эксцентриситет – , наклонение – , начальное значение аргумента широты перигея – . Начальное значение долготы восходящего узла – . Значения моментов инерции спутника кгм, кгм. Момент сил светового давления рассчитывался для четырех вариантов модели паруса, описанных в разделе 1.

Примем, что в начальный момент времени ось  спутника направлена по орту , ось  лежит в плоскости  (положение оси  при  может быть произвольным, указанное положение выбрано из соображений удобства). В этом случае начальные значения элементов  имеют вид

,  ,  ,  ,  ,  ,

,  ,  ,  .

Начальные значения угловой скорости: , , . Выписанные соотношения задают начальные условия для системы (1). Начальные условия для систем (6) и (6) находились из соотношений

   ,     ,       (8)

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

          Сначала опишем результаты интегрирования систем (1), (6) и (6) с согласованными начальными условиями, чтобы оценить точность систем (6) и (6) Вычисляемые решения будем представлять графиками зависимости от времени величин , ,  и . В случае системы (1) эти величины находились из соотношений (8) на каждом шаге численного интегрирования. Полученные результаты представлены на рис. 1 – 8. На рис. 1 – 4 в одних и тех же координатных осях изображены графики решений системы (6) и графики, построенные по решениям системы (1). Как видно из этих графиков, решения обеих систем достаточно точно согласуются между собой на интервале времени около полутора суток. Точность согласования не зависит от использованной модели паруса.

На рис. 5 – 8 в одних и тех же координатных осях приведены графики решений систем (6) и (6) на интервале времени 2 – 3 недели. Эти решения достаточно хорошо совпадают по переменным ,  и вековому изменению амплитуд углов  и , но фазы двух последних углов в конце интервала интегрирования в этих решениях заметно отличаются. Полученные результаты согласуются с типичными оценками погрешности метода усреднения. Анализ рис. 5 – 8 свидетельствует о приемлемой точности систем (6) и (6) на интервалах времени в несколько десятков суток. Особенно удобна в этом отношении система (6), которая позволяет вычислять движение спутника с весьма малыми затратами процессорного времени и оперативно проводить анализ большого числа вариантов движения.

Рис. 1 и 5 иллюстрируют сильноое влияние пропеллерного эффекта на движение спутника. Чтобы этот эффект не очень сильно проявлялся на длительных интервалах времени, в первой модели паруса на длительных интервалах принималось . Рис. 2 и 6 показвают, что сбалансирован- ный поворот лепестков паруса (вариант 2) не вызывает пропеллерного эффекта.

Перейдем к результатам исследования движения спутника на длительных интервалах времени. Представление решений системы (6) на интервале времени в 60 сут не требует каких-либо дополнительных усилий. С решениями системы (1) дело обстоит сложнее. Решения системы (1) на продолжительных интервалах времени будем представлять следующим образом. Введем обозначения:  – драконический период спутника (промежуток времени между двумя его последовательными прохождениями восходящего узла орбиты, в данном случае const),   – временные отрезки, соответствующие орбитальным виткам спутника,  – угол между осью  и ортом ,  – угол между ортами  и . Введенные углы можно найти по формулам

,     .

На решении системы (1) определим последовательности чисел

,     ,     ,     ,

,      ,     ,     ,

,     ,           ().

Члены этих последовательностей суть минимальные и максимальные на каждом орбитальном витке значения некоторых переменных и функций. Введенные последовательности будем представлять графически ломаными. Например, последовательности  и  будем изображать ломаными, проходящими соответственно через точки  и   в плоскости  и т. д. Такие ломаные наглядно и сжато характеризуют границы изменения выбранных величин на каждом витке и позволяют получить представление о движении спутника.

Примеры указанных ломаных приведены на рис. 9, 11, 13 и 15. На рис. 10, 12, 14 и 16 представлены графики описывающих те же движения спутника решений системы (6). Для удобства сравнения разных способов описания движения спутника при использовании одной и той же модели паруса рисунки разбиты на пары. На одной странице рядом помещены ломаные, построенные по решениям системы (1), и графики решений системы (6). Налицо достаточно точное (но несколько хуже, чем на рис. 1 – 9) совпадение результатов, полученных с использованием разных систем. На некоторых графиках ломаные сливаются друг с другом и выглядят как гладкие кривые. В решениях на рис. 11 и 13 последовательности  и  при  заметно изменяют свое поведение. Это связано с появлением на орбите спутника теневого участка. Заметим, что с увеличением угловой скорости спутника эффективность (процессорное время, точность и т. п.) использования системы (1) уменьшается, а эффективность использования системы (6) увеличивается.

                Полученные результаты показывают, что использование эволюционных уравнений для исследования движения спутника в режиме одноосной солнечной ориентации приемлемо по точности. Показано также, что усложнение модели паруса не приводит к каким-либо новым эффектам в динамике этого режима и выводы работы [1] сохраняют свою силу.

          Данная работа выполнена при финансовой поддержке РФФИ (грант 02-01-00323).

Литература

1.     Лихачев В.Н., Сазонов В.В., Ульяшин А.И. Режим одноосной солнечной ориентации искусственного спутника Земли. Препринт ИПМ им. М.В.Келдыша РАН, № 15, 2001.

2.     Комаров М.М., Сазонов В.В. Расчет сил и моментов светового давления, действующих на астероид произвольной формы. Астрономический вестник, 1994, т. 28, № 1, с. 21- 30.

3.     Сазонов В.В. Движение астероида относительно центра масс под действием момента сил светового давления. Астрономический вестник, 1994, т. 28, № 2, с. 95- 107

4.     Maude A.D. Interpolation – mainly for graph plotters. Computer Journal, 1973, v. 16, № 1, p. 64-65.

5.     Меес Ж. Астрономические формулы для калькуляторов. М., Мир, 1988.

6.     Белецкий В.В. Движение искусственного спутника относительно центра масс. М., Наука, 1965.

7.     Черноусько Ф.Л. О движении спутника относительно центра масс под действием гравитационных моментов. Прикладная математика и механика. 1963, т. 27, № 3, с. 474-483.

8.     Коробов Н.М. Теоретико-числовые методы в приближенном анализе. М., Физматгиз, 1963.