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

МАТЕМАТИЧЕСКАЯ ФИЗИКА, МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

MS С 76W05

ШЕСТИСТРУЙНАЯ КИНЕМАТИЧЕСКАЯ МОДЕЛЬ ГЕОДИНАМО

Г.М. Водинчар, Л.К. Фещенко

Институт космофизических исследований и распространения радиоволн ДВО РАН, ул. Мирная, 7, Паратунка, Елизово, 684034, Россия, e-mail: ikir@ikir.ru

Аннотация. Построена кинематическая модель геодинамо, структура нолей скоростей которой согласована с данными о splitting-функциях собственных колебаний в жидком ядре Земли. В модели используются три компоненты скорости, являющиеся аппроксимациями мод собственных колебаний жидкох'о ядра, две из которых определяют 6-етруйную конвекцию в ядре Земли, а третья соответствует опережающему вращению твердoi'o ядра.

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

Введение. Процесс формирования магнитных полой планет и звезд на качественном уровне успешно объясняется теорией гидромагнитного динамо. Разработанные модели конвекции в жидких ядрах планет земного типа, газовых гигантах, конвективных зонах звезд позволяют получать течения, которые могут формировать магнитные ноля, близкие но своей топологии к наблюдаемым |1|- |4|,

Возможности вычислительных систем не позволяют вести прямое численное моделирование трехмерных задач планетарного динамо на геологических временных масштабах. Отметим, что известные теоремы запрета определяют принципиальную трехмерность задачи динамо |2|, В связи с этим численные модели либо воспроизводят МГД-течения с хорошим разрешением но пространству на относительно небольших временных масштабах, порядка десятков тысяч .нет, либо дают возможность просчитывать длительную эволюцию только крупномасштабных пространственных структур. Дня моделей первого тина геометрическая структура течений просчитывается в процессе моделирования, а для моделей второго типа геометрическую крупномасштабную структуру конвекции надо задавать. Возникает ключевой вопрос о том, какова реальная крупномасштабная структура конвекции.

Косвенную информацию об этой структуре можно получить из данных о неоднородностях в плотности жидкого ядра. В статье |5| были проанализированы результаты ряда работ но эрПШ^-функциям собственных колебаний Земли и получены срезы распределения плотности на различных глубинах. Вариации плотности соответствующие эрПШ^-функции п£4, имеющей максимум на глубинах жидкого ядра, представлены па рис. 1. Здесь прослеживается четкая 12-зониая шахматная структура, которой в первом приближении соответствует тессеральная сферическая гармоника У42. Автором [5]

Работа выполнена при поддержке ДВО РАН (проект 10-III-B-07-158) и Минобрнауки России по Программе стратегического развития КамГУ им. Витуса Беринга на 2012-2016 гг.

была высказана гипотеза о соответствующей структуре конвекции, где в шести областях материал ядра «тонет», а в шести - «всплывает». Возможность существования подобной конвективной структуры и генерации в пей магнитного ноля диполыюго тина рассматривалась в работах авторов |6|- |7|. Были найдены квазистационарные решения в которых величина характерной скорости конвекции составляла ~ 10-4 м/с, совпадая с известными оценками дня ядра Земли |8|. При этом величина дииольпой составляющей модельного ноля в пересчете па гауссовский гармонический коэффициент была ~ 105 нТл, при реальном значении в 3 • 104 нТл [9].

Описание геометрии течений только гармоникой У42 не учитывает некоторую «скошенность» распределения па рис. 1. В настоящей работе исследуется более точное воспроизведение наблюдаемых конвективных структур в рамках кинематической модели.

