Исследование движения центра масс малого космического аппарата (МКА)

Размещено на http://www.allbest.ru/

Оглавление

Введение

1. Краткие сведения об орбите

2. Характеристика орбиты

2.1 Связь МКА с наземными пунктами управления

2.2 Выведение на рабочую орбиту

2.3 Исходные данные и цели работы

2.3.1 Исходные данные

2.3.2 Цели работы

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

2.4.1 Уравнения движения МКА

2.4.2 Возмущающие ускорения, действующие на МКА

2.4.3 Расчет параметров текущей орбиты МКА

2.5 Коррекция траектории МКА

2.5.1 Коррекция приведения

2.5.2 Расчет потребного топлива

2.5.3 Коррекция поддержания

2.6 Движение МКА относительно центра масс

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

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

3. Организационно-экономическая часть

3.1 Организация и планирование выполнения темы

3.2 Определение затрат труда

3.3 Расчет сметы затрат на разработку программного продукта

4. Промышленная экология и безопасность

4.1 Введение

4.2 Анализ вредных факторов

4.3 Требования к видеотерминальным устройствам

4.4 Расчет вредных излучений

4.5 Рациональная организация рабочего места

4.6 Рекомендации по снижению утомляемости

4.7 Защита от напряжения прикосновения. Зануление

4.8 Пожарная безопасность

Список литературы

Приложение. Тексты Программ для Borland C++ и Matlab 4.0 for windows

1. Основной программный модуль main.cpp

2. Подпрограмма расчета возмущающих ускорений, параметров орбиты и коррекции sfun.cpp

3. Файл начальной инициализации init.h

4. Файл описания переменных def.h

5. Файл sfun.h

6. Программа построения временных диаграмм

Введение

В данной работе проводится исследование движения центра масс МКА под действием различных возмущающих ускорений (от нецентральности гравитационного поля Земли, сопротивления атмосферы, притяжения Солнца и Луны, из-за давления солнечных лучей) и создание математической модели движения ЦМ МКА, позволяющей учесть при интегрировании уравнений движения ЦМ МКА эволюцию орбиты МКА.

В работе разрабатывается алгоритм коррекции, ликвидирующий ошибки выведения МКА и рассчитывается масса топлива, необходимая для проведения коррекции, необходимой из-за эволюции параметров орбиты и из-за ошибок выведения МКА на рабочую орбиту.

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

В работе приводятся программы, реализующие интегрирование уравнений движения ЦМ МКА, процесс осуществления коррекции и расчет топлива для коррекции.

1. Краткие сведения об орбите

Основными показателями эффективности космической группировки, являются:

- предельная производительность МКА в сутки на освещенной стороне Земли не менее 400-500 объектов.

- периодичность наблюдения районов съемки не реже одного раза в сутки.

Расположение плоскости орбиты по отношению к Солнцу выбрано таким образом, чтобы угол между линией узлов и следом терминатора на плоскости экватора Земли составлял D_т = 30°. При этом северный полувиток орбиты должен проходить над освещенной частью земной поверхности. Для определенности углу D_т приписывается знак «+» в том случае, если восходящий узел орбиты находится над освещенной частью Земли, и знак «-», если ВУ находится над неосвещенной частью. При выборе баллистического построения оперируют углом D, однозначно определяющимся прямым восхождением Солнца a₀ и долготой восходящего узла орбиты в абсолютном пространстве W: D = a₀ - W. Соотношение между углом D_т и углом D: D є D_т - 90°.

2. Характеристика орбиты

Для решения задач наблюдения Земли из космоса с хорошим разрешением при жестких ограничениях на массу КА и минимизации затрат на выведение целесообразно использовать низкие круговые орбиты. В этом классе орбит выделяют солнечно-синхронные орбиты со следующими свойствами:

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

- близость наклонения плоскости орбиты к полярному, что обеспечивает глобальность накрытия полюсами обзора поверхности Земли.

- возможность наблюдения районов на поверхности Земли примерно в одно и то же местное время при незначительном изменении углов места Солнца в точке наблюдения.

Всем этим условиям удовлетворяют солнечно-синхронные орбиты с высотами от нескольких сот до полутора тысяч километров. На больших высотах наклонение солнечно-синхронной орбиты отличается от полярного, и глобальность накрытия поверхности Земли не обеспечивается. Для повышения эффективности наблюдения целесообразно выбрать орбиты с изомаршрутной трассой, у которых следы орбит ежесуточно проходят на одними и теми же районами Земли, что позволяет обеспечивать периодичность наблюдения одного и того же объекта, как минимум, раз в сутки с одного КА.

Предварительные расчеты показали, что целесообразно использовать орбиту с высотой Н = 574 км и наклонением плоскости орбиты к плоскости экватора Земли i = 97,6°.

