Содержание к диссертации
Введение
1 Высокоточное численное моделирование движения систем ИСЗ с использованием параллельных вычислений 19
1.1 Дифференциальные уравнения движения ИСЗ 19
1.2 Описание математической модели сил, действующих на ИСЗ
1.2.1 Возмущения от несферичности геопотенциала 19
1.2.2 Возмущения от приливных деформаций центрального тела 21
1.2.3 Возмущения от третьего тела 27
1.2.4 Возмущения от светового давления 28
1.2.5 Релятивистские эффекты 34
1.2.6 Возмущения от сопротивления атмосферы 1.3 Метод Гаусса-Эверхарта высокого порядка 36
1.4 Оценка эффективности разработанного для кластера ТГУ комплекса программ при решении задач динамики систем ИСЗ 39
1.5 Тестирование комплекса программ на радарных наблюдениях фрагментов распада космических аппаратов 41
2 Комплекс алгоритмов и программ для исследования хаотичности в динамике околоземных объектов 46
2.1 Характеристики хаотичности движения 46
2.2 Алгоритм вычисления параметра MEGNO в задачах динамики ИСЗ 48
2.3 Описание программного комплекса и результатов его тестирования 53
3 Комплекс алгоритмов и программ для решения с применением параллельных вычислений обратных задач динамики исз по данным измерений 55
3.1 Алгоритм решения задачи определения параметров движения по данным измерений55
3.2 Особенности реализации алгоритма в среде параллельных вычислений 56
3.3 Описание программного комплекса и результатов его тестирования 60
4 Численное моделирование влияния слабых возмущений с помощью комплекса программ «численная модель движения систем ИСЗ» 64
4.1 Учет влияния высоких гармоник геопотенциала 64
4.2 Оценка влияния приливных деформаций с использованием различных приливных моделей 67
4.3 Оценка влияния преломления солнечных лучей в земной атмосфере и сжатия Земли на величину возмущений от светового давления 71
4.4 Исследование влияния эффекта Пойнтинга-Робертсона на орбитальную эволюцию объектов геостационарной зоны 72
5 Численное моделирование в решении обратных задач динамики ИСЗ по данным измерений 76
5.1 Определение параметров модели силы светового давления для объектов ГНСС 76
5.2 Результаты исследования орбитального движения фрагмента космического мусора на GEO по данным четырёхлетних оптических наблюдений 80
5.3 Исследование алгоритма быстрого численного оценивания вероятности столкновения двух объектов в околоземном пространстве 85
6 Численное моделирование в задаче исследования динамической структуры околоземного орбитального пространства 92
6.1 Выявление пространственной плотности околоземных объектов, при которой возможно возникновение каскадного эффекта столкновений 92
6.2 Метод ляпуновских характеристик в анализе структуры орбитального пространства зоны GEO 97
6.3 Численно-аналитическая методика выявления вековых резонансов 101
6.4 Исследование динамики ИСЗ Эталон-1 и -2 103
Заключение 108
Список использованных источников
- Возмущения от приливных деформаций центрального тела
- Алгоритм вычисления параметра MEGNO в задачах динамики ИСЗ
- Особенности реализации алгоритма в среде параллельных вычислений
- Результаты исследования орбитального движения фрагмента космического мусора на GEO по данным четырёхлетних оптических наблюдений
Возмущения от приливных деформаций центрального тела
В 1973 г. Э. Эверхарт (Everhart, 1974) предложил интегратор, разработанный специально для численного исследования орбит, и продемонстрировал его высокую эффективность в задачах кометной динамики. Обнаружив принадлежность своего интегратора к семейству интеграторов типа Батчера, Эверхарт обобщил разработанный алгоритм для численного решения любых обыкновенных дифференциальных уравнений первого и второго порядков (Everhart, 1974). Так Э. Эверхарт расширил область применения своего интегратора, который до сих пор не потерял своей популярности в решении задач небесной механики.
В работе (Butcher, 1964) показано, что интегратор Эверхарта (RA15) наследует все замечательные свойства неявных методов Рунге-Кутты типа Батчера, так как он основан на видоизмененных формулах этих методов.
Кроме того, благодаря своей оригинальной схеме, с точки зрения численного интегрирования интегратор имеет следующие преимущества: Алгоритм интегрирования универсален для любого порядка. Интегратор имеет простой критерий для выбора шага интегрирования. В интеграторе реализован довольно точный предиктор решения. Это позволяет сократить количество итераций на шаге в процессе численного интегрирования до двух.
Несмотря на это, код RA15 (Everhart, 1985) и его модификация типа RADAU_27 существенно ограничивают возможности интегратора. По этой причине для численного исследования динамики околоземных объектов мы использовали новый код интегратора Эверхарта, разработанный Авдюшевым В.А. (2006) и названный им GAUSS_15.
В отличие от более ранних версий (RA15 и RADAU_27), код GAUSS_15 позволил устранить следующие трудности:
1. Громоздкий и трудночитаемый код. Использование возможностей современного Фортрана позволило сократить программный код почти в два раза.
2. Устранены все связанные с порядком метода константы (большое количество таких констант затрудняло обобщение кода на другие порядки).
3. Интегратор реализован только для нечетных порядков с разбиением Гаусса-Радо, хотя известно, что неявные методы Рунге-Кутты, построенные на симметричных разбиениях Гаусса-Лобатто и Гаусса-Лежандра, обладают геометрическими свойствами (Hairer, Lubich, Wanner, 2002). 4. Стартовый шаг интегрирования в режиме переменного шага выбирается независимо от дифференциальных уравнений, поэтому не всегда оптимален.
5. Ограничения на величину выбираемого переменного шага заданы в интеграторе из эмпирических соображений и они, кроме того, не связаны с порядком интегратора.
В связи с этим В.А. Авдюшевым был написан новый код интегратора Эверхарта GAUSS34, который разрешает названные выше трудности: 1. Интегрирование на шаге до полной сходимости итерационного процесса. 2. Запоминание величины предпоследнего шага после выполнения процедуры интегрирования, что весьма полезно при многократном использовании программного кода в режиме переменного шага. 3. Быстрый выбор стартового шага, требуемый лишь для первого обращения к интегратору (при повторном обращении используется запоминаемый шаг предыдущего обращения).
Общую теорию интегратора Эверхарта, а также сам программный код GAUSS_15, который был использован для проведения исследований динамики ИСЗ, представленных в данной диссертации, можно найти в работе (Авдюшев, 2006; Авдюшев, 2010).
Используя интегратор Гаусса-Эверхарта, была сделана оценка его возможностей и эффективности в сравнении со старым интегратором RADAU27 (7, 11, 15, 19, 23 и 27 порядков) в возмущенной задаче, где в структуре возмущений учитывались влияния от несферичности Земли (до гармоник 15 порядка) и от притяжения Луны и Солнца. В качестве объекта был рассмотрен ИСЗ системы ГЛОНАСС.
Интегрирование выполнялось на временном интервале 3 месяца с переменным шагом в компьютерной арифметике с четверной точностью на кластере «СКИФ Cyberia». Решение системы нелинейных уравнений на шаге находилось за 2 итерации. После многочисленных вычислений при различных параметрах eto/ были получены характеристики точность-быстродействие. Точность Ах определялась путем сравнения спутниковых эфемерид прямого и обратного интегрирования, а в качестве характеристики быстродействия было принято число обращений к процедуре вычисления правых частей дифференциальных уравнений NF. Результаты для 15-го и 19-го порядков приведены на рисунке 1.6.
Как видно из рисунка 1.6, новая версия интегратора Гаусса-Эверхарта более эффективна, что можно объяснить главным образом использованием в интеграторе модифицированного алгоритма выбора шага. Следует заметить, что именно благодаря повышению вычислительной точности оказывается возможным достичь очень высокой точности интегрирования. При вычислениях же с двойной точностью ошибки интегрирования не удается уменьшить до уровня
Как уже отмечалось, интегратор четного порядка, основанный на разбиении Гаусса-Лобатто, является симметричным и потому обладает геометрическими свойствами, отвечающими свойству обратимости по времени дифференциальных уравнений, к которым также относятся и используемые нами динамические уравнения ИСЗ. Геометрические свойства интеграторов позволяют улучшить поведение глобальной ошибки численного интегрирования. Если для обычных методов ошибка в векторе положения со временем ведет себя квадратично, то при использовании геометрических методов ошибка возрастает линейно, т.е. медленнее.
На рисунке 1.7 показано поведение ошибок интегрирования для четного и нечетного порядков (eto/ = 106 км). В начале интегрирования за счет более высокого порядка результаты интегратора на разбиении Гаусса-Радо значительно точнее, нежели интегратора на разбиении Гаусса-Лобатто. Однако далее, уже после 10 оборотов объекта, за счет лучшего поведения методической ошибки точность интегратора четного порядка становится все более высокой в сравнении с точностью аналога нечетного порядка. Таким образом, с увеличением интервала интегрирования сравнительная эффективность интегратора четного порядка будет только возрастать.
Алгоритм вычисления параметра MEGNO в задачах динамики ИСЗ
В задачах динамики ИСЗ наиболее интересными с точки зрения исследования хаотичности является геостационарная зона (GEO), которая представляет собой область пространства вокруг экваториальной орбиты со средним радиусом A = 42164 км, шириной 150 км вдоль радиуса орбиты и протяженностью ±15 по широте. Наличие либрационного резонанса в области GEO порождает существование хаотического слоя в окрестности сепаратрисы, разделяющей фазовое пространство на области с различным поведением траекторий. Наложение резонансов в нелинейной динамике считается существенным источником появления хаотичности в поведении систем. Поэтому поискам слабых и вторичных резонансов, посвящено в современной небесной механике значительное число работ (Lemaitre et al, 2009; Delsate, et al, 2008; Breiter, 1999; Breiter, 2001). Исследования хаотичности движения с момента формулирования этой проблемы в конце 60-ых годов прошлого столетия и до начала этого столетия проводились аналитическими методами. С начала этого столетия начали применять симплектические интеграторы и высокоточные численные модели движения. Среди них можно выделить ряд наиболее интересных работ. В работе С. Брейтера и др. (Breiter, 2005) исследование стохастичности зоны GEO сделано в характеристиках MEGNO симплектическим интегратором Йошида (Yoshida, 1993) на интервале 40 лет. В работах С. Волк и др.(Valk et al, 2008a; Valk et al, 2008b; Valk et al, 2009) использованы как численно-аналитический, так и численный подходы, причем в численно-аналитических методах для интегрирования усредненных уравнений использовался метод Рунге-Кутты четвертого порядка, а численное интегрирование проводилось методом Адамса-Мультона-Коулла двенадцатого порядка. В работах Э.Д. Кузнецова и др. (2003, 2007, 2008, 2009) исследование эволюции объектов геостационарной зоны выполнено на основе «Численной модели движения ИСЗ», разработанной в свое время сотрудниками НИИ ПММ ТГУ (Бордовицына и др., 2007b). В качестве основных характеристик хаотичности движения в настоящее время используются ляпуновское время и усредненный параметр MEGNO (Mean Exponential Growth of Nearby Orbit) (Cinсotta et al, 2003), который представляет собой взвешенную по времени интегральную форму ляпуновского характеристического числа (LCN).
Пусть динамическая система определяется системой уравнений (1.1), и ф(ґ) = ф(ґ,х0,ґ0) есть решение этой системы при начальных условиях (t0, х0) , тогда LCN определится как 5 Д) A, = lim-ln \Ъф(0) где Ьф(і), так называемый касательный вектор, который измеряет эволюцию начального бесконечно малого отклонения Ъф (ґ0) = 80 между решением ф(ґ) и очень близкой орбитой. Эта эволюция с точностью до бесконечно малых первого порядка может быть описана вариационным уравнением вида 5 , = —5,() = J(f (ф()))5,(), J(f(Ф())) = — (Ф()), (2.2) где і(ф(ґ)) есть матрица Якоби системы дифференциальных уравнений. Заметим, что параметр LCN может быть представлен также в интегральной форме A,= lim-f )fc, (2.3) /- » f J0 8Д.5) причем 5 = Шф , Ъф = 5ф 5ф / Ъф. Как было отмечено выше, параметр MEGNO Уф (ґ), то есть среднее экспоненциальное расхождение двух близких орбит, представляет собой взвешенную по времени интегральную форму LCN: r {,)=ibw) , (24) а средняя величина Уф (ґ) получается как Уф(і) = - Уф(8 (2.5) Эволюция Уф (t) во времени позволяет выявить различный характер орбит. Так, например, известно, что для квазипериодических (регулярных) орбит Уф (ґ) осциллирует около 2. Более того, как показано в (Cincotta et al, 2003), для квазипериодических орбит Уф (ґ) всегда стремиться к 2, а для устойчивых орбит типа гармонического осциллятора Уф (ґ) = 0. Эти особенности Уф (ґ) будут нами использованы как при тестировании программы вычисления MEGNO, так и при анализе получаемых с помощью этого параметра результатов. В задачах численного моделирования целесообразно заменить (Valk et al, 2009) интегральные соотношения (2.4) и (2.5) дифференциальными уравнениями и интегрировать совместно с уравнениями движения (1.1) и уравнениями в вариациях (2.2) еще два уравнения
MEGNO подход, в отличие от общего вариационного метода Ляпунова, дает более полную динамическую информацию об орбитах и эволюции их касательного вектора (Valk et al, 2009). Этот подход является наиболее эффективным для исследования динамики орбит путем разделения регулярных и стохастических режимов. Рекомендуется совместное интегрирование уравнений (1.1) и (2.8) осуществлять в прямоугольных координатах и скоростях, что позволяет исследовать орбиты с любыми эксцентриситетами и наклонениями. Для интегрирования следует использовать эффективный метод высокого порядка и разрядную сетку ЭВМ, гарантирующие, что ошибки интегрирования не попадут в вычисленные величины MEGNO.
Очевидно, что выбор величины начального вектора 80 может оказывать влияние на характер получаемой эволюции движения, поэтому рекомендуется подбор величины вектора проводить на заведомо регулярных орбитах, для которых Y± (ґ) 2 . Все эти рекомендации были нами учтены при создании нашего алгоритмического и программного комплекса.
Как было отмечено выше, алгоритм вычисления параметров MEGNO строится таким образом, что уравнения движения (1.1) интегрируются совместно с уравнениями в вариациях (2.2) и уравнениями (2.6) или (2.8). В нашем случае 80 есть малая вариация начальных условий (1.2), 5Дґ) есть вариация текущих параметров, а матрица Якоби состоит из четырех блоков J-I 0А I где блок О - нулевой, а матрица А является единичной, а матрица В имеет вид В формуле (2.10) для краткости через j(dV/dx ) обозначена матрица частных производных второго порядка от геопотенциала, а через j(dRM/dx ) и J (dRs /их ) соответственно матрицы частных производных второго порядка от возмущающих функций Луны и Солнца, j(dPk/dx. ) - матрица частных производных второго порядка от радиационного давления.
Вычисление вторых частных производных от потенциала притяжения Земли может производиться двумя способами. В основу этих способов вычисления положено представление потенциала Земли и его производных в виде где функции V и их частные производные вычисляются в прямоугольной системе координат жестко связанной с Землей; GHM@ — соответственно гравитационная постоянная и масса Земли; Яф — ее средний экваториальный радиус; Сиот и Snm — нормированные численные коэффициенты, характеризующие структуру гравитационного поля Земли; Vnm — нормированные шаровые функции; i = J-\; символ Real определяет действительную часть соответствующих функций. Связь между нормированными и ненормированными функциями и постоянными определяется формулами
Особенности реализации алгоритма в среде параллельных вычислений
Прежде чем приступить к реализации комплекса программ для решения обратных задач динамики ИСЗ нужно разработать алгоритм работы программного комплекса. Алгоритм должен обеспечивать эффективное использование параллельной вычислительной системы. Для равномерности и сбалансированности загрузки ядер кластера было рассмотрено три подхода к реализации программного комплекса, отличающиеся между собой количеством используемых ядер (Чувашов, 2013 c).
Программный комплекс, использующий два ядра кластера. В этом случае программный комплекс разделяет интервал наблюдения на две части: наблюдения, сделанные до эпохи улучшения орбиты, и наблюдения - после эпохи улучшения. При этом первое ядро интегрирует уравнения движения с отрицательным шагом интегрирования, а второе с положительным шагом (рисунок 3.1).
Программный комплекс, использующий три ядра кластера. В этом случае два ядра, как и в первом случае, используются для интегрирования уравнения движения. А третье ядро выполняет второстепенные функции: чтение из файлов начальных данных, запись в файл результатов представления наблюдений. Этот программный комплекс является промежуточным этапом между программными комплексами, использующими два или четыре ядра кластера (рисунок 3.2). Нужно заметить, что работа третьего ядра тут происходит последовательно относительно двух других ядер. а) два ядра б) три ядра
Программный комплекс, использующий четыре ядра кластера. В этом случае два нечетных ядра, как и в первом случае, используются для интегрирования уравнения движения. А два четных ядра используются для определения отдельных возмущений, входящих в правые части дифференциальных уравнений движения, на каждом шаге интегрирования (рисунок 3.3).
Кроме реализации этих трех подходов, процесс разработки параллельного алгоритма был разбит на четыре этапа.
Информационный граф программного комплекса, использующий четыре ядра кластера. Этап I, II –обмен данных в функции правых частей; III, IV – обмен входными и выходными данных Декомпиляция. На данном этапе задача разработки программного комплекса анализировалась, оценивалась возможность ее распараллеливания. Программный комплекс, адаптированный на кластере «СКИФ Cyberia», должен был учитывать слабые возмущения, влияющие на движение ИСЗ (раздел 4), обладать увеличенным быстродействием программного комплекса и повышенной точностью определения ковариационной матрицы.
Данные, приведенные в таблице 3.1, показывают, как зависит время работы программного комплекса от количества задействованных ядер.
Сравнение времени работы программного комплекса в зависимости от количества задействованных ядер
Количество задействованных ядер Время выполнения задачи улучшения орбиты Лагеоса на семисуточном интервале, сек
Проектирование коммуникаций между задачами. На этом этапе определялись коммуникации, необходимые для пересылки исходных данных, промежуточных результатов выполнения подзадач, а также коммуникации, необходимые для управления работой подзадач. Результаты этого анализа представлены на рисунках 3.1 – 3.3.
Укрупнение. Для повышения эффективности алгоритма и снижения трудоемкости разработки для третьего подхода было выполнено укрупление, которое позволило объединить логику работы четных ядер в одну процедуру. Это повысило масштабируемость программного комплекса. Для первого и второго похода укрупнение не разрабатывалось.
Планирование вычислений. На этом заключительном этапе производилось распределение функций между процессорами и формирование набора задач, в который включаются вновь созданные задачи. Стратегия размещения задач на процессорах представляет собой компромисс между требованием максимальной независимости выполняющихся задач и глобальным учетом состояния вычисления. Для третьего подхода нами была выбрана стратегия «хозяин/ работник» (Немнюгин, Стесик, 2002) (рисунок 3.4).
Нужно заметить, что программный комплекс, использующий четыре ядра, позволяет представлять наблюдения n объектов. В этом случае для каждого объекта задается свой входной файл, и количество ядер становится равно 4n . Информационный граф в этом случае (рисунок 3.5) представляет собой n раз скопированный узел, изображенный на рисунке 3.4.
Информационный граф программного комплекса, который позволяет улучшать орбиты n объектов. Этап I – пересылка начальных данных ядрам, которые управляют локальной группой ядер В таком программном комплексе нужно строго определить функции каждого ядра. Первое ядро является главным ядром для всего программного комплекса и имеет следующие функции: – функция инициализации – пересылка начальных данных ядрам, которые управляют своими узлами; – функция управления ядрами, то есть отслеживание работы ядер, управляющих своими узлами, обработка ошибок, приходящих от этих узлов и прерывание работы всего программного комплекса, из-за ошибок, возникших в программном комплексе; – функция управления узлом, так как первое ядро является главным ядром в узле, то на него ложатся функции по управлению своим узлом; – интегрирование дифференциальных уравнений движения объектов. Функции ядер, управляющих своими узлами, следующие: – функция по обработке пакетов данных, полученных на узле; – функция управления – отслеживание работы ядер, входящих в узел, обработка ошибок, которые возникли в узле, и пересылка этих ошибок главному ядру; – интегрирование дифференциальных уравнений движения объектов.
Эффективность разработанного программного комплекса может быть определена законом Амдала (Немнюгин, Стесик, 2002), который показывает прирост производительности программы, выполняющейся на параллельных процессорах (таблица 3.2). Таблица 3.2 — Эффективность программных комплексов
Исследование эффективности разработанного комплекса программ для решения обратных задач динамики систем ИСЗ и КА проводилось по двум направлениям. Оценивалось повышение быстродействия за счет распараллеливания процесса вычислений и точность вычислений при использовании 128-ми битной разрядной сетки.
Программный комплекс для решения обратных задач тестировался на лазерных наблюдениях геодинамических ИСЗ Лагеос и ИСЗ Эталона. Обработка наблюдений для ИСЗ Лагеоса проводилась на 7-суточной орбитальной дуге, что соответствует 60 оборотам спутника, а для ИСЗ Эталон - на 30-суточной орбитальной дуге. В таблице 3.5 и 3.8 приведены среднеквадратичные ошибки единицы веса для улучшения орбиты ИСЗ Лагеос и ИСЗ Эталон в зависимости от количества отбраковок, которая определялась по правилу «3 j».
Результаты исследования орбитального движения фрагмента космического мусора на GEO по данным четырёхлетних оптических наблюдений
На рисунке 4.3 помимо оценки влияния составляющих прилива, приведены две численные ошибки интегрирования, соответствующие 64 битной и 128 битной разрядным сеткам. Как видно из графиков, влияние значительной части составляющих прилива больше чем численные ошибки интегрирования. Однако оценку влияния возмущения от деформации коэффициента C2,2 можно учесть только на сетке соответствующей 32 десятичным разрядам.
В процессе выполнения программы TOPEX/Poseidon было создано несколько моделей воздействия океанических приливов на коэффициенты разложения гравитационного потенциала Земли. Краткое описание первой аналитической модели TEG-2B, а также всех моделей CSR дано в разделе 1.4.2. Здесь мы приведем оценки вклада в движение спутника Лагеос влияния океанических приливов, представленных вышеперечисленными моделями, для разных интервалов времени (таблицы 4.2 - 4.4). В Таблицах также представлена разность влияния каждой океанической модели относительно друг друга.
Таким образом, исследование показало, что использование 128 битной разрядной сетки позволяет учитывать такие слабые возмущения, которые невозможно отследить на 64 битной разрядной сетке, а это в свою очередь расширяет наши возможности в решении прямых и обратных задач движения ИСЗ, а также в прогнозировании опасных столкновений фрагментов космического мусора с функционирующими и вновь выводимыми космическими аппаратами. 4.3 Оценка влияния преломления солнечных лучей в земной атмосфере и сжатия Земли на величину возмущений от светового давления
Используя формулы (1.22) – (1.28), можно исследовать поведение функции тени на примере спутника Лагеос (рисунок 4.4), где SLCYLIN – «цилиндрическая» модель, SLCONE – «конусная» модель, SLATM – модель, учитывающая атмосферу и сжатие Земли. Видно, что при
использовании полной модели поток радиационного давления, приходящего с Солнца, начинает уменьшаться раньше по причине перекрытия диска Солнца земной атмосферой, и затмение продолжается дольше из-за двойного преломления солнечных лучей в атмосфере. При таком подходе происходит в два раза больше обращений к правым частям дифференциальных уравнений в численной модели (рисунок 4.5), что приводит к падению быстродействия. Связано это с тем, что интегратор, используемый в численной модели, пересчитывает шаг интегрирования при быстром изменении правых частей. Так как процесс длится не продолжительное время - примерно 30 секунд, то и падение быстродействие наблюдается только на этом интервале времени.
В таблице 4.4 показан вклад в движение каждой модели светового давления. Эксперимент проводился на 60 оборотах ИСЗ Лагеоса, на каждом обороте объект пересекал тень от Земли. В эксперименте проводилась оценка вклада модели светового давления в движение спутника. Видно, что модель, учитывающая атмосферу и сжатие Земли, дает значительный вклад в движение объекта в отличие от двух первых моделей. Это отличие объясняется тем, что затмение для последней модели длится в полтора раза дольше, чем для «конусной» модели (рисунок 4.4). Поэтому интегратор успевает почувствовать быстрое изменение правых частей дифференциальных уравнений движения искусственного спутника Земли (рисунок 4.5) и уменьшить шаг интегрирования, что приводит к увеличению точности, но к уменьшению быстродействия.
Точно такой же вклад вносят модели светового давления, когда мы пытаемся с помощью численной модели движения определить координаты спутника Лагеос на интервале времени 60 оборотов (таблица 4.4).
Приведенные результаты показывают, что среди рассмотренных факторов более существенным следует считать учет конусности тени Земли, который дает повышение точности моделирования движения на 16 %, учет преломления светового луча в атмосфере дает 12 %, что в задачах высокоточного прогнозирования также может быть существенно.
Как известно, эффект Пойнтинга—Робертсона заставляет пылевые частицы терять момент количества движения и объект по спирали падает на центральное притягивающее тело. Этот процесс можно наблюдать сегодня чуть ниже геостационарной зоны, куда микрочастицы эволюционировали с геостационарной зоны (Valk et al, 2009).
Исследование влияния эффекта Пойнтинга—Робертсона на динамическую эволюцию пылевых частиц, геостационарной зоны проводилось с помощью разработанного ранее программного комплекса для одновременного численного интегрирования систем уравнений n объектов (Бордовицына и др, 2009). При моделировании движения геостационарных объектов учитывались следующие возмущающие факторы: несферичность Земли до гармоник 10-го порядка и степени, приливные деформации, лунно-солнечные возмущения, эффект Пойнтинга—Робертсона и релятивистские эффекты.
В нашем исследовании мы рассматривали металлические пылевые частицы, имеющие размер с 1 нанометра до 10 микрометров. Зная размер и считая, что частицы имеют форму шара, мы легко определили массу и площадь миделевого сечения частицы. Нами было выбрано семь разных типов движения частиц: устойчивое движение с малой амплитудой либрации относительно точки 75 градусов или точки 255 градусов, движение с большой амплитудой либрации относительно точки 75 градусов или точки 255 градусов, либрационное движение относительно двух устойчивых точек 75 градусов и 255 градусов, круговое движение и неустойчивое движение.
На рисунке 4.6 показана эволюция пылевых частиц размером 1 нанометр в зависимости от типа движения. Так как эволюция частиц мало чем отличается друг от друга, мы увеличили график на небольшом временном интервале, чтобы показать близкие траектории движения этих частиц. Небольшие колебания связаны с периодом вращения объектов. Как видно из графика, характер движения не влияет на динамическую эволюцию движения частиц, а влияние эффекта Пойнтинга—Робертсона приводит к стремительной потере момента количества движения и падению частиц на Землю за 13 лет.
Так как поведение орбитальных характеристик для разных физических параметров объектов будет отличаться друг от друга, то время падения частиц на Землю будет увеличиваться с увеличением массы и, следовательно, размера частиц (рисунок 4.7). Одновременное уменьшение значений большой полуоси и эксцентриситета говорит о том, что объекты снижаются по спирали.