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

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

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

2015. Т. 12. С. 94-114

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

И З В Е С Т И Я

УДК 550.362

Моделирование конвективного теплообмена электропроводной жидкости в шаровой полости. Часть I

С. В. Соловьев ([email protected])

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

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

Введение

Вопросы о тепловом состоянии Земли, распределении источников тепла в ее глубинах, геомагнитном поле Земли и его вариациях имеют фундаментальное значение для любых гипотез о строении и эволюции Земли [1]. Что касается ядра Земли, то распределение температуры в нем еще недостаточно изучено. Интерес к ядру Земли возрос в связи с проблемой происхождения земного магнитного поля, которое оказывает существенное влияние на многие процессы глобального характера, происходящие в Земле. Считается, что магнитное поле Земли создается при магнитогидродинамических течениях в земном ядре [1; 2; 5]. Современные теории геомагнетизма исходят из предположения, что магнитное поле Земли создается и поддерживается за счет динамо-механизма [5]. Математическое исследование гидромагнитного динамо в общем виде еще нуждается в продолжении [1; 2; 5], так как решение уравнений магнитной гидродинамики связано с трудностями. Часто теория гидромагнитного динамо исследуется с применением кинематических моделей, в которых поле скорости жидкости задается, а рассчитывается только тепловое и магнитное поле.

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

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

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

1

Г 6

Эу Эа Эу Эа а Эу п Эу

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

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

_

_ Йё ОГ 1 ЭЗ

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

+--+ ^-г + -

■ + -

Эг 2 5

Йе2 г э6 Йе

Б„ Э2Б.

Б

г Эг Э2 Б6

Э62

Э6 Г2 8т2 6

Эг 2

+ 2

В.. ЭБ6 ЭВ. ЭБ,

г Эг

2

+ -

6 + Б6 ЭБг

Эг Эг г Эг

1 ЭБг ЭБ. + ^дБ +12Б6 ЭБ6

г Э.Э6 . Эг Э6 г Э.Э6 . Э. Э6

Э6

Б6 Э2Б

_6__г_

.2 Э62

1_ ЭА

.2 Э6 Э6

(1)

Э у 1 Э у с tg 6 Эу . , ■ + —--;---°--_ - а . 81И 6

ЭГ2 .2 Э62

Э6

(2)

1

ЭуЭЗ Эу ЭЗ

. 81п6\Э6 Э. ЭГ Э6 с tg 6ЭЗ

1 Г Э2З 2 ЭЗ 1 Э2З

