Разработка численного подхода в теории физической либрации в рамках "главной проблемы"

Загидуллин Артур Александрович; Петрова Наталья Константиновна; Усанин Владимир Сергеевич; Нефедьев Юрий Анатольевич; Глушков Максим Вадимович

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА.

_ СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

2017, Т. 159, кн. 4 С. 529-546

ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)

УДК 521.14

РАЗРАБОТКА ЧИСЛЕННОГО ПОДХОДА В ТЕОРИИ ФИЗИЧЕСКОЙ ЛИБРАЦИИ В РАМКАХ «ГЛАВНОЙ ПРОБЛЕМЫ»

А.А. Загидуллин, Н.К. Петрова, B.C. Усанин, Ю.А. Нефедьев, М.В. Глушков

Казанский (Приволжский) федеральный университет, г. Казань, 420008, Россия

Аннотация

В работе проведено построение численной теории вращения Луны. Математическая модель вращения Луны рассмотрена в рамках «главной проблемы». Уравнения вращения получены на основе гамильтонова подхода. Полученные дифференциальные уравнения решались с помощью метода Рунге-Кутта 10-го порядка точности. Результаты проанализированы на основе остаточных разностей (между численным и аналитическим решениями). Амплитуда остаточной разности по долготе не превосходит по модулю 1.8 угл. сек., а по широте - 0.9 угл. сек. Это достаточно большое расхождение обусловлено неточностью начальных условий, приводящих к появлению ложных гармоник с большими амплитудами.

Ключевые слова: теория физической либрации Луны, главная проблема, уравнение Гамильтона, собственное значение, произвольные периоды, метод Рунге-Кутты, остаточная разность, резонансная частота

Введение

Актуальность исследования физической либрации Луны (ФЛЛ) связана с планируемыми в ближайшие десятилетия международными программами освоения Луны. Они включают в себя решение ряда фундаментальных проблем в изучении Луны, связанных с выявлением особенностей внутреннего строения спутника Земли, нахождение аргументов, подтверждающих существование жидкого ядра Луны, изучение его химического состава и температуры, уточнение размеров и степени сплюснутости (flattening) слоев ядра, изучение процессов тепловой и динамической эволюции Луны, в том числе поиск ответа на фундаментальный вопрос, почему наблюдаемое вращательно-приливное сжатие Луны больше, чем в теоретических расчетах (на основе современных данных о параметрах внутреннего строения Луны) эволюционных переходов к кассиневу вращению и т. д. Решение многих из перечисленных проблем, а также обеспечение навигационных задач в окололунном пространстве могут быть получены путем построения как можно более точной теории вращения Луны.

Настоящая работа посвящена построению численного решения задачи ФЛЛ в рамках задачи трех тел: Луна, Земля и Солнце. При разработке теории вращения Луны важно учитывать резонансный характер спин-орбитального движения Луны. В первом приближении такой тип движения Луны описывается эмпирическими законами Кассини:

1) наклон экватора Луны к эклиптике постоянный;

2) полюс эклиптики лежит между полюсами экватора и орбиты Луны;

3) Луна вращается с постоянной угловой скоростью в прямом направлении вокруг оси, обратно прецессирующей относительно полюса эклиптики, при этом период обращения относительно восходящего узла равен драконическому периоду.

Тип движения по Кассини - это результат длительной приливной эволюции Луны, приведшей к появлению устойчивого вращения с минимальной потерей вращательной энергии. По сравнению с Землей Луна вращается медленно, совершая один оборот вокруг оси за 27.3 сут. Такой тип вращения называется либрацион-ным [1]. В отличие от Земли, для которой характерен ротационный тип вращения, кинетическая энергия Луны достаточно мала по сравнению с работой моментов внешних сил, вследствие чего Луна испытывает лишь небольшие периодические колебания (либрации) относительно устойчивого равномерного вращения. Либрации вызваны рядом причин: это и несферичность лунного тела, и неравномерность распределения масс в нем, и наличие в центре Луны небольшого ядра, отличающегося от мантийного вещества агрегатным состоянием, химическим составом и скоростью вращения. Зависимость либраций от характеристик внутреннего строения Луны обеспечивает возможность изучения этих характеристик путем астрономических наблюдений за особенностями лунного вращения без привлечения дорогостоящих космических экспериментов. Такие эксперименты, помимо всего прочего, предполагают прямое изучение строения Луны с помощью сейсмических, магнитометрических и прочих прямых измерений параметров внутренней структуры Луны. Относительная дешевизна и простота астрономических методов наблюдений объясняет постоянный интерес к теории ФЛЛ, чем и обусловлена актуальность настоящего исследования.

В настоящей работе рассматривается приближенная модель ФЛЛ, называемая в литературе «главной проблемой». Суть её в следующем:

1) тело Луны есть абсолютно твердое;

