Development of a new more precise algorithm for computing tidal Love numbers

全文:

Земля находится под постоянным приливным воздействием со стороны Солнца и Луны. Приливные силы вызывают деформацию фигуры и возмущение гравитационного поля планеты, величины которых могут быть измерены с помощью спутников. В [1] предложены безразмерные геофизические параметры, называемые числами Лява, которые описывают эти изменения: число k пропорционально возмущению гравитационного поля, а числа h и l пропорциональны соответственно радиальному и тангенциальному смещениям поверхности планеты. Эти числа являются интегральными величинами, зависящими от распределения плотности, модуля сдвига и модуля сжатия в недрах планеты.

Приливные числа Лява особенно важны при исследовании тел, для которых еще не проведены сейсмические эксперименты. Вычисляя числа Лява для различных моделей некоторой планеты и сравнивая модельные значения с измеренными, можно получить ценную информацию об ее недрах. С их помощью изучают, например, внутреннее строение Венеры [2, 3], Меркурия [4] и Марса [5].

Для корректного применения чисел Лява в исследовании недр разных тел необходимо уметь их вычислить с высокой точностью для всех рассмотренных моделей. В литературе встречаются различные подходы вычислений, каждый со своими недостатками. В данной работе обсуждаются существующие алгоритмы, на их основе мы разработали и подробно описали новый, более точный и корректный алгоритм, позволяющий вычислить числа Лява для сферически симметричных моделей различных космических тел.

СИСТЕМЫ УРАВНЕНИЙ

В [6] выведена система уравнений для расчета приливной деформации по радиусу r, исходя из уравнений движения упругой среды, уравнения Пуассона для гравитационного потенциала и условия гидростатического равновесия. Рассмотрим сначала случай твердых слоев (например, кора, мантия и внутреннее ядро).

Вводятся 6 переменных где i = 1, 2,…,6, а n обозначает порядок (в [6] все рассмотренные величины разлагаются по сферическим функциям): радиальное и тангенциальное смещения y_1n, y_3n; нормальное и сдвиговое напряжения y_2n, y_4n; гравитационный потенциал y_5n и модифицированная функция гравитации y_6n. Чаще всего используются числа Лява второго порядка и, для простоты, индекс n далее опускается. Переменные y₂, y₄ и y₆ связаны с остальными следующим образом:

,

,

, (1)

,

где X – дилатация, λ и μ – параметры Ламе.

В коре, мантии и во внутреннем ядре Земли будем интегрировать приливную деформацию с помощью дифференциальных уравнений (2), выведенных в [6] без каких-либо изменений:

 

 

 (2)

 

 

 

где χ – частота прилива (равна обратному периоду, умноженному на 2π), а ρ₀ и g₀ – невозмущенные значения плотности и гравитационного ускорения на расстоянии r от центра планеты.

Система (2) изначально была создана для моделирования собственных колебаний Земли, но ее также можно применить к расчету чисел Лява. Периоды приливов как минимум на порядок больше, чем периоды собственных колебаний, и поэтому при расчете чисел Лява часто пренебрегают слагаемыми, пропорциональными χ² (например, в [7]). При частоте полусуточного лунного прилива M₂ (12 ч 25 мин) мы оценили, что эти слагаемые в сотни раз меньше остальных. Действительно, наши расчеты показывают, что без этих слагаемых числа Лява уменьшаются примерно на 0.3%. Мы считаем, что нецелесообразно ими пренебрегать в уравнениях (2) по двум причинам: это создает значительную дополнительную погрешность в числах Лява, но при этом никак не уменьшает время расчета.

Во внешнем ядре (или в любой среде с μ = 0) система (2) неопределенная. В [6] была предложена другая система для интегрирования возмущения Земли через внешнее ядро. В ней тангенциальное смещение пропорционально χ^–2 и при достаточно больших периодах y₃ может принимать нефизически большие значения. Эти уравнения подходят для расчета собственных колебаний, но не для вычисления чисел Лява (особенно для длинных приливов).

В [6–8] предложены разные подходы интегрирования приливного возмущения Земли через внешнее ядро. В этих работах пренебрегают слагаемыми с χ² в жидких слоях, и принимаемые условия в ядре различны. В [8] при выводе уравнений считают, что выполняется условие Адамса–Вильямсона: =, где точка обозначает производную по r. Однако в моделях Земли (например, PREM) это условие не удовлетворяется на протяжении всего внешнего ядра. На рис. 1 изображены значения величин и согласно PREM и видно, что они заметно отличаются у вершины внутреннего ядра и у подошвы мантии. Действительно, современные исследования указывают на присутствие химической стратификации [10] в жидком ядре, что противоречит условию Адамса–Вильямсона.

 

Рис. 1. Сравнение величин и во внешнем ядре. Условие Адамса–Вильямсона не выполняется у границы внутреннего ядра и у границы мантии.

 

