CC BY
46
Небесная механика и астрономия
УДК 523.642
ИНТЕГРИРОВАНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ МАЛЫХ ТЕЛ СОЛНЕЧНОЙ СИСТЕМЫ МЕТОДОМ ОСКУЛИРУЮЩИХ ЭЛЕМЕНТОВ
А. Ф. Заусаев, Д. А. Заусаев
Самарский государственный технический университет,
443100, Самара, ул. Молодогвардейская, 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 г ^
<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>
Как показано в работе [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Г = ^Гл-------2\3 С°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 —
&(р — г)’ х 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 +
dt2 ^2 _|_ 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, 443100.
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. & Math.), Prof., Dept. of Applied Mathematics and, Computer Science.
Zausaev Dmitriy Anatol ’evich, Student, Dept. of Applied Mathematics and Computer Science.