2) орбитальное и вращательное движение Луны не зависимы друг от друга, что позволяет использовать уже готовые эфемериды Луны, построенные ранее другими исследователями;

3) источниками возмущения равномерного вращения (либраций) являются только Земля и Солнце, причем размеры их не учитываются.

Такая модель является базовой, поскольку она обеспечивает описание вращения Луны с достаточно высокой точностью - примерно 10 мс дуги, что, с одной стороны, существенно превышает точность наблюдения ФЛЛ наземными методами, а с другой - отвечает точности первых лазерных измерений ФЛЛ, начавшихся в 70-е годы прошлого столетия. Современные лазерные данные обеспечивают точность в определении ФЛЛ порядка 0.1 мс дуги, что и предъявляет повышенные требования к теоретическому описанию ФЛЛ. Но в любом случае дальнейшее улучшение теории невозможно без решения «главной проблемы».

Перечисленные выше существенные положения модели «главной проблемы» позволяют упростить получение уравнений ФЛЛ и их решение. Важно и то, что в данной модели мы можем рассматривать систему Земля - Луна как консервативную, тем самым не учитывать небольшие потери механической энергии (диссипации) вследствие приливного трения из-за вязкоупругих свойств лунного тела и взаимодействия лунного ядра с мантией. Хотя, конечно, диссипация лунного вращения происходит (о чем неопровержимо свидетельствуют наблюдения), в частности, вследствие приливной эволюции Луна отдаляется от Земли примерно на 4 см/год. Все прочие, более тонкие эффекты, вызывающие либрации на уровне 1 мс, в рамках главной проблемы можно не рассматривать. Полученное решение -первый, необходимый этап в построении высокоточной теории ФЛЛ. Мы предлагаем собственный метод решения этой задачи для того, чтобы, во-первых, иметь

собственный инструмент исследования вращения и внутреннего строения Луны и, во-вторых, приблизиться к объективной оценке точности теоретического описания вращательной динамики Луны, что, как правило, достигается путем сравнения решений, полученных на основе разных подходов, между собой.

Таким образом, точная теория ФЛЛ, как и других тел Солнечной системы, является, с одной стороны, одним из эффективных инструментов изучения внутреннего строения удаленных небесных объектов. С другой стороны, такая теория может быть использована для решения навигационных задач в окололунном пространстве.

Применение модели твердого тела с однородной плотностью для Луны оправдано, так как, во-первых, значение среднего момента инерции Луны близко к 0.4 -моменту инерции твёрдого однородного шара; во-вторых, гипотеза магматического океана [2] позволяет считать, что плотность мантии постоянна и вся неоднородность по плотности находится в тонком слое коры и небольшом ядре. В-третьих, время релаксации Максвелла (т = п/л, где п - вязкость, а л - модуль сдвига, при этом средние значения п = Ю25 пуаз, ¡л = 250 ГПа) для мантии составляет несколько тысячелетий, а значит, при наблюдениях в течение нескольких веков тело Луны будет вести себя как твердое тело, и вязкостью можно смело пренебречь [3].

Следовательно, раз тело Луны в большей степени является твёрдым, то главные моменты инерции Луны жёстко связаны с ним, и это, в свою очередь, даёт нам основание использовать триэдр главных осей как систему координат, положение которой в пространстве и будет описывать вращение Луны. Систему осей инерции Луны с началом в центре масс называем динамической системой координат (ДСК). В качестве инерциальной системы координат (ИСК) выбираем эклиптическую систему координат эпохи Л2000 (рис. 1). Для того чтобы выразить положение тела Луны в пространстве, мы используем матрицы вращения (П (д) - поворот относительно оси г на угол д) на углы л, V, п при переходе от ИСК к ДСК:

Уравнение (1) определяет положение тела Луны относительно инерциальной системы координат через углы либрации л, V, п. Эти углы можно рассматривать как разновидность самолётных углов. Наша задача состоит в нахождении зависимости этих углов от времени.

Для математического описания физической либрации используем метод Гамильтона. Для консервативных систем гамильтониан не зависит явно от времени и является по существу обобщенной энергией системы, то есть суммой потенциальной и кинетической энергий системы, выраженных в виде функций координат и импульсов. Примем в качестве обобщенных канонических координат д = (д\, д2, дз) углы либрации: д\ = ¡л, д2 = V, дз = п. Обобщённый импульс, канонически сопряжённый координате д^, определяется через функцию Лагранжа Ь данной механической системы р^ = дЬ/дд^. Функция Лагранжа задается по формуле Ь = = Т — и, где Т и и - кинетическая и потенциальная энергии. Соответственно, функция Гамильтона Н связана с функцией Лагранжа через преобразование Ле-жандра следующим соотношением:

1. Математическое описание ФЛЛ

(1)

Н = р ■ д — Ь.

(2)

Z (Z)

х (у)