Масса МКА может составить от 500 до 800 кг (что зависит от вида целевой аппаратуры, устанавливаемой на борту МКА). Для выведения МКА на орбиту используется РН СС-19 («Рокот») с разгонным блоком «Бриз».

2.1 Связь МКА с наземными пунктами управления

Управление МКА осуществляется с наземных пунктов управления на территории России. Их количество и место расположения выбирается таким образом, чтобы на любом витке можно было организовать сеанс связи с МКА хотя бы с одного пункта управления. Угол возывшения МКА над горизонтом наземного пункта управления должен быть не менее 7°, а дальность до МКА не должна превышать 2200 км.

В расчете зон связи были использованы следующие исходные данные:

- высота орбиты - 574 км.

- наклонение орбиты - 97,6°.

- географическая долгота восходящего узла первого витка - 4° в.д.

- минимальный угол возвышения МКА над местным горизонтом - 7°.

Из рассматривавшихся возможных наземных пунктов управления (Москва, Новосибирск, Хабаровск, Мурманск, Калининград, Диксон, Комсомольск-на-Амуре, Петропавловск-Камчатский), было выбрано три (Москва, Диксон, Петропавловск-Камчатский), обеспечивающие возможности связи с МКА на любом витке орбиты. При этом зоны связи с МКА составляют от 3 до 9 минут на витке.

Интергральные характеристики возможности связи с МКА:

- высота орбиты - 574 км.

- число витков, видимых из Москвы, вит/сутки - 6.

- суммарное время видимости из Москвы, мин - 41.

- суммарное время видимости с трех пунктов, мин - 153.

- максимальное время видимости одного витка, мин - 9,1.

2.2 Выведение на рабочую орбиту

Выведение МКА на орбиту с наклонением i = 97,6° и высотой Н = 574 км осуществляется ракетой-носителем «Рокот» с разгонным блоком «Бриз». При выведении для каждой отделяющейся части РН (отработанная первая ступень, обтекатель, отработанная вторая ступень) существует свой район падения.

Возможные варианты старта:

1. Полигон Байконур.

Из-за отсутствия зон падения отделяющихся частей возможно сформировать опорную орбиту с наклонением i порядка 65°. Для формирования опорной орбиты с наклонением близким полярному при использовании трассы с азимутом стрельбы более 180° (направление стрельбы на юг) - первая ступень падает в районе Ашхабада, обтекатель сбрасывается на высоте Н порядка 100 км, вторая ступень падает за Аравийским полуостровом. С точки зрения энергетики, выведение осуществляется не по оптимальной схеме, в результате чего на круговую орбиту высотой Н порядка 700 км выводится МКА массой менее 600 кг.

2. Полигон Ледяная (Свободный).

Из-за отсутствия зон падения отделяющихся частей возможно сформировать опорную орбиту с наклонением i порядка 54° и 65°. При северном запуске РН первая ступень падает в районе заповедника в устье реки Олейма (приток Лены).

3. Космодром Плесецк.

Азимуты пуска с космодрома Плесецк обеспечивают наклонения орбит i от 72° до 93°. Формирование требуемого наклонения i = 97,6° осуществляется с помощью разгонного блока «Бриз».

В результате работы двух ступеней РН формируется баллистическая траектория с наклонением i = 93°. Высота в момент окончания работы двигателя второй ступени составляет Н = 190 км, наклонная дальность L = 300 км. Приблизительно через 1,2 секунды после прохождения команды на выключение двигателя второй ступени проходит команда на запуск ДУ РБ. После выключения двигателя второй ступени РН происходит отделение от ракеты связки РБ с КА. Время расцепки t = 318 секунд. Абсолютная скорость в момент отделения V = 5550 м/с. Отделяемая масса 6700 кг.

Двигательная установка РБ «Бриз» выполняет задачу доразгона КА при формировании опорной орбиты.

Характеристики двигателя РБ «Бриз»:

- тяга R, кг - 2000.

- удельный импульс R_уд, сек - 324.

- количество включений, р - 7.

- интервал между включениями, сек - 20.

- время функционирования, час - 7.

В результате работы двигателя РБ «Бриз» при первом включении происходит увеличение высоты баллистической траектории с Н = 190 км до Н = 270 км и к моменту окончания работы двигателя (t = 905,5 сек) в точке с аргументом широты u = 104,1° формируется опорная эллиптическая орбита с параметрами:

- высота в перигее Н_п, км - 190.

- высота в апогее Н_а, км - 574.

- большая полуось орбиты а, км - 6747.

- эксцентриситет e - 0,02548

- наклонение i, ° - 93,4.

- период обращения Т, час - 1,53.

- аргумент перигея w, ° - 128,38.

- долгота восходящего узла в гринвичской СК, фиксированной на момент старта W_г, ° - 48,37.

