CC BY
7УДК 629.785:523.31-852
Влияние возмущающих аэродинамических сил на эволюцию орбиты
космического аппарата
А.В. Ващенко
В работе исследуется возмущенное движение космического аппарата (КА). Описывается влияние возмущающих аэродинамических сил на эволюцию орбиты КА. В работе представлены результаты расчетов эволюции околокруговых орбит для космических аппаратов различной характеристической площади, массы. Также в работе описывается влияние уровня солнечной активности на эволюцию орбиты КА и зависимость величины возмущающего аэродинамического ускорения от высоты орбиты.
На КА, движущийся на высоте 150-1500 км, заметное влияние оказывает сопротивление атмосферы. Аэродинамическая сила от сопротивления земной атмосферы, действующая на КА имеет следующий вид:
Сх - коэффициент лобового сопротивления КА, Sm - характерная площадь КА, m - масса
При расчетах принимается следующее допущение, что &х принимается постоянным по времени. Это упрощение связанно с наличием солнечных батарей и их переменной ориентацией по Солнцу.
Проекции скоростей имеют вид:
Fa = Sk pV2,
а проекции возмущающих ускорений записываются в виде:
AS = -spVoVos AT = -sxpV0Va
КА.
V -(1 + e cos(u -w)) -w3r cos(i)
Уок = соЗг sin(/') cos(w)
у = 1у 2 + у 2 + у 2
'0 V os ^ ' ot ^ ' 0№
2
В формулах для расчета аэродинамической силы и проекций возмущающих ускорений одним из элементов является р - плотность земной атмосферы. Плотность атмосферы вычисляют по ГОСТ 25645.101-83. Данный ГОСТ устанавливает модель плотности, методику расчета и значения средней плотности верхней атмосферы Земли и её предельных отклонений для различных уровней солнечной активности при неопределенности даты и времени запуска искусственного спутника Земли. Модель плотности верхней атмосферы Земли по ГОСТ 25645.101-83 представляет собой среднегодовую плотность атмосферы как функцию высоты для десяти фиксированных значений индекса солнечной активности: F0=65; 75; 100; 125; 150; 175; 200; 225; 250; 275 (в 10*(-22) Вт/(м2*Гц)). В работе расчет производиться только для минимального - 65*(10*(-22) Вт/(м2*Гц)), среднего - 150*(10*(-22) Вт/(м2*Гц)) и максимального - 275*(10*(-22) Вт/(м2*Гц)) индексов солнечной активности.
Модельная плотность атмосферы вычисляют по формуле:
где р - модельная плотность атмосферы, кг/мЗ;
а0, аisa,1, аisa,2, аisa,3 - коэффициенты модели, используемые для расчета плотности атмосферы при различных значениях F0;
F0 - фиксированное значение индекса солнечной активности F10.7 за рассматриваемый период времени;
F10.7 - индекс солнечной активности равный плотности потока радиоизлучения Солнца на длине волны 10,7 см (на частоте 2800 МГц), выраженной в солнечных единицах потока: 10(-22)
h - геометрическая высота над поверхностью общего земного эллипсоида, км.
В приведенных ниже рисунках 1 и 2 показана зависимость плотности земной от высоты орбиты и уровня солнечной активности (данные показаны для минимального, среднего и максимального уровней солнечной активности: 65,150 и 275 *10(-22)Вт/м2Гц).
ВТ/(м2*Гц);
Зависимость п
Рисунок 1 Зависимость плотности земной атмосферы от уровня солнечной активности и высоты орбиты (от 200 до 400 км.).
Зависимость ги
Рисунок 2 Зависимость плотности земной атмосферы от уровня солнечной активности и высоты орбиты (от 450 до 700 км.).
В программе баллистических расчетов задается условная величина индекса солнечной активности (isa=1,2,3; что соответствует: isa=1 - минимуму солнечной активности F0 =65*10(-22) Вт/м2*Гц, isa=2 - номиналу солнечной активности F0 =150*10(-22) Вт/м2*Гц и isa=3 - максимуму солнечной активности F0 =275*10(-22) Вт/м2*Гц.
В соответствии с этими величинами программа производит присвоение соответствующих значений аisa,1, аisa,2, аisa,3 величинам из матрицы: Ж—15.77005 0.78319 70.58367 Ц
a =
-18.70041 0.57145 110.48925 -20.35393 0.42793 135.74445
- данные по ГОСТ 25645.101-83.
ш
Значение а0 принимается равным 9.8067, так как оно не зависит от уровня солнечной активности (ГОСТ 25645.101-83).
Расчет производится с помощью программы баллистических расчетов (на базе математической системы компьютерной алгебры Maple 7).
В расчетах используется система дифференциальных уравнений в оскулирующих элементах с модификацией - вместо времени прохождения перицентра используется аргумент широты. Данная система может использоваться для анализа движения любого типа (эллипс, парабола, гипербола), но она имеет особенности при е=0 и i=0. Для исключения особенности е=0 используются две новых формальных переменных (составляющие вектора Лапласа): l(t)=e(t)sin(ro(t)) и q(t)=e(t)cos(ro(t)).
Метод оскулирующих элементов сводится к тому, что исследование возмущенной траектории КА может быть сведено к анализу совокупности невозмущенных траекторий, соответствующих каждому моменту времени. Т.е., при анализе возмущенного движения мы можем считать, что в любой момент времени КА находится на дуге конического сечения с определенными значениями элементов орбиты, а коническое сечение при движении КА изменяется. При этом изменяются элементы орбиты [1]:
p=p(t); e=e(t); ro=ro(t); Q=Q(t); i=i(t); u=u(t).
Элементы орбиты, рассматриваемые как функции времени, через которые координаты и составляющие скорости в возмущенном движении выражаются теми же формулами, что и в невозмущенном движении, называются оскулирующими элементами [1].
Рисунок 3 Элементы возмущенного движения
Обозначения: р^) - параметр, е^) - эксцентриситет, ю^) - агрумент широты перицентра, О^) -долгота восходящего узла, ОД - наклонение, и(1) - аргумент широты, A(t) - большая полуось, ОД -радиус, ОД - истинная аномалия, _S(t), _Т^), W(t) проекции возмущающих ускорений на оси орбитальной системы координат (по радиус-вектору, по нормали к нему в плоскости орбиты и по нормали к плоскости орбиты).
После всех преобразований система дифференциальных уравнений движения в оскулирующих элементах при исключении особенностей е^-0 принимает вид:
д ^ _Т(г)Р(г)
зг01 := |- р( t ) = 2 У ^
дгГУ 7 1 + q(t)cos(u(t)) + 1(t) sin(u(t))
st02 := ^ 1(г) = г) ^ cos(u(г))
Т(г) (—а(г) sin(u(г)) ^(и(г)) — 2 1(г) - 2 sm(u(г)) + 1(г) ^(и(г))2)
Лт -_
1 + г) cos(u(г)) + 1(г) sin(u(г))
я(г) _W(г) cot(i(г)) sin(u(г)) ' р(г)
_т
1 + г) г)) + 1(г) sin(u(г))
зг03 := ^ я(г) = г) д/лГ sin(u(г))
(я(г) cos(u(г))2 + я(г) + 1(г) sin(u(г)) г)) + 2 г))) рт) _Т(г) + 1 + я(г) cos(u(г)) + 1(г) sin(u(г))
1(г) _W(г) cot(i(г)) sin(u(г + 1 + q(г)cos(u(г)) + 1(г) sin(u(г))
= а лурпТ _W(г) г))
т04 := — г) = (1 + я(г) cos(u(г)) + 1(г) эт^г))) sin(i(г))
= а _ ^рПг г) cos(u(г))
: ^г1(г1 + г) cos(u(г)) + 1(г) г))
д
згОб := дг u( г) =
_W(г) р(г)2 cot(i(г)) г)) (1 + я(г) г)) + 1(г) sin(u(г)))3лт
Лт I 1— ■ ; Л/Лз . I (1 + я(г)cos(u(г)) + 1(г) sm(u(г)))2
(з/2)
р( г)
д М РтГ А(')2 (1 + г г)) + 1(г) sin(u(г)))_Т(г)
иг07 := — А(г)= 2 V рП-—-
дг р(г)
2 /\Jpmr А(г)2 (sin(u(г)) я(г)— cos(u(г)) 1(г)) _Б(г) + р(0
м08 := — г) = 0 дг
Ш09 := дг »( г) = 0
игою := — Г10( г) = 0 дг
Данная система дифференциальных уравнений используется при исследовании эволюции орбиты КА в программе баллистических расчетов.
Результаты баллистических расчетов
В приведенных ниже рисунках 4 и 5 показана зависимость величины возмущающего аэродинамического ускорения от высоты орбиты.
Характеристики космического аппарата следующие: Масса КА на орбите Наклонение орбиты Характерная площадь аппарата
1950 кг. 75° 20 м2
Рисунок 4 Зависимость величины возмущающего аэродинамического ускорения для орбит высотой от 180 до 650 км.
Рисунок 5 Зависимость величины возмущающего аэродинамического ускорения для орбит высотой от 350 до 650 км.
На рисунках 6 и 7 показаны зависимости изменения большой полуоси и долготы восходящего узла от уровня солнечной активности и высоты орбиты.
Рисунок 6 Зависимость изменения большой полуоси от уровня солнечной активности и высоты орбиты.
Зависимость из* 15 <
1Г ' - £ — □
¿ввш □ □
Рисунок 7 Зависимость изменения долготы восходящего узла от уровня солнечной активности и высоты орбиты.
На рисунках 8 и 9 показаны зависимости изменения большой полуоси и долготы восходящего узла от характеристической площади КА и наклонения орбиты.
Рисунок 8 Зависимость изменения большой полуоси от характеристической площади и наклонения орбиты
Зависимо характе|
Рисунок 9 Зависимость изменения долготы восходящего узла от характеристической площади и наклонения орбиты
На рисунках 10 и 11 показаны зависимости изменения большой полуоси и долготы восходящего узла от массы КА и наклонения орбиты.
Зависимость изм<
Рисунок 10 Зависимость изменения большой полуоси от массы КА и наклонения орбиты
Рисунок 11 Зависимость изменения долготы восходящего узла от массы КА и наклонения орбиты
Список используемой литературы:
1. Дубошин Т.Н. Справочное руководство по небесной механике. М.: Изд.«Наука», 1971.-431 с
2. Дубошин Т.Н. Небесная механика (основные задачи и методы). М.: Изд.«Наука», 1975.-385 с
3. Атмосфера Земли верхняя. Модель плотности для проектных баллистических расчетов искусственных спутников Земли. ГОСТ 25645.101-86. Москва. 1983.-168 с.
4. Гушин В.Н. Основы устройства космических аппаратов. М.:Машиностроение, 2003.-272 с.
5. Константинов М.С., Каменков Е.Ф., Перелыгин Б.П., Безвербый В.К. Механика космического полета. М.:Машиностроение, 1989.-408 с.
Сведения об авторах:
Ващенко Алексей Викторович, аспирант кафедры «Космические системы и ракетостроение» Московского авиационного института (государственного технического университета); email: AlekseyVash@yandex.ru