Рис. 1. Селеноцентрические системы координат: ХУZ - эклиптическая система координат, ось X направлена в точку весны (7), Z - к полюсу эклиптики; (X, У, Z) -равномерно вращающаяся эклиптическая система координат со скоростью п = йЬ/№; Ь - средняя долгота Луны. Триэдр (х, у, г) - динамическая система координат: ось х ориентирована в направлении наименьшего момента инерции А; г - в направлении наибольшего момента С; направление оси у выбрано так, что образует правостороннюю декартову систему координат; /, V, п - углы либрации

Поскольку потенциальная энергия не зависит от обобщенных скоростей, обобщенные импульсы можно вычислить по формуле pi = dT/dc¡i.

2. Выражение для кинетической энергии системы

Начнем построение с выражения для кинетической энергии. Абсолютная угловая скорость вращения Луны wabc (относительно ИСК) состоит из суммы двух типов скоростей: а) переносная скорость ^пер соответствует равномерному вращению, синхронизированному с орбитальным движением; б) относительная скорость Н соответствует либрационному вращению относительно равномерного вращения. Получим формулы, которые определяют кинематику вращательного движения Луны, то есть зависимость угловой скорости вращения wabc от углов либрации и их скоростей (аналог кинематических уравнений Эйлера).

Для вывода уравнений либрации необходимо найти проекции вектора угловой скорости вращения Н на оси ДСК (рис. 1). Анализ геометрического взаиморасположения вектора Н и осей ДСК позволил получить следующие соотношения:

НХ = —¡1 sin v — П,

Ну = —¡¡cos v sin n — v cos n, (3)

Hz = ¡¡cos v cos n + v sin n.

Тогда кинетическую энергию T можно выразить через компоненты углового вектора по формуле:

T = 2 (АНХ + ВНу + CQl). (4)

При переходе от абсолютной скорости к относительной и с учетом особенностей переносного вращения из выражения для T следует вычесть переносной импульс Cn, тогда оставшийся импульс pi будет соответствовать только либрационным скоростям по долготе. В дальнейшем принимаем значение главного момента инерции C равным единице. Это означает, что все обобщенные импульсы измеряются в долях главного момента инерции.

В результате алгебраических преобразований из уравнения (4) с учетом (3) получаем следующее выражение для кинетической энергии:

T = 2((1 + &2)pÍ + (1 + ki)(p2 cos 92 - sec 92 (n + Pi - Рз sin 92))2)+

+ 2(p2 sin 93 + cos 93((n + pi) sec 92 - рз tg 92))2. (5)

Здесь ki = C/B — 1, k2 = C/A — 1, где A, B, C - главные моменты инерции Луны.

Выражение (5) для кинетической энергии позволяет нам построить уравнения Гамильтона и получить собственные частоты вращения Луны.

3. Выражение для потенциальной энергии системы

Гамильтониан (2) после преобразований обобщенных координат и импульсов может быть записан в виде H = T(9, р) — U(9), где U - селенопотенциал Луны.

Применим ставший классическим метод разложения потенциала Луны в ряд по сферическим функциям [4]. В силу простоты сферических функций этот метод удобен для аналитических и численных вычислений значения потенциала притяжения небесного тела для любой точки, имеющей массу m. Положение точки определяется её расстоянием от центра масс тела и направляющими косинусами (ui, U2, U3) радиуса-вектора данной точки. Согласно этому разложению потенциал представляется в виде суммы

U = U2 + U2S + U3 + • • • . (6)

Поскольку начало ДСК выбрано в центре инерции Луны, мы не учитываем первую гармонику селенопотенциала, а так как оси ДСК совпадают с главными осями инерции Луны, то значения коэффициентов Стокса есть Сц = 0 и 5ц = 0. Каждое слагаемое в выражении (6) соответствует степени описания формы тела и распределения масс внутри него: гармоники U2 (гравитационное взаимодействие Луны с Землей) и U2S (гравитационное взаимодействие Луны с Солнцем) описывают эллипсоидальную форму лунной динамической фигуры через отношения главных моментов инерции Луны или соответствующие им коэффициенты Стокса 2-го порядка C2m, S2m; гармоника U3 соответствует более сложной, похожей на грушевидную, форме лунного тела, распределение масс в таком теле определяется коэффициентами Стокса 3-го порядка C3m, Ss3m .

В рамках точности, определяемой предположением 3 «главной проблемы», ограничиваемся рассмотрением второй и третьей гармоник селенопотенциала при взаимодействии с Землей и второй гармоникой при взаимодействии с Солнцем.

Выражения для указанных гармоник селенопотенциала можно записать через направляющие косинусы тела в ДСК следующим образом:

U2 = 2 xGME (a )3((C - A)u2 - (B - C )u2), (7)

U3 = -2 ^ЗЕ MmR R (a) VC3iui + S3iU2 + C30U3 + 1OS33u32 2 a3 a \r/ V

