Моделирование нестационарного теплообмена в шаровом слое Земли Текст научной статьи по специальности «Физика»

^отве*,

Серия «Науки о Земле»

И З В Е С Т И Я

« О

2015. Т. 11. С. 106-116

Иркутского государственного университета

Онлайн-доступ к журналу: http://isu.ru/izvestia

УДК 550.362

Моделирование нестационарного теплообмена в шаровом слое Земли

С. В. Соловьев (solovjovsv@rambler.ru)

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

Ключевые слова: математическое моделирование, шаровая полость, конвективный теплообмен, электропроводная жидкость, магнитная гидродинамика, джоулева диссипация.

Введение

В настоящее время принято, что тепловая конвекция в земном ядре является причиной, приводящей к созданию геомагнитного поля Земли [4]. Теория геомагнитного поля получила название гидромагнитного динамо (ГД). Математическое исследование гидромагнитного динамо в общем виде еще далеко от завершения [1; 4], так как решение уравнений магнитной гидродинамики связано с большими трудностями. В этой связи теория ГД часто исследуется с применением кинематических моделей, в которых скорость движения жидкости считается заданной, а определяется только тепловое и магнитное поле.

В настоящей работе, в отличие от кинематических моделей, рассматриваются совместно уравнения: энергии, с учетом внутренних источников тепла и джоулевой диссипации; движения, с учетом магнитных, инерционных, вязких и подъемных сил; магнитной индукции, неразрывности для скорости и магнитной индукции. Используется приближение Буссинеска. Ускорение свободного падения направлено к центру сферической прослойки.

Описание модели

Математическая постановка задачи в переменных вихрь - функция тока - температура в безразмерной форме в сферической системе координат с учетом симметрии по долготе имеет вид [2]:

1 Эа

1

Но Эх Г 6

Эу Эа Эу Эа а Эу _ Эу

--------+ а с tg 6-

Э6 Эг Эг Э6 г Э6 Эг

Яё

Эа 2 Эа 1 Эа с tg 6 Эа

+--+ ^-т + -

ог 1 Э +

Эг2 г Эг Б

Э62

Яе2 г Э6 Яе„

Вг

Э2В6 „ В ЭВ6

6 + 2 Г 6

Эг 2

Э6 г 2 8Ш2 6

, + ЭВг ЭВ6 + В6 ЭВг

г Эг Эг Эг г Эг

_ (1)

В„ Э2В„ 1 ЭВ„ ЭВ„

+ В6 Э2Вв + 1 ЭВв ЭВ6 2В6 ЭВ6

г ЭгЭ6 г Эг Э6 г ЭгЭ6 г Эг Э6

Э6

В6 Э2Вг г2 Э62

ЭВ_ ЭВ^ Г2 Э6 Э6

