Исследование движения центра масс межпланетных космических аппаратов
Страница 5
| rc - r| = Ö((xc-x)2+(yc-y)2+(zc-z)2)
где xc, yc, zc - проекции радиуса-вектора Солнца на оси абсолютной системы координат.
Моделирование движения Солнца проводилось следующим образом: за некоторый промежуток времени t Солнце относительно Земли сместится на угол J = Jн + wct,
где Jн = W + (90 - D) - начальное положение Солнца в эклиптической системе координат.
W = 28,1° - долгота восходящего узла первого витка КА.
D = 30° - угол между восходящим узлом орбиты КА и терминатором.
wc - угловая скорость Солнца относительно Земли.
wc = 2p/T = 2p/365,2422´24´3600 = 1,991´10-7 рад/c = 1,14´10-5 °/c
Таким образом, в эклиптической системе координат проекции составляют:
xce = rccosJ
yce = rcsinJ
zce = 0
rc = 1,496´1011 м (1 астрономическая единица) - расстояние от Земли до Солнца
Плоскость эклиптики наклонена к плоскости экватора на угол e = 23,45°, проекции rc на оси абсолютной системы координат можно найти как
xc = xce = rccosJ
yce = ycecose = rccosJcose
zce = rcsinJsine
Таким образом, проекции возмущающего ускорения на оси абсолютной системы координат:
axc = - mcx/(Ö((xc-x)2+(yc-y)2+(zc-z)2))3
ayc = - mcy/(Ö((xc-x)2+(yc-y)2+(zc-z)2))3
azc = - mcz/(Ö((xc-x)2+(yc-y)2+(zc-z)2))3
С учетом солнечного давления
axc = - (mc-Dmc)x/(Ö((xc-x)2+(yc-y)2+(zc-z)2))3
ayc = - (mc-Dmc)y/(Ö((xc-x)2+(yc-y)2+(zc-z)2))3
azc = - (mc-Dmc)z/(Ö((xc-x)2+(yc-y)2+(zc-z)2))3
5) Возмущающее ускорение, возникающее из-за влияния Луны.
Уравнение движения КА в абсолютной системе координат OXYZ относительно Земли при воздействии Луны:
где mл = 4,902´106 м3/c2- постоянная тяготения Луны.
rл - радиус-вектор от Земли до Луны.
Таким образом, возмущающее ускорение, возникающее из-за влияния Луны:
Так как rл>>r, то в первом слагаемом можно пренебречь r. Следовательно
|rл - r| = Ö((xл-x)2+(yл-y)2+(zл-z)2)
где xл, yл, zл - проекции радиуса-вектора Луны на оси абсолютной системы координат.
Движение Луны учитывается следующим образом: положение Луны в каждый момент времени рассчитывается в соответствии с данными астрономического ежегодника. Все данные заносятся в массив, и далее этот массив считается программой моделирования движения КА. В первом приближении принимается:
- орбита Луны - круговая.
- угол наклона плоскости орбиты Луны к плоскости эклиптики i = 5,15°.
- период обращения линии пересечения плоскостей лунной орбиты и эклиптики (по ходу часовой стрелки, если смотреть с северного полюса) = 18,6 года.
Угол между плоскостями экватора Земли и орбиты Луны можно найти по формуле
cos(hл) = cos(e)cos(i) - sin(e)sin(i)cos(Wл)
где Wл - долгота восходящего узла лунной орбиты, отсчитывается от направления на точку весеннего равноденствия.
e - угол между плоскостями эклиптики и экватора Земли.
Величина hл колеблется с периодом 18,6 лет между минимумом при hл = e - i = 18°18’ и максимумом при hл = e + i = 28°36’ при W = 0.
Долгота восходящего узла лунной орбиты Wл изменяется с течением времени t на величину Wл = t´360/18,6´365,2422´24´3600.
Положение Луны на орбите во время t определяется углом
J л = t´360/27,32´24´3600.
По формулам перехода найдем проекции вектора положения Луны на оси абсолютной системы координат:
xл = rл(cosJлcosWл - coshлsinJлsinWл)
yл = rл(cosJлsinWл + coshлsinJлcosWл)
zл = rлsinhлsinJл
rл = 3,844´108 м - среднее расстояние от Земли до Луны
Таким образом, проекции возмущающего ускорения на оси абсолютной системы координат:
axл = - mлx/(Ö((xл!-x)2+(yл-y)2+(zл-z)2))3
ayл = - mлy/(Ö((xл!-x)2+(yл-y)2+(zл-z)2))3
azл = - mлz/(Ö((xл!-x)2+(yл-y)2+(zл-z)2))3
Уравнения возмущенного движения при действии корректирующего ускорения имеют вид:
или
d2x/dt2 = - (mz/r2)x + axu + axa + axc + axл + axк
d2y/dt2 = - (mz/r2)y + ayu + aya + ayc + ayл + ayк
d2z/dt2 = - (mz/r2)z + azu + aza + azc + azл + azк
2.4.3. РАСЧЕТ ПАРАМЕТРОВ ТЕКУЩЕЙ ОРБИТЫ КА
Полученная система уравнений движения ЦМ КА интегрируется методом Рунге-Кутта 5-го порядка с переменным шагом. Начальные условия x0, y0, z0, Vx0, Vy0, Vz0 - в абсолютной системе координат, соответствуют начальной точке вывода при учете ошибок выведения. После интегрирования мы получаем вектор состояния КА (x, y, z, Vx, Vy, Vz) в любой момент времени.
По вектору состояния можно рассчитать параметры орбиты. соответствующие этому вектору состояния.
а) Фокальный параметр - р.
р = C2/mz, где С - интеграл площадей.
C = r ´ V, |C| = C = Ö(Cx2+Cy2+Cz2)
Cx = yVz - zVy
Cy = zVx - xVz - проекции на оси абсолютной СК
Cz = xVy - yVx
б) Эксцентриситет - е.
e = f/mz, где f - вектор Лапласа
f = V ´ C - mzr/r, |f| = f = Ö(fx2+fy2+fz2)
fx = VyCz - VzCy - mzx/r
fy = VzCx - VxCz - mzy/r - проекции на оси абсолютной СК
fz = VxCy - VyCx - mzz/r
в) Большая полуось орбиты.
a = p/(1 - e2)
г) Наклонение орбиты - i.
Cx = Csin(i)sinW
Cy = - Csin(i)cosW
Cz = Ccos(i)
можно найти наклонение i = arccos(Cz/C)
д) Долгота восходящего узла - W.
Из предыдущей системы можно найти
sinW = Cx/Csin(i)
cosW = - Cy/Csin(i)
Так как наклонение орбиты изменяется несильно в районе i = 97,6°, мы имеем право делить на sin(i).
Если sinW => 0, W = arccos (-Cy/Csin(i))
Если sinW < 0, W = 360 - arccos (-Cy/Csin(i))
е) Аргумент перицентра - w.
fx = f(coswcosW - sinwsinWcos(i))
fy = f(coswsinW + sinwcosWcos(i))
fz = fsinwsin(i)
Отсюда найдем
cosw = fxcosW/f + fysinW/f
sinw = fz/fsin(i)
Если sinw > 0, w = arccos (fxcosW/f + fysinW/f)
Если sinw < 0, w = 360 - arccos (fxcosW/f + fysinW/f)
ж) Период обращения - Т.
T = 2pÖ(a3/mz)
Графики изменения элементов орбиты при действии всех, рассмотренных выше, возмущающих ускорений в течение 2-х периодов (Т = 5765 с) приведены на рис. 1-12.
Графики изменения во времени возмущающих ускорений приведены на рис. 13-18.
2.5. ПРОВЕДЕНИЕ КОРРЕКЦИИ ТРАЕКТОРИИ МКА
Существующие ограничения на точки старта РН и зоны падения отработавших ступеней РН, а также ошибки выведения не позволяют сразу же после пуска реализовать рабочую орбиту. Кроме того, эволюция параметров орбит под действием возмущающих ускорений в процессе полета МКА приводит к отклонению параметров орбиты КА от требуемых значений. Для компенсации воздействия указанных факторов осуществляется коррекция орбиты с помощью корректирующей двигательной установки (КДУ), которая располагается на борту МКА.
В данной работе проведена разработка алгоритма коррекции, моделирование процесса коррекции и расчет топлива, необходимого для проведения коррекции.
Из-за различных причин возникновения отклонений элементов орбиты проводится:
- коррекция приведения - ликвидация ошибок выведения и приведение фактической орбиты к номинальной с заданной точностью.
- коррекция поддержания - ликвидация отклонений параметров орбиты от номинальных, возникающих из-за действия возмущающих ускорений в процессе полета.
Для того, чтобы орбита отвечала заданным требованиям, отклонения параметров задаются следующим образом:
- максимальное отклонение наклонения орбиты Di = 0,1°
- предельное суточное смещение КА по долготе Dl = 0,1°
Следовательно, максимальное отклонение периода орбиты DT = 1,6 сек.
Алгоритм коррекции следующий:
1) Коррекция приведения.
2) Коррекция поддержания.