5

- 3 C30% - 10C33u - 5C3iuiu3 - 10C32uiu3 - 20S32uu2u3 - 5S33^u3+

+ 30C33uiu2 - 30S33u2u + 10C32u2u3 ), (8)

U2S =3 GM 3((C — A)M?S — (В — C)U2s). (9)

2 a s r § /

В уравнениях (7)—(9) использованы следующие обозначения: а, as - средние расстояния от Земли до Луны и от Земли до Солнца; Mm , Me , Ms - массы Луны, Земли и Солнца; R - средний радиус Луны; А - множитель, учитывающий поправку к третьему закону Кеплера; G - гравитационная постоянная; ui, и, из - направляющие косинусы радиуса-вектора Луны во вращающейся эклиптической системе координат; uis, U2S, U3S - направляющие косинусы радиуса-вектора Солнца во вращающейся эклиптической системе координат.

Третью гармонику для описания взаимодействия Луны с Солнцем не учитываем, так как согласно оценкам работы [5], её вклад в потенциал по сравнению с U2 в миллион раз меньше. Оба возмущающих тела, Земля и Солнце, вызывающие либрацию Луны, рассматриваются как точечные массы. В формулах (7)-(9) зависимость U от углов либрации и неявно от времени учитывается через направляющие косинусы и радиусы-векторы Земли и Солнца в ДСК: (ui, U2, u3)T = (x,y,z)T/r в уравнении (1), для (a/r), (as/rs) , L, b, Ls, bs - на основании геометрических построений для радиусов-векторов Земли и Солнца и разложений этих величин в ряды по аргументам Делоне, задаваемых орбитальной теорией [16]

Построение уравнений Гамильтона

В случае консервативной системы гамильтониан H = T — U не зависит от времени явно, это позволяет записать для выбранных нами обобщенных координат и импульсов канонические уравнения в виде:

dqi dH dpi dH . 123 (10)

dt dpi dt dqi ' ' '

В правых частях уравнений (10) стоят функции только обобщенных координат и импульсов. Поэтому система уравнений (10) замкнута относительно этих переменных и представляет собой систему шести дифференциальных уравнений первого порядка, которые полностью определяют изменение во времени обобщенных координат q = (qi,q2,q3) и соответствующих им импульсов p = (р1,Р2,Рз), если заданы начальные условия, то есть значения координат и импульсов в момент времени t = 0.

Используя полученные выше выражения (5) для кинетической и (6)-(9) потенциальной энергии, запишем систему Гамильтона в явном виде:

dqi 2 pi cos2 q3 p2 cos q3 sin q3 sin2 q3 sin2 q3

— = n tg q2 +--2--ki--+ kin-2--+ pi-2--+

dt cos2 q2 cos q2 cos2 q2 cos2 q2

2 2 2 2 sin q3 cos2 q3 tg q2 sin q3 tg q^ sin q3 tg q2

+ kipi-2--p3--p3--kip3-,

cos2 q3 cos q2 cos q2 cos q2

dq2 2 2 — = (1 + ki)p2 cos q3 — ki cos q3 sec q2(n + pi — p3 sin q2)sin q3 + p2 sin q3,

= 1 (2(1 + k2)p3 + 2(1 + ki) sinq3(p2 cosq3 — secq2(n + pi — — p3 sin q2) sin q3) tg q2) — cos q3 tg q2 (p2 sin q3 + cos q3 ((n + pi) sec q2 — p3 tg q2)),

dpi = dU (11)

dt dqi

dp2 dU 2

—f = d92 - (1 + ki) sec 92(P3 - (n + pi) sin 92) sin 9з(р2 cos 93- sec 92 (n + pi - p3 sin 92) sin 93) - cos 93 sec 92 (-p3 sec 92 + (n + pi) tg 92 )p sin 93+

+ cos 93((n + pi) sec 92 - p3 tg 92)),

dp3 dU

— = ---+ ki(p2 cos 93 - sec 92 (n + pi - p3 sin 92) sin 93) x

dt d93

X (p2 sin 93 + cos 93((n + pi) sec 92 - p3 tg 92)).

Уравнения системы (11) являются нелинейными относительно неизвестных переменных 9 = (9i, 92, 93) и p = (pi,p2,p3). Получить их точное решение в аналитическом виде не представляется возможным. Поэтому для решения подобных систем применяются либо приближенные аналитические методы [5-7], когда в разложениях в ряды учитывается ограниченное число слагаемых, либо, как используем мы, численные.

4. Линеаризованная система Гамильтона. Собственные значения

Для понимания корректной оценки и анализа результатов нам необходимо знать собственные частоты описываемого процесса либрации. Для их приближенного вычисления линеаризуем систему (11). Для этого, учитывая, что значения обобщённых координат малы (91 ~ 10-4, 92 ~ Яз ~ 10-2), разложим тригонометрические функции в правых частях этих уравнений по степеням обобщенных координат и выделим члены, линейные относительно 9 и р .

