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

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

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

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

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

[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. ПРОВЕДЕНИЕ КОРРЕКЦИИ ТРАЕКТОРИИ МКА

Существующие ограничения на точки старта РН и зоны падения отработавших

ступеней РН, а также ошибки выведения не позволяют сразу же после пуска

реализовать рабочую орбиту. Кроме того, эволюция параметров орбит под

действием возмущающих ускорений в процессе полета МКА приводит к отклонению

параметров орбиты КА от требуемых значений. Для компенсации воздействия

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

двигательной установки (КДУ), которая располагается на борту МКА.

В данной работе проведена разработка алгоритма коррекции, моделирование

процесса коррекции и расчет топлива, необходимого для проведения коррекции.

Из-за различных причин возникновения отклонений элементов орбиты

проводится:

- коррекция приведения - ликвидация ошибок выведения и приведение

фактической орбиты к номинальной с заданной точностью.

- коррекция поддержания - ликвидация отклонений параметров орбиты от

номинальных, возникающих из-за действия возмущающих ускорений в процессе

полета.

Для того, чтобы орбита отвечала заданным требованиям, отклонения

параметров задаются следующим образом:

- максимальное отклонение наклонения орбиты (i = 0,1(

- предельное суточное смещение КА по долготе (( = 0,1(

Следовательно, максимальное отклонение периода орбиты (T = 1,6 сек.

Алгоритм коррекции следующий:

1) Коррекция приведения.

2) Коррекция поддержания.

2.5.1. КОРРЕКЦИЯ ПРИВЕДЕНИЯ

После окончания процесса выведения МКА, проводятся внешне-траекторные

измерения (ВТИ). Эти измерения обеспечивают, по баллистическим расчетам,

знание вектора состояния с требуемой точностью через 2 суток. После этого

начинается коррекция приведения.

Предложена следующая схема проведения коррекции:

а) Коррекция периода.

б) Коррекция наклонения.

Корректирующий импульс прикладывается в апсидальных точках, либо на линии

узлов в течение 20 сек и происходит исправление одного параметра орбиты.

Таким образом используется однопараметрическая, непрерывная коррекция.

а) Коррекция периода.

Осуществляется в два этапа:

- коррекция перицентра

- коррекция апоцентра

Сначала осуществляется коррекция перицентра - приведение текущего

расстояния до перицентра r( к номинальному радиусу rн = 6952137 м. После

измерения вектора состояния рассчитываются параметры орбиты. Далее

определяется нужный корректирующий импульс (Vк. Направление импульса

(тормозящий или разгоняющий) зависит от взаимного расположения перицентра

орбиты и радиуса номинальной орбиты. Для этого вычисляется (r( = r( - rн.

Возможны ситуации:

1) (r( < 0 - прикладывается разгоняющий импульс

2) (r( > 0 - прикладывается тормозящий импульс

КА долетает до апоцентра и в апоцентре прикладывается корректирующий

импульс. Время работы КДУ - 20 сек.

Так как время работы КДУ ограничено, а (Vк может быть большим,

следовательно, далее рассчитывается максимальный импульс скорости (Vmax за

20 сек работы двигателя:

(Vmax = Pt/m = 25(20/597 = 0,8375 м/с

Если (Vк > (Vmax в апоцентре прикладывается импульс (Vк = (Vmax. В

результате этого r( немного корректируется. На следующем витке опять

рассчитывается (Vк, и если на этот раз (Vк < (Vmax, в апоцентре

прикладывается импульс (Vк. КДУ включается не на полную мощность P =

((Vк/(Vmax)Pmax.

Время включения = 20 сек.

Это происходит до тех пор, пока не приблизится к r( с заданной

точностью.

После того, как скорректирован перицентр, начинается коррекция апоцентра.

Рассчитываются параметры орбиты и нужный корректирующий импульс, такой,

чтобы r( = rн = 6952137 м. Направление корректирующего импульса также

зависит от величин r( и rн.

Вычисляется (r( = r( - rн.

Возможна ситуация:

(r( > 0 - в перицентре прикладывается тормозящий импульс.

КА долетает до перицентра и в перицентре прикладывается корректирующий

импульс. Время работы КДУ - 20 сек.

Так как время работы КДУ ограничено, а (Vк может быть большим,

следовательно, далее рассчитывается максимальный импульс скорости (Vmax за

20 сек работы двигателя:

(Vmax = Pt/m = 25(20/597 = 0,8375 м/с

Если (Vк > (Vmax, в перицентре прикладывается импульс (Vк = (Vmax. В

результате этого немного корректируется r(. На следующем витке опять

рассчитывается (Vк, и если на этот раз (Vк < (Vmax, в перицентре

прикладывается импульс (Vк. КДУ включается не на полную мощность P =

((Vк/(Vmax)Pmax.

Время включения = 20 сек.

Это происходит до тех пор, пока r( не приблизится к rн с заданной

точностью.

Таким образом осуществляется коррекция перехода.

б) Коррекция наклонения.

После коррекции периода проводятся внешне-траекторные измерения и

получают вектор состояния КА. Если снова необходима коррекция периода ее

проводят еще раз и снова измеряют вектор состояния КА.

Далее проводится коррекция наклонения по такой же схеме. Коррекция

производится в точке пересечения орбиты КА с линией узлов.

После того, как рассчитаны корректирующие импульсы скорости, по формулам

перехода проекции вектора на оси абсолютной системы координат. Далее

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

движения центра масс КА. После этого уравнения интегрируются методом Рунге-

Кутта 5-го порядка с переменным шагом.

Графики изменения элементов орбиты в процессе коррекции приведения

приведены на рис.19-30.

2.5.2. РАСЧЕТ ПОТРЕБНОГО ТОПЛИВА

Масса топлива, необходимого для проведения коррекции траектории

рассчитывается по формуле Циолковского:

m = m0(1 - e-(Vк/W)

m0 = 597 кг - начальная масса МКА (кг)

W = 2200 м/с - скорость истечения газов из сопла двигателя.

Результаты проведения коррекции приведения:

| |tн, с |tк, с |(t, |(Vк, |Имп. |m, кг|

| | | |с |м/c | | |

|Коррекция периода |176242 |262592 |300 |12,1 |15 |3,26 |

|Коррекция |273984 |432298 |580 |24,11|29 |6,48 |

|наклонения | | | | | | |

2.5.3.КОРРЕКЦИЯ ПОДДЕРЖАНИЯ

Основная задача МКА - проведение съемки определенных районов Земли по

крайней мере один раз в сутки, т.е. трасса КА должна проходить над заданным

районом каждые сутки.

Требования для проведения коррекции:

- предельное суточное смещение орбиты по долготе (i = 0,1(

- предельное отклонение наклонения (( = 0,1(.

В пересчете отклонения (( на отклонение по периоду получим:

(T = 1,597 сек. - максимальное отклонение по периоду.

При помощи программы моделирования было просчитано 3 месяца и получено,

что средний период изменился на 3,2 сек, а наклонение - на 0,001(.

Таким образом, коррекцию периода надо делать примерно 1 раз в 1,5 мес.

Нужный импульс скорости - 1 м/с за время активного существования - 5 лет

- коррекцию периода надо провести 40 раз, (V = 40 м/с, масса топлива = 10,8

кг.

За 5 лет (i = 0,02( - коррекцию наклонения проводить не надо.

Графики изменения элементов орбиты за 3 месяца приведены на рис.31-42.

2.6. ДВИЖЕНИЕ МКА ОТНОСИТЕЛЬНО ЦЕНТРА МАСС

2.6.1. УРАВНЕНИЯ ДВИЖЕНИЯ ОТНОСИТЕЛЬНО ЦМ КА

При рассмотрении движения относительно ЦМ КА используют уравнения Эйлера:

Jx(x + (Jz-Jy)(y(z = Mxy + Mxв

Jy(y + (Jx-Jz)(x(z = Myy + Myв

Jz(z + (Jy-Jx)(y(x = Mzy + Mzв

где Jx, Jy, Jz - главные моменты инерции,

My - управляющий момент,

Mв - возмущающий момент.

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

произведением угловых скоростей, значит, уравнения Эйлера имеют вид:

Jx(x = Mxy + Mxв

Jy(y = Myy + Myв

Jz(z = Mzy + Mzв

Главные моменты инерции:

Jx = 532 кг(м2, Jy = 563 кг(м2, Jz = 697 кг(м2.

Центробежные моменты инерции принимаются равными 0.

Возмущающий момент Mв возникает из-за того, что двигатель коррекции

расположен не в центре масс КА, и реактивная тяга, линия действия которой

находится на удалении (плече) l от центра масс КА, создает паразитный

крутящий момент Mв.

Mв = P(l,

где P = 25 H - тяга корректирующего двигателя,

l = 4 мм - плечо.

Таким образом, Mв = 25(0,0004 = 0,1 Нм.

2.6.2. СТАБИЛИЗАЦИЯ УГЛОВОГО ПОЛОЖЕНИЯ ПРИ КОРРЕКЦИИ

Основное требование, предъявляемое в этом режиме:

- точность поддержания направления импульса коррекции - не хуже 1

угл.мин.

Целью данной главы является исследование динамики системы при

стабилизации углового положения при коррекции.

Функциональная схема МКА состоит из следующих эелементов:

1) МКА - малый космический аппарат.

МКА описывается как абсолютно твердое тело.

2) ДУС - датчик угловой скорости.

В качестве ДУС используется командный гироскопический прибор. Он

описывается колебательным звеном с параметрами T = 1/30 c-1 и e = 0,7, а

также нелинейным звеном с насыщением 2(/сек.

3) АЦП - аналогово-цифровой преобразователь.

Преобразует аналоговый сигнал с ДУС в цифровой сигнал.

4) ЦАП - цифро-аналоговый преобразователь.

Преобразует цифровой сигнал с ЦВМ в аналоговый.

5) ШИМ - широтно-импульсный модулятор.

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

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

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

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