Моделирование движения центра масс МКА

Моделирование движения центра масс МКА. УРАВНЕНИЕ ДВИЖЕНИЯ КА Рассмотрим невозмущенное движение материальных точек М и m в некоторой инерциальной системе координат.

Движение совершается под действием силы притяжения Fz. Сила Fz для материальной точки m определяется формулой, где - постоянная притяжения, ro - единичный вектор, направленный от М к m где - радиус-вектор, проведенный из т. М до т.m. r - относительное расстояние от М до m. На точку М действует сила Fz, равная по величине и направленная в противоположную сторону.

На основе второго закона Ньютона уравнения движения материальных точек М и m имеют вид 1 , 2 или 3 , 4 где p1 - радиус-вектор, проведенный из начала инерциальной системы координат в точку m. p2 - радиус-вектор, проведенный из начала инерциальной системы координат в точку М Вычитая из уравнения 3 уравнение 4 , получим уравнение движения материальной точки m относительно притягивающего центра М Так как m М, следовательно, можно пренебречь ускорением, которое КА с массой m сообщает притягивающему центру М. Тогда можно совместить начало инерциальной системы координат с притягивающим центром М. Следовательно Таким образом, уравнение невозмущенного движения КА относительно притягивающего центра М в инерциальной системе координат, центр которой находится в М, имеет вид, где fM - гравитационная постоянная Земли. Рассмотрим возмущенное движение КА в геоцентрической экваториальной абсолютной системе координат OXYZ - начало О - в центре масс Земли ось X направлена в точку весеннего равноденствия ось Z совпадает с осью вращения Земли и направлена на Северный полюс Земли ось Y дополняет систему до правой.

Движение КА в абсолютной системе координат OXYZ происходит под действием центральной силы притяжения Земли Fz, а также под действием возмущающих сил Fв. Уравнение движения имеет вид или где m 597 кг - масса КА. В проекциях на оси абсолютной системы координат OXYZ получим или или или где axв, ayв, azв - возмущающиеся ускорения. Основные возмущающиеся ускорения вызываются следующими причинами - нецентральностью поля притяжения Земли сопротивлением атмосферы Земли влиянием Солнца влиянием Луны давлением солнечного света. 2.4.2. ВОЗМУЩАЮЩИЕ УСКОРЕНИЯ, ДЕЙСТВУЮЩИЕ НА МКА 1 Возмущающееся ускорение, вызванное нецентральностью гравитационного поля Земли. Рассмотрим потенциал поля притяжения Земли. При точном расчете параметров орбиты спутников, в качестве хорошего приближения к действительной поверхности Земли принимают геоид.

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

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

Во всех точках геоида потенциал притяжения имеет одно и то же значение. Потенциал притяжения Земли можно представить в виде разложения по сферическим функциям. где 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 находится по следующей формуле Pn0 z 1 2nn! dn z2 - 1 n dzn Присоединенная функция Лежандра порядка n и степени m находится по следующей формуле Pnm z 1-z2 m 2dmPn0 z dzm Возмущающая часть гравитационного потенциала Земли равна Uв U U U - z r U где U - потенциал аномалий силы тяготения Земли. U - часть потенциала Земли, которая учитывает несферичность Земли. Следовательно, Первая зональная гармоника в разложении потенциала учитывает полярное сжатие Земли. Зональные гармоники нечетного порядка и тессеральные гармоники, где 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,0574, d22 - 0,0158. 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 2d2P20 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 . Значит Чтобы найти возмущающее ускорение от нецентральности поля тяготения Земли в проекциях на оси абсолютной системы координат OXYZ, надо взять производные от возмущающего потенциала Uв по координатам X, Y, Z, причем r x2 y2 z2 . Следовательно, 2 Возмущающее ускорение, вызванное сопротивлением атмосферы.

При движении в атмосфере на КА действует сила аэродинамического ускорения Rx, направленная против вектора скорости КА относительно атмосферы где Cx 2 - коэффициент аэродинамического сопротивления.