Такое выделение возможно, так как производные ¿и/¿дп, п = 1, 2, 3, после разложения тригонометрических функций в ряд можно записать в общем виде следующем образом:

Фп(Ч) = ^ = Е К + + кёП) V1 я2-€ЯГ"1 - ¿пЯп. (12)

Щп

Здесь функции Qi2k получены на основе формул (6)—(9) для и с использованием направляющих косинусов щ, щз. Эти функции являются тригонометрическими рядами, реализующими зависимость от времени через аргументы Делоне, и могут быть записаны следующим образом:

Q = ¿ Ar в, Y, C3n, S3n, a, COS (kríl + kr2Í' + kr3F + kr^Dj, (13)

где Pr - коэффициенты, рассчитываемые на основе принятой модели гравитационного поля Луны. В нашем случае это модель LURE2 [8]. Числовые значения коэффициентов Ar и индексов kri, i = 1, 2, 3, получены из аналитического описания орбитального движения Луны [9]. Выполнение этих операций реализовано с помощью специализированного программного комплекса, позволяющего проводить аналитические операции с рядами.

После проведения несложных, но громоздких алгебраических преобразований система (11) может быть переписана в общем виде следующим образом:

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

qi = pi + ii (q,p),

q2 = (1 + ki)p2 - kinq3 + F2(q,p), q3 = (1 + k2)p3 - nq2 + F3(q,p),

р1 = + Q1(q),

Р2 = -п2Ч2 + прз + F5(q,p) + Q2 (д),

рз = -^п29з + к1 ПР2 + F6(q,p) + Qз(q).

(14)

Здесь Fi(q,p), г = 1, . . . , 6, - нелинейные функции, которые содержат члены, полученные путем дифференцирования выражения кинетической энергии (см. (11)), за вычетом линейных членов. Функции Фп(д), п = 1, 2, 3, также нелинейны и содержат члены, найденные путем дифференцирования выражения потенциальной энергии (12), за вычетом выделенных линейных членов с числовыми множителями ¿1, ¿2, ¿з. Эти коэффициенты образованы из постоянных членов в рядах (13). Их значения: ¿1 = -3.53 • 10~5, ¿2 = -9.85 • 10~5, ¿3 = -3.24 • 10~7 - рассчитаны на основе коэффициентов рядов [9], описывающих движение Луны и констант динамической модели ЬиИ,Е2, значения которых приведены в табл. 1.

Следуя методике решения подобного типа уравнений, выделим линейные относительно ц и р члены в уравнениях (14) и построим для них характеристические уравнения. В силу того, что кинетическая энергия Т явно не зависит от д1 (либрации по долготе), линейная часть системы распадается на две независимых подсистемы:

-X 1 -¿1 —X

-X -к1п 1 + к1 0

-п -X 0 1 + к2

п - ¿2 0 -X п

0 -к1п2 + ¿з к1 п -X

(15)

Раскрывая определитель и решая квадратное и биквадратное уравнения, получим корни, определяющие следующие собственные значения:

Х^2 = ± \id1i,

Хз_4 = ^2 = V7п2 - ¿2 - ¿3,

Х5,6 = ±изг,из = л/ккп2-^2^3.

(16)

(17)

(18)

Полученные собственные значения есть собственные частоты либрации Луны, так как выделенная нами линейная часть уравнений ФЛЛ описывает невозмущенное вращение Луны. Как видим, частота либраций по долготе и определяется только коэффициентом ¿1, полученным из возмущающего потенциала и. Луна «затянута» в резонанс с орбитальным движением, и если бы не было внешних возмущений, то либрации по долготе не было бы: «носик» Луны не отклонялся бы от среднего направления на Землю, задаваемого осью X .

Частоты и и из определяют движение полюса Луны, задаваемого главным моментом инерции С, относительно полюса эклиптики. Значения этих частот есть и2 « 0.230185, и2 « 2.31188 • 10-4.

Каждой частоте соответствует свой период: Т = 1057.48 сут (около 2.9 года) для либрации по долготе; Т2 = 27.2963 сут, что примерно соответствует сидерическому месяцу; Тз = 74.4 года. Все значения получены с учетом влияния постоянных членов от потенциала - коэффициентов ¿1, ¿2, ¿з в (16)-(17). Если бы Луна вращалась без учета притяжения Земли, значения частот и соответствующих им периодов были бы другими, особенно существенно изменился бы период Тз : его значение, соответствующее свободному вращению, было бы в два раза больше. Поэтому корректно называть значения и, и, из частотами не свободной, а произвольной либрации [10, 11]. Колебания на этих частотах могут возникать при одномоментном возмущении Луны, например при падении крупного метеорита.

0

Табл. 1

Значения констант, описывающих динамическую фигуру и орбиту Луны

7 0.00022737 /с 1.627905081537 С30 -1.04Е-05

в 0.00063126 /с 0.2308957197 С31 2.86Е-05

1о 2.355555743493 1о 6.240060126913 С32 4.8Е-06

1 0.2280271436 1 0.0172019696 С33 2.7Е-06

к1 0.000403798188 Б 5.1984665886 $31 8.8Е-06

к2 0.000631403620 В 0.21276871 $ 32 1.7Е-06

п 0.2299708345 Сгс -0.0002027 $33 -1.1Е-06

Па 0.017202423838 С22 2.23Е-05

Табл. 2

Начальные значения на момент Л2000 для интегрирования

910 0.001057405207433 920 -0.0269194915987523 930 0.00182899555868

Р1С -6.0963407117Е-05 Р20 3.63692296085Е-04 Р30 -9.7541461736Е-07

Среди гармоник вынужденной либрации имеются частоты, близкие к ш\, ш2, поэтому трудно выделить малые амплитуды произвольных либраций на этих частотах. Тем не менее долговременные лазерные наблюдения позволили определить амплитуды этих либраций с точностью в несколько миллисекунд [12]. Частота шз, с которой полюс С осуществляет ретроградное движение по эллипсу с осями 3 х 8 угл. сек. дуги, соответствует так называемым чандлеро-подобным качаниям, которые вызваны несовпадением главной оси инерции Луны с её осью вращения [13]. Обнаружение либраций с частотами, равными собственным, несет важную информацию о причинах поддержания и возбуждения произвольных колебаний.

5. Анализ результатов

Решение системы Гамильтона (11) проводилось с помощью метода Рунге-Кутты 10-го порядка точности. Интервал интегрирования составил 150 лет, а шаг интегрирования полагался равным 0.1 сут. Значения констант (в радианах), которые использовались в работе, приведены в табл. 1. В табл. 2 приведены начальные значения на момент Л2000 для численного интегрирования на основе вычислений по аналитическим рядам Н.К. Петровой [5].

Чтобы выявить точность построенного решения, мы вычислили остаточные разности между полученным нами численным решением и аналитическим решением Н.К. Петровой, перерасчитанным для модели гравитационного поля ЬиИ,Е2. Точность характеристик, полученных согласно теории Н.К. Петровой, составляет 0.01 угл. сек. на каждую гармонику.

Решение системы (11) имеет периодический характер: на рис. 2 и 3 представлены результаты на интервале 800 сут.

На рис. 2 приведено решение для либрации по долготе, отличительной чертой является наличие постоянного смещения, равного ~ 250 угл. сек., которое обусловлено грушевидностью тела Луны, описываемой селенопотенциалом 3-го порядка.

На рис. 3 показано решение ФЛЛ по широте. Графики сдвинуты по фазе относительно друг друга на 90° и имеют основную частоту, обусловленную обратным движением лунного узла с периодом, равным драконическому месяцу. Амплитуда остаточной разности соответствует изменению угла наклона между экватором Луны и эклиптикой « ±1.5° = ±5400".

50-

—I— 100

—I— 200

—I— 700

—I 800

300

400 Сутки

500

600

Рис. 2. ФЛЛ по долготе ^

Рис. 3. ФЛЛ по широте: слева для V, справа для п

На рис. 4 приведены графики остаточных разностей решений (между численным и аналитическим) на интервале примерно 150 лет.

На всем интервале интегрирования значения остаточных разностей не превышают по модулю 1.8 угл. сек. в либрации по долготе и 0.9 угл. сек. по широте. Аналогичная работа была проделана Г.И. Ерошкиным [14], который сравнивал на интервале 80 лет собственное численное решение с аналитическими решениями Д.Г. Экхарда и М. Мунс, на рис. 5 и 6 приведены остаточные разности только для долготы из работ [6, 7, 14].

Величину амплитуды остаточных разностей и наличие ряда характерных частот (движение лунного узла и собственные частоты системы) Г.И. Ерошкин объясняет [14] следующим образом: аналитические и полуаналитические решения имеют выражения для амплитуд малый знаменатель, который при близости вынужденных частот к собственным стремится к нулю, поэтому величина амплитуды остаточных разностей зависит от того, каким образом авторы аналитических решений находили эти решения в случае частот, близких к резонансным. Вероятно,

1.5 1.0 2 0.5

8 00 ш

ш -0,5

О

с;

* -1.0 -1,5 -2,0

о.е-

0,43

0,2'

8 о.о-

<0

3

О -0,2>: ■ -0,4-0,6-0,8-

10000 20000 30000 40000 50000 Сутки

20000 30000 Сутки

1,0-1 0,80,60,40,20,0-0,2-0,4-0,6-0,8-

Сутки

Рис. 4. Остаточные разности в углах ФЛЛ дг = д2 = V, дз = п

Рис. 5. Остаточные разности между решением Г.И. Ерошкина и Д.Г. Экхарда по долготе

аналогичная ситуация возникает и при интерпретации остаточных разностей при сравнении с результатами Н.К. Петровой. Чтобы проверить это, мы, используя частотный анализ соответствующих остаточных разностей, выделили частоты с наибольшими амплитудами. Построенные периодограммы Шустера [15] представляют собой квадрат спектра мощности (рис. 7).

На диаграммах отчетливо видны все свободные частоты, как и в работе Г.И. Ерошкина: первый свободный период составил 1059 сут, второй период -68 лет, третий период - 27.3 сут. Полученные резонансные периоды зависят от используемой модели гравитационного поля, то есть очень чувствительны к выбору параметров системы и, соответственно, будут различаться в зависимости от модели. Третий свободный период в 68 лет отличается от рассчитанного тео-

1."00

00 Сутки 20 000 30000

Рис. 6. Остаточные разности между решением Г.И. Ерошкина и М. Мунс по долготе

Т, -1059.57 дней

6 8 10 12 14 16 18 20 Частота [1/супси х Е-004]

2,01,51.00,50,0-

-1

Т2=27.3 дня

1 2 3

Частота [1/сутки х Е-002]

1,41.21,00,80,60,40,20,0-0,2 — -1

1 2 3

Частота [1/сутки х Е-002]

Рис. 7. Частотная зависимость для остаточных разностей

ретически и нуждается в дальнейшем улучшении. Возможно, что этот период в 68 лет является фиктивной ложной гармоникой, обусловленной неточностью выбора начальных условий для интегрирования дифференциальных уравнений. В целом анализ нашего решения совпадает с анализом Г.И. Ерошкина, что дополнительно свидетельствует о корректности полученного нами численного решения. Однако определить более строго точность полученного решения пока не удается, так как при вычислении остаточных разностей возникает сильный шум на резонансных частотах. Точность метода Рунге-Кутты 10-го порядка можно оценить путём сравнения решений для трех различных шагов 0.1, 0.05, 0.01 сут (рис. 8). Слева на рисунке показана разность решений для шагов интегрирования в 0.05 и 0.01 сут, а справа - разность решений для шагов в 0.05 и 0.1 сут. Анализируя полу-

Рис. 8. Точность метода Рунге-Кутта 10-го порядка для различных шагов. Слева разница решений между Ь = 0.05 и Ь = 0.01 сут, справа между Ь = 0.05 и Ь = 0.1 сут

ченные графики, можно сказать, что шаг, равный 0.1 сут, является недостаточным для современных данных, но достаточен для сравнения с теорией Н.К. Петровой. Для более точного решения мы использовали шаг, равный 0.05 сут, который обеспечивает внутреннюю точность метода в 10_ 10 угловой секунды на интервале 27 лет.

Нам представляется интересным анализ остаточных разностей при сравнении нашего решения с результатами аналитического представления ФЛЛ в работе Н. Рамбо и Дж.Г. Уильямса [12], поскольку данная теория является на сегодняшний день наиболее точным аналитическим описанием либрации.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Поскольку теория Н. Рамбо и Дж.Г. Уильямса строилась с использованием углов Эйлера, а не самолётных углов, мы имеем возможность сравнить остаточные разности только для направляющих косинусов полюса эклиптики (рис. 9).

Анализируя связь между углами Эйлера, в которых строилась теория Н. Рамбо и Дж.Г. Уильямса, и самолетными углами, используемыми в нашей работе, мы получили нелинейную связь между этими углами. Чтобы уменьшить ожидаемое расхождение между теориями, мы сравнивали направляющие косинусы полюса эклиптики. На рис. 9 показана остаточная разность на интервале примерно 150 лет.

Из рис. 9 видно, что амплитуда полученных нами остаточных разностей значительно больше, чем при сравнении с теорией Н.К. Петровой. Более того,

100

90

80

8

0) 40

