Моделирование трехмерной конвекции в мантии Земли с применением неявного метода расщепления по физическим процессам

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

Том 11, № 4, 2006

МОДЕЛИРОВАНИЕ ТРЕХМЕРНОЙ КОНВЕКЦИИ В МАНТИИ ЗЕМЛИ С ПРИМЕНЕНИЕМ НЕЯВНОГО МЕТОДА РАСЩЕПЛЕНИЯ ПО ФИЗИЧЕСКИМ ПРОЦЕССАМ

В. В. Червов

Институт геологии и минералогии СО РАН, Новосибирск, Россия

e-mail: chervov@ist.nsc.ru

Based on the method of implicit splitting over the physical processes a numerical model of 3-D convection in the Earth mantle has been constructed. The results of some numerical experiments are presented.

Введение

Исследование конвекции в недрах Земли является одной из центральных задач геофизики. Работы в этом направлении, выполненные в последние годы (см., например, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]), значительно расширили наши представления о строении и составе недр. С развитием возможностей вычислительной техники в 90-х годах стали появляться трехмерные математические модели тепловой конвекции [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Так, можно выделить работы В.П. Трубицына с соавторами [10, 14, 15, 16], охватывающие широкий класс прикладных проблем геотектоники. При решении трехмерных задач гидродинамики используется целый ряд подходов, основанных на применении как уравнений в естественных переменных, так и в векторных переменных ф, ш и V, ш [24, 25, 26, 27, 28, 29, 30].

Подобные подходы могут быть обобщены и на случай задач конвекции. Подход с использованием трехмерных векторов завихренности и потенциала успешно применен, в частности, в работе [24], где с применением уравнений Навье — Стокса в приближении Обербека — Буссинеска рассматривалось течение, вызванное градиентами температуры и концентрации в поле силы тяжести в прямоугольном параллелепипеде. Для решения применялись разнесенные сетки с размещением физических величин в различных местах вычислительной ячейки, что позволило обеспечить консервативность для завихренности на дискретном уровне.

В работе [17] приведены тесты трехмерной конвекции, представленные данные получены в результате применения различных методов. В частности, Кристенсен (Christensen) применил гибридный спектрально-конечно-разностный метод и использовал скалярный потенциал для описания течения; Огава (Ogava) применил конечно-разностный метод для переменных скорость — давление; температура вычислялась неявным конечно-разностным

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

методом. Результаты трехмерного тестирования в переменных ф, ш можно найти в работах [18, 19, 20, 21, 22, 23]. Однако постановка граничных условий в переменных ф, ш для задач, связанных с движением вещества через границы (например, для задачи протекания) сопряжена со значительными трудностями, преодолевая которые приходится отказываться от тождественного соблюдения уравнения неразрывности [27, 28, 29, 30, 31], что сводит на нет преимущества такого подхода.

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

Весьма популярным является в настоящее время метод расщепления по физическим процессам [25, 26, 27, 28, 29].

В настоящей работе продолжено тестирование задачи конвекции, предложенной авторами статьи [17]. В работах [18, 19, 20, 21, 22, 23] с использованием подхода в переменных вектор завихренности — векторный потенциал продемонстрировано хорошее согласие с тестами в [17]. С применением метода расщепления по физическим процессам получены результаты, не уступающие по точности найденным ранее с применением ф, ш-подхода. Решена модельная задача протекания мантийного вещества через вертикальные границы.

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

Для описания течений в верхней мантии Земли привлекается хорошо известная математическая модель, включающая в себя обезразмеренные уравнения [5]:

УУ = 0,

Ур = Е + ИаТ g, дТ

— + УУТ = У2Т. дЬ

Здесь р — давление; Т — температура; Ь — время; И,а

ард;

Я3АТ

число Рэлея; V

ПоХ

вектор скорости; g = (0,0, — д*), д* — ускорение силы тяжести; d — вертикальный размер конвектирующей области; АТ = Ттах — Тт-1П; х — температуропроводность; а — коэффициент теплового расширения; р, ' — характерные плотность и динамическая вязкость; Е — вектор скорости:

^ ^ д ди д /'ди д^ д / ди дт х дх'дх дуП \ду дх) дг'\дг дх ^ д /ди д^ ^ д ду д /ду дт\

