Численная модель трехмерной конвекции в мантии Земли Текст научной статьи по специальности «Физика»

Вычислительные технологии Том 11, часть 2, Специальный выпуск, 2006

ЧИСЛЕННАЯ МОДЕЛЬ ТРЕХМЕРНОЙ КОНВЕКЦИИ В МАНТИИ ЗЕМЛИ*

С. А. Тычков, В. В. Червов Институт геологии и минералогии СО РАН, Новосибирск, Россия

Г. Г. Черных

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: chernykh@ict.nsc.ru

A 3d numerical model of convection processes in Earth mantle is constructed on the basis of vorticity — vector potentional variables and the method of fractional steps. The detailed testing of the model was carried out.

Введение

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

При решении трехмерных задач гидродинамики достаточно хорошо зарекомендовали себя переменные завихренность — векторный потенциал [8, 9]. Анализ известных работ показывает, что при численном моделировании конвективных процессов в мантии Земли указанным переменным уделено недостаточное внимание. Так, например, в статье [10] выполнено решение только для постоянной вязкости.

В настоящей работе предпринята попытка построения численной модели конвективных процессов в мантии Земли, основанной на вышеупомянутых переменных и методе дробных шагов [11]. Осуществлено сопоставление с результатами работы [1]. Работа является продолжением и развитием [12-14].

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

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

ди дь дь) дх ду дг

*Работа выполнена при финансовой поддержке СО РАН (интеграционный проект № 116). © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

др ^ д ди ^ д {ди ^ ду \ ^ д ^ди ^ дт\ ^

дх дх дх ду \ду дх) дг \дг дх) '

др д (ди ду \ ^ д ду д / ду дт \ , ,

ду дх ^ у ду дх) ду^ ду дг ^ у дг ду ) '

др д / ди дт \ д / ду дт \ ^ д дт ^ дг дх ^ у дг дх ) ду ^ у дг ду ) дг ^ дг

дт дт дт дт ^2гт1 „

— + и— + у— + т— = У2Т + д. 5

дЬ дх ду дг

Здесь и,У,Ю — компоненты вектора скорости; р — давление; Т — температура; — источник тепловыделения; Яа — число Рэлея; п — динамическая вязкость.

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

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

ш = 1ШХ + ]иу + ки* = V х V, ф = \фх + jфy + кф*. (6)

В результате система уравнений (1) (1) переходит в следующую:

V2 ф = -ш; (7)

V2 [пих] = ПххШх + ( Пхихх + Пуиух + пуиХ) + + яаТу; (8)

V2[пиу] = Пууиу + ( ПхиУх + ПуиУу + пхи'у) + ^2 - яаТх; (9)

V2[пиу] = Пууиу + ( Пхих + ПуиУу + пуиУу) + Р3, (10)

где

= 2(пхх Юу - Пуу у у) + 2Пуу (Уу - Юу) - Пху (иу + -Юх) + Пхг (Ух - щ), 1

^2 = 2(ПххПу - Пуу Юх) + 2Цху (Юу - их) - Пуу (Ух + Пу) + Пху(Юу - Уу), > (11)

Fз = 2(Пуу Ух - ПххЩ) + 2 Пху (их - Уу) - Пху (Юу + Уу) + Пуг(иу - Юх). )

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

Граничные условия проскальзывания для функций ш, ф: на поверхностях х = 0 х = X (0 < у < У, 0 < г < 1):

дфх дих дУ

^Г=ФУ = Ф* = 0, ^=иу = 0, = дх дх дх

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

д-^ = Фх = Ф* = о, = ^ =

ду ду ду

Условия прилипания для функций ш, ф: на поверхностях х = 0 х = X (0 < у < У, 0 < г < 1):

дфу дфх дт Зу

ф* = -£-=ф» = -£-=ф'=0, сох = 0, ^ = - —, и/ = —;

дх дх дх дх

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

дГ =Г = Г=д-^=Г = 0, = ^ = и/ ди•

ду ду ■> х ду-> 1 ду1

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

д41=Г = д41=Фу = Г = 0, 0,* = -^, = и/ = 0.

дг дг дг дг

2. Алгоритм расчета

Конечно-разностный алгоритм решения задачи основан на применении метода дробных шагов [11]. Уравнения (5), (7)-(Ю) интегрируются с помощью схемы стабилизирующей поправки, являющейся итерационной для уравнений (7)—(10). Алгоритм расчета включает на каждом временном слое (до выхода на стационарный режим) следующие этапы:

1. При известном распределении температуры вычисляются завихренность, векторный потенциал и компоненты скорости (до сходимости итерационных процессов).

2, Вычисляется поле температуры.

На примере уравнения переноса тепла (5) схема стабилизирующей поправки выглядит следующим образом:

Тп+1/3 _ тг'