Ре {Эг 2 + . Эг + г2 Э62

Э6

+а

3 (ЭБ 1 1 ЭБг + — Б --

Ре I Эг

г Э6

_ 0:

(3)

0 _-

1

Г 8Ш 6

Б6 Э2Ч 1 ЭБ6 ЭЧ „ Э2Ч

—-т +--1-+ Б,-+. Э62 , Э6 Э6

ЭБ, ЭЧ

Э.Э6 Э6 Эг Э2Б., 2 ЭБ., 1 Э2Б г с tg 6ЭБг

-L +--!- +--L +-Е--Г--

Э.2 Г Эг г2 Э62 г2 Э6

- Щ_ _ 2Б6с ^ 6 - 2_

- ~ Г2 Г ~86

(4)

т

+

0 = ■

1

Г 6

1ЭЧ ЭБ6 г Э6 Эг

- Б,

Э2Ч ЭБ„ ЭЧ

Эг Эг Эг

Яет

Б6 Э2Ч

г ЭгЭ 6

2

АЭЧ.

г2 Э6

- +

Э2Б6 + 2 ЭБ6 + 1 Э2Б6 + с гв 6 ЭБ6

"ЭГ2

Б6

Г2 8Ш2 6

г Эг + 2 Э2Б,

г2 Э62

Э6

г" Э6

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

Qv = ЧГ-- безразмерный внутренний источник теплоты; д - размерный

Ч "

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

яРЧг14

Яе = и^

ре=и^

а, =-

V а \\ Пт р0 и0

мерные числа Рейнольдса, Пекле, Грасгофа, магнитное число Рейнольдса, параметр магнитного взаимодействия. Остальные обозначения общепринятые. Постоянная величина J, входящая в уравнение энергии (3), определяет величину джоулевой диссипации и в зависимости от типа граничных условий для температуры принимает различные значения.

В данной работе при проведении вычислительного эксперимента для температуры задавались следующие граничные условия: на внутренней поверхности полости Г1 ( г = 1) граничное условие второго рода (тепловой поток по закону Фурье) [2], а на внешней поверхности Г2 ( г = я 2) - граничное условие первого рода (постоянное значение температуры): Э9

Яет = и* т В

^ = а Бо Г1 - безраз-

Э п

= 1; 9 г = 0.

Для граничных условий второго рода на внутренней поверхности полости постоянная величина J вычислялась по формуле:

В г

J = А. 4л д

Г

Э

Граничные условия для температуры на оси симметрии имели вид: = 0.

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

Y г» =Y

, ЭБ

_ I _ А , Г

е=о,p = we=o,p = 0 ' Э

ЭБа

e=o,i

эе

= 0;

е=о,т

Бг\Г1 = Бг\Гг = 0; Бе|г1 = -0,01sinе; Бе|г— = 0,01sinе .

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

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

дЗ

эз

Nu1 =--

Эг

Г1

Nu2 =- R2-

Эг

Осредненные числа Нуссельта вычислялись по формулам:

_ 1 p

NU1 =-—j 2 0

эз

Эг

R—

sin еёе, Nu2 =—2-j 2 0

эз

Эг

sin её е.

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

Численное решение задачи осуществлялось с помощью метода конечных элементов на равномерной сетке г = 1 + / -Дг; 9у = у -Д9; / = 0, N г;

у = 0, Ы9 с количеством узлов в области Ъ = Nг - N9; Дг = ——1, Д9 = -Р- - шаг

N.

N

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

w(r, е) = t waNa (г, е); Y(r, е) = t YaNa (г, е); 3(г, е) = t заК (г, е);

Бг (г, е) = ±БгaNa (г, е); Бе (г, е) = £БеаN (г, е);

а=1

Z

£

а=1

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

Дискретный аналог уравнения (1) имеет следующий вид:

£ £ [^у^в^йв + ввьвв + в^в^у о;

b=1 g=1

е=0,р

где

1 Г2 Л

Т;р = -

!\_С 1 0

N

ЭЖ„ ЭЖ„ ЭЖ„ 1 ЭЫ„ эЖ

Эг

--Г-

Эг Эг г Э6 Э6

+N;

г Э6

- N.

N

Г 81П

ё6ёг.

\л N

ТИУ _ Г Г ;

"ру _ J J • г

10 Г 81п(

ЭNУ ЭNP ЭNУ ЭNP 1 Э^

-у--Р---у--^-^к-^ + с 6М, у

Э6 Эг Эг Э6 г р Э6 р

Э^у

'р Эг

й6йг;

/

^ Яе2 ■ ле 1 0

Яет Ю Эг р Э6

ТБГБ6 _ _

5 Г2 „ ЭN, ЭNa „ ЭN, _ ЭNр 1 Э^ Э^ Л

Яе 1

т 1

N« N р-

т 10 у

ТБ6 Б6 ^«ру

Эг

^ р

Эг Эг

■+ N«

—Р N, N,

Эг у г Э6 Э6 у

ё6йг.

/

5

Яет 1

-2 " и

-N. N

ЭNу ЭN _

ч г ; р Э6 Эг р Э6 0

ё 6ёг .

Дискретный аналог уравнения (2) имеет следующий вид:

£ м>„+£ м>„ = 0;

где М «р =- Я

n; + г

Эг

Эг г Э6 Э6

Э6

ё6ёг;

М;р=шр.Я/ Г2N.N 81П 66^ .

Дискретный аналог уравнения (3) имеет следующий вид:

£ 9р КР+££

Р=1 Р=1 у=1

9рУу К;9ру + Б6Р Б6у ^ +

+ Б Б ТБ6Б, + Б Б ТБ6Б6

,р 6у «Ру 6р 6у «Ру

= К;;

1 Г2 Л

где К;р=- Р-ц 1 0

# Гсгв6ЭNР

Эг

Эг Эг г Э6 Э6

Э6

ё 6йг,

Г2 Л лг ¡9 Г Г

к ;ру = я я -

10 Г 81п(

ЭNу ЭNр ЭNу ЭNр

^ Э6 Эг Эг Э6 0

(¡6ёг;

¡^Б6 Б6 л;Ру

Г Г2 л

— Л ГN; Ре Ю ;

ЭNР 1 Эг г

ЭN 1

^ + ^у

Эг г

ё6ёг ;

Т Г2 Л

кб6б, = —— Я я ^

«ру п^ ^^ 1 0

Т Г2 Л

КБГБ6 =-— \ Я ^

«ру р„ JJ ; ^^ 1 0

ЭNр 1

+ - Nр

Эг г 1 Э^

1 ЭN,

г Э6

ё6ёг;

г Э6

ЭN 1

+ ^у

Эг г

ё 6йг ;

J r2 Р

KBrBr =-J [ ) rN ^abg pJ J a

re 1 0

1 SNñ

r SQ

1 SNg

r SQ

dQdr , Ka= Щ) rNadQdr . Pe 1 0

Дискретный аналог уравнения (4) имеет следующий вид:

£ Bqb SBb + ± Br, SBb +± ± [BepyY SBQg + S*^ 0;

где

1 r2 Р

SB = - ~RT))

1 0

Na

b=1 g=1

SNR SNa SNь ■-r-

Sr Sr Sr

1 SNa SNb + c tg e N SNl

r Se Se

-2a NaNb

Se r sin e

dedr,

1 r2 Р

SB = - rT"))

1 0

r2 p ( 2 SN ( л

-N„

sq

dedr, ^ = -)) Jí-^^rnr: abg 1J0 r sin e Sr se

r2 p N SN SN ^Я ^^^dQdr .

1 0 sin Q Sr Sr Дискретный аналог уравнения (5) имеет следующий вид:

±B RBQ +±B RBr +±±[B У RBQy + B У RBry] = 0-

b=1 b=1 b=1 g=1

где

r2 Р

RBr =-— 2 )

л*ь Re J J

1 0

rSK Snl-n 1SKSN.-n Snl+2nn

Sr Sr

Sr r se se

se r

ab

dQdr,

r2 Р

RBe =— ) ) л*ь Re JJ

1 0

A 2c tg e 2 SN„Л

^^ N a N p+ ± N a—Qr

v r r se 0

dQdr;

RaBíMJ

1 0

fr 1 SNa лт cose лт „^SN. л Nb--— Na N

r sin e se

r sin

ab

j Sr y

dQdr;

SK.N - c0se NN

2 • Q b 2 • 2n alyb

r sin e se r sin e

SN g

0 se 0

dQdr.

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

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

индукции для расчета функции тока из уравнения Пуассона (2) осуществляется «внутренний» итерационный цикл до достижения заданной степени точности (10-9).

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

Результаты

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

На рисунке 1 приведены результаты расчетов поля температуры.

а б в г д е ж

Рис. 1. Поле температуры: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил, 3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

На рисунке 1, а представлено поле температуры для неэлектропроводной жидкости, изотермы которого - концентрические окружности, т. е. перенос энергии в полости осуществляется теплопроводностью. Максимальное значение температуры в полости 7шах = 5,208. Такая же ситуация сохраняется и для электропроводной жидкости (рис. 1, б) без учета джоуле-вой диссипации и внутренних источников тепла. Максимальное значение температуры в полости принимает то же значение, что и для результата рис. 1, а. Учет джоулевой диссипации (рис. 1, в) изменят механизм переноса энергии в полости с кондуктивного на конвективный. Изотермы уже отличаются от концентрических окружностей. Максимальное значение температуры в полости Тшах = 6,917. Учет внутренних источников (Qv = 1) и стоков

(Qv = -1) тепла (рис. 1, г, д), совместно с учетом джоулевой диссипации, также приводит к конвективному механизму теплообмена. Максимальное значение температуры в полости для результата рис. 1, г Ттах = 10,177, а для результата рис. 1, д диапазон изменения температуры [-0,498; 3,662]. Неучет джоулевой диссипации, несмотря на наличие внутренних источников (стоков) тепла (рис. 1, е, ж), изменяет механизм переноса энергии в полости с конвективного на кондуктивный. Изотермы представляют собой концентрические окружности. Максимальное значение температуры в полости для результата рис. 1, е Ттах = 8,436, а для результата рис. 1, ж диапазон изменения температуры [-0,655; 1,981].

На рисунке 2 приведены результаты расчетов локальных чисел Нус-сельта на внутренней и внешней поверхности шаровой полости.

Рис. 2. Распределение локальных чисел Нуссельта: а - без учета магнитных сил, 3 = 0, Qy = 0; б - с учетом магнитных сил,

3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных

сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом

магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

На рисунке 2 линии с номером 1 соответствует распределение локальных чисел Нуссельта на внутренней поверхности полости, причем это распределение совпадает с осредненным и принимает постоянное значение

Ыы1 = Ыы1 = 10 (здесь и далее) согласно заданному граничному условию для температуры. Поэтому далее будут анализироваться распределения локальных чисел Нуссельта только на внешней поверхности полости (кривая с номером 2). Распределения локальных чисел Нуссельта на внешней поверхности, приведенные на рис. 2, а, б, совпадают с осредненными значениями и принимают постоянные значения Ыы2 = Л^м21 = 5,319. Это согласуется с распределением температуры (рис. 1, а, б), когда теплообмен осуществляется теплопроводностью. Учет джоулевой диссипации (рис. 2, в) значительно изменяет распределение локальных чисел Нуссельта по сравнению с предыдущими результатами. Такое распределение характерно для конвективного теплообмена и поля температуры, представленного на рис. 1, в. При значении полярного угла 9 »р /2 (экваториальная плоскость) распределение локальных чисел Нуссельта имеет ярко выраженный максимум (рис. 2, в). Значение ос-редненного и диапазон изменения локальных чисел Нуссельта следующие: Ыы2 = 9,646; 5,894 < Ыи2 < 12,386. Учет внутренних источников = 1) и стоков (Qv = -1) тепла (рис. 2, г, д) совместно с учетом джоулевой диссипации, также приводят к распределениям локальных чисел Нуссельта, характерных для конвективного теплообмена. При наличии внутренних источников тепла значения локальных чисел Нуссельта на внешней поверхности больше значений числа Нуссельта на внутренней поверхности (рис. 2, г). Для стоков тепла (рис. 2, д) ситуация обратная. Причем локальные числа Нуссельта принимают отрицательные значения. Для значения полярного угла 9»р /2 распределения локальных чисел Нуссельта имеют максимум (рис. 2, г, д; кривая 2). Значение осредненного и диапазон изменения локальных чисел

