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

DOI: 10.18454/2079-6641-2015-11-2-45-54

математическое моделирование

УДК 517.9

инверсии в модели геодинамо, управляемой 6-ЯЧЕЙКОВОЙ

конвекцией

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

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

2 Камчатский государственный университет имени Витуса Беринга, 683032, г. Петропавловск-Камчатский, ул. Пограничная, 4

E-mail: gvodinehar@gmail.com

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

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

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

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

mathematical modeling

MSC 97А80

reversals in the 6-cells convection driven

G.M. Vodinehar1, 2, L.K. Feshenko2

1 Institute of Cosmophysical Researches and Radio Wave Propagation Far-Eastern Branch, Russian Academy of Sciences, 684034, Kamchatskiy Kray, Paratunka, Mirnaya st., 7, Russia

2 Vitus Bering Kamchatka State University, 683031, Petropavlovsk-Kamchatsky, Pogranichnaya st., 4, Russia

E-mail: gvodinchar@gmail.com

We describe the large-scale model geodynamo, which based on indirect data of inhomogeneities in the density of the Earth’s core. Convection structure is associated with spherical harmonic Y42, which defines the basic poloidal component of velocity. Coriolis drift of this mode determines the toroidal component of velocity. Thus, 6 convective cells are formed. The model takes into account the feedback effect of the magnetic field on convection.It was ascertained that the model contains stable regimes of field generation. The velocity of convection and the dipole component of the magnetic field are close to the observed ones.

Key wards: geodynamo, geomagnetic field reversals, convection

(c) Vodinchar G.M., Feshenko L.K., 2015

Введение

Процесс формирования магнитных полей планет и звєзд успешно объясняется теорией гидродмагнитного динамо. Существующие модели конвекции в жидких ядрах планет и конвективных зонах звезд допускают структуры течений, способные генерировать магнитные поля, подобные реальным [1, 2, 3].

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

Поэтому существующие численные модели либо воспроизводят МГД-тєчєния с хорошим разрешением по пространству на относительно небольших (по геологическим масштабам) промежутках времени (~ 104 лет), либо воспроизводят длительную эволюцию наиболее крупномасштабных пространственных структур. В моделях первого типа пространственная структура течений и поля просчитывается, а для моделей второго типа эта структура должны изначально задаваться. При этом возникает вопрос о том, какова реальная крупномасштабная структура конвекции.

Косвенные данные об этой структуре для ядра Земли можно получить из данных о пространственной неоднородности плотности жидкого ядра. В работе [4] проанализированы результаты исследований функций расщепления собственных колебаний Земли после сильнейшего Чилийского землетрясения. Распределение этих функций характеризует неоднородность ядра Земли на различных глубинах. Проанализировав эти функции, автор [4] пришел к выводу, что в жидком ядре существует 12-зонная шахматного типа неоднородность, которая естественным образом описывается в первом приближении сферической гармоникой Y42. Им же была выдвинута гипотеза о том, что в ядре существует конвекция с 6 областями подъема вещества и 6 областями опускания, т.е. о 6-ячєйковой конвекции. В такой конвективной структуре в каждой из полярных областей расположены по 2 ячейки и еще 2 лежат в экваториальной области.

Отметим также, что 6-ячейковая структура возникает в результате прямого численного моделирования конвекции в ядре Земли при некоторых значениях параметров ядра [5]. В этих расчетах в каждом полушарии расположены по 3 ячейки. Такая структура описывается сферической гармоникой Y43.

Известным свойством реальных космических динамо-систем является наличие инверсий - смены знака дипольной составляющей поля. При этом Большой интерес представляют инверсии не связанные с перестройкой конвекции.

В настоящей работе исследуется вопрос о генерации поля 6-ячєйковой структурой, описываемой гармоникой Y42 с учетом турбулентного а-эффекта, и об инверсиях в такой системе.

Основные уравнения

