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

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

Том 14, № 4, 2009

Маломодовая модель конвекции во вращающемся шаровом слое вязкой жидкости

Г. М. ВодинчАР, Б.М. Шевцов Учреждение Российской академии наук Институт космофизических исследований и распространения радиоволн ДВО РАН, п. Паратунка Камчатского края e-mail: gvodinchar@yandex.ru, bshev@ikir.kamchatka.ru

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

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

Введение

Конвективные процессы во вращающихся шаровых оболочках — основные формы движения жидкостей и газов в астро- и геофизических объектах — играют важнейшую роль в генерации магнитных полей планет и звезд [1-4].

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

1. Уравнения конвекции

Рассмотрим шаровой слой О вязкой жидкости толщиной к, вращающийся с постоянной угловой скоростью О вокруг оси с направляющим ортом ех. Будем использовать сферическую систему координат с началом в центре О слоя, совместив ось вращения слоя с осью в = 0. Температуру на внутренней и внешней границах слоя считаем постоянной, причем на внутренней границе она выше, чем на внешней. Ускорение свободного падения постоянно по величине и равно g = — дег, где ег — радиальный орт.

Используем в качестве единиц измерения (характерных масштабов): по длине — толщину слоя к, по температуре — разность температур на границах слоя 8Т, по времени —

*Работа выполнена по программе фундаментальных исследований Президиума РАН. 2008. № 16. ч. 3. © ИВТ СО РАН, 2009.

характерное время вязкой диссипации t0 = h2/v, где v — коэффициент кинематической вязкости. За единицу давления примем p0v2/h2, где р0 — средняя плотность в слое.

Уравнения конвекции в приближении Обербека—Буссинеска запишем в принятой в работах по геодинамо инвариантной форме [10, 11]:

d v 1

— + (vV) v = Av — Vp + Ra Pr-1Ter - т (ez x v) ,

dT

— + (vV) T = Pr-1 AT, (1)

Vv = 0,

где v, T и p — возмущения полей скорости, температуры и давления соответственно. Управляющими параметрами модели являются: число Рэлея Ra = STgh3/3/(vk), число Прандтля Pr = v/k, число Кориолиса т = 2h2Q/v = 4nt0/T0. Здесь в — коэффициент объемного расширения, k — коэффициент температуропроводности, T0 = 2п/П — период вращения слоя.

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

T U = Tb T U = T1 - 1, v |r=ri = v|r=r2 = °

где r1 — внутренний радиус, а r2 = r1 + 1 — внешний радиус слоя.

Перейдем в системе (1) к однородным граничным условиям для температуры, введя отклонение в температуры от линейного по радиусу профиля, и получим выражение T (r) = —r + r1 + T1 + в (r). Затем возьмем ротор от обеих частей первого уравнения для исключения поля давления, учитывая, что (vV) v = (1/2) grad v2 — v x rot v. Получим систему

dv _■■

rot — — rot (v x rot v) = rot Av + RaPr rot (вег) — тrot (ez x v) ,

дв + (vV) в — (vV) r = Pr-1 f—2 + Ав^| , dt V r ) (2)

Vv = 0,

e|r=ri = в|г=г2 = ° v|r=ri = v|r=r2 =

Следуя идее метода Галеркина, далее разложим поля температуры и скорости на базисные моды с нестационарными амплитудами [2, 5].

2. Разложение полей скорости и температуры

Будем использовать сферическую систему координат с началом в точке О, совместив ось вращения слоя с осью в = 0.

Рассмотрим пространство Н скалярных полей в слое Д со скалярным произведением <ф*>■ = <интегрирование ведется по объему ^удовлетворяющих

условию Ф|дд = 0. Для разложения температуры используем собственные поля оператора Лапласа в пространстве Н1. Их можно представить в виде

©пш = Яп (г, А) Упш (в,р), п = 0,1,..., т = -п,...,п. (3)

