|Vx, м/с |-206,3 |
|Vy, м/с |-1252,03 |
|Vz, м/с |7477,65 |
|(, ( |28,1 |
Параметры орбиты с учетом ошибок выведения:
|(, ( |28,13 |
|T, c |5795,7 |
|(, ( |28,13 |
|p, км |6973,5 |
|а, км |6973,6 |
|e |0,00314 |
|i, ( |97,637 |
2.3.2. ЦЕЛИ РАБОТЫ
1) Исследование и моделирование движения ЦМ МКА при воздействии на КА
возмущающих ускорений.
2) Разработка алгоритмов проведения коррекции траектории МКА,
моделирования процесса, и расчет потребного топлива для проведения
коррекции траектории.
3) Исследование динамики системы коррекции траектории при стабилизации
углового положения в процессе проведения коррекции траектории МКА.
2.4. МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ЦЕНТРА МАСС МКА
2.4.1.УРАВНЕНИЕ ДВИЖЕНИЯ КА
Рассмотрим невозмущенное движение материальных точек М и m в некоторой
инерциальной системе координат. Движение совершается под действием силы
притяжения Fz. Сила Fz для материальной точки m определяется формулой:
[pic],
где ( - постоянная притяжения,
ro - единичный вектор, направленный от М к m,
[pic],
где [pic]- радиус-вектор, проведенный из т.М до т.m.
r - относительное расстояние от М до m.
На точку М действует сила Fz, равная по величине и направленная в
противоположную сторону.
На основе второго закона Ньютона уравнения движения материальных точек М
и m имеют вид:
[pic](1), [pic] (2)
или
[pic](3), [pic] (4)
где p1 - радиус-вектор, проведенный из начала инерциальной системы
координат в точку m.
p2 - радиус-вектор, проведенный из начала инерциальной системы координат
в точку М.
[pic].
Вычитая из уравнения (3) уравнение (4), получим уравнение движения
материальной точки m относительно притягивающего центра М:
[pic][pic]
Так как m<<М, следовательно, можно пренебречь ускорением, которое КА с
массой m сообщает притягивающему центру М. Тогда можно совместить начало
инерциальной системы координат с притягивающим центром М. Следовательно,
[pic].
Таким образом, уравнение невозмущенного движения КА относительно
притягивающего центра М в инерциальной системе координат, центр которой
находится в М, имеет вид
[pic],[pic]
где ( = fM - гравитационная постоянная Земли.
Рассмотрим возмущенное движение КА в геоцентрической экваториальной
(абсолютной) системе координат OXYZ:
- начало О - в центре масс Земли.
- ось X направлена в точку весеннего равноденствия (.
- ось Z совпадает с осью вращения Земли и направлена на Северный полюс
Земли.
- ось Y дополняет систему до правой.
Движение КА в абсолютной системе координат OXYZ происходит под действием
центральной силы притяжения Земли Fz, а также под действием возмущающих сил
Fв. Уравнение движения имеет вид
[pic] или [pic]
где m = 597 кг - масса КА.
В проекциях на оси абсолютной системы координат OXYZ получим
[pic] или [pic]
[pic] или [pic]
[pic] или [pic]
где axв, ayв, azв - возмущающиеся ускорения.
Основные возмущающиеся ускорения вызываются следующими причинами:
- нецентральностью поля притяжения Земли.
- сопротивлением атмосферы Земли.
- влиянием Солнца.
- влиянием Луны.
- давлением солнечного света.
2.4.2. ВОЗМУЩАЮЩИЕ УСКОРЕНИЯ, ДЕЙСТВУЮЩИЕ НА МКА
1) Возмущающееся ускорение, вызванное нецентральностью гравитационного
поля Земли.
Рассмотрим потенциал поля притяжения Земли. При точном расчете параметров
орбиты спутников, в качестве хорошего приближения к действительной
поверхности Земли принимают геоид. Геоид - это гипотетическая уровенная
поверхность, совпадающая с поверхностью спокойного океана и продолженная
под материком.
Иногда в баллистике под геоидом понимают не поверхность, а тело, которое
ограничено поверхностью мирового океана при некотором среднем уровне воды,
свободной от возмущений. Во всех точках геоида потенциал притяжения имеет
одно и то же значение.
Потенциал притяжения Земли можно представить в виде разложения по
сферическим функциям.
[pic]
где (z = fMz - гравитационная постоянная Земли.
r0 - средний экваториальный радиус Земли.
сnm, dnm - коэффициенты, определяемые из гравиметрических данных, а также
по наблюдениям за движением ИСЗ.
L - долгота притягивающей точки.
( - широта притягивающей точки.
Pnm(sin() - присоединенные функции Лежандра степени m и порядка n (при m
( 0).
Pnm(sin() - многочлен Лежандра порядка n (при m = 0).
Составляющие типа ((z/r)(r0/r)ncn0Pn0(sin() - называют зональными
гармониками n-порядка. Т.к. полином Лежандра n-го порядка имеет n
действительных корней, функция Pn0(sin() будет менять знак на n широтах,
сфера делится на n+1 широтную зону, где эти составляющие имеют попеременно
«+» или «-» значения. Поэтому их называют зональными гармониками.
Составляющие типа
((z/r)(r0/r)ncnmcos(mL)Pnm(sin() и ((z/r)(r0/r)ndnmsin(mL)Pnm(sin()
- называют тессеральными гармониками n-порядка и степени m. Они
обращаются в 0 на 2m меридианах, где cos(mL) = 0 и sin(mL) = 0 и на n-m
параллелях, где Pnm(sin() = 0 или dmPnm(sin()/d(sin()m = 0, сфера делится
на n+m+1 трапецию, где эти составляющие сохраняют знак.
Составляющие типа и
((z/r)(r0/r)ncnncos(nL)Pnn(sin() и ((z/r)(r0/r)ndnnsin(nL)Pnn(sin()
- называют секториальными гармониками n-порядка и степени m. Эти
составляющие меняю знак только на меридианах, cos(nL) = 0 и sin(nL) = 0, на
сфере выделяют 2n меридиональных секторов, где эти составляющие сохраняют
знак.
Многочлен Лежандра степени n находится по следующей формуле:[pic]
Pn0(z) = 1/(2nn!)((dn(z2 - 1)n/dzn)
Присоединенная функция Лежандра порядка n и степени m находится по
следующей формуле:
Pnm(z) = (1-z2)m/2(dmPn0(z)/dzm
Возмущающая часть гравитационного потенциала Земли равна
Uв = U’ + (U’ = (U - (z/r) + (U’
где (U’ - потенциал аномалий силы тяготения Земли.
U’ - часть потенциала Земли, которая учитывает несферичность Земли.
Следовательно,
[pic]
Первая зональная гармоника в разложении потенциала учитывает полярное
сжатие Земли.
Зональные гармоники нечетного порядка и тессеральные гармоники, где n-m
нечетное число - учитывают ассиметрию Земли относительно плоскости
экватора.
Секториальные и тессеральные гармоники - учитывают ассиметрию Земли
относительно оси вращения.
Первая зональная гармоника имеет порядок 10-3, а все остальные - порядок
10-6 и выше. Поэтому будем учитывать в разложении потенциала притяжения
только зональную гармонику (n=2, m=0) и секторальную гармонику (n=2, m=2).
Также не будем учитывать потенциал аномалий силы тяготения Земли (U’.
Таким образом,
Uв = ((z/r)(r0/r)2[c20P20(sin() + (c22cos(2L) + d22sin(2L))P22(sin()],
где c20 = - 0,00109808,
c22 = 0,00000574,
d22 = - 0,00000158.
P20(x) = 1/222!(d2(x2 - 1)2/dx2.
Следовательно P20(x) = (3x2 - 1)/2.
Так как sin( = z/r, следовательно P20(sin() = (3(z/r)2 - 1)/2.
P22(x) = (1 - x2)2/2(d2P20(x)/dx2 = 1/2((1 - x2)(d2(3x2 - 1)/dx2
Следовательно P22(x) = 3(1 - x2).
Так как sin( = z/r, следовательно P22(sin() = 3(1 - (z/r)2).
Значит
[pic]
Чтобы найти возмущающее ускорение от нецентральности поля тяготения Земли
в проекциях на оси абсолютной системы координат OXYZ, надо взять
производные от возмущающего потенциала Uв по координатам X, Y, Z, причем r
= ((x2 + y2 + z2).
Следовательно,
[pic][pic][pic][pic][pic][pic]
2) Возмущающее ускорение, вызванное сопротивлением атмосферы.
При движении в атмосфере на КА действует сила аэродинамического ускорения
Rx, направленная против вектора скорости КА относительно атмосферы:
[pic]
где Cx = 2 - коэффициент аэродинамического сопротивления.
Sм = 2,5 м2 - площадь миделевого сечения - проекция КА на плоскость,
перпендикулярную направлению скорости полета.
V - скорость КА.
( - плотность атмосферы в рассматриваемой точке орбиты.
Так как исследуемая орбита - круговая с высотой Н = 574 км, будем
считать, что плотность атмосферы одинакова во всех точках орбиты и равна
плотности атмосферы на высоте 574 км. Из таблицы стандартной атмосферы
находим плотность наиболее близкую к высоте Н = 574 км. Для высоты Н = 580
км ( = 5,098(10-13 кг/м3.
Сила аэродинамического ускорения создает возмущающее касательное
ускорение aa:
[pic]
Найдем проекции аэродинамического ускорения на оси абсолютной системы
координат axa, aya, aza:
aa направлено против скорости КА, следовательно единичный вектор
направления имеет вид
ea = [Vx/|V|, Vy|V|, Vz/|V|], |V| = ((Vx2+Vy2 +Vz2)
Таким образом,
[pic]
Значит
[pic], [pic], [pic]
3) Возмущающее ускорение, вызванное давлением солнечного света.
Давление солнечного света учитывается как добавок к постоянной тяготения
Солнца - ((c. Эта величина вычисляется следующим образом:
((c = pSмA2/m
где p = 4,64(10-6 Н/м2 - давление солнечного света на расстоянии в одну
астрономическую единицу А.
A = 1,496(1011 м - 1 астрономическая единица.
m - масса КА.
Sм = 8 м2 - площадь миделевого сечения - проекция КА на плоскость,
перпендикулярную направления солнечных лучей.
Таким образом,
((c = 1,39154(1015 м3/c2.
4) Возмущающее ускорение, возникающее из-за влияния Солнца.
Уравнение движения КА в абсолютной системе координат OXYZ относительно
Земли при воздействии Солнца:
[pic]
где (z - постоянная тяготения Земли.
(c - постоянная тяготения Солнца.
r - радиус-вектор от Земли до КА.
rc - радиус-вектор от Земли до Солнца.
Таким образом, возмущающее ускорение, возникающее из-за влияния Солнца:
[pic].
Здесь первое слагаемое есть ускорение, которое получил бы КА, если он был
непритягивающим, а Земля отсутствовала.
Второе слагаемое есть ускорение, которое сообщает Солнце Земле, как
непритягивающему телу.
Следовательно, возмущающее ускорение, которое получает КА при движении
относительно Земли - это разность двух слагаемых.
[pic]
Так как rc>>r, то в первом слагаемом можно пренебречь r. Следовательно
[pic]
| rc - r| = (((xc-x)2+(yc-y)2+(zc-z)2)
где xc, yc, zc - проекции радиуса-вектора Солнца на оси абсолютной
системы координат.
Моделирование движения Солнца проводилось следующим образом: за некоторый
промежуток времени t Солнце относительно Земли сместится на угол ( = (н +
(ct,
где (н = ( + (90 - () - начальное положение Солнца в эклиптической
системе координат.
( = 28,1( - долгота восходящего узла первого витка КА.
( = 30( - угол между восходящим узлом орбиты КА и терминатором.
(c - угловая скорость Солнца относительно Земли.
(c = 2(/T = 2(/365,2422(24(3600 = 1,991(10-7 рад/c = 1,14(10-5 (/c
Таким образом, в эклиптической системе координат проекции составляют:
xce = rccos(
yce = rcsin(
zce = 0
rc = 1,496(1011 м (1 астрономическая единица) - расстояние от Земли до
Солнца
Плоскость эклиптики наклонена к плоскости экватора на угол ( = 23,45(,
проекции rc на оси абсолютной системы координат можно найти как
xc = xce = rccos(
yce = ycecos( = rccos(cos(
zce = rcsin(sin(
Таким образом, проекции возмущающего ускорения на оси абсолютной системы
координат:
axc = - (cx/((((xc-x)2+(yc-y)2+(zc-z)2))3
ayc = - (cy/((((xc-x)2+(yc-y)2+(zc-z)2))3
azc = - (cz/((((xc-x)2+(yc-y)2+(zc-z)2))3
С учетом солнечного давления
axc = - ((c-((c)x/((((xc-x)2+(yc-y)2+(zc-z)2))3
ayc = - ((c-((c)y/((((xc-x)2+(yc-y)2+(zc-z)2))3
azc = - ((c-((c)z/((((xc-x)2+(yc-y)2+(zc-z)2))3
5) Возмущающее ускорение, возникающее из-за влияния Луны.
Уравнение движения КА в абсолютной системе координат OXYZ относительно
Земли при воздействии Луны:
[pic]
где (л = 4,902(106 м3/c2- постоянная тяготения Луны.
rл - радиус-вектор от Земли до Луны.
Таким образом, возмущающее ускорение, возникающее из-за влияния Луны:
[pic]
Так как rл>>r, то в первом слагаемом можно пренебречь r. Следовательно
[pic]
|rл - r| = (((xл-x)2+(yл-y)2+(zл-z)2)
где xл, yл, zл - проекции радиуса-вектора Луны на оси абсолютной системы
координат.
Движение Луны учитывается следующим образом: положение Луны в каждый
момент времени рассчитывается в соответствии с данными астрономического
ежегодника. Все данные заносятся в массив, и далее этот массив считается
программой моделирования движения КА. В первом приближении принимается:
- орбита Луны - круговая.
- угол наклона плоскости орбиты Луны к плоскости эклиптики i = 5,15(.
- период обращения линии пересечения плоскостей лунной орбиты и эклиптики
(по ходу часовой стрелки, если смотреть с северного полюса) = 18,6 года.
Угол между плоскостями экватора Земли и орбиты Луны можно найти по
формуле
cos((л) = cos(()cos(i) - sin(()sin(i)cos((л)
где (л - долгота восходящего узла лунной орбиты, отсчитывается от
направления на точку весеннего равноденствия.
( - угол между плоскостями эклиптики и экватора Земли.
Величина (л колеблется с периодом 18,6 лет между минимумом при (л = ( - i
= 18(18’ и максимумом при (л = ( + i = 28(36’ при ( = 0.
Долгота восходящего узла лунной орбиты (л изменяется с течением времени t
на величину (л = t(360/18,6(365,2422(24(3600.
Положение Луны на орбите во время t определяется углом
( л = t(360/27,32(24(3600.
По формулам перехода найдем проекции вектора положения Луны на оси
абсолютной системы координат:
xл = rл(cos(лcos(л - cos(лsin(лsin(л)
yл = rл(cos(лsin(л + cos(лsin(лcos(л)
zл = rлsin(лsin(л
rл = 3,844(108 м - среднее расстояние от Земли до Луны
Таким образом, проекции возмущающего ускорения на оси абсолютной системы
координат:
axл = - (лx/((((xл!-x)2+(yл-y)2+(zл-z)2))3
ayл = - (лy/((((xл!-x)2+(yл-y)2+(zл-z)2))3
azл = - (лz/((((xл!-x)2+(yл-y)2+(zл-z)2))3
Уравнения возмущенного движения при действии корректирующего ускорения
имеют вид:
[pic]
или
d2x/dt2 = - ((z/r2)x + axu + axa + axc + axл + axк
d2y/dt2 = - ((z/r2)y + ayu + aya + ayc + ayл + ayк
d2z/dt2 = - ((z/r2)z + azu + aza + azc + azл + azк
2.4.3. РАСЧЕТ ПАРАМЕТРОВ ТЕКУЩЕЙ ОРБИТЫ КА
Полученная система уравнений движения ЦМ КА интегрируется методом Рунге-
Кутта 5-го порядка с переменным шагом. Начальные условия x0, y0, z0, Vx0,
Vy0, Vz0 - в абсолютной системе координат, соответствуют начальной точке
вывода при учете ошибок выведения. После интегрирования мы получаем вектор
состояния КА (x, y, z, Vx, Vy, Vz) в любой момент времени.
По вектору состояния можно рассчитать параметры орбиты. соответствующие
этому вектору состояния.
а) Фокальный параметр - р.
р = C2/(z, где С - интеграл площадей.
C = r ( V, |C| = C = ((Cx2+Cy2+Cz2)
Cx = yVz - zVy
Cy = zVx - xVz - проекции на оси абсолютной СК
Cz = xVy - yVx
б) Эксцентриситет - е.
e = f/(z, где f - вектор Лапласа
f = V ( C - (zr/r, |f| = f = ((fx2+fy2+fz2)
fx = VyCz - VzCy - (zx/r
fy = VzCx - VxCz - (zy/r - проекции на оси абсолютной СК
fz = VxCy - VyCx - (zz/r
в) Большая полуось орбиты.
a = p/(1 - e2)
г) Наклонение орбиты - i.
Cx = Csin(i)sin(
Cy = - Csin(i)cos(
Cz = Ccos(i)
можно найти наклонение i = arccos(Cz/C)
д) Долгота восходящего узла - (.
Из предыдущей системы можно найти
sin( = Cx/Csin(i)
cos( = - Cy/Csin(i)
Так как наклонение орбиты изменяется несильно в районе i = 97,6(, мы
имеем право делить на sin(i).
Если sin( => 0, ( = arccos (-Cy/Csin(i))
Если sin( < 0, ( = 360 - arccos (-Cy/Csin(i))
е) Аргумент перицентра - (.
fx = f(cos(cos( - sin(sin(cos(i))
fy = f(cos(sin( + sin(cos(cos(i))
fz = fsin(sin(i)
Отсюда найдем
cos( = fxcos(/f + fysin(/f
sin( = fz/fsin(i)
Если sin( > 0, ( = arccos (fxcos(/f + fysin(/f)
Если sin( < 0, ( = 360 - arccos (fxcos(/f + fysin(/f)
ж) Период обращения - Т.
T = 2(((a3/(z)
Графики изменения элементов орбиты при действии всех, рассмотренных выше,
возмущающих ускорений в течение 2-х периодов (Т = 5765 с) приведены на рис.
1-12.
Графики изменения во времени возмущающих ускорений приведены на рис. 13-
18.
2.5. ПРОВЕДЕНИЕ КОРРЕКЦИИ ТРАЕКТОРИИ МКА
Существующие ограничения на точки старта РН и зоны падения отработавших
ступеней РН, а также ошибки выведения не позволяют сразу же после пуска
реализовать рабочую орбиту. Кроме того, эволюция параметров орбит под
действием возмущающих ускорений в процессе полета МКА приводит к отклонению
параметров орбиты КА от требуемых значений. Для компенсации воздействия
указанных факторов осуществляется коррекция орбиты с помощью корректирующей
двигательной установки (КДУ), которая располагается на борту МКА.
В данной работе проведена разработка алгоритма коррекции, моделирование