Рассматриваем сферическую оболочку вязкой несжимаемой жидкости (жидкое ядре), вращающуюся вокруг z-оси с постоянной угловой скоростью ■. В сферической системе координат (r, в, ) внутрення граница r = r и внешняя граница r = ro. Температура на внутренней и внешней границах постоянна и равна, соответственно, Ti и To.

Уравненя геодинамо, учитывающие а-эффект в приближении Буссинеска имеют вид:

Е д + v ■ V ) v = EV2v — Vp — 2ez x v + RaPmTr + rotB x B,

Pm у dt d- + (v ■ V)(t + Ts ) = P!m v2t ,

= rot (v x B) + Rarot (aB) + V2B, ^ ^

о t

V ■ v = 0,

V ■ B = 0.

Основными безразмерными параметрами (параметрами подобия) являются число Экмана Е = v/О-r^, число Релея Ra = вg08Tro/vQ., число Прандтля Pr = v/к, магнитное число Прандтля Pm = v/ц и амплитуда а-эффекта Ra. В этих выражения v обозначает кинематическую вязкость, в - коэффициент объемного расширения, go -ускорение свободного падения на внешней границе, 8T = T — To - разность значений температуры на внутренней и внешней границах, к - коэффициент температуропроводности и п - коэффициент омической диссипации. Такой набор параметров возникает, если ro, r^/ц, 8T и л/ЩрШц, где р - плотность и р - магнитная проницаемость, являются единицей длины, единицей времени, единицей температуры и единицей магнитного поля, соответственно. Радиус внутреннего ядра, выраженный в такой единице длины, составляет r = 0.35.

В уравнениях (1) переменная T обозначает отклонение температуры от равновес-

„ 1/r — 1

ного гиперболического r-профиля Ts = —----- + Ті — 1.

1/ri — 1

Мы считаем, что турбулентность в ядре изотропная и используем скалярную параметризацию а-эффекта в виде a = a (r) cos 0, где max |a(r)| = 1.

Для скорости v принимаются граничные условия прилипания. Также мы считаем, что магнитная проницаемость внутреннего и внешнего ядра одинаковы и среда вне ядра (r > ro) непроводящая. Таким образом, для магнитного поля принимаются вакуумные граничные условия на внешней границе и условия ограниченности в центре Земли. Наконец, отклонения температуры T на границах внешнего ядра считаются нулевыми.

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

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

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

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

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

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

Е 2

—— иv + EV2v — 2ez x v — Vp = 0,

Pm z

V ■ v = 0, (2)

v(r = r) = v(r = ro) = 0.

Для удобства мы использовали те же безразмерные параметры, что и в уравнениях (1).

Нєвязкий аналог этой задачи известен как задача Пуанкаре и является косоэрмитовым [7]. Далее мы будем называть решения задачи (2) модами Пуанкаре. В вязком случае оператор задачи не обладает эрмитовой или косоэрмитовой симметрией, что создает значительные трудности. Точные решения в настоящее время не известны. Есть только некоторые оценки спектра, доказана полнота системы собственных и присоединенных полей этой задачи [8]

Мы будем использовать простую аппроксимацию решений задачи (2), основанную на хорошо известных модах свободных колебаний невращающейся оболочки. Формально они определяются как решения спектральной задачи

^—uv + V2v — V p = 0, V ■ v = 0,

Pm^ y ’ ’ (3)

v(r = r) = v(r = Го) = 0.

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

Тороидальные и полоидальные моды, соответственно, имеют следующий вид:

vk,n,m Rk,n(r)

vp,n,m = n(n + 1)

d

d

ee^ — e^) ОТ, Ф),

sine едф vdej —P (r)

Ynm(e, ф)er+

+ (dr+г) -P,n(r) (eede+Sine ЄфІф) c (e, ф),

k = 0,1,..., n = 1,2,..., m = —n,..., n,

1

(4)

где -Tn(r) и -pn(r) некоторые комбинации степенных функций и сферических функций Бесселя, а Ynm(e, ф) - сферические гармоники.

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

HT = { vT1,1,0, vT3,3,0: 1 vp4,4,0, . . . } ,

