Структура нелинейных движений в слое с концентрационными источниками тепла Текст научной статьи по специальности «Физика»

СТРУКТУРА НЕЛИНЕЙНЫХ ДВИЖЕНИЙ В СЛОЕ С КОНЦЕНТРАЦИОННЫМИ ИСТОЧНИКАМИ ТЕПЛА

Е.С. Братчикова

Пермский государственный педагогический университет, 614990, Пермь, Сибирская, 24

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

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

ВВЕДЕНИЕ

Исследованию термокапиллярной неустойчивости плоского слоя бинарной жидкости с внутренними источниками тепла посвящены работы [1-3]. В работах [1, 2] рассмотрены равномерно распределенные источники тепла. Однако в некоторых физических процессах (радиоактивный распад, селективное светопоглощение, протекание химических реакций) среда является многокомпонентной и распределение источников тепла неоднородно. Для таких сред интенсивность тепловыделения Q можно положить пропорциональной концентрации активной компоненты. В [3] рассмотрены внутренние источники тепла, тепловыделение которых зависит от концентрации тепловыделяющей компоненты, но пренебрегается возможным изменением тепловых источников в ходе процесса.

© Братчикова Е.С., 2007

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

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим в условиях невесомости механическое равновесие плоского слоя жидкой несжимаемой бинарной смеси с концентрационными источниками тепла. Форма области и граничные условия показаны на рис. 1.

Будем считать, что мощность тепловыделения Q пропорциональна концентрации активной компоненты С . Пусть слой ограничен твердой границей г = 0 и свободной недеформируемой поверхностью г = И . Твердая граница слоя изотермическая. Концентрация тепловыделяющей компоненты на твердой границе постоянна, а на свободной границе отсутствует ее поток. На этой границе выполняется ньютоновский закон теплоотдачи, коэффициент поверхностного натяжения зависит от температуры, и имеет место равенство вязких и термокапиллярных сил.

Запишем уравнения плоского конвективного движения, введя функцию тока ¥ и вихрь скорости ю, которые связаны с компонентами скорости следующим образом:

В безразмерных переменных уравнения плоского конвективного движения имеют вид:

Ввиду того, что одна из компонент смеси активна (в жидкости происходит процесс тепловыделения), в уравнении теплопроводности (1.2) появляется добавочный член, описывающий разогрев жидкости. В уравнение (1.3) для концентрации активной компоненты добавлено слагаемое, описывающее убыль ее вследствие выгорания, пропорциональное концентрации - N 2С . Коэффициент N2 определяет скорость выгорания.

В качестве единицы длины выбрана толщина слоя И , единицы времени - И2/у , скорости - %/И , единицы функции тока - %, вихря - %/И2, концентрации - С0, температуры - цк2, где

Ч = Q0C0/Р%ср , % - и коэффициент температуропроводности, V -

кинематическая вязкость смеси.

Уравнения содержат следующие безразмерные параметры: параметр выгорания N, число Прандтля Рг = у/% , диффузионное

д¥

и = Ух = —; V = V г

д¥ д2¥ д2¥

дх ’ дх2 дг2

(1.1)

д/ дг дх дх дг

(1.2)

Рг дС+ГдС-^¥ дС

а д/ Ье I дг дх дх дг

= ДС - N2С,

(1.3)

д2¥ д2¥

(1.4)

ю = -

дх2 дг:

22

число Прандтля Pгd = у/D , число Льюиса Ье = ^% = Pг|Pгd . Здесь D - коэффициент диффузии.

В соответствии со сделанными выше предположениями граничные условия запишем в виде:

д¥

г = 0: ¥ = 0, — = 0, T = 0, С = 1. (1.5)

дг

г = 1: — = -ЫТ , ^ = Ма—, — = 0. (1.6)

дг дх дг дг

Здесь Ы = аИ/% - число Био, Ма = о£0С0Иъ / срур2% - модифицированное число Марангони.

Краевая задача (1.1) - (1.6) при малых значениях числа Маран-гони имеет стационарное решение, соответствующее состоянию механического равновесия:

С0 = сИ(N - МО , (1.7)

cИN

Т„ =—^_ (| -_1_, + ±^ - сИ(N - Щ. (18)

0 N 2(Ы +1) сИК N2 N cИN

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

Результаты линейного анализа устойчивости состояния (1.7) -(1.8) для числа Льюиса Ье = 1 и N = 1 показаны на рис. 2. В области положительных значений числа Марангони наиболее опасными для монотонной неустойчивости являются возмущения с длиной волны 1 = 1.17, для колебательной неустойчивости - с длиной волны 1 = 2.24, в области отрицательных чисел Марангони для монотонной неустойчивости - возмущения с длиной волны 1 = 4.27 . Поэтому в нелинейных расчетах исследовались движения, развивающиеся в ячейках слоя с длиной волны 1 = 1, 1 = 2 и 1 = 4. Вертикальными отрезками на рис. 2, а отмечены параметры, для которых проведены расчеты нелинейных режимов.

Рис. 2. Нейтральные кривые (а) и зависимость частоты нейтральных колебаний от волнового числа (Ь) для Ьв = 1, N = 1, Ы = 1, Рг = 1

2. МЕТОД РЕШЕНИЯ

Для численного построения решения выделим на плоскости х, г прямоугольную область интегрирования 0 < х < 1, 0 < г < 1. Условия на горизонтальных участках г = 0 и г = 1 сформулированы выше. На вертикальных границах поставим условия такие же, как на вертикальных границах ячейки нормальных возмущений с волновым числом 1, кратным I: будем предполагать их свободными и теплоизолированными:

х = 0, х = I : ¥ = 0,

Г* = 0, * = 0, дС = 0.

дх дх дх

Расчеты нелинейных режимов конвекции проводились методом конечных разностей по явной схеме с чередующимися направлениями на равномерной квадратной сетке. При составлении конечноразностного аналога уравнений пространственные производные аппроксимировались центральными разностями, а производные по времени - односторонними. Оператор Лапласа расписывался по пятиточечной схеме. Значения вихря на твердой границе находились по формулам Тома. Для аппроксимации условий тепло изолированности и запрета потока концентрации тепловыделяющей компоненты использовался трехточечный шаблон. Уравнение Пуассона решалось методом верхней релаксации.

Ь

Стационарное решение находилось методом установления. В начальный момент было задано возмущение вихря на фоне равновесных распределений температуры и концентрации. Для получения решения в квадратной ячейке использовалась сетка 55 х 55, в прямоугольной - сетки 110 х 55 и 160 х 40. Шаг по времени, в соответствии с требованием устойчивости счета, выбирался по формуле At = к1 /40 .

Числа Прандтля, Био, Льюиса и параметр выгорания в расчетах фиксировались: Pr = 1, Bi = 1, Le = 1, N = 1.

В ходе расчетов определялись поля Y , T , C , строились изолинии функции тока, температуры и концентрации, а также определялись максимальное значение функции тока Ym (для стационарной конвекции), амплитуда A максимального значения функции тока -для колебательных режимов конвекции.

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Диаграмма режимов конвекции, полученная для расчетной области 2:1, показана на рис. 3. Сплошные линии изображают зависимость от числа Марангони максимального значения функции тока, штриховые линии - зависимость периода регулярных колебательных режимов. На оси чисел Марангони звездочками указаны значения критических чисел Марангони для монотонного возмущения с длиной волны 1 = 1 и для колебательных возмущений с 1 = 2 и 1 = 4 , штрихпунктирные линии разделяют области с разными режимами конвекции. В случае Ma меньше некоторого критического значения движение в ячейке отсутствует, все начальные возмущения со временем затухают.

Нелинейный режим, который согласно линейной теории должен возникнуть с 1 = 1 при Ma* = 4952 , получен в расчетной области с l = 1 при Ma* = 4736 . При Ma* от равновесия ответвляется стационарное движение, характеризующееся двухвихревой структурой, симметричные вихри которой имеют одинаковую интенсивность и локализованы в верхней области ячейки (режим I). По мере повышения надкритичности увеличивается максимальное значение функции тока вихрей. В поле температуры вблизи поверхности, где происходит подъем жидкости к свободной границе, наблюдаются замкнутые изотермы, имеющие более высокую температуру. В этих же точках концентрация тепловыделяющей компоненты больше, чем на периферии ячейки. В рассматриваемом диапазоне чисел Ма-

рангони этот режим конвекции устойчив, при Ма = 12800 он переходит в режим периодических колебаний. При обратной смене режимов обнаружена небольшая область гистерезиса.

Для анализа других возможных нелинейных режимов длина расчетной ячейки была увеличена в два раза (2:1). В расчетах при Ма** > 4740 в ячейке также возникала стационарная конвекция с

1 = 1.

Рис. 3. Диаграмма режимов движения жидкости в ячейке 2:1 при Ьв = 1, Ы = 1, N = 1. Зависимость амплитуды (сплошные кривые) и зависимость периода колебаний от параметра Марангони (штриховые кривые)

С увеличением числа Марангони при Ма > 6900 стационарная конвекция становится неустойчивой по отношению к колебательному режиму. В расчетной области наблюдалось движение, представляющее две длины волны. В некоторые моменты времени большей интенсивностью обладают вихри одного знака, в другие моменты времени - вихри противоположного знака (режим II). Анализ структуры стационарной и колебательной конвекции показывает, что данный вид колебаний представляет собой пульсации около стационарного состояния. Период колебаний несколько увеличивается с ростом числа Марангони. На рис. 4 штрихпунктирной линией представлена зависимость максимальной функции тока от времени для Ма = 7000. Точками на кривых отмечены моменты времени, для которых приведены карты изолиний функции тока на рис. 5.

Рис. 4. Зависимость максимальной функции тока от безразмерного времени для различных чисел Марангони

Рис. 5. Карты изолиний функции тока для Ма = 7000

Режим регулярных колебаний II при Ма > 7460 становится неустойчивым и переходит в режим нерегулярных колебаний. Нерегулярные колебания наблюдались до Ма » 8000 . Этот режим получается как методом изменения (повышения) числа Марангони, так и наложением малых возмущений на состояние покоя.

При увеличении числа Марангони от Ма = 8000 из состояния равновесия с малыми начальными возмущениями удалось обнаружить колебательный режим (III), представляющий собой сложные периодические колебания, максимальные отрицательная и положительная функции тока которых изменяются не одинаково.

Период колебаний заметно уменьшается с ростом числа Марангони. Анализ карт изолиний функции тока, температуры и концентрации показывает, что в некоторые моменты периода происходит перестройка от четырехвихревой структуры к двухвихревой и обратно. Движение достаточно интенсивное, и на свободную границу из глубины слоя выносится жидкость с большей концентрацией тепловыделяющей компоненты. Изменение максимальной функции тока со временем в случае Ма = 8000 приведено на рис. 4 пунктирной линией, изолинии функции тока показаны на рис. 6.

0 00 0 20 0 40 0 60 0 80 1 00 1 20 1 40 1 60 1 80 2 0

0 00 0 20 0 40 0 60 0 80 1 00 1 20 1 40 1 60 1 80 2 00

Рис. 6. Карты изолиний функции тока для Ма = 8000

1

При дальнейшем увеличении числа Марангони в области 9060 > Ма > 10000 наблюдается новый режим сложных колебаний (IV). Колебания вблизи Ма = 9100 похожи на колебания жидкости в режиме III, при этом происходит небольшая деформация каждого второго периода. На рис. 4 сплошной линией представлена зависимость максимальной функции тока со временем для Ма = 9750 , карты изолиний функции тока для моментов времени, отмеченных точками на этой кривой, показаны на рис. 7.

Рис. 7. Карты изолиний функции тока для Ма = 9750

Период колебаний при переходе от режима III удваивается (рис. 3) - происходит бифуркация удвоения периода. Для иллюстрации на рис. 8 приведены фазовые портреты установившихся колебаний (связь температуры и функции тока) в точке расчетной области с координатами х = 5/22 и г = 8/11 при Ма = 9060 (а) и Ма = 9100 (Ь). При увеличении параметра Марангони период колебаний режима IV постепенно уменьшается.

Рассмотрим изменение режимов конвекции при уменьшении числа Марангони. Режим IV переходит в режим III при Ма = 9070 . Области с одновременным существованием двух режимов не обнаружены.

Уменьшая число Марангони и выбирая за начальное состояние поля функции тока, вихря, температуры и концентрации при установившихся колебаниях с Ма = 8000 , методом продолжения по

параметру удается продвинуться до Ма = 7300 . При этом значении Ма периоды колебаний режимов II и III практически равны.

0.204

0.2

0.196

а Ь

Рис. 8. Фазовый портрет (связь температуры и функции тока) в точке с заданными координатами

Т

Х.2

0.192

00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.

Уд —

00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.

Рис. 9 Изолинии функции тока (а), температуры (Ь) и концентрации (с) при Ма = -14500 и X = 4

На границе колебательной (режим II) и монотонной конвекции (режим I) имеется область чисел Марангони, в которой могут существовать оба режима конвекции. При понижении числа Марангони от развитого колебательного режима (Ма » 6550 ) происходит сме-

а

Ь

с

на формы движения: колебательный режим становится неустойчивым, и в слое возникает стационарная конвекция.

В численном эксперименте для отрицательных значений числа Марангони (Ма = -14500) получен стационарный режим с длиной волны 1 = 4 . В ячейке образуются два вихря разного знака, с повышением числа Марангони по абсолютной величине интенсивность вихрей возрастает. На рис. 9 представлены изолинии функции тока, температуры и концентрации этого движения. Интересно отметить, что в области восходящего потока на поверхности возникает "холодное пятно", хотя в этом месте концентрация активной компоненты наибольшая. Кажущееся противоречие связано с тем, что параметр Марангони определен через температурный коэффициент поверхностного натяжения и удельную мощность тепловыделения. Для отрицательного значения числа Марангони одна из этих величин должна быть отрицательной.

Заключение. В работе рассмотрены надкритические режимы движения бинарной смеси с концентрационными источниками тепла, которые возникают в ячейке слоя со строго вертикальными границами 2:1 в отсутствие гравитации. Обнаружено, что в такой системе при Ма/Ма* < 2 существует стационарная конвекция с длиной волны 1 = 1 и три различных колебательных режима движения жидкости. Изучено поведение жидкости в процессе колебаний, определены диапазоны чисел Марангони, при которых существуют колебательные режимы.

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Андреев В.К., Родионов А.А., Рябицкий ЕА. Возникновение термокапиллярной конвекции в жидком цилиндре, цилиндрическом и плоском слоях под действием внутренних источников тепла // ПМТФ. 1989. № 2. С. 101-108.

2. Брискман В.А., Якушин В.И. Термокапиллярная неустойчивость в слое с внутренними источниками тепла // Термо- и концентрационно-капиллярные эффекты в сложных системах. Екатеринбург: УрО РАН, 2003. С. 34-43.

3. Якушин В.И. Термокапиллярная неустойчивость плоского слоя жидкости с концентрационными источниками тепла // Конвективные течения... / Пермь: Перм. пед. ун-т, 2003. С. 62-74.

4. Yakushin V.I., Bratchikova E.S. Thermocapillary instability and finite amplitude convective flows in a plane liquid layer with concentrated heat sources // Proc. Intern. Conf. "Advanced problems in thermal convection". Perm, Russia, 24-27 November 2003 / Perm State University, 2004. P. 239-242.

5. Братчикова Е.С. О термокапиллярной конвекции в слое бинарной смеси с концентрационными источниками тепла // VII Все-рос. конф. молодых ученых по мат. моделированию и информационным технологиям (с участием иностранных ученых). http://www.ict.ru/ws/YM2006/.

STRUCTURE NONLINEAR MOTIONS IN THE PLANE LAYER OF LIQUID WITH CONCENTRATION SOURCES OF HEAT

E.S. Bratchikova

Abstract. Thermocapillary convection of binary mixture in plane horizontal layer under weightness is studied. We consider the case when one component of mixture is burning off with heat generation. Numerical research shows excitation of stationary convection at critical Marangoni number. In overcritical domain different regimes of oscillating convection are found out and studied. The structure of fluid motion is studied versus the di-mensionless parameters and amplitude curves are built.

Key words: plane layer, concentration sources of heat, Marangoni convection, overcritical regimes.