CC BY
27
Вычислительные технологии
Том 13, № 4, 2008
Неустойчивость двухслойной системы при наличии объемных источников тепла*
В. Б. БЕКЕЖАНОВА Институт вычислительного моделирования СО РАН, Сибирский федеральный университет, Красноярск, Россия e-mail: bekezhanova@mail.ru
An initial stage of the penetrative convection in a system of two immiscible liquids is investigated numerically. Densities of liquids are supposed to depend linearly on temperature and pressure. The heat source function is taken into account in the free convection equations. It is shown that
Введение
Исследования физико-химических и биологических характеристик “возрастного” состава вод в оз. Байкал [1, 2] свидетельствуют о наличии в озере механизма глубокого перемешивания, переносящего поверхностные воды в придонные области. Важными факторами, формирующими конвекцию в естественных водоемах, являются непрерывно идущие процессы обмена теплом между водной массой и окружающей средой. Для возникновения конвекции необходимы такие условия на поверхности водоема, чтобы поверхностные слои становились тяжелее нижележащих. Подобные условия создаются весной и летом, когда прогревание поверхностных вод приводит к увеличению их плотности, летом — во время ночного охлаждения, осенью — когда тепловой баланс поверхности становится отрицательным и поверхностные воды, имеющие температуру выше 4 °С, становятся тяжелее нижележащих [3].
Тепловой баланс и режим поверхности планеты формируется под действием такого мощного природного фактора, как солнечное излучение. Наиболее активные термодинамические процессы имеют место при резких сезонных изменениях потока инсоляции. Вода представляет собой светорассеивающую и светопоглощающую среду в достаточно широком диапазоне спектра солнечного излучения, поэтому необходим учет ее оптических свойств. В работе предпринята попытка оценить влияние процесса “переработки” потока излучения на гидродинамическую устойчивость байкальских вод.
Для описания движения жидкости используются уравнения Обербека—Буссинеска. В уравнение энергии добавлена общая энергетическая функция теплового источника Fw, описывающая поглощение солнечной радиации. Значение Fw определяется оптическими параметрами и моделью распространения излучения в неоднородной среде.
Другим важным фактором, определяющим интенсивность вертикального водообмена и связанных с этим процессов переноса примесей, является плотностная стратифи-
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, (гранты № 08-01-00762, НШ 2260.2008.1) и комплексного интеграционного проекта 2.15 СО РАН.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
кация вод озера. Если толщина слоя невелика, то отклонения плотности, вызванные изменениями давления, можно не учитывать. При изучении же процессов, происходящих в глубоководных водоемах, следует принять во внимание, что возникающие перепады давления могут существенно влиять на распределение плотности, а следовательно, и на конвективные процессы. В глубоких водоемах возникает эффект сжимаемости воды при высоких давлениях. Уравнение состояния жидкостей примем в виде
где индекс ] нумерует жидкости, первая жидкость расположена внизу; р0 — некоторое постоянное значение плотности (для воды характерное значение р0 = 999.972 кг/м3); в? — коэффициент теплового расширения жидкости; 0? — температура; 7? — барический коэффициент расширения; 00?, р0? — постоянные средние значения температуры и давления. Такая зависимость плотностей от температуры и давления с малыми коэффициентами в? и 7? позволяет считать слои жидкостей слабо сжимаемыми средами [4].
1. Постановка задачи
Рассмотрим гравитационную тепловую конвекцию в системе двух несмешивающихся жидкостей с общей границей раздела, ограниченных снизу твердой стенкой и сверху свободной поверхностью (рис. 1). Оси х и у находятся в плоскости верхней границы слоя, а ось г направлена вертикально вниз. Точка г = 0 соответствует свободной поверхности, г = г* — поверхности раздела, г = I — нижней границе слоя (твердая стенка). Поверхности Г и Г определяются уравнениями /1,2(х, ¿) = 0, где х = (х,у,г); в частности, для рассматриваемого ниже покоя /1 = г — г*, /2 = г — I.
В каждой области П? справедлива система уравнений Обербека—Буссинеска
(1)
ди '
аіу ^ = о, -д£- + V,- ■ ▽в, = хз △ из + Е*>з¿),
(2)
/
ад)
Рис. 1. Схема течения
где V = (uj,Vj ,Wj) — вектор скорости і-й жидкости, Xj — коэффициент температуропроводности, Fwj (г, і) — мощность внутренних сил. В случае объемного поглощения (проникновения радиации в среду) функция теплового источника описывает внутреннее тепловыделение и определяется по закону [5]
^(-М) = а«2Яехр(—«2,г), ^і(г,і) = ^2(2*, і) ехр(—«іг).
(3)
Здесь о — числовой параметр, представляющий собой отношение интенсивности солнечной радиации к радиационному балансу на свободной поверхности, 0 < о < 1. Далее,
— показатель ослабления солнечной радиации в жидкости (показатель поглощения),
— радиационный баланс поверхности, ^ — динамический коэффициент вязкости, : (0, 0,$), д — ускорение свободного падения. Плотности pj• имеют вид (1).
На твердой стенке задаются температура и условие прилипания
Я
g
01 = 01, VI
На поверхности раздела Г1 при г = г*:
0 при г = /.
VI = V2, 01 = 02, РіП = Р*2П, VI ■ п = V
, д0і _ . 502 = к2тт~ д п дп
На свободной поверхности Г2 при г = 0:
V ■ П = Уп,
Р2 ■ П + Pga
д02
П
0 к2^ + Ь(02 — ^аэ) = ^,
дП
(4)
(5)
(6)
где п — нормаль к поверхности Г, — скорость Г в направлении нормали,
р = — pj + 2^-Dj — тензор напряжений в жидкости, Dj — тензор скоростей деформации векторного поля Vj, р§а8 — давление газа, kj — коэффициент теплопроводности жидкости, Ь — коэффициент межфазного теплообмена, 0ёав — температура газа, Q — заданный поток тепла через свободную поверхность. Учитывая, что на прогрев водяной толщи расходуется энергия, равная Я(1 — а), получим Q = Я(1 — а) — — ^, —
потери тепла на испарение (конденсацию), ^ — потери на конвективный теплообмен.
2. Равновесное состояние
В состоянии механического равновесия vj• = 0 и производные по времени Qj^ = pjt = 0. Из уравнения энергии следует, что 6j (г) имеют следующий вид:
01 (г) = — Аехр(—«іг) + сііг + С21, А = аЯ^2 ехр(—«2^*);
02(г) = —
аЯ
«2X2
«1X1
ехр ( «2 г) + Сі2^ + С22.
Постоянные c1j и с^- (і = 1, 2) определяются из граничных условий:
С21 = 01 — Сц/ + А ехр( —«і/), С12
D — с22
г* + В
аЯ
Х2
(7)
(8)
Сіі = ^ (— ехр(—«2^*) + С12 ) — А«1 ехр(—«іг*), к1 V Х2
D = 0і + — ехр(—«2^*) (--------В ) + Аехр(—«і/) + Аехр(—«іг*)(«і(/ — г*) — 1),
с22 =
В + г*
Ь(В + г*) + &2
О к2 , к2
В = — I — — г* к1 к1
„ . / оД
ф + к2 I----------------+
Б
+
оД
Ь
\Х2 В + г*) «2X2
Рассматривалась зависимость распределения температуры при близких значениях показателей поглощения «1 и «2 от положения поверхности раздела. Вообще говоря, показатель поглощения « не является постоянной величиной, а зависит от природы, состояния вещества, от длины волны проходящего излучения и практически не зависит от мощности светового пучка. Поскольку задача рассматривается применительно к естественным глубоким водоемам, то учитываем коротковолновую радиацию и пренебрегаем длинноволновой, основная часть которой поглощается атмосферой. В таких условиях можно считать «j постоянными и «j ~ 0.2 — 0.4. На рис. 2 даны графики функции температуры при I = 730 м (средняя глубина Байкала). За счет светопоглощающих свойств верхней жидкости влияние солнечного излучения на температурный режим нижнего слоя тем слабее, чем больше толщина верхнего слоя.
Полученные распределения температуры можно сравнить с результатами экспериментальных наблюдений, проводимых на протяжении многих лет на оз. Байкал. По приведенным в [3, 6] значениям 61, Д, ф, 6(г) для отдельных месяцев и зон Байкала восстанавливаются значения «j, г* и о, при которых можно получить истинный профиль температуры.
Из уравнения импульса находятся функции давления pj:
Р1(г)
Род — а2 (с21 — 60 + с11г)
а1
+ Ро +
а2сп\
)
+
+
Аа2
ехр((а1 — «1)г) + ехр(—а1г)
Р2(г)
/ а4 оД
^ — «1
Род — а4 (с22 — 60 + с12г) аз
+ Ро +
а4с12 \
а3 )
+
+
\Х2 «2 (а3 — «2 )
Здесь а1 = Ро71$, а2 = РоАд, аз = Ро72д, а4 = роАд.
ехр((аз — «2)г) + ехр(—азг).
100 200 300 400 500 600 700 _ „ 100 200 300 400 500 600 700 ^ ^
а б
Рис. 2. Графики функций температуры: а — г* = 200 м, б — г* = 50 м; кривая 1 — *1 > *2, кривая 2 — *1 < *2
Постоянные ^1, ^2 находятся из граничных условий и равны
Ро£ — а4 (с22 — 0о + с12^*) + р + а4с12^\
аз
/ Ч / ^4^22 '•'О I '-'12~*У . . ^4ч-'12 \ .
ехр(азг*) ( -------------------а-------+ ро +—азз— ] +
+ (----------------г ехр((аз - «2)2*) + ^
\Х2«2(аз - «2) У
^2 = — е2)
ехр((а1 - аз)г*) - б1,
где
, Ч I Ро^ - а2(с21 - 0о + С11 £*) а2СП\ . Аа2 ,, Ч Ч
61 = ехр(а!,г*)----------------------------------------+ Ро +-------^ +-----------------ехр((а1 - «1 )^*),
а1 а| ) а1 - «1
Ро£ - а4(с22 - ^ . . а4с12 , а4аД
62 =---------------:--------------+ Ро +------------^ +
аз а3 Х2«2(аз - «2)
Функции р, выпуклы вниз и близки к линейным.
Итак, получено стационарное решение р,, 0, краевой задачи (2)-(6), соответствующее состоянию механического равновесия V, = 0.
3. Задача о малых возмущениях равновесия
Сформулируем задачу об устойчивости механического равновесия относительно малых возмущений. Для этого введем определяющие безразмерные параметры. В качестве характерного масштаба длины выберем толщину I* = г* верхнего слоя, в качестве масштаба температуры — разность 0 = 0о - 01. За масштаб скорости примем скорость конвективного всплытия нагретой частицы жидкости V* = \/^1*в0. Для плотности и давления используем масштабы ро и рог'2 соответственно. Температуру будем отсчитывать от температуры нижней границы 01, а давление — от гидростатического. Кроме того, определим модифицированную функцию теплового источника Fj* = v*l*FWj/к,.
Введем безразмерные переменные £ = (С,П,С), т такие, что
х = (х, у, г) = £1*, Ь = — т, р, = роV3pj, у, = v*vj, 0, = 00'.
V*
Здесь р,, V,, 0, — безразмерные функции давления, скорости и температуры соответственно.
Примем в качестве единиц измерения коэффициентов V, х, в, 7, «, к их средние арифметические значения
V! + Х1 + Х2 п в1 + в2 71 + 72 «1 + «2 , к1 + к2
^ = —2— , Х* =--------2-----, в* =-2— , 7* = —2— , «* =--------2---, = —2— .
Тогда безразмерные параметры V,, х,, в,, 7,, «,', к, вычисляются по формулам
/ / Х, /Э/ в, / / «, 7 /
, = — , Х, = — , в = 7^ , 7, = — , «, = — , к = т-
V* Х* в* 7* «* к*
Задача (2)-(6) характеризуется безразмерными параметрами: е# = в*0, = 7*Р0v3,
Ре = 0о/0, = pо/(роV3), ¿* = х*/(l*v*) — число Фурье, /* = 12аЯ«*/(0к*) — параметр
тепловыделения, р* = V*/(/*'У*) — параметр кинематической вязкости (обратное к числу Рейнольдса), Ба = 1/(р*5*) — число Рэлея.
ПУсть (£,г) = у(£,т) + 5*У,(£,т), р?(£,г) = р(£,г) + рАЯ?(£,г), 0ф-(£,т) =
= 0,(£,т) + Т?(£,т), где У, = (и,, V,, Ж?),Р?, Т? — возмущения, а V,,р?, 0? — основное решение. Вид функций УфРф0ф описывающих возмущенное движение, выбран для удобства последующих преобразований. Линеаризуя полную задачу, получим для возмущений скорости, температуры и давления в каждой из жидкостей следующую краевую задачу (“штрихи” опущены):
Ц?? + + Ж? = 0, Тфт + 5*Л?Ж? = 5* А Т?,
~~ = -Р;?? + Аи?, ? = -Рэт + А^?, (9)
р* р*
+ АИф - Ба Т? - ^ Р?,
Р*
где ^ ^ (С) = д0?/дС;
— на твердой стенке
С = ///* : и = V = Ж1 = 0, Т = 0; (10)
— на поверхности раздела
С = 1 : ^1 = и2, VI = V, Ж1 = Ж2, Д1т = 5*ЖЬ
Р^1(^1С + Ж1?) = Р2^(и2С + Ж2?), Р^1(^1С + Ж1п) = Р2^(У2с + Ж2п),
Т - ^ Д1 = Т2 - ^ Яь (11)
дС д(
, дТ1 , дТ2 + ( д202 , д20Л р =0
к1Ж - к2 аГ + (*2 ас2 - к1 ас^у*Й1 = 0
Р1 - Р2 + 2(р21“гИ2( - Р^И-'^) = Ба (рк - р2< - ^ (р2 - Р1 Й1,
где Д1 = Д1(^, п, т) — локальное отклонение поверхности раздела от ее невозмущенного состояния по нормали;
— на свободной границе
С = 0 : и~2£ + Ж2? = 0, + Ж^ = 0, ^2т = 5*^2,
Р2 - 2^2С = Ба (^ + Р2^ Д2, (12)
дТ2 + ^72 Д2 + Б1(Т2 - №) = 0,
дг д(2
где Д2 = Р2(£, П, Т) — возмущение свободной границы, ВІ = 6/*/&2 — число Био.
Рассмотрим нормальные возмущения, пропорциональные ехр[г(аі£ + а2п — Ст)], где аі,а2 — волновые числа вдоль осей х и у соответственно, С = Сг + ¿С — комплексный декремент. Для амплитуд нормальных возмущений получим спектральную краевую задачу, к которой применимо преобразование Сквайра: ^ = га1Ц?- + га2У/. После преобразований система (9) примет вид
Z, + Ж, = 0, -¿СТ,- + Ж, = 5* (Т" — а2Т,),
—с
Т = а2Р, + Z,:| — а2^, (13)
— Ж, = — р + Ж'' — а2 Ж, — Ия Т, — ^ Р,
где а2 = а2 + а2 — квадрат модифицированного волнового числа, а “штрих” обозначает дифференцирование по переменной £ •
Граничные условия (10)—(12) перейдут в следующие:
С = ///* : Ті = 0, Жі = 0, Ті = 0; (14)
С = 1 : Ті = Т2, Жі = Ж2, Р2^2 — а2 Ж) = рі^і^і — а2Жі),
Р1 — Р2 + 2(Р2^2Ж2 — Рі^іЖі) = Иа ^рі — р2 — ~ (Р2 — Рі^ ДЪ
г? _ *тл/ ^ д0і„ д02„ п,ч
Ді = "с Жі, Ті — "д^^і = — "д^ (15)
^2^2 — кіТ + (М2' — кі^'^Ді = 0;
/О
с = 0 : Т2 — а2Ж2 = 0, Д = ^ Ж,
с
Р2 — 2Ж2 = Иа (р2 + Р2/^“)Д2, (16)
Т2 + ^2' Д + ВІ(Т2 — ^2^2) = 0.
Краевая задача (13)—(16) является задачей на собственные значения относительно комплексного декремента С. Для устойчивости равновесного состояния р,,0, по отношению к малым возмущениям вида необходимо и достаточно, чтобы у всех собственных значений С мнимая часть Сі была отрицательной.
4. Численное решение
Полученная спектральная задача (13)—(16) решается методом ортогонализации [7].
Исследовалась устойчивость системы горизонтальных слоев слабо сжимаемых жидкостей с общей поверхностью раздела при следующих значениях параметров: вёав = 12 °С, р§аз = 101300 Па, V* = 1.57 ■ 10-6 м2/с, х* = 1-323 ■ 10-7 м2/с, в* = 8.57 ■ 10-6 К, к* = 0.559 Вт/(м-К), ж* = 0.36 м-1, 7* = 1.40 ■ 10-8 Па-1. Приведенные данные соответствуют значениям параметров для воды в оз. Байкал. При указанных значениях физических параметров найдена зависимость С = 1т С от волнового числа а.
Расчеты проводились для средних глубин Южного (/ = 810 м), Центрального (/ = 803 м) и Северного Байкала (/ = 564 м). При этом учитывались средние значения потоков тепла Q и радиационного баланса Л, характерных для безледного периода указанных зон Байкала.
На рис. 3 изображены кривые СДа), полученные при различных значениях г*. Укажем критические значения волновых чисел а*, при которых наступает неустойчивость, и соответствующие им критические длины волн в размерных переменных А*, а также числа Бк
а) г* = 50 м
Северный Байкал — а* = 1.01, А* = 311.6 м, Б1 = 0.08;
Рис. 3. Инкременты нарастания возмущений Сг(а): а — при г* = 50 м, б — при г* кривая 1 — Северный Байкал, 2 — Южный Байкал, 3 — Центральный Байкал
а
Южный Байкал — а* = 1.10, А* = 284.2 м, Б1 = 0.10;
Центральный Байкал — а* = 1.19, А* = 261.7 м, Б1 = 0.11;
б) —г* = 200 м
Северный Байкал — а* = 3.71, А* = 338.4 м, Б1 = 0.31; Южный Байкал — а* = 4.19, А* = 300.0 м, Б1 = 0.40;
Центральный Байкал — а* = 4.47, А* = 281.1 м, Б1 = 0.42.
= 200 м;
Из полученных результатов видно стабилизирующее влияние теплообмена на устойчивость равновесия. Однако, чем меньше толщина верхнего слоя, тем активнее термодинамические процессы, протекающие в нем, и это оказывает дестабилизирующее влияние на устойчивость.
Граница устойчивости определяется из соотношения СДИа) = 0. Случаю Сі = 0 соответствуют нейтральные возмущения. Полагая в задаче (13)—(16) С = 0, получим нейтральные кривые устойчивости. При расчетах изменялось значение чисел Био, характеризующее интенсивность теплообмена с окружающей средой, положение поверхности раздела г*, а значение I полагалось равным 730 м (средняя глубина Байкала) во всех случаях.
На рис. 4 представлены зависимости чисел Рэлея от волнового числа (нейтральные кривые). Для различных значений чисел Био приведем критические числа Рэлея Иа*, которые являются минимальными значениями на соответствующих нейтральных кривых, и критические волновые числа а*, при которых достигаются значения Иа*:
а) г* = 50 м
Ві = 0.2, Иа* = 2819.16, а* = 0.74;
Ві = 1, Иа* = 3474.20, а* = 0.97;
Ві = 2, Иа* = 3921.71, а* = 1.18;
б) г* = 200 м
Ві = 0.2, Иа* = 3004.27, а* = 1.06;
Ві = 1, Иа* = 3561.67, а* = 1.45;
Ві = 2, Иа* = 3994.19, а* = 1.70.
Видно, что при уменьшении числа Био критические числа Рэлея убывают, а область неустойчивости смещается в область меньших по значению волновых чисел. Полученные результаты позволяют говорить о стабилизирующем влиянии теплообмена на свободной границе. Кроме того, с ростом толщины верхнего слоя также увеличивается зона устойчивости. На рис. 5 показана зависимость критического волнового числа а* от координаты поверхности раздела г*.
Рис. 4. Нейтральные кривые: а — при г* = 50 м, б — при г* = 200 м; кривая 1 — Ві = 2, 2 — Ві = 1, 3 — Ві = 0.2
а
а. а
3 ■■
2 ■
1 ■
О
—------■----■----■----■----■----■--►
100 200 300 400 500 600 700 г„ м
Рис. 5. Зависимость критического волнового числа а* от положения поверхности раздела z* при Bi = 0.2
Полученное решение сравнивалось с решением аналогичной задачи для слоя жидкости толщины l. В двухслойной системе неустойчивость наступает при меньших волновых числах. Укажем критические числа Рэлея при l = 730 м:
Bi = 0.2, Ra* = 3106.18, а* = 1.24;
Bi = 1, Ra* = 3634.87, а* = 1.71;
Bi = 2, Ra* = 4004.73, а* = 1.97.
Заключение
Найдено решение, описывающее состояние покоя системы. На основе численных расчетов показано, что состояние механического равновесия двухслойной системы со свободной границей при наличии объемных источников тепла является неустойчивым. Увеличение толщины верхнего слоя и интенсивность теплообмена на свободной поверхности имеют стабилизирующее влияние. Построены нейтральные кривые и найдены критические числа Рэлея.
Автор выражает благодарность В.К. Андрееву за полезные замечания.
Список литературы
[1] Гранин Н.Г., Шимараев М.Н. К вопросу о стратификации и механизме конвекции в Байкале // Докл. РАН. 1991. Т. 321, № 2. С. 381-385.
[2] Шимараев М.Н., Грачев М.А. и др. Международный гидрофизический эксперимент
на Байкале: процессы обновления глубинных вод в весенний период // Докл. РАН. 1995.
Т. 343, № 6. С. 824-827.
[3] Верволов В.И., Сокольников М.Н., Шимараев М.Н. Гидрометеорологический режим и тепловой балаес озера Байкал. М.; Л.: Наука, 1965. 373 с.
[4] Мосеенков В.Б. Качественные методы исследования задач конвекции вязкой слабо сжимаемой жидкости. Киев: Ин-т математики НАН Украины, 1998. 280 с.
[5] Красс М.С., Мерзликин В.Г. Радиационная теплофизика снега и льда. Л.: Гидрометео-издат, 1990. 261 с.
[6] Shimaraey M.N., Verboloy V.I., Granin N.G., Sherstyankin P.P. Physical limnology of Lake Baikal: a review. Irkutsk—Okayama. 1994. 80 p.
[7] Годунов С.К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений // Успехи мат. наук. 1961. Т. 16, № 3(99). C. 171-174.
Поступила в редакцию 13 июля 2007 г.