Другой подход применяется в [7]: вместо условия Адамса–Вильямсона считается, что дилатация нулевая во внешнем ядре. Это дает дополнительное условие (), которое должно удовлетворяться на границе внутреннего ядра и на нижней границе мантии. Из-за этого задача становится переопределенной и в [7] вводится разрыв радиального смещения на этих границах для замыкания системы. Алгоритм решения задачи при таком подходе изложен в [11], и он действительно дает хорошие результаты (достаточно близкие к результатам, полученным при использовании разработанного нами алгоритма). Тем не менее физического обоснования для разрыва радиального смещения нет, и мы предпочитаем отказаться от такого подхода.

Для упрощения записи граничных условий переменные y_i в жидком ядре далее будут обозначаться как z_i. В работе [9] предложен подход, в котором не используется условие Адамса– Вильямсона и не пренебрегается дилатацией во внешнем ядре, а вводится новая переменная z₇:

. (3)

Обратим внимание на то, что в [9] используется другое определение “шестой переменной” по сравнению с [6]. Для простоты мы предпочитаем определить z₆ так же, как в [6] и в уравнениях (1) и (2) для твердой среды:

. (4)

Из уравнений (3) и (4) получим более удобное определение для переменной z₇:

. (5)

Следуя [9], будем использовать следующие дифференциальные уравнения для интегрирования приливной деформации через жидкое внешнее ядро:

,

. (6)

Будем также использовать следующие соотношения, которые верны во всем внешнем ядре (см. [9]):

, (7)

НАЧАЛЬНЫЕ УСЛОВИЯ В ЦЕНТРЕ ПЛАНЕТЫ

В центре планеты нет смещений, и изменение потенциала нулевое: , и . Существуют тогда три независимых начальных условия для интегрирования от r = 0 до границы внутреннего ядра с помощью системы (2):

 (8)

В системе (2), однако, имеется сингулярность в нуле. Эту проблему можно обойти, если начать интегрирование векторами (8) не от нуля, а от маленького , как предложено в [11]. Используя равным одной тысячной радиуса Земли (примерно 6 км), получаем достаточно точные результаты для чисел Лява.

В нашей программе используем функцию solve_ivp библиотеки scipy [12] и сравниваем два метода интегрирования: явный метод Рунге–Кутты 8-го порядка (DOP853) и неявный метод 5-го порядка (BDF). Описанный подход отхода от нуля имеет огромный недостаток – конечные результаты сильно зависят от используемого метода и от шага интегрирования. Это объясняется тем, что, подставляя (8) в (2), на самом первом шаге интегрирования невязка пропорциональна и может привести к большим погрешностям.

Для того чтобы увеличить точность при отходе от нуля, можно разложить переменные по r. Такой подход был предложен в [13], на основе этой статьи мы разработали свою собственную аппроксимацию – , и аппроксимируются около нуля следующим образом:

,

, (9)

,

где и A, B, C, b, c – неизвестные коэффициенты. Подставляя (9) в (1), можно выразить , и через эти же коэффициенты:

 

,

, (10)

,

здесь считается, что , и – постоянные около нуля.

В системе (2) три уравнения (1-е, 3-е и 5-е) следуют из определений , и и выполняются тождественно при (9) и (10). Подставляя (9) и (10) в оставшиеся три уравнения, получим систему алгебраических уравнений. Все слагаемые, пропорциональные и могут быть исключены, если:

, (11)

b = A∙ J, c = A ∙ nJ.

Таким образом, невязка в системе (2) на первом шаге пропорциональна и не зависит от коэффициента A. Решение теперь зависит от трех неизвестных коэффициентов – A, B, C.

Из (9)–(11) получим три независимых вектора, которые с высокой точностью аппроксимируют решение около нуля:

 

 (12)

 

а само решение при равно

 (13)

где Α, Β, C – неизвестные коэффициенты.

Расчеты показывают, что решения, в которых используются начальные условия (12) вместо (8), значительно более устойчивы к изменению метода и шага интегрирования. Поэтому мы считаем, что при вычислении чисел Лява любых планетных тел с внутренним ядром следует применять полученные нами начальные условия (12).

АЛГОРИТМ ИНТЕГРИРОВАНИЯ

Первый шаг – интегрирование векторов , и от до границы внутреннего ядра (ICB), используя систему (2). Получаем тогда три вектора , и , а полное решение под этой границей имеет вид

 (14)

Решение над этой границей (уже в жидком внешнем ядре) обозначим как . При переходе от твердой к жидкой среде требуем непрерывности всех компонент за исключение третьей [6]:

,

 

,

, (15)

 

 

где – плотность в жидком внешнем ядре непосредственно над границей ICB, а – ускорение силы тяжести на этой же границе. В (15) используются соотношения (7).

