Применение R-функций для решения задачи конвекции в мантии Земли Текст научной статьи по специальности «Математика»

СТАТЬИ

ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА

Применение R-функций для решения задачи конвекции

в мантии Земли

А. Н. Боголюбов1,0, А.Н. Грушинский2, М. И. Светкин1,6

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

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

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

E-mail: a [email protected], 6 [email protected] Статья поступила 06.11.2016, подписана в печать 28.11.2016.

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

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

УДК: 519.63; 536.49. PACS: 02.70.Dh, 91.35.Dc.

Введение

Изучение внутреннего строения Земли и происходящих внутри нее процессов — один из важнейших вопросов геофизики. Численное моделирование позволило за последнее время многое узнать об этих процессах и их влиянии на движение тектонических плит.

Существует большое количество работ (например, [1-4]), посвященных численному моделированию процессов конвекции в мантии в двух и трехмерных случаях в декартовых и сферических координатах. Для решения данного класса задач конвекции хорошо зарекомендовало себя использование переменных типа «завихренности» и векторного потенциала.

В настоящей работе для решения возникающих краевых задач был применен метод Я-функ-ций В. Л. Рвачева [5]. Достоинства этого метода заключается в том, что он позволяет точно аппроксимировать граничные условия различных типов в задачах со сложной геометрией. Идея применения этого аппарата состоит в построении множества (пространства) функций, точно удовлетворяющих нужным краевым условиям. В дальнейшем решение задачи ищется в этом множестве функций методом Ритца, Галеркина, коллокаций или другими аналогичными методами. Показано, что данный аппарат позволяет учесть широкий класс граничных условий в задачах с нестандартной геометрией и может быть применен для решения краевых задач, возникающих в различных областях физики [6-9].

1. Постановка задачи

В Институте физики Земли имени О. Ю. Шмидта РАН возникла гипотеза о процессах, происходящих в мантии Земли в районе моря Скоша, находящемся на границе Южного и Атлантического океанов между Южной Америкой и Антарктидой [10]. Она заключается в том, что аномалия формы геоида в этой области возникла из-за мощного плюма — восходящего потока разогретого вещества из нижней мантии и с границы ядро-мантия. Существующий характер поля высот геоида наталкивает на предположение, что плюм, который сейчас расположен под тройным сочленением Буве [11], мигрировал туда из района западной части моря Скоша, где ранее активный рифтогенез сейчас затухает. При этом активные процессы в недрах этого региона протекают теперь в менее глубоких слоях недр Земли, а именно в земной коре и верхней мантии, лишь слегка захватывая нижнюю мантию.

По-видимому, в мантии происходят и более быстрые процессы, связанные с перестройкой структуры мантийной конвекции. Эта перестройка может быть связана, например, с перестройкой конвекции в жидком ядре или с проскальзыванием мантии относительно ядра Земли.

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

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

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

2. Математическая модель

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

