О движении внутреннего ядра Земли Текст научной статьи по специальности «Математика»

УДК 51-72

о движении внутреннего ядра земли

С.И. Исаева

Сибирский федеральный университет, Красноярск, Россия, sv.is@mail.ru

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

Ключевые слова: внутреннее ядро Земли, математическая модель, метод Мерсона

A mathematical model of the motion of the solid inner core of the Earth is discussed, based on the gravitational interaction of the Sun, the Moon, the Earth, and the core as a separate celestial body, on which additional force of the solid core' interaction with the other layers of the Earth. The choice of the numerical methods for solving the described system of ordinary differential equations.

Key words: Earth's inner core, mathematical model, the Merson's method

ВВЕДЕНИЕ

В течение последнего десятилетия проявился всплеск интереса к исследованиям внутреннего ядра Земли. Это связано,

• во-первых, с участием внутреннего ядра в процессе поддержания конвекции в жидком ядре, что влияет на магнитное поле Земли;

• во-вторых, достигнут значительный прогресс в сейсмологии, в частности, создание цифровой наблюдательной сети, обеспечившей проведение прецизионных измерений с применением средств вычислительной техники;

• в-третьих, большие успехи достигнуты в лабораторных испытаниях материалов при температуре и давлении, достигаемых в условиях внутреннего ядра Земли;

• в-четвертых, развиты вычислительные методы, позволяющие рассчитывать свойства материалов, слагающих недра Земли, при сверхвысоких давлениях.

В 1996 году была опубликована работа авторов Song X., Richard P.G, в которой приводилось сейсмологическое доказательство дифференциального вращения твердого внутреннего ядра Земли, незначительно отличающегося от вращения Земли в целом. На сегодняшний день изучение движения внутреннего ядра Земли осуществляется тремя способами. Первый основан на обработке сейсмических данных, второй - на анализе результатов лабораторных экспериментов, а третий базируется на методах математического моделирования.

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

Данная работа посвящена построению математической модели движения внутреннего ядра Земли как отдельного «небесного тела», выбору численного метода для решения полученной дифференциальной задачи, а также программной реализации выбранного метода.

математическая модель

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

Поскольку расчеты ведутся в гелиоцентрической системе координат, то Солнце считается неподвижной точкой с массой т8 , оказывающей гравитационное воздействие на три остальных объекта. Луна представляет собой точку с массой тм, двигающуюся в поле тяготения Солнца и Земли в целом. Земля представляет собой точку с массой тЕ, двигающуюся в поле тяготения Солнца и Луны.

Твердое ядро моделируется шаром фиксированного радиуса К1 со сферически симметричным распределением плотности р1 и массой т1. Его вращение вокруг собственной оси во внимание не принимается. Допускается смещение ядра от центра Земли, точнее, от центра сферы радиуса ^, ограничивающей поверхность Земли. Учитывается также его вращение вокруг оси суточного вращения Земли в целом с угловой скоростью (0Е. Поскольку движение твердого ядра осуществляется в

жидкой среде с плотностью ро, то принимается во внимание сопротивление окружающей жидкости

с коэффициентом вязкости ^о. Обозначим ГЕ, гм,

Г1 - радиус-векторы между Солнцем и центрами масс Земли, Луны и твердого ядра Земли, соответственно.

В рамках предложенной модели, на твердое внутреннее ядро Земли, как на объект Солнечной системы, действуют гравитационные силы Солнца и Луны:

F = -Gm m ——

1 s KJirigirij

- r — r F =—Gm m —_—

1 M ^'"m"^ I |3

V, — r..\ I M

I

Сумма гравитационных сил и сил гидростатического давления остальной части Земли на твердое ядро равна:

Г \

К = - 3 пОт1 Ро

1 -Р

(—1 - го )

V р1 У

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

К = -У0ЖК21 (¿1 - ¿Е - [¥Е , —1 - —Е ]) .

Сумма центростремительной и центробежной сил, обусловленных вращением внутреннего ядра вокруг оси Земли:

Ка = ¥ ■ т1

О/¥ в |2}

первого порядка. Для задания движения Луны и Земли в целом необходимы 12 начальных условий: 3 координаты и 3 компоненты скорости для каждого объекта, которые взяты нами из астрономических данных.

Для определения траектории движения ядра необходимы еще 6 начальных условий: 3 координаты и 3 компоненты скорости. В координатах задается начальное отклонение ядра от центра Земли - Аг :

г, = г„ + Аг

1 В

(3)

Компоненты скорости должны учитывать как его орбитальное движение вокруг Солнца, так и составляющую от вращения твердого ядра вокруг оси суточного вращения Земли ввиду его отклонения от условного центра Земли, т.е. начальные скорости движения ядра и геометрического центра Земли связаны соотношением