у дхП \ ду дх) дуПду дгП \дг ду) , ^ д /ди + дт\ + д /ду + д^\ + ^ д дт * дх'\дг дх) ду'\дг ду) дг'дг,

где и, у, т — компоненты вектора скорости.

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

Система уравнений устроена так, что в начальный момент времени £ = т0 задаются начальные условия лишь для температуры: Т(ж, у, г, то) = Т0(ж,у,г) [32]. Простейшей областью интегрирования является параллелепипед (рис. 1).

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

— на поверхностях ж = 0, ж = X (0 ^ у ^ У, 0 ^ г ^ 1):

дт

и = = я- =0; дж дж

— на поверхностях у = 0, у = У (0 ^ ж ^ X, 0 ^ г ^ 1):

ди дт 0 ду ду

— на поверхностях г = 0, г =1 (0 ^ ж ^ X, 0 ^ у ^ У):

и = V = т = 0.

Для температуры, как и в [17], на боковых гранях ставятся условия теплоизоляции (адиабатическая стенка), т. е. первые производные на вертикальных стенках равны нулю. На верхней и нижней гранях ставятся условия Дирихле: нулевая температура на верхней и некоторая фиксированная температура (в безразмерных уравнениях равная единице) на нижней:

дТ

— на поверхностях ж = 0, ж = X (0 ^ у ^ У, 0 ^ г ^ 1): = 0;

дджТ

— на поверхностях у = 0, у = У (0 ^ ж ^ X, 0 ^ г ^ 1): = 0;

— на поверхности г = 0 (0 ^ ж ^ X, 0 ^ у ^ У): Т = 0;

— на поверхности г = 1 (0 ^ ж ^ X, 0 ^ у ^ У): Т = 1.

Численная модель трехмерных течений в естественных переменных

Для построения численной модели в случае естественных переменных применялся метод установления с использованием неявной схемы расщепления по физическим процессам [25, 26, 27, 28] и метода дробных шагов [33]:

V - V"

-= Е -Урга + Иа Т^;

т

V2

т

где ¿р = рга+1 — р";

уп+1 — V

-= —V

т

дТ

+ V

Численная реализация (1)-(4) включает в себя следующие этапы.

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

2. Из векторного уравнения (1) находится промежуточное поле скорости "V.

3. Рассчитываются поле разности давлений и давление по (2).

4. Из (3) получается окончательное поле скорости.

5. Путем решения (4) вычисляется поле температуры. Процесс повторяется до некоторого значения £ = пт.

В задаче использовалась разнесенная сетка, в которой давление, вязкость и температура определялись в центре элементарного объема, компоненты вектора скорости — в центрах плоскостей ячеек: и — в центре плоскости УХ, V — в центре плоскости XX, т — в центре плоскости XY. В декартовой системе координат равномерную сетку можно задать в виде

ж/ =

гк = (К — 1/2)Дг, Дг > 0, К = 0,1,...,Ь, (Ь — 1)Дг = X,

где Дж, Ду, Дг — размеры шагов сетки; М, N Ь — число ячеек сетки в направлениях ж, у и г.

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

ига+1/3 — и"

- = ЬииП+1/3 + Ь22и" + Ьззи" + Ь21V™ + Ьз1^га — Ь1РП,

т

ип+2/3 — ип+1/3 , / Ч

и-и- = Ь22

ип+1 — и"+2/3

'>22