Sм 2,5 м2 - площадь миделевого сечения - проекция КА на плоскость, перпендикулярную направлению скорости полета.

V - скорость КА плотность атмосферы в рассматриваемой точке орбиты. Так как исследуемая орбита - круговая с высотой Н 574 км, будем считать, что плотность атмосферы одинакова во всех точках орбиты и равна плотности атмосферы на высоте 574 км. Из таблицы стандартной атмосферы находим плотность наиболее близкую к высоте Н 574 км. Для высоты Н 580 км 5,09810-13 кг м3. Сила аэродинамического ускорения создает возмущающее касательное ускорение aa Найдем проекции аэродинамического ускорения на оси абсолютной системы координат axa, aya, aza aa направлено против скорости КА, следовательно единичный вектор направления имеет вид ea Vx V , Vy V , Vz V , V Vx2 Vy2 Vz2 Таким образом, Значит 3 Возмущающее ускорение, вызванное давлением солнечного света.

Давление солнечного света учитывается как добавок к постоянной тяготения Солнца - c. Эта величина вычисляется следующим образом c pSмA2 m где p 4,6410-6 Н м2 - давление солнечного света на расстоянии в одну астрономическую единицу А. A 1,4961011 м - 1 астрономическая единица. m - масса КА. Sм 8 м2 - площадь миделевого сечения - проекция КА на плоскость, перпендикулярную направления солнечных лучей.

Таким образом, c 1,391541015 м3 c2. 4 Возмущающее ускорение, возникающее из-за влияния Солнца.

Уравнение движения КА в абсолютной системе координат OXYZ относительно Земли при воздействии Солнца где z - постоянная тяготения Земли. c - постоянная тяготения Солнца. r - радиус-вектор от Земли до КА. rc - радиус-вектор от Земли до Солнца.

Таким образом, возмущающее ускорение, возникающее из-за влияния Солнца. Здесь первое слагаемое есть ускорение, которое получил бы КА, если он был непритягивающим, а Земля отсутствовала.

Второе слагаемое есть ускорение, которое сообщает Солнце Земле, как непритягивающему телу. Следовательно, возмущающее ускорение, которое получает КА при движении относительно Земли - это разность двух слагаемых. Так как rc r, то в первом слагаемом можно пренебречь r. Следовательно 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,2422243600 1,99110-7 рад c 1,1410-5 c Таким образом, в эклиптической системе координат проекции составляют xce rccos yce rcsin zce 0 rc 1,4961011 м 1 астрономическая единица - расстояние от Земли до Солнца Плоскость эклиптики наклонена к плоскости экватора на угол 23,45, проекции rc на оси абсолютной системы координат можно найти как xc xce rccos yce ycecos rccoscos zce rcsinsin Таким образом, проекции возмущающего ускорения на оси абсолютной системы координат 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 относительно Земли при воздействии Луны где л 4,902106 м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 л cos cos i - sin sin i cos л где л - долгота восходящего узла лунной орбиты, отсчитывается от направления на точку весеннего равноденствия угол между плоскостями эклиптики и экватора Земли. Величина л колеблется с периодом 18,6 лет между минимумом при л - i 1818 и максимумом при л i 2836 при 0. Долгота восходящего узла лунной орбиты л изменяется с течением времени t на величину л t360 18,6365,2422243600. Положение Луны на орбите во время t определяется углом л t360 27,32243600. По формулам перехода найдем проекции вектора положения Луны на оси абсолютной системы координат xл rл cosлcosл - cosлsinлsinл yл rл cosлsinл cosлsinлcosл zл rлsinлsinл rл 3,844108 м - среднее расстояние от Земли до Луны Таким образом, проекции возмущающего ускорения на оси абсолютной системы координат 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 Уравнения возмущенного движения при действии корректирующего ускорения имеют вид или 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 coscos - sinsincos i fy f cossin sincoscos i fz fsinsin 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 2520 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 2520 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к, м c Имп. m, кг Коррекция периода 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.