Интегрирование уравнений движения малых тел Солнечной системы методом оскулирующих элементов

Тип работы:
Реферат
Предмет:
Физико-математические науки


Узнать стоимость

Детальная информация о работе

Выдержка из работы

Небесная механика и астрономия
УДК 523. 642
ИНТЕГРИРОВАНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ МАЛЫХ ТЕЛ СОЛНЕЧНОЙ СИСТЕМЫ МЕТОДОМ ОСКУЛИРУЮЩИХ ЭЛЕМЕНТОВ
А. Ф. Заусаев, Д. А. Заусаев
Самарский государственный технический университет,
443 100, Самара, ул. Молодогвардейская, 244.
E-mails: Zausaev_AFamail. ru
Предложен алгоритм метода оскулирующих элементов для численного интегрирования уравнений движения малых тел Солнечной системы.
Ключевые слова: интерполяция, численное интегрирование, дифференциальные уравнения движения, метод оскулирующих элементов.
Дифференциальные уравнения движения. В связи с решением проблемы астероидной опасности исследование эволюции орбит малых тел Солнечной системы (астероидов, комет, крупных фрагментов в метеорных потоках) является актуальной задачей.
Большое значение имеет выбор математической модели, описывающей движение небесных тел. Учет лишь гравитационных взаимодействий, который заложен в стандартных уравнениях движения задачи n тел, оказывается недостаточным. Это приводит к более сложной форме дифференциальных уравнений движения, где наряду с гравитационными эффектами учитываются релятивистские эффекты, фигура планет и Солнца, реактивные силы и т. д.
В работе [1] предложен интерполяционный метод оскулирующих элементов, используемый для получения координат и элементов орбит больших планет и малых тел Солнечной системы. Этот метод предназначен для работы с банком данных координат и элементов орбит малых тел Солнечной системы и обеспечивает требуемую точность на ограниченном интервале времени — порядка 100 дней. Данное ограничение связано с тем, что смещение исследуемого объекта на каждом шаге интегрирования проводилось по невозмущенным кепплеровским орбитам. Возмущение от больших планет учитывалось только в конце шага интегрирования. Однако при использовании данного метода для численного интегрирования уравнений движения на интервале времени порядка нескольких сотен лет следует учитывать изменение элементов орбит
Заусаев Анатолий Федорович — профессор кафедры прикладной математики и информатики- д.ф. -м.н., проф.
Заусаев Дмитрий Анатольевич — студент специальности «Прикладная математика и информатика».
на каждом шаге интегрирования. Для решения поставленной задачи необходимо из уравнений движения, приведенных в работе [1], получить дифференциальные уравнения движения для задачи двух тел и проинтегрировать их.
Рассмотрим сначала абсолютное движение двух тел относительно произвольной инерционной системы отсчета. Обозначим через р0 и р векторы тел 5 и М, определяющие их положения относительно начала координат О выбранной инерционной системы координат. Положение тела М относительно 5 будет определяться вектором г = р — ро. Тогда дифференциальные уравнения движения, приведенные в работе [1] в векторной форме, определяющие движение тел 5 и М, запишутся так:
(I2 ро 3асутт2т г
М2 г2 + г у/г3 — г3т + {/(г3 — г3т)2 г'
, (1)
й2 р ЗаазТ2
____ _ ______________^03 '- 05__________
(И2 г2 + Г у/г3 — г33 + у/(г3 — г33)2
где тот, тоз — эффективные радиусы тел М и 5- аот, аоз — соответствующие ускорения на расстоянии тот и то3, т — расстояние между центрами 5 и М.
Вычитая из второго уравнения (1) первое, получим уравнение, определяющее движение тела М относительно 5:
й2 г 3 аозТ^
(Ц2 ^ г2 _|_ г з/гз _ гз^ _|_ зД^з~^Г^Гу2
I____________Заотг2т_______________ г, .
г2 + г у/г3 — г3т + у/{г3 — г3т)2) г
Полагая, что тот и аот тела М пренебрежимо малы по сравнению с тоз и аоз тела 5, что справедливо, если в качестве тела 5 рассматривать Солнце или большую планету, а в качестве тела М — астероид или комету, получим уравнение движения тела М относительно 5 в следующей форме:
с12г 3 а03г23 г ^
& lt-2 Г2 + Г у/Г3 — Г33 + у/(Г3 — Г33)2 Г
«3
Ограничиваясь членами первого порядка относительно % в разложениях
1 — 3 и1 — 3 в бесконечный ряд, уравнения (3) можно представить
Гр з
— X Г
в виде
(12г аоГ203 Г (г3аз
1 + ^ • (4)
2 т3 3т3 /
3
Полагая, что а0г23 = ц ^ = е, запишем уравнение (4) следующим образом:
(12г _^г/ е
Д2 — 7? (.! + рО • (5& gt-
Как показано в работе [3], из уравнений (5) следует, что возмущения в большой полуоси, эксцентриситете и в движении линии апсид можно вычислить по следующим формулам:
1а 2ее 3
¦ 8Ш р (1 + е 008 р) ,
1р а2(1 — е2)5
de е. -з. ч
¦ 81Пр (1 + в СОвр), (6)
1р а3(1 — е2)3
1ш е, 3
7 Г =Гл-------23 С°8^(1 +еС08(^) ,
ар еа3(1 — е2)3
где, а — большая полуось, е — эксцентриситет, р — истинная аномалия возмущаемого тела, ш — аргумент перигелия.
Изменение элементов орбит за время одного полного оборота находится интегрированием по полярному углу р (истинная аномалия).
Легко убедиться в том, что большая полуось и эксцентриситет не имеют вековых изменений, в то время как линия апсид в течении каждого оборота поворачивается в прямом направлении на угол
е г2п
= еа3(1 _ е2)3 ] С08р (1+еС08р)3^, (7)
величина которого с учётом членов первого порядка относительно эксцентриситета такова:
. 3пе
Аш =
а3(1 — е2)3'-
При соответствующем выборе постоянной е можно получить удовлетворительное количественное объяснение невязок в движении линии апсид планетных орбит.
Таким образом, эффект смещения линии апсид в прямом направлении, как следует из (7), должен наблюдаться в невозмущенной задаче двух тел. При численном интегрировании уравнений движения малых тел Солнечной системы методом оскулирующих элементов начальными данными являются координаты х, у, г и скорости X, у, г исследуемого объекта, а также координаты и скорости больших планет хг, у г, гг] X г, у г, гг в экваториальной системе координат на момент времени Ьо.
Нахождение элементов орбит по координатам и скоростям. По известным координатам и скоростям находим элементы орбит по формулам [4]
1
а
2 V2
(8)
где, а — большая полуось орбиты- г = л/ х2 + у2 + г2, V2 = х2 + у2 + г2, ц = 2
= ао г^.
Эксцентриситет е и эксцентрическая аномалия Е на момент времени Ьо вычисляются из соотношений
гг г
евтЕ =--, есовЕ=1-, (9)
а а
где rr = xx + yy + zz. Затем с помощью уравнения Кеплера
M0 = E — e sin E (10)
вычисляется средняя аномалия в эпоху to.
Из соотношений
л/]мр sin i sin Q = yz — zy,
sin г sin cos Q = xz — zx, (11)
¦sjJlpcosi = xy — yx,
находятся наклонение i и долгота восходящего узла Q, отнесенные к экватору, а также параметр p.
Аргумент перигелия и находится по формулам
Jprf zcos г
tgv = -rf-----tg-u —
& amp-(р — г)' х cosQ + у sinQ'
ш = и — V, ш = ш + Аш. (12)
Вычисление координат и скоростей по элементам орбит. Экваториальные координаты х, у, г и скорости X, у, і на следующем шаге Л, = і - іо вычисляются с помощью следующих формул [4]:
М = п (і - і0) + М0, (13)
Е — е sin Е = Ж, (14)
v 1 + e E
Ч2 = Уі^-еЧ2'- (15)
г = (16)
1 + e cos v
и = v + и, (17)
x = r (cos u cos Q — sin u sin Q cos i),
y = r (cos u sin Q +sin u cosQcos i), (18)
z = r sin u sin i.
Здесь M — средняя аномалия- n = -среднесуточное движение, to-на-
чальный момент времени (эпоха), Mo — средняя аномалия в эпоху, E — эксцентрическая аномалия, v — истинная аномалия, u — аргумент широты.
Пусть V — скорость, Vr -радиальная скорость, Vk -трансверсальная скорость. Тогда
У2 = !л (19)
r а)
Vr =. -es'-mv, (20)
V p
Vn = , — (1 + ecosv). (21)
p
Для нахождения экваториальных координат скорости продифференцируем по времени формулы (18), получим
x
х = - Vr + (- sin и cos Q — cos и sin Q cos i) Vn, r
y
у = - Vr + (- sin и sin Q + cos и cos Q cos i) Vn, (22)
z
z = - Vr + cosusmiVn. r
Для учёта возмущающего действия от больших планет уравнения движения возмущаемого тела следует представить в гелиоцентрической системе координат. При переходе от барицентрической системы координат к гелиоцентрической уравнения движения возмущаемого тела, приведенные в работе [1], с учётом формулы (3) запишутся в виде
d2x x + x,
dt2 rp 2 _|_ rp f 3 — r|s + /(r3 r3) 2 1 os) r
d2y _ 3ci0o^ os У + Y
dt2 rp2 _|_ rp ^/^3 -rls+V (r* r3) 2 — os r
d2z За00т os z +
dt22 _|_ rp ^/^3 — r|s + /(r3 r3) 2 — os r
где
V _ (Х^-Х _ ________ 3ао^Го^ I х±________ 3ао^Го^
г Аг Д?+Дг ^/(А|-гЗ.)2 п г2+г. 3/г3_г3.+ 3/(г3_г3.)2
ЛЛ _ |. _________3ао1Го1____________I Уг __________3ао1Го1_________
г ^ Д- Д2+Д. 3/Д3_г3.+ 3/(Д3_Г3.)2 п г2+п з/гз_гз.+ 3/(г3_г3.)2
7 = '- .(У^У. __________ За°^., ^_________ 3ао^ы,
Д2+Д. 3/Д3_г3.+ 3/(Д3_Г3.)2 Г- г2+г. ЗДЗ_Г3.+ 3/(Г3_Г3.)2
А2 = (хг — х)2 + (уг — у)2 + (гг — г)2- тог — эффективный радиус г-того тела- а0г — соответствующее ускорение для г-того тела на расстоянии т0г от центра массы- X, У, 2 — возмущающие ускорения от больших планет.
Исправляя найденные координаты и скорости за счет возмущающего действия от больших планет, окончательно получим координаты и скорости на следующем шаге интегрирования:
Л2 Л2 Л2
хп+1 = х'-п+1 + -X, уп+1 = у'-п+1 + у У, гп+1 = г'-п+1 + у г, (24)
Хп+1 = х? П+1 + ЛХ, у/п+1 = УП+1 + ЛУ, гп+1 = гП+1 + (25)
где хП+1, уП+1, ^П+1, ХП+1, у/П+1, -гП+1 -координаты и компоненты вектора скорости возмущаемого тела на п + 1 шаге, вычисленные без учёта возмущающего действия от больших планет- хп+1, уп+1, 2П+1, Хп+1, уп+1, гп+1 — координаты и компоненты вектора скорости возмущаемого тела на п+1 шаге, полученные с учётом возмущающего действия от больших планет и Луны.
Таким образом, координаты радиус-вектора и скорости возмущаемого тела в точке Ьп+1, если они известны в точке Ьп, методом оскулирующих элементов вычисляются следующим образом:
а) по формулам (8)-(12) находятся элементы орбит исследуемого объекта на момент времени tn-
б) с использованием формул (13)-(22) вычисляются координаты и скорости данного объекта на момент времени tn+i, при этом исправляется величина аргументы перигелия по формулам (12) —
в) с помощью формулы (24) определяются окончательные координаты и компоненты скорости объекта на момент времени tn+i с учётом возмущающего действия от больших планет и Луны.
Эффективность предложенного метода можно существенно повысить, если для координат и скоростей больших планет использовать заранее полученный банк данных.
Работа выполнена при финансовой поддержке Федерального агентства по образованию (РНП. 2.1. 1/745).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Заусаев А. Ф., Заусаев Д. А. Эволюция методы интерполяции, используемые для получения координат и элементов орбит больших планет и малых тел Солненой системы // Вестн. Сам. гос. техн. ун-та. Сер. Физ. -мат. науки, 2008. — № 2(17). — C. 231−238.
2. Standish E. M. JPL. Planetary and Lunar Ephemerides, DE405/LE / In: Jet Lab Technical Report. IOM 321. F-048, 1998. — P. 1−7.
3. Богородский А. Ф. Всемирное тяготение. — Киев: Наукова думка, 1971. — 353 с.
4. Заусаев А. Ф., Заусаев А. А. Математическое моделирование орбитальной эволюции малых тел Солнечной системы. — М.: Машиностроение-1, 2008. — 250 с.
Поступила в редакцию 08/XII/2008- в окончательном варианте — 09/II/2009.
MSC: 85−08
INTEGRATION OF EQUATIONS OF SOLAR SYSTEM SMALL BODIES MOTION WITH OSCULATING ELEMENTS METHOD
A. F. Zausaev, D. A. Zausaev
Samara State Technical University,
244, Molodogvardeyskaya str., Samara, 443 100.
E-mails: Zausaev_AFamail. ru
Algorithm of osculating elements method for numerical integration of equations of Solar system small bodies ' motion is introduced.
Key words: interpolation, numerical integration, differential equations of motion, method of osculating elements.
Original article submitted 08/XII/2008- revision submitted 09/II/2009.
Zausaev Anatoliy Fedorovich, Ph. D. (Phys. & amp- Math.), Prof., Dept. of Applied Mathematics and, Computer Science.
Zausaev Dmitriy Anatol 'evich, Student, Dept. of Applied Mathematics and Computer Science.

ПоказатьСвернуть
Заполнить форму текущей работой