9(Ж 6(^ зоы о

ЗОБ 608 908

0 60Е 120Е 180 120\\^60\¥ 0

Рис. 1. Портрет splitting-фyнкции для моды 11 £4 собственных колебаний Земли из работы [5].

Черный цвет плотность вещества на 0.2% выше средней, белый плотность на 0.2% ниже средней. По горизонтальной оси отмечены градусы долх'оты, но вертикальной широты.

1. Выбор структуры течений. Изображенное па рис. 1 распределение раскладывалась по сферическим гармоникам У™ до п = 6. Вычислялись коэффициенты разложения Сп'±т и амплитуды Ап^т = ^с2 т + с^г _т. В таблице 1 приведены результаты разложения, упорядоченные но убыванию амплитуд. Ряд гармоник оборван но признаку резкого падения амплитуды Л4>4.

Таблица 1

п т С-п,т С>г,—т л ' 11/.П/

4 2 0.059 -0.073 0.094

2 0 0.042 0.042

4 0 0.032 0.03

4 3 0.011 -0.02 0.023

2 2 -0.016 0.006 0.017

4 1 -0.01 -0.013 0.016

2 1 -0.003 0.016 0.016

4 4 -0.003 -0.001 0.003

Дня анализа содержания таблицы 1 рассмотрим сначала представление ноля скорости в задаче геодинамо, использовавшееся авторами в |7|, Все построения ведутся в геоцентрической системе координат, где жидкому ядру соответствует сферическая оболочка г1 < г < г2, а за единицу длины принята толщина жидкого ядра. Тогда г1 = 0.664 и г2 = 1.664.

Материал ядра считается несжимаемым и па границах контактирующим с твердым телом. Тогда Vv = 0 и скорость нулевая на границах.

Соленоидальное поле скорости V раскладывается в сумму тороидальной V71 и поло-идальной V5, составляющих, для которых в свою очередь используется разложения

где 9 и ф - коширота и долгота, соответственно.

Базисные поля V)7 п т(г, 9, ф) и V)7 п т(г, 9, ф) являются собственными модами затухания дня уравнения Навье-Стокса, т.е. являются решениями спектральной задачи

^ + Av — Vp = 0

в пространстве соленоидальных полей, нулевых при г = г1 , 2. В равносильной формулировке эта задача распадается на задачи

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

Сами базисные поля имеют вид

где функции ЯТп(г) и Я^п(г) определяются го задач (2), а собственные значения ^[п и ^^п входят в выражения для этих функций. Уравнения на собственные значения и схема расчета функций ЯТп(г) и Я^п(г) описаны в [7].

Система собственных функций {V)7 п т, V5п т} обладает свойством ортогональности

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

Наблюдаемая в эрПШ^-функциях неоднородность в распределении плотности по сфере соответствует сферическому распределению вертикального переноса вещества ядра, т.е. радиальной составляющей ноля скорости. Тороидальные ноля не имеют радиальной составляющей и не определяют конвекцию как таковую. Собственно они возникают в конвекции в результате корио.нисова сноса нолондальпой части скорости. Соответственно, тороидальные компоненты нельзя увидеть в рис. 1. Радиальная нроек-

„ 5 п(п + 1) , \-trm/п \

ЦИЯ полоидалыюи МОДЫ пт имеет ВИД -----------Пкп{Г)Уп (6, ф), что позволяет связать

данные табл. 1 с вкладом различных V5 п т в поле скорости.

(1)

+ AvT = 0, + rotAvS = 0

(2)

'^'Т, п,т = ^ (ДТп(г)Упт(9, ф)г) > ^1,п,т = Г°^ {^кп(г)УГ(9, ф)Г) > (3)

и полноты относительно скалярного произведения

(4)

Учтем также следующие обстоятельства. В классе еолеиоидальиых нолей линейные оболочки Нт собственных ПОлей V7 п ±т, V5 п ±т устойчивы относительно кориолисова сноса |10|. Точнее говоря, вращение приводит к расщеплению, например, нолоидаль-ной моды V5 п т на бесконечную цепочку чередующихся полоидальных и тороидальных мод вида vfn _т, vTn±1 ±т, vfn±2 ±т, vTn±3 ±т, ... [7]. Аналогична схема расщепления тороидальных мод.

В связи с этим в табл. 1 строки с индексами (4,2) и (2,2) соответствуют одной устойчивой группе мод. Аналогично можно сгруппировать (2,0) и (4,0), (2,1) и (4,1).

На рис. 2 приведены построенные но данным табл. 1 распределения с накоплением числа используемых групп мод. Видно, что различия между рис. 2а и рис. 2Ь практически пет, т.е. группу содержащую индексы (2,1) и (4,1) можно не добавлять.

а

ь

с

С1

Рис. 2. Распределения комбинаций функций из таблицы 1. а - группа У4±2 и У2±2:

Ь - добавлены к предыдущему У20 и У40; с - добавлены к предыдущему У4±3: с! - добавлены к предыдущему У±* и У2±1.

