ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ, 2014, том 40, № 6, с. 426-432
УДК 52-17
ИСПОЛЬЗОВАНИЕ ПЕРЕХОДА К ПЕРЕМЕННЫМ ДЕЙСТВИЯ НЬЮТОНОВОЙ ЗАДАЧИ В ЧИСЛЕННОМ РЕШЕНИИ ЗАДАЧИ N ТЕЛ
© 2014 г. К. В. Лежнин*, С. А. Чернягин**
Московский физико-технический институт, Долгопрудный Поступила в редакцию 21.10.2013 г.
Предложен подход для преодоления проблемы близких прохождений в столкновительных системах — шаровых и рассеяных звездных скоплениях. Как известно, шаг численного интегрирования в таких системах, например, при образовании тесных двойных звезд, начинает дробиться, и скорость расчетов падает вплоть до полной остановки. Показано, что с помощью теории возмущений в предложенном подходе можно изолировать сингулярность и существенно увеличить шаг интегрирования, не теряя физических эффектов, существенно влияющих на эволюцию звездных скоплений.
Ключевые слова: численное моделирование, задача N тел, звездная динамика. DOI: 10.7868/80320010814060035
ВВЕДЕНИЕ
Решение задачи многих тел в создаваемом ими гравитационном поле имеет большое значение для астрофизических приложений. Именно задача многих тел является основополагающей при решении различных задач звездной динамики, динамики и эволюции звездных скоплений, галактик и галактических скоплений. Как известно, задача многих тел не имеет точного аналитического решения, если число тел превышает 2. Во всех этих случаях возможно лишь приближенное численное решение. Рост числа тел ведет к существенному увеличению требований к вычислительным ресурсам. Известно два основных подхода к решению задачи многих тел: прямое решение дифференциальных уравнений движения и статистический метод. В рамках первого метода наибольшие трудности возникают при моделировании близких прохождений тел и возможных образованиях тесных гравитационно-связанных пар (двойных систем). При расчете динамики плотных шаровых скоплений учет образования таких пар необходим ввиду непосредственного влияния на физику процесса: такие двойные звезды при рассеивании на них остальных звезд в большом числе случаев становятся плотнее, передавая часть энергии рассеивающейся звезде. Именно за счет тесных пар происходит так называемый гравитермальный коллапс скоплений
Электронный адрес: klezknin@yandex.ru
Электронный адрес: chernyagin71@mail.ru
(Линден-Бэлл и др., 1968). Метод преодоления вычислительных трудностей, связанных с моделированием близких прохождений и образования тесных пар, рассматривается в данной работе.
Задача двух тел в создаваемом ими гравитационном поле имеет точное решение, восходящее к Ньютону, и рассматривается в любом курсе механики, (см. Ландау, Лившиц, 2007). Если количество тел превышает 2, приходится пользоваться вычислительными методами, все многообразие которых представлено в обзорных статьях (см., например, Денен, Рид, 2011). В настоящей работе мы будем рассматривать метод прямого решения дифференциальных уравнений движения.
Полная функция Лагранжа Ь для системы N гравитационно-взаимодействующих тел имеет следующий вид:
Ь = ^Ьг. (1)
Здесь Ьг — функция Лагранжа г-й частицы:
2 п N Ьг = тщ+С £ (2)
2 2 . , . / . Гг3 3=М = 3
где уг — модуль скорости г-го тела, тг — масса г-го тела, С — гравитационная постоянная, ггз - модуль расстояния между г-м и ]-м телами. Такая функция Лагранжа приводит к следующим уравнениям дви-
жения:
d2
dt2
= - G
j4j
N
j=1,i = j
r 3 ij
(3)
dt = \ Tj~—777
öiäi + ä2
ä i a i + ä2'
где аг — модуль ускорения г-го тела, точка обозначает производную по времени, п — безразмерный параметр, обычно выбираемый ^0.02 (см. Денен, Рид, 2011). Особое внимание в моделировании столкновительной динамики уделяется регуляризации близких прохождений. Понятно, что в случае точного равенства координат гравитационная сила обращается в бесконечность. Обращение с "большими" числами, тем более с бесконечностями, вызывает проблемы у вычислительной техники, поэтому для обхода этой проблемы часто используют "обрезание потенциала" (Арсет, Фолл, 1980; Хокни, Иствуд, 1987) либо так называемую регуляризацию (Орлов, Рубинов, 2008), устраняющую особенность в уравнениях движения. В данной же работе предлагается иной метод для обработки тесных прохождений — переход к переменным действия кеплеровой задачи (Бэйтс, 1992). Его суть заключается в рассмотрении движения тесной пары тел как аналитически известного решения задачи двух тел, взамодействие же с остальными телами рассматривается в качестве малого возмущения.
Рассмотрим гамильтониан системы точечных частиц, взаимодействующих гравитационно по Ньютону:
H
где гг — радиус-вектор г-го тела, г у — вектор, проведенный от г-го тела к j-му телу. При фиксированных начальных координатах и скоростях всех тел задача сводится к задаче Коши и имеет единственное решение. Итак, решение задачи N тел сводится к численному решению системы дифференциальных уравнений с фиксированными начальными условиями.
МЕТОДИКА РЕШЕНИЯ
Прямое решение данной задачи классическими методами сталкивается с большими вычислительными трудностями. Первым упрощением для решения этой задачи является индивидуальный выбор шага интегрирования для каждого тела (Денен, Рид, 2011). Этот прием имеет смысл тогда, когда моделируемая звездная система обладает сильной неоднородностью по плотности. В данном решении шаг выбирается согласно эмпирической формуле С Арсета:
N N
i=1
Vi
2mi
G 2
N N
i=1,j=1,i =
— rjг
(5)
Далее рассмотрим конкретное р-е тело и допустим, что ближайшее к нему имеет номер к. Выделим часть гамильтониана, содержащую особенность и назовем ее Н0:
Hn
Р2р + Р1
2m,
p
Gmpmk
2 тк |rp-rfc|'
(6)
С точки зрения аналитической механики, наиболее просто решение кеплеровой задачи выглядит в переменных действия (Кордани, 2000). Введем следующее действие:
S = jr Ur + je Щ + ЗфПф,
где jr, je и jф имеют вид
= ¿г /PrdT = 1 / / GMm\ L2,
(7)
(8)
(4)
Jd = ifPdde = iHL2-s^(e)
i
(9)
U = f Рф(1ф = Lz-
(10)
Сделаем следующее линейное преобразование:
jE = jr + je,
jL = je, (11)
jLz = ]ф-
Тогда действие представится в виде
S = jE ue + je ue + ,]фиф, (12)
где ul = ur — ue, ue = ur, uLz = иф; можем записать: ua = oüat, (а = E, L, Lz), oje = ыl =
= ulz =0. Связь эксцентрической аномалии £ с фазой ue принимает следующий вид (для эллиптической орбиты):
иЕ = £ — e sin(£).
(13)
mimj
j
Если записать уравнение Гамильтона—Якоби для гиперболической траектории в этих же переменных, то действие будет также линейно расти со временем, при этом фаза пе связана с эксцентрической аномалией для гиперболической орбиты следующим образом:
пе = С - е ).
(14)
Теперь допустим, что пара, образованная к-м и р-м телами, достаточно тесная (приливные силы от взаимодействия с другими телами много меньше силы взаимодействия двух выбранных тел). Тогда можно рассмотреть точные траектории движения этих двух тел как невозмущенное движение, влияние же окружающих тел учтем как малую поправку. При этом шаг интегрирования будет определяться возмущением V = Н — Н0 и может быть выбран существенно большим, нежели шаг при обычной методике разложения в ряд и применении формулы (4) для вычисления шага. Если сравнивать данный метод с КЗ-регуляризацией, то главным преимуществом рассматриваемого подхода является тот факт, что шаг интегрирования двойной системы при V — 0 стремится к бесконечности (определяется медленной подсистемой), чего нет при использовании КЗ-регуляризации (Орлов, Рубинов, 2008). Теперь перейдем к подробному описанию алгоритма. Реализованный код позволяет решать столкновительную задачу N тел, ниже для простоты опишем метод, реализованный для решения задачи трех тел.
АЛГОРИТМ Обозначения
Введем обозначения, которые будут использоваться в математической формулировке обозначенного выше алгоритма. Пусть ТгтеБ£ер — массив шагов интегрирования.
ТгтеБ£ерЕе8гГие — массив времен, содержащих разность между концом времени интегрирования данного тела и текущим временем программы. Будем называть передвигаемое за данный шаг тело г-м, доставляющее ему наибольшее ускорение — ]-м, в случае интегрирования в переменных действия оба тела будут перемещаться на данном шаге. X, V, Ш, Ш1, Ш2, Ш3 — массивы координат тел и производных координат по времени (1—5 производные).
Под собственным временем программы £ подразумевается следующая величина: в начальный момент времени £ = 0, на каждом шаге интегрирования она увеличивается на величину тгп(ТгтеБ£ерКе8гГие).
АйуапсеСоогйгпа1е(г, д£) и AdvanceVelocгty(г, в£) — процедуры, производящие вычисление координат и скоростей г-го тела через время dt при помощи формулы Тейлора:
Хг:=Хг + УгсИ + ]¥г—+ (15)
«г-, г;3 г;4 + \¥1г— + \¥2г—, 6 24
г;2 г;3
Уг:=Уг + \¥гсИ + \¥1г— + \¥2г—. (16)
2 6
Для учета возмущения от внешнего тела при интегрировании в переменных действия вводится функция Рет£итЬа£гопСоотГгпа£е(г,Г£), которая изменяет координату следующим образом:
Г£2 Г£3 Г£4
Хг := Хг + + \VlA- + (17)
2 6 24
Подробно ее использование будет описано ниже.
Инициализация переменных
Задаются начальные данные и параметры задачи — массы, координаты и скорости конфигурации тел, притом координаты выражаются в астрономических единицах, скорости — в скоростях вращения Земли вокруг Солнца, массы — в солнечных массах. Для такого выбора системы единиц гравитационная постоянная равна единице. Далее производится переход в сопутствующую систему координат — систему, в которой ХР = 0. Затем выполняется прямое вычисление ускорения и его производных по времени для всех тел (массивов Ш, Ш1, Ш2, Ш3). Оба массива времени — ТгтеБ£ер и ТгтеБ£ерКе8гГие инициализуются значением шага, вычисленного по упомянутой выше формуле (4).
Шаг интегрирования
Ищем тело с наименьшим ТгтеБ£ерКе8гГие, определяем его номер в массиве, именно это тело будем называть г-м. Находим Шг, Ш 1г, Ш2г, Ш3г — ускорения и его производные по времени, доставляемые г-му телу всеми прочими; среди них находим тело, которое доставляет г-му наибольшее по модулю ускорение — это ]-е тело; ускорение и его производные, доставляемые ]-м телом г-му, обозначаются тз, w1з, т23- и тЗз. Вычисляем временной шаг для г-го тела д
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.