Э2у 1 Э2у с tg 6 Эу . , -+ —--;---;--= - а г 81И (

Эг 2

г2 Э62

Э6

(2)

1 ЭЭ

1

ЭУЭЭ ЭУЭЭ

Но Эх г 8ш 6\Э6 Эг Эг Э6

1 ( Э2Э 2 ЭЭ 1 Э2Э

■ + -

Ре \ Эг2 г Эг г2 Э62

+сс^Э а

г2 Э6 ^

3 (ЭВ,+1В _ I ЭВг

2

Ре I Эг

г Э6

= 0;

(3)

ЭВГ_ Но Эх

1

Г 2 8Ш 6 1

В6 Э2 Ч

1 ЭВ6 ЭЧ

+—

г Э62 г Э6 Э6

+ В.,

Э2Ч ЭВ ЭЧ -+■

ЭгЭ6 Э6 Эг

Яе„

Э2 Вг

2 ЭВ,, 1 Э2В,, с tg 6ЭВГ +--- + ^-Г + —

Эг2 г Эг г2 Э62

2ВГ 2В6с tg 6 2 ЭВ6

Э6

Э6

1 ЭВ6

1

Но Эх Г 8Ш 6

1 ЭЧ ЭВ6

_ В

г Э6 Эг

Э2 Ч Эг 2

Э2 В6 Эг 2

В6

ЭВГ ЭЧ В6

Э2 Ч

_^___^_+ В6 ЭЧ

Эг Эг г ЭгЭ6 г2 Э6

2 ЭВ

г Эг г"

. + . 2 Э2Вг

.1 +с tg 6ЭВ6

Э62

Э6

Э6

(4)

(5)

г2$1п2 6

При записи системы дифференциальных уравнений использованы следующие обозначения: Э = (Т _Т2)1(Т1_Т2), В = В' / В0,Ч,а - безразмерные температура, магнитная индукция, функция тока и напряженность вихря; х = t /10 - безразмерное время; р0, и0, t0, В0 - характерные масшта-

бы; В' - размерная магнитная индукция; Qv - безразмерный внутренний источник теплоты; г = г' / г1 - безразмерный текущий радиус; г' - размерный текущий радиус; 6 - полярный угол; От - коэффициент магнитной вязкости (диффузии); о - электрическая проводимость жидкости.

Но = ^

Яе =

ип г

Ре =

и0 г

Ог =

я Р (Т - т2М3

Яет =

и0 г1 От

8 =

о Во2 г

Ро ио

- безразмерные числа гомохронности, Рейнольдса, Пекле, Грасгофа, магнитное число Рейнольдса, параметр магнитного взаимодействия. Остальные обозначения общепринятые.

Постоянная величина J, входящая в уравнение энергии (3), определяет величину джоулевой диссипации. И в зависимости от типа граничных условий для температуры принимает различные значения. Внутренний источник теплоты Qv задавался как равномерно распределенный по объему шаровой полости источник тепла постоянной мощности [1].

В данной работе при проведении вычислительного эксперимента для температуры задавались граничные условия первого типа, т. е. задавались значения температуры на внутренней и внешней поверхности слоя, которые принимали постоянные значения: Э| г = 1; Э| г = 0. Для граничных условий

J вычислялась по формуле:

первого

J =

От

типа

2

постоянная величина

вид:

4я 1(Т -Т2) Граничные

аэ

условия для температуры на оси симметрии имели

аб

= 0.

6 = 0,7

Граничные условия для функции тока, напряженности вихря и магнитной индукции имели следующий вид:

^ Г, , =1

6=0,Р=Ч=0,Р=°;

авг

аб

ав6

6 = 0,7

аб

В = В = 0; В6

- , - , = -0,0Ыи6; В6

г |Г. г |г, '6 |г. ' ' 6 |г,

= 0;

6=0,7

= 0,0Ыи6.

Граничные условия для вихря на границах сферического слоя предполагают линейное изменение его по нормали [3].

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

Локальные числа Нуссельта на поверхности внутренней и наружной сферы рассчитывались по формулам:

аэ аг

Ыи2 = -Я2

аэ аг

2

а

V

V

Для расчета осредненных чисел Нуссельта (среднеинтегральные значения) производилось осреднение локальных чисел Нуссельта по поверхностям г = 1 и г = Я2:

Метод решения

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

Алгоритм решения задачи следующий:

1. Разбиение и нумерация узлов области.

2. Расчет безразмерных чисел (критериев) подобия.

3. Выбор граничных условий.

4. Задание начальных условий (в стационарном случае начальных приближений).

5. Интегрирование невязок и построение системы нелинейных уравнений.

6. Проведение одной итерации решения системы. Получение новых значений расчетных полей: температуры, вихря, функции тока, магнитной индукции.

Далее возможны два случая:

- стационарный режим:

7. Сравнение величины полученной погрешности с заданной и, если необходимо, возврат к п. 6. При достижении заданной величины погрешности для всех расчетных полей переход к п. 8.

- нестационарный режим

7. Переход к следующему шагу по времени и возврат к п. 6. При достижении заданной продолжительности времени расчета переход к п. 8.

8. Окончание решения.

Необходимо заметить, что при расчете на каждой «внешней» итерации новых значений температуры, вихря и магнитной индукции для расчета функции тока из уравнения Пуассона (2) осуществляется «внутренний» итерационный цикл до достижений заданной степени точности. Это осуществляется независимо от расчета стационарного или нестационарного режима.

_ 1р аэ _

Ыи1 =—| - 6ё6, Ыи.

201_аг _

2 = - ^.г ^ 81П 6^6.

2 01_аг]г2

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

Результаты

На рисунках 1-6 приведены результаты нестационарных расчетов полей температуры, функции тока, напряженности вихря, радиальной и меридиональной составляющей магнитной индукции, а также локальных чисел Нуссельта на внутренней и внешней границе сферической прослойки при следующих значениях безразмерных чисел подобия: Ог = 104;Яе = 10;5 = 10-5;Яе = 1;Ре = 10;г2 /г = 2,5;О = Но = 1.

' ' ' т ' '21''

На рисунке 1 приведены результаты расчетов поля температуры для следующих моментов безразмерного времени т: 0,1; 0,5; 0,7; 1,0; 1,5; 2,0; 2,5; 3,0; 3,5; 5,0. Дальнейшее увеличение времени расчета не приводило к качественным изменениям поля температуры. При этом имели место небольшие количественные различия значений температуры внутри расчетной области по сравнению со стационарным режимом. Приведенные выше моменты времени соответствуют позициям а-к на рис. 1.

е ж з и к

Рис. 1. Эволюция поля температуры

На начальной стадии процесса (т = 0,1) небольшое изменение температуры происходит лишь в области границ слоя (рис. 1, а). С увеличением времени (т = 0,5; 0,7; рис. 1, б, в) изменение температуры происходит уже в экваториальной области. Причем поле температуры приобретает вид, характерный для конвективного теплообмена (линии изотерм далеки от концентрических окружностей, характерных для механизма теплопроводности). По

мере прогрева жидкости в слое поле температуры напоминает форму ушной раковины (рис. 1, г, д) без «перемычки» в области экватора, которая («перемычка») в дальнейшем «срастается» с внешней границей слоя (рис. 1, е), а затем исчезает (рис. 1, ж-к), достигая при т = 13 стационарного режима.

На рисунке 2 приведены результаты расчетов поля функции тока для следующих моментов безразмерного времени т: 0,1; 0,5; 0,7; 1,0; 1,5; 2,0; 2,5; 3,0; 3,5. Дальнейшее увеличение времени расчета не приводило к качественным изменениям поля функции тока. При этом имели место незначительные различия значений функции тока внутри расчетной области по сравнению со стационарным режимом.

В шаровой полости образуются две конвективные ячейки (рис. 2, а-и). Оказалось, что в северном полушарии жидкость движется по часовой стрелке (значения функции тока отрицательные), а в южном - против часовой стрелки (значения функции тока положительные).

е ж з и

Рис. 2. Эволюция поля функции тока

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

В шаровой полости образуются два вихря (рис. 3, а-и). В северном полушарии жидкость так же, как и в случае поля функции тока, движется по часовой стрелке (значения напряженности вихря отрицательные), а в южном - против часовой стрелки (значения напряженности вихря положительные). Как видно из результатов, представленных на рис. 3, на всем времен-

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

а 6 в г д

е ж з и

Рис. 3. Эволюция поля напряженности вихря

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

Дальнейшее увеличение времени расчета не приводило к качественным изменениям поля радиальной составляющей магнитной индукции. При этом имели место незначительные различия значений радиальной составляющей магнитной индукции внутри расчетной области по сравнению со значениями для стационарного режима.

Для всего времени расчета радиальная составляющая магнитной индукции (рис. 4) в северном полушарии принимает отрицательные значения, за исключением небольшой области вблизи внутренней поверхности сферического слоя, а в южном - положительные, за исключением небольшой области вблизи внутренней поверхности сферического слоя. Причем форма этих небольших областей, вблизи внутренней поверхности сферического слоя, в начальный период времени вытянута вдоль самой поверхности (рис. 4, а, б). С течением времени эта область трансформируется в два небольших вихря (рис. 4, в-и), интенсивность которых растет. В северном полушарии значения радиальной составляющей магнитной индукции положительные, а в южном - отрицательные.

а 6 в г д

е ж з и

Рис. 4. Эволюция поля радиальной составляющей магнитной индукции

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

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

На всем временном интервале расчета меридиональная составляющая магнитной индукции (рис. 5) принимает положительные значения практически во всей расчетной области, за исключением небольшой области вблизи внутренней поверхности сферического слоя, в которой значения меридиональной составляющей магнитной индукции отрицательные. С течением времени форма отрицательной области видоизменяется, приближаясь к форме «треугольника» с вершиной на экваторе вблизи внешней поверхности слоя (рис. 5, в-и). Область положительных значений также изменяется (рис. 5, в-и).

На рисунке 6 приведены результаты расчетов локальных чисел Нус-сельта на внутренней и внешней поверхности сферической прослойки для следующих моментов безразмерного времени т: 0,1; 0,5; 0,7; 1,0; 1,5; 2,0; 3,5 (моменты времени соответствуют позициям а-ж на рис. 6).

е ж з и

Рис. 5. Эволюция поля меридиональной составляющей магнитной индукции

1

Рис. 6. Эволюция распределения локальных чисел Нуссельта

Дальнейшее увеличение времени расчета не привело к качественным изменениям профилей локальных чисел Нуссельта. Небольшие количественные различия имели место при стремлении расчетных результатов к их стационарному значению. На рисунке 6 кривой с номером 1 соответствует распределение локальных чисел Нуссельта на внутренней поверхности сферического слоя, а кривой с номером 2 - распределение локальных чисел Нуссельта на внешней поверхности.

Начальная стадия процесса (т = 0,1) характеризуется высокими значениями локальных чисел Нуссельта на всей внутренней поверхности сферического слоя (от северного полюса до южного) по сравнению со значениями локальных чисел Нуссельта на внешней поверхности (рис. 6, а). При значении полярного угла 6 »я/2 (экваториальная плоскость) распределение локальных чисел Нуссельта на внутренней поверхности имеет слабо выраженный минимум (рис. 6, а; кривая 1), а распределение локальных чисел Нуссельта на внешней поверхности - ярко выраженный максимум (рис. 6, а; кривая 2).

Причем значения локальных чисел Нуссельта на внешней границе в области полюсов близки к нулю. Локальные числа Нуссельта изменяются соответственно на внутренней и на внешней поверхности расчетной области

в интервале: 7,40 < Ыи1 < 7,80; Щ= 7,58; 0,19 < Ыи2 < 4,81; ~Ш2 = 2,67.

Далее ситуация изменяется на противоположную с качественно одинаковым характером распределения локальных чисел Нуссельта (рис. 6, б-ж) как на внутренней поверхности, так и на внешней. Здесь при значении полярного угла 6 »я /2 распределение локальных чисел Нуссельта на внутренней поверхности имеет ярко выраженный минимум (рис. 6, б-ж; кривая 1), а распределение локальных чисел Нуссельта на внешней поверхности -ярко выраженный максимум (рис. 6, б-ж; кривая 2). Для этих результатов значения локальных чисел Нуссельта на внутренней поверхности в области полюсов близки к своим максимальным значениям, а значения локальных чисел Нуссельта на внешней поверхности - к минимальным.

Локальные числа Нуссельта изменяются соответственно на внутренней и на внешней поверхности расчетной области в интервале:

- для результатов рис. 6, б:

1,78 < Щ < 3,92; Ыи1= 3,19; 0,49 < Ыы2 < 4,11; ~ый2 = 2,61;

- для результатов рис. 6, ж:

1,39 < Щ < 7,47; Ыи1 = 4,73; 0,67 < Ыы2 < 12,57;Щ, = 8,39.

Выводы

Предложена математическая модель, позволяющая моделировать нестационарный конвективный теплообмен и магнитную гидродинамику электропроводной жидкости (с учетом внутренних источников тепла, джо-улевой диссипации, инерционных и подъемных сил) в шаровом слое. Из анализа полученных результатов можно сделать следующие выводы:

• для рассмотренных граничных условий и значений безразмерных критериев подобия в расчетной области имеет место двухъячеистое течение (т. е. две конвективные ячейки);

• с течением времени поля температуры, функции тока, напряженности вихря, радиальной и меридиональной составляющей магнитной индукции и локальных чисел Нуссельта претерпевают значительные как качественные, так и количественные изменения;

• математическая модель и полученные результаты могут быть полезными для исследования тепловых и гидродинамических процессов, происходящих в недрах Земли и других планет, а также проследить и объяснить их тепловую и гидродинамическую эволюцию.

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

1. Жарков В. Н. Физика Земли и планет. Фигуры и внутреннее строение / В. Н. Жарков, В. П. Трубицын, Л. В. Самсоненко. - М. : Наука, 1971. - 384 с.

2. Соловьев С. В. Моделирование конвективного теплообмена в жидком ядре Земли // Вестн. Тихоокеан. гос. ун-та. - 2012. - № 2 (25). - С. 73-82.

3. Численные методы исследования течений вязкой жидкости / А. Д. Госмен, В. М. Пан, А. К. Ранчел, Д. Б. Сполдинг., М. Вольфштейн. - М. : Мир, 1972. - 320 с.

4. Яновский Б. М. Земной магнетизм / Б. М. Яновский. - Л. : Изд-во Ленингр. ун-та, 1978. - 592 с.

Simulation of Unsteady Heat Transfer Electroconductive Fluid in a Spherical Cavity

S. V. Solovjov

Abstract. In this paper the system of equations describing the convective heat electroconductive liquid with regard to the internal heat sources and Joule dissipation in a spherical cavity simulating liquid core of the Earth is considered. The evolution of the flow structure, temperature field, magnetic field, and the Nusselt numbers is investigated.

Keywords: mathematical modeling, spherical cavity, convective heat transfer, electro-conductive liquid, magnetic hydrodynamics, Joule dissipation.

Соловьев Сергей Викторович доктор физико-математических наук профессор

Тихоокеанский государственный университет 680035, г. Хабаровск, ул. Тихоокеанская, 136 тел. (4212) 37-51-88

Solovjov Sergei Viktorovich

Doctor of Sciences (Physics and

Mathematics), Professor

Pacific National University

136, Tikhookeanskaya st., Khabarovsk,

680035

tel.: (4212) 37-51-88