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

Математическая модель движения конвективных ячеек в мантии под действием источника на границе мантия-ядро

А.Н. Боголюбов1,0, А.Н. Грушинский2, Ю. В. Мухартова1, О. Е. Беленькая1

1 Московский государственный университет имени М.В. Ломоносова, физический факультет, кафедра математики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.

2 Институт физики Земли имени О. Ю. Шмидта РАН.

Россия, 123995, Москва, Большая Грузинская ул., д. 10.

E-mail: abogan7@yandex.ru

Статья поступила 30.06.2014, подписана в печать 23.08.2014.

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

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

УДК: 519.63; 536.49. PACS: 02.30.Jr.

Введение

Математическое моделирование конвективного движения вещества в мантии Земли является одной из важных геофизических задач. В настоящее время имеется много работ, посвященных моделированию конвекции в мантии [1-11]. Однако эти работы посвящены квазистационарным процессам (точнее процессам эволюции в установившейся конвективной структуре), которые происходят в постоянных или медленно (даже в геологическом масштабе времени) изменяющихся условиях. Например, известно, что неподвижный континент сначала подавляет мантийную конвекцию под собой и расширяет конвективную ячейку, а затем через определенное время после прогрева субконтинентальной мантии под континентом возникает горячий восходящий мантийный поток. Но, по-видимому, в мантии происходят и более быстрые процессы, связанные с перестройкой структуры мантийной конвекции. Эта перестройка может быть связана, например, с перестройкой конвекции в жидком ядре или проскальзыванием мантии относительно ядра Земли. С такими процессами, по-видимому, связаны перемещения (перескоки) осей спрединга, например перескок оси спредин-га в Атлантическом океане из Лабрадорского моря на его современное положение (Срединный Атлантический хребет в районе Исландия — хребет Ян-Майен).

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

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

жидкостью с постоянной вязкостью, находящейся в вытянутой прямоугольной области толщиной В и длиной Ь с аспектным отношением Ь : В = 10 : 1. Температура на верхней границе фиксирована и равна 0° С, т. е. рассматривается процесс нагрева мантийного вещества снизу. Боковые стенки области полагаются теплоизолированными, поток вещества через боковую границу также отсутствует.

1. Физическая постановка задачи

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

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

2. Математическая модель движения конвективных ячеек в мантии

2.1. Система уравнений, описывающая движение мантии

Тепловая конвекция вязкой мантии описывается распределением скорости V, температуры Т и давления р. Искомые функции подчиняются системе уравнений гид-

родинамики

dVi

dp dSij

= " дХ + дЫ - pgSt3' дР + £ (V) - 0,

dt dx

и уравнению теплопроводности

HL + VidL -JL i^I i -12

dt +V dxi - Oxi[k i - 1,2

С" VkdxJ '

(1) (2)

(3)

где р — плотность вещества, § — ускорение силы тяжести, £ — коэффициент теплопроводности мантии, $¡1 — символ Кронекера, — тензор вязких напряжений:

с _ (V+

= ^ 8x1 + дх'7 , ц — коэффициент вязкости.

Преобразуем рассматриваемую систему уравнений к безразмерному виду. Сделаем линейную замену переменных: х ' = X, ^ ' = £, ^' = -¡£21, где Б — толщина мантии. В новых координатах х-компонента уравнения (1) примет вид

р£ ^ = _± дР + _1 + + д2Уг \

ц Б2 Л ' Ди дх' + ¡2 V дх'2 + дг '2 + дх'дг'), а г-компонента уравнения (1) примет вид р£ 1 с1Уг = Б2 =

1 др дг'

1

+ D

d2Vz J)2VZ z + 2-Z +

d2Vx

Pg V '

х 2 г 2 х г

Для параметров мантии имеет место соотношение « 10~23, поэтому выражениями в левых частях полученных уравнений можно пренебречь.

Будем считать, что плотность вещества мантии имеет вид [10]