Ь22 (иП+1 — и"

Здесь использованы стандартные аппроксимации [36]:

Ьц ~ Ь

-'21

д д д д д д П^Г, Ь22 ~ тг-Ьзз ~ —'—,

дх дх ду ду дг дг

д д д д д

Ьз1 ~ , Ь1 ~ •

ду дх дг дх дх

Для у и т конечно-разностное представление аналогично (6).

Для наглядности численной схемы в работе применены индексы, удобные для последующего переноса в компьютерную программу. В центре элементарного объема, где определяются скалярные поля температуры Т, вязкости ' и давления р, индексы — прописные латинские буквы: I, .],К. На гранях — индексы г,], к (рис. 2). Таким образом, между I и I + 1 находится индекс г + 1, между ] и ] + 1 находится индекс ., а между К — 1 и К — индекс к. Например, вместо общепринятой записи 2 ^ к-1 становится удобной запись и^ ,з,к-1, так как в общепринятой записи точка г,], к находится в центре ячейки (рис. 3).

Граничные условия для вычислений разности давлений являются следствием условий для уравнения (9). Так как граничные условия для Уга+1 и "У совпадают, решается задача Неймана с однородными граничными условиями. Метод решения задачи Неймана основан

Рис. 2. Местонахождение индексов по направлению X. В нижней части рисунка — шкала расстояний.

Рис. 3. Сеточный шаблон для метода расщепления. В центре ячейки, кроме давления и температуры, вычисляется и вязкость П1, з,к.

на приеме, предложенном в книге [34]:

О/ДК

/ДК = О/ДК —

МЖЬ

где ^ — область определения индексов; О = ¿р — разность давлений; 5 — индекс внутренних итераций; О5+1 — промежуточное значение разности давлений; произведение МЖЬ — количество ячеек в области.

Численная реализация (3) выполнялась следующим образом:

г

™+1 = -п+1 — тО/^К—

игД,К = "гДК т

„,га+1 = ,,~1га+1 — т

= — т

О/,J,K — О/,J-1,К

^У

к

О/ ДК — О/ ДК-1

Индексы над знаком "минус" в формулах демонстрируют связь "больших" и "малых" индексов.

После вычислений скорости следует расчет температуры либо по схеме стабилизирующей поправки, либо по схеме предиктор-корректор, основанной на схеме стабилизирующей поправки:

Т"+1/6 _ Т"

' +Л22Т + Л33

= ЛТ "+1/6 + Л22Т" + Л33Т

т

Т"+2/6 _ Т"+1/6

-= Л22(Т "+2/6 — Т"),

т

Т га+3/6 — Т га+2/6

Л33(Т "+3/6 — Т"),

т

Т"+1 _ Т"

\ИТ ' ' +Л22Т ' ' +Л33

т

Л11Т"+1/6 + Л22Т"+2/6 + Л33Т"+3/6.

Здесь

д2Т д д2Т д д2Т д

Л11Т ~ ^г — 7Г(и"Т); Л22Т - — — — рт); Л33Т - — — — (т"Т). дж2 дж ду2 ду дг2 дг

Численная реализация вычислений температуры осуществлялась так:

грга+1/3 . т4 / грга+1/3 грга+1/3\ т4 /грга+1/3 г)грга+1/3 . грга+1/3 \

Т/ДК + ^иг+1ДКТ г+1,^К — игДКТгДК ^ —

= Т/ДК — Т^+1,К — ^¿.К^Г^К) + (Т/Д+1,К — 2Т/ДК + Т/Д-1,К) —

у у

— ^ Дк+1Т/Дк+1 — ,J,kТ/Дк) + (Т/ДК+1 — 2Т/ДК + Т/ДК-1),

