Журнал вычислительной математики и математической физики, 2020, T. 60, № 10, стр. 1777-1786

ВВЕДЕНИЕ

Исследование движения спутника в центральном ньютоновом силовом поле по круговой орбите под действием гравитационного и аэродинамического моментов представляет значительный практической интерес для создания систем управления ориентацией искусственных спутников Земли. Принцип работы гравитационных систем ориентации основан на том, что в центральном ньютоновом поле сил спутник с неравными главными центральными моментами инерции имеет на круговой орбите 24 положения равновесия, четыре из которых являются устойчивыми [1]–[3]. Подробное рассмотрение динамики спутников с гравитационными системами ориентации представлено в [4]. Важным свойством гравитационных систем ориентации является возможность функционировать на орбите продолжительное время без расходования энергии и (или) рабочего тела. На орбитах в диапазоне высот от 250 до 700 км на ориентацию спутника, помимо гравитационного момента, существенное влияние также оказывает и аэродинамический момент. Поэтому необходимо исследовать совместное влияние гравитационного и аэродинамического момента на движение спутника относительно центра масс, в частности, существующие в данном случае положения равновесия и условия их устойчивости. Обзор основных результатов исследования движения спутников с аэродинамической системой ориентации представлен в [4]–[6].

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

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

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

Для записи уравнений движения введем две правые прямоугольные системы координат с началом в центре масс $O$ спутника. В орбитальной системе координат $OXYZ$ ось $OZ$ направлена вдоль радиуса-вектора, соединяющего центры масс Земли и спутника, ось $OX$ направлена вдоль вектора линейной скорости центра масс $O$ спутника. Тогда ось $OY$ будет направлена вдоль нормали к плоскости орбиты. В связанной со спутником системе координат $Oxyz$ оси $Ox$, $Oy$, $Oz$ – главные центральные оси инерции спутника.

Определим ориентацию системы координат $Oxyz$ относительно орбитальной системы координат с использованием самолетных углов тангажа $\alpha $, рыскания $\beta $ и крена $\gamma $. Направляющие косинусы осей $Ox,\;\,Oy,\;\,Oz$ в орбитальной системе координат выражаются через самолетные углы с помощью соотношений [4]:

При малых колебаниях спутника углу тангажа соответствует поворот вокруг оси $OY$, углу рыскания – поворот вокруг оси $OZ$, углу крена – поворот вокруг оси $OX$.

При выводе уравнений движения примем следующие предположения [4].

1) Действие атмосферы на спутник сводится к силе сопротивления, приложенной в центре давления и направленной против скорости центра масс спутника относительно воздуха; центр давления расположен на оси $Ox$ спутника. Данное предположение достаточно точно выполняется для формы спутника, близкой к сферической.

2) Влияние атмосферы на поступательное движение спутника пренебрежимо мало.

3) Увлечением атмосферы вращающейся Землей пренебрегается.

Пусть на спутник, кроме аэродинамического момента, действуют моменты активного демпфирования, проекции суммарного вектора которых на оси $Ox$, $Oy$, $Oz$ соответственно равны M_x = ${{\bar {k}}_{1}}{{p}_{1}},$ ${{M}_{y}} = {{\bar {k}}_{2}}({{q}_{1}} - {{\omega }_{0}})$ ${{M}_{z}} = {{\bar {k}}_{3}}{{r}_{1}}$. Здесь ${{\bar {k}}_{1}}$, ${{\bar {k}}_{2}}$, ${{\bar {k}}_{3}}$ – коэффициенты демпфирования, $p{}_{1}$, ${{q}_{1}}$, ${{r}_{1}}$ – проекции угловой скорости спутника на оси $Ox$, $Oy$, $Oz$; ${{\omega }_{0}}$ – угловая скорость движения центра масс спутника по круговой орбите. Тогда уравнения движения спутника относительно центра масс запишутся в виде

В уравнениях (2), (3) величины $A$, $B$, $C$ – главные центральные моменты инерции спутника; ${{H}_{1}} = - {{Qa} \mathord{\left/ {\vphantom {{Qa} {\omega _{0}^{2}}}} \right. \kern-0em} {\omega _{0}^{2}}},$ $Q$ – действующая на спутник сила сопротивления; ($a$, $0$, $0$) – координата центра давления спутника в системе координат $Oxyz.$ Для аэродинамически устойчивой конструкции спутника центр давления лежит за его центром тяжести и, следовательно, $a < 0$. Штрихом обозначено дифференцирование по времени $t$.