Величина импульса характеристической скорости, отрабатываемого при первом включении ДУ РБ dV₁ = 2,36 км/с, время работы порядка 600 сек.

Работа двигателя при первом включении происходит вне зоны видимости НПУ на территории России. Географические координаты, соответствующие этому моменту:

- широта j » 76°.

- долгота l » 238°.

В момент прохождения МКА перигея опорной эллиптической орбиты (t = 1231 сек) географические координаты составляют:

- широта j » 53°.

- долгота l » 227°.

На опорной эллиптической орбите МКА совершает пассивный полет до апогея. В районе апогея (t = 1,12 час) осуществляется второе включение ДУ РБ.

В результате приложения второго компланарного импульса характеристической скорости dV₂ = 0,12 км/с, при втором включении (время работы 20 сек) формируется круговая орбита с параметрами:

- высота Н, км - 574.

- наклонение i, ° - 93,4.

- период обращения Т, час - 1,6.

Работа двигателя при втором включении происходит вне зоны видимости НПУ на территории России. Географические координаты, соответствующие этому моменту:

- широта j » 1,5°.

- долгота l » 35,8°.

Для создания круговой, солнечно-синхронной орбиты необходимо изменить наклонение до i = 97,6°. С этой целью осуществляется третье включение ДУ РБ при первом прохождении восходящего узла орбиты (t = 1,3 час).

В результате приложения ортогонального импульса характеристической скорости dV₃ = 0,62 км/с, при третьем включении (время работы 90 сек) формируется солнечно-синхронная круговая орбита с параметрами:

- высота Н, км - 574.

- наклонение i, ° - 97,6.

- период обращения Т, час - 1,6.

- число оборотов в сутки N - 15.

Работа двигателя при третьем включении происходит вне зоны видимости НПУ на территории России. Географические координаты, соответствующие этому моменту:

- широта j » 0°.

- долгота l » 28,1°.

После выключения двигателя при третьем запуске происходит отделение МКА от РБ «Бриз».

Кинематические параметры в гринвичской СК, фиксированной на момент старта РН и оскулирующие элементы орбиты на момент отделения от РБ:

Параметр	Значение
t, сек	4946,5
X, м	4638800
Y, м	5120280
Z, м	689680
Vx, м/с	241,23
Vy, м/с	-1233
Vz, м/с	7473,5
l, °	28,1
T, c	5761,67
e	0,0009
i, °	97,595
Ra, м	6940000
Rп, м	6952000

2.3 Исходные данные и цели работы

2.3.1 Исходные данные

Номинальная орбита, необходимая для выполнения задач МКА, имеет следующие параметры:

- круговая, e = 0.

- солнечно-синхронная, скорость прецессии линии узлов орбиты W равна скорости обращения Солнца относительно Земли

W = 2p / 365,2422 = 0,0172 рад/сут = 0,98 °/сут.

- изомаршрутная, за сутки МКА совершает целое количество оборотов (n = 15).

Это обеспечивает прохождение МКА над одними и теми же районами в одно и тоже местное время.

- период Т = 5765 с.

- высота орбиты Н = 574 км.

- наклонение орбиты i = 97,6°.

- географическая долгота восходящего узла орбиты l_э = 28,1°.

Долгота восходящего узла в геоцентрической экваториальной (абсолютной) системе координат OXYZ определяется как разность

l_э - s₀,

где s₀ - часовой угол, отсчитывающийся от гринвичского меридиана до оси X, направленной в точку весеннего равноденствия.

Часовой угол зависит от даты старта и выбирается из астрономического ежегодника. В данной задаче для моделирования выбран часовой угол = 0.

Следовательно долгота восходящего узла орбиты W = l_э = 28,1°.

Исходя из ТЗ, начальная точка выведения имеет следующие координаты в гринвичской системе координат, фиксированной на момент старта РН:



Параметр	Значение
t, сек	4946.5
X, м	4638800
Y, м	5120280
Z, м	689506,95
Vx, м/с	241,23
Vy, м/с	-1233
Vz, м/с	7472,65

Элементы орбиты:

l, °	28,1
T, c	5761,67
e	0,0009
i, °	97,595
Ra, м	6940000
Rп, м	6952000

Кинематические параметры в геоцентрической экваториальной системе координат:

t, сек	4946.5
X, м	6137262,9
Y, м	3171846,1
Z, м	689506,95
Vx, м/с	-201,3
Vy, м/с	-1247,03
Vz, м/с	7472,65
l, °	28,1

Точность выведения:

- предельная ошибка по координате (3s) - 7 км.

- предельная ошибка по скорости (3s) - 5 м/с.

Пересчитав ошибку по координате на ошибку по периоду выведения орбиты получим предельную ошибку по периоду DT - 10 сек.

Корреляционная матрица ошибок выведения на момент выведения составляет.

