Использование собственных мод колебаний вязкой вращающейся жидкости в задаче крупномасштабного динамо Текст научной статьи по специальности «Физика»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 517.958+537.84

ИСПОЛЬЗОВАНИЕ СОБСТВЕННЫХ МОД КОЛЕБАНИЙ ВЯЗКОЙ ВРАЩАЮЩЕЙСЯ ЖИДКОСТИ В ЗАДАЧЕ КРУПНОМАСШТАБНОГО

ДИНАМО *

Г.М. Водинчар1, 2

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

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

E-mail: gvodinchar@ikir.ru

Рассмотрена модель крупномасштабной МГД-конвекциии во вращающейся сферической оболочке. Скорость представлена модой собственных колебаний вязкой проводящей жидкости в этой оболочке. Найдены стационарные режимы конвекции. Установлена возможность работы динамо в этой модели.

Ключевые слова: динамо, МГД-конвекция, сферическая оболочка, собственные колебания

(с) Водинчар Г.М., 2013

MATHEMATICAL SIMULATION

MSC 37N10+85A30

USING MODES OF FREE OSCILLATION OF A ROTATING VISCOUS FLUID

IN THE LARGE-SCALE DINAMO

G.M. Vodinchar1, 2

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@ikir.ru

Considered a model of large-scale magneto-convection in a roteting spherical shell. Velocity expressed of the one of the modes of free oscillation of a viscous conduction fluid in the shell.

The stationary modes of convection were found. An ability of the dynamo was approved in this model.

Key words: dynamo, magneto-convection, spherical shell, free oscillation

(c) Vodinchar G.M., 2013

*Работа выполнена при поддержке Минобрнауки России по Программе стратегического развития КамГУ им. Витуса Беринга на 2012-2016 гг.

Введение

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

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

При разработке крупномасштабных моделей возникает необходимость задавать из каких либо соображений геометрическую структуру конвекции, поэтому интересно рассмотрение простейших моделей с «базисных» в некотором смысле геометрическими структурами течений. В качестве важнейшего примера надо назвать классическую систему Лоренца [4], которая сыграла огромную роль в понимании многих свойств конвекции и детерминированного хаоса. Отметим, что система Лоренца получена именно путем выделения «базисных» течений. В задаче конвекции в плоском слое поле скорости разложено по собственным модам свободных затухающих колебаний вязкой жидкости в слое, а температура разложена по модам, пространственно согласованным с модами скорости. Далее проведено предельное усечение, сохраняющее нелинейность - одна мода скорости, согласованная с ней температурная мода и еще одна температурная мода, однородная по граничной поверхности. Затем стандартная галеркинская процедура дает динамическую систему для амплитуд -систему Лоренца.

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

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

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

Уравнения модели

Прежде всего уточним постановку задачи динамо, которую далее будем обсуждать.

Рассматривается конвекция проводящей замагниченной жидкости в сферической оболочке Б с твердыми границами. Жидкость находится в центральном поле тяготения вида g = — (^2/^2) г, где радиус-вектор отложен от центра оболочки О, а g2 - ускорение свободного падения на внешней границе г = Г2. Снаружи оболочки среда непроводящая, внутри - проводящее твердое тело. Магнитные проницаемости оболочки и внутреннего тела постоянны и совпадают. Оболочка и внутреннее тело вращаются вокруг оси, проходящей через точку О с постоянной угловой скоростью а. На границах оболочки температуры постоянны, причем на внутренней выше, чем на внешней.

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

д v T 2

Pr-1 — + (vv)v = Av — vp' + Ra—г--ez x v + rotB x B,

д t v 7 r2 E z

^ + (vv)(T + T,) = AT,

^ = rot (v x B) + Ro-1 A B, ()

О t

Vv = 0, vB = 0,

где Pr - число Прандтля, Ra - число Релея, E - число Экмана, Ro - число Робертса, поле T задает отклонение температуры от стационарного гиперболического профиля Ts, а р' - редуцированное давление, собравшее в себя все потенциальные слагаемые из правой части уравнения Навье-Стокса.

Система (1) замыкается нулевыми условиями для скорости и температуры на границах оболочки и вакуумными граничными условиями для магнитного поля [8].

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

-1 / 2

Pr 1 д v + Av-VP - Eez x v, (2)

vv = 0,

c нулевыми условиями на границах оболочки.