HP = Г P T P 1 Vk1,1,0, Vk2,2,0, vk3,3,0, vW...} ,

HT = m= Г vT , vP 1 Vko,m,±m, Vk1 ,m+1 ,±m' vT > Vk2,m+2,±m' vP > Tk3,m+3,±m

HP = m vT

Г vP , vT , 1 vk0,m,±m, Vk1 ,m+1 ,±m, vP , vk2,m+2,±m, vk3,m+3,±m'

ki = 0,1,2,..., m = 1,2,3,...

(5)

Каждое из этих подпространств инвариантно относительно оператора задачи (2). Это легко проверить, если вычислить матрицу оператора в базисе (4). Поэтому каждая из мод Пуанкаре полностью лежит в одном из подпространств (5).

Для представления 6-ячєйковой конвекции важно подпространство Н2р, содержащее полоидальную моду vp4±2. Из уравнений (4) хорошо видно, что радиальная компонента этой моды в (0, ф)-направлении определяется сферической гармоникой

Y ±2

Y 4 .

В нашей модели мы используем только одну моду скорости Vo, являющуюся простейшей аппроксимацией мод Пуанкаре из подпространства Н2р. Мы определим vo в виде

V0 = вvf,3,-2 + №1,3,2 + &<4,-2 + e4Vp>,4,2 + №1,5,

■ + №f

0,5,2’

(6)

где коэффициенты ві находятся из задачи (2) методом Галеркина.

Ясно, что радиальная компонента vo имеет в (0, ф)-направлении структуру гармоники Y4±2.

Теперь рассмотрим собственные моды тєпловой диффузии во внешнем ядре. Они задаются в виде Tk,n,m = Xk,n(r)Y,m, где Xk,n(r) являются линейными комбинациями сферических функций Бесселя первого и второго родов.

В нашей модели главная температурная мода To должна быть согласована с радиальной компонентой скорости, Поэтом определим ее в выражением

T0 = X0,4(r) [ftY4 2 + №4] = в3Т0,4,-2 +

(7)

где в3 и в4 взяты из формулы (6).

Следующая температурная мода Т1 выбирается в виде Т1 = Т1,0,0 = X1,0Y00 подобном моде из системы Лоренца.

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

v(r, t) = u(t) V0 (r), T (r, t) = 00(t)T)( r) + 01(t )T1 (r). (8)

Магнитная часть модели

Для представления магнитного поля мы используем некоторые моды омичєской диссипации, т.е. решения спектральной задачи

ПВ + V2B = 0, V ■ B = 0,

B(r = 0) = ~, BT (r > 1)= 0, rotBP(r > 1) = 0,

(9)

где BT и BP тороидальная и полоидальная составляющие магнитного поля, соответственно.

Тороидальные и полидальные решения этой задачи

ть=aL roM

¥“(Q, V )r

?

= apnrotrot jn

¥n(Q, V)r

(10)

где jn(-) сферические функции Бесселя первого рода. Собственные значения ПП and ПІП определяются из граничных условий, а aTn и aPn являются нормирующими ко-

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

Для отбора магнитных мод мы использовали следующую схему, предложенную в работе [9].

Прежде всего, мы упорядочим моды по возрастанию собственных значений, т.е. по скорости диссипации. Получим следующую последовательность, где прямоугольниками выделены дублеты - пары мод с одинаковыми собственными значениями: 1..1 т~1-І т>—2..2 гр—2..2 и—3..3 и—1..1

p

01

01

p

2..2 02

т

2..2 02 ,

p

03

p

11

т

3..3

03

p

-4..4

04

т

1..1

11

p

2..2 12

гр—4..4 т>—5..5

т04 , p05

Затем, выберем несколько мод с наименьшими собственными значениями, обозначив их B|(r).

Теперь мы запишем магнитное поле в виде B(r,t) = ^gk(t)B|(r) и подставим

k

это выражение в уравнение индукции (трете из уравнений (1)). Скорость задаем в виде (8), где временно считаем, что u(t) = u0 = const > 0. Моду скорости v0 считаем нормированной так, чтобы ее среднеквадратическое значение равнялось единице. Тогда U0 можно интерпретировать как магнитное число Рейнольдса Rem.