Члены, стоящие на главной диагонали представляют собой квадраты предельных ошибок - (3s)².

K₁₁ = K₂₂ = K₃₃ = (3s)² = 7² = 49 км.

K₄₄ = K₅₅ = K₆₆ = (3s)² = 5² = 25 м/с.

Остальные члены представляют собой вторые смешанные моменты K_ij = K_ji = r_ijs_is_j или K_ij = K_ji = r_jj(3s_i)(3s_j), где r_jj - коэффициенты связи величин i и j. В данном случае вторые смешанные моменты K_ij = K_ji = 0.

Кинематические параметры в геоцентрической экваториальной системе координат на момент выведения с учетом ошибок выведения:

t, сек	4946.5
X, м	6144262,9
Y, м	3178846,1
Z, м	696506,95
Vx, м/с	-206,3
Vy, м/с	-1252,03
Vz, м/с	7477,65
l, °	28,1

Параметры орбиты с учетом ошибок выведения:

l, °	28,13
T, c	5795,7
W, °	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 в некоторой инерциальной системе координат. Движение совершается под действием силы притяжения F_z. Сила F_z для материальной точки m

На основе второго закона Ньютона уравнения движения материальных точек М и m имеют вид:

(1),  (2)

или

(3),  (4)

где p₁ - радиус-вектор, проведенный из начала инерциальной системы координат в точку m.

p₂ - радиус-вектор, проведенный из начала инерциальной системы координат в точку М.

.

Вычитая из уравнения (3) уравнение (4), получим уравнение движения материальной точки m относительно притягивающего центра М:

Так как m<<М, следовательно, можно пренебречь ускорением, которое КА с массой m сообщает притягивающему центру М. Тогда можно совместить начало инерциальной системы координат с притягивающим центром М. Следовательно.

Таким образом, уравнение невозмущенного движения КА относительно притягивающего центра М в инерциальной системе координат, центр которой находится в М,

Рассмотрим возмущенное движение КА в геоцентрической экваториальной (абсолютной) системе координат OXYZ:

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

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

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

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

Движение КА в абсолютной системе координат OXYZ происходит под действием центральной силы притяжения Земли F_z, а также под действием возмущающих сил F_в. Уравнение движения имеет вид

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

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

где ax_в, ay_в, az_в - возмущающиеся ускорения.

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

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

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

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

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

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

2.4.2 Возмущающие ускорения, действующие на МКА

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

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

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

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

Составляющие типа (m_z/r)(r₀/r)ⁿc_n0P_n0(sinj) - называют зональными гармониками n-порядка. Т.к. полином Лежандра n-го порядка имеет n действительных корней, функция P_n0(sinj) будет менять знак на n широтах, сфера делится на n+1 широтную зону, где эти составляющие имеют попеременно «+» или «-» значения. Поэтому их называют зональными гармониками.

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

(m_z/r)(r₀/r)ⁿc_nmcos(mL)P_nm(sinj) и (m_z/r)(r₀/r)ⁿd_nmsin(mL)P_nm(sinj)

- называют тессеральными гармониками n-порядка и степени m. Они обращаются в 0 на 2m меридианах, где cos(mL) = 0 и sin(mL) = 0 и на n-m параллелях, где P_nm(sinj) = 0 или d^mP_nm(sinj)/d(sinj)^m = 0, сфера делится на n+m+1 трапецию, где эти составляющие сохраняют знак.

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

(m_z/r)(r₀/r)ⁿc_nncos(nL)P_nn(sinj) и (m_z/r)(r₀/r)ⁿd_nnsin(nL)P_nn(sinj)

- называют секториальными гармониками n-порядка и степени m. Эти составляющие меняю знак только на меридианах, cos(nL) = 0 и sin(nL) = 0, на сфере выделяют 2n меридиональных секторов, где эти составляющие сохраняют знак.

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

P_n0(z) = 1/(2ⁿn!)ґ(dⁿ(z² - 1)ⁿ/dzⁿ)

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

P_nm(z) = (1-z²)^m/2ґd^mP_n0(z)/dz^m

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

U_в = U' + DU' = (U - m_z/r) + DU'

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

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

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

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

Зональные гармоники нечетного порядка и тессеральные гармоники, где n-m нечетное число - учитывают ассиметрию Земли относительно плоскости экватора.

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

Первая зональная гармоника имеет порядок 10^-3, а все остальные - порядок 10^-6 и выше. Поэтому будем учитывать в разложении потенциала притяжения только зональную гармонику (n=2, m=0) и секторальную гармонику (n=2, m=2). Также не будем учитывать потенциал аномалий силы тяготения Земли DU'.

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

U_в = (m_z/r)(r₀/r)²[c₂₀P₂₀(sinj) + (c₂₂cos(2L) + d₂₂sin(2L))P₂₂(sinj)],