р = ро(1 - аТ), где а — некоторая постоянная, а температуру будем искать в виде Т = 5Т•Т', где 5Т — разность температур на нижней и верхней границах расчетной области. Введем обозначения

, Б2

p - kkVp

и представим давление р' в виде суммы р' = р'0+р', где р'0 определяется гидростатическим распределением

kv dPO

D3 dz'

- Pog-

а рр представляет собой отклонение давления от гидростатического. Так как р'0 не зависит от переменной х', то в новых переменных систему уравнений (1)-(3) можно переписать в виде

dx' + dx'2

+

d2Vx + d2V'

dp

z

d2Vz

+2

x

dz '2 d2VZ dz2

+

dr + rdL + rL -

dt ' + x dx' + z dz' Vx Vz

—x +—- - o,

x z

x z

_d!VL

x z 2L

+

= 0,

+ RaT' 2L

= 0,

(4)

dx '2 dz'2'

где Ra - p0agD3ST/kv — число Рэлея.

2.2. Граничные условия

2.2.1. Граничные условия для температуры. Будем считать боковые границы расчетной области теплонепроницаемыми. Это позволяет сформулировать на вертикальных границах условия равенства нулю теплового потока

Т

x

= 0.

х'=0, х' = 10

На верхней свободной границе мантии будем считать температуру равной нулю:

T' I , , - 0.

Iz'-1

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

kdT

z

- A(T0 - T)|z-0 + q,

z=0

где А — коэффициент теплообмена на границе мантия-ядро, или же, в новых переменных:

dT' z

- А '(1 - T')|z'-0 + q'(x', t'),

(5)

z' =0

где

= D , q = D q.

А- k •А- q- klL •q-

Если А' достаточно велико, то при отсутствии источника граничное условие (5) близко к условию Дирихле

Т' |г'=0 = 1,

используемому в классической постановке задачи.

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

Vz Vx

Vz|z'-0 - 0- Szx|z'-0 - ^ + d#

откуда в силу равенства (6) получаем Vx I

- 0,

z'-0

z

- 0.

(6)

(7)

z'-0

Условия (6), (7) также имеют место и на верхней границе, соответствующей г = 1 . На боковых границах области получаем аналогичные условия:

Vx I - 0,

xlx'-0, x'-10

dV

dx'

- 0.

(8)

x'-0, x'-10

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

ш - rot V'.

Поскольку мы рассматриваем двумерную задачу, то ш = еуш, где

dv' dv:

Oj__л __л

dz' дх''

В рассматриваемом приближении несжимаемой жидкости справедливо равенство

div V' = О,

поэтому можно ввести векторный потенциал скоростей ф, так что V' = rot ф. В случае двумерной задачи ф = еуф, причем

dz>' * дхг Подставляя эти соотношения в выражение для функции lu , получаем уравнение

Аф = -ш.

Исключая избыточное давление р' из системы уравнений (4), получаем

аг

Aco = Ra

дхг

z'l

0.8 0.6 0.4 0.2 О

Температура Т'

i L

0 1 2 3 4 5 6 7 8 9 y'10 Положение источника

■:; \.......-

I |

■ г

Уравнение теплопроводности может быть записано в виде

ат' дф ат' дф ат' _ АТ,

дР дг' дх' дх' дг' Граничные условия для функций ф и ш получим из условий (6)-(8):

ф\?'=о, 2'=\ = о, ш\?>=0, 2>=\ = о,

Ф|а-'=0, а-' = Ю = 0, ш\х'=0, а-' = Ю = 0.

3. Численный эксперимент

Расчеты конвекции в мантии были проведены для нескольких значений величины мощности источни-

Температура Т'

1.5 1

¡0.5

о

Положение источника