Получаем тогда 5 уравнений с 6 неизвестными: A, B, C, , и . Из первых четырех уравнений можно выразить B и C через A следующим образом:

 

, (16)

 

Решение над ICB равно:

 (17)

где

 

 (18)

и коэффициент A пока еще неизвестен.

Для интегрирования через жидкое внешнее ядро нам нужны только две компоненты – пятая и седьмая. С ними можно составить следующий вектор:

 (19)

и его следует проинтегрировать до границы внешнего ядра с мантией (CMB) с помощью уравнений (6). Результат обозначим как

 (20)

а полное решение под CMB есть

 (21)

Для того чтобы продолжить интегрирование из CMB до поверхности, необходимо извлечь из (20) значения второй и шестой компонент, используя соотношения (5) и (7):

,

 (22)

где – плотность в жидком внешнем ядре непосредственно под CMB.

Подставляя (21) в (22), можно выразить и через два неизвестных – Α и :

,

(23)

В CMB требуем непрерывности всех компонент, кроме третьей. Решение непосредственно над CMB (в подошве нижней мантии) имеет вид

 

 (24)

и зависит от трех неизвестных – A, и . Для простоты переобозначим и как D и E соответственно, а три вектора из (24) обозначим как , и . Тогда решение над CMB можно записать в виде:

. (25)

Следующий шаг – проинтегрировать векторы , и от CMB до поверхности планеты, используя уравнения (1). Получаем три вектора, которые обозначим как , и и полное решение на поверхности:

 (26)

Пренебрегая влиянием атмосферы в приливной деформации Земли, получаем следующие граничные условия на поверхности (r = a) [6, 11]:

 (27)

Это приводит к системе 5 уравнений с 5 неизвестными (A, D, E, и ):

,

,

, (28)

 

,

решение которой получается из

. (29)

Подставляя полученные значения коэффициентов A, D и E в (26), получим все компоненты . Приливные числа Лява , , определяются как

 (30)

Описанные уравнения необходимо сначала обезразмерить. Все интегрирование выполняется с помощью функции solve_ivp библиотеки scipy [12] на языке программирования Python. Методы RK45, DOP853 и BDF дают одинаковые результаты для выбранной точности вычисления. На первом шаге используется .

ОБСУЖДЕНИЕ И ВЫВОДЫ

В данной работе алгоритм применяется к вычислению приливных чисел Лява Земли, но его можно использовать без изменений для любого тела с похожей внутренней структурой – внутреннее ядро, внешнее жидкое ядро, мантия и твердая кора. Для работы с другими телами необходимо сделать некоторые изменения в алгоритме. Если, например, ядро находится полностью в жидком состоянии, следует применить начальные условия, описанные в [9], и пропустить первый этап интегрирования через внутреннее ядро (т.е. начать с (19)). Если же присутствуют другие жидкие слои, достаточно в них повторить описанный выше способ интегрирования приливной деформации через внешнее ядро.

Алгоритм разработан для расчета чисел Лява в сферически симметричных моделях. Ожидаемая погрешность будет порядка эллиптичности исследуемого тела, в случае Земли она составляет 0.3%, что заметно меньше точности измерения чисел Лява. По этой причине в литературе встречаются только сферически симметричные алгоритмы вычисления чисел Лява.

Применение данного алгоритма к упругой модели Земли PREM при частоте полусуточного лунного прилива дает следующие значения чисел Лява:

 

 

что хорошо согласуется, например, с [13]. Однако планеты не являются идеально упругими телами, и при расчете чисел Лява необходимо учитывать неупругость путем применения некоторой реологии, как, например, реология Андраде. Значение модуля сдвига в каждом слое меняется по заданному закону и преобразуется в комплексную величину, но алгоритм при этом никак не меняется. В связи с этим надо использовать методы интегрирования, которые работают в комплексной области как RK45, DOP853 и BDF из библиотеки solve_ivp.

Описанный выше алгоритм является не только более точным и устойчивым для вычисления приливных чисел Лява, но и более корректным с точки зрения физики по сравнению с теми, которые часто встречаются в литературе. По этой причине именно этот алгоритм следует применять для расчета приливных чисел Лява. В случае если внутренняя структура исследуемого тела сильно отличается от земной, метод расчета легко адаптируется.

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Исследование выполнено при поддержке Российского научного фонда (проект № 23-22-00074), https://rscf.ru/project/23-22-00074.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы данной работы заявляют, что у них нет конфликта интересов.

作者简介

D. Amorim

Moscow Institute of Physics and Technology

编辑信件的主要联系方式.
Email: amorim.dargilan@gmail.com
俄罗斯联邦, Dolgoprudny, Moscow Region

T. Gudkova

Sсhmidt Institute of Physics of the Earth of the Russian Academy of Sciences

Email: gudkova@ifz.ru
俄罗斯联邦, Moscow

参考

补充文件