где c₂₀ = - 0,00109808,

c₂₂ = 0,00000574,

d₂₂ = - 0,00000158.

P₂₀(x) = 1/2²2!ґd²(x² - 1)²/dx².

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

P₂₀(x) = (3x² - 1)/2.

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

P₂₂(x) = (1 - x²)^2/2ґd²P₂₀(x)/dx² = 1/2ґ(1 - x²)ґd²(3x² - 1)/dx²

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

P₂₂(x) = 3(1 - x²).

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

Значит

Чтобы найти возмущающее ускорение от нецентральности поля тяготения Земли в проекциях на оси абсолютной системы координат OXYZ, надо взять производные от возмущающего потенциала U_в по координатам X, Y, Z, причем

r = Ц(x² + y² + z²).

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

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

При движении в атмосфере на КА действует сила аэродинамического ускорения R_x, направленная против вектора скорости КА относительно атмосферы:

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

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

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

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

Так как исследуемая орбита - круговая с высотой Н = 574 км, будем считать, что плотность атмосферы одинакова во всех точках орбиты и равна плотности атмосферы на высоте 574 км. Из таблицы стандартной атмосферы находим плотность наиболее близкую к высоте Н = 574 км. Для высоты Н = 580 км r = 5,098ґ10^-13 кг/м³.

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

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

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

e_a = [V_x/|V|, V_y|V|, V_z/|V|], |V| = Ц(V_x²+V_y² +V_z²)

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

Давление солнечного света учитывается как добавок к постоянной тяготения Солнца - Dm_c. Эта величина вычисляется следующим образом:

Dm_c = pS_мA²/m

где p = 4,64ґ10^-6 Н/м² - давление солнечного света на расстоянии в одну астрономическую единицу А.

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

m - масса КА.

S_м = 8 м² - площадь миделевого сечения - проекция КА на плоскость, перпендикулярную направления солнечных лучей.

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

Dm_c = 1,39154ґ10¹⁵ м³/c².

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

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

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

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

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

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

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

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

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

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

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

| r_c - r| = Ц((x_c-x)²+(y_c-y)²+(z_c-z)²)

где x_c, y_c, z_c - проекции радиуса-вектора Солнца на оси абсолютной системы координат.

Моделирование движения Солнца проводилось следующим образом: за некоторый промежуток времени t Солнце относительно Земли сместится на угол

J = J_н + w_ct,

где J_н = W + (90 - D) - начальное положение Солнца в эклиптической системе координат.

W = 28,1° - долгота восходящего узла первого витка КА.

D = 30° - угол между восходящим узлом орбиты КА и терминатором.

w_c - угловая скорость Солнца относительно Земли.

w_c = 2p/T = 2p/365,2422ґ24ґ3600 = 1,991ґ10^-7 рад/c = 1,14ґ10^-5 °/c

Таким образом, в эклиптической системе координат проекции составляют:

x_c^e = r_ccosJ

y_c^e = r_csinJ

z_c^e = 0

r_c = 1,496ґ10¹¹ м (1 астрономическая единица) - расстояние от Земли до Солнца

Плоскость эклиптики наклонена к плоскости экватора на угол e = 23,45°, проекции r_c на оси абсолютной системы координат можно найти как

x_c = x_c^e = r_ccosJ

y_c^e = y_c^ecose = r_ccosJcose

z_c^e = r_csinJsine

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