Нуссельта для результата рис. 2, г: Ыы2 = 21,030; 17,042 < Ыы2 < 23,929; а

для результата рис. 2, д: Ыы2 = -1,739; -5,283 < Ыы2 < 0,863. Неучет джоулевой диссипации, несмотря на наличие внутренних источников (стоков) тепла (рис. 2, е, ж), приводит к распределению локальных чисел Нуссельта, характерному для режима теплопроводности (рис. 1, е, ж). Распределения локальных чисел Нуссельта совпадают с осредненными значениями и принимают постоянные значения для результата рис. 2, е: Ыы2 = Ыы2 = 16,704; для

результата рис. 2, ж: Ыи2 = Ыы2 = -6,065.

На рисунке 3 приведены результаты расчетов поля функции тока.

а б в г д е ж

Рис. 3. Поле функции тока: а - без учета магнитных сил, J = 0, Qv = 0; б - с учетом магнитных сил, J = 0, Qv = 0; в - с учетом магнитных сил, J ^ 0, Qv = 0; г - с учетом магнитных сил, J ^ 0, Qv = 1; д - с учетом магнитных сил, J ^ 0, Qy = -1; е - с учетом магнитных сил, J = 0, Qy = +1; ж - с учетом магнитных сил, J = 0, Qv = -1