На следующем шаге применим к уравнению индукции метод Галеркина и получим систему

dgk

dt

Rem ^ Wkigi + Ra ^ Akigi i i

Wki = J rot (V0 X Bi) BkdV,

Aki = J rot (a(r) cos QB;) BkdV,

(11)

где Пі - это собственное значение для Bi.

Фактически, мы получим маломодовую аппроксимацию для задачи кинематического динамо с шестью конвективными ячейками. Пусть Я обозначает собственное значение матрицы системы (11). Кинематическое динамо генерирует поле, если

maxReA; > 0. Мода с таким собственным значением является лидирующей, растущей быстрее других. Если собственное значение лидирующей моды комплексное, то магнитное поле будет осциллировать.

Далее, мы будем постепенно наращивать число магнитных мод, пока не получим осциллирующее динамо.

Такова схема отбора магнитных мод.

Разумеется собственные значения зависят от параметров Rem и Ra, а также от формы a(r) радиального профиля а-эффекта. Мы провели серию расчетов, используя два выражения для профиля a(r) = 1 и a(r) = r. Оказалось, что результаты отличаются незначительной этих расчетах параметры Rem и Ra равномерно варьировались в логарифмической шкале в диапазоне [10-1;103]. После получения осциллирующего динамо мы отбрасывали некоторые моды, если это не меняло качественно результата или не меняло существенно порог генерации.

В результате мы остановились на магнитных модах P01, P±32, P^, T±42, P±52.

Полученная схема генерации магнитных мод крупномасштабной конвекцией и а-эффектом показана на рис. 1, а области осциллирующего и неосциллирующего динамо на плоскости параметров (Rem, Ra) показаны на рис. 2.

Рис. 1. Схема динамо. Сплошные стрелки - крупномасштабный генератор; пунктирные стрелки - а-генератор.

Рис. 2. Области генерации поля. Красные точки - неосциллирующее динамо; зеленые точки - осциллирующее динамо. a(r) = 1 (вверху) и a(r) = r (внизу).

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

Нелинейная модель МГД-конвєкции

В вышеописанном кинематическом динамо возможен только неограниченный рост поля. Для получения устойчивой генерации необходимо водить механизм подавления. В связи с этим обратим внимание на области на рис. 2 выделенные эллипсами. Зафиксируем Ra ~ 15. Если точка (Rem, Ra) изменяется вдоль Большой оси эллипса, то значение Rem также варьируется. При этом в крайних положениях динамо не работает. Так можно получить генерацию поля ограниченной величины.

Напомним, что в нашей конструкции Rem является фиксированным уровнем амплитуды u(t) моды скорости vo. Поэтому мы можем составить теперь комбинированную модель МГД-конвєкции с переменной амплитудой скорости. Для ее получения подставим выражения для скорости, температуры и магнитного поля в первые три уравнения (1) и применим метод Галеркина.

Получим нелинейную динамическую систему, описывающую 6-ячейковую МГД-конвєкция в ядре Земли:

Е du(t)

Pm dt

d Qo(t) dt

d 61 (t) dt

dgk (t) dt

k = 1,..., 8,

= -Едu(t) + RaPmSflo(t)+ £ Lijgi(t)gj(t),

i,j=1 Pm

= -Fu(t)01 (t) + Hu(t) - — ZoOo(t),

Pr

Pm

= Fu(t )do(t) - — C1fl1(t ), 88

= u(t) £ Wkigi(t) + Ra £Akigi(t) - nigi(t),

i=1

i=1

(12)

где Zi собственные значеия температурных мод Ті, а д > o, S, Lij, F, H являются коэффициентами Галеркина.

Результаты вычисленных экспериментов с моделью

