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

Вычислительные технологии

Том 16, № 2, 2011

Маломодовая модель геодинамо*

Г. М. Водинчар1,2, Л. К. Крутьева1 1 Институт, космофизических исследований и распространения радиоволн ДВО РАН,

п. Паратунка Камчатского края, Россия 2 Камчатский государственный университет им. Витуса Беринга, Петропавловск-Камчатский, Россия e-mail: gvodinchar@yandex. ru, kruteva_lu@mail. ru

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

Ключевые слова: конвекция, тороидальные и полоидальные поля, ядро Земли, геодинамо.

Введение

Исследование вопросов формирования геомагнитного поля является одним из наиболее интенсивно развивающихся направлений геофизики, тем более, что работы в этой области имеют пересечение с проблемами космического магнетизма, задачами коемо-и астрофизики. Обзоры постоянно увеличивающегося числа исследований данного направления приведены, например, в [1, 2]. В настоящее время практически общепризнанной для геомагнетизма и космического магнетизма является теория динамо. В этой теории достигнут огромный прогресс, однако нельзя считать, что задача формирования и поддержания геомагнитного поля полностью решена. На сегодня нет модели, которая объясняла бы все наблюдаемые свойства поля.

Одним из ключевых для геодинамо является вопрос о структуре конвективных течений в жидком ядре. Косвенную информацию об этой структуре можно получить из данных о неоднородноетях в плотности жидкого ядра. В [3] проанализированы результаты ряда работ по врНШг^-функциям собственных колебаний Земли, в которых получены срезы распределения плотности на различных глубинах. Вариации плотности на глубине 3900 км, соответствующие врНШг^-функции жидкого ядра, приведенные в [3], представлены на рис. 1. Здесь прослеживается четкая 12-зонная шахматная структура. Автором работы [3] на основе этих данных была высказана гипотеза о соответствующей структуре конвекции, где в шести областях материал ядра "тонет", а в шести — "всплывает".

*Работа выполнена при финансовой поддержке ДВО РАН (проект 10-Ш-В-07-158).