Для всех случаев в шаровой полости образуются две конвективные ячейки. На рисунке 3, а представлено поле функции тока для неэлектропроводной жидкости, которая в северном полушарии полости движется по часовой стрелке (значения функции тока отрицательные), а в южном - против часовой стрелки (значения функции тока положительные). Интенсивность движения жидкости в ячейках незначительная (режим теплопроводности), максимальное значение функции тока |^тах| = 3,05 10-8. В случае электропроводной жидкости (рис. 3, б) без учета джоулевой диссипации и внутренних источников тепла направление течения жидкости в ячейках меняется на противоположное по сравнению с результатом рис. 3, а. В северном полушарии жидкость движется против часовой стрелки (значения функции тока положительные), а в южном - по часовой стрелке (значения функции тока отрицательные). Изменение направления циркуляции жидкости в конвективных ячейках можно интерпретировать как изменение магнитных полюсов в шаровой полости. Максимальное значение функции тока незначительно возрастает |^тах| = 2,6210-7. Учет джоулевой диссипации (рис. 3, в), так же как и учет внутренних источников (Qv = 1) и стоков (Qv = -1) тепла (рис.

3, г, д), совместно с учетом джоулевой диссипации, по сравнению с результатом рис. 3, б изменяет направление циркуляции жидкости в конвективных ячейках на противоположное. Максимальное значение функции тока возрастает и для результатов рис. 3, в-д принимает следующие значения: |^тах| = 7,2310-3; 7,7610-3; 6,7710-3. Неучет джоулевой диссипации, несмотря на наличие внутренних источников (стоков) тепла (рис. 3, е, ж), приводит к изменению направления циркуляции жидкости в ячейках на противоположное, по сравнению с результатами рис. 3, в-д. Максимальное значение функции тока для результатов рис. 3, е, ж принимает соответственно следующее значение: |Утах| = 2,78-10-7; 2,47-10-7.

