|
|
|
|
Астрономия>>Исследование движения центра масс межпланетных космических аппаратов 
1. Оглавление 1
2. Исследовательская часть 3
2.1. Введение 3
2.2. Краткие сведения об орбите 4
2.2.1. Характеристика орбиты 4
2.2.2. Связь МКА с наземными пунктами управления 5
2.2.3. Выведение на рабочую орбиту 6
2.3. Исходные данные и цели работы 10
2.3.1. Исходные данные 10
2.3.2. Цели работы 12
2.4. Моделирование движения центра масс МКА 13
2.4.1. Уравнения движения МКА 13
2.4.2. Возмущающие ускорения, действующие на МКА 15
2.4.3. Расчет параметров текущей орбиты МКА 22
2.5. Коррекция траектории МКА 24
2.5.1. Коррекция приведения 24
2.5.2. Расчет потребного топлива 26
2.5.3. Коррекция поддержания 27
2.6. Движение МКА относительно центра масс 28
2.6.1. Уравнения движения относительно центра масс МКА 28
2.6.2. Стабилизация углового положения при коррекции 28
3. ОРГАНИЗАЦИОННО-ЭКОНОМИЧЕСКАЯ ЧАСТЬ 31
3.1. Организация и планирование выполнения темы 31
3.2. Определение затрат труда 31
3.3. Расчет сметы затрат на разработку программного продукта 35
4. ПроМЫШЛЕННАЯ ЭКОЛОГИЯ И БЕЗОПАСНОСТЬ 39
4.1. Введение 39
4.2. Анализ вредных факторов 39
4.3. Требования к видеотерминальным устройствам 44
4.4. Расчет вредных излучений 46
4.5. Рациональная организация рабочего места 46
4.6. Рекомендации по снижению утомляемости 47
4.7. Защита от напряжения прикосновения. Зануление 48
4.8. Пожарная безопасность 49
5. Список литературы 53
6. Приложение. ТекстЫ Программ для Borland C++ и Matlab 4.0 for windows
54
6.1. Основной программный модуль main.cpp 54
6.2. Подпрограмма расчета возмущающих ускорений, параметров орбиты и
коррекции sfun.cpp 57
6.3. Файл начальной инициализации init.h 77
6.4. Файл описания переменных def.h 79
6.5. Файл sfun.h 81
6.6. Файл rk5.h 81
6.7. Программа построения временных диаграмм 82
2. ИССЛЕДОВАТЕЛЬСКАЯ ЧАСТЬ
2.1. ВВЕДЕНИЕ
В данной работе проводится исследование движения центра масс МКА под
действием различных возмущающих ускорений (от нецентральности
гравитационного поля Земли, сопротивления атмосферы, притяжения Солнца и
Луны, из-за давления солнечных лучей) и создание математической модели
движения ЦМ МКА, позволяющей учесть при интегрировании уравнений движения
ЦМ МКА эволюцию орбиты МКА.
В работе разрабатывается алгоритм коррекции, ликвидирующий ошибки
выведения МКА и рассчитывается масса топлива, необходимая для проведения
коррекции, необходимой из-за эволюции параметров орбиты и из-за ошибок
выведения МКА на рабочую орбиту.
Точность проведения коррекции зависит от точности направления
корректирующего импульса, заданной в ТЗ. Было проведено моделирование
системы коррекции в режиме стабилизации углового положения при работе
корректирующей двигательной установки.
В работе приводятся программы, реализующие интегрирование уравнений
движения ЦМ МКА, процесс осуществления коррекции и расчет топлива для
коррекции.
2.2. КРАТКИЕ СВЕДЕНИЯ ОБ ОРБИТЕ
Основными показателями эффективности космической группировки, являются:
- предельная производительность МКА в сутки на освещенной стороне Земли
не менее 400-500 объектов.
- периодичность наблюдения районов съемки не реже одного раза в сутки.
Расположение плоскости орбиты по отношению к Солнцу выбрано таким
образом, чтобы угол между линией узлов и следом терминатора на плоскости
экватора Земли составлял (т = 30(. При этом северный полувиток орбиты
должен проходить над освещенной частью земной поверхности. Для
определенности углу (т приписывается знак «+» в том случае, если восходящий
узел орбиты находится над освещенной частью Земли, и знак «-», если ВУ
находится над неосвещенной частью. При выборе баллистического построения
оперируют углом (, однозначно определяющимся прямым восхождением Солнца (0
и долготой восходящего узла орбиты в абсолютном пространстве (: ( = (0 - (.
Соотношение между углом (т и углом (: ( ( (т - 90(.
2.2.1. ХАРАКТЕРИСТИКА ОРБИТЫ
Для решения задач наблюдения Земли из космоса с хорошим разрешением при
жестких ограничениях на массу КА и минимизации затрат на выведение
целесообразно использовать низкие круговые орбиты. В этом классе орбит
выделяют солнечно-синхронные орбиты со следующими свойствами:
- скорость прецессии плоскости орбиты в пространстве составляет примерно
1( в сутки, что практически обеспечивает постоянство ориентации ее
относительно терминатора Земли в течении всего срока активного
существования КА.
- близость наклонения плоскости орбиты к полярному, что обеспечивает
глобальность накрытия полюсами обзора поверхности Земли.
- возможность наблюдения районов на поверхности Земли примерно в одно и
то же местное время при незначительном изменении углов места Солнца в точке
наблюдения.
Всем этим условиям удовлетворяют солнечно-синхронные орбиты с высотами от
нескольких сот до полутора тысяч километров. На больших высотах наклонение
солнечно-синхронной орбиты отличается от полярного, и глобальность накрытия
поверхности Земли не обеспечивается. Для повышения эффективности наблюдения
целесообразно выбрать орбиты с изомаршрутной трассой, у которых следы орбит
ежесуточно проходят на одними и теми же районами Земли, что позволяет
обеспечивать периодичность наблюдения одного и того же объекта, как
минимум, раз в сутки с одного КА.
Предварительные расчеты показали, что целесообразно использовать орбиту с
высотой Н = 574 км и наклонением плоскости орбиты к плоскости экватора
Земли i = 97,6(.
Масса МКА может составить от 500 до 800 кг (что зависит от вида целевой
аппаратуры, устанавливаемой на борту МКА). Для выведения МКА на орбиту
используется РН СС-19 («Рокот») с разгонным блоком «Бриз».
2.2.3. СВЯЗЬ МКА С НАЗЕМНЫМИ ПУНКТАМИ УПРАВЛЕНИЯ
Управление МКА осуществляется с наземных пунктов управления на территории
России. Их количество и место расположения выбирается таким образом, чтобы
на любом витке можно было организовать сеанс связи с МКА хотя бы с одного
пункта управления. Угол возывшения МКА над горизонтом наземного пункта
управления должен быть не менее 7(, а дальность до МКА не должна превышать
2200 км.
В расчете зон связи были использованы следующие исходные данные:
- высота орбиты - 574 км.
- наклонение орбиты - 97,6(.
- географическая долгота восходящего узла первого витка - 4( в.д.
- минимальный угол возвышения МКА над местным горизонтом - 7(.
Из рассматривавшихся возможных наземных пунктов управления (Москва,
Новосибирск, Хабаровск, Мурманск, Калининград, Диксон, Комсомольск-на-
Амуре, Петропавловск-Камчатский), было выбрано три (Москва, Диксон,
Петропавловск-Камчатский), обеспечивающие возможности связи с МКА на любом
витке орбиты. При этом зоны связи с МКА составляют от 3 до 9 минут на
витке.
Интергральные характеристики возможности связи с МКА:
- высота орбиты - 574 км.
- число витков, видимых из Москвы, вит/сутки - 6.
- суммарное время видимости из Москвы, мин - 41.
- суммарное время видимости с трех пунктов, мин - 153.
- максимальное время видимости одного витка, мин - 9,1.
2.2.4. ВЫВЕДЕНИЕ МКА НА РАБОЧУЮ ОРБИТУ
Выведение МКА на орбиту с наклонением 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.
- эксцентриситет ( - 0,02548
- наклонение i, ( - 93,4.
- период обращения Т, час - 1,53.
- аргумент перигея w, ( - 128,38.
- долгота восходящего узла в гринвичской СК, фиксированной на момент
старта (г, ( - 48,37.
Величина импульса характеристической скорости, отрабатываемого при первом
включении ДУ РБ dV1 = 2,36 км/с, время работы порядка 600 сек.
Работа двигателя при первом включении происходит вне зоны видимости НПУ
на территории России. Географические координаты, соответствующие этому
моменту:
- широта ( ( 76(.
- долгота ( ( 238(.
В момент прохождения МКА перигея опорной эллиптической орбиты (t = 1231
сек) географические координаты составляют:
- широта ( ( 53(.
- долгота ( ( 227(.
На опорной эллиптической орбите МКА совершает пассивный полет до апогея.
В районе апогея (t = 1,12 час) осуществляется второе включение ДУ РБ.
В результате приложения второго компланарного импульса характеристической
скорости dV2 = 0,12 км/с, при втором включении (время работы 20 сек)
формируется круговая орбита с параметрами:
- высота Н, км - 574.
- наклонение i, ( - 93,4.
- период обращения Т, час - 1,6.
Работа двигателя при втором включении происходит вне зоны видимости НПУ
на территории России. Географические координаты, соответствующие этому
моменту:
- широта ( ( 1,5(.
- долгота ( ( 35,8(.
Для создания круговой, солнечно-синхронной орбиты необходимо изменить
наклонение до i = 97,6(. С этой целью осуществляется третье включение ДУ РБ
при первом прохождении восходящего узла орбиты (t = 1,3 час).
В результате приложения ортогонального импульса характеристической
скорости dV3 = 0,62 км/с, при третьем включении (время работы 90 сек)
формируется солнечно-синхронная круговая орбита с параметрами:
- высота Н, км - 574.
- наклонение i, ( - 97,6.
- период обращения Т, час - 1,6.
- число оборотов в сутки N - 15.
Работа двигателя при третьем включении происходит вне зоны видимости НПУ
на территории России. Географические координаты, соответствующие этому
моменту:
- широта ( ( 0(.
- долгота ( ( 28,1(.
После выключения двигателя при третьем запуске происходит отделение МКА
от РБ «Бриз».
Кинематические параметры в гринвичской СК, фиксированной на момент старта
РН и оскулирующие элементы орбиты на момент отделения от РБ:
|Параметр |Значение |
|t, сек |4946,5 |
|X, м |4638800 |
|Y, м |5120280 |
|Z, м |689680 |
|Vx, м/с |241,23 |
|Vy, м/с |-1233 |
|Vz, м/с |7473,5 |
|(, ( |28,1 |
|T, c |5761,67 |
|e |0,0009 |
|i, ( |97,595 |
|Ra, м |6940000 |
|Rп, м |6952000 |
2.3. ИСХОДНЫЕ ДАННЫЕ И ЦЕЛИ РАБОТЫ
2.3.1. ИСХОДНЫЕ ДАННЫЕ
Номинальная орбита, необходимая для выполнения задач МКА, имеет следующие
параметры:
- круговая, e = 0.
- солнечно-синхронная, скорость прецессии линии узлов орбиты ( равна
скорости обращения Солнца относительно Земли
( = 2( / 365,2422 = 0,0172 рад/сут = 0,98 (/сут.
- изомаршрутная, за сутки МКА совершает целое количество оборотов (n =
15).
Это обеспечивает прохождение МКА над одними и теми же районами в одно и
тоже местное время.
- период Т = 5765 с.
- высота орбиты Н = 574 км.
- наклонение орбиты i = 97,6(.
- географическая долгота восходящего узла орбиты (э = 28,1(.
Долгота восходящего узла в геоцентрической экваториальной (абсолютной)
системе координат OXYZ определяется как разность
(э - s0,
где s0 - часовой угол, отсчитывающийся от гринвичского меридиана до оси
X, направленной в точку весеннего равноденствия.
Часовой угол зависит от даты старта и выбирается из астрономического
ежегодника. В данной задаче для моделирования выбран часовой угол = 0.
Следовательно долгота восходящего узла орбиты ( = (э = 28,1(.
Исходя из ТЗ, начальная точка выведения имеет следующие координаты в
гринвичской системе координат, фиксированной на момент старта РН:
|Параметр |Значение |
|t, сек |4946.5 |
|X, м |4638800 |
|Y, м |5120280 |
|Z, м |689506,95 |
|Vx, м/с |241,23 |
|Vy, м/с |-1233 |
|Vz, м/с |7472,65 |
Элементы орбиты:
|(, ( |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 |
|(, ( |28,1 |
Точность выведения:
- предельная ошибка по координате (3() - 7 км.
- предельная ошибка по скорости (3() - 5 м/с.
Пересчитав ошибку по координате на ошибку по периоду выведения орбиты
получим предельную ошибку по периоду (T - 10 сек.
Корреляционная матрица ошибок выведения на момент выведения составляет:
[pic]
Члены, стоящие на главной диагонали представляют собой квадраты
предельных ошибок - (3()2.
K11 = K22 = K33 = (3()2 = 72 = 49 км.
K44 = K55 = K66 = (3()2 = 52 = 25 м/с.
Остальные члены представляют собой вторые смешанные моменты Kij = Kji =
rij(i(j или Kij = Kji = rjj(3(i)(3(j), где rjj - коэффициенты связи величин
i и j. В данном случае вторые смешанные моменты Kij = Kji = 0.
Кинематические параметры в геоцентрической экваториальной системе
координат на момент выведения с учетом ошибок выведения:
|t, сек |4946.5 |
|X, м |6144262,9 |
|Y, м |3178846,1 |
|Z, м |696506,95 |
|Vx, м/с |-206,3 |
|Vy, м/с |-1252,03 |
|Vz, м/с |7477,65 |
|(, ( |28,1 |
Параметры орбиты с учетом ошибок выведения:
|(, ( |28,13 |
|T, c |5795,7 |
|(, ( |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 в некоторой
инерциальной системе координат. Движение совершается под действием силы
притяжения Fz. Сила Fz для материальной точки m определяется формулой:
[pic],
где ( - постоянная притяжения,
ro - единичный вектор, направленный от М к m,
[pic],
где [pic]- радиус-вектор, проведенный из т.М до т.m.
r - относительное расстояние от М до m.
На точку М действует сила Fz, равная по величине и направленная в
противоположную сторону.
На основе второго закона Ньютона уравнения движения материальных точек М
и m имеют вид:
[pic](1), [pic] (2)
или
[pic](3), [pic] (4)
где p1 - радиус-вектор, проведенный из начала инерциальной системы
координат в точку m.
p2 - радиус-вектор, проведенный из начала инерциальной системы координат
в точку М.
[pic].
Вычитая из уравнения (3) уравнение (4), получим уравнение движения
материальной точки m относительно притягивающего центра М:
[pic][pic]
Так как m>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, ( = arccos (fxcos(/f + fysin(/f)
Если sin( 0 - прикладывается тормозящий импульс
КА долетает до апоцентра и в апоцентре прикладывается корректирующий
импульс. Время работы КДУ - 20 сек.
Так как время работы КДУ ограничено, а (Vк может быть большим,
следовательно, далее рассчитывается максимальный импульс скорости (Vmax за
20 сек работы двигателя:
(Vmax = Pt/m = 25(20/597 = 0,8375 м/с
Если (Vк > (Vmax в апоцентре прикладывается импульс (Vк = (Vmax. В
результате этого r( немного корректируется. На следующем витке опять
рассчитывается (Vк, и если на этот раз (Vк 0 - в перицентре прикладывается тормозящий импульс.
КА долетает до перицентра и в перицентре прикладывается корректирующий
импульс. Время работы КДУ - 20 сек.
Так как время работы КДУ ограничено, а (Vк может быть большим,
следовательно, далее рассчитывается максимальный импульс скорости (Vmax за
20 сек работы двигателя:
(Vmax = Pt/m = 25(20/597 = 0,8375 м/с
Если (Vк > (Vmax, в перицентре прикладывается импульс (Vк = (Vmax. В
результате этого немного корректируется r(. На следующем витке опять
рассчитывается (Vк, и если на этот раз (Vк
#include
#include
#include
#include "rk5.h"
#include "sfun.h"
#include "init.h"
#include
typedef long double real;
const float g_r = M_PI/180.;
const float r_g = 180./M_PI;
real t_beg;
real t_end;
real dt;
real toler;
int Np;
int Curp;
real dTp;
real mu_z;
real mu_s;
real mu_l;
real m;
real m_t;
real W;
real w_s;
real w_z;
real w_l;
real ww_l;
real xs,ys,zs;
real xl,yl,zl;
real Fz,Fs,Fl,Fa,U20;
real J1,J2,J3;
int nomin;
real par[8];
real parn[8];
real a_p,e_p,p_p,Om_p,i_p,om_p,Rp_p,Ra_p;
real y_main[6];
real prmt[5];
int Fl_u;
real u_last;
int Fl_ka;
int Fl_kp;
int Fl_ki;
int Fl_i;
int Fl_p;
int Fl_a;
int Fl_lu;
int Fl_pkT;
real dl;
real T_vd;
real dRa;
real dRp;
int Sig;
int Sig_a;
real Tkor;
real Tkore;
real Vkor[3];
real akor[3];
int Fl_l0;
int Fl_l1;
int Fl_pki;
real dV_ps;
real dV_as;
real dV_is;
real dV_ss;
ofstream m_y ("m_y.dat");
ofstream m_f ("m_f.dat");
ofstream m_s ("m_s.dat");
ofstream m_l ("m_l.dat");
ofstream m_par ("m_par.dat");
ofstream u_f ("u_f.dat");
ofstream u_par ("u_par.dat");
ofstream k_par ("k_par.dat");
void out_p(real,real *,real*,int,int,real*);
void main()
{
clrscr();
init_m();
real dery[]={ .167, .167, .167, .167, .166, .166};
int ihlf;
int ndim = 6;
Drkgs(prmt,y_main,dery,ndim,ihlf,fct,out_p);
clrscr();
if (ihlf= (dTp*Curp))
{
Curp++;
gotoxy(1,20);
cout parn[7]))
{
Fl_u = 0;
dl = -(w_z-w_s)*(par[6]-parn[6]);
u_par 79000) && (x 300*g_r)
Fl_u = 1;
u_last = par[7];
}
void korr(real& t, real *f, real *)
{
if (t > (Tkor+172800.))
{
if ((fabs(dl) > 0.1*g_r) && (!Fl_ka) && (!Fl_kp) && (!Fl_ki))
{
Fl_kp = 1;
Fl_ka = 0;
Fl_ki = 0;
cout par[7])
da = fabs(par[5]-par[7]-M_PI);
else
da = fabs(par[5]-par[7]+M_PI);
if (da = (T_vd +20))
{
T_vd = 0;
akor[0] = 0;
akor[1] = 0;
akor[2] = 0;
cout 0)
Sig_a = -1;
else
Sig_a = 1;
cout R_n)
{
Rp_p = R_n;
Ra_p = R_t;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1-e_p);
}
else
{
Rp_p = R_t;
Ra_p = R_n;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1+e_p);
}
real dV = Vk-V_t;
real dVmax = 20*25./m;
cout dVmax)
{
akor[0] = Sig_a*(25./m)*f[3]/V_t;
akor[1] = Sig_a*(25./m)*f[4]/V_t;
akor[2] = Sig_a*(25./m)*f[5]/V_t;
cout fabs(dV))
{
dVmax = fabs(dV);
real Vk_r = Sig_a*dVmax+V_t;
real Ra_r = R_t;
real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;
real a_r = Ra_r/(1+e_r);
real p_r = a_r*(1-e_r*e_r);
real Rp_r = a_r*(1-e_r);
cout 0)
Sig_a = -1;
else
Sig_a = 1;
cout R_n)
{
Rp_p = R_n;
Ra_p = R_t;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1-e_p);
}
else
{
Rp_p = R_t;
Ra_p = R_n;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1+e_p);
}
real dV = Vk-V_t;
real dVmax = 20*25./m;
cout dVmax)
{
akor[0] = Sig_a*(25./m)*f[3]/V_t;
akor[1] = Sig_a*(25./m)*f[4]/V_t;
akor[2] = Sig_a*(25./m)*f[5]/V_t;
cout fabs(dV))
{
dVmax = fabs(dV);
real Vk_r = Sig_a*dVmax+V_t;
real Ra_r = R_t;
real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;
real a_r = Ra_r/(1+e_r);
real p_r = a_r*(1-e_r*e_r);
real Rp_r = a_r*(1-e_r);
cout par[5])
teta = 2*M_PI+par[7]-par[5];
else
teta = par[7]-par[5];
real Vr_i = sqrt(mu_z/par[0])*par[1]*sin(teta);
real Vn_i = sqrt(mu_z/par[0])*(1+par[1]*cos(teta));
V_t = sqrt(f[3]*f[3]+f[4]*f[4]+f[5]*f[5]);
Vk_x = -Vn_i*cos(parn[4])*sin(par[3])+Vr_i*cos(par[3]);
Vk_y = Vn_i*cos(parn[4])*cos(par[3])+Vr_i*sin(par[3]);
Vk_z = Vn_i*sin(parn[4]);
Vk = sqrt(Vk_x*Vk_x+Vk_y*Vk_y+Vk_z*Vk_z);
real dV_x = Vk_x-f[3];
real dV_y = Vk_y-f[4];
real dV_z = Vk_z-f[5];
real dV = sqrt(dV_x*dV_x+dV_y*dV_y+dV_z*dV_z);
real dVmax = 20*25./m;
cout dVmax)
{
akor[0] = (25./m)*dV_x/dV;
akor[1] = (25./m)*dV_y/dV;
akor[2] = (25./m)*dV_z/dV;
cout = 0)
par[3] = acos(c_Om);
else
par[3] = 2*M_PI-acos(c_Om);
real c_om = (f1*cos(par[3])+f2*sin(par[3]))/F;
real s_om = f3/(F*sin(par[4]));
if (s_om > 0)
par[5] = acos(c_om);
else
par[5] = 2*M_PI - acos(c_om);
if (par[2] 0)
par[7] = acos(c_u);
else
par[7] = 2*M_PI - acos(c_u);
}
#include "rk5.h"
#include
void Drkgs(real *prmt,real *y,real *dery,int ndim,int& ihlf,
void (*fct)(real &,real*,real*),
void (*out_p)(real,real*,real*,int,int,real*))
{
static real a[] = { 0.5, 0.292893218811345248, 1.70710678118665475,
0.16666666666666667 };
static real b[] = { 2.0, 1.0, 1.0, 2.0 };
static real c[] = { 0.5, 0.292893218811345248, 1.70710678118665475, 0.5 };
real *aux[8];
int i,j,imod,itest,irec,istep,iend;
real delt,aj,bj,cj,r,r1,r2,x,xend,h;
for (i=0; i= 0.0)
{
iend = 1;
if (r > 0.0) h = xend-x;
}
out_p(x,y,dery,irec,ndim,prmt);
if (prmt[4] != 0.0) goto l40;
itest = 0;
l9: istep++;
j = 0;
l10: aj = a[j];
bj =b[j];
cj = c[j];
for (i=0; i 0.0)
{
if (ihlf-10 >= 0)
{
ihlf = 11;
fct(x,y,dery);
goto l39;
}
for (i=0; i 0) goto l39;
ihlf--;
istep = istep/2;
h = h+h;
if (ihlf 0.0)) goto l4;
ihlf--;
istep = istep/2;
h = h+h;
goto l4;
l39: out_p(x,y,dery,ihlf,ndim,prmt);
l40: for (i=0; i
#include
ifstream if_init;
void nex_ln (void);
void init_m()
{
Np = 150;
t_beg = 0;
t_end = 8000000;
dt = 2;
toler = .05;
dTp = (t_end-t_beg)/float(Np);
Curp = 0;
J1 = 532;
J2 = 563;
J3 = 697;
m = 597.;
W = 2200;
mu_z = 3.9858e14;
mu_s = 1.3249e20;
mu_l = 4.9027e12;
w_s = 2*M_PI/(365.2422*24*3600);
w_z = 2*M_PI/(24*3600);
w_l = 2*M_PI/(27.32*24*3600);
ww_l = 2*M_PI/(18.6*365.2422*24*3600);
parn[0] = 6952137.;
parn[1] = 0;
parn[2] = 6952137;
parn[3] = 28.1*g_r;
parn[4] = 97.6*g_r;
parn[5] = 63.1968*g_r;
parn[6] = 5769.;
parn[7] = 5.751*g_r;
Fl_u = 1;
u_last = parn[7];
Fl_ka = 0;
Fl_kp = 0;
Fl_ki = 0;
Fl_p = 0;
Fl_a = 0;
Fl_i = 0;
Fl_pkT = 0;
Tkor = 0;
T_vd = 0;
akor[0] = 0;
akor[1] = 0;
akor[2] = 0;
dV_ps = 0;
dV_as = 0;
dV_is = 0;
dV_ss = 0;
Fl_l0 = 0;
Fl_l1 = 0;
Fl_pki = 0;
real x0 = 6137262.9+7000;
real y0 = 3171846.1+7000;
real z0 = 689506.95+7000;
real Vx0 = -201.288+5;
real Vy0 = -1247.027+5;
real Vz0 = 7472.65+5;
prmt[0] = t_beg;
prmt[1] = t_end;
prmt[2] = dt;
prmt[3] = toler;
prmt[4] = 0.0;
y_main[0] = x0;
y_main[1] = y0;
y_main[2] = z0;
y_main[3] = Vx0;
y_main[4] = Vy0;
y_main[5] = Vz0;
}
void nex_ln (void)
{
char ch;
if_init.get(ch);
while (ch != '\n')
if_init.get(ch);
}
#endif
6.4 ФАЙЛ ОПИСАНИЯ ПЕРЕМЕННЫХ DEF.H
#ifndef _DEFH
#define _DEFH
#include
typedef long double real;
extern const float g_r;
extern const float r_g;
extern int Np;
extern int Curp;
extern real dTp;
extern real t_beg;
extern real t_end;
extern real dt;
extern real toler;
extern real J1,J2,J3;
extern real mu_z;
extern real mu_s;
extern real mu_l;
extern real m;
extern real m_t;
extern real W;
extern real w_s;
extern real w_z;
extern real w_l;
extern real ww_l;
extern real xs,ys,zs;
extern real xl,yl,zl;
extern real Fz,Fs,Fl,Fa,U20;
extern int nomin;
extern real par[8];
extern real parn[8];
extern real a_p,e_p,p_p,Om_p,i_p,om_p,Rp_p,Ra_p;
extern real y_main[6];
extern real prmt[5];
extern int Fl_u;
extern real u_last;
extern int Fl_ka;
extern int Fl_kp;
extern int Fl_ki;
extern int Fl_i;
extern int Fl_p;
extern int Fl_a;
extern int Fl_lu;
extern int Fl_pkT;
extern real dl;
extern real T_vd;
extern real dRa;
extern real dRp;
extern int Sig;
extern int Sig_a;
extern real Vkor[3];
extern real akor[3];
extern real Tkor;
extern real Tkore;
extern real dV_ps;
extern real dV_as;
extern real dV_is;
extern real dV_ss;
extern int Fl_l0;
extern int Fl_l1;
extern int Fl_pki;
#endif
6.5 ФАЙЛ SFUN.H
#ifndef _SFUN
#define _SFUN
#include "def.h"
#include
#include
#include
void out_p(real x,real *y,real*,int,int,real *);
real interpl(real*,real*,int,real);
void fct(real& ,real *y,real *dery);
void par_or(real *,real *);
#endif
6.5 ФАЙЛ RK5.H
#ifndef _RK5
#define _RK5
#include "def.h"
#include
#include
#include "sfun.h"
void Drkgs(real *prmt,real *y,real *dery,int ndim,int& ihlf,
void (*fct)(real&,real*,real*),
void (*out_p)(real,real*,real*,int,int,real*));
#endif
6.6 ПРОГРАММА ПОСТРОЕНИЯ ВРЕМЕННЫХ ДИАГРАММ
clc
g_r = pi/180;
r_g = 180/pi;
load m_y.dat
t = m_y(:,1);
x = m_y(:,2);
y = m_y(:,3);
z = m_y(:,4);
Vx = m_y(:,5);
Vy = m_y(:,6);
Vz = m_y(:,7);
clear m_y;
s_tmp = size(t);
s_m = s_tmp(1);
clear s_tmp;
load m_f.dat
Fz = m_f(:,2);
Fs = m_f(:,3);
Fl = m_f(:,4);
Fa = m_f(:,5);
U20 = m_f(:,6);
clear m_f;
load m_s.dat
xs = m_s(:,2);
ys = m_s(:,3);
zs = m_s(:,4);
clear m_s;
load m_par.dat
p = m_par(:,2);
e = m_par(:,3);
a = m_par(:,4);
Om = m_par(:,5);
i = m_par(:,6);
omg = m_par(:,7);
T = m_par(:,8);
u = m_par(:,9);
clear m_par;
p_n = 6952137.;
e_n = 0;
a_n = 6952137.;
Om_n0 = 28.1*g_r;
i_n = 97.6*g_r;
omg_n = 346.725*g_r;
T_n = 5765;
ws = 2*pi/(365.2422*24*3600);
for j = 1:s_m, tmp(j) = Om_n0+ws*t(j);
end
Om_n = tmp';
clear tmp;
map = [1,1,1];
colormap(map);
plot(t,p,'y-',[min(t) max(t)],[p_n p_n],'r-'), title (' Фокальный параметр
'), grid on;
print -dwin;
pause;
plot(t,p-p_n,'y-'), title (' dp '), grid on;
print -dwin;
pause;
plot(t,e,'y-',[min(t) max(t)],[e_n e_n],'r-'), title (' Эксцентриситет '),
grid on;
print -dwin;
pause;
plot(t,e-e_n,'y-'), title (' de '), grid on;
print -dwin;
pause;
plot(t,a,'y-',[min(t) max(t)],[a_n a_n],'r-'), title (' Большая полуось
орбиты '), grid on;
print -dwin;
pause;
plot(t,a-a_n,'y-'), title (' da '), grid on;
print -dwin;
pause;
plot(t,Om*r_g,'y-',t,Om_n*r_g,'r-'), title (' Долгота восходящего узла '),
grid on;
print -dwin;
pause;
plot(t,Om*r_g-Om_n*r_g,'y-'), title (' dOm '), grid on;
print -dwin;
pause;
plot(t,i*r_g,'y-',[min(t) max(t)],[i_n*r_g i_n*r_g],'r-'), title ('
Наклонение '), grid on;
print -dwin;
pause;
plot(t,i*r_g-i_n*r_g,'y-'), title (' di '), grid on;
print -dwin;
pause;
plot(t,T,'y-',[min(t) max(t)],[T_n T_n], 'r-'), title (' Период '), grid
on;
print -dwin;
pause;
plot(t,T-T_n,'y-'), title (' dT '), grid on;
print -dwin;
pause;
plot3(x,y,z,'b')
axis([min(x) max(x) min(y) max(y) min(z) max(z)])
set(gca,'box','on')
title (' Положение МКА ')
hold on
plt = plot3(0,0,0,'.','erasemode','xor','markersize',24);
dk = ceil(length(y)/2500);
for k = 1:dk:length(y)
set(plt,'xdata',x(k),'ydata',y(k),'zdata',z(k))
drawnow
end
hold off
pause;
plot(t,Fz,'y-'), title (' Гравитация Земли ' ), grid on;
print -dwin;
pause;
plot(t,Fs,'y-'), title (' Гравитация Солнца и солнечное давление '), grid
on;
print -dwin;
pause;
plot(t,Fl,'y-'), title (' Гравитация Луны '), grid on;
print -dwin;
pause;
plot(t,Fa,'y-'), title (' Сопротивление атмосферы '), grid on;
print -dwin;
pause;
plot(t,U20,'y-'), title (' Нецентральность гравитационного поля Земли '),
grid on;
print -dwin;
pause;
plot(t,Fz+Fs+Fl+Fa+U20,'y-'), title (' Суммарное возмущающее ускорение '),
grid on;
print -dwin;
pause;
clear all
clc
g_r = pi/180;
r_g = 180/pi;
p_n = 6952137.;
e_n = 0;
a_n = 6952137.;
Om_n0 = 28.1*g_r;
i_n = 97.6*g_r;
omg_n = 346.725*g_r;
T_n = 5765;
load u_par.dat
t_u = u_par(:,1);
p_u = u_par(:,2);
e_u = u_par(:,3);
a_u = u_par(:,4);
Om_u = u_par(:,5);
i_u = u_par(:,6);
omg_u = u_par(:,7);
T_u = u_par(:,8);
u_u = u_par(:,9);
clear u_par;
load u_f.dat;
Fz_u = u_f(:,2);
Fs_u = u_f(:,3);
Fl_u = u_f(:,4);
Fa_u = u_f(:,5);
U20_u = u_f(:,6);
clear u_f;
s_tmp = size(t_u);
s_m_u = s_tmp(1);
clear s_tmp;
ws = 2*pi/(365.2422*24*3600);
for j = 1:s_m_u, tmp(j) = Om_n0+ws*t_u(j);
end
Om_n_u = tmp';
clear tmp;
plot(t_u,p_u,'y-',[min(t_u) max(t_u)],[p_n p_n],'r-'), title (' Фокальный
параметр '), grid on;
print -dwin;
pause;
plot(t_u,p_u-p_n,'y-'), title (' dp '), grid on;
print -dwin;
pause;
plot(t_u,e_u,'y-',[min(t_u) max(t_u)],[e_n e_n],'r-'), title ('
Эксцентриситет '), grid on;
print -dwin;
pause;
plot(t_u,e_u-e_n,'y-'), title (' de '), grid on;
print -dwin;
pause;
plot(t_u,a_u,'y-',[min(t_u) max(t_u)],[a_n a_n],'r-'), title (' Большая
полуось орбиты '), grid on;
print -dwin;
pause;
plot(t_u,a_u-a_n,'y-'), title (' da '), grid on;
print -dwin;
pause;
plot(t_u,Om_u*r_g,'y-',t_u,Om_n_u*r_g,'r-'), title (' Долгота восходящего
узла '), grid on;
print -dwin;
pause;
plot(t_u,Om_u*r_g-Om_n_u*r_g,'y-'), title (' dOm '), grid on;
print -dwin;
pause;
plot(t_u,i_u*r_g,'y-',[min(t_u) max(t_u)],[i_n*r_g i_n*r_g],'r-'), title ('
Наклонение '), grid on;
print -dwin;
pause;
plot(t_u,i_u*r_g-i_n*r_g,'y-'), title (' di '), grid on;
print -dwin;
pause;
plot(t_u,T_u,'y-',[min(t_u) max(t_u)],[T_n T_n], 'r-'), title (' Период '),
grid on;
print -dwin;
pause;
plot(t_u,T_u-T_n,'y-'), title (' dT '), grid on;
print -dwin;
pause;
clear all
| Для добавления страницы "Исследование движения центра масс межпланетных космических аппаратов" в избранное нажмине Ctrl+D |
|
|
|
|
|
|
