Санкт-Петербургский государственный университет
Механика и математическое моделирование
Физическая механика
Кремнев Иван Александрович
Математическое моделирование процесса
образования облака из техногенных
микрочастиц, непрерывно инжектируемых
точечным источником, движущимся в
околоземном космическом пространстве по
заданной орбите
Бакалаврская работа
Научный руководитель:
д. ф.-м. н., профессор Е. К. Колесников
Рецензент:
к. ф.-м. н. С. В. Чернов
Санкт-Петербург
2016
SAINT-PETERSBURG STATE UNIVERSITY
Mechanics and mathematical modelling
Physical mechanics
Ivan Kremnev
Mathematical modelling of technogenic
microparticle cloud formation process with
continuous injection from a point source
orbiting in near-earth space along selected
trajectory
Bachelor’s Thesis
Scientific supervisor:
professor Evgeniy Kolesnikov
Reviewer:
PhD Sergey Chernov
Saint-Petersburg
2016
Оглавление
Введение
4
1. Физическая постановка задачи
5
2. Силы, определяющие динамику источника и частиц облака
7
3. Математическая модель
9
4. Особенности динамики облака МЧ на различных орби10
тах
5. Анализ облака МЧ на орбите «Молния»
16
Заключение
19
Список литературы
21
3
Введение
Ввиду постоянного роста числа пассивных орбитальных объектов
техногенного происхождения (космического мусора) в околоземном космическом пространстве (ОКП) проблема наблюдения и предсказания
эволюции облака космического мусора в ОКП приобретает всё большую значимость.
Можно выделить несколько направлений исследований в данной области: изучение динамики крупных объектов космического мусора, таких как вышедшие из строя КА и их обломки; исследование динамики
мелких частиц размера порядка нескольких миллиметров; изучение динамики микрочастиц (МЧ) размера менее миллиметра.
В частности, важное место занимает вопрос о динамике техногенных МЧ размера от 1 µm до 400 µm. Такими частицами являются, в
частности, продукты деградации материала поверхности КА. Наблюдения показывают, что на некоторых околоземных орбитах могут существовать облака МЧ, называемые астрозолями, которые представляют
угрозу для функционирования КА на данных орбитах [7]. В области
низких высот (до 2000 km) максимальная плотность потока частиц в
астрозольных облаках может достигать значений порядка 12 m−2 s−1
при средней фоновой плотности потока порядка 1.7 × 10−3 m2 s−1 . Механизм появления таких облаков до сих пор не изучен, однако предполагается, что одним из возможных источников техногенного загрязнения
ОКП в области низких орбит могут служить объекты, движущиеся по
высоким эллиптическим орбитам с низким перигеем типа орбиты спутника «Молния» [4].
В данной работе ставится задача построения математической модели процесса загрязнения области низких орбит продуктами деградации поверхности КА в условиях их непрерывной инжекции. Основное
внимание уделено исследованию особенности динамики облака техногенных МЧ, инжектируемых точечным источником, движущимся по
орбите типа «Молния».
4
1. Физическая постановка задачи
Исследуется эволюция облака техногенных МЧ, непрерывно инжектируемых на выбранных околоземных орбитах точечным источником.
Рассматриваются орбиты:
• Высокие эллиптические орбиты (ВЭО) типа «Молния»
• Низкие круговые орбиты
• Высокие круговые орбиты
Параметры орбит представлены в Таблице 1.
Орбита
Наклонение Высота перигея Высота апогея
«Молния»
62.8°
40 000 km
500 km
Орбита спутника
28.5°
473 km
483 km
LDEF
Геостационарная
0°
35 786 km
35 786 km
орбита
Таблица 1: Элементы избранных орбит
Движение источника МЧ определяется силой притяжения к Земле и
описывается классическим уравнением движения материальной точки
с радиус-вектором ⃗rs и вектором скорости ⃗vs в геоцентрической инерциальной системе координат (ECI):
m
d⃗
vs
= F⃗g + F⃗gdist
dt
(1)
В формуле (1) F⃗g — сила тяготения в центральном гравитационном
поле Земли, F⃗gdist — возмущение силы тяготения, обусловленное отклонением гравитационного поля Земли от центрального.
Облако МЧ описывается отдельными частицами. Динамика каждой МЧ определяется наборами индивидуальных параметров частицы
и сил, действующих на неё со стороны окружающей среды, и также описывается законом Ньютона. Рассмотрена следующая физическая мо-
5
дель эволюции шарообразной частицы радиуса R и массы m с радиусвектором ⃗r и вектором скорости ⃗v в системе координат ECI:
d⃗v
= F⃗g + F⃗gdist + F⃗rad + F⃗drag
(2)
dt
Кроме упомянутых сил, здесь представлены: F⃗rad — сила солнечного
давления, F⃗drag — сила сопротивления электронейтральной компоненты
m
атмосферы Земли.
Сила Лоренца в данной постановке не принимается в расчёт, так
как для МЧ размера более 0.1 µm ее относительный вклад в динамику
МЧ слишком мал [5]. Гравитационные возмущения, связанные с притяжением Солнца и Луны также не учитывались при решении задачи на
временах до месяца, будучи малыми относительно других возмущений
(порядка 6 mm s−2 и 0.03 mm s−2 соответственно).
6
2. Силы, определяющие динамику источника и частиц облака
Сила F⃗g центрального гравитационного поля Земли определяется
формулой
F⃗g = −m∇V0
V0 = −µ/r
(3)
где V0 — потенциал центрального гравитационного поля Земли,
µ = GME — гравитационный параметр Земли, равный 3.986 × 1014 m3 s−2 .
Возмущение центрального гравитационного поля Земли F⃗ dist опиg
сывается рядом Гаусса в соответствии с формулой
Vdist
F⃗gdist = −m∇Vdist
∞
n
µ ∑ ∑ ( r0 )n
=−
(cnm cos(mϕ) + dnm sin(mϕ))Pnm (sin ψ)
r n=1 m=0 r
Здесь постоянная r0 = 6.363 553 × 106 m, Pnm — присоединенные многочлены Лежандра, а коэффициенты cnm , dnm определяются эмпирически. Аргументами являются географическая долгота ϕ и географическая широта ψ [5]. Основной вклад в возмущение гравитационного поля
Земли вносит вторая гармоника, отвечающая за полярное сжатие Земли и вычисляемая по формуле
V20
F⃗gdist = −m∇V20
2
= µc2r203r0 (3 sin2 ψ − 1)
(4)
Силу солнечного давления F⃗rad , природа которой заключается в передаче МЧ импульса фотонов солнечного света при их поглощении и
рассеянии веществом МЧ, можно выразить соотношением
Qpr N πR2
⃗
⃗ι
(5)
Frad =
c
Коэффициент Qpr является материалозависимой эффективностью
солнечного давления на МЧ, усредненной по всему спектру солнечного
7
излучения, и вычисляется согласно формулам теории Ми [5]. Значения коэффициента для различных материалов и радиусов МЧ взяты
из материалов работы [6]. Величина N отвечает за плотность потока
солнечной энергии на орбите Земли и равна 1366 W m−2 , R — радиус сферической микрочастицы, c = 2.997 924 58 × 108 m s−1 — скорость
света в вакууме. Вектор ⃗ι есть единичный вектор, направленный от
Солнца к Земле. Годичным движением Солнца на рассмотренных временах порядка десятка дней можно пренебречь, положив этот вектор
постоянным. Следует отметить, что эффект Пойнтинга-Робертсона не
принимается в расчёт в виду малости отношения vc .
Последним среди учтённых возмущений является сила сопротивления движению МЧ со стороны атмосферы Земли. Эффект данной силы
сказывается на высотах менее 1000 km и вычисляется по формуле:
2
πR
F⃗drag = −cx
ρa v⃗v
(6)
2
Здесь коэффициент лобового сопротивления шара cx = 2 соответствует модели неупругого столкновения молекул фонового газа с МЧ,
а ρa — общая высотная плотность атмосферы. Значения ρa рассчитаны
на основе модели атмосферы NRLMSISE-00 с усреднением по сезонным
и суточно-широтным колебаниям и значениям параметров в периоды
низкого, среднего и высокого уровней солнечной и геомагнитной активностей для выбранных высот [3]. Между выбранными значениями высот плотность атмосферы аппроксимируется кусочно-экспоненциально.
8
3. Математическая модель
Уравнение движения точечного источника МЧ (1) и N уравнений
движения N частиц вида (2) могут быть записаны в виде N + 1 системы из 6-ти обыкновенных дифференциальных уравнений 1-го порядка.
Численное решение задачи Коши для каждой системы использует метод Рунге-Кутты 4-го порядка точности с постоянным шагом h = 1 min.
Погрешность решения на всём времени существования МЧ оценивается
по правилу Рунге и не превышает 40 km.
Для моделирования непрерывной инжекции МЧ через равные промежутки времени τ = 10 min создается 1 частица выбранного радиуса
со случайной начальной скоростью ⃗v0 относительно источника, распределение которой равномерно по величине от 0 до 10 cm s−1 и по направлению в пространстве относительно источника [1], [2]. Выбор определенного размера частицы обусловлен рассмотрением динамики соответствующей фракции облака МЧ. Эта частица добавляется к уже существующему массиву частиц, тем самым добавляя одну систему дифференциальных уравнений типа (2) к уже имеющимся. Падающие или
покидающие ОКП частицы удаляются из облака, тем самым сокращая
количество имеющихся систем типа (2) на одну.
Таким образом, число решаемых систем уравнений меняется каждый шаг численного метода. Для уравнения движения источника МЧ (1)
начальными данными служат параметры запуска КА (⃗rs0 , ⃗vs0 ), для уравнений движения МЧ (2) — текущая координата источника и случайно
заданная скорость частицы (⃗rs , ⃗vs + ⃗v0 ).
9
4. Особенности динамики облака МЧ на различных орбитах
-5
-5
-5.5
-5.5
1000 km
1000 km
Запуск симуляции осуществляется посредством задания начальных
значений положения и скорости источника в географических координатах: высоты, широты и долготы точки запуска, величины и углов
азимута и подъема скорости.
На высотах менее 1000 km особенности динамики МЧ любого размера определяет сила атмосферного сопротивления F⃗drag . Для низких
круговых орбит типа орбиты спутника LDEF это сказывается в быстром торможении инжектируемых МЧ и их падении на Землю в течение
2-14 часов с момента инжекции. Устойчивого плотного облака МЧ не
образуется (Рис. 1). Здесь и далее при изображении облака МЧ синем
цветом отражены МЧ, красным — источник, бирюзовым — контур Земли; масштаб указан рядом с осями.
-6
-6
-6.5
-6.5
-7
-7
-3
-2.5
-2
-1.5
1000 km
-1
-0.5
0
-3
(a) Частицы радиуса 2 µm
-2.5
-2
-1.5
1000 km
-1
-0.5
0
(b) Частицы радиуса 10 µm
Рис. 1: Облако МЧ спустя 100 часов с момента запуска источника на
низкой круговой орбите
На высоких круговых орбитах для МЧ различных размеров основным возмущением является сила солнечного давления F⃗rad . Хотя в орбитальном масштабе облако держится компактно в течение длительного времени, за счёт ускорения солнечным давлением оно получает
своеобразную зигзагообразную форму (Рис. 2,3), которая сжимается и
10
35
35
30
30
25
25
1000 km
1000 km
расширяется в поперечном направлении при движении по орбите. Подобная форма объясняется действием солнечного давления на МЧ, которое с течением времени спускает частицы, инжектированные на высоких круговых орбитах, на ВЭО с низким перигеем. Так как каждая
частица инжектируется источником в свой момент времени и в определённом месте материнской орбиты, сход частиц на ВЭО выглядит как
непрерывное от частицы к частице изменение параметров их орбиты.
Такой источник также не создаёт плотных потоков частиц в низкоорбитальной области.
20
15
20
15
10
10
5
5
0
0
-10
0
10
20
1000 km
30
40
-10
(a) Частицы радиуса 2 µm,
крупный масштаб
0
10
1000 km
20
30
(b) Частицы радиуса 10 µm,
крупный масштаб
38
37
36
36.5
1000 km
1000 km
34
32
30
28
36
35.5
26
35
24
10
15
20
1000 km
25
18.5
(c) Частицы радиуса 2 µm,
вблизи источника
19
19.5
20
20.5
1000 km
21
21.5
(d) Частицы радиуса 10 µm,
вблизи источника
Рис. 2: Облако МЧ спустя 100 часов с момента запуска источника на
геостационарной орбите
11
0
0
-5
1000 km
1000 km
-10
-20
-10
-15
-20
-30
-25
-30
-40
-10
0
10
20
1000 km
30
40
50
0
(a) Частицы радиуса 2 µm,
крупный масштаб
10
20
1000 km
30
40
(b) Частицы радиуса 10 µm,
крупный масштаб
-22
-25
1000 km
1000 km
-24
-30
-35
-26
-28
-40
-30
-45
15
20
25
30
1000 km
35
40
28
(c) Частицы радиуса 2 µm,
вблизи источника
30
32
34
1000 km
36
38
40
(d) Частицы радиуса 10 µm,
вблизи источника
Рис. 3: Облако МЧ спустя 330 часа с момента запуска источника на
геостационарной орбите
На высоких эллиптических орбитах значительную роль играет возмущение силы притяжения F⃗gdist . Именно оно определяет особенности
динамики облака для частиц радиуса 10 µm и больше. Все частицы данных размеров остаются вблизи материнской орбиты, которая прецессирует с течением времени под действием F⃗gdist . Наименее подвержены
прецессии ВЭО с наклонением около 63°, что соответствует наклонению орбиты «Молния», в то время как на орбитах с малым наклонением при участии солнечного давления происходит расслоение облака
в масштабах нескольких десятков километров (Рис. 4).
12
6
60
4
2
10 km
1000 km
40
20
0
-2
0
-4
-20
-6
-10
0
10
20
30
-10
0
1000 km
10
20
30
1000 km
(a) Плоскость материнской орбиты
(b) Перпендикулярная плоскость
Рис. 4: Облако МЧ радиуса 10 µm спустя 83 часа с момента запуска
источника на ВЭО с малым наклонением
На орбитах «Молния», тем не менее, большое значение вновь приобретает сила солнечного давления. За счёт неё на начальных этапах
частицы идут вблизи материнской орбиты с опережением источника
(Рис. 5). Спустя некоторое время с момента запуска источника более
старые частицы начинают уходить с материнской орбиты, уменьшая
высоту апогея своей орбиты. Со временем облако заполняет всю орбиту (Рис. 6, 7, 8, 9).
-11
10
-12
5
1000 km
1000 km
-13
-14
0
-5
-10
-15
-15
-16
-20
-17
-38
-36
-34
1000 km
-32
-30
(a) Вблизи источника
-20
-10
1000 km
0
(b) Крупный масштаб
Рис. 5: Облако МЧ радиуса 10 µm спустя 80 часов с момента запуска
источника на орбите «Молния»
13
20
15
15
10
10
5
5
1000 km
1000 km
20
0
0
-5
-5
-10
-10
-15
-15
-40
-30
-20
1000 km
-10
0
-40
(a) Спустя 33 часа
-30
-20
1000 km
-10
0
(b) Спустя 330 часов
20
20
15
15
10
10
5
5
1000 km
1000 km
Рис. 6: Облако частиц радиуса 2 µm от источника на орбите «Молния»
0
0
-5
-5
-10
-10
-15
-15
-40
-30
-20
1000 km
-10
0
-40
(a) Спустя 78 часов
-30
-20
1000 km
-10
0
(b) Спустя 330 часов
20
20
15
15
10
10
5
5
1000 km
1000 km
Рис. 7: Облако частиц радиуса 5 µm от источника на орбите «Молния»
0
0
-5
-5
-10
-10
-15
-15
-40
-30
-20
1000 km
-10
0
-40
(a) Спустя 108 часов
-30
-20
1000 km
-10
0
(b) Спустя 330 часов
Рис. 8: Облако частиц радиуса 7 µm от источника на орбите «Молния»
14
15
10
10
5
5
1000 km
1000 km
15
0
0
-5
-5
-10
-10
-15
-15
-40
-30
-20
-10
1000 km
0
-40
(a) Спустя 150 часов
-30
-20
-10
1000 km
0
(b) Спустя 330 часов
Рис. 9: Облако МЧ радиуса 10 µm от источника на орбите «Молния»
15
5. Анализ облака МЧ на орбите «Молния»
Подробные количественные исследования проведены для ВЭО «Молния», представляющей наибольший интерес в вопросе происхождения
облака техногенных МЧ в области низких околоземных орбит. Исследована динамика отдельных фракций частиц размеров от 1 до 50 µm в
течение 19 суток с момента запуска.
При симуляции образования облака скорость инжекции считалась
постоянной для частиц всех размеров: τ −1 = 0.1 min−1 . Чтобы корректно интерпретировать полученные значения потока МЧ, необходимо знать действительную скорость деградации материала поверхности
КА и распределение по размерам инжектируемых частиц. Распределение было принято в соответствии с данными наблюдений оптического
окружения ходе миссии КА MSX [1], орбита которого была круговой
с высотой 900 km. По имеющимся данным, наибольшую часть облака
составляют частицы радиуса 10 µm и 20 µm. Частота инжекции МЧ
была рассчитана исходя из значения массовой потери краски КА Space
Shuttle в ходе эксперимента EOIM-3 на низкой орбите, составившей 700
gm−2 y−1 , что является оценкой сверху массовой потери материала покрашенной поверхности КА, движущегося по высокой эллиптической
орбите.
Как было продемонстрировано на Рис. 6,7,8,9, все фракции облака имеют однотипную динамику, отличающуюся скоростью схода МЧ
с материнской орбиты и темпом уменьшения высоты апогея их орбиты. Перигеи орбит всех отдельных потоков облака почти совпадают с
перигеем материнской орбиты, что позволяет делать заключение о наивысшей плотности потока частиц вблизи него. Поэтому область перигея
материнской орбиты является наиболее вероятной областью нахождения астрозолей.
Проявление эффекта солнечного давления падает с ростом размера частиц, что сказывается в уменьшении скорости их схода с материнской орбиты, однако его наличие определяет локальную структуру
облака. В силу непрерывности инжекции МЧ источником появляются
16
группы МЧ, на которые солнечное давление успевает оказать примерно
одинаковое воздействие на начальных этапах их существования. Такие
группы выглядят как непрерывное образование в пространстве (Рис.
10), которое обладает большей плотностью для частиц большего размера. При прохождении источника через область низких орбит под действием атмосферного сопротивления группа отрезается от источника и
перестает формироваться. Тем не менее, она продолжает двигаться по
орбите в компактном виде в течение некоторого времени, о чем свидетельствует наблюдение подобных групп МЧ на некотором отдалении от
перигея материнской орбиты. Астрозоли в низкоорбитальной области
могут быть представлены, в частности, именно такими группами МЧ.
0
-2
-4
-6
10 km
-8
-10
-12
-14
-16
-18
-20
-42
-40
-38
-36
-34
-32
-30
-28
-26
-24
-22
1000 km
Рис. 10: Группы МЧ радиуса 10 µm. Плоскость, перпендикулярная орбите
При прохождении такой группы МЧ через район перигея происходит ее уплотнение, так как все частицы облака, независимо от своей
текущей орбиты, имеют общий район перигея. При данных обстоятельствах сечение области пространства, занимаемой такой группой, может
иметь площадь 108 m2 при поперечных размерах порядка около 900 m,
а плотность потока МЧ в ней достигает 9 m−2 s−1 (Рис. 11). Протяженность подобной группы МЧ составляет 530 km, и время ее прохода
17
через район перигея приближенно равняется 13 секундам.
×105
×104
4
-0.8
2
m
m
-1
0
-1.2
-2
-1.4
-4
-1.6
5
5.5
6
m
6.5
6.76
×106
6.78
6.8
m
(a) Все частицы
6.82
×106
(b) Две близкие группы
Рис. 11: МЧ радиуса 10 µm, проходящие район перигея в момент времени 452 часа 20 минут с момента запуска источника
Так как скорости падения на Землю частиц каждой фракции облака
отличны, стационарное состояние облака устанавливается по истечении
среднего времени жизни на орбите наиболее крупных частиц облака, в
данном случае 50 µm. Это соответствует 450 часам после запуска источника. По истечении этого времени можно делать выводы о средней
плотности потока частиц в районе перигея орбиты «Молния». На Рис.
12 изображен такой средний поток. На занимаемой им площади около
6.6 × 109 m2 его величина составляет 1.8 × 10−3 m−2 s−1 . Здесь же видна одна плотная группа МЧ, движущаяся на отдалении от основного
потока, с плотностью 8 m−2 s−1 (помечена розовым цветом).
Полученное значение среднего фонового потока МЧ в районе перигея ВЭО «Молния» совпадает с регистрированными в ходе миссии КА
LDEF [4]. Значения потока плотной группы МЧ отличается на порядок
от значения потока частиц в астрозоле. Исходя из этих соображений,
можно сделать вывод, что деградация поверхностных материалов КА,
движущихся на ВЭО типа «Молния» может быть источником астрозолей в низкоорбитальной области.
18
×105
1
m
0
-1
-2
-3
-4
5.5
6
m
6.5
×106
Рис. 12: МЧ, проходящие район перигея в момент времени 479 часов 20
минут с момента запуска источника
Заключение
На основе построенной математической модели динамики облака
техногенных МЧ в условиях их непрерывной инжекции точечным источником, движущимся по выбранной орбите в ходе численных экспериментов были получены следующие результаты:
• Низкие круговые орбиты не являются источником астрозолей на
низких орбитах. Среднее время жизни частиц составляет от 2 до
14 часов для МЧ различных фракций вследствие быстрого торможения частиц атмосферой.
• Высокие круговые орбиты также не являются источником астрозолей на низких орбитах. Образующиеся облака рассеиваются под
действием солнечного давления, которое также медленно сводит
частицы на ВЭО с низким перигеем.
• Высокие эллиптические орбиты с низким перигеем представляют потенциальный источник потоков МЧ в низкоорбитальной области со средними значениями порядка 10−3 m−2 s−1 и пиковыми
19
порядка 9 m−2 s−1 . Вероятно, появление астрозолей с плотностью
потока частиц около 12 m−2 s−1 может быть объяснено деградацией поверхности материала КА, движущегося по данным орбитам.
В качестве дальнейших исследований следует разобрать вопрос динамики частиц земного происхождения, выбрасываемых с поверхности
КА при выходе на ВЭО. Также следует подробнее изучить процесс деградации других материалов поверхности КА в условиях преобладания
безатмосферной среды.
20
Список литературы
[1] Green B.D, Galica G.E. et al. Optical environment surrounding the
MSX spacecraft // Proceedings of the 7th International Symposium on
Materials in Space Environment, Touluse, France. –– 1997.
[2] Klinkrad H. Space Debris. Models and Risk Analysis / European Space
Agency. –– Springer-Praxis, 2006. –– ISBN: 3-540-25448-X.
[3] NRLMSISE-00 empirical model of the atmosphere: Statistical
comparison and scientific issues / J.M. Picone, A.E. Hedin, D.P. Drop,
A.C. Atkin // Geophys. Res. –– 2002. –– Vol. 107, no. A12. –– P. SIA
15–1 — SIA 15–16.
[4] Singer S.F., Mulholland J.D. et al. LDEF Interplanetary Dust
Experiment: Techniques for the Identification and Study of Long-Lived
Orbital Debris Clouds. –– 1991. –– IAF PAPER 91-285.
[5] Колесников Е.К. Динамические модели процессов распространения
потоков заряженных частиц в космической плазме : Диссертация на
соискание степени доктора ф.-м. наук / Е.К. Колесников ; СПбГУ. ––
СПб, 1998.
[6] Исследование астрозолей в околоземном космическом пространстве
с использованием результатов бортовых измерений и математического моделирования. Анализ воздействия потоков астрозолей на
элементы конструкции космических аппаратов : Отчет : МНТЦ
3412 / Научно-исследовательский институт математики и механики им. академика В.И. Смирнова Санкт-Петербургского государственного университета ; исполн.: Е.К. Колесников и др. –– 198504,
Россия, Санкт-Петербург, Петродворец, Университетский пр., 28 :
2010.
[7] Колесников Е.К., Чернов С.В. О возможности длительного орбитального существования субмикронных частиц, инжектируемых в
21
околоземное космическое пространство на вытянутых эллиптических орбитах с низким перигеем // Космические исследования. ––
2013. –– Т. 51 №4. –– С. 287–293.
22
Отзывы:
Авторизуйтесь, чтобы оставить отзыв