На рисунке 4 приведены результаты расчетов поля напряженности вихря.

б

д

Рис. 4. Поле напряженности вихря: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил, 3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv =-1

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

= 2,48-10

,63-10 ; 9,27-10 ; 8,08-10"

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

составляет величину |йшах| = 3,8910-6; 4,0410-6; 3,7610-6 соответственно для

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

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

а

в

г

е

а б

Рис. 5. Поле магнитной индукции: а - радиальной составляющей; б - меридиональной составляющей

Максимальная интенсивность радиальной составляющей магнитной индукции не превышает величину |Вгтах| = 4,83-10-4, а меридиональной -

|£етах| = 1,0010-2.

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

На рисунках 6-10 приведены результаты для Ог = 103. Значения остальных безразмерных критериев подобия те же, что и выше.

а б в г д еж

Рис. 6. Поле температуры: а - без учета магнитных сил, J = 0, Qv = 0; б - с учетом магнитных сил, J = 0, Qv = 0; в - с учетом магнитных сил, J ^ 0, Qv = 0; г - с учетом магнитных сил, J ^ 0, Qv = 1; д - с учетом магнитных сил, J ^ 0, Qv = -1; е - с учетом магнитных сил, J = 0, Qy = +1; ж - с учетом магнитных сил, J = 0, Qy = -1

На рисунке 6 приведены результаты расчетов поля температуры.

На рисунке 6, а представлено поле температуры для неэлектропроводной жидкости, изотермы которого уже отличны от концентрических окружностей. Перенос энергии в полости осуществляется конвекцией. Максимальное значение температуры в полости Ттах = 5,904. Для электропроводной жидкости (рис. 6, б) поле температуры характерно для режима теплопроводности: изотермы - практически концентрические окружности. Максимальное значение температуры в полости незначительно уменьшается Ттах = 5,276. То есть учет магнитных сил подавляет конвекцию, сводя конвективный теплообмен к режиму кондукции. Учет джоулевой диссипации (рис. 6, в) усиливает конвекцию в полости. Максимальное значение температуры в полости Ттах = 6,956. Учет внутренних источников (Qv = 1) и стоков (Qv = -1) тепла (рис. 6, г, д) также приводит к конвективному механизму теплообмена. Максимальное значение температуры в полости для результата рис. 6, г Ттах = 9,370, а для результата рис. 6, д диапазон изменения температуры [-0,751; 3,936]. Неучет джоулевой диссипации (рис. 6, е) также сохраняет конвективный механизм теплообмена при учете внутренних источников тепла, а при учете внутренних стоков тепла (рис. 6, ж) механизм переноса энергии в полости изменяется на кондуктивный. Изотермы представляют собой концентрические окружности. Максимальное значение температуры в полости для результата рис. 6, е Ттах = 8,795, а для результата рис. 6, ж диапазон изменения температуры [-0,635; 2,064].