Полоидальные моды, определяемые сферическим гармониками У±2 и У4±3 задают каждая конвективную структуру в ядре, содержащую 6 ячеек. Они интересны тем, что при существующей относительной толщине земного ядра в 6-ячейковой (6-струйпой) конвекции горизонтальные и вертикальные масштабы ячеек наиболее близки, а такие структуры устойчивы. Отметим, что 6-ячейковая конвекция не реализуется никакими другими комбинациями сферических индексов.

Всего предлагается использовать три группы мод скорости. Две группы соответствуют 6-струйной конвекции. Первая включает в себя полоидальную V5 4 ±2, задающую 6-струйную конвективную структуру, И МОДЫ V5 2 ±2, V7 3 ±2, V7 5 ±2 возникающие

при ее расщеплении за счет вращения. Вторая содержит основную, также задающую 6-струйную конвекцию, моду V5 4 ±3 и ее расщепления V7 3 ±3, vT5 ±3 .

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

Третья группа мод включает полоидальпые составляющие со сферическими индексами (2,0) и (4,0), и ее появление можно объяснить эффектом опережающего вращения твердого ядра. Имеется ряд данных о том, что твердое ядро Земли имеет угловую скорость вращения несколько больше, чем Земля в целом, что объясняют сохранением углового момента Земли при конвективных движениях 1121. В используемом представлении дня поля скорости этот эффект можно учесть,если поменять в выражении (3) тороидальной моды V71 0 радиальную функцию ЯКг) на некоторую убывающую функцию Я(г), удовлетворяющую условиям Я(г1) = 1, Я(г2) = 0. А для того, чтобы эта особая мода скорости не была подвержена вязкой диссипации надо положить

2 2 3

г гг

Я(г) = аг + Ъ/г2, где а = —^—-—3 и Ъ = —3 1 2 3 . Тогда, с учетом (1), варьируя величи-

г1 — г3 г1 — г3

ну соответствующей амплитуды в71 0 можно задать нужную величину опережающего вращения. Кориолисово расщепление этой моды и приводит к возникновению нолои-дальпых мод со сферическими индексами (2,0) и (4,0), а также тороидальной (3,0). Поэтому в состав третьей группы моды входят модифицированная vT1 0, а также V5 2 0,

М , 3 , 0; у0 , 4 , 0-

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

Известно решение задачи Пуанкаре о собственных колебаниях вращающейся идеальной жидкости гAv + 2к х V + Ур = 0, удовлетворяющей условию непроницания на твердой границе [10]. Здесь к - орт оси вращения. Для вязкой жидкости эта задача переопределена, ввиду увеличения чисна граничных условий.

Аналогом задачи Пуанкаре дня вязкой жидкости является спектральная задача Av + бAv + 2к х V + Ур = 0 ОДе е - число Экмана, характеризующее отношение периода вращения ко времени вязкой диссипации. Будем называть ее далее вязкой задачей Пуанкаре.

Эта задача изучалась в работе 1111 дня случая сферической оболочки. Установлены такие важные важные свойства как дискретность спектра, полнота системы собственных функций, получены оценки границ спектра. Однако, точное решение этой задачи по-видимому неизвестно.

Поскольку система полей {V)7 п т, V5п т| полна, можно аппроксимировать вязкие моды Пуанкаре с помощью этих нолей. В частности, дня мод образующих рассмотренные выше три группы можно определить коэффициенты так, чтобы соответствующая группа паилучшим образом аппроксимирована вязкую Пуанкаре моду. Такая аппроксимация в метрике скалярного произведения (4) была построена.

Соответствующие аппроксимации имеют следующий вид:

V! = 1 ,0) — 0.13(8ІП'01^’ 2 ,0) + 0- 1659(^08^ vTз)o) — 0.1374(8^0^^ ,0),

V2 = 0.1577(008^^ 2 , _2 + 2,2) + 0.2575(8^^^ 3 , _2 — со8"02^ 3 , 2) +

+0.5323(cОS^2VS) 4 , _2 + 8ІП^2VS)4,2) + 0.4254( —8ІП^2vT 5 , _2 + C08^2vT 5 , 2), (5)

Vз = 0.3649(8Іп^Т 3 , _3 + сов^з^д 3) + 0.4914(со8^^4 , _3 — 8Іп^^4 , 3) —

— 0.3542(8ІП^Т 5 , _3 + со8^Т 5 , 3),

где ^ - произвольные углы.

Окончательно, скорость в модели принимается в виде

V = а (^іVі + V2 + А^з), (6)

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

