CC BY
47
Вычислительные технологии
Том 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"
-= Е -Урга + Иа Т^; (1)
т
V2 (¿р) = —, (2)
т
где ¿р = рга+1 — р";
уп+1 — V
-= —V (¿р); (3)
т
дТ
+ V (У^Т) = V2T. (4)
Численная реализация (1)-(4) включает в себя следующие этапы.
1. В исследуемой области задается начальное распределение температуры, удовлетворяющее граничным условиям. Компоненты скорости полагаются нулевыми.
2. Из векторного уравнения (1) находится промежуточное поле скорости "V.
3. Рассчитываются поле разности давлений и давление по (2).
4. Из (3) получается окончательное поле скорости.
5. Путем решения (4) вычисляется поле температуры. Процесс повторяется до некоторого значения £ = пт.
В задаче использовалась разнесенная сетка, в которой давление, вязкость и температура определялись в центре элементарного объема, компоненты вектора скорости — в центрах плоскостей ячеек: и — в центре плоскости УХ, V — в центре плоскости XX, т — в центре плоскости XY. В декартовой системе координат равномерную сетку можно задать в виде
ж/ = (I — 1/2)Дж, Дж > 0, I = 0,1,...,М, (М — 1)Дж = X, у-1 = (3 — 1/2)Ду, Ду > 0, 3 = 0,1,...,Ж, (Ж — 1)Ду = У, (5)
гк = (К — 1/2)Дг, Дг > 0, К = 0,1,...,Ь, (Ь — 1)Дг = X,
где Дж, Ду, Дг — размеры шагов сетки; М, N Ь — число ячеек сетки в направлениях ж, у и г.
Для реализации (1) использовалась схема стабилизирующей поправки. Для вектора скорости на примере ж-компоненты и-схема выглядит так:
ига+1/3 — и"
- = ЬииП+1/3 + Ь22и" + Ьззи" + Ь21V™ + Ьз1^га — Ь1РП,
т
ип+2/3 — ип+1/3 , / Ч
и-и- = Ь22 (и"+2/3 — и") , (6)
ип+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,J,K — 2Т/ДК + -1,J,K)
= Т/ДК — Т^+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 (1IJ+1K - i,J,k + 1i,j—1,k)
рП+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 г.