На рисунке 7 приведены распределения локальных чисел Нуссельта. Характер изменения чисел Нуссельта связан с полем температуры. В режиме теплопроводности распределение локальных чисел Нуссельта - это прямая линия 1 (Ии1 = Ии1 = 10), а в режиме конвекции - кривая линия 2.

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

локальных чисел Нуссельта следующие: Ыи2 = 5,319; 4,592 < Ыи2 < 7,982. Учет магнитных сил изменяет характер распределения локальных чисел Нуссельта (рис. 7, б; линия 2). Распределение локальных чисел Нуссельта практически совпадает с осредненным значением (режим теплопроводности). Значение осредненного и диапазон изменения локальных чисел Нуссельта

следующие: Ыи2 = 5,405; 5,387 < Ыи2 < 5,443. Учет джоулевой диссипации (рис. 7, в) значительно изменяет распределение локальных чисел Нуссельта. При значении полярного угла в»п /2 оно имеет ярко выраженный максимум (рис. 7, в). Значение осредненного и диапазон изменения локальных чисел Нуссельта следующие: Ыи2 = 9,702; 2,475 < Ыи2 < 17,926. Учет источников = 1) и стоков = -1) тепла (рис. 7, г, д) также приводит к распределению локальных чисел Нуссельта, характерных для конвективного

теплообмена. Оба распределения имеют максимум при в»п /2. Однако при наличии источников тепла линия 2 (рис. 7, г) лежит выше прямой 1, а при наличии стоков (рис. 7, д) - ниже. Значение осредненного и диапазон изменения локальных чисел Нуссельта для результата рис. 7, г: Ыы2 = 21,012; 10,667 < Ыи2 < 35,530; а для результата рис. 7, д: Шы2 = -1,608; -6,139 < Ыи2 < 2,228. Неучет джоулевой диссипации (рис. 7, е) изменяет распределение локальных чисел Нуссельта, кривая 2 располагается выше прямой 1 и имеет один максимум и два минимума. Значение осредненного и

диапазон изменения локальных чисел Нуссельта следующие: Ыи2 = 16,714; 11,700 < Ыи2 < 23,823. При учете стоков тепла (рис. 7, ж; кривая 2) распределение локальных чисел Нуссельта лежит в области отрицательных значений и соответствует режиму теплопроводности: Ыи2 = Ыи2 = -5,903.

д е ж

Рис. 7. Распределение локальных чисел Нуссельта: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил,

3 = 0, Qy = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных

сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом

магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

На рисунке 8 приведены результаты расчетов поля функции тока.

б в где ж

Рис. 8. Поле функции тока: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил, 3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

Для всех случаев, за исключением режима, представленного на рис. 8, е, в шаровой полости образуются две конвективные ячейки. Для результатов рис. 8, а, б, ж в северном полушарии полости жидкость движется против часовой стрелки (значения функции тока положительные), а в южном - по часовой стрелке (значения отрицательные). Учет джоулевой диссипации, как с источниками (стоками) тепла, так и без них (рис. 8, в-д), изменяет направление циркуляции жидкости в конвективных ячейках. Для режима, представленного на рис. 8, е, в области полюсов образуются еще две мелкомасштабные конвективные ячейки малой интенсивности. Максимальное значение функции тока для результатов рис. 8, а-ж следующие: |Чшах| = 3,8610-2; 8,40 10-4; 2,12-10-1; 3,0010-1; 1,22-10-1; 1,7510-1; 4,07 10-7.

На рисунке 9 приведены результаты расчетов поля напряженности вихря.

Г 1Г 15

а б в г д еж

Рис. 9. Поле напряженности вихря: а - без учета магнитных сил, 3 = 0, Qy = 0; б - с учетом магнитных сил, 3 = 0, Qy = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