Мантия моделируется несжимаемой жидкостью с постоянной вязкостью и коэффициентом теплопроводности, находящейся в вытянутой прямоугольной области толщиной В и длиной Ь. Температура на нижней и верхней границах фиксирована и соответственно равна Т = Т0 и Т = 0, что соответствует нагреву снизу. Боковые стенки полагаются теплоизолированными. Континент выбирается в виде легкой, твердой, толстой, теплопроводной прямоугольной плиты длиной I и толщиной ( + (0, плавающей в мантии, где ( — глубина погружения в мантию и (0 — высота континента над мантией (рис. 1).

D

4l d

L х

Рис. 1. Геометрия модели

Запишем исходную систему уравнений гидродинамики для свободной конвекции в несжимаемой жидкости, состоящую из уравнений несжимаемости, Навье-Стокса и переноса тепла [2, 13]:

div V = 0,

0 = -Vp + AV + Ra ezT,

dL + (V, VT ) = kAT,

dt

(1) (2)

(3)

где V, T, p — неизвестные поля скоростей, температуры и давления соответственно, Яа = = gвT0D3/(kv) — число Рэлея [13], V — кинематиче-

ская вязкость мантии, р — плотность мантии, в — коэффициент температурного расширения мантии, k — коэффициент температуропроводности мантии, g — ускорение свободного падения.

Для описания конвекции хорошо зарекомендовало себя использование переменных вектора «завихренности» w = rot V и векторного потенциала ф , V = rot ф [1]. В таких переменных уравнение несжимаемости (1) выполняется автоматически, а из второго можно исключить неизвестное давление p. Поскольку рассматривается двумерная постановка задачи, то в этих векторах можно выбрать только одну ненулевую компоненту: w = wey, ф = феу. В новых переменных система (1)-(3) примет вид

Аф = -w,

г, dT

Aw = Ra -7—, ox

dL - d + d () = AT

dt dz V dx / dx '

(ф£)

дф VX — >

dz

w=

dz

Vz = ^, dx

dx

(4)

(5)

(6)

(7)

(8)

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

1 £+/

,,du Mdt =

[p\x=z - p\x=+] dz +

dz

z=1-d

1-d

dt = U.

dx,

(9) (10)

Температура внутри континента Tc описывается уравнением теплопроводности

dTr dTc

-dt + U dx = KATc,

(11)

где К — коэффициент температуропроводности континента.

Данная система дифференциальных уравнений в частных производных (4)-(6), (9), (10) дополняется соответствующим количеством начальных и граничных условий.

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

Т|г=1 = 0, Т|г=0 = q(t, X) (12)

— на нижней границе мантии г = 0 возможно наличие некоторого нестационарного источника тепла q(t, х);

тт=о

dx

условие теплоизоляция торцов мантии;

т = т — = K9^ dn dn

(13)

— условия непрерывности температуры и потока тепла на границе мантия-континент, п — нормаль к этой границе;

Т = 0 (15)

на границе континента, не соприкасающейся с мантией.

Для скорости вещества в мантии ставятся условия непротекания и непроскальзывания через границу, а на границе с континентом — непротекания и прилипания. Условия для векторного потенциала и вектора «завихренности» получаются из их определения (7)-(8) и условия неразрывности (1):

ф = 0, ш = 0 (16)

на границе мантии, не соприкасающейся с континентом. На торцах континента условия несколько сложнее:

ф = и(1 - г), ш = -0. (17)

На нижней границе из условий для скорости получаются два условия на векторный потенциал:

Ф = и(1 - й), дф = -и, ш = -0. (18)

На границе мантия-континент для вектора «завихренности» не получается поставить какие-либо конкретные условия, поэтому для него граничные условия ставятся из соображений продолжения оператора задачи (4) на границу.

В начальный момент времени континент неподвижен, его температура равна 0, в мантии отсутствуют потоки, а температура линейно растет с повышением глубины:

Т(х, г, t=0) = 1 - г, ТаЦ=0)= 0, V(t=0) = 0.

(19)

В итоге математическая задача сводится к следующему. Имеются шесть неизвестных функций. Три функции координат и времени для мантийной конвекции: векторный потенциал ф(х, г, 0, вектор завихренности ш(х,г, распределение температуры Т(х, г, три неизвестные функции для континента: распределение температуры Тс(х, г, мгновенная скорость поступательного движения континента как целого ии координата его левого края Для нахождения этих функций имеется

замкнутая система шести взаимосвязанных уравнений (4)-(6), (9)—(11). Эта система дополняется граничными и начальными условиями (12)-(19), позволяющими определить ее решение.

3. Метод R-функций

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

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

В 1968 г. В. Л. Рвачев в своей книге «Геометрические приложения алгебры логики» [14] предложил конструктивное решение этой задачи с помощью Я-функций.

Введем понятие Я-функции. Пусть задано некое разбиение 5 числовой прямой К на £ подмножеств 5: К ^ В£, В£ = {0,1,2,..., £-1}. Эти подмножества будем называть качественными градациями или качествами. Введем также отображение 5": X" ^ Б", 5"(х) = (5(х1),..., 5(хп)), х = (х1,..., х„).

Функция [ называется Я-функцией, если существует функция Р: Б" ^ Вт такая, что 5т о / = Р о 5", т.е. качества значений функции могут быть однозначно определены через качества аргументов с помощью функции £-значной логики Р, называемой сопровождающей функцией.

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

^апё(х, У)= х + у - Vх2 + У2, (20)

Яог(х, у)= х + у +^/х2 + у2, (21)

Лп<й(х) = -х. (22)

Как отмечалось выше, можно использовать другие варианты вместо функций (20)-(22), чтобы получить какие-либо дополнительные свойства (например, гладкость).

В [5] показано, что с помощью данных функций можно строить уравнения ш(х, у) ограниченных областей с кусочно-гладкой границей, являющихся пересечением и/или объединением простых областей. Причем ш(х, у) > 0 во внутренней части области П, ш(х, у) < 0 во внешней части области П и ш(х, у) = 0 на ее границе дП. Функция ш(х, у) не единственная, описывающая данную область, поскольку для произвольной гладкой функции Ф(х, у) > 0 функция а(х, у) = ш(х, у)Ф(х, у) также будет уравнением области, удовлетворяющим тем же условиям. Этой свободой можно воспользоваться, чтобы получить уравнение области, удовлетворяющее некоторым

дополнительным условиям, например условию нор-мализованности до k-го порядка:

дш дn

дП

= 1,

д"

д^

дП

= 0, т = 2,...,

(23)

(26)

и1хедП = м(х), ^П

хедП

= Х(х),

(28)

и(х) =

N

£

¿=1

и (х) ш.' (х)

N

¿=1 ш'(х)

(29)

1

где . — порядок дифференциального оператора В'. Эта структура определена всюду внутри рассматриваемой области.

Чтобы найти неизвестную компоненту Ф, представим ее в виде разложения

N

где п — вектор внешней нормали к границе области дП.

Рассмотрим краевую задачу

Аи = I, х е П, (24)

¿¿и = щ, х е дП', дП ^у дП'. (25)

'

Тогда, зная уравнение области ш(х, у) в представленном выше виде, можно построить оператор В, результат действия которого на произвольную достаточно гладкую функцию Ф будет удовлетворять всем или части граничных условий (25). В первом случае такой оператор называется полной структурой решения, а во втором — неполной. Таким образом, если сделать переход к новой неизвестной функции Ф как и = В[Ф], то все или часть граничных условий будут выполнены и останется произвести выбор функции Ф, чтобы удовлетворялось основное уравнение (24) и, если есть, часть неучтенных граничных условий [5, 8].

Приведем операторы для граничных условий, возникающих в рассматриваемой задаче (4)-(6), (9)-(19). Для граничных условий Дирихле структура выглядит проще всего:

и(х)1хедВ = щ(х), В[Ф] = щ(х)+ ш(х)Ф(х).

Оператор для задачи Неймана требует дополнительного условия нормализованности до первого порядка на границе области дП: ди

ШхедВ = ^ (27)

В[Ф] = Ф + — (УФ, Уш)).

Также в рассматриваемой задаче есть граничное условие вида (18), для него также можно получить структуру решения:

Ф(х) « ФN(х) = Спфп(х)

(30)

п=1

В[Ф] = щ + ш(\ - (Ущ, Уш)) + ш2Ф.

Если граница области дП состоит из гладких частей дП', для которых стоят различные типы граничных условий и построены соответствующие операторы и' = В' [Ф] и уравнения элементов границ Ш', то можно получить общую структуру решения с помощью формулы склейки:

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

Подставив разложение (30) в построенную структуру, получим представление решения в виде и = В(х, у, С1,..., Су). Неизвестные коэффициенты разложения ищутся из условия ортогональности невязки уравнения (24) 5(х, у, С1,..., Су) = = АВ[и] — [ к системе функций {фп}:

(5(х, Сь ..., Су), фп(х)) =

5(х, С1,..., Су )фп(х) (х = 0. (31)

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

Применим полученные методы для построения структуры решения исходной задачи (4)-(6), (9)-(19). Поскольку геометрия задачи не фиксирована, то искомая структура будет зависеть от положения континента и его скорости движения. Помимо этой неявной зависимости от времени, из-за нестационарного граничного условия подвижного источника тепла на нижней границе (12) структура будет зависеть от времени. Таким образом, с помощью методов и операторов (26)-(29), описанных в предыдущем пункте, было построено трехпарамет-рическое семейство структур решений для неизвестных полей:

Т = Вт(£ и, 0[Фт], ф = Вф(£, и, t)[Ф^] & = Вш(£, и, 0[ФШ].

(32)

Для решения начально-краевой задачи была введена равномерная сетка по времени. Определение неизвестных величин происходило в три этапа. Сначала из уравнения (10) находилось новое положение континента. Затем вычислялись распределения неизвестных величин с помощью построенных операторов (32), в которые подставлялись разложения неизвестных компонент (30) по системе двумерных финитных сплайнов и определения коэффициентов разложения с помощью метода Га-леркина (31). По найденным значениям векторного потенциала и вектора «завихренности» из их определения (7), (8) восстанавливались значения скоростей потоков, а затем и давления из уравнения (2). Получив значения всех полей, вычисляли новую скорость континента из уравнения (9).

ш

х

Рис. 2. Начальное состояние мантии и континента

х

Рис. 3. Состояние системы в момент времени t = 10т. Появление конвекционной структуры с восходящим потоком под континентом

х

Рис. 4. Состояние системы в момент времени t = 40т

х

Рис. 5. Состояние системы в момент времени t = 70т. Сформировавшаяся конвекционная структура

Начальное состояние для задачи с континентом представлено на рис. 2. В мантии отсутствуют потоки, распределение температуры равномерное. Континент не прогрет Tc = 0. Вначале под континентом возникает нисходящий поток, затем, спустя некоторое число шагов, континент прогревается. Теплопроводность континента выбрана существенно меньше теплопроводности мантии, поэтому под ним возникает восходящий поток, видный на рис. 3.

Этот поток порождает две первичные конвекционные ячейки слева и справа от континента, затем в силу вязкости мантии возникают вторичные конвекционные ячейки (рис. 4, 5).

Заключение

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

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

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

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

3. Трубицын В.П. // Физика Земли. 2008. № 12. С. 83.

4. Черных Г.Г. и др. // Вычислит. технологии. 2014. 19, № 5. С. 101.

5. Рвачев В.Л. Теория R-функций и некоторые ее приложения. Киев: Наукова думка, 1982.

6. Кравченко В.Ф., Рвачев В.Л. Алгебра логики, атомарные функции и вейвлеты в физических приложениях. М.: Физматлит, 2006.

7. Кравченко В.Ф., Басараб М.А. Булева алгебра и методы аппроксимации в краевых задачах электродинамики. М.: Физматлит, 2004.

8. Максименко-Шейко К.В. R-функции в математическом моделировании геометрических объектов и физических полей. Харьков: ИПМаш НАН Украины, 2009.

9. Басараб М.А. // Математическое моделирование и численные методы. 2014. 1, № 1. С. 18.

10. Удинцев Г.Б., Береснев А.Ф., Кольцова А.В. и др. // Украинский антарктический журнал. 2011-2012. № 10-11. С. 26.

11. Разницин Ю.Н. и др. // ДАН. 1995. 342, № 3. С. 354.

12. Хаин В.Е., Ломизе М.Г. Геотектоника с основами геодинамики. М.: КДУ, 2005.

13. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. 6. Гидродинамика. М.: Наука, 1986.

14. Рвачев В.Л. Геометрические приложения алгебры логики. Киев: Техника, 1967.

Application of R-functions for the solution of the problem of convection in the mantle of the Earth A.N. Bogolyubovu, A.N. Groushinsky2, M.I. Svetkin1b

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

2 The Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences. Moscow 123242, Russia.

E-mail: a [email protected], b [email protected].

In this paper, the process of convection in the Earth's mantle in the presence of a floating continent is considered. The model is a two-dimensional rectangular region of viscous thermally conducting fluid that obeys the equations of hydrodynamics. The method of Rvachev R-functions is used for the description of the problem geometry and the boundary conditions.

Keywords: Earth's mantle, convection, vorticity vector, vector potential, R-functions, Galerkin method. PACS: 02.70.Dh, 91.35.Dc.

Received 6 November 2016.

English version: Moscow University Physics Bulletin. 2017. 72, No. 6. Pp. 507-512.

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

1. Боголюбов Александр Николаевич — доктор физ.-мат. наук, профессор; e-mail: [email protected].

2. Грушинский Андрей Николаевич — канд. физ.-мат. наук, ст. науч. сотрудник; e-mail: [email protected].

3. Светкин Михаил Игоревич — аспирант; e-mail: [email protected].