Таким образом, итоговая система дифференциальных уравнений, описывающая в рамках предложенной модели движение Земли, Луны и твердого ядра Земли, выглядит следующим образом:

d2 г

dt

= -От.

d 2 гм dt2

d2 г.

= -От.

г ' Е - Отм г - Е г 'м

1 |3 г Е г - Е — |3 - г м

г 'м - ОтЕ г м -г ' Е

1 |3 1 .3

г м г м - г Е

dt2

— (+ Км + 4 + К + Ка ). т7

(1)

Система дифференциальных уравнений второго порядка (1) сводится к системе уравнений первого порядка путем введения дополнительных векторов

скорости Земли, Луны и ядра - УЕ, Ум , V1 .

dt

=

—^ = -От —— - От

dt \г

г - г

' Е 'м

г - г 'в 'м

м

dt

= V,

г - г

'м ' Е

ШЕ, Я,

= -От ^ - ОтР -

1. ^ К |3 Е I— —

Ш \г г - г

м м Е

Шг1 ШХ

Шг

= V,

- (+ Км + + К + Ка ).

Ш т1 (2) Соотношение (2) представляет систему из 18-ти обыкновенных дифференциальных уравнений

У1 = УЕ + [¥Е , —1 - —Е ]•

(4)

ВЫБОР ЧИСЛЕННОГО МЕТОДА

Для выбора численного метода была исследована матрица Якоби задачи трех тел: Солнце, Земля, Луна:

Ш 2е

Ш % ШХ 2

= -От.

= -От.

г ' Е Отм г - Е г 'м

| ,3 — |3

г Е г - Е - г м

г 'м - ОтЕ г м -г ' Е

1 ,3 1

г м г м - г Е

(5)

Получены приближенные значения собственных чисел матрицы Якоби дифференциальной системы (5). Структура матрицы Якоби и оценки собственных чисел изложены в работе Исаева С.И., Киреев И.В., Новиков Е.А. Оказалось, что среди собственных значений матрицы Якоби системы имеются как вещественные отрицательные и положительные, так и чисто мнимые. Поэтому метод интегрирования должен быть устойчивым в левой полуплоскости комплексной плоскости, неустойчивым в правой полуплоскости, а модуль функции устойчивости должен быть равен единице на мнимой оси. Такой метод будет правильно по крайней мере с качественной точки зрения реализовывать физические последствия каждой из этих групп собственных чисел, в том числе осцилляции и возможную неустойчивость. Этими свойствами обладает известный пятистадийный метод Рунге-Кутты-Мерсона четвертого порядка точности.

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

у = f (х, у), у(х0 ) = у0, х < х < , (6)

где у и f - вещественные N -мерные вектор-

функции, независимая переменная. Метод Мерсо-на численного решения задачи Коши имеет вид

1 2 1

Уп+1 = Уп +7 к + 7 к4 к5' (7)

6 3 6

к1 = ¥ , Уп),

1 1

к2 = ¥ ($я +- И, Уп + - кх),

1 1 1

к3 = ¥(*п +-И Уп +-к1 +-к2Х

3 6 6

1 1 3

кА = И/(^ + - И, Уп +- кх + - кз), 4 2 8 1 8 3

1 3

кз = И/(^ + И,Уп +-к - -к3 + 2к4). 5 п п 2 1 2 3 4

Пятое вычисление правой части правой части задачи (6) не дает увеличение порядка точности до пятого, однако позволяет расширить интервал устойчивости примерно до величины 3.5 по вещественной оси и оценить величину локальной ошибки 8п 4 с помощью стадий к, 1<г<5, то есть

5пА = (2к - 9к3 + 8к4 -кз)/30. (8)

Предложен эффективный способ учета накопления глобальной ошибки из локальных ошибок (Новиков Е.А., 1997). В результате для контроля точности метода Мерсона можно применять неравенство

||2к - 9к3 + 8к4 - к5||< 150£125 (9)

где е - требуемая точность интегрирования.

Перейдем к построению неравенства для контроля устойчивости. Устойчивость численных методов обычно исследуется на линейном уравнении у=Ху, где X интерпретируется как некоторое собственное число матрицы Якоби задачи (6). Одна из причин успеха метода Мерсона заключается в области устойчивости, которая достаточно велика как по вещественной, так и по мнимой оси комплексной плоскости 2=кк. Это позволяет использовать его при решении задач, собственные числа матрицы Якоби которых имеют достаточно большую мнимую часть. Функция устойчивости Q(z) схемы (7) имеет вид

0( z) = 1 + z + z2/2 + z 3/6 + z4/24 + z5/144

Неравенство для контроля устойчивости получим следующим образом. Применяя к разности (к3-к2) формулу Тейлора с остаточным членом в форме Ла-гранжа первого порядка, будем иметь

к3 - к2 = 16И ^дУу ) (к2 - к ),