+ иЬхТп+1^ - ЬххТп+1'3 = —уЬуТп + ЬууТп - %юЬхТп + ЬххТп

п

Т'п+2/3 _ тп+1/3

+ уЬуТп+2/3 _ ЬууТп+2/3 = уЬуТп _ ЬууТ

П

Тп+1 _ Тп+2/3

--Ь тЬхТп+1 - ЬххТп+1 = шЬхТп - ЬХХТ

п

(12)

где ^ — величина шага по времени, п — помер временного слоя.

В схеме (12) трехточечные разностные операторы Ьхх, Ьуу, Ьххъ Ьх,Ьу,Ьх аппроксимируют дифференциальные следующим стандартным образом (оставлен лишь один индекс):

д2/

хх

дх2

ду2 д2/

Ьуу /

/г+1 -2/ + /.-1 кх д£ х = Ьх/ = /г+1 — /г—1

- 2/ + Л-1 = Ьу / = Л+1 ~~ Л -1

ч у 2НУ

Л+1 -2/ + Д-1 д£ = Ьг / = Л+1 ~ Л-1

г

г2

Для реализации каждого дробного шага схемы (12) применялись скалярные трехточечные прогонки.

Граничные условия для завихренности реализовывались с применением формул Тома с дополнительным усреднением [8, 17].

3. Тестирование на модельной трехмерной задаче конвекции в мантии Земли

Тестирование численной модели осуществлялось путем решения модельной задачи [1] (аналогичные результаты получены и при решении более простой задачи с постоянной вязкостью). Решение отыскивалось в единичном кубе. При этом задавались следующие параметры: масштабный множитель при вязкости

По = 1.20165 ■ 1024, п(Т) = ехр [в/ (Т + О)] - [в/ (0.5 + О)],

причем

в = [225/ 1п(г)] - 0.25 1п (г), О = 15/ 1п (г) - 0.5,

г=ч\т=а/ч\т^=Ж, Ка = адгр£АТ = ^ ^ ^

ПоХ

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

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

\

(и2 + v2 + w2)dxdydz /> ,

XYZ

A

где A — объем параллелепипеда со сторонами X, Y и Z; — число Нуссельта (Nu) по формуле [3]

ЯдТ —dxdy

Stop

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

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

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

Y

[ дТ

— интегральный параметр, вычисляемый по формуле г (х, ;) / — с/у вдоль линии,

J dz о

Y

(XZ )-плоскости;

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

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

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

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

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

В качестве начального распределения температуры выбиралось распределение

Т(х, у, г, ¿0) = Т0(х, у, г) = (1 — г) + 0.2(еов(пх/Х) + еов(пу/У)) вт(пг).

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

СТС — СЬг

Err

Chr

■ 100 %.

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

Хорошо известным эффективным подходом к решению задач математической физики является метод последовательности сеток [18]. Результаты его применения в настоящей работе иллюстрируются табл. 2.

При непосредственных вычислениях на сетке 48x48x48 потребовалось 660 мин 24 с на персональном компьютере с процессором Athlon 1000. Таким образом, выигрыш во времени счета на последовательности сеток более чем восьмикратный.

Таблица 1. Результаты расчетов авторов и Кристенсена [1]

Наименование Результаты Крис- Результаты Относительная

параметра тенсена на сетке авторов на сетке ошибка (Err), %

32x32x64 32x32x32

Nu 3.0393 3.0400 0.0240

Vrms 35.132 35.180 0.1366

w( 1,1,0.5) -58.230 -58.230 0.0000

Т( 1,1,0.5) 0.2393 0.2379 0.5643

0.7684 0.7731 0.6117

т(1,1) -0.1388 -0.1356 2.3055

тт(ощ 0.5816 0.5799 0.2889

U)z(0.75,0.25,0.75) -11.125 -11.25 1.1236

Таблица 2. Расчеты авторов на последовательности сеток

Наименование Результаты Крис-

параметра тенсена на сетке Результаты авторов на сетках

32x32x64

12x12x12 24x24x24 48x48x48

Nu 3.0393 3.2450 3.0397 3.040

Vrms 35.132 37.970 35.360 35.07

W( 1,1,0.5) -58.230 -60.68 -58.262 -58.42

Г(1,1,0.5) 0.2393 0.2326 0.2377 0.239

#(1,1) 0.7684 0.8017 0.7804 0.773

т(1, 0.25) -0.1388 -0.0676 -0.1341 -0.140

Tm(0.50) 0.5816 0.5844 0,5792 0.582

ujz(0.75,0.25,0.75) -11.125 -11.030 -11.203 -11.35

Время счета 65 с 119с 4811с

Существенный выигрыш во времени счета может быть достигнут также с применением экстраполяции по Ричардсону [19]. В частности, на рассмотренной модельной задаче [1] дополнительный выигрыш (для сеток 12x12x12, 24x24x24, 48x48x48) составил «24 раза. При этом рассчитанные на первых двух сетках решения комбинировались и сопоставлялись с расчетом на третьей сетке:

Ф

т+1,т ,3,к

(4 ■ Фт+и-^-: - ф^Ь) /3;

значения индексов г,], к изменяются от 1 до своих максимальных значений, а индекс т = 1, 2, 3 означает номер сетки в последовательности; Ф — одна из искомых функций.