Здесь Ym (0, ф) — нормированные сферические гармоники, радиальные функции Rn (r, Л) = Anjn{v7—Ar) + Bnyn{^f—Xr), где jn (•) и yn (•) — сферические функции Бесселя первого и второго родов, а Л < 0 — собственные значения [12]. Выделяя поля, удовлетворяющие граничному условию OnmldD = 0, получим уравнение на собственные значения:

jn(V—Лп) yn(V—Xr^ — j,n(V—Лг2) yn(V—Xr^ = 0, (4)

которое для каждого n имеет счетное семейство отрицательных решений Xkn.

Таким образом, система собственных полей имеет вид k©m = Rkn (r) Ym (0,ф), k = 0,1,..., n = 0,1,..., m = n,..., n, где Rkn (r) = Aknjn {V—XTnr) + Bknyn{V—Xknr). Эта система ортогональна относительно скалярного произведения (-, •)1, причем каждому собственному значению Xkn соответствуют 2n +1 собственное поле при m = —n,..., n. В дальнейшем будем считать, что коэффициенты Akn и Bkn выбраны так, что функции Rkn нормированы на отрезке [r1; r2] с весом r2. Тогда система полей k©m ортонормальна. Получаем разложение температуры в системе (1):

<х n

© = Е Ekam(t)k©m. (5)

k,n=0 m=-n

Для разложения скорости введем пространство H2 = {u | Vu = 0, u|9D = 0} со скалярным умножением (u, v)2 = J uvdw. Это пространство разлагается в прямую

D

сумму ортогональных подпространств Hp Ф Hp, где Hp и Hp — пространства тороидальных и полоидальных полей соответственно [1]. Тороидальное поле определим как ТФ = rot (Фг), а полоидальное — как Рф = rot rot (Фг), где r — радиус-вектор. При этом [1]

rot Тф = Рф, rot РФ = Т-дф, ДТф = ТДф, ДРф = РДф. (6)

Таким образом, Hp и Hp — инвариантные подпространства оператора Лапласа. Тогда естественно пытаться использовать для разложения скорости собственные поля оператора Лапласа в подпространствах Hp и Hp отдельно.

Для того чтобы поля Тф и Рф были собственными для векторного оператора Лапласа, необходимо и достаточно, чтобы их производящие функции Ф и Ф были собственными для скалярного оператора Лапласа с теми же собственными значениями. Рассмотрим тороидальные собственные поля:

1 dYm dYm т - ЗУГ'

rot (RknY^r) = Rkn~—77 г™ ee — Rkn n ev = Rkn——77 Y- me<? — Rkn n! (7)

sin в дв sin в дв

Поскольку функции Rkn нулевые на границах слоя, то поля (7) принадлежат пространству H2. Введем для них обозначение k T^-

Будем использовать далее разложение скалярного оператора Лапласа на радиальную и угловую части: А = Ar + (1/r2) As, где

1 d ( 2 д \ 1 d ( . В \ 1 д2

Ar = r2— , As = —sin в— +

r r2 dr \ dr j s sin в дв \ дв) sin2 в д' При этом известно, что AsY,¡m = —n (n + 1) Y^•

Полоидальные собственные поля оператора Лапласа

rot rot (Rknvm г):

= _ Rkn Д Vme r

+

dRkn + Rkn

dr r

dYm

n

de

ee +

dRkn + Rkn

dr

dYm

n ,

r J sin e dp

R

kn

n (n +1) Ynmer +

dRkn + Rkn

dYm

n

+

dR

kn

+

R

kn

m

Y-

(8)

r ' ' " \ dr ' r ) дв " ' \ dr ' r ) sin в

не принадлежат пространству H2, поскольку функции Rkn не удовлетворяют граничному условию

dRkn , Rkn

dr

+

r

0.

(9)

r=ri ,Г2

Таким образом, среди собственных полей оператора Лапласа в пространстве H2 нет полоидальных, и использовать их для разложения в Hp невозможно.

Для выбора полоидальных мод скорости заметим, что оператор Лапласа входит в первое уравнение системы (1) в комбинации rot А и поля дT^ удовлетворяют соотношению rot Ад.Tm = Afcnrot kTm. Тогда полоидальные компоненты скорости естественно искать как решения спектральной задачи

rot ДРФ = ^ rot P

ф

(10)

в пространстве Н2.

Покажем, что поле Рф является решением задачи (10) тогда и только тогда, когда функция Ф является решением задачи

Д2Ф = ^ДФ.

(11)

Действительно, пусть Ф удовлетворяет равенству(11). Отметим, что для соленоидаль-

ных

полей Д = -rot2. Тогда, с учетом (6), получим rot ДР^

rot Р

дф

T

T

-Д(ДФ)

-^дф

-^Тдф = -^ДТф = ^rot2Тф = ^rot Р

ф.

Пусть теперь Рф удовлетворяет (10). Тогда rot Рдф_мф = 0. В то же время rot РдФ_мФ = = rot го^дф-^ф = -АТдф-мф = -Тд(дф-^ф) = -rot ((А2Ф - ^АФ) г). Получаем, что (А2Ф — ^АФ) г = VU, где U — произвольное скалярное поле. При этом видно, что U зависит только от радиальной координаты и

. 2 т . т 1dU Д2Ф - ^ДФ = ——.

r dr

(12)

Общее решение неоднородного уравнения (12) представим в виде общего решения Фо однородного уравнения Д2Ф — ^ДФ = 0 и частного решения Ф* уравнения (12), которое можно выбрать зависящим только от г. Таким образом, Ф = Ф0 + Ф* (г), но очевидно, что Рф = Рф0+ф», поэтому можно считать, что Ф = Ф0, т.е. Ф удовлетворяет (11).

Решения задачи (11) имеют вид: ^ (г, А) УП^ (0,ф) , п = 0,1,..., т = —п,...,п, где [13]

К = С1п3п( у/—Цг) + СП уп( у/—Цгг) + Спгп + СП г-(п+1). (13)

Из формулы (8) получаем краевые условия для функции R,

p .

n •

RP (ri) = 0, RP (r2) = 0

dR

dr

dR

r=ri

dr

(14)

Г=Г2

1

p

0

0

образующие однородную линейную систему уравнений относительно коэффициентов СП. Эта система имеет ненулевые решения при условии

jn (v/—ar1)

jn (v/—ar2) d .

dr

dtjn

— ar

r=r 2

d dr

d dr

yn (v/—ar1) yn (V/—ar2)

yn

—ar

yn

—ar

rin rn r-(n+1) ' 1 r-(n+1) '2

nrn-1 - (n +1) r- (n+2)

nrn-1 - (n + 1) r- (n+2)

r=r2

(15)

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

полей. Итак, для разложения скорости в Hp будем использовать поля rot rot (RpnY^7 обозначаемые далее k P^.

Объединяя компоненты разложения в пространствах Ht и Hp, получим представление скорости в системе (1)

v ^ Е (kem (t) kтт+fcYm (t) kpm).

(16)

k,n=0 m=-n

Подставив разложения (5) и (16) в первое и второе уравнения системы (1), получим систему уравнений для амплитуд к(¿), квт ктЛТ Граничные условия и уравнение неразрывности будут обеспечены самим разложением. Отметим, что ортогональны не только тороидальные и полоидальные поля между собой, но и однотипные поля, отличающиеся какими-либо сферическими индексами.

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

0

Г=Г1

3. Динамическая система для амплитуд температуры и скорости

Будем считать, что в разложениях (5) и (16) оставлены какие-либо N членов для скорости и М членов для температуры. Для упрощения записей используем одноиндексные обозначения

М N

е = £ а (*)е,-, V = £ в (*) (17)

j=1 i=1

где vi — поля вида кТ^ и/или к Р^, еj• — поля вида к е^, а (¿) и вi (^ — соответствующие амплитуды. Подставив эти разложения в первое уравнение системы (1), получим

N N N

У^ rot Vi = ввrot (vi x rot Vj) + faairot Vi+ i=1 i,j=1 i=1

M N

+ RaPr-1 rot (Tj er) -T^$rot (ez x v), (18)

j=1 i=1 где ^i — это собственное значение моды vi.

Умножим полученное уравнение (18) скалярно на компоненту завихренности rot Vk и проинтегрируем по объему слоя. Выполняя это для различных k, получим систему обыкновенных дифференциальных уравнений

N N N М N

Е Adt = Е Bjвв + ЕМТЗ + RaPr-^Djaj + т £Ek/Зг, k = 1...N. (19)

i=1 i,j=1 г=1 j=1 г=1

Здесь Ai = (rot vk, rot Vj)2, Bj = (rot vk, rot (v¿ x rot Vj))2, Df = (rot vk, rot 0jer)2, Ef = = — (rot vi, rot (ez x v¿))2. Все эти скалярные произведения несложно могут быть вычислены при заданных радиусах слоя D, причем интегрирование по в и ^ проводится аналитически, а по r — численно.

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

d/ N N М N

-/т = Е B%Згвj + Е Игвг + RaPr-1 J] Dfaj + т J] Ef/г, k = 1...N. (20) i,j=1 г=1 j=1 г=1

Подставив разложение для скорости и температуры во второе уравнение системы (1), умножив почленно на Os и проинтегрировав по слою, получим систему

- N,M N

= Е FjЗга3 + ¿ Я/Зг + Pr-1 (Is + Asas) , s = 1...M. (21)

i,j=1 г=1

Здесь Fj = — (Os, (víV) Qj)1, Я/ = (0s, vгr), Is = — (Os, 2/r)1, As — собственное значение температурной моды Qs.

Объединяя системы (20) и (21), получим замкнутую динамическую систему для амплитуд компонент температуры и скорости:

-в N N M N

-Цт = Е ЩвгРз + Е АР + RaPr-^ DDkaj + т ^ ЕгкР, k = 1-N

i>J=1 i=1 3=1 i=1 (22) - N,M N v 7

--ti = E FijРгаз + ¿ ЯР + Pr-1 (Is + As«s) , s = 1...M.

ij=1 г=1

При выборе мод температуры и скорости в разложениях (17), кроме физических соображений, надо учитывать, что многие из коэффициентов системы (22) могут оказаться нулевыми. Если не принимать это во внимание, то из системы могут исчезнуть физически важные члены. Так, если все коэффициенты Ek = 0, то исчезает кориолисов член, т. е. теряется информация о вращении слоя. Если все Dk = 0, пропадает число Рэлея, содержащее основные физические параметры слоя. К исчезновению гидродинамической нелинейности приводит равенство нулю всех коэффициентов Bj и Fj. Таким

образом, не должны быть полностью нулевыми матрицы Ек и , а также оба массива Вк и Е*

и F . ij ij

4. Конвекция в жидком ядре Земли

Применим рассмотренную модель для описания крупномасштабной конвекции в жидком ядре Земли, пренебрегая его проводимостью и дифференциальным вращением внутреннего ядра. У жидкого ядра т2 = 1391 км и т2 = 3486 км, т.е. к = 2095 км [9]. Тогда безразмерные т2 = 0.664 и т2 = 1.664.

Уравнение (4) решалось численно для п = 0,..., 10, а уравнение (15) — для п = 1,..., 10. При каждом таком п определялись 11 первых корней уравнений, соответствующих к = 0,..., 10. Предварительно корни отделялись графически. Найденные собственные значения Акп и для нескольких низших мод температуры и скорости приведены в таблице.

На рис. 1 и 2 приведены графики нескольких радиальных функций Ккп (т) и Крп (т) соответственно. Видно, что индекс к равен числу нулей функции внутри отрезка [тх; т2].

На рис. 3 изображены несколько линий тока полоидальных мод 0Р3, 1Р1 и 2Р2. Видно, что индекс к соответствует числу линий (к+1), расположенных друг над другом.

Оставим в разложении (17) для скорости моды, отражающие характерные черты движения жидкости в ядре Земли, известные по данным наблюдений. В [6] проанализированы результаты ряда работ по распределению по глубине эрНШ^-функций собственных колебаний Земли. Имеющая максимум на глубинах внешнего ядра эрНШ^-функция сфероидальной моды собственных колебаний ц$4 дает 12 областей. В шести

Собственные значения компонент температуры и скорости

к п ^кп №кп

0 0 —9,8696044010893586188 —

1 0 —39,478417604357434475 —

2 0 —88,826439609804227570 —

3 0 —157, 91367041742973790 —

4 0 —246,74011002723396547 —

0 1 —11,455033134142749864 —38,075147418695209057

1 1 —41,219754268073013367 —81,223365350294919345

2 1 —90,604584355210226459 —156,51998443665601732

3 1 —159,70550740406327457 —239,18054557535862006

4 1 —48,53843939295002972 —353,91379393950900647

0 2 —14,572638397647346606 —36,799344517235429183

1 2 —44,709185001509418657 —82,408073214774901134

2 2 —94,169928226108868868 —155,38531057169492076

3 2 —163,29593636082305466 —240,36767426569833768

4 2 —252,13998293473545468 —352,80364379213322178

0 3 —19,124182167488560428 —37,178684354127621538

1 3 —49,953157609504853354 —84,687302622568106189

2 3 —99,539433622402301062 —155,88308926984687364

3 3 —168,69831014448673504 —242,62778958826604838

4 3 —257, 55451138120985648 —353,31565374005247297

0 4 —24,982958261151985468 —39,867570352259379156

1 4 —56,944257920355348944 —88,414875292654741893

2 4 —106,73522997909339802 —158,49789806819266478

3 4 —175, 93216974699327099 —246,27311073091802174

4 4 —264,79666461286846322 —355,89268966307817916

-1.0

-1.5

-0.5

0.5

1.0

1.5

0

Рис. 1. Графики радиальных функций R¿.n тороидальных мод скорости и температуры: 0 — Roo; 1 — Rio; 2 — R20

Рис. 2. Графики радиальных функций RPn полоидальных мод скорости: 1 — R(1; 2 — R^;

из них плотность вещества выше средней, а в шести других — ниже средней (рис. 4). В качестве возможной интерпретации этого автор работы [6] предлагает 12-ячеистую конвекцию, где в шести областях вещество "тонет", а в шести — "всплывает". Такая конвекция естественным образом связывается с тессеральной сферической гармоникой, имеющей 12-ячеистую структуру максимумов и минимумов.

В нашей модели за перенос вещества в радиальном направлении отвечают радиальные проекции полоидальных мод. Крупномасштабная вертикальная структура конвекции описывается компонентами о Р^, которые могут обеспечить перенос вещества от нижней границы слоя к верхней. Построение карт радиальных проекций 0Рт в плоскости (в, <р) для п =1, 2,10 показало, что искомой 12-ячеистой структурой обладают моды 0Р± (рис. 5). На этом же рисунке изображены линии тока моды 0Р^.

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

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

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

Таким образом, в разложении (17) оставляем две моды скорости: VI — 0Р4 , "2 — 0Р2. По аналогии с параметризацией, используемой при составлении системы Лоренца, для температуры возьмем три моды, обозначив в0 — 0©0, ©1 — 0©4"2, ©2 — 0©2. Компоненты ©1 и ©2 обеспечат разницу температур восходящих и нисходящих потоков, ©0 дает среднее отклонение температуры от линейного по радиусу профиля.

Рассмотрим коэффициенты систем (19) и (21) для отобранных мод. Непосредственное аналитическое интегрирование по в и ^ дало следующие выражения:

A1

A2

r 2

2V Г2

ri

-40r2 R

2 D P

04

d R04

dr2

+ 2r

dR04 ^

dr

+ 4r'

j rP j2-DP dR04 d R04

dr dr2

+ r

, d2Rp4 \

dr2

rlfíP

- 80r<^Rr04 + 400 (RP4)2

dr,

Г2

Di = D2 = -20

R

04 / 2 d2R(°4

r

r

dr2

+ 2r-R04 - 20R04 I dr,

Ei = -E2 = A}/10

(23)

Fi°i = F22 = --0 J R00 (R04RP4 + rdr (R04RP4^ dr

ri

Г2

H° = H22 = 2H rR04 R04dr,

•Ю4 R04d

ri

Г2

I0 = - Vn rR00dr.

ri

Остальные коэффициенты равны нулю. Отметим также, что поля v° и v2 имеют одно и то же собственное значение ^ = -39.8675703523, поля в° и в2 — собственное значение А = -24.9829582612, поле в0 — собственное значение А0 = -9.8696044011. Система (22) в рассматриваемом случае примет вид

= + RaPr-1 (D 1/A 1) + ^

de2 = ^2 + RaPr- 1 (D 1/A 1) «2 - 1T3

da „0 / „ „ ч „ _ 1 (r0

F° (A«1 + в2«2) + Pr-1 (/° + А0 «0)

dt

^ = FoA «0 + Нв + Pr-1Aa1, dt

^ = Fo^ «0 + Н1в2 + Pr-1A«2. dt

Выражение для дивергенции поля скоростей в фазовом пространстве этой системы имеет вид + Рг-1 (А0 + 2А). Оно очевидно отрицательно, т. е. все фазовые траектории сходятся к предельному многообразию, размерность которого меньше 5. Найдем точки покоя системы (24), т.е. решения системы:

+ ИаРг-1 (£1/41) «1 + 10в2 = 0, ^2 + ИаРг-1 (£1/41) «2 — ЮА = 0,

^101 (в1«1 + в2«2) + Рг-1 (I0 + Аоао) = 0, (25)

0^0 + Н1в1 + Рг-1А«1 = 0, ^0в2а0 + Н\^2 + Рг-1А«2 = 0.

Из первых двух уравнений получим

10И,аРг-1£{ (10^ — та2)

в ' А1 (100^2 + т2)

ЮЯа

в = —

ЮИаРг-1^ (10^а2 + та1)

А1 (100^2 + т2)

Подставив эти выражения в четвертое и пятое уравнения системы (25), получим однородную линейную систему уравнений относительно а1 и а2 с основной матрицей:

ШИаЯ^ (^0а0 + Н1) ЮИа^т (^0а0 + Н1)

А — -

А1 (100^2 + т2) А1 (100^2 + т2)

А

ЮИа^т (^¿«0 + Н1) ШИа^ (^а0 + Н1)

А1 (100^2 + т2) А1 (100^2 + т2)

Эта система может иметь ненулевые решения только в том случае, если одновременно

ШИарУ (^¿^0 + Н1) =0 ЮИаР^т (^>0 + Н1) =0 А1 (100^2 + т2) = , А1 (100^2 + т2) = , но эти равенства очевидно несовместны. Таким образом, а1 = а2 = 0, значит, и в1 = в2 = 0. Теперь из третьего уравнения системы (25) получим а0 = —10/А0.

Итак, при любых значениях управляющих параметров система (24) имеет одну точку покоя, соответствующую отсутствию конвекции (V = 0). То, что у системы нет точек покоя с ненулевыми амплитудами компонент скорости, означает, что в данной модели невозможен режим стационарной конвекции.

Рассмотрим вопрос об устойчивости найденной точки. Линеаризуем систему (24) в окрестности точки покоя:

^ = + йаРг-1 И/А!) «1 + Юв2, ^ = ^2 + ИаРг-1 (£1/А1) «2 — ЮА, ^ = Рг-1А0 «0, (26)

^ = (—П1°/А0 + Н1) 01 + Рг-1 Аа1, ^ = (—П^/А, + Н1) 02 + Рг-1 Аа2.

Характеристический многочлен системы (26) имеет вид (Ao/Pr — z) Q (z), где

0

RaPr-1DÍ/ A} 0 '

Pr-1A — z

Аналитически можно определить только одно собственное значение z1 = A/Pr < 0. Дальнейшее исследование устойчивости проведем численно. Вычисление интегралов по формулам (23) дает: A} = 22139.1982314796, D} = 423.6650368945, FÍ0 = 3.7866074003, H} = 16.6933578027, Io = —6'3830764864.

Примем следующие значения для параметров земного ядра [8, 14]: v = 10-6 м2/с, k = 10-5 м2/с, ¿T = 103 К. Известно также, что h =2.1 ■ 106 ми Q = 7.3 ■ 10-5 рад/с. Тогда Pr = 10-1, т = 6, 4 ■ 1014. Значение объемного расширения оставим свободным с ограничением 1 — в^Т >> 0, т. е. в << 2 ■ 10-4 К-1.

Результаты численного исследования корней многочлена Q (z) показывают, что точка покоя теряет устойчивость при превышении коэффициентом в критического значения всг = 10-8'1. Соответствующее критическое значение числа Рэлея Racr = 1026.

Выводы

В работе изучена модель конвекции во вращающемся шаровом слое вязкой жидкости в приближении Буссинеска. Гидродинамические токи разложены на собственные для оператора Лапласа тороидальные компоненты и собственные полоидальные компоненты задачи rot ДР = ^ rot P, что позволяет осуществлять селекцию мод и изучать структуру конвекции на разных пространственных масштабах. Получена квадратично-нелинейная динамическая система для амплитуд компонент скорости и температуры, являющаяся аналогом классической системы Лоренца [2, 5] маломодовой конвекции.

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

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

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

[1] CHANDRASEKHAR S. Hydrodynamics and Hydromagnetic Stability. N.Y.: Dover Publ. Inc., 1981. 654 p.

[2] Монин А.С. Теоретические основы геофизической гидродинамики. Л.: Гидрометеоиздат, 1988. 422 с.

[3] Зельдович Я.Б., РузмАЙкин А.А., Соколов Д.Д. Магнитные поля в астрофизике. М.—Ижевск: НИЦ "Регулярная и хаотическая динамика", 2006. 384 с.

Q (z)

^ — z —т/10

H1 — FWAo 0

т/10 ^ — z 0

H1 — FWAo

RaPr-1 0

Pr-1A — z 0

D1/A1

[4] ПЕдлоски Дж. Геофизическая гидродинамика: Пер. с англ. М.: Мир, 1984. Т. 1. 811 с.

[5] Гледзер Е.Б., ДолжАнский Ф.В., Обухов А.М. Системы гидродинамического типа и их применение. М.: Наука, 1981. 366 с.

[6] Кузнецов В.В. Анизотропия свойств внутреннего ядра Земли // УФН. 1997. Т. 167, № 9. C. 1002-1012.

[7] МолодЕнский С.М. Приливы, нутация и внутреннее строение Земли. М.: Наука, 1984. 215 с.

[8] Джекобс Дж. Земное ядро: Пер. с англ. М.: Мир, 1979. 305 с.

[9] КозЕнко А.В. Земля // Физическая энциклопедия. М.: Советская энциклопедия, 1990. Т. 2. C. 78-80.

[10] Harder H., Hansen U. A finite-volume solution method for thermal convection and dynamo problems in spherical shells // Geophys. J. Int. 2005. Vol. 161. P. 522-532.

[11] Reshetnyak M., Steffen B. A dynamo model in a spherical shell // Numerical Methods and Programming. 2005. Vol. 6. P. 27-32.

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

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

[14] ПАРкинсон У. Введение в геомагнетизм. М.: Мир, 1986. 525 с.

Поступила в редакцию 6 августа 2008 г., в переработанном виде —15 февраля 2009 г.