Если ввести безразмерные параметры ${{\theta }_{A}} = A{\text{/}}B$, ${{\theta }_{C}} = C{\text{/}}B$, $p = {{p}_{1}}{\text{/}}{{\omega }_{0}}$, $q = {{q}_{1}}{\text{/}}{{\omega }_{0}}$, $r = {{r}_{1}}{\text{/}}{{\omega }_{0}}$, ${{k}_{1}} = {{\bar {k}}_{1}}{\text{/}}\omega {}_{0}B$, ${{k}_{2}} = {{\bar {k}}_{2}}{\text{/}}\omega {}_{0}B$, ${{k}_{3}} = {{\bar {k}}_{3}}{\text{/}}\omega {}_{0}B$, ${{h}_{1}} = {{H}_{1}}{\text{/}}B$ и от времени $t$ перейти к безразмерной независимой переменной $\tau = {{\omega }_{0}}t$, то система (2), (3) примет вид

В уравнениях (4), (5) точкой обозначено дифференцирование по переменной $\tau $.

2. ПОЛОЖЕНИЯ РАВНОВЕСИЯ СПУТНИКА

Положив в (4) и (5) $\alpha = {{\alpha }_{0}} = {\text{const}}$, $\beta = {{\beta }_{0}} = {\text{const}}$, $\gamma = {{\gamma }_{0}} = {\text{const}}$, получим при $A \ne B \ne C$ уравнения