2. Кинематическая модель геодинамо. В кинематическом приближении задача динамо сводится к решению уравнения индукции дня магнитного поля

дв

— = го! (V х В) + ЇІе^ДВ, ^

▽в = 0.

при заданной скорости V. Уравнение записано в безразмерном виде. За единицу времени принята величина к/"У0, оде к = 2.1 • 106 м - толщина, жидкого ядра, а «0 = 10_4 м/с - характерная величина скорости конвекции. Основным параметром задачи является магнитное число Рейнольдса И,ет = к^0/^т, оде ^т - магнитная вязкость ядра. Оно дает отношение характерного времени диссипации магнитного поля к2/^т ко времени конвективного цикла к/зд ,

К уравнениям (7) добавляются вакуумные граничные условия дня магнитного ноля 1131. Физически эти условия означают, что среда за пределами ядра непроводящая и токи в пей отсутствуют. Такой допуск оправдан тем, что временные масштабы процессов в ядре составляют от 10з лет и более, т.е. на порядки превосходят характерные времена атмосферных и магпитосферпых процессов. Математически вакуумные условия требуют непрерывного перехода поля В в потенциальное поле при г = г2.

Учитывая описанное выше сиектралыюе представление дня поля скорости, ноле в также представим в виде комбинации собственных тороидальных и нолоидальпых полей оператора Лапласа