-30-

1-'-1-'-1-'-1-'-1-1-1-

0 10000 20000 30000 40000 50000

Сутки

Рис. 9. Остаточные разности в направляющих косинусах полюса эклиптики

очевидно, что имеется систематический сдвиг между решениями. Эти расхождения естественны и обусловлены тем, что в теории Н. Рамбо и Дж.Г. Уильямса [12] учтено существенно больше эффектов в динамике и внутреннем строении Луны, чем в модели главной проблемы, а именно:

1) улучшена точность всех используемых констант, таких как средний момент инерции и др.;

2) учтено большее количество гармоник в разложении селенопотенциала (до шестой включительно);

3) учтены многие особенности внутренней структуры Луны: вязко-упругость мантии и жидкое ядро;

4) учтены прямые и косвенные возмущения от всех планет солнечной системы и ряд других факторов.

Все перечисленные причины расхождения нашей теории с результатами Н. Рамбо и Дж.Г. Уильямса определяют план дальнейшего развития представленного нами подхода в построении численной теории ФЛЛ.

В работе построены нелинейные дифференциальные уравнения, описывающие вращение Луны, с применением гамильтонова подхода. Модель спин-орбитального движения Луны построена в рамках «главной проблемы»: вращательное движение изучено отдельно от орбитального, тело Луны рассмотрено как абсолютно твёрдое и в выражении селенопотенциала учтены лишь вторая и третья гармоники при взаимодействии с Землей и вторая гармоника при взаимодействии Солнцем.

