CC BY
32
Вычислительные технологии
Том 7, № 3, 2002
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ ЗАДАЧ КОНВЕКЦИИ В МАНТИИ ЗЕМЛИ С ПРИМЕНЕНИЕМ ПОСЛЕДОВАТЕЛЬНОСТИ СЕТОК *
В. В. Червов
Объединенный институт геологии, геофизики и минералогии СО РАН
Новосибирск, Россия e-mail: [email protected]
The results of numerical simulation of model convection problem in Earth mantle based on using of the series of grids are presented. The results show the essential gain in CPU time.
Введение
Результаты численного решения задач математической физики с применением последовательности сеток приведены в [1-4]. В работах Е. Л. Тарунина [2, 3] показано, что для двумерных задач конвекции объем вычислений методом установления на последовательности сеток дает значительную экономию (а = 3 ^ 5) затрат машинного времени. Цель настоящей работы — демонстрация эффективности такого подхода при решении трехмерных задач тепловой конвекции в мантии Земли (на модельной задаче).
1. Задача конвективного тепломассопереноса в мантии Земли
Постановка задачи подробно приведена в [5], ниже излагается лишь ее краткое описание для рассматриваемой модельной ситуации с использованием переменных ф, ш: ф = \фх + jфy + кф*, ш = 1шх + jшy + кш*.
В области П (единичный куб) (0 < х < 1, 0 < у < 1, 0 < г < 1) отыскивается решение системы дифференциальных уравнений:
У2фх = -шх, У2фу = -шу, У2ф* = -ш*, (1)
У2Сх = Пххшх + (ПхшХ + Пу шХ + П* шХ) + + Я&Ту, У2Су = Пуушу + (пхшх + Пушу + П*ш*) + ¿2 - ИаТх,
* Работа выполнена при поддержке Интеграционных проектов СО РАН №27 и 30. © В. В. Червов, 2002.
= 'Цгг Ч + (Пх<Х + Пу <1 + Пг <) + ¿3,
дт дт дт дт 2 — + и— + V— + ш— = V т.
дЬ дх ду дг
(2) (3)
На верхней и нижней гранях куба ставились условия прилипания для скорости течения жидкости:
и = V = т = 0, г = 1, 0 < х < 1, 0 < у < 1, г = 0, 0 < х < 1, 0 < у < 1
или, в переменных ф, и:
дфх фХ дфу
фу = фг
0,
их = -
дv дг
и
ди дг
и
0.
дг дг
Температура на верхней грани приравнивалась нулю:
т = 0, г =1, 0 < х < 1, 0 < у < 1. На подошве куба температура равнялась единице:
т = 1, г = 0, 0 < х < 1, 0 < у < 1.
На боковых границах выделенного объема ставились симметричные условия как для температуры, так и для вектора скорости:
дт дv дт
дх дх дх
0,
х = 0, 0 < у < 1, 0 < г < 1, х =1, 0 < у < 1, 0 < г < 1,
или в переменных ф, и:
дфх
дх
фу = фг
дт ди дт 0
ду ду ду
или в переменных ф, и:
дфу ду
фх = фг = 0,
дих
дх у
у
диу
и
0,
и
дv дх
0, 0 < х < 1,0 < г < 1,
1, 0 < х < 1,0 < г < 1,
и
0,
ди
и = --
ду ду
Техника вычислений граничных условий в переменных ф, и описана, например, в работах [2, 6-8].
При Ь = 0 задаются лишь начальные условия для температуры [9]: т(х,у,г, 0) = т(х, у, г). В качестве т в настоящей работе, как и в [5, 10], выбиралась функция т(х, у, г) = (1 - г) + 0.2(сов(пх/Х) + соъ(пу/У)) вт(пг).
В уравнениях (1)-(3), граничных и начальных условиях приняты обозначения: и — вектор скорости с компонентами и = и(х,у,г,Ь), V = v(x,y, г,Ь), т = т(х,у, г,Ь); И = го^; ш = гс^И; т(х,у,г,Ь) — температура жидкости; ц(х,у, г,Ь) > 0 — коэффициент динамической вязкости;
V
д д д 1— + ^ + к
дх ду
дг
С =
= 2(Пгг ту - Пуу Vz) + 2^ ^у - ) - Пху (иг + Шх) + Пхг ^х - Щ); = 2(Пххиг - Пгг Шх) + 2Цхг (Ш - их) - Пуг ^х + Щ) + Пху (ту - Vх);
0
¿3 = 2(ПууУх - ПххПу) + 2Пху(их - %) - Пх*(^у + V*) + Пу*(и* - ^х).
Определяющим параметром задач тепломассопереноса в мантии Земли является число Рэлея:
И,а = ад* рй3ДТ/п0Х,
где й — расстояние от поверхности Земли до ядромантийной границы; ДТ — разность температуры на подошве мантии и температуры на поверхности; п0 — масштабный множитель для динамической вязкости; д* — г-компонента вектора силы тяжести g = (0,0, -д*), т.е. ускорение свободного падения; а — коэффициент теплового расширения жидкости.
2. Результаты расчетов
Алгоритм решения и его численная реализация подробно изложены в [5]. Алгоритм сводится к последовательному интегрированию уравнений Пуассона (1), (2) методом установления с применением схемы стабилизирующей поправки [11]. Вычисления температуры завершают расчеты на текущем шаге по времени. Как ив [5], в настоящей статье решалась тестовая модельная задача Busse et al. [10]. Вычислялись следующие параметры: — среднеквадратичная скорость
Vrm
\
^ XYZ J J J (u2 + v2 + w2)dxdydz
где A — объем параллелепипеда со сторонами X = 1, Y = 1 и Z = 1; — число Нуссельта по формуле [Blankenbach, 12]
Nu = -(XY) / [ Tzdxdy,
Stop
dT
где Tz = ——, Stop — верхняя поверхность куба; dz
— значение вертикальной компоненты скорости w и температуры T в угловых точках среднего сечения конвективного слоя;
dT
— значение теплового потока $ = —— в угловых точках верхней поверхности куба;
dz
— интегральный параметр, вычисляемый по формуле
т (x,z) = / jzdy
о
вдоль линии, параллельной оси Y, начинающейся точками (0, 0.25), (1/2, 0.25), (1, 0.25) фронтальной (XZ)-плоскости;
— средняя температура Tm = f f Tdxdy, вычисляемая на горизонтальных сечениях
области Sz=0.75 и Sz=0.50, на глубинах z = 0.75 и z = 0.50;
— значение вертикальной компоненты uiz вектора завихренности в точке (0.75, 0.25, 0.75), которая выражается через компоненты вектора скорости следующим образом: uiz = dv du
dx dy
В работе Busse [10] наиболее полно представлены результаты Кристенсена (Chr), которые можно принять за эталонные. Сопоставление результатов расчетов Кристенсена и автора (Che) представлено в табл. 1-4. Расчеты выполнялись на ПЭВМ Athlon 1000 MGh с двойной точностью.
Таблица 1
Параметр Результаты Кристенсена Результаты автора на сетках
на сетке 32 х 32 х 64 12 х 12 х 12 24 х 24 х 24 48 х 48 х 48
Nu 3.03927 3.245 3.0397 3.040
^rms 35.132 37.97 35.360 35.07
w(1,1,0.5) -58.23 -60.68 -58.262 -58.42
T(1,1,0.5) 0.23925 0.2326 0.2377 0.2393
0 (1,1) 0.7684 0.80172 0.7804 0.7726
т(1, 0.25) -0.1388 -0.06761 -0.1341 -0.140
Tm(0.75) 0.56593 0.5764 0.5638 0.5644
Tm(0.50) 0.58158 0.5844 0.5792 0.5815
Q(0.75,0.25,0.75) -11.125 -11.03 -11.203 -11.35
Время счета 65 с 119 с 4811 с
При непосредственных вычислениях на сетке 48 х 48 х 48 потребовалось 660 мин 24 с. Таким образом, выигрыш во времени счета на последовательности сеток более чем восьмикратный.
В табл. 2 сведены для сравнения относительные ошибки, которые вычислялась по формуле
Err
Che - Chr
Chr
100%.
Анализируя результаты, приведенные в табл. 2, можно сделать вывод о том, что дальнейшее измельчение сетки не приведет к существенному улучшению решения, так как полученное решение на сетке 48 х 48 х 48 является "точным". Более того, уже на сетке 24 х 24 х 24 относительная ошибка практически по всем параметрам не превышает двух процентов. Исключение составляют только ш(0,0,0.5) и интегральный параметр т(1,0.25), относительная ошибка для которых, впрочем, составляет вполне допустимую величину в пределах 4 %.
Таблица 2
Параметр Относительная ошибка, %
на сетке 12 х 12 х 12 на сетке 24 х 24 х 24 на сетке 48 х 48 х 48
Nu 6.347 0.015 0.028
^rms 7.476 0.650 0.180
w(1,1,0.5) 4.034 0.054 0.326
T(1,1,0.5) 2.847 0.670 0.013
0(1,1) 4.1555 1.540 0.547
т (1,0.25) 105.3 3.500 1.038
Tm(0.75) 1.814 0.364 0.265
Tm(0.50) 0.4785 0.412 0.017
Q(0.75,0.25,0.75) 0.8848 0.699 2.02
Следует заметить, что ухудшение некоторых результатов при переходе от сетки 24 х 24 х 24 к сетке 48 х 48 х 48 имеет, по-видимому, случайный характер, если ошибка менее одного процента. В случае, если ошибка превышает 2 %, могут потребоваться дополнительные исследования на других последовательностях сеток. Например, на последовательности сеток 16 х 16 х 16 - 32 х 32 х 32 - 64 х 64 х 64 были получены результаты, представленные в табл. 3. В табл. 4 сведены для сравнения результаты, полученные на двух последовательностях.
Таблица 3
Параметр Результаты Кристенсена Результаты автора на сетках
на сетке 32 х 32 х 64 16 х 16 х 16 32 х 32 х 32 64 х 64 х 64
3.03927 3.084 3.037 3.041
^гтБ 35.132 36.44 35.23 35.12
Ц1,1,0.5) -58.23 -58.79 -58.29 -58.45
Т(1,1,0.5) 0.23925 0.2342 0.2386 0.2391
0(1,1) 0.7684 0.7855 0.77576 0.77023
т (1,0.25) -0.1388 -0.1155 -0.1414 -0.1377
Тт(0.75) 0.56593 0.5661 0.5636 0.5648
Тт(0.50) 0.58158 0.5777 0.5802 0.5820
0(0.75,0.25,0.75) -11.125 -11.19 -11.40 -11.30
Время счета 2 мин 16 мин 1064 мин
Прямое вычисление на сетке 64 х 64 х 64 займет, по приближенным оценкам, 6 суток, т. е. 144 ч. Выигрыш составляет порядка восьми раз.
Таблица 4
Параметр Относительная ошибка на сетках, %
12 х 12 х 12 16 х 16 х 16 24 х 24 х 24 32 х 32 х 32 48 х 48 х 48 64 х 64 х 64
6.3470 1.4650 0.0150 0.0775 0.0280 0.1047
^гтБ 7.4760 3.5960 0.6500 0.2895 0.1800 0.0653
Ц1,1,0.5) 4.0340 0.9488 0.0540 0.1870 0.3260 0.3740
Т(1,1,0.5) 2.8470 2.1510 0.6700 0.2561 0.0130 0.0590
0(1,1) 4.1555 2.1742 1.5400 0.9485 0.5470 0.2370
т (1,0.25) 105.30 20.140 3.5000 1.8500 1.0380 0.8011
Тт(0.75) 1.8140 0.0290 0.3640 0.4163 0.2650 0.2041
Тт(0.50) 0.4785 0.6728 0.4120 0.2439 0.0170 0.0698
0(0.75,0.25,0.75) 0.8848 0.5622 0.6990 2.4040 2.0225 1.5530
На рис. 1-4 изображены графики изменения параметров решения в зависимости от номера итерации (п), которые получены при решении задачи на сетке 64 х 64 х 64. Начальное распределение температуры было пересчитано с поля температур, которое вычислено на сетке 32 х 32 х 32. Безусловная сходимость к предписанным значениям очевидна. Сравнение результатов, полученных на сетках 48 х 48 х 48 и 64 х 64 х 64, не выявило заметного уменьшения относительной ошибки. Вполне возможно, что и на более мелких сетках результаты окажутся не лучше полученных на сетке 48 х 48 х 48. Главным, впрочем, является то, что для данного класса задач вполне оправдан метод последовательных сеток для нахождения решения.
Рис. 1. Рис. 2.
Рис. 3. Рис. 4.
3. Модельная задача. Уравнение Пуассона
Задача численного интегрирования уравнений (1)-(3) достаточно сложна, поэтому был рассмотрен существенно более простой пример:
V2 Ф(х, у, г) = ^ (х, у, г) = в [6 + 4в (х2 + у2 + г 2)]ев(х2+у2+г2) (4)
в кубе со стороной, равной единице. На верхней грани при г = 1, 0 < х < 1, 0 < у < 1 ставилось условие Неймана:
д Ф
— (х,у, 1) = 2вФ (х,у, 1). (5)
На нижней грани при г = 0, 0 < х < 1, 0 < у < 1 ставилось условие Дирихле
Ф(х,у, 0) = ев (х2+у2). (6)
На боковых границах выделенного объема также ставились условия Неймана
дф(0,у,г) = 0, 0 < у < 1,0 < г < 1;
дг
дф (1, у, г) = 2вФ(1,у,г), 0 < у < 1,0 < г < 1; дг
^(х,0,г) = 0, 0 < х < 1, 0 < г < 1;
дг
дф (х, 1, г) = 2вФ (х, 1, г), 0 < х < 1, 0 < г < 1. дг
(7)
Эта задача имеет точное аналитическое решение Фа = Фа(х, у, г) = ев(х2+у2+г2). Численное интегрирование осуществлялось, как и уравнений (1), (2), с применением итерационной схемы стабилизирующей поправки и последовательности сеток. Параметр в варьировался от -0.2 до -12.0. На последовательности сеток {16 х 16 х 16- 32 х 32 х 32-64x64 х 64} для в = -2.0 время счета 2.03 с + 19.28 с +277.38 с = 298.69 с. На сетке 64 х 64 х 64 время счета 541.40 с.
Ошибка на аналитическом решении вычислялась по формуле (тах Фа = 1)
5т = тах \Ф8^к - (Фа)^| . Критерий выхода из итерационного процесса:
(9)
em = max |Ф?,, - ФГЛI < 5 ■ 10-8.
i,j,k
i,j,k |
Здесь и в (9) индекс т = 1, 2, 3 обозначает номер сетки в последовательности; в итерации; Фа — аналитическое решение.
Результаты решения модельной задачи (4) - (8) приведены в табл. 5 и 6.
номер
Таблица 5
Таблица 6
Сетка 1 (16 х 16 х 16) Сетка 2 (32 х 32 х 32) Сетка 3 Сетка 1 (64 х 64 х 64) (16 х 16 х 16) Сетка 2 (32 х 32 х 32) Сетка 3 (64 х 64 х 64)
ti = 2.03 с ¿2 = 19.28 с ¿3 = 277.38 с ti = 2.30 с t2 = 24.01 с ¿3 = 383.82 с
й1 = 3.44 ■ 10-3 й2 = 7.78 ■ 10-4 й3 = 1.85 ■ 10-4 й1 = 3.44 ■ 10-3 й2 = 7.78 ■ 10-4 й3 = 1.83 ■ 10-4
— й1 — = 44 й2 й2 V =4.2 — й1 — = 44 й2 й2 = 4.2
Выигрыш составляет 1.8 раз. В работе использовались интерполяция сплайнами и линейная интерполяция. Как и в работах [2, 3] для двумерных задач конвекции, так и в результате исследований данной работы было показано, что лучшие результаты дает все же вариант линейной интерполяции.
При повышении на порядок критерия выхода из итерационного процесса, т. е. при
£m = max |Ф j - Ф?"к1 < 5 ■ 10-9,
i,j,k
' i,j,k I
выигрыш составляет 1.58 раз. Время счета на сетке 64 х 64 х 64 равно 647.63 с. Времена на последовательности сеток приведены в табл. 6.
Итак, показано, что при повышении точности выигрыш во времени для модельной задачи не превышает значения 1.6 раз, что вполне согласуется с данными, полученными в результате исследований [1]. Уменьшение ошибки на последовательности сеток также вполне согласуется с порядком аппроксимации итерационной схемы стабилизирующей поправки.
Поскольку при решении задач на основе рассмотренного выше подхода получается последовательность сеточных решений, возникает вопрос о повышении точности с применением экстраполяции Ричардсона [13]. Выполненный автором предварительный анализ показал эффективность этого подхода. Более детально результаты использования экстраполяции Ричардсона автор предполагает изложить в следующей работе.
Таким образом, в работе осуществлено численное решение модельной задачи конвекции в мантии Земли с применением последовательности сеток, демонстрирующее возможность
существенного сокращения времени вычислений за счет использования более совершенной
методики.
Автор благодарит Г. Г. Черных за помощь в постановке задачи и внимание к работе.
Список литературы
[1] Коновалов А. Н. Численное решение задачи теории упругости: лекции для студентов НГУ. Новосибирск: НГУ, 1968. 128 с.
[2] Тарунин Е. Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во ИГУ, 1990. 226 с.
[3] Тарунин Е. Л. Метод последовательности сеток для задач свободной конвекции // Журн. вычисл. матаматики и мат. физики. 1975. Т. 15, №2. C. 436-445.
[4] Федоренко Р. П. О скорости сходимости одного итерационного процесса // Журн. вычисл. математики и мат. физики. 1964. Т. 4, №3. C. 559-564.
[5] Червов В. В. Численное моделирование трехмерных задач конвекции в мантии Земли с применением завихренности и векторного потенциала // Вычисл. технологии. 2002. Т. 7, №1. С. 114-125.
[6] Полежаев В. И., Грязнов В. Л. Метод расчета граничных условий для уравнений Навье — Стокса в переменных вихрь, функция тока // Докл. АН СССР. 1974. Т. 219, №2.
[7] Роуч Х. Вычислительная гидродинамика. М.: Мир, 1980.
[8] Том А., Эйплт К. Числовые расчеты полей в технике и физике. М.: Энергия, 1964.
[9] Федорюк М. В. Характеристики течений несжимаемой жидкости в гравитационном поле // Мат. сб. 1988. №137(179), 4(12).
[10] 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.
[11] ЯнЕнко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ние, 1967.
[12] Blankenbach В., Busse F. et al. A benchmark comparison for mantle convection codes // Geophys. J. Int. 1989. Vol. 98. P. 23-38.
[13] Марчук Г. И., Шайдуров В. В. Повышение точности решения разностных схем. М.: Наука, 1979. 320 с.
Поступила в редакцию 31 января 2002 г.