В = X! ,тСО В1п,т(Г ^ ,тСО В,п, т(Г ^ ^ (8)

к ,п, т к ,п , т

Проведем усечение рядов (8), оставим в них три составляющие динольпой части геомагнитного поля В^ 1 _1, В^ 1 0 и В^ 1 1; а также моды которые структурно связаны

с компонентами скорости, т.е. «зацепляются» с ними в слагаемом rot (v х B) и регенерируют поло. Если использовать одпоипдекспые обозначения для магнитных мод, то интеграл по объему ядра

дает генерацию магнитной моды Вк при нелинейном взаимодействии компоненты скорости V*, заданной одним из выражений (5) и магнитной В.,-, т.е. собственно работу механизма динамо. При этом в рассматриваемой простой модели пас интересовала прежде всего регенерация динольпых компонент.

Расчет интегралов Ж*к для различных комбинаций мод показал, что для этого надо оставить следующие магнитные моды (всего 17): В7 2 0, В7 2 ±1,В^ 0, В7 4 ±2,В^ ±3, В5 3 0,В5 3 ±1, В5 3 ±2,В5 5 ±2, В5 5 ±3, Всего, вместе с дипольными, получается 20 магнитных мод. Отметим, что поскольку эти моды отличаются какими-либо сферическим индексами, то они ортогональны па сфере, а значит и в объеме ядра.

Далее дня удобства будем использовать одпоипдекспые обозначения дня магнитных компонент.

Представляя магнитную индукцию линейными комбинациями вышеуказанных В.,-, а скорость - разложением (6), получаем галеркипские приближения дня уравнения и и, лук ции (7)

Это система линейных уравнений с постоянными коэффициентами дня амплитуд магнитных мод и изменение величины магнитного ноля со временем тогда можно охарактеризовать числом Л = тахШеА*, оде А* собственное значение матрицы системы (9),

*

Это число определяет скорость изменения величины магнитного ноля.

Произведение ЛИ,ет можно интерпретировать как отношение двух времен - характерного времени диссипации магнитного ноля и характерного времени изменения ноля в рассматриваемой модели. В связи с этим предлагается следующий критерий воспроизведения магнитного поля динамо-механизмом в модели - Л > 0 или |Л|Ие ^ 1. В нервом случае поле нарастает, во втором затухает, по скорость его затухания много меньше скорости омического затухания в отсутствие работы динамо.

Ясно, что значение Л зависит от числа И,ет и от весовых коэффициентов к1; к2 в формулах (6), При этом вычислительные эксперименты показали, что собственные значения матрицы системы (9) не зависят от значений углов -0* в (5),

Была проведена серия расчетов с перебором значений к* от 0.1 до 10 и И.ет от 5 до 500 в логарифмической шкале. Для каждой такой комбинации определялись соответствующие Л и |Л|И,е,

l = 1,..., 20

(9)

собственное значение магнитной

моды Bj.

Re m=12.56

10

ГО -I

v: I

0.1

10

■ f + + + + ■¥■ + +

0.1 10

m л I

0.1

ki

Re m=50

10

Re m=31.55

ПП л

I

0.1

- Ч + + + ' + + + + + ■ + + + + + + + + + + 4- + + + + + + + + + + + + + ;

i + + + • t: + + + + + +- + + + -i- + + + *■ *■ , * +

o!i 1 10

ki

10 0.1

1

ki

Re m=125.6

10

Рис. 3. Области поддерживания магнитного поля в плоскости параметров (к\,к3)

при различных Rem,

Диапазон значений для Rem выбирался го следующих соображений, Поскольку h = 2.1 • 106 м и характерная скорость конвекции v0 по оценкам составляет 10-4 м/с, то Rem = 210/vm, Молекулярное значение магнитной вязкости составляет = 1м2/с [13], а турбулентное значение = 20 м2/с [14], Поэтому значения Rem варьировались так, чтобы перекрыть диапазон значений от турбулентного до молекулярного.

Результаты расчетов для некоторых значений Rem представлены на рис. 3. В плоскости (к\,к3) крестиками обозначены пары коэффициентов, в которых в модели поддерживается магнитное поле. Видно, что с ростом Rem области поддержания поля увеличиваются.

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

что режим усиления магнитного ноля реализуется в широком диапазоне значений магнитного числа Рейнольдса от турбулентного до молекулярного значения.

Литература

1. Копо М., Roberts Р.Н. Recent geodvnamo simulations and observations of the field /7

Reviews of Geophysics. 2002. 40, №10. P.B1-B41.

2. Зельдович Я.Б., Рузмайкин А.А., Соколов Д.Д. Магнитные ноля в астрофизике / М.-Ижевск: НИЦ «РХД», 2006.

3. Cupal I., Hejda P., Reshetnvak M. Dynamo model with termal convection and with the free-

rotating inner core /7 Planetary Physics Sciences. 2002. 50. P.1117-1122.

4. Ustvugov S. Three-Dimensional Numerical MHD Simulation of Solar Convection /7 Hyperbolic

Problems: Theory, Numerics, Applications. 2008. 4. P.1061-1068.

5. Кузнецов В.В. Анизотропия свойств внутренних) ядра Земли // Успехи физ. наук.

1997. 169, №9. С.1001-1012.

6. Водинчар Г.М., Шевцов Б.М. Маломодовая модель конвекции во вращающемся шаровом

слое вязкой жидкости /7 Вычислительные технологии. 2009. 14, №4. С.З 15.

7. Водинчар Г.М., Крутьева Л.К. Маломодовая модель конвекции во вращающемся шаровом слое вязкой жидкости /7 Вычислительные технологии. 2011. 16, № 2. С.35 44.

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

9. International Geomagnetic Reference Field. http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html.

10. Гриненен X. Теория вращающихся жидкостей / Л.: Гидрметеоиздат, 1975.

11. Резников Е.Л., Розенкнон Л.М. О собственных колебаниях вращающейся вязкой жид-

кости во внешнем ядре Земли // Вопросы геодинамики и сейсмологии (Вычислительная сейсмология. Выи. 30) / М.: Геое, 1998. С.121-132.

12. Джекобе Дж. Земное ядро / М.: Мир, 1979.

13. Merril R.T., McElhinnv M.W., McFadden P.L. The Magnetic Field of the Earth / N.-Y.: Acad. Press, 1996.

14. Frick P., Reshethvak М., Sokoloff D. Combined grid-shell approach fo convection in a rotating

spherical layer /7 EuroPhvs Letters. 2002. 59, №2. P.212-217.

6-JET KINEMATIC MODEL OF GEODYNAMO G.M. Vodinchar, L.K. Feshchenko

Institute of Cosmophysical Research and Radio Wave Propagation,

Far Eastern Branch of the Russian Academy of Sciences,

Mirnaya Str., 7, Paratunka, Elizovskiy district, Kamchatka region, 684034, Russia, e-mail: ikir@ikir.ru

Abstract. Geodvnamo kinematic model is built. The field velocity structure which is compatible with data over the splitting-funetions of oscillations in the Earth liquid core. In the model three velocity components are used which are approximate oscillations modes of the liquid core. Two of them define 6-jet convection in the Earth’s core, and the third corresponds to anticipatory rotation of the solid core.

Key words: geodvnamo, kinematic dynamo, Earth’s core, Poincare operator.