: Исследование движения центра масс межпланетных космических аппаратов

имеют вид:

(1), (2)

или

(3), (4)

где p1 - радиус-вектор, проведенный из начала инерциальной

сис­темы координат в точку m.

p2 - радиус-вектор, проведенный из начала инерциальной

сис­темы координат в точку М.

.

Вычитая из уравнения (3) уравнение (4), получим уравнение дви­жения

мате­риальной точки m относительно притягивающего цен­тра М:

Так как m<<М, следовательно, можно пренебречь ускорением, которое КА с

массой m сообщает притягивающему центру М. То­гда можно совместить начало

инерциальной системы координат с при­тягивающим центром М. Следовательно,

.

Таким образом, уравнение невозмущенного движения КА отно­сительно

притягивающего центра М в инерциальной системе коор­динат, центр которой

находится в М, имеет вид

,

где m = fM - гравитационная постоянная Земли.

Рассмотрим возмущенное движение КА в геоцентрической эква­ториальной

(абсолютной) системе координат OXYZ:

- начало О - в центре масс Земли.

- ось X направлена в точку весеннего равноденствия g.

- ось Z совпадает с осью вращения Земли и направлена на Север­ный полюс Земли.

- ось Y дополняет систему до правой.

Движение КА в абсолютной системе координат OXYZ происхо­дит под действием

центральной силы притяжения Земли Fz, а также под действием

возмущающих сил Fв. Уравнение движения имеет вид

или

где m = 597 кг - масса КА.

В проекциях на оси абсолютной системы координат OXYZ полу­чим

или

или

или

где axв, ayв, azв - возмущающиеся ускорения.

Основные возмущающиеся ускорения вызываются следующими причинами:

- нецентральностью поля притяжения Земли.

- сопротивлением атмосферы Земли.

- влиянием Солнца.

- влиянием Луны.

- давлением солнечного света.

2.4.2. ВОЗМУЩАЮЩИЕ УСКОРЕНИЯ, ДЕЙСТВУЮЩИЕ НА МКА

1) Возмущающееся ускорение, вызванное нецентральностью гра­витационного поля

Земли.

Рассмотрим потенциал поля притяжения Земли. При точном рас­чете параметров

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

поверхности Земли принимают геоид. Геоид - это гипотетическая уровенная

поверхность, совпадающая с поверх­ностью спокойного океана и продолженная под

материком.

Иногда в баллистике под геоидом понимают не поверхность, а тело, которое

ограничено поверхностью мирового океана при не­ко­тором среднем уровне воды,

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

одно и то же значение.

Потенциал притяжения Земли можно представить в виде разло­же­ния по

сферическим функциям.

где mz = fMz - гравитационная постоянная Земли.

r0 - средний экваториальный радиус Земли.

сnm, dnm - коэффициенты, определяемые из

гравиметрических дан­ных, а также по наблюдениям за движением ИСЗ.

L - долгота притягивающей точки.

j - широта притягивающей точки.

Pnm(sinj) - присоединенные функции Лежандра степени m и

по­рядка n (при m ¹ 0).

Pnm(sinj) - многочлен Лежандра порядка n (при m = 0).

Составляющие типа (mz/r)(r0/r)ncn0

Pn0(sinj) - называют зональ­ными гармониками n-по­рядка. Т.к.

полином Лежандра n-го по­рядка имеет n действительных корней, функция P

n0(sinj) будет ме­нять знак на n широтах, сфера делится на n+1

широтную зону, где эти составляю­щие имеют попеременно «+» или «-» значения.

По­этому их называют зональными гармониками.

Составляющие типа

(mz/r)(r0/r)ncnmcos(mL)Pnm(sinj) и (mz/r)(r0/r)ndnmsin(mL)Pnm(sinj)

- называют тессеральными гармониками n-порядка и степени m. Они обращаются в 0

на 2m меридианах, где cos(mL) = 0 и sin(mL) = 0 и на n-m

параллелях, где Pnm(sinj) = 0 или dmPnm

(sinj)/d(sinj)m = 0, сфера делится на n+m+1 трапецию, где эти

составляющие сохра­няют знак.

Составляющие типа и

(mz/r)(r0/r)ncnncos(nL)Pnn(sinj) и (mz/r)(r0/r)ndnnsin(nL)Pnn(sinj)

- называют секториальными гармониками n-порядка и степени m. Эти составляющие

меняю знак только на меридианах, cos(nL) = 0 и sin(nL) = 0, на

сфере выделяют 2n меридиональных секторов, где эти составляющие со­храняют

знак.

Многочлен Лежандра степени n находится по следующей фор­муле:

Pn0(z) = 1/(2nn!)´(dn(z2 - 1)n/dzn)

Присоединенная функция Лежандра порядка n и степени m нахо­дится по следующей

формуле:

Pnm(z) = (1-z2)m/2´dmPn0(z)/dzm

Возмущающая часть гравитационного потенциала Земли равна

Uв = U’ + DU’ = (U - mz/r) + DU’

где DU’ - потенциал аномалий силы тяготения Земли.

U’ - часть потенциала Земли, которая учитывает несферичность Земли.

Следовательно,

Первая зональная гармоника в разложении потенциала учиты­вает полярное сжатие

Земли.

Зональные гармоники нечетного порядка и тессеральные гармо­ники, где n-m

нечетное число - учитывают ассиметрию Земли отно­сительно плоскости

экватора.

Секториальные и тессеральные гармоники - учитывают ассимет­рию Земли