где вектор £,п вычислен в окрестности решения у(^). Разлагая стадии к1 и к2 в ряды Тейлора по степеням к нетрудно видеть, что имеет место соотношение

к2 - к = к2// / 3 + 0(к3).

С применением степенного метода оценки максимального собственного числа матрицы Якоби системы (6) для контроля устойчивости численной формулы (7) можно использовать неравенство

6 тах

1<г< N

(к3 - к2 )г

(к2 - к )г

< 3.5.

(10)

где числу 3.5 примерно равна длина интервала устойчивости метода Мерсона.

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

Сформулируем алгоритм переменного шага на основе метода (7) с контролем точности вычислений и устойчивости численной схемы. Введем обозначения

Сп = тах1 (2к - 9к3 + 8к4 - кз),. | /(I уп | + г),

¥„ = 6тах|(к -к2), |/|(к2 -к1)1 |.

1<г< N

Тогда для контроля точности схемы (7) можно применять неравенство вида (9), то есть С<£1,25, а для контроля устойчивости (10) - У<Р. Постоянную Р можно выбрать равной числу 3.5, то есть примерно равной длине интервала устойчивости метода (7). В случае выполнения неравенства Уп'|<г по г-й компоненте решения контролируется абсолютная ошибка ге, в противном случае контролируется относительная ошибка е.

Пусть приближенное решение уп в точке 1п вычислено с шагом кп. Тогда алгоритм интегрирования переменного шага на основе схемы (7) с контролем точности и устойчивости формулируется следующим образом.

Шаг 1. Вычисляется стадия к1.

Шаг 2. Вычисляются стадии к., 2<<5, по формулам (7).

Шаг 3. Вычисляется оценка ошибки С .

Шаг 4. Вычисляется значение параметра q по формуле q4 Сп=е1,25.

Шаг 5. Если q<1, и происходит повторное вычисление с шагом qh, то есть происходит переход на шаг 2 (стадия к1 зависит от шага линейно).

Шаг 6. Перевычисляется значение параметра q по формуле q5Cn=e1•25.

Шаг 7. Вычисляется приближение к решению уп+1 по схеме (7).

Шаг 8. Вычисляется значение параметра г по формуле гУ^3.5.

Шаг 9. Вычисляется новый шаг интегрирования йп+1 по формуле

кп+1 = тах(кп, тт(г д)кп). (11)

Шаг 10. Выполняется следующий шаг интегрирования по схеме (7), то есть осуществляется переход на шаг 1.

Отметим, что формула (11) применяется для прогноза величины шага интегрирования к после успешного вычисления решения с предыдущим шагом кп, и поэтому фактически не приводит к увеличе-

нию вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа матрицы Якоби системы (6). Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг кп. В результате для выбора шага и предлагается формула (11), которая позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость.

ЗАКЛЮЧЕНИЕ

Для решения полученной задачи Коши (2), (3), (4) был создан пакет программ Core. В состав комплекса программ включена программная реализация метода Рунге-Кутты-Мерсона с переменным шагом и контролем точности. Программный комплекс реализован на языке программирования С++. Задача характеризуется сильной разномасштабностью, поэтому наряду со стандартными библиотеками iostream, onio.h, math.h, stdio.h была задействована математическая библиотека qd/dd_real.h, позволяющая увеличивать точность математических вычислений. Расчеты проводились с точностью double double - два двойных слова.

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

1. Траектория движения - эллипс, полуоси которого зависят от начального смещения ядра.

2. Расстояние от центра ядра до геометрического центра Земли зависит только от первоначального смещения ядра.

3. Серия расчетов позволила увидеть следующую закономерность - период колебаний внутреннего ядра относительно центра Земли не меняется в зависимости от начального отклонения, и равен 0.076 суток, т.е. 110 минут.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

Исаева С.И., Шайдуров В.В. Математическая модель движения твердого ядра Земли // Известия высших учебных заведений. Поволжский регион. - №1(17). - 2011.

- C. 40-46.

Исаева С.И., Киреев И.В., Новиков Е.А. Сравнение некоторых численных методов на гравитационной задаче трех тел: Земля, Луна, Солнце // Вестник Сибирского государственного аэрокосмического университета -№6(39). - 2011. - C. 20-24. Исаева С.И., Новиков Е.А. Численное сравнение методов Мерсона и Эверхарта на гравитационной задаче двух тел // Системы управления и информационные технологии. - №1(43). - 2011. - C. 25-29. Новиков Е.А. Явные методы для жестких систем. - Новосибирск, Наука, 1997, 197 с. Song X., Richard P.G. Seismological evidence for differential rotation of the Eart's inner core/ Nature. - 1996. - V.382.

- P.221-224

Поступила в редакцию 01 ноября 2012 г. Принято к печати 07 декабря 2012 г.