Опишем теперь некоторые результаты вычислительных экспериментов с моделью (12). Мы брали турбулентные значения для диссипативных коэффициентов V = 1o2 м2/с, к = Ш-2 м2/с, п = 2o м2/с [10]. Внешний радиус ядра ro = 348o км а угловая скорость П = 7.29 х 1o-5 рад/с. Соответствующие значения параметров модели: E = 1o-7, Pr = 1o4, Pm = 5. Также мы использовали Ra ~ 15, следуя оценкам из работы [11] и различные значения числа Релея из диапазона Ra = 1o2 ^ 1o4.

Система уравнений (12) численно решалась с применением стандартных процедур пакета SeiLab.

Были найдены устойчивые режимы динамо с инверсиями в магнитном поле для Ra ~ 1o3. На рис. 3 показаны графики амплитуд скорости и дипольной части магнитного поля для Ra = 5ooo.

t, in units of L In

0.5

t, in units of L In

1.5

0

2

Рис. 3. Амплитуды мод скорости и вертикального диполя для Ra = 5000

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

Разумеется, к конкретным численным значениям, полученным в такой простой модели следует относиться с осторожностью. Тем не менее отметим, что характерное значение скорости на рис. 3 составляет 5 х 10-4 м/с. Это хорошо согласуется с известными оценками реальной скорости конвекции ~ 10-4 м/м [1]. Характерное значение поля в модели 5 х 10-4 Тл также хорошо согласуется с экстраполяцией реальной величины геомагнитного поля на границу ядро-мантия [12]. Временной интервал между инверсиями в модели 3200 лет, что коррелирует со средней величиной интервала полярности для палеомагнитных шкал инверсий [12].

Выводы

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

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

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

Характерные значения скорости и поля в модели совпадают по порядку величины с известными оценками реальных скорости и поля.

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

Библиографический список

1. Jones С.А. Convection-driven geodynamo models // Phil. Trans. R. Soc. bond. A. 2000. VoI. 358. pp. 873-897.

2. Kono M., Roberts P.H. Recent geodynamo simulations and observations of the field // Reviews of Geophysics. 2002. Vol. 40. pp. Б1-Б41.

3. Соколов Д.Д., Степанов P.A., Фрик П.Г. // Динамо: на пути от астрофизических моделей к лабораторному эксперименту // Успехи. физ. наук. 2014. Т. 184. № 3. С. 313-335.

4. Kuznetsov V.V. The anisotropy of properties of the Earth’s inner core // Physics-Uspekhi. 1997. Vol. 40. pp. 951-961.

5. Гореликов А.Б., Ряховский А.Б., Фокин A.C. Численное исследование некоторых нестационарных режимов естественной конвекции во вращающемся сферическом слое // Вычислительная механика сплошной среды. 2012. Т. 5. № 2. С. 184-192.

6. Lorenz E.N. Deterministic nonperiodic flow // J. Atmos. Sci. 1963. Vol. 20. pp. 130-141.

7. Greenspan H.P. The Theory of rotating Fluids. New York: Cambridge University Press, 1968.

8. Резников Е.Л., Розєнкноп Л.М. О собственных колебаниях вращающейся вязкой жидкости во внешнем ядре Земли // Вопросы геодинамики и сейсмологии (Вычислительная сейсмология. Бып. 30). М.: Геос, 1998. С. 121-132.

9. Соколов Д.Д., Нефедов С.Н. Маломодовое приближение в задаче звездного динамо // Вычислительные методы и программирование. 2007. Т. 8. С. 195-204.

10. Решетияк М.Ю. Моделирование процессов динамо в геофизике / Дисс. на соиск. уч. степ. д-ра физ.-мат. наук, Москва: ОИФЗ РАН, 2003.

11. Ануфриев А.П., Решетияка М.Ю., Соколов Д.Д. Оценка динамо-числа в модели турбулентного a-эффекта для жидкого ядра Земли // Геомагнетизм и аэрономия. 1997. Т. 37. С. 141-146.

12. Merril R.T., McElhinny M.W., McFadden P.L. The Magnetic Field of the Earth: Paleomagnetism, the Core, and the Deep Mantle. london: Academic Press, 1996.

Поступила в редакцию / Original article submitted: 17.11.2015