относительно оси вращения.

Первая зональная гармоника имеет порядок 10-3, а все остальные -

порядок 10-6 и выше. Поэтому будем учитывать в разложении

по­тен­циала притяжения только зональную гармонику (n=2, m=0) и секторальную

гармонику (n=2, m=2). Также не будем учитывать потенциал аномалий силы

тяго­тения Земли DU’.

Таким образом,

Uв = (mz/r)(r0/r)2[c20P20(sinj) + (c22cos(2L) + d22sin(2L))P22(sinj)],

где c20 = - 0,00109808,

c22 = 0,00000574,

d22 = - 0,00000158.

P20(x) = 1/222!´d2(x2 - 1)2/dx2.

Следовательно P20(x) = (3x2 - 1)/2.

Так как sinj = z/r, следовательно P20(sinj) = (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).

Так как sinj = z/r, следовательно P22(sinj) = 3(1 - (z/r)2).

Значит

Чтобы найти возмущающее ускорение от нецентральности поля тяготения Земли в

проекциях на оси абсолютной системы коорди­нат OXYZ, надо взять производные от

возмущающего потенциала Uв по координатам X, Y, Z, причем

r = Ö(x2 + y2 + z2).

Следовательно,

2) Возмущающее ускорение, вызванное сопротивлением атмо­сферы.

При движении в атмосфере на КА действует сила аэродинамиче­ского ускорения R

x, направленная против вектора скорости КА от­но­сительно атмосферы:

где Cx = 2 - коэффициент аэродинамического сопротивления.

Sм = 2,5 м2 - площадь миделевого сечения - проекция

КА на плос­кость, пер­пендикулярную направлению скорости полета.

V - скорость КА.

r - плотность атмосферы в рассматриваемой точке орбиты.

Так как исследуемая орбита - круговая с высотой Н = 574 км, бу­дем считать, что

плотность атмосферы одинакова во всех точках ор­биты и равна плотности

атмосферы на высоте 574 км. Из таб­лицы стандартной атмосферы находим плотность

наиболее близ­кую к вы­соте Н = 574 км. Для высоты Н = 580 км r =

5,098´10-13 кг/м3.

Сила аэродинамического ускорения создает возмущающее каса­тельное ускорение aa:

Найдем проекции аэродинамического ускорения на оси абсолют­ной системы координат

aa направлено против скорости КА, следовательно единичный

век­тор направления имеет вид

ea = [Vx/|V|, Vy|V|, Vz/|V|], |V| = Ö(Vx2+Vy2 +Vz2)

Таким образом,

Значит

, ,

3) Возмущающее ускорение, вызванное давлением солнечного света.

Давление солнечного света учитывается как добавок к постоян­ной тяготения Солнца

- Dmc. Эта величина вычисляется следующим об­разом:

Dmc = pSмA2/m

где p = 4,64´10-6 Н/м2 - давление солнечного

света на расстоянии в одну астрономи­ческую единицу А.

A = 1,496´1011 м - 1 астрономическая единица.

m - масса КА.

Sм = 8 м2 - площадь миделевого сечения - проекция

КА на плос­кость, пер­пендикулярную направления солнечных лучей.

Таким образом,

Dmc = 1,39154´1015 м3/c2.

4) Возмущающее ускорение, возникающее из-за влияния Солнца.

Уравнение движения КА в абсолютной системе координат OXYZ относительно Земли

при воздействии Солнца:

где mz - постоянная тяготения Земли.

mc - постоянная тяготения Солнца.

r - радиус-вектор от Земли до КА.

rc - радиус-вектор от Земли до Солнца.

Таким образом, возмущающее ускорение, возникающее из-за влияния Солнца:

.

Здесь первое слагаемое есть ускорение, которое полу­чил бы КА, если он был

непритягиваю­щим, а Земля отсутствовала.

Второе слагаемое есть ускорение, которое сообщает Солнце Земле, как

непритягиваю­щему телу.

Следовательно, возмущающее ускорение, которое получает КА при движении

относительно Земли - это разность двух слагаемых.

Так как rc>>r, то в первом слагаемом можно пренебречь r. Следо­ва­тельно

| 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

rc = 1,496´1011 м (1 астрономическая единица) - расстояние от Земли до Солнца

Плоскость эклиптики наклонена к плоскости экватора на угол e = 23,45°,

проекции rc на оси абсолютной системы координат можно найти

как

xc = xce = rccosJ

zce = rcsinJsine

Таким образом, проекции возмущающего ускорения на оси абсо­лютной системы

координат:

axc = - mcx/(Ö((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

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л)

rл = 3,844´108 м - среднее расстояние от Земли до Луны

Таким образом, проекции возмущающего ускорения на оси абсо­лютной системы

координат:

axл = - mлx/(Ö((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к

d2z/dt2 = - (mz/r2)z + azu + aza + azc + azл + azк

2.4.3. РАСЧЕТ ПАРАМЕТРОВ ТЕКУЩЕЙ ОРБИТЫ КА

Полученная система уравнений движения ЦМ КА интегрируется методом Рунге-Кутта

5-го порядка с переменным шагом. Началь­ные условия x0, y0

, z0, Vx0, Vy0, Vz0 - в абсолютной

системе коорди­нат, соответствуют началь­ной точке вывода при учете ошибок

вы­ведения. После интегриро­вания мы получаем вектор состояния КА (x, y, z, V

x, 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 - вектор Лапласа

Страницы: 1, 2, 3, 4, 5, 6