Для решения полученной системы уравнений построен интегратор для решения дифференциальных уравнений с заданной точностью, основанный на методе Рунге-Кутты 10-го порядка. Полученное решение сравнивалось с аналогичным аналитическим решением Н.К. Петровой [5]. Анализ остаточных разностей показал, что разности амплитуд не превышают 2 угл. сек., хотя точность теории Н.К. Петровой составляет 0.01 угловых секунды на одну гармонику. Основная причина обнаруженных расхождений вызвана, по-видимому, недостаточно точным расчетом резонансных членов в теории Н.К. Петровой. Поведение решения, полученного в нашей работе, близко к результатам поведения решения, полученного на основе численной теории Г.И. Ерошкина.

Заключение

Сравнение с наиболее точным на сегодняшний день аналитическим решением

ФЛЛ, полученным Н. Рамбо и Дж.Г. Уильямса, позволило определить пути дальнейшего развития нашего подхода в построении численной теории ФЛЛ.

Литература

1. Белецкий В.В. Движение искусственного спутника относительно центра масс. - М.: Наука, 1965. - 415 с.

2. Жарков В.Н. Внутреннее строение Земли и планет. Элементарное введение в планетную и спутниковую геофизику. - М.: ООО «Наука и образование», 2013. - 414 с.:

3. Жарков В.Н., Паньков В.Л., Калачников А.А., Оснач А.И. Введение в физику Луны. - М.: Наука, 1969. - 310 с.

4. Дубошин Г.Н. Небесная механика. Основные задачи и методы. - М.: Наука, 1976. -799 с.

5. Petrova N. Analytical extension of lunar libration tables // Earth Moon Planet. - 1996. -V. 73, No 1. - P. 71-99. - doi: 10.1007/BF00058046.

6. Eckhardt D.H. Theory of the libration of the Moon // Earth Moon Planet. - 1981. -V. 25, No 1. - P. 3-49. - doi: 10.1007/BF00911807.

7. Moons M. Analytical theory of the libration of the Moon // Moon Planets. - 1982. -V. 27, No 3. - P. 257-284. - doi: 10.1007/BF00929297.

8. King R.W., Counselman C.C.,III, Shapiro I.I. Lunar dynamics and selenodesy: Results from analysis of VLBI and laser data // J. Geophys. Res. - 1976. - V. 81, No 35. -P. 6251-6256.

9. Gutzwiller M.C., Schmidt D.S. The Motion of the Moon as Computed by the Method of Hill, Brown, and Eckert // Astronomical papers prepared for the use of the American ephemeris and nautical almanac. V. 23, pt. 1. - Washington: G.P.O., 1986. - 272 p.

10. Хабибуллин Ш.Т., Чиканов Ю.А. О произвольной либрации Луны в эйлеровском движении ее полюсов // Труды Казан. гор. астрон. обсерватории. - 1969. - № 36. -С. 49-60.

11. Rambaux N., Williams J.G. The Moon's physical librations and determination of their free modes // Celest. Mech. Dyn. Astron. - 2011. - V. 109, No 1. - P. 85-100. - doi: 10.1007/s10569-010-9314-2.

12. Petrova N., Gusev A., Kawano N., Hanada H. Free librations of the two-layer Moon and the possibilities of their detection // Adv. Space Res. - 2008. - V. 42, No 8. -P. 1398-1404. - doi: 10.1016/j.asr.2008.02.017.

13. Гусев А.В., Петрова Н.К., Ханада X. Вращение, физическая либрация и внутреннее строение многослойной Луны. - Казань: Изд-во Казан. ун-та, 2015. - 323 с.

14. Eroshkin G.I. Comparison of a Numerical Model of the Physical Libration of the Moon with two Semi-Analytical Ones // Figure and Dynamics of the Earth, Moon and Planets: Proc. Int. Symposium / Ed. by P. Holota. - Prague: Czech. Acad. Sci., 1987. - P. 685-696.

15. Теребиж В.Ю. Анализ временных рядов в астрофизике. - М.: Наука, 1992. - 392 c.

Поступила в редакцию 27.03.17

Загидуллин Артур Александрович, аспирант Института физики Казанский (Приволжский) федеральный университет