а

В полости (рис. 9, а, б, ж) образуются два крупномасштабных вихря, в которых жидкость движется против часовой стрелки (значения положительные) в северном полушарии, и по часовой стрелке (отрицательные значения) - в южном. Максимальная интенсивность этих вихрей |ютах| = 6,2710-1;

1,1210-2; 5,2210-6 соответственно для результатов рис. 9, а, б, ж.

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

Максимальная интенсивность этих вихрей составляет величину |®тах| = 2,68; 3,78; 1,53 соответственно для результатов рис. 9, в-д. Неучет

джоулевой диссипации при наличии внутренних источников тепла в электропроводной жидкости приводит к образованию в полости четырех мелкомасштабных вихрей (рис. 9, е). Направление циркуляции жидкости в этих вихрях меняется четыре раза. Максимальная интенсивность вихрей не превосходит величину |ютах| = 2,37.

На рисунке 10 приведены результаты расчетов поля радиальной и меридиональной составляющей магнитной индукции.

а б

Рис. 10. Поле магнитной индукции: а - радиальной составляющей; б - меридиональной составляющей

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

ставляющей магнитной индукции не превышает величину |Вгтах| = 4,93 10-4, а меридиональной - |В0тах| = 1,00 10-2. Можно констатировать, что поле магнитной индукции аналогично результатам, полученным для числа Грасгофа Ог = 102.

Результаты, иллюстрирующие дальнейшее увеличение интенсивности теплообмена в полости при значении числа Грасгофа Ог = 104, представлены на рис. 11-15. Значения остальных безразмерных критериев подобия не изменялись.

На рисунке 11 приведены результаты расчетов поля температуры.

а б в где ж

Рис. 11. Поле температуры: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил, 3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

Для всех режимов теплообмен в полости происходит путем конвекции. Для результатов рис. 11, а, б изменение температуры происходит в области полюсов и экватора. Затем (рис. 11, в, г) основное изменение температуры происходит в области экватора. И для результатов, представленных на рис. 11, д-ж, изменение температуры опять наблюдается в области полюсов и экватора. Максимальное значение температуры Ттах в полости и интервал изменения температуры в полости для результатов, представленных на рис. 11, д, е, следующие: 3,630; 3,630; 3,597; 4 638; [-0,897; 2,283]; 4,617; [-0,934; 2,103] соответственно для результатов, приведенных на рис. 11, а-ж.

На рисунке 12 приведены распределения локальных чисел Нуссельта на внутренней и внешней поверхности шаровой полости.

Для всех режимов (рис. 12; линия 2) распределение локальных чисел Нуссельта носит волновой характер и соответствует конвективному теплообмену. Распределение локальных чисел Нуссельта имеет один или два минимума и один максимум. Значение осредненного и интервалы изменения локальных чисел Нуссельта на наружной поверхности шаровой полости следующие: 'ый2 = 5,319; 5,319; 9,728; 21,166; -1,697; 16,704; -6,065;

Nu2 е [0,590; 17,227], [0,590; 17,227], [0,212; 19,905], [4,738; 37,226], [-5,712; 1,848], [5,266; 34,606], [-6,316; -5,152] соответственно для результатов рис. 12, а-ж.

Рис. 12. Распределение локальных чисел Нуссельта: а - без учета магнитных сил, 3 = 0, Qy = 0; б - с учетом магнитных сил,

3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных

сил, 3 ^ 0, Qy = 1; д - с учетом магнитных сил, 3 ^ 0, Qy = -1; е - с учетом

магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

На рисунке 13 приведены результаты расчетов поля функции тока.

а б в г д е ж

Рис. 13. Поле функции тока: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил, 3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qy = 0; г - с учетом магнитных сил, 3 ^ 0, Qy = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv = -1

В шаровой полости, в зависимости от режима, образуются четыре (рис. 13, а, б, е, ж) и две конвективные ячейки (рис. 13, в-д). Максимальное значение функции тока: |Ушах| = 8,0810-1; 8,08-10-1; 1,15; 1,46; 7,5510-1; 1,15;

3,73 10-1 соответственно для результатов рис. 13, а-ж.