0.8| 0.( О/ 0.2| о!

2 3 4 5 6 7 8 Векторный потенциал скорости

Векторный потенциал скорости

I_

1

Температура Т'

|| щ |Г к

'0 1 2 3 4 5 6 7 8 9 10 Положение источника

г , , j V.....■

'0 1 2 3 4 5 6 7 8 9 10 Векторный потенциал скорости

6 1 4 1 2 " п И

I

1

1.5

J0.5

■

" I

и.

|20 10

z'l

0.8 0.6 0.4 0.2 О

I

10000 5000 Ч' О

z' 1

0.8

0.6 0.4 0.2

Температура Т'

ь.

012345678 Положение источника 9 Л.' 10

, j .....

01234 5 678 Векторный потенциал скорости 9 Ю

1 II •

| J:

1.5

■ 0.5

I

1

1 1с

Г"

1-20

Рис. 1. Эволюция структуры конвективных ячеек в случае /?а = 10000, с = 0.22- 10 8 м/с

ка дтах (^тах ~ 0.4 Вт/м2 и дтах ~ 0.2 Вт/м2) и для нескольких значений скорости с его движения (0.22- Ю-8, 0.55-Ю-8, 1.Ы0-8 и 5.5 • 10~8 м/с). Основные параметры системы: В = 29 • 105 м — толщина мантии, £ = 32- 10~3 Вт/(м-К) — коэффициент теплопроводности мантии, 6Т = 3700 К — разность температур на нижней и верхней границах мантии. При этом нормировочные множители, связывающие время I, компоненты скорости и мощность ц источника тепла на границе с соответствующими безразмерными величинами, равны

t Г)2

- = — = 26.28125 • 1013,

V к

/г ¿7

= 40.82- 10

-6

Температура Т'

Задача рассматривалась для двух значений числа Рэлея Яа = 10000 и 30000.

Коэффициент теплообмена на границе мантия-ядро был взят равным А' = 10000, или же в размерных единицах А = « 1.1 • 10~4 Вт/(м-К). Размер расчетной области был выбран равным Ь х Я = 10£> х В = 29000 х 2900 км.

В зависимости от величины числа Рэлея Яа изменяется горизонтальный размер конвекционной ячейки. При значении Яа = 10000 в расчетной области помещается 5 ячеек, т. е. горизонтальный размер ячейки равен 5800 км, а при значении =30000 — 8 ячеек, что соответствует среднему горизонтальному размеру ячейки, равному 3625 км. По имеющимся оценкам, горизонтальный размер ячейки реально равен примерно 2000 км.

Поставленная задача решалась численно с помощью

Температура Т'

3 4 5 6 7

Положение источника

3 4 5 6 7 8 Положение источника

2 3 4 5 6 7

Векторный потенциал скорости

0.8 0.6 0.4 0.2 о _

□ 0 1

2 3 4 5 6 7 8 Векторный потенциал скорости

0.8 0.6 0.4 0 2 0.

I

10

У:

|20 10 о -10 ■20

Темпепятот!я Т'

Темпепатути Т'

Рис. 2. Эволюция структуры конвективных ячеек в случае /?а = 10000, с = 0.55- 10 8 м/с

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

Мощность ц'{х', Р) движущегося источника была взята вида

' V

A(t') = ^Ao~, если г е (0' ^тах)'

*о(0 =

Ио.

если t' > t„

если

Xmin + ^start) ' С'>

start/

