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

---

Рвссмотрим движение космического аппарата (КА) в центральном гравитационном поле под действием внешних возмущающих сил \( {\bf{S}} \), \( {\bf{T}} \), \( {\bf{W}} \) (Рисунок 1). Положение КА может быть задано в инерциальной декартовой \(OXYZ\) и оскулирующей системах координат. Точка \(O\) совпадает с центром Земли, ось \(OX\) направлена в точку весеннего равноденствия и лежит в экваториальной плоскости; ось \(OZ\) направлена вдоль оси вращения Земли, \(OY\) дополняет систему до правой. Оскулирующая система координат полностью характеризует орбиту КА. Элементами орбиты называют следующие параметры [1 стр 79]:

• наклонение орбиты \( i \) (inclination),
• долгота восходящего узла орбиты \( \Omega \) (longitude of the ascending node),
• аргумент перицентра \( \omega \) (argument of the pericenter),
• большая полуось \( a \) (semi-major axis),
• эксцентриситет \( e \) (eccentricity),
• время прохождения КА через перицентр \( \tau \) (time of periapses passage).

Величины \( \Omega \) и \( i \) характеризуют положение плоскости орбиты в пространстве, а величины \( a \) и \( e \) её размер и форму. Величина \( \omega \) представляет угловое расстояние от восходящего узла \( \Omega \) до перицентра орбиты П, то есть характеризует положение эллипса в пространстве. Для эллиптической орбиты вместо времени прохождения КА через перицентр может быть использована величина \( M_0 \), называемая средней аномалией эпохи (initial mean anomaly) [2 раздел 11.2.3].

\begin{equation} {M_0} = - n\tau, \label{eq:eq1} \end{equation} где \(n = \sqrt {\mu {a^{ - 3}}} \) - среднее движение (mean angular motion), \( \mu \) - гравитационный параметр. Уравнение Кеплера устанавливает связь между средней аномалией \( M=n(t-\tau) \) и эксцентрической аномалией \( E \) (eccentric anomaly) [3 раздел 122] \begin{equation} M = E - e\sin E, \label{eq:eq2} \end{equation} которая в свою очередь связана с углом истинной аномалии \( f \) (true anomaly) \begin{equation} \tan \frac{f}{2} = \sqrt {\frac{{1 + e}}{{1 - e}}} \tan \frac{E}{2} \label{eq:eq3} \end{equation} Показанный на рисунке 1 угол \(\theta = \omega + f\) называется аргументом широты (true latitude angle).
Рисунок 1 - Оскулирующая система координат [1 рис. 2.32]

Уравнения движения в абсолютной прямоугольной системе координат

Дифференциальные уравнения движения КА в абсолютной прямоугольной системе координат имеют вид [1 стр 107, формула (4.1)]: \begin{equation} \begin{array}{l} \ddot X = - \frac{{\mu X}}{{{r^3}}} + \frac{{{F_X}}}{m},\\ \ddot Y = - \frac{{\mu Y}}{{{r^3}}} + \frac{{{F_Y}}}{m},\\ \ddot Z = - \frac{{\mu Z}}{{{r^3}}} + \frac{{{F_Z}}}{m}, \end{array} \label{eq:eq4} \end{equation} где X, Y, Z - координаты КА в OXYZ, m - масса КА, F_X, F_Y, F_Z - проекции возмущающей силы на оси инерциальной системы координат OXYZ, \(r = \sqrt {{X^2} + {Y^2} + {Z^2}} \). Если возмущения заданы в орбитальной системе координат \(C \tau b n\) (Рисунок 1), то для получения F_X, F_Y, F_Z необходимо использовать матрицу перехода \({A_{Eo}}\) [1 стр 83] \begin{equation} \left[ {\begin{array}{*{20}{c}} {{F_X}}\\ {{F_Y}}\\ {{F_Z}} \end{array}} \right] = {A_{Eo}}\left[ {\begin{array}{*{20}{c}} T\\ W\\ S \end{array}} \right] \label{eq:eq5} \end{equation} \begin{equation} {A_{Eo}} = \left[ {\begin{array}{*{20}{c}} {\frac{{{C_2}Z - {C_3}Y}}{{Cr}}}&{\frac{{{C_1}}}{C}}&{\frac{X}{r}}\\ {\frac{{{C_3}X - {C_1}Z}}{{Cr}}}&{\frac{{{C_2}}}{C}}&{\frac{Y}{r}}\\ {\frac{{{C_1}Y - {C_2}X}}{{Cr}}}&{\frac{{{C_3}}}{C}}&{\frac{Z}{r}} \end{array}} \right] \label{eq:eq6} \end{equation}

Здесь компоненты вектора момента количества движения КА относительно центра Земли определяются как \begin{equation} \begin{array}{l} {C_1} = Y{V_Z} - Z{V_Y},\\ {C_2} = Z{V_X} - X{V_Z},\\ {C_3} = X{V_Y} - Y{V_X},\\ C = \sqrt {C_1^2 + C_2^2 + C_3^2} . \end{array} \label{eq:eq7} \end{equation}

Уравнения движения в оскулирующих элементах

Используя элементы орбиты в качестве независимых переменных запишем систему диффеенциальных уравнений, описывающую возмущенное движение КА \begin{equation} \begin{array}{l} \frac{{da}}{{dt}} = \frac{{2a^2}}{{mh}}\left( {eS\sin f\, + \frac{{pT}}{r}} \right),\\ \frac{{de}}{{dt}} = \frac{1}{{mh}}\left( {pS\sin f + ((p + r)\cos f + re)T} \right),\\ \frac{{di}}{{dt}} = \frac{{rW\cos \theta }}{{mh}},\\ \frac{{d\Omega }}{{dt}} = \frac{{rW\sin \theta }}{{mh\sin i}},\\ \frac{{d\omega }}{{dt}} = \frac{{ - pS\cos f + (p + r)T\sin f}}{{mhe}} - \frac{{rW\sin \theta \cos i}}{{mh\sin i}},\\ \frac{{dM}}{{dt}} = n + \frac{{b\left( {(p\cos f - 2re)S - (p + r)T\sin f} \right)}}{{mahe}}. \end{array} \label{eq:eq8} \end{equation}

Здесь \( b = a\sqrt {1 - e_{}^2} \) - малая полуось (semi-minor axis), \( h = \sqrt {\mu p} \). Вместо последнего уравнения системы (8) удобнее использовать \begin{equation} \frac{{df}}{{dt}} = \frac{h}{{r^2}} + \frac{{pS\cos f - (p + r)T\sin f}}{{mhe}} \label{eq:eq9} \end{equation}

Моделирование орбитального движения КА в Matlab

Для проверки системы (8) разработана программа в Matlab, сравнивающая результаты интегирования системы (4) (файл правых частей F_xyz) с системой (8) (файл правых частей F_oe) и с системой (8), в которой последнее уравнение заменено на (9) (файл F_oe2). Результаты совпали. Файл с моделями можно скачать здесь.

Список литературы

1. Нариманов Г.С. Основы теории полета космических аппаратов// М.:Машиностроение, 1972. - 697 с.
2. Schaub H., Junkins J. L. Analytical mechanics of space systems. – 2002, 588 p.
3. Маркеев А. П. Теоретическая механика: учебник для университетов // Ижевск: НИЦ "Регулярная и хаотическая динамика", 2001. - 592 с.