CC BY
39
Вычислительные технологии
Том 14, № 3, 2009
Моделирование трехмерной конвекции в мантии Земли с применением неявного метода слабой сжимаемости*
В. В. ЧЕРВОВ Учреждение Российской академии наук Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, Россия e-mail: elixirexpo@yandex.ru
Представлены результаты трехмерного моделирования конвекции в мантии Земли. Численная модель основана на применении неявного метода слабой сжимаемости.
Ключевые слова: тепловая конвекция в мантии Земли, метод искусственной сжимаемости, численное моделирование.
Введение
В настоящей работе продолжено тестирование задачи конвекции, предложенной в [1]. В работах [2-6] с использованием переменных "вектор завихренности — векторный потенциал метода дробных шагов и последовательности сеток" продемонстрировано хорошее согласие с результатами [1]. В [7] численное решение модельной задачи трехмерной теплогравитационной конвекции осуществлялось с использованием неявного метода расщепления по физическим процессам. Достаточно хорошо известным при решении задач гидродинамики несжимаемых жидкостей является метод слабой сжимаемости [8-11]. В настоящей работе с применением метода слабой сжимаемости получены результаты, не уступающие по точности результатам работ [1-7].
1. Математическая постановка задачи
Для описания течений в верхней мантии Земли привлекается хорошо известная математическая модель, включающая в себя обезразмеренные уравнения [12]:
* Работа выполнена при финансовой поддержке СО РАН (интеграционные проекты № 116 и № 44) и Российского фонда фундаментальных исследований (грант № 08-05-00276-а).
V- v = 0,
(1)
dT
— + v ■ VT = V2T.
(3)
© ИВТ СО РАН, 2009.
и = V = IV = 0
ди дн' ду ду
Рис. 1. Граничные условия для вектора скорости в параллелепипеде: на вертикальных гранях — условия проскальзывания; на горизонтальных плоскостях — условия прилипания
а рд.г //3 А У
Здесь р — давление, Т — температура, £ — время, Ба =--число Рэлея,
ПоХ
V = (и^,т) — вектор скорости, е = (0, 0,1), дх — ускорение силы тяжести, / — вертикальный размер конвектирующей области, АТ = Ттах — Тт;п, х — коэффициент температуропроводности, а — коэффициент теплового расширения, р, п — характерные плотность и динамическая вязкость.
Размерные значения (в системе СИ), которые были использованы в [1] и в настоящей работе, принимались следующими: / = 2.7 • 106 м, АТ = 3700 °С, х = 10-6 м2/с, а = 10-5 °С, р = 3300 кг/м3, дг = 10 м/с2.
Простейшей областью интегрирования является параллелепипед [0 : X], [0 : У] [0 : X]. Для вектора скорости на боковых гранях задаются условия проскальзывания, а на нижней и верхней гранях — условия прилипания (рис. 1):
— на поверхностях х = 0; х = X (0 < у < У, 0 < г < 1): и = д^/дж = дт/дх = 0;
— на поверхностях у = 0; у = У (0 < х < X, 0 < г < 1): ди/ду = V = дт/ду = 0;
— на поверхностях г = 0; г =1(0 < х < X, 0 < у < У): и = V = т = 0.
Для температуры, как и в [1], на боковых гранях ставятся условия теплоизоляции (адиабатическая стенка), т. е. первые производные по нормали на вертикальных стенках равны нулю. На верхней и нижней гранях ставятся условия Дирихле: нулевая температура на верхней и некоторая фиксированная температура (в безразмерных уравнениях равная единице) на нижней:
— на поверхностях х = 0; х = X (0 < у < У, 0 < г < 1): дТ/дж = 0;
— на поверхностях у = 0; у = У (0 < х < X, 0 < г < 1): дТ/ду = 0;
— на поверхности г = 0(0 < х < X, 0 < у < У): Т = 1;
— на поверхности г =1 (0 < х < X, 0 < у < У): Т = 0.
Система уравнений (1)-(3) устроена так [13], что в начальный момент времени £ = £0 задаются начальные условия лишь для температуры. В качестве начального распределения температуры, как ив [1], выбиралось следующее:
Т(х, у, г, £) = Т0(х, у, г) = (1 — г) + 0.2(cos(пx/X) + сов(пу/У))в1п(пг).
2. Численная модель трехмерных течений в естественных переменных, основанная на методе слабой сжимаемости
Как уже отмечалось выше, хорошо известным подходом к решению задач гидродинамики вязкой несжимаемой жидкости является метод искусственной сжимаемости. Для построения численной модели в случае естественных переменных применялся метод установления с использованием неявной схемы искусственной сжимаемости [8-11, 14] и метода дробных шагов [14]:
^ + с2 К+1) = 0, (4)
дт
дМга+1 + (др\га+1 = + дМга+1 + Ка. Тп ■ е-
дт) \ дх-) дхк V дхк дх-)
(5)
дТ
+ V (у^1 ■ Т) = У2Т. (6)
Параметр т является общим итерационным параметром для уравнений (4) и (5); с2 и т играют роль релаксационных параметров. Параметры с и т выбирались из условия устойчивости, предложенного в [15]: т < шт(Дх, Ду, Дг)/с.
Для задачи тестирования с переменной вязкостью результаты получены при с2 = 219.
Численная реализация (4)-(6) включает в себя следующие этапы.
1. В исследуемой области задается начальное распределение температуры, удовлетворяющее граничным условиям. Компоненты скорости полагаются нулевыми.
2. Из векторного уравнения (5) находится поле скорости уга+1; при этом рга+1 определяется из (4) и подставляется в (5).
3. После вычисления всех компонент скорости следует расчет давления из (4):
рга+1 = рп - тс2(Ь1 ип+1 + Ь2УП+1 + Ьз^га+1). (7)
4. Путем решения (6) вычисляется поле температуры.
Процесс повторяется до некоторого значения
гм = N ■ Дг.
Для реализации (5) использовалась схема стабилизирующей поправки. Для вектора скорости, на примере х-компоненты и, схема выглядит так:
—п+1/3 — —п
- = Ьц„П+1/3 + Ь22ПП + Ьззпп + ¿21 ^ + Ьз1ШП - Ь^, (8)
т
„п+2/3 - „п+1/3
---- = Ь22(„п+2/3 - „п), (9)
и
п+1 — —п+2/3
= ¿33(-п+1 - -п). (10)
т
В (7)—(10) использованы стандартные аппроксимации [16]:
д д д д д д д д д Ь1 = тг' Ь = тг' Ьз = тг' Ьп = ; = ; Ьзз = ;
дх ду дг дх дх ду ду дг дг
Ь д д ^ д д 21 дуП дх' 31 дг^ дх
Рис. 2. Двумерная разнесенная сетка
Рис. 3. Сеточный шаблон: в центре ячейки, кроме давления и температуры, вычисляется и вязкость П/ J к
Для V и т конечно-разностное представление аналогично (8)-(10). В задаче использовались как разнесенная, так и неразнесенная сетки. На разнесенной сетке давление, вязкость и температура определялись в центре элементарного объема, компоненты вектора скорости — в центрах плоскостей ячеек: и — в центре плоскости у0г, V — в центре плоскости х0г, т — в центре плоскости х0у. На рис. 2 для простоты и наглядности изображена двумерная разнесенная сетка, сеточный трехмерный шаблон представлен на рис. 3.
3. Результаты расчетов
Тестирование численной модели осуществлялось решением модельной задачи [1]. Решение для переменной вязкости отыскивалось в единичном кубе (X = У = Z =1).
При этом задавались следующие параметры: масштабный множитель при вязкости
По = 1.20165 ■ 1024; ) = exp[9/(T + 0) - 0/(0.5 + 0)]; 9 = 225/ln(r) - 0.25ln(r); 0 = 15/ln(r) - 0.5; r = n |T=0 /n |T=i= 20; Ra = 20 000.
Для решения задачи вводилась равномерная в каждом направлении сетка. Вычислялись следующие параметры:
(1) среднеквадратичная скорость:
где Stop — верхняя поверхность параллелепипеда;
(iii) значение вертикальной компоненты скорости w и температуры T в угловых точках среднего сечения конвективного слоя;
(iv) средняя температура:
вычисляемая на горизонтальных сечениях области Sz=0.75 и Sz=0.50, на глубинах z = 0.75 и z = 0.50;
(v) значения zi — отклонения свободной поверхности Земли от геоида в угловых точках верхней поверхности куба.
Интегралы вычислялись с применением квадратурной формулы трапеций.
Результаты расчетов автора (Che) сопоставлены с данными У. Кристенсена (Chr) как наиболее полными из имеющихся в статье [1] и представлены в табл. 1 и 2. Согласие достаточно хорошее. Изложенный выше алгоритм был с успехом применен к решению модельной задачи о протекании мантийной жидкости [7].
где А — объем параллелепипеда со сторонами X,Y и Z; (ii) число Нуссельта (Nu) по формуле [17]:
Таблица 1. Переменная вязкость. Разнесенная сетка. Естественные переменные. Искусственная сжимаемость
№ Параметр Сетка 32 х 32 х 64 Егг, %
п/п СЬг СЬе
1 3.03927 2.99409 1.51
2 Vттв 35.132 35.159 0.08
3 Ц0,0,1/2) 165.91 165.967 0.04
4 w(0,Y, 1/2) -26.72 -26.977 0.95
5 1/2) -58.23 -58.755 0.89
6 Т(0, 0,1/2) 0.90529 0.9052 0.01
7 Т (0,Г, 1/2) 0.49565 0.4955 0.02
8 Т(X, Г, 1/2) 0.23925 0.2388 0.20
9 Тт (3/4) 0.56593 0.5654 0.09
10 Тт (1/2) 0.58158 0.5890 1.26
11 21(0, 0,1/2) 10869.0 10389.0 4.62
12 *1(0,Г, 1/2) -4145.00 -4271.0 2.95
13 21 (X, Г, 1/2) -12811.0 -12761.6 0.39
Таблица 2. Переменная вязкость. Неразнесенная сетка. Естественные переменные. Искусственная сжимаемость
№ Параметр Сетка 32 х 32 х 64 Егг, %
п/п СЬг СЬе
1 3.03927 3.04469 0.18
2 ^ттв 35.132 35.0888 0.12
3 w(0, 0,1/2) 165.91 165.170 0.44
4 w(0,Y, 1/2) -26.72 -26.601 0.45
5 w(X, Г, 1/2) -58.23 -58.523 0.50
6 Т(0, 0,1/2) 0.90529 0.9062 0.11
7 Т(0, Г, 1/2) 0.49565 0.4996 0.80
8 Т(X, Г, 1/2) 0.23925 0.2409 0.70
9 Тт(3/4) 0.56593 0.5649 0.17
10 Тт(1/2) 0.58158 0.5831 0.26
11 21(0,0,1/2) 10869.0 11070.8 1.82
12 21(0,У,1/2) -4145.00 -4230.92 2.03
13 21(Х,У,1/2) -12811.0 -13035.0 1.72
Относительная ошибка вычислялась по формуле:
Егг
СЬе - СЬг СЬг
■ 100 %.
Заключение
Основные результаты работы сводятся к следующему. Построена трехмерная численная модель конвекции в мантии Земли, основанная на методе слабой сжимаемости. Продемонстрирована ее достаточно высокая эффективность.
Автор выражает благодарность профессору, д-ру физ.-мат. наук Геннадию Георгиевичу Черных за помощь при постановке задачи и постоянное внимание к работе.
Список литературы
[1] Busse F.H., Christensen U., Clever R. et al. 3D Convection at infinite Prandtl number in cartesian geometry — a benchmark comparison // Geophiys. Astrophys. Fluid Dynamics. 1993. Vol. 75. P. 39-59.
[2] Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычисл. технологии. 2002. Т. 7, № 1. С. 114-125.
[3] Червов В.В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением последовательности сеток // Вычисл. технологии. 2002. Т. 7, № 3. C. 85-92.
[4] Тычков С.А., Червов В.В., Черных Г.Г. О численном моделировании тепловой конвекции в мантии Земли // Докл. РАН. 2005. Т. 402, № 2. С. 248-254.
[5] Tyohkov S.A., Chervov V.V., Chernykh G.G. Numerical modeling of 3D convection in the Earth mantle // Russ. J. Numer. Anal. Modell. 2005. Vol. 20, N. 5. Р. 483-500.
[6] Тычков С.А., Червов В.В., Черных Г.Г. Численная модель трехмерной конвекции в верхней мантии Земли // Физика Земли. 2005. № 5. С. 48-64.
[7] Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода расщепления по физическим процессам // Вычисл. технологии. 2006. Т. 11, № 4. C. 73-86.
[8] Владимирова Н.Н., Кузнецов Б.Г., Яненко Н.Н. Численный расчет симметричного обтекания пластинки плоским потоком несжимаемой жидкости // Некоторые вопросы прикладной и вычислительной математики. Новосибирск,1966. C. 186-192.
[9] Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л.: Гидрометеоиздат, 1986.
[10] Флетчер К. Вычислительные методы в динамике жидкостей: М.: Мир,1991. Т. 1, 2.
[11] Черный С.Г., Чирков Д.В., Лапин В.Н. и др. Численное моделирование течений в турбомашинах. Новосибирск: Наука, 2006. 202 с.
[12] Доврецов Н.Л., Кирдяшкин А.Г. Глубинная геодинамика. Новосибирск: НИЦ ОИГГМ СО РАН, 1994.
[13] Федорюк М.В. Характеристики течений несжимаемой жидкости в гравитационном поле / Матем. c6. 1988. Т. 137(179), № 4(12). С. 483-499.
[14] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ние, 1967.
[15] Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во ИГУ, 1990.
[16] Самарский А.А. Теория разностных схем. М.: Наука, 1977.
[17] Blankenbach В., Busse F. et al. A benchmark comparison for mantle convection codes // Geophys J. 1989. Int. 98. P. 23-38.
Поступила в редакцию 19 ноября 2008 г., в переработанном виде — 24 декабря 2008 г.