rHn+2/3 . Tt ( rpn+2/3 гтП+2/з\ Tt (rrn+2/3 nrrn+2/3 . ^n+2/3 ч

1i,j,k + [vi,j+i,KTi,j+1,K - vi,j,KTi,j,K ) - У2

рП+1/3 . Tt ( rpn rpn A Tt (rpn nrpn , rpn \

[I,J,K + h \VI,j+1,K1I,j+1,K - VI,j,K1I,j,K) - fr2 (1I,J+1,K - 21I,J,K + 1I,J-1,K J

rrrn+1 I 4- j rrm+1 rjin+1\ 4 (rpn+1 nrrin+1 I rpn+1 \ _

1 I,J,K + \WI'J'k+11 I,J,k+1 - WI,J,k 1I,J,k) - (1 I,J,K+1 - 21I,J,K + 1 I,J,K —1J =

= Tl+K + JZ {WI,J,k+11I,J,k+1 - WI,J,kTJ) - ^2 (TI',J,K+1 - 21?,.J,K + T'^JK-

Так как температура вычислялась в центрах ячеек (I, К), для интерполяции в точки (г, К) находящихся в центрах граней ячеек привлекались следующие трехточечные выражения (по направлению X):

= т1,з,к + 1(ЗТ/_1д# — Пг^к — ),

о

= Т1,.1,к + 1(ЗТ/+1дк — 2Тгдк — Тг_1^к)-

о

Для направлений У и ^ — аналогичные представления.

На разнесенной сетке, горизонтальный разрез которой представлен на рис. 4, граничные условия прилипания вычислялись следующим образом. Пусть для вектора скорости V на границе Б1 — Б2 заданы условия прилипания, т. е. и = V = = 0. Тогда щ^к = —Щ,о,к и = . Для компоненты V вычисления на границе не производятся, а переносятся

в первую (] = 1) точку массива, определенного в точке как v(I,j,K).

Для представления условий Неймана на границе Б1 — Б 2 можно воспользоваться следующими выражениями. Пусть на Б1 — Б2 поставлены условия ди/ду = дv/дy = в, тогда

ди ду

ui,1,K - ui,0,K 0 ,

— = в ^ Ui,1,K = Ui,0,K + phy,

y=0 hy

dv дУ

-3vi,1,k + 4vi,2,k - VI,3,K = p

2h

y=0 2hy

Определение массива v(:,:,:) в ФОРТРАН-процедуре выглядит так: real, dimension (0:M, 1:N, 0:L) :: v.

Область определения Y-компоненты скорости v:

xI = (I - 1/2)Ax, Ax> 0, I = 0,1,...,M, (M - 1)Ax = X, П H yj = (j - 1)Ay, Ay> 0, j = 1, 2,..., N, (N - 1JAy = Y, zK = (K - 1/2)Az, Az> 0, K = 0,1, ...,L, (L - 1)Az = Z.

Массив u(:,:, :J, в свою очередь, определен внутренними точками и точками, лежащими на границе в направлении X, а также внутренними точками (1 : N -1,1 : K -1) и точками, лежащими за пределами границ в направлениях Y и Z (0, N и 0, K).

Массив u(i,J,K) в ФОРТРАН-процедуре: real, dimension (1:M, 0:N, 0:L) :: u. Область определения X-компоненты скорости:

Xi = (i - 1)Ax, Ax > 0, i = 1, 2,...,M, (M - 1)Ax = X, П H yJ =(J - 1/2)Ay, Ay> 0, j = 0,1,..., N, (N - 1)Ay = Y, zk = (K - 1/2)Az, Az > 0, K = 0,1,...,L, (L - 1)Az = Z.

Рис. 4. Двумерная разнесенная сетка.

Массив w(:,:,:) определен внутренними точками и точками, лежащими на границе в направлении Z, а также внутренними точками (1: M — 1,1: N — 1) и точками, лежащими за пределами границ в направлениях X и Y (0, M и 0, N).

Массив w(I,J,k) в ФОРТРАН-процедуре: real, dimension (0:M, 0:N, 1:L) :: w.

Область определения Z-компоненты скорости:

( xI = (I — 1/2)Дх, Дх > 0, I = 0,1,...,M, (M — 1)Дх = X, Qw =\ yj =(J — 1/2)Ду, Ду> 0, J = 0,1,..., N, (N — 1)Ду = Y, [ zk = (k — 1)Дг, Дг> 0, k = 1, 2,...,L, (L — 1)Дг = Z.

Область определения температуры, давления и вязкости на разнесенной сетке Q приведена выше.

Массивы T(I,J,K), p(I,J,K) и n(I,J,K) в ФОРТРАН-процедуре могут выглядеть так: real, dimension (0:M, 0:N, 0:L) :: Therm, Press, Eta.

2. Результаты расчетов

2.1. Модельная трехмерная задача конвекции в мантии Земли

Тестирование численной модели осуществлялось путем решения модельной задачи [17]. Для постоянной вязкости расчеты проводились при следующих значениях параметров: длина области X = 1.0079, ширина У = 0.6283, высота X = 1.000. Масштабный множитель при вязкости По = 8.0198 • 1023, число Рэлея Иа = (адгр^3ДТ)/(поХ) = 3 • 104.

Решение для переменной вязкости отыскивалось в единичном кубе. При этом задавались следующие параметры: масштабный множитель при вязкости По = 1.201651024, П(Т) = ехр[в/ (Т + ©)] - [в/ (0.5 + 0)], в = [225/ 1п (г)] - 0.251п (г), 0 = 15/ 1п (г) - 0.5;

г = п1т=о / п1т=1 = 20; Иа = (адгР^ДТ)/(пох) = 2 • 104.

Для решения задачи вводилась равномерная в каждом направлении сетка. Вычислялись следующие параметры:

— среднеквадратичная скорость

VTi

\

(ХУZ JJJ (u2 + v2 + w2)dxdydz

A

где A — объем параллелепипеда со сторонами X, У и Z;

— число Нуссельта [35]

Nu = — (XY)-1 УI ddzdxdy,

Stop

где Stop — верхняя поверхность параллелепипеда;

— значения вертикальной компоненты скорости w и температуры T в угловых точках среднего сечения конвективного слоя;

— значения теплового потока = — dT/dz в угловых точках верхней поверхности куба;

Y

Г dT

— интегральный параметр, вычисляемый по формуле т (x,z) = dy вдоль ли-

J dz

о

нии, параллельной оси у, начинающейся точками (0,0.25), (0.5, 0.25), (1, 0.25) фронтальной (^)-плоскости;

— средняя температура Tm = ff Tdxdy, вычисляемая на горизонтальных сечениях

области Sz=0.75 и Sz=0.50, на глубинах z = 0.75 и z = 0.5;

— значение вертикальной компоненты вектора завихренности uiz в точке (0.75, 0.25, 0.75).

— значения z — отклонения свободной поверхности Земли от геоида в угловых точках верхней поверхности куба.

Интегралы вычислялись с применением квадратурной формулы трапеций. Размерные значения (в системе СИ), которые использованы в [17] и в настоящей работе, принимались следующими:

d = 2 700 000, AT = 3700, х = 10-6, a = 10-5, р = 3300, gz = 10.

В качестве начального распределения температуры выбиралось: T(x, у, z, t0) = T0(x, у, z) = (1 — z) + 0.2(cos(nx/Xmax) + cos(ny/ymax)) sin(nz).

Результаты расчетов автора (Che) сопоставлялись с данными Кристенсена (Chr) как наиболее полными из имеющихся в статье [17] (табл. 1 и 2). Относительная ошибка вычислялась по формуле

Err

Che - Chr Chr

100%.

Таблица 1. Постоянная вязкость. Разнесенная сетка. Расщепление в естественных переменных. Консервативная запись всех уравнений. Сравнение решений на двух сетках: 32 х 32 х 31 и 40 х 40 х 39

Параметр Chr на сетке Che на сетке Err, % Che на сетке Err, %

32 x 32 x 64 32 x 32 x 31 40 x 40 x 39

Nu 3.544 3.61954 2.09408 3.59152 1.33025

Vrms 41.00 42.0322 2.45759 41.7962 1.90679

T(0, 0,1/2) 0.8013 0.804283 0.370930 0.803227 0.239961

T(0, Y, 1/2) 0.6188 0.615735 0.492899 0.616842 0.312575

w(0, 0,1/2) 116.6 119.297 2.23934 118.474 1.56013

w(0,Y, 1/2) 40.50 40.2506 0.617252 40.5227 0.05840

т(0, 0) -0.3647 -0.376761 3.19340 -0.374129 2.51234

т (X/2, 0) -0.1292 -0.133626 3.33575 -0.133331 3.12150

т (X, 0) -0.1104 -0.114614 3.69378 -0.114572 3.65906

Tm(3/4) 0.5215 0.522466 0.190724 0.522106 0.121788

Таблица 2. Переменная вязкость. Разнесенная сетка. Расщепление в естественных переменных. Консервативная запись всех уравнений

№ Параметр Chr на сетке 32 x 32 x 64 Che на сетке 40x40 x 39 Err, %

1 Nu 3.03927 3.06762 0.923070

2 Vrms 35.132 35.7431 1.71517

3 w(0, 0,1/2) 165.91 168.643 1.62675

4 w(0, Y, 1/2) -26.72 -26.8058 0.320045

5 w(X,Y, 1/2) -58.23 -59.4940 2.12457

6 T(0, 0,1/2) 0.90529 0.906780 0.164324

7 T(0, Y, 1/2) 0.49565 0.499292 0.729343

8 T(X, Y, 1/2) 0.23925 0.242436 1.31425

9 0(0,0) 5.83390 5.83800 0.0702432

10 0(0,Y) 1.71360 1.72633 0.737409

11 0(X,Y) 0.7684 0.773319 0.636027

12 т (0,1/4) -0.5059 -0.504619 0.253836

13 т (X/2,1/4) -0.1921 -0.196761 2.36874

14 т (X, 1/4) -0.1388 -0.140683 1.33857

15 Tm(3/4) 0.56593 0.570469 0.795663

16 Tm(1/2) 0.58158 0.583488 0.326939

17 z(0, 0,1/2) 10869.0 10577.7 2.75416

18 z(0,Y, 1/2) -4145.00 -3918.71 5.77470

19 z(X,Y, 1/2) -12811.0 -12345.0 3.77509

20 (3/4,1/4, 3/4) -11.125 -10.9630 1.47801

Можно видеть, что результаты расчетов настоящей работы достаточно близки к результатам [17], что свидетельствует о высокой эффективности численной модели.

2.2. Задача протекания

Рассмотрена задача протекания верхнемантийного вещества в кубической области, размеры которой 700 х 700 х 700 км. Число Рэлея, характеризующее режим конвекции, выбрано как И,а = 2.037-105, что также отвечает современным представлениям об условиях в недрах Земли. Основные параметры задачи в системе СИ, пригодные для верхней мантии, полагались следующими:

& = 700 000 м, ДТ = 1800 °С, х = 10-6 м2/с, а = 10-5 °С-1, р = 3300 кг/м3, дх = 10 м/с2, По = 31021 кг/(м • с). Зависимость вязкости от температуры и глубины выражена формулой

Здесь параметры а = 3.89 и Ь = 5.84 обеспечивают перепад вязкости от 20 до 200, что присуще верхнемантийным характеристикам течений: г = Птах/Птт — 20... 200. Скорость на граничной плоскости УХ в точках, где х = 0, равнялась 2 см/год: и|1п = и0 =

Рис. 5. Профиль скорости и = Ух через плоскость у г (х = X = 700 км). На вертикальной оси приведены безразмерные значения скорости.

2 см/год = 312.5 безразмерных единиц; v|in = w|in = 0. На выходной границе YZ в точках, где х = 700 км, задавались "мягкие" условия "свободного" протекания [27, 28]:

= 0.

out

На остальных гранях для вектора скорости задавались нулевые условия прилипания. Условия по температуре (как и в задаче тестирования) в безразмерных единицах:

— на боковых границах дТ/дх = дТ/ду = 0;

— на подошве верхней мантии Т =1;

— на дневной поверхности Т = 0.

Профиль горизонтальной скорости и на выходной границе через t = 50т представлен на рис. 5.

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

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

[1] Алексеев А.С., Лаврентьев М.М., Мухометов Р.Г. и др. Численный метод определения структуры верхней мантии Земли // Мат. пробл. геофизики. 1971. Вып. 2. Новосибирск: ВЦ СО АН СССР. С. 143-165.

[2] Алексеев А.С., Рявой В.З. Модель строения верхней мантии по объемным сейсмическим волнам // Строение земной коры и верхней мантии по данным сейсмических исследований. Киев: Наук. думка, 1977. С. 67-82.

[3] Доврецов Н.Л., Кирдяшкин А.Г. Глубинная геодинамика. Новосибирск: НИЦ ОИГГМ СО РАН, 1994.

[4] Доврецов Н.Л. Пермско-триасовый магматизм и осадконакопление в Евразии как отражение суперплюма // Докл. РАН. 1997. Т. 354. С. 220-223.

[5] Доврецов Н.Л., Кирдяшкин А.Г., Кирдяшкин A.A. Глубинная геодинамика. 2-е изд. Новосибирск: Изд-во СО РАН, 2001. 409 с.

[6] Трувицын В.П., Фрадков А.С. Конвекция под континентами и океанами // Физика Земли. 1985. № 7. С. 3-13.

[7] Трувицын В.П., Бовров А.М., Кувышкин В.В. Влияние континентальной литосферы на структуру мантийной тепловой конвекции // Физика Земли. 1993. № 5. C. 3-11.

[8] Трувицын В.П., Рыков В.В., Трувицын А.П. Конвекция и распределение вязкости в мантии // Физика Земли. 1997. № 3. C. 3-10.

[9] Трувицын В.П., Рыков В.В. Механизм формирования наклонных зон субдукции // Физика Земли. 1997. № 6. C. 3-14.

[10] Трувицын В.П. Основы тектоники плавающих континентов // Физика Земли. 2000. № 9. C. 3-40.

[11] DüBUFFET F., RABINOWICZ M., Monnereaü M. Multiple scales in mantle convection // Earth Planet. Sci. Lett. 2000. Vol. 178. P. 351-366.

ди дх

out

дь дт ду дг

out

д 2v дх2

out

д 2w дх2

[12] Ratcliff J.T., Tackley P.J., Schubert G., Zebib A. Transitions in thermal convection with strongly variable viscosity // Phys. Earth Planet. Intern. 1997. Vol. 102. P. 201-212.

[13] Tackley P.J. Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles //J. Geophys. Res. 1996. Vol. 101. P. 3311-3332.

[14] Рыков В.В., Трубицын В.П. Численное моделирование трехмерной мантийной конвекции и тектоника литосферных плит // Вычисл. сейсмология. 1994. Вып. 26. C. 94-102.

[15] Рыков В.В., Трубицын В.П. Трехмерная модель мантийной конвекции с движущимися континентами // Вычисл. сейсмология. 1994. Вып. 27. C. 21-41.

[16] Trubitsyn V.P., Rykov V.V. A 3D numerical model of the Wilson cycle // J. Geodynam. 1996. Vol. 20. P. 63-75.

[17] Busse F.H., Christensen U., Clever R. et al. 3D convection at infinite Prandl number in Cartesian geometry — a benchmark comparison // Geophiys. Astrophys. Fluid Dynamics. 1993. Vol. 75. P. 39-59.

[18] Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычисл. технологии. 2002. Т. 7, № 1. C. 114-125.

[19] Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением последовательности сеток. // Вычисл. технологии. 2002. Т. 7, № 3. C. 85-92.

[20] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в верхней мантии Земли // Selected Papers of the Intern. Conf. "Fluxes and Structures in Fluids". St. Petersburg, Russia, June 23-26, 2003. M.: IPM RAS, 2004. P. 238—241.

[21] Тычков С.А., Червов В.В., Черных Г.Г. О численном моделировании тепловой конвекции в мантии Земли // Докл. РАН. 2005. Т. 402, № 2. С. 248-254.

[22] Tychkov S.A., Chervov V.V., Chernykh G.G. Numerical modeling of 3D convection in the Earth mantle // Russ. J. Numer. Anal. Modelling. 2005. Vol. 20, N 5. P. 483-500.

[23] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в верхней мантии Земли // Физика Земли. 2005. № 5. С. 48-64.

[24] Бессонов О.А., Брлйловскля В.А., Полежаев В.И. Пространственные эффекты конвекции в расплавах: концентрационные неоднородности, возникновение несимметрии и колебания // Изв. РАН. Механика жидкости и газа. 1997. № 3. С. 74-82.

[25] Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Наука, 1984.

[26] Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л.: Гид-рометеоиздат, 1986.

[27] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. В 2 т. М.: Мир, 1990.

[28] Флетчер К. Вычислительные методы в динамике жидкостей. В 2 т. М.: Мир, 1991.

[29] Voevodin A.F., Goncharova O.N. Method of splitting for physical processes for computation of convection problems in closed domains // Intern. Conf. on Comp. Math. ICCM-2004. 21-25 June, 2004, Novosibirsk, Russia. Books of Proc. P. 942-947.

[30] Aziz K., Hellums J.D. Numerical solution of the three-dimensional equations of motion for laminar natural convection // Phys. Fluids. 1967. Vol. 10, N 2. P. 314-324.

[31] Hirasaki G.J., Hellums J.D. Boundary condition on the vector and scalar potential in viscous three-dimensional hydrodynamics // Quarterly of Appl. Math. 1970. July. P. 293-296.

[32] ФЕДОРЮК М.В. Характеристики течений несжимаемой жидкости в гравитационном поле // Мат. сб. 1988. Т. 137(179), № 4(12).

[33] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[34] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1988.

[35] Blankenbach B., Busse F. et al. A benchmark comparison for mantle convection codes // Geophys. J. 1989. Vol. 98. P. 23-38.

[36] Самарский А.А. Теория разностных схем. М.: Наука, 1977.

Поступила в редакцию 21 февраля 2006 г.