Будем искать ее приближенные решения в виде комбинаций тороидальных и поло-идальных собственных мод колебаний невращающейся оболочки уТпт и УРпт. Явные выражения для них и уравнения на собственные значения пт и дкПтописаны в [9]. Можно показать, что при фиксированном т > 0 линейные оболочки Нт всех тороидальных и полоидальных полей уТп(±)т и уРп(±)т образуют инвариантные подпространства оператора задачи (2). Поэтому задачу можно решать отдельно в каждом таком подпространстве.

Выберем теперь в одном из этих подпространсв некоторую конечную систему из N полей из числа уТп(±)т и ^рП(±)т, которые будем обозначать далее как у,.

N

Ищем решения задачи (2) в виде ^Д'уг, используя метод Галеркина. Известно,

'=1

что применение галеркинской процедуры при нулевых граничных условиях исключает из этой задачи давление [10].

Несложно показать, что тогда столбцы коэффициентов в, будут собственными векторами матрицы с (і,;)-элементами, равными

| ^ (ег X V;) у,^ ,

где интегрирование ведется по объему оболочки, а 8,; - символ Кронекера.

Отберем пару сопряженных собственных значений этой матрицы и соответствующие им сопряженные векторы (наборы коэффициентов Д'). Каждый такой набор определит аппроксимацию одной собственных мод задачи (2). Теперь сложим эту пару аппроксимаций. В результате получим действительное стационарное поле и (г) = <^1 у і (г) + ... + ^vN(г), которое и будем использовать для одномодового приближения скорости в задач динамо, предварительно отнормировав его по объему оболочки.

Далее необходимо составить моду температуры То (г), структурно согласованную со скоростью. Это предполагает, что в тех частях оболочки где положительна (отрицательна) радиальная составляющая поля и (г), должна быть положительна (отрицательна) и То (г). Поскольку тороидальные компоненты не имеют радиальной составляющей, то все определяется только полидальными модами среди у,.

^Р (г)

Пусть, например, у1 = уРпт. Тогда ее радиальная проекция - п(п +1) ——^”(0, Ф)

и с ней согласуется необходимым образом собственная мода скалярного оператора Лапласа в оболочке с нулевыми граничными условиям Ккп(г)УПт(0, ф). Здесь радиальные функции ^Рп(г) и ^кп(г) определяются из соответствующих спектральных задач, обе нулевые на границах и имеют к нулей в интервале (п;г2). При этом соответствующие нули функций практически совпадают. Все детали расчета таких функций описаны в работе автора [9].

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

Также включим в модель еще однородную на сфере моду Т1 (г) с одним нулем на интервале (п;г2). Это мода ^ю(г)ГоО(0,ф). Она вводится в модель по аналогии с построением системы Лоренца.

Тогда подстановка разложений у (г, г) = в (г)у (г) и Т (г, г) = ао(г)То (г) + а1(г)Т1 (г) в первые два уравнения системы (1) и применение галеркинской процедуры даст

систему

Pr 1 = —aß + RaCa + í u(rotB x B)dV,

—,— = Fojß«l + H ß — dt

(3)

Постоянные коэффициенты этой системы

a = Ml^l2 + ... Mn kN > o, dV,

Fo = — / To (uv) TldV, Fl = — / Tl (uv) TodV,

Jd Jd

H = — / To (uv) TsdV, Ai = —Ro—1 / T Л TdV > o

dV,

Также при выводе системы учитывались ортогональность сферических гармоник и равенство нулю интеграла / и (иу) udV, легко устанавливаемое с помощью теоре-

мы Остроградского-Гаусса, для соленоидального нулевого на границе оболочки поля и (г). Кроме этого несложно установить, что всегда ^1 = — Fo.

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

Отметим также , что при нулевом магнитном поле и подходящем перемасштаби-ровании амплитуд и времени система (3) превращается в систему Лоренца.

Представим теперь и магнитное поле В(г,t) в виде комбинации K стационарных мод:

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

Подставим разложения (4) в уравнение индукции (третье уравнение системы (1) и применим галеркинскую процедуру. Получим набор уравнений:

K

Bfot) = £ Ys(t

(4)

s=1

s = 1,..., K.

Подставив также разложения (4) в интеграл в первом из уравнений (3) и собрав уравнения (3) и (5) в одну систему получим:

dß к

Pr 1 ~г~ = — + RaCa + £ LijУгУj,

dt i, j=1

— = Foßa1 + H e — dt

d«l p o a (6)

— = Flßao — Al al, w

к d в к к

I ßsi^7 = I Wsi-вТї- — Ro—1 £ MsíYí, i=l dt i=l i=l

s = 1,..., к,

где Ьг-7- = J и (го1Вг- х Вj) dV .

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

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

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

Уравнения крупномасштабного динамо

Рассмотрим теперь случай, когда используются всего две магнитные моды. Первая В1 - полоидальная, например вертикальный диполь, а вторая В2 - тороидальная, согласованная структурно с модой скорости. При этом будем считать выполненными вышеприведенные требования относительно коэффициентов Wsi и Lij.

Поскольку подпространства тороидальных и полоидальных полей инвариантны относительно оператора Лапласа и ортогональны на сфере, то матрица Q единичная, а М диагональная с положительными Иц. Также несложно показать, что диагональные элементы матриц W и Ь нулевые.

З8

Система (6) тогда примет вид

dß

Pr 1 ~dt = —uß + RaCa + hiWі, d ao dt d a1 dt

dßl

dt

dß2 dt

= Foßal + H ß —

= F1ßao — A1 a1,

= Wl2ßY2 — Ro !Міі7і, = WZlßYl — Ro ^M22Y2.

(7)

Поскольку одновременнная смена знаков коэфициентов ^ равносильна смене знака амплитуды «1, можно считать, что ^0 < 0 и > 0. Смена знака ^2 эквивалентна смене знака у одной из амплитуд $ и у обоих коэффициентов И12 и И2ь Поэтому, не теряя общности, можно считать, что Ь12 < 0.

Выполним в системе (7) замену времени и амплитуд по следующим формулам:

t = — T ß = Ao

= Ao T, ß = </№1

öl, Yl =

м, ao =

Ao

x/PrlLlzWzll

Ao u

RaCv^FIF Bl,

öo,

(8)

Рг|^сЬі2І^і В2 В новых переменных она примет вид:

— = о (öo — м) — BlB2, -т

dö

- T dö1

- т -B1 d t -B2

d T

— —м$1 + гм — öo, = м öo — b öl,

= p^B2 — cB1,

= є ^Bl — /B2,

(9)

где а = 7^Рг > 0, г = Яа > 0, Ь = А1/А0 > 0, р = ^^^|^211 , е = sgnW2l, с = > 0,

Ао -0 д |ЗД Яо-0

М22

7=яоА0 > 0

Видно, что если положить в этой системе Ві = В2 = 0, то она превратиться в систему Лоренца. Если же оставить в ней первое, четвертое и пятое уравнения и считать в первом 0 постоянной, то возникает система, похожая на систему двухдискового динамо Рикитаки с трением, когда трение в обоих дисках совпадает [11, 12].

Далее в работе анализируются стационарные режимы, возникающие в системе (9) при варьировании параметров.

Стационарные режимы в модели крупномасштабного динамо

Рассмотрим состояния равновесия, возникающие в системе (9).

Прежде всего отметим, что у системы (9) есть нулевая точка покоя .Оо, устойчивая при г < 1, и теряющая устойчивость при г = 1. При г > 1 возникают еще точки

Оц,ь2 = (±^Ь(г — 1);±^Ь(г — 1);г — 1;0;0^. Это все «лоренцевские» точки.

Наличие точек покоя с ненулевыми значениями магнитного поля существенно зависит от знаков параметров р и е.

Случай ре > 0.

В этом случае к точкам могут добавиться еще четыре с ненулевыми зна-

чениями магнитного поля («магнитные» точки):

fc. rb\Jfc\p\. rfc. la f IrbI p I - d, e joc\rb\p\-d\

I pi; d ; d ;V d ; V I p I d

fc ; rb\Jfc i pI ; rfc jaf I rb I p I - d ; e la c I rb I p I - d I

p d d

fc rbVfc i p yrfc i p I

p d

d

d

a f rb p - d ac rb p - d

(10)

;e

p d

fc . - rbVfc I p I ; laf I rb I p I - ~d~i ;-e la c I rb I p I - d I

d

d

p d

где d = /о + Ь | р |.

Очевидно, что все точки покоя разбиваются на три группы Оо, О^—¿2 и Ом1—м4, в пределах каждой из которых точки переходят друг в друга при симметриях (отражениях в фазовом пространстве) системы (9). Поэтому геометрии фазовых траекторий в окрестностях точек одной группы одинаковы.

При р > 0 условие возникновения точек Ом1—М2 имеет вид принимает вид

r > 1 + т-, bp

(11)

fc

т.е. в этом случае с ростом числа Релея сначала, при 1 < г < 1 + -—, возникает кон-

Ьр

векция недостаточная для поддержания магнитного поля, а затем, при выполнении (11), появляется режим стационарной МГД-конвекции.

Разумеется, говорить о физических стационарных режимах можно только в том случае, если соответствующие точки покоя устойчивы. Аналитическое исследование на устойчивость ненулевых точек крайне затруднительно. Можно лишь сказать, что в случае (11) «лоренцевские» точки неустойчивы в плоскости магнитных переменных.

Было проведено численное исследование исследование устойчивости «магнитных» точек по первому приближению при случайном варьировании параметров системы (9) в следующих диапазонах:

10—2 < а < 105, 10—1 < Ь < 102,

10—4 < р < 106, 102 < о < 106,

102 < / < 106, 1 < г < 103,

При имитации значений параметров их распределение было равномерным в логарифмическом масштабе. Всего был проведен 1 млн. имитаций. Каждый раз проверялась устойчивость «лоренцевских» точек, а если выполнялось (11) , то и «магнитных».

Во время имитации возникали случаи как устойчивых, так и неустойчивых «магнитных» режимов. При этом «лоренцевские» точки были как устойчивы, так и неустойчивы в пространстве амплитуд скорости и температур независимо от устойчивости «магнитных» точек.

Поэтому, в рассматриваемом случае, в модели могут реализовыватьс различные режимы как гидродинамической, так и МГД конвекции.

/о

В случае р < 0 точки Ом1-м2 являются при г < 1-,—г, т.е. «магнитные» точки

Ь| р1

/о

могут возникнуть при как угодно малых значениях числа Релея, затем, при 1-,—г <

Ь| р|

г < 1, ненулевых точек покоя нет, а при г > 1 появляются только «лоренцевские». Физически это выглядит неправдоподобно.

Было проведено численное исследование на устойчивость этих точек при тех же варьированиях параметров а, Ь, о, /, что и выше. Параметр р варьировался в диапазоне [—106; —10—4] равномерно в логарифмическом масштабе, а параметр г был равномерно распределен на отрезке [0;шт{0,1 — /о/(Ь|р|)}]. Всего было проведено 10 млн. имитаций и во всех «магнитные» точки оказывались неустойчивыми.

Таким образом, при р < 0 физических стационарных режимов динамо в рассматриваемой модели не существует.

Случай ре < 0.

В такой ситуации система (9) имеет только «лоренцевские» точки покоя, поэтому он далее не рассматривается.

Заключение

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

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

Фактически, рассмотренная модель объединяет в себе идеи модели Лоренца и динамо Рикитаки с трением.

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

1. Зельдович Я.Б., Рузмайкин А.А., Соколов Д.Д. Магнитные поля в астрофизике. М. : Ижевск: Регулярная и хаотическая динамика, 2006. 384 с.

2. Glatzmaier G., Roberts P. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal // Nature. 1995. Vol. 377. P. 203-209.

3. Kuang W., Bloxman J. An Earth-like numerical dynamo model // Nature. 1997. Vol. 389. P. 371-374.

4. Lorenz E. Deterministic nonperiodic flow // Journal Atmos. Sci. 1963. Vol. 20. P. 130-141.

5. Гринспен Х. Теория вращающихся жидкостей. Л.: Гидрометеоиздат, 1975. 304 с.

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

7. Рихтмайер Р. Принципы современной математической физики. Т. 1. М.: Мир, 1982. 486 с.

8. Merril R., McElhinny M., McFadden P. The Magnetic Field of the Earth. N.Y.: Acad. Press, 1996. 532 p.

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

10. Ладыженская О.Н. Математические воросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 232 с.

11. Ershov S.V., Malinetskii G.G., Rusmaikin A.A. A generalized two-disk dynamo model // Geophys. Astrophys. Fluid Dynam. 1989. V. 47. P. 251-277.

12. Потапов В.И. Визуализация фазовых траекторий динамической системы Рикитаки // Нелинейная динамика. 2010. Т. 6. № 2. C. 255-265.

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