90n 6(Ш

з<ж о

ЗОБ 608 908

о 60е 120е 180 120\¥ 60 w о

Рис. 1. Портрет 8рНШ^-функции для моды ц£4 собственных колебаний Земли из работы [3]. Черный цвет — плотность вещества на 0.2% выше средней, белый — на 0.2% ниже средней. По горизонтальной оси — градусы долготы, по вертикальной — широты

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

1. Формулировка краевой задачи геодинамо

Рассмотрим вращающуюся вместе с Землей с угловой скоростью Q систему координат, начало которой расположено в центре Земли, а ось Oz проходит через Северный полюс. Обозначим через v, В и P поля скорости, магнитной индукции и давления. Поле температуры внешнего ядра представим в виде Т + Ts, где Ts — стационарное распределение температуры, соответствующее теплопередаче в виде чистой теплопроводности, T

Будем использовать следующие упрощающие предположения: вещество внешнего ядра несжимаемое, относительная магнитная проницаемость всего пространства / = 1, вариации плотности внешнего ядра относительно среднего значения р0 малы, кинематическая вязкость v и температуропроводность k внешнего ядра, а также магнитная вязкость vm всего ядра постоянны. Среду за пределами ядра считаем непроводящей,

В

Т\ и внешней т2 = т\ + h границах жидкого ядра сохраняет постоянные значения Т и Т2 = T — ST. Эти предположения являются обычными при постановке задач геодина-

v

т. е. в рассматриваемой модели отсутствует явление супервращения внутреннего ядра. Тогда уравнения динамо записываются в приближении Буссинеска в виде [1]

d v 1 g2 1

— + (vy)v = i/Av--V Р + (3—Т Г -2(ilxv) +-rot В х В,

dt ро Т 2 /0р0

дТ д В

— + (vV)(T + Ts) = kAT, — = rot (v x В) + um Д В,

▽v = 0, уВ = 0. (1)

Здесь g2 _ ускорение свободного падения на границе ядра, в _ коэффициент объемного теплового расширения внешнего ядра, _ магнитная постоянная.

Известно, что стационарное решение уравнения теплопроводности в сферической оболочке с граничными условиями Ts(r = ri) = и Ts(r = r2) = T2 имеет гиперболи/ ri \

ческий no pa. uivcv профиль Ts = —-— (--1) + T\.

h V r /

Если в качестве единиц измерения расстояния, скорости, времени, давления, температуры и магнитной индукции принять величины h, v/h, h2/и, p^v2/h2, 5T и v^/iopo/h соответственно, то уравнения (1) в безразмерных переменных запишутся в следующем виде (обозначения переменных сохранены):

dv Tr

— + (vy)v = Av - \/P + RaPr"1—er - r (ez x v) + rotB x B, dt r2

dT v dB

— + (vy)T - n r24 = Pr"1 AT, — = rot (v x B) + q-lYi~l Д B, dt r2 dt

Vv = 0, yB = 0. (2)

Управляющими параметрами модели являются: число Релея Ra = STg2h3в/(vk), число Прапдтля Pr = v/k, Кориолиса т = 2h2Q/v, Робертса q = k/vm.

Для исключения поля давления возьмем ротор от обеих частей первого уравнения системы (2), учитывая, что (vy)v = (1/2)gradv2 — v х rot v. Получим

dv Ra /Tr \ .

rot—--rot (v x rotv) = rot Л v + —rot —er — r rot (ez x v) + rot (rotB x B),

dt Pr \r2 J

Г\ГТ1 ЛТЭ

— + (vy)T - пг2Ц = Pr"1 AT, — = rot (v x B) + q-'Pг"1 Д B, dt r2 dt

Vv = 0, VB = 0.

Систему (1) дополним однородными граничными условиями для температуры и условиями прилипания для скорости на внутренней и внешней границах жидкого ядра: T (r = r1) = T (r = r2) = 0, v(r = r1) = v(r = r2) = 0,

Рассмотрим теперь граничные условия для магнитного поля. Поскольку при r > r2 поле потенциально, то в этой области B = —gradU, где потенциал раскладывается в ряд по сферическим гармоникам

/ r \ -(П+1) n

[/ = r2E(fj Е AnWn(e^)- (4)

n=1 ^ ' m=-n

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

n / / \ -(n+1) \

В = ЕЕ A-^TOtTOt[ri[72) Ynm(e^> . (5)

n=1 m=-n \ J

В справедливости этого разложения можно убедиться непосредственным вычислением двойного ротора в формуле (5) и градиента выражения (4), Соленоидальное поле B

деляютея соответственно как гс^(Фг) и rot гс^(Фг), где Ф и Ф — некоторые скалярные

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

( ( \ -(ra+1) \ ВГ™ = rot rot ( ^ J Y?(d,<p)г J . (6)

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

ВПт = rot (Rm(T,t)Yram(^)r) И В^ = rot rot (R^T,t) 1^(9,^) . (7)

Вид функций R^m(т, t) и Rpm(r,t) конкретизируем позже. Тогда для обеспечения непрерывности магнитного поля при переходе через границу ядра т = т2 краевые условия примут следующий вид:

ВПт( Т = Т2) = 0, rot В^т( Т = Т2) = О, ВРт( Т = Т2Ж( Т = Т2), rot В^( Т = Т2) = 0. (8)

2. Спектральное разложение полей

Для построения маломодовой модели геодинамо будем раскладывать поля температуры, скорости и индукции по собственным полям спектральных задач, связанных с оператором Лапласа, Температуру представим в виде Т = ^^ kanm(t) k@nm(r, 9, p), где

k,n,m

k@nm _ собственные функции оператора Лапласа, пулевые при r = т1,2. Тороидальную составляющую поля скорости запишем в виде v = ^^ keTm(t) kvTm(r, 9, p), где kv'Tm =

( ) k,n,m

rot (Rfcn(r)Ynm(9, p)r) — собственные поля векторного оператора Лапласа, удовлетворяющие условию прилипания при r = т1,2. Наконец, полоидальную часть скорости представим в виде vpm = Y^ kAL(t) kvpm(г-Д^, где kvpm = rotrot (Rpn(r)Ynm(9, P)r) -

k,n,m

собственные поля спектральной задачи rot △ P + / rotP = 0 в пространстве полоидальных полей, нулевых при r = т1,2.

Строение функций k Onm описано в [6], построение полей kv^m, kvpm выполнено в [4]. Компоненты (7) магнитного поля также будем раскладывать по полям задачи

rot △ Вnm + V rotВnm = 0 (9)

с соответствующими краевыми условиями. Верхний индекс из формулы (7) опускается, поскольку пока рассуждения носят общий характер для полей обоих типов. Разделяя радиальную и временную переменные, представим поля В^ в виде

Вnm = Y kYnm(t) rot (kRnm(0^(9, p)r) = kYnm(t) kВ„т. kk

Здесь индекс k = 0,1, 2,... соответствует дискретизации спектра задачи (9) по радиальной переменной, В [4] показано, что функции kRnm с учетом их ограниченности

в центре Земли имеют вид (у^г) + У™, где ]п{') ~ сферические функ-

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

Далее рассмотрим поля разных типов отдельно. Краевые условия для тороидальных пт дают систему уравнений

нолей кВТ

^«т (г2)

и

ак Дпт

Иг

0

(10)

Акп Вкп

системы определяет для каждого п счетное множество собственных значений пТп как решений уравнения

Зп {ф]ыГ2) ПГ2 - г2

Фп

Иг

0.

Г=Г 2

После нахождения коэффициентов Пы подставляем их в одно из уравнений снсте-

Акп Вкп

множителя. Условие нормировки, соответствующее единичной безразмерной энергии, выделяемой полем кВ^т во всем пространстве, выберем в виде

у (кВ4тГ ИУ = ] (кВптГ ИУ = 1.

В3 Т<Т2

Для полоидальных полей краевое условие го1 Врт(г = г2) = 0 дает уравнение

И2 2 И п(п + 1)'

Иг2 г Иг

дВР\ Дкп 1

кп |г=Г2

0,

а условие Врт(г = г2)

Вптт(г = г2) приводит к формуле

А

Иг

п + 1

кп |г=Г2

0.

(Н)

(12)

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

(кВРт)2 ИУ

г>ВР Дкп

п2(п + 1)2 Иг +

В3

ВР Дкп

ИДкп

Иг

г2п(п + 1) Иг = 1. (13)

Как и в предыдущем случае, это условие приводит к единичной энергии, выделяемой полем во всем пространстве, В левом интеграле формулы (13) интегрирование ведется по всему пространству, а правая часть формулы представляет собой результат аналитического интегрирования левой части по угловым переменным. При интегрировании по промежутку г > г2 в (13) предполагается, что в этом промежутке

ДВРР (г)

ДьГ Ы ( -

кп г2

-(п+1)

Такое выражение накладывается непрерывным переходом Врт в ВПт» Все вышеописанные расчеты базисных магнитных мод, связанные с решением уравнений на собственные значения, систем для коэффициентов, нормировками, проводились в пакете МАРЬЕ 12. При этом интегрирование по угловым переменным выполнялось аналитически, а по радиальной — численно. Эти расчеты были проведены для к = 0,1, 2 и п = 1,..., 10.

3. Построение маломодовой модели

Выполним отбор мод скорости, которые определят описанную во введении структуру течений с 12 чередующимися зонами поднятия и опускания вещества. В работе [4] показано, что подобная структура вертикальных течений описывается полоидальными компонентами кvp±2. При этом крупномасштабная конвекция с транспортом материа-

к=0

компоненты 0и некоторые ее линии тока приведены на рис. 2. Линейными комбинациями двух таких мод можно обеспечить любой фазовый сдвиг 12-зонной картины по углу р. Поскольку выбор начала отсчета угла р произволен, можно ограничиться только модой 0 2.

Чтобы учесть кориолисов снос основной компоненты скорости, направление еъ х 0действующей на нее силы Кориолиса аппроксимировалось другими компонентами скорости. В качестве критерия приближения использовалась минимизация невязки

£

к , п , т

(кЯ.пт кVnm

+ кёпт кУпт) ег Х

0 у4 , 2

в метрике скалярного произведения (Р, Р) = J (Р Р) ёУ, где интегрирование ведется по объему жидкого ядра. Результаты расчетов показали, что отличными от нуля являются только коэффициенты кд32, кд5,2 и 0в4 -2. Представим выражение е2 х к в порядке убывания коэффициентов:

е2 х кv4P2 » 0.410vn2- 0.34!vп;2- 0.250Vзп2- 0.17!- 0.10<_2- 0.062уП,2 + 0.042у3т2 + ...

180

Рис. 2. Линии тока моды 0V442 (а) и ее радиальная компонента (б). Черный цвет — течение снизу вверх, белый — сверху вниз

С учетом этого разложения для аппроксимации кориолнсова сноса основной моды скорости в модели использовались две тороидальные компоненты и ^Т2-

Для представления температуры были оставлены две моды: 1в00 и 0©42. Первая дает равномерное отклонение по радиусу от стационарного профиля (используется по аналогии с моделью Лоренца маломодовой конвекции в плоском слое [7]), вторая "запускает" основную конвективную моду ^р2.

Магнитное поле представим модами 0Вр0, 0Вр±1, описывающими дипольную часть,

а также пространственно (по индексам) связанными с компонентами скорости модами ВТ ВР ВТ ВР ВТ ВР

0В5,±2> 0В5,±2> 0В3 ,±2> 0В3 ,±2> 0В4 ,±2> 0В4,±2"

Таким образом, первоначально в модели будем использовать три компоненты скорости, две — температуры, 15 — магнитной индукции, В дальнейшем число магнитных мод можно сократить.

Для удобства в соответствии с табл. 1 перейдем к одноиндекеным обозначениям. Итак, принимаются разложения полей:

2 14

Т = «0(^)©0 + «1(^)©1, V = В = 7^)В». (14)

г=0 г=0

Собственные значения мод температуры, скорости и индукции обозначим соответственно через Аг, пг.

Таблица 1. Одноиндексные обозначения

Трехиндексные комбинации и типы полей Один

к п т Тип поля индекс

Моды температуры

1 0 0 0

0 4 2 1

Моды скорости

1 3 2 1юг 0

0 4 2 ро1 1

0 5 2 1юг 2

Моды индукции

0 1 -1 ро1 0

0 1 0 ро1 1

0 1 1 ро1 2

0 4 -2 1юг 3

0 4 2 1юг 4

0 4 -2 ро1 5

0 4 2 ро1 6

1 3 -2 1юг 7

1 3 2 1юг 8

1 3 -2 ро1 9

1 3 2 ро1 10

0 5 -2 1юг 11

0 5 2 1юг 12

1 5 -2 ро1 13

0 5 2 ро1 14

Разложения (14) подставим в систему (1), предварительно взяв в ней ротор третьего уравнения. Следуя идее метода Галеркипа [8], получим квадратичную систему обыкновенных дифференциальных уравнений с постоянными коэффициентами для амплитуд

¿Ев + Е LkjYiYj, k = 0,1, 2,

2 3 3 1

Е А^ = Е - Е + RaPr"1 Е С-«Г

i=0 i,j=0 i=0 j=0 2 14

^k q | \ л г k

i=0 i,j=0

d 2>1 2

-Jr = E + E - Pr_1^> «= o, i,

i,j=0 i=0

14 , 2,14 14

lr E ^ = E - 9"lpr_1 E / = 0,..., 14. (15)

i=0 i,j=0 i=0

Здесь Ak = (rotvfc, rotvi), Bk- = (rotvfc, rot (vi x rotvj)), Cj = (rotvfc, rot (Tjrer/f2)), Ej0 = -(rotvfc,rot (e* x vi)) Lk- = (rotvfc,rot (rotBi x Bj)), Fj = -(©?, (viV) @j) H? = (6?, (vi)r(r1r2/r2)) Qi = (rotBi, rotBi) Wj = (rotBi, rot rot (vi x Вг)). Все эти скалярные произведения являются интегралами по объему либо всего ядра (коэффициенты Qi), либо жидкого ядра (все остальные коэффициенты) от известных базисных полей. Они были вычислены в системе MAPLE 12, причем аналитическое интегрирование по в и показало, что многие коэффициенты пулевые, В частности, равны нулю все коэффициенты Wj для l = 0, 2, 5,..., 8,11,12, а матрицы с элеме нтами Ak и Qi — диагон^ьные. Тогда из системы (15) видно, что амплитуды магнитных мод Bi

l

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

)d£о dt

= -А'фо/Зо + тЕ^г + l7l74 + (L% + J ЪЪ,

= -A\fiifji + ^Ciai + г {El?о + ¿ЭД + {b{3 + Ь\л)

+ (L1,10 + L10,0 71710 + (L1,14 + ^4,1) 71714,

^li = + rE\/3, + L\l7l74 + (L2;13 + L\3>1) 7l7l3,

dQ!° =F1°1/31Q!l + Pr-1AoQ!o,

dt

= ^10^1^0 + - Pr-1A1«1,

dt

= <9/3o79 + <3^73 + winatio + <14/51714+

+W2,13^2713 - q-1Pr-1Q1n1Y1,

Ql^ = <ЛЪ ~ q-'Pr-'Ql^, Qu^f = w^jMi -

Qlt^f = Wl- g-1Pr"1Qi^i47i4. (16)

Примем следующие значения для параметров земного ядра [9]: v = 10-6 м2/е, vm = 10-6 м2/с, k = 10-5 м2/с, iT = 103 К, в = 10-4 К-1, g2 = 7 м/с2. Известно также, что h = 2.1 ■ 106 м, tt = 7.3 ■ 10-6 рад/с, ri = 1391 км.

Система (16) при любых значениях управляющих параметров имеет нулевую точку покоя, соответствующую отсутствию конвекции и магнитного поля. При вычислении ненулевых точек покоя учтем, что система (2) обладает симметрией относительно замены знака магнитного поля, а система (16) обладает симметрией относительно смены знаков амплитуд скоростных мод и температурной амплитуды а1 при сохранении знака а0, Эта вторая симметрия специфична для рассматриваемой маломодовой модели.

Численно, с использованием пакета MAPLE 12, были найдены три несимметричные ненулевые точки покоя. Координаты в1 и 71 этих точек, определяющих характерные скорость конвекции и интенсивность основного диполя, приведены в табл. 2. С учетом вышеуказанных снмметрпй первая точка таблицы расщеляется на две, а каждая из двух оставшихся — на четыре. Первая точка соответствует непроводящему материалу ядра, т. е. в контексте настоящей работы интереса не представляет.

Пересчитывая безразмерную скорость в систему Си, получим, что в1 = 2.46 ■ 108 ~

10-4 м/с, в1 = 2.77 ■ 107 ~ 10-5 м/с. Имеющиеся оценки реальной характерной скорости

10-4

Вне ядра дипольной моде B1 = 0Bp0 соответствует безразмерная компонента потенциала Rqi(^2) ( — ) которая на поверхности Земли имеет вид 2.93- Ю-2 cos 9. Тогда

r2 10 стационарным амплитудам будут соответствовать значения потенциала 2.01 ■ 1010 cos 9

и 1.76 ■ 1011 cos 9. Аналогичная безразмерная компонента потенциала в модели IGR.F [10]

для поверхности Земли равна в безразмерном виде 1.53 ■ 109 cos 9.

Таким образом, точки покоя с координатами в1 = ±2.46 ■ 108, y1 = -6.87 ■ 1011 дают стационарные решения модели, совпадающие по порядку величин с имеющимися оценками скорости конвекции и на порядок отличающиеся от наблюдаемой величины основного диполя. Учитывая большую неопределенность в оценках таких параметров ядра как коэффициенты температуропроводности, объемного расширения и особенно

Таблица2

Номер точки покоя а 71

1 1.54 • 102 0

2 2.46 • 108 6.87- 1011

3 2.77 • 107 6.02 • 1012

вязкости, можно говорить, что эти точки покоя дают стационарные режимы, близкие к наблюдаемым.

Заключение

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

Список литературы

fi] Kono M., Roberts P.H. Recent geodvnamo simulations and observations of the field // Rev. Geophysics. 2002. Vol. 40, No. 10. P. m B!l.

[2] Jones С.A. Convection-driven geodvnamo models // Phil. Trans. R. Soc. Lond. A. 2000. Vol. 358. P. 873-897.

[3] Кузнецов B.B. Анизотропия свойств внутреннего ядра Земли // Успехи физ. наук. 1997. Т. 169, № 9. С. 1001-1012.

[4] Водинчар Г.М., Шевцов Б.М. Маломодовая модель конвекции во вращающемся шаровом слое вязкой жидкости // Вычисл. технологии. 2009. Т. 14, № 4. С. 3-15.

[5] Голицын Г.С. Режимы конвекции на различных вращающихся геофизических и астрофизических объектах // Изв. АН СССР. Физика атмосферы и океана. 1991. Т. 27, № 1. С. 20-31.

[6] Тихонов А.Н., Самарский A.A. Уравнения математической физики. М.: Наука, 1977. 735 с.

[7] Заславский Г.М., Сагдеев Р.З. Введение в нелинейную физику: От маятника до турбулентности и хаоса. М.: Наука, 1988. 368 с.

[8] Ладыженская O.A. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 232 с.

[91 Merril R.T., McElhinny M.W., McFadden P.L. The Magnetic Field of the Earth. N.Y.: Acad. Press, 1996. 532 p.

[10] International Geomagnetic Reference Field.

http ://www.ngdc.noaa.gov/IAGA/vmod/igrf.html

Поступила в редакцию 5 апреля 2010 г., с доработки — 22 апреля 2010 г.