Приведенные выше результаты численных экспериментов свидетельствуют о достаточно высокой эффективности и конкурентоспособности построенной численной модели. Отметим, что в [17] при решении двумерных гидродинамических задач конвекции на последовательности сеток получен выигрыш в 3-5 раз.

В процессе расчетов дополнительно осуществлялся контроль закона сохранения тепла. Следствием исходных дифференциальных уравнений, начального и граничных условий является закон сохранения тепла

1///™ =

А (•*=!)

^ЬоМот (• — О)

Для его проверки осуществлялась последовательность вычислений.

п

1. На каждом шаге по времени Ьп = тк вычисляются интегралы

к—О

1п = ^Т(1А, 1пМр =Л {^^(Ьйу, /га>Ьойот = УУ |

А 5гор(г—1) Яо-иот^—О)

2. Сравнивается величина 1п с величиной

'дТ

дг

dxdy.

10 (^к^ор — 1к,Ьойош )тк, Еп = !п —

к—О

10 + Е (1к ,top 1k,bottom)тk к—О

Результаты контроля закона сохранения тепла, полученные на сетке 32x32x32, представлены в табл. 3. Для сравнения в последнем столбце помещены результаты расчетов на сетке 48x48x48.

п

п

Таблица 3. Анализ закона сохранения тепла

млн лет \Еп\ \Еп/10\ | {1п,Ьор Лг,Ьойот) /Е \ на сетке 32x32x32 1 (1п,top Е,bottom ) /10| на сетке 48x48x48

51.7 1.17Т0-5 2.13Т0"5 2.20Т0"5 3.66Т0"6'

460.9 1.56Т0"4 2.83Т0"4 8.80Т0"6 1.42Т0-6

918.5 2.99Т0"4 5.44Т0"4 7.32Т0"6 1.20Т0-7

1364.5 4.32Т0"4 7.86Т0"4 7.01Т0"6 1.18Т0-7

1711.2 5.35Т0"4 9.73Т0"4 6.83Т0"6 1.13Т0"7

1991.6 6.17Т0-4 1.12Т0-3 6.65Т0"6 1.10Т0"7

2104.1 6.51Т0"4 1.18Т0"3 6.56Т0"6 1.05Т0-7

Выход на стационарное решение в интегральном смысле характеризуется поведением (/n,top — /n,bottom) //о. Интегралы вычислялись по-прежнему с помощью квадратурных формул трапеций,

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

Заключение

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

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

[1] Busse F.H., Christensen и., Clever R. ет al. 3D Convection at infinite Prandl number in cartesian geometry — a benchmark comparison // Geophys. Astrophys. Fluid Dynamics. 1993. Vol. 75. P. 39-59.

[2] Gable W., O'Connell J. Convection in three dimensions with surface plates: generation of toroidal flow //J. Geophys. Research. 1991. Vol. 96, N B5. P. 8391-8405.

[3] Glatzmaier G.A., Schubert G. Three-dimensional spherical model of layered and whole mantle convection //J. Geophys. Research. 1993. Vol. 98, N B12. P. 21969-21976.

[4] Machetel P., Thoravan C., Brunet D. Spectral and geophysical of 3-D spherical mantle convection with an endothermic phase change at the 670 km discontinuity // Phys. Earth and Planetary Interiors. 1995. Vol. 88. P. 43-51.

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

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

[7] Zhong S., Zuber М. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection //J. Geophys. Research. 2000. Vol. 105, N B5. P. 11063-11082.

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

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

[10] Sparks D.W., Parmentier Е.М., Morgan J.P. Three-dimensional mantle convection beneath spreading centers implications for along-axis variations in crustal thickness and gravity // J. Geophys. Research. 1993. Vol. 98, N B12. P. 21977-21995.

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

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

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

[14] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в верхней мантии Земли // 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.

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

[16] Федорюк М.В. Характеристики течений несжимаемой жидкости в гравитационном поле // Мат. сб. 1988. Т. 137(179), № 4(12). С. 483-499.

[17] Тарунин E.JI. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во Иркут. ун-та, 1990.

[18] Федоренко Р.П. О скорости сходимости одного итерационного процесса // Журн. вычисл. математики и мат. физики. 1964. Т. 4, № 3. С. 559-564.

[19] Марчук Г.И., Шайдуров В.В. Повышение точности разностных схем. М.: Наука, 1979.

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

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

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