если ^tart < t' < £itart + — x'min)/c' если t' ^ ^tart + (x'mas - Xmin)/c',

Температура T'

гДе ^тах ~~ момент времени в новых переменных, когда мощность источника достигает своего максимального значения Ло, — момент начала движения источника от точки с координатой х'т-т до точки с координатой х'тт со скоростью с', также выраженной в новых переменных.

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

Опишем некоторые из получившихся результатов. Первый случай: = 10000, дтах « 0.4 Вт/м2, с = 0.22-10~8 м/с. При таких параметрах более разогретый восходящий поток как бы приклеивается к движущемуся источнику и движется вместе с ним. При этом левая ячейка как бы растягивается, что поддерживается дополнительным разогревом от источника (нисходя-

Температура Т'

ill

Положение источника

Векторный потенциал скорости

Температура Т'

Темпепятлтя Т'

3 4 5 6 7

Положение источника

0.8 0.6 0.4 0.2 0

Векторный потенциал скорости

1 П ш

Рис. 3. Эволюция структуры конвективных ячеек в случае Яа= 10000, с = 1.1 • 10 8 м/с

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

Второй случай: Яа = 10000, дтах ~ 0.4 Вт/м2, с = 0.55 • 10~8 м/с. При такой скорости движения и мощности источника, сначала источник тянет за собой восходящий поток, но затем отрывается, успев только лишь немного расширить левую от него конвекционную ячейку и слегка сузить правую. Затем структура конвекционных ячеек по мере остывания почти возвращается к исходной (рис. 2).

Третий случай: Яа = 10000, дтах ~ 0.4 Вт/м2, с = 1.1 • 10~8м/с. При такой скорости картина более сложная. Сначала восходящий поток следует за источ-

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

Четвертый случай: Яа = 10000, дтах « 0.4 Вт/м2, с = 5.5 • 10~8 м/с. При такой скорости движение источника практически не влияет на структуру конвекции.

Рассмотрим теперь вариант с числом Рэлея Яа = 30000. Здесь ситуация гораздо сложнее. Перестройка структуры конвекции оказывается гораздо более существенной. Если для числа Рэлея Яа = 10000 перестройка структуры конвекции в основном сводится

т„, Г"

Рис. 4. Эволюция структуры конвективных ячеек в случае йа = 30000, с = 0.22- 10 8 м/с

к перемещениям восходящих и нисходящих потоков, а количество конвекционных ячеек по большей части остается неизменным, то для Яа = 30000 происходит перестройка с уменьшением числа конвекционных ячеек. Это, по-видимому, связано с большей неустойчивостью структуры и склонностью к возникновению турбулентностей.

Пятый случай: = 30000, дтах « 0.4 Вт/м2, с = 0.22 • 10~8 м/с. При таких параметрах в самом начале движения источника происходит разогрев с возникновением турбулентностей в конвективных ячейках, прилегающих к восходящему потоку, связанному с источником, из-за чего начинается расширение этих ячеек. В результате расширяющиеся ячейки поглощают ближайшие к ним восходящие потоки, и в конвективной структуре остается пять ячеек вместо восьми (рис. 4).

Температура Т'

0.8 0.6 0.4 0.2 о

1 2 3 4 5 6 7 8 Положение источника

г'-1 г 0.8

0.6

0.4

0.2

О -

1 2 3 4 5 6 7 8 9 10 Векторный потенциал скорости

Шестой случай: Яа = 30000, дтах ~ 0.4 Вт/м2, с = 0.55 • 10~8 м/с. При этой скорости источника сначала картина подобна предыдущему случаю, затем источник пытается потянуть за собой восходящий поток, но его голова остается на месте, из-за чего искажается форма восходящего потока. Из-за дополнительного разогрева у его основания возникают зародыши дополнительных восходящих потоков, т. е. сильно искажается форма этих ячеек и распределение тепла внутри них. Две ячейки, испытывающие деформацию, в итоге занимают три четверти площади расчетной области, заставляя оставшиеся три ячейки сузиться, а восходящий поток деформирующихся ячеек сильно смещается вправо (рис. 5).

Седьмой случай: Яа = 30000, дтах ~ 0.4 Вт/м2, с = 1.1 • 10~8 м/с. При такой скорости движения источника сначала все происходит так же, как в предыдущем

0.81 0.6 од| 0.2 0

Температура Т'

1 Л 1Х>

0 1 2 3 4 5 6 7 8 Щ х> 10 Положение источника

. . 7 \

1 2 3 4 5 6 7 а

ЕЖШ1!_

123456789 Векторный потенциал скорости

Температура Т'

Т^л лтг р» г> ат\ го а Т"

3 4 5 6 7

Положение источника

0.8 0.6 0.4 0.2 0

2 3 4 5 6 7 8 9 т' 10 Векторный потенциал скорости

□ и\

Рис. 5. Эволюция структуры конвективных ячеек в случае /?а = 30000, с = 0.55- 10 8 м/с

Температура Т

Температура Т'

2 3 4 5 6 7

Положение источника

0.8 0.6 0.4 0.2

2 3 4 5 6 7 8

Векторный потенциал скорости

ГОТГО

И*

9.2 МС

U:

|40 20 0 -20 -40

_

Положение источника

Векторный потенциал скорости

MINI

Темпепатупа Т'

Темпепят\тя Т'

Рис. 6. Эволюция структуры конвективных ячеек в случае /?а = 30000, с = 1.1 • 10 8 м/с

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

Восьмой случай: Яа = 30000, дтах ~ 0.4 Вт/м2, с = 5.5 • 10~8 м/с. При таких параметрах конвективная структура эволюционирует аналогично предыдущему варианту.

Заключение

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

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

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

1. Christensen и. // ТескторЬузюэ. 1983. 95. Р. 1.

2. Christensen U.R. // Geophys. J. Astr. Soc. 1984. 77. P. 343.

3. Christensen U.R., Yuen D.A. // J. Geophys. Res. 1984. 89. P. 4389.

4. Glatzmaier G.A. // Geophys. astrophys. Fluid Dyn. 1988. 43. P. 223.

5. Lenardic A.,Kaula W.M. // J. Geophys. Res. 1994. 99. P. 15697.

6. MachetelP., Weber P. // Nature. 1991. 350. P. 55.

7. Parmentier E.M. Sotin C., Travis B.J. // Geophys. J. Int. 1994. 116. P. 241.

8. Solheim L.P., Peltier W.R. // J. Geophys. Res. 1994. 99. P. 15861.

9. Tackley P.J., Stevenson D.J., Glatzmaier G.A., Schubert G. // J. Geophys. Res. 1994. 99. P. 15877.

10. Трубицын В.П., Рыков В.В. // Росс. журн. наук о Земле. 1998. 1, №1. С. 1.

11. Червов В.В. // Вычисл. технологии. 2002. 7, № 1. С. 114.

12. Самарский А.А. Теория разностных схем. М., 1983.

A mathematical model of convective cell motion in the mantle under exposure to a heat source at the mantle-core boundary

A.N. Bogolyubov1,a, A.N. Grushinsky2, Yu. V. Mukhartova1, O.E. Belenkaya1

1 Department of Mathematics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

2 O. Yu. Shmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow 123995, Russia. E-mail: a bogan7@yandex.ru.

An algorithm for numerically solving the vector problem on the modeling of the convective motion of matter in the Earth's mantle was built and implemented. A 2D model of convective cell formation in the mantle was considered with a moving heat source at the boundary with a core that generates an upwelling flow of the mantle substance. The motion of the heat source is possible only along the core-mantle boundary. If the heat source is absent, the system is conservative. The mantle was modeled as an incompressible fluid with constant viscosity.

Keywords: mantle of the Earth, convective flow of the mantle substance, moving heat source, hydrodynamics equations, vector potential, vorticity vector, finite difference method. PACS: 02.30.Jr. Received 30 June 2014.

English version: Moscow University Physics Bulletin 6(2014).

Сведения об авторах

1. Боголюбов Александр Николаевич — докт. физ.-мат. наук, профессор; тел.: (495) 939-10-33; e-mail: bogan7@yandex.ru. 3. Грушинский Андрей Николаевич — канд. физ.-мат. наук, ст. науч. сотрудник; тел.: (495) 939-10-33; e-mail: a.grushinsky@mail.ru.

3. Мухартова Юлия Вячеславовна — канд. физ.-мат. наук, ст. науч. сотрудник; e-mail: muhartova@yandex.ru.

4. Беленькая Ольга Евгеньевна — студентка; e-mail: waity@mail.ru.