a_xc = - m_cx/(Ц((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_yc = - m_cy/(Ц((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_zc = - m_cz/(Ц((x_c-x)²+(y_c-y)²+(z_c-z)²))³

С учетом солнечного давления

a_xc = - (m_c-Dm_c)x/(Ц((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_yc = - (m_c-Dm_c)y/(Ц((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_zc = - (m_c-Dm_c)z/(Ц((x_c-x)²+(y_c-y)²+(z_c-z)²))³

5) Возмущающее ускорение, возникающее из-за влияния Луны.

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

где m_л = 4,902ґ10⁶ м³/c²- постоянная тяготения Луны.

r_л - радиус-вектор от Земли до Луны.

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

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

|r_л - r| = Ц((x_л-x)²+(y_л-y)²+(z_л-z)²)

где 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ґ10⁸ м - среднее расстояние от Земли до Луны

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

a_xл = - m_лx/(Ц((x_л!-x)²+(y_л-y)²+(z_л-z)²))³

a_yл = - m_лy/(Ц((x_л!-x)²+(y_л-y)²+(z_л-z)²))³

a_zл = - m_лz/(Ц((x_л!-x)²+(y_л-y)²+(z_л-z)²))³

Уравнения возмущенного движения при действии корректирующего ускорения имеют вид:

Или

d²x/dt² = - (m_z/r²)x + a_xu + a_xa + a_xc + a_xл + a_xк

d²y/dt² = - (m_z/r²)y + a_yu + a_ya + a_yc + a_yл + a_yк

d²z/dt² = - (m_z/r²)z + a_zu + a_za + a_zc + a_zл + a_zк

2.4.3 Расчет параметров текущей орбиты МКА

Полученная система уравнений движения ЦМ КА интегрируется методом Рунге-Кутта 5-го порядка с переменным шагом. Начальные условия x₀, y₀, z₀, V_x0, V_y0, V_z0 - в абсолютной системе координат, соответствуют начальной точке вывода при учете ошибок выведения. После интегрирования мы получаем вектор состояния КА (x, y, z, V_x, V_y, V_z) в любой момент времени.

По вектору состояния можно рассчитать параметры орбиты. соответствующие этому вектору состояния.

а) Фокальный параметр - р.ъ

р = C²/m_z, где С - интеграл площадей.

C = r ґ V, |C| = C = Ц(C_x²+C_y²+C_z²)

C_x = yV_z - zV_y

C_y = zV_x - xV_z - проекции на оси абсолютной СК

C_z = xV_y - yV_x

б) Эксцентриситет - е.

e = f/m_z, где f - вектор Лапласа

f = V ґ C - m_zr/r, |f| = f = Ц(f_x²+f_y²+f_z²)

f_x = V_yC_z - V_zC_y - m_zx/r

f_y = V_zC_x - V_xC_z - m_zy/r - проекции на оси абсолютной СК

f_z = V_xC_y - V_yC_x - m_zz/r

в) Большая полуось орбиты.

a = p/(1 - e²)

г) Наклонение орбиты - i.

C_x = Csin(i)sinW

C_y = - Csin(i)cosW

C_z = Ccos(i)

можно найти наклонение i = arccos(C_z/C)

д) Долгота восходящего узла - W.

Из предыдущей системы можно найти

sinW = C_x/Csin(i)

cosW = - C_y/Csin(i)

Так как наклонение орбиты изменяется несильно в районе i = 97,6°, мы имеем право делить на sin(i).

Если sinW => 0,

W = arccos (-C_y/Csin(i))

Если sinW < 0,

W = 360 - arccos (-C_y/Csin(i))

е) Аргумент перицентра - w.

f_x = f(coswcosW - sinwsinWcos(i))

f_y = f(coswsinW + sinwcosWcos(i))

f_z = fsinwsin(i)

Отсюда найдем

cosw = f_xcosW/f + f_ysinW/f

sinw = f_z/fsin(i)

Если sinw > 0, w = arccos (f_xcosW/f + f_ysinW/f)

Если sinw < 0, w = 360 - arccos (f_xcosW/f + f_ysinW/f)

ж) Период обращения - Т.

T = 2pЦ(a³/m_z)

Графики изменения элементов орбиты при действии всех, рассмотренных выше, возмущающих ускорений в течение 2-х периодов (Т = 5765 с) приведены на рис. 1-12.

Графики изменения во времени возмущающих ускорений приведены на рис. 13-18.

2.5 Коррекция траектории МКА

Существующие ограничения на точки старта РН и зоны падения отработавших ступеней РН, а также ошибки выведения не позволяют сразу же после пуска реализовать рабочую орбиту. Кроме того, эволюция параметров орбит под действием возмущающих ускорений в процессе полета МКА приводит к отклонению параметров орбиты КА от требуемых значений. Для компенсации воздействия указанных факторов осуществляется коррекция орбиты с помощью корректирующей двигательной установки (КДУ), которая располагается на борту МКА.

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

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

- коррекция приведения - ликвидация ошибок выведения и приведение фактической орбиты к номинальной с заданной точностью.

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

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

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

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

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

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

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

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

2.5.1 Коррекция приведения

После окончания процесса выведения МКА, проводятся внешне-траекторные измерения (ВТИ). Эти измерения обеспечивают, по баллистическим расчетам, знание вектора состояния с требуемой точностью через 2 суток. После этого начинается коррекция приведения.

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

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

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

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

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

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

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

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

Сначала осуществляется коррекция перицентра - приведение текущего расстояния до перицентра r_p к номинальному радиусу r_н = 6952137 м. После измерения вектора состояния рассчитываются параметры орбиты. Далее определяется нужный корректирующий импульс DV_к. Направление импульса (тормозящий или разгоняющий) зависит от взаимного расположения перицентра орбиты и радиуса номинальной орбиты. Для этого вычисляется

Dr_p = r_p - r_н.

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

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

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

КА долетает до апоцентра и в апоцентре прикладывается корректирующий импульс. Время работы КДУ - 20 сек.

Так как время работы КДУ ограничено, а DV_к может быть большим, следовательно, далее рассчитывается максимальный импульс скорости DV_max за 20 сек работы двигателя:

DV_max = Pt/m = 25ґ20/597 = 0,8375 м/с

Если DV_к > DV_max в апоцентре прикладывается импульс DV_к = DV_max. В результате этого r_p немного корректируется. На следующем витке опять рассчитывается DV_к, и если на этот раз DV_к < DV_max, в апоцентре прикладывается импульс DV_к. КДУ включается не на полную мощность

P = (DV_к/DV_max)P_max.

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

Это происходит до тех пор, пока не приблизится к r_p с заданной точностью.

После того, как скорректирован перицентр, начинается коррекция апоцентра. Рассчитываются параметры орбиты и нужный корректирующий импульс, такой, чтобы r_a = r_н = 6952137 м. Направление корректирующего импульса также зависит от величин r_a и r_н.

Вычисляется

Dr_a = r_a - r_н.

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

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

КА долетает до перицентра и в перицентре прикладывается корректирующий импульс. Время работы КДУ - 20 сек.

Так как время работы КДУ ограничено, а DV_к может быть большим, следовательно, далее рассчитывается максимальный импульс скорости DV_max за 20 сек работы двигателя:

DV_max = Pt/m = 25ґ20/597 = 0,8375 м/с

Если DV_к > DV_max, в перицентре прикладывается импульс DV_к = DV_max. В результате этого немного корректируется r_a. На следующем витке опять рассчитывается DV_к, и если на этот раз DV_к < DV_max, в перицентре прикладывается импульс DV_к. КДУ включается не на полную мощность

P = (DV_к/DV_max)P_max.

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

Это происходит до тех пор, пока r_a не приблизится к r_н с заданной точностью.

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

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

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

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

После того, как рассчитаны корректирующие импульсы скорости, по формулам перехода проекции вектора на оси абсолютной системы координат. Далее рассчитывается корректирующее ускорение и подставляется в уравнения движения центра масс КА. После этого уравнения интегрируются методом Рунге-Кутта 5-го порядка с переменным шагом.

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

2.5.2 Расчет потребного топлива

Масса топлива, необходимого для проведения коррекции траектории рассчитывается по формуле Циолковского:

m = m₀(1 - e^-DVк/W)

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

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

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

	tн, с	tк, с	Dt, с	DVк, м/c	Имп.	m, кг
Коррекция периода	176242	262592	300	12,1	15	3,26
Коррекция наклонения	273984	432298	580	24,11	29	6,48

2.5.3 Коррекция поддержания

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

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

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

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

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

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

При помощи программы моделирования было просчитано 3 месяца и получено, что средний период изменился на 3,2 сек, а наклонение - на 0,001°.

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

Нужный импульс скорости - 1 м/с за время активного существования - 5 лет - коррекцию периода надо провести 40 раз, DV = 40 м/с, масса топлива = 10,8 кг.

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

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

2.6 Движение МКА относительно центра масс

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

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

J_xw_x + (J_z-J_y)w_yw_z = M_xy + M_xв

J_yw_y + (J_x-J_z)w_xw_z = M_yy + M_yв

J_zw_z + (J_y-J_x)w_yw_x = M_zy + M_zв

где J_x, J_y, J_z - главные моменты инерции,

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

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

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

J_xw_x = M_xy + M_xв

J_yw_y = M_yy + M_yв

J_zw_z = M_zy + M_zв

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

J_x = 532 кгґм², J_y = 563 кгґм², J_z = 697 кгґм².

Центробежные моменты инерции принимаются равными 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) ШИМ - широтно-импульсный модулятор.

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

Так как динамика ЦАП, АЦП, ШИМ как электронных аналоговых приборов оказывает на систему незначительное влияние по сравнению с динамикой механических (ДУС, двигатели) динамические звенья, описывающие эти элементы, можно заменить соответствующими коэффициентами усиления. В первом приближении значения коэффициентов не принципиально.

6) Двигатель стабилизации.

Двигатель описывается нелинейностью с насыщением 0,127 Нм и звеном запаздывания с Т_д = 0,05 сек.

Тяга двигателя 0,1 Н

7) ЦВМ.

В ЦВМ формируется управление по углу и угловой скорости. Закон управления имеет вид:

e = K(K₁j +K₂j), К = 1, К₁ = 550, К₂ = 430.

Эти коэффициенты подбирались на модели, исходя из требований точности поддержания направления корректирующего импульса, а также длительности переходного процесса.

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

Для разомкнутой системы были построены ЛАЧХ и ФЧХ. Эти графики представлены на рис.43.

Результаты моделирования замкнутой системы представлены на рис.44-46.

Таким образом, в результате моделирования получено, что процесс стабилизации углового положения происходит примерно за 15 сек., статическая точность поддержания углового положения - 0,62 угл.мин., что полностью удовлетворяет требованиям технического задания.

3. Организационно-экономическая часть

3.1 Организация и планирование выполнения темы

Сроки выполнения и затраты на исследования в большой мере зависят от организационных условий выполнения исследовательских работ Поэтому необходимо в первую очередь определить, хотя бы в общем виде, порядок и организацию проведения дипломной работы по заданной теме.

Организация дипломной работы по любой теме складывается из определённых этапов и подэтапов, каждый из которых хотя и может иметь разное содержание, однако структурно занимает равное положение для всех дипломных работ, выполняемых в данной отрасли.

Таким образом, структура дипломной работы может быть сформирована по типовой схеме, упорядоченной в соответствии с конкретным видом исследования. Состав дипломной работы по заданной теме, а также потребные категории исследований по этапам и подэтапам представлены в табл.1.

Этапы	Содержание	Исполнители
1.	Техническое задание	Постановка задачи. Определение состава программного продукта.	Руководитель Разработчик
2.	Эскизный проект	Разработка общего описания программного продукта.	Руководитель Разработчик
3.	Технический проект	Разработка структуры программного продукта.	Разработчик
4.	Рабочий проект	Программирования и отладка программы. Проверка результатов и внесение корректив в программу.	Руководитель Разработчик
5.	Внедрение	Оформление необходимой документации.	Разработчик

3.2 Определение затрат труда

Первым шагом при определении себестоимости программного комплекса является расчет трудоемкости создания и внедрения. Расчет производится по методике, приведенной в документе «Типовые нормы времени на программирование задач для ЭВМ». Типовые нормы времени предназначены для определения затрат времени на разработку программных средств вычислительной техники (ПСВТ).

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

Количество разновидностей форм входной информации - 2,

в том числе:

информации, получаемой от решения смежных задач - 1,

справочной, условно постоянной информации (файл инициализации) 1;

Количество разновидностей форм выходной информации - 2,

в том числе:

печатных документов (временные диаграммы) - 1,

информации, наносимой на магнитные носители (файл инициализации) -1;

Степень новизны комплекса задач - Г (разработка программной продукции, основанной на привязке типовых проектных решений).

Сложность алгоритма - 3 (реализуются стандартные методы решения, не предусмотрено применение сложных численных и логических методов).

Вид используемой информации:

количество разновидностей форм переменной информации (ПИ) - 1, в том числе: информации, получаемой от решения смежных задач - 1;

количество разновидностей форм нормативно-справочной информации (НСИ) (файл инициализации) - 1;

Язык программирования - Borland С++.

Вид представления исходной информации - группа 11 (требуется учитывать взаимовлияние различных показателей).

Вид представления выходной информации - группа 22 (вывод информационных массивов на машинные носители).

Трудоемкость разработки программного продукта t_пп может быть определена как сумма величин трудоемкостей выполнения отдельных стадий разработки программного продукта.

t_пп = t_тз + t_эп + t_тп + t_рп + t_в,

где t_тз - трудоемкость разработки технического задания на создание программного продукта,

t_эп - трудоемкость разработки эскизного проекта программного продукта,

t_тп - трудоемкость разработки технического проекта программного продукта,

t_рп - трудоемкость разработки рабочего проекта программного продукта,

t_в - трудоемкость внедрения программного продукта.

Трудоемкость разработки технического задания рассчитывается по формуле

t_тз = Т^з_рз + Т^з_рп,

где Т^з_рз - затраты времени разработчика постановки задач на разработку ТЗ, чел.-дни,

Т^з_рп - затраты времени разработчика программного обеспечения на разработку ТЗ, чел.-дни.

Значения Т^з_рз и Т^з_рп рассчитываются по формуле

Т^з_рз = t_зК^з_рз, Т^з_рп = t_зК^з_рп,



где t_з - норма времени на разработку ТЗ для программного продукта в зависимости от функционального назначения и степени новизны разрабатываемового программного продукта, чел.-дни (t_з = 29),

К^з_рз - коэффициент, учитывающий удельный вес трудоемкости работ, выполняемых разработчиком постановки задач на стадии ТЗ (К^з_рз = 0,65),

К^з_рп - коэффициент, учитывающий удельный вес трудоемкости работ, выполняемых разработчиком программного обеспечения на стадии ТЗ (К^з_рп = 0,35).

Т^з_рз = 29*0,65 = 18,85 чел.-дней.

Т^з_рз = 29*0,35 = 10,15 чел.-дней.

t_тз = Т^з_рз + Т^з_рп =18,85 + 10,15 = 29 чел.-дней.

Трудоемкость разработки эскизного проекта рассчитывается по формуле

t_эп = Т^э_рз + Т^э_рп,

где Т^э_рз - затраты времени разработчика постановки задач на разработку ЭП, чел.-дни,

Т^э_рп - затраты времени разработчика программного обеспечения на разработку ЭП, чел.-дни.