позволяющие определить положения равновесия спутника в орбитальной системе координат. Здесь ${{\tilde {k}}_{1}} = {{{{k}_{1}}} \mathord{\left/ {\vphantom {{{{k}_{1}}} {({{\theta }_{C}} - 1}}} \right. \kern-0em} {({{\theta }_{C}} - 1}})$, ${{\tilde {k}}_{2}} = {{{{k}_{2}}} \mathord{\left/ {\vphantom {{{{k}_{2}}} {({{\theta }_{A}} - {{\theta }_{C}}}}} \right. \kern-0em} {({{\theta }_{A}} - {{\theta }_{C}}}})$, ${{\tilde {k}}_{3}} = {{{{k}_{3}}} \mathord{\left/ {\vphantom {{{{k}_{3}}} {(1 - {{\theta }_{A}}}}} \right. \kern-0em} {(1 - {{\theta }_{A}}}})$ и $h = {{h}_{1}}{\text{/}}(1 - {{\theta }_{C}})$, ${v} = {{(1 - {{\theta }_{A}})} \mathord{\left/ {\vphantom {{(1 - {{\theta }_{A}})} {(1 - {{\theta }_{C}})}}} \right. \kern-0em} {(1 - {{\theta }_{C}})}}$.

Систему (6) можно рассматривать как систему трех уравнений с неизвестными ${{\alpha }_{0}}$, ${{\beta }_{0}}$, ${{\gamma }_{0}}.$ Другой способ замыкания уравнений (6) заключается в добавлении трех условий ортогональности направляющих косинусов:

Уравнения (6) и (7) образуют замкнутую алгебраическую систему уравнений относительно шести направляющих косинусов, определяющих положения равновесия спутника. Для этой системы уравнений ставится следующая задача: при заданных $h$, ${v}$, ${{\tilde {k}}_{1}}$, ${{\tilde {k}}_{2}}$, ${{\tilde {k}}_{3}}$ требуется определить девять направляющих косинусов, т.е. все положения равновесия спутника в орбитальной системе координат. После нахождения шести направляющих косинусов ${{a}_{{21}}}$, ${{a}_{{22}}}$, ${{a}_{{23}}}$, ${{a}_{{31}}}$, ${{a}_{{32}}}$, ${{a}_{{33}}}$ оставшиеся направляющие косинусы ${{a}_{{11}}}$, ${{a}_{{12}}}$, ${{a}_{{13}}}$ можно определить из условий ортогональности.

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

Для частного случая, когда ${{k}_{1}} = {{k}_{2}} = {{k}_{3}} = k$, исследование условий существования различного числа положений равновесия спутника при действии аэродинамического момента можно провести аналитически, применяя алгоритм построения базисов Гребнера [7] с использованием систем компьютерной алгебры. Метод построения базисов Гребнера представляет собой алгоритмическую процедуру для полного приведения задачи в случае полиномов многих переменных к рассмотрению полинома одной переменной. В этом случае можно построить базис Гребнера полиномов (6) и (7) относительно лексикографического упорядочения переменных. В полученном базисе Гребнера существует биквадратный полином, зависящий только от одной переменной ${{a}_{{22}}},$ корни которого позволяют определить все равновесия спутника для данного частного случая. Оставшиеся направляющие косинусы системы (6), (7) ${{a}_{{21}}}$, ${{a}_{{23}}}$, ${{a}_{{31}}}$, ${{a}_{{32}}}$, ${{a}_{{33}}}$ можно получить, приравняв соответствующие полиномы, входящие в базис Гребнера, нулю и найти корни получившихся алгебраических уравнений. Коэффициенты этих уравнений зависят от трех параметров ${v}$, $k$, $h$ и имеют громоздкий вид [8].

В настоящей работе главное внимание уделено исследованию условий асимптотической устойчивости нулевого положения равновесия ${{\alpha }_{0}} = 0$, ${{\beta }_{0}} = 0$, ${{\gamma }_{0}} = 0$, когда оси связанной со спутником системы координат совпадают с осями орбитальной системы координат.

3. ИССЛЕДОВАНИЕ НЕОБХОДИМЫХ И ДОСТАТОЧНЫХ УСЛОВИЙ АСИМПТОТИЧЕСКОЙ УСТОЙЧИВОСТИ ПОЛОЖЕНИЙ РАВНОВЕСИЯ СПУТНИКА

Для исследования необходимых и достаточных условий асимптотической устойчивости положений равновесия системы (6), (7) проведем линеаризацию уравнений системы (4), (5) в окрестности равновесной ориентации спутника $\alpha = {{\alpha }_{0}}$, $\beta = {{\beta }_{0}}$, $\gamma = {{\gamma }_{0}}$, удовлетворяющей системе (6), (7). Представим углы $\alpha $, $\beta $, $\gamma $ как $\alpha = {{\alpha }_{0}} + \bar {\alpha }$, $\beta = {{\beta }_{0}} + \bar {\beta }$, $\gamma = {{\gamma }_{0}} + \bar {\gamma }$ где $\bar {\alpha }$, $\bar {\beta }$, $\bar {\gamma }$ – малые отклонения от положения равновесия спутника. Тогда линеаризованные в окрестности положения равновесия $\alpha = {{\alpha }_{0}}$, $\beta = {{\beta }_{0}}$, $\gamma = {{\gamma }_{0}}$ уравнения системы (4), (5) могут быть записаны в виде

Рассмотрим малые колебания спутника в окрестности нулевого решения

Тогда из уравнений (8) получим систему Характеристическое уравнение системы (10) представляет собой произведение полиномов второй и четвертой степени относительно $\lambda $где

Необходимые и достаточные условия асимптотической устойчивости (условия Рауса–Гурвица) нулевого решения (9) имеют следующий вид:

Рассмотрим частный случай, когда ${{k}_{1}} = {{k}_{2}} = {{k}_{3}} = k$. Тогда условия (12) примут более простой вид:

Исследование эволюции областей, где выполняются необходимые и достаточные условия устойчивости (13), проводилось в плоскости $({{\theta }_{A}},\;{{\theta }_{C}})$ в зависимости от значений параметров $k$ и ${{h}_{1}}$. Необходимо принять во внимание неравенства треугольника для реальных тел, которым должны удовлетворять параметры ${{\theta }_{A}}$ и ${{\theta }_{C}}$: $1 + {{\theta }_{A}} \geqslant {{\theta }_{C}}$, $1 + {{\theta }_{C}} \geqslant {{\theta }_{A}}$, ${{\theta }_{A}} + {{\theta }_{C}} \geqslant 1$. Условия из неравенства треугольника выделяют в плоскости $({{\theta }_{A}},\;{{\theta }_{C}})$ бесконечную полуполосу. В силу неравенства ${{h}_{1}} + 3({{\theta }_{A}} - {{\theta }_{C}}) > 0$ из условий (13) размер полуполосы изменяется в зависимости от величины ${{h}_{1}}$ и при ${{h}_{1}} = 3$ будет совпадать с первым условием неравенства треугольника. Таким образом, исследование существования выполнения условий (13) проводилось в полуполосе, ограниченной прямыми

На фиг. 1 показано изменение границ областей, где выполняются условия асимптотической устойчивости (13) в зависимости от значений параметров $k$ и ${{h}_{1}}$. Области, где выполняются необходимые и достаточные условия устойчивости (13), выделены серым цветом.

При малых значениях $k$ и ${{h}_{1}}$, (фиг. 1а, $k = 0.2$, ${{h}_{1}} = 0.1$), область устойчивости приближается к области, где выполняются необходимые и достаточные условия устойчивости для спутника–твердого тела при $k = 0$, ${{h}_{1}} = 0$. Эта область ограничена прямыми ${{\theta }_{C}} = 1 - {{\theta }_{A}}$, ${{\theta }_{C}} = {{\theta }_{A}}$, ${{\theta }_{A}} = 1$ (см. [1]).

На фиг. 1б при $k = 0.5$, ${{h}_{1}} = 0.5$ подробно показаны кривые из условий (13), которые образуют границы области асимптотической устойчивости решения (9). Кривые ${{\Delta }_{4}} = 0$ представляют собой гиперболы и условие ${{\Delta }_{4}} > 0$ выполняется в области, расположенной между ветвями этой гиперболы.

При возрастании значений параметра ${{h}_{1}}$ происходит качественное изменение области выполнения условий устойчивости. Область выполнения условий (13) распадается на две подобласти (фиг. 1в, $k = 0.5$, ${{h}_{1}} = 1.5$). При ${{h}_{1}} = 3$, $k = 1.0$ (фиг. 1г) область выполнения условий (13) представляет собой два криволинейных треугольника, ограниченных прямыми ${{\theta }_{C}} = 1 - {{\theta }_{A}}$, ${{\theta }_{C}} = {{\theta }_{A}} - 1$ и верхней ветвью гиперболы ${{\Delta }_{4}} = 0$, прямой ${{\theta }_{C}} = {{\theta }_{A}} - 1$ и кривыми ${{\Delta }_{3}} = 0$ и ${{\Delta }_{4}} = 0.$ При дальнейшем увеличении параметра ${{h}_{1}}$ вид областей выполнения необходимых и достаточных условий асимптотической устойчивости (13) решения (9) практически не изменяется и остается близким к форме описанных выше криволинейных треугольников.

4. АНАЛИЗ ПЕРЕХОДНЫХ ПРОЦЕССОВ

Проведем теперь анализ переходных процессов, описываемых системой дифференциальных уравнений (4), (5) при различных значениях аэродинамического параметра и коэффициентов демпфирования. Перепишем систему (4), (5) в стандартном виде, удобном для численного интегрирования при условии, когда ${{k}_{1}} = {{k}_{2}} = {{k}_{3}} = k$:

Численное интегрирование системы (15) проводилось при фиксированных значениях параметров, удовлетворяющих условиям (13). Значения коэффициента демпфирования $k$ изменялись в диапазоне от 0.5 до 2.0, а значения аэродинамического параметра ${{h}_{1}}$ – в диапазоне от 1.0 до 50.0. Исследование решений системы (15) проводилось для случаев: ${{\theta }_{A}} = 0.8$, ${{\theta }_{C}} = 0.4$ и ${{\theta }_{A}} = 0.24$, ${{\theta }_{C}} = 0.95$. При выборе второго варианта условий для инерциальных параметров принималось во внимание, что указанные параметры являлись оптимальными для спутника с аэрогироскопической системой стабилизации при изучении динамики этой системы в работе [5]. На фиг. 2–4 приведены примеры переходных процессов, описываемых системой (15) для значений параметров из диапазонов, указанных выше. Начальные значения переменных при расчетах принимались равными 0.001.

Из фигур видно, что при малых величинах коэффициента демпфирования и аэродинамического момента переходный процесс системы происходит медленно, величина $\tau $ при этом превышает 20 (фиг. 2а, $k = 0.5$, ${{h}_{1}} = 1.0$). При увеличении значения коэффициента демпфирования время переходного процесса системы в положение равновесия уменьшается. При $k = 1.0$, ${{h}_{1}} = 25.0$ время переходного процесса меньше 10 (фиг. 3а), а при $k = 2.0$, ${{h}_{1}} = 25.0$ время переходного процесса меньше 6 (фиг. 3б). С ростом величины аэродинамического момента происходит увеличение частоты колебаний спутника по углам $\alpha $, $\beta $ (фиг. 2б, $k = 0.5$, ${{h}_{1}} = 25.0$; фиг. 4б, $k = 1.0$, ${{h}_{1}} = 50.0$). При больших величинах коэффициента демпфирования и аэродинамического момента переходный процесс системы происходит за один оборот спутника по орбите (фиг. 3б, $k = 2.0$, ${{h}_{1}} = 25.0$, $\tau \approx 6$).

Для случаев, близких к осесимметричному (${{\theta }_{A}} = 0.24$, ${{\theta }_{C}} = 0.95$), время переходного процесса по углу крена $\gamma $ превышает времена переходных процессов по углам $\alpha $, $\beta $ и значительно возрастает при увеличении значений аэродинамического момента (фиг. 4а, $k = 1.0$, ${{h}_{1}} = 5.0$; фиг. 4б, $k = 1.0$, ${{h}_{1}} = 50.0$).

5. ЗАКЛЮЧЕНИЕ

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

Для нулевого положения равновесия, когда оси связанной со спутником системы координат совпадают с осями орбитальной системы координат, получены необходимые и достаточные условия асимптотической устойчивости. Проведен анализ характера изменений областей, где выполняются условия асимптотической устойчивости положения равновесия в плоскости параметров $({{\theta }_{A}},\;{{\theta }_{C}})$ для различных значений аэродинамического параметра и коэффициентов демпфирования. Показано, что при росте значения ${{h}_{1}}$ от 0 до 3 происходит увеличение области устойчивости. При ${{h}_{1}} > 3$ размер области устойчивости практически не изменяется.

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

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