На рисунке 14 приведены результаты расчетов поля напряженности вихря.

а б в г д еж

Рис. 14. Поле напряженности вихря: а - без учета магнитных сил, 3 = 0, Qv = 0; б - с учетом магнитных сил, 3 = 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 0; г - с учетом магнитных сил, 3 ^ 0, Qv = 1; д - с учетом магнитных сил, 3 ^ 0, Qv = -1; е - с учетом магнитных сил, 3 = 0, Qv = +1; ж - с учетом магнитных сил, 3 = 0, Qv =-1

В полости, в зависимости от режима, образуются четыре (рис. 14, а, б, е, ж) и два вихря (рис. 14, в, г, д). Максимальное значение вихря |штах| = 12,0; 12,0; 13,7; 16,3; 10,4; 15,3; 7,2 соответственно для результатов рис. 14, а-ж.

На рисунке 15 приведены результаты расчетов поля радиальной и меридиональной составляющей магнитной индукции.

а б в г д еж

Рис. 15. Поле магнитной индукции. Радиальная составляющая: а - с учетом магнитных сил, 3 = 0, Qy = 0; б - с учетом магнитных сил, 3 ^ 0, Qv = 0; в - с учетом магнитных сил, 3 ^ 0, Qv = 1; г - с учетом магнитных сил, 3 ^ 0, Qv = -1; д - с учетом магнитных сил, 3 = 0, Qv = +1; е - с учетом магнитных сил, 3 = 0, Qv = -1; ж - меридиональная составляющая магнитной

индукции

Структура поля радиальной составляющей магнитной индукции претерпела значительные изменения по сравнению с результатами, полученными для значений числа Грасгофа 100 и 1000. Для всех случаев радиальная составляющая магнитной индукции в северном полушарии принимает отрицательные значения, за исключением небольшой области вблизи внутренней поверхности полости, а в южном - положительные, за исключением небольшой области вблизи внутренней поверхности полости. Максимальная интенсивность радиальной составляющей магнитной индукции: |5rmax| = 7,1210-4; 6,87-10-4; 7,33■ 10-4;

5,6410-4; 8,3310-4; 4,11 10-4 соответственно для результатов рис. 15, а-е. Структура поля меридиональной составляющей магнитной индукции практически не изменилась по сравнению с аналогичными результатами, представленными выше (рис. 15, ж).

В данной работе все результаты были получены для одной толщины шаровой полости R2 = 2,0. Определенный интерес вызывает исследование влияния толщины шаровой полости на теплообмен и структуру течения жидкости в полости. Результаты, позволяющие проследить это влияние, будут представлены во второй части статьи.

Выводы

Анализ результатов, полученных по предложенной математической модели, позволяет сделать следующие выводы:

- при нагреве шаровой полости снизу безразмерные критерии подобия в расчетной области образуют две или четыре конвективные ячейки;

- для всех значений чисел Грасгофа учет джоулевой диссипации изменяет направление циркуляции жидкости в конвективных ячейках и увеличивает ее интенсивность;

- при Gr = 104 в полости имеет место развитая конвекция, а поле радиальной составляющей магнитной индукции значительно отличается от соответствующих полей при Gr = 102;103;

- математическая модель и полученные результаты могут быть полезными при исследовании тепловых и гидродинамических процессов в недрах Земли и других планет.

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

1. Жарков В. Н. Внутреннее строение Земли и планет / В. Н. Жарков. - М. : Наука. 1983. - 414 с.

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

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

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

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

Modeling Convective Heat Transfer of Electroconductive Fluid in a Spherical Cavity. PART I

S. V. Solovyov

Abstract. In this paper, based on mathematical modeling, the convective heat transfer electroconductive liquids with regard to the internal sources of heat and Joule dissipation in a spherical cavity with heat from the bottom, simulating liquid core of the Earth. The structure of flow, temperature field, magnetic field, as well as the distribution of local Nusselt numbers is investigated.

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

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

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

Solovyov Sergei Viktorovich

Doctor of Sciences (Physics and

Mathematics), Professor

Pacific National University

136, Tikhookeanskaja st., Khabarovsk,

680035

tel.: 8(4212) 37-51-88